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 849 by mmeineke, Fri Oct 31 21:06:47 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 28 | Line 30 | NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff):
30    Integrator( theInfo, the_ff )
31   {
32    GenericData* data;
31 –  DoubleData * chiValue;
32 –  DoubleData * integralOfChidtValue;
33    
34 –  chiValue = NULL;
35 –  integralOfChidtValue = NULL;
36 –  
34    chi = 0.0;
35    integralOfChidt = 0.0;
36    have_tau_thermostat = 0;
# Line 46 | Line 43 | NPT::NPT ( SimInfo *theInfo, ForceFields* the_ff):
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 186 | Line 177 | void NPT::moveB( void ){
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 | void NPT::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 | void NPT::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 274 | Line 287 | void NPT::evolveChiA() {
287   }
288  
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   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   bool NPT::chiConverged() {

Diff Legend

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