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