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, 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 <sys/types.h>
53 #include "gromacs/fileio/gmxfio.h"
54 #include "gromacs/fileio/warninp.h"
55 #include "gromacs/gmxpreprocess/gmxcpp.h"
56 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
57 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
58 #include "gromacs/gmxpreprocess/grompp-impl.h"
59 #include "gromacs/gmxpreprocess/readir.h"
60 #include "gromacs/gmxpreprocess/topdirs.h"
61 #include "gromacs/gmxpreprocess/toppush.h"
62 #include "gromacs/gmxpreprocess/topshake.h"
63 #include "gromacs/gmxpreprocess/toputil.h"
64 #include "gromacs/gmxpreprocess/vsite_parm.h"
65 #include "gromacs/math/units.h"
66 #include "gromacs/math/utilities.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/pbcutil/pbc.h"
70 #include "gromacs/topology/block.h"
71 #include "gromacs/topology/ifunc.h"
72 #include "gromacs/topology/symtab.h"
73 #include "gromacs/topology/topology.h"
74 #include "gromacs/utility/cstringutil.h"
75 #include "gromacs/utility/fatalerror.h"
76 #include "gromacs/utility/futil.h"
77 #include "gromacs/utility/gmxassert.h"
78 #include "gromacs/utility/pleasecite.h"
79 #include "gromacs/utility/smalloc.h"
81 #define OPENDIR '[' /* starting sign for directive */
82 #define CLOSEDIR ']' /* ending sign for directive */
84 static void free_nbparam(t_nbparam **param, int nr)
89 for (i = 0; i < nr; i++)
97 static int copy_nbparams(t_nbparam **param, int ftype, t_params *plist, int nr)
105 for (i = 0; i < nr; i++)
107 for (j = 0; j <= i; j++)
110 if (param[i][j].bSet)
112 for (f = 0; f < nrfp; f++)
114 plist->param[nr*i+j].c[f] = param[i][j].c[f];
115 plist->param[nr*j+i].c[f] = param[i][j].c[f];
125 static void gen_pairs(t_params *nbs, t_params *pairs, real fudge, int comb)
127 int i, j, ntp, nrfp, nrfpA, nrfpB, nnn;
130 nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
131 GMX_ASSERT(nnn * nnn == ntp, "Number of pairs of generated non-bonded parameters should be a perfect square");
133 nrfpA = interaction_function[F_LJ14].nrfpA;
134 nrfpB = interaction_function[F_LJ14].nrfpB;
137 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
139 gmx_incons("Number of force parameters in gen_pairs wrong");
142 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
143 snew(pairs->param, pairs->nr);
144 for (i = 0; (i < ntp); i++)
147 pairs->param[i].a[0] = i / nnn;
148 pairs->param[i].a[1] = i % nnn;
149 /* Copy normal and FEP parameters and multiply by fudge factor */
153 for (j = 0; (j < nrfp); j++)
155 /* If we are using sigma/epsilon values, only the epsilon values
156 * should be scaled, but not sigma.
157 * The sigma values have even indices 0,2, etc.
159 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2 == 0))
168 pairs->param[i].c[j] = scaling*nbs->param[i].c[j];
169 /* NOTE: this should be cleat to the compiler, but some gcc 5.2 versions
170 * issue false positive warnings for the pairs->param.c[] indexing below.
172 assert(2*nrfp <= MAXFORCEPARAM);
173 pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j];
178 double check_mol(const gmx_mtop_t *mtop, warninp_t wi)
185 /* Check mass and charge */
188 for (const gmx_molblock_t &molb : mtop->molblock)
190 const t_atoms *atoms = &mtop->moltype[molb.type].atoms;
191 for (i = 0; (i < atoms->nr); i++)
193 q += molb.nmol*atoms->atom[i].q;
194 m = atoms->atom[i].m;
195 mB = atoms->atom[i].mB;
196 pt = atoms->atom[i].ptype;
197 /* If the particle is an atom or a nucleus it must have a mass,
198 * else, if it is a shell, a vsite or a bondshell it can have mass zero
200 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
202 ri = atoms->atom[i].resind;
203 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
204 *(atoms->atomname[i]),
205 *(atoms->resinfo[ri].name),
206 atoms->resinfo[ri].nr,
208 warning_error(wi, buf);
211 if (((m != 0) || (mB != 0)) && (pt == eptVSite))
213 ri = atoms->atom[i].resind;
214 sprintf(buf, "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state B)\n"
215 " Check your topology.\n",
216 *(atoms->atomname[i]),
217 *(atoms->resinfo[ri].name),
218 atoms->resinfo[ri].nr,
220 warning_error(wi, buf);
221 /* The following statements make LINCS break! */
222 /* atoms->atom[i].m=0; */
229 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
231 * The results of this routine are only used for checking and for
232 * printing warning messages. Thus we can assume that charges of molecules
233 * should be integer. If the user wanted non-integer molecular charge,
234 * an undesired warning is printed and the user should use grompp -maxwarn 1.
236 * \param qMol The total, unrounded, charge of the molecule
237 * \param sumAbsQ The sum of absolute values of the charges, used for determining the tolerance for the rounding.
239 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
241 /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
242 * of the charges for ascii float truncation in the topology files.
243 * Although the summation here uses double precision, the charges
244 * are read and stored in single precision when real=float. This can
245 * lead to rounding errors of half the least significant bit.
246 * Note that, unfortunately, we can not assume addition of random
247 * rounding errors. It is not entirely unlikely that many charges
248 * have a near half-bit rounding error with the same sign.
250 double tolAbs = 1e-6;
251 double tol = std::max(tolAbs, 0.5*GMX_REAL_EPS*sumAbsQ);
252 double qRound = std::round(qMol);
253 if (std::abs(qMol - qRound) <= tol)
263 static void sum_q(const t_atoms *atoms, int numMols,
264 double *qTotA, double *qTotB)
271 for (int i = 0; i < atoms->nr; i++)
273 qmolA += atoms->atom[i].q;
274 qmolB += atoms->atom[i].qB;
275 sumAbsQA += std::abs(atoms->atom[i].q);
276 sumAbsQB += std::abs(atoms->atom[i].qB);
279 *qTotA += numMols*roundedMoleculeCharge(qmolA, sumAbsQA);
280 *qTotB += numMols*roundedMoleculeCharge(qmolB, sumAbsQB);
283 static void get_nbparm(char *nb_str, char *comb_str, int *nb, int *comb,
287 char warn_buf[STRLEN];
290 for (i = 1; (i < eNBF_NR); i++)
292 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
299 *nb = strtol(nb_str, nullptr, 10);
301 if ((*nb < 1) || (*nb >= eNBF_NR))
303 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s",
304 nb_str, enbf_names[1]);
305 warning_error(wi, warn_buf);
309 for (i = 1; (i < eCOMB_NR); i++)
311 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
318 *comb = strtol(comb_str, nullptr, 10);
320 if ((*comb < 1) || (*comb >= eCOMB_NR))
322 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s",
323 comb_str, ecomb_names[1]);
324 warning_error(wi, warn_buf);
329 static char ** cpp_opts(const char *define, const char *include,
334 const char *cppadds[2];
335 char **cppopts = nullptr;
336 const char *option[2] = { "-D", "-I" };
337 const char *nopt[2] = { "define", "include" };
341 char warn_buf[STRLEN];
344 cppadds[1] = include;
345 for (n = 0; (n < 2); n++)
352 while ((*ptr != '\0') && isspace(*ptr))
357 while ((*rptr != '\0') && !isspace(*rptr))
365 strncpy(buf, ptr, len);
366 if (strstr(ptr, option[n]) != ptr)
368 set_warning_line(wi, "mdp file", -1);
369 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
370 warning(wi, warn_buf);
374 srenew(cppopts, ++ncppopts);
375 cppopts[ncppopts-1] = gmx_strdup(buf);
383 srenew(cppopts, ++ncppopts);
384 cppopts[ncppopts-1] = nullptr;
390 static void make_atoms_sys(const std::vector<gmx_molblock_t> &molblock,
391 const t_molinfo *molinfo,
395 atoms->atom = nullptr;
397 for (const gmx_molblock_t &molb : molblock)
399 const t_atoms &mol_atoms = molinfo[molb.type].atoms;
401 srenew(atoms->atom, atoms->nr + molb.nmol*mol_atoms.nr);
403 for (int m = 0; m < molb.nmol; m++)
405 for (int a = 0; a < mol_atoms.nr; a++)
407 atoms->atom[atoms->nr++] = mol_atoms.atom[a];
414 static char **read_topol(const char *infile, const char *outfile,
415 const char *define, const char *include,
417 gpp_atomtype_t atype,
420 t_molinfo **intermolecular_interactions,
422 int *combination_rule,
426 std::vector<gmx_molblock_t> *molblock,
429 bool usingFullRangeElectrostatics,
434 char *pline = nullptr, **title = nullptr;
435 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
437 char *dirstr, *dummy2;
438 int nrcopies, nmol, nscan, ncombs, ncopy;
439 double fLJ, fQQ, fPOW;
440 t_molinfo *mi0 = nullptr;
443 t_nbparam **nbparam, **pair;
445 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
446 bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
447 double qt = 0, qBt = 0; /* total charge */
448 t_bond_atomtype batype;
450 int dcatt = -1, nmol_couple;
451 /* File handling variables */
455 char *tmp_line = nullptr;
456 char warn_buf[STRLEN];
457 const char *floating_point_arithmetic_tip =
458 "Total charge should normally be an integer. See\n"
459 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
460 "for discussion on how close it should be to an integer.\n";
461 /* We need to open the output file before opening the input file,
462 * because cpp_open_file can change the current working directory.
466 out = gmx_fio_fopen(outfile, "w");
473 /* open input file */
474 auto cpp_opts_return = cpp_opts(define, include, wi);
475 status = cpp_open_file(infile, &handle, cpp_opts_return);
478 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
481 /* some local variables */
482 DS_Init(&DS); /* directive stack */
483 nmol = 0; /* no molecules yet... */
484 d = d_invalid; /* first thing should be a directive */
485 nbparam = nullptr; /* The temporary non-bonded matrix */
486 pair = nullptr; /* The temporary pair interaction matrix */
487 block2 = nullptr; /* the extra exclusions */
490 *reppow = 12.0; /* Default value for repulsion power */
492 *intermolecular_interactions = nullptr;
494 /* Init the number of CMAP torsion angles and grid spacing */
495 plist[F_CMAP].grid_spacing = 0;
496 plist[F_CMAP].nc = 0;
498 bWarn_copy_A_B = bFEP;
500 batype = init_bond_atomtype();
501 /* parse the actual file */
502 bReadDefaults = FALSE;
504 bReadMolType = FALSE;
509 status = cpp_read_line(&handle, STRLEN, line);
510 done = (status == eCPP_EOF);
513 if (status != eCPP_OK)
515 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
519 fprintf(out, "%s\n", line);
522 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
524 pline = gmx_strdup(line);
526 /* Strip trailing '\' from pline, if it exists */
528 if ((sl > 0) && (pline[sl-1] == CONTINUE))
533 /* build one long line from several fragments - necessary for CMAP */
534 while (continuing(line))
536 status = cpp_read_line(&handle, STRLEN, line);
537 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
539 /* Since we depend on the '\' being present to continue to read, we copy line
540 * to a tmp string, strip the '\' from that string, and cat it to pline
542 tmp_line = gmx_strdup(line);
544 sl = strlen(tmp_line);
545 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE))
547 tmp_line[sl-1] = ' ';
550 done = (status == eCPP_EOF);
553 if (status != eCPP_OK)
555 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
559 fprintf(out, "%s\n", line);
563 srenew(pline, strlen(pline)+strlen(tmp_line)+1);
564 strcat(pline, tmp_line);
568 /* skip trailing and leading spaces and comment text */
569 strip_comment (pline);
572 /* if there is something left... */
573 if (static_cast<int>(strlen(pline)) > 0)
575 if (pline[0] == OPENDIR)
577 /* A directive on this line: copy the directive
578 * without the brackets into dirstr, then
579 * skip spaces and tabs on either side of directive
581 dirstr = gmx_strdup((pline+1));
582 if ((dummy2 = strchr (dirstr, CLOSEDIR)) != nullptr)
588 if ((newd = str2dir(dirstr)) == d_invalid)
590 sprintf(errbuf, "Invalid directive %s", dirstr);
591 warning_error(wi, errbuf);
595 /* Directive found */
596 if (DS_Check_Order (DS, newd))
603 /* we should print here which directives should have
604 been present, and which actually are */
605 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
606 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
610 if (d == d_intermolecular_interactions)
612 if (*intermolecular_interactions == nullptr)
614 /* We (mis)use the moleculetype processing
615 * to process the intermolecular interactions
616 * by making a "molecule" of the size of the system.
618 snew(*intermolecular_interactions, 1);
619 init_molinfo(*intermolecular_interactions);
620 mi0 = *intermolecular_interactions;
621 make_atoms_sys(*molblock, *molinfo,
628 else if (d != d_invalid)
630 /* Not a directive, just a plain string
631 * use a gigantic switch to decode,
632 * if there is a valid directive!
639 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
640 cpp_error(&handle, eCPP_SYNTAX));
642 bReadDefaults = TRUE;
643 nscan = sscanf(pline, "%s%s%s%lf%lf%lf",
644 nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
655 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
658 bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
659 if (nb_funct != eNBF_LJ && bGenPairs)
661 gmx_fatal(FARGS, "Generating pair parameters is only supported with LJ non-bonded interactions");
677 nb_funct = ifunc_index(d_nonbond_params, nb_funct);
681 push_at(symtab, atype, batype, pline, nb_funct,
682 &nbparam, bGenPairs ? &pair : nullptr, wi);
686 push_bt(d, plist, 2, nullptr, batype, pline, wi);
688 case d_constrainttypes:
689 push_bt(d, plist, 2, nullptr, batype, pline, wi);
694 push_nbt(d, pair, atype, pline, F_LJ14, wi);
698 push_bt(d, plist, 2, atype, nullptr, pline, wi);
702 push_bt(d, plist, 3, nullptr, batype, pline, wi);
704 case d_dihedraltypes:
705 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
706 push_dihedraltype(d, plist, batype, pline, wi);
709 case d_nonbond_params:
710 push_nbt(d, nbparam, atype, pline, nb_funct, wi);
715 srenew(block,nblock);
716 srenew(blockinfo,nblock);
717 blk0=&(block[nblock-1]);
718 bi0=&(blockinfo[nblock-1]);
721 push_molt(symtab,bi0,pline);
725 case d_implicit_genborn_params:
726 // Skip this line, so old topologies with
727 // GB parameters can be read.
730 case d_implicit_surface_params:
731 // Skip this line, so that any topologies
732 // with surface parameters can be read
733 // (even though these were never formally
738 push_cmaptype(d, plist, 5, atype, batype, pline, wi);
746 if (opts->couple_moltype != nullptr &&
747 (opts->couple_lam0 == ecouplamNONE ||
748 opts->couple_lam0 == ecouplamQ ||
749 opts->couple_lam1 == ecouplamNONE ||
750 opts->couple_lam1 == ecouplamQ))
752 dcatt = add_atomtype_decoupled(symtab, atype,
753 &nbparam, bGenPairs ? &pair : nullptr);
755 ntype = get_atomtype_ntypes(atype);
756 ncombs = (ntype*(ntype+1))/2;
757 generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi);
758 ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
760 fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
761 free_nbparam(nbparam, ntype);
764 gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule);
765 ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
767 fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
768 free_nbparam(pair, ntype);
770 /* Copy GBSA parameters to atomtype array? */
775 push_molt(symtab, &nmol, molinfo, pline, wi);
776 srenew(block2, nmol);
777 block2[nmol-1].nr = 0;
778 mi0 = &((*molinfo)[nmol-1]);
779 mi0->atoms.haveMass = TRUE;
780 mi0->atoms.haveCharge = TRUE;
781 mi0->atoms.haveType = TRUE;
782 mi0->atoms.haveBState = TRUE;
783 mi0->atoms.havePdbInfo = FALSE;
787 push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi);
791 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
792 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
795 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
796 FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
806 case d_position_restraints:
807 case d_angle_restraints:
808 case d_angle_restraints_z:
809 case d_distance_restraints:
810 case d_orientation_restraints:
811 case d_dihedral_restraints:
814 case d_water_polarization:
815 case d_thole_polarization:
816 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE,
817 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
820 push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi);
824 push_vsitesn(d, mi0->plist, &(mi0->atoms), pline, wi);
827 GMX_ASSERT(block2, "block2 must always be allocated so exclusions can be processed");
828 if (!block2[nmol-1].nr)
830 init_block2(&(block2[nmol-1]), mi0->atoms.nr);
832 push_excl(pline, &(block2[nmol-1]), wi);
836 title = put_symtab(symtab, pline);
843 push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
844 mi0 = &((*molinfo)[whichmol]);
845 molblock->resize(molblock->size() + 1);
846 molblock->back().type = whichmol;
847 molblock->back().nmol = nrcopies;
849 bCouple = (opts->couple_moltype != nullptr &&
850 (gmx_strcasecmp("system", opts->couple_moltype) == 0 ||
851 strcmp(*(mi0->name), opts->couple_moltype) == 0));
854 nmol_couple += nrcopies;
857 if (mi0->atoms.nr == 0)
859 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms",
863 "Excluding %d bonded neighbours molecule type '%s'\n",
864 mi0->nrexcl, *mi0->name);
865 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
866 if (!mi0->bProcessed)
869 generate_excl(mi0->nrexcl,
874 merge_excl(&(mi0->excls), &(block2[whichmol]), wi);
875 done_block2(&(block2[whichmol]));
876 make_shake(mi0->plist, &mi0->atoms, opts->nshake);
884 convert_moltype_couple(mi0, dcatt, *fudgeQQ,
885 opts->couple_lam0, opts->couple_lam1,
887 nb_funct, &(plist[nb_funct]), wi);
889 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
890 mi0->bProcessed = TRUE;
895 fprintf (stderr, "case: %d\n", static_cast<int>(d));
896 gmx_incons("unknown directive");
905 sfree(cpp_opts_return);
912 if (opts->couple_moltype)
914 if (nmol_couple == 0)
916 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling",
917 opts->couple_moltype);
919 fprintf(stderr, "Coupling %d copies of molecule type '%s'\n",
920 nmol_couple, opts->couple_moltype);
923 /* this is not very clean, but fixes core dump on empty system name */
926 title = put_symtab(symtab, "");
931 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
932 warning_note(wi, warn_buf);
934 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
936 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
937 warning_note(wi, warn_buf);
939 if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
941 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.");
942 please_cite(stdout, "Hub2014a");
946 for (i = 0; i < nmol; i++)
948 done_block2(&(block2[i]));
952 done_bond_atomtype(&batype);
954 if (*intermolecular_interactions != nullptr)
956 sfree(mi0->atoms.atom);
964 char **do_top(bool bVerbose,
966 const char *topppfile,
971 int *combination_rule,
972 double *repulsion_power,
974 gpp_atomtype_t atype,
977 t_molinfo **intermolecular_interactions,
978 const t_inputrec *ir,
979 std::vector<gmx_molblock_t> *molblock,
982 /* Tmpfile might contain a long path */
997 printf("processing topology...\n");
999 title = read_topol(topfile, tmpfile, opts->define, opts->include,
1001 nrmols, molinfo, intermolecular_interactions,
1002 plist, combination_rule, repulsion_power,
1003 opts, fudgeQQ, molblock,
1004 ir->efep != efepNO, bZero,
1005 EEL_FULL(ir->coulombtype), wi);
1007 if ((*combination_rule != eCOMB_GEOMETRIC) &&
1008 (ir->vdwtype == evdwUSER))
1010 warning(wi, "Using sigma/epsilon based combination rules with"
1011 " user supplied potential function may produce unwanted"
1019 static void generate_qmexcl_moltype(gmx_moltype_t *molt, const unsigned char *grpnr,
1020 t_inputrec *ir, warninp_t wi)
1022 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1024 /* generates the exclusions between the individual QM atoms, as
1025 * these interactions should be handled by the QM subroutines and
1026 * not by the gromacs routines
1028 int qm_max = 0, qm_nr = 0, link_nr = 0, link_max = 0;
1029 int *qm_arr = nullptr, *link_arr = nullptr;
1030 bool *bQMMM, *blink;
1032 /* First we search and select the QM atoms in an qm_arr array that
1033 * we use to create the exclusions.
1035 * we take the possibility into account that a user has defined more
1036 * than one QM group:
1038 * for that we also need to do this an ugly work-about just in case
1039 * the QM group contains the entire system...
1042 /* we first search for all the QM atoms and put them in an array
1044 for (int j = 0; j < ir->opts.ngQM; j++)
1046 for (int i = 0; i < molt->atoms.nr; i++)
1048 if (qm_nr >= qm_max)
1051 srenew(qm_arr, qm_max);
1053 if ((grpnr ? grpnr[i] : 0) == j)
1055 qm_arr[qm_nr++] = i;
1059 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1060 * QM/not QM. We first set all elements to false. Afterwards we use
1061 * the qm_arr to change the elements corresponding to the QM atoms
1064 snew(bQMMM, molt->atoms.nr);
1065 for (int i = 0; i < molt->atoms.nr; i++)
1069 for (int i = 0; i < qm_nr; i++)
1071 bQMMM[qm_arr[i]] = TRUE;
1074 /* We remove all bonded interactions (i.e. bonds,
1075 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1076 * are removed is as follows: if the interaction invloves 2 atoms,
1077 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1078 * it is removed if at least two of the atoms are QM atoms, if the
1079 * interaction involves 4 atoms, it is removed if there are at least
1080 * 2 QM atoms. Since this routine is called once before any forces
1081 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1082 * be rewritten at this poitn without any problem. 25-9-2002 */
1084 /* first check whether we already have CONNBONDS.
1085 * Note that if we don't, we don't add a param entry and set ftype=0,
1086 * which is ok, since CONNBONDS does not use parameters.
1088 int ftype_connbond = 0;
1089 int ind_connbond = 0;
1090 if (molt->ilist[F_CONNBONDS].nr != 0)
1092 fprintf(stderr, "nr. of CONNBONDS present already: %d\n",
1093 molt->ilist[F_CONNBONDS].nr/3);
1094 ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1095 ind_connbond = molt->ilist[F_CONNBONDS].nr;
1097 /* now we delete all bonded interactions, except the ones describing
1098 * a chemical bond. These are converted to CONNBONDS
1100 for (int ftype = 0; ftype < F_NRE; ftype++)
1102 if (!(interaction_function[ftype].flags & IF_BOND) ||
1103 ftype == F_CONNBONDS)
1107 int nratoms = interaction_function[ftype].nratoms;
1109 while (j < molt->ilist[ftype].nr)
1115 /* Remove an interaction between two atoms when both are
1116 * in the QM region. Note that we don't have to worry about
1117 * link atoms here, as they won't have 2-atom interactions.
1119 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1120 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1121 bexcl = (bQMMM[a1] && bQMMM[a2]);
1122 /* A chemical bond between two QM atoms will be copied to
1123 * the F_CONNBONDS list, for reasons mentioned above.
1125 if (bexcl && IS_CHEMBOND(ftype))
1127 srenew(molt->ilist[F_CONNBONDS].iatoms, ind_connbond + 3);
1128 molt->ilist[F_CONNBONDS].nr += 3;
1129 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = ftype_connbond;
1130 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a1;
1131 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a2;
1136 /* MM interactions have to be excluded if they are included
1137 * in the QM already. Because we use a link atom (H atom)
1138 * when the QM/MM boundary runs through a chemical bond, this
1139 * means that as long as one atom is MM, we still exclude,
1140 * as the interaction is included in the QM via:
1141 * QMatom1-QMatom2-QMatom-3-Linkatom.
1144 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1146 if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1151 bexcl = (numQmAtoms >= nratoms - 1);
1153 if (bexcl && ftype == F_SETTLE)
1155 gmx_fatal(FARGS, "Can not apply QM to molecules with SETTLE, replace the moleculetype using QM and SETTLE by one without SETTLE");
1160 /* since the interaction involves QM atoms, these should be
1161 * removed from the MM ilist
1163 molt->ilist[ftype].nr -= (nratoms+1);
1164 for (int l = j; l < molt->ilist[ftype].nr; l++)
1166 molt->ilist[ftype].iatoms[l] = molt->ilist[ftype].iatoms[l+(nratoms+1)];
1171 j += nratoms+1; /* the +1 is for the functype */
1175 /* Now, we search for atoms bonded to a QM atom because we also want
1176 * to exclude their nonbonded interactions with the QM atoms. The
1177 * reason for this is that this interaction is accounted for in the
1178 * linkatoms interaction with the QMatoms and would be counted
1181 for (int i = 0; i < F_NRE; i++)
1186 while (j < molt->ilist[i].nr)
1188 int a1 = molt->ilist[i].iatoms[j+1];
1189 int a2 = molt->ilist[i].iatoms[j+2];
1190 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1192 if (link_nr >= link_max)
1195 srenew(link_arr, link_max);
1199 link_arr[link_nr++] = a2;
1203 link_arr[link_nr++] = a1;
1210 snew(blink, molt->atoms.nr);
1211 for (int i = 0; i < molt->atoms.nr; i++)
1215 for (int i = 0; i < link_nr; i++)
1217 blink[link_arr[i]] = TRUE;
1219 /* creating the exclusion block for the QM atoms. Each QM atom has
1220 * as excluded elements all the other QMatoms (and itself).
1223 qmexcl.nr = molt->atoms.nr;
1224 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1225 snew(qmexcl.index, qmexcl.nr+1);
1226 snew(qmexcl.a, qmexcl.nra);
1228 for (int i = 0; i < qmexcl.nr; i++)
1230 qmexcl.index[i] = j;
1233 for (int k = 0; k < qm_nr; k++)
1235 qmexcl.a[k+j] = qm_arr[k];
1237 for (int k = 0; k < link_nr; k++)
1239 qmexcl.a[qm_nr+k+j] = link_arr[k];
1241 j += (qm_nr+link_nr);
1245 for (int k = 0; k < qm_nr; k++)
1247 qmexcl.a[k+j] = qm_arr[k];
1252 qmexcl.index[qmexcl.nr] = j;
1254 /* and merging with the exclusions already present in sys.
1258 init_block2(&qmexcl2, molt->atoms.nr);
1259 b_to_b2(&qmexcl, &qmexcl2);
1260 merge_excl(&(molt->excls), &qmexcl2, wi);
1261 done_block2(&qmexcl2);
1263 /* Finally, we also need to get rid of the pair interactions of the
1264 * classical atom bonded to the boundary QM atoms with the QMatoms,
1265 * as this interaction is already accounted for by the QM, so also
1266 * here we run the risk of double counting! We proceed in a similar
1267 * way as we did above for the other bonded interactions: */
1268 for (int i = F_LJ14; i < F_COUL14; i++)
1270 int nratoms = interaction_function[i].nratoms;
1272 while (j < molt->ilist[i].nr)
1274 int a1 = molt->ilist[i].iatoms[j+1];
1275 int a2 = molt->ilist[i].iatoms[j+2];
1276 bool bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1277 (blink[a1] && bQMMM[a2]) ||
1278 (bQMMM[a1] && blink[a2]));
1281 /* since the interaction involves QM atoms, these should be
1282 * removed from the MM ilist
1284 molt->ilist[i].nr -= (nratoms+1);
1285 for (int k = j; k < molt->ilist[i].nr; k++)
1287 molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1292 j += nratoms+1; /* the +1 is for the functype */
1301 } /* generate_qmexcl */
1303 void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi)
1305 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1308 unsigned char *grpnr;
1309 int mol, nat_mol, nr_mol_with_qm_atoms = 0;
1310 gmx_molblock_t *molb;
1313 grpnr = sys->groups.grpnr[egcQMMM];
1315 for (size_t mb = 0; mb < sys->molblock.size(); mb++)
1317 molb = &sys->molblock[mb];
1318 nat_mol = sys->moltype[molb->type].atoms.nr;
1319 for (mol = 0; mol < molb->nmol; mol++)
1322 for (int i = 0; i < nat_mol; i++)
1324 if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM)
1331 nr_mol_with_qm_atoms++;
1334 /* We need to split this molblock */
1337 /* Split the molblock at this molecule */
1338 auto pos = sys->molblock.begin() + mb + 1;
1339 sys->molblock.insert(pos, sys->molblock[mb]);
1340 sys->molblock[mb ].nmol = mol;
1341 sys->molblock[mb+1].nmol -= mol;
1343 molb = &sys->molblock[mb];
1347 /* Split the molblock after this molecule */
1348 auto pos = sys->molblock.begin() + mb + 1;
1349 sys->molblock.insert(pos, sys->molblock[mb]);
1350 molb = &sys->molblock[mb];
1351 sys->molblock[mb ].nmol = 1;
1352 sys->molblock[mb+1].nmol -= 1;
1355 /* Add a moltype for the QMMM molecule */
1356 sys->moltype.push_back(sys->moltype[molb->type]);
1357 /* Copy the exclusions to a new array, since this is the only
1358 * thing that needs to be modified for QMMM.
1360 copy_blocka(&sys->moltype[molb->type ].excls,
1361 &sys->moltype.back().excls);
1362 /* Set the molecule type for the QMMM molblock */
1363 molb->type = sys->moltype.size() - 1;
1365 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi);
1373 if (nr_mol_with_qm_atoms > 1)
1375 /* generate a warning is there are QM atoms in different
1376 * topologies. In this case it is not possible at this stage to
1377 * mutualy exclude the non-bonded interactions via the
1378 * exclusions (AFAIK). Instead, the user is advised to use the
1379 * energy group exclusions in the mdp file
1382 "\nThe QM subsystem is divided over multiple topologies. "
1383 "The mutual non-bonded interactions cannot be excluded. "
1384 "There are two ways to achieve this:\n\n"
1385 "1) merge the topologies, such that the atoms of the QM "
1386 "subsystem are all present in one single topology file. "
1387 "In this case this warning will dissappear\n\n"
1388 "2) exclude the non-bonded interactions explicitly via the "
1389 "energygrp-excl option in the mdp file. if this is the case "
1390 "this warning may be ignored"