Dear Alex:
I want to use lammps to caculate the stress between two groups, and I use the command “compute group/group”. However, when I run the script there is error show that “Pair hybrid sub style does not support the single call”.
I look the error of the lammps Manual, and it says this is because I have used the 3-body ptential.
So, there are some other ways to solve this problem.
My script is below:
#lattice topology(include SiC)===========================================
lattice custom {aa} origin 0.1 0.25 0 a1 1 0 0 a2 0 1 0 a3 0 0 1 basis 0 0 0 basis 0.5 0.5 0.0 basis 0.0 0.5 0.5 basis 0.5 0.0 0.5 basis 0.25 0.25 0.25 basis 0.75 0.75 0.25 basis 0.25 0.75 0.75 basis 0.75 0.25 0.75 orient x 1 0 0 orient y 0 1 0 orient z 0 0 1
region SiC-Cu block 0 {l1} {l5} {l4} 0 {l6}
create_box 3 SiC-Cu
region SiC block INF INF {a} INF INF INF
create_atoms 2 region SiC basis 1 1 basis 2 1 basis 3 1 basis 4 1 basis 5 2 basis 6 2 basis 7 2 basis 8 2
group SiC region SiC
#lattice topology(include Cu )===========================================
lattice none
lattice fcc {bb} origin 0.15 0.35 0.1 orient x 1 0 0 orient y 0 -1 0 orient z 0 0 -1
region Cu block INF INF INF {a} INF INF
create_atoms 3 region Cu
group Cu region Cu
#the mass of atom(include Cu and SiC)=====================================
mass 1 12.01
mass 2 28.086
mass 3 63.55
#the part I want to compute=====================================
region displace1 block INF INF {a} {l74} INF INF
group displace1 region displace1
region displace2 block INF INF {l84} {a} INF INF
group displace2 region displace2
region SiC11 block INF INF {yhigh} INF INF INF
group SiC11 region SiC11
region Cu1 block INF INF INF {ylow} INF INF
group Cu1 region Cu1
group mix subtract all Cu1 SiC11
#the mass of atom(include Cu and SiC)=====================================
group Cu type 3
group SiC type 1 2
#potential===========================================
pair_style hybrid tersoff morse/opt 2.0 eam/alloy/opt
pair_coeff * * tersoff SiC.tersoff C Si NULL
pair_coeff * * eam/alloy/opt Cu01.eam.alloy NULL NULL Cu
pair_coeff 1*2 3 none
pair_coeff 1 3 morse/opt
pair_coeff 2 3 morse/opt #(shun-Fa Huang) Takatoshi Kato
neighbor 2 bin
neigh_modify every 1 delay 0 check yes
#compute the intial temp====================================================
compute 3dtemp all temp
compute 3dpress all pressure 3dtemp
#caculate atoms displacement====================================================
variable v equal vol
variable v1 equal {v}
compute disp all displace/atom
compute Eke all ke/atom
compute Epe all pe/atom
compute strs all stress/atom
compute strs2 all stress/atom virial
compute diffu all msd
compute Ek all reduce ave c_Eke
compute Energy all reduce ave c_Eke c_Epe
compute SiC all reduce sum c_strs[1] c_strs[2] c_strs[3] c_strs[4] c_strs[5] c_strs[6] c_strs2[2]
variable SiCmis equal "sqrt(((c_SiC[1]-c_SiC[2])^2+(c_SiC[2]-c_SiC[3])^2+(c_SiC[3]-c_SiC[1])^2+6.0*((c_SiC[4])^2+(c_SiC[5])^2+(c_SiC[6])^2))/2.0)"
compute mine displace1 group/group displace2
variable mutual equal "sqrt((c_mine[1])^2+(c_mine[2])^2+(c_mine[3])^2)"
#zhuanhuadanwei====================================================
variable factor equal {v}/count(all)*10000
variable strss3 atom c_strs2[2]/{factor}
variable s11 atom c_strs[1]/{factor}
variable s22 atom c_strs[2]/{factor}
variable s33 atom c_strs[3]/{factor}
variable s12 atom c_strs[4]/{factor}
variable s13 atom c_strs[5]/{factor}
variable s23 atom c_strs[6]/{factor}
variable mises atom "sqrt(((v_s11-v_s22)^2+(v_s22-v_s33)^2+(v_s33-v_s11)^2+6.0*((v_s12)^2+(v_s13)^2+(v_s23)^2))/2.0)"
variable Fx equal c_SiC[1]/({v}*10000)
variable Fy equal c_SiC[2]/{v}/10000
variable Fz equal c_SiC[3]/{v}/10000
#minimization=================================================
thermo_style custom step temp c_3dtemp etotal c_Ek v_Fx v_Fy v_Fz lx ly lz v_v
thermo_modify lost warn line one norm yes flush yes
timestep 0.001
thermo 200
min_style cg
min_modify dmax 0.1 line backtrack
minimize 1.0e-10 1.0e-12 10000 5000
fix 6 all box/relax x 0.0 z 0.0 vmax 0.001
unfix 6
reset_timestep 0
#initial velocity & equilibrate======================
velocity all create {Tem1} 5812775 dist gaussian temp 3dtemp units box
fix 10 all ave/time 1 {t4} {t3} c_diffu[1] c_diffu[2] c_diffu[3] c_diffu[4] file SiC_{l0}.2.msd
fix 15 all ave/time {t4} 1 {t3} c_SiC[2] v_v c_Energy[1] c_Energy[2] file SiC_{l0}.2.stress
fix 16 displace1 ave/time {t4} 1 {t3} c_mine[1] c_mine[2] c_mine[3] v_mutual file SiC_{l0}.3.stress
fix 1 mix nvt temp {Tem1} {Tem1} 0.05 drag 1.0
fix 2 Cu1 nve
fix 3 Cu1 setforce 0 0 0
fix 4 SiC11 nve
fix 5 SiC11 setforce 0 0 0
thermo 1000
run $r
unfix 15
unfix 10
#unfix 16
reset_timestep 0
write_restart sic-cu-shear.restart
Load Process ======================================================
#tension strain rate=10e9/s
#tensile loading along y axis
0.2% strain in each step ==========================================
variable tl equal 200
variable delta equal 0.04
fix 20 displace1 ave/time {t4} 1 {t3} c_SiC[2] v_v c_mine[1] c_mine[2] c_mine[3] v_mutual c_Energy[1] c_Energy[2] file SiC_{l0}.stress
fix 22 all ave/time 1 {te} {te} c_diffu[1] c_diffu[2] c_diffu[3] c_diffu[4] file SiC_{l0}.msd
fix 13 displace1 ave/time {t4} 1 {t3} c_mine[1] c_mine[2] c_mine[3] v_mutual file SiC_{l0}.4.stress
dump mydump1 all cfg {te} SiC_{l0}.{l9}.*.cfg id type xsu ysu zsu vx vy vz c_disp[4] c_Eke c_Epe v_mises v_s12 v_s22 v_s13
dump_modify mydump1 element C Si Cu
label tloop
variable t loop {tl}
#fix 2 all deform 1 yz erate 0.001 remap v units box
displace_atoms mix ramp x {a} {delta} y {ylow}{bb} {yhigh}{aa} units box
displace_atoms SiC11 move {delta} {a} {a} units box
thermo ${te}
run n
#run {n1}
variable i equal $t
next t
jump sic-cu-shear.in tloop
variable t delete
variable strain21 delete
variable tl delete
variable delta delete
variable i delete
unfix 20
unfix 22
undump mydump1
unfix 1
unfix 2
unfix 3
unfix 4
unfix 5
u
nfix 20
unfix 22
clear
Simen Zhou
School of Aeronautics sicence and engineering
Beihang university
Xueyuan road No 31. Haidian district. Beijing