Problem in constructing Fe-C alloy system

Dear Lammps users

I am trying to build up bcc-structured Fe(100) doped with carbon atoms. I plan to place the carbon atoms in the octahedral sites. Please see scripts of the input code and log files below:

input code:

Introduction of defects into the lattice

include X.txt
include Y.txt

include Z.txt
include ID.txt

units metal
dimension 3
boundary p p s
atom_style atomic
neighbor 2.0 bin
neigh_modify delay 5
thermo 1000

----------------------------------------Define variables-----------------------------------------

variable Temp index 1 # temperature for indentation tests
timestep 0.001
variable myseed equal 482748 # the value seed for the velocity
variable atomrate equal 5000 # the rate in timesteps that the atoms are dumped
variable time_eq equal 20000 # time steps for the equilibrium process
variable time_indent equal 20000 # time steps for indentation

variable tdamp equal “v_tstep100" # damping time for thermal equilibrium
variable pdamp equal "v_tstep
1000”

variable a equal 2.866 # lattice parameter for the metal
variable size equal 40 # length of the simulation box
variable height equal 20 # height of the simulation box
variable P_in equal ${size}*$a/2 # position of indentation on the X axis
variable m equal 55.85 # atomic weight gram per mole
variable r equal 50 # Tip radius in angstrom
variable k equal 0.05 # loading rate
variable s equal 10 # stiffness of the virtual indenter

print “-------------------------------Temperature=${Temp} K---------”

-------------------------------------Atom definition--------------------------------

lattice bcc a region total block 0 {size} 0 {size} 0 {height}
region body block 0 {size} 0 {size} 3 {height} region bottom block 0 {size} 0 ${size} 0 2

create_box 2 total
create_atoms 1 region total

variable i loop 10

label loopa

region {ID} block {X} {X}+1 {Y} {Y}+1 {Z} {Z}+1 #delete_atoms region {ID}
create_atoms 2 single {X}+0.5 {Y}+0.5 ${Z}

next X
next Y
next Z
next ID
next i

jump Fe100_C_T.txt loopa

mass 1 $m
mass 2 12

group koerper region body
group boden region bottom

#-----------------------------------Force field--------------------------------

pair_style eam/fs
pair_coeff * * FeC.eam.fs Fe C

#pair_style tersoff/zbl
#pair_coeff * * FeC.tersoff.zbl Fe C

-------------------------------------settings--------------------------------

compute myTemp all temp
compute new koerper temp
compute disp all displace/atom
compute pot all pe/atom
compute ke all ke/atom
compute s1 all stress/atom virial

variable depth1 equal $r+a*{height}+0.5-(step-${time_eq})dt$k # loading part

#------------------------------------------Equilibrium------------------------------------

#velocity all create $x 825577 dist gaussian

minimize 1.0e-6 1.0e-12 10000 100000

In the log file, I found the minimization failed, and the energy goes to infinity:

#velocity all create $x 825577 dist gaussian

minimize 1.0e-6 1.0e-12 10000 100000
WARNING: Resetting reneighboring criteria during minimization (…/min.cpp:173)
Memory usage per processor = 4.21951 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
0 0 inf 0 inf -nan 753469
1 0 inf 0 inf -nan 753469
Loop time of 6.66577 on 45 procs for 1 steps with 65610 atoms

Minimization stats:
Stopping criterion = linesearch alpha is zero
Energy initial, next-to-last, final =
inf inf inf
Force two-norm initial, final = -nan -nan
Force max component initial, final = 0.832085 0.832085
Final line search alpha, max atom move = 0 0
Iterations, force evaluations = 1 1021

Pair time () = 4.78796 (71.8291) Neigh time () = 0 (0)
Comm time () = 1.75883 (26.386) Outpt time () = 0 (0)
Other time (%) = 0.118975 (1.78487)

Nlocal: 1458 ave 1514 max 1352 min
Histogram: 5 0 4 6 0 0 10 0 0 20
Nghost: 3618.69 ave 4414 max 3222 min
Histogram: 30 0 0 0 0 0 0 0 0 15
Neighs: 94537.1 ave 103422 max 85049 min
Histogram: 5 5 5 4 1 9 1 5 0 10

Total # of neighbors = 4254171
Ave neighs/atom = 64.8403
Neighbor list builds = 0
Dangerous builds = 0
velocity boden set 0.0 0.0 0.0 units box
velocity koerper create 300 ${myseed} temp new
velocity koerper create 300 482748 temp new

fix 1 boden setforce 0.0 0.0 0.0
fix 2 all temp/rescale 1000 {Temp} {Temp} 1.0 1.0
fix 2 all temp/rescale 1000 1 ${Temp} 1.0 1.0
fix 2 all temp/rescale 1000 1 1 1.0 1.0

fix 3 all npt temp 300 300 0.8 x 0.0 0.0 0.1 y 0 0 0.1 drag 1.0

dump 2 all custom 50000 initial_${Temp}K.lammpstrj id type x y z
dump 2 all custom 50000 initial_1K.lammpstrj id type x y z

run ${time_eq}
run 20000
Memory usage per processor = 3.53286 Mbytes
Step Temp E_pair E_mol TotEng Press Volume
1 256.09901 inf 0 inf -nan 753469
1000 -nan 0 0 -nan -nan -nan
2000 -nan 0 0 -nan -nan -nan

Could you please give me some clues about this problem? Thank you so much!

Best regards

Dong Wu

Graduate Research Assistant
Ph. D Student
Department of Material Science and Engineering
The University of Tennessee, Knoxville

log.lammps (17.2 KB)

Your script is long and complicated. You need
to debug it. Start with building one or 2 unit

cells. Dump the coords. Check the file to see

if you built the lattice structure you want. Repeat
for a larger system and visualize it. If you
do a run 0, you should get the same energy (per atom)
as for the a few unit cells. Then proceed.

Steve