[gmx-users] Wall potential for a membrane-simulation

[Marianne Schulte]

Hi everyone!







I'm new to Gromacs and trying to
simulate a membrane system with two walls, one at the bottom of my
box at z=0 and one at the top, using the gromos53a6 forcefield
(GROMACS version 4.5.5).


My testing system consists of a
membrane in the middle, water and sodium ions (40) above the membrane
and water and chloride ions (also 40) beneath it.


First, I built up the system, did an
energy minimization and then implemented the walls for the
nvt-equilibration. I didn't change anything in the topology.  (Do I
have to define the wall-atem-type there?  If yes, what exactly do I
have to put in the topology?) 



My nvt.mdp file looks like that:







title		= NVT equilibration
for KALP15-DPPC 



define		= -DSTRONG_POSRES
-DPOSRES_LIPID   	; position restrain the lipid


; Run parameters


integrator	= md                                                        		;
leap-frog integrator


nsteps		= 50000                                                         		; 2 *
50000 =  100 ps


dt		= 0.002                                                                  ; 2 fs  



; Output control


nstxout		= 100                                                            ; save
coordinates every 0.2 ps


nstvout		= 100                                                            ; save
velocities every 0.2 ps


nstenergy	= 100                                                        ; save
energies every 0.2 ps


nstlog		= 100                                                              ; update
log file every 0.2 ps


; Bond parameters


continuation	= yes                                                    ;
first dynamics run


constraint_algorithm =
lincs                                     ; holonomic constraints 



constraints	= all-bonds                                             ;
all bonds (even heavy atom-H bonds) constrained


lincs_iter	= 1                                                             ; accuracy
of LINCS


lincs_order	= 4                                                          ; also
related to accuracy


; Neighborsearching


ns_type		= grid                                                          		; search
neighboring grid cels


nstlist		= 5	                                                                	; 10 fs


rlist		= 1.2                                                                 		;
short-range neighborlist cutoff (in nm)


rcoulomb	= 1.2                                                        		;
short-range electrostatic cutoff (in nm)


rvdw		= 1.2                                                                		; short-range
van der Waals cutoff (in nm)


; Electrostatics


coulombtype	= PME                                               		;
Particle Mesh Ewald for long-range electrostatics


pme_order	= 4                                                        		; cubic
interpolation


fourierspacing	= 0.16                                             		;
grid spacing for FFT


; Temperature coupling is
on


tcoupl		= V-rescale                                                  	;
modified Berendsen thermostat


tc-grps		= lipid
SOL_NA_CL                                    ; two coupling groups - more accurate


tau_t		= 0.1	0.1                                                            	;
time constant, in ps


ref_t		= 303 	303                                                          	;
reference temperature, one for each group, in K


; Pressure coupling is off


pcoupl		= no                                                             		; no
pressure coupling in NVT


; Periodic boundary
conditions


;wall


nwall = 2


wall_type = 10-4


wall_density = 100
100		


wall_atomtype = C C


wall_r_linpot = 1 1	                                    	;
with -1 I got a fatal error, because atoms were beyond the wall


wall_ewald_zfac = 3


ewald_geometry= 3dc


pbc = xy                                                   			; Periodic
Boundary Conditions in x/y direction


; Dispersion correction


DispCorr	= EnerPres                                             	;
account for cut-off vdW scheme


; Velocity generation


gen_vel		= yes	                                                        	; assign
velocities from Maxwell distribution


gen_temp	= 303                                                     		;
temperature for Maxwell distribution


gen_seed	= -1                                                        		; generate
a random seed


; COM motion removal


; These options remove
motion of the protein/bilayer relative to the solvent/ions


nstcomm	= 10


comm-mode	= Linear


comm-grps	= lipid
SOL_NA_CL







What I see in VMD is the
following: as expected some NA-ions are going through the membrane to
reach the chloride/water side. But at the same time water molecules
(and one chloride-ion) at both sides disappear in z-direction and
never come back. In the first 10 ps (with a tiny time-step of 0.1 fs)
there were only a few molecules that went outside the box, but after
10 ps it increased dramatically and at some point the equilibration
crashed.


I already played around
with time-steps, atom_types and density, but couldn't find anything
that works to keep the molecules inside the box..







Anyone has an idea what I
did wrong and how I can correct it?  Already read the manual and
mailing-list archive, but could not find anything that really helped
me....












Thank you so much in
advance!!!!


Cheers, 



Marianne