理论计算常用脚本

发表于 更新于 分类于 ， 本文字数： 12k 阅读时长 ≈ 11 分钟

记录平时在Windows和CentOS系统中使用Gaussian、Orca、Multiwfn和VMD等软件进行理论计算所使用的脚本。
CentOS系统已安装PBS队列系统、Multiwfn软件。

脚本使用方法

将下列命令保存为pbs，并赋予执行权限，将pbs所在文件夹的路径添加至系统路径中，如export PATH=$PATH:/home/script；按照需要将本文中其他脚本保存至pbs所在文件夹，赋予脚本执行权限；使用脚本时，直接输入pbs命令按照要求再输入对应的数字即可。

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
#! /bin/bash
echo "-----------------选项-----------------"
 
echo "please input the number for selecting country:
1. 提交当前文件夹中的gjf文件(Gaussian>log)
2. 汇总当前文件夹中的吸收数据(Gaussian)
3. 汇总当前文件夹中的HOMO-LUMO数据(Gaussian)
4. 提交当前文件夹中的inp文件(orca>out)
5. 提交当前文件夹中的同名dal/mol文件(dalton>out)
6. 输出当前文件夹中的HOMO-LUMO-cube文件(Gaussian)"

value=0;
read -p "input: " value
case $value in
    1) echo "Gaussian任务已提交至PBS，具体查看pbs.log"
	gaussian-pbs
        ;;
    2) abs-read
	echo "吸收数据已汇总至当前文件夹abs-read.txt文件中"
        ;;
    3) HOMO-LUMO
	echo "HOMO-LUMO数据已汇总至当前文件夹HOMO-LUMO.txt文件中"
        ;;
    4)	orca-pbs
	echo "orca计算任务已提交至PBS，具体查看pbs.log"
        ;;
    5) 
	read -p "请输入所需CPU核心数（1~48）：" valuecpu
	read -p "请输入每个核心所需内存（MB）：" valuemem
	dalton-pbs $valuecpu $valuemem
	echo "Dalton计算任务已提交至PBS，具体查看pbs.log"
        ;;
    6) mo-cube
        ;;
    *) echo "You select number is not in the menu"
        exit 1
        ;;
esac

PBS批量提交Gaussian任务

该脚本会自动识别当前文件夹中的gjf文件，并根据log文件判断该gjf文件是否被计算过；任务提交和计算完成时会在pbs.log文件中添加记录，计算完成后会自动转换生成fchk文件。

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
#!/bin/bash
for inf in *.gjf
do
path=$(cd "$(dirname "$inf")";pwd)
if [ ! -f ${inf//gjf/log} ];  ##判断gjf文件是否存在对应的log文件
then
cat > pbs.sh << "EOF"
#PBS -N gaussian
#PBS -l nodes=1:ppn=12
#PBS -l mem=24000mb
#PBS -q normal
#PBS -j oe
cd $PBS_O_WORKDIR
INPUT_NAME=1       ##这里输入文件名，不含拓展名，不需要修改##
g09 $INPUT_NAME.gjf   ##对应G09软件
wait
formchk $INPUT_NAME.chk
echo `date +"%Y-%m-%d %H:%M:%S"` job任务已经计算完成！ >> /home/pbs.log
EOF
name="INPUT_NAME="${inf%.*}
filename="$path/$inf"
sed -i -r "s/(INPUT_NAME=)[^#]*/$name    /" pbs.sh
sed -i -r "s?job?$filename?" pbs.sh   ##替换pbs.sh中文件路径
qsub pbs.sh
rm -f pbs.sh
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经提交至PBS队列！ >> /home/pbs.log
else
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经计算过了！ >> /home/pbs.log
fi
done

PBS批量提交Orca任务

该脚本会自动识别当前文件夹中的inp文件，并根据out文件判断该inp文件是否被计算过；任务提交和计算完成时会在pbs.log文件中添加记录，计算完成后会自动转换生成molden.input文件。需要根据实际安装情况修改orca的安装路径。

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
#!/bin/bash
for inf in *.inp
do
path=$(cd "$(dirname "$inf")";pwd)
if [ ! -f ${inf//inp/out} ];  ##判断inp文件是否存在对应的out文件
then
cat > pbs.sh << "EOF"
#PBS -S /bin/bash
#PBS -N orca
#PBS -l nodes=1:ppn=12
#PBS -l mem=80000mb
#PBS -q single
#PBS -j oe
cd $PBS_O_WORKDIR
INPUT_NAME=1       ##这里输入文件名，不含拓展名##
export OMP_NUM_THREADS=12
/home/orcauser/software/orca502-static/orca $INPUT_NAME.inp > $INPUT_NAME.out
wait
/home/orcauser/software/orca502-static/orca_2mkl $INPUT_NAME -molden
echo `date +"%Y-%m-%d %H:%M:%S"` job任务已经计算完成！来自Orca！ >> /home/pbs.log
EOF
name="INPUT_NAME="${inf%.*}
filename="$path/$inf"
sed -i -r "s/(INPUT_NAME=)[^#]*/$name    /" pbs.sh
sed -i -r "s?job?$filename?" pbs.sh   ##替换pbs.sh中文件路径
qsub pbs.sh
rm -f pbs.sh
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经提交至PBS队列！来自Orca！ >> /home/pbs.log
else
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经计算过了！来自Orca！ >> /home/pbs.log
fi
done

激发态能量(吸收值)汇总

读取当前文件夹中所有的out或log文件，将激发态能量汇总至abs-read.txt文件。

1
2
3
4
5
6
7
#!/bin/bash
for inf in *.chk
do
echo "当前数据为：${inf%.*}" >> abs-read.txt
grep -E 'Excited State\ \ ' ${inf//chk/out}  >> abs-read.txt
grep -E 'Excited State\ \ ' ${inf//chk/log}  >> abs-read.txt
done

HOMO-LUMO能量汇总(Gaussian)

读取当前文件夹中所有的fchk文件，将HOMO-LUMO能量以及能量差汇总至HOMO-LUMO.txt文件，需要提前安装Multiwfn。

1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/bash
cat > script.txt << EOF 
0
EOF
for inf in *.fchk
do
multiwfn $inf < script.txt > ${inf//fchk/txt}   #Multiwfn按照实际情况修改大小写
echo "当前数据为：${inf%.*}" >> HOMO-LUMO.txt
grep -E 'HOMO|LUMO' ${inf//fchk/txt}  >> HOMO-LUMO.txt
rm -f ${inf//fchk/txt}
done
rm -f script.txt

HOMO-LUMO能量汇总(Orca)

读取当前文件夹中所有的molden.input文件，将HOMO-LUMO能量以及能量差汇总至HOMO-LUMO.txt文件，需要提前安装Multiwfn。

1
2
3
4
5
6
7
8
9
10
11
12
#!/bin/bash
cat > script.txt << EOF 
0
EOF
for inf in *.molden.input
do
multiwfn $inf < script.txt > ${inf//molden.input/txt}
echo "当前数据为：${inf%.*.*}" >> HOMO-LUMO.txt
grep -E 'HOMO|LUMO' ${inf//molden.input/txt}  >> HOMO-LUMO.txt
rm -f ${inf//molden.input/txt}
done
rm -f script.txt

批量输出HOMO-LUMO cube文件(Gaussian)

读取当前文件夹中所有的fchk文件，使用Multifn软件输出HOMO与LUMO的cube文件并命名。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#!/bin/bash
cat > script.txt << EOF 
200
3
h
3
3
l
3
EOF
for inf in *.fchk
do
multiwfn $inf < script.txt
name1=${inf%.*}"-homo.cube"
name2=${inf%.*}"-lumo.cube"
mv h.cub $name1
mv l.cub $name2
done
rm -f script.txt

批量输出HOMO-LUMO cube文件(Orca)

读取当前文件夹中所有的molden.input文件，使用Multifn软件输出HOMO与LUMO的cube文件并命名。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
#!/bin/bash
cat > script.txt << EOF 
200
3
h
3
3
l
3
EOF
for inf in *.molden.input
do
multiwfn $inf < script.txt 
name1=${inf%.*}"-homo.cube"
name2=${inf%.*}"-lumo.cube"
mv h.cub $name1
mv l.cub $name2
done
rm -f script.txt

pbs提交dalton任务

该脚本会自动识别当前文件夹中的dal文件，并根据out文件判断该dal文件是否被计算过；任务提交和计算完成时会在pbs.log文件中添加记录。需要根据实际安装情况修改orca的安装路径。使用pbs执行该脚本时会要求输入CPU核心数和每个核心所使用的最大内存数。

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
#!/bin/bash
for inf in *.dal
do
path=$(cd "$(dirname "$inf")";pwd)
if [ ! -f ${inf//dal/out} ];  ##判断dal文件是否存在对应的out文件
then
cat > pbs.sh << "EOF"
#PBS -N dalton
#PBS -q single
#PBS -j oe
source /opt/intel/oneapi/setvars.sh
cd $PBS_O_WORKDIR
INPUT_NAME=1    ##这里输入文件名，不含拓展名##
cpu=1       	##这里输入CPU核心数##
mem=500       	##这里输入内存##
dalton -N $cpu -mb $mem $INPUT_NAME $INPUT_NAME
wait
echo `date +"%Y-%m-%d %H:%M:%S"` job任务已经计算完成！ >> /home/pbs.log
EOF
name="INPUT_NAME="${inf%.*}
cpuinput="cpu="$1
meminput="mem="$2
filename="$path/$inf"
sed -i -r "s/(INPUT_NAME=)[^#]*/$name    /" pbs.sh
sed -i -r "s/(cpu=)[^#]*/$cpuinput    /" pbs.sh
sed -i -r "s/(mem=)[^#]*/$meminput    /" pbs.sh
sed -i -r "s?job?$filename?" pbs.sh   ##替换pbs.sh中文件路径
qsub pbs.sh
rm -f pbs.sh
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经提交至PBS队列！ >> /home/pbs.log
else
echo `date +"%Y-%m-%d %H:%M:%S"` ${path}/${inf}任务已经计算过了！ >> /home/pbs.log
fi
done

fchk-gjf批量转换(保留溶剂，保留赝势)

该脚本将读取fchk文件中的结构，并按照特定需要转换为gjf任务文件，例如结构优化完计算吸收光谱，溶剂与优化结构时选用的溶剂一致。

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
#!/bin/bash
cat > template.gjf <<\EOF 
%nprocshared=12
%mem=24GB
%chk=[name].chk
#p td=(triplets,nstates=3) mn15/genecp scrf=(solvent=acetonitrile) g09default

gaussian

1 1
[geometry]
EOF
for inf in *.fchk
do
oldname=${inf%.*}".gjf"
name=${inf%.*}"-ph.gjf"
cat > script.txt << EOF 
100
2
10
$name
n
q
EOF
multiwfn $inf < script.txt
rm -f script.txt
#使用原计算文件溶剂
oldstring='grep -o "scrf=(solvent=.*)" $name'
newstring='grep -o "scrf=(solvent=.*)" $oldname'
sed -i "s/$oldstring/$newstring/" $name
#通过genecp判断赝势，使用原计算文件赝势，原gjf文件空两行
if cat $name |grep 'genecp' > /dev/null
then
tail -n 11 $oldname >> $name
else
echo "" >> $name
fi
#修改第8行电荷与自旋多重度
sed -i '8c 1 1' $name
done
rm -f template.gjf

fchk-inp批量转换(Orca)

本脚本为bat批处理文件，将下列代码保存为bat文件并执行，将会读取当前文件夹中所有fchk文件，通过调用multiwfn（软件按照实际情况使用，或者直接使用绝对路径）输出Orca计算的输入文件，goto :eof之后的文字为Orca计算所使用的关键词。

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
del soc-multiwfn.txt
echo 100 >> soc-multiwfn.txt
echo 2 >> soc-multiwfn.txt
echo 12 >> soc-multiwfn.txt
echo. >> soc-multiwfn.txt
echo 1 >> soc-multiwfn.txt
for %%i in (*.fchk) do (
multiwfn %%i < soc-multiwfn.txt
rename %%~ni.inp %%~ni-soc.inp
)
del soc-multiwfn.txt
:: 注意把要加的内容写在第18(即代码中more +18的那个数)行之下
for %%i in (*.inp) do (
    more +18 "%~0" > "%%i.tmp"
    more +3 "%%i" >> "%%i.tmp"
    move /y "%%i.tmp" "%%i"
)
goto :eof
! B3LYP/G def2-svp miniprint tightSCF
%pal nprocs 24 end
%maxcore 2500
%tddft
nroots 20
dosoc true
tda false
printlevel 3
end

提取柔性扫描轨迹结构

参考链接：批量提取柔性扫描的能量及结构 (123qwertybobo)
如下处理的柔性扫描任务是计算改变D(2,1,53,40)二面角时的激发态能量，该bash脚本能够处理柔性扫描log文件得到轨迹文件和不同二面角是的激发态能量。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
#!/bin/bash
i=0
for lod in *.log
do
icc=0
loc=${lod%.log}
mkdir -p ${loc}
for i in `sed -n -e '/Stationary point found/=' ${loc}.log`
do
((icc++))
dir0=`head -n ${i} ${loc}.log | sed -n -e '/Standard orientation:/=' | tail -n 1`
dir2=`sed -n ''${dir0}','${i}'{/Excited State   1:      Triplet-A/p}' ${loc}.log | tail -n 1 | awk -F " " '{print $5}'`
dir3=`sed -n ''${dir0}','${i}'{/Excited State   1:      Triplet-A/p}' ${loc}.log | tail -n 1 | awk -F " " '{print $7}'`
dir1=`sed -n ''${i}','$((${i}+1000))'{/D(2,1,53,40)/p}' ${loc}.log | tail -n 1 | awk -F " " '{print $4}'`
printf {"%-4d %10.4f %10.4f eV %10.2f nm \n",${icc},${dir1},${dir2},${dir3}} >> ${loc}/${loc}-abs.txt
NATOM=`grep "NAtoms=" ${loc}.log | tail -1 | awk '{print $2}'`
echo -e ${NATOM} >> ${loc}/${loc}.xyz
echo -e ${icc} >> ${loc}/${loc}.xyz
sed -n ''$((${dir0}+5))','$((${dir0}+${NATOM}+4))'p' ${loc}.log | tail -$NATOM | awk '{gsub(17,"Cl",$2); gsub(16,"S",$2); gsub(6,"C",$2); gsub(7,"N",$2); gsub(8,"O",$2);gsub(15,"P",$2); gsub(33,"As",$2);gsub(51,"Sb",$2);gsub(16,"S",$2);gsub(34,"Se",$2); gsub(52,"Te",$2);gsub(52,"Te",$2);gsub(9,"F",$2);gsub(17,"Cl",$2);gsub(33,"Br",$2);gsub(53,"I",$2);gsub("1","H",$2);printf "%2s %12.7f %12.7f %12.7f\n", $2, $4, $5, $6}' >> ${loc}/${loc}.xyz
done
done

SOC数据后处理(Orca)

本脚本为python脚本，将下列代码保存并执行，将会读取当前文件夹中Orca计算旋轨耦合矩阵元的结果文件（out文件），并输出单线态和三线态之间的旋轨耦合矩阵元。脚本会输出两个文件：一个是计算好的SOC文件，第一和二列分别是三线态和单线态序号，第三列是SOC值；第二个文件是提取出来的原始CALCULATED REDUCED SOCME BETWEEN TRIPLETS AND SINGLETS数据。

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
#!/usr/bin/python
#coding=utf-8
import re
import sys
import ast
import math
import time
import os
import glob

def get_info(text):
  head='T      S           MS= 0                  -1                    \+1'
  tail='CALCULATED REDUCED SOCME BETWEEN TRIPLETS AND SINGLETS'
  data_content=re.findall(head+'(.*)'+tail,text,re.S)
  return data_content

start = time.time()
path_vec = glob.glob(r'*.out')
for k in range(0,len(path_vec)):
    print("源文件：" + path_vec[k])
    portion = os.path.splitext(path_vec[k])#将文件名拆成名字和后缀
    newname = os.path.join(portion[0] + '-SOC.out') 
    newname1 = os.path.join(portion[0] + '-SOC-raw.out') 
    print("数据汇总：" + newname)
    with open(path_vec[k], 'r') as f:
      text=f.read()
      content=get_info(text)
      str2 = ''.join(content)
      str1 = re.findall("[-+]?[0-9]*\.?[0-9]+",str2)
      str3=[]
      for x in range(0,len(str1),8):
          i=str1[x]
          j=str1[x+1]
          i = int (i)
          j = int (j)
          b = 0
          for c in range(6):
            a = str1[x+c+2]
            a =float (a)
            b = a*a+b
          g = math.sqrt(b)
          str4=[i,j,g]
          str3.append(str4)
      str3 = ''.join('%s' %id for id in str3)
      str3 = str3.replace('][','\n')
      str3 = str3.strip('[')
      str3 = str3.strip(']')
      with open(newname,"w") as f:
              f.write(str3)
      str2.split('\n')
      with open(newname1,"w") as f:
              f.write(str2+'\n')    
print ("第一列为三线态，第二列为单线态！")     
end = time.time()
time = end-start
print ("elapsed time(s):"+ str(time))