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.
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_bond_atomtype.h"
58 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
59 #include "gromacs/gmxpreprocess/grompp-impl.h"
60 #include "gromacs/gmxpreprocess/readir.h"
61 #include "gromacs/gmxpreprocess/topdirs.h"
62 #include "gromacs/gmxpreprocess/toppush.h"
63 #include "gromacs/gmxpreprocess/topshake.h"
64 #include "gromacs/gmxpreprocess/toputil.h"
65 #include "gromacs/gmxpreprocess/vsite_parm.h"
66 #include "gromacs/math/units.h"
67 #include "gromacs/math/utilities.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/pbcutil/pbc.h"
71 #include "gromacs/topology/block.h"
72 #include "gromacs/topology/ifunc.h"
73 #include "gromacs/topology/symtab.h"
74 #include "gromacs/topology/topology.h"
75 #include "gromacs/utility/cstringutil.h"
76 #include "gromacs/utility/fatalerror.h"
77 #include "gromacs/utility/futil.h"
78 #include "gromacs/utility/gmxassert.h"
79 #include "gromacs/utility/pleasecite.h"
80 #include "gromacs/utility/smalloc.h"
82 #define OPENDIR '[' /* starting sign for directive */
83 #define CLOSEDIR ']' /* ending sign for directive */
85 static void free_nbparam(t_nbparam **param, int nr)
90 for (i = 0; i < nr; i++)
98 static int copy_nbparams(t_nbparam **param, int ftype, t_params *plist, int nr)
106 for (i = 0; i < nr; i++)
108 for (j = 0; j <= i; j++)
111 if (param[i][j].bSet)
113 for (f = 0; f < nrfp; f++)
115 plist->param[nr*i+j].c[f] = param[i][j].c[f];
116 plist->param[nr*j+i].c[f] = param[i][j].c[f];
126 static void gen_pairs(t_params *nbs, t_params *pairs, real fudge, int comb)
128 int i, j, ntp, nrfp, nrfpA, nrfpB, nnn;
131 nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
132 GMX_ASSERT(nnn * nnn == ntp, "Number of pairs of generated non-bonded parameters should be a perfect square");
134 nrfpA = interaction_function[F_LJ14].nrfpA;
135 nrfpB = interaction_function[F_LJ14].nrfpB;
138 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
140 gmx_incons("Number of force parameters in gen_pairs wrong");
143 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
146 fprintf(debug, "Fudge factor for 1-4 interactions: %g\n", fudge);
147 fprintf(debug, "Holy Cow! there are %d types\n", ntp);
149 snew(pairs->param, pairs->nr);
150 for (i = 0; (i < ntp); i++)
153 pairs->param[i].a[0] = i / nnn;
154 pairs->param[i].a[1] = i % nnn;
155 /* Copy normal and FEP parameters and multiply by fudge factor */
159 for (j = 0; (j < nrfp); j++)
161 /* If we are using sigma/epsilon values, only the epsilon values
162 * should be scaled, but not sigma.
163 * The sigma values have even indices 0,2, etc.
165 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2 == 0))
174 pairs->param[i].c[j] = scaling*nbs->param[i].c[j];
175 /* NOTE: this should be cleat to the compiler, but some gcc 5.2 versions
176 * issue false positive warnings for the pairs->param.c[] indexing below.
178 assert(2*nrfp <= MAXFORCEPARAM);
179 pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j];
184 double check_mol(gmx_mtop_t *mtop, warninp_t wi)
187 int i, mb, nmol, ri, pt;
192 /* Check mass and charge */
195 for (mb = 0; mb < mtop->nmolblock; mb++)
197 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
198 nmol = mtop->molblock[mb].nmol;
199 for (i = 0; (i < atoms->nr); i++)
201 q += nmol*atoms->atom[i].q;
202 m = atoms->atom[i].m;
203 mB = atoms->atom[i].mB;
204 pt = atoms->atom[i].ptype;
205 /* If the particle is an atom or a nucleus it must have a mass,
206 * else, if it is a shell, a vsite or a bondshell it can have mass zero
208 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
210 ri = atoms->atom[i].resind;
211 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
212 *(atoms->atomname[i]),
213 *(atoms->resinfo[ri].name),
214 atoms->resinfo[ri].nr,
216 warning_error(wi, buf);
219 if (((m != 0) || (mB != 0)) && (pt == eptVSite))
221 ri = atoms->atom[i].resind;
222 sprintf(buf, "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state B)\n"
223 " Check your topology.\n",
224 *(atoms->atomname[i]),
225 *(atoms->resinfo[ri].name),
226 atoms->resinfo[ri].nr,
228 warning_error(wi, buf);
229 /* The following statements make LINCS break! */
230 /* atoms->atom[i].m=0; */
237 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
239 * The results of this routine are only used for checking and for
240 * printing warning messages. Thus we can assume that charges of molecules
241 * should be integer. If the user wanted non-integer molecular charge,
242 * an undesired warning is printed and the user should use grompp -maxwarn 1.
244 * \param qMol The total, unrounded, charge of the molecule
245 * \param sumAbsQ The sum of absolute values of the charges, used for determining the tolerance for the rounding.
247 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
249 /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
250 * of the charges for ascii float truncation in the topology files.
251 * Although the summation here uses double precision, the charges
252 * are read and stored in single precision when real=float. This can
253 * lead to rounding errors of half the least significant bit.
254 * Note that, unfortunately, we can not assume addition of random
255 * rounding errors. It is not entirely unlikely that many charges
256 * have a near half-bit rounding error with the same sign.
258 double tolAbs = 1e-6;
259 double tol = std::max(tolAbs, 0.5*GMX_REAL_EPS*sumAbsQ);
260 double qRound = std::round(qMol);
261 if (std::abs(qMol - qRound) <= tol)
271 static void sum_q(const t_atoms *atoms, int numMols,
272 double *qTotA, double *qTotB)
279 for (int i = 0; i < atoms->nr; i++)
281 qmolA += atoms->atom[i].q;
282 qmolB += atoms->atom[i].qB;
283 sumAbsQA += std::abs(atoms->atom[i].q);
284 sumAbsQB += std::abs(atoms->atom[i].qB);
287 *qTotA += numMols*roundedMoleculeCharge(qmolA, sumAbsQA);
288 *qTotB += numMols*roundedMoleculeCharge(qmolB, sumAbsQB);
291 static void get_nbparm(char *nb_str, char *comb_str, int *nb, int *comb,
295 char warn_buf[STRLEN];
298 for (i = 1; (i < eNBF_NR); i++)
300 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
307 *nb = strtol(nb_str, nullptr, 10);
309 if ((*nb < 1) || (*nb >= eNBF_NR))
311 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s",
312 nb_str, enbf_names[1]);
313 warning_error(wi, warn_buf);
317 for (i = 1; (i < eCOMB_NR); i++)
319 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
326 *comb = strtol(comb_str, nullptr, 10);
328 if ((*comb < 1) || (*comb >= eCOMB_NR))
330 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s",
331 comb_str, ecomb_names[1]);
332 warning_error(wi, warn_buf);
337 static char ** cpp_opts(const char *define, const char *include,
342 const char *cppadds[2];
343 char **cppopts = nullptr;
344 const char *option[2] = { "-D", "-I" };
345 const char *nopt[2] = { "define", "include" };
349 char warn_buf[STRLEN];
352 cppadds[1] = include;
353 for (n = 0; (n < 2); n++)
360 while ((*ptr != '\0') && isspace(*ptr))
365 while ((*rptr != '\0') && !isspace(*rptr))
373 strncpy(buf, ptr, len);
374 if (strstr(ptr, option[n]) != ptr)
376 set_warning_line(wi, "mdp file", -1);
377 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
378 warning(wi, warn_buf);
382 srenew(cppopts, ++ncppopts);
383 cppopts[ncppopts-1] = gmx_strdup(buf);
391 srenew(cppopts, ++ncppopts);
392 cppopts[ncppopts-1] = nullptr;
398 static void make_atoms_sys(int nmolb, const gmx_molblock_t *molb,
399 const t_molinfo *molinfo,
403 const t_atoms *mol_atoms;
406 atoms->atom = nullptr;
408 for (mb = 0; mb < nmolb; mb++)
411 mol_atoms = &molinfo[molb[mb].type].atoms;
413 srenew(atoms->atom, atoms->nr + molb[mb].nmol*mol_atoms->nr);
415 for (m = 0; m < molb[mb].nmol; m++)
417 for (a = 0; a < mol_atoms->nr; a++)
419 atoms->atom[atoms->nr++] = mol_atoms->atom[a];
426 static char **read_topol(const char *infile, const char *outfile,
427 const char *define, const char *include,
429 gpp_atomtype_t atype,
432 t_molinfo **intermolecular_interactions,
434 int *combination_rule,
439 gmx_molblock_t **molblock,
442 gmx_bool usingFullRangeElectrostatics,
447 char *pline = nullptr, **title = nullptr;
448 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
450 char *dirstr, *dummy2;
451 int nrcopies, nmol, nmolb = 0, nscan, ncombs, ncopy;
452 double fLJ, fQQ, fPOW;
453 gmx_molblock_t *molb = nullptr;
454 t_molinfo *mi0 = nullptr;
457 t_nbparam **nbparam, **pair;
459 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
460 gmx_bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
461 double qt = 0, qBt = 0; /* total charge */
462 t_bond_atomtype batype;
464 int dcatt = -1, nmol_couple;
465 /* File handling variables */
468 char *tmp_line = nullptr;
469 char warn_buf[STRLEN];
470 const char *floating_point_arithmetic_tip =
471 "Total charge should normally be an integer. See\n"
472 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
473 "for discussion on how close it should be to an integer.\n";
474 /* We need to open the output file before opening the input file,
475 * because cpp_open_file can change the current working directory.
479 out = gmx_fio_fopen(outfile, "w");
486 /* open input file */
487 status = cpp_open_file(infile, &handle, cpp_opts(define, include, wi));
490 gmx_fatal(FARGS, cpp_error(&handle, status));
493 /* some local variables */
494 DS_Init(&DS); /* directive stack */
495 nmol = 0; /* no molecules yet... */
496 d = d_invalid; /* first thing should be a directive */
497 nbparam = nullptr; /* The temporary non-bonded matrix */
498 pair = nullptr; /* The temporary pair interaction matrix */
499 block2 = nullptr; /* the extra exclusions */
502 *reppow = 12.0; /* Default value for repulsion power */
504 *intermolecular_interactions = nullptr;
506 /* Init the number of CMAP torsion angles and grid spacing */
507 plist[F_CMAP].grid_spacing = 0;
508 plist[F_CMAP].nc = 0;
510 bWarn_copy_A_B = bFEP;
512 batype = init_bond_atomtype();
513 /* parse the actual file */
514 bReadDefaults = FALSE;
516 bReadMolType = FALSE;
521 status = cpp_read_line(&handle, STRLEN, line);
522 done = (status == eCPP_EOF);
525 if (status != eCPP_OK)
527 gmx_fatal(FARGS, cpp_error(&handle, status));
531 fprintf(out, "%s\n", line);
534 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
536 pline = gmx_strdup(line);
538 /* Strip trailing '\' from pline, if it exists */
540 if ((sl > 0) && (pline[sl-1] == CONTINUE))
545 /* build one long line from several fragments - necessary for CMAP */
546 while (continuing(line))
548 status = cpp_read_line(&handle, STRLEN, line);
549 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
551 /* Since we depend on the '\' being present to continue to read, we copy line
552 * to a tmp string, strip the '\' from that string, and cat it to pline
554 tmp_line = gmx_strdup(line);
556 sl = strlen(tmp_line);
557 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE))
559 tmp_line[sl-1] = ' ';
562 done = (status == eCPP_EOF);
565 if (status != eCPP_OK)
567 gmx_fatal(FARGS, cpp_error(&handle, status));
571 fprintf(out, "%s\n", line);
575 srenew(pline, strlen(pline)+strlen(tmp_line)+1);
576 strcat(pline, tmp_line);
580 /* skip trailing and leading spaces and comment text */
581 strip_comment (pline);
584 /* if there is something left... */
585 if ((int)strlen(pline) > 0)
587 if (pline[0] == OPENDIR)
589 /* A directive on this line: copy the directive
590 * without the brackets into dirstr, then
591 * skip spaces and tabs on either side of directive
593 dirstr = gmx_strdup((pline+1));
594 if ((dummy2 = strchr (dirstr, CLOSEDIR)) != nullptr)
600 if ((newd = str2dir(dirstr)) == d_invalid)
602 sprintf(errbuf, "Invalid directive %s", dirstr);
603 warning_error(wi, errbuf);
607 /* Directive found */
610 fprintf(debug, "found directive '%s'\n", dir2str(newd));
612 if (DS_Check_Order (DS, newd))
619 /* we should print here which directives should have
620 been present, and which actually are */
621 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
622 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
626 if (d == d_intermolecular_interactions)
628 if (*intermolecular_interactions == nullptr)
630 /* We (mis)use the moleculetype processing
631 * to process the intermolecular interactions
632 * by making a "molecule" of the size of the system.
634 snew(*intermolecular_interactions, 1);
635 init_molinfo(*intermolecular_interactions);
636 mi0 = *intermolecular_interactions;
637 make_atoms_sys(nmolb, molb, *molinfo,
644 else if (d != d_invalid)
646 /* Not a directive, just a plain string
647 * use a gigantic switch to decode,
648 * if there is a valid directive!
655 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
656 cpp_error(&handle, eCPP_SYNTAX));
658 bReadDefaults = TRUE;
659 nscan = sscanf(pline, "%s%s%s%lf%lf%lf",
660 nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
671 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
674 bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
675 if (nb_funct != eNBF_LJ && bGenPairs)
677 gmx_fatal(FARGS, "Generating pair parameters is only supported with LJ non-bonded interactions");
693 nb_funct = ifunc_index(d_nonbond_params, nb_funct);
697 push_at(symtab, atype, batype, pline, nb_funct,
698 &nbparam, bGenPairs ? &pair : nullptr, wi);
702 push_bt(d, plist, 2, nullptr, batype, pline, wi);
704 case d_constrainttypes:
705 push_bt(d, plist, 2, nullptr, batype, pline, wi);
710 push_nbt(d, pair, atype, pline, F_LJ14, wi);
714 push_bt(d, plist, 2, atype, nullptr, pline, wi);
718 push_bt(d, plist, 3, nullptr, batype, pline, wi);
720 case d_dihedraltypes:
721 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
722 push_dihedraltype(d, plist, batype, pline, wi);
725 case d_nonbond_params:
726 push_nbt(d, nbparam, atype, pline, nb_funct, wi);
731 srenew(block,nblock);
732 srenew(blockinfo,nblock);
733 blk0=&(block[nblock-1]);
734 bi0=&(blockinfo[nblock-1]);
737 push_molt(symtab,bi0,pline);
741 case d_implicit_genborn_params:
742 // Skip this line, so old topologies with
743 // GB parameters can be read.
746 case d_implicit_surface_params:
747 // Skip this line, so that any topologies
748 // with surface parameters can be read
749 // (even though these were never formally
754 push_cmaptype(d, plist, 5, atype, batype, pline, wi);
762 if (opts->couple_moltype != nullptr &&
763 (opts->couple_lam0 == ecouplamNONE ||
764 opts->couple_lam0 == ecouplamQ ||
765 opts->couple_lam1 == ecouplamNONE ||
766 opts->couple_lam1 == ecouplamQ))
768 dcatt = add_atomtype_decoupled(symtab, atype,
769 &nbparam, bGenPairs ? &pair : nullptr);
771 ntype = get_atomtype_ntypes(atype);
772 ncombs = (ntype*(ntype+1))/2;
773 generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi);
774 ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
776 fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
777 free_nbparam(nbparam, ntype);
780 gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule);
781 ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
783 fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
784 free_nbparam(pair, ntype);
786 /* Copy GBSA parameters to atomtype array? */
791 push_molt(symtab, &nmol, molinfo, pline, wi);
792 srenew(block2, nmol);
793 block2[nmol-1].nr = 0;
794 mi0 = &((*molinfo)[nmol-1]);
795 mi0->atoms.haveMass = TRUE;
796 mi0->atoms.haveCharge = TRUE;
797 mi0->atoms.haveType = TRUE;
798 mi0->atoms.haveBState = TRUE;
799 mi0->atoms.havePdbInfo = FALSE;
803 push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi);
807 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
808 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
811 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
812 FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
822 case d_position_restraints:
823 case d_angle_restraints:
824 case d_angle_restraints_z:
825 case d_distance_restraints:
826 case d_orientation_restraints:
827 case d_dihedral_restraints:
830 case d_water_polarization:
831 case d_thole_polarization:
832 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE,
833 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
836 push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi);
840 push_vsitesn(d, mi0->plist, &(mi0->atoms), pline, wi);
843 GMX_ASSERT(block2, "block2 must always be allocated so exclusions can be processed");
844 if (!block2[nmol-1].nr)
846 init_block2(&(block2[nmol-1]), mi0->atoms.nr);
848 push_excl(pline, &(block2[nmol-1]), wi);
852 title = put_symtab(symtab, pline);
859 push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
860 mi0 = &((*molinfo)[whichmol]);
861 srenew(molb, nmolb+1);
862 molb[nmolb].type = whichmol;
863 molb[nmolb].nmol = nrcopies;
866 bCouple = (opts->couple_moltype != nullptr &&
867 (gmx_strcasecmp("system", opts->couple_moltype) == 0 ||
868 strcmp(*(mi0->name), opts->couple_moltype) == 0));
871 nmol_couple += nrcopies;
874 if (mi0->atoms.nr == 0)
876 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms",
880 "Excluding %d bonded neighbours molecule type '%s'\n",
881 mi0->nrexcl, *mi0->name);
882 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
883 if (!mi0->bProcessed)
886 generate_excl(mi0->nrexcl,
891 merge_excl(&(mi0->excls), &(block2[whichmol]), wi);
892 done_block2(&(block2[whichmol]));
893 make_shake(mi0->plist, &mi0->atoms, opts->nshake);
901 convert_moltype_couple(mi0, dcatt, *fudgeQQ,
902 opts->couple_lam0, opts->couple_lam1,
904 nb_funct, &(plist[nb_funct]), wi);
906 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
907 mi0->bProcessed = TRUE;
912 fprintf (stderr, "case: %d\n", (int)d);
913 gmx_incons("unknown directive");
922 status = cpp_close_file(&handle);
923 if (status != eCPP_OK)
925 gmx_fatal(FARGS, cpp_error(&handle, status));
933 if (opts->couple_moltype)
935 if (nmol_couple == 0)
937 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling",
938 opts->couple_moltype);
940 fprintf(stderr, "Coupling %d copies of molecule type '%s'\n",
941 nmol_couple, opts->couple_moltype);
944 /* this is not very clean, but fixes core dump on empty system name */
947 title = put_symtab(symtab, "");
952 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
953 warning_note(wi, warn_buf);
955 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
957 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
958 warning_note(wi, warn_buf);
960 if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
962 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.");
963 please_cite(stdout, "Hub2014a");
967 for (i = 0; i < nmol; i++)
969 done_block2(&(block2[i]));
973 done_bond_atomtype(&batype);
975 if (*intermolecular_interactions != nullptr)
977 sfree(mi0->atoms.atom);
988 char **do_top(gmx_bool bVerbose,
990 const char *topppfile,
995 int *combination_rule,
996 double *repulsion_power,
998 gpp_atomtype_t atype,
1000 t_molinfo **molinfo,
1001 t_molinfo **intermolecular_interactions,
1002 const t_inputrec *ir,
1004 gmx_molblock_t **molblock,
1007 /* Tmpfile might contain a long path */
1008 const char *tmpfile;
1013 tmpfile = topppfile;
1022 printf("processing topology...\n");
1024 title = read_topol(topfile, tmpfile, opts->define, opts->include,
1026 nrmols, molinfo, intermolecular_interactions,
1027 plist, combination_rule, repulsion_power,
1028 opts, fudgeQQ, nmolblock, molblock,
1029 ir->efep != efepNO, bZero,
1030 EEL_FULL(ir->coulombtype), wi);
1032 if ((*combination_rule != eCOMB_GEOMETRIC) &&
1033 (ir->vdwtype == evdwUSER))
1035 warning(wi, "Using sigma/epsilon based combination rules with"
1036 " user supplied potential function may produce unwanted"
1044 static void generate_qmexcl_moltype(gmx_moltype_t *molt, unsigned char *grpnr,
1045 t_inputrec *ir, warninp_t wi)
1047 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1049 /* generates the exclusions between the individual QM atoms, as
1050 * these interactions should be handled by the QM subroutines and
1051 * not by the gromacs routines
1053 int qm_max = 0, qm_nr = 0, link_nr = 0, link_max = 0;
1054 int *qm_arr = nullptr, *link_arr = nullptr;
1055 gmx_bool *bQMMM, *blink;
1057 /* First we search and select the QM atoms in an qm_arr array that
1058 * we use to create the exclusions.
1060 * we take the possibility into account that a user has defined more
1061 * than one QM group:
1063 * for that we also need to do this an ugly work-about just in case
1064 * the QM group contains the entire system...
1067 /* we first search for all the QM atoms and put them in an array
1069 for (int j = 0; j < ir->opts.ngQM; j++)
1071 for (int i = 0; i < molt->atoms.nr; i++)
1073 if (qm_nr >= qm_max)
1076 srenew(qm_arr, qm_max);
1078 if ((grpnr ? grpnr[i] : 0) == j)
1080 qm_arr[qm_nr++] = i;
1084 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1085 * QM/not QM. We first set all elements to false. Afterwards we use
1086 * the qm_arr to change the elements corresponding to the QM atoms
1089 snew(bQMMM, molt->atoms.nr);
1090 for (int i = 0; i < molt->atoms.nr; i++)
1094 for (int i = 0; i < qm_nr; i++)
1096 bQMMM[qm_arr[i]] = TRUE;
1099 /* We remove all bonded interactions (i.e. bonds,
1100 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1101 * are removed is as follows: if the interaction invloves 2 atoms,
1102 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1103 * it is removed if at least two of the atoms are QM atoms, if the
1104 * interaction involves 4 atoms, it is removed if there are at least
1105 * 2 QM atoms. Since this routine is called once before any forces
1106 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1107 * be rewritten at this poitn without any problem. 25-9-2002 */
1109 /* first check whether we already have CONNBONDS.
1110 * Note that if we don't, we don't add a param entry and set ftype=0,
1111 * which is ok, since CONNBONDS does not use parameters.
1113 int ftype_connbond = 0;
1114 int ind_connbond = 0;
1115 if (molt->ilist[F_CONNBONDS].nr != 0)
1117 fprintf(stderr, "nr. of CONNBONDS present already: %d\n",
1118 molt->ilist[F_CONNBONDS].nr/3);
1119 ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1120 ind_connbond = molt->ilist[F_CONNBONDS].nr;
1122 /* now we delete all bonded interactions, except the ones describing
1123 * a chemical bond. These are converted to CONNBONDS
1125 for (int ftype = 0; ftype < F_NRE; ftype++)
1127 if (!(interaction_function[ftype].flags & IF_BOND) ||
1128 ftype == F_CONNBONDS)
1132 int nratoms = interaction_function[ftype].nratoms;
1134 while (j < molt->ilist[ftype].nr)
1140 /* Remove an interaction between two atoms when both are
1141 * in the QM region. Note that we don't have to worry about
1142 * link atoms here, as they won't have 2-atom interactions.
1144 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1145 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1146 bexcl = (bQMMM[a1] && bQMMM[a2]);
1147 /* A chemical bond between two QM atoms will be copied to
1148 * the F_CONNBONDS list, for reasons mentioned above.
1150 if (bexcl && IS_CHEMBOND(ftype))
1152 srenew(molt->ilist[F_CONNBONDS].iatoms, ind_connbond + 3);
1153 molt->ilist[F_CONNBONDS].nr += 3;
1154 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = ftype_connbond;
1155 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a1;
1156 molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a2;
1161 /* MM interactions have to be excluded if they are included
1162 * in the QM already. Because we use a link atom (H atom)
1163 * when the QM/MM boundary runs through a chemical bond, this
1164 * means that as long as one atom is MM, we still exclude,
1165 * as the interaction is included in the QM via:
1166 * QMatom1-QMatom2-QMatom-3-Linkatom.
1169 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1171 if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1176 bexcl = (numQmAtoms >= nratoms - 1);
1178 if (bexcl && ftype == F_SETTLE)
1180 gmx_fatal(FARGS, "Can not apply QM to molecules with SETTLE, replace the moleculetype using QM and SETTLE by one without SETTLE");
1185 /* since the interaction involves QM atoms, these should be
1186 * removed from the MM ilist
1188 molt->ilist[ftype].nr -= (nratoms+1);
1189 for (int l = j; l < molt->ilist[ftype].nr; l++)
1191 molt->ilist[ftype].iatoms[l] = molt->ilist[ftype].iatoms[l+(nratoms+1)];
1196 j += nratoms+1; /* the +1 is for the functype */
1200 /* Now, we search for atoms bonded to a QM atom because we also want
1201 * to exclude their nonbonded interactions with the QM atoms. The
1202 * reason for this is that this interaction is accounted for in the
1203 * linkatoms interaction with the QMatoms and would be counted
1206 for (int i = 0; i < F_NRE; i++)
1211 while (j < molt->ilist[i].nr)
1213 int a1 = molt->ilist[i].iatoms[j+1];
1214 int a2 = molt->ilist[i].iatoms[j+2];
1215 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1217 if (link_nr >= link_max)
1220 srenew(link_arr, link_max);
1224 link_arr[link_nr++] = a2;
1228 link_arr[link_nr++] = a1;
1235 snew(blink, molt->atoms.nr);
1236 for (int i = 0; i < molt->atoms.nr; i++)
1240 for (int i = 0; i < link_nr; i++)
1242 blink[link_arr[i]] = TRUE;
1244 /* creating the exclusion block for the QM atoms. Each QM atom has
1245 * as excluded elements all the other QMatoms (and itself).
1248 qmexcl.nr = molt->atoms.nr;
1249 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1250 snew(qmexcl.index, qmexcl.nr+1);
1251 snew(qmexcl.a, qmexcl.nra);
1253 for (int i = 0; i < qmexcl.nr; i++)
1255 qmexcl.index[i] = j;
1258 for (int k = 0; k < qm_nr; k++)
1260 qmexcl.a[k+j] = qm_arr[k];
1262 for (int k = 0; k < link_nr; k++)
1264 qmexcl.a[qm_nr+k+j] = link_arr[k];
1266 j += (qm_nr+link_nr);
1270 for (int k = 0; k < qm_nr; k++)
1272 qmexcl.a[k+j] = qm_arr[k];
1277 qmexcl.index[qmexcl.nr] = j;
1279 /* and merging with the exclusions already present in sys.
1283 init_block2(&qmexcl2, molt->atoms.nr);
1284 b_to_b2(&qmexcl, &qmexcl2);
1285 merge_excl(&(molt->excls), &qmexcl2, wi);
1286 done_block2(&qmexcl2);
1288 /* Finally, we also need to get rid of the pair interactions of the
1289 * classical atom bonded to the boundary QM atoms with the QMatoms,
1290 * as this interaction is already accounted for by the QM, so also
1291 * here we run the risk of double counting! We proceed in a similar
1292 * way as we did above for the other bonded interactions: */
1293 for (int i = F_LJ14; i < F_COUL14; i++)
1295 int nratoms = interaction_function[i].nratoms;
1297 while (j < molt->ilist[i].nr)
1299 int a1 = molt->ilist[i].iatoms[j+1];
1300 int a2 = molt->ilist[i].iatoms[j+2];
1301 bool bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1302 (blink[a1] && bQMMM[a2]) ||
1303 (bQMMM[a1] && blink[a2]));
1306 /* since the interaction involves QM atoms, these should be
1307 * removed from the MM ilist
1309 molt->ilist[i].nr -= (nratoms+1);
1310 for (int k = j; k < molt->ilist[i].nr; k++)
1312 molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1317 j += nratoms+1; /* the +1 is for the functype */
1326 } /* generate_qmexcl */
1328 void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi)
1330 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1333 unsigned char *grpnr;
1334 int mb, mol, nat_mol, i, nr_mol_with_qm_atoms = 0;
1335 gmx_molblock_t *molb;
1338 grpnr = sys->groups.grpnr[egcQMMM];
1340 for (mb = 0; mb < sys->nmolblock; mb++)
1342 molb = &sys->molblock[mb];
1343 nat_mol = sys->moltype[molb->type].atoms.nr;
1344 for (mol = 0; mol < molb->nmol; mol++)
1347 for (i = 0; i < nat_mol; i++)
1349 if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM)
1356 nr_mol_with_qm_atoms++;
1359 /* We need to split this molblock */
1362 /* Split the molblock at this molecule */
1364 srenew(sys->molblock, sys->nmolblock);
1365 for (i = sys->nmolblock-2; i >= mb; i--)
1367 sys->molblock[i+1] = sys->molblock[i];
1369 sys->molblock[mb ].nmol = mol;
1370 sys->molblock[mb+1].nmol -= mol;
1372 molb = &sys->molblock[mb];
1376 /* Split the molblock after this molecule */
1378 srenew(sys->molblock, sys->nmolblock);
1379 molb = &sys->molblock[mb];
1380 for (i = sys->nmolblock-2; i >= mb; i--)
1382 sys->molblock[i+1] = sys->molblock[i];
1384 sys->molblock[mb ].nmol = 1;
1385 sys->molblock[mb+1].nmol -= 1;
1388 /* Add a moltype for the QMMM molecule */
1390 srenew(sys->moltype, sys->nmoltype);
1391 /* Copy the moltype struct */
1392 sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
1393 /* Copy the exclusions to a new array, since this is the only
1394 * thing that needs to be modified for QMMM.
1396 copy_blocka(&sys->moltype[molb->type ].excls,
1397 &sys->moltype[sys->nmoltype-1].excls);
1398 /* Set the molecule type for the QMMM molblock */
1399 molb->type = sys->nmoltype - 1;
1401 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, wi);
1409 if (nr_mol_with_qm_atoms > 1)
1411 /* generate a warning is there are QM atoms in different
1412 * topologies. In this case it is not possible at this stage to
1413 * mutualy exclude the non-bonded interactions via the
1414 * exclusions (AFAIK). Instead, the user is advised to use the
1415 * energy group exclusions in the mdp file
1418 "\nThe QM subsystem is divided over multiple topologies. "
1419 "The mutual non-bonded interactions cannot be excluded. "
1420 "There are two ways to achieve this:\n\n"
1421 "1) merge the topologies, such that the atoms of the QM "
1422 "subsystem are all present in one single topology file. "
1423 "In this case this warning will dissappear\n\n"
1424 "2) exclude the non-bonded interactions explicitly via the "
1425 "energygrp-excl option in the mdp file. if this is the case "
1426 "this warning may be ignored"