xtb-Molclus分子构象采样

分子模拟 > 构象搜索分子模拟构象搜索 xtb Molclus

发布时间 : 2025-01-26 01:22

阅读 :

基本流程

使用xtb跑分子动力学
使用Molclus调用xtb分别在gfn0和gfn2下优化结构
利用编写好的脚本使用Gaussian和ORCA计算分子的单点能
利用Molclus求得玻尔兹曼分布
使用得到的结构继续后续的操作

xtb分子动力学模拟

直接在xtb工作目录中进行，需要：
1. 分子的xyz结构文件
2. md.inp配置文件

md.inp配置文件

$md
   temp= 400  # in K
   time= 500.0  # in ps
   dump= 50.0  # in fs
   step=  1.0  # in fs
   hmass=1 # do not change mass of hydrogens
   shake=1 # constraint X-H bonds
$end
复制代码

一般只需要调整温度和时长即可，要根据具体的分子结构来调整参数。比如说有的分子在高温下容易分解，这时候可以延长模拟时间。在运行的过程中，将产生的xtb.trj文件复制出来，后缀改为xyz，使用VMD或者pymol查看分子动力学运行的情况。

一行命令运行xtb

如果要采样，不需要对用到的指令有大幅度的变动。唯一要注意的是分子的电荷。
使用下面的指令即可使用xtb进行分子动力学模拟：

xtb 分子名称.xyz --input md.inp --omd --gfn 0 --chrg 0
bash复制代码

Caution! :

将分子名称.xyz改为分子结构文件的名称
保证md.inp文件和你的xyz文件在同一目录下
--chrg 0一定要根据你的分子进行修改！！！否则在运行的时候分子有可能会碎掉，而且结果不可用。

xtb运行时可能会用到的脚本

要使用以下的脚本有几点需要注意：

脚本文件的编码格式，如果直接用记事本编辑有可能无法在Linux系统下正确识别运行
创建好脚本文件之后，运行chmod +x YourScript.sh为脚本文件添加运行权限

run.sh

#!/bin/bash

cd /home/chain/xtb_work/XDS1-noZn
xtb XDS1-noZn.xyz --input md.inp --omd --gfn 0 --chrg 0
echo "XDS1-noZn have been done"

cd /home/chain/xtb_work/XDS1-Zn
xtb XDS1-Zn.xyz --input md.inp --omd --gfn 0 --chrg 2
echo "XDS1-Zn have been done"

cd /home/chain/xtb_work
bash复制代码

当你有多个分子，并且每个分子的电荷不同时，可以使用该脚本。
一个分子指定一个目录，例如：
cd /home/chain/xtb_work/XDS1-noZn

下一行写具体的运行指令，例如：
xtb XDS1-noZn.xyz --input md.inp --omd --gfn 0 --chrg 0

echo行可写可不写

run.sh的运行

./run.sh
bash复制代码

runall.sh

遍历当前路径下的每个文件夹全都进行xtb计算

#!/bin/bash

# 遍历当前路径下的每个文件夹全都进行xtb计算
for dir in */; do
    # 进入文件夹
    cd "$dir" || continue
    
    # 检查是否存在 .xyz 文件
    if ls *.xyz 1> /dev/null 2>&1; then
        # 遍历每个 .xyz 文件
        for file in *.xyz; do
            echo "Processing $file in $dir..."
            # 运行 xtb 命令
            xtb "$file" --input ../md.inp --omd --gfn 0 --chrg 0
        done
    else
        echo "No .xyz files found in $dir."
    fi

    # 返回上一级目录
    cd ..
done
bash复制代码

当你有多个分子，并且每个分子的电荷相同时，可以使用该脚本。
该脚本将当前目录下的文件夹遍历一遍，并且执行其中的xtb计算任务。
默认每个分子的电荷为0，这在脚本中的xtb "$file" --input ../md.inp --omd --gfn 0 --chrg 0中有体现，修改--chrg参数即可修改计算使用的电荷。

runall.sh的运行

./runall.sh
bash复制代码

Molclus调用xtb进行结构优化

使用gfn 0进行粗优化

该过程要在Molclus目录下进行，需要遵循以下步骤：

将上一步xtb分子动力学产生的xtb.trj文件改为traj.xyz

cp xtb.trj traj.xyz
复制代码

将traj.xyz复制到Molclus目录下

cp traj.xyz /home/chain/molclus_1.12_Linux/traj.xyz
复制代码

修改setting.ini参数，有几点参数需要注意：
1. iprog= 4表示调用xtb计算。
2. ngeom= 0一般而言设置为0，但需要根据分子动力学具体结果调整。比如，分子动力学模拟后半程分子碎裂，那么为了节约时间，就可以将ngeom设置为感兴趣的范围。
3. itask= 0表示进行优化结构任务。
4. xtb_arg= "--gfn 0 --chrg 2 --uhf 0"，其中--chrg参数要根据分子电荷进行调整。
运行./molclus即可开始计算。
运行结束生成最终的isomers.xyz，其中包含着使用gfn 0优化过后的每一帧结构。
运行./isostat，对isomers.xyz中的结构去重，依次执行下面操作：（参数依然需要根据具体体系调整，比如柔性分子可以适当增大两个阈值）
1. 【回车】
2. 0.5（能量去重阈值）
3. 0.5（结构去重阈值）
4. 【回车】
运行结束生成cluster.xyz文件，其中包含去重之后的结构。

使用gfn 2进一步优化

同样在Molclus目录下进行：

将之前的traj.xyz文件删除

rm traj.xyz
复制代码

将cluster.xyz复制为traj.xyz

cp cluster.xyz traj.xyz
复制代码

修改setting.ini文件参数，除了xtb_arg，其余同上：
1. xtb_arg= "--gfn 2 --gbsa h2o --chrg 0 --uhf 0"
其余步骤同[[#使用gfn 0进行粗优化]]。
最后得到cluster.xyz文件。

Gaussian-ORCA联用优化并计算单点能（该步骤可根据需求跳到[[#不使用Molclus调用Gaussian-ORCA计算]]）

删除traj.xyz文件

rm traj.xyz
复制代码

将cluster.xyz改为traj.xyz

mv cluster.xyz traj.xyz
复制代码

修改setting.ini文件参数：
1. iprog= 1，表示调用Gaussian。
2. ngeom= n，表示计算traj.xyz前n个构象。
3. itask= 3，表示任务类型是优化+振动分析获得自由能
检查template.gjf文件参数，此文件用于自动生成Gaussian输入文件，进行结构优化。
检查template_SP.inp文件参数，此文件用于自动生成ORCA输入文件，进行单点能计算。
运行./molclus | tee out.txt，可以在out.txt文件中查看屏幕上输出的信息。
运行结束生成isomers.xyz，可以查看其中每个结构的自由能，虚频数
运行./isostat，对isomers.xyz文件中的结构进行去重
1. 【回车】
2. 0.5（能量去重阈值）
3. 0.5（结构去重阈值）
4. 298.15
屏幕打印298.15K下，各结构的玻尔兹曼分布情况。根据结果选取结构即可。同时生成cluster.xyz文件。

一键脚本

gfn0-gfn2连算脚本

经xtb分子动力学模拟的文件traj.xyz置于Molclus目录下，且setting.ini配置文件中的iprog, ngeom, itask均已设置好。
运行以下脚本：

#!/bin/bash

#脚本参数
#molclus目录
molclus=/home/chain/molclus_1.12_Linux
#gfn 0计算参数
gfn0_arg="--gfn 0 --chrg 0 --uhf 0"
#gfn 2计算参数
gfn2_arg="--gfn 2 --gbsa h2o --chrg 0 --uhf 0"
#创建结果存储目录
results=/home/chain/molclus_1.12_Linux/gfn_results
#能量阈值
gfn0_en=0.5
gfn2_en=0.5
#结构阈值
gfn0_st=0.5
gfn2_st=0.5

#开始计时
start=$(date)
start_time=$(date +%s)

#初始化setting.ini
sed -i 's/xtb_arg= "[^"]*"/xtb_arg= "$gfn0_arg"/' "$molclus/setting.ini"

#进入molclus目录
cd $molclus

#gfn 0结构优化
./molclus
echo -e "\n$gfn0_en\n$gfn0_st\n\n" | ./isostat

#备份gfn 0计算文件
mv traj.xyz traj.xyz.gfn0.bak
mv isomers.xyz isomers.xyz.gfn0.bak
cp cluster.xyz cluster.xyz.gfn0.bak
mv cluster.xyz traj.xyz

#修改setting.ini中配置
sed -i 's/xtb_arg= "$gfn0_arg"/xtb_arg= "$gfn2_arg"/' "setting.ini"
echo "配置文件已更新。"

#gfn 2结构优化
./molclus
echo -e "\n$gfn2_en\n$gfn2_st\n\n" | ./isostat

#备份gfn 2计算文件
mv traj.xyz traj.xyz.gfn2.bak
mv isomers.xyz isomers.xyz.gfn2.bak
mv cluster.xyz cluster.xyz.gfn2.bak

#存储结果
mkdir -p $results
mv traj.xyz.gfn0.bak $results/traj.xyz.gfn0.bak
mv isomers.xyz.gfn0.bak $results/isomers.xyz.gfn0.bak
mv cluster.xyz.gfn0.bak $results/cluster.xyz.gfn0.bak
mv traj.xyz.gfn2.bak $results/traj.xyz.gfn2.bak
mv isomers.xyz.gfn2.bak $results/isomers.xyz.gfn2.bak
mv cluster.xyz.gfn2.bak $results/cluster.xyz.gfn2.bak

#结束计时
end=$(date)
end_time=$(date +%s)

#计算耗时
elapsed_time=$((end_time - start_time))

echo "gfn 0-gfn 2结构优化任务已完成。"
echo "开始时间：$start"
echo "结束时间：$end"
echo "用时：$elapsed_time 秒"
bash复制代码

xtb分子动力学-gfn优化连算脚本

将分子xyz结构文件命名为traj.xyz置于xtb_work文件夹下，保证同一目录下有md.inp文件。

使用前需要确认md.inp文件中参数、分子电荷

#!/bin/bash
set -e

##脚本参数
#molclus目录
molclus=/home/chain/molclus_1.12_Linux
#分子电荷
charge=0
#gfn 0计算参数
gfn0_arg="--gfn 0 --chrg 0 --uhf 0"
#gfn 2计算参数
gfn2_arg="--gfn 2 --gbsa h2o --chrg 0 --uhf 0"
#创建结果存储目录
results=/home/chain/molclus_1.12_Linux/xtb_gfn_results
#能量阈值
gfn0_en=0.5
gfn2_en=0.5
#结构阈值
gfn0_st=0.5
gfn2_st=0.5

#开始计时
start_0=$(date)
start_time=$(date +%s)

#调用xtb
xtb traj.xyz --input md.inp --omd --gfn 0 --chrg $charge
#复制文件至molclus目录下
cp xtb.trj $molclus/traj.xyz

#初始化setting.ini
sed -i 's/xtb_arg= "[^"]*"/xtb_arg= "$gfn0_arg"/' "$molclus/setting.ini"

#进入molclus目录
cd $molclus

#gfn 0结构优化
./molclus
echo -e "\n$gfn0_en\n$gfn0_st\n\n" | ./isostat

#备份gfn 0计算文件
mv traj.xyz traj.xyz.gfn0.bak
mv isomers.xyz isomers.xyz.gfn0.bak
cp cluster.xyz cluster.xyz.gfn0.bak
mv cluster.xyz traj.xyz

#修改setting.ini中配置
sed -i 's/xtb_arg= "$gfn0_arg"/xtb_arg= "$gfn2_arg"/' "setting.ini"
echo "配置文件已更新。"

#gfn 2结构优化
./molclus
echo -e "\n$gfn2_en\n$gfn2_st\n\n" | ./isostat

#备份gfn 2计算文件
mv traj.xyz traj.xyz.gfn2.bak
mv isomers.xyz isomers.xyz.gfn2.bak
mv cluster.xyz cluster.xyz.gfn2.bak

#存储结果
mkdir -p $results
mv traj.xyz.gfn0.bak $results/traj.xyz.gfn0.bak
mv isomers.xyz.gfn0.bak $results/isomers.xyz.gfn0.bak
mv cluster.xyz.gfn0.bak $results/cluster.xyz.gfn0.bak
mv traj.xyz.gfn2.bak $results/traj.xyz.gfn2.bak
mv isomers.xyz.gfn2.bak $results/isomers.xyz.gfn2.bak
mv cluster.xyz.gfn2.bak $results/cluster.xyz.gfn2.bak

#结束计时
end_0=$(date)
end_time=$(date +%s)

#计算耗时
elapsed_time=$((end_time - start_time))

echo "xtb-gfn 0-gfn 2任务已完成。"
echo "开始时间：$start_0"
echo "结束时间：$end_0"
echo "用时：$elapsed_time 秒"
bash复制代码

不使用Molclus调用Gaussian-ORCA计算

操作流程

将gfn 2优化去重产生的cluster.xyz文件改名为traj.xyz。
将traj.xyz放进traj2gjf文件夹下。
配置好template.gjf模板文件，稍后生成的输入文件全都是基于此配置文件生成的。
运行./traj2gjf.sh [arguments]，其中[arguments]为感兴趣的结构范围。例如：

./traj2gjf.sh 1-3 5 7-9
复制代码

表示生成traj.xyz文件中1-3 5 7-9的结构的输入文件。
5. 生成gau*.gjf文件，例如：gau1.gjf。可将gau*.gjf文件放在服务器上运行。
6. 运行产生gau*.out文件。
7. 将gau*.out文件全部放入out2orca_SP文件夹下。
8. 配置好template_SP.inp模板文件，稍后生成的输入文件全都是基于此配置文件生成的。
9. 运行./gau2orca_SP.sh。
10. 生成gau*.xyz文件和orca*.inp文件，可将orca*.inp文件放在服务器上运行。
11. 运行产生orca*.out文件。
12. 将gau*.xyz, gau*.out, orca*.out文件全部放入gene_isomers文件夹下。
13. 运行./gene_isomers，生成isomers.xyz文件。
14. 将isomers.xyz文件放入Molclus目录下
15. 运行./isostat，对isomers.xyz文件中的结构进行去重
1. 【回车】
2. 0.5（能量去重阈值）
3. 0.5（结构去重阈值）
4. 298.15
16. 屏幕打印298.15K下，各结构的玻尔兹曼分布情况。根据结果选取结构即可。同时生成cluster.xyz文件。

使用到的脚本及使用注意事项

脚本文件夹配置

需要脚本文件（将三个脚本分别放入不同文件夹）：

traj2gjf，生成Gaussian计算输入文件。所需文件如下：
1. traj2gjf.sh
2. template.gjf
3. traj.xyz
out2orca_SP，利用Gaussian输出文件生成ORCA单点能计算输入文件。所需文件如下：
1. gau2orca_SP.sh
2. template_SP.inp
3. gau*.out
gene_isomers，生成isomers.xyz文件，用于结构去重找到最后的稳定结构。所需文件如下：
1. gene_isomers.sh
2. gau*.out
3. gau*.xyz
4. orca*.out

traj2gjf.sh

#!/bin/bash

# Script: traj2gjf.sh
# Description: Extracts geometries from traj.xyz and generates .gjf files based on template.gjf
# Usage: ./traj2gjf.sh 1-3 5 7-9

# Check if the required files exist
if [[ ! -f "traj.xyz" ]]; then
    echo "Error: traj.xyz file not found."
    exit 1
fi

if [[ ! -f "template.gjf" ]]; then
    echo "Error: template.gjf file not found."
    exit 1
fi

# Function to expand ranges into a list of numbers
expand_range() {
    local input="$1"
    local result=()
    IFS=',' read -ra tokens <<< "$input"
    for token in "${tokens[@]}"; do
        if [[ $token =~ ^[0-9]+$ ]]; then
            result+=("$token")
        elif [[ $token =~ ^([0-9]+)-([0-9]+)$ ]]; then
            start=${BASH_REMATCH[1]}
            end=${BASH_REMATCH[2]}
            start=${token%-*}
            end=${token#*-}
            if ! [[ $start =~ ^[0-9]+$ && $end =~ ^[0-9]+$ ]]; then
                echo "Invalid range: $token"
                exit 1
            fi
            if (( start > end )); then
                echo "Invalid range (start > end): $token"
                exit 1
            fi
            for ((i=start; i<=end; i++)); do
                result+=("$i")
            done
        else
            echo "Invalid token: $token"
            exit 1
        fi
    done
    echo "${result[@]}"
}

# Parse input ranges to generate a list of cluster numbers
clusters=()
for arg in "$@"; do
    cluster_numbers=$(expand_range "$arg")
    clusters+=($cluster_numbers)
done

# Remove duplicates and sort the cluster numbers
clusters=($(echo "${clusters[@]}" | tr ' ' '\n' | sort -n | uniq))

if [[ ${#clusters[@]} -eq 0 ]]; then
    echo "No valid cluster numbers provided."
    exit 1
fi

# Read template.gjf content
template_content=$(cat template.gjf)
if [[ -z "$template_content" ]]; then
    echo "Error: template.gjf is empty."
    exit 1
fi

# Initialize variables
current_cluster=0
read_geometry=false
geometry_lines=()
cluster_number=0

# Read traj.xyz and process clusters
line_number=0

while IFS= read -r line || [[ -n "$line" ]]; do
    ((line_number++))

    # Trim leading and trailing whitespace
    trimmed_line="${line#"${line%%[![:space:]]*}"}"
    trimmed_line="${trimmed_line%"${trimmed_line##*[![:space:]]}"}"

    # Skip empty lines
    if [[ -z "$trimmed_line" ]]; then
        continue
    fi

    # Check if the line indicates the number of atoms
    if [[ $trimmed_line =~ ^([0-9]+)$ ]]; then
        num_atoms=${BASH_REMATCH[1]}
        read_geometry=false
        geometry_lines=()
        continue
    fi

    # Check if the line contains the cluster number
    if [[ $trimmed_line =~ ^Energy.*#Cluster:[[:space:]]*([0-9]+) ]]; then
        cluster_number=${BASH_REMATCH[1]}
        if [[ " ${clusters[@]} " =~ " ${cluster_number} " ]]; then
            read_geometry=true
            atom_count=0
        else
            read_geometry=false
        fi
        continue
    fi

    # Collect geometry lines
    if $read_geometry; then
        if [[ -n "$trimmed_line" ]]; then
            geometry_lines+=("$trimmed_line")
            ((atom_count++))
        fi
        if [[ $atom_count -eq $num_atoms ]]; then
            # Replace [GEOMETRY] in template and write to output file
            geometry_text=$(printf "%s\n" "${geometry_lines[@]}")
            output_content="${template_content/\[GEOMETRY\]/$geometry_text}"
            output_content="${output_content/\%chk=file.chk/\%chk=gau${cluster_number}.chk}"
            output_file="gau${cluster_number}.gjf"
            echo "$output_content" > "$output_file"
            echo "Generated $output_file"
            read_geometry=false
        fi
    fi

done < "traj.xyz"

# Check if any clusters were generated
if [[ ! -f "gau${clusters[0]}.gjf" ]]; then
    echo "No clusters were generated. Please check if the cluster numbers exist in traj.xyz."
    exit 1
fi

bash复制代码

out2orca_SP.sh

使用该脚本前一定要安装好Multiwfn，并配置好环境变量。由于我使用的是WSL，所以在这里用的是Multiwfn_noGUI。

#converte Gaussian ouput file to ORCA inp file based on template_SP.inp. (Modified from out2gjf.sh of Sobereva)
#!/bin/bash
icc=0
nfile=`ls *.out|wc -l`
for inf in *.out
do
((icc++))
echo Converting ${inf} to ${inf//out/xyz} ... \($icc of $nfile\)
Multiwfn_noGUI ${inf} << EOF > /dev/null
100
2
2
${inf//out/xyz}
0
q
EOF
done
# Now, it has converted the all .out file to .xyz file in this directory. 
# Then, we need to extract the geometry from .xyz file and replace into template_SP.inp
icc=0
nfile=`ls *.out|wc -l`
template=$(cat template_SP.inp)
for xyz in *.xyz
do
((icc++))
output_file1=${xyz//gau/orca}
output_file=${output_file1//xyz/inp}
echo Converting ${xyz} to ${output_file} \($icc of $nfile\)
geometry=$(tail -n +3 ${xyz})
output_content=${template//\[GEOMETRY\]/${geometry}}
echo "$output_content" > "$output_file"
done
bash复制代码

gene_isomers.sh

#!/bin/bash
icc=0
nfile=`ls gau* | wc -l`
output_file="isomers.xyz"
for inf in gau*.out
do
    ((icc++))
    gau_file="gau${icc}.out"
    orca_file="orca${icc}.out"
    xyz_file="gau${icc}.xyz"
    gcorr=$(grep "Thermal correction to Gibbs Free Energy=" ${gau_file} | awk '{print $NF}')
    single_point=$(grep "FINAL SINGLE POINT ENERGY" ${orca_file} | awk '{print $NF}')
    gibbs=$(echo "$gcorr + $single_point" | bc)
    NImag=$(grep "Frequencies --" ${gau_file} | awk '{for (i=3; i<=NF; i++) if ($i < 0) count++} END {print count+0}')
    NAtoms=$(grep "NAtoms=" ${gau_file} | awk 'NR==1 {print $2}')
    geometry=$(tail -n +3 ${xyz_file})
    echo "${NAtoms}" >> ${output_file}
    echo "G = ${gibbs} Gcorr = ${gcorr} a.u. Nimag = ${NImag} Frame: ${icc} Step: ${icc} Geom: ${icc}" >> ${output_file}
    echo "${geometry}" >> ${output_file}
done
bash复制代码

template.gjf

%nprocshared=64
%mem=150GB
# B3LYP/def2svp em=GD3BJ opt freq int=fine

Template file

0 1
[GEOMETRY]


复制代码

template_SP.inp

! wB97M-V def2-TZVP def2/J RIJCOSX strongSCF noautostart miniprint nopop
%maxcore     950
%pal nprocs   8 end
%cpcm
smd true
SMDsolvent "water"
end
* xyz   0   1
[GEOMETRY]
 *
复制代码

转载请注明来源，欢迎对文章中的引用来源进行考证，欢迎指出任何有错误或不够清晰的表达。可以在下面评论区评论，也可以邮件至 1971945941@qq.com