lammps实例(2)

表面与界面能

铜的表面能

当物体形成表面时，表面上的原子键发生断裂，接近表面的几层原子不再如之前处于平衡状态，从而导致能量的升高，升高的温度便是物体的表面能。

利用LAMMPS 做出 20*20*40 fcc 的盒子，删去边缘的原子制造出一段真空层；算出此时体系的总能量E0

，

然后从中间把盒子切成两半并移至足够远的距离，此时的体系总能量为

Efinal

，

从而表面能：                      

γsurface=

Efinal−E0

2A

A为表面的面积

(100) 面与 (111) 面

如下是输入文件in.surface_Cu_100

# LAMMPS Cu _Surface_100 units metal boundary p p p atom_style atomic lattice fcc 3.61 region box block 0 20 0 20 0 40 create_box 1 box create_atoms 1 box timestep 0.005 thermo 5

pair_style eam/alloy pair_coeff * * jin_copper_lammps.setfl Cu

region boundary1 block INF INF INF INF 29.9 INF region boundary2 block INF INF INF INF INF 9.9 group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2

delete_atoms group boundary

neighbor 0.6 bin

neigh_modify every 5 delay 0 check yes

compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0

dump 1 all custom 100 dump.atom id xs ys zs c_3 c_4 c_5

dump_modify 1 format \"%d %16.9g %16.9g %16.9g %16.9g %16.9g %g\"

min_style sd

minimize 1.0e-30 1.0e-15 1000 10000 variable E equal pe print \"---------------------- E=$E---------------------------\" run 0

region down block INF INF INF INF INF 19.94 region up block INF INF INF INF 19.95 INF group up region up group down region down

displace_box all z delta 0 40 units lattice remap none displace_atoms up move 0 0 40 units lattice minimize 1.0e-30 1.0e-20 10000 100000 print \"-----SURFACE----------E=$E--------------------------------\"

Plane (100) (111) Surface energy(mJ/m2)

1330 1228

Plane (111)

Plane (100)

镁的hcp/fcc界面能

计算镁的hcp/fcc界面能比较有趣。正如我们所知，镁是hcp结构的金属，它是按

ABABABABAB… 所排列的；当我们进行一系列移动之后得到ABABCABCAB… ，从而产生镁的hcp/fcc界面。

如下是输入文件 in.hcp_fcc_vac

# LAMMPS hcp_fcc_Mg

units metal boundary p p p atom_style atomic lattice custom 4.5198 a1 1.2247448713915049 0 0 a2 0 0.70710678118675244 0 a3 0 0 1.1476 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.33333333333333333333 0.0 0.5 basis 0.83333333333333333333 0.5 0.5 region box block 0 6 0 6 0 36 create_box 1 box create_atoms 1 box pair_style eam/alloy

pair_coeff * * mg.eam.alloy Mg neighbor 0.8 bin

neigh_modify every 5 delay 0 check yes thermo 10 variable N equal atoms variable E equal pe

region boundary1 block INF INF INF INF 26.9 INF region boundary2 block INF INF INF INF INF 8.9 group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2 delete_atoms group boundary

label dj

variable x index 17.9 18.9 19.9 20.9 21.9 22.9 23.9 24.9 25.9 region up_$x block INF INF INF INF $x INF group up_$x region up_$x displace_atoms up_$x move

-0.3333333333333333333333333333333333333333333333333333333 0 0 units lattice next x jump in.hcp_fcc_vac dj

timestep 0.01

region part block INF INF INF INF 12 24 group part region part compute 3 all pe/atom compute 4 all ke/atom compute 5 all coord/atom 3.0 dump 1 part custom 100 dump.atom id xs ys zs c_3 c_4 c_5

dump_modify 1 format \"%d %16.9g %16.9g %16.9g %16.9g %16.9g %g\" min_style sd

minimize 1.0e-30 1.0e-15 10000 100000 shell mkdir dump_fcc_hcp shell mv dump.atom dump_fcc_hcp/dump.atom print \"----------fcc_hcp---------------N=$N E=$E------------------------\"

delete_atoms group all lattice none

lattice custom 4.5198 a1 1.2247448713915049 0 0 a2 0 0.70710678118675244 0 a3 0 0 1.1476 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.33333333333333333333 0.0 0.5 basis 0.8333333333333333333333333333333333333333333333333333 0.5 0.5 create_atoms 1 box group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2 delete_atoms group boundary undump 1 dump 1 all custom 100 dump_hcp.atom id xs ys zs c_3 c_4 c_5 dump_modify 1 format \"%d %16.9g %16.9g %16.9g %16.9g %16.9g %g\" min_style sd

minimize 1.0e-30 1.0e-15 10000 10000 print \"----------hcp---------------N=$N E=$E------------------------\" shell mkdir dump_hcp shell mv dump_hcp.atom dump_hcp/dump.atom

delete_atoms group all lattice none

lattice custom 4.5198 a1 1.2247448713915049 0 0 a2 0 0.70710678118675244 0 a3 0 0 1.1476 basis 0.0 0.0 0.0 basis 0.5 0.5 0.0 basis 0.33333333333333333333 0.0 0.5 basis 0.83333333333333333333 0.5 0.5 create_atoms 1 box group boundary1 region boundary1 group boundary2 region boundary2 group boundary union boundary1 boundary2 delete_atoms group boundary

label ddj

variable y index 9.9 10.9 11.9 12.9 13.9 14.9 15.9 16.9 17.9 18.9 19.9 20.9 21.9 22.9 23.9 24.9 25.9 region up_fcc_$y block INF INF INF INF $y INF group up_fcc_$y region up_fcc_$y displace_atoms up_fcc_$y move

-0.333333333333333333333333333333333333333333333333333333333 0 0 units lattice next y jump in.hcp_fcc_vac ddj undump 1 dump 1 all custom 1000 dump_fc.atom id xs ys zs c_3 c_4 c_5 dump_modify 1 format \"%d %16.9g %16.9g %16.9g %16.9g %16.9g %g\"

min_style sd

minimize 1.0e-30 1.0e-15 10000 100000

shell mkdir dump_fcc shell mv dump_fcc.atom dump_fcc/dump.atom print \"-----------fcc--------------N=$N E=$E------------------------\"

根据公式：            

γ=

Et−

1

(E+Efcc)2hcp

2A

可以得到：γ=27.058(mJ/m) ，这与之前所得到结果比较吻合。

2