Simulating Melting-Solidification Cycle

Dear LAMMPS Users,

I am trying to set up a dynamic simulation such that two regions are
thermostatted regions within a periodic cell. One of the regions is kept
above the melting point temperature and the other below the melting point
temperature such that when atoms are constantly pulled (using fix move)
across these regions, they undergo melting and solidification
respectively.

I understand that using a static region in LAMMPS means assigns the initial
group of atoms in the geometric space to a region and those atoms remain in
that region even they are out of the geometric space. Because of this, I
have decided to use dynamic groups so that atoms are periodically assigned
to the group. However, I did not observe a cycle of melting and
solidification as expected.

I have inserted my input file below:

# my silicon #change_box by 0.473% # minimize # box

units metal
dimension 3
boundary p p p
atom_style atomic
newton on

lattice diamond 5.4307 orient x 1 0 -1 orient y -1 2 -1 orient z 1 1 1
region simbox block 0 5 0 5 0 10
create_box 1 simbox
create_atoms 1 box

change_box all x scale 1.00485 y scale 1.00485 z volume remap

mass 1 28.0855

pair_style sw
pair_coeff * * si.sw Si

neighbor 2.0 bin
neigh_modify every 1 delay 0 check yes

#Define Region
region 1 block INF INF INF INF 0 5.0
region 2 block INF INF INF INF 5.0 10
group cold region 1
group hot region 2

#Define Smaller Region
region 3 block INF INF INF INF 2 3.25
region 4 block INF INF INF INF 3.25 6.75
region 5 block INF INF INF INF 6.75 8

compute Tcold all temp/region 3
compute Thot all temp/region 5

group coldS dynamic all region 3 every 1 # updating group every timestep
group middle dynamic all region 4 every 1
group hotL dynamic all region 5 every 1

#compute and conversions
compute x all property/atom x
compute k all ke/atom
compute cc all chunk/atom bin/1d z lower 9.3735584 units box
variable kb equal 8.617339e-5
variable temp atom c_k/(1.5*v_kb)

timestep 0.001

#outputs
dump id all custom 10000 move_02.d id x y z vx vy vz
dump_modify id sort id
log move_02.log
thermo_style custom step temp press ke pe density enthalpy lx ly lz

fix AVE all ave/chunk 500 1000 500000 cc c_k temp density/mass file
move_02.profile #
restart 10000 restart_move_02a restart_move_02b

#Initialization @ 1000K
minimize 1.0e-4 1.0e-6 100 1000
velocity all create 1000 987456 dist gaussian
fix one all npt temp 1000 1000 0.1 z 0 0 1
run 50000
unfix one

## raise temperature to 1670K
fix 1 all npt temp 1000 1670 0.1 z 0 0 1
thermo 5000
run 50000
unfix 1

##Equilibrate at 1670 K
fix 2 all npt temp 1670 1670 0.1 z 0 0 1
thermo 5000
run 100000
unfix 2

# Half liquid and Half solid
fix 3 cold npt temp 1670 1670 0.1 z 0 0 1
fix 4 hot npt temp 1670 3000 0.1 z 0 0 1
thermo 5000
run 50000
unfix 4

#equilibrate melt region
fix 5 hot npt temp 3000 3000 0.1 z 0 0 1
thermo 5000
run 100000
unfix 5

##quench hot region
fix 6 hot npt temp 3000 1700 0.1 z 0 0 1
thermo 5000
run 50000
unfix 6

##equilibrate hot side
fix 7 hot npt temp 1700 1700 0.1 z 0 0 1
thermo 5000
run 100000
unfix 3
unfix 7

## fix smaller region and constant pulling
fix 8 coldS langevin 1650 1650 1.00 698748 tally yes
fix 9 hotL langevin 1700 1700 1.00 698748 tally yes
fix pull all move linear NULL NULL -0.005 ##constant
pulling

run 3000000

  I would be glad to receive useful insights.

Regards,
Victor Fabiyi
Email: [email protected]

I forgot to mention, I am using January 2017 LAMMPS version

I doubt that anyone will make the effort to understand
a quite complicated input script and simulation protocol.

So I’ll just give generic advice. Start with as simple
a script as possible and verify it does what you
expect (thermo output, computed diagnostics, visualization).

Then add complexity one step at a time, and repeat
the verification.

The last thing I would add is the pulling via fix move.

Steve