| 5 |
|
#include "Thermo.hpp" |
| 6 |
|
#include "ReadWrite.hpp" |
| 7 |
|
#include "ForceFields.hpp" |
| 8 |
+ |
#include "ExtendedSystem.hpp" |
| 9 |
|
#include "simError.h" |
| 10 |
|
|
| 11 |
|
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 |
|
molecules = entry_plug->molecules; |
| 51 |
|
mass = new double[entry_plug->n_atoms]; |
| 52 |
|
for(int i = 0; i < entry_plug->n_atoms; i++){ |
| 53 |
|
mass[i] = entry_plug->atoms[i]->getMass(); |
| 54 |
< |
} |
| 52 |
< |
|
| 53 |
< |
|
| 54 |
> |
} |
| 55 |
|
|
| 56 |
|
// check for constraints |
| 57 |
|
|
| 184 |
|
int status_n = (int)( statusTime / dt ); |
| 185 |
|
int vel_n = (int)( thermalTime / dt ); |
| 186 |
|
|
| 187 |
< |
int calcPot; |
| 187 |
> |
int calcPot, calcStress; |
| 188 |
|
|
| 189 |
|
Thermo *tStats = new Thermo( entry_plug ); |
| 190 |
|
|
| 197 |
|
|
| 198 |
|
// initialize the forces the before the first step |
| 199 |
|
|
| 200 |
< |
myFF->doForces(1,0); |
| 200 |
> |
myFF->doForces(1,1); |
| 201 |
|
|
| 202 |
|
if( entry_plug->setTemp ){ |
| 203 |
|
|
| 325 |
|
|
| 326 |
|
// calculate the forces |
| 327 |
|
|
| 328 |
< |
myFF->doForces(calcPot, 0); |
| 328 |
> |
myFF->doForces(calcPot, calcStress); |
| 329 |
|
|
| 330 |
|
// move b |
| 331 |
|
|
| 396 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 397 |
|
} |
| 398 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 399 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 400 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 399 |
> |
if( !((time+1) % status_n) ) { |
| 400 |
> |
calcPot = 1; |
| 401 |
> |
calcStress = 1; |
| 402 |
> |
} |
| 403 |
> |
if( !(time % status_n) ){ |
| 404 |
> |
e_out->writeStat( time * dt ); |
| 405 |
> |
calcPot = 0; |
| 406 |
> |
calcStress = 0; |
| 407 |
> |
} |
| 408 |
|
} |
| 409 |
|
} |
| 410 |
|
else{ |
| 502 |
|
|
| 503 |
|
// calculate the forces |
| 504 |
|
|
| 505 |
< |
myFF->doForces(calcPot,0); |
| 505 |
> |
myFF->doForces(calcPot,calcStress); |
| 506 |
|
|
| 507 |
|
for( i=0; i< nAtoms; i++ ){ |
| 508 |
|
|
| 561 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 562 |
|
} |
| 563 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 564 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 565 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 564 |
> |
if( !((time+1) % status_n) ) { |
| 565 |
> |
calcPot = 1; |
| 566 |
> |
calcStress = 1; |
| 567 |
> |
} |
| 568 |
> |
if( !(time % status_n) ){ |
| 569 |
> |
e_out->writeStat( time * dt ); |
| 570 |
> |
calcPot = 0; |
| 571 |
> |
calcStress = 0; |
| 572 |
> |
} |
| 573 |
|
} |
| 574 |
|
} |
| 575 |
|
|