2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018, 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);
143 /* Get a list of all the molecule types that are present in a group of atoms. *
144 * Because all interaction within the group to embed are removed on the topology *
145 * level, if the same molecule type is found in another part of the system, these *
146 * would also be affected. Therefore we have to check if the embedded and rest group *
147 * share common molecule types. If so, membed will stop with an error. */
148 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
151 int type = 0, block = 0;
155 snew(tlist->index, at->nr);
156 for (i = 0; i < at->nr; i++)
159 get_mol_id(at->index[i], mtop, &type, &block);
160 for (j = 0; j < nr; j++)
162 if (tlist->index[j] == type)
170 tlist->index[nr] = type;
174 srenew(tlist->index, nr);
178 /* Do the actual check of the molecule types between embedded and rest group */
179 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
181 t_block *ins_mtype, *rest_mtype;
186 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
187 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
189 for (i = 0; i < ins_mtype->nr; i++)
191 for (j = 0; j < rest_mtype->nr; j++)
193 if (ins_mtype->index[i] == rest_mtype->index[j])
195 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
196 "1. Your *.ndx and *.top do not match\n"
197 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
198 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
199 "we need to exclude all interactions between the atoms in the group to\n"
200 "insert, the same moleculetype can not be used in both groups. Change the\n"
201 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
202 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
203 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
208 done_block(ins_mtype);
209 done_block(rest_mtype);
214 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
215 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
216 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
219 std::vector <t_inpfile> inp;
221 wi = init_warning(TRUE, 0);
224 gmx::TextInputFile stream(membed_input);
225 inp = read_inpfile(&stream, membed_input, wi);
228 *it_xy = get_eint(&inp, "nxy", 1000, wi);
229 *it_z = get_eint(&inp, "nz", 0, wi);
230 *xy_fac = get_ereal(&inp, "xyinit", 0.5, wi);
231 *xy_max = get_ereal(&inp, "xyend", 1.0, wi);
232 *z_fac = get_ereal(&inp, "zinit", 1.0, wi);
233 *z_max = get_ereal(&inp, "zend", 1.0, wi);
234 *probe_rad = get_ereal(&inp, "rad", 0.22, wi);
235 *low_up_rm = get_eint(&inp, "ndiff", 0, wi);
236 *maxwarn = get_eint(&inp, "maxwarn", 0, wi);
237 *pieces = get_eint(&inp, "pieces", 1, wi);
238 *bALLOW_ASYMMETRY = get_eeenum(&inp, "asymmetry", yesno_names, wi);
240 check_warning_error(wi, FARGS);
242 gmx::TextOutputFile stream(membed_input);
243 write_inpfile(&stream, membed_input, &inp, FALSE, WriteMdpHeader::yes, wi);
246 done_warning(wi, FARGS);
249 /* Obtain the maximum and minimum coordinates of the group to be embedded */
250 static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins,
251 gmx_groups_t *groups, int ins_grp_id, real xy_max)
254 real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax;
255 const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
256 * determine the overlap between molecule to embed and membrane */
257 const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
258 snew(rest_at->index, state->natoms);
260 xmin = xmax = state->x[ins_at->index[0]][XX];
261 ymin = ymax = state->x[ins_at->index[0]][YY];
262 zmin = zmax = state->x[ins_at->index[0]][ZZ];
264 for (i = 0; i < state->natoms; i++)
266 gid = groups->grpnr[egcFREEZE][i];
267 if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
299 rest_at->index[c] = i;
305 srenew(rest_at->index, c);
307 if (xy_max > fac_inp_size)
309 pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
310 pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
312 pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
313 pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
317 pos_ins->xmin[XX] = xmin;
318 pos_ins->xmin[YY] = ymin;
320 pos_ins->xmax[XX] = xmax;
321 pos_ins->xmax[YY] = ymax;
324 if ( (zmax-zmin) < min_memthick)
326 pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
327 pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
331 pos_ins->xmin[ZZ] = zmin;
332 pos_ins->xmax[ZZ] = zmax;
338 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
339 * xy-plane and counting the number of occupied grid points */
340 static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
342 real x, y, dx = 0.15, dy = 0.15, area = 0.0;
343 real add, memmin, memmax;
346 /* min and max membrane coordinate are altered to reduce the influence of the *
348 memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
349 memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
351 //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
352 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
354 //NOLINTNEXTLINE(clang-analyzer-security.FloatLoopCounter)
355 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
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) )
370 while ( (c < ins_at->nr) && (add < 0.5) );
379 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
381 int i, j, at, mol, nmol, nmolbox, count;
383 real z, zmin, zmax, mem_area;
386 int type = 0, block = 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++)
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]) )
400 mol = get_mol_id(at, mtop, &type, &block);
402 for (j = 0; j < nmol; j++)
404 if (mol == mol_id[j])
432 srenew(mol_id, nmol);
433 mem_p->mol_id = mol_id;
435 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
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);
444 mem_p->zmed = (zmax-zmin)/2+zmin;
446 /*number of membrane molecules in protein box*/
447 nmolbox = count/mtop->moltype[mtop->molblock[block].type].atoms.nr;
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;
453 return mem_p->mem_at.nr;
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)
459 int i, j, at, c, outsidesum, gctr = 0;
463 for (i = 0; i < pos_ins->pieces; i++)
465 idxsum += pos_ins->nidx[i];
468 if (idxsum != ins_at->nr)
470 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
473 snew(pos_ins->geom_cent, pos_ins->pieces);
474 for (i = 0; i < pos_ins->pieces; i++)
478 for (j = 0; j < DIM; j++)
480 pos_ins->geom_cent[i][j] = 0;
483 for (j = 0; j < pos_ins->nidx[i]; j++)
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) )
489 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
501 svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
504 if (!bALLOW_ASYMMETRY)
506 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
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]);
512 fprintf(stderr, "\n");
515 /* resize performed in the md loop */
516 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, const rvec fac)
518 int i, j, k, at, c = 0;
519 for (k = 0; k < pos_ins->pieces; k++)
521 for (i = 0; i < pos_ins->nidx[k]; i++)
523 at = pos_ins->subindex[k][i];
524 for (j = 0; j < DIM; j++)
526 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
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)
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;
548 r_min_rad = probe_rad*probe_rad;
549 gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
550 snew(rm_p->block, molecules.numBlocks());
553 for (i = 0; i < ins_at->nr; i++)
555 at = ins_at->index[i];
556 for (j = 0; j < rest_at->nr; j++)
558 at2 = rest_at->index[j];
559 pbc_dx(pbc, r[at], r[at2], dr);
561 if (norm2(dr) < r_min_rad)
563 mol_id = get_mol_id(at2, mtop, &type, &block);
565 for (l = 0; l < nrm; l++)
567 if (rm_p->mol[l] == mol_id)
575 rm_p->mol[nrm] = mol_id;
576 rm_p->block[nrm] = block;
579 for (l = 0; l < mem_p->nmol; l++)
581 if (mol_id == mem_p->mol_id[l])
583 for (int k : molecules.block(mol_id))
587 z_lip /= molecules.block(mol_id).size();
588 if (z_lip < mem_p->zmed)
603 /*make sure equal number of lipids from upper and lower layer are removed */
604 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
606 snew(dist, mem_p->nmol);
607 snew(order, mem_p->nmol);
608 for (i = 0; i < mem_p->nmol; i++)
610 at = molecules.block(mem_p->mol_id[i]).begin();
611 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
612 if (pos_ins->pieces > 1)
615 min_norm = norm2(dr);
616 for (k = 1; k < pos_ins->pieces; k++)
618 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
619 if (norm2(dr_tmp) < min_norm)
621 min_norm = norm2(dr_tmp);
622 copy_rvec(dr_tmp, dr);
626 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
628 while (j >= 0 && dist[i] < dist[order[j]])
630 order[j+1] = order[j];
637 while (nupper != nlower)
639 mol_id = mem_p->mol_id[order[i]];
640 block = get_molblock(mol_id, mtop->molblock);
642 for (l = 0; l < nrm; l++)
644 if (rm_p->mol[l] == mol_id)
653 for (int k : molecules.block(mol_id))
657 z_lip /= molecules.block(mol_id).size();
658 if (nupper > nlower && z_lip < mem_p->zmed)
660 rm_p->mol[nrm] = mol_id;
661 rm_p->block[nrm] = block;
665 else if (nupper < nlower && z_lip > mem_p->zmed)
667 rm_p->mol[nrm] = mol_id;
668 rm_p->block[nrm] = block;
677 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
685 srenew(rm_p->mol, nrm);
686 srenew(rm_p->block, nrm);
688 return nupper+nlower;
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)
695 int j, k, n, rm, mol_id, at, block;
698 unsigned char *new_egrp[egcNR];
702 /* Construct the molecule range information */
703 gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
705 snew(list, state->natoms);
707 for (int i = 0; i < rm_p->nr; i++)
709 mol_id = rm_p->mol[i];
710 at = molecules.block(mol_id).size();
711 block = rm_p->block[i];
712 mtop->molblock[block].nmol--;
713 for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++)
722 snew(x_tmp, state->natoms);
723 snew(v_tmp, state->natoms);
725 for (int i = 0; i < egcNR; i++)
727 if (groups->grpnr[i] != nullptr)
729 groups->ngrpnr[i] = state->natoms;
730 snew(new_egrp[i], state->natoms);
735 for (int i = 0; i < state->natoms+n; i++)
738 for (j = 0; j < n; j++)
749 for (j = 0; j < egcNR; j++)
751 if (groups->grpnr[j] != nullptr)
753 new_egrp[j][i-rm] = groups->grpnr[j][i];
756 copy_rvec(state->x[i], x_tmp[i-rm]);
757 copy_rvec(state->v[i], v_tmp[i-rm]);
758 for (j = 0; j < ins_at->nr; j++)
760 if (i == ins_at->index[j])
762 ins_at->index[j] = i-rm;
766 for (j = 0; j < pos_ins->pieces; j++)
768 for (k = 0; k < pos_ins->nidx[j]; k++)
770 if (i == pos_ins->subindex[j][k])
772 pos_ins->subindex[j][k] = i-rm;
778 for (int i = 0; i < state->natoms; i++)
780 copy_rvec(x_tmp[i], state->x[i]);
783 for (int i = 0; i < state->natoms; i++)
785 copy_rvec(v_tmp[i], state->v[i]);
789 for (int i = 0; i < egcNR; i++)
791 if (groups->grpnr[i] != nullptr)
793 sfree(groups->grpnr[i]);
794 groups->grpnr[i] = new_egrp[i];
798 /* remove empty molblocks */
800 for (size_t i = 0; i < mtop->molblock.size(); i++)
802 if (mtop->molblock[i].nmol == 0)
808 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
811 mtop->molblock.resize(mtop->molblock.size() - RMmolblock);
814 /* remove al bonded interactions from mtop for the molecule to be embedded */
815 static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
818 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
820 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
821 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
822 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
823 * ins_at index group. MGWolf 050710 */
826 snew(bRM, mtop->moltype.size());
827 for (size_t i = 0; i < mtop->moltype.size(); i++)
832 for (size_t i = 0; i < mtop->molblock.size(); i++)
834 /*loop over molecule blocks*/
835 type = mtop->molblock[i].type;
836 natom = mtop->moltype[type].atoms.nr;
837 nmol = mtop->molblock[i].nmol;
839 for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)
841 /*loop over atoms in the block*/
842 at = j+atom1; /*atom index = block index + offset*/
845 for (m = 0; (m < ins_at->nr) && (bINS == FALSE); m++)
847 /*loop over atoms in insertion index group to determine if we're inserting one*/
848 if (at == ins_at->index[m])
855 atom1 += natom*nmol; /*update offset*/
858 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
862 for (size_t i = 0; i < mtop->moltype.size(); i++)
866 for (j = 0; j < F_LJ; j++)
868 mtop->moltype[i].ilist[j].nr = 0;
871 for (j = F_POSRES; j <= F_VSITEN; j++)
873 mtop->moltype[i].ilist[j].nr = 0;
882 /* Write a topology where the number of molecules is correct for the system after embedding */
883 static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
887 char buf[STRLEN], buf2[STRLEN], *temp;
888 int i, *nmol_rm, nmol, line;
889 char temporary_filename[STRLEN];
891 fpin = gmx_ffopen(topfile, "r");
892 strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
893 gmx_tmpnam(temporary_filename);
894 fpout = gmx_ffopen(temporary_filename, "w");
896 snew(nmol_rm, mtop->moltype.size());
897 for (i = 0; i < rm_p->nr; i++)
899 nmol_rm[rm_p->block[i]]++;
903 while (fgets(buf, STRLEN, fpin))
909 if ((temp = strchr(buf2, '\n')) != nullptr)
917 if ((temp = strchr(buf2, '\n')) != nullptr)
922 if (buf2[strlen(buf2)-1] == ']')
924 buf2[strlen(buf2)-1] = '\0';
927 if (gmx_strcasecmp(buf2, "molecules") == 0)
932 fprintf(fpout, "%s", buf);
934 else if (bMolecules == 1)
936 for (size_t i = 0; i < mtop->molblock.size(); i++)
938 nmol = mtop->molblock[i].nmol;
939 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
940 fprintf(fpout, "%s", buf);
944 else if (bMolecules == 2)
950 fprintf(fpout, "%s", buf);
955 fprintf(fpout, "%s", buf);
960 /* use gmx_ffopen to generate backup of topinout */
961 fpout = gmx_ffopen(topfile, "w");
963 rename(temporary_filename, topfile);
966 void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x)
968 /* Set new positions for the group to embed */
969 if (step_rel <= membed->it_xy)
971 membed->fac[0] += membed->xy_step;
972 membed->fac[1] += membed->xy_step;
974 else if (step_rel <= (membed->it_xy+membed->it_z))
976 membed->fac[2] += membed->z_step;
978 resize(membed->r_ins, x, membed->pos_ins, membed->fac);
981 /* We would like gn to be const as well, but C doesn't allow this */
982 /* TODO this is utility functionality (search for the index of a
983 string in a collection), so should be refactored and located more
984 centrally. Also, it nearly duplicates the same string in readir.c */
985 static int search_string(const char *s, int ng, char *gn[])
989 for (i = 0; (i < ng); i++)
991 if (gmx_strcasecmp(s, gn[i]) == 0)
998 "Group %s selected for embedding was not found in the index file.\n"
999 "Group names must match either [moleculetype] names or custom index group\n"
1000 "names, in which case you must supply an index file to the '-n' option\n"
1007 gmx_membed_t *init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
1008 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
1010 char *ins, **gnames;
1011 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
1012 int ng, j, max_lip_rm, ins_grp_id, ntype, lip_rm;
1014 rvec *r_ins = nullptr;
1015 t_block *ins_at, *rest_at;
1019 gmx_groups_t *groups;
1020 gmx_bool bExcl = FALSE;
1023 char **piecename = nullptr;
1024 gmx_membed_t *membed = nullptr;
1026 /* input variables */
1033 real probe_rad = 0.22;
1037 gmx_bool bALLOW_ASYMMETRY = FALSE;
1039 /* sanity check constants */ /* Issue a warning when: */
1040 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1041 * and rest smaller than this value is probably too small */
1042 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1043 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1044 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1045 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1046 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1054 /* get input data out membed file */
1057 get_input(opt2fn("-membed", nfile, fnm),
1058 &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1059 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1061 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1063 if (!EI_DYNAMICS(inputrec->eI) )
1065 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1070 gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1075 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1078 groups = &(mtop->groups);
1079 snew(gnames, groups->ngrpname);
1080 for (i = 0; i < groups->ngrpname; i++)
1082 gnames[i] = *(groups->grpname[i]);
1085 atoms = gmx_mtop_global_atoms(mtop);
1087 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1088 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1089 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1090 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1091 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1093 pos_ins->pieces = pieces;
1094 snew(pos_ins->nidx, pieces);
1095 snew(pos_ins->subindex, pieces);
1096 snew(piecename, pieces);
1099 fprintf(stderr, "\nSelect pieces to embed:\n");
1100 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1104 /*use whole embedded group*/
1105 snew(pos_ins->nidx, 1);
1106 snew(pos_ins->subindex, 1);
1107 pos_ins->nidx[0] = ins_at->nr;
1108 pos_ins->subindex[0] = ins_at->index;
1111 if (probe_rad < min_probe_rad)
1114 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1115 "in overlap between waters and the group to embed, which will result "
1116 "in Lincs errors etc.\n\n", warn);
1119 if (xy_fac < min_xy_init)
1122 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too small.\n\n", warn, ins);
1125 if (it_xy < min_it_xy)
1128 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1129 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1132 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1135 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1136 " is probably too small.\nIncrease -nz or the maxwarn setting in the membed input file.\n\n", warn, ins, it_z);
1139 if (it_xy+it_z > inputrec->nsteps)
1142 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1143 "number of steps in the tpr.\n"
1144 "(increase maxwarn in the membed input file to override)\n\n", warn);
1148 if (inputrec->opts.ngfrz == 1)
1150 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1153 for (i = 0; i < inputrec->opts.ngfrz; i++)
1155 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1156 if (ins_grp_id == tmp_id)
1165 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1168 for (i = 0; i < DIM; i++)
1170 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1172 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1176 ng = groups->grps[egcENER].nr;
1179 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1182 for (i = 0; i < ng; i++)
1184 for (j = 0; j < ng; j++)
1186 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1189 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1190 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1192 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1193 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1194 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1202 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1203 "the freeze group");
1206 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1208 init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1209 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1210 check_types(ins_at, rest_at, mtop);
1212 init_mem_at(mem_p, mtop, as_rvec_array(state->x.data()), state->box, pos_ins);
1214 prot_area = est_prot_area(pos_ins, as_rvec_array(state->x.data()), ins_at, mem_p);
1215 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) )
1218 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1219 "This might cause pressure problems during the growth phase. Just try with\n"
1220 "current setup and increase 'maxwarn' in your membed settings file, but lower the\n"
1221 "compressibility in the mdp-file or disable pressure coupling if problems occur.\n\n", warn);
1226 gmx_fatal(FARGS, "Too many warnings (override by setting maxwarn in the membed input file)\n");
1229 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1230 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1231 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1233 /* Maximum number of lipids to be removed*/
1234 max_lip_rm = (int)(2*prot_area/mem_p->lip_area);
1235 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1237 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1238 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1239 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1240 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1242 /* resize the protein by xy and by z if necessary*/
1243 snew(r_ins, ins_at->nr);
1244 init_resize(ins_at, r_ins, pos_ins, mem_p, as_rvec_array(state->x.data()), bALLOW_ASYMMETRY);
1245 membed->fac[0] = membed->fac[1] = xy_fac;
1246 membed->fac[2] = z_fac;
1248 membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
1249 membed->z_step = (z_max-z_fac)/(double)(it_z-1);
1251 resize(r_ins, as_rvec_array(state->x.data()), pos_ins, membed->fac);
1253 /* remove overlapping lipids and water from the membrane box*/
1254 /*mark molecules to be removed*/
1256 set_pbc(pbc, inputrec->ePBC, state->box);
1259 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, as_rvec_array(state->x.data()), mem_p, pos_ins,
1260 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1261 lip_rm -= low_up_rm;
1265 for (i = 0; i < rm_p->nr; i++)
1267 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1271 for (size_t i = 0; i < mtop->molblock.size(); i++)
1274 for (j = 0; j < rm_p->nr; j++)
1276 if (rm_p->block[j] == static_cast<int>(i))
1281 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1284 if (lip_rm > max_lip_rm)
1287 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1288 "protein area\nTry making the -xyinit resize factor smaller or increase "
1289 "maxwarn in the membed input file.\n\n", warn);
1292 /*remove all lipids and waters overlapping and update all important structures*/
1293 rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
1295 rm_bonded_at = rm_bonded(ins_at, mtop);
1296 if (rm_bonded_at != ins_at->nr)
1298 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1299 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1300 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1305 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless,\n"
1306 "you can increase the maxwarn setting in the membed input file.");
1309 // Re-establish the invariants of the derived values within
1311 gmx_mtop_finalize(mtop);
1313 if (ftp2bSet(efTOP, nfile, fnm))
1315 top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
1325 membed->it_xy = it_xy;
1326 membed->it_z = it_z;
1327 membed->pos_ins = pos_ins;
1328 membed->r_ins = r_ins;
1334 void free_membed(gmx_membed_t *membed)