OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
Mie.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/Mie.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <cstring>
50
51#include "types/MieInteractionType.hpp"
52#include "utils/simError.h"
53
54using namespace std;
55
56namespace OpenMD {
57
58 Mie::Mie() : initialized_(false), forceField_(NULL), name_("Mie") {}
59
60 void Mie::initialize() {
61 MieTypes.clear();
62 MieTids.clear();
63 MixingMap.clear();
64 MieTids.resize(forceField_->getNAtomType(), -1);
65
66 ForceField::NonBondedInteractionTypeContainer* nbiTypes =
67 forceField_->getNonBondedInteractionTypes();
68 ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
69 ForceField::NonBondedInteractionTypeContainer::KeyType keys;
70 NonBondedInteractionType* nbt;
71 int mietid1, mietid2;
72
73 for (nbt = nbiTypes->beginType(j); nbt != NULL;
74 nbt = nbiTypes->nextType(j)) {
75 if (nbt->isMie()) {
76 keys = nbiTypes->getKeys(j);
77 AtomType* at1 = forceField_->getAtomType(keys[0]);
78 if (at1 == NULL) {
79 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
80 "Mie::initialize could not find AtomType %s\n"
81 "\tto for for %s - %s interaction.\n",
82 keys[0].c_str(), keys[0].c_str(), keys[1].c_str());
83 painCave.severity = OPENMD_ERROR;
84 painCave.isFatal = 1;
85 simError();
86 }
87
88 AtomType* at2 = forceField_->getAtomType(keys[1]);
89 if (at2 == NULL) {
90 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
91 "Mie::initialize could not find AtomType %s\n"
92 "\tfor %s - %s nonbonded interaction.\n",
93 keys[1].c_str(), keys[0].c_str(), keys[1].c_str());
94 painCave.severity = OPENMD_ERROR;
95 painCave.isFatal = 1;
96 simError();
97 }
98
99 int atid1 = at1->getIdent();
100 if (MieTids[atid1] == -1) {
101 mietid1 = MieTypes.size();
102 MieTypes.insert(atid1);
103 MieTids[atid1] = mietid1;
104 }
105 int atid2 = at2->getIdent();
106 if (MieTids[atid2] == -1) {
107 mietid2 = MieTypes.size();
108 MieTypes.insert(atid2);
109 MieTids[atid2] = mietid2;
110 }
111
112 MieInteractionType* mit = dynamic_cast<MieInteractionType*>(nbt);
113 if (mit == NULL) {
114 snprintf(
115 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
116 "Mie::initialize could not convert NonBondedInteractionType\n"
117 "\tto MieInteractionType for %s - %s interaction.\n",
118 at1->getName().c_str(), at2->getName().c_str());
119 painCave.severity = OPENMD_ERROR;
120 painCave.isFatal = 1;
121 simError();
122 }
123
124 RealType sigma = mit->getSigma();
125 RealType epsilon = mit->getEpsilon();
126 int nRep = mit->getNrep();
127 int mAtt = mit->getMatt();
128
129 addExplicitInteraction(at1, at2, sigma, epsilon, nRep, mAtt);
130 }
131 }
132 initialized_ = true;
133 }
134
135 void Mie::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
136 RealType sigma, RealType epsilon, int nRep,
137 int mAtt) {
138 MieInteractionData mixer;
139 mixer.sigma = sigma;
140 mixer.epsilon = epsilon;
141 mixer.sigmai = 1.0 / mixer.sigma;
142 mixer.nRep = nRep;
143 mixer.mAtt = mAtt;
144
145 RealType n = RealType(nRep);
146 RealType m = RealType(mAtt);
147
148 mixer.nmScale = n * pow(n / m, m / (n - m)) / (n - m);
149
150 int nMie = MieTypes.size();
151 int atid1 = atype1->getIdent();
152 int atid2 = atype2->getIdent();
153 int mietid1, mietid2;
154
155 pair<set<int>::iterator, bool> ret;
156 ret = MieTypes.insert(atid1);
157 if (ret.second == false) {
158 // already had this type in the MieMap, just get the mietid:
159 mietid1 = MieTids[atid1];
160 } else {
161 // didn't already have it, so make a new one and assign it:
162 mietid1 = nMie;
163 MieTids[atid1] = nMie;
164 nMie++;
165 }
166
167 ret = MieTypes.insert(atid2);
168 if (ret.second == false) {
169 // already had this type in the MieMap, just get the mietid:
170 mietid2 = MieTids[atid2];
171 } else {
172 // didn't already have it, so make a new one and assign it:
173 mietid2 = nMie;
174 MieTids[atid2] = nMie;
175 nMie++;
176 }
177
178 MixingMap.resize(nMie);
179 MixingMap[mietid1].resize(nMie);
180
181 MixingMap[mietid1][mietid2] = mixer;
182 if (mietid2 != mietid1) {
183 MixingMap[mietid2].resize(nMie);
184 MixingMap[mietid2][mietid1] = mixer;
185 }
186 }
187
188 void Mie::calcForce(InteractionData& idat) {
189 if (!initialized_) initialize();
190
191 MieInteractionData& mixer =
192 MixingMap[MieTids[idat.atid1]][MieTids[idat.atid2]];
193 RealType sigmai = mixer.sigmai;
194 RealType epsilon = mixer.epsilon;
195 int nRep = mixer.nRep;
196 int mAtt = mixer.mAtt;
197 RealType nmScale = mixer.nmScale;
198
199 RealType ros;
200 RealType rcos;
201 RealType myPot = 0.0;
202 RealType myPotC = 0.0;
203 RealType myDeriv = 0.0;
204 RealType myDerivC = 0.0;
205
206 ros = idat.rij * sigmai;
207
208 getMieFunc(ros, nRep, mAtt, myPot, myDeriv);
209
210 if (idat.shiftedPot) {
211 rcos = idat.rcut * sigmai;
212 getMieFunc(rcos, nRep, mAtt, myPotC, myDerivC);
213 myDerivC = 0.0;
214 } else if (idat.shiftedForce) {
215 rcos = idat.rcut * sigmai;
216 getMieFunc(rcos, nRep, mAtt, myPotC, myDerivC);
217 myPotC = myPotC + myDerivC * (idat.rij - idat.rcut) * sigmai;
218 } else {
219 myPotC = 0.0;
220 myDerivC = 0.0;
221 }
222
223 RealType pot_temp = idat.vdwMult * nmScale * epsilon * (myPot - myPotC);
224 idat.vpair += pot_temp;
225
226 RealType dudr = idat.sw * idat.vdwMult * nmScale * epsilon *
227 (myDeriv - myDerivC) * sigmai;
228
229 idat.pot[VANDERWAALS_FAMILY] += idat.sw * pot_temp;
230 if (idat.isSelected) idat.selePot[VANDERWAALS_FAMILY] += idat.sw * pot_temp;
231
232 idat.f1 += idat.d * dudr / idat.rij;
233
234 return;
235 }
236
237 void Mie::getMieFunc(const RealType& r, int& n, int& m, RealType& pot,
238 RealType& deriv) {
239 RealType ri = 1.0 / r;
240 RealType rin = pow(ri, n);
241 RealType rim = pow(ri, m);
242 RealType rin1 = rin * ri;
243 RealType rim1 = rim * ri;
244
245 pot = rin - rim;
246 deriv = -n * rin1 + m * rim1;
247 return;
248 }
249
250 RealType Mie::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
251 if (!initialized_) initialize();
252
253 int atid1 = atypes.first->getIdent();
254 int atid2 = atypes.second->getIdent();
255 int mietid1 = MieTids[atid1];
256 int mietid2 = MieTids[atid2];
257
258 if (mietid1 == -1 || mietid2 == -1)
259 return 0.0;
260 else {
261 MieInteractionData mixer = MixingMap[mietid1][mietid2];
262 return 2.5 * mixer.sigma;
263 }
264 }
265} // 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.