# This is the control script for LAMMPS. Run LAMMPS on GPU
echo both
log 1.1_Initialize.out
#-------------------------------------------------------------------------------
# Stage 1.1: Initialize LAMMPS run for 3-d periodic
#-------------------------------------------------------------------------------
units real
boundary p p p
atom_style full
atom_modify map array
pair_style lj/class2/coul/long 9.5
pair_modify mix sixthpower
pair_modify tail yes
bond_style class2
angle_style class2
dihedral_style class2
improper_style class2
special_bonds lj 0.0 0.0 1.0 coul 0.0 0.0 1.0
box tilt large
read_data structure.dat
include parameters.dat
neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes
kspace_style pppm 0.00001
variable R equal 0.00198722
variable sysvol equal vol
variable sysmass equal mass(all)/6.0221367e+23
variable sysdensity equal v_sysmass/v_sysvol/1.0e-24
variable coulomb equal ecoul+elong
variable etotal equal etotal
variable pe equal pe
variable ke equal ke
variable evdwl equal evdwl
variable epair equal epair
variable ebond equal ebond
variable eangle equal eangle
variable edihed equal edihed
variable eimp equal eimp
variable lx equal lx
variable ly equal ly
variable lz equal lz
variable Nthermo equal 0
variable cella equal lx
variable cellb equal sqrt(ly*ly+xy*xy)
variable cellc equal sqrt(lz*lz+xz*xz+yz*yz)
variable cellalpha equal acos((xy*xz+ly*yz)/(v_cellb*v_cellc))
variable cellbeta equal acos(xz/v_cellc)
variable cellgamma equal acos(xy/v_cellb)
variable p equal press
variable pxx equal pxx
variable pyy equal pyy
variable pzz equal pzz
variable pyz equal pyz
variable pxz equal pxz
variable pxy equal pxy
variable sxx equal -pxx
variable syy equal -pyy
variable szz equal -pzz
variable syz equal -pyz
variable sxz equal -pxz
variable sxy equal -pxy
variable fmax equal fmax
variable fnorm equal fnorm
variable time equal step*dt+0.000001
variable surfacetension equal 0.5*v_lz*(0.5*(v_sxx+v_syy)-v_szz)
thermo_style custom step v_time press vol v_sysdensity temp ebond eangle edihed eimp evdwl ecoul etail elong pe ke
thermo_modify flush yes
#
# Set up the fixed and movable groups
#
group movable union all
group fixed subtract all movable
compute sum_f1 movable reduce sum fx fy fz
variable sum_fx equal c_sum_f1[1]
variable sum_fy equal c_sum_f1[2]
variable sum_fz equal c_sum_f1[3]
log 1.2_NVT.out
#-------------------------------------------------------------------------------
# Stage 1.2: NVT integration for 10 ps with a timestep of 1 fs
# Initial temperature 30 °C, final 300 °C
#-------------------------------------------------------------------------------
reset_timestep 0
thermo_style custom step v_time press vol v_sysdensity temp ebond eangle edihed eimp evdwl ecoul etail elong pe ke
thermo ${Nthermo}
fix 1 movable nvt temp 303.15 573.15 100 drag 0.0
fix 2 movable ave/time 1 1000 1000 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 1.2_averages.txt off 1
fix 3 movable ave/time 1 1 1 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 1.2_instantaneous.txt
restart 10000 1.2.restart
dump sci all custom 1000 1.2.xyz id mol type q xs ys zs
timestep 1
run 10000
undump sci
restart 0
dump sci all custom 10000 1.2.xyz id mol type q xs ys zs
run 0
undump sci
unfix 1
unfix 2
unfix 3
log 1.3_NVT.out
#-------------------------------------------------------------------------------
# Stage 1.3: NVT integration for 10 ps with a timestep of 1 fs
# Initial temperature 300 °C, final 30 °C
#-------------------------------------------------------------------------------
reset_timestep 0
thermo_style custom step v_time press vol v_sysdensity temp ebond eangle edihed eimp evdwl ecoul etail elong pe ke
thermo ${Nthermo}
fix 1 movable nvt temp 573.15 303.15 100 drag 0.0
fix 2 movable ave/time 1 1000 1000 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 1.3_averages.txt off 1
fix 3 movable ave/time 1 1 1 v_time c_thermo_temp c_thermo_press v_sysvol v_sysdensity v_etotal v_pe v_ke v_evdwl v_coulomb v_sxx v_syy v_szz v_syz v_sxz v_sxy file 1.3_instantaneous.txt
restart 10000 1.3.restart
dump trj all custom 1000 1.3.Trajectory.xyz id mol type q xs ys zs
fix trjE all ave/time 1000 1 1000 v_etotal v_pe v_ke file 1.3.energies.txt
dump sci all custom 1000 1.3.xyz id mol type q xs ys zs
timestep 1
run 10000
undump trj
unfix trjE
undump sci
restart 0
dump sci all custom 10000 1.3.xyz id mol type q xs ys zs
run 0
undump sci
unfix 1
unfix 2
unfix 3
log 1.4_Custom.out
#-------------------------------------------------------------------------------
# Stage 1.4: dipolemoment
#-------------------------------------------------------------------------------
# ----- !ìÃx -----
timestep 1.0 # fs
thermo_style custom step temp press etotal pe ke vol
thermo 1000
fix nvt all nvt temp 300.0 300.0 100.0
# sa 100 ps (100000 steps with 1 fs timestep)
# " 500 ps (500000 steps)
# ----- ©=vu M(t) = sum q_i * r_i -----
variable qx atom q*x
variable qy atom q*y
variable qz atom q*z
compute sumx all reduce sum v_qx
compute sumy all reduce sum v_qy
compute sumz all reduce sum v_qz
variable Mx equal c_sumx
variable My equal c_sumy
variable Mz equal c_sumz
variable Mabs equal sqrt(v_Mx*v_Mx + v_My*v_My + v_Mz*v_Mz)
# ----- Ïe8ú M(t) -----
fix fM all ave/time 1 1 1 v_Mx v_My v_Mz v_Mabs file dipole_Mt.dat
# ----- êøÜGreenKubo ( -----
# 5000 ÞøÜ 5 ps ã1 fs timestep
fix fcorr all ave/correlate 1 1 5000 v_Mx v_My v_Mz type auto file dipole_corr.dat ave running
# ----- ÷L -----
run 100000 # equilibration 100 ps
unfix fcorr
fix fcorr all ave/correlate 1 1 5000 v_Mx v_My v_Mz type auto file dipole_corr.dat ave running
run 500000 # production 500 ps
unfix fM
unfix fcorr
unfix nvt