[lammps-users] NPT cannot hold the pressure

LAMMPS LAMMPS Mailing List Mirror
1

Hi all,

Hope you are doing well.

The NPT ensemble doesn’t work well for my system. I’m trying to run a simulation for a pre-compressed(say, 10% strain) defective domain using NPT with fixed normal pressures(from the simulation for a bulk domain at the same 10% strain state) and fixed temperature. However the normal pressures start to decrease at the beginning of the simulation and never rise to the given values. I thought it might be because the normal pressures don’t match the relevant strain state, thus I ran a simulation with NVT first and used the final pressures from the NVT simulation to apply the NPT. Nonetheless, the problem resists. Going through the mail list archive and LAMMPS documentation I didn’t find something helpful. Could anyone please give me some ideas or references that I can turn to?

Sincerely,
Joe

2

Since fix npt does work well for “normal” systems and when used correctly, the problem is most likely caused by how you have set up your system or incorrect settings in your input.
Axel

3

Dr. Axel,

Thanks for your reply. I agree with you that something is wrong.

  1. In terms of the domain, I pre-compressed my domain using the info. from the bulk simulation including strain, stress. And to get the velocities, I relaxed the defective domain at 300K and 0 pressure using NPT, then I got the atomic velocities. Checking the potential energy of each atom at the 0ps, I don’t observe a concentration of the potential energy which may indicate the domain is stable.

  2. There maybe some stupid mistakes in my input file. However, I can’t find them out. If you can, could you please give it a glance?

INITIALIZATION

units real
dimension 3
processors * * *
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic

box tilt large

read_data data.txt

pair_style born/coul/long 15.0
kspace_style pppm 1.0e-4
include pair_born.dat
special_bonds lj/coul 0.0 0.0 1.0 dihedral yes

suffix off
kspace_style pppm 1.0e-4

neighbor 2.0 nsq
neigh_modify delay 0 every 1 check yes one 100000 page 1000000

compute 1 all temp
compute PE all pe/atom
compute KE all ke/atom
compute eng all pe/atom
compute eatoms all reduce sum c_eng
compute stress all stress/atom NULL

variable totaleng equal “c_eatoms”
variable p1 equal “-pxx0.000101325"
variable p2 equal "-pyy0.000101325”
variable p3 equal “-pzz0.000101325"
variable p4 equal "-pxy0.000101325”
variable p5 equal “-pxz0.000101325"
variable p6 equal "-pyz0.000101325”

suffix gpu

pair_style buck/coul/long 15.0
include pair_buck.dat

reset_timestep 0

variable stress_xx equal “-2.6176e+04”
variable stress_yy equal “-5.6426e+04”
variable stress_zz equal “-7.3006e+04”

fix 1 all npt temp 300 300 1000 x {stress_xx} {stress_xx} 100000 y {stress_yy} {stress_yy} 100000 z {stress_zz} {stress_zz} 100000

thermo 100
thermo_style custom step pe ke lx ly lz press v_p1 v_p2 v_p3 v_p4 v_p5 v_p6 vol v_beta temp

dump 1 all custom 1000 deform.* id mol type q xu yu zu c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] vx vy vz ix iy iz c_PE c_KE

run 100000

undump 1
unfix 1
unfix 2

Sincerely,
Joe

4

Dr. Axel,

Thanks for your reply. I agree with you that something is wrong.

  1. In terms of the domain, I pre-compressed my domain using the info. from the bulk simulation including strain, stress. And to get the velocities, I relaxed the defective domain at 300K and 0 pressure using NPT, then I got the atomic velocities. Checking the potential energy of each atom at the 0ps, I don’t observe a concentration of the potential energy which may indicate the domain is stable.

  2. There maybe some stupid mistakes in my input file. However, I can’t find them out. If you can, could you please give it a glance?

this is not very helpful since it is missing files that are read by this input. having the log file would be more helpful.
from looking through the input, there are a few things that stand out. not all of them related to your issue. please see my comments following those lines.

INITIALIZATION

units real
dimension 3
processors * * *
boundary p p p
atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style harmonic
improper_style harmonic

box tilt large

is this really needed? it is not a good idea to do an NPT simulation of a strongly tilted simulation cell. the larger the tilt the more likely that there will be a tendency to “flatten” the cell (as that will become a favorable configuration, especially with long-range electrostatics). in general, it may be better to construct a less tilted or even orthogonal supercell.

read_data data.txt

pair_style born/coul/long 15.0
kspace_style pppm 1.0e-4

1e-4 is a very low convergence. while it is acceptable for constant volume simulations, because there is a lot of error cancellation, it will introduce significant errors for variable cell calculations since stress components are not so easily converged. something like 1e-6 would be better.

include pair_born.dat
special_bonds lj/coul 0.0 0.0 1.0 dihedral yes

why “dihedral yes”?

suffix off
kspace_style pppm 1.0e-4

neighbor 2.0 nsq

why “nsq”? this is rarely an efficient choice unless your system is tiny. and then you may be better off using ewald instead of pppm.

neigh_modify delay 0 every 1 check yes one 100000 page 1000000

why so large values for “one” and “page”? they are like 10x larger than the default? for a normal density system, this should not be needed.

since you seem to be using the GPU package, that setting is irrelevant anyway.

compute 1 all temp
compute PE all pe/atom
compute KE all ke/atom
compute eng all pe/atom
compute eatoms all reduce sum c_eng
compute stress all stress/atom NULL

variable totaleng equal “c_eatoms”
variable p1 equal “-pxx0.000101325"
variable p2 equal "-pyy0.000101325”
variable p3 equal “-pzz0.000101325"
variable p4 equal "-pxy0.000101325”
variable p5 equal “-pxz0.000101325"
variable p6 equal "-pyz0.000101325”

suffix gpu

pair_style buck/coul/long 15.0
include pair_buck.dat

reset_timestep 0

variable stress_xx equal “-2.6176e+04”
variable stress_yy equal “-5.6426e+04”
variable stress_zz equal “-7.3006e+04”

fix 1 all npt temp 300 300 1000 x {stress_xx} {stress_xx} 100000 y {stress_yy} {stress_yy} 100000 z {stress_zz} {stress_zz} 100000

thermo 100
thermo_style custom step pe ke lx ly lz press v_p1 v_p2 v_p3 v_p4 v_p5 v_p6 vol v_beta temp

dump 1 all custom 1000 deform.* id mol type q xu yu zu c_stress[1] c_stress[2] c_stress[3] c_stress[4] c_stress[5] c_stress[6] vx vy vz ix iy iz c_PE c_KE

run 100000

this run is far too short to have an equilibrated system given the time constants of your fix npt command.

axel.

5

please note that your final fix npt command after expanding all variables looks like this:

fix 1 all npt temp 300 300 1000 x -47838.4766090731 -47838.4766090731 100000 y -103121.901165786 -103121.901165786 100000 z -133423.242804235 -133423.242804235 100000

this means you have a very large negative pressure in all three directions, which translates into the system wanting to expand forever.
so whatever you did to determine those pressure values, it was wrong.
given the size of your system it is not clear to me why you cannot use a supercell and an orthogonal cell. your fix npt command is chosen to preserve the tilt anyway, so using an orthogonal cell will make things much easier and also numerically more stable and faster, as using a triclinic cell requires that during the run the system is converted multiple times from regular cartesian coordinates to fractional coordinates and back.

in short, LAMMPS is doing what you are asking for.

axel.

6

Dr. Axel,

What would be your suggestion for the run time? Say, 1000000?

you have to run an equilibration simulation for as long as it takes until the system is equilibrated.
there is no single reliable indicator since nobody can know the shape of your potential hypersurface and how long it will take to explore all relevant parts of it sufficiently well.
if your system can have “rare events” it can take a long time. typically for something like your simulation one would look at running averages (with a suitable but not too short and not too long window) for properties like total energy or volume. or structural properties like radial distribution functions. once they are converged to a sufficient degree, you could stop your calculation.

axel.

7

Dr. Axel,

I appreciate your patience to look into the log file and explain to me. However there are still a few questions remaining.

  1. I do have very high pressures in all three directions. But I got them from a simulation using NVT to hold the temperature and keep the volume. So in my opinion, theoretically applying those pressures and corresponding strain, the NPT ensemble should be able to hold the pressures. Therefore I don’t quite understand why “LAMMPS is doing what you are asking for.” Could you please give me a bit more explanation? If it is wrong, how can I get them correctly?

  2. The reason I use a triclinic cell is that 1) the system I’m simulating is triclinic with a large tilt factor; 2) LAMMPS documentation says that my tilt factor should not be larger than the x-dimension of my domain; 3) the x-dimension of my system is small. So that it seems that the system is almost orthogonal. But LAMMPS documentation does say that I can set all tilt factors as 0 initially, I may try it.

Sincerely,
Joe

8

Dr. Axel,

I appreciate your patience to look into the log file and explain to me. However there are still a few questions remaining.

  1. I do have very high pressures in all three directions. But I got them from a simulation using NVT to hold the temperature and keep the volume. So in my opinion, theoretically applying those pressures and corresponding strain, the NPT ensemble should be able to hold the pressures. Therefore I don’t quite understand why “LAMMPS is doing what you are asking for.” Could you please give me a bit more explanation? If it is wrong, how can I get them correctly?

you are confusing “internal” and “external” pressure and need to refresh your understanding of some basic physics. when your system has a negative pressure it “wants to shrink”. if it has a positive pressure it “wants to expand”. however, the pressure you apply with fix npt is an “external” pressure. if that pressure is positive, it will try to compress the system, if it is negative, it will expand the system. since the pressure at inifinite dilution is zero, a negative external pressure will always expand your system. even though your initial geometry has a negative pressure, as your system relaxes and equilibrates it should relax to become less negative. thus you can have a balance between internal and external pressure only for positive pressure values.

the fact that you have a very high initial pressure and also that a rather large amount of kinetic energy is released in the first part of your simulation (the temperature rises quite significantly) are both an indication that your initial geometry is not very consistent with the force field that you are using. it is not possible to say whether this is because of a bad geometry or bad potential parameters or an incorrect topology.

  1. The reason I use a triclinic cell is that 1) the system I’m simulating is triclinic with a large tilt factor; 2) LAMMPS documentation says that my tilt factor should not be larger than the x-dimension of my domain; 3) the x-dimension of my system is small. So that it seems that the system is almost orthogonal. But LAMMPS documentation does say that I can set all tilt factors as 0 initially, I may try it.

any triclinic system can be represented as an orthogonal supercell. depending on the tilt that supercell may be very large, though. but at some point the errors will be small compared to the error from using empirical classical potentials instead of more accurate force computation methods.

axel.

9

Dr Axel,

Thanks for your explanation. It is clear to me now. Appreciate a lot!

Sincerely,
Joe

10

Dr. Kohlmeyer,

According to our previous communication, I improved and modified the setting of my simulation. However, except temperature and Kinetic energy, other properties of my system cannot equilibrate. And I’ve run it for a reasonably long time(800000 timesteps, 800 picoseconds). Does that mean my system won’t be able to equilibrate?

Sincerely,
Joe

11

Dr. Kohlmeyer,

According to our previous communication, I improved and modified the setting of my simulation. However, except temperature and Kinetic energy, other properties of my system cannot equilibrate. And I’ve run it for a reasonably long time(800000 timesteps, 800 picoseconds). Does that mean my system won’t be able to equilibrate?

sorry, but i do not own a crystal ball, thus I cannot answer that kind of question.

axel.

12

Dr. Kohlmeyer,

My system is too big to be sent by email. Still appreciate your help.

Sincerely,
Joe

13

Dr. Kohlmeyer,

My system is too big to be sent by email. Still appreciate your help.

i could not tell if it was smaller from just looking at files and i don’t have the time (and interest) to figure this out on your behalf. sorry.
this is something that has next to nothing to do with LAMMPS but rather is a matter of general MD expertise. you need to talk with your adviser/supervisor or experienced colleagues or collaborators for guidance.

axel.

14

Dr. Kohlmeyer,

Thanks for your suggestion.

Sincerely,
Joe