ComputeSphericalBC.C

Go to the documentation of this file.
00001 
00007 #include "ComputeSphericalBC.h"
00008 #include "Node.h"
00009 #include "SimParameters.h"
00010 #include "Patch.h"
00011 
00012 /************************************************************************/
00013 /*                                                                      */
00014 /*                      FUNCTION ComputeSphericalBC                     */
00015 /*                                                                      */
00016 /*      This is the constructor for the ComputeSphericalBC 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 spherical       */
00021 /*   boundaries.                                                        */
00022 /*                                                                      */
00023 /************************************************************************/
00024 
00025 ComputeSphericalBC::ComputeSphericalBC(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         r1 = simParams->sphericalBCr1;
00034         r2 = simParams->sphericalBCr2;
00035         r1_2 = r1*r1;
00036         r2_2 = r2*r2;
00037         k1 = simParams->sphericalBCk1;
00038         k2 = simParams->sphericalBCk2;
00039         exp1 = simParams->sphericalBCexp1;
00040         exp2 = simParams->sphericalBCexp2;
00041 
00042         //  Check to see if this is one set of parameters or two
00043         if (r2 > -1.0)
00044         {
00045                 twoForces = TRUE;
00046         }
00047         else
00048         {
00049                 twoForces = FALSE;
00050         }
00051 
00052         center = simParams->sphericalCenter;
00053 
00054 }
00055 /*                      END OF FUNCTION ComputeSphericalBC              */
00056 
00057 /************************************************************************/
00058 /*                                                                      */
00059 /*                      FUNCTION ~ComputeSphericalBC                    */
00060 /*                                                                      */
00061 /*      This is the destructor for the ComputeSphericalBC force object. */
00062 /*   It currently does ABSOLUTELY NOTHING!!                             */
00063 /*                                                                      */
00064 /************************************************************************/
00065 
00066 ComputeSphericalBC::~ComputeSphericalBC()
00067 
00068 {
00069         delete reduction;
00070 }
00071 /*                      END OF FUNCTION ~ComputeSphericalBC             */
00072 
00073 /************************************************************************/
00074 /*                                                                      */
00075 /*                              FUNCTION force                          */
00076 /*                                                                      */
00077 /*   INPUTS:                                                            */
00078 /*      numAtoms - Number of coordinates being passed                   */
00079 /*      x - Array of atom coordinates                                   */
00080 /*      forces - Array of force vectors                                 */
00081 /*                                                                      */
00082 /*      This function calculates the force and energy for the shereical */
00083 /*   boundary conditions for this patch.                                */
00084 /*                                                                      */
00085 /************************************************************************/
00086 void ComputeSphericalBC::doForce(FullAtom* p, Results* r)
00087 {
00088         Vector diff;            //  Distance from atom to center of sphere
00089         Vector f;               //  Calculated force vector
00090         int i, j;               //  Loop counters
00091         BigReal dist, dist_2;   //  Distance from atom to center, and this
00092                                 //  distance squared
00093         BigReal rval;           //  Difference between distance from atom
00094                                 //  to center and radius of sphere
00095         BigReal eval;           //  Energy value for this atom
00096         BigReal fval;           //  Force magnitude for this atom
00097 
00098         // aliases to work with old code
00099         FullAtom *x = p;
00100         Force *forces = r->f[Results::normal];
00101         BigReal energy = 0;
00102 
00103         //  Loop through and check each atom
00104         for (i=0; i<numAtoms; i++)
00105         {
00106                 //  Calculate the vector from the atom to the center of the
00107                 //  sphere
00108                 diff.x = x[i].position.x - center.x;
00109                 diff.y = x[i].position.y - center.y;
00110                 diff.z = x[i].position.z - center.z;
00111                 
00112                 //  Calculate the distance squared
00113                 dist_2 = diff.x*diff.x + diff.y*diff.y + diff.z*diff.z;
00114 
00115                 //  Look to see if we are outside either radius
00116                 if ( (dist_2 > r1_2) || (twoForces && (dist_2 > r2_2)) )
00117                 {
00118                         //  Calculate the distance to the center
00119                         dist = sqrt(dist_2);
00120 
00121                         //  Normalize the direction vector
00122                         diff /= dist;
00123 
00124                         //  Check to see if we are outside radius 1
00125                         if (dist > r1)
00126                         {
00127                                 //  Assign the force vector to the
00128                                 //  unit direction vector
00129                                 f.x = diff.x;
00130                                 f.y = diff.y;
00131                                 f.z = diff.z;
00132 
00133                                 //  Calculate the energy which is
00134                                 //  e = k1*(r_i-r_center)^exp1
00135                                 eval = k1;
00136                                 rval = fabs(dist - r1);
00137 
00138                                 for (j=0; j<exp1; j++)
00139                                 {
00140                                         eval *= rval;
00141                                 }
00142 
00143                                 energy += eval;
00144 
00145                                 //  Now calculate the force which is
00146                                 //  e = -k1*exp1*(r_i-r_center1)^(exp1-1)
00147                                 fval = -exp1*k1;
00148 
00149                                 for (j=0; j<exp1-1; j++)
00150                                 {
00151                                         fval *= rval;
00152                                 }
00153 
00154                                 //  Multiply the force magnitude to the
00155                                 //  unit direction vector to get the
00156                                 //  resulting force
00157                                 f *= fval;
00158 
00159                                 //  Add the force to the force vectors
00160                                 forces[i].x += f.x;
00161                                 forces[i].y += f.y;
00162                                 forces[i].z += f.z;
00163                         }
00164 
00165                         //  Check to see if there is a second radius
00166                         //  and if we are outside of it
00167                         if (twoForces && (dist > r2) )
00168                         {
00169                                 //  Assign the force vector to the
00170                                 //  unit direction vector
00171                                 f.x = diff.x;
00172                                 f.y = diff.y;
00173                                 f.z = diff.z;
00174 
00175                                 //  Calculate the energy which is
00176                                 //  e = k2*(r_i-r_center2)^exp2
00177                                 eval = k2;
00178                                 rval = fabs(dist - r2);
00179 
00180                                 for (j=0; j<exp2; j++)
00181                                 {
00182                                         eval *= rval;
00183                                 }
00184 
00185                                 energy += eval;
00186 
00187                                 //  Now calculate the force which is
00188                                 //  e = -k2*exp2*(r_i-r_center2)^(exp2-1)
00189                                 fval = -exp2*k2;
00190 
00191                                 for (j=0; j<exp2-1; j++)
00192                                 {
00193                                         fval *= rval;
00194                                 }
00195 
00196                                 //  Multiply the force magnitude to the
00197                                 //  unit direction vector to get the
00198                                 //  resulting force
00199                                 f *= fval;
00200 
00201                                 //  Add the force to the force vectors
00202                                 forces[i].x += f.x;
00203                                 forces[i].y += f.y;
00204                                 forces[i].z += f.z;
00205                         }
00206                 }
00207         }
00208 
00209     reduction->item(REDUCTION_BC_ENERGY) += energy;
00210     reduction->submit();
00211 
00212 }
00213 /*                      END OF FUNCTION force                           */

Generated on Sun Nov 19 01:17:12 2017 for NAMD by  doxygen 1.4.7