src/programs/mdrun/membed.c

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