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 static void print_nbt (FILE *out,char *title,gpp_atomtype_t at,
176 int ftype,t_params *nbt)
178 int f,i,j,k,l,nrfp,ntype;
188 fprintf (out,"; %s\n",title);
189 fprintf (out,"[ %s ]\n",dir2str(d_nonbond_params));
190 fprintf (out,"; %6s %8s","ai","aj");
191 fprintf (out," %8s","funct");
192 for (j=0; (j<nrfp); j++)
193 fprintf (out," %11c%1d",'c',j);
196 /* print non-bondtypes */
197 ntype = get_atomtype_ntypes(at);
198 for (i=k=0; (i<ntype); i++)
199 for(j=0; (j<ntype); j++,k++) {
200 fprintf (out,"%8s %8s %8d",
201 get_atomtype_name(i,at),get_atomtype_name(f,at),f);
202 for(l=0; (l<nrfp); l++)
203 fprintf (out," %12.5e",nbt->param[k].c[l]);
210 void print_bt(FILE *out, directive d, gpp_atomtype_t at,
211 int ftype,int fsubtype,t_params plist[],
214 /* This dihp is a DIRTY patch because the dih-types do not use
215 * all four atoms to determine the type.
217 const int dihp[2][2] = { { 1,2 }, { 0,3 } };
220 bool bDih=FALSE,bSwapParity;
247 case F_CROSS_BOND_ANGLES:
250 case F_CROSS_BOND_BONDS:
296 fprintf(out,"[ %s ]\n",dir2str(d));
299 fprintf (out,"%3s %4s","ai","aj");
300 for (j=2; (j<nral); j++)
301 fprintf (out," %3c%c",'a','i'+j);
304 for (j=0; (j<2); j++)
305 fprintf (out,"%3c%c",'a','i'+dihp[f][j]);
307 fprintf (out," funct");
308 for (j=0; (j<nrfp); j++)
309 fprintf (out," %12c%1d",'c',j);
312 /* print bondtypes */
313 for (i=0; (i<bt->nr); i++) {
314 bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
316 for (j=0; (j<nral); j++)
317 fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[j],at));
320 fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[dihp[f][j]],at));
321 fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
323 if (bt->param[i].s[0])
324 fprintf(out," %s",bt->param[i].s);
326 for (j=0; (j<nrfp && (bt->param[i].c[j] != NOTSET)); j++)
327 fprintf (out,"%13.6e ",bt->param[i].c[j]);
335 void print_blocka(FILE *out, const char *szName,
336 const char *szIndex, const char *szA,
341 fprintf (out,"; %s\n",szName);
342 fprintf (out,"; %4s %s\n",szIndex,szA);
343 for (i=0; (i < block->nr); i++) {
344 for (i=0; (i < block->nr); i++) {
345 fprintf (out,"%6d",i+1);
346 for (j=block->index[i]; (j < ((int)block->index[i+1])); j++)
347 fprintf (out,"%5d",block->a[j]+1);
354 void print_excl(FILE *out, int natoms, t_excls excls[])
361 for(i=0; i<natoms && !have_excl; i++)
362 have_excl = (excls[i].nr > 0);
365 fprintf (out,"[ %s ]\n",dir2str(d_exclusions));
366 fprintf (out,"; %4s %s\n","i","excluded from i");
367 for(i=0; i<natoms; i++)
368 if (excls[i].nr > 0) {
369 fprintf (out,"%6d ",i+1);
370 for(j=0; j<excls[i].nr; j++)
371 fprintf (out," %5d",excls[i].e[j]+1);
379 static double get_residue_charge(const t_atoms *atoms,int at)
384 ri = atoms->atom[at].resind;
386 while (at < atoms->nr && atoms->atom[at].resind == ri) {
387 q += atoms->atom[at].q;
394 void print_atoms(FILE *out,gpp_atomtype_t atype,t_atoms *at,int *cgnr,
404 fprintf(out,"[ %s ]\n",as);
405 fprintf(out,"; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
406 "nr","type","resnr","residue","atom","cgnr","charge","mass","typeB","chargeB","massB");
411 fprintf(debug,"This molecule has %d atoms and %d residues\n",
415 /* if the information is present... */
416 for (i=0; (i < at->nr); i++) {
417 ri = at->atom[i].resind;
418 if ((i == 0 || ri != at->atom[i-1].resind) &&
419 at->resinfo[ri].rtp != NULL) {
420 qres = get_residue_charge(at,i);
421 fprintf(out,"; residue %3d %-3s rtp %-4s q ",
423 *at->resinfo[ri].name,
424 *at->resinfo[ri].rtp);
425 if (fabs(qres) < 0.001) {
426 fprintf(out," %s","0.0");
428 fprintf(out,"%+3.1f",qres);
432 tpA = at->atom[i].type;
433 if ((tpnmA = get_atomtype_name(tpA,atype)) == NULL)
434 gmx_fatal(FARGS,"tpA = %d, i= %d in print_atoms",tpA,i);
436 fprintf(out,"%6d %10s %6d%c %5s %6s %6d %10g %10g",
441 *(at->resinfo[at->atom[i].resind].rtp) :
442 *(at->resinfo[at->atom[i].resind].name),
443 *(at->atomname[i]),cgnr[i],
444 at->atom[i].q,at->atom[i].m);
445 if (PERTURBED(at->atom[i])) {
446 tpB = at->atom[i].typeB;
447 if ((tpnmB = get_atomtype_name(tpB,atype)) == NULL)
448 gmx_fatal(FARGS,"tpB = %d, i= %d in print_atoms",tpB,i);
449 fprintf(out," %6s %10g %10g",
450 tpnmB,at->atom[i].qB,at->atom[i].mB);
452 qtot += (double)at->atom[i].q;
453 if ( fabs(qtot) < 4*GMX_REAL_EPS )
455 fprintf(out," ; qtot %.4g\n",qtot);
462 void print_bondeds(FILE *out,int natoms,directive d,
463 int ftype,int fsubtype,t_params plist[])
466 gpp_atomtype_t atype;
471 atype = init_atomtype();
475 for (i=0; (i < natoms); i++) {
477 sprintf(buf,"%4d",(i+1));
478 add_atomtype(atype,&stab,a,buf,param,0,0,0,0,0,0,0);
480 print_bt(out,d,atype,ftype,fsubtype,plist,TRUE);
485 done_atomtype(atype);