%%writefile in.single_comp_diffusion

#---------------Initialize Simulation -----------------------#
units lj                # 单位
dimension 2             # 维度
boundary p p p          # 边界条件
atom_style atomic       # 原子类型

#-------------- Create Atoms Initial Conditions------------- #
lattice hex 1.0                        # 晶格类型
region box block 0 20 0 10 -0.1 0.1    # 空间区域
create_box 1 box                       # 模拟盒子
region 2 block 5 15 0 10 -0.1 0.1
create_atoms 1 region 2                # 基于晶格点阵创建原子
mass 1 1.0                             # 质量
velocity all create 2.5 87287          # 速度

#---------------- Define Interatomic Potential --------------#
pair_style lj/cut 2.5                  # 原子相互作用势
pair_coeff 1 1 1.0 1.0 2.5             # 对势参数
neighbor 0.3 bin                       # 近邻列表
neigh_modify every 20 delay 0 check no # 近邻算法
fix 1 all nvt temp 0.5 0.5 0.01        # 时间步操作
fix 2 all enforce2d

#--------------- Run MD Simulation --------------------------#
dump 1 all custom 100 toEquil.lammpstrj id type x y z vx vy vz    # 输出分子动力学模拟信息
thermo 500                                                        # 输出热力学信息间隔
run 10000                                                         # 运行指定步数分子动力学模拟

# !lmp -i in.single_comp_diffusion

%%writefile in.multi_comp_diffusion

#---------------Initialize Simulation -----------------------#
units lj
dimension 2
boundary p p p
atom_style atomic
variable t equal 0.5             # 定义系统温度

#-------------- Create Atoms Initial Conditions------------- #
lattice sq 1.0
region box block 0 100 0 100 -0.5 0.5
create_box 2 box
create_atoms 1 random 2500 12345 box
create_atoms 2 random 2500 54321 box        # 随机创建原子
mass 1 1.0
mass 2 1.0

#---------------- Define Interatomic Potential --------------#
pair_style hybrid lj/cut 2.5 soft 5.0
pair_coeff 1 1 lj/cut 1.0 1.0 2.5
pair_coeff 2 2 lj/cut 1.0 1.0 2.5
pair_coeff 1 2 soft 5.0

#--------------- Run MD Simulation --------------------------#
compute eng all pe/atom                                        # 计算每个原子势能
compute eatoms all reduce sum c_eng                            # 计算所有原子势能
thermo_style custom step temp epair etotal press c_eatoms      # 定义输出的热力学量
thermo 1000
dump id all atom 100 dump.lammpstrj
minimize 1e-4 1e-6 1000 10000                                  # 防止原子距离过近
velocity all create $t 87287
fix nvt all nvt temp $t $t 0.01
run 50000

%%capture
# !lmp -i in.multi_comp_diffusion

%%writefile in.Ar_single_point

units lj
boundary p p p
atom_style atomic
lattice fcc 1.00
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 1.0
pair_style lj/cut 8.0
pair_coeff 1 1 1.0 1.0 8.0
variable P equal pe
variable L equal (count(all)/1.00)^(1/3)
run 0
print “Lattice constant: $L"
print "Cohesive Energy of Ar: $P"

%%writefile in.Ar_loop

units lj
boundary p p p
atom_style atomic

label loop_i
variable i  loop  50
variable x  equal  1.02+0.002*$i

lattice fcc $x
region     box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 1.0

pair_style  lj/cut  4.0
pair_coeff  1  1  1.0  1.0  4.0

variable P equal pe
variable L equal (count(all)/$x)^(1/3)

run 0
print  "Cohesive Energy of Ar a = $L E = $P"
clear
next i
jump SELF loop_i

# Read the file and search for lines containing the pattern "Cohesive Energy of Ar a = * E = *
import re

# Initialize list to hold the matching lines
matching_lines = []

# Define the regex pattern for capturing a and E values
pattern = re.compile(r"Cohesive Energy of Ar a = ([\d\.-]+)\s+E = ([\d\.-]+)")

# Read the file and search for the pattern
with open('./log.lammps', 'r') as file:
    for line in file:
        match = pattern.search(line)
        if match:
            a_value = float(match.group(1))
            E_value = float(match.group(2))
            matching_lines.append((a_value, E_value))

# Write the extracted a and E values to a CSV file
with open('Ar_fcc.csv', 'w') as f:
    f.write("a,E\n")  # Write header
    for a, E in matching_lines:
        f.write(f"{a},{E}\n")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


filename = "Ar_fcc.csv"
DataFrame = pd.read_csv(filename)
lat = np.array(DataFrame['a'])
ene = np.array(DataFrame['E'])

a, b, c = np.polyfit(lat, ene, 2)
mesh = np.linspace(np.min(lat), np.max(lat), 1000)
opt_lat = -b / (2 * a)
Emin = np.polyval([a, b, c], opt_lat)
print("a,b,c =", a, b, c)
print("Emin =", Emin, "opt_lat =", opt_lat)

plt.scatter(lat, ene, 10)
plt.plot(mesh, a * mesh ** 2 + b * mesh + c)
plt.xlabel('Lattice parameter')
plt.ylabel('Energy')
plt.show()

%%writefile in.Al_loop

units metal
boundary p p p
atom_style atomic
variable i loop 30
variable x equal 3.90+0.01*$i

# Define a lattice for use by other commands.
lattice fcc $x
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box

pair_style eam/fs
pair_coeff * * Al_mm.eam.fs Al

variable n equal count(all)
variable P equal pe/$n

run 0
print "Cohesive Energy of Al a = $x E = $P"
clear 
next i
jump SELF

%%writefile in.Al_relax

units metal
dimension 3
boundary p p p
atom_style atomic

lattice fcc 4.00
region box block 0 1 0 1 0 1 units lattice
create_box 1 box
create_atoms 1 box

pair_style eam/fs
pair_coeff * * Al_mm.eam.fs Al

neighbor 2.0 bin
neigh_modify delay 10 check yes

compute eng all pe/atom
compute eatoms all reduce sum c_eng
dump 1 all atom 1 relax.lammpstrj

#--------------Run minimization--------------#
reset_timestep 0
fix 1 all box/relax iso 0.0 vmax 0.001
thermo 10
thermo_style custom step pe lx ly lz press pxx pyy pzz c_eatoms

min_style cg
minimize 1e-25 1e-25 5000 10000

variable natoms equal count(all)
variable teng equal "c_eatoms"
variable a equal lx
variable ecoh equal "v_teng/v_natoms“

print "Total energy (eV)= ${teng};" 
print "Number of atoms = ${natoms};"
print "Lattice constant (Angstroms) = ${a};" 
print "Cohesive energy (eV/atom) = ${ecoh};" 

%%writefile in.fcc_bulk_modulus

units lj
boundary p p p
atom_style atomic
label LOOP

# The volume change is controlled by lattice constant
variable i loop 40
variable x equal 1.05+0.002*$i

lattice fcc $x
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box

pair_style lj/cut 8.50
pair_coeff * * 1.0 1.0
mass 1 1.0

thermo_style custom etotal
variable P equal pe                  # Potential energy per atom
variable v equal (1.0/$x)            # vol in LJ unit is N/x, vol per atom in LJ units is 1/x

dump 1 all custom 10 trj.lammpstrj id type x y z vx vy vz
run 0
print "FCC lattice rho = $v E = $P"
clear
next i
jump SELF LOOP

# Read the file and search for lines containing the pattern "FCC lattice rho = * E = *"
import re

# Initialize list to hold the matching lines
matching_lines = []

# Define the regex pattern for capturing a and E values
pattern = re.compile(r"FCC lattice rho = ([\d\.-]+)\s+E = ([\d\.-]+)")

# Read the file and search for the pattern
with open('./log.lammps', 'r') as file:
    for line in file:
        match = pattern.search(line)
        if match:
            v_value = float(match.group(1))
            E_value = float(match.group(2))
            matching_lines.append((v_value, E_value))

# Write the extracted a and E values to a CSV file
with open('fcc_bulk_modulus.csv', 'w') as f:
    f.write("V,E\n")  # Write header
    for a, E in matching_lines:
        f.write(f"{V},{E}\n")

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt


filename = "fcc_bulk_modulus.csv"
DataFrame = pd.read_csv(filename)
vol = np.array(DataFrame['V'])
ene = np.array(DataFrame['E'])

a, b, c = np.polyfit(vol, ene, 2)
mesh = np.linspace(np.min(vol), np.max(vol), 1000)
opt_vol = -b / (2 * a)
Emin = np.polyval([a, b, c], opt_vol)
print("a,b,c =", a, b, c)
print("Emin =", Emin, "opt_vol =", opt_vol)

plt.scatter(vol, ene, 10)
plt.plot(mesh, a * mesh ** 2 + b * mesh + c)
plt.xlabel('Volume')
plt.ylabel('Energy')

plt.show()

import pandas as pd
import numpy as np


df = pd.read_csv('fcc_bulk_modulus.csv)
v = np.array(df['V'])
e = np.array(df['E'])

a, b, c = np.polyfit(v, e, 2)
v0 = -b / (2 * a)
e0 = a * v0 ** 2 + b * v0 + c
b0 = 2 * a * v0
bP = 3.5
x0 = [e0, b0, bP, v0]


# Birch–Murnaghan equation of state
def Murnaghan(parameters, vol):
    E0 = parameters[0]
    B0 = parameters[1]
    BP = parameters[2]
    V0 = parameters[3]
    E = E0 + B0 * vol / BP * (((V0 / vol) ** BP) / (BP - 1) + 1) - V0 * B0 / (BP - 1.)

    return E


# error
def residual(pars, y, x):
    err = y - Murnaghan(pars, x)

    return err


from scipy.optimize import leastsq

murnpars, ier = leastsq(residual, x0, args=(e, v))

print('Bulk Modulus:' + str(murnpars[1]))
print('lattice constant:', murnpars[3] ** (1 / 3))

from matplotlib import pyplot as plt

v_mesh = np.linspace(np.min(v), np.max(v), 1000)
plt.scatter(v, e, 10)
plt.plot(v_mesh, Murnaghan(murnpars, v_mesh))
plt.xlabel('Volume')
plt.ylabel('Energy')

plt.show()

%%writefile in.Cu_lj_mudulus

units metal
boundary p p p
atom_style atomic
#------------- setup loop -------------#
variable i loop 40
variable x equal 3.40+0.01*$i
lattice fcc $x
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 64

#------------- Cu LJ parameter -------------#
pair_style lj/cut 10.0
pair_coeff 1 1 0.40933 2.338
variable v equal ($x)^3
variable n equal count(all)
variable P equal pe

#------------- run -------------#
run 0
print "Cohesive Energy of Cu v = $v x = $x E = $P "
clear 
next i
jump SELF

%%writefile in.Cu_eam_mudulus

units metal
boundary p p p
atom_style atomic
#------------- setup loop -------------#
variable i loop 40
variable x equal 3.40+0.01*$i
lattice fcc $x
region box block 0 1 0 1 0 1
create_box 1 box
create_atoms 1 box
mass 1 64

#------------- Cu EAM parameter -------------#
pair_style eam
pair_coeff * * Cu_u3.eam
variable v equal ($x)^3
variable n equal count(all)
variable P equal pe

#------------- run -------------#
run 0
print "Cohesive Energy of Cu v = $v x = $x E = $P "
clear 
next i
jump SELF

%%writefile in.graphene

# Parameters
variable        Lx      equal   20
variable        Ly      equal   10
# variable      Lz      equal   0.8

# Unit Size
variable        x         equal    ${Lx}
variable        y         equal    ${Ly}
# variable        z         equal    ${Lz}
variable        xbox         equal    round(v_x)
variable        ybox         equal    round(v_y)
# variable        zbox         equal    round(v_z)

# Initialization
units                real
dimension            3
boundary             p p p
neighbor             2.0 bin
neigh_modify         every 10 delay 0 check yes
timestep             0.001
atom_style           charge

# Modeling
lattice custom  1.421 a1  3  0  0   a2   0 1.732  0 a3 0 0 2.357   &
        basis   0        0  0  &
        basis   0.333    0  0  &
        basis   0.5     0.5 0  &
        basis   0.833   0.5 0
region                box block 0 ${xbox} 0 ${ybox} -5.0 5.0
create_box            1 box
region                graphene block 0 ${xbox} 0 ${ybox} -0.1 0.1
create_atoms          1 region graphene
mass                  * 12.011150
write_data            graphene.dat

%%writefile in.tensile

# UNITS
units real
timestep 1
variable fpress equal 0.000101325        # atm -> GPa
variable fenergy equal 0.043             # kcal/mole -> eV

# GRAPHENE SHEET
dimension 3
boundary p p f
atom_style charge
read_data graphene.dat
pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.reax.cho C

# TEMPERATURE SETTINGS
variable temp equal 300
variable seed equal 1717
velocity all create ${temp} ${seed} dist uniform

# EQUILIBRATION
fix fnpt all npt temp ${temp} ${temp} 10 x 0 0 500 y 0 0 500
thermo 100
run 1000

# OUTPUT
# assume that the thickness of monolayer graphene is 3.35A.
variable tmp equal lx
variable lx0 equal ${tmp}
variable tmp equal ly
variable ly0 equal ${tmp}
variable Eavg equal etotal/atoms*${fenergy}                # eV/atom
variable pe equal pe/atoms*${fenergy}                      # eV/atom
variable ke equal ke/atoms*${fenergy}                      # eV/atom
variable strainx equal (lx-${lx0})/${lx0}
variable strainy equal (ly-${ly0})/${ly0}
variable stressx equal -pxx*(lz/3.35)*${fpress}            # GPa
variable stressy equal -pyy*(lz/3.35)*${fpress}            # GPa
thermo_style custom step time temp etotal press v_Eavg v_pe v_ke v_strainx v_stressx v_strainy v_stressy

# DEFORMATION
fix boxdeform all deform 1 x scale 2 remap x
fix fnpt all npt temp ${temp} ${temp} 10 y 0 0 500
fix output all ave/time 1 100 100 v_strainx v_stressx v_strainy v_stressy file stress_strain.txt
dump 1 all atom 10 dump.lammpstrj
fix stop all halt 100 v_strainx > 0.7  error  continue
thermo 100
run 10000

%%writefile in.thermal_expansion_coefficent

variable T index 300 400 500 600 700 800 900 1000
label T_loop

units metal
boundary p p p
atom_style atomic

# Create lattice
lattice fcc 3.62
region box block 0 8 0 8 0 8
create_box 1 box
create_atoms 1 box

# Set interatomic potential
pair_style eam
pair_coeff 1 1 Cu_u3.eam

# Reset timestep
reset_timestep 0

# Initialize velocity
velocity all create ${T} 87287 dist gaussian

# Equilibrate using NPT ensemble
fix 1 all npt temp ${T} ${T} $(100*dt) iso 0 0 1

# Define thermo output
thermo_style custom step temp epair press lx ly lz
thermo 1000

# Define temperature, length and volume computation
compute actual_T all temp
variable Lx equal lx
variable V equal vol

# Average properties over time and write to a single file
fix 2 all ave/time 100 10 10000 c_actual_T v_Lx v_V file thermal_expansion_data.${T}.txt

# Run simulation
run 10000

# Unfix the NPT ensemble
unfix 1

clear
next T
jump SELF T_loop

%%writefile in.Cv

units metal 
boundary p p p
atom_style atomic

variable Tstart equal 10.0
variable Tstop equal 1000
variable Tdamp equal 0.2

lattice fcc 3.62
region simbox block 0 8 0 8 0 8
create_box 1 simbox
create_atoms 1 box

pair_style eam
pair_coeff 1 1 Cu_u3.eam

velocity all create ${Tstart} 825577 dist gaussian

# Equilibrium step
fix 1 all nvt temp ${Tstart} ${Tstart} ${Tdamp}
thermo 100
run 5000
unfix 1
reset_timestep 0
thermo 1000
thermo_style custom step temp pe etotal

# Heating step
print "Production steps started!"
fix 2 all nvt temp ${Tstart} ${Tstop} ${Tdamp}
run 120000 
print "Simulation complete!"

%%writefile in.fusion

units lj
boundary p p p
atom_style atomic

lattice fcc 1.073
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0

pair_style lj/cut 5.0
pair_coeff 1 1 1.0 1.0 5.0
velocity all create 0.01 87287

timestep 0.005

# Equilibrium step
fix 1 all nvt temp 0.01 0.01 0.2
run 50000
unfix 1

# Heating step
reset_timestep 0
thermo 1000
thermo_style custom step temp pe etotal press vol
fix 1 all npt temp 0.01 0.85 1.0 iso 0 0 1.0
dump 1 all atom 10000 fusion.lammpstrj
run 1000000

%%writefile in.quench

units lj
boundary p p p
atom_style atomic

lattice fcc 1.073
region box block 0 10 0 10 0 10
create_box 1 box
create_atoms 1 box
mass 1 1.0

pair_style lj/cut 5.0
pair_coeff 1 1 1.0 1.0 5.0
velocity all create 0.01 87287

timestep 0.005

# Equilibrium step
fix 1 all nvt temp 0.85 0.85 0.2
run 50000
unfix 1

# Heating step
reset_timestep 0
thermo 1000
thermo_style custom step temp pe etotal press vol
fix 1 all npt temp 0.85 0.01 1.0 iso 0 0 1.0
dump 1 all atom 10000 quench.lammpstrj
run 1000000

%%writefile in.Cu_fusion

units metal
boundary p p p
atom_style atomic

lattice fcc 3.62
region box block 0 8 0 8 0 8
create_box 1 box
create_atoms 1 box

pair_style eam
pair_coeff 1 1 Cu_u3.eam

thermo 1000
thermo_style custom step temp pe etotal press vol

velocity all create 100.0 4928459 dist gaussian
fix 1 all npt temp 100.0 100.0 0.1 iso 0.0 0.0 1.0
run 10000
unfix 1

reset_timestep 0
compute 1 all msd
fix msd all ave/time 10 100 1000 c_1[1] c_1[2] c_1[3] c_1[4] file msd.dat     # Output msd
fix 1 all npt temp 100 2000 0.1 iso 0 0 1.0
dump 1 all atom 1000 dump.lammpstrj
run 120000

from ase import io
import numpy as np


def lindemann(filename, n_start=None, n_end=None):
    """
    Calculate the Lindemann index for atomic vibrations using the provided trajectory.

    Parameters:
    - filename: Name of the file containing atomic positions (typically a LAMMPS dump file).
    - n_start: Starting index for the trajectory.
    - n_end: Ending index for the trajectory.

    Returns:
    - Lindemann index value.
    """

    # Construct the slice string based on provided indices
    if n_start is None and n_end is None:
        slice_str = ":"
    elif n_start is None:
        slice_str = f":{n_end}"
    elif n_end is None:
        slice_str = f"{n_start}:"
    else:
        slice_str = f"{n_start}:{n_end}"

    # Read atomic positions from the file for the given range
    atoms = io.read(filename, slice_str, 'lammps-dump-text')

    # Total number of atoms and total number of timesteps
    N = atoms[0].positions.shape[0]
    N_t = len(atoms)
    print('number of atoms in each frame =', N, 'number of frames =', N_t)

    # Initialize arrays to store distances between atoms
    r_square = np.zeros((N_t, N, N))

    # Loop through each timestep
    for t in range(N_t):
        G = np.dot(atoms[t].positions, atoms[t].positions.T)
        H = np.tile(np.diag(G), (N, 1))

        # Compute the squared distance between each pair of atoms at time t
        r_square[t, :, :] = H + H.T - 2 * G

    # Compute average distance and squared average distance over all timesteps
    r_ave = np.mean(np.sqrt(r_square), axis=0)
    print(r_square.shape)
    print(r_ave.shape)
    print('r_ave_max =', np.max(r_ave), 'r_ave_min =', np.min(r_ave))

    r2_ave = np.mean(r_square, axis=0)
    print('r2_ave =', np.max(r2_ave))

    # Extract upper triangular part of matrices
    r_ave_upper = r_ave[np.triu_indices(N, k=1)]
    print('r_ave_upper_max =', np.max(r_ave_upper), 'r_ave_upper_min =', np.min(r_ave_upper))
    r2_ave_upper = r2_ave[np.triu_indices(N, k=1)]

    # Calculate Lindemann criterion for upper triangular elements
    value_upper = np.sqrt(r2_ave_upper - np.power(r_ave_upper, 2)) / r_ave_upper

    # Sum over all unique atom pairs
    total_value = np.sum(value_upper)

    # Normalize by the total number of atom pairs
    return total_value * 2 / (N * (N - 1))

%%writefile in.nanoparticles_fusion

# Initialization
units metal
boundary p p p
atom_style atomic

# Generate NiCu nanoparticles
variable A0 equal 3.589
lattice fcc ${A0}
region mybox block 0 40 0 40 0 40
create_box 2 mybox                         # 2 types of atoms
region CuNi_nano sphere 20 20 20 2
create_atoms 1 region CuNi_nano
set type 1 type/fraction 2 0.5 7777
mass 1 63.54600000                         # Cu
mass 2 58.69340000                         # Ni

# Using EAM potential
pair_style eam/alloy
pair_coeff * * CuNi.eam.alloy Cu Ni
write_data nanoparticle.cif
run 0

# Thermal equilibrium step
thermo 1000
variable j loop 0 20
label loop_j
variable temperature equal 900+10*$j
variable T equal temp
variable Eatom equal etotal/atoms
fix 1 all nvt temp ${temperature} ${temperature} 0.1
run 50000
unfix 1

# Statistical step
fix 2 all ave/time 100 5 1000 v_T v_Eatom file data_ave${temperature}.txt
dump 1 all atom 5000 fusion_${temperature}.atom
fix 1 all nvt temp ${temperature} ${temperature} 0.1
run 500000

unfix 1
undump 1
next j
jump SELF loop_j

%%writefile in.nanoparticles_sintering

# Initialization
units metal
boundary p p p
atom_style atomic
timestep 0.001

# Get three Pt nanoparticles 
variable A0 equal 3.9239
lattice fcc ${A0}
region mybox block 0 20 0 20 0 20
region sphere_pt1 sphere 7 10 8 3
region sphere_pt2 sphere 10 10 13.4 3
region sphere_pt3 sphere 13 10 8 3
create_box 1 mybox
create_atoms 1 region sphere_pt1
create_atoms 1 region sphere_pt2
create_atoms 1 region sphere_pt3

# Using EAM potential
pair_style eam
pair_coeff 1 1 Pt_u3.eam

# Output initial structure
dump 1 all cfg 1 coord.*.cfg mass type xs ys zs
run 0
undump 1

# Define triple_neck region
region triple_neck block 7.5 12.5 5 15 9 13
group 1 dynamic all region triple_neck every 1000

# Output number of atoms in triple neck region
variable N equal step
variable T equal temp
variable Natom equal count(all)
variable V equal vol/v_Natom
variable sinter_atom equal count(1)
dump 1 all xyz 1000 melt.xyz
thermo 1000
fix extra all print 1000 "${N} ${T} ${sinter_atom}" file data.txt

# Run in 500k, 1000k, 1400k, respectively
fix 1 all npt temp 500 500 0.1 iso 1 1 1
run 10000
unfix 1
fix 1 all npt temp 1000 1000 0.1 iso 1 1 1
run 20000
unfix 1
fix 1 all npt temp 1400 1400 0.1 iso 1 1 1
run 30000
unfix 1

%%writefile in.thermal_expansion

units real
boundary p p p
atom_style charge

read_data Li2S.dat

pair_style reax/c NULL checkqeq no
pair_coeff * * ffield.reax.LiS S Li

# Fix 1 all qeq/reax 1 0.0 10.0 1e-6 reax/c
neighbor 3.0 bin
min_style cg
minimize 1e-4 1e-6 100 1000

# Stablize at room Temp
velocity all create 300 96588
thermo 1000
thermo_style custom step temp press density vol
compute myRDF all rdf 50 1 1 1 2 2 1 2 2
fix 2 all ave/time 100 10 10000 c_myRDF[*] file initial.rdf mode vector
fix 3 all nvt temp 300 300 10
run 10000
unfix 2
unfix 3

# Heating slowly
fix 4 all npt temp 300.0 1800.0 100 iso 1.0 1.0 1000 drag 1.0
run 50000
unfix 4

# Keep melting
fix 5 all npt temp 1800.0 1800.0 100 iso 1.0 1.0 1000 drag 1.0
fix 6 all ave/time 100 10 5000 c_myRDF[*] file melted.rdf mode vector
run 5000
unfix 5
unfix 6

# Rapidly quenching
fix 7 all npt temp 1800.0 300.0 1 iso 1.0 1.0 1000 drag 0.1
thermo 10
run 100
unfix 7

# Stablize the amorphous LixS at RT
fix 8 all ave/time 100 10 5000 c_myRDF[*] file amorphous.rdf mode vector
fix 9 all npt temp 300 300 100 iso 1.0 1.0 1000 drag 1.0
thermo 1000

%%writefile in.hydrogen_storage

units metal
dimension 3
atom_style atomic
boundary p p p
read_data MgH.data
neighbor 0.5 bin
neigh_modify every 1 delay 0 check yes

region Mg block INF INF INF INF -0.01 5 units box
group Mg region Mg
velocity Mg set 0 0 0
fix 01 Mg setforce 0 0 0
region H2 block INF INF INF INF 6 INF units box
group H2 region H2

mass 1 24.3050
mass 2 1

pair_style adp
pair_coeff * * Mg_H.adp.alloy.txt Mg H

region 1H block INF INF INF 9 INF units box
group mobile1 dynamic all region 1H every 100
variable number1 equal count(mobile1)

region 2H block INF INF INF 5 9 units box
group mobile2 dynamic all region 2H every 100
variable number2 equal count(mobile2)

region 3H block INF INF INF -0.001 5 units box
group mobile3 dynamic all region 3H every 100
variable number3 equal count(mobile3)

# compute 1 Mg temp
compute 2 H2 temp
timestep 0.001

variable T equal TEMP
thermo 1000
thermo_style custom step temp pe etotal press vol c_2

min_style cg
minimize 1e-5 1e-5 10000 10000
run 0

velocity H2 create ${T} 71717
dump 1 all atom 100 dump.xyz
fix 1 H2 nvt temp ${T} ${T} 0.1
thermo_modify lost ignore
fix extra all print 100 "${N} ${P} ${number1} ${number2} ${number3}" file energy.data
run 200000

%%writefile in.langevin

units lj
atom_style atomic
lattice fcc 0.6
region box block 0 10 0 10 0 20
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 1.35 71717
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0
neighbor 0.3 bin

region hot block INF INF INF INF 0 1
region cold block INF INF INF INF 10 11

compute Thot all temp/region hot
compute Tcold all temp/region cold

# 1st Equilibration run
fix 1 all nvt temp 1.35 1.35 0.5
thermo 100
run 1000

velocity all scale 1.35
unfix 1

# 2nd Equilibration run
fix 1 all nve
fix hot all heat 1 100.0 region hot
fix cold all heat 1 -100.0 region cold
thermo_style custom step temp c_Thot c_Tcold
thermo 1000
run 10000

compute ke all ke/atom
variable temp atom c_ke/1.5

compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.heat
variable tdiff equal f_2[11][3]-f_2[1][3]

fix ave all ave/time 1 1 1000 v_tdiff ave running start 13000
thermo_style custom step temp c_Thot c_Tcold v_tdiff f_ave
run 20000

%%writefile in.mp

units lj
atom_style atomic
lattice fcc 0.6
region box block 0 10 0 10 0 20
create_box 1 box
create_atoms 1 box
mass 1 1.0
velocity all create 1.35 71717
pair_style lj/cut 2.5
pair_coeff 1 1 1.0 1.0
neighbor 0.3 bin

# 1st Equilibration run
fix 1 all nvt temp 1.35 1.35 0.5
thermo 100
run 1000
velocity all scale 1.35
unfix 1

# 2nd Equilibration run
compute ke all ke/atom
variable temp atom c_ke/1.5
fix 1 all nve
compute layers all chunk/atom bin/1d z lower 0.05 units reduced
fix 2 all ave/chunk 10 100 1000 layers v_temp file profile.mp
fix 3 all thermal/conductivity 10 z 20
variable tdiff equal f_2[11][3]-f_2[1][3]
thermo_style custom step temp epair etotal f_3 v_tdiff
thermo 1000
run 20000

# Thermal conductivity calculation
fix 3 all thermal/conductivity 10 z 20
fix ave all ave/time 1 1 1000 v_tdiff ave running
thermo_style custom step temp epair etotal f_3 v_tdiff f_ave
run 20000