Acetone Solvent Box Simulation

User discussions
1

Dear Sir,
Following energy minimization, NVT, and NPT simulations, I visualized the npt.gro file in VMD. Notably, an empty portion resembling a hole is evident within the acetone solvent box. Seeking guidance on addressing this issue.

image
image770×560 155 KB

em file
; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions

NVT
title = OPLS Lysozyme NVT equilibration
define = ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000000 ; 2 * 50000000 = 100000 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
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
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 = 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 off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed

NPT
title = OPLS Lysozyme NPT equilibration
define = ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10000 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
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
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 0.5 ; short-range electrostatic cutoff (in nm)
rvdw = 0.5 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 = 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 = Berendsen
pcoupltype = isotropic
nstpcouple = -1
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility = 2.31943e-09 ; isothermal compressibility of acetone, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off

Best regards,
LUKHMANUL HAKEEM K.

2

Hi,
Nanobubbles should not appear in NPT simulations of liquids. My best guess is that the system hasn’t equilibrated yet. Did you have a look at the box volume over time? The isothermal compressibility is very small, you can increase it (even setting it to 4.5e-5 bar^-1, like the one of water) to speed-up the relaxation of the box volume. You can set it back to its physically accurate value once the system has reached equilibrium. Also, consider using c-rescale instead of berendsen for pressure coupling. Moreover, rvdw is quite small; it’s usually about 0.9 or 1.0 nm. That could also mess up things in systems with interfaces.

3

Sir,
I am uncertain about the proper methods to confirm whether my system has equilibrated effectively. While I did not monitor the box volume over time. I was did the NPT simulation of acetone box by giving the 4.5e-5 bar^-1 as isothermal compressibility value . Then by using this files I was calculated by ‘gmx energy’(i got small value). Then again I did NPT simulation by using this isothermal compressibility value. This are the steps i Followed in my previous simulation. Is it correct or not?.
Could you please explain your message little bit more?.
So my next step is do the simulation of perylene derived dimer in this acetone solvent box. In experimental observations, it has been noted that the residues of this dimer tend to approach each other in polar solvents, while exhibiting a tendency to remain apart in non-polar solvents. To provide empirical confirmation of these findings, I am keen to conduct MD simulations. I seek your expertise on the necessary adjustments to be made in the energy minimization, NVT, NPT, and production run files to achieve accurate results. Your insights and guidance on these matters would be immensely valuable to me.
I eagerly await your response.

Warm regards,
LUKHMANUL HAKEEM K.

4

This is expected mostly for the polar solvents, equilibrate the system for longer time ~10-20 ns.

5

And, as Michele pointed out above, make sure that you use simulation settings that are suitable for the force field you are using. I cannot think of any circumstance where switching rvdw and rcoul from 1.2 nm to 0.5 nm, when going from NVT to NPT is motivated.

6

Ok, Thank you sir for your reply.

7

And I believe that the low compressibility is also a problem. I think it’s risky basing it on simulations, especially if those simulations might not be correct. I would suggest that you use reported experimental values instead. Unless I’m wrong, acetone should have almost three times higher compressibility than water - not more than a factor 1000 lower.

8

Thank You sir for your reply.

9

Dear Sir,
I hope this message finds you well. I wanted to share with you the results of the recent simulation I conducted. After making modifications to the EM, NVT, and NPT files, I ran the simulation and believe that the system has now reached equilibrium. But Upon reviewing the NPT.gro file in VMD, I noticed that certain acetone molecules appear to extend beyond the simulation box boundaries.

Attached are the plots depicting temperature, energy, and density variations throughout the simulation and also em.mdp,nvt.mdp and npt.mdp files,NPT.gro files. Kindly review them at your convenience and let me know if you have any insights or further suggestions.

; minim.mdp - used as input into grompp to generate em.tpr
; Parameters describing what to do, when to stop and what to save
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Minimization step size
nsteps = 5000000 ; Maximum number of (minimization) steps to perform

; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
vdw-type = cutoff
vdw-modifier= force-switch
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.2 ; Short-range electrostatic cut-off
rvdw = 1.2 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions in all 3 dimensions

; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10000 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
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
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 = 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 off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed

; Run parameters
integrator = md ; leap-frog integrator
nsteps = 5000000 ; 2 * 5000000 = 10000 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = yes ; Restarting after NVT
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
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet scheme
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; 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 = system ; two coupling groups - more accurate
tau_t = 2 ; time constant, in ps
ref_t = 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = c-rescale
pcoupltype = isotropic
nstpcouple = -1
tau_p = 4.0 ; time constant, in ps
ref_p = 1.0 ; reference pressure, in bar
compressibility =1.26e-4 bar^-1 ; isothermal compressibility of acetone, bar^-1
refcoord_scaling = com
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = no ; Velocity generation is off
image

energy_plot
energy_plot1200×900 57 KB

temp
temp1200×900 59 KB

dens
dens1200×900 87.8 KB

I eagerly await your response.

Warm regards,
LUKHMANUL HAKEEM K.

10

That molecules cross the periodic boundary (or if they get broken by it) is not a problem: https://manual.gromacs.org/2024.1/user-guide/terminology.html#periodic-boundary-conditions.

From your plot your system definitely looks equilibrated and the density looks correct as well. It’s difficult to say exactly from the plot, but it looks reasonably close to the experimental value ~780 kg/m3 at 300K. This is an indication that the force field parameters are good.

11

Dear Sir
Thank you for your reply. So my next step is to solvate perylene dimer using this acetone solvent box, followed by simulating and visualizing the geometry of the perylene molecule. Our experimental observations indicate close proximity of the perylene molecules in polar solvents, and we intend to validate this through MD simulations.
However, a concern has arisen regarding the utilization of two different force fields: one for the perylene dimer, derived from the AMBER Antechamber tool and converted to GROMACS format using ACPYPE, and another for acetone, created through CHARMM-GUI. Could this discrepancy in force fields pose a problem for our simulations? Additionally, could you advise on any necessary adjustments to the energy minimization (EM), NVT ensemble, and NPT ensemble parameters to ensure results align more closely with experimental observations?. Please help me.

I eagerly await your response.

Warm regards,
LUKHMANUL HAKEEM K.

12

I would recommend against mixing force fields. They recommend different water models (which might not be a problem in your case) and different nonbonded interaction parameters, and potentially also different dispersion correction settings. I’d suggest using either Amber or Charmm based force fields.

I think your simulations parameters look OK at a glance (but I won’t take any responsibility if you’re publishing the data - you have to go through your settings yourself). You might want to run at 298K instead of 300K and you could probably use a tau_p of 10 ps. I would also recommend writing less often, unless you really know you need all that data:

nstxout-compressed = 50000 ; save compressed coordinates every 100 ps
nstvout = 0 ; Do not save velocities
nstenergy = 100000 ; save energies every 200 ps
nstlog = 100000 ; update log file every 200 ps

Should be enough for most purposes (not necessarily all). Especially if you intend to store that simulation results for a long time, it’s often a waste of disk space if you write too often. I think 100-500 compressed coordinate frames per simulation is usually enough.

13

Ok, Thank you sir for your reply.