ComputeNonbondedFEP.C File Reference

#include "ComputeNonbondedInl.h"
#include "ComputeNonbondedBase.h"

Go to the source code of this file.

Defines

#define FEPFLAG
#define CALCENERGY
#define NBTYPE   NBPAIR
#define FULLELECT
#define MERGEELECT
#define SLOWONLY
#define NBTYPE   NBSELF
#define FULLELECT
#define MERGEELECT
#define SLOWONLY

Functions

void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, BigReal cutoff2, BigReal myVdwLambda, BigReal myVdwLambda2, Bool Fep_WCA_repuOn, Bool Fep_WCA_dispOn, bool Fep_Wham, BigReal WCA_rcut1, BigReal WCA_rcut2, BigReal WCA_rcut3, BigReal switchfactor, Bool vdwForceSwitching, Bool LJcorrection, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_energy_2, BigReal *alch_vdw_energy_2_Left)


Define Documentation

#define CALCENERGY

Definition at line 188 of file ComputeNonbondedFEP.C.

#define FEPFLAG

Definition at line 187 of file ComputeNonbondedFEP.C.

#define FULLELECT

Definition at line 205 of file ComputeNonbondedFEP.C.

#define FULLELECT

Definition at line 205 of file ComputeNonbondedFEP.C.

#define MERGEELECT

Definition at line 207 of file ComputeNonbondedFEP.C.

#define MERGEELECT

Definition at line 207 of file ComputeNonbondedFEP.C.

#define NBTYPE   NBSELF

Definition at line 203 of file ComputeNonbondedFEP.C.

#define NBTYPE   NBPAIR

Definition at line 203 of file ComputeNonbondedFEP.C.

#define SLOWONLY

Definition at line 210 of file ComputeNonbondedFEP.C.

#define SLOWONLY

Definition at line 210 of file ComputeNonbondedFEP.C.


Function Documentation

void fep_vdw_forceandenergies ( BigReal  A,
BigReal  B,
BigReal  r2,
BigReal  myVdwShift,
BigReal  myVdwShift2,
BigReal  switchdist2,
BigReal  cutoff2,
BigReal  myVdwLambda,
BigReal  myVdwLambda2,
Bool  Fep_WCA_repuOn,
Bool  Fep_WCA_dispOn,
bool  Fep_Wham,
BigReal  WCA_rcut1,
BigReal  WCA_rcut2,
BigReal  WCA_rcut3,
BigReal  switchfactor,
Bool  vdwForceSwitching,
Bool  LJcorrection,
BigReal alch_vdw_energy,
BigReal alch_vdw_force,
BigReal alch_vdw_energy_2,
BigReal alch_vdw_energy_2_Left 
) [inline]

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 19 of file ComputeNonbondedFEP.C.

References f.

00026                                    {
00027   // switching function (this is correct whether switching is active or not)
00028   const BigReal switchmul = (r2 > switchdist2 ? switchfactor*(cutoff2 - r2) \
00029                              *(cutoff2 - r2)                            \
00030                              *(cutoff2 - 3.*switchdist2 + 2.*r2) : 1.);
00031   const BigReal switchmul2 = (r2 > switchdist2 ?                        \
00032                               12.*switchfactor*(cutoff2 - r2)           \
00033                               *(r2 - switchdist2) : 0.);
00034   
00035   *alch_vdw_energy_2_Left = 0.;
00036 
00037   if(Fep_WCA_repuOn){
00038     if(B != 0.0) {// Excluded when some atoms like H, drude, lone pair are involved
00039       const BigReal Emin = B*B/(4.0*A);
00040       const BigReal Rmin_SQ = powf(2.0*A/B, 1.f/3);
00041       const BigReal r2_1 = r2 + (1.-WCA_rcut1)*(1.-WCA_rcut1)*Rmin_SQ;
00042       const BigReal r2_2 = r2 + (1.-WCA_rcut2)*(1.-WCA_rcut2)*Rmin_SQ;
00043       const BigReal r2_3 = r2 + (1.-WCA_rcut3)*(1.-WCA_rcut3)*Rmin_SQ;
00044       const BigReal r6_1 = r2_1*r2_1*r2_1;
00045       const BigReal r6_2 = r2_2*r2_2*r2_2;
00046       const BigReal r6_3 = r2_3*r2_3*r2_3;
00047       const BigReal WCA_rcut1_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut1) \
00048                                                        *(1.0-WCA_rcut1)) ? \
00049                                         A/(r6_1*r6_1) - B/r6_1 + Emin : 0.);
00050       const BigReal WCA_rcut2_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut2) \
00051                                                        *(1.0-WCA_rcut2)) ? \
00052                                         A/(r6_2*r6_2) - B/r6_2 + Emin : 0.);
00053       const BigReal WCA_rcut3_energy = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut3) \
00054                                                        *(1.0-WCA_rcut3)) ? \
00055                                         A/(r6_3*r6_3) - B/r6_3 + Emin : 0.);
00056       const BigReal WCA_rcut1_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut1) \
00057                                                       *(1.0-WCA_rcut1)) ? \
00058                                        (12.*(WCA_rcut1_energy)          \
00059                                         + 6.*B/r6_1 - 12.0 * Emin )/r2_1: 0.);
00060       const BigReal WCA_rcut2_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut2) \
00061                                                       *(1.0-WCA_rcut2))? \
00062                                        (12.*(WCA_rcut2_energy)          \
00063                                         + 6.*B/r6_2 - 12.0 * Emin )/r2_2: 0.);
00064       const BigReal WCA_rcut3_force = (r2 <= Rmin_SQ*(1.0 - (1.0-WCA_rcut3) \
00065                                                       *(1.0-WCA_rcut3)) ? \
00066                                        (12.*(WCA_rcut3_energy)          \
00067                                         + 6.*B/r6_3 - 12.0 * Emin )/r2_3: 0.);
00068       // separation-shifted repulsion force and energy
00069       *alch_vdw_energy = WCA_rcut2_energy; 
00070       *alch_vdw_force = WCA_rcut2_force;
00071       if(WCA_rcut1 < WCA_rcut2) {
00072         *alch_vdw_energy_2_Left = *alch_vdw_energy + WCA_rcut2_energy - WCA_rcut1_energy; 
00073       }
00074       if(WCA_rcut2 < WCA_rcut3) {
00075         *alch_vdw_energy_2 = *alch_vdw_energy + WCA_rcut3_energy - WCA_rcut2_energy; 
00076       }
00077     }
00078     else {
00079       *alch_vdw_energy = 0.0;
00080       *alch_vdw_force = 0.0;
00081       *alch_vdw_energy_2_Left = 0.0;
00082       *alch_vdw_energy_2 = 0.0;
00083     }
00084   }
00085   else if(Fep_WCA_dispOn) {
00086     // separation-shifted dispersion force and energy
00087     if(B == 0.0) {      // some atoms like H, drude, lone pair are involved
00088       *alch_vdw_energy = 0.0;
00089       *alch_vdw_force = 0.0;
00090       *alch_vdw_energy_2 = 0.0;
00091     }
00092     else {
00093       const BigReal Emin = B*B/(4.0*A);
00094       const BigReal Rmin_SQ = powf(2.0*A/B, 1.f/3);
00095       const BigReal r2_1 = r2;
00096       const BigReal r2_2 = r2;
00097       const BigReal r6_1 = r2_1*r2_1*r2_1;
00098       const BigReal r6_2 = r2_2*r2_2*r2_2;
00099       *alch_vdw_energy = r2 > Rmin_SQ? \
00100                          myVdwLambda*(A/(r6_1*r6_1) - B/r6_1): \
00101                          A/(r6_1*r6_1) - B/r6_1 + (1.-myVdwLambda)* Emin;
00102       *alch_vdw_force =  r2 > Rmin_SQ? \
00103                          (12.*(*alch_vdw_energy) + 6.*myVdwLambda*B/r6_1)/r2_1 * switchmul \
00104                          + (*alch_vdw_energy) * switchmul2:(12.*(*alch_vdw_energy) \
00105                          + 6.*B/r6_1 - 12.0*(1.-myVdwLambda)* Emin )/r2_1;
00106       *alch_vdw_energy *= switchmul; 
00107                         if(!Fep_Wham){ 
00108         *alch_vdw_energy_2 = r2 > Rmin_SQ? \
00109                              myVdwLambda2*switchmul*(A/(r6_1*r6_1) - B/r6_1): \
00110                              A/(r6_1*r6_1) - B/r6_1 + (1.-myVdwLambda2)* Emin;
00111                         }
00112                         else{
00113                                 *alch_vdw_energy_2 = r2 > Rmin_SQ? \
00114                                                            switchmul*(A/(r6_1*r6_1) - B/r6_1): - Emin;
00115         *alch_vdw_energy_2 += *alch_vdw_energy;
00116                         }
00117     }
00118   } 
00119   else {
00120     //myVdwShift already multplied by relevant (1-vdwLambda)
00121     const BigReal r2_1 = 1./(r2 + myVdwShift);
00122     const BigReal r2_2 = 1./(r2 + myVdwShift2);
00123     const BigReal r6_1 = r2_1*r2_1*r2_1;
00124     const BigReal r6_2 = r2_2*r2_2*r2_2;
00125     // separation-shifted vdW force and energy
00126     const BigReal U1 = A*r6_1*r6_1 - B*r6_1; // NB: unscaled, shorthand only!
00127     const BigReal U2 = A*r6_2*r6_2 - B*r6_2;
00128     *alch_vdw_energy = myVdwLambda*switchmul*U1;
00129     *alch_vdw_energy_2 = myVdwLambda2*switchmul*U2;
00130     *alch_vdw_force = (myVdwLambda*(switchmul*(12.*U1 + 6.*B*r6_1)*r2_1 \
00131                                     + switchmul2*U1));
00132     // BKR - separation-shifted vdW force switching and potential shifting
00133     if(vdwForceSwitching){ // add potential shifts term
00134       const BigReal cutoff6 = cutoff2*cutoff2*cutoff2;
00135       const BigReal cutoff3 = cutoff2*sqrt(cutoff2);
00136 
00137       const BigReal shifted_switchdist2 = switchdist2 + myVdwShift;
00138       const BigReal shifted_switchdist = sqrt(shifted_switchdist2);
00139       const BigReal shifted_switchdist6 = \
00140           shifted_switchdist2*shifted_switchdist2*shifted_switchdist2;
00141       const BigReal shifted_switchdist3 = shifted_switchdist2*shifted_switchdist;
00142      
00143       const BigReal shifted_switchdist2_2 = switchdist2 + myVdwShift2;
00144       const BigReal shifted_switchdist_2 = sqrt(shifted_switchdist2_2);
00145       const BigReal shifted_switchdist6_2 = \
00146           shifted_switchdist2_2*shifted_switchdist2_2*shifted_switchdist2_2;
00147       const BigReal shifted_switchdist3_2 = shifted_switchdist2_2*shifted_switchdist_2;
00148 
00149       const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
00150       const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
00151       const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
00152  
00153       const BigReal v_vdwa2 = -A / (cutoff6*shifted_switchdist6_2);
00154       const BigReal v_vdwb2 = -B / (cutoff3*shifted_switchdist3_2);
00155       const BigReal dU2 = v_vdwa2 - v_vdwb2; //deltaV2 from Steinbach & Brooks
00156 
00157       if(r2 > switchdist2) {
00158         const BigReal k_vdwa = A*cutoff6 / (cutoff6 - shifted_switchdist6);
00159         const BigReal k_vdwb = B*cutoff3 / (cutoff3 - shifted_switchdist3);
00160         const BigReal r_1 = sqrt(r2_1);
00161         const BigReal tmpa = r6_1 - (1./cutoff6);
00162         const BigReal tmpb = r2_1*r_1 - (1./cutoff3);
00163 
00164         const BigReal k_vdwa2 = A*cutoff6 / (cutoff6 - shifted_switchdist6_2);
00165         const BigReal k_vdwb2 = B*cutoff3 / (cutoff3 - shifted_switchdist3_2);
00166         const BigReal r_2 = sqrt(r2_2);
00167         const BigReal tmpa2 = r6_2 - (1./cutoff6);
00168         const BigReal tmpb2 = r2_2*r_2 - (1./cutoff3);
00169 
00170         *alch_vdw_energy = (myVdwLambda*(k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb
00171                                          - (LJcorrection ? dU : 0.)));
00172         *alch_vdw_energy_2 = (myVdwLambda2*(k_vdwa2*tmpa2*tmpa2 \
00173                                             - k_vdwb2*tmpb2*tmpb2
00174                                             - (LJcorrection ? dU2 : 0.)));
00175         *alch_vdw_force = (myVdwLambda*r2_1*(12.*k_vdwa*tmpa*r6_1
00176                                              - 6.*k_vdwb*tmpb*r2_1*r_1));
00177       }else{
00178         if(!LJcorrection) {
00179           *alch_vdw_energy += myVdwLambda*dU;
00180           *alch_vdw_energy_2 += myVdwLambda2*dU2;
00181         }
00182       }
00183     }
00184   }
00185 }


Generated on Fri Sep 22 01:17:15 2017 for NAMD by  doxygen 1.4.7