1 /* -*- c++ -*- */ 2 #ifndef ANGLESYSTEMFORCE_H 3 #define ANGLESYSTEMFORCE_H 4 5 #include "SystemForce.h" 6 #include "AngleSystemForceBase.h" 7 #include "ScalarStructure.h" 8 #include "Parallel.h" 9 10 namespace ProtoMol { 11 //_________________________________________________________________ AngleSystemForce 12 template<class TBoundaryConditions> 13 class AngleSystemForce : public SystemForce, private AngleSystemForceBase { 14 15 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 16 // Constructors, destructors, assignment 17 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 18 public: 19 20 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 21 // New methods of class AngleSystemForce 22 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 23 void calcAngle(const TBoundaryConditions &boundary, 24 const Angle& currentAngle, 25 const Vector3DBlock* positions, 26 Vector3DBlock* forces, 27 ScalarStructure* energies); 28 29 Real calcAngleEnergy(const TBoundaryConditions &boundary, 30 const Angle& currentAngle, 31 const Vector3DBlock* positions); 32 33 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 34 // From class SystemForce 35 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 36 public: 37 virtual void evaluate(const GenericTopology* topo, 38 const Vector3DBlock* positions, 39 Vector3DBlock* forces, 40 ScalarStructure* energies); 41 42 virtual void parallelEvaluate(const GenericTopology* topo, 43 const Vector3DBlock* positions, 44 Vector3DBlock* forces, 45 ScalarStructure* energies); 46 47 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 48 // From class Force 49 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 50 public: getKeyword()51 virtual std::string getKeyword() const{return keyword;} 52 virtual unsigned int numberOfBlocks(const GenericTopology* topo, 53 const Vector3DBlock* pos); 54 private: doMake(std::string &,std::vector<Value>)55 virtual Force* doMake(std::string&, std::vector<Value>) const { 56 return (new AngleSystemForce()); 57 } 58 59 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 60 // From class Makeable 61 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 62 public: getIdNoAlias()63 virtual std::string getIdNoAlias() const{return keyword;} getParameterSize()64 virtual unsigned int getParameterSize() const{return 0;} getParameters(std::vector<Parameter> &)65 virtual void getParameters(std::vector<Parameter>&) const {} 66 private: doSetParameters(std::string &,std::vector<Value>)67 virtual void doSetParameters(std::string&, std::vector<Value>){} 68 69 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 70 // My data members 71 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 72 private: 73 74 }; 75 76 //______________________________________________________________________ INLINES 77 78 template<class TBoundaryConditions> evaluate(const GenericTopology * topo,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)79 inline void AngleSystemForce<TBoundaryConditions>::evaluate(const GenericTopology* topo, 80 const Vector3DBlock* positions, 81 Vector3DBlock* forces, 82 ScalarStructure* energies) { 83 84 const TBoundaryConditions &boundary = 85 (dynamic_cast<const SemiGenericTopology<TBoundaryConditions>& >(*topo)).boundaryConditions; 86 for (unsigned int i = 0; i < topo->angles.size(); i++) 87 calcAngle(boundary, topo->angles[i], positions, forces, energies); 88 } 89 90 template<class TBoundaryConditions> calcAngle(const TBoundaryConditions & boundary,const Angle & currentAngle,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)91 inline void AngleSystemForce<TBoundaryConditions>::calcAngle(const TBoundaryConditions &boundary, 92 const Angle& currentAngle, 93 const Vector3DBlock* positions, 94 Vector3DBlock* forces, 95 ScalarStructure* energies){ 96 97 int a1 = currentAngle.atom1; 98 int a2 = currentAngle.atom2; 99 int a3 = currentAngle.atom3; 100 101 Real restAngle = currentAngle.restAngle; 102 Real forceConstant = currentAngle.forceConstant; 103 Real ureyBradleyConstant = currentAngle.ureyBradleyConstant; 104 Real ureyBradleyRestLength = currentAngle.ureyBradleyRestLength; 105 // Report::report << restAngle <<","<<forceConstant<<","<<ureyBradleyConstant<<","<<ureyBradleyRestLength<<Report::endr; 106 107 Vector3D atom1((*positions)[a1]); 108 Vector3D atom2((*positions)[a2]); 109 Vector3D atom3((*positions)[a3]); 110 111 //Vector3D r12 = atom1 - atom2; // Vector from atom 1 to atom 2. 112 Vector3D r12(boundary.minimalDifference(atom2,atom1)); 113 //Vector3D r32 = atom3 - atom2; // Vector from atom 3 to atom 2. 114 Vector3D r32(boundary.minimalDifference(atom2,atom3)); 115 //Vector3D r13 = atom1 - atom3; // Vector from atom 1 to atom 3. 116 Vector3D r13(boundary.minimalDifference(atom3,atom1)); 117 Real d12 = r12.norm(); // Distance between atom 1 and 2. 118 Real d32 = r32.norm(); // Distance between atom 3 and 2. 119 Real d13 = r13.norm(); // Distance between atom 1 and 3. 120 121 122 // Calculate theta. 123 Real theta = atan2((r12.cross(r32)).norm(),r12.dot(r32)); 124 Real sinTheta = sin(theta); 125 Real cosTheta = cos(theta); 126 127 // Calculate dpot/dtheta 128 Real dpotdtheta = 2.0 * forceConstant * ( theta - restAngle); 129 130 // Calculate dr/dx, dr/dy, dr/dz. 131 Vector3D dr12(r12/d12); 132 Vector3D dr32(r32/d32); 133 Vector3D dr13(r13/d13); 134 // Calulate dtheta/dx, dtheta/dy, dtheta/dz. 135 Vector3D dtheta1((dr12 * cosTheta - dr32)/(sinTheta*d12)); // atom1 136 Vector3D dtheta3((dr32 * cosTheta - dr12)/(sinTheta*d32)); // atom3 137 138 // Calculate Urey Bradley force. 139 Vector3D ureyBradleyforce1(dr13 * (2.0 * ureyBradleyConstant * (d13 - ureyBradleyRestLength))); 140 Vector3D ureyBradleyforce3(-ureyBradleyforce1); 141 142 // Calculate force on atom1 due to atom 2 and 3. 143 Vector3D force1(-dtheta1 * dpotdtheta - ureyBradleyforce1); 144 145 // Calculate force on atom3 due to atom 1 and 2. 146 Vector3D force3(-dtheta3 * dpotdtheta - ureyBradleyforce3); 147 148 // Calculate force on atom2 due to atom 1 and 3. 149 Vector3D force2(-force1 - force3); 150 151 // Add to the total force. 152 (*forces)[a1] += force1; 153 (*forces)[a2] += force2; 154 (*forces)[a3] += force3; 155 156 // Calculate Energy. 157 Real eHarmonic = forceConstant * (theta - restAngle)*(theta - restAngle); 158 Real eUreyBradley = ureyBradleyConstant * (d13 - ureyBradleyRestLength)*(d13 - ureyBradleyRestLength); 159 // Add Energy 160 (*energies)[ScalarStructure::ANGLE] += eHarmonic + eUreyBradley; 161 162 // Add virial 163 if(energies->virial()){ 164 Real xy = force1.x * r12.y + force3.x * r32.y; 165 Real xz = force1.x * r12.z + force3.x * r32.z; 166 Real yz = force1.y * r12.z + force3.y * r32.z; 167 (*energies)[ScalarStructure::VIRIALXX] += force1.x * r12.x + force3.x * r32.x; 168 (*energies)[ScalarStructure::VIRIALXY] += xy; 169 (*energies)[ScalarStructure::VIRIALXZ] += xz; 170 (*energies)[ScalarStructure::VIRIALYX] += xy; 171 (*energies)[ScalarStructure::VIRIALYY] += force1.y * r12.y + force3.y * r32.y; 172 (*energies)[ScalarStructure::VIRIALYZ] += yz; 173 (*energies)[ScalarStructure::VIRIALZX] += xz; 174 (*energies)[ScalarStructure::VIRIALZY] += yz; 175 (*energies)[ScalarStructure::VIRIALZZ] += force1.z * r12.z + force3.z * r32.z; 176 } 177 } 178 179 template<class TBoundaryConditions> calcAngleEnergy(const TBoundaryConditions & boundary,const Angle & currentAngle,const Vector3DBlock * positions)180 inline Real AngleSystemForce<TBoundaryConditions>::calcAngleEnergy(const TBoundaryConditions &boundary, 181 const Angle& currentAngle, 182 const Vector3DBlock* positions) 183 { 184 int a1, a2, a3; 185 Real restAngle, forceConstant, ureyBradleyConstant, ureyBradleyRestLength; 186 Vector3D atom1, atom2, atom3, r12, r32, r13; 187 Real d13, theta, eHarmonic, eUreyBradley; 188 189 a1 = currentAngle.atom1; 190 a2 = currentAngle.atom2; 191 a3 = currentAngle.atom3; 192 restAngle = currentAngle.restAngle; 193 forceConstant = currentAngle.forceConstant; 194 ureyBradleyConstant = currentAngle.ureyBradleyConstant; 195 ureyBradleyRestLength = currentAngle.ureyBradleyRestLength; 196 197 atom1 = (*positions)[a1]; 198 atom2 = (*positions)[a2]; 199 atom3 = (*positions)[a3]; 200 201 //r12 = atom1 - atom2; // Vector from atom 1 to atom 2. 202 r12 = boundary.minimalDifference(atom2,atom1); 203 //r32 = atom3 - atom2; // Vector from atom 3 to atom 2. 204 r32 = boundary.minimalDifference(atom2,atom3); 205 //r13 = atom1 - atom3; // Vector from atom 1 to atom 3. 206 r13 = boundary.minimalDifference(atom3,atom1); 207 d13 = r13.norm(); // Distance between atom 1 and 3. 208 209 // Calculate theta. 210 theta = atan2((r12.cross(r32)).norm(),r12.dot(r32)); 211 212 // Calculate Energy. 213 eHarmonic = forceConstant * (theta - restAngle)*(theta - restAngle); 214 eUreyBradley = ureyBradleyConstant * (d13 - ureyBradleyRestLength) 215 *(d13 - ureyBradleyRestLength); 216 217 return (eHarmonic + eUreyBradley); 218 } 219 template<class TBoundaryConditions> parallelEvaluate(const GenericTopology * topo,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)220 inline void AngleSystemForce<TBoundaryConditions>::parallelEvaluate(const GenericTopology* topo, 221 const Vector3DBlock* positions, 222 Vector3DBlock* forces, 223 ScalarStructure* energies) 224 { 225 const TBoundaryConditions &boundary = 226 (dynamic_cast<const SemiGenericTopology<TBoundaryConditions>& >(*topo)).boundaryConditions; 227 228 unsigned int n = topo->angles.size(); 229 unsigned int count = numberOfBlocks(topo,positions); 230 231 for(unsigned int i = 0;i<count;i++){ 232 if(Parallel::next()){ 233 int to = (n*(i+1))/count; 234 if(to > static_cast<int>(n)) 235 to = n; 236 int from = (n*i)/count; 237 for (int j = from; j < to; j++) 238 calcAngle(boundary, topo->angles[j], positions, forces, energies); 239 } 240 } 241 } 242 243 template<class TBoundaryConditions> numberOfBlocks(const GenericTopology * topo,const Vector3DBlock *)244 inline unsigned int AngleSystemForce<TBoundaryConditions>::numberOfBlocks(const GenericTopology* topo, 245 const Vector3DBlock*){ 246 return std::min(Parallel::getAvailableNum(),static_cast<int>(topo->angles.size())); 247 } 248 249 } 250 #endif /* ANGLESYSTEMFORCE_H */ 251