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

Comparing branches/development/src/nonbonded/GB.cpp (file contents):
Revision 1483 by gezelter, Tue Jul 27 21:17:31 2010 UTC vs.
Revision 1674 by gezelter, Thu Feb 16 15:59:20 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 44 | Line 45
45  
46   #include <cmath>
47   #include "nonbonded/GB.hpp"
47 – #include "nonbonded/LJ.hpp"
48   #include "utils/simError.h"
49  
50   using namespace std;
51   namespace OpenMD {
52  
53 <  bool GB::initialized_ = false;
54 <  RealType GB::mu_ = 2.0;
55 <  RealType GB::nu_ = 1.0;
56 <  ForceField* GB::forceField_ = NULL;
57 <  map<int, AtomType*> GB::GBMap;
58 <  map<pair<AtomType*, AtomType*>, GBInteractionData> GB::MixingMap;
53 >  GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
54    
60 –  GB* GB::_instance = NULL;
61 –
62 –  GB* GB::Instance() {
63 –    if (!_instance) {
64 –      _instance = new GB();
65 –    }
66 –    return _instance;
67 –  }
68 –
55    GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
56      
57      // Do sanity checking on the AtomType we were passed before
# Line 105 | Line 91 | namespace OpenMD {
91      return gbData->getData();
92    }
93  
94 +  LJParam GB::getLJParam(AtomType* atomType) {
95 +    
96 +    // Do sanity checking on the AtomType we were passed before
97 +    // building any data structures:
98 +    if (!atomType->isLennardJones()) {
99 +      sprintf( painCave.errMsg,
100 +               "GB::getLJParam was passed an atomType (%s) that does not\n"
101 +               "\tappear to be a Lennard-Jones atom.\n",
102 +               atomType->getName().c_str());
103 +      painCave.severity = OPENMD_ERROR;
104 +      painCave.isFatal = 1;
105 +      simError();
106 +    }
107 +    
108 +    GenericData* data = atomType->getPropertyByName("LennardJones");
109 +    if (data == NULL) {
110 +      sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
111 +               "\tparameters for atomType %s.\n", atomType->getName().c_str());
112 +      painCave.severity = OPENMD_ERROR;
113 +      painCave.isFatal = 1;
114 +      simError();
115 +    }
116 +    
117 +    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
118 +    if (ljData == NULL) {
119 +      sprintf( painCave.errMsg,
120 +               "GB::getLJParam could not convert GenericData to LJParam for\n"
121 +               "\tatom type %s\n", atomType->getName().c_str());
122 +      painCave.severity = OPENMD_ERROR;
123 +      painCave.isFatal = 1;
124 +      simError();          
125 +    }
126 +    
127 +    return ljData->getData();
128 +  }
129 +  
130 +  RealType GB::getLJEpsilon(AtomType* atomType) {    
131 +    LJParam ljParam = getLJParam(atomType);
132 +    return ljParam.epsilon;
133 +  }
134 +  RealType GB::getLJSigma(AtomType* atomType) {    
135 +    LJParam ljParam = getLJParam(atomType);
136 +    return ljParam.sigma;
137 +  }
138 +  
139    void GB::initialize() {    
140 +    
141 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
142 +    mu_ = fopts.getGayBerneMu();
143 +    nu_ = fopts.getGayBerneNu();
144      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
145      ForceField::AtomTypeContainer::MapTypeIterator i;
146      AtomType* at;
# Line 148 | Line 183 | namespace OpenMD {
183        er1 = gb1.GB_eps_ratio;
184        dw1 = gb1.GB_dw;
185      } else if (atomType->isLennardJones()) {
186 <      d1 = LJ::Instance()->getSigma(atomType) / sqrt(2.0);
187 <      e1 = LJ::Instance()->getEpsilon(atomType);
186 >      d1 = getLJSigma(atomType) / sqrt(2.0);
187 >      e1 = getLJEpsilon(atomType);
188        l1 = d1;
189        er1 = 1.0;
190        dw1 = 1.0;      
# Line 181 | Line 216 | namespace OpenMD {
216          er2 = gb2.GB_eps_ratio;
217          dw2 = gb2.GB_dw;
218        } else if (atype2->isLennardJones()) {
219 <        d2 = LJ::Instance()->getSigma(atype2) / sqrt(2.0);
220 <        e2 = LJ::Instance()->getEpsilon(atype2);
219 >        d2 = getLJSigma(atype2) / sqrt(2.0);
220 >        e2 = getLJEpsilon(atype2);
221          l2 = d2;
222          er2 = 1.0;
223          dw2 = 1.0;
224        }
225                        
226 <      GBInteractionData mixer;        
226 >      GBInteractionData mixer1, mixer2;    
227        
228        //  Cleaver paper uses sqrt of squares to get sigma0 for
229        //  mixed interactions.
230              
231 <      mixer.sigma0 = sqrt(d1*d1 + d2*d2);
232 <      mixer.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
233 <      mixer.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
234 <      mixer.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
231 >      mixer1.sigma0 = sqrt(d1*d1 + d2*d2);
232 >      mixer1.xa2 = (l1*l1 - d1*d1)/(l1*l1 + d2*d2);
233 >      mixer1.xai2 = (l2*l2 - d2*d2)/(l2*l2 + d1*d1);
234 >      mixer1.x2 = (l1*l1 - d1*d1) * (l2*l2 - d2*d2) /
235          ((l2*l2 + d1*d1) * (l1*l1 + d2*d2));
236 +
237 +      mixer2.sigma0 = mixer1.sigma0;
238 +      // xa2 and xai2 for j-i pairs are reversed from the same i-j pairing.
239 +      // Swapping the particles reverses the anisotropy parameters:
240 +      mixer2.xa2 = mixer1.xai2;
241 +      mixer2.xai2 = mixer1.xa2;
242 +      mixer2.x2 = mixer1.x2;
243  
244        // assumed LB mixing rules for now:
245  
246 <      mixer.dw = 0.5 * (dw1 + dw2);
247 <      mixer.eps0 = sqrt(e1 * e2);
246 >      mixer1.dw = 0.5 * (dw1 + dw2);
247 >      mixer1.eps0 = sqrt(e1 * e2);
248 >
249 >      mixer2.dw = mixer1.dw;
250 >      mixer2.eps0 = mixer1.eps0;
251        
252        RealType er = sqrt(er1 * er2);
253 <      RealType ermu = pow(er,(1.0 / mu_));
253 >      RealType ermu = pow(er, (RealType(1.0) / mu_));
254        RealType xp = (1.0 - ermu) / (1.0 + ermu);
255        RealType ap2 = 1.0 / (1.0 + ermu);
256        
257 <      mixer.xp2 = xp * xp;
258 <      mixer.xpap2 = xp * ap2;
259 <      mixer.xpapi2 = xp / ap2;
257 >      mixer1.xp2 = xp * xp;
258 >      mixer1.xpap2 = xp * ap2;
259 >      mixer1.xpapi2 = xp / ap2;
260  
261 +      mixer2.xp2 = mixer1.xp2;
262 +      mixer2.xpap2 = mixer1.xpap2;
263 +      mixer2.xpapi2 = mixer1.xpapi2;
264 +
265        // only add this pairing if at least one of the atoms is a Gay-Berne atom
266  
267        if (atomType->isGayBerne() || atype2->isGayBerne()) {
# Line 221 | Line 270 | namespace OpenMD {
270          key1 = make_pair(atomType, atype2);
271          key2 = make_pair(atype2, atomType);
272          
273 <        MixingMap[key1] = mixer;
273 >        MixingMap[key1] = mixer1;
274          if (key2 != key1) {
275 <          MixingMap[key2] = mixer;
275 >          MixingMap[key2] = mixer2;
276          }
277        }
278      }      
279    }
280 <  
281 <
233 <  RealType GB::getGayBerneCut(int atid) {
234 <    if (!initialized_) initialize();
235 <    std::map<int, AtomType*> :: const_iterator it;
236 <    it = GBMap.find(atid);
237 <    if (it == GBMap.end()) {
238 <      sprintf( painCave.errMsg,
239 <               "GB::getGayBerneCut could not find atid %d in GBMap\n",
240 <               (atid));
241 <      painCave.severity = OPENMD_ERROR;
242 <      painCave.isFatal = 1;
243 <      simError();          
244 <    }
245 <
246 <    AtomType* atype = it->second;
247 <
248 <    RealType gbCut;
249 <    
250 <    if (atype->isGayBerne()) {
251 <      GayBerneParam gb = getGayBerneParam(atype);
252 <
253 <      // sigma is actually sqrt(2) * l for prolate ellipsoids
254 <      gbCut = 2.5 * sqrt(2.0) * max(gb.GB_l, gb.GB_d);
255 <
256 <    } else if (atype->isLennardJones()) {
257 <      gbCut = 2.5 * LJ::Instance()->getSigma(atype);
258 <    }
259 <    
260 <    return gbCut;
261 <  }
262 <
263 <
264 <  void GB::calcForce(AtomType* at1, AtomType* at2, Vector3d d,
265 <                     RealType r, RealType r2, RealType sw,
266 <                     RealType vdwMult, RealType &vpair, RealType &pot,
267 <                     RotMat3x3d A1, RotMat3x3d A2, Vector3d &f1,
268 <                     Vector3d &t1, Vector3d &t2) {
280 >  
281 >  void GB::calcForce(InteractionData &idat) {
282  
283      if (!initialized_) initialize();
284      
285 <    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
273 <    GBInteractionData mixer = MixingMap[key];
285 >    GBInteractionData mixer = MixingMap[idat.atypes];
286  
287      RealType sigma0 = mixer.sigma0;
288      RealType dw     = mixer.dw;
# Line 282 | Line 294 | namespace OpenMD {
294      RealType xpap2  = mixer.xpap2;
295      RealType xpapi2 = mixer.xpapi2;
296  
297 <    Vector3d ul1 = A1.getColumn(2);
298 <    Vector3d ul2 = A2.getColumn(2);
297 >    Vector3d ul1 = idat.A1->getRow(2);
298 >    Vector3d ul2 = idat.A2->getRow(2);
299  
300      RealType a, b, g;
301  
302 <    bool i_is_LJ = at1->isLennardJones();
303 <    bool j_is_LJ = at2->isLennardJones();
304 <
302 >    bool i_is_LJ = idat.atypes.first->isLennardJones();
303 >    bool j_is_LJ = idat.atypes.second->isLennardJones();
304 >    
305      if (i_is_LJ) {
306        a = 0.0;
307        ul1 = V3Zero;
308      } else {
309 <      a = dot(d, ul1);
309 >      a = dot(*(idat.d), ul1);
310      }
311  
312      if (j_is_LJ) {
313        b = 0.0;
314        ul2 = V3Zero;
315      } else {
316 <      b = dot(d, ul2);
316 >      b = dot(*(idat.d), ul2);
317      }
318  
319      if (i_is_LJ || j_is_LJ)
# Line 309 | Line 321 | namespace OpenMD {
321      else
322        g = dot(ul1, ul2);
323  
324 <    RealType au = a / r;
325 <    RealType bu = b / r;
324 >    RealType au = a / *(idat.rij);
325 >    RealType bu = b / *(idat.rij);
326      
327      RealType au2 = au * au;
328      RealType bu2 = bu * bu;
# Line 323 | Line 335 | namespace OpenMD {
335      RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
336      RealType e2 = 1.0 - Hp;
337      RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
338 <    RealType BigR = dw*sigma0 / (r - sigma + dw*sigma0);
338 >    RealType BigR = dw*sigma0 / (*(idat.rij) - sigma + dw*sigma0);
339      
340      RealType R3 = BigR*BigR*BigR;
341      RealType R6 = R3*R3;
# Line 331 | Line 343 | namespace OpenMD {
343      RealType R12 = R6*R6;
344      RealType R13 = R6*R7;
345  
346 <    RealType U = vdwMult * 4.0 * eps * (R12 - R6);
346 >    RealType U = *(idat.vdwMult) * 4.0 * eps * (R12 - R6);
347  
348      RealType s3 = sigma*sigma*sigma;
349      RealType s03 = sigma0*sigma0*sigma0;
350  
351 <    RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);
351 >    RealType pref1 = - *(idat.vdwMult) * 8.0 * eps * mu_ * (R12 - R6) /
352 >      (e2 * *(idat.rij));
353  
354 <    RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*r*s03);
354 >    RealType pref2 = *(idat.vdwMult) * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /
355 >      (dw*  *(idat.rij) * s03);
356  
357 <    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H));
357 >    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *  
358 >                                             *(idat.rij) / s3 + H));
359      
360      RealType dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0 - xp2 * g2)
361        + pref2 * (xa2 * au - x2 *bu*g) / (1.0 - x2 * g2);
# Line 352 | Line 367 | namespace OpenMD {
367        + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) /
368        (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
369        (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
355 –    
370  
371 <    Vector3d rhat = d / r;  
372 <    Vector3d rxu1 = cross(d, ul1);
373 <    Vector3d rxu2 = cross(d, ul2);
371 >    Vector3d rhat = *(idat.d) / *(idat.rij);  
372 >    Vector3d rxu1 = cross(*(idat.d), ul1);
373 >    Vector3d rxu2 = cross(*(idat.d), ul2);
374      Vector3d uxu = cross(ul1, ul2);
361 –    
362 –    pot += U*sw;
363 –    f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
364 –    t1 += dUda * rxu1 - dUdg * uxu;
365 –    t2 += dUdb * rxu2 - dUdg * uxu;
366 –    vpair += U*sw;
375  
376 +    (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
377 +    *(idat.f1) += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
378 +    *(idat.t1) += dUda * rxu1 - dUdg * uxu;
379 +    *(idat.t2) += dUdb * rxu2 + dUdg * uxu;
380 +    *(idat.vpair) += U * *(idat.sw);
381 +
382      return;
383  
384    }
385  
386 <  void GB::do_gb_pair(int *atid1, int *atid2, RealType *d, RealType *r,
387 <                      RealType *r2, RealType *sw, RealType *vdwMult,
374 <                      RealType *vpair, RealType *pot, RealType *A1,
375 <                      RealType *A2, RealType *f1, RealType *t1, RealType *t2) {
376 <    
377 <    if (!initialized_) initialize();
378 <    
379 <    AtomType* atype1 = GBMap[*atid1];
380 <    AtomType* atype2 = GBMap[*atid2];
381 <    
382 <    Vector3d disp(d);
383 <    Vector3d frc(f1);
384 <    Vector3d trq1(t1);
385 <    Vector3d trq2(t2);
386 <    RotMat3x3d Ai(A1);
387 <    RotMat3x3d Aj(A2);
388 <  
389 <    // Fortran has the opposite matrix ordering from c++, so we'll use
390 <    // transpose here.  When we finish the conversion to C++, this wrapper
391 <    // will disappear, as will the transpose below:
386 >  RealType GB::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
387 >    if (!initialized_) initialize();  
388  
389 <    calcForce(atype1, atype2, disp, *r, *r2, *sw, *vdwMult, *vpair, *pot,
394 <              Ai.transpose(), Aj.transpose(), frc, trq1, trq1);
395 <      
396 <    f1[0] = frc.x();
397 <    f1[1] = frc.y();
398 <    f1[2] = frc.z();
389 >    RealType cut = 0.0;
390  
391 <    t1[0] = trq1.x();
392 <    t1[1] = trq1.y();
393 <    t1[2] = trq1.z();
391 >    if (atypes.first->isGayBerne()) {
392 >      GayBerneParam gb1 = getGayBerneParam(atypes.first);
393 >      RealType d1 = gb1.GB_d;
394 >      RealType l1 = gb1.GB_l;
395 >      // sigma is actually sqrt(2)*l  for prolate ellipsoids
396 >      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
397 >    } else if (atypes.first->isLennardJones()) {
398 >      cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
399 >    }
400  
401 <    t2[0] = trq2.x();
402 <    t2[1] = trq2.y();
403 <    t2[2] = trq2.z();
404 <
405 <    return;    
401 >    if (atypes.second->isGayBerne()) {
402 >      GayBerneParam gb2 = getGayBerneParam(atypes.second);
403 >      RealType d2 = gb2.GB_d;
404 >      RealType l2 = gb2.GB_l;
405 >      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
406 >    } else if (atypes.second->isLennardJones()) {
407 >      cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
408 >    }
409 >  
410 >    return cut;
411    }
412   }
413  
412 – extern "C" {
413 –  
414 – #define fortranGetGayBerneCut FC_FUNC(getgaybernecut, GETGAYBERNECUT)
415 – #define fortranDoGBPair FC_FUNC(do_gb_pair, DO_GB_PAIR)
416 –  
417 –  RealType fortranGetGayBerneCut(int* atid) {
418 –    return OpenMD::GB::Instance()->getGayBerneCut(*atid);
419 –  }
420 –
421 –  void fortranDoGBPair(int *atid1, int *atid2, RealType *d, RealType *r,
422 –                       RealType *r2, RealType *sw, RealType *vdwMult,
423 –                       RealType *vpair, RealType *pot, RealType *A1,
424 –                       RealType *A2, RealType *f1, RealType *t1, RealType *t2){
425 –    
426 –    return OpenMD::GB::Instance()->do_gb_pair(atid1, atid2, d, r, r2, sw,
427 –                                              vdwMult, vpair, pot, A1, A2, f1,
428 –                                              t1, t2);
429 –  }
430 – }

Diff Legend

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