OpenMD 3.0
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
EAM.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/EAM.hpp"
46
47#include <cmath>
48#include <cstdio>
49#include <cstring>
50
51#include "types/EAMInteractionType.hpp"
52#include "types/FluctuatingChargeAdapter.hpp"
53#include "utils/simError.h"
54
55namespace OpenMD {
56
57 EAM::EAM() :
58 initialized_(false), haveCutoffRadius_(false), forceField_(NULL),
59 electrostatic_(NULL), eamRcut_(0.0), mixMeth_(eamJohnson), name_("EAM") {
60 // This prefactor converts charge-charge interactions into kcal /
61 // mol assuming distances are measured in angstroms and charges
62 // are measured in electrons. Matches the value in
63 // Electrostatics.cpp
64 pre11_ = 332.0637778;
65 }
66
67 RealType EAM::fastPower(RealType x, int y) {
68 RealType temp;
69 if (y == 0) return 1;
70 temp = fastPower(x, y / 2);
71 if (y % 2 == 0)
72 return temp * temp;
73 else {
74 if (y > 0)
75 return x * temp * temp;
76 else
77 return (temp * temp) / x;
78 }
79 }
80
81 // Zhou functional form with switching functions
82
83 // RealType EAM::ZhouPhiCoreCore(RealType r, RealType re, RealType A,
84 // RealType alpha, RealType kappa) {
85 // return (A * exp(-alpha * (r / re - 1.0))) /
86 // (1.0 + fastPower(r / re - kappa, 20));
87 // }
88 //
89 // RealType EAM::ZhouPhiCoreValence(RealType r, RealType re, RealType B,
90 // RealType beta, RealType lambda) {
91 // return -(B * exp(-beta * (r / re - 1.0))) /
92 // (1.0 + fastPower(r / re - lambda, 20));
93 // }
94 //
95 // RealType EAM::ZhouPhi(RealType r, RealType re, RealType A, RealType B,
96 // RealType alpha, RealType beta, RealType kappa,
97 // RealType lambda) {
98 // return ZhouPhiCoreCore(r, re, A, alpha, kappa) +
99 // ZhouPhiCoreValence(r, re, B, beta, lambda);
100 // }
101 //
102 // RealType EAM::ZhouRho(RealType r, RealType re, RealType fe, RealType beta,
103 // RealType lambda) {
104 // return (fe * exp(-beta * (r / re - 1.0))) /
105 // (1.0 + fastPower(r / re - lambda, 20));
106 // }
107
108 // Zhou functional form without switching functions
109
110 RealType EAM::PhiCoreCore(RealType r, RealType re, RealType A, RealType alpha,
111 RealType kappa) {
112 if (mixMeth_ == eamDream2) {
113 RealType a = 0.02;
114 RealType sigma = 1;
115 return (A * exp(-alpha * (r / re - 1.0)) +
116 a * fastPower(sigma / (r + 1e-5), 30));
117 }
118
119 else if ((mixMeth_ == eamDream1) || (mixMeth_ == eamJohnson)) {
120 return (A * exp(-alpha * (r / re - 1.0))) /
121 (1.0 + fastPower(r / re - kappa, 20));
122 }
123
124 else
125 return 0;
126 }
127
128 RealType EAM::PhiCoreValence(RealType r, RealType re, RealType B,
129 RealType beta, RealType lambda) {
130 if (mixMeth_ == eamDream2) {
131 return -(B * exp(-beta * (r / re - 1.0)));
132 }
133
134 else if ((mixMeth_ == eamDream1) || (mixMeth_ == eamJohnson)) {
135 return -(B * exp(-beta * (r / re - 1.0))) /
136 (1.0 + fastPower(r / re - lambda, 20));
137 }
138
139 else
140 return 0;
141 }
142
143 RealType EAM::Phi(RealType r, RealType re, RealType A, RealType B,
144 RealType alpha, RealType beta, RealType kappa,
145 RealType lambda) {
146 return PhiCoreCore(r, re, A, alpha, kappa) +
147 PhiCoreValence(r, re, B, beta, lambda);
148 }
149
150 RealType EAM::Rho(RealType r, RealType re, RealType fe, RealType beta,
151 RealType lambda) {
152 if (mixMeth_ == eamDream2) {
153 return (fe * exp(-beta * (r / re - 1.0)));
154 }
155
156 else if ((mixMeth_ == eamDream1) || (mixMeth_ == eamJohnson)) {
157 return (fe * exp(-beta * (r / re - 1.0))) /
158 (1.0 + fastPower(r / re - lambda, 20));
159 }
160
161 else
162 return 0;
163 }
164
165 RealType EAM::gFunc(RealType q, RealType nV, RealType nM) {
166 if (q >= nV) { return 0.0; }
167 if (q <= -nM) { return (nM + nV) / nV; }
168
169 return ((q - nV) * (q - nV) * (nM + q + nV)) / (nV * nV * (nM + nV));
170 }
171
172 RealType EAM::gPrime(RealType q, RealType nV, RealType nM) {
173 if (q >= nV) { return 0.0; }
174 if (q <= -nM) { return 0.0; }
175
176 return ((q - nV) * (2 * nM + 3 * q + nV)) / (nV * nV * (nM + nV));
177 }
178
179 RealType EAM::Zhou2001Functional(RealType rho, RealType rhoe,
180 std::vector<RealType> Fn,
181 std::vector<RealType> F, RealType Fe,
182 RealType eta) {
183 RealType rhon = 0.85 * rhoe;
184 RealType rho0 = 1.15 * rhoe;
185 RealType functional(0.0);
186 RealType t;
187 if (rho < rhon) {
188 t = rho / rhon - 1.0;
189 functional = Fn.at(0) + t * (Fn.at(1) + t * (Fn.at(2) + t * Fn.at(3)));
190 } else if (rho >= rhon && rho < rho0) {
191 t = rho / rhoe - 1.0;
192 functional = F.at(0) + t * (F.at(1) + t * (F.at(2) + t * F.at(3)));
193 } else if (rho >= rho0) {
194 t = rho / rhoe;
195 functional = Fe * (1.0 - log(pow(t, eta))) * pow(t, eta);
196 }
197 return functional;
198 }
199
200 RealType EAM::Zhou2004Functional(RealType rho, RealType rhoe, RealType rhos,
201 std::vector<RealType> Fn,
202 std::vector<RealType> F, RealType Fe,
203 RealType eta, RealType rhol, RealType rhoh) {
204 RealType rhon = rhol * rhoe;
205 RealType rho0 = rhoh * rhoe;
206 RealType functional(0.0);
207 RealType t;
208 if (rho < rhon) {
209 t = rho / rhon - 1.0;
210 functional = Fn.at(0) + t * (Fn.at(1) + t * (Fn.at(2) + t * Fn.at(3)));
211 } else if (rhon <= rho && rho < rho0) {
212 t = rho / rhoe - 1.0;
213 functional = F.at(0) + t * (F.at(1) + t * (F.at(2) + t * F.at(3)));
214 } else if (rho0 <= rho) {
215 t = rho / rhos;
216 functional = Fe * (1.0 - log(pow(t, eta))) * pow(t, eta);
217 }
218 return functional;
219 }
220
221 RealType EAM::Zhou2005Functional(RealType rho, RealType rhoe, RealType rhos,
222 std::vector<RealType> Fn,
223 std::vector<RealType> F, RealType F3plus,
224 RealType F3minus, RealType Fe,
225 RealType eta) {
226 RealType rhon = 0.85 * rhoe;
227 RealType rho0 = 1.15 * rhoe;
228 RealType functional(0.0);
229 RealType t;
230 if (rho < rhon) {
231 t = rho / rhon - 1.0;
232 functional = Fn.at(0) + t * (Fn.at(1) + t * (Fn.at(2) + t * Fn.at(3)));
233 } else if (rho >= rhon && rho < rhoe) {
234 t = rho / rhoe - 1.0;
235 functional = F.at(0) + t * (F.at(1) + t * (F.at(2) + t * F3minus));
236 } else if (rho >= rhoe && rho < rho0) {
237 t = rho / rhoe - 1.0;
238 functional = F.at(0) + t * (F.at(1) + t * (F.at(2) + t * F3plus));
239 } else if (rho >= rho0) {
240 t = rho / rhos;
241 functional = Fe * (1.0 - log(pow(t, eta))) * pow(t, eta);
242 }
243 return functional;
244 }
245
246 RealType EAM::Zhou2005OxygenFunctional(
247 RealType rho, std::vector<RealType> OrhoLimits,
248 std::vector<RealType> OrhoE, std::vector<std::vector<RealType>> OF) {
249 RealType functional(0.0);
250 RealType t;
251
252 if (rho < OrhoLimits[1]) {
253 t = rho / OrhoE[0] - 1.0;
254 functional = OF[0][0] + t * (OF[0][1] + t * (OF[0][2] + t * OF[0][3]));
255 } else if (rho >= OrhoLimits[1] && rho < OrhoLimits[2]) {
256 t = rho / OrhoE[1] - 1.0;
257 functional = OF[1][0] + t * (OF[1][1] + t * OF[1][2]);
258 } else if (rho >= OrhoLimits[2] && rho < OrhoLimits[3]) {
259 t = rho / OrhoE[2] - 1.0;
260 functional = OF[2][0] + t * (OF[2][1] + t * OF[2][2]);
261 } else if (rho >= OrhoLimits[3] && rho < OrhoLimits[4]) {
262 t = rho / OrhoE[3] - 1.0;
263 functional = OF[3][0] + t * (OF[3][1] + t * OF[3][2]);
264 } else if (rho >= OrhoLimits[4]) {
265 t = rho / OrhoE[4] - 1.0;
266 functional = OF[4][0] + t * (OF[4][1] + t * OF[4][2]);
267 }
268 return functional;
269 }
270
271 RealType EAM::RoseFunctional(RealType rho, RealType rhoe, RealType F0) {
272 RealType t = rho / rhoe;
273 RealType functional(0.0);
274 if (t > 0.0) { functional = -F0 * log(t) * t; }
275 return functional;
276 }
277
278 CubicSplinePtr EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {
279 EAMAdapter ea1 = EAMAdapter(atomType1);
280 EAMAdapter ea2 = EAMAdapter(atomType2);
281
282 CubicSplinePtr cs {std::make_shared<CubicSpline>()};
283
284 RealType rha(0.0), rhb(0.0), pha(0.0), phb(0.0), phab(0.0);
285 vector<RealType> rvals;
286 vector<RealType> phivals;
287 RealType rmax(0.0), dr, r;
288 int nr;
289
290 if (ea1.getEAMType() == eamFuncfl && ea2.getEAMType() == eamFuncfl) {
291 CubicSplinePtr z1 = ea1.getZSpline();
292 CubicSplinePtr z2 = ea2.getZSpline();
293 CubicSplinePtr rho1 = ea1.getRhoSpline();
294 CubicSplinePtr rho2 = ea2.getRhoSpline();
295
296 // make the r grid:
297
298 // we need phi out to the largest value we'll encounter in the
299 // radial space;
300
301 rmax = max(rmax, ea1.getRcut());
302 rmax = max(rmax, ea1.getNr() * ea1.getDr());
303
304 rmax = max(rmax, ea2.getRcut());
305 rmax = max(rmax, ea2.getNr() * ea2.getDr());
306
307 // use the smallest dr (finest grid) to build our grid:
308
309 dr = min(ea1.getDr(), ea2.getDr());
310 nr = int(rmax / dr + 0.5);
311
312 for (int i = 0; i < nr; i++)
313 rvals.push_back(RealType(i * dr));
314
315 // construct the pair potential:
316
317 RealType za, zb;
318
319 phivals.push_back(0.0);
320
321 for (unsigned int i = 1; i < rvals.size(); i++) {
322 r = rvals[i];
323
324 // only use z(r) if we're inside this atom's cutoff radius,
325 // otherwise, we'll use zero for the charge. This effectively
326 // means that our phi grid goes out beyond the cutoff of the
327 // pair potential
328
329 rha = r <= ea1.getRcut() ? rho1->getValueAt(r) : 0.0;
330 rhb = r <= ea2.getRcut() ? rho2->getValueAt(r) : 0.0;
331
332 za = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
333 zb = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
334
335 switch (mixMeth_) {
336 case eamJohnson:
337 case eamDream1:
338 case eamDream2:
339 phab = 0.0;
340 if (rha > 0.0) {
341 pha = pre11_ * (za * za) / r;
342 phab = phab + 0.5 * (rhb / rha) * pha;
343 }
344 if (rhb > 0.0) {
345 phb = pre11_ * (zb * zb) / r;
346 phab = phab + 0.5 * (rha / rhb) * phb;
347 }
348
349 break;
350
351 case eamDaw:
352 if (r <= ea1.getRcut() && r <= ea2.getRcut()) {
353 phab = pre11_ * (za * zb) / r;
354 } else {
355 phab = 0.0;
356 }
357
358 break;
359 case eamUnknownMix:
360 default:
361
362 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
363 "EAM::getPhi hit a mixing method it doesn't know about!\n");
364 painCave.severity = OPENMD_ERROR;
365 painCave.isFatal = 1;
366 simError();
367 }
368 phivals.push_back(phab);
369 }
370 cs->addPoints(rvals, phivals);
371 } else {
372 EAMAtomData& data1 = EAMdata[EAMtids[atomType1->getIdent()]];
373 EAMAtomData& data2 = EAMdata[EAMtids[atomType2->getIdent()]];
374
375 RealType re1, A1, B1, alpha1, beta1, kappa1, lambda1;
376 RealType re2, A2, B2, alpha2, beta2, kappa2, lambda2;
377
378 re1 = ea1.getRe();
379 A1 = ea1.getA();
380 B1 = ea1.getB();
381 alpha1 = ea1.getAlpha();
382 beta1 = ea1.getBeta();
383 kappa1 = ea1.getKappa();
384 lambda1 = ea1.getLambda();
385
386 re2 = ea2.getRe();
387 A2 = ea2.getA();
388 B2 = ea2.getB();
389 alpha2 = ea2.getAlpha();
390 beta2 = ea2.getBeta();
391 kappa2 = ea2.getKappa();
392 lambda2 = ea2.getLambda();
393
394 RealType rmax = 0.0;
395 rmax = max(rmax, data1.rcut);
396 rmax = max(rmax, data2.rcut);
397 rmax = max(rmax, eamRcut_);
398
399 // use the smallest dr (finest grid) to build our grid:
400
401 RealType dr = min(data1.rho->getSpacing(), data2.rho->getSpacing());
402
403 int nr = int(rmax / dr + 1);
404
405 for (int i = 0; i < nr; i++)
406 rvals.push_back(RealType(i * dr));
407
408 vector<RealType> phivals;
409 RealType r;
410
411 for (unsigned int i = 0; i < rvals.size(); i++) {
412 r = rvals[i];
413 rha = 0.0;
414 rhb = 0.0;
415 pha = 0.0;
416 phb = 0.0;
417 phab = 0.0;
418
419 // rcut values are derived parameters if no splines are pre-set:
420
421 if (r < data1.rcut) {
422 rha = data1.rho->getValueAt(r);
423
424 // pha = ZhouPhi(r, re1, A1, B1, alpha1, beta1, kappa1, lambda1);
425 pha = Phi(r, re1, A1, B1, alpha1, beta1, kappa1, lambda1);
426 }
427 if (r < data2.rcut) {
428 rhb = data2.rho->getValueAt(r);
429 // phb = ZhouPhi(r, re2, A2, B2, alpha2, beta2, kappa2, lambda2);
430 phb = Phi(r, re2, A2, B2, alpha2, beta2, kappa2, lambda2);
431 }
432
433 if (r < data1.rcut) phab = phab + 0.5 * (rhb / rha) * pha;
434 if (r < data2.rcut) phab = phab + 0.5 * (rha / rhb) * phb;
435
436 phivals.push_back(phab);
437 }
438 cs->addPoints(rvals, phivals);
439 }
440
441 return cs;
442 }
443
444 void EAM::setCutoffRadius(RealType rCut) {
445 eamRcut_ = rCut;
446 haveCutoffRadius_ = true;
447 }
448
449 void EAM::initialize() {
450 // set up the mixing method:
451 ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
452 string EAMMixMeth = fopts.getEAMMixingMethod();
453 toUpper(EAMMixMeth);
454 if (EAMMixMeth == "JOHNSON")
455 mixMeth_ = eamJohnson;
456 else if (EAMMixMeth == "DREAM1")
457 mixMeth_ = eamDream1;
458 else if (EAMMixMeth == "DREAM2")
459 mixMeth_ = eamDream2;
460 else if (EAMMixMeth == "DAW")
461 mixMeth_ = eamDaw;
462 else
463 mixMeth_ = eamUnknownMix;
464
465 // find all of the EAM atom Types:
466 EAMtypes.clear();
467 EAMtids.clear();
468 EAMdata.clear();
469 MixingMap.clear();
470 nEAM_ = 0;
471
472 EAMtids.resize(forceField_->getNAtomType(), -1);
473
474 AtomTypeSet::iterator at;
475 for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
476 if ((*at)->isEAM()) nEAM_++;
477 }
478 EAMdata.resize(nEAM_);
479 MixingMap.resize(nEAM_);
480
481 for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
482 if ((*at)->isEAM()) addType(*at);
483 }
484
485 // find all of the explicit EAM interactions:
487 forceField_->getNonBondedInteractionTypes();
488 ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
489 ForceField::NonBondedInteractionTypeContainer::KeyType keys;
491
492 for (nbt = nbiTypes->beginType(j); nbt != NULL;
493 nbt = nbiTypes->nextType(j)) {
494 if (nbt->isEAMTable() || nbt->isEAMZhou() || nbt->isEAMOxides()) {
495 keys = nbiTypes->getKeys(j);
496 AtomType* at1 = forceField_->getAtomType(keys[0]);
497 if (at1 == NULL) {
498 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
499 "EAM::initialize could not find AtomType %s\n"
500 "\tfor %s - %s explicit interaction.\n",
501 keys[0].c_str(), keys[0].c_str(), keys[1].c_str());
502 painCave.severity = OPENMD_ERROR;
503 painCave.isFatal = 1;
504 simError();
505 }
506
507 AtomType* at2 = forceField_->getAtomType(keys[1]);
508 if (at2 == NULL) {
509 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
510 "EAM::initialize could not find AtomType %s\n"
511 "\tfor %s - %s explicit interaction.\n",
512 keys[1].c_str(), keys[0].c_str(), keys[1].c_str());
513 painCave.severity = OPENMD_ERROR;
514 painCave.isFatal = 1;
515 simError();
516 }
517
518 EAMInteractionType* eamit = dynamic_cast<EAMInteractionType*>(nbt);
519
520 if (eamit == NULL) {
521 snprintf(
522 painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
523 "EAM::initialize could not convert NonBondedInteractionType\n"
524 "\tto EAMInteractionType for %s - %s interaction.\n",
525 at1->getName().c_str(), at2->getName().c_str());
526 painCave.severity = OPENMD_ERROR;
527 painCave.isFatal = 1;
528 simError();
529 }
530 if (nbt->isEAMTable()) {
531 vector<RealType> phiAB = eamit->getPhi();
532 RealType dr = eamit->getDr();
533 int nr = eamit->getNr();
534
535 addExplicitInteraction(at1, at2, dr, nr, phiAB);
536
537 } else if (nbt->isEAMZhou()) {
538 RealType re = eamit->getRe();
539 RealType alpha = eamit->getAlpha();
540 RealType beta = eamit->getBeta();
541 RealType A = eamit->getA();
542 RealType B = eamit->getB();
543 RealType kappa = eamit->getKappa();
544 RealType lambda = eamit->getLambda();
545
546 addExplicitInteraction(at1, at2, re, alpha, beta, A, B, kappa,
547 lambda);
548
549 } else if (nbt->isEAMOxides()) {
550 RealType re = eamit->getRe();
551 RealType alpha = eamit->getAlpha();
552 RealType A = eamit->getA();
553 RealType Ci = eamit->getCi();
554 RealType Cj = eamit->getCj();
555 addExplicitInteraction(at1, at2, re, alpha, A, Ci, Cj);
556 }
557 }
558 }
559 initialized_ = true;
560 }
561
562 void EAM::addType(AtomType* atomType) {
563 EAMAdapter ea = EAMAdapter(atomType);
564 EAMAtomData eamAtomData;
565 EAMType et = ea.getEAMType();
566
567 switch (et) {
568 case eamFuncfl: {
569 eamAtomData.rho = ea.getRhoSpline();
570 eamAtomData.F = ea.getFSpline();
571 eamAtomData.Z = ea.getZSpline();
572 eamAtomData.rcut = ea.getRcut();
573 break;
574 }
575 case eamZhou2001:
576 case eamZhou2004:
577 case eamZhou2005:
578 case eamZhouRose: {
579 RealType re = ea.getRe();
580 RealType fe = ea.get_fe();
581 RealType rhoe = ea.getRhoe();
582 RealType A = ea.getA();
583 RealType B = ea.getB();
584 RealType alpha = ea.getAlpha();
585 RealType beta = ea.getBeta();
586 RealType kappa = ea.getKappa();
587 RealType lambda = ea.getLambda();
588 std::vector<RealType> Fn = ea.getFn();
589 std::vector<RealType> F = ea.getF();
590 RealType Fe = ea.getFe();
591 RealType eta = ea.getEta();
592 RealType F0 = ea.getF0();
593 // RealType latticeConstant = ea.getLatticeConstant();
594
595 int Nr = 2000;
596 // eamAtomData.rcut = latticeConstant * sqrt(10.0) / 2.0;
597 // eamAtomData.rcut = re * (pow(10.0, 0.3) + lambda);
598 eamAtomData.rcut = std::min(12.0, std::max(9.0, 4.0 * re));
599 // eamAtomData.rcut = re + (20.0 * re / beta);
600 RealType dr = eamAtomData.rcut / (RealType)(Nr - 1);
601 RealType r;
602
603 int Nrho = 3000;
604 RealType rhomax =
605 max(900.0, max(2.0 * rhoe, Rho(0.0, re, fe, beta, lambda)));
606 RealType drho = rhomax / (RealType)(Nrho - 1);
607 RealType rho;
608 RealType phiCC, phiCV;
609
610 std::vector<RealType> rvals;
611 std::vector<RealType> zvals;
612 std::vector<RealType> ccvals;
613 std::vector<RealType> cvvals;
614
615 std::vector<RealType> rhovals;
616 std::vector<RealType> funcvals;
617
618 for (int i = 0; i < Nr; i++) {
619 r = RealType(i) * dr;
620 rvals.push_back(r);
621 rhovals.push_back(Rho(r, re, fe, beta, lambda));
622 phiCC = PhiCoreCore(r, re, A, alpha, kappa);
623 phiCV = PhiCoreValence(r, re, B, beta, lambda);
624 ccvals.push_back(phiCC);
625 cvvals.push_back(phiCV);
626 }
627 eamAtomData.rho = std::make_shared<CubicSpline>();
628 eamAtomData.rho->addPoints(rvals, rhovals);
629 eamAtomData.phiCC = std::make_shared<CubicSpline>();
630 eamAtomData.phiCC->addPoints(rvals, ccvals);
631 eamAtomData.phiCV = std::make_shared<CubicSpline>();
632 eamAtomData.phiCV->addPoints(rvals, cvvals);
633
634 rhovals.clear();
635 if (et == eamZhou2001) {
636 for (int i = 0; i < Nrho; i++) {
637 rho = RealType(i) * drho;
638 rhovals.push_back(rho);
639 funcvals.push_back(Zhou2001Functional(rho, rhoe, Fn, F, Fe, eta));
640 }
641 } else if (et == eamZhou2004) {
642 RealType rhos = ea.getRhos();
643 RealType rhol = ea.getRhol();
644 RealType rhoh = ea.getRhoh();
645 for (int i = 0; i < Nrho; i++) {
646 rho = RealType(i) * drho;
647 rhovals.push_back(rho);
648 funcvals.push_back(
649 Zhou2004Functional(rho, rhoe, rhos, Fn, F, Fe, eta, rhol, rhoh));
650 }
651 } else if (et == eamZhou2005) {
652 RealType rhos = ea.getRhos();
653 RealType F3plus = ea.getF3plus();
654 RealType F3minus = ea.getF3minus();
655 for (int i = 0; i < Nrho; i++) {
656 rho = RealType(i) * drho;
657 rhovals.push_back(rho);
658 funcvals.push_back(Zhou2005Functional(rho, rhoe, rhos, Fn, F, F3plus,
659 F3minus, Fe, eta));
660 }
661 } else if (et == eamZhouRose) {
662 for (int i = 0; i < Nrho; i++) {
663 rho = RealType(i) * drho;
664 rhovals.push_back(rho);
665 funcvals.push_back(RoseFunctional(rho, rhoe, F0));
666 }
667 }
668
669 eamAtomData.F = std::make_shared<CubicSpline>();
670 eamAtomData.F->addPoints(rhovals, funcvals);
671 break;
672 }
673 case eamOxygenFuncfl: {
674 RealType re = ea.getRe();
675 RealType fe = ea.get_fe();
676
677 RealType A = ea.getA();
678 RealType B = ea.getB();
679 RealType alpha = ea.getAlpha();
680 RealType beta = ea.getBeta();
681 RealType kappa = ea.getKappa();
682 RealType lambda = ea.getLambda();
683
684 // RealType latticeConstant = ea.getLatticeConstant();
685
686 int Nr = 2000;
687 // eamAtomData.rcut = latticeConstant * sqrt(10.0) / 2.0;
688 // eamAtomData.rcut = re * (pow(10.0, 0.3) + lambda);
689 eamAtomData.rcut = std::min(12.0, std::max(9.0, 4.0 * re));
690 // eamAtomData.rcut = re + (20.0 * re / beta);
691 RealType dr = eamAtomData.rcut / (RealType)(Nr - 1);
692 RealType r;
693 RealType phiCC, phiCV;
694
695 std::vector<RealType> rvals;
696 std::vector<RealType> zvals;
697 std::vector<RealType> ccvals;
698 std::vector<RealType> cvvals;
699
700 std::vector<RealType> rhovals;
701
702 for (int i = 0; i < Nr; i++) {
703 r = RealType(i) * dr;
704 rvals.push_back(r);
705 rhovals.push_back(Rho(r, re, fe, beta, lambda));
706 phiCC = PhiCoreCore(r, re, A, alpha, kappa);
707 phiCV = PhiCoreValence(r, re, B, beta, lambda);
708 ccvals.push_back(phiCC);
709 cvvals.push_back(phiCV);
710 }
711 eamAtomData.rho = std::make_shared<CubicSpline>();
712 eamAtomData.rho->addPoints(rvals, rhovals);
713 eamAtomData.phiCC = std::make_shared<CubicSpline>();
714 eamAtomData.phiCC->addPoints(rvals, ccvals);
715 eamAtomData.phiCV = std::make_shared<CubicSpline>();
716 eamAtomData.phiCV->addPoints(rvals, cvvals);
717
718 eamAtomData.F = ea.getFSpline();
719 break;
720 }
721 case eamZhou2005Oxygen: {
722 RealType re = ea.getRe();
723 RealType fe = ea.get_fe();
724 RealType A = ea.getA();
725 RealType B = ea.getB();
726 RealType alpha = ea.getAlpha();
727 RealType beta = ea.getBeta();
728 RealType kappa = ea.getKappa();
729 RealType lambda = ea.getLambda();
730 RealType gamma = ea.getGamma();
731 RealType nu = ea.getNu();
732 std::vector<RealType> OrhoLimits = ea.getOrhoLimits();
733 std::vector<RealType> OrhoE = ea.getOrhoE();
734 std::vector<std::vector<RealType>> OF = ea.getOF();
735
736 int Nr = 2000;
737 // eamAtomData.rcut = 6.0;
738 // eamAtomData.rcut = re * (pow(10.0, 3.0/10.0) + nu );
739
740 eamAtomData.rcut = std::min(12.0, std::max(9.0, 4.0 * re));
741 // eamAtomData.rcut = re + (20.0 * re / beta);
742 RealType dr = eamAtomData.rcut / (RealType)(Nr - 1);
743 RealType r;
744
745 int Nrho = 3000;
746 // RealType rhomax = max(800.0, ZhouRho(0.0, re, fe, gamma, nu));
747 RealType rhomax = max(800.0, Rho(0.0, re, fe, gamma, nu));
748 RealType drho = rhomax / (RealType)(Nrho - 1);
749 RealType rho, phiCC, phiCV;
750
751 std::vector<RealType> rvals;
752 std::vector<RealType> zvals;
753 std::vector<RealType> ccvals;
754 std::vector<RealType> cvvals;
755
756 std::vector<RealType> rhovals;
757 std::vector<RealType> funcvals;
758
759 for (int i = 0; i < Nr; i++) {
760 r = RealType(i) * dr;
761 rvals.push_back(r);
762 rhovals.push_back(Rho(r, re, fe, gamma, nu));
763 phiCC = PhiCoreCore(r, re, A, alpha, kappa);
764 phiCV = PhiCoreValence(r, re, B, beta, lambda);
765 ccvals.push_back(phiCC);
766 cvvals.push_back(phiCV);
767 }
768 eamAtomData.rho = std::make_shared<CubicSpline>();
769 eamAtomData.rho->addPoints(rvals, rhovals);
770 eamAtomData.phiCC = std::make_shared<CubicSpline>();
771 eamAtomData.phiCC->addPoints(rvals, ccvals);
772 eamAtomData.phiCV = std::make_shared<CubicSpline>();
773 eamAtomData.phiCV->addPoints(rvals, cvvals);
774
775 rhovals.clear();
776 for (int i = 0; i < Nrho; i++) {
777 rho = RealType(i) * drho;
778 rhovals.push_back(rho);
779 funcvals.push_back(
780 Zhou2005OxygenFunctional(rho, OrhoLimits, OrhoE, OF));
781 }
782
783 eamAtomData.F = std::make_shared<CubicSpline>();
784 eamAtomData.F->addPoints(rhovals, funcvals);
785 break;
786 }
787
788 case eamUnknown: {
789 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
790 "EAM::addType found an unknown EAM type for atomType %s\n",
791 atomType->getName().c_str());
792 painCave.severity = OPENMD_ERROR;
793 painCave.isFatal = 1;
794 simError();
795 break;
796 }
797 }
798
800 if (fqa.isFluctuatingCharge()) {
801 if (fqa.isMetallic()) {
802 eamAtomData.isFluctuatingCharge = true;
803 eamAtomData.nValence = fqa.getNValence();
804 eamAtomData.nMobile = fqa.getNMobile();
805 } else {
806 eamAtomData.isFluctuatingCharge = false;
807 }
808 } else {
809 eamAtomData.isFluctuatingCharge = false;
810 }
811
812 // add it to the map:
813 int atid = atomType->getIdent();
814 int eamtid = EAMtypes.size();
815
816 pair<set<int>::iterator, bool> ret;
817 ret = EAMtypes.insert(atid);
818 if (ret.second == false) {
819 snprintf(painCave.errMsg, MAX_SIM_ERROR_MSG_LENGTH,
820 "EAM already had a previous entry with ident %d\n", atid);
821 painCave.severity = OPENMD_INFO;
822 painCave.isFatal = 0;
823 simError();
824 }
825
826 EAMtids[atid] = eamtid;
827 EAMdata[eamtid] = eamAtomData;
828 MixingMap[eamtid].resize(nEAM_);
829
830 // Now, iterate over all known types and add to the mixing map:
831
832 std::set<int>::iterator it;
833 for (it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
834 int eamtid2 = EAMtids[(*it)];
835 AtomType* atype2 = forceField_->getAtomType((*it));
836
837 EAMInteractionData mixer;
838 mixer.phi = getPhi(atomType, atype2);
839 mixer.rcut = mixer.phi->getLimits().second;
840 mixer.explicitlySet = false;
841
842 MixingMap[eamtid2].resize(nEAM_);
843
844 MixingMap[eamtid][eamtid2] = mixer;
845 if (eamtid2 != eamtid) { MixingMap[eamtid2][eamtid] = mixer; }
846 }
847 return;
848 }
849
850 void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
851 RealType dr, int nr,
852 vector<RealType> phiVals) {
853 EAMInteractionData mixer;
854 CubicSplinePtr cs {std::make_shared<CubicSpline>()};
855 vector<RealType> rVals;
856
857 for (int i = 0; i < nr; i++)
858 rVals.push_back(i * dr);
859
860 cs->addPoints(rVals, phiVals);
861 mixer.phi = cs;
862 mixer.rcut = mixer.phi->getLimits().second;
863 mixer.explicitlySet = true;
864
865 int atid1 = atype1->getIdent();
866 int atid2 = atype2->getIdent();
867
868 int eamtid1, eamtid2;
869
870 pair<set<int>::iterator, bool> ret;
871 ret = EAMtypes.insert(atid1);
872 if (ret.second == false) {
873 // already had this type in the EAMMap, just get the eamtid:
874 eamtid1 = EAMtids[atid1];
875 } else {
876 // didn't already have it, so make a new one and assign it:
877 eamtid1 = nEAM_;
878 EAMtids[atid1] = nEAM_;
879 nEAM_++;
880 }
881
882 ret = EAMtypes.insert(atid2);
883 if (ret.second == false) {
884 // already had this type in the EAMMap, just get the eamtid:
885 eamtid2 = EAMtids[atid2];
886 } else {
887 // didn't already have it, so make a new one and assign it:
888 eamtid2 = nEAM_;
889 EAMtids[atid2] = nEAM_;
890 nEAM_++;
891 }
892
893 MixingMap.resize(nEAM_);
894 MixingMap[eamtid1].resize(nEAM_);
895 MixingMap[eamtid1][eamtid2] = mixer;
896
897 if (eamtid2 != eamtid1) {
898 MixingMap[eamtid2].resize(nEAM_);
899 MixingMap[eamtid2][eamtid1] = mixer;
900 }
901 return;
902 }
903
904 void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
905 RealType re, RealType alpha, RealType beta,
906 RealType A, RealType B, RealType kappa,
907 RealType lambda) {
908 CubicSplinePtr cs {std::make_shared<CubicSpline>()};
909
910 EAMInteractionData mixer;
911 std::vector<RealType> rVals;
912 std::vector<RealType> phiVals;
913 std::vector<RealType> jVals;
914
915 int Nr = 2000;
916 RealType r;
917
918 // default to FCC if we don't specify HCP or BCC:
919 // RealType rcut = sqrt(5.0) * re;
920 RealType rcut = std::min(12.0, std::max(9.0, 4.0 * re));
921 RealType dr = rcut / (RealType)(Nr - 1);
922
923 for (int i = 0; i < Nr; i++) {
924 r = RealType(i) * dr;
925 rVals.push_back(r);
926 phiVals.push_back(Phi(r, re, A, B, alpha, beta, kappa, lambda));
927 }
928 cs->addPoints(rVals, phiVals);
929 mixer.phi = cs;
930 mixer.rcut = mixer.phi->getLimits().second;
931 mixer.explicitlySet = true;
932
933 int atid1 = atype1->getIdent();
934 int atid2 = atype2->getIdent();
935
936 int eamtid1, eamtid2;
937
938 pair<set<int>::iterator, bool> ret;
939 ret = EAMtypes.insert(atid1);
940 if (ret.second == false) {
941 // already had this type in the EAMMap, just get the eamtid:
942 eamtid1 = EAMtids[atid1];
943 } else {
944 // didn't already have it, so make a new one and assign it:
945 eamtid1 = nEAM_;
946 EAMtids[atid1] = nEAM_;
947 nEAM_++;
948 }
949
950 ret = EAMtypes.insert(atid2);
951 if (ret.second == false) {
952 // already had this type in the EAMMap, just get the eamtid:
953 eamtid2 = EAMtids[atid2];
954 } else {
955 // didn't already have it, so make a new one and assign it:
956 eamtid2 = nEAM_;
957 EAMtids[atid2] = nEAM_;
958 nEAM_++;
959 }
960
961 MixingMap.resize(nEAM_);
962 MixingMap[eamtid1].resize(nEAM_);
963 MixingMap[eamtid1][eamtid2] = mixer;
964 if (eamtid2 != eamtid1) {
965 MixingMap[eamtid2].resize(nEAM_);
966 MixingMap[eamtid2][eamtid1] = mixer;
967 }
968 return;
969 }
970
971 void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
972 RealType re, RealType alpha, RealType A,
973 RealType Ci, RealType Cj) {
974 CubicSplinePtr cs {std::make_shared<CubicSpline>()};
975
976 EAMInteractionData mixer;
977 std::vector<RealType> rVals;
978 std::vector<RealType> phiVals;
979
980 int Nr = 2000;
981 RealType r;
982 RealType phiCC;
983
984 RealType rcut = std::min(12.0, std::max(9.0, 4.0 * re));
985 RealType dr = rcut / (RealType)(Nr - 1);
986
987 for (int i = 0; i < Nr; i++) {
988 r = RealType(i) * dr;
989 rVals.push_back(r);
990 // for Dream2 and EAMOxide explicit potential we dont care about kappa
991 // but it is needed while calling PhiCoreCore. So, lets define it as 0
992 // its value however is not used in calculation
993 RealType kappa = 0;
994 phiCC = PhiCoreCore(r, re, A, alpha, kappa);
995 phiVals.push_back(phiCC);
996 }
997
998 cs->addPoints(rVals, phiVals);
999 // this is the repulsive piece in explicit EAM Oxide potential
1000 // this is phiCC for the explicit EAM oxide potential for convinence, it
1001 // is represented as phi
1002 mixer.phi = cs;
1003 mixer.rcut = mixer.phi->getLimits().second;
1004
1005 mixer.Ci = Ci;
1006 mixer.Cj = Cj;
1007
1008 mixer.explicitlySet = true;
1009
1010 int atid1 = atype1->getIdent();
1011 int atid2 = atype2->getIdent();
1012
1013 int eamtid1, eamtid2;
1014
1015 pair<set<int>::iterator, bool> ret;
1016 ret = EAMtypes.insert(atid1);
1017 if (ret.second == false) {
1018 // already had this type in the EAMMap, just get the eamtid:
1019 eamtid1 = EAMtids[atid1];
1020 } else {
1021 // didn't already have it, so make a new one and assign it:
1022 eamtid1 = nEAM_;
1023 EAMtids[atid1] = nEAM_;
1024 nEAM_++;
1025 }
1026
1027 ret = EAMtypes.insert(atid2);
1028 if (ret.second == false) {
1029 // already had this type in the EAMMap, just get the eamtid:
1030 eamtid2 = EAMtids[atid2];
1031 } else {
1032 // didn't already have it, so make a new one and assign it:
1033 eamtid2 = nEAM_;
1034 EAMtids[atid2] = nEAM_;
1035 nEAM_++;
1036 }
1037
1038 MixingMap.resize(nEAM_);
1039 MixingMap[eamtid1].resize(nEAM_);
1040 MixingMap[eamtid1][eamtid2] = mixer;
1041 if (eamtid2 != eamtid1) {
1042 MixingMap[eamtid2].resize(nEAM_);
1043 MixingMap[eamtid2][eamtid1] = mixer;
1044 }
1045 return;
1046 }
1047
1048 void EAM::calcDensity(InteractionData& idat) {
1049 if (!initialized_) initialize();
1050
1051 EAMAtomData& data1 = EAMdata[EAMtids[idat.atid1]];
1052 EAMAtomData& data2 = EAMdata[EAMtids[idat.atid2]];
1053 RealType s;
1054
1055 if (haveCutoffRadius_)
1056 if (idat.rij > eamRcut_) return;
1057
1058 if (idat.rij < data1.rcut) {
1059 s = 1.0;
1060 if (data1.isFluctuatingCharge) {
1061 if (mixMeth_ == eamDream2)
1062 s = gFunc(idat.flucQ1, data1.nValence, data1.nMobile);
1063 else
1064 s = (data1.nValence - idat.flucQ1) / (data1.nValence);
1065 }
1066 idat.rho2 += s * data1.rho->getValueAt(idat.rij);
1067 }
1068
1069 if (idat.rij < data2.rcut) {
1070 s = 1.0;
1071 if (data2.isFluctuatingCharge) {
1072 if (mixMeth_ == eamDream2)
1073 s = gFunc(idat.flucQ2, data2.nValence, data2.nMobile);
1074 else
1075 s = (data2.nValence - idat.flucQ2) / (data2.nValence);
1076 }
1077 idat.rho1 += s * data2.rho->getValueAt(idat.rij);
1078 }
1079
1080 return;
1081 }
1082
1083 void EAM::calcFunctional(SelfData& sdat) {
1084 if (!initialized_) initialize();
1085 EAMAtomData& data1 = EAMdata[EAMtids[sdat.atid]];
1086
1087 data1.F->getValueAndDerivativeAt(sdat.rho, sdat.frho, sdat.dfrhodrho);
1088
1090
1091 if (sdat.isSelected) sdat.selePot[METALLIC_EMBEDDING_FAMILY] += sdat.frho;
1092
1093 if (sdat.doParticlePot) sdat.particlePot += sdat.frho;
1094
1095 return;
1096 }
1097
1098 void EAM::calcForce(InteractionData& idat) {
1099 if (!initialized_) initialize();
1100
1101 if (haveCutoffRadius_)
1102 if (idat.rij > eamRcut_) return;
1103
1104 int eamtid1 = EAMtids[idat.atid1];
1105 int eamtid2 = EAMtids[idat.atid2];
1106 EAMAtomData& data1 = EAMdata[eamtid1];
1107 EAMAtomData& data2 = EAMdata[eamtid2];
1108
1109 // get type-specific cutoff radii
1110
1111 RealType rci = data1.rcut;
1112 RealType rcj = data2.rcut;
1113 RealType rcij = MixingMap[eamtid1][eamtid2].rcut;
1114
1115 RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
1116 RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
1117
1118 RealType phab(0.0), dvpdr(0.0);
1119 RealType drhoidr(0.0), drhojdr(0.0), dudr(0.0);
1120 // For DREAM, we need nValence (Va, Vb), nMobile (Ma, Mb), for
1121 // each type, and density scaling variables (si, sj) for each atom:
1122 RealType Va(0.0), Vb(0.0);
1123 RealType Ma(0.0), Mb(0.0);
1124 RealType si(1.0), sj(1.0);
1125 RealType sip(0.0), sjp(0.0);
1126 RealType Ci(1.0), Cj(1.0);
1127
1128 rhat = idat.d / idat.rij;
1129 if (idat.rij < rci && idat.rij < rcij) {
1130 data1.rho->getValueAndDerivativeAt(idat.rij, rha, drha);
1131 CubicSplinePtr phi = MixingMap[eamtid1][eamtid1].phi;
1132 phi->getValueAndDerivativeAt(idat.rij, pha, dpha);
1133 }
1134
1135 if (idat.rij < rcj && idat.rij < rcij) {
1136 data2.rho->getValueAndDerivativeAt(idat.rij, rhb, drhb);
1137 CubicSplinePtr phi = MixingMap[eamtid2][eamtid2].phi;
1138 phi->getValueAndDerivativeAt(idat.rij, phb, dphb);
1139 }
1140
1141 bool hasFlucQ = data1.isFluctuatingCharge || data2.isFluctuatingCharge;
1142 bool isExplicit = MixingMap[eamtid1][eamtid2].explicitlySet;
1143 if (hasFlucQ) {
1144 if (data1.isFluctuatingCharge) {
1145 Va = data1.nValence;
1146 Ma = data1.nMobile;
1147 if (mixMeth_ == eamDream2)
1148 si = gFunc(idat.flucQ1, Va, Ma);
1149 else
1150 si = (Va - idat.flucQ1) / Va;
1151 }
1152 if (data2.isFluctuatingCharge) {
1153 Vb = data2.nValence;
1154 Mb = data2.nMobile;
1155 if (mixMeth_ == eamDream2)
1156 sj = gFunc(idat.flucQ2, Vb, Mb);
1157 else
1158 sj = (Vb - idat.flucQ2) / Vb;
1159 }
1160
1161 if (mixMeth_ == eamJohnson || mixMeth_ == eamDream1) {
1162 if (idat.rij < rci && idat.rij < rcij) {
1163 phab = phab + 0.5 * (sj / si) * (rhb / rha) * pha;
1164 dvpdr = dvpdr + 0.5 * (sj / si) *
1165 ((rhb / rha) * dpha +
1166 pha * ((drhb / rha) - (rhb * drha / rha / rha)));
1167
1168 if (data1.isFluctuatingCharge) {
1169 idat.dVdFQ1 += 0.5 * (rhb * sj * pha) / (rha * Ma * si * si);
1170 }
1171 if (data2.isFluctuatingCharge) {
1172 idat.dVdFQ2 -= 0.5 * (rhb * pha) / (rha * Mb * si);
1173 }
1174 }
1175
1176 if (idat.rij < rcj && idat.rij < rcij) {
1177 phab = phab + 0.5 * (si / sj) * (rha / rhb) * phb;
1178 dvpdr = dvpdr + 0.5 * (si / sj) *
1179 ((rha / rhb) * dphb +
1180 phb * ((drha / rhb) - (rha * drhb / rhb / rhb)));
1181
1182 if (data1.isFluctuatingCharge) {
1183 idat.dVdFQ1 -= 0.5 * (rha * phb) / (rhb * Ma * sj);
1184 }
1185 if (data2.isFluctuatingCharge) {
1186 idat.dVdFQ2 += 0.5 * (rha * si * phb) / (rhb * Mb * sj * sj);
1187 }
1188 }
1189 } else {
1190 if (isExplicit) {
1191 // phi is total potential for EAMTable and EAMZhou but CC interaction
1192 // for EAMOxide potential
1193 bool insideCutOff = (idat.rij < rcij);
1194 if (insideCutOff) {
1195 CubicSplinePtr phi = MixingMap[eamtid1][eamtid2].phi;
1196 phi->getValueAndDerivativeAt(idat.rij, phab, dvpdr);
1197 }
1198
1199 } else {
1200 // Core-Core part first - no fluctuating charge, just Johnson mixing:
1201
1202 if (idat.rij < rci && idat.rij < rcij) {
1203 CubicSplinePtr phiACC = data1.phiCC;
1204 phiACC->getValueAndDerivativeAt(idat.rij, pha, dpha);
1205 phab += 0.5 * (rhb / rha) * pha;
1206 dvpdr += 0.5 * ((rhb / rha) * dpha +
1207 pha * ((drhb / rha) - (rhb * drha / rha / rha)));
1208 }
1209 if (idat.rij < rcj && idat.rij < rcij) {
1210 CubicSplinePtr phiBCC = data2.phiCC;
1211 phiBCC->getValueAndDerivativeAt(idat.rij, phb, dphb);
1212 phab += 0.5 * (rha / rhb) * phb;
1213 dvpdr += 0.5 * ((rha / rhb) * dphb +
1214 phb * ((drha / rhb) - (rha * drhb / rhb / rhb)));
1215 }
1216 }
1217
1218 if (isExplicit) {
1219 // Ci, Cj are 0 for EAMTable and EAMZhou. So, there is no CV
1220 // interaction
1221 Ci = MixingMap[eamtid1][eamtid2].Ci;
1222 Cj = MixingMap[eamtid1][eamtid2].Cj;
1223 }
1224 // Core-Valence next, this does include fluctuating charges:
1225
1226 if (data1.isFluctuatingCharge) { sip = gPrime(idat.flucQ1, Va, Ma); }
1227 if (data2.isFluctuatingCharge) { sjp = gPrime(idat.flucQ2, Vb, Mb); }
1228
1229 if (idat.rij < rci && idat.rij < rcij) {
1230 CubicSplinePtr phiACV = data1.phiCV;
1231 phiACV->getValueAndDerivativeAt(idat.rij, pha, dpha);
1232
1233 phab += 0.5 * sj * Ci * (rhb / rha) * pha;
1234 dvpdr += 0.5 * sj * Ci *
1235 ((rhb / rha) * dpha +
1236 pha * ((drhb / rha) - (rhb * drha / rha / rha)));
1237
1238 if (data2.isFluctuatingCharge) {
1239 idat.dVdFQ2 += 0.5 * sjp * Ci * (rhb / rha) * pha;
1240 }
1241 }
1242 if (idat.rij < rcj && idat.rij < rcij) {
1243 CubicSplinePtr phiBCV = data2.phiCV;
1244 phiBCV->getValueAndDerivativeAt(idat.rij, phb, dphb);
1245
1246 phab += 0.5 * si * Cj * (rha / rhb) * phb;
1247 dvpdr += 0.5 * si * Cj *
1248 ((rha / rhb) * dphb +
1249 phb * ((drha / rhb) - (rha * drhb / rhb / rhb)));
1250
1251 if (data1.isFluctuatingCharge) {
1252 idat.dVdFQ1 += 0.5 * sip * Cj * (rha / rhb) * phb;
1253 }
1254 }
1255 }
1256 } else {
1257 if (idat.rij < rcij) {
1258 CubicSplinePtr phi = MixingMap[eamtid1][eamtid2].phi;
1259 phi->getValueAndDerivativeAt(idat.rij, phab, dvpdr);
1260 }
1261 }
1262 drhoidr = si * drha;
1263 drhojdr = sj * drhb;
1264 dudr = drhojdr * idat.dfrho1 + drhoidr * idat.dfrho2 + dvpdr;
1265
1266 idat.f1 += rhat * dudr;
1267
1268 if (idat.doParticlePot) {
1269 // particlePot is the difference between the full potential and
1270 // the full potential without the presence of a particular
1271 // particle (atom1).
1272 //
1273 // This reduces the density at other particle locations, so we
1274 // need to recompute the density at atom2 assuming atom1 didn't
1275 // contribute. This then requires recomputing the density
1276 // functional for atom2 as well.
1277
1278 idat.particlePot1 += data2.F->getValueAt(idat.rho2 - rha) - idat.frho2;
1279
1280 idat.particlePot2 += data1.F->getValueAt(idat.rho1 - rhb) - idat.frho1;
1281 }
1282
1283 idat.pot[METALLIC_PAIR_FAMILY] += phab;
1284 if (idat.isSelected) idat.selePot[METALLIC_PAIR_FAMILY] += phab;
1285
1286 idat.vpair += phab;
1287
1288 // When densities are fluctuating, the functional depends on the
1289 // fluctuating densities from other sites:
1290 if (data1.isFluctuatingCharge) {
1291 if (mixMeth_ == eamDream2)
1292 idat.dVdFQ1 += idat.dfrho2 * rha * sip;
1293 else
1294 idat.dVdFQ1 -= idat.dfrho2 * rha / Va;
1295 }
1296 if (data2.isFluctuatingCharge) {
1297 if (mixMeth_ == eamDream2)
1298 idat.dVdFQ2 += idat.dfrho1 * rhb * sjp;
1299 else
1300 idat.dVdFQ2 -= idat.dfrho1 * rhb / Vb;
1301 }
1302
1303 return;
1304 }
1305
1306 RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
1307 if (!initialized_) initialize();
1308
1309 RealType cut = 0.0;
1310
1311 int atid1 = atypes.first->getIdent();
1312 int atid2 = atypes.second->getIdent();
1313 int eamtid1 = EAMtids[atid1];
1314 int eamtid2 = EAMtids[atid2];
1315
1316 if (eamtid1 != -1) {
1317 EAMAtomData data1 = EAMdata[eamtid1];
1318 cut = data1.rcut;
1319 }
1320
1321 if (eamtid2 != -1) {
1322 EAMAtomData data2 = EAMdata[eamtid2];
1323 if (data2.rcut > cut) cut = data2.rcut;
1324 }
1325
1326 return cut;
1327 }
1328} // namespace OpenMD
AtomType is what OpenMD looks to for unchanging data about an atom.
Definition AtomType.hpp:66
EAMInteractionType is one of the basic metallic interactions for representing the bonding in metallic...
AtomType * getAtomType(const std::string &at)
getAtomType by string
NonBondedInteractionType class is responsible for keeping track of static (unchanging) parameters for...
This basic Periodic Table class was originally taken from the data.cpp file in OpenBabel.
@ METALLIC_EMBEDDING_FAMILY
Transition metal interactions involving electron density.
@ METALLIC_PAIR_FAMILY
Transition metal interactions involving pair potentials.
The InteractionData struct.
Vector3d d
interatomic vector (already wrapped into box)
RealType dVdFQ1
fluctuating charge force on atom1
RealType dVdFQ2
fluctuating charge force on atom2
RealType particlePot2
particle potential for atom2
int atid1
atomType ident for atom 1
potVec selePot
potential energies of the selected sites
RealType flucQ2
fluctuating charge on atom2
RealType rho2
total electron density at second atom
bool isSelected
one of the particles has been selected for selection potential
RealType flucQ1
fluctuating charge on atom1
RealType vpair
pair potential
int atid2
atomType ident for atom 2
RealType dfrho2
derivative of functional for atom 2
RealType frho2
density functional at second atom
bool doParticlePot
should we bother with the particle pot?
RealType dfrho1
derivative of functional for atom 1
RealType frho1
density functional at first atom
Vector3d f1
force between the two atoms
RealType rij
interatomic separation
RealType particlePot1
particle potential for atom1
RealType rho1
total electron density at first atom
The SelfData struct.
potVec selePot
potential energy of the selected site
RealType frho
value of density functional for atom
RealType dfrhodrho
derivative of density functional for atom
potVec selfPot
total potential (including embedding energy)
RealType particlePot
contribution to potential from this particle
RealType rho
electron density
bool isSelected
this site has been selected for selection potential
bool doParticlePot
should we bother with the particle pot?
int atid
atomType ident for the atom