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