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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /* This file is completely threadsafe - keep it that way! */
62 #include "gpp_nextnb.h"
63 #include "gpp_atomtype.h"
65 #include "fflibutil.h"
67 static void rd_nm2type_file(const char *fn,int *nnm,t_nm2type **nmp)
71 char libfilename[128];
72 char format[128],f1[128];
73 char buf[1024],elem[16],type[16],nbbuf[16],**newbuf;
80 gmx_fatal(FARGS,"Can not find %s in library directory",fn);
85 /* Read a line from the file */
86 bCont = (fgets2(buf,1023,fp) != NULL);
91 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);
97 strcpy(format,"%*s%*s%*s%*s%*s");
98 for(i=0; (i<nb); i++) {
99 /* Complicated format statement */
102 if (sscanf(buf,f1,nbbuf,&(nm2t[nnnm].blen[i])) != 2)
103 gmx_fatal(FARGS,"Error on line %d of %s",line,libfilename);
104 newbuf[i] = strdup(nbbuf);
105 strcat(format,"%*s%*s");
110 nm2t[nnnm].elem = strdup(elem);
111 nm2t[nnnm].type = strdup(type);
114 nm2t[nnnm].nbonds = nb;
115 nm2t[nnnm].bond = newbuf;
127 t_nm2type *rd_nm2type(const char *ffdir,int *nnm)
133 nff = fflib_search_file_end(ffdir,".n2t",FALSE,&ff);
136 for(f=0; f<nff; f++) {
137 rd_nm2type_file(ff[f],nnm,&nm);
145 void dump_nm2type(FILE *fp,int nnm,t_nm2type nm2t[])
149 fprintf(fp,"; nm2type database\n");
150 for(i=0; (i<nnm); i++) {
151 fprintf(fp,"%-8s %-8s %8.4f %8.4f %-4d",
152 nm2t[i].elem,nm2t[i].type,
153 nm2t[i].q,nm2t[i].m,nm2t[i].nbonds);
154 for(j=0; (j<nm2t[i].nbonds); j++)
155 fprintf(fp," %-5s %6.4f",nm2t[i].bond[j],nm2t[i].blen[j]);
160 enum { ematchNone, ematchWild, ematchElem, ematchExact, ematchNR };
162 static int match_str(const char *atom,const char *template_string)
164 if (!atom || !template_string)
166 else if (gmx_strcasecmp(atom,template_string) == 0)
168 else if (atom[0] == template_string[0])
170 else if (strcmp(template_string,"*") == 0)
176 int nm2type(int nnm,t_nm2type nm2t[],t_symtab *tab,t_atoms *atoms,
177 gpp_atomtype_t atype,int *nbonds,t_params *bonds)
181 int i,j,k,m,n,nresolved,nb,maxbond,ai,aj,best,im,nqual[2][ematchNR];
182 int *bbb,*n_mask,*m_mask,**match,**quality;
183 char *aname_i,*aname_m,*aname_n,*type;
189 for(i=0; (i<atoms->nr); i++)
190 maxbond = max(maxbond,nbonds[i]);
192 fprintf(debug,"Max number of bonds per atom is %d\n",maxbond);
194 snew(n_mask,maxbond);
195 snew(m_mask,maxbond);
197 for(i=0; (i<maxbond); i++)
198 snew(match[i],maxbond);
201 for(i=0; (i<atoms->nr); i++) {
202 aname_i = *atoms->atomname[i];
204 for(j=0; (j<bonds->nr); j++) {
205 ai = bonds->param[j].AI;
206 aj = bonds->param[j].AJ;
213 gmx_fatal(FARGS,"Counting number of bonds nb = %d, nbonds[%d] = %d",
216 fprintf(debug,"%4s has bonds to",aname_i);
217 for(j=0; (j<nb); j++)
218 fprintf(debug," %4s",*atoms->atomname[bbb[j]]);
222 for(k=0; (k<ematchNR); k++)
225 /* First check for names */
226 for(k=0; (k<nnm); k++) {
227 if (nm2t[k].nbonds == nb) {
228 im = match_str(*atoms->atomname[i],nm2t[k].elem);
229 if (im > ematchWild) {
230 for(j=0; (j<ematchNR); j++)
233 /* Fill a matrix with matching quality */
234 for(m=0; (m<nb); m++) {
235 aname_m = *atoms->atomname[bbb[m]];
236 for(n=0; (n<nb); n++) {
237 aname_n = nm2t[k].bond[n];
238 match[m][n] = match_str(aname_m,aname_n);
241 /* Now pick the best matches */
242 for(m=0; (m<nb); m++) {
246 for(j=ematchNR-1; (j>0); j--) {
247 for(m=0; (m<nb); m++) {
248 for(n=0; (n<nb); n++) {
249 if ((n_mask[n] == 0) &&
251 (match[m][n] == j)) {
259 if ((nqual[cur][ematchExact]+
260 nqual[cur][ematchElem]+
261 nqual[cur][ematchWild]) == nb) {
262 if ((nqual[cur][ematchExact] > nqual[prev][ematchExact]) ||
264 ((nqual[cur][ematchExact] == nqual[prev][ematchExact]) &&
265 (nqual[cur][ematchElem] > nqual[prev][ematchElem])) ||
267 ((nqual[cur][ematchExact] == nqual[prev][ematchExact]) &&
268 (nqual[cur][ematchElem] == nqual[prev][ematchElem]) &&
269 (nqual[cur][ematchWild] > nqual[prev][ematchWild]))) {
283 type = nm2t[best].type;
285 if ((k = get_atomtype_type(type,atype)) == NOTSET) {
286 atoms->atom[i].qB = alpha;
287 atoms->atom[i].m = atoms->atom[i].mB = mm;
288 k = add_atomtype(atype,tab,&(atoms->atom[i]),type,param,
289 atoms->atom[i].type,0,0,0,atomnr,0,0);
291 atoms->atom[i].type = k;
292 atoms->atom[i].typeB = k;
293 atoms->atom[i].q = qq;
294 atoms->atom[i].qB = qq;
295 atoms->atom[i].m = mm;
296 atoms->atom[i].mB = mm;
300 fprintf(stderr,"Can not find forcefield for atom %s-%d with %d bonds\n",
301 *atoms->atomname[i],i+1,nb);