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