NAMD
ComputeCylindricalBC.C
Go to the documentation of this file.
1 
7 #include "ComputeCylindricalBC.h"
8 #include "Node.h"
9 #include "SimParameters.h"
10 #include "Patch.h"
11 
12 /************************************************************************/
13 /* */
14 /* FUNCTION ComputeCylindricalBC */
15 /* */
16 /* This is the constructor for the ComputeCylindricalBC force object.*/
17 /* It is responsible for getting all the parameters from the */
18 /* SimParameters object and then determining if this object needs */
19 /* to perform any computation. It only needs to do so if there is */
20 /* some portion of the patch that lays outside of the cylindrical */
21 /* boundaries. */
22 /* */
23 /************************************************************************/
24 
25 ComputeCylindricalBC::ComputeCylindricalBC(ComputeID c, PatchID pid)
26  : ComputeHomePatch(c,pid)
27 {
28  reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
29 
30  SimParameters *simParams = Node::Object()->simParameters;
31 
32  // Get parameters from the SimParameters object
33  axis = simParams->cylindricalBCAxis;
34  r1 = simParams->cylindricalBCr1;
35  r2 = simParams->cylindricalBCr2;
36  r1_2 = r1*r1;
37  r2_2 = r2*r2;
38  k1 = simParams->cylindricalBCk1;
39  k2 = simParams->cylindricalBCk2;
40  exp1 = simParams->cylindricalBCexp1;
41  exp2 = simParams->cylindricalBCexp2;
42  //Additions for ends
43  l1 = simParams->cylindricalBCl1;
44  l2 = simParams->cylindricalBCl2;
45  l1_2 = l1*l1;
46  l2_2 = l2*l2;
47 
48  // Check to see if this is one set of parameters or two
49  if (r2 > -1.0)
50  {
51  twoForces = TRUE;
52  }
53  else
54  {
55  twoForces = FALSE;
56  }
57 
58  center = simParams->cylindricalCenter;
59 
60 }
61 /* END OF FUNCTION ComputeCylindricalBC */
62 
63 /************************************************************************/
64 /* */
65 /* FUNCTION ~ComputeCylindricalBC */
66 /* */
67 /* This is the destructor for the ComputeCylindricalBC force object. */
68 /* It currently does ABSOLUTELY NOTHING!! */
69 /* */
70 /************************************************************************/
71 
72 ComputeCylindricalBC::~ComputeCylindricalBC()
73 
74 {
75  delete reduction;
76 }
77 /* END OF FUNCTION ~ComputeCylindricalBC */
78 
79 /************************************************************************/
80 /* */
81 /* FUNCTION force */
82 /* */
83 /* INPUTS: */
84 /* numAtoms - Number of coordinates being passed */
85 /* x - Array of atom coordinates */
86 /* forces - Array of force vectors */
87 /* */
88 /* This function calculates the force and energy for the cylindri. */
89 /* boundary conditions for this patch. */
90 /* */
91 /************************************************************************/
92 
93 void ComputeCylindricalBC::doForce(FullAtom* p, Results* r)
94 {
95  Vector diff; // Distance from atom to center of cylinder
96  Vector f; // Calculated force vector
97  int i, j; // Loop counters
98  BigReal dist, dist_2; // Distance from atom to center, and this
99  // distance squared
100  BigReal rval; // Difference between distance from atom
101  // to center and radius of cylinder
102  BigReal eval; // Energy value for this atom
103  BigReal fval; // Force magnitude for this atom
104 
105  // aliases to work with old code
106  FullAtom *x = p;
107  Force *forces = r->f[Results::normal];
108  BigReal energy = 0;
109 
110  // Loop through and check each atom
111  for (i=0; i<numAtoms; i++)
112  {
113  // Calculate the vector from the atom to the center of the
114  // cylinder
115  diff.x = ( axis == 'x' ? 0.0 : x[i].position.x - center.x );
116  diff.y = ( axis == 'y' ? 0.0 : x[i].position.y - center.y );
117  diff.z = ( axis == 'z' ? 0.0 : x[i].position.z - center.z );
118 
119  // Calculate the distance squared
120  dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
121 
122  // Look to see if we are outside either radius
123  if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
124  {
125  // Calculate the distance to the center
126  dist = sqrt(dist_2);
127 
128  // Normalize the direction vector
129  diff /= dist;
130 
131  // Check to see if we are outside radius 1
132  if (dist > r1)
133  {
134  // Assign the force vector to the
135  // unit direction vector
136  f.x = diff.x;
137  f.y = diff.y;
138  f.z = diff.z;
139 
140  // Calculate the energy which is
141  // e = k1*(r_i-r_center)^exp1
142  eval = k1;
143  rval = fabs(dist - r1);
144 
145  for (j=0; j<exp1; j++)
146  {
147  eval *= rval;
148  }
149 
150  energy += eval;
151 
152  // Now calculate the force which is
153  // e = -k1*exp1*(r_i-r_center1)^(exp1-1)
154  fval = -exp1*k1;
155 
156  for (j=0; j<exp1-1; j++)
157  {
158  fval *= rval;
159  }
160 
161  // Multiply the force magnitude to the
162  // unit direction vector to get the
163  // resulting force
164  f *= fval;
165 
166  // Add the force to the force vectors
167  forces[i].x += f.x;
168  forces[i].y += f.y;
169  forces[i].z += f.z;
170  }
171 
172  // Check to see if there is a second radius
173  // and if we are outside of it
174  if (twoForces && (dist > r2) )
175  {
176  // Assign the force vector to the
177  // unit direction vector
178  f.x = diff.x;
179  f.y = diff.y;
180  f.z = diff.z;
181 
182  // Calculate the energy which is
183  // e = k2*(r_i-r_center2)^exp2
184  eval = k2;
185  rval = fabs(dist - r2);
186 
187  for (j=0; j<exp2; j++)
188  {
189  eval *= rval;
190  }
191 
192  energy += eval;
193 
194  // Now calculate the force which is
195  // e = -k2*exp2*(r_i-r_center2)^(exp2-1)
196  fval = -exp2*k2;
197 
198  for (j=0; j<exp2-1; j++)
199  {
200  fval *= rval;
201  }
202 
203  // Multiply the force magnitude to the
204  // unit direction vector to get the
205  // resulting force
206  f *= fval;
207 
208  // Add the force to the force vectors
209  forces[i].x += f.x;
210  forces[i].y += f.y;
211  forces[i].z += f.z;
212  }
213  }
214  }
215  // Loop through and check each atom
216  for (i=0; i<numAtoms; i++)
217  {
218  // Calculate the vector from the atom to the center of the
219  // cylinder
220  diff.x = ( axis != 'x' ? 0.0 : x[i].position.x - center.x );
221  diff.y = ( axis != 'y' ? 0.0 : x[i].position.y - center.y );
222  diff.z = ( axis != 'z' ? 0.0 : x[i].position.z - center.z );
223 
224  // Calculate the distance squared
225  dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
226 
227  // Look to see if we are outside either radius
228  if ( (dist_2 > l1_2) || (twoForces && (dist_2 > l2_2)) )
229  {
230  // Calculate the distance to the center
231  dist = sqrt(dist_2);
232 
233  // Normalize the direction vector
234  diff /= dist;
235 
236  // Check to see if we are outside radius 1
237  if (dist > l1)
238  {
239 //printf ("Do faces one force z=%f\n", dist);
240  // Assign the force vector to the
241  // unit direction vector
242  f.x = diff.x;
243  f.y = diff.y;
244  f.z = diff.z;
245 
246  // Calculate the energy which is
247  // e = k1*(r_i-r_center)^exp1
248  eval = k1;
249  rval = fabs(dist - l1);
250 
251  for (j=0; j<exp1; j++)
252  {
253  eval *= rval;
254  }
255 
256  energy += eval;
257 
258  // Now calculate the force which is
259  // e = -k1*exp1*(r_i-r_center1)^(exp1-1)
260  fval = -exp1*k1;
261 
262  for (j=0; j<exp1-1; j++)
263  {
264  fval *= rval;
265  }
266 
267  // Multiply the force magnitude to the
268  // unit direction vector to get the
269  // resulting force
270  f *= fval;
271 
272  // Add the force to the force vectors
273  forces[i].x += f.x;
274  forces[i].y += f.y;
275  forces[i].z += f.z;
276  }
277 
278  // Check to see if there is a second radius
279  // and if we are outside of it
280  if (twoForces && (dist > l2) )
281  {
282 //printf ("Do faces two force z=%f\n", dist);
283  // Assign the force vector to the
284  // unit direction vector
285  f.x = diff.x;
286  f.y = diff.y;
287  f.z = diff.z;
288 
289  // Calculate the energy which is
290  // e = k2*(r_i-r_center2)^exp2
291  eval = k2;
292  rval = fabs(dist - l2);
293 
294  for (j=0; j<exp2; j++)
295  {
296  eval *= rval;
297  }
298 
299  energy += eval;
300 
301  // Now calculate the force which is
302  // e = -k2*exp2*(r_i-r_center2)^(exp2-1)
303  fval = -exp2*k2;
304 
305  for (j=0; j<exp2-1; j++)
306  {
307  fval *= rval;
308  }
309 
310  // Multiply the force magnitude to the
311  // unit direction vector to get the
312  // resulting force
313  f *= fval;
314 
315  // Add the force to the force vectors
316  forces[i].x += f.x;
317  forces[i].y += f.y;
318  forces[i].z += f.z;
319  }
320  }
321  }
322 
323  reduction->item(REDUCTION_BC_ENERGY) += energy;
324  reduction->submit();
325 
326 }
327 /* END OF FUNCTION force */
328 
Node::Object
static Node * Object()
Definition: Node.h:86
ComputeCylindricalBC::reduction
SubmitReduction * reduction
Definition: ComputeCylindricalBC.h:37
ComputeID
int32 ComputeID
Definition: NamdTypes.h:278
Vector
Definition: Vector.h:72
Node::simParameters
SimParameters * simParameters
Definition: Node.h:181
SubmitReduction::item
BigReal & item(int i)
Definition: ReductionMgr.h:313
ComputeCylindricalBC::doForce
virtual void doForce(FullAtom *p, Results *r)
Definition: ComputeCylindricalBC.C:93
Vector::z
BigReal z
Definition: Vector.h:74
FALSE
#define FALSE
Definition: common.h:127
CompAtom::position
Position position
Definition: NamdTypes.h:77
ReductionMgr::willSubmit
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:366
ReductionMgr::Object
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:279
Results::f
Force * f[maxNumForces]
Definition: PatchTypes.h:146
Vector::x
BigReal x
Definition: Vector.h:74
ComputeCylindricalBC::ComputeCylindricalBC
ComputeCylindricalBC(ComputeID c, PatchID pid)
Definition: ComputeCylindricalBC.C:25
simParams
#define simParams
Definition: Output.C:129
Vector::y
BigReal y
Definition: Vector.h:74
SubmitReduction::submit
void submit(void)
Definition: ReductionMgr.h:324
PatchID
int32 PatchID
Definition: NamdTypes.h:277
TRUE
#define TRUE
Definition: common.h:128
BigReal
double BigReal
Definition: common.h:123
Generated by   doxygen 1.8.14