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 1928 by gezelter, Sat Aug 17 13:03:17 2013 UTC vs.
Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 UTC

# Line 51 | Line 51 | namespace OpenMD {
51  
52   namespace OpenMD {
53  
54 <  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
55 <               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
54 >  EAM::EAM() : initialized_(false), haveCutoffRadius_(false),
55 >               forceField_(NULL), eamRcut_(0.0), mixMeth_(eamJohnson),
56 >               name_("EAM") {}
57    
58    CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {  
59      EAMAdapter ea1 = EAMAdapter(atomType1);
# Line 107 | Line 108 | namespace OpenMD {
108        zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
109  
110        phi = pre11_ * (zi * zj) / r;
111 <
111 >      
112        phivals.push_back(phi);
113      }
113 –      
114      CubicSpline* cs = new CubicSpline();
115      cs->addPoints(rvals, phivals);
116      return cs;
# Line 212 | Line 212 | namespace OpenMD {
212      eamAtomData.F = ea.getF();
213      eamAtomData.Z = ea.getZ();
214      eamAtomData.rcut = ea.getRcut();
215 –    eamAtomData.isFluctuating = atomType->isFluctuatingCharge();
215        
216      // add it to the map:
217      int atid = atomType->getIdent();
# Line 227 | Line 226 | namespace OpenMD {
226        painCave.severity = OPENMD_INFO;
227        painCave.isFatal = 0;
228        simError();        
230 –    }
231 –
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;
229      }
230  
231  
# Line 314 | Line 296 | namespace OpenMD {
296        if ( *(idat.rij) > eamRcut_) return;
297      
298      if ( *(idat.rij) < data1.rcut) {
299 <      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 <      }
299 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
300      }
301        
302      if ( *(idat.rij) < data2.rcut) {
303 <      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 <      }
303 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
304      }
305      
306      return;  
# Line 337 | Line 309 | namespace OpenMD {
309    void EAM::calcFunctional(SelfData &sdat) {
310      
311      if (!initialized_) initialize();
340 –
312      EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
313              
314 <    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
314 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho),
315 >                                      *(sdat.dfrhodrho) );
316  
317      (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
318      if (sdat.doParticlePot) {
# Line 361 | Line 333 | namespace OpenMD {
333  
334      int eamtid1 = EAMtids[idat.atid1];
335      int eamtid2 = EAMtids[idat.atid2];
364 –    
336      EAMAtomData &data1 = EAMdata[eamtid1];
337      EAMAtomData &data2 = EAMdata[eamtid2];
338      
# Line 369 | Line 340 | namespace OpenMD {
340      
341      RealType rci = data1.rcut;
342      RealType rcj = data2.rcut;
343 +
344      
345      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
346      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
# Line 379 | Line 351 | namespace OpenMD {
351        data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
352        CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
353        phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
382 –      if (data1.isFluctuating) {
383 –        *(idat.dVdFQ1) -= *(idat.dfrho2) * rha * data1.qToRhoScaling;
384 –      }
354      }
355      
356      if ( *(idat.rij) < rcj) {
357        data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
358        CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
359        phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
391 –      if (data2.isFluctuating) {
392 –        *(idat.dVdFQ2) -= *(idat.dfrho1) * rhb * data2.qToRhoScaling;
393 –      }
360      }
395 –
361      switch(mixMeth_) {
362      case eamJohnson:
398 –      
363        if ( *(idat.rij) < rci) {
364          phab = phab + 0.5 * (rhb / rha) * pha;
365          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
366                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
367        }
368 <      
405 <      
406 <      
368 >                
369        if ( *(idat.rij) < rcj) {
370          phab = phab + 0.5 * (rha / rhb) * phb;
371          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
372                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
373        }
412 –      
374        break;
414 –      
375      case eamDaw:
416 –      
376        if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
377          MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
378                                                                    phab, dvpdr);
379        }
421 –      
380        break;
381      case eamUnknown:
382      default:
# Line 439 | Line 397 | namespace OpenMD {
397      
398      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
399  
442 –        
400      if (idat.doParticlePot) {
401        // particlePot is the difference between the full potential and
402        // the full potential without the presence of a particular
# Line 457 | Line 414 | namespace OpenMD {
414          - *(idat.frho1);
415      }
416      
417 <    (*(idat.pot))[METALLIC_FAMILY] += phab;
418 <    
419 <    *(idat.vpair) += phab;
463 <  
464 <    return;
465 <    
417 >    (*(idat.pot))[METALLIC_FAMILY] += phab;    
418 >    *(idat.vpair) += phab;  
419 >    return;    
420    }
421  
422    RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {

Diff Legend

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