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) 2012,2014,2015,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
49 #include "gromacs/gmxpreprocess/notset.h"
50 #include "gromacs/gmxpreprocess/topdirs.h"
51 #include "gromacs/topology/block.h"
52 #include "gromacs/topology/ifunc.h"
53 #include "gromacs/topology/symtab.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/gmxassert.h"
56 #include "gromacs/utility/smalloc.h"
60 void set_p_string(t_param *p, const char *s)
64 if (strlen(s) < sizeof(p->s)-1)
66 strncpy(p->s, s, sizeof(p->s));
70 gmx_fatal(FARGS, "Increase MAXSLEN in the grompp code to at least %d,"
71 " or shorten your definition of bonds like %s to at most %d",
72 strlen(s)+1, s, MAXSLEN-1);
81 void pr_alloc (int extra, t_params *pr)
85 /* get new space for arrays */
88 gmx_fatal(FARGS, "Trying to make array smaller.\n");
94 GMX_ASSERT(pr->nr != 0 || pr->param == NULL, "Invalid t_params object");
95 if (pr->nr+extra > pr->maxnr)
97 pr->maxnr = std::max(static_cast<int>(1.2*pr->maxnr), pr->maxnr + extra);
98 srenew(pr->param, pr->maxnr);
99 for (i = pr->nr; (i < pr->maxnr); i++)
101 for (j = 0; (j < MAXATOMLIST); j++)
103 pr->param[i].a[j] = 0;
105 for (j = 0; (j < MAXFORCEPARAM); j++)
107 pr->param[i].c[j] = 0;
109 set_p_string(&(pr->param[i]), "");
114 void init_plist(t_params plist[])
118 for (i = 0; (i < F_NRE); i++)
122 plist[i].param = nullptr;
126 plist[i].cmap = nullptr;
127 plist[i].grid_spacing = 0;
130 plist[i].cmap_types = nullptr;
134 void done_plist(t_params *plist)
136 for (int i = 0; i < F_NRE; i++)
138 t_params *pl = &plist[i];
141 sfree(pl->cmap_types);
145 void cp_param(t_param *dest, t_param *src)
149 for (j = 0; (j < MAXATOMLIST); j++)
151 dest->a[j] = src->a[j];
153 for (j = 0; (j < MAXFORCEPARAM); j++)
155 dest->c[j] = src->c[j];
157 strncpy(dest->s, src->s, sizeof(dest->s));
160 void add_param_to_list(t_params *list, t_param *b)
164 /* allocate one position extra */
167 /* fill the arrays */
168 for (j = 0; (j < MAXFORCEPARAM); j++)
170 list->param[list->nr].c[j] = b->c[j];
172 for (j = 0; (j < MAXATOMLIST); j++)
174 list->param[list->nr].a[j] = b->a[j];
176 memset(list->param[list->nr].s, 0, sizeof(list->param[list->nr].s));
182 void init_molinfo(t_molinfo *mol)
185 mol->excl_set = FALSE;
186 mol->bProcessed = FALSE;
187 init_plist(mol->plist);
188 init_block(&mol->cgs);
189 init_block(&mol->mols);
190 init_blocka(&mol->excls);
191 init_atom(&mol->atoms);
196 void done_mi(t_molinfo *mi)
198 done_atom (&(mi->atoms));
199 done_block(&(mi->cgs));
200 done_block(&(mi->mols));
201 done_plist(mi->plist);
204 /* PRINTING STRUCTURES */
206 static void print_bt(FILE *out, directive d, gpp_atomtype_t at,
207 int ftype, int fsubtype, t_params plist[],
210 /* This dihp is a DIRTY patch because the dih-types do not use
211 * all four atoms to determine the type.
213 const int dihp[2][2] = { { 1, 2 }, { 0, 3 } };
215 int i, j, f, nral, nrfp;
216 gmx_bool bDih = FALSE, bSwapParity;
218 bt = &(plist[ftype]);
246 case F_CROSS_BOND_ANGLES:
249 case F_CROSS_BOND_BONDS:
299 fprintf(out, "[ %s ]\n", dir2str(d));
303 fprintf (out, "%3s %4s", "ai", "aj");
304 for (j = 2; (j < nral); j++)
306 fprintf (out, " %3c%c", 'a', 'i'+j);
311 for (j = 0; (j < 2); j++)
313 fprintf (out, "%3c%c", 'a', 'i'+dihp[f][j]);
317 fprintf (out, " funct");
318 for (j = 0; (j < nrfp); j++)
320 fprintf (out, " %12c%1d", 'c', j);
324 /* print bondtypes */
325 for (i = 0; (i < bt->nr); i++)
327 bSwapParity = (bt->param[i].c0() == NOTSET) && (bt->param[i].c1() == -1);
330 for (j = 0; (j < nral); j++)
332 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[j], at));
337 for (j = 0; (j < 2); j++)
339 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[dihp[f][j]], at));
342 fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1);
344 if (bt->param[i].s[0])
346 fprintf(out, " %s", bt->param[i].s);
350 for (j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++)
352 fprintf (out, "%13.6e ", bt->param[i].c[j]);
362 void print_blocka(FILE *out, const char *szName,
363 const char *szIndex, const char *szA,
368 fprintf (out, "; %s\n", szName);
369 fprintf (out, "; %4s %s\n", szIndex, szA);
370 for (i = 0; (i < block->nr); i++)
372 for (i = 0; (i < block->nr); i++)
374 fprintf (out, "%6d", i+1);
375 for (j = block->index[i]; (j < ((int)block->index[i+1])); j++)
377 fprintf (out, "%5d", block->a[j]+1);
385 void print_excl(FILE *out, int natoms, t_excls excls[])
392 for (i = 0; i < natoms && !have_excl; i++)
394 have_excl = (excls[i].nr > 0);
399 fprintf (out, "[ %s ]\n", dir2str(d_exclusions));
400 fprintf (out, "; %4s %s\n", "i", "excluded from i");
401 for (i = 0; i < natoms; i++)
405 fprintf (out, "%6d ", i+1);
406 for (j = 0; j < excls[i].nr; j++)
408 fprintf (out, " %5d", excls[i].e[j]+1);
418 static double get_residue_charge(const t_atoms *atoms, int at)
423 ri = atoms->atom[at].resind;
425 while (at < atoms->nr && atoms->atom[at].resind == ri)
427 q += atoms->atom[at].q;
434 void print_atoms(FILE *out, gpp_atomtype_t atype, t_atoms *at, int *cgnr,
435 gmx_bool bRTPresname)
443 as = dir2str(d_atoms);
444 fprintf(out, "[ %s ]\n", as);
445 fprintf(out, "; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
446 "nr", "type", "resnr", "residue", "atom", "cgnr", "charge", "mass", "typeB", "chargeB", "massB");
452 fprintf(debug, "This molecule has %d atoms and %d residues\n",
458 /* if the information is present... */
459 for (i = 0; (i < at->nr); i++)
461 ri = at->atom[i].resind;
462 if ((i == 0 || ri != at->atom[i-1].resind) &&
463 at->resinfo[ri].rtp != nullptr)
465 qres = get_residue_charge(at, i);
466 fprintf(out, "; residue %3d %-3s rtp %-4s q ",
468 *at->resinfo[ri].name,
469 *at->resinfo[ri].rtp);
470 if (fabs(qres) < 0.001)
472 fprintf(out, " %s", "0.0");
476 fprintf(out, "%+3.1f", qres);
480 tpA = at->atom[i].type;
481 if ((tpnmA = get_atomtype_name(tpA, atype)) == nullptr)
483 gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i);
486 /* This is true by construction, but static analysers don't know */
487 GMX_ASSERT(!bRTPresname || at->resinfo[at->atom[i].resind].rtp, "-rtpres did not have residue name available");
488 fprintf(out, "%6d %10s %6d%c %5s %6s %6d %10g %10g",
493 *(at->resinfo[at->atom[i].resind].rtp) :
494 *(at->resinfo[at->atom[i].resind].name),
495 *(at->atomname[i]), cgnr[i],
496 at->atom[i].q, at->atom[i].m);
497 if (PERTURBED(at->atom[i]))
499 tpB = at->atom[i].typeB;
500 if ((tpnmB = get_atomtype_name(tpB, atype)) == nullptr)
502 gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i);
504 fprintf(out, " %6s %10g %10g",
505 tpnmB, at->atom[i].qB, at->atom[i].mB);
507 qtot += (double)at->atom[i].q;
508 if (fabs(qtot) < 4*GMX_REAL_EPS)
512 fprintf(out, " ; qtot %.4g\n", qtot);
519 void print_bondeds(FILE *out, int natoms, directive d,
520 int ftype, int fsubtype, t_params plist[])
523 gpp_atomtype_t atype;
528 atype = init_atomtype();
532 for (i = 0; (i < natoms); i++)
535 sprintf(buf, "%4d", (i+1));
536 add_atomtype(atype, &stab, a, buf, param, 0, 0);
538 print_bt(out, d, atype, ftype, fsubtype, plist, TRUE);
543 done_atomtype(atype);