+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile.in
-Makefile
-config.h
-stamp-h.in
--- /dev/null
+/*
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gromacs Runs On Most of All Computer Systems
+ */
+
+#ifndef _gmx_membed_h
+#define _gmx_membed_h
+
+#include "typedefs.h"
+
+/* information about scaling center */
+typedef struct {
+ rvec xmin;
+ rvec xmax;
+ rvec *geom_cent;
+ int pieces;
+ int *nidx;
+ atom_id **subindex;
+} pos_ins_t;
+
+/* variables needed in do_md */
+typedef struct {
+ int it_xy;
+ int it_z;
+ real xy_step;
+ real z_step;
+ rvec fac; /* scaling factor */
+ rvec *r_ins; /* target coordinates */
+ pos_ins_t *pos_ins;
+} gmx_membed_t;
+
+/* initialisation of membed code */
+extern void init_membed(FILE *fplog, gmx_membed_t *membed, int nfile, const t_filenm fnm[],
+ gmx_mtop_t *mtop, t_inputrec *inputrec, t_state *state, t_commrec *cr,
+ real *cpt);
+
+/* rescaling the coordinates voor de membed code */
+extern void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x);
+#endif
#include "vsite.h"
#include "pull.h"
#include "update.h"
+#include "gmx_membed.h"
#ifdef __cplusplus
extern "C" {
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
#define unit_density_SI unit_mass_SI "/" unit_length_SI "^3"
#define unit_invvisc_SI unit_length_SI " " unit_time_SI "/" unit_mass_SI
+ /* The routines below can be used for converting units from or to GROMACS
+ internal units. */
+ enum { eg2cAngstrom, eg2cNm, eg2cBohr, eg2cKcal_Mole,
+ eg2cHartree, eg2cHartree_e, eg2cAngstrom3, eg2cCoulomb,
+ eg2cDebye, eg2cElectron, eg2cBuckingham, eg2cNR };
+
+ /* Convert value x to GROMACS units. Energy -> Energy, Length -> Length etc.
+ The type of x is deduced from unit,
+ which should be taken from the enum above. */
+ extern double convert2gmx(double x,int unit);
+
+ /* Convert value x from GROMACS units to the desired one.
+ The type of return value is deduced from unit, see above */
+ extern double gmx2convert(double x,int unit);
+
+ /* Convert the string to one of the units supported. Returns -1 if not found. */
+ extern int string2unit(char *string);
+
+ /* Convert the unit to a strong. Return NULL when unit is out of range. */
+ extern const char *unit2string(int unit);
+
#ifdef __cplusplus
}
#endif
*/
-char **split(char sep,char *str);
+char **split(char sep,const char *str);
/* Implementation of the well-known Perl function split */
gmx_large_int_t str_to_large_int_t(const char *str, char **endptr);
+++ /dev/null
-Makefile.in
-Makefile
\ No newline at end of file
return (m[XX][XX]+m[YY][YY]+m[ZZ][ZZ]);
}
-static gmx_inline real _divide(real a,real b,const char *file,int line)
+static gmx_inline real _divide_err(real a,real b,const char *file,int line)
{
if (fabs(b) <= GMX_REAL_MIN)
gmx_fatal(FARGS,"Dividing by zero, file %s, line %d",file,line);
-box[ZZ][XX]*box[ZZ][XX]-box[ZZ][YY]*box[ZZ][YY]);
}
-#define divide(a,b) _divide((a),(b),__FILE__,__LINE__)
+#define divide_err(a,b) _divide_err((a),(b),__FILE__,__LINE__)
#define mod(a,b) _mod((a),(b),__FILE__,__LINE__)
#ifdef __cplusplus
+++ /dev/null
-Makefile.in
-Makefile
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
--- /dev/null
+GMXRC
+GMXRC.*
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
--- /dev/null
+template
+gromacs
+Makefile.*-*
+++ /dev/null
-Makefile
-Makefile.in
-363
+364
��ߦ��ߨ���߬�����������ߦ��߷���ߋ��߶�ߦ�����������������߶��
�߳���ߋ��߳���߫���ߨ���߫������������
�����߽����߲�߷������ײ��������
�����،ߑ�߈��߆��ߜ��ߍ���ߐ�ߞ�ߚ���������������߸��������
�ߓ���ߋ�߈����ߋ���߶ߙ���ߓ���߶ߛ�ߌ�������������������
���߶ߗ���ߚ���������ߓ�����ߋ���ߚ���������ߚ�������߯������
+�����ߐ�ߜ�����������ߞ��ߋ��ߌ������߲�ߌ������ߖ�߈�������ߑ�����ߜ��������ײ�������
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
+++ /dev/null
-Makefile
-Makefile.in
-stamp-h.in
-stamp-h
-config.h
-config.h.in
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
main.c maths.c matio.c mshift.c \
mtop_util.c mtxio.c mvdata.c names.c \
network.c nrama.c nrjac.c nrnb.c \
- pargs.c pbc.c pdbio.c princ.c \
+ pargs.c pbc.c pdbio.c physics.c princ.c \
rando.c random.c gmx_random.c rbin.c \
readinp.c replace.c rmpbc.c shift_util.c \
sortwater.c smalloc.c statutil.c sfactor.c \
FORCE:
endif
+LDADD = ../gmxlib/libgmx@LIBSUFFIX@.la
+
+EXTRA_PROGRAMS = physics_test
EXTRA_DIST = version.h libgmx.pc.cmakein
--- /dev/null
+/* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
+ * $Id: molprop_util.c,v 1.51 2009/06/01 06:13:18 spoel Exp $
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 4.0.99
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2008, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Groningen Machine for Chemical Simulation
+ */
+#include <stdio.h>
+#include "string2.h"
+#include "physics.h"
+
+double convert2gmx(double x,int unit)
+{
+ switch (unit)
+ {
+ case eg2cAngstrom:
+ return x*A2NM;
+ case eg2cNm:
+ return x;
+ case eg2cBohr:
+ return x*BOHR2NM;
+ case eg2cKcal_Mole:
+ return x/CAL2JOULE;
+ case eg2cHartree:
+ return x*ONE_4PI_EPS0/BOHR2NM;
+ case eg2cHartree_e:
+ return x*ONE_4PI_EPS0/BOHR2NM;
+ case eg2cAngstrom3:
+ return x*A2NM*A2NM*A2NM;
+ case eg2cCoulomb:
+ return x/E_CHARGE;
+ case eg2cDebye:
+ return x*DEBYE2ENM;
+ case eg2cElectron:
+ return x;
+ case eg2cBuckingham:
+ return x*A2NM*DEBYE2ENM;
+ default:
+ fprintf(stderr,"Unknown unit %d, not converting.\n",unit);
+ }
+ return x;
+}
+
+double gmx2convert(double x,int unit)
+{
+ switch (unit)
+ {
+ case eg2cAngstrom:
+ return x/A2NM;
+ case eg2cNm:
+ return x;
+ case eg2cBohr:
+ return x/BOHR2NM;
+ case eg2cKcal_Mole:
+ return x*CAL2JOULE;
+ case eg2cHartree:
+ return x/(ONE_4PI_EPS0/BOHR2NM);
+ case eg2cHartree_e:
+ return x/(ONE_4PI_EPS0/BOHR2NM);
+ case eg2cAngstrom3:
+ return x/(A2NM*A2NM*A2NM);
+ case eg2cCoulomb:
+ return x*E_CHARGE;
+ case eg2cDebye:
+ return x/DEBYE2ENM;
+ case eg2cElectron:
+ return x;
+ case eg2cBuckingham:
+ return x/(A2NM*DEBYE2ENM);
+ default:
+ fprintf(stderr,"Unknown unit %d, not converting.\n",unit);
+ }
+ return x;
+}
+
+/* This has to have the same order as the enums. */
+static const char *eg2c_names[eg2cNR] = {
+ "Angstrom", "Nm", "Bohr", "Kcal_Mole",
+ "Hartree", "Hartree_e", "Angstrom3", "Coulomb",
+ "Debye", "Electron", "Buckingham"
+};
+
+int string2unit(char *string)
+{
+ int i;
+
+ for(i=0; (i<eg2cNR); i++)
+ if (strcasecmp(string,eg2c_names[i]) == 0)
+ return i;
+ return -1;
+}
+
+const char *unit2string(int unit)
+{
+ if ((unit >= 0) && (unit < eg2cNR))
+ return eg2c_names[unit];
+
+ return NULL;
+}
--- /dev/null
+#include <stdio.h>
+#include "physics.h"
+
+int main(int argc,char *argv[])
+{
+ int i;
+ double x,y,z;
+
+ x = 3.25;
+ for(i=0; (i<eg2cNR); i++) {
+ y = gmx2convert(x,i);
+ z = convert2gmx(y,i);
+ printf("Converted %g [type %d] to %g and back to %g. Diff %g\n",
+ x,i,y,z,x-z);
+ }
+}
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
return b2;
}
-char **split(char sep,char *str)
+char **split(char sep,const char *str)
{
char **ptr = NULL;
int n,nn,nptr = 0;
+++ /dev/null
-Makefile
-Makefile.in
-.libs
-.deps
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
.deps
.libs
+g_luck
+g_protonate
+g_x2top
+gmxcheck
+gmxdump
+grompp
+mdrun
+pdb2gmx
+tpbconv
g_x2top_SOURCES = g_x2top.c nm2type.c g_x2top.h
mdrun_SOURCES = \
- gctio.c \
+ gctio.c gmx_membed.c gmx_membed.h \
ionize.c ionize.h xmdrun.h \
do_gct.c repl_ex.c repl_ex.h \
xutils.c runner.c md.c mdrun.c \
--- /dev/null
+/*
+ * $Id: mdrun.c,v 1.139.2.9 2009/05/04 16:13:29 hess Exp $
+ *
+ * This source code is part of
+ *
+ * G R O M A C S
+ *
+ * GROningen MAchine for Chemical Simulations
+ *
+ * VERSION 3.2.0
+ * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team,
+ * check out http://www.gromacs.org for more information.
+
+ * This program is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU General Public License
+ * as published by the Free Software Foundation; either version 2
+ * of the License, or (at your option) any later version.
+ *
+ * If you want to redistribute modifications, please consider that
+ * scientific software is very special. Version control is crucial -
+ * bugs must be traceable. We will be happy to consider code for
+ * inclusion in the official distribution, but derived work must not
+ * be called official GROMACS. Details are found in the README & COPYING
+ * files - if they are missing, get the official version at www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the papers on the package - you can find them in the top README file.
+ *
+ * For more info, check our website at http://www.gromacs.org
+ *
+ * And Hey:
+ * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
+ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <signal.h>
+#include <stdlib.h>
+#include "typedefs.h"
+#include "smalloc.h"
+#include "sysstuff.h"
+#include "vec.h"
+#include "statutil.h"
+#include "macros.h"
+#include "copyrite.h"
+#include "main.h"
+#include "futil.h"
+#include "edsam.h"
+#include "index.h"
+#include "physics.h"
+#include "names.h"
+#include "mtop_util.h"
+#include "tpxio.h"
+#include "string2.h"
+#include "gmx_membed.h"
+#include "pbc.h"
+#include "readinp.h"
+
+typedef struct {
+ int id;
+ char *name;
+ int nr;
+ int natoms; /*nr of atoms per lipid*/
+ int mol1; /*id of the first lipid molecule*/
+ real area;
+} lip_t;
+
+typedef struct {
+ char *name;
+ t_block mem_at;
+ int *mol_id;
+ int nmol;
+ real lip_area;
+ real zmin;
+ real zmax;
+ real zmed;
+} mem_t;
+
+typedef struct {
+ int *mol;
+ int *block;
+ int nr;
+} rm_t;
+
+int search_string(char *s,int ng,char ***gn)
+{
+ int i;
+
+ for(i=0; (i<ng); i++)
+ if (gmx_strcasecmp(s,*gn[i]) == 0)
+ return i;
+
+ gmx_fatal(FARGS,"Group %s not found in indexfile.\nMaybe you have non-default groups in your mdp file, while not using the '-n' option of grompp.\nIn that case use the '-n' option.\n",s);
+
+ return -1;
+}
+
+int get_mol_id(int at,int nmblock,gmx_molblock_t *mblock, int *type, int *block)
+{
+ int mol_id=0;
+ int i;
+
+ for(i=0;i<nmblock;i++)
+ {
+ if(at<(mblock[i].nmol*mblock[i].natoms_mol))
+ {
+ mol_id+=at/mblock[i].natoms_mol;
+ *type = mblock[i].type;
+ *block = i;
+ return mol_id;
+ } else {
+ at-= mblock[i].nmol*mblock[i].natoms_mol;
+ mol_id+=mblock[i].nmol;
+ }
+ }
+
+ gmx_fatal(FARGS,"Something is wrong in mol ids, at %d, mol_id %d",at,mol_id);
+
+ return -1;
+}
+
+int get_block(int mol_id,int nmblock,gmx_molblock_t *mblock)
+{
+ int i;
+ int nmol=0;
+
+ for(i=0;i<nmblock;i++)
+ {
+ nmol+=mblock[i].nmol;
+ if(mol_id<nmol)
+ return i;
+ }
+
+ gmx_fatal(FARGS,"mol_id %d larger than total number of molecules %d.\n",mol_id,nmol);
+
+ return -1;
+}
+
+int get_tpr_version(const char *infile)
+{
+ char buf[STRLEN];
+ gmx_bool bDouble;
+ int precision,fver;
+ t_fileio *fio;
+
+ fio = open_tpx(infile,"r");
+ gmx_fio_checktype(fio);
+
+ precision = sizeof(real);
+
+ gmx_fio_do_string(fio,buf);
+ if (strncmp(buf,"VERSION",7))
+ gmx_fatal(FARGS,"Can not read file %s,\n"
+ " this file is from a Gromacs version which is older than 2.0\n"
+ " Make a new one with grompp or use a gro or pdb file, if possible",
+ gmx_fio_getname(fio));
+ gmx_fio_do_int(fio,precision);
+ bDouble = (precision == sizeof(double));
+ if ((precision != sizeof(float)) && !bDouble)
+ gmx_fatal(FARGS,"Unknown precision in file %s: real is %d bytes "
+ "instead of %d or %d",
+ gmx_fio_getname(fio),precision,sizeof(float),sizeof(double));
+ gmx_fio_setprecision(fio,bDouble);
+ fprintf(stderr,"Reading file %s, %s (%s precision)\n",
+ gmx_fio_getname(fio),buf,bDouble ? "double" : "single");
+
+ gmx_fio_do_int(fio,fver);
+
+ close_tpx(fio);
+
+ return fver;
+}
+
+int get_mtype_list(t_block *at, gmx_mtop_t *mtop, t_block *tlist)
+{
+ int i,j,nr,mol_id;
+ int type=0,block=0;
+ gmx_bool bNEW;
+
+ nr=0;
+ snew(tlist->index,at->nr);
+ for (i=0;i<at->nr;i++)
+ {
+ bNEW=TRUE;
+ mol_id = get_mol_id(at->index[i],mtop->nmolblock,mtop->molblock,&type,&block);
+ for(j=0;j<nr;j++)
+ {
+ if(tlist->index[j]==type)
+ bNEW=FALSE;
+ }
+ if(bNEW==TRUE)
+ {
+ tlist->index[nr]=type;
+ nr++;
+ }
+ }
+
+ srenew(tlist->index,nr);
+ return nr;
+}
+
+void check_types(t_block *ins_at,t_block *rest_at,gmx_mtop_t *mtop)
+{
+ t_block *ins_mtype,*rest_mtype;
+ int i,j;
+
+ snew(ins_mtype,1);
+ snew(rest_mtype,1);
+ ins_mtype->nr = get_mtype_list(ins_at , mtop, ins_mtype );
+ rest_mtype->nr = get_mtype_list(rest_at, mtop, rest_mtype);
+
+ for(i=0;i<ins_mtype->nr;i++)
+ {
+ for(j=0;j<rest_mtype->nr;j++)
+ {
+ if(ins_mtype->index[i]==rest_mtype->index[j])
+ gmx_fatal(FARGS,"Moleculetype %s is found both in the group to insert and the rest of the system.\n"
+ "1. Your *.ndx and *.top do not match\n"
+ "2. You are inserting some molecules of type %s (for example xray-solvent), while\n"
+ "the same moleculetype is also used in the rest of the system (solvent box). Because\n"
+ "we need to exclude all interactions between the atoms in the group to\n"
+ "insert, the same moleculetype can not be used in both groups. Change the\n"
+ "moleculetype of the molecules %s in the inserted group. Do not forget to provide\n"
+ "an appropriate *.itp file",*(mtop->moltype[rest_mtype->index[j]].name),
+ *(mtop->moltype[rest_mtype->index[j]].name),*(mtop->moltype[rest_mtype->index[j]].name));
+ }
+ }
+
+ sfree(ins_mtype->index);
+ sfree(rest_mtype->index);
+ sfree(ins_mtype);
+ sfree(rest_mtype);
+}
+
+void get_input(const char *membed_input, real *xy_fac, real *xy_max, real *z_fac, real *z_max,
+ int *it_xy, int *it_z, real *probe_rad, int *low_up_rm, int *maxwarn,
+ int *pieces, gmx_bool *bALLOW_ASYMMETRY)
+{
+ warninp_t wi;
+ t_inpfile *inp;
+ int ninp;
+
+ wi = init_warning(TRUE,0);
+
+ inp = read_inpfile(membed_input, &ninp, NULL, wi);
+ ITYPE ("nxy", *it_xy, 1000);
+ ITYPE ("nz", *it_z, 0);
+ RTYPE ("xyinit", *xy_fac, 0.5);
+ RTYPE ("xyend", *xy_max, 1.0);
+ RTYPE ("zinit", *z_fac, 1.0);
+ RTYPE ("zend", *z_max, 1.0);
+ RTYPE ("rad", *probe_rad, 0.22);
+ ITYPE ("ndiff", *low_up_rm, 0);
+ ITYPE ("maxwarn", *maxwarn, 0);
+ ITYPE ("pieces", *pieces, 1);
+ EETYPE("asymmetry", *bALLOW_ASYMMETRY, yesno_names);
+
+ write_inpfile(membed_input,ninp,inp,FALSE,wi);
+}
+
+int init_ins_at(t_block *ins_at,t_block *rest_at,t_state *state, pos_ins_t *pos_ins,gmx_groups_t *groups,int ins_grp_id, real xy_max)
+{
+ int i,gid,c=0;
+ real x,xmin,xmax,y,ymin,ymax,z,zmin,zmax;
+
+ snew(rest_at->index,state->natoms);
+
+ xmin=xmax=state->x[ins_at->index[0]][XX];
+ ymin=ymax=state->x[ins_at->index[0]][YY];
+ zmin=zmax=state->x[ins_at->index[0]][ZZ];
+
+ for(i=0;i<state->natoms;i++)
+ {
+ gid = groups->grpnr[egcFREEZE][i];
+ if(groups->grps[egcFREEZE].nm_ind[gid]==ins_grp_id)
+ {
+ x=state->x[i][XX];
+ if (x<xmin) xmin=x;
+ if (x>xmax) xmax=x;
+ y=state->x[i][YY];
+ if (y<ymin) ymin=y;
+ if (y>ymax) ymax=y;
+ z=state->x[i][ZZ];
+ if (z<zmin) zmin=z;
+ if (z>zmax) zmax=z;
+ } else {
+ rest_at->index[c]=i;
+ c++;
+ }
+ }
+
+ rest_at->nr=c;
+ srenew(rest_at->index,c);
+
+ if(xy_max>1.000001)
+ {
+ pos_ins->xmin[XX]=xmin-((xmax-xmin)*xy_max-(xmax-xmin))/2;
+ pos_ins->xmin[YY]=ymin-((ymax-ymin)*xy_max-(ymax-ymin))/2;
+
+ pos_ins->xmax[XX]=xmax+((xmax-xmin)*xy_max-(xmax-xmin))/2;
+ pos_ins->xmax[YY]=ymax+((ymax-ymin)*xy_max-(ymax-ymin))/2;
+ } else {
+ pos_ins->xmin[XX]=xmin;
+ pos_ins->xmin[YY]=ymin;
+
+ pos_ins->xmax[XX]=xmax;
+ pos_ins->xmax[YY]=ymax;
+ }
+
+ /* 6.0 is estimated thickness of bilayer */
+ if( (zmax-zmin) < 6.0 )
+ {
+ pos_ins->xmin[ZZ]=zmin+(zmax-zmin)/2.0-3.0;
+ pos_ins->xmax[ZZ]=zmin+(zmax-zmin)/2.0+3.0;
+ } else {
+ pos_ins->xmin[ZZ]=zmin;
+ pos_ins->xmax[ZZ]=zmax;
+ }
+
+ return c;
+}
+
+real est_prot_area(pos_ins_t *pos_ins,rvec *r,t_block *ins_at, mem_t *mem_p)
+{
+ real x,y,dx=0.15,dy=0.15,area=0.0;
+ real add;
+ int c,at;
+
+ for(x=pos_ins->xmin[XX];x<pos_ins->xmax[XX];x+=dx)
+ {
+ for(y=pos_ins->xmin[YY];y<pos_ins->xmax[YY];y+=dy)
+ {
+ c=0;
+ add=0.0;
+ do
+ {
+ at=ins_at->index[c];
+ if ( (r[at][XX]>=x) && (r[at][XX]<x+dx) &&
+ (r[at][YY]>=y) && (r[at][YY]<y+dy) &&
+ (r[at][ZZ]>mem_p->zmin+1.0) && (r[at][ZZ]<mem_p->zmax-1.0) )
+ add=1.0;
+ c++;
+ } while ( (c<ins_at->nr) && (add<0.5) );
+ area+=add;
+ }
+ }
+ area=area*dx*dy;
+
+ return area;
+}
+
+void init_lip(matrix box, gmx_mtop_t *mtop, lip_t *lip)
+{
+ int i;
+ real mem_area;
+ int mol1=0;
+
+ mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
+ for(i=0;i<mtop->nmolblock;i++)
+ {
+ if(mtop->molblock[i].type == lip->id)
+ {
+ lip->nr=mtop->molblock[i].nmol;
+ lip->natoms=mtop->molblock[i].natoms_mol;
+ }
+ }
+ lip->area=2.0*mem_area/(double)lip->nr;
+
+ for (i=0;i<lip->id;i++)
+ mol1+=mtop->molblock[i].nmol;
+ lip->mol1=mol1;
+}
+
+int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_ins_t *pos_ins)
+{
+ int i,j,at,mol,nmol,nmolbox,count;
+ t_block *mem_a;
+ real z,zmin,zmax,mem_area;
+ gmx_bool bNew;
+ atom_id *mol_id;
+ int type=0,block=0;
+
+ nmol=count=0;
+ mem_a=&(mem_p->mem_at);
+ snew(mol_id,mem_a->nr);
+/* snew(index,mem_a->nr); */
+ zmin=pos_ins->xmax[ZZ];
+ zmax=pos_ins->xmin[ZZ];
+ for(i=0;i<mem_a->nr;i++)
+ {
+ at=mem_a->index[i];
+ if( (r[at][XX]>pos_ins->xmin[XX]) && (r[at][XX]<pos_ins->xmax[XX]) &&
+ (r[at][YY]>pos_ins->xmin[YY]) && (r[at][YY]<pos_ins->xmax[YY]) &&
+ (r[at][ZZ]>pos_ins->xmin[ZZ]) && (r[at][ZZ]<pos_ins->xmax[ZZ]) )
+ {
+ mol = get_mol_id(at,mtop->nmolblock,mtop->molblock,&type,&block);
+
+ bNew=TRUE;
+ for(j=0;j<nmol;j++)
+ if(mol == mol_id[j])
+ bNew=FALSE;
+
+ if(bNew)
+ {
+ mol_id[nmol]=mol;
+ nmol++;
+ }
+
+ z=r[at][ZZ];
+ if(z<zmin) zmin=z;
+ if(z>zmax) zmax=z;
+
+/* index[count]=at;*/
+ count++;
+ }
+ }
+
+ mem_p->nmol=nmol;
+ srenew(mol_id,nmol);
+ mem_p->mol_id=mol_id;
+/* srenew(index,count);*/
+/* mem_p->mem_at.nr=count;*/
+/* sfree(mem_p->mem_at.index);*/
+/* mem_p->mem_at.index=index;*/
+
+ if((zmax-zmin)>(box[ZZ][ZZ]-0.5))
+ gmx_fatal(FARGS,"Something is wrong with your membrane. Max and min z values are %f and %f.\n"
+ "Maybe your membrane is not centered in the box, but located at the box edge in the z-direction,\n"
+ "so that one membrane is distributed over two periodic box images. Another possibility is that\n"
+ "your water layer is not thick enough.\n",zmax,zmin);
+ mem_p->zmin=zmin;
+ mem_p->zmax=zmax;
+ mem_p->zmed=(zmax-zmin)/2+zmin;
+
+ /*number of membrane molecules in protein box*/
+ nmolbox = count/mtop->molblock[block].natoms_mol;
+ /*mem_area = box[XX][XX]*box[YY][YY]-box[XX][YY]*box[YY][XX];
+ mem_p->lip_area = 2.0*mem_area/(double)mem_p->nmol;*/
+ mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
+ mem_p->lip_area = 2.0*mem_area/(double)nmolbox;
+
+ return mem_p->mem_at.nr;
+}
+
+void init_resize(t_block *ins_at,rvec *r_ins,pos_ins_t *pos_ins,mem_t *mem_p,rvec *r, gmx_bool bALLOW_ASYMMETRY)
+{
+ int i,j,at,c,outsidesum,gctr=0;
+ int idxsum=0;
+
+ /*sanity check*/
+ for (i=0;i<pos_ins->pieces;i++)
+ idxsum+=pos_ins->nidx[i];
+ if (idxsum!=ins_at->nr)
+ gmx_fatal(FARGS,"Piecewise sum of inserted atoms not same as size of group selected to insert.");
+
+ snew(pos_ins->geom_cent,pos_ins->pieces);
+ for (i=0;i<pos_ins->pieces;i++)
+ {
+ c=0;
+ outsidesum=0;
+ for(j=0;j<DIM;j++)
+ pos_ins->geom_cent[i][j]=0;
+
+ for(j=0;j<DIM;j++)
+ pos_ins->geom_cent[i][j]=0;
+ for (j=0;j<pos_ins->nidx[i];j++)
+ {
+ at=pos_ins->subindex[i][j];
+ copy_rvec(r[at],r_ins[gctr]);
+ if( (r_ins[gctr][ZZ]<mem_p->zmax) && (r_ins[gctr][ZZ]>mem_p->zmin) )
+ {
+ rvec_inc(pos_ins->geom_cent[i],r_ins[gctr]);
+ c++;
+ }
+ else
+ outsidesum++;
+ gctr++;
+ }
+ if (c>0)
+ svmul(1/(double)c,pos_ins->geom_cent[i],pos_ins->geom_cent[i]);
+ if (!bALLOW_ASYMMETRY)
+ pos_ins->geom_cent[i][ZZ]=mem_p->zmed;
+
+ fprintf(stderr,"Embedding piece %d with center of geometry: %f %f %f\n",i,pos_ins->geom_cent[i][XX],pos_ins->geom_cent[i][YY],pos_ins->geom_cent[i][ZZ]);
+ }
+ fprintf(stderr,"\n");
+}
+
+void resize(rvec *r_ins, rvec *r, pos_ins_t *pos_ins,rvec fac)
+{
+ int i,j,k,at,c=0;
+ for (k=0;k<pos_ins->pieces;k++)
+ for(i=0;i<pos_ins->nidx[k];i++)
+ {
+ at=pos_ins->subindex[k][i];
+ for(j=0;j<DIM;j++)
+ r[at][j]=pos_ins->geom_cent[k][j]+fac[j]*(r_ins[c][j]-pos_ins->geom_cent[k][j]);
+ c++;
+ }
+}
+
+int gen_rm_list(rm_t *rm_p,t_block *ins_at,t_block *rest_at,t_pbc *pbc, gmx_mtop_t *mtop,
+ rvec *r, rvec *r_ins, mem_t *mem_p, pos_ins_t *pos_ins, real probe_rad, int low_up_rm, gmx_bool bALLOW_ASYMMETRY)
+{
+ int i,j,k,l,at,at2,mol_id;
+ int type=0,block=0;
+ int nrm,nupper,nlower;
+ real r_min_rad,z_lip,min_norm;
+ gmx_bool bRM;
+ rvec dr,dr_tmp;
+ real *dist;
+ int *order;
+
+ r_min_rad=probe_rad*probe_rad;
+ snew(rm_p->mol,mtop->mols.nr);
+ snew(rm_p->block,mtop->mols.nr);
+ nrm=nupper=0;
+ nlower=low_up_rm;
+ for(i=0;i<ins_at->nr;i++)
+ {
+ at=ins_at->index[i];
+ for(j=0;j<rest_at->nr;j++)
+ {
+ at2=rest_at->index[j];
+ pbc_dx(pbc,r[at],r[at2],dr);
+
+ if(norm2(dr)<r_min_rad)
+ {
+ mol_id = get_mol_id(at2,mtop->nmolblock,mtop->molblock,&type,&block);
+ bRM=TRUE;
+ for(l=0;l<nrm;l++)
+ if(rm_p->mol[l]==mol_id)
+ bRM=FALSE;
+ if(bRM)
+ {
+ /*fprintf(stderr,"%d wordt toegevoegd\n",mol_id);*/
+ rm_p->mol[nrm]=mol_id;
+ rm_p->block[nrm]=block;
+ nrm++;
+ z_lip=0.0;
+ for(l=0;l<mem_p->nmol;l++)
+ {
+ if(mol_id==mem_p->mol_id[l])
+ {
+ for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
+ z_lip+=r[k][ZZ];
+ z_lip/=mtop->molblock[block].natoms_mol;
+ if(z_lip<mem_p->zmed)
+ nlower++;
+ else
+ nupper++;
+ }
+ }
+ }
+ }
+ }
+ }
+
+ /*make sure equal number of lipids from upper and lower layer are removed */
+ if( (nupper!=nlower) && (!bALLOW_ASYMMETRY) )
+ {
+ snew(dist,mem_p->nmol);
+ snew(order,mem_p->nmol);
+ for(i=0;i<mem_p->nmol;i++)
+ {
+ at = mtop->mols.index[mem_p->mol_id[i]];
+ pbc_dx(pbc,r[at],pos_ins->geom_cent[0],dr);
+ if (pos_ins->pieces>1)
+ {
+ /*minimum dr value*/
+ min_norm=norm2(dr);
+ for (k=1;k<pos_ins->pieces;k++)
+ {
+ pbc_dx(pbc,r[at],pos_ins->geom_cent[k],dr_tmp);
+ if (norm2(dr_tmp) < min_norm)
+ {
+ min_norm=norm2(dr_tmp);
+ copy_rvec(dr_tmp,dr);
+ }
+ }
+ }
+ dist[i]=dr[XX]*dr[XX]+dr[YY]*dr[YY];
+ j=i-1;
+ while (j>=0 && dist[i]<dist[order[j]])
+ {
+ order[j+1]=order[j];
+ j--;
+ }
+ order[j+1]=i;
+ }
+
+ i=0;
+ while(nupper!=nlower)
+ {
+ mol_id=mem_p->mol_id[order[i]];
+ block=get_block(mol_id,mtop->nmolblock,mtop->molblock);
+
+ bRM=TRUE;
+ for(l=0;l<nrm;l++)
+ if(rm_p->mol[l]==mol_id)
+ bRM=FALSE;
+ if(bRM)
+ {
+ z_lip=0;
+ for(k=mtop->mols.index[mol_id];k<mtop->mols.index[mol_id+1];k++)
+ z_lip+=r[k][ZZ];
+ z_lip/=mtop->molblock[block].natoms_mol;
+ if(nupper>nlower && z_lip<mem_p->zmed)
+ {
+ rm_p->mol[nrm]=mol_id;
+ rm_p->block[nrm]=block;
+ nrm++;
+ nlower++;
+ }
+ else if (nupper<nlower && z_lip>mem_p->zmed)
+ {
+ rm_p->mol[nrm]=mol_id;
+ rm_p->block[nrm]=block;
+ nrm++;
+ nupper++;
+ }
+ }
+ i++;
+
+ if(i>mem_p->nmol)
+ gmx_fatal(FARGS,"Trying to remove more lipid molecules than there are in the membrane");
+ }
+ sfree(dist);
+ sfree(order);
+ }
+
+ rm_p->nr=nrm;
+ srenew(rm_p->mol,nrm);
+ srenew(rm_p->block,nrm);
+
+ return nupper+nlower;
+}
+
+void rm_group(t_inputrec *ir, gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state *state, t_block *ins_at, pos_ins_t *pos_ins)
+{
+ int i,j,k,n,rm,mol_id,at,block;
+ rvec *x_tmp,*v_tmp;
+ atom_id *list,*new_mols;
+ unsigned char *new_egrp[egcNR];
+ gmx_bool bRM;
+
+ snew(list,state->natoms);
+ n=0;
+ for(i=0;i<rm_p->nr;i++)
+ {
+ mol_id=rm_p->mol[i];
+ at=mtop->mols.index[mol_id];
+ block =rm_p->block[i];
+ mtop->molblock[block].nmol--;
+ for(j=0;j<mtop->molblock[block].natoms_mol;j++)
+ {
+ list[n]=at+j;
+ n++;
+ }
+
+ mtop->mols.index[mol_id]=-1;
+ }
+
+ mtop->mols.nr-=rm_p->nr;
+ mtop->mols.nalloc_index-=rm_p->nr;
+ snew(new_mols,mtop->mols.nr);
+ for(i=0;i<mtop->mols.nr+rm_p->nr;i++)
+ {
+ j=0;
+ if(mtop->mols.index[i]!=-1)
+ {
+ new_mols[j]=mtop->mols.index[i];
+ j++;
+ }
+ }
+ sfree(mtop->mols.index);
+ mtop->mols.index=new_mols;
+
+
+ mtop->natoms-=n;
+ state->natoms-=n;
+ state->nalloc=state->natoms;
+ snew(x_tmp,state->nalloc);
+ snew(v_tmp,state->nalloc);
+
+ for(i=0;i<egcNR;i++)
+ {
+ if(groups->grpnr[i]!=NULL)
+ {
+ groups->ngrpnr[i]=state->natoms;
+ snew(new_egrp[i],state->natoms);
+ }
+ }
+
+ rm=0;
+ for (i=0;i<state->natoms+n;i++)
+ {
+ bRM=FALSE;
+ for(j=0;j<n;j++)
+ {
+ if(i==list[j])
+ {
+ bRM=TRUE;
+ rm++;
+ }
+ }
+
+ if(!bRM)
+ {
+ for(j=0;j<egcNR;j++)
+ {
+ if(groups->grpnr[j]!=NULL)
+ {
+ new_egrp[j][i-rm]=groups->grpnr[j][i];
+ }
+ }
+ copy_rvec(state->x[i],x_tmp[i-rm]);
+ copy_rvec(state->v[i],v_tmp[i-rm]);
+ for(j=0;j<ins_at->nr;j++)
+ {
+ if (i==ins_at->index[j])
+ ins_at->index[j]=i-rm;
+ }
+ for(j=0;j<pos_ins->pieces;j++)
+ {
+ for(k=0;k<pos_ins->nidx[j];k++)
+ {
+ if (i==pos_ins->subindex[j][k])
+ pos_ins->subindex[j][k]=i-rm;
+ }
+ }
+ }
+ }
+ sfree(state->x);
+ state->x=x_tmp;
+ sfree(state->v);
+ state->v=v_tmp;
+
+ for(i=0;i<egcNR;i++)
+ {
+ if(groups->grpnr[i]!=NULL)
+ {
+ sfree(groups->grpnr[i]);
+ groups->grpnr[i]=new_egrp[i];
+ }
+ }
+}
+
+int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
+{
+ int i,j,m;
+ int type,natom,nmol,at,atom1=0,rm_at=0;
+ gmx_bool *bRM,bINS;
+ /*this routine lives dangerously by assuming that all molecules of a given type are in order in the structure*/
+ /*this routine does not live as dangerously as it seems. There is namely a check in mdrunner_membed to make
+ *sure that g_membed exits with a warning when there are molecules of the same type not in the
+ *ins_at index group. MGWolf 050710 */
+
+
+ snew(bRM,mtop->nmoltype);
+ for (i=0;i<mtop->nmoltype;i++)
+ {
+ bRM[i]=TRUE;
+ }
+
+ for (i=0;i<mtop->nmolblock;i++)
+ {
+ /*loop over molecule blocks*/
+ type =mtop->molblock[i].type;
+ natom =mtop->molblock[i].natoms_mol;
+ nmol =mtop->molblock[i].nmol;
+
+ for(j=0;j<natom*nmol && bRM[type]==TRUE;j++)
+ {
+ /*loop over atoms in the block*/
+ at=j+atom1; /*atom index = block index + offset*/
+ bINS=FALSE;
+
+ for (m=0;(m<ins_at->nr) && (bINS==FALSE);m++)
+ {
+ /*loop over atoms in insertion index group to determine if we're inserting one*/
+ if(at==ins_at->index[m])
+ {
+ bINS=TRUE;
+ }
+ }
+ bRM[type]=bINS;
+ }
+ atom1+=natom*nmol; /*update offset*/
+ if(bRM[type])
+ {
+ rm_at+=natom*nmol; /*increment bonded removal counter by # atoms in block*/
+ }
+ }
+
+ for(i=0;i<mtop->nmoltype;i++)
+ {
+ if(bRM[i])
+ {
+ for(j=0;j<F_LJ;j++)
+ {
+ mtop->moltype[i].ilist[j].nr=0;
+ }
+ for(j=F_POSRES;j<=F_VSITEN;j++)
+ {
+ mtop->moltype[i].ilist[j].nr=0;
+ }
+ }
+ }
+ sfree(bRM);
+
+ return rm_at;
+}
+
+void top_update(const char *topfile, char *ins, rm_t *rm_p, gmx_mtop_t *mtop)
+{
+#define TEMP_FILENM "temp.top"
+ int bMolecules=0;
+ FILE *fpin,*fpout;
+ char buf[STRLEN],buf2[STRLEN],*temp;
+ int i,*nmol_rm,nmol,line;
+
+ fpin = ffopen(topfile,"r");
+ fpout = ffopen(TEMP_FILENM,"w");
+
+ snew(nmol_rm,mtop->nmoltype);
+ for(i=0;i<rm_p->nr;i++)
+ nmol_rm[rm_p->block[i]]++;
+
+ line=0;
+ while(fgets(buf,STRLEN,fpin))
+ {
+ line++;
+ if(buf[0]!=';')
+ {
+ strcpy(buf2,buf);
+ if ((temp=strchr(buf2,'\n')) != NULL)
+ temp[0]='\0';
+ ltrim(buf2);
+
+ if (buf2[0]=='[')
+ {
+ buf2[0]=' ';
+ if ((temp=strchr(buf2,'\n')) != NULL)
+ temp[0]='\0';
+ rtrim(buf2);
+ if (buf2[strlen(buf2)-1]==']')
+ {
+ buf2[strlen(buf2)-1]='\0';
+ ltrim(buf2);
+ rtrim(buf2);
+ if (gmx_strcasecmp(buf2,"molecules")==0)
+ bMolecules=1;
+ }
+ fprintf(fpout,"%s",buf);
+ } else if (bMolecules==1)
+ {
+ for(i=0;i<mtop->nmolblock;i++)
+ {
+ nmol=mtop->molblock[i].nmol;
+ sprintf(buf,"%-15s %5d\n",*(mtop->moltype[mtop->molblock[i].type].name),nmol);
+ fprintf(fpout,"%s",buf);
+ }
+ bMolecules=2;
+ } else if (bMolecules==2)
+ {
+ /* print nothing */
+ } else
+ {
+ fprintf(fpout,"%s",buf);
+ }
+ } else
+ {
+ fprintf(fpout,"%s",buf);
+ }
+ }
+
+ fclose(fpout);
+ /* use ffopen to generate backup of topinout */
+ fpout=ffopen(topfile,"w");
+ fclose(fpout);
+ rename(TEMP_FILENM,topfile);
+#undef TEMP_FILENM
+}
+
+void rescale_membed(int step_rel, gmx_membed_t *membed, rvec *x)
+{
+ /* Set new positions for the group to embed */
+ if(step_rel<=membed->it_xy)
+ {
+ membed->fac[0]+=membed->xy_step;
+ membed->fac[1]+=membed->xy_step;
+ } else if (step_rel<=(membed->it_xy+membed->it_z))
+ {
+ membed->fac[2]+=membed->z_step;
+ }
+ resize(membed->r_ins,x,membed->pos_ins,membed->fac);
+}
+
+void init_membed(FILE *fplog, gmx_membed_t *membed, int nfile, const t_filenm fnm[], gmx_mtop_t *mtop, t_inputrec *inputrec, t_state *state, t_commrec *cr,real *cpt)
+{
+ char *ins;
+ int i,rm_bonded_at,fr_id,fr_i=0,tmp_id,warn=0;
+ int ng,j,max_lip_rm,ins_grp_id,ins_nat,mem_nat,ntype,lip_rm,tpr_version;
+ real prot_area;
+ rvec *r_ins=NULL;
+ t_block *ins_at,*rest_at;
+ pos_ins_t *pos_ins;
+ mem_t *mem_p;
+ rm_t *rm_p;
+ gmx_groups_t *groups;
+ gmx_bool bExcl=FALSE;
+ t_atoms atoms;
+ t_pbc *pbc;
+ char **piecename=NULL;
+
+ /* input variables */
+ const char *membed_input;
+ real xy_fac = 0.5;
+ real xy_max = 1.0;
+ real z_fac = 1.0;
+ real z_max = 1.0;
+ int it_xy = 1000;
+ int it_z = 0;
+ real probe_rad = 0.22;
+ int low_up_rm = 0;
+ int maxwarn=0;
+ int pieces=1;
+ gmx_bool bALLOW_ASYMMETRY=FALSE;
+
+ snew(ins_at,1);
+ snew(pos_ins,1);
+
+ if(MASTER(cr))
+ {
+ /* get input data out membed file */
+ membed_input = opt2fn("-membed",nfile,fnm);
+ get_input(membed_input,&xy_fac,&xy_max,&z_fac,&z_max,&it_xy,&it_z,&probe_rad,&low_up_rm,&maxwarn,&pieces,&bALLOW_ASYMMETRY);
+
+ tpr_version = get_tpr_version(ftp2fn(efTPX,nfile,fnm));
+ if (tpr_version<58)
+ gmx_fatal(FARGS,"Version of *.tpr file to old (%d). Rerun grompp with gromacs VERSION 4.0.3 or newer.\n",tpr_version);
+
+ if( !EI_DYNAMICS(inputrec->eI) )
+ gmx_input("Change integrator to a dynamics integrator in mdp file (e.g. md or sd).");
+
+ if(PAR(cr))
+ gmx_input("Sorry, parallel g_membed is not yet fully functional.");
+
+#ifdef GMX_OPENMM
+ gmx_input("Sorry, g_membed does not work with openmm.");
+#endif
+
+ if(*cpt>=0)
+ {
+ fprintf(stderr,"\nSetting -cpt to -1, because embedding cannot be restarted from cpt-files.\n");
+ *cpt=-1;
+ }
+ groups=&(mtop->groups);
+
+ atoms=gmx_mtop_global_atoms(mtop);
+ snew(mem_p,1);
+ fprintf(stderr,"\nSelect a group to embed in the membrane:\n");
+ get_index(&atoms,opt2fn_null("-mn",nfile,fnm),1,&(ins_at->nr),&(ins_at->index),&ins);
+ ins_grp_id = search_string(ins,groups->ngrpname,(groups->grpname));
+ fprintf(stderr,"\nSelect a group to embed %s into (e.g. the membrane):\n",ins);
+ get_index(&atoms,opt2fn_null("-mn",nfile,fnm),1,&(mem_p->mem_at.nr),&(mem_p->mem_at.index),&(mem_p->name));
+
+ pos_ins->pieces=pieces;
+ snew(pos_ins->nidx,pieces);
+ snew(pos_ins->subindex,pieces);
+ snew(piecename,pieces);
+ if (pieces>1)
+ {
+ fprintf(stderr,"\nSelect pieces to embed:\n");
+ get_index(&atoms,opt2fn_null("-mn",nfile,fnm),pieces,pos_ins->nidx,pos_ins->subindex,piecename);
+ }
+ else
+ {
+ /*use whole embedded group*/
+ snew(pos_ins->nidx,1);
+ snew(pos_ins->subindex,1);
+ pos_ins->nidx[0]=ins_at->nr;
+ pos_ins->subindex[0]=ins_at->index;
+ }
+
+ if(probe_rad<0.2199999)
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d:\nA probe radius (-rad) smaller than 0.2 can result in overlap between waters "
+ "and the group to embed, which will result in Lincs errors etc.\nIf you are sure, you can increase maxwarn.\n\n",warn);
+ }
+
+ if(xy_fac<0.09999999)
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d:\nThe initial size of %s is probably too smal.\n"
+ "If you are sure, you can increase maxwarn.\n\n",warn,ins);
+ }
+
+ if(it_xy<1000)
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the xy-coordinates of %s (%d) is probably too small.\n"
+ "Increase -nxy or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_xy);
+ }
+
+ if( (it_z<100) && ( z_fac<0.99999999 || z_fac>1.0000001) )
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d;\nThe number of steps used to grow the z-coordinate of %s (%d) is probably too small.\n"
+ "Increase -nz or, if you are sure, you can increase maxwarn.\n\n",warn,ins,it_z);
+ }
+
+ if(it_xy+it_z>inputrec->nsteps)
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d:\nThe number of growth steps (-nxy + -nz) is larger than the number of steps in the tpr.\n"
+ "If you are sure, you can increase maxwarn.\n\n",warn);
+ }
+
+ fr_id=-1;
+ if( inputrec->opts.ngfrz==1)
+ gmx_fatal(FARGS,"You did not specify \"%s\" as a freezegroup.",ins);
+ for(i=0;i<inputrec->opts.ngfrz;i++)
+ {
+ tmp_id = mtop->groups.grps[egcFREEZE].nm_ind[i];
+ if(ins_grp_id==tmp_id)
+ {
+ fr_id=tmp_id;
+ fr_i=i;
+ }
+ }
+ if (fr_id == -1 )
+ gmx_fatal(FARGS,"\"%s\" not as freezegroup defined in the mdp-file.",ins);
+
+ for(i=0;i<DIM;i++)
+ if( inputrec->opts.nFreeze[fr_i][i] != 1)
+ gmx_fatal(FARGS,"freeze dimensions for %s are not Y Y Y\n",ins);
+
+ ng = groups->grps[egcENER].nr;
+ if (ng == 1)
+ gmx_input("No energy groups defined. This is necessary for energy exclusion in the freeze group");
+
+ for(i=0;i<ng;i++)
+ {
+ for(j=0;j<ng;j++)
+ {
+ if (inputrec->opts.egp_flags[ng*i+j] == EGP_EXCL)
+ {
+ bExcl = TRUE;
+ if ( (groups->grps[egcENER].nm_ind[i] != ins_grp_id) || (groups->grps[egcENER].nm_ind[j] != ins_grp_id) )
+ gmx_fatal(FARGS,"Energy exclusions \"%s\" and \"%s\" do not match the group to embed \"%s\"",
+ *groups->grpname[groups->grps[egcENER].nm_ind[i]],
+ *groups->grpname[groups->grps[egcENER].nm_ind[j]],ins);
+ }
+ }
+ }
+ if (!bExcl)
+ gmx_input("No energy exclusion groups defined. This is necessary for energy exclusion in the freeze group");
+
+ /* Guess the area the protein will occupy in the membrane plane Calculate area per lipid*/
+ snew(rest_at,1);
+ ins_nat = init_ins_at(ins_at,rest_at,state,pos_ins,groups,ins_grp_id,xy_max);
+ /* Check moleculetypes in insertion group */
+ check_types(ins_at,rest_at,mtop);
+
+ mem_nat = init_mem_at(mem_p,mtop,state->x,state->box,pos_ins);
+
+ prot_area = est_prot_area(pos_ins,state->x,ins_at,mem_p);
+ if ( (prot_area>7.5) && ( (state->box[XX][XX]*state->box[YY][YY]-state->box[XX][YY]*state->box[YY][XX])<50) )
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d:\nThe xy-area is very small compared to the area of the protein.\n"
+ "This might cause pressure problems during the growth phase. Just try with\n"
+ "current setup (-maxwarn + 1), but if pressure problems occur, lower the\n"
+ "compressibility in the mdp-file or use no pressure coupling at all.\n\n",warn);
+ }
+ if(warn>maxwarn)
+ gmx_fatal(FARGS,"Too many warnings.\n");
+
+ printf("The estimated area of the protein in the membrane is %.3f nm^2\n",prot_area);
+ printf("\nThere are %d lipids in the membrane part that overlaps the protein.\nThe area per lipid is %.4f nm^2.\n",mem_p->nmol,mem_p->lip_area);
+
+ /* Maximum number of lipids to be removed*/
+ max_lip_rm=(int)(2*prot_area/mem_p->lip_area);
+ printf("Maximum number of lipids that will be removed is %d.\n",max_lip_rm);
+
+ printf("\nWill resize the protein by a factor of %.3f in the xy plane and %.3f in the z direction.\n"
+ "This resizing will be done with respect to the geometrical center of all protein atoms\n"
+ "that span the membrane region, i.e. z between %.3f and %.3f\n\n",xy_fac,z_fac,mem_p->zmin,mem_p->zmax);
+
+ /* resize the protein by xy and by z if necessary*/
+ snew(r_ins,ins_at->nr);
+ init_resize(ins_at,r_ins,pos_ins,mem_p,state->x,bALLOW_ASYMMETRY);
+ membed->fac[0]=membed->fac[1]=xy_fac;
+ membed->fac[2]=z_fac;
+
+ membed->xy_step =(xy_max-xy_fac)/(double)(it_xy);
+ membed->z_step =(z_max-z_fac)/(double)(it_z-1);
+
+ resize(r_ins,state->x,pos_ins,membed->fac);
+
+ /* remove overlapping lipids and water from the membrane box*/
+ /*mark molecules to be removed*/
+ snew(pbc,1);
+ set_pbc(pbc,inputrec->ePBC,state->box);
+
+ snew(rm_p,1);
+ lip_rm = gen_rm_list(rm_p,ins_at,rest_at,pbc,mtop,state->x, r_ins, mem_p,pos_ins,probe_rad,low_up_rm,bALLOW_ASYMMETRY);
+ lip_rm -= low_up_rm;
+
+ if(fplog)
+ for(i=0;i<rm_p->nr;i++)
+ fprintf(fplog,"rm mol %d\n",rm_p->mol[i]);
+
+ for(i=0;i<mtop->nmolblock;i++)
+ {
+ ntype=0;
+ for(j=0;j<rm_p->nr;j++)
+ if(rm_p->block[j]==i)
+ ntype++;
+ printf("Will remove %d %s molecules\n",ntype,*(mtop->moltype[mtop->molblock[i].type].name));
+ }
+
+ if(lip_rm>max_lip_rm)
+ {
+ warn++;
+ fprintf(stderr,"\nWarning %d:\nTrying to remove a larger lipid area than the estimated protein area\n"
+ "Try making the -xyinit resize factor smaller.\n\n",warn);
+ }
+
+ /*remove all lipids and waters overlapping and update all important structures*/
+ rm_group(inputrec,groups,mtop,rm_p,state,ins_at,pos_ins);
+
+ rm_bonded_at = rm_bonded(ins_at,mtop);
+ if (rm_bonded_at != ins_at->nr)
+ {
+ fprintf(stderr,"Warning: The number of atoms for which the bonded interactions are removed is %d, "
+ "while %d atoms are embedded. Make sure that the atoms to be embedded are not in the same"
+ "molecule type as atoms that are not to be embedded.\n",rm_bonded_at,ins_at->nr);
+ }
+
+ if(warn>maxwarn)
+ gmx_fatal(FARGS,"Too many warnings.\nIf you are sure these warnings are harmless, you can increase -maxwarn");
+
+ if (ftp2bSet(efTOP,nfile,fnm))
+ top_update(opt2fn("-p",nfile,fnm),ins,rm_p,mtop);
+
+ sfree(pbc);
+ sfree(rest_at);
+ if (pieces>1) { sfree(piecename); }
+
+ membed->it_xy=it_xy;
+ membed->it_z=it_z;
+ membed->pos_ins=pos_ins;
+ membed->r_ins=r_ins;
+ }
+}
#include "checkpoint.h"
#include "mtop_util.h"
#include "sighandler.h"
+#include "gmx_membed.h"
#ifdef GMX_LIB_MPI
#include <mpi.h>
t_mdatoms *mdatoms,
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,t_forcerec *fr,
- int repl_ex_nst,int repl_ex_seed,
+ int repl_ex_nst,int repl_ex_seed,gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
}
/* ####### END SET VARIABLES FOR NEXT ITERATION ###### */
+
+ if ( (membed!=NULL) && (!bLastStep) )
+ rescale_membed(step_rel,membed,state_global->x);
if (bRerunMD)
{
#include "genborn.h"
#include "string2.h"
#include "copyrite.h"
+#include "gmx_membed.h"
#ifdef GMX_THREADS
#include "tmpi.h"
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
t_nrnb *nrnb,gmx_wallcycle_t wcycle,
gmx_edsam_t ed,t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
"appropriate options have been given. Currently under",
"investigation are: polarizability, and X-Ray bombardments.",
"[PAR]",
+ "The option [TT]-membed[dd] does what used to be g_membed, i.e. embed",
+ "a protein into a membrane. The data file should contain the options",
+ "that where passed to g_membed before. The [TT]-mn[tt] and [TT]-mp[tt]",
+ "both apply to this as well.",
+ "[PAR]",
"The option [TT]-pforce[tt] is useful when you suspect a simulation",
"crashes due to too large forces. With this option coordinates and",
"forces of atoms with a force larger than a certain value will",
{ efXVG, "-px", "pullx", ffOPTWR },
{ efXVG, "-pf", "pullf", ffOPTWR },
{ efMTX, "-mtx", "nm", ffOPTWR },
- { efNDX, "-dn", "dipole", ffOPTWR }
+ { efNDX, "-dn", "dipole", ffOPTWR },
+ { efDAT, "-membed", "membed", ffOPTRD },
+ { efTOP, "-mp", "membed", ffOPTRD },
+ { efNDX, "-mn", "membed", ffOPTRD }
};
#define NFILE asize(fnm)
#include "sighandler.h"
#include "tpxio.h"
#include "txtdump.h"
+#include "gmx_membed.h"
#include "md_openmm.h"
gmx_edsam_t ed=NULL;
t_commrec *cr_old=cr;
int nthreads=1;
+ gmx_membed_t *membed=NULL;
/* CAUTION: threads may be started later on in this function, so
cr doesn't reflect the final parallel state right now */
}
/* END OF CAUTION: cr is now reliable */
+ /* g_membed initialisation *
+ * Because we change the mtop, init_membed is called before the init_parallel *
+ * (in case we ever want to make it run in parallel) */
+ if (opt2bSet("-membed",nfile,fnm))
+ {
+ fprintf(stderr,"Entering membed code");
+ snew(membed,1);
+ init_membed(fplog,membed,nfile,fnm,mtop,inputrec,state,cr,&cpt_period);
+ }
+
if (PAR(cr))
{
/* now broadcast everything to the non-master nodes/threads: */
fcd,state,
mdatoms,nrnb,wcycle,ed,fr,
repl_ex_nst,repl_ex_seed,
+ membed,
cpt_period,max_hours,
deviceOptions,
Flags,
finish_run(fplog,cr,ftp2fn(efSTO,nfile,fnm),
inputrec,nrnb,wcycle,&runtime,
EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
+
+ if (opt2bSet("-membed",nfile,fnm))
+ {
+ sfree(membed);
+ }
/* Does what it says */
print_date_and_time(fplog,cr->nodeid,"Finished mdrun",&runtime);
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
#include "mtop_util.h"
#include "gmxfio.h"
#include "pme.h"
+#include "gmx_membed.h"
typedef struct {
t_state s;
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
{
for(m=0; (m<DIM); m++)
{
- b_inv[m] = divide(1.0,box_size[m]);
+ b_inv[m] = divide_err(1.0,box_size[m]);
}
bTriclinic = TRICLINIC(box);
}
gmx_edsam_t ed,
t_forcerec *fr,
int repl_ex_nst,int repl_ex_seed,
+ gmx_membed_t *membed,
real cpt_period,real max_hours,
const char *deviceOptions,
unsigned long Flags,
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
.deps
.libs
+g_highway
+g_logo
+g_xrama
+ngmx
+++ /dev/null
-Makefile
-Makefile.in
-.deps
-.libs
\ No newline at end of file
.deps
.libs
+do_dssp
+editconf
+eneconv
+g_anadock
+g_anaeig
+g_analyze
+g_angle
+g_bar
+g_bond
+g_bundle
+g_chi
+g_cluster
+g_clustsize
+g_confrms
+g_covar
+g_current
+g_density
+g_densmap
+g_dielectric
+g_dih
+g_dipoles
+g_disre
+g_dist
+g_dyndom
+g_enemat
+g_energy
+g_filter
+g_gyrate
+g_h2order
+g_hbond
+g_helix
+g_helixorient
+g_kinetics
+g_lie
+g_mdmat
+g_membed
+g_mindist
+g_morph
+g_msd
+g_nmeig
+g_nmens
+g_nmtraj
+g_order
+g_polystat
+g_potential
+g_principal
+g_rama
+g_rdf
+g_rms
+g_rmsdist
+g_rmsf
+g_rotacf
+g_rotmat
+g_saltbr
+g_sas
+g_sdf
+g_select
+g_sgangle
+g_sham
+g_sigeps
+g_sorient
+g_spatial
+g_spol
+g_tcaf
+g_traj
+g_tune_pme
+g_vanhove
+g_velacc
+g_wham
+g_wheel
+genbox
+genconf
+genion
+genrestr
+make_edi
+make_ndx
+mk_angndx
+trjcat
+trjconv
+trjorder
+xpm2ps
gmx_trjconv.c gmx_trjcat.c gmx_trjorder.c gmx_xpm2ps.c \
gmx_editconf.c gmx_genbox.c gmx_genion.c gmx_genconf.c \
gmx_genpr.c gmx_eneconv.c gmx_vanhove.c gmx_wheel.c \
- addconf.c addconf.h gmx_tune_pme.c gmx_membed.c \
+ addconf.c addconf.h gmx_tune_pme.c \
calcpot.c calcpot.h edittop.c
bin_PROGRAMS = \
g_sham g_sorient g_spol \
g_spatial g_pme_error \
g_tcaf g_traj g_tune_pme \
- g_vanhove g_velacc g_membed \
+ g_vanhove g_velacc \
g_clustsize g_mdmat g_wham \
g_sigeps
#include <signal.h>
#include <stdlib.h>
#include "typedefs.h"
-#include "smalloc.h"
#include "sysstuff.h"
-#include "vec.h"
#include "statutil.h"
#include "macros.h"
#include "copyrite.h"
#include "main.h"
-#include "futil.h"
-#include "edsam.h"
-#include "checkpoint.h"
-#include "vcm.h"
-#include "mdebin.h"
-#include "nrnb.h"
-#include "calcmu.h"
-#include "index.h"
-#include "vsite.h"
-#include "update.h"
-#include "ns.h"
-#include "trnio.h"
-#include "xtcio.h"
-#include "mdrun.h"
-#include "confio.h"
-#include "network.h"
-#include "pull.h"
-#include "xvgr.h"
-#include "physics.h"
-#include "names.h"
-#include "disre.h"
-#include "orires.h"
-#include "dihre.h"
-#include "pppm.h"
-#include "pme.h"
-#include "mdatoms.h"
-#include "qmmm.h"
-#include "mpelogging.h"
-#include "domdec.h"
-#include "partdec.h"
-#include "topsort.h"
-#include "coulomb.h"
-#include "constr.h"
-#include "shellfc.h"
-#include "mvdata.h"
-#include "checkpoint.h"
-#include "mtop_util.h"
-#include "tpxio.h"
-#include "string2.h"
-#include "sighandler.h"
#include "gmx_ana.h"
#ifdef GMX_LIB_MPI
return rc;
}
+
int gmx_membed(int argc,char *argv[])
{
const char *desc[] = {
" - It is recommended to perform a short equilibration run after the embedding",
"(see Wolf et al, J Comp Chem 31 (2010) 2169-2174, to re-equilibrate the membrane. Clearly",
"protein equilibration might require longer.\n",
+ " - It is now also possible to use the g_membed functionality with mdrun. You should than pass",
+ "a data file containing the command line options of g_membed following the -membed option, for",
+ "example mdrun -s into_mem.tpr -membed membed.dat.",
"\n"
};
- t_commrec *cr;
t_filenm fnm[] = {
{ efTPX, "-f", "into_mem", ffREAD },
{ efNDX, "-n", "index", ffOPTRD },
{ efTOP, "-p", "topol", ffOPTRW },
{ efTRN, "-o", NULL, ffWRITE },
{ efXTC, "-x", NULL, ffOPTWR },
- { efCPT, "-cpi", NULL, ffOPTRD },
- { efCPT, "-cpo", NULL, ffOPTWR },
{ efSTO, "-c", "membedded", ffWRITE },
{ efEDR, "-e", "ener", ffWRITE },
- { efLOG, "-g", "md", ffWRITE },
- { efEDI, "-ei", "sam", ffOPTRD },
- { efTRX, "-rerun", "rerun", ffOPTRD },
- { efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tablep", "tablep", ffOPTRD },
- { efXVG, "-tableb", "table", ffOPTRD },
- { efXVG, "-dhdl", "dhdl", ffOPTWR },
- { efXVG, "-field", "field", ffOPTWR },
- { efXVG, "-table", "table", ffOPTRD },
- { efXVG, "-tablep", "tablep", ffOPTRD },
- { efXVG, "-tableb", "table", ffOPTRD },
- { efTRX, "-rerun", "rerun", ffOPTRD },
- { efXVG, "-tpi", "tpi", ffOPTWR },
- { efXVG, "-tpid", "tpidist", ffOPTWR },
- { efEDI, "-ei", "sam", ffOPTRD },
- { efEDO, "-eo", "sam", ffOPTWR },
- { efGCT, "-j", "wham", ffOPTRD },
- { efGCT, "-jo", "bam", ffOPTWR },
- { efXVG, "-ffout", "gct", ffOPTWR },
- { efXVG, "-devout", "deviatie", ffOPTWR },
- { efXVG, "-runav", "runaver", ffOPTWR },
- { efXVG, "-px", "pullx", ffOPTWR },
- { efXVG, "-pf", "pullf", ffOPTWR },
- { efMTX, "-mtx", "nm", ffOPTWR },
- { efNDX, "-dn", "dipole", ffOPTWR }
+ { efDAT, "-dat", "membed", ffWRITE }
};
#define NFILE asize(fnm)
/* Command line options ! */
- gmx_bool bCart = FALSE;
- gmx_bool bPPPME = FALSE;
- gmx_bool bPartDec = FALSE;
- gmx_bool bDDBondCheck = TRUE;
- gmx_bool bDDBondComm = TRUE;
- gmx_bool bVerbose = FALSE;
- gmx_bool bCompact = TRUE;
- gmx_bool bSepPot = FALSE;
- gmx_bool bRerunVSite = FALSE;
- gmx_bool bIonize = FALSE;
- gmx_bool bConfout = TRUE;
- gmx_bool bReproducible = FALSE;
-
- int npme=-1;
- int nmultisim=0;
- int nstglobalcomm=-1;
- int repl_ex_nst=0;
- int repl_ex_seed=-1;
- int nstepout=100;
- int nthreads=0; /* set to determine # of threads automatically */
- int resetstep=-1;
-
- rvec realddxyz={0,0,0};
- const char *ddno_opt[ddnoNR+1] =
- { NULL, "interleave", "pp_pme", "cartesian", NULL };
- const char *dddlb_opt[] =
- { NULL, "auto", "no", "yes", NULL };
- real rdd=0.0,rconstr=0.0,dlb_scale=0.8,pforce=-1;
- char *ddcsx=NULL,*ddcsy=NULL,*ddcsz=NULL;
- real cpt_period=15.0,max_hours=-1;
- gmx_bool bAppendFiles=TRUE,bAddPart=TRUE;
- gmx_bool bResetCountersHalfWay=FALSE;
- output_env_t oenv=NULL;
- const char *deviceOptions = "";
-
real xy_fac = 0.5;
real xy_max = 1.0;
real z_fac = 1.0;
int low_up_rm = 0;
int maxwarn=0;
int pieces=1;
- gmx_bool bALLOW_ASYMMETRY=FALSE;
-
+ gmx_bool bALLOW_ASYMMETRY=FALSE;
+ int nstepout=100;
+ gmx_bool bVerbose=FALSE;
+ char *mdrun_path=NULL;
/* arguments relevant to OPENMM only*/
#ifdef GMX_OPENMM
#endif
t_pargs pa[] = {
- { "-xyinit", FALSE, etREAL, {&xy_fac}, "Resize factor for the protein in the xy dimension before starting embedding" },
- { "-xyend", FALSE, etREAL, {&xy_max}, "Final resize factor in the xy dimension" },
- { "-zinit", FALSE, etREAL, {&z_fac}, "Resize factor for the protein in the z dimension before starting embedding" },
- { "-zend", FALSE, etREAL, {&z_max}, "Final resize faction in the z dimension" },
- { "-nxy", FALSE, etINT, {&it_xy}, "Number of iteration for the xy dimension" },
- { "-nz", FALSE, etINT, {&it_z}, "Number of iterations for the z dimension" },
- { "-rad", FALSE, etREAL, {&probe_rad}, "Probe radius to check for overlap between the group to embed and the membrane"},
- { "-pieces", FALSE, etINT, {&pieces}, "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
- { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY}, "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
- { "-ndiff" , FALSE, etINT, {&low_up_rm}, "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
- { "-maxwarn", FALSE, etINT, {&maxwarn}, "Maximum number of warning allowed" },
- { "-pd", FALSE, etBOOL,{&bPartDec},
- "HIDDENUse particle decompostion" },
- { "-dd", FALSE, etRVEC,{&realddxyz},
- "HIDDENDomain decomposition grid, 0 is optimize" },
- { "-nt", FALSE, etINT, {&nthreads},
- "HIDDENNumber of threads to start (0 is guess)" },
- { "-npme", FALSE, etINT, {&npme},
- "HIDDENNumber of separate nodes to be used for PME, -1 is guess" },
- { "-ddorder", FALSE, etENUM, {ddno_opt},
- "HIDDENDD node order" },
- { "-ddcheck", FALSE, etBOOL, {&bDDBondCheck},
- "HIDDENCheck for all bonded interactions with DD" },
- { "-ddbondcomm", FALSE, etBOOL, {&bDDBondComm},
- "HIDDENUse special bonded atom communication when -rdd > cut-off" },
- { "-rdd", FALSE, etREAL, {&rdd},
- "HIDDENThe maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates" },
- { "-rcon", FALSE, etREAL, {&rconstr},
- "HIDDENMaximum distance for P-LINCS (nm), 0 is estimate" },
- { "-dlb", FALSE, etENUM, {dddlb_opt},
- "HIDDENDynamic load balancing (with DD)" },
- { "-dds", FALSE, etREAL, {&dlb_scale},
- "HIDDENMinimum allowed dlb scaling of the DD cell size" },
- { "-ddcsx", FALSE, etSTR, {&ddcsx},
- "HIDDENThe DD cell sizes in x" },
- { "-ddcsy", FALSE, etSTR, {&ddcsy},
- "HIDDENThe DD cell sizes in y" },
- { "-ddcsz", FALSE, etSTR, {&ddcsz},
- "HIDDENThe DD cell sizes in z" },
- { "-gcom", FALSE, etINT,{&nstglobalcomm},
- "HIDDENGlobal communication frequency" },
- { "-compact", FALSE, etBOOL,{&bCompact},
- "Write a compact log file" },
- { "-seppot", FALSE, etBOOL, {&bSepPot},
- "HIDDENWrite separate V and dVdl terms for each interaction type and node to the log file(s)" },
- { "-pforce", FALSE, etREAL, {&pforce},
- "HIDDENPrint all forces larger than this (kJ/mol nm)" },
- { "-reprod", FALSE, etBOOL,{&bReproducible},
- "HIDDENTry to avoid optimizations that affect binary reproducibility" },
- { "-multi", FALSE, etINT,{&nmultisim},
- "HIDDENDo multiple simulations in parallel" },
- { "-replex", FALSE, etINT, {&repl_ex_nst},
- "HIDDENAttempt replica exchange every # steps" },
- { "-reseed", FALSE, etINT, {&repl_ex_seed},
- "HIDDENSeed for replica exchange, -1 is generate a seed" },
- { "-rerunvsite", FALSE, etBOOL, {&bRerunVSite},
- "HIDDENRecalculate virtual site coordinates with -rerun" },
- { "-ionize", FALSE, etBOOL,{&bIonize},
- "HIDDENDo a simulation including the effect of an X-Ray bombardment on your system" },
- { "-confout", TRUE, etBOOL, {&bConfout},
- "HIDDENWrite the last configuration with -c and force checkpointing at the last step" },
- { "-stepout", FALSE, etINT, {&nstepout},
- "HIDDENFrequency of writing the remaining runtime" },
- { "-resetstep", FALSE, etINT, {&resetstep},
- "HIDDENReset cycle counters after these many time steps" },
- { "-resethway", FALSE, etBOOL, {&bResetCountersHalfWay},
- "HIDDENReset the cycle counters after half the number of steps or halfway -maxh" },
- { "-v", FALSE, etBOOL,{&bVerbose},
- "Be loud and noisy" },
- { "-maxh", FALSE, etREAL, {&max_hours},
- "HIDDENTerminate after 0.99 times this time (hours)" },
- { "-cpt", FALSE, etREAL, {&cpt_period},
- "HIDDENCheckpoint interval (minutes)" },
- { "-append", FALSE, etBOOL, {&bAppendFiles},
- "HIDDENAppend to previous output files when continuing from checkpoint" },
- { "-addpart", FALSE, etBOOL, {&bAddPart},
- "HIDDENAdd the simulation part number to all output files when continuing from checkpoint" },
+ { "-xyinit", FALSE, etREAL, {&xy_fac},
+ "Resize factor for the protein in the xy dimension before starting embedding" },
+ { "-xyend", FALSE, etREAL, {&xy_max},
+ "Final resize factor in the xy dimension" },
+ { "-zinit", FALSE, etREAL, {&z_fac},
+ "Resize factor for the protein in the z dimension before starting embedding" },
+ { "-zend", FALSE, etREAL, {&z_max},
+ "Final resize faction in the z dimension" },
+ { "-nxy", FALSE, etINT, {&it_xy},
+ "Number of iteration for the xy dimension" },
+ { "-nz", FALSE, etINT, {&it_z},
+ "Number of iterations for the z dimension" },
+ { "-rad", FALSE, etREAL, {&probe_rad},
+ "Probe radius to check for overlap between the group to embed and the membrane"},
+ { "-pieces", FALSE, etINT, {&pieces},
+ "Perform piecewise resize. Select parts of the group to insert and resize these with respect to their own geometrical center." },
+ { "-asymmetry",FALSE, etBOOL,{&bALLOW_ASYMMETRY},
+ "Allow asymmetric insertion, i.e. the number of lipids removed from the upper and lower leaflet will not be checked." },
+ { "-ndiff" , FALSE, etINT, {&low_up_rm},
+ "Number of lipids that will additionally be removed from the lower (negative number) or upper (positive number) membrane leaflet." },
+ { "-maxwarn", FALSE, etINT, {&maxwarn},
+ "Maximum number of warning allowed" },
+ { "-stepout", FALSE, etINT, {&nstepout},
+ "HIDDENFrequency of writing the remaining runtime" },
+ { "-v", FALSE, etBOOL,{&bVerbose},
+ "Be loud and noisy" },
+ { "-mdrun_path", FALSE, etSTR, {&mdrun_path},
+ "Path to the mdrun executable compiled with this g_membed version" }
};
- gmx_edsam_t ed;
- unsigned long Flags, PCA_Flags;
- ivec ddxyz;
- int dd_node_order;
- gmx_bool HaveCheckpoint;
- FILE *fplog,*fptest;
- int sim_part,sim_part_fn;
- const char *part_suffix=".part";
- char suffix[STRLEN];
- int rc;
-
-
- cr = init_par(&argc,&argv);
-
- PCA_Flags = (PCA_KEEP_ARGS | PCA_NOEXIT_ON_ARGS | PCA_CAN_SET_DEFFNM
- | (MASTER(cr) ? 0 : PCA_QUIET));
-
-
- /* Comment this in to do fexist calls only on master
- * works not with rerun or tables at the moment
- * also comment out the version of init_forcerec in md.c
- * with NULL instead of opt2fn
- */
- /*
- if (!MASTER(cr))
- {
- PCA_Flags |= PCA_NOT_READ_NODE;
- }
- */
-
- parse_common_args(&argc,argv,PCA_Flags, NFILE,fnm,asize(pa),pa,
- asize(desc),desc,0,NULL, &oenv);
-
- /* we set these early because they might be used in init_multisystem()
- Note that there is the potential for npme>nnodes until the number of
- threads is set later on, if there's thread parallelization. That shouldn't
- lead to problems. */
- dd_node_order = nenum(ddno_opt);
- cr->npmenodes = npme;
-#ifdef GMX_THREADS
- /* now determine the number of threads automatically. The threads are
- only started at mdrunner_threads, though. */
- if (nthreads<1)
- {
- nthreads=tMPI_Get_recommended_nthreads();
- }
-#else
- nthreads=1;
-#endif
+ FILE *data_out;
+ output_env_t oenv;
+ char buf[256],buf2[64];
- if (repl_ex_nst != 0 && nmultisim < 2)
- gmx_fatal(FARGS,"Need at least two replicas for replica exchange (option -multi)");
+ parse_common_args(&argc,argv,0, NFILE,fnm,asize(pa),pa,
+ asize(desc),desc,0,NULL, &oenv);
- if (nmultisim > 1) {
-#ifndef GMX_THREADS
- init_multisystem(cr,nmultisim,NFILE,fnm,TRUE);
-#else
- gmx_fatal(FARGS,"mdrun -multi is not supported with the thread library.Please compile GROMACS with MPI support");
-#endif
- }
+ data_out = ffopen(opt2fn("-dat",NFILE,fnm),"w");
+ fprintf(data_out,"nxy = %d\nnz = %d\nxyinit = %f\nxyend = %f\nzinit = %f\nzend = %f\n"
+ "rad = %f\npieces = %d\nasymmetry = %s\nndiff = %d\nmaxwarn = %d\n",
+ it_xy,it_z,xy_fac,xy_max,z_fac,z_max,probe_rad,pieces,
+ bALLOW_ASYMMETRY ? "yes" : "no",low_up_rm,maxwarn);
+ fclose(data_out);
- /* Check if there is ANY checkpoint file available */
- sim_part = 1;
- sim_part_fn = sim_part;
- if (opt2bSet("-cpi",NFILE,fnm))
- {
- bAppendFiles =
- read_checkpoint_simulation_part(opt2fn_master("-cpi", NFILE,fnm,cr),
- &sim_part_fn,NULL,cr,
- bAppendFiles,NFILE,fnm,
- part_suffix,&bAddPart);
- if (sim_part_fn==0 && MASTER(cr))
- {
- fprintf(stdout,"No previous checkpoint file present, assuming this is a new run.\n");
- }
- else
- {
- sim_part = sim_part_fn + 1;
- }
- }
- else
+ sprintf(buf,"%s -s %s -membed %s -o %s -c %s -e %s -nt 1 -cpt -1",
+ (mdrun_path==NULL) ? "mdrun" : mdrun_path,
+ opt2fn("-f",NFILE,fnm),opt2fn("-dat",NFILE,fnm),opt2fn("-o",NFILE,fnm),
+ opt2fn("-c",NFILE,fnm),opt2fn("-e",NFILE,fnm));
+ if (opt2bSet("-n",NFILE,fnm))
{
- bAppendFiles = FALSE;
- }
-
- if (!bAppendFiles)
- {
- sim_part_fn = sim_part;
- }
-
- if (bAddPart && sim_part_fn > 1)
- {
- /* This is a continuation run, rename trajectory output files
- (except checkpoint files) */
- /* create new part name first (zero-filled) */
- sprintf(suffix,"%s%04d",part_suffix,sim_part_fn);
-
- add_suffix_to_output_names(fnm,NFILE,suffix);
- fprintf(stdout,"Checkpoint file is from part %d, new output files will be suffixed '%s'.\n",sim_part-1,suffix);
- }
-
- Flags = opt2bSet("-rerun",NFILE,fnm) ? MD_RERUN : 0;
- Flags = Flags | (bSepPot ? MD_SEPPOT : 0);
- Flags = Flags | (bIonize ? MD_IONIZE : 0);
- Flags = Flags | (bPartDec ? MD_PARTDEC : 0);
- Flags = Flags | (bDDBondCheck ? MD_DDBONDCHECK : 0);
- Flags = Flags | (bDDBondComm ? MD_DDBONDCOMM : 0);
- Flags = Flags | (bConfout ? MD_CONFOUT : 0);
- Flags = Flags | (bRerunVSite ? MD_RERUN_VSITE : 0);
- Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
- Flags = Flags | (bAppendFiles ? MD_APPENDFILES : 0);
- Flags = Flags | (sim_part>1 ? MD_STARTFROMCPT : 0);
- Flags = Flags | (bResetCountersHalfWay ? MD_RESETCOUNTERSHALFWAY : 0);
-
-
- /* We postpone opening the log file if we are appending, so we can
- first truncate the old log file and append to the correct position
- there instead. */
- if ((MASTER(cr) || bSepPot) && !bAppendFiles)
+ sprintf(buf2," -mn %s",opt2fn("-n",NFILE,fnm));
+ strcat(buf,buf2);
+ }
+ if (opt2bSet("-x",NFILE,fnm))
{
- gmx_log_open(ftp2fn(efLOG,NFILE,fnm),cr,!bSepPot,Flags,&fplog);
- CopyRight(fplog,argv[0]);
- please_cite(fplog,"Hess2008b");
- please_cite(fplog,"Spoel2005a");
- please_cite(fplog,"Lindahl2001a");
- please_cite(fplog,"Berendsen95a");
+ sprintf(buf2," -x %s",opt2fn("-x",NFILE,fnm));
+ strcat(buf,buf2);
}
- else
+ if (opt2bSet("-p",NFILE,fnm))
+ {
+ sprintf(buf2," -mp %s",opt2fn("-p",NFILE,fnm));
+ strcat(buf,buf2);
+ }
+ if (bVerbose)
{
- fplog = NULL;
+ sprintf(buf2," -v -stepout %d",nstepout);
+ strcat(buf,buf2);
}
- ddxyz[XX] = (int)(realddxyz[XX] + 0.5);
- ddxyz[YY] = (int)(realddxyz[YY] + 0.5);
- ddxyz[ZZ] = (int)(realddxyz[ZZ] + 0.5);
-
- /* even if nthreads = 1, we still call this one */
-
- rc = mdrunner_membed(fplog, cr, NFILE, fnm, oenv, bVerbose, bCompact,
- nstglobalcomm,
- ddxyz, dd_node_order, rdd, rconstr, dddlb_opt[0], dlb_scale,
- ddcsx, ddcsy, ddcsz, nstepout, resetstep, nmultisim, repl_ex_nst,
- repl_ex_seed, pforce, cpt_period, max_hours, deviceOptions, Flags,
- xy_fac,xy_max,z_fac,z_max,
- it_xy,it_z,probe_rad,low_up_rm,
- pieces,bALLOW_ASYMMETRY,maxwarn);
+ printf("%s\n",buf);
+ system(buf);
- if (gmx_parallel_env_initialized())
- gmx_finalize();
-
- if (MULTIMASTER(cr)) {
- thanx(stderr);
- }
-
- /* Log file has to be closed in mdrunner if we are appending to it
- (fplog not set here) */
- fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
-
- if (MASTER(cr) && !bAppendFiles)
- {
- gmx_log_close(fplog);
- }
+ fprintf(stderr,"Please cite:\nWolf et al, J Comp Chem 31 (2010) 2169-2174.\n");
- return rc;
+ return 0;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff