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 |
< |
|
40 |
< |
|
39 |
> |
myES = the_es; |
40 |
> |
|
41 |
|
c_natoms = info.n_atoms; |
42 |
|
c_atoms = info.atoms; |
43 |
< |
c_sr_interactions = info.sr_interactions; |
44 |
< |
c_n_SRI = info.n_SRI; |
43 |
> |
nMols = info.n_mol; |
44 |
> |
molecules = info.molecules; |
45 |
|
c_is_constrained = 0; |
46 |
|
c_box_x = info.box_x; |
47 |
|
c_box_y = info.box_y; |
68 |
|
|
69 |
|
Constraint *temp_con; |
70 |
|
Constraint *dummy_plug; |
71 |
< |
temp_con = new Constraint[c_n_SRI]; |
71 |
> |
temp_con = new Constraint[info.n_SRI]; |
72 |
|
|
73 |
|
c_n_constrained = 0; |
74 |
|
int constrained = 0; |
75 |
< |
|
76 |
< |
for(int i = 0; i < c_n_SRI; i++){ |
75 |
> |
SRI** theArray; |
76 |
> |
for(int i = 0; i < nMols; i++){ |
77 |
|
|
78 |
< |
constrained = c_sr_interactions[i]->is_constrained(); |
78 |
> |
theArray = (SRI**) molecules[i].getMyBonds(); |
79 |
> |
for(int j=0; j<molecules[i].getNBonds(); j++){ |
80 |
> |
|
81 |
> |
constrained = theArray[j]->is_constrained(); |
82 |
> |
|
83 |
> |
if(constrained){ |
84 |
> |
|
85 |
> |
dummy_plug = theArray[j]->get_constraint(); |
86 |
> |
temp_con[c_n_constrained].set_a( dummy_plug->get_a() ); |
87 |
> |
temp_con[c_n_constrained].set_b( dummy_plug->get_b() ); |
88 |
> |
temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() ); |
89 |
> |
|
90 |
> |
c_n_constrained++; |
91 |
> |
constrained = 0; |
92 |
> |
} |
93 |
> |
} |
94 |
|
|
95 |
< |
if(constrained){ |
95 |
> |
theArray = (SRI**) molecules[i].getMyBends(); |
96 |
> |
for(int j=0; j<molecules[i].getNBends(); j++){ |
97 |
|
|
98 |
< |
dummy_plug = c_sr_interactions[i]->get_constraint(); |
99 |
< |
temp_con[c_n_constrained].set_a( dummy_plug->get_a() ); |
100 |
< |
temp_con[c_n_constrained].set_b( dummy_plug->get_b() ); |
101 |
< |
temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() ); |
98 |
> |
constrained = theArray[j]->is_constrained(); |
99 |
> |
|
100 |
> |
if(constrained){ |
101 |
> |
|
102 |
> |
dummy_plug = theArray[j]->get_constraint(); |
103 |
> |
temp_con[c_n_constrained].set_a( dummy_plug->get_a() ); |
104 |
> |
temp_con[c_n_constrained].set_b( dummy_plug->get_b() ); |
105 |
> |
temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() ); |
106 |
> |
|
107 |
> |
c_n_constrained++; |
108 |
> |
constrained = 0; |
109 |
> |
} |
110 |
> |
} |
111 |
|
|
112 |
< |
c_n_constrained++; |
113 |
< |
constrained = 0; |
112 |
> |
theArray = (SRI**) molecules[i].getMyTorsions(); |
113 |
> |
for(int j=0; j<molecules[i].getNTorsions(); j++){ |
114 |
> |
|
115 |
> |
constrained = theArray[j]->is_constrained(); |
116 |
> |
|
117 |
> |
if(constrained){ |
118 |
> |
|
119 |
> |
dummy_plug = theArray[j]->get_constraint(); |
120 |
> |
temp_con[c_n_constrained].set_a( dummy_plug->get_a() ); |
121 |
> |
temp_con[c_n_constrained].set_b( dummy_plug->get_b() ); |
122 |
> |
temp_con[c_n_constrained].set_dsqr( dummy_plug->get_dsqr() ); |
123 |
> |
|
124 |
> |
c_n_constrained++; |
125 |
> |
constrained = 0; |
126 |
> |
} |
127 |
|
} |
128 |
+ |
|
129 |
+ |
|
130 |
|
} |
131 |
|
|
132 |
|
if(c_n_constrained > 0){ |
135 |
|
c_constrained_i = new int[c_n_constrained]; |
136 |
|
c_constrained_j = new int[c_n_constrained]; |
137 |
|
c_constrained_dsqr = new double[c_n_constrained]; |
138 |
< |
|
138 |
> |
|
139 |
|
for( int i = 0; i < c_n_constrained; i++){ |
140 |
|
|
141 |
|
/* add 1 to the index for the fortran arrays. */ |
142 |
< |
|
142 |
> |
|
143 |
|
c_constrained_i[i] = temp_con[i].get_a() + 1; |
144 |
|
c_constrained_j[i] = temp_con[i].get_b() + 1; |
145 |
|
c_constrained_dsqr[i] = temp_con[i].get_dsqr(); |
167 |
|
void Verlet::integrate( void ){ |
168 |
|
|
169 |
|
int i, j; /* loop counters */ |
170 |
< |
int calcPot; |
170 |
> |
int calcPot, calcStress; |
171 |
|
|
172 |
|
double kE; |
173 |
|
|
184 |
|
double *Fz = new double[c_natoms]; |
185 |
|
|
186 |
|
int time; |
187 |
+ |
|
188 |
+ |
double press[9]; |
189 |
|
|
190 |
|
double dt = entry_plug->dt; |
191 |
|
double runTime = entry_plug->run_time; |
204 |
|
DumpWriter* dump_out = new DumpWriter( entry_plug ); |
205 |
|
|
206 |
|
// the first time integrate is called, the forces need to be initialized |
164 |
– |
|
207 |
|
|
208 |
< |
myFF->doForces(1,0); |
208 |
> |
myFF->doForces(1,1); |
209 |
|
|
210 |
|
if( entry_plug->setTemp ){ |
211 |
|
tStats->velocitize(); |
217 |
|
|
218 |
|
calcPot = 0; |
219 |
|
|
220 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT")) { |
221 |
+ |
calcStress = 1; |
222 |
+ |
} else { |
223 |
+ |
calcStress = 0; |
224 |
+ |
} |
225 |
+ |
|
226 |
|
if( c_is_constrained ){ |
227 |
|
for(i = 0; i < n_loops; i++){ |
228 |
|
|
229 |
+ |
|
230 |
|
// fill R, V, and F arrays and RATTLE in fortran |
231 |
< |
|
231 |
> |
|
232 |
|
for( j=0; j<c_natoms; j++ ){ |
233 |
< |
|
233 |
> |
|
234 |
|
Rx[j] = c_atoms[j]->getX(); |
235 |
|
Ry[j] = c_atoms[j]->getY(); |
236 |
|
Rz[j] = c_atoms[j]->getZ(); |
262 |
|
c_atoms[j]->set_vz(Vz[j]); |
263 |
|
} |
264 |
|
|
265 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
266 |
+ |
myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() ); |
267 |
+ |
|
268 |
|
// calculate the forces |
269 |
|
|
270 |
< |
myFF->doForces(calcPot,0); |
270 |
> |
myFF->doForces(calcPot,calcStress); |
271 |
|
|
272 |
|
// finish the constrain move ( same as above. ) |
273 |
|
|
286 |
|
Fz[j] = c_atoms[j]->getFz(); |
287 |
|
} |
288 |
|
|
289 |
+ |
|
290 |
|
v_constrain_b_( dt, c_natoms, c_mass, Rx, Ry, Rz, Vx, Vy, Vz, |
291 |
|
Fx, Fy, Fz, |
292 |
|
kE, c_n_constrained, c_constrained_dsqr, |
304 |
|
c_atoms[j]->set_vz(Vz[j]); |
305 |
|
} |
306 |
|
|
307 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
308 |
+ |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
309 |
+ |
|
310 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
311 |
+ |
tStats->getPressureTensor(press); |
312 |
+ |
myES->NoseHooverAndersonNPT( dt, |
313 |
+ |
tStats->getKinetic(), |
314 |
+ |
press); |
315 |
+ |
} |
316 |
+ |
|
317 |
|
time = i + 1; |
318 |
|
|
319 |
|
if( entry_plug->setTemp ){ |
320 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
321 |
|
} |
322 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
323 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
324 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
323 |
> |
|
324 |
> |
if( !((time+1) % status_n) ) { |
325 |
> |
calcPot = 1; |
326 |
> |
calcStress = 1; |
327 |
> |
} |
328 |
> |
if( !(time % status_n) ){ |
329 |
> |
e_out->writeStat( time * dt ); |
330 |
> |
calcPot = 0; |
331 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
332 |
> |
else calcStress = 0; |
333 |
> |
} |
334 |
|
} |
335 |
|
} |
336 |
|
else{ |
337 |
|
for(i = 0; i < n_loops; i++){ |
338 |
< |
|
338 |
> |
|
339 |
> |
|
340 |
|
move_a( dt ); |
341 |
+ |
|
342 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
343 |
+ |
myES->NoseHooverNVT( dt / 2.0, tStats->getKinetic() ); |
344 |
|
|
345 |
|
// calculate the forces |
346 |
|
|
347 |
< |
myFF->doForces(calcPot,0); |
347 |
> |
myFF->doForces(calcPot,calcStress); |
348 |
|
|
349 |
|
// complete the verlet move |
350 |
|
|
351 |
|
move_b( dt ); |
352 |
|
|
353 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NVT")) |
354 |
+ |
myES->NoseHooverNVT( dt / 2.0 , tStats->getKinetic() ); |
355 |
+ |
|
356 |
+ |
if (!strcasecmp( entry_plug->ensemble, "NPT") ) { |
357 |
+ |
tStats->getPressureTensor(press); |
358 |
+ |
myES->NoseHooverAndersonNPT( dt, |
359 |
+ |
tStats->getKinetic(), |
360 |
+ |
press); |
361 |
+ |
} |
362 |
+ |
|
363 |
|
time = i + 1; |
364 |
|
|
365 |
|
if( entry_plug->setTemp ){ |
366 |
|
if( !(time % vel_n) ) tStats->velocitize(); |
367 |
|
} |
368 |
|
if( !(time % sample_n) ) dump_out->writeDump( time * dt ); |
369 |
< |
if( !((time+1) % status_n) ) calcPot = 1; |
370 |
< |
if( !(time % status_n) ){ e_out->writeStat( time * dt ); calcPot = 0; } |
369 |
> |
if( !((time+1) % status_n) ) { |
370 |
> |
calcPot = 1; |
371 |
> |
calcStress = 1; |
372 |
> |
} |
373 |
> |
if( !(time % status_n) ){ |
374 |
> |
e_out->writeStat( time * dt ); |
375 |
> |
calcPot = 0; |
376 |
> |
if (!strcasecmp(entry_plug->ensemble, "NPT")) calcStress = 1; |
377 |
> |
else calcStress = 0; |
378 |
> |
} |
379 |
|
} |
380 |
|
} |
381 |
|
|