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 = C11*e+0.5*C111*e*e+… where e is the strain defined as e=(lz-108.45)/108.45 in this calculation.