From 675bfc76ab7febe514af377490f1368e94b79c56 Mon Sep 17 00:00:00 2001 From: Carsten Kutzner Date: Tue, 1 Mar 2016 15:43:13 +0100 Subject: [PATCH] Moved essential dynamics initialization call into do_md() - init_edsam is now called in do_md() in the same code block where the membed and the swap code gets initialized - the ed struct is removed from the integrator calls - removed a TODO, which was actually already done (closing of edo file) - Fixes one of the TODOs left from gerrit.gromacs.org/#/c/5660/ Change-Id: I396cf943fe34ef1683204effd4d5e935e6132dd7 --- src/gromacs/essentialdynamics/edsam.cpp | 63 ++++++++++++++++++++++----------- src/gromacs/essentialdynamics/edsam.h | 50 ++++++++++++-------------- src/gromacs/mdlib/constr.cpp | 22 ++++++------ src/gromacs/mdlib/constr.h | 7 ++-- src/gromacs/mdlib/integrator.h | 2 -- src/gromacs/mdlib/minimize.cpp | 4 --- src/gromacs/mdlib/tpi.cpp | 1 - src/programs/mdrun/md.cpp | 18 ++++++++-- src/programs/mdrun/membed.cpp | 1 - src/programs/mdrun/runner.cpp | 23 ++++-------- 10 files changed, 103 insertions(+), 88 deletions(-) diff --git a/src/gromacs/essentialdynamics/edsam.cpp b/src/gromacs/essentialdynamics/edsam.cpp index 8abd08b2e3..911c648456 100644 --- a/src/gromacs/essentialdynamics/edsam.cpp +++ b/src/gromacs/essentialdynamics/edsam.cpp @@ -54,6 +54,7 @@ #include "gromacs/math/vec.h" #include "gromacs/math/vectypes.h" #include "gromacs/mdlib/broadcaststructs.h" +#include "gromacs/mdlib/constr.h" #include "gromacs/mdlib/groupcoord.h" #include "gromacs/mdlib/mdrun.h" #include "gromacs/mdlib/sim_util.h" @@ -1182,12 +1183,23 @@ static void get_flood_energies(t_edpar *edi, real Vfl[], int nnames) #endif -gmx_edsam_t ed_open(int natoms, ObservablesHistory *oh, int nfile, const t_filenm fnm[], unsigned long Flags, const gmx_output_env_t *oenv, t_commrec *cr) +/* This function opens the ED input and output files, reads in all datasets it finds + * in the input file, and cross-checks whether the .edi file information is consistent + * with the essential dynamics data found in the checkpoint file (if present). + * gmx make_edi can be used to create an .edi input file. + */ +static gmx_edsam_t ed_open( + int natoms, + ObservablesHistory *oh, + const char *ediFileName, + const char *edoFileName, + gmx_bool bAppend, + const gmx_output_env_t *oenv, + t_commrec *cr) { gmx_edsam_t ed; int nED; - /* Allocate space for the ED data structure */ snew(ed, 1); @@ -1196,7 +1208,6 @@ gmx_edsam_t ed_open(int natoms, ObservablesHistory *oh, int nfile, const t_filen if (MASTER(cr)) { - fprintf(stderr, "ED sampling will be performed!\n"); snew(ed->edpar, 1); // If we start from a checkpoint file, we already have an edsamHistory struct @@ -1207,7 +1218,7 @@ gmx_edsam_t ed_open(int natoms, ObservablesHistory *oh, int nfile, const t_filen edsamhistory_t *EDstate = oh->edsamHistory.get(); /* Read the edi input file: */ - nED = read_edi_file(ftp2fn(efEDI, nfile, fnm), ed->edpar, natoms); + nED = read_edi_file(ediFileName, ed->edpar, natoms); /* Make sure the checkpoint was produced in a run using this .edi file */ if (EDstate->bFromCpt) @@ -1221,14 +1232,13 @@ gmx_edsam_t ed_open(int natoms, ObservablesHistory *oh, int nfile, const t_filen init_edsamstate(ed, EDstate); /* The master opens the ED output file */ - /* TODO This file is never closed... */ - if (Flags & MD_APPENDFILES) + if (bAppend) { - ed->edo = gmx_fio_fopen(opt2fn("-eo", nfile, fnm), "a+"); + ed->edo = gmx_fio_fopen(edoFileName, "a+"); } else { - ed->edo = xvgropen(opt2fn("-eo", nfile, fnm), + ed->edo = xvgropen(edoFileName, "Essential dynamics / flooding output", "Time (ps)", "RMSDs (nm), projections on EVs (nm), ...", oenv); @@ -2648,13 +2658,18 @@ static void write_edo_legend(gmx_edsam_t ed, int nED, const gmx_output_env_t *oe /* Init routine for ED and flooding. Calls init_edi in a loop for every .edi-cycle * contained in the input file, creates a NULL terminated list of t_edpar structures */ -void init_edsam(const gmx_mtop_t *mtop, - const t_inputrec *ir, - t_commrec *cr, - gmx_edsam_t ed, - rvec x[], - matrix box, - edsamhistory_t *EDstate) // Used on MASTER only +gmx_edsam_t init_edsam( + const char *ediFileName, + const char *edoFileName, + const gmx_mtop_t *mtop, + const t_inputrec *ir, + t_commrec *cr, + gmx_constr *constr, + rvec x[], + matrix box, + ObservablesHistory *oh, + const gmx_output_env_t *oenv, + gmx_bool bAppend) { t_edpar *edi = nullptr; /* points to a single edi data set */ int i, avindex; @@ -2665,17 +2680,22 @@ void init_edsam(const gmx_mtop_t *mtop, matrix fit_rotmat; /* ... and rotation from fit to reference structure */ rvec *ref_x_old = nullptr; /* helper pointer */ + if (MASTER(cr)) { fprintf(stderr, "ED: Initializing essential dynamics constraints.\n"); - if (nullptr == ed) + if (oh->edsamHistory != nullptr && oh->edsamHistory->bFromCpt && !gmx_fexist(ediFileName) ) { - gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or\n" - "flooding simulation. Please also provide the correct .edi file with -ei.\n"); + gmx_fatal(FARGS, "The checkpoint file you provided is from an essential dynamics or flooding\n" + "simulation. Please also set the .edi file on the command line with -ei.\n"); } } + /* Open input and output files, allocate space for ED data structure */ + gmx_edsam_t ed = ed_open(mtop->natoms, oh, ediFileName, edoFileName, bAppend, oenv, cr); + saveEdsamPointer(constr, ed); + /* Needed for initializing radacc radius in do_edsam */ ed->bFirst = TRUE; @@ -2701,6 +2721,8 @@ void init_edsam(const gmx_mtop_t *mtop, * not before dd_partition_system which is called after init_edsam */ if (MASTER(cr)) { + edsamhistory_t *EDstate = oh->edsamHistory.get(); + if (!EDstate->bFromCpt) { /* Remove PBC, make molecule(s) subject to ED whole. */ @@ -2917,10 +2939,9 @@ void init_edsam(const gmx_mtop_t *mtop, { /* In the single-CPU case, point the local atom numbers pointers to the global * one, so that we can use the same notation in serial and parallel case: */ - /* Loop over all ED data sets (usually only one, though) */ edi = ed->edpar; - for (int nr_edi = 1; nr_edi <= EDstate->nED; nr_edi++) + for (int nr_edi = 1; nr_edi <= nED; nr_edi++) { edi->sref.anrs_loc = edi->sref.anrs; edi->sav.anrs_loc = edi->sav.anrs; @@ -2997,6 +3018,8 @@ void init_edsam(const gmx_mtop_t *mtop, { fflush(ed->edo); } + + return ed; } diff --git a/src/gromacs/essentialdynamics/edsam.h b/src/gromacs/essentialdynamics/edsam.h index bd3a6b3e0c..5e882083c5 100644 --- a/src/gromacs/essentialdynamics/edsam.h +++ b/src/gromacs/essentialdynamics/edsam.h @@ -54,14 +54,14 @@ /*! \brief Abstract type for essential dynamics * - * The main type is defined only in edsam.c + * The main type is defined only in edsam.cpp */ typedef struct gmx_edsam *gmx_edsam_t; -struct edsamhistory_t; struct gmx_domdec_t; struct gmx_mtop_t; struct gmx_output_env_t; +struct ObservablesHistory; struct t_commrec; struct t_filenm; struct t_inputrec; @@ -80,40 +80,34 @@ void do_edsam(const t_inputrec *ir, gmx_int64_t step, t_commrec *cr, rvec xs[], rvec v[], matrix box, gmx_edsam_t ed); -/*! \brief Reads in the .edi file containing the essential dynamics (ED) and flooding data. - * - * This function opens the ED input and output files, reads in all datasets it finds - * in the input file, and cross-checks whether the .edi file information is consistent - * with the ED data found in the checkpoint file (if present). - * gmx make_edi can be used to create an .edi input file. - * - * \param natoms Number of atoms of the whole MD system. - * \param oh Observables history, contains ED observables history - * \param nfile Number of entries (files) in the fnm structure. - * \param fnm The filenames struct; it contains also the names of the - * essential dynamics and flooding in + output files. - * \param Flags Flags passed over from main, used to determine - * whether we are appending. - * \param oenv Needed to open the output xvgr file. - * \param cr Data needed for MPI communication. - * \returns Pointer to the initialized essential dynamics / flooding data. - */ -gmx_edsam_t ed_open(int natoms, struct ObservablesHistory *oh, int nfile, const t_filenm fnm[], - unsigned long Flags, const gmx_output_env_t *oenv, t_commrec *cr); - /*! \brief Initializes the essential dynamics and flooding module. * + * \param ediFileName Essential dynamics input file. + * \param edoFileName Output file for essential dynamics data. * \param mtop Molecular topology. * \param ir MD input parameter record. * \param cr Data needed for MPI communication. - * \param ed The essential dynamics data. + * \param constr Data structure keeping the constraint information. * \param x Positions of the whole MD system. * \param box The simulation box. - * \param EDstate ED data stored in the checkpoint file. + * \param oh The observables history container. + * \param oenv The output environment information. + * \param bAppend Append to existing output files? + * + * \returns A pointer to the ED data structure. */ -void init_edsam(const gmx_mtop_t *mtop, const t_inputrec *ir, t_commrec *cr, - gmx_edsam_t ed, rvec x[], matrix box, edsamhistory_t *EDstate); - +gmx_edsam_t init_edsam( + const char *ediFileName, + const char *edoFileName, + const gmx_mtop_t *mtop, + const t_inputrec *ir, + t_commrec *cr, + struct gmx_constr *constr, + rvec x[], + matrix box, + ObservablesHistory *oh, + const gmx_output_env_t *oenv, + gmx_bool bAppend); /*! \brief Make a selection of the home atoms for the ED groups. * diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 377b4f733e..452c156ff0 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -629,7 +629,7 @@ gmx_bool constrain(FILE *fplog, gmx_bool bLog, gmx_bool bEner, } if (constr->ed && delta_step > 0) { - /* apply the essential dynamcs constraints here */ + /* apply the essential dynamics constraints here */ do_edsam(ir, step, cr, xprime, v, box, constr->ed); } } @@ -1148,7 +1148,7 @@ real constr_r_max(FILE *fplog, const gmx_mtop_t *mtop, const t_inputrec *ir) gmx_constr_t init_constraints(FILE *fplog, const gmx_mtop_t *mtop, const t_inputrec *ir, - gmx_edsam_t ed, edsamhistory_t *edsamHistory, t_state *state, + bool doEssentialDynamics, t_commrec *cr) { int nconstraints = @@ -1161,7 +1161,7 @@ gmx_constr_t init_constraints(FILE *fplog, if (nconstraints + nsettles == 0 && !(ir->bPull && pull_have_constraint(ir->pull_work)) && - ed == nullptr) + !doEssentialDynamics) { return nullptr; } @@ -1300,19 +1300,17 @@ gmx_constr_t init_constraints(FILE *fplog, constr->warncount_lincs = 0; constr->warncount_settle = 0; - /* Initialize the essential dynamics sampling. - * Put the pointer to the ED struct in constr */ - constr->ed = ed; - if (ed != nullptr || edsamHistory != nullptr) - { - init_edsam(mtop, ir, cr, ed, as_rvec_array(state->x.data()), state->box, edsamHistory); - } - - constr->warn_mtop = mtop; + constr->warn_mtop = mtop; return constr; } +/* Put a pointer to the essential dynamics constraints into the constr struct */ +void saveEdsamPointer(gmx_constr_t constr, gmx_edsam_t ed) +{ + constr->ed = ed; +} + const t_blocka *atom2constraints_moltype(gmx_constr_t constr) { return constr->at2con_mt; diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index e0b8a0cd87..fb3c8e8e4f 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -38,7 +38,6 @@ #ifndef GMX_MBLIB_CONSTR_H #define GMX_MBLIB_CONSTR_H -#include "gromacs/essentialdynamics/edsam.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" @@ -206,10 +205,14 @@ gmx_bool constrain(FILE *log, gmx_bool bLog, gmx_bool bEner, gmx_constr_t init_constraints(FILE *log, const gmx_mtop_t *mtop, const t_inputrec *ir, - gmx_edsam_t ed, edsamhistory_t *edsamHistory, t_state *state, + bool doEssentialDynamics, struct t_commrec *cr); /* Initialize constraints stuff */ +void saveEdsamPointer(gmx_constr_t constr, + struct gmx_edsam *ed); +/* Put a pointer to the essential dynamics constraints into the constr struct */ + void set_constraints(gmx_constr_t constr, gmx_localtop_t *top, const t_inputrec *ir, diff --git a/src/gromacs/mdlib/integrator.h b/src/gromacs/mdlib/integrator.h index f9b12ef952..f48e6e910b 100644 --- a/src/gromacs/mdlib/integrator.h +++ b/src/gromacs/mdlib/integrator.h @@ -91,7 +91,6 @@ class MDLogger; * \param[in] mdatoms Structure containing atom information * \param[in] nrnb Accounting for floating point operations * \param[in] wcycle Wall cycle timing information - * \param[in] ed Essential dynamics sampling information * \param[in] fr Force record with cut-off information and more * \param[in] repl_ex_nst How often we do replica exchange (in steps) * \param[in] repl_ex_nex How many replicas we have @@ -116,7 +115,6 @@ typedef double integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &md ObservablesHistory *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t ed, t_forcerec *fr, int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t gmx_unused * membed, diff --git a/src/gromacs/mdlib/minimize.cpp b/src/gromacs/mdlib/minimize.cpp index 55937087f7..10a5b19191 100644 --- a/src/gromacs/mdlib/minimize.cpp +++ b/src/gromacs/mdlib/minimize.cpp @@ -998,7 +998,6 @@ double do_cg(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog, ObservablesHistory *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t gmx_unused ed, t_forcerec *fr, int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed, gmx_membed_t gmx_unused *membed, @@ -1651,7 +1650,6 @@ double do_lbfgs(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlo ObservablesHistory *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t gmx_unused ed, t_forcerec *fr, int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed, gmx_membed_t gmx_unused *membed, @@ -2424,7 +2422,6 @@ double do_steep(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlo ObservablesHistory *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t gmx_unused ed, t_forcerec *fr, int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed, gmx_membed_t gmx_unused *membed, @@ -2695,7 +2692,6 @@ double do_nm(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog, ObservablesHistory gmx_unused *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t gmx_unused ed, t_forcerec *fr, int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed, gmx_membed_t gmx_unused *membed, diff --git a/src/gromacs/mdlib/tpi.cpp b/src/gromacs/mdlib/tpi.cpp index b365d98270..8ab7bb3221 100644 --- a/src/gromacs/mdlib/tpi.cpp +++ b/src/gromacs/mdlib/tpi.cpp @@ -160,7 +160,6 @@ double do_tpi(FILE *fplog, t_commrec *cr, const gmx::MDLogger gmx_unused &mdlog, ObservablesHistory gmx_unused *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t gmx_unused ed, t_forcerec *fr, int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed, gmx_membed_t gmx_unused *membed, diff --git a/src/programs/mdrun/md.cpp b/src/programs/mdrun/md.cpp index 5c110ab45a..b73092eaf8 100644 --- a/src/programs/mdrun/md.cpp +++ b/src/programs/mdrun/md.cpp @@ -53,6 +53,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/domdec_struct.h" +#include "gromacs/essentialdynamics/edsam.h" #include "gromacs/ewald/pme.h" #include "gromacs/ewald/pme-load-balancing.h" #include "gromacs/fileio/trxio.h" @@ -204,7 +205,6 @@ static void reset_all_counters(FILE *fplog, const gmx::MDLogger &mdlog, t_commre t_state *state_global, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t ed, t_forcerec *fr, int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, real cpt_period, real max_hours, @@ -225,7 +225,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog, ObservablesHistory *observablesHistory, t_mdatoms *mdatoms, t_nrnb *nrnb, gmx_wallcycle_t wcycle, - gmx_edsam_t ed, t_forcerec *fr, + t_forcerec *fr, int repl_ex_nst, int repl_ex_nex, int repl_ex_seed, gmx_membed_t *membed, real cpt_period, real max_hours, @@ -347,6 +347,17 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog, } groups = &top_global->groups; + gmx_edsam *ed = nullptr; + if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr) + { + /* Initialize essential dynamics sampling */ + ed = init_edsam(opt2fn_null("-ei", nfile, fnm), opt2fn("-eo", nfile, fnm), + top_global, ir, cr, constr, + as_rvec_array(state_global->x.data()), + state_global->box, observablesHistory, + oenv, Flags & MD_APPENDFILES); + } + if (ir->eSwapCoords != eswapNO) { /* Initialize ion swapping code */ @@ -1858,6 +1869,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog, finish_swapcoords(ir->swap); } + /* Do essential dynamics cleanup if needed. Close .edo file */ + done_ed(&ed); + /* IMD cleanup, if bIMD is TRUE. */ IMD_finalize(ir->bIMD, ir->imd); diff --git a/src/programs/mdrun/membed.cpp b/src/programs/mdrun/membed.cpp index 4a71a749f0..6a10a22149 100644 --- a/src/programs/mdrun/membed.cpp +++ b/src/programs/mdrun/membed.cpp @@ -40,7 +40,6 @@ #include #include "gromacs/commandline/filenm.h" -#include "gromacs/essentialdynamics/edsam.h" #include "gromacs/fileio/readinp.h" #include "gromacs/fileio/warninp.h" #include "gromacs/gmxlib/network.h" diff --git a/src/programs/mdrun/runner.cpp b/src/programs/mdrun/runner.cpp index b9984d63e7..ab159f87c0 100644 --- a/src/programs/mdrun/runner.cpp +++ b/src/programs/mdrun/runner.cpp @@ -57,7 +57,6 @@ #include "gromacs/commandline/filenm.h" #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" -#include "gromacs/essentialdynamics/edsam.h" #include "gromacs/ewald/pme.h" #include "gromacs/fileio/checkpoint.h" #include "gromacs/fileio/oenv.h" @@ -728,7 +727,6 @@ int mdrunner(gmx_hw_opt_t *hw_opt, gmx_walltime_accounting_t walltime_accounting = nullptr; int rc; gmx_int64_t reset_counters; - gmx_edsam_t ed = nullptr; int nthreads_pme = 1; gmx_membed_t * membed = nullptr; gmx_hw_info_t *hwinfo = nullptr; @@ -1050,16 +1048,6 @@ int mdrunner(gmx_hw_opt_t *hw_opt, gmx_bcast(sizeof(box), box, cr); } - // TODO This should move to do_md(), because it only makes sense - // with dynamical integrators, but there is no test coverage and - // it interacts with constraints, somehow. - /* Essential dynamics */ - if (opt2bSet("-ei", nfile, fnm)) - { - /* Open input and output files, allocate space for ED data structure */ - ed = ed_open(mtop->natoms, &observablesHistory, nfile, fnm, Flags, oenv, cr); - } - if (PAR(cr) && !(EI_TPI(inputrec->eI) || inputrec->eI == eiNM)) { @@ -1369,7 +1357,12 @@ int mdrunner(gmx_hw_opt_t *hw_opt, bVerbose, Flags); } - constr = init_constraints(fplog, mtop, inputrec, ed, observablesHistory.edsamHistory.get(), state, cr); + /* Let init_constraints know whether we have essential dynamics constraints. + * TODO: inputrec should tell us whether we use an algorithm, not a file option or the checkpoint + */ + bool doEdsam = (opt2fn_null("-ei", nfile, fnm) != nullptr || observablesHistory.edsamHistory); + + constr = init_constraints(fplog, mtop, inputrec, doEdsam, cr); if (DOMAINDECOMP(cr)) { @@ -1389,7 +1382,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt, nstepout, mdModules.outputProvider(), inputrec, mtop, fcd, state, &observablesHistory, - mdatoms, nrnb, wcycle, ed, fr, + mdatoms, nrnb, wcycle, fr, repl_ex_nst, repl_ex_nex, repl_ex_seed, membed, cpt_period, max_hours, @@ -1455,8 +1448,6 @@ int mdrunner(gmx_hw_opt_t *hw_opt, rc = (int)gmx_get_stop_condition(); - done_ed(&ed); - #if GMX_THREAD_MPI /* we need to join all threads. The sub-threads join when they exit this function, but the master thread needs to be told to -- 2.11.4.GIT