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 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC

# Line 44 | Line 44
44  
45   #include <cmath>
46   #include "nonbonded/GB.hpp"
47 – #include "nonbonded/LJ.hpp"
47   #include "utils/simError.h"
48  
49   using namespace std;
50   namespace OpenMD {
51  
52 <  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;
52 >  GB::GB() : name_("GB"), initialized_(false), mu_(2.0), nu_(1.0), forceField_(NULL) {}
53    
60 –  GB* GB::_instance = NULL;
61 –
62 –  GB* GB::Instance() {
63 –    if (!_instance) {
64 –      _instance = new GB();
65 –    }
66 –    return _instance;
67 –  }
68 –
54    GayBerneParam GB::getGayBerneParam(AtomType* atomType) {
55      
56      // Do sanity checking on the AtomType we were passed before
# Line 105 | Line 90 | namespace OpenMD {
90      return gbData->getData();
91    }
92  
93 +  LJParam GB::getLJParam(AtomType* atomType) {
94 +    
95 +    // Do sanity checking on the AtomType we were passed before
96 +    // building any data structures:
97 +    if (!atomType->isLennardJones()) {
98 +      sprintf( painCave.errMsg,
99 +               "GB::getLJParam was passed an atomType (%s) that does not\n"
100 +               "\tappear to be a Lennard-Jones atom.\n",
101 +               atomType->getName().c_str());
102 +      painCave.severity = OPENMD_ERROR;
103 +      painCave.isFatal = 1;
104 +      simError();
105 +    }
106 +    
107 +    GenericData* data = atomType->getPropertyByName("LennardJones");
108 +    if (data == NULL) {
109 +      sprintf( painCave.errMsg, "GB::getLJParam could not find Lennard-Jones\n"
110 +               "\tparameters for atomType %s.\n", atomType->getName().c_str());
111 +      painCave.severity = OPENMD_ERROR;
112 +      painCave.isFatal = 1;
113 +      simError();
114 +    }
115 +    
116 +    LJParamGenericData* ljData = dynamic_cast<LJParamGenericData*>(data);
117 +    if (ljData == NULL) {
118 +      sprintf( painCave.errMsg,
119 +               "GB::getLJParam could not convert GenericData to LJParam for\n"
120 +               "\tatom type %s\n", atomType->getName().c_str());
121 +      painCave.severity = OPENMD_ERROR;
122 +      painCave.isFatal = 1;
123 +      simError();          
124 +    }
125 +    
126 +    return ljData->getData();
127 +  }
128 +  
129 +  RealType GB::getLJEpsilon(AtomType* atomType) {    
130 +    LJParam ljParam = getLJParam(atomType);
131 +    return ljParam.epsilon;
132 +  }
133 +  RealType GB::getLJSigma(AtomType* atomType) {    
134 +    LJParam ljParam = getLJParam(atomType);
135 +    return ljParam.sigma;
136 +  }
137 +  
138    void GB::initialize() {    
139 +    
140 +    ForceFieldOptions& fopts = forceField_->getForceFieldOptions();
141 +    mu_ = fopts.getGayBerneMu();
142 +    nu_ = fopts.getGayBerneNu();
143      ForceField::AtomTypeContainer* atomTypes = forceField_->getAtomTypes();
144      ForceField::AtomTypeContainer::MapTypeIterator i;
145      AtomType* at;
# Line 148 | Line 182 | namespace OpenMD {
182        er1 = gb1.GB_eps_ratio;
183        dw1 = gb1.GB_dw;
184      } else if (atomType->isLennardJones()) {
185 <      d1 = LJ::Instance()->getSigma(atomType) / sqrt(2.0);
186 <      e1 = LJ::Instance()->getEpsilon(atomType);
185 >      d1 = getLJSigma(atomType) / sqrt(2.0);
186 >      e1 = getLJEpsilon(atomType);
187        l1 = d1;
188        er1 = 1.0;
189        dw1 = 1.0;      
# Line 181 | Line 215 | namespace OpenMD {
215          er2 = gb2.GB_eps_ratio;
216          dw2 = gb2.GB_dw;
217        } else if (atype2->isLennardJones()) {
218 <        d2 = LJ::Instance()->getSigma(atype2) / sqrt(2.0);
219 <        e2 = LJ::Instance()->getEpsilon(atype2);
218 >        d2 = getLJSigma(atype2) / sqrt(2.0);
219 >        e2 = getLJEpsilon(atype2);
220          l2 = d2;
221          er2 = 1.0;
222          dw2 = 1.0;
# Line 228 | Line 262 | namespace OpenMD {
262        }
263      }      
264    }
265 <  
265 >  
266 >  void GB::calcForce(InteractionData idat) {
267  
233 –  RealType GB::getGayBerneCut(int atid) {
268      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;
269      
270 <    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) {
269 <
270 <    if (!initialized_) initialize();
271 <    
272 <    pair<AtomType*, AtomType*> key = make_pair(at1, at2);
270 >    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
271      GBInteractionData mixer = MixingMap[key];
272  
273      RealType sigma0 = mixer.sigma0;
# Line 282 | Line 280 | namespace OpenMD {
280      RealType xpap2  = mixer.xpap2;
281      RealType xpapi2 = mixer.xpapi2;
282  
283 <    Vector3d ul1 = A1.getColumn(2);
284 <    Vector3d ul2 = A2.getColumn(2);
283 >    Vector3d ul1 = idat.A1.getRow(2);
284 >    Vector3d ul2 = idat.A2.getRow(2);
285  
286      RealType a, b, g;
287  
288 <    bool i_is_LJ = at1->isLennardJones();
289 <    bool j_is_LJ = at2->isLennardJones();
288 >    bool i_is_LJ = idat.atype1->isLennardJones();
289 >    bool j_is_LJ = idat.atype2->isLennardJones();
290  
291      if (i_is_LJ) {
292        a = 0.0;
293        ul1 = V3Zero;
294      } else {
295 <      a = dot(d, ul1);
295 >      a = dot(idat.d, ul1);
296      }
297  
298      if (j_is_LJ) {
299        b = 0.0;
300        ul2 = V3Zero;
301      } else {
302 <      b = dot(d, ul2);
302 >      b = dot(idat.d, ul2);
303      }
304  
305      if (i_is_LJ || j_is_LJ)
# Line 309 | Line 307 | namespace OpenMD {
307      else
308        g = dot(ul1, ul2);
309  
310 <    RealType au = a / r;
311 <    RealType bu = b / r;
310 >    RealType au = a / idat.rij;
311 >    RealType bu = b / idat.rij;
312      
313      RealType au2 = au * au;
314      RealType bu2 = bu * bu;
# Line 323 | Line 321 | namespace OpenMD {
321      RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
322      RealType e2 = 1.0 - Hp;
323      RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
324 <    RealType BigR = dw*sigma0 / (r - sigma + dw*sigma0);
324 >    RealType BigR = dw*sigma0 / (idat.rij - sigma + dw*sigma0);
325      
326      RealType R3 = BigR*BigR*BigR;
327      RealType R6 = R3*R3;
# Line 331 | Line 329 | namespace OpenMD {
329      RealType R12 = R6*R6;
330      RealType R13 = R6*R7;
331  
332 <    RealType U = vdwMult * 4.0 * eps * (R12 - R6);
332 >    RealType U = idat.vdwMult * 4.0 * eps * (R12 - R6);
333  
334      RealType s3 = sigma*sigma*sigma;
335      RealType s03 = sigma0*sigma0*sigma0;
336  
337 <    RealType pref1 = - vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * r);
337 >    RealType pref1 = - idat.vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * idat.rij);
338  
339 <    RealType pref2 = vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*r*s03);
339 >    RealType pref2 = idat.vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*idat.rij*s03);
340  
341 <    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*r/s3 + H));
341 >    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*idat.rij/s3 + H));
342      
343      RealType dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0 - xp2 * g2)
344        + pref2 * (xa2 * au - x2 *bu*g) / (1.0 - x2 * g2);
# Line 354 | Line 352 | namespace OpenMD {
352        (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
353      
354  
355 <    Vector3d rhat = d / r;  
356 <    Vector3d rxu1 = cross(d, ul1);
357 <    Vector3d rxu2 = cross(d, ul2);
355 >    Vector3d rhat = idat.d / idat.rij;  
356 >    Vector3d rxu1 = cross(idat.d, ul1);
357 >    Vector3d rxu2 = cross(idat.d, ul2);
358      Vector3d uxu = cross(ul1, ul2);
359      
360 <    pot += U*sw;
361 <    f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
362 <    t1 += dUda * rxu1 - dUdg * uxu;
363 <    t2 += dUdb * rxu2 - dUdg * uxu;
364 <    vpair += U*sw;
360 >    idat.pot += U*idat.sw;
361 >    idat.f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
362 >    idat.t1 += dUda * rxu1 - dUdg * uxu;
363 >    idat.t2 += dUdb * rxu2 - dUdg * uxu;
364 >    idat.vpair += U*idat.sw;
365  
366      return;
367  
368    }
369  
370 <  void GB::do_gb_pair(int *atid1, int *atid2, RealType *d, RealType *r,
371 <                      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:
370 >  RealType GB::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
371 >    if (!initialized_) initialize();  
372  
373 <    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();
373 >    RealType cut = 0.0;
374  
375 <    t1[0] = trq1.x();
376 <    t1[1] = trq1.y();
377 <    t1[2] = trq1.z();
375 >    if (at1->isGayBerne()) {
376 >      GayBerneParam gb1 = getGayBerneParam(at1);
377 >      RealType d1 = gb1.GB_d;
378 >      RealType l1 = gb1.GB_l;
379 >      // sigma is actually sqrt(2)*l  for prolate ellipsoids
380 >      cut = max(cut, 2.5 * sqrt(2.0) * max(d1, l1));
381 >    } else if (at1->isLennardJones()) {
382 >      cut = max(cut, 2.5 * getLJSigma(at1));
383 >    }
384  
385 <    t2[0] = trq2.x();
386 <    t2[1] = trq2.y();
387 <    t2[2] = trq2.z();
388 <
389 <    return;    
385 >    if (at2->isGayBerne()) {
386 >      GayBerneParam gb2 = getGayBerneParam(at2);
387 >      RealType d2 = gb2.GB_d;
388 >      RealType l2 = gb2.GB_l;
389 >      cut = max(cut, 2.5 * sqrt(2.0) * max(d2, l2));
390 >    } else if (at1->isLennardJones()) {
391 >      cut = max(cut, 2.5 * getLJSigma(at2));
392 >    }
393 >  
394 >    return cut;
395    }
396   }
397  
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