NAMD
Public Types | Public Member Functions | Static Public Member Functions | Public Attributes | Static Public Attributes | List of all members
CrosstermElem Class Reference

#include <ComputeCrossterms.h>

Public Types

enum  { size = 8 }
 
enum  {
  crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2,
  TENSOR =(virialIndex), reductionDataSize
}
 
enum  { reductionChecksumLabel = REDUCTION_CROSSTERM_CHECKSUM }
 

Public Member Functions

int hash () const
 
 CrosstermElem ()
 
 CrosstermElem (AtomID atom0, const TupleSignature *sig, const CrosstermValue *v)
 
 CrosstermElem (const Crossterm *a, const CrosstermValue *v)
 
 CrosstermElem (AtomID atom0, AtomID atom1, AtomID atom2, AtomID atom3, AtomID atom4, AtomID atom5, AtomID atom6, AtomID atom7)
 
 ~CrosstermElem ()
 
int operator== (const CrosstermElem &a) const
 
int operator< (const CrosstermElem &a) const
 

Static Public Member Functions

static void computeForce (CrosstermElem *, int, BigReal *, BigReal *)
 
static void getMoleculePointers (Molecule *, int *, int32 ***, Crossterm **)
 
static void getParameterPointers (Parameters *, const CrosstermValue **)
 
static void getTupleInfo (AtomSignature *sig, int *count, TupleSignature **t)
 
static void submitReductionData (BigReal *, SubmitReduction *)
 

Public Attributes

AtomID atomID [size]
 
int localIndex [size]
 
TuplePatchElem * p [size]
 
Real scale
 
const CrosstermValue * value
 

Static Public Attributes

static int pressureProfileSlabs = 0
 
static int pressureProfileAtomTypes = 1
 
static BigReal pressureProfileThickness = 0
 
static BigReal pressureProfileMin = 0
 

Detailed Description

Definition at line 17 of file ComputeCrossterms.h.

Member Enumeration Documentation

◆ anonymous enum

anonymous enum
Enumerator
size 

Definition at line 20 of file ComputeCrossterms.h.

20 { size = 8 };

◆ anonymous enum

anonymous enum
Enumerator
crosstermEnergyIndex 
crosstermEnergyIndex_f 
crosstermEnergyIndex_ti_1 
crosstermEnergyIndex_ti_2 
TENSOR 
reductionDataSize 

Definition at line 46 of file ComputeCrossterms.h.

46  { crosstermEnergyIndex, crosstermEnergyIndex_f,
47  crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2,
48  TENSOR(virialIndex), reductionDataSize };

◆ anonymous enum

anonymous enum
Enumerator
reductionChecksumLabel 

Definition at line 49 of file ComputeCrossterms.h.

49 { reductionChecksumLabel = REDUCTION_CROSSTERM_CHECKSUM };

Constructor & Destructor Documentation

◆ CrosstermElem() [1/4]

CrosstermElem::CrosstermElem ( )
inline

Definition at line 52 of file ComputeCrossterms.h.

52 { ; }

◆ CrosstermElem() [2/4]

CrosstermElem::CrosstermElem ( AtomID  atom0,
const TupleSignature *  sig,
const CrosstermValue *  v 
)
inline

Definition at line 54 of file ComputeCrossterms.h.

References atomID, TupleSignature::offset, TupleSignature::tupleParamType, and value.

54  {
55  atomID[0] = atom0;
56  atomID[1] = atom0 + sig->offset[0];
57  atomID[2] = atom0 + sig->offset[1];
58  atomID[3] = atom0 + sig->offset[2];
59  atomID[4] = atom0 + sig->offset[3];
60  atomID[5] = atom0 + sig->offset[4];
61  atomID[6] = atom0 + sig->offset[5];
62  atomID[7] = atom0 + sig->offset[6];
63  value = &v[sig->tupleParamType];
64 
65  }
CrosstermElem::value
const CrosstermValue * value
Definition: ComputeCrossterms.h:41
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21
TupleSignature::tupleParamType
Index tupleParamType
Definition: structures.h:204

◆ CrosstermElem() [3/4]

CrosstermElem::CrosstermElem ( const Crossterm *  a,
const CrosstermValue *  v 
)
inline

Definition at line 67 of file ComputeCrossterms.h.

References crossterm::atom1, crossterm::atom2, crossterm::atom3, crossterm::atom4, crossterm::atom5, crossterm::atom6, crossterm::atom7, crossterm::atom8, atomID, crossterm::crossterm_type, and value.

67  {
68  atomID[0] = a->atom1;
69  atomID[1] = a->atom2;
70  atomID[2] = a->atom3;
71  atomID[3] = a->atom4;
72  atomID[4] = a->atom5;
73  atomID[5] = a->atom6;
74  atomID[6] = a->atom7;
75  atomID[7] = a->atom8;
76  value = &v[a->crossterm_type];
77  }
crossterm::crossterm_type
Index crossterm_type
Definition: structures.h:91
CrosstermElem::value
const CrosstermValue * value
Definition: ComputeCrossterms.h:41
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21
crossterm::atom5
int32 atom5
Definition: structures.h:87
crossterm::atom8
int32 atom8
Definition: structures.h:90
crossterm::atom1
int32 atom1
Definition: structures.h:83
crossterm::atom4
int32 atom4
Definition: structures.h:86
crossterm::atom3
int32 atom3
Definition: structures.h:85
crossterm::atom2
int32 atom2
Definition: structures.h:84
crossterm::atom7
int32 atom7
Definition: structures.h:89
crossterm::atom6
int32 atom6
Definition: structures.h:88

◆ CrosstermElem() [4/4]

CrosstermElem::CrosstermElem ( AtomID  atom0,
AtomID  atom1,
AtomID  atom2,
AtomID  atom3,
AtomID  atom4,
AtomID  atom5,
AtomID  atom6,
AtomID  atom7 
)
inline

Definition at line 79 of file ComputeCrossterms.h.

References atomID.

80  {
81  atomID[0] = atom0;
82  atomID[1] = atom1;
83  atomID[2] = atom2;
84  atomID[3] = atom3;
85  atomID[4] = atom4;
86  atomID[5] = atom5;
87  atomID[6] = atom6;
88  atomID[7] = atom7;
89  }
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21

◆ ~CrosstermElem()

CrosstermElem::~CrosstermElem ( )
inline

Definition at line 90 of file ComputeCrossterms.h.

90 {};

Member Function Documentation

◆ computeForce()

void CrosstermElem::computeForce ( CrosstermElem *  tuples,
int  ntuple,
BigReal *  reduction,
BigReal *  pressureProfileData 
)
static

Definition at line 53 of file ComputeCrossterms.C.

References TuplePatchElem::af, atomID, CrosstermValue::c, CMAP_PHI_0, CMAP_PSI_0, CMAP_SETUP_DIM, CMAP_SPACING, CMAP_TABLE_DIM, crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2, CrosstermData::d00, CrosstermData::d01, CrosstermData::d10, CrosstermData::d11, DebugM, Lattice::delta(), TuplePatchElem::f, Patch::flags, Molecule::get_fep_type(), INDEX, Patch::lattice, Vector::length(), localIndex, Node::molecule, Node::Object(), order, p, TuplePatchElem::p, CompAtom::partition, PI, CompAtom::position, pp_clamp(), pp_reduction(), pressureProfileAtomTypes, pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, scale, Node::simParameters, simParams, size, Flags::step, value, TuplePatchElem::x, Vector::x, Vector::y, and Vector::z.

55 {
56  const Lattice & lattice = tuples[0].p[0]->p->lattice;
57  //fepb BKR
58  SimParameters *const simParams = Node::Object()->simParameters;
59  const int step = tuples[0].p[0]->p->flags.step;
60  const BigReal alchLambda = simParams->getCurrentLambda(step);
61  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
62  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
63  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
64  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
65  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
66  Molecule *const mol = Node::Object()->molecule;
67  //fepe
68 
69  for ( int ituple=0; ituple<ntuple; ++ituple ) {
70  const CrosstermElem &tup = tuples[ituple];
71  enum { size = 8 };
72  const AtomID (&atomID)[size](tup.atomID);
73  const int (&localIndex)[size](tup.localIndex);
74  TuplePatchElem * const(&p)[size](tup.p);
75  const Real (&scale)(tup.scale);
76  const CrosstermValue * const(&value)(tup.value);
77 
78  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
79  << localIndex[1] << " " << localIndex[2] << " " <<
80  localIndex[3] << std::endl);
81 
82  // Vector r12, r23, r34; // vector between atoms
83  Vector A,B,C,D,E,F; // cross products
84  BigReal rA,rB,rC,rD,rE,rF; // length of vectors
85  BigReal energy=0; // energy from the angle
86  BigReal phi,psi; // angle between the plans
87  double cos_phi,cos_psi; // cos(phi)
88  double sin_phi,sin_psi; // sin(phi)
89  Vector dcosdA,dcosdD; // Derivative d(cos(phi))/dA
90  Vector dcosdB,dcosdE; // Derivative d(cos(phi))/dB
91  Vector dsindC,dsindF; // Derivative d(sin(phi))/dC
92  Vector dsindB,dsindE; // Derivative d(sin(phi))/dB
93  BigReal U,U_phi,U_psi; // energy constants
94  BigReal diff; // for periodicity
95  Force f1(0,0,0),f2(0,0,0),f3(0,0,0); // force components
96  Force f4(0,0,0),f5(0,0,0),f6(0,0,0); // force components
97 
98  //DebugM(3, "::computeForce() -- starting with crossterm type " << crosstermType << std::endl);
99 
100  // Calculate the vectors between atoms
101  const Position & pos0 = p[0]->x[localIndex[0]].position;
102  const Position & pos1 = p[1]->x[localIndex[1]].position;
103  const Position & pos2 = p[2]->x[localIndex[2]].position;
104  const Position & pos3 = p[3]->x[localIndex[3]].position;
105  const Position & pos4 = p[4]->x[localIndex[4]].position;
106  const Position & pos5 = p[5]->x[localIndex[5]].position;
107  const Position & pos6 = p[6]->x[localIndex[6]].position;
108  const Position & pos7 = p[7]->x[localIndex[7]].position;
109  const Vector r12 = lattice.delta(pos0,pos1);
110  const Vector r23 = lattice.delta(pos1,pos2);
111  const Vector r34 = lattice.delta(pos2,pos3);
112  const Vector r56 = lattice.delta(pos4,pos5);
113  const Vector r67 = lattice.delta(pos5,pos6);
114  const Vector r78 = lattice.delta(pos6,pos7);
115 
116  // Calculate the cross products
117  A = cross(r12,r23);
118  B = cross(r23,r34);
119  C = cross(r23,A);
120  D = cross(r56,r67);
121  E = cross(r67,r78);
122  F = cross(r67,D);
123 
124  // Calculate the distances
125  rA = A.length();
126  rB = B.length();
127  rC = C.length();
128  rD = D.length();
129  rE = E.length();
130  rF = F.length();
131 
132  // Calculate the sin and cos
133  cos_phi = A*B/(rA*rB);
134  sin_phi = C*B/(rC*rB);
135  cos_psi = D*E/(rD*rE);
136  sin_psi = F*E/(rF*rE);
137 
138  // Normalize B
139  rB = 1.0/rB;
140  B *= rB;
141  rE = 1.0/rE;
142  E *= rE;
143 
144  phi= -atan2(sin_phi,cos_phi);
145  psi= -atan2(sin_psi,cos_psi);
146 
147  if (fabs(sin_phi) > 0.1) {
148  // Normalize A
149  rA = 1.0/rA;
150  A *= rA;
151  dcosdA = rA*(cos_phi*A-B);
152  dcosdB = rB*(cos_phi*B-A);
153  } else {
154  // Normalize C
155  rC = 1.0/rC;
156  C *= rC;
157  dsindC = rC*(sin_phi*C-B);
158  dsindB = rB*(sin_phi*B-C);
159  }
160 
161  if (fabs(sin_psi) > 0.1) {
162  // Normalize A
163  rD = 1.0/rD;
164  D *= rD;
165  dcosdD = rD*(cos_psi*D-E);
166  dcosdE = rE*(cos_psi*E-D);
167  } else {
168  // Normalize C
169  rF = 1.0/rF;
170  F *= rF;
171  dsindF = rF*(sin_psi*F-E);
172  dsindE = rE*(sin_psi*E-F);
173  }
174 
175  // Calculate the energy
176 {
177  const double h = CMAP_SPACING * PI / 180.0;
178  const double h_1 = 1.0 / h;
179 #ifdef FASTER
180  const double six_h = 6.0 * h_1;
181 #endif
182  const double phi_0 = CMAP_PHI_0 * PI / 180.0;
183  const double psi_0 = CMAP_PSI_0 * PI / 180.0;
184 
185  enum { D = CMAP_TABLE_DIM };
186  double xa[2], xb[2], dxa[2], dxb[2];
187  double ya[2], yb[2], dya[2], dyb[2];
188  double t, dx_h, dy_h;
189 #ifdef FASTER
190  double s1, s2, s3, s4, s5;
191 #endif
192  double f = 0, fx = 0, fy = 0;
193  const CrosstermData *table = &value->c[0][0];
194  int i, j, ilo, jlo, ij;
195 
196  /* distance measured in grid points between angle and smallest value */
197  dx_h = (phi - phi_0) * h_1;
198  dy_h = (psi - psi_0) * h_1;
199 
200  /* find smallest numbered grid point in stencil */
201  ilo = (int) floor(dx_h);
202  if ( ilo < 0 ) ilo = 0;
203  if ( ilo >= CMAP_SETUP_DIM ) ilo = CMAP_SETUP_DIM - 1;
204  jlo = (int) floor(dy_h);
205  if ( jlo < 0 ) jlo = 0;
206  if ( jlo >= CMAP_SETUP_DIM ) jlo = CMAP_SETUP_DIM - 1;
207 
208 #if !defined(FASTER)
209 
210  /* find t for x-dimension and compute xa, xb, dxa, dxb */
211  t = dx_h - (double) ilo;
212  xa[0] = (1 - t) * (1 - t) * (1 + 2*t);
213  xb[0] = h * t * (1 - t) * (1 - t);
214  dxa[0] = -6 * t * (1 - t) * h_1;
215  dxb[0] = (1 - t) * (1 - 3*t);
216  t--;
217  xa[1] = (1 + t) * (1 + t) * (1 - 2*t);
218  xb[1] = h * t * (1 + t) * (1 + t);
219  dxa[1] = -6 * t * (1 + t) * h_1;
220  dxb[1] = (1 + t) * (1 + 3*t);
221 
222  /* find t for y-dimension and compute ya, yb, dya, dyb */
223  t = dy_h - (double) jlo;
224  ya[0] = (1 - t) * (1 - t) * (1 + 2*t);
225  yb[0] = h * t * (1 - t) * (1 - t);
226  dya[0] = -6 * t * (1 - t) * h_1;
227  dyb[0] = (1 - t) * (1 - 3*t);
228  t--;
229  ya[1] = (1 + t) * (1 + t) * (1 - 2*t);
230  yb[1] = h * t * (1 + t) * (1 + t);
231  dya[1] = -6 * t * (1 + t) * h_1;
232  dyb[1] = (1 + t) * (1 + 3*t);
233 
234 #else
235 
236  /* find t for x-dimension and compute xa, xb, dxa, dxb */
237  t = dx_h - (double) ilo;
238  s1 = 1-t;
239  s2 = 2*t;
240  s3 = 3*t;
241  s4 = t*s1;
242  s5 = h*s4;
243  xa[0] = s1*s1*(1+s2);
244  xa[1] = t*t*(3-s2);
245  xb[0] = s5*s1;
246  xb[1] = -s5*t;
247  dxa[0] = -six_h*s4;
248  dxa[1] = -dxa[0];
249  dxb[0] = s1*(1-s3);
250  dxb[1] = t*(-2+s3);
251 
252  /* find t for y-dimension and compute ya, yb, dya, dyb */
253  t = dy_h - (double) jlo;
254  s1 = 1-t;
255  s2 = 2*t;
256  s3 = 3*t;
257  s4 = t*s1;
258  s5 = h*s4;
259  ya[0] = s1*s1*(1+s2);
260  ya[1] = t*t*(3-s2);
261  yb[0] = s5*s1;
262  yb[1] = -s5*t;
263  dya[0] = -six_h*s4;
264  dya[1] = -dya[0];
265  dyb[0] = s1*(1-s3);
266  dyb[1] = t*(-2+s3);
267 
268 #endif
269 
270  for (i = 0; i < 2; i++) {
271  for (j = 0; j < 2; j++) {
272  ij = INDEX(D,i+ilo,j+jlo);
273 
274 #if !defined(FASTER)
275 
276  f += xa[i] * ya[j] * table[ij].d00
277  + xb[i] * ya[j] * table[ij].d10
278  + xa[i] * yb[j] * table[ij].d01
279  + xb[i] * yb[j] * table[ij].d11;
280 
281  fx += dxa[i] * ya[j] * table[ij].d00
282  + dxb[i] * ya[j] * table[ij].d10
283  + dxa[i] * yb[j] * table[ij].d01
284  + dxb[i] * yb[j] * table[ij].d11;
285 
286  fy += xa[i] * dya[j] * table[ij].d00
287  + xb[i] * dya[j] * table[ij].d10
288  + xa[i] * dyb[j] * table[ij].d01
289  + xb[i] * dyb[j] * table[ij].d11;
290 
291 #else
292 
293  s1=ya[j]*table[ij].d00;
294  s2=yb[j]*table[ij].d01;
295  s3=ya[j]*table[ij].d10;
296  s4=yb[j]*table[ij].d11;
297 
298  f+=xa[i]*(s1+s2)+xb[i]*(s3+s4);
299  fx+=dxa[i]*(s1+s2)+dxb[i]*(s3+s4);
300  fy+=xa[i]*(dya[j]*table[ij].d00+dyb[j]*table[ij].d01)
301  +xb[i]*(dya[j]*table[ij].d10+dyb[j]*table[ij].d11);
302 
303 #endif
304  }
305  }
306 
307  /* return accumulated values */
308  U = f * scale;
309  U_phi = fx * scale;
310  U_psi = fy * scale;
311 
312 /*
313 CkPrintf("crossterm %d-%d-%d-%d %d-%d-%d-%d %lf %lf %d %d %lf %lf %lf\n",
314  atomID[0], atomID[1], atomID[2], atomID[3],
315  atomID[4], atomID[5], atomID[6], atomID[7],
316  phi, psi, ilo, jlo, U, U_phi, U_psi);
317 CkPrintf("%d %d-%d-%d-%d %d-%d-%d-%d\n", CkMyPe(),
318  p[0]->patchID, p[1]->patchID, p[2]->patchID, p[3]->patchID,
319  p[4]->patchID, p[5]->patchID, p[6]->patchID, p[7]->patchID);
320 */
321 
322 }
323 
324  // Add the energy from this crossterm to the total energy
325  energy += U;
326 
327  // Next, we want to calculate the forces. In order
328  // to do that, we first need to figure out whether the
329  // sin or cos form will be more stable. For this,
330  // just look at the value of phi
331  if (fabs(sin_phi) > 0.1)
332  {
333  // use the sin version to avoid 1/cos terms
334  U_phi = U_phi/sin_phi;
335 
336  f1.x += U_phi*(r23.y*dcosdA.z - r23.z*dcosdA.y);
337  f1.y += U_phi*(r23.z*dcosdA.x - r23.x*dcosdA.z);
338  f1.z += U_phi*(r23.x*dcosdA.y - r23.y*dcosdA.x);
339 
340  f3.x += U_phi*(r23.z*dcosdB.y - r23.y*dcosdB.z);
341  f3.y += U_phi*(r23.x*dcosdB.z - r23.z*dcosdB.x);
342  f3.z += U_phi*(r23.y*dcosdB.x - r23.x*dcosdB.y);
343 
344  f2.x += U_phi*(r12.z*dcosdA.y - r12.y*dcosdA.z
345  + r34.y*dcosdB.z - r34.z*dcosdB.y);
346  f2.y += U_phi*(r12.x*dcosdA.z - r12.z*dcosdA.x
347  + r34.z*dcosdB.x - r34.x*dcosdB.z);
348  f2.z += U_phi*(r12.y*dcosdA.x - r12.x*dcosdA.y
349  + r34.x*dcosdB.y - r34.y*dcosdB.x);
350  }
351  else
352  {
353  // This angle is closer to 0 or 180 than it is to
354  // 90, so use the cos version to avoid 1/sin terms
355  U_phi = -U_phi/cos_phi;
356 
357  f1.x += U_phi*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
358  - r23.x*r23.y*dsindC.y
359  - r23.x*r23.z*dsindC.z);
360  f1.y += U_phi*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
361  - r23.y*r23.z*dsindC.z
362  - r23.y*r23.x*dsindC.x);
363  f1.z += U_phi*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
364  - r23.z*r23.x*dsindC.x
365  - r23.z*r23.y*dsindC.y);
366 
367  f3 += cross(U_phi,dsindB,r23);
368 
369  f2.x += U_phi*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
370  +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
371  +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
372  +dsindB.z*r34.y - dsindB.y*r34.z);
373  f2.y += U_phi*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
374  +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
375  +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
376  +dsindB.x*r34.z - dsindB.z*r34.x);
377  f2.z += U_phi*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
378  +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
379  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
380  +dsindB.y*r34.x - dsindB.x*r34.y);
381  }
382 
383  if (fabs(sin_psi) > 0.1)
384  {
385  // use the sin version to avoid 1/cos terms
386  U_psi = U_psi/sin_psi;
387 
388  f4.x += U_psi*(r67.y*dcosdD.z - r67.z*dcosdD.y);
389  f4.y += U_psi*(r67.z*dcosdD.x - r67.x*dcosdD.z);
390  f4.z += U_psi*(r67.x*dcosdD.y - r67.y*dcosdD.x);
391 
392  f6.x += U_psi*(r67.z*dcosdE.y - r67.y*dcosdE.z);
393  f6.y += U_psi*(r67.x*dcosdE.z - r67.z*dcosdE.x);
394  f6.z += U_psi*(r67.y*dcosdE.x - r67.x*dcosdE.y);
395 
396  f5.x += U_psi*(r56.z*dcosdD.y - r56.y*dcosdD.z
397  + r78.y*dcosdE.z - r78.z*dcosdE.y);
398  f5.y += U_psi*(r56.x*dcosdD.z - r56.z*dcosdD.x
399  + r78.z*dcosdE.x - r78.x*dcosdE.z);
400  f5.z += U_psi*(r56.y*dcosdD.x - r56.x*dcosdD.y
401  + r78.x*dcosdE.y - r78.y*dcosdE.x);
402  }
403  else
404  {
405  // This angle is closer to 0 or 180 than it is to
406  // 90, so use the cos version to avoid 1/sin terms
407  U_psi = -U_psi/cos_psi;
408 
409  f4.x += U_psi*((r67.y*r67.y + r67.z*r67.z)*dsindF.x
410  - r67.x*r67.y*dsindF.y
411  - r67.x*r67.z*dsindF.z);
412  f4.y += U_psi*((r67.z*r67.z + r67.x*r67.x)*dsindF.y
413  - r67.y*r67.z*dsindF.z
414  - r67.y*r67.x*dsindF.x);
415  f4.z += U_psi*((r67.x*r67.x + r67.y*r67.y)*dsindF.z
416  - r67.z*r67.x*dsindF.x
417  - r67.z*r67.y*dsindF.y);
418 
419  f6 += cross(U_psi,dsindE,r67);
420 
421  f5.x += U_psi*(-(r67.y*r56.y + r67.z*r56.z)*dsindF.x
422  +(2.0*r67.x*r56.y - r56.x*r67.y)*dsindF.y
423  +(2.0*r67.x*r56.z - r56.x*r67.z)*dsindF.z
424  +dsindE.z*r78.y - dsindE.y*r78.z);
425  f5.y += U_psi*(-(r67.z*r56.z + r67.x*r56.x)*dsindF.y
426  +(2.0*r67.y*r56.z - r56.y*r67.z)*dsindF.z
427  +(2.0*r67.y*r56.x - r56.y*r67.x)*dsindF.x
428  +dsindE.x*r78.z - dsindE.z*r78.x);
429  f5.z += U_psi*(-(r67.x*r56.x + r67.y*r56.y)*dsindF.z
430  +(2.0*r67.z*r56.x - r56.z*r67.x)*dsindF.x
431  +(2.0*r67.z*r56.y - r56.z*r67.y)*dsindF.y
432  +dsindE.y*r78.x - dsindE.x*r78.y);
433  }
434 
435  //fepb - BKR scaling of alchemical bonded terms
436  // NB: TI derivative is the _unscaled_ energy.
437  if ( simParams->alchOn ) {
438  int typeSum1, typeSum2;
439  typeSum1 = typeSum2 = 0;
440  for (int i=0; i < 4; ++i ) {
441  typeSum1 += (mol->get_fep_type(atomID[i]) == 2 ? -1 :\
442  mol->get_fep_type(atomID[i]));
443  typeSum2 += (mol->get_fep_type(atomID[i+4]) == 2 ? -1 :\
444  mol->get_fep_type(atomID[i+4]));
445  }
446  int order = (simParams->alchBondDecouple ? 5 : 4);
447  if ( (0 < typeSum1 && typeSum1 < order) ||
448  (0 < typeSum2 && typeSum2 < order) ) {
449  reduction[crosstermEnergyIndex_ti_1] += energy;
450  reduction[crosstermEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy;
451  energy *= bond_lambda_1;
452  f1 *= bond_lambda_1;
453  f2 *= bond_lambda_1;
454  f3 *= bond_lambda_1;
455  f4 *= bond_lambda_1;
456  f5 *= bond_lambda_1;
457  f6 *= bond_lambda_1;
458  } else if ( (0 > typeSum1 && typeSum1 > -order) ||
459  (0 > typeSum2 && typeSum2 > -order) ) {
460  reduction[crosstermEnergyIndex_ti_2] += energy;
461  reduction[crosstermEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
462  energy *= bond_lambda_2;
463  f1 *= bond_lambda_2;
464  f2 *= bond_lambda_2;
465  f3 *= bond_lambda_2;
466  f4 *= bond_lambda_2;
467  f5 *= bond_lambda_2;
468  f6 *= bond_lambda_2;
469  }
470  }
471  //fepe
472 
473  /* store the forces */
474  p[0]->f[localIndex[0]] += f1;
475  p[1]->f[localIndex[1]] += f2 - f1;
476  p[2]->f[localIndex[2]] += f3 - f2;
477  p[3]->f[localIndex[3]] += -f3;
478  p[4]->f[localIndex[4]] += f4;
479  p[5]->f[localIndex[5]] += f5 - f4;
480  p[6]->f[localIndex[6]] += f6 - f5;
481  p[7]->f[localIndex[7]] += -f6;
482 
483  if ( p[0]->af ) {
484  p[0]->af[localIndex[0]] += f1;
485  p[1]->af[localIndex[1]] += f2 - f1;
486  p[2]->af[localIndex[2]] += f3 - f2;
487  p[3]->af[localIndex[3]] += -f3;
488  p[4]->af[localIndex[4]] += f4;
489  p[5]->af[localIndex[5]] += f5 - f4;
490  p[6]->af[localIndex[6]] += f6 - f5;
491  p[7]->af[localIndex[7]] += -f6;
492  }
493 
494  DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
495  reduction[crosstermEnergyIndex] += energy;
496  reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
497  reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
498  reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
499  reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
500  reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
501  reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
502  reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
503  reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
504  reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
505 
506  reduction[virialIndex_XX] += ( f4.x * r56.x + f5.x * r67.x + f6.x * r78.x );
507  reduction[virialIndex_XY] += ( f4.x * r56.y + f5.x * r67.y + f6.x * r78.y );
508  reduction[virialIndex_XZ] += ( f4.x * r56.z + f5.x * r67.z + f6.x * r78.z );
509  reduction[virialIndex_YX] += ( f4.y * r56.x + f5.y * r67.x + f6.y * r78.x );
510  reduction[virialIndex_YY] += ( f4.y * r56.y + f5.y * r67.y + f6.y * r78.y );
511  reduction[virialIndex_YZ] += ( f4.y * r56.z + f5.y * r67.z + f6.y * r78.z );
512  reduction[virialIndex_ZX] += ( f4.z * r56.x + f5.z * r67.x + f6.z * r78.x );
513  reduction[virialIndex_ZY] += ( f4.z * r56.y + f5.z * r67.y + f6.z * r78.y );
514  reduction[virialIndex_ZZ] += ( f4.z * r56.z + f5.z * r67.z + f6.z * r78.z );
515 
516  if (pressureProfileData) {
517  BigReal z1 = p[0]->x[localIndex[0]].position.z;
518  BigReal z2 = p[1]->x[localIndex[1]].position.z;
519  BigReal z3 = p[2]->x[localIndex[2]].position.z;
520  BigReal z4 = p[3]->x[localIndex[3]].position.z;
521  BigReal z5 = p[4]->x[localIndex[4]].position.z;
522  BigReal z6 = p[5]->x[localIndex[5]].position.z;
523  BigReal z7 = p[6]->x[localIndex[6]].position.z;
524  BigReal z8 = p[7]->x[localIndex[7]].position.z;
525  int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
526  int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
527  int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
528  int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
529  int n5 = (int)floor((z5-pressureProfileMin)/pressureProfileThickness);
530  int n6 = (int)floor((z6-pressureProfileMin)/pressureProfileThickness);
531  int n7 = (int)floor((z7-pressureProfileMin)/pressureProfileThickness);
532  int n8 = (int)floor((z8-pressureProfileMin)/pressureProfileThickness);
533  pp_clamp(n1, pressureProfileSlabs);
534  pp_clamp(n2, pressureProfileSlabs);
535  pp_clamp(n3, pressureProfileSlabs);
536  pp_clamp(n4, pressureProfileSlabs);
537  pp_clamp(n5, pressureProfileSlabs);
538  pp_clamp(n6, pressureProfileSlabs);
539  pp_clamp(n7, pressureProfileSlabs);
540  pp_clamp(n8, pressureProfileSlabs);
541  int p1 = p[0]->x[localIndex[0]].partition;
542  int p2 = p[1]->x[localIndex[1]].partition;
543  int p3 = p[2]->x[localIndex[2]].partition;
544  int p4 = p[3]->x[localIndex[3]].partition;
545  int p5 = p[4]->x[localIndex[4]].partition;
546  int p6 = p[5]->x[localIndex[5]].partition;
547  int p7 = p[6]->x[localIndex[6]].partition;
548  int p8 = p[7]->x[localIndex[7]].partition;
549  int pn = pressureProfileAtomTypes;
550  pp_reduction(pressureProfileSlabs, n1, n2, p1, p2, pn,
551  f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
552  pressureProfileData);
553  pp_reduction(pressureProfileSlabs, n2, n3, p2, p3, pn,
554  f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
555  pressureProfileData);
556  pp_reduction(pressureProfileSlabs, n3, n4, p3, p4, pn,
557  f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
558  pressureProfileData);
559  pp_reduction(pressureProfileSlabs, n5, n6, p5, p6, pn,
560  f4.x * r56.x, f4.y * r56.y, f4.z * r56.z,
561  pressureProfileData);
562  pp_reduction(pressureProfileSlabs, n6, n7, p6, p7, pn,
563  f5.x * r67.x, f5.y * r67.y, f5.z * r67.z,
564  pressureProfileData);
565  pp_reduction(pressureProfileSlabs, n7, n8, p7, p8, pn,
566  f6.x * r78.x, f6.y * r78.y, f6.z * r78.z,
567  pressureProfileData);
568  }
569 
570  }
571 }
Node::Object
static Node * Object()
Definition: Node.h:86
CrosstermElem::value
const CrosstermValue * value
Definition: ComputeCrossterms.h:41
CrosstermValue::c
CrosstermData c[dim][dim]
Definition: Parameters.h:145
Patch::lattice
Lattice & lattice
Definition: Patch.h:127
pp_reduction
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
Definition: PressureProfile.h:22
Vector
Definition: Vector.h:72
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21
INDEX
#define INDEX(ncols, i, j)
Definition: ComputeCrossterms.C:27
Real
float Real
Definition: common.h:118
DebugM
#define DebugM(x, y)
Definition: Debug.h:75
Vector::z
BigReal z
Definition: Vector.h:74
CompAtom::position
Position position
Definition: NamdTypes.h:77
CrosstermElem::pressureProfileSlabs
static int pressureProfileSlabs
Definition: ComputeCrossterms.h:35
Molecule
Molecule stores the structural information for the system.
Definition: Molecule.h:175
CrosstermElem::pressureProfileAtomTypes
static int pressureProfileAtomTypes
Definition: ComputeCrossterms.h:36
Patch::flags
Flags flags
Definition: Patch.h:128
CrosstermData::d11
BigReal d11
Definition: Parameters.h:137
PI
#define PI
Definition: common.h:92
CrosstermElem::p
TuplePatchElem * p[size]
Definition: ComputeCrossterms.h:23
pp_clamp
void pp_clamp(int &n, int nslabs)
Definition: PressureProfile.h:15
order
#define order
Definition: PmeRealSpace.C:235
Vector::length
NAMD_HOST_DEVICE BigReal length(void) const
Definition: Vector.h:202
CrosstermElem::pressureProfileMin
static BigReal pressureProfileMin
Definition: ComputeCrossterms.h:38
CrosstermElem::localIndex
int localIndex[size]
Definition: ComputeCrossterms.h:22
CrosstermData::d00
BigReal d00
Definition: Parameters.h:137
CompAtom::partition
uint8 partition
Definition: NamdTypes.h:80
Vector::x
BigReal x
Definition: Vector.h:74
Molecule::get_fep_type
unsigned char get_fep_type(int anum) const
Definition: Molecule.h:1432
simParams
#define simParams
Definition: Output.C:129
CrosstermData::d01
BigReal d01
Definition: Parameters.h:137
AtomID
int32 AtomID
Definition: NamdTypes.h:35
Vector::y
BigReal y
Definition: Vector.h:74
Node::molecule
Molecule * molecule
Definition: Node.h:179
CrosstermData::d10
BigReal d10
Definition: Parameters.h:137
BigReal
double BigReal
Definition: common.h:123
CrosstermElem::pressureProfileThickness
static BigReal pressureProfileThickness
Definition: ComputeCrossterms.h:37
Lattice::delta
NAMD_HOST_DEVICE Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:149
Flags::step
int step
Definition: PatchTypes.h:16

◆ getMoleculePointers()

void CrosstermElem::getMoleculePointers ( Molecule *  mol,
int *  count,
int32 ***  byatom,
Crossterm **  structarray 
)
static

Definition at line 38 of file ComputeCrossterms.C.

References NAMD_die(), and Molecule::numCrossterms.

39 {
40 #ifdef MEM_OPT_VERSION
41  NAMD_die("Should not be called in CrosstermElem::getMoleculePointers in memory optimized version!");
42 #else
43  *count = mol->numCrossterms;
44  *byatom = mol->crosstermsByAtom;
45  *structarray = mol->crossterms;
46 #endif
47 }
Molecule::numCrossterms
int numCrossterms
Definition: Molecule.h:596
NAMD_die
void NAMD_die(const char *err_msg)
Definition: common.C:147

◆ getParameterPointers()

void CrosstermElem::getParameterPointers ( Parameters *  p,
const CrosstermValue **  v 
)
static

Definition at line 49 of file ComputeCrossterms.C.

References p.

49  {
50  *v = p->crossterm_array;
51 }
CrosstermElem::p
TuplePatchElem * p[size]
Definition: ComputeCrossterms.h:23

◆ getTupleInfo()

static void CrosstermElem::getTupleInfo ( AtomSignature *  sig,
int *  count,
TupleSignature **  t 
)
inlinestatic

Definition at line 29 of file ComputeCrossterms.h.

References AtomSignature::crosstermCnt, and AtomSignature::crosstermSigs.

29  {
30  *count = sig->crosstermCnt;
31  *t = sig->crosstermSigs;
32  }
AtomSignature::crosstermSigs
TupleSignature * crosstermSigs
Definition: structures.h:339

◆ hash()

int CrosstermElem::hash ( void  ) const
inline

Definition at line 43 of file ComputeCrossterms.h.

References atomID.

43  {
44  return 0x7FFFFFFF &((atomID[0]<<24) + (atomID[1]<<16) + (atomID[2]<<8) + atomID[3]);
45  }
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21

◆ operator<()

int CrosstermElem::operator< ( const CrosstermElem &  a) const
inline

Definition at line 99 of file ComputeCrossterms.h.

References atomID.

99  {
100  return (atomID[0] < a.atomID[0] ||
101  (atomID[0] == a.atomID[0] &&
102  (atomID[1] < a.atomID[1] ||
103  (atomID[1] == a.atomID[1] &&
104  (atomID[2] < a.atomID[2] ||
105  (atomID[2] == a.atomID[2] &&
106  (atomID[3] < a.atomID[3] ||
107  (atomID[3] == a.atomID[3] &&
108  (atomID[4] < a.atomID[4] ||
109  (atomID[4] == a.atomID[4] &&
110  (atomID[5] < a.atomID[5] ||
111  (atomID[5] == a.atomID[5] &&
112  (atomID[6] < a.atomID[6] ||
113  (atomID[6] == a.atomID[6] &&
114  atomID[7] < a.atomID[7]
115  ))))))))))))));
116  }
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21

◆ operator==()

int CrosstermElem::operator== ( const CrosstermElem &  a) const
inline

Definition at line 92 of file ComputeCrossterms.h.

References atomID.

92  {
93  return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
94  a.atomID[2] == atomID[2] && a.atomID[3] == atomID[3] &&
95  a.atomID[4] == atomID[4] && a.atomID[5] == atomID[5] &&
96  a.atomID[6] == atomID[6] && a.atomID[7] == atomID[7]);
97  }
CrosstermElem::atomID
AtomID atomID[size]
Definition: ComputeCrossterms.h:21

◆ submitReductionData()

void CrosstermElem::submitReductionData ( BigReal *  data,
SubmitReduction *  reduction 
)
static

Definition at line 573 of file ComputeCrossterms.C.

References ADD_TENSOR, crosstermEnergyIndex, crosstermEnergyIndex_f, crosstermEnergyIndex_ti_1, crosstermEnergyIndex_ti_2, SubmitReduction::item(), REDUCTION_BONDED_ENERGY_F, REDUCTION_BONDED_ENERGY_TI_1, REDUCTION_BONDED_ENERGY_TI_2, and REDUCTION_CROSSTERM_ENERGY.

574 {
575  reduction->item(REDUCTION_CROSSTERM_ENERGY) += data[crosstermEnergyIndex];
576  reduction->item(REDUCTION_BONDED_ENERGY_F) += data[crosstermEnergyIndex_f];
577  reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[crosstermEnergyIndex_ti_1];
578  reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[crosstermEnergyIndex_ti_2];
579  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
580  ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
581 }
ADD_TENSOR
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:33
SubmitReduction::item
BigReal & item(int i)
Definition: ReductionMgr.h:313

Member Data Documentation

◆ atomID

AtomID CrosstermElem::atomID[size]

Definition at line 21 of file ComputeCrossterms.h.

Referenced by computeForce(), CrosstermElem(), hash(), operator<(), and operator==().

◆ localIndex

int CrosstermElem::localIndex[size]

Definition at line 22 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ p

TuplePatchElem* CrosstermElem::p[size]

Definition at line 23 of file ComputeCrossterms.h.

Referenced by computeForce(), and getParameterPointers().

◆ pressureProfileAtomTypes

int CrosstermElem::pressureProfileAtomTypes = 1
static

Definition at line 36 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ pressureProfileMin

BigReal CrosstermElem::pressureProfileMin = 0
static

Definition at line 38 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ pressureProfileSlabs

int CrosstermElem::pressureProfileSlabs = 0
static

Definition at line 35 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ pressureProfileThickness

BigReal CrosstermElem::pressureProfileThickness = 0
static

Definition at line 37 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ scale

Real CrosstermElem::scale

Definition at line 24 of file ComputeCrossterms.h.

Referenced by computeForce().

◆ value

const CrosstermValue* CrosstermElem::value

Definition at line 41 of file ComputeCrossterms.h.

Referenced by computeForce(), and CrosstermElem().

The documentation for this class was generated from the following files:
Generated by   doxygen 1.8.14