ComputeDihedrals.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeDihedrals.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "Node.h"
00012 #include "ReductionMgr.h"
00013 #include "Lattice.h"
00014 #include "PressureProfile.h"
00015 #include "Debug.h"
00016 
00017 
00018 // static initialization
00019 int DihedralElem::pressureProfileSlabs = 0;
00020 int DihedralElem::pressureProfileAtomTypes = 1;
00021 BigReal DihedralElem::pressureProfileThickness = 0;
00022 BigReal DihedralElem::pressureProfileMin = 0;
00023 
00024 void DihedralElem::getMoleculePointers
00025     (Molecule* mol, int* count, int32*** byatom, Dihedral** structarray)
00026 {
00027 #ifdef MEM_OPT_VERSION
00028   NAMD_die("Should not be called in DihedralElem::getMoleculePointers in memory optimized version!");
00029 #else
00030   *count = mol->numDihedrals;
00031   *byatom = mol->dihedralsByAtom;
00032   *structarray = mol->dihedrals;
00033 #endif
00034 }
00035 
00036 void DihedralElem::getParameterPointers(Parameters *p, const DihedralValue **v) {
00037   *v = p->dihedral_array;
00038 }
00039 
00040 void DihedralElem::computeForce(DihedralElem *tuples, int ntuple, BigReal *reduction, 
00041                                 BigReal *pressureProfileData)
00042 {
00043  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00044 
00045  //fepb BKR
00046  SimParameters *const simParams = Node::Object()->simParameters;
00047  const int step = tuples[0].p[0]->p->flags.step;
00048  const BigReal alchLambda = simParams->getCurrentLambda(step);
00049  const BigReal alchLambda2 = simParams->alchLambda2;
00050  const BigReal bond_lambda_1 = simParams->getBondLambda(alchLambda);
00051  const BigReal bond_lambda_2 = simParams->getBondLambda(1-alchLambda);
00052  const BigReal bond_lambda_12 = simParams->getBondLambda(alchLambda2);
00053  const BigReal bond_lambda_22 = simParams->getBondLambda(1-alchLambda2);
00054  Molecule *const mol = Node::Object()->molecule;
00055  //fepe
00056 
00057  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00058   const DihedralElem &tup = tuples[ituple];
00059   enum { size = 4 };
00060   const AtomID (&atomID)[size](tup.atomID);
00061   const int    (&localIndex)[size](tup.localIndex);
00062   TuplePatchElem * const(&p)[size](tup.p);
00063   const Real (&scale)(tup.scale);
00064   const DihedralValue * const(&value)(tup.value);
00065 
00066   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00067                << localIndex[1] << " " << localIndex[2] << std::endl);
00068 
00069   //  Calculate the vectors between atoms
00070   const Position & pos0 = p[0]->x[localIndex[0]].position;
00071   const Position & pos1 = p[1]->x[localIndex[1]].position;
00072   const Vector r12 = lattice.delta(pos0,pos1);
00073   const Position & pos2 = p[2]->x[localIndex[2]].position;
00074   const Vector r23 = lattice.delta(pos1,pos2);
00075   const Position & pos3 = p[3]->x[localIndex[3]].position;
00076   const Vector r34 = lattice.delta(pos2,pos3);
00077 
00078   //  Calculate the cross products and distances
00079   Vector A = cross(r12,r23);
00080   register  BigReal rAinv = A.rlength();
00081   Vector B = cross(r23,r34);
00082   register  BigReal rBinv = B.rlength();
00083   Vector C = cross(r23,A);
00084   register  BigReal rCinv = C.rlength();
00085 
00086   //  Calculate the sin and cos
00087   BigReal cos_phi = (A*B)*(rAinv*rBinv);
00088   BigReal sin_phi = (C*B)*(rCinv*rBinv);
00089 
00090   BigReal phi= -atan2(sin_phi,cos_phi);
00091 
00092   BigReal K=0;          // energy
00093   BigReal K1=0;         // force
00094 
00095   // get the dihedral information
00096   int multiplicity = value->multiplicity;
00097 
00098   //  Loop through the multiple parameter sets for this
00099   //  bond.  We will only loop more than once if this
00100   //  has multiple parameter sets from Charmm22
00101   for (int mult_num=0; mult_num<multiplicity; mult_num++)
00102   {
00103     /* get angle information */
00104     Real k = value->values[mult_num].k * scale;
00105     Real delta = value->values[mult_num].delta;
00106     int n = value->values[mult_num].n;
00107 
00108     //  Calculate the energy
00109     if (n)
00110     {
00111       //  Periodicity is greater than 0, so use cos form
00112       K += k*(1+cos(n*phi - delta));
00113       K1 += -n*k*sin(n*phi - delta);
00114     }
00115     else
00116     {
00117       //  Periodicity is 0, so just use the harmonic form
00118       BigReal diff = phi-delta;
00119       if (diff < -PI)           diff += TWOPI;
00120       else if (diff > PI)       diff -= TWOPI;
00121 
00122       K += k*diff*diff;
00123       K1 += 2.0*k*diff;
00124     }
00125   } /* for multiplicity */
00126 
00127   //fepb - BKR scaling of alchemical bonded terms
00128   //       NB: TI derivative is the _unscaled_ energy.
00129   if ( simParams->alchOn ) {
00130     switch ( mol->get_fep_bonded_type(atomID, 4) ) {
00131     case 1:
00132       reduction[dihedralEnergyIndex_ti_1] += K;
00133       reduction[dihedralEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1)*K;
00134       K *= bond_lambda_1;
00135       K1 *= bond_lambda_1;
00136       break;
00137     case 2:
00138       reduction[dihedralEnergyIndex_ti_2] += K;
00139       reduction[dihedralEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2)*K;
00140       K *= bond_lambda_2;
00141       K1 *= bond_lambda_2;
00142       break;
00143     }
00144   }
00145   //fepe
00146 
00147   Force f1,f2,f3;
00148 
00149   //  Normalize B
00150   //rB = 1.0/rB;
00151   B *= rBinv;
00152 
00153     //  Next, we want to calculate the forces.  In order
00154     //  to do that, we first need to figure out whether the
00155     //  sin or cos form will be more stable.  For this,
00156     //  just look at the value of phi
00157     if (fabs(sin_phi) > 0.1)
00158     {
00159       //  use the sin version to avoid 1/cos terms
00160 
00161       //  Normalize A
00162       A *= rAinv;
00163       Vector dcosdA;
00164       Vector dcosdB;
00165 
00166       dcosdA.x = rAinv*(cos_phi*A.x-B.x);
00167       dcosdA.y = rAinv*(cos_phi*A.y-B.y);
00168       dcosdA.z = rAinv*(cos_phi*A.z-B.z);
00169             
00170       dcosdB.x = rBinv*(cos_phi*B.x-A.x);
00171       dcosdB.y = rBinv*(cos_phi*B.y-A.y);
00172       dcosdB.z = rBinv*(cos_phi*B.z-A.z);
00173 
00174       K1 = K1/sin_phi;
00175 
00176       f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
00177       f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
00178       f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
00179                                               
00180       f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
00181       f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
00182       f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
00183                                               
00184       f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
00185                + r34.y*dcosdB.z - r34.z*dcosdB.y);
00186       f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
00187                + r34.z*dcosdB.x - r34.x*dcosdB.z);
00188       f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
00189                + r34.x*dcosdB.y - r34.y*dcosdB.x);
00190     }
00191     else
00192     {
00193       //  This angle is closer to 0 or 180 than it is to
00194       //  90, so use the cos version to avoid 1/sin terms
00195 
00196       //  Normalize C
00197       //      rC = 1.0/rC;
00198       C *= rCinv;
00199       
00200       Vector dsindC;
00201       Vector dsindB;
00202 
00203       dsindC.x = rCinv*(sin_phi*C.x-B.x);
00204       dsindC.y = rCinv*(sin_phi*C.y-B.y);
00205       dsindC.z = rCinv*(sin_phi*C.z-B.z);
00206 
00207       dsindB.x = rBinv*(sin_phi*B.x-C.x);
00208       dsindB.y = rBinv*(sin_phi*B.y-C.y);
00209       dsindB.z = rBinv*(sin_phi*B.z-C.z);
00210 
00211       K1 = -K1/cos_phi;
00212 
00213       f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00214                 - r23.x*r23.y*dsindC.y
00215                 - r23.x*r23.z*dsindC.z);
00216       f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00217                 - r23.y*r23.z*dsindC.z
00218                 - r23.y*r23.x*dsindC.x);
00219       f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00220                 - r23.z*r23.x*dsindC.x
00221                 - r23.z*r23.y*dsindC.y);
00222 
00223       f3 = cross(K1,dsindB,r23);
00224 
00225       f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00226              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00227              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00228              +dsindB.z*r34.y - dsindB.y*r34.z);
00229       f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00230              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00231              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00232              +dsindB.x*r34.z - dsindB.z*r34.x);
00233       f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00234              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00235              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00236              +dsindB.y*r34.x - dsindB.x*r34.y);
00237     }
00238 
00239   /* store the forces */
00240   //  p[0]->f[localIndex[0]] += f1;
00241   //  p[1]->f[localIndex[1]] += f2 - f1;
00242   //  p[2]->f[localIndex[2]] += f3 - f2;
00243   //  p[3]->f[localIndex[3]] += -f3;
00244 
00245   p[0]->f[localIndex[0]].x += f1.x;
00246   p[0]->f[localIndex[0]].y += f1.y;
00247   p[0]->f[localIndex[0]].z += f1.z;
00248 
00249   p[1]->f[localIndex[1]].x += f2.x - f1.x;
00250   p[1]->f[localIndex[1]].y += f2.y - f1.y;
00251   p[1]->f[localIndex[1]].z += f2.z - f1.z;
00252 
00253   p[2]->f[localIndex[2]].x += f3.x - f2.x;
00254   p[2]->f[localIndex[2]].y += f3.y - f2.y;
00255   p[2]->f[localIndex[2]].z += f3.z - f2.z;
00256 
00257   p[3]->f[localIndex[3]].x += -f3.x;
00258   p[3]->f[localIndex[3]].y += -f3.y;
00259   p[3]->f[localIndex[3]].z += -f3.z;  
00260 
00261     /* store the force for dihedral-only accelMD */
00262   if ( p[0]->af ) {
00263     p[0]->af[localIndex[0]].x += f1.x;
00264     p[0]->af[localIndex[0]].y += f1.y;
00265     p[0]->af[localIndex[0]].z += f1.z;
00266 
00267     p[1]->af[localIndex[1]].x += f2.x - f1.x;
00268     p[1]->af[localIndex[1]].y += f2.y - f1.y;
00269     p[1]->af[localIndex[1]].z += f2.z - f1.z;
00270 
00271     p[2]->af[localIndex[2]].x += f3.x - f2.x;
00272     p[2]->af[localIndex[2]].y += f3.y - f2.y;
00273     p[2]->af[localIndex[2]].z += f3.z - f2.z;
00274 
00275     p[3]->af[localIndex[3]].x += -f3.x;
00276     p[3]->af[localIndex[3]].y += -f3.y;
00277     p[3]->af[localIndex[3]].z += -f3.z;
00278   }
00279 
00280   DebugM(3, "::computeForce() -- ending with delta energy " << K << std::endl);
00281   reduction[dihedralEnergyIndex] += K;
00282   reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
00283   reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
00284   reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
00285   reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
00286   reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
00287   reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
00288   reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
00289   reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
00290   reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
00291 
00292   if (pressureProfileData) {
00293     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00294     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00295     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00296     BigReal z4 = p[3]->x[localIndex[3]].position.z;
00297     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00298     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00299     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00300     int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
00301     pp_clamp(n1, pressureProfileSlabs);
00302     pp_clamp(n2, pressureProfileSlabs);
00303     pp_clamp(n3, pressureProfileSlabs);
00304     pp_clamp(n4, pressureProfileSlabs);
00305     int p1 = p[0]->x[localIndex[0]].partition;
00306     int p2 = p[1]->x[localIndex[1]].partition;
00307     int p3 = p[2]->x[localIndex[2]].partition;
00308     int p4 = p[3]->x[localIndex[3]].partition;
00309     int pn = pressureProfileAtomTypes;
00310     pp_reduction(pressureProfileSlabs, n1, n2,
00311                 p1, p2, pn,
00312                 f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
00313                 pressureProfileData);
00314     pp_reduction(pressureProfileSlabs, n2, n3,
00315                 p2, p3, pn,
00316                 f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
00317                 pressureProfileData);
00318     pp_reduction(pressureProfileSlabs, n3, n4,
00319                 p3, p4, pn,
00320                 f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
00321                 pressureProfileData);
00322   }
00323 
00324  }
00325 }
00326 
00327 
00328 void DihedralElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00329 {
00330   reduction->item(REDUCTION_DIHEDRAL_ENERGY) += data[dihedralEnergyIndex];
00331   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[dihedralEnergyIndex_f];
00332   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[dihedralEnergyIndex_ti_1];
00333   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[dihedralEnergyIndex_ti_2];
00334   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00335   ADD_TENSOR(reduction,REDUCTION_VIRIAL_AMD_DIHE,data,virialIndex);
00336 }
00337 

Generated on Sat Sep 23 01:17:11 2017 for NAMD by  doxygen 1.4.7