src/kernel/membed.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-2012, The GROMACS development team,
   6  * check out http://www.gromacs.org for more information.
   7  * Copyright (c) 2012,2013, 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 <signal.h>
  43 #include <stdlib.h>
  44 #include "typedefs.h"
  45 #include "smalloc.h"
  46 #include "sysstuff.h"
  47 #include "vec.h"
  48 #include "statutil.h"
  49 #include "macros.h"
  50 #include "copyrite.h"
  51 #include "main.h"
  52 #include "futil.h"
  53 #include "edsam.h"
  54 #include "index.h"
  55 #include "physics.h"
  56 #include "names.h"
  57 #include "mtop_util.h"
  58 #include "tpxio.h"
  59 #include "string2.h"
  60 #include "membed.h"
  61 #include "pbc.h"
  62 #include "readinp.h"
  63 #include "readir.h"
  64
  65 /* information about scaling center */
  66 typedef struct {
  67     rvec    xmin;         /* smallest coordinates of all embedded molecules */
  68     rvec    xmax;         /* largest coordinates of all embedded molecules */
  69     rvec    *geom_cent;   /* scaling center of each independent molecule to embed */
  70     int     pieces;       /* number of molecules to embed independently */
  71     int     *nidx;        /* n atoms for every independent embedded molecule (index in subindex) */
  72     atom_id **subindex;   /* atomids for independent molecule *
  73                            * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
  74 } pos_ins_t;
  75
  76 /* variables needed in do_md */
  77 struct membed {
  78     int   it_xy;          /* number of iterations (steps) used to grow something in the xy-plane */
  79     int   it_z;           /* same, but for z */
  80     real  xy_step;        /* stepsize used during resize in xy-plane */
  81     real  z_step;         /* same, but in z */
  82     rvec  fac;            /* initial scaling of the molecule to grow into the membrane */
  83     rvec  *r_ins;         /* final coordinates of the molecule to grow  */
  84     pos_ins_t *pos_ins;   /* scaling center for each piece to embed */
  85 };
  86
  87 /* membrane related variables */
  88 typedef struct {
  89     char      *name;     /* name of index group to embed molecule into (usually membrane) */
  90     t_block   mem_at;    /* list all atoms in membrane */
  91     int       nmol;      /* number of membrane molecules overlapping with the molecule to embed */
  92     int       *mol_id;   /* list of molecules in membrane that overlap with the molecule to embed */
  93     real      lip_area;  /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
  94     real      zmin;      /* minimum z coordinate of membrane */
  95     real      zmax;      /* maximum z coordinate of membrane */
  96     real      zmed;      /* median z coordinate of membrane */
  97 } mem_t;
  98
  99 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
 100  * These will then be removed from the system */
 101 typedef struct {
 102     int  nr;      /* number of molecules to remove */
 103     int  *mol;    /* list of molecule ids to remove */
 104     int  *block;  /* id of the molblock that the molecule to remove is part of */
 105 } rm_t;
 106
 107 /* Get the global molecule id, and the corresponding molecule type and id of the *
 108  * molblock from the global atom nr. */
 109 static int get_mol_id(int at, gmx_mtop_t  *mtop, int *type, int *block)
 110 {
 111     int mol_id=0;
 112     int i;
 113     int atnr_mol;
 114     gmx_mtop_atomlookup_t alook;
 115
 116     alook = gmx_mtop_atomlookup_settle_init(mtop);
 117     gmx_mtop_atomnr_to_molblock_ind(alook,at,block,&mol_id,&atnr_mol);
 118     for(i=0;i<*block;i++)
 119     {
 120         mol_id += mtop->molblock[i].nmol;
 121     }
 122     *type = mtop->molblock[*block].type;
 123
 124     gmx_mtop_atomlookup_destroy(alook);
 125
 126     return mol_id;
 127 }
 128
 129 /* Get the id of the molblock from a global molecule id */
 130 static int get_molblock(int mol_id,int nmblock,gmx_molblock_t *mblock)
 131 {
 132     int i;
 133     int nmol=0;
 134
 135     for(i=0;i<nmblock;i++)
 136     {
 137         nmol+=mblock[i].nmol;
 138         if(mol_id<nmol)
 139         {
 140             return i;
 141         }
 142     }
 143
 144     gmx_fatal(FARGS,"mol_id %d larger than total number of molecules %d.\n",mol_id,nmol);
 145
 146     return -1;
 147 }
 148
 149 static int get_tpr_version(const char *infile)
 150 {
 151     t_tpxheader  header;
 152     int          version,generation;
 153
 154     read_tpxheader(infile,&header,TRUE,&version,&generation);
 155
 156     return version;
 157 }
 158
 159 /* Get a list of all the molecule types that are present in a group of atoms. *
 160  * Because all interaction within the group to embed are removed on the topology *
 161  * level, if the same molecule type is found in another part of the system, these *
 162  * would also be affected. Therefore we have to check if the embedded and rest group *
 163  * share common molecule types. If so, membed will stop with an error. */
 164 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
 165 {
 166     int i,j,nr,mol_id;
 167     int type=0,block=0;
 168     gmx_bool bNEW;
 169
 170     nr=0;
 171     snew(tlist->index,at->nr);
 172     for (i=0;i<at->nr;i++)
 173     {
 174         bNEW=TRUE;
 175         mol_id = get_mol_id(at->index[i],mtop,&type,&block);
 176         for(j=0;j<nr;j++)
 177         {
 178             if(tlist->index[j]==type)
 179             {
 180                 bNEW=FALSE;
 181             }
 182         }
 183
 184         if(bNEW==TRUE)
 185         {
 186             tlist->index[nr]=type;
 187             nr++;
 188         }
 189     }
 190     srenew(tlist->index,nr);
 191     return nr;
 192 }
 193
 194 /* Do the actual check of the molecule types between embedded and rest group */
 195 static void check_types(t_block *ins_at,t_block *rest_at,gmx_mtop_t *mtop)
 196 {
 197     t_block        *ins_mtype,*rest_mtype;
 198     int            i,j;
 199
 200     snew(ins_mtype,1);
 201     snew(rest_mtype,1);
 202     ins_mtype->nr  = get_mtype_list(ins_at , mtop, ins_mtype );
 203     rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
 204
 205     for(i=0;i<ins_mtype->nr;i++)
 206     {
 207         for(j=0;j<rest_mtype->nr;j++)
 208         {
 209             if(ins_mtype->index[i]==rest_mtype->index[j])
 210             {
 211                 gmx_fatal(FARGS,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
 212                           "1. Your *.ndx and *.top do not match\n"
 213                           "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
 214                           "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
 215                           "we need to exclude all interactions between the atoms in the group to\n"
 216                           "insert, the same moleculetype can not be used in both groups. Change the\n"
 217                           "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
 218                           "an appropriate *.itp file",*(mtop->moltype[rest_mtype->index[j]].name),
 219                           *(mtop->moltype[rest_mtype->index[j]].name),*(mtop->moltype[rest_mtype->index[j]].name));
 220             }
 221         }
 222     }
 223
 224     sfree(ins_mtype->index);
 225     sfree(rest_mtype->index);
 226     sfree(ins_mtype);
 227     sfree(rest_mtype);
 228 }
 229
 230 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
 231                       int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
 232                       int *pieces, gmx_bool *bALLOW_ASYMMETRY)
 233 {
 234     warninp_t wi;
 235     t_inpfile *inp;
 236     int       ninp;
 237
 238     wi = init_warning(TRUE,0);
 239
 240     inp = read_inpfile(membed_input, &ninp, NULL, wi);
 241     ITYPE ("nxy", *it_xy, 1000);
 242     ITYPE ("nz", *it_z, 0);
 243     RTYPE ("xyinit", *xy_fac, 0.5);
 244     RTYPE ("xyend", *xy_max, 1.0);
 245     RTYPE ("zinit", *z_fac, 1.0);
 246     RTYPE ("zend", *z_max, 1.0);
 247     RTYPE ("rad", *probe_rad, 0.22);
 248     ITYPE ("ndiff", *low_up_rm, 0);
 249     ITYPE ("maxwarn", *maxwarn, 0);
 250     ITYPE ("pieces", *pieces, 1);
 251     EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
 252
 253     write_inpfile(membed_input,ninp,inp,FALSE,wi);
 254 }
 255
 256 /* Obtain the maximum and minimum coordinates of the group to be embedded */
 257 static int init_ins_at(t_block *ins_at,t_block *rest_at,t_state *state, pos_ins_t *pos_ins,
 258                        gmx_groups_t *groups,int ins_grp_id, real xy_max)
 259 {
 260     int i,gid,c=0;
 261     real x,xmin,xmax,y,ymin,ymax,z,zmin,zmax;
 262     const real min_memthick=6.0;      /* minimum thickness of the bilayer that will be used to *
 263                                        * determine the overlap between molecule to embed and membrane */
 264     const real fac_inp_size=1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
 265     snew(rest_at->index,state->natoms);
 266
 267     xmin=xmax=state->x[ins_at->index[0]][XX];
 268     ymin=ymax=state->x[ins_at->index[0]][YY];
 269     zmin=zmax=state->x[ins_at->index[0]][ZZ];
 270
 271     for(i=0;i<state->natoms;i++)
 272     {
 273         gid = groups->grpnr[egcFREEZE][i];
 274         if(groups->grps[egcFREEZE].nm_ind[gid]==ins_grp_id)
 275         {
 276             x=state->x[i][XX];
 277             if (x<xmin)                xmin=x;
 278             if (x>xmax)                xmax=x;
 279             y=state->x[i][YY];
 280             if (y<ymin)                ymin=y;
 281             if (y>ymax)                ymax=y;
 282             z=state->x[i][ZZ];
 283             if (z<zmin)                zmin=z;
 284             if (z>zmax)                zmax=z;
 285         }
 286         else
 287         {
 288             rest_at->index[c]=i;
 289             c++;
 290         }
 291     }
 292
 293     rest_at->nr=c;
 294     srenew(rest_at->index,c);
 295
 296     if(xy_max>fac_inp_size)
 297     {
 298         pos_ins->xmin[XX]=xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
 299         pos_ins->xmin[YY]=ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
 300
 301         pos_ins->xmax[XX]=xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
 302         pos_ins->xmax[YY]=ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
 303     }
 304     else
 305     {
 306         pos_ins->xmin[XX]=xmin;
 307         pos_ins->xmin[YY]=ymin;
 308
 309         pos_ins->xmax[XX]=xmax;
 310         pos_ins->xmax[YY]=ymax;
 311     }
 312
 313     if( (zmax-zmin) < min_memthick )
 314     {
 315         pos_ins->xmin[ZZ]=zmin+(zmax-zmin)/2.0-0.5*min_memthick;
 316         pos_ins->xmax[ZZ]=zmin+(zmax-zmin)/2.0+0.5*min_memthick;
 317     }
 318     else
 319     {
 320         pos_ins->xmin[ZZ]=zmin;
 321         pos_ins->xmax[ZZ]=zmax;
 322     }
 323
 324     return c;
 325 }
 326
 327 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
 328  * xy-plane and counting the number of occupied grid points */
 329 static real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
 330 {
 331     real x,y,dx=0.15,dy=0.15,area=0.0;
 332     real add,memmin,memmax;
 333     int c,at;
 334
 335     /* min and max membrane coordinate are altered to reduce the influence of the *
 336      * boundary region */
 337     memmin=mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
 338     memmax=mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
 339
 340     for(x=pos_ins->xmin[XX];x<pos_ins->xmax[XX];x+=dx)
 341     {
 342         for(y=pos_ins->xmin[YY];y<pos_ins->xmax[YY];y+=dy)
 343         {
 344             c=0;
 345             add=0.0;
 346             do
 347             {
 348                 at=ins_at->index[c];
 349                 if ( (r[at][XX]>=x) && (r[at][XX]<x+dx) &&
 350                      (r[at][YY]>=y) && (r[at][YY]<y+dy) &&
 351                      (r[at][ZZ]>memmin) && (r[at][ZZ]<memmax) )
 352                 {
 353                     add=1.0;
 354                 }
 355                 c++;
 356             } while ( (c<ins_at->nr) && (add<0.5) );
 357             area+=add;
 358         }
 359     }
 360     area=area*dx*dy;
 361
 362     return area;
 363 }
 364
 365 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
 366 {
 367     int i,j,at,mol,nmol,nmolbox,count;
 368     t_block *mem_a;
 369     real z,zmin,zmax,mem_area;
 370     gmx_bool bNew;
 371     atom_id *mol_id;
 372     int type=0,block=0;
 373
 374     nmol=count=0;
 375     mem_a=&(mem_p->mem_at);
 376     snew(mol_id,mem_a->nr);
 377     zmin=pos_ins->xmax[ZZ];
 378     zmax=pos_ins->xmin[ZZ];
 379     for(i=0;i<mem_a->nr;i++)
 380     {
 381         at=mem_a->index[i];
 382         if( (r[at][XX]>pos_ins->xmin[XX]) && (r[at][XX]<pos_ins->xmax[XX]) &&
 383             (r[at][YY]>pos_ins->xmin[YY]) && (r[at][YY]<pos_ins->xmax[YY]) &&
 384             (r[at][ZZ]>pos_ins->xmin[ZZ]) && (r[at][ZZ]<pos_ins->xmax[ZZ]) )
 385         {
 386             mol = get_mol_id(at,mtop,&type,&block);
 387             bNew=TRUE;
 388             for(j=0;j<nmol;j++)
 389             {
 390                 if(mol == mol_id[j])
 391                 {
 392                     bNew=FALSE;
 393                 }
 394             }
 395
 396             if(bNew)
 397             {
 398                 mol_id[nmol]=mol;
 399                 nmol++;
 400             }
 401
 402             z=r[at][ZZ];
 403             if(z<zmin)
 404             {
 405                 zmin=z;
 406             }
 407
 408             if(z>zmax)
 409             {
 410                 zmax=z;
 411             }
 412
 413             count++;
 414         }
 415     }
 416
 417     mem_p->nmol=nmol;
 418     srenew(mol_id,nmol);
 419     mem_p->mol_id=mol_id;
 420
 421     if((zmax-zmin)>(box[ZZ][ZZ]-0.5))
 422     {
 423         gmx_fatal(FARGS,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
 424                   "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
 425                   "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
 426                   "your water layer is not thick enough.\n",zmax,zmin);
 427     }
 428     mem_p->zmin=zmin;
 429     mem_p->zmax=zmax;
 430     mem_p->zmed=(zmax-zmin)/2+zmin;
 431
 432     /*number of membrane molecules in protein box*/
 433     nmolbox = count/mtop->molblock[block].natoms_mol;
 434     /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
 435     mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
 436     /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
 437     mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
 438
 439     return mem_p->mem_at.nr;
 440 }
 441
 442 static void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r,
 443                         gmx_bool bALLOW_ASYMMETRY)
 444 {
 445     int i,j,at,c,outsidesum,gctr=0;
 446     int idxsum=0;
 447
 448     /*sanity check*/
 449     for (i=0;i<pos_ins->pieces;i++)
 450     {
 451         idxsum+=pos_ins->nidx[i];
 452     }
 453
 454     if (idxsum!=ins_at->nr)
 455     {
 456         gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
 457     }
 458
 459     snew(pos_ins->geom_cent,pos_ins->pieces);
 460     for (i=0;i<pos_ins->pieces;i++)
 461     {
 462         c=0;
 463         outsidesum=0;
 464         for(j=0;j<DIM;j++)
 465         {
 466             pos_ins->geom_cent[i][j]=0;
 467         }
 468
 469         for (j=0;j<pos_ins->nidx[i];j++)
 470         {
 471             at=pos_ins->subindex[i][j];
 472             copy_rvec(r[at],r_ins[gctr]);
 473             if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )
 474             {
 475                 rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);
 476                 c++;
 477             }
 478             else
 479             {
 480                 outsidesum++;
 481             }
 482             gctr++;
 483         }
 484
 485         if (c>0)
 486         {
 487             svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);
 488         }
 489
 490         if (!bALLOW_ASYMMETRY)
 491         {
 492             pos_ins->geom_cent[i][ZZ]=mem_p->zmed;
 493         }
 494
 495         fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f\n",
 496                 i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);
 497     }
 498     fprintf(stderr,"\n");
 499 }
 500
 501 /* resize performed in the md loop */
 502 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins,rvec fac)
 503 {
 504     int i,j,k,at,c=0;
 505     for (k=0;k<pos_ins->pieces;k++)
 506     {
 507         for(i=0;i<pos_ins->nidx[k];i++)
 508         {
 509             at=pos_ins->subindex[k][i];
 510             for(j=0;j<DIM;j++)
 511             {
 512                 r[at][j]=pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
 513             }
 514             c++;
 515         }
 516     }
 517 }
 518
 519 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
 520  * The molecule to be embedded is already reduced in size. */
 521 static int gen_rm_list(rm_t *rm_p,t_block *ins_at,t_block *rest_at,t_pbc *pbc, gmx_mtop_t *mtop,
 522                        rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
 523                        int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
 524 {
 525     int i,j,k,l,at,at2,mol_id;
 526     int type=0,block=0;
 527     int nrm,nupper,nlower;
 528     real r_min_rad,z_lip,min_norm;
 529     gmx_bool bRM;
 530     rvec dr,dr_tmp;
 531     real *dist;
 532     int *order;
 533
 534     r_min_rad=probe_rad*probe_rad;
 535     snew(rm_p->mol,mtop->mols.nr);
 536     snew(rm_p->block,mtop->mols.nr);
 537     nrm=nupper=0;
 538     nlower=low_up_rm;
 539     for(i=0;i<ins_at->nr;i++)
 540     {
 541         at=ins_at->index[i];
 542         for(j=0;j<rest_at->nr;j++)
 543         {
 544             at2=rest_at->index[j];
 545             pbc_dx(pbc,r[at],r[at2],dr);
 546
 547             if(norm2(dr)<r_min_rad)
 548             {
 549                 mol_id = get_mol_id(at2,mtop,&type,&block);
 550                 bRM=TRUE;
 551                 for(l=0;l<nrm;l++)
 552                 {
 553                     if(rm_p->mol[l]==mol_id)
 554                     {
 555                         bRM=FALSE;
 556                     }
 557                 }
 558
 559                 if(bRM)
 560                 {
 561                     rm_p->mol[nrm]=mol_id;
 562                     rm_p->block[nrm]=block;
 563                     nrm++;
 564                     z_lip=0.0;
 565                     for(l=0;l<mem_p->nmol;l++)
 566                     {
 567                         if(mol_id==mem_p->mol_id[l])
 568                         {
 569                             for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
 570                             {
 571                                 z_lip+=r[k][ZZ];
 572                             }
 573                             z_lip/=mtop->molblock[block].natoms_mol;
 574                             if(z_lip<mem_p->zmed)
 575                             {
 576                                 nlower++;
 577                             }
 578                             else
 579                             {
 580                                 nupper++;
 581                             }
 582                         }
 583                     }
 584                 }
 585             }
 586         }
 587     }
 588
 589     /*make sure equal number of lipids from upper and lower layer are removed */
 590     if( (nupper!=nlower) && (!bALLOW_ASYMMETRY) )
 591     {
 592         snew(dist,mem_p->nmol);
 593         snew(order,mem_p->nmol);
 594         for(i=0;i<mem_p->nmol;i++)
 595         {
 596             at = mtop->mols.index[mem_p->mol_id[i]];
 597             pbc_dx(pbc,r[at],pos_ins->geom_cent[0],dr);
 598             if (pos_ins->pieces>1)
 599             {
 600                 /*minimum dr value*/
 601                 min_norm=norm2(dr);
 602                 for (k=1;k<pos_ins->pieces;k++)
 603                 {
 604                     pbc_dx(pbc,r[at],pos_ins->geom_cent[k],dr_tmp);
 605                     if (norm2(dr_tmp) < min_norm)
 606                     {
 607                         min_norm=norm2(dr_tmp);
 608                         copy_rvec(dr_tmp,dr);
 609                     }
 610                 }
 611             }
 612             dist[i]=dr[XX]*dr[XX]+dr[YY]*dr[YY];
 613             j=i-1;
 614             while (j>=0 && dist[i]<dist[order[j]])
 615             {
 616                 order[j+1]=order[j];
 617                 j--;
 618             }
 619             order[j+1]=i;
 620         }
 621
 622         i=0;
 623         while(nupper!=nlower)
 624         {
 625             mol_id=mem_p->mol_id[order[i]];
 626             block=get_molblock(mol_id,mtop->nmolblock,mtop->molblock);
 627             bRM=TRUE;
 628             for(l=0;l<nrm;l++)
 629             {
 630                 if(rm_p->mol[l]==mol_id)
 631                 {
 632                     bRM=FALSE;
 633                 }
 634             }
 635
 636             if(bRM)
 637             {
 638                 z_lip=0;
 639                 for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
 640                 {
 641                     z_lip+=r[k][ZZ];
 642                 }
 643                 z_lip/=mtop->molblock[block].natoms_mol;
 644                 if(nupper>nlower && z_lip<mem_p->zmed)
 645                 {
 646                     rm_p->mol[nrm]=mol_id;
 647                     rm_p->block[nrm]=block;
 648                     nrm++;
 649                     nlower++;
 650                 }
 651                 else if (nupper<nlower && z_lip>mem_p->zmed)
 652                 {
 653                     rm_p->mol[nrm]=mol_id;
 654                     rm_p->block[nrm]=block;
 655                     nrm++;
 656                     nupper++;
 657                 }
 658             }
 659             i++;
 660
 661             if(i>mem_p->nmol)
 662             {
 663                 gmx_fatal(FARGS,"Trying to remove more lipid molecules than there are in the membrane");
 664             }
 665         }
 666         sfree(dist);
 667         sfree(order);
 668     }
 669
 670     rm_p->nr=nrm;
 671     srenew(rm_p->mol,nrm);
 672     srenew(rm_p->block,nrm);
 673
 674     return nupper+nlower;
 675 }
 676
 677 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
 678 static void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
 679                      t_block *ins_at, pos_ins_t *pos_ins)
 680 {
 681     int i,j,k,n,rm,mol_id,at,block;
 682     rvec *x_tmp,*v_tmp;
 683     atom_id *list,*new_mols;
 684     unsigned char  *new_egrp[egcNR];
 685     gmx_bool bRM;
 686     int RMmolblock;
 687
 688     snew(list,state->natoms);
 689     n=0;
 690     for(i=0;i<rm_p->nr;i++)
 691     {
 692         mol_id=rm_p->mol[i];
 693         at=mtop->mols.index[mol_id];
 694         block =rm_p->block[i];
 695         mtop->molblock[block].nmol--;
 696         for(j=0;j<mtop->molblock[block].natoms_mol;j++)
 697         {
 698             list[n]=at+j;
 699             n++;
 700         }
 701         mtop->mols.index[mol_id]=-1;
 702     }
 703
 704     mtop->mols.nr-=rm_p->nr;
 705     mtop->mols.nalloc_index-=rm_p->nr;
 706     snew(new_mols,mtop->mols.nr);
 707     for(i=0;i<mtop->mols.nr+rm_p->nr;i++)
 708     {
 709         j=0;
 710         if(mtop->mols.index[i]!=-1)
 711         {
 712             new_mols[j]=mtop->mols.index[i];
 713             j++;
 714         }
 715     }
 716     sfree(mtop->mols.index);
 717     mtop->mols.index=new_mols;
 718     mtop->natoms-=n;
 719     state->natoms-=n;
 720     state->nalloc=state->natoms;
 721     snew(x_tmp,state->nalloc);
 722     snew(v_tmp,state->nalloc);
 723
 724     for(i=0;i<egcNR;i++)
 725     {
 726         if(groups->grpnr[i]!=NULL)
 727         {
 728             groups->ngrpnr[i]=state->natoms;
 729             snew(new_egrp[i],state->natoms);
 730         }
 731     }
 732
 733     rm=0;
 734     for (i=0;i<state->natoms+n;i++)
 735     {
 736         bRM=FALSE;
 737         for(j=0;j<n;j++)
 738         {
 739             if(i==list[j])
 740             {
 741                 bRM=TRUE;
 742                 rm++;
 743             }
 744         }
 745
 746         if(!bRM)
 747         {
 748             for(j=0;j<egcNR;j++)
 749             {
 750                 if(groups->grpnr[j]!=NULL)
 751                 {
 752                     new_egrp[j][i-rm]=groups->grpnr[j][i];
 753                 }
 754             }
 755             copy_rvec(state->x[i],x_tmp[i-rm]);
 756             copy_rvec(state->v[i],v_tmp[i-rm]);
 757             for(j=0;j<ins_at->nr;j++)
 758             {
 759                 if (i==ins_at->index[j])
 760                 {
 761                     ins_at->index[j]=i-rm;
 762                 }
 763             }
 764
 765             for(j=0;j<pos_ins->pieces;j++)
 766             {
 767                 for(k=0;k<pos_ins->nidx[j];k++)
 768                 {
 769                     if (i==pos_ins->subindex[j][k])
 770                     {
 771                         pos_ins->subindex[j][k]=i-rm;
 772                     }
 773                 }
 774             }
 775         }
 776     }
 777     sfree(state->x);
 778     state->x=x_tmp;
 779     sfree(state->v);
 780     state->v=v_tmp;
 781
 782     for(i=0;i<egcNR;i++)
 783     {
 784         if(groups->grpnr[i]!=NULL)
 785         {
 786             sfree(groups->grpnr[i]);
 787             groups->grpnr[i]=new_egrp[i];
 788         }
 789     }
 790
 791     /* remove empty molblocks */
 792     RMmolblock=0;
 793     for (i=0;i<mtop->nmolblock;i++)
 794     {
 795         if(mtop->molblock[i].nmol==0)
 796         {
 797             RMmolblock++;
 798         }
 799         else
 800         {
 801             mtop->molblock[i-RMmolblock]=mtop->molblock[i];
 802         }
 803     }
 804     mtop->nmolblock-=RMmolblock;
 805 }
 806
 807 /* remove al bonded interactions from mtop for the molecule to be embedded */
 808 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
 809 {
 810     int i,j,m;
 811     int type,natom,nmol,at,atom1=0,rm_at=0;
 812     gmx_bool *bRM,bINS;
 813     /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
 814     /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
 815      *sure that g_membed exits with a warning when there are molecules of the same type not in the *
 816      *ins_at index group. MGWolf 050710 */
 817
 818
 819     snew(bRM,mtop->nmoltype);
 820     for (i=0;i<mtop->nmoltype;i++)
 821     {
 822         bRM[i]=TRUE;
 823     }
 824
 825     for (i=0;i<mtop->nmolblock;i++)
 826     {
 827         /*loop over molecule blocks*/
 828         type        =mtop->molblock[i].type;
 829         natom        =mtop->molblock[i].natoms_mol;
 830         nmol        =mtop->molblock[i].nmol;
 831
 832         for(j=0;j<natom*nmol && bRM[type]==TRUE;j++)
 833         {
 834             /*loop over atoms in the block*/
 835             at=j+atom1; /*atom index = block index + offset*/
 836             bINS=FALSE;
 837
 838             for (m=0;(m<ins_at->nr) && (bINS==FALSE);m++)
 839             {
 840                 /*loop over atoms in insertion index group to determine if we're inserting one*/
 841                 if(at==ins_at->index[m])
 842                 {
 843                     bINS=TRUE;
 844                 }
 845             }
 846             bRM[type]=bINS;
 847         }
 848         atom1+=natom*nmol; /*update offset*/
 849         if(bRM[type])
 850         {
 851             rm_at+=natom*nmol; /*increment bonded removal counter by # atoms in block*/
 852         }
 853     }
 854
 855     for(i=0;i<mtop->nmoltype;i++)
 856     {
 857         if(bRM[i])
 858         {
 859             for(j=0;j<F_LJ;j++)
 860             {
 861                 mtop->moltype[i].ilist[j].nr=0;
 862             }
 863
 864             for(j=F_POSRES;j<=F_VSITEN;j++)
 865             {
 866                 mtop->moltype[i].ilist[j].nr=0;
 867             }
 868         }
 869     }
 870     sfree(bRM);
 871
 872     return rm_at;
 873 }
 874
 875 /* Write a topology where the number of molecules is correct for the system after embedding */
 876 static void top_update(const char *topfile, char *ins, rm_t *rm_p, gmx_mtop_t *mtop)
 877 {
 878 #define TEMP_FILENM "temp.top"
 879     int    bMolecules=0;
 880     FILE    *fpin,*fpout;
 881     char    buf[STRLEN],buf2[STRLEN],*temp;
 882     int        i,*nmol_rm,nmol,line;
 883
 884     fpin  = ffopen(topfile,"r");
 885     fpout = ffopen(TEMP_FILENM,"w");
 886
 887     snew(nmol_rm,mtop->nmoltype);
 888     for(i=0;i<rm_p->nr;i++)
 889     {
 890         nmol_rm[rm_p->block[i]]++;
 891     }
 892
 893     line=0;
 894     while(fgets(buf,STRLEN,fpin))
 895     {
 896         line++;
 897         if(buf[0]!=';')
 898         {
 899             strcpy(buf2,buf);
 900             if ((temp=strchr(buf2,'\n')) != NULL)
 901             {
 902                 temp[0]='\0';
 903             }
 904             ltrim(buf2);
 905             if (buf2[0]=='[')
 906             {
 907                 buf2[0]=' ';
 908                 if ((temp=strchr(buf2,'\n')) != NULL)
 909                 {
 910                     temp[0]='\0';
 911                 }
 912                 rtrim(buf2);
 913                 if (buf2[strlen(buf2)-1]==']')
 914                 {
 915                     buf2[strlen(buf2)-1]='\0';
 916                     ltrim(buf2);
 917                     rtrim(buf2);
 918                     if (gmx_strcasecmp(buf2,"molecules")==0)
 919                     {
 920                         bMolecules=1;
 921                     }
 922                 }
 923                 fprintf(fpout,"%s",buf);
 924             }
 925             else if (bMolecules==1)
 926             {
 927                 for(i=0;i<mtop->nmolblock;i++)
 928                 {
 929                     nmol=mtop->molblock[i].nmol;
 930                     sprintf(buf,"%-15s %5d\n",*(mtop->moltype[mtop->molblock[i].type].name),nmol);
 931                     fprintf(fpout,"%s",buf);
 932                 }
 933                 bMolecules=2;
 934             }
 935             else if (bMolecules==2)
 936             {
 937                 /* print nothing */
 938             }
 939             else
 940             {
 941                 fprintf(fpout,"%s",buf);
 942             }
 943         }
 944         else
 945         {
 946             fprintf(fpout,"%s",buf);
 947         }
 948     }
 949
 950     ffclose(fpout);
 951     /* use ffopen to generate backup of topinout */
 952     fpout=ffopen(topfile,"w");
 953     ffclose(fpout);
 954     rename(TEMP_FILENM,topfile);
 955 #undef TEMP_FILENM
 956 }
 957
 958 void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x)
 959 {
 960     /* Set new positions for the group to embed */
 961     if(step_rel<=membed->it_xy)
 962     {
 963         membed->fac[0]+=membed->xy_step;
 964         membed->fac[1]+=membed->xy_step;
 965     }
 966     else if (step_rel<=(membed->it_xy+membed->it_z))
 967     {
 968         membed->fac[2]+=membed->z_step;
 969     }
 970     resize(membed->r_ins,x,membed->pos_ins,membed->fac);
 971 }
 972
 973 gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
 974                  t_inputrec *inputrec, t_state *state, t_commrec *cr,real *cpt)
 975 {
 976     char                    *ins,**gnames;
 977     int                     i,rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
 978     int                     ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
 979     real                    prot_area;
 980     rvec                    *r_ins=NULL;
 981     t_block                 *ins_at,*rest_at;
 982     pos_ins_t               *pos_ins;
 983     mem_t                   *mem_p;
 984     rm_t                    *rm_p;
 985     gmx_groups_t            *groups;
 986     gmx_bool                bExcl=FALSE;
 987     t_atoms                 atoms;
 988     t_pbc                   *pbc;
 989     char                    **piecename=NULL;
 990     gmx_membed_t            membed=NULL;
 991
 992     /* input variables */
 993     const char *membed_input;
 994     real xy_fac = 0.5;
 995     real xy_max = 1.0;
 996     real z_fac = 1.0;
 997     real z_max = 1.0;
 998     int it_xy = 1000;
 999     int it_z = 0;
1000     real probe_rad = 0.22;
1001     int low_up_rm = 0;
1002     int maxwarn=0;
1003     int pieces=1;
1004     gmx_bool bALLOW_ASYMMETRY=FALSE;
1005
1006     /* sanity check constants */         /* Issue a warning when: */
1007     const int membed_version=58;         /* tpr version is smaller */
1008     const real min_probe_rad=0.2199999;  /* A probe radius for overlap between embedded molecule *
1009                                           * and rest smaller than this value is probably too small */
1010     const real min_xy_init=0.0999999;    /* the initial shrinking of the molecule to embed is smaller */
1011     const int min_it_xy=1000;            /* the number of steps to embed in xy-plane is smaller */
1012     const int min_it_z=100;              /* the number of steps to embed in z is smaller */
1013     const real prot_vs_box=7.5;          /* molecule to embed is large (more then prot_vs_box) with respect */
1014     const real box_vs_prot=50;           /* to the box size (less than box_vs_prot) */
1015
1016     snew(membed,1);
1017     snew(ins_at,1);
1018     snew(pos_ins,1);
1019
1020     if(MASTER(cr))
1021     {
1022         /* get input data out membed file */
1023         membed_input = opt2fn("-membed",nfile,fnm);
1024         get_input(membed_input,&xy_fac,&xy_max,&z_fac,&z_max,&it_xy,&it_z,&probe_rad,&low_up_rm,
1025                   &maxwarn,&pieces,&bALLOW_ASYMMETRY);
1026
1027         tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
1028         if (tpr_version<membed_version)
1029         {
1030             gmx_fatal(FARGS,"Version of *.tpr file to old (%d). "
1031                             "Rerun grompp with GROMACS version 4.0.3 or newer.\n",tpr_version);
1032         }
1033
1034         if( !EI_DYNAMICS(inputrec->eI) )
1035         {
1036             gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1037         }
1038
1039         if(PAR(cr))
1040         {
1041             gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1042         }
1043
1044         if(*cpt>=0)
1045         {
1046             fprintf(stderr,"\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1047              *cpt=-1;
1048         }
1049         groups=&(mtop->groups);
1050         snew(gnames,groups->ngrpname);
1051         for (i=0; i<groups->ngrpname; i++)
1052         {
1053             gnames[i] = *(groups->grpname[i]);
1054         }
1055
1056         atoms=gmx_mtop_global_atoms(mtop);
1057         snew(mem_p,1);
1058         fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
1059         get_index(&atoms,opt2fn_null("-mn",nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
1060         ins_grp_id = search_string(ins,groups->ngrpname,gnames);
1061         fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
1062         get_index(&atoms,opt2fn_null("-mn",nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
1063
1064         pos_ins->pieces=pieces;
1065         snew(pos_ins->nidx,pieces);
1066         snew(pos_ins->subindex,pieces);
1067         snew(piecename,pieces);
1068         if (pieces>1)
1069         {
1070             fprintf(stderr,"\nSelect pieces to embed:\n");
1071             get_index(&atoms,opt2fn_null("-mn",nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
1072         }
1073         else
1074         {
1075             /*use whole embedded group*/
1076             snew(pos_ins->nidx,1);
1077             snew(pos_ins->subindex,1);
1078             pos_ins->nidx[0]=ins_at->nr;
1079             pos_ins->subindex[0]=ins_at->index;
1080         }
1081
1082         if(probe_rad<min_probe_rad)
1083         {
1084             warn++;
1085             fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1086                            "in overlap between waters and the group to embed, which will result "
1087                            "in Lincs errors etc.\n\n",warn);
1088         }
1089
1090         if(xy_fac<min_xy_init)
1091         {
1092             warn++;
1093             fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n\n",warn,ins);
1094         }
1095
1096         if(it_xy<min_it_xy)
1097         {
1098             warn++;
1099             fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1100                            " is probably too small.\nIncrease -nxy or.\n\n",warn,ins,it_xy);
1101         }
1102
1103         if( (it_z<min_it_z) && ( z_fac<0.99999999 || z_fac>1.0000001) )
1104         {
1105             warn++;
1106             fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1107                            " is probably too small.\nIncrease -nz or maxwarn.\n\n",warn,ins,it_z);
1108         }
1109
1110         if(it_xy+it_z>inputrec->nsteps)
1111         {
1112             warn++;
1113             fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1114                            "number of steps in the tpr.\n\n",warn);
1115         }
1116
1117         fr_id=-1;
1118         if( inputrec->opts.ngfrz==1)
1119         {
1120             gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
1121         }
1122
1123         for(i=0;i<inputrec->opts.ngfrz;i++)
1124         {
1125             tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1126             if(ins_grp_id==tmp_id)
1127             {
1128                 fr_id=tmp_id;
1129                 fr_i=i;
1130             }
1131         }
1132
1133         if (fr_id == -1 )
1134         {
1135             gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
1136         }
1137
1138         for(i=0;i<DIM;i++)
1139         {
1140             if( inputrec->opts.nFreeze[fr_i][i] != 1)
1141             {
1142                 gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
1143             }
1144         }
1145
1146         ng = groups->grps[egcENER].nr;
1147         if (ng == 1)
1148         {
1149             gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1150         }
1151
1152         for(i=0;i<ng;i++)
1153         {
1154             for(j=0;j<ng;j++)
1155             {
1156                 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1157                 {
1158                     bExcl = TRUE;
1159                     if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1160                          (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1161                     {
1162                         gmx_fatal(FARGS,"Energy exclusions \"%s\" and  \"%s\" do not match the group "
1163                                   "to embed \"%s\"",*groups->grpname[groups->grps[egcENER].nm_ind[i]],
1164                                   *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
1165                     }
1166                 }
1167             }
1168         }
1169
1170         if (!bExcl) {
1171             gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1172                       "the freeze group");
1173         }
1174
1175         /* Obtain the maximum and minimum coordinates of the group to be embedded */
1176         snew(rest_at,1);
1177         ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
1178         /* Check that moleculetypes in insertion group are not part of the rest of the system */
1179         check_types(ins_at,rest_at,mtop);
1180
1181         mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
1182
1183         prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
1184         if ( (prot_area>prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<box_vs_prot) )
1185         {
1186             warn++;
1187             fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1188                     "This might cause pressure problems during the growth phase. Just try with\n"
1189                     "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
1190                     "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
1191         }
1192
1193         if(warn>maxwarn)
1194         {
1195             gmx_fatal(FARGS,"Too many warnings.\n");
1196         }
1197
1198         printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
1199         printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1200                "The area per lipid is %.4f nm^2.\n",mem_p->nmol,mem_p->lip_area);
1201
1202         /* Maximum number of lipids to be removed*/
1203         max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
1204         printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
1205
1206         printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1207                "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1208                "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1209                xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
1210
1211         /* resize the protein by xy and by z if necessary*/
1212         snew(r_ins,ins_at->nr);
1213         init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
1214         membed->fac[0]=membed->fac[1]=xy_fac;
1215         membed->fac[2]=z_fac;
1216
1217         membed->xy_step =(xy_max-xy_fac)/(double)(it_xy);
1218         membed->z_step  =(z_max-z_fac)/(double)(it_z-1);
1219
1220         resize(r_ins,state->x,pos_ins,membed->fac);
1221
1222         /* remove overlapping lipids and water from the membrane box*/
1223         /*mark molecules to be removed*/
1224         snew(pbc,1);
1225         set_pbc(pbc,inputrec->ePBC,state->box);
1226
1227         snew(rm_p,1);
1228         lip_rm = gen_rm_list(rm_p,ins_at,rest_at,pbc,mtop,state->x, r_ins, mem_p,pos_ins,
1229                              probe_rad,low_up_rm,bALLOW_ASYMMETRY);
1230         lip_rm -= low_up_rm;
1231
1232         if(fplog)
1233         {
1234             for(i=0;i<rm_p->nr;i++)
1235             {
1236                 fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
1237             }
1238         }
1239
1240         for(i=0;i<mtop->nmolblock;i++)
1241         {
1242             ntype=0;
1243             for(j=0;j<rm_p->nr;j++)
1244             {
1245                 if(rm_p->block[j]==i)
1246                 {
1247                     ntype++;
1248                 }
1249             }
1250             printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
1251         }
1252
1253         if(lip_rm>max_lip_rm)
1254         {
1255             warn++;
1256             fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1257                            "protein area\nTry making the -xyinit resize factor smaller or increase "
1258                            "maxwarn.\n\n",warn);
1259         }
1260
1261         /*remove all lipids and waters overlapping and update all important structures*/
1262         rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
1263
1264         rm_bonded_at = rm_bonded(ins_at,mtop);
1265         if (rm_bonded_at != ins_at->nr)
1266         {
1267             fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
1268                     "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1269                     "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
1270         }
1271
1272         if(warn>maxwarn)
1273         {
1274             gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, "
1275                             "you can increase -maxwarn");
1276         }
1277
1278         if (ftp2bSet(efTOP,nfile,fnm))
1279         {
1280             top_update(opt2fn("-mp",nfile,fnm),ins,rm_p,mtop);
1281         }
1282
1283         sfree(pbc);
1284         sfree(rest_at);
1285         if (pieces>1)
1286         {
1287             sfree(piecename);
1288         }
1289
1290         membed->it_xy=it_xy;
1291         membed->it_z=it_z;
1292         membed->pos_ins=pos_ins;
1293         membed->r_ins=r_ins;
1294     }
1295
1296     return membed;
1297 }