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.
42 #include "types/commrec.h"
43 #include "gromacs/utility/smalloc.h"
44 #include "gromacs/math/vec.h"
46 #include "gromacs/utility/futil.h"
47 #include "gromacs/essentialdynamics/edsam.h"
48 #include "gromacs/topology/index.h"
50 #include "gromacs/topology/mtop_util.h"
51 #include "gromacs/fileio/tpxio.h"
52 #include "gromacs/utility/cstringutil.h"
54 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/gmxpreprocess/readir.h"
58 /* information about scaling center */
60 rvec xmin; /* smallest coordinates of all embedded molecules */
61 rvec xmax; /* largest coordinates of all embedded molecules */
62 rvec *geom_cent; /* scaling center of each independent molecule to embed */
63 int pieces; /* number of molecules to embed independently */
64 int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
65 atom_id **subindex; /* atomids for independent molecule *
66 * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
69 /* variables needed in do_md */
71 int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
72 int it_z; /* same, but for z */
73 real xy_step; /* stepsize used during resize in xy-plane */
74 real z_step; /* same, but in z */
75 rvec fac; /* initial scaling of the molecule to grow into the membrane */
76 rvec *r_ins; /* final coordinates of the molecule to grow */
77 pos_ins_t *pos_ins; /* scaling center for each piece to embed */
80 /* membrane related variables */
82 char *name; /* name of index group to embed molecule into (usually membrane) */
83 t_block mem_at; /* list all atoms in membrane */
84 int nmol; /* number of membrane molecules overlapping with the molecule to embed */
85 int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
86 real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
87 real zmin; /* minimum z coordinate of membrane */
88 real zmax; /* maximum z coordinate of membrane */
89 real zmed; /* median z coordinate of membrane */
92 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
93 * These will then be removed from the system */
95 int nr; /* number of molecules to remove */
96 int *mol; /* list of molecule ids to remove */
97 int *block; /* id of the molblock that the molecule to remove is part of */
100 /* Get the global molecule id, and the corresponding molecule type and id of the *
101 * molblock from the global atom nr. */
102 static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
107 gmx_mtop_atomlookup_t alook;
109 alook = gmx_mtop_atomlookup_settle_init(mtop);
110 gmx_mtop_atomnr_to_molblock_ind(alook, at, block, &mol_id, &atnr_mol);
111 for (i = 0; i < *block; i++)
113 mol_id += mtop->molblock[i].nmol;
115 *type = mtop->molblock[*block].type;
117 gmx_mtop_atomlookup_destroy(alook);
122 /* Get the id of the molblock from a global molecule id */
123 static int get_molblock(int mol_id, int nmblock, gmx_molblock_t *mblock)
128 for (i = 0; i < nmblock; i++)
130 nmol += mblock[i].nmol;
137 gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
142 static int get_tpr_version(const char *infile)
145 int version, generation;
147 read_tpxheader(infile, &header, TRUE, &version, &generation);
152 /* Get a list of all the molecule types that are present in a group of atoms. *
153 * Because all interaction within the group to embed are removed on the topology *
154 * level, if the same molecule type is found in another part of the system, these *
155 * would also be affected. Therefore we have to check if the embedded and rest group *
156 * share common molecule types. If so, membed will stop with an error. */
157 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
159 int i, j, nr, mol_id;
160 int type = 0, block = 0;
164 snew(tlist->index, at->nr);
165 for (i = 0; i < at->nr; i++)
168 mol_id = get_mol_id(at->index[i], mtop, &type, &block);
169 for (j = 0; j < nr; j++)
171 if (tlist->index[j] == type)
179 tlist->index[nr] = type;
183 srenew(tlist->index, nr);
187 /* Do the actual check of the molecule types between embedded and rest group */
188 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
190 t_block *ins_mtype, *rest_mtype;
195 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
196 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
198 for (i = 0; i < ins_mtype->nr; i++)
200 for (j = 0; j < rest_mtype->nr; j++)
202 if (ins_mtype->index[i] == rest_mtype->index[j])
204 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
205 "1. Your *.ndx and *.top do not match\n"
206 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
207 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
208 "we need to exclude all interactions between the atoms in the group to\n"
209 "insert, the same moleculetype can not be used in both groups. Change the\n"
210 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
211 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
212 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
217 sfree(ins_mtype->index);
218 sfree(rest_mtype->index);
223 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
224 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
225 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
231 wi = init_warning(TRUE, 0);
233 inp = read_inpfile(membed_input, &ninp, wi);
234 ITYPE ("nxy", *it_xy, 1000);
235 ITYPE ("nz", *it_z, 0);
236 RTYPE ("xyinit", *xy_fac, 0.5);
237 RTYPE ("xyend", *xy_max, 1.0);
238 RTYPE ("zinit", *z_fac, 1.0);
239 RTYPE ("zend", *z_max, 1.0);
240 RTYPE ("rad", *probe_rad, 0.22);
241 ITYPE ("ndiff", *low_up_rm, 0);
242 ITYPE ("maxwarn", *maxwarn, 0);
243 ITYPE ("pieces", *pieces, 1);
244 EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
246 write_inpfile(membed_input, ninp, inp, FALSE, wi);
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 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
353 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
359 at = ins_at->index[c];
360 if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
361 (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
362 (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
368 while ( (c < ins_at->nr) && (add < 0.5) );
377 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
379 int i, j, at, mol, nmol, nmolbox, count;
381 real z, zmin, zmax, mem_area;
384 int type = 0, block = 0;
387 mem_a = &(mem_p->mem_at);
388 snew(mol_id, mem_a->nr);
389 zmin = pos_ins->xmax[ZZ];
390 zmax = pos_ins->xmin[ZZ];
391 for (i = 0; i < mem_a->nr; i++)
393 at = mem_a->index[i];
394 if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
395 (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
396 (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
398 mol = get_mol_id(at, mtop, &type, &block);
400 for (j = 0; j < nmol; j++)
402 if (mol == mol_id[j])
430 srenew(mol_id, nmol);
431 mem_p->mol_id = mol_id;
433 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
435 gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
436 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
437 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
438 "your water layer is not thick enough.\n", zmax, zmin);
442 mem_p->zmed = (zmax-zmin)/2+zmin;
444 /*number of membrane molecules in protein box*/
445 nmolbox = count/mtop->molblock[block].natoms_mol;
446 /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
447 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
448 /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
449 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
451 return mem_p->mem_at.nr;
454 static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
455 gmx_bool bALLOW_ASYMMETRY)
457 int i, j, at, c, outsidesum, gctr = 0;
461 for (i = 0; i < pos_ins->pieces; i++)
463 idxsum += pos_ins->nidx[i];
466 if (idxsum != ins_at->nr)
468 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
471 snew(pos_ins->geom_cent, pos_ins->pieces);
472 for (i = 0; i < pos_ins->pieces; i++)
476 for (j = 0; j < DIM; j++)
478 pos_ins->geom_cent[i][j] = 0;
481 for (j = 0; j < pos_ins->nidx[i]; j++)
483 at = pos_ins->subindex[i][j];
484 copy_rvec(r[at], r_ins[gctr]);
485 if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
487 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
499 svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
502 if (!bALLOW_ASYMMETRY)
504 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
507 fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
508 i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
510 fprintf(stderr, "\n");
513 /* resize performed in the md loop */
514 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, rvec fac)
516 int i, j, k, at, c = 0;
517 for (k = 0; k < pos_ins->pieces; k++)
519 for (i = 0; i < pos_ins->nidx[k]; i++)
521 at = pos_ins->subindex[k][i];
522 for (j = 0; j < DIM; j++)
524 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
531 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
532 * The molecule to be embedded is already reduced in size. */
533 static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
534 rvec *r, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
535 int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
537 int i, j, k, l, at, at2, mol_id;
538 int type = 0, block = 0;
539 int nrm, nupper, nlower;
540 real r_min_rad, z_lip, min_norm;
546 r_min_rad = probe_rad*probe_rad;
547 snew(rm_p->mol, mtop->mols.nr);
548 snew(rm_p->block, mtop->mols.nr);
551 for (i = 0; i < ins_at->nr; i++)
553 at = ins_at->index[i];
554 for (j = 0; j < rest_at->nr; j++)
556 at2 = rest_at->index[j];
557 pbc_dx(pbc, r[at], r[at2], dr);
559 if (norm2(dr) < r_min_rad)
561 mol_id = get_mol_id(at2, mtop, &type, &block);
563 for (l = 0; l < nrm; l++)
565 if (rm_p->mol[l] == mol_id)
573 rm_p->mol[nrm] = mol_id;
574 rm_p->block[nrm] = block;
577 for (l = 0; l < mem_p->nmol; l++)
579 if (mol_id == mem_p->mol_id[l])
581 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
585 z_lip /= mtop->molblock[block].natoms_mol;
586 if (z_lip < mem_p->zmed)
601 /*make sure equal number of lipids from upper and lower layer are removed */
602 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
604 snew(dist, mem_p->nmol);
605 snew(order, mem_p->nmol);
606 for (i = 0; i < mem_p->nmol; i++)
608 at = mtop->mols.index[mem_p->mol_id[i]];
609 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
610 if (pos_ins->pieces > 1)
613 min_norm = norm2(dr);
614 for (k = 1; k < pos_ins->pieces; k++)
616 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
617 if (norm2(dr_tmp) < min_norm)
619 min_norm = norm2(dr_tmp);
620 copy_rvec(dr_tmp, dr);
624 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
626 while (j >= 0 && dist[i] < dist[order[j]])
628 order[j+1] = order[j];
635 while (nupper != nlower)
637 mol_id = mem_p->mol_id[order[i]];
638 block = get_molblock(mol_id, mtop->nmolblock, mtop->molblock);
640 for (l = 0; l < nrm; l++)
642 if (rm_p->mol[l] == mol_id)
651 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
655 z_lip /= mtop->molblock[block].natoms_mol;
656 if (nupper > nlower && z_lip < mem_p->zmed)
658 rm_p->mol[nrm] = mol_id;
659 rm_p->block[nrm] = block;
663 else if (nupper < nlower && z_lip > mem_p->zmed)
665 rm_p->mol[nrm] = mol_id;
666 rm_p->block[nrm] = block;
675 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
683 srenew(rm_p->mol, nrm);
684 srenew(rm_p->block, nrm);
686 return nupper+nlower;
689 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
690 static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
691 t_block *ins_at, pos_ins_t *pos_ins)
693 int i, j, k, n, rm, mol_id, at, block;
695 atom_id *list, *new_mols;
696 unsigned char *new_egrp[egcNR];
700 snew(list, state->natoms);
702 for (i = 0; i < rm_p->nr; i++)
704 mol_id = rm_p->mol[i];
705 at = mtop->mols.index[mol_id];
706 block = rm_p->block[i];
707 mtop->molblock[block].nmol--;
708 for (j = 0; j < mtop->molblock[block].natoms_mol; j++)
713 mtop->mols.index[mol_id] = -1;
716 mtop->mols.nr -= rm_p->nr;
717 mtop->mols.nalloc_index -= rm_p->nr;
718 snew(new_mols, mtop->mols.nr);
719 for (i = 0; i < mtop->mols.nr+rm_p->nr; i++)
722 if (mtop->mols.index[i] != -1)
724 new_mols[j] = mtop->mols.index[i];
728 sfree(mtop->mols.index);
729 mtop->mols.index = new_mols;
732 state->nalloc = state->natoms;
733 snew(x_tmp, state->nalloc);
734 snew(v_tmp, state->nalloc);
736 for (i = 0; i < egcNR; i++)
738 if (groups->grpnr[i] != NULL)
740 groups->ngrpnr[i] = state->natoms;
741 snew(new_egrp[i], state->natoms);
746 for (i = 0; i < state->natoms+n; i++)
749 for (j = 0; j < n; j++)
760 for (j = 0; j < egcNR; j++)
762 if (groups->grpnr[j] != NULL)
764 new_egrp[j][i-rm] = groups->grpnr[j][i];
767 copy_rvec(state->x[i], x_tmp[i-rm]);
768 copy_rvec(state->v[i], v_tmp[i-rm]);
769 for (j = 0; j < ins_at->nr; j++)
771 if (i == ins_at->index[j])
773 ins_at->index[j] = i-rm;
777 for (j = 0; j < pos_ins->pieces; j++)
779 for (k = 0; k < pos_ins->nidx[j]; k++)
781 if (i == pos_ins->subindex[j][k])
783 pos_ins->subindex[j][k] = i-rm;
794 for (i = 0; i < egcNR; i++)
796 if (groups->grpnr[i] != NULL)
798 sfree(groups->grpnr[i]);
799 groups->grpnr[i] = new_egrp[i];
803 /* remove empty molblocks */
805 for (i = 0; i < mtop->nmolblock; i++)
807 if (mtop->molblock[i].nmol == 0)
813 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
816 mtop->nmolblock -= RMmolblock;
819 /* remove al bonded interactions from mtop for the molecule to be embedded */
820 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
823 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
825 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
826 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
827 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
828 * ins_at index group. MGWolf 050710 */
831 snew(bRM, mtop->nmoltype);
832 for (i = 0; i < mtop->nmoltype; i++)
837 for (i = 0; i < mtop->nmolblock; i++)
839 /*loop over molecule blocks*/
840 type = mtop->molblock[i].type;
841 natom = mtop->molblock[i].natoms_mol;
842 nmol = mtop->molblock[i].nmol;
844 for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)
846 /*loop over atoms in the block*/
847 at = j+atom1; /*atom index = block index + offset*/
850 for (m = 0; (m < ins_at->nr) && (bINS == FALSE); m++)
852 /*loop over atoms in insertion index group to determine if we're inserting one*/
853 if (at == ins_at->index[m])
860 atom1 += natom*nmol; /*update offset*/
863 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
867 for (i = 0; i < mtop->nmoltype; i++)
871 for (j = 0; j < F_LJ; j++)
873 mtop->moltype[i].ilist[j].nr = 0;
876 for (j = F_POSRES; j <= F_VSITEN; j++)
878 mtop->moltype[i].ilist[j].nr = 0;
887 /* Write a topology where the number of molecules is correct for the system after embedding */
888 static void top_update(const char *topfile, rm_t *rm_p, gmx_mtop_t *mtop)
890 #define TEMP_FILENM "temp.top"
893 char buf[STRLEN], buf2[STRLEN], *temp;
894 int i, *nmol_rm, nmol, line;
896 fpin = gmx_ffopen(topfile, "r");
897 fpout = gmx_ffopen(TEMP_FILENM, "w");
899 snew(nmol_rm, mtop->nmoltype);
900 for (i = 0; i < rm_p->nr; i++)
902 nmol_rm[rm_p->block[i]]++;
906 while (fgets(buf, STRLEN, fpin))
912 if ((temp = strchr(buf2, '\n')) != NULL)
920 if ((temp = strchr(buf2, '\n')) != NULL)
925 if (buf2[strlen(buf2)-1] == ']')
927 buf2[strlen(buf2)-1] = '\0';
930 if (gmx_strcasecmp(buf2, "molecules") == 0)
935 fprintf(fpout, "%s", buf);
937 else if (bMolecules == 1)
939 for (i = 0; i < mtop->nmolblock; i++)
941 nmol = mtop->molblock[i].nmol;
942 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
943 fprintf(fpout, "%s", buf);
947 else if (bMolecules == 2)
953 fprintf(fpout, "%s", buf);
958 fprintf(fpout, "%s", buf);
963 /* use gmx_ffopen to generate backup of topinout */
964 fpout = gmx_ffopen(topfile, "w");
966 rename(TEMP_FILENM, 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 gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
986 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
989 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
990 int ng, j, max_lip_rm, ins_grp_id, ins_nat, mem_nat, ntype, lip_rm, tpr_version;
993 t_block *ins_at, *rest_at;
997 gmx_groups_t *groups;
998 gmx_bool bExcl = FALSE;
1001 char **piecename = NULL;
1002 gmx_membed_t membed = NULL;
1004 /* input variables */
1005 const char *membed_input;
1012 real probe_rad = 0.22;
1016 gmx_bool bALLOW_ASYMMETRY = FALSE;
1018 /* sanity check constants */ /* Issue a warning when: */
1019 const int membed_version = 58; /* tpr version is smaller */
1020 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1021 * and rest smaller than this value is probably too small */
1022 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1023 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1024 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1025 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1026 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1034 /* get input data out membed file */
1035 membed_input = opt2fn("-membed", nfile, fnm);
1036 get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1037 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1039 tpr_version = get_tpr_version(ftp2fn(efTPX, nfile, fnm));
1040 if (tpr_version < membed_version)
1042 gmx_fatal(FARGS, "Version of *.tpr file to old (%d). "
1043 "Rerun grompp with GROMACS version 4.0.3 or newer.\n", tpr_version);
1046 if (!EI_DYNAMICS(inputrec->eI) )
1048 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1053 gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1058 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1061 groups = &(mtop->groups);
1062 snew(gnames, groups->ngrpname);
1063 for (i = 0; i < groups->ngrpname; i++)
1065 gnames[i] = *(groups->grpname[i]);
1068 atoms = gmx_mtop_global_atoms(mtop);
1070 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1071 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1072 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1073 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1074 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1076 pos_ins->pieces = pieces;
1077 snew(pos_ins->nidx, pieces);
1078 snew(pos_ins->subindex, pieces);
1079 snew(piecename, pieces);
1082 fprintf(stderr, "\nSelect pieces to embed:\n");
1083 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1087 /*use whole embedded group*/
1088 snew(pos_ins->nidx, 1);
1089 snew(pos_ins->subindex, 1);
1090 pos_ins->nidx[0] = ins_at->nr;
1091 pos_ins->subindex[0] = ins_at->index;
1094 if (probe_rad < min_probe_rad)
1097 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1098 "in overlap between waters and the group to embed, which will result "
1099 "in Lincs errors etc.\n\n", warn);
1102 if (xy_fac < min_xy_init)
1105 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too smal.\n\n", warn, ins);
1108 if (it_xy < min_it_xy)
1111 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1112 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1115 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1118 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1119 " is probably too small.\nIncrease -nz or maxwarn.\n\n", warn, ins, it_z);
1122 if (it_xy+it_z > inputrec->nsteps)
1125 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1126 "number of steps in the tpr.\n\n", warn);
1130 if (inputrec->opts.ngfrz == 1)
1132 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1135 for (i = 0; i < inputrec->opts.ngfrz; i++)
1137 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1138 if (ins_grp_id == tmp_id)
1147 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1150 for (i = 0; i < DIM; i++)
1152 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1154 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1158 ng = groups->grps[egcENER].nr;
1161 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1164 for (i = 0; i < ng; i++)
1166 for (j = 0; j < ng; j++)
1168 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1171 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1172 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1174 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1175 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1176 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1184 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1185 "the freeze group");
1188 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1190 ins_nat = init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1191 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1192 check_types(ins_at, rest_at, mtop);
1194 mem_nat = init_mem_at(mem_p, mtop, state->x, state->box, pos_ins);
1196 prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p);
1197 if ( (prot_area > prot_vs_box) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX]) < box_vs_prot) )
1200 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1201 "This might cause pressure problems during the growth phase. Just try with\n"
1202 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
1203 "compressibility in the mdp-file or use no pressure coupling at all.\n\n", warn);
1208 gmx_fatal(FARGS, "Too many warnings.\n");
1211 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1212 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1213 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1215 /* Maximum number of lipids to be removed*/
1216 max_lip_rm = (int)(2*prot_area/mem_p->lip_area);
1217 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1219 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1220 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1221 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1222 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1224 /* resize the protein by xy and by z if necessary*/
1225 snew(r_ins, ins_at->nr);
1226 init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY);
1227 membed->fac[0] = membed->fac[1] = xy_fac;
1228 membed->fac[2] = z_fac;
1230 membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
1231 membed->z_step = (z_max-z_fac)/(double)(it_z-1);
1233 resize(r_ins, state->x, pos_ins, membed->fac);
1235 /* remove overlapping lipids and water from the membrane box*/
1236 /*mark molecules to be removed*/
1238 set_pbc(pbc, inputrec->ePBC, state->box);
1241 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, mem_p, pos_ins,
1242 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1243 lip_rm -= low_up_rm;
1247 for (i = 0; i < rm_p->nr; i++)
1249 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1253 for (i = 0; i < mtop->nmolblock; i++)
1256 for (j = 0; j < rm_p->nr; j++)
1258 if (rm_p->block[j] == i)
1263 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1266 if (lip_rm > max_lip_rm)
1269 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1270 "protein area\nTry making the -xyinit resize factor smaller or increase "
1271 "maxwarn.\n\n", warn);
1274 /*remove all lipids and waters overlapping and update all important structures*/
1275 rm_group(groups, mtop, rm_p, state, ins_at, pos_ins);
1277 rm_bonded_at = rm_bonded(ins_at, mtop);
1278 if (rm_bonded_at != ins_at->nr)
1280 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1281 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1282 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1287 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless, "
1288 "you can increase -maxwarn");
1291 if (ftp2bSet(efTOP, nfile, fnm))
1293 top_update(opt2fn("-mp", nfile, fnm), rm_p, mtop);
1303 membed->it_xy = it_xy;
1304 membed->it_z = it_z;
1305 membed->pos_ins = pos_ins;
1306 membed->r_ins = r_ins;