%
% This file is part of the GROMACS molecular simulation package.
%
-% Copyright (c) 2013,2014, by the GROMACS development team, led by
+% Copyright (c) 2013,2014,2015, 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.
them in {\tt atomtypes.atp} as well.
+\section{Molecule definition\index{molecule definition}}
+
+\subsection{Moleculetype entries}
+An organizational structure that usually corresponds to molecules is
+the {\tt [ moleculetype ]} entry. This entry serves two main purposes. One is
+to give structure to the topology file(s), usually corresponding to real
+molecules. This makes the topology easier to read and writing it less labor
+intensive. A second purpose is computational efficiency. The system definition
+that is kept in memory is proportional in size of the {\tt moleculetype}
+definitions. If a molecule is present in 100000 copies, this saves a factor
+of 100000 in memory, which means the system usually fits in cache, which
+can improve performance tremendously. Interactions that correspond to chemical
+bonds, that generate exclusions, can only be defined between atoms within
+a {\tt moleculetype}. It is allowed to have multiple molecules which are
+not covalently bonded in one {\tt moleculetype} definition. Molecules can
+be made infinitely long by connecting to themselves over periodic boundaries.
+When such periodic molecules are present, an option in the {\tt mdp} file
+needs to be set to tell {\gromacs} not to attempt to make molecules
+that are broken over periodic boundaries whole again.
+
+\subsection{Intermolecular interactions\index{intermolecular interaction}}
+In some cases, one would like atoms in different molecules to also interact
+with other interactions than the usual non-bonded interactions. This is often
+the case in binding studies. When the molecules are covalently bound, e.g.
+a ligand binding covalently to a protein, they are effectively one molecule
+and they should be defined in one {\tt [ moleculetype ]} entry. Note that
+{\tt pdb2gmx} has an option to put two or more molecules in one
+{\tt [ moleculetype ]} entry. When molecules are not covalently bound,
+it is much more convenient to use separate {\tt moleculetype} definitions
+and specify the intermolecular interactions in the
+{\tt [ intermolecular\_interactions] } section. In this section, which is
+placed at the end of the topology (see \tabref{topfile1}), normal bonded
+interactions can be specified using global atom indices. The only restrictions
+are that no interactions can be used that generates exclusions and no
+constraints can be used.
\subsection{Intramolecular pair interactions\index{intramolecular pair interaction}}
\label{subsec:pairinteractions}
in separate ``LJ-14'' and ``Coulomb-14'' entries per energy group pair.
Energies for {\tt [~pairs_nb~]} are added to the ``LJ-(SR)'' and ``Coulomb-(SR)'' terms.
-\subsection{Implicit solvation parameters\index{implicit solvation parameters}}
+
+\subsection{Exclusions}
+\label{sec:excl}
+The \normindex{exclusions} for non-bonded interactions are generated by {\tt
+grompp} for neighboring atoms up to a certain number of bonds away, as
+defined in the {\tt [~moleculetype~]} section in the topology file
+(see \ssecref{topfile}). Particles are considered bonded when they are
+connected by ``chemical'' bonds ({\tt [~bonds~]} types 1 to 5, 7 or 8)
+or constraints ({\tt [~constraints~]} type 1).
+Type 5 {\tt [~bonds~]} can be used to create a \normindex{connection}
+between two atoms without creating an interaction.
+There is a \normindex{harmonic interaction}
+({\tt [~bonds~]} type 6) that does not connect the atoms by a chemical bond.
+There is also a second constraint type ({\tt [~constraints~]} type 2)
+that fixes the distance, but does not connect
+the atoms by a chemical bond.
+For a complete list of all these interactions, see \tabref{topfile2}.
+
+Extra exclusions within a molecule can be added manually
+in a {\tt [~exclusions~]} section. Each line should start with one
+atom index, followed by one or more atom indices. All non-bonded
+interactions between the first atom and the other atoms will be excluded.
+
+When all non-bonded interactions within or between groups of atoms need
+to be excluded, is it more convenient and much more efficient to use
+energy monitor group exclusions (see \secref{groupconcept}).
+
+
+\section{Implicit solvation parameters\index{implicit solvation parameters}}
Starting with {\gromacs} 4.5, implicit solvent is supported. A section in the
topology has been introduced to list those parameters:
basis.
-
-\section{Exclusions}
-\label{sec:excl}
-The \normindex{exclusions} for non-bonded interactions are generated by {\tt
-grompp} for neighboring atoms up to a certain number of bonds away, as
-defined in the {\tt [~moleculetype~]} section in the topology file
-(see \ssecref{topfile}). Particles are considered bonded when they are
-connected by ``chemical'' bonds ({\tt [~bonds~]} types 1 to 5, 7 or 8)
-or constraints ({\tt [~constraints~]} type 1).
-Type 5 {\tt [~bonds~]} can be used to create a \normindex{connection}
-between two atoms without creating an interaction.
-There is a \normindex{harmonic interaction}
-({\tt [~bonds~]} type 6) that does not connect the atoms by a chemical bond.
-There is also a second constraint type ({\tt [~constraints~]} type 2)
-that fixes the distance, but does not connect
-the atoms by a chemical bond.
-For a complete list of all these interactions, see \tabref{topfile2}.
-
-Extra exclusions within a molecule can be added manually
-in a {\tt [~exclusions~]} section. Each line should start with one
-atom index, followed by one or more atom indices. All non-bonded
-interactions between the first atom and the other atoms will be excluded.
-
-When all non-bonded interactions within or between groups of atoms need
-to be excluded, is it more convenient and much more efficient to use
-energy monitor group exclusions (see \secref{groupconcept}).
-
\section{Constraint algorithms\index{constraint algorithms}}
\label{sec:constraints}
Constraints are defined in the {\tt [~constraints~]} section.
{\em mandatory} & {\tts molecules} & & & \multicolumn{2}{l|}{molecule name; number of molecules} \\
\dline
\multicolumn{6}{c}{~} \\
+\multicolumn{6}{c}{\bf \large Inter-molecular interactions} \\
+\dline
+{\em optional} & \multicolumn{4}{l}{\tts intermolecular_interactions} & \\
+\hline
+\multicolumn{6}{|c|}{one or more bonded interactions as described in \tabref{topfile2}, with two or more atoms,} \\
+\multicolumn{6}{|c|}{no interactions that generate exclusions, no constraints, use global atom numbers} \\
+\dline
+\multicolumn{6}{c}{~} \\
\multicolumn{6}{l}{`\# at' is the required number of atom type indices for this directive} \\
\multicolumn{6}{l}{`f. tp' is the value used to select this function type} \\
-\multicolumn{6}{l}{`F. E.' indicates which of the parameters for this interaction can be} \\
-\multicolumn{6}{l}{\phantom{`F. E.'} interpolated during free energy calculations} \\
+\multicolumn{6}{l}{`F. E.' indicates which of the parameters can be interpolated in free energy calculations} \\
\multicolumn{6}{l}{~$^{(cr)}$ the combination rule determines the type of LJ parameters, see~\ssecref{nbpar}}\\
\multicolumn{6}{l}{~$^{(*)}$ for {\tts dihedraltypes} one can specify 4 atoms or the inner (outer for improper) 2 atoms}\\
\multicolumn{6}{l}{~$^{(nrexcl)}$ exclude neighbors $n_{ex}$ bonds away for non-bonded interactions}\\
% The footnotetext fields only work inside the body, and not from a
% column with ``p'' formatting'!
bond & {\tts bonds}\fnm{4},\fnm{5} & 2 & 1 & $b_0$ (nm); $k_b$ (\kJmolnm{-2}) & all & \ssecref{harmonicbond}
-\label{tab:topfile2} \footnotetext[1]{The required number of atom indices for this directive}\footnotetext[2]{The index to use to select this function type}\footnotetext[3]{Indicates which of the parameters for this interaction can be interpolated during free energy calculations}\footnotetext[4]{This interaction type will be used by {{\tts grompp}} for generating exclusions}\footnotetext[5]{This interaction type can be converted to constraints by {{\tts grompp}}}\footnotetext[7]{The combination rule determines the type of LJ parameters, see~\ssecref{nbpar}}\footnotetext[6]{No connection, and so no exclusions, are generated for this interaction}
+\label{tab:topfile2} \footnotetext[1]{The required number of atom indices for this directive}\footnotetext[2]{The index to use to select this function type}\footnotetext[3]{Indicates which of the parameters can be interpolated in free energy calculations}\footnotetext[4]{This interaction type will be used by {{\tts grompp}} for generating exclusions}\footnotetext[5]{This interaction type can be converted to constraints by {{\tts grompp}}}\footnotetext[7]{The combination rule determines the type of LJ parameters, see~\ssecref{nbpar}}\footnotetext[6]{No connection, and so no exclusions, are generated for this interaction}
\\
G96 bond & {\tts bonds}\fnm{4},\fnm{5} & 2 & 2 & $b_0$ (nm); $k_b$ (\kJmolnm{-4}) & all & \ssecref{G96bond} \\
Morse & {\tts bonds}\fnm{4},\fnm{5} & 2 & 3 & $b_0$ (nm); $D$ (\kJmol); $\beta$ (nm$^{-1}$) & all & \ssecref{Morsebond} \\
comm->bCGs = (ncg_mtop(mtop) < mtop->natoms);
- comm->bInterCGBondeds = (ncg_mtop(mtop) > mtop->mols.nr);
+ comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) ||
+ mtop->bIntermolecularInteractions);
if (comm->bInterCGBondeds)
{
comm->bInterCGMultiBody = (multi_body_bondeds_count(mtop) > 0);
#include "gmxpre.h"
+#include <assert.h>
#include <stdlib.h>
#include <string.h>
/*! \brief Struct for the reverse topology: links bonded interactions to atomsx */
typedef struct gmx_reverse_top {
//! @cond Doxygen_Suppress
- gmx_bool bExclRequired; /**< Do we require all exclusions to be assigned? */
- int n_excl_at_max; /**< The maximum number of exclusions one atom can have */
- gmx_bool bConstr; /**< Are there constraints in this revserse top? */
- gmx_bool bSettle; /**< Are there settles in this revserse top? */
- gmx_bool bBCheck; /**< All bonded interactions have to be assigned? */
- gmx_bool bMultiCGmols; /**< Are the multi charge-group molecules? */
- reverse_ilist_t *ril_mt; /**< Reverse ilist for all moltypes */
- int ril_mt_tot_size; /**< The size of ril_mt[?].index summed over all entries */
- int ilsort; /**< The sorting state of bondeds for free energy */
- molblock_ind_t *mbi; /**< molblock to global atom index for quick lookup of molblocks on atom index */
- int nmolblock; /**< The number of molblocks, size of \p mbi */
+ gmx_bool bExclRequired; /**< Do we require all exclusions to be assigned? */
+ int n_excl_at_max; /**< The maximum number of exclusions one atom can have */
+ gmx_bool bConstr; /**< Are there constraints in this revserse top? */
+ gmx_bool bSettle; /**< Are there settles in this revserse top? */
+ gmx_bool bBCheck; /**< All bonded interactions have to be assigned? */
+ gmx_bool bMultiCGmols; /**< Are the multi charge-group molecules? */
+ reverse_ilist_t *ril_mt; /**< Reverse ilist for all moltypes */
+ int ril_mt_tot_size; /**< The size of ril_mt[?].index summed over all entries */
+ int ilsort; /**< The sorting state of bondeds for free energy */
+ molblock_ind_t *mbi; /**< molblock to global atom index for quick lookup of molblocks on atom index */
+ int nmolblock; /**< The number of molblocks, size of \p mbi */
+
+ gmx_bool bIntermolecularInteractions; /**< Do we have intermolecular interactions? */
+ reverse_ilist_t ril_intermol; /**< Intermolecular reverse ilist */
/* Work data structures for multi-threading */
int nthread; /**< The number of threads to be used */
}
/*! \brief Run the reverse ilist generation and store it when \p bAssign = TRUE */
-static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom,
- int **vsite_pbc,
+static int low_make_reverse_ilist(const t_ilist *il_mt,
+ const t_atom *atom,
+ int **vsite_pbc, /* should be const */
int *count,
gmx_bool bConstr, gmx_bool bSettle,
gmx_bool bBCheck,
gmx_bool bLinkToAllAtoms,
gmx_bool bAssign)
{
- int ftype, nral, i, j, nlink, link;
- t_ilist *il;
- t_iatom *ia;
- atom_id a;
- int nint;
- gmx_bool bVSite;
+ int ftype, nral, i, j, nlink, link;
+ const t_ilist *il;
+ t_iatom *ia;
+ atom_id a;
+ int nint;
+ gmx_bool bVSite;
nint = 0;
for (ftype = 0; ftype < F_NRE; ftype++)
}
/*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
-static int make_reverse_ilist(gmx_moltype_t *molt,
- int **vsite_pbc,
+static int make_reverse_ilist(const t_ilist *ilist,
+ const t_atoms *atoms,
+ int **vsite_pbc, /* should be const (C issue) */
gmx_bool bConstr, gmx_bool bSettle,
gmx_bool bBCheck,
gmx_bool bLinkToAllAtoms,
int nat_mt, *count, i, nint_mt;
/* Count the interactions */
- nat_mt = molt->atoms.nr;
+ nat_mt = atoms->nr;
snew(count, nat_mt);
- low_make_reverse_ilist(molt->ilist, molt->atoms.atom, vsite_pbc,
+ low_make_reverse_ilist(ilist, atoms->atom, vsite_pbc,
count,
bConstr, bSettle, bBCheck, NULL, NULL,
bLinkToAllAtoms, FALSE);
/* Store the interactions */
nint_mt =
- low_make_reverse_ilist(molt->ilist, molt->atoms.atom, vsite_pbc,
+ low_make_reverse_ilist(ilist, atoms->atom, vsite_pbc,
count,
bConstr, bSettle, bBCheck,
ril_mt->index, ril_mt->il,
/* Make the atom to interaction list for this molecule type */
nint_mt[mt] =
- make_reverse_ilist(molt, vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
+ make_reverse_ilist(molt->ilist, &molt->atoms,
+ vsite_pbc_molt ? vsite_pbc_molt[mt] : NULL,
rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
&rt->ril_mt[mt]);
}
sfree(nint_mt);
+ /* Make an intermolecular reverse top, if necessary */
+ rt->bIntermolecularInteractions = mtop->bIntermolecularInteractions;
+ if (rt->bIntermolecularInteractions)
+ {
+ t_atoms atoms_global;
+
+ rt->ril_intermol.index = NULL;
+ rt->ril_intermol.il = NULL;
+
+ atoms_global.nr = mtop->natoms;
+ atoms_global.atom = NULL; /* Only used with virtual sites */
+
+ *nint +=
+ make_reverse_ilist(mtop->intermolecular_ilist, &atoms_global,
+ NULL,
+ rt->bConstr, rt->bSettle, rt->bBCheck, FALSE,
+ &rt->ril_intermol);
+ }
+
if (bFE && gmx_mtop_bondeds_free_energy(mtop))
{
rt->ilsort = ilsortFE_UNSORTED;
/* If normal and/or settle constraints act only within charge groups,
* we can store them in the reverse top and simply assign them to domains.
* Otherwise we need to assign them to multiple domains and set up
- * the parallel version constraint algoirthm(s).
+ * the parallel version constraint algorithm(s).
*/
dd->reverse_top = make_reverse_top(mtop, ir->efep != efepNO,
}
/*! \brief Store a virtual site interaction, complex because of PBC and recursion */
-static void add_vsite(gmx_ga2la_t ga2la, int *index, int *rtil,
+static void add_vsite(gmx_ga2la_t ga2la, const int *index, const int *rtil,
int ftype, int nral,
gmx_bool bHomeA, int a, int a_gl, int a_mol,
- t_iatom *iatoms,
+ const t_iatom *iatoms,
t_idef *idef, int **vsite_pbc, int *vsite_pbc_nalloc)
{
int k, vsi, pbc_a_mol;
}
}
-/*! \brief This function looks up and assigns bonded interactions for zone iz.
- *
- * With thread parallelizing each thread acts on a different atom range:
- * at_start to at_end.
+/*! \brief Check and when available assign bonded interactions for local atom i
*/
-static int make_bondeds_zone(gmx_domdec_t *dd,
- const gmx_domdec_zones_t *zones,
- const gmx_molblock_t *molb,
- gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
- real rc2,
- int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
- const t_iparams *ip_in,
- t_idef *idef,
- int **vsite_pbc,
- int *vsite_pbc_nalloc,
- int iz, int nzone,
- int at_start, int at_end)
+static gmx_inline void
+check_assign_interactions_atom(int i, int i_gl,
+ int mol, int i_mol,
+ const int *index, const int *rtil,
+ gmx_bool bInterMolInteractions,
+ int ind_start, int ind_end,
+ const gmx_domdec_t *dd,
+ const gmx_domdec_zones_t *zones,
+ const gmx_molblock_t *molb,
+ gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
+ real rc2,
+ int *la2lc,
+ t_pbc *pbc_null,
+ rvec *cg_cm,
+ const t_iparams *ip_in,
+ t_idef *idef,
+ int **vsite_pbc, int *vsite_pbc_nalloc,
+ int iz,
+ gmx_bool bBCheck,
+ int *nbonded_local)
{
- int i, i_gl, mb, mt, mol, i_mol, j, ftype, nral, d, k;
- int *index, *rtil;
- t_iatom *iatoms, tiatoms[1+MAXATOMLIST];
- gmx_bool bBCheck, bUse, bLocal;
- ivec k_zero, k_plus;
- gmx_ga2la_t ga2la;
- int a_loc;
- int kz;
- int nizone;
- const gmx_domdec_ns_ranges_t *izone;
- gmx_reverse_top_t *rt;
- int nbonded_local;
-
- nizone = zones->nizone;
- izone = zones->izone;
-
- rt = dd->reverse_top;
+ int j;
- bBCheck = rt->bBCheck;
+ j = ind_start;
+ while (j < ind_end)
+ {
+ int ftype;
+ const t_iatom *iatoms;
+ int nral;
+ t_iatom tiatoms[1 + MAXATOMLIST];
- nbonded_local = 0;
+ ftype = rtil[j++];
+ iatoms = rtil + j;
+ nral = NRAL(ftype);
+ if (ftype == F_SETTLE)
+ {
+ /* Settles are only in the reverse top when they
+ * operate within a charge group. So we can assign
+ * them without checks. We do this only for performance
+ * reasons; it could be handled by the code below.
+ */
+ if (iz == 0)
+ {
+ /* Home zone: add this settle to the local topology */
+ tiatoms[0] = iatoms[0];
+ tiatoms[1] = i;
+ tiatoms[2] = i + iatoms[2] - iatoms[1];
+ tiatoms[3] = i + iatoms[3] - iatoms[1];
+ add_ifunc(nral, tiatoms, &idef->il[ftype]);
+ (*nbonded_local)++;
+ }
+ j += 1 + nral;
+ }
+ else if (interaction_function[ftype].flags & IF_VSITE)
+ {
+ assert(!bInterMolInteractions);
+ /* The vsite construction goes where the vsite itself is */
+ if (iz == 0)
+ {
+ add_vsite(dd->ga2la, index, rtil, ftype, nral,
+ TRUE, i, i_gl, i_mol,
+ iatoms, idef, vsite_pbc, vsite_pbc_nalloc);
+ }
+ j += 1 + nral + 2;
+ }
+ else
+ {
+ gmx_bool bUse;
- ga2la = dd->ga2la;
+ /* Copy the type */
+ tiatoms[0] = iatoms[0];
- for (i = at_start; i < at_end; i++)
- {
- /* Get the global atom number */
- i_gl = dd->gatindex[i];
- global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
- /* Check all interactions assigned to this atom */
- index = rt->ril_mt[mt].index;
- rtil = rt->ril_mt[mt].il;
- j = index[i_mol];
- while (j < index[i_mol+1])
- {
- ftype = rtil[j++];
- iatoms = rtil + j;
- nral = NRAL(ftype);
- if (ftype == F_SETTLE)
+ if (nral == 1)
{
- /* Settles are only in the reverse top when they
- * operate within a charge group. So we can assign
- * them without checks. We do this only for performance
- * reasons; it could be handled by the code below.
- */
+ assert(!bInterMolInteractions);
+ /* Assign single-body interactions to the home zone */
if (iz == 0)
{
- /* Home zone: add this settle to the local topology */
- tiatoms[0] = iatoms[0];
+ bUse = TRUE;
tiatoms[1] = i;
- tiatoms[2] = i + iatoms[2] - iatoms[1];
- tiatoms[3] = i + iatoms[3] - iatoms[1];
- add_ifunc(nral, tiatoms, &idef->il[ftype]);
- nbonded_local++;
+ if (ftype == F_POSRES)
+ {
+ add_posres(mol, i_mol, molb, tiatoms, ip_in,
+ idef);
+ }
+ else if (ftype == F_FBPOSRES)
+ {
+ add_fbposres(mol, i_mol, molb, tiatoms, ip_in,
+ idef);
+ }
}
- j += 1 + nral;
- }
- else if (interaction_function[ftype].flags & IF_VSITE)
- {
- /* The vsite construction goes where the vsite itself is */
- if (iz == 0)
+ else
{
- add_vsite(dd->ga2la, index, rtil, ftype, nral,
- TRUE, i, i_gl, i_mol,
- iatoms, idef, vsite_pbc, vsite_pbc_nalloc);
+ bUse = FALSE;
}
- j += 1 + nral + 2;
}
- else
+ else if (nral == 2)
{
- /* Copy the type */
- tiatoms[0] = iatoms[0];
+ /* This is a two-body interaction, we can assign
+ * analogous to the non-bonded assignments.
+ */
+ int k_gl, a_loc, kz;
- if (nral == 1)
+ if (!bInterMolInteractions)
+ {
+ /* Get the global index using the offset in the molecule */
+ k_gl = i_gl + iatoms[2] - i_mol;
+ }
+ else
+ {
+ k_gl = iatoms[2];
+ }
+ if (!ga2la_get(dd->ga2la, k_gl, &a_loc, &kz))
{
- /* Assign single-body interactions to the home zone */
- if (iz == 0)
+ bUse = FALSE;
+ }
+ else
+ {
+ if (kz >= zones->n)
+ {
+ kz -= zones->n;
+ }
+ /* Check zone interaction assignments */
+ bUse = ((iz < zones->nizone &&
+ iz <= kz &&
+ kz >= zones->izone[iz].j0 &&
+ kz < zones->izone[iz].j1) ||
+ (kz < zones->nizone &&
+ iz > kz &&
+ iz >= zones->izone[kz].j0 &&
+ iz < zones->izone[kz].j1));
+ if (bUse)
{
- bUse = TRUE;
tiatoms[1] = i;
- if (ftype == F_POSRES)
+ tiatoms[2] = a_loc;
+ /* If necessary check the cgcm distance */
+ if (bRCheck2B &&
+ dd_dist2(pbc_null, cg_cm, la2lc,
+ tiatoms[1], tiatoms[2]) >= rc2)
{
- add_posres(mol, i_mol, &molb[mb], tiatoms, ip_in,
- idef);
- }
- else if (ftype == F_FBPOSRES)
- {
- add_fbposres(mol, i_mol, &molb[mb], tiatoms, ip_in,
- idef);
+ bUse = FALSE;
}
}
+ }
+ }
+ else
+ {
+ /* Assign this multi-body bonded interaction to
+ * the local node if we have all the atoms involved
+ * (local or communicated) and the minimum zone shift
+ * in each dimension is zero, for dimensions
+ * with 2 DD cells an extra check may be necessary.
+ */
+ ivec k_zero, k_plus;
+ int k;
+
+ bUse = TRUE;
+ clear_ivec(k_zero);
+ clear_ivec(k_plus);
+ for (k = 1; k <= nral && bUse; k++)
+ {
+ gmx_bool bLocal;
+ int k_gl, a_loc;
+ int kz;
+
+ if (!bInterMolInteractions)
+ {
+ /* Get the global index using the offset in the molecule */
+ k_gl = i_gl + iatoms[k] - i_mol;
+ }
else
{
- bUse = FALSE;
+ k_gl = iatoms[k];
}
- }
- else if (nral == 2)
- {
- /* This is a two-body interaction, we can assign
- * analogous to the non-bonded assignments.
- */
- if (!ga2la_get(ga2la, i_gl+iatoms[2]-i_mol, &a_loc, &kz))
+ bLocal = ga2la_get(dd->ga2la, k_gl, &a_loc, &kz);
+ if (!bLocal || kz >= zones->n)
{
+ /* We do not have this atom of this interaction
+ * locally, or it comes from more than one cell
+ * away.
+ */
bUse = FALSE;
}
else
{
- if (kz >= nzone)
- {
- kz -= nzone;
- }
- /* Check zone interaction assignments */
- bUse = ((iz < nizone && iz <= kz &&
- izone[iz].j0 <= kz && kz < izone[iz].j1) ||
- (kz < nizone &&iz > kz &&
- izone[kz].j0 <= iz && iz < izone[kz].j1));
- if (bUse)
+ int d;
+
+ tiatoms[k] = a_loc;
+ for (d = 0; d < DIM; d++)
{
- tiatoms[1] = i;
- tiatoms[2] = a_loc;
- /* If necessary check the cgcm distance */
- if (bRCheck2B &&
- dd_dist2(pbc_null, cg_cm, la2lc,
- tiatoms[1], tiatoms[2]) >= rc2)
+ if (zones->shift[kz][d] == 0)
{
- bUse = FALSE;
+ k_zero[d] = k;
+ }
+ else
+ {
+ k_plus[d] = k;
}
}
}
}
- else
+ bUse = (bUse &&
+ k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
+ if (bRCheckMB)
{
- /* Assign this multi-body bonded interaction to
- * the local node if we have all the atoms involved
- * (local or communicated) and the minimum zone shift
- * in each dimension is zero, for dimensions
- * with 2 DD cells an extra check may be necessary.
- */
- bUse = TRUE;
- clear_ivec(k_zero);
- clear_ivec(k_plus);
- for (k = 1; k <= nral && bUse; k++)
+ int d;
+
+ for (d = 0; (d < DIM && bUse); d++)
{
- bLocal = ga2la_get(ga2la, i_gl+iatoms[k]-i_mol,
- &a_loc, &kz);
- if (!bLocal || kz >= zones->n)
+ /* Check if the cg_cm distance falls within
+ * the cut-off to avoid possible multiple
+ * assignments of bonded interactions.
+ */
+ if (rcheck[d] &&
+ k_plus[d] &&
+ dd_dist2(pbc_null, cg_cm, la2lc,
+ tiatoms[k_zero[d]], tiatoms[k_plus[d]]) >= rc2)
{
- /* We do not have this atom of this interaction
- * locally, or it comes from more than one cell
- * away.
- */
bUse = FALSE;
}
- else
- {
- tiatoms[k] = a_loc;
- for (d = 0; d < DIM; d++)
- {
- if (zones->shift[kz][d] == 0)
- {
- k_zero[d] = k;
- }
- else
- {
- k_plus[d] = k;
- }
- }
- }
- }
- bUse = (bUse &&
- k_zero[XX] && k_zero[YY] && k_zero[ZZ]);
- if (bRCheckMB)
- {
- for (d = 0; (d < DIM && bUse); d++)
- {
- /* Check if the cg_cm distance falls within
- * the cut-off to avoid possible multiple
- * assignments of bonded interactions.
- */
- if (rcheck[d] &&
- k_plus[d] &&
- dd_dist2(pbc_null, cg_cm, la2lc,
- tiatoms[k_zero[d]], tiatoms[k_plus[d]]) >= rc2)
- {
- bUse = FALSE;
- }
- }
}
}
- if (bUse)
+ }
+ if (bUse)
+ {
+ /* Add this interaction to the local topology */
+ add_ifunc(nral, tiatoms, &idef->il[ftype]);
+ /* Sum so we can check in global_stat
+ * if we have everything.
+ */
+ if (bBCheck ||
+ !(interaction_function[ftype].flags & IF_LIMZERO))
{
- /* Add this interaction to the local topology */
- add_ifunc(nral, tiatoms, &idef->il[ftype]);
- /* Sum so we can check in global_stat
- * if we have everything.
- */
- if (bBCheck ||
- !(interaction_function[ftype].flags & IF_LIMZERO))
- {
- nbonded_local++;
- }
+ (*nbonded_local)++;
}
- j += 1 + nral;
}
+ j += 1 + nral;
+ }
+ }
+}
+
+/*! \brief This function looks up and assigns bonded interactions for zone iz.
+ *
+ * With thread parallelizing each thread acts on a different atom range:
+ * at_start to at_end.
+ */
+static int make_bondeds_zone(gmx_domdec_t *dd,
+ const gmx_domdec_zones_t *zones,
+ const gmx_molblock_t *molb,
+ gmx_bool bRCheckMB, ivec rcheck, gmx_bool bRCheck2B,
+ real rc2,
+ int *la2lc, t_pbc *pbc_null, rvec *cg_cm,
+ const t_iparams *ip_in,
+ t_idef *idef,
+ int **vsite_pbc,
+ int *vsite_pbc_nalloc,
+ int izone,
+ int at_start, int at_end)
+{
+ int i, i_gl, mb, mt, mol, i_mol;
+ int *index, *rtil;
+ gmx_bool bBCheck;
+ gmx_reverse_top_t *rt;
+ int nbonded_local;
+
+ rt = dd->reverse_top;
+
+ bBCheck = rt->bBCheck;
+
+ nbonded_local = 0;
+
+ for (i = at_start; i < at_end; i++)
+ {
+ /* Get the global atom number */
+ i_gl = dd->gatindex[i];
+ global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
+ /* Check all intramolecular interactions assigned to this atom */
+ index = rt->ril_mt[mt].index;
+ rtil = rt->ril_mt[mt].il;
+
+ check_assign_interactions_atom(i, i_gl, mol, i_mol,
+ index, rtil, FALSE,
+ index[i_mol], index[i_mol+1],
+ dd, zones,
+ &molb[mb],
+ bRCheckMB, rcheck, bRCheck2B, rc2,
+ la2lc,
+ pbc_null,
+ cg_cm,
+ ip_in,
+ idef, vsite_pbc, vsite_pbc_nalloc,
+ izone,
+ bBCheck,
+ &nbonded_local);
+
+
+ if (rt->bIntermolecularInteractions)
+ {
+ /* Check all intermolecular interactions assigned to this atom */
+ index = rt->ril_intermol.index;
+ rtil = rt->ril_intermol.il;
+
+ check_assign_interactions_atom(i, i_gl, mol, i_mol,
+ index, rtil, TRUE,
+ index[i_gl], index[i_gl + 1],
+ dd, zones,
+ &molb[mb],
+ bRCheckMB, rcheck, bRCheck2B, rc2,
+ la2lc,
+ pbc_null,
+ cg_cm,
+ ip_in,
+ idef, vsite_pbc, vsite_pbc_nalloc,
+ izone,
+ bBCheck,
+ &nbonded_local);
}
}
t_blocka *lexcls, int *excl_count)
{
int nzone_bondeds, nzone_excl;
- int iz, cg0, cg1;
+ int izone, cg0, cg1;
real rc2;
int nbonded_local;
int thread;
lexcls->nra = 0;
*excl_count = 0;
- for (iz = 0; iz < nzone_bondeds; iz++)
+ for (izone = 0; izone < nzone_bondeds; izone++)
{
- cg0 = zones->cg_range[iz];
- cg1 = zones->cg_range[iz+1];
+ cg0 = zones->cg_range[izone];
+ cg1 = zones->cg_range[izone + 1];
#pragma omp parallel for num_threads(rt->nthread) schedule(static)
for (thread = 0; thread < rt->nthread; thread++)
la2lc, pbc_null, cg_cm, idef->iparams,
idef_t,
vsite_pbc, vsite_pbc_nalloc,
- iz, zones->n,
+ izone,
dd->cgindex[cg0t], dd->cgindex[cg1t]);
- if (iz < nzone_excl)
+ if (izone < nzone_excl)
{
if (thread == 0)
{
make_exclusions_zone(dd, zones,
mtop->moltype, cginfo,
excl_t,
- iz,
+ izone,
cg0t, cg1t);
}
else
mtop->moltype, bRCheck2B, rc2,
la2lc, pbc_null, cg_cm, cginfo,
excl_t,
- iz,
+ izone,
cg0t, cg1t);
}
}
nbonded_local += rt->th_work[thread].nbonded;
}
- if (iz < nzone_excl)
+ if (izone < nzone_excl)
{
if (rt->nthread > 1)
{
/* Some zones might not have exclusions, but some code still needs to
* loop over the index, so we set the indices here.
*/
- for (iz = nzone_excl; iz < zones->nizone; iz++)
+ for (izone = nzone_excl; izone < zones->nizone; izone++)
{
- set_no_exclusions_zone(dd, zones, iz, lexcls);
+ set_no_exclusions_zone(dd, zones, izone, lexcls);
}
finish_local_exclusions(dd, zones, lexcls);
t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd,
cginfo_mb_t *cginfo_mb)
{
- gmx_reverse_top_t *rt;
+ gmx_bool bExclRequired;
int mb, cg_offset, cg, cg_gl, a, aj, i, j, ftype, nral, nlink_mol, mol, ncgi;
gmx_molblock_t *molb;
gmx_moltype_t *molt;
t_block *cgs;
t_blocka *excls;
int *a2c;
- reverse_ilist_t ril;
+ reverse_ilist_t ril, ril_intermol;
t_blocka *link;
cginfo_mb_t *cgi_mb;
* which are also stored in reverse_top.
*/
- rt = dd->reverse_top;
+ bExclRequired = dd->reverse_top->bExclRequired;
+
+ if (mtop->bIntermolecularInteractions)
+ {
+ t_atoms atoms;
+
+ atoms.nr = mtop->natoms;
+ atoms.atom = NULL;
+
+ make_reverse_ilist(mtop->intermolecular_ilist, &atoms,
+ NULL, FALSE, FALSE, FALSE, TRUE, &ril_intermol);
+ }
snew(link, 1);
snew(link->index, ncg_mtop(mtop)+1);
* to all atoms, not only the first atom as in gmx_reverse_top.
* The constraints are discarded here.
*/
- make_reverse_ilist(molt, NULL, FALSE, FALSE, FALSE, TRUE, &ril);
+ make_reverse_ilist(molt->ilist, &molt->atoms,
+ NULL, FALSE, FALSE, FALSE, TRUE, &ril);
cgi_mb = &cginfo_mb[mb];
- for (cg = 0; cg < cgs->nr; cg++)
+ for (mol = 0; mol < (mtop->bIntermolecularInteractions ? molb->nmol : 1); mol++)
{
- cg_gl = cg_offset + cg;
- link->index[cg_gl+1] = link->index[cg_gl];
- for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
+ for (cg = 0; cg < cgs->nr; cg++)
{
- i = ril.index[a];
- while (i < ril.index[a+1])
+ cg_gl = cg_offset + cg;
+ link->index[cg_gl+1] = link->index[cg_gl];
+ for (a = cgs->index[cg]; a < cgs->index[cg+1]; a++)
{
- ftype = ril.il[i++];
- nral = NRAL(ftype);
- /* Skip the ifunc index */
- i++;
- for (j = 0; j < nral; j++)
+ i = ril.index[a];
+ while (i < ril.index[a+1])
{
- aj = ril.il[i+j];
- if (a2c[aj] != cg)
+ ftype = ril.il[i++];
+ nral = NRAL(ftype);
+ /* Skip the ifunc index */
+ i++;
+ for (j = 0; j < nral; j++)
{
- check_link(link, cg_gl, cg_offset+a2c[aj]);
+ aj = ril.il[i+j];
+ if (a2c[aj] != cg)
+ {
+ check_link(link, cg_gl, cg_offset+a2c[aj]);
+ }
}
+ i += nral_rt(ftype);
}
- i += nral_rt(ftype);
- }
- if (rt->bExclRequired)
- {
- /* Exclusions always go both ways */
- for (j = excls->index[a]; j < excls->index[a+1]; j++)
+ if (bExclRequired)
+ {
+ /* Exclusions always go both ways */
+ for (j = excls->index[a]; j < excls->index[a+1]; j++)
+ {
+ aj = excls->a[j];
+ if (a2c[aj] != cg)
+ {
+ check_link(link, cg_gl, cg_offset+a2c[aj]);
+ }
+ }
+ }
+
+ if (mtop->bIntermolecularInteractions)
{
- aj = excls->a[j];
- if (a2c[aj] != cg)
+ i = ril_intermol.index[a];
+ while (i < ril.index[a+1])
{
- check_link(link, cg_gl, cg_offset+a2c[aj]);
+ ftype = ril_intermol.il[i++];
+ nral = NRAL(ftype);
+ /* Skip the ifunc index */
+ i++;
+ for (j = 0; j < nral; j++)
+ {
+ aj = ril_intermol.il[i+j];
+ if (a2c[aj] != cg_offset + cg)
+ {
+ check_link(link, cg_gl, a2c[aj]);
+ }
+ }
+ i += nral_rt(ftype);
}
}
}
+ if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
+ {
+ SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
+ ncgi++;
+ }
}
- if (link->index[cg_gl+1] - link->index[cg_gl] > 0)
- {
- SET_CGINFO_BOND_INTER(cgi_mb->cginfo[cg]);
- ncgi++;
- }
- }
- nlink_mol = link->index[cg_offset+cgs->nr] - link->index[cg_offset];
- cg_offset += cgs->nr;
+ cg_offset += cgs->nr;
+ }
+ nlink_mol = link->index[cg_offset] - link->index[cg_offset-cgs->nr];
destroy_reverse_ilist(&ril);
sfree(a2c);
fprintf(debug, "molecule type '%s' %d cgs has %d cg links through bonded interac.\n", *molt->name, cgs->nr, nlink_mol);
}
- if (molb->nmol > 1)
+ if (molb->nmol > mol)
{
/* Copy the data for the rest of the molecules in this block */
- link->nalloc_a += (molb->nmol - 1)*nlink_mol;
+ link->nalloc_a += (molb->nmol - mol)*nlink_mol;
srenew(link->a, link->nalloc_a);
- for (mol = 1; mol < molb->nmol; mol++)
+ for (; mol < molb->nmol; mol++)
{
for (cg = 0; cg < cgs->nr; cg++)
{
}
}
+ if (mtop->bIntermolecularInteractions)
+ {
+ destroy_reverse_ilist(&ril_intermol);
+ }
+
if (debug)
{
fprintf(debug, "Of the %d charge groups %d are linked via bonded interactions\n", ncg_mtop(mtop), ncgi);
tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
tpxv_RemoveObsoleteParameters1, /**< remove optimize_fft, dihre_fc, nstcheckpoint */
tpxv_PullCoordTypeGeom, /**< add pull type and geometry per group and flat-bottom */
- tpxv_PullGeomDirRel /**< add pull geometry direction-relative */
+ tpxv_PullGeomDirRel, /**< add pull geometry direction-relative */
+ tpxv_IntermolecularBondeds /**< permit inter-molecular bonded interactions in the topology */
};
/*! \brief Version number of the file format written to run input
*
* When developing a feature branch that needs to change the run input
* file format, change tpx_tag instead. */
-static const int tpx_version = tpxv_PullGeomDirRel;
+static const int tpx_version = tpxv_IntermolecularBondeds;
/* This number should only be increased when you edit the TOPOLOGY section
static void set_disres_npair(gmx_mtop_t *mtop)
{
- int mt, i, npair;
- t_iparams *ip;
- t_ilist *il;
- t_iatom *a;
+ t_iparams *ip;
+ gmx_mtop_ilistloop_t iloop;
+ t_ilist *ilist, *il;
+ int nmol, i, npair;
+ t_iatom *a;
ip = mtop->ffparams.iparams;
- for (mt = 0; mt < mtop->nmoltype; mt++)
+ iloop = gmx_mtop_ilistloop_init(mtop);
+ while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
{
- il = &mtop->moltype[mt].ilist[F_DISRES];
+ il = &ilist[F_DISRES];
+
if (il->nr > 0)
{
a = il->iatoms;
mtop->molblock[0].nposres_xB = 0;
}
+ if (file_version >= tpxv_IntermolecularBondeds)
+ {
+ gmx_fio_do_gmx_bool(fio, mtop->bIntermolecularInteractions);
+ if (mtop->bIntermolecularInteractions)
+ {
+ if (bRead)
+ {
+ snew(mtop->intermolecular_ilist, F_NRE);
+ }
+ do_ilists(fio, mtop->intermolecular_ilist, bRead, file_version);
+ }
+ }
+ else
+ {
+ mtop->bIntermolecularInteractions = FALSE;
+ }
+
do_atomtypes (fio, &(mtop->atomtypes), bRead, file_version);
if (bRead && debug)
{
void pr_mtop(FILE *fp, int indent, const char *title, gmx_mtop_t *mtop,
gmx_bool bShowNumbers)
{
- int mt, mb;
+ int mt, mb, j;
if (available(fp, mtop, indent, title))
{
{
pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
}
+ pr_str(fp, indent, "bIntermolecularInteractions", EBOOL(mtop->bIntermolecularInteractions));
+ if (mtop->bIntermolecularInteractions)
+ {
+ for (j = 0; (j < F_NRE); j++)
+ {
+ pr_ilist(fp, indent, interaction_function[j].longname,
+ mtop->ffparams.functype,
+ &mtop->intermolecular_ilist[j], bShowNumbers);
+ }
+ }
pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
for (mt = 0; mt < mtop->nmoltype; mt++)
pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
+ pr_str(fp, indent, "bIntermolecularInteractions", EBOOL(top->bIntermolecularInteractions));
pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
pr_idef(fp, indent, "idef", &top->idef, bShowNumbers);
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, 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.
}
void convert_params(int atnr, t_params nbtypes[],
- t_molinfo *mi, int comb, double reppow, real fudgeQQ,
+ t_molinfo *mi, t_molinfo *intermolecular_interactions,
+ int comb, double reppow, real fudgeQQ,
gmx_mtop_t *mtop)
{
int i, j, maxtypes, mt;
}
}
}
+
+ mtop->bIntermolecularInteractions = FALSE;
+ if (intermolecular_interactions != NULL)
+ {
+ /* Process the intermolecular interaction list */
+ snew(mtop->intermolecular_ilist, F_NRE);
+
+ for (i = 0; (i < F_NRE); i++)
+ {
+ mtop->intermolecular_ilist[i].nr = 0;
+ mtop->intermolecular_ilist[i].iatoms = NULL;
+
+ plist = intermolecular_interactions->plist;
+
+ if (plist[i].nr > 0)
+ {
+ flags = interaction_function[i].flags;
+ /* For intermolecular interactions we (currently)
+ * only support potentials.
+ * Constraints and virtual sites would be possible,
+ * but require a lot of extra (bug-prone) code.
+ */
+ if (!(flags & IF_BOND))
+ {
+ gmx_fatal(FARGS, "The intermolecular_interaction section may only contain bonded potentials");
+ }
+ else if (NRAL(i) == 1) /* e.g. position restraints */
+ {
+ gmx_fatal(FARGS, "Single atom interactions don't make sense in the intermolecular_interaction section, you can put them in the moleculetype section");
+ }
+ else if (flags & IF_CHEMBOND)
+ {
+ gmx_fatal(FARGS, "The intermolecular_interaction can not contain chemically bonding interactions");
+ }
+ else
+ {
+ enter_function(&(plist[i]), (t_functype)i, comb, reppow,
+ ffp, &mtop->intermolecular_ilist[i],
+ &maxtypes, FALSE, FALSE);
+
+ mtop->bIntermolecularInteractions = TRUE;
+ }
+ }
+ }
+
+ if (!mtop->bIntermolecularInteractions)
+ {
+ sfree(mtop->intermolecular_ilist);
+ }
+ }
+
if (debug)
{
fprintf(debug, "%s, line %d: There are %d functypes in idef\n",
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, 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.
#endif
void convert_params(int atnr, t_params nbtypes[],
- t_molinfo *mi, int comb, double reppow, real fudgeQQ,
+ t_molinfo *mi,
+ t_molinfo *intermolecular_interactions,
+ int comb, double reppow, real fudgeQQ,
gmx_mtop_t *mtop);
#ifdef __cplusplus
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, 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.
d_orientation_restraints,
d_dihedral_restraints,
d_cmap,
+ d_intermolecular_interactions,
d_maxdir,
d_invalid,
d_none
t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero,
gmx_bool bGenVel, gmx_bool bVerbose, t_state *state,
gpp_atomtype_t atype, gmx_mtop_t *sys,
- int *nmi, t_molinfo **mi, t_params plist[],
+ int *nmi, t_molinfo **mi, t_molinfo **intermolecular_interactions,
+ t_params plist[],
int *comb, double *reppow, real *fudgeQQ,
gmx_bool bMorse,
warninp_t wi)
/* TOPOLOGY processing */
sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab),
plist, comb, reppow, fudgeQQ,
- atype, &nrmols, &molinfo, ir,
+ atype, &nrmols, &molinfo, intermolecular_interactions,
+ ir,
&nmolblock, &molblock, bGB,
wi);
t_gromppopts *opts;
gmx_mtop_t *sys;
int nmi;
- t_molinfo *mi;
+ t_molinfo *mi, *intermolecular_interactions;
gpp_atomtype_t atype;
t_inputrec *ir;
int natoms, nvsite, comb, mt;
snew(state, 1);
new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
opts, ir, bZero, bGenVel, bVerbose, state,
- atype, sys, &nmi, &mi, plist, &comb, &reppow, &fudgeQQ,
+ atype, sys, &nmi, &mi, &intermolecular_interactions,
+ plist, &comb, &reppow, &fudgeQQ,
opts->bMorse,
wi);
}
ntype = get_atomtype_ntypes(atype);
- convert_params(ntype, plist, mi, comb, reppow, fudgeQQ, sys);
+ convert_params(ntype, plist, mi, intermolecular_interactions,
+ comb, reppow, fudgeQQ, sys);
if (debug)
{
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015, 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.
"orientation_restraints",
"dihedral_restraints",
"cmap",
+ "intermolecular_interactions",
"invalid"
};
set_nec(&(necessary[d_orientation_restraints]), d_atoms, d_none);
set_nec(&(necessary[d_dihedral_restraints]), d_atoms, d_none);
set_nec(&(necessary[d_cmap]), d_atoms, d_none);
+ set_nec(&(necessary[d_intermolecular_interactions]), d_molecules, d_none);
for (i = 0; (i < d_maxdir); i++)
{
}
+static void make_atoms_sys(int nmolb, const gmx_molblock_t *molb,
+ const t_molinfo *molinfo,
+ t_atoms *atoms)
+{
+ int mb, m, a;
+ const t_atoms *mol_atoms;
+
+ atoms->nr = 0;
+ atoms->atom = NULL;
+
+ for (mb = 0; mb < nmolb; mb++)
+ {
+ mol_atoms = &molinfo[molb[mb].type].atoms;
+
+ srenew(atoms->atom, atoms->nr + molb[mb].nmol*mol_atoms->nr);
+
+ for (m = 0; m < molb[mb].nmol; m++)
+ {
+ for (a = 0; a < mol_atoms->nr; a++)
+ {
+ atoms->atom[atoms->nr++] = mol_atoms->atom[a];
+ }
+ }
+ }
+}
static char **read_topol(const char *infile, const char *outfile,
gpp_atomtype_t atype,
int *nrmols,
t_molinfo **molinfo,
+ t_molinfo **intermolecular_interactions,
t_params plist[],
int *combination_rule,
double *reppow,
DirStack *DS;
directive d, newd;
t_nbparam **nbparam, **pair;
+ gmx_bool bIntermolecularInteractions;
t_block2 *block2;
real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
gmx_bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
pair = NULL; /* The temporary pair interaction matrix */
block2 = NULL; /* the extra exclusions */
nb_funct = F_LJ;
+
*reppow = 12.0; /* Default value for repulsion power */
- comb = 0;
+ *intermolecular_interactions = NULL;
/* Init the number of CMAP torsion angles and grid spacing */
plist[F_CMAP].grid_spacing = 0;
cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
/* d = d_invalid; */
}
+
+ if (d == d_intermolecular_interactions)
+ {
+ if (*intermolecular_interactions == NULL)
+ {
+ /* We (mis)use the moleculetype processing
+ * to process the intermolecular interactions
+ * by making a "molecule" of the size of the system.
+ */
+ snew(*intermolecular_interactions, 1);
+ init_molinfo(*intermolecular_interactions);
+ mi0 = *intermolecular_interactions;
+ make_atoms_sys(nmolb, molb, *molinfo,
+ &mi0->atoms);
+ }
+ }
}
sfree(dirstr);
}
fudgeLJ = 1.0;
*fudgeQQ = 1.0;
- get_nbparm(nb_str, comb_str, &nb_funct, &comb, wi);
- *combination_rule = comb;
+ get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
if (nscan >= 3)
{
bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
}
ntype = get_atomtype_ntypes(atype);
ncombs = (ntype*(ntype+1))/2;
- generate_nbparams(comb, nb_funct, &(plist[nb_funct]), atype, wi);
+ generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi);
ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
ntype);
fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
free_nbparam(nbparam, ntype);
if (bGenPairs)
{
- gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, comb);
+ gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule);
ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
ntype);
fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
done_bond_atomtype(&batype);
+ if (*intermolecular_interactions != NULL)
+ {
+ sfree(mi0->atoms.atom);
+ }
+
*nrmols = nmol;
*nmolblock = nmolb;
gpp_atomtype_t atype,
int *nrmols,
t_molinfo **molinfo,
+ t_molinfo **intermolecular_interactions,
t_inputrec *ir,
int *nmolblock,
gmx_molblock_t **molblock,
printf("processing topology...\n");
}
title = read_topol(topfile, tmpfile, opts->define, opts->include,
- symtab, atype, nrmols, molinfo,
+ symtab, atype,
+ nrmols, molinfo, intermolecular_interactions,
plist, combination_rule, repulsion_power,
opts, fudgeQQ, nmolblock, molblock,
ir->efep != efepNO, bGenborn, bZero, wi);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015, 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.
gpp_atomtype_t atype,
int *nrmols,
t_molinfo **molinfo,
+ t_molinfo **intermolecular_interactions,
t_inputrec *ir,
int *nmolblock,
gmx_molblock_t **molblock,
gmx_bool bGB,
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);
bc_moltype(cr, &mtop->symtab, &mtop->moltype[i]);
}
+ block_bc(cr, mtop->bIntermolecularInteractions);
+ if (mtop->bIntermolecularInteractions)
+ {
+ snew_bc(cr, mtop->intermolecular_ilist, F_NRE);
+ bc_ilists(cr, mtop->intermolecular_ilist);
+ }
+
block_bc(cr, mtop->nmolblock);
snew_bc(cr, mtop->molblock, mtop->nmolblock);
for (i = 0; i < mtop->nmolblock; i++)
{
if (!DOMAINDECOMP(cr))
{
+ gmx_bool bSHAKE;
+
+ bSHAKE = (ir->eConstrAlg == econtSHAKE &&
+ (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 ||
+ gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0));
+
/* The group cut-off scheme and SHAKE assume charge groups
* are whole, but not using molpbc is faster in most cases.
+ * With intermolecular interactions we need PBC for calculating
+ * distances between atoms in different molecules.
*/
- if (fr->cutoff_scheme == ecutsGROUP ||
- (ir->eConstrAlg == econtSHAKE &&
- (gmx_mtop_ftype_count(mtop, F_CONSTR) > 0 ||
- gmx_mtop_ftype_count(mtop, F_CONSTRNC) > 0)))
+ if ((fr->cutoff_scheme == ecutsGROUP || bSHAKE) &&
+ !mtop->bIntermolecularInteractions)
{
fr->bMolPBC = ir->bPeriodicMols;
+
+ if (bSHAKE && fr->bMolPBC)
+ {
+ gmx_fatal(FARGS, "SHAKE is not supported with periodic molecules");
+ }
}
else
{
fr->bMolPBC = TRUE;
+
if (getenv("GMX_USE_GRAPH") != NULL)
{
fr->bMolPBC = FALSE;
if (fp)
{
- fprintf(fp, "\nGMX_MOLPBC is set, using the graph for bonded interactions\n\n");
+ md_print_warn(cr, fp, "GMX_USE_GRAPH is set, using the graph for bonded interactions\n");
}
+
+ if (mtop->bIntermolecularInteractions)
+ {
+ md_print_warn(cr, fp, "WARNING: Molecules linked by intermolecular interactions have to reside in the same periodic image, otherwise artifacts will occur!\n");
+ }
+ }
+
+ if (bSHAKE && fr->bMolPBC)
+ {
+ gmx_fatal(FARGS, "SHAKE is not properly supported with intermolecular interactions. For short simulations where linked molecules remain in the same periodic image, the environment variable GMX_USE_GRAPH can be used to override this check.\n");
}
}
}
cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
cmp_block(fp, &t1->mols, &t2->mols, "mols");
+ cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
}
else
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2008,2009,2010,2012,2013,2014, by the GROMACS development team, led by
+ * Copyright (c) 2008,2009,2010,2012,2013,2014,2015, 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.
}
iloop->mblock++;
- if (iloop->mblock == iloop->mtop->nmolblock)
+ if (iloop->mblock >= iloop->mtop->nmolblock)
{
+ if (iloop->mblock == iloop->mtop->nmolblock &&
+ iloop->mtop->bIntermolecularInteractions)
+ {
+ *ilist_mol = iloop->mtop->intermolecular_ilist;
+ *nmol = 1;
+ return TRUE;
+ }
+
gmx_mtop_ilistloop_destroy(iloop);
return FALSE;
}
iloop->mol++;
- if (iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
+ /* Inter-molecular interactions, if present, are indexed with
+ * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
+ * check for this value in this conditional.
+ */
+ if (iloop->mblock == iloop->mtop->nmolblock ||
+ iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
{
iloop->mblock++;
iloop->mol = 0;
- if (iloop->mblock == iloop->mtop->nmolblock)
+ if (iloop->mblock >= iloop->mtop->nmolblock)
{
+ if (iloop->mblock == iloop->mtop->nmolblock &&
+ iloop->mtop->bIntermolecularInteractions)
+ {
+ *ilist_mol = iloop->mtop->intermolecular_ilist;
+ *atnr_offset = 0;
+ return TRUE;
+ }
+
gmx_mtop_ilistloop_all_destroy(iloop);
return FALSE;
}
n += nmol*il[ftype].nr/(1+NRAL(ftype));
}
+ if (mtop->bIntermolecularInteractions)
+ {
+ n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
+ }
+
return n;
}
natoms += molb->nmol*srcnr;
}
+ if (mtop->bIntermolecularInteractions)
+ {
+ for (ftype = 0; ftype < F_NRE; ftype++)
+ {
+ ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
+ 1, 0, mtop->natoms);
+ }
+ }
+
if (ir == NULL)
{
top->idef.ilsort = ilsortUNKNOWN;
gen_local_top(mtop, NULL, FALSE, <op);
- top.name = mtop->name;
- top.idef = ltop.idef;
- top.atomtypes = ltop.atomtypes;
- top.cgs = ltop.cgs;
- top.excls = ltop.excls;
- top.atoms = gmx_mtop_global_atoms(mtop);
- top.mols = mtop->mols;
- top.symtab = mtop->symtab;
+ top.name = mtop->name;
+ top.idef = ltop.idef;
+ top.atomtypes = ltop.atomtypes;
+ top.cgs = ltop.cgs;
+ top.excls = ltop.excls;
+ top.atoms = gmx_mtop_global_atoms(mtop);
+ top.mols = mtop->mols;
+ top.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
+ top.symtab = mtop->symtab;
/* We only need to free the moltype and molblock data,
* all other pointers have been copied to top.
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011,2014, by the GROMACS development team, led by
+ * Copyright (c) 2011,2014,2015, 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.
typedef struct gmx_moltype_t
{
- char **name; /* Name of the molecule type */
- t_atoms atoms; /* The atoms */
- t_ilist ilist[F_NRE];
- t_block cgs; /* The charge groups */
- t_blocka excls; /* The exclusions */
+ char **name; /* Name of the molecule type */
+ t_atoms atoms; /* The atoms in this molecule */
+ t_ilist ilist[F_NRE]; /* Interaction list with local indices */
+ t_block cgs; /* The charge groups */
+ t_blocka excls; /* The exclusions */
} gmx_moltype_t;
typedef struct gmx_molblock_t
gmx_moltype_t *moltype;
int nmolblock;
gmx_molblock_t *molblock;
+ gmx_bool bIntermolecularInteractions; /* Are there intermolecular
+ * interactions? */
+ t_ilist *intermolecular_ilist; /* List of intermolecular interactions
+ * using system wide atom indices,
+ * either NULL or size F_NRE */
int natoms;
- int maxres_renum; /* Parameter for residue numbering */
- int maxresnr; /* The maximum residue number in moltype */
- t_atomtypes atomtypes; /* Atomtype properties */
- t_block mols; /* The molecules */
+ int maxres_renum; /* Parameter for residue numbering */
+ int maxresnr; /* The maximum residue number in moltype */
+ t_atomtypes atomtypes; /* Atomtype properties */
+ t_block mols; /* The molecules */
gmx_groups_t groups;
- t_symtab symtab; /* The symbol table */
+ t_symtab symtab; /* The symbol table */
} gmx_mtop_t;
/* The mdrun node-local topology struct, completely written out */
/* The old topology struct, completely written out, used in analysis tools */
typedef struct t_topology
{
- char **name; /* Name of the topology */
- t_idef idef; /* The interaction function definition */
- t_atoms atoms; /* The atoms */
- t_atomtypes atomtypes; /* Atomtype properties */
- t_block cgs; /* The charge groups */
- t_block mols; /* The molecules */
- t_blocka excls; /* The exclusions */
- t_symtab symtab; /* The symbol table */
+ char **name; /* Name of the topology */
+ t_idef idef; /* The interaction function definition */
+ t_atoms atoms; /* The atoms */
+ t_atomtypes atomtypes; /* Atomtype properties */
+ t_block cgs; /* The charge groups */
+ t_block mols; /* The molecules */
+ gmx_bool bIntermolecularInteractions; /* Inter.mol. int. ? */
+ t_blocka excls; /* The exclusions */
+ t_symtab symtab; /* The symbol table */
} t_topology;
void init_mtop(gmx_mtop_t *mtop);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff