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