OpenMD 3.1
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"
62#include "primitives/Molecule.hpp"
63#include "types/FixedChargeAdapter.hpp"
64#include "types/FluctuatingChargeAdapter.hpp"
65#include "types/MultipoleAdapter.hpp"
66#include "types/NonBondedInteractionType.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
692 FixedChargeAdapter fca = FixedChargeAdapter(atomType);
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
711 FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atomType);
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(
811 ElectrostaticSummationMethod esm) {
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;
1313 ElectrostaticAtomData data;
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();
1345 info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
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();
1391 info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
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
1436 Mat3x3d hmat = info_->getSnapshotManager()->getCurrentSnapshot()->getHmat();
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();
1502 info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(r);
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;
1577 ElectrostaticAtomData data;
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();
1777 info_->getSnapshotManager()->getCurrentSnapshot()->wrapVector(d);
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
OpenMD::Atom::getAtomType
AtomType * getAtomType()
Returns the AtomType of this Atom.
Definition A.hpp:93
OpenMD::AtomType
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
OpenMD::SimInfo::beginMolecule
Molecule * beginMolecule(MoleculeIterator &i)
Returns the first molecule in this SimInfo and intialize the iterator.
Definition SimInfo.cpp:240
OpenMD::SimInfo::getNAtoms
unsigned int getNAtoms()
Returns the number of local atoms.
Definition SimInfo.hpp:172
OpenMD::SimInfo::nextMolecule
Molecule * nextMolecule(MoleculeIterator &i)
Returns the next avaliable Molecule based on the iterator.
Definition SimInfo.cpp:245
OpenMD::SimInfo::getSnapshotManager
SnapshotManager * getSnapshotManager()
Returns the snapshot manager.
Definition SimInfo.hpp:248
OpenMD::Snapshot::getHmat
Mat3x3d getHmat()
Returns the H-Matrix.
Definition Snapshot.cpp:214
OpenMD::Snapshot::wrapVector
void wrapVector(Vector3d &v)
Wrapping the vector according to periodic boundary condition.
Definition Snapshot.cpp:337
OpenMD::SnapshotManager::getCurrentSnapshot
Snapshot * getCurrentSnapshot()
Returns the pointer of current snapshot.
Definition SnapshotManager.hpp:89
OpenMD::SquareMatrix3::trace
Real trace() const
Returns the trace of this matrix.
Definition SquareMatrix3.hpp:412
OpenMD::SquareMatrix::diagonals
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,...
Definition SquareMatrix.hpp:212
OpenMD::StuntDouble::getQuadrupole
Mat3x3d getQuadrupole()
Returns the current quadrupole tensor of this stuntDouble.
Definition StuntDouble.hpp:627
OpenMD::StuntDouble::getPos
Vector3d getPos()
Returns the current position of this stuntDouble.
Definition StuntDouble.hpp:196
OpenMD::StuntDouble::getFlucQPos
RealType getFlucQPos()
Returns the current fluctuating charge of this stuntDouble.
Definition StuntDouble.hpp:930
OpenMD::StuntDouble::getDipole
Vector3d getDipole()
Returns the current dipole vector of this stuntDouble.
Definition StuntDouble.hpp:598
OpenMD::Vector3::z
Real & z()
Returns reference of the third element of Vector3.
Definition Vector3.hpp:120
OpenMD::Vector3::x
Real & x()
Returns reference of the first element of Vector3.
Definition Vector3.hpp:96
OpenMD::Vector3::y
Real & y()
Returns reference of the second element of Vector3.
Definition Vector3.hpp:108
OpenMD::Vector::zero
void zero()
Zeros out the values in this vector in place.
Definition Vector.hpp:194
OpenMD::Vector::length
Real length()
Returns the length of this vector.
Definition Vector.hpp:393
OpenMD::Vector::Vmul
void Vmul(const Vector< Real, Dim > &v1, const Vector< Real, Dim > &v2)
Sets the elements of this vector to the multiplication of elements of two other vectors.
Definition Vector.hpp:284
OpenMD::Vector::lengthSquare
Real lengthSquare()
Returns the squared length of this vector.
Definition Vector.hpp:399
OpenMD::Vector::Vdiv
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
OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
Definition ActionCorrFunc.cpp:60
OpenMD::cross
Vector3< Real > cross(const Vector3< Real > &v1, const Vector3< Real > &v2)
Returns the cross product of two Vectors.
Definition Vector3.hpp:136
OpenMD::dot
Real dot(const DynamicVector< Real > &v1, const DynamicVector< Real > &v2)
Returns the dot product of two DynamicVectors.
Definition DynamicVector.hpp:438
OpenMD::ElectrostaticSummationMethod
ElectrostaticSummationMethod
Definition Electrostatic.hpp:74
OpenMD::esm_EWALD_SPME
@ esm_EWALD_SPME
SPME Ewald methods aren't supported yet.
Definition Electrostatic.hpp:83
OpenMD::esm_EWALD_PME
@ esm_EWALD_PME
PME Ewald methods aren't supported yet.
Definition Electrostatic.hpp:82
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:636
OpenMD::ELECTROSTATIC_FAMILY
@ ELECTROSTATIC_FAMILY
Coulombic and point-multipole interactions.
Definition NonBondedInteraction.hpp:66
OpenMD::InteractionData
The InteractionData struct.
Definition NonBondedInteraction.hpp:103
OpenMD::InteractionData::d
Vector3d d
interatomic vector (already wrapped into box)
Definition NonBondedInteraction.hpp:106
OpenMD::InteractionData::sPot2
RealType sPot2
site potential on second atom
Definition NonBondedInteraction.hpp:162
OpenMD::InteractionData::D_1
Vector3d D_1
dipole vector of first atom
Definition NonBondedInteraction.hpp:139
OpenMD::InteractionData::skippedCharge1
RealType skippedCharge1
charge skipped on atom1 in pairwise interaction loop with atom2
Definition NonBondedInteraction.hpp:157
OpenMD::InteractionData::dVdFQ1
RealType dVdFQ1
fluctuating charge force on atom1
Definition NonBondedInteraction.hpp:153
OpenMD::InteractionData::dVdFQ2
RealType dVdFQ2
fluctuating charge force on atom2
Definition NonBondedInteraction.hpp:154
OpenMD::InteractionData::electroMult
RealType electroMult
multiplier for electrostatic interactions
Definition NonBondedInteraction.hpp:121
OpenMD::InteractionData::atid1
int atid1
atomType ident for atom 1
Definition NonBondedInteraction.hpp:104
OpenMD::InteractionData::selePot
potVec selePot
potential energies of the selected sites
Definition NonBondedInteraction.hpp:133
OpenMD::InteractionData::flucQ2
RealType flucQ2
fluctuating charge on atom2
Definition NonBondedInteraction.hpp:152
OpenMD::InteractionData::excluded
bool excluded
is this excluded from direct pairwise interactions? (some indirect interactions may still apply)
Definition NonBondedInteraction.hpp:115
OpenMD::InteractionData::sw
RealType sw
switching function value at rij
Definition NonBondedInteraction.hpp:113
OpenMD::InteractionData::eField1
Vector3d eField1
electric field on first atom
Definition NonBondedInteraction.hpp:155
OpenMD::InteractionData::isSelected
bool isSelected
one of the particles has been selected for selection potential
Definition NonBondedInteraction.hpp:131
OpenMD::InteractionData::doElectricField
bool doElectricField
should we bother with the electric field?
Definition NonBondedInteraction.hpp:127
OpenMD::InteractionData::flucQ1
RealType flucQ1
fluctuating charge on atom1
Definition NonBondedInteraction.hpp:151
OpenMD::InteractionData::Q_2
Mat3x3d Q_2
quadrupole tensor of first atom
Definition NonBondedInteraction.hpp:142
OpenMD::InteractionData::eField2
Vector3d eField2
electric field on second atom
Definition NonBondedInteraction.hpp:156
OpenMD::InteractionData::t1
Vector3d t1
torque on first atom
Definition NonBondedInteraction.hpp:143
OpenMD::InteractionData::t2
Vector3d t2
torque on second atom
Definition NonBondedInteraction.hpp:144
OpenMD::InteractionData::vpair
RealType vpair
pair potential
Definition NonBondedInteraction.hpp:125
OpenMD::InteractionData::sPot1
RealType sPot1
site potential on first atom
Definition NonBondedInteraction.hpp:161
OpenMD::InteractionData::atid2
int atid2
atomType ident for atom 2
Definition NonBondedInteraction.hpp:105
OpenMD::InteractionData::doSitePotential
bool doSitePotential
should we bother with electrostatic site potential
Definition NonBondedInteraction.hpp:129
OpenMD::InteractionData::D_2
Vector3d D_2
dipole vector of first atom
Definition NonBondedInteraction.hpp:140
OpenMD::InteractionData::skippedCharge2
RealType skippedCharge2
charge skipped on atom2 in pairwise interaction loop with atom1
Definition NonBondedInteraction.hpp:159
OpenMD::InteractionData::excludedPot
potVec excludedPot
potential energy excluded from the overall calculation
Definition NonBondedInteraction.hpp:123
OpenMD::InteractionData::Q_1
Mat3x3d Q_1
quadrupole tensor of first atom
Definition NonBondedInteraction.hpp:141
OpenMD::InteractionData::f1
Vector3d f1
force between the two atoms
Definition NonBondedInteraction.hpp:136
OpenMD::InteractionData::rij
RealType rij
interatomic separation
Definition NonBondedInteraction.hpp:107
OpenMD::SelfData
The SelfData struct.
Definition NonBondedInteraction.hpp:176
OpenMD::SelfData::selePot
potVec selePot
potential energy of the selected site
Definition NonBondedInteraction.hpp:191
OpenMD::SelfData::flucQfrc
RealType flucQfrc
fluctuating charge derivative
Definition NonBondedInteraction.hpp:188
OpenMD::SelfData::selfPot
potVec selfPot
total potential (including embedding energy)
Definition NonBondedInteraction.hpp:180
OpenMD::SelfData::particlePot
RealType particlePot
contribution to potential from this particle
Definition NonBondedInteraction.hpp:183
OpenMD::SelfData::isSelected
bool isSelected
this site has been selected for selection potential
Definition NonBondedInteraction.hpp:189
OpenMD::SelfData::flucQ
RealType flucQ
current value of atom's fluctuating charge
Definition NonBondedInteraction.hpp:187
OpenMD::SelfData::skippedCharge
RealType skippedCharge
charge skipped in pairwise interaction loop
Definition NonBondedInteraction.hpp:179
OpenMD::SelfData::doParticlePot
bool doParticlePot
should we bother with the particle pot?
Definition NonBondedInteraction.hpp:181
OpenMD::SelfData::atid
int atid
atomType ident for the atom
Definition NonBondedInteraction.hpp:177