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