ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE/libmdtools/Symplectic.cpp
(Generate patch)

Comparing trunk/OOPSE/libmdtools/Symplectic.cpp (file contents):
Revision 423 by mmeineke, Thu Mar 27 20:12:15 2003 UTC vs.
Revision 477 by gezelter, Tue Apr 8 14:34:30 2003 UTC

# Line 5 | Line 5
5   #include "Thermo.hpp"
6   #include "ReadWrite.hpp"
7   #include "ForceFields.hpp"
8 + #include "ExtendedSystem.hpp"
9   #include "simError.h"
10  
11   extern "C"{
# Line 30 | Line 31 | extern "C"{
31  
32  
33  
34 < Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff ){
34 > Symplectic::Symplectic( SimInfo* the_entry_plug, ForceFields* the_ff,
35 >                        ExtendedSystem* the_es ){
36    entry_plug = the_entry_plug;
37    myFF = the_ff;
38 +  myES = the_es;
39    isFirst = 1;
40  
41 +  std::cerr<< "calling symplectic constructor\n";
42 +
43    molecules = entry_plug->molecules;
44    nMols = entry_plug->n_mol;
45  
# Line 48 | Line 53 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
53    mass = new double[entry_plug->n_atoms];
54    for(int i = 0; i < entry_plug->n_atoms; i++){
55      mass[i] = entry_plug->atoms[i]->getMass();
56 <  }
52 <
53 <  
56 >  }  
57  
58    // check for constraints
59  
# Line 65 | Line 68 | Symplectic::Symplectic( SimInfo* the_entry_plug, Force
68    SRI** theArray;
69    for(int i = 0; i < nMols; i++){
70      
71 <    theArray = molecules[i].getMyBonds();
72 <    for(int j=0; j<molecules[i].getNbonds(); j++){
71 >    theArray = (SRI**) molecules[i].getMyBonds();
72 >    for(int j=0; j<molecules[i].getNBonds(); j++){
73        
74        constrained = theArray[j]->is_constrained();
75        
76        if(constrained){
77          
78          dummy_plug = theArray[j]->get_constraint();
79 <        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
80 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
81 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
79 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
80 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
81 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
82          
83 <        c_n_constrained++;
83 >        n_constrained++;
84          constrained = 0;
85        }
86      }
87  
88 <    theArray = molecules[i].getMyBends();
89 <    for(int j=0; j<molecules[i].getNbends(); j++){
88 >    theArray = (SRI**) molecules[i].getMyBends();
89 >    for(int j=0; j<molecules[i].getNBends(); j++){
90        
91        constrained = theArray[j]->is_constrained();
92        
93        if(constrained){
94          
95          dummy_plug = theArray[j]->get_constraint();
96 <        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
97 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
98 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
96 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
97 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
98 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
99          
100 <        c_n_constrained++;
100 >        n_constrained++;
101          constrained = 0;
102        }
103      }
104  
105 <    theArray = molecules[i].getMyTorsions();
106 <    for(int j=0; j<molecules[i].getNtorsions(); j++){
105 >    theArray = (SRI**) molecules[i].getMyTorsions();
106 >    for(int j=0; j<molecules[i].getNTorsions(); j++){
107        
108        constrained = theArray[j]->is_constrained();
109        
110        if(constrained){
111          
112          dummy_plug = theArray[j]->get_constraint();
113 <        temp_con[c_n_constrained].set_a( dummy_plug->get_a() );
114 <        temp_con[c_n_constrained].set_b( dummy_plug->get_b() );
115 <        temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() );
113 >        temp_con[n_constrained].set_a( dummy_plug->get_a() );
114 >        temp_con[n_constrained].set_b( dummy_plug->get_b() );
115 >        temp_con[n_constrained].set_dsqr( dummy_plug->get_dsqr() );
116          
117 <        c_n_constrained++;
117 >        n_constrained++;
118          constrained = 0;
119        }
120      }
# Line 183 | Line 186 | void Symplectic::integrate( void ){
186    int status_n = (int)( statusTime / dt );
187    int vel_n    = (int)( thermalTime / dt );
188  
189 <  int calcPot;
189 >  int calcPot, calcStress;
190  
191 <   Thermo *tStats = new Thermo( entry_plug );
191 >  Thermo *tStats;
192 >  StatWriter*  e_out;
193 >  DumpWriter*  dump_out;
194  
195 <  StatWriter*  e_out    = new StatWriter( entry_plug );
191 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
195 >  std::cerr << "about to call new thermo\n";
196  
197 +  tStats   = new Thermo( entry_plug );
198 +  e_out    = new StatWriter( entry_plug );
199 +
200 +  std::cerr << "calling dumpWriter \n";
201 +  dump_out = new DumpWriter( entry_plug );
202 +  std::cerr << "called dumpWriter \n";
203 +
204    Atom** atoms = entry_plug->atoms;
205    DirectionalAtom* dAtom;
206    dt2 = 0.5 * dt;
207  
208    // initialize the forces the before the first step
209  
210 <  myFF->doForces(1,0);
210 >  myFF->doForces(1,1);
211    
212    if( entry_plug->setTemp ){
213      
# Line 208 | Line 219 | void Symplectic::integrate( void ){
219    
220    calcPot = 0;
221  
222 +  if (!strcasecmp( entry_plug->ensemble, "NPT")) {
223 +    calcStress = 1;
224 +  } else {
225 +    calcStress = 0;
226 +  }
227 +
228    if( n_constrained ){
229  
230      double *Rx = new double[nAtoms];
# Line 224 | Line 241 | void Symplectic::integrate( void ){
241      
242  
243      for( tl=0; tl < n_loops; tl++ ){
244 +
245 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
246 +        myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() );
247        
248        for( j=0; j<nAtoms; j++ ){
249  
# Line 261 | Line 281 | void Symplectic::integrate( void ){
281        }
282  
283  
284 <      for( i=0; i<nAtoms; i++ ){
285 <        if( atoms[i]->isDirectional() ){
284 > //       for( i=0; i<nAtoms; i++ ){
285 > // //   if( atoms[i]->isDirectional() ){
286                    
287 <          dAtom = (DirectionalAtom *)atoms[i];
287 > // //     dAtom = (DirectionalAtom *)atoms[i];
288            
289 <          // get and convert the torque to body frame
289 > // //     // get and convert the torque to body frame
290            
291 <          Tb[0] = dAtom->getTx();
292 <          Tb[1] = dAtom->getTy();
293 <          Tb[2] = dAtom->getTz();
291 > // //     Tb[0] = dAtom->getTx();
292 > // //     Tb[1] = dAtom->getTy();
293 > // //     Tb[2] = dAtom->getTz();
294            
295 <          dAtom->lab2Body( Tb );
295 > // //     dAtom->lab2Body( Tb );
296            
297 <          // get the angular momentum, and propagate a half step
297 > // //     // get the angular momentum, and propagate a half step
298            
299 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
299 > // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
300 > // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
301 > // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
302            
303 <          // get the atom's rotation matrix
303 > // //     // get the atom's rotation matrix
304            
305 <          A[0][0] = dAtom->getAxx();
306 <          A[0][1] = dAtom->getAxy();
307 <          A[0][2] = dAtom->getAxz();
305 > // //     A[0][0] = dAtom->getAxx();
306 > // //     A[0][1] = dAtom->getAxy();
307 > // //     A[0][2] = dAtom->getAxz();
308            
309 <          A[1][0] = dAtom->getAyx();
310 <          A[1][1] = dAtom->getAyy();
311 <          A[1][2] = dAtom->getAyz();
309 > // //     A[1][0] = dAtom->getAyx();
310 > // //     A[1][1] = dAtom->getAyy();
311 > // //     A[1][2] = dAtom->getAyz();
312            
313 <          A[2][0] = dAtom->getAzx();
314 <          A[2][1] = dAtom->getAzy();
315 <          A[2][2] = dAtom->getAzz();
296 <          
297 <          
298 <          // use the angular velocities to propagate the rotation matrix a
299 <          // full time step
313 > // //     A[2][0] = dAtom->getAzx();
314 > // //     A[2][1] = dAtom->getAzy();
315 > // //     A[2][2] = dAtom->getAzz();
316            
317            
318 <          angle = dt2 * ji[0] / dAtom->getIxx();
319 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
318 > // //     // use the angular velocities to propagate the rotation matrix a
319 > // //     // full time step
320            
305 –          angle = dt2 * ji[1] / dAtom->getIyy();
306 –          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
321            
322 <          angle = dt * ji[2] / dAtom->getIzz();
323 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
322 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
323 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
324            
325 <          angle = dt2 * ji[1] / dAtom->getIyy();
326 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
325 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
326 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
327            
328 <          angle = dt2 * ji[0] / dAtom->getIxx();
329 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
328 > // //     angle = dt * ji[2] / dAtom->getIzz();
329 > // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
330            
331 + // //     angle = dt2 * ji[1] / dAtom->getIyy();
332 + // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
333            
334 <          dAtom->setA( A );
335 <          dAtom->setJx( ji[0] );
336 <          dAtom->setJy( ji[1] );
337 <          dAtom->setJz( ji[2] );
338 <        }
339 <      }
334 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
335 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
336 >          
337 >          
338 > // //     dAtom->setA( A );
339 > // //     dAtom->setJx( ji[0] );
340 > // //     dAtom->setJy( ji[1] );
341 > // //     dAtom->setJz( ji[2] );
342 > // //   }
343 > //       }
344        
345        // calculate the forces
346        
347 <      myFF->doForces(calcPot, 0);
347 >      myFF->doForces(calcPot, calcStress);
348        
349        // move b
350  
# Line 362 | Line 382 | void Symplectic::integrate( void ){
382          atoms[j]->set_vz(Vz[j]);
383        }
384        
385 <      for( i=0; i< nAtoms; i++ ){
385 > //       for( i=0; i< nAtoms; i++ ){
386  
387 <        if( atoms[i]->isDirectional() ){
387 > //      if( atoms[i]->isDirectional() ){
388  
389 <          dAtom = (DirectionalAtom *)atoms[i];
389 > //        dAtom = (DirectionalAtom *)atoms[i];
390            
391 <          // get and convert the torque to body frame
391 > //        // get and convert the torque to body frame
392            
393 <          Tb[0] = dAtom->getTx();
394 <          Tb[1] = dAtom->getTy();
395 <          Tb[2] = dAtom->getTz();
393 > //        Tb[0] = dAtom->getTx();
394 > //        Tb[1] = dAtom->getTy();
395 > //        Tb[2] = dAtom->getTz();
396            
397 <          dAtom->lab2Body( Tb );
397 > //        dAtom->lab2Body( Tb );
398            
399 <          // get the angular momentum, and complete the angular momentum
400 <          // half step
399 > //        // get the angular momentum, and complete the angular momentum
400 > //        // half step
401            
402 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
402 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
403 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
404 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
405            
406 <          dAtom->setJx( ji[0] );
407 <          dAtom->setJy( ji[1] );
408 <          dAtom->setJz( ji[2] );
409 <        }
410 <      }
406 > //        dAtom->setJx( ji[0] );
407 > //        dAtom->setJy( ji[1] );
408 > //        dAtom->setJz( ji[2] );
409 > //      }
410 > //       }
411      
412 +
413 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
414 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
415 +
416 +      if (!strcasecmp( entry_plug->ensemble, "NPT") )
417 +        myES->NoseHooverAndersonNPT( dt,
418 +                                     tStats->getKinetic(),
419 +                                     tStats->getPressure());
420 +
421        time = tl + 1;
422        
423        if( entry_plug->setTemp ){
424          if( !(time % vel_n) ) tStats->velocitize();
425        }
426        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
427 <      if( !((time+1) % status_n) ) calcPot = 1;
428 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
427 >      if( !((time+1) % status_n) ) {
428 >        calcPot = 1;
429 >        // bitwise masking in case we need it for NPT
430 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
431 >      }
432 >      if( !(time % status_n) ){
433 >        e_out->writeStat( time * dt );
434 >        calcPot = 0;
435 >        // bitwise masking in case we need it for NPT
436 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
437 >      }
438      }
439    }
440    else{
# Line 406 | Line 444 | void Symplectic::integrate( void ){
444        kE = 0.0;
445        rot_kE= 0.0;
446        trans_kE = 0.0;
447 +
448 +      if (!strcasecmp( entry_plug->ensemble, "NVT"))
449 +        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
450        
451        for( i=0; i<nAtoms; i++ ){
452          
# Line 432 | Line 473 | void Symplectic::integrate( void ){
473          atoms[i]->set_vy( vy );
474          atoms[i]->set_vz( vz );
475          
476 <        if( atoms[i]->isDirectional() ){
476 > //      if( atoms[i]->isDirectional() ){
477  
478 <          dAtom = (DirectionalAtom *)atoms[i];
478 > //        dAtom = (DirectionalAtom *)atoms[i];
479            
480 <          // get and convert the torque to body frame
480 > //        // get and convert the torque to body frame
481            
482 <          Tb[0] = dAtom->getTx();
483 <          Tb[1] = dAtom->getTy();
484 <          Tb[2] = dAtom->getTz();
482 > //        Tb[0] = dAtom->getTx();
483 > //        Tb[1] = dAtom->getTy();
484 > //        Tb[2] = dAtom->getTz();
485            
486 <          dAtom->lab2Body( Tb );
486 > //        dAtom->lab2Body( Tb );
487            
488 <          // get the angular momentum, and propagate a half step
488 > //        // get the angular momentum, and propagate a half step
489            
490 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
490 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
491 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
492 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
493            
494 <          // get the atom's rotation matrix
494 > //        // get the atom's rotation matrix
495            
496 <          A[0][0] = dAtom->getAxx();
497 <          A[0][1] = dAtom->getAxy();
498 <          A[0][2] = dAtom->getAxz();
496 > //        A[0][0] = dAtom->getAxx();
497 > //        A[0][1] = dAtom->getAxy();
498 > //        A[0][2] = dAtom->getAxz();
499            
500 <          A[1][0] = dAtom->getAyx();
501 <          A[1][1] = dAtom->getAyy();
502 <          A[1][2] = dAtom->getAyz();
500 > //        A[1][0] = dAtom->getAyx();
501 > //        A[1][1] = dAtom->getAyy();
502 > //        A[1][2] = dAtom->getAyz();
503            
504 <          A[2][0] = dAtom->getAzx();
505 <          A[2][1] = dAtom->getAzy();
506 <          A[2][2] = dAtom->getAzz();
504 > //        A[2][0] = dAtom->getAzx();
505 > //        A[2][1] = dAtom->getAzy();
506 > //        A[2][2] = dAtom->getAzz();
507            
508            
509 <          // use the angular velocities to propagate the rotation matrix a
510 <          // full time step
509 > //        // use the angular velocities to propagate the rotation matrix a
510 > //        // full time step
511            
512            
513 <          angle = dt2 * ji[0] / dAtom->getIxx();
514 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
513 > //        angle = dt2 * ji[0] / dAtom->getIxx();
514 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
515            
516 <          angle = dt2 * ji[1] / dAtom->getIyy();
517 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
516 > //        angle = dt2 * ji[1] / dAtom->getIyy();
517 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
518            
519 <          angle = dt * ji[2] / dAtom->getIzz();
520 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
519 > //        angle = dt * ji[2] / dAtom->getIzz();
520 > //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
521            
522 <          angle = dt2 * ji[1] / dAtom->getIyy();
523 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
522 > //        angle = dt2 * ji[1] / dAtom->getIyy();
523 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
524            
525 <          angle = dt2 * ji[0] / dAtom->getIxx();
526 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
525 > //        angle = dt2 * ji[0] / dAtom->getIxx();
526 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
527            
528            
529 <          dAtom->setA( A );
530 <          dAtom->setJx( ji[0] );
531 <          dAtom->setJy( ji[1] );
532 <          dAtom->setJz( ji[2] );
533 <        }
529 > //        dAtom->setA( A );
530 > //        dAtom->setJx( ji[0] );
531 > //        dAtom->setJy( ji[1] );
532 > //        dAtom->setJz( ji[2] );
533 > //      }
534        }
535        
536        // calculate the forces
537        
538 <      myFF->doForces(calcPot,0);
538 >      myFF->doForces(calcPot,calcStress);
539        
540        for( i=0; i< nAtoms; i++ ){
541          
# Line 515 | Line 556 | void Symplectic::integrate( void ){
556   //      vy2 = vy * vy;
557   //      vz2 = vz * vz;
558          
559 <        if( atoms[i]->isDirectional() ){
559 > //      if( atoms[i]->isDirectional() ){
560  
561 <          dAtom = (DirectionalAtom *)atoms[i];
561 > //        dAtom = (DirectionalAtom *)atoms[i];
562            
563 <          // get and convert the torque to body frame
563 > //        // get and convert the torque to body frame
564            
565 <          Tb[0] = dAtom->getTx();
566 <          Tb[1] = dAtom->getTy();
567 <          Tb[2] = dAtom->getTz();
565 > //        Tb[0] = dAtom->getTx();
566 > //        Tb[1] = dAtom->getTy();
567 > //        Tb[2] = dAtom->getTz();
568            
569 <          dAtom->lab2Body( Tb );
569 > //        dAtom->lab2Body( Tb );
570            
571 <          // get the angular momentum, and complete the angular momentum
572 <          // half step
571 > //        // get the angular momentum, and complete the angular momentum
572 > //        // half step
573            
574 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
574 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
575 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
576 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
577            
578 <          jx2 = ji[0] * ji[0];
579 <          jy2 = ji[1] * ji[1];
580 <          jz2 = ji[2] * ji[2];
578 > //        jx2 = ji[0] * ji[0];
579 > //        jy2 = ji[1] * ji[1];
580 > //        jz2 = ji[2] * ji[2];
581            
582 <          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583 <            + (jz2 / dAtom->getIzz());
582 > //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
583 > //          + (jz2 / dAtom->getIzz());
584            
585 <          dAtom->setJx( ji[0] );
586 <          dAtom->setJy( ji[1] );
587 <          dAtom->setJz( ji[2] );
588 <        }
585 > //        dAtom->setJx( ji[0] );
586 > //        dAtom->setJy( ji[1] );
587 > //        dAtom->setJz( ji[2] );
588 > //      }
589        }
590 <      
590 >    
591 >      if (!strcasecmp( entry_plug->ensemble, "NVT"))
592 >        myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() );
593 >
594 >      if (!strcasecmp( entry_plug->ensemble, "NPT") )
595 >        myES->NoseHooverAndersonNPT( dt,
596 >                                     tStats->getKinetic(),
597 >                                     tStats->getPressure());
598 >  
599        time = tl + 1;
600        
601        if( entry_plug->setTemp ){
602          if( !(time % vel_n) ) tStats->velocitize();
603        }
604        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
605 <      if( !((time+1) % status_n) ) calcPot = 1;
606 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
605 >      if( !((time+1) % status_n) ) {
606 >        calcPot = 1;
607 >        // bitwise masking in case we need it for NPT
608 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 1;
609 >      }
610 >      if( !(time % status_n) ){
611 >        e_out->writeStat( time * dt );
612 >        calcPot = 0;
613 >        // bitwise masking in case we need it for NPT
614 >        calcStress = (!strcasecmp(entry_plug->ensemble,"NPT")) && 0;
615 >      }      
616      }
617    }
618  
# Line 581 | Line 639 | void Symplectic::rotate( int axes1, int axes2, double
639  
640    for(i=0; i<3; i++){
641      for(j=0; j<3; j++){
642 <      tempA[i][j] = A[i][j];
642 >      tempA[j][i] = A[i][j];
643      }
644    }
645  
# Line 640 | Line 698 | void Symplectic::rotate( int axes1, int axes2, double
698      for(j=0; j<3; j++){
699        A[j][i] = 0.0;
700        for(k=0; k<3; k++){
701 <        A[j][i] += tempA[k][i] * rot[j][k];
701 >        A[j][i] += tempA[i][k] * rot[j][k];
702        }
703      }
704    }

Diff Legend

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