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