src/kernel/gpp_atomtype.c

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, by the GROMACS development team, led by
   8  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   9  * others, as listed in the AUTHORS file in the top-level source
  10  * directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include <math.h>
  43 #include "smalloc.h"
  44 #include "sysstuff.h"
  45 #include "macros.h"
  46 #include "string2.h"
  47 #include "topdirs.h"
  48 #include "toputil.h"
  49 #include "topdirs.h"
  50 #include "toputil.h"
  51 #include "symtab.h"
  52 #include "gmx_fatal.h"
  53 #include "txtdump.h"
  54 #include "gpp_atomtype.h"
  55
  56 typedef struct gpp_atomtype {
  57   int           nr;             /* The number of atomtypes              */
  58   t_atom        *atom;          /* Array of atoms                       */
  59   char          ***atomname;    /* Names of the atomtypes               */
  60   t_param       *nb;            /* Nonbonded force default params       */
  61   int           *bondatomtype;  /* The bond_atomtype for each atomtype  */
  62   real          *radius;        /* Radius for GBSA stuff                */
  63   real          *vol;           /* Effective volume for GBSA            */
  64   real          *surftens;      /* Surface tension with water, for GBSA */
  65   real          *gb_radius;     /* Radius for Still model               */
  66   real          *S_hct;         /* Overlap factor for HCT model         */
  67   int           *atomnumber;    /* Atomic number, used for QM/MM        */
  68 } t_gpp_atomtype;
  69
  70 int get_atomtype_type(const char *str,gpp_atomtype_t ga)
  71 {
  72   int i;
  73
  74   /* Atom types are always case sensitive */
  75   for (i=0; (i<ga->nr); i++)
  76     if (strcmp(str,*(ga->atomname[i])) == 0)
  77       return i;
  78
  79   return NOTSET;
  80 }
  81
  82 int get_atomtype_ntypes(gpp_atomtype_t ga)
  83 {
  84   return ga->nr;
  85 }
  86
  87 char *get_atomtype_name(int nt, gpp_atomtype_t ga)
  88 {
  89   if ((nt < 0) || (nt >= ga->nr))
  90     return NULL;
  91
  92   return *(ga->atomname[nt]);
  93 }
  94
  95 real get_atomtype_massA(int nt,gpp_atomtype_t ga)
  96 {
  97   if ((nt < 0) || (nt >= ga->nr))
  98     return NOTSET;
  99
 100   return ga->atom[nt].m;
 101 }
 102
 103 real get_atomtype_massB(int nt,gpp_atomtype_t ga)
 104 {
 105   if ((nt < 0) || (nt >= ga->nr))
 106     return NOTSET;
 107
 108   return ga->atom[nt].mB;
 109 }
 110
 111 real get_atomtype_qA(int nt,gpp_atomtype_t ga)
 112 {
 113   if ((nt < 0) || (nt >= ga->nr))
 114     return NOTSET;
 115
 116   return ga->atom[nt].q;
 117 }
 118
 119 real get_atomtype_qB(int nt,gpp_atomtype_t ga)
 120 {
 121   if ((nt < 0) || (nt >= ga->nr))
 122     return NOTSET;
 123
 124   return ga->atom[nt].qB;
 125 }
 126
 127 int get_atomtype_ptype(int nt,gpp_atomtype_t ga)
 128 {
 129   if ((nt < 0) || (nt >= ga->nr))
 130     return NOTSET;
 131
 132   return ga->atom[nt].ptype;
 133 }
 134
 135 int get_atomtype_batype(int nt,gpp_atomtype_t ga)
 136 {
 137   if ((nt < 0) || (nt >= ga->nr))
 138     return NOTSET;
 139
 140   return ga->bondatomtype[nt];
 141 }
 142
 143 int get_atomtype_atomnumber(int nt,gpp_atomtype_t ga)
 144 {
 145   if ((nt < 0) || (nt >= ga->nr))
 146     return NOTSET;
 147
 148   return ga->atomnumber[nt];
 149 }
 150
 151 real get_atomtype_radius(int nt,gpp_atomtype_t ga)
 152 {
 153   if ((nt < 0) || (nt >= ga->nr))
 154     return NOTSET;
 155
 156   return ga->radius[nt];
 157 }
 158
 159 real get_atomtype_vol(int nt,gpp_atomtype_t ga)
 160 {
 161   if ((nt < 0) || (nt >= ga->nr))
 162     return NOTSET;
 163
 164   return ga->vol[nt];
 165 }
 166
 167 real get_atomtype_surftens(int nt,gpp_atomtype_t ga)
 168 {
 169   if ((nt < 0) || (nt >= ga->nr))
 170     return NOTSET;
 171
 172   return ga->surftens[nt];
 173 }
 174
 175 real get_atomtype_gb_radius(int nt,gpp_atomtype_t ga)
 176 {
 177   if ((nt < 0) || (nt >= ga->nr))
 178     return NOTSET;
 179
 180   return ga->gb_radius[nt];
 181 }
 182
 183 real get_atomtype_S_hct(int nt,gpp_atomtype_t ga)
 184 {
 185   if ((nt < 0) || (nt >= ga->nr))
 186     return NOTSET;
 187
 188   return ga->S_hct[nt];
 189 }
 190
 191 real get_atomtype_nbparam(int nt,int param,gpp_atomtype_t ga)
 192 {
 193   if ((nt < 0) || (nt >= ga->nr))
 194     return NOTSET;
 195   if ((param < 0) || (param >= MAXFORCEPARAM))
 196     return NOTSET;
 197   return ga->nb[nt].c[param];
 198 }
 199
 200 gpp_atomtype_t init_atomtype(void)
 201 {
 202   gpp_atomtype_t ga;
 203
 204   snew(ga,1);
 205
 206   ga->nr           = 0;
 207   ga->atom         = NULL;
 208   ga->atomname     = NULL;
 209   ga->nb           = NULL;
 210   ga->bondatomtype = NULL;
 211   ga->radius       = NULL;
 212   ga->vol          = NULL;
 213   ga->surftens     = NULL;
 214   ga->atomnumber   = NULL;
 215   ga->gb_radius    = NULL;
 216   ga->S_hct        = NULL;
 217
 218   return ga;
 219 }
 220
 221 int
 222 set_atomtype_gbparam(gpp_atomtype_t ga, int i,
 223                      real radius,real vol,real surftens,
 224                      real gb_radius, real S_hct)
 225 {
 226   if ( (i < 0) || (i>= ga->nr))
 227     return NOTSET;
 228
 229   ga->radius[i]    = radius;
 230   ga->vol[i]       = vol;
 231   ga->surftens[i]  = surftens;
 232   ga->gb_radius[i] = gb_radius;
 233   ga->S_hct[i]     = S_hct;
 234
 235   return i;
 236 }
 237
 238
 239 int set_atomtype(int nt,gpp_atomtype_t ga,t_symtab *tab,
 240                  t_atom *a,const char *name,t_param *nb,
 241                  int bondatomtype,
 242                  real radius,real vol,real surftens,int atomnumber,
 243                  real gb_radius, real S_hct)
 244 {
 245   if ((nt < 0) || (nt >= ga->nr))
 246     return NOTSET;
 247
 248   ga->atom[nt] = *a;
 249   ga->atomname[nt] = put_symtab(tab,name);
 250   ga->nb[nt] = *nb;
 251   ga->bondatomtype[nt] = bondatomtype;
 252   ga->radius[nt] = radius;
 253   ga->vol[nt] = vol;
 254   ga->surftens[nt] = surftens;
 255   ga->atomnumber[nt] = atomnumber;
 256   ga->gb_radius[nt] = gb_radius;
 257   ga->S_hct[nt]     = S_hct;
 258
 259   return nt;
 260 }
 261
 262 int add_atomtype(gpp_atomtype_t ga,t_symtab *tab,
 263                  t_atom *a,const char *name,t_param *nb,
 264                  int bondatomtype,
 265                  real radius,real vol,real surftens,real atomnumber,
 266                  real gb_radius, real S_hct)
 267 {
 268   int i;
 269
 270   for(i=0; (i<ga->nr); i++) {
 271     if (strcmp(*ga->atomname[i],name) == 0) {
 272       if (NULL != debug)
 273         fprintf(debug,"Trying to add atomtype %s again. Skipping it.\n",name);
 274       break;
 275     }
 276   }
 277   if (i == ga->nr) {
 278     ga->nr++;
 279     srenew(ga->atom,ga->nr);
 280     srenew(ga->atomname,ga->nr);
 281     srenew(ga->nb,ga->nr);
 282     srenew(ga->bondatomtype,ga->nr);
 283     srenew(ga->radius,ga->nr);
 284     srenew(ga->vol,ga->nr);
 285     srenew(ga->surftens,ga->nr);
 286     srenew(ga->atomnumber,ga->nr);
 287     srenew(ga->gb_radius,ga->nr);
 288     srenew(ga->S_hct,ga->nr);
 289
 290     return set_atomtype(ga->nr-1,ga,tab,a,name,nb,bondatomtype,radius,
 291                         vol,surftens,atomnumber,gb_radius,S_hct);
 292   }
 293   else
 294     return i;
 295 }
 296
 297 void print_at (FILE * out, gpp_atomtype_t ga)
 298 {
 299   int i;
 300   t_atom     *atom = ga->atom;
 301   t_param    *nb   = ga->nb;
 302
 303   fprintf (out,"[ %s ]\n",dir2str(d_atomtypes));
 304   fprintf (out,"; %6s  %8s  %8s  %8s  %12s  %12s\n",
 305            "type","mass","charge","particle","c6","c12");
 306   for (i=0; (i<ga->nr); i++)
 307     fprintf(out,"%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
 308             *(ga->atomname[i]),atom[i].m,atom[i].q,"A",
 309             nb[i].C0,nb[i].C1);
 310
 311   fprintf (out,"\n");
 312 }
 313
 314 void done_atomtype(gpp_atomtype_t ga)
 315 {
 316   sfree(ga->atom);
 317   sfree(ga->atomname);
 318   sfree(ga->nb);
 319   sfree(ga->bondatomtype);
 320   sfree(ga->radius);
 321   sfree(ga->vol);
 322   sfree(ga->gb_radius);
 323   sfree(ga->S_hct);
 324   sfree(ga->surftens);
 325   sfree(ga->atomnumber);
 326   ga->nr = 0;
 327   sfree(ga);
 328   ga = NULL;
 329 }
 330
 331 static int search_atomtypes(gpp_atomtype_t ga,int *n,int typelist[],
 332                             int thistype,
 333                             t_param param[],int ftype)
 334 {
 335   int i,nn,nrfp,j,k,ntype,tli;
 336   gmx_bool bFound=FALSE;
 337
 338   nn    = *n;
 339   nrfp  = NRFP(ftype);
 340   ntype = get_atomtype_ntypes(ga);
 341
 342   for(i=0; (i<nn); i++)
 343   {
 344     if (typelist[i] == thistype)
 345     {
 346       /* This type number has already been added */
 347       break;
 348     }
 349
 350     /* Otherwise, check if the parameters are identical to any previously added type */
 351
 352     bFound=TRUE;
 353     for(j=0;j<ntype && bFound;j++)
 354     {
 355       /* Check nonbonded parameters */
 356       for(k=0;k<nrfp && bFound;k++)
 357       {
 358         bFound=(param[ntype*typelist[i]+j].c[k]==param[ntype*thistype+j].c[k]);
 359       }
 360
 361       /* Check radius, volume, surftens */
 362       tli = typelist[i];
 363       bFound = bFound &&
 364         (get_atomtype_radius(tli,ga) == get_atomtype_radius(thistype,ga)) &&
 365         (get_atomtype_vol(tli,ga) == get_atomtype_vol(thistype,ga)) &&
 366         (get_atomtype_surftens(tli,ga) == get_atomtype_surftens(thistype,ga)) &&
 367         (get_atomtype_atomnumber(tli,ga) == get_atomtype_atomnumber(thistype,ga)) &&
 368         (get_atomtype_gb_radius(tli,ga) == get_atomtype_gb_radius(thistype,ga)) &&
 369         (get_atomtype_S_hct(tli,ga) == get_atomtype_S_hct(thistype,ga));
 370     }
 371     if (bFound)
 372     {
 373       break;
 374     }
 375   }
 376
 377   if (i == nn) {
 378     if (debug)
 379       fprintf(debug,"Renumbering atomtype %d to %d\n",thistype,nn);
 380     if (nn == ntype)
 381       gmx_fatal(FARGS,"Atomtype horror n = %d, %s, %d",nn,__FILE__,__LINE__);
 382     typelist[nn]=thistype;
 383     nn++;
 384   }
 385   *n = nn;
 386
 387   return i;
 388 }
 389
 390 void renum_atype(t_params plist[],gmx_mtop_t *mtop,
 391                 int *wall_atomtype,
 392                 gpp_atomtype_t ga,gmx_bool bVerbose)
 393 {
 394   int      i,j,k,l,molt,mi,mj,nat,nrfp,ftype,ntype;
 395   t_atoms  *atoms;
 396   t_param  *nbsnew;
 397   int      *typelist;
 398   real     *new_radius;
 399   real     *new_vol;
 400   real     *new_surftens;
 401   real     *new_gb_radius;
 402   real     *new_S_hct;
 403   int      *new_atomnumber;
 404   char     ***new_atomname;
 405
 406   ntype = get_atomtype_ntypes(ga);
 407   snew(typelist,ntype);
 408
 409   if (bVerbose)
 410     fprintf(stderr,"renumbering atomtypes...\n");
 411
 412   /* Since the bonded interactions have been assigned now,
 413    * we want to reduce the number of atom types by merging
 414    * ones with identical nonbonded interactions, in addition
 415    * to removing unused ones.
 416    *
 417    * With Generalized-Born electrostatics, or implicit solvent
 418    * we also check that the atomtype radius, effective_volume
 419    * and surface tension match.
 420    *
 421    * With QM/MM we also check that the atom numbers match
 422    */
 423
 424   /* Get nonbonded interaction type */
 425   if (plist[F_LJ].nr > 0)
 426     ftype=F_LJ;
 427   else
 428     ftype=F_BHAM;
 429
 430   /* Renumber atomtypes by first making a list of which ones are actually used.
 431    * We provide the list of nonbonded parameters so search_atomtypes
 432    * can determine if two types should be merged.
 433    */
 434   nat=0;
 435   for(molt=0; molt<mtop->nmoltype; molt++) {
 436     atoms = &mtop->moltype[molt].atoms;
 437     for(i=0; (i<atoms->nr); i++) {
 438       atoms->atom[i].type =
 439         search_atomtypes(ga,&nat,typelist,atoms->atom[i].type,
 440                          plist[ftype].param,ftype);
 441       atoms->atom[i].typeB=
 442         search_atomtypes(ga,&nat,typelist,atoms->atom[i].typeB,
 443                          plist[ftype].param,ftype);
 444     }
 445   }
 446
 447   for(i=0; i<2; i++) {
 448     if (wall_atomtype[i] >= 0)
 449       wall_atomtype[i] = search_atomtypes(ga,&nat,typelist,wall_atomtype[i],
 450                                           plist[ftype].param,ftype);
 451   }
 452
 453   snew(new_radius,nat);
 454   snew(new_vol,nat);
 455   snew(new_surftens,nat);
 456   snew(new_atomnumber,nat);
 457   snew(new_gb_radius,nat);
 458   snew(new_S_hct,nat);
 459   snew(new_atomname,nat);
 460
 461   /* We now have a list of unique atomtypes in typelist */
 462
 463   if (debug)
 464     pr_ivec(debug,0,"typelist",typelist,nat,TRUE);
 465
 466   /* Renumber nlist */
 467   nbsnew = NULL;
 468   snew(nbsnew,plist[ftype].nr);
 469
 470   nrfp  = NRFP(ftype);
 471
 472   for(i=k=0; (i<nat); i++)
 473   {
 474     mi=typelist[i];
 475     for(j=0; (j<nat); j++,k++)
 476     {
 477       mj=typelist[j];
 478       for(l=0; (l<nrfp); l++)
 479       {
 480         nbsnew[k].c[l]=plist[ftype].param[ntype*mi+mj].c[l];
 481       }
 482     }
 483     new_radius[i]     = get_atomtype_radius(mi,ga);
 484     new_vol[i]        = get_atomtype_vol(mi,ga);
 485     new_surftens[i]   = get_atomtype_surftens(mi,ga);
 486     new_atomnumber[i] = get_atomtype_atomnumber(mi,ga);
 487     new_gb_radius[i]  = get_atomtype_gb_radius(mi,ga);
 488     new_S_hct[i]      = get_atomtype_S_hct(mi,ga);
 489     new_atomname[i]   = ga->atomname[mi];
 490   }
 491
 492   for(i=0; (i<nat*nat); i++) {
 493     for(l=0; (l<nrfp); l++)
 494       plist[ftype].param[i].c[l]=nbsnew[i].c[l];
 495   }
 496   plist[ftype].nr=i;
 497   mtop->ffparams.atnr = nat;
 498
 499   sfree(ga->radius);
 500   sfree(ga->vol);
 501   sfree(ga->surftens);
 502   sfree(ga->atomnumber);
 503   sfree(ga->gb_radius);
 504   sfree(ga->S_hct);
 505   /* Dangling atomname pointers ? */
 506   sfree(ga->atomname);
 507
 508   ga->radius     = new_radius;
 509   ga->vol        = new_vol;
 510   ga->surftens   = new_surftens;
 511   ga->atomnumber = new_atomnumber;
 512   ga->gb_radius  = new_gb_radius;
 513   ga->S_hct      = new_S_hct;
 514   ga->atomname   = new_atomname;
 515
 516   ga->nr=nat;
 517
 518   sfree(nbsnew);
 519   sfree(typelist);
 520 }
 521
 522 void copy_atomtype_atomtypes(gpp_atomtype_t ga,t_atomtypes *atomtypes)
 523 {
 524   int i,ntype;
 525
 526   /* Copy the atomtype data to the topology atomtype list */
 527   ntype = get_atomtype_ntypes(ga);
 528   atomtypes->nr=ntype;
 529   snew(atomtypes->radius,ntype);
 530   snew(atomtypes->vol,ntype);
 531   snew(atomtypes->surftens,ntype);
 532   snew(atomtypes->atomnumber,ntype);
 533   snew(atomtypes->gb_radius,ntype);
 534   snew(atomtypes->S_hct,ntype);
 535
 536   for(i=0; i<ntype; i++) {
 537     atomtypes->radius[i]     = ga->radius[i];
 538     atomtypes->vol[i]        = ga->vol[i];
 539     atomtypes->surftens[i]   = ga->surftens[i];
 540     atomtypes->atomnumber[i] = ga->atomnumber[i];
 541     atomtypes->gb_radius[i]  = ga->gb_radius[i];
 542     atomtypes->S_hct[i]      = ga->S_hct[i];
 543   }
 544 }
 545