ComputeNonbondedFEP.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 lambda2 energy for FEP */
00018 /* ******************************************** */
00019 inline void fep_vdw_forceandenergies (BigReal A, BigReal B, BigReal r2, 
00020   BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, 
00021   BigReal cutoff2, BigReal myVdwLambda, BigReal myVdwLambda2, 
00022   Bool Fep_WCA_repuOn, Bool Fep_WCA_dispOn, bool Fep_Wham, BigReal WCA_rcut1, 
00023   BigReal WCA_rcut2, BigReal WCA_rcut3, BigReal switchfactor, 
00024   Bool vdwForceSwitching, Bool LJcorrection, BigReal* alch_vdw_energy,
00025   BigReal* alch_vdw_force, BigReal* alch_vdw_energy_2,
00026   BigReal* alch_vdw_energy_2_Left) {
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 }
00186 
00187 #define FEPFLAG
00188 #define CALCENERGY
00189 
00190 #define NBTYPE NBPAIR
00191 #include "ComputeNonbondedBase.h"
00192 #define FULLELECT
00193 #include "ComputeNonbondedBase.h"
00194 #define MERGEELECT
00195 #include "ComputeNonbondedBase.h"
00196 #undef MERGEELECT
00197 #define SLOWONLY
00198 #include "ComputeNonbondedBase.h"
00199 #undef SLOWONLY
00200 #undef FULLELECT
00201 #undef  NBTYPE
00202 
00203 #define NBTYPE NBSELF
00204 #include "ComputeNonbondedBase.h"
00205 #define FULLELECT
00206 #include "ComputeNonbondedBase.h"
00207 #define MERGEELECT
00208 #include "ComputeNonbondedBase.h"
00209 #undef MERGEELECT
00210 #define SLOWONLY
00211 #include "ComputeNonbondedBase.h"
00212 #undef SLOWONLY
00213 #undef FULLELECT
00214 #undef  NBTYPE
00215 
00216 #undef CALCENERGY
00217 #undef FEPFLAG
00218 
00219 

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