src/kernel/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 static void print_nbt (FILE *out,char *title,gpp_atomtype_t at,
 176                        int ftype,t_params *nbt)
 177 {
 178   int f,i,j,k,l,nrfp,ntype;
 179
 180   if (ftype == F_LJ)
 181     f=1;
 182   else
 183     f=2;
 184   nrfp=NRFP(ftype);
 185
 186   if (nbt->nr) {
 187     /* header */
 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);
 194     fprintf (out,"\n");
 195
 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]);
 204         fprintf (out,"\n");
 205       }
 206   }
 207   fprintf (out,"\n");
 208 }
 209
 210 void print_bt(FILE *out, directive d, gpp_atomtype_t at,
 211               int ftype,int fsubtype,t_params plist[],
 212               bool bFullDih)
 213 {
 214   /* This dihp is a DIRTY patch because the dih-types do not use
 215    * all four atoms to determine the type.
 216    */
 217   const int dihp[2][2] = { { 1,2 }, { 0,3 } };
 218   t_params *bt;
 219   int      i,j,f,nral,nrfp;
 220   bool     bDih=FALSE,bSwapParity;
 221
 222   bt=&(plist[ftype]);
 223
 224   if (!bt->nr)
 225     return;
 226
 227   f = 0;
 228   switch(ftype) {
 229   case F_G96ANGLES:
 230     f = 1;
 231     break;
 232   case F_G96BONDS:
 233     f = 1;
 234     break;
 235   case F_MORSE:
 236     f = 2;
 237     break;
 238   case F_CUBICBONDS:
 239     f = 3;
 240     break;
 241   case F_CONNBONDS:
 242     f = 4;
 243     break;
 244   case F_HARMONIC:
 245     f = 5;
 246     break;
 247   case F_CROSS_BOND_ANGLES:
 248     f = 2;
 249     break;
 250   case F_CROSS_BOND_BONDS:
 251     f = 3;
 252     break;
 253   case F_UREY_BRADLEY:
 254     f = 4;
 255     break;
 256   case F_PDIHS:
 257   case F_RBDIHS:
 258   case F_FOURDIHS:
 259     bDih=TRUE;
 260     break;
 261   case F_IDIHS:
 262     f=1;
 263     bDih=TRUE;
 264     break;
 265   case F_CONSTRNC:
 266     f=1;
 267     break;
 268   case F_VSITE3FD:
 269     f = 1;
 270     break;
 271   case F_VSITE3FAD:
 272     f = 2;
 273     break;
 274   case F_VSITE3OUT:
 275     f = 3;
 276     break;
 277   case F_VSITE4FDN:
 278       f = 1;
 279     break;
 280   case F_CMAP:
 281           f = 1;
 282         break;
 283
 284   default:
 285     bDih=FALSE;
 286   }
 287   if (bFullDih)
 288     bDih=FALSE;
 289   if (fsubtype)
 290     f = fsubtype-1;
 291
 292   nral = NRAL(ftype);
 293   nrfp = NRFP(ftype);
 294
 295   /* header */
 296   fprintf(out,"[ %s ]\n",dir2str(d));
 297   fprintf(out,"; ");
 298   if (!bDih) {
 299     fprintf (out,"%3s  %4s","ai","aj");
 300     for (j=2; (j<nral); j++)
 301       fprintf (out,"  %3c%c",'a','i'+j);
 302   }
 303   else
 304     for (j=0; (j<2); j++)
 305       fprintf (out,"%3c%c",'a','i'+dihp[f][j]);
 306
 307   fprintf (out," funct");
 308   for (j=0; (j<nrfp); j++)
 309     fprintf (out," %12c%1d",'c',j);
 310   fprintf (out,"\n");
 311
 312   /* print bondtypes */
 313   for (i=0; (i<bt->nr); i++) {
 314     bSwapParity = (bt->param[i].C0==NOTSET) && (bt->param[i].C1==-1);
 315     if (!bDih)
 316       for (j=0; (j<nral); j++)
 317         fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[j],at));
 318     else
 319       for(j=0; (j<2); j++)
 320         fprintf (out,"%5s ",get_atomtype_name(bt->param[i].a[dihp[f][j]],at));
 321     fprintf (out,"%5d ", bSwapParity ? -f-1 : f+1);
 322
 323     if (bt->param[i].s[0])
 324       fprintf(out,"   %s",bt->param[i].s);
 325     else
 326       for (j=0; (j<nrfp && (bt->param[i].c[j] != NOTSET)); j++)
 327         fprintf (out,"%13.6e ",bt->param[i].c[j]);
 328
 329     fprintf (out,"\n");
 330   }
 331   fprintf (out,"\n");
 332   fflush (out);
 333 }
 334
 335 void print_blocka(FILE *out, const char *szName,
 336                   const char *szIndex, const char *szA,
 337                   t_blocka *block)
 338 {
 339   int i,j;
 340
 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);
 348       fprintf (out,"\n");
 349     }
 350     fprintf (out,"\n");
 351   }
 352 }
 353
 354 void print_excl(FILE *out, int natoms, t_excls excls[])
 355 {
 356   atom_id i;
 357   bool    have_excl;
 358   int     j;
 359
 360   have_excl=FALSE;
 361   for(i=0; i<natoms && !have_excl; i++)
 362     have_excl = (excls[i].nr > 0);
 363
 364   if (have_excl) {
 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);
 372         fprintf (out,"\n");
 373       }
 374     fprintf (out,"\n");
 375     fflush(out);
 376   }
 377 }
 378
 379 static double get_residue_charge(const t_atoms *atoms,int at)
 380 {
 381   int ri;
 382   double q;
 383
 384   ri = atoms->atom[at].resind;
 385   q = 0;
 386   while (at < atoms->nr && atoms->atom[at].resind == ri) {
 387     q += atoms->atom[at].q;
 388     at++;
 389   }
 390
 391   return q;
 392 }
 393
 394 void print_atoms(FILE *out,gpp_atomtype_t atype,t_atoms *at,int *cgnr,
 395                  bool bRTPresname)
 396 {
 397   int  i,ri;
 398   int  tpA,tpB;
 399   const char *as;
 400   char *tpnmA,*tpnmB;
 401   double qres,qtot;
 402
 403   as=dir2str(d_atoms);
 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");
 407
 408   qtot  = 0;
 409
 410   if (debug)
 411     fprintf(debug,"This molecule has %d atoms and %d residues\n",
 412             at->nr,at->nres);
 413
 414   if (at->nres) {
 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 ",
 422                 at->resinfo[ri].nr,
 423                 *at->resinfo[ri].name,
 424                 *at->resinfo[ri].rtp);
 425         if (fabs(qres) < 0.001) {
 426           fprintf(out," %s","0.0");
 427         } else {
 428           fprintf(out,"%+3.1f",qres);
 429         }
 430         fprintf(out,"\n");
 431       }
 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);
 435
 436       fprintf(out,"%6d %10s %6d%c %5s %6s %6d %10g %10g",
 437               i+1,tpnmA,
 438               at->resinfo[ri].nr,
 439               at->resinfo[ri].ic,
 440               bRTPresname ?
 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);
 451       }
 452       qtot += (double)at->atom[i].q;
 453       if ( fabs(qtot) < 4*GMX_REAL_EPS )
 454         qtot=0;
 455       fprintf(out,"   ; qtot %.4g\n",qtot);
 456     }
 457   }
 458   fprintf(out,"\n");
 459   fflush(out);
 460 }
 461
 462 void print_bondeds(FILE *out,int natoms,directive d,
 463                    int ftype,int fsubtype,t_params plist[])
 464 {
 465   t_symtab   stab;
 466   gpp_atomtype_t atype;
 467   t_param    *param;
 468   t_atom     *a;
 469   int i;
 470
 471   atype = init_atomtype();
 472   snew(a,1);
 473   snew(param,1);
 474   open_symtab(&stab);
 475   for (i=0; (i < natoms); i++) {
 476     char buf[12];
 477     sprintf(buf,"%4d",(i+1));
 478     add_atomtype(atype,&stab,a,buf,param,0,0,0,0,0,0,0);
 479   }
 480   print_bt(out,d,atype,ftype,fsubtype,plist,TRUE);
 481
 482   done_symtab(&stab);
 483   sfree(a);
 484   sfree(param);
 485   done_atomtype(atype);
 486 }
 487