Stress concentrations in elongated sample in tension (p p p)

Hi Liam, first, let me clarify that I only guessed something about the definition of the virial for many-body potential.

If I do what you are supposed to do, i.e. look up the doc page:
I find the following:

Details of how LAMMPS computes the virial for individual atoms for either pairwise or many-body potentials, and including the effects of periodic boundary conditions is discussed in (Thompson). The basic idea for many-body potentials is to treat each component of the force computation between a small cluster of atoms in the same manner as in the formula above for bond, angle, dihedral, etc interactions. Namely the quantity R dot F is summed over the atoms in the interaction, with the R vectors unwrapped by periodic boundaries so that the cluster of atoms is close together. The total contribution for the cluster interaction is divided evenly among those atoms.

So my guess was inaccurate after all, and you should probably read the above article for more information pertinent to your model. Since in this case the code implementation is well described in a paper, that should be the reference more than the mailing list.

I wrote |\vec{r}| in TeX notation; if you are not familiar, paste the code between $'s in one of the many webpages like this:


Thanks for passing the reference along. I do not invoke the virial key word, so I am also including the KE terms. I can try with the virial keyword.

Based on what I am finding, this appears to mean that within a periodic box the average stress in a specific region can be different than in another region. I am discovering higher stresses near the root as opposed to the middle. This surprised me as I used periodic to avoid free surface and edge effects. I have not seen any reference to this so I wanted to be sure it was not an error in my code or within LAMMPS.


I don’t think there should be different stresses in different parts of a homogeneous system (if I understand your statements correctly). The different modulus you observed might be caused by insufficient time average. What are the uncertainties (statistical errors) in your data?


Giacomo Fiorin <giacomo.fiorin@…29…> 于 2018年12月6日周四 21:15写道:


I did not expect it either, leading me to believe it is likely as issue with something in my procedure.

I am equilibrating for 20 ps, and the structure appears well equilibrated. I then strain at a strain rate of 0.01 per picosecond, a rate I have seen in literature. My modulus in the middle 80% is about 140 while the root 10% is about 170. What is odd though is that if I make the root area even smaller (5%) I see even higher concentrations at the root, lending itself to the conclusion of concentrations near the root. Overall though, the average modulus for the three regions is very close to the modulus predicted when using the total stress (sum of the per atom stresses divied by volume). My input script is attached!

Thanks for your help.


saintvenant.txt (4.28 KB)

Have you printed out the volumes of the different groups and checked that they are what you intended?

There should be no problem from the lammps side. It must be related to some details of your simulation.



Morrissey, Liam S. <lsm088@…7202…> 于 2018年12月6日周四 21:40写道:


To make sure the volumes were ok I actually just defined them as constants which were 20% of the total volume, 80% and then 20% for root, middle, tip volumes respectively. During straining the length will change but volume of each region will not so I think this is ok. To make sure my numbers were correct I also compared the total pressure in the x (pxx) with the total stress in the x as calculated from: ( sum the stress per atom in the [1] direction)/volume. These two values were always equal, as they should be based on the manual.

During straining at each step I then calculated the stress in each group by:

  • summing the stress per atom in the [1] direction for each group (root, middle, tip). For this compute reduce the groups were defined from a region definition pre straining that was 0-20, 20-180, and then 180-200 of the length.
  • and dividing this total stress by the volume of the respective group
  • finally, for each step I could check the total pressure in the x (pxx) against the average stress using the stress from each region: (.2stressroot+.8stressmiddle+.2*stresstip). These values were within 2% error at each step. This makes me think that the numbers are ok.

For example, in the first step of strain pxx was -0.13, stress root was -0.061, stress middle was -0.167, stress tip was 0.067.

When I then get the modulus for each region I am see a higher modulus at the root than the middle.

I am really having trouble seeing an error in my methodology here and am confused why the stress at the root, middle and tip when divided by their respective volumes are not the same.

All this above makes me feel like my methods are good, but I am confused as the manual and common sense says I should not see stress concentrations at boundaries when periodic.

Thanks and sorry for the long explanation.

I can only imagine two psssible problems, one is incorrect volume, the other is insufficient statistics. It seems you are sure the volumes are correct. Then the second question: did you get similar results in different independent runs ( with different random velocity seeds)?


Morrissey, Liam S. <lsm088@…7202…> 于 2018年12月7日周五 17:14写道:

Being that these values for the stress are contrary to what you were previously reporting, I’d bet these are due to fluctuations - be it the actual value of the stress or even number of particles in your volumes. Its also consistent with your observation that using a smaller volume produces ‘larger concentrations’ -> fluctuations.

On a more personal level. From the small amount of observations you’ve provided, the values in the tip and root are also fairly close and substantially different from the middle. I’d probably also check that there are no low frequency waves in the system.

Hi Eric,

Thanks for the response. I should have clarified, initially I was reporting the elastic modulus from each section, in my most recent email I simplified things and just reported the stresses at each section for a strain increment. Ovearll the elastic modulus for the root volume is higher than the middle. However, after your suggestion to see if the number of atoms in each section are fluctuating I re-ran the simulation while tracking n atoms in each region:
1981 15780 1979

1912 15679 1978

1861 15653 1978

1836 15679 1998

My simulation began with 20000 atoms but for some reason I am losing atoms during the strain (I did not get a lost atoms error). I think this may be the source of the error. However, even when I look at the stresses pre straining and after equilibrium they are not evenly distributed.

However, pre my NPT equilibration the stresses were fairly even, and became successively worse after NPT to 300K, making me think the way I am defining regions/volumes is the problem. Perhaps the atoms are moving from their region and because they region is not dynamic it is making it appear like a concentration?

Stress root Middle Tip (Pa)

2264.645 2205.3 2280.68
-73.097 -2415.64 265.139
-237.312 412.0085 -250.507
318.2054 1029.441 112.1347
111.66 2872.638 -504.999
154.2957 1733.8 -345.347
70.912 -997.353 260.7136
444.6972 -2071.95 890.4407
237.9774 -1416.84 631.7219
162.3972 1121.202 -45.7307
-48.245 3171.495 -449.219
14.55666 -3271.69 91.44401
-148.958 -2843.78 384.6902
-196.273 2100.708 -253.064
89.01935 655.7885 -404.937
-16.892 3877.771 -672.246
-277.099 -40.0287 -505.547
-332.326 -2333.35 -3.34882
162.2108 -1102.77 640.3921
268.0092 -1657.26 615.848
-45.0259 1994.259 -363.071

I have attached the log file to see the actual numbers.

log.lammps (11.7 KB)

Sounds like you’re on the right track then. :slight_smile:


I am still a little surprised that the number of atoms in each group is changing.I was under the impression that when a group is defined those atoms wouldn’t change as groups are not dynamic unless I specify the dynamic key word. So to see number of atoms in each group changing during straining is not really following what is indicated by the manual.

Also, the fact that the stresses are fairly even pre npt indicates a lot of atom motion between groups during temp equil. Is there anywhere in my script that seems like an error to explain this?

Perhaps I just shouldn’t be defining groups before straining and then attempting to look at stresses of these individual groups because a lot is going on within the simulation box during the deformation. All the values match when I look at the whole box, making it seem like it really is a problem with groups changing during deformation.

As a brief follow up, I am now just focusing on the equil, to simplify further and not worry about straining. When I start my equil there is 2100 atoms in the root, 16000 in the middle, and 1900 in the tip, totally 20000 as per the log file. After 1000 steps it is 1944 in the root, 14643 at the middle, at 1801 at the tip. I have done other studies using the group function and once assigned to the group the number shouldn’t change unless it leaves the simulation box? Why are they changing and why is there a net reduction in atoms?

saintvenant.txt (4.36 KB)

I can’t speak for the pressure issues you have, but the number of atoms remains at 20000 based on your logfile:
Loop time of 130.477 on 1 procs for 2000 steps with 20000 atoms (Line 223, end of the simulation)

Your method of computing number of atoms in your group is wrong. As stated in the documentation, the group assignments do not change as atoms move out of/into the region. So as the simulation progresses, atoms from outside the region come in, displacing atoms in the region. This doesn’t change the number of atoms in the group. However, your method of counting (e.g., count(root,root) ) counts the number of atoms that are in the group root AND in the region root. So that count will go down as the system relaxes.