src/gromacs/gmxpreprocess/gpp_atomtype.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 "txtdump.h"
  51 #include "gpp_atomtype.h"
  52
  53 typedef struct gpp_atomtype {
  54   int           nr;             /* The number of atomtypes              */
  55   t_atom        *atom;          /* Array of atoms                       */
  56   char          ***atomname;    /* Names of the atomtypes               */
  57   t_param       *nb;            /* Nonbonded force default params       */
  58   int           *bondatomtype;  /* The bond_atomtype for each atomtype  */
  59   real          *radius;        /* Radius for GBSA stuff                */
  60   real          *vol;           /* Effective volume for GBSA            */
  61   real          *surftens;      /* Surface tension with water, for GBSA */
  62   real          *gb_radius;     /* Radius for Still model               */
  63   real          *S_hct;         /* Overlap factor for HCT model         */
  64   int           *atomnumber;    /* Atomic number, used for QM/MM        */
  65 } t_gpp_atomtype;
  66
  67 int get_atomtype_type(const char *str,gpp_atomtype_t ga)
  68 {
  69   int i;
  70
  71   /* Atom types are always case sensitive */
  72   for (i=0; (i<ga->nr); i++)
  73     if (strcmp(str,*(ga->atomname[i])) == 0)
  74       return i;
  75
  76   return NOTSET;
  77 }
  78
  79 int get_atomtype_ntypes(gpp_atomtype_t ga)
  80 {
  81   return ga->nr;
  82 }
  83
  84 char *get_atomtype_name(int nt, gpp_atomtype_t ga)
  85 {
  86   if ((nt < 0) || (nt >= ga->nr))
  87     return NULL;
  88
  89   return *(ga->atomname[nt]);
  90 }
  91
  92 real get_atomtype_massA(int nt,gpp_atomtype_t ga)
  93 {
  94   if ((nt < 0) || (nt >= ga->nr))
  95     return NOTSET;
  96
  97   return ga->atom[nt].m;
  98 }
  99
 100 real get_atomtype_massB(int nt,gpp_atomtype_t ga)
 101 {
 102   if ((nt < 0) || (nt >= ga->nr))
 103     return NOTSET;
 104
 105   return ga->atom[nt].mB;
 106 }
 107
 108 real get_atomtype_qA(int nt,gpp_atomtype_t ga)
 109 {
 110   if ((nt < 0) || (nt >= ga->nr))
 111     return NOTSET;
 112
 113   return ga->atom[nt].q;
 114 }
 115
 116 real get_atomtype_qB(int nt,gpp_atomtype_t ga)
 117 {
 118   if ((nt < 0) || (nt >= ga->nr))
 119     return NOTSET;
 120
 121   return ga->atom[nt].qB;
 122 }
 123
 124 int get_atomtype_ptype(int nt,gpp_atomtype_t ga)
 125 {
 126   if ((nt < 0) || (nt >= ga->nr))
 127     return NOTSET;
 128
 129   return ga->atom[nt].ptype;
 130 }
 131
 132 int get_atomtype_batype(int nt,gpp_atomtype_t ga)
 133 {
 134   if ((nt < 0) || (nt >= ga->nr))
 135     return NOTSET;
 136
 137   return ga->bondatomtype[nt];
 138 }
 139
 140 int get_atomtype_atomnumber(int nt,gpp_atomtype_t ga)
 141 {
 142   if ((nt < 0) || (nt >= ga->nr))
 143     return NOTSET;
 144
 145   return ga->atomnumber[nt];
 146 }
 147
 148 real get_atomtype_radius(int nt,gpp_atomtype_t ga)
 149 {
 150   if ((nt < 0) || (nt >= ga->nr))
 151     return NOTSET;
 152
 153   return ga->radius[nt];
 154 }
 155
 156 real get_atomtype_vol(int nt,gpp_atomtype_t ga)
 157 {
 158   if ((nt < 0) || (nt >= ga->nr))
 159     return NOTSET;
 160
 161   return ga->vol[nt];
 162 }
 163
 164 real get_atomtype_surftens(int nt,gpp_atomtype_t ga)
 165 {
 166   if ((nt < 0) || (nt >= ga->nr))
 167     return NOTSET;
 168
 169   return ga->surftens[nt];
 170 }
 171
 172 real get_atomtype_gb_radius(int nt,gpp_atomtype_t ga)
 173 {
 174   if ((nt < 0) || (nt >= ga->nr))
 175     return NOTSET;
 176
 177   return ga->gb_radius[nt];
 178 }
 179
 180 real get_atomtype_S_hct(int nt,gpp_atomtype_t ga)
 181 {
 182   if ((nt < 0) || (nt >= ga->nr))
 183     return NOTSET;
 184
 185   return ga->S_hct[nt];
 186 }
 187
 188 real get_atomtype_nbparam(int nt,int param,gpp_atomtype_t ga)
 189 {
 190   if ((nt < 0) || (nt >= ga->nr))
 191     return NOTSET;
 192   if ((param < 0) || (param >= MAXFORCEPARAM))
 193     return NOTSET;
 194   return ga->nb[nt].c[param];
 195 }
 196
 197 gpp_atomtype_t init_atomtype(void)
 198 {
 199   gpp_atomtype_t ga;
 200
 201   snew(ga,1);
 202
 203   ga->nr           = 0;
 204   ga->atom         = NULL;
 205   ga->atomname     = NULL;
 206   ga->nb           = NULL;
 207   ga->bondatomtype = NULL;
 208   ga->radius       = NULL;
 209   ga->vol          = NULL;
 210   ga->surftens     = NULL;
 211   ga->atomnumber   = NULL;
 212   ga->gb_radius    = NULL;
 213   ga->S_hct        = NULL;
 214
 215   return ga;
 216 }
 217
 218 int
 219 set_atomtype_gbparam(gpp_atomtype_t ga, int i,
 220                      real radius,real vol,real surftens,
 221                      real gb_radius, real S_hct)
 222 {
 223   if ( (i < 0) || (i>= ga->nr))
 224     return NOTSET;
 225
 226   ga->radius[i]    = radius;
 227   ga->vol[i]       = vol;
 228   ga->surftens[i]  = surftens;
 229   ga->gb_radius[i] = gb_radius;
 230   ga->S_hct[i]     = S_hct;
 231
 232   return i;
 233 }
 234
 235
 236 int set_atomtype(int nt,gpp_atomtype_t ga,t_symtab *tab,
 237                  t_atom *a,const char *name,t_param *nb,
 238                  int bondatomtype,
 239                  real radius,real vol,real surftens,int atomnumber,
 240                  real gb_radius, real S_hct)
 241 {
 242   if ((nt < 0) || (nt >= ga->nr))
 243     return NOTSET;
 244
 245   ga->atom[nt] = *a;
 246   ga->atomname[nt] = put_symtab(tab,name);
 247   ga->nb[nt] = *nb;
 248   ga->bondatomtype[nt] = bondatomtype;
 249   ga->radius[nt] = radius;
 250   ga->vol[nt] = vol;
 251   ga->surftens[nt] = surftens;
 252   ga->atomnumber[nt] = atomnumber;
 253   ga->gb_radius[nt] = gb_radius;
 254   ga->S_hct[nt]     = S_hct;
 255
 256   return nt;
 257 }
 258
 259 int add_atomtype(gpp_atomtype_t ga,t_symtab *tab,
 260                  t_atom *a,const char *name,t_param *nb,
 261                  int bondatomtype,
 262                  real radius,real vol,real surftens,real atomnumber,
 263                  real gb_radius, real S_hct)
 264 {
 265   int i;
 266
 267   for(i=0; (i<ga->nr); i++) {
 268     if (strcmp(*ga->atomname[i],name) == 0) {
 269       if (NULL != debug)
 270         fprintf(debug,"Trying to add atomtype %s again. Skipping it.\n",name);
 271       break;
 272     }
 273   }
 274   if (i == ga->nr) {
 275     ga->nr++;
 276     srenew(ga->atom,ga->nr);
 277     srenew(ga->atomname,ga->nr);
 278     srenew(ga->nb,ga->nr);
 279     srenew(ga->bondatomtype,ga->nr);
 280     srenew(ga->radius,ga->nr);
 281     srenew(ga->vol,ga->nr);
 282     srenew(ga->surftens,ga->nr);
 283     srenew(ga->atomnumber,ga->nr);
 284     srenew(ga->gb_radius,ga->nr);
 285     srenew(ga->S_hct,ga->nr);
 286
 287     return set_atomtype(ga->nr-1,ga,tab,a,name,nb,bondatomtype,radius,
 288                         vol,surftens,atomnumber,gb_radius,S_hct);
 289   }
 290   else
 291     return i;
 292 }
 293
 294 void print_at (FILE * out, gpp_atomtype_t ga)
 295 {
 296   int i;
 297   t_atom     *atom = ga->atom;
 298   t_param    *nb   = ga->nb;
 299
 300   fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
 301   fprintf (out,"; %6s  %8s  %8s  %8s  %12s  %12s\n",
 302            "type","mass","charge","particle","c6","c12");
 303   for (i=0; (i<ga->nr); i++)
 304     fprintf(out,"%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
 305             *(ga->atomname[i]),atom[i].m,atom[i].q,"A",
 306             nb[i].C0,nb[i].C1);
 307
 308   fprintf (out,"\n");
 309 }
 310
 311 void done_atomtype(gpp_atomtype_t ga)
 312 {
 313   sfree(ga->atom);
 314   sfree(ga->atomname);
 315   sfree(ga->nb);
 316   sfree(ga->bondatomtype);
 317   sfree(ga->radius);
 318   sfree(ga->vol);
 319   sfree(ga->gb_radius);
 320   sfree(ga->S_hct);
 321   sfree(ga->surftens);
 322   sfree(ga->atomnumber);
 323   ga->nr = 0;
 324   sfree(ga);
 325   ga = NULL;
 326 }
 327
 328 static int search_atomtypes(gpp_atomtype_t ga,int *n,int typelist[],
 329                             int thistype,
 330                             t_param param[],int ftype)
 331 {
 332   int i,nn,nrfp,j,k,ntype,tli;
 333   gmx_bool bFound=FALSE;
 334
 335   nn    = *n;
 336   nrfp  = NRFP(ftype);
 337   ntype = get_atomtype_ntypes(ga);
 338
 339   for(i=0; (i<nn); i++)
 340   {
 341     if (typelist[i] == thistype)
 342     {
 343       /* This type number has already been added */
 344       break;
 345     }
 346
 347     /* Otherwise, check if the parameters are identical to any previously added type */
 348
 349     bFound=TRUE;
 350     for(j=0;j<ntype && bFound;j++)
 351     {
 352       /* Check nonbonded parameters */
 353       for(k=0;k<nrfp && bFound;k++)
 354       {
 355         bFound=(param[ntype*typelist[i]+j].c[k]==param[ntype*thistype+j].c[k]);
 356       }
 357
 358       /* Check radius, volume, surftens */
 359       tli = typelist[i];
 360       bFound = bFound &&
 361         (get_atomtype_radius(tli,ga) == get_atomtype_radius(thistype,ga)) &&
 362         (get_atomtype_vol(tli,ga) == get_atomtype_vol(thistype,ga)) &&
 363         (get_atomtype_surftens(tli,ga) == get_atomtype_surftens(thistype,ga)) &&
 364         (get_atomtype_atomnumber(tli,ga) == get_atomtype_atomnumber(thistype,ga)) &&
 365         (get_atomtype_gb_radius(tli,ga) == get_atomtype_gb_radius(thistype,ga)) &&
 366         (get_atomtype_S_hct(tli,ga) == get_atomtype_S_hct(thistype,ga));
 367     }
 368     if (bFound)
 369     {
 370       break;
 371     }
 372   }
 373
 374   if (i == nn) {
 375     if (debug)
 376       fprintf(debug,"Renumbering atomtype %d to %d\n",thistype,nn);
 377     if (nn == ntype)
 378       gmx_fatal(FARGS,"Atomtype horror n = %d, %s, %d",nn,__FILE__,__LINE__);
 379     typelist[nn]=thistype;
 380     nn++;
 381   }
 382   *n = nn;
 383
 384   return i;
 385 }
 386
 387 void renum_atype(t_params plist[],gmx_mtop_t *mtop,
 388                 int *wall_atomtype,
 389                 gpp_atomtype_t ga,gmx_bool bVerbose)
 390 {
 391   int      i,j,k,l,molt,mi,mj,nat,nrfp,ftype,ntype;
 392   t_atoms  *atoms;
 393   t_param  *nbsnew;
 394   int      *typelist;
 395   real     *new_radius;
 396   real     *new_vol;
 397   real     *new_surftens;
 398   real     *new_gb_radius;
 399   real     *new_S_hct;
 400   int      *new_atomnumber;
 401   char     ***new_atomname;
 402
 403   ntype = get_atomtype_ntypes(ga);
 404   snew(typelist,ntype);
 405
 406   if (bVerbose)
 407     fprintf(stderr,"renumbering atomtypes...\n");
 408
 409   /* Since the bonded interactions have been assigned now,
 410    * we want to reduce the number of atom types by merging
 411    * ones with identical nonbonded interactions, in addition
 412    * to removing unused ones.
 413    *
 414    * With Generalized-Born electrostatics, or implicit solvent
 415    * we also check that the atomtype radius, effective_volume
 416    * and surface tension match.
 417    *
 418    * With QM/MM we also check that the atom numbers match
 419    */
 420
 421   /* Get nonbonded interaction type */
 422   if (plist[F_LJ].nr > 0)
 423     ftype=F_LJ;
 424   else
 425     ftype=F_BHAM;
 426
 427   /* Renumber atomtypes by first making a list of which ones are actually used.
 428    * We provide the list of nonbonded parameters so search_atomtypes
 429    * can determine if two types should be merged.
 430    */
 431   nat=0;
 432   for(molt=0; molt<mtop->nmoltype; molt++) {
 433     atoms = &mtop->moltype[molt].atoms;
 434     for(i=0; (i<atoms->nr); i++) {
 435       atoms->atom[i].type =
 436         search_atomtypes(ga,&nat,typelist,atoms->atom[i].type,
 437                          plist[ftype].param,ftype);
 438       atoms->atom[i].typeB=
 439         search_atomtypes(ga,&nat,typelist,atoms->atom[i].typeB,
 440                          plist[ftype].param,ftype);
 441     }
 442   }
 443
 444   for(i=0; i<2; i++) {
 445     if (wall_atomtype[i] >= 0)
 446       wall_atomtype[i] = search_atomtypes(ga,&nat,typelist,wall_atomtype[i],
 447                                           plist[ftype].param,ftype);
 448   }
 449
 450   snew(new_radius,nat);
 451   snew(new_vol,nat);
 452   snew(new_surftens,nat);
 453   snew(new_atomnumber,nat);
 454   snew(new_gb_radius,nat);
 455   snew(new_S_hct,nat);
 456   snew(new_atomname,nat);
 457
 458   /* We now have a list of unique atomtypes in typelist */
 459
 460   if (debug)
 461     pr_ivec(debug,0,"typelist",typelist,nat,TRUE);
 462
 463   /* Renumber nlist */
 464   nbsnew = NULL;
 465   snew(nbsnew,plist[ftype].nr);
 466
 467   nrfp  = NRFP(ftype);
 468
 469   for(i=k=0; (i<nat); i++)
 470   {
 471     mi=typelist[i];
 472     for(j=0; (j<nat); j++,k++)
 473     {
 474       mj=typelist[j];
 475       for(l=0; (l<nrfp); l++)
 476       {
 477         nbsnew[k].c[l]=plist[ftype].param[ntype*mi+mj].c[l];
 478       }
 479     }
 480     new_radius[i]     = get_atomtype_radius(mi,ga);
 481     new_vol[i]        = get_atomtype_vol(mi,ga);
 482     new_surftens[i]   = get_atomtype_surftens(mi,ga);
 483     new_atomnumber[i] = get_atomtype_atomnumber(mi,ga);
 484     new_gb_radius[i]  = get_atomtype_gb_radius(mi,ga);
 485     new_S_hct[i]      = get_atomtype_S_hct(mi,ga);
 486     new_atomname[i]   = ga->atomname[mi];
 487   }
 488
 489   for(i=0; (i<nat*nat); i++) {
 490     for(l=0; (l<nrfp); l++)
 491       plist[ftype].param[i].c[l]=nbsnew[i].c[l];
 492   }
 493   plist[ftype].nr=i;
 494   mtop->ffparams.atnr = nat;
 495
 496   sfree(ga->radius);
 497   sfree(ga->vol);
 498   sfree(ga->surftens);
 499   sfree(ga->atomnumber);
 500   sfree(ga->gb_radius);
 501   sfree(ga->S_hct);
 502   /* Dangling atomname pointers ? */
 503   sfree(ga->atomname);
 504
 505   ga->radius     = new_radius;
 506   ga->vol        = new_vol;
 507   ga->surftens   = new_surftens;
 508   ga->atomnumber = new_atomnumber;
 509   ga->gb_radius  = new_gb_radius;
 510   ga->S_hct      = new_S_hct;
 511   ga->atomname   = new_atomname;
 512
 513   ga->nr=nat;
 514
 515   sfree(nbsnew);
 516   sfree(typelist);
 517 }
 518
 519 void copy_atomtype_atomtypes(gpp_atomtype_t ga,t_atomtypes *atomtypes)
 520 {
 521   int i,ntype;
 522
 523   /* Copy the atomtype data to the topology atomtype list */
 524   ntype = get_atomtype_ntypes(ga);
 525   atomtypes->nr=ntype;
 526   snew(atomtypes->radius,ntype);
 527   snew(atomtypes->vol,ntype);
 528   snew(atomtypes->surftens,ntype);
 529   snew(atomtypes->atomnumber,ntype);
 530   snew(atomtypes->gb_radius,ntype);
 531   snew(atomtypes->S_hct,ntype);
 532
 533   for(i=0; i<ntype; i++) {
 534     atomtypes->radius[i]     = ga->radius[i];
 535     atomtypes->vol[i]        = ga->vol[i];
 536     atomtypes->surftens[i]   = ga->surftens[i];
 537     atomtypes->atomnumber[i] = ga->atomnumber[i];
 538     atomtypes->gb_radius[i]  = ga->gb_radius[i];
 539     atomtypes->S_hct[i]      = ga->S_hct[i];
 540   }
 541 }
 542