Uniaxial behavior change between LAMMPS Sep 2021 and Jun 2022

Hi LAMMPS users,

I noticed a change in behavior between a data set run with the stable 29 Sep 2021 build of LAMMPS and the stable 23 Jun 2022 build, with supporting images in the pdf.

LAMMPS 2021 2022 Stress Strain Behavior Change.pdf (318.4 KB)

The set up is a graphene sheet with a hole in it, shown in Figure 1, sampled at 100 points in time after a 5 ns equilibration. Figure 2 contains the first ten restarts for each set of damping parameter combinations, with Figures 2B and 2C being the same input (below) on the same data using the June 2022 and Sep 2021 LAMMPS builds, respectively. Figure 2A is an NVT run and Figure 2D an NPT with the standard ‘rule of thumb’ damping parameter values for reference. The behavior of the Jun 2022 build is consistent under slightly different (but effectively the same) damping parameters, with the transverse stress stable at some value around 0 GPa (see extra image below), compared to the relatively volatile transverse values returned by the Sep 2021 build.

NPT 100 1001 temp 2022

Figure 3 plots both behaviors as failure probabilities, with the transverse behaviors in Figure 3B sampled at the ultimate axial stress value used in Figure 3A. These distributions use all 100 restarts generated from the graphene sheet.

I’ve looked through the change logs and related src files on github to see if there were any changes that might lead to this difference in behavior, but I could not find anything that stood out. All restarts, starting data, and the REBO potential used were either developed using, or from, the 29 Sept 2021 build, so it could be a backwards compatibility issue on my part that I missed when looking at the documentation.

Does anyone have any ideas as to what might be causing this?

Thanks,
Tj

Edit: Added a restart file and REBO potential
LuoArmchairSingleCrystal-a0-10-r02_1.rst (8.5 MB)
CH_Luo.rebo (905.0 KB)

Tensile code for reference, with command line input:

# Command line input values
mpirun -np 6 --bind-to socket -x OMP_NUM_THREADS=2 lmp -sf hybrid gpu omp -pk gpu 1 gpuID 0 -in 20220517_Graphene-Tensile.in


#-----------------------------------------------------------
# Simualtion Variables
#-----------------------------------------------------------
# Temperature
variable T1 equal 300

# Timestep
variable tms equal 0.001 # ps
variable tdamp equal ${tms}*100.0
variable pdamp equal ${tms}*999.0
#-----------------------------------------------------------------
# Graphene tensile
#-----------------------------------------------------------------
units		metal
atom_style	molecular

variable fname string 'LuoArmchairSingleCrystal-a0-10-r02'

label equi_timesteps
variable equi_loop loop 1 100

    variable loop_val equal ${equi_loop}
    read_restart   ${fname}_${loop_val}.rst


    pair_style rebo #1.92
    pair_coeff * * CH_Luo.rebo C

    mass * 12.0107 #g/mol

    newton on

    #-----------------------------------------------------------
    # Set up computes
    #-----------------------------------------------------------

    compute temps all temp
    # Stress/ Atom Method
    compute SAperatom all stress/atom NULL #virial
    # average per atom. (Summing and dividing over total volume) and (averaging peratom stress then dividing by single volume) are the same
    compute SApress all reduce ave c_SAperatom[1] c_SAperatom[2] c_SAperatom[3] c_SAperatom[4] c_SAperatom[5] c_SAperatom[6]

    #-----------------------------------------------------------
    # Run Equilibration
    #-----------------------------------------------------------
    # fix	1 all npt temp ${T1} ${T1} ${tdamp} x 1.01325 1.01325 ${pdamp} z 0.0 0.0 ${pdamp}
    #
    # thermo_style custom step temp
    # thermo          1000
    # timestep ${tms}
    # run 100000

    #-----------------------------------------------------------
    # Run Tensile
    #-----------------------------------------------------------
    # Set up
    variable temps equal c_temps

    variable tmpx equal "lx"
    variable tmpy equal "ly"
    variable Lx0 equal ${tmpx}
    variable Ly0 equal ${tmpy}

    # Strain
    variable strainx equal "(lx - v_Lx0)/v_Lx0"
    variable strainy equal "(ly - v_Ly0)/v_Ly0"

    variable truestrainx equal "ln(lx/v_Lx0)"
    variable truestrainy equal "ln(ly/v_Ly0)"

    # graphene volume
    variable xyarea equal "lx*ly"
    variable graph_thickness equal 3.4 # taken from Luo Paper
    variable totalbeads equal count(all)
    variable graph_vol equal "v_xyarea*v_graph_thickness/v_totalbeads"
    # area ratio
    variable aratio equal "v_xyarea/(v_Lx0*v_Ly0)"

    # sampling rate
    variable sample_rate equal 10

    # stress/atom is units pressure*volume, so divide by volume, and convert bar to Pa
    variable SA_stressxx equal c_SApress[1]/v_graph_vol*1E5
    variable SA_stressyy equal c_SApress[2]/v_graph_vol*1E5
    variable SA_stresszz equal c_SApress[3]/v_graph_vol*1E5
    variable SA_stressxy equal c_SApress[4]/v_graph_vol*1E5
    variable SA_stressxz equal c_SApress[5]/v_graph_vol*1E5
    variable SA_stressyz equal c_SApress[6]/v_graph_vol*1E5

    # out
    fix def21 all print ${sample_rate} "${strainy} ${truestrainy} ${SA_stressxx} ${SA_stressyy} ${SA_stresszz} ${SA_stressxy} ${SA_stressxz} ${SA_stressyz}" file Graphene_NPT_REBO_Strain.SA${sample_rate}_${loop_val}.txt screen no
    fix def22 all print ${sample_rate} "${temps} ${strainy} ${aratio} ${SA_stressxx} ${SA_stressyy} ${SA_stresszz}" file Graphene_NPT_TEMP_REBO_Strain.SA${sample_rate}_${loop_val}.txt screen no

    dump            14 all custom 100 Graphene_NPT_REBO_Strain.dump id mol type x y z c_SAperatom[*]

    # fixes and run

    # fix	1 all nvt temp ${T1} ${T1} ${tdamp}
    fix	1 all npt temp ${T1} ${T1} ${tdamp} x 1.01325 1.01325 ${pdamp} z 0.0 0.0 ${pdamp}
    fix 2 all deform 1 y erate ${tms}

    thermo_style custom step temp v_strainy v_SA_stressyy v_aratio
    thermo          100
    timestep ${tms}
    run 40000

    clear
next equi_loop
jump SELF equi_timesteps

@tjbarrett I don’t know of any significant changes to the AIREBO code between the last two stable releases. All I have are a few notes:

  • there is no GPU accelerated version of this, so trying to use the gpu suffix is a waste. There is a small probability that this causes a problem, since there is a /gpu version of fix npt (which does not do any acceleration, but expects to be used with /gpu pair styles)
  • similarly, you should check if there are differences in behavior with using /omp versus not using it (and between 1 and more OpenMP threads as well). This could help to narrow down where the difference is coming from.
  • if you want apples to apples comparisons, you should not use binary restart files, but convert them to data files and then run with either LAMMPS version starting from the data file.
  • you may want to check how sensitive these calculations are to small changes in the initial data by adding a small random displacement with displace_atoms to your calculations

@akohlmey Thanks Dr. Kohlmeyer. I’ll test the OMP settings, change to data files, and update on any differences I find.

So from some preliminary tests, it appears the fluctuations found in the Sept 2021 build results are self inflicted, the results of the gpu inclusion in the command line. When run with just OMP, the results are similar to the June 2022 build, as shown below.

No GPU

I’ll run some other combinations to see if anything else comes up, but for now that seems to be the culprit