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

Comparing:
branches/development/src/nonbonded/SC.cpp (file contents), Revision 1489 by gezelter, Tue Aug 10 18:34:59 2010 UTC vs.
trunk/src/nonbonded/SC.cpp (file contents), Revision 1895 by gezelter, Mon Jul 1 21:09:37 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 49 | Line 50 | namespace OpenMD {
50  
51   namespace OpenMD {
52  
52 –  bool SC::initialized_ = false;
53 –  RealType SC::scRcut_ = 0.0;
54 –  int SC::np_ = 3000;
55 –  ForceField* SC::forceField_ = NULL;
56 –  map<int, AtomType*> SC::SClist;
57 –  map<AtomType*, SCAtomData> SC::SCMap;
58 –  map<pair<AtomType*, AtomType*>, SCInteractionData> SC::MixingMap;
59 –  
60 –  SC* SC::_instance = NULL;
53  
54 <  SC* SC::Instance() {
55 <    if (!_instance) {
64 <      _instance = new SC();
65 <    }
66 <    return _instance;
67 <  }
54 >  SC::SC() : name_("SC"), initialized_(false), forceField_(NULL),
55 >             scRcut_(0.0), np_(3000) {}
56    
57 <  SCParam SC::getSCParam(AtomType* atomType) {
58 <    
71 <    // Do sanity checking on the AtomType we were passed before
72 <    // building any data structures:
73 <    if (!atomType->isSC()) {
74 <      sprintf( painCave.errMsg,
75 <               "SC::getSCParam was passed an atomType (%s) that does not\n"
76 <               "\tappear to be a Sutton-Chen (SC) atom.\n",
77 <               atomType->getName().c_str());
78 <      painCave.severity = OPENMD_ERROR;
79 <      painCave.isFatal = 1;
80 <      simError();
81 <    }
82 <    
83 <    GenericData* data = atomType->getPropertyByName("SC");
84 <    if (data == NULL) {
85 <      sprintf( painCave.errMsg, "SC::getSCParam could not find SC\n"
86 <               "\tparameters for atomType %s.\n",
87 <               atomType->getName().c_str());
88 <      painCave.severity = OPENMD_ERROR;
89 <      painCave.isFatal = 1;
90 <      simError();
91 <    }
92 <    
93 <    SCParamGenericData* scData = dynamic_cast<SCParamGenericData*>(data);
94 <    if (scData == NULL) {
95 <      sprintf( painCave.errMsg,
96 <               "SC::getSCParam could not convert GenericData to SCParamGenericData\n"
97 <               "\tfor atom type %s\n", atomType->getName().c_str());
98 <      painCave.severity = OPENMD_ERROR;
99 <      painCave.isFatal = 1;
100 <      simError();          
101 <    }
102 <    
103 <    return scData->getData();
104 <  }
57 >  SC::~SC() {
58 >    initialized_ = false;
59  
60 <  RealType SC::getC(AtomType* atomType) {    
61 <    SCParam scParam = getSCParam(atomType);
62 <    return scParam.c;
60 >    MixingMap.clear();
61 >    SCtypes.clear();
62 >    SCdata.clear();
63 >    SCtids.clear();
64    }
65 <
111 <  RealType SC::getM(AtomType* atomType) {    
112 <    SCParam scParam = getSCParam(atomType);
113 <    return scParam.m;
114 <  }
115 <
65 >        
66    RealType SC::getM(AtomType* atomType1, AtomType* atomType2) {    
67 <    RealType m1 = getM(atomType1);
68 <    RealType m2 = getM(atomType2);
67 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
68 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
69 >    RealType m1 = sca1.getM();
70 >    RealType m2 = sca2.getM();
71      return 0.5 * (m1 + m2);
72    }
73  
122 –  RealType SC::getN(AtomType* atomType) {    
123 –    SCParam scParam = getSCParam(atomType);
124 –    return scParam.n;
125 –  }
126 –
74    RealType SC::getN(AtomType* atomType1, AtomType* atomType2) {    
75 <    RealType n1 = getN(atomType1);
76 <    RealType n2 = getN(atomType2);
75 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
76 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
77 >    RealType n1 = sca1.getN();
78 >    RealType n2 = sca2.getN();
79      return 0.5 * (n1 + n2);
80    }
81  
133 –  RealType SC::getAlpha(AtomType* atomType) {    
134 –    SCParam scParam = getSCParam(atomType);
135 –    return scParam.alpha;
136 –  }
137 –
82    RealType SC::getAlpha(AtomType* atomType1, AtomType* atomType2) {    
83 <    RealType alpha1 = getAlpha(atomType1);
84 <    RealType alpha2 = getAlpha(atomType2);
83 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
84 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
85 >    RealType alpha1 = sca1.getAlpha();
86 >    RealType alpha2 = sca2.getAlpha();
87  
88      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
89      std::string DistanceMix = fopts.getDistanceMixingRule();
# Line 149 | Line 95 | namespace OpenMD {
95        return 0.5 * (alpha1 + alpha2);
96    }
97  
98 <  RealType SC::getEpsilon(AtomType* atomType) {    
99 <    SCParam scParam = getSCParam(atomType);
100 <    return scParam.epsilon;
101 <  }
102 <
157 <  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
158 <    RealType epsilon1 = getEpsilon(atomType1);
159 <    RealType epsilon2 = getEpsilon(atomType2);
98 >  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {  
99 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
100 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
101 >    RealType epsilon1 = sca1.getEpsilon();
102 >    RealType epsilon2 = sca2.getEpsilon();
103      return sqrt(epsilon1 * epsilon2);
104    }
105  
106 <  void SC::initialize() {
107 <    // find all of the SC atom Types:
108 <    ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
109 <    ForceField::AtomTypeContainer::MapTypeIterator i;
110 <    AtomType* at;
106 >  void SC::initialize() {      
107 >    // find all of the SC atom Types:
108 >    SCtypes.clear();
109 >    SCtids.clear();
110 >    SCdata.clear();
111 >    MixingMap.clear();
112 >    nSC_ = 0;
113  
114 <    for (at = atomTypes->beginType(i); at != NULL;
115 <         at = atomTypes->nextType(i)) {
116 <      if (at->isSC())
117 <        addType(at);
118 <    }    
114 >    SCtids.resize( forceField_->getNAtomType(), -1);
115 >
116 >    set<AtomType*>::iterator at;
117 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
118 >      if ((*at)->isSC()) nSC_++;
119 >    }
120 >    SCdata.resize(nSC_);
121 >    MixingMap.resize(nSC_);
122 >    for (at = simTypes_.begin(); at != simTypes_.end(); ++at) {
123 >      if ((*at)->isSC()) addType((*at));
124 >    }
125      initialized_ = true;
126    }
127    
# Line 178 | Line 129 | namespace OpenMD {
129  
130    void SC::addType(AtomType* atomType){
131  
132 +    SuttonChenAdapter sca = SuttonChenAdapter(atomType);
133      SCAtomData scAtomData;
134      
135 <    scAtomData.c = getC(atomType);
136 <    scAtomData.m = getM(atomType);
137 <    scAtomData.n = getN(atomType);
138 <    scAtomData.alpha = getAlpha(atomType);
139 <    scAtomData.epsilon = getEpsilon(atomType);
135 >    scAtomData.c = sca.getC();
136 >    scAtomData.m = sca.getM();
137 >    scAtomData.n = sca.getN();
138 >    scAtomData.alpha = sca.getAlpha();
139 >    scAtomData.epsilon = sca.getEpsilon();
140      scAtomData.rCut = 2.0 * scAtomData.alpha;
141 <
141 >
142      // add it to the map:
143 <    AtomTypeProperties atp = atomType->getATP();    
143 >    int atid = atomType->getIdent();
144 >    int sctid = SCtypes.size();
145  
146 <    pair<map<int,AtomType*>::iterator,bool> ret;    
147 <    ret = SClist.insert( pair<int, AtomType*>(atp.ident, atomType) );
146 >    pair<set<int>::iterator,bool> ret;    
147 >    ret = SCtypes.insert( atid );
148      if (ret.second == false) {
149        sprintf( painCave.errMsg,
150                 "SC already had a previous entry with ident %d\n",
151 <               atp.ident);
151 >               atid );
152        painCave.severity = OPENMD_INFO;
153        painCave.isFatal = 0;
154        simError();        
155      }
203 –
204 –    SCMap[atomType] = scAtomData;
156      
157 +    SCtids[atid] = sctid;
158 +    SCdata[sctid] = scAtomData;
159 +    MixingMap[sctid].resize(nSC_);
160 +    
161      // Now, iterate over all known types and add to the mixing map:
162      
163 <    map<AtomType*, SCAtomData>::iterator it;
164 <    for( it = SCMap.begin(); it != SCMap.end(); ++it) {
163 >    std::set<int>::iterator it;
164 >    for( it = SCtypes.begin(); it != SCtypes.end(); ++it) {
165        
166 <      AtomType* atype2 = (*it).first;
166 >      int sctid2 = SCtids[ (*it) ];
167 >      AtomType* atype2 = forceField_->getAtomType( (*it) );
168        
169        SCInteractionData mixer;
170  
# Line 247 | Line 203 | namespace OpenMD {
203  
204        mixer.explicitlySet = false;
205  
206 <      pair<AtomType*, AtomType*> key1, key2;
251 <      key1 = make_pair(atomType, atype2);
252 <      key2 = make_pair(atype2, atomType);
206 >      MixingMap[sctid2].resize( nSC_ );
207        
208 <      MixingMap[key1] = mixer;
209 <      if (key2 != key1) {
210 <        MixingMap[key2] = mixer;
208 >      MixingMap[sctid][sctid2] = mixer;
209 >      if (sctid2 != sctid) {
210 >        MixingMap[sctid2][sctid] = mixer;
211        }
212      }      
213      return;
# Line 304 | Line 258 | namespace OpenMD {
258      
259      mixer.explicitlySet = true;
260  
261 <    pair<AtomType*, AtomType*> key1, key2;
262 <    key1 = make_pair(atype1, atype2);
263 <    key2 = make_pair(atype2, atype1);
264 <    
265 <    MixingMap[key1] = mixer;
266 <    if (key2 != key1) {
313 <      MixingMap[key2] = mixer;
261 >    int sctid1 = SCtids[ atype1->getIdent() ];
262 >    int sctid2 = SCtids[ atype2->getIdent() ];
263 >
264 >    MixingMap[sctid1][sctid2] = mixer;
265 >    if (sctid2 != sctid1) {
266 >      MixingMap[sctid2][sctid1] = mixer;
267      }    
268      return;
269    }
270  
271 <  void SC::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
319 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
271 >  void SC::calcDensity(InteractionData &idat) {
272      
273      if (!initialized_) initialize();
274 +    int sctid1 = SCtids[idat.atid1];
275 +    int sctid2 = SCtids[idat.atid2];
276      
277 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
277 >    SCInteractionData &mixer = MixingMap[sctid1][sctid2];
278  
279 <    rho_i_at_j = mixer.phi->getValueAt(rij);
326 <    rho_j_at_i = rho_i_at_j;
279 >    RealType rcij = mixer.rCut;
280  
281 +    if ( *(idat.rij)  < rcij) {
282 +      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
283 +      *(idat.rho1) += rho;
284 +      *(idat.rho2) += rho;
285 +    }
286 +    
287      return;
288    }
289  
290 <  void SC::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
332 <                          RealType &dfrhodrho) {
290 >  void SC::calcFunctional(SelfData &sdat) {
291  
292      if (!initialized_) initialize();
293  
294 <    SCAtomData data1 = SCMap[at1];
295 <    
296 <    frho = - data1.c * data1.epsilon * sqrt(rho);
297 <    dfrhodrho = 0.5 * frho / rho;
298 <    
341 <    return;
342 <  }
294 >    SCAtomData &data1 = SCdata[SCtids[sdat.atid]];
295 >  
296 >    RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
297 >    *(sdat.frho) = u;
298 >    *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
299  
300 <
301 <  void SC::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
302 <                     RealType rij, RealType r2, RealType sw,
303 <                     RealType &vpair, RealType &pot, Vector3d &f1,
348 <                     RealType rho_i, RealType rho_j,
349 <                     RealType dfrhodrho_i, RealType dfrhodrho_j,
350 <                     RealType &fshift_i, RealType &fshift_j) {
351 <    
352 <    if (!initialized_) initialize();
353 <    
354 <    SCAtomData data1 = SCMap[at1];
355 <    SCAtomData data2 = SCMap[at1];
300 >    (*(sdat.pot))[METALLIC_FAMILY] += u;
301 >    if (sdat.doParticlePot) {
302 >      *(sdat.particlePot) += u;
303 >    }
304  
305 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
358 <
359 <    RealType rcij = mixer.rCut;
360 <    RealType vcij = mixer.vCut;
361 <
362 <    pair<RealType, RealType> res;
363 <
364 <    res = mixer.phi->getValueAndDerivativeAt(rij);
365 <    RealType rhtmp = res.first;
366 <    RealType drhodr = res.second;
367 <
368 <    res = mixer.V->getValueAndDerivativeAt(rij);
369 <    RealType vptmp = res.first;
370 <    RealType dvpdr = res.second;
371 <
372 <    RealType pot_temp = vptmp - vcij;
373 <    vpair += pot_temp;
374 <
375 <    RealType dudr = drhodr * (dfrhodrho_i + dfrhodrho_j) + dvpdr;
376 <
377 <    f1 += d * dudr / rij;
378 <        
379 <    // particle_pot is the difference between the full potential
380 <    // and the full potential without the presence of a particular
381 <    // particle (atom1).
382 <    //
383 <    // This reduces the density at other particle locations, so
384 <    // we need to recompute the density at atom2 assuming atom1
385 <    // didn't contribute.  This then requires recomputing the
386 <    // density functional for atom2 as well.
387 <    //
388 <    // Most of the particle_pot heavy lifting comes from the
389 <    // pair interaction, and will be handled by vpair.
390 <    
391 <    fshift_i = - data1.c * data1.epsilon * sqrt(rho_i - rhtmp);
392 <    fshift_j = - data2.c * data2.epsilon * sqrt(rho_j - rhtmp);
393 <    
394 <    pot += pot_temp;
395 <
396 <    return;    
305 >    return;
306    }
398 –
399 –
400 –  void SC::calc_sc_prepair_rho(int *atid1, int *atid2, RealType *rij,
401 –                               RealType* rho_i_at_j, RealType* rho_j_at_i){
402 –    
403 –    if (!initialized_) initialize();
404 –    
405 –    AtomType* atype1 = SClist[*atid1];
406 –    AtomType* atype2 = SClist[*atid2];
407 –    
408 –    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
409 –    
410 –    return;    
411 –  }
307    
308 <  void SC::calc_sc_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
309 <                                 RealType *dfrhodrho) {
308 >
309 >  void SC::calcForce(InteractionData &idat) {
310      
311      if (!initialized_) initialize();
312      
313 <    AtomType* atype1 = SClist[*atid1];  
313 >    int &sctid1 = SCtids[idat.atid1];
314 >    int &sctid2 = SCtids[idat.atid2];
315  
316 <    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
317 <    
422 <    return;    
423 <  }
424 <  
425 <  RealType SC::getSCcut(int *atid1) {
316 >    SCAtomData &data1 = SCdata[sctid1];
317 >    SCAtomData &data2 = SCdata[sctid2];
318  
319 <    if (!initialized_) initialize();
428 <    
429 <    AtomType* atype1 = SClist[*atid1];  
430 <      
431 <    return 2.0 * getAlpha(atype1);
432 <  }
319 >    SCInteractionData &mixer = MixingMap[sctid1][sctid2];
320  
321 <  void SC::do_sc_pair(int *atid1, int *atid2, RealType *d, RealType *rij,
435 <                        RealType *r2, RealType *sw, RealType *vpair,
436 <                        RealType *pot, RealType *f1, RealType *rho1,
437 <                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
438 <                        RealType *fshift1, RealType *fshift2) {
321 >    RealType rcij = mixer.rCut;
322  
323 <    if (!initialized_) initialize();
324 <    
325 <    AtomType* atype1 = SClist[*atid1];
443 <    AtomType* atype2 = SClist[*atid2];
444 <    
445 <    Vector3d disp(d[0], d[1], d[2]);
446 <    Vector3d frc(f1[0], f1[1], f1[2]);
447 <    
448 <    calcForce(atype1, atype2, disp, *rij, *r2, *sw, *vpair,  *pot, frc,
449 <              *rho1, *rho2, *dfrho1, *dfrho2, *fshift1, *fshift2);
323 >    if ( *(idat.rij)  < rcij) {
324 >      RealType vcij = mixer.vCut;
325 >      RealType rhtmp, drhodr, vptmp, dvpdr;
326        
327 <    f1[0] = frc.x();
328 <    f1[1] = frc.y();
329 <    f1[2] = frc.z();
327 >      mixer.phi->getValueAndDerivativeAt( *(idat.rij), rhtmp, drhodr );      
328 >      mixer.V->getValueAndDerivativeAt( *(idat.rij), vptmp, dvpdr);
329 >      
330 >      RealType pot_temp = vptmp - vcij;
331 >      *(idat.vpair) += pot_temp;
332 >      
333 >      RealType dudr = drhodr * ( *(idat.dfrho1) + *(idat.dfrho2) ) + dvpdr;
334 >      
335 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
336 >        
337 >      if (idat.doParticlePot) {
338 >        // particlePot is the difference between the full potential and
339 >        // the full potential without the presence of a particular
340 >        // particle (atom1).
341 >        //
342 >        // This reduces the density at other particle locations, so we
343 >        // need to recompute the density at atom2 assuming atom1 didn't
344 >        // contribute.  This then requires recomputing the density
345 >        // functional for atom2 as well.
346 >        
347 >        *(idat.particlePot1) -= data2.c * data2.epsilon *
348 >          sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
349  
350 +        *(idat.particlePot2) -= data1.c * data1.epsilon *
351 +          sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
352 +      }
353 +      
354 +      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
355 +    }
356 +      
357      return;    
358    }
457 –  
458 –  void SC::setCutoffSC(RealType *thisRcut) {
459 –    scRcut_ = *thisRcut;
460 –  }
461 – }
359  
360 < extern "C" {
361 <  
362 < #define fortranCalcDensity FC_FUNC(calc_sc_prepair_rho, CALC_SC_PREPAIR_RHO)
363 < #define fortranCalcFunctional FC_FUNC(calc_sc_preforce_frho, CALC_SC_PREFORCE_FRHO)
364 < #define fortranCalcForce FC_FUNC(do_sc_pair, DO_SC_PAIR)
365 < #define fortranSetCutoffSC FC_FUNC(setcutoffsc, SETCUTOFFSC)
366 < #define fortranGetSCcut FC_FUNC(getsccut, GETSCCUT)
470 <  
471 <  
472 <  void fortranCalcDensity(int *atid1, int *atid2, RealType *rij,
473 <                          RealType *rho_i_at_j, RealType *rho_j_at_i) {
360 >  RealType SC::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
361 >    if (!initialized_) initialize();  
362 >
363 >    int atid1 = atypes.first->getIdent();
364 >    int atid2 = atypes.second->getIdent();
365 >    int &sctid1 = SCtids[atid1];
366 >    int &sctid2 = SCtids[atid2];
367      
368 <    return OpenMD::SC::Instance()->calc_sc_prepair_rho(atid1, atid2, rij,
369 <                                                       rho_i_at_j,  
370 <                                                       rho_j_at_i);
368 >    if (sctid1 == -1 || sctid2 == -1) {
369 >      return 0.0;
370 >    } else {
371 >      return MixingMap[sctid1][sctid2].rCut;
372 >    }
373    }
479 –  void fortranCalcFunctional(int *atid1, RealType *rho, RealType *frho,
480 –                             RealType *dfrhodrho) {  
481 –    
482 –    return OpenMD::SC::Instance()->calc_sc_preforce_Frho(atid1, rho, frho,
483 –                                                         dfrhodrho);
484 –    
485 –  }
486 –  void fortranSetCutoffSC(RealType *rcut) {
487 –    return OpenMD::SC::Instance()->setCutoffSC(rcut);
488 –  }
489 –  void fortranCalcForce(int *atid1, int *atid2, RealType *d, RealType *rij,
490 –                        RealType *r2, RealType *sw, RealType *vpair,
491 –                        RealType *pot, RealType *f1, RealType *rho1,
492 –                        RealType *rho2, RealType *dfrho1, RealType *dfrho2,
493 –                        RealType *fshift1, RealType *fshift2){
494 –    
495 –    return OpenMD::SC::Instance()->do_sc_pair(atid1, atid2, d, rij,
496 –                                              r2, sw, vpair,
497 –                                              pot, f1, rho1,
498 –                                              rho2, dfrho1, dfrho2,
499 –                                              fshift1,  fshift2);
500 –  }
501 –  RealType fortranGetSCcut(int* atid) {
502 –    return OpenMD::SC::Instance()->getSCcut(atid);
503 –  }
504 –  
374   }

Diff Legend

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