OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Electrostatic.cpp
1/*
2 * Copyright (c) 2004-present, The University of Notre Dame. All rights
3 * reserved.
4 *
5 * Redistribution and use in source and binary forms, with or without
6 * modification, are permitted provided that the following conditions are met:
7 *
8 * 1. Redistributions of source code must retain the above copyright notice,
9 * this list of conditions and the following disclaimer.
10 *
11 * 2. Redistributions in binary form must reproduce the above copyright notice,
12 * this list of conditions and the following disclaimer in the documentation
13 * and/or other materials provided with the distribution.
14 *
15 * 3. Neither the name of the copyright holder nor the names of its
16 * contributors may be used to endorse or promote products derived from
17 * this software without specific prior written permission.
18 *
19 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
20 * AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
21 * IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
22 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
23 * LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
24 * CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
25 * SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
26 * INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
27 * CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
28 * ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
29 * POSSIBILITY OF SUCH DAMAGE.
30 *
31 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
32 * research, please cite the appropriate papers when you publish your
33 * work. Good starting points are:
34 *
35 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
36 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
37 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
38 * [4] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
39 * [5] Kuang & Gezelter, Mol. Phys., 110, 691-701 (2012).
40 * [6] Lamichhane, Gezelter & Newman, J. Chem. Phys. 141, 134109 (2014).
41 * [7] Lamichhane, Newman & Gezelter, J. Chem. Phys. 141, 134110 (2014).
42 * [8] Bhattarai, Newman & Gezelter, Phys. Rev. B 99, 094106 (2019).
43 */
44
45#include "nonbonded/Electrostatic.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <cstring>
50#include <memory>
51#include <numeric>
52
53#ifdef IS_MPI
54#include <mpi.h>
55#endif
56
57#include "flucq/FluctuatingChargeForces.hpp"
58#include "io/Globals.hpp"
59#include "math/SquareMatrix.hpp"
60#include "math/erfc.hpp"
61#include "nonbonded/SlaterIntegrals.hpp"
63#include "types/FixedChargeAdapter.hpp"
64#include "types/FluctuatingChargeAdapter.hpp"
65#include "types/MultipoleAdapter.hpp"
67#include "utils/Constants.hpp"
68#include "utils/simError.h"
69
70namespace OpenMD {
71
72 Electrostatic::Electrostatic() :
73 name_("Electrostatic"), initialized_(false), haveCutoffRadius_(false),
74 haveDampingAlpha_(false), haveDielectric_(false),
75 haveElectroSplines_(false), info_(NULL), forceField_(NULL)
76
77 {
78 flucQ_ = new FluctuatingChargeForces(info_);
79 }
80
81 Electrostatic::~Electrostatic() { delete flucQ_; }
82
83 void Electrostatic::setForceField(ForceField* ff) {
84 forceField_ = ff;
85 flucQ_->setForceField(forceField_);
86 }
87
88 void Electrostatic::setSimulatedAtomTypes(AtomTypeSet& simtypes) {
89 simTypes_ = simtypes;
90 flucQ_->setSimulatedAtomTypes(simTypes_);
91 }
92
93 void Electrostatic::initialize() {
94 Globals* simParams_ = info_->getSimParams();
95
96 summationMap_["HARD"] = esm_HARD;
97 summationMap_["NONE"] = esm_HARD;
98 summationMap_["SWITCHING_FUNCTION"] = esm_SWITCHING_FUNCTION;
99 summationMap_["SHIFTED_POTENTIAL"] = esm_SHIFTED_POTENTIAL;
100 summationMap_["SHIFTED_FORCE"] = esm_SHIFTED_FORCE;
101 summationMap_["TAYLOR_SHIFTED"] = esm_TAYLOR_SHIFTED;
102 summationMap_["REACTION_FIELD"] = esm_REACTION_FIELD;
103 summationMap_["EWALD_FULL"] = esm_EWALD_FULL;
104 summationMap_["EWALD_PME"] = esm_EWALD_PME;
105 summationMap_["EWALD_SPME"] = esm_EWALD_SPME;
106 screeningMap_["DAMPED"] = DAMPED;
107 screeningMap_["UNDAMPED"] = UNDAMPED;
108
109 // these prefactors convert the multipole interactions into kcal / mol
110 // all were computed assuming distances are measured in angstroms
111 // Charge-Charge, assuming charges are measured in electrons
112 pre11_ = 332.0637778;
113 // Charge-Dipole, assuming charges are measured in electrons, and
114 // dipoles are measured in debyes
115 pre12_ = 69.13373;
116 // Dipole-Dipole, assuming dipoles are measured in Debye
117 pre22_ = 14.39325;
118 // Charge-Quadrupole, assuming charges are measured in electrons, and
119 // quadrupoles are measured in 10^-26 esu cm^2
120 // This unit is also known affectionately as an esu centibarn.
121 pre14_ = 69.13373;
122 // Dipole-Quadrupole, assuming dipoles are measured in debyes and
123 // quadrupoles in esu centibarns:
124 pre24_ = 14.39325;
125 // Quadrupole-Quadrupole, assuming esu centibarns:
126 pre44_ = 14.39325;
127
128 // conversions for the simulation box dipole moment
129 chargeToC_ = 1.60217733e-19;
130 angstromToM_ = 1.0e-10;
131 debyeToCm_ = 3.33564095198e-30;
132
133 // Default number of points for electrostatic splines
134 np_ = 100;
135
136 // variables to handle different summation methods for long-range
137 // electrostatics:
138 summationMethod_ = esm_HARD;
139 screeningMethod_ = UNDAMPED;
140 dielectric_ = 1.0;
141
142 // check the summation method:
143 if (simParams_->haveElectrostaticSummationMethod()) {
144 string myMethod = simParams_->getElectrostaticSummationMethod();
145 toUpper(myMethod);
146 map<string, ElectrostaticSummationMethod>::iterator i;
147 i = summationMap_.find(myMethod);
148 if (i != summationMap_.end()) {
149 summationMethod_ = (*i).second;
150 } else {
151 // throw error
152 snprintf(
153 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
154 "Electrostatic::initialize: Unknown electrostaticSummationMethod.\n"
155 "\t(Input file specified %s .)\n"
156 "\telectrostaticSummationMethod must be one of: \"hard\",\n"
157 "\t\"shifted_potential\", \"shifted_force\",\n"
158 "\t\"taylor_shifted\", or \"reaction_field\".\n",
159 myMethod.c_str());
160 painCave.isFatal = 1;
161 simError();
162 }
163 } else {
164 // set ElectrostaticSummationMethod to the cutoffMethod:
165 if (simParams_->haveCutoffMethod()) {
166 string myMethod = simParams_->getCutoffMethod();
167 toUpper(myMethod);
168 map<string, ElectrostaticSummationMethod>::iterator i;
169 i = summationMap_.find(myMethod);
170 if (i != summationMap_.end()) { summationMethod_ = (*i).second; }
171 }
172 }
173
174 if (summationMethod_ == esm_REACTION_FIELD) {
175 if (!simParams_->haveDielectric()) {
176 // throw warning
177 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
178 "SimInfo warning: dielectric was not specified in the input "
179 "file\n\tfor "
180 "the reaction field correction method.\n"
181 "\tA default value of %f will be used for the dielectric.\n",
182 dielectric_);
183 painCave.isFatal = 0;
184 painCave.severity = OPENMD_INFO;
185 simError();
186 } else {
187 dielectric_ = simParams_->getDielectric();
188 }
189 haveDielectric_ = true;
190 }
191
192 if (simParams_->haveElectrostaticScreeningMethod()) {
193 string myScreen = simParams_->getElectrostaticScreeningMethod();
194 toUpper(myScreen);
195 map<string, ElectrostaticScreeningMethod>::iterator i;
196 i = screeningMap_.find(myScreen);
197 if (i != screeningMap_.end()) {
198 screeningMethod_ = (*i).second;
199 } else {
200 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
201 "SimInfo error: Unknown electrostaticScreeningMethod.\n"
202 "\t(Input file specified %s .)\n"
203 "\telectrostaticScreeningMethod must be one of: \"undamped\"\n"
204 "or \"damped\".\n",
205 myScreen.c_str());
206 painCave.isFatal = 1;
207 simError();
208 }
209 }
210
211 // check to make sure a cutoff value has been set:
212 if (!haveCutoffRadius_) {
213 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
214 "Electrostatic::initialize has no Default "
215 "Cutoff value!\n");
216 painCave.severity = OPENMD_ERROR;
217 painCave.isFatal = 1;
218 simError();
219 }
220
221 if (screeningMethod_ == DAMPED || summationMethod_ == esm_EWALD_FULL) {
222 if (!simParams_->haveDampingAlpha()) {
223 haveDampingAlpha_ = false;
224 // first compute a cutoff dependent alpha value
225 // we assume alpha depends linearly with rcut from 0 to 20.5 ang
226 dampingAlpha_ = 0.425 - cutoffRadius_ * 0.02;
227 if (dampingAlpha_ < 0.0) {
228 screeningMethod_ = UNDAMPED;
229 dampingAlpha_ = 0.0;
230 // throw warning
231 snprintf(
232 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
233 "Electrostatic::initialize: dampingAlpha was not specified in "
234 "the\n"
235 "\tinput file, but the computed value would be 0.0 with a\n"
236 "\tcutoff of %f (ang). Switching to UNDAMPED electrostatics.\n",
237 cutoffRadius_);
238 painCave.severity = OPENMD_INFO;
239 painCave.isFatal = 0;
240 simError();
241 } else {
242 // throw warning
243 snprintf(
244 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
245 "Electrostatic::initialize: dampingAlpha was not specified in "
246 "the\n"
247 "\tinput file. A default value of %f (1/ang) will be used for "
248 "the\n"
249 "\tcutoff of %f (ang).\n",
250 dampingAlpha_, cutoffRadius_);
251 painCave.severity = OPENMD_INFO;
252 painCave.isFatal = 0;
253 simError();
254 haveDampingAlpha_ = true;
255 }
256 } else {
257 dampingAlpha_ = simParams_->getDampingAlpha();
258 haveDampingAlpha_ = true;
259 }
260 }
261
262 Etypes.clear();
263 Etids.clear();
264 FQtypes.clear();
265 FQtids.clear();
266 ElectrostaticMap.clear();
267 Jij.clear();
268 nElectro_ = 0;
269 nFlucq_ = 0;
270
271 Etids.resize(forceField_->getNAtomType(), -1);
272 FQtids.resize(forceField_->getNAtomType(), -1);
273
274 AtomTypeSet::iterator at;
275 for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
276 if ((*at)->isElectrostatic()) nElectro_++;
277 if ((*at)->isFluctuatingCharge()) nFlucq_++;
278 }
279
280 Jij.resize(nFlucq_);
281
282 for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
283 if ((*at)->isElectrostatic()) addType(*at);
284 }
285
286 if (summationMethod_ == esm_REACTION_FIELD) {
287 preRF_ = (dielectric_ - 1.0) /
288 ((2.0 * dielectric_ + 1.0) * pow(cutoffRadius_, 3));
289 }
290
291 RealType b0c, b1c, b2c, b3c, b4c, b5c;
292 RealType db0c_1, db0c_2, db0c_3, db0c_4, db0c_5;
293 RealType a2, expTerm, invArootPi(0.0);
294
295 RealType r = cutoffRadius_;
296 RealType r2 = r * r;
297 RealType ric = 1.0 / r;
298 RealType ric2 = ric * ric;
299
300 if (screeningMethod_ == DAMPED) {
301 a2 = dampingAlpha_ * dampingAlpha_;
302 invArootPi = 1.0 / (dampingAlpha_ * sqrt(Constants::PI));
303 expTerm = exp(-a2 * r2);
304 // values of Smith's B_l functions at the cutoff radius:
305 b0c = erfc(dampingAlpha_ * r) / r;
306 b1c = (b0c + 2.0 * a2 * expTerm * invArootPi) / r2;
307 b2c = (3.0 * b1c + pow(2.0 * a2, 2) * expTerm * invArootPi) / r2;
308 b3c = (5.0 * b2c + pow(2.0 * a2, 3) * expTerm * invArootPi) / r2;
309 b4c = (7.0 * b3c + pow(2.0 * a2, 4) * expTerm * invArootPi) / r2;
310 b5c = (9.0 * b4c + pow(2.0 * a2, 5) * expTerm * invArootPi) / r2;
311 // Half the Smith self piece:
312 selfMult1_ = -a2 * invArootPi;
313 selfMult2_ = -2.0 * a2 * a2 * invArootPi / 3.0;
314 selfMult4_ = -4.0 * a2 * a2 * a2 * invArootPi / 5.0;
315 } else {
316 a2 = 0.0;
317 b0c = 1.0 / r;
318 b1c = (b0c) / r2;
319 b2c = (3.0 * b1c) / r2;
320 b3c = (5.0 * b2c) / r2;
321 b4c = (7.0 * b3c) / r2;
322 b5c = (9.0 * b4c) / r2;
323 selfMult1_ = 0.0;
324 selfMult2_ = 0.0;
325 selfMult4_ = 0.0;
326 }
327
328 // higher derivatives of B_0 at the cutoff radius:
329 db0c_1 = -r * b1c;
330 db0c_2 = -b1c + r2 * b2c;
331 db0c_3 = 3.0 * r * b2c - r2 * r * b3c;
332 db0c_4 = 3.0 * b2c - 6.0 * r2 * b3c + r2 * r2 * b4c;
333 db0c_5 = -15.0 * r * b3c + 10.0 * r2 * r * b4c - r2 * r2 * r * b5c;
334
335 if (summationMethod_ != esm_EWALD_FULL) {
336 selfMult1_ -= b0c;
337 selfMult2_ += (db0c_2 + 2.0 * db0c_1 * ric) / 3.0;
338 selfMult4_ -= (db0c_4 + 4.0 * db0c_3 * ric) / 15.0;
339 }
340
341 // working variables for the splines:
342 RealType ri, ri2;
343 RealType b0, b1, b2, b3, b4, b5;
344 RealType db0_1, db0_2, db0_3, db0_4, db0_5;
345 RealType f, fc, f0;
346 RealType g, gc, g0, g1, g2, g3, g4;
347 RealType h, hc, h1, h2, h3, h4;
348 RealType s, sc, s2, s3, s4;
349 RealType t, tc, t3, t4;
350 RealType u, uc, u4;
351
352 // working variables for Taylor expansion:
353 RealType rmRc, rmRc2, rmRc3, rmRc4;
354
355 // Approximate using splines using a maximum of 0.1 Angstroms
356 // between points.
357 int nptest = int((cutoffRadius_ + 2.0) / 0.1);
358 np_ = (np_ > nptest) ? np_ : nptest;
359
360 // Add a 2 angstrom safety window to deal with cutoffGroups that
361 // have charged atoms longer than the cutoffRadius away from each
362 // other. Splining is almost certainly the best choice here.
363 // Direct calls to erfc would be preferrable if it is a very fast
364 // implementation.
365
366 RealType dx = (cutoffRadius_ + 2.0) / RealType(np_);
367
368 // Storage vectors for the computed functions
369 vector<RealType> rv;
370 vector<RealType> v01v;
371 vector<RealType> v11v;
372 vector<RealType> v21v, v22v;
373 vector<RealType> v31v, v32v;
374 vector<RealType> v41v, v42v, v43v;
375
376 for (int i = 1; i < np_ + 1; i++) {
377 r = RealType(i) * dx;
378 rv.push_back(r);
379
380 ri = 1.0 / r;
381 ri2 = ri * ri;
382
383 r2 = r * r;
384 expTerm = exp(-a2 * r2);
385
386 // Taylor expansion factors (no need for factorials this way):
387 rmRc = r - cutoffRadius_;
388 rmRc2 = rmRc * rmRc / 2.0;
389 rmRc3 = rmRc2 * rmRc / 3.0;
390 rmRc4 = rmRc3 * rmRc / 4.0;
391
392 // values of Smith's B_l functions at r:
393 if (screeningMethod_ == DAMPED) {
394 b0 = erfc(dampingAlpha_ * r) * ri;
395 b1 = (b0 + 2.0 * a2 * expTerm * invArootPi) * ri2;
396 b2 = (3.0 * b1 + pow(2.0 * a2, 2) * expTerm * invArootPi) * ri2;
397 b3 = (5.0 * b2 + pow(2.0 * a2, 3) * expTerm * invArootPi) * ri2;
398 b4 = (7.0 * b3 + pow(2.0 * a2, 4) * expTerm * invArootPi) * ri2;
399 b5 = (9.0 * b4 + pow(2.0 * a2, 5) * expTerm * invArootPi) * ri2;
400 } else {
401 b0 = ri;
402 b1 = (b0)*ri2;
403 b2 = (3.0 * b1) * ri2;
404 b3 = (5.0 * b2) * ri2;
405 b4 = (7.0 * b3) * ri2;
406 b5 = (9.0 * b4) * ri2;
407 }
408
409 // higher derivatives of B_0 at r:
410 db0_1 = -r * b1;
411 db0_2 = -b1 + r2 * b2;
412 db0_3 = 3.0 * r * b2 - r2 * r * b3;
413 db0_4 = 3.0 * b2 - 6.0 * r2 * b3 + r2 * r2 * b4;
414 db0_5 = -15.0 * r * b3 + 10.0 * r2 * r * b4 - r2 * r2 * r * b5;
415
416 f = b0;
417 fc = b0c;
418 f0 = f - fc - rmRc * db0c_1;
419
420 g = db0_1;
421 gc = db0c_1;
422 g0 = g - gc;
423 g1 = g0 - rmRc * db0c_2;
424 g2 = g1 - rmRc2 * db0c_3;
425 g3 = g2 - rmRc3 * db0c_4;
426 g4 = g3 - rmRc4 * db0c_5;
427
428 h = db0_2;
429 hc = db0c_2;
430 h1 = h - hc;
431 h2 = h1 - rmRc * db0c_3;
432 h3 = h2 - rmRc2 * db0c_4;
433 h4 = h3 - rmRc3 * db0c_5;
434
435 s = db0_3;
436 sc = db0c_3;
437 s2 = s - sc;
438 s3 = s2 - rmRc * db0c_4;
439 s4 = s3 - rmRc2 * db0c_5;
440
441 t = db0_4;
442 tc = db0c_4;
443 t3 = t - tc;
444 t4 = t3 - rmRc * db0c_5;
445
446 u = db0_5;
447 uc = db0c_5;
448 u4 = u - uc;
449
450 // in what follows below, the various v functions are used for
451 // potentials and torques, while the w functions show up in the
452 // forces.
453
454 switch (summationMethod_) {
455 case esm_SHIFTED_FORCE:
456
457 v01 = f - fc - rmRc * gc;
458 v11 = g - gc - rmRc * hc;
459 v21 = g * ri - gc * ric - rmRc * (hc - gc * ric) * ric;
460 v22 =
461 h - g * ri - (hc - gc * ric) - rmRc * (sc - (hc - gc * ric) * ric);
462 v31 = (h - g * ri) * ri - (hc - gc * ric) * ric -
463 rmRc * (sc - 2.0 * (hc - gc * ric) * ric) * ric;
464 v32 = (s - 3.0 * (h - g * ri) * ri) -
465 (sc - 3.0 * (hc - gc * ric) * ric) -
466 rmRc * (tc - 3.0 * (sc - 2.0 * (hc - gc * ric) * ric) * ric);
467 v41 = (h - g * ri) * ri2 - (hc - gc * ric) * ric2 -
468 rmRc * (sc - 3.0 * (hc - gc * ric) * ric) * ric2;
469 v42 =
470 (s - 3.0 * (h - g * ri) * ri) * ri -
471 (sc - 3.0 * (hc - gc * ric) * ric) * ric -
472 rmRc * (tc - (4.0 * sc - 9.0 * (hc - gc * ric) * ric) * ric) * ric;
473
474 v43 =
475 (t - 3.0 * (2.0 * s - 5.0 * (h - g * ri) * ri) * ri) -
476 (tc - 3.0 * (2.0 * sc - 5.0 * (hc - gc * ric) * ric) * ric) -
477 rmRc * (uc - 3.0 *
478 (2.0 * tc -
479 (7.0 * sc - 15.0 * (hc - gc * ric) * ric) * ric) *
480 ric);
481
482 dv01 = g - gc;
483 dv11 = h - hc;
484 dv21 = (h - g * ri) * ri - (hc - gc * ric) * ric;
485 dv22 = (s - (h - g * ri) * ri) - (sc - (hc - gc * ric) * ric);
486 dv31 = (s - 2.0 * (h - g * ri) * ri) * ri -
487 (sc - 2.0 * (hc - gc * ric) * ric) * ric;
488 dv32 = (t - 3.0 * (s - 2.0 * (h - g * ri) * ri) * ri) -
489 (tc - 3.0 * (sc - 2.0 * (hc - gc * ric) * ric) * ric);
490 dv41 = (s - 3.0 * (h - g * ri) * ri) * ri2 -
491 (sc - 3.0 * (hc - gc * ric) * ric) * ric2;
492 dv42 = (t - (4.0 * s - 9.0 * (h - g * ri) * ri) * ri) * ri -
493 (tc - (4.0 * sc - 9.0 * (hc - gc * ric) * ric) * ric) * ric;
494 dv43 =
495 (u -
496 3.0 * (2.0 * t - (7.0 * s - 15.0 * (h - g * ri) * ri) * ri) * ri) -
497 (uc -
498 3.0 *
499 (2.0 * tc - (7.0 * sc - 15.0 * (hc - gc * ric) * ric) * ric) *
500 ric);
501
502 break;
503
504 case esm_TAYLOR_SHIFTED:
505
506 v01 = f0;
507 v11 = g1;
508 v21 = g2 * ri;
509 v22 = h2 - v21;
510 v31 = (h3 - g3 * ri) * ri;
511 v32 = s3 - 3.0 * v31;
512 v41 = (h4 - g4 * ri) * ri2;
513 v42 = s4 * ri - 3.0 * v41;
514 v43 = t4 - 6.0 * v42 - 3.0 * v41;
515
516 dv01 = g0;
517 dv11 = h1;
518 dv21 = (h2 - g2 * ri) * ri;
519 dv22 = (s2 - (h2 - g2 * ri) * ri);
520 dv31 = (s3 - 2.0 * (h3 - g3 * ri) * ri) * ri;
521 dv32 = (t3 - 3.0 * (s3 - 2.0 * (h3 - g3 * ri) * ri) * ri);
522 dv41 = (s4 - 3.0 * (h4 - g4 * ri) * ri) * ri2;
523 dv42 = (t4 - (4.0 * s4 - 9.0 * (h4 - g4 * ri) * ri) * ri) * ri;
524 dv43 =
525 (u4 -
526 3.0 * (2.0 * t4 - (7.0 * s4 - 15.0 * (h4 - g4 * ri) * ri) * ri) *
527 ri);
528
529 break;
530
531 case esm_SHIFTED_POTENTIAL:
532
533 v01 = f - fc;
534 v11 = g - gc;
535 v21 = g * ri - gc * ric;
536 v22 = h - g * ri - (hc - gc * ric);
537 v31 = (h - g * ri) * ri - (hc - gc * ric) * ric;
538 v32 =
539 (s - 3.0 * (h - g * ri) * ri) - (sc - 3.0 * (hc - gc * ric) * ric);
540 v41 = (h - g * ri) * ri2 - (hc - gc * ric) * ric2;
541 v42 = (s - 3.0 * (h - g * ri) * ri) * ri -
542 (sc - 3.0 * (hc - gc * ric) * ric) * ric;
543 v43 = (t - 3.0 * (2.0 * s - 5.0 * (h - g * ri) * ri) * ri) -
544 (tc - 3.0 * (2.0 * sc - 5.0 * (hc - gc * ric) * ric) * ric);
545
546 dv01 = g;
547 dv11 = h;
548 dv21 = (h - g * ri) * ri;
549 dv22 = (s - (h - g * ri) * ri);
550 dv31 = (s - 2.0 * (h - g * ri) * ri) * ri;
551 dv32 = (t - 3.0 * (s - 2.0 * (h - g * ri) * ri) * ri);
552 dv41 = (s - 3.0 * (h - g * ri) * ri) * ri2;
553 dv42 = (t - (4.0 * s - 9.0 * (h - g * ri) * ri) * ri) * ri;
554 dv43 =
555 (u -
556 3.0 * (2.0 * t - (7.0 * s - 15.0 * (h - g * ri) * ri) * ri) * ri);
557
558 break;
559
560 case esm_SWITCHING_FUNCTION:
561 case esm_HARD:
562 case esm_EWALD_FULL:
563
564 v01 = f;
565 v11 = g;
566 v21 = g * ri;
567 v22 = h - g * ri;
568 v31 = (h - g * ri) * ri;
569 v32 = (s - 3.0 * (h - g * ri) * ri);
570 v41 = (h - g * ri) * ri2;
571 v42 = (s - 3.0 * (h - g * ri) * ri) * ri;
572 v43 = (t - 3.0 * (2.0 * s - 5.0 * (h - g * ri) * ri) * ri);
573
574 dv01 = g;
575 dv11 = h;
576 dv21 = (h - g * ri) * ri;
577 dv22 = (s - (h - g * ri) * ri);
578 dv31 = (s - 2.0 * (h - g * ri) * ri) * ri;
579 dv32 = (t - 3.0 * (s - 2.0 * (h - g * ri) * ri) * ri);
580 dv41 = (s - 3.0 * (h - g * ri) * ri) * ri2;
581 dv42 = (t - (4.0 * s - 9.0 * (h - g * ri) * ri) * ri) * ri;
582 dv43 =
583 (u -
584 3.0 * (2.0 * t - (7.0 * s - 15.0 * (h - g * ri) * ri) * ri) * ri);
585
586 break;
587
588 case esm_REACTION_FIELD:
589
590 // following DL_POLY's lead for shifting the image charge potential:
591 f = b0 + preRF_ * r2;
592 fc = b0c + preRF_ * cutoffRadius_ * cutoffRadius_;
593
594 g = db0_1 + preRF_ * 2.0 * r;
595 gc = db0c_1 + preRF_ * 2.0 * cutoffRadius_;
596
597 h = db0_2 + preRF_ * 2.0;
598 hc = db0c_2 + preRF_ * 2.0;
599
600 v01 = f - fc;
601 v11 = g - gc;
602 v21 = g * ri - gc * ric;
603 v22 = h - g * ri - (hc - gc * ric);
604 v31 = (h - g * ri) * ri - (hc - gc * ric) * ric;
605 v32 =
606 (s - 3.0 * (h - g * ri) * ri) - (sc - 3.0 * (hc - gc * ric) * ric);
607 v41 = (h - g * ri) * ri2 - (hc - gc * ric) * ric2;
608 v42 = (s - 3.0 * (h - g * ri) * ri) * ri -
609 (sc - 3.0 * (hc - gc * ric) * ric) * ric;
610 v43 = (t - 3.0 * (2.0 * s - 5.0 * (h - g * ri) * ri) * ri) -
611 (tc - 3.0 * (2.0 * sc - 5.0 * (hc - gc * ric) * ric) * ric);
612
613 dv01 = g;
614 dv11 = h;
615 dv21 = (h - g * ri) * ri;
616 dv22 = (s - (h - g * ri) * ri);
617 dv31 = (s - 2.0 * (h - g * ri) * ri) * ri;
618 dv32 = (t - 3.0 * (s - 2.0 * (h - g * ri) * ri) * ri);
619 dv41 = (s - 3.0 * (h - g * ri) * ri) * ri2;
620 dv42 = (t - (4.0 * s - 9.0 * (h - g * ri) * ri) * ri) * ri;
621 dv43 =
622 (u -
623 3.0 * (2.0 * t - (7.0 * s - 15.0 * (h - g * ri) * ri) * ri) * ri);
624
625 break;
626
627 case esm_EWALD_PME:
628 case esm_EWALD_SPME:
629 default:
630 map<string, ElectrostaticSummationMethod>::iterator i;
631 std::string meth;
632 for (i = summationMap_.begin(); i != summationMap_.end(); ++i) {
633 if ((*i).second == summationMethod_) meth = (*i).first;
634 }
635 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
636 "Electrostatic::initialize: electrostaticSummationMethod %s \n"
637 "\thas not been implemented yet. Please select one of:\n"
638 "\t\"hard\", \"shifted_potential\", or \"shifted_force\"\n",
639 meth.c_str());
640 painCave.isFatal = 1;
641 simError();
642 break;
643 }
644
645 // Add these computed values to the storage vectors for spline creation:
646 v01v.push_back(v01);
647 v11v.push_back(v11);
648 v21v.push_back(v21);
649 v22v.push_back(v22);
650 v31v.push_back(v31);
651 v32v.push_back(v32);
652 v41v.push_back(v41);
653 v42v.push_back(v42);
654 v43v.push_back(v43);
655 }
656
657 // construct the spline structures and fill them with the values we've
658 // computed:
659
660 v01s = std::make_shared<CubicSpline>();
661 v01s->addPoints(rv, v01v);
662 v11s = std::make_shared<CubicSpline>();
663 v11s->addPoints(rv, v11v);
664 v21s = std::make_shared<CubicSpline>();
665 v21s->addPoints(rv, v21v);
666 v22s = std::make_shared<CubicSpline>();
667 v22s->addPoints(rv, v22v);
668 v31s = std::make_shared<CubicSpline>();
669 v31s->addPoints(rv, v31v);
670 v32s = std::make_shared<CubicSpline>();
671 v32s->addPoints(rv, v32v);
672 v41s = std::make_shared<CubicSpline>();
673 v41s->addPoints(rv, v41v);
674 v42s = std::make_shared<CubicSpline>();
675 v42s->addPoints(rv, v42v);
676 v43s = std::make_shared<CubicSpline>();
677 v43s->addPoints(rv, v43v);
678
679 haveElectroSplines_ = true;
680
681 initialized_ = true;
682 }
683
684 void Electrostatic::addType(AtomType* atomType) {
685 ElectrostaticAtomData electrostaticAtomData;
686 electrostaticAtomData.is_Charge = false;
687 electrostaticAtomData.is_Dipole = false;
688 electrostaticAtomData.is_Quadrupole = false;
689 electrostaticAtomData.is_Fluctuating = false;
690 electrostaticAtomData.uses_SlaterIntramolecular = false;
691
693
694 if (fca.isFixedCharge()) {
695 electrostaticAtomData.is_Charge = true;
696 electrostaticAtomData.fixedCharge = fca.getCharge();
697 }
698
699 MultipoleAdapter ma = MultipoleAdapter(atomType);
700 if (ma.isMultipole()) {
701 if (ma.isDipole()) {
702 electrostaticAtomData.is_Dipole = true;
703 electrostaticAtomData.dipole = ma.getDipole();
704 }
705 if (ma.isQuadrupole()) {
706 electrostaticAtomData.is_Quadrupole = true;
707 electrostaticAtomData.quadrupole = ma.getQuadrupole();
708 }
709 }
710
712
713 if (fqa.isFluctuatingCharge()) {
714 electrostaticAtomData.is_Fluctuating = true;
715 electrostaticAtomData.uses_SlaterIntramolecular =
716 fqa.usesSlaterIntramolecular();
717 electrostaticAtomData.electronegativity = fqa.getElectronegativity();
718 electrostaticAtomData.hardness = fqa.getHardness();
719 electrostaticAtomData.slaterN = fqa.getSlaterN();
720 electrostaticAtomData.slaterZeta = fqa.getSlaterZeta();
721 }
722
723 int atid = atomType->getIdent();
724 int etid = Etypes.size();
725 int fqtid = FQtypes.size();
726
727 pair<set<int>::iterator, bool> ret;
728 ret = Etypes.insert(atid);
729 if (ret.second == false) {
730 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
731 "Electrostatic already had a previous entry with ident %d\n",
732 atid);
733 painCave.severity = OPENMD_INFO;
734 painCave.isFatal = 0;
735 simError();
736 }
737
738 Etids[atid] = etid;
739 ElectrostaticMap.push_back(electrostaticAtomData);
740
741 if (electrostaticAtomData.is_Fluctuating) {
742 ret = FQtypes.insert(atid);
743 if (ret.second == false) {
744 snprintf(
745 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
746 "Electrostatic already had a previous fluctuating charge entry "
747 "with ident %d\n",
748 atid);
749 painCave.severity = OPENMD_INFO;
750 painCave.isFatal = 0;
751 simError();
752 }
753 FQtids[atid] = fqtid;
754 Jij[fqtid].resize(nFlucq_);
755
756 // Now, iterate over all known fluctuating and add to the
757 // coulomb integral map:
758
759 std::set<int>::iterator it;
760 for (it = FQtypes.begin(); it != FQtypes.end(); ++it) {
761 int etid2 = Etids[(*it)];
762 int fqtid2 = FQtids[(*it)];
763 ElectrostaticAtomData eaData2 = ElectrostaticMap[etid2];
764 RealType a = electrostaticAtomData.slaterZeta;
765 RealType b = eaData2.slaterZeta;
766 int m = electrostaticAtomData.slaterN;
767 int n = eaData2.slaterN;
768 CubicSplinePtr J {std::make_shared<CubicSpline>()};
769
770 // do both types actually use Slater orbitals?
771
772 if ((electrostaticAtomData.uses_SlaterIntramolecular &&
773 eaData2.uses_SlaterIntramolecular)) {
774 // Create the spline of the coulombic integral for s-type
775 // Slater orbitals. Add a 2 angstrom safety window to deal
776 // with cutoffGroups that have charged atoms longer than the
777 // cutoffRadius away from each other.
778
779 RealType rval;
780 RealType dr = (cutoffRadius_ + 2.0) / RealType(np_ - 1);
781 vector<RealType> rvals;
782 vector<RealType> Jvals;
783 // RealType j0, j0c, j1c;
784 // don't start at i = 0, as rval = 0 is undefined for the
785 // slater overlap integrals.
786 for (int i = 1; i < np_ + 1; i++) {
787 rval = RealType(i) * dr;
788 rvals.push_back(rval);
789
790 Jvals.push_back(
791 sSTOCoulInt(a, b, m, n, rval * Constants::angstromToBohr) *
792 Constants::hartreeToKcal);
793 }
794
795 J->addPoints(rvals, Jvals);
796 }
797 Jij[fqtid][fqtid2] = J;
798 Jij[fqtid2].resize(nFlucq_);
799 Jij[fqtid2][fqtid] = J;
800 }
801 }
802 return;
803 }
804
805 void Electrostatic::setCutoffRadius(RealType rCut) {
806 cutoffRadius_ = rCut;
807 haveCutoffRadius_ = true;
808 }
809
810 void Electrostatic::setElectrostaticSummationMethod(
812 summationMethod_ = esm;
813 }
814 void Electrostatic::setElectrostaticScreeningMethod(
815 ElectrostaticScreeningMethod sm) {
816 screeningMethod_ = sm;
817 }
818 void Electrostatic::setDampingAlpha(RealType alpha) {
819 dampingAlpha_ = alpha;
820 haveDampingAlpha_ = true;
821 }
822 void Electrostatic::setReactionFieldDielectric(RealType dielectric) {
823 dielectric_ = dielectric;
824 haveDielectric_ = true;
825 }
826
827 void Electrostatic::calcForce(InteractionData& idat) {
828 if (!initialized_) initialize();
829
830 if (Etids[idat.atid1] != -1) {
831 data1 = ElectrostaticMap[Etids[idat.atid1]];
832 a_is_Charge = data1.is_Charge;
833 a_is_Dipole = data1.is_Dipole;
834 a_is_Quadrupole = data1.is_Quadrupole;
835 a_is_Fluctuating = data1.is_Fluctuating;
836 a_uses_SlaterIntra = data1.uses_SlaterIntramolecular;
837
838 } else {
839 a_is_Charge = false;
840 a_is_Dipole = false;
841 a_is_Quadrupole = false;
842 a_is_Fluctuating = false;
843 a_uses_SlaterIntra = false;
844 }
845 if (Etids[idat.atid2] != -1) {
846 data2 = ElectrostaticMap[Etids[idat.atid2]];
847 b_is_Charge = data2.is_Charge;
848 b_is_Dipole = data2.is_Dipole;
849 b_is_Quadrupole = data2.is_Quadrupole;
850 b_is_Fluctuating = data2.is_Fluctuating;
851 b_uses_SlaterIntra = data2.uses_SlaterIntramolecular;
852
853 } else {
854 b_is_Charge = false;
855 b_is_Dipole = false;
856 b_is_Quadrupole = false;
857 b_is_Fluctuating = false;
858 b_uses_SlaterIntra = false;
859 }
860
861 U = 0.0; // Potential
862 F.zero(); // Force
863 Ta.zero(); // Torque on site a
864 Tb.zero(); // Torque on site b
865 Ea.zero(); // Electric field at site a
866 Eb.zero(); // Electric field at site b
867 Pa = 0.0; // Site potential at site a
868 Pb = 0.0; // Site potential at site b
869 dUdCa = 0.0; // fluctuating charge force at site a
870 dUdCb = 0.0; // fluctuating charge force at site a
871
872 // Indirect interactions mediated by the reaction field.
873 indirect_Pot = 0.0; // Potential
874 indirect_F.zero(); // Force
875 indirect_Ta.zero(); // Torque on site a
876 indirect_Tb.zero(); // Torque on site b
877
878 // Excluded potential that is still computed for fluctuating charges
879 excluded_Pot = 0.0;
880
881 // some variables we'll need independent of electrostatic type:
882
883 ri = 1.0 / idat.rij;
884 rhat = idat.d * ri;
885
886 // Obtain all of the required radial function values from the
887 // spline structures:
888
889 if (((a_uses_SlaterIntra || b_uses_SlaterIntra) && idat.excluded)) {
890 J = Jij[FQtids[idat.atid1]][FQtids[idat.atid2]];
891 }
892
893 // needed for fields (and forces):
894 if (a_is_Charge || b_is_Charge) {
895 v01s->getValueAndDerivativeAt(idat.rij, v01, dv01);
896 }
897 if (a_is_Dipole || b_is_Dipole) {
898 v11s->getValueAndDerivativeAt(idat.rij, v11, dv11);
899 v11or = ri * v11;
900 }
901 if (a_is_Quadrupole || b_is_Quadrupole || (a_is_Dipole && b_is_Dipole)) {
902 v21s->getValueAndDerivativeAt(idat.rij, v21, dv21);
903 v22s->getValueAndDerivativeAt(idat.rij, v22, dv22);
904 v22or = ri * v22;
905 }
906
907 // needed for potentials (and forces and torques):
908 if ((a_is_Dipole && b_is_Quadrupole) || (b_is_Dipole && a_is_Quadrupole)) {
909 v31s->getValueAndDerivativeAt(idat.rij, v31, dv31);
910 v32s->getValueAndDerivativeAt(idat.rij, v32, dv32);
911 v31or = v31 * ri;
912 v32or = v32 * ri;
913 }
914 if (a_is_Quadrupole && b_is_Quadrupole) {
915 v41s->getValueAndDerivativeAt(idat.rij, v41, dv41);
916 v42s->getValueAndDerivativeAt(idat.rij, v42, dv42);
917 v43s->getValueAndDerivativeAt(idat.rij, v43, dv43);
918 v42or = v42 * ri;
919 v43or = v43 * ri;
920 }
921
922 // calculate the single-site contributions (fields, etc).
923
924 if (a_is_Charge) {
925 C_a = data1.fixedCharge;
926
927 if (a_is_Fluctuating) { C_a += idat.flucQ1; }
928
929 if (idat.excluded) {
930 idat.skippedCharge2 += C_a;
931 } else {
932 // only do the field and site potentials if we're not excluded:
933 Eb -= C_a * pre11_ * dv01 * rhat;
934 Pb += C_a * pre11_ * v01;
935 }
936 }
937
938 if (a_is_Dipole) {
939 rdDa = dot(rhat, idat.D_1);
940 rxDa = cross(rhat, idat.D_1);
941 if (!idat.excluded) {
942 Eb -= pre12_ * ((dv11 - v11or) * rdDa * rhat + v11or * idat.D_1);
943 Pb += pre12_ * v11 * rdDa;
944 }
945 }
946
947 if (a_is_Quadrupole) {
948 trQa = idat.Q_1.trace();
949 Qar = idat.Q_1 * rhat;
950 rQa = rhat * idat.Q_1;
951 rdQar = dot(rhat, Qar);
952 rxQar = cross(rhat, Qar);
953 if (!idat.excluded) {
954 Eb -= pre14_ * (trQa * rhat * dv21 + 2.0 * Qar * v22or +
955 rdQar * rhat * (dv22 - 2.0 * v22or));
956 Pb += pre14_ * (v21 * trQa + v22 * rdQar);
957 }
958 }
959
960 if (b_is_Charge) {
961 C_b = data2.fixedCharge;
962
963 if (b_is_Fluctuating) { C_b += idat.flucQ2; }
964
965 if (idat.excluded) {
966 idat.skippedCharge1 += C_b;
967 } else {
968 // only do the field if we're not excluded:
969 Ea += C_b * pre11_ * dv01 * rhat;
970 Pa += C_b * pre11_ * v01;
971 }
972 }
973
974 if (b_is_Dipole) {
975 rdDb = dot(rhat, idat.D_2);
976 rxDb = cross(rhat, idat.D_2);
977 if (!idat.excluded) {
978 Ea += pre12_ * ((dv11 - v11or) * rdDb * rhat + v11or * idat.D_2);
979 Pa += pre12_ * v11 * rdDb;
980 }
981 }
982
983 if (b_is_Quadrupole) {
984 trQb = idat.Q_2.trace();
985 Qbr = idat.Q_2 * rhat;
986 rQb = rhat * idat.Q_2;
987 rdQbr = dot(rhat, Qbr);
988 rxQbr = cross(rhat, Qbr);
989 if (!idat.excluded) {
990 Ea += pre14_ * (trQb * rhat * dv21 + 2.0 * Qbr * v22or +
991 rdQbr * rhat * (dv22 - 2.0 * v22or));
992 Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
993 }
994 }
995
996 if (a_is_Charge) {
997 if (b_is_Charge) {
998 pref = pre11_ * idat.electroMult;
999 U += C_a * C_b * pref * v01;
1000 F += C_a * C_b * pref * dv01 * rhat;
1001
1002 // If this is an excluded pair, there are still indirect
1003 // interactions via the reaction field we must worry about:
1004
1005 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1006 rfContrib = preRF_ * pref * C_a * C_b * idat.r2;
1007 indirect_Pot += rfContrib;
1008 indirect_F += rfContrib * 2.0 * ri * rhat;
1009 }
1010
1011 // Fluctuating charge forces are handled via Coulomb integrals
1012 // for excluded pairs (i.e. those connected via bonds) and
1013 // with the standard charge-charge interaction otherwise.
1014
1015 if (idat.excluded) {
1016 if (a_uses_SlaterIntra || b_uses_SlaterIntra) {
1017 coulInt = J->getValueAt(idat.rij);
1018 excluded_Pot += C_a * C_b * coulInt;
1019 if (a_is_Fluctuating) dUdCa += C_b * coulInt;
1020 if (b_is_Fluctuating) dUdCb += C_a * coulInt;
1021 }
1022 } else {
1023 if (a_is_Fluctuating) dUdCa += C_b * pref * v01;
1024 if (b_is_Fluctuating) dUdCb += C_a * pref * v01;
1025 }
1026 }
1027
1028 if (b_is_Dipole) {
1029 pref = pre12_ * idat.electroMult;
1030 U += C_a * pref * v11 * rdDb;
1031 F += C_a * pref * ((dv11 - v11or) * rdDb * rhat + v11or * idat.D_2);
1032 Tb += C_a * pref * v11 * rxDb;
1033
1034 if (a_is_Fluctuating) dUdCa += pref * v11 * rdDb;
1035
1036 // Even if we excluded this pair from direct interactions, we
1037 // still have the reaction-field-mediated charge-dipole
1038 // interaction:
1039
1040 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1041 rfContrib = C_a * pref * preRF_ * 2.0 * idat.rij;
1042 indirect_Pot += rfContrib * rdDb;
1043 indirect_F += rfContrib * idat.D_2 / idat.rij;
1044 indirect_Tb += C_a * pref * preRF_ * rxDb;
1045 }
1046 }
1047
1048 if (b_is_Quadrupole) {
1049 pref = pre14_ * idat.electroMult;
1050 U += C_a * pref * (v21 * trQb + v22 * rdQbr);
1051 F += C_a * pref * (trQb * dv21 * rhat + 2.0 * Qbr * v22or);
1052 F += C_a * pref * rdQbr * rhat * (dv22 - 2.0 * v22or);
1053 Tb += C_a * pref * 2.0 * rxQbr * v22;
1054
1055 if (a_is_Fluctuating) dUdCa += pref * (v21 * trQb + v22 * rdQbr);
1056 }
1057 }
1058
1059 if (a_is_Dipole) {
1060 if (b_is_Charge) {
1061 pref = pre12_ * idat.electroMult;
1062
1063 U -= C_b * pref * v11 * rdDa;
1064 F -= C_b * pref * ((dv11 - v11or) * rdDa * rhat + v11or * idat.D_1);
1065 Ta -= C_b * pref * v11 * rxDa;
1066
1067 if (b_is_Fluctuating) dUdCb -= pref * v11 * rdDa;
1068
1069 // Even if we excluded this pair from direct interactions,
1070 // we still have the reaction-field-mediated charge-dipole
1071 // interaction:
1072 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1073 rfContrib = C_b * pref * preRF_ * 2.0 * idat.rij;
1074 indirect_Pot -= rfContrib * rdDa;
1075 indirect_F -= rfContrib * idat.D_1 / idat.rij;
1076 indirect_Ta -= C_b * pref * preRF_ * rxDa;
1077 }
1078 }
1079
1080 if (b_is_Dipole) {
1081 pref = pre22_ * idat.electroMult;
1082 DadDb = dot(idat.D_1, idat.D_2);
1083 DaxDb = cross(idat.D_1, idat.D_2);
1084
1085 U -= pref * (DadDb * v21 + rdDa * rdDb * v22);
1086 F -= pref * (dv21 * DadDb * rhat +
1087 v22or * (rdDb * idat.D_1 + rdDa * idat.D_2));
1088 F -= pref * (rdDa * rdDb) * (dv22 - 2.0 * v22or) * rhat;
1089 Ta += pref * (v21 * DaxDb - v22 * rdDb * rxDa);
1090 Tb += pref * (-v21 * DaxDb - v22 * rdDa * rxDb);
1091 // Even if we excluded this pair from direct interactions, we
1092 // still have the reaction-field-mediated dipole-dipole
1093 // interaction:
1094 if (summationMethod_ == esm_REACTION_FIELD && idat.excluded) {
1095 rfContrib = -pref * preRF_ * 2.0;
1096 indirect_Pot += rfContrib * DadDb;
1097 indirect_Ta += rfContrib * DaxDb;
1098 indirect_Tb -= rfContrib * DaxDb;
1099 }
1100 }
1101
1102 if (b_is_Quadrupole) {
1103 pref = pre24_ * idat.electroMult;
1104 DadQb = idat.D_1 * idat.Q_2;
1105 DadQbr = dot(idat.D_1, Qbr);
1106 DaxQbr = cross(idat.D_1, Qbr);
1107
1108 U -= pref * ((trQb * rdDa + 2.0 * DadQbr) * v31 + rdDa * rdQbr * v32);
1109 F -= pref * (trQb * idat.D_1 + 2.0 * DadQb) * v31or;
1110 F -= pref * (trQb * rdDa + 2.0 * DadQbr) * (dv31 - v31or) * rhat;
1111 F -= pref * (idat.D_1 * rdQbr + 2.0 * rdDa * rQb) * v32or;
1112 F -= pref * (rdDa * rdQbr * rhat * (dv32 - 3.0 * v32or));
1113 Ta += pref * ((-trQb * rxDa + 2.0 * DaxQbr) * v31 - rxDa * rdQbr * v32);
1114 Tb += pref * ((2.0 * cross(DadQb, rhat) - 2.0 * DaxQbr) * v31 -
1115 2.0 * rdDa * rxQbr * v32);
1116 }
1117 }
1118
1119 if (a_is_Quadrupole) {
1120 if (b_is_Charge) {
1121 pref = pre14_ * idat.electroMult;
1122 U += C_b * pref * (v21 * trQa + v22 * rdQar);
1123 F += C_b * pref * (trQa * rhat * dv21 + 2.0 * Qar * v22or);
1124 F += C_b * pref * rdQar * rhat * (dv22 - 2.0 * v22or);
1125 Ta += C_b * pref * 2.0 * rxQar * v22;
1126
1127 if (b_is_Fluctuating) dUdCb += pref * (v21 * trQa + v22 * rdQar);
1128 }
1129 if (b_is_Dipole) {
1130 pref = pre24_ * idat.electroMult;
1131 DbdQa = idat.D_2 * idat.Q_1;
1132 DbdQar = dot(idat.D_2, Qar);
1133 DbxQar = cross(idat.D_2, Qar);
1134
1135 U += pref * ((trQa * rdDb + 2.0 * DbdQar) * v31 + rdDb * rdQar * v32);
1136 F += pref * (trQa * idat.D_2 + 2.0 * DbdQa) * v31or;
1137 F += pref * (trQa * rdDb + 2.0 * DbdQar) * (dv31 - v31or) * rhat;
1138 F += pref * (idat.D_2 * rdQar + 2.0 * rdDb * rQa) * v32or;
1139 F += pref * (rdDb * rdQar * rhat * (dv32 - 3.0 * v32or));
1140 Ta += pref * ((-2.0 * cross(DbdQa, rhat) + 2.0 * DbxQar) * v31 +
1141 2.0 * rdDb * rxQar * v32);
1142 Tb += pref * ((trQa * rxDb - 2.0 * DbxQar) * v31 + rxDb * rdQar * v32);
1143 }
1144 if (b_is_Quadrupole) {
1145 pref = pre44_ * idat.electroMult; // yes
1146 QaQb = idat.Q_1 * idat.Q_2;
1147 trQaQb = QaQb.trace();
1148 rQaQb = rhat * QaQb;
1149 QaQbr = QaQb * rhat;
1150 QaxQb = mCross(idat.Q_1, idat.Q_2);
1151 rQaQbr = dot(rQa, Qbr);
1152 rQaxQbr = cross(rQa, Qbr);
1153
1154 U += pref * (trQa * trQb + 2.0 * trQaQb) * v41;
1155 U += pref * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) * v42;
1156 U += pref * (rdQar * rdQbr) * v43;
1157
1158 F += pref * rhat * (trQa * trQb + 2.0 * trQaQb) * dv41;
1159 F += pref * rhat * (trQa * rdQbr + trQb * rdQar + 4.0 * rQaQbr) *
1160 (dv42 - 2.0 * v42or);
1161 F += pref * rhat * (rdQar * rdQbr) * (dv43 - 4.0 * v43or);
1162
1163 F += pref * 2.0 * (trQb * rQa + trQa * rQb) * v42or;
1164 F += pref * 4.0 * (rQaQb + QaQbr) * v42or;
1165 F += pref * 2.0 * (rQa * rdQbr + rdQar * rQb) * v43or;
1166
1167 Ta += pref * (-4.0 * QaxQb * v41);
1168 Ta += pref *
1169 (-2.0 * trQb * cross(rQa, rhat) + 4.0 * cross(rhat, QaQbr) -
1170 4.0 * rQaxQbr) *
1171 v42;
1172 Ta += pref * 2.0 * cross(rhat, Qar) * rdQbr * v43;
1173
1174 Tb += pref * (+4.0 * QaxQb * v41);
1175 Tb += pref *
1176 (-2.0 * trQa * cross(rQb, rhat) - 4.0 * cross(rQaQb, rhat) +
1177 4.0 * rQaxQbr) *
1178 v42;
1179 // Possible replacement for line 2 above:
1180 // + 4.0 * cross(rhat, QbQar)
1181
1182 Tb += pref * 2.0 * cross(rhat, Qbr) * rdQar * v43;
1183 }
1184 }
1185
1186 if (idat.doElectricField) {
1187 idat.eField1 += Ea * idat.electroMult;
1188 idat.eField2 += Eb * idat.electroMult;
1189 }
1190
1191 if (idat.doSitePotential) {
1192 idat.sPot1 += Pa * idat.electroMult;
1193 idat.sPot2 += Pb * idat.electroMult;
1194 }
1195
1196 if (a_is_Fluctuating) idat.dVdFQ1 += dUdCa * idat.sw;
1197 if (b_is_Fluctuating) idat.dVdFQ2 += dUdCb * idat.sw;
1198
1199 if (!idat.excluded) {
1200 idat.vpair += U;
1201 idat.pot[ELECTROSTATIC_FAMILY] += U * idat.sw;
1202 if (idat.isSelected) idat.selePot[ELECTROSTATIC_FAMILY] += U * idat.sw;
1203
1204 idat.f1 += F * idat.sw;
1205
1206 if (a_is_Dipole || a_is_Quadrupole) idat.t1 += Ta * idat.sw;
1207
1208 if (b_is_Dipole || b_is_Quadrupole) idat.t2 += Tb * idat.sw;
1209
1210 } else {
1211 // only accumulate the forces and torques resulting from the
1212 // indirect reaction field terms.
1213
1214 idat.vpair += indirect_Pot;
1215 idat.excludedPot[ELECTROSTATIC_FAMILY] += excluded_Pot;
1216 idat.pot[ELECTROSTATIC_FAMILY] += idat.sw * indirect_Pot;
1217 if (idat.isSelected)
1218 idat.selePot[ELECTROSTATIC_FAMILY] += idat.sw * indirect_Pot;
1219
1220 idat.f1 += idat.sw * indirect_F;
1221
1222 if (a_is_Dipole || a_is_Quadrupole) idat.t1 += idat.sw * indirect_Ta;
1223
1224 if (b_is_Dipole || b_is_Quadrupole) idat.t2 += idat.sw * indirect_Tb;
1225 }
1226 return;
1227 }
1228
1229 void Electrostatic::calcSelfCorrection(SelfData& sdat) {
1230 if (!initialized_) initialize();
1231
1232 ElectrostaticAtomData data = ElectrostaticMap[Etids[sdat.atid]];
1233
1234 // logicals
1235 bool i_is_Charge = data.is_Charge;
1236 bool i_is_Dipole = data.is_Dipole;
1237 bool i_is_Quadrupole = data.is_Quadrupole;
1238 bool i_is_Fluctuating = data.is_Fluctuating;
1239 RealType C_a = data.fixedCharge;
1240 RealType selfPot(0.0), fqf(0.0), preVal, DdD(0.0), trQ, trQQ;
1241
1242 if (i_is_Dipole) { DdD = data.dipole.lengthSquare(); }
1243
1244 if (i_is_Fluctuating) {
1245 // We're now doing all of the self pieces for fluctuating charges in
1246 // explicit self interactions.
1247 C_a += sdat.flucQ;
1248 flucQ_->getSelfInteraction(sdat.atid, sdat.flucQ, selfPot, fqf);
1249 }
1250
1251 switch (summationMethod_) {
1252 case esm_REACTION_FIELD:
1253
1254 if (i_is_Charge) {
1255 // Self potential [see Wang and Hermans, "Reaction Field
1256 // Molecular Dynamics Simulation with Friedman’s Image Charge
1257 // Method," J. Phys. Chem. 99, 12001-12007 (1995).]
1258 preVal = pre11_ * preRF_ * C_a * C_a;
1259 selfPot -= 0.5 * preVal / cutoffRadius_;
1260 // if (i_is_Fluctuating) {
1261 // fqf += pre11_ * preRF_ * C_a / cutoffRadius_;
1262 //}
1263 }
1264
1265 if (i_is_Dipole) { selfPot -= pre22_ * preRF_ * DdD; }
1266
1267 break;
1268
1269 case esm_SHIFTED_FORCE:
1270 case esm_SHIFTED_POTENTIAL:
1271 case esm_TAYLOR_SHIFTED:
1272 case esm_EWALD_FULL:
1273 if (i_is_Charge) {
1274 selfPot += selfMult1_ * pre11_ * pow(C_a + sdat.skippedCharge, 2);
1275 if (i_is_Fluctuating) {
1276 fqf -= selfMult1_ * pre11_ * 2.0 * (C_a + sdat.skippedCharge);
1277 }
1278 }
1279 if (i_is_Dipole) selfPot += selfMult2_ * pre22_ * DdD;
1280 if (i_is_Quadrupole) {
1281 trQ = data.quadrupole.trace();
1282 trQQ = (data.quadrupole * data.quadrupole).trace();
1283 selfPot += selfMult4_ * pre44_ * (2.0 * trQQ + trQ * trQ);
1284 if (i_is_Charge) {
1285 selfPot -= selfMult2_ * pre14_ * 2.0 * C_a * trQ;
1286 if (i_is_Fluctuating) { fqf += selfMult2_ * pre14_ * 2.0 * trQ; }
1287 }
1288 }
1289 break;
1290 default:
1291 break;
1292 }
1293
1294 sdat.selfPot[ELECTROSTATIC_FAMILY] += selfPot;
1295
1296 if (sdat.isSelected) sdat.selePot[ELECTROSTATIC_FAMILY] += selfPot;
1297
1298 if (sdat.doParticlePot) { sdat.particlePot += selfPot; }
1299
1300 if (i_is_Fluctuating) sdat.flucQfrc += fqf;
1301 }
1302
1303 void Electrostatic::calcSurfaceTerm(bool slabGeometry, int axis,
1304 RealType& pot) {
1305 SimInfo::MoleculeIterator mi;
1306 Molecule::AtomIterator ai;
1307 RealType C;
1308 Vector3d r;
1309 Vector3d D;
1310 Vector3d t;
1311 Vector3d netDipole(0.0);
1312 int atid;
1314
1315 const RealType mPoleConverter = 0.20819434; // converts from the
1316 // internal units of
1317 // Debye (for dipoles)
1318 // or Debye-angstroms
1319 // (for quadrupoles) to
1320 // electron angstroms or
1321 // electron-angstroms^2
1322
1323 const RealType eConverter = 332.0637778; // convert the
1324 // Charge-Charge
1325 // electrostatic
1326 // interactions into kcal /
1327 // mol assuming distances
1328 // are measured in
1329 // angstroms.
1330
1331 if (!initialized_) initialize();
1332
1333 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1334 mol = info_->nextMolecule(mi)) {
1335 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1336 atom = mol->nextAtom(ai)) {
1337 atid = atom->getAtomType()->getIdent();
1338 data = ElectrostaticMap[Etids[atid]];
1339
1340 if (data.is_Charge) {
1341 C = data.fixedCharge;
1342 if (data.is_Fluctuating) C += atom->getFlucQPos();
1343
1344 r = atom->getPos();
1346
1347 netDipole += C * r;
1348 }
1349
1350 if (data.is_Dipole) {
1351 D = atom->getDipole() * mPoleConverter;
1352 netDipole += D;
1353 }
1354 }
1355 }
1356
1357#ifdef IS_MPI
1358 MPI_Allreduce(MPI_IN_PLACE, netDipole.getArrayPointer(), 3, MPI_REALTYPE,
1359 MPI_SUM, MPI_COMM_WORLD);
1360#endif
1361
1362 RealType V = info_->getSnapshotManager()->getCurrentSnapshot()->getVolume();
1363 RealType prefactor;
1364
1365 if (slabGeometry) {
1366 prefactor = 2.0 * Constants::PI / V;
1367 // Compute complementary axes to the privileged axis
1368 int axis1 = (axis + 1) % 3;
1369 int axis2 = (axis + 2) % 3;
1370 netDipole[axis1] = 0.0;
1371 netDipole[axis2] = 0.0;
1372 } else {
1373 prefactor = 2.0 * Constants::PI / (3.0 * V);
1374 }
1375
1376 pot += eConverter * prefactor * netDipole.lengthSquare();
1377
1378 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1379 mol = info_->nextMolecule(mi)) {
1380 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1381 atom = mol->nextAtom(ai)) {
1382 atom->addElectricField(-eConverter * prefactor * 2.0 * netDipole);
1383
1384 atid = atom->getAtomType()->getIdent();
1385 data = ElectrostaticMap[Etids[atid]];
1386
1387 if (data.is_Charge) {
1388 C = data.fixedCharge;
1389 if (data.is_Fluctuating) {
1390 r = atom->getPos();
1392 atom->addFlucQFrc(-eConverter * prefactor * 2.0 *
1393 dot(r, netDipole));
1394 }
1395 atom->addFrc(-eConverter * prefactor * 2.0 * C * netDipole);
1396 }
1397
1398 if (data.is_Dipole) {
1399 D = atom->getDipole() * mPoleConverter;
1400 t = -eConverter * prefactor * 2.0 * cross(D, netDipole);
1401 atom->addTrq(t);
1402 }
1403 }
1404 }
1405 }
1406
1407 RealType Electrostatic::getSuggestedCutoffRadius(pair<AtomType*, AtomType*>) {
1408 // This seems to work moderately well as a default. There's no
1409 // inherent scale for 1/r interactions that we can standardize.
1410 // 12 angstroms seems to be a reasonably good guess for most
1411 // cases.
1412 return 12.0;
1413 }
1414
1415 void Electrostatic::ReciprocalSpaceSum(RealType& pot) {
1416 RealType kPot = 0.0;
1417 RealType kVir = 0.0;
1418 // Mat3x3d kVirTens(0.0);
1419
1420 const RealType mPoleConverter = 0.20819434; // converts from the
1421 // internal units of
1422 // Debye (for dipoles)
1423 // or Debye-angstroms
1424 // (for quadrupoles) to
1425 // electron angstroms or
1426 // electron-angstroms^2
1427
1428 const RealType eConverter = 332.0637778; // convert the
1429 // Charge-Charge
1430 // electrostatic
1431 // interactions into kcal /
1432 // mol assuming distances
1433 // are measured in
1434 // angstroms.
1435
1437 Vector3d box = hmat.diagonals();
1438 RealType boxMax = box.max();
1439
1440 // int kMax = int(2.0 * Constants::PI / (pow(dampingAlpha_,2)*cutoffRadius_
1441 // * boxMax) );
1442 int kMax = 7;
1443 int kSqMax = kMax * kMax + 2;
1444
1445 int kLimit = kMax + 1;
1446 int kLim2 = 2 * kMax + 1;
1447 int kSqLim = kSqMax;
1448
1449 vector<RealType> AK(kSqLim + 1, 0.0);
1450 RealType xcl = 2.0 * Constants::PI / box.x();
1451 RealType ycl = 2.0 * Constants::PI / box.y();
1452 RealType zcl = 2.0 * Constants::PI / box.z();
1453 RealType rcl = 2.0 * Constants::PI / boxMax;
1454 RealType rvol = 2.0 * Constants::PI / (box.x() * box.y() * box.z());
1455
1456 if (dampingAlpha_ < 1.0e-12) return;
1457
1458 RealType ralph = -0.25 / pow(dampingAlpha_, 2);
1459
1460 // Calculate and store exponential factors
1461
1462 vector<vector<RealType>> elc;
1463 vector<vector<RealType>> emc;
1464 vector<vector<RealType>> enc;
1465 vector<vector<RealType>> els;
1466 vector<vector<RealType>> ems;
1467 vector<vector<RealType>> ens;
1468
1469 int nMax = info_->getNAtoms();
1470
1471 elc.resize(kLimit + 1);
1472 emc.resize(kLimit + 1);
1473 enc.resize(kLimit + 1);
1474 els.resize(kLimit + 1);
1475 ems.resize(kLimit + 1);
1476 ens.resize(kLimit + 1);
1477
1478 for (int j = 0; j < kLimit + 1; j++) {
1479 elc[j].resize(nMax);
1480 emc[j].resize(nMax);
1481 enc[j].resize(nMax);
1482 els[j].resize(nMax);
1483 ems[j].resize(nMax);
1484 ens[j].resize(nMax);
1485 }
1486
1487 Vector3d t(2.0 * Constants::PI);
1488 t.Vdiv(t, box);
1489
1490 SimInfo::MoleculeIterator mi;
1491 Molecule::AtomIterator ai;
1492 int i;
1493 Vector3d r;
1494 Vector3d tt;
1495
1496 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1497 mol = info_->nextMolecule(mi)) {
1498 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1499 atom = mol->nextAtom(ai)) {
1500 i = atom->getLocalIndex();
1501 r = atom->getPos();
1503
1504 tt.Vmul(t, r);
1505
1506 elc[1][i] = 1.0;
1507 emc[1][i] = 1.0;
1508 enc[1][i] = 1.0;
1509 els[1][i] = 0.0;
1510 ems[1][i] = 0.0;
1511 ens[1][i] = 0.0;
1512
1513 elc[2][i] = cos(tt.x());
1514 emc[2][i] = cos(tt.y());
1515 enc[2][i] = cos(tt.z());
1516 els[2][i] = sin(tt.x());
1517 ems[2][i] = sin(tt.y());
1518 ens[2][i] = sin(tt.z());
1519
1520 for (int l = 3; l <= kLimit; l++) {
1521 elc[l][i] = elc[l - 1][i] * elc[2][i] - els[l - 1][i] * els[2][i];
1522 emc[l][i] = emc[l - 1][i] * emc[2][i] - ems[l - 1][i] * ems[2][i];
1523 enc[l][i] = enc[l - 1][i] * enc[2][i] - ens[l - 1][i] * ens[2][i];
1524 els[l][i] = els[l - 1][i] * elc[2][i] + elc[l - 1][i] * els[2][i];
1525 ems[l][i] = ems[l - 1][i] * emc[2][i] + emc[l - 1][i] * ems[2][i];
1526 ens[l][i] = ens[l - 1][i] * enc[2][i] + enc[l - 1][i] * ens[2][i];
1527 }
1528 }
1529 }
1530
1531 // Calculate and store AK coefficients:
1532
1533 RealType eksq = 1.0;
1534 RealType expf = 0.0;
1535 if (ralph < 0.0) expf = exp(ralph * rcl * rcl);
1536 for (i = 1; i <= kSqLim; i++) {
1537 RealType rksq = float(i) * rcl * rcl;
1538 eksq = expf * eksq;
1539 AK[i] = eConverter * eksq / rksq;
1540 }
1541
1542 /*
1543 * Loop over all k vectors k = 2 pi (ll/Lx, mm/Ly, nn/Lz)
1544 * the values of ll, mm and nn are selected so that the symmetry of
1545 * reciprocal lattice is taken into account i.e. the following
1546 * rules apply.
1547 *
1548 * ll ranges over the values 0 to kMax only.
1549 *
1550 * mm ranges over 0 to kMax when ll=0 and over
1551 * -kMax to kMax otherwise.
1552 * nn ranges over 1 to kMax when ll=mm=0 and over
1553 * -kMax to kMax otherwise.
1554 *
1555 * Hence the result of the summation must be doubled at the end.
1556 */
1557
1558 std::vector<RealType> clm(nMax, 0.0);
1559 std::vector<RealType> slm(nMax, 0.0);
1560 std::vector<RealType> ckr(nMax, 0.0);
1561 std::vector<RealType> skr(nMax, 0.0);
1562 std::vector<RealType> ckc(nMax, 0.0);
1563 std::vector<RealType> cks(nMax, 0.0);
1564 std::vector<RealType> dkc(nMax, 0.0);
1565 std::vector<RealType> dks(nMax, 0.0);
1566 std::vector<RealType> qkc(nMax, 0.0);
1567 std::vector<RealType> qks(nMax, 0.0);
1568 std::vector<Vector3d> dxk(nMax, V3Zero);
1569 std::vector<Vector3d> qxk(nMax, V3Zero);
1570 RealType rl, rm, rn;
1571 Vector3d kVec;
1572 Vector3d Qk;
1573 // RealType k2;
1574 Mat3x3d Kmat;
1575 RealType ckcs, ckss, dkcs, dkss, qkcs, qkss;
1576 int atid;
1578 RealType C, dk, qk;
1579 Vector3d D;
1580 Mat3x3d Q;
1581
1582 int mMin = kLimit;
1583 int nMin = kLimit + 1;
1584 for (int l = 1; l <= kLimit; l++) {
1585 int ll = l - 1;
1586 rl = xcl * float(ll);
1587 for (int mmm = mMin; mmm <= kLim2; mmm++) {
1588 int mm = mmm - kLimit;
1589 int m = abs(mm) + 1;
1590 rm = ycl * float(mm);
1591 // Set temporary products of exponential terms
1592 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1593 mol = info_->nextMolecule(mi)) {
1594 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1595 atom = mol->nextAtom(ai)) {
1596 i = atom->getLocalIndex();
1597 if (mm < 0) {
1598 clm[i] = elc[l][i] * emc[m][i] + els[l][i] * ems[m][i];
1599 slm[i] = els[l][i] * emc[m][i] - ems[m][i] * elc[l][i];
1600 } else {
1601 clm[i] = elc[l][i] * emc[m][i] - els[l][i] * ems[m][i];
1602 slm[i] = els[l][i] * emc[m][i] + ems[m][i] * elc[l][i];
1603 }
1604 }
1605 }
1606 for (int nnn = nMin; nnn <= kLim2; nnn++) {
1607 int nn = nnn - kLimit;
1608 int n = abs(nn) + 1;
1609 rn = zcl * float(nn);
1610 // Test on magnitude of k vector:
1611 int kk = ll * ll + mm * mm + nn * nn;
1612 if (kk <= kSqLim) {
1613 kVec = Vector3d(rl, rm, rn);
1614 // k2 = dot(kVec, kVec); // length^2 of kVec
1615 Kmat = outProduct(kVec, kVec); // kMatrix
1616 // Calculate exp(ikr) terms
1617 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1618 mol = info_->nextMolecule(mi)) {
1619 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1620 atom = mol->nextAtom(ai)) {
1621 i = atom->getLocalIndex();
1622
1623 if (nn < 0) {
1624 ckr[i] = clm[i] * enc[n][i] + slm[i] * ens[n][i];
1625 skr[i] = slm[i] * enc[n][i] - clm[i] * ens[n][i];
1626
1627 } else {
1628 ckr[i] = clm[i] * enc[n][i] - slm[i] * ens[n][i];
1629 skr[i] = slm[i] * enc[n][i] + clm[i] * ens[n][i];
1630 }
1631 }
1632 }
1633
1634 // Calculate scalar and vector products for each site:
1635
1636 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1637 mol = info_->nextMolecule(mi)) {
1638 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1639 atom = mol->nextAtom(ai)) {
1640 i = atom->getLocalIndex();
1641 int atid = atom->getAtomType()->getIdent();
1642 data = ElectrostaticMap[Etids[atid]];
1643
1644 if (data.is_Charge) {
1645 C = data.fixedCharge;
1646 if (data.is_Fluctuating) C += atom->getFlucQPos();
1647 ckc[i] = C * ckr[i];
1648 cks[i] = C * skr[i];
1649 }
1650
1651 if (data.is_Dipole) {
1652 D = atom->getDipole() * mPoleConverter;
1653 dk = dot(D, kVec);
1654 dxk[i] = cross(D, kVec);
1655 dkc[i] = dk * ckr[i];
1656 dks[i] = dk * skr[i];
1657 }
1658 if (data.is_Quadrupole) {
1659 Q = atom->getQuadrupole() * mPoleConverter;
1660 Qk = Q * kVec;
1661 qk = dot(kVec, Qk);
1662 qxk[i] = -cross(kVec, Qk);
1663 qkc[i] = qk * ckr[i];
1664 qks[i] = qk * skr[i];
1665 }
1666 }
1667 }
1668
1669 // calculate vector sums
1670
1671 ckcs = std::accumulate(ckc.begin(), ckc.end(), 0.0);
1672 ckss = std::accumulate(cks.begin(), cks.end(), 0.0);
1673 dkcs = std::accumulate(dkc.begin(), dkc.end(), 0.0);
1674 dkss = std::accumulate(dks.begin(), dks.end(), 0.0);
1675 qkcs = std::accumulate(qkc.begin(), qkc.end(), 0.0);
1676 qkss = std::accumulate(qks.begin(), qks.end(), 0.0);
1677
1678#ifdef IS_MPI
1679 MPI_Allreduce(MPI_IN_PLACE, &ckcs, 1, MPI_REALTYPE, MPI_SUM,
1680 MPI_COMM_WORLD);
1681 MPI_Allreduce(MPI_IN_PLACE, &ckss, 1, MPI_REALTYPE, MPI_SUM,
1682 MPI_COMM_WORLD);
1683 MPI_Allreduce(MPI_IN_PLACE, &dkcs, 1, MPI_REALTYPE, MPI_SUM,
1684 MPI_COMM_WORLD);
1685 MPI_Allreduce(MPI_IN_PLACE, &dkss, 1, MPI_REALTYPE, MPI_SUM,
1686 MPI_COMM_WORLD);
1687 MPI_Allreduce(MPI_IN_PLACE, &qkcs, 1, MPI_REALTYPE, MPI_SUM,
1688 MPI_COMM_WORLD);
1689 MPI_Allreduce(MPI_IN_PLACE, &qkss, 1, MPI_REALTYPE, MPI_SUM,
1690 MPI_COMM_WORLD);
1691#endif
1692
1693 // Accumulate potential energy and virial contribution:
1694
1695 kPot += 2.0 * rvol * AK[kk] *
1696 ((ckss + dkcs - qkss) * (ckss + dkcs - qkss) +
1697 (ckcs - dkss - qkcs) * (ckcs - dkss - qkcs));
1698
1699 kVir +=
1700 2.0 * rvol * AK[kk] *
1701 (ckcs * ckcs + ckss * ckss + 4.0 * (ckss * dkcs - ckcs * dkss) +
1702 3.0 * (dkcs * dkcs + dkss * dkss) -
1703 6.0 * (ckss * qkss + ckcs * qkcs) +
1704 8.0 * (dkss * qkcs - dkcs * qkss) +
1705 5.0 * (qkss * qkss + qkcs * qkcs));
1706
1707 // kVirTens += 2 * rvol * AK[kk] *
1708 // ( Mat3x3d::identity() - 2.0*( 1.0 / k2 - ralph) * Kmat ) *
1709 // (ckcs * ckcs + ckss * ckss + 4.0 * (ckss * dkcs - ckcs * dkss)
1710 // +
1711 // 3.0 * (dkcs * dkcs + dkss * dkss) -
1712 // 6.0 * (ckss * qkss + ckcs * qkcs) +
1713 // 8.0 * (dkss * qkcs - dkcs * qkss) +
1714 // 5.0 * (qkss * qkss + qkcs * qkcs));
1715
1716 // Calculate force and torque for each site:
1717
1718 for (Molecule* mol = info_->beginMolecule(mi); mol != NULL;
1719 mol = info_->nextMolecule(mi)) {
1720 for (Atom* atom = mol->beginAtom(ai); atom != NULL;
1721 atom = mol->nextAtom(ai)) {
1722 i = atom->getLocalIndex();
1723 atid = atom->getAtomType()->getIdent();
1724 data = ElectrostaticMap[Etids[atid]];
1725
1726 RealType qfrc =
1727 AK[kk] *
1728 ((cks[i] + dkc[i] - qks[i]) * (ckcs - dkss - qkcs) -
1729 (ckc[i] - dks[i] - qkc[i]) * (ckss + dkcs - qkss));
1730 RealType qtrq1 = AK[kk] * (skr[i] * (ckcs - dkss - qkcs) -
1731 ckr[i] * (ckss + dkcs - qkss));
1732 RealType qtrq2 = 2.0 * AK[kk] *
1733 (ckr[i] * (ckcs - dkss - qkcs) +
1734 skr[i] * (ckss + dkcs - qkss));
1735
1736 atom->addFrc(4.0 * rvol * qfrc * kVec);
1737
1738 if (data.is_Fluctuating) {
1739 atom->addFlucQFrc(-2.0 * rvol * qtrq2);
1740 }
1741
1742 if (data.is_Dipole) {
1743 atom->addTrq(4.0 * rvol * qtrq1 * dxk[i]);
1744 }
1745 if (data.is_Quadrupole) {
1746 atom->addTrq(4.0 * rvol * qtrq2 * qxk[i]);
1747 }
1748 }
1749 }
1750 }
1751 }
1752 nMin = 1;
1753 }
1754 mMin = 1;
1755 }
1756 pot += kPot;
1757 // virialTensor += kVirTens;
1758 }
1759
1760 void Electrostatic::getSitePotentials(Atom* a1, Atom* a2, bool excluded,
1761 RealType& spot1, RealType& spot2) {
1762 if (!initialized_) { initialize(); }
1763
1764 const RealType mPoleConverter = 0.20819434;
1765
1766 AtomType* atype1 = a1->getAtomType();
1767 AtomType* atype2 = a2->getAtomType();
1768 int atid1 = atype1->getIdent();
1769 int atid2 = atype2->getIdent();
1770 data1 = ElectrostaticMap[Etids[atid1]];
1771 data2 = ElectrostaticMap[Etids[atid2]];
1772
1773 Pa = 0.0; // Site potential at site a
1774 Pb = 0.0; // Site potential at site b
1775
1776 Vector3d d = a2->getPos() - a1->getPos();
1778 RealType rij = d.length();
1779 // some variables we'll need independent of electrostatic type:
1780
1781 RealType ri = 1.0 / rij;
1782 rhat = d * ri;
1783
1784 if ((rij >= cutoffRadius_) || excluded) {
1785 spot1 = 0.0;
1786 spot2 = 0.0;
1787 return;
1788 }
1789
1790 // logicals
1791
1792 a_is_Charge = data1.is_Charge;
1793 a_is_Dipole = data1.is_Dipole;
1794 a_is_Quadrupole = data1.is_Quadrupole;
1795 a_is_Fluctuating = data1.is_Fluctuating;
1796
1797 b_is_Charge = data2.is_Charge;
1798 b_is_Dipole = data2.is_Dipole;
1799 b_is_Quadrupole = data2.is_Quadrupole;
1800 b_is_Fluctuating = data2.is_Fluctuating;
1801
1802 // Obtain all of the required radial function values from the
1803 // spline structures:
1804
1805 if (a_is_Charge || b_is_Charge) { v01 = v01s->getValueAt(rij); }
1806 if (a_is_Dipole || b_is_Dipole) {
1807 v11 = v11s->getValueAt(rij);
1808 v11or = ri * v11;
1809 }
1810 if (a_is_Quadrupole || b_is_Quadrupole) {
1811 v21 = v21s->getValueAt(rij);
1812 v22 = v22s->getValueAt(rij);
1813 v22or = ri * v22;
1814 }
1815
1816 if (a_is_Charge) {
1817 C_a = data1.fixedCharge;
1818
1819 if (a_is_Fluctuating) { C_a += a1->getFlucQPos(); }
1820
1821 Pb += C_a * pre11_ * v01;
1822 }
1823
1824 if (a_is_Dipole) {
1825 D_a = a1->getDipole() * mPoleConverter;
1826 rdDa = dot(rhat, D_a);
1827 Pb += pre12_ * v11 * rdDa;
1828 }
1829
1830 if (a_is_Quadrupole) {
1831 Q_a = a1->getQuadrupole() * mPoleConverter;
1832 trQa = Q_a.trace();
1833 Qar = Q_a * rhat;
1834 rdQar = dot(rhat, Qar);
1835 Pb += pre14_ * (v21 * trQa + v22 * rdQar);
1836 }
1837
1838 if (b_is_Charge) {
1839 C_b = data2.fixedCharge;
1840
1841 if (b_is_Fluctuating) C_b += a2->getFlucQPos();
1842
1843 Pa += C_b * pre11_ * v01;
1844 }
1845
1846 if (b_is_Dipole) {
1847 D_b = a2->getDipole() * mPoleConverter;
1848 rdDb = dot(rhat, D_b);
1849 Pa += pre12_ * v11 * rdDb;
1850 }
1851
1852 if (b_is_Quadrupole) {
1853 Q_a = a2->getQuadrupole() * mPoleConverter;
1854 trQb = Q_b.trace();
1855 Qbr = Q_b * rhat;
1856 rdQbr = dot(rhat, Qbr);
1857 Pa += pre14_ * (v21 * trQb + v22 * rdQbr);
1858 }
1859
1860 spot1 = Pa;
1861 spot2 = Pb;
1862 }
1863
1864 RealType Electrostatic::getFieldFunction(RealType r) {
1865 if (!initialized_) { initialize(); }
1866 RealType v01, dv01;
1867 v01s->getValueAndDerivativeAt(r, v01, dv01);
1868 return dv01 * pre11_;
1869 }
1870
1871} // namespace OpenMD
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
unsigned int getNAtoms()
Returns the number of local atoms.
Definition SimInfo.hpp:172
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
Definition Snapshot.cpp:337
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Real trace() const
Returns the trace of this matrix.
Vector< Real, Dim > diagonals() const
Returns a column vector that contains the elements from the diagonal of m in the order R(0) = m(0,...
void zero()
Zeros out the values in this vector in place.
Definition Vector.hpp:194
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
Real lengthSquare()
Returns the squared length of this vector.
Definition Vector.hpp:399
void Vdiv(const Vector< Real, Dim > &v1, const Vector< Real, Dim > &v2)
Sets the elements of this vector to the division of elements of two other vectors.
Definition Vector.hpp:336
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
ElectrostaticSummationMethod
@ esm_EWALD_SPME
SPME Ewald methods aren't supported yet.
@ esm_EWALD_PME
PME Ewald methods aren't supported yet.
Vector< Real, Row > mCross(const RectMatrix< Real, Row, Col > &t1, const RectMatrix< Real, Row, Col > &t2)
Returns the vector (cross) product of two matrices.
@ ELECTROSTATIC_FAMILY
Coulombic and point-multipole interactions.
The InteractionData struct.
Vector3d d
interatomic vector (already wrapped into box)
RealType sPot2
site potential on second atom
Vector3d D_1
dipole vector of first atom
RealType skippedCharge1
charge skipped on atom1 in pairwise interaction loop with atom2
RealType dVdFQ1
fluctuating charge force on atom1
RealType dVdFQ2
fluctuating charge force on atom2
RealType electroMult
multiplier for electrostatic interactions
int atid1
atomType ident for atom 1
potVec selePot
potential energies of the selected sites
RealType flucQ2
fluctuating charge on atom2
bool excluded
is this excluded from direct pairwise interactions? (some indirect interactions may still apply)
RealType sw
switching function value at rij
Vector3d eField1
electric field on first atom
bool isSelected
one of the particles has been selected for selection potential
bool doElectricField
should we bother with the electric field?
RealType flucQ1
fluctuating charge on atom1
Mat3x3d Q_2
quadrupole tensor of first atom
Vector3d eField2
electric field on second atom
Vector3d t1
torque on first atom
Vector3d t2
torque on second atom
RealType vpair
pair potential
RealType sPot1
site potential on first atom
int atid2
atomType ident for atom 2
bool doSitePotential
should we bother with electrostatic site potential
Vector3d D_2
dipole vector of first atom
RealType skippedCharge2
charge skipped on atom2 in pairwise interaction loop with atom1
potVec excludedPot
potential energy excluded from the overall calculation
Mat3x3d Q_1
quadrupole tensor of first atom
Vector3d f1
force between the two atoms
RealType rij
interatomic separation
The SelfData struct.
potVec selePot
potential energy of the selected site
RealType flucQfrc
fluctuating charge derivative
potVec selfPot
total potential (including embedding energy)
RealType particlePot
contribution to potential from this particle
bool isSelected
this site has been selected for selection potential
RealType flucQ
current value of atom's fluctuating charge
RealType skippedCharge
charge skipped in pairwise interaction loop
bool doParticlePot
should we bother with the particle pot?
int atid
atomType ident for the atom