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 1760 by gezelter, Thu Jun 21 19:26:46 2012 UTC

# Line 48 | Line 48 | namespace OpenMD {
48   namespace OpenMD {
49  
50    FluctuatingChargeNVT::FluctuatingChargeNVT(SimInfo* info) :
51 <    FluctuatingChargePropagator(info), chiTolerance_ (1e-6), maxIterNum_(4),
52 <    thermo(info),
51 >    FluctuatingChargePropagator(info), chiTolerance_ (1e-6),
52 >    maxIterNum_(4), thermo(info),
53      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 <        }
54 >    
55 >    if (hasFlucQ_) {    
56 >      if (info_->getSimParams()->haveDt()) {
57 >        dt_ = info_->getSimParams()->getDt();
58 >        dt2_ = dt_ * 0.5;
59 >      } else {
60 >        sprintf(painCave.errMsg,
61 >                "FluctuatingChargeNVT Error: dt is not set\n");
62 >        painCave.isFatal = 1;
63 >        simError();
64 >      }
65 >      
66 >      if (!info_->getSimParams()->getUseIntialExtendedSystemState()) {
67 >        currentSnapshot_->setChiElectronic(0.0);
68 >        currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
69 >      }
70 >      
71 >      if (!fqParams_->haveTargetTemp()) {
72 >        sprintf(painCave.errMsg, "You can't use the FluctuatingChargeNVT "
73 >                "propagator without a flucQ.targetTemp!\n");
74 >        painCave.isFatal = 1;
75 >        painCave.severity = OPENMD_ERROR;
76 >        simError();
77 >      } else {
78 >        targetTemp_ = fqParams_->getTargetTemp();
79 >      }
80 >      
81 >      // We must set tauThermostat.
82 >      
83 >      if (!fqParams_->haveTauThermostat()) {
84 >        sprintf(painCave.errMsg, "If you use the FluctuatingChargeNVT\n"
85 >                "\tpropagator, you must set flucQ.tauThermostat .\n");
86          
87 <        if (!simParams->getUseIntialExtendedSystemState()) {
88 <          currentSnapshot_->setChiElectronic(0.0);
89 <          currentSnapshot_->setIntegralOfChiElectronicDt(0.0);
90 <        }
91 <        
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();
87 >        painCave.severity = OPENMD_ERROR;
88 >        painCave.isFatal = 1;
89 >        simError();
90 >      } else {
91 >        tauThermostat_ = fqParams_->getTauThermostat();
92        }
93 <    }      
93 >      updateSizes();
94 >    }
95    }
96  
97    void FluctuatingChargeNVT::initialize() {
98 <
105 <    if (!hasFlucQ_) return;
106 <
107 <    SimInfo::MoleculeIterator i;
108 <    Molecule::FluctuatingChargeIterator  j;
109 <    Molecule* mol;
110 <    Atom* atom;
111 <    
112 <    for (mol = info_->beginMolecule(i); mol != NULL;
113 <         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);
118 <      }
119 <    }
120 <    
121 <    cerr << "Yeah, you should probably implement this\n";
98 >    FluctuatingChargePropagator::initialize();
99    }
100  
101 +
102    void FluctuatingChargeNVT::moveA() {
103  
104      if (!hasFlucQ_) return;
# Line 135 | Line 113 | namespace OpenMD {
113      RealType integralOfChidt = currentSnapshot_->getIntegralOfChiElectronicDt();
114      RealType instTemp = thermo.getElectronicTemperature();
115  
138 –    cerr << "why are we here?\n";
116  
117      for (mol = info_->beginMolecule(i); mol != NULL;
118           mol = info_->nextMolecule(i)) {
# Line 147 | Line 124 | namespace OpenMD {
124          cfrc = atom->getFlucQFrc();
125          cmass = atom->getChargeMass();
126          
127 <        // velocity half step
128 <        cvel += dt2_ *PhysicalConstants::energyConvert/cmass*cfrc - dt2_*chi*cvel;                    
127 >        cerr << atom->getType() << "\n";
128 >        cerr << "Before\n";
129 >        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\tcfrc: " << cfrc << "\tcmass: " << cmass << "\n";
130 >        
131 >        // velocity half step
132 >        cvel += dt2_ * cfrc / cmass - dt2_*chi*cvel;                    
133          // position whole step
134          cpos += dt_ * cvel;
135 <        
135 >        
136 >        cerr << "After\n";
137 >        cerr << "cvel: " << cvel << "\tcpos: " << cpos << "\n\n";        
138 >
139          atom->setFlucQVel(cvel);
140          atom->setFlucQPos(cpos);
141        }
# Line 161 | Line 145 | namespace OpenMD {
145        (tauThermostat_ * tauThermostat_);
146  
147      integralOfChidt += chi * dt2_;
148 +    cerr << "Move A instTemp: " << instTemp << "\n";
149      currentSnapshot_->setChiElectronic(chi);
150      currentSnapshot_->setIntegralOfChiElectronicDt(integralOfChidt);
151        
# Line 201 | Line 186 | namespace OpenMD {
186      for(int k = 0; k < maxIterNum_; k++) {
187        index = 0;
188        instTemp = thermo.getElectronicTemperature();
189 <
189 >      cerr << "MoveB instTemp: " << instTemp << "\n";
190        // evolve chi another half step using the temperature at t + dt/2
191        prevChi = chi;
192        chi = oldChi + dt2_ * (instTemp / targetTemp_ - 1.0) /
# Line 217 | Line 202 | namespace OpenMD {
202            cmass = atom->getChargeMass();
203            
204            // velocity half step
205 <          cvel = oldVel_[index] + dt2_/cmass*PhysicalConstants::energyConvert * cfrc - dt2_*chi*oldVel_[index];
205 >          cvel = oldVel_[index] + dt2_ * cfrc / cmass - dt2_*chi*oldVel_[index];
206            
207            atom->setFlucQVel(cvel);      
208            ++index;          

Diff Legend

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