ON
GMX_MPI)
mark_as_advanced(GMX_MPI_IN_PLACE)
+
+option(GMX_MIMIC "Enable MiMiC QM/MM interface (CPMD is required)" OFF)
+
option(GMX_FAHCORE "Build a library with mdrun functionality" OFF)
mark_as_advanced(GMX_FAHCORE)
--- /dev/null
+#
+# This file is part of the GROMACS molecular simulation package.
+#
+# Copyright (c) 2017,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.
+
+if(GMX_MIMIC)
+ find_package(MimicCommLib REQUIRED)
+ include_directories(SYSTEM ${MimicCommLib_INCLUDE_DIRS})
+ list(APPEND GMX_COMMON_LIBRARIES Upstream::mimiccomm)
+endif()
/* Build using clang analyzer */
#cmakedefine01 GMX_CLANG_ANALYZER
+/* Use MiMiC QM/MM interface */
+#cmakedefine01 GMX_MIMIC
+
/*! \endcond */
#endif
tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */
tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */
tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */
+ tpxv_MimicQMMM, /**< Inroduced support for MiMiC QM/MM interface */
tpxv_Count /**< the total number of tpxv versions */
};
snew(ir->opts.SAoff, ir->opts.ngQM);
snew(ir->opts.SAsteps, ir->opts.ngQM);
}
- if (ir->opts.ngQM > 0)
+ if (ir->opts.ngQM > 0 && ir->bQMMM)
{
gmx_fio_ndo_int(fio, ir->opts.QMmethod, ir->opts.ngQM);
gmx_fio_ndo_int(fio, ir->opts.QMbasis, ir->opts.ngQM);
}
else
{
- generate_qmexcl(&sys, ir, wi);
+ generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL);
}
}
+ if (ir->eI == eiMimic)
+ {
+ generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC);
+ }
+
if (ftp2bSet(efTRN, NFILE, fnm))
{
if (bVerbose)
auto qmGroupNames = gmx::splitString(is->QMMM);
auto qmMethods = gmx::splitString(is->QMmethod);
auto qmBasisSets = gmx::splitString(is->QMbasis);
- if (qmMethods.size() != qmGroupNames.size() ||
- qmBasisSets.size() != qmGroupNames.size())
+ if (ir->eI != eiMimic)
+ {
+ if (qmMethods.size() != qmGroupNames.size() ||
+ qmBasisSets.size() != qmGroupNames.size())
+ {
+ gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
+ " and %zu methods\n", qmGroupNames.size(),
+ qmBasisSets.size(), qmMethods.size());
+ }
+ /* group rest, if any, is always MM! */
+ do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM,
+ restnm, egrptpALL_GENREST, bVerbose, wi);
+ nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
+ ir->opts.ngQM = qmGroupNames.size();
+ snew(ir->opts.QMmethod, nr);
+ snew(ir->opts.QMbasis, nr);
+ for (i = 0; i < nr; i++)
+ {
+ /* input consists of strings: RHF CASSCF PM3 .. These need to be
+ * converted to the corresponding enum in names.c
+ */
+ ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(),
+ eQMmethodNR,
+ eQMmethod_names);
+ ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(),
+ eQMbasisNR,
+ eQMbasis_names);
+
+ }
+ auto qmMultiplicities = gmx::splitString(is->QMmult);
+ auto qmCharges = gmx::splitString(is->QMcharge);
+ auto qmbSH = gmx::splitString(is->bSH);
+ snew(ir->opts.QMmult, nr);
+ snew(ir->opts.QMcharge, nr);
+ snew(ir->opts.bSH, nr);
+ convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
+ convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
+ convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
+
+ auto CASelectrons = gmx::splitString(is->CASelectrons);
+ auto CASorbitals = gmx::splitString(is->CASorbitals);
+ snew(ir->opts.CASelectrons, nr);
+ snew(ir->opts.CASorbitals, nr);
+ convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
+ convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
+
+ auto SAon = gmx::splitString(is->SAon);
+ auto SAoff = gmx::splitString(is->SAoff);
+ auto SAsteps = gmx::splitString(is->SAsteps);
+ snew(ir->opts.SAon, nr);
+ snew(ir->opts.SAoff, nr);
+ snew(ir->opts.SAsteps, nr);
+ convertInts(wi, SAon, "SAon", ir->opts.SAon);
+ convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
+ convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
+ }
+ else
{
- gmx_fatal(FARGS, "Invalid QMMM input: %zu groups %zu basissets"
- " and %zu methods\n", qmGroupNames.size(),
- qmBasisSets.size(), qmMethods.size());
+ /* MiMiC */
+ if (qmGroupNames.size() > 1)
+ {
+ gmx_fatal(FARGS, "Currently, having more than one QM group in MiMiC is not supported");
+ }
+ /* group rest, if any, is always MM! */
+ do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM,
+ restnm, egrptpALL_GENREST, bVerbose, wi);
+
+ ir->opts.ngQM = qmGroupNames.size();
}
- /* group rest, if any, is always MM! */
- do_numbering(natoms, groups, qmGroupNames, grps, gnames, egcQMMM,
- restnm, egrptpALL_GENREST, bVerbose, wi);
- nr = qmGroupNames.size(); /*atoms->grps[egcQMMM].nr;*/
- ir->opts.ngQM = qmGroupNames.size();
- snew(ir->opts.QMmethod, nr);
- snew(ir->opts.QMbasis, nr);
- for (i = 0; i < nr; i++)
- {
- /* input consists of strings: RHF CASSCF PM3 .. These need to be
- * converted to the corresponding enum in names.c
- */
- ir->opts.QMmethod[i] = search_QMstring(qmMethods[i].c_str(), eQMmethodNR,
- eQMmethod_names);
- ir->opts.QMbasis[i] = search_QMstring(qmBasisSets[i].c_str(), eQMbasisNR,
- eQMbasis_names);
-
- }
- auto qmMultiplicities = gmx::splitString(is->QMmult);
- auto qmCharges = gmx::splitString(is->QMcharge);
- auto qmbSH = gmx::splitString(is->bSH);
- snew(ir->opts.QMmult, nr);
- snew(ir->opts.QMcharge, nr);
- snew(ir->opts.bSH, nr);
- convertInts(wi, qmMultiplicities, "QMmult", ir->opts.QMmult);
- convertInts(wi, qmCharges, "QMcharge", ir->opts.QMcharge);
- convertYesNos(wi, qmbSH, "bSH", ir->opts.bSH);
-
- auto CASelectrons = gmx::splitString(is->CASelectrons);
- auto CASorbitals = gmx::splitString(is->CASorbitals);
- snew(ir->opts.CASelectrons, nr);
- snew(ir->opts.CASorbitals, nr);
- convertInts(wi, CASelectrons, "CASelectrons", ir->opts.CASelectrons);
- convertInts(wi, CASorbitals, "CASOrbitals", ir->opts.CASorbitals);
-
- auto SAon = gmx::splitString(is->SAon);
- auto SAoff = gmx::splitString(is->SAoff);
- auto SAsteps = gmx::splitString(is->SAsteps);
- snew(ir->opts.SAon, nr);
- snew(ir->opts.SAoff, nr);
- snew(ir->opts.SAsteps, nr);
- convertInts(wi, SAon, "SAon", ir->opts.SAon);
- convertInts(wi, SAoff, "SAoff", ir->opts.SAoff);
- convertInts(wi, SAsteps, "SAsteps", ir->opts.SAsteps);
/* end of QMMM input */
EXPECT_DEATH_IF_SUPPORTED(runTest(inputMdpFile), "Invalid enum");
}
+TEST_F(GetIrTest, HandlesMimic)
+{
+ const char *inputMdpFile[] = {"integrator = mimic", "QMMM-grps = QMatoms"};
+ runTest(joinStrings(inputMdpFile, "\n"));
+}
+
} // namespace test
} // namespace gmx
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Bool Name="Error parsing mdp file">false</Bool>
+ <String Name="OutputMdpFile">
+; VARIOUS PREPROCESSING OPTIONS
+; Preprocessor information: use cpp syntax.
+; e.g.: -I/home/joe/doe -I/home/mary/roe
+include =
+; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)
+define =
+
+; RUN CONTROL PARAMETERS
+integrator = mimic
+; Start time and timestep in ps
+tinit = 0
+dt = 0.001
+nsteps = 0
+; For exact run continuation or redoing part of a run
+init-step = 0
+; Part index is updated automatically on checkpointing (keeps files separate)
+simulation-part = 1
+; mode for center of mass motion removal
+comm-mode = Linear
+; number of steps for center of mass motion removal
+nstcomm = 100
+; group(s) for center of mass motion removal
+comm-grps =
+
+; LANGEVIN DYNAMICS OPTIONS
+; Friction coefficient (amu/ps) and random seed
+bd-fric = 0
+ld-seed = -1
+
+; ENERGY MINIMIZATION OPTIONS
+; Force tolerance and initial step-size
+emtol = 10
+emstep = 0.01
+; Max number of iterations in relax-shells
+niter = 20
+; Step size (ps^2) for minimization of flexible constraints
+fcstep = 0
+; Frequency of steepest descents steps when doing CG
+nstcgsteep = 1000
+nbfgscorr = 10
+
+; TEST PARTICLE INSERTION OPTIONS
+rtpi = 0.05
+
+; OUTPUT CONTROL OPTIONS
+; Output frequency for coords (x), velocities (v) and forces (f)
+nstxout = 0
+nstvout = 0
+nstfout = 0
+; Output frequency for energies to log file and energy file
+nstlog = 1000
+nstcalcenergy = 100
+nstenergy = 1000
+; Output frequency and precision for .xtc file
+nstxout-compressed = 0
+compressed-x-precision = 1000
+; This selects the subset of atoms for the compressed
+; trajectory file. You can select multiple groups. By
+; default, all atoms will be written.
+compressed-x-grps =
+; Selection of energy groups
+energygrps =
+
+; NEIGHBORSEARCHING PARAMETERS
+; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+cutoff-scheme = Verlet
+; nblist update frequency
+nstlist = 10
+; ns algorithm (simple or grid)
+ns-type = Grid
+; Periodic boundary conditions: xyz, no, xy
+pbc = xyz
+periodic-molecules = no
+; Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,
+; a value of -1 means: use rlist
+verlet-buffer-tolerance = 0.005
+; nblist cut-off
+rlist = 1
+; long-range cut-off for switched potentials
+
+; OPTIONS FOR ELECTROSTATICS AND VDW
+; Method for doing electrostatics
+coulombtype = Cut-off
+coulomb-modifier = Potential-shift-Verlet
+rcoulomb-switch = 0
+rcoulomb = 1
+; Relative dielectric constant for the medium and the reaction field
+epsilon-r = 1
+epsilon-rf = 0
+; Method for doing Van der Waals
+vdw-type = Cut-off
+vdw-modifier = Potential-shift-Verlet
+; cut-off lengths
+rvdw-switch = 0
+rvdw = 1
+; Apply long range dispersion corrections for Energy and Pressure
+DispCorr = No
+; Extension of the potential lookup tables beyond the cut-off
+table-extension = 1
+; Separate tables between energy group pairs
+energygrp-table =
+; Spacing for the PME/PPPM FFT grid
+fourierspacing = 0.12
+; FFT grid size, when a value is 0 fourierspacing will be used
+fourier-nx = 0
+fourier-ny = 0
+fourier-nz = 0
+; EWALD/PME/PPPM parameters
+pme-order = 4
+ewald-rtol = 1e-05
+ewald-rtol-lj = 0.001
+lj-pme-comb-rule = Geometric
+ewald-geometry = 3d
+epsilon-surface = 0
+implicit-solvent = no
+
+; OPTIONS FOR WEAK COUPLING ALGORITHMS
+; Temperature coupling
+tcoupl = No
+nsttcouple = -1
+nh-chain-length = 10
+print-nose-hoover-chain-variables = no
+; Groups to couple separately
+tc-grps =
+; Time constant (ps) and reference temperature (K)
+tau-t =
+ref-t =
+; pressure coupling
+pcoupl = No
+pcoupltype = Isotropic
+nstpcouple = -1
+; Time constant (ps), compressibility (1/bar) and reference P (bar)
+tau-p = 1
+compressibility =
+ref-p =
+; Scaling of reference coordinates, No, All or COM
+refcoord-scaling = No
+
+; OPTIONS FOR QMMM calculations
+QMMM = no
+; Groups treated Quantum Mechanically
+QMMM-grps = QMatoms
+; QM method
+QMmethod =
+; QMMM scheme
+QMMMscheme = normal
+; QM basisset
+QMbasis =
+; QM charge
+QMcharge =
+; QM multiplicity
+QMmult =
+; Surface Hopping
+SH =
+; CAS space options
+CASorbitals =
+CASelectrons =
+SAon =
+SAoff =
+SAsteps =
+; Scale factor for MM charges
+MMChargeScaleFactor = 1
+
+; SIMULATED ANNEALING
+; Type of annealing for each temperature group (no/single/periodic)
+annealing =
+; Number of time points to use for specifying annealing in each group
+annealing-npoints =
+; List of times at the annealing points for each group
+annealing-time =
+; Temp. at each annealing point, for each group.
+annealing-temp =
+
+; GENERATE VELOCITIES FOR STARTUP RUN
+gen-vel = no
+gen-temp = 300
+gen-seed = -1
+
+; OPTIONS FOR BONDS
+constraints = none
+; Type of constraint algorithm
+constraint-algorithm = Lincs
+; Do not constrain the start configuration
+continuation = no
+; Use successive overrelaxation to reduce the number of shake iterations
+Shake-SOR = no
+; Relative tolerance of shake
+shake-tol = 0.0001
+; Highest order in the expansion of the constraint coupling matrix
+lincs-order = 4
+; Number of iterations in the final step of LINCS. 1 is fine for
+; normal simulations, but use 2 to conserve energy in NVE runs.
+; For energy minimization with constraints it should be 4 to 8.
+lincs-iter = 1
+; Lincs will write a warning to the stderr if in one step a bond
+; rotates over more degrees than
+lincs-warnangle = 30
+; Convert harmonic bonds to morse potentials
+morse = no
+
+; ENERGY GROUP EXCLUSIONS
+; Pairs of energy groups for which all non-bonded interactions are excluded
+energygrp-excl =
+
+; WALLS
+; Number of walls, type, atom types, densities and box-z scale factor for Ewald
+nwall = 0
+wall-type = 9-3
+wall-r-linpot = -1
+wall-atomtype =
+wall-density =
+wall-ewald-zfac = 3
+
+; COM PULLING
+pull = no
+
+; AWH biasing
+awh = no
+
+; ENFORCED ROTATION
+; Enforced rotation: No or Yes
+rotation = no
+
+; Group to display and/or manipulate in interactive MD session
+IMD-group =
+
+; NMR refinement stuff
+; Distance restraints type: No, Simple or Ensemble
+disre = No
+; Force weighting of pairs in one distance restraint: Conservative or Equal
+disre-weighting = Conservative
+; Use sqrt of the time averaged times the instantaneous violation
+disre-mixed = no
+disre-fc = 1000
+disre-tau = 0
+; Output frequency for pair distances to energy file
+nstdisreout = 100
+; Orientation restraints: No or Yes
+orire = no
+; Orientation restraints force constant and tau for time averaging
+orire-fc = 0
+orire-tau = 0
+orire-fitgrp =
+; Output frequency for trace(SD) and S to energy file
+nstorireout = 100
+
+; Free energy variables
+free-energy = no
+couple-moltype =
+couple-lambda0 = vdw-q
+couple-lambda1 = vdw-q
+couple-intramol = no
+init-lambda = -1
+init-lambda-state = -1
+delta-lambda = 0
+nstdhdl = 50
+fep-lambdas =
+mass-lambdas =
+coul-lambdas =
+vdw-lambdas =
+bonded-lambdas =
+restraint-lambdas =
+temperature-lambdas =
+calc-lambda-neighbors = 1
+init-lambda-weights =
+dhdl-print-energy = no
+sc-alpha = 0
+sc-power = 1
+sc-r-power = 6
+sc-sigma = 0.3
+sc-coul = no
+separate-dhdl-file = yes
+dhdl-derivatives = yes
+dh_hist_size = 0
+dh_hist_spacing = 0.1
+
+; Non-equilibrium MD stuff
+acc-grps =
+accelerate =
+freezegrps =
+freezedim =
+cos-acceleration = 0
+deform =
+
+; simulated tempering variables
+simulated-tempering = no
+simulated-tempering-scaling = geometric
+sim-temp-low = 300
+sim-temp-high = 300
+
+; Ion/water position swapping for computational electrophysiology setups
+; Swap positions along direction: no, X, Y, Z
+swapcoords = no
+adress = no
+
+; User defined thingies
+user1-grps =
+user2-grps =
+userint1 = 0
+userint2 = 0
+userint3 = 0
+userint4 = 0
+userreal1 = 0
+userreal2 = 0
+userreal3 = 0
+userreal4 = 0
+; Electric fields
+; Format for electric-field-x, etc. is: four real variables:
+; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
+; and sigma (ps) width of the pulse. Omega = 0 means static field,
+; sigma = 0 means no pulse, leaving the field to be a cosine function.
+electric-field-x = 0 0 0 0
+electric-field-y = 0 0 0 0
+electric-field-z = 0 0 0 0
+</String>
+</ReferenceData>
return title;
}
-
+/*! \brief
+ * Generate exclusion lists for QM/MM.
+ *
+ * This routine updates the exclusion lists for QM atoms in order to include all other QM
+ * atoms of this molecule. Moreover, this routine replaces bonds between QM atoms with
+ * CONNBOND and, when MiMiC is not used, removes bonded interactions between QM and link atoms.
+ * Finally, in case if MiMiC QM/MM is used - charges of QM atoms are set to 0
+ *
+ * @param molt molecule type with QM atoms
+ * @param grpnr group informatio
+ * @param ir input record
+ * @param wi warning handler
+ * @param qmmmMode QM/MM mode switch: original/MiMiC
+ */
static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *grpnr,
- t_inputrec *ir, warninp_t wi)
+ t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode)
{
/* This routine expects molt->ilist to be of size F_NRE and ordered. */
}
if ((grpnr ? grpnr[i] : 0) == j)
{
- qm_arr[qm_nr++] = i;
+ qm_arr[qm_nr++] = i;
+ molt->atoms.atom[i].q = 0.0;
+ molt->atoms.atom[i].qB = 0.0;
}
}
}
numQmAtoms++;
}
}
- bexcl = (numQmAtoms >= nratoms - 1);
+
+ /* MiMiC treats link atoms as quantum atoms - therefore
+ * we do not need do additional exclusions here */
+ if (qmmmMode == GmxQmmmMode::GMX_QMMM_MIMIC)
+ {
+ bexcl = numQmAtoms == nratoms;
+ }
+ else
+ {
+ bexcl = (numQmAtoms >= nratoms - 1);
+ }
if (bexcl && ftype == F_SETTLE)
{
* linkatoms interaction with the QMatoms and would be counted
* twice. */
- for (int i = 0; i < F_NRE; i++)
+ if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
{
- if (IS_CHEMBOND(i))
+ for (int i = 0; i < F_NRE; i++)
{
- int j = 0;
- while (j < molt->ilist[i].size())
+ if (IS_CHEMBOND(i))
{
- int a1 = molt->ilist[i].iatoms[j+1];
- int a2 = molt->ilist[i].iatoms[j+2];
- if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
+ int j = 0;
+ while (j < molt->ilist[i].size())
{
- if (link_nr >= link_max)
- {
- link_max += 10;
- srenew(link_arr, link_max);
- }
- if (bQMMM[a1])
+ int a1 = molt->ilist[i].iatoms[j + 1];
+ int a2 = molt->ilist[i].iatoms[j + 2];
+ if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
{
- link_arr[link_nr++] = a2;
- }
- else
- {
- link_arr[link_nr++] = a1;
+ if (link_nr >= link_max)
+ {
+ link_max += 10;
+ srenew(link_arr, link_max);
+ }
+ if (bQMMM[a1])
+ {
+ link_arr[link_nr++] = a2;
+ }
+ else
+ {
+ link_arr[link_nr++] = a1;
+ }
}
+ j += 3;
}
- j += 3;
}
}
}
{
blink[i] = FALSE;
}
- for (int i = 0; i < link_nr; i++)
+
+ if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
{
- blink[link_arr[i]] = TRUE;
+ for (int i = 0; i < link_nr; i++)
+ {
+ blink[link_arr[i]] = TRUE;
+ }
}
/* creating the exclusion block for the QM atoms. Each QM atom has
* as excluded elements all the other QMatoms (and itself).
free(blink);
} /* generate_qmexcl */
-void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi)
+void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode)
{
/* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
*/
- unsigned char *grpnr;
- int mol, nat_mol, nr_mol_with_qm_atoms = 0;
- gmx_molblock_t *molb;
- bool bQMMM;
+ unsigned char *grpnr;
+ int mol, nat_mol, nr_mol_with_qm_atoms = 0;
+ gmx_molblock_t *molb;
+ bool bQMMM;
+ int index_offset = 0;
+ int qm_nr = 0;
grpnr = sys->groups.grpnr[egcQMMM];
bQMMM = FALSE;
for (int i = 0; i < nat_mol; i++)
{
- if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM)
+ if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM))
{
- bQMMM = TRUE;
+ bQMMM = TRUE;
+ qm_nr++;
}
}
+
if (bQMMM)
{
nr_mol_with_qm_atoms++;
sys->molblock[mb+1].nmol -= 1;
}
- /* Add a moltype for the QMMM molecule */
- sys->moltype.push_back(sys->moltype[molb->type]);
+ /* Create a copy of a moltype for a molecule
+ * containing QM atoms and append it in the end of the list
+ */
+ std::vector<gmx_moltype_t> temp(sys->moltype.size());
+ for (size_t i = 0; i < sys->moltype.size(); ++i)
+ {
+ copy_moltype(&sys->moltype[i], &temp[i]);
+ }
+ sys->moltype.resize(sys->moltype.size() + 1);
+ for (size_t i = 0; i < temp.size(); ++i)
+ {
+ copy_moltype(&temp[i], &sys->moltype[i]);
+ }
+ copy_moltype(&sys->moltype[molb->type], &sys->moltype.back());
/* Copy the exclusions to a new array, since this is the only
* thing that needs to be modified for QMMM.
*/
- copy_blocka(&sys->moltype[molb->type ].excls,
+ copy_blocka(&sys->moltype[molb->type].excls,
&sys->moltype.back().excls);
/* Set the molecule type for the QMMM molblock */
molb->type = sys->moltype.size() - 1;
}
- generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi);
+ generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi, qmmmMode);
}
if (grpnr)
{
grpnr += nat_mol;
}
+ index_offset += nat_mol;
}
}
- if (nr_mol_with_qm_atoms > 1)
+ if (qmmmMode == GmxQmmmMode::GMX_QMMM_ORIGINAL &&
+ nr_mol_with_qm_atoms > 1)
{
/* generate a warning is there are QM atoms in different
* topologies. In this case it is not possible at this stage to
struct t_gromppopts;
struct t_inputrec;
struct warninp;
+enum struct GmxQmmmMode;
typedef warninp *warninp_t;
double check_mol(const gmx_mtop_t *mtop, warninp_t wi);
warninp_t wi);
/* This routine expects sys->molt[m].ilist to be of size F_NRE and ordered. */
-void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi);
+void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode);
#endif
"no", "yes", nullptr
};
-
-
const char *ens_names[ensNR+1] =
{
"Grid", "Simple", nullptr
const char *ei_names[eiNR+1] =
{
- "md", "steep", "cg", "bd", "sd2 - removed", "nm", "l-bfgs", "tpi", "tpic", "sd", "md-vv", "md-vv-avek", nullptr
+ "md", "steep", "cg", "bd", "sd2 - removed", "nm", "l-bfgs", "tpi", "tpic", "sd", "md-vv", "md-vv-avek", "mimic", nullptr
};
const char *ecutscheme_names[ecutsNR+1] = {
* and the half step kinetic energy for temperature control
*/
enum {
- eiMD, eiSteep, eiCG, eiBD, eiSD2_REMOVED, eiNM, eiLBFGS, eiTPI, eiTPIC, eiSD1, eiVV, eiVVAK, eiNR
+ eiMD, eiSteep, eiCG, eiBD, eiSD2_REMOVED, eiNM, eiLBFGS, eiTPI, eiTPIC, eiSD1, eiVV, eiVVAK, eiMimic, eiNR
};
//! Name of the integrator algorithm
extern const char *ei_names[eiNR+1];
//! Macro returning integrator string
#define EI(e) enum_name(e, eiNR, ei_names)
+//! Do we use MiMiC QM/MM?
+#define EI_MIMIC(e) ((e) == eiMimic)
//! Do we use velocity Verlet
#define EI_VV(e) ((e) == eiVV || (e) == eiVVAK)
//! Do we use molecular dynamics
-#define EI_MD(e) ((e) == eiMD || EI_VV(e))
+#define EI_MD(e) ((e) == eiMD || EI_VV(e) || EI_MIMIC(e))
//! Do we use stochastic dynamics
#define EI_SD(e) ((e) == eiSD1)
//! Do we use any stochastic integrator
//! String corresponding to interactions in neighborlist code
extern const char *gmx_nblist_interaction_names[GMX_NBLIST_INTERACTION_NR+1];
+
+//! \brief QM/MM mode
+enum struct GmxQmmmMode {
+ GMX_QMMM_ORIGINAL,
+ GMX_QMMM_MIMIC
+};
#endif /* GMX_MDTYPES_MD_ENUMS_H */
dst->atomtypeB[i] = src->atomtypeB[i];
}
}
- dst->nres = src->nres;
+ dst->haveBState = src->haveBState;
+ dst->haveCharge = src->haveCharge;
+ dst->haveMass = src->haveMass;
+ dst->havePdbInfo = src->havePdbInfo;
+ dst->haveType = src->haveType;
+ dst->nres = src->nres;
for (i = 0; (i < src->nres); i++)
{
dst->resinfo[i] = src->resinfo[i];
}
}
}
+
+void copy_block(const t_block *src, t_block *dst)
+{
+ dst->nr = src->nr;
+ /* Workaround for inconsistent handling of nalloc_index in
+ * other parts of the code. Often nalloc_index and nalloc_a
+ * are not set.
+ */
+ dst->nalloc_index = std::max(src->nalloc_index, dst->nr + 1);
+ snew(dst->index, dst->nalloc_index);
+ for (int i = 0; i < dst->nr+1; ++i)
+ {
+ dst->index[i] = src->index[i];
+ }
+}
void copy_blocka(const t_blocka *src, t_blocka *dest);
+void copy_block(const t_block *src, t_block *dst);
+
void stupid_fill_block(t_block *grp, int natom, gmx_bool bOneIndexGroup);
/* Fill a block structure with numbers identical to the index
* (0, 1, 2, .. natom-1)
delete idef->cmap_grid;
init_idef(idef);
}
+
+void copy_ilist(const t_ilist *src, t_ilist *dst)
+{
+ dst->nr = src->nr;
+ dst->nr_nonperturbed = src->nr_nonperturbed;
+ dst->nalloc = src->nalloc;
+
+ snew(dst->iatoms, dst->nr);
+ for (int i = 0; i < dst->nr; ++i)
+ {
+ dst->iatoms[i] = src->iatoms[i];
+ }
+}
*/
void done_idef(t_idef *idef);
+void copy_ilist(const t_ilist *src, t_ilist *dst);
+
#endif
{
return (group->grpnr[type] ? group->grpnr[type][atom] : 0);
}
+
+void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst)
+{
+ dst->name = src->name;
+ copy_blocka(&src->excls, &dst->excls);
+ copy_block(&src->cgs, &dst->cgs);
+ t_atoms *atomsCopy = copy_t_atoms(&src->atoms);
+ dst->atoms = *atomsCopy;
+ sfree(atomsCopy);
+
+ for (int i = 0; i < F_NRE; ++i)
+ {
+ dst->ilist[i] = src->ilist[i];
+ }
+}
//! Deleter for gmx_localtop_t, needed until it has a proper destructor.
using ExpandedTopologyPtr = gmx::unique_cptr<gmx_localtop_t, done_and_sfree_localtop>;
+void copy_moltype(const gmx_moltype_t *src, gmx_moltype_t *dst);
+
#endif
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff