########################################################################
include(gmxManageMPI)
+########################################################################
+#Process MiMiC settings
+########################################################################
+include(gmxManageMimic)
########################################################################
#Process shared/static library settings
gpu_utils -> hardware
topology -> listed-forces
listed-forces -> mdlib
+topology -> gmxpreprocess
\ No newline at end of file
* ``-DGMX_FFT_LIBRARY=xxx`` to select whether to use ``fftw3``, ``mkl`` or ``fftpack`` libraries for `FFT support`_
* ``-DCMAKE_BUILD_TYPE=Debug`` to build |Gromacs| in debug mode
+Building with MiMiC QM/MM support
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+
+MiMiC QM/MM interface integration will require linking against MiMiC
+communication library, that establishes the communication channel between
+|Gromacs| and CPMD. Check that the installation folder of the library
+is added to CMAKE_PREFIX_PATH if it is installed in non-standard location.
+Building QM/MM-capable version requires double-precision version of |Gromacs|
+compiled with MPI support:
+
+* ``-DGMX_DOUBLE=ON -DGMX_MPI -DGMX_MIMIC=ON``
+
Building older versions
^^^^^^^^^^^^^^^^^^^^^^^
:mdp-value:`integrator=tpic` gives identical results to
single-rank :mdp-value:`integrator=tpic`.
+ .. mdp-value:: mimic
+
+ Enable MiMiC QM/MM coupling to run hybrid molecular dynamics.
+ Keey in mind that its required to launch CPMD compiled with MiMiC as well.
+ In this mode all options regarding integration (T-coupling, P-coupling,
+ timestep and number of steps) are ignored as CPMD will do the integration
+ instead. Options related to forces computation (cutoffs, PME parameters,
+ etc.) are working as usual. Atom selection to define QM atoms is read
+ from :mdp:`QMMM-grps`
+
.. mdp:: tinit
(0) [ps]
.. mdp:: QMMM-grps
- groups to be descibed at the QM level
+ groups to be descibed at the QM level (works also in case of MiMiC QM/MM)
.. mdp:: QMMMscheme
add_subdirectory(simd)
add_subdirectory(imd)
add_subdirectory(compat)
+add_subdirectory(mimic)
if (NOT GMX_BUILD_MDRUN_ONLY)
add_subdirectory(gmxana)
add_subdirectory(gmxpreprocess)
}
/*! \brief Set the exclusion data for i-zone \p iz */
-static void make_exclusions_zone(gmx_domdec_t *dd,
- gmx_domdec_zones_t *zones,
+static void make_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
const std::vector<gmx_moltype_t> &moltype,
- const int *cginfo,
- t_blocka *lexcls,
- int iz,
- int at_start, int at_end)
+ const int *cginfo, t_blocka *lexcls, int iz,
+ int at_start, int at_end,
+ const gmx::ArrayRef<const int> intermolecularExclusionGroup)
{
int n_excl_at_max, n, at;
lexcls->nalloc_a = over_alloc_large(n + 1000);
srenew(lexcls->a, lexcls->nalloc_a);
}
+
if (GET_CGINFO_EXCL_INTER(cginfo[at]))
{
int a_gl, mb, mt, mol, a_mol, j;
/* We don't need exclusions for this atom */
lexcls->index[at] = n;
}
+
+ bool isExcludedAtom = !intermolecularExclusionGroup.empty() &&
+ std::find(intermolecularExclusionGroup.begin(),
+ intermolecularExclusionGroup.end(),
+ dd->globalAtomIndices[at]) !=
+ intermolecularExclusionGroup.end();
+
+ if (isExcludedAtom)
+ {
+ if (n + intermolecularExclusionGroup.size() > lexcls->nalloc_a)
+ {
+ lexcls->nalloc_a =
+ over_alloc_large(n + intermolecularExclusionGroup.size());
+ srenew(lexcls->a, lexcls->nalloc_a);
+ }
+ for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
+ {
+ if (const auto *entry = dd->ga2la->find(qmAtomGlobalIndex))
+ {
+ lexcls->a[n++] = entry->la;
+ }
+ }
+ }
}
lexcls->index[lexcls->nr] = n;
!rt->bExclRequired)
{
/* No charge groups and no distance check required */
- make_exclusions_zone(dd, zones,
- mtop->moltype, cginfo,
- excl_t,
- izone,
- cg0t, cg1t);
+ make_exclusions_zone(dd, zones, mtop->moltype, cginfo,
+ excl_t, izone, cg0t,
+ cg1t,
+ mtop->intermolecularExclusions);
}
else
{
multisim.cpp
replicaexchange.cpp
rerun.cpp
+ mimic.cpp
runner.cpp
simulationcontext.cpp
tpi.cpp
do_md();
}
break;
+ case eiMimic:
+ if (doRerun)
+ {
+ do_rerun();
+ }
+ else
+ {
+ do_mimic();
+ }
+ break;
case eiSteep:
do_steep();
break;
IntegratorFunctionType do_nm;
//! Implements test particle insertion
IntegratorFunctionType do_tpi;
+ //! Implements MiMiC QM/MM workflow
+ IntegratorFunctionType do_mimic;
/*! \brief Function to run the correct IntegratorFunctionType,
* based on the .mdp integrator field. */
void run(unsigned int ei, bool doRerun);
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief Declares the loop for MiMiC QM/MM
+ *
+ * \author Viacheslav Bolnykh <v.bolnykh@hpc-leap.eu>
+ * \ingroup module_mdrun
+ */
+#include "gmxpre.h"
+
+#include <cinttypes>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+
+#include <algorithm>
+#include <memory>
+
+#include "gromacs/awh/awh.h"
+#include "gromacs/commandline/filenm.h"
+#include "gromacs/compat/make_unique.h"
+#include "gromacs/domdec/collect.h"
+#include "gromacs/domdec/domdec.h"
+#include "gromacs/domdec/domdec_network.h"
+#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/domdec/partition.h"
+#include "gromacs/essentialdynamics/edsam.h"
+#include "gromacs/ewald/pme.h"
+#include "gromacs/ewald/pme-load-balancing.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/gmxlib/nrnb.h"
+#include "gromacs/gpu_utils/gpu_utils.h"
+#include "gromacs/imd/imd.h"
+#include "gromacs/listed-forces/manage-threading.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/checkpointhandler.h"
+#include "gromacs/mdlib/compute_io.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/mdlib/ebin.h"
+#include "gromacs/mdlib/expanded.h"
+#include "gromacs/mdlib/force.h"
+#include "gromacs/mdlib/force_flags.h"
+#include "gromacs/mdlib/forcerec.h"
+#include "gromacs/mdlib/md_support.h"
+#include "gromacs/mdlib/mdatoms.h"
+#include "gromacs/mdlib/mdebin.h"
+#include "gromacs/mdlib/mdoutf.h"
+#include "gromacs/mdlib/mdrun.h"
+#include "gromacs/mdlib/mdsetup.h"
+#include "gromacs/mdlib/membed.h"
+#include "gromacs/mdlib/nb_verlet.h"
+#include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
+#include "gromacs/mdlib/ns.h"
+#include "gromacs/mdlib/resethandler.h"
+#include "gromacs/mdlib/shellfc.h"
+#include "gromacs/mdlib/sighandler.h"
+#include "gromacs/mdlib/sim_util.h"
+#include "gromacs/mdlib/simulationsignal.h"
+#include "gromacs/mdlib/stophandler.h"
+#include "gromacs/mdlib/tgroup.h"
+#include "gromacs/mdlib/trajectory_writing.h"
+#include "gromacs/mdlib/update.h"
+#include "gromacs/mdlib/vcm.h"
+#include "gromacs/mdlib/vsite.h"
+#include "gromacs/mdtypes/awh-history.h"
+#include "gromacs/mdtypes/awh-params.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/df_history.h"
+#include "gromacs/mdtypes/energyhistory.h"
+#include "gromacs/mdtypes/fcdata.h"
+#include "gromacs/mdtypes/forcerec.h"
+#include "gromacs/mdtypes/group.h"
+#include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/mdtypes/interaction_const.h"
+#include "gromacs/mdtypes/md_enums.h"
+#include "gromacs/mdtypes/mdatom.h"
+#include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/state.h"
+#include "gromacs/mimic/MimicCommunicator.h"
+#include "gromacs/mimic/MimicUtils.h"
+#include "gromacs/pbcutil/mshift.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pulling/pull.h"
+#include "gromacs/swap/swapcoords.h"
+#include "gromacs/timing/wallcycle.h"
+#include "gromacs/timing/walltime_accounting.h"
+#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/basedefinitions.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/logger.h"
+#include "gromacs/utility/real.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "integrator.h"
+#include "replicaexchange.h"
+
+using gmx::SimulationSignaller;
+
+void gmx::Integrator::do_mimic()
+{
+ t_inputrec *ir = inputrec;
+ gmx_mdoutf *outf = nullptr;
+ int64_t step, step_rel;
+ double t, lam0[efptNR];
+ bool isLastStep = false;
+ bool doFreeEnergyPerturbation = false;
+ int force_flags;
+ tensor force_vir, shake_vir, total_vir, pres;
+ rvec mu_tot;
+ gmx_localtop_t *top;
+ t_mdebin *mdebin = nullptr;
+ gmx_enerdata_t *enerd;
+ PaddedVector<gmx::RVec> f {};
+ gmx_global_stat_t gstat;
+ t_graph *graph = nullptr;
+ gmx_groups_t *groups;
+ gmx_ekindata_t *ekind;
+ gmx_shellfc_t *shellfc;
+
+ double cycles;
+
+ /* Domain decomposition could incorrectly miss a bonded
+ interaction, but checking for that requires a global
+ communication stage, which does not otherwise happen in DD
+ code. So we do that alongside the first global energy reduction
+ after a new DD is made. These variables handle whether the
+ check happens, and the result it returns. */
+ bool shouldCheckNumberOfBondedInteractions = false;
+ int totalNumberOfBondedInteractions = -1;
+
+ SimulationSignals signals;
+ // Most global communnication stages don't propagate mdrun
+ // signals, and will use this object to achieve that.
+ SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
+
+ if (ir->bExpanded)
+ {
+ gmx_fatal(FARGS, "Expanded ensemble not supported by MiMiC.");
+ }
+ if (ir->bSimTemp)
+ {
+ gmx_fatal(FARGS, "Simulated tempering not supported by MiMiC.");
+ }
+ if (ir->bDoAwh)
+ {
+ gmx_fatal(FARGS, "AWH not supported by MiMiC.");
+ }
+ if (replExParams.exchangeInterval > 0)
+ {
+ gmx_fatal(FARGS, "Replica exchange not supported by MiMiC.");
+ }
+ if (opt2bSet("-ei", nfile, fnm) || observablesHistory->edsamHistory != nullptr)
+ {
+ gmx_fatal(FARGS, "Essential dynamics not supported by MiMiC.");
+ }
+ if (ir->bIMD)
+ {
+ gmx_fatal(FARGS, "Interactive MD not supported by MiMiC.");
+ }
+ if (isMultiSim(ms))
+ {
+ gmx_fatal(FARGS, "Multiple simulations not supported by MiMiC.");
+ }
+ if (std::any_of(ir->opts.annealing, ir->opts.annealing + ir->opts.ngtc,
+ [](int i){return i != eannNO; }))
+ {
+ gmx_fatal(FARGS, "Simulated annealing not supported by MiMiC.");
+ }
+
+ /* Settings for rerun */
+ ir->nstlist = 1;
+ ir->nstcalcenergy = 1;
+ int nstglobalcomm = 1;
+ const bool bNS = true;
+
+ // Communicator to interact with MiMiC
+ MimicCommunicator mimicCommunicator {};
+ if (MASTER(cr))
+ {
+ mimicCommunicator.init();
+ mimicCommunicator.sendInitData(top_global, state_global->x);
+ ir->nsteps = mimicCommunicator.getStepNumber();
+ }
+
+ ir->nstxout_compressed = 0;
+ groups = &top_global->groups;
+ top_global->intermolecularExclusions = genQmmmIndices(*top_global);
+
+ /* Initial values */
+ init_rerun(fplog, cr, outputProvider, ir, oenv, mdrunOptions,
+ state_global, lam0, nrnb, top_global,
+ nfile, fnm, &outf, &mdebin, wcycle);
+
+ /* Energy terms and groups */
+ snew(enerd, 1);
+ init_enerdata(top_global->groups.grps[egcENER].nr, ir->fepvals->n_lambda,
+ enerd);
+
+ /* Kinetic energy data */
+ snew(ekind, 1);
+ init_ekindata(fplog, top_global, &(ir->opts), ekind);
+ /* Copy the cos acceleration to the groups struct */
+ ekind->cosacc.cos_accel = ir->cos_accel;
+
+ gstat = global_stat_init(ir);
+
+ /* Check for polarizable models and flexible constraints */
+ shellfc = init_shell_flexcon(fplog,
+ top_global, constr ? constr->numFlexibleConstraints() : 0,
+ ir->nstcalcenergy, DOMAINDECOMP(cr));
+
+ {
+ double io = compute_io(ir, top_global->natoms, groups, mdebin->ebin->nener, 1);
+ if ((io > 2000) && MASTER(cr))
+ {
+ fprintf(stderr,
+ "\nWARNING: This run will generate roughly %.0f Mb of data\n\n",
+ io);
+ }
+ }
+
+ // Local state only becomes valid now.
+ std::unique_ptr<t_state> stateInstance;
+ t_state * state;
+
+ if (DOMAINDECOMP(cr))
+ {
+ top = dd_init_local_top(top_global);
+
+ stateInstance = compat::make_unique<t_state>();
+ state = stateInstance.get();
+ dd_init_local_state(cr->dd, state_global, state);
+
+ /* Distribute the charge groups over the nodes from the master node */
+ dd_partition_system(fplog, mdlog, ir->init_step, cr, TRUE, 1,
+ state_global, top_global, ir,
+ state, &f, mdAtoms, top, fr,
+ vsite, constr,
+ nrnb, nullptr, FALSE);
+ shouldCheckNumberOfBondedInteractions = true;
+ gmx_bcast(sizeof(ir->nsteps), &ir->nsteps, cr);
+ }
+ else
+ {
+ state_change_natoms(state_global, state_global->natoms);
+ /* We need to allocate one element extra, since we might use
+ * (unaligned) 4-wide SIMD loads to access rvec entries.
+ */
+ f.resizeWithPadding(state_global->natoms);
+ /* Copy the pointer to the global state */
+ state = state_global;
+
+ snew(top, 1);
+ mdAlgorithmsSetupAtomData(cr, ir, top_global, top, fr,
+ &graph, mdAtoms, constr, vsite, shellfc);
+ }
+
+ auto mdatoms = mdAtoms->mdatoms();
+
+ // NOTE: The global state is no longer used at this point.
+ // But state_global is still used as temporary storage space for writing
+ // the global state to file and potentially for replica exchange.
+ // (Global topology should persist.)
+
+ update_mdatoms(mdatoms, state->lambda[efptMASS]);
+
+ if (ir->efep != efepNO && ir->fepvals->nstdhdl != 0)
+ {
+ doFreeEnergyPerturbation = true;
+ }
+
+ {
+ int cglo_flags = (CGLO_INITIALIZATION | CGLO_GSTAT |
+ (shouldCheckNumberOfBondedInteractions ?
+ CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
+ bool bSumEkinhOld = false;
+ t_vcm *vcm = nullptr;
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ nullptr, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+ constr, &nullSignaller, state->box,
+ &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags);
+ }
+ checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
+
+ if (MASTER(cr))
+ {
+ fprintf(stderr, "starting MiMiC MD run '%s'\n\n",
+ *(top_global->name));
+ if (mdrunOptions.verbose)
+ {
+ fprintf(stderr, "Calculated time to finish depends on nsteps from "
+ "run input file,\nwhich may not correspond to the time "
+ "needed to process input trajectory.\n\n");
+ }
+ fprintf(fplog, "\n");
+ }
+
+ walltime_accounting_start_time(walltime_accounting);
+ wallcycle_start(wcycle, ewcRUN);
+ print_start(fplog, cr, walltime_accounting, "mdrun");
+
+ /***********************************************************
+ *
+ * Loop over MD steps
+ *
+ ************************************************************/
+
+ if (constr)
+ {
+ GMX_LOG(mdlog.info).asParagraph().
+ appendText("Simulations has constraints. Constraints will "
+ "be handled by CPMD.");
+ }
+
+ GMX_LOG(mdlog.info).asParagraph().
+ appendText("MiMiC does not report kinetic energy, total energy, temperature, virial and pressure.");
+
+ auto stopHandler = stopHandlerBuilder->getStopHandlerMD(
+ compat::not_null<SimulationSignal*>(&signals[eglsSTOPCOND]), false,
+ MASTER(cr), ir->nstlist, mdrunOptions.reproducible, nstglobalcomm,
+ mdrunOptions.maximumHoursToRun, ir->nstlist == 0, fplog, step, bNS, walltime_accounting);
+
+ // we don't do counter resetting in rerun - finish will always be valid
+ walltime_accounting_set_valid_finish(walltime_accounting);
+
+ DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion =
+ (DOMAINDECOMP(cr) ?
+ DdOpenBalanceRegionBeforeForceComputation::yes :
+ DdOpenBalanceRegionBeforeForceComputation::no);
+ DdCloseBalanceRegionAfterForceComputation ddCloseBalanceRegion =
+ (DOMAINDECOMP(cr) ?
+ DdCloseBalanceRegionAfterForceComputation::yes :
+ DdCloseBalanceRegionAfterForceComputation::no);
+
+ step = ir->init_step;
+ step_rel = 0;
+
+ /* and stop now if we should */
+ isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
+ while (!isLastStep)
+ {
+ isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel == ir->nsteps));
+ wallcycle_start(wcycle, ewcSTEP);
+
+ t = step;
+
+ if (MASTER(cr))
+ {
+ mimicCommunicator.getCoords(&state_global->x, state_global->natoms);
+ for (int i = 0; i < state_global->natoms; i++)
+ {
+ copy_rvec(state_global->x[i], state->x[i]);
+ }
+ }
+
+ if (ir->efep != efepNO)
+ {
+ setCurrentLambdasLocal(step, ir->fepvals, lam0, state);
+ }
+
+ if (MASTER(cr))
+ {
+ const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites);
+ if (constructVsites && DOMAINDECOMP(cr))
+ {
+ gmx_fatal(FARGS, "Vsite recalculation with -rerun is not implemented with domain decomposition, "
+ "use a single rank");
+ }
+ }
+
+ if (DOMAINDECOMP(cr))
+ {
+ /* Repartition the domain decomposition */
+ const bool bMasterState = true;
+ dd_partition_system(fplog, mdlog, step, cr,
+ bMasterState, nstglobalcomm,
+ state_global, top_global, ir,
+ state, &f, mdAtoms, top, fr,
+ vsite, constr,
+ nrnb, wcycle,
+ mdrunOptions.verbose);
+ shouldCheckNumberOfBondedInteractions = true;
+ }
+
+ if (MASTER(cr))
+ {
+ print_ebin_header(fplog, step, t); /* can we improve the information printed here? */
+ }
+
+ if (ir->efep != efepNO)
+ {
+ update_mdatoms(mdatoms, state->lambda[efptMASS]);
+ }
+
+ force_flags = (GMX_FORCE_STATECHANGED |
+ GMX_FORCE_DYNAMICBOX |
+ GMX_FORCE_ALLFORCES |
+ GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved
+ GMX_FORCE_ENERGY |
+ (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
+
+ if (shellfc)
+ {
+ /* Now is the time to relax the shells */
+ relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose,
+ enforcedRotation, step,
+ ir, bNS, force_flags, top,
+ constr, enerd, fcd,
+ state, f, force_vir, mdatoms,
+ nrnb, wcycle, graph, groups,
+ shellfc, fr, t, mu_tot,
+ vsite,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
+ }
+ else
+ {
+ /* The coordinates (x) are shifted (to get whole molecules)
+ * in do_force.
+ * This is parallellized as well, and does communication too.
+ * Check comments in sim_util.c
+ */
+ Awh *awh = nullptr;
+ gmx_edsam *ed = nullptr;
+ do_force(fplog, cr, ms, ir, awh, enforcedRotation,
+ step, nrnb, wcycle, top, groups,
+ state->box, state->x, &state->hist,
+ f, force_vir, mdatoms, enerd, fcd,
+ state->lambda, graph,
+ fr, vsite, mu_tot, t, ed,
+ GMX_FORCE_NS | force_flags,
+ ddOpenBalanceRegion, ddCloseBalanceRegion);
+ }
+
+ /* Now we have the energies and forces corresponding to the
+ * coordinates at time t.
+ */
+ {
+ const bool isCheckpointingStep = false;
+ const bool doRerun = false;
+ const bool bSumEkinhOld = false;
+ do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t,
+ ir, state, state_global, observablesHistory,
+ top_global, fr,
+ outf, mdebin, ekind, f,
+ isCheckpointingStep, doRerun, isLastStep,
+ mdrunOptions.writeConfout,
+ bSumEkinhOld);
+ }
+
+ stopHandler->setSignal();
+
+ if (graph)
+ {
+ /* Need to unshift here */
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
+ }
+
+ if (vsite != nullptr)
+ {
+ wallcycle_start(wcycle, ewcVSITECONSTR);
+ if (graph != nullptr)
+ {
+ shift_self(graph, state->box, as_rvec_array(state->x.data()));
+ }
+ construct_vsites(vsite, as_rvec_array(state->x.data()), ir->delta_t, as_rvec_array(state->v.data()),
+ top->idef.iparams, top->idef.il,
+ fr->ePBC, fr->bMolPBC, cr, state->box);
+
+ if (graph != nullptr)
+ {
+ unshift_self(graph, state->box, as_rvec_array(state->x.data()));
+ }
+ wallcycle_stop(wcycle, ewcVSITECONSTR);
+ }
+
+ {
+ const bool doInterSimSignal = false;
+ const bool doIntraSimSignal = true;
+ bool bSumEkinhOld = false;
+ t_vcm *vcm = nullptr;
+ SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
+
+ compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+ wcycle, enerd, nullptr, nullptr, nullptr, nullptr, mu_tot,
+ constr, &signaller,
+ state->box,
+ &totalNumberOfBondedInteractions, &bSumEkinhOld,
+ CGLO_GSTAT | CGLO_ENERGY
+ | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
+ );
+ checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
+ top_global, top, state,
+ &shouldCheckNumberOfBondedInteractions);
+ }
+
+ {
+ gmx::HostVector<gmx::RVec> fglobal(top_global->natoms);
+ gmx::ArrayRef<gmx::RVec> ftemp;
+ gmx::ArrayRef<const gmx::RVec> flocal = gmx::makeArrayRef(f);
+ if (DOMAINDECOMP(cr))
+ {
+ ftemp = gmx::makeArrayRef(fglobal);
+ dd_collect_vec(cr->dd, state, flocal, ftemp);
+ }
+ else
+ {
+ ftemp = gmx::makeArrayRef(f);
+ }
+
+ if (MASTER(cr))
+ {
+ mimicCommunicator.sendEnergies(enerd->term[F_EPOT]);
+ mimicCommunicator.sendForces(ftemp, state_global->natoms);
+ }
+ }
+
+
+
+ /* Note: this is OK, but there are some numerical precision issues with using the convergence of
+ the virial that should probably be addressed eventually. state->veta has better properies,
+ but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
+ generate the new shake_vir, but test the veta value for convergence. This will take some thought. */
+
+ if (ir->efep != efepNO)
+ {
+ /* Sum up the foreign energy and dhdl terms for md and sd.
+ Currently done every step so that dhdl is correct in the .edr */
+ sum_dhdl(enerd, state->lambda, ir->fepvals);
+ }
+
+ /* Output stuff */
+ if (MASTER(cr))
+ {
+ const bool bCalcEnerStep = true;
+ upd_mdebin(mdebin, doFreeEnergyPerturbation, bCalcEnerStep,
+ t, mdatoms->tmass, enerd, state,
+ ir->fepvals, ir->expandedvals, state->box,
+ shake_vir, force_vir, total_vir, pres,
+ ekind, mu_tot, constr);
+
+ const bool do_ene = true;
+ const bool do_log = true;
+ Awh *awh = nullptr;
+ const bool do_dr = ir->nstdisreout != 0;
+ const bool do_or = ir->nstorireout != 0;
+
+ print_ebin(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or, do_log ? fplog : nullptr,
+ step, t,
+ eprNORMAL, mdebin, fcd, groups, &(ir->opts), awh);
+
+ if (do_per_step(step, ir->nstlog))
+ {
+ if (fflush(fplog) != 0)
+ {
+ gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
+ }
+ }
+ }
+
+ /* Print the remaining wall clock time for the run */
+ if (isMasterSimMasterRank(ms, cr) &&
+ (mdrunOptions.verbose || gmx_got_usr_signal()))
+ {
+ if (shellfc)
+ {
+ fprintf(stderr, "\n");
+ }
+ print_time(stderr, walltime_accounting, step, ir, cr);
+ }
+
+ cycles = wallcycle_stop(wcycle, ewcSTEP);
+ if (DOMAINDECOMP(cr) && wcycle)
+ {
+ dd_cycles_add(cr->dd, cycles, ddCyclStep);
+ }
+
+ /* increase the MD step number */
+ step++;
+ step_rel++;
+ }
+ /* End of main MD loop */
+
+ /* Closing TNG files can include compressing data. Therefore it is good to do that
+ * before stopping the time measurements. */
+ mdoutf_tng_close(outf);
+
+ /* Stop measuring walltime */
+ walltime_accounting_end_time(walltime_accounting);
+
+ if (MASTER(cr))
+ {
+ mimicCommunicator.finalize();
+ }
+
+ if (!thisRankHasDuty(cr, DUTY_PME))
+ {
+ /* Tell the PME only node to finish */
+ gmx_pme_send_finish(cr);
+ }
+
+ done_mdebin(mdebin);
+ done_mdoutf(outf);
+
+ done_shellfc(fplog, shellfc, step_rel);
+
+ // Clean up swapcoords
+ if (ir->eSwapCoords != eswapNO)
+ {
+ finish_swapcoords(ir->swap);
+ }
+
+ walltime_accounting_set_nsteps_done(walltime_accounting, step_rel);
+
+ destroy_enerdata(enerd);
+ sfree(enerd);
+ sfree(top);
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ *
+ * \brief Declares the loop for MiMiC QM/MM
+ *
+ * \author Viacheslav Bolnykh <v.bolnykh@hpc-leap.eu>
+ * \ingroup module_mdrun
+ */
+#ifndef GMX_MIMIC_MIMIC_H
+#define GMX_MIMIC_MIMIC_H
+#include "integrator.h"
+
+namespace gmx
+{
+
+//! MD simulations
+integrator_t do_mimic;
+
+} // namespace gmx
+#endif //GMX_MIMIC_MIMIC_H
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "gromacs/mdtypes/state.h"
+#include "gromacs/mimic/MimicUtils.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
ir->nstxout_compressed = 0;
groups = &top_global->groups;
+ if (ir->eI == eiMimic)
+ {
+ top_global->intermolecularExclusions = genQmmmIndices(*top_global);
+ }
/* Initial values */
init_rerun(fplog, cr, outputProvider, ir, oenv, mdrunOptions,
--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2018, by the GROMACS development team, led by
+# Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+# and including many others, as listed in the AUTHORS file in the
+# top-level source directory and at http://www.gromacs.org.
+#
+# GROMACS is free software; you can redistribute it and/or
+# modify it under the terms of the GNU Lesser General Public License
+# as published by the Free Software Foundation; either version 2.1
+# of the License, or (at your option) any later version.
+#
+# GROMACS is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+# Lesser General Public License for more details.
+#
+# You should have received a copy of the GNU Lesser General Public
+# License along with GROMACS; if not, see
+# http://www.gnu.org/licenses, or write to the Free Software Foundation,
+# Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+# If you want to redistribute modifications to GROMACS, please
+# consider that scientific software is very special. Version
+# control is crucial - bugs must be traceable. We will be happy to
+# consider code for inclusion in the official distribution, but
+# derived work must not be called official GROMACS. Details are found
+# in the README & COPYING files - if they are missing, get the
+# official version at http://www.gromacs.org.
+#
+# To help us fund GROMACS development, we humbly ask that you cite
+# the research papers on the package. Check out http://www.gromacs.org.
+gmx_add_libgromacs_sources(
+ MimicCommunicator.cpp
+ MimicUtils.cpp
+)
+
+gmx_install_headers()
\ No newline at end of file
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include "MimicCommunicator.h"
+
+#include "config.h"
+
+#include <unordered_set>
+
+#include "gromacs/math/units.h"
+#include "gromacs/utility/fatalerror.h"
+
+#if GMX_MIMIC
+#include <DataTypes.h>
+#include <MessageApi.h>
+#endif
+
+#if !GMX_MIMIC
+//! \brief Definitions to stub the ones defined in DataTypes.h
+constexpr int TYPE_INT = 0, TYPE_DOUBLE = 0;
+
+/*! \brief Stub communication library function to call in case if
+ * GROMACS is compiled without MiMiC. Calling causes GROMACS to exit!
+ */
+static void MCL_init_client(const char *) // NOLINT(readability-named-parameter)
+{
+ gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+}
+
+/*! \brief Stub communication library function to call in case if
+ * GROMACS is compiled without MiMiC. Calling causes GROMACS to exit!
+ */
+static void MCL_send(void *, int, int, int) // NOLINT(readability-named-parameter)
+{
+ gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+}
+
+/*! \brief Stub communication library function to call in case if
+ * GROMACS is compiled without MiMiC. Calling causes GROMACS to exit!
+ */
+static void MCL_receive(void *, int, int, int) // NOLINT(readability-named-parameter)
+{
+ gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+}
+
+/*! \brief Stub communication library function to call in case if
+ * GROMACS is compiled without MiMiC. Calling causes GROMACS to exit!
+ */
+static void MCL_destroy() // NOLINT(readability-named-parameter)
+{
+ gmx_fatal(FARGS, "GROMACS is compiled without MiMiC support! Please, recompile with -DGMX_MIMIC=ON");
+}
+#endif
+
+void gmx::MimicCommunicator::init()
+{
+ char path[GMX_PATH_MAX];
+ gmx_getcwd(path, GMX_PATH_MAX);
+ return MCL_init_client(path);
+}
+
+void gmx::MimicCommunicator::sendInitData(gmx_mtop_t *mtop,
+ HostVector<gmx::RVec> coords)
+{
+ MCL_send(&mtop->natoms, 1, TYPE_INT, 0);
+ MCL_send(&mtop->atomtypes.nr, 1, TYPE_INT, 0);
+
+ std::vector<int> atomTypes;
+ std::vector<int> nAtomsMol;
+ std::vector<int> idOrder;
+ std::vector<double> charges;
+ std::vector<double> masses(mtop->atomtypes.nr, -1);
+ std::vector<int> elements(mtop->atomtypes.nr, -1);
+ std::vector<int> bonds;
+ std::vector<double> bondLengths;
+ std::unordered_set<int> existingTypes;
+
+ atomTypes.reserve(static_cast<size_t>(mtop->natoms));
+ charges.reserve(static_cast<size_t>(mtop->natoms));
+
+ int offset = 0;
+ for (const gmx_molblock_t &molblock : mtop->molblock)
+ {
+ gmx_moltype_t *type = &mtop->moltype[molblock.type];
+ for (int mol = 0; mol < molblock.nmol; ++mol)
+ {
+ int nconstr = type->ilist[F_CONSTR].size() / 3;
+ int nconstrc = type->ilist[F_CONSTRNC].size() / 3;
+ int nsettle = type->ilist[F_SETTLE].size() / 4;
+
+ for (int ncon = 0; ncon < nconstr + nconstrc; ++ncon)
+ {
+ int contype = type->ilist[F_CONSTR].iatoms[0];
+ int at1 = type->ilist[F_CONSTR].iatoms[1];
+ int at2 = type->ilist[F_CONSTR].iatoms[2];
+ bonds.push_back(offset + at1 + 1);
+ bonds.push_back(offset + at2 + 1);
+ bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA) / BOHR2NM);
+ }
+
+ for (int ncon = 0; ncon < nsettle; ++ncon)
+ {
+ t_iatom ox;
+ t_iatom h1;
+ t_iatom h2;
+
+ int contype = type->ilist[F_SETTLE].iatoms[0];
+
+ ox = type->ilist[F_SETTLE].iatoms[1];
+ h1 = type->ilist[F_SETTLE].iatoms[2];
+ h2 = type->ilist[F_SETTLE].iatoms[3];
+
+ bonds.push_back(offset + ox + 1);
+ bonds.push_back(offset + h1 + 1);
+
+ bonds.push_back(offset + ox + 1);
+ bonds.push_back(offset + h2 + 1);
+
+ bonds.push_back(offset + h1 + 1);
+ bonds.push_back(offset + h2 + 1);
+ bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA) / BOHR2NM);
+ bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA) / BOHR2NM);
+ bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dB) / BOHR2NM);
+ }
+
+ nAtomsMol.push_back(type->atoms.nr);
+ for (int at = 0; at < type->atoms.nr; ++at)
+ {
+ int atomtype = type->atoms.atom[at].type;
+ auto charge = static_cast<double>(type->atoms.atom[at].q);
+ idOrder.push_back(offset + 1);
+ offset++;
+ atomTypes.push_back(atomtype + 1);
+ charges.push_back(charge);
+ if (existingTypes.insert(atomtype).second)
+ {
+ masses[atomtype] = type->atoms.atom[at].m;
+ elements[atomtype] = type->atoms.atom[at].atomnumber;
+ }
+ }
+ }
+ }
+ // sending atom types
+ MCL_send(&*atomTypes.begin(), mtop->natoms, TYPE_INT, 0);
+
+ int max_multipole_order = 0;
+ // sending multipole orders
+ MCL_send(&max_multipole_order, 1, TYPE_INT, 0);
+
+ int nMolecules = nAtomsMol.size();
+ // sending molecule number
+ MCL_send(&nMolecules, 1, TYPE_INT, 0);
+
+ // sending number of atoms per molecules
+ MCL_send(&*nAtomsMol.begin(), nAtomsMol.size(), TYPE_INT, 0);
+
+ int nBonds = bonds.size() / 2;
+ // sending number of bond constraints
+ MCL_send(&nBonds, 1, TYPE_INT, 0);
+
+ // sending number of angle constraints
+ MCL_send(&max_multipole_order, 1, TYPE_INT, 0);
+
+ if (nBonds > 0)
+ {
+ // sending bonded atoms indices
+ MCL_send(&*bonds.begin(), bonds.size(), TYPE_INT, 0);
+
+ // sending bond lengths
+ MCL_send(&*bondLengths.begin(), bondLengths.size(), TYPE_DOUBLE, 0);
+ }
+
+ // sending array of atomic charges
+ MCL_send(&*charges.begin(), mtop->natoms, TYPE_DOUBLE, 0);
+
+ // sending array of atomic masses
+ MCL_send(&*masses.begin(), mtop->atomtypes.nr, TYPE_DOUBLE, 0);
+
+ // sending ids of atoms per molecule
+ MCL_send(&*idOrder.begin(), idOrder.size(), TYPE_INT, 0);
+
+ // sending list of elements
+ MCL_send(&*elements.begin(), mtop->atomtypes.nr, TYPE_INT, 0);
+
+ std::vector<double> convertedCoords;
+ for (auto &coord : coords)
+ {
+ convertedCoords.push_back(static_cast<double>(coord[0]) / BOHR2NM);
+ convertedCoords.push_back(static_cast<double>(coord[1]) / BOHR2NM);
+ convertedCoords.push_back(static_cast<double>(coord[2]) / BOHR2NM);
+ }
+
+ // sending array of coordinates
+ MCL_send(&*convertedCoords.begin(), 3 * mtop->natoms, TYPE_DOUBLE, 0);
+}
+
+int64_t gmx::MimicCommunicator::getStepNumber()
+{
+ int steps;
+ MCL_receive(&steps, 1, TYPE_INT, 0);
+ return steps;
+}
+
+void gmx::MimicCommunicator::getCoords(HostVector<RVec> *x, const int natoms)
+{
+ std::vector<double> coords(natoms * 3);
+ MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, 0);
+ for (int j = 0; j < natoms; ++j)
+ {
+ (*x)[j][0] = static_cast<real>(coords[j * 3] * BOHR2NM);
+ (*x)[j][1] = static_cast<real>(coords[j * 3 + 1] * BOHR2NM);
+ (*x)[j][2] = static_cast<real>(coords[j * 3 + 2] * BOHR2NM);
+ }
+}
+
+void gmx::MimicCommunicator::sendEnergies(real energy)
+{
+ double convertedEnergy = energy / (HARTREE2KJ * AVOGADRO);
+ MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, 0);
+}
+
+void gmx::MimicCommunicator::sendForces(gmx::ArrayRef<gmx::RVec> forces, int natoms)
+{
+ std::vector<double> convertedForce;
+ for (int j = 0; j < natoms; ++j)
+ {
+ convertedForce.push_back(static_cast<real>(forces[j][0]) / HARTREE_BOHR2MD);
+ convertedForce.push_back(static_cast<real>(forces[j][1]) / HARTREE_BOHR2MD);
+ convertedForce.push_back(static_cast<real>(forces[j][2]) / HARTREE_BOHR2MD);
+ }
+ MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, 0);
+}
+
+void gmx::MimicCommunicator::finalize()
+{
+ MCL_destroy();
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+#ifndef GMX_MIMIC_COMMUNICATOR_H
+#define GMX_MIMIC_COMMUNICATOR_H
+
+#include "gromacs/gpu_utils/hostallocator.h"
+#include "gromacs/math/paddedvector.h"
+#include "gromacs/mdlib/constr.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/futil.h"
+
+namespace gmx
+{
+/**
+ * \inlibraryapi
+ * \internal \brief
+ * Class-wrapper around MiMiC communication library
+ * It uses GROMACS' unit conversion to switch from GROMACS' units to a.u.
+ *
+ * \author Viacheslav Bolnykh <v.bolnykh@hpc-leap.eu>
+ * \ingroup module_mimic
+ */
+class MimicCommunicator
+{
+
+ public:
+ /*! \brief
+ * Initializes the communicator
+ */
+ void init();
+
+ /*! \brief
+ * Sends the data needed for MiMiC initialization
+ *
+ * That includes number of atoms, element numbers, charges, masses,
+ * maximal order of multipoles (0 for point-charge forcefields),
+ * number of molecules, number of atoms per each molecule,
+ * bond constraints data
+ *
+ * @param mtop global topology data
+ * @param coords coordinates of all atoms
+ */
+ void sendInitData(gmx_mtop_t *mtop,
+ HostVector<gmx::RVec> coords);
+
+ /*! \brief
+ * Gets the number of MD steps to perform from MiMiC
+ *
+ * @return nsteps the number of MD steps to perform
+ */
+ int64_t getStepNumber();
+
+ /*! \brief
+ * Receive and array of updated atomic coordinates from MiMiC
+ *
+ * @param x array of coordinates to fill
+ * @param natoms number of atoms in the system
+ */
+ void getCoords(HostVector<RVec> *x, int natoms);
+
+ /*! \brief
+ * Send the potential energy value to MiMiC
+ *
+ * @param energy energy value to send
+ */
+ void sendEnergies(real energy);
+
+ /*! \brief
+ * Send classical forces acting on all atoms in the system
+ * to MiMiC.
+ *
+ * @param forces array of forces to send
+ * @param natoms number of atoms in the system
+ */
+ void sendForces(ArrayRef<gmx::RVec> forces, int natoms);
+
+ /*! \brief
+ * Finish communications and disconnect from the server
+ */
+ void finalize();
+
+};
+
+} // namespace gmx
+
+#endif //GMX_MIMIC_COMMUNICATOR_H
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+#include "gmxpre.h"
+
+#include "MimicUtils.h"
+
+std::vector<int> genQmmmIndices(const gmx_mtop_t &mtop)
+{
+ std::vector<int> output;
+ int global_at = 0;
+ unsigned char *grpnr = mtop.groups.grpnr[egcQMMM];
+ for (const gmx_molblock_t &molb : mtop.molblock)
+ {
+ for (int mol = 0; mol < molb.nmol; ++mol)
+ {
+ for (int n_atom = 0; n_atom < mtop.moltype[molb.type].atoms.nr; ++n_atom)
+ {
+ if (!grpnr || !grpnr[global_at])
+ {
+ output.push_back(global_at);
+ }
+ ++global_at;
+ }
+ }
+ }
+ return output;
+}
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ * \brief Provides utility functions for MiMiC QM/MM
+ * \inlibraryapi
+ *
+ * \author Viacheslav Bolnykh <v.bolnykh@hpc-leap.eu>
+ * \ingroup module_mimic
+ */
+#ifndef GMX_MIMIC_MIMICUTILS_H
+#define GMX_MIMIC_MIMICUTILS_H
+
+#include <vector>
+
+#include "gromacs/gmxpreprocess/toppush.h"
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/topology.h"
+
+/*! \brief Generates the list of QM atoms
+ *
+ * This function generates vector containing global IDs of QM atoms
+ * based on the information stored in the QM/MM group (egcQMMM)
+ *
+ * \param[in] mtop Global topology object
+ * \return The list of global IDs of QM atoms
+ */
+std::vector<int> genQmmmIndices(const gmx_mtop_t &mtop);
+
+#endif //GMX_MIMIC_MIMICUTILS_H
#include <cstdlib>
#include <cstring>
+#include "gromacs/gmxpreprocess/toppush.h"
#include "gromacs/math/vectypes.h"
#include "gromacs/topology/atoms.h"
#include "gromacs/topology/block.h"
dest->nr += src->nr;
}
dest->index[dest->nr] = dest->nra;
+ dest->nalloc_index = dest->nr;
+ dest->nalloc_a = dest->nra;
}
static void ilistcat(int ftype,
}
}
-static void gen_local_top(const gmx_mtop_t &mtop,
- bool freeEnergyInteractionsAtEnd,
- bool bMergeConstr,
- gmx_localtop_t *top)
+/*! \brief Updates inter-molecular exclusion lists
+ *
+ * This function updates inter-molecular exclusions to exclude all
+ * non-bonded interactions between QM atoms
+ *
+ * \param[inout] excls existing exclusions in local topology
+ * \param[in] ids list of global QM atoms IDs
+ */
+static void addMimicExclusions(t_blocka *excls,
+ const gmx::ArrayRef<const int> ids)
+{
+ t_blocka inter_excl {};
+ init_blocka(&inter_excl);
+ size_t n_q = ids.size();
+
+ inter_excl.nr = excls->nr;
+ inter_excl.nra = n_q * n_q;
+
+ size_t total_nra = n_q * n_q;
+
+ snew(inter_excl.index, excls->nr + 1);
+ snew(inter_excl.a, total_nra);
+
+ for (int i = 0; i < excls->nr; ++i)
+ {
+ inter_excl.index[i] = 0;
+ }
+
+ /* Here we loop over the list of QM atom ids
+ * and create exclusions between all of them resulting in
+ * n_q * n_q sized exclusion list
+ */
+ int prev_index = 0;
+ for (int k = 0; k < inter_excl.nr; ++k)
+ {
+ inter_excl.index[k] = prev_index;
+ for (long i = 0; i < ids.size(); i++)
+ {
+ if (k != ids[i])
+ {
+ continue;
+ }
+ size_t index = n_q * i;
+ inter_excl.index[ids[i]] = index;
+ prev_index = index + n_q;
+ for (size_t j = 0; j < n_q; ++j)
+ {
+ inter_excl.a[n_q * i + j] = ids[j];
+ }
+ }
+ }
+ inter_excl.index[ids[n_q - 1] + 1] = n_q * n_q;
+
+ inter_excl.index[inter_excl.nr] = n_q * n_q;
+
+ t_block2 qmexcl2 {};
+ init_block2(&qmexcl2, excls->nr);
+ b_to_b2(&inter_excl, &qmexcl2);
+
+ // Merge the created exclusion list with the existing one
+ merge_excl(excls, &qmexcl2, nullptr);
+ done_block2(&qmexcl2);
+}
+
+static void gen_local_top(const gmx_mtop_t &mtop,
+ bool freeEnergyInteractionsAtEnd,
+ bool bMergeConstr,
+ gmx_localtop_t *top)
{
copyAtomtypesFromMtop(mtop, &top->atomtypes);
copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
copyCgsFromMtop(mtop, &top->cgs);
copyExclsFromMtop(mtop, &top->excls);
+ if (!mtop.intermolecularExclusions.empty())
+ {
+ addMimicExclusions(&top->excls,
+ mtop.intermolecularExclusions);
+ }
}
gmx_localtop_t *
struct t_atoms;
struct t_block;
struct t_symtab;
+enum struct GmxQmmmMode;
// TODO All of the functions taking a const gmx_mtop * are deprecated
// and should be replaced by versions taking const gmx_mtop & when
t_symtab symtab;
//! Tells whether we have valid molecule indices
bool haveMoleculeIndices = false;
+ /*! \brief List of global atom indices of atoms between which
+ * non-bonded interactions must be excluded.
+ */
+ std::vector<int> intermolecularExclusions;
/* Derived data below */
//! Indices for each molblock entry for fast lookup of atom properties
--- /dev/null
+[ System ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ Water ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ SOL ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ QMatoms ]
+ 4 5 6
--- /dev/null
+[ System ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ Water ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ SOL ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ QMatoms ]
+ 7 8 9 10 11 12
--- /dev/null
+Generated by gmx solvate
+ 12
+ 1SOL OW 1 0.768 1.144 1.023
+ 1SOL HW1 2 0.690 1.161 1.083
+ 1SOL HW2 3 0.802 1.231 0.987
+ 2SOL OW 4 0.830 1.273 1.422
+ 2SOL HW1 5 0.825 1.306 1.517
+ 2SOL HW2 6 0.744 1.292 1.376
+ 3SOL OW 7 0.822 1.002 1.372
+ 3SOL HW1 8 0.862 1.001 1.280
+ 3SOL HW2 9 0.832 1.094 1.412
+ 4SOL OW 10 0.859 0.956 0.861
+ 4SOL HW1 11 0.913 0.887 0.909
+ 4SOL HW2 12 0.827 1.025 0.927
+ 2.50000 2.50000 2.50000
--- /dev/null
+[ defaults ]
+
+; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
+1 2 yes 0.5 0.8333
+
+[ atomtypes ]
+; name at.num mass charge ptype sigma epsilon
+HW 1 1.008 0.0000 A 0.00000e+00 0.00000e+00
+OW 8 16.00 0.0000 A 3.15061e-01 6.36386e-01
+OW_spc 8 15.9994 0.0000 A 3.16557e-01 6.50629e-01
+HW_spc 1 1.0080 0.0000 A 0.00000e+00 0.00000e+00
+[ bondtypes ]
+; i j func b0 kb
+ OW HW 1 0.09572 462750.4 ; P water
+ HW HW 1 0.15136 462750.4 ; P water
+
+[ angletypes ]
+; i j k func th0 cth
+HW OW HW 1 104.520 836.800 ; TIP3P water
+HW HW OW 1 127.740 0.000 ; (found in crystallographic water with 3 bonds)
+
+; Include water topology
+[ moleculetype ]
+; molname nrexcl
+SOL 2
+
+[ atoms ]
+; id at type res nr res name at name cg nr charge mass
+ 1 OW_spc 1 SOL OW 1 -0.82 15.99940
+ 2 HW_spc 1 SOL HW1 1 0.41 1.00800
+ 3 HW_spc 1 SOL HW2 1 0.41 1.00800
+
+
+[ settles ]
+; OW funct doh dhh
+1 1 0.1 0.16330
+
+[ exclusions ]
+1 2 3
+2 1 3
+3 1 2
+
+[ system ]
+; Name
+Generated by gmx solvate
+
+[ molecules ]
+; Compound #mols
+SOL 4
\ No newline at end of file
tabulated_bonded_interactions.cpp
termination.cpp
trajectory_writing.cpp
+ mimic.cpp
# pseudo-library for code for testing mdrun
$<TARGET_OBJECTS:mdrun_test_objlib>
# pseudo-library for code for mdrun
# files with code for tests
domain_decomposition.cpp
minimize.cpp
+ mimic.cpp
multisim.cpp
multisimtest.cpp
pmetest.cpp
--- /dev/null
+Alanine dipeptide in water
+ 23
+ 1ALA N 1 2.367 0.023 0.092 -0.4721 -0.2209 -0.2428
+ 1ALA H1 2 2.342 0.122 0.074 3.8920 2.2407 2.1908
+ 1ALA H2 3 2.357 -0.028 0.008 0.1636 -1.6509 0.8158
+ 1ALA H3 4 2.311 -0.013 0.166 -0.2998 -0.4202 -0.5465
+ 1ALA CA 5 2.509 0.017 0.132 0.1102 0.3694 0.0487
+ 1ALA HA 6 2.516 -0.027 0.232 0.7745 -0.0467 -0.0173
+ 1ALA CB 7 2.589 -0.069 0.033 -0.0111 -0.1315 -0.1634
+ 1ALA HB1 8 2.548 -0.170 0.028 -0.0695 -0.0101 0.4311
+ 1ALA HB2 9 2.587 -0.025 -0.065 -0.0267 0.1572 0.5044
+ 1ALA HB3 10 2.693 -0.076 0.065 -0.2808 0.1818 0.4398
+ 1ALA C 11 2.567 0.160 0.139 0.5705 0.3686 0.1845
+ 1ALA O 12 2.496 0.257 0.112 -0.5448 -0.5389 -0.3907
+ 2ALA N 13 2.692 0.179 0.172 -0.4313 -0.2012 -0.0255
+ 2ALA H 14 2.757 0.107 0.201 1.4888 1.2812 -0.2476
+ 2ALA CA 15 2.747 0.315 0.177 0.4311 -0.1270 0.0374
+ 2ALA HA 16 2.742 0.360 0.079 0.6671 0.4015 0.4405
+ 2ALA CB 17 2.668 0.402 0.276 -0.0301 0.2343 0.3454
+ 2ALA HB1 18 2.566 0.417 0.256 0.5589 3.2748 2.1666
+ 2ALA HB2 19 2.676 0.360 0.375 0.5922 -0.3322 -0.4112
+ 2ALA HB3 20 2.708 0.502 0.276 0.7373 -0.5381 -0.4917
+ 2ALA C 21 2.897 0.313 0.218 1.0552 1.8657 0.1344
+ 2ALA OC1 22 2.959 0.421 0.229 -0.6645 -0.9445 -0.1410
+ 2ALA OC2 23 2.951 0.200 0.241 -0.1119 -0.3693 -0.0061
+ 2.50000 2.50000 2.50000
--- /dev/null
+[ System ]
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+ 16 17 18 19 20 21 22 23
+[ Protein ]
+ 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
+ 16 17 18 19 20 21 22 23
+[ Protein-H ]
+ 1 5 7 11 12 13 15 17 21 22 23
+[ C-alpha ]
+ 5 15
+[ Backbone ]
+ 1 5 11 13 15 21
+[ MainChain ]
+ 1 5 11 12 13 15 21 22 23
+[ MainChain+Cb ]
+ 1 5 7 11 12 13 15 17 21 22 23
+[ MainChain+H ]
+ 1 2 3 4 5 11 12 13 14 15 21 22 23
+[ SideChain ]
+ 6 7 8 9 10 16 17 18 19 20
+[ SideChain-H ]
+ 7 17
+[ QMatoms ]
+ 13 14 15 16 17 18 19 20 21 22 23
--- /dev/null
+;
+; File 'topol.top' was generated
+; By user: chaosit (1000)
+; On host: localhost.localdomain
+; At date: Tue May 16 18:44:37 2017
+;
+; This is a standalone topology file
+;
+; Created by:
+; :-) GROMACS - gmx pdb2gmx, 2017-dev-20170512-bd242e8-dirty-unknown (-:
+;
+; Executable: /home/chaosit/gromacs/bin/gmx_mpi
+; Data prefix: /home/chaosit/gromacs
+; Working dir: /home/chaosit/models/alala
+; Command line:
+; gmx_mpi pdb2gmx -f alaala.pdb
+; Force field was read from the standard GROMACS share directory.
+;
+
+; Include forcefield parameters
+#include "amber03.ff/forcefield.itp"
+
+[ moleculetype ]
+; Name nrexcl
+Protein_chain_A 3
+
+[ atoms ]
+; nr type resnr residue atom cgnr charge mass typeB chargeB massB
+; residue 1 ALA rtp NALA q +1.0
+ 1 N3 1 ALA N 1 -0.589266 14.01 ; qtot -0.5893
+ 2 H 1 ALA H1 2 0.446422 1.008 ; qtot -0.1428
+ 3 H 1 ALA H2 3 0.446422 1.008 ; qtot 0.3036
+ 4 H 1 ALA H3 4 0.446422 1.008 ; qtot 0.75
+ 5 CT 1 ALA CA 5 0.113871 12.01 ; qtot 0.8639
+ 6 HP 1 ALA HA 6 0.06715 1.008 ; qtot 0.931
+ 7 CT 1 ALA CB 7 -0.204113 12.01 ; qtot 0.7269
+ 8 HC 1 ALA HB1 8 0.063056 1.008 ; qtot 0.79
+ 9 HC 1 ALA HB2 9 0.063056 1.008 ; qtot 0.853
+ 10 HC 1 ALA HB3 10 0.063056 1.008 ; qtot 0.9161
+ 11 C 1 ALA C 11 0.676687 12.01 ; qtot 1.593
+ 12 O 1 ALA O 12 -0.592764 16 ; qtot 1
+; residue 2 ALA rtp CALA q -1.0
+ 13 N 2 ALA N 13 -0.622202 14.01 ; qtot 0.3778
+ 14 H 2 ALA H 14 0.359784 1.008 ; qtot 0.7376
+ 15 CT 2 ALA CA 15 -0.054102 12.01 ; qtot 0.6835
+ 16 H1 2 ALA HA 16 0.123283 1.008 ; qtot 0.8068
+ 17 CT 2 ALA CB 17 -0.185718 12.01 ; qtot 0.621
+ 18 HC 2 ALA HB1 18 0.069652 1.008 ; qtot 0.6907
+ 19 HC 2 ALA HB2 19 0.069652 1.008 ; qtot 0.7603
+ 20 HC 2 ALA HB3 20 0.069652 1.008 ; qtot 0.83
+ 21 C 2 ALA C 21 0.664014 12.01 ; qtot 1.494
+ 22 O2 2 ALA OC1 22 -0.747007 16 ; qtot 0.747
+ 23 O2 2 ALA OC2 23 -0.747007 16 ; qtot 0
+
+[ bonds ]
+; ai aj funct c0 c1 c2 c3
+ 1 2 1
+ 1 3 1
+ 1 4 1
+ 1 5 1
+ 5 6 1
+ 5 7 1
+ 5 11 1
+ 7 8 1
+ 7 9 1
+ 7 10 1
+ 11 12 1
+ 11 13 1
+ 13 14 1
+ 13 15 1
+ 15 16 1
+ 15 17 1
+ 15 21 1
+ 17 18 1
+ 17 19 1
+ 17 20 1
+ 21 22 1
+ 21 23 1
+
+[ pairs ]
+; ai aj funct c0 c1 c2 c3
+ 1 8 1
+ 1 9 1
+ 1 10 1
+ 1 12 1
+ 1 13 1
+ 2 6 1
+ 2 7 1
+ 2 11 1
+ 3 6 1
+ 3 7 1
+ 3 11 1
+ 4 6 1
+ 4 7 1
+ 4 11 1
+ 5 14 1
+ 5 15 1
+ 6 8 1
+ 6 9 1
+ 6 10 1
+ 6 12 1
+ 6 13 1
+ 7 12 1
+ 7 13 1
+ 8 11 1
+ 9 11 1
+ 10 11 1
+ 11 16 1
+ 11 17 1
+ 11 21 1
+ 12 14 1
+ 12 15 1
+ 13 18 1
+ 13 19 1
+ 13 20 1
+ 13 22 1
+ 13 23 1
+ 14 16 1
+ 14 17 1
+ 14 21 1
+ 16 18 1
+ 16 19 1
+ 16 20 1
+ 16 22 1
+ 16 23 1
+ 17 22 1
+ 17 23 1
+ 18 21 1
+ 19 21 1
+ 20 21 1
+
+[ angles ]
+; ai aj ak funct c0 c1 c2 c3
+ 2 1 3 1
+ 2 1 4 1
+ 2 1 5 1
+ 3 1 4 1
+ 3 1 5 1
+ 4 1 5 1
+ 1 5 6 1
+ 1 5 7 1
+ 1 5 11 1
+ 6 5 7 1
+ 6 5 11 1
+ 7 5 11 1
+ 5 7 8 1
+ 5 7 9 1
+ 5 7 10 1
+ 8 7 9 1
+ 8 7 10 1
+ 9 7 10 1
+ 5 11 12 1
+ 5 11 13 1
+ 12 11 13 1
+ 11 13 14 1
+ 11 13 15 1
+ 14 13 15 1
+ 13 15 16 1
+ 13 15 17 1
+ 13 15 21 1
+ 16 15 17 1
+ 16 15 21 1
+ 17 15 21 1
+ 15 17 18 1
+ 15 17 19 1
+ 15 17 20 1
+ 18 17 19 1
+ 18 17 20 1
+ 19 17 20 1
+ 15 21 22 1
+ 15 21 23 1
+ 22 21 23 1
+
+[ dihedrals ]
+; ai aj ak al funct c0 c1 c2 c3 c4 c5
+ 2 1 5 6 9
+ 2 1 5 7 9
+ 2 1 5 11 9
+ 3 1 5 6 9
+ 3 1 5 7 9
+ 3 1 5 11 9
+ 4 1 5 6 9
+ 4 1 5 7 9
+ 4 1 5 11 9
+ 1 5 7 8 9
+ 1 5 7 9 9
+ 1 5 7 10 9
+ 6 5 7 8 9
+ 6 5 7 9 9
+ 6 5 7 10 9
+ 11 5 7 8 9
+ 11 5 7 9 9
+ 11 5 7 10 9
+ 1 5 11 12 9
+ 1 5 11 13 9
+ 6 5 11 12 9
+ 6 5 11 13 9
+ 7 5 11 12 9
+ 7 5 11 13 9
+ 5 11 13 14 9
+ 5 11 13 15 9
+ 12 11 13 14 9
+ 12 11 13 15 9
+ 11 13 15 16 9
+ 11 13 15 17 9
+ 11 13 15 21 9
+ 14 13 15 16 9
+ 14 13 15 17 9
+ 14 13 15 21 9
+ 13 15 17 18 9
+ 13 15 17 19 9
+ 13 15 17 20 9
+ 16 15 17 18 9
+ 16 15 17 19 9
+ 16 15 17 20 9
+ 21 15 17 18 9
+ 21 15 17 19 9
+ 21 15 17 20 9
+ 13 15 21 22 9
+ 13 15 21 23 9
+ 16 15 21 22 9
+ 16 15 21 23 9
+ 17 15 21 22 9
+ 17 15 21 23 9
+
+[ dihedrals ]
+; ai aj ak al funct c0 c1 c2 c3
+ 5 13 11 12 4
+ 11 15 13 14 4
+ 15 22 21 23 4
+
+; Include Position restraint file
+#ifdef POSRES
+#include "posre.itp"
+#endif
+
+; Include water topology
+#include "amber03.ff/tip3p.itp"
+
+#ifdef POSRES_WATER
+; Position restraint for each water oxygen
+[ position_restraints ]
+; i funct fcx fcy fcz
+ 1 1 1000 1000 1000
+#endif
+
+; Include topology for ions
+#include "amber03.ff/ions.itp"
+
+[ system ]
+; Name
+UNNAMED in water
+
+[ molecules ]
+; Compound #mols
+Protein_chain_A 1
--- /dev/null
+[ System ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ Water ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ SOL ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
+[ QMatoms ]
+ 1 2 3 4 5 6 7 8 9 10 11 12
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2018, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \internal \file
+ * \brief
+ * Tests for MiMiC forces computation
+ *
+ * \author Viacheslav Bolnykh <v.bolnykh@hpc-leap.eu>
+ * \ingroup module_mdrun_integration_tests
+ */
+#include "gmxpre.h"
+
+#include <gtest/gtest.h>
+
+#include "gromacs/topology/ifunc.h"
+#include "gromacs/trajectory/energyframe.h"
+#include "gromacs/utility/basenetwork.h"
+#include "gromacs/utility/stringutil.h"
+
+#include "testutils/cmdlinetest.h"
+#include "testutils/refdata.h"
+
+#include "energycomparison.h"
+#include "energyreader.h"
+#include "moduletest.h"
+#include "simulationdatabase.h"
+
+namespace gmx
+{
+namespace test
+{
+
+//! Test fixture for bonded interactions
+class MimicTest : public gmx::test::MdrunTestFixture
+{
+ public:
+ //! Execute the trajectory writing test
+ void setupGrompp(const char *index_file, const char *top_file, const char *gro_file)
+ {
+ runner_.topFileName_ = TestFileManager::getInputFilePath(top_file);
+ runner_.groFileName_ = TestFileManager::getInputFilePath(gro_file);
+ runner_.ndxFileName_ = TestFileManager::getInputFilePath(index_file);
+ runner_.useStringAsMdpFile("integrator = mimic\n"
+ "QMMM-grps = QMatoms");
+ }
+ //! Prepare an mdrun caller
+ CommandLine setupMdrun()
+ {
+ CommandLine rerunCaller;
+ rerunCaller.append("mdrun");
+ rerunCaller.addOption("-rerun", runner_.groFileName_);
+ runner_.edrFileName_ = fileManager_.getTemporaryFilePath(".edr");
+ return rerunCaller;
+ }
+ //! Check the output of mdrun
+ void checkMdrun()
+ {
+ EnergyTolerances energiesToMatch
+ {{
+ {
+ interaction_function[F_EPOT].longname, relativeToleranceAsFloatingPoint(-20.1, 1e-4)
+ },
+ }};
+
+ TestReferenceData refData;
+ auto checker = refData.rootChecker();
+ checkEnergiesAgainstReferenceData(runner_.edrFileName_,
+ energiesToMatch,
+ &checker);
+ }
+};
+
+// This test checks if the energies produced with one quantum molecule are reasonable
+TEST_F(MimicTest, OneQuantumMol)
+{
+ setupGrompp("1quantum.ndx", "4water.top", "4water.gro");
+ if (gmx_node_rank() == 0)
+ {
+ EXPECT_EQ(0, runner_.callGromppOnThisRank());
+ }
+
+ test::CommandLine rerunCaller = setupMdrun();
+
+ int dummy = 0;
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ if (gmx_node_rank() == 0)
+ {
+ checkMdrun();
+ }
+ // Stupid way of synchronizining ranks but simplest API does not allow us to have any other
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+}
+
+// This test checks if the energies produced with all quantum molecules are reasonable (0)
+TEST_F(MimicTest, AllQuantumMol)
+{
+ setupGrompp("allquantum.ndx", "4water.top", "4water.gro");
+ if (gmx_node_rank() == 0)
+ {
+ EXPECT_EQ(0, runner_.callGromppOnThisRank());
+ }
+
+ test::CommandLine rerunCaller = setupMdrun();
+ int dummy = 0;
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ if (gmx_node_rank() == 0)
+ {
+ checkMdrun();
+ }
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+}
+
+// This test checks if the energies produced with two quantum molecules are reasonable
+// Needed to check the LJ intermolecular exclusions
+TEST_F(MimicTest, TwoQuantumMol)
+{
+ setupGrompp("2quantum.ndx", "4water.top", "4water.gro");
+ if (gmx_node_rank() == 0)
+ {
+ EXPECT_EQ(0, runner_.callGromppOnThisRank());
+ }
+
+ test::CommandLine rerunCaller = setupMdrun();
+ int dummy = 0;
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ if (gmx_node_rank() == 0)
+ {
+ checkMdrun();
+ }
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+}
+
+// This test checks if the energies produced with QM/MM boundary cutting the bond are ok
+TEST_F(MimicTest, BondCuts)
+{
+ setupGrompp("ala.ndx", "ala.top", "ala.gro");
+ if (gmx_node_rank() == 0)
+ {
+ EXPECT_EQ(0, runner_.callGromppOnThisRank());
+ }
+
+ test::CommandLine rerunCaller = setupMdrun();
+ int dummy = 0;
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+ ASSERT_EQ(0, runner_.callMdrun(rerunCaller));
+ if (gmx_node_rank() == 0)
+ {
+ checkMdrun();
+ }
+ if (gmx_mpi_initialized())
+ {
+ gmx_broadcast_world(sizeof(int), &dummy);
+ }
+}
+
+} // namespace test
+
+} // namespace gmx
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Energy Name="Potential">
+ <Real Name="Time 0.000000 Step 0 in frame 0">0.000000</Real>
+ </Energy>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Energy Name="Potential">
+ <Real Name="Time 0.000000 Step 0 in frame 0">153.635712</Real>
+ </Energy>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Energy Name="Potential">
+ <Real Name="Time 0.000000 Step 0 in frame 0">-21.228647</Real>
+ </Energy>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Energy Name="Potential">
+ <Real Name="Time 0.000000 Step 0 in frame 0">21.707289</Real>
+ </Energy>
+</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff