src/programs/mdrun/membed.c

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