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

Comparing trunk/src/nonbonded/EAM.cpp (file contents):
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC vs.
Revision 2031 by jmichalk, Fri Oct 31 18:40:40 2014 UTC

# Line 60 | Line 60 | namespace OpenMD {
60      CubicSpline* z1 = ea1.getZ();
61      CubicSpline* z2 = ea2.getZ();
62  
63 +    // Thise prefactors convert the charge-charge interactions into
64 +    // kcal / mol all were computed assuming distances are measured in
65 +    // angstroms Charge-Charge, assuming charges are measured in
66 +    // electrons.  Matches value in Electrostatics.cpp
67 +    pre11_ = 332.0637778;
68 +
69      // make the r grid:
70  
65 –
71      // we need phi out to the largest value we'll encounter in the radial space;
72      
73      RealType rmax = 0.0;
# Line 101 | Line 106 | namespace OpenMD {
106        zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
107        zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
108  
109 <      phi = 331.999296 * (zi * zj) / r;
110 <
109 >      phi = pre11_ * (zi * zj) / r;
110 >      
111        phivals.push_back(phi);
112      }
108 –      
113      CubicSpline* cs = new CubicSpline();
114      cs->addPoints(rvals, phivals);
115      return cs;
# Line 117 | Line 121 | namespace OpenMD {
121    }
122  
123    void EAM::initialize() {
120 –
124      // set up the mixing method:
125      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
126      string EAMMixMeth = fopts.getEAMMixingMethod();
# Line 131 | Line 134 | namespace OpenMD {
134        mixMeth_ = eamUnknown;
135        
136      // find all of the EAM atom Types:
137 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
138 <    ForceField::AtomTypeContainer::MapTypeIterator i;
139 <    AtomType* at;
137 >    EAMtypes.clear();
138 >    EAMtids.clear();
139 >    EAMdata.clear();
140 >    MixingMap.clear();
141 >    nEAM_ = 0;
142 >    
143 >    EAMtids.resize( forceField_->getNAtomType(), -1);
144  
145 <    for (at = atomTypes->beginType(i); at != NULL;
146 <         at = atomTypes->nextType(i)) {
147 <      
141 <      if (at->isEAM())
142 <        addType(at);
145 >    set<AtomType*>::iterator at;
146 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
147 >      if ((*at)->isEAM()) nEAM_++;
148      }
149 +    EAMdata.resize(nEAM_);
150 +    MixingMap.resize(nEAM_);
151 +
152 +    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
153 +      if ((*at)->isEAM()) addType(*at);
154 +    }
155      
156      // find all of the explicit EAM interactions (setfl):
157      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
# Line 200 | Line 211 | namespace OpenMD {
211      eamAtomData.F = ea.getF();
212      eamAtomData.Z = ea.getZ();
213      eamAtomData.rcut = ea.getRcut();
214 <
214 >      
215      // add it to the map:
216 +    int atid = atomType->getIdent();
217 +    int eamtid = EAMtypes.size();
218  
219 <    pair<map<int,AtomType*>::iterator,bool> ret;    
220 <    ret = EAMlist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
219 >    pair<set<int>::iterator,bool> ret;    
220 >    ret = EAMtypes.insert( atid );
221      if (ret.second == false) {
222        sprintf( painCave.errMsg,
223                 "EAM already had a previous entry with ident %d\n",
224 <               atomType->getIdent());
224 >               atid);
225        painCave.severity = OPENMD_INFO;
226        painCave.isFatal = 0;
227        simError();        
228      }
229  
230 <    EAMMap[atomType] = eamAtomData;
230 >
231 >    EAMtids[atid] = eamtid;
232 >    EAMdata[eamtid] = eamAtomData;
233 >    MixingMap[eamtid].resize(nEAM_);
234      
235      // Now, iterate over all known types and add to the mixing map:
236      
237 <    map<AtomType*, EAMAtomData>::iterator it;
238 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
237 >    std::set<int>::iterator it;
238 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
239        
240 <      AtomType* atype2 = (*it).first;
240 >      int eamtid2 = EAMtids[ (*it) ];
241 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
242  
243        EAMInteractionData mixer;
244        mixer.phi = getPhi(atomType, atype2);
245 +      mixer.rcut = mixer.phi->getLimits().second;
246        mixer.explicitlySet = false;
247  
248 <      pair<AtomType*, AtomType*> key1, key2;
231 <      key1 = make_pair(atomType, atype2);
232 <      key2 = make_pair(atype2, atomType);
248 >      MixingMap[eamtid2].resize( nEAM_ );
249        
250 <      MixingMap[key1] = mixer;
251 <      if (key2 != key1) {
252 <        MixingMap[key2] = mixer;
250 >      MixingMap[eamtid][eamtid2] = mixer;
251 >      if (eamtid2 != eamtid) {
252 >        MixingMap[eamtid2][eamtid] = mixer;
253        }
254      }      
255      return;
# Line 255 | Line 271 | namespace OpenMD {
271  
272      cs->addPoints(rVals, phiVals);
273      mixer.phi = cs;
274 +    mixer.rcut = mixer.phi->getLimits().second;
275      mixer.explicitlySet = true;
276  
277 <    pair<AtomType*, AtomType*> key1, key2;
278 <    key1 = make_pair(atype1, atype2);
262 <    key2 = make_pair(atype2, atype1);
277 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
278 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
279      
280 <    MixingMap[key1] = mixer;
281 <    if (key2 != key1) {
282 <      MixingMap[key2] = mixer;
280 >    MixingMap[eamtid1][eamtid2] = mixer;
281 >    if (eamtid2 != eamtid1) {
282 >      MixingMap[eamtid2][eamtid1] = mixer;
283      }    
284      return;
285    }
# Line 272 | Line 288 | namespace OpenMD {
288      
289      if (!initialized_) initialize();
290      
291 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
292 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
293 <    
291 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
292 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
293 >
294      if (haveCutoffRadius_)
295        if ( *(idat.rij) > eamRcut_) return;
296      
297 <    if ( *(idat.rij) < data1.rcut)
298 <      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
299 <    
297 >    if ( *(idat.rij) < data1.rcut) {
298 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
299 >    }
300        
301 <    if ( *(idat.rij) < data2.rcut)
302 <      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
301 >    if ( *(idat.rij) < data2.rcut) {
302 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
303 >    }
304      
305      return;  
306    }
# Line 291 | Line 308 | namespace OpenMD {
308    void EAM::calcFunctional(SelfData &sdat) {
309      
310      if (!initialized_) initialize();
311 +    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
312 +            
313 +    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
314  
315 <    EAMAtomData data1 = EAMMap[ sdat.atype ];
296 <        
297 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
298 <
299 <    *(sdat.frho) = result.first;
300 <    *(sdat.dfrhodrho) = result.second;
301 <
302 <    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
315 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
316      if (sdat.doParticlePot) {
317 <      *(sdat.particlePot) += result.first;
317 >      *(sdat.particlePot) += *(sdat.frho);
318      }
319  
320      return;
# Line 315 | Line 328 | namespace OpenMD {
328      if (haveCutoffRadius_)
329        if ( *(idat.rij) > eamRcut_) return;
330    
331 <    pair<RealType, RealType> res;
331 >
332 >    int eamtid1 = EAMtids[idat.atid1];
333 >    int eamtid2 = EAMtids[idat.atid2];
334 >    EAMAtomData &data1 = EAMdata[eamtid1];
335 >    EAMAtomData &data2 = EAMdata[eamtid2];
336      
320 –    EAMAtomData data1 = EAMMap[idat.atypes.first];
321 –    EAMAtomData data2 = EAMMap[idat.atypes.second];
322 –    
337      // get type-specific cutoff radii
338      
339      RealType rci = data1.rcut;
340      RealType rcj = data2.rcut;
341 +
342      
343      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
344      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
# Line 331 | Line 346 | namespace OpenMD {
346      RealType drhoidr, drhojdr, dudr;
347      
348      if ( *(idat.rij) < rci) {
349 <      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
350 <      rha = res.first;
351 <      drha = res.second;
337 <      
338 <      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
339 <      pha = res.first;
340 <      dpha = res.second;
349 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
350 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
351 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
352      }
353      
354      if ( *(idat.rij) < rcj) {
355 <      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
356 <      rhb = res.first;
357 <      drhb = res.second;
347 <      
348 <      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
349 <      phb = res.first;
350 <      dphb = res.second;
355 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
356 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
357 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
358      }
352 –
359      switch(mixMeth_) {
360      case eamJohnson:
355 –      
361        if ( *(idat.rij) < rci) {
362          phab = phab + 0.5 * (rhb / rha) * pha;
363          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
# Line 366 | Line 371 | namespace OpenMD {
371          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
372                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
373        }
369 –      
374        break;
371 –      
375      case eamDaw:
376 <      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
377 <      phab = res.first;
378 <      dvpdr = res.second;
379 <      
376 >      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
377 >        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
378 >                                                                  phab, dvpdr);
379 >      }
380        break;
381      case eamUnknown:
382      default:
# Line 393 | Line 396 | namespace OpenMD {
396      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
397      
398      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
399 +
400          
401      if (idat.doParticlePot) {
402        // particlePot is the difference between the full potential and
# Line 424 | Line 428 | namespace OpenMD {
428  
429      RealType cut = 0.0;
430  
431 <    map<AtomType*, EAMAtomData>::iterator it;
432 <
433 <    it = EAMMap.find(atypes.first);
434 <    if (it != EAMMap.end()) {
435 <      EAMAtomData data1 = (*it).second;
431 >    int atid1 = atypes.first->getIdent();
432 >    int atid2 = atypes.second->getIdent();
433 >    int eamtid1 = EAMtids[atid1];
434 >    int eamtid2 = EAMtids[atid2];
435 >    
436 >    if (eamtid1 != -1) {
437 >      EAMAtomData data1 = EAMdata[eamtid1];
438        cut = data1.rcut;
439      }
440  
441 <    it = EAMMap.find(atypes.second);
442 <    if (it != EAMMap.end()) {
437 <      EAMAtomData data2 = (*it).second;
441 >    if (eamtid2 != -1) {
442 >      EAMAtomData data2 = EAMdata[eamtid2];
443        if (data2.rcut > cut)
444          cut = data2.rcut;
445      }
446 <
446 >    
447      return cut;
448    }
449   }

Diff Legend

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