ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/branches/new-templateless/OOPSE/libmdtools/NPT.cpp
(Generate patch)

Comparing branches/new-templateless/OOPSE/libmdtools/NPT.cpp (file contents):
Revision 848, Fri Oct 31 18:28:53 2003 UTC vs.
Revision 851 by mmeineke, Wed Nov 5 19:18:17 2003 UTC

# Line 1 | Line 1
1 + #include <stdlib.h>
2   #include <math.h>
3 +
4   #include "Atom.hpp"
5   #include "SRI.hpp"
6   #include "AbstractClasses.hpp"
# Line 24 | Line 26
26   //
27   //    Hoover, W. G., 1986, Phys. Rev. A, 34, 2499.
28  
29 < template<typename T> NPT<T>::NPT ( SimInfo *theInfo, ForceFields* the_ff):
30 <  T( theInfo, the_ff )
29 > NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff):
30 >  Integrator( theInfo, the_ff )
31   {
32    GenericData* data;
33 <  DoubleData * chiValue;
32 <  DoubleData * integralOfChidtValue;
33 <
34 <  chiValue = NULL;
35 <  integralOfChidtValue = NULL;
36 <
33 >  
34    chi = 0.0;
35    integralOfChidt = 0.0;
36    have_tau_thermostat = 0;
# Line 46 | Line 43 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
43  
44    // retrieve chi and integralOfChidt from simInfo
45    data = info->getProperty(CHIVALUE_ID);
46 <  if(data){
47 <    chiValue = dynamic_cast<DoubleData*>(data);
46 >  if(data != NULL ){
47 >    chi = data->getDval();
48    }
49  
50    data = info->getProperty(INTEGRALOFCHIDT_ID);
51 <  if(data){
52 <    integralOfChidtValue = dynamic_cast<DoubleData*>(data);
51 >  if(data != NULL ){
52 >    integralOfChidt = data->getDval();
53    }
54  
58 –  // chi and integralOfChidt should appear by pair
59 –  if(chiValue && integralOfChidtValue){
60 –    chi = chiValue->getData();
61 –    integralOfChidt = integralOfChidtValue->getData();
62 –  }
63 –
55    oldPos = new double[3*nAtoms];
56    oldVel = new double[3*nAtoms];
57    oldJi = new double[3*nAtoms];
# Line 72 | Line 63 | template<typename T> NPT<T>::NPT ( SimInfo *theInfo, F
63  
64   }
65  
66 < template<typename T> NPT<T>::~NPT() {
66 > NPT::~NPT() {
67    delete[] oldPos;
68    delete[] oldVel;
69    delete[] oldJi;
70   }
71  
72 < template<typename T> void NPT<T>::moveA() {
72 > void NPT::moveA() {
73  
74    //new version of NPT
75    int i, j, k;
# Line 180 | Line 171 | template<typename T> void NPT<T>::moveA() {
171    this->scaleSimBox();
172   }
173  
174 < template<typename T> void NPT<T>::moveB( void ){
174 > void NPT::moveB( void ){
175  
176    //new version of NPT
177    int i, j, k;
178    DirectionalAtom* dAtom;
179    double Tb[3], ji[3], sc[3];
180 <  double vel[3], frc[3];
180 >  double vel[3], frc[3], qVel[3];
181    double mass;
182  
183    // Set things up for the iteration:
# Line 216 | Line 207 | template<typename T> void NPT<T>::moveB( void ){
207  
208    for (k=0; k < 4; k++) {
209  
210 +    atoms[0]->getVel(vel);
211 +
212 + //    std::cerr << "vel atom0 = "
213 + //            << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";
214 +    
215      instaTemp = tStats->getTemperature();
216 +    
217 +    std::cerr << "instaTemp MoveB = " << instaTemp << "\n";
218 +    
219      instaPress = tStats->getPressure();
220  
221      // evolve chi another half step using the temperature at t + dt/2
# Line 233 | Line 232 | template<typename T> void NPT<T>::moveB( void ){
232  
233        getVelScaleB( sc, i );
234  
235 +      for(j=0;j<3;j++)
236 +        qVel[j] = vel[j];
237 +
238        // velocity half step
239        for (j=0; j < 3; j++)
240          vel[j] = oldVel[3*i+j] + dt2 * ((frc[j] / mass ) * eConvert - sc[j]);
241  
242 + //       std::cerr << " vel of " << i << " => "
243 + //              << vel[0] << ", " << vel[1] << ", " << vel[2] << "\n";
244 +
245 +      double delVel = 0.0;
246 +      for(j=0;j<3;j++)
247 +        delVel += pow( (qVel[j] - vel[j]), 2 );
248 +      
249 +      delVel = sqrt(delVel);
250 +      
251 +      std::cerr << "delVel[" << i << "] = " << delVel << "\n";
252 +
253        atoms[i]->setVel( vel );
254  
255        if( atoms[i]->isDirectional() ){
# Line 268 | Line 281 | template<typename T> void NPT<T>::moveB( void ){
281  
282   }
283  
284 < template<typename T> void NPT<T>::resetIntegrator() {
284 > void NPT::resetIntegrator() {
285    chi = 0.0;
286 <  T::resetIntegrator();
286 >  Integrator::resetIntegrator();
287   }
288  
289 < template<typename T> void NPT<T>::evolveChiA() {
289 > void NPT::evolveChiA() {
290 >
291 >  std::cerr << "\n"
292 >            << "evolveChiA before:\n"
293 >            << "chi = " << chi << "\n"
294 >            << "oldChi = " << oldChi << "\n"
295 >            << "\n"
296 >            << "instaTemp = " << instaTemp << "\n"
297 >            << "targetTemp = " << targetTemp << "\n"
298 >            << "dt2 = " << dt2 << "\n"
299 >            << "tt2 = " << tt2 << "\n";
300 >
301    chi += dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
302    oldChi = chi;
303 +
304 +  std::cerr << "\n"
305 +            << "evolveChiA after:\n"
306 +            << "chi = " << chi << "\n"
307 +            << "oldChi = " << oldChi << "\n"
308 +            << "\n"
309 +            << "instaTemp = " << instaTemp << "\n"
310 +            << "targetTemp = " << targetTemp << "\n"
311 +            << "dt2 = " << dt2 << "\n"
312 +            << "tt2 = " << tt2 << "\n";
313 +
314   }
315  
316 < template<typename T> void NPT<T>::evolveChiB() {
316 > void NPT::evolveChiB() {
317  
318 +  std::cerr << "\n"
319 +            << "evolveChiB before:\n"
320 +            << "chi = " << chi << "\n"
321 +            << "oldChi = " << oldChi << "\n"
322 +            << "\n"
323 +            << "instaTemp = " << instaTemp << "\n"
324 +            << "targetTemp = " << targetTemp << "\n"
325 +            << "dt2 = " << dt2 << "\n"
326 +            << "tt2 = " << tt2 << "\n";
327 +
328    prevChi = chi;
329    chi = oldChi + dt2 * ( instaTemp / targetTemp - 1.0) / tt2;
330 +
331 +  std::cerr << "\n"
332 +            << "evolveChiB after:\n"
333 +            << "chi = " << chi << "\n"
334 +            << "oldChi = " << oldChi << "\n"
335 +            << "\n"
336 +            << "instaTemp = " << instaTemp << "\n"
337 +            << "targetTemp = " << targetTemp << "\n"
338 +            << "dt2 = " << dt2 << "\n"
339 +            << "tt2 = " << tt2 << "\n";
340   }
341  
342 < template<typename T> bool NPT<T>::chiConverged() {
342 > bool NPT::chiConverged() {
343  
344    return ( fabs( prevChi - chi ) <= chiTolerance );
345   }
346  
347 < template<typename T> int NPT<T>::readyCheck() {
347 > int NPT::readyCheck() {
348  
349    //check parent's readyCheck() first
350 <  if (T::readyCheck() == -1)
350 >  if (Integrator::readyCheck() == -1)
351      return -1;
352  
353    // First check to see if we have a target temperature.

Diff Legend

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