src/kernel/grompp.c

   1 /*  -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
   2  *
   3  *
   4  *                This source code is part of
   5  *
   6  *                 G   R   O   M   A   C   S
   7  *
   8  *          GROningen MAchine for Chemical Simulations
   9  *
  10  *                        VERSION 3.2.03
  11  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  12  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  13  * Copyright (c) 2001-2004, The GROMACS development team,
  14  * check out http://www.gromacs.org for more information.
  15
  16  * This program is free software; you can redistribute it and/or
  17  * modify it under the terms of the GNU General Public License
  18  * as published by the Free Software Foundation; either version 2
  19  * of the License, or (at your option) any later version.
  20  *
  21  * If you want to redistribute modifications, please consider that
  22  * scientific software is very special. Version control is crucial -
  23  * bugs must be traceable. We will be happy to consider code for
  24  * inclusion in the official distribution, but derived work must not
  25  * be called official GROMACS. Details are found in the README & COPYING
  26  * files - if they are missing, get the official version at www.gromacs.org.
  27  *
  28  * To help us fund GROMACS development, we humbly ask that you cite
  29  * the papers on the package - you can find them in the top README file.
  30  *
  31  * For more info, check our website at http://www.gromacs.org
  32  *
  33  * And Hey:
  34  * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
  35  */
  36 #ifdef HAVE_CONFIG_H
  37 #include <config.h>
  38 #endif
  39
  40 #include <sys/types.h>
  41 #include <math.h>
  42 #include <string.h>
  43 #include <errno.h>
  44 #include <limits.h>
  45
  46 #include "sysstuff.h"
  47 #include "smalloc.h"
  48 #include "macros.h"
  49 #include "string2.h"
  50 #include "readir.h"
  51 #include "toputil.h"
  52 #include "topio.h"
  53 #include "confio.h"
  54 #include "copyrite.h"
  55 #include "readir.h"
  56 #include "symtab.h"
  57 #include "names.h"
  58 #include "grompp.h"
  59 #include "random.h"
  60 #include "vec.h"
  61 #include "futil.h"
  62 #include "statutil.h"
  63 #include "splitter.h"
  64 #include "sortwater.h"
  65 #include "convparm.h"
  66 #include "gmx_fatal.h"
  67 #include "warninp.h"
  68 #include "index.h"
  69 #include "gmxfio.h"
  70 #include "trnio.h"
  71 #include "tpxio.h"
  72 #include "vsite_parm.h"
  73 #include "txtdump.h"
  74 #include "calcgrid.h"
  75 #include "add_par.h"
  76 #include "enxio.h"
  77 #include "perf_est.h"
  78 #include "compute_io.h"
  79 #include "gpp_atomtype.h"
  80 #include "gpp_tomorse.h"
  81 #include "mtop_util.h"
  82 #include "genborn.h"
  83
  84 static int rm_interactions(int ifunc,int nrmols,t_molinfo mols[])
  85 {
  86   int  i,n;
  87
  88   n=0;
  89   /* For all the molecule types */
  90   for(i=0; i<nrmols; i++) {
  91     n += mols[i].plist[ifunc].nr;
  92     mols[i].plist[ifunc].nr=0;
  93   }
  94   return n;
  95 }
  96
  97 static int check_atom_names(const char *fn1, const char *fn2,
  98                             gmx_mtop_t *mtop, t_atoms *at)
  99 {
 100   int mb,m,i,j,nmismatch;
 101   t_atoms *tat;
 102 #define MAXMISMATCH 20
 103
 104   if (mtop->natoms != at->nr)
 105     gmx_incons("comparing atom names");
 106
 107   nmismatch=0;
 108   i = 0;
 109   for(mb=0; mb<mtop->nmolblock; mb++) {
 110     tat = &mtop->moltype[mtop->molblock[mb].type].atoms;
 111     for(m=0; m<mtop->molblock[mb].nmol; m++) {
 112       for(j=0; j < tat->nr; j++) {
 113         if (strcmp( *(tat->atomname[j]) , *(at->atomname[i]) ) != 0) {
 114           if (nmismatch < MAXMISMATCH) {
 115             fprintf(stderr,
 116                     "Warning: atom name %d in %s and %s does not match (%s - %s)\n",
 117                     i+1, fn1, fn2, *(tat->atomname[j]), *(at->atomname[i]));
 118           } else if (nmismatch == MAXMISMATCH) {
 119             fprintf(stderr,"(more than %d non-matching atom names)\n",MAXMISMATCH);
 120           }
 121           nmismatch++;
 122         }
 123         i++;
 124       }
 125     }
 126   }
 127
 128   return nmismatch;
 129 }
 130
 131 static void check_eg_vs_cg(gmx_mtop_t *mtop)
 132 {
 133   int astart,mb,m,cg,j,firstj;
 134   unsigned char firsteg,eg;
 135   gmx_moltype_t *molt;
 136
 137   /* Go through all the charge groups and make sure all their
 138    * atoms are in the same energy group.
 139    */
 140
 141   astart = 0;
 142   for(mb=0; mb<mtop->nmolblock; mb++) {
 143     molt = &mtop->moltype[mtop->molblock[mb].type];
 144     for(m=0; m<mtop->molblock[mb].nmol; m++) {
 145       for(cg=0; cg<molt->cgs.nr;cg++) {
 146         /* Get the energy group of the first atom in this charge group */
 147         firstj = astart + molt->cgs.index[cg];
 148         firsteg = ggrpnr(&mtop->groups,egcENER,firstj);
 149         for(j=molt->cgs.index[cg]+1;j<molt->cgs.index[cg+1];j++) {
 150           eg = ggrpnr(&mtop->groups,egcENER,astart+j);
 151           if(eg != firsteg) {
 152             gmx_fatal(FARGS,"atoms %d and %d in charge group %d of molecule type '%s' are in different energy groups",
 153                       firstj+1,astart+j+1,cg+1,*molt->name);
 154           }
 155         }
 156       }
 157       astart += molt->atoms.nr;
 158     }
 159   }
 160 }
 161
 162 static void check_cg_sizes(const char *topfn,t_block *cgs,warninp_t wi)
 163 {
 164     int  maxsize,cg;
 165     char warn_buf[STRLEN];
 166
 167     maxsize = 0;
 168     for(cg=0; cg<cgs->nr; cg++)
 169     {
 170         maxsize = max(maxsize,cgs->index[cg+1]-cgs->index[cg]);
 171     }
 172
 173     if (maxsize > MAX_CHARGEGROUP_SIZE)
 174     {
 175         gmx_fatal(FARGS,"The largest charge group contains %d atoms. The maximum is %d.",maxsize,MAX_CHARGEGROUP_SIZE);
 176     }
 177     else if (maxsize > 10)
 178     {
 179         set_warning_line(wi,topfn,-1);
 180         sprintf(warn_buf,
 181                 "The largest charge group contains %d atoms.\n"
 182                 "Since atoms only see each other when the centers of geometry of the charge groups they belong to are within the cut-off distance, too large charge groups can lead to serious cut-off artifacts.\n"
 183                 "For efficiency and accuracy, charge group should consist of a few atoms.\n"
 184                 "For all-atom force fields use: CH3, CH2, CH, NH2, NH, OH, CO2, CO, etc.",
 185                 maxsize);
 186         warning_note(wi,warn_buf);
 187     }
 188 }
 189
 190 static void check_bonds_timestep(gmx_mtop_t *mtop,double dt,warninp_t wi)
 191 {
 192     /* This check is not intended to ensure accurate integration,
 193      * rather it is to signal mistakes in the mdp settings.
 194      * A common mistake is to forget to turn on constraints
 195      * for MD after energy minimization with flexible bonds.
 196      * This check can also detect too large time steps for flexible water
 197      * models, but such errors will often be masked by the constraints
 198      * mdp options, which turns flexible water into water with bond constraints,
 199      * but without an angle constraint. Unfortunately such incorrect use
 200      * of water models can not easily be detected without checking
 201      * for specific model names.
 202      *
 203      * The stability limit of leap-frog or velocity verlet is 4.44 steps
 204      * per oscillational period.
 205      * But accurate bonds distributions are lost far before that limit.
 206      * To allow relatively common schemes (although not common with Gromacs)
 207      * of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
 208      * we set the note limit to 10.
 209      */
 210     int       min_steps_warn=5;
 211     int       min_steps_note=10;
 212     t_iparams *ip;
 213     int       molt;
 214     gmx_moltype_t *moltype,*w_moltype;
 215     t_atom    *atom;
 216     t_ilist   *ilist,*ilb,*ilc,*ils;
 217     int       ftype;
 218     int       i,a1,a2,w_a1,w_a2,j;
 219     real      twopi2,limit2,fc,re,m1,m2,period2,w_period2;
 220     gmx_bool  bFound,bWater,bWarn;
 221     char      warn_buf[STRLEN];
 222
 223     ip = mtop->ffparams.iparams;
 224
 225     twopi2 = sqr(2*M_PI);
 226
 227     limit2 = sqr(min_steps_note*dt);
 228
 229     w_a1 = w_a2 = -1;
 230     w_period2 = -1.0;
 231
 232     w_moltype = NULL;
 233     for(molt=0; molt<mtop->nmoltype; molt++)
 234     {
 235         moltype = &mtop->moltype[molt];
 236         atom  = moltype->atoms.atom;
 237         ilist = moltype->ilist;
 238         ilc = &ilist[F_CONSTR];
 239         ils = &ilist[F_SETTLE];
 240         for(ftype=0; ftype<F_NRE; ftype++)
 241         {
 242             if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
 243             {
 244                 continue;
 245             }
 246
 247             ilb = &ilist[ftype];
 248             for(i=0; i<ilb->nr; i+=3)
 249             {
 250                 fc = ip[ilb->iatoms[i]].harmonic.krA;
 251                 re = ip[ilb->iatoms[i]].harmonic.rA;
 252                 if (ftype == F_G96BONDS)
 253                 {
 254                     /* Convert squared sqaure fc to harmonic fc */
 255                     fc = 2*fc*re;
 256                 }
 257                 a1 = ilb->iatoms[i+1];
 258                 a2 = ilb->iatoms[i+2];
 259                 m1 = atom[a1].m;
 260                 m2 = atom[a2].m;
 261                 if (fc > 0 && m1 > 0 && m2 > 0)
 262                 {
 263                     period2 = twopi2*m1*m2/((m1 + m2)*fc);
 264                 }
 265                 else
 266                 {
 267                     period2 = GMX_FLOAT_MAX;
 268                 }
 269                 if (debug)
 270                 {
 271                     fprintf(debug,"fc %g m1 %g m2 %g period %g\n",
 272                             fc,m1,m2,sqrt(period2));
 273                 }
 274                 if (period2 < limit2)
 275                 {
 276                     bFound = FALSE;
 277                     for(j=0; j<ilc->nr; j+=3)
 278                     {
 279                         if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) ||
 280                             (ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1))
 281                             {
 282                                 bFound = TRUE;
 283                             }
 284                         }
 285                     for(j=0; j<ils->nr; j+=2)
 286                     {
 287                         if ((a1 >= ils->iatoms[j+1] && a1 < ils->iatoms[j+1]+3) &&
 288                             (a2 >= ils->iatoms[j+1] && a2 < ils->iatoms[j+1]+3))
 289                         {
 290                             bFound = TRUE;
 291                         }
 292                     }
 293                     if (!bFound &&
 294                         (w_moltype == NULL || period2 < w_period2))
 295                     {
 296                         w_moltype = moltype;
 297                         w_a1      = a1;
 298                         w_a2      = a2;
 299                         w_period2 = period2;
 300                     }
 301                 }
 302             }
 303         }
 304     }
 305
 306     if (w_moltype != NULL)
 307     {
 308         bWarn = (w_period2 < sqr(min_steps_warn*dt));
 309         /* A check that would recognize most water models */
 310         bWater = ((*w_moltype->atoms.atomname[0])[0] == 'O' &&
 311                   w_moltype->atoms.nr <= 5);
 312         sprintf(warn_buf,"The bond in molecule-type %s between atoms %d %s and %d %s has an estimated oscillational period of %.1e ps, which is less than %d times the time step of %.1e ps.\n"
 313                 "%s",
 314                 *w_moltype->name,
 315                 w_a1+1,*w_moltype->atoms.atomname[w_a1],
 316                 w_a2+1,*w_moltype->atoms.atomname[w_a2],
 317                 sqrt(w_period2),bWarn ? min_steps_warn : min_steps_note,dt,
 318                 bWater ?
 319                 "Maybe you asked for fexible water." :
 320                 "Maybe you forgot to change the constraints mdp option.");
 321         if (bWarn)
 322         {
 323             warning(wi,warn_buf);
 324         }
 325         else
 326         {
 327             warning_note(wi,warn_buf);
 328         }
 329     }
 330 }
 331
 332 static void check_vel(gmx_mtop_t *mtop,rvec v[])
 333 {
 334   gmx_mtop_atomloop_all_t aloop;
 335   t_atom *atom;
 336   int a;
 337
 338   aloop = gmx_mtop_atomloop_all_init(mtop);
 339   while (gmx_mtop_atomloop_all_next(aloop,&a,&atom)) {
 340     if (atom->ptype == eptShell ||
 341         atom->ptype == eptBond  ||
 342         atom->ptype == eptVSite) {
 343       clear_rvec(v[a]);
 344     }
 345   }
 346 }
 347
 348 static gmx_bool nint_ftype(gmx_mtop_t *mtop,t_molinfo *mi,int ftype)
 349 {
 350   int nint,mb;
 351
 352   nint = 0;
 353   for(mb=0; mb<mtop->nmolblock; mb++) {
 354     nint += mtop->molblock[mb].nmol*mi[mtop->molblock[mb].type].plist[ftype].nr;
 355   }
 356
 357   return nint;
 358 }
 359
 360 /* This routine reorders the molecule type array
 361  * in the order of use in the molblocks,
 362  * unused molecule types are deleted.
 363  */
 364 static void renumber_moltypes(gmx_mtop_t *sys,
 365                               int *nmolinfo,t_molinfo **molinfo)
 366 {
 367   int *order,norder,i;
 368   int mb,mi;
 369   t_molinfo *minew;
 370
 371   snew(order,*nmolinfo);
 372   norder = 0;
 373   for(mb=0; mb<sys->nmolblock; mb++) {
 374     for(i=0; i<norder; i++) {
 375       if (order[i] == sys->molblock[mb].type) {
 376         break;
 377       }
 378     }
 379     if (i == norder) {
 380       /* This type did not occur yet, add it */
 381       order[norder] = sys->molblock[mb].type;
 382       /* Renumber the moltype in the topology */
 383       norder++;
 384     }
 385     sys->molblock[mb].type = i;
 386   }
 387
 388   /* We still need to reorder the molinfo structs */
 389   snew(minew,norder);
 390   for(mi=0; mi<*nmolinfo; mi++) {
 391     for(i=0; i<norder; i++) {
 392       if (order[i] == mi) {
 393         break;
 394       }
 395     }
 396     if (i == norder) {
 397       done_mi(&(*molinfo)[mi]);
 398     } else {
 399       minew[i] = (*molinfo)[mi];
 400     }
 401   }
 402   sfree(*molinfo);
 403
 404   *nmolinfo = norder;
 405   *molinfo  = minew;
 406 }
 407
 408 static void molinfo2mtop(int nmi,t_molinfo *mi,gmx_mtop_t *mtop)
 409 {
 410   int m;
 411   gmx_moltype_t *molt;
 412
 413   mtop->nmoltype = nmi;
 414   snew(mtop->moltype,nmi);
 415   for(m=0; m<nmi; m++) {
 416     molt = &mtop->moltype[m];
 417     molt->name  = mi[m].name;
 418     molt->atoms = mi[m].atoms;
 419     /* ilists are copied later */
 420     molt->cgs   = mi[m].cgs;
 421     molt->excls = mi[m].excls;
 422   }
 423 }
 424
 425 static void
 426 new_status(const char *topfile,const char *topppfile,const char *confin,
 427            t_gromppopts *opts,t_inputrec *ir,gmx_bool bZero,
 428            gmx_bool bGenVel,gmx_bool bVerbose,t_state *state,
 429            gpp_atomtype_t atype,gmx_mtop_t *sys,
 430            int *nmi,t_molinfo **mi,t_params plist[],
 431            int *comb,double *reppow,real *fudgeQQ,
 432            gmx_bool bMorse,
 433            warninp_t wi)
 434 {
 435   t_molinfo   *molinfo=NULL;
 436   int         nmolblock;
 437   gmx_molblock_t *molblock,*molbs;
 438   t_atoms     *confat;
 439   int         mb,i,nrmols,nmismatch;
 440   char        buf[STRLEN];
 441   gmx_bool        bGB=FALSE;
 442   char        warn_buf[STRLEN];
 443
 444   init_mtop(sys);
 445
 446   /* Set gmx_boolean for GB */
 447   if(ir->implicit_solvent)
 448     bGB=TRUE;
 449
 450   /* TOPOLOGY processing */
 451   sys->name = do_top(bVerbose,topfile,topppfile,opts,bZero,&(sys->symtab),
 452                      plist,comb,reppow,fudgeQQ,
 453                      atype,&nrmols,&molinfo,ir,
 454                      &nmolblock,&molblock,bGB,
 455                      wi);
 456
 457   sys->nmolblock = 0;
 458   snew(sys->molblock,nmolblock);
 459
 460   sys->natoms = 0;
 461   for(mb=0; mb<nmolblock; mb++) {
 462     if (sys->nmolblock > 0 &&
 463         molblock[mb].type == sys->molblock[sys->nmolblock-1].type) {
 464       /* Merge consecutive blocks with the same molecule type */
 465       sys->molblock[sys->nmolblock-1].nmol += molblock[mb].nmol;
 466       sys->natoms += molblock[mb].nmol*sys->molblock[sys->nmolblock-1].natoms_mol;
 467     } else if (molblock[mb].nmol > 0) {
 468       /* Add a new molblock to the topology */
 469       molbs = &sys->molblock[sys->nmolblock];
 470       *molbs = molblock[mb];
 471       molbs->natoms_mol = molinfo[molbs->type].atoms.nr;
 472       molbs->nposres_xA = 0;
 473       molbs->nposres_xB = 0;
 474       sys->natoms += molbs->nmol*molbs->natoms_mol;
 475       sys->nmolblock++;
 476     }
 477   }
 478   if (sys->nmolblock == 0) {
 479     gmx_fatal(FARGS,"No molecules were defined in the system");
 480   }
 481
 482   renumber_moltypes(sys,&nrmols,&molinfo);
 483
 484   if (bMorse)
 485     convert_harmonics(nrmols,molinfo,atype);
 486
 487   if (ir->eDisre == edrNone) {
 488     i = rm_interactions(F_DISRES,nrmols,molinfo);
 489     if (i > 0) {
 490       set_warning_line(wi,"unknown",-1);
 491       sprintf(warn_buf,"disre = no, removed %d distance restraints",i);
 492       warning_note(wi,warn_buf);
 493     }
 494   }
 495   if (opts->bOrire == FALSE) {
 496     i = rm_interactions(F_ORIRES,nrmols,molinfo);
 497     if (i > 0) {
 498       set_warning_line(wi,"unknown",-1);
 499       sprintf(warn_buf,"orire = no, removed %d orientation restraints",i);
 500       warning_note(wi,warn_buf);
 501     }
 502   }
 503   if (opts->bDihre == FALSE) {
 504     i = rm_interactions(F_DIHRES,nrmols,molinfo);
 505     if (i > 0) {
 506       set_warning_line(wi,"unknown",-1);
 507       sprintf(warn_buf,"dihre = no, removed %d dihedral restraints",i);
 508       warning_note(wi,warn_buf);
 509     }
 510   }
 511
 512   /* Copy structures from msys to sys */
 513   molinfo2mtop(nrmols,molinfo,sys);
 514
 515   gmx_mtop_finalize(sys);
 516
 517   /* COORDINATE file processing */
 518   if (bVerbose)
 519     fprintf(stderr,"processing coordinates...\n");
 520
 521   get_stx_coordnum(confin,&state->natoms);
 522   if (state->natoms != sys->natoms)
 523     gmx_fatal(FARGS,"number of coordinates in coordinate file (%s, %d)\n"
 524                 "             does not match topology (%s, %d)",
 525               confin,state->natoms,topfile,sys->natoms);
 526   else {
 527     /* make space for coordinates and velocities */
 528     char title[STRLEN];
 529     snew(confat,1);
 530     init_t_atoms(confat,state->natoms,FALSE);
 531     init_state(state,state->natoms,0,0,0);
 532     read_stx_conf(confin,title,confat,state->x,state->v,NULL,state->box);
 533     /* This call fixes the box shape for runs with pressure scaling */
 534     set_box_rel(ir,state);
 535
 536     nmismatch = check_atom_names(topfile, confin, sys, confat);
 537     free_t_atoms(confat,TRUE);
 538     sfree(confat);
 539
 540     if (nmismatch) {
 541       sprintf(buf,"%d non-matching atom name%s\n"
 542               "atom names from %s will be used\n"
 543               "atom names from %s will be ignored\n",
 544               nmismatch,(nmismatch == 1) ? "" : "s",topfile,confin);
 545       warning(wi,buf);
 546     }
 547     if (bVerbose)
 548       fprintf(stderr,"double-checking input for internal consistency...\n");
 549     double_check(ir,state->box,nint_ftype(sys,molinfo,F_CONSTR),wi);
 550   }
 551
 552   if (bGenVel) {
 553     real *mass;
 554     gmx_mtop_atomloop_all_t aloop;
 555     t_atom *atom;
 556
 557     snew(mass,state->natoms);
 558     aloop = gmx_mtop_atomloop_all_init(sys);
 559     while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
 560       mass[i] = atom->m;
 561     }
 562
 563     if (opts->seed == -1) {
 564       opts->seed = make_seed();
 565       fprintf(stderr,"Setting gen_seed to %d\n",opts->seed);
 566     }
 567     maxwell_speed(opts->tempi,opts->seed,sys,state->v);
 568
 569     stop_cm(stdout,state->natoms,mass,state->x,state->v);
 570     sfree(mass);
 571   }
 572
 573   *nmi = nrmols;
 574   *mi  = molinfo;
 575 }
 576
 577 static void copy_state(const char *slog,t_trxframe *fr,
 578                        gmx_bool bReadVel,t_state *state,
 579                        double *use_time)
 580 {
 581     int i;
 582
 583     if (fr->not_ok & FRAME_NOT_OK)
 584     {
 585         gmx_fatal(FARGS,"Can not start from an incomplete frame");
 586     }
 587     if (!fr->bX)
 588     {
 589         gmx_fatal(FARGS,"Did not find a frame with coordinates in file %s",
 590                   slog);
 591     }
 592
 593     for(i=0; i<state->natoms; i++)
 594     {
 595         copy_rvec(fr->x[i],state->x[i]);
 596     }
 597     if (bReadVel)
 598     {
 599         if (!fr->bV)
 600         {
 601             gmx_incons("Trajecory frame unexpectedly does not contain velocities");
 602         }
 603         for(i=0; i<state->natoms; i++)
 604         {
 605             copy_rvec(fr->v[i],state->v[i]);
 606         }
 607     }
 608     if (fr->bBox)
 609     {
 610         copy_mat(fr->box,state->box);
 611     }
 612
 613     *use_time = fr->time;
 614 }
 615
 616 static void cont_status(const char *slog,const char *ener,
 617                         gmx_bool bNeedVel,gmx_bool bGenVel, real fr_time,
 618                         t_inputrec *ir,t_state *state,
 619                         gmx_mtop_t *sys,
 620                         const output_env_t oenv)
 621      /* If fr_time == -1 read the last frame available which is complete */
 622 {
 623     gmx_bool bReadVel;
 624     t_trxframe  fr;
 625     t_trxstatus *fp;
 626     int i;
 627     double use_time;
 628
 629     bReadVel = (bNeedVel && !bGenVel);
 630
 631     fprintf(stderr,
 632             "Reading Coordinates%s and Box size from old trajectory\n",
 633             bReadVel ? ", Velocities" : "");
 634     if (fr_time == -1)
 635     {
 636         fprintf(stderr,"Will read whole trajectory\n");
 637     }
 638     else
 639     {
 640         fprintf(stderr,"Will read till time %g\n",fr_time);
 641     }
 642     if (!bReadVel)
 643     {
 644         if (bGenVel)
 645         {
 646             fprintf(stderr,"Velocities generated: "
 647                     "ignoring velocities in input trajectory\n");
 648         }
 649         read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X);
 650     }
 651     else
 652     {
 653         read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X | TRX_NEED_V);
 654
 655         if (!fr.bV)
 656         {
 657             fprintf(stderr,
 658                     "\n"
 659                     "WARNING: Did not find a frame with velocities in file %s,\n"
 660                     "         all velocities will be set to zero!\n\n",slog);
 661             for(i=0; i<sys->natoms; i++)
 662             {
 663                 clear_rvec(state->v[i]);
 664             }
 665             close_trj(fp);
 666             /* Search for a frame without velocities */
 667             bReadVel = FALSE;
 668             read_first_frame(oenv,&fp,slog,&fr,TRX_NEED_X);
 669         }
 670     }
 671
 672     state->natoms = fr.natoms;
 673
 674     if (sys->natoms != state->natoms)
 675     {
 676         gmx_fatal(FARGS,"Number of atoms in Topology "
 677                   "is not the same as in Trajectory");
 678     }
 679     copy_state(slog,&fr,bReadVel,state,&use_time);
 680
 681     /* Find the appropriate frame */
 682     while ((fr_time == -1 || fr.time < fr_time) &&
 683            read_next_frame(oenv,fp,&fr))
 684     {
 685         copy_state(slog,&fr,bReadVel,state,&use_time);
 686     }
 687
 688     close_trj(fp);
 689
 690     /* Set the relative box lengths for preserving the box shape.
 691      * Note that this call can lead to differences in the last bit
 692      * with respect to using tpbconv to create a [TT].tpx[tt] file.
 693      */
 694     set_box_rel(ir,state);
 695
 696     fprintf(stderr,"Using frame at t = %g ps\n",use_time);
 697     fprintf(stderr,"Starting time for run is %g ps\n",ir->init_t);
 698
 699     if ((ir->epc != epcNO  || ir->etc ==etcNOSEHOOVER) && ener)
 700     {
 701         get_enx_state(ener,use_time,&sys->groups,ir,state);
 702         preserve_box_shape(ir,state->box_rel,state->boxv);
 703     }
 704 }
 705
 706 static void read_posres(gmx_mtop_t *mtop,t_molinfo *molinfo,gmx_bool bTopB,
 707                         char *fn,
 708                         int rc_scaling, int ePBC,
 709                         rvec com,
 710                         warninp_t wi)
 711 {
 712   gmx_bool   bFirst = TRUE;
 713   rvec   *x,*v,*xp;
 714   dvec   sum;
 715   double totmass;
 716   t_atoms dumat;
 717   matrix box,invbox;
 718   int    natoms,npbcdim=0;
 719   char   warn_buf[STRLEN],title[STRLEN];
 720   int    a,i,ai,j,k,mb,nat_molb;
 721   gmx_molblock_t *molb;
 722   t_params *pr;
 723   t_atom *atom;
 724
 725   get_stx_coordnum(fn,&natoms);
 726   if (natoms != mtop->natoms) {
 727     sprintf(warn_buf,"The number of atoms in %s (%d) does not match the number of atoms in the topology (%d). Will assume that the first %d atoms in the topology and %s match.",fn,natoms,mtop->natoms,min(mtop->natoms,natoms),fn);
 728     warning(wi,warn_buf);
 729   }
 730   snew(x,natoms);
 731   snew(v,natoms);
 732   init_t_atoms(&dumat,natoms,FALSE);
 733   read_stx_conf(fn,title,&dumat,x,v,NULL,box);
 734
 735   npbcdim = ePBC2npbcdim(ePBC);
 736   clear_rvec(com);
 737   if (rc_scaling != erscNO) {
 738     copy_mat(box,invbox);
 739     for(j=npbcdim; j<DIM; j++) {
 740       clear_rvec(invbox[j]);
 741       invbox[j][j] = 1;
 742     }
 743     m_inv_ur0(invbox,invbox);
 744   }
 745
 746   /* Copy the reference coordinates to mtop */
 747   clear_dvec(sum);
 748   totmass = 0;
 749   a = 0;
 750   for(mb=0; mb<mtop->nmolblock; mb++) {
 751     molb = &mtop->molblock[mb];
 752     nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
 753     pr = &(molinfo[molb->type].plist[F_POSRES]);
 754     if (pr->nr > 0) {
 755       atom = mtop->moltype[molb->type].atoms.atom;
 756       for(i=0; (i<pr->nr); i++) {
 757         ai=pr->param[i].AI;
 758         if (ai >= natoms) {
 759           gmx_fatal(FARGS,"Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n",
 760                     ai+1,*molinfo[molb->type].name,fn,natoms);
 761         }
 762         if (rc_scaling == erscCOM) {
 763           /* Determine the center of mass of the posres reference coordinates */
 764           for(j=0; j<npbcdim; j++) {
 765             sum[j] += atom[ai].m*x[a+ai][j];
 766           }
 767           totmass  += atom[ai].m;
 768         }
 769       }
 770       if (!bTopB) {
 771         molb->nposres_xA = nat_molb;
 772         snew(molb->posres_xA,molb->nposres_xA);
 773         for(i=0; i<nat_molb; i++) {
 774           copy_rvec(x[a+i],molb->posres_xA[i]);
 775         }
 776       } else {
 777         molb->nposres_xB = nat_molb;
 778         snew(molb->posres_xB,molb->nposres_xB);
 779         for(i=0; i<nat_molb; i++) {
 780           copy_rvec(x[a+i],molb->posres_xB[i]);
 781         }
 782       }
 783     }
 784     a += nat_molb;
 785   }
 786   if (rc_scaling == erscCOM) {
 787     if (totmass == 0)
 788       gmx_fatal(FARGS,"The total mass of the position restraint atoms is 0");
 789     for(j=0; j<npbcdim; j++)
 790       com[j] = sum[j]/totmass;
 791     fprintf(stderr,"The center of mass of the position restraint coord's is %6.3f %6.3f %6.3f\n",com[XX],com[YY],com[ZZ]);
 792   }
 793
 794   if (rc_scaling != erscNO) {
 795     for(mb=0; mb<mtop->nmolblock; mb++) {
 796       molb = &mtop->molblock[mb];
 797       nat_molb = molb->nmol*mtop->moltype[molb->type].atoms.nr;
 798       if (molb->nposres_xA > 0 || molb->nposres_xB > 0) {
 799         xp = (!bTopB ? molb->posres_xA : molb->posres_xB);
 800         for(i=0; i<nat_molb; i++) {
 801           for(j=0; j<npbcdim; j++) {
 802             if (rc_scaling == erscALL) {
 803               /* Convert from Cartesian to crystal coordinates */
 804               xp[i][j] *= invbox[j][j];
 805               for(k=j+1; k<npbcdim; k++) {
 806                 xp[i][j] += invbox[k][j]*xp[i][k];
 807               }
 808             } else if (rc_scaling == erscCOM) {
 809               /* Subtract the center of mass */
 810               xp[i][j] -= com[j];
 811             }
 812           }
 813         }
 814       }
 815     }
 816
 817     if (rc_scaling == erscCOM) {
 818       /* Convert the COM from Cartesian to crystal coordinates */
 819       for(j=0; j<npbcdim; j++) {
 820         com[j] *= invbox[j][j];
 821         for(k=j+1; k<npbcdim; k++) {
 822           com[j] += invbox[k][j]*com[k];
 823         }
 824       }
 825     }
 826   }
 827
 828   free_t_atoms(&dumat,TRUE);
 829   sfree(x);
 830   sfree(v);
 831 }
 832
 833 static void gen_posres(gmx_mtop_t *mtop,t_molinfo *mi,
 834                        char *fnA, char *fnB,
 835                        int rc_scaling, int ePBC,
 836                        rvec com, rvec comB,
 837                        warninp_t wi)
 838 {
 839   int i,j;
 840
 841   read_posres  (mtop,mi,FALSE,fnA,rc_scaling,ePBC,com,wi);
 842   if (strcmp(fnA,fnB) != 0) {
 843       read_posres(mtop,mi,TRUE ,fnB,rc_scaling,ePBC,comB,wi);
 844   }
 845 }
 846
 847 static void set_wall_atomtype(gpp_atomtype_t at,t_gromppopts *opts,
 848                               t_inputrec *ir)
 849 {
 850   int i;
 851
 852   if (ir->nwall > 0)
 853     fprintf(stderr,"Searching the wall atom type(s)\n");
 854   for(i=0; i<ir->nwall; i++)
 855     ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i],at);
 856 }
 857
 858 static int nrdf_internal(t_atoms *atoms)
 859 {
 860   int i,nmass,nrdf;
 861
 862   nmass = 0;
 863   for(i=0; i<atoms->nr; i++) {
 864     /* Vsite ptype might not be set here yet, so also check the mass */
 865     if ((atoms->atom[i].ptype == eptAtom ||
 866          atoms->atom[i].ptype == eptNucleus)
 867         && atoms->atom[i].m > 0) {
 868       nmass++;
 869     }
 870   }
 871   switch (nmass) {
 872   case 0:  nrdf = 0; break;
 873   case 1:  nrdf = 0; break;
 874   case 2:  nrdf = 1; break;
 875   default: nrdf = nmass*3 - 6; break;
 876   }
 877
 878   return nrdf;
 879 }
 880
 881 void
 882 spline1d( double        dx,
 883                  double *      y,
 884                  int           n,
 885                  double *      u,
 886                  double *      y2 )
 887 {
 888     int i;
 889     double p,q;
 890
 891     y2[0] = 0.0;
 892     u[0]  = 0.0;
 893
 894     for(i=1;i<n-1;i++)
 895     {
 896                 p = 0.5*y2[i-1]+2.0;
 897         y2[i] = -0.5/p;
 898         q = (y[i+1]-2.0*y[i]+y[i-1])/dx;
 899                 u[i] = (3.0*q/dx-0.5*u[i-1])/p;
 900     }
 901
 902     y2[n-1] = 0.0;
 903
 904     for(i=n-2;i>=0;i--)
 905     {
 906         y2[i] = y2[i]*y2[i+1]+u[i];
 907     }
 908 }
 909
 910
 911 void
 912 interpolate1d( double     xmin,
 913                           double     dx,
 914                           double *   ya,
 915                           double *   y2a,
 916                           double     x,
 917                           double *   y,
 918                           double *   y1)
 919 {
 920     int ix;
 921     double a,b;
 922
 923     ix = (x-xmin)/dx;
 924
 925     a = (xmin+(ix+1)*dx-x)/dx;
 926     b = (x-xmin-ix*dx)/dx;
 927
 928     *y  = a*ya[ix]+b*ya[ix+1]+((a*a*a-a)*y2a[ix]+(b*b*b-b)*y2a[ix+1])*(dx*dx)/6.0;
 929     *y1 = (ya[ix+1]-ya[ix])/dx-(3.0*a*a-1.0)/6.0*dx*y2a[ix]+(3.0*b*b-1.0)/6.0*dx*y2a[ix+1];
 930 }
 931
 932
 933 void
 934 setup_cmap (int              grid_spacing,
 935                         int              nc,
 936                         real *           grid ,
 937                         gmx_cmap_t *     cmap_grid)
 938 {
 939         double *tmp_u,*tmp_u2,*tmp_yy,*tmp_y1,*tmp_t2,*tmp_grid;
 940
 941     int    i,j,k,ii,jj,kk,idx;
 942         int    offset;
 943     double dx,xmin,v,v1,v2,v12;
 944     double phi,psi;
 945
 946         snew(tmp_u,2*grid_spacing);
 947         snew(tmp_u2,2*grid_spacing);
 948         snew(tmp_yy,2*grid_spacing);
 949         snew(tmp_y1,2*grid_spacing);
 950         snew(tmp_t2,2*grid_spacing*2*grid_spacing);
 951         snew(tmp_grid,2*grid_spacing*2*grid_spacing);
 952
 953     dx = 360.0/grid_spacing;
 954     xmin = -180.0-dx*grid_spacing/2;
 955
 956         for(kk=0;kk<nc;kk++)
 957         {
 958                 /* Compute an offset depending on which cmap we are using
 959                  * Offset will be the map number multiplied with the grid_spacing * grid_spacing * 2
 960                  */
 961                 offset = kk * grid_spacing * grid_spacing * 2;
 962
 963                 for(i=0;i<2*grid_spacing;i++)
 964                 {
 965                         ii=(i+grid_spacing-grid_spacing/2)%grid_spacing;
 966
 967                         for(j=0;j<2*grid_spacing;j++)
 968                         {
 969                                 jj=(j+grid_spacing-grid_spacing/2)%grid_spacing;
 970                                 tmp_grid[i*grid_spacing*2+j] = grid[offset+ii*grid_spacing+jj];
 971                         }
 972                 }
 973
 974                 for(i=0;i<2*grid_spacing;i++)
 975                 {
 976                         spline1d(dx,&(tmp_grid[2*grid_spacing*i]),2*grid_spacing,tmp_u,&(tmp_t2[2*grid_spacing*i]));
 977                 }
 978
 979                 for(i=grid_spacing/2;i<grid_spacing+grid_spacing/2;i++)
 980                 {
 981                         ii = i-grid_spacing/2;
 982                         phi = ii*dx-180.0;
 983
 984                         for(j=grid_spacing/2;j<grid_spacing+grid_spacing/2;j++)
 985                         {
 986                                 jj = j-grid_spacing/2;
 987                                 psi = jj*dx-180.0;
 988
 989                                 for(k=0;k<2*grid_spacing;k++)
 990                                 {
 991                                         interpolate1d(xmin,dx,&(tmp_grid[2*grid_spacing*k]),
 992                                                                   &(tmp_t2[2*grid_spacing*k]),psi,&tmp_yy[k],&tmp_y1[k]);
 993                                 }
 994
 995                                 spline1d(dx,tmp_yy,2*grid_spacing,tmp_u,tmp_u2);
 996                                 interpolate1d(xmin,dx,tmp_yy,tmp_u2,phi,&v,&v1);
 997                                 spline1d(dx,tmp_y1,2*grid_spacing,tmp_u,tmp_u2);
 998                                 interpolate1d(xmin,dx,tmp_y1,tmp_u2,phi,&v2,&v12);
 999
1000                                 idx = ii*grid_spacing+jj;
1001                                 cmap_grid->cmapdata[kk].cmap[idx*4] = grid[offset+ii*grid_spacing+jj];
1002                                 cmap_grid->cmapdata[kk].cmap[idx*4+1] = v1;
1003                                 cmap_grid->cmapdata[kk].cmap[idx*4+2] = v2;
1004                                 cmap_grid->cmapdata[kk].cmap[idx*4+3] = v12;
1005                         }
1006                 }
1007         }
1008 }
1009
1010 void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
1011 {
1012         int i,k,nelem;
1013
1014         cmap_grid->ngrid        = ngrid;
1015         cmap_grid->grid_spacing = grid_spacing;
1016         nelem                   = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
1017
1018         snew(cmap_grid->cmapdata,ngrid);
1019
1020         for(i=0;i<cmap_grid->ngrid;i++)
1021         {
1022                 snew(cmap_grid->cmapdata[i].cmap,4*nelem);
1023         }
1024 }
1025
1026
1027 static int count_constraints(gmx_mtop_t *mtop,t_molinfo *mi,warninp_t wi)
1028 {
1029   int count,count_mol,i,mb;
1030   gmx_molblock_t *molb;
1031   t_params *plist;
1032   char buf[STRLEN];
1033
1034   count = 0;
1035   for(mb=0; mb<mtop->nmolblock; mb++) {
1036     count_mol = 0;
1037     molb  = &mtop->molblock[mb];
1038     plist = mi[molb->type].plist;
1039
1040     for(i=0; i<F_NRE; i++) {
1041       if (i == F_SETTLE)
1042         count_mol += 3*plist[i].nr;
1043       else if (interaction_function[i].flags & IF_CONSTRAINT)
1044         count_mol += plist[i].nr;
1045     }
1046
1047     if (count_mol > nrdf_internal(&mi[molb->type].atoms)) {
1048       sprintf(buf,
1049               "Molecule type '%s' has %d constraints.\n"
1050               "For stability and efficiency there should not be more constraints than internal number of degrees of freedom: %d.\n",
1051               *mi[molb->type].name,count_mol,
1052               nrdf_internal(&mi[molb->type].atoms));
1053       warning(wi,buf);
1054     }
1055     count += molb->nmol*count_mol;
1056   }
1057
1058   return count;
1059 }
1060
1061 static void check_gbsa_params_charged(gmx_mtop_t *sys, gpp_atomtype_t atype)
1062 {
1063     int i,nmiss,natoms,mt;
1064     real q;
1065     const t_atoms *atoms;
1066
1067     nmiss = 0;
1068     for(mt=0;mt<sys->nmoltype;mt++)
1069     {
1070         atoms  = &sys->moltype[mt].atoms;
1071         natoms = atoms->nr;
1072
1073         for(i=0;i<natoms;i++)
1074         {
1075             q = atoms->atom[i].q;
1076             if ((get_atomtype_radius(atoms->atom[i].type,atype)    == 0  ||
1077                  get_atomtype_vol(atoms->atom[i].type,atype)       == 0  ||
1078                  get_atomtype_surftens(atoms->atom[i].type,atype)  == 0  ||
1079                  get_atomtype_gb_radius(atoms->atom[i].type,atype) == 0  ||
1080                  get_atomtype_S_hct(atoms->atom[i].type,atype)     == 0) &&
1081                 q != 0)
1082             {
1083                 fprintf(stderr,"\nGB parameter(s) zero for atom type '%s' while charge is %g\n",
1084                         get_atomtype_name(atoms->atom[i].type,atype),q);
1085                 nmiss++;
1086             }
1087         }
1088     }
1089
1090     if (nmiss > 0)
1091     {
1092         gmx_fatal(FARGS,"Can't do GB electrostatics; the implicit_genborn_params section of the forcefield has parameters with value zero for %d atomtypes that occur as charged atoms.",nmiss);
1093     }
1094 }
1095
1096
1097 static void check_gbsa_params(t_inputrec *ir,gpp_atomtype_t atype)
1098 {
1099     int  nmiss,i;
1100
1101     /* If we are doing GBSA, check that we got the parameters we need
1102      * This checking is to see if there are GBSA paratmeters for all
1103      * atoms in the force field. To go around this for testing purposes
1104      * comment out the nerror++ counter temporarily
1105      */
1106     nmiss = 0;
1107     for(i=0;i<get_atomtype_ntypes(atype);i++)
1108     {
1109         if (get_atomtype_radius(i,atype)    < 0 ||
1110             get_atomtype_vol(i,atype)       < 0 ||
1111             get_atomtype_surftens(i,atype)  < 0 ||
1112             get_atomtype_gb_radius(i,atype) < 0 ||
1113             get_atomtype_S_hct(i,atype)     < 0)
1114         {
1115             fprintf(stderr,"\nGB parameter(s) missing or negative for atom type '%s'\n",
1116                     get_atomtype_name(i,atype));
1117             nmiss++;
1118         }
1119     }
1120
1121     if (nmiss > 0)
1122     {
1123         gmx_fatal(FARGS,"Can't do GB electrostatics; the implicit_genborn_params section of the forcefield is missing parameters for %d atomtypes or they might be negative.",nmiss);
1124     }
1125
1126 }
1127
1128 static void check_settle(gmx_mtop_t   *sys)
1129 {
1130     int i,j,cgj1,nra;
1131
1132     nra = interaction_function[F_SETTLE].nratoms;
1133     for(i=0; (i<sys->nmoltype); i++)
1134     {
1135         for(j=0; (j<sys->moltype[i].ilist[F_SETTLE].nr); j+=nra+1)
1136         {
1137             cgj1 = sys->moltype[i].cgs.index[j+1];
1138             if (j+2 >= cgj1)
1139                 gmx_fatal(FARGS,"For SETTLE you need to have all atoms involved in one charge group. Please fix your topology.");
1140         }
1141     }
1142 }
1143
1144 int main (int argc, char *argv[])
1145 {
1146   static const char *desc[] = {
1147     "The gromacs preprocessor",
1148     "reads a molecular topology file, checks the validity of the",
1149     "file, expands the topology from a molecular description to an atomic",
1150     "description. The topology file contains information about",
1151     "molecule types and the number of molecules, the preprocessor",
1152     "copies each molecule as needed. ",
1153     "There is no limitation on the number of molecule types. ",
1154     "Bonds and bond-angles can be converted into constraints, separately",
1155     "for hydrogens and heavy atoms.",
1156     "Then a coordinate file is read and velocities can be generated",
1157     "from a Maxwellian distribution if requested.",
1158     "[TT]grompp[tt] also reads parameters for the [TT]mdrun[tt] ",
1159     "(eg. number of MD steps, time step, cut-off), and others such as",
1160     "NEMD parameters, which are corrected so that the net acceleration",
1161     "is zero.",
1162     "Eventually a binary file is produced that can serve as the sole input",
1163     "file for the MD program.[PAR]",
1164
1165     "[TT]grompp[tt] uses the atom names from the topology file. The atom names",
1166     "in the coordinate file (option [TT]-c[tt]) are only read to generate",
1167     "warnings when they do not match the atom names in the topology.",
1168     "Note that the atom names are irrelevant for the simulation as",
1169     "only the atom types are used for generating interaction parameters.[PAR]",
1170
1171     "[TT]grompp[tt] uses a built-in preprocessor to resolve includes, macros, ",
1172     "etc. The preprocessor supports the following keywords:[PAR]",
1173     "#ifdef VARIABLE[BR]",
1174     "#ifndef VARIABLE[BR]",
1175     "#else[BR]",
1176     "#endif[BR]",
1177     "#define VARIABLE[BR]",
1178     "#undef VARIABLE[BR]"
1179     "#include \"filename\"[BR]",
1180     "#include <filename>[PAR]",
1181     "The functioning of these statements in your topology may be modulated by",
1182     "using the following two flags in your [TT].mdp[tt] file:[PAR]",
1183     "[TT]define = -DVARIABLE1 -DVARIABLE2[BR]",
1184     "include = -I/home/john/doe[tt][BR]",
1185     "For further information a C-programming textbook may help you out.",
1186     "Specifying the [TT]-pp[tt] flag will get the pre-processed",
1187     "topology file written out so that you can verify its contents.[PAR]",
1188
1189     /* cpp has been unnecessary for some time, hasn't it?
1190         "If your system does not have a C-preprocessor, you can still",
1191         "use [TT]grompp[tt], but you do not have access to the features ",
1192         "from the cpp. Command line options to the C-preprocessor can be given",
1193         "in the [TT].mdp[tt] file. See your local manual (man cpp).[PAR]",
1194     */
1195
1196     "When using position restraints a file with restraint coordinates",
1197     "can be supplied with [TT]-r[tt], otherwise restraining will be done",
1198     "with respect to the conformation from the [TT]-c[tt] option.",
1199     "For free energy calculation the the coordinates for the B topology",
1200     "can be supplied with [TT]-rb[tt], otherwise they will be equal to",
1201     "those of the A topology.[PAR]",
1202
1203     "Starting coordinates can be read from trajectory with [TT]-t[tt].",
1204     "The last frame with coordinates and velocities will be read,",
1205     "unless the [TT]-time[tt] option is used. Only if this information",
1206     "is absent will the coordinates in the [TT]-c[tt] file be used.",
1207     "Note that these velocities will not be used when [TT]gen_vel = yes[tt]",
1208     "in your [TT].mdp[tt] file. An energy file can be supplied with",
1209     "[TT]-e[tt] to read Nose-Hoover and/or Parrinello-Rahman coupling",
1210     "variables.[PAR]",
1211
1212     "[TT]grompp[tt] can be used to restart simulations (preserving",
1213     "continuity) by supplying just a checkpoint file with [TT]-t[tt].",
1214     "However, for simply changing the number of run steps to extend",
1215     "a run, using [TT]tpbconv[tt] is more convenient than [TT]grompp[tt].",
1216     "You then supply the old checkpoint file directly to [TT]mdrun[tt]",
1217     "with [TT]-cpi[tt]. If you wish to change the ensemble or things",
1218     "like output frequency, then supplying the checkpoint file to",
1219     "[TT]grompp[tt] with [TT]-t[tt] along with a new [TT].mdp[tt] file",
1220     "with [TT]-f[tt] is the recommended procedure.[PAR]",
1221
1222     "By default, all bonded interactions which have constant energy due to",
1223     "virtual site constructions will be removed. If this constant energy is",
1224     "not zero, this will result in a shift in the total energy. All bonded",
1225     "interactions can be kept by turning off [TT]-rmvsbds[tt]. Additionally,",
1226     "all constraints for distances which will be constant anyway because",
1227     "of virtual site constructions will be removed. If any constraints remain",
1228     "which involve virtual sites, a fatal error will result.[PAR]"
1229
1230     "To verify your run input file, please take note of all warnings",
1231     "on the screen, and correct where necessary. Do also look at the contents",
1232     "of the [TT]mdout.mdp[tt] file; this contains comment lines, as well as",
1233     "the input that [TT]grompp[tt] has read. If in doubt, you can start [TT]grompp[tt]",
1234     "with the [TT]-debug[tt] option which will give you more information",
1235     "in a file called [TT]grompp.log[tt] (along with real debug info). You",
1236     "can see the contents of the run input file with the [TT]gmxdump[tt]",
1237     "program. [TT]gmxcheck[tt] can be used to compare the contents of two",
1238     "run input files.[PAR]"
1239
1240     "The [TT]-maxwarn[tt] option can be used to override warnings printed",
1241     "by [TT]grompp[tt] that otherwise halt output. In some cases, warnings are",
1242     "harmless, but usually they are not. The user is advised to carefully",
1243     "interpret the output messages before attempting to bypass them with",
1244     "this option."
1245   };
1246   t_gromppopts *opts;
1247   gmx_mtop_t   *sys;
1248   int          nmi;
1249   t_molinfo    *mi;
1250   gpp_atomtype_t atype;
1251   t_inputrec   *ir;
1252   int          natoms,nvsite,comb,mt;
1253   t_params     *plist;
1254   t_state      state;
1255   matrix       box;
1256   real         max_spacing,fudgeQQ;
1257   double       reppow;
1258   char         fn[STRLEN],fnB[STRLEN];
1259   const char   *mdparin;
1260   int          ntype;
1261   gmx_bool         bNeedVel,bGenVel;
1262   gmx_bool         have_atomnumber;
1263   int              n12,n13,n14;
1264   t_params     *gb_plist = NULL;
1265   gmx_genborn_t *born = NULL;
1266   output_env_t oenv;
1267   gmx_bool         bVerbose = FALSE;
1268   warninp_t    wi;
1269   char         warn_buf[STRLEN];
1270
1271   t_filenm fnm[] = {
1272     { efMDP, NULL,  NULL,        ffREAD  },
1273     { efMDP, "-po", "mdout",     ffWRITE },
1274     { efSTX, "-c",  NULL,        ffREAD  },
1275     { efSTX, "-r",  NULL,        ffOPTRD },
1276     { efSTX, "-rb", NULL,        ffOPTRD },
1277     { efNDX, NULL,  NULL,        ffOPTRD },
1278     { efTOP, NULL,  NULL,        ffREAD  },
1279     { efTOP, "-pp", "processed", ffOPTWR },
1280     { efTPX, "-o",  NULL,        ffWRITE },
1281     { efTRN, "-t",  NULL,        ffOPTRD },
1282     { efEDR, "-e",  NULL,        ffOPTRD },
1283     { efTRN, "-ref","rotref",    ffOPTRW }
1284   };
1285 #define NFILE asize(fnm)
1286
1287   /* Command line options */
1288   static gmx_bool bRenum=TRUE;
1289   static gmx_bool bRmVSBds=TRUE,bZero=FALSE;
1290   static int  i,maxwarn=0;
1291   static real fr_time=-1;
1292   t_pargs pa[] = {
1293     { "-v",       FALSE, etBOOL,{&bVerbose},
1294       "Be loud and noisy" },
1295     { "-time",    FALSE, etREAL, {&fr_time},
1296       "Take frame at or first after this time." },
1297     { "-rmvsbds",FALSE, etBOOL, {&bRmVSBds},
1298       "Remove constant bonded interactions with virtual sites" },
1299     { "-maxwarn", FALSE, etINT,  {&maxwarn},
1300       "Number of allowed warnings during input processing. Not for normal use and may generate unstable systems" },
1301     { "-zero",    FALSE, etBOOL, {&bZero},
1302       "Set parameters for bonded interactions without defaults to zero instead of generating an error" },
1303     { "-renum",   FALSE, etBOOL, {&bRenum},
1304       "Renumber atomtypes and minimize number of atomtypes" }
1305   };
1306
1307   CopyRight(stdout,argv[0]);
1308
1309   /* Initiate some variables */
1310   snew(ir,1);
1311   snew(opts,1);
1312   init_ir(ir,opts);
1313
1314   /* Parse the command line */
1315   parse_common_args(&argc,argv,0,NFILE,fnm,asize(pa),pa,
1316                     asize(desc),desc,0,NULL,&oenv);
1317
1318   wi = init_warning(TRUE,maxwarn);
1319
1320   /* PARAMETER file processing */
1321   mdparin = opt2fn("-f",NFILE,fnm);
1322   set_warning_line(wi,mdparin,-1);
1323   get_ir(mdparin,opt2fn("-po",NFILE,fnm),ir,opts,wi);
1324
1325   if (bVerbose)
1326     fprintf(stderr,"checking input for internal consistency...\n");
1327   check_ir(mdparin,ir,opts,wi);
1328
1329   if (ir->ld_seed == -1) {
1330     ir->ld_seed = make_seed();
1331     fprintf(stderr,"Setting the LD random seed to %d\n",ir->ld_seed);
1332   }
1333
1334   bNeedVel = EI_STATE_VELOCITY(ir->eI);
1335   bGenVel  = (bNeedVel && opts->bGenVel);
1336
1337   snew(plist,F_NRE);
1338   init_plist(plist);
1339   snew(sys,1);
1340   atype = init_atomtype();
1341   if (debug)
1342     pr_symtab(debug,0,"Just opened",&sys->symtab);
1343
1344   strcpy(fn,ftp2fn(efTOP,NFILE,fnm));
1345   if (!gmx_fexist(fn))
1346     gmx_fatal(FARGS,"%s does not exist",fn);
1347   new_status(fn,opt2fn_null("-pp",NFILE,fnm),opt2fn("-c",NFILE,fnm),
1348              opts,ir,bZero,bGenVel,bVerbose,&state,
1349              atype,sys,&nmi,&mi,plist,&comb,&reppow,&fudgeQQ,
1350              opts->bMorse,
1351              wi);
1352
1353   if (debug)
1354     pr_symtab(debug,0,"After new_status",&sys->symtab);
1355
1356   if (count_constraints(sys,mi,wi) && (ir->eConstrAlg == econtSHAKE)) {
1357     if (ir->eI == eiCG || ir->eI == eiLBFGS) {
1358         sprintf(warn_buf,"Can not do %s with %s, use %s",
1359                 EI(ir->eI),econstr_names[econtSHAKE],econstr_names[econtLINCS]);
1360         warning_error(wi,warn_buf);
1361     }
1362     if (ir->bPeriodicMols) {
1363         sprintf(warn_buf,"Can not do periodic molecules with %s, use %s",
1364                 econstr_names[econtSHAKE],econstr_names[econtLINCS]);
1365         warning_error(wi,warn_buf);
1366     }
1367   }
1368
1369   /* If we are doing QM/MM, check that we got the atom numbers */
1370   have_atomnumber = TRUE;
1371   for (i=0; i<get_atomtype_ntypes(atype); i++) {
1372     have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i,atype) >= 0);
1373   }
1374   if (!have_atomnumber && ir->bQMMM)
1375   {
1376       warning_error(wi,
1377                     "\n"
1378                     "It appears as if you are trying to run a QM/MM calculation, but the force\n"
1379                     "field you are using does not contain atom numbers fields. This is an\n"
1380                     "optional field (introduced in Gromacs 3.3) for general runs, but mandatory\n"
1381                     "for QM/MM. The good news is that it is easy to add - put the atom number as\n"
1382                     "an integer just before the mass column in ffXXXnb.itp.\n"
1383                     "NB: United atoms have the same atom numbers as normal ones.\n\n");
1384   }
1385
1386   if (ir->bAdress) {
1387     if ((ir->adress->const_wf>1) || (ir->adress->const_wf<0)) {
1388       warning_error(wi,"AdResS contant weighting function should be between 0 and 1\n\n");
1389     }
1390     /** \TODO check size of ex+hy width against box size */
1391   }
1392
1393   /* Check for errors in the input now, since they might cause problems
1394    * during processing further down.
1395    */
1396   check_warning_error(wi,FARGS);
1397
1398   if (opt2bSet("-r",NFILE,fnm))
1399     sprintf(fn,"%s",opt2fn("-r",NFILE,fnm));
1400   else
1401     sprintf(fn,"%s",opt2fn("-c",NFILE,fnm));
1402   if (opt2bSet("-rb",NFILE,fnm))
1403     sprintf(fnB,"%s",opt2fn("-rb",NFILE,fnm));
1404   else
1405     strcpy(fnB,fn);
1406
1407     if (nint_ftype(sys,mi,F_POSRES) > 0)
1408     {
1409         if (bVerbose)
1410         {
1411             fprintf(stderr,"Reading position restraint coords from %s",fn);
1412             if (strcmp(fn,fnB) == 0)
1413             {
1414                 fprintf(stderr,"\n");
1415             }
1416             else
1417             {
1418                 fprintf(stderr," and %s\n",fnB);
1419                 if (ir->efep != efepNO && ir->n_flambda > 0)
1420                 {
1421                     warning_error(wi,"Can not change the position restraint reference coordinates with lambda togther with foreign lambda calculation.");
1422                 }
1423             }
1424         }
1425         gen_posres(sys,mi,fn,fnB,
1426                    ir->refcoord_scaling,ir->ePBC,
1427                    ir->posres_com,ir->posres_comB,
1428                    wi);
1429     }
1430
1431   nvsite = 0;
1432   /* set parameters for virtual site construction (not for vsiten) */
1433   for(mt=0; mt<sys->nmoltype; mt++) {
1434     nvsite +=
1435       set_vsites(bVerbose, &sys->moltype[mt].atoms, atype, mi[mt].plist);
1436   }
1437   /* now throw away all obsolete bonds, angles and dihedrals: */
1438   /* note: constraints are ALWAYS removed */
1439   if (nvsite) {
1440     for(mt=0; mt<sys->nmoltype; mt++) {
1441       clean_vsite_bondeds(mi[mt].plist,sys->moltype[mt].atoms.nr,bRmVSBds);
1442     }
1443   }
1444
1445         /* If we are using CMAP, setup the pre-interpolation grid */
1446         if(plist->ncmap>0)
1447         {
1448                 init_cmap_grid(&sys->ffparams.cmap_grid, plist->nc, plist->grid_spacing);
1449                 setup_cmap(plist->grid_spacing, plist->nc, plist->cmap,&sys->ffparams.cmap_grid);
1450         }
1451
1452   set_wall_atomtype(atype,opts,ir);
1453   if (bRenum) {
1454     renum_atype(plist, sys, ir->wall_atomtype, atype, bVerbose);
1455     ntype = get_atomtype_ntypes(atype);
1456   }
1457
1458     if (ir->implicit_solvent != eisNO)
1459     {
1460         /* Now we have renumbered the atom types, we can check the GBSA params */
1461         check_gbsa_params(ir,atype);
1462
1463       /* Check that all atoms that have charge and/or LJ-parameters also have
1464        * sensible GB-parameters
1465        */
1466       check_gbsa_params_charged(sys,atype);
1467     }
1468
1469         /* PELA: Copy the atomtype data to the topology atomtype list */
1470         copy_atomtype_atomtypes(atype,&(sys->atomtypes));
1471
1472         if (debug)
1473     pr_symtab(debug,0,"After renum_atype",&sys->symtab);
1474
1475   if (bVerbose)
1476     fprintf(stderr,"converting bonded parameters...\n");
1477
1478   ntype = get_atomtype_ntypes(atype);
1479   convert_params(ntype, plist, mi, comb, reppow, fudgeQQ, sys);
1480
1481   if (debug)
1482     pr_symtab(debug,0,"After convert_params",&sys->symtab);
1483
1484   /* set ptype to VSite for virtual sites */
1485   for(mt=0; mt<sys->nmoltype; mt++) {
1486     set_vsites_ptype(FALSE,&sys->moltype[mt]);
1487   }
1488   if (debug) {
1489     pr_symtab(debug,0,"After virtual sites",&sys->symtab);
1490   }
1491   /* Check velocity for virtual sites and shells */
1492   if (bGenVel) {
1493     check_vel(sys,state.v);
1494   }
1495
1496   /* check for charge groups in settles */
1497   check_settle(sys);
1498
1499   /* check masses */
1500   check_mol(sys,wi);
1501
1502   for(i=0; i<sys->nmoltype; i++) {
1503       check_cg_sizes(ftp2fn(efTOP,NFILE,fnm),&sys->moltype[i].cgs,wi);
1504   }
1505
1506   if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)
1507   {
1508       check_bonds_timestep(sys,ir->delta_t,wi);
1509   }
1510
1511   if (EI_ENERGY_MINIMIZATION(ir->eI) && 0 == ir->nsteps)
1512   {
1513       warning_note(wi,"Zero-step energy minimization will alter the coordinates before calculating the energy. If you just want the energy of a single point, try zero-step MD (with unconstrained_start = yes). To do multiple single-point energy evaluations of different configurations of the same topology, use mdrun -rerun.");
1514   }
1515
1516   check_warning_error(wi,FARGS);
1517
1518   if (bVerbose)
1519     fprintf(stderr,"initialising group options...\n");
1520   do_index(mdparin,ftp2fn_null(efNDX,NFILE,fnm),
1521            sys,bVerbose,ir,
1522            bGenVel ? state.v : NULL,
1523            wi);
1524
1525   /* Init the temperature coupling state */
1526   init_gtc_state(&state,ir->opts.ngtc,0,ir->opts.nhchainlength);
1527
1528   if (bVerbose)
1529     fprintf(stderr,"Checking consistency between energy and charge groups...\n");
1530   check_eg_vs_cg(sys);
1531
1532   if (debug)
1533     pr_symtab(debug,0,"After index",&sys->symtab);
1534   triple_check(mdparin,ir,sys,wi);
1535   close_symtab(&sys->symtab);
1536   if (debug)
1537     pr_symtab(debug,0,"After close",&sys->symtab);
1538
1539   /* make exclusions between QM atoms */
1540   if (ir->bQMMM) {
1541     generate_qmexcl(sys,ir);
1542   }
1543
1544   if (ftp2bSet(efTRN,NFILE,fnm)) {
1545     if (bVerbose)
1546       fprintf(stderr,"getting data from old trajectory ...\n");
1547     cont_status(ftp2fn(efTRN,NFILE,fnm),ftp2fn_null(efEDR,NFILE,fnm),
1548                 bNeedVel,bGenVel,fr_time,ir,&state,sys,oenv);
1549   }
1550
1551     if (ir->ePBC==epbcXY && ir->nwall!=2)
1552     {
1553         clear_rvec(state.box[ZZ]);
1554     }
1555
1556     if (ir->rlist > 0)
1557     {
1558         set_warning_line(wi,mdparin,-1);
1559         check_chargegroup_radii(sys,ir,state.x,wi);
1560     }
1561
1562   if (EEL_FULL(ir->coulombtype)) {
1563     /* Calculate the optimal grid dimensions */
1564     copy_mat(state.box,box);
1565     if (ir->ePBC==epbcXY && ir->nwall==2)
1566       svmul(ir->wall_ewald_zfac,box[ZZ],box[ZZ]);
1567     max_spacing = calc_grid(stdout,box,opts->fourierspacing,
1568                             &(ir->nkx),&(ir->nky),&(ir->nkz));
1569     if ((ir->coulombtype == eelPPPM) && (max_spacing > 0.1)) {
1570         set_warning_line(wi,mdparin,-1);
1571         warning_note(wi,"Grid spacing larger then 0.1 while using PPPM.");
1572     }
1573   }
1574
1575   if (ir->ePull != epullNO)
1576     set_pull_init(ir,sys,state.x,state.box,oenv,opts->pull_start);
1577
1578   if (ir->bRot)
1579   {
1580       set_reference_positions(ir->rot,sys,state.x,state.box,
1581                               opt2fn("-ref",NFILE,fnm),opt2bSet("-ref",NFILE,fnm),
1582                               wi);
1583   }
1584
1585   /*  reset_multinr(sys); */
1586
1587   if (EEL_PME(ir->coulombtype)) {
1588     float ratio = pme_load_estimate(sys,ir,state.box);
1589     fprintf(stderr,"Estimate for the relative computational load of the PME mesh part: %.2f\n",ratio);
1590     /* With free energy we might need to do PME both for the A and B state
1591      * charges. This will double the cost, but the optimal performance will
1592      * then probably be at a slightly larger cut-off and grid spacing.
1593      */
1594     if ((ir->efep == efepNO && ratio > 1.0/2.0) ||
1595         (ir->efep != efepNO && ratio > 2.0/3.0)) {
1596         warning_note(wi,
1597                      "The optimal PME mesh load for parallel simulations is below 0.5\n"
1598                    "and for highly parallel simulations between 0.25 and 0.33,\n"
1599                    "for higher performance, increase the cut-off and the PME grid spacing");
1600     }
1601   }
1602
1603     {
1604         char warn_buf[STRLEN];
1605         double cio = compute_io(ir,sys->natoms,&sys->groups,F_NRE,1);
1606         sprintf(warn_buf,"This run will generate roughly %.0f Mb of data",cio);
1607         if (cio > 2000) {
1608             set_warning_line(wi,mdparin,-1);
1609             warning_note(wi,warn_buf);
1610         } else {
1611             printf("%s\n",warn_buf);
1612         }
1613     }
1614
1615   if (bVerbose)
1616     fprintf(stderr,"writing run input file...\n");
1617
1618   done_warning(wi,FARGS);
1619
1620   state.lambda = ir->init_lambda;
1621   write_tpx_state(ftp2fn(efTPX,NFILE,fnm),ir,&state,sys);
1622
1623   thanx(stderr);
1624
1625   return 0;
1626 }