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