File: | programs/mdrun/md.c |
Location: | line 1378, column 22 |
Description: | The right operand of '>=' is a garbage value |
1 | /* | |||
2 | * This file is part of the GROMACS molecular simulation package. | |||
3 | * | |||
4 | * Copyright (c) 1991-2000, University of Groningen, The Netherlands. | |||
5 | * Copyright (c) 2001-2004, The GROMACS development team. | |||
6 | * Copyright (c) 2011,2012,2013,2014, by the GROMACS development team, led by | |||
7 | * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, | |||
8 | * and including many others, as listed in the AUTHORS file in the | |||
9 | * top-level source directory and at http://www.gromacs.org. | |||
10 | * | |||
11 | * GROMACS is free software; you can redistribute it and/or | |||
12 | * modify it under the terms of the GNU Lesser General Public License | |||
13 | * as published by the Free Software Foundation; either version 2.1 | |||
14 | * of the License, or (at your option) any later version. | |||
15 | * | |||
16 | * GROMACS is distributed in the hope that it will be useful, | |||
17 | * but WITHOUT ANY WARRANTY; without even the implied warranty of | |||
18 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU | |||
19 | * Lesser General Public License for more details. | |||
20 | * | |||
21 | * You should have received a copy of the GNU Lesser General Public | |||
22 | * License along with GROMACS; if not, see | |||
23 | * http://www.gnu.org/licenses, or write to the Free Software Foundation, | |||
24 | * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. | |||
25 | * | |||
26 | * If you want to redistribute modifications to GROMACS, please | |||
27 | * consider that scientific software is very special. Version | |||
28 | * control is crucial - bugs must be traceable. We will be happy to | |||
29 | * consider code for inclusion in the official distribution, but | |||
30 | * derived work must not be called official GROMACS. Details are found | |||
31 | * in the README & COPYING files - if they are missing, get the | |||
32 | * official version at http://www.gromacs.org. | |||
33 | * | |||
34 | * To help us fund GROMACS development, we humbly ask that you cite | |||
35 | * the research papers on the package. Check out http://www.gromacs.org. | |||
36 | */ | |||
37 | #ifdef HAVE_CONFIG_H1 | |||
38 | #include <config.h> | |||
39 | #endif | |||
40 | ||||
41 | #include <stdlib.h> | |||
42 | ||||
43 | #include "typedefs.h" | |||
44 | #include "gromacs/utility/smalloc.h" | |||
45 | #include "gromacs/math/vec.h" | |||
46 | #include "vcm.h" | |||
47 | #include "mdebin.h" | |||
48 | #include "nrnb.h" | |||
49 | #include "calcmu.h" | |||
50 | #include "index.h" | |||
51 | #include "vsite.h" | |||
52 | #include "update.h" | |||
53 | #include "ns.h" | |||
54 | #include "mdrun.h" | |||
55 | #include "md_support.h" | |||
56 | #include "md_logging.h" | |||
57 | #include "network.h" | |||
58 | #include "physics.h" | |||
59 | #include "names.h" | |||
60 | #include "force.h" | |||
61 | #include "disre.h" | |||
62 | #include "orires.h" | |||
63 | #include "pme.h" | |||
64 | #include "mdatoms.h" | |||
65 | #include "repl_ex.h" | |||
66 | #include "deform.h" | |||
67 | #include "qmmm.h" | |||
68 | #include "domdec.h" | |||
69 | #include "domdec_network.h" | |||
70 | #include "gromacs/gmxlib/topsort.h" | |||
71 | #include "coulomb.h" | |||
72 | #include "constr.h" | |||
73 | #include "shellfc.h" | |||
74 | #include "gromacs/gmxpreprocess/compute_io.h" | |||
75 | #include "checkpoint.h" | |||
76 | #include "mtop_util.h" | |||
77 | #include "sighandler.h" | |||
78 | #include "txtdump.h" | |||
79 | #include "gromacs/utility/cstringutil.h" | |||
80 | #include "pme_loadbal.h" | |||
81 | #include "bondf.h" | |||
82 | #include "membed.h" | |||
83 | #include "types/nlistheuristics.h" | |||
84 | #include "types/iteratedconstraints.h" | |||
85 | #include "nbnxn_cuda_data_mgmt.h" | |||
86 | ||||
87 | #include "gromacs/utility/gmxmpi.h" | |||
88 | #include "gromacs/fileio/confio.h" | |||
89 | #include "gromacs/fileio/trajectory_writing.h" | |||
90 | #include "gromacs/fileio/trnio.h" | |||
91 | #include "gromacs/fileio/trxio.h" | |||
92 | #include "gromacs/fileio/xtcio.h" | |||
93 | #include "gromacs/timing/wallcycle.h" | |||
94 | #include "gromacs/timing/walltime_accounting.h" | |||
95 | #include "gromacs/pulling/pull.h" | |||
96 | #include "gromacs/swap/swapcoords.h" | |||
97 | #include "gromacs/imd/imd.h" | |||
98 | ||||
99 | #ifdef GMX_FAHCORE | |||
100 | #include "corewrap.h" | |||
101 | #endif | |||
102 | ||||
103 | static void reset_all_counters(FILE *fplog, t_commrec *cr, | |||
104 | gmx_int64_t step, | |||
105 | gmx_int64_t *step_rel, t_inputrec *ir, | |||
106 | gmx_wallcycle_t wcycle, t_nrnb *nrnb, | |||
107 | gmx_walltime_accounting_t walltime_accounting, | |||
108 | nbnxn_cuda_ptr_t cu_nbv) | |||
109 | { | |||
110 | char sbuf[STEPSTRSIZE22]; | |||
111 | ||||
112 | /* Reset all the counters related to performance over the run */ | |||
113 | md_print_warn(cr, fplog, "step %s: resetting all time and cycle counters\n", | |||
114 | gmx_step_str(step, sbuf)); | |||
115 | ||||
116 | if (cu_nbv) | |||
117 | { | |||
118 | nbnxn_cuda_reset_timings(cu_nbv); | |||
119 | } | |||
120 | ||||
121 | wallcycle_stop(wcycle, ewcRUN); | |||
122 | wallcycle_reset_all(wcycle); | |||
123 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
124 | { | |||
125 | reset_dd_statistics_counters(cr->dd); | |||
126 | } | |||
127 | init_nrnb(nrnb); | |||
128 | ir->init_step += *step_rel; | |||
129 | ir->nsteps -= *step_rel; | |||
130 | *step_rel = 0; | |||
131 | wallcycle_start(wcycle, ewcRUN); | |||
132 | walltime_accounting_start(walltime_accounting); | |||
133 | print_date_and_time(fplog, cr->nodeid, "Restarted time", gmx_gettime()); | |||
134 | } | |||
135 | ||||
136 | double do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[], | |||
137 | const output_env_t oenv, gmx_bool bVerbose, gmx_bool bCompact, | |||
138 | int nstglobalcomm, | |||
139 | gmx_vsite_t *vsite, gmx_constr_t constr, | |||
140 | int stepout, t_inputrec *ir, | |||
141 | gmx_mtop_t *top_global, | |||
142 | t_fcdata *fcd, | |||
143 | t_state *state_global, | |||
144 | t_mdatoms *mdatoms, | |||
145 | t_nrnb *nrnb, gmx_wallcycle_t wcycle, | |||
146 | gmx_edsam_t ed, t_forcerec *fr, | |||
147 | int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t membed, | |||
148 | real cpt_period, real max_hours, | |||
149 | const char gmx_unused__attribute__ ((unused)) *deviceOptions, | |||
150 | int imdport, | |||
151 | unsigned long Flags, | |||
152 | gmx_walltime_accounting_t walltime_accounting) | |||
153 | { | |||
154 | gmx_mdoutf_t outf = NULL((void*)0); | |||
155 | gmx_int64_t step, step_rel; | |||
156 | double elapsed_time; | |||
157 | double t, t0, lam0[efptNR]; | |||
158 | gmx_bool bGStatEveryStep, bGStat, bCalcVir, bCalcEner; | |||
159 | gmx_bool bNS, bNStList, bSimAnn, bStopCM, bRerunMD, bNotLastFrame = FALSE0, | |||
160 | bFirstStep, bStateFromCP, bStateFromTPX, bInitStep, bLastStep, | |||
161 | bBornRadii, bStartingFromCpt; | |||
162 | gmx_bool bDoDHDL = FALSE0, bDoFEP = FALSE0, bDoExpanded = FALSE0; | |||
163 | gmx_bool do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE1, | |||
164 | bForceUpdate = FALSE0, bCPT; | |||
165 | gmx_bool bMasterState; | |||
166 | int force_flags, cglo_flags; | |||
167 | tensor force_vir, shake_vir, total_vir, tmp_vir, pres; | |||
168 | int i, m; | |||
169 | t_trxstatus *status; | |||
170 | rvec mu_tot; | |||
171 | t_vcm *vcm; | |||
172 | t_state *bufstate = NULL((void*)0); | |||
173 | matrix *scale_tot, pcoupl_mu, M, ebox; | |||
174 | gmx_nlheur_t nlh; | |||
175 | t_trxframe rerun_fr; | |||
176 | gmx_repl_ex_t repl_ex = NULL((void*)0); | |||
177 | int nchkpt = 1; | |||
178 | gmx_localtop_t *top; | |||
179 | t_mdebin *mdebin = NULL((void*)0); | |||
180 | t_state *state = NULL((void*)0); | |||
181 | rvec *f_global = NULL((void*)0); | |||
182 | gmx_enerdata_t *enerd; | |||
183 | rvec *f = NULL((void*)0); | |||
184 | gmx_global_stat_t gstat; | |||
185 | gmx_update_t upd = NULL((void*)0); | |||
186 | t_graph *graph = NULL((void*)0); | |||
187 | globsig_t gs; | |||
188 | gmx_groups_t *groups; | |||
189 | gmx_ekindata_t *ekind, *ekind_save; | |||
190 | gmx_shellfc_t shellfc; | |||
191 | int count, nconverged = 0; | |||
192 | real timestep = 0; | |||
193 | double tcount = 0; | |||
194 | gmx_bool bConverged = TRUE1, bOK, bSumEkinhOld, bExchanged, bNeedRepartition; | |||
195 | gmx_bool bAppend; | |||
196 | gmx_bool bResetCountersHalfMaxH = FALSE0; | |||
197 | gmx_bool bVV, bIterativeCase, bFirstIterate, bTemp, bPres, bTrotter; | |||
198 | gmx_bool bUpdateDoLR; | |||
199 | real dvdl_constr; | |||
200 | rvec *cbuf = NULL((void*)0); | |||
201 | matrix lastbox; | |||
202 | real veta_save, scalevir, tracevir; | |||
203 | real vetanew = 0; | |||
204 | int lamnew = 0; | |||
205 | /* for FEP */ | |||
206 | int nstfep; | |||
207 | double cycles; | |||
208 | real saved_conserved_quantity = 0; | |||
209 | real last_ekin = 0; | |||
210 | int iter_i; | |||
211 | t_extmass MassQ; | |||
212 | int **trotter_seq; | |||
213 | char sbuf[STEPSTRSIZE22], sbuf2[STEPSTRSIZE22]; | |||
214 | int handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/ | |||
215 | gmx_iterate_t iterate; | |||
216 | gmx_int64_t multisim_nsteps = -1; /* number of steps to do before first multisim | |||
217 | simulation stops. If equal to zero, don't | |||
218 | communicate any more between multisims.*/ | |||
219 | /* PME load balancing data for GPU kernels */ | |||
220 | pme_load_balancing_t pme_loadbal = NULL((void*)0); | |||
221 | double cycles_pmes; | |||
222 | gmx_bool bPMETuneTry = FALSE0, bPMETuneRunning = FALSE0; | |||
223 | ||||
224 | /* Interactive MD */ | |||
225 | gmx_bool bIMDstep = FALSE0; | |||
226 | ||||
227 | #ifdef GMX_FAHCORE | |||
228 | /* Temporary addition for FAHCORE checkpointing */ | |||
229 | int chkpt_ret; | |||
230 | #endif | |||
231 | ||||
232 | /* Check for special mdrun options */ | |||
233 | bRerunMD = (Flags & MD_RERUN(1<<4)); | |||
234 | bAppend = (Flags & MD_APPENDFILES(1<<15)); | |||
235 | if (Flags & MD_RESETCOUNTERSHALFWAY(1<<19)) | |||
| ||||
236 | { | |||
237 | if (ir->nsteps > 0) | |||
238 | { | |||
239 | /* Signal to reset the counters half the simulation steps. */ | |||
240 | wcycle_set_reset_counters(wcycle, ir->nsteps/2); | |||
241 | } | |||
242 | /* Signal to reset the counters halfway the simulation time. */ | |||
243 | bResetCountersHalfMaxH = (max_hours > 0); | |||
244 | } | |||
245 | ||||
246 | /* md-vv uses averaged full step velocities for T-control | |||
247 | md-vv-avek uses averaged half step velocities for T-control (but full step ekin for P control) | |||
248 | md uses averaged half step kinetic energies to determine temperature unless defined otherwise by GMX_EKIN_AVE_VEL; */ | |||
249 | bVV = EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK); | |||
250 | if (bVV) /* to store the initial velocities while computing virial */ | |||
251 | { | |||
252 | snew(cbuf, top_global->natoms)(cbuf) = save_calloc("cbuf", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 252, (top_global->natoms), sizeof(*(cbuf))); | |||
253 | } | |||
254 | /* all the iteratative cases - only if there are constraints */ | |||
255 | bIterativeCase = ((IR_NPH_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && (!(((ir)->etc == etcNOSEHOOVER ))))) || IR_NPT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && ((ir)->etc == etcNOSEHOOVER )))) && (constr) && (!bRerunMD)); | |||
256 | gmx_iterate_init(&iterate, FALSE0); /* The default value of iterate->bIterationActive is set to | |||
257 | false in this step. The correct value, true or false, | |||
258 | is set at each step, as it depends on the frequency of temperature | |||
259 | and pressure control.*/ | |||
260 | bTrotter = (bVV && (IR_NPT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && ((ir)->etc == etcNOSEHOOVER ))) || IR_NPH_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && (((ir)->epc == epcMTTK) && (!(((ir)->etc == etcNOSEHOOVER ))))) || IR_NVT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER ))))); | |||
261 | ||||
262 | if (bRerunMD) | |||
263 | { | |||
264 | /* Since we don't know if the frames read are related in any way, | |||
265 | * rebuild the neighborlist at every step. | |||
266 | */ | |||
267 | ir->nstlist = 1; | |||
268 | ir->nstcalcenergy = 1; | |||
269 | nstglobalcomm = 1; | |||
270 | } | |||
271 | ||||
272 | check_ir_old_tpx_versions(cr, fplog, ir, top_global); | |||
273 | ||||
274 | nstglobalcomm = check_nstglobalcomm(fplog, cr, nstglobalcomm, ir); | |||
275 | bGStatEveryStep = (nstglobalcomm == 1); | |||
276 | ||||
277 | if (!bGStatEveryStep && ir->nstlist == -1 && fplog != NULL((void*)0)) | |||
278 | { | |||
279 | fprintf(fplog, | |||
280 | "To reduce the energy communication with nstlist = -1\n" | |||
281 | "the neighbor list validity should not be checked at every step,\n" | |||
282 | "this means that exact integration is not guaranteed.\n" | |||
283 | "The neighbor list validity is checked after:\n" | |||
284 | " <n.list life time> - 2*std.dev.(n.list life time) steps.\n" | |||
285 | "In most cases this will result in exact integration.\n" | |||
286 | "This reduces the energy communication by a factor of 2 to 3.\n" | |||
287 | "If you want less energy communication, set nstlist > 3.\n\n"); | |||
288 | } | |||
289 | ||||
290 | if (bRerunMD) | |||
291 | { | |||
292 | ir->nstxout_compressed = 0; | |||
293 | } | |||
294 | groups = &top_global->groups; | |||
295 | ||||
296 | /* Initial values */ | |||
297 | init_md(fplog, cr, ir, oenv, &t, &t0, state_global->lambda, | |||
298 | &(state_global->fep_state), lam0, | |||
299 | nrnb, top_global, &upd, | |||
300 | nfile, fnm, &outf, &mdebin, | |||
301 | force_vir, shake_vir, mu_tot, &bSimAnn, &vcm, Flags); | |||
302 | ||||
303 | clear_mat(total_vir); | |||
304 | clear_mat(pres); | |||
305 | /* Energy terms and groups */ | |||
306 | snew(enerd, 1)(enerd) = save_calloc("enerd", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 306, (1), sizeof(*(enerd))); | |||
307 | init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda, | |||
308 | enerd); | |||
309 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
310 | { | |||
311 | f = NULL((void*)0); | |||
312 | } | |||
313 | else | |||
314 | { | |||
315 | snew(f, top_global->natoms)(f) = save_calloc("f", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 315, (top_global->natoms), sizeof(*(f))); | |||
316 | } | |||
317 | ||||
318 | /* Kinetic energy data */ | |||
319 | snew(ekind, 1)(ekind) = save_calloc("ekind", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 319, (1), sizeof(*(ekind))); | |||
320 | init_ekindata(fplog, top_global, &(ir->opts), ekind); | |||
321 | /* needed for iteration of constraints */ | |||
322 | snew(ekind_save, 1)(ekind_save) = save_calloc("ekind_save", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 322, (1), sizeof(*(ekind_save))); | |||
323 | init_ekindata(fplog, top_global, &(ir->opts), ekind_save); | |||
324 | /* Copy the cos acceleration to the groups struct */ | |||
325 | ekind->cosacc.cos_accel = ir->cos_accel; | |||
326 | ||||
327 | gstat = global_stat_init(ir); | |||
328 | debug_gmx(); | |||
329 | ||||
330 | /* Check for polarizable models and flexible constraints */ | |||
331 | shellfc = init_shell_flexcon(fplog, fr->cutoff_scheme == ecutsVERLET, | |||
332 | top_global, n_flexible_constraints(constr), | |||
333 | (ir->bContinuation || | |||
334 | (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) && !MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)))) ? | |||
335 | NULL((void*)0) : state_global->x); | |||
336 | ||||
337 | if (DEFORM(*ir)((*ir).deform[0][0] != 0 || (*ir).deform[1][1] != 0 || (*ir). deform[2][2] != 0 || (*ir).deform[1][0] != 0 || (*ir).deform[ 2][0] != 0 || (*ir).deform[2][1] != 0)) | |||
338 | { | |||
339 | tMPI_Thread_mutex_lock(&deform_init_box_mutex); | |||
340 | set_deform_reference_box(upd, | |||
341 | deform_init_init_step_tpx, | |||
342 | deform_init_box_tpx); | |||
343 | tMPI_Thread_mutex_unlock(&deform_init_box_mutex); | |||
344 | } | |||
345 | ||||
346 | { | |||
347 | double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1); | |||
348 | if ((io > 2000) && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
349 | { | |||
350 | fprintf(stderrstderr, | |||
351 | "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", | |||
352 | io); | |||
353 | } | |||
354 | } | |||
355 | ||||
356 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
357 | { | |||
358 | top = dd_init_local_top(top_global); | |||
359 | ||||
360 | snew(state, 1)(state) = save_calloc("state", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 360, (1), sizeof(*(state))); | |||
361 | dd_init_local_state(cr->dd, state_global, state); | |||
362 | ||||
363 | if (DDMASTER(cr->dd)((cr->dd)->rank == (cr->dd)->masterrank) && ir->nstfout) | |||
364 | { | |||
365 | snew(f_global, state_global->natoms)(f_global) = save_calloc("f_global", "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c" , 365, (state_global->natoms), sizeof(*(f_global))); | |||
366 | } | |||
367 | } | |||
368 | else | |||
369 | { | |||
370 | top = gmx_mtop_generate_local_top(top_global, ir); | |||
371 | ||||
372 | forcerec_set_excl_load(fr, top); | |||
373 | ||||
374 | state = serial_init_local_state(state_global); | |||
375 | f_global = f; | |||
376 | ||||
377 | atoms2md(top_global, ir, 0, NULL((void*)0), top_global->natoms, mdatoms); | |||
378 | ||||
379 | if (vsite) | |||
380 | { | |||
381 | set_vsite_top(vsite, top, mdatoms, cr); | |||
382 | } | |||
383 | ||||
384 | if (ir->ePBC != epbcNONE && !fr->bMolPBC) | |||
385 | { | |||
386 | graph = mk_graph(fplog, &(top->idef), 0, top_global->natoms, FALSE0, FALSE0); | |||
387 | } | |||
388 | ||||
389 | if (shellfc) | |||
390 | { | |||
391 | make_local_shells(cr, mdatoms, shellfc); | |||
392 | } | |||
393 | ||||
394 | setup_bonded_threading(fr, &top->idef); | |||
395 | } | |||
396 | ||||
397 | /* Set up interactive MD (IMD) */ | |||
398 | init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, state_global->x, | |||
399 | nfile, fnm, oenv, imdport, Flags); | |||
400 | ||||
401 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
402 | { | |||
403 | /* Distribute the charge groups over the nodes from the master node */ | |||
404 | dd_partition_system(fplog, ir->init_step, cr, TRUE1, 1, | |||
405 | state_global, top_global, ir, | |||
406 | state, &f, mdatoms, top, fr, | |||
407 | vsite, shellfc, constr, | |||
408 | nrnb, wcycle, FALSE0); | |||
409 | ||||
410 | } | |||
411 | ||||
412 | update_mdatoms(mdatoms, state->lambda[efptMASS]); | |||
413 | ||||
414 | if (opt2bSet("-cpi", nfile, fnm)) | |||
415 | { | |||
416 | bStateFromCP = gmx_fexist_master(opt2fn_master("-cpi", nfile, fnm, cr), cr); | |||
417 | } | |||
418 | else | |||
419 | { | |||
420 | bStateFromCP = FALSE0; | |||
421 | } | |||
422 | ||||
423 | if (ir->bExpanded) | |||
424 | { | |||
425 | init_expanded_ensemble(bStateFromCP, ir, &state->dfhist); | |||
426 | } | |||
427 | ||||
428 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
429 | { | |||
430 | if (bStateFromCP) | |||
431 | { | |||
432 | /* Update mdebin with energy history if appending to output files */ | |||
433 | if (Flags & MD_APPENDFILES(1<<15)) | |||
434 | { | |||
435 | restore_energyhistory_from_state(mdebin, &state_global->enerhist); | |||
436 | } | |||
437 | else | |||
438 | { | |||
439 | /* We might have read an energy history from checkpoint, | |||
440 | * free the allocated memory and reset the counts. | |||
441 | */ | |||
442 | done_energyhistory(&state_global->enerhist); | |||
443 | init_energyhistory(&state_global->enerhist); | |||
444 | } | |||
445 | } | |||
446 | /* Set the initial energy history in state by updating once */ | |||
447 | update_energyhistory(&state_global->enerhist, mdebin); | |||
448 | } | |||
449 | ||||
450 | /* Initialize constraints */ | |||
451 | if (constr && !DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
452 | { | |||
453 | set_constraints(constr, top, ir, mdatoms, cr); | |||
454 | } | |||
455 | ||||
456 | if (repl_ex_nst > 0) | |||
457 | { | |||
458 | /* We need to be sure replica exchange can only occur | |||
459 | * when the energies are current */ | |||
460 | check_nst_param(fplog, cr, "nstcalcenergy", ir->nstcalcenergy, | |||
461 | "repl_ex_nst", &repl_ex_nst); | |||
462 | /* This check needs to happen before inter-simulation | |||
463 | * signals are initialized, too */ | |||
464 | } | |||
465 | if (repl_ex_nst > 0 && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
466 | { | |||
467 | repl_ex = init_replica_exchange(fplog, cr->ms, state_global, ir, | |||
468 | repl_ex_nst, repl_ex_nex, repl_ex_seed); | |||
469 | } | |||
470 | ||||
471 | /* PME tuning is only supported with GPUs or PME nodes and not with rerun. | |||
472 | * PME tuning is not supported with PME only for LJ and not for Coulomb. | |||
473 | */ | |||
474 | if ((Flags & MD_TUNEPME(1<<20)) && | |||
475 | EEL_PME(fr->eeltype)((fr->eeltype) == eelPME || (fr->eeltype) == eelPMESWITCH || (fr->eeltype) == eelPMEUSER || (fr->eeltype) == eelPMEUSERSWITCH || (fr->eeltype) == eelP3M_AD) && | |||
476 | ( (fr->cutoff_scheme == ecutsVERLET && fr->nbv->bUseGPU) || !(cr->duty & DUTY_PME(1<<1))) && | |||
477 | !bRerunMD) | |||
478 | { | |||
479 | pme_loadbal_init(&pme_loadbal, ir, state->box, fr->ic, fr->pmedata); | |||
480 | cycles_pmes = 0; | |||
481 | if (cr->duty & DUTY_PME(1<<1)) | |||
482 | { | |||
483 | /* Start tuning right away, as we can't measure the load */ | |||
484 | bPMETuneRunning = TRUE1; | |||
485 | } | |||
486 | else | |||
487 | { | |||
488 | /* Separate PME nodes, we can measure the PP/PME load balance */ | |||
489 | bPMETuneTry = TRUE1; | |||
490 | } | |||
491 | } | |||
492 | ||||
493 | if (!ir->bContinuation && !bRerunMD) | |||
494 | { | |||
495 | if (mdatoms->cFREEZE && (state->flags & (1<<estV))) | |||
496 | { | |||
497 | /* Set the velocities of frozen particles to zero */ | |||
498 | for (i = 0; i < mdatoms->homenr; i++) | |||
499 | { | |||
500 | for (m = 0; m < DIM3; m++) | |||
501 | { | |||
502 | if (ir->opts.nFreeze[mdatoms->cFREEZE[i]][m]) | |||
503 | { | |||
504 | state->v[i][m] = 0; | |||
505 | } | |||
506 | } | |||
507 | } | |||
508 | } | |||
509 | ||||
510 | if (constr) | |||
511 | { | |||
512 | /* Constrain the initial coordinates and velocities */ | |||
513 | do_constrain_first(fplog, constr, ir, mdatoms, state, | |||
514 | cr, nrnb, fr, top); | |||
515 | } | |||
516 | if (vsite) | |||
517 | { | |||
518 | /* Construct the virtual sites for the initial configuration */ | |||
519 | construct_vsites(vsite, state->x, ir->delta_t, NULL((void*)0), | |||
520 | top->idef.iparams, top->idef.il, | |||
521 | fr->ePBC, fr->bMolPBC, cr, state->box); | |||
522 | } | |||
523 | } | |||
524 | ||||
525 | debug_gmx(); | |||
526 | ||||
527 | /* set free energy calculation frequency as the minimum | |||
528 | greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/ | |||
529 | nstfep = ir->fepvals->nstdhdl; | |||
530 | if (ir->bExpanded) | |||
531 | { | |||
532 | nstfep = gmx_greatest_common_divisor(ir->fepvals->nstdhdl, nstfep); | |||
533 | } | |||
534 | if (repl_ex_nst > 0) | |||
535 | { | |||
536 | nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep); | |||
537 | } | |||
538 | ||||
539 | /* I'm assuming we need global communication the first time! MRS */ | |||
540 | cglo_flags = (CGLO_TEMPERATURE(1<<7) | CGLO_GSTAT(1<<4) | |||
541 | | ((ir->comm_mode != ecmNO) ? CGLO_STOPCM(1<<3) : 0) | |||
542 | | (bVV ? CGLO_PRESSURE(1<<8) : 0) | |||
543 | | (bVV ? CGLO_CONSTRAINT(1<<9) : 0) | |||
544 | | (bRerunMD ? CGLO_RERUNMD(1<<1) : 0) | |||
545 | | ((Flags & MD_READ_EKIN(1<<17)) ? CGLO_READEKIN(1<<12) : 0)); | |||
546 | ||||
547 | bSumEkinhOld = FALSE0; | |||
548 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
549 | NULL((void*)0), enerd, force_vir, shake_vir, total_vir, pres, mu_tot, | |||
550 | constr, NULL((void*)0), FALSE0, state->box, | |||
551 | top_global, &bSumEkinhOld, cglo_flags); | |||
552 | if (ir->eI == eiVVAK) | |||
553 | { | |||
554 | /* a second call to get the half step temperature initialized as well */ | |||
555 | /* we do the same call as above, but turn the pressure off -- internally to | |||
556 | compute_globals, this is recognized as a velocity verlet half-step | |||
557 | kinetic energy calculation. This minimized excess variables, but | |||
558 | perhaps loses some logic?*/ | |||
559 | ||||
560 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
561 | NULL((void*)0), enerd, force_vir, shake_vir, total_vir, pres, mu_tot, | |||
562 | constr, NULL((void*)0), FALSE0, state->box, | |||
563 | top_global, &bSumEkinhOld, | |||
564 | cglo_flags &~(CGLO_STOPCM(1<<3) | CGLO_PRESSURE(1<<8))); | |||
565 | } | |||
566 | ||||
567 | /* Calculate the initial half step temperature, and save the ekinh_old */ | |||
568 | if (!(Flags & MD_STARTFROMCPT(1<<18))) | |||
569 | { | |||
570 | for (i = 0; (i < ir->opts.ngtc); i++) | |||
571 | { | |||
572 | copy_mat(ekind->tcstat[i].ekinh, ekind->tcstat[i].ekinh_old); | |||
573 | } | |||
574 | } | |||
575 | if (ir->eI != eiVV) | |||
576 | { | |||
577 | enerd->term[F_TEMP] *= 2; /* result of averages being done over previous and current step, | |||
578 | and there is no previous step */ | |||
579 | } | |||
580 | ||||
581 | /* if using an iterative algorithm, we need to create a working directory for the state. */ | |||
582 | if (bIterativeCase) | |||
583 | { | |||
584 | bufstate = init_bufstate(state); | |||
585 | } | |||
586 | ||||
587 | /* need to make an initiation call to get the Trotter variables set, as well as other constants for non-trotter | |||
588 | temperature control */ | |||
589 | trotter_seq = init_npt_vars(ir, state, &MassQ, bTrotter); | |||
590 | ||||
591 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
592 | { | |||
593 | if (constr && !ir->bContinuation && ir->eConstrAlg == econtLINCS) | |||
594 | { | |||
595 | fprintf(fplog, | |||
596 | "RMS relative constraint deviation after constraining: %.2e\n", | |||
597 | constr_rmsd(constr, FALSE0)); | |||
598 | } | |||
599 | if (EI_STATE_VELOCITY(ir->eI)(((ir->eI) == eiMD || ((ir->eI) == eiVV || (ir->eI) == eiVVAK)) || ((ir->eI) == eiSD1 || (ir->eI) == eiSD2))) | |||
600 | { | |||
601 | fprintf(fplog, "Initial temperature: %g K\n", enerd->term[F_TEMP]); | |||
602 | } | |||
603 | if (bRerunMD) | |||
604 | { | |||
605 | fprintf(stderrstderr, "starting md rerun '%s', reading coordinates from" | |||
606 | " input trajectory '%s'\n\n", | |||
607 | *(top_global->name), opt2fn("-rerun", nfile, fnm)); | |||
608 | if (bVerbose) | |||
609 | { | |||
610 | fprintf(stderrstderr, "Calculated time to finish depends on nsteps from " | |||
611 | "run input file,\nwhich may not correspond to the time " | |||
612 | "needed to process input trajectory.\n\n"); | |||
613 | } | |||
614 | } | |||
615 | else | |||
616 | { | |||
617 | char tbuf[20]; | |||
618 | fprintf(stderrstderr, "starting mdrun '%s'\n", | |||
619 | *(top_global->name)); | |||
620 | if (ir->nsteps >= 0) | |||
621 | { | |||
622 | sprintf(tbuf, "%8.1f", (ir->init_step+ir->nsteps)*ir->delta_t); | |||
623 | } | |||
624 | else | |||
625 | { | |||
626 | sprintf(tbuf, "%s", "infinite"); | |||
627 | } | |||
628 | if (ir->init_step > 0) | |||
629 | { | |||
630 | fprintf(stderrstderr, "%s steps, %s ps (continuing from step %s, %8.1f ps).\n", | |||
631 | gmx_step_str(ir->init_step+ir->nsteps, sbuf), tbuf, | |||
632 | gmx_step_str(ir->init_step, sbuf2), | |||
633 | ir->init_step*ir->delta_t); | |||
634 | } | |||
635 | else | |||
636 | { | |||
637 | fprintf(stderrstderr, "%s steps, %s ps.\n", | |||
638 | gmx_step_str(ir->nsteps, sbuf), tbuf); | |||
639 | } | |||
640 | } | |||
641 | fprintf(fplog, "\n"); | |||
642 | } | |||
643 | ||||
644 | walltime_accounting_start(walltime_accounting); | |||
645 | wallcycle_start(wcycle, ewcRUN); | |||
646 | print_start(fplog, cr, walltime_accounting, "mdrun"); | |||
647 | ||||
648 | /* safest point to do file checkpointing is here. More general point would be immediately before integrator call */ | |||
649 | #ifdef GMX_FAHCORE | |||
650 | chkpt_ret = fcCheckPointParallel( cr->nodeid, | |||
651 | NULL((void*)0), 0); | |||
652 | if (chkpt_ret == 0) | |||
653 | { | |||
654 | gmx_fatal( 3, __FILE__"/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", __LINE__654, "Checkpoint error on step %d\n", 0 ); | |||
655 | } | |||
656 | #endif | |||
657 | ||||
658 | debug_gmx(); | |||
659 | /*********************************************************** | |||
660 | * | |||
661 | * Loop over MD steps | |||
662 | * | |||
663 | ************************************************************/ | |||
664 | ||||
665 | /* if rerunMD then read coordinates and velocities from input trajectory */ | |||
666 | if (bRerunMD) | |||
667 | { | |||
668 | if (getenv("GMX_FORCE_UPDATE")) | |||
669 | { | |||
670 | bForceUpdate = TRUE1; | |||
671 | } | |||
672 | ||||
673 | rerun_fr.natoms = 0; | |||
674 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
675 | { | |||
676 | bNotLastFrame = read_first_frame(oenv, &status, | |||
677 | opt2fn("-rerun", nfile, fnm), | |||
678 | &rerun_fr, TRX_NEED_X(1<<1) | TRX_READ_V(1<<2)); | |||
679 | if (rerun_fr.natoms != top_global->natoms) | |||
680 | { | |||
681 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 681, | |||
682 | "Number of atoms in trajectory (%d) does not match the " | |||
683 | "run input file (%d)\n", | |||
684 | rerun_fr.natoms, top_global->natoms); | |||
685 | } | |||
686 | if (ir->ePBC != epbcNONE) | |||
687 | { | |||
688 | if (!rerun_fr.bBox) | |||
689 | { | |||
690 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 690, "Rerun trajectory frame step %d time %f does not contain a box, while pbc is used", rerun_fr.step, rerun_fr.time); | |||
691 | } | |||
692 | if (max_cutoff2(ir->ePBC, rerun_fr.box) < sqr(fr->rlistlong)) | |||
693 | { | |||
694 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 694, "Rerun trajectory frame step %d time %f has too small box dimensions", rerun_fr.step, rerun_fr.time); | |||
695 | } | |||
696 | } | |||
697 | } | |||
698 | ||||
699 | if (PAR(cr)((cr)->nnodes > 1)) | |||
700 | { | |||
701 | rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame); | |||
702 | } | |||
703 | ||||
704 | if (ir->ePBC != epbcNONE) | |||
705 | { | |||
706 | /* Set the shift vectors. | |||
707 | * Necessary here when have a static box different from the tpr box. | |||
708 | */ | |||
709 | calc_shifts(rerun_fr.box, fr->shift_vec); | |||
710 | } | |||
711 | } | |||
712 | ||||
713 | /* loop over MD steps or if rerunMD to end of input trajectory */ | |||
714 | bFirstStep = TRUE1; | |||
715 | /* Skip the first Nose-Hoover integration when we get the state from tpx */ | |||
716 | bStateFromTPX = !bStateFromCP; | |||
717 | bInitStep = bFirstStep && (bStateFromTPX || bVV); | |||
718 | bStartingFromCpt = (Flags & MD_STARTFROMCPT(1<<18)) && bInitStep; | |||
719 | bLastStep = FALSE0; | |||
720 | bSumEkinhOld = FALSE0; | |||
721 | bExchanged = FALSE0; | |||
722 | bNeedRepartition = FALSE0; | |||
723 | ||||
724 | init_global_signals(&gs, cr, ir, repl_ex_nst); | |||
725 | ||||
726 | step = ir->init_step; | |||
727 | step_rel = 0; | |||
728 | ||||
729 | if (ir->nstlist == -1) | |||
730 | { | |||
731 | init_nlistheuristics(&nlh, bGStatEveryStep, step); | |||
732 | } | |||
733 | ||||
734 | if (MULTISIM(cr)((cr)->ms) && (repl_ex_nst <= 0 )) | |||
735 | { | |||
736 | /* check how many steps are left in other sims */ | |||
737 | multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps); | |||
738 | } | |||
739 | ||||
740 | ||||
741 | /* and stop now if we should */ | |||
742 | bLastStep = (bRerunMD || (ir->nsteps >= 0 && step_rel > ir->nsteps) || | |||
743 | ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps ))); | |||
744 | while (!bLastStep || (bRerunMD && bNotLastFrame)) | |||
745 | { | |||
746 | ||||
747 | wallcycle_start(wcycle, ewcSTEP); | |||
748 | ||||
749 | if (bRerunMD) | |||
750 | { | |||
751 | if (rerun_fr.bStep) | |||
752 | { | |||
753 | step = rerun_fr.step; | |||
754 | step_rel = step - ir->init_step; | |||
755 | } | |||
756 | if (rerun_fr.bTime) | |||
757 | { | |||
758 | t = rerun_fr.time; | |||
759 | } | |||
760 | else | |||
761 | { | |||
762 | t = step; | |||
763 | } | |||
764 | } | |||
765 | else | |||
766 | { | |||
767 | bLastStep = (step_rel == ir->nsteps); | |||
768 | t = t0 + step*ir->delta_t; | |||
769 | } | |||
770 | ||||
771 | if (ir->efep != efepNO || ir->bSimTemp) | |||
772 | { | |||
773 | /* find and set the current lambdas. If rerunning, we either read in a state, or a lambda value, | |||
774 | requiring different logic. */ | |||
775 | ||||
776 | set_current_lambdas(step, ir->fepvals, bRerunMD, &rerun_fr, state_global, state, lam0); | |||
777 | bDoDHDL = do_per_step(step, ir->fepvals->nstdhdl); | |||
778 | bDoFEP = (do_per_step(step, nstfep) && (ir->efep != efepNO)); | |||
779 | bDoExpanded = (do_per_step(step, ir->expandedvals->nstexpanded) | |||
780 | && (ir->bExpanded) && (step > 0) && (!bStartingFromCpt)); | |||
781 | } | |||
782 | ||||
783 | if (bSimAnn) | |||
784 | { | |||
785 | update_annealing_target_temp(&(ir->opts), t); | |||
786 | } | |||
787 | ||||
788 | if (bRerunMD) | |||
789 | { | |||
790 | if (!DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) || MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
791 | { | |||
792 | for (i = 0; i < state_global->natoms; i++) | |||
793 | { | |||
794 | copy_rvec(rerun_fr.x[i], state_global->x[i]); | |||
795 | } | |||
796 | if (rerun_fr.bV) | |||
797 | { | |||
798 | for (i = 0; i < state_global->natoms; i++) | |||
799 | { | |||
800 | copy_rvec(rerun_fr.v[i], state_global->v[i]); | |||
801 | } | |||
802 | } | |||
803 | else | |||
804 | { | |||
805 | for (i = 0; i < state_global->natoms; i++) | |||
806 | { | |||
807 | clear_rvec(state_global->v[i]); | |||
808 | } | |||
809 | if (bRerunWarnNoV) | |||
810 | { | |||
811 | fprintf(stderrstderr, "\nWARNING: Some frames do not contain velocities.\n" | |||
812 | " Ekin, temperature and pressure are incorrect,\n" | |||
813 | " the virial will be incorrect when constraints are present.\n" | |||
814 | "\n"); | |||
815 | bRerunWarnNoV = FALSE0; | |||
816 | } | |||
817 | } | |||
818 | } | |||
819 | copy_mat(rerun_fr.box, state_global->box); | |||
820 | copy_mat(state_global->box, state->box); | |||
821 | ||||
822 | if (vsite && (Flags & MD_RERUN_VSITE(1<<5))) | |||
823 | { | |||
824 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
825 | { | |||
826 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 826, "Vsite recalculation with -rerun is not implemented with domain decomposition, use a single rank"); | |||
827 | } | |||
828 | if (graph) | |||
829 | { | |||
830 | /* Following is necessary because the graph may get out of sync | |||
831 | * with the coordinates if we only have every N'th coordinate set | |||
832 | */ | |||
833 | mk_mshift(fplog, graph, fr->ePBC, state->box, state->x); | |||
834 | shift_self(graph, state->box, state->x); | |||
835 | } | |||
836 | construct_vsites(vsite, state->x, ir->delta_t, state->v, | |||
837 | top->idef.iparams, top->idef.il, | |||
838 | fr->ePBC, fr->bMolPBC, cr, state->box); | |||
839 | if (graph) | |||
840 | { | |||
841 | unshift_self(graph, state->box, state->x); | |||
842 | } | |||
843 | } | |||
844 | } | |||
845 | ||||
846 | /* Stop Center of Mass motion */ | |||
847 | bStopCM = (ir->comm_mode != ecmNO && do_per_step(step, ir->nstcomm)); | |||
848 | ||||
849 | if (bRerunMD) | |||
850 | { | |||
851 | /* for rerun MD always do Neighbour Searching */ | |||
852 | bNS = (bFirstStep || ir->nstlist != 0); | |||
853 | bNStList = bNS; | |||
854 | } | |||
855 | else | |||
856 | { | |||
857 | /* Determine whether or not to do Neighbour Searching and LR */ | |||
858 | bNStList = (ir->nstlist > 0 && step % ir->nstlist == 0); | |||
859 | ||||
860 | bNS = (bFirstStep || bExchanged || bNeedRepartition || bNStList || bDoFEP || | |||
861 | (ir->nstlist == -1 && nlh.nabnsb > 0)); | |||
862 | ||||
863 | if (bNS && ir->nstlist == -1) | |||
864 | { | |||
865 | set_nlistheuristics(&nlh, bFirstStep || bExchanged || bNeedRepartition || bDoFEP, step); | |||
866 | } | |||
867 | } | |||
868 | ||||
869 | /* check whether we should stop because another simulation has | |||
870 | stopped. */ | |||
871 | if (MULTISIM(cr)((cr)->ms)) | |||
872 | { | |||
873 | if ( (multisim_nsteps >= 0) && (step_rel >= multisim_nsteps) && | |||
874 | (multisim_nsteps != ir->nsteps) ) | |||
875 | { | |||
876 | if (bNS) | |||
877 | { | |||
878 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
879 | { | |||
880 | fprintf(stderrstderr, | |||
881 | "Stopping simulation %d because another one has finished\n", | |||
882 | cr->ms->sim); | |||
883 | } | |||
884 | bLastStep = TRUE1; | |||
885 | gs.sig[eglsCHKPT] = 1; | |||
886 | } | |||
887 | } | |||
888 | } | |||
889 | ||||
890 | /* < 0 means stop at next step, > 0 means stop at next NS step */ | |||
891 | if ( (gs.set[eglsSTOPCOND] < 0) || | |||
892 | ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) ) | |||
893 | { | |||
894 | bLastStep = TRUE1; | |||
895 | } | |||
896 | ||||
897 | /* Determine whether or not to update the Born radii if doing GB */ | |||
898 | bBornRadii = bFirstStep; | |||
899 | if (ir->implicit_solvent && (step % ir->nstgbradii == 0)) | |||
900 | { | |||
901 | bBornRadii = TRUE1; | |||
902 | } | |||
903 | ||||
904 | do_log = do_per_step(step, ir->nstlog) || bFirstStep || bLastStep; | |||
905 | do_verbose = bVerbose && | |||
906 | (step % stepout == 0 || bFirstStep || bLastStep); | |||
907 | ||||
908 | if (bNS && !(bFirstStep && ir->bContinuation && !bRerunMD)) | |||
909 | { | |||
910 | if (bRerunMD) | |||
911 | { | |||
912 | bMasterState = TRUE1; | |||
913 | } | |||
914 | else | |||
915 | { | |||
916 | bMasterState = FALSE0; | |||
917 | /* Correct the new box if it is too skewed */ | |||
918 | if (DYNAMIC_BOX(*ir)((*ir).epc != epcNO || (*ir).eI == eiTPI || ((*ir).deform[0][ 0] != 0 || (*ir).deform[1][1] != 0 || (*ir).deform[2][2] != 0 || (*ir).deform[1][0] != 0 || (*ir).deform[2][0] != 0 || (*ir ).deform[2][1] != 0))) | |||
919 | { | |||
920 | if (correct_box(fplog, step, state->box, graph)) | |||
921 | { | |||
922 | bMasterState = TRUE1; | |||
923 | } | |||
924 | } | |||
925 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) && bMasterState) | |||
926 | { | |||
927 | dd_collect_state(cr->dd, state, state_global); | |||
928 | } | |||
929 | } | |||
930 | ||||
931 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
932 | { | |||
933 | /* Repartition the domain decomposition */ | |||
934 | wallcycle_start(wcycle, ewcDOMDEC); | |||
935 | dd_partition_system(fplog, step, cr, | |||
936 | bMasterState, nstglobalcomm, | |||
937 | state_global, top_global, ir, | |||
938 | state, &f, mdatoms, top, fr, | |||
939 | vsite, shellfc, constr, | |||
940 | nrnb, wcycle, | |||
941 | do_verbose && !bPMETuneRunning); | |||
942 | wallcycle_stop(wcycle, ewcDOMDEC); | |||
943 | /* If using an iterative integrator, reallocate space to match the decomposition */ | |||
944 | } | |||
945 | } | |||
946 | ||||
947 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && do_log) | |||
948 | { | |||
949 | print_ebin_header(fplog, step, t, state->lambda[efptFEP]); /* can we improve the information printed here? */ | |||
950 | } | |||
951 | ||||
952 | if (ir->efep != efepNO) | |||
953 | { | |||
954 | update_mdatoms(mdatoms, state->lambda[efptMASS]); | |||
955 | } | |||
956 | ||||
957 | if ((bRerunMD && rerun_fr.bV) || bExchanged) | |||
958 | { | |||
959 | ||||
960 | /* We need the kinetic energy at minus the half step for determining | |||
961 | * the full step kinetic energy and possibly for T-coupling.*/ | |||
962 | /* This may not be quite working correctly yet . . . . */ | |||
963 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
964 | wcycle, enerd, NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0), mu_tot, | |||
965 | constr, NULL((void*)0), FALSE0, state->box, | |||
966 | top_global, &bSumEkinhOld, | |||
967 | CGLO_RERUNMD(1<<1) | CGLO_GSTAT(1<<4) | CGLO_TEMPERATURE(1<<7)); | |||
968 | } | |||
969 | clear_mat(force_vir); | |||
970 | ||||
971 | /* We write a checkpoint at this MD step when: | |||
972 | * either at an NS step when we signalled through gs, | |||
973 | * or at the last step (but not when we do not want confout), | |||
974 | * but never at the first step or with rerun. | |||
975 | */ | |||
976 | bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) || | |||
977 | (bLastStep && (Flags & MD_CONFOUT(1<<12)))) && | |||
978 | step > ir->init_step && !bRerunMD); | |||
979 | if (bCPT) | |||
980 | { | |||
981 | gs.set[eglsCHKPT] = 0; | |||
982 | } | |||
983 | ||||
984 | /* Determine the energy and pressure: | |||
985 | * at nstcalcenergy steps and at energy output steps (set below). | |||
986 | */ | |||
987 | if (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) && (!bInitStep)) | |||
988 | { | |||
989 | /* for vv, the first half of the integration actually corresponds | |||
990 | to the previous step. bCalcEner is only required to be evaluated on the 'next' step, | |||
991 | but the virial needs to be calculated on both the current step and the 'next' step. Future | |||
992 | reorganization may be able to get rid of one of the bCalcVir=TRUE steps. */ | |||
993 | ||||
994 | bCalcEner = do_per_step(step-1, ir->nstcalcenergy); | |||
995 | bCalcVir = bCalcEner || | |||
996 | (ir->epc != epcNO && (do_per_step(step, ir->nstpcouple) || do_per_step(step-1, ir->nstpcouple))); | |||
997 | } | |||
998 | else | |||
999 | { | |||
1000 | bCalcEner = do_per_step(step, ir->nstcalcenergy); | |||
1001 | bCalcVir = bCalcEner || | |||
1002 | (ir->epc != epcNO && do_per_step(step, ir->nstpcouple)); | |||
1003 | } | |||
1004 | ||||
1005 | /* Do we need global communication ? */ | |||
1006 | bGStat = (bCalcVir || bCalcEner || bStopCM || | |||
1007 | do_per_step(step, nstglobalcomm) || (bVV && IR_NVT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER ))) && do_per_step(step-1, nstglobalcomm)) || | |||
1008 | (ir->nstlist == -1 && !bRerunMD && step >= nlh.step_nscheck)); | |||
1009 | ||||
1010 | do_ene = (do_per_step(step, ir->nstenergy) || bLastStep); | |||
1011 | ||||
1012 | if (do_ene || do_log) | |||
1013 | { | |||
1014 | bCalcVir = TRUE1; | |||
1015 | bCalcEner = TRUE1; | |||
1016 | bGStat = TRUE1; | |||
1017 | } | |||
1018 | ||||
1019 | /* these CGLO_ options remain the same throughout the iteration */ | |||
1020 | cglo_flags = ((bRerunMD ? CGLO_RERUNMD(1<<1) : 0) | | |||
1021 | (bGStat ? CGLO_GSTAT(1<<4) : 0) | |||
1022 | ); | |||
1023 | ||||
1024 | force_flags = (GMX_FORCE_STATECHANGED(1<<0) | | |||
1025 | ((DYNAMIC_BOX(*ir)((*ir).epc != epcNO || (*ir).eI == eiTPI || ((*ir).deform[0][ 0] != 0 || (*ir).deform[1][1] != 0 || (*ir).deform[2][2] != 0 || (*ir).deform[1][0] != 0 || (*ir).deform[2][0] != 0 || (*ir ).deform[2][1] != 0)) || bRerunMD) ? GMX_FORCE_DYNAMICBOX(1<<1) : 0) | | |||
1026 | GMX_FORCE_ALLFORCES((1<<4) | (1<<6) | (1<<7)) | | |||
1027 | GMX_FORCE_SEPLRF(1<<5) | | |||
1028 | (bCalcVir ? GMX_FORCE_VIRIAL(1<<8) : 0) | | |||
1029 | (bCalcEner ? GMX_FORCE_ENERGY(1<<9) : 0) | | |||
1030 | (bDoFEP ? GMX_FORCE_DHDL(1<<10) : 0) | |||
1031 | ); | |||
1032 | ||||
1033 | if (fr->bTwinRange) | |||
1034 | { | |||
1035 | if (do_per_step(step, ir->nstcalclr)) | |||
1036 | { | |||
1037 | force_flags |= GMX_FORCE_DO_LR(1<<11); | |||
1038 | } | |||
1039 | } | |||
1040 | ||||
1041 | if (shellfc) | |||
1042 | { | |||
1043 | /* Now is the time to relax the shells */ | |||
1044 | count = relax_shell_flexcon(fplog, cr, bVerbose, step, | |||
1045 | ir, bNS, force_flags, | |||
1046 | top, | |||
1047 | constr, enerd, fcd, | |||
1048 | state, f, force_vir, mdatoms, | |||
1049 | nrnb, wcycle, graph, groups, | |||
1050 | shellfc, fr, bBornRadii, t, mu_tot, | |||
1051 | &bConverged, vsite, | |||
1052 | mdoutf_get_fp_field(outf)); | |||
1053 | tcount += count; | |||
1054 | ||||
1055 | if (bConverged) | |||
1056 | { | |||
1057 | nconverged++; | |||
1058 | } | |||
1059 | } | |||
1060 | else | |||
1061 | { | |||
1062 | /* The coordinates (x) are shifted (to get whole molecules) | |||
1063 | * in do_force. | |||
1064 | * This is parallellized as well, and does communication too. | |||
1065 | * Check comments in sim_util.c | |||
1066 | */ | |||
1067 | do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups, | |||
1068 | state->box, state->x, &state->hist, | |||
1069 | f, force_vir, mdatoms, enerd, fcd, | |||
1070 | state->lambda, graph, | |||
1071 | fr, vsite, mu_tot, t, mdoutf_get_fp_field(outf), ed, bBornRadii, | |||
1072 | (bNS ? GMX_FORCE_NS(1<<2) : 0) | force_flags); | |||
1073 | } | |||
1074 | ||||
1075 | if (bVV && !bStartingFromCpt && !bRerunMD) | |||
1076 | /* ############### START FIRST UPDATE HALF-STEP FOR VV METHODS############### */ | |||
1077 | { | |||
1078 | if (ir->eI == eiVV && bInitStep) | |||
1079 | { | |||
1080 | /* if using velocity verlet with full time step Ekin, | |||
1081 | * take the first half step only to compute the | |||
1082 | * virial for the first step. From there, | |||
1083 | * revert back to the initial coordinates | |||
1084 | * so that the input is actually the initial step. | |||
1085 | */ | |||
1086 | copy_rvecn(state->v, cbuf, 0, state->natoms); /* should make this better for parallelizing? */ | |||
1087 | } | |||
1088 | else | |||
1089 | { | |||
1090 | /* this is for NHC in the Ekin(t+dt/2) version of vv */ | |||
1091 | trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ1); | |||
1092 | } | |||
1093 | ||||
1094 | /* If we are using twin-range interactions where the long-range component | |||
1095 | * is only evaluated every nstcalclr>1 steps, we should do a special update | |||
1096 | * step to combine the long-range forces on these steps. | |||
1097 | * For nstcalclr=1 this is not done, since the forces would have been added | |||
1098 | * directly to the short-range forces already. | |||
1099 | */ | |||
1100 | bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); | |||
1101 | ||||
1102 | update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, | |||
1103 | f, bUpdateDoLR, fr->f_twin, fcd, | |||
1104 | ekind, M, upd, bInitStep, etrtVELOCITY1, | |||
1105 | cr, nrnb, constr, &top->idef); | |||
1106 | ||||
1107 | if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep) | |||
1108 | { | |||
1109 | gmx_iterate_init(&iterate, TRUE1); | |||
1110 | } | |||
1111 | /* for iterations, we save these vectors, as we will be self-consistently iterating | |||
1112 | the calculations */ | |||
1113 | ||||
1114 | /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */ | |||
1115 | ||||
1116 | /* save the state */ | |||
1117 | if (iterate.bIterationActive) | |||
1118 | { | |||
1119 | copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts)); | |||
1120 | } | |||
1121 | ||||
1122 | bFirstIterate = TRUE1; | |||
1123 | while (bFirstIterate || iterate.bIterationActive) | |||
1124 | { | |||
1125 | if (iterate.bIterationActive) | |||
1126 | { | |||
1127 | copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts)); | |||
1128 | if (bFirstIterate && bTrotter) | |||
1129 | { | |||
1130 | /* The first time through, we need a decent first estimate | |||
1131 | of veta(t+dt) to compute the constraints. Do | |||
1132 | this by computing the box volume part of the | |||
1133 | trotter integration at this time. Nothing else | |||
1134 | should be changed by this routine here. If | |||
1135 | !(first time), we start with the previous value | |||
1136 | of veta. */ | |||
1137 | ||||
1138 | veta_save = state->veta; | |||
1139 | trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0); | |||
1140 | vetanew = state->veta; | |||
1141 | state->veta = veta_save; | |||
1142 | } | |||
1143 | } | |||
1144 | ||||
1145 | bOK = TRUE1; | |||
1146 | if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */ | |||
1147 | { | |||
1148 | update_constraints(fplog, step, NULL((void*)0), ir, ekind, mdatoms, | |||
1149 | state, fr->bMolPBC, graph, f, | |||
1150 | &top->idef, shake_vir, | |||
1151 | cr, nrnb, wcycle, upd, constr, | |||
1152 | TRUE1, bCalcVir, vetanew); | |||
1153 | ||||
1154 | if (!bOK) | |||
1155 | { | |||
1156 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 1156, "Constraint error: Shake, Lincs or Settle could not solve the constrains"); | |||
1157 | } | |||
1158 | ||||
1159 | } | |||
1160 | else if (graph) | |||
1161 | { | |||
1162 | /* Need to unshift here if a do_force has been | |||
1163 | called in the previous step */ | |||
1164 | unshift_self(graph, state->box, state->x); | |||
1165 | } | |||
1166 | ||||
1167 | /* if VV, compute the pressure and constraints */ | |||
1168 | /* For VV2, we strictly only need this if using pressure | |||
1169 | * control, but we really would like to have accurate pressures | |||
1170 | * printed out. | |||
1171 | * Think about ways around this in the future? | |||
1172 | * For now, keep this choice in comments. | |||
1173 | */ | |||
1174 | /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */ | |||
1175 | /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/ | |||
1176 | bPres = TRUE1; | |||
1177 | bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK)); | |||
1178 | if (bCalcEner && ir->eI == eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/ | |||
1179 | { | |||
1180 | bSumEkinhOld = TRUE1; | |||
1181 | } | |||
1182 | /* for vv, the first half of the integration actually corresponds to the previous step. | |||
1183 | So we need information from the last step in the first half of the integration */ | |||
1184 | if (bGStat || do_per_step(step-1, nstglobalcomm)) | |||
1185 | { | |||
1186 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
1187 | wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot, | |||
1188 | constr, NULL((void*)0), FALSE0, state->box, | |||
1189 | top_global, &bSumEkinhOld, | |||
1190 | cglo_flags | |||
1191 | | CGLO_ENERGY(1<<6) | |||
1192 | | (bTemp ? CGLO_TEMPERATURE(1<<7) : 0) | |||
1193 | | (bPres ? CGLO_PRESSURE(1<<8) : 0) | |||
1194 | | (bPres ? CGLO_CONSTRAINT(1<<9) : 0) | |||
1195 | | ((iterate.bIterationActive) ? CGLO_ITERATE(1<<10) : 0) | |||
1196 | | (bFirstIterate ? CGLO_FIRSTITERATE(1<<11) : 0) | |||
1197 | | CGLO_SCALEEKIN(1<<13) | |||
1198 | ); | |||
1199 | /* explanation of above: | |||
1200 | a) We compute Ekin at the full time step | |||
1201 | if 1) we are using the AveVel Ekin, and it's not the | |||
1202 | initial step, or 2) if we are using AveEkin, but need the full | |||
1203 | time step kinetic energy for the pressure (always true now, since we want accurate statistics). | |||
1204 | b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in | |||
1205 | EkinAveVel because it's needed for the pressure */ | |||
1206 | } | |||
1207 | /* temperature scaling and pressure scaling to produce the extended variables at t+dt */ | |||
1208 | if (!bInitStep) | |||
1209 | { | |||
1210 | if (bTrotter) | |||
1211 | { | |||
1212 | m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */ | |||
1213 | trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2); | |||
1214 | } | |||
1215 | else | |||
1216 | { | |||
1217 | if (bExchanged) | |||
1218 | { | |||
1219 | ||||
1220 | /* We need the kinetic energy at minus the half step for determining | |||
1221 | * the full step kinetic energy and possibly for T-coupling.*/ | |||
1222 | /* This may not be quite working correctly yet . . . . */ | |||
1223 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
1224 | wcycle, enerd, NULL((void*)0), NULL((void*)0), NULL((void*)0), NULL((void*)0), mu_tot, | |||
1225 | constr, NULL((void*)0), FALSE0, state->box, | |||
1226 | top_global, &bSumEkinhOld, | |||
1227 | CGLO_RERUNMD(1<<1) | CGLO_GSTAT(1<<4) | CGLO_TEMPERATURE(1<<7)); | |||
1228 | } | |||
1229 | } | |||
1230 | } | |||
1231 | ||||
1232 | if (iterate.bIterationActive && | |||
1233 | done_iterating(cr, fplog, step, &iterate, bFirstIterate, | |||
1234 | state->veta, &vetanew)) | |||
1235 | { | |||
1236 | break; | |||
1237 | } | |||
1238 | bFirstIterate = FALSE0; | |||
1239 | } | |||
1240 | ||||
1241 | if (bTrotter && !bInitStep) | |||
1242 | { | |||
1243 | copy_mat(shake_vir, state->svir_prev); | |||
1244 | copy_mat(force_vir, state->fvir_prev); | |||
1245 | if (IR_NVT_TROTTER(ir)((((ir)->eI == eiVV) || ((ir)->eI == eiVVAK)) && ((!((ir)->epc == epcMTTK)) && ((ir)->etc == etcNOSEHOOVER ))) && ir->eI == eiVV) | |||
1246 | { | |||
1247 | /* update temperature and kinetic energy now that step is over - this is the v(t+dt) point */ | |||
1248 | enerd->term[F_TEMP] = sum_ekin(&(ir->opts), ekind, NULL((void*)0), (ir->eI == eiVV), FALSE0); | |||
1249 | enerd->term[F_EKIN] = trace(ekind->ekin); | |||
1250 | } | |||
1251 | } | |||
1252 | /* if it's the initial step, we performed this first step just to get the constraint virial */ | |||
1253 | if (bInitStep && ir->eI == eiVV) | |||
1254 | { | |||
1255 | copy_rvecn(cbuf, state->v, 0, state->natoms); | |||
1256 | } | |||
1257 | } | |||
1258 | ||||
1259 | /* MRS -- now done iterating -- compute the conserved quantity */ | |||
1260 | if (bVV) | |||
1261 | { | |||
1262 | saved_conserved_quantity = compute_conserved_from_auxiliary(ir, state, &MassQ); | |||
1263 | if (ir->eI == eiVV) | |||
1264 | { | |||
1265 | last_ekin = enerd->term[F_EKIN]; | |||
1266 | } | |||
1267 | if ((ir->eDispCorr != edispcEnerPres) && (ir->eDispCorr != edispcAllEnerPres)) | |||
1268 | { | |||
1269 | saved_conserved_quantity -= enerd->term[F_DISPCORR]; | |||
1270 | } | |||
1271 | /* sum up the foreign energy and dhdl terms for vv. currently done every step so that dhdl is correct in the .edr */ | |||
1272 | if (!bRerunMD) | |||
1273 | { | |||
1274 | sum_dhdl(enerd, state->lambda, ir->fepvals); | |||
1275 | } | |||
1276 | } | |||
1277 | ||||
1278 | /* ######## END FIRST UPDATE STEP ############## */ | |||
1279 | /* ######## If doing VV, we now have v(dt) ###### */ | |||
1280 | if (bDoExpanded) | |||
1281 | { | |||
1282 | /* perform extended ensemble sampling in lambda - we don't | |||
1283 | actually move to the new state before outputting | |||
1284 | statistics, but if performing simulated tempering, we | |||
1285 | do update the velocities and the tau_t. */ | |||
1286 | ||||
1287 | lamnew = ExpandedEnsembleDynamics(fplog, ir, enerd, state, &MassQ, state->fep_state, &state->dfhist, step, state->v, mdatoms); | |||
1288 | /* history is maintained in state->dfhist, but state_global is what is sent to trajectory and log output */ | |||
1289 | copy_df_history(&state_global->dfhist, &state->dfhist); | |||
1290 | } | |||
1291 | ||||
1292 | /* Now we have the energies and forces corresponding to the | |||
1293 | * coordinates at time t. We must output all of this before | |||
1294 | * the update. | |||
1295 | */ | |||
1296 | do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, | |||
1297 | ir, state, state_global, top_global, fr, | |||
1298 | outf, mdebin, ekind, f, f_global, | |||
1299 | wcycle, &nchkpt, | |||
1300 | bCPT, bRerunMD, bLastStep, (Flags & MD_CONFOUT(1<<12)), | |||
1301 | bSumEkinhOld); | |||
1302 | /* Check if IMD step and do IMD communication, if bIMD is TRUE. */ | |||
1303 | bIMDstep = do_IMD(ir->bIMD, step, cr, bNS, state->box, state->x, ir, t, wcycle); | |||
1304 | ||||
1305 | /* kludge -- virial is lost with restart for NPT control. Must restart */ | |||
1306 | if (bStartingFromCpt && bVV) | |||
1307 | { | |||
1308 | copy_mat(state->svir_prev, shake_vir); | |||
1309 | copy_mat(state->fvir_prev, force_vir); | |||
1310 | } | |||
1311 | ||||
1312 | elapsed_time = walltime_accounting_get_current_elapsed_time(walltime_accounting); | |||
1313 | ||||
1314 | /* Check whether everything is still allright */ | |||
1315 | if (((int)gmx_get_stop_condition() > handled_stop_condition) | |||
1316 | #ifdef GMX_THREAD_MPI | |||
1317 | && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) | |||
1318 | #endif | |||
1319 | ) | |||
1320 | { | |||
1321 | /* this is just make gs.sig compatible with the hack | |||
1322 | of sending signals around by MPI_Reduce with together with | |||
1323 | other floats */ | |||
1324 | if (gmx_get_stop_condition() == gmx_stop_cond_next_ns) | |||
1325 | { | |||
1326 | gs.sig[eglsSTOPCOND] = 1; | |||
1327 | } | |||
1328 | if (gmx_get_stop_condition() == gmx_stop_cond_next) | |||
1329 | { | |||
1330 | gs.sig[eglsSTOPCOND] = -1; | |||
1331 | } | |||
1332 | /* < 0 means stop at next step, > 0 means stop at next NS step */ | |||
1333 | if (fplog) | |||
1334 | { | |||
1335 | fprintf(fplog, | |||
1336 | "\n\nReceived the %s signal, stopping at the next %sstep\n\n", | |||
1337 | gmx_get_signal_name(), | |||
1338 | gs.sig[eglsSTOPCOND] == 1 ? "NS " : ""); | |||
1339 | fflush(fplog); | |||
1340 | } | |||
1341 | fprintf(stderrstderr, | |||
1342 | "\n\nReceived the %s signal, stopping at the next %sstep\n\n", | |||
1343 | gmx_get_signal_name(), | |||
1344 | gs.sig[eglsSTOPCOND] == 1 ? "NS " : ""); | |||
1345 | fflush(stderrstderr); | |||
1346 | handled_stop_condition = (int)gmx_get_stop_condition(); | |||
1347 | } | |||
1348 | else if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && (bNS || ir->nstlist <= 0) && | |||
1349 | (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) && | |||
1350 | gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0) | |||
1351 | { | |||
1352 | /* Signal to terminate the run */ | |||
1353 | gs.sig[eglsSTOPCOND] = 1; | |||
1354 | if (fplog) | |||
1355 | { | |||
1356 | fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99); | |||
1357 | } | |||
1358 | fprintf(stderrstderr, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99); | |||
1359 | } | |||
1360 | ||||
1361 | if (bResetCountersHalfMaxH && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && | |||
1362 | elapsed_time > max_hours*60.0*60.0*0.495) | |||
1363 | { | |||
1364 | gs.sig[eglsRESETCOUNTERS] = 1; | |||
1365 | } | |||
1366 | ||||
1367 | if (ir->nstlist == -1 && !bRerunMD) | |||
1368 | { | |||
1369 | /* When bGStatEveryStep=FALSE, global_stat is only called | |||
1370 | * when we check the atom displacements, not at NS steps. | |||
1371 | * This means that also the bonded interaction count check is not | |||
1372 | * performed immediately after NS. Therefore a few MD steps could | |||
1373 | * be performed with missing interactions. | |||
1374 | * But wrong energies are never written to file, | |||
1375 | * since energies are only written after global_stat | |||
1376 | * has been called. | |||
1377 | */ | |||
1378 | if (step >= nlh.step_nscheck) | |||
| ||||
1379 | { | |||
1380 | nlh.nabnsb = natoms_beyond_ns_buffer(ir, fr, &top->cgs, | |||
1381 | nlh.scale_tot, state->x); | |||
1382 | } | |||
1383 | else | |||
1384 | { | |||
1385 | /* This is not necessarily true, | |||
1386 | * but step_nscheck is determined quite conservatively. | |||
1387 | */ | |||
1388 | nlh.nabnsb = 0; | |||
1389 | } | |||
1390 | } | |||
1391 | ||||
1392 | /* In parallel we only have to check for checkpointing in steps | |||
1393 | * where we do global communication, | |||
1394 | * otherwise the other nodes don't know. | |||
1395 | */ | |||
1396 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && ((bGStat || !PAR(cr)((cr)->nnodes > 1)) && | |||
1397 | cpt_period >= 0 && | |||
1398 | (cpt_period == 0 || | |||
1399 | elapsed_time >= nchkpt*cpt_period*60.0)) && | |||
1400 | gs.set[eglsCHKPT] == 0) | |||
1401 | { | |||
1402 | gs.sig[eglsCHKPT] = 1; | |||
1403 | } | |||
1404 | ||||
1405 | /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */ | |||
1406 | if (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK)) | |||
1407 | { | |||
1408 | if (!bInitStep) | |||
1409 | { | |||
1410 | update_tcouple(step, ir, state, ekind, &MassQ, mdatoms); | |||
1411 | } | |||
1412 | if (ETC_ANDERSEN(ir->etc)(((ir->etc) == etcANDERSENMASSIVE) || ((ir->etc) == etcANDERSEN ))) /* keep this outside of update_tcouple because of the extra info required to pass */ | |||
1413 | { | |||
1414 | gmx_bool bIfRandomize; | |||
1415 | bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr); | |||
1416 | /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */ | |||
1417 | if (constr && bIfRandomize) | |||
1418 | { | |||
1419 | update_constraints(fplog, step, NULL((void*)0), ir, ekind, mdatoms, | |||
1420 | state, fr->bMolPBC, graph, f, | |||
1421 | &top->idef, tmp_vir, | |||
1422 | cr, nrnb, wcycle, upd, constr, | |||
1423 | TRUE1, bCalcVir, vetanew); | |||
1424 | } | |||
1425 | } | |||
1426 | } | |||
1427 | ||||
1428 | if (bIterativeCase && do_per_step(step, ir->nstpcouple)) | |||
1429 | { | |||
1430 | gmx_iterate_init(&iterate, TRUE1); | |||
1431 | /* for iterations, we save these vectors, as we will be redoing the calculations */ | |||
1432 | copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts)); | |||
1433 | } | |||
1434 | ||||
1435 | bFirstIterate = TRUE1; | |||
1436 | while (bFirstIterate || iterate.bIterationActive) | |||
1437 | { | |||
1438 | /* We now restore these vectors to redo the calculation with improved extended variables */ | |||
1439 | if (iterate.bIterationActive) | |||
1440 | { | |||
1441 | copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts)); | |||
1442 | } | |||
1443 | ||||
1444 | /* We make the decision to break or not -after- the calculation of Ekin and Pressure, | |||
1445 | so scroll down for that logic */ | |||
1446 | ||||
1447 | /* ######### START SECOND UPDATE STEP ################# */ | |||
1448 | /* Box is changed in update() when we do pressure coupling, | |||
1449 | * but we should still use the old box for energy corrections and when | |||
1450 | * writing it to the energy file, so it matches the trajectory files for | |||
1451 | * the same timestep above. Make a copy in a separate array. | |||
1452 | */ | |||
1453 | copy_mat(state->box, lastbox); | |||
1454 | ||||
1455 | bOK = TRUE1; | |||
1456 | dvdl_constr = 0; | |||
1457 | ||||
1458 | if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate)) | |||
1459 | { | |||
1460 | wallcycle_start(wcycle, ewcUPDATE); | |||
1461 | /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */ | |||
1462 | if (bTrotter) | |||
1463 | { | |||
1464 | if (iterate.bIterationActive) | |||
1465 | { | |||
1466 | if (bFirstIterate) | |||
1467 | { | |||
1468 | scalevir = 1; | |||
1469 | } | |||
1470 | else | |||
1471 | { | |||
1472 | /* we use a new value of scalevir to converge the iterations faster */ | |||
1473 | scalevir = tracevir/trace(shake_vir); | |||
1474 | } | |||
1475 | msmul(shake_vir, scalevir, shake_vir); | |||
1476 | m_add(force_vir, shake_vir, total_vir); | |||
1477 | clear_mat(shake_vir); | |||
1478 | } | |||
1479 | trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3); | |||
1480 | /* We can only do Berendsen coupling after we have summed | |||
1481 | * the kinetic energy or virial. Since the happens | |||
1482 | * in global_state after update, we should only do it at | |||
1483 | * step % nstlist = 1 with bGStatEveryStep=FALSE. | |||
1484 | */ | |||
1485 | } | |||
1486 | else | |||
1487 | { | |||
1488 | update_tcouple(step, ir, state, ekind, &MassQ, mdatoms); | |||
1489 | update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep); | |||
1490 | } | |||
1491 | ||||
1492 | if (bVV) | |||
1493 | { | |||
1494 | bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); | |||
1495 | ||||
1496 | /* velocity half-step update */ | |||
1497 | update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, | |||
1498 | bUpdateDoLR, fr->f_twin, fcd, | |||
1499 | ekind, M, upd, FALSE0, etrtVELOCITY2, | |||
1500 | cr, nrnb, constr, &top->idef); | |||
1501 | } | |||
1502 | ||||
1503 | /* Above, initialize just copies ekinh into ekin, | |||
1504 | * it doesn't copy position (for VV), | |||
1505 | * and entire integrator for MD. | |||
1506 | */ | |||
1507 | ||||
1508 | if (ir->eI == eiVVAK) | |||
1509 | { | |||
1510 | copy_rvecn(state->x, cbuf, 0, state->natoms); | |||
1511 | } | |||
1512 | bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); | |||
1513 | ||||
1514 | update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, | |||
1515 | bUpdateDoLR, fr->f_twin, fcd, | |||
1516 | ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef); | |||
1517 | wallcycle_stop(wcycle, ewcUPDATE); | |||
1518 | ||||
1519 | update_constraints(fplog, step, &dvdl_constr, ir, ekind, mdatoms, state, | |||
1520 | fr->bMolPBC, graph, f, | |||
1521 | &top->idef, shake_vir, | |||
1522 | cr, nrnb, wcycle, upd, constr, | |||
1523 | FALSE0, bCalcVir, state->veta); | |||
1524 | ||||
1525 | if (ir->eI == eiVVAK) | |||
1526 | { | |||
1527 | /* erase F_EKIN and F_TEMP here? */ | |||
1528 | /* just compute the kinetic energy at the half step to perform a trotter step */ | |||
1529 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
1530 | wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot, | |||
1531 | constr, NULL((void*)0), FALSE0, lastbox, | |||
1532 | top_global, &bSumEkinhOld, | |||
1533 | cglo_flags | CGLO_TEMPERATURE(1<<7) | |||
1534 | ); | |||
1535 | wallcycle_start(wcycle, ewcUPDATE); | |||
1536 | trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4); | |||
1537 | /* now we know the scaling, we can compute the positions again again */ | |||
1538 | copy_rvecn(cbuf, state->x, 0, state->natoms); | |||
1539 | ||||
1540 | bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); | |||
1541 | ||||
1542 | update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, | |||
1543 | bUpdateDoLR, fr->f_twin, fcd, | |||
1544 | ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef); | |||
1545 | wallcycle_stop(wcycle, ewcUPDATE); | |||
1546 | ||||
1547 | /* do we need an extra constraint here? just need to copy out of state->v to upd->xp? */ | |||
1548 | /* are the small terms in the shake_vir here due | |||
1549 | * to numerical errors, or are they important | |||
1550 | * physically? I'm thinking they are just errors, but not completely sure. | |||
1551 | * For now, will call without actually constraining, constr=NULL*/ | |||
1552 | update_constraints(fplog, step, NULL((void*)0), ir, ekind, mdatoms, | |||
1553 | state, fr->bMolPBC, graph, f, | |||
1554 | &top->idef, tmp_vir, | |||
1555 | cr, nrnb, wcycle, upd, NULL((void*)0), | |||
1556 | FALSE0, bCalcVir, | |||
1557 | state->veta); | |||
1558 | } | |||
1559 | if (!bOK) | |||
1560 | { | |||
1561 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 1561, "Constraint error: Shake, Lincs or Settle could not solve the constrains"); | |||
1562 | } | |||
1563 | ||||
1564 | if (fr->bSepDVDL && fplog && do_log) | |||
1565 | { | |||
1566 | gmx_print_sepdvdl(fplog, "Constraint dV/dl", 0.0, dvdl_constr); | |||
1567 | } | |||
1568 | if (bVV) | |||
1569 | { | |||
1570 | /* this factor or 2 correction is necessary | |||
1571 | because half of the constraint force is removed | |||
1572 | in the vv step, so we have to double it. See | |||
1573 | the Redmine issue #1255. It is not yet clear | |||
1574 | if the factor of 2 is exact, or just a very | |||
1575 | good approximation, and this will be | |||
1576 | investigated. The next step is to see if this | |||
1577 | can be done adding a dhdl contribution from the | |||
1578 | rattle step, but this is somewhat more | |||
1579 | complicated with the current code. Will be | |||
1580 | investigated, hopefully for 4.6.3. However, | |||
1581 | this current solution is much better than | |||
1582 | having it completely wrong. | |||
1583 | */ | |||
1584 | enerd->term[F_DVDL_CONSTR] += 2*dvdl_constr; | |||
1585 | } | |||
1586 | else | |||
1587 | { | |||
1588 | enerd->term[F_DVDL_CONSTR] += dvdl_constr; | |||
1589 | } | |||
1590 | } | |||
1591 | else if (graph) | |||
1592 | { | |||
1593 | /* Need to unshift here */ | |||
1594 | unshift_self(graph, state->box, state->x); | |||
1595 | } | |||
1596 | ||||
1597 | if (vsite != NULL((void*)0)) | |||
1598 | { | |||
1599 | wallcycle_start(wcycle, ewcVSITECONSTR); | |||
1600 | if (graph != NULL((void*)0)) | |||
1601 | { | |||
1602 | shift_self(graph, state->box, state->x); | |||
1603 | } | |||
1604 | construct_vsites(vsite, state->x, ir->delta_t, state->v, | |||
1605 | top->idef.iparams, top->idef.il, | |||
1606 | fr->ePBC, fr->bMolPBC, cr, state->box); | |||
1607 | ||||
1608 | if (graph != NULL((void*)0)) | |||
1609 | { | |||
1610 | unshift_self(graph, state->box, state->x); | |||
1611 | } | |||
1612 | wallcycle_stop(wcycle, ewcVSITECONSTR); | |||
1613 | } | |||
1614 | ||||
1615 | /* ############## IF NOT VV, Calculate globals HERE, also iterate constraints ############ */ | |||
1616 | /* With Leap-Frog we can skip compute_globals at | |||
1617 | * non-communication steps, but we need to calculate | |||
1618 | * the kinetic energy one step before communication. | |||
1619 | */ | |||
1620 | if (bGStat || (!EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) && do_per_step(step+1, nstglobalcomm))) | |||
1621 | { | |||
1622 | if (ir->nstlist == -1 && bFirstIterate) | |||
1623 | { | |||
1624 | gs.sig[eglsNABNSB] = nlh.nabnsb; | |||
1625 | } | |||
1626 | compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, | |||
1627 | wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot, | |||
1628 | constr, | |||
1629 | bFirstIterate ? &gs : NULL((void*)0), | |||
1630 | (step_rel % gs.nstms == 0) && | |||
1631 | (multisim_nsteps < 0 || (step_rel < multisim_nsteps)), | |||
1632 | lastbox, | |||
1633 | top_global, &bSumEkinhOld, | |||
1634 | cglo_flags | |||
1635 | | (!EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) || bRerunMD ? CGLO_ENERGY(1<<6) : 0) | |||
1636 | | (!EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) && bStopCM ? CGLO_STOPCM(1<<3) : 0) | |||
1637 | | (!EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) ? CGLO_TEMPERATURE(1<<7) : 0) | |||
1638 | | (!EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK) || bRerunMD ? CGLO_PRESSURE(1<<8) : 0) | |||
1639 | | (iterate.bIterationActive ? CGLO_ITERATE(1<<10) : 0) | |||
1640 | | (bFirstIterate ? CGLO_FIRSTITERATE(1<<11) : 0) | |||
1641 | | CGLO_CONSTRAINT(1<<9) | |||
1642 | ); | |||
1643 | if (ir->nstlist == -1 && bFirstIterate) | |||
1644 | { | |||
1645 | nlh.nabnsb = gs.set[eglsNABNSB]; | |||
1646 | gs.set[eglsNABNSB] = 0; | |||
1647 | } | |||
1648 | } | |||
1649 | /* bIterate is set to keep it from eliminating the old ekin kinetic energy terms */ | |||
1650 | /* ############# END CALC EKIN AND PRESSURE ################# */ | |||
1651 | ||||
1652 | /* Note: this is OK, but there are some numerical precision issues with using the convergence of | |||
1653 | the virial that should probably be addressed eventually. state->veta has better properies, | |||
1654 | but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could | |||
1655 | generate the new shake_vir, but test the veta value for convergence. This will take some thought. */ | |||
1656 | ||||
1657 | if (iterate.bIterationActive && | |||
1658 | done_iterating(cr, fplog, step, &iterate, bFirstIterate, | |||
1659 | trace(shake_vir), &tracevir)) | |||
1660 | { | |||
1661 | break; | |||
1662 | } | |||
1663 | bFirstIterate = FALSE0; | |||
1664 | } | |||
1665 | ||||
1666 | if (!bVV || bRerunMD) | |||
1667 | { | |||
1668 | /* sum up the foreign energy and dhdl terms for md and sd. currently done every step so that dhdl is correct in the .edr */ | |||
1669 | sum_dhdl(enerd, state->lambda, ir->fepvals); | |||
1670 | } | |||
1671 | update_box(fplog, step, ir, mdatoms, state, f, | |||
1672 | ir->nstlist == -1 ? &nlh.scale_tot : NULL((void*)0), pcoupl_mu, nrnb, upd); | |||
1673 | ||||
1674 | /* ################# END UPDATE STEP 2 ################# */ | |||
1675 | /* #### We now have r(t+dt) and v(t+dt/2) ############# */ | |||
1676 | ||||
1677 | /* The coordinates (x) were unshifted in update */ | |||
1678 | if (!bGStat) | |||
1679 | { | |||
1680 | /* We will not sum ekinh_old, | |||
1681 | * so signal that we still have to do it. | |||
1682 | */ | |||
1683 | bSumEkinhOld = TRUE1; | |||
1684 | } | |||
1685 | ||||
1686 | /* ######### BEGIN PREPARING EDR OUTPUT ########### */ | |||
1687 | ||||
1688 | /* use the directly determined last velocity, not actually the averaged half steps */ | |||
1689 | if (bTrotter && ir->eI == eiVV) | |||
1690 | { | |||
1691 | enerd->term[F_EKIN] = last_ekin; | |||
1692 | } | |||
1693 | enerd->term[F_ETOT] = enerd->term[F_EPOT] + enerd->term[F_EKIN]; | |||
1694 | ||||
1695 | if (bVV) | |||
1696 | { | |||
1697 | enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + saved_conserved_quantity; | |||
1698 | } | |||
1699 | else | |||
1700 | { | |||
1701 | enerd->term[F_ECONSERVED] = enerd->term[F_ETOT] + compute_conserved_from_auxiliary(ir, state, &MassQ); | |||
1702 | } | |||
1703 | /* ######### END PREPARING EDR OUTPUT ########### */ | |||
1704 | ||||
1705 | /* Output stuff */ | |||
1706 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
1707 | { | |||
1708 | gmx_bool do_dr, do_or; | |||
1709 | ||||
1710 | if (fplog && do_log && bDoExpanded) | |||
1711 | { | |||
1712 | /* only needed if doing expanded ensemble */ | |||
1713 | PrintFreeEnergyInfoToFile(fplog, ir->fepvals, ir->expandedvals, ir->bSimTemp ? ir->simtempvals : NULL((void*)0), | |||
1714 | &state_global->dfhist, state->fep_state, ir->nstlog, step); | |||
1715 | } | |||
1716 | if (!(bStartingFromCpt && (EI_VV(ir->eI)((ir->eI) == eiVV || (ir->eI) == eiVVAK)))) | |||
1717 | { | |||
1718 | if (bCalcEner) | |||
1719 | { | |||
1720 | upd_mdebin(mdebin, bDoDHDL, TRUE1, | |||
1721 | t, mdatoms->tmass, enerd, state, | |||
1722 | ir->fepvals, ir->expandedvals, lastbox, | |||
1723 | shake_vir, force_vir, total_vir, pres, | |||
1724 | ekind, mu_tot, constr); | |||
1725 | } | |||
1726 | else | |||
1727 | { | |||
1728 | upd_mdebin_step(mdebin); | |||
1729 | } | |||
1730 | ||||
1731 | do_dr = do_per_step(step, ir->nstdisreout); | |||
1732 | do_or = do_per_step(step, ir->nstorireout); | |||
1733 | ||||
1734 | print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : NULL((void*)0), | |||
1735 | step, t, | |||
1736 | eprNORMAL, bCompact, mdebin, fcd, groups, &(ir->opts)); | |||
1737 | } | |||
1738 | if (ir->ePull != epullNO) | |||
1739 | { | |||
1740 | pull_print_output(ir->pull, step, t); | |||
1741 | } | |||
1742 | ||||
1743 | if (do_per_step(step, ir->nstlog)) | |||
1744 | { | |||
1745 | if (fflush(fplog) != 0) | |||
1746 | { | |||
1747 | gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/programs/mdrun/md.c", 1747, "Cannot flush logfile - maybe you are out of disk space?"); | |||
1748 | } | |||
1749 | } | |||
1750 | } | |||
1751 | if (bDoExpanded) | |||
1752 | { | |||
1753 | /* Have to do this part _after_ outputting the logfile and the edr file */ | |||
1754 | /* Gets written into the state at the beginning of next loop*/ | |||
1755 | state->fep_state = lamnew; | |||
1756 | } | |||
1757 | /* Print the remaining wall clock time for the run */ | |||
1758 | if (MULTIMASTER(cr)((((((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && ((cr)->duty & (1<<0))) || !((cr)->nnodes > 1)) && (!((cr)->ms) || (((cr)->ms)->sim == 0 ))) && (do_verbose || gmx_got_usr_signal()) && !bPMETuneRunning) | |||
1759 | { | |||
1760 | if (shellfc) | |||
1761 | { | |||
1762 | fprintf(stderrstderr, "\n"); | |||
1763 | } | |||
1764 | print_time(stderrstderr, walltime_accounting, step, ir, cr); | |||
1765 | } | |||
1766 | ||||
1767 | /* Ion/water position swapping. | |||
1768 | * Not done in last step since trajectory writing happens before this call | |||
1769 | * in the MD loop and exchanges would be lost anyway. */ | |||
1770 | bNeedRepartition = FALSE0; | |||
1771 | if ((ir->eSwapCoords != eswapNO) && (step > 0) && !bLastStep && | |||
1772 | do_per_step(step, ir->swap->nstswap)) | |||
1773 | { | |||
1774 | bNeedRepartition = do_swapcoords(cr, step, t, ir, wcycle, | |||
1775 | bRerunMD ? rerun_fr.x : state->x, | |||
1776 | bRerunMD ? rerun_fr.box : state->box, | |||
1777 | top_global, MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)) && bVerbose, bRerunMD); | |||
1778 | ||||
1779 | if (bNeedRepartition && DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1))) | |||
1780 | { | |||
1781 | dd_collect_state(cr->dd, state, state_global); | |||
1782 | } | |||
1783 | } | |||
1784 | ||||
1785 | /* Replica exchange */ | |||
1786 | bExchanged = FALSE0; | |||
1787 | if ((repl_ex_nst > 0) && (step > 0) && !bLastStep && | |||
1788 | do_per_step(step, repl_ex_nst)) | |||
1789 | { | |||
1790 | bExchanged = replica_exchange(fplog, cr, repl_ex, | |||
1791 | state_global, enerd, | |||
1792 | state, step, t); | |||
1793 | } | |||
1794 | ||||
1795 | if ( (bExchanged || bNeedRepartition) && DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) ) | |||
1796 | { | |||
1797 | dd_partition_system(fplog, step, cr, TRUE1, 1, | |||
1798 | state_global, top_global, ir, | |||
1799 | state, &f, mdatoms, top, fr, | |||
1800 | vsite, shellfc, constr, | |||
1801 | nrnb, wcycle, FALSE0); | |||
1802 | } | |||
1803 | ||||
1804 | bFirstStep = FALSE0; | |||
1805 | bInitStep = FALSE0; | |||
1806 | bStartingFromCpt = FALSE0; | |||
1807 | ||||
1808 | /* ####### SET VARIABLES FOR NEXT ITERATION IF THEY STILL NEED IT ###### */ | |||
1809 | /* With all integrators, except VV, we need to retain the pressure | |||
1810 | * at the current step for coupling at the next step. | |||
1811 | */ | |||
1812 | if ((state->flags & (1<<estPRES_PREV)) && | |||
1813 | (bGStatEveryStep || | |||
1814 | (ir->nstpcouple > 0 && step % ir->nstpcouple == 0))) | |||
1815 | { | |||
1816 | /* Store the pressure in t_state for pressure coupling | |||
1817 | * at the next MD step. | |||
1818 | */ | |||
1819 | copy_mat(pres, state->pres_prev); | |||
1820 | } | |||
1821 | ||||
1822 | /* ####### END SET VARIABLES FOR NEXT ITERATION ###### */ | |||
1823 | ||||
1824 | if ( (membed != NULL((void*)0)) && (!bLastStep) ) | |||
1825 | { | |||
1826 | rescale_membed(step_rel, membed, state_global->x); | |||
1827 | } | |||
1828 | ||||
1829 | if (bRerunMD) | |||
1830 | { | |||
1831 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
1832 | { | |||
1833 | /* read next frame from input trajectory */ | |||
1834 | bNotLastFrame = read_next_frame(oenv, status, &rerun_fr); | |||
1835 | } | |||
1836 | ||||
1837 | if (PAR(cr)((cr)->nnodes > 1)) | |||
1838 | { | |||
1839 | rerun_parallel_comm(cr, &rerun_fr, &bNotLastFrame); | |||
1840 | } | |||
1841 | } | |||
1842 | ||||
1843 | if (!bRerunMD || !rerun_fr.bStep) | |||
1844 | { | |||
1845 | /* increase the MD step number */ | |||
1846 | step++; | |||
1847 | step_rel++; | |||
1848 | } | |||
1849 | ||||
1850 | cycles = wallcycle_stop(wcycle, ewcSTEP); | |||
1851 | if (DOMAINDECOMP(cr)(((cr)->dd != ((void*)0)) && ((cr)->nnodes > 1)) && wcycle) | |||
1852 | { | |||
1853 | dd_cycles_add(cr->dd, cycles, ddCyclStep); | |||
1854 | } | |||
1855 | ||||
1856 | if (bPMETuneRunning || bPMETuneTry) | |||
1857 | { | |||
1858 | /* PME grid + cut-off optimization with GPUs or PME nodes */ | |||
1859 | ||||
1860 | /* Count the total cycles over the last steps */ | |||
1861 | cycles_pmes += cycles; | |||
1862 | ||||
1863 | /* We can only switch cut-off at NS steps */ | |||
1864 | if (step % ir->nstlist == 0) | |||
1865 | { | |||
1866 | /* PME grid + cut-off optimization with GPUs or PME nodes */ | |||
1867 | if (bPMETuneTry) | |||
1868 | { | |||
1869 | if (DDMASTER(cr->dd)((cr->dd)->rank == (cr->dd)->masterrank)) | |||
1870 | { | |||
1871 | /* PME node load is too high, start tuning */ | |||
1872 | bPMETuneRunning = (dd_pme_f_ratio(cr->dd) >= 1.05); | |||
1873 | } | |||
1874 | dd_bcast(cr->dd, sizeof(gmx_bool), &bPMETuneRunning); | |||
1875 | ||||
1876 | if (bPMETuneRunning || step_rel > ir->nstlist*50) | |||
1877 | { | |||
1878 | bPMETuneTry = FALSE0; | |||
1879 | } | |||
1880 | } | |||
1881 | if (bPMETuneRunning) | |||
1882 | { | |||
1883 | /* init_step might not be a multiple of nstlist, | |||
1884 | * but the first cycle is always skipped anyhow. | |||
1885 | */ | |||
1886 | bPMETuneRunning = | |||
1887 | pme_load_balance(pme_loadbal, cr, | |||
1888 | (bVerbose && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) ? stderrstderr : NULL((void*)0), | |||
1889 | fplog, | |||
1890 | ir, state, cycles_pmes, | |||
1891 | fr->ic, fr->nbv, &fr->pmedata, | |||
1892 | step); | |||
1893 | ||||
1894 | /* Update constants in forcerec/inputrec to keep them in sync with fr->ic */ | |||
1895 | fr->ewaldcoeff_q = fr->ic->ewaldcoeff_q; | |||
1896 | fr->ewaldcoeff_lj = fr->ic->ewaldcoeff_lj; | |||
1897 | fr->rlist = fr->ic->rlist; | |||
1898 | fr->rlistlong = fr->ic->rlistlong; | |||
1899 | fr->rcoulomb = fr->ic->rcoulomb; | |||
1900 | fr->rvdw = fr->ic->rvdw; | |||
1901 | ||||
1902 | if (ir->eDispCorr != edispcNO) | |||
1903 | { | |||
1904 | calc_enervirdiff(NULL((void*)0), ir->eDispCorr, fr); | |||
1905 | } | |||
1906 | } | |||
1907 | cycles_pmes = 0; | |||
1908 | } | |||
1909 | } | |||
1910 | ||||
1911 | if (step_rel == wcycle_get_reset_counters(wcycle) || | |||
1912 | gs.set[eglsRESETCOUNTERS] != 0) | |||
1913 | { | |||
1914 | /* Reset all the counters related to performance over the run */ | |||
1915 | reset_all_counters(fplog, cr, step, &step_rel, ir, wcycle, nrnb, walltime_accounting, | |||
1916 | fr->nbv != NULL((void*)0) && fr->nbv->bUseGPU ? fr->nbv->cu_nbv : NULL((void*)0)); | |||
1917 | wcycle_set_reset_counters(wcycle, -1); | |||
1918 | if (!(cr->duty & DUTY_PME(1<<1))) | |||
1919 | { | |||
1920 | /* Tell our PME node to reset its counters */ | |||
1921 | gmx_pme_send_resetcounters(cr, step); | |||
1922 | } | |||
1923 | /* Correct max_hours for the elapsed time */ | |||
1924 | max_hours -= elapsed_time/(60.0*60.0); | |||
1925 | bResetCountersHalfMaxH = FALSE0; | |||
1926 | gs.set[eglsRESETCOUNTERS] = 0; | |||
1927 | } | |||
1928 | ||||
1929 | /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */ | |||
1930 | IMD_prep_energies_send_positions(ir->bIMD && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1)), bIMDstep, ir->imd, enerd, step, bCalcEner, wcycle); | |||
1931 | ||||
1932 | } | |||
1933 | /* End of main MD loop */ | |||
1934 | debug_gmx(); | |||
1935 | ||||
1936 | /* Stop measuring walltime */ | |||
1937 | walltime_accounting_end(walltime_accounting); | |||
1938 | ||||
1939 | if (bRerunMD && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
1940 | { | |||
1941 | close_trj(status); | |||
1942 | } | |||
1943 | ||||
1944 | if (!(cr->duty & DUTY_PME(1<<1))) | |||
1945 | { | |||
1946 | /* Tell the PME only node to finish */ | |||
1947 | gmx_pme_send_finish(cr); | |||
1948 | } | |||
1949 | ||||
1950 | if (MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
1951 | { | |||
1952 | if (ir->nstcalcenergy > 0 && !bRerunMD) | |||
1953 | { | |||
1954 | print_ebin(mdoutf_get_fp_ene(outf), FALSE0, FALSE0, FALSE0, fplog, step, t, | |||
1955 | eprAVER, FALSE0, mdebin, fcd, groups, &(ir->opts)); | |||
1956 | } | |||
1957 | } | |||
1958 | ||||
1959 | done_mdoutf(outf); | |||
1960 | debug_gmx(); | |||
1961 | ||||
1962 | if (ir->nstlist == -1 && nlh.nns > 0 && fplog) | |||
1963 | { | |||
1964 | fprintf(fplog, "Average neighborlist lifetime: %.1f steps, std.dev.: %.1f steps\n", nlh.s1/nlh.nns, sqrt(nlh.s2/nlh.nns - sqr(nlh.s1/nlh.nns))); | |||
1965 | fprintf(fplog, "Average number of atoms that crossed the half buffer length: %.1f\n\n", nlh.ab/nlh.nns); | |||
1966 | } | |||
1967 | ||||
1968 | if (pme_loadbal != NULL((void*)0)) | |||
1969 | { | |||
1970 | pme_loadbal_done(pme_loadbal, cr, fplog, | |||
1971 | fr->nbv != NULL((void*)0) && fr->nbv->bUseGPU); | |||
1972 | } | |||
1973 | ||||
1974 | if (shellfc && fplog) | |||
1975 | { | |||
1976 | fprintf(fplog, "Fraction of iterations that converged: %.2f %%\n", | |||
1977 | (nconverged*100.0)/step_rel); | |||
1978 | fprintf(fplog, "Average number of force evaluations per MD step: %.2f\n\n", | |||
1979 | tcount/step_rel); | |||
1980 | } | |||
1981 | ||||
1982 | if (repl_ex_nst > 0 && MASTER(cr)(((cr)->nodeid == 0) || !((cr)->nnodes > 1))) | |||
1983 | { | |||
1984 | print_replica_exchange_statistics(fplog, repl_ex); | |||
1985 | } | |||
1986 | ||||
1987 | /* IMD cleanup, if bIMD is TRUE. */ | |||
1988 | IMD_finalize(ir->bIMD, ir->imd); | |||
1989 | ||||
1990 | walltime_accounting_set_nsteps_done(walltime_accounting, step_rel); | |||
1991 | ||||
1992 | return 0; | |||
1993 | } |