| 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, 24107 (2008).           | 
| 39 | 
< | 
 * [4]  Vardeman & Gezelter, in progress (2009).                         | 
| 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 | 
  | 
#include <cmath> | 
| 46 | 
  | 
#include "utils/simError.h" | 
| 47 | 
  | 
#include "utils/PhysicalConstants.hpp" | 
| 48 | 
  | 
#include "utils/StringUtils.hpp" | 
| 49 | 
+ | 
#ifdef IS_MPI | 
| 50 | 
+ | 
#include <mpi.h> | 
| 51 | 
+ | 
#endif | 
| 52 | 
+ | 
 | 
| 53 | 
  | 
namespace OpenMD { | 
| 54 | 
  | 
  ZconstraintForceManager::ZconstraintForceManager(SimInfo* info): ForceManager(info), infiniteTime(1e31) { | 
| 55 | 
  | 
    currSnapshot_ = info_->getSnapshotManager()->getCurrentSnapshot(); | 
| 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 | 
  | 
 | 
| 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; | 
| 405 | 
  | 
 | 
| 406 | 
  | 
    // apply negative to fixed z-constrained molecues; | 
| 407 | 
  | 
    for ( i = fixedZMols_.begin(); i !=  fixedZMols_.end(); ++i) { | 
| 408 | 
+ | 
 | 
| 409 | 
  | 
      mol = i->mol; | 
| 395 | 
– | 
      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  | 
| 396 | 
– | 
          integrableObject = mol->nextIntegrableObject(ii)) { | 
| 410 | 
  | 
 | 
| 411 | 
< | 
        force[whichDirection] = -getZFOfFixedZMols(mol, integrableObject, i->fz); | 
| 412 | 
< | 
        integrableObject->addFrc(force); | 
| 411 | 
> | 
      for(sd = mol->beginIntegrableObject(ii); sd != NULL;  | 
| 412 | 
> | 
          sd = mol->nextIntegrableObject(ii)) { | 
| 413 | 
> | 
 | 
| 414 | 
> | 
        force[whichDirection] = -getZFOfFixedZMols(mol, sd, i->fz); | 
| 415 | 
> | 
        sd->addFrc(force); | 
| 416 | 
  | 
      } | 
| 417 | 
  | 
    } | 
| 418 | 
  | 
 | 
| 419 | 
  | 
    //modify the forces of moving z-constrained molecules | 
| 420 | 
  | 
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) { | 
| 421 | 
+ | 
 | 
| 422 | 
  | 
      mol = i->mol; | 
| 406 | 
– | 
      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  | 
| 407 | 
– | 
          integrableObject = mol->nextIntegrableObject(ii)) { | 
| 423 | 
  | 
 | 
| 424 | 
+ | 
      for(sd = mol->beginIntegrableObject(ii); sd != NULL;  | 
| 425 | 
+ | 
          sd = mol->nextIntegrableObject(ii)) { | 
| 426 | 
+ | 
 | 
| 427 | 
  | 
        force[whichDirection] = -getZFOfMovingMols(mol,totalFZ); | 
| 428 | 
< | 
        integrableObject->addFrc(force); | 
| 428 | 
> | 
        sd->addFrc(force); | 
| 429 | 
  | 
      } | 
| 430 | 
  | 
    } | 
| 431 | 
  | 
 | 
| 432 | 
  | 
    //modify the forces of unconstrained molecules | 
| 433 | 
  | 
    std::vector<Molecule*>::iterator j; | 
| 434 | 
  | 
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) { | 
| 435 | 
+ | 
 | 
| 436 | 
  | 
      mol =*j; | 
| 418 | 
– | 
      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  | 
| 419 | 
– | 
          integrableObject = mol->nextIntegrableObject(ii)) { | 
| 437 | 
  | 
 | 
| 438 | 
+ | 
      for(sd = mol->beginIntegrableObject(ii); sd != NULL;  | 
| 439 | 
+ | 
          sd = mol->nextIntegrableObject(ii)) { | 
| 440 | 
+ | 
 | 
| 441 | 
  | 
        force[whichDirection] = -getZFOfMovingMols(mol, totalFZ); | 
| 442 | 
< | 
        integrableObject->addFrc(force); | 
| 442 | 
> | 
        sd->addFrc(force); | 
| 443 | 
  | 
      } | 
| 444 | 
  | 
    } | 
| 445 | 
  | 
 | 
| 451 | 
  | 
    Vector3d force(0.0); | 
| 452 | 
  | 
    Vector3d com; | 
| 453 | 
  | 
    RealType totalFZ_local = 0; | 
| 454 | 
+ | 
    RealType lrPot; | 
| 455 | 
  | 
    std::list<ZconstraintMol>::iterator i; | 
| 456 | 
< | 
    StuntDouble* integrableObject; | 
| 456 | 
> | 
    StuntDouble* sd; | 
| 457 | 
  | 
    Molecule::IntegrableObjectIterator ii; | 
| 458 | 
  | 
    Molecule* mol; | 
| 459 | 
  | 
    for ( i = movingZMols_.begin(); i !=  movingZMols_.end(); ++i) { | 
| 462 | 
  | 
      RealType resPos = usingSMD_? i->cantPos : i->param.zTargetPos; | 
| 463 | 
  | 
      RealType diff = com[whichDirection] - resPos;  | 
| 464 | 
  | 
      RealType harmonicU = 0.5 * i->param.kz * diff * diff; | 
| 465 | 
< | 
      currSnapshot_->statData[Stats::LONG_RANGE_POTENTIAL] += harmonicU; | 
| 465 | 
> | 
      lrPot = currSnapshot_->getLongRangePotential(); | 
| 466 | 
> | 
      lrPot += harmonicU; | 
| 467 | 
> | 
      currSnapshot_->setLongRangePotential(lrPot); | 
| 468 | 
  | 
      RealType harmonicF = -i->param.kz * diff; | 
| 469 | 
  | 
      totalFZ_local += harmonicF; | 
| 470 | 
  | 
 | 
| 471 | 
  | 
      //adjust force | 
| 472 | 
< | 
      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  | 
| 473 | 
< | 
          integrableObject = mol->nextIntegrableObject(ii)) { | 
| 472 | 
> | 
      for(sd = mol->beginIntegrableObject(ii); sd != NULL;  | 
| 473 | 
> | 
          sd = mol->nextIntegrableObject(ii)) { | 
| 474 | 
  | 
 | 
| 475 | 
< | 
        force[whichDirection] = getHFOfFixedZMols(mol, integrableObject, harmonicF); | 
| 476 | 
< | 
        integrableObject->addFrc(force);             | 
| 475 | 
> | 
        force[whichDirection] = getHFOfFixedZMols(mol, sd, harmonicF); | 
| 476 | 
> | 
        sd->addFrc(force);             | 
| 477 | 
  | 
      } | 
| 478 | 
  | 
    } | 
| 479 | 
  | 
 | 
| 486 | 
  | 
    //modify the forces of unconstrained molecules | 
| 487 | 
  | 
    std::vector<Molecule*>::iterator j; | 
| 488 | 
  | 
    for ( j = unzconsMols_.begin(); j !=  unzconsMols_.end(); ++j) { | 
| 489 | 
+ | 
 | 
| 490 | 
  | 
      mol = *j; | 
| 467 | 
– | 
      for(integrableObject = mol->beginIntegrableObject(ii); integrableObject != NULL;  | 
| 468 | 
– | 
          integrableObject = mol->nextIntegrableObject(ii)) { | 
| 491 | 
  | 
 | 
| 492 | 
+ | 
      for(sd = mol->beginIntegrableObject(ii); sd != NULL;  | 
| 493 | 
+ | 
          sd = mol->nextIntegrableObject(ii)) { | 
| 494 | 
+ | 
 | 
| 495 | 
  | 
        force[whichDirection] = getHFOfUnconsMols(mol, totalFZ); | 
| 496 | 
< | 
        integrableObject->addFrc(force);             | 
| 496 | 
> | 
        sd->addFrc(force);             | 
| 497 | 
  | 
      } | 
| 498 | 
  | 
    } | 
| 499 | 
  | 
 |