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 1929 by gezelter, Mon Aug 19 13:12:00 2013 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;
109 >      phi = pre11_ * (zi * zj) / r;
110  
111        phivals.push_back(phi);
112      }
# Line 117 | Line 122 | namespace OpenMD {
122    }
123  
124    void EAM::initialize() {
120 –
125      // set up the mixing method:
126      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
127      string EAMMixMeth = fopts.getEAMMixingMethod();
# Line 131 | Line 135 | namespace OpenMD {
135        mixMeth_ = eamUnknown;
136        
137      // find all of the EAM atom Types:
138 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
139 <    ForceField::AtomTypeContainer::MapTypeIterator i;
140 <    AtomType* at;
138 >    EAMtypes.clear();
139 >    EAMtids.clear();
140 >    EAMdata.clear();
141 >    MixingMap.clear();
142 >    nEAM_ = 0;
143 >    
144 >    EAMtids.resize( forceField_->getNAtomType(), -1);
145  
146 <    for (at = atomTypes->beginType(i); at != NULL;
147 <         at = atomTypes->nextType(i)) {
148 <      
141 <      if (at->isEAM())
142 <        addType(at);
146 >    set<AtomType*>::iterator at;
147 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
148 >      if ((*at)->isEAM()) nEAM_++;
149      }
150 +    EAMdata.resize(nEAM_);
151 +    MixingMap.resize(nEAM_);
152 +
153 +    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
154 +      if ((*at)->isEAM()) addType(*at);
155 +    }
156      
157      // find all of the explicit EAM interactions (setfl):
158      ForceField::NonBondedInteractionTypeContainer* nbiTypes = forceField_->getNonBondedInteractionTypes();
# Line 200 | Line 212 | namespace OpenMD {
212      eamAtomData.F = ea.getF();
213      eamAtomData.Z = ea.getZ();
214      eamAtomData.rcut = ea.getRcut();
215 <
215 >      
216      // add it to the map:
217 +    int atid = atomType->getIdent();
218 +    int eamtid = EAMtypes.size();
219  
220 <    pair<map<int,AtomType*>::iterator,bool> ret;    
221 <    ret = EAMlist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
220 >    pair<set<int>::iterator,bool> ret;    
221 >    ret = EAMtypes.insert( atid );
222      if (ret.second == false) {
223        sprintf( painCave.errMsg,
224                 "EAM already had a previous entry with ident %d\n",
225 <               atomType->getIdent());
225 >               atid);
226        painCave.severity = OPENMD_INFO;
227        painCave.isFatal = 0;
228        simError();        
229      }
230  
231 <    EAMMap[atomType] = eamAtomData;
231 >
232 >    EAMtids[atid] = eamtid;
233 >    EAMdata[eamtid] = eamAtomData;
234 >    MixingMap[eamtid].resize(nEAM_);
235      
236      // Now, iterate over all known types and add to the mixing map:
237      
238 <    map<AtomType*, EAMAtomData>::iterator it;
239 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
238 >    std::set<int>::iterator it;
239 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
240        
241 <      AtomType* atype2 = (*it).first;
241 >      int eamtid2 = EAMtids[ (*it) ];
242 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
243  
244        EAMInteractionData mixer;
245        mixer.phi = getPhi(atomType, atype2);
246 +      mixer.rcut = mixer.phi->getLimits().second;
247        mixer.explicitlySet = false;
248  
249 <      pair<AtomType*, AtomType*> key1, key2;
231 <      key1 = make_pair(atomType, atype2);
232 <      key2 = make_pair(atype2, atomType);
249 >      MixingMap[eamtid2].resize( nEAM_ );
250        
251 <      MixingMap[key1] = mixer;
252 <      if (key2 != key1) {
253 <        MixingMap[key2] = mixer;
251 >      MixingMap[eamtid][eamtid2] = mixer;
252 >      if (eamtid2 != eamtid) {
253 >        MixingMap[eamtid2][eamtid] = mixer;
254        }
255      }      
256      return;
# Line 255 | Line 272 | namespace OpenMD {
272  
273      cs->addPoints(rVals, phiVals);
274      mixer.phi = cs;
275 +    mixer.rcut = mixer.phi->getLimits().second;
276      mixer.explicitlySet = true;
277  
278 <    pair<AtomType*, AtomType*> key1, key2;
279 <    key1 = make_pair(atype1, atype2);
262 <    key2 = make_pair(atype2, atype1);
278 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
279 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
280      
281 <    MixingMap[key1] = mixer;
282 <    if (key2 != key1) {
283 <      MixingMap[key2] = mixer;
281 >    MixingMap[eamtid1][eamtid2] = mixer;
282 >    if (eamtid2 != eamtid1) {
283 >      MixingMap[eamtid2][eamtid1] = mixer;
284      }    
285      return;
286    }
# Line 272 | Line 289 | namespace OpenMD {
289      
290      if (!initialized_) initialize();
291      
292 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
293 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
294 <    
292 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
293 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
294 >
295      if (haveCutoffRadius_)
296        if ( *(idat.rij) > eamRcut_) return;
297      
298 <    if ( *(idat.rij) < data1.rcut)
299 <      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
300 <    
298 >    if ( *(idat.rij) < data1.rcut) {
299 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
300 >    }
301        
302 <    if ( *(idat.rij) < data2.rcut)
303 <      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
302 >    if ( *(idat.rij) < data2.rcut) {
303 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
304 >    }
305      
306      return;  
307    }
# Line 292 | Line 310 | namespace OpenMD {
310      
311      if (!initialized_) initialize();
312  
313 <    EAMAtomData data1 = EAMMap[ sdat.atype ];
314 <        
315 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
313 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
314 >            
315 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
316  
317 <    *(sdat.frho) = result.first;
300 <    *(sdat.dfrhodrho) = result.second;
301 <
302 <    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
317 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
318      if (sdat.doParticlePot) {
319 <      *(sdat.particlePot) += result.first;
319 >      *(sdat.particlePot) += *(sdat.frho);
320      }
321  
322      return;
# Line 315 | Line 330 | namespace OpenMD {
330      if (haveCutoffRadius_)
331        if ( *(idat.rij) > eamRcut_) return;
332    
333 <    pair<RealType, RealType> res;
333 >
334 >    int eamtid1 = EAMtids[idat.atid1];
335 >    int eamtid2 = EAMtids[idat.atid2];
336      
337 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
338 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
337 >    EAMAtomData &data1 = EAMdata[eamtid1];
338 >    EAMAtomData &data2 = EAMdata[eamtid2];
339      
340      // get type-specific cutoff radii
341      
# Line 331 | Line 348 | namespace OpenMD {
348      RealType drhoidr, drhojdr, dudr;
349      
350      if ( *(idat.rij) < rci) {
351 <      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
352 <      rha = res.first;
353 <      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;
351 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
352 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
353 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
354      }
355      
356      if ( *(idat.rij) < rcj) {
357 <      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
358 <      rhb = res.first;
359 <      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;
357 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
358 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
359 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
360      }
361  
362      switch(mixMeth_) {
# Line 370 | Line 379 | namespace OpenMD {
379        break;
380        
381      case eamDaw:
373 –      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
374 –      phab = res.first;
375 –      dvpdr = res.second;
382        
383 +      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
384 +        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
385 +                                                                  phab, dvpdr);
386 +      }
387 +      
388        break;
389      case eamUnknown:
390      default:
# Line 393 | Line 404 | namespace OpenMD {
404      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
405      
406      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
407 +
408          
409      if (idat.doParticlePot) {
410        // particlePot is the difference between the full potential and
# Line 424 | Line 436 | namespace OpenMD {
436  
437      RealType cut = 0.0;
438  
439 <    map<AtomType*, EAMAtomData>::iterator it;
440 <
441 <    it = EAMMap.find(atypes.first);
442 <    if (it != EAMMap.end()) {
443 <      EAMAtomData data1 = (*it).second;
439 >    int atid1 = atypes.first->getIdent();
440 >    int atid2 = atypes.second->getIdent();
441 >    int eamtid1 = EAMtids[atid1];
442 >    int eamtid2 = EAMtids[atid2];
443 >    
444 >    if (eamtid1 != -1) {
445 >      EAMAtomData data1 = EAMdata[eamtid1];
446        cut = data1.rcut;
447      }
448  
449 <    it = EAMMap.find(atypes.second);
450 <    if (it != EAMMap.end()) {
437 <      EAMAtomData data2 = (*it).second;
449 >    if (eamtid2 != -1) {
450 >      EAMAtomData data2 = EAMdata[eamtid2];
451        if (data2.rcut > cut)
452          cut = data2.rcut;
453      }
454 <
454 >    
455      return cut;
456    }
457   }

Diff Legend

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