High negative pressure during simulation

LAMMPS LAMMPS Beginners
, ,
1

Respected peers,
I am trying to simulate the lignin molecule with a mixture of ozone and water with reaxff potential (I used the reaxff potential available on the LAMMPS example.). I am getting some questionable values for pressure. While running for NVE simulation on the following system: lignin + 500 ozone + 3500 water for 1 ns and 0K starting temperature, the system recorded a very high negative pressure of negative tens of thousand. Later this value saturated to around -2000 as the simulation reaches the end of 1 ns (the temperature at the end of 1 ns is around 473 K). According to LAMMPS documentation, the pressure is supposed to be atm units (I used ‘real’ units). Considering the negative pressure, I incremented water molecules to 5000 but the pressure trend remained the same (I simulated this particular system for 10 ns). I am pretty new with LAMMPS and I will be grateful if someone will help me understand these results
Below is my input file (for lignin + 500 ozone + 3500 water system).

units real
boundary p p p

atom_style charge
read_data data.lignin+ozone+water.charge

pair_style reaxff lmp_control
pair_coeff * * ffield.reax.cho C H O

group lig id 1:201
group ozn id 202:1701
group water id >= 1702
group lig+ozn id 1:1701

compute reax all pair reaxff

#compute pressure all pressure thermo_temp

variable eb equal c_reax[1]
variable ea equal c_reax[2]
variable elp equal c_reax[3]
variable emol equal c_reax[4]
variable ev equal c_reax[5]
variable epen equal c_reax[6]
variable ecoa equal c_reax[7]
variable ehb equal c_reax[8]
variable et equal c_reax[9]
variable eco equal c_reax[10]
variable ew equal c_reax[11]
variable ep equal c_reax[12]
variable efi equal c_reax[13]
variable eqeq equal c_reax[14]

thermo 1

neighbor 2.5 bin
neigh_modify every 10 delay 0 check yes

#velocity all create 300.0 4928459 mom yes rot yes dist gaussian

#fix 1 all nve temp 300.0 300.0 0.1
fix 1 all nve
fix 2 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff
fix 4 all reaxff/bonds 20 bonds.reaxff.lignin+ozone+water_nve300_1ns
fix 5 lig reaxff/bonds 20 bonds.reaxff.lignin+ozone+water(lignin)_nve300_1ns
fix 6 ozn reaxff/bonds 20 bonds.reaxff.lignin+ozone+water(ozone)_nve300_1ns
fix 7 lig+ozn reaxff/bonds 20 bonds.reaxff.lignin+ozone+water(lignin+ozone)_nve300_1ns
#fix 8 all press/berendsen iso 0.0 0.0 1000.0
#fix 5 all bond/break 5 1 1.9

minimize 1.0e-6 1.0e-9 10000 100000

log log.lignin+ozone+water.nve300_1ns

variable nqeq equal f_2

thermo_style custom step temp epair etotal press &
v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa &
v_ehb v_et v_eco v_ew v_ep v_efi v_eqeq v_nqeq

timestep 0.1

dump 1 all atom 20 dump.reaxff.lignin+ozone+water_nve300_1ns
dump 2 lig atom 20 dump.reaxff.lignin+ozone+water(lig)_nve300_1ns
dump 3 ozn atom 20 dump.reaxff.lignin+ozone+water(ozn)_nve300_1ns
dump 4 lig+ozn atom 20 dump.reaxff.lignin+ozone+water(lig+ozn)_nve300_1ns

restart 10000 restart.lignin+ozone+water_nve300_1ns

run 10000

Please let me know if the above-mentioned information is enough to pinpoint the issue. My apologies as being a new user, I couldn’t upload my log file.
Thank you in advance.

2

a) you are using a combustion parameter set for a molecular system. That is a very bad idea.
b) a negative pressure is an indication that your system “wants” to shrink. Does your system have a proper density? You may need to shrink the box. But also you need to equilibrate your system. Also, condensed systems are not very compressible; so small changes in density will have large pressure changes.

3

If you are also new to MD, then it is not a good idea to start with such a complex system and an advanced force field as ReaxFF. You are likely to make many simple “procedural” mistakes that you can avoid by practicing your MD skils with reproducing results from simpler systems with simpler force fields.

4

Hi,
Thank you for your reply.

  1. Can you help me locate that ‘combustion parameter’?
  2. I took a lignin + 500 ozone molecule system and solvated using gromacs to produce an input pdb file. This added around 3500 water molecules to the box. On finding the negative pressure, I manually solvated the system with 5000 water molecules instead, hoping an increased density will make the pressure positive. But the pressure trend remained the same.
    From equilibrium, do you mean energy minimization? Because I did that before NVE.

I have some prior practice of MD with GROMACS but I wanted to perform a ReaxFF simulation so I decided to go with LAMMPS.

5

Please look up the publication describing the parameter set you are using (the reference should be in the file). It will tell you that this is not a parameter set for your kind of system. ReaxFF parameter set are not very transferable. You need to look for some publication describing a parameter set that is suitable for your sytstem.

You are not answering my question and are using a very questionable procedure. You don’t just change the system until you have the desired pressure, but you check whether temperature, pressure, and density conform with corresponding experimental data (or suitable estimates). If then, while reproducing the desired conditions, you cannot reproduce all at the same time, you have an indication that your force field is not able to reproduce the expected state and thus is not suitable. That assumes, of course, that there are no mistakes in your input.

This whole procedure is independent from the MD simulation software and thus definitely not a LAMMPS issue and therefore off-topic for this category. Rather it should be a topic for discussion with your adviser/tutor or experienced colleagues. It is the job of your adviser to arrange for you being properly trained for your work.