1
2 #ifdef HAVE_CONFIG_H
3 #include <scconfig.h>
4 #endif
5
6 #include <fstream>
7
8 #include <util/keyval/keyval.h>
9 #include <math/isosurf/shape.h>
10 #include <chemistry/qc/wfn/solvent.h>
11 #include <chemistry/molecule/formula.h>
12
13 #ifdef USING_NAMESPACE_STD
14 using namespace std;
15 #endif
16 using namespace sc;
17
18 static inline double
get_ki(int z)19 get_ki(int z)
20 {
21 // The ki values (used in the computation of the dispersion coefficients)
22 // for H, C, and N were taken from Vigne-Maeder and Claverie, JACS 1987, v109, pp24-28
23 // and the value for O from Huron and Claverie, J. Phys. Chem. 1974, v78, p1862
24
25 double ki;
26
27 if (z <= 0) {
28 ExEnv::errn() << "Non-positive nuclear charge encountered in computation of"
29 << " dispersion coefficient" << endl;
30 abort();
31 }
32 else if (z == 1) ki = 1.0;
33 else if (z == 6) ki = 1.0;
34 else if (z == 7) ki = 1.18;
35 else if (z == 8) ki = 1.36; // from Huron & Claverie, J.Phys.Chem v78, 1974, p1862
36 else if (z > 1 && z < 6) {
37 ki = 1.0;
38 ExEnv::out0() << "Warning: No d6 dispersion coefficient available for atomic number " <<
39 z << "; using value for carbon instead" << endl;
40 }
41 else {
42 ki = 1.18;
43 ExEnv::out0() << "Warning: No d6 dispersion coefficient available for atomic number " <<
44 z << "; using value for nitrogen instead" << endl;
45 }
46
47 return ki;
48 }
49
50 static inline double
get_d6ii(int z,double r_vdw)51 get_d6ii(int z, double r_vdw)
52 {
53 // The dispersion coefficient d6 for a pair of atoms ij can be computed
54 // from the dispersion coefficient d6ii for atom pair ii and d6jj for
55 // atom pair jj by the formula: d6 = sqrt(d6ii*d6jj).
56 // The dispersion coefficients d8 and d10 can be obtained from d6.
57 // The d6ii values given below were taken from: Vigne-Maeder and Claverie
58 // JACS 1987, v. 109, pp. 24-28.
59
60 const double a6 = 0.143; // [kcal/mol]
61 double d6ii;
62 double ki;
63
64 Ref<Units> unit = new Units("kcal/mol");
65
66 ki = get_ki(z);
67 d6ii = ki*ki*a6*pow(4*r_vdw*r_vdw,3.0); // units of (kcal mol^-1)*bohr^6
68 d6ii *= unit->to_atomic_units(); // convert to atomic units
69 return d6ii;
70 }
71
72 static inline double
get_d8ii(double d6ii,double r_vdw)73 get_d8ii(double d6ii, double r_vdw)
74 {
75 // The value of c8 was taken from Vigne-Maeder and Claverie, JACS 1987,
76 // v. 109, pp 24-28 and is here obtained in atomic units by using
77 // atomic units for d6ii and r_vdw
78
79 double d8ii;
80 const double c8 = 0.26626;
81
82 d8ii = d6ii*c8*4*pow(r_vdw,2.0);
83
84 return d8ii;
85 }
86
87 static inline double
get_d10ii(double d6ii,double r_vdw)88 get_d10ii(double d6ii, double r_vdw)
89 {
90 // The value of c10 was taken from Vigne-Maeder and Claverie, JACS 1987,
91 // v. 109, pp 24-28 and is here obtained in atomic units by using
92 // atomic units for d6ii and r_vdw
93
94 double d10ii;
95 const double c10 = 0.095467;
96
97 d10ii = d6ii*c10*16*pow(r_vdw,4.0);
98
99 return d10ii;
100 }
101
102 // For debugging compute 6, 8, and 10 contributions separately
103 static inline double
disp6_contrib(double rasnorm,double d6)104 disp6_contrib(double rasnorm, double d6)
105 {
106 double edisp6_contrib;
107
108 edisp6_contrib = d6/(3*pow(rasnorm,6.0)); // atomic units
109
110 return edisp6_contrib;
111 }
112
113 static inline double
disp8_contrib(double rasnorm,double d8)114 disp8_contrib(double rasnorm, double d8)
115 {
116 double edisp8_contrib;
117
118 edisp8_contrib = d8/(5*pow(rasnorm,8.0)); // atomic units
119
120 return edisp8_contrib;
121 }
122
123 static inline double
disp10_contrib(double rasnorm,double d10)124 disp10_contrib(double rasnorm, double d10)
125 {
126 double edisp10_contrib;
127
128 edisp10_contrib = d10/(7*pow(rasnorm,10.0)); // atomic units
129
130 return edisp10_contrib;
131 }
132
133 static inline double
disp_contrib(double rasnorm,double d6,double d8,double d10)134 disp_contrib(double rasnorm, double d6, double d8, double d10)
135 {
136 double edisp_contrib;
137
138 edisp_contrib = d6/(3*pow(rasnorm,6.0)) + d8/(5*pow(rasnorm,8.0))
139 + d10/(7*pow(rasnorm,10.0));
140
141 return edisp_contrib;
142 }
143
144 static inline double
rep_contrib(double rasnorm,double ri_vdw,double rj_vdw,double ki,double kj,double kcalpermol_to_hartree)145 rep_contrib(double rasnorm, double ri_vdw, double rj_vdw, double ki, double kj,
146 double kcalpermol_to_hartree)
147 {
148 // The expression and the parameters used for the repulsion energy
149 // were taken from Vigne-Maeder and Claverie, JACS 1987, v109, pp24-28
150 // NB: We have omitted the factor Gij
151
152 const double c = 90000; // [kcal/mol]
153 const double gamma = 12.35;
154 double erep_contrib;
155 double tmp;
156
157 tmp = gamma*rasnorm/(2.0*sqrt(ri_vdw*rj_vdw));
158
159 erep_contrib = -ki*kj*c*(1.0/tmp + 2.0/(tmp*tmp) + 2.0/(tmp*tmp*tmp))*exp(-tmp);
160 erep_contrib *= kcalpermol_to_hartree; // convert from kcal/mol to atomic units
161
162 return erep_contrib;
163 }
164
165 double
disprep()166 BEMSolvent::disprep()
167 {
168 double edisprep = 0.0;
169 double edisprep_contrib;
170 double edisp6_contrib, edisp8_contrib, edisp10_contrib; // for debugging
171 double erep_contrib;
172 double edisp6 = 0.0; // for debugging
173 double edisp8 = 0.0; // for debugging
174 double edisp10 = 0.0; // for debugging
175 double erep = 0.0;
176 double proberadius;
177 double radius;
178 double rasnorm;
179 double weight;
180 double d6, d8, d10; // dispersion coefficients
181 double d6aa, d8aa, d10aa; // dispersion coefficients for atom pair aa
182 double d6ss, d8ss, d10ss; // dispersion coefficients for atom pair ss
183 int i, iloop, isolute;
184 int natomtypes;
185 int z_solvent_atom;
186
187 Ref<Units> unit = new Units("kcal/mol");
188 double kcalpermol_to_hartree = unit->to_atomic_units();
189
190 Ref<AtomInfo> atominfo = solute_->atominfo();
191 Ref<AtomInfo> solventatominfo = solvent_->atominfo();
192 MolecularFormula formula(solvent_);
193
194 // Compute number of different atom types in solvent molecule
195 natomtypes = formula.natomtypes();
196
197 double *solute_d6ii = new double[solute_->natom()];
198 double *solute_d8ii = new double[solute_->natom()];
199 double *solute_d10ii = new double[solute_->natom()];
200 double *solute_ki = new double[solute_->natom()];
201
202 for (isolute=0; isolute<solute_->natom(); isolute++) {
203 int Z_solute = solute_->Z(isolute);
204 double radius = atominfo->vdw_radius(Z_solute);
205 solute_d6ii[isolute] = get_d6ii(Z_solute,radius);
206 solute_d8ii[isolute] = get_d8ii(solute_d6ii[isolute],radius);
207 solute_d10ii[isolute] = get_d10ii(solute_d6ii[isolute],radius);
208 solute_ki[isolute] = get_ki(Z_solute);
209 }
210
211 // Loop over atom types in solvent molecule
212 for (iloop=0; iloop<natomtypes; iloop++) {
213
214 // define the shape of the surface for current atom type
215 Ref<UnionShape> us = new UnionShape;
216 z_solvent_atom = formula.Z(iloop);
217 proberadius = solventatominfo->vdw_radius(z_solvent_atom);
218 for (i=0; i<solute_->natom(); i++) {
219 us->add_shape(new SphereShape(solute_->r(i),
220 atominfo->vdw_radius(solute_->Z(i))+proberadius));
221 }
222
223 // triangulate the surface
224 Ref<AssignedKeyVal> keyval = new AssignedKeyVal;
225 keyval->assign("volume", us.pointer());
226 keyval->assign("order", 2);
227 keyval->assign("remove_short_edges", 1);
228 keyval->assign("remove_small_triangles", 1);
229 keyval->assign("remove_slender_triangles", 1);
230 keyval->assign("short_edge_factor", 0.8);
231 keyval->assign("small_triangle_factor", 0.8);
232 keyval->assign("slender_triangle_factor", 0.8);
233 Ref<TriangulatedImplicitSurface> ts = new TriangulatedImplicitSurface(keyval.pointer());
234 ts->init();
235
236 // Debug print: check the triangulated surface
237 // if (iloop == 0) {
238 // ofstream geomviewfile("geomview.input");
239 // ts->print_geomview_format(geomviewfile);
240 // }
241
242 ExEnv::out0().setf(ios::scientific,ios::floatfield); // use scientific format
243 ExEnv::out0() << "Area of disp-rep surface generated with atom number "
244 << setw(3) << setfill(' ') << z_solvent_atom
245 << " as probe: " << setprecision(4) << ts->area()
246 << " bohr^2" << endl;
247
248 edisprep_contrib = 0.0;
249 edisp6_contrib = 0.0; // for debugging
250 edisp8_contrib = 0.0; // for debugging
251 edisp10_contrib = 0.0; // for debugging
252 erep_contrib = 0.0;
253 TriangulatedSurfaceIntegrator triint(ts.pointer());
254
255 double solvent_ki = get_ki(z_solvent_atom);
256 d6ss = get_d6ii(z_solvent_atom,proberadius);
257 d8ss = get_d8ii(d6ss, proberadius);
258 d10ss = get_d10ii(d6ss, proberadius);
259
260 // integrate the surface
261 for (triint=0; triint.update(); triint++) {
262 SCVector3 dA = triint.dA();
263 SCVector3 location = triint.current()->point();
264 weight = triint.weight();
265
266 //Loop over atoms in solute
267 for (isolute=0; isolute<solute_->natom(); isolute++) {
268
269 SCVector3 atom(solute_->r(isolute));
270 SCVector3 ras = location - atom;
271 rasnorm = ras.norm();
272 radius = atominfo->vdw_radius(solute_->Z(isolute));
273 d6aa = solute_d6ii[isolute];
274 d8aa = solute_d8ii[isolute];
275 d10aa = solute_d10ii[isolute];
276 d6 = sqrt(d6aa*d6ss);
277 d8 = sqrt(d8aa*d8ss);
278 d10 = sqrt(d10aa*d10ss);
279
280 double f = ras.dot(dA)*weight;
281 double tdisp6 = f*disp6_contrib(rasnorm,d6);
282 double tdisp8 = f*disp8_contrib(rasnorm,d8);
283 double tdisp10 = f*disp10_contrib(rasnorm,d10);
284 double trep = f*rep_contrib(rasnorm,radius,proberadius,
285 solute_ki[isolute],solvent_ki,
286 kcalpermol_to_hartree);
287 double tdisp = tdisp6+tdisp8+tdisp10;
288
289 // add in contributions to various energies; the minus sign
290 // is there to get the normal pointing into the cavity
291 edisprep_contrib -= tdisp+trep;
292 edisp6_contrib -= tdisp6;
293 edisp8_contrib -= tdisp8;
294 edisp10_contrib -= tdisp10;
295 erep_contrib -= trep;
296
297 }
298 }
299
300 edisprep += edisprep_contrib*formula.nZ(iloop);
301 edisp6 += edisp6_contrib*formula.nZ(iloop);
302 edisp8 += edisp8_contrib*formula.nZ(iloop);
303 edisp10 += edisp10_contrib*formula.nZ(iloop);
304 erep += erep_contrib*formula.nZ(iloop);
305 }
306
307 delete[] solute_d6ii;
308 delete[] solute_d8ii;
309 delete[] solute_d10ii;
310 delete[] solute_ki;
311
312 // Multiply energies by number density of solvent
313 // Print out individual energy contributions in kcal/mol
314
315 ExEnv::out0().setf(ios::scientific,ios::floatfield); // use scientific format
316 ExEnv::out0().precision(5);
317 ExEnv::out0() << "Edisp6: " << edisp6*solvent_density_*unit->from_atomic_units()
318 << " kcal/mol" << endl;
319 ExEnv::out0() << "Edisp8: " << edisp8*solvent_density_*unit->from_atomic_units()
320 << " kcal/mol" << endl;
321 ExEnv::out0() << "Edisp10: " << edisp10*solvent_density_*unit->from_atomic_units()
322 << " kcal/mol" << endl;
323
324
325 ExEnv::out0() << "Total dispersion energy: "
326 << (edisp6 + edisp8 + edisp10)*solvent_density_*unit->from_atomic_units()
327 << " kcal/mol" << endl;
328 ExEnv::out0() << "Repulsion energy: " << setw(12) << setfill(' ')
329 << erep*solvent_density_*unit->from_atomic_units() << " kcal/mol" << endl;
330
331 return edisprep*solvent_density_; // atomic units
332
333 }
334
335
336