Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | File List | Class Members | File Members

ComputeCylindricalBC.C

Go to the documentation of this file.
00001 
00007 #include "ComputeCylindricalBC.h"
00008 #include "Node.h"
00009 #include "SimParameters.h"
00010 #include "Patch.h"
00011 
00012 /************************************************************************/
00013 /*                                                                      */
00014 /*                      FUNCTION ComputeCylindricalBC                   */
00015 /*                                                                      */
00016 /*      This is the constructor for the ComputeCylindricalBC force object.*/
00017 /*   It is responsible for getting all the parameters from the          */
00018 /*   SimParameters object and then determining if this object needs     */
00019 /*   to perform any computation.  It only needs to do so if there is    */
00020 /*   some portion of the patch that lays outside of the cylindrical     */
00021 /*   boundaries.                                                        */
00022 /*                                                                      */
00023 /************************************************************************/
00024 
00025 ComputeCylindricalBC::ComputeCylindricalBC(ComputeID c, PatchID pid)
00026   : ComputePatch(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 }
00061 /*                      END OF FUNCTION ComputeCylindricalBC            */
00062 
00063 /************************************************************************/
00064 /*                                                                      */
00065 /*                      FUNCTION ~ComputeCylindricalBC                  */
00066 /*                                                                      */
00067 /*      This is the destructor for the ComputeCylindricalBC force object.       */
00068 /*   It currently does ABSOLUTELY NOTHING!!                             */
00069 /*                                                                      */
00070 /************************************************************************/
00071 
00072 ComputeCylindricalBC::~ComputeCylindricalBC()
00073 
00074 {
00075         delete reduction;
00076 }
00077 /*                      END OF FUNCTION ~ComputeCylindricalBC           */
00078 
00079 /************************************************************************/
00080 /*                                                                      */
00081 /*                              FUNCTION force                          */
00082 /*                                                                      */
00083 /*   INPUTS:                                                            */
00084 /*      numAtoms - Number of coordinates being passed                   */
00085 /*      x - Array of atom coordinates                                   */
00086 /*      forces - Array of force vectors                                 */
00087 /*                                                                      */
00088 /*      This function calculates the force and energy for the cylindri. */
00089 /*   boundary conditions for this patch.                                */
00090 /*                                                                      */
00091 /************************************************************************/
00092 
00093 #ifdef MEM_OPT_VERSION
00094 void ComputeCylindricalBC::doForce(CompAtom* p, CompAtomExt* pExt, Results* r)
00095 #else
00096 void ComputeCylindricalBC::doForce(CompAtom* p, Results* r)
00097 #endif
00098 {
00099         Vector diff;            //  Distance from atom to center of cylinder
00100         Vector f;               //  Calculated force vector
00101         int i, j;               //  Loop counters
00102         BigReal dist, dist_2;   //  Distance from atom to center, and this
00103                                 //  distance squared
00104         BigReal rval;           //  Difference between distance from atom
00105                                 //  to center and radius of cylinder
00106         BigReal eval;           //  Energy value for this atom
00107         BigReal fval;           //  Force magnitude for this atom
00108 
00109         // aliases to work with old code
00110         CompAtom *x = p;
00111         Force *forces = r->f[Results::normal];
00112         BigReal energy = 0;
00113 
00114         //  Loop through and check each atom
00115         for (i=0; i<numAtoms; i++)
00116         {
00117                 //  Calculate the vector from the atom to the center of the
00118                 //  cylinder
00119                 diff.x = ( axis == 'x' ? 0.0 : x[i].position.x - center.x );
00120                 diff.y = ( axis == 'y' ? 0.0 : x[i].position.y - center.y );
00121                 diff.z = ( axis == 'z' ? 0.0 : x[i].position.z - center.z );
00122                 
00123                 //  Calculate the distance squared
00124                 dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
00125 
00126                 //  Look to see if we are outside either radius
00127                 if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
00128                 {
00129                         //  Calculate the distance to the center
00130                         dist = sqrt(dist_2);
00131 
00132                         //  Normalize the direction vector
00133                         diff /= dist;
00134 
00135                         //  Check to see if we are outside radius 1
00136                         if (dist > r1)
00137                         {
00138                                 //  Assign the force vector to the
00139                                 //  unit direction vector
00140                                 f.x = diff.x;
00141                                 f.y = diff.y;
00142                                 f.z = diff.z;
00143 
00144                                 //  Calculate the energy which is
00145                                 //  e = k1*(r_i-r_center)^exp1
00146                                 eval = k1;
00147                                 rval = fabs(dist - r1);
00148 
00149                                 for (j=0; j<exp1; j++)
00150                                 {
00151                                         eval *= rval;
00152                                 }
00153 
00154                                 energy += eval;
00155 
00156                                 //  Now calculate the force which is
00157                                 //  e = -k1*exp1*(r_i-r_center1)^(exp1-1)
00158                                 fval = -exp1*k1;
00159 
00160                                 for (j=0; j<exp1-1; j++)
00161                                 {
00162                                         fval *= rval;
00163                                 }
00164 
00165                                 //  Multiply the force magnitude to the
00166                                 //  unit direction vector to get the
00167                                 //  resulting force
00168                                 f *= fval;
00169 
00170                                 //  Add the force to the force vectors
00171                                 forces[i].x += f.x;
00172                                 forces[i].y += f.y;
00173                                 forces[i].z += f.z;
00174                         }
00175 
00176                         //  Check to see if there is a second radius
00177                         //  and if we are outside of it
00178                         if (twoForces && (dist > r2) )
00179                         {
00180                                 //  Assign the force vector to the
00181                                 //  unit direction vector
00182                                 f.x = diff.x;
00183                                 f.y = diff.y;
00184                                 f.z = diff.z;
00185 
00186                                 //  Calculate the energy which is
00187                                 //  e = k2*(r_i-r_center2)^exp2
00188                                 eval = k2;
00189                                 rval = fabs(dist - r2);
00190 
00191                                 for (j=0; j<exp2; j++)
00192                                 {
00193                                         eval *= rval;
00194                                 }
00195 
00196                                 energy += eval;
00197 
00198                                 //  Now calculate the force which is
00199                                 //  e = -k2*exp2*(r_i-r_center2)^(exp2-1)
00200                                 fval = -exp2*k2;
00201 
00202                                 for (j=0; j<exp2-1; j++)
00203                                 {
00204                                         fval *= rval;
00205                                 }
00206 
00207                                 //  Multiply the force magnitude to the
00208                                 //  unit direction vector to get the
00209                                 //  resulting force
00210                                 f *= fval;
00211 
00212                                 //  Add the force to the force vectors
00213                                 forces[i].x += f.x;
00214                                 forces[i].y += f.y;
00215                                 forces[i].z += f.z;
00216                         }
00217                 }
00218         }
00219        //  Loop through and check each atom
00220         for (i=0; i<numAtoms; i++)
00221         {
00222                 //  Calculate the vector from the atom to the center of the
00223                 //  cylinder
00224                 diff.x = ( axis != 'x' ? 0.0 : x[i].position.x - center.x );
00225                 diff.y = ( axis != 'y' ? 0.0 : x[i].position.y - center.y );
00226                 diff.z = ( axis != 'z' ? 0.0 : x[i].position.z - center.z );
00227 
00228                 //  Calculate the distance squared
00229                 dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
00230 
00231                 //  Look to see if we are outside either radius
00232                 if ( (dist_2 > l1_2)  || (twoForces && (dist_2 > l2_2)) )
00233                 {
00234                         //  Calculate the distance to the center
00235                         dist = sqrt(dist_2);
00236 
00237                         //  Normalize the direction vector
00238                         diff /= dist;
00239 
00240                         //  Check to see if we are outside radius 1
00241                         if (dist > l1)
00242                         {
00243 //printf ("Do faces one force z=%f\n", dist);
00244                                 //  Assign the force vector to the
00245                                 //  unit direction vector
00246                                 f.x = diff.x;
00247                                 f.y = diff.y;
00248                                 f.z = diff.z;
00249 
00250                                 //  Calculate the energy which is
00251                                 //  e = k1*(r_i-r_center)^exp1
00252                                 eval = k1;
00253                                 rval = fabs(dist - l1);
00254 
00255                                 for (j=0; j<exp1; j++)
00256                                 {
00257                                         eval *= rval;
00258                                 }
00259 
00260                                 energy += eval;
00261 
00262                                 //  Now calculate the force which is
00263                                 //  e = -k1*exp1*(r_i-r_center1)^(exp1-1)
00264                                 fval = -exp1*k1;
00265 
00266                                 for (j=0; j<exp1-1; j++)
00267                                 {
00268                                         fval *= rval;
00269                                 }
00270 
00271                                 //  Multiply the force magnitude to the
00272                                 //  unit direction vector to get the
00273                                 //  resulting force
00274                                 f *= fval;
00275 
00276                                 //  Add the force to the force vectors
00277                                 forces[i].x += f.x;
00278                                 forces[i].y += f.y;
00279                                 forces[i].z += f.z;
00280                         }
00281 
00282                         //  Check to see if there is a second radius
00283                         //  and if we are outside of it
00284                         if (twoForces && (dist > l2) )
00285                         {
00286 //printf ("Do faces two force z=%f\n", dist);
00287                                 //  Assign the force vector to the
00288                                 //  unit direction vector
00289                                 f.x = diff.x;
00290                                 f.y = diff.y;
00291                                 f.z = diff.z;
00292 
00293                                 //  Calculate the energy which is
00294                                 //  e = k2*(r_i-r_center2)^exp2
00295                                 eval = k2;
00296                                 rval = fabs(dist - l2);
00297 
00298                                 for (j=0; j<exp2; j++)
00299                                 {
00300                                         eval *= rval;
00301                                 }
00302 
00303                                 energy += eval;
00304 
00305                                 //  Now calculate the force which is
00306                                 //  e = -k2*exp2*(r_i-r_center2)^(exp2-1)
00307                                 fval = -exp2*k2;
00308 
00309                                 for (j=0; j<exp2-1; j++)
00310                                 {
00311                                         fval *= rval;
00312                                 }
00313 
00314                                 //  Multiply the force magnitude to the
00315                                 //  unit direction vector to get the
00316                                 //  resulting force
00317                                 f *= fval;
00318 
00319                                 //  Add the force to the force vectors
00320                                 forces[i].x += f.x;
00321                                 forces[i].y += f.y;
00322                                 forces[i].z += f.z;
00323                         }
00324                 }
00325         }
00326 
00327     reduction->item(REDUCTION_BC_ENERGY) += energy;
00328     reduction->submit();
00329 
00330 }
00331 /*                      END OF FUNCTION force                           */
00332 

Generated on Sat Sep 6 04:07:40 2008 for NAMD by  doxygen 1.3.9.1