NAMD
ComputeNonbondedFEP.C
Go to the documentation of this file.
1 
7 /*
8  Common operations for ComputeNonbonded classes
9 */
10 
11 #include "ComputeNonbondedInl.h"
12 #include "ComputeNonbondedAlch.h"
13 
14 /* ******************************************** */
15 /* vdW energy, force and lambda2 energy for FEP */
16 /* ******************************************** */
18  BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2,
19  BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching,
20  BigReal myVdwLambda, BigReal myVdwLambda2, Bool alchWCAOn,
21  BigReal myRepLambda, BigReal myRepLambda2, BigReal* alch_vdw_energy,
22  BigReal* alch_vdw_force, BigReal* alch_vdw_energy_2) {
23  BigReal U, F, U_2, dU, dU_2, switchmul, switchmul2;
24  /*
25  * The variable naming here is unavoidably awful. A number immediately after
26  * a variable implies a distance to that power. The suffix "_2" indicates an
27  * evaluation at alchLambda2. The prefix "inv" means the inverse (previously
28  * this also used underscores - it was awful).
29  */
30  if (alchWCAOn) {
31  // WCA-on, auxilliary, lambda-dependent cutoff based on Rmin
32  //
33  // Avoid divide by zero - corectly zeroes interaction below.
34  const BigReal Rmin2 = (B <= 0.0 ? 0.0 : powf(2.0*A/B, 1.f/3));
35  if (myRepLambda < 1.0) {
36  // modified repulsive soft-core
37  // All of the scaling is baked into the shift and cutoff values. There
38  // are no linear coupling terms out front.
39  const BigReal WCAshift = Rmin2*(1 - myRepLambda)*(1 - myRepLambda);
40  if (r2 <= Rmin2 - WCAshift) {
41  const BigReal epsilon = B*B/(4.0*A);
42  vdw_forceandenergy(A, B, r2 + WCAshift, &U, &F);
43  *alch_vdw_energy = U + epsilon;
44  *alch_vdw_force = F;
45  } else {
46  *alch_vdw_energy = 0.0;
47  *alch_vdw_force = 0.0;
48  }
49  } else { // myRepLambda == 1.0
50  if (vdwForceSwitching) {
51  // force and potential switching
52  if (r2 <= Rmin2) {
53  // normal LJ, potential w/two shifts
54  const BigReal epsilon = B*B/(4.0*A);
55  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
56  vdw_forceandenergy(A, B, r2, &U, &F);
57  *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon + myVdwLambda*dU;
58  *alch_vdw_force = F;
59  } else if (r2 <= switchdist2) {
60  // normal LJ potential w/shift
61  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
62  vdw_forceandenergy(A, B, r2, &U, &F);
63  *alch_vdw_energy = myVdwLambda*(U + dU);
64  *alch_vdw_force = myVdwLambda*F;
65  } else { // r2 > switchdist
66  // normal LJ potential with linear coupling
67  vdw_fswitch_forceandenergy(A, B, r2, switchdist2, cutoff2, &U, &F);
68  *alch_vdw_energy = myVdwLambda*U;
69  *alch_vdw_force = myVdwLambda*F;
70  }
71  } else {
72  // potential switching - also correct for no switching
73  if (r2 <= Rmin2) {
74  // normal LJ potential w/shift
75  const BigReal epsilon = B*B/(4.0*A);
76  vdw_forceandenergy(A, B, r2, &U, &F);
77  *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon;
78  *alch_vdw_force = F;
79  } else { // r2 > Rmin2
80  // normal LJ potential with linear coupling
81  vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
82  &switchmul2);
83  vdw_forceandenergy(A, B, r2, &U, &F);
84  *alch_vdw_energy = myVdwLambda*switchmul*U;
85  *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
86  }
87  }
88  }
89  if (myRepLambda2 < 1.0) {
90  // Same as above, but energy only. This is done separately bc the
91  // cutoffs are lambda dependent.
92  const BigReal WCAshift_2 = Rmin2*(1 - myRepLambda2)*(1 - myRepLambda2);
93  if (r2 <= Rmin2 - WCAshift_2) {
94  const BigReal epsilon = B*B/(4.0*A);
95  vdw_energy(A, B, r2 + WCAshift_2, &U_2);
96  *alch_vdw_energy_2 = U_2 + epsilon;
97  } else {
98  *alch_vdw_energy_2 = 0.0;
99  }
100  } else { // myRepLambda2 == 1.0
101  if (vdwForceSwitching) {
102  // force and potential switching
103  if (r2 <= Rmin2) {
104  // normal LJ, potential w/two shifts
105  const BigReal epsilon = B*B/(4.0*A);
106  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU_2);
107  vdw_energy(A, B, r2, &U_2);
108  *alch_vdw_energy_2 = \
109  U_2 + (1 - myVdwLambda2)*epsilon + myVdwLambda2*dU_2;
110  } else if (r2 <= switchdist2) {
111  // normal LJ potential w/shift
112  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU_2);
113  vdw_energy(A, B, r2, &U_2);
114  *alch_vdw_energy_2 = myVdwLambda2*(U_2 + dU_2);
115  } else { // r2 > switchdist
116  vdw_fswitch_energy(A, B, r2, switchdist2, cutoff2, &U_2);
117  *alch_vdw_energy_2 = myVdwLambda2*U_2;
118  }
119  } else {
120  // potential switching - also correct for no switching
121  if (r2 <= Rmin2) {
122  // normal LJ potential w/shift
123  const BigReal epsilon = B*B/(4.0*A);
124  vdw_energy(A, B, r2, &U_2);
125  *alch_vdw_energy_2 = U_2 + (1 - myVdwLambda2)*epsilon;
126  } else { // r2 > Rmin2
127  // normal LJ potential with linear coupling
128  vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
129  &switchmul2);
130  vdw_energy(A, B, r2, &U_2);
131  *alch_vdw_energy_2 = myVdwLambda2*switchmul*U_2;
132  }
133  }
134  }
135  } else { // WCA-off
136  if (vdwForceSwitching) {
137  // force and potential switching
138  if (r2 <= switchdist2) {
139  // normal LJ potential w/shift
140  vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
141  vdw_fswitch_shift(A, B, switchdist2 + myVdwShift, cutoff2, &dU);
142  vdw_energy(A, B, r2 + myVdwShift2, &U_2);
143  vdw_fswitch_shift(A, B, switchdist2 + myVdwShift2, cutoff2, &dU_2);
144  *alch_vdw_energy = myVdwLambda*(U + dU);
145  *alch_vdw_energy_2 = myVdwLambda2*(U_2 + dU_2);
146  *alch_vdw_force = myVdwLambda*F;
147  } else { // r2 > switchdist2
148  vdw_fswitch_forceandenergy(A, B, r2 + myVdwShift, \
149  switchdist2 + myVdwShift, cutoff2, &U, &F);
150  vdw_fswitch_energy(A, B, r2 + myVdwShift2, switchdist2 + myVdwShift2, \
151  cutoff2, &U_2);
152  *alch_vdw_energy = myVdwLambda*U;
153  *alch_vdw_energy_2 = myVdwLambda2*U_2;
154  *alch_vdw_force = myVdwLambda*F;
155  }
156  } else {
157  // potential switching - also correct for no switching
158  vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
159  &switchmul2);
160  vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
161  vdw_energy(A, B, r2 + myVdwShift2, &U_2);
162  *alch_vdw_energy = myVdwLambda*switchmul*U;
163  *alch_vdw_energy_2 = myVdwLambda2*switchmul*U_2;
164  *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
165  }
166  }
167 }
168 
169 #define FEPFLAG
170 #define CALCENERGY
171 
172 #define NBTYPE NBPAIR
173 #include "ComputeNonbondedBase.h"
174 #define FULLELECT
175 #include "ComputeNonbondedBase.h"
176 #define MERGEELECT
177 #include "ComputeNonbondedBase.h"
178 #undef MERGEELECT
179 #define SLOWONLY
180 #include "ComputeNonbondedBase.h"
181 #undef SLOWONLY
182 #undef FULLELECT
183 #undef NBTYPE
184 
185 #define NBTYPE NBSELF
186 #include "ComputeNonbondedBase.h"
187 #define FULLELECT
188 #include "ComputeNonbondedBase.h"
189 #define MERGEELECT
190 #include "ComputeNonbondedBase.h"
191 #undef MERGEELECT
192 #define SLOWONLY
193 #include "ComputeNonbondedBase.h"
194 #undef SLOWONLY
195 #undef FULLELECT
196 #undef NBTYPE
197 
198 #undef CALCENERGY
199 #undef FEPFLAG
200 
201 
void vdw_fswitch_energy(const BigReal A, const BigReal B, const BigReal r2, const BigReal switchdist2, const BigReal cutoff2, BigReal *U)
void vdw_forceandenergy(const BigReal A, const BigReal B, const BigReal r2, BigReal *U, BigReal *F)
const BigReal A
void vdw_switch(const BigReal r2, const BigReal switchdist2, const BigReal cutoff2, const BigReal switchfactor, BigReal *switchmul, BigReal *switchmul2)
void vdw_fswitch_forceandenergy(const BigReal A, const BigReal B, const BigReal r2, const BigReal switchdist2, const BigReal cutoff2, BigReal *U, BigReal *F)
void fep_vdw_forceandenergies(BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal myVdwShift2, BigReal switchdist2, BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda, BigReal myVdwLambda2, Bool alchWCAOn, BigReal myRepLambda, BigReal myRepLambda2, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_energy_2)
void vdw_energy(const BigReal A, const BigReal B, const BigReal r2, BigReal *U)
int Bool
Definition: common.h:133
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ cudaTextureObject_t const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cutoff2
const BigReal B
void vdw_fswitch_shift(const BigReal A, const BigReal B, const BigReal switchdist2, const BigReal cutoff2, BigReal *dU)
double BigReal
Definition: common.h:114