2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
53 #include <unordered_set>
54 #include <sys/types.h>
56 #include "gromacs/fileio/gmxfio.h"
57 #include "gromacs/fileio/warninp.h"
58 #include "gromacs/gmxpreprocess/gmxcpp.h"
59 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
60 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
61 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
62 #include "gromacs/gmxpreprocess/grompp_impl.h"
63 #include "gromacs/gmxpreprocess/readir.h"
64 #include "gromacs/gmxpreprocess/topdirs.h"
65 #include "gromacs/gmxpreprocess/toppush.h"
66 #include "gromacs/gmxpreprocess/topshake.h"
67 #include "gromacs/gmxpreprocess/toputil.h"
68 #include "gromacs/gmxpreprocess/vsite_parm.h"
69 #include "gromacs/math/units.h"
70 #include "gromacs/math/utilities.h"
71 #include "gromacs/mdtypes/inputrec.h"
72 #include "gromacs/mdtypes/md_enums.h"
73 #include "gromacs/pbcutil/pbc.h"
74 #include "gromacs/topology/block.h"
75 #include "gromacs/topology/exclusionblocks.h"
76 #include "gromacs/topology/ifunc.h"
77 #include "gromacs/topology/symtab.h"
78 #include "gromacs/topology/topology.h"
79 #include "gromacs/utility/cstringutil.h"
80 #include "gromacs/utility/fatalerror.h"
81 #include "gromacs/utility/futil.h"
82 #include "gromacs/utility/gmxassert.h"
83 #include "gromacs/utility/pleasecite.h"
84 #include "gromacs/utility/smalloc.h"
86 #define OPENDIR '[' /* starting sign for directive */
87 #define CLOSEDIR ']' /* ending sign for directive */
89 static void gen_pairs(const InteractionsOfType& nbs, InteractionsOfType* pairs, real fudge, int comb)
93 int nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
94 GMX_ASSERT(nnn * nnn == ntp,
95 "Number of pairs of generated non-bonded parameters should be a perfect square");
96 int nrfp = NRFP(F_LJ);
97 int nrfpA = interaction_function[F_LJ14].nrfpA;
98 int nrfpB = interaction_function[F_LJ14].nrfpB;
100 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
102 gmx_incons("Number of force parameters in gen_pairs wrong");
105 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
106 pairs->interactionTypes.clear();
108 std::array<int, 2> atomNumbers;
109 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
110 for (const auto& type : nbs.interactionTypes)
112 /* Copy type.atoms */
113 atomNumbers = { i / nnn, i % nnn };
114 /* Copy normal and FEP parameters and multiply by fudge factor */
115 gmx::ArrayRef<const real> existingParam = type.forceParam();
116 GMX_RELEASE_ASSERT(2 * nrfp <= MAXFORCEPARAM,
117 "Can't have more parameters than half of maximum p arameter number");
118 for (int j = 0; j < nrfp; j++)
120 /* If we are using sigma/epsilon values, only the epsilon values
121 * should be scaled, but not sigma.
122 * The sigma values have even indices 0,2, etc.
124 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j % 2 == 0))
133 forceParam[j] = scaling * existingParam[j];
134 forceParam[nrfp + j] = scaling * existingParam[j];
136 pairs->interactionTypes.emplace_back(InteractionOfType(atomNumbers, forceParam));
141 double check_mol(const gmx_mtop_t* mtop, warninp* wi)
148 /* Check mass and charge */
151 for (const gmx_molblock_t& molb : mtop->molblock)
153 const t_atoms* atoms = &mtop->moltype[molb.type].atoms;
154 for (i = 0; (i < atoms->nr); i++)
156 q += molb.nmol * atoms->atom[i].q;
157 m = atoms->atom[i].m;
158 mB = atoms->atom[i].mB;
159 pt = atoms->atom[i].ptype;
160 /* If the particle is an atom or a nucleus it must have a mass,
161 * else, if it is a shell, a vsite or a bondshell it can have mass zero
163 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
165 ri = atoms->atom[i].resind;
166 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
167 *(atoms->atomname[i]), *(atoms->resinfo[ri].name), atoms->resinfo[ri].nr, m, mB);
168 warning_error(wi, buf);
170 else if (((m != 0) || (mB != 0)) && (pt == eptVSite))
172 ri = atoms->atom[i].resind;
174 "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state "
176 " Check your topology.\n",
177 *(atoms->atomname[i]), *(atoms->resinfo[ri].name), atoms->resinfo[ri].nr, m, mB);
178 warning_error(wi, buf);
179 /* The following statements make LINCS break! */
180 /* atoms->atom[i].m=0; */
187 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
189 * The results of this routine are only used for checking and for
190 * printing warning messages. Thus we can assume that charges of molecules
191 * should be integer. If the user wanted non-integer molecular charge,
192 * an undesired warning is printed and the user should use grompp -maxwarn 1.
194 * \param qMol The total, unrounded, charge of the molecule
195 * \param sumAbsQ The sum of absolute values of the charges, used for determining the tolerance for the rounding.
197 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
199 /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
200 * of the charges for ascii float truncation in the topology files.
201 * Although the summation here uses double precision, the charges
202 * are read and stored in single precision when real=float. This can
203 * lead to rounding errors of half the least significant bit.
204 * Note that, unfortunately, we can not assume addition of random
205 * rounding errors. It is not entirely unlikely that many charges
206 * have a near half-bit rounding error with the same sign.
208 double tolAbs = 1e-6;
209 double tol = std::max(tolAbs, 0.5 * GMX_REAL_EPS * sumAbsQ);
210 double qRound = std::round(qMol);
211 if (std::abs(qMol - qRound) <= tol)
221 static void sum_q(const t_atoms* atoms, int numMols, double* qTotA, double* qTotB)
228 for (int i = 0; i < atoms->nr; i++)
230 qmolA += atoms->atom[i].q;
231 qmolB += atoms->atom[i].qB;
232 sumAbsQA += std::abs(atoms->atom[i].q);
233 sumAbsQB += std::abs(atoms->atom[i].qB);
236 *qTotA += numMols * roundedMoleculeCharge(qmolA, sumAbsQA);
237 *qTotB += numMols * roundedMoleculeCharge(qmolB, sumAbsQB);
240 static void get_nbparm(char* nb_str, char* comb_str, int* nb, int* comb, warninp* wi)
243 char warn_buf[STRLEN];
246 for (i = 1; (i < eNBF_NR); i++)
248 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
255 *nb = strtol(nb_str, nullptr, 10);
257 if ((*nb < 1) || (*nb >= eNBF_NR))
259 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s", nb_str, enbf_names[1]);
260 warning_error(wi, warn_buf);
264 for (i = 1; (i < eCOMB_NR); i++)
266 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
273 *comb = strtol(comb_str, nullptr, 10);
275 if ((*comb < 1) || (*comb >= eCOMB_NR))
277 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s", comb_str, ecomb_names[1]);
278 warning_error(wi, warn_buf);
283 static char** cpp_opts(const char* define, const char* include, warninp* wi)
287 const char* cppadds[2];
288 char** cppopts = nullptr;
289 const char* option[2] = { "-D", "-I" };
290 const char* nopt[2] = { "define", "include" };
294 char warn_buf[STRLEN];
297 cppadds[1] = include;
298 for (n = 0; (n < 2); n++)
305 while ((*ptr != '\0') && isspace(*ptr))
310 while ((*rptr != '\0') && !isspace(*rptr))
317 snew(buf, (len + 1));
318 strncpy(buf, ptr, len);
319 if (strstr(ptr, option[n]) != ptr)
321 set_warning_line(wi, "mdp file", -1);
322 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
323 warning(wi, warn_buf);
327 srenew(cppopts, ++ncppopts);
328 cppopts[ncppopts - 1] = gmx_strdup(buf);
336 srenew(cppopts, ++ncppopts);
337 cppopts[ncppopts - 1] = nullptr;
343 static void make_atoms_sys(gmx::ArrayRef<const gmx_molblock_t> molblock,
344 gmx::ArrayRef<const MoleculeInformation> molinfo,
348 atoms->atom = nullptr;
350 for (const gmx_molblock_t& molb : molblock)
352 const t_atoms& mol_atoms = molinfo[molb.type].atoms;
354 srenew(atoms->atom, atoms->nr + molb.nmol * mol_atoms.nr);
356 for (int m = 0; m < molb.nmol; m++)
358 for (int a = 0; a < mol_atoms.nr; a++)
360 atoms->atom[atoms->nr++] = mol_atoms.atom[a];
367 static char** read_topol(const char* infile,
372 PreprocessingAtomTypes* atypes,
373 std::vector<MoleculeInformation>* molinfo,
374 std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
375 gmx::ArrayRef<InteractionsOfType> interactions,
376 int* combination_rule,
380 std::vector<gmx_molblock_t>* molblock,
381 bool* ffParametrizedWithHBondConstraints,
384 bool usingFullRangeElectrostatics,
389 char * pline = nullptr, **title = nullptr;
390 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
392 char * dirstr, *dummy2;
393 int nrcopies, nscan, ncombs, ncopy;
394 double fLJ, fQQ, fPOW;
395 MoleculeInformation* mi0 = nullptr;
398 t_nbparam ** nbparam, **pair;
399 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
400 bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
401 double qt = 0, qBt = 0; /* total charge */
402 int dcatt = -1, nmol_couple;
403 /* File handling variables */
407 char* tmp_line = nullptr;
408 char warn_buf[STRLEN];
409 const char* floating_point_arithmetic_tip =
410 "Total charge should normally be an integer. See\n"
411 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
412 "for discussion on how close it should be to an integer.\n";
413 /* We need to open the output file before opening the input file,
414 * because cpp_open_file can change the current working directory.
418 out = gmx_fio_fopen(outfile, "w");
425 /* open input file */
426 auto cpp_opts_return = cpp_opts(define, include, wi);
427 status = cpp_open_file(infile, &handle, cpp_opts_return);
430 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
433 /* some local variables */
434 DS_Init(&DS); /* directive stack */
435 d = Directive::d_invalid; /* first thing should be a directive */
436 nbparam = nullptr; /* The temporary non-bonded matrix */
437 pair = nullptr; /* The temporary pair interaction matrix */
438 std::vector<std::vector<gmx::ExclusionBlock>> exclusionBlocks;
441 *reppow = 12.0; /* Default value for repulsion power */
443 /* Init the number of CMAP torsion angles and grid spacing */
444 interactions[F_CMAP].cmakeGridSpacing = 0;
445 interactions[F_CMAP].cmapAngles = 0;
447 bWarn_copy_A_B = bFEP;
449 PreprocessingBondAtomType bondAtomType;
450 /* parse the actual file */
451 bReadDefaults = FALSE;
453 bReadMolType = FALSE;
458 status = cpp_read_line(&handle, STRLEN, line);
459 done = (status == eCPP_EOF);
462 if (status != eCPP_OK)
464 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
468 fprintf(out, "%s\n", line);
471 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
473 pline = gmx_strdup(line);
475 /* Strip trailing '\' from pline, if it exists */
477 if ((sl > 0) && (pline[sl - 1] == CONTINUE))
482 /* build one long line from several fragments - necessary for CMAP */
483 while (continuing(line))
485 status = cpp_read_line(&handle, STRLEN, line);
486 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
488 /* Since we depend on the '\' being present to continue to read, we copy line
489 * to a tmp string, strip the '\' from that string, and cat it to pline
491 tmp_line = gmx_strdup(line);
493 sl = strlen(tmp_line);
494 if ((sl > 0) && (tmp_line[sl - 1] == CONTINUE))
496 tmp_line[sl - 1] = ' ';
499 done = (status == eCPP_EOF);
502 if (status != eCPP_OK)
504 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
508 fprintf(out, "%s\n", line);
512 srenew(pline, strlen(pline) + strlen(tmp_line) + 1);
513 strcat(pline, tmp_line);
517 /* skip trailing and leading spaces and comment text */
518 strip_comment(pline);
521 /* if there is something left... */
522 if (static_cast<int>(strlen(pline)) > 0)
524 if (pline[0] == OPENDIR)
526 /* A directive on this line: copy the directive
527 * without the brackets into dirstr, then
528 * skip spaces and tabs on either side of directive
530 dirstr = gmx_strdup((pline + 1));
531 if ((dummy2 = strchr(dirstr, CLOSEDIR)) != nullptr)
537 if ((newd = str2dir(dirstr)) == Directive::d_invalid)
539 sprintf(errbuf, "Invalid directive %s", dirstr);
540 warning_error(wi, errbuf);
544 /* Directive found */
545 if (DS_Check_Order(DS, newd))
552 /* we should print here which directives should have
553 been present, and which actually are */
554 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
555 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
556 /* d = Directive::d_invalid; */
559 if (d == Directive::d_intermolecular_interactions)
561 if (*intermolecular_interactions == nullptr)
563 /* We (mis)use the moleculetype processing
564 * to process the intermolecular interactions
565 * by making a "molecule" of the size of the system.
567 *intermolecular_interactions = std::make_unique<MoleculeInformation>();
568 mi0 = intermolecular_interactions->get();
570 make_atoms_sys(*molblock, *molinfo, &mi0->atoms);
576 else if (d != Directive::d_invalid)
578 /* Not a directive, just a plain string
579 * use a gigantic switch to decode,
580 * if there is a valid directive!
584 case Directive::d_defaults:
587 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
588 cpp_error(&handle, eCPP_SYNTAX));
590 bReadDefaults = TRUE;
591 nscan = sscanf(pline, "%s%s%s%lf%lf%lf", nb_str, comb_str, genpairs,
603 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
606 bGenPairs = (gmx::equalCaseInsensitive(genpairs, "Y", 1));
607 if (nb_funct != eNBF_LJ && bGenPairs)
610 "Generating pair parameters is only supported "
611 "with LJ non-bonded interactions");
627 nb_funct = ifunc_index(Directive::d_nonbond_params, nb_funct);
630 case Directive::d_atomtypes:
631 push_at(symtab, atypes, &bondAtomType, pline, nb_funct, &nbparam,
632 bGenPairs ? &pair : nullptr, wi);
635 case Directive::d_bondtypes:
636 push_bt(d, interactions, 2, nullptr, &bondAtomType, pline, wi);
638 case Directive::d_constrainttypes:
639 push_bt(d, interactions, 2, nullptr, &bondAtomType, pline, wi);
641 case Directive::d_pairtypes:
644 push_nbt(d, pair, atypes, pline, F_LJ14, wi);
648 push_bt(d, interactions, 2, atypes, nullptr, pline, wi);
651 case Directive::d_angletypes:
652 push_bt(d, interactions, 3, nullptr, &bondAtomType, pline, wi);
654 case Directive::d_dihedraltypes:
655 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
656 push_dihedraltype(d, interactions, &bondAtomType, pline, wi);
659 case Directive::d_nonbond_params:
660 push_nbt(d, nbparam, atypes, pline, nb_funct, wi);
663 case Directive::d_implicit_genborn_params:
664 // Skip this line, so old topologies with
665 // GB parameters can be read.
668 case Directive::d_implicit_surface_params:
669 // Skip this line, so that any topologies
670 // with surface parameters can be read
671 // (even though these were never formally
675 case Directive::d_cmaptypes:
676 push_cmaptype(d, interactions, 5, atypes, &bondAtomType, pline, wi);
679 case Directive::d_moleculetype:
684 if (opts->couple_moltype != nullptr
685 && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam0 == ecouplamQ
686 || opts->couple_lam1 == ecouplamNONE
687 || opts->couple_lam1 == ecouplamQ))
689 dcatt = add_atomtype_decoupled(symtab, atypes, &nbparam,
690 bGenPairs ? &pair : nullptr);
692 ntype = atypes->size();
693 ncombs = (ntype * (ntype + 1)) / 2;
694 generate_nbparams(*combination_rule, nb_funct,
695 &(interactions[nb_funct]), atypes, wi);
696 ncopy = copy_nbparams(nbparam, nb_funct, &(interactions[nb_funct]), ntype);
698 "Generated %d of the %d non-bonded parameter "
700 ncombs - ncopy, ncombs);
701 free_nbparam(nbparam, ntype);
704 gen_pairs((interactions[nb_funct]), &(interactions[F_LJ14]),
705 fudgeLJ, *combination_rule);
706 ncopy = copy_nbparams(pair, nb_funct, &(interactions[F_LJ14]), ntype);
708 "Generated %d of the %d 1-4 parameter combinations\n",
709 ncombs - ncopy, ncombs);
710 free_nbparam(pair, ntype);
712 /* Copy GBSA parameters to atomtype array? */
717 push_molt(symtab, molinfo, pline, wi);
718 exclusionBlocks.emplace_back();
719 mi0 = &molinfo->back();
720 mi0->atoms.haveMass = TRUE;
721 mi0->atoms.haveCharge = TRUE;
722 mi0->atoms.haveType = TRUE;
723 mi0->atoms.haveBState = TRUE;
724 mi0->atoms.havePdbInfo = FALSE;
727 case Directive::d_atoms:
728 push_atom(symtab, &(mi0->atoms), atypes, pline, wi);
731 case Directive::d_pairs:
734 "Need to have a valid MoleculeInformation object to work on");
735 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
736 pline, FALSE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
738 case Directive::d_pairs_nb:
741 "Need to have a valid MoleculeInformation object to work on");
742 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
743 pline, FALSE, FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
746 case Directive::d_vsites2:
747 case Directive::d_vsites3:
748 case Directive::d_vsites4:
749 case Directive::d_bonds:
750 case Directive::d_angles:
751 case Directive::d_constraints:
752 case Directive::d_settles:
753 case Directive::d_position_restraints:
754 case Directive::d_angle_restraints:
755 case Directive::d_angle_restraints_z:
756 case Directive::d_distance_restraints:
757 case Directive::d_orientation_restraints:
758 case Directive::d_dihedral_restraints:
759 case Directive::d_dihedrals:
760 case Directive::d_polarization:
761 case Directive::d_water_polarization:
762 case Directive::d_thole_polarization:
765 "Need to have a valid MoleculeInformation object to work on");
766 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
767 pline, TRUE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
769 case Directive::d_cmap:
772 "Need to have a valid MoleculeInformation object to work on");
773 push_cmap(d, interactions, mi0->interactions, &(mi0->atoms), atypes, pline, wi);
776 case Directive::d_vsitesn:
779 "Need to have a valid MoleculeInformation object to work on");
780 push_vsitesn(d, mi0->interactions, &(mi0->atoms), pline, wi);
782 case Directive::d_exclusions:
783 GMX_ASSERT(!exclusionBlocks.empty(),
784 "exclusionBlocks must always be allocated so exclusions can "
786 if (exclusionBlocks.back().empty())
788 GMX_RELEASE_ASSERT(mi0,
789 "Need to have a valid MoleculeInformation "
790 "object to work on");
791 exclusionBlocks.back().resize(mi0->atoms.nr);
793 push_excl(pline, exclusionBlocks.back(), wi);
795 case Directive::d_system:
797 title = put_symtab(symtab, pline);
799 case Directive::d_molecules:
804 push_mol(*molinfo, pline, &whichmol, &nrcopies, wi);
805 mi0 = &((*molinfo)[whichmol]);
806 molblock->resize(molblock->size() + 1);
807 molblock->back().type = whichmol;
808 molblock->back().nmol = nrcopies;
810 bCouple = (opts->couple_moltype != nullptr
811 && (gmx_strcasecmp("system", opts->couple_moltype) == 0
812 || strcmp(*(mi0->name), opts->couple_moltype) == 0));
815 nmol_couple += nrcopies;
818 if (mi0->atoms.nr == 0)
820 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms", *mi0->name);
822 fprintf(stderr, "Excluding %d bonded neighbours molecule type '%s'\n",
823 mi0->nrexcl, *mi0->name);
824 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
825 if (!mi0->bProcessed)
827 generate_excl(mi0->nrexcl, mi0->atoms.nr, mi0->interactions, &(mi0->excls));
828 gmx::mergeExclusions(&(mi0->excls), exclusionBlocks[whichmol]);
829 make_shake(mi0->interactions, &mi0->atoms, opts->nshake);
833 convert_moltype_couple(mi0, dcatt, *fudgeQQ, opts->couple_lam0,
834 opts->couple_lam1, opts->bCoupleIntra,
835 nb_funct, &(interactions[nb_funct]), wi);
837 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
838 mi0->bProcessed = TRUE;
843 fprintf(stderr, "case: %d\n", static_cast<int>(d));
844 gmx_incons("unknown directive");
853 // Check that all strings defined with -D were used when processing topology
854 std::string unusedDefineWarning = checkAndWarnForUnusedDefines(*handle);
855 if (!unusedDefineWarning.empty())
857 warning(wi, unusedDefineWarning);
860 sfree(cpp_opts_return);
867 /* List of GROMACS define names for force fields that have been
868 * parametrized using constraints involving hydrogens only.
870 * We should avoid hardcoded names, but this is hopefully only
871 * needed temparorily for discouraging use of constraints=all-bonds.
873 const std::array<std::string, 3> ffDefines = { "_FF_AMBER", "_FF_CHARMM", "_FF_OPLSAA" };
874 *ffParametrizedWithHBondConstraints = false;
875 for (const std::string& ffDefine : ffDefines)
877 if (cpp_find_define(&handle, ffDefine))
879 *ffParametrizedWithHBondConstraints = true;
883 if (cpp_find_define(&handle, "_FF_GROMOS96") != nullptr)
886 "The GROMOS force fields have been parametrized with a physically incorrect "
887 "multiple-time-stepping scheme for a twin-range cut-off. When used with "
888 "a single-range cut-off (or a correct Trotter multiple-time-stepping scheme), "
889 "physical properties, such as the density, might differ from the intended values. "
890 "Since there are researchers actively working on validating GROMOS with modern "
891 "integrators we have not yet removed the GROMOS force fields, but you should be "
892 "aware of these issues and check if molecules in your system are affected before "
894 "Further information is available at https://redmine.gromacs.org/issues/2884 , "
895 "and a longer explanation of our decision to remove physically incorrect "
897 "can be found at https://doi.org/10.26434/chemrxiv.11474583.v1 .");
902 if (opts->couple_moltype)
904 if (nmol_couple == 0)
906 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling", opts->couple_moltype);
908 fprintf(stderr, "Coupling %d copies of molecule type '%s'\n", nmol_couple, opts->couple_moltype);
911 /* this is not very clean, but fixes core dump on empty system name */
914 title = put_symtab(symtab, "");
919 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
920 warning_note(wi, warn_buf);
922 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
924 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt,
925 floating_point_arithmetic_tip);
926 warning_note(wi, warn_buf);
928 if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
931 "You are using Ewald electrostatics in a system with net charge. This can lead to "
932 "severe artifacts, such as ions moving into regions with low dielectric, due to "
933 "the uniform background charge. We suggest to neutralize your system with counter "
934 "ions, possibly in combination with a physiological salt concentration.");
935 please_cite(stdout, "Hub2014a");
940 if (*intermolecular_interactions != nullptr)
942 sfree(intermolecular_interactions->get()->atoms.atom);
948 char** do_top(bool bVerbose,
950 const char* topppfile,
954 gmx::ArrayRef<InteractionsOfType> interactions,
955 int* combination_rule,
956 double* repulsion_power,
958 PreprocessingAtomTypes* atypes,
959 std::vector<MoleculeInformation>* molinfo,
960 std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
961 const t_inputrec* ir,
962 std::vector<gmx_molblock_t>* molblock,
963 bool* ffParametrizedWithHBondConstraints,
966 /* Tmpfile might contain a long path */
981 printf("processing topology...\n");
983 title = read_topol(topfile, tmpfile, opts->define, opts->include, symtab, atypes, molinfo,
984 intermolecular_interactions, interactions, combination_rule, repulsion_power,
985 opts, fudgeQQ, molblock, ffParametrizedWithHBondConstraints,
986 ir->efep != efepNO, bZero, EEL_FULL(ir->coulombtype), wi);
988 if ((*combination_rule != eCOMB_GEOMETRIC) && (ir->vdwtype == evdwUSER))
991 "Using sigma/epsilon based combination rules with"
992 " user supplied potential function may produce unwanted"
1000 * Generate exclusion lists for QM/MM.
1002 * This routine updates the exclusion lists for QM atoms in order to include all other QM
1003 * atoms of this molecule. Moreover, this routine replaces bonds between QM atoms with
1004 * CONNBOND and, when MiMiC is not used, removes bonded interactions between QM and link atoms.
1005 * Finally, in case if MiMiC QM/MM is used - charges of QM atoms are set to 0
1007 * @param molt molecule type with QM atoms
1008 * @param grpnr group informatio
1009 * @param ir input record
1010 * @param qmmmMode QM/MM mode switch: original/MiMiC
1012 static void generate_qmexcl_moltype(gmx_moltype_t* molt, const unsigned char* grpnr, t_inputrec* ir, GmxQmmmMode qmmmMode)
1014 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1016 /* generates the exclusions between the individual QM atoms, as
1017 * these interactions should be handled by the QM subroutines and
1018 * not by the gromacs routines
1020 int qm_max = 0, qm_nr = 0, link_nr = 0, link_max = 0;
1021 int * qm_arr = nullptr, *link_arr = nullptr;
1022 bool *bQMMM, *blink;
1024 /* First we search and select the QM atoms in an qm_arr array that
1025 * we use to create the exclusions.
1027 * we take the possibility into account that a user has defined more
1028 * than one QM group:
1030 * for that we also need to do this an ugly work-about just in case
1031 * the QM group contains the entire system...
1034 /* we first search for all the QM atoms and put them in an array
1036 for (int j = 0; j < ir->opts.ngQM; j++)
1038 for (int i = 0; i < molt->atoms.nr; i++)
1040 if (qm_nr >= qm_max)
1043 srenew(qm_arr, qm_max);
1045 if ((grpnr ? grpnr[i] : 0) == j)
1047 qm_arr[qm_nr++] = i;
1048 molt->atoms.atom[i].q = 0.0;
1049 molt->atoms.atom[i].qB = 0.0;
1053 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1054 * QM/not QM. We first set all elements to false. Afterwards we use
1055 * the qm_arr to change the elements corresponding to the QM atoms
1058 snew(bQMMM, molt->atoms.nr);
1059 for (int i = 0; i < molt->atoms.nr; i++)
1063 for (int i = 0; i < qm_nr; i++)
1065 bQMMM[qm_arr[i]] = TRUE;
1068 /* We remove all bonded interactions (i.e. bonds,
1069 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1070 * are removed is as follows: if the interaction invloves 2 atoms,
1071 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1072 * it is removed if at least two of the atoms are QM atoms, if the
1073 * interaction involves 4 atoms, it is removed if there are at least
1074 * 2 QM atoms. Since this routine is called once before any forces
1075 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1076 * be rewritten at this poitn without any problem. 25-9-2002 */
1078 /* first check whether we already have CONNBONDS.
1079 * Note that if we don't, we don't add a param entry and set ftype=0,
1080 * which is ok, since CONNBONDS does not use parameters.
1082 int ftype_connbond = 0;
1083 int ind_connbond = 0;
1084 if (molt->ilist[F_CONNBONDS].size() != 0)
1086 fprintf(stderr, "nr. of CONNBONDS present already: %d\n", molt->ilist[F_CONNBONDS].size() / 3);
1087 ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1088 ind_connbond = molt->ilist[F_CONNBONDS].size();
1090 /* now we delete all bonded interactions, except the ones describing
1091 * a chemical bond. These are converted to CONNBONDS
1093 for (int ftype = 0; ftype < F_NRE; ftype++)
1095 if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_CONNBONDS)
1099 int nratoms = interaction_function[ftype].nratoms;
1101 while (j < molt->ilist[ftype].size())
1107 /* Remove an interaction between two atoms when both are
1108 * in the QM region. Note that we don't have to worry about
1109 * link atoms here, as they won't have 2-atom interactions.
1111 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1112 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1113 bexcl = (bQMMM[a1] && bQMMM[a2]);
1114 /* A chemical bond between two QM atoms will be copied to
1115 * the F_CONNBONDS list, for reasons mentioned above.
1117 if (bexcl && IS_CHEMBOND(ftype))
1119 InteractionList& ilist = molt->ilist[F_CONNBONDS];
1120 ilist.iatoms.resize(ind_connbond + 3);
1121 ilist.iatoms[ind_connbond++] = ftype_connbond;
1122 ilist.iatoms[ind_connbond++] = a1;
1123 ilist.iatoms[ind_connbond++] = a2;
1128 /* MM interactions have to be excluded if they are included
1129 * in the QM already. Because we use a link atom (H atom)
1130 * when the QM/MM boundary runs through a chemical bond, this
1131 * means that as long as one atom is MM, we still exclude,
1132 * as the interaction is included in the QM via:
1133 * QMatom1-QMatom2-QMatom-3-Linkatom.
1136 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1138 if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1144 /* MiMiC treats link atoms as quantum atoms - therefore
1145 * we do not need do additional exclusions here */
1146 if (qmmmMode == GmxQmmmMode::GMX_QMMM_MIMIC)
1148 bexcl = numQmAtoms == nratoms;
1152 bexcl = (numQmAtoms >= nratoms - 1);
1155 if (bexcl && ftype == F_SETTLE)
1158 "Can not apply QM to molecules with SETTLE, replace the moleculetype "
1159 "using QM and SETTLE by one without SETTLE");
1164 /* since the interaction involves QM atoms, these should be
1165 * removed from the MM ilist
1167 InteractionList& ilist = molt->ilist[ftype];
1168 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1170 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1172 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1176 j += nratoms + 1; /* the +1 is for the functype */
1180 /* Now, we search for atoms bonded to a QM atom because we also want
1181 * to exclude their nonbonded interactions with the QM atoms. The
1182 * reason for this is that this interaction is accounted for in the
1183 * linkatoms interaction with the QMatoms and would be counted
1186 if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
1188 for (int i = 0; i < F_NRE; i++)
1193 while (j < molt->ilist[i].size())
1195 int a1 = molt->ilist[i].iatoms[j + 1];
1196 int a2 = molt->ilist[i].iatoms[j + 2];
1197 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1199 if (link_nr >= link_max)
1202 srenew(link_arr, link_max);
1206 link_arr[link_nr++] = a2;
1210 link_arr[link_nr++] = a1;
1218 snew(blink, molt->atoms.nr);
1219 for (int i = 0; i < molt->atoms.nr; i++)
1224 if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
1226 for (int i = 0; i < link_nr; i++)
1228 blink[link_arr[i]] = TRUE;
1231 /* creating the exclusion block for the QM atoms. Each QM atom has
1232 * as excluded elements all the other QMatoms (and itself).
1235 qmexcl.nr = molt->atoms.nr;
1236 qmexcl.nra = qm_nr * (qm_nr + link_nr) + link_nr * qm_nr;
1237 snew(qmexcl.index, qmexcl.nr + 1);
1238 snew(qmexcl.a, qmexcl.nra);
1240 for (int i = 0; i < qmexcl.nr; i++)
1242 qmexcl.index[i] = j;
1245 for (int k = 0; k < qm_nr; k++)
1247 qmexcl.a[k + j] = qm_arr[k];
1249 for (int k = 0; k < link_nr; k++)
1251 qmexcl.a[qm_nr + k + j] = link_arr[k];
1253 j += (qm_nr + link_nr);
1257 for (int k = 0; k < qm_nr; k++)
1259 qmexcl.a[k + j] = qm_arr[k];
1264 qmexcl.index[qmexcl.nr] = j;
1266 /* and merging with the exclusions already present in sys.
1269 std::vector<gmx::ExclusionBlock> qmexcl2(molt->atoms.nr);
1270 gmx::blockaToExclusionBlocks(&qmexcl, qmexcl2);
1271 gmx::mergeExclusions(&(molt->excls), qmexcl2);
1273 /* Finally, we also need to get rid of the pair interactions of the
1274 * classical atom bonded to the boundary QM atoms with the QMatoms,
1275 * as this interaction is already accounted for by the QM, so also
1276 * here we run the risk of double counting! We proceed in a similar
1277 * way as we did above for the other bonded interactions: */
1278 for (int i = F_LJ14; i < F_COUL14; i++)
1280 int nratoms = interaction_function[i].nratoms;
1282 while (j < molt->ilist[i].size())
1284 int a1 = molt->ilist[i].iatoms[j + 1];
1285 int a2 = molt->ilist[i].iatoms[j + 2];
1287 ((bQMMM[a1] && bQMMM[a2]) || (blink[a1] && bQMMM[a2]) || (bQMMM[a1] && blink[a2]));
1290 /* since the interaction involves QM atoms, these should be
1291 * removed from the MM ilist
1293 InteractionList& ilist = molt->ilist[i];
1294 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1296 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1298 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1302 j += nratoms + 1; /* the +1 is for the functype */
1311 } /* generate_qmexcl */
1313 void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, warninp* wi, GmxQmmmMode qmmmMode)
1315 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1318 unsigned char* grpnr;
1319 int mol, nat_mol, nr_mol_with_qm_atoms = 0;
1320 gmx_molblock_t* molb;
1322 int index_offset = 0;
1325 grpnr = sys->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].data();
1327 for (size_t mb = 0; mb < sys->molblock.size(); mb++)
1329 molb = &sys->molblock[mb];
1330 nat_mol = sys->moltype[molb->type].atoms.nr;
1331 for (mol = 0; mol < molb->nmol; mol++)
1334 for (int i = 0; i < nat_mol; i++)
1336 if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM))
1345 nr_mol_with_qm_atoms++;
1348 /* We need to split this molblock */
1351 /* Split the molblock at this molecule */
1352 auto pos = sys->molblock.begin() + mb + 1;
1353 sys->molblock.insert(pos, sys->molblock[mb]);
1354 sys->molblock[mb].nmol = mol;
1355 sys->molblock[mb + 1].nmol -= mol;
1357 molb = &sys->molblock[mb];
1361 /* Split the molblock after this molecule */
1362 auto pos = sys->molblock.begin() + mb + 1;
1363 sys->molblock.insert(pos, sys->molblock[mb]);
1364 molb = &sys->molblock[mb];
1365 sys->molblock[mb].nmol = 1;
1366 sys->molblock[mb + 1].nmol -= 1;
1369 /* Create a copy of a moltype for a molecule
1370 * containing QM atoms and append it in the end of the list
1372 std::vector<gmx_moltype_t> temp(sys->moltype.size());
1373 for (size_t i = 0; i < sys->moltype.size(); ++i)
1375 copy_moltype(&sys->moltype[i], &temp[i]);
1377 sys->moltype.resize(sys->moltype.size() + 1);
1378 for (size_t i = 0; i < temp.size(); ++i)
1380 copy_moltype(&temp[i], &sys->moltype[i]);
1382 copy_moltype(&sys->moltype[molb->type], &sys->moltype.back());
1383 /* Copy the exclusions to a new array, since this is the only
1384 * thing that needs to be modified for QMMM.
1386 sys->moltype.back().excls = sys->moltype[molb->type].excls;
1387 /* Set the molecule type for the QMMM molblock */
1388 molb->type = sys->moltype.size() - 1;
1390 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, qmmmMode);
1396 index_offset += nat_mol;
1399 if (qmmmMode == GmxQmmmMode::GMX_QMMM_ORIGINAL && nr_mol_with_qm_atoms > 1)
1401 /* generate a warning is there are QM atoms in different
1402 * topologies. In this case it is not possible at this stage to
1403 * mutualy exclude the non-bonded interactions via the
1404 * exclusions (AFAIK). Instead, the user is advised to use the
1405 * energy group exclusions in the mdp file
1408 "\nThe QM subsystem is divided over multiple topologies. "
1409 "The mutual non-bonded interactions cannot be excluded. "
1410 "There are two ways to achieve this:\n\n"
1411 "1) merge the topologies, such that the atoms of the QM "
1412 "subsystem are all present in one single topology file. "
1413 "In this case this warning will dissappear\n\n"
1414 "2) exclude the non-bonded interactions explicitly via the "
1415 "energygrp-excl option in the mdp file. if this is the case "
1416 "this warning may be ignored"