src/programs/mdrun/membed.c

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