3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2008, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
35 /* This file is completely threadsafe - keep it that way! */
59 #include "gpp_nextnb.h"
60 #include "gpp_atomtype.h"
62 #include "fflibutil.h"
64 static void rd_nm2type_file(const char *fn,int *nnm,t_nm2type **nmp)
68 char libfilename[128];
69 char format[128],f1[128];
70 char buf[1024],elem[16],type[16],nbbuf[16],**newbuf;
77 gmx_fatal(FARGS,"Can not find %s in library directory",fn);
82 /* Read a line from the file */
83 bCont = (fgets2(buf,1023,fp) != NULL);
88 if (sscanf(buf,"%s%s%lf%lf%d",elem,type,&qq,&mm,&nb) == 5) {
89 /* If we can read the first four, there probably is more */
91 snew(nm2t[nnnm].blen,nb);
94 strcpy(format,"%*s%*s%*s%*s%*s");
95 for(i=0; (i<nb); i++) {
96 /* Complicated format statement */
99 if (sscanf(buf,f1,nbbuf,&(nm2t[nnnm].blen[i])) != 2)
100 gmx_fatal(FARGS,"Error on line %d of %s",line,libfilename);
101 newbuf[i] = strdup(nbbuf);
102 strcat(format,"%*s%*s");
107 nm2t[nnnm].elem = strdup(elem);
108 nm2t[nnnm].type = strdup(type);
111 nm2t[nnnm].nbonds = nb;
112 nm2t[nnnm].bond = newbuf;
124 t_nm2type *rd_nm2type(const char *ffdir,int *nnm)
130 nff = fflib_search_file_end(ffdir,".n2t",FALSE,&ff);
133 for(f=0; f<nff; f++) {
134 rd_nm2type_file(ff[f],nnm,&nm);
142 void dump_nm2type(FILE *fp,int nnm,t_nm2type nm2t[])
146 fprintf(fp,"; nm2type database\n");
147 for(i=0; (i<nnm); i++) {
148 fprintf(fp,"%-8s %-8s %8.4f %8.4f %-4d",
149 nm2t[i].elem,nm2t[i].type,
150 nm2t[i].q,nm2t[i].m,nm2t[i].nbonds);
151 for(j=0; (j<nm2t[i].nbonds); j++)
152 fprintf(fp," %-5s %6.4f",nm2t[i].bond[j],nm2t[i].blen[j]);
157 enum { ematchNone, ematchWild, ematchElem, ematchExact, ematchNR };
159 static int match_str(const char *atom,const char *template_string)
161 if (!atom || !template_string)
163 else if (gmx_strcasecmp(atom,template_string) == 0)
165 else if (atom[0] == template_string[0])
167 else if (strcmp(template_string,"*") == 0)
173 int nm2type(int nnm,t_nm2type nm2t[],t_symtab *tab,t_atoms *atoms,
174 gpp_atomtype_t atype,int *nbonds,t_params *bonds)
178 int i,j,k,m,n,nresolved,nb,maxbond,ai,aj,best,im,nqual[2][ematchNR];
179 int *bbb,*n_mask,*m_mask,**match,**quality;
180 char *aname_i,*aname_m,*aname_n,*type;
186 for(i=0; (i<atoms->nr); i++)
187 maxbond = max(maxbond,nbonds[i]);
189 fprintf(debug,"Max number of bonds per atom is %d\n",maxbond);
191 snew(n_mask,maxbond);
192 snew(m_mask,maxbond);
194 for(i=0; (i<maxbond); i++)
195 snew(match[i],maxbond);
198 for(i=0; (i<atoms->nr); i++) {
199 aname_i = *atoms->atomname[i];
201 for(j=0; (j<bonds->nr); j++) {
202 ai = bonds->param[j].AI;
203 aj = bonds->param[j].AJ;
210 gmx_fatal(FARGS,"Counting number of bonds nb = %d, nbonds[%d] = %d",
213 fprintf(debug,"%4s has bonds to",aname_i);
214 for(j=0; (j<nb); j++)
215 fprintf(debug," %4s",*atoms->atomname[bbb[j]]);
219 for(k=0; (k<ematchNR); k++)
222 /* First check for names */
223 for(k=0; (k<nnm); k++) {
224 if (nm2t[k].nbonds == nb) {
225 im = match_str(*atoms->atomname[i],nm2t[k].elem);
226 if (im > ematchWild) {
227 for(j=0; (j<ematchNR); j++)
230 /* Fill a matrix with matching quality */
231 for(m=0; (m<nb); m++) {
232 aname_m = *atoms->atomname[bbb[m]];
233 for(n=0; (n<nb); n++) {
234 aname_n = nm2t[k].bond[n];
235 match[m][n] = match_str(aname_m,aname_n);
238 /* Now pick the best matches */
239 for(m=0; (m<nb); m++) {
243 for(j=ematchNR-1; (j>0); j--) {
244 for(m=0; (m<nb); m++) {
245 for(n=0; (n<nb); n++) {
246 if ((n_mask[n] == 0) &&
248 (match[m][n] == j)) {
256 if ((nqual[cur][ematchExact]+
257 nqual[cur][ematchElem]+
258 nqual[cur][ematchWild]) == nb) {
259 if ((nqual[cur][ematchExact] > nqual[prev][ematchExact]) ||
261 ((nqual[cur][ematchExact] == nqual[prev][ematchExact]) &&
262 (nqual[cur][ematchElem] > nqual[prev][ematchElem])) ||
264 ((nqual[cur][ematchExact] == nqual[prev][ematchExact]) &&
265 (nqual[cur][ematchElem] == nqual[prev][ematchElem]) &&
266 (nqual[cur][ematchWild] > nqual[prev][ematchWild]))) {
280 type = nm2t[best].type;
282 if ((k = get_atomtype_type(type,atype)) == NOTSET) {
283 atoms->atom[i].qB = alpha;
284 atoms->atom[i].m = atoms->atom[i].mB = mm;
285 k = add_atomtype(atype,tab,&(atoms->atom[i]),type,param,
286 atoms->atom[i].type,0,0,0,atomnr,0,0);
288 atoms->atom[i].type = k;
289 atoms->atom[i].typeB = k;
290 atoms->atom[i].q = qq;
291 atoms->atom[i].qB = qq;
292 atoms->atom[i].m = mm;
293 atoms->atom[i].mB = mm;
297 fprintf(stderr,"Can not find forcefield for atom %s-%d with %d bonds\n",
298 *atoms->atomname[i],i+1,nb);