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-2004, 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 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
49 #include "gmx_fatal.h"
51 #include "gpp_atomtype.h"
53 typedef struct gpp_atomtype {
54 int nr; /* The number of atomtypes */
55 t_atom *atom; /* Array of atoms */
56 char ***atomname; /* Names of the atomtypes */
57 t_param *nb; /* Nonbonded force default params */
58 int *bondatomtype; /* The bond_atomtype for each atomtype */
59 real *radius; /* Radius for GBSA stuff */
60 real *vol; /* Effective volume for GBSA */
61 real *surftens; /* Surface tension with water, for GBSA */
62 real *gb_radius; /* Radius for Still model */
63 real *S_hct; /* Overlap factor for HCT model */
64 int *atomnumber; /* Atomic number, used for QM/MM */
67 int get_atomtype_type(const char *str,gpp_atomtype_t ga)
71 /* Atom types are always case sensitive */
72 for (i=0; (i<ga->nr); i++)
73 if (strcmp(str,*(ga->atomname[i])) == 0)
79 int get_atomtype_ntypes(gpp_atomtype_t ga)
84 char *get_atomtype_name(int nt, gpp_atomtype_t ga)
86 if ((nt < 0) || (nt >= ga->nr))
89 return *(ga->atomname[nt]);
92 real get_atomtype_massA(int nt,gpp_atomtype_t ga)
94 if ((nt < 0) || (nt >= ga->nr))
97 return ga->atom[nt].m;
100 real get_atomtype_massB(int nt,gpp_atomtype_t ga)
102 if ((nt < 0) || (nt >= ga->nr))
105 return ga->atom[nt].mB;
108 real get_atomtype_qA(int nt,gpp_atomtype_t ga)
110 if ((nt < 0) || (nt >= ga->nr))
113 return ga->atom[nt].q;
116 real get_atomtype_qB(int nt,gpp_atomtype_t ga)
118 if ((nt < 0) || (nt >= ga->nr))
121 return ga->atom[nt].qB;
124 int get_atomtype_ptype(int nt,gpp_atomtype_t ga)
126 if ((nt < 0) || (nt >= ga->nr))
129 return ga->atom[nt].ptype;
132 int get_atomtype_batype(int nt,gpp_atomtype_t ga)
134 if ((nt < 0) || (nt >= ga->nr))
137 return ga->bondatomtype[nt];
140 int get_atomtype_atomnumber(int nt,gpp_atomtype_t ga)
142 if ((nt < 0) || (nt >= ga->nr))
145 return ga->atomnumber[nt];
148 real get_atomtype_radius(int nt,gpp_atomtype_t ga)
150 if ((nt < 0) || (nt >= ga->nr))
153 return ga->radius[nt];
156 real get_atomtype_vol(int nt,gpp_atomtype_t ga)
158 if ((nt < 0) || (nt >= ga->nr))
164 real get_atomtype_surftens(int nt,gpp_atomtype_t ga)
166 if ((nt < 0) || (nt >= ga->nr))
169 return ga->surftens[nt];
172 real get_atomtype_gb_radius(int nt,gpp_atomtype_t ga)
174 if ((nt < 0) || (nt >= ga->nr))
177 return ga->gb_radius[nt];
180 real get_atomtype_S_hct(int nt,gpp_atomtype_t ga)
182 if ((nt < 0) || (nt >= ga->nr))
185 return ga->S_hct[nt];
188 real get_atomtype_nbparam(int nt,int param,gpp_atomtype_t ga)
190 if ((nt < 0) || (nt >= ga->nr))
192 if ((param < 0) || (param >= MAXFORCEPARAM))
194 return ga->nb[nt].c[param];
197 gpp_atomtype_t init_atomtype(void)
207 ga->bondatomtype = NULL;
211 ga->atomnumber = NULL;
212 ga->gb_radius = NULL;
219 set_atomtype_gbparam(gpp_atomtype_t ga, int i,
220 real radius,real vol,real surftens,
221 real gb_radius, real S_hct)
223 if ( (i < 0) || (i>= ga->nr))
226 ga->radius[i] = radius;
228 ga->surftens[i] = surftens;
229 ga->gb_radius[i] = gb_radius;
230 ga->S_hct[i] = S_hct;
236 int set_atomtype(int nt,gpp_atomtype_t ga,t_symtab *tab,
237 t_atom *a,const char *name,t_param *nb,
239 real radius,real vol,real surftens,int atomnumber,
240 real gb_radius, real S_hct)
242 if ((nt < 0) || (nt >= ga->nr))
246 ga->atomname[nt] = put_symtab(tab,name);
248 ga->bondatomtype[nt] = bondatomtype;
249 ga->radius[nt] = radius;
251 ga->surftens[nt] = surftens;
252 ga->atomnumber[nt] = atomnumber;
253 ga->gb_radius[nt] = gb_radius;
254 ga->S_hct[nt] = S_hct;
259 int add_atomtype(gpp_atomtype_t ga,t_symtab *tab,
260 t_atom *a,const char *name,t_param *nb,
262 real radius,real vol,real surftens,real atomnumber,
263 real gb_radius, real S_hct)
267 for(i=0; (i<ga->nr); i++) {
268 if (strcmp(*ga->atomname[i],name) == 0) {
270 fprintf(debug,"Trying to add atomtype %s again. Skipping it.\n",name);
276 srenew(ga->atom,ga->nr);
277 srenew(ga->atomname,ga->nr);
278 srenew(ga->nb,ga->nr);
279 srenew(ga->bondatomtype,ga->nr);
280 srenew(ga->radius,ga->nr);
281 srenew(ga->vol,ga->nr);
282 srenew(ga->surftens,ga->nr);
283 srenew(ga->atomnumber,ga->nr);
284 srenew(ga->gb_radius,ga->nr);
285 srenew(ga->S_hct,ga->nr);
287 return set_atomtype(ga->nr-1,ga,tab,a,name,nb,bondatomtype,radius,
288 vol,surftens,atomnumber,gb_radius,S_hct);
294 void print_at (FILE * out, gpp_atomtype_t ga)
297 t_atom *atom = ga->atom;
298 t_param *nb = ga->nb;
300 fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
301 fprintf (out,"; %6s %8s %8s %8s %12s %12s\n",
302 "type","mass","charge","particle","c6","c12");
303 for (i=0; (i<ga->nr); i++)
304 fprintf(out,"%8s %8.3f %8.3f %8s %12e %12e\n",
305 *(ga->atomname[i]),atom[i].m,atom[i].q,"A",
311 void done_atomtype(gpp_atomtype_t ga)
316 sfree(ga->bondatomtype);
319 sfree(ga->gb_radius);
322 sfree(ga->atomnumber);
328 static int search_atomtypes(gpp_atomtype_t ga,int *n,int typelist[],
330 t_param param[],int ftype)
332 int i,nn,nrfp,j,k,ntype,tli;
333 gmx_bool bFound=FALSE;
337 ntype = get_atomtype_ntypes(ga);
339 for(i=0; (i<nn); i++)
341 if (typelist[i] == thistype)
343 /* This type number has already been added */
347 /* Otherwise, check if the parameters are identical to any previously added type */
350 for(j=0;j<ntype && bFound;j++)
352 /* Check nonbonded parameters */
353 for(k=0;k<nrfp && bFound;k++)
355 bFound=(param[ntype*typelist[i]+j].c[k]==param[ntype*thistype+j].c[k]);
358 /* Check radius, volume, surftens */
361 (get_atomtype_radius(tli,ga) == get_atomtype_radius(thistype,ga)) &&
362 (get_atomtype_vol(tli,ga) == get_atomtype_vol(thistype,ga)) &&
363 (get_atomtype_surftens(tli,ga) == get_atomtype_surftens(thistype,ga)) &&
364 (get_atomtype_atomnumber(tli,ga) == get_atomtype_atomnumber(thistype,ga)) &&
365 (get_atomtype_gb_radius(tli,ga) == get_atomtype_gb_radius(thistype,ga)) &&
366 (get_atomtype_S_hct(tli,ga) == get_atomtype_S_hct(thistype,ga));
376 fprintf(debug,"Renumbering atomtype %d to %d\n",thistype,nn);
378 gmx_fatal(FARGS,"Atomtype horror n = %d, %s, %d",nn,__FILE__,__LINE__);
379 typelist[nn]=thistype;
387 void renum_atype(t_params plist[],gmx_mtop_t *mtop,
389 gpp_atomtype_t ga,gmx_bool bVerbose)
391 int i,j,k,l,molt,mi,mj,nat,nrfp,ftype,ntype;
401 char ***new_atomname;
403 ntype = get_atomtype_ntypes(ga);
404 snew(typelist,ntype);
407 fprintf(stderr,"renumbering atomtypes...\n");
409 /* Since the bonded interactions have been assigned now,
410 * we want to reduce the number of atom types by merging
411 * ones with identical nonbonded interactions, in addition
412 * to removing unused ones.
414 * With Generalized-Born electrostatics, or implicit solvent
415 * we also check that the atomtype radius, effective_volume
416 * and surface tension match.
418 * With QM/MM we also check that the atom numbers match
421 /* Get nonbonded interaction type */
422 if (plist[F_LJ].nr > 0)
427 /* Renumber atomtypes by first making a list of which ones are actually used.
428 * We provide the list of nonbonded parameters so search_atomtypes
429 * can determine if two types should be merged.
432 for(molt=0; molt<mtop->nmoltype; molt++) {
433 atoms = &mtop->moltype[molt].atoms;
434 for(i=0; (i<atoms->nr); i++) {
435 atoms->atom[i].type =
436 search_atomtypes(ga,&nat,typelist,atoms->atom[i].type,
437 plist[ftype].param,ftype);
438 atoms->atom[i].typeB=
439 search_atomtypes(ga,&nat,typelist,atoms->atom[i].typeB,
440 plist[ftype].param,ftype);
445 if (wall_atomtype[i] >= 0)
446 wall_atomtype[i] = search_atomtypes(ga,&nat,typelist,wall_atomtype[i],
447 plist[ftype].param,ftype);
450 snew(new_radius,nat);
452 snew(new_surftens,nat);
453 snew(new_atomnumber,nat);
454 snew(new_gb_radius,nat);
456 snew(new_atomname,nat);
458 /* We now have a list of unique atomtypes in typelist */
461 pr_ivec(debug,0,"typelist",typelist,nat,TRUE);
465 snew(nbsnew,plist[ftype].nr);
469 for(i=k=0; (i<nat); i++)
472 for(j=0; (j<nat); j++,k++)
475 for(l=0; (l<nrfp); l++)
477 nbsnew[k].c[l]=plist[ftype].param[ntype*mi+mj].c[l];
480 new_radius[i] = get_atomtype_radius(mi,ga);
481 new_vol[i] = get_atomtype_vol(mi,ga);
482 new_surftens[i] = get_atomtype_surftens(mi,ga);
483 new_atomnumber[i] = get_atomtype_atomnumber(mi,ga);
484 new_gb_radius[i] = get_atomtype_gb_radius(mi,ga);
485 new_S_hct[i] = get_atomtype_S_hct(mi,ga);
486 new_atomname[i] = ga->atomname[mi];
489 for(i=0; (i<nat*nat); i++) {
490 for(l=0; (l<nrfp); l++)
491 plist[ftype].param[i].c[l]=nbsnew[i].c[l];
494 mtop->ffparams.atnr = nat;
499 sfree(ga->atomnumber);
500 sfree(ga->gb_radius);
502 /* Dangling atomname pointers ? */
505 ga->radius = new_radius;
507 ga->surftens = new_surftens;
508 ga->atomnumber = new_atomnumber;
509 ga->gb_radius = new_gb_radius;
510 ga->S_hct = new_S_hct;
511 ga->atomname = new_atomname;
519 void copy_atomtype_atomtypes(gpp_atomtype_t ga,t_atomtypes *atomtypes)
523 /* Copy the atomtype data to the topology atomtype list */
524 ntype = get_atomtype_ntypes(ga);
526 snew(atomtypes->radius,ntype);
527 snew(atomtypes->vol,ntype);
528 snew(atomtypes->surftens,ntype);
529 snew(atomtypes->atomnumber,ntype);
530 snew(atomtypes->gb_radius,ntype);
531 snew(atomtypes->S_hct,ntype);
533 for(i=0; i<ntype; i++) {
534 atomtypes->radius[i] = ga->radius[i];
535 atomtypes->vol[i] = ga->vol[i];
536 atomtypes->surftens[i] = ga->surftens[i];
537 atomtypes->atomnumber[i] = ga->atomnumber[i];
538 atomtypes->gb_radius[i] = ga->gb_radius[i];
539 atomtypes->S_hct[i] = ga->S_hct[i];