2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014,2015, 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/essentialdynamics/edsam.h"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/legacyheaders/macros.h"
45 #include "gromacs/legacyheaders/names.h"
46 #include "gromacs/legacyheaders/readinp.h"
47 #include "gromacs/legacyheaders/typedefs.h"
48 #include "gromacs/legacyheaders/types/commrec.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/pbcutil/pbc.h"
51 #include "gromacs/topology/index.h"
52 #include "gromacs/topology/mtop_util.h"
53 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/smalloc.h"
57 /* information about scaling center */
59 rvec xmin; /* smallest coordinates of all embedded molecules */
60 rvec xmax; /* largest coordinates of all embedded molecules */
61 rvec *geom_cent; /* scaling center of each independent molecule to embed */
62 int pieces; /* number of molecules to embed independently */
63 int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
64 atom_id **subindex; /* atomids for independent molecule *
65 * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
68 /* variables needed in do_md */
70 int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
71 int it_z; /* same, but for z */
72 real xy_step; /* stepsize used during resize in xy-plane */
73 real z_step; /* same, but in z */
74 rvec fac; /* initial scaling of the molecule to grow into the membrane */
75 rvec *r_ins; /* final coordinates of the molecule to grow */
76 pos_ins_t *pos_ins; /* scaling center for each piece to embed */
79 /* membrane related variables */
81 char *name; /* name of index group to embed molecule into (usually membrane) */
82 t_block mem_at; /* list all atoms in membrane */
83 int nmol; /* number of membrane molecules overlapping with the molecule to embed */
84 int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
85 real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
86 real zmin; /* minimum z coordinate of membrane */
87 real zmax; /* maximum z coordinate of membrane */
88 real zmed; /* median z coordinate of membrane */
91 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
92 * These will then be removed from the system */
94 int nr; /* number of molecules to remove */
95 int *mol; /* list of molecule ids to remove */
96 int *block; /* id of the molblock that the molecule to remove is part of */
99 /* Get the global molecule id, and the corresponding molecule type and id of the *
100 * molblock from the global atom nr. */
101 static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
106 gmx_mtop_atomlookup_t alook;
108 alook = gmx_mtop_atomlookup_settle_init(mtop);
109 gmx_mtop_atomnr_to_molblock_ind(alook, at, block, &mol_id, &atnr_mol);
110 for (i = 0; i < *block; i++)
112 mol_id += mtop->molblock[i].nmol;
114 *type = mtop->molblock[*block].type;
116 gmx_mtop_atomlookup_destroy(alook);
121 /* Get the id of the molblock from a global molecule id */
122 static int get_molblock(int mol_id, int nmblock, gmx_molblock_t *mblock)
127 for (i = 0; i < nmblock; i++)
129 nmol += mblock[i].nmol;
136 gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
141 static int get_tpr_version(const char *infile)
144 int version, generation;
146 read_tpxheader(infile, &header, TRUE, &version, &generation);
151 /* Get a list of all the molecule types that are present in a group of atoms. *
152 * Because all interaction within the group to embed are removed on the topology *
153 * level, if the same molecule type is found in another part of the system, these *
154 * would also be affected. Therefore we have to check if the embedded and rest group *
155 * share common molecule types. If so, membed will stop with an error. */
156 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
158 int i, j, nr, mol_id;
159 int type = 0, block = 0;
163 snew(tlist->index, at->nr);
164 for (i = 0; i < at->nr; i++)
167 mol_id = get_mol_id(at->index[i], mtop, &type, &block);
168 for (j = 0; j < nr; j++)
170 if (tlist->index[j] == type)
178 tlist->index[nr] = type;
182 srenew(tlist->index, nr);
186 /* Do the actual check of the molecule types between embedded and rest group */
187 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
189 t_block *ins_mtype, *rest_mtype;
194 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
195 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
197 for (i = 0; i < ins_mtype->nr; i++)
199 for (j = 0; j < rest_mtype->nr; j++)
201 if (ins_mtype->index[i] == rest_mtype->index[j])
203 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
204 "1. Your *.ndx and *.top do not match\n"
205 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
206 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
207 "we need to exclude all interactions between the atoms in the group to\n"
208 "insert, the same moleculetype can not be used in both groups. Change the\n"
209 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
210 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
211 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
216 sfree(ins_mtype->index);
217 sfree(rest_mtype->index);
222 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
223 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
224 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
230 wi = init_warning(TRUE, 0);
232 inp = read_inpfile(membed_input, &ninp, wi);
233 ITYPE ("nxy", *it_xy, 1000);
234 ITYPE ("nz", *it_z, 0);
235 RTYPE ("xyinit", *xy_fac, 0.5);
236 RTYPE ("xyend", *xy_max, 1.0);
237 RTYPE ("zinit", *z_fac, 1.0);
238 RTYPE ("zend", *z_max, 1.0);
239 RTYPE ("rad", *probe_rad, 0.22);
240 ITYPE ("ndiff", *low_up_rm, 0);
241 ITYPE ("maxwarn", *maxwarn, 0);
242 ITYPE ("pieces", *pieces, 1);
243 EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
245 write_inpfile(membed_input, ninp, inp, FALSE, wi);
248 /* Obtain the maximum and minimum coordinates of the group to be embedded */
249 static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins,
250 gmx_groups_t *groups, int ins_grp_id, real xy_max)
253 real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax;
254 const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
255 * determine the overlap between molecule to embed and membrane */
256 const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
257 snew(rest_at->index, state->natoms);
259 xmin = xmax = state->x[ins_at->index[0]][XX];
260 ymin = ymax = state->x[ins_at->index[0]][YY];
261 zmin = zmax = state->x[ins_at->index[0]][ZZ];
263 for (i = 0; i < state->natoms; i++)
265 gid = groups->grpnr[egcFREEZE][i];
266 if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
298 rest_at->index[c] = i;
304 srenew(rest_at->index, c);
306 if (xy_max > fac_inp_size)
308 pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
309 pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
311 pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
312 pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
316 pos_ins->xmin[XX] = xmin;
317 pos_ins->xmin[YY] = ymin;
319 pos_ins->xmax[XX] = xmax;
320 pos_ins->xmax[YY] = ymax;
323 if ( (zmax-zmin) < min_memthick)
325 pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
326 pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
330 pos_ins->xmin[ZZ] = zmin;
331 pos_ins->xmax[ZZ] = zmax;
337 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
338 * xy-plane and counting the number of occupied grid points */
339 static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
341 real x, y, dx = 0.15, dy = 0.15, area = 0.0;
342 real add, memmin, memmax;
345 /* min and max membrane coordinate are altered to reduce the influence of the *
347 memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
348 memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
350 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
352 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
358 at = ins_at->index[c];
359 if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
360 (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
361 (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
367 while ( (c < ins_at->nr) && (add < 0.5) );
376 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
378 int i, j, at, mol, nmol, nmolbox, count;
380 real z, zmin, zmax, mem_area;
383 int type = 0, block = 0;
386 mem_a = &(mem_p->mem_at);
387 snew(mol_id, mem_a->nr);
388 zmin = pos_ins->xmax[ZZ];
389 zmax = pos_ins->xmin[ZZ];
390 for (i = 0; i < mem_a->nr; i++)
392 at = mem_a->index[i];
393 if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
394 (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
395 (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
397 mol = get_mol_id(at, mtop, &type, &block);
399 for (j = 0; j < nmol; j++)
401 if (mol == mol_id[j])
429 srenew(mol_id, nmol);
430 mem_p->mol_id = mol_id;
432 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
434 gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
435 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
436 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
437 "your water layer is not thick enough.\n", zmax, zmin);
441 mem_p->zmed = (zmax-zmin)/2+zmin;
443 /*number of membrane molecules in protein box*/
444 nmolbox = count/mtop->molblock[block].natoms_mol;
445 /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
446 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
447 /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
448 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
450 return mem_p->mem_at.nr;
453 static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
454 gmx_bool bALLOW_ASYMMETRY)
456 int i, j, at, c, outsidesum, gctr = 0;
460 for (i = 0; i < pos_ins->pieces; i++)
462 idxsum += pos_ins->nidx[i];
465 if (idxsum != ins_at->nr)
467 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
470 snew(pos_ins->geom_cent, pos_ins->pieces);
471 for (i = 0; i < pos_ins->pieces; i++)
475 for (j = 0; j < DIM; j++)
477 pos_ins->geom_cent[i][j] = 0;
480 for (j = 0; j < pos_ins->nidx[i]; j++)
482 at = pos_ins->subindex[i][j];
483 copy_rvec(r[at], r_ins[gctr]);
484 if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
486 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
498 svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
501 if (!bALLOW_ASYMMETRY)
503 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
506 fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
507 i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
509 fprintf(stderr, "\n");
512 /* resize performed in the md loop */
513 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, rvec fac)
515 int i, j, k, at, c = 0;
516 for (k = 0; k < pos_ins->pieces; k++)
518 for (i = 0; i < pos_ins->nidx[k]; i++)
520 at = pos_ins->subindex[k][i];
521 for (j = 0; j < DIM; j++)
523 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
530 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
531 * The molecule to be embedded is already reduced in size. */
532 static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
533 rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
534 int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
536 int i, j, k, l, at, at2, mol_id;
537 int type = 0, block = 0;
538 int nrm, nupper, nlower;
539 real r_min_rad, z_lip, min_norm;
545 r_min_rad = probe_rad*probe_rad;
546 snew(rm_p->mol, mtop->mols.nr);
547 snew(rm_p->block, mtop->mols.nr);
550 for (i = 0; i < ins_at->nr; i++)
552 at = ins_at->index[i];
553 for (j = 0; j < rest_at->nr; j++)
555 at2 = rest_at->index[j];
556 pbc_dx(pbc, r[at], r[at2], dr);
558 if (norm2(dr) < r_min_rad)
560 mol_id = get_mol_id(at2, mtop, &type, &block);
562 for (l = 0; l < nrm; l++)
564 if (rm_p->mol[l] == mol_id)
572 rm_p->mol[nrm] = mol_id;
573 rm_p->block[nrm] = block;
576 for (l = 0; l < mem_p->nmol; l++)
578 if (mol_id == mem_p->mol_id[l])
580 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
584 z_lip /= mtop->molblock[block].natoms_mol;
585 if (z_lip < mem_p->zmed)
600 /*make sure equal number of lipids from upper and lower layer are removed */
601 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
603 snew(dist, mem_p->nmol);
604 snew(order, mem_p->nmol);
605 for (i = 0; i < mem_p->nmol; i++)
607 at = mtop->mols.index[mem_p->mol_id[i]];
608 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
609 if (pos_ins->pieces > 1)
612 min_norm = norm2(dr);
613 for (k = 1; k < pos_ins->pieces; k++)
615 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
616 if (norm2(dr_tmp) < min_norm)
618 min_norm = norm2(dr_tmp);
619 copy_rvec(dr_tmp, dr);
623 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
625 while (j >= 0 && dist[i] < dist[order[j]])
627 order[j+1] = order[j];
634 while (nupper != nlower)
636 mol_id = mem_p->mol_id[order[i]];
637 block = get_molblock(mol_id, mtop->nmolblock, mtop->molblock);
639 for (l = 0; l < nrm; l++)
641 if (rm_p->mol[l] == mol_id)
650 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
654 z_lip /= mtop->molblock[block].natoms_mol;
655 if (nupper > nlower && z_lip < mem_p->zmed)
657 rm_p->mol[nrm] = mol_id;
658 rm_p->block[nrm] = block;
662 else if (nupper < nlower && z_lip > mem_p->zmed)
664 rm_p->mol[nrm] = mol_id;
665 rm_p->block[nrm] = block;
674 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
682 srenew(rm_p->mol, nrm);
683 srenew(rm_p->block, nrm);
685 return nupper+nlower;
688 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
689 static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
690 t_block *ins_at, pos_ins_t *pos_ins)
692 int i, j, k, n, rm, mol_id, at, block;
694 atom_id *list, *new_mols;
695 unsigned char *new_egrp[egcNR];
699 snew(list, state->natoms);
701 for (i = 0; i < rm_p->nr; i++)
703 mol_id = rm_p->mol[i];
704 at = mtop->mols.index[mol_id];
705 block = rm_p->block[i];
706 mtop->molblock[block].nmol--;
707 for (j = 0; j < mtop->molblock[block].natoms_mol; j++)
712 mtop->mols.index[mol_id] = -1;
715 mtop->mols.nr -= rm_p->nr;
716 mtop->mols.nalloc_index -= rm_p->nr;
717 snew(new_mols, mtop->mols.nr);
718 for (i = 0; i < mtop->mols.nr+rm_p->nr; i++)
721 if (mtop->mols.index[i] != -1)
723 new_mols[j] = mtop->mols.index[i];
727 sfree(mtop->mols.index);
728 mtop->mols.index = new_mols;
731 state->nalloc = state->natoms;
732 snew(x_tmp, state->nalloc);
733 snew(v_tmp, state->nalloc);
735 for (i = 0; i < egcNR; i++)
737 if (groups->grpnr[i] != NULL)
739 groups->ngrpnr[i] = state->natoms;
740 snew(new_egrp[i], state->natoms);
745 for (i = 0; i < state->natoms+n; i++)
748 for (j = 0; j < n; j++)
759 for (j = 0; j < egcNR; j++)
761 if (groups->grpnr[j] != NULL)
763 new_egrp[j][i-rm] = groups->grpnr[j][i];
766 copy_rvec(state->x[i], x_tmp[i-rm]);
767 copy_rvec(state->v[i], v_tmp[i-rm]);
768 for (j = 0; j < ins_at->nr; j++)
770 if (i == ins_at->index[j])
772 ins_at->index[j] = i-rm;
776 for (j = 0; j < pos_ins->pieces; j++)
778 for (k = 0; k < pos_ins->nidx[j]; k++)
780 if (i == pos_ins->subindex[j][k])
782 pos_ins->subindex[j][k] = i-rm;
793 for (i = 0; i < egcNR; i++)
795 if (groups->grpnr[i] != NULL)
797 sfree(groups->grpnr[i]);
798 groups->grpnr[i] = new_egrp[i];
802 /* remove empty molblocks */
804 for (i = 0; i < mtop->nmolblock; i++)
806 if (mtop->molblock[i].nmol == 0)
812 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
815 mtop->nmolblock -= RMmolblock;
818 /* remove al bonded interactions from mtop for the molecule to be embedded */
819 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
822 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
824 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
825 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
826 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
827 * ins_at index group. MGWolf 050710 */
830 snew(bRM, mtop->nmoltype);
831 for (i = 0; i < mtop->nmoltype; i++)
836 for (i = 0; i < mtop->nmolblock; i++)
838 /*loop over molecule blocks*/
839 type = mtop->molblock[i].type;
840 natom = mtop->molblock[i].natoms_mol;
841 nmol = mtop->molblock[i].nmol;
843 for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)
845 /*loop over atoms in the block*/
846 at = j+atom1; /*atom index = block index + offset*/
849 for (m = 0; (m < ins_at->nr) && (bINS == FALSE); m++)
851 /*loop over atoms in insertion index group to determine if we're inserting one*/
852 if (at == ins_at->index[m])
859 atom1 += natom*nmol; /*update offset*/
862 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
866 for (i = 0; i < mtop->nmoltype; i++)
870 for (j = 0; j < F_LJ; j++)
872 mtop->moltype[i].ilist[j].nr = 0;
875 for (j = F_POSRES; j <= F_VSITEN; j++)
877 mtop->moltype[i].ilist[j].nr = 0;
886 /* Write a topology where the number of molecules is correct for the system after embedding */
887 static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
891 char buf[STRLEN], buf2[STRLEN], *temp;
892 int i, *nmol_rm, nmol, line;
893 char temporary_filename[STRLEN];
895 fpin = gmx_ffopen(topfile, "r");
896 strncpy(temporary_filename, "temp.topXXXXXX", STRLEN);
897 gmx_tmpnam(temporary_filename);
898 fpout = gmx_ffopen(temporary_filename, "w");
900 snew(nmol_rm, mtop->nmoltype);
901 for (i = 0; i < rm_p->nr; i++)
903 nmol_rm[rm_p->block[i]]++;
907 while (fgets(buf, STRLEN, fpin))
913 if ((temp = strchr(buf2, '\n')) != NULL)
921 if ((temp = strchr(buf2, '\n')) != NULL)
926 if (buf2[strlen(buf2)-1] == ']')
928 buf2[strlen(buf2)-1] = '\0';
931 if (gmx_strcasecmp(buf2, "molecules") == 0)
936 fprintf(fpout, "%s", buf);
938 else if (bMolecules == 1)
940 for (i = 0; i < mtop->nmolblock; i++)
942 nmol = mtop->molblock[i].nmol;
943 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
944 fprintf(fpout, "%s", buf);
948 else if (bMolecules == 2)
954 fprintf(fpout, "%s", buf);
959 fprintf(fpout, "%s", buf);
964 /* use gmx_ffopen to generate backup of topinout */
965 fpout = gmx_ffopen(topfile, "w");
967 rename(temporary_filename, topfile);
970 void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x)
972 /* Set new positions for the group to embed */
973 if (step_rel <= membed->it_xy)
975 membed->fac[0] += membed->xy_step;
976 membed->fac[1] += membed->xy_step;
978 else if (step_rel <= (membed->it_xy+membed->it_z))
980 membed->fac[2] += membed->z_step;
982 resize(membed->r_ins, x, membed->pos_ins, membed->fac);
985 /* We would like gn to be const as well, but C doesn't allow this */
986 /* TODO this is utility functionality (search for the index of a
987 string in a collection), so should be refactored and located more
988 centrally. Also, it nearly duplicates the same string in readir.c */
989 static int search_string(const char *s, int ng, char *gn[])
993 for (i = 0; (i < ng); i++)
995 if (gmx_strcasecmp(s, gn[i]) == 0)
1002 "Group %s selected for embedding was not found in the index file.\n"
1003 "Group names must match either [moleculetype] names or custom index group\n"
1004 "names, in which case you must supply an index file to the '-n' option\n"
1011 gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
1012 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
1014 char *ins, **gnames;
1015 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
1016 int ng, j, max_lip_rm, ins_grp_id, ins_nat, mem_nat, ntype, lip_rm, tpr_version;
1019 t_block *ins_at, *rest_at;
1023 gmx_groups_t *groups;
1024 gmx_bool bExcl = FALSE;
1027 char **piecename = NULL;
1028 gmx_membed_t membed = NULL;
1030 /* input variables */
1031 const char *membed_input;
1038 real probe_rad = 0.22;
1042 gmx_bool bALLOW_ASYMMETRY = FALSE;
1044 /* sanity check constants */ /* Issue a warning when: */
1045 const int membed_version = 58; /* tpr version is smaller */
1046 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1047 * and rest smaller than this value is probably too small */
1048 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1049 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1050 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1051 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1052 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1060 /* get input data out membed file */
1061 membed_input = opt2fn("-membed", nfile, fnm);
1062 get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1063 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1065 tpr_version = get_tpr_version(ftp2fn(efTPR, nfile, fnm));
1066 if (tpr_version < membed_version)
1068 gmx_fatal(FARGS, "Version of *.tpr file to old (%d). "
1069 "Rerun grompp with GROMACS version 4.0.3 or newer.\n", tpr_version);
1072 if (!EI_DYNAMICS(inputrec->eI) )
1074 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1079 gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1084 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1087 groups = &(mtop->groups);
1088 snew(gnames, groups->ngrpname);
1089 for (i = 0; i < groups->ngrpname; i++)
1091 gnames[i] = *(groups->grpname[i]);
1094 atoms = gmx_mtop_global_atoms(mtop);
1096 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1097 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1098 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1099 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1100 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1102 pos_ins->pieces = pieces;
1103 snew(pos_ins->nidx, pieces);
1104 snew(pos_ins->subindex, pieces);
1105 snew(piecename, pieces);
1108 fprintf(stderr, "\nSelect pieces to embed:\n");
1109 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1113 /*use whole embedded group*/
1114 snew(pos_ins->nidx, 1);
1115 snew(pos_ins->subindex, 1);
1116 pos_ins->nidx[0] = ins_at->nr;
1117 pos_ins->subindex[0] = ins_at->index;
1120 if (probe_rad < min_probe_rad)
1123 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1124 "in overlap between waters and the group to embed, which will result "
1125 "in Lincs errors etc.\n\n", warn);
1128 if (xy_fac < min_xy_init)
1131 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too smal.\n\n", warn, ins);
1134 if (it_xy < min_it_xy)
1137 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1138 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1141 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1144 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1145 " is probably too small.\nIncrease -nz or maxwarn.\n\n", warn, ins, it_z);
1148 if (it_xy+it_z > inputrec->nsteps)
1151 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1152 "number of steps in the tpr.\n\n", warn);
1156 if (inputrec->opts.ngfrz == 1)
1158 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1161 for (i = 0; i < inputrec->opts.ngfrz; i++)
1163 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1164 if (ins_grp_id == tmp_id)
1173 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1176 for (i = 0; i < DIM; i++)
1178 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1180 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1184 ng = groups->grps[egcENER].nr;
1187 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1190 for (i = 0; i < ng; i++)
1192 for (j = 0; j < ng; j++)
1194 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1197 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1198 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1200 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1201 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1202 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1210 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1211 "the freeze group");
1214 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1216 ins_nat = init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1217 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1218 check_types(ins_at, rest_at, mtop);
1220 mem_nat = init_mem_at(mem_p, mtop, state->x, state->box, pos_ins);
1222 prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p);
1223 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) )
1226 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1227 "This might cause pressure problems during the growth phase. Just try with\n"
1228 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
1229 "compressibility in the mdp-file or use no pressure coupling at all.\n\n", warn);
1234 gmx_fatal(FARGS, "Too many warnings.\n");
1237 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1238 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1239 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1241 /* Maximum number of lipids to be removed*/
1242 max_lip_rm = (int)(2*prot_area/mem_p->lip_area);
1243 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1245 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1246 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1247 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1248 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1250 /* resize the protein by xy and by z if necessary*/
1251 snew(r_ins, ins_at->nr);
1252 init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY);
1253 membed->fac[0] = membed->fac[1] = xy_fac;
1254 membed->fac[2] = z_fac;
1256 membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
1257 membed->z_step = (z_max-z_fac)/(double)(it_z-1);
1259 resize(r_ins, state->x, pos_ins, membed->fac);
1261 /* remove overlapping lipids and water from the membrane box*/
1262 /*mark molecules to be removed*/
1264 set_pbc(pbc, inputrec->ePBC, state->box);
1267 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, mem_p, pos_ins,
1268 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1269 lip_rm -= low_up_rm;
1273 for (i = 0; i < rm_p->nr; i++)
1275 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1279 for (i = 0; i < mtop->nmolblock; i++)
1282 for (j = 0; j < rm_p->nr; j++)
1284 if (rm_p->block[j] == i)
1289 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1292 if (lip_rm > max_lip_rm)
1295 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1296 "protein area\nTry making the -xyinit resize factor smaller or increase "
1297 "maxwarn.\n\n", warn);
1300 /*remove all lipids and waters overlapping and update all important structures*/
1301 rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
1303 rm_bonded_at = rm_bonded(ins_at, mtop);
1304 if (rm_bonded_at != ins_at->nr)
1306 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1307 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1308 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1313 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless, "
1314 "you can increase -maxwarn");
1317 if (ftp2bSet(efTOP, nfile, fnm))
1319 top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
1329 membed->it_xy = it_xy;
1330 membed->it_z = it_z;
1331 membed->pos_ins = pos_ins;
1332 membed->r_ins = r_ins;