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

ComputeAngles.C

Go to the documentation of this file.
00001 
00007 /*
00008    Methods for ComputeAngles.  Main code is for
00009    loading in the AngleElem information and
00010    for computing forces and energies for all angles on node's.
00011    HomePatch(es)
00012 */
00013 
00014 #include "InfoStream.h"
00015 #include "ComputeAngles.h"
00016 #include "Molecule.h"
00017 #include "Parameters.h"
00018 #include "Node.h"
00019 #include "ReductionMgr.h"
00020 #include "Lattice.h"
00021 #include "PressureProfile.h"
00022 #include "Debug.h"
00023 
00024 #if 0
00025 void AngleElem::loadTuplesForAtom
00026   (void *voidlist, AtomID atomID, Molecule *molecule)
00027 {
00028       DebugM(1, "::loadTuplesForAtom - atomID " << atomID << std::endl );
00029       UniqueSet<AngleElem> &angleList =
00030                   *( (UniqueSet<AngleElem>*) voidlist );
00031 
00032       DebugM(1, "::loadTuplesForAtom - current list size " << angleList.size() << std::endl );
00033 
00034       /* get list of all angles for the atom */
00035       int *angles = molecule->get_angles_for_atom(atomID);
00036       DebugM(1, "::loadTuplesForAtom - atomID " << atomID << std::endl );
00037 
00038       /* cycle through each angle */
00039       int angleNum = *angles;
00040       while(angleNum != -1)
00041       {
00042         /* store angle in the list */
00043         DebugM(1, "::loadTuplesForAtom - loading angle " << angleNum << std::endl );
00044         angleList.add(AngleElem(molecule->get_angle(angleNum)));
00045         angleNum = *(++angles);
00046       }
00047 }
00048 #endif
00049 
00050 // static initialization
00051 int AngleElem::pressureProfileSlabs = 0;
00052 int AngleElem::pressureProfileAtomTypes = 1;
00053 BigReal AngleElem::pressureProfileThickness = 0;
00054 BigReal AngleElem::pressureProfileMin = 0;
00055 
00056 void AngleElem::getMoleculePointers
00057     (Molecule* mol, int* count, int32*** byatom, Angle** structarray)
00058 {
00059 #ifdef MEM_OPT_VERSION
00060   NAMD_die("Should not be called in AngleElem::getMoleculePointers in memory optimized version!");
00061 #else
00062   *count = mol->numAngles;
00063   *byatom = mol->anglesByAtom;
00064   *structarray = mol->angles;
00065 #endif
00066 }
00067 
00068 void AngleElem::getParameterPointers(Parameters *p, const AngleValue **v) {
00069   *v = p->angle_array;
00070 }
00071 
00072 void AngleElem::computeForce(BigReal *reduction, BigReal *pressureProfileData)
00073 {
00074   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00075                << localIndex[1] << " " << localIndex[2] << std::endl);
00076 
00077   const Position & pos1 = p[0]->x[localIndex[0]].position;
00078   const Lattice & lattice = p[0]->p->lattice;
00079   const Position & pos2 = p[1]->x[localIndex[1]].position;
00080   Vector r12 = lattice.delta(pos1,pos2);
00081   register BigReal d12inv = r12.rlength();
00082   const Position & pos3 = p[2]->x[localIndex[2]].position;
00083   Vector r32 = lattice.delta(pos3,pos2);
00084   register BigReal d32inv = r32.rlength();
00085   int normal = value->normal;
00086 
00087   BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
00088   //  Make sure that the cosine value is acceptable.  With roundoff, you
00089   //  can get values like 1.0+2e-16, which makes acos puke.  So instead,
00090   //  just set these kinds of values to exactly 1.0
00091   if (cos_theta > 1.0) cos_theta = 1.0;
00092   else if (cos_theta < -1.0) cos_theta = -1.0;
00093 
00094   BigReal k = value->k * scale;
00095   BigReal theta0 = value->theta0;
00096 
00097   //  Get theta
00098   BigReal theta = acos(cos_theta);
00099 
00100   //  Compare it to the rest angle
00101   BigReal diff;
00102 
00103   if (normal == 1) {
00104     diff = theta - theta0;
00105   } else {
00106     diff = cos_theta - cos(theta0);
00107   }
00108   
00109   //  Add the energy from this angle to the total energy
00110   BigReal energy = k *diff*diff;
00111 
00112   //  Normalize vector r12 and r32
00113   //BigReal d12inv = 1. / d12;
00114   //BigReal d32inv = 1. / d32;
00115 
00116 
00117   //  Calculate constant factor 2k(theta-theta0)/sin(theta)
00118   BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00119   if (normal != 1) {
00120     diff *= (2.0* k);
00121   } else if ( sin_theta < 1.e-6 ) {
00122     // catch case where bonds are parallel
00123     // this gets the force approximately right for theta0 of 0 or pi
00124     // and at least avoids small division for other values
00125     if ( diff < 0. ) diff = 2.0 * k;
00126     else diff = -2.0 * k;
00127   } else { 
00128     diff *= (-2.0* k) / sin_theta; 
00129   }
00130   BigReal c1 = diff * d12inv;
00131   BigReal c2 = diff * d32inv;
00132 
00133   //  Calculate the actual forces
00134   Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00135   Force force2 = force1;
00136   Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00137   force2 += force3;  force2 *= -1;
00138 
00139   //  Check to see if we need to do the Urey-Bradley term
00140   if (value->k_ub)
00141   {
00142         //  Non-zero k_ub value, so calculate the harmonic
00143         //  potential between the 1-3 atoms
00144 
00145   if (normal != 1) {
00146     NAMD_die("ERROR: Can't use cosAngles with Urey-Bradley angles");
00147   }
00148         BigReal k_ub = value->k_ub;
00149         BigReal r_ub = value->r_ub;
00150         Vector r13 = r12 - r32;
00151         BigReal d13 = r13.length();
00152         diff = d13- r_ub;
00153 
00154         energy += k_ub *diff*diff;
00155 
00156         diff *= -2.0*k_ub / d13;
00157         r13 *= diff;
00158 
00159         force1 += r13;
00160         force3 -= r13;
00161   }
00162 
00163   p[0]->f[localIndex[0]].x += force1.x;
00164   p[0]->f[localIndex[0]].y += force1.y;
00165   p[0]->f[localIndex[0]].z += force1.z;
00166 
00167   p[1]->f[localIndex[1]].x += force2.x;
00168   p[1]->f[localIndex[1]].y += force2.y;
00169   p[1]->f[localIndex[1]].z += force2.z;
00170 
00171   p[2]->f[localIndex[2]].x += force3.x;
00172   p[2]->f[localIndex[2]].y += force3.y;
00173   p[2]->f[localIndex[2]].z += force3.z;
00174 
00175   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00176 
00177   reduction[angleEnergyIndex] += energy;
00178   reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
00179   reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
00180   reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
00181   reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
00182   reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
00183   reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
00184   reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
00185   reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
00186   reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
00187 
00188   if (pressureProfileData) {
00189     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00190     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00191     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00192     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00193     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00194     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00195     pp_clamp(n1, pressureProfileSlabs);
00196     pp_clamp(n2, pressureProfileSlabs);
00197     pp_clamp(n3, pressureProfileSlabs);
00198     int p1 = p[0]->x[localIndex[0]].partition;
00199     int p2 = p[1]->x[localIndex[1]].partition;
00200     int p3 = p[2]->x[localIndex[2]].partition;
00201     int pn = pressureProfileAtomTypes;
00202     pp_reduction(pressureProfileSlabs, n1, n2, 
00203                 p1, p2, pn,
00204                 force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
00205                 pressureProfileData);
00206     pp_reduction(pressureProfileSlabs, n3, n2, 
00207                 p3, p2, pn,
00208                 force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
00209                 pressureProfileData);
00210   }
00211 }
00212 
00213 
00214 void AngleElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00215 {
00216   reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
00217   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00218 }
00219 

Generated on Tue Nov 24 04:07:42 2009 for NAMD by  doxygen 1.3.9.1