src/gromacs/gmxpreprocess/toputil.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  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.
  14
  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.
  19  *
  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.
  26  *
  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.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "smalloc.h"
  41 #include "sysstuff.h"
  42 #include "macros.h"
  43 #include "string2.h"
  44 #include "topdirs.h"
  45 #include "toputil.h"
  46 #include "topdirs.h"
  47 #include "toputil.h"
  48 #include "symtab.h"
  49 #include "gmx_fatal.h"
  50 #include "gpp_atomtype.h"
  51
  52 /* UTILITIES */
  53
  54 void set_p_string(t_param *p,const char *s)
  55 {
  56   if (s) {
  57     if (strlen(s) < sizeof(p->s)-1)
  58       strncpy(p->s,s,sizeof(p->s));
  59     else
  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);
  63   }
  64   else
  65     strcpy(p->s,"");
  66 }
  67
  68 void pr_alloc (int extra, t_params *pr)
  69 {
  70   int i,j;
  71
  72   /* get new space for arrays */
  73   if (extra < 0)
  74     gmx_fatal(FARGS,"Trying to make array smaller.\n");
  75   if (extra == 0)
  76     return;
  77   if ((pr->nr == 0) && (pr->param != NULL)) {
  78     fprintf(stderr,"Warning: dangling pointer at %lx\n",
  79             (unsigned long)pr->param);
  80     pr->param = NULL;
  81   }
  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++)
  87         pr->param[i].a[j]=0;
  88       for(j=0; (j<MAXFORCEPARAM); j++)
  89         pr->param[i].c[j]=0;
  90       set_p_string(&(pr->param[i]),"");
  91     }
  92   }
  93 }
  94
  95 void init_plist(t_params plist[])
  96 {
  97   int i;
  98
  99   for(i=0; (i<F_NRE); i++) {
 100     plist[i].nr    = 0;
 101     plist[i].maxnr = 0;
 102     plist[i].param = NULL;
 103
 104     /* CMAP */
 105     plist[i].ncmap=0;
 106     plist[i].cmap=NULL;
 107     plist[i].grid_spacing = 0;
 108     plist[i].nc = 0;
 109     plist[i].nct        = 0;
 110     plist[i].cmap_types = NULL;
 111   }
 112 }
 113
 114 void cp_param(t_param *dest,t_param *src)
 115 {
 116   int j;
 117
 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));
 123 }
 124
 125 void add_param_to_list(t_params *list, t_param *b)
 126 {
 127   int j;
 128
 129   /* allocate one position extra */
 130   pr_alloc (1,list);
 131
 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));
 138
 139   list->nr++;
 140 }
 141
 142
 143 void init_molinfo(t_molinfo *mol)
 144 {
 145   mol->nrexcl = 0;
 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);
 153 }
 154
 155 /* FREEING MEMORY */
 156
 157 void done_bt (t_params *pl)
 158 {
 159   sfree(pl->param);
 160 }
 161
 162 void done_mi(t_molinfo *mi)
 163 {
 164   int i;
 165
 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]));
 171 }
 172
 173 /* PRINTING STRUCTURES */
 174
 175 void print_bt(FILE *out, directive d, gpp_atomtype_t at,
 176               int ftype,int fsubtype,t_params plist[],
 177               gmx_bool bFullDih)
 178 {
 179   /* This dihp is a DIRTY patch because the dih-types do not use
 180    * all four atoms to determine the type.
 181    */
 182   const int dihp[2][2] = { { 1,2 }, { 0,3 } };
 183   t_params *bt;
 184   int      i,j,f,nral,nrfp;
 185   gmx_bool     bDih=FALSE,bSwapParity;
 186
 187   bt=&(plist[ftype]);
 188
 189   if (!bt->nr)
 190     return;
 191
 192   f = 0;
 193   switch(ftype) {
 194   case F_G96ANGLES:
 195     f = 1;
 196     break;
 197   case F_G96BONDS:
 198     f = 1;
 199     break;
 200   case F_MORSE:
 201     f = 2;
 202     break;
 203   case F_CUBICBONDS:
 204     f = 3;
 205     break;
 206   case F_CONNBONDS:
 207     f = 4;
 208     break;
 209   case F_HARMONIC:
 210     f = 5;
 211     break;
 212   case F_CROSS_BOND_ANGLES:
 213     f = 2;
 214     break;
 215   case F_CROSS_BOND_BONDS:
 216     f = 3;
 217     break;
 218   case F_UREY_BRADLEY:
 219     f = 4;
 220     break;
 221   case F_PDIHS:
 222   case F_RBDIHS:
 223   case F_FOURDIHS:
 224     bDih=TRUE;
 225     break;
 226   case F_IDIHS:
 227     f=1;
 228     bDih=TRUE;
 229     break;
 230   case F_CONSTRNC:
 231     f=1;
 232     break;
 233   case F_VSITE3FD:
 234     f = 1;
 235     break;
 236   case F_VSITE3FAD:
 237     f = 2;
 238     break;
 239   case F_VSITE3OUT:
 240     f = 3;
 241     break;
 242   case F_VSITE4FDN:
 243       f = 1;
 244     break;
 245   case F_CMAP:
 246           f = 1;
 247         break;
 248
 249   default:
 250     bDih=FALSE;
 251   }
 252   if (bFullDih)
 253     bDih=FALSE;
 254   if (fsubtype)
 255     f = fsubtype-1;
 256
 257   nral = NRAL(ftype);
 258   nrfp = NRFP(ftype);
 259
 260   /* header */
 261   fprintf(out,"[ %s ]\n",dir2str(d));
 262   fprintf(out,"; ");
 263   if (!bDih) {
 264     fprintf (out,"%3s  %4s","ai","aj");
 265     for (j=2; (j<nral); j++)
 266       fprintf (out,"  %3c%c",'a','i'+j);
 267   }
 268   else
 269     for (j=0; (j<2); j++)
 270       fprintf (out,"%3c%c",'a','i'+dihp[f][j]);
 271
 272   fprintf (out," funct");
 273   for (j=0; (j<nrfp); j++)
 274     fprintf (out," %12c%1d",'c',j);
 275   fprintf (out,"\n");
 276
 277   /* print bondtypes */
 278   for (i=0; (i<bt->nr); i++) {
 279     bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
 280     if (!bDih)
 281       for (j=0; (j<nral); j++)
 282         fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[j],at));
 283     else
 284       for(j=0; (j<2); j++)
 285         fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[dihp[f][j]],at));
 286     fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
 287
 288     if (bt->param[i].s[0])
 289       fprintf(out,"   %s",bt->param[i].s);
 290     else
 291       for (j=0; (j<nrfp && (bt->param[i].c[j] != NOTSET)); j++)
 292         fprintf (out,"%13.6e ",bt->param[i].c[j]);
 293
 294     fprintf (out,"\n");
 295   }
 296   fprintf (out,"\n");
 297   fflush (out);
 298 }
 299
 300 void print_blocka(FILE *out, const char *szName,
 301                   const char *szIndex, const char *szA,
 302                   t_blocka *block)
 303 {
 304   int i,j;
 305
 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);
 313       fprintf (out,"\n");
 314     }
 315     fprintf (out,"\n");
 316   }
 317 }
 318
 319 void print_excl(FILE *out, int natoms, t_excls excls[])
 320 {
 321   atom_id i;
 322   gmx_bool    have_excl;
 323   int     j;
 324
 325   have_excl=FALSE;
 326   for(i=0; i<natoms && !have_excl; i++)
 327     have_excl = (excls[i].nr > 0);
 328
 329   if (have_excl) {
 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);
 337         fprintf (out,"\n");
 338       }
 339     fprintf (out,"\n");
 340     fflush(out);
 341   }
 342 }
 343
 344 static double get_residue_charge(const t_atoms *atoms,int at)
 345 {
 346   int ri;
 347   double q;
 348
 349   ri = atoms->atom[at].resind;
 350   q = 0;
 351   while (at < atoms->nr && atoms->atom[at].resind == ri) {
 352     q += atoms->atom[at].q;
 353     at++;
 354   }
 355
 356   return q;
 357 }
 358
 359 void print_atoms(FILE *out,gpp_atomtype_t atype,t_atoms *at,int *cgnr,
 360                  gmx_bool bRTPresname)
 361 {
 362   int  i,ri;
 363   int  tpA,tpB;
 364   const char *as;
 365   char *tpnmA,*tpnmB;
 366   double qres,qtot;
 367
 368   as=dir2str(d_atoms);
 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");
 372
 373   qtot  = 0;
 374
 375   if (debug)
 376     fprintf(debug,"This molecule has %d atoms and %d residues\n",
 377             at->nr,at->nres);
 378
 379   if (at->nres) {
 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 ",
 387                 at->resinfo[ri].nr,
 388                 *at->resinfo[ri].name,
 389                 *at->resinfo[ri].rtp);
 390         if (fabs(qres) < 0.001) {
 391           fprintf(out," %s","0.0");
 392         } else {
 393           fprintf(out,"%+3.1f",qres);
 394         }
 395         fprintf(out,"\n");
 396       }
 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);
 400
 401       fprintf(out,"%6d %10s %6d%c %5s %6s %6d %10g %10g",
 402               i+1,tpnmA,
 403               at->resinfo[ri].nr,
 404               at->resinfo[ri].ic,
 405               bRTPresname ?
 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);
 416       }
 417       qtot += (double)at->atom[i].q;
 418       if ( fabs(qtot) < 4*GMX_REAL_EPS )
 419         qtot=0;
 420       fprintf(out,"   ; qtot %.4g\n",qtot);
 421     }
 422   }
 423   fprintf(out,"\n");
 424   fflush(out);
 425 }
 426
 427 void print_bondeds(FILE *out,int natoms,directive d,
 428                    int ftype,int fsubtype,t_params plist[])
 429 {
 430   t_symtab   stab;
 431   gpp_atomtype_t atype;
 432   t_param    *param;
 433   t_atom     *a;
 434   int i;
 435
 436   atype = init_atomtype();
 437   snew(a,1);
 438   snew(param,1);
 439   open_symtab(&stab);
 440   for (i=0; (i < natoms); i++) {
 441     char buf[12];
 442     sprintf(buf,"%4d",(i+1));
 443     add_atomtype(atype,&stab,a,buf,param,0,0,0,0,0,0,0);
 444   }
 445   print_bt(out,d,atype,ftype,fsubtype,plist,TRUE);
 446
 447   done_symtab(&stab);
 448   sfree(a);
 449   sfree(param);
 450   done_atomtype(atype);
 451 }
 452