From: Berk Hess Date: Wed, 21 Aug 2013 15:30:36 +0000 (+0200) Subject: Intermolecular bonded interaction support added X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=6ab8ccf204441da0dcb138f55201027ac10a4091;p=alexxy%2Fgromacs.git Intermolecular bonded interaction support added The .top file can now contain an [ intermolecular_interactions ] directive, after which bonded interactions can be entered using global atom indices. Added a molecule defition section to the topologies chapter in the manual. A description of the moleculetype directive was missing. This section has a subsection on intermolecular interactions. Change-Id: I383287dd0729fef1c54f27a4fe9f5d628445549c --- diff --git a/docs/manual/topology.tex b/docs/manual/topology.tex index 9646061882..2f5b102015 100644 --- a/docs/manual/topology.tex +++ b/docs/manual/topology.tex @@ -1,7 +1,7 @@ % % 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. @@ -409,6 +409,41 @@ want to include parameters for new atom types, make sure you define 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} @@ -477,7 +512,35 @@ Energies for types 1 and 2 are written to the energy and log file 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: @@ -508,33 +571,6 @@ as proposed by Hawkins {\em et al.}~\cite{Truhlar96}, but on a per-atom (rather 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. @@ -1169,10 +1205,17 @@ in \tabref{topfile2}} \\ {\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}\\ @@ -1202,7 +1245,7 @@ Name of interaction & Topology file directive & num. & fu % 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} \\ diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 2cb8e7a64d..4f31d1fe54 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -6740,7 +6740,8 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr, 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); diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index 6795e28145..f0863e061c 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -45,6 +45,7 @@ #include "gmxpre.h" +#include #include #include @@ -103,17 +104,20 @@ typedef struct { /*! \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 */ @@ -503,8 +507,9 @@ static void count_excls(const t_block *cgs, const t_blocka *excls, } /*! \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, @@ -512,12 +517,12 @@ static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom, 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++) @@ -606,8 +611,9 @@ static int low_make_reverse_ilist(t_ilist *il_mt, t_atom *atom, } /*! \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, @@ -616,9 +622,9 @@ static int make_reverse_ilist(gmx_moltype_t *molt, 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); @@ -634,7 +640,7 @@ static int make_reverse_ilist(gmx_moltype_t *molt, /* 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, @@ -685,7 +691,8 @@ static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE, /* 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]); @@ -703,6 +710,25 @@ static gmx_reverse_top_t *make_reverse_top(gmx_mtop_t *mtop, gmx_bool bFE, } 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; @@ -752,7 +778,7 @@ void dd_make_reverse_top(FILE *fplog, /* 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, @@ -1000,10 +1026,10 @@ static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb, } /*! \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; @@ -1269,228 +1295,324 @@ static void combine_idef(t_idef *dest, const thread_work_t *src, int nsrc, } } -/*! \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); } } @@ -1811,7 +1933,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, 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; @@ -1854,10 +1976,10 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, 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++) @@ -1907,10 +2029,10 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, 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) { @@ -1930,7 +2052,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, make_exclusions_zone(dd, zones, mtop->moltype, cginfo, excl_t, - iz, + izone, cg0t, cg1t); } else @@ -1940,7 +2062,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, mtop->moltype, bRCheck2B, rc2, la2lc, pbc_null, cg_cm, cginfo, excl_t, - iz, + izone, cg0t, cg1t); } } @@ -1956,7 +2078,7 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, nbonded_local += rt->th_work[thread].nbonded; } - if (iz < nzone_excl) + if (izone < nzone_excl) { if (rt->nthread > 1) { @@ -1973,9 +2095,9 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd, /* 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); @@ -2203,14 +2325,14 @@ static int *make_at2cg(t_block *cgs) 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; @@ -2219,7 +2341,18 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, * 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); @@ -2244,55 +2377,80 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, * 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); @@ -2302,12 +2460,12 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, 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++) { @@ -2330,6 +2488,11 @@ t_blocka *make_charge_group_links(gmx_mtop_t *mtop, gmx_domdec_t *dd, } } + 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); diff --git a/src/gromacs/fileio/tpxio.c b/src/gromacs/fileio/tpxio.c index cf77d95e6b..26aed407ee 100644 --- a/src/gromacs/fileio/tpxio.c +++ b/src/gromacs/fileio/tpxio.c @@ -92,7 +92,8 @@ enum tpxv { 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 @@ -106,7 +107,7 @@ enum tpxv { * * 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 @@ -3024,16 +3025,19 @@ static void add_posres_molblock(gmx_mtop_t *mtop) 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; @@ -3122,6 +3126,23 @@ static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead, 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) { diff --git a/src/gromacs/gmxlib/txtdump.c b/src/gromacs/gmxlib/txtdump.c index 4d04ff8759..606a8f4c07 100644 --- a/src/gromacs/gmxlib/txtdump.c +++ b/src/gromacs/gmxlib/txtdump.c @@ -1882,7 +1882,7 @@ static void pr_molblock(FILE *fp, int indent, const char *title, 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)) { @@ -1895,6 +1895,16 @@ void pr_mtop(FILE *fp, int indent, const char *title, gmx_mtop_t *mtop, { 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++) @@ -1917,6 +1927,7 @@ void pr_top(FILE *fp, int indent, const char *title, t_topology *top, gmx_bool b 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); } diff --git a/src/gromacs/gmxpreprocess/convparm.c b/src/gromacs/gmxpreprocess/convparm.c index e3bac95963..93f1603c8e 100644 --- a/src/gromacs/gmxpreprocess/convparm.c +++ b/src/gromacs/gmxpreprocess/convparm.c @@ -3,7 +3,7 @@ * * 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. @@ -567,7 +567,8 @@ static void enter_function(t_params *p, t_functype ftype, int comb, real reppow, } 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; @@ -611,6 +612,57 @@ void convert_params(int atnr, t_params nbtypes[], } } } + + 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", diff --git a/src/gromacs/gmxpreprocess/convparm.h b/src/gromacs/gmxpreprocess/convparm.h index d8c69a250d..0505967510 100644 --- a/src/gromacs/gmxpreprocess/convparm.h +++ b/src/gromacs/gmxpreprocess/convparm.h @@ -3,7 +3,7 @@ * * 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. @@ -46,7 +46,9 @@ extern "C" { #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 diff --git a/src/gromacs/gmxpreprocess/grompp-impl.h b/src/gromacs/gmxpreprocess/grompp-impl.h index 6c1093c50f..20cafb0cd7 100644 --- a/src/gromacs/gmxpreprocess/grompp-impl.h +++ b/src/gromacs/gmxpreprocess/grompp-impl.h @@ -3,7 +3,7 @@ * * 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. @@ -165,6 +165,7 @@ typedef enum { d_orientation_restraints, d_dihedral_restraints, d_cmap, + d_intermolecular_interactions, d_maxdir, d_invalid, d_none diff --git a/src/gromacs/gmxpreprocess/grompp.c b/src/gromacs/gmxpreprocess/grompp.c index 30881cbf29..2c9ed586d7 100644 --- a/src/gromacs/gmxpreprocess/grompp.c +++ b/src/gromacs/gmxpreprocess/grompp.c @@ -499,7 +499,8 @@ new_status(const char *topfile, const char *topppfile, const char *confin, 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) @@ -524,7 +525,8 @@ new_status(const char *topfile, const char *topppfile, const char *confin, /* 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); @@ -1510,7 +1512,7 @@ int gmx_grompp(int argc, char *argv[]) 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; @@ -1637,7 +1639,8 @@ int gmx_grompp(int argc, char *argv[]) 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); @@ -1806,7 +1809,8 @@ int gmx_grompp(int argc, char *argv[]) } 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) { diff --git a/src/gromacs/gmxpreprocess/topdirs.c b/src/gromacs/gmxpreprocess/topdirs.c index 1c42dfab2d..a1e14b270d 100644 --- a/src/gromacs/gmxpreprocess/topdirs.c +++ b/src/gromacs/gmxpreprocess/topdirs.c @@ -3,7 +3,7 @@ * * 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. @@ -86,6 +86,7 @@ static const char *directive_names[d_maxdir+1] = { "orientation_restraints", "dihedral_restraints", "cmap", + "intermolecular_interactions", "invalid" }; @@ -396,6 +397,7 @@ void DS_Init(DirStack **DS) 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++) { diff --git a/src/gromacs/gmxpreprocess/topio.c b/src/gromacs/gmxpreprocess/topio.c index 7a3fe0033f..3478619757 100644 --- a/src/gromacs/gmxpreprocess/topio.c +++ b/src/gromacs/gmxpreprocess/topio.c @@ -528,6 +528,31 @@ generate_gb_exclusion_interactions(t_molinfo *mi, gpp_atomtype_t atype, t_nextnb } +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, @@ -536,6 +561,7 @@ 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, @@ -562,6 +588,7 @@ static char **read_topol(const char *infile, const char *outfile, 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; @@ -605,9 +632,10 @@ static char **read_topol(const char *infile, const char *outfile, 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; @@ -728,6 +756,22 @@ static char **read_topol(const char *infile, const char *outfile, 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); } @@ -758,8 +802,7 @@ static char **read_topol(const char *infile, const char *outfile, 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); @@ -857,14 +900,14 @@ static char **read_topol(const char *infile, const char *outfile, } 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); @@ -1055,6 +1098,11 @@ static char **read_topol(const char *infile, const char *outfile, done_bond_atomtype(&batype); + if (*intermolecular_interactions != NULL) + { + sfree(mi0->atoms.atom); + } + *nrmols = nmol; *nmolblock = nmolb; @@ -1076,6 +1124,7 @@ char **do_top(gmx_bool bVerbose, gpp_atomtype_t atype, int *nrmols, t_molinfo **molinfo, + t_molinfo **intermolecular_interactions, t_inputrec *ir, int *nmolblock, gmx_molblock_t **molblock, @@ -1100,7 +1149,8 @@ char **do_top(gmx_bool bVerbose, 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); diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h index e03b102aa4..744f31202a 100644 --- a/src/gromacs/gmxpreprocess/topio.h +++ b/src/gromacs/gmxpreprocess/topio.h @@ -3,7 +3,7 @@ * * 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. @@ -63,13 +63,13 @@ char **do_top(gmx_bool bVerbose, 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); diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp index 049fc4cba7..caadeb67a3 100644 --- a/src/gromacs/mdlib/broadcaststructs.cpp +++ b/src/gromacs/mdlib/broadcaststructs.cpp @@ -813,6 +813,13 @@ void bcast_ir_mtop(const t_commrec *cr, t_inputrec *inputrec, gmx_mtop_t *mtop) 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++) diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index 1333573e30..8edb8d5879 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -2524,26 +2524,48 @@ void init_forcerec(FILE *fp, { 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"); } } } diff --git a/src/gromacs/tools/compare.c b/src/gromacs/tools/compare.c index fd73cd8c97..975b62f536 100644 --- a/src/gromacs/tools/compare.c +++ b/src/gromacs/tools/compare.c @@ -451,6 +451,7 @@ static void cmp_top(FILE *fp, t_topology *t1, t_topology *t2, real ftol, real ab 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 diff --git a/src/gromacs/topology/mtop_util.c b/src/gromacs/topology/mtop_util.c index 24eb127af2..771ccd68ad 100644 --- a/src/gromacs/topology/mtop_util.c +++ b/src/gromacs/topology/mtop_util.c @@ -1,7 +1,7 @@ /* * 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. @@ -625,8 +625,16 @@ gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, } 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; } @@ -683,12 +691,25 @@ gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, 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; } @@ -716,6 +737,11 @@ int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype) n += nmol*il[ftype].nr/(1+NRAL(ftype)); } + if (mtop->bIntermolecularInteractions) + { + n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype)); + } + return n; } @@ -1112,6 +1138,15 @@ static void gen_local_top(const gmx_mtop_t *mtop, const t_inputrec *ir, 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; @@ -1159,14 +1194,15 @@ t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop) 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. diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index 6f4e4792d4..85d32c5cd1 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -3,7 +3,7 @@ * * 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. @@ -56,11 +56,11 @@ enum { 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 @@ -99,13 +99,18 @@ typedef struct gmx_mtop_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 */ @@ -120,14 +125,15 @@ typedef struct gmx_localtop_t /* 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);