How to solve the problem of energy and temperature deviation in the two-dimensional configuration of Kob-Anderson LJ

[jun_wu]

The above is the configuration details of the paper, and the following is the lammps script I wrote：
units lj
dimension 2
atom_style atomic

boundary p p p

timestep 0.005

region myreg block 0 40.8248 0 40.8248 -0.5 0.5
#egion myreg block 0 65.6062 0 65.6062 -0.5 0.5

create_box 2 myreg
create_atoms 1 random 1300 31235 myreg
create_atoms 2 random 700 53341 myreg

mass 1 1
mass 2 1

pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2

#min_style cg
minimize 1.0e-4 1.0e-4 10000 10000 #

#The call terminal (shell) creates a folder to store the computed data
shell mkdir ./configs_log_e
shell mkdir ./configs_quench_e
shell mkdir ./configs_equ_e
log log_e.lammps

#generates an ensemble of velocities
velocity all create 4 34343 dist gaussian
reset_timestep 0

#sets parameters for the building of pairwise neighbor lists
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

compute and thermo
compute peratom all pe/atom
compute pe all reduce sum c_peratom

compute 1 all temp/com
compute myTemp mobile temp/com
#compute_modify myTemp dynamic/dof yes

thermo 1000
thermo_style custom step atoms temp press pe c_pe ke etotal
thermo_modify flush yes

#--------melting at a high temperature---------
fix 1 all nvt temp 4 4 $(100*dt)
fix f2d1 all enforce2d
run 1000000
unfix 1
unfix f2d1

#-------------------quench---------------------
fix 2 all nvt temp 4 0.4 $(100*dt)
fix f2d2 all enforce2d
reset_timestep 0

dump 2 all custom 10000 ./configs_quench_e/heat.* id type xu yu x y
dump_modify 2 sort id

run 9000000
unfix 2
unfix f2d2
undump 2

#------------------equ----------------------
dump files at fixed intervals
variable f file merged_powers.txt
variable s equal “0+next(f)”

fix 3 all nvt temp 0.4 0.4 $(100*dt)
fix f2d3 all enforce2d
run 2000000

reset_timestep 0

dump 3 all custom 10000000 ./configs_log_e/heat.* id type xu yu x y
dump_modify 3 sort id
#dump_modify 3 format 3 .15lf format 4 .15lf #format 7 .15lf format 8 .15lf #Determine output accuracy
dump_modify 3 every v_s first yes
dump 4 all custom 10000 ./configs_equ_e/heat.* id type xu yu x y
dump_modify 4 sort id
run 5000000
unfix 3
undump 3
undump 4

write_data kalj_e.data nofix nocoeff
------------------------------end------------------------------
However, I found that the Potential energy per particle in the cooling process is inconsistent with the paper：

Did I ignore some commands in the process of writing the IN script, which caused the energy to deviate？I would be grateful if you could answer this question

[akohlmey]

The quoted paper says “shitfted force LJ” potential, but you are using a plain LJ with cutoff only.

[jun_wu]

Thank you very much for your answer. I thought it was because I ignored some commands. I will correct the error here. Again, thank you very much.

[akohlmey]

There could be other issues. I did not look any closer after I saw the mismatched pair style.

[jun_wu]

I’m sorry to bother you again. After receiving your reply, I modified my script as follows：
units lj
dimension 2
atom_style atomic

boundary p p p

timestep 0.005

region myreg block 0 40.8248 0 40.8248 -0.5 0.5
#egion myreg block 0 65.6062 0 65.6062 -0.5 0.5

create_box 2 myreg
create_atoms 1 random 1300 31235 myreg
create_atoms 2 random 700 53341 myreg

mass 1 1
mass 2 1

pair_style lj/smooth/linear 2.5
pair_coeff 1 1 1.0 1.0 2.5
pair_coeff 1 2 1.5 0.8 2.0
pair_coeff 2 2 0.5 0.88 2.2

#min_style cg
minimize 1.0e-4 1.0e-4 10000 10000 #

#The call terminal (shell) creates a folder to store the computed data
shell mkdir ./configs_log_e
shell mkdir ./configs_quench_e
shell mkdir ./configs_equ_e
log log_e.lammps

#generates an ensemble of velocities
velocity all create 4 34343 dist gaussian
reset_timestep 0

#sets parameters for the building of pairwise neighbor lists
neighbor 0.3 bin
neigh_modify every 1 delay 0 check yes

compute potential energy per particle and thermo
compute peratom all pe/atom
compute pe all reduce sum c_peratom

thermo 10000
thermo_style custom step atoms temp press pe c_pe ke etotal
thermo_modify flush yes

#--------melting at a high temperature---------
fix 1 all nvt temp 4 4 $(100*dt)
fix f2d1 all enforce2d
run 1000000
unfix 1
unfix f2d1

#-------------------quench---------------------
fix 2 all nvt temp 4 0.05 $(100*dt)
fix f2d2 all enforce2d
reset_timestep 0

dump 2 all custom 10000 ./configs_quench_e/heat.* id type xu yu x y
dump_modify 2 sort id

run 9500000
unfix 2
unfix f2d2
undump 2

#------------------equ----------------------
dump files at fixed intervals
variable f file merged_powers.txt
variable s equal “0+next(f)”

fix 3 all nvt temp 0.05 0.05 $(100*dt)
fix f2d3 all enforce2d
run 2000000

reset_timestep 0

dump 3 all custom 10000000 ./configs_log_e/heat.* id type xu yu x y
dump_modify 3 sort id
dump_modify 3 every v_s first yes
dump 4 all custom 10000 ./configs_equ_e/heat.* id type xu yu x y
dump_modify 4 sort id
run 5000000
unfix 3
undump 3
undump 4

write_data kalj_e.data
However, after recalculating the calculation, I found that the potential energy per particle in the inherent structure still could not correspond to the text. Theoretically, the curve should be flat during the cooling process, but it did not. Could you please take a look at my script again to see if there are any mistakes I made？Thank you very much.

[akohlmey]

Please look at the post. Can you read all of it correctly? I cannot. If you quote text, you have to do it properly to avoid the markup processing of the forum mess things up. This is explained in the guidelines post.

@jun_wu

Could you please take a look at my script again to see if there are any mistakes I made？

Sorry, but I am not your adviser or tutor and I am not familiar with your research and all details related to it. I cannot see any general bad choices or syntax errors. That is all that can be done from the LAMMPS perspective. Anything else is really more a question of the science and how to translate it to LAMMPS commands. That kind of question is beyond the scope of this forum.

[jun_wu]

Thank you for your reply. I will discuss with my tutor and check if there are any details wrong here

[srtee]

As a quick observation, in the paper screenshot, the average potential energy ranges between -3.34 and -3.42. In your results the range is much larger, which would seem far more typical of a less dense system. I can’t see any obvious issue with your input, but you should visualise your trajectories and see whether the particle density created is correct.