57 for (
int ituple=0; ituple<ntuple; ++ituple ) {
66 #if defined(DEBUG_PROTOCELL) 82 CkPrintf(
"%11d %11d %11d %11d mult=%d", i, j, k, l, mult);
83 for (
int m=0; m < mult; m++) {
87 CkPrintf(
" k=%g delta=%g n=%d\n", kdelta, delta, n);
105 Vector A = cross(r12,r23);
107 Vector B = cross(r23,r34);
113 BigReal cos_phi = (A*B)*(rAinv*rBinv);
114 BigReal sin_phi = (C*B)*(rCinv*rBinv);
116 BigReal phi= -atan2(sin_phi,cos_phi);
127 for (
int mult_num=0; mult_num<multiplicity; mult_num++)
138 K += k*(1+cos(n*phi - delta));
139 K1 += -n*k*sin(n*phi - delta);
146 else if (diff >
PI) diff -=
TWOPI;
183 if (fabs(sin_phi) > 0.1)
192 dcosdA.
x = rAinv*(cos_phi*A.
x-B.
x);
193 dcosdA.
y = rAinv*(cos_phi*A.
y-B.
y);
194 dcosdA.
z = rAinv*(cos_phi*A.
z-B.
z);
196 dcosdB.
x = rBinv*(cos_phi*B.
x-A.
x);
197 dcosdB.
y = rBinv*(cos_phi*B.
y-A.
y);
198 dcosdB.
z = rBinv*(cos_phi*B.
z-A.
z);
202 f1.
x = K1*(r23.
y*dcosdA.
z - r23.
z*dcosdA.
y);
203 f1.
y = K1*(r23.
z*dcosdA.
x - r23.
x*dcosdA.
z);
204 f1.
z = K1*(r23.
x*dcosdA.
y - r23.
y*dcosdA.
x);
206 f3.
x = K1*(r23.
z*dcosdB.
y - r23.
y*dcosdB.
z);
207 f3.
y = K1*(r23.
x*dcosdB.
z - r23.
z*dcosdB.
x);
208 f3.
z = K1*(r23.
y*dcosdB.
x - r23.
x*dcosdB.
y);
210 f2.
x = K1*(r12.
z*dcosdA.
y - r12.
y*dcosdA.
z 211 + r34.
y*dcosdB.
z - r34.
z*dcosdB.
y);
212 f2.
y = K1*(r12.
x*dcosdA.
z - r12.
z*dcosdA.
x 213 + r34.
z*dcosdB.
x - r34.
x*dcosdB.
z);
214 f2.
z = K1*(r12.
y*dcosdA.
x - r12.
x*dcosdA.
y 215 + r34.
x*dcosdB.
y - r34.
y*dcosdB.
x);
229 dsindC.
x = rCinv*(sin_phi*C.
x-B.
x);
230 dsindC.
y = rCinv*(sin_phi*C.
y-B.
y);
231 dsindC.
z = rCinv*(sin_phi*C.
z-B.
z);
233 dsindB.
x = rBinv*(sin_phi*B.
x-C.
x);
234 dsindB.
y = rBinv*(sin_phi*B.
y-C.
y);
235 dsindB.
z = rBinv*(sin_phi*B.
z-C.
z);
239 f1.
x = K1*((r23.
y*r23.
y + r23.
z*r23.
z)*dsindC.
x 240 - r23.
x*r23.
y*dsindC.
y 241 - r23.
x*r23.
z*dsindC.
z);
242 f1.
y = K1*((r23.
z*r23.
z + r23.
x*r23.
x)*dsindC.y
243 - r23.
y*r23.
z*dsindC.z
244 - r23.
y*r23.
x*dsindC.x);
245 f1.
z = K1*((r23.
x*r23.
x + r23.
y*r23.
y)*dsindC.z
246 - r23.
z*r23.
x*dsindC.x
247 - r23.
z*r23.
y*dsindC.y);
249 f3 = cross(K1,dsindB,r23);
251 f2.
x = K1*(-(r23.
y*r12.
y + r23.
z*r12.
z)*dsindC.x
252 +(2.0*r23.
x*r12.
y - r12.
x*r23.
y)*dsindC.y
253 +(2.0*r23.
x*r12.
z - r12.
x*r23.
z)*dsindC.z
254 +dsindB.z*r34.
y - dsindB.y*r34.
z);
255 f2.y = K1*(-(r23.
z*r12.
z + r23.
x*r12.
x)*dsindC.y
256 +(2.0*r23.
y*r12.
z - r12.
y*r23.
z)*dsindC.z
257 +(2.0*r23.
y*r12.
x - r12.
y*r23.
x)*dsindC.x
258 +dsindB.x*r34.
z - dsindB.z*r34.
x);
259 f2.z = K1*(-(r23.
x*r12.
x + r23.
y*r12.
y)*dsindC.z
260 +(2.0*r23.
z*r12.
x - r12.
z*r23.
x)*dsindC.x
261 +(2.0*r23.
z*r12.
y - r12.
z*r23.
y)*dsindC.y
262 +dsindB.y*r34.
x - dsindB.x*r34.
y);
306 DebugM(3,
"::computeForce() -- ending with delta energy " << K << std::endl);
308 reduction[virialIndex_XX] += ( f1.x * r12.
x + f2.x * r23.
x + f3.x * r34.
x );
309 reduction[virialIndex_XY] += ( f1.x * r12.
y + f2.x * r23.
y + f3.x * r34.
y );
310 reduction[virialIndex_XZ] += ( f1.x * r12.
z + f2.x * r23.
z + f3.x * r34.
z );
311 reduction[virialIndex_YX] += ( f1.y * r12.
x + f2.y * r23.
x + f3.y * r34.
x );
312 reduction[virialIndex_YY] += ( f1.y * r12.
y + f2.y * r23.
y + f3.y * r34.
y );
313 reduction[virialIndex_YZ] += ( f1.y * r12.
z + f2.y * r23.
z + f3.y * r34.
z );
314 reduction[virialIndex_ZX] += ( f1.z * r12.
x + f2.z * r23.
x + f3.z * r34.
x );
315 reduction[virialIndex_ZY] += ( f1.z * r12.
y + f2.z * r23.
y + f3.z * r34.
y );
316 reduction[virialIndex_ZZ] += ( f1.z * r12.
z + f2.z * r23.
z + f3.z * r34.
z );
318 if (pressureProfileData) {
338 f1.x * r12.
x, f1.y * r12.
y, f1.z * r12.
z,
339 pressureProfileData);
342 f2.x * r23.
x, f2.y * r23.
y, f2.z * r23.
z,
343 pressureProfileData);
346 f3.x * r34.
x, f3.y * r34.
y, f3.z * r34.
z,
347 pressureProfileData);
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
SimParameters * simParameters
static BigReal pressureProfileMin
Molecule stores the structural information for the system.
FourBodyConsts values[MAX_MULTIPLICITY]
static int pressureProfileAtomTypes
void pp_clamp(int &n, int nslabs)
static BigReal pressureProfileThickness
const DihedralValue * value
int get_fep_bonded_type(const int *atomID, unsigned int order) const
static int pressureProfileSlabs
NAMD_HOST_DEVICE BigReal rlength(void) const
NAMD_HOST_DEVICE Vector delta(const Position &pos1, const Position &pos2) const