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 1489 by gezelter, Tue Aug 10 18:34:59 2010 UTC vs.
Revision 1868 by gezelter, Tue Apr 30 15:56:54 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  
53 <  bool SC::initialized_ = false;
54 <  RealType SC::scRcut_ = 0.0;
55 <  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;
53 >
54 >  SC::SC() : name_("SC"), initialized_(false), forceField_(NULL),
55 >             scRcut_(0.0), np_(3000) {}
56    
57 <  SC* SC::_instance = NULL;
57 >  SC::~SC() {
58 >    initialized_ = false;
59  
60 <  SC* SC::Instance() {
61 <    if (!_instance) {
62 <      _instance = new SC();
60 >    map<pair<AtomType*, AtomType*>, SCInteractionData>::iterator it;
61 >    for (it = MixingMap.begin(); it != MixingMap.end(); ++it) {
62 >      SCInteractionData mixer = (*it).second;
63 >      delete mixer.V;
64 >      delete mixer.phi;
65      }
66 –    return _instance;
67 –  }
68 –  
69 –  SCParam SC::getSCParam(AtomType* atomType) {
70 –    
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 –  }
66  
67 <  RealType SC::getC(AtomType* atomType) {    
68 <    SCParam scParam = getSCParam(atomType);
69 <    return scParam.c;
67 >    MixingMap.clear();
68 >    SCMap.clear();
69 >    SClist.clear();
70    }
71 <
111 <  RealType SC::getM(AtomType* atomType) {    
112 <    SCParam scParam = getSCParam(atomType);
113 <    return scParam.m;
114 <  }
115 <
71 >        
72    RealType SC::getM(AtomType* atomType1, AtomType* atomType2) {    
73 <    RealType m1 = getM(atomType1);
74 <    RealType m2 = getM(atomType2);
73 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
74 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
75 >    RealType m1 = sca1.getM();
76 >    RealType m2 = sca2.getM();
77      return 0.5 * (m1 + m2);
78    }
79  
122 –  RealType SC::getN(AtomType* atomType) {    
123 –    SCParam scParam = getSCParam(atomType);
124 –    return scParam.n;
125 –  }
126 –
80    RealType SC::getN(AtomType* atomType1, AtomType* atomType2) {    
81 <    RealType n1 = getN(atomType1);
82 <    RealType n2 = getN(atomType2);
81 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
82 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
83 >    RealType n1 = sca1.getN();
84 >    RealType n2 = sca2.getN();
85      return 0.5 * (n1 + n2);
86    }
87  
133 –  RealType SC::getAlpha(AtomType* atomType) {    
134 –    SCParam scParam = getSCParam(atomType);
135 –    return scParam.alpha;
136 –  }
137 –
88    RealType SC::getAlpha(AtomType* atomType1, AtomType* atomType2) {    
89 <    RealType alpha1 = getAlpha(atomType1);
90 <    RealType alpha2 = getAlpha(atomType2);
89 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
90 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
91 >    RealType alpha1 = sca1.getAlpha();
92 >    RealType alpha2 = sca2.getAlpha();
93  
94      ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
95      std::string DistanceMix = fopts.getDistanceMixingRule();
# Line 149 | Line 101 | namespace OpenMD {
101        return 0.5 * (alpha1 + alpha2);
102    }
103  
104 <  RealType SC::getEpsilon(AtomType* atomType) {    
105 <    SCParam scParam = getSCParam(atomType);
106 <    return scParam.epsilon;
107 <  }
108 <
157 <  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {    
158 <    RealType epsilon1 = getEpsilon(atomType1);
159 <    RealType epsilon2 = getEpsilon(atomType2);
104 >  RealType SC::getEpsilon(AtomType* atomType1, AtomType* atomType2) {  
105 >    SuttonChenAdapter sca1 = SuttonChenAdapter(atomType1);
106 >    SuttonChenAdapter sca2 = SuttonChenAdapter(atomType2);
107 >    RealType epsilon1 = sca1.getEpsilon();
108 >    RealType epsilon2 = sca2.getEpsilon();
109      return sqrt(epsilon1 * epsilon2);
110    }
111  
# Line 168 | Line 117 | namespace OpenMD {
117  
118      for (at = atomTypes->beginType(i); at != NULL;
119           at = atomTypes->nextType(i)) {
120 <      if (at->isSC())
120 >      SuttonChenAdapter sca = SuttonChenAdapter(at);
121 >      if (sca.isSuttonChen())
122          addType(at);
123      }    
124      initialized_ = true;
# Line 178 | Line 128 | namespace OpenMD {
128  
129    void SC::addType(AtomType* atomType){
130  
131 +    SuttonChenAdapter sca = SuttonChenAdapter(atomType);
132      SCAtomData scAtomData;
133      
134 <    scAtomData.c = getC(atomType);
135 <    scAtomData.m = getM(atomType);
136 <    scAtomData.n = getN(atomType);
137 <    scAtomData.alpha = getAlpha(atomType);
138 <    scAtomData.epsilon = getEpsilon(atomType);
134 >    scAtomData.c = sca.getC();
135 >    scAtomData.m = sca.getM();
136 >    scAtomData.n = sca.getN();
137 >    scAtomData.alpha = sca.getAlpha();
138 >    scAtomData.epsilon = sca.getEpsilon();
139      scAtomData.rCut = 2.0 * scAtomData.alpha;
140  
141      // add it to the map:
191 –    AtomTypeProperties atp = atomType->getATP();    
142  
143      pair<map<int,AtomType*>::iterator,bool> ret;    
144 <    ret = SClist.insert( pair<int, AtomType*>(atp.ident, atomType) );
144 >    ret = SClist.insert( pair<int, AtomType*>(atomType->getIdent(), atomType) );
145      if (ret.second == false) {
146        sprintf( painCave.errMsg,
147                 "SC already had a previous entry with ident %d\n",
148 <               atp.ident);
148 >               atomType->getIdent() );
149        painCave.severity = OPENMD_INFO;
150        painCave.isFatal = 0;
151        simError();        
# Line 315 | Line 265 | namespace OpenMD {
265      return;
266    }
267  
268 <  void SC::calcDensity(AtomType* at1, AtomType* at2, const RealType rij,
319 <                        RealType &rho_i_at_j, RealType &rho_j_at_i) {
268 >  void SC::calcDensity(InteractionData &idat) {
269      
270      if (!initialized_) initialize();
271      
272 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
272 >    SCInteractionData mixer = MixingMap[ idat.atypes ];
273  
274 <    rho_i_at_j = mixer.phi->getValueAt(rij);
326 <    rho_j_at_i = rho_i_at_j;
274 >    RealType rcij = mixer.rCut;
275  
276 +    if ( *(idat.rij)  < rcij) {
277 +      RealType rho = mixer.phi->getValueAt( *(idat.rij) );
278 +      *(idat.rho1) += rho;
279 +      *(idat.rho2) += rho;
280 +    }
281 +    
282      return;
283    }
284  
285 <  void SC::calcFunctional(AtomType* at1, RealType rho, RealType &frho,
332 <                          RealType &dfrhodrho) {
285 >  void SC::calcFunctional(SelfData &sdat) {
286  
287      if (!initialized_) initialize();
288  
289 <    SCAtomData data1 = SCMap[at1];
290 <    
291 <    frho = - data1.c * data1.epsilon * sqrt(rho);
292 <    dfrhodrho = 0.5 * frho / rho;
293 <    
289 >    SCAtomData data1 = SCMap[sdat.atype];
290 >  
291 >    RealType u = - data1.c * data1.epsilon * sqrt( *(sdat.rho) );
292 >    *(sdat.frho) = u;
293 >    *(sdat.dfrhodrho) = 0.5 * *(sdat.frho) / *(sdat.rho);
294 >
295 >    (*(sdat.pot))[METALLIC_FAMILY] += u;
296 >    if (sdat.doParticlePot) {
297 >      *(sdat.particlePot) += u;
298 >    }
299 >
300      return;
301    }
302 <
302 >  
303  
304 <  void SC::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
346 <                     RealType rij, RealType r2, RealType sw,
347 <                     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) {
304 >  void SC::calcForce(InteractionData &idat) {
305      
306      if (!initialized_) initialize();
307      
308 <    SCAtomData data1 = SCMap[at1];
309 <    SCAtomData data2 = SCMap[at1];
308 >    SCAtomData data1 = SCMap[idat.atypes.first];
309 >    SCAtomData data2 = SCMap[idat.atypes.second];
310  
311 <    SCInteractionData mixer = MixingMap[make_pair(at1, at2)];
311 >    SCInteractionData mixer = MixingMap[idat.atypes];
312  
313      RealType rcij = mixer.rCut;
360 –    RealType vcij = mixer.vCut;
314  
315 <    pair<RealType, RealType> res;
316 <
317 <    res = mixer.phi->getValueAndDerivativeAt(rij);
318 <    RealType rhtmp = res.first;
319 <    RealType drhodr = res.second;
320 <
321 <    res = mixer.V->getValueAndDerivativeAt(rij);
322 <    RealType vptmp = res.first;
323 <    RealType dvpdr = res.second;
324 <
325 <    RealType pot_temp = vptmp - vcij;
326 <    vpair += pot_temp;
327 <
328 <    RealType dudr = drhodr * (dfrhodrho_i + dfrhodrho_j) + dvpdr;
329 <
330 <    f1 += d * dudr / rij;
315 >    if ( *(idat.rij)  < rcij) {
316 >      RealType vcij = mixer.vCut;
317 >      
318 >      pair<RealType, RealType> res;
319 >      
320 >      res = mixer.phi->getValueAndDerivativeAt( *(idat.rij) );
321 >      RealType rhtmp = res.first;
322 >      RealType drhodr = res.second;
323 >      
324 >      res = mixer.V->getValueAndDerivativeAt( *(idat.rij) );
325 >      RealType vptmp = res.first;
326 >      RealType dvpdr = res.second;
327 >      
328 >      RealType pot_temp = vptmp - vcij;
329 >      *(idat.vpair) += pot_temp;
330 >      
331 >      RealType dudr = drhodr * ( *(idat.dfrho1) + *(idat.dfrho2) ) + dvpdr;
332 >      
333 >      *(idat.f1) += *(idat.d) * dudr / *(idat.rij) ;
334          
335 <    // particle_pot is the difference between the full potential
336 <    // and the full potential without the presence of a particular
337 <    // particle (atom1).
338 <    //
339 <    // This reduces the density at other particle locations, so
340 <    // we need to recompute the density at atom2 assuming atom1
341 <    // didn't contribute.  This then requires recomputing the
342 <    // density functional for atom2 as well.
343 <    //
344 <    // Most of the particle_pot heavy lifting comes from the
345 <    // pair interaction, and will be handled by vpair.
346 <    
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;
335 >      if (idat.doParticlePot) {
336 >        // particlePot is the difference between the full potential and
337 >        // the full potential without the presence of a particular
338 >        // particle (atom1).
339 >        //
340 >        // This reduces the density at other particle locations, so we
341 >        // need to recompute the density at atom2 assuming atom1 didn't
342 >        // contribute.  This then requires recomputing the density
343 >        // functional for atom2 as well.
344 >        
345 >        *(idat.particlePot1) -= data2.c * data2.epsilon *
346 >          sqrt( *(idat.rho2) - rhtmp) + *(idat.frho2);
347  
348 +        *(idat.particlePot2) -= data1.c * data1.epsilon *
349 +          sqrt( *(idat.rho1) - rhtmp) + *(idat.frho1);
350 +      }
351 +      
352 +      (*(idat.pot))[METALLIC_FAMILY] += pot_temp;
353 +    }
354 +      
355      return;    
356    }
357  
358 +  RealType SC::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
359 +    if (!initialized_) initialize();  
360  
361 <  void SC::calc_sc_prepair_rho(int *atid1, int *atid2, RealType *rij,
362 <                               RealType* rho_i_at_j, RealType* rho_j_at_i){
363 <    
364 <    if (!initialized_) initialize();
365 <    
366 <    AtomType* atype1 = SClist[*atid1];
367 <    AtomType* atype2 = SClist[*atid2];
368 <    
408 <    calcDensity(atype1, atype2, *rij, *rho_i_at_j, *rho_j_at_i);
409 <    
410 <    return;    
361 >    map<pair<AtomType*, AtomType*>, SCInteractionData>::iterator it;
362 >    it = MixingMap.find(atypes);
363 >    if (it == MixingMap.end())
364 >      return 0.0;
365 >    else  {
366 >      SCInteractionData mixer = (*it).second;
367 >      return mixer.rCut;
368 >    }
369    }
412 –  
413 –  void SC::calc_sc_preforce_Frho(int *atid1, RealType *rho, RealType *frho,
414 –                                 RealType *dfrhodrho) {
415 –    
416 –    if (!initialized_) initialize();
417 –    
418 –    AtomType* atype1 = SClist[*atid1];  
419 –
420 –    calcFunctional(atype1, *rho, *frho, *dfrhodrho);
421 –    
422 –    return;    
423 –  }
424 –  
425 –  RealType SC::getSCcut(int *atid1) {
426 –
427 –    if (!initialized_) initialize();
428 –    
429 –    AtomType* atype1 = SClist[*atid1];  
430 –      
431 –    return 2.0 * getAlpha(atype1);
432 –  }
433 –
434 –  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) {
439 –
440 –    if (!initialized_) initialize();
441 –    
442 –    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);
450 –      
451 –    f1[0] = frc.x();
452 –    f1[1] = frc.y();
453 –    f1[2] = frc.z();
454 –
455 –    return;    
456 –  }
457 –  
458 –  void SC::setCutoffSC(RealType *thisRcut) {
459 –    scRcut_ = *thisRcut;
460 –  }
370   }
462 –
463 – extern "C" {
464 –  
465 – #define fortranCalcDensity FC_FUNC(calc_sc_prepair_rho, CALC_SC_PREPAIR_RHO)
466 – #define fortranCalcFunctional FC_FUNC(calc_sc_preforce_frho, CALC_SC_PREFORCE_FRHO)
467 – #define fortranCalcForce FC_FUNC(do_sc_pair, DO_SC_PAIR)
468 – #define fortranSetCutoffSC FC_FUNC(setcutoffsc, SETCUTOFFSC)
469 – #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) {
474 –    
475 –    return OpenMD::SC::Instance()->calc_sc_prepair_rho(atid1, atid2, rij,
476 –                                                       rho_i_at_j,  
477 –                                                       rho_j_at_i);
478 –  }
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 –  
505 – }

Diff Legend

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