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