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, 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.
49 #include <sys/types.h>
51 #include "gromacs/fileio/gmxfio.h"
52 #include "gromacs/gmxpreprocess/gmxcpp.h"
53 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
54 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
55 #include "gromacs/gmxpreprocess/grompp-impl.h"
56 #include "gromacs/gmxpreprocess/topdirs.h"
57 #include "gromacs/gmxpreprocess/toppush.h"
58 #include "gromacs/gmxpreprocess/topshake.h"
59 #include "gromacs/gmxpreprocess/toputil.h"
60 #include "gromacs/gmxpreprocess/vsite_parm.h"
61 #include "gromacs/legacyheaders/genborn.h"
62 #include "gromacs/legacyheaders/macros.h"
63 #include "gromacs/legacyheaders/names.h"
64 #include "gromacs/legacyheaders/txtdump.h"
65 #include "gromacs/legacyheaders/typedefs.h"
66 #include "gromacs/legacyheaders/warninp.h"
67 #include "gromacs/math/units.h"
68 #include "gromacs/math/utilities.h"
69 #include "gromacs/topology/block.h"
70 #include "gromacs/topology/symtab.h"
71 #include "gromacs/topology/topology.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/futil.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/smalloc.h"
78 #define OPENDIR '[' /* starting sign for directive */
79 #define CLOSEDIR ']' /* ending sign for directive */
81 static void free_nbparam(t_nbparam **param, int nr)
85 for (i = 0; i < nr; i++)
92 static int copy_nbparams(t_nbparam **param, int ftype, t_params *plist, int nr)
100 for (i = 0; i < nr; i++)
102 for (j = 0; j <= i; j++)
104 if (param[i][j].bSet)
106 for (f = 0; f < nrfp; f++)
108 plist->param[nr*i+j].c[f] = param[i][j].c[f];
109 plist->param[nr*j+i].c[f] = param[i][j].c[f];
119 static void gen_pairs(t_params *nbs, t_params *pairs, real fudge, int comb)
121 int i, j, ntp, nrfp, nrfpA, nrfpB, nnn;
124 nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
125 GMX_ASSERT(nnn * nnn == ntp, "Number of pairs of generated non-bonded parameters should be a perfect square");
127 nrfpA = interaction_function[F_LJ14].nrfpA;
128 nrfpB = interaction_function[F_LJ14].nrfpB;
131 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
133 gmx_incons("Number of force parameters in gen_pairs wrong");
136 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
139 fprintf(debug, "Fudge factor for 1-4 interactions: %g\n", fudge);
140 fprintf(debug, "Holy Cow! there are %d types\n", ntp);
142 snew(pairs->param, pairs->nr);
143 for (i = 0; (i < ntp); i++)
146 pairs->param[i].a[0] = i / nnn;
147 pairs->param[i].a[1] = i % nnn;
148 /* Copy normal and FEP parameters and multiply by fudge factor */
152 for (j = 0; (j < nrfp); j++)
154 /* If we are using sigma/epsilon values, only the epsilon values
155 * should be scaled, but not sigma.
156 * The sigma values have even indices 0,2, etc.
158 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2 == 0))
167 pairs->param[i].c[j] = scaling*nbs->param[i].c[j];
168 pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j];
173 double check_mol(gmx_mtop_t *mtop, warninp_t wi)
176 int i, mb, nmol, ri, pt;
181 /* Check mass and charge */
184 for (mb = 0; mb < mtop->nmoltype; mb++)
186 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
187 nmol = mtop->molblock[mb].nmol;
188 for (i = 0; (i < atoms->nr); i++)
190 q += nmol*atoms->atom[i].q;
191 m = atoms->atom[i].m;
192 mB = atoms->atom[i].mB;
193 pt = atoms->atom[i].ptype;
194 /* If the particle is an atom or a nucleus it must have a mass,
195 * else, if it is a shell, a vsite or a bondshell it can have mass zero
197 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
199 ri = atoms->atom[i].resind;
200 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
201 *(atoms->atomname[i]),
202 *(atoms->resinfo[ri].name),
203 atoms->resinfo[ri].nr,
205 warning_error(wi, buf);
208 if (((m != 0) || (mB != 0)) && (pt == eptVSite))
210 ri = atoms->atom[i].resind;
211 sprintf(buf, "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state B)\n"
212 " Check your topology.\n",
213 *(atoms->atomname[i]),
214 *(atoms->resinfo[ri].name),
215 atoms->resinfo[ri].nr,
217 warning_error(wi, buf);
218 /* The following statements make LINCS break! */
219 /* atoms->atom[i].m=0; */
226 static void sum_q(t_atoms *atoms, int n, double *qt, double *qBt)
234 for (i = 0; i < atoms->nr; i++)
236 qmolA += atoms->atom[i].q;
237 qmolB += atoms->atom[i].qB;
239 /* Unfortunately an absolute comparison,
240 * but this avoids unnecessary warnings and gmx-users mails.
242 if (fabs(qmolA) >= 1e-6 || fabs(qmolB) >= 1e-6)
249 static void get_nbparm(char *nb_str, char *comb_str, int *nb, int *comb,
253 char warn_buf[STRLEN];
256 for (i = 1; (i < eNBF_NR); i++)
258 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
265 *nb = strtol(nb_str, NULL, 10);
267 if ((*nb < 1) || (*nb >= eNBF_NR))
269 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s",
270 nb_str, enbf_names[1]);
271 warning_error(wi, warn_buf);
275 for (i = 1; (i < eCOMB_NR); i++)
277 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
284 *comb = strtol(comb_str, NULL, 10);
286 if ((*comb < 1) || (*comb >= eCOMB_NR))
288 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s",
289 comb_str, ecomb_names[1]);
290 warning_error(wi, warn_buf);
295 static char ** cpp_opts(const char *define, const char *include,
300 const char *cppadds[2];
301 char **cppopts = NULL;
302 const char *option[2] = { "-D", "-I" };
303 const char *nopt[2] = { "define", "include" };
307 char warn_buf[STRLEN];
310 cppadds[1] = include;
311 for (n = 0; (n < 2); n++)
318 while ((*ptr != '\0') && isspace(*ptr))
323 while ((*rptr != '\0') && !isspace(*rptr))
331 strncpy(buf, ptr, len);
332 if (strstr(ptr, option[n]) != ptr)
334 set_warning_line(wi, "mdp file", -1);
335 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
336 warning(wi, warn_buf);
340 srenew(cppopts, ++ncppopts);
341 cppopts[ncppopts-1] = gmx_strdup(buf);
349 srenew(cppopts, ++ncppopts);
350 cppopts[ncppopts-1] = NULL;
357 find_gb_bondlength(t_params *plist, int ai, int aj, real *length)
364 for (i = 0; i < F_NRE && !found; i++)
368 for (j = 0; j < plist[i].nr; j++)
370 a1 = plist[i].param[j].a[0];
371 a2 = plist[i].param[j].a[1];
373 if ( (a1 == ai && a2 == aj) || (a1 == aj && a2 == ai))
375 /* Equilibrium bond distance */
376 *length = plist[i].param[j].c[0];
389 find_gb_anglelength(t_params *plist, int ai, int ak, real *length)
391 int i, j, a1, a2, a3;
394 int status, status1, status2;
398 for (i = 0; i < F_NRE && !found; i++)
402 for (j = 0; j < plist[i].nr; j++)
404 a1 = plist[i].param[j].a[0];
405 a2 = plist[i].param[j].a[1];
406 a3 = plist[i].param[j].a[2];
408 /* We dont care what the middle atom is, but use it below */
409 if ( (a1 == ai && a3 == ak) || (a1 == ak && a3 == ai) )
411 /* Equilibrium bond distance */
412 a123 = plist[i].param[j].c[0];
413 /* Use middle atom to find reference distances r12 and r23 */
414 status1 = find_gb_bondlength(plist, a1, a2, &r12);
415 status2 = find_gb_bondlength(plist, a2, a3, &r23);
417 if (status1 == 0 && status2 == 0)
419 /* cosine theorem to get r13 */
420 *length = std::sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG)));
433 generate_gb_exclusion_interactions(t_molinfo *mi, gpp_atomtype_t atype, t_nextnb *nnb)
435 int j, n, ai, aj, ti, tj;
440 real radiusi, radiusj;
441 real gb_radiusi, gb_radiusj;
442 real param_c2, param_c4;
448 for (n = 1; n <= nnb->nrex; n++)
463 /* Put all higher-order exclusions into 1,4 list so we dont miss them */
470 for (ai = 0; ai < nnb->nr; ai++)
472 ti = at->atom[ai].type;
473 radiusi = get_atomtype_radius(ti, atype);
474 gb_radiusi = get_atomtype_gb_radius(ti, atype);
476 for (j = 0; j < nnb->nrexcl[ai][n]; j++)
478 aj = nnb->a[ai][n][j];
480 /* Only add the interactions once */
483 tj = at->atom[aj].type;
484 radiusj = get_atomtype_radius(tj, atype);
485 gb_radiusj = get_atomtype_gb_radius(tj, atype);
487 /* There is an exclusion of type "ftype" between atoms ai and aj */
491 /* Reference distance, not used for 1-4 interactions */
495 if (find_gb_bondlength(plist, ai, aj, &distance) != 0)
497 gmx_fatal(FARGS, "Cannot find bond length for atoms %d-%d", ai, aj);
501 if (find_gb_anglelength(plist, ai, aj, &distance) != 0)
503 gmx_fatal(FARGS, "Cannot find length for atoms %d-%d involved in angle", ai, aj);
510 /* Assign GB parameters */
512 param.c[0] = radiusi+radiusj;
513 /* Reference distance distance */
514 param.c[1] = distance;
515 /* Still parameter */
516 param.c[2] = param_c2;
518 param.c[3] = gb_radiusi+gb_radiusj;
520 param.c[4] = param_c4;
522 /* Add it to the parameter list */
523 add_param_to_list(&plist[ftype], ¶m);
532 static void make_atoms_sys(int nmolb, const gmx_molblock_t *molb,
533 const t_molinfo *molinfo,
537 const t_atoms *mol_atoms;
542 for (mb = 0; mb < nmolb; mb++)
544 mol_atoms = &molinfo[molb[mb].type].atoms;
546 srenew(atoms->atom, atoms->nr + molb[mb].nmol*mol_atoms->nr);
548 for (m = 0; m < molb[mb].nmol; m++)
550 for (a = 0; a < mol_atoms->nr; a++)
552 atoms->atom[atoms->nr++] = mol_atoms->atom[a];
559 static char **read_topol(const char *infile, const char *outfile,
560 const char *define, const char *include,
562 gpp_atomtype_t atype,
565 t_molinfo **intermolecular_interactions,
567 int *combination_rule,
572 gmx_molblock_t **molblock,
580 char *pline = NULL, **title = NULL;
581 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
583 char *dirstr, *dummy2;
584 int nrcopies, nmol, nmolb = 0, nscan, ncombs, ncopy;
585 double fLJ, fQQ, fPOW;
586 gmx_molblock_t *molb = NULL;
587 t_molinfo *mi0 = NULL;
590 t_nbparam **nbparam, **pair;
592 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
593 gmx_bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
594 double qt = 0, qBt = 0; /* total charge */
595 t_bond_atomtype batype;
597 int dcatt = -1, nmol_couple;
598 /* File handling variables */
601 char *tmp_line = NULL;
602 char warn_buf[STRLEN];
603 const char *floating_point_arithmetic_tip =
604 "Total charge should normally be an integer. See\n"
605 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
606 "for discussion on how close it should be to an integer.\n";
607 /* We need to open the output file before opening the input file,
608 * because cpp_open_file can change the current working directory.
612 out = gmx_fio_fopen(outfile, "w");
619 /* open input file */
620 status = cpp_open_file(infile, &handle, cpp_opts(define, include, wi));
623 gmx_fatal(FARGS, cpp_error(&handle, status));
626 /* some local variables */
627 DS_Init(&DS); /* directive stack */
628 nmol = 0; /* no molecules yet... */
629 d = d_invalid; /* first thing should be a directive */
630 nbparam = NULL; /* The temporary non-bonded matrix */
631 pair = NULL; /* The temporary pair interaction matrix */
632 block2 = NULL; /* the extra exclusions */
635 *reppow = 12.0; /* Default value for repulsion power */
637 *intermolecular_interactions = NULL;
639 /* Init the number of CMAP torsion angles and grid spacing */
640 plist[F_CMAP].grid_spacing = 0;
641 plist[F_CMAP].nc = 0;
643 bWarn_copy_A_B = bFEP;
645 batype = init_bond_atomtype();
646 /* parse the actual file */
647 bReadDefaults = FALSE;
649 bReadMolType = FALSE;
654 status = cpp_read_line(&handle, STRLEN, line);
655 done = (status == eCPP_EOF);
658 if (status != eCPP_OK)
660 gmx_fatal(FARGS, cpp_error(&handle, status));
664 fprintf(out, "%s\n", line);
667 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
669 pline = gmx_strdup(line);
671 /* Strip trailing '\' from pline, if it exists */
673 if ((sl > 0) && (pline[sl-1] == CONTINUE))
678 /* build one long line from several fragments - necessary for CMAP */
679 while (continuing(line))
681 status = cpp_read_line(&handle, STRLEN, line);
682 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
684 /* Since we depend on the '\' being present to continue to read, we copy line
685 * to a tmp string, strip the '\' from that string, and cat it to pline
687 tmp_line = gmx_strdup(line);
689 sl = strlen(tmp_line);
690 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE))
692 tmp_line[sl-1] = ' ';
695 done = (status == eCPP_EOF);
698 if (status != eCPP_OK)
700 gmx_fatal(FARGS, cpp_error(&handle, status));
704 fprintf(out, "%s\n", line);
708 srenew(pline, strlen(pline)+strlen(tmp_line)+1);
709 strcat(pline, tmp_line);
713 /* skip trailing and leading spaces and comment text */
714 strip_comment (pline);
717 /* if there is something left... */
718 if ((int)strlen(pline) > 0)
720 if (pline[0] == OPENDIR)
722 /* A directive on this line: copy the directive
723 * without the brackets into dirstr, then
724 * skip spaces and tabs on either side of directive
726 dirstr = gmx_strdup((pline+1));
727 if ((dummy2 = strchr (dirstr, CLOSEDIR)) != NULL)
733 if ((newd = str2dir(dirstr)) == d_invalid)
735 sprintf(errbuf, "Invalid directive %s", dirstr);
736 warning_error(wi, errbuf);
740 /* Directive found */
743 fprintf(debug, "found directive '%s'\n", dir2str(newd));
745 if (DS_Check_Order (DS, newd))
752 /* we should print here which directives should have
753 been present, and which actually are */
754 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
755 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
759 if (d == d_intermolecular_interactions)
761 if (*intermolecular_interactions == NULL)
763 /* We (mis)use the moleculetype processing
764 * to process the intermolecular interactions
765 * by making a "molecule" of the size of the system.
767 snew(*intermolecular_interactions, 1);
768 init_molinfo(*intermolecular_interactions);
769 mi0 = *intermolecular_interactions;
770 make_atoms_sys(nmolb, molb, *molinfo,
777 else if (d != d_invalid)
779 /* Not a directive, just a plain string
780 * use a gigantic switch to decode,
781 * if there is a valid directive!
788 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
789 cpp_error(&handle, eCPP_SYNTAX));
791 bReadDefaults = TRUE;
792 nscan = sscanf(pline, "%s%s%s%lf%lf%lf",
793 nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
804 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
807 bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
808 if (nb_funct != eNBF_LJ && bGenPairs)
810 gmx_fatal(FARGS, "Generating pair parameters is only supported with LJ non-bonded interactions");
826 nb_funct = ifunc_index(d_nonbond_params, nb_funct);
830 push_at(symtab, atype, batype, pline, nb_funct,
831 &nbparam, bGenPairs ? &pair : NULL, wi);
835 push_bt(d, plist, 2, NULL, batype, pline, wi);
837 case d_constrainttypes:
838 push_bt(d, plist, 2, NULL, batype, pline, wi);
843 push_nbt(d, pair, atype, pline, F_LJ14, wi);
847 push_bt(d, plist, 2, atype, NULL, pline, wi);
851 push_bt(d, plist, 3, NULL, batype, pline, wi);
853 case d_dihedraltypes:
854 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
855 push_dihedraltype(d, plist, batype, pline, wi);
858 case d_nonbond_params:
859 push_nbt(d, nbparam, atype, pline, nb_funct, wi);
864 srenew(block,nblock);
865 srenew(blockinfo,nblock);
866 blk0=&(block[nblock-1]);
867 bi0=&(blockinfo[nblock-1]);
870 push_molt(symtab,bi0,pline);
874 case d_implicit_genborn_params:
875 push_gb_params(atype, pline, wi);
878 case d_implicit_surface_params:
879 gmx_fatal(FARGS, "Implicit surface directive not supported yet.");
883 push_cmaptype(d, plist, 5, atype, batype, pline, wi);
891 if (opts->couple_moltype != NULL &&
892 (opts->couple_lam0 == ecouplamNONE ||
893 opts->couple_lam0 == ecouplamQ ||
894 opts->couple_lam1 == ecouplamNONE ||
895 opts->couple_lam1 == ecouplamQ))
897 dcatt = add_atomtype_decoupled(symtab, atype,
898 &nbparam, bGenPairs ? &pair : NULL);
900 ntype = get_atomtype_ntypes(atype);
901 ncombs = (ntype*(ntype+1))/2;
902 generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi);
903 ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
905 fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
906 free_nbparam(nbparam, ntype);
909 gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule);
910 ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
912 fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
913 free_nbparam(pair, ntype);
915 /* Copy GBSA parameters to atomtype array? */
920 push_molt(symtab, &nmol, molinfo, pline, wi);
921 srenew(block2, nmol);
922 block2[nmol-1].nr = 0;
923 mi0 = &((*molinfo)[nmol-1]);
927 push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi);
931 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
932 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
935 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE,
936 FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
946 case d_position_restraints:
947 case d_angle_restraints:
948 case d_angle_restraints_z:
949 case d_distance_restraints:
950 case d_orientation_restraints:
951 case d_dihedral_restraints:
954 case d_water_polarization:
955 case d_thole_polarization:
956 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE,
957 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
960 push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi);
964 push_vsitesn(d, mi0->plist, &(mi0->atoms), pline, wi);
967 GMX_ASSERT(block2, "block2 must always be allocated so exclusions can be processed");
968 if (!block2[nmol-1].nr)
970 init_block2(&(block2[nmol-1]), mi0->atoms.nr);
972 push_excl(pline, &(block2[nmol-1]));
976 title = put_symtab(symtab, pline);
983 push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
984 mi0 = &((*molinfo)[whichmol]);
985 srenew(molb, nmolb+1);
986 molb[nmolb].type = whichmol;
987 molb[nmolb].nmol = nrcopies;
990 bCouple = (opts->couple_moltype != NULL &&
991 (gmx_strcasecmp("system", opts->couple_moltype) == 0 ||
992 gmx_strcasecmp(*(mi0->name), opts->couple_moltype) == 0));
995 nmol_couple += nrcopies;
998 if (mi0->atoms.nr == 0)
1000 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms",
1004 "Excluding %d bonded neighbours molecule type '%s'\n",
1005 mi0->nrexcl, *mi0->name);
1006 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
1007 if (!mi0->bProcessed)
1010 generate_excl(mi0->nrexcl,
1015 merge_excl(&(mi0->excls), &(block2[whichmol]));
1016 done_block2(&(block2[whichmol]));
1017 make_shake(mi0->plist, &mi0->atoms, opts->nshake);
1021 /* nnb contains information about first,2nd,3rd bonded neighbors.
1022 * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
1024 if (bGenborn == TRUE)
1026 generate_gb_exclusion_interactions(mi0, atype, &nnb);
1033 convert_moltype_couple(mi0, dcatt, *fudgeQQ,
1034 opts->couple_lam0, opts->couple_lam1,
1036 nb_funct, &(plist[nb_funct]));
1038 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
1039 mi0->bProcessed = TRUE;
1044 fprintf (stderr, "case: %d\n", (int)d);
1045 gmx_incons("unknown directive");
1054 status = cpp_close_file(&handle);
1055 if (status != eCPP_OK)
1057 gmx_fatal(FARGS, cpp_error(&handle, status));
1062 gmx_fio_fclose(out);
1065 if (opts->couple_moltype)
1067 if (nmol_couple == 0)
1069 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling",
1070 opts->couple_moltype);
1072 fprintf(stderr, "Coupling %d copies of molecule type '%s'\n",
1073 nmol_couple, opts->couple_moltype);
1076 /* this is not very clean, but fixes core dump on empty system name */
1079 title = put_symtab(symtab, "");
1081 if (fabs(qt) > 1e-4)
1083 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
1084 warning_note(wi, warn_buf);
1086 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
1088 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
1089 warning_note(wi, warn_buf);
1092 for (i = 0; i < nmol; i++)
1094 done_block2(&(block2[i]));
1098 done_bond_atomtype(&batype);
1100 if (*intermolecular_interactions != NULL)
1102 sfree(mi0->atoms.atom);
1113 char **do_top(gmx_bool bVerbose,
1114 const char *topfile,
1115 const char *topppfile,
1120 int *combination_rule,
1121 double *repulsion_power,
1123 gpp_atomtype_t atype,
1125 t_molinfo **molinfo,
1126 t_molinfo **intermolecular_interactions,
1129 gmx_molblock_t **molblock,
1133 /* Tmpfile might contain a long path */
1134 const char *tmpfile;
1139 tmpfile = topppfile;
1148 printf("processing topology...\n");
1150 title = read_topol(topfile, tmpfile, opts->define, opts->include,
1152 nrmols, molinfo, intermolecular_interactions,
1153 plist, combination_rule, repulsion_power,
1154 opts, fudgeQQ, nmolblock, molblock,
1155 ir->efep != efepNO, bGenborn, bZero, wi);
1156 if ((*combination_rule != eCOMB_GEOMETRIC) &&
1157 (ir->vdwtype == evdwUSER))
1159 warning(wi, "Using sigma/epsilon based combination rules with"
1160 " user supplied potential function may produce unwanted"
1168 static void generate_qmexcl_moltype(gmx_moltype_t *molt, unsigned char *grpnr,
1171 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1173 /* generates the exclusions between the individual QM atoms, as
1174 * these interactions should be handled by the QM subroutines and
1175 * not by the gromacs routines
1178 i, j, l, k = 0, jmax, qm_max = 0, qm_nr = 0, nratoms = 0, link_nr = 0, link_max = 0;
1180 *qm_arr = NULL, *link_arr = NULL, a1, a2, a3, a4, ftype = 0;
1186 *bQMMM, *blink, bexcl;
1188 /* First we search and select the QM atoms in an qm_arr array that
1189 * we use to create the exclusions.
1191 * we take the possibility into account that a user has defined more
1192 * than one QM group:
1194 * for that we also need to do this an ugly work-about just in case
1195 * the QM group contains the entire system...
1197 jmax = ir->opts.ngQM;
1199 /* we first search for all the QM atoms and put them in an array
1201 for (j = 0; j < jmax; j++)
1203 for (i = 0; i < molt->atoms.nr; i++)
1205 if (qm_nr >= qm_max)
1208 srenew(qm_arr, qm_max);
1210 if ((grpnr ? grpnr[i] : 0) == j)
1212 qm_arr[qm_nr++] = i;
1216 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1217 * QM/not QM. We first set all elements to false. Afterwards we use
1218 * the qm_arr to change the elements corresponding to the QM atoms
1221 snew(bQMMM, molt->atoms.nr);
1222 for (i = 0; i < molt->atoms.nr; i++)
1226 for (i = 0; i < qm_nr; i++)
1228 bQMMM[qm_arr[i]] = TRUE;
1231 /* We remove all bonded interactions (i.e. bonds,
1232 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1233 * are removed is as follows: if the interaction invloves 2 atoms,
1234 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1235 * it is removed if at least two of the atoms are QM atoms, if the
1236 * interaction involves 4 atoms, it is removed if there are at least
1237 * 2 QM atoms. Since this routine is called once before any forces
1238 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1239 * be rewritten at this poitn without any problem. 25-9-2002 */
1241 /* first check weter we already have CONNBONDS: */
1242 if (molt->ilist[F_CONNBONDS].nr != 0)
1244 fprintf(stderr, "nr. of CONNBONDS present already: %d\n",
1245 molt->ilist[F_CONNBONDS].nr/3);
1246 ftype = molt->ilist[F_CONNBONDS].iatoms[0];
1247 k = molt->ilist[F_CONNBONDS].nr;
1249 /* now we delete all bonded interactions, except the ones describing
1250 * a chemical bond. These are converted to CONNBONDS
1252 for (i = 0; i < F_LJ; i++)
1254 if (i == F_CONNBONDS)
1258 nratoms = interaction_function[i].nratoms;
1260 while (j < molt->ilist[i].nr)
1265 a1 = molt->ilist[i].iatoms[j+1];
1266 a2 = molt->ilist[i].iatoms[j+2];
1267 bexcl = (bQMMM[a1] && bQMMM[a2]);
1268 /* a bonded beteen two QM atoms will be copied to the
1269 * CONNBONDS list, for reasons mentioned above
1271 if (bexcl && i < F_ANGLES)
1273 srenew(molt->ilist[F_CONNBONDS].iatoms, k+3);
1274 molt->ilist[F_CONNBONDS].nr += 3;
1275 molt->ilist[F_CONNBONDS].iatoms[k++] = ftype;
1276 molt->ilist[F_CONNBONDS].iatoms[k++] = a1;
1277 molt->ilist[F_CONNBONDS].iatoms[k++] = a2;
1281 a1 = molt->ilist[i].iatoms[j+1];
1282 a2 = molt->ilist[i].iatoms[j+2];
1283 a3 = molt->ilist[i].iatoms[j+3];
1284 bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1285 (bQMMM[a1] && bQMMM[a3]) ||
1286 (bQMMM[a2] && bQMMM[a3]));
1289 a1 = molt->ilist[i].iatoms[j+1];
1290 a2 = molt->ilist[i].iatoms[j+2];
1291 a3 = molt->ilist[i].iatoms[j+3];
1292 a4 = molt->ilist[i].iatoms[j+4];
1293 bexcl = ((bQMMM[a1] && bQMMM[a2] && bQMMM[a3]) ||
1294 (bQMMM[a1] && bQMMM[a2] && bQMMM[a4]) ||
1295 (bQMMM[a1] && bQMMM[a3] && bQMMM[a4]) ||
1296 (bQMMM[a2] && bQMMM[a3] && bQMMM[a4]));
1299 gmx_fatal(FARGS, "no such bonded interactions with %d atoms\n", nratoms);
1304 /* since the interaction involves QM atoms, these should be
1305 * removed from the MM ilist
1307 molt->ilist[i].nr -= (nratoms+1);
1308 for (l = j; l < molt->ilist[i].nr; l++)
1310 molt->ilist[i].iatoms[l] = molt->ilist[i].iatoms[l+(nratoms+1)];
1315 j += nratoms+1; /* the +1 is for the functype */
1319 /* Now, we search for atoms bonded to a QM atom because we also want
1320 * to exclude their nonbonded interactions with the QM atoms. The
1321 * reason for this is that this interaction is accounted for in the
1322 * linkatoms interaction with the QMatoms and would be counted
1325 for (i = 0; i < F_NRE; i++)
1330 while (j < molt->ilist[i].nr)
1332 a1 = molt->ilist[i].iatoms[j+1];
1333 a2 = molt->ilist[i].iatoms[j+2];
1334 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1336 if (link_nr >= link_max)
1339 srenew(link_arr, link_max);
1343 link_arr[link_nr++] = a2;
1347 link_arr[link_nr++] = a1;
1354 snew(blink, molt->atoms.nr);
1355 for (i = 0; i < molt->atoms.nr; i++)
1359 for (i = 0; i < link_nr; i++)
1361 blink[link_arr[i]] = TRUE;
1363 /* creating the exclusion block for the QM atoms. Each QM atom has
1364 * as excluded elements all the other QMatoms (and itself).
1366 qmexcl.nr = molt->atoms.nr;
1367 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1368 snew(qmexcl.index, qmexcl.nr+1);
1369 snew(qmexcl.a, qmexcl.nra);
1371 for (i = 0; i < qmexcl.nr; i++)
1373 qmexcl.index[i] = j;
1376 for (k = 0; k < qm_nr; k++)
1378 qmexcl.a[k+j] = qm_arr[k];
1380 for (k = 0; k < link_nr; k++)
1382 qmexcl.a[qm_nr+k+j] = link_arr[k];
1384 j += (qm_nr+link_nr);
1388 for (k = 0; k < qm_nr; k++)
1390 qmexcl.a[k+j] = qm_arr[k];
1395 qmexcl.index[qmexcl.nr] = j;
1397 /* and merging with the exclusions already present in sys.
1400 init_block2(&qmexcl2, molt->atoms.nr);
1401 b_to_b2(&qmexcl, &qmexcl2);
1402 merge_excl(&(molt->excls), &qmexcl2);
1403 done_block2(&qmexcl2);
1405 /* Finally, we also need to get rid of the pair interactions of the
1406 * classical atom bonded to the boundary QM atoms with the QMatoms,
1407 * as this interaction is already accounted for by the QM, so also
1408 * here we run the risk of double counting! We proceed in a similar
1409 * way as we did above for the other bonded interactions: */
1410 for (i = F_LJ14; i < F_COUL14; i++)
1412 nratoms = interaction_function[i].nratoms;
1414 while (j < molt->ilist[i].nr)
1416 a1 = molt->ilist[i].iatoms[j+1];
1417 a2 = molt->ilist[i].iatoms[j+2];
1418 bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1419 (blink[a1] && bQMMM[a2]) ||
1420 (bQMMM[a1] && blink[a2]));
1423 /* since the interaction involves QM atoms, these should be
1424 * removed from the MM ilist
1426 molt->ilist[i].nr -= (nratoms+1);
1427 for (k = j; k < molt->ilist[i].nr; k++)
1429 molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1434 j += nratoms+1; /* the +1 is for the functype */
1443 } /* generate_qmexcl */
1445 void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi)
1447 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1450 unsigned char *grpnr;
1451 int mb, mol, nat_mol, i, nr_mol_with_qm_atoms = 0;
1452 gmx_molblock_t *molb;
1455 grpnr = sys->groups.grpnr[egcQMMM];
1457 for (mb = 0; mb < sys->nmolblock; mb++)
1459 molb = &sys->molblock[mb];
1460 nat_mol = sys->moltype[molb->type].atoms.nr;
1461 for (mol = 0; mol < molb->nmol; mol++)
1464 for (i = 0; i < nat_mol; i++)
1466 if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM)
1473 nr_mol_with_qm_atoms++;
1476 /* We need to split this molblock */
1479 /* Split the molblock at this molecule */
1481 srenew(sys->molblock, sys->nmolblock);
1482 for (i = sys->nmolblock-2; i >= mb; i--)
1484 sys->molblock[i+1] = sys->molblock[i];
1486 sys->molblock[mb ].nmol = mol;
1487 sys->molblock[mb+1].nmol -= mol;
1489 molb = &sys->molblock[mb];
1493 /* Split the molblock after this molecule */
1495 srenew(sys->molblock, sys->nmolblock);
1496 molb = &sys->molblock[mb];
1497 for (i = sys->nmolblock-2; i >= mb; i--)
1499 sys->molblock[i+1] = sys->molblock[i];
1501 sys->molblock[mb ].nmol = 1;
1502 sys->molblock[mb+1].nmol -= 1;
1505 /* Add a moltype for the QMMM molecule */
1507 srenew(sys->moltype, sys->nmoltype);
1508 /* Copy the moltype struct */
1509 sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
1510 /* Copy the exclusions to a new array, since this is the only
1511 * thing that needs to be modified for QMMM.
1513 copy_blocka(&sys->moltype[molb->type ].excls,
1514 &sys->moltype[sys->nmoltype-1].excls);
1515 /* Set the molecule type for the QMMM molblock */
1516 molb->type = sys->nmoltype - 1;
1518 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir);
1526 if (nr_mol_with_qm_atoms > 1)
1528 /* generate a warning is there are QM atoms in different
1529 * topologies. In this case it is not possible at this stage to
1530 * mutualy exclude the non-bonded interactions via the
1531 * exclusions (AFAIK). Instead, the user is advised to use the
1532 * energy group exclusions in the mdp file
1535 "\nThe QM subsystem is divided over multiple topologies. "
1536 "The mutual non-bonded interactions cannot be excluded. "
1537 "There are two ways to achieve this:\n\n"
1538 "1) merge the topologies, such that the atoms of the QM "
1539 "subsystem are all present in one single topology file. "
1540 "In this case this warning will dissappear\n\n"
1541 "2) exclude the non-bonded interactions explicitly via the "
1542 "energygrp-excl option in the mdp file. if this is the case "
1543 "this warning may be ignored"