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