src/gromacs/gmxpreprocess/toputil.c

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