2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2010,2011,2012,2013,2014, 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.
39 #include "gromacs/legacyheaders/typedefs.h"
40 #include "gromacs/legacyheaders/types/commrec.h"
41 #include "gromacs/utility/smalloc.h"
42 #include "gromacs/math/vec.h"
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/utility/futil.h"
45 #include "gromacs/essentialdynamics/edsam.h"
46 #include "gromacs/topology/index.h"
47 #include "gromacs/legacyheaders/names.h"
48 #include "gromacs/topology/mtop_util.h"
49 #include "gromacs/fileio/tpxio.h"
50 #include "gromacs/utility/cstringutil.h"
52 #include "gromacs/pbcutil/pbc.h"
53 #include "gromacs/legacyheaders/readinp.h"
54 #include "gromacs/gmxpreprocess/readir.h"
56 /* information about scaling center */
58 rvec xmin; /* smallest coordinates of all embedded molecules */
59 rvec xmax; /* largest coordinates of all embedded molecules */
60 rvec *geom_cent; /* scaling center of each independent molecule to embed */
61 int pieces; /* number of molecules to embed independently */
62 int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
63 atom_id **subindex; /* atomids for independent molecule *
64 * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
67 /* variables needed in do_md */
69 int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
70 int it_z; /* same, but for z */
71 real xy_step; /* stepsize used during resize in xy-plane */
72 real z_step; /* same, but in z */
73 rvec fac; /* initial scaling of the molecule to grow into the membrane */
74 rvec *r_ins; /* final coordinates of the molecule to grow */
75 pos_ins_t *pos_ins; /* scaling center for each piece to embed */
78 /* membrane related variables */
80 char *name; /* name of index group to embed molecule into (usually membrane) */
81 t_block mem_at; /* list all atoms in membrane */
82 int nmol; /* number of membrane molecules overlapping with the molecule to embed */
83 int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
84 real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
85 real zmin; /* minimum z coordinate of membrane */
86 real zmax; /* maximum z coordinate of membrane */
87 real zmed; /* median z coordinate of membrane */
90 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
91 * These will then be removed from the system */
93 int nr; /* number of molecules to remove */
94 int *mol; /* list of molecule ids to remove */
95 int *block; /* id of the molblock that the molecule to remove is part of */
98 /* Get the global molecule id, and the corresponding molecule type and id of the *
99 * molblock from the global atom nr. */
100 static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
105 gmx_mtop_atomlookup_t alook;
107 alook = gmx_mtop_atomlookup_settle_init(mtop);
108 gmx_mtop_atomnr_to_molblock_ind(alook, at, block, &mol_id, &atnr_mol);
109 for (i = 0; i < *block; i++)
111 mol_id += mtop->molblock[i].nmol;
113 *type = mtop->molblock[*block].type;
115 gmx_mtop_atomlookup_destroy(alook);
120 /* Get the id of the molblock from a global molecule id */
121 static int get_molblock(int mol_id, int nmblock, gmx_molblock_t *mblock)
126 for (i = 0; i < nmblock; i++)
128 nmol += mblock[i].nmol;
135 gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
140 static int get_tpr_version(const char *infile)
143 int version, generation;
145 read_tpxheader(infile, &header, TRUE, &version, &generation);
150 /* Get a list of all the molecule types that are present in a group of atoms. *
151 * Because all interaction within the group to embed are removed on the topology *
152 * level, if the same molecule type is found in another part of the system, these *
153 * would also be affected. Therefore we have to check if the embedded and rest group *
154 * share common molecule types. If so, membed will stop with an error. */
155 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
157 int i, j, nr, mol_id;
158 int type = 0, block = 0;
162 snew(tlist->index, at->nr);
163 for (i = 0; i < at->nr; i++)
166 mol_id = get_mol_id(at->index[i], mtop, &type, &block);
167 for (j = 0; j < nr; j++)
169 if (tlist->index[j] == type)
177 tlist->index[nr] = type;
181 srenew(tlist->index, nr);
185 /* Do the actual check of the molecule types between embedded and rest group */
186 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
188 t_block *ins_mtype, *rest_mtype;
193 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
194 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
196 for (i = 0; i < ins_mtype->nr; i++)
198 for (j = 0; j < rest_mtype->nr; j++)
200 if (ins_mtype->index[i] == rest_mtype->index[j])
202 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
203 "1. Your *.ndx and *.top do not match\n"
204 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
205 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
206 "we need to exclude all interactions between the atoms in the group to\n"
207 "insert, the same moleculetype can not be used in both groups. Change the\n"
208 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
209 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
210 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
215 sfree(ins_mtype->index);
216 sfree(rest_mtype->index);
221 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
222 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
223 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
229 wi = init_warning(TRUE, 0);
231 inp = read_inpfile(membed_input, &ninp, wi);
232 ITYPE ("nxy", *it_xy, 1000);
233 ITYPE ("nz", *it_z, 0);
234 RTYPE ("xyinit", *xy_fac, 0.5);
235 RTYPE ("xyend", *xy_max, 1.0);
236 RTYPE ("zinit", *z_fac, 1.0);
237 RTYPE ("zend", *z_max, 1.0);
238 RTYPE ("rad", *probe_rad, 0.22);
239 ITYPE ("ndiff", *low_up_rm, 0);
240 ITYPE ("maxwarn", *maxwarn, 0);
241 ITYPE ("pieces", *pieces, 1);
242 EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
244 write_inpfile(membed_input, ninp, inp, FALSE, wi);
247 /* Obtain the maximum and minimum coordinates of the group to be embedded */
248 static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins,
249 gmx_groups_t *groups, int ins_grp_id, real xy_max)
252 real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax;
253 const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
254 * determine the overlap between molecule to embed and membrane */
255 const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
256 snew(rest_at->index, state->natoms);
258 xmin = xmax = state->x[ins_at->index[0]][XX];
259 ymin = ymax = state->x[ins_at->index[0]][YY];
260 zmin = zmax = state->x[ins_at->index[0]][ZZ];
262 for (i = 0; i < state->natoms; i++)
264 gid = groups->grpnr[egcFREEZE][i];
265 if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
297 rest_at->index[c] = i;
303 srenew(rest_at->index, c);
305 if (xy_max > fac_inp_size)
307 pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
308 pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
310 pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
311 pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
315 pos_ins->xmin[XX] = xmin;
316 pos_ins->xmin[YY] = ymin;
318 pos_ins->xmax[XX] = xmax;
319 pos_ins->xmax[YY] = ymax;
322 if ( (zmax-zmin) < min_memthick)
324 pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
325 pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
329 pos_ins->xmin[ZZ] = zmin;
330 pos_ins->xmax[ZZ] = zmax;
336 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
337 * xy-plane and counting the number of occupied grid points */
338 static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
340 real x, y, dx = 0.15, dy = 0.15, area = 0.0;
341 real add, memmin, memmax;
344 /* min and max membrane coordinate are altered to reduce the influence of the *
346 memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
347 memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
349 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
351 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
357 at = ins_at->index[c];
358 if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
359 (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
360 (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
366 while ( (c < ins_at->nr) && (add < 0.5) );
375 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
377 int i, j, at, mol, nmol, nmolbox, count;
379 real z, zmin, zmax, mem_area;
382 int type = 0, block = 0;
385 mem_a = &(mem_p->mem_at);
386 snew(mol_id, mem_a->nr);
387 zmin = pos_ins->xmax[ZZ];
388 zmax = pos_ins->xmin[ZZ];
389 for (i = 0; i < mem_a->nr; i++)
391 at = mem_a->index[i];
392 if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
393 (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
394 (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
396 mol = get_mol_id(at, mtop, &type, &block);
398 for (j = 0; j < nmol; j++)
400 if (mol == mol_id[j])
428 srenew(mol_id, nmol);
429 mem_p->mol_id = mol_id;
431 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
433 gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
434 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
435 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
436 "your water layer is not thick enough.\n", zmax, zmin);
440 mem_p->zmed = (zmax-zmin)/2+zmin;
442 /*number of membrane molecules in protein box*/
443 nmolbox = count/mtop->molblock[block].natoms_mol;
444 /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
445 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
446 /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
447 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
449 return mem_p->mem_at.nr;
452 static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
453 gmx_bool bALLOW_ASYMMETRY)
455 int i, j, at, c, outsidesum, gctr = 0;
459 for (i = 0; i < pos_ins->pieces; i++)
461 idxsum += pos_ins->nidx[i];
464 if (idxsum != ins_at->nr)
466 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
469 snew(pos_ins->geom_cent, pos_ins->pieces);
470 for (i = 0; i < pos_ins->pieces; i++)
474 for (j = 0; j < DIM; j++)
476 pos_ins->geom_cent[i][j] = 0;
479 for (j = 0; j < pos_ins->nidx[i]; j++)
481 at = pos_ins->subindex[i][j];
482 copy_rvec(r[at], r_ins[gctr]);
483 if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
485 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
497 svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
500 if (!bALLOW_ASYMMETRY)
502 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
505 fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
506 i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
508 fprintf(stderr, "\n");
511 /* resize performed in the md loop */
512 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, rvec fac)
514 int i, j, k, at, c = 0;
515 for (k = 0; k < pos_ins->pieces; k++)
517 for (i = 0; i < pos_ins->nidx[k]; i++)
519 at = pos_ins->subindex[k][i];
520 for (j = 0; j < DIM; j++)
522 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
529 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
530 * The molecule to be embedded is already reduced in size. */
531 static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
532 rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
533 int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
535 int i, j, k, l, at, at2, mol_id;
536 int type = 0, block = 0;
537 int nrm, nupper, nlower;
538 real r_min_rad, z_lip, min_norm;
544 r_min_rad = probe_rad*probe_rad;
545 snew(rm_p->mol, mtop->mols.nr);
546 snew(rm_p->block, mtop->mols.nr);
549 for (i = 0; i < ins_at->nr; i++)
551 at = ins_at->index[i];
552 for (j = 0; j < rest_at->nr; j++)
554 at2 = rest_at->index[j];
555 pbc_dx(pbc, r[at], r[at2], dr);
557 if (norm2(dr) < r_min_rad)
559 mol_id = get_mol_id(at2, mtop, &type, &block);
561 for (l = 0; l < nrm; l++)
563 if (rm_p->mol[l] == mol_id)
571 rm_p->mol[nrm] = mol_id;
572 rm_p->block[nrm] = block;
575 for (l = 0; l < mem_p->nmol; l++)
577 if (mol_id == mem_p->mol_id[l])
579 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
583 z_lip /= mtop->molblock[block].natoms_mol;
584 if (z_lip < mem_p->zmed)
599 /*make sure equal number of lipids from upper and lower layer are removed */
600 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
602 snew(dist, mem_p->nmol);
603 snew(order, mem_p->nmol);
604 for (i = 0; i < mem_p->nmol; i++)
606 at = mtop->mols.index[mem_p->mol_id[i]];
607 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
608 if (pos_ins->pieces > 1)
611 min_norm = norm2(dr);
612 for (k = 1; k < pos_ins->pieces; k++)
614 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
615 if (norm2(dr_tmp) < min_norm)
617 min_norm = norm2(dr_tmp);
618 copy_rvec(dr_tmp, dr);
622 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
624 while (j >= 0 && dist[i] < dist[order[j]])
626 order[j+1] = order[j];
633 while (nupper != nlower)
635 mol_id = mem_p->mol_id[order[i]];
636 block = get_molblock(mol_id, mtop->nmolblock, mtop->molblock);
638 for (l = 0; l < nrm; l++)
640 if (rm_p->mol[l] == mol_id)
649 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
653 z_lip /= mtop->molblock[block].natoms_mol;
654 if (nupper > nlower && z_lip < mem_p->zmed)
656 rm_p->mol[nrm] = mol_id;
657 rm_p->block[nrm] = block;
661 else if (nupper < nlower && z_lip > mem_p->zmed)
663 rm_p->mol[nrm] = mol_id;
664 rm_p->block[nrm] = block;
673 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
681 srenew(rm_p->mol, nrm);
682 srenew(rm_p->block, nrm);
684 return nupper+nlower;
687 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
688 static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
689 t_block *ins_at, pos_ins_t *pos_ins)
691 int i, j, k, n, rm, mol_id, at, block;
693 atom_id *list, *new_mols;
694 unsigned char *new_egrp[egcNR];
698 snew(list, state->natoms);
700 for (i = 0; i < rm_p->nr; i++)
702 mol_id = rm_p->mol[i];
703 at = mtop->mols.index[mol_id];
704 block = rm_p->block[i];
705 mtop->molblock[block].nmol--;
706 for (j = 0; j < mtop->molblock[block].natoms_mol; j++)
711 mtop->mols.index[mol_id] = -1;
714 mtop->mols.nr -= rm_p->nr;
715 mtop->mols.nalloc_index -= rm_p->nr;
716 snew(new_mols, mtop->mols.nr);
717 for (i = 0; i < mtop->mols.nr+rm_p->nr; i++)
720 if (mtop->mols.index[i] != -1)
722 new_mols[j] = mtop->mols.index[i];
726 sfree(mtop->mols.index);
727 mtop->mols.index = new_mols;
730 state->nalloc = state->natoms;
731 snew(x_tmp, state->nalloc);
732 snew(v_tmp, state->nalloc);
734 for (i = 0; i < egcNR; i++)
736 if (groups->grpnr[i] != NULL)
738 groups->ngrpnr[i] = state->natoms;
739 snew(new_egrp[i], state->natoms);
744 for (i = 0; i < state->natoms+n; i++)
747 for (j = 0; j < n; j++)
758 for (j = 0; j < egcNR; j++)
760 if (groups->grpnr[j] != NULL)
762 new_egrp[j][i-rm] = groups->grpnr[j][i];
765 copy_rvec(state->x[i], x_tmp[i-rm]);
766 copy_rvec(state->v[i], v_tmp[i-rm]);
767 for (j = 0; j < ins_at->nr; j++)
769 if (i == ins_at->index[j])
771 ins_at->index[j] = i-rm;
775 for (j = 0; j < pos_ins->pieces; j++)
777 for (k = 0; k < pos_ins->nidx[j]; k++)
779 if (i == pos_ins->subindex[j][k])
781 pos_ins->subindex[j][k] = i-rm;
792 for (i = 0; i < egcNR; i++)
794 if (groups->grpnr[i] != NULL)
796 sfree(groups->grpnr[i]);
797 groups->grpnr[i] = new_egrp[i];
801 /* remove empty molblocks */
803 for (i = 0; i < mtop->nmolblock; i++)
805 if (mtop->molblock[i].nmol == 0)
811 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
814 mtop->nmolblock -= RMmolblock;
817 /* remove al bonded interactions from mtop for the molecule to be embedded */
818 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
821 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
823 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
824 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
825 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
826 * ins_at index group. MGWolf 050710 */
829 snew(bRM, mtop->nmoltype);
830 for (i = 0; i < mtop->nmoltype; i++)
835 for (i = 0; i < mtop->nmolblock; i++)
837 /*loop over molecule blocks*/
838 type = mtop->molblock[i].type;
839 natom = mtop->molblock[i].natoms_mol;
840 nmol = mtop->molblock[i].nmol;
842 for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)
844 /*loop over atoms in the block*/
845 at = j+atom1; /*atom index = block index + offset*/
848 for (m = 0; (m < ins_at->nr) && (bINS == FALSE); m++)
850 /*loop over atoms in insertion index group to determine if we're inserting one*/
851 if (at == ins_at->index[m])
858 atom1 += natom*nmol; /*update offset*/
861 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
865 for (i = 0; i < mtop->nmoltype; i++)
869 for (j = 0; j < F_LJ; j++)
871 mtop->moltype[i].ilist[j].nr = 0;
874 for (j = F_POSRES; j <= F_VSITEN; j++)
876 mtop->moltype[i].ilist[j].nr = 0;
885 /* Write a topology where the number of molecules is correct for the system after embedding */
886 static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
888 #define TEMP_FILENM "temp.top"
891 char buf[STRLEN], buf2[STRLEN], *temp;
892 int i, *nmol_rm, nmol, line;
894 fpin = gmx_ffopen(topfile, "r");
895 fpout = gmx_ffopen(TEMP_FILENM, "w");
897 snew(nmol_rm, mtop->nmoltype);
898 for (i = 0; i < rm_p->nr; i++)
900 nmol_rm[rm_p->block[i]]++;
904 while (fgets(buf, STRLEN, fpin))
910 if ((temp = strchr(buf2, '\n')) != NULL)
918 if ((temp = strchr(buf2, '\n')) != NULL)
923 if (buf2[strlen(buf2)-1] == ']')
925 buf2[strlen(buf2)-1] = '\0';
928 if (gmx_strcasecmp(buf2, "molecules") == 0)
933 fprintf(fpout, "%s", buf);
935 else if (bMolecules == 1)
937 for (i = 0; i < mtop->nmolblock; i++)
939 nmol = mtop->molblock[i].nmol;
940 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
941 fprintf(fpout, "%s", buf);
945 else if (bMolecules == 2)
951 fprintf(fpout, "%s", buf);
956 fprintf(fpout, "%s", buf);
961 /* use gmx_ffopen to generate backup of topinout */
962 fpout = gmx_ffopen(topfile, "w");
964 rename(TEMP_FILENM, topfile);
968 void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x)
970 /* Set new positions for the group to embed */
971 if (step_rel <= membed->it_xy)
973 membed->fac[0] += membed->xy_step;
974 membed->fac[1] += membed->xy_step;
976 else if (step_rel <= (membed->it_xy+membed->it_z))
978 membed->fac[2] += membed->z_step;
980 resize(membed->r_ins, x, membed->pos_ins, membed->fac);
983 gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
984 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
987 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
988 int ng, j, max_lip_rm, ins_grp_id, ins_nat, mem_nat, ntype, lip_rm, tpr_version;
991 t_block *ins_at, *rest_at;
995 gmx_groups_t *groups;
996 gmx_bool bExcl = FALSE;
999 char **piecename = NULL;
1000 gmx_membed_t membed = NULL;
1002 /* input variables */
1003 const char *membed_input;
1010 real probe_rad = 0.22;
1014 gmx_bool bALLOW_ASYMMETRY = FALSE;
1016 /* sanity check constants */ /* Issue a warning when: */
1017 const int membed_version = 58; /* tpr version is smaller */
1018 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1019 * and rest smaller than this value is probably too small */
1020 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1021 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1022 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1023 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1024 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1032 /* get input data out membed file */
1033 membed_input = opt2fn("-membed", nfile, fnm);
1034 get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1035 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1037 tpr_version = get_tpr_version(ftp2fn(efTPX, nfile, fnm));
1038 if (tpr_version < membed_version)
1040 gmx_fatal(FARGS, "Version of *.tpr file to old (%d). "
1041 "Rerun grompp with GROMACS version 4.0.3 or newer.\n", tpr_version);
1044 if (!EI_DYNAMICS(inputrec->eI) )
1046 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1051 gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1056 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1059 groups = &(mtop->groups);
1060 snew(gnames, groups->ngrpname);
1061 for (i = 0; i < groups->ngrpname; i++)
1063 gnames[i] = *(groups->grpname[i]);
1066 atoms = gmx_mtop_global_atoms(mtop);
1068 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1069 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1070 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1071 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1072 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1074 pos_ins->pieces = pieces;
1075 snew(pos_ins->nidx, pieces);
1076 snew(pos_ins->subindex, pieces);
1077 snew(piecename, pieces);
1080 fprintf(stderr, "\nSelect pieces to embed:\n");
1081 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1085 /*use whole embedded group*/
1086 snew(pos_ins->nidx, 1);
1087 snew(pos_ins->subindex, 1);
1088 pos_ins->nidx[0] = ins_at->nr;
1089 pos_ins->subindex[0] = ins_at->index;
1092 if (probe_rad < min_probe_rad)
1095 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1096 "in overlap between waters and the group to embed, which will result "
1097 "in Lincs errors etc.\n\n", warn);
1100 if (xy_fac < min_xy_init)
1103 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too smal.\n\n", warn, ins);
1106 if (it_xy < min_it_xy)
1109 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1110 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1113 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1116 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1117 " is probably too small.\nIncrease -nz or maxwarn.\n\n", warn, ins, it_z);
1120 if (it_xy+it_z > inputrec->nsteps)
1123 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1124 "number of steps in the tpr.\n\n", warn);
1128 if (inputrec->opts.ngfrz == 1)
1130 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1133 for (i = 0; i < inputrec->opts.ngfrz; i++)
1135 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1136 if (ins_grp_id == tmp_id)
1145 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1148 for (i = 0; i < DIM; i++)
1150 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1152 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1156 ng = groups->grps[egcENER].nr;
1159 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1162 for (i = 0; i < ng; i++)
1164 for (j = 0; j < ng; j++)
1166 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1169 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1170 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1172 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1173 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1174 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1182 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1183 "the freeze group");
1186 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1188 ins_nat = init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1189 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1190 check_types(ins_at, rest_at, mtop);
1192 mem_nat = init_mem_at(mem_p, mtop, state->x, state->box, pos_ins);
1194 prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p);
1195 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) )
1198 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1199 "This might cause pressure problems during the growth phase. Just try with\n"
1200 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
1201 "compressibility in the mdp-file or use no pressure coupling at all.\n\n", warn);
1206 gmx_fatal(FARGS, "Too many warnings.\n");
1209 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1210 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1211 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1213 /* Maximum number of lipids to be removed*/
1214 max_lip_rm = (int)(2*prot_area/mem_p->lip_area);
1215 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1217 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1218 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1219 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1220 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1222 /* resize the protein by xy and by z if necessary*/
1223 snew(r_ins, ins_at->nr);
1224 init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY);
1225 membed->fac[0] = membed->fac[1] = xy_fac;
1226 membed->fac[2] = z_fac;
1228 membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
1229 membed->z_step = (z_max-z_fac)/(double)(it_z-1);
1231 resize(r_ins, state->x, pos_ins, membed->fac);
1233 /* remove overlapping lipids and water from the membrane box*/
1234 /*mark molecules to be removed*/
1236 set_pbc(pbc, inputrec->ePBC, state->box);
1239 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, mem_p, pos_ins,
1240 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1241 lip_rm -= low_up_rm;
1245 for (i = 0; i < rm_p->nr; i++)
1247 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1251 for (i = 0; i < mtop->nmolblock; i++)
1254 for (j = 0; j < rm_p->nr; j++)
1256 if (rm_p->block[j] == i)
1261 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1264 if (lip_rm > max_lip_rm)
1267 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1268 "protein area\nTry making the -xyinit resize factor smaller or increase "
1269 "maxwarn.\n\n", warn);
1272 /*remove all lipids and waters overlapping and update all important structures*/
1273 rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
1275 rm_bonded_at = rm_bonded(ins_at, mtop);
1276 if (rm_bonded_at != ins_at->nr)
1278 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1279 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1280 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1285 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless, "
1286 "you can increase -maxwarn");
1289 if (ftp2bSet(efTOP, nfile, fnm))
1291 top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
1301 membed->it_xy = it_xy;
1302 membed->it_z = it_z;
1303 membed->pos_ins = pos_ins;
1304 membed->r_ins = r_ins;