LAMMPS code debugging

Hello! I am trying to run a uniaxial compression simulation for a High Entropy Alloy. The results are unexpected & I have found one possible reason from another thread. But I think there might be some syntax or other problems too. I’d be grateful if someone could help me check if the code is working as intended. Thanks

# ------------------------ INITIALIZATION ----------------------------
units 		metal
dimension	3
boundary 	p	p	p
atom_style	atomic
variable latparam equal 3.404
variable 	fHf equal 20
variable	fNb equal 20
variable 	fTa equal 20
variable 	fTi equal 20
variable 	fZr equal 20

# ----------------------- ATOM DEFINITION ----------------------------
lattice		bcc ${latparam}
region		HEA block 0 10 0 10 0 10
create_box	5 HEA

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

create_atoms 1 region HEA
set type 1 type/fraction 2 $((v_fNb+v_fTa+v_fTi+v_fZr)/100.0) 1238
set type 2 type/fraction 3 $((v_fTa+v_fTi+v_fZr)/(v_fNb+v_fTa+v_fTi+v_fZr)) 1750
set type 3 type/fraction 4 $((v_fTi+v_fZr)/(v_fTa+v_fTi+v_fZr)) 6648
set type 4 type/fraction 5 $((v_fZr)/(v_fTi+v_fZr)) 1441



# ------------------------ FORCE FIELDS ------------------------------
pair_style	meam
pair_coeff	* * library.meam Hf Nb Ta Ti Zr HfNbTaTiZr.meam Hf Nb Ta Ti 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

minimize 1e-4 1e- 100 1000

velocity all create 300 12345 mom yes rot no
fix 1 all npt temp 300 300 1 iso 0 0 1 drag 1 

# Set thermo output
thermo 500
thermo_style custom step lx ly lz press pxx pyy pzz pe temp

dump 		1 all cfg 250 dump.hnttz_T1w_velocity.comp_*.cfg mass type xs ys zs c_csym c_potatom fx fy fz

# Run for at least 10 picosecond (assuming 1 fs timestep)
run 20000
unfix 1

# Store final cell length for strain calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print " Initial Length, L0: ${L0} "

######################################
# DEFORMATION
reset_timestep	0

fix		1 all npt temp 300 300 1 y 0 0 1 z 0 0 1 drag 1
variable srate equal 1.0e10
variable srate1 equal "-v_srate / 1.0e12"
fix		2 all deform 1 x erate ${srate1} units box remap x

# Output strain and stress info to file
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# p2, p3, p4 are in GPa
variable strain equal " (lx - v_L0)/v_L0 "
variable p1 equal "v_strain"
variable p2 equal "-pxx/10000"
variable p3 equal "-pyy/10000"
variable p4 equal "-pzz/10000"
fix def1 all print 100 " ${p1} ${p2} ${p3} ${p4} " file hnttz_T1w_comp.def1.txt screen no

# Use cfg for AtomEye
dump 		2 all cfg 250 dump.hnttz_T1w_strain.comp_*.cfg mass type xs ys zs c_csym c_potatom fx fy fz

# Display thermo
thermo 	500
thermo_style	custom step v_strain temp v_p2 v_p3 v_p4 ke pe press

run		50000


######################################
# SIMULATION DONE

NB: Some words & symbols have been deformed due to direct copy-pasting. I would request to not bother with that & just help me with the logistic side. Thanks in advance

What is the error message you receive?

Btw, I doubt anybody will be able to help you without knowing what “the unexpected results are”.

Hi @Shrabon_B,

You cannot reasonably ask for debugging and mention that your copy-pasting altered the text. I see no reason for that, and if you are referring to the format, this is because the forum uses markdown for message formatting (you won’t regret taking time to learn a bit of it).

If you want to maintain the plain text/code formatting, encapsulate your pasted code with three back ticks (the ` character). This will make your input file way more readable.

But as @simongravelle mentioned, you also need to be more specific about your problem and the discrepancy between your output and the expected results.

3 Likes

@simongravelle @Germain Hey guys, sorry for being unclear. The problems are: i) The Stress-Strain curve doesn’t really look like regular Stress-Strain curves. It has two peaks at the initial stage. I haven’t let it run until fracture due to that. ii) Even though I specified it as BCC structure with lattice constant of 3.404A, the OVITO tool fails to detect the particles as BCC. Initially, 90% of the atoms are shown to have ‘Other’ structure type & just 10% as BCC. After 200 timesteps, almost 95% of those are ‘Other’ type, some are HCP, some are BCC & a few are FCC

When looking at your input, I dont see anything that seems wrong. I can’t say for sure, but it could be that the issue comes from the HfNbTaTiZr.meam file.

You are right to worry about the system not behaving as a BCC. In your place, I would try first to get the system behaving ‘well’ at equilibrium before launching such stress-strain measurement.

1 Like

I obtained the meam file from a fellow undergrad student through a middleman, so it isn’t completely reliable. In that case, what is the best way for getting this type of potential files which aren’t available by default on lammps? I haven’t found anything for this specific alloy on Github. Or do you recommend any other potential type instead of meam ? I just chose it because it was available

The system you are trying to simulate is rather complicated. The best place to look is always the literature. As you state that:

I obtained the meam file from a fellow undergrad student through a middleman

I cannot be sure but given the specificity of the elements and the potential form, I highly suspect that the files come from the provided material of this paper. You should start by being sure that this is the case, and if so go on with a proper reading to understand the expected behaviour of the potential with regard to the structure and mechanical properties and compare with what you obtained after each step of your simulation procedure (even equilibration as stated by Simon).

Keep in mind that the best person to ask for advices on your project might be a direct advisor or colleague.

edit: After quickly doing the maths from your strain parameters, it seems reasonable to me that you loose the structure at some point given the final expected x box length. Also consider the sign of your strain rate. As far as I can tell, there is no reason to consider any bugs here. Large strains will result in out of equilibrium structure.

3 Likes