src/gromacs/gmxpreprocess/gpp_atomtype.c

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