How to solve this energy minimization issue for high entropy alloys in lammps?

I am trying to determine some mechanical properties of Hf-Nb-Ta-Zr high entropy alloy. But I’m not being able to relax it properly & it’s causing absurd values of elastic moduli, double-necking during tensile tests etc. issues. Same codes work for single metals but I don’t know what to do for HEAs. The potential energy values do converge in thermo outputs… I can’t figure out what the exact issue is. The minimization part of the code is:

#minimization

fix rel all box/relax x 0 y 0 z 0 couple none

minimize 1e-6 1e-8 10000 100000

min_style cg

unfix rel

reset_timestep 0

#NVT_thermalization

velocity all create 1.0e-5 123456 mom yes rot yes dist gaussian

fix NVT all nvt temp 1.0e-5 ${STemp} 0.05 tchain 10

run ${EQRUN}

unfix NVT

#Annealing

fix equil all langevin ${STemp} 1e-5 0.05 123456 tally yes zero yes

fix NVE all nve

run ${EQRUN}

unfix equil

unfix NVE

#Pressure Relaxation

fix NPT all npt temp 1.0e-5 1.0e-5 100 x 0.0 0.0 1000 y 0.0 0.0 1000

run ${EQRUN}

unfix NPT

#NVE Equilbrium

fix NVE all nve

run ${EQRUN}

unfix NVE

You may be able to use the protocols in this paper: Full article: Molecular dynamics study on dynamic mechanical behaviour of FeCoCrCuNi high entropy alloy the material has high entropy (as the name says!) and so having a properly mixed initial condition will be crucial, by random generation on a lattice and/or Monte Carlo atom swaps.

1 Like

What makes you think that the issue is with the minimization protocol and not the system setup or the choice of potentials?

I have run the same code for for Al and Si; both yielded satisfactory results. Those two only required the direct minimization command. For the alloy, just the necessary create_box and create_atoms command were changed and it caused weird results even after nvt & npt. So the problem is either the potential or my relaxation procedure. Since the potential was procured from the reference of a journal paper, I think my code is the issue here

Yes, but it is more likely your code for creating the system is at fault for issues like you describe than the minimization protocol. This is supported by the fact, that you get expected results for single element systems.

Wow so that means the problem lies before minimization. I mistook the initialization as too simple to have a problem there, sorry. Can you please let me know where the issue is? This is the initial geometry creation upto minimization:

#-------variables-----------
variable STemp equal 300.0
variable EQRUN equal 50000
variable srate equal 1.0e10

------------------------ INITIALIZATION ----------------------------

units metal
dimension 3
boundary p p p
atom_style atomic
variable latparam equal 3.63
variable fHf equal 25
variable fNb equal 25
variable fTa equal 25
variable fZr equal 25

----------------------- ATOM DEFINITION ----------------------------

lattice bcc ${latparam}
region HEA block 0 100 0 20 0 20
create_box 4 HEA

------------------ Create different types of Atoms --------------------

create_atoms 1 region HEA
set type 1 type/fraction 2 ((v_fNb+v_fTa+v_fZr)/100.0) 1738 set type 2 type/fraction 3 ((v_fTa+v_fZr)/(v_fNb+v_fTa+v_fZr)) 1755
set type 3 type/fraction 4 $(v_fZr/(v_fTa+v_fZr)) 6688

------------------------ FORCE FIELDS ------------------------------

pair_style eam/alloy
pair_coeff * * Hf-Nb-Ta-Zr.eam.alloy Hf Nb Ta Zr
neighbor 2.0 bin
neigh_modify delay 10 check yes

------------------------- SETTINGS ---------------------------------

compute csym all centro/atom bcc
compute potatom all pe/atom

######################################

EQUILIBRATION

reset_timestep 0
timestep 0.001
#minimization

Please wrap your input script as a code block using triple backticks,

```
like this,
```

so we can read it without formatting problems. Thanks :slight_smile:

#-------variables-----------
variable STemp equal 300.0
variable EQRUN equal 50000
variable srate equal 1.0e10

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

units metal
dimension 3
boundary p p p
atom_style atomic
variable latparam equal 3.63
variable fHf equal 25
variable fNb equal 25
variable fTa equal 25
variable fZr equal 25

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

lattice bcc ${latparam}
region HEA block 0 100 0 20 0 20
create_box 4 HEA

# ------------------ Create different types of Atoms --------------------

create_atoms 1 region HEA
set type 1 type/fraction 2 ((v_fNb+v_fTa+v_fZr)/100.0) 1738 set type 2 type/fraction 3 ((vfNb+vfTa+vfZr)/100.0)1738settype2type/fraction3((v_fNb+v_fTa+v_fZr)/100.0) 1738 set type 2 type/fraction 3 ((v_fTa+v_fZr)/(v_fNb+v_fTa+v_fZr)) 1755
set type 3 type/fraction 4 $(v_fZr/(v_fTa+v_fZr)) 6688

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

pair_style eam/alloy
pair_coeff * * Hf-Nb-Ta-Zr.eam.alloy Hf Nb Ta Zr
neighbor 2.0 bin
neigh_modify delay 10 check yes

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

compute csym all centro/atom bcc
compute potatom all pe/atom

######################################

# EQUILIBRATION

reset_timestep 0
timestep 0.001
#minimization