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_t_atoms(&mol
->atoms
, 0, FALSE
);
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
);