| 51 |
|
|
| 52 |
|
namespace OpenMD { |
| 53 |
|
|
| 54 |
< |
bool Morse::initialized_ = false; |
| 55 |
< |
bool Morse::shiftedPot_ = false; |
| 56 |
< |
bool Morse::shiftedFrc_ = false; |
| 57 |
< |
ForceField* Morse::forceField_ = NULL; |
| 58 |
< |
map<int, AtomType*> Morse::MorseMap; |
| 59 |
< |
map<pair<AtomType*, AtomType*>, MorseInteractionData> Morse::MixingMap; |
| 60 |
< |
map<string, MorseInteractionType> Morse::stringToEnumMap_; |
| 54 |
> |
Morse::Morse() : name_("Morse"), initialized_(false), forceField_(NULL), |
| 55 |
> |
shiftedPot_(false), shiftedFrc_(false) {} |
| 56 |
|
|
| 62 |
– |
Morse* Morse::_instance = NULL; |
| 63 |
– |
|
| 64 |
– |
Morse* Morse::Instance() { |
| 65 |
– |
if (!_instance) { |
| 66 |
– |
_instance = new Morse(); |
| 67 |
– |
} |
| 68 |
– |
return _instance; |
| 69 |
– |
} |
| 70 |
– |
|
| 57 |
|
void Morse::initialize() { |
| 58 |
|
|
| 59 |
|
stringToEnumMap_["shiftedMorse"] = shiftedMorse; |
| 126 |
|
RealType De, RealType Re, RealType beta, |
| 127 |
|
MorseInteractionType mit) { |
| 128 |
|
|
| 143 |
– |
AtomTypeProperties atp = atype1->getATP(); |
| 144 |
– |
MorseMap.insert( pair<int, AtomType*>(atp.ident, atype1) ); |
| 145 |
– |
|
| 146 |
– |
atp = atype2->getATP(); |
| 147 |
– |
MorseMap.insert( pair<int, AtomType*>(atp.ident, atype2) ); |
| 148 |
– |
|
| 129 |
|
MorseInteractionData mixer; |
| 130 |
|
mixer.De = De; |
| 131 |
|
mixer.Re = Re; |
| 142 |
|
} |
| 143 |
|
} |
| 144 |
|
|
| 145 |
< |
void Morse::calcForce(AtomType* at1, AtomType* at2, Vector3d d, |
| 166 |
< |
RealType r, RealType r2, RealType rcut, RealType sw, |
| 167 |
< |
RealType vdwMult, RealType &vpair, RealType &pot, |
| 168 |
< |
Vector3d &f1) { |
| 145 |
> |
void Morse::calcForce(InteractionData idat) { |
| 146 |
|
|
| 147 |
|
if (!initialized_) initialize(); |
| 148 |
|
|
| 151 |
|
RealType myDeriv = 0.0; |
| 152 |
|
RealType myDerivC = 0.0; |
| 153 |
|
|
| 154 |
< |
pair<AtomType*, AtomType*> key = make_pair(at1, at2); |
| 154 |
> |
pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2); |
| 155 |
|
MorseInteractionData mixer = MixingMap[key]; |
| 156 |
|
|
| 157 |
|
RealType De = mixer.De; |
| 161 |
|
|
| 162 |
|
// V(r) = D_e exp(-a(r-re)(exp(-a(r-re))-2) |
| 163 |
|
|
| 164 |
< |
RealType expt = -beta*(r - Re); |
| 164 |
> |
RealType expt = -beta*(idat.rij - Re); |
| 165 |
|
RealType expfnc = exp(expt); |
| 166 |
|
RealType expfnc2 = expfnc*expfnc; |
| 167 |
|
|
| 170 |
|
RealType expfnc2C = 0.0; |
| 171 |
|
|
| 172 |
|
if (Morse::shiftedPot_ || Morse::shiftedFrc_) { |
| 173 |
< |
exptC = -beta*(rcut - Re); |
| 173 |
> |
exptC = -beta*(idat.rcut - Re); |
| 174 |
|
expfncC = exp(exptC); |
| 175 |
|
expfnc2C = expfncC*expfncC; |
| 176 |
|
} |
| 188 |
|
} else if (Morse::shiftedFrc_) { |
| 189 |
|
myPotC = De * (expfnc2C - 2.0 * expfncC); |
| 190 |
|
myDerivC = 2.0 * De * beta * (expfnc2C - expfnc2C); |
| 191 |
< |
myPotC += myDerivC * (r - rcut); |
| 191 |
> |
myPotC += myDerivC * (idat.rij - idat.rcut); |
| 192 |
|
} else { |
| 193 |
|
myPotC = 0.0; |
| 194 |
|
myDerivC = 0.0; |
| 207 |
|
} else if (Morse::shiftedFrc_) { |
| 208 |
|
myPotC = De * expfnc2C; |
| 209 |
|
myDerivC = -2.0 * De * beta * expfnc2C; |
| 210 |
< |
myPotC += myDerivC * (r - rcut); |
| 210 |
> |
myPotC += myDerivC * (idat.rij - idat.rcut); |
| 211 |
|
} else { |
| 212 |
|
myPotC = 0.0; |
| 213 |
|
myDerivC = 0.0; |
| 217 |
|
} |
| 218 |
|
} |
| 219 |
|
|
| 220 |
< |
RealType pot_temp = vdwMult * (myPot - myPotC); |
| 221 |
< |
vpair += pot_temp; |
| 220 |
> |
RealType pot_temp = idat.vdwMult * (myPot - myPotC); |
| 221 |
> |
idat.vpair += pot_temp; |
| 222 |
|
|
| 223 |
< |
RealType dudr = sw * vdwMult * (myDeriv - myDerivC); |
| 223 |
> |
RealType dudr = idat.sw * idat.vdwMult * (myDeriv - myDerivC); |
| 224 |
|
|
| 225 |
< |
pot += sw * pot_temp; |
| 226 |
< |
f1 = d * dudr / r; |
| 225 |
> |
idat.pot += idat.sw * pot_temp; |
| 226 |
> |
idat.f1 = idat.d * dudr / idat.rij; |
| 227 |
|
|
| 228 |
|
return; |
| 229 |
|
|
| 230 |
|
} |
| 254 |
– |
|
| 255 |
– |
void Morse::do_morse_pair(int *atid1, int *atid2, RealType *d, |
| 256 |
– |
RealType *rij, RealType *r2, RealType *rcut, |
| 257 |
– |
RealType *sw, RealType *vdwMult, |
| 258 |
– |
RealType *vpair, RealType *pot, RealType *f1) { |
| 259 |
– |
|
| 260 |
– |
if (!initialized_) initialize(); |
| 231 |
|
|
| 262 |
– |
AtomType* atype1 = MorseMap[*atid1]; |
| 263 |
– |
AtomType* atype2 = MorseMap[*atid2]; |
| 264 |
– |
|
| 265 |
– |
Vector3d disp(d[0], d[1], d[2]); |
| 266 |
– |
Vector3d frc(f1[0], f1[1], f1[2]); |
| 267 |
– |
|
| 268 |
– |
calcForce(atype1, atype2, disp, *rij, *r2, *rcut, *sw, *vdwMult, *vpair, |
| 269 |
– |
*pot, frc); |
| 270 |
– |
|
| 271 |
– |
f1[0] = frc.x(); |
| 272 |
– |
f1[1] = frc.y(); |
| 273 |
– |
f1[2] = frc.z(); |
| 274 |
– |
|
| 275 |
– |
return; |
| 276 |
– |
} |
| 277 |
– |
|
| 278 |
– |
void Morse::setMorseDefaultCutoff(RealType *thisRcut, int *sP, int *sF) { |
| 279 |
– |
shiftedPot_ = (bool)(*sP); |
| 280 |
– |
shiftedFrc_ = (bool)(*sF); |
| 281 |
– |
} |
| 232 |
|
} |
| 233 |
|
|
| 284 |
– |
extern "C" { |
| 285 |
– |
|
| 286 |
– |
#define fortranSetMorseCutoff FC_FUNC(setmorsedefaultcutoff, SETMORSEDEFAULTCUTOFF) |
| 287 |
– |
#define fortranDoMorsePair FC_FUNC(do_morse_pair, DO_MORSE_PAIR) |
| 288 |
– |
|
| 289 |
– |
void fortranSetMorseCutoff(RealType *rcut, int *shiftedPot, int *shiftedFrc) { |
| 290 |
– |
return OpenMD::Morse::Instance()->setMorseDefaultCutoff(rcut, shiftedPot, |
| 291 |
– |
shiftedFrc); |
| 292 |
– |
} |
| 293 |
– |
void fortranDoMorsePair(int *atid1, int *atid2, RealType *d, RealType *rij, |
| 294 |
– |
RealType *r2, RealType *rcut, RealType *sw, |
| 295 |
– |
RealType *vdwMult, RealType* vpair, RealType* pot, |
| 296 |
– |
RealType *f1){ |
| 297 |
– |
|
| 298 |
– |
return OpenMD::Morse::Instance()->do_morse_pair(atid1, atid2, d, rij, r2, |
| 299 |
– |
rcut, sw, vdwMult, vpair, |
| 300 |
– |
pot, f1); |
| 301 |
– |
} |
| 302 |
– |
} |
