1 /* -*- c++ -*- */ 2 #ifndef COULOMBFORCE_H 3 #define COULOMBFORCE_H 4 5 #include "GenericTopology.h" 6 #include "ScalarStructure.h" 7 #include "ExclusionTable.h" 8 #include "Parameter.h" 9 #include <string> 10 11 namespace ProtoMol { 12 13 //_________________________________________________________________ CoulombForce 14 class CoulombForce { 15 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 16 // This uses the weighted charges on each atom, so the Coulomb 17 // constant here is one. 18 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 19 20 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 21 // Constructors, destructors, assignment 22 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 23 public: 24 25 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 26 // New methods of class CoulombForce 27 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 28 public: operator()29 void operator()(Real &energy, 30 Real &force, 31 Real /*distSquared*/, 32 Real rDistSquared, 33 const Vector3D &, 34 const GenericTopology* topo, 35 int atom1, int atom2, 36 ExclusionClass excl) const { 37 energy = topo->atoms[atom1].scaledCharge * 38 topo->atoms[atom2].scaledCharge * 39 sqrt(rDistSquared); 40 if (topo->coulombScalingFactor != 1.0 && excl == EXCLUSION_MODIFIED) 41 energy *= topo->coulombScalingFactor; 42 43 force = energy * rDistSquared; 44 } 45 accumulateEnergy(ScalarStructure * energies,Real energy)46 static void accumulateEnergy(ScalarStructure* energies, Real energy) { 47 (*energies)[ScalarStructure::COULOMB] += energy; 48 } 49 getEnergy(const ScalarStructure * energies)50 static Real getEnergy(const ScalarStructure* energies) {return (*energies)[ScalarStructure::COULOMB];} 51 52 // Parsing getId()53 static std::string getId() {return keyword;} getParameterSize()54 static unsigned int getParameterSize() {return 0;} getParameters(std::vector<Parameter> &)55 void getParameters(std::vector<Parameter>&) const{} make(std::string &,const std::vector<Value> &)56 static CoulombForce make(std::string&, const std::vector<Value>&) { 57 return CoulombForce(); 58 } 59 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 60 // Sub Classes 61 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 62 public: 63 class C1 { 64 public: kernel(Real r)65 static Real kernel(Real r) {return (1.0/r);} dKernel(Real r)66 static Real dKernel(Real r) {r=1.0/r;return (-r*r);} kernelR(Real rr)67 static Real kernelR(Real rr) {return (rr);} dKernelR(Real rr)68 static Real dKernelR(Real rr) {return (-rr*rr);} smooth(Real r,Real,Real cr)69 static Real smooth(Real r,Real /*c*/,Real cr) {return (cr*(1.5-0.5*r*r*cr*cr));} smooth0(Real,Real cr)70 static Real smooth0(Real /*c*/,Real cr) {return (1.5*cr);} dSmooth(Real r,Real,Real cr)71 static Real dSmooth(Real r,Real /*c*/,Real cr) {return (-r*cr*cr*cr);} smoothKernel(Real r,Real c,Real cr)72 static Real smoothKernel(Real r,Real c,Real cr) {return (r<c ? smooth(r,c,cr) : kernel(r));} dSmoothKernel(Real r,Real c,Real cr)73 static Real dSmoothKernel(Real r,Real c,Real cr){return (r<c ? dSmooth(r,c,cr) : dKernel(r));} getKeyword()74 static std::string getKeyword() {return keyword;} getForceKeyword()75 static std::string getForceKeyword() {return CoulombForce::keyword;} 76 public: 77 static const std::string keyword; 78 }; 79 public: 80 class C2 { 81 public: kernel(Real r)82 static Real kernel(Real r) {return (1.0/r);} dKernel(Real r)83 static Real dKernel(Real r) {r=1.0/r;return (-r*r);} kernelR(Real rr)84 static Real kernelR(Real rr) {return (rr);} dKernelR(Real rr)85 static Real dKernelR(Real rr) {return (-rr*rr);} smooth(Real r,Real,Real cr)86 static Real smooth(Real r,Real /*c*/,Real cr) {r=r*r*cr*cr;return (cr*(1.875-r*(1.25-0.375*r)));} smooth0(Real,Real cr)87 static Real smooth0(Real /*c*/,Real cr) {return (1.875*cr);} dSmooth(Real r,Real c,Real cr)88 static Real dSmooth(Real r,Real c,Real cr) {c=r*cr*cr;return (cr*c*(1.5*r*c-2.5));} smoothKernel(Real r,Real c,Real cr)89 static Real smoothKernel(Real r,Real c,Real cr) {return (r<c ? smooth(r,c,cr) : kernel(r));} dSmoothKernel(Real r,Real c,Real cr)90 static Real dSmoothKernel(Real r,Real c,Real cr){return (r<c ? dSmooth(r,c,cr) : dKernel(r));} getKeyword()91 static std::string getKeyword() {return keyword;} getForceKeyword()92 static std::string getForceKeyword() {return CoulombForce::keyword;} 93 public: 94 static const std::string keyword; 95 }; 96 public: 97 class C3 { 98 public: kernel(Real r)99 static Real kernel(Real r) {return (1.0/r);} dKernel(Real r)100 static Real dKernel(Real r) {r=1.0/r;return (-r*r);} kernelR(Real rr)101 static Real kernelR(Real rr) {return (rr);} dKernelR(Real rr)102 static Real dKernelR(Real rr) {return (-rr*rr);} smooth(Real r,Real,Real cr)103 static Real smooth(Real r,Real /*c*/,Real cr) {r=r*r*cr*cr;return (0.0625*cr*(35.0-r*(35.0-r*(21.0-5.0*r))));} smooth0(Real,Real cr)104 static Real smooth0(Real /*c*/,Real cr) {return (2.1875*cr);} dSmooth(Real r,Real c,Real cr)105 static Real dSmooth(Real r,Real c,Real cr) {c=r*r*cr*cr;return (r*cr*cr*cr*(-4.375 + c * (5.25 - 1.875 * c)));} smoothKernel(Real r,Real c,Real cr)106 static Real smoothKernel(Real r,Real c,Real cr) {return (r<c ? smooth(r,c,cr) : kernel(r));} dSmoothKernel(Real r,Real c,Real cr)107 static Real dSmoothKernel(Real r,Real c,Real cr){return (r<c ? dSmooth(r,c,cr) : dKernel(r));} getKeyword()108 static std::string getKeyword() {return keyword;} getForceKeyword()109 static std::string getForceKeyword() {return CoulombForce::keyword;} 110 public: 111 static const std::string keyword; 112 }; 113 public: 114 class C4 { 115 public: kernel(Real r)116 static Real kernel(Real r) {return (1.0/r);} dKernel(Real r)117 static Real dKernel(Real r) {r=1.0/r;return (-r*r);} kernelR(Real rr)118 static Real kernelR(Real rr) {return (rr);} dKernelR(Real rr)119 static Real dKernelR(Real rr) {return (-rr*rr);} smooth(Real r,Real,Real cr)120 static Real smooth(Real r,Real /*c*/,Real cr) {r=r*r*cr*cr;return (0.0078125*cr*(315.0-r*(420.0-r*(378.0-r*(180.0-r*35.0)))));} smooth0(Real,Real cr)121 static Real smooth0(Real /*c*/,Real cr) {return (2.4609375*cr);} dSmooth(Real r,Real c,Real cr)122 static Real dSmooth(Real r,Real c,Real cr) {c=r*r*cr*cr;return(-r*cr*cr*cr*(6.5625-c*(11.8125-c*(8.4375-c*2.1875)))); 123 } smoothKernel(Real r,Real c,Real cr)124 static Real smoothKernel(Real r,Real c,Real cr) {return (r<c ? smooth(r,c,cr) : kernel(r));} dSmoothKernel(Real r,Real c,Real cr)125 static Real dSmoothKernel(Real r,Real c,Real cr){return (r<c ? dSmooth(r,c,cr) : dKernel(r));} getKeyword()126 static std::string getKeyword() {return keyword;} getForceKeyword()127 static std::string getForceKeyword() {return CoulombForce::keyword;} 128 public: 129 static const std::string keyword; 130 }; 131 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 132 // My data members 133 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 134 public: 135 static const std::string keyword; 136 private: 137 138 }; 139 140 //______________________________________________________________________ INLINES 141 } 142 #endif /* COULOMBFORCE_H */ 143