本文以及所有腳本版權都歸李雲海博士，裡面提到的腳本可以免費使用！
第一種實現方法腳本:以石墨烯為例，其晶格常數約為2.47埃。我們在2.30埃和2.70埃之間，每隔0.05埃取一個點。腳本中主循環如下(opt_1.sh)：
for a in $(seq 2.30 0.05 2.70); do
geninp
mpirun -np 20 vasp_std &> log.$a
wait
mv OUTCAR OUTCAR.$a
done
每一個晶格常數對應的POSCAR通過geninp函數生成：
#! /bin/bash
function geninp {
local a1x=$(awk 'BEGIN {print '$a' * 0.5}')
local a1y=$(awk 'BEGIN {print '$a' * 0.5 * sqrt(3)}')
local a2x=$(awk 'BEGIN {print '$a' * -0.5}')
local a2y=$(awk 'BEGIN {print '$a' * 0.5 * sqrt(3)}')
cat > POSCAR << EOF
Graphene
1.0
$a1x $a1y 0.0
$a2x $a2y 0.0
0.0 0.0 15.0
C
2
Direct
0.0 0.0 0.5
0.333333333 0.333333333 0.5
EOF
}
需要注意的幾個地方：
① 因為變量a是全局變量，所以在函數geninp中可以直接使用
② 因為Shell不支持浮點計算，所以要用awk來計算原胞基矢的各個分量，在awk命令中引用變量需要加單引號； ③ cat和EOF之間的部分相當於模板，將模板中的變量用相應的值替換，就生成了POSACR；
如果POSCAR中原子較多，再把模板放到腳本中，腳本就顯得十分冗長。這時可以用宏展開程序m4生成POSCAR。首先創建模板POSCAR.m4，內容如下：
Graphene
1.0
A1X A1Y 0.0
A2X A2Y 0.0
0.0 0.0 15.0
C
2
Direct
0.0 0.0 0.5
0.333333333 0.333333333 0.5
其中A1X, A1Y, A2X, A2Y都是待展開的宏。這些宏的定義在geninp中給出(opt_2.sh)：
#! /bin/bash
function geninp {
local a1x=$(awk 'BEGIN {print '$a' * 0.5}')
local a1y=$(awk 'BEGIN {print '$a' * 0.5 * sqrt(3)}')
local a2x=$(awk 'BEGIN {print '$a' * -0.5}')
local a2y=$(awk 'BEGIN {print '$a' * 0.5 * sqrt(3)}')
m4 -DA1X=$a1x -DA1Y=$a1y -DA2X=$a2x -DA2Y=$a2y POSCAR.m4 > POSCAR
}
for a in $(seq 2.30 0.05 2.70); do
geninp
mpirun -np 20 vasp_std &> log.$a
wait
mv OUTCAR OUTCAR.$a
done
-DFOO=bar的意思是將宏FOO定義為bar。因為m4將處理後的文本直接寫入到標準輸出，所以末尾還要重定向到POSCAR。
第二種實現方法腳本：由於宏處理器m4比較艱深晦澀，對數學計算的支持也不是很完美。因此我決定用Python自己寫一個宏處理器。核心算法是字符串遞歸替換，因此主程序（mace.py）非常簡單：
"""
Library for generating files using MACro Expansion.
You need to call the main() function, which accepts the following arguments:
macro: dictionary, macro definitions
template: string, filename of the template
output: string, filename of the output file
An auxiliary function include() is provided for including a file and defining
a macro from it.
"""
def expand_line(macro, line):
flag = False
for mac_name, mac_text in macro.items():
word = "<%s>" % mac_name
if line.find(word) is not -1:
line = line.replace(word, str(mac_text).lstrip("\n").rstrip("\n"))
flag = True
return line, flag
def include(filename, nl0=None, nl1=None):
with open(filename, "r") as infile:
content = infile.readlines()
nl_start = nl0 if nl0 is not None else 1
nl_end = nl1 if nl1 is not None else len(content)
longline = "".join(content[(nl_start-1):nl_end])
return longline
def main(macro, template, output):
with open(template, "r") as in_file:
content_raw = in_file.readlines()
content_expanded = []
for line in content_raw:
flag = True
while flag is True:
line, flag = expand_line(macro, line)
content_expanded.append(line)
with open(output, "w") as out_file:
for line in content_expanded:
out_file.write(line)
使用時只需編寫一個腳本（runmace.py）調用main函數即可。調用main函數需要三個參數：存儲宏定義的字典macro，模板文件名template和生成的輸入文件名output。除main函數外，還提供了一個include函數，用於把指定文件中的內容包含進模板。include函數除文件名filename外，另有兩個可選參數nl0和nl1用於指定行號，即可以只包含從開始行到結束行的內容（含這兩行）。行號從1開始，nl0默認為第一行，nl1默認為最後一行。runmace.py示例如下：
#! /usr/bin/env python
from mace import main
from math import sqrt
m = dict()
for a in [2.46, 2.47, 2.48]:
m["ax"] = a * 0.5
m["ay"] = a * 0.5 * sqrt(3.0)
m["bx"] = -m["ax"]
m["by"] = m["ay"]
for key, val in m.items():
m[key] = str("%14.9f" % val)
main(m, "POSCAR.tpl", "POSCAR.%4.2f" % a)
接下來給一個實用性較強的例子，首先是生成不同晶格常數下的石墨烯POSCAR。這種操作在優化晶體結構時很常見，所用模板（POSCAR.tpl）如下：
Graphene
1.0
<ax><ay> 0.000000000
<bx><by> 0.000000000
0.000000000 0.000000000 15.000000000
C
2
Direct
0.000000000 0.000000000 0.000000000
0.333333333 0.333333333 0.333333333
運行上面的run_mace.py腳本即可產生POSCAR.2.46、POSCAR.2.47和POSCAR.2.48，其中一個內容如下：
Graphene
1.0
1.230000000 2.130422493 0.000000000
-1.230000000 2.130422493 0.000000000
0.000000000 0.000000000 15.000000000
C
2
Direct
0.000000000 0.000000000 0.000000000
0.333333333 0.333333333 0.333333333
李雲海博士科學網博客地址：
http://blog.sciencenet.cn/blog-2909108-1176842.html
李雲海博士程序源碼分享網站：
https://gitee.com/yhli/misc
分享一個算平均電勢和面內積分電荷密度的程序
MobaXterm詳細使用教程系列二