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 1505 by gezelter, Sun Oct 3 22:18:59 2010 UTC vs.
Revision 1668 by gezelter, Fri Jan 6 19:03:05 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 239 | Line 240 | namespace OpenMD {
240        mixer.eps0 = sqrt(e1 * e2);
241        
242        RealType er = sqrt(er1 * er2);
243 <      RealType ermu = pow(er,(1.0 / mu_));
243 >      RealType ermu = pow(er, (RealType(1.0) / mu_));
244        RealType xp = (1.0 - ermu) / (1.0 + ermu);
245        RealType ap2 = 1.0 / (1.0 + ermu);
246        
247        mixer.xp2 = xp * xp;
248        mixer.xpap2 = xp * ap2;
249        mixer.xpapi2 = xp / ap2;
250 +
251 +      cerr << "mixer" << er1 << " " << er2 << " " << mu_ << " " << ermu << " " << xp <<" " << ap2 << "\n";
252  
253        // only add this pairing if at least one of the atoms is a Gay-Berne atom
254  
# Line 263 | Line 266 | namespace OpenMD {
266      }      
267    }
268    
269 <  void GB::calcForce(InteractionData idat) {
269 >  void GB::calcForce(InteractionData &idat) {
270  
271      if (!initialized_) initialize();
272      
273 <    pair<AtomType*, AtomType*> key = make_pair(idat.atype1, idat.atype2);
271 <    GBInteractionData mixer = MixingMap[key];
273 >    GBInteractionData mixer = MixingMap[idat.atypes];
274  
275      RealType sigma0 = mixer.sigma0;
276      RealType dw     = mixer.dw;
# Line 280 | Line 282 | namespace OpenMD {
282      RealType xpap2  = mixer.xpap2;
283      RealType xpapi2 = mixer.xpapi2;
284  
285 <    Vector3d ul1 = idat.A1.getRow(2);
286 <    Vector3d ul2 = idat.A2.getRow(2);
285 >    Vector3d ul1 = idat.A1->getRow(2);
286 >    Vector3d ul2 = idat.A2->getRow(2);
287  
288      RealType a, b, g;
289  
290 <    bool i_is_LJ = idat.atype1->isLennardJones();
291 <    bool j_is_LJ = idat.atype2->isLennardJones();
292 <
290 >    bool i_is_LJ = idat.atypes.first->isLennardJones();
291 >    bool j_is_LJ = idat.atypes.second->isLennardJones();
292 >    
293      if (i_is_LJ) {
294        a = 0.0;
295        ul1 = V3Zero;
296      } else {
297 <      a = dot(idat.d, ul1);
297 >      a = dot(*(idat.d), ul1);
298      }
299  
300      if (j_is_LJ) {
301        b = 0.0;
302        ul2 = V3Zero;
303      } else {
304 <      b = dot(idat.d, ul2);
304 >      b = dot(*(idat.d), ul2);
305      }
306  
307      if (i_is_LJ || j_is_LJ)
# Line 307 | Line 309 | namespace OpenMD {
309      else
310        g = dot(ul1, ul2);
311  
312 <    RealType au = a / idat.rij;
313 <    RealType bu = b / idat.rij;
312 >    cerr << "in GB, d = " << *(idat.d) << "\n";
313 >    cerr << "abg = " << a << " " << b << " " << g <<"\n";
314 >
315 >    RealType au = a / *(idat.rij);
316 >    RealType bu = b / *(idat.rij);
317      
318      RealType au2 = au * au;
319      RealType bu2 = bu * bu;
# Line 316 | Line 321 | namespace OpenMD {
321      
322      RealType H  = (xa2 * au2 + xai2 * bu2 - 2.0*x2*au*bu*g)  / (1.0 - x2*g2);
323      RealType Hp = (xpap2*au2 + xpapi2*bu2 - 2.0*xp2*au*bu*g) / (1.0 - xp2*g2);
324 +    cerr << "xa2, xai2 " << xa2 << " " << xai2 << "\n";
325 +    cerr << "xpap2, xpapi2 " << xpap2 << " " << xpapi2 << "\n";
326 +    cerr << "H Hp = " << H <<  " " << Hp << "\n";
327  
328      RealType sigma = sigma0 / sqrt(1.0 - H);
329      RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
330      RealType e2 = 1.0 - Hp;
331      RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
332 <    RealType BigR = dw*sigma0 / (idat.rij - sigma + dw*sigma0);
332 >    cerr << "eps = " << eps0 << " " << e1 << " " << nu_ << " " << e2 << " " << mu_ << "\n";
333 >    RealType BigR = dw*sigma0 / (*(idat.rij) - sigma + dw*sigma0);
334      
335      RealType R3 = BigR*BigR*BigR;
336      RealType R6 = R3*R3;
# Line 329 | Line 338 | namespace OpenMD {
338      RealType R12 = R6*R6;
339      RealType R13 = R6*R7;
340  
341 <    RealType U = idat.vdwMult * 4.0 * eps * (R12 - R6);
341 >    RealType U = *(idat.vdwMult) * 4.0 * eps * (R12 - R6);
342  
343 +    cerr << "R12, R6, eps = " << R12 << " " << R6 << " " << eps << " " <<  *(idat.vdwMult) << "\n";
344 +
345      RealType s3 = sigma*sigma*sigma;
346      RealType s03 = sigma0*sigma0*sigma0;
347  
348 <    RealType pref1 = - idat.vdwMult * 8.0 * eps * mu_ * (R12 - R6) / (e2 * idat.rij);
348 >    RealType pref1 = - *(idat.vdwMult) * 8.0 * eps * mu_ * (R12 - R6) /
349 >      (e2 * *(idat.rij));
350  
351 <    RealType pref2 = idat.vdwMult * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /(dw*idat.rij*s03);
351 >    RealType pref2 = *(idat.vdwMult) * 8.0 * eps * s3 * (6.0*R13 - 3.0*R7) /
352 >      (dw*  *(idat.rij) * s03);
353  
354 <    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0*sigma0*idat.rij/s3 + H));
354 >    RealType dUdr = - (pref1 * Hp + pref2 * (sigma0 * sigma0 *  
355 >                                             *(idat.rij) / s3 + H));
356      
357      RealType dUda = pref1 * (xpap2*au - xp2*bu*g) / (1.0 - xp2 * g2)
358        + pref2 * (xa2 * au - x2 *bu*g) / (1.0 - x2 * g2);
# Line 350 | Line 364 | namespace OpenMD {
364        + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) /
365        (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
366        (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
353 –    
367  
368 <    Vector3d rhat = idat.d / idat.rij;  
369 <    Vector3d rxu1 = cross(idat.d, ul1);
370 <    Vector3d rxu2 = cross(idat.d, ul2);
368 >    cerr << pref1 << " " << pref2 << " " << dUdr <<" " << dUda << " " << dUdb << dUdg << "\n";
369 >
370 >    Vector3d rhat = *(idat.d) / *(idat.rij);  
371 >    Vector3d rxu1 = cross(*(idat.d), ul1);
372 >    Vector3d rxu2 = cross(*(idat.d), ul2);
373      Vector3d uxu = cross(ul1, ul2);
374 +
375 +    cerr << "U = " << U << "\n";
376 +    cerr << "f1 = " << dUdr * rhat + dUda * ul1 + dUdb * ul2 << "\n";
377 +    cerr << "t1 = " << dUda * rxu1 - dUdg * uxu << "\n";
378 +    cerr << "t2 = " << dUdb * rxu2 - dUdg * uxu << "\n";
379 +
380      
381 <    idat.pot += U*idat.sw;
382 <    idat.f1 += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
383 <    idat.t1 += dUda * rxu1 - dUdg * uxu;
384 <    idat.t2 += dUdb * rxu2 - dUdg * uxu;
385 <    idat.vpair += U*idat.sw;
381 >    (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
382 >    *(idat.f1) += dUdr * rhat + dUda * ul1 + dUdb * ul2;    
383 >    *(idat.t1) += dUda * rxu1 - dUdg * uxu;
384 >    *(idat.t2) += dUdb * rxu2 - dUdg * uxu;
385 >    *(idat.vpair) += U * *(idat.sw);
386  
387      return;
388  
389    }
390  
391 <  RealType GB::getSuggestedCutoffRadius(AtomType* at1, AtomType* at2) {
391 >  RealType GB::getSuggestedCutoffRadius(pair<AtomType*, AtomType*> atypes) {
392      if (!initialized_) initialize();  
393  
394      RealType cut = 0.0;
395  
396 <    if (at1->isGayBerne()) {
397 <      GayBerneParam gb1 = getGayBerneParam(at1);
396 >    if (atypes.first->isGayBerne()) {
397 >      GayBerneParam gb1 = getGayBerneParam(atypes.first);
398        RealType d1 = gb1.GB_d;
399        RealType l1 = gb1.GB_l;
400        // sigma is actually sqrt(2)*l  for prolate ellipsoids
401 <      cut = max(cut, 2.5 * sqrt(2.0) * max(d1, l1));
402 <    } else if (at1->isLennardJones()) {
403 <      cut = max(cut, 2.5 * getLJSigma(at1));
401 >      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d1, l1));
402 >    } else if (atypes.first->isLennardJones()) {
403 >      cut = max(cut, RealType(2.5) * getLJSigma(atypes.first));
404      }
405  
406 <    if (at2->isGayBerne()) {
407 <      GayBerneParam gb2 = getGayBerneParam(at2);
406 >    if (atypes.second->isGayBerne()) {
407 >      GayBerneParam gb2 = getGayBerneParam(atypes.second);
408        RealType d2 = gb2.GB_d;
409        RealType l2 = gb2.GB_l;
410 <      cut = max(cut, 2.5 * sqrt(2.0) * max(d2, l2));
411 <    } else if (at1->isLennardJones()) {
412 <      cut = max(cut, 2.5 * getLJSigma(at2));
410 >      cut = max(cut, RealType(2.5) * sqrt(RealType(2.0)) * max(d2, l2));
411 >    } else if (atypes.second->isLennardJones()) {
412 >      cut = max(cut, RealType(2.5) * getLJSigma(atypes.second));
413      }
414    
415      return cut;

Diff Legend

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