| 30 |
|
} |
| 31 |
|
|
| 32 |
|
|
| 33 |
< |
Verlet::Verlet( SimInfo &info ){ |
| 33 |
> |
Verlet::Verlet( SimInfo &info, ForceFields* the_ff ){ |
| 34 |
|
|
| 35 |
|
// get what information we need from the SimInfo object |
| 36 |
|
|
| 37 |
|
entry_plug = &info; |
| 38 |
+ |
myFF = the_ff; |
| 39 |
|
|
| 40 |
+ |
|
| 41 |
|
c_natoms = info.n_atoms; |
| 42 |
|
c_atoms = info.atoms; |
| 43 |
|
c_sr_interactions = info.sr_interactions; |
| 42 |
– |
longRange = info.longRange; |
| 44 |
|
c_n_SRI = info.n_SRI; |
| 45 |
|
c_is_constrained = 0; |
| 46 |
|
c_box_x = info.box_x; |
| 124 |
|
} |
| 125 |
|
|
| 126 |
|
|
| 126 |
– |
|
| 127 |
– |
void Verlet::integrate_b( double time_length, double dt, |
| 128 |
– |
int n_bond_0, int n_bond_f, |
| 129 |
– |
int n_bend_0, int n_bend_f, |
| 130 |
– |
int n_torsion_0, int n_torsion_f, |
| 131 |
– |
bool do_bonds, bool do_bends, bool do_torsions, |
| 132 |
– |
bool do_LRI ){ |
| 133 |
– |
|
| 134 |
– |
// double percent_tolerance = 0.001; |
| 135 |
– |
// int max_iterations = 10000; |
| 136 |
– |
|
| 137 |
– |
int i, j; /* loop counters */ |
| 138 |
– |
double n_loops = time_length / dt; |
| 139 |
– |
|
| 140 |
– |
// the first time integrate is called, the forces need to be initialized |
| 141 |
– |
|
| 142 |
– |
if(is_first){ |
| 143 |
– |
is_first = 0; |
| 144 |
– |
|
| 145 |
– |
for(i = 0; i < c_natoms; i++){ |
| 146 |
– |
c_atoms[i]->zeroForces(); |
| 147 |
– |
} |
| 148 |
– |
|
| 149 |
– |
if( do_bonds ){ |
| 150 |
– |
for(i = n_bond_0; i <= n_bond_f; i++){ |
| 151 |
– |
c_sr_interactions[i]->calc_forces(); |
| 152 |
– |
} |
| 153 |
– |
} |
| 154 |
– |
|
| 155 |
– |
if( do_bends ){ |
| 156 |
– |
for(i = n_bend_0; i <= n_bend_f; i++){ |
| 157 |
– |
c_sr_interactions[i]->calc_forces(); |
| 158 |
– |
} |
| 159 |
– |
} |
| 160 |
– |
|
| 161 |
– |
if( do_torsions ){ |
| 162 |
– |
for(i = n_torsion_0; i <= n_torsion_f; i++){ |
| 163 |
– |
c_sr_interactions[i]->calc_forces(); |
| 164 |
– |
} |
| 165 |
– |
} |
| 166 |
– |
|
| 167 |
– |
if( do_LRI ) longRange->calc_forces(); |
| 168 |
– |
} |
| 169 |
– |
|
| 170 |
– |
for(i = 0; i < n_loops; i++){ |
| 171 |
– |
|
| 172 |
– |
move_a( dt ); |
| 173 |
– |
|
| 174 |
– |
// calculate the forces |
| 175 |
– |
|
| 176 |
– |
for(j = 0; j < c_natoms; j++){ |
| 177 |
– |
c_atoms[j]->zeroForces(); |
| 178 |
– |
} |
| 179 |
– |
|
| 180 |
– |
|
| 181 |
– |
if( do_bonds ){ |
| 182 |
– |
for(i = n_bond_0; i <= n_bond_f; i++){ |
| 183 |
– |
c_sr_interactions[i]->calc_forces(); |
| 184 |
– |
} |
| 185 |
– |
} |
| 186 |
– |
|
| 187 |
– |
if( do_bends ){ |
| 188 |
– |
for(i = n_bend_0; i <= n_bend_f; i++){ |
| 189 |
– |
c_sr_interactions[i]->calc_forces(); |
| 190 |
– |
} |
| 191 |
– |
} |
| 192 |
– |
|
| 193 |
– |
if( do_torsions ){ |
| 194 |
– |
for(i = n_torsion_0; i <= n_torsion_f; i++){ |
| 195 |
– |
c_sr_interactions[i]->calc_forces(); |
| 196 |
– |
} |
| 197 |
– |
} |
| 198 |
– |
|
| 199 |
– |
if( do_LRI ) longRange->calc_forces(); |
| 200 |
– |
|
| 201 |
– |
|
| 202 |
– |
// complete the verlet move |
| 203 |
– |
|
| 204 |
– |
move_b( dt ); |
| 205 |
– |
} |
| 206 |
– |
} |
| 207 |
– |
|
| 208 |
– |
|
| 127 |
|
void Verlet::integrate( void ){ |
| 128 |
|
|
| 129 |
|
int i, j; /* loop counters */ |
| 130 |
< |
|
| 130 |
> |
int calcPot; |
| 131 |
> |
|
| 132 |
|
double kE; |
| 133 |
|
|
| 134 |
|
double *Rx = new double[c_natoms]; |
| 164 |
|
// the first time integrate is called, the forces need to be initialized |
| 165 |
|
|
| 166 |
|
|
| 167 |
< |
for(i = 0; i < c_natoms; i++){ |
| 249 |
< |
c_atoms[i]->zeroForces(); |
| 250 |
< |
} |
| 167 |
> |
myFF->doForces(1); |
| 168 |
|
|
| 252 |
– |
for(i = 0; i < c_n_SRI; i++){ |
| 253 |
– |
c_sr_interactions[i]->calc_forces(); |
| 254 |
– |
} |
| 255 |
– |
|
| 256 |
– |
longRange->calc_forces(); |
| 257 |
– |
|
| 169 |
|
if( entry_plug->setTemp ){ |
| 170 |
|
tStats->velocitize(); |
| 171 |
|
} |
| 172 |
|
|
| 173 |
|
dump_out->writeDump( 0.0 ); |
| 174 |
+ |
|
| 175 |
|
e_out->writeStat( 0.0 ); |
| 176 |
|
|
| 177 |
+ |
calcPot = 0; |
| 178 |
+ |
|
| 179 |
|
if( c_is_constrained ){ |
| 180 |
|
for(i = 0; i < n_loops; i++){ |
| 181 |
|
|
| 216 |
|
|
| 217 |
|
// calculate the forces |
| 218 |
|
|
| 219 |
< |
for(j = 0; j < c_natoms; j++){ |
| 306 |
< |
c_atoms[j]->zeroForces(); |
| 307 |
< |
} |
| 219 |
> |
myFF->doForces(calcPot); |
| 220 |
|
|
| 309 |
– |
for(j = 0; j < c_n_SRI; j++){ |
| 310 |
– |
c_sr_interactions[j]->calc_forces(); |
| 311 |
– |
} |
| 312 |
– |
|
| 313 |
– |
longRange->calc_forces(); |
| 314 |
– |
|
| 221 |
|
// finish the constrain move ( same as above. ) |
| 222 |
|
|
| 223 |
|
for( j=0; j<c_natoms; j++ ){ |
| 258 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 259 |
|
} |
| 260 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 261 |
< |
if( !(time % status_n) ) e_out->writeStat( time * dt ); |
| 261 |
> |
if( !((time+1) % status_n) ) calcPot = 1; |
| 262 |
> |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 263 |
|
} |
| 264 |
|
} |
| 265 |
|
else{ |
| 269 |
|
|
| 270 |
|
// calculate the forces |
| 271 |
|
|
| 272 |
< |
for(j = 0; j < c_natoms; j++){ |
| 366 |
< |
c_atoms[j]->zeroForces(); |
| 367 |
< |
} |
| 368 |
< |
|
| 369 |
< |
for(j = 0; j < c_n_SRI; j++){ |
| 370 |
< |
c_sr_interactions[j]->calc_forces(); |
| 371 |
< |
} |
| 372 |
< |
|
| 373 |
< |
longRange->calc_forces(); |
| 272 |
> |
myFF->doForces(calcPot); |
| 273 |
|
|
| 274 |
|
// complete the verlet move |
| 275 |
|
|
| 280 |
|
if( entry_plug->setTemp ){ |
| 281 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
| 282 |
|
} |
| 283 |
< |
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 284 |
< |
if( !(time % status_n) ) e_out->writeStat( time * dt ); |
| 283 |
> |
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
| 284 |
> |
if( !((time+1) % status_n) ) calcPot = 1; |
| 285 |
> |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
| 286 |
|
} |
| 287 |
|
} |
| 288 |
|
|
| 335 |
|
const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 |
| 336 |
|
|
| 337 |
|
double vx, vy, vz; |
| 438 |
– |
double v_sqr; |
| 338 |
|
int mb; |
| 339 |
|
double h_dt = 0.5 * dt; |
| 340 |
|
|
