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 |
|
|