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

Comparing branches/development/src/nonbonded/SC.cpp (file contents):
Revision 1502 by gezelter, Sat Oct 2 19:53:32 2010 UTC vs.
Revision 1710 by gezelter, Fri May 18 21:44:02 2012 UTC

# Line 36 | Line 36
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).                        
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 53 | Line 54 | namespace OpenMD {
54    SC::SC() : name_("SC"), initialized_(false), forceField_(NULL),
55               scRcut_(0.0), np_(3000) {}
56    
56 –  SCParam SC::getSCParam(AtomType* atomType) {
57 –    
58 –    // Do sanity checking on the AtomType we were passed before
59 –    // building any data structures:
60 –    if (!atomType->isSC()) {
61 –      sprintf( painCave.errMsg,
62 –               "SC::getSCParam was passed an atomType (%s) that does not\n"
63 –               "\tappear to be a Sutton-Chen (SC) atom.\n",
64 –               atomType->getName().c_str());
65 –      painCave.severity = OPENMD_ERROR;
66 –      painCave.isFatal = 1;
67 –      simError();
68 –    }
69 –    
70 –    GenericData* data = atomType->getPropertyByName("SC");
71 –    if (data == NULL) {
72 –      sprintf( painCave.errMsg, "SC::getSCParam could not find SC\n"
73 –               "\tparameters for atomType %s.\n",
74 –               atomType->getName().c_str());
75 –      painCave.severity = OPENMD_ERROR;
76 –      painCave.isFatal = 1;
77 –      simError();
78 –    }
79 –    
80 –    SCParamGenericData* scData = dynamic_cast<SCParamGenericData*>(data);
81 –    if (scData == NULL) {
82 –      sprintf( painCave.errMsg,
83 –               "SC::getSCParam could not convert GenericData to SCParamGenericData\n"
84 –               "\tfor atom type %s\n", atomType->getName().c_str());
85 –      painCave.severity = OPENMD_ERROR;
86 –      painCave.isFatal = 1;
87 –      simError();          
88 –    }
89 –    
90 –    return scData->getData();
91 –  }
92 –
93 –  RealType SC::getC(AtomType* atomType) {    
94 –    SCParam scParam = getSCParam(atomType);
95 –    return scParam.c;
96 –  }
97 –
98 –  RealType SC::getM(AtomType* atomType) {    
99 –    SCParam scParam = getSCParam(atomType);
100 –    return scParam.m;
101 –  }
102 –
57    RealType SC::getM(AtomType* atomType1, AtomType* atomType2) {    
58 <    RealType m1 = getM(atomType1);
59 <    RealType m2 = getM(atomType2);
58 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
59 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
60 >    RealType m1 = sca1.getM();
61 >    RealType m2 = sca2.getM();
62      return 0.5 * (m1 + m2);
63    }
64  
109 –  RealType SC::getN(AtomType* atomType) {    
110 –    SCParam scParam = getSCParam(atomType);
111 –    return scParam.n;
112 –  }
113 –
65    RealType SC::getN(AtomType* atomType1, AtomType* atomType2) {    
66 <    RealType n1 = getN(atomType1);
67 <    RealType n2 = getN(atomType2);
66 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
67 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
68 >    RealType n1 = sca1.getN();
69 >    RealType n2 = sca2.getN();
70      return 0.5 * (n1 + n2);
71    }
72  
120 –  RealType SC::getAlpha(AtomType* atomType) {    
121 –    SCParam scParam = getSCParam(atomType);
122 –    return scParam.alpha;
123 –  }
124 –
73    RealType SC::getAlpha(AtomType* atomType1, AtomType* atomType2) {    
74 <    RealType alpha1 = getAlpha(atomType1);
75 <    RealType alpha2 = getAlpha(atomType2);
74 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
75 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
76 >    RealType alpha1 = sca1.getAlpha();
77 >    RealType alpha2 = sca2.getAlpha();
78  
79      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
80      std::string DistanceMix = fopts.getDistanceMixingRule();
# Line 136 | Line 86 | namespace OpenMD {
86        return 0.5 * (alpha1 + alpha2);
87    }
88  
89 <  RealType SC::getEpsilon(AtomType* atomType) {    
90 <    SCParam scParam = getSCParam(atomType);
91 <    return scParam.epsilon;
92 <  }
93 <
144 <  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
145 <    RealType epsilon1 = getEpsilon(atomType1);
146 <    RealType epsilon2 = getEpsilon(atomType2);
89 >  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {  
90 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
91 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
92 >    RealType epsilon1 = sca1.getEpsilon();
93 >    RealType epsilon2 = sca2.getEpsilon();
94      return sqrt(epsilon1 * epsilon2);
95    }
96  
# Line 155 | Line 102 | namespace OpenMD {
102  
103      for (at = atomTypes->beginType(i); at != NULL;
104           at = atomTypes->nextType(i)) {
105 <      if (at->isSC())
105 >      SuttonChenAdapter sca = SuttonChenAdapter(at);
106 >      if (sca.isSuttonChen())
107          addType(at);
108      }    
109      initialized_ = true;
# Line 165 | Line 113 | namespace OpenMD {
113  
114    void SC::addType(AtomType* atomType){
115  
116 +    SuttonChenAdapter sca = SuttonChenAdapter(atomType);
117      SCAtomData scAtomData;
118      
119 <    scAtomData.c = getC(atomType);
120 <    scAtomData.m = getM(atomType);
121 <    scAtomData.n = getN(atomType);
122 <    scAtomData.alpha = getAlpha(atomType);
123 <    scAtomData.epsilon = getEpsilon(atomType);
119 >    scAtomData.c = sca.getC();
120 >    scAtomData.m = sca.getM();
121 >    scAtomData.n = sca.getN();
122 >    scAtomData.alpha = sca.getAlpha();
123 >    scAtomData.epsilon = sca.getEpsilon();
124      scAtomData.rCut = 2.0 * scAtomData.alpha;
125  
126      // add it to the map:
178 –    AtomTypeProperties atp = atomType->getATP();    
127  
128      pair<map<int,AtomType*>::iterator,bool> ret;    
129 <    ret = SClist.insert( pair<int, AtomType*>(atp.ident, atomType) );
129 >    ret = SClist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
130      if (ret.second == false) {
131        sprintf( painCave.errMsg,
132                 "SC already had a previous entry with ident %d\n",
133 <               atp.ident);
133 >               atomType->getIdent() );
134        painCave.severity = OPENMD_INFO;
135        painCave.isFatal = 0;
136        simError();        
# Line 302 | Line 250 | namespace OpenMD {
250      return;
251    }
252  
253 <  void SC::calcDensity(DensityData ddat) {
253 >  void SC::calcDensity(InteractionData &idat) {
254      
255      if (!initialized_) initialize();
256      
257 <    SCInteractionData mixer = MixingMap[make_pair(ddat.atype1, ddat.atype2)];
257 >    SCInteractionData mixer = MixingMap[ idat.atypes ];
258  
259      RealType rcij = mixer.rCut;
260  
261 <    if (ddat.rij < rcij) {
262 <      ddat.rho_i_at_j = mixer.phi->getValueAt(ddat.rij);
263 <      ddat.rho_j_at_i = ddat.rho_i_at_j;
264 <    } else {
265 <      ddat.rho_i_at_j = 0.0;
266 <      ddat.rho_j_at_i = 0.0;
319 <    }
320 <
261 >    if ( *(idat.rij)  < rcij) {
262 >      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
263 >      *(idat.rho1) += rho;
264 >      *(idat.rho2) += rho;
265 >    }
266 >    
267      return;
268    }
269  
270 <  void SC::calcFunctional(FunctionalData fdat) {
270 >  void SC::calcFunctional(SelfData &sdat) {
271  
272      if (!initialized_) initialize();
273  
274 <    SCAtomData data1 = SCMap[fdat.atype];
274 >    SCAtomData data1 = SCMap[sdat.atype];
275 >  
276 >    RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
277 >    *(sdat.frho) = u;
278 >    *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
279 >
280 >    (*(sdat.pot))[METALLIC_FAMILY] += u;
281 >    *(sdat.particlePot) += u;
282      
330 –    fdat.frho = - data1.c * data1.epsilon * sqrt(fdat.rho);
331 –    fdat.dfrhodrho = 0.5 * fdat.frho / fdat.rho;
332 –    
283      return;
284    }
285    
286  
287 <  void SC::calcForce(InteractionData idat) {
287 >  void SC::calcForce(InteractionData &idat) {
288      
289      if (!initialized_) initialize();
290      
291 <    SCAtomData data1 = SCMap[idat.atype1];
292 <    SCAtomData data2 = SCMap[idat.atype2];
291 >    SCAtomData data1 = SCMap[idat.atypes.first];
292 >    SCAtomData data2 = SCMap[idat.atypes.second];
293  
294 <    SCInteractionData mixer = MixingMap[make_pair(idat.atype1, idat.atype2)];
294 >    SCInteractionData mixer = MixingMap[idat.atypes];
295  
296      RealType rcij = mixer.rCut;
297  
298 <    if (idat.rij < rcij) {
298 >    if ( *(idat.rij)  < rcij) {
299        RealType vcij = mixer.vCut;
300        
301        pair<RealType, RealType> res;
302        
303 <      res = mixer.phi->getValueAndDerivativeAt(idat.rij);
303 >      res = mixer.phi->getValueAndDerivativeAt( *(idat.rij) );
304        RealType rhtmp = res.first;
305        RealType drhodr = res.second;
306        
307 <      res = mixer.V->getValueAndDerivativeAt(idat.rij);
307 >      res = mixer.V->getValueAndDerivativeAt( *(idat.rij) );
308        RealType vptmp = res.first;
309        RealType dvpdr = res.second;
310        
311        RealType pot_temp = vptmp - vcij;
312 <      idat.vpair += pot_temp;
312 >      *(idat.vpair) += pot_temp;
313        
314 <      RealType dudr = drhodr * (idat.dfrho1 + idat.dfrho2) + dvpdr;
314 >      RealType dudr = drhodr * ( *(idat.dfrho1) + *(idat.dfrho2) ) + dvpdr;
315        
316 <      idat.f1 += idat.d * dudr / idat.rij;
316 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
317          
318 <      // particle_pot is the difference between the full potential
319 <      // and the full potential without the presence of a particular
318 >      // particlePot is the difference between the full potential and
319 >      // the full potential without the presence of a particular
320        // particle (atom1).
321        //
322 <      // This reduces the density at other particle locations, so
323 <      // we need to recompute the density at atom2 assuming atom1
324 <      // didn't contribute.  This then requires recomputing the
325 <      // density functional for atom2 as well.
326 <      //
327 <      // Most of the particle_pot heavy lifting comes from the
328 <      // pair interaction, and will be handled by vpair.
322 >      // This reduces the density at other particle locations, so we
323 >      // need to recompute the density at atom2 assuming atom1 didn't
324 >      // contribute.  This then requires recomputing the density
325 >      // functional for atom2 as well.
326 >          
327 >      *(idat.particlePot1) -= data2.c * data2.epsilon *
328 >        sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
329 >
330 >      *(idat.particlePot2) -= data1.c * data1.epsilon *
331 >        sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
332        
333 <      idat.fshift1 = - data1.c * data1.epsilon * sqrt(idat.rho1 - rhtmp);
381 <      idat.fshift2 = - data2.c * data2.epsilon * sqrt(idat.rho2 - rhtmp);
382 <      
383 <      idat.pot += pot_temp;
333 >      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
334      }
335        
336      return;    
337    }
338 +
339 +  RealType SC::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
340 +    if (!initialized_) initialize();  
341 +
342 +    map<pair<AtomType*, AtomType*>, SCInteractionData>::iterator it;
343 +    it = MixingMap.find(atypes);
344 +    if (it == MixingMap.end())
345 +      return 0.0;
346 +    else  {
347 +      SCInteractionData mixer = (*it).second;
348 +      return mixer.rCut;
349 +    }
350 +  }
351   }

Diff Legend

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