ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new_design/OOPSE-3.0/src/integrators/NVT.cpp
(Generate patch)

Comparing branches/new_design/OOPSE-3.0/src/integrators/NVT.cpp (file contents):
Revision 1762 by tim, Fri Nov 19 21:38:22 2004 UTC vs.
Revision 1774 by tim, Tue Nov 23 23:12:23 2004 UTC

# Line 1 | Line 1
1 < #include <math.h>
1 > #inlcude "integrators/NVT.hpp"
2  
3 < #include "primitives/Atom.hpp"
4 < #include "primitives/SRI.hpp"
5 < #include "primitives/AbstractClasses.hpp"
6 < #include "brains/SimInfo.hpp"
7 < #include "UseTheForce/ForceFields.hpp"
8 < #include "brains/Thermo.hpp"
9 < #include "io/ReadWrite.hpp"
10 < #include "integrators/Integrator.hpp"
11 < #include "utils/simError.h"
3 > namespace oopse {
4  
5 < // Basic thermostating via Hoover, Phys.Rev.A, 1985, Vol. 31 (5) 1695-1697
5 > NVT::NVT(SimInfo* info) : VelocityVerletIntegrator(info), chiTolerance_ (1e-6) {
6  
7 < template <typename T>NVT<T>::NVT(SimInfo *theInfo, ForceFields *the_ff) :
16 <    T(theInfo, the_ff) {
17 <    GenericData *data;
18 <    DoubleGenericData *chiValue;
19 <    DoubleGenericData *integralOfChidtValue;
7 >    Globals* globals = info_->getGlobals();
8  
9 <    chiValue = NULL;
10 <    integralOfChidtValue = NULL;
9 >    if (globals->getUseInitXSstate()) {
10 >        Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
11 >        currSnapshot->setChi(0.0);
12 >        currSnapshot->setIntegralOfChiDt(0.0);
13 >    }
14 >    
15 >    if (!globals->haveTargetTemp()) {
16 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
17 >        painCave.isFatal = 1;
18 >        painCave.severity = OOPSE_ERROR;
19 >        simError();
20 >    } else {
21 >        targetTemp_ = globals->getTargetTemp();
22 >    }
23  
24 <    chi = 0.0;
25 <    have_tau_thermostat = 0;
26 <    have_target_temp = 0;
27 <    have_chi_tolerance = 0;
28 <    integralOfChidt = 0.0;
24 >    // We must set tauThermostat_.
25  
26 <    if (theInfo->useInitXSstate) {
26 >    if (!globals->haveTauThermostat()) {
27 >        sprintf(painCave.errMsg, "If you use the constant temperature\n"
28 >                                     "\tintegrator, you must set tauThermostat_.\n");
29  
30 <        // retrieve chi and integralOfChidt from simInfo
31 <        data = info->getPropertyByName(CHIVALUE_ID);
32 <
33 <        if (data) {
34 <            chiValue = dynamic_cast<DoubleGenericData *>(data);
37 <        }
38 <
39 <        data = info->getPropertyByName(INTEGRALOFCHIDT_ID);
40 <
41 <        if (data) {
42 <            integralOfChidtValue = dynamic_cast<DoubleGenericData *>(data);
43 <        }
44 <
45 <        // chi and integralOfChidt should appear by pair
46 <        if (chiValue && integralOfChidtValue) {
47 <            chi = chiValue->getData();
48 <            integralOfChidt = integralOfChidtValue->getData();
49 <        }
30 >        painCave.severity = OOPSE_ERROR;
31 >        painCave.isFatal = 1;
32 >        simError();
33 >    } else {
34 >        tauThermostat_ = globals->getTauThermostat();
35      }
36  
37 <    oldVel = new double[3 * integrableObjects.size()];
53 <    oldJi = new double[3 * integrableObjects.size()];
37 >    update();
38   }
39  
40 < template <typename T>NVT<T>::~NVT() {
57 <    delete [] oldVel;
58 <    delete [] oldJi;
40 > NVT::~NVT() {
41   }
42  
43 < template <typename T>
44 < void NVT<T>::moveA() {
45 <    int i, j;
46 <    DirectionalAtom *dAtom;
43 > void NVT::update() {
44 >    oldVel_.resize(info_->getNIntegrableObjects());
45 >    oldJi_.resize(info_->getNIntegrableObjects());    
46 > }
47 > void NVT::moveA() {
48 >    typename SimInfo::MoleculeIterator i;
49 >    typename Molecule::IntegrableObjectIterator  j;
50 >    Molecule* mol;
51 >    StuntDouble* integrableObject;
52      Vector3d Tb;
53      Vector3d ji;
54      double mass;
# Line 75 | Line 62 | void NVT<T>::moveA() {
62  
63      instTemp = tStats->getTemperature();
64  
65 <    for(i = 0; i < integrableObjects.size(); i++) {
66 <        vel = integrableObjects[i]->getVel();
67 <        pos = integrableObjects[i]->getPos();
81 <        integrableObjects[i]->getFrc(frc);
65 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
66 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
67 >               integrableObject = mol->nextIntegrableObject(j)) {
68  
69 <        mass = integrableObjects[i]->getMass();
69 >        vel = integrableObject->getVel();
70 >        pos = integrableObject->getPos();
71 >        frc = integrableObject->getFrc();
72  
73 <        for(j = 0; j < 3; j++) {
86 <            // velocity half step  (use chi from previous step here):
87 <            vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
73 >        mass = integrableObject->getMass();
74  
75 <            // position whole step
76 <            pos[j] += dt * vel[j];
77 <        }
75 >        // velocity half step  (use chi from previous step here):
76 >        //vel[j] += dt2 * ((frc[j] / mass ) * eConvert - vel[j]*chi);
77 >        vel = dt2 *eConvert/mass*frc - dt2*chi*vel;
78 >        
79 >        // position whole step
80 >        //pos[j] += dt * vel[j];
81 >        pos += dt * vel;
82  
83 <        integrableObjects[i]->setVel(vel);
84 <        integrableObjects[i]->setPos(pos);
83 >        integrableObject->setVel(vel);
84 >        integrableObject->setPos(pos);
85  
86 <        if (integrableObjects[i]->isDirectional()) {
86 >        if (integrableObject->isDirectional()) {
87  
88              // get and convert the torque to body frame
89  
90 <            Tb = integrableObjects[i]->getTrq();
91 <            integrableObjects[i]->lab2Body(Tb);
90 >            Tb = integrableObject->getTrq();
91 >            integrableObject->lab2Body(Tb);
92  
93              // get the angular momentum, and propagate a half step
94  
95 <            ji = integrableObjects[i]->getJ();
95 >            ji = integrableObject->getJ();
96  
97 <            for(j = 0; j < 3; j++)
98 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
97 >            //ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
98 >            ji += dt2*eConvert*Tb - dt2*ch *ji;
99 >            this->rotationPropagation(integrableObject, ji);
100  
101 <            this->rotationPropagation(integrableObjects[i], ji);
111 <
112 <            integrableObjects[i]->setJ(ji);
101 >            integrableObject->setJ(ji);
102          }
103      }
104  
105 <    if (nConstrained)
106 <        constrainA();
105 >    }
106 >    
107 >    //constraintAlgorithm->doConstrainA();
108  
109      // Finally, evolve chi a half step (just like a velocity) using
110      // temperature at time t, not time t+dt/2
111  
112 <    //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
113 <
114 <    chi += dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
112 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
113 >    double chi = currSnapshot->getChi();
114 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
115 >    
116 >    chi += dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
117      integralOfChidt += chi * dt2;
118 +
119 +    currSnapshot->setChi(chi);
120 +    currSnapshot->setIntegralOfChiDt(integralOfChidt);
121   }
122  
123 < template <typename T>
124 < void NVT<T>::moveB(void) {
125 <    int i, j, k;
126 <    double Tb[3], ji[3];
127 <    double vel[3], frc[3];
123 > void NVT::moveB() {
124 >    typename SimInfo::MoleculeIterator i;
125 >    typename Molecule::IntegrableObjectIterator  j;
126 >    Molecule* mol;
127 >    StuntDouble* integrableObject;
128 >    
129 >    Vector3d Tb;
130 >    Vector3d ji;    
131 >    Vector3d vel;
132 >    Vector3d frc;
133      double mass;
134      double instTemp;
135      double oldChi, prevChi;
136 <
136 >    int index;
137      // Set things up for the iteration:
138  
139 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
140 +    double chi = currSnapshot->getChi();
141 +    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
142      oldChi = chi;
143  
144 <    for(i = 0; i < integrableObjects.size(); i++) {
145 <        vel = integrableObjects[i]->getVel();
146 <
147 <        for(j = 0; j < 3; j++)
148 <    oldVel[3 * i + j] = vel[j];
149 <
147 <        if (integrableObjects[i]->isDirectional()) {
148 <            ji = integrableObjects[i]->getJ();
149 <
150 <            for(j = 0; j < 3; j++)
151 <        oldJi[3 * i + j] = ji[j];
144 >    index = 0;
145 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
146 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
147 >               integrableObject = mol->nextIntegrableObject(j)) {
148 >                oldVel_[index] = integrableObject->getVel();
149 >                oldJi_[index] = integrableObject->getJ();                
150          }
151 +        ++index;              
152      }
153  
154      // do the iteration:
155  
156 <    for(k = 0; k < 4; k++) {
156 >    for(int k = 0; k < maxIterNum_; k++) {
157 >        index = 0;
158          instTemp = tStats->getTemperature();
159  
160          // evolve chi another half step using the temperature at t + dt/2
161  
162          prevChi = chi;
163 <        chi = oldChi + dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
163 >        chi = oldChi + dt2 * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
164  
165 <        for(i = 0; i < integrableObjects.size(); i++) {
166 <            integrableObjects[i]->getFrc(frc);
167 <            vel = integrableObjects[i]->getVel();
165 >        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
166 >            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
167 >                   integrableObject = mol->nextIntegrableObject(j)) {
168  
169 <            mass = integrableObjects[i]->getMass();
169 >                frc = integrableObject->getFrc();
170 >                vel = integrableObject->getVel();
171  
172 <            // velocity half step
172 <            for(j = 0; j < 3; j++)
173 <        vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel[3*i + j]*chi);
172 >                mass = integrableObject->getMass();
173  
174 <            integrableObjects[i]->setVel(vel);
174 >                // velocity half step
175 >                //for(j = 0; j < 3; j++)
176 >                //    vel[j] = oldVel_[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
177 >                vel = oldVel_ + dt2/mass*eConvert * frc - dt2*chi*oldVel_[index];
178 >            
179 >                integrableObject->setVel(vel);
180  
181 <            if (integrableObjects[i]->isDirectional()) {
181 >                if (integrableObject->isDirectional()) {
182  
183 <                // get and convert the torque to body frame
183 >                    // get and convert the torque to body frame
184  
185 <                Tb = integrableObjects[i]->getTrq();
186 <                integrableObjects[i]->lab2Body(Tb);
185 >                    Tb = integrableObject->getTrq();
186 >                    integrableObject->lab2Body(Tb);
187  
188 <                for(j = 0; j < 3; j++)
189 <            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
188 >                    //for(j = 0; j < 3; j++)
189 >                    //    ji[j] = oldJi_[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
190 >                    ji += dt2*eConvert*Tb - dt2*ch *oldJi_[index];
191  
192 <                integrableObjects[i]->setJ(ji);
192 >                    integrableObject->setJ(ji);
193 >                }
194              }
195          }
196 +    
197  
198 <        if (nConstrained)
192 <            constrainB();
198 >        //constraintAlgorithm->doConstrainB();
199  
200          if (fabs(prevChi - chi) <= chiTolerance)
201              break;
202 +
203 +        ++index;
204      }
205  
206      integralOfChidt += dt2 * chi;
199 – }
207  
208 < template <typename T>
209 < void NVT<T>::resetIntegrator(void) {
203 <    chi = 0.0;
204 <    integralOfChidt = 0.0;
208 >    currSnapshot->setChi(chi);
209 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);
210   }
211  
207 – template <typename T>
208 – int NVT<T>::readyCheck() {
212  
213 <    //check parent's readyCheck() first
211 <    if (T::readyCheck() == -1)
212 <        return -1;
213 <
214 <    // First check to see if we have a target temperature.
215 <    // Not having one is fatal.
216 <
217 <    if (!have_target_temp) {
218 <        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
219 <        painCave.isFatal = 1;
220 <        painCave.severity = OOPSE_ERROR;
221 <        simError();
222 <        return -1;
223 <    }
224 <
225 <    // We must set tauThermostat.
226 <
227 <    if (!have_tau_thermostat) {
228 <        sprintf(painCave.errMsg, "If you use the constant temperature\n"
229 <                                     "\tintegrator, you must set tauThermostat.\n");
230 <
231 <        painCave.severity = OOPSE_ERROR;
232 <        painCave.isFatal = 1;
233 <        simError();
234 <        return -1;
235 <    }
236 <
237 <    if (!have_chi_tolerance) {
238 <        sprintf(painCave.errMsg, "In NVT integrator: setting chi tolerance to 1e-6\n");
239 <        chiTolerance = 1e - 6;
240 <        have_chi_tolerance = 1;
241 <        painCave.severity = OOPSE_INFO;
242 <        painCave.isFatal = 0;
243 <        simError();
244 <    }
245 <
246 <    return 1;
247 < }
248 <
249 < template <typename T>
250 < double NVT<T>::getConservedQuantity(void) {
213 > double NVT::calcConservedQuantity() {
214      double conservedQuantity;
215      double fkBT;
216      double Energy;
217      double thermostat_kinetic;
218      double thermostat_potential;
219  
220 <    fkBT = (double)(info->ndf) *kB *targetTemp;
220 >    fkBT = info_->getNdf() *kB *targetTemp_;
221  
222      Energy = tStats->getTotalE();
223  
224 <    thermostat_kinetic = fkBT * tauThermostat * tauThermostat * chi * chi / (2.0 * eConvert);
224 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
225  
226      thermostat_potential = fkBT * integralOfChidt / eConvert;
227  
# Line 267 | Line 230 | template <typename T>
230      return conservedQuantity;
231   }
232  
270 – template <typename T>
271 – string NVT<T>::getAdditionalParameters(void) {
272 –    string parameters;
273 –    const int BUFFERSIZE = 2000; // size of the read buffer
274 –    char buffer[BUFFERSIZE];
233  
234 <    sprintf(buffer, "\t%G\t%G;", chi, integralOfChidt);
277 <    parameters += buffer;
278 <
279 <    return parameters;
280 < }
234 > }//end namespace oopse

Diff Legend

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