| 6 |  | * redistribute this software in source and binary code form, provided | 
| 7 |  | * that the following conditions are met: | 
| 8 |  | * | 
| 9 | < | * 1. Acknowledgement of the program authors must be made in any | 
| 10 | < | *    publication of scientific results based in part on use of the | 
| 11 | < | *    program.  An acceptable form of acknowledgement is citation of | 
| 12 | < | *    the article in which the program was described (Matthew | 
| 13 | < | *    A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher | 
| 14 | < | *    J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented | 
| 15 | < | *    Parallel Simulation Engine for Molecular Dynamics," | 
| 16 | < | *    J. Comput. Chem. 26, pp. 252-271 (2005)) | 
| 17 | < | * | 
| 18 | < | * 2. Redistributions of source code must retain the above copyright | 
| 9 | > | * 1. Redistributions of source code must retain the above copyright | 
| 10 |  | *    notice, this list of conditions and the following disclaimer. | 
| 11 |  | * | 
| 12 | < | * 3. Redistributions in binary form must reproduce the above copyright | 
| 12 | > | * 2. Redistributions in binary form must reproduce the above copyright | 
| 13 |  | *    notice, this list of conditions and the following disclaimer in the | 
| 14 |  | *    documentation and/or other materials provided with the | 
| 15 |  | *    distribution. | 
| 28 |  | * arising out of the use of or inability to use software, even if the | 
| 29 |  | * University of Notre Dame has been advised of the possibility of | 
| 30 |  | * such damages. | 
| 31 | + | * | 
| 32 | + | * SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your | 
| 33 | + | * research, please cite the appropriate papers when you publish your | 
| 34 | + | * work.  Good starting points are: | 
| 35 | + | * | 
| 36 | + | * [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). | 
| 37 | + | * [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). | 
| 38 | + | * [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). | 
| 39 | + | * [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010). | 
| 40 | + | * [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). | 
| 41 |  | */ | 
| 42 |  |  | 
| 43 | + | #ifdef IS_MPI | 
| 44 | + | #include <mpi.h> | 
| 45 | + | #endif | 
| 46 |  | #include <cmath> | 
| 47 |  | #include "constraints/ZconstraintForceManager.hpp" | 
| 48 |  | #include "integrators/Integrator.hpp" | 
| 49 |  | #include "utils/simError.h" | 
| 50 | < | #include "utils/OOPSEConstant.hpp" | 
| 50 | > | #include "utils/PhysicalConstants.hpp" | 
| 51 |  | #include "utils/StringUtils.hpp" | 
| 52 | < | namespace oopse { | 
| 52 | > |  | 
| 53 | > | namespace OpenMD { | 
| 54 |  | ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) { | 
| 55 |  | currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 56 |  | Globals* simParam = info_->getSimParams(); | 
| 82 |  | zconsTol_ = 0.01; | 
| 83 |  | sprintf(painCave.errMsg, | 
| 84 |  | "ZConstraint Warning: Tolerance for z-constraint method is not specified.\n" | 
| 85 | < | "\tOOPSE will use a default value of %f.\n" | 
| 85 | > | "\tOpenMD will use a default value of %f.\n" | 
| 86 |  | "\tTo set the tolerance, use the zconsTol variable.\n", | 
| 87 |  | zconsTol_); | 
| 88 |  | painCave.isFatal = 0; | 
| 123 |  | } else { | 
| 124 |  | targetTemp = 298.0; | 
| 125 |  | } | 
| 126 | < | RealType zforceConstant = OOPSEConstant::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox); | 
| 126 | > | RealType zforceConstant = PhysicalConstants::kb * targetTemp / (halfOfLargestBox * halfOfLargestBox); | 
| 127 |  |  | 
| 128 |  | int nZconstraints = simParam->getNZconsStamps(); | 
| 129 |  | std::vector<ZConsStamp*> stamp = simParam->getZconsStamps(); | 
| 161 |  | #ifndef IS_MPI | 
| 162 |  | totMassUnconsMols_ = totMassUnconsMols_local; | 
| 163 |  | #else | 
| 164 | < | MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, MPI_REALTYPE, | 
| 165 | < | MPI_SUM, MPI_COMM_WORLD); | 
| 164 | > | MPI_Allreduce(&totMassUnconsMols_local, &totMassUnconsMols_, 1, | 
| 165 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 166 |  | #endif | 
| 167 |  |  | 
| 168 |  | // creat zconsWriter | 
| 245 |  | currZconsTime_ = currSnapshot_->getTime(); | 
| 246 |  | } | 
| 247 |  |  | 
| 248 | < | void ZconstraintForceManager::calcForces(bool needPotential, bool needStress){ | 
| 249 | < | ForceManager::calcForces(needPotential, needStress); | 
| 248 | > | void ZconstraintForceManager::calcForces(){ | 
| 249 | > | ForceManager::calcForces(); | 
| 250 |  |  | 
| 251 |  | if (usingZconsGap_){ | 
| 252 |  | updateZPos(); | 
| 287 |  | Vector3d vel; | 
| 288 |  | std::list<ZconstraintMol>::iterator i; | 
| 289 |  | Molecule* mol; | 
| 290 | < | StuntDouble* integrableObject; | 
| 290 | > | StuntDouble* sd; | 
| 291 |  | Molecule::IntegrableObjectIterator ii; | 
| 292 |  |  | 
| 293 |  | //zero out the velocities of center of mass of fixed z-constrained molecules | 
| 294 |  | for(i = fixedZMols_.begin(); i != fixedZMols_.end(); ++i) { | 
| 295 | + |  | 
| 296 |  | mol = i->mol; | 
| 297 |  | comVel = mol->getComVel(); | 
| 298 | < | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 299 | < | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 300 | < | vel = integrableObject->getVel(); | 
| 298 | > |  | 
| 299 | > | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 300 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 301 | > |  | 
| 302 | > | vel = sd->getVel(); | 
| 303 |  | vel[whichDirection] -= comVel[whichDirection]; | 
| 304 | < | integrableObject->setVel(vel); | 
| 304 | > | sd->setVel(vel); | 
| 305 |  | } | 
| 306 |  | } | 
| 307 |  |  | 
| 326 |  | #ifndef IS_MPI | 
| 327 |  | pzMovingMols = pzMovingMols_local; | 
| 328 |  | #else | 
| 329 | < | MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, MPI_REALTYPE, | 
| 330 | < | MPI_SUM, MPI_COMM_WORLD); | 
| 329 | > | MPI_Allreduce(&pzMovingMols_local, &pzMovingMols, 1, | 
| 330 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 331 |  | #endif | 
| 332 |  |  | 
| 333 |  | RealType vzMovingMols = pzMovingMols / (totMassMovingZMols_ + totMassUnconsMols_); | 
| 334 |  |  | 
| 335 |  | //modify the velocities of moving z-constrained molecuels | 
| 336 |  | for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) { | 
| 337 | + |  | 
| 338 |  | mol = i->mol; | 
| 330 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 331 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 339 |  |  | 
| 340 | < | vel = integrableObject->getVel(); | 
| 340 | > | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 341 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 342 | > |  | 
| 343 | > | vel = sd->getVel(); | 
| 344 |  | vel[whichDirection] -= vzMovingMols; | 
| 345 | < | integrableObject->setVel(vel); | 
| 345 | > | sd->setVel(vel); | 
| 346 |  | } | 
| 347 |  | } | 
| 348 |  |  | 
| 349 |  | //modify the velocites of unconstrained molecules | 
| 350 |  | for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) { | 
| 351 | + |  | 
| 352 |  | mol =*j; | 
| 342 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 343 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 353 |  |  | 
| 354 | < | vel = integrableObject->getVel(); | 
| 354 | > | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 355 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 356 | > |  | 
| 357 | > | vel = sd->getVel(); | 
| 358 |  | vel[whichDirection] -= vzMovingMols; | 
| 359 | < | integrableObject->setVel(vel); | 
| 359 | > | sd->setVel(vel); | 
| 360 |  | } | 
| 361 |  | } | 
| 362 |  |  | 
| 378 |  | //calculate the total z-contrained force of fixed z-contrained molecules | 
| 379 |  | std::list<ZconstraintMol>::iterator i; | 
| 380 |  | Molecule* mol; | 
| 381 | < | StuntDouble* integrableObject; | 
| 381 | > | StuntDouble* sd; | 
| 382 |  | Molecule::IntegrableObjectIterator ii; | 
| 383 |  |  | 
| 384 |  | for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) { | 
| 385 | + |  | 
| 386 |  | mol = i->mol; | 
| 387 |  | i->fz = 0.0; | 
| 375 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 376 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 388 |  |  | 
| 389 | < | force = integrableObject->getFrc(); | 
| 389 | > | for( sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 390 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 391 | > |  | 
| 392 | > | force = sd->getFrc(); | 
| 393 |  | i->fz += force[whichDirection]; | 
| 394 |  | } | 
| 395 |  | totalFZ_local += i->fz; | 
| 397 |  |  | 
| 398 |  | //calculate total z-constraint force | 
| 399 |  | #ifdef IS_MPI | 
| 400 | < | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 400 | > | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, | 
| 401 | > | MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 402 |  | #else | 
| 403 |  | totalFZ = totalFZ_local; | 
| 404 |  | #endif | 
| 406 |  |  | 
| 407 |  | // apply negative to fixed z-constrained molecues; | 
| 408 |  | for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) { | 
| 409 | + |  | 
| 410 |  | mol = i->mol; | 
| 395 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 396 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 411 |  |  | 
| 412 | < | force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz); | 
| 413 | < | integrableObject->addFrc(force); | 
| 412 | > | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 413 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 414 | > |  | 
| 415 | > | force[whichDirection] = -getZFOfFixedZMols(mol, sd, i->fz); | 
| 416 | > | sd->addFrc(force); | 
| 417 |  | } | 
| 418 |  | } | 
| 419 |  |  | 
| 420 |  | //modify the forces of moving z-constrained molecules | 
| 421 |  | for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) { | 
| 422 | + |  | 
| 423 |  | mol = i->mol; | 
| 406 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 407 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 424 |  |  | 
| 425 | + | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 426 | + | sd = mol->nextIntegrableObject(ii)) { | 
| 427 | + |  | 
| 428 |  | force[whichDirection] = -getZFOfMovingMols(mol,totalFZ); | 
| 429 | < | integrableObject->addFrc(force); | 
| 429 | > | sd->addFrc(force); | 
| 430 |  | } | 
| 431 |  | } | 
| 432 |  |  | 
| 433 |  | //modify the forces of unconstrained molecules | 
| 434 |  | std::vector<Molecule*>::iterator j; | 
| 435 |  | for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) { | 
| 436 | + |  | 
| 437 |  | mol =*j; | 
| 418 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 419 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 438 |  |  | 
| 439 | + | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 440 | + | sd = mol->nextIntegrableObject(ii)) { | 
| 441 | + |  | 
| 442 |  | force[whichDirection] = -getZFOfMovingMols(mol, totalFZ); | 
| 443 | < | integrableObject->addFrc(force); | 
| 443 | > | sd->addFrc(force); | 
| 444 |  | } | 
| 445 |  | } | 
| 446 |  |  | 
| 452 |  | Vector3d force(0.0); | 
| 453 |  | Vector3d com; | 
| 454 |  | RealType totalFZ_local = 0; | 
| 455 | + | RealType lrPot; | 
| 456 |  | std::list<ZconstraintMol>::iterator i; | 
| 457 | < | StuntDouble* integrableObject; | 
| 457 | > | StuntDouble* sd; | 
| 458 |  | Molecule::IntegrableObjectIterator ii; | 
| 459 |  | Molecule* mol; | 
| 460 |  | for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) { | 
| 463 |  | RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos; | 
| 464 |  | RealType diff = com[whichDirection] - resPos; | 
| 465 |  | RealType harmonicU = 0.5 * i->param.kz * diff * diff; | 
| 466 | < | currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU; | 
| 466 | > | lrPot = currSnapshot_->getLongRangePotential(); | 
| 467 | > | lrPot += harmonicU; | 
| 468 | > | currSnapshot_->setLongRangePotential(lrPot); | 
| 469 |  | RealType harmonicF = -i->param.kz * diff; | 
| 470 |  | totalFZ_local += harmonicF; | 
| 471 |  |  | 
| 472 |  | //adjust force | 
| 473 | < | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 474 | < | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 473 | > | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 474 | > | sd = mol->nextIntegrableObject(ii)) { | 
| 475 |  |  | 
| 476 | < | force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF); | 
| 477 | < | integrableObject->addFrc(force); | 
| 476 | > | force[whichDirection] = getHFOfFixedZMols(mol, sd, harmonicF); | 
| 477 | > | sd->addFrc(force); | 
| 478 |  | } | 
| 479 |  | } | 
| 480 |  |  | 
| 481 |  | #ifndef IS_MPI | 
| 482 |  | totalFZ = totalFZ_local; | 
| 483 |  | #else | 
| 484 | < | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 484 | > | MPI_Allreduce(&totalFZ_local, &totalFZ, 1, MPI_REALTYPE, | 
| 485 | > | MPI_SUM, MPI_COMM_WORLD); | 
| 486 |  | #endif | 
| 487 |  |  | 
| 488 |  | //modify the forces of unconstrained molecules | 
| 489 |  | std::vector<Molecule*>::iterator j; | 
| 490 |  | for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) { | 
| 491 | + |  | 
| 492 |  | mol = *j; | 
| 467 | – | for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL; | 
| 468 | – | integrableObject = mol->nextIntegrableObject(ii)) { | 
| 493 |  |  | 
| 494 | + | for(sd = mol->beginIntegrableObject(ii); sd != NULL; | 
| 495 | + | sd = mol->nextIntegrableObject(ii)) { | 
| 496 | + |  | 
| 497 |  | force[whichDirection] = getHFOfUnconsMols(mol, totalFZ); | 
| 498 | < | integrableObject->addFrc(force); | 
| 498 | > | sd->addFrc(force); | 
| 499 |  | } | 
| 500 |  | } | 
| 501 |  |  | 
| 567 |  | #ifndef IS_MPI | 
| 568 |  | havingFixed = havingFixed_local; | 
| 569 |  | #else | 
| 570 | < | MPI_Allreduce(&havingFixed_local, &havingFixed, 1, MPI_INT, MPI_SUM, | 
| 571 | < | MPI_COMM_WORLD); | 
| 570 | > | MPI_Allreduce(&havingFixed_local, &havingFixed, 1, | 
| 571 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 572 |  | #endif | 
| 573 |  |  | 
| 574 |  | return havingFixed > 0; | 
| 584 |  | #ifndef IS_MPI | 
| 585 |  | havingMoving = havingMoving_local; | 
| 586 |  | #else | 
| 587 | < | MPI_Allreduce(&havingMoving_local, &havingMoving, 1, MPI_INT, MPI_SUM, | 
| 588 | < | MPI_COMM_WORLD); | 
| 587 | > | MPI_Allreduce(&havingMoving_local, &havingMoving, 1, | 
| 588 | > | MPI_INT, MPI_SUM, MPI_COMM_WORLD); | 
| 589 |  | #endif | 
| 590 |  |  | 
| 591 |  | return havingMoving > 0; | 
| 600 |  | } | 
| 601 |  |  | 
| 602 |  | #ifdef IS_MPI | 
| 603 | < | MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, 1, MPI_REALTYPE, | 
| 604 | < | MPI_SUM, MPI_COMM_WORLD); | 
| 603 | > | MPI_Allreduce(&totMassMovingZMols_local, &totMassMovingZMols_, | 
| 604 | > | 1, MPI_REALTYPE, MPI_SUM, MPI_COMM_WORLD); | 
| 605 |  | #else | 
| 606 |  | totMassMovingZMols_ = totMassMovingZMols_local; | 
| 607 |  | #endif |