ComputeImpropers.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeImpropers.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 ImproperElem::pressureProfileSlabs = 0;
00020 int ImproperElem::pressureProfileAtomTypes = 1;
00021 BigReal ImproperElem::pressureProfileThickness = 0;
00022 BigReal ImproperElem::pressureProfileMin = 0;
00023 
00024 void ImproperElem::getMoleculePointers
00025     (Molecule* mol, int* count, int32*** byatom, Improper** structarray)
00026 {
00027 #ifdef MEM_OPT_VERSION
00028   NAMD_die("Should not be called in ImproperElem::getMoleculePointers in memory optimized version!");
00029 #else
00030   *count = mol->numImpropers;
00031   *byatom = mol->impropersByAtom;
00032   *structarray = mol->impropers;
00033 #endif
00034 }
00035 
00036 void ImproperElem::getParameterPointers(Parameters *p, const ImproperValue **v) {
00037   *v = p->improper_array;
00038 }
00039 
00040 void ImproperElem::computeForce(ImproperElem *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 ImproperElem &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 ImproperValue * const(&value)(tup.value);
00065 
00066   DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00067                << localIndex[1] << " " << localIndex[2] << " " <<
00068                localIndex[3] << std::endl);
00069 
00070   // Vector r12, r23, r34;      // vector between atoms
00071   Vector A,B,C;         // cross products
00072   BigReal rA, rB, rC;   // length of vectors A, B, and C
00073   BigReal energy=0;     // energy from the angle
00074   BigReal phi;          // angle between the plans
00075   double cos_phi;       // cos(phi)
00076   double sin_phi;       // sin(phi)
00077   Vector dcosdA;        // Derivative d(cos(phi))/dA
00078   Vector dcosdB;        // Derivative d(cos(phi))/dB
00079   Vector dsindC;        // Derivative d(sin(phi))/dC
00080   Vector dsindB;        // Derivative d(sin(phi))/dB
00081   BigReal K,K1;         // energy constants
00082   BigReal diff;         // for periodicity
00083   Force f1(0,0,0),f2(0,0,0),f3(0,0,0);  // force components
00084 
00085   //DebugM(3, "::computeForce() -- starting with improper type " << improperType << std::endl);
00086 
00087   // get the improper information
00088   int multiplicity = value->multiplicity;
00089 
00090   //  Calculate the vectors between atoms
00091   const Position & pos0 = p[0]->x[localIndex[0]].position;
00092   const Position & pos1 = p[1]->x[localIndex[1]].position;
00093   const Position & pos2 = p[2]->x[localIndex[2]].position;
00094   const Position & pos3 = p[3]->x[localIndex[3]].position;
00095   const Vector r12 = lattice.delta(pos0,pos1);
00096   const Vector r23 = lattice.delta(pos1,pos2);
00097   const Vector r34 = lattice.delta(pos2,pos3);
00098 
00099   //  Calculate the cross products
00100   A = cross(r12,r23);
00101   B = cross(r23,r34);
00102   C = cross(r23,A);
00103 
00104   //  Calculate the distances
00105   rA = A.length();
00106   rB = B.length();
00107   rC = C.length();
00108 
00109   //  Calculate the sin and cos
00110   cos_phi = A*B/(rA*rB);
00111   sin_phi = C*B/(rC*rB);
00112 
00113   //  Normalize B
00114   rB = 1.0/rB;
00115   B *= rB;
00116 
00117   phi= -atan2(sin_phi,cos_phi);
00118 
00119   if (fabs(sin_phi) > 0.1)
00120   {
00121     //  Normalize A
00122     rA = 1.0/rA;
00123     A *= rA;
00124     dcosdA = rA*(cos_phi*A-B);
00125     dcosdB = rB*(cos_phi*B-A);
00126   }
00127   else
00128   {
00129     //  Normalize C
00130     rC = 1.0/rC;
00131     C *= rC;
00132     dsindC = rC*(sin_phi*C-B);
00133     dsindB = rB*(sin_phi*B-C);
00134   }
00135 
00136   //  Loop through the multiple parameter sets for this
00137   //  bond.  We will only loop more than once if this
00138   //  has multiple parameter sets from Charmm22
00139   for (int mult_num=0; mult_num<multiplicity; mult_num++)
00140   {
00141     /* get angle information */
00142     Real k = value->values[mult_num].k * scale;
00143     Real delta = value->values[mult_num].delta;
00144     int n = value->values[mult_num].n;
00145 
00146     //  Calculate the energy
00147     if (n)
00148     {
00149       //  Periodicity is greater than 0, so use cos form
00150       K = k*(1+cos(n*phi - delta));
00151       K1 = -n*k*sin(n*phi - delta);
00152     }
00153     else
00154     {
00155       //  Periodicity is 0, so just use the harmonic form
00156       diff = phi-delta;
00157       if (diff < -PI)           diff += TWOPI;
00158       else if (diff > PI)       diff -= TWOPI;
00159 
00160       K = k*diff*diff;
00161       K1 = 2.0*k*diff;
00162     }
00163 
00164     //  Add the energy from this improper to the total energy
00165     energy += K;
00166 
00167     //  Next, we want to calculate the forces.  In order
00168     //  to do that, we first need to figure out whether the
00169     //  sin or cos form will be more stable.  For this,
00170     //  just look at the value of phi
00171     if (fabs(sin_phi) > 0.1)
00172     {
00173       //  use the sin version to avoid 1/cos terms
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       K1 = -K1/cos_phi;
00196 
00197       f1.x += K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00198                 - r23.x*r23.y*dsindC.y
00199                 - r23.x*r23.z*dsindC.z);
00200       f1.y += K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00201                 - r23.y*r23.z*dsindC.z
00202                 - r23.y*r23.x*dsindC.x);
00203       f1.z += K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00204                 - r23.z*r23.x*dsindC.x
00205                 - r23.z*r23.y*dsindC.y);
00206 
00207       f3 += cross(K1,dsindB,r23);
00208 
00209       f2.x += K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00210              +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00211              +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00212              +dsindB.z*r34.y - dsindB.y*r34.z);
00213       f2.y += K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00214              +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00215              +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00216              +dsindB.x*r34.z - dsindB.z*r34.x);
00217       f2.z += K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00218              +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00219              +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00220              +dsindB.y*r34.x - dsindB.x*r34.y);
00221     }
00222   } /* for multiplicity */
00223 
00224   //fepb - BKR scaling of alchemical bonded terms
00225   //       NB: TI derivative is the _unscaled_ energy.
00226   if ( simParams->alchOn ) {
00227     switch ( mol->get_fep_bonded_type(atomID, 4) ) {
00228     case 1:
00229       reduction[improperEnergyIndex_ti_1] += energy;
00230       reduction[improperEnergyIndex_f] += (bond_lambda_12 - bond_lambda_1) * 
00231                                            energy;
00232       energy *= bond_lambda_1;
00233       f1 *= bond_lambda_1;
00234       f2 *= bond_lambda_1;
00235       f3 *= bond_lambda_1;
00236       break;
00237     case 2:
00238       reduction[improperEnergyIndex_ti_2] += energy;
00239       reduction[improperEnergyIndex_f] += (bond_lambda_22 - bond_lambda_2) *
00240                                            energy;
00241       energy *= bond_lambda_2;
00242       f1 *= bond_lambda_2;
00243       f2 *= bond_lambda_2;
00244       f3 *= bond_lambda_2;
00245       break;
00246     }
00247   }
00248   //fepe
00249 
00250   /* store the forces */
00251   p[0]->f[localIndex[0]] += f1;
00252   p[1]->f[localIndex[1]] += f2 - f1;
00253   p[2]->f[localIndex[2]] += f3 - f2;
00254   p[3]->f[localIndex[3]] += -f3;
00255 
00256   DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00257   reduction[improperEnergyIndex] += energy;
00258   reduction[virialIndex_XX] += ( f1.x * r12.x + f2.x * r23.x + f3.x * r34.x );
00259   reduction[virialIndex_XY] += ( f1.x * r12.y + f2.x * r23.y + f3.x * r34.y );
00260   reduction[virialIndex_XZ] += ( f1.x * r12.z + f2.x * r23.z + f3.x * r34.z );
00261   reduction[virialIndex_YX] += ( f1.y * r12.x + f2.y * r23.x + f3.y * r34.x );
00262   reduction[virialIndex_YY] += ( f1.y * r12.y + f2.y * r23.y + f3.y * r34.y );
00263   reduction[virialIndex_YZ] += ( f1.y * r12.z + f2.y * r23.z + f3.y * r34.z );
00264   reduction[virialIndex_ZX] += ( f1.z * r12.x + f2.z * r23.x + f3.z * r34.x );
00265   reduction[virialIndex_ZY] += ( f1.z * r12.y + f2.z * r23.y + f3.z * r34.y );
00266   reduction[virialIndex_ZZ] += ( f1.z * r12.z + f2.z * r23.z + f3.z * r34.z );
00267 
00268   if (pressureProfileData) {
00269     BigReal z1 = p[0]->x[localIndex[0]].position.z;
00270     BigReal z2 = p[1]->x[localIndex[1]].position.z;
00271     BigReal z3 = p[2]->x[localIndex[2]].position.z;
00272     BigReal z4 = p[3]->x[localIndex[3]].position.z;
00273     int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00274     int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00275     int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00276     int n4 = (int)floor((z4-pressureProfileMin)/pressureProfileThickness);
00277     pp_clamp(n1, pressureProfileSlabs);
00278     pp_clamp(n2, pressureProfileSlabs);
00279     pp_clamp(n3, pressureProfileSlabs);
00280     pp_clamp(n4, pressureProfileSlabs);
00281     int p1 = p[0]->x[localIndex[0]].partition;
00282     int p2 = p[1]->x[localIndex[1]].partition;
00283     int p3 = p[2]->x[localIndex[2]].partition;
00284     int p4 = p[3]->x[localIndex[3]].partition;
00285     int pn = pressureProfileAtomTypes;
00286     pp_reduction(pressureProfileSlabs, n1, n2,
00287                 p1, p2, pn,
00288                 f1.x * r12.x, f1.y * r12.y, f1.z * r12.z,
00289                 pressureProfileData);
00290     pp_reduction(pressureProfileSlabs, n2, n3,
00291                 p2, p3, pn,
00292                 f2.x * r23.x, f2.y * r23.y, f2.z * r23.z,
00293                 pressureProfileData);
00294     pp_reduction(pressureProfileSlabs, n3, n4,
00295                 p3, p4, pn,
00296                 f3.x * r34.x, f3.y * r34.y, f3.z * r34.z,
00297                 pressureProfileData);
00298   }
00299 
00300  }
00301 }
00302 
00303 void ImproperElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00304 {
00305   reduction->item(REDUCTION_IMPROPER_ENERGY) += data[improperEnergyIndex];
00306   reduction->item(REDUCTION_BONDED_ENERGY_F) += data[improperEnergyIndex_f];
00307   reduction->item(REDUCTION_BONDED_ENERGY_TI_1) += data[improperEnergyIndex_ti_1];
00308   reduction->item(REDUCTION_BONDED_ENERGY_TI_2) += data[improperEnergyIndex_ti_2];
00309   ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00310 }
00311 

Generated on Sat Nov 18 01:17:12 2017 for NAMD by  doxygen 1.4.7