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 
00025 // static initialization
00026 int AngleElem::pressureProfileSlabs = 0;
00027 int AngleElem::pressureProfileAtomTypes = 1;
00028 BigReal AngleElem::pressureProfileThickness = 0;
00029 BigReal AngleElem::pressureProfileMin = 0;
00030 
00031 void AngleElem::getMoleculePointers
00032     (Molecule* mol, int* count, int32*** byatom, Angle** structarray)
00033 {
00034 #ifdef MEM_OPT_VERSION
00035   NAMD_die("Should not be called in AngleElem::getMoleculePointers in memory optimized version!");
00036 #else
00037   *count = mol->numAngles;
00038   *byatom = mol->anglesByAtom;
00039   *structarray = mol->angles;
00040 #endif
00041 }
00042 
00043 void AngleElem::getParameterPointers(Parameters *p, const AngleValue **v) {
00044   *v = p->angle_array;
00045 }
00046 
00047 void AngleElem::computeForce(AngleElem *tuples, int ntuple, BigReal *reduction, BigReal *pressureProfileData)
00048 {
00049  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00050 
00051  //fepb BKR
00052  SimParameters *const simParams = Node::Object()->simParameters;
00053  const int step = tuples[0].p[0]->p->flags.step;
00054  const BigReal alchLambda = simParams->getCurrentLambda(step);
00055  const BigReal alchLambda2 = simParams->alchLambda2;
00056  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00057  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00058  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00059  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00060  Molecule *const mol = Node::Object()->molecule;
00061  //fepe
00062 
00063  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00064   const AngleElem &tup = tuples[ituple];
00065   enum { size = 3 };
00066   const AtomID (&atomID)[size](tup.atomID);
00067   const int    (&localIndex)[size](tup.localIndex);
00068   TuplePatchElem * const(&p)[size](tup.p);
00069   const Real (&scale)(tup.scale);
00070   const AngleValue * const(&value)(tup.value);
00071 
00072   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00073                << localIndex[1] << " " << localIndex[2] << std::endl);
00074 
00075   const Position & pos1 = p[0]->x[localIndex[0]].position;
00076   const Position & pos2 = p[1]->x[localIndex[1]].position;
00077   Vector r12 = lattice.delta(pos1,pos2);
00078   register BigReal d12inv = r12.rlength();
00079   const Position & pos3 = p[2]->x[localIndex[2]].position;
00080   Vector r32 = lattice.delta(pos3,pos2);
00081   register BigReal d32inv = r32.rlength();
00082   int normal = value->normal;
00083 
00084   BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
00085   //  Make sure that the cosine value is acceptable.  With roundoff, you
00086   //  can get values like 1.0+2e-16, which makes acos puke.  So instead,
00087   //  just set these kinds of values to exactly 1.0
00088   if (cos_theta > 1.0) cos_theta = 1.0;
00089   else if (cos_theta < -1.0) cos_theta = -1.0;
00090 
00091   BigReal k = value->k * scale;
00092   BigReal theta0 = value->theta0;
00093 
00094   //  Get theta
00095   BigReal theta = acos(cos_theta);
00096 
00097   //  Compare it to the rest angle
00098   BigReal diff;
00099 
00100   if (normal == 1) {
00101     diff = theta - theta0;
00102   } else {
00103     diff = cos_theta - cos(theta0);
00104   }
00105   
00106   //  Add the energy from this angle to the total energy
00107   BigReal energy = k *diff*diff;
00108 
00109   //  Normalize vector r12 and r32
00110   //BigReal d12inv = 1. / d12;
00111   //BigReal d32inv = 1. / d32;
00112 
00113 
00114   //  Calculate constant factor 2k(theta-theta0)/sin(theta)
00115   BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00116   if (normal != 1) {
00117     diff *= (2.0* k);
00118   } else if ( sin_theta < 1.e-6 ) {
00119     // catch case where bonds are parallel
00120     // this gets the force approximately right for theta0 of 0 or pi
00121     // and at least avoids small division for other values
00122     if ( diff < 0. ) diff = 2.0 * k;
00123     else diff = -2.0 * k;
00124   } else { 
00125     diff *= (-2.0* k) / sin_theta; 
00126   }
00127   BigReal c1 = diff * d12inv;
00128   BigReal c2 = diff * d32inv;
00129 
00130   //  Calculate the actual forces
00131   Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00132   Force force2 = force1;
00133   Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00134   force2 += force3;  force2 *= -1;
00135 
00136   //  Check to see if we need to do the Urey-Bradley term
00137   if (value->k_ub)
00138   {
00139         //  Non-zero k_ub value, so calculate the harmonic
00140         //  potential between the 1-3 atoms
00141 
00142   if (normal != 1) {
00143     NAMD_die("ERROR: Can't use cosAngles with Urey-Bradley angles");
00144   }
00145         BigReal k_ub = value->k_ub;
00146         BigReal r_ub = value->r_ub;
00147         Vector r13 = r12 - r32;
00148         BigReal d13 = r13.length();
00149         diff = d13- r_ub;
00150 
00151         energy += k_ub *diff*diff;
00152 
00153         diff *= -2.0*k_ub / d13;
00154         r13 *= diff;
00155 
00156         force1 += r13;
00157         force3 -= r13;
00158   }
00159 
00160   //fepb - BKR scaling of alchemical bonded terms
00161   //       NB: TI derivative is the _unscaled_ energy.
00162   if ( simParams->alchOn ) {
00163     switch ( mol->get_fep_bonded_type(atomID, 3) ) {
00164     case 1:
00165       reduction[angleEnergyIndex_ti_1] += energy;
00166       reduction[angleEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*energy; 
00167       energy *= bond_lambda_1;
00168       force1 *= bond_lambda_1;
00169       force2 *= bond_lambda_1;
00170       force3 *= bond_lambda_1;
00171       break;
00172     case 2:
00173       reduction[angleEnergyIndex_ti_2] += energy;
00174       reduction[angleEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*energy;
00175       energy *= bond_lambda_2;
00176       force1 *= bond_lambda_2;
00177       force2 *= bond_lambda_2;
00178       force3 *= bond_lambda_2;
00179       break;
00180     }
00181   }
00182   //fepe
00183 
00184   p[0]->f[localIndex[0]].x += force1.x;
00185   p[0]->f[localIndex[0]].y += force1.y;
00186   p[0]->f[localIndex[0]].z += force1.z;
00187 
00188   p[1]->f[localIndex[1]].x += force2.x;
00189   p[1]->f[localIndex[1]].y += force2.y;
00190   p[1]->f[localIndex[1]].z += force2.z;
00191 
00192   p[2]->f[localIndex[2]].x += force3.x;
00193   p[2]->f[localIndex[2]].y += force3.y;
00194   p[2]->f[localIndex[2]].z += force3.z;
00195 
00196   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00197 
00198   reduction[angleEnergyIndex] += energy;
00199   reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
00200   reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
00201   reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
00202   reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
00203   reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
00204   reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
00205   reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
00206   reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
00207   reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
00208 
00209   if (pressureProfileData) {
00210     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00211     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00212     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00213     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00214     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00215     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00216     pp_clamp(n1, pressureProfileSlabs);
00217     pp_clamp(n2, pressureProfileSlabs);
00218     pp_clamp(n3, pressureProfileSlabs);
00219     int p1 = p[0]->x[localIndex[0]].partition;
00220     int p2 = p[1]->x[localIndex[1]].partition;
00221     int p3 = p[2]->x[localIndex[2]].partition;
00222     int pn = pressureProfileAtomTypes;
00223     pp_reduction(pressureProfileSlabs, n1, n2, 
00224                 p1, p2, pn,
00225                 force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
00226                 pressureProfileData);
00227     pp_reduction(pressureProfileSlabs, n3, n2, 
00228                 p3, p2, pn,
00229                 force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
00230                 pressureProfileData);
00231   }
00232 
00233  }
00234 }
00235 
00236 
00237 void AngleElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00238 {
00239   reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
00240   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[angleEnergyIndex_f];
00241   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[angleEnergyIndex_ti_1];
00242   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[angleEnergyIndex_ti_2];
00243   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00244 }
00245 

Generated on Thu Nov 23 01:17:10 2017 for NAMD by  doxygen 1.4.7