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 1927 by gezelter, Wed Aug 14 20:19:19 2013 UTC vs.
Revision 2033 by gezelter, Sat Nov 1 14:12:16 2014 UTC

# Line 107 | Line 107 | namespace OpenMD {
107        zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
108  
109        phi = pre11_ * (zi * zj) / r;
110 <
110 >      
111        phivals.push_back(phi);
112      }
113 –      
113      CubicSpline* cs = new CubicSpline();
114      cs->addPoints(rvals, phivals);
115      return cs;
# Line 212 | Line 211 | namespace OpenMD {
211      eamAtomData.F = ea.getF();
212      eamAtomData.Z = ea.getZ();
213      eamAtomData.rcut = ea.getRcut();
215 –    eamAtomData.isFluctuating = atomType->isFluctuatingCharge();
214        
215      // add it to the map:
216      int atid = atomType->getIdent();
# Line 227 | Line 225 | namespace OpenMD {
225        painCave.severity = OPENMD_INFO;
226        painCave.isFatal = 0;
227        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;
228      }
229  
230  
# Line 261 | 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 289 | 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 312 | Line 295 | namespace OpenMD {
295        if ( *(idat.rij) > eamRcut_) return;
296      
297      if ( *(idat.rij) < data1.rcut) {
298 <      if (data1.isFluctuating) {
316 <        *(idat.rho2) += (1.0 -  *(idat.flucQ1) * data1.qToRhoScaling ) *
317 <          data1.rho->getValueAt( *(idat.rij) );
318 <      } else {
319 <        *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
320 <      }
298 >      *(idat.rho2) += data1.rho->getValueAt( *(idat.rij));
299      }
300        
301      if ( *(idat.rij) < data2.rcut) {
302 <      if (data2.isFluctuating) {
325 <        *(idat.rho1) += (1.0 -  *(idat.flucQ2) * data2.qToRhoScaling ) *
326 <          data2.rho->getValueAt( *(idat.rij) );
327 <      } else {
328 <        *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
329 <      }
302 >      *(idat.rho1) += data2.rho->getValueAt( *(idat.rij));
303      }
304      
305      return;  
# Line 335 | Line 308 | namespace OpenMD {
308    void EAM::calcFunctional(SelfData &sdat) {
309      
310      if (!initialized_) initialize();
338 –
311      EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
312              
313 <    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
313 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho),
314 >                                      *(sdat.dfrhodrho) );
315  
316      (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
317      if (sdat.doParticlePot) {
# Line 359 | Line 332 | namespace OpenMD {
332  
333      int eamtid1 = EAMtids[idat.atid1];
334      int eamtid2 = EAMtids[idat.atid2];
362 –    
335      EAMAtomData &data1 = EAMdata[eamtid1];
336      EAMAtomData &data2 = EAMdata[eamtid2];
337      
# Line 367 | Line 339 | namespace OpenMD {
339      
340      RealType rci = data1.rcut;
341      RealType rcj = data2.rcut;
342 +
343      
344      RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
345      RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
# Line 377 | Line 350 | namespace OpenMD {
350        data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
351        CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
352        phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
380 –      if (data1.isFluctuating) {
381 –        *(idat.dVdFQ1) -= *(idat.dfrho2) * rha * data1.qToRhoScaling;
382 –      }
353      }
354      
355      if ( *(idat.rij) < rcj) {
356        data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
357        CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
358        phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
389 –      if (data2.isFluctuating) {
390 –        *(idat.dVdFQ2) -= *(idat.dfrho1) * rhb * data2.qToRhoScaling;
391 –      }
359      }
393 –
360      switch(mixMeth_) {
361      case eamJohnson:
396 –      
362        if ( *(idat.rij) < rci) {
363          phab = phab + 0.5 * (rhb / rha) * pha;
364          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
365                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
366        }
367 <      
403 <      
404 <      
367 >                
368        if ( *(idat.rij) < rcj) {
369          phab = phab + 0.5 * (rha / rhb) * phb;
370          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
371                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
372        }
410 –      
373        break;
412 –      
374      case eamDaw:
375 <      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
376 <      
375 >      if ( *(idat.rij) <  MixingMap[eamtid1][eamtid2].rcut) {
376 >        MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij),
377 >                                                                  phab, dvpdr);
378 >      }
379        break;
380      case eamUnknown:
381      default:
# Line 433 | Line 396 | namespace OpenMD {
396      
397      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
398  
436 –        
399      if (idat.doParticlePot) {
400        // particlePot is the difference between the full potential and
401        // the full potential without the presence of a particular
# Line 451 | Line 413 | namespace OpenMD {
413          - *(idat.frho1);
414      }
415      
416 <    (*(idat.pot))[METALLIC_FAMILY] += phab;
417 <    
418 <    *(idat.vpair) += phab;
457 <  
458 <    return;
459 <    
416 >    (*(idat.pot))[METALLIC_FAMILY] += phab;    
417 >    *(idat.vpair) += phab;  
418 >    return;    
419    }
420  
421    RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {

Diff Legend

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