ComputeNonbondedTI.C

Go to the documentation of this file.
00001 
00007 /*
00008    Common operations for ComputeNonbonded classes
00009 */
00010 
00011 #include "ComputeNonbondedInl.h"
00012 
00013 // 3 inline functions to handle explicit calculation of separation-shifted
00014 // vdW for FEP and TI, and a shifted electrostatics potential for decoupling
00015 
00016 /* ********************************** */
00017 /* vdW energy, force and dU/dl for TI */
00018 /* ********************************** */
00019 inline void ti_vdw_force_energy_dUdl (BigReal A, BigReal B, BigReal r2, 
00020   BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2, 
00021   BigReal myVdwLambda, BigReal alchVdwShiftCoeff, BigReal switchfactor, 
00022   Bool vdwForceSwitching, Bool LJcorrection, BigReal* alch_vdw_energy,
00023   BigReal* alch_vdw_force, BigReal* alch_vdw_dUdl) {
00024   //myVdwShift already multplied by relevant (1-vdwLambda)
00025   const BigReal r2_1 = 1./(r2 + myVdwShift);
00026   const BigReal r6_1 = r2_1*r2_1*r2_1;
00027     
00028   // switching function (this is correct whether switching is active or not)
00029   const BigReal switchmul = (r2 > switchdist2 ? switchfactor*(cutoff2 - r2) \
00030                              *(cutoff2 - r2) \
00031                              *(cutoff2 - 3.*switchdist2 + 2.*r2) : 1.);
00032   const BigReal switchmul2 = (r2 > switchdist2 ?                      \
00033                               12.*switchfactor*(cutoff2 - r2)       \
00034                               *(r2 - switchdist2) : 0.);
00035 
00036   // separation-shifted vdW force and energy
00037   const BigReal U = A*r6_1*r6_1 - B*r6_1; // NB: unscaled! for shorthand only!
00038   *alch_vdw_energy = myVdwLambda*switchmul*U;
00039   *alch_vdw_force = (myVdwLambda*(switchmul*(12.*U + 6.*B*r6_1)*r2_1 \
00040                                   + switchmul2*U));
00041   *alch_vdw_dUdl = (switchmul*(U + myVdwLambda*alchVdwShiftCoeff \
00042                                    *(6.*U + 3.*B*r6_1)*r2_1));
00043 
00044   // BKR - separation-shifted vdW force switching and potential shifting
00045   if(vdwForceSwitching){ // add potential shifts and additional dU/dl terms
00046     const BigReal cutoff6 = cutoff2*cutoff2*cutoff2;
00047     const BigReal cutoff3 = cutoff2*sqrt(cutoff2);
00048     const BigReal shifted_switchdist2 = switchdist2 + myVdwShift;
00049     const BigReal shifted_switchdist = sqrt(shifted_switchdist2);
00050     const BigReal shifted_switchdist6 = \
00051         shifted_switchdist2*shifted_switchdist2*shifted_switchdist2;
00052     const BigReal shifted_switchdist3 = shifted_switchdist2*shifted_switchdist;
00053 
00054     const BigReal v_vdwa = -A / (cutoff6*shifted_switchdist6);
00055     const BigReal v_vdwb = -B / (cutoff3*shifted_switchdist3);
00056     const BigReal dU = v_vdwa - v_vdwb; //deltaV2 from Steinbach & Brooks
00057 
00058     if(r2 > switchdist2) {
00059       const BigReal k_vdwa = A*cutoff6 / (cutoff6 - shifted_switchdist6);
00060       const BigReal k_vdwb = B*cutoff3 / (cutoff3 - shifted_switchdist3);
00061       const BigReal tmpa = r6_1 - (1./cutoff6);
00062       const BigReal r_1 = sqrt(r2_1);
00063       const BigReal tmpb = r2_1*r_1 - (1./cutoff3);
00064       const BigReal Uh = (k_vdwa*tmpa*tmpa - k_vdwb*tmpb*tmpb \
00065                           - (LJcorrection ? dU : 0.));
00066       *alch_vdw_energy = myVdwLambda*Uh;
00067       *alch_vdw_force = (myVdwLambda*r2_1*(12.*k_vdwa*tmpa*r6_1
00068                                            - 6.*k_vdwb*tmpb*r2_1*r_1));
00069       *alch_vdw_dUdl = (Uh + alchVdwShiftCoeff*myVdwLambda*r2_1 \
00070                         *(6.*k_vdwa*tmpa*r6_1                   \
00071                           - 3.*k_vdwb*tmpb*r2_1*r_1));
00072     }else{
00073       if(!LJcorrection) {
00074         *alch_vdw_energy += myVdwLambda*dU;
00075         *alch_vdw_dUdl += dU;
00076       }
00077     }
00078   }
00079 }
00080 
00081 
00082 /*************THERMODYNAMIC INTEGRATION*************/
00083 #define TIFLAG
00084 #define CALCENERGY
00085 
00086 #define NBTYPE NBPAIR
00087 #include "ComputeNonbondedBase.h"
00088 #define FULLELECT
00089 #include "ComputeNonbondedBase.h"
00090 #define MERGEELECT
00091 #include "ComputeNonbondedBase.h"
00092 #undef MERGEELECT
00093 #define SLOWONLY
00094 #include "ComputeNonbondedBase.h"
00095 #undef SLOWONLY
00096 #undef FULLELECT
00097 #undef  NBTYPE
00098 
00099 #define NBTYPE NBSELF
00100 #include "ComputeNonbondedBase.h"
00101 #define FULLELECT
00102 #include "ComputeNonbondedBase.h"
00103 #define MERGEELECT
00104 #include "ComputeNonbondedBase.h"
00105 #undef MERGEELECT
00106 #define SLOWONLY
00107 #include "ComputeNonbondedBase.h"
00108 #undef SLOWONLY
00109 #undef FULLELECT
00110 #undef  NBTYPE
00111 
00112 #undef CALCENERGY
00113 
00114 #define NBTYPE NBPAIR
00115 #include "ComputeNonbondedBase.h"
00116 #define FULLELECT
00117 #include "ComputeNonbondedBase.h"
00118 #define MERGEELECT
00119 #include "ComputeNonbondedBase.h"
00120 #undef MERGEELECT
00121 #define SLOWONLY
00122 #include "ComputeNonbondedBase.h"
00123 #undef SLOWONLY
00124 #undef FULLELECT
00125 #undef  NBTYPE
00126 
00127 #define NBTYPE NBSELF
00128 #include "ComputeNonbondedBase.h"
00129 #define FULLELECT
00130 #include "ComputeNonbondedBase.h"
00131 #define MERGEELECT
00132 #include "ComputeNonbondedBase.h"
00133 #undef MERGEELECT
00134 #define SLOWONLY
00135 #include "ComputeNonbondedBase.h"
00136 #undef SLOWONLY
00137 #undef FULLELECT
00138 #undef  NBTYPE
00139 
00140 #undef TIFLAG
00141 

Generated on Thu Sep 21 01:17:11 2017 for NAMD by  doxygen 1.4.7