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, 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! */
44 #include "gromacs/math/utilities.h"
47 #include "gromacs/utility/cstringutil.h"
48 #include "gromacs/utility/smalloc.h"
49 #include "gromacs/fileio/confio.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/math/3dview.h"
57 #include "gpp_nextnb.h"
58 #include "gpp_atomtype.h"
59 #include "fflibutil.h"
63 #include "gromacs/utility/fatalerror.h"
65 static void rd_nm2type_file(const char *fn
, int *nnm
, t_nm2type
**nmp
)
69 char libfilename
[128];
70 char format
[128], f1
[128];
71 char buf
[1024], elem
[16], type
[16], nbbuf
[16], **newbuf
;
72 int i
, nb
, nnnm
, line
= 1;
74 t_nm2type
*nm2t
= NULL
;
79 gmx_fatal(FARGS
, "Can not find %s in library directory", fn
);
86 /* Read a line from the file */
87 bCont
= (fgets2(buf
, 1023, fp
) != NULL
);
93 if (sscanf(buf
, "%s%s%lf%lf%d", elem
, type
, &qq
, &mm
, &nb
) == 5)
95 /* If we can read the first four, there probably is more */
97 snew(nm2t
[nnnm
].blen
, nb
);
101 strcpy(format
, "%*s%*s%*s%*s%*s");
102 for (i
= 0; (i
< nb
); i
++)
104 /* Complicated format statement */
107 if (sscanf(buf
, f1
, nbbuf
, &(nm2t
[nnnm
].blen
[i
])) != 2)
109 gmx_fatal(FARGS
, "Error on line %d of %s", line
, libfilename
);
111 newbuf
[i
] = strdup(nbbuf
);
112 strcat(format
, "%*s%*s");
119 nm2t
[nnnm
].elem
= strdup(elem
);
120 nm2t
[nnnm
].type
= strdup(type
);
123 nm2t
[nnnm
].nbonds
= nb
;
124 nm2t
[nnnm
].bond
= newbuf
;
137 t_nm2type
*rd_nm2type(const char *ffdir
, int *nnm
)
143 nff
= fflib_search_file_end(ffdir
, ".n2t", FALSE
, &ff
);
146 for (f
= 0; f
< nff
; f
++)
148 rd_nm2type_file(ff
[f
], nnm
, &nm
);
156 void dump_nm2type(FILE *fp
, int nnm
, t_nm2type nm2t
[])
160 fprintf(fp
, "; nm2type database\n");
161 for (i
= 0; (i
< nnm
); i
++)
163 fprintf(fp
, "%-8s %-8s %8.4f %8.4f %-4d",
164 nm2t
[i
].elem
, nm2t
[i
].type
,
165 nm2t
[i
].q
, nm2t
[i
].m
, nm2t
[i
].nbonds
);
166 for (j
= 0; (j
< nm2t
[i
].nbonds
); j
++)
168 fprintf(fp
, " %-5s %6.4f", nm2t
[i
].bond
[j
], nm2t
[i
].blen
[j
]);
175 ematchNone
, ematchWild
, ematchElem
, ematchExact
, ematchNR
178 static int match_str(const char *atom
, const char *template_string
)
180 if (!atom
|| !template_string
)
184 else if (gmx_strcasecmp(atom
, template_string
) == 0)
188 else if (atom
[0] == template_string
[0])
192 else if (strcmp(template_string
, "*") == 0)
202 int nm2type(int nnm
, t_nm2type nm2t
[], t_symtab
*tab
, t_atoms
*atoms
,
203 gpp_atomtype_t atype
, int *nbonds
, t_params
*bonds
)
207 int i
, j
, k
, m
, n
, nresolved
, nb
, maxbond
, ai
, aj
, best
, im
, nqual
[2][ematchNR
];
208 int *bbb
, *n_mask
, *m_mask
, **match
, **quality
;
209 char *aname_i
, *aname_m
, *aname_n
, *type
;
215 for (i
= 0; (i
< atoms
->nr
); i
++)
217 maxbond
= max(maxbond
, nbonds
[i
]);
221 fprintf(debug
, "Max number of bonds per atom is %d\n", maxbond
);
224 snew(n_mask
, maxbond
);
225 snew(m_mask
, maxbond
);
226 snew(match
, maxbond
);
227 for (i
= 0; (i
< maxbond
); i
++)
229 snew(match
[i
], maxbond
);
233 for (i
= 0; (i
< atoms
->nr
); i
++)
235 aname_i
= *atoms
->atomname
[i
];
237 for (j
= 0; (j
< bonds
->nr
); j
++)
239 ai
= bonds
->param
[j
].AI
;
240 aj
= bonds
->param
[j
].AJ
;
252 gmx_fatal(FARGS
, "Counting number of bonds nb = %d, nbonds[%d] = %d",
257 fprintf(debug
, "%4s has bonds to", aname_i
);
258 for (j
= 0; (j
< nb
); j
++)
260 fprintf(debug
, " %4s", *atoms
->atomname
[bbb
[j
]]);
262 fprintf(debug
, "\n");
265 for (k
= 0; (k
< ematchNR
); k
++)
270 /* First check for names */
271 for (k
= 0; (k
< nnm
); k
++)
273 if (nm2t
[k
].nbonds
== nb
)
275 im
= match_str(*atoms
->atomname
[i
], nm2t
[k
].elem
);
278 for (j
= 0; (j
< ematchNR
); j
++)
283 /* Fill a matrix with matching quality */
284 for (m
= 0; (m
< nb
); m
++)
286 aname_m
= *atoms
->atomname
[bbb
[m
]];
287 for (n
= 0; (n
< nb
); n
++)
289 aname_n
= nm2t
[k
].bond
[n
];
290 match
[m
][n
] = match_str(aname_m
, aname_n
);
293 /* Now pick the best matches */
294 for (m
= 0; (m
< nb
); m
++)
299 for (j
= ematchNR
-1; (j
> 0); j
--)
301 for (m
= 0; (m
< nb
); m
++)
303 for (n
= 0; (n
< nb
); n
++)
305 if ((n_mask
[n
] == 0) &&
316 if ((nqual
[cur
][ematchExact
]+
317 nqual
[cur
][ematchElem
]+
318 nqual
[cur
][ematchWild
]) == nb
)
320 if ((nqual
[cur
][ematchExact
] > nqual
[prev
][ematchExact
]) ||
322 ((nqual
[cur
][ematchExact
] == nqual
[prev
][ematchExact
]) &&
323 (nqual
[cur
][ematchElem
] > nqual
[prev
][ematchElem
])) ||
325 ((nqual
[cur
][ematchExact
] == nqual
[prev
][ematchExact
]) &&
326 (nqual
[cur
][ematchElem
] == nqual
[prev
][ematchElem
]) &&
327 (nqual
[cur
][ematchWild
] > nqual
[prev
][ematchWild
])))
343 type
= nm2t
[best
].type
;
345 if ((k
= get_atomtype_type(type
, atype
)) == NOTSET
)
347 atoms
->atom
[i
].qB
= alpha
;
348 atoms
->atom
[i
].m
= atoms
->atom
[i
].mB
= mm
;
349 k
= add_atomtype(atype
, tab
, &(atoms
->atom
[i
]), type
, param
,
350 atoms
->atom
[i
].type
, 0, 0, 0, atomnr
, 0, 0);
352 atoms
->atom
[i
].type
= k
;
353 atoms
->atom
[i
].typeB
= k
;
354 atoms
->atom
[i
].q
= qq
;
355 atoms
->atom
[i
].qB
= qq
;
356 atoms
->atom
[i
].m
= mm
;
357 atoms
->atom
[i
].mB
= mm
;
362 fprintf(stderr
, "Can not find forcefield for atom %s-%d with %d bonds\n",
363 *atoms
->atomname
[i
], i
+1, nb
);