Hello,

I am trying to find the Young Modulus values of the titanium crystal structure at 0 Kelvin. For stress calculation I use Clausius virial theorem or BDT (atomic) stress approach, both give the same result. In addition, since I work around 0 Kelvin in the stress calculation, the main effect comes from the potential part and there is no stress effect from the kinetic part, and when I include the temperature effect, there is no difference in the results.

When I use nvt or nve ensemble in the analysis, the results are in full agreement with the results in the literature, while when I use npt ensemble, there is a 35% decrease in the stress value in the direction I apply load and therefore in the young modulus value. Stress value appears to be zero because I reset the stress value in other directions.

I relax until the temperature and pressure are balanced. I take the barostat damping value as 1000xdelta(t) and the thermostat damping value as 100xdelta(t). If I increase the damping values, the system stabilizes more slowly and I have to run for longer. Thus, I reset the residual stresses, if any, on the part before starting the deformation. When calculating stress and strain in the analysis, I take the time average, the results do not change when I change the statistical averages.

Although I have tried the analyzes for different statistical situations and titanium crystal structures, I get correct results with the nve or nvt ensemble, and quite different results with the npt ensemble.

My intention was to simulate testing actually done in vacuum or atmospheric pressure, but for weeks now I can’t find the reason. I’ve included the script I used below, I’d be very grateful if you could help.

Kind regards,

Emel

# Input file for uniaxial tensile loading of single crystal Titanium

#--------------- INITIALIZATION ---------------

units metal

dimension 3

boundary p p p

atom_style full

newton on

neighbor 0.1 bin

neigh_modify every 1 check yes

##--------------- ATOM DEFINITION ---------------

read_data Ti_alpha.data

##--------------- FORCE FIELDS ---------------

pair_style meam/spline

pair_coeff * * Ti.meam.spline Ti

#--------------- SETTINGS ---------------

timestep 0.001

variable ts equal 0.001

velocity all create 0.5 4928459 rot yes dist gaussian

#--------------- EQUILIBRATION & MINIMIZATION ---------------

dump 1 all atom 1 minimize.lammpstrj

minimize 1.0e-8 1.0e-8 10000 100000

undump 1

#-------------- RELAXATION ---------------

fix 1 all npt temp 0.5 0.5 0.1 iso 0 0 1 #or aniso

dump 1 all atom 500 relaxation_test_npt.lammpstrj

thermo 500

run 50000

undump 1

#--------------- COMPUTATION ---------------

compute 1 all stress/atom NULL

compute 2 all reduce sum c_1[1] c_1[2] c_1[3] c_1[4] c_1[5] c_1[6]

#--------------- VARIABLES ----------------

reset_timestep 0

variable Lx equal lx

variable tmp equal “lx”

variable L0 equal ${tmp}

variable Ly equal ly

variable Lz equal lz

variable Vol equal vol

variable up equal 2e-3

#--------------- DEFORMATION -----------------

unfix 1

reset_timestep 0

fix 1 all npt temp 0.5 0.5 0.1 y 0 0 1 z 0 0 1

fix 2 all ave/time 1 100 100 c_2[1] c_2[2] c_2[3]

fix per all ave/atom 1 100 100 c_1[1] c_1[2] c_1[3]

variable delta equal {up}*{tmp}

fix 3 all deform 1 x delta 0 ${delta} units box remap x

#--------------- Output strain and stress info to file for units metal

variable strain equal “(lx - v_L0)/v_L0”

variable STR equal “v_strain”

variable PXX equal f_2[1]/10000/v_Vol

variable PYY equal f_2[2]/10000/v_Vol

variable PZZ equal f_2[3]/10000/v_Vol

fix def1 all print 100 “{STR} {PXX} {PYY} {PZZ}” file strain_stress.def1.txt screen no

thermo 500

thermo_style custom step temp v_STR v_PXX v_PYY v_PZZ ke pe press lx ly lz vol

dump impactor all custom 500 dump.Ti.* id type x y z fx fy fz f_per[1] f_per[2] f_per[3]

dump 1 all atom 500 tensile_npt.lammpstrj

run 50000

# SIMULATION DONE

print “All done”