1 /* -*- c++ -*- */ 2 #ifndef ISGANGLESYSTEMFORCE_H 3 #define ISGANGLESYSTEMFORCE_H 4 5 #include "SystemForce.h" 6 #include "iSGAngleSystemForceBase.h" 7 #include "ScalarStructure.h" 8 #include "Parallel.h" 9 10 namespace ProtoMol { 11 //_________________________________________________________________ iSGAngleSystemForce 12 template<class TBoundaryConditions> 13 class iSGAngleSystemForce : public SystemForce, private iSGAngleSystemForceBase { 14 15 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 16 // Constructors, destructors, assignment 17 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 18 public: 19 20 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 21 // New methods of class iSGAngleSystemForce 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 iSGAngleSystemForce()); 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 67 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 68 // My data members 69 //~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 70 private: 71 72 }; 73 74 //______________________________________________________________________ INLINES 75 76 template<class TBoundaryConditions> evaluate(const GenericTopology * topo,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)77 inline void iSGAngleSystemForce<TBoundaryConditions>::evaluate(const GenericTopology* topo, 78 const Vector3DBlock* positions, 79 Vector3DBlock* forces, 80 ScalarStructure* energies) { 81 82 const TBoundaryConditions &boundary = 83 (dynamic_cast<const SemiGenericTopology<TBoundaryConditions>& >(*topo)).boundaryConditions; 84 for (unsigned int i = 0; i < topo->angles.size(); i++) 85 calcAngle(boundary, topo->angles[i], positions, forces, energies); 86 } 87 88 template<class TBoundaryConditions> calcAngle(const TBoundaryConditions & boundary,const Angle & currentAngle,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)89 inline void iSGAngleSystemForce<TBoundaryConditions>::calcAngle(const TBoundaryConditions &boundary, 90 const Angle& currentAngle, 91 const Vector3DBlock* positions, 92 Vector3DBlock* forces, 93 ScalarStructure* energies){ 94 95 int a1 = currentAngle.atom1; 96 int a2 = currentAngle.atom2; 97 int a3 = currentAngle.atom3; 98 99 Real restAngle = currentAngle.restAngle; 100 Real forceConstant = currentAngle.forceConstant; 101 Real ureyBradleyConstant = currentAngle.ureyBradleyConstant; 102 Real ureyBradleyRestLength = currentAngle.ureyBradleyRestLength; 103 104 // Variables for computation of the chemical potential difference 105 Real DeltaK = currentAngle.DeltaK; 106 Real DeltaTheta0 = currentAngle.DeltaTheta0; 107 Real Delta_ubK = currentAngle.Delta_ubK; 108 Real Delta_ubR0 = currentAngle.Delta_ubR0; 109 110 Vector3D atom1((*positions)[a1]); 111 Vector3D atom2((*positions)[a2]); 112 Vector3D atom3((*positions)[a3]); 113 114 //Vector3D r12 = atom1 - atom2; // Vector from atom 1 to atom 2. 115 Vector3D r12(boundary.minimalDifference(atom2,atom1)); 116 //Vector3D r32 = atom3 - atom2; // Vector from atom 3 to atom 2. 117 Vector3D r32(boundary.minimalDifference(atom2,atom3)); 118 //Vector3D r13 = atom1 - atom3; // Vector from atom 1 to atom 3. 119 Vector3D r13(boundary.minimalDifference(atom3,atom1)); 120 Real d12 = r12.norm(); // Distance between atom 1 and 2. 121 Real d32 = r32.norm(); // Distance between atom 3 and 2. 122 Real d13 = r13.norm(); // Distance between atom 1 and 3. 123 124 125 // Calculate theta. 126 Real theta = atan2((r12.cross(r32)).norm(),r12.dot(r32)); 127 Real sinTheta = sin(theta); 128 Real cosTheta = cos(theta); 129 130 // Calculate dpot/dtheta 131 Real dpotdtheta = 2.0 * forceConstant * (theta - restAngle); 132 133 // Calculate dr/dx, dr/dy, dr/dz. 134 Vector3D dr12(r12/d12); 135 Vector3D dr32(r32/d32); 136 Vector3D dr13(r13/d13); 137 // Calulate dtheta/dx, dtheta/dy, dtheta/dz. 138 Vector3D dtheta1((dr12 * cosTheta - dr32)/(sinTheta*d12)); // atom1 139 Vector3D dtheta3((dr32 * cosTheta - dr12)/(sinTheta*d32)); // atom3 140 141 // Calculate Urey Bradley force. 142 Vector3D ureyBradleyforce1(dr13 * (2.0 * ureyBradleyConstant * (d13 - ureyBradleyRestLength))); 143 Vector3D ureyBradleyforce3(-ureyBradleyforce1); 144 145 // Calculate force on atom1 due to atom 2 and 3. 146 Vector3D force1(-dtheta1 * dpotdtheta - ureyBradleyforce1); 147 148 // Calculate force on atom3 due to atom 1 and 2. 149 Vector3D force3(-dtheta3 * dpotdtheta - ureyBradleyforce3); 150 151 // Calculate force on atom2 due to atom 1 and 3. 152 Vector3D force2(-force1 - force3); 153 154 // Add to the total force. 155 (*forces)[a1] += force1; 156 (*forces)[a2] += force2; 157 (*forces)[a3] += force3; 158 159 // Calculate Energy. 160 Real eHarmonic = forceConstant * (theta - restAngle)*(theta - restAngle); 161 Real eUreyBradley = ureyBradleyConstant * (d13 - ureyBradleyRestLength)*(d13 - ureyBradleyRestLength); 162 163 // Calculate Chemical potential difference 164 if (DeltaK != 0.0 && DeltaTheta0 != 0.0) { 165 Real DeltaMuHarmonic = eHarmonic * DeltaK / forceConstant - dpotdtheta * DeltaTheta0; 166 Real DeltaMuUreyBradley = 0.0; 167 if (ureyBradleyConstant != 0) DeltaMuUreyBradley = eUreyBradley * Delta_ubK / ureyBradleyConstant 168 - 2.0 * ureyBradleyConstant * (d13 - ureyBradleyRestLength) * Delta_ubR0; 169 170 // add the chemical potential difference from this angle to the total 171 (*energies)[ScalarStructure::ANGLE_DELTAMU] += DeltaMuHarmonic + DeltaMuUreyBradley; 172 } 173 174 // Add Energy 175 (*energies)[ScalarStructure::ANGLE] += eHarmonic + eUreyBradley; 176 177 //#define ISG_DEBUG_ANGLE 178 #ifdef ISG_DEBUG_ANGLE 179 Report::report.precision(8); 180 Report::report << "Angle atoms: " << a1 << ", " << a2 << ", " << a3 << Report::endr; 181 Report::report << "k = " << forceConstant << ", theta0 = " << restAngle << Report::endr; 182 Report::report << "DeltaK = " << DeltaK << ", DeltaTheta0 = " << DeltaTheta0 << Report::endr; 183 Report::report << "theta = " << theta << Report::endr; 184 Report::report << "current DeltaMu (angle) = " << (*energies)[ScalarStructure::ANGLE_DELTAMU] << Report::endr; 185 Report::report << "current energy (angle) = " << (*energies)[ScalarStructure::ANGLE] << Report::endr; 186 #endif 187 188 // Add virial 189 if(energies->virial()) { 190 Real xy = force1.x * r12.y + force3.x * r32.y; 191 Real xz = force1.x * r12.z + force3.x * r32.z; 192 Real yz = force1.y * r12.z + force3.y * r32.z; 193 (*energies)[ScalarStructure::VIRIALXX] += force1.x * r12.x + force3.x * r32.x; 194 (*energies)[ScalarStructure::VIRIALXY] += xy; 195 (*energies)[ScalarStructure::VIRIALXZ] += xz; 196 (*energies)[ScalarStructure::VIRIALYX] += xy; 197 (*energies)[ScalarStructure::VIRIALYY] += force1.y * r12.y + force3.y * r32.y; 198 (*energies)[ScalarStructure::VIRIALYZ] += yz; 199 (*energies)[ScalarStructure::VIRIALZX] += xz; 200 (*energies)[ScalarStructure::VIRIALZY] += yz; 201 (*energies)[ScalarStructure::VIRIALZZ] += force1.z * r12.z + force3.z * r32.z; 202 } 203 } 204 205 template<class TBoundaryConditions> calcAngleEnergy(const TBoundaryConditions & boundary,const Angle & currentAngle,const Vector3DBlock * positions)206 inline Real iSGAngleSystemForce<TBoundaryConditions>::calcAngleEnergy(const TBoundaryConditions &boundary, 207 const Angle& currentAngle, 208 const Vector3DBlock* positions) 209 { 210 int a1, a2, a3; 211 Real restAngle, forceConstant, ureyBradleyConstant, ureyBradleyRestLength; 212 Vector3D atom1, atom2, atom3, r12, r32, r13; 213 Real d13, theta, eHarmonic, eUreyBradley; 214 215 a1 = currentAngle.atom1; 216 a2 = currentAngle.atom2; 217 a3 = currentAngle.atom3; 218 restAngle = currentAngle.restAngle; 219 forceConstant = currentAngle.forceConstant; 220 ureyBradleyConstant = currentAngle.ureyBradleyConstant; 221 ureyBradleyRestLength = currentAngle.ureyBradleyRestLength; 222 223 atom1 = (*positions)[a1]; 224 atom2 = (*positions)[a2]; 225 atom3 = (*positions)[a3]; 226 227 //r12 = atom1 - atom2; // Vector from atom 1 to atom 2. 228 r12 = boundary.minimalDifference(atom2,atom1); 229 //r32 = atom3 - atom2; // Vector from atom 3 to atom 2. 230 r32 = boundary.minimalDifference(atom2,atom3); 231 //r13 = atom1 - atom3; // Vector from atom 1 to atom 3. 232 r13 = boundary.minimalDifference(atom3,atom1); 233 d13 = r13.norm(); // Distance between atom 1 and 3. 234 235 // Calculate theta. 236 theta = atan2((r12.cross(r32)).norm(),r12.dot(r32)); 237 238 // Calculate Energy. 239 eHarmonic = forceConstant * (theta - restAngle)*(theta - restAngle); 240 eUreyBradley = ureyBradleyConstant * (d13 - ureyBradleyRestLength) 241 *(d13 - ureyBradleyRestLength); 242 243 return (eHarmonic + eUreyBradley); 244 } 245 246 template<class TBoundaryConditions> parallelEvaluate(const GenericTopology * topo,const Vector3DBlock * positions,Vector3DBlock * forces,ScalarStructure * energies)247 inline void iSGAngleSystemForce<TBoundaryConditions>::parallelEvaluate(const GenericTopology* topo, 248 const Vector3DBlock* positions, 249 Vector3DBlock* forces, 250 ScalarStructure* energies) 251 { 252 const TBoundaryConditions &boundary = 253 (dynamic_cast<const SemiGenericTopology<TBoundaryConditions>& >(*topo)).boundaryConditions; 254 255 unsigned int n = topo->angles.size(); 256 unsigned int count = numberOfBlocks(topo,positions); 257 258 for(unsigned int i = 0;i<count;i++){ 259 if(Parallel::next()){ 260 int to = (n*(i+1))/count; 261 if(to > static_cast<int>(n)) 262 to = n; 263 int from = (n*i)/count; 264 for (int j = from; j < to; j++) 265 calcAngle(boundary, topo->angles[j], positions, forces, energies); 266 } 267 } 268 } 269 270 template<class TBoundaryConditions> numberOfBlocks(const GenericTopology * topo,const Vector3DBlock *)271 inline unsigned int iSGAngleSystemForce<TBoundaryConditions>::numberOfBlocks(const GenericTopology* topo, 272 const Vector3DBlock*){ 273 return std::min(Parallel::getAvailableNum(),static_cast<int>(topo->angles.size())); 274 } 275 276 } 277 #endif /* ISGANGLESYSTEMFORCE_H */ 278