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) 2011,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gpp_atomtype.h"
46 #include "gromacs/gmxpreprocess/notset.h"
47 #include "gromacs/gmxpreprocess/topdirs.h"
48 #include "gromacs/gmxpreprocess/toputil.h"
49 #include "gromacs/math/vecdump.h"
50 #include "gromacs/topology/ifunc.h"
51 #include "gromacs/topology/symtab.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/smalloc.h"
56 typedef struct gpp_atomtype
{
57 int nr
; /* The number of atomtypes */
58 t_atom
*atom
; /* Array of atoms */
59 char ***atomname
; /* Names of the atomtypes */
60 t_param
*nb
; /* Nonbonded force default params */
61 int *bondatomtype
; /* The bond_atomtype for each atomtype */
62 real
*radius
; /* Radius for GBSA stuff */
63 real
*vol
; /* Effective volume for GBSA */
64 real
*surftens
; /* Surface tension with water, for GBSA */
65 real
*gb_radius
; /* Radius for Still model */
66 real
*S_hct
; /* Overlap factor for HCT model */
67 int *atomnumber
; /* Atomic number, used for QM/MM */
70 int get_atomtype_type(const char *str
, gpp_atomtype_t ga
)
74 /* Atom types are always case sensitive */
75 for (i
= 0; (i
< ga
->nr
); i
++)
77 if (strcmp(str
, *(ga
->atomname
[i
])) == 0)
86 int get_atomtype_ntypes(gpp_atomtype_t ga
)
91 char *get_atomtype_name(int nt
, gpp_atomtype_t ga
)
93 if ((nt
< 0) || (nt
>= ga
->nr
))
98 return *(ga
->atomname
[nt
]);
101 real
get_atomtype_massA(int nt
, gpp_atomtype_t ga
)
103 if ((nt
< 0) || (nt
>= ga
->nr
))
108 return ga
->atom
[nt
].m
;
111 real
get_atomtype_massB(int nt
, gpp_atomtype_t ga
)
113 if ((nt
< 0) || (nt
>= ga
->nr
))
118 return ga
->atom
[nt
].mB
;
121 real
get_atomtype_qA(int nt
, gpp_atomtype_t ga
)
123 if ((nt
< 0) || (nt
>= ga
->nr
))
128 return ga
->atom
[nt
].q
;
131 real
get_atomtype_qB(int nt
, gpp_atomtype_t ga
)
133 if ((nt
< 0) || (nt
>= ga
->nr
))
138 return ga
->atom
[nt
].qB
;
141 int get_atomtype_ptype(int nt
, gpp_atomtype_t ga
)
143 if ((nt
< 0) || (nt
>= ga
->nr
))
148 return ga
->atom
[nt
].ptype
;
151 int get_atomtype_batype(int nt
, gpp_atomtype_t ga
)
153 if ((nt
< 0) || (nt
>= ga
->nr
))
158 return ga
->bondatomtype
[nt
];
161 int get_atomtype_atomnumber(int nt
, gpp_atomtype_t ga
)
163 if ((nt
< 0) || (nt
>= ga
->nr
))
168 return ga
->atomnumber
[nt
];
171 real
get_atomtype_radius(int nt
, gpp_atomtype_t ga
)
173 if ((nt
< 0) || (nt
>= ga
->nr
))
178 return ga
->radius
[nt
];
181 real
get_atomtype_vol(int nt
, gpp_atomtype_t ga
)
183 if ((nt
< 0) || (nt
>= ga
->nr
))
191 real
get_atomtype_surftens(int nt
, gpp_atomtype_t ga
)
193 if ((nt
< 0) || (nt
>= ga
->nr
))
198 return ga
->surftens
[nt
];
201 real
get_atomtype_gb_radius(int nt
, gpp_atomtype_t ga
)
203 if ((nt
< 0) || (nt
>= ga
->nr
))
208 return ga
->gb_radius
[nt
];
211 real
get_atomtype_S_hct(int nt
, gpp_atomtype_t ga
)
213 if ((nt
< 0) || (nt
>= ga
->nr
))
218 return ga
->S_hct
[nt
];
221 real
get_atomtype_nbparam(int nt
, int param
, gpp_atomtype_t ga
)
223 if ((nt
< 0) || (nt
>= ga
->nr
))
227 if ((param
< 0) || (param
>= MAXFORCEPARAM
))
231 return ga
->nb
[nt
].c
[param
];
234 gpp_atomtype_t
init_atomtype(void)
244 ga
->bondatomtype
= NULL
;
248 ga
->atomnumber
= NULL
;
249 ga
->gb_radius
= NULL
;
256 set_atomtype_gbparam(gpp_atomtype_t ga
, int i
,
257 real radius
, real vol
, real surftens
,
258 real gb_radius
, real S_hct
)
260 if ( (i
< 0) || (i
>= ga
->nr
))
265 ga
->radius
[i
] = radius
;
267 ga
->surftens
[i
] = surftens
;
268 ga
->gb_radius
[i
] = gb_radius
;
269 ga
->S_hct
[i
] = S_hct
;
275 int set_atomtype(int nt
, gpp_atomtype_t ga
, t_symtab
*tab
,
276 t_atom
*a
, const char *name
, t_param
*nb
,
278 real radius
, real vol
, real surftens
, int atomnumber
,
279 real gb_radius
, real S_hct
)
281 if ((nt
< 0) || (nt
>= ga
->nr
))
287 ga
->atomname
[nt
] = put_symtab(tab
, name
);
289 ga
->bondatomtype
[nt
] = bondatomtype
;
290 ga
->radius
[nt
] = radius
;
292 ga
->surftens
[nt
] = surftens
;
293 ga
->atomnumber
[nt
] = atomnumber
;
294 ga
->gb_radius
[nt
] = gb_radius
;
295 ga
->S_hct
[nt
] = S_hct
;
300 int add_atomtype(gpp_atomtype_t ga
, t_symtab
*tab
,
301 t_atom
*a
, const char *name
, t_param
*nb
,
303 real radius
, real vol
, real surftens
, int atomnumber
,
304 real gb_radius
, real S_hct
)
308 for (i
= 0; (i
< ga
->nr
); i
++)
310 if (strcmp(*ga
->atomname
[i
], name
) == 0)
314 fprintf(debug
, "Trying to add atomtype %s again. Skipping it.\n", name
);
322 srenew(ga
->atom
, ga
->nr
);
323 srenew(ga
->atomname
, ga
->nr
);
324 srenew(ga
->nb
, ga
->nr
);
325 srenew(ga
->bondatomtype
, ga
->nr
);
326 srenew(ga
->radius
, ga
->nr
);
327 srenew(ga
->vol
, ga
->nr
);
328 srenew(ga
->surftens
, ga
->nr
);
329 srenew(ga
->atomnumber
, ga
->nr
);
330 srenew(ga
->gb_radius
, ga
->nr
);
331 srenew(ga
->S_hct
, ga
->nr
);
333 return set_atomtype(ga
->nr
-1, ga
, tab
, a
, name
, nb
, bondatomtype
, radius
,
334 vol
, surftens
, atomnumber
, gb_radius
, S_hct
);
342 void print_at (FILE * out
, gpp_atomtype_t ga
)
345 t_atom
*atom
= ga
->atom
;
346 t_param
*nb
= ga
->nb
;
348 fprintf (out
, "[ %s ]\n", dir2str(d_atomtypes
));
349 fprintf (out
, "; %6s %8s %8s %8s %12s %12s\n",
350 "type", "mass", "charge", "particle", "c6", "c12");
351 for (i
= 0; (i
< ga
->nr
); i
++)
353 fprintf(out
, "%8s %8.3f %8.3f %8s %12e %12e\n",
354 *(ga
->atomname
[i
]), atom
[i
].m
, atom
[i
].q
, "A",
355 nb
[i
].c0(), nb
[i
].c1());
361 void done_atomtype(gpp_atomtype_t ga
)
366 sfree(ga
->bondatomtype
);
369 sfree(ga
->gb_radius
);
372 sfree(ga
->atomnumber
);
377 static int search_atomtypes(gpp_atomtype_t ga
, int *n
, int typelist
[],
379 t_param param
[], int ftype
)
381 int i
, nn
, nrfp
, j
, k
, ntype
, tli
;
382 gmx_bool bFound
= FALSE
;
386 ntype
= get_atomtype_ntypes(ga
);
388 for (i
= 0; (i
< nn
); i
++)
390 if (typelist
[i
] == thistype
)
392 /* This type number has already been added */
396 /* Otherwise, check if the parameters are identical to any previously added type */
399 for (j
= 0; j
< ntype
&& bFound
; j
++)
401 /* Check nonbonded parameters */
402 for (k
= 0; k
< nrfp
&& bFound
; k
++)
404 bFound
= (param
[ntype
*typelist
[i
]+j
].c
[k
] == param
[ntype
*thistype
+j
].c
[k
]);
407 /* Check radius, volume, surftens */
410 (get_atomtype_radius(tli
, ga
) == get_atomtype_radius(thistype
, ga
)) &&
411 (get_atomtype_vol(tli
, ga
) == get_atomtype_vol(thistype
, ga
)) &&
412 (get_atomtype_surftens(tli
, ga
) == get_atomtype_surftens(thistype
, ga
)) &&
413 (get_atomtype_atomnumber(tli
, ga
) == get_atomtype_atomnumber(thistype
, ga
)) &&
414 (get_atomtype_gb_radius(tli
, ga
) == get_atomtype_gb_radius(thistype
, ga
)) &&
415 (get_atomtype_S_hct(tli
, ga
) == get_atomtype_S_hct(thistype
, ga
));
427 fprintf(debug
, "Renumbering atomtype %d to %d\n", thistype
, nn
);
431 gmx_fatal(FARGS
, "Atomtype horror n = %d, %s, %d", nn
, __FILE__
, __LINE__
);
433 typelist
[nn
] = thistype
;
441 void renum_atype(t_params plist
[], gmx_mtop_t
*mtop
,
443 gpp_atomtype_t ga
, gmx_bool bVerbose
)
445 int i
, j
, k
, l
, molt
, mi
, mj
, nat
, nrfp
, ftype
, ntype
;
455 char ***new_atomname
;
457 ntype
= get_atomtype_ntypes(ga
);
458 snew(typelist
, ntype
);
462 fprintf(stderr
, "renumbering atomtypes...\n");
465 /* Since the bonded interactions have been assigned now,
466 * we want to reduce the number of atom types by merging
467 * ones with identical nonbonded interactions, in addition
468 * to removing unused ones.
470 * With Generalized-Born electrostatics, or implicit solvent
471 * we also check that the atomtype radius, effective_volume
472 * and surface tension match.
474 * With QM/MM we also check that the atom numbers match
477 /* Get nonbonded interaction type */
478 if (plist
[F_LJ
].nr
> 0)
487 /* Renumber atomtypes by first making a list of which ones are actually used.
488 * We provide the list of nonbonded parameters so search_atomtypes
489 * can determine if two types should be merged.
492 for (molt
= 0; molt
< mtop
->nmoltype
; molt
++)
494 atoms
= &mtop
->moltype
[molt
].atoms
;
495 for (i
= 0; (i
< atoms
->nr
); i
++)
497 atoms
->atom
[i
].type
=
498 search_atomtypes(ga
, &nat
, typelist
, atoms
->atom
[i
].type
,
499 plist
[ftype
].param
, ftype
);
500 atoms
->atom
[i
].typeB
=
501 search_atomtypes(ga
, &nat
, typelist
, atoms
->atom
[i
].typeB
,
502 plist
[ftype
].param
, ftype
);
506 for (i
= 0; i
< 2; i
++)
508 if (wall_atomtype
[i
] >= 0)
510 wall_atomtype
[i
] = search_atomtypes(ga
, &nat
, typelist
, wall_atomtype
[i
],
511 plist
[ftype
].param
, ftype
);
515 snew(new_radius
, nat
);
517 snew(new_surftens
, nat
);
518 snew(new_atomnumber
, nat
);
519 snew(new_gb_radius
, nat
);
520 snew(new_S_hct
, nat
);
521 snew(new_atomname
, nat
);
523 /* We now have a list of unique atomtypes in typelist */
527 pr_ivec(debug
, 0, "typelist", typelist
, nat
, TRUE
);
532 snew(nbsnew
, plist
[ftype
].nr
);
536 for (i
= k
= 0; (i
< nat
); i
++)
539 for (j
= 0; (j
< nat
); j
++, k
++)
542 for (l
= 0; (l
< nrfp
); l
++)
544 nbsnew
[k
].c
[l
] = plist
[ftype
].param
[ntype
*mi
+mj
].c
[l
];
547 new_radius
[i
] = get_atomtype_radius(mi
, ga
);
548 new_vol
[i
] = get_atomtype_vol(mi
, ga
);
549 new_surftens
[i
] = get_atomtype_surftens(mi
, ga
);
550 new_atomnumber
[i
] = get_atomtype_atomnumber(mi
, ga
);
551 new_gb_radius
[i
] = get_atomtype_gb_radius(mi
, ga
);
552 new_S_hct
[i
] = get_atomtype_S_hct(mi
, ga
);
553 new_atomname
[i
] = ga
->atomname
[mi
];
556 for (i
= 0; (i
< nat
*nat
); i
++)
558 for (l
= 0; (l
< nrfp
); l
++)
560 plist
[ftype
].param
[i
].c
[l
] = nbsnew
[i
].c
[l
];
564 mtop
->ffparams
.atnr
= nat
;
569 sfree(ga
->atomnumber
);
570 sfree(ga
->gb_radius
);
572 /* Dangling atomname pointers ? */
575 ga
->radius
= new_radius
;
577 ga
->surftens
= new_surftens
;
578 ga
->atomnumber
= new_atomnumber
;
579 ga
->gb_radius
= new_gb_radius
;
580 ga
->S_hct
= new_S_hct
;
581 ga
->atomname
= new_atomname
;
589 void copy_atomtype_atomtypes(gpp_atomtype_t ga
, t_atomtypes
*atomtypes
)
593 /* Copy the atomtype data to the topology atomtype list */
594 ntype
= get_atomtype_ntypes(ga
);
595 atomtypes
->nr
= ntype
;
596 snew(atomtypes
->radius
, ntype
);
597 snew(atomtypes
->vol
, ntype
);
598 snew(atomtypes
->surftens
, ntype
);
599 snew(atomtypes
->atomnumber
, ntype
);
600 snew(atomtypes
->gb_radius
, ntype
);
601 snew(atomtypes
->S_hct
, ntype
);
603 for (i
= 0; i
< ntype
; i
++)
605 atomtypes
->radius
[i
] = ga
->radius
[i
];
606 atomtypes
->vol
[i
] = ga
->vol
[i
];
607 atomtypes
->surftens
[i
] = ga
->surftens
[i
];
608 atomtypes
->atomnumber
[i
] = ga
->atomnumber
[i
];
609 atomtypes
->gb_radius
[i
] = ga
->gb_radius
[i
];
610 atomtypes
->S_hct
[i
] = ga
->S_hct
[i
];