57 |
|
#include "primitives/Torsion.hpp" |
58 |
|
#include "primitives/Inversion.hpp" |
59 |
|
#include "nonbonded/NonBondedInteraction.hpp" |
60 |
< |
#include "perturbations/ElectricField.hpp" |
60 |
> |
#include "perturbations/UniformField.hpp" |
61 |
> |
#include "perturbations/UniformGradient.hpp" |
62 |
|
#include "parallel/ForceMatrixDecomposition.hpp" |
63 |
|
|
64 |
|
#include <cstdio> |
395 |
|
switcher_->setSwitch(rSwitch_, rCut_); |
396 |
|
} |
397 |
|
|
397 |
– |
|
398 |
– |
|
399 |
– |
|
398 |
|
void ForceManager::initialize() { |
399 |
|
|
400 |
|
if (!info_->isTopologyDone()) { |
403 |
|
interactionMan_->setSimInfo(info_); |
404 |
|
interactionMan_->initialize(); |
405 |
|
|
406 |
< |
// We want to delay the cutoffs until after the interaction |
407 |
< |
// manager has set up the atom-atom interactions so that we can |
408 |
< |
// query them for suggested cutoff values |
406 |
> |
//! We want to delay the cutoffs until after the interaction |
407 |
> |
//! manager has set up the atom-atom interactions so that we can |
408 |
> |
//! query them for suggested cutoff values |
409 |
|
setupCutoffs(); |
410 |
|
|
411 |
|
info_->prepareTopology(); |
421 |
|
|
422 |
|
ForceFieldOptions& fopts = forceField_->getForceFieldOptions(); |
423 |
|
|
424 |
< |
// Force fields can set options on how to scale van der Waals and |
425 |
< |
// electrostatic interactions for atoms connected via bonds, bends |
426 |
< |
// and torsions in this case the topological distance between |
427 |
< |
// atoms is: |
428 |
< |
// 0 = topologically unconnected |
429 |
< |
// 1 = bonded together |
430 |
< |
// 2 = connected via a bend |
431 |
< |
// 3 = connected via a torsion |
424 |
> |
//! Force fields can set options on how to scale van der Waals and |
425 |
> |
//! electrostatic interactions for atoms connected via bonds, bends |
426 |
> |
//! and torsions in this case the topological distance between |
427 |
> |
//! atoms is: |
428 |
> |
//! 0 = topologically unconnected |
429 |
> |
//! 1 = bonded together |
430 |
> |
//! 2 = connected via a bend |
431 |
> |
//! 3 = connected via a torsion |
432 |
|
|
433 |
|
vdwScale_.reserve(4); |
434 |
|
fill(vdwScale_.begin(), vdwScale_.end(), 0.0); |
446 |
|
electrostaticScale_[2] = fopts.getelectrostatic13scale(); |
447 |
|
electrostaticScale_[3] = fopts.getelectrostatic14scale(); |
448 |
|
|
449 |
< |
if (info_->getSimParams()->haveElectricField()) { |
450 |
< |
ElectricField* eField = new ElectricField(info_); |
449 |
> |
if (info_->getSimParams()->haveUniformField()) { |
450 |
> |
UniformField* eField = new UniformField(info_); |
451 |
|
perturbations_.push_back(eField); |
452 |
|
} |
453 |
< |
|
453 |
> |
if (info_->getSimParams()->haveUniformGradientStrength() || |
454 |
> |
info_->getSimParams()->haveUniformGradientDirection1() || |
455 |
> |
info_->getSimParams()->haveUniformGradientDirection2() ) { |
456 |
> |
UniformGradient* eGrad = new UniformGradient(info_); |
457 |
> |
perturbations_.push_back(eGrad); |
458 |
> |
} |
459 |
> |
|
460 |
|
usePeriodicBoundaryConditions_ = info_->getSimParams()->getUsePeriodicBoundaryConditions(); |
461 |
|
|
462 |
|
fDecomp_->distributeInitialData(); |
745 |
|
idat.f1 = &f1; |
746 |
|
idat.sw = &sw; |
747 |
|
idat.shiftedPot = (cutoffMethod_ == SHIFTED_POTENTIAL) ? true : false; |
748 |
< |
idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; |
748 |
> |
idat.shiftedForce = (cutoffMethod_ == SHIFTED_FORCE || |
749 |
> |
cutoffMethod_ == TAYLOR_SHIFTED) ? true : false; |
750 |
|
idat.doParticlePot = doParticlePot_; |
751 |
|
idat.doElectricField = doElectricField_; |
752 |
|
idat.doSitePotential = doSitePotential_; |
770 |
|
} |
771 |
|
|
772 |
|
for (vector<pair<int, int> >::iterator it = neighborList_.begin(); |
773 |
< |
it != neighborList_.end(); ++it) { |
773 |
> |
it != neighborList_.end(); ++it) { |
774 |
|
|
775 |
|
cg1 = (*it).first; |
776 |
|
cg2 = (*it).second; |
963 |
|
curSnapshot->setLongRangePotential(longRangePotential); |
964 |
|
|
965 |
|
curSnapshot->setExcludedPotentials(*(fDecomp_->getExcludedSelfPotential()) + |
966 |
< |
*(fDecomp_->getExcludedPotential())); |
966 |
> |
*(fDecomp_->getExcludedPotential())); |
967 |
|
|
968 |
|
} |
969 |
|
|
999 |
|
|
1000 |
|
if (info_->getSimParams()->getUseLongRangeCorrections()) { |
1001 |
|
/* |
1002 |
< |
RealType vol = curSnapshot->getVolume(); |
1003 |
< |
RealType Elrc(0.0); |
1004 |
< |
RealType Wlrc(0.0); |
1002 |
> |
RealType vol = curSnapshot->getVolume(); |
1003 |
> |
RealType Elrc(0.0); |
1004 |
> |
RealType Wlrc(0.0); |
1005 |
|
|
1006 |
< |
set<AtomType*>::iterator i; |
1007 |
< |
set<AtomType*>::iterator j; |
1006 |
> |
set<AtomType*>::iterator i; |
1007 |
> |
set<AtomType*>::iterator j; |
1008 |
|
|
1009 |
< |
RealType n_i, n_j; |
1010 |
< |
RealType rho_i, rho_j; |
1011 |
< |
pair<RealType, RealType> LRI; |
1009 |
> |
RealType n_i, n_j; |
1010 |
> |
RealType rho_i, rho_j; |
1011 |
> |
pair<RealType, RealType> LRI; |
1012 |
|
|
1013 |
< |
for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { |
1013 |
> |
for (i = atomTypes_.begin(); i != atomTypes_.end(); ++i) { |
1014 |
|
n_i = RealType(info_->getGlobalCountOfType(*i)); |
1015 |
|
rho_i = n_i / vol; |
1016 |
|
for (j = atomTypes_.begin(); j != atomTypes_.end(); ++j) { |
1017 |
< |
n_j = RealType(info_->getGlobalCountOfType(*j)); |
1018 |
< |
rho_j = n_j / vol; |
1017 |
> |
n_j = RealType(info_->getGlobalCountOfType(*j)); |
1018 |
> |
rho_j = n_j / vol; |
1019 |
|
|
1020 |
< |
LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); |
1020 |
> |
LRI = interactionMan_->getLongRangeIntegrals( (*i), (*j) ); |
1021 |
|
|
1022 |
< |
Elrc += n_i * rho_j * LRI.first; |
1023 |
< |
Wlrc -= rho_i * rho_j * LRI.second; |
1022 |
> |
Elrc += n_i * rho_j * LRI.first; |
1023 |
> |
Wlrc -= rho_i * rho_j * LRI.second; |
1024 |
|
} |
1025 |
< |
} |
1026 |
< |
Elrc *= 2.0 * NumericConstant::PI; |
1027 |
< |
Wlrc *= 2.0 * NumericConstant::PI; |
1025 |
> |
} |
1026 |
> |
Elrc *= 2.0 * NumericConstant::PI; |
1027 |
> |
Wlrc *= 2.0 * NumericConstant::PI; |
1028 |
|
|
1029 |
< |
RealType lrp = curSnapshot->getLongRangePotential(); |
1030 |
< |
curSnapshot->setLongRangePotential(lrp + Elrc); |
1031 |
< |
stressTensor += Wlrc * SquareMatrix3<RealType>::identity(); |
1032 |
< |
curSnapshot->setStressTensor(stressTensor); |
1029 |
> |
RealType lrp = curSnapshot->getLongRangePotential(); |
1030 |
> |
curSnapshot->setLongRangePotential(lrp + Elrc); |
1031 |
> |
stressTensor += Wlrc * SquareMatrix3<RealType>::identity(); |
1032 |
> |
curSnapshot->setStressTensor(stressTensor); |
1033 |
|
*/ |
1034 |
|
|
1035 |
|
} |