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