27 #ifdef MEM_OPT_VERSION
28 NAMD_die(
"Should not be called in ImproperElem::getMoleculePointers in memory optimized version!");
31 *byatom = mol->impropersByAtom;
32 *structarray = mol->impropers;
57 for (
int ituple=0; ituple<ntuple; ++ituple ) {
83 Force f1(0,0,0),f2(0,0,0),f3(0,0,0);
110 cos_phi = A*B/(rA*rB);
111 sin_phi = C*B/(rC*rB);
117 phi= -atan2(sin_phi,cos_phi);
119 if (fabs(sin_phi) > 0.1)
124 dcosdA = rA*(cos_phi*A-
B);
125 dcosdB = rB*(cos_phi*B-
A);
132 dsindC = rC*(sin_phi*C-
B);
133 dsindB = rB*(sin_phi*B-C);
139 for (
int mult_num=0; mult_num<multiplicity; mult_num++)
150 K = k*(1+cos(n*phi - delta));
151 K1 = -n*k*sin(n*phi - delta);
158 else if (diff >
PI) diff -=
TWOPI;
171 if (fabs(sin_phi) > 0.1)
176 f1.x += K1*(r23.
y*dcosdA.
z - r23.
z*dcosdA.
y);
177 f1.y += K1*(r23.
z*dcosdA.
x - r23.
x*dcosdA.
z);
178 f1.z += K1*(r23.
x*dcosdA.
y - r23.
y*dcosdA.
x);
180 f3.
x += K1*(r23.
z*dcosdB.
y - r23.
y*dcosdB.
z);
181 f3.
y += K1*(r23.
x*dcosdB.
z - r23.
z*dcosdB.
x);
182 f3.
z += K1*(r23.
y*dcosdB.
x - r23.
x*dcosdB.
y);
184 f2.x += K1*(r12.
z*dcosdA.
y - r12.
y*dcosdA.
z
185 + r34.
y*dcosdB.
z - r34.
z*dcosdB.
y);
186 f2.y += K1*(r12.
x*dcosdA.
z - r12.
z*dcosdA.
x
187 + r34.
z*dcosdB.
x - r34.
x*dcosdB.
z);
188 f2.z += K1*(r12.
y*dcosdA.
x - r12.
x*dcosdA.
y
189 + r34.
x*dcosdB.
y - r34.
y*dcosdB.
x);
197 f1.x += K1*((r23.
y*r23.
y + r23.
z*r23.
z)*dsindC.
x
198 - r23.
x*r23.
y*dsindC.
y
199 - r23.
x*r23.
z*dsindC.
z);
200 f1.y += K1*((r23.
z*r23.
z + r23.
x*r23.
x)*dsindC.
y
201 - r23.
y*r23.
z*dsindC.
z
202 - r23.
y*r23.
x*dsindC.
x);
203 f1.z += K1*((r23.
x*r23.
x + r23.
y*r23.
y)*dsindC.
z
204 - r23.
z*r23.
x*dsindC.
x
205 - r23.
z*r23.
y*dsindC.
y);
207 f3 +=
cross(K1,dsindB,r23);
209 f2.x += K1*(-(r23.
y*r12.
y + r23.
z*r12.
z)*dsindC.
x
210 +(2.0*r23.
x*r12.
y - r12.
x*r23.
y)*dsindC.
y
211 +(2.0*r23.
x*r12.
z - r12.
x*r23.
z)*dsindC.
z
212 +dsindB.
z*r34.
y - dsindB.
y*r34.
z);
213 f2.y += K1*(-(r23.
z*r12.
z + r23.
x*r12.
x)*dsindC.
y
214 +(2.0*r23.
y*r12.
z - r12.
y*r23.
z)*dsindC.
z
215 +(2.0*r23.
y*r12.
x - r12.
y*r23.
x)*dsindC.
x
216 +dsindB.
x*r34.
z - dsindB.
z*r34.
x);
217 f2.z += K1*(-(r23.
x*r12.
x + r23.
y*r12.
y)*dsindC.
z
218 +(2.0*r23.
z*r12.
x - r12.
z*r23.
x)*dsindC.
x
219 +(2.0*r23.
z*r12.
y - r12.
z*r23.
y)*dsindC.
y
220 +dsindB.
y*r34.
x - dsindB.
x*r34.
y);
232 energy *= bond_lambda_1;
241 energy *= bond_lambda_2;
256 DebugM(3,
"::computeForce() -- ending with delta energy " << energy << std::endl);
258 reduction[virialIndex_XX] += ( f1.x * r12.
x + f2.x * r23.
x + f3.
x * r34.
x );
259 reduction[virialIndex_XY] += ( f1.x * r12.
y + f2.x * r23.
y + f3.
x * r34.
y );
260 reduction[virialIndex_XZ] += ( f1.x * r12.
z + f2.x * r23.
z + f3.
x * r34.
z );
261 reduction[virialIndex_YX] += ( f1.y * r12.
x + f2.y * r23.
x + f3.
y * r34.
x );
262 reduction[virialIndex_YY] += ( f1.y * r12.
y + f2.y * r23.
y + f3.
y * r34.
y );
263 reduction[virialIndex_YZ] += ( f1.y * r12.
z + f2.y * r23.
z + f3.
y * r34.
z );
264 reduction[virialIndex_ZX] += ( f1.z * r12.
x + f2.z * r23.
x + f3.
z * r34.
x );
265 reduction[virialIndex_ZY] += ( f1.z * r12.
y + f2.z * r23.
y + f3.
z * r34.
y );
266 reduction[virialIndex_ZZ] += ( f1.z * r12.
z + f2.z * r23.
z + f3.
z * r34.
z );
268 if (pressureProfileData) {
288 f1.x * r12.
x, f1.y * r12.
y, f1.z * r12.
z,
289 pressureProfileData);
292 f2.x * r23.
x, f2.y * r23.
y, f2.z * r23.
z,
293 pressureProfileData);
296 f3.
x * r34.
x, f3.
y * r34.
y, f3.
z * r34.
z,
297 pressureProfileData);
309 ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
static int pressureProfileAtomTypes
static void submitReductionData(BigReal *, SubmitReduction *)
static void getMoleculePointers(Molecule *, int *, int32 ***, Improper **)
const ImproperValue * value
#define ADD_TENSOR(R, RL, D, DL)
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
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
int get_fep_bonded_type(const int *atomID, unsigned int order) const
BigReal length(void) const
static void getParameterPointers(Parameters *, const ImproperValue **)
static BigReal pressureProfileMin
void pp_clamp(int &n, int nslabs)
static int pressureProfileSlabs
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
ImproperValue * improper_array
void NAMD_die(const char *err_msg)
FourBodyConsts values[MAX_MULTIPLICITY]
static void computeForce(ImproperElem *, int, BigReal *, BigReal *)
BigReal getCurrentLambda(const int)
static BigReal pressureProfileThickness