ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/flucq/FluctuatingChargeNVT.cpp
(Generate patch)

Comparing:
branches/development/src/integrators/FluctuatingChargeNVT.cpp (file contents), Revision 1716 by gezelter, Wed May 23 01:26:15 2012 UTC vs.
branches/development/src/flucq/FluctuatingChargeNVT.cpp (file contents), Revision 1874 by gezelter, Wed May 15 15:09:35 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 45 | Line 45
45   #include "utils/simError.h"
46   #include "utils/PhysicalConstants.hpp"
47  
48 +
49   namespace OpenMD {
50  
51    FluctuatingChargeNVT::FluctuatingChargeNVT(SimInfo* info) :
52 <    FluctuatingChargePropagator(info), chiTolerance_ (1e-6), maxIterNum_(4),
53 <    thermo(info),
54 <    currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) {
54 <
55 <    if (info_->usesFluctuatingCharges()) {
56 <      if (info_->getNFluctuatingCharges() > 0) {
57 <
58 <        hasFlucQ_ = true;
59 <        Globals* simParams = info_->getSimParams();
60 <
61 <        if (simParams->haveDt()) {
62 <          dt_ = simParams->getDt();
63 <          dt2_ = dt_ * 0.5;
64 <        } else {
65 <          sprintf(painCave.errMsg,
66 <                  "FluctuatingChargeNVT Error: dt is not set\n");
67 <          painCave.isFatal = 1;
68 <          simError();
69 <        }
70 <        
71 <        if (!simParams->getUseIntialExtendedSystemState()) {
72 <          currentSnapshot_->setChiElectronic(0.0);
73 <          currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
74 <        }
75 <        
76 <        if (!simParams->haveFlucQTargetTemp()) {
77 <          sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
78 <                  "propagator without a flucQ.targetTemp!\n");
79 <          painCave.isFatal = 1;
80 <          painCave.severity = OPENMD_ERROR;
81 <          simError();
82 <        } else {
83 <          targetTemp_ = simParams->getFlucQTargetTemp();
84 <        }
85 <        
86 <        // We must set tauThermostat.
87 <        
88 <        if (!simParams->haveFlucQtauThermostat()) {
89 <          sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
90 <                  "\tpropagator, you must set flucQ.tauThermostat .\n");
91 <          
92 <          painCave.severity = OPENMD_ERROR;
93 <          painCave.isFatal = 1;
94 <          simError();
95 <        } else {
96 <          tauThermostat_ = simParams->getFlucQtauThermostat();
97 <        }
98 <        updateSizes();
99 <      }
100 <    }      
52 >    FluctuatingChargePropagator(info), chiTolerance_ (1e-6),
53 >    maxIterNum_(4), thermo(info),
54 >    snap(info->getSnapshotManager()->getCurrentSnapshot()) {  
55    }
56  
57    void FluctuatingChargeNVT::initialize() {
58 <
59 <    if (!hasFlucQ_) return;
60 <
61 <    SimInfo::MoleculeIterator i;
62 <    Molecule::FluctuatingChargeIterator  j;
63 <    Molecule* mol;
64 <    Atom* atom;
65 <    
66 <    for (mol = info_->beginMolecule(i); mol != NULL;
67 <         mol = info_->nextMolecule(i)) {
114 <      for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
115 <           atom = mol->nextFluctuatingCharge(j)) {
116 <        atom->setFlucQPos(0.0);
117 <        atom->setFlucQVel(0.0);
58 >    FluctuatingChargePropagator::initialize();  
59 >    if (hasFlucQ_) {
60 >      if (info_->getSimParams()->haveDt()) {
61 >        dt_ = info_->getSimParams()->getDt();
62 >        dt2_ = dt_ * 0.5;
63 >      } else {
64 >        sprintf(painCave.errMsg,
65 >                "FluctuatingChargeNVT Error: dt is not set\n");
66 >        painCave.isFatal = 1;
67 >        simError();
68        }
69 +      
70 +      if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
71 +        snap->setElectronicThermostat(make_pair(0.0, 0.0));
72 +      }
73 +      
74 +      if (!fqParams_->haveTargetTemp()) {
75 +        sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
76 +                "propagator without a flucQ.targetTemp!\n");
77 +        painCave.isFatal = 1;
78 +        painCave.severity = OPENMD_ERROR;
79 +        simError();
80 +      } else {
81 +        targetTemp_ = fqParams_->getTargetTemp();
82 +      }
83 +      
84 +      // We must set tauThermostat.
85 +      
86 +      if (!fqParams_->haveTauThermostat()) {
87 +        sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
88 +                "\tpropagator, you must set flucQ.tauThermostat .\n");
89 +        
90 +        painCave.severity = OPENMD_ERROR;
91 +        painCave.isFatal = 1;
92 +        simError();
93 +      } else {
94 +        tauThermostat_ = fqParams_->getTauThermostat();
95 +      }
96 +      updateSizes();
97      }
120 –    
121 –    cerr << "Yeah, you should probably implement this\n";
98    }
99  
100 +
101    void FluctuatingChargeNVT::moveA() {
102  
103      if (!hasFlucQ_) return;
# Line 131 | Line 108 | namespace OpenMD {
108      Atom* atom;
109      RealType cvel, cpos, cfrc, cmass;
110  
111 <    RealType chi = currentSnapshot_->getChiElectronic();
112 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
111 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
112 >    RealType chi = thermostat.first;
113 >    RealType integralOfChidt = thermostat.second;
114      RealType instTemp = thermo.getElectronicTemperature();
115  
138 –    cerr << "why are we here?\n";
139 –
116      for (mol = info_->beginMolecule(i); mol != NULL;
117           mol = info_->nextMolecule(i)) {
118        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 145 | Line 121 | namespace OpenMD {
121          cvel = atom->getFlucQVel();
122          cpos = atom->getFlucQPos();
123          cfrc = atom->getFlucQFrc();
124 <        cmass = atom->getChargeMass();
125 <        
126 <        // velocity half step
127 <        cvel += dt2_ *PhysicalConstants::energyConvert/cmass*cfrc - dt2_*chi*cvel;                    
124 >        cmass = atom->getChargeMass();      
125 >
126 >        // velocity half step
127 >        cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel;                    
128          // position whole step
129          cpos += dt_ * cvel;
130 <        
130 >        
131          atom->setFlucQVel(cvel);
132          atom->setFlucQPos(cpos);
133        }
# Line 161 | Line 137 | namespace OpenMD {
137        (tauThermostat_ * tauThermostat_);
138  
139      integralOfChidt += chi * dt2_;
140 <    currentSnapshot_->setChiElectronic(chi);
165 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
166 <      
140 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
141    }
142  
143    void FluctuatingChargeNVT::updateSizes() {
170 –    if (!hasFlucQ_) return;
144      oldVel_.resize(info_->getNFluctuatingCharges());
145    }
146  
# Line 178 | Line 151 | namespace OpenMD {
151      Molecule* mol;
152      Atom* atom;
153      RealType instTemp;
154 <    RealType chi = currentSnapshot_->getChiElectronic();
154 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
155 >    RealType chi = thermostat.first;
156      RealType oldChi = chi;
157      RealType prevChi;
158 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
158 >    RealType integralOfChidt = thermostat.second;
159      int index;
160      RealType cfrc, cvel, cmass;
161  
# Line 201 | Line 175 | namespace OpenMD {
175      for(int k = 0; k < maxIterNum_; k++) {
176        index = 0;
177        instTemp = thermo.getElectronicTemperature();
204 –
178        // evolve chi another half step using the temperature at t + dt/2
179        prevChi = chi;
180        chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) /
# Line 213 | Line 186 | namespace OpenMD {
186               atom = mol->nextFluctuatingCharge(j)) {
187  
188            cfrc = atom->getFlucQFrc();
216 –          cvel =atom->getFlucQVel();
189            cmass = atom->getChargeMass();
190            
191            // velocity half step
192 <          cvel = oldVel_[index] + dt2_/cmass*PhysicalConstants::energyConvert * cfrc - dt2_*chi*oldVel_[index];
221 <          
192 >          cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index];
193            atom->setFlucQVel(cvel);      
194            ++index;          
195          }
# Line 227 | Line 198 | namespace OpenMD {
198          break;
199      }
200      integralOfChidt += dt2_ * chi;
201 <    currentSnapshot_->setChiElectronic(chi);
231 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
201 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
202    }
203    
204    void FluctuatingChargeNVT::resetPropagator() {
205      if (!hasFlucQ_) return;
206 <    currentSnapshot_->setChiElectronic(0.0);
237 <    currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
206 >    snap->setElectronicThermostat(make_pair(0.0, 0.0));
207    }
208    
209    RealType FluctuatingChargeNVT::calcConservedQuantity() {
210      if (!hasFlucQ_) return 0.0;
211 <    RealType chi = currentSnapshot_->getChiElectronic();
212 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
211 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
212 >    RealType chi = thermostat.first;
213 >    RealType integralOfChidt = thermostat.second;
214      RealType fkBT = info_->getNFluctuatingCharges() *
215        PhysicalConstants::kB *targetTemp_;
216  

Diff Legend

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