A problem in setting or choosing an integrator in the uniaxial strain of 2D structures

When subjecting a two-dimensional structure to uniaxial strain in some direction, I have had problems.

The problem is to strain one direction of a 2D structure. I used a fix nvt to perform this action, together with a fix deform in the x direction, for example. The result is exciting, and the structure breaks as expected. However, when checking the stress-strain curve, in the y direction, there is a stress similar to the tension in the strain direction. This is because the box does not adjust in the y direction, and there is stress in both directions. When studying a “more brittle” structure, I noticed that it would break in y when stretched in x because of this stress in the direction perpendicular to the applied strain.

To solve this, I started using an NPT integrator instead of the NVT. When using the fix npt, I realized that the box collapses, depending on the choice I allow regarding the directions of the npt. For example, if I use fix npt … x 0 0 0 {pdamp}, with deform x, it gives an error. If I use fix npt … y 0 0 pdamp z 0 0 pdamp, the box collapses at z. Switching to fix npt … y 0 0 pdamp, the box collapses rapidly at y, causing the structure to lose its two-dimensional topology. I tried adding a drug 1.0 to NPT, but that didn’t solve the problem either.

Note: I have always used a pdamp of 1000*dt. I don’t know if this influences much.

I’d fix the z direction at some distance and only barostat the y dimension if it is a 2D material in xy. That should prevent the z dimension from collapsing and any periodic behavior that may arise from too small of a box in that direction.

NVT will create an off axis stress if the 2D material is periodic in 2D. NPT will not. Pdamp of NPT will eventually influence the behavior of the simulation, but only at very high values, by effectively converging on the behavior of NVT. See our preprint for more details on that
https://papers.ssrn.com/sol3/papers.cfm?abstract_id=4289377

In my 2D system, the sizes are approx. 250x250 Ang^2, periodic in x and y. The box size in z I left at 200 Ang. Could this size generate the collapse? Still, if the thermostat is only in y, the box does not collapse in z. But, I’m not sure what to do to avoid the collapse itself in y.

Thanks for the preprint. I will look at it carefully, but I can already see that it describes many important aspects and will definitely be quoted in my future work.

Unless your nonbonded cutoff exceeds 200 Ang, I don’t see why that z dimension length would cause any issues. The collapse in the z if barostating it would be the barostating function trying to match the prescribed value. If it is fixed and not used, it shouldn’t cause any issues ( boundary p p f ).

For an xy periodic 2D system, the integrators could be one of the following

boundary p p f

# NVT y
fix 1 all nvt temp ${T1} ${T1} ${tdamp}
fix 2 all deform 1 y erate ${rate}

# NPT y
fix 1 all npt temp ${T1} ${T1} ${t_tau} x 0.0 0.0 ${p_tau}
fix 2 all deform 1 y erate ${rate}

# NPH y
fix 1 all nph x 0.0 0.0 ${p_tau}
fix 2 all deform 1 y erate ${rate}

# NVE y
fix 1 all nve
fix 2 all deform 1 y erate ${rate}

For the collapsing in y, what is your input deck? Would you mind posting it here?

What happens if you don’t use fix deform, but you use this fix npt in the y-direction only? Are you sure that your system is stable at 0 pressure in that direction?

Antes de fazer o fix deform eu realizei 250k passos de npt iso, e também outros 250k passos de nvt.

Then I agree with TJ, we’ll need to see your full input deck to get any further.

Basically, that’s it!

fix     npt_structure all npt temp 300.0 300.0 10 y 0.0 0.0 100
fix     strain all deform 1 x erate 1e-6 remap x
# -----------------------------------------------------
run     7000000
# -----------------------------------------------------
unfix   strain
unfix   npt_structure

What does your data look like? If y 0.0 0.0 100 is compressing in the y direction, then there probably aren’t bonded / non-bonded forces resisting the dimension change in compression

edit: And by data I mean the atomic structure, bonds, angles, etc.

Look at this preprint (https://arxiv.org/pdf/2007.07168.pdf) please. In it, in figure 04, the same problem occurred. The more porous system allows for even more collapse in the direction perpendicular to the strain.

When we ask for an input deck, we mean the full input, the data files needed to run the input, and the terminal command used to execute the program. The log file is also useful.

This is important because it saves us from having to ask many tedious questions like “What units are you using? Are you using fix enforce2d?”, etc.

What do you mean when you say this? Does it gain z-height? Or does it simply break down?

Right. I will prepare the complete deck to send.

As for your last question: Yes, the structure gains height in z, and it is no longer a two-dimensional structure. While in the opposite direction of strain application, the box drastically decreases in size. In figure 4 of the preprint I sent above, I show a result figure with this.

As far as I can tell, LAMMPS is working as intended. As you increase the strain in the y-direction, you observe the stress increase in a predictable way. The only deviation from your expectations is the failure mode, but crumpling like this is a valid result. From here on it is your responsibility to interpret the results and understand them enough to prove to your reviewers that they are valid.

As there is no issue with LAMMPS, but with your expectations, this is now off topic for this forum.

As some general advice, I would be careful of how the simulation box is divided up among the processors. Your timing shows wildly varying times (between 8-800s for pairwise computation in the NVT section!). Check out the balance and processors commands, consider shrinking the z-height of the box, and maybe reduce the number of processors (64 may be unnecessary for 12.7k particles).

1 Like