I am currently trying to set up a simulation of an ion channel with an applied electric field of 60 mV for direction z. I now encounter the problem that the simulation box keeps shrinking in the z-axis (see snaphot), resulting in a very flat and distorted system. However, GMX does not output any error messages or warnings during the run. Our group also encountered this problem with Gromacs 2018 on a different ion channel. Runs without electric field are stable and do not shrink.
Can you recommend any way to handle this problem? I read in the forum that a solution might be to switch to isotropic pressure coupling, but since the system contains a membrane, I am not sure about the validity of this approach. Might a longer equilibration help?
Thanks in advance for your advice!
All the best,
Theres
Here a bit of information about the system:
equilibration: 1 ns nvt, 1 ns npt
electric field: 60 mV in z direction
mdrun: 50 ns, 4 runs - in all runs the box shrinks in z-dimension during the simulation
mdp file:
; Run parameters
integrator = md
nsteps = 25000000
dt = 0.002
; Output control
nstxout-compressed = 10000
nstenergy = 10000
nstlog = 10000
; Bond parameters
continuation = yes
constraint_algorithm = lincs
constraints = h-bonds
lincs_iter = 1
lincs_order = 4
; Neighborsearching
cutoff-scheme = Verlet
ns_type = grid
nstlist = 10
rlist = 1.0
rcoulomb = 1.0
rvdw = 1.0
; Electrostatics
coulombtype = PME
pme_order = 4
fourierspacing = 0.16
; Temperature coupling is on
tcoupl = V-rescale
tc-grps = Protein_MOL_POPC Water_Ion_ETR
tau_t = 0.5 0.5
ref_t = 310 310
; Pressure coupling is on
pcoupl = Parrinello-Rahman
pcoupltype = semiisotropic
tau_p = 2.0
ref_p = 1.0 1.0
compressibility = 4.5e-5 4.5e-5
; Periodic boundary conditions
pbc = xyz
; Dispersion correction
DispCorr = EnerPres
; Velocity generation
gen_vel = no
; COM motion removal
comm-mode = Linear
comm-grps = Protein_MOL_POPC Water_Ion_ETR
;Electric fieldselectric field
Electric-field-z = 0.06 0 0 0 ; electric field 0.06 Volt per nm in z direction
This appears to be a rather obvious bug. Please try isotropic pressure coupling to see if you at least get sensible coordinates. If you do, this points to a problem specifically in semiisotropic coupling, which is important information for correcting the problem.
As already mentioned, the system is perfectly stable in the free MD simulation. Can you recommend any next steps (perhaps check the lipid density of the membrane or equilibrate longer?) Or should I report my problem
By the way, I forgot to mention that I am using amber99sb in combination with Berger lipids.
Hi, in my case, it helped to do the npt with C-rescale barostat (instead of Parrinello-Rahman). You can also try to increase tau_p, sometimes this helps as well…
I am working for few months in setup a electric field in my simulation unsuccessfully. My system consist o two membranes dividing the space in two chambers with different concentration of ions, the membrane protein is embedded in a membrane with a lipid composition resembling a mammalian membrane. The system without an electric field is stable and have a microsec production run. However, when try to do the same with an electric field the result after equilibration is similar to the theresfriesacher described in his last pic after equilibration.
I am using 2022.4 version in a NPT simulation with a Z compressibility of zero . I tried different barostats
and field intensities without relevant differences.
I noticed that the average pressure is far from 1 bar and remain far as soon that switch on the field.
My question is: Can the barostat handle or take in account the electrostatic pressure?
If my calculation using the virial theorem is correct (a big if). This pressure should be proportional to the sum over all atom charge times its coordinate in the direction of the field ( Sum of qi.zi) that only will be zero or small with a charge distribution symmetrical or almost symmetrical with respect to an arbitrary plane perpendicular to the field.
I don’t understand what is happening. If you compressibility along Z is zero then the system can only scale along x and y, right? If your electric field is along Z then I don’t see why the system would expand much along X and Y. The issue mentioned in the beginning of this thread is typically that the system contracts along the direction of the field.
Right, I was referring to the second. pic (sorry). The system scrambles, the contraction on the direction of the field (Z) is avoided by zeroing compresibility (as you suggested in other thread). In the code, there is a mention of the non-conservativiness of the field (However, in the case of a constant electric field it is conservative) and other about the selfconsistence of the electric field that I do not understand (does it need to be selfconsistent? It is an external force). The “misbehavior” of the pressure when the E is on, make me question if electrostatic pressure is handled well by the barostat in the calculations.
Thanks for your prompt answer.
An electric field can put enormous force on the atoms. These will be compensates by inter and intramolecular forces. These, in turn, contribute to the pressure. Thus the pressure can be very high.
