Definition at line 40 of file ComputeDihedrals.C.
References A, TuplePatchElem::af, SimParameters::alchOn, atomID, B, cross(), DebugM, four_body_consts::delta, Lattice::delta(), dihedralEnergyIndex, dihedralEnergyIndex_f, dihedralEnergyIndex_ti_1, dihedralEnergyIndex_ti_2, TuplePatchElem::f, Patch::flags, Molecule::get_fep_bonded_type(), SimParameters::getBondLambda(), SimParameters::getCurrentLambda(), SimParameters::getCurrentLambda2(), four_body_consts::k, Patch::lattice, localIndex, Node::molecule, DihedralValue::multiplicity, four_body_consts::n, Node::Object(), p, TuplePatchElem::p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, Vector::rlength(), scale, Node::simParameters, simParams, SimParameters::singleTopology, size, Flags::step, TWOPI, value, DihedralValue::values, TuplePatchElem::x, Vector::x, Vector::y, and Vector::z.
57 for (
int ituple=0; ituple<ntuple; ++ituple ) {
87 BigReal cos_phi = (A*
B)*(rAinv*rBinv);
88 BigReal sin_phi = (C*
B)*(rCinv*rBinv);
90 BigReal phi= -atan2(sin_phi,cos_phi);
101 for (
int mult_num=0; mult_num<multiplicity; mult_num++)
112 K += k*(1+cos(n*phi - delta));
113 K1 += -n*k*sin(n*phi - delta);
120 else if (diff >
PI) diff -=
TWOPI;
157 if (fabs(sin_phi) > 0.1)
166 dcosdA.
x = rAinv*(cos_phi*A.
x-B.
x);
167 dcosdA.
y = rAinv*(cos_phi*A.
y-B.
y);
168 dcosdA.
z = rAinv*(cos_phi*A.
z-B.
z);
170 dcosdB.
x = rBinv*(cos_phi*B.
x-A.
x);
171 dcosdB.
y = rBinv*(cos_phi*B.
y-A.
y);
172 dcosdB.
z = rBinv*(cos_phi*B.
z-A.
z);
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);
203 dsindC.
x = rCinv*(sin_phi*C.
x-B.
x);
204 dsindC.
y = rCinv*(sin_phi*C.
y-B.
y);
205 dsindC.
z = rCinv*(sin_phi*C.
z-B.
z);
207 dsindB.
x = rBinv*(sin_phi*B.
x-C.
x);
208 dsindB.
y = rBinv*(sin_phi*B.
y-C.
y);
209 dsindB.
z = rBinv*(sin_phi*B.
z-C.
z);
213 f1.
x = K1*((r23.
y*r23.
y + r23.
z*r23.
z)*dsindC.
x
214 - r23.
x*r23.
y*dsindC.
y
215 - r23.
x*r23.
z*dsindC.
z);
216 f1.
y = K1*((r23.
z*r23.
z + r23.
x*r23.
x)*dsindC.y
217 - r23.
y*r23.
z*dsindC.z
218 - r23.
y*r23.
x*dsindC.x);
219 f1.
z = K1*((r23.
x*r23.
x + r23.
y*r23.
y)*dsindC.z
220 - r23.
z*r23.
x*dsindC.x
221 - r23.
z*r23.
y*dsindC.y);
223 f3 =
cross(K1,dsindB,r23);
225 f2.
x = K1*(-(r23.
y*r12.
y + r23.
z*r12.
z)*dsindC.x
226 +(2.0*r23.
x*r12.
y - r12.
x*r23.
y)*dsindC.y
227 +(2.0*r23.
x*r12.
z - r12.
x*r23.
z)*dsindC.z
228 +dsindB.z*r34.
y - dsindB.y*r34.
z);
229 f2.y = K1*(-(r23.
z*r12.
z + r23.
x*r12.
x)*dsindC.y
230 +(2.0*r23.
y*r12.
z - r12.
y*r23.
z)*dsindC.z
231 +(2.0*r23.
y*r12.
x - r12.
y*r23.
x)*dsindC.x
232 +dsindB.x*r34.
z - dsindB.z*r34.
x);
233 f2.z = K1*(-(r23.
x*r12.
x + r23.
y*r12.
y)*dsindC.z
234 +(2.0*r23.
z*r12.
x - r12.
z*r23.
x)*dsindC.x
235 +(2.0*r23.
z*r12.
y - r12.
z*r23.
y)*dsindC.y
236 +dsindB.y*r34.
x - dsindB.x*r34.
y);
280 DebugM(3,
"::computeForce() -- ending with delta energy " << K << std::endl);
282 reduction[virialIndex_XX] += ( f1.x * r12.
x + f2.x * r23.
x + f3.x * r34.
x );
283 reduction[virialIndex_XY] += ( f1.x * r12.
y + f2.x * r23.
y + f3.x * r34.
y );
284 reduction[virialIndex_XZ] += ( f1.x * r12.
z + f2.x * r23.
z + f3.x * r34.
z );
285 reduction[virialIndex_YX] += ( f1.y * r12.
x + f2.y * r23.
x + f3.y * r34.
x );
286 reduction[virialIndex_YY] += ( f1.y * r12.
y + f2.y * r23.
y + f3.y * r34.
y );
287 reduction[virialIndex_YZ] += ( f1.y * r12.
z + f2.y * r23.
z + f3.y * r34.
z );
288 reduction[virialIndex_ZX] += ( f1.z * r12.
x + f2.z * r23.
x + f3.z * r34.
x );
289 reduction[virialIndex_ZY] += ( f1.z * r12.
y + f2.z * r23.
y + f3.z * r34.
y );
290 reduction[virialIndex_ZZ] += ( f1.z * r12.
z + f2.z * r23.
z + f3.z * r34.
z );
292 if (pressureProfileData) {
312 f1.x * r12.
x, f1.y * r12.
y, f1.z * r12.
z,
313 pressureProfileData);
316 f2.x * r23.
x, f2.y * r23.
y, f2.z * r23.
z,
317 pressureProfileData);
320 f3.x * r34.
x, f3.y * r34.
y, f3.z * r34.
z,
321 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
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
int get_fep_bonded_type(const int *atomID, unsigned int order) const
static BigReal pressureProfileMin
FourBodyConsts values[MAX_MULTIPLICITY]
static int pressureProfileAtomTypes
void pp_clamp(int &n, int nslabs)
static BigReal pressureProfileThickness
const DihedralValue * value
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
BigReal getCurrentLambda(const int)
static int pressureProfileSlabs