GROMACS version:
GROMACS modification: Yes/No
Here post your question
Dear users,
I am applying an electric field of 1.5 V/nm along the Z-axis to a system with a liquid-liquid interface. The system has the interface aligned perpendicular to the Z-axis, but after the MD simulation, the interface becomes aligned parallel to the Z-axis. md.mdp is :
title = OPLS Lysozyme NPT equilibration
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 100000000 ; 2 * 500000 = 1000 ps (1 ns)
dt = 0.001 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 10000 ; save energies every 10.0 ps
nstlog = 10000 ; update log file every 10.0 ps
nstxout-compressed = 10000 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = no ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; 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
electric-field-z = 1.5 0 0 0
; Temperature coupling is on
tcoupl = nose-hoover ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; Apply pressure coupling separately in X-Y and Z
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is off
With a liquid system, the electric field in enhanced by a factor of the dieletric permittivity when using PME with epsilon-surface=infinity. You need the scale down the applied field by the dielectric permittivity. Alternatively you can set epsilon-surface=1, but that will suppress the fluctuation of the dipole of the system by the same factor.
Still you might see get the conformation you see. If the gain in electrostatic energy is higher than the cost of the extra interface area, that is the lowest free-energy state.
PS: Your have a too coarse PME grid. You should use fourier-spacing=0.125 with a cut-off 1 and pme-order=4.
Dear @hess
Thank you for your response. I will follow your suggestion.
I have one more concern: do we need to add any electrodes to the system manually, or will the simulation automatically assign the positive and negative electrodes? In several papers, I have seen references to the movement of positive and negative electrodes, but the figures do not explicitly show any added electrodes.
Thank you again for your previous guidance regarding distortion upon applying an electric field.
As a follow-up, I am now studying the effect of an electric field on polymer chains at a liquid–liquid (oil–water) interface. The system consists of multiple polymer chains in the organic phase and counterions in the aqueous phase, with graphene electrodes at the top and bottom of the box (along the z-axis) (inorder to hinder the transfer of polymer cation from one organic phase to the other I used graphene sheets). I applied a field of 1.5 V/nm in the z-direction.
The simulations done in two ways:
NVT → NPT → MD, all with the electric field.
NVT → NPT without EF, then MD with EF.
In both cases, I observed no significant change in polymer conformation or alignment due to the field. However, some polymer chains, which were initially randomly distributed in the organic phase before simulation, are now partially adsorbed near the graphene electrodes after applying the field.
Could this indicate an issue with the force field, system size, or field strength? I’d appreciate any suggestions on how to improve the setup to better observe field-induced effects.
This is the mdp file I used:
title = OPLS Lysozyme NPT equilibration
; Run parameters
integrator = md ; leap-frog integra
define =-DPOSRES
nsteps = 200000000 ; 2 * 500000 = 1000 ps (1 ns)
dt = 0.001 ; 2 fs
; Output control
nstxout = 0 ; suppress bulky .trr file by specifying
nstvout = 0 ; 0 for output frequency of nstxout,
nstfout = 0 ; nstvout, and nstfout
nstenergy = 10000 ; save energies every 10.0 ps
nstlog = 10000 ; update log file every 10.0 ps
nstxout-compressed = 10000 ; save compressed coordinates every 10.0 ps
compressed-x-grps = System ; save the whole system
; Bond parameters
continuation = no ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
comm-mode = linear
comm-grps = System
; Neighborsearching
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
ewald-geometry = 3dc
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
electric-field-z = 1.5 0 0 0
; Temperature coupling is on
tcoupl = nose-hoover ; modified Berendsen thermostat
tc-grps = system ; two coupling groups - more accurate
tau_t = 0.1 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT
pcoupltype = isotropic ; Apply pressure coupling separately in X-Y and Z
tau_p = 5.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; Velocity generation is off
periodic_molecules = yes
freezegrps = RNA
freezedim = Y Y Y
Output configuration is attached below:
As you suggested, instead of using the electric-field-z option, I would like to try creating an electric field by setting charges on two surfaces in my system.
Could you please advise on how to correctly assign and maintain the surface charges in GROMACS? Specifically:
Should I modify the atomic charges of the graphene sheet atoms directly in the topology?
