NAMD
ComputeAniso.C
Go to the documentation of this file.
1 
7 #include "InfoStream.h"
8 #include "ComputeAniso.h"
9 #include "Molecule.h"
10 #include "Parameters.h"
11 #include "Node.h"
12 #include "ReductionMgr.h"
13 #include "Lattice.h"
14 #include "PressureProfile.h"
15 #include "Debug.h"
16 
17 #define CALCULATE_ANISO
18 
19 // static initialization
24 
26  (Molecule* mol, int* count, int32*** byatom, Aniso** structarray)
27 {
28 #ifdef MEM_OPT_VERSION
29  NAMD_die("Should not be called in AnisoElem::getMoleculePointers in memory optimized version!");
30 #else
31  *count = mol->numAnisos;
32  *byatom = mol->anisosByAtom;
33  *structarray = mol->anisos;
34 #endif
35 }
36 
38  *v = NULL; // parameters are stored in the structure
39 }
40 
41 void AnisoElem::computeForce(AnisoElem *tuples, int ntuple, BigReal *reduction,
42  BigReal *pressureProfileData)
43 {
44  const Lattice & lattice = tuples[0].p[0]->p->lattice;
45 
46  //fepb BKR
48  const int step = tuples[0].p[0]->p->flags.step;
49  const BigReal alchLambda = simParams->getCurrentLambda(step);
50  const BigReal alchLambda2 = simParams->getCurrentLambda2(step);
51  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
52  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
53  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
54  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
55  Molecule *const mol = Node::Object()->molecule;
56  //fepe
57 
58  for ( int ituple=0; ituple<ntuple; ++ituple ) {
59  const AnisoElem &tup = tuples[ituple];
60  enum { size = 4 };
61  const AtomID (&atomID)[size](tup.atomID);
62  const int (&localIndex)[size](tup.localIndex);
63  TuplePatchElem * const(&p)[size](tup.p);
64  const Real (&scale)(tup.scale);
65  const AnisoValue * const(&value)(tup.value);
66 
67  DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
68  << localIndex[1] << " " << localIndex[2] << " "
69  << localIndex[3] << std::endl);
70 
71 #ifdef CALCULATE_ANISO
72  // used some comments from Ed Harder's implementation in CHARMM
73 
74  const BigReal kpar0 = 2*value->k11; // force constants
75  const BigReal kperp0 = 2*value->k22;
76  const BigReal kiso0 = 2*value->k33;
77 
78  const Position & ri = p[0]->x[localIndex[0]].position; // atom I
79  const Position & rj = p[0]->x[localIndex[0]+1].position; // atom I's Drude
80  const Position & rl = p[1]->x[localIndex[1]].position; // atom L
81  const Position & rm = p[2]->x[localIndex[2]].position; // atom M
82  const Position & rn = p[3]->x[localIndex[3]].position; // atom N
83 
84  // calculate parallel and perpendicular displacement vectors
85  Vector r_il = lattice.delta(ri,rl); // shortest vector image: ri - rl
86  Vector r_mn = lattice.delta(rm,rn); // shortest vector image: rm - rn
87 
88  BigReal r_il_invlen = r_il.rlength(); // need recip lengths of r_il, r_mn
89  BigReal r_mn_invlen = r_mn.rlength();
90 
91  Vector u1 = r_il * r_il_invlen; // normalize r_il, r_mn
92  Vector u2 = r_mn * r_mn_invlen;
93 
94  Vector dr = rj - ri; // Drude displacement vector (ri, rj are in same patch)
95 
96  BigReal dpar = dr * u1; // parallel displacement
97  BigReal dperp = dr * u2; // perpendicular displacement
98 
99  // aniso spring energy
100  // kpar reduces response along carbonyl vector
101  // kperp reduces response perp to bond vector
102  // (reg in and out of plane response)
103  BigReal eaniso;
104  eaniso = 0.5*kpar0*dpar*dpar + 0.5*kperp0*dperp*dperp + 0.5*kiso0*(dr*dr);
105 
106  // calculate force vectors in one direction only:
107  // fi = -(fj + fl), fn = -fm
108 
109  // force on atom j
110  Vector fj = -kiso0 * dr;
111  fj -= kpar0 * dpar * u1;
112  fj -= kperp0 * dperp * u2;
113 
114  // force on atom l
115  Vector fl = kpar0 * dpar * r_il_invlen * dr;
116  fl -= kpar0 * dpar * dpar * r_il_invlen * u1;
117 
118  // force on atom m
119  Vector fm = kperp0 * dperp * dperp * r_mn_invlen * u2;
120  fm -= kperp0 * dperp * r_mn_invlen * dr;
121 
122  //fepb - BKR scaling of alchemical bonded terms
123  // NB: TI derivative is the _unscaled_ energy.
124  if ( simParams->alchOn && !simParams->singleTopology) {
125  switch ( mol->get_fep_bonded_type(atomID, 4) ) {
126  case 1:
127  reduction[anisoEnergyIndex_ti_1] += eaniso;
128  reduction[anisoEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*eaniso;
129  eaniso *= bond_lambda_1;
130  fj *= bond_lambda_1;
131  fl *= bond_lambda_1;
132  fm *= bond_lambda_1;
133  break;
134  case 2:
135  reduction[anisoEnergyIndex_ti_2] += eaniso;
136  reduction[anisoEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*eaniso;
137  eaniso *= bond_lambda_2;
138  fj *= bond_lambda_2;
139  fl *= bond_lambda_2;
140  fm *= bond_lambda_2;
141  break;
142  }
143  }
144  //fepe
145 
146  // accumulate forces
147  p[0]->f[localIndex[0]] -= (fj + fl);
148  p[0]->f[localIndex[0]+1] += fj;
149  p[1]->f[localIndex[1]] += fl;
150  p[2]->f[localIndex[2]] += fm;
151  p[3]->f[localIndex[3]] -= fm;
152 
153  // update potential
154  reduction[anisoEnergyIndex] += eaniso;
155 
156  // update virial
157  reduction[virialIndex_XX] += fj.x * dr.x - fl.x * r_il.x + fm.x * r_mn.x;
158  reduction[virialIndex_XY] += fj.x * dr.y - fl.x * r_il.y + fm.x * r_mn.y;
159  reduction[virialIndex_XZ] += fj.x * dr.z - fl.x * r_il.z + fm.x * r_mn.z;
160  reduction[virialIndex_YX] += fj.y * dr.x - fl.y * r_il.x + fm.y * r_mn.x;
161  reduction[virialIndex_YY] += fj.y * dr.y - fl.y * r_il.y + fm.y * r_mn.y;
162  reduction[virialIndex_YZ] += fj.y * dr.z - fl.y * r_il.z + fm.y * r_mn.z;
163  reduction[virialIndex_ZX] += fj.z * dr.x - fl.z * r_il.x + fm.z * r_mn.x;
164  reduction[virialIndex_ZY] += fj.z * dr.y - fl.z * r_il.y + fm.z * r_mn.y;
165  reduction[virialIndex_ZZ] += fj.z * dr.z - fl.z * r_il.z + fm.z * r_mn.z;
166 
167  // update pressure profile data
168  if (pressureProfileData) {
169  BigReal zi = p[0]->x[localIndex[0]].position.z;
170  BigReal zj = p[0]->x[localIndex[0]+1].position.z;
171  BigReal zl = p[1]->x[localIndex[1]].position.z;
172  BigReal zm = p[2]->x[localIndex[2]].position.z;
173  BigReal zn = p[3]->x[localIndex[3]].position.z;
174  int ni = (int)floor((zi-pressureProfileMin)/pressureProfileThickness);
175  int nj = (int)floor((zj-pressureProfileMin)/pressureProfileThickness);
176  int nl = (int)floor((zl-pressureProfileMin)/pressureProfileThickness);
177  int nm = (int)floor((zm-pressureProfileMin)/pressureProfileThickness);
178  int nn = (int)floor((zn-pressureProfileMin)/pressureProfileThickness);
184  int pi = p[0]->x[localIndex[0]].partition;
185  int pj = p[0]->x[localIndex[0]+1].partition;
186  int pl = p[1]->x[localIndex[1]].partition;
187  int pm = p[2]->x[localIndex[2]].partition;
188  int pn = p[3]->x[localIndex[3]].partition;
189  int pt = pressureProfileAtomTypes;
191  pj, pi, pt, fj.x * dr.x, fj.y * dr.y, fj.z * dr.z,
192  pressureProfileData);
194  pi, pl, pt, -fl.x * r_il.x, -fl.y * r_il.y, -fl.z * r_il.z,
195  pressureProfileData);
197  pm, pn, pt, fm.x * r_mn.x, fm.y * r_mn.y, fm.z * r_mn.z,
198  pressureProfileData);
199  }
200 #endif
201 
202  }
203 }
204 
205 
207 {
208  reduction->item(REDUCTION_BOND_ENERGY) += data[anisoEnergyIndex];
212  ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
213 }
214 
static Node * Object()
Definition: Node.h:86
unsigned char partition
Definition: NamdTypes.h:56
Real k11
Definition: structures.h:136
short int32
Definition: dumpdcd.c:24
int AtomID
Definition: NamdTypes.h:29
Lattice & lattice
Definition: Patch.h:126
#define ADD_TENSOR(R, RL, D, DL)
Definition: ReductionMgr.h:32
int numAnisos
Number of anisotropic terms.
Definition: Molecule.h:583
void pp_reduction(int nslabs, int n1, int n2, int atype1, int atype2, int numtypes, BigReal vxx, BigReal vyy, BigReal vzz, BigReal *reduction)
Definition: Vector.h:64
static BigReal pressureProfileMin
Definition: ComputeAniso.h:40
SimParameters * simParameters
Definition: Node.h:178
float Real
Definition: common.h:107
BigReal & item(int i)
Definition: ReductionMgr.h:312
#define DebugM(x, y)
Definition: Debug.h:59
BigReal z
Definition: Vector.h:66
Position position
Definition: NamdTypes.h:53
static int pressureProfileSlabs
Definition: ComputeAniso.h:37
int get_fep_bonded_type(const int *atomID, unsigned int order) const
Definition: Molecule.h:1385
Real scale
Definition: ComputeAniso.h:25
static void computeForce(AnisoElem *, int, BigReal *, BigReal *)
Definition: ComputeAniso.C:41
static BigReal pressureProfileThickness
Definition: ComputeAniso.h:39
Flags flags
Definition: Patch.h:127
void pp_clamp(int &n, int nslabs)
BigReal rlength(void)
Definition: Vector.h:177
BigReal getBondLambda(const BigReal)
BigReal getCurrentLambda2(const int)
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
void NAMD_die(const char *err_msg)
Definition: common.C:83
int localIndex[size]
Definition: ComputeAniso.h:23
Real k22
Definition: structures.h:137
Real k33
Definition: structures.h:138
#define simParams
Definition: Output.C:127
static int pressureProfileAtomTypes
Definition: ComputeAniso.h:38
static void submitReductionData(BigReal *, SubmitReduction *)
Definition: ComputeAniso.C:206
BigReal y
Definition: Vector.h:66
BigReal getCurrentLambda(const int)
AtomID atomID[size]
Definition: ComputeAniso.h:22
TuplePatchElem * p[size]
Definition: ComputeAniso.h:24
static void getMoleculePointers(Molecule *, int *, int32 ***, Aniso **)
Definition: ComputeAniso.C:26
const AnisoValue * value
Definition: ComputeAniso.h:43
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:112
int step
Definition: PatchTypes.h:16
static void getParameterPointers(Parameters *, const AnisoValue **)
Definition: ComputeAniso.C:37