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,2019, 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 <climits>
  42 #include <cmath>
  43 #include <cstring>
  44
  45 #include "gromacs/gmxpreprocess/grompp_impl.h"
  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 struct gpp_atomtype
  57 {
  58     int              nr;           /* The number of atomtypes           */
  59     t_atom          *atom;         /* Array of atoms                    */
  60     char          ***atomname;     /* Names of the atomtypes            */
  61     t_param         *nb;           /* Nonbonded force default params    */
  62     int             *bondatomtype; /* The bond_atomtype for each atomtype  */
  63     int             *atomnumber;   /* Atomic number, used for QM/MM        */
  64 };
  65
  66 int get_atomtype_type(const char *str, gpp_atomtype *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 *ga)
  83 {
  84     return ga->nr;
  85 }
  86
  87 char *get_atomtype_name(int nt, gpp_atomtype *ga)
  88 {
  89     if ((nt < 0) || (nt >= ga->nr))
  90     {
  91         return nullptr;
  92     }
  93
  94     return *(ga->atomname[nt]);
  95 }
  96
  97 real get_atomtype_massA(int nt, gpp_atomtype *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 *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 *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 *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 *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, const gpp_atomtype* 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 *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_nbparam(int nt, int param, gpp_atomtype *ga)
 168 {
 169     if ((nt < 0) || (nt >= ga->nr))
 170     {
 171         return NOTSET;
 172     }
 173     if ((param < 0) || (param >= MAXFORCEPARAM))
 174     {
 175         return NOTSET;
 176     }
 177     return ga->nb[nt].c[param];
 178 }
 179
 180 gpp_atomtype *init_atomtype()
 181 {
 182     gpp_atomtype *ga;
 183
 184     snew(ga, 1);
 185
 186     ga->nr           = 0;
 187     ga->atom         = nullptr;
 188     ga->atomname     = nullptr;
 189     ga->nb           = nullptr;
 190     ga->bondatomtype = nullptr;
 191     ga->atomnumber   = nullptr;
 192
 193     return ga;
 194 }
 195
 196 int set_atomtype(int nt, gpp_atomtype *ga, t_symtab *tab,
 197                  t_atom *a, const char *name, t_param *nb,
 198                  int bondatomtype, int atomnumber)
 199 {
 200     if ((nt < 0) || (nt >= ga->nr))
 201     {
 202         return NOTSET;
 203     }
 204
 205     ga->atom[nt]         = *a;
 206     ga->atomname[nt]     = put_symtab(tab, name);
 207     ga->nb[nt]           = *nb;
 208     ga->bondatomtype[nt] = bondatomtype;
 209     ga->atomnumber[nt]   = atomnumber;
 210
 211     return nt;
 212 }
 213
 214 int add_atomtype(gpp_atomtype *ga, t_symtab *tab,
 215                  t_atom *a, const char *name, t_param *nb,
 216                  int bondatomtype, int atomnumber)
 217 {
 218     int i;
 219
 220     for (i = 0; (i < ga->nr); i++)
 221     {
 222         if (strcmp(*ga->atomname[i], name) == 0)
 223         {
 224             break;
 225         }
 226     }
 227     if (i == ga->nr)
 228     {
 229         ga->nr++;
 230         srenew(ga->atom, ga->nr);
 231         srenew(ga->atomname, ga->nr);
 232         srenew(ga->nb, ga->nr);
 233         srenew(ga->bondatomtype, ga->nr);
 234         srenew(ga->atomnumber, ga->nr);
 235
 236         return set_atomtype(ga->nr-1, ga, tab, a, name, nb, bondatomtype, atomnumber);
 237     }
 238     else
 239     {
 240         return i;
 241     }
 242 }
 243
 244 void print_at (FILE * out, gpp_atomtype *ga)
 245 {
 246     int         i;
 247     t_atom     *atom = ga->atom;
 248     t_param    *nb   = ga->nb;
 249
 250     fprintf (out, "[ %s ]\n", dir2str(Directive::d_atomtypes));
 251     fprintf (out, "; %6s  %8s  %8s  %8s  %12s  %12s\n",
 252              "type", "mass", "charge", "particle", "c6", "c12");
 253     for (i = 0; (i < ga->nr); i++)
 254     {
 255         fprintf(out, "%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
 256                 *(ga->atomname[i]), atom[i].m, atom[i].q, "A",
 257                 nb[i].c0(), nb[i].c1());
 258     }
 259
 260     fprintf (out, "\n");
 261 }
 262
 263 void done_atomtype(gpp_atomtype *ga)
 264 {
 265     sfree(ga->atom);
 266     sfree(ga->atomname);
 267     sfree(ga->nb);
 268     sfree(ga->bondatomtype);
 269     sfree(ga->atomnumber);
 270     ga->nr = 0;
 271     sfree(ga);
 272 }
 273
 274 static int search_atomtypes(gpp_atomtype *ga, int *n, int typelist[],
 275                             int thistype,
 276                             t_param param[], int ftype)
 277 {
 278     int      i, nn, nrfp, j, k, ntype, tli;
 279     bool     bFound = FALSE;
 280
 281     nn    = *n;
 282     nrfp  = NRFP(ftype);
 283     ntype = get_atomtype_ntypes(ga);
 284
 285     for (i = 0; (i < nn); i++)
 286     {
 287         if (typelist[i] == thistype)
 288         {
 289             /* This type number has already been added */
 290             break;
 291         }
 292
 293         /* Otherwise, check if the parameters are identical to any previously added type */
 294
 295         bFound = TRUE;
 296         for (j = 0; j < ntype && bFound; j++)
 297         {
 298             /* Check nonbonded parameters */
 299             for (k = 0; k < nrfp && bFound; k++)
 300             {
 301                 bFound = (param[ntype*typelist[i]+j].c[k] == param[ntype*thistype+j].c[k]);
 302             }
 303
 304             /* Check atomnumber */
 305             tli    = typelist[i];
 306             bFound = bFound &&
 307                 (get_atomtype_atomnumber(tli, ga) == get_atomtype_atomnumber(thistype, ga));
 308         }
 309         if (bFound)
 310         {
 311             break;
 312         }
 313     }
 314
 315     if (i == nn)
 316     {
 317         if (nn == ntype)
 318         {
 319             gmx_fatal(FARGS, "Atomtype horror n = %d, %s, %d", nn, __FILE__, __LINE__);
 320         }
 321         typelist[nn] = thistype;
 322         nn++;
 323     }
 324     *n = nn;
 325
 326     return i;
 327 }
 328
 329 void renum_atype(gmx::ArrayRef<InteractionTypeParameters> plist, gmx_mtop_t *mtop,
 330                  int *wall_atomtype,
 331                  gpp_atomtype *ga, bool bVerbose)
 332 {
 333     int         i, j, k, l, mi, mj, nat, nrfp, ftype, ntype;
 334     t_atoms    *atoms;
 335     t_param    *nbsnew;
 336     int        *typelist;
 337     int        *new_atomnumber;
 338     char     ***new_atomname;
 339
 340     ntype = get_atomtype_ntypes(ga);
 341     snew(typelist, ntype);
 342
 343     if (bVerbose)
 344     {
 345         fprintf(stderr, "renumbering atomtypes...\n");
 346     }
 347
 348     /* Since the bonded interactions have been assigned now,
 349      * we want to reduce the number of atom types by merging
 350      * ones with identical nonbonded interactions, in addition
 351      * to removing unused ones.
 352      *
 353      * With QM/MM we also check that the atom numbers match
 354      */
 355
 356     /* Get nonbonded interaction type */
 357     if (plist[F_LJ].nr > 0)
 358     {
 359         ftype = F_LJ;
 360     }
 361     else
 362     {
 363         ftype = F_BHAM;
 364     }
 365
 366     /* Renumber atomtypes by first making a list of which ones are actually used.
 367      * We provide the list of nonbonded parameters so search_atomtypes
 368      * can determine if two types should be merged.
 369      */
 370     nat = 0;
 371     for (gmx_moltype_t &moltype : mtop->moltype)
 372     {
 373         atoms = &moltype.atoms;
 374         for (i = 0; (i < atoms->nr); i++)
 375         {
 376             atoms->atom[i].type =
 377                 search_atomtypes(ga, &nat, typelist, atoms->atom[i].type,
 378                                  plist[ftype].param, ftype);
 379             atoms->atom[i].typeB =
 380                 search_atomtypes(ga, &nat, typelist, atoms->atom[i].typeB,
 381                                  plist[ftype].param, ftype);
 382         }
 383     }
 384
 385     for (i = 0; i < 2; i++)
 386     {
 387         if (wall_atomtype[i] >= 0)
 388         {
 389             wall_atomtype[i] = search_atomtypes(ga, &nat, typelist, wall_atomtype[i],
 390                                                 plist[ftype].param, ftype);
 391         }
 392     }
 393
 394     snew(new_atomnumber, nat);
 395     snew(new_atomname, nat);
 396     /* We now have a list of unique atomtypes in typelist */
 397
 398     /* Renumber nlist */
 399     nbsnew = nullptr;
 400     snew(nbsnew, plist[ftype].nr);
 401
 402     nrfp  = NRFP(ftype);
 403
 404     for (i = k = 0; (i < nat); i++)
 405     {
 406         mi = typelist[i];
 407         for (j = 0; (j < nat); j++, k++)
 408         {
 409             mj = typelist[j];
 410             for (l = 0; (l < nrfp); l++)
 411             {
 412                 nbsnew[k].c[l] = plist[ftype].param[ntype*mi+mj].c[l];
 413             }
 414         }
 415         new_atomnumber[i] = get_atomtype_atomnumber(mi, ga);
 416         new_atomname[i]   = ga->atomname[mi];
 417     }
 418
 419     for (i = 0; (i < nat*nat); i++)
 420     {
 421         for (l = 0; (l < nrfp); l++)
 422         {
 423             plist[ftype].param[i].c[l] = nbsnew[i].c[l];
 424         }
 425     }
 426     plist[ftype].nr     = i;
 427     mtop->ffparams.atnr = nat;
 428
 429     sfree(ga->atomnumber);
 430     /* Dangling atomname pointers ? */
 431     sfree(ga->atomname);
 432
 433     ga->atomnumber = new_atomnumber;
 434     ga->atomname   = new_atomname;
 435
 436     ga->nr = nat;
 437
 438     sfree(nbsnew);
 439     sfree(typelist);
 440 }
 441
 442 void copy_atomtype_atomtypes(gpp_atomtype *ga, t_atomtypes *atomtypes)
 443 {
 444     int i, ntype;
 445
 446     /* Copy the atomtype data to the topology atomtype list */
 447     ntype         = get_atomtype_ntypes(ga);
 448     atomtypes->nr = ntype;
 449     snew(atomtypes->atomnumber, ntype);
 450
 451     for (i = 0; i < ntype; i++)
 452     {
 453         atomtypes->atomnumber[i] = ga->atomnumber[i];
 454     }
 455 }