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-2008, The GROMACS development team.
6 * Copyright (c) 2013,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.
37 /* This file is completely threadsafe - keep it that way! */
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/gmxlib/readinp.h"
48 #include "gromacs/gmxpreprocess/fflibutil.h"
49 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
50 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
51 #include "gromacs/gmxpreprocess/notset.h"
52 #include "gromacs/gmxpreprocess/pdb2top.h"
53 #include "gromacs/gmxpreprocess/toppush.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/mdtypes/md_enums.h"
57 #include "gromacs/utility/cstringutil.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/futil.h"
60 #include "gromacs/utility/smalloc.h"
62 static void rd_nm2type_file(const char *fn
, int *nnm
, t_nm2type
**nmp
)
66 char libfilename
[128];
67 char format
[128], f1
[128];
68 char buf
[1024], elem
[16], type
[16], nbbuf
[16], **newbuf
;
69 int i
, nb
, nnnm
, line
= 1;
71 t_nm2type
*nm2t
= NULL
;
76 gmx_fatal(FARGS
, "Can not find %s in library directory", fn
);
83 /* Read a line from the file */
84 bCont
= (fgets2(buf
, 1023, fp
) != NULL
);
90 if (sscanf(buf
, "%s%s%lf%lf%d", elem
, type
, &qq
, &mm
, &nb
) == 5)
92 /* If we can read the first four, there probably is more */
94 snew(nm2t
[nnnm
].blen
, nb
);
98 strcpy(format
, "%*s%*s%*s%*s%*s");
99 for (i
= 0; (i
< nb
); i
++)
101 /* Complicated format statement */
104 if (sscanf(buf
, f1
, nbbuf
, &(nm2t
[nnnm
].blen
[i
])) != 2)
106 gmx_fatal(FARGS
, "Error on line %d of %s", line
, libfilename
);
108 newbuf
[i
] = gmx_strdup(nbbuf
);
109 strcat(format
, "%*s%*s");
116 nm2t
[nnnm
].elem
= gmx_strdup(elem
);
117 nm2t
[nnnm
].type
= gmx_strdup(type
);
120 nm2t
[nnnm
].nbonds
= nb
;
121 nm2t
[nnnm
].bond
= newbuf
;
134 t_nm2type
*rd_nm2type(const char *ffdir
, int *nnm
)
140 nff
= fflib_search_file_end(ffdir
, ".n2t", FALSE
, &ff
);
143 for (f
= 0; f
< nff
; f
++)
145 rd_nm2type_file(ff
[f
], nnm
, &nm
);
153 void dump_nm2type(FILE *fp
, int nnm
, t_nm2type nm2t
[])
157 fprintf(fp
, "; nm2type database\n");
158 for (i
= 0; (i
< nnm
); i
++)
160 fprintf(fp
, "%-8s %-8s %8.4f %8.4f %-4d",
161 nm2t
[i
].elem
, nm2t
[i
].type
,
162 nm2t
[i
].q
, nm2t
[i
].m
, nm2t
[i
].nbonds
);
163 for (j
= 0; (j
< nm2t
[i
].nbonds
); j
++)
165 fprintf(fp
, " %-5s %6.4f", nm2t
[i
].bond
[j
], nm2t
[i
].blen
[j
]);
172 ematchNone
, ematchWild
, ematchElem
, ematchExact
, ematchNR
175 static int match_str(const char *atom
, const char *template_string
)
177 if (!atom
|| !template_string
)
181 else if (gmx_strcasecmp(atom
, template_string
) == 0)
185 else if (atom
[0] == template_string
[0])
189 else if (strcmp(template_string
, "*") == 0)
199 int nm2type(int nnm
, t_nm2type nm2t
[], struct t_symtab
*tab
, t_atoms
*atoms
,
200 gpp_atomtype_t atype
, int *nbonds
, t_params
*bonds
)
204 int i
, j
, k
, m
, n
, nresolved
, nb
, maxbond
, ai
, aj
, best
, im
, nqual
[2][ematchNR
];
205 int *bbb
, *n_mask
, *m_mask
, **match
;
206 char *aname_i
, *aname_m
, *aname_n
, *type
;
212 for (i
= 0; (i
< atoms
->nr
); i
++)
214 maxbond
= std::max(maxbond
, nbonds
[i
]);
218 fprintf(debug
, "Max number of bonds per atom is %d\n", maxbond
);
221 snew(n_mask
, maxbond
);
222 snew(m_mask
, maxbond
);
223 snew(match
, maxbond
);
224 for (i
= 0; (i
< maxbond
); i
++)
226 snew(match
[i
], maxbond
);
230 for (i
= 0; (i
< atoms
->nr
); i
++)
232 aname_i
= *atoms
->atomname
[i
];
234 for (j
= 0; (j
< bonds
->nr
); j
++)
236 ai
= bonds
->param
[j
].ai();
237 aj
= bonds
->param
[j
].aj();
249 gmx_fatal(FARGS
, "Counting number of bonds nb = %d, nbonds[%d] = %d",
254 fprintf(debug
, "%4s has bonds to", aname_i
);
255 for (j
= 0; (j
< nb
); j
++)
257 fprintf(debug
, " %4s", *atoms
->atomname
[bbb
[j
]]);
259 fprintf(debug
, "\n");
262 for (k
= 0; (k
< ematchNR
); k
++)
267 /* First check for names */
268 for (k
= 0; (k
< nnm
); k
++)
270 if (nm2t
[k
].nbonds
== nb
)
272 im
= match_str(*atoms
->atomname
[i
], nm2t
[k
].elem
);
275 for (j
= 0; (j
< ematchNR
); j
++)
280 /* Fill a matrix with matching quality */
281 for (m
= 0; (m
< nb
); m
++)
283 aname_m
= *atoms
->atomname
[bbb
[m
]];
284 for (n
= 0; (n
< nb
); n
++)
286 aname_n
= nm2t
[k
].bond
[n
];
287 match
[m
][n
] = match_str(aname_m
, aname_n
);
290 /* Now pick the best matches */
291 for (m
= 0; (m
< nb
); m
++)
296 for (j
= ematchNR
-1; (j
> 0); j
--)
298 for (m
= 0; (m
< nb
); m
++)
300 for (n
= 0; (n
< nb
); n
++)
302 if ((n_mask
[n
] == 0) &&
313 if ((nqual
[cur
][ematchExact
]+
314 nqual
[cur
][ematchElem
]+
315 nqual
[cur
][ematchWild
]) == nb
)
317 if ((nqual
[cur
][ematchExact
] > nqual
[prev
][ematchExact
]) ||
319 ((nqual
[cur
][ematchExact
] == nqual
[prev
][ematchExact
]) &&
320 (nqual
[cur
][ematchElem
] > nqual
[prev
][ematchElem
])) ||
322 ((nqual
[cur
][ematchExact
] == nqual
[prev
][ematchExact
]) &&
323 (nqual
[cur
][ematchElem
] == nqual
[prev
][ematchElem
]) &&
324 (nqual
[cur
][ematchWild
] > nqual
[prev
][ematchWild
])))
340 type
= nm2t
[best
].type
;
342 if ((k
= get_atomtype_type(type
, atype
)) == NOTSET
)
344 atoms
->atom
[i
].qB
= alpha
;
345 atoms
->atom
[i
].m
= atoms
->atom
[i
].mB
= mm
;
346 k
= add_atomtype(atype
, tab
, &(atoms
->atom
[i
]), type
, param
,
347 atoms
->atom
[i
].type
, 0, 0, 0, atomnr
, 0, 0);
349 atoms
->atom
[i
].type
= k
;
350 atoms
->atom
[i
].typeB
= k
;
351 atoms
->atom
[i
].q
= qq
;
352 atoms
->atom
[i
].qB
= qq
;
353 atoms
->atom
[i
].m
= mm
;
354 atoms
->atom
[i
].mB
= mm
;
359 fprintf(stderr
, "Can not find forcefield for atom %s-%d with %d bonds\n",
360 *atoms
->atomname
[i
], i
+1, nb
);