src/gromacs/gmxpreprocess/toputil.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) 2012,2014,2015, 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 "toputil.h"
  40
  41 #include <string.h>
  42
  43 #include <cmath>
  44
  45 #include <algorithm>
  46
  47 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  48 #include "gromacs/gmxpreprocess/topdirs.h"
  49 #include "gromacs/legacyheaders/macros.h"
  50 #include "gromacs/topology/block.h"
  51 #include "gromacs/topology/symtab.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/gmxassert.h"
  54 #include "gromacs/utility/smalloc.h"
  55
  56 /* UTILITIES */
  57
  58 void set_p_string(t_param *p, const char *s)
  59 {
  60     if (s)
  61     {
  62         if (strlen(s) < sizeof(p->s)-1)
  63         {
  64             strncpy(p->s, s, sizeof(p->s));
  65         }
  66         else
  67         {
  68             gmx_fatal(FARGS, "Increase MAXSLEN in the grompp code to at least %d,"
  69                       " or shorten your definition of bonds like %s to at most %d",
  70                       strlen(s)+1, s, MAXSLEN-1);
  71         }
  72     }
  73     else
  74     {
  75         strcpy(p->s, "");
  76     }
  77 }
  78
  79 void pr_alloc (int extra, t_params *pr)
  80 {
  81     int i, j;
  82
  83     /* get new space for arrays */
  84     if (extra < 0)
  85     {
  86         gmx_fatal(FARGS, "Trying to make array smaller.\n");
  87     }
  88     if (extra == 0)
  89     {
  90         return;
  91     }
  92     GMX_ASSERT(pr->nr != 0 || pr->param == NULL, "Invalid t_params object");
  93     if (pr->nr+extra > pr->maxnr)
  94     {
  95         pr->maxnr = std::max(static_cast<int>(1.2*pr->maxnr), pr->maxnr + extra);
  96         srenew(pr->param, pr->maxnr);
  97         for (i = pr->nr; (i < pr->maxnr); i++)
  98         {
  99             for (j = 0; (j < MAXATOMLIST); j++)
 100             {
 101                 pr->param[i].a[j] = 0;
 102             }
 103             for (j = 0; (j < MAXFORCEPARAM); j++)
 104             {
 105                 pr->param[i].c[j] = 0;
 106             }
 107             set_p_string(&(pr->param[i]), "");
 108         }
 109     }
 110 }
 111
 112 void init_plist(t_params plist[])
 113 {
 114     int i;
 115
 116     for (i = 0; (i < F_NRE); i++)
 117     {
 118         plist[i].nr    = 0;
 119         plist[i].maxnr = 0;
 120         plist[i].param = NULL;
 121
 122         /* CMAP */
 123         plist[i].ncmap        = 0;
 124         plist[i].cmap         = NULL;
 125         plist[i].grid_spacing = 0;
 126         plist[i].nc           = 0;
 127         plist[i].nct          = 0;
 128         plist[i].cmap_types   = NULL;
 129     }
 130 }
 131
 132 void cp_param(t_param *dest, t_param *src)
 133 {
 134     int j;
 135
 136     for (j = 0; (j < MAXATOMLIST); j++)
 137     {
 138         dest->a[j] = src->a[j];
 139     }
 140     for (j = 0; (j < MAXFORCEPARAM); j++)
 141     {
 142         dest->c[j] = src->c[j];
 143     }
 144     strncpy(dest->s, src->s, sizeof(dest->s));
 145 }
 146
 147 void add_param_to_list(t_params *list, t_param *b)
 148 {
 149     int j;
 150
 151     /* allocate one position extra */
 152     pr_alloc (1, list);
 153
 154     /* fill the arrays */
 155     for (j = 0; (j < MAXFORCEPARAM); j++)
 156     {
 157         list->param[list->nr].c[j]   = b->c[j];
 158     }
 159     for (j = 0; (j < MAXATOMLIST); j++)
 160     {
 161         list->param[list->nr].a[j]   = b->a[j];
 162     }
 163     memset(list->param[list->nr].s, 0, sizeof(list->param[list->nr].s));
 164
 165     list->nr++;
 166 }
 167
 168
 169 void init_molinfo(t_molinfo *mol)
 170 {
 171     mol->nrexcl     = 0;
 172     mol->excl_set   = FALSE;
 173     mol->bProcessed = FALSE;
 174     init_plist(mol->plist);
 175     init_block(&mol->cgs);
 176     init_block(&mol->mols);
 177     init_blocka(&mol->excls);
 178     init_atom(&mol->atoms);
 179 }
 180
 181 /* FREEING MEMORY */
 182
 183 void done_bt (t_params *pl)
 184 {
 185     sfree(pl->param);
 186 }
 187
 188 void done_mi(t_molinfo *mi)
 189 {
 190     int i;
 191
 192     done_atom (&(mi->atoms));
 193     done_block(&(mi->cgs));
 194     done_block(&(mi->mols));
 195     for (i = 0; (i < F_NRE); i++)
 196     {
 197         done_bt(&(mi->plist[i]));
 198     }
 199 }
 200
 201 /* PRINTING STRUCTURES */
 202
 203 void print_bt(FILE *out, directive d, gpp_atomtype_t at,
 204               int ftype, int fsubtype, t_params plist[],
 205               gmx_bool bFullDih)
 206 {
 207     /* This dihp is a DIRTY patch because the dih-types do not use
 208      * all four atoms to determine the type.
 209      */
 210     const int    dihp[2][2] = { { 1, 2 }, { 0, 3 } };
 211     t_params    *bt;
 212     int          i, j, f, nral, nrfp;
 213     gmx_bool     bDih = FALSE, bSwapParity;
 214
 215     bt = &(plist[ftype]);
 216
 217     if (!bt->nr)
 218     {
 219         return;
 220     }
 221
 222     f = 0;
 223     switch (ftype)
 224     {
 225         case F_G96ANGLES:
 226             f = 1;
 227             break;
 228         case F_G96BONDS:
 229             f = 1;
 230             break;
 231         case F_MORSE:
 232             f = 2;
 233             break;
 234         case F_CUBICBONDS:
 235             f = 3;
 236             break;
 237         case F_CONNBONDS:
 238             f = 4;
 239             break;
 240         case F_HARMONIC:
 241             f = 5;
 242             break;
 243         case F_CROSS_BOND_ANGLES:
 244             f = 2;
 245             break;
 246         case F_CROSS_BOND_BONDS:
 247             f = 3;
 248             break;
 249         case F_UREY_BRADLEY:
 250             f = 4;
 251             break;
 252         case F_PDIHS:
 253         case F_RBDIHS:
 254         case F_FOURDIHS:
 255             bDih = TRUE;
 256             break;
 257         case F_IDIHS:
 258             f    = 1;
 259             bDih = TRUE;
 260             break;
 261         case F_CONSTRNC:
 262             f = 1;
 263             break;
 264         case F_VSITE3FD:
 265             f = 1;
 266             break;
 267         case F_VSITE3FAD:
 268             f = 2;
 269             break;
 270         case F_VSITE3OUT:
 271             f = 3;
 272             break;
 273         case F_VSITE4FDN:
 274             f = 1;
 275             break;
 276         case F_CMAP:
 277             f = 1;
 278             break;
 279
 280         default:
 281             bDih = FALSE;
 282     }
 283     if (bFullDih)
 284     {
 285         bDih = FALSE;
 286     }
 287     if (fsubtype)
 288     {
 289         f = fsubtype-1;
 290     }
 291
 292     nral = NRAL(ftype);
 293     nrfp = NRFP(ftype);
 294
 295     /* header */
 296     fprintf(out, "[ %s ]\n", dir2str(d));
 297     fprintf(out, "; ");
 298     if (!bDih)
 299     {
 300         fprintf (out, "%3s  %4s", "ai", "aj");
 301         for (j = 2; (j < nral); j++)
 302         {
 303             fprintf (out, "  %3c%c", 'a', 'i'+j);
 304         }
 305     }
 306     else
 307     {
 308         for (j = 0; (j < 2); j++)
 309         {
 310             fprintf (out, "%3c%c", 'a', 'i'+dihp[f][j]);
 311         }
 312     }
 313
 314     fprintf (out, " funct");
 315     for (j = 0; (j < nrfp); j++)
 316     {
 317         fprintf (out, " %12c%1d", 'c', j);
 318     }
 319     fprintf (out, "\n");
 320
 321     /* print bondtypes */
 322     for (i = 0; (i < bt->nr); i++)
 323     {
 324         bSwapParity = (bt->param[i].C0 == NOTSET) && (bt->param[i].C1 == -1);
 325         if (!bDih)
 326         {
 327             for (j = 0; (j < nral); j++)
 328             {
 329                 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[j], at));
 330             }
 331         }
 332         else
 333         {
 334             for (j = 0; (j < 2); j++)
 335             {
 336                 fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[dihp[f][j]], at));
 337             }
 338         }
 339         fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1);
 340
 341         if (bt->param[i].s[0])
 342         {
 343             fprintf(out, "   %s", bt->param[i].s);
 344         }
 345         else
 346         {
 347             for (j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++)
 348             {
 349                 fprintf (out, "%13.6e ", bt->param[i].c[j]);
 350             }
 351         }
 352
 353         fprintf (out, "\n");
 354     }
 355     fprintf (out, "\n");
 356     fflush (out);
 357 }
 358
 359 void print_blocka(FILE *out, const char *szName,
 360                   const char *szIndex, const char *szA,
 361                   t_blocka *block)
 362 {
 363     int i, j;
 364
 365     fprintf (out, "; %s\n", szName);
 366     fprintf (out, "; %4s    %s\n", szIndex, szA);
 367     for (i = 0; (i < block->nr); i++)
 368     {
 369         for (i = 0; (i < block->nr); i++)
 370         {
 371             fprintf (out, "%6d", i+1);
 372             for (j = block->index[i]; (j < ((int)block->index[i+1])); j++)
 373             {
 374                 fprintf (out, "%5d", block->a[j]+1);
 375             }
 376             fprintf (out, "\n");
 377         }
 378         fprintf (out, "\n");
 379     }
 380 }
 381
 382 void print_excl(FILE *out, int natoms, t_excls excls[])
 383 {
 384     atom_id     i;
 385     gmx_bool    have_excl;
 386     int         j;
 387
 388     have_excl = FALSE;
 389     for (i = 0; i < natoms && !have_excl; i++)
 390     {
 391         have_excl = (excls[i].nr > 0);
 392     }
 393
 394     if (have_excl)
 395     {
 396         fprintf (out, "[ %s ]\n", dir2str(d_exclusions));
 397         fprintf (out, "; %4s    %s\n", "i", "excluded from i");
 398         for (i = 0; i < natoms; i++)
 399         {
 400             if (excls[i].nr > 0)
 401             {
 402                 fprintf (out, "%6d ", i+1);
 403                 for (j = 0; j < excls[i].nr; j++)
 404                 {
 405                     fprintf (out, " %5d", excls[i].e[j]+1);
 406                 }
 407                 fprintf (out, "\n");
 408             }
 409         }
 410         fprintf (out, "\n");
 411         fflush(out);
 412     }
 413 }
 414
 415 static double get_residue_charge(const t_atoms *atoms, int at)
 416 {
 417     int    ri;
 418     double q;
 419
 420     ri = atoms->atom[at].resind;
 421     q  = 0;
 422     while (at < atoms->nr && atoms->atom[at].resind == ri)
 423     {
 424         q += atoms->atom[at].q;
 425         at++;
 426     }
 427
 428     return q;
 429 }
 430
 431 void print_atoms(FILE *out, gpp_atomtype_t atype, t_atoms *at, int *cgnr,
 432                  gmx_bool bRTPresname)
 433 {
 434     int         i, ri;
 435     int         tpA, tpB;
 436     const char *as;
 437     char       *tpnmA, *tpnmB;
 438     double      qres, qtot;
 439
 440     as = dir2str(d_atoms);
 441     fprintf(out, "[ %s ]\n", as);
 442     fprintf(out, "; %4s %10s %6s %7s%6s %6s %10s %10s %6s %10s %10s\n",
 443             "nr", "type", "resnr", "residue", "atom", "cgnr", "charge", "mass", "typeB", "chargeB", "massB");
 444
 445     qtot  = 0;
 446
 447     if (debug)
 448     {
 449         fprintf(debug, "This molecule has %d atoms and %d residues\n",
 450                 at->nr, at->nres);
 451     }
 452
 453     if (at->nres)
 454     {
 455         /* if the information is present... */
 456         for (i = 0; (i < at->nr); i++)
 457         {
 458             ri = at->atom[i].resind;
 459             if ((i == 0 || ri != at->atom[i-1].resind) &&
 460                 at->resinfo[ri].rtp != NULL)
 461             {
 462                 qres = get_residue_charge(at, i);
 463                 fprintf(out, "; residue %3d %-3s rtp %-4s q ",
 464                         at->resinfo[ri].nr,
 465                         *at->resinfo[ri].name,
 466                         *at->resinfo[ri].rtp);
 467                 if (fabs(qres) < 0.001)
 468                 {
 469                     fprintf(out, " %s", "0.0");
 470                 }
 471                 else
 472                 {
 473                     fprintf(out, "%+3.1f", qres);
 474                 }
 475                 fprintf(out, "\n");
 476             }
 477             tpA = at->atom[i].type;
 478             if ((tpnmA = get_atomtype_name(tpA, atype)) == NULL)
 479             {
 480                 gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i);
 481             }
 482
 483             /* This is true by construction, but static analysers don't know */
 484             GMX_ASSERT(!bRTPresname || at->resinfo[at->atom[i].resind].rtp, "-rtpres did not have residue name available");
 485             fprintf(out, "%6d %10s %6d%c %5s %6s %6d %10g %10g",
 486                     i+1, tpnmA,
 487                     at->resinfo[ri].nr,
 488                     at->resinfo[ri].ic,
 489                     bRTPresname ?
 490                     *(at->resinfo[at->atom[i].resind].rtp) :
 491                     *(at->resinfo[at->atom[i].resind].name),
 492                     *(at->atomname[i]), cgnr[i],
 493                     at->atom[i].q, at->atom[i].m);
 494             if (PERTURBED(at->atom[i]))
 495             {
 496                 tpB = at->atom[i].typeB;
 497                 if ((tpnmB = get_atomtype_name(tpB, atype)) == NULL)
 498                 {
 499                     gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i);
 500                 }
 501                 fprintf(out, " %6s %10g %10g",
 502                         tpnmB, at->atom[i].qB, at->atom[i].mB);
 503             }
 504             qtot += (double)at->atom[i].q;
 505             if (fabs(qtot) < 4*GMX_REAL_EPS)
 506             {
 507                 qtot = 0;
 508             }
 509             fprintf(out, "   ; qtot %.4g\n", qtot);
 510         }
 511     }
 512     fprintf(out, "\n");
 513     fflush(out);
 514 }
 515
 516 void print_bondeds(FILE *out, int natoms, directive d,
 517                    int ftype, int fsubtype, t_params plist[])
 518 {
 519     t_symtab       stab;
 520     gpp_atomtype_t atype;
 521     t_param       *param;
 522     t_atom        *a;
 523     int            i;
 524
 525     atype = init_atomtype();
 526     snew(a, 1);
 527     snew(param, 1);
 528     open_symtab(&stab);
 529     for (i = 0; (i < natoms); i++)
 530     {
 531         char buf[12];
 532         sprintf(buf, "%4d", (i+1));
 533         add_atomtype(atype, &stab, a, buf, param, 0, 0, 0, 0, 0, 0, 0);
 534     }
 535     print_bt(out, d, atype, ftype, fsubtype, plist, TRUE);
 536
 537     done_symtab(&stab);
 538     sfree(a);
 539     sfree(param);
 540     done_atomtype(atype);
 541 }