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