NAMD
ComputeRestraints.C
Go to the documentation of this file.
1 
7 #include "ComputeRestraints.h"
8 #include "Node.h"
9 #include "Molecule.h"
10 #include "SimParameters.h"
11 #include "Patch.h"
12 #include "NamdOneTools.h"
13 
14 /************************************************************************/
15 /* */
16 /* FUNCTION ComputeRestraints */
17 /* */
18 /************************************************************************/
19 
21  : ComputeHomePatch(c,pid)
22 {
24 
26 
27  // Get parameters from the SimParameters object
28  consExp = simParams->constraintExp;
29 
30  //****** BEGIN selective restraints (X,Y,Z) changes
31  consSelectOn = simParams->selectConstraintsOn;
32  if (consSelectOn) {
33  consSelectX = simParams->constrXOn;
34  consSelectY = simParams->constrYOn;
35  consSelectZ = simParams->constrZOn;
36  }
37  //****** END selective restraints (X,Y,Z) changes
38 
39  //****** BEGIN moving constraints changes
40  consMoveOn = simParams->movingConstraintsOn;
41  if (consMoveOn) {
42  moveVel = simParams->movingConsVel;
43  }
44  //****** END moving constraints changes
45  //****** BEGIN rotating constraints changes
46  consRotOn = simParams->rotConstraintsOn;
47  if (consRotOn) {
48  rotVel = simParams->rotConsVel;
49  rotAxis = simParams->rotConsAxis;
50  rotPivot = simParams->rotConsPivot;
51  }
52  //****** END rotating constraints changes
53 
54 }
55 /* END OF FUNCTION ComputeRestraints */
56 
57 /************************************************************************/
58 /* */
59 /* FUNCTION ~ComputeRestraints */
60 /* */
61 /* This is the destructor for the ComputeRestraints force object. */
62 /* */
63 /************************************************************************/
64 
66 
67 {
68  delete reduction;
69 }
70 /* END OF FUNCTION ~ComputeRestraints */
71 
72 /************************************************************************/
73 /* */
74 /* FUNCTION force */
75 /* */
76 /************************************************************************/
78 {
79  Molecule *molecule = Node::Object()->molecule;
80  Real k; // Force constant
82  BigReal scaling = simParams->constraintScaling;
83  Vector refPos; // Reference position
84  BigReal r, r2; // r=distance between atom position and the
85  // reference position, r2 = r^2
86  Vector Rij; // vector between current position and reference pos
87  BigReal value; // Current calculated value
88 
89  // aliases to work with old code
90  Force *f = res->f[Results::normal];
91  BigReal energy = 0;
92  BigReal m[9];
93  Tensor virial;
94  Vector netForce = 0;
95 
96  // BEGIN moving and rotating constraint changes ******
97 
98  // This version only allows one atom to be moved
99  // and only ALL ref positions to be rotated
100 
101  int currentTime = patch->flags.step;
102  if (consRotOn) {
103  vec_rotation_matrix(rotVel * currentTime, rotAxis, m);
104  }
105 
106  // END moving and rotating constraint changes ******
107 
108 
109  if (scaling != 0.) for (int localID=0; localID<numAtoms; ++localID)
110  {
111  if (molecule->is_atom_constrained(p[localID].id))
112  {
113 #ifndef MEM_OPT_VERSION
114  molecule->get_cons_params(k, refPos, p[localID].id);
115 #else
116  k = 1.0;
117  refPos = p[localID].fixedPosition;
118 #endif
119 
120  k *= scaling;
121 
122  // BEGIN moving and rotating constraint changes ******
123 
124  if (consMoveOn) {
125  refPos += currentTime * moveVel;
126  }
127  else if(consRotOn) {
128  refPos = mat_multiply_vec(refPos - rotPivot, m) + rotPivot;
129  }
130 
131  // END moving and rotating constraint changes *******
132 
133  if (simParams->sphericalConstraintsOn) {
134  BigReal refRad = (refPos - simParams->sphericalConstrCenter).length();
135  Vector relPos = patch->lattice.delta(p[localID].position, simParams->sphericalConstrCenter);
136  refPos = simParams->sphericalConstrCenter + relPos * (refRad/relPos.length());
137  }
138 
139  Rij = patch->lattice.delta(refPos,p[localID].position);
140  Vector vpos = refPos - Rij;
141 
142  //****** BEGIN selective restraints (X,Y,Z) changes
143  if (consSelectOn) { // check which components we want to restrain:
144  if (!consSelectX) {Rij.x=0.0;} // by setting the appropriate
145  if (!consSelectY) {Rij.y=0.0;} // Cartesian component to zero
146  if (!consSelectZ) {Rij.z=0.0;} // we will avoid a restoring force
147  } // along that component, and the
148  // energy will also be correct
149  //****** END selective restraints (X,Y,Z) changes
150 
151 
152  // Calculate the distance and the distance squared
153  r2 = Rij.length2();
154  r = sqrt(r2);
155 
156  // Only calculate the energy and the force if the distance is
157  // non-zero. Otherwise, you could end up dividing by 0, which
158  // is bad
159  if (r>0.0)
160  {
161  value=k;
162 
163  // Loop through and multiple k by r consExp times.
164  // i.e. calculate kr^e
165  // I know, I could use pow(), but I don't trust it.
166  for (int k=0; k<consExp; ++k)
167  {
168  value *= r;
169  }
170 
171  // Add to the energy total
172  energy += value;
173 
174  // Now calculate the force, which is ekr^(e-1). Also, divide
175  // by another factor of r to normalize the vector before we
176  // multiple it by the magnitude
177  value *= consExp;
178  value /= r2;
179 
180  Rij *= value;
181  //iout << iINFO << "restraining force" << Rij << "\n";
182 
183  f[localID] += Rij;
184  netForce += Rij;
185  virial += outer(Rij,vpos);
186  }
187  }
188  }
189 
190  reduction->item(REDUCTION_BC_ENERGY) += energy;
191  ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
192  ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,netForce);
193  reduction->submit();
194 
195 }
196 /* END OF FUNCTION force */
197 
static Node * Object()
Definition: Node.h:86
zVector rotConsAxis
int ComputeID
Definition: NamdTypes.h:183
Position fixedPosition
Definition: NamdTypes.h:102
Lattice & lattice
Definition: Patch.h:126
zVector sphericalConstrCenter
Definition: Vector.h:64
#define ADD_TENSOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:43
SimParameters * simParameters
Definition: Node.h:178
BigReal constraintScaling
Bool selectConstraintsOn
float Real
Definition: common.h:109
BigReal & item(int i)
Definition: ReductionMgr.h:312
Bool movingConstraintsOn
void vec_rotation_matrix(BigReal angle, Vector v, BigReal m[])
Definition: NamdOneTools.C:132
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
zVector rotConsPivot
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
BigReal length(void) const
Definition: Vector.h:169
Bool sphericalConstraintsOn
Flags flags
Definition: Patch.h:127
BigReal rotConsVel
SubmitReduction * reduction
virtual void doForce(FullAtom *p, Results *r)
Vector mat_multiply_vec(const Vector &v, BigReal m[])
Definition: NamdOneTools.C:232
ComputeRestraints(ComputeID c, PatchID pid)
Force * f[maxNumForces]
Definition: PatchTypes.h:67
virtual ~ComputeRestraints()
Vector delta(const Position &pos1, const Position &pos2) const
Definition: Lattice.h:144
BigReal x
Definition: Vector.h:66
int PatchID
Definition: NamdTypes.h:182
Tensor outer(const Vector &v1, const Vector &v2)
Definition: Tensor.h:241
Bool is_atom_constrained(int atomnum) const
Definition: Molecule.h:1201
void get_cons_params(Real &k, Vector &refPos, int atomnum) const
Definition: Molecule.h:1266
#define simParams
Definition: Output.C:127
Definition: Tensor.h:15
#define ADD_VECTOR_OBJECT(R, RL, D)
Definition: ReductionMgr.h:27
void submit(void)
Definition: ReductionMgr.h:323
Molecule * molecule
Definition: Node.h:176
double BigReal
Definition: common.h:114
zVector movingConsVel
int step
Definition: PatchTypes.h:16