--- branches/development/src/perturbations/ElectricField.cpp 2013/01/25 15:21:04 1839 +++ branches/development/src/perturbations/ElectricField.cpp 2013/05/17 14:41:42 1875 @@ -35,7 +35,7 @@ * * [1] Meineke, et al., J. Comp. Chem. 26, 252-271 (2005). * [2] Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006). - * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 24107 (2008). + * [3] Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008). * [4] Kuang & Gezelter, J. Chem. Phys. 133, 164101 (2010). * [5] Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011). */ @@ -77,30 +77,35 @@ namespace OpenMD { Vector3d trq; Vector3d EFfrc; Vector3d pos; - RealType chrg; - RealType pot, fieldPot; - RealType chrgToKcal = 23.0609; - RealType debyeToKcal = 4.8018969509; - bool isCharge; if (doElectricField) { - fieldPot = 0.0; + const RealType chrgToKcal = 23.0609; + const RealType debyeToKcal = 4.8018969509; + RealType pot; + RealType fieldPot = 0.0; - for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) { + for (mol = info_->beginMolecule(i); mol != NULL; + mol = info_->nextMolecule(i)) { + for (atom = mol->beginAtom(j); atom != NULL; atom = mol->nextAtom(j)) { - isCharge = false; - chrg = 0.0; - atom->addElectricField(EF * chrgToKcal); - - FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType()); + bool isCharge = false; + RealType chrg = 0.0; + + AtomType* atype = atom->getAtomType(); + + if (atype->isElectrostatic()) { + atom->addElectricField(EF * chrgToKcal); + } + + FixedChargeAdapter fca = FixedChargeAdapter(atype); if ( fca.isFixedCharge() ) { isCharge = true; chrg = fca.getCharge(); } - FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType()); + FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atype); if ( fqa.isFluctuatingCharge() ) { isCharge = true; chrg += atom->getFlucQPos(); @@ -119,7 +124,7 @@ namespace OpenMD { fieldPot += pot; } - MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType()); + MultipoleAdapter ma = MultipoleAdapter(atype); if (ma.isDipole() ) { Vector3d dipole = atom->getDipole(); dipole *= debyeToKcal;