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 472 by mmeineke, Mon Apr 7 21:16:35 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 >  std::cerr << "about to call new thermo\n";
196  
197 <  StatWriter*  e_out    = new StatWriter( entry_plug );
198 <  DumpWriter*  dump_out = new DumpWriter( entry_plug );
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 261 | Line 272 | void Symplectic::integrate( void ){
272        }
273  
274  
275 <      for( i=0; i<nAtoms; i++ ){
276 <        if( atoms[i]->isDirectional() ){
277 <                  
278 <          dAtom = (DirectionalAtom *)atoms[i];
275 > //       for( i=0; i<nAtoms; i++ ){
276 > // //   if( atoms[i]->isDirectional() ){
277 >                  
278 > // //     dAtom = (DirectionalAtom *)atoms[i];
279            
280 <          // get and convert the torque to body frame
280 > // //     // get and convert the torque to body frame
281            
282 <          Tb[0] = dAtom->getTx();
283 <          Tb[1] = dAtom->getTy();
284 <          Tb[2] = dAtom->getTz();
282 > // //     Tb[0] = dAtom->getTx();
283 > // //     Tb[1] = dAtom->getTy();
284 > // //     Tb[2] = dAtom->getTz();
285            
286 <          dAtom->lab2Body( Tb );
286 > // //     dAtom->lab2Body( Tb );
287            
288 <          // get the angular momentum, and propagate a half step
288 > // //     // get the angular momentum, and propagate a half step
289            
290 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
291 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
292 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
290 > // //     ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
291 > // //     ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
292 > // //     ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
293            
294 <          // get the atom's rotation matrix
294 > // //     // get the atom's rotation matrix
295            
296 <          A[0][0] = dAtom->getAxx();
297 <          A[0][1] = dAtom->getAxy();
298 <          A[0][2] = dAtom->getAxz();
296 > // //     A[0][0] = dAtom->getAxx();
297 > // //     A[0][1] = dAtom->getAxy();
298 > // //     A[0][2] = dAtom->getAxz();
299            
300 <          A[1][0] = dAtom->getAyx();
301 <          A[1][1] = dAtom->getAyy();
302 <          A[1][2] = dAtom->getAyz();
300 > // //     A[1][0] = dAtom->getAyx();
301 > // //     A[1][1] = dAtom->getAyy();
302 > // //     A[1][2] = dAtom->getAyz();
303            
304 <          A[2][0] = dAtom->getAzx();
305 <          A[2][1] = dAtom->getAzy();
306 <          A[2][2] = dAtom->getAzz();
304 > // //     A[2][0] = dAtom->getAzx();
305 > // //     A[2][1] = dAtom->getAzy();
306 > // //     A[2][2] = dAtom->getAzz();
307            
308            
309 <          // use the angular velocities to propagate the rotation matrix a
310 <          // full time step
309 > // //     // use the angular velocities to propagate the rotation matrix a
310 > // //     // full time step
311            
312            
313 <          angle = dt2 * ji[0] / dAtom->getIxx();
314 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
313 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
314 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
315            
316 <          angle = dt2 * ji[1] / dAtom->getIyy();
317 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
316 > // //     angle = dt2 * ji[1] / dAtom->getIyy();
317 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
318            
319 <          angle = dt * ji[2] / dAtom->getIzz();
320 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
319 > // //     angle = dt * ji[2] / dAtom->getIzz();
320 > // //     this->rotate( 0, 1, angle, ji, A ); // rotate about the z-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[1] / dAtom->getIyy();
323 > // //     this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
324            
325 <          angle = dt2 * ji[0] / dAtom->getIxx();
326 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
325 > // //     angle = dt2 * ji[0] / dAtom->getIxx();
326 > // //     this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
327            
328            
329 <          dAtom->setA( A );
330 <          dAtom->setJx( ji[0] );
331 <          dAtom->setJy( ji[1] );
332 <          dAtom->setJz( ji[2] );
333 <        }
334 <      }
329 > // //     dAtom->setA( A );
330 > // //     dAtom->setJx( ji[0] );
331 > // //     dAtom->setJy( ji[1] );
332 > // //     dAtom->setJz( ji[2] );
333 > // //   }
334 > //       }
335        
336        // calculate the forces
337        
338 <      myFF->doForces(calcPot, 0);
338 >      myFF->doForces(calcPot, calcStress);
339        
340        // move b
341  
# Line 362 | Line 373 | void Symplectic::integrate( void ){
373          atoms[j]->set_vz(Vz[j]);
374        }
375        
376 <      for( i=0; i< nAtoms; i++ ){
376 > //       for( i=0; i< nAtoms; i++ ){
377  
378 <        if( atoms[i]->isDirectional() ){
378 > //      if( atoms[i]->isDirectional() ){
379  
380 <          dAtom = (DirectionalAtom *)atoms[i];
380 > //        dAtom = (DirectionalAtom *)atoms[i];
381            
382 <          // get and convert the torque to body frame
382 > //        // get and convert the torque to body frame
383            
384 <          Tb[0] = dAtom->getTx();
385 <          Tb[1] = dAtom->getTy();
386 <          Tb[2] = dAtom->getTz();
384 > //        Tb[0] = dAtom->getTx();
385 > //        Tb[1] = dAtom->getTy();
386 > //        Tb[2] = dAtom->getTz();
387            
388 <          dAtom->lab2Body( Tb );
388 > //        dAtom->lab2Body( Tb );
389            
390 <          // get the angular momentum, and complete the angular momentum
391 <          // half step
390 > //        // get the angular momentum, and complete the angular momentum
391 > //        // half step
392            
393 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
393 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
394 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
395 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
396            
397 <          dAtom->setJx( ji[0] );
398 <          dAtom->setJy( ji[1] );
399 <          dAtom->setJz( ji[2] );
400 <        }
401 <      }
397 > //        dAtom->setJx( ji[0] );
398 > //        dAtom->setJy( ji[1] );
399 > //        dAtom->setJz( ji[2] );
400 > //      }
401 > //       }
402      
403        time = tl + 1;
404        
# Line 395 | Line 406 | void Symplectic::integrate( void ){
406          if( !(time % vel_n) ) tStats->velocitize();
407        }
408        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
409 <      if( !((time+1) % status_n) ) calcPot = 1;
410 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
409 >      if( !((time+1) % status_n) ) {
410 >        calcPot = 1;
411 >        calcStress = 1;
412 >      }
413 >      if( !(time % status_n) ){
414 >        e_out->writeStat( time * dt );
415 >        calcPot = 0;
416 >        calcStress = 0;
417 >      }
418      }
419    }
420    else{
# Line 432 | Line 450 | void Symplectic::integrate( void ){
450          atoms[i]->set_vy( vy );
451          atoms[i]->set_vz( vz );
452          
453 <        if( atoms[i]->isDirectional() ){
453 > //      if( atoms[i]->isDirectional() ){
454  
455 <          dAtom = (DirectionalAtom *)atoms[i];
455 > //        dAtom = (DirectionalAtom *)atoms[i];
456            
457 <          // get and convert the torque to body frame
457 > //        // get and convert the torque to body frame
458            
459 <          Tb[0] = dAtom->getTx();
460 <          Tb[1] = dAtom->getTy();
461 <          Tb[2] = dAtom->getTz();
444 <          
445 <          dAtom->lab2Body( Tb );
459 > //        Tb[0] = dAtom->getTx();
460 > //        Tb[1] = dAtom->getTy();
461 > //        Tb[2] = dAtom->getTz();
462            
463 <          // get the angular momentum, and propagate a half step
463 > //        dAtom->lab2Body( Tb );
464            
465 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
450 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
451 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
465 > //        // get the angular momentum, and propagate a half step
466            
467 <          // get the atom's rotation matrix
467 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
468 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
469 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
470            
471 <          A[0][0] = dAtom->getAxx();
456 <          A[0][1] = dAtom->getAxy();
457 <          A[0][2] = dAtom->getAxz();
471 > //        // get the atom's rotation matrix
472            
473 <          A[1][0] = dAtom->getAyx();
474 <          A[1][1] = dAtom->getAyy();
475 <          A[1][2] = dAtom->getAyz();
473 > //        A[0][0] = dAtom->getAxx();
474 > //        A[0][1] = dAtom->getAxy();
475 > //        A[0][2] = dAtom->getAxz();
476            
477 <          A[2][0] = dAtom->getAzx();
478 <          A[2][1] = dAtom->getAzy();
479 <          A[2][2] = dAtom->getAzz();
477 > //        A[1][0] = dAtom->getAyx();
478 > //        A[1][1] = dAtom->getAyy();
479 > //        A[1][2] = dAtom->getAyz();
480            
481 + //        A[2][0] = dAtom->getAzx();
482 + //        A[2][1] = dAtom->getAzy();
483 + //        A[2][2] = dAtom->getAzz();
484            
468 –          // use the angular velocities to propagate the rotation matrix a
469 –          // full time step
485            
486 + //        // use the angular velocities to propagate the rotation matrix a
487 + //        // full time step
488            
472 –          angle = dt2 * ji[0] / dAtom->getIxx();
473 –          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
489            
490 <          angle = dt2 * ji[1] / dAtom->getIyy();
491 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
490 > //        angle = dt2 * ji[0] / dAtom->getIxx();
491 > //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
492            
493 <          angle = dt * ji[2] / dAtom->getIzz();
494 <          this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
493 > //        angle = dt2 * ji[1] / dAtom->getIyy();
494 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
495            
496 <          angle = dt2 * ji[1] / dAtom->getIyy();
497 <          this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
496 > //        angle = dt * ji[2] / dAtom->getIzz();
497 > //        this->rotate( 0, 1, angle, ji, A ); // rotate about the z-axis
498            
499 <          angle = dt2 * ji[0] / dAtom->getIxx();
500 <          this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
499 > //        angle = dt2 * ji[1] / dAtom->getIyy();
500 > //        this->rotate( 2, 0, angle, ji, A ); // rotate about the y-axis
501            
502 + //        angle = dt2 * ji[0] / dAtom->getIxx();
503 + //        this->rotate( 1, 2, angle, ji, A ); // rotate about the x-axis
504            
505 <          dAtom->setA( A );
506 <          dAtom->setJx( ji[0] );
507 <          dAtom->setJy( ji[1] );
508 <          dAtom->setJz( ji[2] );
509 <        }
505 >          
506 > //        dAtom->setA( A );
507 > //        dAtom->setJx( ji[0] );
508 > //        dAtom->setJy( ji[1] );
509 > //        dAtom->setJz( ji[2] );
510 > //      }
511        }
512        
513        // calculate the forces
514        
515 <      myFF->doForces(calcPot,0);
515 >      myFF->doForces(calcPot,calcStress);
516        
517        for( i=0; i< nAtoms; i++ ){
518          
# Line 515 | Line 533 | void Symplectic::integrate( void ){
533   //      vy2 = vy * vy;
534   //      vz2 = vz * vz;
535          
536 <        if( atoms[i]->isDirectional() ){
536 > //      if( atoms[i]->isDirectional() ){
537  
538 <          dAtom = (DirectionalAtom *)atoms[i];
538 > //        dAtom = (DirectionalAtom *)atoms[i];
539            
540 <          // get and convert the torque to body frame
540 > //        // get and convert the torque to body frame
541            
542 <          Tb[0] = dAtom->getTx();
543 <          Tb[1] = dAtom->getTy();
544 <          Tb[2] = dAtom->getTz();
542 > //        Tb[0] = dAtom->getTx();
543 > //        Tb[1] = dAtom->getTy();
544 > //        Tb[2] = dAtom->getTz();
545            
546 <          dAtom->lab2Body( Tb );
546 > //        dAtom->lab2Body( Tb );
547            
548 <          // get the angular momentum, and complete the angular momentum
549 <          // half step
548 > //        // get the angular momentum, and complete the angular momentum
549 > //        // half step
550            
551 <          ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 <          ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 <          ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
551 > //        ji[0] = dAtom->getJx() + ( dt2 * Tb[0] ) * e_convert;
552 > //        ji[1] = dAtom->getJy() + ( dt2 * Tb[1] ) * e_convert;
553 > //        ji[2] = dAtom->getJz() + ( dt2 * Tb[2] ) * e_convert;
554            
555 <          jx2 = ji[0] * ji[0];
556 <          jy2 = ji[1] * ji[1];
557 <          jz2 = ji[2] * ji[2];
555 > //        jx2 = ji[0] * ji[0];
556 > //        jy2 = ji[1] * ji[1];
557 > //        jz2 = ji[2] * ji[2];
558            
559 <          rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 <            + (jz2 / dAtom->getIzz());
559 > //        rot_kE += (jx2 / dAtom->getIxx()) + (jy2 / dAtom->getIyy())
560 > //          + (jz2 / dAtom->getIzz());
561            
562 <          dAtom->setJx( ji[0] );
563 <          dAtom->setJy( ji[1] );
564 <          dAtom->setJz( ji[2] );
565 <        }
562 > //        dAtom->setJx( ji[0] );
563 > //        dAtom->setJy( ji[1] );
564 > //        dAtom->setJz( ji[2] );
565 > //      }
566        }
567        
568        time = tl + 1;
# Line 553 | Line 571 | void Symplectic::integrate( void ){
571          if( !(time % vel_n) ) tStats->velocitize();
572        }
573        if( !(time % sample_n) ) dump_out->writeDump( time * dt );
574 <      if( !((time+1) % status_n) ) calcPot = 1;
575 <      if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; }
574 >      if( !((time+1) % status_n) ) {
575 >        calcPot = 1;
576 >        calcStress = 1;
577 >      }
578 >      if( !(time % status_n) ){
579 >        e_out->writeStat( time * dt );
580 >        calcPot = 0;
581 >        calcStress = 0;
582 >      }      
583      }
584    }
585  
# Line 581 | Line 606 | void Symplectic::rotate( int axes1, int axes2, double
606  
607    for(i=0; i<3; i++){
608      for(j=0; j<3; j++){
609 <      tempA[i][j] = A[i][j];
609 >      tempA[j][i] = A[i][j];
610      }
611    }
612  
# Line 640 | Line 665 | void Symplectic::rotate( int axes1, int axes2, double
665      for(j=0; j<3; j++){
666        A[j][i] = 0.0;
667        for(k=0; k<3; k++){
668 <        A[j][i] += tempA[k][i] * rot[j][k];
668 >        A[j][i] += tempA[i][k] * rot[j][k];
669        }
670      }
671    }

Diff Legend

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