I want to calculate the thermal conductivity for my system. Thus, I am trying to run NVE for my system after a 20 ns equilibration using the NVT ensemble at 300 K. I used the Nosé-Hoover thermostat. My box is orthogonal with dimensions x = 57 Å, y = 57 Å, z = 100 Å, and gamma = 120 degrees. Inside the box, there is a solid slab is at the bottom (0–15 Å), above of that there is liquid (15–60 Å), and vacuum from (60–100) Å. When I run NVE, I see a continuous energy drift (straight line type increase).
Initially my timestep was 1 fs, I decreased it to 0.5 fs, but the same error continues.
I also tried to tighten the kspace from 1.0e^-5 to 1.0e^-6 , but again the same issue continues.
I am not able to figure out how to fix it.
For your reference, I am sharing my input file.
units real
dimension 3
boundary p p p
atom_style full
neigh_modify every 1 delay 0 check yes
neighbor 1.0 bin
pair_style lj/class2/coul/long 12
pair_modify shift yes mix arithmetic
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
special_bonds lj/coul 0.0 0.0 1 angle yes dihedral yes # if 1-3 atoms connected by angle, lj yes
read_restart equilibration.rst
kspace_style pppm 1.0e-5
variable dt equal 0.5
timestep ${dt}
thermo 1000
log nve.log
thermo_style custom step temp press density etotal ke pe ebond eangle edihed eimp evdwl ecoul elong etail enthalpy
fix 1 all nve
fix recenter_slab slab recenter INIT INIT INIT
restart 10000 nve1_1.rst nve1_2.rst
dump mydcd all dcd 10000 nve.dcd
dump_modify mydcd unwrap yes
run 200000000 #production run for 100ns
write_data nve.data
For a sample with slab geometry and long-range electrostatics, you should use:
boundary p p f
kspace_modify slab 3.0
To decouple the electrostatic interactions perpendicular to the slab. Also, you should use appropriate wall fixes if there are atoms in the gas phase. I would also increase the skin distance to 2.5, to be on the safe side.
@Ankit from your description of the system, there is a solid slab is at the bottom (0–15 Å), above of that there is liquid (15–60 Å), and vacuum from (60–100) Å, I suppose there is a vacuum-liquid interface, which is stabilized by the thermostat in the preceding simulation. When you switch to nve, there is no such stabilizer anymore. The momenta of the atoms in the restart file, consistent with the canonical ensemble distribution at 300K allow the atoms in the liquid phase to escape to the vacuum. You may see this effect by visualizing the system. The increase in the energy is likely due to the interaction between the atoms raising up.
Yes, by commenting out fix recenter command shows little bit lesser upward drift in Energy compared to the one when it is not commented out. for e.g., Along with the use of Recenter command energy drift is around 5000 for a 5ns run, but the same drift would be around 3000 if i commented out fix recenter command.
But in the preceding I have seen the slab seems drifting away from its original position when no fix rigid command is implemented.
Yes for your reference i am sharing the restart file. Remember i am using 2024 version of lammps. This is the restart file after NVT simulation of 20ns.
@hothello Thanks for your suggestion. But as I said my box is non orthogonal and kspace modify slab 3.0 is not supported for non-orthogonal box.
Increasing the skin distance I could try.
But I suspect the real issue lies with the solid framework/slab , when i ran the NVE simulation for the framework its not conserving energy. Is there anything I could do for that may solve the issue
@ndtrung yes you are correct about the fact that liquid vapour interface and also the solid liquid interface was stabilized by thermostat in the preceding simulation.
When I ran the NVE simulation for liquid only by putting vacuum on both side of liquid, there was also drift in energy which i fixed using fix momentum command of lammps.
But the same thing does not solve the issue for my combined system.
LIke i said I suspect the main issue lies with the slab framework/slab , when I ran NVE simulation for solid framework alone it is also showing drift in energy.
Basically I removed the liquid part of the system (as the poster said it can be reproduced with the solid part only), decreased the time step to 0.01fs (since sometimes this kind of problem is caused by large time steps), then remove the interactions one by one until getting good energy conservation.
The problem is in the dihedral forces. With the dihedral forces enabled I observed abrupt energy changes like this (timestep 0.01fs):
When I replace the dihedral style with dihedral_style zero nocoeff the problem is gone.
I also tried disabling dihedral only while keeping every other force (lj, electrostatic, bond, etc.) untouched; the energy drift is <0.5kcal/mol after 100ps with a 1fs timestep (only the slab; with the liquid a smaller timestep is likely needed due to the hydrogen atoms).
Also as a side note, you can convert a hexagonal lattice to an orthogonal lattice with a lattice transformation, which may speed up the simulation (WARNING: Triclinic box skew is large. LAMMPS will run inefficiently. (src/domain.cpp:221)).
@initialize 1.I tried this part dihedral_style zero nocoeff for only the slab and apply nve integration for that .Although there is no increase in energy for first 100ps but during long run the energy drift slowly.
2. Also I tried to implement the same command for slab+liquid, the energy drifts significantly even for first 10ps.
3. Although when i try for 0.01fs I see better energy conservation,because I need at least 40ns run(since my goal is to calculate thermal conductivity using green kubo) in NVE simulation. Doing that at this small time step is computationally expensive.
Although there is no increase in energy for first 100ps but during long run the energy drift slowly.
That will be your homework I think. You can’t expect others to run lengthy simulations for you.
Also I tried to implement the same command for slab+liquid, the energy drifts significantly even for first 10ps
I would repeat what I said: the energy drift is … 1fs timestep (only the slab; with the liquid a smaller timestep is likely needed due to the hydrogen atoms).
As @initialize pointed out, a 1 x 2 x 1 hexagonal supercell can be readily transformed into an orthogonal one [here a guide].
My guess for the energy drift is that the slab was created from a unit cell that was not relaxed (at NPT conditions) with the force field that you are currently using.
Please note that strictly speaking you cannot expect your simulation to represent a liquid-vapor interface with complete fidelity.
At 300 K, the vapour pressure of water is about 3.5 kPa, which equals an ideal gas volume of about 8,000 L/mol. At the same temperature liquid water occupies about 18 mL/mol (since 1 kg of water is about 55 moles), which is about two-millionths of its molar gas volume.
Thus, if you attempt to simulate a liquid-gas interface of water at 300 K, there are only two possibilities. The first is that you have a simulation which is one-millionth liquid water by volume and have the rest of it be occupied by water vapour. Clearly that is technically unfeasible, but it is theoretically possible. Otherwise, the gaseous portion of your system will be so small that it is no longer anything remotely resembling a thermodynamics limit state, and its small system fluctuations will noticeably deviate from thermodynamic ideal values.
Of course, water vapour is close enough to an ideal gas under those conditions that we don’t care, and your results will be practically usable. But you can’t presume to replicate every little detail of a microcanonical ensemble at that point. And in particular, energy conservation is very difficult in a molecular simulation because when a timestep is small enough to rule out integrator error, you need to have many more timesteps instead – and the more timesteps you have the more floating point error you accumulate overall.
You will simply have to settle for showing that energy is conserved well enough – for example that it doesn’t drift that much over the course of one correlation time period.