ComputeGromacsPair.C

Go to the documentation of this file.
00001 
00007 #include "InfoStream.h"
00008 #include "ComputeGromacsPair.h"
00009 #include "Molecule.h"
00010 #include "Parameters.h"
00011 #include "LJTable.h"
00012 #include "Node.h"
00013 #include "ReductionMgr.h"
00014 #include "Lattice.h"
00015 #include "PressureProfile.h"
00016 #include "Debug.h"
00017 #include <iostream>
00018 
00019 // static initialization
00020 int GromacsPairElem::pressureProfileSlabs = 0;
00021 int GromacsPairElem::pressureProfileAtomTypes = 1;
00022 BigReal GromacsPairElem::pressureProfileThickness = 0;
00023 BigReal GromacsPairElem::pressureProfileMin = 0;
00024 
00025 
00026 void GromacsPairElem::getMoleculePointers
00027     (Molecule* mol, int* count, int32*** byatom, GromacsPair** structarray)
00028 {
00029 #ifdef MEM_OPT_VERSION
00030     NAMD_die("Should not be called in GromacsPairElem::getMoleculePointers in memory optimized version!");
00031 #else
00032     *count = mol->numLJPair;
00033     *byatom = mol->gromacsPairByAtom;
00034     *structarray = mol->gromacsPair;
00035 #endif
00036 }
00037 
00038 //void GromacsPairElem::getParameterPointers(Parameters *p, const GromacsPairValue **v) {
00039 void GromacsPairElem::getParameterPointers(Parameters *p, const GromacsPairValue **v) {
00040     // JLai
00041     *v=p->gromacsPair_array;
00042 }
00043 
00044 void GromacsPairElem::computeForce(GromacsPairElem *tuples, int ntuple, BigReal *reduction, 
00045                             BigReal *pressureProfileData)
00046 {
00047  const Lattice & lattice = tuples[0].p[0]->p->lattice;
00048 
00049  for ( int ituple=0; ituple<ntuple; ++ituple ) {
00050   const GromacsPairElem &tup = tuples[ituple];
00051   enum { size = 2 };
00052   const AtomID (&atomID)[size](tup.atomID);
00053   const int    (&localIndex)[size](tup.localIndex);
00054   TuplePatchElem * const(&p)[size](tup.p);
00055   const Real (&scale)(tup.scale);
00056   const GromacsPairValue * const(&value)(tup.value);
00057 
00058     DebugM(1, "::computeforce() localIndex = " << localIndex [0] << " "
00059            << localIndex[1] << std::endl);
00060 
00061     // compute vectors between atoms and their distances
00062     const Vector r12 = lattice.delta(p[0]->x[localIndex[0]].position,
00063                                      p[1]->x[localIndex[1]].position);
00064 
00065     if ( p[0]->patchID == p[1]->patchID && localIndex[0] == localIndex[1] ) {
00066         continue;
00067     }
00068     SimParameters *simParams = Node::Object()->simParameters;
00069     BigReal cutoff = simParams->cutoff;
00070     BigReal r12_len = r12.length2();
00071     if ( r12_len == 0 ) continue;
00072     //if ( r12_len > cutoff) {
00073         //continue;
00074     //}
00075     BigReal ri2 = 1.0/r12_len;
00076     BigReal ri = sqrt(ri2);
00077     BigReal ri6 = ri2*ri2*ri2;
00078     BigReal ri7 = ri*ri6;
00079     BigReal ri8 = ri2*ri2*ri2*ri2;
00080     BigReal ri12 = ri6*ri6;
00081     BigReal ri13 = ri12*ri;
00082     BigReal ri14 = ri12*ri2;
00083 
00084     BigReal energy = 0;
00085     BigReal diff = 0;
00086     BigReal pairC12 = value->pairC12;
00087     BigReal pairC6 = value->pairC6;
00088 
00089     // Add the energy for the 12-6 LJ interaction to the total energy
00090     energy = (pairC12*ri12) - (pairC6*ri6);
00091     // This is a dirty hack; currently the code is looping over N^2 instead of N(N-1)/2 
00092     // This is happening because the LJ list is 2x long as it has to be
00093     //energy *= 0.5; 
00094 
00095     // Determine the magnitude of the force
00096     diff = ((12*pairC12*ri14) - (6*pairC6*ri8));
00097     // This is a dirty hack; currently the code is looping over N^2 instead of N(N-1)/2 
00098     // This is happening because the LJ list is 2x long as it has to be
00099     //diff *= 0.5; 
00100     //std::cout << "Force: " << diff << " " << pairC12 << " " << pairC6 << " " << r12.length() << "\n";
00101 
00102     //Scale the force vector accordingly
00103     const Force f12 = diff * r12;
00104     //std::cout << "Atoms: " << localIndex[0] << " " << localIndex[1] << " " << f12.length() << " " << f12 << "\n";
00105     //std::cout << "Force2: " << f12 << "\n";
00106 
00107     // Now add the forces to each force vector
00108     p[0]->f[localIndex[0]] += f12;
00109     p[1]->f[localIndex[1]] -= f12;
00110 
00111     DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00112     reduction[gromacsPairEnergyIndex] += energy;
00113     reduction[virialIndex_XX] += f12.x * r12.x;
00114     reduction[virialIndex_XY] += f12.x * r12.y;
00115     reduction[virialIndex_XZ] += f12.x * r12.z;
00116     reduction[virialIndex_YX] += f12.y * r12.x;
00117     reduction[virialIndex_YY] += f12.y * r12.y;
00118     reduction[virialIndex_YZ] += f12.y * r12.z;
00119     reduction[virialIndex_ZX] += f12.z * r12.x;
00120     reduction[virialIndex_ZY] += f12.z * r12.y;
00121     reduction[virialIndex_ZZ] += f12.z * r12.z;
00122 
00123     if (pressureProfileData) {
00124         BigReal z1 = p[0]->x[localIndex[0]].position.z;
00125         BigReal z2 = p[1]->x[localIndex[1]].position.z;
00126         int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00127         int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00128         pp_clamp(n1, pressureProfileSlabs);
00129         pp_clamp(n2, pressureProfileSlabs);
00130         int p1 = p[0]->x[localIndex[0]].partition;
00131         int p2 = p[1]->x[localIndex[1]].partition;
00132         int pn = pressureProfileAtomTypes;
00133         pp_reduction(pressureProfileSlabs,
00134                      n1, n2, 
00135                      p1, p2, pn, 
00136                      f12.x * r12.x, f12.y * r12.y, f12.z * r12.z,
00137                      pressureProfileData);
00138     } 
00139 
00140  }
00141 }
00142 
00143 void GromacsPairElem::submitReductionData(BigReal *data, SubmitReduction *reduction)
00144 {
00145     // JLai
00146     reduction->item(REDUCTION_GRO_LJ_ENERGY) += data[gromacsPairEnergyIndex];
00147     ADD_TENSOR(reduction,REDUCTION_VIRIAL_NBOND,data,virialIndex);
00148     return;
00149 }

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