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   : 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 }
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 void ComputeCylindricalBC::doForce(FullAtom* p, Results* r)
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 }
00327 /*                      END OF FUNCTION force                           */
00328 

Generated on Mon Nov 20 01:17:10 2017 for NAMD by  doxygen 1.4.7