# Problem in equilibrating the system

Dear Lammps users,

I am trying to simulate the indentation process. Because of the nature of the problem, I tried to fix the position of a layer of atoms in the bottom, and fix the temperature of the lateral layers of the box named thermostat layers.

For doing so I modeled a 300200300 box and divided it in in three different groups (base; fixed atoms, mobile; newton atoms and thermo; thermostatic atoms), as follows.

units metal
echo both
atom_style atomic

dimension 3
boundary p s p
region box block 0 300 0 200 0 300 units box
create_box 1 box
lattice fcc 4.08 orient x 1 1 1 orient y 1 -1 0 orient z 1 1 -2
region au block 0 300 0 300 0 300 units box
create_atoms 1 region au units box
timestep 0.0015
mass * 196.96655
pair_style eam/alloy
pair_coeff * * Au_GLJ10_3.eam.alloy Au

variable t equal 300 # Desired Simulation Temperature

# rigid boundry

region 1 block 0 300 0 10 0 300 units box
group base region 1
region 2 block 0 20 10 200 0 300 units box
group thermo2 region 2
region 3 block 280 300 10 200 0 300 units box
group thermo3 region 3
region 4 block 0 300 10 200 0 20 units box
group thermo4 region 4
region 5 block 0 300 10 200 280 300 units box
group thermo5 region 5
region 6 block 20 280 10 30 20 280 units box
group thermo6 region 6

group thermo union thermo2 thermo3 thermo4 thermo5 thermo6
group fixed union base thermo
group mobile subtract all base

Then I minimized the structure;

min_style cg
min_modify line backtrack
minimize 1.0e-6 1.0e-8 3000 100000
undump cg

After that I tried to equilibrate all the atoms in 300K; and remove the temperature fixture from mobile atoms and assign it to thermo atoms after equilibration;

# initial velocities and thermostat layer

compute new mobile temp
velocity mobile create 300 482748 temp new
fix 1 all nve
fix 2 base setforce 0.0 0.0 0.0
fix 3 all temp/rescale 10 300 300 1.0 0.5

thermo 400

thermo_style custom step ke pe etotal temp
thermo_modify lost warn flush yes
dump 1 all custom 2000 CSN.eqil id type x y z c_csym c_pnrg

run 30000

unfix 3
undump 1
fix 3 thermo temp/rescale 10 300 300 1.0 0.5

The problem is that after equilibration state and when I look at my sysstem, I see that there are some dislocation in the system after equilibration (as shown below), while there shouldn’t be any dislocation in the system if it’s well equilibrated.

I have also tried to this with nvt ensemble, but I got the same result;

compute new mobile temp
velocity mobile create 300 482748 temp new
fix 1 all nve
fix 2 base setforce 0.0 0.0 0.0
fix 3 all nvt temp 300.0 300.0 0.5

thermo 400

thermo_style custom step ke pe etotal temp
thermo_modify lost warn flush yes
dump 1 all custom 2000 CSN.eqil id type x y z c_csym c_pnrg

run 30000

unfix 3
undump 1
fix 3 thermo nvt temp 300.0 300.0 0.5

I also tried to with assigning temperature to the mobile group, but still I got the same result;

compute new mobile temp
velocity mobile create 300 482748 temp new
fix 1 all nve
fix 2 base setforce 0.0 0.0 0.0
fix 3 mobile nvt temp 300.0 300.0 0.5

thermo 400

thermo_style custom step ke pe etotal temp
thermo_modify lost warn flush yes
dump 1 all custom 2000 CSN.eqil id type x y z c_csym c_pnrg

run 30000

unfix 3
undump 1
fix 3 thermo nvt temp 300.0 300.0 0.5

After that I tried with ruining the system for 30000 steps after assigning fix to thermo layer, but it didn’t help.
It’s now a long time that I’m trying to solve this problem without a result. I would really appreciate if any one can help me to figure it out.

Thanks

Just quick assessment;

I would advise against the fixed atoms on the lateral sides, only fix the bottom layer opposite of the indenter.

Same goes for the thermostated layers, I suggest only use at the bottom. Or not at all, I don’t see the benefit.

Are you using the lattice constant at 300K? Otherwise the fixed atoms will cause a biaxial compression on the system which in this case I would assume is undesirable. You can minimize the entire system with minimize+box/relax and then fix the bottom later right before dynamics.

If you suspect the EAM force field is bad (do your own testing), you could equilibrate a bulk system (prefect crystal, no fixed regions) at the desired temperature/pressure.

-Mitch

Thanks a lot Mitch;

Yes, you are right. As the box is periodic in lateral directions, the lateral thermostat layers wouldn’t have a significant effect.

No, unfortunately I didn’t use the lattice constant at 300K. So, I’ll go with minimizing entire system and then fixing the bottom later as you mentioned.

The EAM force field is working well as I have previously equilibrated a bulk system to the 300K and there was no defect in the structure after equilibration.

Thanks a lot again,

Mahdi