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 1586 by gezelter, Tue Jun 14 20:41:44 2011 UTC vs.
Revision 1587 by gezelter, Fri Jul 8 20:25:32 2011 UTC

# Line 239 | Line 239 | namespace OpenMD {
239        mixer.eps0 = sqrt(e1 * e2);
240        
241        RealType er = sqrt(er1 * er2);
242 <      RealType ermu = pow(er,(1.0 / mu_));
242 >      RealType ermu = pow(er, (1.0 / mu_));
243        RealType xp = (1.0 - ermu) / (1.0 + ermu);
244        RealType ap2 = 1.0 / (1.0 + ermu);
245        
# Line 247 | Line 247 | namespace OpenMD {
247        mixer.xpap2 = xp * ap2;
248        mixer.xpapi2 = xp / ap2;
249  
250 +      cerr << "mixer" << er1 << " " << er2 << " " << mu_ << " " << ermu << " " << xp <<" " << ap2 << "\n";
251 +
252        // only add this pairing if at least one of the atoms is a Gay-Berne atom
253  
254        if (atomType->isGayBerne() || atype2->isGayBerne()) {
# Line 306 | Line 308 | namespace OpenMD {
308      else
309        g = dot(ul1, ul2);
310  
311 +    cerr << "in GB, d = " << *(idat.d) << "\n";
312 +    cerr << "abg = " << a << " " << b << " " << g <<"\n";
313 +
314      RealType au = a / *(idat.rij);
315      RealType bu = b / *(idat.rij);
316      
# Line 315 | Line 320 | namespace OpenMD {
320      
321      RealType H  = (xa2 * au2 + xai2 * bu2 - 2.0*x2*au*bu*g)  / (1.0 - x2*g2);
322      RealType Hp = (xpap2*au2 + xpapi2*bu2 - 2.0*xp2*au*bu*g) / (1.0 - xp2*g2);
323 +    cerr << "xa2, xai2 " << xa2 << " " << xai2 << "\n";
324 +    cerr << "xpap2, xpapi2 " << xpap2 << " " << xpapi2 << "\n";
325 +    cerr << "H Hp = " << H <<  " " << Hp << "\n";
326  
327      RealType sigma = sigma0 / sqrt(1.0 - H);
328      RealType e1 = 1.0 / sqrt(1.0 - x2*g2);
329      RealType e2 = 1.0 - Hp;
330      RealType eps = eps0 * pow(e1,nu_) * pow(e2,mu_);
331 +    cerr << "eps = " << eps0 << " " << e1 << " " << nu_ << " " << e2 << " " << mu_ << "\n";
332      RealType BigR = dw*sigma0 / (*(idat.rij) - sigma + dw*sigma0);
333      
334      RealType R3 = BigR*BigR*BigR;
# Line 330 | Line 339 | namespace OpenMD {
339  
340      RealType U = *(idat.vdwMult) * 4.0 * eps * (R12 - R6);
341  
342 +    cerr << "R12, R6, eps = " << R12 << " " << R6 << " " << eps << " " <<  *(idat.vdwMult) << "\n";
343 +
344      RealType s3 = sigma*sigma*sigma;
345      RealType s03 = sigma0*sigma0*sigma0;
346  
# Line 352 | Line 363 | namespace OpenMD {
363        + 8.0 * eps * mu_ * (R12 - R6) * (xp2*au*bu - Hp*xp2*g) /
364        (1.0 - xp2 * g2) / e2 + 8.0 * eps * s3 * (3.0 * R7 - 6.0 * R13) *
365        (x2 * au * bu - H * x2 * g) / (1.0 - x2 * g2) / (dw * s03);
355 –    
366  
367 +    cerr << pref1 << " " << pref2 << " " << dUdr <<" " << dUda << " " << dUdb << dUdg << "\n";
368 +
369      Vector3d rhat = *(idat.d) / *(idat.rij);  
370      Vector3d rxu1 = cross(*(idat.d), ul1);
371      Vector3d rxu2 = cross(*(idat.d), ul2);
372      Vector3d uxu = cross(ul1, ul2);
373 +
374 +    cerr << "U = " << U << "\n";
375 +    cerr << "f1 = " << dUdr * rhat + dUda * ul1 + dUdb * ul2 << "\n";
376 +    cerr << "t1 = " << dUda * rxu1 - dUdg * uxu << "\n";
377 +    cerr << "t2 = " << dUdb * rxu2 - dUdg * uxu << "\n";
378 +
379      
380      (*(idat.pot))[VANDERWAALS_FAMILY] += U *  *(idat.sw);
381      *(idat.f1) += dUdr * rhat + dUda * ul1 + dUdb * ul2;    

Diff Legend

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