| 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, 24107 (2008). |
| 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 |
|
*/ |
| 127 |
|
} |
| 128 |
|
|
| 129 |
|
void FluctuatingChargeLangevin::applyConstraints() { |
| 130 |
+ |
|
| 131 |
+ |
if (!hasFlucQ_) return; |
| 132 |
+ |
|
| 133 |
|
SimInfo::MoleculeIterator i; |
| 134 |
|
Molecule::FluctuatingChargeIterator j; |
| 135 |
|
Molecule* mol; |
| 136 |
|
Atom* atom; |
| 137 |
< |
RealType cvel, cpos, cfrc, cmass, randomForce, frictionForce; |
| 137 |
> |
RealType cvel, cfrc, cmass, randomForce, frictionForce; |
| 138 |
|
RealType velStep, oldFF; // used to test for convergence |
| 139 |
|
|
| 140 |
|
for (mol = info_->beginMolecule(i); mol != NULL; |
| 141 |
|
mol = info_->nextMolecule(i)) { |
| 142 |
|
for (atom = mol->beginFluctuatingCharge(j); atom != NULL; |
| 143 |
|
atom = mol->nextFluctuatingCharge(j)) { |
| 144 |
< |
|
| 142 |
< |
cvel = atom->getFlucQVel(); |
| 143 |
< |
cpos = atom->getFlucQPos(); |
| 144 |
< |
cfrc = atom->getFlucQFrc(); |
| 145 |
< |
cmass = atom->getChargeMass(); |
| 146 |
< |
|
| 144 |
> |
|
| 145 |
|
randomForce = randNumGen_.randNorm(0, variance_ ); |
| 146 |
|
atom->addFlucQFrc(randomForce); |
| 147 |
|
|
| 149 |
|
// required is at the full step: v(t + h), while we have |
| 150 |
|
// initially the velocity at the half step: v(t + h/2). We |
| 151 |
|
// need to iterate to converge the friction force vector. |
| 152 |
< |
|
| 152 |
> |
|
| 153 |
|
// this is the velocity at the half-step: |
| 154 |
< |
|
| 154 |
> |
|
| 155 |
|
cvel = atom->getFlucQVel(); |
| 156 |
< |
|
| 156 |
> |
|
| 157 |
|
// estimate velocity at full-step using everything but |
| 158 |
|
// friction forces: |
| 159 |
|
|
| 160 |
|
cfrc = atom->getFlucQFrc(); |
| 161 |
+ |
cmass = atom->getChargeMass(); |
| 162 |
|
velStep = cvel + dt2_ * cfrc / cmass; |
| 163 |
< |
|
| 163 |
> |
|
| 164 |
|
frictionForce = 0.0; |
| 165 |
+ |
|
| 166 |
|
//iteration starts here: |
| 167 |
|
|
| 168 |
|
for (int k = 0; k < maxIterNum_; k++) { |
