Bug Summary

File:gromacs/gmxpreprocess/topio.c
Location:line 1363, column 13
Description:Value stored to 'bexcl' is never read

Annotated Source Code

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, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37#ifdef HAVE_CONFIG_H1
38#include <config.h>
39#endif
40
41#include <assert.h>
42#include <ctype.h>
43#include <errno(*__errno_location ()).h>
44#include <math.h>
45#include <stdio.h>
46#include <stdlib.h>
47#include <string.h>
48
49#include <sys/types.h>
50
51#include "gromacs/utility/futil.h"
52#include "typedefs.h"
53#include "gromacs/utility/smalloc.h"
54#include "macros.h"
55#include "gromacs/fileio/gmxfio.h"
56#include "txtdump.h"
57#include "physics.h"
58#include "macros.h"
59#include "names.h"
60#include "gromacs/utility/cstringutil.h"
61#include "symtab.h"
62#include "gromacs/utility/fatalerror.h"
63#include "warninp.h"
64#include "vsite_parm.h"
65
66#include "grompp-impl.h"
67#include "toputil.h"
68#include "toppush.h"
69#include "topdirs.h"
70#include "gpp_nextnb.h"
71#include "topio.h"
72#include "topshake.h"
73#include "gmxcpp.h"
74#include "gpp_bond_atomtype.h"
75#include "genborn.h"
76#include "gromacs/math/utilities.h"
77
78#define OPENDIR'[' '[' /* starting sign for directive */
79#define CLOSEDIR']' ']' /* ending sign for directive */
80
81static void free_nbparam(t_nbparam **param, int nr)
82{
83 int i;
84
85 for (i = 0; i < nr; i++)
86 {
87 sfree(param[i])save_free("param[i]", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 87, (param[i]));
88 }
89 sfree(param)save_free("param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 89, (param));
90}
91
92static int copy_nbparams(t_nbparam **param, int ftype, t_params *plist, int nr)
93{
94 int i, j, f;
95 int nrfp, ncopy;
96
97 nrfp = NRFP(ftype)((interaction_function[(ftype)].nrfpA)+(interaction_function[
(ftype)].nrfpB));
98
99 ncopy = 0;
100 for (i = 0; i < nr; i++)
101 {
102 for (j = 0; j <= i; j++)
103 {
104 if (param[i][j].bSet)
105 {
106 for (f = 0; f < nrfp; f++)
107 {
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];
110 }
111 ncopy++;
112 }
113 }
114 }
115
116 return ncopy;
117}
118
119static void gen_pairs(t_params *nbs, t_params *pairs, real fudge, int comb)
120{
121 int i, j, ntp, nrfp, nrfpA, nrfpB, nnn;
122 real scaling;
123 ntp = nbs->nr;
124 nnn = sqrt(ntp);
125 nrfp = NRFP(F_LJ)((interaction_function[(F_LJ)].nrfpA)+(interaction_function[(
F_LJ)].nrfpB));
126 nrfpA = interaction_function[F_LJ14].nrfpA;
127 nrfpB = interaction_function[F_LJ14].nrfpB;
128 pairs->nr = ntp;
129
130 if ((nrfp != nrfpA) || (nrfpA != nrfpB))
131 {
132 gmx_incons("Number of force parameters in gen_pairs wrong")_gmx_error("incons", "Number of force parameters in gen_pairs wrong"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 132);
133 }
134
135 fprintf(stderrstderr, "Generating 1-4 interactions: fudge = %g\n", fudge);
136 if (debug)
137 {
138 fprintf(debug, "Fudge factor for 1-4 interactions: %g\n", fudge);
139 fprintf(debug, "Holy Cow! there are %d types\n", ntp);
140 }
141 snew(pairs->param, pairs->nr)(pairs->param) = save_calloc("pairs->param", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 141, (pairs->nr), sizeof(*(pairs->param)));
142 for (i = 0; (i < ntp); i++)
143 {
144 /* Copy param.a */
145 pairs->param[i].a[0] = i / nnn;
146 pairs->param[i].a[1] = i % nnn;
147 /* Copy normal and FEP parameters and multiply by fudge factor */
148
149
150
151 for (j = 0; (j < nrfp); j++)
152 {
153 /* If we are using sigma/epsilon values, only the epsilon values
154 * should be scaled, but not sigma.
155 * The sigma values have even indices 0,2, etc.
156 */
157 if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2 == 0))
158 {
159 scaling = 1.0;
160 }
161 else
162 {
163 scaling = fudge;
164 }
165
166 pairs->param[i].c[j] = scaling*nbs->param[i].c[j];
167 pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j];
168 }
169 }
170}
171
172double check_mol(gmx_mtop_t *mtop, warninp_t wi)
173{
174 char buf[256];
175 int i, mb, nmol, ri, pt;
176 double q;
177 real m;
178 t_atoms *atoms;
179
180 /* Check mass and charge */
181 q = 0.0;
182
183 for (mb = 0; mb < mtop->nmoltype; mb++)
184 {
185 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
186 nmol = mtop->molblock[mb].nmol;
187 for (i = 0; (i < atoms->nr); i++)
188 {
189 q += nmol*atoms->atom[i].q;
190 m = atoms->atom[i].m;
191 pt = atoms->atom[i].ptype;
192 /* If the particle is an atom or a nucleus it must have a mass,
193 * else, if it is a shell, a vsite or a bondshell it can have mass zero
194 */
195 if ((m <= 0.0) && ((pt == eptAtom) || (pt == eptNucleus)))
196 {
197 ri = atoms->atom[i].resind;
198 sprintf(buf, "atom %s (Res %s-%d) has mass %g\n",
199 *(atoms->atomname[i]),
200 *(atoms->resinfo[ri].name),
201 atoms->resinfo[ri].nr,
202 m);
203 warning_error(wi, buf);
204 }
205 else
206 if ((m != 0) && (pt == eptVSite))
207 {
208 ri = atoms->atom[i].resind;
209 sprintf(buf, "virtual site %s (Res %s-%d) has non-zero mass %g\n"
210 " Check your topology.\n",
211 *(atoms->atomname[i]),
212 *(atoms->resinfo[ri].name),
213 atoms->resinfo[ri].nr,
214 m);
215 warning_error(wi, buf);
216 /* The following statements make LINCS break! */
217 /* atoms->atom[i].m=0; */
218 }
219 }
220 }
221 return q;
222}
223
224static void sum_q(t_atoms *atoms, int n, double *qt, double *qBt)
225{
226 double qmolA, qmolB;
227 int i;
228
229 /* sum charge */
230 qmolA = 0;
231 qmolB = 0;
232 for (i = 0; i < atoms->nr; i++)
233 {
234 qmolA += atoms->atom[i].q;
235 qmolB += atoms->atom[i].qB;
236 }
237 /* Unfortunately an absolute comparison,
238 * but this avoids unnecessary warnings and gmx-users mails.
239 */
240 if (fabs(qmolA) >= 1e-6 || fabs(qmolB) >= 1e-6)
241 {
242 *qt += n*qmolA;
243 *qBt += n*qmolB;
244 }
245}
246
247static void get_nbparm(char *nb_str, char *comb_str, int *nb, int *comb,
248 warninp_t wi)
249{
250 int i;
251 char warn_buf[STRLEN4096];
252
253 *nb = -1;
254 for (i = 1; (i < eNBF_NR); i++)
255 {
256 if (gmx_strcasecmp(nb_str, enbf_names[i]) == 0)
257 {
258 *nb = i;
259 }
260 }
261 if (*nb == -1)
262 {
263 *nb = strtol(nb_str, NULL((void*)0), 10);
264 }
265 if ((*nb < 1) || (*nb >= eNBF_NR))
266 {
267 sprintf(warn_buf, "Invalid nonbond function selector '%s' using %s",
268 nb_str, enbf_names[1]);
269 warning_error(wi, warn_buf);
270 *nb = 1;
271 }
272 *comb = -1;
273 for (i = 1; (i < eCOMB_NR); i++)
274 {
275 if (gmx_strcasecmp(comb_str, ecomb_names[i]) == 0)
276 {
277 *comb = i;
278 }
279 }
280 if (*comb == -1)
281 {
282 *comb = strtol(comb_str, NULL((void*)0), 10);
283 }
284 if ((*comb < 1) || (*comb >= eCOMB_NR))
285 {
286 sprintf(warn_buf, "Invalid combination rule selector '%s' using %s",
287 comb_str, ecomb_names[1]);
288 warning_error(wi, warn_buf);
289 *comb = 1;
290 }
291}
292
293static char ** cpp_opts(const char *define, const char *include,
294 warninp_t wi)
295{
296 int n, len;
297 int ncppopts = 0;
298 const char *cppadds[2];
299 char **cppopts = NULL((void*)0);
300 const char *option[2] = { "-D", "-I" };
301 const char *nopt[2] = { "define", "include" };
302 const char *ptr;
303 const char *rptr;
304 char *buf;
305 char warn_buf[STRLEN4096];
306
307 cppadds[0] = define;
308 cppadds[1] = include;
309 for (n = 0; (n < 2); n++)
310 {
311 if (cppadds[n])
312 {
313 ptr = cppadds[n];
314 while (*ptr != '\0')
315 {
316 while ((*ptr != '\0') && isspace(*ptr)((*__ctype_b_loc ())[(int) ((*ptr))] & (unsigned short int
) _ISspace))
317 {
318 ptr++;
319 }
320 rptr = ptr;
321 while ((*rptr != '\0') && !isspace(*rptr)((*__ctype_b_loc ())[(int) ((*rptr))] & (unsigned short int
) _ISspace))
322 {
323 rptr++;
324 }
325 len = (rptr - ptr);
326 if (len > 2)
327 {
328 snew(buf, (len+1))(buf) = save_calloc("buf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 328, ((len+1)), sizeof(*(buf)));
329 strncpy(buf, ptr, len)__builtin_strncpy (buf, ptr, len);
330 if (strstr(ptr, option[n]) != ptr)
331 {
332 set_warning_line(wi, "mdp file", -1);
333 sprintf(warn_buf, "Malformed %s option %s", nopt[n], buf);
334 warning(wi, warn_buf);
335 }
336 else
337 {
338 srenew(cppopts, ++ncppopts)(cppopts) = save_realloc("cppopts", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 338, (cppopts), (++ncppopts), sizeof(*(cppopts)));
339 cppopts[ncppopts-1] = strdup(buf)(__extension__ (__builtin_constant_p (buf) && ((size_t
)(const void *)((buf) + 1) - (size_t)(const void *)(buf) == 1
) ? (((const char *) (buf))[0] == '\0' ? (char *) calloc ((size_t
) 1, (size_t) 1) : ({ size_t __len = strlen (buf) + 1; char *
__retval = (char *) malloc (__len); if (__retval != ((void*)0
)) __retval = (char *) memcpy (__retval, buf, __len); __retval
; })) : __strdup (buf)));
340 }
341 sfree(buf)save_free("buf", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 341, (buf));
342 ptr = rptr;
343 }
344 }
345 }
346 }
347 srenew(cppopts, ++ncppopts)(cppopts) = save_realloc("cppopts", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 347, (cppopts), (++ncppopts), sizeof(*(cppopts)));
348 cppopts[ncppopts-1] = NULL((void*)0);
349
350 return cppopts;
351}
352
353
354int
355find_gb_bondlength(t_params *plist, int ai, int aj, real *length)
356{
357 int i, j, a1, a2;
358
359 int found = 0;
360 int status;
361
362 for (i = 0; i < F_NRE && !found; i++)
363 {
364 if (IS_CHEMBOND(i)(interaction_function[(i)].nratoms == 2 && (interaction_function
[(i)].flags & 1<<3)))
365 {
366 for (j = 0; j < plist[i].nr; j++)
367 {
368 a1 = plist[i].param[j].a[0];
369 a2 = plist[i].param[j].a[1];
370
371 if ( (a1 == ai && a2 == aj) || (a1 == aj && a2 == ai))
372 {
373 /* Equilibrium bond distance */
374 *length = plist[i].param[j].c[0];
375 found = 1;
376 }
377 }
378 }
379 }
380 status = !found;
381
382 return status;
383}
384
385
386int
387find_gb_anglelength(t_params *plist, int ai, int ak, real *length)
388{
389 int i, j, a1, a2, a3;
390 real r12, r23, a123;
391 int found = 0;
392 int status, status1, status2;
393
394 r12 = r23 = 0;
395
396 for (i = 0; i < F_NRE && !found; i++)
397 {
398 if (IS_ANGLE(i)(interaction_function[(i)].nratoms == 3 && (interaction_function
[(i)].flags & 1<<5)))
399 {
400 for (j = 0; j < plist[i].nr; j++)
401 {
402 a1 = plist[i].param[j].a[0];
403 a2 = plist[i].param[j].a[1];
404 a3 = plist[i].param[j].a[2];
405
406 /* We dont care what the middle atom is, but use it below */
407 if ( (a1 == ai && a3 == ak) || (a1 == ak && a3 == ai) )
408 {
409 /* Equilibrium bond distance */
410 a123 = plist[i].param[j].c[0];
411 /* Use middle atom to find reference distances r12 and r23 */
412 status1 = find_gb_bondlength(plist, a1, a2, &r12);
413 status2 = find_gb_bondlength(plist, a2, a3, &r23);
414
415 if (status1 == 0 && status2 == 0)
416 {
417 /* cosine theorem to get r13 */
418 *length = sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG(180.0/3.14159265358979323846))));
419 found = 1;
420 }
421 }
422 }
423 }
424 }
425 status = !found;
426
427 return status;
428}
429
430int
431generate_gb_exclusion_interactions(t_molinfo *mi, gpp_atomtype_t atype, t_nextnb *nnb)
432{
433 int i, j, k, n, ai, aj, ti, tj;
434 int n12, n13, n14;
435 int ftype;
436 t_param param;
437 t_params * plist;
438 t_atoms * at;
439 real radiusi, radiusj;
440 real gb_radiusi, gb_radiusj;
441 real param_c2, param_c4;
442 real distance;
443
444 plist = mi->plist;
445 at = &mi->atoms;
446
447 for (n = 1; n <= nnb->nrex; n++)
448 {
449 switch (n)
450 {
451 case 1:
452 ftype = F_GB12;
453 param_c2 = STILL_P20.921*0.1*(4.184);
454 param_c4 = 0.8875;
455 break;
456 case 2:
457 ftype = F_GB13;
458 param_c2 = STILL_P36.211*0.1*(4.184);
459 param_c4 = 0.3516;
460 break;
461 default:
462 /* Put all higher-order exclusions into 1,4 list so we dont miss them */
463 ftype = F_GB14;
464 param_c2 = STILL_P36.211*0.1*(4.184);
465 param_c4 = 0.3516;
466 break;
467 }
468
469 for (ai = 0; ai < nnb->nr; ai++)
470 {
471 ti = at->atom[ai].type;
472 radiusi = get_atomtype_radius(ti, atype);
473 gb_radiusi = get_atomtype_gb_radius(ti, atype);
474
475 for (j = 0; j < nnb->nrexcl[ai][n]; j++)
476 {
477 aj = nnb->a[ai][n][j];
478
479 /* Only add the interactions once */
480 if (aj > ai)
481 {
482 tj = at->atom[aj].type;
483 radiusj = get_atomtype_radius(tj, atype);
484 gb_radiusj = get_atomtype_gb_radius(tj, atype);
485
486 /* There is an exclusion of type "ftype" between atoms ai and aj */
487 param.a[0] = ai;
488 param.a[1] = aj;
489
490 /* Reference distance, not used for 1-4 interactions */
491 switch (ftype)
492 {
493 case F_GB12:
494 if (find_gb_bondlength(plist, ai, aj, &distance) != 0)
495 {
496 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 496, "Cannot find bond length for atoms %d-%d", ai, aj);
497 }
498 break;
499 case F_GB13:
500 if (find_gb_anglelength(plist, ai, aj, &distance) != 0)
501 {
502 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 502, "Cannot find length for atoms %d-%d involved in angle", ai, aj);
503 }
504 break;
505 default:
506 distance = -1;
507 break;
508 }
509 /* Assign GB parameters */
510 /* Sum of radii */
511 param.c[0] = radiusi+radiusj;
512 /* Reference distance distance */
513 param.c[1] = distance;
514 /* Still parameter */
515 param.c[2] = param_c2;
516 /* GB radius */
517 param.c[3] = gb_radiusi+gb_radiusj;
518 /* Parameter */
519 param.c[4] = param_c4;
520
521 /* Add it to the parameter list */
522 add_param_to_list(&plist[ftype], &param);
523 }
524 }
525 }
526 }
527 return 0;
528}
529
530
531
532
533static char **read_topol(const char *infile, const char *outfile,
534 const char *define, const char *include,
535 t_symtab *symtab,
536 gpp_atomtype_t atype,
537 int *nrmols,
538 t_molinfo **molinfo,
539 t_params plist[],
540 int *combination_rule,
541 double *reppow,
542 t_gromppopts *opts,
543 real *fudgeQQ,
544 int *nmolblock,
545 gmx_molblock_t **molblock,
546 gmx_bool bFEP,
547 gmx_bool bGenborn,
548 gmx_bool bZero,
549 warninp_t wi)
550{
551 FILE *out;
552 int i, sl, nb_funct, comb;
553 char *pline = NULL((void*)0), **title = NULL((void*)0);
554 char line[STRLEN4096], errbuf[256], comb_str[256], nb_str[256];
555 char genpairs[32];
556 char *dirstr, *dummy2;
557 int nrcopies, nmol, nmolb = 0, nscan, ncombs, ncopy;
558 double fLJ, fQQ, fPOW;
559 gmx_molblock_t *molb = NULL((void*)0);
560 t_topology *block = NULL((void*)0);
561 t_molinfo *mi0 = NULL((void*)0);
562 DirStack *DS;
563 directive d, newd;
564 t_nbparam **nbparam, **pair;
565 t_block2 *block2;
566 real fudgeLJ = -1; /* Multiplication factor to generate 1-4 from LJ */
567 gmx_bool bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
568 double qt = 0, qBt = 0; /* total charge */
569 t_bond_atomtype batype;
570 int lastcg = -1;
571 int dcatt = -1, nmol_couple;
572 /* File handling variables */
573 int status, done;
574 gmx_cpp_t handle;
575 char *tmp_line = NULL((void*)0);
576 char warn_buf[STRLEN4096];
577 const char *floating_point_arithmetic_tip =
578 "Total charge should normally be an integer. See\n"
579 "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
580 "for discussion on how close it should be to an integer.\n";
581 /* We need to open the output file before opening the input file,
582 * because cpp_open_file can change the current working directory.
583 */
584 if (outfile)
585 {
586 out = gmx_fio_fopen(outfile, "w");
587 }
588 else
589 {
590 out = NULL((void*)0);
591 }
592
593 /* open input file */
594 status = cpp_open_file(infile, &handle, cpp_opts(define, include, wi));
595 if (status != 0)
596 {
597 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 597, cpp_error(&handle, status));
598 }
599
600 /* some local variables */
601 DS_Init(&DS); /* directive stack */
602 nmol = 0; /* no molecules yet... */
603 d = d_invalid; /* first thing should be a directive */
604 nbparam = NULL((void*)0); /* The temporary non-bonded matrix */
605 pair = NULL((void*)0); /* The temporary pair interaction matrix */
606 block2 = NULL((void*)0); /* the extra exclusions */
607 nb_funct = F_LJ;
608 *reppow = 12.0; /* Default value for repulsion power */
609
610 comb = 0;
611
612 /* Init the number of CMAP torsion angles and grid spacing */
613 plist->grid_spacing = 0;
614 plist->nc = 0;
615
616 bWarn_copy_A_B = bFEP;
617
618 batype = init_bond_atomtype();
619 /* parse the actual file */
620 bReadDefaults = FALSE0;
621 bGenPairs = FALSE0;
622 bReadMolType = FALSE0;
623 nmol_couple = 0;
624
625 do
626 {
627 status = cpp_read_line(&handle, STRLEN4096, line);
628 done = (status == eCPP_EOF);
629 if (!done)
630 {
631 if (status != eCPP_OK)
632 {
633 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 633, cpp_error(&handle, status));
634 }
635 else if (out)
636 {
637 fprintf(out, "%s\n", line);
638 }
639
640 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
641
642 pline = strdup(line)(__extension__ (__builtin_constant_p (line) && ((size_t
)(const void *)((line) + 1) - (size_t)(const void *)(line) ==
1) ? (((const char *) (line))[0] == '\0' ? (char *) calloc (
(size_t) 1, (size_t) 1) : ({ size_t __len = strlen (line) + 1
; char *__retval = (char *) malloc (__len); if (__retval != (
(void*)0)) __retval = (char *) memcpy (__retval, line, __len)
; __retval; })) : __strdup (line)));
643
644 /* Strip trailing '\' from pline, if it exists */
645 sl = strlen(pline);
646 if ((sl > 0) && (pline[sl-1] == CONTINUE'\\'))
647 {
648 pline[sl-1] = ' ';
649 }
650
651 /* build one long line from several fragments - necessary for CMAP */
652 while (continuing(line))
653 {
654 status = cpp_read_line(&handle, STRLEN4096, line);
655 set_warning_line(wi, cpp_cur_file(&handle), cpp_cur_linenr(&handle));
656
657 /* Since we depend on the '\' being present to continue to read, we copy line
658 * to a tmp string, strip the '\' from that string, and cat it to pline
659 */
660 tmp_line = strdup(line)(__extension__ (__builtin_constant_p (line) && ((size_t
)(const void *)((line) + 1) - (size_t)(const void *)(line) ==
1) ? (((const char *) (line))[0] == '\0' ? (char *) calloc (
(size_t) 1, (size_t) 1) : ({ size_t __len = strlen (line) + 1
; char *__retval = (char *) malloc (__len); if (__retval != (
(void*)0)) __retval = (char *) memcpy (__retval, line, __len)
; __retval; })) : __strdup (line)));
661
662 sl = strlen(tmp_line);
663 if ((sl > 0) && (tmp_line[sl-1] == CONTINUE'\\'))
664 {
665 tmp_line[sl-1] = ' ';
666 }
667
668 done = (status == eCPP_EOF);
669 if (!done)
670 {
671 if (status != eCPP_OK)
672 {
673 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 673, cpp_error(&handle, status));
674 }
675 else if (out)
676 {
677 fprintf(out, "%s\n", line);
678 }
679 }
680
681 srenew(pline, strlen(pline)+strlen(tmp_line)+1)(pline) = save_realloc("pline", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 681, (pline), (strlen(pline)+strlen(tmp_line)+1), sizeof(*(
pline)));
682 strcat(pline, tmp_line);
683 sfree(tmp_line)save_free("tmp_line", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 683, (tmp_line));
684 }
685
686 /* skip trailing and leading spaces and comment text */
687 strip_comment (pline);
688 trim (pline);
689
690 /* if there is something left... */
691 if ((int)strlen(pline) > 0)
692 {
693 if (pline[0] == OPENDIR'[')
694 {
695 /* A directive on this line: copy the directive
696 * without the brackets into dirstr, then
697 * skip spaces and tabs on either side of directive
698 */
699 dirstr = strdup((pline+1))(__extension__ (__builtin_constant_p ((pline+1)) && (
(size_t)(const void *)(((pline+1)) + 1) - (size_t)(const void
*)((pline+1)) == 1) ? (((const char *) ((pline+1)))[0] == '\0'
? (char *) calloc ((size_t) 1, (size_t) 1) : ({ size_t __len
= strlen ((pline+1)) + 1; char *__retval = (char *) malloc (
__len); if (__retval != ((void*)0)) __retval = (char *) memcpy
(__retval, (pline+1), __len); __retval; })) : __strdup ((pline
+1))));
700 if ((dummy2 = strchr (dirstr, CLOSEDIR)(__extension__ (__builtin_constant_p (']') && !__builtin_constant_p
(dirstr) && (']') == '\0' ? (char *) __rawmemchr (dirstr
, ']') : __builtin_strchr (dirstr, ']')))) != NULL((void*)0))
701 {
702 (*dummy2) = 0;
703 }
704 trim (dirstr);
705
706 if ((newd = str2dir(dirstr)) == d_invalid)
707 {
708 sprintf(errbuf, "Invalid directive %s", dirstr);
709 warning_error(wi, errbuf);
710 }
711 else
712 {
713 /* Directive found */
714 if (debug)
715 {
716 fprintf(debug, "found directive '%s'\n", dir2str(newd));
717 }
718 if (DS_Check_Order (DS, newd))
719 {
720 DS_Push (&DS, newd);
721 d = newd;
722 }
723 else
724 {
725 /* we should print here which directives should have
726 been present, and which actually are */
727 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 727, "%s\nInvalid order for directive %s",
728 cpp_error(&handle, eCPP_SYNTAX), dir2str(newd));
729 /* d = d_invalid; */
730 }
731 }
732 sfree(dirstr)save_free("dirstr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 732, (dirstr));
733 }
734 else if (d != d_invalid)
735 {
736 /* Not a directive, just a plain string
737 * use a gigantic switch to decode,
738 * if there is a valid directive!
739 */
740 switch (d)
741 {
742 case d_defaults:
743 if (bReadDefaults)
744 {
745 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 745, "%s\nFound a second defaults directive.\n",
746 cpp_error(&handle, eCPP_SYNTAX));
747 }
748 bReadDefaults = TRUE1;
749 nscan = sscanf(pline, "%s%s%s%lf%lf%lf",
750 nb_str, comb_str, genpairs, &fLJ, &fQQ, &fPOW);
751 if (nscan < 2)
752 {
753 too_few(wi)_too_few(wi, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 753);
754 }
755 else
756 {
757 bGenPairs = FALSE0;
758 fudgeLJ = 1.0;
759 *fudgeQQ = 1.0;
760
761 get_nbparm(nb_str, comb_str, &nb_funct, &comb, wi);
762 *combination_rule = comb;
763 if (nscan >= 3)
764 {
765 bGenPairs = (gmx_strncasecmp(genpairs, "Y", 1) == 0);
766 if (nb_funct != eNBF_LJ && bGenPairs)
767 {
768 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 768, "Generating pair parameters is only supported with LJ non-bonded interactions");
769 }
770 }
771 if (nscan >= 4)
772 {
773 fudgeLJ = fLJ;
774 }
775 if (nscan >= 5)
776 {
777 *fudgeQQ = fQQ;
778 }
779 if (nscan >= 6)
780 {
781 *reppow = fPOW;
782 }
783 }
784 nb_funct = ifunc_index(d_nonbond_params, nb_funct);
785
786 break;
787 case d_atomtypes:
788 push_at(symtab, atype, batype, pline, nb_funct,
789 &nbparam, bGenPairs ? &pair : NULL((void*)0), wi);
790 break;
791
792 case d_bondtypes:
793 push_bt(d, plist, 2, NULL((void*)0), batype, pline, wi);
794 break;
795 case d_constrainttypes:
796 push_bt(d, plist, 2, NULL((void*)0), batype, pline, wi);
797 break;
798 case d_pairtypes:
799 if (bGenPairs)
800 {
801 push_nbt(d, pair, atype, pline, F_LJ14, wi);
802 }
803 else
804 {
805 push_bt(d, plist, 2, atype, NULL((void*)0), pline, wi);
806 }
807 break;
808 case d_angletypes:
809 push_bt(d, plist, 3, NULL((void*)0), batype, pline, wi);
810 break;
811 case d_dihedraltypes:
812 /* Special routine that can read both 2 and 4 atom dihedral definitions. */
813 push_dihedraltype(d, plist, batype, pline, wi);
814 break;
815
816 case d_nonbond_params:
817 push_nbt(d, nbparam, atype, pline, nb_funct, wi);
818 break;
819 /*
820 case d_blocktype:
821 nblock++;
822 srenew(block,nblock);
823 srenew(blockinfo,nblock);
824 blk0=&(block[nblock-1]);
825 bi0=&(blockinfo[nblock-1]);
826 init_top(blk0);
827 init_molinfo(bi0);
828 push_molt(symtab,bi0,pline);
829 break;
830 */
831
832 case d_implicit_genborn_params:
833 push_gb_params(atype, pline, wi);
834 break;
835
836 case d_implicit_surface_params:
837 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 837, "Implicit surface directive not supported yet.");
838 break;
839
840 case d_cmaptypes:
841 push_cmaptype(d, plist, 5, atype, batype, pline, wi);
842 break;
843
844 case d_moleculetype:
845 {
846 if (!bReadMolType)
847 {
848 int ntype;
849 if (opts->couple_moltype != NULL((void*)0) &&
850 (opts->couple_lam0 == ecouplamNONE ||
851 opts->couple_lam0 == ecouplamQ ||
852 opts->couple_lam1 == ecouplamNONE ||
853 opts->couple_lam1 == ecouplamQ))
854 {
855 dcatt = add_atomtype_decoupled(symtab, atype,
856 &nbparam, bGenPairs ? &pair : NULL((void*)0));
857 }
858 ntype = get_atomtype_ntypes(atype);
859 ncombs = (ntype*(ntype+1))/2;
860 generate_nbparams(comb, nb_funct, &(plist[nb_funct]), atype, wi);
861 ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]),
862 ntype);
863 fprintf(stderrstderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs);
864 free_nbparam(nbparam, ntype);
865 if (bGenPairs)
866 {
867 gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, comb);
868 ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]),
869 ntype);
870 fprintf(stderrstderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs);
871 free_nbparam(pair, ntype);
872 }
873 /* Copy GBSA parameters to atomtype array? */
874
875 bReadMolType = TRUE1;
876 }
877
878 push_molt(symtab, &nmol, molinfo, pline, wi);
879 srenew(block2, nmol)(block2) = save_realloc("block2", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 879, (block2), (nmol), sizeof(*(block2)));
880 block2[nmol-1].nr = 0;
881 mi0 = &((*molinfo)[nmol-1]);
882 break;
883 }
884 case d_atoms:
885 push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi);
886 break;
887
888 case d_pairs:
889 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE0,
890 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
891 break;
892
893 case d_vsites2:
894 case d_vsites3:
895 case d_vsites4:
896 case d_bonds:
897 case d_angles:
898 case d_constraints:
899 case d_settles:
900 case d_position_restraints:
901 case d_angle_restraints:
902 case d_angle_restraints_z:
903 case d_distance_restraints:
904 case d_orientation_restraints:
905 case d_dihedral_restraints:
906 case d_dihedrals:
907 case d_polarization:
908 case d_water_polarization:
909 case d_thole_polarization:
910 push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE1,
911 bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi);
912 break;
913 case d_cmap:
914 push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi);
915 break;
916
917 case d_vsitesn:
918 push_vsitesn(d, mi0->plist, &(mi0->atoms), pline, wi);
919 break;
920 case d_exclusions:
921 assert(block2)((void) (0));
922 if (!block2[nmol-1].nr)
923 {
924 init_block2(&(block2[nmol-1]), mi0->atoms.nr);
925 }
926 push_excl(pline, &(block2[nmol-1]));
927 break;
928 case d_system:
929 trim(pline);
930 title = put_symtab(symtab, pline);
931 break;
932 case d_molecules:
933 {
934 int whichmol;
935 gmx_bool bCouple;
936
937 push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
938 mi0 = &((*molinfo)[whichmol]);
939 srenew(molb, nmolb+1)(molb) = save_realloc("molb", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 939, (molb), (nmolb+1), sizeof(*(molb)));
940 molb[nmolb].type = whichmol;
941 molb[nmolb].nmol = nrcopies;
942 nmolb++;
943
944 bCouple = (opts->couple_moltype != NULL((void*)0) &&
945 (gmx_strcasecmp("system", opts->couple_moltype) == 0 ||
946 gmx_strcasecmp(*(mi0->name), opts->couple_moltype) == 0));
947 if (bCouple)
948 {
949 nmol_couple += nrcopies;
950 }
951
952 if (mi0->atoms.nr == 0)
953 {
954 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 954, "Molecule type '%s' contains no atoms",
955 *mi0->name);
956 }
957 fprintf(stderrstderr,
958 "Excluding %d bonded neighbours molecule type '%s'\n",
959 mi0->nrexcl, *mi0->name);
960 sum_q(&mi0->atoms, nrcopies, &qt, &qBt);
961 if (!mi0->bProcessed)
962 {
963 t_nextnb nnb;
964 generate_excl(mi0->nrexcl,
965 mi0->atoms.nr,
966 mi0->plist,
967 &nnb,
968 &(mi0->excls));
969 merge_excl(&(mi0->excls), &(block2[whichmol]));
970 done_block2(&(block2[whichmol]));
971 make_shake(mi0->plist, &mi0->atoms, opts->nshake);
972
973
974
975 /* nnb contains information about first,2nd,3rd bonded neighbors.
976 * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
977 */
978 if (bGenborn == TRUE1)
979 {
980 generate_gb_exclusion_interactions(mi0, atype, &nnb);
981 }
982
983 done_nnb(&nnb);
984
985 if (bCouple)
986 {
987 convert_moltype_couple(mi0, dcatt, *fudgeQQ,
988 opts->couple_lam0, opts->couple_lam1,
989 opts->bCoupleIntra,
990 nb_funct, &(plist[nb_funct]));
991 }
992 stupid_fill_block(&mi0->mols, mi0->atoms.nr, TRUE1);
993 mi0->bProcessed = TRUE1;
994 }
995 break;
996 }
997 default:
998 fprintf (stderrstderr, "case: %d\n", d);
999 gmx_incons("unknown directive")_gmx_error("incons", "unknown directive", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 999);
1000 }
1001 }
1002 }
1003 sfree(pline)save_free("pline", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1003, (pline));
1004 pline = NULL((void*)0);
1005 }
1006 }
1007 while (!done);
1008 status = cpp_close_file(&handle);
1009 if (status != eCPP_OK)
1010 {
1011 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1011, cpp_error(&handle, status));
1012 }
1013 cpp_done();
1014 if (out)
1015 {
1016 gmx_fio_fclose(out);
1017 }
1018
1019 if (opts->couple_moltype)
1020 {
1021 if (nmol_couple == 0)
1022 {
1023 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1023, "Did not find any molecules of type '%s' for coupling",
1024 opts->couple_moltype);
1025 }
1026 fprintf(stderrstderr, "Coupling %d copies of molecule type '%s'\n",
1027 nmol_couple, opts->couple_moltype);
1028 }
1029
1030 /* this is not very clean, but fixes core dump on empty system name */
1031 if (!title)
1032 {
1033 title = put_symtab(symtab, "");
1034 }
1035 if (fabs(qt) > 1e-4)
1036 {
1037 sprintf(warn_buf, "System has non-zero total charge: %.6f\n%s\n", qt, floating_point_arithmetic_tip);
1038 warning_note(wi, warn_buf);
1039 }
1040 if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt, qt, 1e-6))
1041 {
1042 sprintf(warn_buf, "State B has non-zero total charge: %.6f\n%s\n", qBt, floating_point_arithmetic_tip);
1043 warning_note(wi, warn_buf);
1044 }
1045 DS_Done (&DS);
1046 for (i = 0; i < nmol; i++)
1047 {
1048 done_block2(&(block2[i]));
1049 }
1050 free(block2);
1051
1052 done_bond_atomtype(&batype);
1053
1054 *nrmols = nmol;
1055
1056 *nmolblock = nmolb;
1057 *molblock = molb;
1058
1059 return title;
1060}
1061
1062char **do_top(gmx_bool bVerbose,
1063 const char *topfile,
1064 const char *topppfile,
1065 t_gromppopts *opts,
1066 gmx_bool bZero,
1067 t_symtab *symtab,
1068 t_params plist[],
1069 int *combination_rule,
1070 double *repulsion_power,
1071 real *fudgeQQ,
1072 gpp_atomtype_t atype,
1073 int *nrmols,
1074 t_molinfo **molinfo,
1075 t_inputrec *ir,
1076 int *nmolblock,
1077 gmx_molblock_t **molblock,
1078 gmx_bool bGenborn,
1079 warninp_t wi)
1080{
1081 /* Tmpfile might contain a long path */
1082 const char *tmpfile;
1083 char **title;
1084
1085 if (topppfile)
1086 {
1087 tmpfile = topppfile;
1088 }
1089 else
1090 {
1091 tmpfile = NULL((void*)0);
1092 }
1093
1094 if (bVerbose)
1095 {
1096 printf("processing topology...\n");
1097 }
1098 title = read_topol(topfile, tmpfile, opts->define, opts->include,
1099 symtab, atype, nrmols, molinfo,
1100 plist, combination_rule, repulsion_power,
1101 opts, fudgeQQ, nmolblock, molblock,
1102 ir->efep != efepNO, bGenborn, bZero, wi);
1103 if ((*combination_rule != eCOMB_GEOMETRIC) &&
1104 (ir->vdwtype == evdwUSER))
1105 {
1106 warning(wi, "Using sigma/epsilon based combination rules with"
1107 " user supplied potential function may produce unwanted"
1108 " results");
1109 }
1110
1111 return title;
1112}
1113
1114
1115static void generate_qmexcl_moltype(gmx_moltype_t *molt, unsigned char *grpnr,
1116 t_inputrec *ir)
1117{
1118 /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1119
1120 /* generates the exclusions between the individual QM atoms, as
1121 * these interactions should be handled by the QM subroutines and
1122 * not by the gromacs routines
1123 */
1124 int
1125 i, j, l, k = 0, jmax, qm_max = 0, qm_nr = 0, nratoms = 0, link_nr = 0, link_max = 0;
1126 atom_id
1127 *qm_arr = NULL((void*)0), *link_arr = NULL((void*)0), a1, a2, a3, a4, ftype = 0;
1128 t_blocka
1129 qmexcl;
1130 t_block2
1131 qmexcl2;
1132 gmx_bool
1133 *bQMMM, *blink, bexcl;
1134
1135 /* First we search and select the QM atoms in an qm_arr array that
1136 * we use to create the exclusions.
1137 *
1138 * we take the possibility into account that a user has defined more
1139 * than one QM group:
1140 *
1141 * for that we also need to do this an ugly work-about just in case
1142 * the QM group contains the entire system...
1143 */
1144 jmax = ir->opts.ngQM;
1145
1146 /* we first search for all the QM atoms and put them in an array
1147 */
1148 for (j = 0; j < jmax; j++)
1149 {
1150 for (i = 0; i < molt->atoms.nr; i++)
1151 {
1152 if (qm_nr >= qm_max)
1153 {
1154 qm_max += 100;
1155 srenew(qm_arr, qm_max)(qm_arr) = save_realloc("qm_arr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1155, (qm_arr), (qm_max), sizeof(*(qm_arr)));
1156 }
1157 if ((grpnr ? grpnr[i] : 0) == j)
1158 {
1159 qm_arr[qm_nr++] = i;
1160 }
1161 }
1162 }
1163 /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1164 * QM/not QM. We first set all elements to false. Afterwards we use
1165 * the qm_arr to change the elements corresponding to the QM atoms
1166 * to TRUE.
1167 */
1168 snew(bQMMM, molt->atoms.nr)(bQMMM) = save_calloc("bQMMM", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1168, (molt->atoms.nr), sizeof(*(bQMMM)));
1169 for (i = 0; i < molt->atoms.nr; i++)
1170 {
1171 bQMMM[i] = FALSE0;
1172 }
1173 for (i = 0; i < qm_nr; i++)
1174 {
1175 bQMMM[qm_arr[i]] = TRUE1;
1176 }
1177
1178 /* We remove all bonded interactions (i.e. bonds,
1179 * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1180 * are removed is as follows: if the interaction invloves 2 atoms,
1181 * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1182 * it is removed if at least two of the atoms are QM atoms, if the
1183 * interaction involves 4 atoms, it is removed if there are at least
1184 * 2 QM atoms. Since this routine is called once before any forces
1185 * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1186 * be rewritten at this poitn without any problem. 25-9-2002 */
1187
1188 /* first check weter we already have CONNBONDS: */
1189 if (molt->ilist[F_CONNBONDS].nr != 0)
1190 {
1191 fprintf(stderrstderr, "nr. of CONNBONDS present already: %d\n",
1192 molt->ilist[F_CONNBONDS].nr/3);
1193 ftype = molt->ilist[F_CONNBONDS].iatoms[0];
1194 k = molt->ilist[F_CONNBONDS].nr;
1195 }
1196 /* now we delete all bonded interactions, except the ones describing
1197 * a chemical bond. These are converted to CONNBONDS
1198 */
1199 for (i = 0; i < F_LJ; i++)
1200 {
1201 if (i == F_CONNBONDS)
1202 {
1203 continue;
1204 }
1205 nratoms = interaction_function[i].nratoms;
1206 j = 0;
1207 while (j < molt->ilist[i].nr)
1208 {
1209 bexcl = FALSE0;
1210 switch (nratoms)
1211 {
1212 case 2:
1213 a1 = molt->ilist[i].iatoms[j+1];
1214 a2 = molt->ilist[i].iatoms[j+2];
1215 bexcl = (bQMMM[a1] && bQMMM[a2]);
1216 /* a bonded beteen two QM atoms will be copied to the
1217 * CONNBONDS list, for reasons mentioned above
1218 */
1219 if (bexcl && i < F_ANGLES)
1220 {
1221 srenew(molt->ilist[F_CONNBONDS].iatoms, k+3)(molt->ilist[F_CONNBONDS].iatoms) = save_realloc("molt->ilist[F_CONNBONDS].iatoms"
, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1221, (molt->ilist[F_CONNBONDS].iatoms), (k+3), sizeof(*
(molt->ilist[F_CONNBONDS].iatoms)));
1222 molt->ilist[F_CONNBONDS].nr += 3;
1223 molt->ilist[F_CONNBONDS].iatoms[k++] = ftype;
1224 molt->ilist[F_CONNBONDS].iatoms[k++] = a1;
1225 molt->ilist[F_CONNBONDS].iatoms[k++] = a2;
1226 }
1227 break;
1228 case 3:
1229 a1 = molt->ilist[i].iatoms[j+1];
1230 a2 = molt->ilist[i].iatoms[j+2];
1231 a3 = molt->ilist[i].iatoms[j+3];
1232 bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1233 (bQMMM[a1] && bQMMM[a3]) ||
1234 (bQMMM[a2] && bQMMM[a3]));
1235 break;
1236 case 4:
1237 a1 = molt->ilist[i].iatoms[j+1];
1238 a2 = molt->ilist[i].iatoms[j+2];
1239 a3 = molt->ilist[i].iatoms[j+3];
1240 a4 = molt->ilist[i].iatoms[j+4];
1241 bexcl = ((bQMMM[a1] && bQMMM[a2] && bQMMM[a3]) ||
1242 (bQMMM[a1] && bQMMM[a2] && bQMMM[a4]) ||
1243 (bQMMM[a1] && bQMMM[a3] && bQMMM[a4]) ||
1244 (bQMMM[a2] && bQMMM[a3] && bQMMM[a4]));
1245 break;
1246 default:
1247 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1247, "no such bonded interactions with %d atoms\n", nratoms);
1248 }
1249 if (bexcl)
1250 {
1251 /* since the interaction involves QM atoms, these should be
1252 * removed from the MM ilist
1253 */
1254 molt->ilist[i].nr -= (nratoms+1);
1255 for (l = j; l < molt->ilist[i].nr; l++)
1256 {
1257 molt->ilist[i].iatoms[l] = molt->ilist[i].iatoms[l+(nratoms+1)];
1258 }
1259 }
1260 else
1261 {
1262 j += nratoms+1; /* the +1 is for the functype */
1263 }
1264 }
1265 }
1266 /* Now, we search for atoms bonded to a QM atom because we also want
1267 * to exclude their nonbonded interactions with the QM atoms. The
1268 * reason for this is that this interaction is accounted for in the
1269 * linkatoms interaction with the QMatoms and would be counted
1270 * twice. */
1271
1272 for (i = 0; i < F_NRE; i++)
1273 {
1274 if (IS_CHEMBOND(i)(interaction_function[(i)].nratoms == 2 && (interaction_function
[(i)].flags & 1<<3)))
1275 {
1276 j = 0;
1277 while (j < molt->ilist[i].nr)
1278 {
1279 a1 = molt->ilist[i].iatoms[j+1];
1280 a2 = molt->ilist[i].iatoms[j+2];
1281 if ((bQMMM[a1] && !bQMMM[a2]) || (!bQMMM[a1] && bQMMM[a2]))
1282 {
1283 if (link_nr >= link_max)
1284 {
1285 link_max += 10;
1286 srenew(link_arr, link_max)(link_arr) = save_realloc("link_arr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1286, (link_arr), (link_max), sizeof(*(link_arr)));
1287 }
1288 if (bQMMM[a1])
1289 {
1290 link_arr[link_nr++] = a2;
1291 }
1292 else
1293 {
1294 link_arr[link_nr++] = a1;
1295 }
1296 }
1297 j += 3;
1298 }
1299 }
1300 }
1301 snew(blink, molt->atoms.nr)(blink) = save_calloc("blink", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1301, (molt->atoms.nr), sizeof(*(blink)));
1302 for (i = 0; i < molt->atoms.nr; i++)
1303 {
1304 blink[i] = FALSE0;
1305 }
1306 for (i = 0; i < link_nr; i++)
1307 {
1308 blink[link_arr[i]] = TRUE1;
1309 }
1310 /* creating the exclusion block for the QM atoms. Each QM atom has
1311 * as excluded elements all the other QMatoms (and itself).
1312 */
1313 qmexcl.nr = molt->atoms.nr;
1314 qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1315 snew(qmexcl.index, qmexcl.nr+1)(qmexcl.index) = save_calloc("qmexcl.index", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1315, (qmexcl.nr+1), sizeof(*(qmexcl.index)));
1316 snew(qmexcl.a, qmexcl.nra)(qmexcl.a) = save_calloc("qmexcl.a", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1316, (qmexcl.nra), sizeof(*(qmexcl.a)));
1317 j = 0;
1318 for (i = 0; i < qmexcl.nr; i++)
1319 {
1320 qmexcl.index[i] = j;
1321 if (bQMMM[i])
1322 {
1323 for (k = 0; k < qm_nr; k++)
1324 {
1325 qmexcl.a[k+j] = qm_arr[k];
1326 }
1327 for (k = 0; k < link_nr; k++)
1328 {
1329 qmexcl.a[qm_nr+k+j] = link_arr[k];
1330 }
1331 j += (qm_nr+link_nr);
1332 }
1333 if (blink[i])
1334 {
1335 for (k = 0; k < qm_nr; k++)
1336 {
1337 qmexcl.a[k+j] = qm_arr[k];
1338 }
1339 j += qm_nr;
1340 }
1341 }
1342 qmexcl.index[qmexcl.nr] = j;
1343
1344 /* and merging with the exclusions already present in sys.
1345 */
1346
1347 init_block2(&qmexcl2, molt->atoms.nr);
1348 b_to_b2(&qmexcl, &qmexcl2);
1349 merge_excl(&(molt->excls), &qmexcl2);
1350 done_block2(&qmexcl2);
1351
1352 /* Finally, we also need to get rid of the pair interactions of the
1353 * classical atom bonded to the boundary QM atoms with the QMatoms,
1354 * as this interaction is already accounted for by the QM, so also
1355 * here we run the risk of double counting! We proceed in a similar
1356 * way as we did above for the other bonded interactions: */
1357 for (i = F_LJ14; i < F_COUL14; i++)
1358 {
1359 nratoms = interaction_function[i].nratoms;
1360 j = 0;
1361 while (j < molt->ilist[i].nr)
1362 {
1363 bexcl = FALSE0;
Value stored to 'bexcl' is never read
1364 a1 = molt->ilist[i].iatoms[j+1];
1365 a2 = molt->ilist[i].iatoms[j+2];
1366 bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1367 (blink[a1] && bQMMM[a2]) ||
1368 (bQMMM[a1] && blink[a2]));
1369 if (bexcl)
1370 {
1371 /* since the interaction involves QM atoms, these should be
1372 * removed from the MM ilist
1373 */
1374 molt->ilist[i].nr -= (nratoms+1);
1375 for (k = j; k < molt->ilist[i].nr; k++)
1376 {
1377 molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1378 }
1379 }
1380 else
1381 {
1382 j += nratoms+1; /* the +1 is for the functype */
1383 }
1384 }
1385 }
1386
1387 free(qm_arr);
1388 free(bQMMM);
1389 free(link_arr);
1390 free(blink);
1391} /* generate_qmexcl */
1392
1393void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi)
1394{
1395 /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1396 */
1397
1398 unsigned char *grpnr;
1399 int mb, mol, nat_mol, i, nr_mol_with_qm_atoms = 0;
1400 gmx_molblock_t *molb;
1401 gmx_bool bQMMM;
1402
1403 grpnr = sys->groups.grpnr[egcQMMM];
1404
1405 for (mb = 0; mb < sys->nmolblock; mb++)
1406 {
1407 molb = &sys->molblock[mb];
1408 nat_mol = sys->moltype[molb->type].atoms.nr;
1409 for (mol = 0; mol < molb->nmol; mol++)
1410 {
1411 bQMMM = FALSE0;
1412 for (i = 0; i < nat_mol; i++)
1413 {
1414 if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM)
1415 {
1416 bQMMM = TRUE1;
1417 }
1418 }
1419 if (bQMMM)
1420 {
1421 nr_mol_with_qm_atoms++;
1422 if (molb->nmol > 1)
1423 {
1424 /* We need to split this molblock */
1425 if (mol > 0)
1426 {
1427 /* Split the molblock at this molecule */
1428 sys->nmolblock++;
1429 srenew(sys->molblock, sys->nmolblock)(sys->molblock) = save_realloc("sys->molblock", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1429, (sys->molblock), (sys->nmolblock), sizeof(*(sys
->molblock)));
1430 for (i = sys->nmolblock-2; i >= mb; i--)
1431 {
1432 sys->molblock[i+1] = sys->molblock[i];
1433 }
1434 sys->molblock[mb ].nmol = mol;
1435 sys->molblock[mb+1].nmol -= mol;
1436 mb++;
1437 molb = &sys->molblock[mb];
1438 }
1439 if (molb->nmol > 1)
1440 {
1441 /* Split the molblock after this molecule */
1442 sys->nmolblock++;
1443 srenew(sys->molblock, sys->nmolblock)(sys->molblock) = save_realloc("sys->molblock", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1443, (sys->molblock), (sys->nmolblock), sizeof(*(sys
->molblock)));
1444 molb = &sys->molblock[mb];
1445 for (i = sys->nmolblock-2; i >= mb; i--)
1446 {
1447 sys->molblock[i+1] = sys->molblock[i];
1448 }
1449 sys->molblock[mb ].nmol = 1;
1450 sys->molblock[mb+1].nmol -= 1;
1451 }
1452
1453 /* Add a moltype for the QMMM molecule */
1454 sys->nmoltype++;
1455 srenew(sys->moltype, sys->nmoltype)(sys->moltype) = save_realloc("sys->moltype", "/home/alexxy/Develop/gromacs/src/gromacs/gmxpreprocess/topio.c"
, 1455, (sys->moltype), (sys->nmoltype), sizeof(*(sys->
moltype)));
1456 /* Copy the moltype struct */
1457 sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
1458 /* Copy the exclusions to a new array, since this is the only
1459 * thing that needs to be modified for QMMM.
1460 */
1461 copy_blocka(&sys->moltype[molb->type ].excls,
1462 &sys->moltype[sys->nmoltype-1].excls);
1463 /* Set the molecule type for the QMMM molblock */
1464 molb->type = sys->nmoltype - 1;
1465 }
1466 generate_qmexcl_moltype(&sys->moltype[molb->type], grpnr, ir);
1467 }
1468 if (grpnr)
1469 {
1470 grpnr += nat_mol;
1471 }
1472 }
1473 }
1474 if (nr_mol_with_qm_atoms > 1)
1475 {
1476 /* generate a warning is there are QM atoms in different
1477 * topologies. In this case it is not possible at this stage to
1478 * mutualy exclude the non-bonded interactions via the
1479 * exclusions (AFAIK). Instead, the user is advised to use the
1480 * energy group exclusions in the mdp file
1481 */
1482 warning_note(wi,
1483 "\nThe QM subsystem is divided over multiple topologies. "
1484 "The mutual non-bonded interactions cannot be excluded. "
1485 "There are two ways to achieve this:\n\n"
1486 "1) merge the topologies, such that the atoms of the QM "
1487 "subsystem are all present in one single topology file. "
1488 "In this case this warning will dissappear\n\n"
1489 "2) exclude the non-bonded interactions explicitly via the "
1490 "energygrp-excl option in the mdp file. if this is the case "
1491 "this warning may be ignored"
1492 "\n\n");
1493 }
1494}