ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/nonbonded/EAM.cpp
Revision: 1582
Committed: Tue Jun 14 20:41:44 2011 UTC (13 years, 10 months ago) by gezelter
File size: 15635 byte(s)
Log Message:
Well, the potential energy values are still garbage, but the LJ sample
runs.

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Redistributions of source code must retain the above copyright
10 * notice, this list of conditions and the following disclaimer.
11 *
12 * 2. Redistributions in binary form must reproduce the above copyright
13 * notice, this list of conditions and the following disclaimer in the
14 * documentation and/or other materials provided with the
15 * distribution.
16 *
17 * This software is provided "AS IS," without a warranty of any
18 * kind. All express or implied conditions, representations and
19 * warranties, including any implied warranty of merchantability,
20 * fitness for a particular purpose or non-infringement, are hereby
21 * excluded. The University of Notre Dame and its licensors shall not
22 * be liable for any damages suffered by licensee as a result of
23 * using, modifying or distributing the software or its
24 * derivatives. In no event will the University of Notre Dame or its
25 * licensors be liable for any lost revenue, profit or data, or for
26 * direct, indirect, special, consequential, incidental or punitive
27 * damages, however caused and regardless of the theory of liability,
28 * arising out of the use of or inability to use software, even if the
29 * University of Notre Dame has been advised of the possibility of
30 * such damages.
31 *
32 * SUPPORT OPEN SCIENCE! If you use OpenMD or its source code in your
33 * research, please cite the appropriate papers when you publish your
34 * work. Good starting points are:
35 *
36 * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37 * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38 * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).
39 * [4] Vardeman & Gezelter, in progress (2009).
40 */
41
42 #include <stdio.h>
43 #include <string.h>
44
45 #include <cmath>
46 #include "nonbonded/EAM.hpp"
47 #include "utils/simError.h"
48 #include "types/NonBondedInteractionType.hpp"
49
50
51 namespace OpenMD {
52
53 EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
54 mixMeth_(eamJohnson), eamRcut_(0.0) {}
55
56 EAMParam EAM::getEAMParam(AtomType* atomType) {
57
58 // Do sanity checking on the AtomType we were passed before
59 // building any data structures:
60 if (!atomType->isEAM()) {
61 sprintf( painCave.errMsg,
62 "EAM::getEAMParam was passed an atomType (%s) that does not\n"
63 "\tappear to be an embedded atom method (EAM) atom.\n",
64 atomType->getName().c_str());
65 painCave.severity = OPENMD_ERROR;
66 painCave.isFatal = 1;
67 simError();
68 }
69
70 GenericData* data = atomType->getPropertyByName("EAM");
71 if (data == NULL) {
72 sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
73 "\tparameters for atomType %s.\n",
74 atomType->getName().c_str());
75 painCave.severity = OPENMD_ERROR;
76 painCave.isFatal = 1;
77 simError();
78 }
79
80 EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
81 if (eamData == NULL) {
82 sprintf( painCave.errMsg,
83 "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
84 "\tatom type %s\n", atomType->getName().c_str());
85 painCave.severity = OPENMD_ERROR;
86 painCave.isFatal = 1;
87 simError();
88 }
89
90 return eamData->getData();
91 }
92
93 CubicSpline* EAM::getZ(AtomType* atomType) {
94 EAMParam eamParam = getEAMParam(atomType);
95 int nr = eamParam.nr;
96 RealType dr = eamParam.dr;
97 vector<RealType> rvals;
98
99 for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
100
101 CubicSpline* cs = new CubicSpline();
102 cs->addPoints(rvals, eamParam.Z);
103 return cs;
104 }
105
106 RealType EAM::getRcut(AtomType* atomType) {
107 EAMParam eamParam = getEAMParam(atomType);
108 return eamParam.rcut;
109 }
110
111 CubicSpline* EAM::getRho(AtomType* atomType) {
112 EAMParam eamParam = getEAMParam(atomType);
113 int nr = eamParam.nr;
114 RealType dr = eamParam.dr;
115 vector<RealType> rvals;
116
117 for (int i = 0; i < nr; i++) rvals.push_back(RealType(i) * dr);
118
119 CubicSpline* cs = new CubicSpline();
120 cs->addPoints(rvals, eamParam.rho);
121 return cs;
122 }
123
124 CubicSpline* EAM::getF(AtomType* atomType) {
125 EAMParam eamParam = getEAMParam(atomType);
126 int nrho = eamParam.nrho;
127 RealType drho = eamParam.drho;
128 vector<RealType> rhovals;
129 vector<RealType> scaledF;
130
131 for (int i = 0; i < nrho; i++) {
132 rhovals.push_back(RealType(i) * drho);
133 scaledF.push_back( eamParam.F[i] * 23.06054 );
134 }
135
136 CubicSpline* cs = new CubicSpline();
137 cs->addPoints(rhovals, scaledF);
138 return cs;
139 }
140
141 CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {
142 EAMParam eamParam1 = getEAMParam(atomType1);
143 EAMParam eamParam2 = getEAMParam(atomType2);
144 CubicSpline* z1 = getZ(atomType1);
145 CubicSpline* z2 = getZ(atomType2);
146
147 // make the r grid:
148
149
150 // we need phi out to the largest value we'll encounter in the radial space;
151
152 RealType rmax = 0.0;
153 rmax = max(rmax, eamParam1.rcut);
154 rmax = max(rmax, eamParam1.nr * eamParam1.dr);
155
156 rmax = max(rmax, eamParam2.rcut);
157 rmax = max(rmax, eamParam2.nr * eamParam2.dr);
158
159 // use the smallest dr (finest grid) to build our grid:
160
161 RealType dr = min(eamParam1.dr, eamParam2.dr);
162
163 int nr = int(rmax/dr + 0.5);
164
165 vector<RealType> rvals;
166 for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
167
168 // construct the pair potential:
169
170 vector<RealType> phivals;
171 RealType phi;
172 RealType r;
173 RealType zi, zj;
174
175 phivals.push_back(0.0);
176
177 for (int i = 1; i < rvals.size(); i++ ) {
178 r = rvals[i];
179
180 // only use z(r) if we're inside this atom's cutoff radius,
181 // otherwise, we'll use zero for the charge. This effectively
182 // means that our phi grid goes out beyond the cutoff of the
183 // pair potential
184
185 zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
186 zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
187
188 phi = 331.999296 * (zi * zj) / r;
189
190 phivals.push_back(phi);
191 }
192
193 CubicSpline* cs = new CubicSpline();
194 cs->addPoints(rvals, phivals);
195 return cs;
196 }
197
198 void EAM::initialize() {
199
200 // set up the mixing method:
201 ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
202 string EAMMixMeth = fopts.getEAMMixingMethod();
203 toUpper(EAMMixMeth);
204
205 if (EAMMixMeth == "JOHNSON")
206 mixMeth_ = eamJohnson;
207 else if (EAMMixMeth == "DAW")
208 mixMeth_ = eamDaw;
209 else
210 mixMeth_ = eamUnknown;
211
212 // find all of the EAM atom Types:
213 ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
214 ForceField::AtomTypeContainer::MapTypeIterator i;
215 AtomType* at;
216
217 for (at = atomTypes->beginType(i); at != NULL;
218 at = atomTypes->nextType(i)) {
219
220 if (at->isEAM())
221 addType(at);
222 }
223
224 // find all of the explicit EAM interactions (setfl):
225 ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
226 ForceField::NonBondedInteractionTypeContainer::MapTypeIterator j;
227 NonBondedInteractionType* nbt;
228
229 for (nbt = nbiTypes->beginType(j); nbt != NULL;
230 nbt = nbiTypes->nextType(j)) {
231
232 if (nbt->isEAM()) {
233
234 pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
235
236 GenericData* data = nbt->getPropertyByName("EAM");
237 if (data == NULL) {
238 sprintf( painCave.errMsg, "EAM::rebuildMixingMap could not find\n"
239 "\tEAM parameters for %s - %s interaction.\n",
240 atypes.first->getName().c_str(),
241 atypes.second->getName().c_str());
242 painCave.severity = OPENMD_ERROR;
243 painCave.isFatal = 1;
244 simError();
245 }
246
247 EAMMixingData* eamData = dynamic_cast<EAMMixingData*>(data);
248 if (eamData == NULL) {
249 sprintf( painCave.errMsg,
250 "EAM::rebuildMixingMap could not convert GenericData to\n"
251 "\tEAMMixingData for %s - %s interaction.\n",
252 atypes.first->getName().c_str(),
253 atypes.second->getName().c_str());
254 painCave.severity = OPENMD_ERROR;
255 painCave.isFatal = 1;
256 simError();
257 }
258
259 EAMMixingParam eamParam = eamData->getData();
260
261 vector<RealType> phiAB = eamParam.phi;
262 RealType dr = eamParam.dr;
263 int nr = eamParam.nr;
264
265 addExplicitInteraction(atypes.first, atypes.second, dr, nr, phiAB);
266 }
267 }
268 initialized_ = true;
269 }
270
271
272
273 void EAM::addType(AtomType* atomType){
274
275 EAMAtomData eamAtomData;
276
277 eamAtomData.rho = getRho(atomType);
278 eamAtomData.F = getF(atomType);
279 eamAtomData.Z = getZ(atomType);
280 eamAtomData.rcut = getRcut(atomType);
281
282 // add it to the map:
283 AtomTypeProperties atp = atomType->getATP();
284
285 pair<map<int,AtomType*>::iterator,bool> ret;
286 ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
287 if (ret.second == false) {
288 sprintf( painCave.errMsg,
289 "EAM already had a previous entry with ident %d\n",
290 atp.ident);
291 painCave.severity = OPENMD_INFO;
292 painCave.isFatal = 0;
293 simError();
294 }
295
296 EAMMap[atomType] = eamAtomData;
297
298 // Now, iterate over all known types and add to the mixing map:
299
300 map<AtomType*, EAMAtomData>::iterator it;
301 for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
302
303 AtomType* atype2 = (*it).first;
304
305 EAMInteractionData mixer;
306 mixer.phi = getPhi(atomType, atype2);
307 mixer.explicitlySet = false;
308
309 pair<AtomType*, AtomType*> key1, key2;
310 key1 = make_pair(atomType, atype2);
311 key2 = make_pair(atype2, atomType);
312
313 MixingMap[key1] = mixer;
314 if (key2 != key1) {
315 MixingMap[key2] = mixer;
316 }
317 }
318 return;
319 }
320
321 void EAM::addExplicitInteraction(AtomType* atype1, AtomType* atype2,
322 RealType dr, int nr,
323 vector<RealType> phiVals) {
324
325 // in case these weren't already in the map
326 addType(atype1);
327 addType(atype2);
328
329 EAMInteractionData mixer;
330 CubicSpline* cs = new CubicSpline();
331 vector<RealType> rVals;
332
333 for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
334
335 cs->addPoints(rVals, phiVals);
336 mixer.phi = cs;
337 mixer.explicitlySet = true;
338
339 pair<AtomType*, AtomType*> key1, key2;
340 key1 = make_pair(atype1, atype2);
341 key2 = make_pair(atype2, atype1);
342
343 MixingMap[key1] = mixer;
344 if (key2 != key1) {
345 MixingMap[key2] = mixer;
346 }
347 return;
348 }
349
350 void EAM::calcDensity(InteractionData &idat) {
351
352 if (!initialized_) initialize();
353
354 EAMAtomData data1 = EAMMap[idat.atypes.first];
355 EAMAtomData data2 = EAMMap[idat.atypes.second];
356
357 if ( *(idat.rij) < data1.rcut)
358 *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
359
360
361 if ( *(idat.rij) < data2.rcut)
362 *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
363
364 return;
365 }
366
367 void EAM::calcFunctional(SelfData &sdat) {
368
369 if (!initialized_) initialize();
370
371 EAMAtomData data1 = EAMMap[ sdat.atype ];
372
373 pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
374
375 *(sdat.frho) = result.first;
376 *(sdat.dfrhodrho) = result.second;
377
378 sdat.pot[METALLIC_FAMILY] += result.first;
379 *(sdat.particlePot) += result.first;
380
381 return;
382 }
383
384
385 void EAM::calcForce(InteractionData &idat) {
386
387 if (!initialized_) initialize();
388
389 pair<RealType, RealType> res;
390
391 if ( *(idat.rij) < eamRcut_) {
392
393 EAMAtomData data1 = EAMMap[idat.atypes.first];
394 EAMAtomData data2 = EAMMap[idat.atypes.second];
395
396 // get type-specific cutoff radii
397
398 RealType rci = data1.rcut;
399 RealType rcj = data2.rcut;
400
401 RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
402 RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
403 RealType phab(0.0), dvpdr(0.0);
404 RealType drhoidr, drhojdr, dudr;
405
406 if ( *(idat.rij) < rci) {
407 res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
408 rha = res.first;
409 drha = res.second;
410
411 res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
412 pha = res.first;
413 dpha = res.second;
414 }
415
416 if ( *(idat.rij) < rcj) {
417 res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
418 rhb = res.first;
419 drhb = res.second;
420
421 res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
422 phb = res.first;
423 dphb = res.second;
424 }
425
426 switch(mixMeth_) {
427 case eamJohnson:
428
429 if ( *(idat.rij) < rci) {
430 phab = phab + 0.5 * (rhb / rha) * pha;
431 dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
432 pha*((drhb/rha) - (rhb*drha/rha/rha)));
433 }
434
435
436
437 if ( *(idat.rij) < rcj) {
438 phab = phab + 0.5 * (rha / rhb) * phb;
439 dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
440 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
441 }
442
443 break;
444
445 case eamDaw:
446 res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
447 phab = res.first;
448 dvpdr = res.second;
449
450 break;
451 case eamUnknown:
452 default:
453
454 sprintf(painCave.errMsg,
455 "EAM::calcForce hit a mixing method it doesn't know about!\n"
456 );
457 painCave.severity = OPENMD_ERROR;
458 painCave.isFatal = 1;
459 simError();
460
461 }
462
463 drhoidr = drha;
464 drhojdr = drhb;
465
466 dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
467
468 *(idat.f1) = *(idat.d) * dudr / *(idat.rij);
469
470 // particlePot is the difference between the full potential and
471 // the full potential without the presence of a particular
472 // particle (atom1).
473 //
474 // This reduces the density at other particle locations, so we
475 // need to recompute the density at atom2 assuming atom1 didn't
476 // contribute. This then requires recomputing the density
477 // functional for atom2 as well.
478
479 *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
480 - *(idat.frho2);
481
482 *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
483 - *(idat.frho1);
484
485 (*(idat.pot))[METALLIC_FAMILY] += phab;
486
487 *(idat.vpair) += phab;
488 }
489
490 return;
491
492 }
493
494 RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
495 if (!initialized_) initialize();
496
497 RealType cut = 0.0;
498
499 map<AtomType*, EAMAtomData>::iterator it;
500
501 it = EAMMap.find(atypes.first);
502 if (it != EAMMap.end()) {
503 EAMAtomData data1 = (*it).second;
504 cut = data1.rcut;
505 }
506
507 it = EAMMap.find(atypes.second);
508 if (it != EAMMap.end()) {
509 EAMAtomData data2 = (*it).second;
510 if (data2.rcut > cut)
511 cut = data2.rcut;
512 }
513
514 return cut;
515 }
516 }
517

Properties

Name Value
svn:eol-style native