src/gromacs/gmxpreprocess/topio.c

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