[gmx-users] Naughty Vacuum Bubble in our Vesicle!

[Björn Sommer]

Dear all,

we are trying to simulate a vesicle in water using united-atoms 
(Gromos96/ffG45a3). The system was modelled with the VesicleBuilder and 
the MembraneEditor. So first the vesicle was built (with 3 components: 2 
PC, 1 Chol), and then it was embedded in a water (spc216) box with 
genbox. The membrane-intersecting water was removed by a custom Python 
script in VMD. After the removel or the intersecting water, the water 
seems to be very well enclosed in the inner membrane, without 
intersecting water atoms and with only a little space between the inner 
head groups and the water.

The system minimization in water (spc216) worked pretty well, but after 
NPT-equillibration I found a vacuum bubble in the intracellular room of 
the vesicle.

I tried to do a NVT-equillibration before the NPT, which ended with the 
same result.

Using different barostats (parrinello-rahman, Berendsen) and 
refcoord-scaling options could'nt change anything, too.

To "repair" the vacuum I tried to manually insert some water using 
pymol, after equillibrating this system again, i got the same result 
with an even bigger bubble in the centre.

Analysing my system with g_energy showed a volume increase about 400nm³, 
which is rawly 4% more than the volume of the starting system. The 
systems energy increased by 20 kJ/mol.

Trying to simulate the vesicle with the vacuum bubble inside resulted in 
a deformed vesicle and an increased distance between the outer and the 
inner lipid layer.

The NPT.mdp file I used is the following:

;**********************************************************
; NEIGHBORSEARCHING PARAMETERS =
; nblist update frequency =
nstlist                  = 5
; ns algorithm (simple or grid) =
ns_type                  = grid
; Periodic boundary conditions: xyz or none =
pbc                      = xyz
; nblist cut-off         =
rlist                    = 1.6

; OPTIONS FOR ELECTROSTATICS AND VDW =
; Method for doing electrostatics =
coulombtype              = PME
rcoulomb_switch          = 0.0
rcoulomb                 = 1.6
; Method for doing Van der Waals =
vdw_type                 = Shift
; cut-off lengths        =
rvdw_switch              = 0.9
rvdw                     = 1.0
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr                 = AllEnerPres

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; Temperature coupling   =
tcoupl                   = Berendsen
; Groups to couple separately =
;TODO: für mehrere Lipidtypen anpassen
tc-grps                     = CHO DPC DPE SOL
; Time constant (ps) and reference temperature (K) =
;TAUT
tau_t                    = 0.1 0.1 0.1 0.1
;REFT
ref_t                    = 300 300 300 300
; Pressure coupling      =
Pcoupl                   = berendsen
Pcoupltype               = isotropic ;semiisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar) =
tau_p                    = 4.0  4.0
compressibility          = 3e-5 3e-5
ref_p                    = 1.0  1.0
refcoord-scaling         = no

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel                  = no
gen_temp                 = 105
gen_seed                 = 473529

; OPTIONS FOR BONDS     =
constraints              = all-bonds
fourierspacing       =
pme_order            =  6
optimize_fft         =  yes
; Type of constraint algorithm =
constraint_algorithm     = Lincs
; Do not constrain the start configuration =
unconstrained_start      = no
; Highest order in the expansion of the constraint coupling matrix =
lincs_order              = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle          = 30

;**********************************************************

Also, I tried to relax my vesicle under NPT in vacuum in order to add 
water in a later step, but the NPT ended with a lot of LINCS warning 
(rotation more than 30 degrees). Okay, we read already that there are 
some information in the gmx-list discussing this problem, but the 
question is if it basically makes sense to follow the idea of a vaccum 
simulation or if we should directly start with solvated system in any case.

Can you suggest us any method to solve this problem or maybe help us to 
improve our .mdp? Would be great!

Best wishes,
Manuel (and Björn)