NAMD
ComputeSphericalBC.C
Go to the documentation of this file.
1 
7 #include "ComputeSphericalBC.h"
8 #include "Node.h"
9 #include "SimParameters.h"
10 #include "Patch.h"
11 
12 /************************************************************************/
13 /* */
14 /* FUNCTION ComputeSphericalBC */
15 /* */
16 /* This is the constructor for the ComputeSphericalBC 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 spherical */
21 /* boundaries. */
22 /* */
23 /************************************************************************/
24 
26  : ComputeHomePatch(c,pid)
27 {
29 
31 
32  // Get parameters from the SimParameters object
33  r1 = simParams->sphericalBCr1;
34  r2 = simParams->sphericalBCr2;
35  r1_2 = r1*r1;
36  r2_2 = r2*r2;
37  k1 = simParams->sphericalBCk1;
38  k2 = simParams->sphericalBCk2;
39  exp1 = simParams->sphericalBCexp1;
40  exp2 = simParams->sphericalBCexp2;
41 
42  // Check to see if this is one set of parameters or two
43  if (r2 > -1.0)
44  {
45  twoForces = TRUE;
46  }
47  else
48  {
49  twoForces = FALSE;
50  }
51 
52  center = simParams->sphericalCenter;
53 
54 }
55 /* END OF FUNCTION ComputeSphericalBC */
56 
57 /************************************************************************/
58 /* */
59 /* FUNCTION ~ComputeSphericalBC */
60 /* */
61 /* This is the destructor for the ComputeSphericalBC force object. */
62 /* It currently does ABSOLUTELY NOTHING!! */
63 /* */
64 /************************************************************************/
65 
67 
68 {
69  delete reduction;
70 }
71 /* END OF FUNCTION ~ComputeSphericalBC */
72 
73 /************************************************************************/
74 /* */
75 /* FUNCTION force */
76 /* */
77 /* INPUTS: */
78 /* numAtoms - Number of coordinates being passed */
79 /* x - Array of atom coordinates */
80 /* forces - Array of force vectors */
81 /* */
82 /* This function calculates the force and energy for the shereical */
83 /* boundary conditions for this patch. */
84 /* */
85 /************************************************************************/
87 {
88  Vector diff; // Distance from atom to center of sphere
89  Vector f; // Calculated force vector
90  int i, j; // Loop counters
91  BigReal dist, dist_2; // Distance from atom to center, and this
92  // distance squared
93  BigReal rval; // Difference between distance from atom
94  // to center and radius of sphere
95  BigReal eval; // Energy value for this atom
96  BigReal fval; // Force magnitude for this atom
97 
98  // aliases to work with old code
99  FullAtom *x = p;
100  Force *forces = r->f[Results::normal];
101  BigReal energy = 0;
102 
103  // Loop through and check each atom
104  for (i=0; i<numAtoms; i++)
105  {
106  // Calculate the vector from the atom to the center of the
107  // sphere
108  diff.x = x[i].position.x - center.x;
109  diff.y = x[i].position.y - center.y;
110  diff.z = x[i].position.z - center.z;
111 
112  // Calculate the distance squared
113  dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
114 
115  // Look to see if we are outside either radius
116  if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
117  {
118  // Calculate the distance to the center
119  dist = sqrt(dist_2);
120 
121  // Normalize the direction vector
122  diff /= dist;
123 
124  // Check to see if we are outside radius 1
125  if (dist > r1)
126  {
127  // Assign the force vector to the
128  // unit direction vector
129  f.x = diff.x;
130  f.y = diff.y;
131  f.z = diff.z;
132 
133  // Calculate the energy which is
134  // e = k1*(r_i-r_center)^exp1
135  eval = k1;
136  rval = fabs(dist - r1);
137 
138  for (j=0; j<exp1; j++)
139  {
140  eval *= rval;
141  }
142 
143  energy += eval;
144 
145  // Now calculate the force which is
146  // e = -k1*exp1*(r_i-r_center1)^(exp1-1)
147  fval = -exp1*k1;
148 
149  for (j=0; j<exp1-1; j++)
150  {
151  fval *= rval;
152  }
153 
154  // Multiply the force magnitude to the
155  // unit direction vector to get the
156  // resulting force
157  f *= fval;
158 
159  // Add the force to the force vectors
160  forces[i].x += f.x;
161  forces[i].y += f.y;
162  forces[i].z += f.z;
163  }
164 
165  // Check to see if there is a second radius
166  // and if we are outside of it
167  if (twoForces && (dist > r2) )
168  {
169  // Assign the force vector to the
170  // unit direction vector
171  f.x = diff.x;
172  f.y = diff.y;
173  f.z = diff.z;
174 
175  // Calculate the energy which is
176  // e = k2*(r_i-r_center2)^exp2
177  eval = k2;
178  rval = fabs(dist - r2);
179 
180  for (j=0; j<exp2; j++)
181  {
182  eval *= rval;
183  }
184 
185  energy += eval;
186 
187  // Now calculate the force which is
188  // e = -k2*exp2*(r_i-r_center2)^(exp2-1)
189  fval = -exp2*k2;
190 
191  for (j=0; j<exp2-1; j++)
192  {
193  fval *= rval;
194  }
195 
196  // Multiply the force magnitude to the
197  // unit direction vector to get the
198  // resulting force
199  f *= fval;
200 
201  // Add the force to the force vectors
202  forces[i].x += f.x;
203  forces[i].y += f.y;
204  forces[i].z += f.z;
205  }
206  }
207  }
208 
209  reduction->item(REDUCTION_BC_ENERGY) += energy;
210  reduction->submit();
211 
212 }
213 /* END OF FUNCTION force */
static Node * Object()
Definition: Node.h:86
zVector sphericalCenter
BigReal sphericalBCr2
int ComputeID
Definition: NamdTypes.h:183
virtual void doForce(FullAtom *p, Results *r)
Definition: Vector.h:64
SimParameters * simParameters
Definition: Node.h:178
static __thread float4 * forces
BigReal sphericalBCk1
BigReal & item(int i)
Definition: ReductionMgr.h:312
BigReal z
Definition: Vector.h:66
#define FALSE
Definition: common.h:116
Position position
Definition: NamdTypes.h:53
SubmitReduction * willSubmit(int setID, int size=-1)
Definition: ReductionMgr.C:365
static ReductionMgr * Object(void)
Definition: ReductionMgr.h:278
ComputeSphericalBC(ComputeID c, PatchID pid)
Force * f[maxNumForces]
Definition: PatchTypes.h:67
BigReal x
Definition: Vector.h:66
int PatchID
Definition: NamdTypes.h:182
BigReal sphericalBCr1
#define simParams
Definition: Output.C:127
BigReal y
Definition: Vector.h:66
void submit(void)
Definition: ReductionMgr.h:323
gridSize x
#define TRUE
Definition: common.h:117
BigReal sphericalBCk2
double BigReal
Definition: common.h:112
SubmitReduction * reduction