From: Christian Blau Date: Fri, 13 Mar 2020 17:02:56 +0000 (+0100) Subject: Remove dysfunctional QMMM interface pt1 X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=e67f9e8c189137e044db1ccedc61801d71c11721;p=alexxy%2Fgromacs.git Remove dysfunctional QMMM interface pt1 Removing some of the by now obsolete QMMM functionality to ease refactoring work on grompp and mtop Change-Id: I4a973e905f63d06028ec5d0e73afbcbafa7f0cf5 --- diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 06185c4d44..0ba7a763c5 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2227,20 +2227,9 @@ int gmx_grompp(int argc, char* argv[]) pr_symtab(debug, 0, "After close", &sys.symtab); } - /* make exclusions between QM atoms and remove charges if needed */ - if (ir->bQMMM) - { - generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_ORIGINAL, logger); - if (ir->QMMMscheme != eQMMMschemeoniom) - { - std::vector qmmmAtoms = qmmmAtomIndices(*ir, sys); - removeQmmmAtomCharges(&sys, qmmmAtoms); - } - } - if (ir->eI == eiMimic) { - generate_qmexcl(&sys, ir, wi, GmxQmmmMode::GMX_QMMM_MIMIC, logger); + generate_qmexcl(&sys, ir, logger); } if (ftp2bSet(efTRN, NFILE, fnm)) diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index 377bb4cb33..c03802eef9 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -1013,23 +1013,20 @@ char** do_top(bool bVerbose, } /*! \brief - * Generate exclusion lists for QM/MM. + * Exclude molecular interactions for QM atoms handled by MiMic. * - * This routine updates the exclusion lists for QM atoms in order to include all other QM - * atoms of this molecule. Moreover, this routine replaces bonds between QM atoms with - * CONNBOND and, when MiMiC is not used, removes bonded interactions between QM and link atoms. - * Finally, in case if MiMiC QM/MM is used - charges of QM atoms are set to 0 + * Update the exclusion lists to include all QM atoms of this molecule, + * replace bonds between QM atoms with CONNBOND and + * set charges of QM atoms to 0. * * \param[in,out] molt molecule type with QM atoms * \param[in] grpnr group informatio * \param[in,out] ir input record - * \param[in,out] qmmmMode QM/MM mode switch: original/MiMiC * \param[in] logger Handle to logging interface. */ static void generate_qmexcl_moltype(gmx_moltype_t* molt, const unsigned char* grpnr, t_inputrec* ir, - GmxQmmmMode qmmmMode, const gmx::MDLogger& logger) { /* This routine expects molt->ilist to be of size F_NRE and ordered. */ @@ -1038,9 +1035,9 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, * these interactions should be handled by the QM subroutines and * not by the gromacs routines */ - int qm_max = 0, qm_nr = 0, link_nr = 0, link_max = 0; + int qm_max = 0, qm_nr = 0, link_nr = 0; int * qm_arr = nullptr, *link_arr = nullptr; - bool *bQMMM, *blink; + bool* bQMMM; /* First we search and select the QM atoms in an qm_arr array that * we use to create the exclusions. @@ -1167,14 +1164,7 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, /* MiMiC treats link atoms as quantum atoms - therefore * we do not need do additional exclusions here */ - if (qmmmMode == GmxQmmmMode::GMX_QMMM_MIMIC) - { - bexcl = numQmAtoms == nratoms; - } - else - { - bexcl = (numQmAtoms >= nratoms - 1); - } + bexcl = numQmAtoms == nratoms; if (bexcl && ftype == F_SETTLE) { @@ -1207,51 +1197,6 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, * linkatoms interaction with the QMatoms and would be counted * twice. */ - if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC) - { - for (int i = 0; i < F_NRE; i++) - { - if (IS_CHEMBOND(i)) - { - int j = 0; - while (j < molt->ilist[i].size()) - { - int a1 = molt->ilist[i].iatoms[j + 1]; - int a2 = molt->ilist[i].iatoms[j + 2]; - if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2])) - { - if (link_nr >= link_max) - { - link_max += 10; - srenew(link_arr, link_max); - } - if (bQMMM[a1]) - { - link_arr[link_nr++] = a2; - } - else - { - link_arr[link_nr++] = a1; - } - } - j += 3; - } - } - } - } - snew(blink, molt->atoms.nr); - for (int i = 0; i < molt->atoms.nr; i++) - { - blink[i] = FALSE; - } - - if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC) - { - for (int i = 0; i < link_nr; i++) - { - blink[link_arr[i]] = TRUE; - } - } /* creating the exclusion block for the QM atoms. Each QM atom has * as excluded elements all the other QMatoms (and itself). */ @@ -1276,14 +1221,6 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, } j += (qm_nr + link_nr); } - if (blink[i]) - { - for (int k = 0; k < qm_nr; k++) - { - qmexcl.a[k + j] = qm_arr[k]; - } - j += qm_nr; - } } qmexcl.index[qmexcl.nr] = j; @@ -1305,10 +1242,9 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, int j = 0; while (j < molt->ilist[i].size()) { - int a1 = molt->ilist[i].iatoms[j + 1]; - int a2 = molt->ilist[i].iatoms[j + 2]; - bool bexcl = - ((bQMMM[a1] && bQMMM[a2]) || (blink[a1] && bQMMM[a2]) || (bQMMM[a1] && blink[a2])); + int a1 = molt->ilist[i].iatoms[j + 1]; + int a2 = molt->ilist[i].iatoms[j + 2]; + bool bexcl = (bQMMM[a1] && bQMMM[a2]); if (bexcl) { /* since the interaction involves QM atoms, these should be @@ -1331,10 +1267,9 @@ static void generate_qmexcl_moltype(gmx_moltype_t* molt, free(qm_arr); free(bQMMM); free(link_arr); - free(blink); } /* generate_qmexcl */ -void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, warninp* wi, GmxQmmmMode qmmmMode, const gmx::MDLogger& logger) +void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, const gmx::MDLogger& logger) { /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered. */ @@ -1411,7 +1346,7 @@ void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, warninp* wi, GmxQmmmMode q /* Set the molecule type for the QMMM molblock */ molb->type = sys->moltype.size() - 1; } - generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, qmmmMode, logger); + generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, logger); } if (grpnr) { @@ -1420,24 +1355,4 @@ void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, warninp* wi, GmxQmmmMode q index_offset += nat_mol; } } - if (qmmmMode == GmxQmmmMode::GMX_QMMM_ORIGINAL && nr_mol_with_qm_atoms > 1) - { - /* generate a warning is there are QM atoms in different - * topologies. In this case it is not possible at this stage to - * mutualy exclude the non-bonded interactions via the - * exclusions (AFAIK). Instead, the user is advised to use the - * energy group exclusions in the mdp file - */ - warning_note(wi, - "\nThe QM subsystem is divided over multiple topologies. " - "The mutual non-bonded interactions cannot be excluded. " - "There are two ways to achieve this:\n\n" - "1) merge the topologies, such that the atoms of the QM " - "subsystem are all present in one single topology file. " - "In this case this warning will dissappear\n\n" - "2) exclude the non-bonded interactions explicitly via the " - "energygrp-excl option in the mdp file. if this is the case " - "this warning may be ignored" - "\n\n"); - } } diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h index d1177212c5..48c7508fd0 100644 --- a/src/gromacs/gmxpreprocess/topio.h +++ b/src/gromacs/gmxpreprocess/topio.h @@ -85,6 +85,6 @@ char** do_top(bool bVerbose, const gmx::MDLogger& logger); /* 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, GmxQmmmMode qmmmMode, const gmx::MDLogger& logger); +void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, const gmx::MDLogger& logger); #endif diff --git a/src/gromacs/mdlib/qmmm.cpp b/src/gromacs/mdlib/qmmm.cpp index 5de20a019e..508f87cef3 100644 --- a/src/gromacs/mdlib/qmmm.cpp +++ b/src/gromacs/mdlib/qmmm.cpp @@ -353,69 +353,6 @@ t_QMMMrec* mk_QMMMrec() } /* mk_QMMMrec */ #endif -std::vector qmmmAtomIndices(const t_inputrec& ir, const gmx_mtop_t& mtop) -{ - const int numQmmmGroups = ir.opts.ngQM; - const SimulationGroups& groups = mtop.groups; - std::vector qmmmAtoms; - for (int i = 0; i < numQmmmGroups; i++) - { - for (const AtomProxy atomP : AtomRange(mtop)) - { - int index = atomP.globalAtomNumber(); - if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, index) == i) - { - qmmmAtoms.push_back(index); - } - } - if (ir.QMMMscheme == eQMMMschemeoniom) - { - /* I assume that users specify the QM groups from small to - * big(ger) in the mdp file - */ - gmx_mtop_ilistloop_all_t iloop = gmx_mtop_ilistloop_all_init(&mtop); - int nral1 = 1 + NRAL(F_VSITE2); - int atomOffset = 0; - while (const InteractionLists* ilists = gmx_mtop_ilistloop_all_next(iloop, &atomOffset)) - { - const InteractionList& ilist = (*ilists)[F_VSITE2]; - for (int j = 0; j < ilist.size(); j += nral1) - { - const int vsite = atomOffset + ilist.iatoms[j]; /* the vsite */ - const int ai = atomOffset + ilist.iatoms[j + 1]; /* constructing atom */ - const int aj = atomOffset + ilist.iatoms[j + 2]; /* constructing atom */ - if (getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite) - == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, ai) - && getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, vsite) - == getGroupType(groups, SimulationAtomGroupType::QuantumMechanics, aj)) - { - /* this dummy link atom needs to be removed from qmmmAtoms - * before making the QMrec of this layer! - */ - qmmmAtoms.erase(std::remove_if(qmmmAtoms.begin(), qmmmAtoms.end(), - [&vsite](int atom) { return atom == vsite; }), - qmmmAtoms.end()); - } - } - } - } - } - return qmmmAtoms; -} - -void removeQmmmAtomCharges(gmx_mtop_t* mtop, gmx::ArrayRef qmmmAtoms) -{ - int molb = 0; - for (gmx::index i = 0; i < qmmmAtoms.ssize(); i++) - { - int indexInMolecule; - mtopGetMolblockIndex(mtop, qmmmAtoms[i], &molb, nullptr, &indexInMolecule); - t_atom* atom = &mtop->moltype[mtop->molblock[molb].type].atoms.atom[indexInMolecule]; - atom->q = 0.0; - atom->qB = 0.0; - } -} - void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* ir, const t_forcerec* fr) { /* we put the atomsnumbers of atoms that belong to the QMMM group in @@ -484,7 +421,7 @@ void init_QMMMrec(const t_commrec* cr, const gmx_mtop_t* mtop, const t_inputrec* * independent neighbourlists. I have to think about * that.. 11-11-2003 */ - std::vector qmmmAtoms = qmmmAtomIndices(*ir, *mtop); + std::vector qmmmAtoms; snew(qr->qm, numQmmmGroups); for (int i = 0; i < numQmmmGroups; i++) { diff --git a/src/gromacs/mdlib/qmmm.h b/src/gromacs/mdlib/qmmm.h index dce3ed8941..0eda21d216 100644 --- a/src/gromacs/mdlib/qmmm.h +++ b/src/gromacs/mdlib/qmmm.h @@ -149,21 +149,4 @@ real calculate_QMMM(const t_commrec* cr, gmx::ForceWithShiftForces* forceWithShi * called by system(). */ -/*! \brief - * Return vector of atom indices for atoms in the QMMM region. - * - * \param[in] mtop Topology to use for populating array. - * \param[in] ir Inputrec used in simulation. - * \returns Vector of atoms. - */ -std::vector qmmmAtomIndices(const t_inputrec& ir, const gmx_mtop_t& mtop); - -/*! \brief - * Remove charges from QMMM atoms. - * - * \param[in] mtop Topology used for removing atoms. - * \param[in] qmmmAtoms ArrayRef to vector conatining qmmm atom indices. - */ -void removeQmmmAtomCharges(gmx_mtop_t* mtop, gmx::ArrayRef qmmmAtoms); - #endif diff --git a/src/gromacs/mdtypes/md_enums.h b/src/gromacs/mdtypes/md_enums.h index bd66df5857..7f9bf5f3d8 100644 --- a/src/gromacs/mdtypes/md_enums.h +++ b/src/gromacs/mdtypes/md_enums.h @@ -891,7 +891,6 @@ extern const char* gmx_nblist_interaction_names[GMX_NBLIST_INTERACTION_NR + 1]; //! \brief QM/MM mode enum struct GmxQmmmMode { - GMX_QMMM_ORIGINAL, GMX_QMMM_MIMIC }; #endif /* GMX_MDTYPES_MD_ENUMS_H */ diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index 13c79f817c..f358e51437 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -406,68 +406,6 @@ typedef struct gmx_mtop_ilistloop_all int a_offset; } t_gmx_mtop_ilist_all; -gmx_mtop_ilistloop_all_t gmx_mtop_ilistloop_all_init(const gmx_mtop_t* mtop) -{ - struct gmx_mtop_ilistloop_all* iloop; - - snew(iloop, 1); - - iloop->mtop = mtop; - iloop->mblock = 0; - iloop->mol = -1; - iloop->a_offset = 0; - - return iloop; -} - -static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop) -{ - sfree(iloop); -} - -const InteractionLists* gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, int* atnr_offset) -{ - - if (iloop == nullptr) - { - gmx_incons( - "gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init"); - } - - if (iloop->mol >= 0) - { - iloop->a_offset += iloop->mtop->moleculeBlockIndices[iloop->mblock].numAtomsPerMolecule; - } - - iloop->mol++; - - /* 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->molblock.size() - || iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol) - { - iloop->mblock++; - iloop->mol = 0; - if (iloop->mblock >= iloop->mtop->molblock.size()) - { - if (iloop->mblock == iloop->mtop->molblock.size() && iloop->mtop->bIntermolecularInteractions) - { - *atnr_offset = 0; - return iloop->mtop->intermolecular_ilist.get(); - } - - gmx_mtop_ilistloop_all_destroy(iloop); - return nullptr; - } - } - - *atnr_offset = iloop->a_offset; - - return &iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist; -} - int gmx_mtop_ftype_count(const gmx_mtop_t* mtop, int ftype) { gmx_mtop_ilistloop_t iloop; diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index 2843f28c14..24fef841e9 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -217,25 +217,6 @@ gmx_mtop_ilistloop_t gmx_mtop_ilistloop_init(const gmx_mtop_t& mtop); */ const InteractionLists* gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, int* nmol); -/* Abstract type for ilist loop over all ilists of all molecules */ -typedef struct gmx_mtop_ilistloop_all* gmx_mtop_ilistloop_all_t; - -/* Initialize an ilist loop over all molecule types in the system. - * Only use this when you really need to loop over all molecules, - * i.e. when you use groups which might differ per molecule, - * otherwise use gmx_mtop_ilistloop. - */ -gmx_mtop_ilistloop_all_t gmx_mtop_ilistloop_all_init(const gmx_mtop_t* mtop); - -/* Loop to the next molecule, - * When not at the end: - * returns a valid pointer to the next array ilist_mol[F_NRE], - * writes the atom offset which should be added to iatoms in atnr_offset. - * When at the end, destroys iloop and returns nullptr. - */ -const InteractionLists* gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, int* atnr_offset); - - /* Returns the total number of interactions in the system of type ftype */ int gmx_mtop_ftype_count(const gmx_mtop_t* mtop, int ftype);