2 * $Id: mdrun.c,v 1.139.2.9 2009/05/04 16:13:29 hess Exp $
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2012, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
54 #include "mtop_util.h"
62 /* information about scaling center */
64 rvec xmin; /* smallest coordinates of all embedded molecules */
65 rvec xmax; /* largest coordinates of all embedded molecules */
66 rvec *geom_cent; /* scaling center of each independent molecule to embed */
67 int pieces; /* number of molecules to embed independently */
68 int *nidx; /* n atoms for every independent embedded molecule (index in subindex) */
69 atom_id **subindex; /* atomids for independent molecule *
70 * atoms of piece i run from subindex[i][0] to subindex[i][nidx[i]] */
73 /* variables needed in do_md */
75 int it_xy; /* number of iterations (steps) used to grow something in the xy-plane */
76 int it_z; /* same, but for z */
77 real xy_step; /* stepsize used during resize in xy-plane */
78 real z_step; /* same, but in z */
79 rvec fac; /* initial scaling of the molecule to grow into the membrane */
80 rvec *r_ins; /* final coordinates of the molecule to grow */
81 pos_ins_t *pos_ins; /* scaling center for each piece to embed */
84 /* membrane related variables */
86 char *name; /* name of index group to embed molecule into (usually membrane) */
87 t_block mem_at; /* list all atoms in membrane */
88 int nmol; /* number of membrane molecules overlapping with the molecule to embed */
89 int *mol_id; /* list of molecules in membrane that overlap with the molecule to embed */
90 real lip_area; /* average area per lipid in membrane (only correct for homogeneous bilayers)*/
91 real zmin; /* minimum z coordinate of membrane */
92 real zmax; /* maximum z coordinate of membrane */
93 real zmed; /* median z coordinate of membrane */
96 /* Lists all molecules in the membrane that overlap with the molecule to be embedded. *
97 * These will then be removed from the system */
99 int nr; /* number of molecules to remove */
100 int *mol; /* list of molecule ids to remove */
101 int *block; /* id of the molblock that the molecule to remove is part of */
104 /* Get the global molecule id, and the corresponding molecule type and id of the *
105 * molblock from the global atom nr. */
106 static int get_mol_id(int at, gmx_mtop_t *mtop, int *type, int *block)
111 gmx_mtop_atomlookup_t alook;
113 alook = gmx_mtop_atomlookup_settle_init(mtop);
114 gmx_mtop_atomnr_to_molblock_ind(alook, at, block, &mol_id, &atnr_mol);
115 for (i = 0; i < *block; i++)
117 mol_id += mtop->molblock[i].nmol;
119 *type = mtop->molblock[*block].type;
121 gmx_mtop_atomlookup_destroy(alook);
126 /* Get the id of the molblock from a global molecule id */
127 static int get_molblock(int mol_id, int nmblock, gmx_molblock_t *mblock)
132 for (i = 0; i < nmblock; i++)
134 nmol += mblock[i].nmol;
141 gmx_fatal(FARGS, "mol_id %d larger than total number of molecules %d.\n", mol_id, nmol);
146 static int get_tpr_version(const char *infile)
149 int version, generation;
151 read_tpxheader(infile, &header, TRUE, &version, &generation);
156 /* Get a list of all the molecule types that are present in a group of atoms. *
157 * Because all interaction within the group to embed are removed on the topology *
158 * level, if the same molecule type is found in another part of the system, these *
159 * would also be affected. Therefore we have to check if the embedded and rest group *
160 * share common molecule types. If so, membed will stop with an error. */
161 static int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
163 int i, j, nr, mol_id;
164 int type = 0, block = 0;
168 snew(tlist->index, at->nr);
169 for (i = 0; i < at->nr; i++)
172 mol_id = get_mol_id(at->index[i], mtop, &type, &block);
173 for (j = 0; j < nr; j++)
175 if (tlist->index[j] == type)
183 tlist->index[nr] = type;
187 srenew(tlist->index, nr);
191 /* Do the actual check of the molecule types between embedded and rest group */
192 static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop)
194 t_block *ins_mtype, *rest_mtype;
199 ins_mtype->nr = get_mtype_list(ins_at, mtop, ins_mtype );
200 rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
202 for (i = 0; i < ins_mtype->nr; i++)
204 for (j = 0; j < rest_mtype->nr; j++)
206 if (ins_mtype->index[i] == rest_mtype->index[j])
208 gmx_fatal(FARGS, "Moleculetype %s is found both in the group to insert and the rest of the system.\n"
209 "1. Your *.ndx and *.top do not match\n"
210 "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
211 "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
212 "we need to exclude all interactions between the atoms in the group to\n"
213 "insert, the same moleculetype can not be used in both groups. Change the\n"
214 "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
215 "an appropriate *.itp file", *(mtop->moltype[rest_mtype->index[j]].name),
216 *(mtop->moltype[rest_mtype->index[j]].name), *(mtop->moltype[rest_mtype->index[j]].name));
221 sfree(ins_mtype->index);
222 sfree(rest_mtype->index);
227 static void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
228 int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
229 int *pieces, gmx_bool *bALLOW_ASYMMETRY)
235 wi = init_warning(TRUE, 0);
237 inp = read_inpfile(membed_input, &ninp, NULL, wi);
238 ITYPE ("nxy", *it_xy, 1000);
239 ITYPE ("nz", *it_z, 0);
240 RTYPE ("xyinit", *xy_fac, 0.5);
241 RTYPE ("xyend", *xy_max, 1.0);
242 RTYPE ("zinit", *z_fac, 1.0);
243 RTYPE ("zend", *z_max, 1.0);
244 RTYPE ("rad", *probe_rad, 0.22);
245 ITYPE ("ndiff", *low_up_rm, 0);
246 ITYPE ("maxwarn", *maxwarn, 0);
247 ITYPE ("pieces", *pieces, 1);
248 EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
250 write_inpfile(membed_input, ninp, inp, FALSE, wi);
253 /* Obtain the maximum and minimum coordinates of the group to be embedded */
254 static int init_ins_at(t_block *ins_at, t_block *rest_at, t_state *state, pos_ins_t *pos_ins,
255 gmx_groups_t *groups, int ins_grp_id, real xy_max)
258 real x, xmin, xmax, y, ymin, ymax, z, zmin, zmax;
259 const real min_memthick = 6.0; /* minimum thickness of the bilayer that will be used to *
260 * determine the overlap between molecule to embed and membrane */
261 const real fac_inp_size = 1.000001; /* scaling factor to obtain input_size + 0.000001 (comparing reals) */
262 snew(rest_at->index, state->natoms);
264 xmin = xmax = state->x[ins_at->index[0]][XX];
265 ymin = ymax = state->x[ins_at->index[0]][YY];
266 zmin = zmax = state->x[ins_at->index[0]][ZZ];
268 for (i = 0; i < state->natoms; i++)
270 gid = groups->grpnr[egcFREEZE][i];
271 if (groups->grps[egcFREEZE].nm_ind[gid] == ins_grp_id)
303 rest_at->index[c] = i;
309 srenew(rest_at->index, c);
311 if (xy_max > fac_inp_size)
313 pos_ins->xmin[XX] = xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
314 pos_ins->xmin[YY] = ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
316 pos_ins->xmax[XX] = xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
317 pos_ins->xmax[YY] = ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
321 pos_ins->xmin[XX] = xmin;
322 pos_ins->xmin[YY] = ymin;
324 pos_ins->xmax[XX] = xmax;
325 pos_ins->xmax[YY] = ymax;
328 if ( (zmax-zmin) < min_memthick)
330 pos_ins->xmin[ZZ] = zmin+(zmax-zmin)/2.0-0.5*min_memthick;
331 pos_ins->xmax[ZZ] = zmin+(zmax-zmin)/2.0+0.5*min_memthick;
335 pos_ins->xmin[ZZ] = zmin;
336 pos_ins->xmax[ZZ] = zmax;
342 /* Estimate the area of the embedded molecule by projecting all coordinates on a grid in the *
343 * xy-plane and counting the number of occupied grid points */
344 static real est_prot_area(pos_ins_t *pos_ins, rvec *r, t_block *ins_at, mem_t *mem_p)
346 real x, y, dx = 0.15, dy = 0.15, area = 0.0;
347 real add, memmin, memmax;
350 /* min and max membrane coordinate are altered to reduce the influence of the *
352 memmin = mem_p->zmin+0.1*(mem_p->zmax-mem_p->zmin);
353 memmax = mem_p->zmax-0.1*(mem_p->zmax-mem_p->zmin);
355 for (x = pos_ins->xmin[XX]; x < pos_ins->xmax[XX]; x += dx)
357 for (y = pos_ins->xmin[YY]; y < pos_ins->xmax[YY]; y += dy)
363 at = ins_at->index[c];
364 if ( (r[at][XX] >= x) && (r[at][XX] < x+dx) &&
365 (r[at][YY] >= y) && (r[at][YY] < y+dy) &&
366 (r[at][ZZ] > memmin) && (r[at][ZZ] < memmax) )
372 while ( (c < ins_at->nr) && (add < 0.5) );
381 static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
383 int i, j, at, mol, nmol, nmolbox, count;
385 real z, zmin, zmax, mem_area;
388 int type = 0, block = 0;
391 mem_a = &(mem_p->mem_at);
392 snew(mol_id, mem_a->nr);
393 zmin = pos_ins->xmax[ZZ];
394 zmax = pos_ins->xmin[ZZ];
395 for (i = 0; i < mem_a->nr; i++)
397 at = mem_a->index[i];
398 if ( (r[at][XX] > pos_ins->xmin[XX]) && (r[at][XX] < pos_ins->xmax[XX]) &&
399 (r[at][YY] > pos_ins->xmin[YY]) && (r[at][YY] < pos_ins->xmax[YY]) &&
400 (r[at][ZZ] > pos_ins->xmin[ZZ]) && (r[at][ZZ] < pos_ins->xmax[ZZ]) )
402 mol = get_mol_id(at, mtop, &type, &block);
404 for (j = 0; j < nmol; j++)
406 if (mol == mol_id[j])
434 srenew(mol_id, nmol);
435 mem_p->mol_id = mol_id;
437 if ((zmax-zmin) > (box[ZZ][ZZ]-0.5))
439 gmx_fatal(FARGS, "Something is wrong with your membrane. Max and min z values are %f and %f.\n"
440 "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
441 "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
442 "your water layer is not thick enough.\n", zmax, zmin);
446 mem_p->zmed = (zmax-zmin)/2+zmin;
448 /*number of membrane molecules in protein box*/
449 nmolbox = count/mtop->molblock[block].natoms_mol;
450 /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
451 mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
452 /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
453 mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
455 return mem_p->mem_at.nr;
458 static void init_resize(t_block *ins_at, rvec *r_ins, pos_ins_t *pos_ins, mem_t *mem_p, rvec *r,
459 gmx_bool bALLOW_ASYMMETRY)
461 int i, j, at, c, outsidesum, gctr = 0;
465 for (i = 0; i < pos_ins->pieces; i++)
467 idxsum += pos_ins->nidx[i];
470 if (idxsum != ins_at->nr)
472 gmx_fatal(FARGS, "Piecewise sum of inserted atoms not same as size of group selected to insert.");
475 snew(pos_ins->geom_cent, pos_ins->pieces);
476 for (i = 0; i < pos_ins->pieces; i++)
480 for (j = 0; j < DIM; j++)
482 pos_ins->geom_cent[i][j] = 0;
485 for (j = 0; j < pos_ins->nidx[i]; j++)
487 at = pos_ins->subindex[i][j];
488 copy_rvec(r[at], r_ins[gctr]);
489 if ( (r_ins[gctr][ZZ] < mem_p->zmax) && (r_ins[gctr][ZZ] > mem_p->zmin) )
491 rvec_inc(pos_ins->geom_cent[i], r_ins[gctr]);
503 svmul(1/(double)c, pos_ins->geom_cent[i], pos_ins->geom_cent[i]);
506 if (!bALLOW_ASYMMETRY)
508 pos_ins->geom_cent[i][ZZ] = mem_p->zmed;
511 fprintf(stderr, "Embedding piece %d with center of geometry: %f %f %f\n",
512 i, pos_ins->geom_cent[i][XX], pos_ins->geom_cent[i][YY], pos_ins->geom_cent[i][ZZ]);
514 fprintf(stderr, "\n");
517 /* resize performed in the md loop */
518 static void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins, rvec fac)
520 int i, j, k, at, c = 0;
521 for (k = 0; k < pos_ins->pieces; k++)
523 for (i = 0; i < pos_ins->nidx[k]; i++)
525 at = pos_ins->subindex[k][i];
526 for (j = 0; j < DIM; j++)
528 r[at][j] = pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
535 /* generate the list of membrane molecules that overlap with the molecule to be embedded. *
536 * The molecule to be embedded is already reduced in size. */
537 static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc, gmx_mtop_t *mtop,
538 rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad,
539 int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
541 int i, j, k, l, at, at2, mol_id;
542 int type = 0, block = 0;
543 int nrm, nupper, nlower;
544 real r_min_rad, z_lip, min_norm;
550 r_min_rad = probe_rad*probe_rad;
551 snew(rm_p->mol, mtop->mols.nr);
552 snew(rm_p->block, mtop->mols.nr);
555 for (i = 0; i < ins_at->nr; i++)
557 at = ins_at->index[i];
558 for (j = 0; j < rest_at->nr; j++)
560 at2 = rest_at->index[j];
561 pbc_dx(pbc, r[at], r[at2], dr);
563 if (norm2(dr) < r_min_rad)
565 mol_id = get_mol_id(at2, mtop, &type, &block);
567 for (l = 0; l < nrm; l++)
569 if (rm_p->mol[l] == mol_id)
577 rm_p->mol[nrm] = mol_id;
578 rm_p->block[nrm] = block;
581 for (l = 0; l < mem_p->nmol; l++)
583 if (mol_id == mem_p->mol_id[l])
585 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
589 z_lip /= mtop->molblock[block].natoms_mol;
590 if (z_lip < mem_p->zmed)
605 /*make sure equal number of lipids from upper and lower layer are removed */
606 if ( (nupper != nlower) && (!bALLOW_ASYMMETRY) )
608 snew(dist, mem_p->nmol);
609 snew(order, mem_p->nmol);
610 for (i = 0; i < mem_p->nmol; i++)
612 at = mtop->mols.index[mem_p->mol_id[i]];
613 pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr);
614 if (pos_ins->pieces > 1)
617 min_norm = norm2(dr);
618 for (k = 1; k < pos_ins->pieces; k++)
620 pbc_dx(pbc, r[at], pos_ins->geom_cent[k], dr_tmp);
621 if (norm2(dr_tmp) < min_norm)
623 min_norm = norm2(dr_tmp);
624 copy_rvec(dr_tmp, dr);
628 dist[i] = dr[XX]*dr[XX]+dr[YY]*dr[YY];
630 while (j >= 0 && dist[i] < dist[order[j]])
632 order[j+1] = order[j];
639 while (nupper != nlower)
641 mol_id = mem_p->mol_id[order[i]];
642 block = get_molblock(mol_id, mtop->nmolblock, mtop->molblock);
644 for (l = 0; l < nrm; l++)
646 if (rm_p->mol[l] == mol_id)
655 for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++)
659 z_lip /= mtop->molblock[block].natoms_mol;
660 if (nupper > nlower && z_lip < mem_p->zmed)
662 rm_p->mol[nrm] = mol_id;
663 rm_p->block[nrm] = block;
667 else if (nupper < nlower && z_lip > mem_p->zmed)
669 rm_p->mol[nrm] = mol_id;
670 rm_p->block[nrm] = block;
679 gmx_fatal(FARGS, "Trying to remove more lipid molecules than there are in the membrane");
687 srenew(rm_p->mol, nrm);
688 srenew(rm_p->block, nrm);
690 return nupper+nlower;
693 /*remove all lipids and waters overlapping and update all important structures (e.g. state and mtop)*/
694 static void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state,
695 t_block *ins_at, pos_ins_t *pos_ins)
697 int i, j, k, n, rm, mol_id, at, block;
699 atom_id *list, *new_mols;
700 unsigned char *new_egrp[egcNR];
704 snew(list, state->natoms);
706 for (i = 0; i < rm_p->nr; i++)
708 mol_id = rm_p->mol[i];
709 at = mtop->mols.index[mol_id];
710 block = rm_p->block[i];
711 mtop->molblock[block].nmol--;
712 for (j = 0; j < mtop->molblock[block].natoms_mol; j++)
717 mtop->mols.index[mol_id] = -1;
720 mtop->mols.nr -= rm_p->nr;
721 mtop->mols.nalloc_index -= rm_p->nr;
722 snew(new_mols, mtop->mols.nr);
723 for (i = 0; i < mtop->mols.nr+rm_p->nr; i++)
726 if (mtop->mols.index[i] != -1)
728 new_mols[j] = mtop->mols.index[i];
732 sfree(mtop->mols.index);
733 mtop->mols.index = new_mols;
736 state->nalloc = state->natoms;
737 snew(x_tmp, state->nalloc);
738 snew(v_tmp, state->nalloc);
740 for (i = 0; i < egcNR; i++)
742 if (groups->grpnr[i] != NULL)
744 groups->ngrpnr[i] = state->natoms;
745 snew(new_egrp[i], state->natoms);
750 for (i = 0; i < state->natoms+n; i++)
753 for (j = 0; j < n; j++)
764 for (j = 0; j < egcNR; j++)
766 if (groups->grpnr[j] != NULL)
768 new_egrp[j][i-rm] = groups->grpnr[j][i];
771 copy_rvec(state->x[i], x_tmp[i-rm]);
772 copy_rvec(state->v[i], v_tmp[i-rm]);
773 for (j = 0; j < ins_at->nr; j++)
775 if (i == ins_at->index[j])
777 ins_at->index[j] = i-rm;
781 for (j = 0; j < pos_ins->pieces; j++)
783 for (k = 0; k < pos_ins->nidx[j]; k++)
785 if (i == pos_ins->subindex[j][k])
787 pos_ins->subindex[j][k] = i-rm;
798 for (i = 0; i < egcNR; i++)
800 if (groups->grpnr[i] != NULL)
802 sfree(groups->grpnr[i]);
803 groups->grpnr[i] = new_egrp[i];
807 /* remove empty molblocks */
809 for (i = 0; i < mtop->nmolblock; i++)
811 if (mtop->molblock[i].nmol == 0)
817 mtop->molblock[i-RMmolblock] = mtop->molblock[i];
820 mtop->nmolblock -= RMmolblock;
823 /* remove al bonded interactions from mtop for the molecule to be embedded */
824 int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
827 int type, natom, nmol, at, atom1 = 0, rm_at = 0;
829 /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
830 /*this routine does not live as dangerously as it seems. There is namely a check in init_membed to make *
831 * sure that g_membed exits with a warning when there are molecules of the same type not in the *
832 * ins_at index group. MGWolf 050710 */
835 snew(bRM, mtop->nmoltype);
836 for (i = 0; i < mtop->nmoltype; i++)
841 for (i = 0; i < mtop->nmolblock; i++)
843 /*loop over molecule blocks*/
844 type = mtop->molblock[i].type;
845 natom = mtop->molblock[i].natoms_mol;
846 nmol = mtop->molblock[i].nmol;
848 for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)
850 /*loop over atoms in the block*/
851 at = j+atom1; /*atom index = block index + offset*/
854 for (m = 0; (m < ins_at->nr) && (bINS == FALSE); m++)
856 /*loop over atoms in insertion index group to determine if we're inserting one*/
857 if (at == ins_at->index[m])
864 atom1 += natom*nmol; /*update offset*/
867 rm_at += natom*nmol; /*increment bonded removal counter by # atoms in block*/
871 for (i = 0; i < mtop->nmoltype; i++)
875 for (j = 0; j < F_LJ; j++)
877 mtop->moltype[i].ilist[j].nr = 0;
880 for (j = F_POSRES; j <= F_VSITEN; j++)
882 mtop->moltype[i].ilist[j].nr = 0;
891 /* Write a topology where the number of molecules is correct for the system after embedding */
892 static void top_update(const char *topfile, char *ins, rm_t *rm_p, gmx_mtop_t *mtop)
894 #define TEMP_FILENM "temp.top"
897 char buf[STRLEN], buf2[STRLEN], *temp;
898 int i, *nmol_rm, nmol, line;
900 fpin = ffopen(topfile, "r");
901 fpout = ffopen(TEMP_FILENM, "w");
903 snew(nmol_rm, mtop->nmoltype);
904 for (i = 0; i < rm_p->nr; i++)
906 nmol_rm[rm_p->block[i]]++;
910 while (fgets(buf, STRLEN, fpin))
916 if ((temp = strchr(buf2, '\n')) != NULL)
924 if ((temp = strchr(buf2, '\n')) != NULL)
929 if (buf2[strlen(buf2)-1] == ']')
931 buf2[strlen(buf2)-1] = '\0';
934 if (gmx_strcasecmp(buf2, "molecules") == 0)
939 fprintf(fpout, "%s", buf);
941 else if (bMolecules == 1)
943 for (i = 0; i < mtop->nmolblock; i++)
945 nmol = mtop->molblock[i].nmol;
946 sprintf(buf, "%-15s %5d\n", *(mtop->moltype[mtop->molblock[i].type].name), nmol);
947 fprintf(fpout, "%s", buf);
951 else if (bMolecules == 2)
957 fprintf(fpout, "%s", buf);
962 fprintf(fpout, "%s", buf);
967 /* use ffopen to generate backup of topinout */
968 fpout = ffopen(topfile, "w");
970 rename(TEMP_FILENM, topfile);
974 void rescale_membed(int step_rel, gmx_membed_t membed, rvec *x)
976 /* Set new positions for the group to embed */
977 if (step_rel <= membed->it_xy)
979 membed->fac[0] += membed->xy_step;
980 membed->fac[1] += membed->xy_step;
982 else if (step_rel <= (membed->it_xy+membed->it_z))
984 membed->fac[2] += membed->z_step;
986 resize(membed->r_ins, x, membed->pos_ins, membed->fac);
989 gmx_membed_t init_membed(FILE *fplog, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop,
990 t_inputrec *inputrec, t_state *state, t_commrec *cr, real *cpt)
993 int i, rm_bonded_at, fr_id, fr_i = 0, tmp_id, warn = 0;
994 int ng, j, max_lip_rm, ins_grp_id, ins_nat, mem_nat, ntype, lip_rm, tpr_version;
997 t_block *ins_at, *rest_at;
1001 gmx_groups_t *groups;
1002 gmx_bool bExcl = FALSE;
1005 char **piecename = NULL;
1006 gmx_membed_t membed = NULL;
1008 /* input variables */
1009 const char *membed_input;
1016 real probe_rad = 0.22;
1020 gmx_bool bALLOW_ASYMMETRY = FALSE;
1022 /* sanity check constants */ /* Issue a warning when: */
1023 const int membed_version = 58; /* tpr version is smaller */
1024 const real min_probe_rad = 0.2199999; /* A probe radius for overlap between embedded molecule *
1025 * and rest smaller than this value is probably too small */
1026 const real min_xy_init = 0.0999999; /* the initial shrinking of the molecule to embed is smaller */
1027 const int min_it_xy = 1000; /* the number of steps to embed in xy-plane is smaller */
1028 const int min_it_z = 100; /* the number of steps to embed in z is smaller */
1029 const real prot_vs_box = 7.5; /* molecule to embed is large (more then prot_vs_box) with respect */
1030 const real box_vs_prot = 50; /* to the box size (less than box_vs_prot) */
1038 /* get input data out membed file */
1039 membed_input = opt2fn("-membed", nfile, fnm);
1040 get_input(membed_input, &xy_fac, &xy_max, &z_fac, &z_max, &it_xy, &it_z, &probe_rad, &low_up_rm,
1041 &maxwarn, &pieces, &bALLOW_ASYMMETRY);
1043 tpr_version = get_tpr_version(ftp2fn(efTPX, nfile, fnm));
1044 if (tpr_version < membed_version)
1046 gmx_fatal(FARGS, "Version of *.tpr file to old (%d). "
1047 "Rerun grompp with GROMACS version 4.0.3 or newer.\n", tpr_version);
1050 if (!EI_DYNAMICS(inputrec->eI) )
1052 gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
1057 gmx_input("Sorry, parallel g_membed is not yet fully functional.");
1062 fprintf(stderr, "\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
1065 groups = &(mtop->groups);
1066 snew(gnames, groups->ngrpname);
1067 for (i = 0; i < groups->ngrpname; i++)
1069 gnames[i] = *(groups->grpname[i]);
1072 atoms = gmx_mtop_global_atoms(mtop);
1074 fprintf(stderr, "\nSelect a group to embed in the membrane:\n");
1075 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(ins_at->nr), &(ins_at->index), &ins);
1076 ins_grp_id = search_string(ins, groups->ngrpname, gnames);
1077 fprintf(stderr, "\nSelect a group to embed %s into (e.g. the membrane):\n", ins);
1078 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), 1, &(mem_p->mem_at.nr), &(mem_p->mem_at.index), &(mem_p->name));
1080 pos_ins->pieces = pieces;
1081 snew(pos_ins->nidx, pieces);
1082 snew(pos_ins->subindex, pieces);
1083 snew(piecename, pieces);
1086 fprintf(stderr, "\nSelect pieces to embed:\n");
1087 get_index(&atoms, opt2fn_null("-mn", nfile, fnm), pieces, pos_ins->nidx, pos_ins->subindex, piecename);
1091 /*use whole embedded group*/
1092 snew(pos_ins->nidx, 1);
1093 snew(pos_ins->subindex, 1);
1094 pos_ins->nidx[0] = ins_at->nr;
1095 pos_ins->subindex[0] = ins_at->index;
1098 if (probe_rad < min_probe_rad)
1101 fprintf(stderr, "\nWarning %d:\nA probe radius (-rad) smaller than 0.2 nm can result "
1102 "in overlap between waters and the group to embed, which will result "
1103 "in Lincs errors etc.\n\n", warn);
1106 if (xy_fac < min_xy_init)
1109 fprintf(stderr, "\nWarning %d:\nThe initial size of %s is probably too smal.\n\n", warn, ins);
1112 if (it_xy < min_it_xy)
1115 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d)"
1116 " is probably too small.\nIncrease -nxy or.\n\n", warn, ins, it_xy);
1119 if ( (it_z < min_it_z) && ( z_fac < 0.99999999 || z_fac > 1.0000001) )
1122 fprintf(stderr, "\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d)"
1123 " is probably too small.\nIncrease -nz or maxwarn.\n\n", warn, ins, it_z);
1126 if (it_xy+it_z > inputrec->nsteps)
1129 fprintf(stderr, "\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the "
1130 "number of steps in the tpr.\n\n", warn);
1134 if (inputrec->opts.ngfrz == 1)
1136 gmx_fatal(FARGS, "You did not specify \"%s\" as a freezegroup.", ins);
1139 for (i = 0; i < inputrec->opts.ngfrz; i++)
1141 tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
1142 if (ins_grp_id == tmp_id)
1151 gmx_fatal(FARGS, "\"%s\" not as freezegroup defined in the mdp-file.", ins);
1154 for (i = 0; i < DIM; i++)
1156 if (inputrec->opts.nFreeze[fr_i][i] != 1)
1158 gmx_fatal(FARGS, "freeze dimensions for %s are not Y Y Y\n", ins);
1162 ng = groups->grps[egcENER].nr;
1165 gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
1168 for (i = 0; i < ng; i++)
1170 for (j = 0; j < ng; j++)
1172 if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
1175 if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) ||
1176 (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
1178 gmx_fatal(FARGS, "Energy exclusions \"%s\" and \"%s\" do not match the group "
1179 "to embed \"%s\"", *groups->grpname[groups->grps[egcENER].nm_ind[i]],
1180 *groups->grpname[groups->grps[egcENER].nm_ind[j]], ins);
1188 gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in "
1189 "the freeze group");
1192 /* Obtain the maximum and minimum coordinates of the group to be embedded */
1194 ins_nat = init_ins_at(ins_at, rest_at, state, pos_ins, groups, ins_grp_id, xy_max);
1195 /* Check that moleculetypes in insertion group are not part of the rest of the system */
1196 check_types(ins_at, rest_at, mtop);
1198 mem_nat = init_mem_at(mem_p, mtop, state->x, state->box, pos_ins);
1200 prot_area = est_prot_area(pos_ins, state->x, ins_at, mem_p);
1201 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) )
1204 fprintf(stderr, "\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
1205 "This might cause pressure problems during the growth phase. Just try with\n"
1206 "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
1207 "compressibility in the mdp-file or use no pressure coupling at all.\n\n", warn);
1212 gmx_fatal(FARGS, "Too many warnings.\n");
1215 printf("The estimated area of the protein in the membrane is %.3f nm^2\n", prot_area);
1216 printf("\nThere are %d lipids in the membrane part that overlaps the protein.\n"
1217 "The area per lipid is %.4f nm^2.\n", mem_p->nmol, mem_p->lip_area);
1219 /* Maximum number of lipids to be removed*/
1220 max_lip_rm = (int)(2*prot_area/mem_p->lip_area);
1221 printf("Maximum number of lipids that will be removed is %d.\n", max_lip_rm);
1223 printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
1224 "This resizing will be done with respect to the geometrical center of all protein atoms\n"
1225 "that span the membrane region, i.e. z between %.3f and %.3f\n\n",
1226 xy_fac, z_fac, mem_p->zmin, mem_p->zmax);
1228 /* resize the protein by xy and by z if necessary*/
1229 snew(r_ins, ins_at->nr);
1230 init_resize(ins_at, r_ins, pos_ins, mem_p, state->x, bALLOW_ASYMMETRY);
1231 membed->fac[0] = membed->fac[1] = xy_fac;
1232 membed->fac[2] = z_fac;
1234 membed->xy_step = (xy_max-xy_fac)/(double)(it_xy);
1235 membed->z_step = (z_max-z_fac)/(double)(it_z-1);
1237 resize(r_ins, state->x, pos_ins, membed->fac);
1239 /* remove overlapping lipids and water from the membrane box*/
1240 /*mark molecules to be removed*/
1242 set_pbc(pbc, inputrec->ePBC, state->box);
1245 lip_rm = gen_rm_list(rm_p, ins_at, rest_at, pbc, mtop, state->x, r_ins, mem_p, pos_ins,
1246 probe_rad, low_up_rm, bALLOW_ASYMMETRY);
1247 lip_rm -= low_up_rm;
1251 for (i = 0; i < rm_p->nr; i++)
1253 fprintf(fplog, "rm mol %d\n", rm_p->mol[i]);
1257 for (i = 0; i < mtop->nmolblock; i++)
1260 for (j = 0; j < rm_p->nr; j++)
1262 if (rm_p->block[j] == i)
1267 printf("Will remove %d %s molecules\n", ntype, *(mtop->moltype[mtop->molblock[i].type].name));
1270 if (lip_rm > max_lip_rm)
1273 fprintf(stderr, "\nWarning %d:\nTrying to remove a larger lipid area than the estimated "
1274 "protein area\nTry making the -xyinit resize factor smaller or increase "
1275 "maxwarn.\n\n", warn);
1278 /*remove all lipids and waters overlapping and update all important structures*/
1279 rm_group(inputrec, groups, mtop, rm_p, state, ins_at, pos_ins);
1281 rm_bonded_at = rm_bonded(ins_at, mtop);
1282 if (rm_bonded_at != ins_at->nr)
1284 fprintf(stderr, "Warning: The number of atoms for which the bonded interactions are removed is %d, "
1285 "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
1286 "molecule type as atoms that are not to be embedded.\n", rm_bonded_at, ins_at->nr);
1291 gmx_fatal(FARGS, "Too many warnings.\nIf you are sure these warnings are harmless, "
1292 "you can increase -maxwarn");
1295 if (ftp2bSet(efTOP, nfile, fnm))
1297 top_update(opt2fn("-mp", nfile, fnm), ins, rm_p, mtop);
1307 membed->it_xy = it_xy;
1308 membed->it_z = it_z;
1309 membed->pos_ins = pos_ins;
1310 membed->r_ins = r_ins;