src/kernel/toputil.c

   1 /*
   2  * $Id$
   3  *
   4  *       This source code is part of
   5  *
   6  *        G   R   O   M   A   C   S
   7  *
   8  * GROningen MAchine for Chemical Simulations
   9  *
  10  *               VERSION 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * GRowing Old MAkes el Chrono Sweat
  28  */
  29 static char *SRCID_toputil_c = "$Id$";
  30
  31 #ifdef HAVE_IDENT
  32 #ident "@(#) toputil.c 1.68 9/30/97"
  33 #endif
  34
  35 #include "assert.h"
  36 #include "smalloc.h"
  37 #include "sysstuff.h"
  38 #include "macros.h"
  39 #include "string2.h"
  40 #include "topdirs.h"
  41 #include "toputil.h"
  42 #include "topdirs.h"
  43 #include "toputil.h"
  44 #include "symtab.h"
  45 #include "fatal.h"
  46
  47 /* UTILITIES */
  48
  49 int at2type(char *str,t_atomtype *at)
  50 {
  51   int i;
  52
  53   for (i=0; (i<at->nr) && strcasecmp(str,*(at->atomname[i])); i++)
  54     ;
  55   if (i == at->nr)
  56     fatal_error(0,"Atomtype '%s' not found!",str);
  57
  58   return i;
  59 }
  60
  61 char *type2nm(int nt, t_atomtype *at)
  62 {
  63   if ((nt < 0) || (nt >= at->nr))
  64     fatal_error(0,"nt out of range: %d",nt);
  65
  66   return *(at->atomname[nt]);
  67 }
  68
  69 void pr_alloc (int extra, t_params *pr)
  70 {
  71   int i,j;
  72
  73   /* get new space for arrays */
  74   if (extra < 0)
  75     fatal_error(0,"Trying to make array < 0 bytes\n");
  76   if (extra == 0)
  77     return;
  78   if ((pr->nr == 0) && (pr->param != NULL)) {
  79     fprintf(stderr,"Warning: dangling pointer at %x\n",
  80             (unsigned int)pr->param);
  81     pr->param = NULL;
  82   }
  83   srenew(pr->param,pr->nr+extra);
  84   for(i=pr->nr; (i<pr->nr+extra); i++) {
  85     for(j=0; (j<MAXATOMLIST); j++)
  86       pr->param[i].a[j]=0;
  87     for(j=0; (j<MAXFORCEPARAM); j++)
  88       pr->param[i].c[j]=0;
  89     pr->param[i].s = NULL;
  90   }
  91 }
  92
  93 /* INIT STRUCTURES */
  94
  95 void init_atomtype (t_atomtype *at)
  96 {
  97   at->nr   = 0;
  98   at->atom     = NULL;
  99   at->atomname = NULL;
 100   at->nb       = NULL;
 101 }
 102
 103 void init_plist(t_params plist[])
 104 {
 105   int i;
 106
 107   for(i=0; (i<F_NRE); i++) {
 108     plist[i].nr    = 0;
 109     plist[i].param = NULL;
 110   }
 111 }
 112
 113 void init_molinfo(t_molinfo *mol)
 114 {
 115   mol->nrexcl = 0;
 116   mol->excl_set = FALSE;
 117   mol->bProcessed = FALSE;
 118   init_plist(mol->plist);
 119   init_block(&mol->cgs);
 120   init_block(&mol->mols);
 121   init_atom(&mol->atoms);
 122 }
 123
 124 /* FREEING MEMORY */
 125
 126 void done_bt (t_params *pl)
 127 {
 128   sfree(pl->param);
 129 }
 130
 131 void done_mi(t_molinfo *mi)
 132 {
 133   int i;
 134
 135   done_atom (&(mi->atoms));
 136   done_block(&(mi->cgs));
 137   done_block(&(mi->mols));
 138   for(i=0; (i<F_NRE); i++)
 139     done_bt(&(mi->plist[i]));
 140 }
 141
 142 /* PRINTING STRUCTURES */
 143
 144 void print_at (FILE * out, t_atomtype *at)
 145 {
 146   int i;
 147   t_atom     *atom = at->atom;
 148   t_param    *nb   = at->nb;
 149
 150   fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
 151   fprintf (out,"; %6s  %8s  %8s  %8s  %12s  %12s\n",
 152            "type","mass","charge","particle","c6","c12");
 153   for (i=0; (i<at->nr); i++)
 154     fprintf(out,"%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
 155             *(at->atomname[i]),atom[i].m,atom[i].q,"A",
 156             nb[i].C0,nb[i].C1);
 157
 158   fprintf (out,"\n");
 159 }
 160
 161 static void print_nbt (FILE *out,char *title,t_atomtype *at,
 162                        int ftype,t_params *nbt)
 163 {
 164   int f,i,j,k,l,nrfp;
 165
 166   if (ftype == F_LJ)
 167     f=1;
 168   else
 169     f=2;
 170   nrfp=NRFP(ftype);
 171
 172   if (nbt->nr) {
 173     /* header */
 174     fprintf (out,"; %s\n",title);
 175     fprintf (out,"[ %s ]\n",dir2str(d_nonbond_params));
 176     fprintf (out,"; %6s  %8s","ai","aj");
 177     fprintf (out,"  %8s","funct");
 178     for (j=0; (j<nrfp); j++)
 179       fprintf (out,"  %11c%1d",'c',j);
 180     fprintf (out,"\n");
 181
 182     /* print non-bondtypes */
 183     for (i=k=0; (i<at->nr); i++)
 184       for(j=0; (j<at->nr); j++,k++) {
 185         fprintf (out,"%8s  %8s  %8d",
 186                  *(at->atomname[i]),*(at->atomname[j]),f);
 187         for(l=0; (l<nrfp); l++)
 188           fprintf (out,"  %12.5e",nbt->param[k].c[l]);
 189         fprintf (out,"\n");
 190       }
 191   }
 192   fprintf (out,"\n");
 193 }
 194
 195 void print_bt(FILE *out, directive d, t_atomtype *at,
 196               int ftype,int fsubtype,t_params plist[],
 197               bool bFullDih)
 198 {
 199   /* This dihp is a DIRTY patch because the dih-types do not use
 200    * all four atoms to determine the type.
 201    */
 202   static int dihp[2][2] = { { 1,2 }, { 0,3 } };
 203   t_params *bt;
 204   int      i,j,f,nral,nrfp;
 205   bool     bDih=FALSE,bSwapParity;
 206
 207   bt=&(plist[ftype]);
 208
 209   if (!bt->nr)
 210     return;
 211
 212   f = 0;
 213   switch(ftype) {
 214   case F_G96ANGLES:
 215   case F_G96BONDS:
 216     f = 1;
 217     f = 1;
 218     break;
 219   case F_MORSE:
 220     f = 2;
 221     break;
 222   case F_PDIHS:
 223   case F_RBDIHS:
 224     bDih=TRUE;
 225     break;
 226   case F_IDIHS:
 227     f=1;
 228     bDih=TRUE;
 229     break;
 230   case F_SHAKENC:
 231     f=1;
 232     break;
 233   case F_DUMMY3FD:
 234     f = 1;
 235     break;
 236   case F_DUMMY3FAD:
 237     f = 2;
 238     break;
 239   case F_DUMMY3OUT:
 240     f = 3;
 241     break;
 242   default:
 243     bDih=FALSE;
 244   }
 245   if (bFullDih)
 246     bDih=FALSE;
 247   if (fsubtype)
 248     f = fsubtype-1;
 249
 250   nral = NRAL(ftype);
 251   nrfp = NRFP(ftype);
 252
 253   /* header */
 254   fprintf(out,"[ %s ]\n",dir2str(d));
 255   fprintf(out,"; ");
 256   if (!bDih) {
 257     fprintf (out,"%3s  %4s","ai","aj");
 258     for (j=2; (j<nral); j++)
 259       fprintf (out,"  %3c%c",'a','i'+j);
 260   }
 261   else
 262     for (j=0; (j<2); j++)
 263       fprintf (out,"%3c%c",'a','i'+dihp[f][j]);
 264
 265   fprintf (out," funct");
 266   for (j=0; (j<nrfp); j++)
 267     fprintf (out," %12c%1d",'c',j);
 268   fprintf (out,"\n");
 269
 270   /* print bondtypes */
 271   for (i=0; (i<bt->nr); i++) {
 272     bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
 273     if (!bDih)
 274       for (j=0; (j<nral); j++)
 275         fprintf (out,"%5s ",*(at->atomname[bt->param[i].a[j]]));
 276     else
 277       for(j=0; (j<2); j++)
 278         fprintf (out,"%5s ",*(at->atomname[bt->param[i].a[dihp[f][j]]]));
 279     fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
 280
 281     if (bt->param[i].s[0])
 282       fprintf(out,"   %s",bt->param[i].s);
 283     else
 284       for (j=0; (j<nrfp && (bt->param[i].c[j] != NOTSET)); j++)
 285         fprintf (out,"%13.6e ",bt->param[i].c[j]);
 286
 287     fprintf (out,"\n");
 288   }
 289   fprintf (out,"\n");
 290   fflush (out);
 291 }
 292
 293 void print_block (FILE *out, char *szName,
 294                   char *szIndex, char *szA,
 295                   t_block *block)
 296 {
 297   int i,j;
 298
 299   fprintf (out,"; %s\n",szName);
 300   fprintf (out,"; %4s    %s\n",szIndex,szA);
 301   for (i=0; (i < block->nr); i++) {
 302     for (i=0; (i < block->nr); i++) {
 303       fprintf (out,"%6d",i+1);
 304       for (j=block->index[i]; (j < ((int)block->index[i+1])); j++)
 305         fprintf (out,"%5u",block->a[j]+1);
 306       fprintf (out,"\n");
 307     }
 308     fprintf (out,"\n");
 309   }
 310 }
 311
 312 void print_excl(FILE *out, t_block *excl)
 313 {
 314   int     i;
 315   bool    have_excl;
 316   atom_id j;
 317
 318   have_excl=FALSE;
 319   if (excl->nr)
 320     for (i=0; (i < excl->nr) && !have_excl; i++)
 321       if (excl->index[i] != excl->index[i+1])
 322         have_excl=TRUE;
 323
 324   if (have_excl) {
 325     fprintf (out,"[ %s ]\n",dir2str(d_exclusions));
 326     fprintf (out,"; %4s    %s\n","i","excluded from i");
 327     for (i=0; (i < excl->nr); i++)
 328       if (excl->index[i] != excl->index[i+1]) {
 329         fprintf (out,"%6d",i+1);
 330         for (j=excl->index[i]; (j < excl->index[i+1]); j++)
 331           fprintf (out,"%5u",excl->a[j]+1);
 332         fprintf (out,"\n");
 333       }
 334     fprintf (out,"\n");
 335     fflush(out);
 336   }
 337 }
 338
 339 void print_atoms(FILE *out,t_atomtype *atype,t_atoms *at,int *cgnr)
 340 {
 341   int  i;
 342   int  itype;
 343   char *as;
 344   real qtot;
 345
 346   as=dir2str(d_atoms);
 347   fprintf(out,"[ %s ]\n",as);
 348   fprintf(out,"; %4s %6s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
 349           "nr","type","resnr","residue","atom","cgnr","charge","mass","typeB","chargeB","massB");
 350
 351   qtot  = 0;
 352
 353   if (debug)
 354     fprintf(debug,"This molecule has %d atoms and %d residues\n",
 355             at->nr,at->nres);
 356
 357   if (at->nres) {
 358     /* if the information is present... */
 359     for (i=0; (i < at->nr); i++) {
 360       itype=at->atom[i].type;
 361       if ((itype < 0) || (itype > atype->nr))
 362         fatal_error(0,"itype = %d, i= %d in print_atoms",itype,i);
 363
 364       fprintf(out,"%6d %6s %6d %6s %6s %6d %10g %10g",
 365               i+1,*(atype->atomname[itype]),
 366               at->atom[i].resnr+1,
 367               *(at->resname[at->atom[i].resnr]),
 368               *(at->atomname[i]),cgnr[i],
 369               at->atom[i].q,at->atom[i].m);
 370       if (PERTURBED(at->atom[i])) {
 371         fprintf(out," %6s %10g %10g",
 372                 *(atype->atomname[at->atom[i].typeB]),
 373                 at->atom[i].qB,at->atom[i].mB);
 374       }
 375       qtot+=at->atom[i].q;
 376       fprintf(out,"   ; qtot %g\n",qtot);
 377     }
 378   }
 379   fprintf(out,"\n");
 380   fflush(out);
 381 }
 382
 383 void print_bondeds(FILE *out,int natoms,directive d,
 384                    int ftype,int fsubtype,t_params plist[])
 385 {
 386   t_symtab   stab;
 387   t_atomtype atype;
 388   int i;
 389
 390   snew(atype.atom,natoms);
 391   snew(atype.atomname,natoms);
 392   open_symtab(&stab);
 393   for (i=0; (i < natoms); i++) {
 394     char buf[12];
 395     sprintf(buf,"%4d",(i+1));
 396     atype.atomname[i]=put_symtab(&stab,buf);
 397   }
 398   print_bt(out,d,&atype,ftype,fsubtype,plist,TRUE);
 399
 400   done_symtab(&stab);
 401   sfree(atype.atom);
 402   sfree(atype.atomname);
 403 }
 404