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/flucq/FluctuatingChargeNVT.cpp (file contents):
Revision 1760 by gezelter, Thu Jun 21 19:26:46 2012 UTC vs.
Revision 1766 by gezelter, Thu Jul 5 17:08:25 2012 UTC

# 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),
53 <    maxIterNum_(4), thermo(info),
54 <    currentSnapshot_(info->getSnapshotManager()->getCurrentSnapshot()) {
55 <    
56 <    if (hasFlucQ_) {    
53 >    maxIterNum_(4), thermo(info),
54 >    snap(info->getSnapshotManager()->getCurrentSnapshot()) {    
55 >  }
56 >
57 >  void FluctuatingChargeNVT::initialize() {
58 >    FluctuatingChargePropagator::initialize();  
59 >    if (hasFlucQ_) {
60        if (info_->getSimParams()->haveDt()) {
61          dt_ = info_->getSimParams()->getDt();
62          dt2_ = dt_ * 0.5;
# Line 64 | Line 68 | namespace OpenMD {
68        }
69        
70        if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
71 <        currentSnapshot_->setChiElectronic(0.0);
68 <        currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
71 >        snap->setElectronicThermostat(make_pair(0.0, 0.0));
72        }
73        
74        if (!fqParams_->haveTargetTemp()) {
# Line 90 | Line 93 | namespace OpenMD {
93        } else {
94          tauThermostat_ = fqParams_->getTauThermostat();
95        }
96 <      updateSizes();
96 >      updateSizes();
97      }
98    }
99  
97 –  void FluctuatingChargeNVT::initialize() {
98 –    FluctuatingChargePropagator::initialize();
99 –  }
100  
101 –
101    void FluctuatingChargeNVT::moveA() {
102  
103      if (!hasFlucQ_) return;
# Line 109 | 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  
116 –
116      for (mol = info_->beginMolecule(i); mol != NULL;
117           mol = info_->nextMolecule(i)) {
118        for (atom = mol->beginFluctuatingCharge(j); atom != NULL;
# Line 122 | Line 121 | namespace OpenMD {
121          cvel = atom->getFlucQVel();
122          cpos = atom->getFlucQPos();
123          cfrc = atom->getFlucQFrc();
124 <        cmass = atom->getChargeMass();
125 <        
126 <        cerr << atom->getType() << "\n";
128 <        cerr << "Before\n";
129 <        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\tcfrc: " << cfrc << "\tcmass: " << cmass << "\n";
130 <        
124 >        cmass = atom->getChargeMass();      
125 >
126 >        cerr << cpos << "\t" << cvel << "\t" << cfrc << "\t" << instTemp << "\t" <<chi << "\t" << integralOfChidt << "\n";
127          // velocity half step
128          cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel;                    
129          // position whole step
130          cpos += dt_ * cvel;
131          
136 –        cerr << "After\n";
137 –        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\n\n";        
138 –
132          atom->setFlucQVel(cvel);
133          atom->setFlucQPos(cpos);
134        }
# Line 145 | Line 138 | namespace OpenMD {
138        (tauThermostat_ * tauThermostat_);
139  
140      integralOfChidt += chi * dt2_;
141 <    cerr << "Move A instTemp: " << instTemp << "\n";
149 <    currentSnapshot_->setChiElectronic(chi);
150 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
151 <      
141 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
142    }
143  
144    void FluctuatingChargeNVT::updateSizes() {
155 –    if (!hasFlucQ_) return;
145      oldVel_.resize(info_->getNFluctuatingCharges());
146    }
147  
# Line 163 | Line 152 | namespace OpenMD {
152      Molecule* mol;
153      Atom* atom;
154      RealType instTemp;
155 <    RealType chi = currentSnapshot_->getChiElectronic();
155 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
156 >    RealType chi = thermostat.first;
157      RealType oldChi = chi;
158      RealType prevChi;
159 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
159 >    RealType integralOfChidt = thermostat.second;
160      int index;
161      RealType cfrc, cvel, cmass;
162  
# Line 186 | Line 176 | namespace OpenMD {
176      for(int k = 0; k < maxIterNum_; k++) {
177        index = 0;
178        instTemp = thermo.getElectronicTemperature();
189 –      cerr << "MoveB instTemp: " << instTemp << "\n";
179        // evolve chi another half step using the temperature at t + dt/2
180        prevChi = chi;
181        chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) /
# Line 203 | Line 192 | namespace OpenMD {
192            
193            // velocity half step
194            cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index];
206 –          
195            atom->setFlucQVel(cvel);      
196            ++index;          
197          }
# Line 212 | Line 200 | namespace OpenMD {
200          break;
201      }
202      integralOfChidt += dt2_ * chi;
203 <    currentSnapshot_->setChiElectronic(chi);
216 <    currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
203 >    snap->setElectronicThermostat(make_pair(chi, integralOfChidt));
204    }
205    
206    void FluctuatingChargeNVT::resetPropagator() {
207      if (!hasFlucQ_) return;
208 <    currentSnapshot_->setChiElectronic(0.0);
222 <    currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
208 >    snap->setElectronicThermostat(make_pair(0.0, 0.0));
209    }
210    
211    RealType FluctuatingChargeNVT::calcConservedQuantity() {
212      if (!hasFlucQ_) return 0.0;
213 <    RealType chi = currentSnapshot_->getChiElectronic();
214 <    RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
213 >    pair<RealType, RealType> thermostat = snap->getElectronicThermostat();
214 >    RealType chi = thermostat.first;
215 >    RealType integralOfChidt = thermostat.second;
216      RealType fkBT = info_->getNFluctuatingCharges() *
217        PhysicalConstants::kB *targetTemp_;
218  

Diff Legend

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