该文章由
D
r
.
B
I
N
B
O
\rm \color{Violet} Dr.BINBO
Dr.BINBO 撰写，主要介绍
l
a
m
m
p
s
\color{Violet} lammps
lammps 实现剪切作用（两种方法）
?
\heartsuit
?
E
m
a
i
l
：
A
m
e
l
i
a
990608
@
163.
c
o
m
\rm Email：Amelia990608@163.com
Email：Amelia990608@163.com
N O T E ? ! ! ! ￣ \rm \underline {NOTE\ !!!} NOTE?!!!? 力 \color{red}{力} 力 学 \color{blue}{学} 学 专 \color{orange}{专} 专 题 \color{green}{题} 题！
感
谢
\rm 感谢
感谢
?
\heartsuit
?
D
r
.
B
I
N
B
O
\rm \color{Violet} Dr.BINBO
Dr.BINBO
?
\heartsuit
?
提
供
的
力
学
相
关
博
文
！
这
里
做
简
单
整
理
为
小
专
题
:
\rm 提供的力学相关博文！这里做简单整理为小专题:
提供的力学相关博文！这里做简单整理为小专题:
《
力
学
专
题
—
纳
米
压
痕
、
单
轴
压
缩
、
剪
切
行
为
的
L
A
M
M
P
S
模
拟
》
￣
\rm \color{red}\underline {《力学专题—纳米压痕、单轴压缩、剪切行为的LAMMPS模拟》}
《力学专题—纳米压痕、单轴压缩、剪切行为的LAMMPS模拟》?:
?
\dag
? 《Lammps计算纳米压痕—包含全部In文件》
?
\dag
?
?
\dag
? 《石墨烯复合材料单轴压缩-应力应变曲线以及弹性模量计算（包含in文件）》
?
\dag
?
?
\dag
?
?
\dag
? 《Lammps剪切作用的两种实现方法及对比—Deform以及Velocity》
Lammps剪切作用的两种实现方法及对比—Deform以及Velocity
一、Fix deform实现剪切作用
1. 初始设置—达到平衡态
############################################################################
################################# Deform method ##########################
############################################################################
#------------------------------Model settings--------------------------
clear
shell mkdir output ./output/minimize ./output/relax ./output/shear ./output/data
#------------------------------Basic settings--------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
variable Temp equal 300 # desired temperature
variable tstp equal 0.001 # timestep
variable tdamp equal 100*${tstp} # damping parameters #
variable pdamp equal 1000*${tstp} # for NPT ensemble #
variable dump_every equal 1000
variable equi_step equal 10000
read_data Ni_supercell.lmp
#----------------------------------------- Force Field ------------------------------------
pair_style eam
pair_coeff * * Ni_u3.eam
#--------------------------------------------- Minimization ---------------------------------------
thermo 100
timestep ${tstp}
thermo_style custom step temp etotal pe ke pxy
# compute specific quantities
compute pot_energy all pe/atom
dump min all custom ${dump_every} ./output/minimize/dump_minimize.*.lammpstrj id type x y z c_pot_energy
min_style cg
minimize 1.0e-6 1.0e-6 10000 10000
write_data ./output/minimize/mini.data
#------------------------------------------ equilibrium -------------------------------
undump min
reset_timestep 0
compute new3d all temp
compute new2d all temp/partial 0 1 1
velocity all create ${Temp} 12345 units box
fix 1 all npt temp ${Temp} ${Temp} ${tdamp} iso 0.0 0.0 ${pdamp}
fix_modify 1 temp new3d
thermo 1000
thermo_modify temp new3d
run ${equi_step}
write_data ./output/relax/equil.data
write_restart ./output/relax/equil.restart
2. Deform剪切实现—xy方向
#------------------------------------------- shear in xy plane -------------------------------------
reset_timestep 0
unfix 1
change_box all triclinic
variable tmpx equal lx
variable tmpy equal ly
variable tmpz equal lz
variable L0x equal ${tmpx}
variable L0y equal ${tmpy}
variable L0z equal ${tmpz}
# In units metal, pressure is in [bars]; sxz is in GPa
variable srate equal 0.01
variable final_strain equal 0.4
variable dyna_step equal ${final_strain}/(${srate}*${tstp})
variable sxy equal "-pxy/10000" # calculate shear stress
variable strain equal ${srate}*step*${tstp} # calculate shear strain
# deformation starts!
#fix 1 all nve
#fix 2 all temp/rescale 1 ${Temp} ${Temp} ${tdamp} 0.5
fix 2 all npt temp ${Temp} ${Temp} ${tdamp} iso 0.0 0.0 ${pdamp}
fix_modify 2 temp new2d
fix 3 all deform 1 xy erate ${srate} remap x
dump result all custom ${dump_every} ./output/shear/dump.*.lammps id type x y z
fix def1 all print 100 "${strain} ${sxy}" file ./output/data/Shear_data.txt screen no
thermo_style custom step press pe pxx pyy pzz pxy pxz pyz v_strain v_sxy
run ${dyna_step}
write_data ./output/shear/shear_end.*.data
write_restart ./output/shear/shear_end.*.restart
print "-------------------------------- All done! -----------------------------------------"
3. Deform剪切实现—应力应变分析
二、Velocity 实现剪切作用
1. 初始设置—达到平衡态
############################################################################
################################# Velocity method ##########################
############################################################################
#------------------------------Model settings--------------------------
clear
shell mkdir output ./output/minimize ./output/relax ./output/shear ./output/data
#------------------------------Basic settings--------------------------
units metal
dimension 3
boundary s s p
atom_style atomic
variable Temp equal 10 # desired temperature
variable tstp equal 0.001 # timestep
variable tdamp equal 100*${tstp} # damping parameters #
variable pdamp equal 1000*${tstp} # for NPT ensemble #
variable dump_every equal 1000
variable equi_step equal 10000
read_data Ni_supercell.lmp
variable size_x equal xhi-xlo
variable size_y equal yhi-ylo
variable size_z equal zhi-zlo
variable layer_thickness equal 0.1*${size_y}
variable upper_layer_limit equal ${size_y}-${layer_thickness}
variable lower_layer_limit equal ${layer_thickness}
# Define regions
region upper block INF INF ${upper_layer_limit} INF INF INF units box
region lower block INF INF INF ${lower_layer_limit} INF INF units box
# Define groups
group upper region upper
group lower region lower
group boundary union upper lower
group mobile subtract all boundary
#----------------------------------------- Force Field ------------------------------------
pair_style eam
pair_coeff * * Ni_u3.eam
#--------------------------------------------- Minimization ---------------------------------------
thermo 100
timestep ${tstp}
thermo_style custom step temp etotal pe ke pxz
# fix boundary group
# compute specific quantities
compute pot_energy all pe/atom
compute stress all stress/atom NULL
dump min all custom ${dump_every} ./output/minimize/dump_minimize.*.lammpstrj id type x y z c_pot_energy
min_style cg
minimize 1.0e-6 1.0e-6 10000 10000
write_data ./output/minimize/mini.data
#------------------------------------------ equilibrium -------------------------------
undump min
reset_timestep 0
compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1
velocity mobile create ${Temp} 12345 units box
#fix 1 all npt temp ${Temp} ${Temp} ${tdamp} aniso 0.0 0.0 ${pdamp}
fix 1 all nve
fix 0 boundary setforce 0.0 0.0 0.0
fix 2 mobile temp/rescale 10 ${Temp} ${Temp} 10.0 1.0
fix_modify 2 temp new3d
#fix_modify 1 temp new3d
thermo 1000
thermo_modify temp new3d
run ${equi_step}
write_data ./output/relax/equil.data
write_restart ./output/relax/equil.restart
2. Velocity—剪切实现
#------------------------------------------- shear in xy plane -------------------------------------
#change_box all boundary s s p
reset_timestep 0
unfix 1
#fix 1 all nvt temp ${Temp} ${Temp} ${tdamp}
#fix_modify 1 temp new3d
#run 10000
#reset_timestep 0
#unfix 1
variable tmpx equal lx
variable tmpy equal ly
variable tmpz equal lz
variable L0x equal ${tmpx}
variable L0y equal ${tmpy}
variable L0z equal ${tmpz}
variable vol_ini equal v_L0x*v_L0y*v_L0z
print "Initial volumn, vol = ${vol_ini}"
# In units metal, pressure is in [bars]; sxx, syy, szz are in GPa
compute 2 mobile stress/atom NULL
compute strsxy mobile reduce sum c_2[4]
variable unit equal 1E-4
#variable sxy equal ${unit}*(c_strsxy/v_vol_ini)
variable deltaL equal (lx-${L0x})
variable sxy equal "-pxy/10000" # calculate shear stress
variable strain equal v_deltaL/ly # calculate shear strain
variable srate equal 0.01
variable vx equal ${srate}*${L0y}
variable final_strain equal 0.3
#variable dyna_step equal ${final_strain}/(${srate}*${tstp})
variable dyna_step equal 15000
# deformation starts!
fix 1 all nve
fix 2 mobile temp/rescale 10 ${Temp} ${Temp} 10.0 1.0
fix_modify 2 temp new2d
velocity upper set 2.0 0.0 0.0 units box
velocity mobile ramp vx 0.0 2.0 y ${lower_layer_limit} ${upper_layer_limit} sum yes units box
dump result all custom ${dump_every} ./output/shear/dump.*.lammps id type x y z
fix def1 all print 1000 "${strain} ${sxy}" file ./output/data/Shear_data.txt screen no
#
fix def1 all ave/time 100 10 1000 v_strain v_sxy file ./output/data/Shear_data.txt
thermo_style custom step press pe pxx pyy pzz pxy pxz pyz v_strain v_sxy
run ${dyna_step}
write_data ./output/shear/shear_end.*.data
write_restart ./output/shear/shear_end.*.restart
print "-------------------------------- All done! -----------------------------------------"
3. Velocity 剪切实现—应力应变分析
三、全部代码-剪切作用
1. Deform 方法
############################################################################
################################# Deform method ##########################
############################################################################
#------------------------------Model settings--------------------------
clear
shell mkdir output ./output/minimize ./output/relax ./output/shear ./output/data
#------------------------------Basic settings--------------------------
units metal
dimension 3
boundary p p p
atom_style atomic
variable Temp equal 300 # desired temperature
variable tstp equal 0.001 # timestep
variable tdamp equal 100*${tstp} # damping parameters #
variable pdamp equal 1000*${tstp} # for NPT ensemble #
variable dump_every equal 1000
variable equi_step equal 10000
read_data Ni_supercell.lmp
#----------------------------------------- Force Field ------------------------------------
pair_style eam
pair_coeff * * Ni_u3.eam
#--------------------------------------------- Minimization ---------------------------------------
thermo 100
timestep ${tstp}
thermo_style custom step temp etotal pe ke pxy
# compute specific quantities
compute pot_energy all pe/atom
dump min all custom ${dump_every} ./output/minimize/dump_minimize.*.lammpstrj id type x y z c_pot_energy
min_style cg
minimize 1.0e-6 1.0e-6 10000 10000
write_data ./output/minimize/mini.data
#------------------------------------------ equilibrium -------------------------------
undump min
reset_timestep 0
compute new3d all temp
compute new2d all temp/partial 0 1 1
velocity all create ${Temp} 12345 units box
fix 1 all npt temp ${Temp} ${Temp} ${tdamp} iso 0.0 0.0 ${pdamp}
fix_modify 1 temp new3d
thermo 1000
thermo_modify temp new3d
run ${equi_step}
write_data ./output/relax/equil.data
write_restart ./output/relax/equil.restart
#------------------------------------------- shear in xy plane -------------------------------------
reset_timestep 0
unfix 1
change_box all triclinic
variable tmpx equal lx
variable tmpy equal ly
variable tmpz equal lz
variable L0x equal ${tmpx}
variable L0y equal ${tmpy}
variable L0z equal ${tmpz}
# In units metal, pressure is in [bars]; sxz is in GPa
variable srate equal 0.01
variable final_strain equal 0.4
variable dyna_step equal ${final_strain}/(${srate}*${tstp})
variable sxy equal "-pxy/10000" # calculate shear stress
variable strain equal ${srate}*step*${tstp} # calculate shear strain
# deformation starts!
#fix 1 all nve
#fix 2 all temp/rescale 1 ${Temp} ${Temp} ${tdamp} 0.5
fix 2 all npt temp ${Temp} ${Temp} ${tdamp} iso 0.0 0.0 ${pdamp}
fix_modify 2 temp new2d
fix 3 all deform 1 xy erate ${srate} remap x
dump result all custom ${dump_every} ./output/shear/dump.*.lammps id type x y z
fix def1 all print 100 "${strain} ${sxy}" file ./output/data/Shear_data.txt screen no
thermo_style custom step press pe pxx pyy pzz pxy pxz pyz v_strain v_sxy
run ${dyna_step}
write_data ./output/shear/shear_end.*.data
write_restart ./output/shear/shear_end.*.restart
print "-------------------------------- All done! -----------------------------------------"
2. Velocity 方法
############################################################################
################################# Velocity method ##########################
############################################################################
#------------------------------Model settings--------------------------
clear
shell mkdir output ./output/minimize ./output/relax ./output/shear ./output/data
#------------------------------Basic settings--------------------------
units metal
dimension 3
boundary s s p
atom_style atomic
variable Temp equal 10 # desired temperature
variable tstp equal 0.001 # timestep
variable tdamp equal 100*${tstp} # damping parameters #
variable pdamp equal 1000*${tstp} # for NPT ensemble #
variable dump_every equal 1000
variable equi_step equal 10000
read_data Ni_supercell.lmp
variable size_x equal xhi-xlo
variable size_y equal yhi-ylo
variable size_z equal zhi-zlo
variable layer_thickness equal 0.1*${size_y}
variable upper_layer_limit equal ${size_y}-${layer_thickness}
variable lower_layer_limit equal ${layer_thickness}
# Define regions
region upper block INF INF ${upper_layer_limit} INF INF INF units box
region lower block INF INF INF ${lower_layer_limit} INF INF units box
# Define groups
group upper region upper
group lower region lower
group boundary union upper lower
group mobile subtract all boundary
#----------------------------------------- Force Field ------------------------------------
pair_style eam
pair_coeff * * Ni_u3.eam
#--------------------------------------------- Minimization ---------------------------------------
thermo 100
timestep ${tstp}
thermo_style custom step temp etotal pe ke pxz
# fix boundary group
# compute specific quantities
compute pot_energy all pe/atom
compute stress all stress/atom NULL
dump min all custom ${dump_every} ./output/minimize/dump_minimize.*.lammpstrj id type x y z c_pot_energy
min_style cg
minimize 1.0e-6 1.0e-6 10000 10000
write_data ./output/minimize/mini.data
#------------------------------------------ equilibrium -------------------------------
undump min
reset_timestep 0
compute new3d mobile temp
compute new2d mobile temp/partial 0 1 1
velocity mobile create ${Temp} 12345 units box
#fix 1 all npt temp ${Temp} ${Temp} ${tdamp} aniso 0.0 0.0 ${pdamp}
fix 1 all nve
fix 0 boundary setforce 0.0 0.0 0.0
fix 2 mobile temp/rescale 10 ${Temp} ${Temp} 10.0 1.0
fix_modify 2 temp new3d
#fix_modify 1 temp new3d
thermo 1000
thermo_modify temp new3d
run ${equi_step}
write_data ./output/relax/equil.data
write_restart ./output/relax/equil.restart
#------------------------------------------- shear in xy plane -------------------------------------
#change_box all boundary s s p
reset_timestep 0
unfix 1
#fix 1 all nvt temp ${Temp} ${Temp} ${tdamp}
#fix_modify 1 temp new3d
#run 10000
#reset_timestep 0
#unfix 1
variable tmpx equal lx
variable tmpy equal ly
variable tmpz equal lz
variable L0x equal ${tmpx}
variable L0y equal ${tmpy}
variable L0z equal ${tmpz}
variable vol_ini equal v_L0x*v_L0y*v_L0z
print "Initial volumn, vol = ${vol_ini}"
# In units metal, pressure is in [bars]; sxx, syy, szz are in GPa
compute 2 mobile stress/atom NULL
compute strsxy mobile reduce sum c_2[4]
variable unit equal 1E-4
#variable sxy equal ${unit}*(c_strsxy/v_vol_ini)
variable deltaL equal (lx-${L0x})
variable sxy equal "-pxy/10000" # calculate shear stress
variable strain equal v_deltaL/ly # calculate shear strain
variable srate equal 0.01
variable vx equal ${srate}*${L0y}
variable final_strain equal 0.3
#variable dyna_step equal ${final_strain}/(${srate}*${tstp})
variable dyna_step equal 15000
# deformation starts!
fix 1 all nve
fix 2 mobile temp/rescale 10 ${Temp} ${Temp} 10.0 1.0
fix_modify 2 temp new2d
velocity upper set 2.0 0.0 0.0 units box
velocity mobile ramp vx 0.0 2.0 y ${lower_layer_limit} ${upper_layer_limit} sum yes units box
dump result all custom ${dump_every} ./output/shear/dump.*.lammps id type x y z
fix def1 all print 1000 "${strain} ${sxy}" file ./output/data/Shear_data.txt screen no
#
fix def1 all ave/time 100 10 1000 v_strain v_sxy file ./output/data/Shear_data.txt
thermo_style custom step press pe pxx pyy pzz pxy pxz pyz v_strain v_sxy
run ${dyna_step}
write_data ./output/shear/shear_end.*.data
write_restart ./output/shear/shear_end.*.restart
print "-------------------------------- All done! -----------------------------------------"
四、Python绘图代码
1. Python 代码
具体请参考 ：《I论文—Python绘图模板及技巧（未完成）》
import os
def make_save_file(filepath,filename,savename):
# filepath 存储路径； filename：创建文件夹的名字 savename：存储图片的名字
save_path = filepath+os.sep+filename+os.sep+savename
all_path = filepath+os.sep+filename
if not os.path.exists(all_path):
os.mkdir(all_path)
plt.savefig(save_path,dpi=150)
filepath = r"PureNi_shear"
font1 = {'family': 'Times New Roman',
'weight': 'bold',
'size': 22,
}
font2 = {'family': 'Times New Roman',
'weight': 'normal',
'size': 11,
}
# Global setting
# ----------------------------------------------------------------------#
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
# coding:utf-8
import pylab
import matplotlib.ticker as plticker
from matplotlib.pyplot import MultipleLocator
import xlrd
import matplotlib.pyplot as plt
from matplotlib.lines import Line2D
import math
from scipy import optimize
plt.style.use('science')
# ----------------------------------------------------------------------#
# ----------------------------------------------------------------------#
def setup(ax,x_label,y_label):
# 边框线的线宽设置
width = 1.5
ax.spines['top'].set_linewidth(width)
ax.spines['bottom'].set_linewidth(width)
ax.spines['left'].set_linewidth(width)
ax.spines['right'].set_linewidth(width)
# 边框上的ticks的出现
ax.tick_params(top='on', bottom='on', left='on', right='on', direction='in')
# 边框上的ticks对应的lable,即数字的尺寸
ax.tick_params(labelsize=20,pad=6)
#ax.yaxis.set_ticks_position('right')
ax.tick_params(which='major', width=1.50, length=6)
ax.tick_params(which='minor', width=0.75, length=0)
# 边框上的ticks的出现的间隔
locx = plticker.MultipleLocator(base=0.1) # this locator puts ticks at regular intervals
ax.xaxis.set_major_locator(locx)
locy = plticker.MultipleLocator(base=10) # this locator puts ticks at regular intervals
ax.yaxis.set_major_locator(locy)
# 边框上的注释labels
labels = ax.get_xticklabels() + ax.get_yticklabels()
[label.set_fontname('Times New Roman') for label in labels]
# 边框上的注释labels距离坐标轴的距离
xAxisLable = x_label
yAxisLable = y_label
ax.set_xlabel(xAxisLable, font1, labelpad=5)
ax.set_ylabel(yAxisLable, font1, labelpad=5)
# ----------------------------------------------------------------------#
# 读取excel
excel_path = r"out.xlsx"
work_sheet = xlrd.open_workbook(excel_path);
# sheet的数量
sheet_num = work_sheet.nsheets;
# 获取sheet name
sheet_name = []
for sheet in work_sheet.sheets():
sheet_name.append(sheet.name)
# ----------------------------------------------------------------------#
# ----------------------------------------------------------------------#
for sheet_num in range(100):
# 选取sheet
now_sheet = work_sheet.sheets()[0]
#now_sheet_ref = work_sheet.sheets()[3]
# sheet的行
row = now_sheet.nrows
# sheet的列
ncol = now_sheet.ncols
#----------------------------------------------------------------------#
#----------------------------------------------------------------------#
font1 = {'family': 'Times New Roman',
'weight': 'bold',
'size': 18,
}
font2 = {'family': 'Times New Roman',
'weight': 'normal',
'size': 11,
}
#----------------------------------------------------------------------#
# temp
deform_strain = now_sheet.col_values(0)
deform_strain_matrix = deform_strain[1:]
deform_stress = now_sheet.col_values(1)
deform_stress_matrix = list(filter(None, deform_stress[1:]))
velocity_strain = now_sheet.col_values(2)
velocity_strain_matrix = list(filter(None, velocity_strain[1:]))
velocity_stress = now_sheet.col_values(3)
velocity_stress_matrix = list(filter(None, velocity_stress[1:]))
markersize = 10
linewidth = 1
markevery =2
alpha =0.5
p_1 = dict(marker='o',color = 'r',linestyle = 'none',markevery=markevery,markerfacecolor='r',markersize=markersize,linewidth=linewidth,alpha=alpha)
p_2 = dict(marker='s',color = 'b',linestyle = 'none',markevery=markevery,markerfacecolor='none',markersize=markersize,linewidth=linewidth,alpha=alpha)
#----------------------------------------------------------------------#
fig, axe = plt.subplots(1, 1, figsize=(5,5))
axe.plot(deform_strain_matrix[0:sheet_num], deform_stress_matrix[0:sheet_num],**p_1,label='Deform method')
#axe.plot(velocity_strain_matrix, velocity_stress_matrix,**p_2,label='Velocity method')
plt.axvline(deform_strain_matrix[sheet_num],color='k',lw=2,linestyle='dashed')
axe.set_ylim(0, 20)
axe.set_xlim(0, 0.2)
setup(axe,r'Sheer strain',r'$\rm Stress\ /Gpa$')
axe.legend(loc='best', frameon=False, \
labelspacing=0.3)
axe.legend(bbox_to_anchor=(0.4, 1), \
loc='upper left', borderaxespad=0.,fontsize='xx-large')
# 图片保存
filename = "deform"
savename = str(sheet_num)+".jpg"
make_save_file(filepath,filename,savename)
# ----------------------------------------------------------------------#
plt.show()
2. 数据下载
链接：（请点击）
提取码：9rpj