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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
51 #include <unordered_set>
52 #include <sys/types.h>
54 #include "gromacs/fileio/gmxfio.h"
55 #include "gromacs/fileio/warninp.h"
56 #include "gromacs/gmxpreprocess/gmxcpp.h"
57 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
58 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
59 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
60 #include "gromacs/gmxpreprocess/grompp-impl.h"
61 #include "gromacs/gmxpreprocess/readir.h"
62 #include "gromacs/gmxpreprocess/topdirs.h"
63 #include "gromacs/gmxpreprocess/toppush.h"
64 #include "gromacs/gmxpreprocess/topshake.h"
65 #include "gromacs/gmxpreprocess/toputil.h"
66 #include "gromacs/gmxpreprocess/vsite_parm.h"
67 #include "gromacs/math/units.h"
68 #include "gromacs/math/utilities.h"
69 #include "gromacs/mdtypes/inputrec.h"
70 #include "gromacs/mdtypes/md_enums.h"
71 #include "gromacs/pbcutil/pbc.h"
72 #include "gromacs/topology/block.h"
73 #include "gromacs/topology/exclusionblocks.h"
74 #include "gromacs/topology/ifunc.h"
75 #include "gromacs/topology/symtab.h"
76 #include "gromacs/topology/topology.h"
77 #include "gromacs/utility/cstringutil.h"
78 #include "gromacs/utility/fatalerror.h"
79 #include "gromacs/utility/futil.h"
80 #include "gromacs/utility/gmxassert.h"
81 #include "gromacs/utility/pleasecite.h"
82 #include "gromacs/utility/smalloc.h"
84 #define OPENDIR '[' /* starting sign for directive */
85 #define CLOSEDIR ']' /* ending sign for directive */
87 static void gen_pairs(t_params *nbs, t_params *pairs, real fudge, int comb)
89 int i, j, ntp, nrfp, nrfpA, nrfpB, nnn;
92 nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
93 GMX_ASSERT(nnn * nnn == ntp, "Number of pairs of generated non-bonded parameters should be a perfect square");
95 nrfpA = interaction_function[F_LJ14].nrfpA;
96 nrfpB = interaction_function[F_LJ14].nrfpB;
99 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
101 gmx_incons("Number of force parameters in gen_pairs wrong");
104 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
105 snew(pairs->param, pairs->nr);
106 for (i = 0; (i < ntp); i++)
109 pairs->param[i].a[0] = i / nnn;
110 pairs->param[i].a[1] = i % nnn;
111 /* Copy normal and FEP parameters and multiply by fudge factor */
115 for (j = 0; (j < nrfp); j++)
117 /* If we are using sigma/epsilon values, only the epsilon values
118 * should be scaled, but not sigma.
119 * The sigma values have even indices 0,2, etc.
121 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2 == 0))
130 pairs->param[i].c[j] = scaling*nbs->param[i].c[j];
131 /* NOTE: this should be clear to the compiler, but some gcc 5.2 versions
132 * issue false positive warnings for the pairs->param.c[] indexing below.
134 assert(2*nrfp <= MAXFORCEPARAM);
135 pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j];
140 double check_mol(const gmx_mtop_t *mtop, warninp *wi)
147 /* Check mass and charge */
150 for (const gmx_molblock_t &molb : mtop->molblock)
152 const t_atoms *atoms = &mtop->moltype[molb.type].atoms;
153 for (i = 0; (i < atoms->nr); i++)
155 q += molb.nmol*atoms->atom[i].q;
156 m = atoms->atom[i].m;
157 mB = atoms->atom[i].mB;
158 pt = atoms->atom[i].ptype;
159 /* If the particle is an atom or a nucleus it must have a mass,
160 * else, if it is a shell, a vsite or a bondshell it can have mass zero
162 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
164 ri = atoms->atom[i].resind;
165 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
166 *(atoms->atomname[i]),
167 *(atoms->resinfo[ri].name),
168 atoms->resinfo[ri].nr,
170 warning_error(wi, buf);
173 if (((m != 0) || (mB != 0)) && (pt == eptVSite))
175 ri = atoms->atom[i].resind;
176 sprintf(buf, "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state B)\n"
177 " Check your topology.\n",
178 *(atoms->atomname[i]),
179 *(atoms->resinfo[ri].name),
180 atoms->resinfo[ri].nr,
182 warning_error(wi, buf);
183 /* The following statements make LINCS break! */
184 /* atoms->atom[i].m=0; */
191 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
193 * The results of this routine are only used for checking and for
194 * printing warning messages. Thus we can assume that charges of molecules
195 * should be integer. If the user wanted non-integer molecular charge,
196 * an undesired warning is printed and the user should use grompp -maxwarn 1.
198 * \param qMol The total, unrounded, charge of the molecule
199 * \param sumAbsQ The sum of absolute values of the charges, used for determining the tolerance for the rounding.
201 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
203 /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
204 * of the charges for ascii float truncation in the topology files.
205 * Although the summation here uses double precision, the charges
206 * are read and stored in single precision when real=float. This can
207 * lead to rounding errors of half the least significant bit.
208 * Note that, unfortunately, we can not assume addition of random
209 * rounding errors. It is not entirely unlikely that many charges
210 * have a near half-bit rounding error with the same sign.
212 double tolAbs = 1e-6;
213 double tol = std::max(tolAbs, 0.5*GMX_REAL_EPS*sumAbsQ);
214 double qRound = std::round(qMol);
215 if (std::abs(qMol - qRound) <= tol)
225 static void sum_q(const t_atoms *atoms, int numMols,
226 double *qTotA, double *qTotB)
233 for (int i = 0; i < atoms->nr; i++)
235 qmolA += atoms->atom[i].q;
236 qmolB += atoms->atom[i].qB;
237 sumAbsQA += std::abs(atoms->atom[i].q);
238 sumAbsQB += std::abs(atoms->atom[i].qB);
241 *qTotA += numMols*roundedMoleculeCharge(qmolA, sumAbsQA);
242 *qTotB += numMols*roundedMoleculeCharge(qmolB, sumAbsQB);
245 static void get_nbparm(char *nb_str, char *comb_str, int *nb, int *comb,
249 char warn_buf[STRLEN];
252 for (i = 1; (i < eNBF_NR); i++)
254 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
261 *nb = strtol(nb_str, nullptr, 10);
263 if ((*nb < 1) || (*nb >= eNBF_NR))
265 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s",
266 nb_str, enbf_names[1]);
267 warning_error(wi, warn_buf);
271 for (i = 1; (i < eCOMB_NR); i++)
273 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
280 *comb = strtol(comb_str, nullptr, 10);
282 if ((*comb < 1) || (*comb >= eCOMB_NR))
284 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s",
285 comb_str, ecomb_names[1]);
286 warning_error(wi, warn_buf);
291 static char ** cpp_opts(const char *define, const char *include,
296 const char *cppadds[2];
297 char **cppopts = nullptr;
298 const char *option[2] = { "-D", "-I" };
299 const char *nopt[2] = { "define", "include" };
303 char warn_buf[STRLEN];
306 cppadds[1] = include;
307 for (n = 0; (n < 2); n++)
314 while ((*ptr != '\0') && isspace(*ptr))
319 while ((*rptr != '\0') && !isspace(*rptr))
327 strncpy(buf, ptr, len);
328 if (strstr(ptr, option[n]) != ptr)
330 set_warning_line(wi, "mdp file", -1);
331 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
332 warning(wi, warn_buf);
336 srenew(cppopts, ++ncppopts);
337 cppopts[ncppopts-1] = gmx_strdup(buf);
345 srenew(cppopts, ++ncppopts);
346 cppopts[ncppopts-1] = nullptr;
352 static void make_atoms_sys(const std::vector<gmx_molblock_t> &molblock,
353 const t_molinfo *molinfo,
357 atoms->atom = nullptr;
359 for (const gmx_molblock_t &molb : molblock)
361 const t_atoms &mol_atoms = molinfo[molb.type].atoms;
363 srenew(atoms->atom, atoms->nr + molb.nmol*mol_atoms.nr);
365 for (int m = 0; m < molb.nmol; m++)
367 for (int a = 0; a < mol_atoms.nr; a++)
369 atoms->atom[atoms->nr++] = mol_atoms.atom[a];
376 static char **read_topol(const char *infile, const char *outfile,
377 const char *define, const char *include,
382 t_molinfo **intermolecular_interactions,
384 int *combination_rule,
388 std::vector<gmx_molblock_t> *molblock,
389 bool *ffParametrizedWithHBondConstraints,
392 bool usingFullRangeElectrostatics,
397 char *pline = nullptr, **title = nullptr;
398 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
400 char *dirstr, *dummy2;
401 int nrcopies, nmol, nscan, ncombs, ncopy;
402 double fLJ, fQQ, fPOW;
403 t_molinfo *mi0 = nullptr;
406 t_nbparam **nbparam, **pair;
407 gmx::ExclusionBlocks *exclusionBlocks;
408 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
409 bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
410 double qt = 0, qBt = 0; /* total charge */
411 gpp_bond_atomtype *batype;
413 int dcatt = -1, nmol_couple;
414 /* File handling variables */
418 char *tmp_line = nullptr;
419 char warn_buf[STRLEN];
420 const char *floating_point_arithmetic_tip =
421 "Total charge should normally be an integer. See\n"
422 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
423 "for discussion on how close it should be to an integer.\n";
424 /* We need to open the output file before opening the input file,
425 * because cpp_open_file can change the current working directory.
429 out = gmx_fio_fopen(outfile, "w");
436 /* open input file */
437 auto cpp_opts_return = cpp_opts(define, include, wi);
438 status = cpp_open_file(infile, &handle, cpp_opts_return);
441 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
444 /* some local variables */
445 DS_Init(&DS); /* directive stack */
446 nmol = 0; /* no molecules yet... */
447 d = Directive::d_invalid; /* first thing should be a directive */
448 nbparam = nullptr; /* The temporary non-bonded matrix */
449 pair = nullptr; /* The temporary pair interaction matrix */
450 exclusionBlocks = nullptr; /* the extra exclusions */
453 *reppow = 12.0; /* Default value for repulsion power */
455 *intermolecular_interactions = nullptr;
457 /* Init the number of CMAP torsion angles and grid spacing */
458 plist[F_CMAP].grid_spacing = 0;
459 plist[F_CMAP].nc = 0;
461 bWarn_copy_A_B = bFEP;
463 batype = init_bond_atomtype();
464 /* parse the actual file */
465 bReadDefaults = FALSE;
467 bReadMolType = FALSE;
472 status = cpp_read_line(&handle, STRLEN, line);
473 done = (status == eCPP_EOF);
476 if (status != eCPP_OK)
478 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
482 fprintf(out, "%s\n", line);
485 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
487 pline = gmx_strdup(line);
489 /* Strip trailing '\' from pline, if it exists */
491 if ((sl > 0) && (pline[sl-1] == CONTINUE))
496 /* build one long line from several fragments - necessary for CMAP */
497 while (continuing(line))
499 status = cpp_read_line(&handle, STRLEN, line);
500 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
502 /* Since we depend on the '\' being present to continue to read, we copy line
503 * to a tmp string, strip the '\' from that string, and cat it to pline
505 tmp_line = gmx_strdup(line);
507 sl = strlen(tmp_line);
508 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE))
510 tmp_line[sl-1] = ' ';
513 done = (status == eCPP_EOF);
516 if (status != eCPP_OK)
518 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
522 fprintf(out, "%s\n", line);
526 srenew(pline, strlen(pline)+strlen(tmp_line)+1);
527 strcat(pline, tmp_line);
531 /* skip trailing and leading spaces and comment text */
532 strip_comment (pline);
535 /* if there is something left... */
536 if (static_cast<int>(strlen(pline)) > 0)
538 if (pline[0] == OPENDIR)
540 /* A directive on this line: copy the directive
541 * without the brackets into dirstr, then
542 * skip spaces and tabs on either side of directive
544 dirstr = gmx_strdup((pline+1));
545 if ((dummy2 = strchr (dirstr, CLOSEDIR)) != nullptr)
551 if ((newd = str2dir(dirstr)) == Directive::d_invalid)
553 sprintf(errbuf, "Invalid directive %s", dirstr);
554 warning_error(wi, errbuf);
558 /* Directive found */
559 if (DS_Check_Order (DS, newd))
566 /* we should print here which directives should have
567 been present, and which actually are */
568 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
569 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
570 /* d = Directive::d_invalid; */
573 if (d == Directive::d_intermolecular_interactions)
575 if (*intermolecular_interactions == nullptr)
577 /* We (mis)use the moleculetype processing
578 * to process the intermolecular interactions
579 * by making a "molecule" of the size of the system.
581 snew(*intermolecular_interactions, 1);
582 init_molinfo(*intermolecular_interactions);
583 mi0 = *intermolecular_interactions;
584 make_atoms_sys(*molblock, *molinfo,
591 else if (d != Directive::d_invalid)
593 /* Not a directive, just a plain string
594 * use a gigantic switch to decode,
595 * if there is a valid directive!
599 case Directive::d_defaults:
602 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
603 cpp_error(&handle, eCPP_SYNTAX));
605 bReadDefaults = TRUE;
606 nscan = sscanf(pline, "%s%s%s%lf%lf%lf",
607 nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
618 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
621 bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
622 if (nb_funct != eNBF_LJ && bGenPairs)
624 gmx_fatal(FARGS, "Generating pair parameters is only supported with LJ non-bonded interactions");
640 nb_funct = ifunc_index(Directive::d_nonbond_params, nb_funct);
643 case Directive::d_atomtypes:
644 push_at(symtab, atype, batype, pline, nb_funct,
645 &nbparam, bGenPairs ? &pair : nullptr, wi);
648 case Directive::d_bondtypes:
649 push_bt(d, plist, 2, nullptr, batype, pline, wi);
651 case Directive::d_constrainttypes:
652 push_bt(d, plist, 2, nullptr, batype, pline, wi);
654 case Directive::d_pairtypes:
657 push_nbt(d, pair, atype, pline, F_LJ14, wi);
661 push_bt(d, plist, 2, atype, nullptr, pline, wi);
664 case Directive::d_angletypes:
665 push_bt(d, plist, 3, nullptr, batype, pline, wi);
667 case Directive::d_dihedraltypes:
668 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
669 push_dihedraltype(d, plist, batype, pline, wi);
672 case Directive::d_nonbond_params:
673 push_nbt(d, nbparam, atype, pline, nb_funct, wi);
676 case Directive::d_implicit_genborn_params:
677 // Skip this line, so old topologies with
678 // GB parameters can be read.
681 case Directive::d_implicit_surface_params:
682 // Skip this line, so that any topologies
683 // with surface parameters can be read
684 // (even though these were never formally
688 case Directive::d_cmaptypes:
689 push_cmaptype(d, plist, 5, atype, batype, pline, wi);
692 case Directive::d_moleculetype:
697 if (opts->couple_moltype != nullptr &&
698 (opts->couple_lam0 == ecouplamNONE ||
699 opts->couple_lam0 == ecouplamQ ||
700 opts->couple_lam1 == ecouplamNONE ||
701 opts->couple_lam1 == ecouplamQ))
703 dcatt = add_atomtype_decoupled(symtab, atype,
704 &nbparam, bGenPairs ? &pair : nullptr);
706 ntype = get_atomtype_ntypes(atype);
707 ncombs = (ntype*(ntype+1))/2;
708 generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi);
709 ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
711 fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
712 free_nbparam(nbparam, ntype);
715 gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule);
716 ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
718 fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
719 free_nbparam(pair, ntype);
721 /* Copy GBSA parameters to atomtype array? */
726 push_molt(symtab, &nmol, molinfo, pline, wi);
727 srenew(exclusionBlocks, nmol);
728 exclusionBlocks[nmol-1].nr = 0;
729 mi0 = &((*molinfo)[nmol-1]);
730 mi0->atoms.haveMass = TRUE;
731 mi0->atoms.haveCharge = TRUE;
732 mi0->atoms.haveType = TRUE;
733 mi0->atoms.haveBState = TRUE;
734 mi0->atoms.havePdbInfo = FALSE;
737 case Directive::d_atoms:
738 push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi);
741 case Directive::d_pairs:
742 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
743 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
745 case Directive::d_pairs_nb:
746 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
747 FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
750 case Directive::d_vsites2:
751 case Directive::d_vsites3:
752 case Directive::d_vsites4:
753 case Directive::d_bonds:
754 case Directive::d_angles:
755 case Directive::d_constraints:
756 case Directive::d_settles:
757 case Directive::d_position_restraints:
758 case Directive::d_angle_restraints:
759 case Directive::d_angle_restraints_z:
760 case Directive::d_distance_restraints:
761 case Directive::d_orientation_restraints:
762 case Directive::d_dihedral_restraints:
763 case Directive::d_dihedrals:
764 case Directive::d_polarization:
765 case Directive::d_water_polarization:
766 case Directive::d_thole_polarization:
767 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE,
768 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
770 case Directive::d_cmap:
771 push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi);
774 case Directive::d_vsitesn:
775 push_vsitesn(d, mi0->plist, &(mi0->atoms), pline, wi);
777 case Directive::d_exclusions:
778 GMX_ASSERT(exclusionBlocks, "exclusionBlocks must always be allocated so exclusions can be processed");
779 if (!exclusionBlocks[nmol-1].nr)
781 initExclusionBlocks(&(exclusionBlocks[nmol-1]), mi0->atoms.nr);
783 push_excl(pline, &(exclusionBlocks[nmol-1]), wi);
785 case Directive::d_system:
787 title = put_symtab(symtab, pline);
789 case Directive::d_molecules:
794 push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
795 mi0 = &((*molinfo)[whichmol]);
796 molblock->resize(molblock->size() + 1);
797 molblock->back().type = whichmol;
798 molblock->back().nmol = nrcopies;
800 bCouple = (opts->couple_moltype != nullptr &&
801 (gmx_strcasecmp("system", opts->couple_moltype) == 0 ||
802 strcmp(*(mi0->name), opts->couple_moltype) == 0));
805 nmol_couple += nrcopies;
808 if (mi0->atoms.nr == 0)
810 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms",
814 "Excluding %d bonded neighbours molecule type '%s'\n",
815 mi0->nrexcl, *mi0->name);
816 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
817 if (!mi0->bProcessed)
820 generate_excl(mi0->nrexcl,
825 gmx::mergeExclusions(&(mi0->excls), &(exclusionBlocks[whichmol]));
826 gmx::doneExclusionBlocks(&(exclusionBlocks[whichmol]));
827 make_shake(mi0->plist, &mi0->atoms, opts->nshake);
835 convert_moltype_couple(mi0, dcatt, *fudgeQQ,
836 opts->couple_lam0, opts->couple_lam1,
838 nb_funct, &(plist[nb_funct]), wi);
840 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
841 mi0->bProcessed = TRUE;
846 fprintf (stderr, "case: %d\n", static_cast<int>(d));
847 gmx_incons("unknown directive");
857 // Check that all strings defined with -D were used when processing topology
858 std::string unusedDefineWarning = checkAndWarnForUnusedDefines(*handle);
859 if (!unusedDefineWarning.empty())
861 warning(wi, unusedDefineWarning);
864 sfree(cpp_opts_return);
871 /* List of GROMACS define names for force fields that have been
872 * parametrized using constraints involving hydrogens only.
874 * We should avoid hardcoded names, but this is hopefully only
875 * needed temparorily for discouraging use of constraints=all-bonds.
877 const std::array<std::string, 3> ffDefines = {
882 *ffParametrizedWithHBondConstraints = false;
883 for (const std::string &ffDefine : ffDefines)
885 if (cpp_find_define(&handle, ffDefine))
887 *ffParametrizedWithHBondConstraints = true;
893 if (opts->couple_moltype)
895 if (nmol_couple == 0)
897 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling",
898 opts->couple_moltype);
900 fprintf(stderr, "Coupling %d copies of molecule type '%s'\n",
901 nmol_couple, opts->couple_moltype);
904 /* this is not very clean, but fixes core dump on empty system name */
907 title = put_symtab(symtab, "");
912 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
913 warning_note(wi, warn_buf);
915 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
917 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
918 warning_note(wi, warn_buf);
920 if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
922 warning(wi, "You are using Ewald electrostatics in a system with net charge. This can lead to severe artifacts, such as ions moving into regions with low dielectric, due to the uniform background charge. We suggest to neutralize your system with counter ions, possibly in combination with a physiological salt concentration.");
923 please_cite(stdout, "Hub2014a");
927 for (i = 0; i < nmol; i++)
929 gmx::doneExclusionBlocks(&(exclusionBlocks[i]));
931 free(exclusionBlocks);
933 done_bond_atomtype(&batype);
935 if (*intermolecular_interactions != nullptr)
937 sfree(mi0->atoms.atom);
945 char **do_top(bool bVerbose,
947 const char *topppfile,
952 int *combination_rule,
953 double *repulsion_power,
958 t_molinfo **intermolecular_interactions,
959 const t_inputrec *ir,
960 std::vector<gmx_molblock_t> *molblock,
961 bool *ffParametrizedWithHBondConstraints,
964 /* Tmpfile might contain a long path */
979 printf("processing topology...\n");
981 title = read_topol(topfile, tmpfile, opts->define, opts->include,
983 nrmols, molinfo, intermolecular_interactions,
984 plist, combination_rule, repulsion_power,
985 opts, fudgeQQ, molblock,
986 ffParametrizedWithHBondConstraints,
987 ir->efep != efepNO, bZero,
988 EEL_FULL(ir->coulombtype), wi);
990 if ((*combination_rule != eCOMB_GEOMETRIC) &&
991 (ir->vdwtype == evdwUSER))
993 warning(wi, "Using sigma/epsilon based combination rules with"
994 " user supplied potential function may produce unwanted"
1002 * Generate exclusion lists for QM/MM.
1004 * This routine updates the exclusion lists for QM atoms in order to include all other QM
1005 * atoms of this molecule. Moreover, this routine replaces bonds between QM atoms with
1006 * CONNBOND and, when MiMiC is not used, removes bonded interactions between QM and link atoms.
1007 * Finally, in case if MiMiC QM/MM is used - charges of QM atoms are set to 0
1009 * @param molt molecule type with QM atoms
1010 * @param grpnr group informatio
1011 * @param ir input record
1012 * @param qmmmMode QM/MM mode switch: original/MiMiC
1014 static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *grpnr,
1015 t_inputrec *ir, GmxQmmmMode qmmmMode)
1017 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1019 /* generates the exclusions between the individual QM atoms, as
1020 * these interactions should be handled by the QM subroutines and
1021 * not by the gromacs routines
1023 int qm_max = 0, qm_nr = 0, link_nr = 0, link_max = 0;
1024 int *qm_arr = nullptr, *link_arr = nullptr;
1025 bool *bQMMM, *blink;
1027 /* First we search and select the QM atoms in an qm_arr array that
1028 * we use to create the exclusions.
1030 * we take the possibility into account that a user has defined more
1031 * than one QM group:
1033 * for that we also need to do this an ugly work-about just in case
1034 * the QM group contains the entire system...
1037 /* we first search for all the QM atoms and put them in an array
1039 for (int j = 0; j < ir->opts.ngQM; j++)
1041 for (int i = 0; i < molt->atoms.nr; i++)
1043 if (qm_nr >= qm_max)
1046 srenew(qm_arr, qm_max);
1048 if ((grpnr ? grpnr[i] : 0) == j)
1050 qm_arr[qm_nr++] = i;
1051 molt->atoms.atom[i].q = 0.0;
1052 molt->atoms.atom[i].qB = 0.0;
1056 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1057 * QM/not QM. We first set all elements to false. Afterwards we use
1058 * the qm_arr to change the elements corresponding to the QM atoms
1061 snew(bQMMM, molt->atoms.nr);
1062 for (int i = 0; i < molt->atoms.nr; i++)
1066 for (int i = 0; i < qm_nr; i++)
1068 bQMMM[qm_arr[i]] = TRUE;
1071 /* We remove all bonded interactions (i.e. bonds,
1072 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1073 * are removed is as follows: if the interaction invloves 2 atoms,
1074 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1075 * it is removed if at least two of the atoms are QM atoms, if the
1076 * interaction involves 4 atoms, it is removed if there are at least
1077 * 2 QM atoms. Since this routine is called once before any forces
1078 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1079 * be rewritten at this poitn without any problem. 25-9-2002 */
1081 /* first check whether we already have CONNBONDS.
1082 * Note that if we don't, we don't add a param entry and set ftype=0,
1083 * which is ok, since CONNBONDS does not use parameters.
1085 int ftype_connbond = 0;
1086 int ind_connbond = 0;
1087 if (molt->ilist[F_CONNBONDS].size() != 0)
1089 fprintf(stderr, "nr. of CONNBONDS present already: %d\n",
1090 molt->ilist[F_CONNBONDS].size()/3);
1091 ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1092 ind_connbond = molt->ilist[F_CONNBONDS].size();
1094 /* now we delete all bonded interactions, except the ones describing
1095 * a chemical bond. These are converted to CONNBONDS
1097 for (int ftype = 0; ftype < F_NRE; ftype++)
1099 if (!(interaction_function[ftype].flags & IF_BOND) ||
1100 ftype == F_CONNBONDS)
1104 int nratoms = interaction_function[ftype].nratoms;
1106 while (j < molt->ilist[ftype].size())
1112 /* Remove an interaction between two atoms when both are
1113 * in the QM region. Note that we don't have to worry about
1114 * link atoms here, as they won't have 2-atom interactions.
1116 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1117 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1118 bexcl = (bQMMM[a1] && bQMMM[a2]);
1119 /* A chemical bond between two QM atoms will be copied to
1120 * the F_CONNBONDS list, for reasons mentioned above.
1122 if (bexcl && IS_CHEMBOND(ftype))
1124 InteractionList &ilist = molt->ilist[F_CONNBONDS];
1125 ilist.iatoms.resize(ind_connbond + 3);
1126 ilist.iatoms[ind_connbond++] = ftype_connbond;
1127 ilist.iatoms[ind_connbond++] = a1;
1128 ilist.iatoms[ind_connbond++] = a2;
1133 /* MM interactions have to be excluded if they are included
1134 * in the QM already. Because we use a link atom (H atom)
1135 * when the QM/MM boundary runs through a chemical bond, this
1136 * means that as long as one atom is MM, we still exclude,
1137 * as the interaction is included in the QM via:
1138 * QMatom1-QMatom2-QMatom-3-Linkatom.
1141 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1143 if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1149 /* MiMiC treats link atoms as quantum atoms - therefore
1150 * we do not need do additional exclusions here */
1151 if (qmmmMode == GmxQmmmMode::GMX_QMMM_MIMIC)
1153 bexcl = numQmAtoms == nratoms;
1157 bexcl = (numQmAtoms >= nratoms - 1);
1160 if (bexcl && ftype == F_SETTLE)
1162 gmx_fatal(FARGS, "Can not apply QM to molecules with SETTLE, replace the moleculetype using QM and SETTLE by one without SETTLE");
1167 /* since the interaction involves QM atoms, these should be
1168 * removed from the MM ilist
1170 InteractionList &ilist = molt->ilist[ftype];
1171 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1173 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1175 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1179 j += nratoms+1; /* the +1 is for the functype */
1183 /* Now, we search for atoms bonded to a QM atom because we also want
1184 * to exclude their nonbonded interactions with the QM atoms. The
1185 * reason for this is that this interaction is accounted for in the
1186 * linkatoms interaction with the QMatoms and would be counted
1189 if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
1191 for (int i = 0; i < F_NRE; i++)
1196 while (j < molt->ilist[i].size())
1198 int a1 = molt->ilist[i].iatoms[j + 1];
1199 int a2 = molt->ilist[i].iatoms[j + 2];
1200 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1202 if (link_nr >= link_max)
1205 srenew(link_arr, link_max);
1209 link_arr[link_nr++] = a2;
1213 link_arr[link_nr++] = a1;
1221 snew(blink, molt->atoms.nr);
1222 for (int i = 0; i < molt->atoms.nr; i++)
1227 if (qmmmMode != GmxQmmmMode::GMX_QMMM_MIMIC)
1229 for (int i = 0; i < link_nr; i++)
1231 blink[link_arr[i]] = TRUE;
1234 /* creating the exclusion block for the QM atoms. Each QM atom has
1235 * as excluded elements all the other QMatoms (and itself).
1238 qmexcl.nr = molt->atoms.nr;
1239 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1240 snew(qmexcl.index, qmexcl.nr+1);
1241 snew(qmexcl.a, qmexcl.nra);
1243 for (int i = 0; i < qmexcl.nr; i++)
1245 qmexcl.index[i] = j;
1248 for (int k = 0; k < qm_nr; k++)
1250 qmexcl.a[k+j] = qm_arr[k];
1252 for (int k = 0; k < link_nr; k++)
1254 qmexcl.a[qm_nr+k+j] = link_arr[k];
1256 j += (qm_nr+link_nr);
1260 for (int k = 0; k < qm_nr; k++)
1262 qmexcl.a[k+j] = qm_arr[k];
1267 qmexcl.index[qmexcl.nr] = j;
1269 /* and merging with the exclusions already present in sys.
1272 gmx::ExclusionBlocks qmexcl2;
1273 initExclusionBlocks(&qmexcl2, molt->atoms.nr);
1274 gmx::blockaToExclusionBlocks(&qmexcl, &qmexcl2);
1275 gmx::mergeExclusions(&(molt->excls), &qmexcl2);
1276 gmx::doneExclusionBlocks(&qmexcl2);
1278 /* Finally, we also need to get rid of the pair interactions of the
1279 * classical atom bonded to the boundary QM atoms with the QMatoms,
1280 * as this interaction is already accounted for by the QM, so also
1281 * here we run the risk of double counting! We proceed in a similar
1282 * way as we did above for the other bonded interactions: */
1283 for (int i = F_LJ14; i < F_COUL14; i++)
1285 int nratoms = interaction_function[i].nratoms;
1287 while (j < molt->ilist[i].size())
1289 int a1 = molt->ilist[i].iatoms[j+1];
1290 int a2 = molt->ilist[i].iatoms[j+2];
1291 bool bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1292 (blink[a1] && bQMMM[a2]) ||
1293 (bQMMM[a1] && blink[a2]));
1296 /* since the interaction involves QM atoms, these should be
1297 * removed from the MM ilist
1299 InteractionList &ilist = molt->ilist[i];
1300 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1302 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1304 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1308 j += nratoms+1; /* the +1 is for the functype */
1317 } /* generate_qmexcl */
1319 void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp *wi, GmxQmmmMode qmmmMode)
1321 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1324 unsigned char *grpnr;
1325 int mol, nat_mol, nr_mol_with_qm_atoms = 0;
1326 gmx_molblock_t *molb;
1328 int index_offset = 0;
1331 grpnr = sys->groups.grpnr[egcQMMM];
1333 for (size_t mb = 0; mb < sys->molblock.size(); mb++)
1335 molb = &sys->molblock[mb];
1336 nat_mol = sys->moltype[molb->type].atoms.nr;
1337 for (mol = 0; mol < molb->nmol; mol++)
1340 for (int i = 0; i < nat_mol; i++)
1342 if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM))
1351 nr_mol_with_qm_atoms++;
1354 /* We need to split this molblock */
1357 /* Split the molblock at this molecule */
1358 auto pos = sys->molblock.begin() + mb + 1;
1359 sys->molblock.insert(pos, sys->molblock[mb]);
1360 sys->molblock[mb ].nmol = mol;
1361 sys->molblock[mb+1].nmol -= mol;
1363 molb = &sys->molblock[mb];
1367 /* Split the molblock after this molecule */
1368 auto pos = sys->molblock.begin() + mb + 1;
1369 sys->molblock.insert(pos, sys->molblock[mb]);
1370 molb = &sys->molblock[mb];
1371 sys->molblock[mb ].nmol = 1;
1372 sys->molblock[mb+1].nmol -= 1;
1375 /* Create a copy of a moltype for a molecule
1376 * containing QM atoms and append it in the end of the list
1378 std::vector<gmx_moltype_t> temp(sys->moltype.size());
1379 for (size_t i = 0; i < sys->moltype.size(); ++i)
1381 copy_moltype(&sys->moltype[i], &temp[i]);
1383 sys->moltype.resize(sys->moltype.size() + 1);
1384 for (size_t i = 0; i < temp.size(); ++i)
1386 copy_moltype(&temp[i], &sys->moltype[i]);
1388 copy_moltype(&sys->moltype[molb->type], &sys->moltype.back());
1389 /* Copy the exclusions to a new array, since this is the only
1390 * thing that needs to be modified for QMMM.
1392 copy_blocka(&sys->moltype[molb->type].excls,
1393 &sys->moltype.back().excls);
1394 /* Set the molecule type for the QMMM molblock */
1395 molb->type = sys->moltype.size() - 1;
1397 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, qmmmMode);
1403 index_offset += nat_mol;
1406 if (qmmmMode == GmxQmmmMode::GMX_QMMM_ORIGINAL &&
1407 nr_mol_with_qm_atoms > 1)
1409 /* generate a warning is there are QM atoms in different
1410 * topologies. In this case it is not possible at this stage to
1411 * mutualy exclude the non-bonded interactions via the
1412 * exclusions (AFAIK). Instead, the user is advised to use the
1413 * energy group exclusions in the mdp file
1416 "\nThe QM subsystem is divided over multiple topologies. "
1417 "The mutual non-bonded interactions cannot be excluded. "
1418 "There are two ways to achieve this:\n\n"
1419 "1) merge the topologies, such that the atoms of the QM "
1420 "subsystem are all present in one single topology file. "
1421 "In this case this warning will dissappear\n\n"
1422 "2) exclude the non-bonded interactions explicitly via the "
1423 "energygrp-excl option in the mdp file. if this is the case "
1424 "this warning may be ignored"