src/gromacs/gmxpreprocess/topio.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.0
  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 <math.h>
  41 #include <sys/types.h>
  42 #include <stdio.h>
  43 #include <string.h>
  44 #include <errno.h>
  45 #include <ctype.h>
  46 #include "futil.h"
  47 #include "sysstuff.h"
  48 #include "typedefs.h"
  49 #include "smalloc.h"
  50 #include "macros.h"
  51 #include "gmxfio.h"
  52 #include "txtdump.h"
  53 #include "physics.h"
  54 #include "macros.h"
  55 #include "names.h"
  56 #include "string2.h"
  57 #include "symtab.h"
  58 #include "gmx_fatal.h"
  59 #include "warninp.h"
  60 #include "vsite_parm.h"
  61
  62 #include "grompp.h"
  63 #include "toputil.h"
  64 #include "toppush.h"
  65 #include "topdirs.h"
  66 #include "gpp_nextnb.h"
  67 #include "topio.h"
  68 #include "topshake.h"
  69 #include "gmxcpp.h"
  70 #include "gpp_bond_atomtype.h"
  71 #include "genborn.h"
  72 #include "maths.h"
  73
  74 #define CPPMARK         '#'     /* mark from cpp                        */
  75 #define OPENDIR         '['     /* starting sign for directive          */
  76 #define CLOSEDIR        ']'     /* ending sign for directive            */
  77
  78 #define PM()  \
  79   printf("line: %d, maxavail: %d\n",__LINE__,maxavail()); \
  80   fflush(stdout)
  81
  82 static void free_nbparam(t_nbparam **param,int nr)
  83 {
  84   int i;
  85
  86   for(i=0; i<nr; i++)
  87     sfree(param[i]);
  88   sfree(param);
  89 }
  90
  91 static int copy_nbparams(t_nbparam **param,int ftype,t_params *plist,int nr)
  92 {
  93   int i,j,f;
  94   int nrfp,ncopy;
  95
  96   nrfp = NRFP(ftype);
  97
  98   ncopy = 0;
  99   for(i=0; i<nr; i++)
 100     for(j=0; j<=i; j++)
 101       if (param[i][j].bSet) {
 102         for(f=0; f<nrfp; f++) {
 103           plist->param[nr*i+j].c[f] = param[i][j].c[f];
 104           plist->param[nr*j+i].c[f] = param[i][j].c[f];
 105         }
 106         ncopy++;
 107       }
 108
 109   return ncopy;
 110 }
 111
 112 static void gen_pairs(t_params *nbs,t_params *pairs,real fudge, int comb, gmx_bool bVerbose)
 113 {
 114     int     i,j,ntp,nrfp,nrfpA,nrfpB,nnn;
 115     real    scaling;
 116     ntp       = nbs->nr;
 117     nnn       = sqrt(ntp);
 118     nrfp      = NRFP(F_LJ);
 119     nrfpA     = interaction_function[F_LJ14].nrfpA;
 120     nrfpB     = interaction_function[F_LJ14].nrfpB;
 121     pairs->nr = ntp;
 122
 123   if ((nrfp  != nrfpA) || (nrfpA != nrfpB))
 124     gmx_incons("Number of force parameters in gen_pairs wrong");
 125
 126   fprintf(stderr,"Generating 1-4 interactions: fudge = %g\n",fudge);
 127   if (debug) {
 128     fprintf(debug,"Fudge factor for 1-4 interactions: %g\n",fudge);
 129     fprintf(debug,"Holy Cow! there are %d types\n",ntp);
 130   }
 131   snew(pairs->param,pairs->nr);
 132   for(i=0; (i<ntp); i++) {
 133     /* Copy param.a */
 134     pairs->param[i].a[0] = i / nnn;
 135     pairs->param[i].a[1] = i % nnn;
 136     /* Copy normal and FEP parameters and multiply by fudge factor */
 137
 138
 139
 140     for(j=0; (j<nrfp); j++) {
 141         /* If we are using sigma/epsilon values, only the epsilon values
 142          * should be scaled, but not sigma.
 143          * The sigma values have even indices 0,2, etc.
 144          */
 145         if ((comb == eCOMB_ARITHMETIC || comb == eCOMB_GEOM_SIG_EPS) && (j%2==0))
 146             scaling = 1.0;
 147         else
 148             scaling = fudge;
 149
 150         pairs->param[i].c[j]=scaling*nbs->param[i].c[j];
 151         pairs->param[i].c[nrfp+j]=scaling*nbs->param[i].c[j];
 152     }
 153   }
 154 }
 155
 156 double check_mol(gmx_mtop_t *mtop,warninp_t wi)
 157 {
 158   char    buf[256];
 159   int     i,mb,nmol,ri,pt;
 160   double  q;
 161   real    m;
 162   t_atoms *atoms;
 163
 164   /* Check mass and charge */
 165   q=0.0;
 166
 167   for(mb=0; mb<mtop->nmoltype; mb++) {
 168     atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
 169     nmol  = mtop->molblock[mb].nmol;
 170     for (i=0; (i<atoms->nr); i++) {
 171       q += nmol*atoms->atom[i].q;
 172       m  = atoms->atom[i].m;
 173       pt = atoms->atom[i].ptype;
 174       /* If the particle is an atom or a nucleus it must have a mass,
 175        * else, if it is a shell, a vsite or a bondshell it can have mass zero
 176        */
 177       if ((m <= 0.0) && ((pt == eptAtom) || (pt == eptNucleus))) {
 178         ri = atoms->atom[i].resind;
 179         sprintf(buf,"atom %s (Res %s-%d) has mass %g\n",
 180                 *(atoms->atomname[i]),
 181                 *(atoms->resinfo[ri].name),
 182                 atoms->resinfo[ri].nr,
 183                 m);
 184         warning_error(wi,buf);
 185       } else
 186         if ((m!=0) && (pt == eptVSite)) {
 187           ri = atoms->atom[i].resind;
 188           sprintf(buf,"virtual site %s (Res %s-%d) has non-zero mass %g\n"
 189                   "     Check your topology.\n",
 190                   *(atoms->atomname[i]),
 191                   *(atoms->resinfo[ri].name),
 192                   atoms->resinfo[ri].nr,
 193                   m);
 194           warning_error(wi,buf);
 195           /* The following statements make LINCS break! */
 196           /* atoms->atom[i].m=0; */
 197         }
 198     }
 199   }
 200   return q;
 201 }
 202
 203 static void sum_q(t_atoms *atoms,int n,double *qt,double *qBt)
 204 {
 205     double qmolA,qmolB;
 206     int     i;
 207
 208     /* sum charge */
 209     qmolA = 0;
 210     qmolB = 0;
 211     for (i=0; i<atoms->nr; i++)
 212     {
 213         qmolA += atoms->atom[i].q;
 214         qmolB += atoms->atom[i].qB;
 215     }
 216     /* Unfortunately an absolute comparison,
 217      * but this avoids unnecessary warnings and gmx-users mails.
 218      */
 219     if (fabs(qmolA) >= 1e-6 || fabs(qmolB) >= 1e-6)
 220     {
 221         *qt  += n*qmolA;
 222         *qBt += n*qmolB;
 223     }
 224 }
 225
 226 static void get_nbparm(char *nb_str,char *comb_str,int *nb,int *comb,
 227                        warninp_t wi)
 228 {
 229   int  i;
 230   char warn_buf[STRLEN];
 231
 232   *nb   = -1;
 233   for(i=1; (i<eNBF_NR); i++)
 234     if (gmx_strcasecmp(nb_str,enbf_names[i]) == 0)
 235       *nb = i;
 236   if (*nb == -1)
 237     *nb = strtol(nb_str,NULL,10);
 238   if ((*nb < 1) || (*nb >= eNBF_NR)) {
 239     sprintf(warn_buf,"Invalid nonbond function selector '%s' using %s",
 240             nb_str,enbf_names[1]);
 241     warning_error(wi,warn_buf);
 242     *nb = 1;
 243   }
 244   *comb = -1;
 245   for(i=1; (i<eCOMB_NR); i++)
 246     if (gmx_strcasecmp(comb_str,ecomb_names[i]) == 0)
 247       *comb = i;
 248   if (*comb == -1)
 249     *comb = strtol(comb_str,NULL,10);
 250   if ((*comb < 1) || (*comb >= eCOMB_NR)) {
 251     sprintf(warn_buf,"Invalid combination rule selector '%s' using %s",
 252             comb_str,ecomb_names[1]);
 253     warning_error(wi,warn_buf);
 254     *comb = 1;
 255   }
 256 }
 257
 258 static char ** cpp_opts(const char *define,const char *include,
 259                         const char *infile,
 260                         warninp_t wi)
 261 {
 262   int  n,len;
 263   int  ncppopts=0;
 264   const char *cppadds[2];
 265   char **cppopts  = NULL;
 266   const char *option[2] = { "-D","-I" };
 267   const char *nopt[2]   = { "define", "include" };
 268   const char *ptr;
 269   const char *rptr;
 270   char *buf;
 271   char warn_buf[STRLEN];
 272
 273   cppadds[0] = define;
 274   cppadds[1] = include;
 275   for(n=0; (n<2); n++) {
 276     if (cppadds[n]) {
 277       ptr = cppadds[n];
 278       while(*ptr != '\0') {
 279         while((*ptr != '\0') && isspace(*ptr))
 280           ptr++;
 281         rptr = ptr;
 282         while((*rptr != '\0') && !isspace(*rptr))
 283           rptr++;
 284         len = (rptr - ptr);
 285         if (len > 2) {
 286           snew(buf,(len+1));
 287           strncpy(buf,ptr,len);
 288           if (strstr(ptr,option[n]) != ptr) {
 289           set_warning_line(wi,"mdp file",-1);
 290           sprintf(warn_buf,"Malformed %s option %s",nopt[n],buf);
 291           warning(wi,warn_buf);
 292           }
 293           else {
 294             srenew(cppopts,++ncppopts);
 295             cppopts[ncppopts-1] = strdup(buf);
 296           }
 297           sfree(buf);
 298           ptr = rptr;
 299         }
 300       }
 301     }
 302   }
 303   srenew(cppopts,++ncppopts);
 304   cppopts[ncppopts-1] = NULL;
 305
 306   return cppopts;
 307 }
 308
 309
 310 int
 311 find_gb_bondlength(t_params *plist,int ai,int aj, real *length)
 312 {
 313     int i,j,a1,a2;
 314
 315     int found=0;
 316     int status;
 317
 318     for(i=0;i<F_NRE && !found;i++)
 319     {
 320         if(IS_CHEMBOND(i))
 321         {
 322             for(j=0;j<plist[i].nr; j++)
 323             {
 324                 a1=plist[i].param[j].a[0];
 325                 a2=plist[i].param[j].a[1];
 326
 327                 if( (a1==ai && a2==aj) || (a1==aj && a2==ai))
 328                 {
 329                     /* Equilibrium bond distance */
 330                     *length=plist[i].param[j].c[0];
 331                     found=1;
 332                 }
 333             }
 334         }
 335     }
 336     status = !found;
 337
 338     return status;
 339 }
 340
 341
 342 int
 343 find_gb_anglelength(t_params *plist,int ai,int ak, real *length)
 344 {
 345     int i,j,a1,a2,a3;
 346     real r12,r23,a123;
 347     int found=0;
 348     int status,status1,status2;
 349
 350     r12 = r23 = 0;
 351
 352     for(i=0;i<F_NRE && !found;i++)
 353     {
 354         if(IS_ANGLE(i))
 355         {
 356             for(j=0;j<plist[i].nr; j++)
 357             {
 358                 a1=plist[i].param[j].a[0];
 359                 a2=plist[i].param[j].a[1];
 360                 a3=plist[i].param[j].a[2];
 361
 362                 /* We dont care what the middle atom is, but use it below */
 363                 if( (a1==ai && a3==ak) || (a1==ak && a3==ai) )
 364                 {
 365                     /* Equilibrium bond distance */
 366                     a123 = plist[i].param[j].c[0];
 367                     /* Use middle atom to find reference distances r12 and r23 */
 368                     status1 = find_gb_bondlength(plist,a1,a2,&r12);
 369                     status2 = find_gb_bondlength(plist,a2,a3,&r23);
 370
 371                     if(status1==0 && status2==0)
 372                     {
 373                         /* cosine theorem to get r13 */
 374                         *length=sqrt(r12*r12+r23*r23-(2*r12*r23*cos(a123/RAD2DEG)));
 375                         found=1;
 376                     }
 377                 }
 378             }
 379         }
 380     }
 381     status = !found;
 382
 383     return status;
 384 }
 385
 386 int
 387 generate_gb_exclusion_interactions(t_molinfo *mi,gpp_atomtype_t atype,t_nextnb *nnb)
 388 {
 389         int          i,j,k,n,ai,aj,ti,tj;
 390     int          n12,n13,n14;
 391     int          ftype;
 392     t_param      param;
 393     t_params *   plist;
 394     t_atoms *    at;
 395     real         radiusi,radiusj;
 396     real         gb_radiusi,gb_radiusj;
 397     real         param_c2,param_c4;
 398     real         distance;
 399
 400     plist = mi->plist;
 401     at    = &mi->atoms;
 402
 403     for(n=1;n<=nnb->nrex;n++)
 404     {
 405         switch(n)
 406         {
 407             case 1:
 408                 ftype=F_GB12;
 409                 param_c2 = STILL_P2;
 410                 param_c4 = 0.8875;
 411                 break;
 412             case 2:
 413                 ftype=F_GB13;
 414                 param_c2 = STILL_P3;
 415                 param_c4 = 0.3516;
 416                 break;
 417             default:
 418                 /* Put all higher-order exclusions into 1,4 list so we dont miss them */
 419                 ftype=F_GB14;
 420                 param_c2 = STILL_P3;
 421                 param_c4 = 0.3516;
 422                 break;
 423         }
 424
 425         for(ai=0;ai<nnb->nr;ai++)
 426         {
 427             ti         = at->atom[ai].type;
 428             radiusi    = get_atomtype_radius(ti,atype);
 429             gb_radiusi = get_atomtype_gb_radius(ti,atype);
 430
 431             for(j=0;j<nnb->nrexcl[ai][n];j++)
 432             {
 433                 aj = nnb->a[ai][n][j];
 434
 435                 /* Only add the interactions once */
 436                 if(aj>ai)
 437                 {
 438                     tj = at->atom[aj].type;
 439                     radiusj    = get_atomtype_radius(tj,atype);
 440                     gb_radiusj = get_atomtype_gb_radius(tj,atype);
 441
 442                    /* There is an exclusion of type "ftype" between atoms ai and aj */
 443                     param.a[0] = ai;
 444                     param.a[1] = aj;
 445
 446                     /* Reference distance, not used for 1-4 interactions */
 447                     switch(ftype)
 448                     {
 449                         case F_GB12:
 450                             if(find_gb_bondlength(plist,ai,aj,&distance)!=0)
 451                             {
 452                                 gmx_fatal(FARGS,"Cannot find bond length for atoms %d-%d",ai,aj);
 453                             }
 454                             break;
 455                         case F_GB13:
 456                             if(find_gb_anglelength(plist,ai,aj,&distance)!=0)
 457                             {
 458                                 gmx_fatal(FARGS,"Cannot find length for atoms %d-%d involved in angle",ai,aj);
 459                             }
 460                             break;
 461                         default:
 462                             distance=-1;
 463                             break;
 464                     }
 465                     /* Assign GB parameters */
 466                     /* Sum of radii */
 467                     param.c[0] = radiusi+radiusj;
 468                     /* Reference distance distance */
 469                     param.c[1] = distance;
 470                     /* Still parameter */
 471                     param.c[2] = param_c2;
 472                     /* GB radius */
 473                     param.c[3] = gb_radiusi+gb_radiusj;
 474                     /* Parameter */
 475                     param.c[4] = param_c4;
 476
 477                     /* Add it to the parameter list */
 478                     add_param_to_list(&plist[ftype],&param);
 479                 }
 480             }
 481         }
 482         }
 483         return 0;
 484 }
 485
 486
 487
 488
 489 static char **read_topol(const char *infile,const char *outfile,
 490                          const char *define,const char *include,
 491                          t_symtab    *symtab,
 492                          gpp_atomtype_t atype,
 493                          int         *nrmols,
 494                          t_molinfo   **molinfo,
 495                          t_params    plist[],
 496                          int         *combination_rule,
 497                          double      *reppow,
 498                          t_gromppopts *opts,
 499                          real        *fudgeQQ,
 500                          int         *nmolblock,
 501                          gmx_molblock_t **molblock,
 502                          gmx_bool        bFEP,
 503                          gmx_bool        bGenborn,
 504                          gmx_bool        bZero,
 505                          gmx_bool        bVerbose,
 506                          warninp_t   wi)
 507 {
 508   FILE       *out;
 509   int        i,sl,nb_funct,comb;
 510   char       *pline=NULL,**title=NULL;
 511   char       line[STRLEN],errbuf[256],comb_str[256],nb_str[256];
 512   char       genpairs[32];
 513   char       *dirstr,*dummy2;
 514   int        nrcopies,nmol,nmolb=0,nscan,ncombs,ncopy;
 515   double     fLJ,fQQ,fPOW;
 516   gmx_molblock_t *molb=NULL;
 517   t_topology *block=NULL;
 518   t_molinfo  *mi0=NULL;
 519   DirStack   *DS;
 520   directive  d,newd;
 521   t_nbparam  **nbparam,**pair;
 522   t_block2   *block2;
 523   real       fudgeLJ=-1;    /* Multiplication factor to generate 1-4 from LJ */
 524   gmx_bool       bReadDefaults,bReadMolType,bGenPairs,bWarn_copy_A_B;
 525   double     qt=0,qBt=0; /* total charge */
 526   t_bond_atomtype batype;
 527   int        lastcg=-1;
 528   int        dcatt=-1,nmol_couple;
 529   /* File handling variables */
 530   int        status,done;
 531   gmx_cpp_t  handle;
 532   char     *tmp_line=NULL;
 533   char       warn_buf[STRLEN];
 534   const char *floating_point_arithmetic_tip =
 535       "Total charge should normally be an integer. See\n"
 536       "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
 537       "for discussion on how close it should be to an integer.\n";
 538   /* We need to open the output file before opening the input file,
 539    * because cpp_open_file can change the current working directory.
 540    */
 541   if (outfile) {
 542     out = gmx_fio_fopen(outfile,"w");
 543   } else {
 544     out = NULL;
 545   }
 546
 547   /* open input file */
 548   status = cpp_open_file(infile,&handle,cpp_opts(define,include,infile,wi));
 549   if (status != 0)
 550     gmx_fatal(FARGS,cpp_error(&handle,status));
 551
 552   /* some local variables */
 553   DS_Init(&DS);                 /* directive stack                       */
 554   nmol     = 0;                 /* no molecules yet...                   */
 555   d        = d_invalid;         /* first thing should be a directive     */
 556   nbparam  = NULL;              /* The temporary non-bonded matrix       */
 557   pair     = NULL;              /* The temporary pair interaction matrix */
 558   block2   = NULL;              /* the extra exclusions                  */
 559   nb_funct = F_LJ;
 560   *reppow  = 12.0;              /* Default value for repulsion power     */
 561
 562   comb     = 0;
 563
 564   /* Init the number of CMAP torsion angles  and grid spacing */
 565   plist->grid_spacing = 0;
 566   plist->nc           = 0;
 567
 568   bWarn_copy_A_B = bFEP;
 569
 570   batype = init_bond_atomtype();
 571   /* parse the actual file */
 572   bReadDefaults = FALSE;
 573   bGenPairs     = FALSE;
 574   bReadMolType  = FALSE;
 575   nmol_couple = 0;
 576
 577   do {
 578     status = cpp_read_line(&handle,STRLEN,line);
 579     done = (status == eCPP_EOF);
 580     if (!done) {
 581       if (status != eCPP_OK)
 582         gmx_fatal(FARGS,cpp_error(&handle,status));
 583       else if (out)
 584         fprintf(out,"%s\n",line);
 585
 586       set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
 587
 588       pline = strdup(line);
 589
 590         /* Strip trailing '\' from pline, if it exists */
 591         sl = strlen(pline);
 592         if ((sl > 0) && (pline[sl-1] == CONTINUE)) {
 593                 pline[sl-1] = ' ';
 594         }
 595
 596                 /* build one long line from several fragments - necessary for CMAP */
 597                 while (continuing(line))
 598                 {
 599                         status = cpp_read_line(&handle,STRLEN,line);
 600             set_warning_line(wi,cpp_cur_file(&handle),cpp_cur_linenr(&handle));
 601
 602                         /* Since we depend on the '\' being present to continue to read, we copy line
 603                          * to a tmp string, strip the '\' from that string, and cat it to pline
 604                          */
 605                         tmp_line=strdup(line);
 606
 607                         sl = strlen(tmp_line);
 608                         if ((sl > 0) && (tmp_line[sl-1] == CONTINUE)) {
 609                                 tmp_line[sl-1] = ' ';
 610                         }
 611
 612                         done = (status == eCPP_EOF);
 613                         if (!done) {
 614                                 if (status != eCPP_OK)
 615                                         gmx_fatal(FARGS,cpp_error(&handle,status));
 616                                 else if (out)
 617                                         fprintf(out,"%s\n",line);
 618                         }
 619
 620                         srenew(pline,strlen(pline)+strlen(tmp_line)+1);
 621                         strcat(pline,tmp_line);
 622             sfree(tmp_line);
 623                 }
 624
 625       /* skip trailing and leading spaces and comment text */
 626       strip_comment (pline);
 627       trim (pline);
 628
 629       /* if there is something left... */
 630       if ((int)strlen(pline) > 0) {
 631         if (pline[0] == OPENDIR) {
 632           /* A directive on this line: copy the directive
 633            * without the brackets into dirstr, then
 634            * skip spaces and tabs on either side of directive
 635            */
 636           dirstr = strdup((pline+1));
 637           if ((dummy2 = strchr (dirstr,CLOSEDIR)) != NULL)
 638             (*dummy2) = 0;
 639           trim (dirstr);
 640
 641           if ((newd = str2dir(dirstr)) == d_invalid) {
 642             sprintf(errbuf,"Invalid directive %s",dirstr);
 643             warning_error(wi,errbuf);
 644           }
 645           else {
 646             /* Directive found */
 647             if (debug)
 648               fprintf(debug,"found directive '%s'\n",dir2str(newd));
 649             if (DS_Check_Order (DS,newd)) {
 650               DS_Push (&DS,newd);
 651                 d = newd;
 652             }
 653             else {
 654               /* we should print here which directives should have
 655                  been present, and which actually are */
 656               gmx_fatal(FARGS,"%s\nInvalid order for directive %s",
 657                         cpp_error(&handle,eCPP_SYNTAX),dir2str(newd));
 658               /* d = d_invalid; */
 659             }
 660           }
 661           sfree(dirstr);
 662         }
 663         else if (d != d_invalid) {
 664           /* Not a directive, just a plain string
 665            * use a gigantic switch to decode,
 666            * if there is a valid directive!
 667            */
 668           switch (d) {
 669           case d_defaults:
 670             if (bReadDefaults)
 671               gmx_fatal(FARGS,"%s\nFound a second defaults directive.\n",
 672                         cpp_error(&handle,eCPP_SYNTAX));
 673             bReadDefaults = TRUE;
 674             nscan = sscanf(pline,"%s%s%s%lf%lf%lf",
 675                            nb_str,comb_str,genpairs,&fLJ,&fQQ,&fPOW);
 676             if (nscan < 2)
 677             too_few(wi);
 678             else {
 679               bGenPairs = FALSE;
 680               fudgeLJ   = 1.0;
 681               *fudgeQQ  = 1.0;
 682
 683               get_nbparm(nb_str,comb_str,&nb_funct,&comb,wi);
 684               *combination_rule = comb;
 685               if (nscan >= 3) {
 686                 bGenPairs = (gmx_strncasecmp(genpairs,"Y",1) == 0);
 687                 if (nb_funct != eNBF_LJ && bGenPairs) {
 688                   gmx_fatal(FARGS,"Generating pair parameters is only supported with LJ non-bonded interactions");
 689                 }
 690               }
 691               if (nscan >= 4)
 692                 fudgeLJ   = fLJ;
 693               if (nscan >= 5)
 694                 *fudgeQQ  = fQQ;
 695               if (nscan >= 6)
 696                 *reppow   = fPOW;
 697             }
 698             nb_funct = ifunc_index(d_nonbond_params,nb_funct);
 699
 700             break;
 701           case d_atomtypes:
 702             push_at(symtab,atype,batype,pline,nb_funct,
 703                 &nbparam,bGenPairs ? &pair : NULL,wi);
 704             break;
 705
 706           case d_bondtypes:
 707           push_bt(d,plist,2,NULL,batype,pline,wi);
 708             break;
 709           case d_constrainttypes:
 710           push_bt(d,plist,2,NULL,batype,pline,wi);
 711             break;
 712           case d_pairtypes:
 713             if (bGenPairs)
 714             push_nbt(d,pair,atype,pline,F_LJ14,wi);
 715             else
 716             push_bt(d,plist,2,atype,NULL,pline,wi);
 717               break;
 718           case d_angletypes:
 719           push_bt(d,plist,3,NULL,batype,pline,wi);
 720             break;
 721           case d_dihedraltypes:
 722             /* Special routine that can read both 2 and 4 atom dihedral definitions. */
 723           push_dihedraltype(d,plist,batype,pline,wi);
 724             break;
 725
 726           case d_nonbond_params:
 727           push_nbt(d,nbparam,atype,pline,nb_funct,wi);
 728             break;
 729             /*
 730               case d_blocktype:
 731               nblock++;
 732               srenew(block,nblock);
 733               srenew(blockinfo,nblock);
 734               blk0=&(block[nblock-1]);
 735               bi0=&(blockinfo[nblock-1]);
 736               init_top(blk0);
 737               init_molinfo(bi0);
 738               push_molt(symtab,bi0,pline);
 739               break;
 740             */
 741
 742           case d_implicit_genborn_params:
 743           push_gb_params(atype,pline,wi);
 744             break;
 745
 746           case d_implicit_surface_params:
 747             gmx_fatal(FARGS,"Implicit surface directive not supported yet.");
 748             break;
 749
 750           case d_cmaptypes:
 751           push_cmaptype(d, plist, 5, atype, batype,pline,wi);
 752                 break;
 753
 754           case d_moleculetype: {
 755             if (!bReadMolType) {
 756               int ntype;
 757           if (opts->couple_moltype != NULL &&
 758               (opts->couple_lam0 == ecouplamNONE ||
 759                opts->couple_lam0 == ecouplamQ ||
 760                opts->couple_lam1 == ecouplamNONE ||
 761                opts->couple_lam1 == ecouplamQ))
 762           {
 763               dcatt = add_atomtype_decoupled(symtab,atype,
 764                                              &nbparam,bGenPairs?&pair:NULL);
 765           }
 766               ntype = get_atomtype_ntypes(atype);
 767               ncombs = (ntype*(ntype+1))/2;
 768               generate_nbparams(comb,nb_funct,&(plist[nb_funct]),atype,wi);
 769               ncopy = copy_nbparams(nbparam,nb_funct,&(plist[nb_funct]),
 770                                     ntype);
 771               fprintf(stderr,"Generated %d of the %d non-bonded parameter combinations\n",ncombs-ncopy,ncombs);
 772               free_nbparam(nbparam,ntype);
 773               if (bGenPairs) {
 774                 gen_pairs(&(plist[nb_funct]),&(plist[F_LJ14]),fudgeLJ,comb,bVerbose);
 775                 ncopy = copy_nbparams(pair,nb_funct,&(plist[F_LJ14]),
 776                                       ntype);
 777                 fprintf(stderr,"Generated %d of the %d 1-4 parameter combinations\n",ncombs-ncopy,ncombs);
 778                 free_nbparam(pair,ntype);
 779               }
 780               /* Copy GBSA parameters to atomtype array? */
 781
 782               bReadMolType = TRUE;
 783             }
 784
 785             push_molt(symtab,&nmol,molinfo,pline,wi);
 786             srenew(block2,nmol);
 787             block2[nmol-1].nr=0;
 788             mi0=&((*molinfo)[nmol-1]);
 789             break;
 790           }
 791           case d_atoms:
 792           push_atom(symtab,&(mi0->cgs),&(mi0->atoms),atype,pline,&lastcg,wi);
 793             break;
 794
 795           case d_pairs:
 796           push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,FALSE,
 797                     bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
 798             break;
 799
 800           case d_vsites2:
 801           case d_vsites3:
 802           case d_vsites4:
 803           case d_bonds:
 804           case d_angles:
 805           case d_constraints:
 806           case d_settles:
 807           case d_position_restraints:
 808           case d_angle_restraints:
 809           case d_angle_restraints_z:
 810           case d_distance_restraints:
 811           case d_orientation_restraints:
 812           case d_dihedral_restraints:
 813           case d_dihedrals:
 814           case d_polarization:
 815           case d_water_polarization:
 816           case d_thole_polarization:
 817             push_bond(d,plist,mi0->plist,&(mi0->atoms),atype,pline,TRUE,
 818                   bGenPairs,*fudgeQQ,bZero,&bWarn_copy_A_B,wi);
 819             break;
 820           case d_cmap:
 821           push_cmap(d,plist,mi0->plist,&(mi0->atoms),atype,pline,
 822                     &bWarn_copy_A_B,wi);
 823                 break;
 824
 825           case d_vsitesn:
 826           push_vsitesn(d,plist,mi0->plist,&(mi0->atoms),atype,pline,wi);
 827             break;
 828           case d_exclusions:
 829             if (!block2[nmol-1].nr)
 830               init_block2(&(block2[nmol-1]),mi0->atoms.nr);
 831             push_excl(pline,&(block2[nmol-1]));
 832             break;
 833           case d_system:
 834             trim(pline);
 835             title=put_symtab(symtab,pline);
 836             break;
 837           case d_molecules: {
 838             int  whichmol;
 839             gmx_bool bCouple;
 840
 841             push_mol(nmol,*molinfo,pline,&whichmol,&nrcopies,wi);
 842             mi0=&((*molinfo)[whichmol]);
 843             srenew(molb,nmolb+1);
 844             molb[nmolb].type = whichmol;
 845             molb[nmolb].nmol = nrcopies;
 846             nmolb++;
 847
 848         bCouple = (opts->couple_moltype != NULL &&
 849                    (gmx_strcasecmp("system"    ,opts->couple_moltype) == 0 ||
 850                     gmx_strcasecmp(*(mi0->name),opts->couple_moltype) == 0));
 851         if (bCouple) {
 852             nmol_couple += nrcopies;
 853         }
 854
 855             if (mi0->atoms.nr == 0) {
 856               gmx_fatal(FARGS,"Molecule type '%s' contains no atoms",
 857                         *mi0->name);
 858             }
 859             fprintf(stderr,
 860                     "Excluding %d bonded neighbours molecule type '%s'\n",
 861                     mi0->nrexcl,*mi0->name);
 862             sum_q(&mi0->atoms,nrcopies,&qt,&qBt);
 863             if (!mi0->bProcessed)
 864               {
 865             t_nextnb nnb;
 866             generate_excl(mi0->nrexcl,
 867                           mi0->atoms.nr,
 868                           mi0->plist,
 869                           &nnb,
 870                           &(mi0->excls));
 871             merge_excl(&(mi0->excls),&(block2[whichmol]));
 872             done_block2(&(block2[whichmol]));
 873             make_shake(mi0->plist,&mi0->atoms,atype,opts->nshake);
 874
 875
 876
 877             /* nnb contains information about first,2nd,3rd bonded neighbors.
 878              * Use this to generate GB 1-2,1-3,1-4 interactions when necessary.
 879              */
 880             if (bGenborn==TRUE)
 881             {
 882                 generate_gb_exclusion_interactions(mi0,atype,&nnb);
 883             }
 884
 885             done_nnb(&nnb);
 886
 887             if (bCouple)
 888             {
 889                 convert_moltype_couple(mi0,dcatt,*fudgeQQ,
 890                                        opts->couple_lam0,opts->couple_lam1,
 891                                        opts->bCoupleIntra,
 892                                        nb_funct,&(plist[nb_funct]));
 893             }
 894             stupid_fill_block(&mi0->mols,mi0->atoms.nr,TRUE);
 895             mi0->bProcessed=TRUE;
 896             }
 897             break;
 898           }
 899           default:
 900             fprintf (stderr,"case: %d\n",d);
 901             invalid_case();
 902           }
 903         }
 904       }
 905       sfree(pline);
 906       pline=NULL;
 907     }
 908   } while (!done);
 909   status = cpp_close_file(&handle);
 910   if (status != eCPP_OK)
 911     gmx_fatal(FARGS,cpp_error(&handle,status));
 912   if (out)
 913     gmx_fio_fclose(out);
 914
 915   if (opts->couple_moltype) {
 916     if (nmol_couple == 0) {
 917       gmx_fatal(FARGS,"Did not find any molecules of type '%s' for coupling",
 918                 opts->couple_moltype);
 919     }
 920     fprintf(stderr,"Coupling %d copies of molecule type '%s'\n",
 921             nmol_couple,opts->couple_moltype);
 922   }
 923
 924   /* this is not very clean, but fixes core dump on empty system name */
 925   if(!title)
 926     title=put_symtab(symtab,"");
 927   if (fabs(qt) > 1e-4) {
 928     sprintf(warn_buf,"System has non-zero total charge: %.6f\n%s\n",qt,floating_point_arithmetic_tip);
 929     warning_note(wi,warn_buf);
 930   }
 931   if (fabs(qBt) > 1e-4 && !gmx_within_tol(qBt,qt,1e-6)) {
 932     sprintf(warn_buf,"State B has non-zero total charge: %.6f\n%s\n",qBt,floating_point_arithmetic_tip);
 933     warning_note(wi,warn_buf);
 934   }
 935   DS_Done (&DS);
 936   for(i=0; i<nmol; i++)
 937     done_block2(&(block2[i]));
 938   free(block2);
 939
 940   done_bond_atomtype(&batype);
 941
 942   *nrmols=nmol;
 943
 944   *nmolblock = nmolb;
 945   *molblock  = molb;
 946
 947   return title;
 948 }
 949
 950 char **do_top(gmx_bool         bVerbose,
 951               const char   *topfile,
 952               const char   *topppfile,
 953               t_gromppopts *opts,
 954               gmx_bool         bZero,
 955               t_symtab     *symtab,
 956               t_params     plist[],
 957               int          *combination_rule,
 958               double       *repulsion_power,
 959               real         *fudgeQQ,
 960               gpp_atomtype_t atype,
 961               int          *nrmols,
 962               t_molinfo    **molinfo,
 963               t_inputrec   *ir,
 964               int          *nmolblock,
 965               gmx_molblock_t **molblock,
 966               gmx_bool          bGenborn,
 967               warninp_t     wi)
 968 {
 969     /* Tmpfile might contain a long path */
 970     const char *tmpfile;
 971     char **title;
 972
 973     if (topppfile)
 974     {
 975         tmpfile = topppfile;
 976     }
 977     else
 978     {
 979         tmpfile = NULL;
 980     }
 981
 982     if (bVerbose)
 983     {
 984         printf("processing topology...\n");
 985     }
 986     title = read_topol(topfile,tmpfile,opts->define,opts->include,
 987                        symtab,atype,nrmols,molinfo,
 988                        plist,combination_rule,repulsion_power,
 989                        opts,fudgeQQ,nmolblock,molblock,
 990                        ir->efep!=efepNO,bGenborn,bZero,bVerbose,
 991                        wi);
 992     if ((*combination_rule != eCOMB_GEOMETRIC) &&
 993         (ir->vdwtype == evdwUSER))
 994     {
 995         warning(wi,"Using sigma/epsilon based combination rules with"
 996                 " user supplied potential function may produce unwanted"
 997                 " results");
 998     }
 999
1000     return title;
1001 }
1002
1003
1004 static void generate_qmexcl_moltype(gmx_moltype_t *molt,unsigned char *grpnr,
1005                                     t_inputrec *ir)
1006 {
1007   /* This routine expects molt->ilist to be of size F_NRE and ordered. */
1008
1009   /* generates the exclusions between the individual QM atoms, as
1010    * these interactions should be handled by the QM subroutines and
1011    * not by the gromacs routines
1012    */
1013   int
1014     i,j,l,k=0,jmax,qm_max=0,qm_nr=0,nratoms=0,link_nr=0,link_max=0;
1015   atom_id
1016     *qm_arr=NULL,*link_arr=NULL,a1,a2,a3,a4,ftype=0;
1017   t_blocka
1018     qmexcl;
1019   t_block2
1020     qmexcl2;
1021   gmx_bool
1022     *bQMMM,*blink,bexcl;
1023
1024   /* First we search and select the QM atoms in an qm_arr array that
1025    * we use to create the exclusions.
1026    *
1027    * we take the possibility into account that a user has defined more
1028    * than one QM group:
1029    *
1030    * for that we also need to do this an ugly work-about just in case
1031    * the QM group contains the entire system...
1032    */
1033   fprintf(stderr,"excluding classical QM-QM interactions...\n");
1034
1035   jmax = ir->opts.ngQM;
1036
1037   /* we first search for all the QM atoms and put them in an array
1038    */
1039   for(j=0;j<jmax;j++){
1040     for(i=0;i<molt->atoms.nr;i++){
1041       if(qm_nr>=qm_max){
1042         qm_max += 100;
1043         srenew(qm_arr,qm_max);
1044       }
1045       if ((grpnr ? grpnr[i] : 0) == j){
1046         qm_arr[qm_nr++] = i;
1047       }
1048     }
1049   }
1050   /* bQMMM[..] is an array containin TRUE/FALSE for atoms that are
1051    * QM/not QM. We first set all elements to false. Afterwards we use
1052    * the qm_arr to change the elements corresponding to the QM atoms
1053    * to TRUE.
1054    */
1055   snew(bQMMM,molt->atoms.nr);
1056   for (i=0;i<molt->atoms.nr;i++)
1057     bQMMM[i]=FALSE;
1058   for (i=0;i<qm_nr;i++)
1059     bQMMM[qm_arr[i]]=TRUE;
1060
1061   /* We remove all bonded interactions (i.e. bonds,
1062    * angles, dihedrals, 1-4's), involving the QM atoms. The way they
1063    * are removed is as follows: if the interaction invloves 2 atoms,
1064    * it is removed if both atoms are QMatoms. If it involves 3 atoms,
1065    * it is removed if at least two of the atoms are QM atoms, if the
1066    * interaction involves 4 atoms, it is removed if there are at least
1067    * 2 QM atoms.  Since this routine is called once before any forces
1068    * are computed, the top->idef.il[N].iatom[] array (see idef.h) can
1069    * be rewritten at this poitn without any problem. 25-9-2002 */
1070
1071   /* first check weter we already have CONNBONDS: */
1072   if (molt->ilist[F_CONNBONDS].nr != 0){
1073     fprintf(stderr,"nr. of CONNBONDS present already: %d\n",
1074             molt->ilist[F_CONNBONDS].nr/3);
1075     ftype = molt->ilist[F_CONNBONDS].iatoms[0];
1076     k = molt->ilist[F_CONNBONDS].nr;
1077   }
1078   /* now we delete all bonded interactions, except the ones describing
1079    * a chemical bond. These are converted to CONNBONDS
1080    */
1081   for (i=0;i<F_LJ;i++){
1082     if(i==F_CONNBONDS)
1083       continue;
1084     nratoms = interaction_function[i].nratoms;
1085     j = 0;
1086     while (j<molt->ilist[i].nr){
1087       bexcl = FALSE;
1088       switch (nratoms){
1089       case 2:
1090         a1 = molt->ilist[i].iatoms[j+1];
1091         a2 = molt->ilist[i].iatoms[j+2];
1092         bexcl = (bQMMM[a1] && bQMMM[a2]);
1093         /* a bonded beteen two QM atoms will be copied to the
1094          * CONNBONDS list, for reasons mentioned above
1095          */
1096         if(bexcl && i<F_ANGLES){
1097           srenew(molt->ilist[F_CONNBONDS].iatoms,k+3);
1098           molt->ilist[F_CONNBONDS].nr         += 3;
1099           molt->ilist[F_CONNBONDS].iatoms[k++] = ftype;
1100           molt->ilist[F_CONNBONDS].iatoms[k++] = a1;
1101           molt->ilist[F_CONNBONDS].iatoms[k++] = a2;
1102         }
1103         break;
1104       case 3:
1105         a1 = molt->ilist[i].iatoms[j+1];
1106         a2 = molt->ilist[i].iatoms[j+2];
1107         a3 = molt->ilist[i].iatoms[j+3];
1108         bexcl = ((bQMMM[a1] && bQMMM[a2]) ||
1109                  (bQMMM[a1] && bQMMM[a3]) ||
1110                  (bQMMM[a2] && bQMMM[a3]));
1111         break;
1112       case 4:
1113         a1 = molt->ilist[i].iatoms[j+1];
1114         a2 = molt->ilist[i].iatoms[j+2];
1115         a3 = molt->ilist[i].iatoms[j+3];
1116         a4 = molt->ilist[i].iatoms[j+4];
1117         bexcl = ((bQMMM[a1] && bQMMM[a2] && bQMMM[a3]) ||
1118                  (bQMMM[a1] && bQMMM[a2] && bQMMM[a4]) ||
1119                  (bQMMM[a1] && bQMMM[a3] && bQMMM[a4]) ||
1120                  (bQMMM[a2] && bQMMM[a3] && bQMMM[a4]));
1121         break;
1122       default:
1123         gmx_fatal(FARGS,"no such bonded interactions with %d atoms\n",nratoms);
1124       }
1125       if (bexcl){
1126         /* since the interaction involves QM atoms, these should be
1127          * removed from the MM ilist
1128          */
1129         molt->ilist[i].nr -= (nratoms+1);
1130         for (l=j;l<molt->ilist[i].nr;l++)
1131           molt->ilist[i].iatoms[l] = molt->ilist[i].iatoms[l+(nratoms+1)];
1132       } else {
1133         j += nratoms+1; /* the +1 is for the functype */
1134       }
1135     }
1136   }
1137   /* Now, we search for atoms bonded to a QM atom because we also want
1138    * to exclude their nonbonded interactions with the QM atoms. The
1139    * reason for this is that this interaction is accounted for in the
1140    * linkatoms interaction with the QMatoms and would be counted
1141    * twice.  */
1142
1143   for(i=0;i<F_NRE;i++){
1144     if(IS_CHEMBOND(i)){
1145       j=0;
1146       while(j<molt->ilist[i].nr){
1147         a1 = molt->ilist[i].iatoms[j+1];
1148         a2 = molt->ilist[i].iatoms[j+2];
1149         if((bQMMM[a1] && !bQMMM[a2])||(!bQMMM[a1] && bQMMM[a2])){
1150           if(link_nr>=link_max){
1151             link_max += 10;
1152             srenew(link_arr,link_max);
1153           }
1154           if(bQMMM[a1]){
1155             link_arr[link_nr++] = a2;
1156           } else {
1157             link_arr[link_nr++] = a1;
1158           }
1159         }
1160         j+=3;
1161       }
1162     }
1163   }
1164   snew(blink,molt->atoms.nr);
1165   for (i=0;i<molt->atoms.nr;i++)
1166     blink[i]=FALSE;
1167   for (i=0;i<link_nr;i++)
1168     blink[link_arr[i]]=TRUE;
1169   /* creating the exclusion block for the QM atoms. Each QM atom has
1170    * as excluded elements all the other QMatoms (and itself).
1171    */
1172   qmexcl.nr = molt->atoms.nr;
1173   qmexcl.nra = qm_nr*(qm_nr+link_nr)+link_nr*qm_nr;
1174   snew(qmexcl.index,qmexcl.nr+1);
1175   snew(qmexcl.a,qmexcl.nra);
1176   j=0;
1177   for(i=0;i<qmexcl.nr;i++){
1178     qmexcl.index[i]=j;
1179     if(bQMMM[i]){
1180       for(k=0;k<qm_nr;k++){
1181         qmexcl.a[k+j] = qm_arr[k];
1182       }
1183       for(k=0;k<link_nr;k++){
1184         qmexcl.a[qm_nr+k+j] = link_arr[k];
1185       }
1186       j+=(qm_nr+link_nr);
1187     }
1188     if(blink[i]){
1189       for(k=0;k<qm_nr;k++){
1190         qmexcl.a[k+j] = qm_arr[k];
1191       }
1192       j+=qm_nr;
1193     }
1194   }
1195   qmexcl.index[qmexcl.nr]=j;
1196
1197   /* and merging with the exclusions already present in sys.
1198    */
1199
1200   init_block2(&qmexcl2,molt->atoms.nr);
1201   b_to_b2(&qmexcl, &qmexcl2);
1202   merge_excl(&(molt->excls),&qmexcl2);
1203   done_block2(&qmexcl2);
1204
1205   /* Finally, we also need to get rid of the pair interactions of the
1206    * classical atom bonded to the boundary QM atoms with the QMatoms,
1207    * as this interaction is already accounted for by the QM, so also
1208    * here we run the risk of double counting! We proceed in a similar
1209    * way as we did above for the other bonded interactions: */
1210   for (i=F_LJ14;i<F_COUL14;i++){
1211     nratoms = interaction_function[i].nratoms;
1212     j = 0;
1213     while (j<molt->ilist[i].nr){
1214       bexcl = FALSE;
1215       a1 = molt->ilist[i].iatoms[j+1];
1216       a2 = molt->ilist[i].iatoms[j+2];
1217       bexcl = ((bQMMM[a1] && bQMMM[a2])||
1218                (blink[a1] && bQMMM[a2])||
1219                (bQMMM[a1] && blink[a2]));
1220       if (bexcl){
1221         /* since the interaction involves QM atoms, these should be
1222          * removed from the MM ilist
1223          */
1224         molt->ilist[i].nr -= (nratoms+1);
1225         for (k=j;k<molt->ilist[i].nr;k++)
1226           molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
1227       } else {
1228         j += nratoms+1; /* the +1 is for the functype */
1229       }
1230     }
1231   }
1232
1233   free(qm_arr);
1234   free(bQMMM);
1235   free(link_arr);
1236   free(blink);
1237 } /* generate_qmexcl */
1238
1239 void generate_qmexcl(gmx_mtop_t *sys,t_inputrec *ir)
1240 {
1241   /* This routine expects molt->molt[m].ilist to be of size F_NRE and ordered.
1242    */
1243
1244   unsigned char *grpnr;
1245   int mb,mol,nat_mol,i;
1246   gmx_molblock_t *molb;
1247   gmx_bool bQMMM;
1248
1249   grpnr = sys->groups.grpnr[egcQMMM];
1250
1251   for(mb=0; mb<sys->nmolblock; mb++) {
1252     molb = &sys->molblock[mb];
1253     nat_mol = sys->moltype[molb->type].atoms.nr;
1254     for(mol=0; mol<molb->nmol; mol++) {
1255       bQMMM = FALSE;
1256       for(i=0; i<nat_mol; i++) {
1257         if ((grpnr ? grpnr[i] : 0) < ir->opts.ngQM) {
1258           bQMMM = TRUE;
1259         }
1260       }
1261       if (bQMMM) {
1262         if (molb->nmol > 1) {
1263           /* We need to split this molblock */
1264           if (mol > 0) {
1265             /* Split the molblock at this molecule */
1266             sys->nmolblock++;
1267             srenew(sys->molblock,sys->nmolblock);
1268             for(i=mb; i<sys->nmolblock-1; i++) {
1269               sys->molblock[i+1] = sys->molblock[i];
1270             }
1271             sys->molblock[mb  ].nmol  = mol;
1272             sys->molblock[mb+1].nmol -= mol;
1273             mb++;
1274             molb = &sys->molblock[mb];
1275           }
1276           if (molb->nmol > 1) {
1277             /* Split the molblock after this molecule */
1278             sys->nmolblock++;
1279             srenew(sys->molblock,sys->nmolblock);
1280             for(i=mb; i<sys->nmolblock-1; i++) {
1281               sys->molblock[i+1] = sys->molblock[i];
1282             }
1283             sys->molblock[mb  ].nmol  = 1;
1284             sys->molblock[mb+1].nmol -= 1;
1285           }
1286
1287           /* Add a moltype for the QMMM molecule */
1288           sys->nmoltype++;
1289           srenew(sys->moltype,sys->nmoltype);
1290           /* Copy the moltype struct */
1291           sys->moltype[sys->nmoltype-1] = sys->moltype[molb->type];
1292           /* Copy the exclusions to a new array, since this is the only
1293            * thing that needs to be modified for QMMM.
1294            */
1295           copy_blocka(&sys->moltype[molb->type     ].excls,
1296                       &sys->moltype[sys->nmoltype-1].excls);
1297           /* Set the molecule type for the QMMM molblock */
1298           molb->type = sys->nmoltype - 1;
1299         }
1300
1301         generate_qmexcl_moltype(&sys->moltype[molb->type],grpnr,ir);
1302       }
1303       if (grpnr) {
1304         grpnr += nat_mol;
1305       }
1306     }
1307   }
1308 }