ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/Electrostatic.cpp
Revision: 1893
Committed: Wed Jun 19 17:19:07 2013 UTC (12 years, 1 month ago) by gezelter
File size: 41580 byte(s)
Log Message:
Some performance improvements

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39 * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010).
40 * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41 */
42
43 #include <stdio.h>
44 #include <string.h>
45
46 #include <cmath>
47 #include "nonbonded/Electrostatic.hpp"
48 #include "utils/simError.h"
49 #include "types/NonBondedInteractionType.hpp"
50 #include "types/FixedChargeAdapter.hpp"
51 #include "types/FluctuatingChargeAdapter.hpp"
52 #include "types/MultipoleAdapter.hpp"
53 #include "io/Globals.hpp"
54 #include "nonbonded/SlaterIntegrals.hpp"
55 #include "utils/PhysicalConstants.hpp"
56 #include "math/erfc.hpp"
57 #include "math/SquareMatrix.hpp"
58
59 namespace OpenMD {
60
61 Electrostatic::Electrostatic(): name_("Electrostatic"), initialized_(false),
62 forceField_(NULL), info_(NULL),
63 haveCutoffRadius_(false),
64 haveDampingAlpha_(false),
65 haveDielectric_(false),
66 haveElectroSplines_(false)
67 {}
68
69 void Electrostatic::initialize() {
70
71 Globals* simParams_ = info_->getSimParams();
72
73 summationMap_["HARD"] = esm_HARD;
74 summationMap_["NONE"] = esm_HARD;
75 summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
76 summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
77 summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
78 summationMap_["TAYLOR_SHIFTED"] = esm_TAYLOR_SHIFTED;
79 summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD;
80 summationMap_["EWALD_FULL"] = esm_EWALD_FULL;
81 summationMap_["EWALD_PME"] = esm_EWALD_PME;
82 summationMap_["EWALD_SPME"] = esm_EWALD_SPME;
83 screeningMap_["DAMPED"] = DAMPED;
84 screeningMap_["UNDAMPED"] = UNDAMPED;
85
86 // these prefactors convert the multipole interactions into kcal / mol
87 // all were computed assuming distances are measured in angstroms
88 // Charge-Charge, assuming charges are measured in electrons
89 pre11_ = 332.0637778;
90 // Charge-Dipole, assuming charges are measured in electrons, and
91 // dipoles are measured in debyes
92 pre12_ = 69.13373;
93 // Dipole-Dipole, assuming dipoles are measured in Debye
94 pre22_ = 14.39325;
95 // Charge-Quadrupole, assuming charges are measured in electrons, and
96 // quadrupoles are measured in 10^-26 esu cm^2
97 // This unit is also known affectionately as an esu centibarn.
98 pre14_ = 69.13373;
99 // Dipole-Quadrupole, assuming dipoles are measured in debyes and
100 // quadrupoles in esu centibarns:
101 pre24_ = 14.39325;
102 // Quadrupole-Quadrupole, assuming esu centibarns:
103 pre44_ = 14.39325;
104
105 // conversions for the simulation box dipole moment
106 chargeToC_ = 1.60217733e-19;
107 angstromToM_ = 1.0e-10;
108 debyeToCm_ = 3.33564095198e-30;
109
110 // Default number of points for electrostatic splines
111 np_ = 100;
112
113 // variables to handle different summation methods for long-range
114 // electrostatics:
115 summationMethod_ = esm_HARD;
116 screeningMethod_ = UNDAMPED;
117 dielectric_ = 1.0;
118
119 // check the summation method:
120 if (simParams_->haveElectrostaticSummationMethod()) {
121 string myMethod = simParams_->getElectrostaticSummationMethod();
122 toUpper(myMethod);
123 map<string, ElectrostaticSummationMethod>::iterator i;
124 i = summationMap_.find(myMethod);
125 if ( i != summationMap_.end() ) {
126 summationMethod_ = (*i).second;
127 } else {
128 // throw error
129 sprintf( painCave.errMsg,
130 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
131 "\t(Input file specified %s .)\n"
132 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
133 "\t\"shifted_potential\", \"shifted_force\",\n"
134 "\t\"taylor_shifted\", or \"reaction_field\".\n",
135 myMethod.c_str() );
136 painCave.isFatal = 1;
137 simError();
138 }
139 } else {
140 // set ElectrostaticSummationMethod to the cutoffMethod:
141 if (simParams_->haveCutoffMethod()){
142 string myMethod = simParams_->getCutoffMethod();
143 toUpper(myMethod);
144 map<string, ElectrostaticSummationMethod>::iterator i;
145 i = summationMap_.find(myMethod);
146 if ( i != summationMap_.end() ) {
147 summationMethod_ = (*i).second;
148 }
149 }
150 }
151
152 if (summationMethod_ == esm_REACTION_FIELD) {
153 if (!simParams_->haveDielectric()) {
154 // throw warning
155 sprintf( painCave.errMsg,
156 "SimInfo warning: dielectric was not specified in the input file\n\tfor the reaction field correction method.\n"
157 "\tA default value of %f will be used for the dielectric.\n", dielectric_);
158 painCave.isFatal = 0;
159 painCave.severity = OPENMD_INFO;
160 simError();
161 } else {
162 dielectric_ = simParams_->getDielectric();
163 }
164 haveDielectric_ = true;
165 }
166
167 if (simParams_->haveElectrostaticScreeningMethod()) {
168 string myScreen = simParams_->getElectrostaticScreeningMethod();
169 toUpper(myScreen);
170 map<string, ElectrostaticScreeningMethod>::iterator i;
171 i = screeningMap_.find(myScreen);
172 if ( i != screeningMap_.end()) {
173 screeningMethod_ = (*i).second;
174 } else {
175 sprintf( painCave.errMsg,
176 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
177 "\t(Input file specified %s .)\n"
178 "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
179 "or \"damped\".\n", myScreen.c_str() );
180 painCave.isFatal = 1;
181 simError();
182 }
183 }
184
185 // check to make sure a cutoff value has been set:
186 if (!haveCutoffRadius_) {
187 sprintf( painCave.errMsg, "Electrostatic::initialize has no Default "
188 "Cutoff value!\n");
189 painCave.severity = OPENMD_ERROR;
190 painCave.isFatal = 1;
191 simError();
192 }
193
194 if (screeningMethod_ == DAMPED) {
195 if (!simParams_->haveDampingAlpha()) {
196 // first set a cutoff dependent alpha value
197 // we assume alpha depends linearly with rcut from 0 to 20.5 ang
198 dampingAlpha_ = 0.425 - cutoffRadius_* 0.02;
199 if (dampingAlpha_ < 0.0) dampingAlpha_ = 0.0;
200
201 // throw warning
202 sprintf( painCave.errMsg,
203 "Electrostatic::initialize: dampingAlpha was not specified in the\n"
204 "\tinput file. A default value of %f (1/ang) will be used for the\n"
205 "\tcutoff of %f (ang).\n",
206 dampingAlpha_, cutoffRadius_);
207 painCave.severity = OPENMD_INFO;
208 painCave.isFatal = 0;
209 simError();
210 } else {
211 dampingAlpha_ = simParams_->getDampingAlpha();
212 }
213 haveDampingAlpha_ = true;
214 }
215
216 // find all of the Electrostatic atom Types:
217 ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
218 ForceField::AtomTypeContainer::MapTypeIterator i;
219 AtomType* at;
220
221 for (at = atomTypes->beginType(i); at != NULL;
222 at = atomTypes->nextType(i)) {
223
224 if (at->isElectrostatic())
225 addType(at);
226 }
227
228 if (summationMethod_ == esm_REACTION_FIELD) {
229 preRF_ = (dielectric_ - 1.0) /
230 ((2.0 * dielectric_ + 1.0) * pow(cutoffRadius_,3) );
231 }
232
233 RealType b0c, b1c, b2c, b3c, b4c, b5c;
234 RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
235 RealType a2, expTerm, invArootPi;
236
237 RealType r = cutoffRadius_;
238 RealType r2 = r * r;
239 RealType ric = 1.0 / r;
240 RealType ric2 = ric * ric;
241
242 if (screeningMethod_ == DAMPED) {
243 a2 = dampingAlpha_ * dampingAlpha_;
244 invArootPi = 1.0 / (dampingAlpha_ * sqrt(M_PI));
245 expTerm = exp(-a2 * r2);
246 // values of Smith's B_l functions at the cutoff radius:
247 b0c = erfc(dampingAlpha_ * r) / r;
248 b1c = ( b0c + 2.0*a2 * expTerm * invArootPi) / r2;
249 b2c = (3.0 * b1c + pow(2.0*a2, 2) * expTerm * invArootPi) / r2;
250 b3c = (5.0 * b2c + pow(2.0*a2, 3) * expTerm * invArootPi) / r2;
251 b4c = (7.0 * b3c + pow(2.0*a2, 4) * expTerm * invArootPi) / r2;
252 b5c = (9.0 * b4c + pow(2.0*a2, 5) * expTerm * invArootPi) / r2;
253 selfMult_ = b0c + a2 * invArootPi;
254 } else {
255 a2 = 0.0;
256 b0c = 1.0 / r;
257 b1c = ( b0c) / r2;
258 b2c = (3.0 * b1c) / r2;
259 b3c = (5.0 * b2c) / r2;
260 b4c = (7.0 * b3c) / r2;
261 b5c = (9.0 * b4c) / r2;
262 selfMult_ = b0c;
263 }
264
265 // higher derivatives of B_0 at the cutoff radius:
266 db0c_1 = -r * b1c;
267 db0c_2 = -b1c + r2 * b2c;
268 db0c_3 = 3.0*r*b2c - r2*r*b3c;
269 db0c_4 = 3.0*b2c - 6.0*r2*b3c + r2*r2*b4c;
270 db0c_5 = -15.0*r*b3c + 10.0*r2*r*b4c - r2*r2*r*b5c;
271
272
273 // working variables for the splines:
274 RealType ri, ri2;
275 RealType b0, b1, b2, b3, b4, b5;
276 RealType db0_1, db0_2, db0_3, db0_4, db0_5;
277 RealType f, fc, f0;
278 RealType g, gc, g0, g1, g2, g3, g4;
279 RealType h, hc, h1, h2, h3, h4;
280 RealType s, sc, s2, s3, s4;
281 RealType t, tc, t3, t4;
282 RealType u, uc, u4;
283
284 // working variables for Taylor expansion:
285 RealType rmRc, rmRc2, rmRc3, rmRc4;
286
287 // Approximate using splines using a maximum of 0.1 Angstroms
288 // between points.
289 int nptest = int((cutoffRadius_ + 2.0) / 0.1);
290 np_ = (np_ > nptest) ? np_ : nptest;
291
292 // Add a 2 angstrom safety window to deal with cutoffGroups that
293 // have charged atoms longer than the cutoffRadius away from each
294 // other. Splining is almost certainly the best choice here.
295 // Direct calls to erfc would be preferrable if it is a very fast
296 // implementation.
297
298 RealType dx = (cutoffRadius_ + 2.0) / RealType(np_);
299
300 // Storage vectors for the computed functions
301 vector<RealType> rv;
302 vector<RealType> v01v;
303 vector<RealType> v11v;
304 vector<RealType> v21v, v22v;
305 vector<RealType> v31v, v32v;
306 vector<RealType> v41v, v42v, v43v;
307
308 /*
309 vector<RealType> dv01v;
310 vector<RealType> dv11v;
311 vector<RealType> dv21v, dv22v;
312 vector<RealType> dv31v, dv32v;
313 vector<RealType> dv41v, dv42v, dv43v;
314 */
315
316 for (int i = 1; i < np_ + 1; i++) {
317 r = RealType(i) * dx;
318 rv.push_back(r);
319
320 ri = 1.0 / r;
321 ri2 = ri * ri;
322
323 r2 = r * r;
324 expTerm = exp(-a2 * r2);
325
326 // Taylor expansion factors (no need for factorials this way):
327 rmRc = r - cutoffRadius_;
328 rmRc2 = rmRc * rmRc / 2.0;
329 rmRc3 = rmRc2 * rmRc / 3.0;
330 rmRc4 = rmRc3 * rmRc / 4.0;
331
332 // values of Smith's B_l functions at r:
333 if (screeningMethod_ == DAMPED) {
334 b0 = erfc(dampingAlpha_ * r) * ri;
335 b1 = ( b0 + 2.0*a2 * expTerm * invArootPi) * ri2;
336 b2 = (3.0 * b1 + pow(2.0*a2, 2) * expTerm * invArootPi) * ri2;
337 b3 = (5.0 * b2 + pow(2.0*a2, 3) * expTerm * invArootPi) * ri2;
338 b4 = (7.0 * b3 + pow(2.0*a2, 4) * expTerm * invArootPi) * ri2;
339 b5 = (9.0 * b4 + pow(2.0*a2, 5) * expTerm * invArootPi) * ri2;
340 } else {
341 b0 = ri;
342 b1 = ( b0) * ri2;
343 b2 = (3.0 * b1) * ri2;
344 b3 = (5.0 * b2) * ri2;
345 b4 = (7.0 * b3) * ri2;
346 b5 = (9.0 * b4) * ri2;
347 }
348
349 // higher derivatives of B_0 at r:
350 db0_1 = -r * b1;
351 db0_2 = -b1 + r2 * b2;
352 db0_3 = 3.0*r*b2 - r2*r*b3;
353 db0_4 = 3.0*b2 - 6.0*r2*b3 + r2*r2*b4;
354 db0_5 = -15.0*r*b3 + 10.0*r2*r*b4 - r2*r2*r*b5;
355
356 f = b0;
357 fc = b0c;
358 f0 = f - fc - rmRc*db0c_1;
359
360 g = db0_1;
361 gc = db0c_1;
362 g0 = g - gc;
363 g1 = g0 - rmRc *db0c_2;
364 g2 = g1 - rmRc2*db0c_3;
365 g3 = g2 - rmRc3*db0c_4;
366 g4 = g3 - rmRc4*db0c_5;
367
368 h = db0_2;
369 hc = db0c_2;
370 h1 = h - hc;
371 h2 = h1 - rmRc *db0c_3;
372 h3 = h2 - rmRc2*db0c_4;
373 h4 = h3 - rmRc3*db0c_5;
374
375 s = db0_3;
376 sc = db0c_3;
377 s2 = s - sc;
378 s3 = s2 - rmRc *db0c_4;
379 s4 = s3 - rmRc2*db0c_5;
380
381 t = db0_4;
382 tc = db0c_4;
383 t3 = t - tc;
384 t4 = t3 - rmRc *db0c_5;
385
386 u = db0_5;
387 uc = db0c_5;
388 u4 = u - uc;
389
390 // in what follows below, the various v functions are used for
391 // potentials and torques, while the w functions show up in the
392 // forces.
393
394 switch (summationMethod_) {
395 case esm_SHIFTED_FORCE:
396
397 v01 = f - fc - rmRc*gc;
398 v11 = g - gc - rmRc*hc;
399 v21 = g*ri - gc*ric - rmRc*(hc - gc*ric)*ric;
400 v22 = h - g*ri - (hc - gc*ric) - rmRc*(sc - (hc - gc*ric)*ric);
401 v31 = (h-g*ri)*ri - (hc-g*ric)*ric - rmRc*(sc-2.0*(hc-gc*ric)*ric)*ric;
402 v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric)
403 - rmRc*(tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
404 v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2
405 - rmRc*(sc - 3.0*(hc-gc*ric)*ric)*ric2;
406 v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric
407 - rmRc*(tc - (4.0*sc - 9.0*(hc - gc*ric)*ric)*ric)*ric;
408
409 v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
410 - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric)
411 - rmRc*(uc-3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
412
413 dv01 = g - gc;
414 dv11 = h - hc;
415 dv21 = (h - g*ri)*ri - (hc - gc*ric)*ric;
416 dv22 = (s - (h - g*ri)*ri) - (sc - (hc - gc*ric)*ric);
417 dv31 = (s - 2.0*(h-g*ri)*ri)*ri - (sc - 2.0*(hc-gc*ric)*ric)*ric;
418 dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri)
419 - (tc - 3.0*(sc-2.0*(hc-gc*ric)*ric)*ric);
420 dv41 = (s - 3.0*(h - g*ri)*ri)*ri2 - (sc - 3.0*(hc - gc*ric)*ric)*ric2;
421 dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri
422 - (tc - (4.0*sc - 9.0*(hc-gc*ric)*ric)*ric)*ric;
423 dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri)
424 - (uc - 3.0*(2.0*tc - (7.0*sc - 15.0*(hc - gc*ric)*ric)*ric)*ric);
425
426 break;
427
428 case esm_TAYLOR_SHIFTED:
429
430 v01 = f0;
431 v11 = g1;
432 v21 = g2 * ri;
433 v22 = h2 - v21;
434 v31 = (h3 - g3 * ri) * ri;
435 v32 = s3 - 3.0*v31;
436 v41 = (h4 - g4 * ri) * ri2;
437 v42 = s4 * ri - 3.0*v41;
438 v43 = t4 - 6.0*v42 - 3.0*v41;
439
440 dv01 = g0;
441 dv11 = h1;
442 dv21 = (h2 - g2*ri)*ri;
443 dv22 = (s2 - (h2 - g2*ri)*ri);
444 dv31 = (s3 - 2.0*(h3-g3*ri)*ri)*ri;
445 dv32 = (t3 - 3.0*(s3-2.0*(h3-g3*ri)*ri)*ri);
446 dv41 = (s4 - 3.0*(h4 - g4*ri)*ri)*ri2;
447 dv42 = (t4 - (4.0*s4 - 9.0*(h4-g4*ri)*ri)*ri)*ri;
448 dv43 = (u4 - 3.0*(2.0*t4 - (7.0*s4 - 15.0*(h4 - g4*ri)*ri)*ri)*ri);
449
450 break;
451
452 case esm_SHIFTED_POTENTIAL:
453
454 v01 = f - fc;
455 v11 = g - gc;
456 v21 = g*ri - gc*ric;
457 v22 = h - g*ri - (hc - gc*ric);
458 v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
459 v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
460 v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
461 v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
462 v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
463 - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
464
465 dv01 = g;
466 dv11 = h;
467 dv21 = (h - g*ri)*ri;
468 dv22 = (s - (h - g*ri)*ri);
469 dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
470 dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
471 dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
472 dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
473 dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
474
475 break;
476
477 case esm_SWITCHING_FUNCTION:
478 case esm_HARD:
479
480 v01 = f;
481 v11 = g;
482 v21 = g*ri;
483 v22 = h - g*ri;
484 v31 = (h-g*ri)*ri;
485 v32 = (s - 3.0*(h-g*ri)*ri);
486 v41 = (h - g*ri)*ri2;
487 v42 = (s-3.0*(h-g*ri)*ri)*ri;
488 v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri);
489
490 dv01 = g;
491 dv11 = h;
492 dv21 = (h - g*ri)*ri;
493 dv22 = (s - (h - g*ri)*ri);
494 dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
495 dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
496 dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
497 dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
498 dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
499
500 break;
501
502 case esm_REACTION_FIELD:
503
504 // following DL_POLY's lead for shifting the image charge potential:
505 f = b0 + preRF_ * r2;
506 fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_;
507
508 g = db0_1 + preRF_ * 2.0 * r;
509 gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_;
510
511 h = db0_2 + preRF_ * 2.0;
512 hc = db0c_2 + preRF_ * 2.0;
513
514 v01 = f - fc;
515 v11 = g - gc;
516 v21 = g*ri - gc*ric;
517 v22 = h - g*ri - (hc - gc*ric);
518 v31 = (h-g*ri)*ri - (hc-g*ric)*ric;
519 v32 = (s - 3.0*(h-g*ri)*ri) - (sc - 3.0*(hc-gc*ric)*ric);
520 v41 = (h - g*ri)*ri2 - (hc - gc*ric)*ric2;
521 v42 = (s-3.0*(h-g*ri)*ri)*ri - (sc-3.0*(hc-gc*ric)*ric)*ric;
522 v43 = (t - 3.0*(2.0*s - 5.0*(h - g*ri)*ri)*ri)
523 - (tc - 3.0*(2.0*sc - 5.0*(hc - gc*ric)*ric)*ric);
524
525 dv01 = g;
526 dv11 = h;
527 dv21 = (h - g*ri)*ri;
528 dv22 = (s - (h - g*ri)*ri);
529 dv31 = (s - 2.0*(h-g*ri)*ri)*ri;
530 dv32 = (t - 3.0*(s-2.0*(h-g*ri)*ri)*ri);
531 dv41 = (s - 3.0*(h - g*ri)*ri)*ri2;
532 dv42 = (t - (4.0*s - 9.0*(h-g*ri)*ri)*ri)*ri;
533 dv43 = (u - 3.0*(2.0*t - (7.0*s - 15.0*(h - g*ri)*ri)*ri)*ri);
534
535 break;
536
537 case esm_EWALD_FULL:
538 case esm_EWALD_PME:
539 case esm_EWALD_SPME:
540 default :
541 map<string, ElectrostaticSummationMethod>::iterator i;
542 std::string meth;
543 for (i = summationMap_.begin(); i != summationMap_.end(); ++i) {
544 if ((*i).second == summationMethod_) meth = (*i).first;
545 }
546 sprintf( painCave.errMsg,
547 "Electrostatic::initialize: electrostaticSummationMethod %s \n"
548 "\thas not been implemented yet. Please select one of:\n"
549 "\t\"hard\", \"shifted_potential\", or \"shifted_force\"\n",
550 meth.c_str() );
551 painCave.isFatal = 1;
552 simError();
553 break;
554 }
555
556 // Add these computed values to the storage vectors for spline creation:
557 v01v.push_back(v01);
558 v11v.push_back(v11);
559 v21v.push_back(v21);
560 v22v.push_back(v22);
561 v31v.push_back(v31);
562 v32v.push_back(v32);
563 v41v.push_back(v41);
564 v42v.push_back(v42);
565 v43v.push_back(v43);
566 /*
567 dv01v.push_back(dv01);
568 dv11v.push_back(dv11);
569 dv21v.push_back(dv21);
570 dv22v.push_back(dv22);
571 dv31v.push_back(dv31);
572 dv32v.push_back(dv32);
573 dv41v.push_back(dv41);
574 dv42v.push_back(dv42);
575 dv43v.push_back(dv43);
576 */
577 }
578
579 // construct the spline structures and fill them with the values we've
580 // computed:
581
582 v01s = new CubicSpline();
583 v01s->addPoints(rv, v01v);
584 v11s = new CubicSpline();
585 v11s->addPoints(rv, v11v);
586 v21s = new CubicSpline();
587 v21s->addPoints(rv, v21v);
588 v22s = new CubicSpline();
589 v22s->addPoints(rv, v22v);
590 v31s = new CubicSpline();
591 v31s->addPoints(rv, v31v);
592 v32s = new CubicSpline();
593 v32s->addPoints(rv, v32v);
594 v41s = new CubicSpline();
595 v41s->addPoints(rv, v41v);
596 v42s = new CubicSpline();
597 v42s->addPoints(rv, v42v);
598 v43s = new CubicSpline();
599 v43s->addPoints(rv, v43v);
600
601 /*
602 dv01s = new CubicSpline();
603 dv01s->addPoints(rv, dv01v);
604 dv11s = new CubicSpline();
605 dv11s->addPoints(rv, dv11v);
606 dv21s = new CubicSpline();
607 dv21s->addPoints(rv, dv21v);
608 dv22s = new CubicSpline();
609 dv22s->addPoints(rv, dv22v);
610 dv31s = new CubicSpline();
611 dv31s->addPoints(rv, dv31v);
612 dv32s = new CubicSpline();
613 dv32s->addPoints(rv, dv32v);
614 dv41s = new CubicSpline();
615 dv41s->addPoints(rv, dv41v);
616 dv42s = new CubicSpline();
617 dv42s->addPoints(rv, dv42v);
618 dv43s = new CubicSpline();
619 dv43s->addPoints(rv, dv43v);
620 */
621
622 haveElectroSplines_ = true;
623
624 initialized_ = true;
625 }
626
627 void Electrostatic::addType(AtomType* atomType){
628
629 ElectrostaticAtomData electrostaticAtomData;
630 electrostaticAtomData.is_Charge = false;
631 electrostaticAtomData.is_Dipole = false;
632 electrostaticAtomData.is_Quadrupole = false;
633 electrostaticAtomData.is_Fluctuating = false;
634
635 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
636
637 if (fca.isFixedCharge()) {
638 electrostaticAtomData.is_Charge = true;
639 electrostaticAtomData.fixedCharge = fca.getCharge();
640 }
641
642 MultipoleAdapter ma = MultipoleAdapter(atomType);
643 if (ma.isMultipole()) {
644 if (ma.isDipole()) {
645 electrostaticAtomData.is_Dipole = true;
646 electrostaticAtomData.dipole = ma.getDipole();
647 }
648 if (ma.isQuadrupole()) {
649 electrostaticAtomData.is_Quadrupole = true;
650 electrostaticAtomData.quadrupole = ma.getQuadrupole();
651 }
652 }
653
654 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
655
656 if (fqa.isFluctuatingCharge()) {
657 electrostaticAtomData.is_Fluctuating = true;
658 electrostaticAtomData.electronegativity = fqa.getElectronegativity();
659 electrostaticAtomData.hardness = fqa.getHardness();
660 electrostaticAtomData.slaterN = fqa.getSlaterN();
661 electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
662 }
663
664 pair<map<int,AtomType*>::iterator,bool> ret;
665 ret = ElectrostaticList.insert( pair<int,AtomType*>(atomType->getIdent(),
666 atomType) );
667 if (ret.second == false) {
668 sprintf( painCave.errMsg,
669 "Electrostatic already had a previous entry with ident %d\n",
670 atomType->getIdent() );
671 painCave.severity = OPENMD_INFO;
672 painCave.isFatal = 0;
673 simError();
674 }
675
676 ElectrostaticMap[atomType] = electrostaticAtomData;
677
678 // Now, iterate over all known types and add to the mixing map:
679
680 map<AtomType*, ElectrostaticAtomData>::iterator it;
681 for( it = ElectrostaticMap.begin(); it != ElectrostaticMap.end(); ++it) {
682 AtomType* atype2 = (*it).first;
683 ElectrostaticAtomData eaData2 = (*it).second;
684 if (eaData2.is_Fluctuating && electrostaticAtomData.is_Fluctuating) {
685
686 RealType a = electrostaticAtomData.slaterZeta;
687 RealType b = eaData2.slaterZeta;
688 int m = electrostaticAtomData.slaterN;
689 int n = eaData2.slaterN;
690
691 // Create the spline of the coulombic integral for s-type
692 // Slater orbitals. Add a 2 angstrom safety window to deal
693 // with cutoffGroups that have charged atoms longer than the
694 // cutoffRadius away from each other.
695
696 RealType rval;
697 RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
698 vector<RealType> rvals;
699 vector<RealType> Jvals;
700 // don't start at i = 0, as rval = 0 is undefined for the
701 // slater overlap integrals.
702 for (int i = 1; i < np_+1; i++) {
703 rval = RealType(i) * dr;
704 rvals.push_back(rval);
705 Jvals.push_back(sSTOCoulInt( a, b, m, n, rval *
706 PhysicalConstants::angstromToBohr ) *
707 PhysicalConstants::hartreeToKcal );
708 }
709
710 CubicSpline* J = new CubicSpline();
711 J->addPoints(rvals, Jvals);
712
713 pair<AtomType*, AtomType*> key1, key2;
714 key1 = make_pair(atomType, atype2);
715 key2 = make_pair(atype2, atomType);
716
717 Jij[key1] = J;
718 Jij[key2] = J;
719 }
720 }
721
722 return;
723 }
724
725 void Electrostatic::setCutoffRadius( RealType rCut ) {
726 cutoffRadius_ = rCut;
727 haveCutoffRadius_ = true;
728 }
729
730 void Electrostatic::setElectrostaticSummationMethod( ElectrostaticSummationMethod esm ) {
731 summationMethod_ = esm;
732 }
733 void Electrostatic::setElectrostaticScreeningMethod( ElectrostaticScreeningMethod sm ) {
734 screeningMethod_ = sm;
735 }
736 void Electrostatic::setDampingAlpha( RealType alpha ) {
737 dampingAlpha_ = alpha;
738 haveDampingAlpha_ = true;
739 }
740 void Electrostatic::setReactionFieldDielectric( RealType dielectric ){
741 dielectric_ = dielectric;
742 haveDielectric_ = true;
743 }
744
745 void Electrostatic::calcForce(InteractionData &idat) {
746
747 if (!initialized_) initialize();
748
749 data1 = ElectrostaticMap[idat.atypes.first];
750 data2 = ElectrostaticMap[idat.atypes.second];
751
752 U = 0.0; // Potential
753 F.zero(); // Force
754 Ta.zero(); // Torque on site a
755 Tb.zero(); // Torque on site b
756 Ea.zero(); // Electric field at site a
757 Eb.zero(); // Electric field at site b
758 dUdCa = 0.0; // fluctuating charge force at site a
759 dUdCb = 0.0; // fluctuating charge force at site a
760
761 // Indirect interactions mediated by the reaction field.
762 indirect_Pot = 0.0; // Potential
763 indirect_F.zero(); // Force
764 indirect_Ta.zero(); // Torque on site a
765 indirect_Tb.zero(); // Torque on site b
766
767 // Excluded potential that is still computed for fluctuating charges
768 excluded_Pot= 0.0;
769
770
771 // some variables we'll need independent of electrostatic type:
772
773 ri = 1.0 / *(idat.rij);
774 rhat = *(idat.d) * ri;
775
776 // logicals
777
778 a_is_Charge = data1.is_Charge;
779 a_is_Dipole = data1.is_Dipole;
780 a_is_Quadrupole = data1.is_Quadrupole;
781 a_is_Fluctuating = data1.is_Fluctuating;
782
783 b_is_Charge = data2.is_Charge;
784 b_is_Dipole = data2.is_Dipole;
785 b_is_Quadrupole = data2.is_Quadrupole;
786 b_is_Fluctuating = data2.is_Fluctuating;
787
788 // Obtain all of the required radial function values from the
789 // spline structures:
790
791 // needed for fields (and forces):
792 if (a_is_Charge || b_is_Charge) {
793 v01s->getValueAndDerivativeAt( *(idat.rij), v01, dv01);
794 }
795 if (a_is_Dipole || b_is_Dipole) {
796 v11s->getValueAndDerivativeAt( *(idat.rij), v11, dv11);
797 v11or = ri * v11;
798 }
799 if (a_is_Quadrupole || b_is_Quadrupole || (a_is_Dipole && b_is_Dipole)) {
800 v21s->getValueAndDerivativeAt( *(idat.rij), v21, dv21);
801 v22s->getValueAndDerivativeAt( *(idat.rij), v22, dv22);
802 v22or = ri * v22;
803 }
804
805 // needed for potentials (and forces and torques):
806 if ((a_is_Dipole && b_is_Quadrupole) ||
807 (b_is_Dipole && a_is_Quadrupole)) {
808 v31s->getValueAndDerivativeAt( *(idat.rij), v31, dv31);
809 v32s->getValueAndDerivativeAt( *(idat.rij), v32, dv32);
810 v31or = v31 * ri;
811 v32or = v32 * ri;
812 }
813 if (a_is_Quadrupole && b_is_Quadrupole) {
814 v41s->getValueAndDerivativeAt( *(idat.rij), v41, dv41);
815 v42s->getValueAndDerivativeAt( *(idat.rij), v42, dv42);
816 v43s->getValueAndDerivativeAt( *(idat.rij), v43, dv43);
817 v42or = v42 * ri;
818 v43or = v43 * ri;
819 }
820
821 // calculate the single-site contributions (fields, etc).
822
823 if (a_is_Charge) {
824 C_a = data1.fixedCharge;
825
826 if (a_is_Fluctuating) {
827 C_a += *(idat.flucQ1);
828 }
829
830 if (idat.excluded) {
831 *(idat.skippedCharge2) += C_a;
832 } else {
833 // only do the field if we're not excluded:
834 Eb -= C_a * pre11_ * dv01 * rhat;
835 }
836 }
837
838 if (a_is_Dipole) {
839 D_a = *(idat.dipole1);
840 rdDa = dot(rhat, D_a);
841 rxDa = cross(rhat, D_a);
842 if (!idat.excluded)
843 Eb -= pre12_ * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
844 }
845
846 if (a_is_Quadrupole) {
847 Q_a = *(idat.quadrupole1);
848 trQa = Q_a.trace();
849 Qar = Q_a * rhat;
850 rQa = rhat * Q_a;
851 rdQar = dot(rhat, Qar);
852 rxQar = cross(rhat, Qar);
853 if (!idat.excluded)
854 Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or
855 + rdQar * rhat * (dv22 - 2.0*v22or));
856 }
857
858 if (b_is_Charge) {
859 C_b = data2.fixedCharge;
860
861 if (b_is_Fluctuating)
862 C_b += *(idat.flucQ2);
863
864 if (idat.excluded) {
865 *(idat.skippedCharge1) += C_b;
866 } else {
867 // only do the field if we're not excluded:
868 Ea += C_b * pre11_ * dv01 * rhat;
869 }
870 }
871
872 if (b_is_Dipole) {
873 D_b = *(idat.dipole2);
874 rdDb = dot(rhat, D_b);
875 rxDb = cross(rhat, D_b);
876 if (!idat.excluded)
877 Ea += pre12_ * ((dv11-v11or) * rdDb * rhat + v11or * D_b);
878 }
879
880 if (b_is_Quadrupole) {
881 Q_b = *(idat.quadrupole2);
882 trQb = Q_b.trace();
883 Qbr = Q_b * rhat;
884 rQb = rhat * Q_b;
885 rdQbr = dot(rhat, Qbr);
886 rxQbr = cross(rhat, Qbr);
887 if (!idat.excluded)
888 Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or
889 + rdQbr * rhat * (dv22 - 2.0*v22or));
890 }
891
892 if ((a_is_Fluctuating || b_is_Fluctuating) && idat.excluded) {
893 J = Jij[idat.atypes];
894 }
895
896 if (a_is_Charge) {
897
898 if (b_is_Charge) {
899 pref = pre11_ * *(idat.electroMult);
900 U += C_a * C_b * pref * v01;
901 F += C_a * C_b * pref * dv01 * rhat;
902
903 // If this is an excluded pair, there are still indirect
904 // interactions via the reaction field we must worry about:
905
906 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
907 rfContrib = preRF_ * pref * C_a * C_b * *(idat.r2);
908 indirect_Pot += rfContrib;
909 indirect_F += rfContrib * 2.0 * ri * rhat;
910 }
911
912 // Fluctuating charge forces are handled via Coulomb integrals
913 // for excluded pairs (i.e. those connected via bonds) and
914 // with the standard charge-charge interaction otherwise.
915
916 if (idat.excluded) {
917 if (a_is_Fluctuating || b_is_Fluctuating) {
918 coulInt = J->getValueAt( *(idat.rij) );
919 if (a_is_Fluctuating) dUdCa += coulInt * C_b;
920 if (b_is_Fluctuating) dUdCb += coulInt * C_a;
921 excluded_Pot += C_a * C_b * coulInt;
922 }
923 } else {
924 if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
925 if (a_is_Fluctuating) dUdCb += C_a * pref * v01;
926 }
927 }
928
929 if (b_is_Dipole) {
930 pref = pre12_ * *(idat.electroMult);
931 U += C_a * pref * v11 * rdDb;
932 F += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * D_b);
933 Tb += C_a * pref * v11 * rxDb;
934
935 if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb;
936
937 // Even if we excluded this pair from direct interactions, we
938 // still have the reaction-field-mediated charge-dipole
939 // interaction:
940
941 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
942 rfContrib = C_a * pref * preRF_ * 2.0 * *(idat.rij);
943 indirect_Pot += rfContrib * rdDb;
944 indirect_F += rfContrib * D_b / (*idat.rij);
945 indirect_Tb += C_a * pref * preRF_ * rxDb;
946 }
947 }
948
949 if (b_is_Quadrupole) {
950 pref = pre14_ * *(idat.electroMult);
951 U += C_a * pref * (v21 * trQb + v22 * rdQbr);
952 F += C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or);
953 F += C_a * pref * rdQbr * rhat * (dv22 - 2.0*v22or);
954 Tb += C_a * pref * 2.0 * rxQbr * v22;
955
956 if (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr);
957 }
958 }
959
960 if (a_is_Dipole) {
961
962 if (b_is_Charge) {
963 pref = pre12_ * *(idat.electroMult);
964
965 U -= C_b * pref * v11 * rdDa;
966 F -= C_b * pref * ((dv11-v11or) * rdDa * rhat + v11or * D_a);
967 Ta -= C_b * pref * v11 * rxDa;
968
969 if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa;
970
971 // Even if we excluded this pair from direct interactions,
972 // we still have the reaction-field-mediated charge-dipole
973 // interaction:
974 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
975 rfContrib = C_b * pref * preRF_ * 2.0 * *(idat.rij);
976 indirect_Pot -= rfContrib * rdDa;
977 indirect_F -= rfContrib * D_a / (*idat.rij);
978 indirect_Ta -= C_b * pref * preRF_ * rxDa;
979 }
980 }
981
982 if (b_is_Dipole) {
983 pref = pre22_ * *(idat.electroMult);
984 DadDb = dot(D_a, D_b);
985 DaxDb = cross(D_a, D_b);
986
987 U -= pref * (DadDb * v21 + rdDa * rdDb * v22);
988 F -= pref * (dv21 * DadDb * rhat + v22or * (rdDb * D_a + rdDa * D_b));
989 F -= pref * (rdDa * rdDb) * (dv22 - 2.0*v22or) * rhat;
990 Ta += pref * ( v21 * DaxDb - v22 * rdDb * rxDa);
991 Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
992
993 // Even if we excluded this pair from direct interactions, we
994 // still have the reaction-field-mediated dipole-dipole
995 // interaction:
996 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
997 rfContrib = -pref * preRF_ * 2.0;
998 indirect_Pot += rfContrib * DadDb;
999 indirect_Ta += rfContrib * DaxDb;
1000 indirect_Tb -= rfContrib * DaxDb;
1001 }
1002 }
1003
1004 if (b_is_Quadrupole) {
1005 pref = pre24_ * *(idat.electroMult);
1006 DadQb = D_a * Q_b;
1007 DadQbr = dot(D_a, Qbr);
1008 DaxQbr = cross(D_a, Qbr);
1009
1010 U -= pref * ((trQb*rdDa + 2.0*DadQbr)*v31 + rdDa*rdQbr*v32);
1011 F -= pref * (trQb*D_a + 2.0*DadQb) * v31or;
1012 F -= pref * (trQb*rdDa + 2.0*DadQbr) * (dv31-v31or) * rhat;
1013 F -= pref * (D_a*rdQbr + 2.0*rdDa*rQb) * v32or;
1014 F -= pref * (rdDa * rdQbr * rhat * (dv32-3.0*v32or));
1015 Ta += pref * ((-trQb*rxDa + 2.0 * DaxQbr)*v31 - rxDa*rdQbr*v32);
1016 Tb += pref * ((2.0*cross(DadQb, rhat) - 2.0*DaxQbr)*v31
1017 - 2.0*rdDa*rxQbr*v32);
1018 }
1019 }
1020
1021 if (a_is_Quadrupole) {
1022 if (b_is_Charge) {
1023 pref = pre14_ * *(idat.electroMult);
1024 U += C_b * pref * (v21 * trQa + v22 * rdQar);
1025 F += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or);
1026 F += C_b * pref * rdQar * rhat * (dv22 - 2.0*v22or);
1027 Ta += C_b * pref * 2.0 * rxQar * v22;
1028
1029 if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar);
1030 }
1031 if (b_is_Dipole) {
1032 pref = pre24_ * *(idat.electroMult);
1033 DbdQa = D_b * Q_a;
1034 DbdQar = dot(D_b, Qar);
1035 DbxQar = cross(D_b, Qar);
1036
1037 U += pref * ((trQa*rdDb + 2.0*DbdQar)*v31 + rdDb*rdQar*v32);
1038 F += pref * (trQa*D_b + 2.0*DbdQa) * v31or;
1039 F += pref * (trQa*rdDb + 2.0*DbdQar) * (dv31-v31or) * rhat;
1040 F += pref * (D_b*rdQar + 2.0*rdDb*rQa) * v32or;
1041 F += pref * (rdDb * rdQar * rhat * (dv32-3.0*v32or));
1042 Ta += pref * ((-2.0*cross(DbdQa, rhat) + 2.0*DbxQar)*v31
1043 + 2.0*rdDb*rxQar*v32);
1044 Tb += pref * ((trQa*rxDb - 2.0 * DbxQar)*v31 + rxDb*rdQar*v32);
1045 }
1046 if (b_is_Quadrupole) {
1047 pref = pre44_ * *(idat.electroMult); // yes
1048 QaQb = Q_a * Q_b;
1049 trQaQb = QaQb.trace();
1050 rQaQb = rhat * QaQb;
1051 QaQbr = QaQb * rhat;
1052 QaxQb = cross(Q_a, Q_b);
1053 rQaQbr = dot(rQa, Qbr);
1054 rQaxQbr = cross(rQa, Qbr);
1055
1056 U += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1057 U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42;
1058 U += pref * (rdQar * rdQbr) * v43;
1059
1060 F += pref * rhat * (trQa * trQb + 2.0 * trQaQb)*dv41;
1061 F += pref*rhat*(trQa*rdQbr + trQb*rdQar + 4.0*rQaQbr)*(dv42-2.0*v42or);
1062 F += pref * rhat * (rdQar * rdQbr)*(dv43 - 4.0*v43or);
1063
1064 F += pref * 2.0 * (trQb*rQa + trQa*rQb) * v42or;
1065 F += pref * 4.0 * (rQaQb + QaQbr) * v42or;
1066 F += pref * 2.0 * (rQa*rdQbr + rdQar*rQb) * v43or;
1067
1068 Ta += pref * (- 4.0 * QaxQb * v41);
1069 Ta += pref * (- 2.0 * trQb * cross(rQa, rhat)
1070 + 4.0 * cross(rhat, QaQbr)
1071 - 4.0 * rQaxQbr) * v42;
1072 Ta += pref * 2.0 * cross(rhat,Qar) * rdQbr * v43;
1073
1074
1075 Tb += pref * (+ 4.0 * QaxQb * v41);
1076 Tb += pref * (- 2.0 * trQa * cross(rQb, rhat)
1077 - 4.0 * cross(rQaQb, rhat)
1078 + 4.0 * rQaxQbr) * v42;
1079 // Possible replacement for line 2 above:
1080 // + 4.0 * cross(rhat, QbQar)
1081
1082 Tb += pref * 2.0 * cross(rhat,Qbr) * rdQar * v43;
1083
1084 // cerr << " tsum = " << Ta + Tb - cross( *(idat.d) , F ) << "\n";
1085 }
1086 }
1087
1088 if (idat.doElectricField) {
1089 *(idat.eField1) += Ea * *(idat.electroMult);
1090 *(idat.eField2) += Eb * *(idat.electroMult);
1091 }
1092
1093 if (a_is_Fluctuating) *(idat.dVdFQ1) += dUdCa * *(idat.sw);
1094 if (b_is_Fluctuating) *(idat.dVdFQ2) += dUdCb * *(idat.sw);
1095
1096 if (!idat.excluded) {
1097
1098 *(idat.vpair) += U;
1099 (*(idat.pot))[ELECTROSTATIC_FAMILY] += U * *(idat.sw);
1100 *(idat.f1) += F * *(idat.sw);
1101
1102 if (a_is_Dipole || a_is_Quadrupole)
1103 *(idat.t1) += Ta * *(idat.sw);
1104
1105 if (b_is_Dipole || b_is_Quadrupole)
1106 *(idat.t2) += Tb * *(idat.sw);
1107
1108 } else {
1109
1110 // only accumulate the forces and torques resulting from the
1111 // indirect reaction field terms.
1112
1113 *(idat.vpair) += indirect_Pot;
1114 (*(idat.excludedPot))[ELECTROSTATIC_FAMILY] += excluded_Pot;
1115 (*(idat.pot))[ELECTROSTATIC_FAMILY] += *(idat.sw) * indirect_Pot;
1116 *(idat.f1) += *(idat.sw) * indirect_F;
1117
1118 if (a_is_Dipole || a_is_Quadrupole)
1119 *(idat.t1) += *(idat.sw) * indirect_Ta;
1120
1121 if (b_is_Dipole || b_is_Quadrupole)
1122 *(idat.t2) += *(idat.sw) * indirect_Tb;
1123 }
1124 return;
1125 }
1126
1127 void Electrostatic::calcSelfCorrection(SelfData &sdat) {
1128
1129 if (!initialized_) initialize();
1130
1131 ElectrostaticAtomData data = ElectrostaticMap[sdat.atype];
1132
1133 // logicals
1134 bool i_is_Charge = data.is_Charge;
1135 bool i_is_Dipole = data.is_Dipole;
1136 bool i_is_Fluctuating = data.is_Fluctuating;
1137 RealType C_a = data.fixedCharge;
1138 RealType self, preVal, DadDa;
1139
1140 if (i_is_Fluctuating) {
1141 C_a += *(sdat.flucQ);
1142 // dVdFQ is really a force, so this is negative the derivative
1143 *(sdat.dVdFQ) -= *(sdat.flucQ) * data.hardness + data.electronegativity;
1144 (*(sdat.excludedPot))[ELECTROSTATIC_FAMILY] += (*sdat.flucQ) *
1145 (*(sdat.flucQ) * data.hardness * 0.5 + data.electronegativity);
1146 }
1147
1148 switch (summationMethod_) {
1149 case esm_REACTION_FIELD:
1150
1151 if (i_is_Charge) {
1152 // Self potential [see Wang and Hermans, "Reaction Field
1153 // Molecular Dynamics Simulation with Friedman’s Image Charge
1154 // Method," J. Phys. Chem. 99, 12001-12007 (1995).]
1155 preVal = pre11_ * preRF_ * C_a * C_a;
1156 (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= 0.5 * preVal / cutoffRadius_;
1157 }
1158
1159 if (i_is_Dipole) {
1160 DadDa = data.dipole.lengthSquare();
1161 (*(sdat.pot))[ELECTROSTATIC_FAMILY] -= pre22_ * preRF_ * DadDa;
1162 }
1163
1164 break;
1165
1166 case esm_SHIFTED_FORCE:
1167 case esm_SHIFTED_POTENTIAL:
1168 if (i_is_Charge) {
1169 self = - selfMult_ * C_a * (C_a + *(sdat.skippedCharge)) * pre11_;
1170 (*(sdat.pot))[ELECTROSTATIC_FAMILY] += self;
1171 }
1172 break;
1173 default:
1174 break;
1175 }
1176 }
1177
1178 RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1179 // This seems to work moderately well as a default. There's no
1180 // inherent scale for 1/r interactions that we can standardize.
1181 // 12 angstroms seems to be a reasonably good guess for most
1182 // cases.
1183 return 12.0;
1184 }
1185 }

Properties

Name Value
svn:eol-style native