src/gromacs/gmxpreprocess/gpp_atomtype.cpp

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