From: Mark Abraham Date: Wed, 14 Oct 2015 15:10:46 +0000 (+0200) Subject: Merge release-5-0 into release-5-1 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=b7fb47a196bd1d9236fe823efde8a968d93602f4;p=alexxy%2Fgromacs.git Merge release-5-0 into release-5-1 Change-Id: I6cdeb2f67b4a0f6c84db9ca086b5600c47320cf0 --- b7fb47a196bd1d9236fe823efde8a968d93602f4 diff --cc src/programs/mdrun/md.cpp index 6c0dc10aa4,0b12964ea5..c9e0a847d7 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@@ -509,37 -536,16 +509,38 @@@ double do_md(FILE *fplog, t_commrec *cr debug_gmx(); - /* set free energy calculation frequency as the minimum - greatest common denominator of nstdhdl, nstexpanded, and repl_ex_nst*/ - nstfep = ir->fepvals->nstdhdl; - if (ir->bExpanded) + if (IR_TWINRANGE(*ir) && repl_ex_nst % ir->nstcalclr != 0) { - nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep); + /* We should exchange at nstcalclr steps to get correct integration */ + gmx_fatal(FARGS, "The replica exchange period (%d) is not divisible by nstcalclr (%d)", repl_ex_nst, ir->nstcalclr); } - if (repl_ex_nst > 0) + + if (ir->efep != efepNO) { - nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep); + /* Set free energy calculation frequency as the greatest common + * denominator of nstdhdl and repl_ex_nst. + * Check for nstcalclr with twin-range, since we need the long-range + * contribution to the free-energy at the correct (nstcalclr) steps. + */ + nstfep = ir->fepvals->nstdhdl; + if (ir->bExpanded) + { + if (IR_TWINRANGE(*ir) && + ir->expandedvals->nstexpanded % ir->nstcalclr != 0) + { + gmx_fatal(FARGS, "nstexpanded should be divisible by nstcalclr"); + } ++ nstfep = gmx_greatest_common_divisor(ir->expandedvals->nstexpanded, nstfep); + } + if (repl_ex_nst > 0) + { + nstfep = gmx_greatest_common_divisor(repl_ex_nst, nstfep); + } + /* We checked divisibility of repl_ex_nst and nstcalclr above */ + if (IR_TWINRANGE(*ir) && nstfep % ir->nstcalclr != 0) + { + gmx_incons("nstfep not divisible by nstcalclr"); + } } /* Be REALLY careful about what flags you set here. You CANNOT assume @@@ -1117,92 -1133,152 +1118,93 @@@ ekind, M, upd, bInitStep, etrtVELOCITY1, cr, nrnb, constr, &top->idef); - if (bIterativeCase && do_per_step(step-1, ir->nstpcouple) && !bInitStep) + if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */ { - gmx_iterate_init(&iterate, TRUE); + wallcycle_stop(wcycle, ewcUPDATE); + update_constraints(fplog, step, NULL, ir, mdatoms, + state, fr->bMolPBC, graph, f, + &top->idef, shake_vir, + cr, nrnb, wcycle, upd, constr, + TRUE, bCalcVir); + wallcycle_start(wcycle, ewcUPDATE); + if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1) + { + /* Correct the virial for multiple time stepping */ + m_sub(shake_vir, fr->vir_twin_constr, shake_vir); + } } - /* for iterations, we save these vectors, as we will be self-consistently iterating - the calculations */ - - /*#### UPDATE EXTENDED VARIABLES IN TROTTER FORMULATION */ - - /* save the state */ - if (iterate.bIterationActive) + else if (graph) { - copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts)); + /* Need to unshift here if a do_force has been + called in the previous step */ + unshift_self(graph, state->box, state->x); } - - bFirstIterate = TRUE; - while (bFirstIterate || iterate.bIterationActive) + /* if VV, compute the pressure and constraints */ + /* For VV2, we strictly only need this if using pressure + * control, but we really would like to have accurate pressures + * printed out. + * Think about ways around this in the future? + * For now, keep this choice in comments. + */ + /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */ + /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/ + bPres = TRUE; + bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK)); + if (bCalcEner && ir->eI == eiVVAK) { - if (iterate.bIterationActive) - { - copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts)); - if (bFirstIterate && bTrotter) - { - /* The first time through, we need a decent first estimate - of veta(t+dt) to compute the constraints. Do - this by computing the box volume part of the - trotter integration at this time. Nothing else - should be changed by this routine here. If - !(first time), we start with the previous value - of veta. */ - - veta_save = state->veta; - trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ0); - vetanew = state->veta; - state->veta = veta_save; - } - } - - bOK = TRUE; - if (!bRerunMD || rerun_fr.bV || bForceUpdate) /* Why is rerun_fr.bV here? Unclear. */ - { - wallcycle_stop(wcycle, ewcUPDATE); - update_constraints(fplog, step, NULL, ir, ekind, mdatoms, - state, fr->bMolPBC, graph, f, - &top->idef, shake_vir, - cr, nrnb, wcycle, upd, constr, - TRUE, bCalcVir, vetanew); - wallcycle_start(wcycle, ewcUPDATE); - - if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1) - { - /* Correct the virial for multiple time stepping */ - m_sub(shake_vir, fr->vir_twin_constr, shake_vir); - } - - if (!bOK) - { - gmx_fatal(FARGS, "Constraint error: Shake, Lincs or Settle could not solve the constrains"); - } - - } - else if (graph) - { - /* Need to unshift here if a do_force has been - called in the previous step */ - unshift_self(graph, state->box, state->x); - } - - /* if VV, compute the pressure and constraints */ - /* For VV2, we strictly only need this if using pressure - * control, but we really would like to have accurate pressures - * printed out. - * Think about ways around this in the future? - * For now, keep this choice in comments. - */ - /*bPres = (ir->eI==eiVV || IR_NPT_TROTTER(ir)); */ - /*bTemp = ((ir->eI==eiVV &&(!bInitStep)) || (ir->eI==eiVVAK && IR_NPT_TROTTER(ir)));*/ - bPres = TRUE; - bTemp = ((ir->eI == eiVV && (!bInitStep)) || (ir->eI == eiVVAK)); - if (bCalcEner && ir->eI == eiVVAK) /*MRS: 7/9/2010 -- this still doesn't fix it?*/ - { - bSumEkinhOld = TRUE; - } - /* for vv, the first half of the integration actually corresponds to the previous step. - So we need information from the last step in the first half of the integration */ - if (bGStat || do_per_step(step-1, nstglobalcomm)) + bSumEkinhOld = TRUE; + } + /* for vv, the first half of the integration actually corresponds to the previous step. + So we need information from the last step in the first half of the integration */ + if (bGStat || do_per_step(step-1, nstglobalcomm)) + { + wallcycle_stop(wcycle, ewcUPDATE); + compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, + wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot, + constr, NULL, FALSE, state->box, + top_global, &bSumEkinhOld, + cglo_flags + | CGLO_ENERGY + | (bTemp ? CGLO_TEMPERATURE : 0) + | (bPres ? CGLO_PRESSURE : 0) + | (bPres ? CGLO_CONSTRAINT : 0) ++ | (bStopCM ? CGLO_STOPCM : 0) + | CGLO_SCALEEKIN + ); + /* explanation of above: + a) We compute Ekin at the full time step + if 1) we are using the AveVel Ekin, and it's not the + initial step, or 2) if we are using AveEkin, but need the full + time step kinetic energy for the pressure (always true now, since we want accurate statistics). + b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in + EkinAveVel because it's needed for the pressure */ + wallcycle_start(wcycle, ewcUPDATE); + } + /* temperature scaling and pressure scaling to produce the extended variables at t+dt */ + if (!bInitStep) + { + if (bTrotter) { - wallcycle_stop(wcycle, ewcUPDATE); - compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, - wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot, - constr, NULL, FALSE, state->box, - top_global, &bSumEkinhOld, - cglo_flags - | CGLO_ENERGY - | (bTemp ? CGLO_TEMPERATURE : 0) - | (bPres ? CGLO_PRESSURE : 0) - | (bPres ? CGLO_CONSTRAINT : 0) - | (bStopCM ? CGLO_STOPCM : 0) - | ((iterate.bIterationActive) ? CGLO_ITERATE : 0) - | (bFirstIterate ? CGLO_FIRSTITERATE : 0) - | CGLO_SCALEEKIN - ); - /* explanation of above: - a) We compute Ekin at the full time step - if 1) we are using the AveVel Ekin, and it's not the - initial step, or 2) if we are using AveEkin, but need the full - time step kinetic energy for the pressure (always true now, since we want accurate statistics). - b) If we are using EkinAveEkin for the kinetic energy for the temperature control, we still feed in - EkinAveVel because it's needed for the pressure */ - wallcycle_start(wcycle, ewcUPDATE); + m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */ + trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2); } - /* temperature scaling and pressure scaling to produce the extended variables at t+dt */ - if (!bInitStep) + else { - if (bTrotter) - { - m_add(force_vir, shake_vir, total_vir); /* we need the un-dispersion corrected total vir here */ - trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ2); - } - else + if (bExchanged) { - if (bExchanged) - { - wallcycle_stop(wcycle, ewcUPDATE); - /* We need the kinetic energy at minus the half step for determining - * the full step kinetic energy and possibly for T-coupling.*/ - /* This may not be quite working correctly yet . . . . */ - compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, - wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot, - constr, NULL, FALSE, state->box, - top_global, &bSumEkinhOld, - CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE); - wallcycle_start(wcycle, ewcUPDATE); - } + wallcycle_stop(wcycle, ewcUPDATE); + /* We need the kinetic energy at minus the half step for determining + * the full step kinetic energy and possibly for T-coupling.*/ + /* This may not be quite working correctly yet . . . . */ + compute_globals(fplog, gstat, cr, ir, fr, ekind, state, state_global, mdatoms, nrnb, vcm, + wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot, + constr, NULL, FALSE, state->box, + top_global, &bSumEkinhOld, + CGLO_RERUNMD | CGLO_GSTAT | CGLO_TEMPERATURE); + wallcycle_start(wcycle, ewcUPDATE); } } - - if (iterate.bIterationActive && - done_iterating(cr, fplog, step, &iterate, bFirstIterate, - state->veta, &vetanew)) - { - break; - } - bFirstIterate = FALSE; } - if (bTrotter && !bInitStep) { copy_mat(shake_vir, state->svir_prev);