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 428 by mmeineke, Thu Mar 27 21:07:14 2003 UTC vs.
Revision 475 by gezelter, Tue Apr 8 12:44:18 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 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 , 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;
282 <          
283 <          // get the atom's rotation matrix
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 <          A[0][0] = dAtom->getAxx();
286 <          A[0][1] = dAtom->getAxy();
287 <          A[0][2] = dAtom->getAxz();
303 > // //     // get the atom's rotation matrix
304            
305 <          A[1][0] = dAtom->getAyx();
306 <          A[1][1] = dAtom->getAyy();
307 <          A[1][2] = dAtom->getAyz();
305 > // //     A[0][0] = dAtom->getAxx();
306 > // //     A[0][1] = dAtom->getAxy();
307 > // //     A[0][2] = dAtom->getAxz();
308            
309 <          A[2][0] = dAtom->getAzx();
310 <          A[2][1] = dAtom->getAzy();
311 <          A[2][2] = dAtom->getAzz();
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();
316            
298 –          // use the angular velocities to propagate the rotation matrix a
299 –          // full time step
317            
318 + // //     // use the angular velocities to propagate the rotation matrix a
319 + // //     // full time step
320            
302 –          angle = dt2 * ji[0] / dAtom->getIxx();
303 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
321            
322 <          angle = dt2 * ji[1] / dAtom->getIyy();
323 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
322 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
323 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
324            
325 <          angle = dt * ji[2] / dAtom->getIzz();
326 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-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[1] / dAtom->getIyy();
329 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-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[0] / dAtom->getIxx();
332 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
331 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
332 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
333 >          
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 <      }
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 , 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 , 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 , 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