OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Buckingham.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/Buckingham.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <cstring>
50
51#include "types/BuckinghamInteractionType.hpp"
52#include "utils/simError.h"
53
54using namespace std;
55
56namespace OpenMD {
57
58 Buckingham::Buckingham() :
59 initialized_(false), forceField_(NULL), name_("Buckingham") {}
60
61 void Buckingham::initialize() {
62 Btypes.clear();
63 Btids.clear();
64 MixingMap.clear();
65 Btids.resize(forceField_->getNAtomType(), -1);
66
67 ForceField::NonBondedInteractionTypeContainer* nbiTypes =
68 forceField_->getNonBondedInteractionTypes();
69 ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
70 ForceField::NonBondedInteractionTypeContainer::KeyType keys;
71 NonBondedInteractionType* nbt;
72 int btid1, btid2;
73
74 for (nbt = nbiTypes->beginType(j); nbt != NULL;
75 nbt = nbiTypes->nextType(j)) {
76 if (nbt->isBuckingham()) {
77 keys = nbiTypes->getKeys(j);
78 AtomType* at1 = forceField_->getAtomType(keys[0]);
79 if (at1 == NULL) {
80 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
81 "Buckingham::initialize could not find AtomType %s\n"
82 "\tto for for %s - %s interaction.\n",
83 keys[0].c_str(), keys[0].c_str(), keys[1].c_str());
84 painCave.severity = OPENMD_ERROR;
85 painCave.isFatal = 1;
86 simError();
87 }
88
89 AtomType* at2 = forceField_->getAtomType(keys[1]);
90 if (at2 == NULL) {
91 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
92 "Buckingham::initialize could not find AtomType %s\n"
93 "\tfor %s - %s nonbonded interaction.\n",
94 keys[1].c_str(), keys[0].c_str(), keys[1].c_str());
95 painCave.severity = OPENMD_ERROR;
96 painCave.isFatal = 1;
97 simError();
98 }
99
100 int atid1 = at1->getIdent();
101 if (Btids[atid1] == -1) {
102 btid1 = Btypes.size();
103 Btypes.insert(atid1);
104 Btids[atid1] = btid1;
105 }
106 int atid2 = at2->getIdent();
107 if (Btids[atid2] == -1) {
108 btid2 = Btypes.size();
109 Btypes.insert(atid2);
110 Btids[atid2] = btid2;
111 }
112
113 BuckinghamInteractionType* bit =
114 dynamic_cast<BuckinghamInteractionType*>(nbt);
115
116 if (bit == NULL) {
117 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
118 "Buckingham::initialize could not convert "
119 "NonBondedInteractionType\n"
120 "\tto BuckinghamInteractionType for %s - %s interaction.\n",
121 at1->getName().c_str(), at2->getName().c_str());
122 painCave.severity = OPENMD_ERROR;
123 painCave.isFatal = 1;
124 simError();
125 }
126
127 RealType A = bit->getA();
128 RealType B = bit->getB();
129 RealType C = bit->getC();
130
131 BuckinghamType variant = bit->getInteractionType();
132 addExplicitInteraction(at1, at2, A, B, C, variant);
133 }
134 }
135 initialized_ = true;
136 }
137
138 void Buckingham::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
139 RealType A, RealType B, RealType C,
140 BuckinghamType bt) {
141 BuckinghamInteractionData mixer;
142 mixer.A = A;
143 mixer.B = B;
144 mixer.C = C;
145 mixer.variant = bt;
146
147 int btid1 = Btids[atype1->getIdent()];
148 int btid2 = Btids[atype2->getIdent()];
149 int nB = Btypes.size();
150
151 MixingMap.resize(nB);
152 MixingMap[btid1].resize(nB);
153
154 MixingMap[btid1][btid2] = mixer;
155 if (btid2 != btid1) {
156 MixingMap[btid2].resize(nB);
157 MixingMap[btid2][btid1] = mixer;
158 }
159 }
160
161 void Buckingham::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
162 RealType A, RealType B, RealType C,
163 RealType sigma, RealType epsilon,
164 BuckinghamType bt) {
165 BuckinghamInteractionData mixer;
166 mixer.A = A;
167 mixer.B = B;
168 mixer.C = C;
169 mixer.sigma = sigma;
170 mixer.epsilon = epsilon;
171 mixer.variant = bt;
172
173 int btid1 = Btids[atype1->getIdent()];
174 int btid2 = Btids[atype2->getIdent()];
175 int nB = Btypes.size();
176
177 MixingMap.resize(nB);
178 MixingMap[btid1].resize(nB);
179
180 MixingMap[btid1][btid2] = mixer;
181 if (btid2 != btid1) {
182 MixingMap[btid2].resize(nB);
183 MixingMap[btid2][btid1] = mixer;
184 }
185 }
186
187 void Buckingham::calcForce(InteractionData& idat) {
188 if (!initialized_) initialize();
189
190 BuckinghamInteractionData& mixer =
191 MixingMap[Btids[idat.atid1]][Btids[idat.atid2]];
192
193 RealType myPot = 0.0;
194 RealType myPotC = 0.0;
195 RealType myDeriv = 0.0;
196 RealType myDerivC = 0.0;
197
198 RealType A = mixer.A;
199 RealType B = mixer.B;
200 RealType C = mixer.C;
201 RealType sigma = mixer.sigma;
202 RealType epsilon = mixer.epsilon;
203 BuckinghamType variant = mixer.variant;
204
205 RealType expt = -B * idat.rij;
206 RealType expfnc = exp(expt);
207 RealType fnc6 = 1.0 / pow(idat.rij, 6);
208 RealType fnc7 = fnc6 / idat.rij;
209
210 RealType exptC = 0.0;
211 RealType expfncC = 0.0;
212 RealType fnc6C = 0.0;
213 RealType fnc7C = 0.0;
214
215 if (idat.shiftedPot || idat.shiftedForce) {
216 exptC = -B * idat.rcut;
217 expfncC = exp(exptC);
218 fnc6C = 1.0 / pow(idat.rcut, 6);
219 fnc7C = fnc6C / idat.rcut;
220 }
221
222 switch (variant) {
223 case btTraditional: {
224 // V(r) = A exp(-B*r) - C/r^6
225 myPot = A * expfnc - C * fnc6;
226 myDeriv = -A * B * expfnc + C * fnc7;
227
228 if (idat.shiftedPot) {
229 myPotC = A * expfncC - C * fnc6C;
230 myDerivC = 0.0;
231 } else if (idat.shiftedForce) {
232 myPotC = A * expfncC - C * fnc6C;
233 myDerivC = -A * B * expfncC + C * fnc7C;
234 myPotC += myDerivC * (idat.rij - idat.rcut);
235 } else {
236 myPotC = 0.0;
237 myDerivC = 0.0;
238 }
239 break;
240 }
241 case btModified: {
242 RealType s6 = pow(sigma, 6);
243 RealType s7 = pow(sigma, 7);
244 RealType fnc30 = pow(sigma / idat.rij, 30);
245 RealType fnc31 = fnc30 * sigma / idat.rij;
246 RealType fnc30C = 0.0;
247 RealType fnc31C = 0.0;
248
249 if (idat.shiftedPot || idat.shiftedForce) {
250 fnc30C = pow(sigma / idat.rcut, 30);
251 fnc31C = fnc30C * sigma / idat.rcut;
252 }
253
254 // V(r) = A exp(-B*r) - C/r^6 + 4 epsilon ((sigma/r)^30 - (sigma/r)^6)
255 myPot = A * expfnc - C * fnc6 + 4.0 * epsilon * (fnc30 - s6 * fnc6);
256 myDeriv = -A * B * expfnc + C * fnc7 +
257 4.0 * epsilon * (-30.0 * fnc31 + 6.0 * s7 * fnc7) / sigma;
258
259 if (idat.shiftedPot) {
260 myPotC =
261 A * expfncC - C * fnc6C + 4.0 * epsilon * (fnc30C - s6 * fnc6C);
262 myDerivC = 0.0;
263 } else if (idat.shiftedForce) {
264 myPotC =
265 A * expfncC - C * fnc6C + 4.0 * epsilon * (fnc30C - s6 * fnc6C);
266 myDeriv = -A * B * expfncC + C * fnc7C +
267 4.0 * epsilon * (-30.0 * fnc31C + 6.0 * s7 * fnc7C) / sigma;
268 myPotC += myDerivC * (idat.rij - idat.rcut);
269 } else {
270 myPotC = 0.0;
271 myDerivC = 0.0;
272 }
273
274 break;
275 }
276 case btUnknown: {
277 // don't know what to do so don't do anything
278 break;
279 }
280 }
281
282 RealType pot_temp = idat.vdwMult * (myPot - myPotC);
283 idat.vpair += pot_temp;
284
285 RealType dudr = idat.sw * idat.vdwMult * (myDeriv - myDerivC);
286
287 idat.pot[VANDERWAALS_FAMILY] += idat.sw * pot_temp;
288 if (idat.isSelected) idat.selePot[VANDERWAALS_FAMILY] += idat.sw * pot_temp;
289
290 idat.f1 += idat.d * dudr / idat.rij;
291
292 return;
293 }
294
295 RealType Buckingham::getSuggestedCutoffRadius(
296 pair<AtomType*, AtomType*> atypes) {
297 if (!initialized_) initialize();
298
299 int atid1 = atypes.first->getIdent();
300 int atid2 = atypes.second->getIdent();
301 int btid1 = Btids[atid1];
302 int btid2 = Btids[atid2];
303
304 if (btid1 == -1 || btid2 == -1)
305 return 0.0;
306 else {
307 // Uncomment if we ever want to query the simulated atoms types
308 // for a suggested cutoff:
309 //
310 // BuckinghamInteractionData mixer = MixingMap[btid1][btid2];
311 //
312 // suggested cutoff for most implementations of the BKS potential are
313 // around 1 nm (10 angstroms):
314 return 10.0;
315 }
316 }
317} // namespace OpenMD
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
@ VANDERWAALS_FAMILY
Long-range dispersion and short-range pauli repulsion.