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