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"
50 #include "gpp_atomtype.h"
54 void set_p_string(t_param *p,const char *s)
57 if (strlen(s) < sizeof(p->s)-1)
58 strncpy(p->s,s,sizeof(p->s));
60 gmx_fatal(FARGS,"Increase MAXSLEN in include/grompp.h to at least %d,"
61 " or shorten your definition of bonds like %s to at most %d",
62 strlen(s)+1,s,MAXSLEN-1);
68 void pr_alloc (int extra, t_params *pr)
72 /* get new space for arrays */
74 gmx_fatal(FARGS,"Trying to make array smaller.\n");
77 if ((pr->nr == 0) && (pr->param != NULL)) {
78 fprintf(stderr,"Warning: dangling pointer at %lx\n",
79 (unsigned long)pr->param);
82 if (pr->nr+extra > pr->maxnr) {
83 pr->maxnr = max(1.2*pr->maxnr,pr->maxnr + extra);
84 srenew(pr->param,pr->maxnr);
85 for(i=pr->nr; (i<pr->maxnr); i++) {
86 for(j=0; (j<MAXATOMLIST); j++)
88 for(j=0; (j<MAXFORCEPARAM); j++)
90 set_p_string(&(pr->param[i]),"");
95 void init_plist(t_params plist[])
99 for(i=0; (i<F_NRE); i++) {
102 plist[i].param = NULL;
107 plist[i].grid_spacing = 0;
110 plist[i].cmap_types = NULL;
114 void cp_param(t_param *dest,t_param *src)
118 for(j=0; (j<MAXATOMLIST); j++)
119 dest->a[j] = src->a[j];
120 for(j=0; (j<MAXFORCEPARAM); j++)
121 dest->c[j] = src->c[j];
122 strncpy(dest->s,src->s,sizeof(dest->s));
125 void add_param_to_list(t_params *list, t_param *b)
129 /* allocate one position extra */
132 /* fill the arrays */
133 for (j=0; (j < MAXFORCEPARAM); j++)
134 list->param[list->nr].c[j] = b->c[j];
135 for (j=0; (j < MAXATOMLIST); j++)
136 list->param[list->nr].a[j] = b->a[j];
137 memset(list->param[list->nr].s,0,sizeof(list->param[list->nr].s));
143 void init_molinfo(t_molinfo *mol)
146 mol->excl_set = FALSE;
147 mol->bProcessed = FALSE;
148 init_plist(mol->plist);
149 init_block(&mol->cgs);
150 init_block(&mol->mols);
151 init_blocka(&mol->excls);
152 init_atom(&mol->atoms);
157 void done_bt (t_params *pl)
162 void done_mi(t_molinfo *mi)
166 done_atom (&(mi->atoms));
167 done_block(&(mi->cgs));
168 done_block(&(mi->mols));
169 for(i=0; (i<F_NRE); i++)
170 done_bt(&(mi->plist[i]));
173 /* PRINTING STRUCTURES */
175 void print_bt(FILE *out, directive d, gpp_atomtype_t at,
176 int ftype,int fsubtype,t_params plist[],
179 /* This dihp is a DIRTY patch because the dih-types do not use
180 * all four atoms to determine the type.
182 const int dihp[2][2] = { { 1,2 }, { 0,3 } };
185 gmx_bool bDih=FALSE,bSwapParity;
212 case F_CROSS_BOND_ANGLES:
215 case F_CROSS_BOND_BONDS:
261 fprintf(out,"[ %s ]\n",dir2str(d));
264 fprintf (out,"%3s %4s","ai","aj");
265 for (j=2; (j<nral); j++)
266 fprintf (out," %3c%c",'a','i'+j);
269 for (j=0; (j<2); j++)
270 fprintf (out,"%3c%c",'a','i'+dihp[f][j]);
272 fprintf (out," funct");
273 for (j=0; (j<nrfp); j++)
274 fprintf (out," %12c%1d",'c',j);
277 /* print bondtypes */
278 for (i=0; (i<bt->nr); i++) {
279 bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
281 for (j=0; (j<nral); j++)
282 fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[j],at));
285 fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[dihp[f][j]],at));
286 fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
288 if (bt->param[i].s[0])
289 fprintf(out," %s",bt->param[i].s);
291 for (j=0; (j<nrfp && (bt->param[i].c[j] != NOTSET)); j++)
292 fprintf (out,"%13.6e ",bt->param[i].c[j]);
300 void print_blocka(FILE *out, const char *szName,
301 const char *szIndex, const char *szA,
306 fprintf (out,"; %s\n",szName);
307 fprintf (out,"; %4s %s\n",szIndex,szA);
308 for (i=0; (i < block->nr); i++) {
309 for (i=0; (i < block->nr); i++) {
310 fprintf (out,"%6d",i+1);
311 for (j=block->index[i]; (j < ((int)block->index[i+1])); j++)
312 fprintf (out,"%5d",block->a[j]+1);
319 void print_excl(FILE *out, int natoms, t_excls excls[])
326 for(i=0; i<natoms && !have_excl; i++)
327 have_excl = (excls[i].nr > 0);
330 fprintf (out,"[ %s ]\n",dir2str(d_exclusions));
331 fprintf (out,"; %4s %s\n","i","excluded from i");
332 for(i=0; i<natoms; i++)
333 if (excls[i].nr > 0) {
334 fprintf (out,"%6d ",i+1);
335 for(j=0; j<excls[i].nr; j++)
336 fprintf (out," %5d",excls[i].e[j]+1);
344 static double get_residue_charge(const t_atoms *atoms,int at)
349 ri = atoms->atom[at].resind;
351 while (at < atoms->nr && atoms->atom[at].resind == ri) {
352 q += atoms->atom[at].q;
359 void print_atoms(FILE *out,gpp_atomtype_t atype,t_atoms *at,int *cgnr,
360 gmx_bool bRTPresname)
369 fprintf(out,"[ %s ]\n",as);
370 fprintf(out,"; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
371 "nr","type","resnr","residue","atom","cgnr","charge","mass","typeB","chargeB","massB");
376 fprintf(debug,"This molecule has %d atoms and %d residues\n",
380 /* if the information is present... */
381 for (i=0; (i < at->nr); i++) {
382 ri = at->atom[i].resind;
383 if ((i == 0 || ri != at->atom[i-1].resind) &&
384 at->resinfo[ri].rtp != NULL) {
385 qres = get_residue_charge(at,i);
386 fprintf(out,"; residue %3d %-3s rtp %-4s q ",
388 *at->resinfo[ri].name,
389 *at->resinfo[ri].rtp);
390 if (fabs(qres) < 0.001) {
391 fprintf(out," %s","0.0");
393 fprintf(out,"%+3.1f",qres);
397 tpA = at->atom[i].type;
398 if ((tpnmA = get_atomtype_name(tpA,atype)) == NULL)
399 gmx_fatal(FARGS,"tpA = %d, i= %d in print_atoms",tpA,i);
401 fprintf(out,"%6d %10s %6d%c %5s %6s %6d %10g %10g",
406 *(at->resinfo[at->atom[i].resind].rtp) :
407 *(at->resinfo[at->atom[i].resind].name),
408 *(at->atomname[i]),cgnr[i],
409 at->atom[i].q,at->atom[i].m);
410 if (PERTURBED(at->atom[i])) {
411 tpB = at->atom[i].typeB;
412 if ((tpnmB = get_atomtype_name(tpB,atype)) == NULL)
413 gmx_fatal(FARGS,"tpB = %d, i= %d in print_atoms",tpB,i);
414 fprintf(out," %6s %10g %10g",
415 tpnmB,at->atom[i].qB,at->atom[i].mB);
417 qtot += (double)at->atom[i].q;
418 if ( fabs(qtot) < 4*GMX_REAL_EPS )
420 fprintf(out," ; qtot %.4g\n",qtot);
427 void print_bondeds(FILE *out,int natoms,directive d,
428 int ftype,int fsubtype,t_params plist[])
431 gpp_atomtype_t atype;
436 atype = init_atomtype();
440 for (i=0; (i < natoms); i++) {
442 sprintf(buf,"%4d",(i+1));
443 add_atomtype(atype,&stab,a,buf,param,0,0,0,0,0,0,0);
445 print_bt(out,d,atype,ftype,fsubtype,plist,TRUE);
450 done_atomtype(atype);