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