src/kernel/topio.c

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