[gmx-users] walls combining with the PME and pressure coupling

[彭春望]

Dear all,


I have sent this e-mail to the gmx-users mailing list for a few days, but without any reply until now. I'm not sure if it has something to do with the format, because the link can’t be opened when I search the mailing list. So this time I use the plain text and send it once more, would anyone please help me with this?


I’m working on some simulations about the adsorption of protein on solid surfaces, which have slab geometry in the x-y plane. In order to reduce the amount of water molecules and at the same time to decrease the unphysical Coulomb interaction between periodic images in the z-direction, an empty layer should be added in the z-direction. To prevent the water molecules from evaporating into the vacuum layer, I’m going to use the wall option. However, there are some details that I’m not pretty sure about.


What are your suggestions about the choice of wall_atomtype and wall_type? Is it necessary to leave some room for the two walls (or it may lead to high interaction energies between the walls and the system) and how to determine its height (in my simulations, I leave about 1.5 Angstrom, respectively)? Does it have something to do with the wall_atomtype and wall_type?


When choosing nwall=2, pressure coupling and Ewald summation can be used (it is usually best to use semi-isotropic pressure coupling with x/y compressibility set to 0). It means the wall can move in the z-direction? Can the pressure coupling be used in system with fixed atoms and how can we control/calculate the pressure of the mobile phase (as it has been discussed in the paper [Biointerphases 5, 85 (2010)] using the CHARMM package)? 


When combining walls with the PME method, it is suggested the eward_geometry be set to 3dc and the wall_eward_zfac be 3. Does this mean there will be an empty layer whose height is 3 times the slab height added to increase the z-dimension of the box? And I’m not sure about the exact meaning of the “slab height”; it seems to be the length/width of the slab (as described in the paper [J. Chem. Phys. 111, 3155 (1999)]).


I’m sorry for troubling you with so many questions. But I really need your help. Anyone could help me? Thank you!!


Below is an attachment of mdp file used in my simulation work. If there is anything wrong, please be kind to point it out. Thanks again!


md.mdp
------------------------------------------------------------------------------
title= cyt-c on Au MD 
; Run parameters
integrator= md; leap-frog integrator
nsteps= 10000000; 2 * 10000000 = 20000 ps, 20 ns
dt= 0.002; 2 fs
comm-mode     = Linear       ; mode for center of mass motion removal


; Output control
nstxout= 10000; save coordinates every 20 ps
nstvout= 10000; save velocities every 20 ps
nstenergy= 10000; save energies every 20 ps
nstlog= 10000; update log file every 20 ps


; Selection of energy groups
energygrps    = Protein_HEM  GLD ;Group(s) to write to energy file


; Bond parameters
continuation= yes; Restarting after NPT 
constraint_algorithm = lincs; holonomic constraints, GolP has been tested with lincs only 
constraints= hbonds; bonds with H-atoms constrained
lincs_iter= 1; accuracy of LINCS
lincs_order= 4; also related to accuracy


; Neighborsearching
ns_type= grid; search neighboring grid cells
nstlist= 10; 20 fs
rlist= 1.1; short-range neighborlist cutoff (in nm)
; Periodic boundary conditions
pbc= xy; 2-D PBC


; Method for doing Van der Waals
vdw-type      = switch
rvdw-switch   = 0.9
rvdw= 1.0; short-range van der Waals cutoff (in nm)


; Electrostatics
coulombtype= PME; Particle Mesh Ewald for long-range electrostatics
rcoulomb= 1.1; short-range electrostatic cutoff (in nm)
pme_order= 4; cubic interpolation
fourierspacing= 0.12; grid spacing for FFT
ewald_rtol     = 1e-5
ewald_geometry  = 3dc


; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx               = 0
fourier_ny               = 0
fourier_nz               = 0
optimize_fft             = yes


; Temperature coupling is on
tcoupl= Nose-Hoover; Nose-Hoover thermostat
tc-grps= Protein_HEM GLDWater_and_ions  ; three coupling groups - more accurate
tau_t= 0.50.50.5; time constant, in ps
ref_t= 300 300 300; reference temperature, one for each group, in K


; Pressure coupling is on
pcoupl= Parrinello-Rahman; Pressure coupling on in NPT
pcoupltype= semiisotropic; nonuniform scaling of box vectors
tau_p= 1.01.0; time constant, in ps
ref_p= 1.01.0; reference pressure, in bar
compressibility = 04.5e-5; isothermal compressibility of water, bar^-1


; Velocity generation
gen_vel= no; Velocity generation is off 


; Non-equilibrium MD stuff
freezegrps    = LOCK
freezedim     = Y  Y  Y


; WALLS                
; Number of walls, type, atom types, densities and box-z scale factor for Ewald
nwall                    = 2
wall_type                = 9-3
wall_r_linpot            = -1  -1
wall_atomtype            = OWT3  OWT3     ; oxygen atom of TIP3P water in charmm27.ff
wall_density             = 33.4  33.4
wall_ewald_zfac          = 3
---------------------------------------------------------------------------


Yours, sincerely!
Chunwang Peng


Chemistry & Chemical Engineering，
South China University of Technology，
Tianhe District, Guangzhou, China.