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
];