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 1927 by gezelter, Wed Aug 14 20:19:19 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.explicitlySet = false;
265  
266 <      pair<AtomType*, AtomType*> key1, key2;
231 <      key1 = make_pair(atomType, atype2);
232 <      key2 = make_pair(atype2, atomType);
266 >      MixingMap[eamtid2].resize( nEAM_ );
267        
268 <      MixingMap[key1] = mixer;
269 <      if (key2 != key1) {
270 <        MixingMap[key2] = mixer;
268 >      MixingMap[eamtid][eamtid2] = mixer;
269 >      if (eamtid2 != eamtid) {
270 >        MixingMap[eamtid2][eamtid] = mixer;
271        }
272      }      
273      return;
# Line 257 | Line 291 | namespace OpenMD {
291      mixer.phi = cs;
292      mixer.explicitlySet = true;
293  
294 <    pair<AtomType*, AtomType*> key1, key2;
295 <    key1 = make_pair(atype1, atype2);
262 <    key2 = make_pair(atype2, atype1);
294 >    int eamtid1 = EAMtids[ atype1->getIdent() ];
295 >    int eamtid2 = EAMtids[ atype2->getIdent() ];
296      
297 <    MixingMap[key1] = mixer;
298 <    if (key2 != key1) {
299 <      MixingMap[key2] = mixer;
297 >    MixingMap[eamtid1][eamtid2] = mixer;
298 >    if (eamtid2 != eamtid1) {
299 >      MixingMap[eamtid2][eamtid1] = mixer;
300      }    
301      return;
302    }
# Line 272 | Line 305 | namespace OpenMD {
305      
306      if (!initialized_) initialize();
307      
308 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
309 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
310 <    
308 >    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
309 >    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
310 >
311      if (haveCutoffRadius_)
312        if ( *(idat.rij) > eamRcut_) return;
313      
314 <    if ( *(idat.rij) < data1.rcut)
315 <      *(idat.rho1) += data1.rho->getValueAt( *(idat.rij));
316 <    
314 >    if ( *(idat.rij) < data1.rcut) {
315 >      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 >      }
321 >    }
322        
323 <    if ( *(idat.rij) < data2.rcut)
324 <      *(idat.rho2) += data2.rho->getValueAt( *(idat.rij));
323 >    if ( *(idat.rij) < data2.rcut) {
324 >      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 >      }
330 >    }
331      
332      return;  
333    }
# Line 292 | Line 336 | namespace OpenMD {
336      
337      if (!initialized_) initialize();
338  
339 <    EAMAtomData data1 = EAMMap[ sdat.atype ];
340 <        
341 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt( *(sdat.rho) );
339 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
340 >            
341 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
342  
343 <    *(sdat.frho) = result.first;
300 <    *(sdat.dfrhodrho) = result.second;
301 <
302 <    (*(sdat.pot))[METALLIC_FAMILY] += result.first;
343 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
344      if (sdat.doParticlePot) {
345 <      *(sdat.particlePot) += result.first;
345 >      *(sdat.particlePot) += *(sdat.frho);
346      }
347  
348      return;
# Line 315 | Line 356 | namespace OpenMD {
356      if (haveCutoffRadius_)
357        if ( *(idat.rij) > eamRcut_) return;
358    
359 <    pair<RealType, RealType> res;
359 >
360 >    int eamtid1 = EAMtids[idat.atid1];
361 >    int eamtid2 = EAMtids[idat.atid2];
362      
363 <    EAMAtomData data1 = EAMMap[idat.atypes.first];
364 <    EAMAtomData data2 = EAMMap[idat.atypes.second];
363 >    EAMAtomData &data1 = EAMdata[eamtid1];
364 >    EAMAtomData &data2 = EAMdata[eamtid2];
365      
366      // get type-specific cutoff radii
367      
# Line 331 | Line 374 | namespace OpenMD {
374      RealType drhoidr, drhojdr, dudr;
375      
376      if ( *(idat.rij) < rci) {
377 <      res = data1.rho->getValueAndDerivativeAt( *(idat.rij));
378 <      rha = res.first;
379 <      drha = res.second;
380 <      
381 <      res = MixingMap[make_pair(idat.atypes.first, idat.atypes.first)].phi->getValueAndDerivativeAt( *(idat.rij) );
382 <      pha = res.first;
340 <      dpha = res.second;
377 >      data1.rho->getValueAndDerivativeAt( *(idat.rij), rha, drha);
378 >      CubicSpline* phi = MixingMap[eamtid1][eamtid1].phi;
379 >      phi->getValueAndDerivativeAt( *(idat.rij), pha, dpha);
380 >      if (data1.isFluctuating) {
381 >        *(idat.dVdFQ1) -= *(idat.dfrho2) * rha * data1.qToRhoScaling;
382 >      }
383      }
384      
385      if ( *(idat.rij) < rcj) {
386 <      res = data2.rho->getValueAndDerivativeAt( *(idat.rij) );
387 <      rhb = res.first;
388 <      drhb = res.second;
389 <      
390 <      res = MixingMap[make_pair(idat.atypes.second, idat.atypes.second)].phi->getValueAndDerivativeAt( *(idat.rij) );
391 <      phb = res.first;
350 <      dphb = res.second;
386 >      data2.rho->getValueAndDerivativeAt( *(idat.rij), rhb, drhb );
387 >      CubicSpline* phi = MixingMap[eamtid2][eamtid2].phi;
388 >      phi->getValueAndDerivativeAt( *(idat.rij), phb, dphb);
389 >      if (data2.isFluctuating) {
390 >        *(idat.dVdFQ2) -= *(idat.dfrho1) * rhb * data2.qToRhoScaling;
391 >      }
392      }
393  
394      switch(mixMeth_) {
# Line 370 | Line 411 | namespace OpenMD {
411        break;
412        
413      case eamDaw:
414 <      res = MixingMap[idat.atypes].phi->getValueAndDerivativeAt( *(idat.rij));
374 <      phab = res.first;
375 <      dvpdr = res.second;
414 >      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
415        
416        break;
417      case eamUnknown:
# Line 393 | Line 432 | namespace OpenMD {
432      dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
433      
434      *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
435 +
436          
437      if (idat.doParticlePot) {
438        // particlePot is the difference between the full potential and
# Line 424 | Line 464 | namespace OpenMD {
464  
465      RealType cut = 0.0;
466  
467 <    map<AtomType*, EAMAtomData>::iterator it;
468 <
469 <    it = EAMMap.find(atypes.first);
470 <    if (it != EAMMap.end()) {
471 <      EAMAtomData data1 = (*it).second;
467 >    int atid1 = atypes.first->getIdent();
468 >    int atid2 = atypes.second->getIdent();
469 >    int eamtid1 = EAMtids[atid1];
470 >    int eamtid2 = EAMtids[atid2];
471 >    
472 >    if (eamtid1 != -1) {
473 >      EAMAtomData data1 = EAMdata[eamtid1];
474        cut = data1.rcut;
475      }
476  
477 <    it = EAMMap.find(atypes.second);
478 <    if (it != EAMMap.end()) {
437 <      EAMAtomData data2 = (*it).second;
477 >    if (eamtid2 != -1) {
478 >      EAMAtomData data2 = EAMdata[eamtid2];
479        if (data2.rcut > cut)
480          cut = data2.rcut;
481      }
482 <
482 >    
483      return cut;
484    }
485   }

Diff Legend

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