1 #include "topologyutilities.h" 2 3 #include "Vector3DBlock.h" 4 #include "GenericTopology.h" 5 #include "Report.h" 6 #include "Topology.h" 7 #include "mathutilities.h" 8 #include "pmconstants.h" 9 #include "stringutilities.h" 10 #include "ScalarStructure.h" 11 12 #include <vector> 13 #include <algorithm> 14 #include <set> 15 16 using namespace ProtoMol::Report; 17 using std::string; 18 using std::vector; 19 using std::set; 20 21 namespace ProtoMol { 22 23 //________________________________________________________________randomVelocity randomVelocity(Real temperature,const GenericTopology * topology,Vector3DBlock * velocities,unsigned int seed)24 void randomVelocity(Real temperature, 25 const GenericTopology* topology, 26 Vector3DBlock* velocities, 27 unsigned int seed){ 28 29 // Argument tests 30 if (temperature < 0){ 31 report << error << "[randomVelocity] : " 32 << "Kelvin temperature cannot be negative. " 33 << "Aborting ...."<<endr; 34 return; 35 } 36 37 // Short cuts 38 const Real kbT = temperature*Constant::BOLTZMANN; 39 const unsigned int nAtoms = topology->atoms.size(); 40 41 // Make sure that the velocitie array has the right size ... 42 velocities->resize(nAtoms); 43 44 // Assign the random velocity to each atom 45 for (unsigned int i = 0; i < nAtoms; i++){ 46 Real kbToverM = sqrt(kbT/topology->atoms[i].scaledMass); 47 48 // Generates a random number between -6.0 and 6.0 49 (*velocities)[i].x = kbToverM * randomGaussian(6.0, seed); 50 (*velocities)[i].y = kbToverM * randomGaussian(6.0, seed); 51 (*velocities)[i].z = kbToverM * randomGaussian(6.0, seed); 52 } 53 } 54 55 //________________________________________________________________randomVelocity randomVelocity(Real temperatureFrom,Real temperatureTo,const GenericTopology * topology,Vector3DBlock * velocities,unsigned int seed)56 void randomVelocity(Real temperatureFrom, 57 Real temperatureTo, 58 const GenericTopology* topology, 59 Vector3DBlock* velocities, 60 unsigned int seed){ 61 // Make sure temperatureFrom <= temperatureTo 62 if(temperatureFrom > temperatureTo) 63 std::swap(temperatureFrom,temperatureTo); 64 65 // Argument tests 66 if (temperatureFrom < 0){ 67 report << error << "[randomVelocity] : " 68 << "Kelvin temperature cannot be negative. " 69 << "Aborting ...."<<endr; 70 return; 71 } 72 73 // Random velocities 74 randomVelocity((temperatureFrom+temperatureTo)*0.5,topology,velocities,seed); 75 Real actual = temperature(topology,velocities); 76 77 if(actual <= 0.0 && temperatureFrom > 0.0){ 78 report << error << "[randomVelocity] : " 79 << "Actual temperature zero, can not rescale. " 80 << "Aborting ...."<<endr; 81 return; 82 83 } 84 85 // rescale if necessary 86 if(actual < temperatureFrom || actual > temperatureTo){ 87 Real target = temperatureFrom + (temperatureTo-temperatureFrom) * randomNumber(); 88 Real scale = sqrt( target / actual); 89 const unsigned int nAtoms = topology->atoms.size(); 90 for (unsigned int i = 0; i < nAtoms; i++){ 91 (*velocities)[i] *= scale; 92 } 93 } 94 95 96 } 97 //___________________________________________________________temperature temperature(const GenericTopology * topology,const Vector3DBlock * velocities)98 Real temperature(const GenericTopology* topology, 99 const Vector3DBlock* velocities){ 100 return temperature(kineticEnergy(topology,velocities),topology->degreesOfFreedom); 101 } 102 103 //___________________________________________________________temperature temperature(Real kineticEnergy,unsigned int degreesOfFreedom)104 Real temperature(Real kineticEnergy, 105 unsigned int degreesOfFreedom){ 106 return (2.0 * kineticEnergy / (Constant::BOLTZMANN * degreesOfFreedom)); 107 } 108 109 //___________________________________________________________temperatureForAtomType temperatureForAtomType(const GenericTopology * topology,const Vector3DBlock * velocities,int atomType,waterOption option)110 Real temperatureForAtomType(const GenericTopology* topology, 111 const Vector3DBlock* velocities, 112 int atomType, 113 waterOption option) { 114 // The number of atoms 115 int count; 116 // The total KE 117 Real KE = kineticEnergyForAtomType(topology,velocities,atomType,option,count); 118 119 return temperature(KE,3*count); 120 } 121 122 //___________________________________________________________temperatureForWater temperatureForWater(const GenericTopology * topology,const Vector3DBlock * velocities)123 Real temperatureForWater(const GenericTopology* topology, 124 const Vector3DBlock* velocities) { 125 // The number of waters 126 int count; 127 // The total water KE 128 Real KE = kineticEnergyForWater(topology,velocities,count); 129 130 return temperature(KE,3*count); 131 } 132 133 //___________________________________________________________temperatureForNonWater temperatureForNonWater(const GenericTopology * topology,const Vector3DBlock * velocities)134 Real temperatureForNonWater(const GenericTopology* topology, 135 const Vector3DBlock* velocities) { 136 // The number of non-water molecules 137 int count; 138 // The total non-water KE 139 Real KE = kineticEnergyForNonWater(topology,velocities,count); 140 141 return temperature(KE,3*count); 142 } 143 144 //___________________________________________________________kineticEnergy kineticEnergy(const GenericTopology * topology,const Vector3DBlock * velocities)145 Real kineticEnergy(const GenericTopology* topology, 146 const Vector3DBlock* velocities){ 147 Real kineticEnergy = 0.0; 148 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 149 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 150 } 151 return (kineticEnergy * 0.5); 152 } 153 154 //___________________________________________________________kineticEnergyForAtomType kineticEnergyForAtomType(const GenericTopology * topology,const Vector3DBlock * velocities,int atomType,waterOption option)155 Real kineticEnergyForAtomType(const GenericTopology* topology, 156 const Vector3DBlock* velocities, 157 int atomType, 158 waterOption option){ 159 Real kineticEnergy = 0.0; 160 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 161 if (topology->atoms[i].type == atomType) { 162 bool good = false; 163 if ((option == IGNORE_WATER) && (topology->molecules[topology->atoms[i].molecule].water == false)) { 164 good = true; 165 } 166 else if ((option == ONLY_WATER) && (topology->molecules[topology->atoms[i].molecule].water == true)) { 167 good = true; 168 } 169 else if (option == ALL) { 170 good = true; 171 } 172 173 if (good == true) { 174 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 175 } 176 } 177 } 178 return (kineticEnergy * 0.5); 179 } 180 kineticEnergyForAtomType(const GenericTopology * topology,const Vector3DBlock * velocities,int atomType,waterOption option,int & atomCount)181 Real kineticEnergyForAtomType(const GenericTopology* topology, 182 const Vector3DBlock* velocities, 183 int atomType, 184 waterOption option, 185 int & atomCount){ 186 Real kineticEnergy = 0.0; 187 atomCount = 0; 188 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 189 if (topology->atoms[i].type == atomType) { 190 bool good = false; 191 if ((option == IGNORE_WATER) && (topology->molecules[topology->atoms[i].molecule].water == false)) { 192 good = true; 193 } 194 else if ((option == ONLY_WATER) && (topology->molecules[topology->atoms[i].molecule].water == true)) { 195 good = true; 196 } 197 else if (option == ALL) { 198 good = true; 199 } 200 201 if (good == true) { 202 atomCount++; 203 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 204 } 205 } 206 } 207 return (kineticEnergy * 0.5); 208 } 209 210 //___________________________________________________________kineticEnergyForWater kineticEnergyForWater(const GenericTopology * topology,const Vector3DBlock * velocities)211 Real kineticEnergyForWater(const GenericTopology* topology, 212 const Vector3DBlock* velocities) { 213 Real kineticEnergy = 0.0; 214 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 215 if (topology->molecules[topology->atoms[i].molecule].water == true) { 216 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 217 } 218 } 219 return (kineticEnergy * 0.5); 220 } 221 kineticEnergyForWater(const GenericTopology * topology,const Vector3DBlock * velocities,int & waterCount)222 Real kineticEnergyForWater(const GenericTopology* topology, 223 const Vector3DBlock* velocities, 224 int & waterCount) { 225 Real kineticEnergy = 0.0; 226 waterCount = 0; 227 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 228 if (topology->molecules[topology->atoms[i].molecule].water == true) { 229 waterCount++; 230 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 231 } 232 } 233 return (kineticEnergy * 0.5); 234 } 235 236 //___________________________________________________________kineticEnergyForNonWater kineticEnergyForNonWater(const GenericTopology * topology,const Vector3DBlock * velocities)237 Real kineticEnergyForNonWater(const GenericTopology* topology, 238 const Vector3DBlock* velocities) { 239 Real kineticEnergy = 0.0; 240 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 241 if (topology->molecules[topology->atoms[i].molecule].water == false) { 242 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 243 } 244 } 245 return (kineticEnergy * 0.5); 246 } 247 kineticEnergyForNonWater(const GenericTopology * topology,const Vector3DBlock * velocities,int & nonWaterCount)248 Real kineticEnergyForNonWater(const GenericTopology* topology, 249 const Vector3DBlock* velocities, 250 int & nonWaterCount) { 251 Real kineticEnergy = 0.0; 252 nonWaterCount = 0; 253 for (unsigned int i = 0; i < topology->atoms.size(); i++) { 254 if (topology->molecules[topology->atoms[i].molecule].water == false) { 255 nonWaterCount++; 256 kineticEnergy += topology->atoms[i].scaledMass *((*velocities)[i]).normSquared(); 257 } 258 } 259 return (kineticEnergy * 0.5); 260 } 261 262 //___________________________________________________________molecularKineticEnergy molecularKineticEnergy(const GenericTopology * topology,const Vector3DBlock * velocities)263 Real molecularKineticEnergy(const GenericTopology* topology, 264 const Vector3DBlock* velocities){ 265 266 Real kineticEnergy = 0.0; 267 for (unsigned int i = 0; i < topology->molecules.size(); i++) { 268 269 // compute the COM momentum 270 Vector3D tempMolMom = molecularMomentum(topology->molecules[i].atoms , velocities, topology); 271 272 // compute twice the kinetic energy 273 kineticEnergy += tempMolMom.dot( tempMolMom ) / topology->molecules[i].mass; 274 } 275 276 // multiply by 1/2 to get the kinetic energy 277 return (kineticEnergy * 0.5); 278 } 279 280 //___________________________________________________________velocityVirial velocityVirial(const GenericTopology * topology,const Vector3DBlock * velocities)281 ScalarStructure velocityVirial(const GenericTopology* topology, 282 const Vector3DBlock* velocities){ 283 ScalarStructure res; 284 res.clear(); 285 addVelocityVirial(&res,topology,velocities); 286 return res; 287 } 288 289 //___________________________________________________________addVelocityVirial addVelocityVirial(ScalarStructure * energies,const GenericTopology * topology,const Vector3DBlock * velocities)290 void addVelocityVirial(ScalarStructure* energies, 291 const GenericTopology* topology, 292 const Vector3DBlock* velocities){ 293 for (unsigned int i = 0; i < topology->atoms.size(); i++) 294 energies->addVirial((*velocities)[i],(*velocities)[i] * topology->atoms[i].scaledMass); 295 } 296 297 //___________________________________________________________atomTypeToSymbolName atomTypeToSymbolName(const string & type)298 string atomTypeToSymbolName(const string& type) { 299 if (equalBeginNocase(type, "SOD")) 300 return string("NA"); 301 else if (equalBeginNocase(type, "POT")) 302 return string("K"); 303 else if (equalBeginNocase(type, "CLA")) 304 return string("CL"); 305 else if (equalBeginNocase(type, "CAL")) 306 return string("CA"); 307 else if (equalBeginNocase(type, "CES")) 308 return string("CS"); 309 else if (equalBeginNocase(type, "FE")) 310 return string("FE"); 311 else if (equalBeginNocase(type, "HE")) 312 return string("HE"); 313 else if (equalBeginNocase(type, "NE")) 314 return string("NE"); 315 else if (equalBeginNocase(type, "H")) 316 return string("H"); 317 else if (equalBeginNocase(type, "C")) 318 return string("C"); 319 else if (equalBeginNocase(type, "O")) 320 return string("O"); 321 else if (equalBeginNocase(type, "F")) 322 return string("F"); 323 else if (equalBeginNocase(type, "N")) 324 return string("N"); 325 else if (equalBeginNocase(type, "P")) 326 return string("P"); 327 else if (equalBeginNocase(type, "S")) 328 return string("S"); 329 else if (equalBeginNocase(type, "MG")) 330 return string("MG"); 331 else if (equalBeginNocase(type, "ZN")) 332 return string("ZN"); 333 else 334 return type; 335 } 336 337 338 //___________________________________________________________linearMomentum linearMomentum(const Vector3DBlock * velocities,const GenericTopology * topo)339 Vector3D linearMomentum(const Vector3DBlock *velocities, 340 const GenericTopology *topo){ 341 if(velocities->empty()) 342 return Vector3D(0.0,0.0,0.0); 343 344 // Loop through all the atoms and remove the motion of center of mass 345 // using an advanced method -- Kahan's magic addition algorithm to get 346 // rid of round-off errors: Scientific Computing pp34. 347 348 Vector3D sumMomentum((*velocities)[0]*topo->atoms[0].scaledMass); 349 Vector3D tempC(0.0,0.0,0.0); 350 for (unsigned int i=1; i<topo->atoms.size(); i++) { 351 Vector3D tempX((*velocities)[i]*topo->atoms[i].scaledMass); 352 Vector3D tempY(tempX - tempC); 353 Vector3D tempT(sumMomentum + tempY); 354 tempC = (tempT - sumMomentum) - tempY; 355 sumMomentum = tempT; 356 } 357 return sumMomentum; 358 } 359 360 //___________________________________________________________removeLinearMomentum removeLinearMomentum(Vector3DBlock * velocities,const GenericTopology * topo)361 Vector3D removeLinearMomentum(Vector3DBlock *velocities, 362 const GenericTopology *topo){ 363 if(velocities->empty()) 364 return Vector3D(0.0,0.0,0.0); 365 366 Vector3D momentum = linearMomentum(velocities,topo); 367 Vector3D avgMomentum = momentum/velocities->size(); 368 for (unsigned int i=0; i< topo->atoms.size(); i++){ 369 (*velocities)[i] -= avgMomentum/topo->atoms[i].scaledMass; 370 } 371 372 return momentum; 373 } 374 375 //___________________________________________________________centerOfMass centerOfMass(const Vector3DBlock * positions,const GenericTopology * topo)376 Vector3D centerOfMass(const Vector3DBlock *positions, 377 const GenericTopology *topo){ 378 379 if(positions->empty()) 380 return Vector3D(0.0,0.0,0.0); 381 382 // Loop through all the atoms and compute center of mass 383 // using an advanced method -- Kahan's magic addition algorithm to get 384 // rid of round-off errors: Scientific Computing pp34. 385 // Remove angular momentum 386 387 Real sumM = 0.0; 388 Real m = topo->atoms[0].scaledMass; 389 sumM += m; 390 Vector3D centerOfMass((*positions)[0] * m); 391 Vector3D tempC(0.0,0.0,0.0); 392 for (unsigned int i=1; i<topo->atoms.size(); i++) { 393 m = topo->atoms[i].scaledMass; 394 sumM += m; 395 Vector3D tempX((*positions)[i] * m); 396 Vector3D tempY(tempX - tempC); 397 Vector3D tempT(centerOfMass + tempY); 398 tempC = (tempT - centerOfMass) - tempY; 399 centerOfMass = tempT; 400 } 401 return centerOfMass/sumM; 402 403 404 } 405 406 //___________________________________________________________angularMomentum angularMomentum(const Vector3DBlock * positions,const Vector3DBlock * velocities,const GenericTopology * topo)407 Vector3D angularMomentum(const Vector3DBlock *positions, 408 const Vector3DBlock *velocities, 409 const GenericTopology *topo){ 410 411 if(velocities->empty()) 412 return Vector3D(0.0,0.0,0.0); 413 414 return angularMomentum(positions,velocities,topo,centerOfMass(positions,topo)); 415 } 416 417 //___________________________________________________________angularMomentum angularMomentum(const Vector3DBlock * positions,const Vector3DBlock * velocities,const GenericTopology * topo,const Vector3D & centerOfMass)418 Vector3D angularMomentum(const Vector3DBlock *positions, 419 const Vector3DBlock *velocities, 420 const GenericTopology *topo, 421 const Vector3D& centerOfMass){ 422 423 if(velocities->empty()) 424 return Vector3D(0.0,0.0,0.0); 425 426 // Loop through all the atoms and compute center of mass 427 // using an advanced method -- Kahan's magic addition algorithm to get 428 // rid of round-off errors: Scientific Computing pp34. 429 // Remove angular momentum 430 431 Vector3D momentum(((*positions)[0]-centerOfMass).cross((*velocities)[0]) * topo->atoms[0].scaledMass); 432 Vector3D tempC = Vector3D(0.0,0.0,0.0); 433 for (unsigned int i=1; i<topo->atoms.size(); i++) { 434 Vector3D tempX(((*positions)[i]-centerOfMass).cross((*velocities)[i]) * topo->atoms[i].scaledMass); 435 Vector3D tempY(tempX - tempC); 436 Vector3D tempT(momentum + tempY); 437 tempC = (tempT - momentum) - tempY; 438 momentum = tempT; 439 } 440 return momentum; 441 } 442 443 //___________________________________________________________inertiaMomentum inertiaMomentum(const Vector3DBlock * positions,const GenericTopology * topo,const Vector3D & centerOfMass)444 Matrix3by3 inertiaMomentum(const Vector3DBlock *positions, 445 const GenericTopology *topo, 446 const Vector3D& centerOfMass){ 447 448 Matrix3by3 inertia; 449 inertia.zeroMatrix(); 450 for (unsigned int i=0; i<topo->atoms.size(); i++){ 451 Vector3D r((*positions)[i]-centerOfMass); 452 Real x = r.x; 453 Real y = r.y; 454 Real z = r.z; 455 inertia += Matrix3by3(y*y+z*z, -x*y, -x*z, 456 -x*y, x*x+z*z, -y*z, 457 -x*z, -y*z, x*x+y*y)*topo->atoms[i].scaledMass; 458 } 459 return inertia; 460 } 461 462 //___________________________________________________________removeAngularMomentum removeAngularMomentum(const Vector3DBlock * positions,Vector3DBlock * velocities,const GenericTopology * topo)463 Vector3D removeAngularMomentum(const Vector3DBlock *positions, 464 Vector3DBlock *velocities, 465 const GenericTopology *topo){ 466 467 if(velocities->empty()) 468 return Vector3D(0.0,0.0,0.0); 469 470 // Center of mass 471 Vector3D center = centerOfMass(positions,topo); 472 473 // Angular momentum 474 Vector3D momentum = angularMomentum(positions,velocities,topo,center); 475 if(momentum.normSquared() == 0){ 476 return Vector3D(0.0,0.0,0.0); 477 } 478 479 // Inertia moment 480 Matrix3by3 inertia = inertiaMomentum(positions,topo,center); 481 482 if(!inertia.invert()){ 483 return Vector3D(0.0,0.0,0.0); 484 } 485 // Remove angular momentum 486 Vector3D w(inertia*momentum); 487 for (unsigned int i=0; i< topo->atoms.size(); i++){ 488 Vector3D r((*positions)[i]-center); 489 (*velocities)[i] -= w.cross(r); 490 } 491 492 return momentum; 493 } 494 495 496 //___________________________________________________________computePressure computePressure(const GenericTopology * topology,const Vector3DBlock * positions,const Vector3DBlock * velocities,const ScalarStructure * energies)497 Real computePressure(const GenericTopology* topology, 498 const Vector3DBlock* positions, 499 const Vector3DBlock* velocities, 500 const ScalarStructure* energies){ 501 return computePressure(energies,topology->getVolume(*positions),kineticEnergy(topology,velocities)); 502 } 503 504 //___________________________________________________________computePressure computePressure(const ScalarStructure * energies,Real volume,Real kineticEnergy)505 Real computePressure(const ScalarStructure* energies, 506 Real volume, 507 Real kineticEnergy){ 508 return (energies->pressure(volume) + kineticEnergy*2.0/3.0/volume*Constant::PRESSUREFACTOR); 509 } 510 511 //___________________________________________________________computeMolecularPressure computeMolecularPressure(const ScalarStructure * energies,Real volume,Real kineticEnergy)512 Real computeMolecularPressure(const ScalarStructure* energies, 513 Real volume, 514 Real kineticEnergy){ 515 return (energies->molecularPressure(volume) + kineticEnergy*2.0/3.0/volume*Constant::PRESSUREFACTOR); 516 } 517 518 519 //___________________________________________________________molecularMomentum molecularMomentum(const vector<int> & atomList,const Vector3DBlock * velocities,const GenericTopology * topo)520 Vector3D molecularMomentum(const vector<int> &atomList, 521 const Vector3DBlock *velocities, 522 const GenericTopology *topo){ 523 // Loop through all the atoms and remove the motion of center of mass 524 // using an advanced method -- Kahan's magic addition algorithm to get 525 // rid of round-off errors: Scientific Computing pp34. 526 527 Vector3D sumMomentum((*velocities)[atomList[0]]*topo->atoms[atomList[0]].scaledMass); 528 Vector3D tempC(0.0,0.0,0.0); 529 for (unsigned int i=1; i<atomList.size(); i++) { 530 Vector3D tempX((*velocities)[atomList[i]]*topo->atoms[atomList[i]].scaledMass); 531 Vector3D tempY(tempX - tempC); 532 Vector3D tempT(sumMomentum + tempY); 533 tempC = (tempT - sumMomentum) - tempY; 534 sumMomentum = tempT; 535 } 536 return sumMomentum; 537 } 538 //___________________________________________________________molecularCenterOfMass molecularCenterOfMass(const vector<int> & atomList,const Vector3DBlock * positions,const GenericTopology * topo)539 Vector3D molecularCenterOfMass(const vector<int> &atomList, 540 const Vector3DBlock *positions, 541 const GenericTopology *topo ) { 542 543 if(atomList.empty()) 544 return Vector3D(0.0,0.0,0.0); 545 546 // Loop through all the atoms and compute center of mass using an advanced 547 // method -- Kahan's magic addition algorithm to get rid of round-off 548 // errors: Heath, "Scientific Computing", pp48. 549 550 const int numAtoms = atomList.size(); 551 552 Real sumMass = 0.0; 553 Real mass = topo->atoms[ atomList[0] ].scaledMass; 554 Vector3D center( (*positions)[ atomList[0] ] * mass ); 555 Vector3D tempC( 0.0, 0.0, 0.0 ); 556 sumMass += mass; 557 558 for( int i = 1; i < numAtoms; i++ ) { 559 mass = topo->atoms[ atomList[i] ].scaledMass; 560 sumMass += mass; 561 Vector3D tempX( (*positions)[ atomList[i] ] * mass ); 562 Vector3D tempY( tempX - tempC ); 563 Vector3D tempT( center + tempY ); 564 tempC = ( tempT - center ) - tempY; 565 center = tempT; 566 } 567 568 return( center / sumMass ); 569 570 } 571 572 //___________________________________________________________buildMolecularCenterOfMass buildMolecularCenterOfMass(const Vector3DBlock * positions,GenericTopology * topo)573 void buildMolecularCenterOfMass(const Vector3DBlock *positions, 574 GenericTopology *topo ) { 575 for(unsigned int i=0;i<topo->molecules.size();++i){ 576 topo->molecules[i].position = molecularCenterOfMass(topo->molecules[i].atoms,positions,topo); 577 } 578 } 579 580 //___________________________________________________________buildMolecularMomentum buildMolecularMomentum(const Vector3DBlock * velocities,GenericTopology * topo)581 void buildMolecularMomentum(const Vector3DBlock *velocities, 582 GenericTopology *topo ) { 583 584 for(unsigned int i=0;i<topo->molecules.size();++i){ 585 topo->molecules[i].momentum = molecularMomentum(topo->molecules[i].atoms,velocities,topo); 586 } 587 } 588 589 buildRattleShakeBondConstraintList(GenericTopology * topology,vector<Bond::Constraint> & bondConstraints)590 void buildRattleShakeBondConstraintList(GenericTopology * topology, vector<Bond::Constraint>& bondConstraints){ 591 // here we go through the bond list first, then the angle list to extract the possible 592 // third pair if they are both hydrogen. Thus, all bond lengths plus the third pair in 593 // the waters will be constrained. The angles rather than H-(heavy)-H are left free 594 // of vibrations 595 596 bondConstraints.clear(); 597 598 for (unsigned int i = 0; i < topology->bonds.size(); i++){ 599 bondConstraints.push_back(Bond::Constraint(topology->bonds[i].atom1,topology->bonds[i].atom2,topology->bonds[i].restLength)); 600 } 601 602 // now we have used all the spaces allocated for bondConstraints. 603 // we will use push_back method for more constraints. 604 for (unsigned int i = 0; i < topology->angles.size(); i++){ 605 int a1 = topology->angles[i].atom1; 606 int a2 = topology->angles[i].atom2; 607 int a3 = topology->angles[i].atom3; 608 609 PairIntSorted p1(a1,a2); 610 PairIntSorted p2(a2,a3); 611 612 Real M1 = topology->atoms[a1].scaledMass; 613 Real M3 = topology->atoms[a3].scaledMass; 614 // For regular water M3 = M1 = 1.008, and for heavy waters, 615 // M1 or M3 should be 2 ~ 3 times heavier 616 if((M1<5) && (M3 <5)){ 617 // this exclude (heavy atom)-H pairs and (heavy)-(heavy) pairs 618 // sqrt( 2 * 0.957 * 0.957 * ( 1 - cos( 104.52 * 0.017453 ) ) ); 619 //bondConstraints.push_back(Bond::Constraint(a1,a3,1.51356642665)); 620 621 // ... properly retrieve the right length! 622 int b1 = -1; 623 int b2 = -1; 624 const std::vector< int >& mybonds1 = topology->atoms[a1].mybonds; 625 const std::vector< int >& mybonds3 = topology->atoms[a3].mybonds; 626 627 for (unsigned int j = 0; j < mybonds1.size(); j++){ 628 PairIntSorted s1(topology->bonds[mybonds1[j]].atom1,topology->bonds[mybonds1[j]].atom2); 629 if(p1 == s1){ 630 for (unsigned int k = 0; k < mybonds3.size(); k++){ 631 PairIntSorted s2(topology->bonds[mybonds3[k]].atom1,topology->bonds[mybonds3[k]].atom2); 632 if(p2 == s2){ 633 if(b1 > -1 || b2 > -1){ 634 report << warning << "Found already two bonds ("<<b1<<","<<b2<<") for angle "<<i<<"."<<endr; 635 } 636 else { 637 b1 = mybonds1[j]; 638 b2 = mybonds3[k]; 639 } 640 641 } 642 } 643 } 644 } 645 if(b1 > -1 && b2 > -1){ 646 Real b = topology->bonds[b1].restLength; 647 Real c = topology->bonds[b2].restLength; 648 Real alpha = topology->angles[i].restAngle; 649 bondConstraints.push_back(Bond::Constraint(a1,a3,sqrt(b*b+c*c-2*b*c*cos(alpha)))); 650 report << debug(10) <<i<<" \t :b="<<b<<", c="<<c<<", alpha="<<alpha<<", H-H r="<<sqrt(b*b+c*c-2*b*c*cos(alpha))<<endr; 651 } 652 else { 653 report << debug(10) << "Could not find two matching bonds for angle "<<i<<"."<<endr; 654 } 655 } 656 } 657 if(bondConstraints.size() == topology->bonds.size()) 658 report << hint << "No additional H-X-H SHAKE/RATTLE constraint contributions."<<endr; 659 } 660 661 //___________________________________________________________getAtomsBondedtoDihedral 662 // this function gets all the atoms bonded to ONE side of the dihedral getAtomsBondedtoDihedral(const GenericTopology * topology,set<int,std::less<int>> * atomSet,const int atomID,const int inAtomID,const int outAtomID,const int exclAtomID)663 void getAtomsBondedtoDihedral(const GenericTopology* topology, 664 set< int, std::less< int > >* atomSet, 665 const int atomID, 666 const int inAtomID, 667 const int outAtomID, 668 const int exclAtomID){ 669 670 atomSet->insert(atomID); 671 672 for (unsigned int j = 0; j < topology->atoms[atomID].mybonds.size(); j++) 673 { 674 int bondindex = topology->atoms[atomID].mybonds[j]; 675 int atom1 = topology->bonds[bondindex].atom1; 676 int atom2 = topology->bonds[bondindex].atom2; 677 678 //MUST IMPLEMENT ERROR CASE FOR A BOND LOOP THROUGH THE DIHEDRAL 679 680 if (!((atom1 == inAtomID && atom2 == outAtomID) || (atom1 == outAtomID && atom2 == inAtomID))) 681 { 682 if (atom1 == atomID){ 683 if (atomSet->find(atom2) == atomSet->end()){ 684 getAtomsBondedtoDihedral(topology, atomSet, atom2, inAtomID, outAtomID, exclAtomID); 685 } 686 } 687 else{ 688 if (atomSet->find(atom1) == atomSet->end()){ 689 getAtomsBondedtoDihedral(topology, atomSet, atom1, inAtomID, outAtomID, exclAtomID); 690 } 691 } 692 } 693 694 } 695 } 696 rotateDihedral(const GenericTopology * topology,Vector3DBlock * positions,Vector3DBlock * velocities,const int dihedralID,Real angle)697 void rotateDihedral(const GenericTopology* topology, 698 Vector3DBlock* positions, 699 Vector3DBlock* velocities, 700 const int dihedralID,Real angle){ 701 int exclusionAtomID = topology->dihedrals[dihedralID].atom1; 702 int innerAtomID1 = topology->dihedrals[dihedralID].atom2; 703 int innerAtomID2 = topology->dihedrals[dihedralID].atom3; 704 int atomID = topology->dihedrals[dihedralID].atom4; 705 vector<AngleInfo> angles; 706 707 build_angle_list(topology,atomID,innerAtomID1,innerAtomID2,exclusionAtomID,angle, &angles); 708 general_rotation(innerAtomID1,innerAtomID2,positions,velocities,&angles); 709 710 } 711 rotateDihedral(const GenericTopology * topology,Vector3DBlock * positions,const int dihedralID,Real angle)712 void rotateDihedral(const GenericTopology* topology, 713 Vector3DBlock* positions, 714 const int dihedralID,Real angle){ 715 int exclusionAtomID = topology->dihedrals[dihedralID].atom1; 716 int innerAtomID1 = topology->dihedrals[dihedralID].atom2; 717 int innerAtomID2 = topology->dihedrals[dihedralID].atom3; 718 int atomID = topology->dihedrals[dihedralID].atom4; 719 vector<AngleInfo> angles; 720 721 build_angle_list(topology,atomID,innerAtomID1,innerAtomID2,exclusionAtomID,angle, &angles); 722 general_rotation(innerAtomID1,innerAtomID2,positions,&angles); 723 724 } 725 set_angles(Stack<unsigned int> * nodeStack,vector<AngleInfo> * angles,bool lastIsInnerAtom,Real wholeAngle)726 void set_angles (Stack<unsigned int> *nodeStack,vector<AngleInfo> *angles, bool lastIsInnerAtom, Real wholeAngle) { 727 Real startAngle = 0.0; 728 Real finishAngle = 0.0; 729 bool done = false; 730 unsigned int atomID; 731 unsigned int pos; 732 unsigned int count; 733 734 count = nodeStack->getNumElements(); 735 finishAngle = (*angles)[nodeStack->getElement(count-1)].getAngle(); 736 if (lastIsInnerAtom) { 737 finishAngle = wholeAngle; 738 } 739 pos = nodeStack->getNumElements()-2; 740 while (!done) { 741 atomID = nodeStack->getElement(pos); 742 if ((*angles)[atomID].getAngleType() == AngleInfo::ANGLE_POINTER) { 743 pos = pos - 1; 744 } 745 else if ((*angles)[atomID].getAngleType() == AngleInfo::ANGLE_VALUE) { 746 startAngle = (*angles)[atomID].getAngle(); 747 done = true; 748 } 749 else if (pos > nodeStack->getNumElements()-2) { 750 report << "Error: pos out of range!!!" << endr; 751 } 752 else { 753 report << "Error: No specified Angle Type!!!" << endr; 754 } 755 } 756 757 count = count - pos - 1; 758 finishAngle = startAngle - finishAngle; 759 760 finishAngle = finishAngle/count; 761 762 for (unsigned int x = pos;x < nodeStack->getNumElements(); x++) { 763 atomID = nodeStack->getElement(x); 764 (*angles)[atomID].setAngle(startAngle); 765 startAngle = startAngle - finishAngle; 766 } 767 768 } 769 770 build_angle_list(const GenericTopology * topo,const unsigned int atomID,const unsigned int inAtomID,const unsigned int outAtomID,const unsigned int exclAtomID,Real rotAngle,vector<AngleInfo> * angles)771 void build_angle_list (const GenericTopology *topo, const unsigned int atomID,const unsigned int inAtomID, const unsigned int outAtomID, const unsigned int exclAtomID, Real rotAngle, vector<AngleInfo> *angles) { 772 Stack<unsigned int> nodeStack, indexStack; 773 unsigned int curElement; 774 unsigned int index; 775 unsigned int tempBond1,tempBond2; 776 bool done; 777 AngleInfo temp; 778 779 for (unsigned int x = 0; x < topo->atoms.size(); x++) { 780 angles->push_back(AngleInfo()); 781 } 782 783 for (unsigned int x = 0; x < topo->bonds.size(); x++) { 784 tempBond1 = topo->bonds[x].atom1; 785 tempBond2 = topo->bonds[x].atom2; 786 if (tempBond2 != atomID) { 787 (*angles)[tempBond1].addBond(tempBond2); 788 } 789 if (tempBond1 != atomID) { 790 (*angles)[tempBond2].addBond(tempBond1); 791 } 792 } 793 (*angles)[exclAtomID].setExclusionAtom(); 794 (*angles)[exclAtomID].setAngle(0); 795 (*angles)[inAtomID].setInnerAtom(); 796 (*angles)[inAtomID].setVisited(); 797 (*angles)[outAtomID].setInnerAtom(); 798 (*angles)[outAtomID].setVisited(); 799 800 nodeStack.addElement(atomID); 801 (*angles)[atomID].setAngle(rotAngle); 802 curElement = atomID; 803 index = 0; 804 while (nodeStack.getNumElements() > 0) { 805 done = false; 806 (*angles)[curElement].setVisited(); 807 while (!done) { 808 if (index < ((*angles)[curElement].numBonds())) { 809 if ((*angles)[(*angles)[curElement].getBond(index)].isInnerAtom()) { 810 if (curElement != atomID) { 811 nodeStack.addElement((*angles)[curElement].getBond(index)); 812 set_angles(&nodeStack,angles,true,rotAngle); 813 nodeStack.popElement(); 814 } 815 index++; 816 } 817 else if ((*angles)[(*angles)[curElement].getBond(index)].isExclusionAtom()) { 818 nodeStack.addElement((*angles)[curElement].getBond(index)); 819 set_angles(&nodeStack,angles,false,rotAngle); 820 nodeStack.popElement(); 821 index++; 822 } 823 else if ((*angles)[(*angles)[curElement].getBond(index)].isVisited()) { 824 if ((*angles)[curElement].getBond(index) != nodeStack.getElement(nodeStack.getNumElements()-2)) { 825 nodeStack.addElement((*angles)[curElement].getBond(index)); 826 set_angles(&nodeStack,angles,false,rotAngle); 827 nodeStack.popElement(); 828 } 829 index++; 830 } 831 else { 832 indexStack.addElement(index); 833 (*angles)[(*angles)[curElement].getBond(index)].setPointer(curElement); 834 curElement = (*angles)[curElement].getBond(index); 835 nodeStack.addElement(curElement); 836 done = true; 837 index = 0; 838 } 839 } 840 if ((index > ((*angles)[curElement].numBonds()-1)) || ((*angles)[curElement].numBonds() == 0)) { 841 index = 0; 842 nodeStack.popElement(); 843 if (indexStack.getNumElements() > 0) { 844 index = indexStack.popElement() + 1; 845 } 846 curElement = nodeStack.getElement(nodeStack.getNumElements()-1); 847 done = true; 848 } 849 } 850 } 851 852 done = false; 853 nodeStack.reset(false); 854 index = 0; 855 rotAngle = -1.0; 856 while (!done) { 857 if ((*angles)[index].getAngleType() == AngleInfo::ANGLE_POINTER) { 858 nodeStack.addElement(index); 859 curElement = (*angles)[index].getPointer(); 860 while (nodeStack.getNumElements() > 0) { 861 if ((*angles)[curElement].getAngleType() == AngleInfo::ANGLE_VALUE) { 862 rotAngle = (*angles)[curElement].getAngle(); 863 curElement = nodeStack.popElement(); 864 } 865 else if (((*angles)[curElement].getAngleType() == AngleInfo::ANGLE_POINTER) && (rotAngle != -1.0)) { 866 (*angles)[curElement].setAngle(rotAngle); 867 curElement = nodeStack.popElement(); 868 } 869 else { 870 nodeStack.addElement(curElement); 871 curElement = (*angles)[curElement].getPointer(); 872 } 873 } 874 (*angles)[curElement].setAngle(rotAngle); 875 } 876 index++; 877 if (index >= angles->size()) { 878 done = true; 879 } 880 } 881 882 } 883 884 //___________________________________________________________ computePhiDihedral(const GenericTopology * topo,const Vector3DBlock * positions,int index)885 Real computePhiDihedral(const GenericTopology *topo, const Vector3DBlock *positions, int index){ 886 int a1 = topo->dihedrals[index].atom1; 887 int a2 = topo->dihedrals[index].atom2; 888 int a3 = topo->dihedrals[index].atom3; 889 int a4 = topo->dihedrals[index].atom4; 890 891 //Vector3D r12 = (*positions)[a1] - (*positions)[a2]; // Vector from atom 1 to atom 2 892 Vector3D r12 = topo->minimalDifference((*positions)[a2],(*positions)[a1]); 893 //Vector3D r23 = (*positions)[a2] - (*positions)[a3]; // Vector from atom 2 to atom 3 894 Vector3D r23 = topo->minimalDifference((*positions)[a3],(*positions)[a2]); 895 //Vector3D r34 = (*positions)[a3] - (*positions)[a4]; // Vector from atom 3 to atom 4 896 Vector3D r34 = topo->minimalDifference((*positions)[a4],(*positions)[a3]); 897 898 // Cross product of r12 and r23, represents the plane shared by these two vectors 899 Vector3D a = r12.cross(r23); 900 // Cross product of r12 and r23, represents the plane shared by these two vectors 901 Vector3D b = r23.cross(r34); 902 903 Vector3D c = r23.cross(a); 904 905 // Calculate phi. 906 Real cosPhi = a.dot(b)/(a.norm()*b.norm()); 907 Real sinPhi = c.dot(b)/(c.norm()*b.norm()); 908 909 return -atan2(sinPhi,cosPhi); 910 } 911 912 //___________________________________________________________ computePhiDihedralEnergy(const GenericTopology * topo,int index,Real phi)913 Real computePhiDihedralEnergy(const GenericTopology *topo, int index, Real phi){ 914 Torsion currTorsion = topo->dihedrals[index]; 915 Real energy = 0.0; 916 for( int i = 0; i < currTorsion.multiplicity; i++ ) { 917 if( currTorsion.periodicity[i] > 0 ) { 918 919 // Add energy 920 energy += currTorsion.forceConstant[i] * ( 1.0 + cos( 921 currTorsion.periodicity[i] * phi 922 + currTorsion.phaseShift[i] ) ); 923 //report << "force constant = " << currTorsion.forceConstant[i] << endr; 924 //report << "periodicity = " << currTorsion.periodicity[i] << endr; 925 //report << "trial angle = " << phi << endr; 926 //report << "phase shift = " << currTorsion.phaseShift[i] << endr; 927 } 928 929 else { 930 Real diff = phi - currTorsion.phaseShift[i]; 931 if( diff < -M_PI ) 932 diff += 2 * M_PI; 933 else if( diff > M_PI ) 934 diff -= 2 * M_PI; 935 936 // Add energy 937 energy += currTorsion.forceConstant[i] * diff * diff; 938 } 939 } 940 941 return energy; 942 943 } 944 general_rotation(unsigned int innerAtom1,unsigned int innerAtom2,Vector3DBlock * positions,vector<AngleInfo> * angles)945 void general_rotation(unsigned int innerAtom1, unsigned int innerAtom2, Vector3DBlock *positions, vector<AngleInfo> *angles) { 946 unsigned int numAtoms; 947 Vector3D a; 948 Real norm_a; 949 Real d; 950 Real cosTheta; 951 Real sinTheta; 952 Real temp[9]; 953 Real angle; 954 955 numAtoms = angles->size(); 956 957 a.x = (*positions)[innerAtom2].x - (*positions)[innerAtom1].x; 958 a.y = (*positions)[innerAtom2].y - (*positions)[innerAtom1].y; 959 a.z = (*positions)[innerAtom2].z - (*positions)[innerAtom1].z; 960 norm_a = sqrt(a.x*a.x + a.y*a.y + a.z*a.z); 961 a.x = a.x/norm_a; 962 a.y = a.y/norm_a; 963 a.z = a.z/norm_a; 964 965 d = sqrt(a.y*a.y + a.z*a.z); 966 if (d == 0.0) { 967 for (unsigned int c = 0; c < numAtoms; c++) { 968 angle = (*angles)[c].getAngle(); 969 cosTheta = cos(angle); 970 sinTheta = sin(angle); 971 972 if (0.0 < angle) { 973 temp[1] = (*positions)[c].y - (*positions)[innerAtom1].y; 974 temp[2] = (*positions)[c].z - (*positions)[innerAtom1].z; 975 976 (*positions)[c].y = temp[1]*cosTheta - temp[2]*sinTheta + (*positions)[innerAtom1].y; 977 (*positions)[c].z = temp[1]*sinTheta + temp[2]*cosTheta + (*positions)[innerAtom1].z; 978 } 979 } 980 } 981 else { 982 for (unsigned int c = 0; c < numAtoms; c++) { 983 984 angle = (*angles)[c].getAngle(); 985 cosTheta = cos(angle); 986 sinTheta = sin(angle); 987 988 if (0.0 < angle) { 989 cosTheta = cos(angle); 990 sinTheta = sin(angle); 991 temp[0] = (*positions)[c].x - (*positions)[innerAtom1].x; 992 temp[1] = (*positions)[c].y - (*positions)[innerAtom1].y; 993 temp[8] = (*positions)[c].z - (*positions)[innerAtom1].z;; 994 temp[2] = d*temp[0] + -a.x*((temp[1]*a.y + temp[8]*a.z)/d); 995 temp[3] = (temp[1]*a.z - temp[8]*a.y)/d; 996 temp[4] = (cosTheta*temp[2] - sinTheta*temp[3]); 997 temp[5] = (a.x*temp[0] + (temp[1]*a.y + temp[8]*a.z)); 998 temp[6] = sinTheta*temp[2] + cosTheta*temp[3]; 999 temp[7] = -a.x*temp[4] + d*temp[5]; 1000 1001 (*positions)[c].x = d*temp[4] + a.x*temp[5] + (*positions)[innerAtom1].x; 1002 (*positions)[c].y = (a.z*(temp[6]) + a.y*(temp[7]))/d + (*positions)[innerAtom1].y; 1003 (*positions)[c].z = (-a.y*(temp[6]) + a.z*(temp[7]))/d + (*positions)[innerAtom1].z; 1004 } 1005 } 1006 1007 } 1008 1009 return; 1010 1011 } 1012 1013 general_rotation(unsigned int innerAtom1,unsigned int innerAtom2,Vector3DBlock * positions,Vector3DBlock * velocities,vector<AngleInfo> * angles)1014 void general_rotation(unsigned int innerAtom1, unsigned int innerAtom2, Vector3DBlock *positions, Vector3DBlock *velocities,vector<AngleInfo> *angles) { 1015 unsigned int numAtoms; 1016 Vector3D a; 1017 Real norm_a; 1018 Real d; 1019 Real cosTheta; 1020 Real sinTheta; 1021 Real temp[9]; 1022 Real tempVelo[9]; 1023 Real angle; 1024 1025 numAtoms = angles->size(); 1026 1027 a.x = (*positions)[innerAtom2].x - (*positions)[innerAtom1].x; 1028 a.y = (*positions)[innerAtom2].y - (*positions)[innerAtom1].y; 1029 a.z = (*positions)[innerAtom2].z - (*positions)[innerAtom1].z; 1030 norm_a = sqrt(a.x*a.x + a.y*a.y + a.z*a.z); 1031 a.x = a.x/norm_a; 1032 a.y = a.y/norm_a; 1033 a.z = a.z/norm_a; 1034 1035 d = sqrt(a.y*a.y + a.z*a.z); 1036 1037 if (d == 0.0) { 1038 for (unsigned int c = 0; c < numAtoms; c++) { 1039 angle = (*angles)[c].getAngle(); 1040 cosTheta = cos(angle); 1041 sinTheta = sin(angle); 1042 1043 if (0.0 < angle) { 1044 temp[1] = (*positions)[c].y - (*positions)[innerAtom1].y; 1045 temp[2] = (*positions)[c].z - (*positions)[innerAtom1].z; 1046 1047 tempVelo[1] = (*velocities)[c].y + temp[1]; 1048 tempVelo[2] = (*velocities)[c].z + temp[2]; 1049 1050 (*positions)[c].y = temp[1]*cosTheta - temp[2]*sinTheta; 1051 (*positions)[c].z = temp[1]*sinTheta + temp[2]*cosTheta; 1052 1053 (*velocities)[c].y = tempVelo[1]*cosTheta - tempVelo[2]*sinTheta - (*positions)[c].y; 1054 (*velocities)[c].z = tempVelo[1]*sinTheta + tempVelo[2]*cosTheta - (*positions)[c].z; 1055 1056 (*positions)[c].y += (*positions)[innerAtom1].y; 1057 (*positions)[c].z += (*positions)[innerAtom1].z; 1058 } 1059 } 1060 } 1061 else { 1062 for (unsigned int c = 0; c < numAtoms; c++) { 1063 1064 angle = (*angles)[c].getAngle(); 1065 cosTheta = cos(angle); 1066 sinTheta = sin(angle); 1067 1068 if (0.0 < angle) { 1069 cosTheta = cos(angle); 1070 sinTheta = sin(angle); 1071 temp[0] = (*positions)[c].x - (*positions)[innerAtom1].x; 1072 temp[1] = (*positions)[c].y - (*positions)[innerAtom1].y; 1073 temp[8] = (*positions)[c].z - (*positions)[innerAtom1].z; 1074 temp[2] = d*temp[0] + -a.x*((temp[1]*a.y + temp[8]*a.z)/d); 1075 temp[3] = (temp[1]*a.z - temp[8]*a.y)/d; 1076 temp[4] = (cosTheta*temp[2] - sinTheta*temp[3]); 1077 temp[5] = (a.x*temp[0] + (temp[1]*a.y + temp[8]*a.z)); 1078 temp[6] = sinTheta*temp[2] + cosTheta*temp[3]; 1079 temp[7] = -a.x*temp[4] + d*temp[5]; 1080 1081 tempVelo[0] = (*velocities)[c].x + temp[0]; 1082 tempVelo[1] = (*velocities)[c].y + temp[1]; 1083 tempVelo[8] = (*velocities)[c].z + temp[8]; 1084 tempVelo[2] = d*tempVelo[0] + -a.x*((tempVelo[1]*a.y + tempVelo[8]*a.z)/d); 1085 tempVelo[3] = (tempVelo[1]*a.z - tempVelo[8]*a.y)/d; 1086 tempVelo[4] = (cosTheta*tempVelo[2] - sinTheta*tempVelo[3]); 1087 tempVelo[5] = (a.x*tempVelo[0] + (tempVelo[1]*a.y + tempVelo[8]*a.z)); 1088 tempVelo[6] = sinTheta*tempVelo[2] + cosTheta*tempVelo[3]; 1089 tempVelo[7] = -a.x*tempVelo[4] + d*tempVelo[5]; 1090 1091 (*positions)[c].x = d*temp[4] + a.x*temp[5]; 1092 (*positions)[c].y = (a.z*(temp[6]) + a.y*(temp[7]))/d; 1093 (*positions)[c].z = (-a.y*(temp[6]) + a.z*(temp[7]))/d; 1094 1095 (*velocities)[c].x = d*tempVelo[4] + a.x*tempVelo[5] - (*positions)[c].x; 1096 (*velocities)[c].y = (a.z*(tempVelo[6]) + a.y*(tempVelo[7]))/d - (*positions)[c].y; 1097 (*velocities)[c].z = (-a.y*(tempVelo[6]) + a.z*(tempVelo[7]))/d - (*positions)[c].z; 1098 1099 (*positions)[c].x += (*positions)[innerAtom1].x; 1100 (*positions)[c].y += (*positions)[innerAtom1].y; 1101 (*positions)[c].z += (*positions)[innerAtom1].z; 1102 } 1103 } 1104 } 1105 1106 return; 1107 1108 } 1109 1110 1111 } 1112