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

Comparing:
branches/development/src/flucq/FluctuatingChargeNVT.cpp (file contents), Revision 1735 by jmichalk, Tue Jun 5 17:50:53 2012 UTC vs.
trunk/src/flucq/FluctuatingChargeNVT.cpp (file contents), Revision 2071 by gezelter, Sat Mar 7 21:41:51 2015 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),
53 <    currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) {
54 <
55 <
56 <    if (info_->usesFluctuatingCharges()) {
57 <      if (info_->getNFluctuatingCharges() > 0) {
58 <
59 <        hasFlucQ_ = true;
60 <        Globals* simParams = info_->getSimParams();
61 <        FluctuatingChargeParameters* fqParams = simParams->getFluctuatingChargeParameters();
62 <
63 <        if (simParams->haveDt()) {
64 <          dt_ = simParams->getDt();
65 <          dt2_ = dt_ * 0.5;
66 <        } else {
67 <          sprintf(painCave.errMsg,
68 <                  "FluctuatingChargeNVT Error: dt is not set\n");
69 <          painCave.isFatal = 1;
70 <          simError();
71 <        }
72 <        
73 <        if (!simParams->getUseIntialExtendedSystemState()) {
74 <          currentSnapshot_->setChiElectronic(0.0);
75 <          currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
76 <        }
77 <        
78 <        if (!fqParams->haveTargetTemp()) {
79 <          sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
80 <                  "propagator without a flucQ.targetTemp!\n");
81 <          painCave.isFatal = 1;
82 <          painCave.severity = OPENMD_ERROR;
83 <          simError();
84 <        } else {
85 <          targetTemp_ = fqParams->getTargetTemp();
86 <        }
87 <        
88 <        // We must set tauThermostat.
89 <        
90 <        if (!fqParams->haveTauThermostat()) {
91 <          sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
92 <                  "\tpropagator, you must set flucQ.tauThermostat .\n");
93 <          
94 <          painCave.severity = OPENMD_ERROR;
95 <          painCave.isFatal = 1;
96 <          simError();
97 <        } else {
98 <          tauThermostat_ = fqParams->getTauThermostat();
99 <        }
100 <        updateSizes();
101 <      }
102 <    }      
52 >    FluctuatingChargePropagator(info), maxIterNum_(4), chiTolerance_ (1e-6),
53 >    snap(info->getSnapshotManager()->getCurrentSnapshot()), thermo(info) {  
54    }
55  
56    void FluctuatingChargeNVT::initialize() {
57 <
58 <    if (!hasFlucQ_) return;
59 <
60 <    SimInfo::MoleculeIterator i;
61 <    Molecule::FluctuatingChargeIterator  j;
62 <    Molecule* mol;
63 <    Atom* atom;
64 <    
65 <    for (mol = info_->beginMolecule(i); mol != NULL;
66 <         mol = info_->nextMolecule(i)) {
116 <      for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
117 <           atom = mol->nextFluctuatingCharge(j)) {
118 <        atom->setFlucQPos(0.0);
119 <        atom->setFlucQVel(0.0);
57 >    FluctuatingChargePropagator::initialize();  
58 >    if (hasFlucQ_) {
59 >      if (info_->getSimParams()->haveDt()) {
60 >        dt_ = info_->getSimParams()->getDt();
61 >        dt2_ = dt_ * 0.5;
62 >      } else {
63 >        sprintf(painCave.errMsg,
64 >                "FluctuatingChargeNVT Error: dt is not set\n");
65 >        painCave.isFatal = 1;
66 >        simError();
67        }
68 +      
69 +      if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
70 +        snap->setElectronicThermostat(make_pair(0.0, 0.0));
71 +      }
72 +      
73 +      if (!fqParams_->haveTargetTemp()) {
74 +        sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
75 +                "propagator without a flucQ.targetTemp!\n");
76 +        painCave.isFatal = 1;
77 +        painCave.severity = OPENMD_ERROR;
78 +        simError();
79 +      } else {
80 +        targetTemp_ = fqParams_->getTargetTemp();
81 +      }
82 +      
83 +      // We must set tauThermostat.
84 +      
85 +      if (!fqParams_->haveTauThermostat()) {
86 +        sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
87 +                "\tpropagator, you must set flucQ.tauThermostat .\n");
88 +        
89 +        painCave.severity = OPENMD_ERROR;
90 +        painCave.isFatal = 1;
91 +        simError();
92 +      } else {
93 +        tauThermostat_ = fqParams_->getTauThermostat();
94 +      }
95 +      updateSizes();
96      }
122 –    
123 –    cerr << "Yeah, you should probably implement this\n";
97    }
98  
99 +
100    void FluctuatingChargeNVT::moveA() {
101  
102      if (!hasFlucQ_) return;
# Line 133 | Line 107 | namespace OpenMD {
107      Atom* atom;
108      RealType cvel, cpos, cfrc, cmass;
109  
110 <    RealType chi = currentSnapshot_->getChiElectronic();
111 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
110 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
111 >    RealType chi = thermostat.first;
112 >    RealType integralOfChidt = thermostat.second;
113      RealType instTemp = thermo.getElectronicTemperature();
114  
140 –
115      for (mol = info_->beginMolecule(i); mol != NULL;
116           mol = info_->nextMolecule(i)) {
117        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 146 | Line 120 | namespace OpenMD {
120          cvel = atom->getFlucQVel();
121          cpos = atom->getFlucQPos();
122          cfrc = atom->getFlucQFrc();
123 <        cmass = atom->getChargeMass();
124 <        
151 <        cerr << atom->getType() << "\n";
152 <        cerr << "Before\n";
153 <        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\tcfrc: " << cfrc << "\tcmass: " << cmass << "\n";
154 <        
123 >        cmass = atom->getChargeMass();      
124 >
125          // velocity half step
126          cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel;                    
127          // position whole step
128          cpos += dt_ * cvel;
129          
160 –        cerr << "After\n";
161 –        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\n\n";        
162 –
130          atom->setFlucQVel(cvel);
131          atom->setFlucQPos(cpos);
132        }
# Line 169 | Line 136 | namespace OpenMD {
136        (tauThermostat_ * tauThermostat_);
137  
138      integralOfChidt += chi * dt2_;
139 <    cerr << "Move A instTemp: " << instTemp << "\n";
173 <    currentSnapshot_->setChiElectronic(chi);
174 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
175 <      
139 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
140    }
141  
142    void FluctuatingChargeNVT::updateSizes() {
179 –    if (!hasFlucQ_) return;
143      oldVel_.resize(info_->getNFluctuatingCharges());
144    }
145  
# Line 187 | Line 150 | namespace OpenMD {
150      Molecule* mol;
151      Atom* atom;
152      RealType instTemp;
153 <    RealType chi = currentSnapshot_->getChiElectronic();
153 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
154 >    RealType chi = thermostat.first;
155      RealType oldChi = chi;
156      RealType prevChi;
157 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
157 >    RealType integralOfChidt = thermostat.second;
158      int index;
159      RealType cfrc, cvel, cmass;
160  
# Line 210 | Line 174 | namespace OpenMD {
174      for(int k = 0; k < maxIterNum_; k++) {
175        index = 0;
176        instTemp = thermo.getElectronicTemperature();
213 –      cerr << "MoveB instTemp: " << instTemp << "\n";
177        // evolve chi another half step using the temperature at t + dt/2
178        prevChi = chi;
179        chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) /
# Line 222 | Line 185 | namespace OpenMD {
185               atom = mol->nextFluctuatingCharge(j)) {
186  
187            cfrc = atom->getFlucQFrc();
225 –          cvel =atom->getFlucQVel();
188            cmass = atom->getChargeMass();
189            
190            // velocity half step
191            cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index];
230 –          
192            atom->setFlucQVel(cvel);      
193            ++index;          
194          }
# Line 236 | Line 197 | namespace OpenMD {
197          break;
198      }
199      integralOfChidt += dt2_ * chi;
200 <    currentSnapshot_->setChiElectronic(chi);
240 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
200 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
201    }
202    
203    void FluctuatingChargeNVT::resetPropagator() {
204      if (!hasFlucQ_) return;
205 <    currentSnapshot_->setChiElectronic(0.0);
246 <    currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
205 >    snap->setElectronicThermostat(make_pair(0.0, 0.0));
206    }
207    
208    RealType FluctuatingChargeNVT::calcConservedQuantity() {
209      if (!hasFlucQ_) return 0.0;
210 <    RealType chi = currentSnapshot_->getChiElectronic();
211 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
210 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
211 >    RealType chi = thermostat.first;
212 >    RealType integralOfChidt = thermostat.second;
213      RealType fkBT = info_->getNFluctuatingCharges() *
214        PhysicalConstants::kB *targetTemp_;
215  

Diff Legend

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