NAMD
ComputeNonbondedTI.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 dU/dl for TI */
16 /* ********************************** */
18  BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2,
19  BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda,
20  BigReal alchVdwShiftCoeff, Bool alchWCAOn, BigReal myRepLambda,
21  BigReal* alch_vdw_energy, BigReal* alch_vdw_force,
22  BigReal* alch_vdw_dUdl) {
23  BigReal U, F, dU, 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  /* THIS BRANCH IS PARTIALLY INCOMPLETE AND BLOCKED INSIDE SimParameters
32  *
33  * The problem is that WCA creates TWO energy terms and thus two different
34  * TI derivatives, but the reduction code is currently only set up for
35  * one derivative. There is no such problem with FEP because the energy
36  * is not decomposed. In principle, the contribution to the free energy
37  * from each derivative is zero while the other is changing (e.g. while
38  * repulsion is changing dispersion is not), but the derivatives are still
39  * non-zero and it is NAMD convention to report them even as such.
40  * However, even if we were willing to break this convention, both
41  * derivatives still compute to the free energy at the boundary where
42  * repulsion is completely on and dispersion is exactly zero - two
43  * derivatives are thus required for a complete working implementation.
44  *
45  */
46 
47  // WCA-on, auxilliary, lambda-dependent cutoff based on Rmin
48  //
49  // Avoid divide by zero - correctly zeroes interaction below.
50  const BigReal Rmin2 = (B <= 0.0 ? 0.0 : powf(2.0*A/B, 1.f/3));
51  if (myRepLambda < 1.0) {
52  // modified repulsive soft-core
53  // All of the scaling is baked into the shift and cutoff values. There
54  // are no linear coupling terms out front.
55  const BigReal WCAshift = Rmin2*(1 - myRepLambda)*(1 - myRepLambda);
56  if (r2 <= Rmin2 - WCAshift) {
57  const BigReal epsilon = B*B/(4.0*A);
58  vdw_forceandenergy(A, B, r2 + WCAshift, &U, &F);
59  *alch_vdw_energy = U + epsilon;
60  *alch_vdw_force = F;
61  *alch_vdw_dUdl = Rmin2*(1 - myRepLambda)*F;
62  } else {
63  *alch_vdw_energy = 0.0;
64  *alch_vdw_force = 0.0;
65  *alch_vdw_dUdl = 0.0;
66  }
67  } else { // myRepLambda == 1.0
68  if (vdwForceSwitching) {
69  // force and potential switching
70  if (r2 <= Rmin2) {
71  // normal LJ, potential w/two shifts
72  const BigReal epsilon = B*B/(4.0*A);
73  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
74  vdw_forceandenergy(A, B, r2, &U, &F);
75  *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon + myVdwLambda*dU;
76  *alch_vdw_force = F;
77  *alch_vdw_dUdl = dU - epsilon;
78  } else if (r2 <= switchdist2) {
79  // normal LJ potential w/shift
80  vdw_fswitch_shift(A, B, switchdist2, cutoff2, &dU);
81  vdw_forceandenergy(A, B, r2, &U, &F);
82  *alch_vdw_energy = myVdwLambda*(U + dU);
83  *alch_vdw_force = myVdwLambda*F;
84  *alch_vdw_dUdl = U + dU;
85  } else { // r2 > switchdist
86  // normal LJ potential with linear coupling
87  vdw_fswitch_forceandenergy(A, B, r2, switchdist2, cutoff2, &U, &F);
88  *alch_vdw_energy = myVdwLambda*U;
89  *alch_vdw_force = myVdwLambda*F;
90  *alch_vdw_dUdl = U;
91  }
92  } else {
93  // potential switching - also correct for no switching
94  if (r2 <= Rmin2) {
95  // normal LJ potential w/shift
96  const BigReal epsilon = B*B/(4.0*A);
97  vdw_forceandenergy(A, B, r2, &U, &F);
98  *alch_vdw_energy = U + (1 - myVdwLambda)*epsilon;
99  *alch_vdw_force = F;
100  *alch_vdw_dUdl = -epsilon;
101  } else { // r2 > Rmin2
102  // normal LJ potential with linear coupling
103  vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
104  &switchmul2);
105  vdw_forceandenergy(A, B, r2, &U, &F);
106  *alch_vdw_energy = myVdwLambda*switchmul*U;
107  *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
108  *alch_vdw_dUdl = switchmul*U;
109  }
110  }
111  }
112  } else { // WCA-off
113  if (vdwForceSwitching) {
114  // force and potential switching
115  if (r2 <= switchdist2) {
116  // normal LJ potential w/shift
117  vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
118  vdw_fswitch_shift(A, B, switchdist2 + myVdwShift, cutoff2, &dU);
119  *alch_vdw_energy = myVdwLambda*(U + dU);
120  *alch_vdw_force = myVdwLambda*F;
121  *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F + dU;
122  } else { // r2 > switchdist2
123  vdw_fswitch_forceandenergy(A, B, r2 + myVdwShift, \
124  switchdist2 + myVdwShift, cutoff2, &U, &F);
125  *alch_vdw_energy = myVdwLambda*U;
126  *alch_vdw_force = myVdwLambda*F;
127  *alch_vdw_dUdl = U + 0.5*myVdwLambda*alchVdwShiftCoeff*F;
128  }
129  } else {
130  // potential switching - also correct for no switching
131  vdw_switch(r2, switchdist2, cutoff2, switchfactor, &switchmul, \
132  &switchmul2);
133  vdw_forceandenergy(A, B, r2 + myVdwShift, &U, &F);
134  *alch_vdw_energy = myVdwLambda*switchmul*U;
135  *alch_vdw_force = myVdwLambda*(switchmul*F + switchmul2*U);
136  *alch_vdw_dUdl = switchmul*(U + 0.5*myVdwLambda*alchVdwShiftCoeff*F);
137  }
138  }
139 }
140 
141 
142 /*************THERMODYNAMIC INTEGRATION*************/
143 #define TIFLAG
144 #define CALCENERGY
145 
146 #define NBTYPE NBPAIR
147 #include "ComputeNonbondedBase.h"
148 #define FULLELECT
149 #include "ComputeNonbondedBase.h"
150 #define MERGEELECT
151 #include "ComputeNonbondedBase.h"
152 #undef MERGEELECT
153 #define SLOWONLY
154 #include "ComputeNonbondedBase.h"
155 #undef SLOWONLY
156 #undef FULLELECT
157 #undef NBTYPE
158 
159 #define NBTYPE NBSELF
160 #include "ComputeNonbondedBase.h"
161 #define FULLELECT
162 #include "ComputeNonbondedBase.h"
163 #define MERGEELECT
164 #include "ComputeNonbondedBase.h"
165 #undef MERGEELECT
166 #define SLOWONLY
167 #include "ComputeNonbondedBase.h"
168 #undef SLOWONLY
169 #undef FULLELECT
170 #undef NBTYPE
171 
172 #undef CALCENERGY
173 
174 #define NBTYPE NBPAIR
175 #include "ComputeNonbondedBase.h"
176 #define FULLELECT
177 #include "ComputeNonbondedBase.h"
178 #define MERGEELECT
179 #include "ComputeNonbondedBase.h"
180 #undef MERGEELECT
181 #define SLOWONLY
182 #include "ComputeNonbondedBase.h"
183 #undef SLOWONLY
184 #undef FULLELECT
185 #undef NBTYPE
186 
187 #define NBTYPE NBSELF
188 #include "ComputeNonbondedBase.h"
189 #define FULLELECT
190 #include "ComputeNonbondedBase.h"
191 #define MERGEELECT
192 #include "ComputeNonbondedBase.h"
193 #undef MERGEELECT
194 #define SLOWONLY
195 #include "ComputeNonbondedBase.h"
196 #undef SLOWONLY
197 #undef FULLELECT
198 #undef NBTYPE
199 
200 #undef TIFLAG
201 
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)
int Bool
Definition: common.h:133
void ti_vdw_force_energy_dUdl(BigReal A, BigReal B, BigReal r2, BigReal myVdwShift, BigReal switchdist2, BigReal cutoff2, BigReal switchfactor, Bool vdwForceSwitching, BigReal myVdwLambda, BigReal alchVdwShiftCoeff, Bool alchWCAOn, BigReal myRepLambda, BigReal *alch_vdw_energy, BigReal *alch_vdw_force, BigReal *alch_vdw_dUdl)
const BigReal B
void vdw_fswitch_shift(const BigReal A, const BigReal B, const BigReal switchdist2, const BigReal cutoff2, BigReal *dU)
__global__ void const int const TileList *__restrict__ TileExcl *__restrict__ const int *__restrict__ const int const float2 *__restrict__ const int *__restrict__ const float3 const float3 const float3 const float4 *__restrict__ const float cutoff2
double BigReal
Definition: common.h:114