30 |
|
} |
31 |
|
|
32 |
|
|
33 |
< |
Verlet::Verlet( SimInfo &info ){ |
33 |
> |
Verlet::Verlet( SimInfo &info, ForceField* 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; |
122 |
|
|
123 |
|
delete[] c_mass; |
124 |
|
c_mass = 0; |
123 |
– |
} |
124 |
– |
|
125 |
– |
|
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 |
– |
} |
125 |
|
} |
126 |
|
|
127 |
|
|
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(); |
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 |
|
} |
213 |
|
|
214 |
|
// calculate the forces |
215 |
|
|
216 |
< |
for(j = 0; j < c_natoms; j++){ |
306 |
< |
c_atoms[j]->zeroForces(); |
307 |
< |
} |
216 |
> |
myFF->doForces(); |
217 |
|
|
309 |
– |
for(j = 0; j < c_n_SRI; j++){ |
310 |
– |
c_sr_interactions[j]->calc_forces(); |
311 |
– |
} |
312 |
– |
|
313 |
– |
longRange->calc_forces(); |
314 |
– |
|
218 |
|
// finish the constrain move ( same as above. ) |
219 |
|
|
220 |
|
for( j=0; j<c_natoms; j++ ){ |
265 |
|
|
266 |
|
// calculate the forces |
267 |
|
|
268 |
< |
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(); |
268 |
> |
myFF->doForces(); |
269 |
|
|
270 |
|
// complete the verlet move |
271 |
|
|
330 |
|
const double e_convert = 4.184e-4; // converts kcal/mol -> amu*A^2/fs^2 |
331 |
|
|
332 |
|
double vx, vy, vz; |
438 |
– |
double v_sqr; |
333 |
|
int mb; |
334 |
|
double h_dt = 0.5 * dt; |
335 |
|
|