Hi,
I used EAM potential to obtain second-order elastic constants (SOECs) C11 and third-order elastic constants (TOECs) C111 of copper crystalline. The value of C11 I got is 169 GPa, which agress well with experimental data . However, the value of C111 is -387 GPa, which is much more smaller than experimental value -1500 GPa (See Ref. PR. 161, 673 (1967)). I believe the wrong value has no relevance to the accuracy of EAM potential, since a previous analystic work based on EAM potential got correct value of C111 (-1271 GPa) (See Ref. PRB. 53, 14080 (1996)). The problem may come from my LAMMPS script. Can anyone help me to find the bug of the script? Thanks very much.
The LAMMPS script is as follows:
##########################################Begin######
#file name: in.C11-100
variable f index Cu-C11-100
log $f.log
variable freq equal 1
units metal
boundary p p p
atom_style atomic
lattice fcc 3.615
region box block 0 15 0 15 0 30 units lattice
create_box 1 box
create_atoms 1 box
pair_style eam/alloy #EAM potential, I am sure it is correct.
pair_coeff * * Cu01.eam.alloy Cu
thermo_style custom step pzz lz
dump 1 all xyz ${freq} $f.xyz
timestep 0.001
min_style cg
thermo ${freq}
fix 1 all deform 1 z delta 0.0 -5.4225 units box #push 5%
run 1 post no
unfix 1
label begin
variable t loop 100 #Loop 100 times
minimize 1.0e-8 1.0e-8 1000 10000
fix 1 all deform 1 z delta 0.0 108.45E-3 units box #0.1% each step
run 1 post no
unfix 1
next t
jump in.C11-100 begin
############################################END######
After calculation, the results is just like:
Step Pzz Lz
2 93233.746 103.0275
4 91131.496 103.13595
6 89046.053 103.2444
8 86976.593 103.35285
…
Then I fit the stress-strain relation as: Pzz = C11e+0.5C111ee+… where e is the strain defined as e=(lz-108.45)/108.45 in this calculation.