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