From: Mark Abraham Date: Tue, 12 Jan 2016 14:18:06 +0000 (+0100) Subject: Merge branch release-5-0 into release-5-1 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=4b9949f8c2fd84359dfe96ecaf966ad78ae5d9c6;p=alexxy%2Fgromacs.git Merge branch release-5-0 into release-5-1 Change-Id: I49ca6d8edc18c764eff368cd396873385f11ea9c --- 4b9949f8c2fd84359dfe96ecaf966ad78ae5d9c6 diff --cc src/programs/mdrun/md.cpp index a34678f561,3de0982b21..a065d8f9ef --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@@ -1347,117 -1445,114 +1347,112 @@@ double do_md(FILE *fplog, t_commrec *cr gs.sig[eglsCHKPT] = 1; } - /* at the start of step, randomize or scale the velocities (trotter done elsewhere) */ - if (EI_VV(ir->eI)) + /* at the start of step, randomize or scale the velocities ((if vv. Restriction of Andersen controlled + in preprocessing */ + + if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */ { - if (!bInitStep) + gmx_bool bIfRandomize; + bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr); + /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */ + if (constr && bIfRandomize) { - update_tcouple(step, ir, state, ekind, &MassQ, mdatoms); - } - if (ETC_ANDERSEN(ir->etc)) /* keep this outside of update_tcouple because of the extra info required to pass */ - { - gmx_bool bIfRandomize; - bIfRandomize = update_randomize_velocities(ir, step, cr, mdatoms, state, upd, constr); - /* if we have constraints, we have to remove the kinetic energy parallel to the bonds */ - if (constr && bIfRandomize) - { - update_constraints(fplog, step, NULL, ir, mdatoms, - state, fr->bMolPBC, graph, f, - &top->idef, tmp_vir, - cr, nrnb, wcycle, upd, constr, - TRUE, bCalcVir); - } - update_constraints(fplog, step, NULL, ir, ekind, mdatoms, ++ update_constraints(fplog, step, NULL, ir, mdatoms, + state, fr->bMolPBC, graph, f, + &top->idef, tmp_vir, + cr, nrnb, wcycle, upd, constr, - TRUE, bCalcVir, vetanew); ++ TRUE, bCalcVir); } } + /* ######### START SECOND UPDATE STEP ################# */ + /* Box is changed in update() when we do pressure coupling, + * but we should still use the old box for energy corrections and when + * writing it to the energy file, so it matches the trajectory files for + * the same timestep above. Make a copy in a separate array. + */ + copy_mat(state->box, lastbox); - if (bIterativeCase && do_per_step(step, ir->nstpcouple)) - { - gmx_iterate_init(&iterate, TRUE); - /* for iterations, we save these vectors, as we will be redoing the calculations */ - copy_coupling_state(state, bufstate, ekind, ekind_save, &(ir->opts)); - } + dvdl_constr = 0; - bFirstIterate = TRUE; - while (bFirstIterate || iterate.bIterationActive) + if (!bRerunMD || rerun_fr.bV || bForceUpdate) { - /* We now restore these vectors to redo the calculation with improved extended variables */ - if (iterate.bIterationActive) + wallcycle_start(wcycle, ewcUPDATE); + /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */ + if (bTrotter) { - copy_coupling_state(bufstate, state, ekind_save, ekind, &(ir->opts)); + trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3); + /* We can only do Berendsen coupling after we have summed + * the kinetic energy or virial. Since the happens + * in global_state after update, we should only do it at + * step % nstlist = 1 with bGStatEveryStep=FALSE. + */ + } + else + { + update_tcouple(step, ir, state, ekind, &MassQ, mdatoms); + update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep); } - /* We make the decision to break or not -after- the calculation of Ekin and Pressure, - so scroll down for that logic */ + if (bVV) + { + bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); - /* ######### START SECOND UPDATE STEP ################# */ - /* Box is changed in update() when we do pressure coupling, - * but we should still use the old box for energy corrections and when - * writing it to the energy file, so it matches the trajectory files for - * the same timestep above. Make a copy in a separate array. - */ - copy_mat(state->box, lastbox); + /* velocity half-step update */ + update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, + bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd, + ekind, M, upd, FALSE, etrtVELOCITY2, + cr, nrnb, constr, &top->idef); + } - bOK = TRUE; - dvdl_constr = 0; + /* Above, initialize just copies ekinh into ekin, + * it doesn't copy position (for VV), + * and entire integrator for MD. + */ - if (!(bRerunMD && !rerun_fr.bV && !bForceUpdate)) + if (ir->eI == eiVVAK) { - wallcycle_start(wcycle, ewcUPDATE); - /* UPDATE PRESSURE VARIABLES IN TROTTER FORMULATION WITH CONSTRAINTS */ - if (bTrotter) + /* We probably only need md->homenr, not state->natoms */ + if (state->natoms > cbuf_nalloc) { - if (iterate.bIterationActive) - { - if (bFirstIterate) - { - scalevir = 1; - } - else - { - /* we use a new value of scalevir to converge the iterations faster */ - scalevir = tracevir/trace(shake_vir); - } - msmul(shake_vir, scalevir, shake_vir); - m_add(force_vir, shake_vir, total_vir); - clear_mat(shake_vir); - } - trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ3); - /* We can only do Berendsen coupling after we have summed - * the kinetic energy or virial. Since the happens - * in global_state after update, we should only do it at - * step % nstlist = 1 with bGStatEveryStep=FALSE. - */ - } - else - { - update_tcouple(step, ir, state, ekind, &MassQ, mdatoms); - update_pcouple(fplog, step, ir, state, pcoupl_mu, M, bInitStep); + cbuf_nalloc = state->natoms; + srenew(cbuf, cbuf_nalloc); } + copy_rvecn(state->x, cbuf, 0, state->natoms); + } + bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); - if (bVV) - { - bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); + update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, + bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd, + ekind, M, upd, bInitStep, etrtPOSITION, cr, nrnb, constr, &top->idef); + wallcycle_stop(wcycle, ewcUPDATE); - /* velocity half-step update */ - update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f, - bUpdateDoLR, fr->f_twin, bCalcVir ? &fr->vir_twin_constr : NULL, fcd, - ekind, M, upd, FALSE, etrtVELOCITY2, - cr, nrnb, constr, &top->idef); - } + update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state, + fr->bMolPBC, graph, f, + &top->idef, shake_vir, + cr, nrnb, wcycle, upd, constr, + FALSE, bCalcVir); - /* Above, initialize just copies ekinh into ekin, - * it doesn't copy position (for VV), - * and entire integrator for MD. - */ + if (bCalcVir && bUpdateDoLR && ir->nstcalclr > 1) + { + /* Correct the virial for multiple time stepping */ + m_sub(shake_vir, fr->vir_twin_constr, shake_vir); + } + + if (ir->eI == eiVVAK) + { + /* erase F_EKIN and F_TEMP here? */ + /* just compute the kinetic energy at the half step to perform a trotter step */ + 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, lastbox, + top_global, &bSumEkinhOld, + cglo_flags | CGLO_TEMPERATURE + ); + wallcycle_start(wcycle, ewcUPDATE); + trotter_update(ir, step, ekind, enerd, state, total_vir, mdatoms, &MassQ, trotter_seq, ettTSEQ4); + /* now we know the scaling, we can compute the positions again again */ + copy_rvecn(cbuf, state->x, 0, state->natoms); - if (ir->eI == eiVVAK) - { - /* We probably only need md->homenr, not state->natoms */ - if (state->natoms > cbuf_nalloc) - { - cbuf_nalloc = state->natoms; - srenew(cbuf, cbuf_nalloc); - } - copy_rvecn(state->x, cbuf, 0, state->natoms); - } bUpdateDoLR = (fr->bTwinRange && do_per_step(step, ir->nstcalclr)); update_coords(fplog, step, ir, mdatoms, state, fr->bMolPBC, f,