ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/nonbonded/EAM.cpp
(Generate patch)

Comparing:
branches/development/src/nonbonded/EAM.cpp (file contents), Revision 1478 by gezelter, Fri Jul 23 20:45:40 2010 UTC vs.
trunk/src/nonbonded/EAM.cpp (file contents), Revision 1927 by gezelter, Wed Aug 14 20:19:19 2013 UTC

# Line 35 | Line 35
35   *                                                                      
36   * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).            
37   * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).          
38 < * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008).          
39 < * [4]  Vardeman & Gezelter, in progress (2009).                        
38 > * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).          
39 > * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40 > * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41   */
42  
43   #include <stdio.h>
# Line 45 | Line 46
46   #include <cmath>
47   #include "nonbonded/EAM.hpp"
48   #include "utils/simError.h"
49 + #include "types/NonBondedInteractionType.hpp"
50  
51  
52   namespace OpenMD {
53  
54 <  bool EAM::initialized_ = false;
55 <  ForceField* EAM::forceField_ = NULL;
54 <  std::map<int, AtomType*> EAM::EAMlist;
55 <  std::map<AtomType*, EAMAtomData> EAM::EAMMap;
56 <  std::map<std::pair<AtomType*, AtomType*>, EAMInteractionData> EAM::MixingMap;
54 >  EAM::EAM() : name_("EAM"), initialized_(false), forceField_(NULL),
55 >               mixMeth_(eamJohnson), eamRcut_(0.0), haveCutoffRadius_(false) {}
56    
57 <  EAM* EAM::_instance = NULL;
57 >  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {  
58 >    EAMAdapter ea1 = EAMAdapter(atomType1);
59 >    EAMAdapter ea2 = EAMAdapter(atomType2);
60 >    CubicSpline* z1 = ea1.getZ();
61 >    CubicSpline* z2 = ea2.getZ();
62  
63 <  EAM* EAM::Instance() {
64 <    if (!_instance) {
65 <      _instance = new EAM();
66 <    }
67 <    return _instance;
65 <  }
66 <  
67 <  EAMParam EAM::getEAMParam(AtomType* atomType) {
68 <    
69 <    // Do sanity checking on the AtomType we were passed before
70 <    // building any data structures:
71 <    if (!atomType->isEAM()) {
72 <      sprintf( painCave.errMsg,
73 <               "EAM::getEAMParam was passed an atomType (%s) that does not\n"
74 <               "\tappear to be an embedded atom method (EAM) atom.\n",
75 <               atomType->getName().c_str());
76 <      painCave.severity = OPENMD_ERROR;
77 <      painCave.isFatal = 1;
78 <      simError();
79 <    }
80 <    
81 <    GenericData* data = atomType->getPropertyByName("EAM");
82 <    if (data == NULL) {
83 <      sprintf( painCave.errMsg, "EAM::getEAMParam could not find EAM\n"
84 <               "\tparameters for atomType %s.\n",
85 <               atomType->getName().c_str());
86 <      painCave.severity = OPENMD_ERROR;
87 <      painCave.isFatal = 1;
88 <      simError();
89 <    }
90 <    
91 <    EAMParamGenericData* eamData = dynamic_cast<EAMParamGenericData*>(data);
92 <    if (eamData == NULL) {
93 <      sprintf( painCave.errMsg,
94 <               "EAM::getEAMParam could not convert GenericData to EAMParam for\n"
95 <               "\tatom type %s\n", atomType->getName().c_str());
96 <      painCave.severity = OPENMD_ERROR;
97 <      painCave.isFatal = 1;
98 <      simError();          
99 <    }
100 <    
101 <    return eamData->getData();
102 <  }
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 <  CubicSpline* EAM::getZ(AtomType* atomType) {    
105 <    EAMParam eamParam = getEAMParam(atomType);
106 <    int nr = eamParam.nr;
107 <    RealType dr = eamParam.dr;
108 <    vector<RealType> rvals;
109 <    
110 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
111 <      
112 <    CubicSpline* cs = new CubicSpline();
113 <    cs->addPoints(rvals, eamParam.Z);
114 <    return cs;
115 <  }
69 >    // make the r grid:
70  
71 <  CubicSpline* EAM::getRho(AtomType* atomType) {    
118 <    EAMParam eamParam = getEAMParam(atomType);
119 <    int nr = eamParam.nr;
120 <    RealType dr = eamParam.dr;
121 <    vector<RealType> rvals;
71 >    // we need phi out to the largest value we'll encounter in the radial space;
72      
73 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
74 <      
75 <    CubicSpline* cs = new CubicSpline();
126 <    cs->addPoints(rvals, eamParam.rho);
127 <    return cs;
128 <  }
73 >    RealType rmax = 0.0;
74 >    rmax = max(rmax, ea1.getRcut());
75 >    rmax = max(rmax, ea1.getNr() * ea1.getDr());
76  
77 <  CubicSpline* EAM::getF(AtomType* atomType) {    
78 <    EAMParam eamParam = getEAMParam(atomType);
132 <    int nrho = eamParam.nrho;
133 <    RealType drho = eamParam.drho;
134 <    vector<RealType> rhovals;
135 <    vector<RealType> scaledF;
136 <    
137 <    for (int i = 0; i < nrho; i++) {
138 <      rhovals.push_back(i * drho);
139 <      scaledF.push_back( eamParam.F[i] * 23.06054 );
140 <    }
141 <      
142 <    CubicSpline* cs = new CubicSpline();
143 <    cs->addPoints(rhovals, eamParam.F);
144 <    return cs;
145 <  }
146 <  
147 <  CubicSpline* EAM::getPhi(AtomType* atomType1, AtomType* atomType2) {    
148 <    EAMParam eamParam1 = getEAMParam(atomType1);
149 <    EAMParam eamParam2 = getEAMParam(atomType2);
150 <    CubicSpline* z1 = getZ(atomType1);
151 <    CubicSpline* z2 = getZ(atomType2);
77 >    rmax = max(rmax, ea2.getRcut());
78 >    rmax = max(rmax, ea2.getNr() * ea2.getDr());
79  
80 <    // make the r grid:
80 >    // use the smallest dr (finest grid) to build our grid:
81  
82 <    // set rcut to be the smaller of the two atomic rcuts
82 >    RealType dr = min(ea1.getDr(), ea2.getDr());
83  
84 <    RealType rcut = eamParam1.rcut < eamParam2.rcut ?
158 <      eamParam1.rcut : eamParam2.rcut;
84 >    int nr = int(rmax/dr + 0.5);
85  
160 –    // use the smallest dr (finest grid) to build our grid:
161 –
162 –    RealType dr = eamParam1.dr < eamParam2.dr ? eamParam1.dr : eamParam2.dr;
163 –    int nr = int(rcut/dr);
86      vector<RealType> rvals;
87 <    for (int i = 0; i < nr; i++) rvals.push_back(i*dr);
87 >    for (int i = 0; i < nr; i++) rvals.push_back(RealType(i*dr));
88  
89      // construct the pair potential:
90  
# Line 173 | Line 95 | namespace OpenMD {
95  
96      phivals.push_back(0.0);
97  
98 <    for (int i = 1; i < rvals.size(); i++ ) {
98 >    for (unsigned int i = 1; i < rvals.size(); i++ ) {
99        r = rvals[i];
178 –      zi = z1->getValueAt(r);
179 –      zj = z2->getValueAt(r);
100  
101 <      phi = 331.999296 * (zi * zj) / r;
101 >      // only use z(r) if we're inside this atom's cutoff radius,
102 >      // otherwise, we'll use zero for the charge.  This effectively
103 >      // means that our phi grid goes out beyond the cutoff of the
104 >      // pair potential
105 >
106 >      zi = r <= ea1.getRcut() ? z1->getValueAt(r) : 0.0;
107 >      zj = r <= ea2.getRcut() ? z2->getValueAt(r) : 0.0;
108 >
109 >      phi = pre11_ * (zi * zj) / r;
110 >
111        phivals.push_back(phi);
112      }
113        
# Line 187 | Line 116 | namespace OpenMD {
116      return cs;
117    }
118  
119 <  void EAM::initialize() {
119 >  void EAM::setCutoffRadius( RealType rCut ) {
120 >    eamRcut_ = rCut;
121 >    haveCutoffRadius_ = true;
122 >  }
123  
124 +  void EAM::initialize() {
125      // set up the mixing method:
126 <    ForceFieldOptions ffo = forceField_->getForceFieldOptions();
127 <    string EAMMixMeth = toUpperCopy(ffo.getEAMMixingMethod());
128 <
126 >    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
127 >    string EAMMixMeth = fopts.getEAMMixingMethod();
128 >    toUpper(EAMMixMeth);
129 >  
130      if (EAMMixMeth == "JOHNSON")
131        mixMeth_ = eamJohnson;    
132      else if (EAMMixMeth == "DAW")
# Line 201 | 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 <      
211 <      if (at->isEAM())
212 <        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 222 | Line 164 | namespace OpenMD {
164        
165        if (nbt->isEAM()) {
166          
167 <        std::pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
167 >        pair<AtomType*, AtomType*> atypes = nbt->getAtomTypes();
168          
169          GenericData* data = nbt->getPropertyByName("EAM");
170          if (data == NULL) {
# Line 247 | Line 189 | namespace OpenMD {
189            simError();          
190          }
191          
192 <        EAMMix eamParam = eamData->getData();
192 >        EAMMixingParam eamParam = eamData->getData();
193  
194 <        vector<RealType> phiAB = eamParam.phiAB;
194 >        vector<RealType> phiAB = eamParam.phi;
195          RealType dr = eamParam.dr;
196          int nr = eamParam.nr;
197  
# Line 263 | Line 205 | namespace OpenMD {
205  
206    void EAM::addType(AtomType* atomType){
207  
208 +    EAMAdapter ea = EAMAdapter(atomType);
209      EAMAtomData eamAtomData;
210  
211 <    eamAtomData.rho = getRho(atomType);
212 <    eamAtomData.F = getF(atomType);
213 <    eamAtomData.Z = getZ(atomType);
214 <    eamAtomData.rcut = getRcut(atomType);
215 <
211 >    eamAtomData.rho = ea.getRho();
212 >    eamAtomData.F = ea.getF();
213 >    eamAtomData.Z = ea.getZ();
214 >    eamAtomData.rcut = ea.getRcut();
215 >    eamAtomData.isFluctuating = atomType->isFluctuatingCharge();
216 >      
217      // add it to the map:
218 <    AtomTypeProperties atp = atomType->getATP();    
218 >    int atid = atomType->getIdent();
219 >    int eamtid = EAMtypes.size();
220  
221 <    std::pair<std::map<int,AtomType*>::iterator,bool> ret;    
222 <    ret = EAMlist.insert( std::pair<int, AtomType*>(atp.ident, 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 <               atp.ident);
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 <    std::map<int, AtomType*>::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).second;
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 <      std::pair<AtomType*, AtomType*> key1, key2;
301 <      key1 = std::make_pair(atomType, atype2);
302 <      key2 = std::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 319 | Line 283 | namespace OpenMD {
283  
284      EAMInteractionData mixer;
285      CubicSpline* cs = new CubicSpline();
286 <    vector<RealType> rvals;
286 >    vector<RealType> rVals;
287  
288 <    for (int i = 0; i < nr; i++) rvals.push_back(i * dr);
288 >    for (int i = 0; i < nr; i++) rVals.push_back(i * dr);
289  
290      cs->addPoints(rVals, phiVals);
291      mixer.phi = cs;
292      mixer.explicitlySet = true;
293  
294 <    std::pair<AtomType*, AtomType*> key1, key2;
295 <    key1 = std::make_pair(atype1, atype2);
332 <    key2 = std::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    }
303  
304 <  void EAM::calcDensity(AtomType* at1, AtomType* at2, Vector3d d,
305 <                        RealType rij, RealType r2, RealType rho_i_at_j,
343 <                        RealType rho_j_at_i) {
344 <
304 >  void EAM::calcDensity(InteractionData &idat) {
305 >    
306      if (!initialized_) initialize();
307 +    
308 +    EAMAtomData &data1 = EAMdata[EAMtids[idat.atid1]];
309 +    EAMAtomData &data2 = EAMdata[EAMtids[idat.atid2]];
310  
311 <    EAMAtomData data1 = EAMMap[at1];
312 <    EAMAtomData data2 = EAMMap[at2];
313 <
314 <    if (rij < data1.rcut) rho_i_at_j = data1.rho->getValueAt(rij);
315 <    if (rij < data2.rcut) rho_j_at_i = data2.rho->getValueAt(rij);
316 <    return;
311 >    if (haveCutoffRadius_)
312 >      if ( *(idat.rij) > eamRcut_) return;
313 >    
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 >      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    }
334 <
335 <  void EAM::calcFunctional(AtomType* at1, RealType rho, RealType frho,
336 <                           RealType dfrhodrho) {
357 <
334 >  
335 >  void EAM::calcFunctional(SelfData &sdat) {
336 >    
337      if (!initialized_) initialize();
338  
339 <    EAMAtomData data1 = EAMMap[at1];
340 <        
341 <    pair<RealType, RealType> result = data1.F->getValueAndDerivativeAt(rho);
339 >    EAMAtomData &data1 = EAMdata[ EAMtids[sdat.atid] ];
340 >            
341 >    data1.F->getValueAndDerivativeAt( *(sdat.rho), *(sdat.frho), *(sdat.dfrhodrho) );
342  
343 <    frho = result.first;
344 <    dfrhodrho = result.second;
343 >    (*(sdat.pot))[METALLIC_FAMILY] += *(sdat.frho);
344 >    if (sdat.doParticlePot) {
345 >      *(sdat.particlePot) += *(sdat.frho);
346 >    }
347 >
348      return;
349    }
350  
351  
352 <  void EAM::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
371 <                      RealType rij, RealType r2, RealType sw,
372 <                      RealType &vpair, RealType &pot, Vector3d &f1,
373 <                      RealType rho1, RealType rho2, RealType dfrho1,
374 <                      RealType dfrho2, RealType fshift1, RealType fshift2) {
352 >  void EAM::calcForce(InteractionData &idat) {
353  
354      if (!initialized_) initialize();
377 –    
378 –    pair<RealType, RealType> res;
379 –    
380 –    if (rij < eamRcut_) {
355  
356 <      EAMAtomData data1 = EAMMap[at1];
357 <      EAMAtomData data2 = EAMMap[at2];
356 >    if (haveCutoffRadius_)
357 >      if ( *(idat.rij) > eamRcut_) return;
358 >  
359  
360 <      // get type-specific cutoff radii
361 <
362 <      RealType rci = data1.rcut;
363 <      RealType rcj = data2.rcut;
364 <      
365 <      RealType rha, drha, rhb, drhb;
366 <      RealType pha, dpha, phb, dphb;
367 <      RealType phab, dvpdr;
368 <      RealType drhoidr, drhojdr, dudr;
369 <      
370 <      if (rij < rci) {
371 <        res = data1.rho->getValueAndDerivativeAt(rij);
372 <        rha = res.first;
373 <        drha = res.second;
374 <
375 <        res = MixingMap[make_pair(at1, at1)].phi->getValueAndDerivativeAt(rij);
376 <        pha = res.first;
377 <        dpha = res.second;
360 >    int eamtid1 = EAMtids[idat.atid1];
361 >    int eamtid2 = EAMtids[idat.atid2];
362 >    
363 >    EAMAtomData &data1 = EAMdata[eamtid1];
364 >    EAMAtomData &data2 = EAMdata[eamtid2];
365 >    
366 >    // get type-specific cutoff radii
367 >    
368 >    RealType rci = data1.rcut;
369 >    RealType rcj = data2.rcut;
370 >    
371 >    RealType rha(0.0), drha(0.0), rhb(0.0), drhb(0.0);
372 >    RealType pha(0.0), dpha(0.0), phb(0.0), dphb(0.0);
373 >    RealType phab(0.0), dvpdr(0.0);
374 >    RealType drhoidr, drhojdr, dudr;
375 >    
376 >    if ( *(idat.rij) < rci) {
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 <      if (rij < rcj) {
385 <        res = data2.rho->getValueAndDerivativeAt(rij);
386 <        rhb = res.first;
387 <        drhb = res.second;
388 <
389 <        res = MixingMap[make_pair(at2, at2)].phi->getValueAndDerivativeAt(rij);
390 <        phb = res.first;
412 <        dphb = res.second;
383 >    }
384 >    
385 >    if ( *(idat.rij) < rcj) {
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 <      phab = 0.0;
395 <      dvpdr = 0.0;
396 <
397 <      switch(mixMeth_) {
398 <      case eamJohnson:
399 <      
400 <        if (rij < rci) {
422 <          phab = phab + 0.5 * (rhb / rha) * pha;
423 <          dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
424 <                               pha*((drhb/rha) - (rhb*drha/rha/rha)));
425 <        }
426 <
427 <        if (rij < rcj) {
428 <          phab = phab + 0.5 * (rha / rhb) * phb;
429 <          dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
430 <                                 phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
431 <        }
432 <
433 <        break;
434 <
435 <      case eamDaw:
436 <                
437 <        res = MixingMap[make_pair(at1,at2)].phi->getValueAndDerivativeAt(rij);
438 <        phab = res.first;
439 <        dvpdr = res.second;
440 <
441 <        break;
442 <      case eamUnknown:
443 <      default:
444 <
445 <        sprintf(painCave.errMsg,
446 <                "EAM::calcForce hit a mixing method it doesn't know about!\n"
447 <                );
448 <        painCave.severity = OPENMD_ERROR;
449 <        painCave.isFatal = 1;
450 <        simError();        
451 <          
394 >    switch(mixMeth_) {
395 >    case eamJohnson:
396 >      
397 >      if ( *(idat.rij) < rci) {
398 >        phab = phab + 0.5 * (rhb / rha) * pha;
399 >        dvpdr = dvpdr + 0.5*((rhb/rha)*dpha +
400 >                             pha*((drhb/rha) - (rhb*drha/rha/rha)));
401        }
402        
403 <      drhoidr = drha;
404 <      drhojdr = drhb;
403 >      
404 >      
405 >      if ( *(idat.rij) < rcj) {
406 >        phab = phab + 0.5 * (rha / rhb) * phb;
407 >        dvpdr = dvpdr + 0.5 * ((rha/rhb)*dphb +
408 >                               phb*((drha/rhb) - (rha*drhb/rhb/rhb)));
409 >      }
410 >      
411 >      break;
412 >      
413 >    case eamDaw:
414 >      MixingMap[eamtid1][eamtid2].phi->getValueAndDerivativeAt( *(idat.rij), phab, dvpdr);
415 >      
416 >      break;
417 >    case eamUnknown:
418 >    default:
419 >      
420 >      sprintf(painCave.errMsg,
421 >              "EAM::calcForce hit a mixing method it doesn't know about!\n"
422 >              );
423 >      painCave.severity = OPENMD_ERROR;
424 >      painCave.isFatal = 1;
425 >      simError();        
426 >      
427 >    }
428 >    
429 >    drhoidr = drha;
430 >    drhojdr = drhb;
431 >    
432 >    dudr = drhojdr* *(idat.dfrho1) + drhoidr* *(idat.dfrho2) + dvpdr;
433 >    
434 >    *(idat.f1) += *(idat.d) * dudr / *(idat.rij);
435  
457 –      dudr = drhojdr*dfrhodrho_i + drhoidr*dfrhodrho_j + dvpdr;
458 –
459 –      f1 = d * dudr / rij;
436          
437 <      // particle_pot is the difference between the full potential
438 <      // and the full potential without the presence of a particular
437 >    if (idat.doParticlePot) {
438 >      // particlePot is the difference between the full potential and
439 >      // the full potential without the presence of a particular
440        // particle (atom1).
441        //
442 <      // This reduces the density at other particle locations, so
443 <      // we need to recompute the density at atom2 assuming atom1
444 <      // didn't contribute.  This then requires recomputing the
445 <      // density functional for atom2 as well.
446 <      //
447 <      // Most of the particle_pot heavy lifting comes from the
448 <      // pair interaction, and will be handled by vpair.
449 <    
450 <      fshift_i = data1.F->getValueAt( rho_i - rhb );
451 <      fshift_j = data1.F->getValueAt( rho_j - rha );
475 <
476 <      pot += phab;
477 <
478 <      vpair += phab;
442 >      // This reduces the density at other particle locations, so we
443 >      // need to recompute the density at atom2 assuming atom1 didn't
444 >      // contribute.  This then requires recomputing the density
445 >      // functional for atom2 as well.
446 >      
447 >      *(idat.particlePot1) += data2.F->getValueAt( *(idat.rho2) - rha )
448 >        - *(idat.frho2);
449 >      
450 >      *(idat.particlePot2) += data1.F->getValueAt( *(idat.rho1) - rhb)
451 >        - *(idat.frho1);
452      }
453 <
453 >    
454 >    (*(idat.pot))[METALLIC_FAMILY] += phab;
455 >    
456 >    *(idat.vpair) += phab;
457 >  
458      return;
459      
460    }
461  
462 +  RealType EAM::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
463 +    if (!initialized_) initialize();  
464  
465 <  void EAM::calc_eam_prepair_rho(int *atid1, int *atid2, RealType *d,
487 <                                 RealType *rij, RealType *r2,
488 <                                 RealType* rho_i_at_j, RealType* rho_j_at_i){
489 <    if (!initialized_) initialize();
465 >    RealType cut = 0.0;
466  
467 <    AtomType* atype1 = EAMlist[*atid1];
468 <    AtomType* atype2 = EAMlist[*atid2];
467 >    int atid1 = atypes.first->getIdent();
468 >    int atid2 = atypes.second->getIdent();
469 >    int eamtid1 = EAMtids[atid1];
470 >    int eamtid2 = EAMtids[atid2];
471      
472 <    Vector3d disp(d[0], d[1], d[2]);
472 >    if (eamtid1 != -1) {
473 >      EAMAtomData data1 = EAMdata[eamtid1];
474 >      cut = data1.rcut;
475 >    }
476  
477 <    calcDensity(atype1, atype2, disp, *rij, *r2, *rho_i_at_j, *rho_j_at_i);
478 <
479 <    return;    
480 <  }
481 <
501 <  void EAM::calc_eam_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
502 <                                   RealType *dfrhodrho) {
503 <
504 <    if (!initialized_) initialize();
505 <
506 <    AtomType* atype1 = EAMlist[*atid1];  
507 <
508 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
477 >    if (eamtid2 != -1) {
478 >      EAMAtomData data2 = EAMdata[eamtid2];
479 >      if (data2.rcut > cut)
480 >        cut = data2.rcut;
481 >    }
482      
483 <    return;    
483 >    return cut;
484    }
512 –
513 –  void EAM::do_eam_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
514 –                        RealType *r2, RealType *sw, RealType *vpair,
515 –                        RealType *pot, RealType *f1, RealType *rho1,
516 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
517 –                        RealType *fshift1, RealType *fshift2) {
518 –
519 –    if (!initialized_) initialize();
520 –    
521 –    AtomType* atype1 = EAMMap[*atid1];
522 –    AtomType* atype2 = EAMMap[*atid2];
523 –    
524 –    Vector3d disp(d[0], d[1], d[2]);
525 –    Vector3d frc(f1[0], f1[1], f1[2]);
526 –    
527 –    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
528 –              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
529 –      
530 –    f1[0] = frc.x();
531 –    f1[1] = frc.y();
532 –    f1[2] = frc.z();
533 –
534 –    return;    
535 –  }
536 –  
537 –  void EAM::setCutoffEAM(RealType *thisRcut) {
538 –    eamRcut_ = thisRcut;
539 –  }
485   }
486  
542 – extern "C" {
543 –  
544 – #define fortranCalcDensity FC_FUNC(calc_eam_prepair_rho, CALC_EAM_PREPAIR_RHO)
545 – #define fortranCalcFunctional FC_FUNC(calc_eam_preforce_frho, CALC_EAM_PREFORCE_FRHO)
546 – #define fortranCalcForce FC_FUNC(do_eam_pair, DO_EAM_PAIR)
547 – #define fortranSetCutoffEAM FC_FUNC(setcutoffeam, SETCUTOFFEAM)
548 –  
549 –  RealType fortranCalcDensity(int *atid1, int *atid2, RealType *d,
550 –                              RealType *rij, RealType *r2,
551 –                              RealType *rho_i_at_j, RealType *rho_j_at_i) {
552 –
553 –    return OpenMD::EAM::Instance()->calc_eam_prepair_rho(*atid1, *atid2, *d,
554 –                                                         *rij, *r2,
555 –                                                         *rho_i_at_j,  
556 –                                                         *rho_j_at_i);
557 –  }
558 –  RealType fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
559 –                                 RealType *dfrhodrho) {  
560 –
561 –    return OpenMD::EAM::Instance()->calc_eam_preforce_Frho(*atid1,
562 –                                                           *rho,
563 –                                                           *frho,
564 –                                                           *dfrhodrho);
565 –
566 –  }
567 –  void fortranSetEAMCutoff(RealType *rcut) {
568 –    return OpenMD::EAM::Instance()->setCutoffEAM(rcut);
569 –  }
570 –  void fortranDoEAMPair(int *atid1, int *atid2, RealType *d, RealType *rij,
571 –                        RealType *r2, RealType *sw, RealType *vpair,
572 –                        RealType *pot, RealType *f1, RealType *rho1,
573 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
574 –                        RealType *fshift1, RealType *fshift2){
575 –    
576 –    return OpenMD::EAM::Instance()->do_eam_pair(*atid1, *atid2, *d, *rij,
577 –                                                *r2, *sw,  *vpair,
578 –                                                *pot, *f1,  *rho1,
579 –                                                *rho2,  *dfrho1,  *dfrho2,
580 –                                                *fshift1,  *fshift2);
581 –  }
582 – }

Diff Legend

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