- Email: [email protected]
Design of reduced molecular models by integral equation theory
Chemical Physics Letters 418 (2006) 202–207 www.elsevier.com/locate/cplett
Design of reduced molecular models by integral equation theory Stefan M. Kast *, Wulf Hauptmann 1, Bernd Schilling Physikalische Chemie, Technische Universita¨t Darmstadt, Petersenstraße 20, 64287 Darmstadt, Germany Received 2 September 2005; in final form 24 October 2005 Available online 18 November 2005
Abstract Molecular integral equation theory is inverted with respect to the calculation of effective potentials for a coarse-grained, reduced molecular model composed of fewer sites than the reference. The methodology is applied to an aqueous solution of cyclohexanol as a model case, the quality of the effective potentials is studied by extensive molecular dynamics simulations. Practical guidelines for the choice of reduced model sites are worked out and perspectives for mesoscale modelling are discussed. 2005 Elsevier B.V. All rights reserved.
1. Introduction The development of methods for bridging length and time scale limitations in molecular simulations is currently a very active area of research. Conceptually, with the present-day computational resources a systematic coarse-graining of detailed atomic models represents a necessary and useful approach for studying systems on the mesoscale. In the soft matter sciences, particularly for problems of colloids and polymers, coarse-graining procedures are well established: based on a theorem proved by Henderson more than 30 years ago [1], equating distribution functions for a reference and a reduced model leads to a unique coarse-grained model effective potential by inversion of the distribution, typically computed from molecular simulations. The inversion process is usually performed by iterative molecular simulations [2–4], by inverse Monte-Carlo in conjunction with dissipative particle dynamics simulations [5], or by inversion of the atomic Ornstein-Zernike integral equation within the hypernetted chain (HNC) approximation [6,7]. For biomolecular problems, however, where strongly directional forces dominate the scenarios, reduced model *
Corresponding author. Fax: +49 6151 164298. E-mail address: [email protected] (S.M. Kast). 1 Present address: Technische Chemie, Technische Universita¨t Darmstadt, Petersenstraße 20, 64287 Darmstadt, Germany. 0009-2614/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.10.125
development is far less common. Procedures typically start heuristically (see e.g. [8,9]) and imply checking/adapting the results by explicit simulations. In the case of very large systems where the target observable cannot even be computed from reference simulations, we need another, systematic approach for effective potential development: molecular integral equation theory could provide the key to these difficulties, since it is possible to apply the uniqueness theorem [1] and to invert the potential directly without simulation effort for various reduced models that can retain enough structural information to represent strongly anisotropic interactions. Of course, the necessarily approximate correlation functions from an integral equation treatment imply errors in the computed effective potentials. As a proof of concept, we present in this Letter an integral equation-based methodology for molecular coarse-graining, i.e. reduction to fewer interaction sites than the reference system, on the example of solute–solvent interactions. We test the resulting strategy with a simple amphiphilic model system, cyclohexanol in aqueous solution by extensive molecular dynamics simulations. To the best of our knowledge no such attempt on the basis of integral equations alone has been described in the past. Since the present work represents only an intermediate, though necessary step towards a general coarse-graining approach, the next stages of development are discussed in the concluding section.
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
solution to the solute–solvent RISM–HNC equations for the reference system, and c2 from analogously requiring according to (2)
2. Theory and methods 2.1. Model reduction by integral equation theory The one-dimensional reference interaction site model (RISM) integral equation [10,11] for fluids composed of molecules with sites a and c that is used in this work reads hvv ¼ xvv cvv xvv þ xvv cvv ðqhvv Þ;
ð1Þ
where h = (hac(r)) is the matrix of the total correlation functions (the radial distribution function is g = h + 1), c = (cac(r)) are the direct correlation functions, x = (xac(r)) are the intramolecular correlation functions (for rigid molecules these are normalised Dirac delta functions constraining the site distance), q is the site density, r is the distance and the star denotes a convolution product. In contrast to the pure solvent case (superscript ÔvvÕ), for a solute system in infinite dilution (ÔuvÕ) the simpler relation qhuv ¼ xuu cuv vvv ;
203
ð2Þ
ðq^h1 ^v1 Þred ¼ ðq^h2 ^v1 Þ;
ð6Þ
where the carets denote Fourier transforms and the subscript ÔredÕ indicates that only those lines in the product matrix are taken into account that correspond to the solute–solvent site pairs of the reduced model. The structural difference between reference and reduced model is reflected by a modified intramolecular correlation matrix, such that we obtain finally from (2) and (6) ^c2 ¼ x ^ 1 ^ 1^c1 Þred . 2 ðx
ð7Þ
Inverse Fourier transform of the latter direct correlation function matrix and inserting in (5) yields a direct, noniterative and numerically stable route to the effective potentials u2 using only structural information (in the case of rigid models) and the direct correlation functions of the reference system.
holds where vvv ¼ qxvv þ q2 hvv ;
ð3Þ
the solvent susceptibility (density correlation function), can be computed beforehand and independently of any solute. These equations must be supplied with a closure relation hac ðrÞ þ 1 ¼ expðbuac ðrÞ þ hac ðrÞ cac ðrÞ þ Bac ðrÞÞ;
ð4Þ
with the pair potential uac and b = 1/kT where k is the Boltzmann constant and T is the absolute temperature; Bac are the bridge functions (for site-averaging theories like the RISM approach these functions have not the same diagrammatic meaning as in orientation-preserving theories). In the HNC approximation that is used exclusively here for both solvent–solvent and solute–solvent problems (extended RISM, RISM–HNC) B identically vanishes. Our goal in this Letter is to demonstrate a proof of concept that a coarse-graining model reduction strategy based on integral equation theory is feasible even within the simple RISM–HNC approximation under well-defined conditions. We focus exemplarily on the solute–solvent case where we intend to derive effective potentials for a reduced solute model composed of fewer sites than a reference model by requiring identity between total solute–solvent correlation functions for reduced sites and corresponding atoms of the full model. Direct inversion of the closure with respect to the potential by inserting a modified direct correlation function from inverting the RISM integral equation is numerically difficult because of the r ! 0 singularity of both potential and ln(h + 1). A more practical route is available by considering the difference between the known full potential, u1, and reduced model effective potential, u2, for the set of reduced model sites Duac ¼ u2;ac u1;ac ¼ b1 ðc2;ac c1;ac Þ;
2.2. Models and numerical methods We tested the methodology on a simple model system: a rigid cyclohexanol model as an example for an amphiphilic solute system containing both hydrophobic and hydrophilic regions, solvated in TIP3P water [12] under ambient conditions. The atom numbers and the location of the reduced sites are shown in Fig. 1, the atomic coordinates and force field parameters for the reference system are given in Table 1. The structure was modelled and optimised by Sybyl [13]. A simple Lennard-Jones/Coulomb reference potential function was applied using Lorentz–Berthelot combination rules. The Lennard-Jones parameters were taken from the Amber4.1 force field [14], partial charges resulted from fitting to the electrostatic potential derived from HF/6-31G* quantum-chemical calculations with GAUSSIAN94 [15]. We explicitly investigated three reduced models for comparison with the full 19-site system: reduced sites R1–R5 (5-site), R1, R2, R5 (3-site), and the 1-site model comprising only R2. The effective potentials derived by 1D-RISM–HNC theory are spherically symmetric and given numerically on a grid that can directly be used within the DL_POLY simulation package [16]. Except for the
ð5Þ
as derived from the closure (4) under the condition of total correlation function identity. We take c1 directly from a
Fig. 1. Atom and reduced site numbers of the cyclohexanol models. The 5-site model consists of sites R1–R5, the 3-site model of R1, R2, R5, the 1site model comprises R2.
204
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
Table 1 Structural data and potential parameters (full model: Coulomb charge q, Lennard-Jones e, r) for the cyclohexanol model ˚ ˚ ˚ Atom x/A y/A z/A q/e e/kJ mol1 C1 (R1) H2 H3 C4 H5 H6 C7 (R2) H8 O9 (R3) H10 (R4) C11 H12 H13 C14 (R5) H15 H16 C17 H18 H19
1.0530 1.5410 1.4340 0.4820 0.7210 0.9620 1.0430 0.5940 2.4630 2.9450 0.6910 1.0760 1.1770 0.8440 1.3260 1.0780 1.4040 0.9750 2.4990
1.1261 1.3741 1.8231 1.2941 2.3231 1.1171 0.2931 0.5131 0.4551 0.2851 1.1639 1.8549 1.4119 1.3309 1.1519 2.3629 0.3309 0.5549 0.4419
0.4050 0.5530 1.1710 0.2640 0.0530 1.2400 0.7800 1.7640 0.9340 0.1310 0.3810 1.1480 0.5770 0.2410 1.2170 0.0710 0.8050 1.7950 0.8740
1-site model all energy functions are long-ranged (see below) which necessitated the application of a shifted-force modification: Reduced potentials to be used in simulations were multiplied by the CHARMM function for long-range ˚. electrostatics [17] (1-(r/rc)2)2, with a cutoff distance of 12 A
0.02990 0.00747 0.01621 0.25237 0.06025 0.07063 0.45409 0.05127 0.76214 0.43158 0.25189 0.06021 0.07053 0.02924 0.00758 0.01649 0.09673 0.03880 0.01887
0.44964 0.06453 0.06453 0.44964 0.06453 0.06453 0.44964 0.06453 0.69871 0.19211 0.44964 0.06453 0.06453 0.44964 0.06453 0.06453 0.44964 0.06453 0.06453
˚ r/A 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953 3.39967 2.64953 3.00000 0.40000 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953
Lennard-Jones interactions were also treated by the appropriate CHARMM shifted-force expression and truncated ˚ . Potentials were handled in the same way in beyond 12 A the simulations and in the integral equation treatment, ˚, except that small Lennard-Jones parameters (rH = 0.4 A
Fig. 2. Full atomic model and reduced model effective potentials between site R2 and water oxygen (top), R2 and water hydrogen (bottom) for various reduction stages.
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
eH = 0.1921 kJ mol1) [18] were attributed to the hydrogen sites in order to avoid numerical singularities within the RISM theory. Molecular dynamics simulations were performed in the isothermal/isobaric ensemble [19] (300 K/1 bar) with a temperature coupling constant of 0.5, and 2 ps for the barostat, using a time step of 2 fs. Long-range electrostatics was treated by Ewald summation [20] with 8 · 8· 8 recipro˚ and a charge cal space vectors, a real-space cutoff of 12 A 1 ˚ smearing parameter of 0.3 A . A single solute molecule was embedded in a pre-equilibrated water box of 2000 molecules, equilibrated for another 230 ps and finally sampled for 1.36 ns (1-site), 1.3 ns (3-site), 1.58 ns (5-site), 1.38 ns (full model). Distribution functions were computed on a ˚ spacing (radial averages) and 1 A ˚ 3 spacing grid with 0.1 A (3D density). The RISM solutions were obtained on a log˚ to a arithmic grid of 512 points ranging from 5.98 · 103 A ˚ maximum distance of 164.02 A at the same average density ˚ 3) and temperature as in the simulations. For (0.03291 A the solvent–solvent case, no phenomenological adjustment of the dielectric constant was applied. The nonlinear system of equations was solved using a variant of the Ômodified direct inversion of iterative subspaceÕ (MDIIS) method [21] for coarse optimisation of the correlation func-
205
tions and the traditional Picard iteration scheme for ultimate convergence to a threshold of max(Dc) < 107(vv) and max(Dc) < 104(uv), where Dc means the difference of the direct correlation functions between two successive iteration steps. 3. Results and discussion The resulting effective potentials are shown in Fig. 2 exemplarily for the site R2 that is preserved in all models. As can be seen, the potentials are smooth enough in any case for a direct numerical tabulation within the simulation programme. Except for the 1-site model that is dominated by the overall dipole of the solute, all potentials are apparently long-ranged as expected, but the asymptotics do not follow a simple power law. A shifted-force treatment using truncation is therefore necessary in this case. Apparently, the effective potentials can substantially change their ÔpolarityÕ, i.e. change from overall repulsive to repulsive/attractive upon switching to a different number of reduced sites. The quality of the resulting models can only be assessed by investigating the corresponding simulated distribution functions. As shown in Fig. 3 again for the case of site R2, the initial peak in the radial distribution functions near
1.4 1.2
g (R2–O)
1.0 0.8 MD, full MD, 5–site
0.6
MD, 3–site 0.4
MD, 1–site RISM–HNC
0.2 0.0 1
2
3
4
5
6
7
8
9
10
11
9
10
11
1.4 1.2
g (R2–H)
1.0 0.8 MD, full 0.6
MD, 5–site MD, 3–site MD, 1–site RISM–HNC
0.4 0.2 0.0 1
2
3
4
5
6 r/Å
7
8
Fig. 3. Full atomic model and reduced model radial distribution functions between site R2 and water oxygen (top), R2 and water hydrogen (bottom) for various reduction stages, simulated data and integral equation results.
206
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
the onset of the reference distribution (full model) is only resolved by the 5-site model that preserves the hydroxy group. The RISM–HNC solution also shows no significantly peaked structure in this range. The 3-site and the 1-site models are not able to resolve the reference structure of the distribution functions adequately. For the 5-site model however, the quality of the correspondence is very good, remaining peak shifts can be correlated with the application of the shifted-force modifications. The same statement holds for all the other sites, as exemplarily shown for the OH group oxygen g functions in Fig. 4. We can correlate the quality of the reduced models with not too strong differences between original and reduced potentials. In a number of further experiments we found that the quality of the 1-site and the 3-site models cannot be improved by using the ÔexactÕ solvent susceptibility (having almost no influence on the results) or bridge functions extracted from the simulation of the reference system (cf. [22]). In the latter case, assuming the same bridge function for corresponding sites in the reference and the reduced models is therefore not sufficient; the local molecular structure has to enter into a suitable bridge model expression if improvements beyond the HNC approximation are sought.
In summary, useful rules for the rational development of reduced molecular models by RISM–HNC inversion are: (1) strongly polar groups should be preserved in the reduced models, weakly polar moieties can safely be compressed into substantially fewer sites and (2) the quality of the resulting model hinges upon the form of the reduced potential that should not deviate too severely (in a sense of changing polarity) from the original potential. In a final check, we calculated an observable that did not enter the original, radial g function-based derivation, namely the full threedimensional density distribution of water oxygen around the entire solute. As shown in Fig. 5 (image created with MOLCAD [23]) we find an overall good quality of the agreement, particularly the solvent cavity size and form is well represented, as well as the distinct coordination structure of the hydroxy function. Though it is clear that the quality strongly depends on the choice of the location of the reduced sites, the deficiencies of the RISM–HNC approximations can be eliminated by the rules formulated above. 4. Concluding remarks We have demonstrated that RISM–HNC theory can be inverted successfully and in a numerically stable way in order
R3–O (MD, full)
2.0
R3–O (MD, 5–site) R3–H (MD, full) 1.5
g
R3–H (MD, 5–site)
1.0
0.5
0.0 0
1
2
3
4
5 r/Å
6
7
8
9
10
Fig. 4. Full atomic model and 5-site reduced model radial distribution functions between oxygen site (R3) and water oxygen/hydrogen.
Fig. 5. Three-dimensional density distribution of water oxygen sites around the full (left) and the 5-site model (right), colour-coded ranging from transparent near a density of zero over blue (dark) for bulk density up to yellow (light) and red (dark, OH-water peaks) for high density [23]. (For interpretation of the references in colour in this figure legend, the reader is referred to the Web version of this article.)
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
to yield effective potentials between a reduced set of atomic sites and the solvent. Results were checked in an unbiased manner by comparison with reference simulations. Rules have been formulated that guide the rational model design in such a way that the resulting models retain as much of the reference distribution function structure as possible. In this way, a bottom-up approach to reduced model development starting from established molecular force fields is possible without referring to iterative simulation-based compression. This proof of principle represents an important step towards our ultimate goal: the derivation of solventmediated effective solute–solute potentials between large and highly anisotropic molecules for studying their behaviour on the mesoscale. To this end the solute–solute RISM equations have to be solved which is still accompanied by severe numerical difficulties at present. The reduced model solute–solvent direct correlation functions as used in this work are needed as input to the solute–solute equation. Work into these directions is currently underway. Acknowledgements We thank the Deutsche Forschungsgemeinschaft, the Adolf-Messer-Stiftung and the Fonds der Chemischen Industrie for financial support. References [1] R.L. Henderson, Phys. Lett. 49A (1974) 197. [2] F. Mu¨ller-Plathe, ChemPhysChem 3 (2002) 754.
207
[3] F. Mu¨ller-Plathe, J. Comput. Chem. 24 (2003) 1624. [4] C.C. Liew, M. Mikami, Chem. Phys. Lett. 368 (2003) 346. [5] A.P. Lyubartsev, M. Karttunen, I. Vattulainen, A. Laaksonen, Soft. Mater. 1 (2002) 121. [6] E. Gua`rdia, J.L. Go´mez-Este´vez, J.A. Padro´, J. Chem. Phys. 86 (1987) 6438. [7] S. Amokrane, M. Bouaskarne, J. Chem. Phys. 112 (2000) 11107. [8] S. Takada, Proteins: Struct. Funct. Gen. 42 (2001) 85. [9] L. Saiz, M.L. Klein, Acc. Chem. Res. 35 (2002) 482. [10] J.-P. Hansen, I.R. McDonald, Theory of simple liquids, second edn., Academic Press, London, 1991. [11] D. Chandler, H.C. Andersen, J. Chem. Phys. 57 (1972) 1930. [12] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [13] SYBYL, Tripos Inc., St. Louis, MO. [14] D. Pearlman et al., Amber4.1, University of California, San Francisco, CA, 1995. [15] M.J. Frisch et al., GAUSSIAN 94, Gaussian, Inc., Pittsburgh, PA, 1995. [16] T.R. Forester, W. Smith, DL_POLY Molecular Dynamics Code, version 2.14, CCP5 of the EPSRC, 1995. [17] P.J. Steinbach, B.R. Brooks, J. Comput. Chem. 15 (1994) 667. [18] W.E. Reiher III, Ph.D. thesis, Harvard University, Cambridge, MA, 1985. [19] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. Di Nola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [20] P.P. Ewald, Ann. Phys. (Leipzig) 64 (1921) 253. [21] A. Kovalenko, S. Ten-no, F. Hirata, J. Comput. Chem. 20 (1999) 928. [22] S.M. Kast, K.F. Schmidt, B. Schilling, Chem. Phys. Lett. 367 (2003) 398. [23] J. Brickmann, T. Goetze, W. Heiden, G. Moeckel, S. Reiling, H. Vollhardt, C.-D. Zachmann, in: J.E. Bowie (Ed.), Data Visualisation in Molecular Science: Tools for Insight and Innovation, Addison-Wesley, Reading, 1995, p. 83; extended by M. Keil.
Design of reduced molecular models by integral equation theory Stefan M. Kast *, Wulf Hauptmann 1, Bernd Schilling Physikalische Chemie, Technische Universita¨t Darmstadt, Petersenstraße 20, 64287 Darmstadt, Germany Received 2 September 2005; in final form 24 October 2005 Available online 18 November 2005
Abstract Molecular integral equation theory is inverted with respect to the calculation of effective potentials for a coarse-grained, reduced molecular model composed of fewer sites than the reference. The methodology is applied to an aqueous solution of cyclohexanol as a model case, the quality of the effective potentials is studied by extensive molecular dynamics simulations. Practical guidelines for the choice of reduced model sites are worked out and perspectives for mesoscale modelling are discussed. 2005 Elsevier B.V. All rights reserved.
1. Introduction The development of methods for bridging length and time scale limitations in molecular simulations is currently a very active area of research. Conceptually, with the present-day computational resources a systematic coarse-graining of detailed atomic models represents a necessary and useful approach for studying systems on the mesoscale. In the soft matter sciences, particularly for problems of colloids and polymers, coarse-graining procedures are well established: based on a theorem proved by Henderson more than 30 years ago [1], equating distribution functions for a reference and a reduced model leads to a unique coarse-grained model effective potential by inversion of the distribution, typically computed from molecular simulations. The inversion process is usually performed by iterative molecular simulations [2–4], by inverse Monte-Carlo in conjunction with dissipative particle dynamics simulations [5], or by inversion of the atomic Ornstein-Zernike integral equation within the hypernetted chain (HNC) approximation [6,7]. For biomolecular problems, however, where strongly directional forces dominate the scenarios, reduced model *
Corresponding author. Fax: +49 6151 164298. E-mail address: [email protected] (S.M. Kast). 1 Present address: Technische Chemie, Technische Universita¨t Darmstadt, Petersenstraße 20, 64287 Darmstadt, Germany. 0009-2614/$ - see front matter 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2005.10.125
development is far less common. Procedures typically start heuristically (see e.g. [8,9]) and imply checking/adapting the results by explicit simulations. In the case of very large systems where the target observable cannot even be computed from reference simulations, we need another, systematic approach for effective potential development: molecular integral equation theory could provide the key to these difficulties, since it is possible to apply the uniqueness theorem [1] and to invert the potential directly without simulation effort for various reduced models that can retain enough structural information to represent strongly anisotropic interactions. Of course, the necessarily approximate correlation functions from an integral equation treatment imply errors in the computed effective potentials. As a proof of concept, we present in this Letter an integral equation-based methodology for molecular coarse-graining, i.e. reduction to fewer interaction sites than the reference system, on the example of solute–solvent interactions. We test the resulting strategy with a simple amphiphilic model system, cyclohexanol in aqueous solution by extensive molecular dynamics simulations. To the best of our knowledge no such attempt on the basis of integral equations alone has been described in the past. Since the present work represents only an intermediate, though necessary step towards a general coarse-graining approach, the next stages of development are discussed in the concluding section.
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
solution to the solute–solvent RISM–HNC equations for the reference system, and c2 from analogously requiring according to (2)
2. Theory and methods 2.1. Model reduction by integral equation theory The one-dimensional reference interaction site model (RISM) integral equation [10,11] for fluids composed of molecules with sites a and c that is used in this work reads hvv ¼ xvv cvv xvv þ xvv cvv ðqhvv Þ;
ð1Þ
where h = (hac(r)) is the matrix of the total correlation functions (the radial distribution function is g = h + 1), c = (cac(r)) are the direct correlation functions, x = (xac(r)) are the intramolecular correlation functions (for rigid molecules these are normalised Dirac delta functions constraining the site distance), q is the site density, r is the distance and the star denotes a convolution product. In contrast to the pure solvent case (superscript ÔvvÕ), for a solute system in infinite dilution (ÔuvÕ) the simpler relation qhuv ¼ xuu cuv vvv ;
203
ð2Þ
ðq^h1 ^v1 Þred ¼ ðq^h2 ^v1 Þ;
ð6Þ
where the carets denote Fourier transforms and the subscript ÔredÕ indicates that only those lines in the product matrix are taken into account that correspond to the solute–solvent site pairs of the reduced model. The structural difference between reference and reduced model is reflected by a modified intramolecular correlation matrix, such that we obtain finally from (2) and (6) ^c2 ¼ x ^ 1 ^ 1^c1 Þred . 2 ðx
ð7Þ
Inverse Fourier transform of the latter direct correlation function matrix and inserting in (5) yields a direct, noniterative and numerically stable route to the effective potentials u2 using only structural information (in the case of rigid models) and the direct correlation functions of the reference system.
holds where vvv ¼ qxvv þ q2 hvv ;
ð3Þ
the solvent susceptibility (density correlation function), can be computed beforehand and independently of any solute. These equations must be supplied with a closure relation hac ðrÞ þ 1 ¼ expðbuac ðrÞ þ hac ðrÞ cac ðrÞ þ Bac ðrÞÞ;
ð4Þ
with the pair potential uac and b = 1/kT where k is the Boltzmann constant and T is the absolute temperature; Bac are the bridge functions (for site-averaging theories like the RISM approach these functions have not the same diagrammatic meaning as in orientation-preserving theories). In the HNC approximation that is used exclusively here for both solvent–solvent and solute–solvent problems (extended RISM, RISM–HNC) B identically vanishes. Our goal in this Letter is to demonstrate a proof of concept that a coarse-graining model reduction strategy based on integral equation theory is feasible even within the simple RISM–HNC approximation under well-defined conditions. We focus exemplarily on the solute–solvent case where we intend to derive effective potentials for a reduced solute model composed of fewer sites than a reference model by requiring identity between total solute–solvent correlation functions for reduced sites and corresponding atoms of the full model. Direct inversion of the closure with respect to the potential by inserting a modified direct correlation function from inverting the RISM integral equation is numerically difficult because of the r ! 0 singularity of both potential and ln(h + 1). A more practical route is available by considering the difference between the known full potential, u1, and reduced model effective potential, u2, for the set of reduced model sites Duac ¼ u2;ac u1;ac ¼ b1 ðc2;ac c1;ac Þ;
2.2. Models and numerical methods We tested the methodology on a simple model system: a rigid cyclohexanol model as an example for an amphiphilic solute system containing both hydrophobic and hydrophilic regions, solvated in TIP3P water [12] under ambient conditions. The atom numbers and the location of the reduced sites are shown in Fig. 1, the atomic coordinates and force field parameters for the reference system are given in Table 1. The structure was modelled and optimised by Sybyl [13]. A simple Lennard-Jones/Coulomb reference potential function was applied using Lorentz–Berthelot combination rules. The Lennard-Jones parameters were taken from the Amber4.1 force field [14], partial charges resulted from fitting to the electrostatic potential derived from HF/6-31G* quantum-chemical calculations with GAUSSIAN94 [15]. We explicitly investigated three reduced models for comparison with the full 19-site system: reduced sites R1–R5 (5-site), R1, R2, R5 (3-site), and the 1-site model comprising only R2. The effective potentials derived by 1D-RISM–HNC theory are spherically symmetric and given numerically on a grid that can directly be used within the DL_POLY simulation package [16]. Except for the
ð5Þ
as derived from the closure (4) under the condition of total correlation function identity. We take c1 directly from a
Fig. 1. Atom and reduced site numbers of the cyclohexanol models. The 5-site model consists of sites R1–R5, the 3-site model of R1, R2, R5, the 1site model comprises R2.
204
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
Table 1 Structural data and potential parameters (full model: Coulomb charge q, Lennard-Jones e, r) for the cyclohexanol model ˚ ˚ ˚ Atom x/A y/A z/A q/e e/kJ mol1 C1 (R1) H2 H3 C4 H5 H6 C7 (R2) H8 O9 (R3) H10 (R4) C11 H12 H13 C14 (R5) H15 H16 C17 H18 H19
1.0530 1.5410 1.4340 0.4820 0.7210 0.9620 1.0430 0.5940 2.4630 2.9450 0.6910 1.0760 1.1770 0.8440 1.3260 1.0780 1.4040 0.9750 2.4990
1.1261 1.3741 1.8231 1.2941 2.3231 1.1171 0.2931 0.5131 0.4551 0.2851 1.1639 1.8549 1.4119 1.3309 1.1519 2.3629 0.3309 0.5549 0.4419
0.4050 0.5530 1.1710 0.2640 0.0530 1.2400 0.7800 1.7640 0.9340 0.1310 0.3810 1.1480 0.5770 0.2410 1.2170 0.0710 0.8050 1.7950 0.8740
1-site model all energy functions are long-ranged (see below) which necessitated the application of a shifted-force modification: Reduced potentials to be used in simulations were multiplied by the CHARMM function for long-range ˚. electrostatics [17] (1-(r/rc)2)2, with a cutoff distance of 12 A
0.02990 0.00747 0.01621 0.25237 0.06025 0.07063 0.45409 0.05127 0.76214 0.43158 0.25189 0.06021 0.07053 0.02924 0.00758 0.01649 0.09673 0.03880 0.01887
0.44964 0.06453 0.06453 0.44964 0.06453 0.06453 0.44964 0.06453 0.69871 0.19211 0.44964 0.06453 0.06453 0.44964 0.06453 0.06453 0.44964 0.06453 0.06453
˚ r/A 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953 3.39967 2.64953 3.00000 0.40000 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953 3.39967 2.64953 2.64953
Lennard-Jones interactions were also treated by the appropriate CHARMM shifted-force expression and truncated ˚ . Potentials were handled in the same way in beyond 12 A the simulations and in the integral equation treatment, ˚, except that small Lennard-Jones parameters (rH = 0.4 A
Fig. 2. Full atomic model and reduced model effective potentials between site R2 and water oxygen (top), R2 and water hydrogen (bottom) for various reduction stages.
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
eH = 0.1921 kJ mol1) [18] were attributed to the hydrogen sites in order to avoid numerical singularities within the RISM theory. Molecular dynamics simulations were performed in the isothermal/isobaric ensemble [19] (300 K/1 bar) with a temperature coupling constant of 0.5, and 2 ps for the barostat, using a time step of 2 fs. Long-range electrostatics was treated by Ewald summation [20] with 8 · 8· 8 recipro˚ and a charge cal space vectors, a real-space cutoff of 12 A 1 ˚ smearing parameter of 0.3 A . A single solute molecule was embedded in a pre-equilibrated water box of 2000 molecules, equilibrated for another 230 ps and finally sampled for 1.36 ns (1-site), 1.3 ns (3-site), 1.58 ns (5-site), 1.38 ns (full model). Distribution functions were computed on a ˚ spacing (radial averages) and 1 A ˚ 3 spacing grid with 0.1 A (3D density). The RISM solutions were obtained on a log˚ to a arithmic grid of 512 points ranging from 5.98 · 103 A ˚ maximum distance of 164.02 A at the same average density ˚ 3) and temperature as in the simulations. For (0.03291 A the solvent–solvent case, no phenomenological adjustment of the dielectric constant was applied. The nonlinear system of equations was solved using a variant of the Ômodified direct inversion of iterative subspaceÕ (MDIIS) method [21] for coarse optimisation of the correlation func-
205
tions and the traditional Picard iteration scheme for ultimate convergence to a threshold of max(Dc) < 107(vv) and max(Dc) < 104(uv), where Dc means the difference of the direct correlation functions between two successive iteration steps. 3. Results and discussion The resulting effective potentials are shown in Fig. 2 exemplarily for the site R2 that is preserved in all models. As can be seen, the potentials are smooth enough in any case for a direct numerical tabulation within the simulation programme. Except for the 1-site model that is dominated by the overall dipole of the solute, all potentials are apparently long-ranged as expected, but the asymptotics do not follow a simple power law. A shifted-force treatment using truncation is therefore necessary in this case. Apparently, the effective potentials can substantially change their ÔpolarityÕ, i.e. change from overall repulsive to repulsive/attractive upon switching to a different number of reduced sites. The quality of the resulting models can only be assessed by investigating the corresponding simulated distribution functions. As shown in Fig. 3 again for the case of site R2, the initial peak in the radial distribution functions near
1.4 1.2
g (R2–O)
1.0 0.8 MD, full MD, 5–site
0.6
MD, 3–site 0.4
MD, 1–site RISM–HNC
0.2 0.0 1
2
3
4
5
6
7
8
9
10
11
9
10
11
1.4 1.2
g (R2–H)
1.0 0.8 MD, full 0.6
MD, 5–site MD, 3–site MD, 1–site RISM–HNC
0.4 0.2 0.0 1
2
3
4
5
6 r/Å
7
8
Fig. 3. Full atomic model and reduced model radial distribution functions between site R2 and water oxygen (top), R2 and water hydrogen (bottom) for various reduction stages, simulated data and integral equation results.
206
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
the onset of the reference distribution (full model) is only resolved by the 5-site model that preserves the hydroxy group. The RISM–HNC solution also shows no significantly peaked structure in this range. The 3-site and the 1-site models are not able to resolve the reference structure of the distribution functions adequately. For the 5-site model however, the quality of the correspondence is very good, remaining peak shifts can be correlated with the application of the shifted-force modifications. The same statement holds for all the other sites, as exemplarily shown for the OH group oxygen g functions in Fig. 4. We can correlate the quality of the reduced models with not too strong differences between original and reduced potentials. In a number of further experiments we found that the quality of the 1-site and the 3-site models cannot be improved by using the ÔexactÕ solvent susceptibility (having almost no influence on the results) or bridge functions extracted from the simulation of the reference system (cf. [22]). In the latter case, assuming the same bridge function for corresponding sites in the reference and the reduced models is therefore not sufficient; the local molecular structure has to enter into a suitable bridge model expression if improvements beyond the HNC approximation are sought.
In summary, useful rules for the rational development of reduced molecular models by RISM–HNC inversion are: (1) strongly polar groups should be preserved in the reduced models, weakly polar moieties can safely be compressed into substantially fewer sites and (2) the quality of the resulting model hinges upon the form of the reduced potential that should not deviate too severely (in a sense of changing polarity) from the original potential. In a final check, we calculated an observable that did not enter the original, radial g function-based derivation, namely the full threedimensional density distribution of water oxygen around the entire solute. As shown in Fig. 5 (image created with MOLCAD [23]) we find an overall good quality of the agreement, particularly the solvent cavity size and form is well represented, as well as the distinct coordination structure of the hydroxy function. Though it is clear that the quality strongly depends on the choice of the location of the reduced sites, the deficiencies of the RISM–HNC approximations can be eliminated by the rules formulated above. 4. Concluding remarks We have demonstrated that RISM–HNC theory can be inverted successfully and in a numerically stable way in order
R3–O (MD, full)
2.0
R3–O (MD, 5–site) R3–H (MD, full) 1.5
g
R3–H (MD, 5–site)
1.0
0.5
0.0 0
1
2
3
4
5 r/Å
6
7
8
9
10
Fig. 4. Full atomic model and 5-site reduced model radial distribution functions between oxygen site (R3) and water oxygen/hydrogen.
Fig. 5. Three-dimensional density distribution of water oxygen sites around the full (left) and the 5-site model (right), colour-coded ranging from transparent near a density of zero over blue (dark) for bulk density up to yellow (light) and red (dark, OH-water peaks) for high density [23]. (For interpretation of the references in colour in this figure legend, the reader is referred to the Web version of this article.)
S.M. Kast et al. / Chemical Physics Letters 418 (2006) 202–207
to yield effective potentials between a reduced set of atomic sites and the solvent. Results were checked in an unbiased manner by comparison with reference simulations. Rules have been formulated that guide the rational model design in such a way that the resulting models retain as much of the reference distribution function structure as possible. In this way, a bottom-up approach to reduced model development starting from established molecular force fields is possible without referring to iterative simulation-based compression. This proof of principle represents an important step towards our ultimate goal: the derivation of solventmediated effective solute–solute potentials between large and highly anisotropic molecules for studying their behaviour on the mesoscale. To this end the solute–solute RISM equations have to be solved which is still accompanied by severe numerical difficulties at present. The reduced model solute–solvent direct correlation functions as used in this work are needed as input to the solute–solute equation. Work into these directions is currently underway. Acknowledgements We thank the Deutsche Forschungsgemeinschaft, the Adolf-Messer-Stiftung and the Fonds der Chemischen Industrie for financial support. References [1] R.L. Henderson, Phys. Lett. 49A (1974) 197. [2] F. Mu¨ller-Plathe, ChemPhysChem 3 (2002) 754.
207
[3] F. Mu¨ller-Plathe, J. Comput. Chem. 24 (2003) 1624. [4] C.C. Liew, M. Mikami, Chem. Phys. Lett. 368 (2003) 346. [5] A.P. Lyubartsev, M. Karttunen, I. Vattulainen, A. Laaksonen, Soft. Mater. 1 (2002) 121. [6] E. Gua`rdia, J.L. Go´mez-Este´vez, J.A. Padro´, J. Chem. Phys. 86 (1987) 6438. [7] S. Amokrane, M. Bouaskarne, J. Chem. Phys. 112 (2000) 11107. [8] S. Takada, Proteins: Struct. Funct. Gen. 42 (2001) 85. [9] L. Saiz, M.L. Klein, Acc. Chem. Res. 35 (2002) 482. [10] J.-P. Hansen, I.R. McDonald, Theory of simple liquids, second edn., Academic Press, London, 1991. [11] D. Chandler, H.C. Andersen, J. Chem. Phys. 57 (1972) 1930. [12] W.L. Jorgensen, J. Chandrasekhar, J.D. Madura, R.W. Impey, M.L. Klein, J. Chem. Phys. 79 (1983) 926. [13] SYBYL, Tripos Inc., St. Louis, MO. [14] D. Pearlman et al., Amber4.1, University of California, San Francisco, CA, 1995. [15] M.J. Frisch et al., GAUSSIAN 94, Gaussian, Inc., Pittsburgh, PA, 1995. [16] T.R. Forester, W. Smith, DL_POLY Molecular Dynamics Code, version 2.14, CCP5 of the EPSRC, 1995. [17] P.J. Steinbach, B.R. Brooks, J. Comput. Chem. 15 (1994) 667. [18] W.E. Reiher III, Ph.D. thesis, Harvard University, Cambridge, MA, 1985. [19] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. Di Nola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684. [20] P.P. Ewald, Ann. Phys. (Leipzig) 64 (1921) 253. [21] A. Kovalenko, S. Ten-no, F. Hirata, J. Comput. Chem. 20 (1999) 928. [22] S.M. Kast, K.F. Schmidt, B. Schilling, Chem. Phys. Lett. 367 (2003) 398. [23] J. Brickmann, T. Goetze, W. Heiden, G. Moeckel, S. Reiling, H. Vollhardt, C.-D. Zachmann, in: J.E. Bowie (Ed.), Data Visualisation in Molecular Science: Tools for Insight and Innovation, Addison-Wesley, Reading, 1995, p. 83; extended by M. Keil.