1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
11 *
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
16 *
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
21 *
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 *
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
34 *
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
37 */
38 #include "gmxpre.h"
39
40 #include "topio.h"
41
42 #include <cassert>
43 #include <cctype>
44 #include <cerrno>
45 #include <cmath>
46 #include <cstdio>
47 #include <cstdlib>
48 #include <cstring>
49
50 #include <algorithm>
51 #include <memory>
52 #include <unordered_set>
53
54 #include <sys/types.h>
55
56 #include "gromacs/fileio/gmxfio.h"
57 #include "gromacs/fileio/warninp.h"
58 #include "gromacs/gmxpreprocess/gmxcpp.h"
59 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
60 #include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
61 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
62 #include "gromacs/gmxpreprocess/grompp_impl.h"
63 #include "gromacs/gmxpreprocess/readir.h"
64 #include "gromacs/gmxpreprocess/topdirs.h"
65 #include "gromacs/gmxpreprocess/toppush.h"
66 #include "gromacs/gmxpreprocess/topshake.h"
67 #include "gromacs/gmxpreprocess/toputil.h"
68 #include "gromacs/gmxpreprocess/vsite_parm.h"
69 #include "gromacs/math/units.h"
70 #include "gromacs/math/utilities.h"
71 #include "gromacs/mdtypes/inputrec.h"
72 #include "gromacs/mdtypes/md_enums.h"
73 #include "gromacs/pbcutil/pbc.h"
74 #include "gromacs/topology/block.h"
75 #include "gromacs/topology/exclusionblocks.h"
76 #include "gromacs/topology/ifunc.h"
77 #include "gromacs/topology/symtab.h"
78 #include "gromacs/topology/topology.h"
79 #include "gromacs/utility/cstringutil.h"
80 #include "gromacs/utility/fatalerror.h"
81 #include "gromacs/utility/futil.h"
82 #include "gromacs/utility/gmxassert.h"
83 #include "gromacs/utility/logger.h"
84 #include "gromacs/utility/pleasecite.h"
85 #include "gromacs/utility/smalloc.h"
86
87 #define OPENDIR '[' /* starting sign for directive */
88 #define CLOSEDIR ']' /* ending sign for directive */
89
gen_pairs(const InteractionsOfType & nbs,InteractionsOfType * pairs,real fudge,int comb)90 static void gen_pairs(const InteractionsOfType& nbs, InteractionsOfType* pairs, real fudge, int comb)
91 {
92 real scaling;
93 int ntp = nbs.size();
94 int nnn = static_cast<int>(std::sqrt(static_cast<double>(ntp)));
95 GMX_ASSERT(nnn * nnn == ntp,
96 "Number of pairs of generated non-bonded parameters should be a perfect square");
97 int nrfp = NRFP(F_LJ);
98 int nrfpA = interaction_function[F_LJ14].nrfpA;
99 int nrfpB = interaction_function[F_LJ14].nrfpB;
100
101 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
102 {
103 gmx_incons("Number of force parameters in gen_pairs wrong");
104 }
105
106 fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
107 pairs->interactionTypes.clear();
108 int i = 0;
109 std::array<int, 2> atomNumbers;
110 std::array<real, MAXFORCEPARAM> forceParam = { NOTSET };
111 for (const auto& type : nbs.interactionTypes)
112 {
113 /* Copy type.atoms */
114 atomNumbers = { i / nnn, i % nnn };
115 /* Copy normal and FEP parameters and multiply by fudge factor */
116 gmx::ArrayRef<const real> existingParam = type.forceParam();
117 GMX_RELEASE_ASSERT(2 * nrfp <= MAXFORCEPARAM,
118 "Can't have more parameters than half of maximum p arameter number");
119 for (int j = 0; j < nrfp; j++)
120 {
121 /* If we are using sigma/epsilon values, only the epsilon values
122 * should be scaled, but not sigma.
123 * The sigma values have even indices 0,2, etc.
124 */
125 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j % 2 == 0))
126 {
127 scaling = 1.0;
128 }
129 else
130 {
131 scaling = fudge;
132 }
133
134 forceParam[j] = scaling * existingParam[j];
135 forceParam[nrfp + j] = scaling * existingParam[j];
136 }
137 pairs->interactionTypes.emplace_back(InteractionOfType(atomNumbers, forceParam));
138 i++;
139 }
140 }
141
check_mol(const gmx_mtop_t * mtop,warninp * wi)142 double check_mol(const gmx_mtop_t* mtop, warninp* wi)
143 {
144 char buf[256];
145 int i, ri, pt;
146 double q;
147 real m, mB;
148
149 /* Check mass and charge */
150 q = 0.0;
151
152 for (const gmx_molblock_t& molb : mtop->molblock)
153 {
154 const t_atoms* atoms = &mtop->moltype[molb.type].atoms;
155 for (i = 0; (i < atoms->nr); i++)
156 {
157 q += molb.nmol * atoms->atom[i].q;
158 m = atoms->atom[i].m;
159 mB = atoms->atom[i].mB;
160 pt = atoms->atom[i].ptype;
161 /* If the particle is an atom or a nucleus it must have a mass,
162 * else, if it is a shell, a vsite or a bondshell it can have mass zero
163 */
164 if (((m <= 0.0) || (mB <= 0.0)) && ((pt == eptAtom) || (pt == eptNucleus)))
165 {
166 ri = atoms->atom[i].resind;
167 sprintf(buf, "atom %s (Res %s-%d) has mass %g (state A) / %g (state B)\n",
168 *(atoms->atomname[i]), *(atoms->resinfo[ri].name), atoms->resinfo[ri].nr, m, mB);
169 warning_error(wi, buf);
170 }
171 else if (((m != 0) || (mB != 0)) && (pt == eptVSite))
172 {
173 ri = atoms->atom[i].resind;
174 sprintf(buf,
175 "virtual site %s (Res %s-%d) has non-zero mass %g (state A) / %g (state "
176 "B)\n"
177 " Check your topology.\n",
178 *(atoms->atomname[i]), *(atoms->resinfo[ri].name), atoms->resinfo[ri].nr, m, mB);
179 warning_error(wi, buf);
180 /* The following statements make LINCS break! */
181 /* atoms->atom[i].m=0; */
182 }
183 }
184 }
185 return q;
186 }
187
188 /*! \brief Returns the rounded charge of a molecule, when close to integer, otherwise returns the original charge.
189 *
190 * The results of this routine are only used for checking and for
191 * printing warning messages. Thus we can assume that charges of molecules
192 * should be integer. If the user wanted non-integer molecular charge,
193 * an undesired warning is printed and the user should use grompp -maxwarn 1.
194 *
195 * \param qMol The total, unrounded, charge of the molecule
196 * \param sumAbsQ The sum of absolute values of the charges, used for determining the tolerance for the rounding.
197 */
roundedMoleculeCharge(double qMol,double sumAbsQ)198 static double roundedMoleculeCharge(double qMol, double sumAbsQ)
199 {
200 /* We use a tolerance of 1e-6 for inaccuracies beyond the 6th decimal
201 * of the charges for ascii float truncation in the topology files.
202 * Although the summation here uses double precision, the charges
203 * are read and stored in single precision when real=float. This can
204 * lead to rounding errors of half the least significant bit.
205 * Note that, unfortunately, we can not assume addition of random
206 * rounding errors. It is not entirely unlikely that many charges
207 * have a near half-bit rounding error with the same sign.
208 */
209 double tolAbs = 1e-6;
210 double tol = std::max(tolAbs, 0.5 * GMX_REAL_EPS * sumAbsQ);
211 double qRound = std::round(qMol);
212 if (std::abs(qMol - qRound) <= tol)
213 {
214 return qRound;
215 }
216 else
217 {
218 return qMol;
219 }
220 }
221
sum_q(const t_atoms * atoms,int numMols,double * qTotA,double * qTotB)222 static void sum_q(const t_atoms* atoms, int numMols, double* qTotA, double* qTotB)
223 {
224 /* sum charge */
225 double qmolA = 0;
226 double qmolB = 0;
227 double sumAbsQA = 0;
228 double sumAbsQB = 0;
229 for (int i = 0; i < atoms->nr; i++)
230 {
231 qmolA += atoms->atom[i].q;
232 qmolB += atoms->atom[i].qB;
233 sumAbsQA += std::abs(atoms->atom[i].q);
234 sumAbsQB += std::abs(atoms->atom[i].qB);
235 }
236
237 *qTotA += numMols * roundedMoleculeCharge(qmolA, sumAbsQA);
238 *qTotB += numMols * roundedMoleculeCharge(qmolB, sumAbsQB);
239 }
240
get_nbparm(char * nb_str,char * comb_str,int * nb,int * comb,warninp * wi)241 static void get_nbparm(char* nb_str, char* comb_str, int* nb, int* comb, warninp* wi)
242 {
243 int i;
244 char warn_buf[STRLEN];
245
246 *nb = -1;
247 for (i = 1; (i < eNBF_NR); i++)
248 {
249 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
250 {
251 *nb = i;
252 }
253 }
254 if (*nb == -1)
255 {
256 *nb = strtol(nb_str, nullptr, 10);
257 }
258 if ((*nb < 1) || (*nb >= eNBF_NR))
259 {
260 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s", nb_str, enbf_names[1]);
261 warning_error(wi, warn_buf);
262 *nb = 1;
263 }
264 *comb = -1;
265 for (i = 1; (i < eCOMB_NR); i++)
266 {
267 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
268 {
269 *comb = i;
270 }
271 }
272 if (*comb == -1)
273 {
274 *comb = strtol(comb_str, nullptr, 10);
275 }
276 if ((*comb < 1) || (*comb >= eCOMB_NR))
277 {
278 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s", comb_str, ecomb_names[1]);
279 warning_error(wi, warn_buf);
280 *comb = 1;
281 }
282 }
283
cpp_opts(const char * define,const char * include,warninp * wi)284 static char** cpp_opts(const char* define, const char* include, warninp* wi)
285 {
286 int n, len;
287 int ncppopts = 0;
288 const char* cppadds[2];
289 char** cppopts = nullptr;
290 const char* option[2] = { "-D", "-I" };
291 const char* nopt[2] = { "define", "include" };
292 const char* ptr;
293 const char* rptr;
294 char* buf;
295 char warn_buf[STRLEN];
296
297 cppadds[0] = define;
298 cppadds[1] = include;
299 for (n = 0; (n < 2); n++)
300 {
301 if (cppadds[n])
302 {
303 ptr = cppadds[n];
304 while (*ptr != '\0')
305 {
306 while ((*ptr != '\0') && isspace(*ptr))
307 {
308 ptr++;
309 }
310 rptr = ptr;
311 while ((*rptr != '\0') && !isspace(*rptr))
312 {
313 rptr++;
314 }
315 len = (rptr - ptr);
316 if (len > 2)
317 {
318 snew(buf, (len + 1));
319 strncpy(buf, ptr, len);
320 if (strstr(ptr, option[n]) != ptr)
321 {
322 set_warning_line(wi, "mdp file", -1);
323 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
324 warning(wi, warn_buf);
325 }
326 else
327 {
328 srenew(cppopts, ++ncppopts);
329 cppopts[ncppopts - 1] = gmx_strdup(buf);
330 }
331 sfree(buf);
332 ptr = rptr;
333 }
334 }
335 }
336 }
337 srenew(cppopts, ++ncppopts);
338 cppopts[ncppopts - 1] = nullptr;
339
340 return cppopts;
341 }
342
343
make_atoms_sys(gmx::ArrayRef<const gmx_molblock_t> molblock,gmx::ArrayRef<const MoleculeInformation> molinfo,t_atoms * atoms)344 static void make_atoms_sys(gmx::ArrayRef<const gmx_molblock_t> molblock,
345 gmx::ArrayRef<const MoleculeInformation> molinfo,
346 t_atoms* atoms)
347 {
348 atoms->nr = 0;
349 atoms->atom = nullptr;
350
351 for (const gmx_molblock_t& molb : molblock)
352 {
353 const t_atoms& mol_atoms = molinfo[molb.type].atoms;
354
355 srenew(atoms->atom, atoms->nr + molb.nmol * mol_atoms.nr);
356
357 for (int m = 0; m < molb.nmol; m++)
358 {
359 for (int a = 0; a < mol_atoms.nr; a++)
360 {
361 atoms->atom[atoms->nr++] = mol_atoms.atom[a];
362 }
363 }
364 }
365 }
366
367
read_topol(const char * infile,const char * outfile,const char * define,const char * include,t_symtab * symtab,PreprocessingAtomTypes * atypes,std::vector<MoleculeInformation> * molinfo,std::unique_ptr<MoleculeInformation> * intermolecular_interactions,gmx::ArrayRef<InteractionsOfType> interactions,int * combination_rule,double * reppow,t_gromppopts * opts,real * fudgeQQ,std::vector<gmx_molblock_t> * molblock,bool * ffParametrizedWithHBondConstraints,bool bFEP,bool bZero,bool usingFullRangeElectrostatics,warninp * wi,const gmx::MDLogger & logger)368 static char** read_topol(const char* infile,
369 const char* outfile,
370 const char* define,
371 const char* include,
372 t_symtab* symtab,
373 PreprocessingAtomTypes* atypes,
374 std::vector<MoleculeInformation>* molinfo,
375 std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
376 gmx::ArrayRef<InteractionsOfType> interactions,
377 int* combination_rule,
378 double* reppow,
379 t_gromppopts* opts,
380 real* fudgeQQ,
381 std::vector<gmx_molblock_t>* molblock,
382 bool* ffParametrizedWithHBondConstraints,
383 bool bFEP,
384 bool bZero,
385 bool usingFullRangeElectrostatics,
386 warninp* wi,
387 const gmx::MDLogger& logger)
388 {
389 FILE* out;
390 int sl, nb_funct;
391 char * pline = nullptr, **title = nullptr;
392 char line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
393 char genpairs[32];
394 char * dirstr, *dummy2;
395 int nrcopies, nscan, ncombs, ncopy;
396 double fLJ, fQQ, fPOW;
397 MoleculeInformation* mi0 = nullptr;
398 DirStack* DS;
399 Directive d, newd;
400 t_nbparam ** nbparam, **pair;
401 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
402 bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
403 double qt = 0, qBt = 0; /* total charge */
404 int dcatt = -1, nmol_couple;
405 /* File handling variables */
406 int status;
407 bool done;
408 gmx_cpp_t handle;
409 char* tmp_line = nullptr;
410 char warn_buf[STRLEN];
411 const char* floating_point_arithmetic_tip =
412 "Total charge should normally be an integer. See\n"
413 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
414 "for discussion on how close it should be to an integer.\n";
415 /* We need to open the output file before opening the input file,
416 * because cpp_open_file can change the current working directory.
417 */
418 if (outfile)
419 {
420 out = gmx_fio_fopen(outfile, "w");
421 }
422 else
423 {
424 out = nullptr;
425 }
426
427 /* open input file */
428 auto cpp_opts_return = cpp_opts(define, include, wi);
429 status = cpp_open_file(infile, &handle, cpp_opts_return);
430 if (status != 0)
431 {
432 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
433 }
434
435 /* some local variables */
436 DS_Init(&DS); /* directive stack */
437 d = Directive::d_invalid; /* first thing should be a directive */
438 nbparam = nullptr; /* The temporary non-bonded matrix */
439 pair = nullptr; /* The temporary pair interaction matrix */
440 std::vector<std::vector<gmx::ExclusionBlock>> exclusionBlocks;
441 nb_funct = F_LJ;
442
443 *reppow = 12.0; /* Default value for repulsion power */
444
445 /* Init the number of CMAP torsion angles and grid spacing */
446 interactions[F_CMAP].cmakeGridSpacing = 0;
447 interactions[F_CMAP].cmapAngles = 0;
448
449 bWarn_copy_A_B = bFEP;
450
451 PreprocessingBondAtomType bondAtomType;
452 /* parse the actual file */
453 bReadDefaults = FALSE;
454 bGenPairs = FALSE;
455 bReadMolType = FALSE;
456 nmol_couple = 0;
457
458 do
459 {
460 status = cpp_read_line(&handle, STRLEN, line);
461 done = (status == eCPP_EOF);
462 if (!done)
463 {
464 if (status != eCPP_OK)
465 {
466 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
467 }
468 else if (out)
469 {
470 fprintf(out, "%s\n", line);
471 }
472
473 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
474
475 pline = gmx_strdup(line);
476
477 /* Strip trailing '\' from pline, if it exists */
478 sl = strlen(pline);
479 if ((sl > 0) && (pline[sl - 1] == CONTINUE))
480 {
481 pline[sl - 1] = ' ';
482 }
483
484 /* build one long line from several fragments - necessary for CMAP */
485 while (continuing(line))
486 {
487 status = cpp_read_line(&handle, STRLEN, line);
488 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
489
490 /* Since we depend on the '\' being present to continue to read, we copy line
491 * to a tmp string, strip the '\' from that string, and cat it to pline
492 */
493 tmp_line = gmx_strdup(line);
494
495 sl = strlen(tmp_line);
496 if ((sl > 0) && (tmp_line[sl - 1] == CONTINUE))
497 {
498 tmp_line[sl - 1] = ' ';
499 }
500
501 done = (status == eCPP_EOF);
502 if (!done)
503 {
504 if (status != eCPP_OK)
505 {
506 gmx_fatal(FARGS, "%s", cpp_error(&handle, status));
507 }
508 else if (out)
509 {
510 fprintf(out, "%s\n", line);
511 }
512 }
513
514 srenew(pline, strlen(pline) + strlen(tmp_line) + 1);
515 strcat(pline, tmp_line);
516 sfree(tmp_line);
517 }
518
519 /* skip trailing and leading spaces and comment text */
520 strip_comment(pline);
521 trim(pline);
522
523 /* if there is something left... */
524 if (static_cast<int>(strlen(pline)) > 0)
525 {
526 if (pline[0] == OPENDIR)
527 {
528 /* A directive on this line: copy the directive
529 * without the brackets into dirstr, then
530 * skip spaces and tabs on either side of directive
531 */
532 dirstr = gmx_strdup((pline + 1));
533 if ((dummy2 = strchr(dirstr, CLOSEDIR)) != nullptr)
534 {
535 (*dummy2) = 0;
536 }
537 trim(dirstr);
538
539 if ((newd = str2dir(dirstr)) == Directive::d_invalid)
540 {
541 sprintf(errbuf, "Invalid directive %s", dirstr);
542 warning_error(wi, errbuf);
543 }
544 else
545 {
546 /* Directive found */
547 if (DS_Check_Order(DS, newd))
548 {
549 DS_Push(&DS, newd);
550 d = newd;
551 }
552 else
553 {
554 /* we should print here which directives should have
555 been present, and which actually are */
556 gmx_fatal(FARGS, "%s\nInvalid order for directive %s",
557 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
558 /* d = Directive::d_invalid; */
559 }
560
561 if (d == Directive::d_intermolecular_interactions)
562 {
563 if (*intermolecular_interactions == nullptr)
564 {
565 /* We (mis)use the moleculetype processing
566 * to process the intermolecular interactions
567 * by making a "molecule" of the size of the system.
568 */
569 *intermolecular_interactions = std::make_unique<MoleculeInformation>();
570 mi0 = intermolecular_interactions->get();
571 mi0->initMolInfo();
572 make_atoms_sys(*molblock, *molinfo, &mi0->atoms);
573 }
574 }
575 }
576 sfree(dirstr);
577 }
578 else if (d != Directive::d_invalid)
579 {
580 /* Not a directive, just a plain string
581 * use a gigantic switch to decode,
582 * if there is a valid directive!
583 */
584 switch (d)
585 {
586 case Directive::d_defaults:
587 if (bReadDefaults)
588 {
589 gmx_fatal(FARGS, "%s\nFound a second defaults directive.\n",
590 cpp_error(&handle, eCPP_SYNTAX));
591 }
592 bReadDefaults = TRUE;
593 nscan = sscanf(pline, "%s%s%s%lf%lf%lf", nb_str, comb_str, genpairs,
594 &fLJ, &fQQ, &fPOW);
595 if (nscan < 2)
596 {
597 too_few(wi);
598 }
599 else
600 {
601 bGenPairs = FALSE;
602 fudgeLJ = 1.0;
603 *fudgeQQ = 1.0;
604
605 get_nbparm(nb_str, comb_str, &nb_funct, combination_rule, wi);
606 if (nscan >= 3)
607 {
608 bGenPairs = (gmx::equalCaseInsensitive(genpairs, "Y", 1));
609 if (nb_funct != eNBF_LJ && bGenPairs)
610 {
611 gmx_fatal(FARGS,
612 "Generating pair parameters is only supported "
613 "with LJ non-bonded interactions");
614 }
615 }
616 if (nscan >= 4)
617 {
618 fudgeLJ = fLJ;
619 }
620 if (nscan >= 5)
621 {
622 *fudgeQQ = fQQ;
623 }
624 if (nscan >= 6)
625 {
626 *reppow = fPOW;
627 }
628 }
629 nb_funct = ifunc_index(Directive::d_nonbond_params, nb_funct);
630
631 break;
632 case Directive::d_atomtypes:
633 push_at(symtab, atypes, &bondAtomType, pline, nb_funct, &nbparam,
634 bGenPairs ? &pair : nullptr, wi);
635 break;
636
637 case Directive::d_bondtypes: // Intended to fall through
638 case Directive::d_constrainttypes:
639 push_bt(d, interactions, 2, nullptr, &bondAtomType, pline, wi);
640 break;
641 case Directive::d_pairtypes:
642 if (bGenPairs)
643 {
644 push_nbt(d, pair, atypes, pline, F_LJ14, wi);
645 }
646 else
647 {
648 push_bt(d, interactions, 2, atypes, nullptr, pline, wi);
649 }
650 break;
651 case Directive::d_angletypes:
652 push_bt(d, interactions, 3, nullptr, &bondAtomType, pline, wi);
653 break;
654 case Directive::d_dihedraltypes:
655 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
656 push_dihedraltype(d, interactions, &bondAtomType, pline, wi);
657 break;
658
659 case Directive::d_nonbond_params:
660 push_nbt(d, nbparam, atypes, pline, nb_funct, wi);
661 break;
662
663 case Directive::d_implicit_genborn_params: // NOLINT bugprone-branch-clone
664 // Skip this line, so old topologies with
665 // GB parameters can be read.
666 break;
667
668 case Directive::d_implicit_surface_params:
669 // Skip this line, so that any topologies
670 // with surface parameters can be read
671 // (even though these were never formally
672 // supported).
673 break;
674
675 case Directive::d_cmaptypes:
676 push_cmaptype(d, interactions, 5, atypes, &bondAtomType, pline, wi);
677 break;
678
679 case Directive::d_moleculetype:
680 {
681 if (!bReadMolType)
682 {
683 int ntype;
684 if (opts->couple_moltype != nullptr
685 && (opts->couple_lam0 == ecouplamNONE || opts->couple_lam0 == ecouplamQ
686 || opts->couple_lam1 == ecouplamNONE
687 || opts->couple_lam1 == ecouplamQ))
688 {
689 dcatt = add_atomtype_decoupled(symtab, atypes, &nbparam,
690 bGenPairs ? &pair : nullptr);
691 }
692 ntype = atypes->size();
693 ncombs = (ntype * (ntype + 1)) / 2;
694 generate_nbparams(*combination_rule, nb_funct,
695 &(interactions[nb_funct]), atypes, wi);
696 ncopy = copy_nbparams(nbparam, nb_funct, &(interactions[nb_funct]), ntype);
697 GMX_LOG(logger.info)
698 .asParagraph()
699 .appendTextFormatted(
700 "Generated %d of the %d non-bonded parameter "
701 "combinations",
702 ncombs - ncopy, ncombs);
703 free_nbparam(nbparam, ntype);
704 if (bGenPairs)
705 {
706 gen_pairs((interactions[nb_funct]), &(interactions[F_LJ14]),
707 fudgeLJ, *combination_rule);
708 ncopy = copy_nbparams(pair, nb_funct, &(interactions[F_LJ14]), ntype);
709 GMX_LOG(logger.info)
710 .asParagraph()
711 .appendTextFormatted(
712 "Generated %d of the %d 1-4 parameter "
713 "combinations",
714 ncombs - ncopy, ncombs);
715 free_nbparam(pair, ntype);
716 }
717 /* Copy GBSA parameters to atomtype array? */
718
719 bReadMolType = TRUE;
720 }
721
722 push_molt(symtab, molinfo, pline, wi);
723 exclusionBlocks.emplace_back();
724 mi0 = &molinfo->back();
725 mi0->atoms.haveMass = TRUE;
726 mi0->atoms.haveCharge = TRUE;
727 mi0->atoms.haveType = TRUE;
728 mi0->atoms.haveBState = TRUE;
729 mi0->atoms.havePdbInfo = FALSE;
730 break;
731 }
732 case Directive::d_atoms:
733 push_atom(symtab, &(mi0->atoms), atypes, pline, wi);
734 break;
735
736 case Directive::d_pairs:
737 GMX_RELEASE_ASSERT(
738 mi0,
739 "Need to have a valid MoleculeInformation object to work on");
740 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
741 pline, FALSE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
742 break;
743 case Directive::d_pairs_nb:
744 GMX_RELEASE_ASSERT(
745 mi0,
746 "Need to have a valid MoleculeInformation object to work on");
747 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
748 pline, FALSE, FALSE, 1.0, bZero, &bWarn_copy_A_B, wi);
749 break;
750
751 case Directive::d_vsites1:
752 case Directive::d_vsites2:
753 case Directive::d_vsites3:
754 case Directive::d_vsites4:
755 case Directive::d_bonds:
756 case Directive::d_angles:
757 case Directive::d_constraints:
758 case Directive::d_settles:
759 case Directive::d_position_restraints:
760 case Directive::d_angle_restraints:
761 case Directive::d_angle_restraints_z:
762 case Directive::d_distance_restraints:
763 case Directive::d_orientation_restraints:
764 case Directive::d_dihedral_restraints:
765 case Directive::d_dihedrals:
766 case Directive::d_polarization:
767 case Directive::d_water_polarization:
768 case Directive::d_thole_polarization:
769 GMX_RELEASE_ASSERT(
770 mi0,
771 "Need to have a valid MoleculeInformation object to work on");
772 push_bond(d, interactions, mi0->interactions, &(mi0->atoms), atypes,
773 pline, TRUE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
774 break;
775 case Directive::d_cmap:
776 GMX_RELEASE_ASSERT(
777 mi0,
778 "Need to have a valid MoleculeInformation object to work on");
779 push_cmap(d, interactions, mi0->interactions, &(mi0->atoms), atypes, pline, wi);
780 break;
781
782 case Directive::d_vsitesn:
783 GMX_RELEASE_ASSERT(
784 mi0,
785 "Need to have a valid MoleculeInformation object to work on");
786 push_vsitesn(d, mi0->interactions, &(mi0->atoms), pline, wi);
787 break;
788 case Directive::d_exclusions:
789 GMX_ASSERT(!exclusionBlocks.empty(),
790 "exclusionBlocks must always be allocated so exclusions can "
791 "be processed");
792 if (exclusionBlocks.back().empty())
793 {
794 GMX_RELEASE_ASSERT(mi0,
795 "Need to have a valid MoleculeInformation "
796 "object to work on");
797 exclusionBlocks.back().resize(mi0->atoms.nr);
798 }
799 push_excl(pline, exclusionBlocks.back(), wi);
800 break;
801 case Directive::d_system:
802 trim(pline);
803 title = put_symtab(symtab, pline);
804 break;
805 case Directive::d_molecules:
806 {
807 int whichmol;
808 bool bCouple;
809
810 push_mol(*molinfo, pline, &whichmol, &nrcopies, wi);
811 mi0 = &((*molinfo)[whichmol]);
812 molblock->resize(molblock->size() + 1);
813 molblock->back().type = whichmol;
814 molblock->back().nmol = nrcopies;
815
816 bCouple = (opts->couple_moltype != nullptr
817 && (gmx_strcasecmp("system", opts->couple_moltype) == 0
818 || strcmp(*(mi0->name), opts->couple_moltype) == 0));
819 if (bCouple)
820 {
821 nmol_couple += nrcopies;
822 }
823
824 if (mi0->atoms.nr == 0)
825 {
826 gmx_fatal(FARGS, "Molecule type '%s' contains no atoms", *mi0->name);
827 }
828 GMX_LOG(logger.info)
829 .asParagraph()
830 .appendTextFormatted(
831 "Excluding %d bonded neighbours molecule type '%s'",
832 mi0->nrexcl, *mi0->name);
833 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
834 if (!mi0->bProcessed)
835 {
836 generate_excl(mi0->nrexcl, mi0->atoms.nr, mi0->interactions, &(mi0->excls));
837 gmx::mergeExclusions(&(mi0->excls), exclusionBlocks[whichmol]);
838 make_shake(mi0->interactions, &mi0->atoms, opts->nshake, logger);
839
840 if (bCouple)
841 {
842 convert_moltype_couple(mi0, dcatt, *fudgeQQ, opts->couple_lam0,
843 opts->couple_lam1, opts->bCoupleIntra,
844 nb_funct, &(interactions[nb_funct]), wi);
845 }
846 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE);
847 mi0->bProcessed = TRUE;
848 }
849 break;
850 }
851 default:
852 GMX_LOG(logger.warning)
853 .asParagraph()
854 .appendTextFormatted("case: %d", static_cast<int>(d));
855 gmx_incons("unknown directive");
856 }
857 }
858 }
859 sfree(pline);
860 pline = nullptr;
861 }
862 } while (!done);
863
864 // Check that all strings defined with -D were used when processing topology
865 std::string unusedDefineWarning = checkAndWarnForUnusedDefines(*handle);
866 if (!unusedDefineWarning.empty())
867 {
868 warning(wi, unusedDefineWarning);
869 }
870
871 sfree(cpp_opts_return);
872
873 if (out)
874 {
875 gmx_fio_fclose(out);
876 }
877
878 /* List of GROMACS define names for force fields that have been
879 * parametrized using constraints involving hydrogens only.
880 *
881 * We should avoid hardcoded names, but this is hopefully only
882 * needed temparorily for discouraging use of constraints=all-bonds.
883 */
884 const std::array<std::string, 3> ffDefines = { "_FF_AMBER", "_FF_CHARMM", "_FF_OPLSAA" };
885 *ffParametrizedWithHBondConstraints = false;
886 for (const std::string& ffDefine : ffDefines)
887 {
888 if (cpp_find_define(&handle, ffDefine))
889 {
890 *ffParametrizedWithHBondConstraints = true;
891 }
892 }
893
894 if (cpp_find_define(&handle, "_FF_GROMOS96") != nullptr)
895 {
896 warning(wi,
897 "The GROMOS force fields have been parametrized with a physically incorrect "
898 "multiple-time-stepping scheme for a twin-range cut-off. When used with "
899 "a single-range cut-off (or a correct Trotter multiple-time-stepping scheme), "
900 "physical properties, such as the density, might differ from the intended values. "
901 "Since there are researchers actively working on validating GROMOS with modern "
902 "integrators we have not yet removed the GROMOS force fields, but you should be "
903 "aware of these issues and check if molecules in your system are affected before "
904 "proceeding. "
905 "Further information is available at https://redmine.gromacs.org/issues/2884 , "
906 "and a longer explanation of our decision to remove physically incorrect "
907 "algorithms "
908 "can be found at https://doi.org/10.26434/chemrxiv.11474583.v1 .");
909 }
910 // TODO: Update URL for Issue #2884 in conjunction with updating grompp.warn in regressiontests.
911
912 cpp_done(handle);
913
914 if (opts->couple_moltype)
915 {
916 if (nmol_couple == 0)
917 {
918 gmx_fatal(FARGS, "Did not find any molecules of type '%s' for coupling", opts->couple_moltype);
919 }
920 GMX_LOG(logger.info)
921 .asParagraph()
922 .appendTextFormatted("Coupling %d copies of molecule type '%s'", nmol_couple,
923 opts->couple_moltype);
924 }
925
926 /* this is not very clean, but fixes core dump on empty system name */
927 if (!title)
928 {
929 title = put_symtab(symtab, "");
930 }
931
932 if (fabs(qt) > 1e-4)
933 {
934 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
935 warning_note(wi, warn_buf);
936 }
937 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
938 {
939 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt,
940 floating_point_arithmetic_tip);
941 warning_note(wi, warn_buf);
942 }
943 if (usingFullRangeElectrostatics && (fabs(qt) > 1e-4 || fabs(qBt) > 1e-4))
944 {
945 warning(wi,
946 "You are using Ewald electrostatics in a system with net charge. This can lead to "
947 "severe artifacts, such as ions moving into regions with low dielectric, due to "
948 "the uniform background charge. We suggest to neutralize your system with counter "
949 "ions, possibly in combination with a physiological salt concentration.");
950 please_cite(stdout, "Hub2014a");
951 }
952
953 DS_Done(&DS);
954
955 if (*intermolecular_interactions != nullptr)
956 {
957 sfree(intermolecular_interactions->get()->atoms.atom);
958 }
959
960 return title;
961 }
962
do_top(bool bVerbose,const char * topfile,const char * topppfile,t_gromppopts * opts,bool bZero,t_symtab * symtab,gmx::ArrayRef<InteractionsOfType> interactions,int * combination_rule,double * repulsion_power,real * fudgeQQ,PreprocessingAtomTypes * atypes,std::vector<MoleculeInformation> * molinfo,std::unique_ptr<MoleculeInformation> * intermolecular_interactions,const t_inputrec * ir,std::vector<gmx_molblock_t> * molblock,bool * ffParametrizedWithHBondConstraints,warninp * wi,const gmx::MDLogger & logger)963 char** do_top(bool bVerbose,
964 const char* topfile,
965 const char* topppfile,
966 t_gromppopts* opts,
967 bool bZero,
968 t_symtab* symtab,
969 gmx::ArrayRef<InteractionsOfType> interactions,
970 int* combination_rule,
971 double* repulsion_power,
972 real* fudgeQQ,
973 PreprocessingAtomTypes* atypes,
974 std::vector<MoleculeInformation>* molinfo,
975 std::unique_ptr<MoleculeInformation>* intermolecular_interactions,
976 const t_inputrec* ir,
977 std::vector<gmx_molblock_t>* molblock,
978 bool* ffParametrizedWithHBondConstraints,
979 warninp* wi,
980 const gmx::MDLogger& logger)
981 {
982 /* Tmpfile might contain a long path */
983 const char* tmpfile;
984 char** title;
985
986 if (topppfile)
987 {
988 tmpfile = topppfile;
989 }
990 else
991 {
992 tmpfile = nullptr;
993 }
994
995 if (bVerbose)
996 {
997 GMX_LOG(logger.info).asParagraph().appendTextFormatted("processing topology...");
998 }
999 title = read_topol(topfile, tmpfile, opts->define, opts->include, symtab, atypes, molinfo,
1000 intermolecular_interactions, interactions, combination_rule, repulsion_power,
1001 opts, fudgeQQ, molblock, ffParametrizedWithHBondConstraints,
1002 ir->efep != efepNO, bZero, EEL_FULL(ir->coulombtype), wi, logger);
1003
1004 if ((*combination_rule != eCOMB_GEOMETRIC) && (ir->vdwtype == evdwUSER))
1005 {
1006 warning(wi,
1007 "Using sigma/epsilon based combination rules with"
1008 " user supplied potential function may produce unwanted"
1009 " results");
1010 }
1011
1012 return title;
1013 }
1014
1015 /*! \brief
1016 * Exclude molecular interactions for QM atoms handled by MiMic.
1017 *
1018 * Update the exclusion lists to include all QM atoms of this molecule,
1019 * replace bonds between QM atoms with CONNBOND and
1020 * set charges of QM atoms to 0.
1021 *
1022 * \param[in,out] molt molecule type with QM atoms
1023 * \param[in] grpnr group informatio
1024 * \param[in,out] ir input record
1025 * \param[in] logger Handle to logging interface.
1026 */
generate_qmexcl_moltype(gmx_moltype_t * molt,const unsigned char * grpnr,t_inputrec * ir,const gmx::MDLogger & logger)1027 static void generate_qmexcl_moltype(gmx_moltype_t* molt,
1028 const unsigned char* grpnr,
1029 t_inputrec* ir,
1030 const gmx::MDLogger& logger)
1031 {
1032 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1033
1034 /* generates the exclusions between the individual QM atoms, as
1035 * these interactions should be handled by the QM subroutines and
1036 * not by the gromacs routines
1037 */
1038 int qm_max = 0, qm_nr = 0, link_nr = 0;
1039 int * qm_arr = nullptr, *link_arr = nullptr;
1040 bool* bQMMM;
1041
1042 /* First we search and select the QM atoms in an qm_arr array that
1043 * we use to create the exclusions.
1044 *
1045 * we take the possibility into account that a user has defined more
1046 * than one QM group:
1047 *
1048 * for that we also need to do this an ugly work-about just in case
1049 * the QM group contains the entire system...
1050 */
1051
1052 /* we first search for all the QM atoms and put them in an array
1053 */
1054 for (int j = 0; j < ir->opts.ngQM; j++)
1055 {
1056 for (int i = 0; i < molt->atoms.nr; i++)
1057 {
1058 if (qm_nr >= qm_max)
1059 {
1060 qm_max += 100;
1061 srenew(qm_arr, qm_max);
1062 }
1063 if ((grpnr ? grpnr[i] : 0) == j)
1064 {
1065 qm_arr[qm_nr++] = i;
1066 molt->atoms.atom[i].q = 0.0;
1067 molt->atoms.atom[i].qB = 0.0;
1068 }
1069 }
1070 }
1071 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1072 * QM/not QM. We first set all elements to false. Afterwards we use
1073 * the qm_arr to change the elements corresponding to the QM atoms
1074 * to TRUE.
1075 */
1076 snew(bQMMM, molt->atoms.nr);
1077 for (int i = 0; i < molt->atoms.nr; i++)
1078 {
1079 bQMMM[i] = FALSE;
1080 }
1081 for (int i = 0; i < qm_nr; i++)
1082 {
1083 bQMMM[qm_arr[i]] = TRUE;
1084 }
1085
1086 /* We remove all bonded interactions (i.e. bonds,
1087 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1088 * are removed is as follows: if the interaction invloves 2 atoms,
1089 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1090 * it is removed if at least two of the atoms are QM atoms, if the
1091 * interaction involves 4 atoms, it is removed if there are at least
1092 * 2 QM atoms. Since this routine is called once before any forces
1093 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1094 * be rewritten at this poitn without any problem. 25-9-2002 */
1095
1096 /* first check whether we already have CONNBONDS.
1097 * Note that if we don't, we don't add a param entry and set ftype=0,
1098 * which is ok, since CONNBONDS does not use parameters.
1099 */
1100 int ftype_connbond = 0;
1101 int ind_connbond = 0;
1102 if (!molt->ilist[F_CONNBONDS].empty())
1103 {
1104 GMX_LOG(logger.info)
1105 .asParagraph()
1106 .appendTextFormatted("nr. of CONNBONDS present already: %d",
1107 molt->ilist[F_CONNBONDS].size() / 3);
1108 ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
1109 ind_connbond = molt->ilist[F_CONNBONDS].size();
1110 }
1111 /* now we delete all bonded interactions, except the ones describing
1112 * a chemical bond. These are converted to CONNBONDS
1113 */
1114 for (int ftype = 0; ftype < F_NRE; ftype++)
1115 {
1116 if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_CONNBONDS)
1117 {
1118 continue;
1119 }
1120 int nratoms = interaction_function[ftype].nratoms;
1121 int j = 0;
1122 while (j < molt->ilist[ftype].size())
1123 {
1124 bool bexcl;
1125
1126 if (nratoms == 2)
1127 {
1128 /* Remove an interaction between two atoms when both are
1129 * in the QM region. Note that we don't have to worry about
1130 * link atoms here, as they won't have 2-atom interactions.
1131 */
1132 int a1 = molt->ilist[ftype].iatoms[1 + j + 0];
1133 int a2 = molt->ilist[ftype].iatoms[1 + j + 1];
1134 bexcl = (bQMMM[a1] && bQMMM[a2]);
1135 /* A chemical bond between two QM atoms will be copied to
1136 * the F_CONNBONDS list, for reasons mentioned above.
1137 */
1138 if (bexcl && IS_CHEMBOND(ftype))
1139 {
1140 InteractionList& ilist = molt->ilist[F_CONNBONDS];
1141 ilist.iatoms.resize(ind_connbond + 3);
1142 ilist.iatoms[ind_connbond++] = ftype_connbond;
1143 ilist.iatoms[ind_connbond++] = a1;
1144 ilist.iatoms[ind_connbond++] = a2;
1145 }
1146 }
1147 else
1148 {
1149 /* MM interactions have to be excluded if they are included
1150 * in the QM already. Because we use a link atom (H atom)
1151 * when the QM/MM boundary runs through a chemical bond, this
1152 * means that as long as one atom is MM, we still exclude,
1153 * as the interaction is included in the QM via:
1154 * QMatom1-QMatom2-QMatom-3-Linkatom.
1155 */
1156 int numQmAtoms = 0;
1157 for (int jj = j + 1; jj < j + 1 + nratoms; jj++)
1158 {
1159 if (bQMMM[molt->ilist[ftype].iatoms[jj]])
1160 {
1161 numQmAtoms++;
1162 }
1163 }
1164
1165 /* MiMiC treats link atoms as quantum atoms - therefore
1166 * we do not need do additional exclusions here */
1167 bexcl = numQmAtoms == nratoms;
1168
1169 if (bexcl && ftype == F_SETTLE)
1170 {
1171 gmx_fatal(FARGS,
1172 "Can not apply QM to molecules with SETTLE, replace the moleculetype "
1173 "using QM and SETTLE by one without SETTLE");
1174 }
1175 }
1176 if (bexcl)
1177 {
1178 /* since the interaction involves QM atoms, these should be
1179 * removed from the MM ilist
1180 */
1181 InteractionList& ilist = molt->ilist[ftype];
1182 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1183 {
1184 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1185 }
1186 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1187 }
1188 else
1189 {
1190 j += nratoms + 1; /* the +1 is for the functype */
1191 }
1192 }
1193 }
1194 /* Now, we search for atoms bonded to a QM atom because we also want
1195 * to exclude their nonbonded interactions with the QM atoms. The
1196 * reason for this is that this interaction is accounted for in the
1197 * linkatoms interaction with the QMatoms and would be counted
1198 * twice. */
1199
1200 /* creating the exclusion block for the QM atoms. Each QM atom has
1201 * as excluded elements all the other QMatoms (and itself).
1202 */
1203 t_blocka qmexcl;
1204 qmexcl.nr = molt->atoms.nr;
1205 qmexcl.nra = qm_nr * (qm_nr + link_nr) + link_nr * qm_nr;
1206 snew(qmexcl.index, qmexcl.nr + 1);
1207 snew(qmexcl.a, qmexcl.nra);
1208 int j = 0;
1209 for (int i = 0; i < qmexcl.nr; i++)
1210 {
1211 qmexcl.index[i] = j;
1212 if (bQMMM[i])
1213 {
1214 for (int k = 0; k < qm_nr; k++)
1215 {
1216 qmexcl.a[k + j] = qm_arr[k];
1217 }
1218 for (int k = 0; k < link_nr; k++)
1219 {
1220 qmexcl.a[qm_nr + k + j] = link_arr[k];
1221 }
1222 j += (qm_nr + link_nr);
1223 }
1224 }
1225 qmexcl.index[qmexcl.nr] = j;
1226
1227 /* and merging with the exclusions already present in sys.
1228 */
1229
1230 std::vector<gmx::ExclusionBlock> qmexcl2(molt->atoms.nr);
1231 gmx::blockaToExclusionBlocks(&qmexcl, qmexcl2);
1232 gmx::mergeExclusions(&(molt->excls), qmexcl2);
1233
1234 /* Finally, we also need to get rid of the pair interactions of the
1235 * classical atom bonded to the boundary QM atoms with the QMatoms,
1236 * as this interaction is already accounted for by the QM, so also
1237 * here we run the risk of double counting! We proceed in a similar
1238 * way as we did above for the other bonded interactions: */
1239 for (int i = F_LJ14; i < F_COUL14; i++)
1240 {
1241 int nratoms = interaction_function[i].nratoms;
1242 int j = 0;
1243 while (j < molt->ilist[i].size())
1244 {
1245 int a1 = molt->ilist[i].iatoms[j + 1];
1246 int a2 = molt->ilist[i].iatoms[j + 2];
1247 bool bexcl = (bQMMM[a1] && bQMMM[a2]);
1248 if (bexcl)
1249 {
1250 /* since the interaction involves QM atoms, these should be
1251 * removed from the MM ilist
1252 */
1253 InteractionList& ilist = molt->ilist[i];
1254 for (int k = j; k < ilist.size() - (nratoms + 1); k++)
1255 {
1256 ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
1257 }
1258 ilist.iatoms.resize(ilist.size() - (nratoms + 1));
1259 }
1260 else
1261 {
1262 j += nratoms + 1; /* the +1 is for the functype */
1263 }
1264 }
1265 }
1266
1267 free(qm_arr);
1268 free(bQMMM);
1269 free(link_arr);
1270 } /* generate_qmexcl */
1271
generate_qmexcl(gmx_mtop_t * sys,t_inputrec * ir,const gmx::MDLogger & logger)1272 void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, const gmx::MDLogger& logger)
1273 {
1274 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1275 */
1276
1277 unsigned char* grpnr;
1278 int mol, nat_mol, nr_mol_with_qm_atoms = 0;
1279 gmx_molblock_t* molb;
1280 bool bQMMM;
1281 int index_offset = 0;
1282 int qm_nr = 0;
1283
1284 grpnr = sys->groups.groupNumbers[SimulationAtomGroupType::QuantumMechanics].data();
1285
1286 for (size_t mb = 0; mb < sys->molblock.size(); mb++)
1287 {
1288 molb = &sys->molblock[mb];
1289 nat_mol = sys->moltype[molb->type].atoms.nr;
1290 for (mol = 0; mol < molb->nmol; mol++)
1291 {
1292 bQMMM = FALSE;
1293 for (int i = 0; i < nat_mol; i++)
1294 {
1295 if ((grpnr ? grpnr[i] : 0) < (ir->opts.ngQM))
1296 {
1297 bQMMM = TRUE;
1298 qm_nr++;
1299 }
1300 }
1301
1302 if (bQMMM)
1303 {
1304 nr_mol_with_qm_atoms++;
1305 if (molb->nmol > 1)
1306 {
1307 /* We need to split this molblock */
1308 if (mol > 0)
1309 {
1310 /* Split the molblock at this molecule */
1311 auto pos = sys->molblock.begin() + mb + 1;
1312 sys->molblock.insert(pos, sys->molblock[mb]);
1313 sys->molblock[mb].nmol = mol;
1314 sys->molblock[mb + 1].nmol -= mol;
1315 mb++;
1316 molb = &sys->molblock[mb];
1317 }
1318 if (molb->nmol > 1)
1319 {
1320 /* Split the molblock after this molecule */
1321 auto pos = sys->molblock.begin() + mb + 1;
1322 sys->molblock.insert(pos, sys->molblock[mb]);
1323 molb = &sys->molblock[mb];
1324 sys->molblock[mb].nmol = 1;
1325 sys->molblock[mb + 1].nmol -= 1;
1326 }
1327
1328 /* Create a copy of a moltype for a molecule
1329 * containing QM atoms and append it in the end of the list
1330 */
1331 std::vector<gmx_moltype_t> temp(sys->moltype.size());
1332 for (size_t i = 0; i < sys->moltype.size(); ++i)
1333 {
1334 copy_moltype(&sys->moltype[i], &temp[i]);
1335 }
1336 sys->moltype.resize(sys->moltype.size() + 1);
1337 for (size_t i = 0; i < temp.size(); ++i)
1338 {
1339 copy_moltype(&temp[i], &sys->moltype[i]);
1340 }
1341 copy_moltype(&sys->moltype[molb->type], &sys->moltype.back());
1342 /* Copy the exclusions to a new array, since this is the only
1343 * thing that needs to be modified for QMMM.
1344 */
1345 sys->moltype.back().excls = sys->moltype[molb->type].excls;
1346 /* Set the molecule type for the QMMM molblock */
1347 molb->type = sys->moltype.size() - 1;
1348 }
1349 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir, logger);
1350 }
1351 if (grpnr)
1352 {
1353 grpnr += nat_mol;
1354 }
1355 index_offset += nat_mol;
1356 }
1357 }
1358 }
1359