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 1928 by gezelter, Sat Aug 17 13:03:17 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 >    eamAtomData.isFluctuating = atomType->isFluctuatingCharge();
216 >      
217      // add it to the map:
218 +    int atid = atomType->getIdent();
219 +    int eamtid = EAMtypes.size();
220  
221 <    pair<map<int,AtomType*>::iterator,bool> ret;    
222 <    ret = EAMlist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
221 >    pair<set<int>::iterator,bool> ret;    
222 >    ret = EAMtypes.insert( atid );
223      if (ret.second == false) {
224        sprintf( painCave.errMsg,
225                 "EAM already had a previous entry with ident %d\n",
226 <               atomType->getIdent());
226 >               atid);
227        painCave.severity = OPENMD_INFO;
228        painCave.isFatal = 0;
229        simError();        
230      }
231  
232 <    EAMMap[atomType] = eamAtomData;
232 >    if (eamAtomData.isFluctuating) {
233 >      // compute charge to rho scaling:
234 >      RealType z0 = eamAtomData.Z->getValueAt(0.0);
235 >      RealType dr = ea.getDr();
236 >      RealType rmax = max(eamAtomData.rcut, ea.getNr() * dr);
237 >      int nr = int(rmax/dr + 0.5);
238 >      RealType r;
239 >      RealType sum(0.0);
240 >
241 >      for (int i = 0; i < nr; i++) {
242 >        r = RealType(i*dr);
243 >        sum += r * r * eamAtomData.rho->getValueAt(r) * dr;      
244 >      }
245 >      sum *= 4.0 * M_PI;
246 >      eamAtomData.qToRhoScaling = sum / z0;
247 >    }
248 >
249 >
250 >    EAMtids[atid] = eamtid;
251 >    EAMdata[eamtid] = eamAtomData;
252 >    MixingMap[eamtid].resize(nEAM_);
253      
254      // Now, iterate over all known types and add to the mixing map:
255      
256 <    map<AtomType*, EAMAtomData>::iterator it;
257 <    for( it = EAMMap.begin(); it != EAMMap.end(); ++it) {
256 >    std::set<int>::iterator it;
257 >    for( it = EAMtypes.begin(); it != EAMtypes.end(); ++it) {
258        
259 <      AtomType* atype2 = (*it).first;
259 >      int eamtid2 = EAMtids[ (*it) ];
260 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
261  
262        EAMInteractionData mixer;
263        mixer.phi = getPhi(atomType, atype2);
264 +      mixer.rcut = mixer.phi->getLimits().second;
265        mixer.explicitlySet = false;
266  
267 <      pair<AtomType*, AtomType*> key1, key2;
231 <      key1 = make_pair(atomType, atype2);
232 <      key2 = make_pair(atype2, atomType);
267 >      MixingMap[eamtid2].resize( nEAM_ );
268        
269 <      MixingMap[key1] = mixer;
270 <      if (key2 != key1) {
271 <        MixingMap[key2] = mixer;
269 >      MixingMap[eamtid][eamtid2] = mixer;
270 >      if (eamtid2 != eamtid) {
271 >        MixingMap[eamtid2][eamtid] = mixer;
272        }
273      }      
274      return;
# Line 255 | Line 290 | namespace OpenMD {
290  
291      cs->addPoints(rVals, phiVals);
292      mixer.phi = cs;
293 +    mixer.rcut = mixer.phi->getLimits().second;
294      mixer.explicitlySet = true;
295  
296 <    pair<AtomType*, AtomType*> key1, key2;
297 <    key1 = make_pair(atype1, atype2);
262 <    key2 = make_pair(atype2, atype1);
296 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
297 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
298      
299 <    MixingMap[key1] = mixer;
300 <    if (key2 != key1) {
301 <      MixingMap[key2] = mixer;
299 >    MixingMap[eamtid1][eamtid2] = mixer;
300 >    if (eamtid2 != eamtid1) {
301 >      MixingMap[eamtid2][eamtid1] = mixer;
302      }    
303      return;
304    }
# Line 272 | Line 307 | namespace OpenMD {
307      
308      if (!initialized_) initialize();
309      
310 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
311 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
312 <    
310 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
311 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
312 >
313      if (haveCutoffRadius_)
314        if ( *(idat.rij) > eamRcut_) return;
315      
316 <    if ( *(idat.rij) < data1.rcut)
317 <      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
318 <    
316 >    if ( *(idat.rij) < data1.rcut) {
317 >      if (data1.isFluctuating) {
318 >        *(idat.rho2) += (1.0 -  *(idat.flucQ1) * data1.qToRhoScaling ) *
319 >          data1.rho->getValueAt( *(idat.rij) );
320 >      } else {
321 >        *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
322 >      }
323 >    }
324        
325 <    if ( *(idat.rij) < data2.rcut)
326 <      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
325 >    if ( *(idat.rij) < data2.rcut) {
326 >      if (data2.isFluctuating) {
327 >        *(idat.rho1) += (1.0 -  *(idat.flucQ2) * data2.qToRhoScaling ) *
328 >          data2.rho->getValueAt( *(idat.rij) );
329 >      } else {
330 >        *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
331 >      }
332 >    }
333      
334      return;  
335    }
# Line 292 | Line 338 | namespace OpenMD {
338      
339      if (!initialized_) initialize();
340  
341 <    EAMAtomData data1 = EAMMap[ sdat.atype ];
342 <        
343 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
341 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
342 >            
343 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
344  
345 <    *(sdat.frho) = result.first;
300 <    *(sdat.dfrhodrho) = result.second;
301 <
302 <    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
345 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
346      if (sdat.doParticlePot) {
347 <      *(sdat.particlePot) += result.first;
347 >      *(sdat.particlePot) += *(sdat.frho);
348      }
349  
350      return;
# Line 315 | Line 358 | namespace OpenMD {
358      if (haveCutoffRadius_)
359        if ( *(idat.rij) > eamRcut_) return;
360    
361 <    pair<RealType, RealType> res;
361 >
362 >    int eamtid1 = EAMtids[idat.atid1];
363 >    int eamtid2 = EAMtids[idat.atid2];
364      
365 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
366 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
365 >    EAMAtomData &data1 = EAMdata[eamtid1];
366 >    EAMAtomData &data2 = EAMdata[eamtid2];
367      
368      // get type-specific cutoff radii
369      
# Line 331 | Line 376 | namespace OpenMD {
376      RealType drhoidr, drhojdr, dudr;
377      
378      if ( *(idat.rij) < rci) {
379 <      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
380 <      rha = res.first;
381 <      drha = res.second;
382 <      
383 <      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
384 <      pha = res.first;
340 <      dpha = res.second;
379 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
380 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
381 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
382 >      if (data1.isFluctuating) {
383 >        *(idat.dVdFQ1) -= *(idat.dfrho2) * rha * data1.qToRhoScaling;
384 >      }
385      }
386      
387      if ( *(idat.rij) < rcj) {
388 <      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
389 <      rhb = res.first;
390 <      drhb = res.second;
391 <      
392 <      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
393 <      phb = res.first;
350 <      dphb = res.second;
388 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
389 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
390 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
391 >      if (data2.isFluctuating) {
392 >        *(idat.dVdFQ2) -= *(idat.dfrho1) * rhb * data2.qToRhoScaling;
393 >      }
394      }
395  
396      switch(mixMeth_) {
# Line 370 | Line 413 | namespace OpenMD {
413        break;
414        
415      case eamDaw:
373 –      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
374 –      phab = res.first;
375 –      dvpdr = res.second;
416        
417 +      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
418 +        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
419 +                                                                  phab, dvpdr);
420 +      }
421 +      
422        break;
423      case eamUnknown:
424      default:
# Line 393 | Line 438 | namespace OpenMD {
438      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
439      
440      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
441 +
442          
443      if (idat.doParticlePot) {
444        // particlePot is the difference between the full potential and
# Line 424 | Line 470 | namespace OpenMD {
470  
471      RealType cut = 0.0;
472  
473 <    map<AtomType*, EAMAtomData>::iterator it;
474 <
475 <    it = EAMMap.find(atypes.first);
476 <    if (it != EAMMap.end()) {
477 <      EAMAtomData data1 = (*it).second;
473 >    int atid1 = atypes.first->getIdent();
474 >    int atid2 = atypes.second->getIdent();
475 >    int eamtid1 = EAMtids[atid1];
476 >    int eamtid2 = EAMtids[atid2];
477 >    
478 >    if (eamtid1 != -1) {
479 >      EAMAtomData data1 = EAMdata[eamtid1];
480        cut = data1.rcut;
481      }
482  
483 <    it = EAMMap.find(atypes.second);
484 <    if (it != EAMMap.end()) {
437 <      EAMAtomData data2 = (*it).second;
483 >    if (eamtid2 != -1) {
484 >      EAMAtomData data2 = EAMdata[eamtid2];
485        if (data2.rcut > cut)
486          cut = data2.rcut;
487      }
488 <
488 >    
489      return cut;
490    }
491   }

Diff Legend

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