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

Comparing branches/new_design/OOPSE-2.0/src/integrators/NVT.cpp (file contents):
Revision 1764 by tim, Fri Nov 19 21:38:22 2004 UTC vs.
Revision 1765 by tim, Mon Nov 22 20:55:52 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
14 <
15 < template <typename T>NVT<T>::NVT(SimInfo *theInfo, ForceFields *the_ff) :
16 <    T(theInfo, the_ff) {
5 > NVT::NVT(SimInfo* Info) : VelocityVerletIntegrator(info){
6      GenericData *data;
7      DoubleGenericData *chiValue;
8      DoubleGenericData *integralOfChidtValue;
# Line 49 | Line 38 | template <typename T>NVT<T>::NVT(SimInfo *theInfo, For
38          }
39      }
40  
41 <    oldVel = new double[3 * integrableObjects.size()];
42 <    oldJi = new double[3 * integrableObjects.size()];
41 >    oldVel_.resize(info_->getNIntegrableObjects());
42 >    oldJi_.resize(info_->getNIntegrableObjects());
43   }
44  
45 < template <typename T>NVT<T>::~NVT() {
57 <    delete [] oldVel;
58 <    delete [] oldJi;
45 > NVT::~NVT() {
46   }
47  
48 < template <typename T>
49 < void NVT<T>::moveA() {
50 <    int i, j;
51 <    DirectionalAtom *dAtom;
48 > void NVT::moveA() {
49 >    typename SimInfo::MoleculeIterator i;
50 >    typename Molecule::IntegrableObjectIterator  j;
51 >    Molecule* mol;
52 >    StuntDouble* integrableObject;
53      Vector3d Tb;
54      Vector3d ji;
55      double mass;
# Line 75 | Line 63 | void NVT<T>::moveA() {
63  
64      instTemp = tStats->getTemperature();
65  
66 <    for(i = 0; i < integrableObjects.size(); i++) {
67 <        vel = integrableObjects[i]->getVel();
68 <        pos = integrableObjects[i]->getPos();
81 <        integrableObjects[i]->getFrc(frc);
66 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
67 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
68 >               integrableObject = mol->nextIntegrableObject(j)) {
69  
70 <        mass = integrableObjects[i]->getMass();
70 >        vel = integrableObject->getVel();
71 >        pos = integrableObject->getPos();
72 >        frc = integrableObject->getFrc();
73  
74 <        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);
74 >        mass = integrableObject->getMass();
75  
76 <            // position whole step
77 <            pos[j] += dt * vel[j];
78 <        }
76 >        // velocity half step  (use chi from previous step here):
77 >        //vel[j] += halfStep * ((frc[j] / mass ) * eConvert - vel[j]*chi);
78 >        vel = halfStep *eConvert/mass*frc - halfStep*chi*vel;
79 >        
80 >        // position whole step
81 >        //pos[j] += dt * vel[j];
82 >        pos += fullStep * vel;
83  
84 <        integrableObjects[i]->setVel(vel);
85 <        integrableObjects[i]->setPos(pos);
84 >        integrableObject->setVel(vel);
85 >        integrableObject->setPos(pos);
86  
87 <        if (integrableObjects[i]->isDirectional()) {
87 >        if (integrableObject->isDirectional()) {
88  
89              // get and convert the torque to body frame
90  
91 <            Tb = integrableObjects[i]->getTrq();
92 <            integrableObjects[i]->lab2Body(Tb);
91 >            Tb = integrableObject->getTrq();
92 >            integrableObject->lab2Body(Tb);
93  
94              // get the angular momentum, and propagate a half step
95  
96 <            ji = integrableObjects[i]->getJ();
96 >            ji = integrableObject->getJ();
97  
98 <            for(j = 0; j < 3; j++)
99 <        ji[j] += dt2 * (Tb[j] * eConvert - ji[j]*chi);
98 >            //ji[j] += halfStep * (Tb[j] * eConvert - ji[j]*chi);
99 >            ji += halfStep*eConvert*Tb - halfStep*ch *ji;
100 >            this->rotationPropagation(integrableObject, ji);
101  
102 <            this->rotationPropagation(integrableObjects[i], ji);
111 <
112 <            integrableObjects[i]->setJ(ji);
102 >            integrableObject->setJ(ji);
103          }
104      }
105  
106 +    }
107 +    
108      if (nConstrained)
109          constrainA();
110  
111      // Finally, evolve chi a half step (just like a velocity) using
112      // temperature at time t, not time t+dt/2
113  
114 <    //std::cerr << "targetTemp = " << targetTemp << " instTemp = " << instTemp << " tauThermostat = " << tauThermostat << " integral of Chi = " << integralOfChidt << "\n";
114 >    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
115 >    double chi = currSnapshot->getChi();
116 >    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
117 >    
118 >    chi += halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
119 >    integralOfChidt += chi * halfStep;
120  
121 <    chi += dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
122 <    integralOfChidt += chi * dt2;
121 >    currSnapshot->setChi(chi);
122 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);
123   }
124  
125 < template <typename T>
126 < void NVT<T>::moveB(void) {
127 <    int i, j, k;
128 <    double Tb[3], ji[3];
129 <    double vel[3], frc[3];
125 > void NVT::moveB() {
126 >    typename SimInfo::MoleculeIterator i;
127 >    typename Molecule::IntegrableObjectIterator  j;
128 >    Molecule* mol;
129 >    StuntDouble* integrableObject;
130 >    
131 >    Vector3d Tb;
132 >    Vector3d ji;    
133 >    Vector3d vel;
134 >    Vector3d frc;
135      double mass;
136      double instTemp;
137      double oldChi, prevChi;
138 <
138 >    int index;
139      // Set things up for the iteration:
140  
141 +    Snapshot* currSnapshot = info_->getSnapshotManager()->getCurrentSnapshot();
142 +    double chi = currSnapshot->getChi();
143 +    double integralOfChidt = currSnapshot->getIntegralOfChiDt();
144      oldChi = chi;
145  
146 <    for(i = 0; i < integrableObjects.size(); i++) {
147 <        vel = integrableObjects[i]->getVel();
148 <
149 <        for(j = 0; j < 3; j++)
150 <    oldVel[3 * i + j] = vel[j];
151 <
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];
146 >    index = 0;
147 >    for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
148 >        for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
149 >               integrableObject = mol->nextIntegrableObject(j)) {
150 >                oldVel_[index] = integrableObject->getVel();
151 >                oldJi_[index] = integrableObject->getJ();                
152          }
153 +        ++index;              
154      }
155  
156      // do the iteration:
157  
158 <    for(k = 0; k < 4; k++) {
158 >    for(int k = 0; k < maxIterNum_; k++) {
159 >        index = 0;
160          instTemp = tStats->getTemperature();
161  
162          // evolve chi another half step using the temperature at t + dt/2
163  
164          prevChi = chi;
165 <        chi = oldChi + dt2 * (instTemp / targetTemp - 1.0) / (tauThermostat * tauThermostat);
165 >        chi = oldChi + halfStep * (instTemp / targetTemp_ - 1.0) / (tauThermostat_ * tauThermostat_);
166  
167 <        for(i = 0; i < integrableObjects.size(); i++) {
168 <            integrableObjects[i]->getFrc(frc);
169 <            vel = integrableObjects[i]->getVel();
167 >        for (mol = info_->beginMolecule(i); mol != NULL; mol = info_->nextMolecule(i)) {
168 >            for (integrableObject = mol->beginIntegrableObject(j); integrableObject != NULL;
169 >                   integrableObject = mol->nextIntegrableObject(j)) {
170  
171 <            mass = integrableObjects[i]->getMass();
171 >                frc = integrableObject->getFrc();
172 >                vel = integrableObject->getVel();
173  
174 <            // 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);
174 >                mass = integrableObject->getMass();
175  
176 <            integrableObjects[i]->setVel(vel);
176 >                // velocity half step
177 >                //for(j = 0; j < 3; j++)
178 >                //    vel[j] = oldVel_[3*i+j] + halfStep * ((frc[j] / mass ) * eConvert - oldVel_[3*i + j]*chi);
179 >                vel = oldVel_ + halfStep/mass*eConvert * frc - halfStep*chi*oldVel_[index];
180 >            
181 >                integrableObject->setVel(vel);
182  
183 <            if (integrableObjects[i]->isDirectional()) {
183 >                if (integrableObject->isDirectional()) {
184  
185 <                // get and convert the torque to body frame
185 >                    // get and convert the torque to body frame
186  
187 <                Tb = integrableObjects[i]->getTrq();
188 <                integrableObjects[i]->lab2Body(Tb);
187 >                    Tb = integrableObject->getTrq();
188 >                    integrableObject->lab2Body(Tb);
189  
190 <                for(j = 0; j < 3; j++)
191 <            ji[j] = oldJi[3*i + j] + dt2 * (Tb[j] * eConvert - oldJi[3*i+j]*chi);
190 >                    //for(j = 0; j < 3; j++)
191 >                    //    ji[j] = oldJi_[3*i + j] + halfStep * (Tb[j] * eConvert - oldJi_[3*i+j]*chi);
192 >                    ji += halfStep*eConvert*Tb - halfStep*ch *oldJi_[index];
193  
194 <                integrableObjects[i]->setJ(ji);
194 >                    integrableObject->setJ(ji);
195 >                }
196              }
197          }
198 +    
199  
200          if (nConstrained)
201              constrainB();
202  
203          if (fabs(prevChi - chi) <= chiTolerance)
204              break;
205 +
206 +        ++index;
207      }
208  
209 <    integralOfChidt += dt2 * chi;
199 < }
209 >    integralOfChidt += halfStep * chi;
210  
211 < template <typename T>
212 < void NVT<T>::resetIntegrator(void) {
203 <    chi = 0.0;
204 <    integralOfChidt = 0.0;
211 >    currSnapshot->setChi(chi);
212 >    currSnapshot->setIntegralOfChiDt(integralOfChidt);
213   }
214  
215   template <typename T>
# Line 215 | Line 223 | int NVT<T>::readyCheck() {
223      // Not having one is fatal.
224  
225      if (!have_target_temp) {
226 <        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp!\n");
226 >        sprintf(painCave.errMsg, "You can't use the NVT integrator without a targetTemp_!\n");
227          painCave.isFatal = 1;
228          painCave.severity = OOPSE_ERROR;
229          simError();
230          return -1;
231      }
232  
233 <    // We must set tauThermostat.
233 >    // We must set tauThermostat_.
234  
235      if (!have_tau_thermostat) {
236          sprintf(painCave.errMsg, "If you use the constant temperature\n"
237 <                                     "\tintegrator, you must set tauThermostat.\n");
237 >                                     "\tintegrator, you must set tauThermostat_.\n");
238  
239          painCave.severity = OOPSE_ERROR;
240          painCave.isFatal = 1;
# Line 254 | Line 262 | double NVT<T>::getConservedQuantity(void) {
262      double thermostat_kinetic;
263      double thermostat_potential;
264  
265 <    fkBT = (double)(info->ndf) *kB *targetTemp;
265 >    fkBT = info_->getNdf() *kB *targetTemp_;
266  
267      Energy = tStats->getTotalE();
268  
269 <    thermostat_kinetic = fkBT * tauThermostat * tauThermostat * chi * chi / (2.0 * eConvert);
269 >    thermostat_kinetic = fkBT * tauThermostat_ * tauThermostat_ * chi * chi / (2.0 * eConvert);
270  
271      thermostat_potential = fkBT * integralOfChidt / eConvert;
272  
# Line 267 | Line 275 | template <typename T>
275      return conservedQuantity;
276   }
277  
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];
278  
279 <    sprintf(buffer, "\t%G\t%G;", chi, integralOfChidt);
277 <    parameters += buffer;
278 <
279 <    return parameters;
280 < }
279 > }//end namespace oopse

Diff Legend

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