| 7 |
|
#include "SimInfo.hpp" |
| 8 |
|
#include "Thermo.hpp" |
| 9 |
|
#include "ReadWrite.hpp" |
| 10 |
+ |
#include "ExtendedSystem.hpp" |
| 11 |
|
|
| 12 |
|
extern "C"{ |
| 13 |
|
|
| 30 |
|
} |
| 31 |
|
|
| 32 |
|
|
| 33 |
< |
Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){ |
| 33 |
> |
Verlet::Verlet( SimInfo &info, ForceFields* the_ff, ExtendedSystem* the_es ){ |
| 34 |
|
|
| 35 |
|
// get what information we need from the SimInfo object |
| 36 |
|
|
| 37 |
|
entry_plug = &info; |
| 38 |
|
myFF = the_ff; |
| 39 |
< |
|
| 39 |
> |
myES = the_es; |
| 40 |
|
|
| 41 |
|
c_natoms = info.n_atoms; |
| 42 |
|
c_atoms = info.atoms; |
| 75 |
|
SRI** theArray; |
| 76 |
|
for(int i = 0; i < nMols; i++){ |
| 77 |
|
|
| 78 |
< |
theArray = molecules[i].getMyBonds(); |
| 79 |
< |
for(int j=0; j<molecules[i].getNbonds(); j++){ |
| 78 |
> |
theArray = (SRI**) molecules[i].getMyBonds(); |
| 79 |
> |
for(int j=0; j<molecules[i].getNBonds(); j++){ |
| 80 |
|
|
| 81 |
|
constrained = theArray[j]->is_constrained(); |
| 82 |
|
|
| 92 |
|
} |
| 93 |
|
} |
| 94 |
|
|
| 95 |
< |
theArray = molecules[i].getMyBends(); |
| 96 |
< |
for(int j=0; j<molecules[i].getNbends(); j++){ |
| 95 |
> |
theArray = (SRI**) molecules[i].getMyBends(); |
| 96 |
> |
for(int j=0; j<molecules[i].getNBends(); j++){ |
| 97 |
|
|
| 98 |
|
constrained = theArray[j]->is_constrained(); |
| 99 |
|
|
| 109 |
|
} |
| 110 |
|
} |
| 111 |
|
|
| 112 |
< |
theArray = molecules[i].getMyTorsions(); |
| 113 |
< |
for(int j=0; j<molecules[i].getNtorsions(); j++){ |
| 112 |
> |
theArray = (SRI**) molecules[i].getMyTorsions(); |
| 113 |
> |
for(int j=0; j<molecules[i].getNTorsions(); j++){ |
| 114 |
|
|
| 115 |
|
constrained = theArray[j]->is_constrained(); |
| 116 |
|
|
| 167 |
|
void Verlet::integrate( void ){ |
| 168 |
|
|
| 169 |
|
int i, j; /* loop counters */ |
| 170 |
< |
int calcPot; |
| 170 |
> |
int calcPot, calcStress; |
| 171 |
|
|
| 172 |
|
double kE; |
| 173 |
|
|
| 204 |
|
// the first time integrate is called, the forces need to be initialized |
| 205 |
|
|
| 206 |
|
|
| 207 |
< |
myFF->doForces(1,0); |
| 207 |
> |
myFF->doForces(1,1); |
| 208 |
|
|
| 209 |
|
if( entry_plug->setTemp ){ |
| 210 |
|
tStats->velocitize(); |
| 215 |
|
e_out->writeStat( 0.0 ); |
| 216 |
|
|
| 217 |
|
calcPot = 0; |
| 218 |
+ |
calcStress = 0; |
| 219 |
|
|
| 220 |
|
if( c_is_constrained ){ |
| 221 |
|
for(i = 0; i < n_loops; i++){ |
| 257 |
|
|
| 258 |
|
// calculate the forces |
| 259 |
|
|
| 260 |
< |
myFF->doForces(calcPot,0); |
| 260 |
> |
myFF->doForces(calcPot,calcStress); |
| 261 |
|
|
| 262 |
|
// finish the constrain move ( same as above. ) |
| 263 |
|
|
| 299 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 300 |
|
} |
| 301 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 302 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 303 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 302 |
> |
if( !((time+1) % status_n) ) { |
| 303 |
> |
calcPot = 1; |
| 304 |
> |
calcStress = 1; |
| 305 |
> |
} |
| 306 |
> |
if( !(time % status_n) ){ |
| 307 |
> |
e_out->writeStat( time * dt ); |
| 308 |
> |
calcPot = 0; |
| 309 |
> |
calcStress = 0; |
| 310 |
> |
} |
| 311 |
|
} |
| 312 |
|
} |
| 313 |
|
else{ |
| 317 |
|
|
| 318 |
|
// calculate the forces |
| 319 |
|
|
| 320 |
< |
myFF->doForces(calcPot,0); |
| 320 |
> |
myFF->doForces(calcPot,calcStress); |
| 321 |
|
|
| 322 |
|
// complete the verlet move |
| 323 |
|
|
| 329 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 330 |
|
} |
| 331 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 332 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
| 333 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 332 |
> |
if( !((time+1) % status_n) ) { |
| 333 |
> |
calcPot = 1; |
| 334 |
> |
calcStress = 1; |
| 335 |
> |
} |
| 336 |
> |
if( !(time % status_n) ){ |
| 337 |
> |
e_out->writeStat( time * dt ); |
| 338 |
> |
calcPot = 0; |
| 339 |
> |
calcStress = 0; |
| 340 |
> |
} |
| 341 |
|
} |
| 342 |
|
} |
| 343 |
|
|
