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 1895 by gezelter, Mon Jul 1 21:09:37 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 207 | 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();
# Line 223 | Line 227 | namespace OpenMD {
227        simError();        
228      }
229  
230 +
231      EAMtids[atid] = eamtid;
232      EAMdata[eamtid] = eamAtomData;
233      MixingMap[eamtid].resize(nEAM_);
# Line 237 | Line 242 | namespace OpenMD {
242  
243        EAMInteractionData mixer;
244        mixer.phi = getPhi(atomType, atype2);
245 +      mixer.rcut = mixer.phi->getLimits().second;
246        mixer.explicitlySet = false;
247  
248        MixingMap[eamtid2].resize( nEAM_ );
# Line 265 | 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      int eamtid1 = EAMtids[ atype1->getIdent() ];
# Line 287 | Line 294 | namespace OpenMD {
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 300 | Line 308 | namespace OpenMD {
308    void EAM::calcFunctional(SelfData &sdat) {
309      
310      if (!initialized_) initialize();
303 –
311      EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
312 <        
312 >            
313      data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
314  
315      (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
# Line 324 | Line 331 | namespace OpenMD {
331  
332      int eamtid1 = EAMtids[idat.atid1];
333      int eamtid2 = EAMtids[idat.atid2];
327 –    
334      EAMAtomData &data1 = EAMdata[eamtid1];
335      EAMAtomData &data2 = EAMdata[eamtid2];
336      
# Line 332 | Line 338 | namespace OpenMD {
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 349 | Line 356 | namespace OpenMD {
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 <      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
377 <      
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 391 | 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

Diff Legend

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