ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/trunk/src/perturbations/UniformField.cpp
(Generate patch)

Comparing:
branches/development/src/perturbations/ElectricField.cpp (file contents), Revision 1780 by jmarr, Mon Aug 20 18:28:22 2012 UTC vs.
trunk/src/perturbations/ElectricField.cpp (file contents), Revision 1908 by gezelter, Fri Jul 19 21:25:45 2013 UTC

# Line 35 | Line 35
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   */
# Line 77 | Line 77 | namespace OpenMD {
77      Vector3d trq;
78      Vector3d EFfrc;                            
79      Vector3d pos;
80 –    RealType chrg;
81 –    RealType pot, fieldPot, moment;
82 –    RealType chrgToKcal = 23.0609;
83 –    RealType debyeToKcal = 4.8018969509;
84 –    bool isCharge;
80  
81      if (doElectricField) {
82 <      fieldPot = 0.0;
82 >      const RealType chrgToKcal = 23.0609;
83 >      const RealType debyeToKcal = 4.8018969509;
84 >      RealType pot;
85 >      RealType fieldPot = 0.0;
86  
87 <      for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {      
87 >      for (mol = info_->beginMolecule(i); mol != NULL;
88 >           mol = info_->nextMolecule(i)) {      
89 >
90          for (atom = mol->beginAtom(j); atom != NULL;
91               atom = mol->nextAtom(j)) {
92 –          isCharge = false;
93 –          chrg = 0.0;
92  
93 <          FixedChargeAdapter fca = FixedChargeAdapter(atom->getAtomType());
93 >          bool isCharge = false;
94 >          RealType chrg = 0.0;
95 >          
96 >          AtomType* atype = atom->getAtomType();
97 >
98 >          // ad-hoc choice of the origin for potential calculation and
99 >          // fluctuating charge force:
100 >          pos = atom->getPos();
101 >          
102 >          if (atype->isElectrostatic()) {
103 >            atom->addElectricField(EF * chrgToKcal);
104 >          }
105 >          
106 >          FixedChargeAdapter fca = FixedChargeAdapter(atype);
107            if ( fca.isFixedCharge() ) {
108              isCharge = true;
109              chrg = fca.getCharge();
110            }
111            
112 <          FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atom->getAtomType());
112 >          FluctuatingChargeAdapter fqa = FluctuatingChargeAdapter(atype);
113            if ( fqa.isFluctuatingCharge() ) {
114              isCharge = true;
115              chrg += atom->getFlucQPos();
116 +            atom->addFlucQFrc( dot(pos,EF) * chrgToKcal );
117            }
118            
119            if (isCharge) {
120              EFfrc = EF*chrg;
121              EFfrc *= chrgToKcal;
122              atom->addFrc(EFfrc);
111 –            // totally ad-hoc choice of the origin for potential calculation
112 –            pos = atom->getPos();
123              pot = -dot(pos, EFfrc);
124              if (doParticlePot) {      
125                atom->addParticlePot(pot);
# Line 117 | Line 127 | namespace OpenMD {
127              fieldPot += pot;
128            }
129              
130 <          MultipoleAdapter ma = MultipoleAdapter(atom->getAtomType());
130 >          MultipoleAdapter ma = MultipoleAdapter(atype);
131            if (ma.isDipole() ) {
132 <            Vector3d u_i = atom->getElectroFrame().getColumn(2);
133 <            moment = ma.getDipoleMoment();
134 <            moment *= debyeToKcal;
135 <            dip = u_i * moment;
126 <            trq = cross(dip, EF);
127 <            //cerr << "dip = " << dip << "\n";
128 <            // cerr << "trq = " << trq << "\n";
132 >            Vector3d dipole = atom->getDipole();
133 >            dipole *= debyeToKcal;
134 >
135 >            trq = cross(dipole, EF);
136              atom->addTrq(trq);
137 <            pot = -dot(dip, EF);
138 <            //cerr << "pot = " << pot << "\n";
137 >
138 >            pot = -dot(dipole, EF);
139              if (doParticlePot) {      
140                atom->addParticlePot(pot);
141              }
# Line 142 | Line 149 | namespace OpenMD {
149   #endif
150        Snapshot* snap = info_->getSnapshotManager()->getCurrentSnapshot();
151        longRangePotential = snap->getLongRangePotentials();
145 –      // << "longRangePotential = " << longRangePotential << "\n";
152        longRangePotential[ELECTROSTATIC_FAMILY] += fieldPot;
147 –      //cerr << "longRangePotential[ELECTROSTATIC_FAMILY] = " << longRangePotential[ELECTROSTATIC_FAMILY] << "\n";
153        snap->setLongRangePotential(longRangePotential);
154      }
155    }
151 –
156   }

Diff Legend

– Removed lines
+ Added lines
< Changed lines
> Changed lines