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<int> 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))
}
/*! \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. */
* 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.
/* 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)
{
* 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).
*/
}
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;
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
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.
*/
/* 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)
{
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");
- }
}
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
} /* mk_QMMMrec */
#endif
-std::vector<int> qmmmAtomIndices(const t_inputrec& ir, const gmx_mtop_t& mtop)
-{
- const int numQmmmGroups = ir.opts.ngQM;
- const SimulationGroups& groups = mtop.groups;
- std::vector<int> 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<const int> 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
* independent neighbourlists. I have to think about
* that.. 11-11-2003
*/
- std::vector<int> qmmmAtoms = qmmmAtomIndices(*ir, *mtop);
+ std::vector<int> qmmmAtoms;
snew(qr->qm, numQmmmGroups);
for (int i = 0; i < numQmmmGroups; i++)
{
* 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<int> 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<const int> qmmmAtoms);
-
#endif
//! \brief QM/MM mode
enum struct GmxQmmmMode
{
- GMX_QMMM_ORIGINAL,
GMX_QMMM_MIMIC
};
#endif /* GMX_MDTYPES_MD_ENUMS_H */
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;
*/
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);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff