ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/EAM.cpp
(Generate patch)

Comparing:
branches/development/src/nonbonded/EAM.cpp (file contents), Revision 1481 by gezelter, Mon Jul 26 21:55:18 2010 UTC vs.
trunk/src/nonbonded/EAM.cpp (file contents), Revision 1895 by gezelter, Mon Jul 1 21:09:37 2013 UTC

# Line 35 | Line 35
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).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 50 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  bool EAM::initialized_ = false;
55 <  RealType EAM::eamRcut_ = 0.0;
55 <  EAMMixingMethod EAM::mixMeth_ = eamJohnson;
56 <  ForceField* EAM::forceField_ = NULL;
57 <  map<int, AtomType*> EAM::EAMlist;
58 <  map<AtomType*, EAMAtomData> EAM::EAMMap;
59 <  map<pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
60 <
54 >  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
55 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
56    
57 <  EAM* EAM::_instance = NULL;
57 >  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {  
58 >    EAMAdapter ea1 = EAMAdapter(atomType1);
59 >    EAMAdapter ea2 = EAMAdapter(atomType2);
60 >    CubicSpline* z1 = ea1.getZ();
61 >    CubicSpline* z2 = ea2.getZ();
62  
64 –  EAM* EAM::Instance() {
65 –    if (!_instance) {
66 –      _instance = new EAM();
67 –    }
68 –    return _instance;
69 –  }
70 –  
71 –  EAMParam EAM::getEAMParam(AtomType* atomType) {
72 –    
73 –    // Do sanity checking on the AtomType we were passed before
74 –    // building any data structures:
75 –    if (!atomType->isEAM()) {
76 –      sprintf( painCave.errMsg,
77 –               "EAM::getEAMParam was passed an atomType (%s) that does not\n"
78 –               "\tappear to be an embedded atom method (EAM) atom.\n",
79 –               atomType->getName().c_str());
80 –      painCave.severity = OPENMD_ERROR;
81 –      painCave.isFatal = 1;
82 –      simError();
83 –    }
84 –    
85 –    GenericData* data = atomType->getPropertyByName("EAM");
86 –    if (data == NULL) {
87 –      sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
88 –               "\tparameters for atomType %s.\n",
89 –               atomType->getName().c_str());
90 –      painCave.severity = OPENMD_ERROR;
91 –      painCave.isFatal = 1;
92 –      simError();
93 –    }
94 –    
95 –    EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
96 –    if (eamData == NULL) {
97 –      sprintf( painCave.errMsg,
98 –               "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
99 –               "\tatom type %s\n", atomType->getName().c_str());
100 –      painCave.severity = OPENMD_ERROR;
101 –      painCave.isFatal = 1;
102 –      simError();          
103 –    }
104 –    
105 –    return eamData->getData();
106 –  }
107 –
108 –  CubicSpline* EAM::getZ(AtomType* atomType) {    
109 –    EAMParam eamParam = getEAMParam(atomType);
110 –    int nr = eamParam.nr;
111 –    RealType dr = eamParam.dr;
112 –    vector<RealType> rvals;
113 –    
114 –    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
115 –      
116 –    CubicSpline* cs = new CubicSpline();
117 –    cs->addPoints(rvals, eamParam.Z);
118 –    return cs;
119 –  }
120 –
121 –  RealType EAM::getRcut(AtomType* atomType) {    
122 –    EAMParam eamParam = getEAMParam(atomType);
123 –    return eamParam.rcut;
124 –  }
125 –
126 –  CubicSpline* EAM::getRho(AtomType* atomType) {    
127 –    EAMParam eamParam = getEAMParam(atomType);
128 –    int nr = eamParam.nr;
129 –    RealType dr = eamParam.dr;
130 –    vector<RealType> rvals;
131 –    
132 –    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
133 –      
134 –    CubicSpline* cs = new CubicSpline();
135 –    cs->addPoints(rvals, eamParam.rho);
136 –    return cs;
137 –  }
138 –
139 –  CubicSpline* EAM::getF(AtomType* atomType) {    
140 –    EAMParam eamParam = getEAMParam(atomType);
141 –    int nrho = eamParam.nrho;
142 –    RealType drho = eamParam.drho;
143 –    vector<RealType> rhovals;
144 –    vector<RealType> scaledF;
145 –    
146 –    for (int i = 0; i < nrho; i++) {
147 –      rhovals.push_back(i * drho);
148 –      scaledF.push_back( eamParam.F[i] * 23.06054 );
149 –    }
150 –      
151 –    CubicSpline* cs = new CubicSpline();
152 –    cs->addPoints(rhovals, eamParam.F);
153 –    return cs;
154 –  }
155 –  
156 –  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {    
157 –    EAMParam eamParam1 = getEAMParam(atomType1);
158 –    EAMParam eamParam2 = getEAMParam(atomType2);
159 –    CubicSpline* z1 = getZ(atomType1);
160 –    CubicSpline* z2 = getZ(atomType2);
161 –
63      // make the r grid:
64  
65  
66      // we need phi out to the largest value we'll encounter in the radial space;
67      
68      RealType rmax = 0.0;
69 <    rmax = max(rmax, eamParam1.rcut);
70 <    rmax = max(rmax, eamParam1.nr * eamParam1.dr);
69 >    rmax = max(rmax, ea1.getRcut());
70 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
71  
72 <    rmax = max(rmax, eamParam2.rcut);
73 <    rmax = max(rmax, eamParam2.nr * eamParam2.dr);
72 >    rmax = max(rmax, ea2.getRcut());
73 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
74  
75      // use the smallest dr (finest grid) to build our grid:
76  
77 <    RealType dr = min(eamParam1.dr, eamParam2.dr);
77 >    RealType dr = min(ea1.getDr(), ea2.getDr());
78  
79      int nr = int(rmax/dr + 0.5);
80  
# Line 189 | Line 90 | namespace OpenMD {
90  
91      phivals.push_back(0.0);
92  
93 <    for (int i = 1; i < rvals.size(); i++ ) {
93 >    for (unsigned int i = 1; i < rvals.size(); i++ ) {
94        r = rvals[i];
95  
96 <      // only use z(r) if we're inside this atoms cutoff radius, otherwise, we'll use zero for the charge.
97 <      // This effectively means that our phi grid goes out beyond the cutoff of the pair potential
96 >      // only use z(r) if we're inside this atom's cutoff radius,
97 >      // otherwise, we'll use zero for the charge.  This effectively
98 >      // means that our phi grid goes out beyond the cutoff of the
99 >      // pair potential
100  
101 <      zi = r <= eamParam1.rcut ? z1->getValueAt(r) : 0.0;
102 <      zj = r <= eamParam2.rcut ? z2->getValueAt(r) : 0.0;
101 >      zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
102 >      zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
103  
104        phi = 331.999296 * (zi * zj) / r;
105  
# Line 208 | Line 111 | namespace OpenMD {
111      return cs;
112    }
113  
114 <  void EAM::initialize() {
114 >  void EAM::setCutoffRadius( RealType rCut ) {
115 >    eamRcut_ = rCut;
116 >    haveCutoffRadius_ = true;
117 >  }
118  
119 +  void EAM::initialize() {
120      // set up the mixing method:
121      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
122      string EAMMixMeth = fopts.getEAMMixingMethod();
# Line 223 | Line 130 | namespace OpenMD {
130        mixMeth_ = eamUnknown;
131        
132      // find all of the EAM atom Types:
133 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
134 <    ForceField::AtomTypeContainer::MapTypeIterator i;
135 <    AtomType* at;
133 >    EAMtypes.clear();
134 >    EAMtids.clear();
135 >    EAMdata.clear();
136 >    MixingMap.clear();
137 >    nEAM_ = 0;
138 >    
139 >    EAMtids.resize( forceField_->getNAtomType(), -1);
140 >
141 >    set<AtomType*>::iterator at;
142 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
143 >      if ((*at)->isEAM()) nEAM_++;
144 >    }
145 >    EAMdata.resize(nEAM_);
146 >    MixingMap.resize(nEAM_);
147  
148 <    for (at = atomTypes->beginType(i); at != NULL;
149 <         at = atomTypes->nextType(i)) {
232 <      
233 <      if (at->isEAM())
234 <        addType(at);
148 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
149 >      if ((*at)->isEAM()) addType(*at);
150      }
151      
152      // find all of the explicit EAM interactions (setfl):
# Line 285 | Line 200 | namespace OpenMD {
200  
201    void EAM::addType(AtomType* atomType){
202  
203 +    EAMAdapter ea = EAMAdapter(atomType);
204      EAMAtomData eamAtomData;
289 –    
290 –    eamAtomData.rho = getRho(atomType);
291 –    eamAtomData.F = getF(atomType);
292 –    eamAtomData.Z = getZ(atomType);
293 –    eamAtomData.rcut = getRcut(atomType);
205  
206 +    eamAtomData.rho = ea.getRho();
207 +    eamAtomData.F = ea.getF();
208 +    eamAtomData.Z = ea.getZ();
209 +    eamAtomData.rcut = ea.getRcut();
210 +
211      // add it to the map:
212 <    AtomTypeProperties atp = atomType->getATP();    
212 >    int atid = atomType->getIdent();
213 >    int eamtid = EAMtypes.size();
214  
215 <    pair<map<int,AtomType*>::iterator,bool> ret;    
216 <    ret = EAMlist.insert( pair<int, AtomType*>(atp.ident, atomType) );
215 >    pair<set<int>::iterator,bool> ret;    
216 >    ret = EAMtypes.insert( atid );
217      if (ret.second == false) {
218        sprintf( painCave.errMsg,
219                 "EAM already had a previous entry with ident %d\n",
220 <               atp.ident);
220 >               atid);
221        painCave.severity = OPENMD_INFO;
222        painCave.isFatal = 0;
223        simError();        
224      }
225  
226 <    EAMMap[atomType] = eamAtomData;
226 >    EAMtids[atid] = eamtid;
227 >    EAMdata[eamtid] = eamAtomData;
228 >    MixingMap[eamtid].resize(nEAM_);
229      
230      // Now, iterate over all known types and add to the mixing map:
231      
232 <    map<AtomType*, EAMAtomData>::iterator it;
233 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
232 >    std::set<int>::iterator it;
233 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
234        
235 <      AtomType* atype2 = (*it).first;
235 >      int eamtid2 = EAMtids[ (*it) ];
236 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
237  
238        EAMInteractionData mixer;
239        mixer.phi = getPhi(atomType, atype2);
240        mixer.explicitlySet = false;
241  
242 <      pair<AtomType*, AtomType*> key1, key2;
323 <      key1 = make_pair(atomType, atype2);
324 <      key2 = make_pair(atype2, atomType);
242 >      MixingMap[eamtid2].resize( nEAM_ );
243        
244 <      MixingMap[key1] = mixer;
245 <      if (key2 != key1) {
246 <        MixingMap[key2] = mixer;
244 >      MixingMap[eamtid][eamtid2] = mixer;
245 >      if (eamtid2 != eamtid) {
246 >        MixingMap[eamtid2][eamtid] = mixer;
247        }
248      }      
249      return;
# Line 349 | Line 267 | namespace OpenMD {
267      mixer.phi = cs;
268      mixer.explicitlySet = true;
269  
270 <    pair<AtomType*, AtomType*> key1, key2;
271 <    key1 = make_pair(atype1, atype2);
354 <    key2 = make_pair(atype2, atype1);
270 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
271 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
272      
273 <    MixingMap[key1] = mixer;
274 <    if (key2 != key1) {
275 <      MixingMap[key2] = mixer;
273 >    MixingMap[eamtid1][eamtid2] = mixer;
274 >    if (eamtid2 != eamtid1) {
275 >      MixingMap[eamtid2][eamtid1] = mixer;
276      }    
277      return;
278    }
279  
280 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
364 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
280 >  void EAM::calcDensity(InteractionData &idat) {
281      
282      if (!initialized_) initialize();
283      
284 <    EAMAtomData data1 = EAMMap[at1];
285 <    EAMAtomData data2 = EAMMap[at2];
284 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
285 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
286  
287 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
288 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
289 <    return;
287 >    if (haveCutoffRadius_)
288 >      if ( *(idat.rij) > eamRcut_) return;
289 >    
290 >    if ( *(idat.rij) < data1.rcut)
291 >      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
292 >    
293 >      
294 >    if ( *(idat.rij) < data2.rcut)
295 >      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
296 >    
297 >    return;  
298    }
299 <
300 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
301 <                           RealType &dfrhodrho) {
378 <
299 >  
300 >  void EAM::calcFunctional(SelfData &sdat) {
301 >    
302      if (!initialized_) initialize();
303  
304 <    EAMAtomData data1 = EAMMap[at1];
304 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
305          
306 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
306 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
307  
308 <    frho = result.first;
309 <    dfrhodrho = result.second;
308 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
309 >    if (sdat.doParticlePot) {
310 >      *(sdat.particlePot) += *(sdat.frho);
311 >    }
312 >
313      return;
314    }
315  
316  
317 <  void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
392 <                      RealType rij, RealType r2, RealType sw,
393 <                      RealType &vpair, RealType &pot, Vector3d &f1,
394 <                      RealType rho_i, RealType rho_j,
395 <                      RealType dfrhodrho_i, RealType dfrhodrho_j,
396 <                      RealType &fshift_i, RealType &fshift_j) {
317 >  void EAM::calcForce(InteractionData &idat) {
318  
319      if (!initialized_) initialize();
320  
321 <    pair<RealType, RealType> res;
322 <    
323 <    if (rij < eamRcut_) {
321 >    if (haveCutoffRadius_)
322 >      if ( *(idat.rij) > eamRcut_) return;
323 >  
324  
325 <      EAMAtomData data1 = EAMMap[at1];
326 <      EAMAtomData data2 = EAMMap[at2];
327 <
328 <      // get type-specific cutoff radii
329 <
330 <      RealType rci = data1.rcut;
331 <      RealType rcj = data2.rcut;
332 <      
333 <      RealType rha, drha, rhb, drhb;
334 <      RealType pha, dpha, phb, dphb;
335 <      RealType phab, dvpdr;
336 <      RealType drhoidr, drhojdr, dudr;
337 <      
338 <      if (rij < rci) {
339 <        res = data1.rho->getValueAndDerivativeAt(rij);
340 <        rha = res.first;
341 <        drha = res.second;
325 >    int eamtid1 = EAMtids[idat.atid1];
326 >    int eamtid2 = EAMtids[idat.atid2];
327 >    
328 >    EAMAtomData &data1 = EAMdata[eamtid1];
329 >    EAMAtomData &data2 = EAMdata[eamtid2];
330 >    
331 >    // get type-specific cutoff radii
332 >    
333 >    RealType rci = data1.rcut;
334 >    RealType rcj = data2.rcut;
335 >    
336 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
337 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
338 >    RealType phab(0.0), dvpdr(0.0);
339 >    RealType drhoidr, drhojdr, dudr;
340 >    
341 >    if ( *(idat.rij) < rci) {
342 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
343 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
344 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
345 >    }
346 >    
347 >    if ( *(idat.rij) < rcj) {
348 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
349 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
350 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
351 >    }
352  
353 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
354 <        pha = res.first;
355 <        dpha = res.second;
353 >    switch(mixMeth_) {
354 >    case eamJohnson:
355 >      
356 >      if ( *(idat.rij) < rci) {
357 >        phab = phab + 0.5 * (rhb / rha) * pha;
358 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
359 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
360        }
361 <
362 <      if (rij < rcj) {
363 <        res = data2.rho->getValueAndDerivativeAt(rij);
364 <        rhb = res.first;
365 <        drhb = res.second;
366 <
367 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
433 <        phb = res.first;
434 <        dphb = res.second;
361 >      
362 >      
363 >      
364 >      if ( *(idat.rij) < rcj) {
365 >        phab = phab + 0.5 * (rha / rhb) * phb;
366 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
367 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
368        }
436 –
437 –      phab = 0.0;
438 –      dvpdr = 0.0;
439 –
440 –      switch(mixMeth_) {
441 –      case eamJohnson:
442 –      
443 –        if (rij < rci) {
444 –          phab = phab + 0.5 * (rhb / rha) * pha;
445 –          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
446 –                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
447 –        }
448 –
449 –        if (rij < rcj) {
450 –          phab = phab + 0.5 * (rha / rhb) * phb;
451 –          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
452 –                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
453 –        }
454 –
455 –        break;
456 –
457 –      case eamDaw:
458 –        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
459 –        phab = res.first;
460 –        dvpdr = res.second;
461 –
462 –        break;
463 –      case eamUnknown:
464 –      default:
465 –
466 –        sprintf(painCave.errMsg,
467 –                "EAM::calcForce hit a mixing method it doesn't know about!\n"
468 –                );
469 –        painCave.severity = OPENMD_ERROR;
470 –        painCave.isFatal = 1;
471 –        simError();        
472 –          
473 –      }
369        
370 <      drhoidr = drha;
371 <      drhojdr = drhb;
372 <
373 <      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
374 <
375 <      f1 = d * dudr / rij;
370 >      break;
371 >      
372 >    case eamDaw:
373 >      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
374 >      
375 >      break;
376 >    case eamUnknown:
377 >    default:
378 >      
379 >      sprintf(painCave.errMsg,
380 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
381 >              );
382 >      painCave.severity = OPENMD_ERROR;
383 >      painCave.isFatal = 1;
384 >      simError();        
385 >      
386 >    }
387 >    
388 >    drhoidr = drha;
389 >    drhojdr = drhb;
390 >    
391 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
392 >    
393 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
394          
395 <      // particle_pot is the difference between the full potential
396 <      // and the full potential without the presence of a particular
395 >    if (idat.doParticlePot) {
396 >      // particlePot is the difference between the full potential and
397 >      // the full potential without the presence of a particular
398        // particle (atom1).
399        //
400 <      // This reduces the density at other particle locations, so
401 <      // we need to recompute the density at atom2 assuming atom1
402 <      // didn't contribute.  This then requires recomputing the
403 <      // density functional for atom2 as well.
404 <      //
405 <      // Most of the particle_pot heavy lifting comes from the
406 <      // pair interaction, and will be handled by vpair.
407 <    
408 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
409 <      fshift_j = data1.F->getValueAt( rho_j - rha );
496 <
497 <      pot += phab;
498 <
499 <      vpair += phab;
400 >      // This reduces the density at other particle locations, so we
401 >      // need to recompute the density at atom2 assuming atom1 didn't
402 >      // contribute.  This then requires recomputing the density
403 >      // functional for atom2 as well.
404 >      
405 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
406 >        - *(idat.frho2);
407 >      
408 >      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
409 >        - *(idat.frho1);
410      }
501 –
502 –    return;
411      
412 <  }
505 <
506 <
507 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *rij,
508 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
509 <
510 <    if (!initialized_) initialize();
412 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
413      
414 <    AtomType* atype1 = EAMlist[*atid1];
415 <    AtomType* atype2 = EAMlist[*atid2];
414 >    *(idat.vpair) += phab;
415 >  
416 >    return;
417      
515 –    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
516 –
517 –    return;    
418    }
419  
420 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
421 <                                   RealType *dfrhodrho) {
420 >  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
421 >    if (!initialized_) initialize();  
422  
423 <    if (!initialized_) initialize();
423 >    RealType cut = 0.0;
424  
425 <    AtomType* atype1 = EAMlist[*atid1];  
426 <
427 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
425 >    int atid1 = atypes.first->getIdent();
426 >    int atid2 = atypes.second->getIdent();
427 >    int eamtid1 = EAMtids[atid1];
428 >    int eamtid2 = EAMtids[atid2];
429      
430 <    return;    
431 <  }
432 <  RealType EAM::getEAMcut(int *atid1) {
430 >    if (eamtid1 != -1) {
431 >      EAMAtomData data1 = EAMdata[eamtid1];
432 >      cut = data1.rcut;
433 >    }
434  
435 <    if (!initialized_) initialize();
435 >    if (eamtid2 != -1) {
436 >      EAMAtomData data2 = EAMdata[eamtid2];
437 >      if (data2.rcut > cut)
438 >        cut = data2.rcut;
439 >    }
440      
441 <    AtomType* atype1 = EAMlist[*atid1];  
536 <      
537 <    return getRcut(atype1);
441 >    return cut;
442    }
539 –
540 –  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
541 –                        RealType *r2, RealType *sw, RealType *vpair,
542 –                        RealType *pot, RealType *f1, RealType *rho1,
543 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
544 –                        RealType *fshift1, RealType *fshift2) {
545 –
546 –    if (!initialized_) initialize();
547 –    
548 –    AtomType* atype1 = EAMlist[*atid1];
549 –    AtomType* atype2 = EAMlist[*atid2];
550 –    
551 –    Vector3d disp(d[0], d[1], d[2]);
552 –    Vector3d frc(f1[0], f1[1], f1[2]);
553 –    
554 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
555 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
556 –      
557 –    f1[0] = frc.x();
558 –    f1[1] = frc.y();
559 –    f1[2] = frc.z();
560 –
561 –    return;    
562 –  }
563 –  
564 –  void EAM::setCutoffEAM(RealType *thisRcut) {
565 –    eamRcut_ = *thisRcut;
566 –  }
443   }
444  
569 – extern "C" {
570 –  
571 – #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
572 – #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
573 – #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
574 – #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
575 – #define fortranGetEAMcut FC_FUNC(geteamcut, GETEAMCUT)
576 –
577 –  
578 –  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
579 –                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
580 –    
581 –    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(atid1, atid2, rij,
582 –                                                         rho_i_at_j,  
583 –                                                         rho_j_at_i);
584 –  }
585 –  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
586 –                             RealType *dfrhodrho) {  
587 –    
588 –    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(atid1, rho, frho,
589 –                                                           dfrhodrho);
590 –    
591 –  }
592 –  void fortranSetCutoffEAM(RealType *rcut) {
593 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
594 –  }
595 –  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
596 –                        RealType *r2, RealType *sw, RealType *vpair,
597 –                        RealType *pot, RealType *f1, RealType *rho1,
598 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
599 –                        RealType *fshift1, RealType *fshift2){
600 –    
601 –    return OpenMD::EAM::Instance()->do_eam_pair(atid1, atid2, d, rij,
602 –                                                r2, sw, vpair,
603 –                                                pot, f1, rho1,
604 –                                                rho2, dfrho1, dfrho2,
605 –                                                fshift1,  fshift2);
606 –  }
607 –  RealType fortranGetEAMcut(int* atid) {
608 –    return OpenMD::EAM::Instance()->getEAMcut(atid);
609 –  }
610 –
611 – }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines