File: | programs/mdrun/md.c |
Location: | line 234, column 5 |
Description: | Value stored to 'bAppend' is never read |
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)); |
Value stored to 'bAppend' is never read | |
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 | } |