Pressure in non periodic boundaries.

LAMMPS LAMMPS Mailing List Mirror
1

Dear Lammps Users,

I am simulating a crystal by using periodic boundary conditions in x and y directions and free surface with a vacuum above and below in z direction. I am using NVT for equilibrating the crystal at 300 K. As there is vacuum above and below the crystal (in Z direction), I thought the crystal will expand (to achieve mechanical equilibration) and pressure will be zero during equilibration. But I get a very high value of pressure. I also see high pressures during heating simulations in NVT and little/no expansion.

1 .Is there a physical reason why crystal is not expanding or is the pressure calculation not correct?

  1. In general, I am wondering what would be the best method and ensemble to simulate melting of nanoparticles and thin films (with free surfaces in one direction).

After reading the lammps mailing list archives, I come to know that under non-periodic boundaries, the volume is effectively infinity and doesn’t change. what LAMMPS uses as volume of the simulation cell is truncated completely arbitrarily for technical and convenience reasons.

  1. But I would like to know why volume should be infinity, given that pressure of crystal should be based on volume of crystal itself and the simulation domain volume is finite

I would be very thankful if someone can help me to clarify my doubt.

Thanks

Nikhil Joshi

2

Here are some responses to your questions:

  1. [Self-plagiarized from recent post on same subject] Consider a single pair of atoms, 1 and 2, in vacuum, interacting via a harmonic spring. The virial stress tensor is proportional to F_{12} x R_{12}. At any instant, this can be positive or negative, depending on whether the spring is in tension or compression. We can only expect that the average value will be close to zero, and even then, there is lots of stat mech fine print, but for large systems, it is true that the non-periodic components of the virial will be approximately zero (when the kinetic contribution is added in). As you average over larger and larger surface areas A, and/or longer and longer times T, the magnitude of the fluctuations in the virial should decrease roughly as 1/(AxT)[End of self-plagiarizing]
    So you should expect to see the Pzz undergo large fluctuations about zero. If this is not happening, something else is going on (external forces, insufficient equilibration, insufficient time averaging)

  2. This is an impossible question to answer. It depends on what you are trying to accomplish in the melting. If you are interested in realistic kinetics, you can not just superheat a perfect crystal, even at a surface. That will greatly underestimate the melting rate, due to absence of nucleation sites. Read the literature on MD simulation of melting, particularly laser surface heating.

  3. This is a matter of interpretation. LAMMPS will use whatever volume is defined by the simulation cell boundaries i.e. lxlylz. It is up to you, as the scientist, to decide what you think is an appropriate volume. The key point here is that pressure is a thermodynamic or continuum mechanical property that requires the definition of a material volume or at least a material surface. In MD, which deals with points masses, there are no surfaces or volumes, except those artificially imposed by the simulation cell. In the case of periodic boundaries, these have an important physical meaning, in the case of a non-periodic system, not so much.

Aidan

3

Thank you for the quick response and helpful comments. I did not decompose the pressure previously. When I decompose the pressure into Pxx, Pyy, Pzz, I observe that Pzz flucutates around zero which is expected as there is vacuum above and below the crystal in z direction. But Pxx and Pyy are as high as 20000 bars when NVT ensemble is used. We did apply periodic B.C. in other two directions (X and Y). Of course, Pxx and Pyy did fluctuate about zero when NPT ensemble is used and when barostat is applied in X and Y directions.

(1) Is this trend for NVT case normal? After melting we did observe that the crystal expanded to fill in the vacuum space and Pxx and Pyy dropped to zero.

(2) We did read a paper titled ‘‘Calculation of pressure in case of periodic boundary conditions’’ by Manuel J. Louwerse , Chemical Physics Letters,April 2006, Pages 138–141, http://dx.doi.org/10.1016/j.cplett.2006.01.087, which says that virial equation state cannot be used for non-pairwise additive force fields and that a correction must be applied. We are using Vashishta’s potential function for alumina (Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina J. Appl. Phys. 103, 083504 (2008); http://dx.doi.org/10.1063/1.2901171). Has LAMMPS applied this correction?

(3) We are bit concerned with high pressures, since we are interested in determining melting points at 0 bar pressure. It appears that NPT is suitable, but we are not sure if the pressure is calculated correctly.

Thanks again for all the help. I greatly appreciate it.

Thanks
Nikhil Joshi

4

Re: (2) LAMMPS has many non-pairwise potentials,
which it computes the viral for. Not sure what additive
means in this context. The vashishita pair style
in LAMMPS is relatively new. You can contact the
authors (top of file) if you have a Q about its viral.

Steve

5

Here are responses to your numbered questions:

  1. Yes, this is normal. The behavior you describe in the first paragraph is exactly what you should expect for NVT MD of solid with free surfaces. Essentially, you have equilibrated the solid at a stress/strain state of exx != 0, eyy 1= 0, Pzz=0, where exx and eyy are the axial strains relative to the the strain states for which Pxx=0, Pyy=0. This corresponds to a combination of hydrostatic compression/tension plus non-zero shear stress. When you melt the crystal, the shear stress is able to fully relax, and you end up with Pxx=Pyy=Pzz=Phydro. If you give the free surface enough vacuum to expand into, Pzz will be able to relax to zero, and so Phydro=0. Even if the free surface is fully consumed by the expansion, so the volume relaxation is incomplete, Phydro will be much smaller than the original Pxx, Pyy.

  2. It always surprises my when people suspect that they have discovered a secret error in LAMMPS, when other explanations are much more likely. I am familiar with the Louwerse paper. In fact, their misguided claim that the virial equation for periodic systems is somehow in need of correction motivated us to write our own account of how LAMMPS evaluations the virial in periodic systems. There is no need for corrections. There is simply one correct answer and multiple ways to calculate it exactly. LAMMPS supports three distinct ways. [Thompson, Plimpton, Mattson, J Chem Phys, 131, 154107 (2009).] Rest assured, what LAMMPS calls “press, pxx, pyy, pzz, pxy, pxz, pyz” in thermo output is the correct pressure and stress, regardless of what pair_style you are using. All this and more is described in great detail on the doc page for compute pressure. I suspect you have not read it.

  3. You can certainly use NPT in x and y to relax Pxx and Pyy, but you should be careful. It is very easy to generate large non-equilibrium volume changes when you run NPT MD starting at a far-from-equilibrium initial state.

I will also mention that you might find it scientifically useful to eliminate the issue of nucleation kinetics by starting the system out with the top half (z>0) of the material melted and the bottom half (z<0) solid. If you adjust the temperature, you will be able to maintain the system in the two-phase state for extended periods of time. The extensive body of work by J.J. (Jeff) Hoyt and collaborators is also very relevant.

Aidan