ComputeRestraints.C

Go to the documentation of this file.
00001 
00007 #include "ComputeRestraints.h"
00008 #include "Node.h"
00009 #include "Molecule.h"
00010 #include "SimParameters.h"
00011 #include "Patch.h"
00012 #include "NamdOneTools.h"
00013 
00014 /************************************************************************/
00015 /*                                                                      */
00016 /*                      FUNCTION ComputeRestraints                      */
00017 /*                                                                      */
00018 /************************************************************************/
00019 
00020 ComputeRestraints::ComputeRestraints(ComputeID c, PatchID pid)
00021   : ComputeHomePatch(c,pid)
00022 {
00023         reduction = ReductionMgr::Object()->willSubmit(REDUCTIONS_BASIC);
00024 
00025         SimParameters *simParams = Node::Object()->simParameters;
00026 
00027         //  Get parameters from the SimParameters object
00028         consExp = simParams->constraintExp;
00029         
00030         //****** BEGIN selective restraints (X,Y,Z) changes 
00031         consSelectOn = simParams->selectConstraintsOn;
00032         if (consSelectOn) {
00033           consSelectX = simParams->constrXOn;
00034           consSelectY = simParams->constrYOn;
00035           consSelectZ = simParams->constrZOn;
00036         }
00037         //****** END selective restraints (X,Y,Z) changes 
00038         
00039         //****** BEGIN moving constraints changes 
00040         consMoveOn = simParams->movingConstraintsOn;
00041         if (consMoveOn) {
00042           moveVel = simParams->movingConsVel;
00043         }
00044         //****** END moving constraints changes 
00045         //****** BEGIN rotating constraints changes 
00046         consRotOn = simParams->rotConstraintsOn;
00047         if (consRotOn) {
00048           rotVel = simParams->rotConsVel;
00049           rotAxis = simParams->rotConsAxis;
00050           rotPivot = simParams->rotConsPivot;
00051         }
00052         //****** END rotating constraints changes 
00053 
00054 }
00055 /*                      END OF FUNCTION ComputeRestraints               */
00056 
00057 /************************************************************************/
00058 /*                                                                      */
00059 /*                      FUNCTION ~ComputeRestraints                     */
00060 /*                                                                      */
00061 /*      This is the destructor for the ComputeRestraints force object.  */
00062 /*                                                                      */
00063 /************************************************************************/
00064 
00065 ComputeRestraints::~ComputeRestraints()
00066 
00067 {
00068         delete reduction;
00069 }
00070 /*                      END OF FUNCTION ~ComputeRestraints              */
00071 
00072 /************************************************************************/
00073 /*                                                                      */
00074 /*                              FUNCTION force                          */
00075 /*                                                                      */
00076 /************************************************************************/
00077 void ComputeRestraints::doForce(FullAtom* p, Results* res)
00078 {
00079         Molecule *molecule = Node::Object()->molecule;
00080         Real k;                 //  Force constant
00081         SimParameters *simParams = Node::Object()->simParameters;
00082         BigReal scaling = simParams->constraintScaling;
00083         Vector refPos;          //  Reference position
00084         BigReal r, r2;  //  r=distance between atom position and the
00085                         //  reference position, r2 = r^2
00086         Vector Rij;     //  vector between current position and reference pos
00087         BigReal value;  //  Current calculated value
00088 
00089         // aliases to work with old code
00090         Force *f = res->f[Results::normal];
00091         BigReal energy = 0;
00092         BigReal m[9];
00093         Tensor virial;
00094         Vector netForce = 0;
00095 
00096         // BEGIN moving and rotating constraint changes ******
00097 
00098         // This version only allows one atom to be moved 
00099         // and only ALL ref positions to be rotated
00100 
00101         int currentTime = patch->flags.step;
00102         if (consRotOn) {
00103           vec_rotation_matrix(rotVel * currentTime, rotAxis, m);
00104         }
00105 
00106         // END moving and rotating constraint changes ******
00107 
00108           
00109         if (scaling != 0.) for (int localID=0; localID<numAtoms; ++localID)
00110         {
00111           if (molecule->is_atom_constrained(p[localID].id))
00112           {
00113 #ifndef MEM_OPT_VERSION
00114             molecule->get_cons_params(k, refPos, p[localID].id);
00115 #else
00116             k = 1.0;
00117             refPos = p[localID].fixedPosition;
00118 #endif
00119 
00120             k *= scaling;
00121 
00122             // BEGIN moving and rotating constraint changes ******
00123             
00124             if (consMoveOn) {
00125               refPos += currentTime * moveVel;
00126             }
00127             else if(consRotOn) {
00128               refPos = mat_multiply_vec(refPos - rotPivot, m) + rotPivot;
00129             }
00130 
00131             // END moving and rotating constraint changes *******
00132 
00133             if (simParams->sphericalConstraintsOn) {
00134               BigReal refRad = (refPos - simParams->sphericalConstrCenter).length();
00135               Vector relPos = patch->lattice.delta(p[localID].position, simParams->sphericalConstrCenter);
00136               refPos = simParams->sphericalConstrCenter + relPos * (refRad/relPos.length());
00137             }
00138 
00139             Rij = patch->lattice.delta(refPos,p[localID].position);
00140             Vector vpos = refPos - Rij;
00141 
00142             //****** BEGIN selective restraints (X,Y,Z) changes 
00143             if (consSelectOn) { // check which components we want to restrain:
00144               if (!consSelectX) {Rij.x=0.0;}  // by setting the appropriate
00145               if (!consSelectY) {Rij.y=0.0;}  // Cartesian component to zero
00146               if (!consSelectZ) {Rij.z=0.0;}  // we will avoid a restoring force
00147             }                                 // along that component, and the
00148                                               // energy will also be correct
00149             //****** END selective restraints (X,Y,Z) changes 
00150 
00151             
00152             //  Calculate the distance and the distance squared
00153             r2 = Rij.length2();
00154             r = sqrt(r2);
00155 
00156             //  Only calculate the energy and the force if the distance is
00157             //  non-zero.   Otherwise, you could end up dividing by 0, which
00158             //  is bad
00159             if (r>0.0)
00160             {
00161               value=k;
00162       
00163               //  Loop through and multiple k by r consExp times.
00164               //  i.e.  calculate kr^e
00165               //  I know, I could use pow(), but I don't trust it.
00166               for (int k=0; k<consExp; ++k)
00167               {
00168                 value *= r;
00169               }
00170 
00171               //  Add to the energy total
00172               energy += value;
00173       
00174               //  Now calculate the force, which is ekr^(e-1).  Also, divide
00175               //  by another factor of r to normalize the vector before we
00176               //  multiple it by the magnitude
00177               value *= consExp;
00178               value /= r2;
00179       
00180               Rij *= value;
00181               //iout << iINFO << "restraining force" << Rij        << "\n"; 
00182 
00183               f[localID] += Rij;
00184               netForce += Rij;
00185               virial += outer(Rij,vpos);
00186             }
00187           }
00188         }
00189 
00190         reduction->item(REDUCTION_BC_ENERGY) += energy;
00191         ADD_TENSOR_OBJECT(reduction,REDUCTION_VIRIAL_NORMAL,virial);
00192         ADD_VECTOR_OBJECT(reduction,REDUCTION_EXT_FORCE_NORMAL,netForce);
00193         reduction->submit();
00194 
00195 }
00196 /*                      END OF FUNCTION force                           */
00197 

Generated on Tue Sep 19 01:17:11 2017 for NAMD by  doxygen 1.4.7