2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018,2019, 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.
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.
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.
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.
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.
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.
42 #include "gromacs/commandline/filenm.h"
43 #include "gromacs/fileio/readinp.h"
44 #include "gromacs/fileio/warninp.h"
45 #include "gromacs/gmxlib/network.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/mdtypes/commrec.h"
48 #include "gromacs/mdtypes/inputrec.h"
49 #include "gromacs/mdtypes/md_enums.h"
50 #include "gromacs/mdtypes/state.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/topology/index.h"
53 #include "gromacs/topology/mtop_lookup.h"
54 #include "gromacs/topology/mtop_util.h"
55 #include "gromacs/topology/topology.h"
56 #include "gromacs/utility/cstringutil.h"
57 #include "gromacs/utility/exceptions.h"
58 #include "gromacs/utility/fatalerror.h"
59 #include "gromacs/utility/filestream.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/smalloc.h"
63 /* information about scaling center */
65 rvec xmin; /* smallest coordinates of all embedded molecules */
66 rvec xmax; /* largest coordinates of all embedded molecules */
67 rvec *geom_cent; /* scaling center of each independent molecule to embed */
68 int pieces; /* number of molecules to embed independently */
69 int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
70 int **subindex; /* atomids for independent molecule *
71 * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
74 /* variables needed in do_md */
76 int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
77 int it_z; /* same, but for z */
78 real xy_step; /* stepsize used during resize in xy-plane */
79 real z_step; /* same, but in z */
80 rvec fac; /* initial scaling of the molecule to grow into the membrane */
81 rvec *r_ins; /* final coordinates of the molecule to grow */
82 pos_ins_t *pos_ins; /* scaling center for each piece to embed */
85 /* membrane related variables */
87 char *name; /* name of index group to embed molecule into (usually membrane) */
88 t_block mem_at; /* list all atoms in membrane */
89 int nmol; /* number of membrane molecules overlapping with the molecule to embed */
90 int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
91 real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
92 real zmin; /* minimum z coordinate of membrane */
93 real zmax; /* maximum z coordinate of membrane */
94 real zmed; /* median z coordinate of membrane */
97 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
98 * These will then be removed from the system */
100 int nr; /* number of molecules to remove */
101 int *mol; /* list of molecule ids to remove */
102 int *block; /* id of the molblock that the molecule to remove is part of */
105 /* Get the global molecule id, and the corresponding molecule type and id of the *
106 * molblock from the global atom nr. */
107 static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
114 mtopGetMolblockIndex(mtop, at, block, &mol_id, &atnr_mol);
115 for (i = 0; i < *block; i++)
117 mol_id += mtop->molblock[i].nmol;
119 *type = mtop->molblock[*block].type;
124 /* Get the id of the molblock from a global molecule id */
125 static int get_molblock(int mol_id, const std::vector<gmx_molblock_t> &mblock)
129 for (size_t i = 0; i < mblock.size(); i++)
131 nmol += mblock[i].nmol;
138 gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
141 /* Get a list of all the molecule types that are present in a group of atoms. *
142 * Because all interaction within the group to embed are removed on the topology *
143 * level, if the same molecule type is found in another part of the system, these *
144 * would also be affected. Therefore we have to check if the embedded and rest group *
145 * share common molecule types. If so, membed will stop with an error. */
146 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
149 int type = 0, block = 0;
153 snew(tlist->index, at->nr);
154 for (i = 0; i < at->nr; i++)
157 get_mol_id(at->index[i], mtop, &type, &block);
158 for (j = 0; j < nr; j++)
160 if (tlist->index[j] == type)
168 tlist->index[nr] = type;
172 srenew(tlist->index, nr);
176 /* Do the actual check of the molecule types between embedded and rest group */
177 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
179 t_block *ins_mtype, *rest_mtype;
184 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
185 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
187 for (i = 0; i < ins_mtype->nr; i++)
189 for (j = 0; j < rest_mtype->nr; j++)
191 if (ins_mtype->index[i] == rest_mtype->index[j])
193 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
194 "1. Your *.ndx and *.top do not match\n"
195 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
196 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
197 "we need to exclude all interactions between the atoms in the group to\n"
198 "insert, the same moleculetype can not be used in both groups. Change the\n"
199 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
200 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
201 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
206 done_block(ins_mtype);
207 done_block(rest_mtype);
212 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
213 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
214 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
217 std::vector <t_inpfile> inp;
219 wi = init_warning(TRUE, 0);
222 gmx::TextInputFile stream(membed_input);
223 inp = read_inpfile(&stream, membed_input, wi);
226 *it_xy = get_eint(&inp, "nxy", 1000, wi);
227 *it_z = get_eint(&inp, "nz", 0, wi);
228 *xy_fac = get_ereal(&inp, "xyinit", 0.5, wi);
229 *xy_max = get_ereal(&inp, "xyend", 1.0, wi);
230 *z_fac = get_ereal(&inp, "zinit", 1.0, wi);
231 *z_max = get_ereal(&inp, "zend", 1.0, wi);
232 *probe_rad = get_ereal(&inp, "rad", 0.22, wi);
233 *low_up_rm = get_eint(&inp, "ndiff", 0, wi);
234 *maxwarn = get_eint(&inp, "maxwarn", 0, wi);
235 *pieces = get_eint(&inp, "pieces", 1, wi);
236 const char *yesno_names[BOOL_NR+1] =
240 *bALLOW_ASYMMETRY = (get_eeenum(&inp, "asymmetry", yesno_names, wi) != 0);
242 check_warning_error(wi, FARGS);
244 gmx::TextOutputFile stream(membed_input);
245 write_inpfile(&stream, membed_input, &inp, FALSE, WriteMdpHeader::yes, wi);
248 done_warning(wi, FARGS);
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)
256 real xmin, xmax, ymin, ymax, 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 auto x = makeArrayRef(state->x);
263 xmin = xmax = x[ins_at->index[0]][XX];
264 ymin = ymax = x[ins_at->index[0]][YY];
265 zmin = zmax = x[ins_at->index[0]][ZZ];
267 for (i = 0; i < state->natoms; i++)
269 gid = groups->grpnr[egcFREEZE][i];
270 if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
272 xmin = std::min(xmin, x[i][XX]);
273 xmax = std::max(xmax, x[i][XX]);
274 ymin = std::min(ymin, x[i][YY]);
275 ymax = std::max(ymax, x[i][YY]);
276 zmin = std::min(zmin, x[i][ZZ]);
277 zmax = std::max(zmax, x[i][ZZ]);
281 rest_at->index[c] = i;
287 srenew(rest_at->index, c);
289 if (xy_max > fac_inp_size)
291 pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
292 pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
294 pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
295 pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
299 pos_ins->xmin[XX] = xmin;
300 pos_ins->xmin[YY] = ymin;
302 pos_ins->xmax[XX] = xmax;
303 pos_ins->xmax[YY] = ymax;
306 if ( (zmax-zmin) < min_memthick)
308 pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
309 pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
313 pos_ins->xmin[ZZ] = zmin;
314 pos_ins->xmax[ZZ] = zmax;
320 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
321 * xy-plane and counting the number of occupied grid points */
322 static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
324 real x, y, dx = 0.15, dy = 0.15, area = 0.0;
325 real add, memmin, memmax;
328 /* min and max membrane coordinate are altered to reduce the influence of the *
330 memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
331 memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
333 //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
334 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
336 //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
337 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
343 at = ins_at->index[c];
344 if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
345 (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
346 (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
352 while ( (c < ins_at->nr) && (add < 0.5) );
361 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
363 int i, j, at, mol, nmol, nmolbox, count;
365 real z, zmin, zmax, mem_area;
368 int type = 0, block = 0;
371 mem_a = &(mem_p->mem_at);
372 snew(mol_id, mem_a->nr);
373 zmin = pos_ins->xmax[ZZ];
374 zmax = pos_ins->xmin[ZZ];
375 for (i = 0; i < mem_a->nr; i++)
377 at = mem_a->index[i];
378 if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
379 (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
380 (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
382 mol = get_mol_id(at, mtop, &type, &block);
384 for (j = 0; j < nmol; j++)
386 if (mol == mol_id[j])
414 srenew(mol_id, nmol);
415 mem_p->mol_id = mol_id;
417 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
419 gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
420 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
421 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
422 "your water layer is not thick enough.\n", zmax, zmin);
426 mem_p->zmed = (zmax-zmin)/2+zmin;
428 /*number of membrane molecules in protein box*/
429 nmolbox = count/mtop->moltype[mtop->molblock[block].type].atoms.nr;
430 /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
431 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
432 /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
433 mem_p->lip_area = 2.0*mem_area/static_cast<double>(nmolbox);
435 return mem_p->mem_at.nr;
438 static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
439 gmx_bool bALLOW_ASYMMETRY)
441 int i, j, at, c, outsidesum, gctr = 0;
445 for (i = 0; i < pos_ins->pieces; i++)
447 idxsum += pos_ins->nidx[i];
450 if (idxsum != ins_at->nr)
452 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
455 snew(pos_ins->geom_cent, pos_ins->pieces);
456 for (i = 0; i < pos_ins->pieces; i++)
460 for (j = 0; j < DIM; j++)
462 pos_ins->geom_cent[i][j] = 0;
465 for (j = 0; j < pos_ins->nidx[i]; j++)
467 at = pos_ins->subindex[i][j];
468 copy_rvec(r[at], r_ins[gctr]);
469 if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
471 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
483 svmul(1/static_cast<double>(c), pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
486 if (!bALLOW_ASYMMETRY)
488 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
491 fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
492 i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
494 fprintf(stderr, "\n");
497 /* resize performed in the md loop */
498 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, const rvec fac)
500 int i, j, k, at, c = 0;
501 for (k = 0; k < pos_ins->pieces; k++)
503 for (i = 0; i < pos_ins->nidx[k]; i++)
505 at = pos_ins->subindex[k][i];
506 for (j = 0; j < DIM; j++)
508 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
515 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
516 * The molecule to be embedded is already reduced in size. */
517 static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
518 rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
519 int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
521 int i, j, k, l, at, at2, mol_id;
522 int type = 0, block = 0;
523 int nrm, nupper, nlower;
524 real r_min_rad, z_lip, min_norm;
530 r_min_rad = probe_rad*probe_rad;
531 gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
532 snew(rm_p->block, molecules.numBlocks());
533 snew(rm_p->mol, molecules.numBlocks());
536 for (i = 0; i < ins_at->nr; i++)
538 at = ins_at->index[i];
539 for (j = 0; j < rest_at->nr; j++)
541 at2 = rest_at->index[j];
542 pbc_dx(pbc, r[at], r[at2], dr);
544 if (norm2(dr) < r_min_rad)
546 mol_id = get_mol_id(at2, mtop, &type, &block);
548 for (l = 0; l < nrm; l++)
550 if (rm_p->mol[l] == mol_id)
558 rm_p->mol[nrm] = mol_id;
559 rm_p->block[nrm] = block;
562 for (l = 0; l < mem_p->nmol; l++)
564 if (mol_id == mem_p->mol_id[l])
566 for (int k : molecules.block(mol_id))
570 z_lip /= molecules.block(mol_id).size();
571 if (z_lip < mem_p->zmed)
586 /*make sure equal number of lipids from upper and lower layer are removed */
587 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
589 snew(dist, mem_p->nmol);
590 snew(order, mem_p->nmol);
591 for (i = 0; i < mem_p->nmol; i++)
593 at = molecules.block(mem_p->mol_id[i]).begin();
594 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
595 if (pos_ins->pieces > 1)
598 min_norm = norm2(dr);
599 for (k = 1; k < pos_ins->pieces; k++)
601 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
602 if (norm2(dr_tmp) < min_norm)
604 min_norm = norm2(dr_tmp);
605 copy_rvec(dr_tmp, dr);
609 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
611 while (j >= 0 && dist[i] < dist[order[j]])
613 order[j+1] = order[j];
620 while (nupper != nlower)
622 mol_id = mem_p->mol_id[order[i]];
623 block = get_molblock(mol_id, mtop->molblock);
625 for (l = 0; l < nrm; l++)
627 if (rm_p->mol[l] == mol_id)
636 for (int k : molecules.block(mol_id))
640 z_lip /= molecules.block(mol_id).size();
641 if (nupper > nlower && z_lip < mem_p->zmed)
643 rm_p->mol[nrm] = mol_id;
644 rm_p->block[nrm] = block;
648 else if (nupper < nlower && z_lip > mem_p->zmed)
650 rm_p->mol[nrm] = mol_id;
651 rm_p->block[nrm] = block;
660 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
668 srenew(rm_p->mol, nrm);
669 srenew(rm_p->block, nrm);
671 return nupper+nlower;
674 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
675 static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
676 t_block *ins_at, pos_ins_t *pos_ins)
678 int j, k, n, rm, mol_id, at, block;
681 unsigned char *new_egrp[egcNR];
685 /* Construct the molecule range information */
686 gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
688 snew(list, state->natoms);
690 for (int i = 0; i < rm_p->nr; i++)
692 mol_id = rm_p->mol[i];
693 at = molecules.block(mol_id).begin();
694 block = rm_p->block[i];
695 mtop->molblock[block].nmol--;
696 for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++)
704 /* We cannot change the size of the state datastructures here
705 * because we still access the coordinate arrays for all positions
706 * before removing the molecules we want to remove.
708 const int newStateAtomNumber = state->natoms - n;
709 snew(x_tmp, newStateAtomNumber);
710 snew(v_tmp, newStateAtomNumber);
712 for (int i = 0; i < egcNR; i++)
714 if (groups->grpnr[i] != nullptr)
716 groups->ngrpnr[i] = newStateAtomNumber;
717 snew(new_egrp[i], newStateAtomNumber);
721 auto x = makeArrayRef(state->x);
722 auto v = makeArrayRef(state->v);
724 for (int i = 0; i < state->natoms; i++)
727 for (j = 0; j < n; j++)
738 for (j = 0; j < egcNR; j++)
740 if (groups->grpnr[j] != nullptr)
742 new_egrp[j][i-rm] = groups->grpnr[j][i];
745 copy_rvec(x[i], x_tmp[i-rm]);
746 copy_rvec(v[i], v_tmp[i-rm]);
747 for (j = 0; j < ins_at->nr; j++)
749 if (i == ins_at->index[j])
751 ins_at->index[j] = i-rm;
755 for (j = 0; j < pos_ins->pieces; j++)
757 for (k = 0; k < pos_ins->nidx[j]; k++)
759 if (i == pos_ins->subindex[j][k])
761 pos_ins->subindex[j][k] = i-rm;
767 state_change_natoms(state, newStateAtomNumber);
768 for (int i = 0; i < state->natoms; i++)
770 copy_rvec(x_tmp[i], x[i]);
773 for (int i = 0; i < state->natoms; i++)
775 copy_rvec(v_tmp[i], v[i]);
779 for (int i = 0; i < egcNR; i++)
781 if (groups->grpnr[i] != nullptr)
783 sfree(groups->grpnr[i]);
784 groups->grpnr[i] = new_egrp[i];
788 /* remove empty molblocks */
790 for (size_t i = 0; i < mtop->molblock.size(); i++)
792 if (mtop->molblock[i].nmol == 0)
798 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
801 mtop->molblock.resize(mtop->molblock.size() - RMmolblock);
804 /* remove al bonded interactions from mtop for the molecule to be embedded */
805 static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
808 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
810 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
811 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
812 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
813 * ins_at index group. MGWolf 050710 */
816 snew(bRM, mtop->moltype.size());
817 for (size_t i = 0; i < mtop->moltype.size(); i++)
822 for (size_t i = 0; i < mtop->molblock.size(); i++)
824 /*loop over molecule blocks*/
825 type = mtop->molblock[i].type;
826 natom = mtop->moltype[type].atoms.nr;
827 nmol = mtop->molblock[i].nmol;
829 for (j = 0; j < natom*nmol && bRM[type]; j++)
831 /*loop over atoms in the block*/
832 at = j+atom1; /*atom index = block index + offset*/
835 for (m = 0; (m < ins_at->nr) && (!bINS); m++)
837 /*loop over atoms in insertion index group to determine if we're inserting one*/
838 if (at == ins_at->index[m])
845 atom1 += natom*nmol; /*update offset*/
848 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
852 for (size_t i = 0; i < mtop->moltype.size(); i++)
856 for (j = 0; j < F_LJ; j++)
858 mtop->moltype[i].ilist[j].iatoms.clear();
861 for (j = F_POSRES; j <= F_VSITEN; j++)
863 mtop->moltype[i].ilist[j].iatoms.clear();
872 /* Write a topology where the number of molecules is correct for the system after embedding */
873 static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
877 char buf[STRLEN], buf2[STRLEN], *temp;
878 int i, *nmol_rm, nmol, line;
879 char temporary_filename[STRLEN];
881 fpin = gmx_ffopen(topfile, "r");
882 strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
883 gmx_tmpnam(temporary_filename);
884 fpout = gmx_ffopen(temporary_filename, "w");
886 snew(nmol_rm, mtop->moltype.size());
887 for (i = 0; i < rm_p->nr; i++)
889 nmol_rm[rm_p->block[i]]++;
893 while (fgets(buf, STRLEN, fpin))
899 if ((temp = strchr(buf2, '\n')) != nullptr)
907 if ((temp = strchr(buf2, '\n')) != nullptr)
912 if (buf2[strlen(buf2)-1] == ']')
914 buf2[strlen(buf2)-1] = '\0';
917 if (gmx_strcasecmp(buf2, "molecules") == 0)
922 fprintf(fpout, "%s", buf);
924 else if (bMolecules == 1)
926 for (size_t i = 0; i < mtop->molblock.size(); i++)
928 nmol = mtop->molblock[i].nmol;
929 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
930 fprintf(fpout, "%s", buf);
934 else if (bMolecules == 2)
940 fprintf(fpout, "%s", buf);
945 fprintf(fpout, "%s", buf);
950 /* use gmx_ffopen to generate backup of topinout */
951 fpout = gmx_ffopen(topfile, "w");
953 rename(temporary_filename, topfile);
956 void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x)
958 /* Set new positions for the group to embed */
959 if (step_rel <= membed->it_xy)
961 membed->fac[0] += membed->xy_step;
962 membed->fac[1] += membed->xy_step;
964 else if (step_rel <= (membed->it_xy+membed->it_z))
966 membed->fac[2] += membed->z_step;
968 resize(membed->r_ins, x, membed->pos_ins, membed->fac);
971 /* We would like gn to be const as well, but C doesn't allow this */
972 /* TODO this is utility functionality (search for the index of a
973 string in a collection), so should be refactored and located more
974 centrally. Also, it nearly duplicates the same string in readir.c */
975 static int search_string(const char *s, int ng, char *gn[])
979 for (i = 0; (i < ng); i++)
981 if (gmx_strcasecmp(s, gn[i]) == 0)
988 "Group %s selected for embedding was not found in the index file.\n"
989 "Group names must match either [moleculetype] names or custom index group\n"
990 "names, in which case you must supply an index file to the '-n' option\n"
995 gmx_membed_t *init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
996 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
999 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
1000 int ng, j, max_lip_rm, ins_grp_id, ntype, lip_rm;
1002 rvec *r_ins = nullptr;
1003 t_block *ins_at, *rest_at;
1007 gmx_groups_t *groups;
1008 gmx_bool bExcl = FALSE;
1011 char **piecename = nullptr;
1012 gmx_membed_t *membed = nullptr;
1014 /* input variables */
1021 real probe_rad = 0.22;
1025 gmx_bool bALLOW_ASYMMETRY = FALSE;
1027 /* sanity check constants */ /* Issue a warning when: */
1028 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1029 * and rest smaller than this value is probably too small */
1030 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1031 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1032 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1033 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1034 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1042 fprintf(fplog, "Note: it is expected that in future gmx mdrun -membed will not be the "
1043 "way to access this feature, perhaps changing to e.g. gmx membed.");
1044 /* get input data out membed file */
1047 get_input(opt2fn("-membed", nfile, fnm),
1048 &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1049 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1051 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1053 if (!EI_DYNAMICS(inputrec->eI) )
1055 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1060 gmx_input("Sorry, parallel membed is not yet fully functional.");
1065 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1068 groups = &(mtop->groups);
1069 snew(gnames, groups->ngrpname);
1070 for (i = 0; i < groups->ngrpname; i++)
1072 gnames[i] = *(groups->grpname[i]);
1075 atoms = gmx_mtop_global_atoms(mtop);
1077 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1078 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1079 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1080 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1081 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1083 pos_ins->pieces = pieces;
1084 snew(pos_ins->nidx, pieces);
1085 snew(pos_ins->subindex, pieces);
1086 snew(piecename, pieces);
1089 fprintf(stderr, "\nSelect pieces to embed:\n");
1090 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1094 /*use whole embedded group*/
1095 snew(pos_ins->nidx, 1);
1096 snew(pos_ins->subindex, 1);
1097 pos_ins->nidx[0] = ins_at->nr;
1098 pos_ins->subindex[0] = ins_at->index;
1101 if (probe_rad < min_probe_rad)
1104 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1105 "in overlap between waters and the group to embed, which will result "
1106 "in Lincs errors etc.\n\n", warn);
1109 if (xy_fac < min_xy_init)
1112 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too small.\n\n", warn, ins);
1115 if (it_xy < min_it_xy)
1118 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1119 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1122 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1125 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1126 " is probably too small.\nIncrease -nz or the maxwarn setting in the membed input file.\n\n", warn, ins, it_z);
1129 if (it_xy+it_z > inputrec->nsteps)
1132 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1133 "number of steps in the tpr.\n"
1134 "(increase maxwarn in the membed input file to override)\n\n", warn);
1138 if (inputrec->opts.ngfrz == 1)
1140 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1143 for (i = 0; i < inputrec->opts.ngfrz; i++)
1145 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1146 if (ins_grp_id == tmp_id)
1155 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1158 for (i = 0; i < DIM; i++)
1160 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1162 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1166 ng = groups->grps[egcENER].nr;
1169 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1172 for (i = 0; i < ng; i++)
1174 for (j = 0; j < ng; j++)
1176 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1179 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1180 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1182 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1183 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1184 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1192 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1193 "the freeze group");
1196 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1198 init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1199 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1200 check_types(ins_at, rest_at, mtop);
1202 init_mem_at(mem_p, mtop, state->x.rvec_array(), state->box, pos_ins);
1204 prot_area = est_prot_area(pos_ins, state->x.rvec_array(), ins_at, mem_p);
1205 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) )
1208 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1209 "This might cause pressure problems during the growth phase. Just try with\n"
1210 "current setup and increase 'maxwarn' in your membed settings file, but lower the\n"
1211 "compressibility in the mdp-file or disable pressure coupling if problems occur.\n\n", warn);
1216 gmx_fatal(FARGS, "Too many warnings (override by setting maxwarn in the membed input file)\n");
1219 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1220 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1221 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1223 /* Maximum number of lipids to be removed*/
1224 max_lip_rm = static_cast<int>(2*prot_area/mem_p->lip_area);
1225 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1227 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1228 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1229 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1230 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1232 /* resize the protein by xy and by z if necessary*/
1233 snew(r_ins, ins_at->nr);
1234 init_resize(ins_at, r_ins, pos_ins, mem_p, state->x.rvec_array(), bALLOW_ASYMMETRY);
1235 membed->fac[0] = membed->fac[1] = xy_fac;
1236 membed->fac[2] = z_fac;
1238 membed->xy_step = (xy_max-xy_fac)/static_cast<double>(it_xy);
1239 membed->z_step = (z_max-z_fac)/static_cast<double>(it_z-1);
1241 resize(r_ins, state->x.rvec_array(), pos_ins, membed->fac);
1243 /* remove overlapping lipids and water from the membrane box*/
1244 /*mark molecules to be removed*/
1246 set_pbc(pbc, inputrec->ePBC, state->box);
1249 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x.rvec_array(), mem_p, pos_ins,
1250 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1251 lip_rm -= low_up_rm;
1255 for (i = 0; i < rm_p->nr; i++)
1257 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1261 for (size_t i = 0; i < mtop->molblock.size(); i++)
1264 for (j = 0; j < rm_p->nr; j++)
1266 if (rm_p->block[j] == static_cast<int>(i))
1271 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1274 if (lip_rm > max_lip_rm)
1277 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1278 "protein area\nTry making the -xyinit resize factor smaller or increase "
1279 "maxwarn in the membed input file.\n\n", warn);
1282 /*remove all lipids and waters overlapping and update all important structures*/
1283 rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
1285 rm_bonded_at = rm_bonded(ins_at, mtop);
1286 if (rm_bonded_at != ins_at->nr)
1288 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1289 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1290 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1295 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless,\n"
1296 "you can increase the maxwarn setting in the membed input file.");
1299 // Re-establish the invariants of the derived values within
1301 gmx_mtop_finalize(mtop);
1303 if (ftp2bSet(efTOP, nfile, fnm))
1305 top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
1315 membed->it_xy = it_xy;
1316 membed->it_z = it_z;
1317 membed->pos_ins = pos_ins;
1318 membed->r_ins = r_ins;
1324 void free_membed(gmx_membed_t *membed)