#include <ComputeCylindricalBC.h>
Inheritance diagram for ComputeCylindricalBC:

Public Member Functions | |
| ComputeCylindricalBC (ComputeID c, PatchID pid) | |
| virtual | ~ComputeCylindricalBC () |
| virtual void | doForce (FullAtom *p, Results *r) |
Public Attributes | |
| SubmitReduction * | reduction |
Definition at line 13 of file ComputeCylindricalBC.h.
|
||||||||||||
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 25 of file ComputeCylindricalBC.C. References SimParameters::cylindricalBCAxis, SimParameters::cylindricalBCexp1, SimParameters::cylindricalBCexp2, SimParameters::cylindricalBCk1, SimParameters::cylindricalBCk2, SimParameters::cylindricalBCl1, SimParameters::cylindricalBCl2, SimParameters::cylindricalBCr1, SimParameters::cylindricalBCr2, SimParameters::cylindricalCenter, Node::Object(), ReductionMgr::Object(), reduction, REDUCTIONS_BASIC, Node::simParameters, simParams, and ReductionMgr::willSubmit(). 00026 : ComputeHomePatch(c,pid) 00027 { 00028 reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC); 00029 00030 SimParameters *simParams = Node::Object()->simParameters; 00031 00032 // Get parameters from the SimParameters object 00033 axis = simParams->cylindricalBCAxis; 00034 r1 = simParams->cylindricalBCr1; 00035 r2 = simParams->cylindricalBCr2; 00036 r1_2 = r1*r1; 00037 r2_2 = r2*r2; 00038 k1 = simParams->cylindricalBCk1; 00039 k2 = simParams->cylindricalBCk2; 00040 exp1 = simParams->cylindricalBCexp1; 00041 exp2 = simParams->cylindricalBCexp2; 00042 //Additions for ends 00043 l1 = simParams->cylindricalBCl1; 00044 l2 = simParams->cylindricalBCl2; 00045 l1_2 = l1*l1; 00046 l2_2 = l2*l2; 00047 00048 // Check to see if this is one set of parameters or two 00049 if (r2 > -1.0) 00050 { 00051 twoForces = TRUE; 00052 } 00053 else 00054 { 00055 twoForces = FALSE; 00056 } 00057 00058 center = simParams->cylindricalCenter; 00059 00060 }
|
|
|
Definition at line 72 of file ComputeCylindricalBC.C. 00074 {
00075 delete reduction;
00076 }
|
|
||||||||||||
|
Implements ComputeHomePatch. Definition at line 93 of file ComputeCylindricalBC.C. References BigReal, Results::f, Force, SubmitReduction::item(), j, CompAtom::position, reduction, REDUCTION_BC_ENERGY, SubmitReduction::submit(), Vector::x, Vector::y, and Vector::z. 00094 {
00095 Vector diff; // Distance from atom to center of cylinder
00096 Vector f; // Calculated force vector
00097 int i, j; // Loop counters
00098 BigReal dist, dist_2; // Distance from atom to center, and this
00099 // distance squared
00100 BigReal rval; // Difference between distance from atom
00101 // to center and radius of cylinder
00102 BigReal eval; // Energy value for this atom
00103 BigReal fval; // Force magnitude for this atom
00104
00105 // aliases to work with old code
00106 FullAtom *x = p;
00107 Force *forces = r->f[Results::normal];
00108 BigReal energy = 0;
00109
00110 // Loop through and check each atom
00111 for (i=0; i<numAtoms; i++)
00112 {
00113 // Calculate the vector from the atom to the center of the
00114 // cylinder
00115 diff.x = ( axis == 'x' ? 0.0 : x[i].position.x - center.x );
00116 diff.y = ( axis == 'y' ? 0.0 : x[i].position.y - center.y );
00117 diff.z = ( axis == 'z' ? 0.0 : x[i].position.z - center.z );
00118
00119 // Calculate the distance squared
00120 dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
00121
00122 // Look to see if we are outside either radius
00123 if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
00124 {
00125 // Calculate the distance to the center
00126 dist = sqrt(dist_2);
00127
00128 // Normalize the direction vector
00129 diff /= dist;
00130
00131 // Check to see if we are outside radius 1
00132 if (dist > r1)
00133 {
00134 // Assign the force vector to the
00135 // unit direction vector
00136 f.x = diff.x;
00137 f.y = diff.y;
00138 f.z = diff.z;
00139
00140 // Calculate the energy which is
00141 // e = k1*(r_i-r_center)^exp1
00142 eval = k1;
00143 rval = fabs(dist - r1);
00144
00145 for (j=0; j<exp1; j++)
00146 {
00147 eval *= rval;
00148 }
00149
00150 energy += eval;
00151
00152 // Now calculate the force which is
00153 // e = -k1*exp1*(r_i-r_center1)^(exp1-1)
00154 fval = -exp1*k1;
00155
00156 for (j=0; j<exp1-1; j++)
00157 {
00158 fval *= rval;
00159 }
00160
00161 // Multiply the force magnitude to the
00162 // unit direction vector to get the
00163 // resulting force
00164 f *= fval;
00165
00166 // Add the force to the force vectors
00167 forces[i].x += f.x;
00168 forces[i].y += f.y;
00169 forces[i].z += f.z;
00170 }
00171
00172 // Check to see if there is a second radius
00173 // and if we are outside of it
00174 if (twoForces && (dist > r2) )
00175 {
00176 // Assign the force vector to the
00177 // unit direction vector
00178 f.x = diff.x;
00179 f.y = diff.y;
00180 f.z = diff.z;
00181
00182 // Calculate the energy which is
00183 // e = k2*(r_i-r_center2)^exp2
00184 eval = k2;
00185 rval = fabs(dist - r2);
00186
00187 for (j=0; j<exp2; j++)
00188 {
00189 eval *= rval;
00190 }
00191
00192 energy += eval;
00193
00194 // Now calculate the force which is
00195 // e = -k2*exp2*(r_i-r_center2)^(exp2-1)
00196 fval = -exp2*k2;
00197
00198 for (j=0; j<exp2-1; j++)
00199 {
00200 fval *= rval;
00201 }
00202
00203 // Multiply the force magnitude to the
00204 // unit direction vector to get the
00205 // resulting force
00206 f *= fval;
00207
00208 // Add the force to the force vectors
00209 forces[i].x += f.x;
00210 forces[i].y += f.y;
00211 forces[i].z += f.z;
00212 }
00213 }
00214 }
00215 // Loop through and check each atom
00216 for (i=0; i<numAtoms; i++)
00217 {
00218 // Calculate the vector from the atom to the center of the
00219 // cylinder
00220 diff.x = ( axis != 'x' ? 0.0 : x[i].position.x - center.x );
00221 diff.y = ( axis != 'y' ? 0.0 : x[i].position.y - center.y );
00222 diff.z = ( axis != 'z' ? 0.0 : x[i].position.z - center.z );
00223
00224 // Calculate the distance squared
00225 dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
00226
00227 // Look to see if we are outside either radius
00228 if ( (dist_2 > l1_2) || (twoForces && (dist_2 > l2_2)) )
00229 {
00230 // Calculate the distance to the center
00231 dist = sqrt(dist_2);
00232
00233 // Normalize the direction vector
00234 diff /= dist;
00235
00236 // Check to see if we are outside radius 1
00237 if (dist > l1)
00238 {
00239 //printf ("Do faces one force z=%f\n", dist);
00240 // Assign the force vector to the
00241 // unit direction vector
00242 f.x = diff.x;
00243 f.y = diff.y;
00244 f.z = diff.z;
00245
00246 // Calculate the energy which is
00247 // e = k1*(r_i-r_center)^exp1
00248 eval = k1;
00249 rval = fabs(dist - l1);
00250
00251 for (j=0; j<exp1; j++)
00252 {
00253 eval *= rval;
00254 }
00255
00256 energy += eval;
00257
00258 // Now calculate the force which is
00259 // e = -k1*exp1*(r_i-r_center1)^(exp1-1)
00260 fval = -exp1*k1;
00261
00262 for (j=0; j<exp1-1; j++)
00263 {
00264 fval *= rval;
00265 }
00266
00267 // Multiply the force magnitude to the
00268 // unit direction vector to get the
00269 // resulting force
00270 f *= fval;
00271
00272 // Add the force to the force vectors
00273 forces[i].x += f.x;
00274 forces[i].y += f.y;
00275 forces[i].z += f.z;
00276 }
00277
00278 // Check to see if there is a second radius
00279 // and if we are outside of it
00280 if (twoForces && (dist > l2) )
00281 {
00282 //printf ("Do faces two force z=%f\n", dist);
00283 // Assign the force vector to the
00284 // unit direction vector
00285 f.x = diff.x;
00286 f.y = diff.y;
00287 f.z = diff.z;
00288
00289 // Calculate the energy which is
00290 // e = k2*(r_i-r_center2)^exp2
00291 eval = k2;
00292 rval = fabs(dist - l2);
00293
00294 for (j=0; j<exp2; j++)
00295 {
00296 eval *= rval;
00297 }
00298
00299 energy += eval;
00300
00301 // Now calculate the force which is
00302 // e = -k2*exp2*(r_i-r_center2)^(exp2-1)
00303 fval = -exp2*k2;
00304
00305 for (j=0; j<exp2-1; j++)
00306 {
00307 fval *= rval;
00308 }
00309
00310 // Multiply the force magnitude to the
00311 // unit direction vector to get the
00312 // resulting force
00313 f *= fval;
00314
00315 // Add the force to the force vectors
00316 forces[i].x += f.x;
00317 forces[i].y += f.y;
00318 forces[i].z += f.z;
00319 }
00320 }
00321 }
00322
00323 reduction->item(REDUCTION_BC_ENERGY) += energy;
00324 reduction->submit();
00325
00326 }
|
|
|
Definition at line 37 of file ComputeCylindricalBC.h. Referenced by ComputeCylindricalBC(), and doForce(). |
1.3.9.1