Questions regarding the boundary conditions and pressure calculations in 2D materials

Hello LAMMPS Community,

I’ve got some questions about 2D materials, specifically crystalline nanosheets with a thickness of up to 3 atomic layers along the Z direction, similar to graphene. I’ve tried searching through the archives, but I couldn’t find answers that directly address the following questions:

  1. In the literature, I’ve noticed that most simulations involving very large sheets use full periodic boundary conditions in the XY plane and add a sufficiently large vacuum space along the Z direction to prevent interactions between periodic images. I’m thinking about trying shrink-wrapped boundary conditions along the Z direction while keeping everything else the same. Would this approach give me a material state and physical properties very similar to the traditional method, except for potentially better control over load imbalance with the ‘boundary p p s’ setting? Can you confirm if my understanding is correct?

  2. I came across Axel’s explanation regarding the fix nve, nvt, and npt commands, which suggests that none of these fixes explicitly enforce the associated ensembles. Instead, they mainly focus on integrating particle velocities and positions. For 2D materials, whether with full periodic boundaries or ‘boundary p p s,’ it seems that none of these fixes can mimic the corresponding ensembles. Is my understanding accurate?

  3. I’m finding it challenging to calculate pressure or stress components in the periodic directions within the XY plane. I’m not even sure if the stresses reported are as meaningful as those in bulk materials. Could you provide guidance on the best practices in LAMMPS for volumetrically relaxing 2D materials in the XY plane under periodic boundary conditions? If I choose ‘fix npt+vacuum,’ I’m concerned about artificially increasing the area for pressure calculations along the X and Y directions. On the other hand, if I choose ‘boundary p p s,’ I’m unsure about the physical reliability of the calculated pressures along X and Y directions when using fix npt for materials with very few atomic layers. Any insights on this matter would be greatly appreciated.

I greatly appreciate your time and assistance in addressing these questions in advance!

I am no expert, but I think the key thing is to think that the fix nvt and fix npt calls some specific equation of motions (e.g. Nose Hoover). If you read the official publication, the Nose Hoover equations of motion were proved to be able to reproduce the ensembles in the sense that they were able to reproduce the probability density function of the given ensemble (if I remember correctly that was it). Other equations of motion, developed by other people, can also ensure a fix NVT condition (or fix NPT also) but should not reproduce any ensemble (for example, it seems (I’ve read) that Langevin is one of these cases, although I never really read the official Langevin publication to be honest).

I think you cannot say for granted you are reproducing an ensemble when you are using these fix nvt, npt, nve commands because you may be adding additional features to the dynamics that would change it compared to using the original underlying equations of motion on their own.

Why do you say that for 2D materials the equations wouldnt reproduce the ensemble? (I am asking just because I am curious)

A not fully periodic system is an “open” system with effectively an infinite volume. For that reason, neither global volume nor global pressure are well defined. LAMMPS well still output the volume of the cell, but that is not the volume of the system, only the volume that LAMMPS is “monitoring”. Since pressure depends on volume, that is also tainted. You can still apply fix nvt, nve, npt, but in the latter case you can only monitor Pxx and Pyy if z is non-periodic and adjust only those components of the cell size. That means using integrators, thermostats, and barostats is still valid but you will not represent the ensembles of their names since the system does not conform to how the ensembles are defined.

ahhh, I see, makes sense. I thought it was not possible to reproduce any ensemble for 2D materials in general (even if it is a case where the “nanophase” is surrounded with another phase instead of being in vacuum or so). Thankssss!

To give my own answer to your three questions:

  1. Shrink-wrapping the z-direction is unlikely to give you a more efficient simulation:

    • If your force field includes partial charges and long-range electrostatics, then you cannot use shrink-wrapping at all (you must either have periodic z, or fixed z with a slab correction).
    • Otherwise, you still need enough distance between your 2D material slabs to remove interactions (since you are simulating a 2D material, not a 3D stack of material). So why not just fix the z-boundary? If you expect particles to “boil off” your material so that you need the z-boundary to change then that is quite an unusual 2D material indeed.
    • An immediate way to control load imbalance is to set processors * * 1. This will create a 2D x-y grid of processors, each of which should have roughly the same number of particles.
  2. It is hard to give specific further insight without knowing which specific instance of “Axel’s explanation” you are referring to (half the site is his explanations after all, given his job). As he has explained, you can measure the xx and yy components of the pressure tensor in a 2D simulation and use them to just change the x and y dimensions of the box (but see my next point). This is not an “NPT ensemble” from the strict point of view of an ensemble in equilibrium with a three-dimensional pressure “reservoir”, but it is certainly a well-defined simulation with a reasonable interpretation of results.

  3. You could certainly do something like

fix npt all npt temp 300 300 $(100*dt) x 1 1 $(1000*dt) y 1 1 $(1000*dt) couple xy

(check documentation for full details). But honestly I would ask someone with DFT expertise to tackle this, unless the force field you use is well-validated enough for you to draw quantitative conclusions (in which case there would be a status quo in terms of the right way to use this force field, and you would be expected to follow it for correct results).

1 Like

@akohlmey @ceciliaalvares

Thank you both for these fruitful discussions. I always learn something new from them.


Thank you for your elaborative comments and suggestions, as always. I’m interested in the mechanical properties of the material, but still struggling with the terminology of stress in 2D materials. In my opinion, the term is somewhat abused because there is no well-defined area normal to the edges of these materials. To partially correct this definition, people use stress with force per length unit.

Regarding LAMMPS, if I want to conduct a thought experiment, I can enlarge the vacuum as much as I want. If I’m not mistaken, per Axel’s comments, the pressure should drop as LAMMPS monitors the predefined volume and area rather than the actual volume of the system. In that sense, do you think there is no difference in using full periodic with vacuum and ‘boundary p p s’ to take care of the lateral stress or pressure components in the XY plane?

And why not just do it? There’s no reason for you to do a “thought experiment”. You and I are not Albert Einstein, who had to use the power of his imagination to come up with special relativity because there’s no such thing in real life as a train that runs at the speed of light.

LAMMPS is right there for you to use. Set up some simulations and run them and see what result you get. In particular, both Axel and I have talked to you about the concept of the pressure tensor and its xx and yy components. This is standard materials mechanics that you can learn about in most textbooks. They are how we can work with the concept of a 2D material and its 2D strain under various 2D stresses.