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 1665 by gezelter, Tue Nov 22 20:38:56 2011 UTC vs.
Revision 1711 by gezelter, Sat May 19 02:58:35 2012 UTC

# Line 54 | Line 54 | namespace OpenMD {
54    SC::SC() : name_("SC"), initialized_(false), forceField_(NULL),
55               scRcut_(0.0), np_(3000) {}
56    
57 –  SCParam SC::getSCParam(AtomType* atomType) {
58 –    
59 –    // Do sanity checking on the AtomType we were passed before
60 –    // building any data structures:
61 –    if (!atomType->isSC()) {
62 –      sprintf( painCave.errMsg,
63 –               "SC::getSCParam was passed an atomType (%s) that does not\n"
64 –               "\tappear to be a Sutton-Chen (SC) atom.\n",
65 –               atomType->getName().c_str());
66 –      painCave.severity = OPENMD_ERROR;
67 –      painCave.isFatal = 1;
68 –      simError();
69 –    }
70 –    
71 –    GenericData* data = atomType->getPropertyByName("SC");
72 –    if (data == NULL) {
73 –      sprintf( painCave.errMsg, "SC::getSCParam could not find SC\n"
74 –               "\tparameters for atomType %s.\n",
75 –               atomType->getName().c_str());
76 –      painCave.severity = OPENMD_ERROR;
77 –      painCave.isFatal = 1;
78 –      simError();
79 –    }
80 –    
81 –    SCParamGenericData* scData = dynamic_cast<SCParamGenericData*>(data);
82 –    if (scData == NULL) {
83 –      sprintf( painCave.errMsg,
84 –               "SC::getSCParam could not convert GenericData to SCParamGenericData\n"
85 –               "\tfor atom type %s\n", atomType->getName().c_str());
86 –      painCave.severity = OPENMD_ERROR;
87 –      painCave.isFatal = 1;
88 –      simError();          
89 –    }
90 –    
91 –    return scData->getData();
92 –  }
93 –
94 –  RealType SC::getC(AtomType* atomType) {    
95 –    SCParam scParam = getSCParam(atomType);
96 –    return scParam.c;
97 –  }
98 –
99 –  RealType SC::getM(AtomType* atomType) {    
100 –    SCParam scParam = getSCParam(atomType);
101 –    return scParam.m;
102 –  }
103 –
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  
110 –  RealType SC::getN(AtomType* atomType) {    
111 –    SCParam scParam = getSCParam(atomType);
112 –    return scParam.n;
113 –  }
114 –
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  
121 –  RealType SC::getAlpha(AtomType* atomType) {    
122 –    SCParam scParam = getSCParam(atomType);
123 –    return scParam.alpha;
124 –  }
125 –
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 137 | 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 <
145 <  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
146 <    RealType epsilon1 = getEpsilon(atomType1);
147 <    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 156 | 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 166 | 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:
179 –    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 331 | Line 278 | namespace OpenMD {
278      *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
279  
280      (*(sdat.pot))[METALLIC_FAMILY] += u;
281 <    *(sdat.particlePot) += u;
282 <    
281 >    if (sdat.doParticlePot) {
282 >      *(sdat.particlePot) += u;
283 >    }
284 >
285      return;
286    }
287    
# Line 368 | Line 317 | namespace OpenMD {
317        
318        *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
319          
320 <      // particlePot is the difference between the full potential and
321 <      // the full potential without the presence of a particular
322 <      // particle (atom1).
323 <      //
324 <      // This reduces the density at other particle locations, so we
325 <      // need to recompute the density at atom2 assuming atom1 didn't
326 <      // contribute.  This then requires recomputing the density
327 <      // functional for atom2 as well.
328 <          
329 <      *(idat.particlePot1) -= data2.c * data2.epsilon *
330 <        sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
320 >      if (idat.doParticlePot) {
321 >        // particlePot is the difference between the full potential and
322 >        // the full potential without the presence of a particular
323 >        // particle (atom1).
324 >        //
325 >        // This reduces the density at other particle locations, so we
326 >        // need to recompute the density at atom2 assuming atom1 didn't
327 >        // contribute.  This then requires recomputing the density
328 >        // functional for atom2 as well.
329 >        
330 >        *(idat.particlePot1) -= data2.c * data2.epsilon *
331 >          sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
332  
333 <      *(idat.particlePot2) -= data1.c * data1.epsilon *
334 <        sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
333 >        *(idat.particlePot2) -= data1.c * data1.epsilon *
334 >          sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
335 >      }
336        
337        (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
338      }

Diff Legend

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