-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathin_read_fix_Siabox.fcc
77 lines (75 loc) · 3.88 KB
/
in_read_fix_Siabox.fcc
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
#xyz文件中只有原子的类型和坐标因此使用read_dump命令重启就一定要提前定义好模拟盒子，势函数，邻近原子列表
#弊端：xyz文件中没有原子速度因此如果没有velocity命令那么原子的速度都为0！
#----------------------initiation---------------------------------------
units metal
boundary p p p
#----------------------atom definition-------------------------------------
lattice fcc ${A}
variable ax equal ${A}
variable strain equal ${Strain}
variable d1 equal "v_ax/sqrt(2)"
region whole block -10 10 -10 10 -10 10
create_box 1 whole
change_box all x scale ${Strain} remap
fix 1 all deform 1 x scale 1 y scale 1 z scale 1
pair_style eam/alloy
pair_coeff * * FeNiCr_Bonny_2013_ptDef.eam.alloy Ni
neighbor 2.0 bin
neigh_modify delay 5 check yes one 100000 page 100000000
read_dump ${Siabox} 0 x y z box no add yes format xyz
#读取xyz文件
#0为读入文件对应的步数，但是因为某些原因只能设为0.
#------------------------------定义间隙--------------------------------
variable dis equal 0.27
variable dx atom "x/v_strain/v_ax - floor(x/v_strain/v_ax)"
variable dy atom "y/v_ax - floor(y/v_ax)"
variable dz atom "z/v_ax - floor(z/v_ax)"
variable SIAs atom " sqrt((v_dx)^2+(v_dy)^2+(v_dz)^2) > v_dis && &
sqrt((v_dx-0.5)^2+(v_dy)^2+(v_dz-0.5)^2) > v_dis && &
sqrt((v_dx-0.5)^2+(v_dy-0.5)^2+(v_dz)^2) > v_dis && &
sqrt((v_dx)^2+(v_dy-0.5)^2+(v_dz-0.5)^2) > v_dis && &
sqrt((v_dx-1)^2+(v_dy-0.5)^2+(v_dz-0.5)^2) > v_dis && &
sqrt((v_dx-0.5)^2+(v_dy-1)^2+(v_dz-0.5)^2) > v_dis && &
sqrt((v_dx-0.5)^2+(v_dy-0.5)^2+(v_dz-1)^2) > v_dis && &
sqrt((v_dx)^2+(v_dy)^2+(v_dz-1)^2) > v_dis && &
sqrt((v_dx)^2+(v_dy-1)^2+(v_dz)^2) > v_dis && &
sqrt((v_dx-1)^2+(v_dy)^2+(v_dz)^2) > v_dis && &
sqrt((v_dx)^2+(v_dy-1)^2+(v_dz-1)^2) > v_dis && &
sqrt((v_dx-1)^2+(v_dy-1)^2+(v_dz)^2) > v_dis && &
sqrt((v_dx-1)^2+(v_dy)^2+(v_dz-1)^2) > v_dis && &
sqrt((v_dx-1)^2+(v_dy-1)^2+(v_dz-1)^2) > v_dis"
#-------------------------------输出间隙--------------------------------
group iniSIAs variable SIAs
group SIAs dynamic all var SIAs every 1
compute clus SIAs cluster/atom ${d1}
variable numSIAs equal count(SIAs)
#-----------------------------MSD-------------------------------
compute msd all msd
#------------------------------Fixed MD-------------------------------
fix Fixed SIAs setforce 0 0 0
group normal subtract all iniSIAs
velocity normal create ${Temp} ${Nseed} rot yes mom yes dist gaussian
velocity all zero linear
velocity all zero angular
dump fixeddump all xyz 1 fixall.*.xyz
fix nvt all nvt temp ${Temp} ${Temp} 1
thermo 1
log log.fix.${Action}.${m}
thermo_style custom step temp etotal pe press v_numSIAs
run ${Fixedstep}
#-------------------------------赋予速度并开始真正MD----------------------------------------
reset_timestep 0
#----- relaxation --------------------------------------------------------------------------
unfix Fixed
undump fixeddump
log log.${Action}.${m}
thermo 100
thermo_style custom step time temp press pxx pyy pzz pxy pxz pyz v_numSIAs c_msd[4] etotal pe vol lx
timestep 0.002
compute 2 all stress/atom NULL
dump 2 all custom 20000 stress.${Action}.temp${Temp}.strain${Strain}.${Nseed}.* x y z c_2[1] c_2[2] c_2[3] c_2[4] c_2[5] c_2[6]
dump dumpposition all xyz ${gap} all.*.xyz
dump dumpsia SIAs custom ${gap} sia.*.xyz x y z c_clus
restart 250000 restart.*
write_restart restart.*
run ${Nstep}