计算化学公社

标题: 【已解决】Gromacs模拟水溶液蒸发的脚本 [打印本页]

作者
Author: Eianghuan 时间: 2025-12-12 22:26
标题: 【已解决】Gromacs模拟水溶液蒸发的脚本
本帖最后由 Eianghuan 于 2025-12-16 17:14 编辑

请教各位老师，我在用gromacs 2023.5做水溶液中水分子蒸发，溶质形成薄膜的模拟，我构建了一个10*10*20的盒子，溶液相位于下侧，真空相位于上侧。用nvt系综进行1ns模拟后，通过select命令：
gmx select -f evap_1.gro -s evap_1.tpr -on index_delMOL.ndx -select "resname SOL and z > 10.5 and z < 14.5"
生成了含有蒸发出去的水分子的信息的index文件，接下来要把这部分水分子从gro中删除。我尝试了命令：
gmx trjconv -f evap_1.gro -s evap_1.tpr -n index_delMOL.ndx -o evap_1_modified.gro
但这个命令是只留下index中记录的原子，请问我应该用什么指令从gro中删除index中的原子?
谢谢各位老师了

作者
Author: SharkYYX2025 时间: 2025-12-13 01:31
本帖最后由 SharkYYX2025 于 2025-12-13 01:41 编辑

你好，也可以考虑用VMD删除。可以按照这个类似格式的语句，在这个REP里面选择你想要的水分子，然后保存坐标，生成gro，再生成新的index。数字我随便输的，可以换成你需要的
(, 下载次数 Times of downloads: 0)

作者
Author: pal 时间: 2025-12-13 09:17
在select的时候选择想要保留的原子不就好了

作者
Author: Eianghuan 时间: 2025-12-13 10:30

SharkYYX2025 发表于 2025-12-13 01:31
你好，也可以考虑用VMD删除。可以按照这个类似格式的语句，在这个REP里面选择你想要的水分子，然后保存坐标 ...

好的，谢谢您，我需要用脚本反复删除续跑删除，所以VMD好像不太行

作者
Author: Eianghuan 时间: 2025-12-13 10:32

pal 发表于 2025-12-13 09:17
在select的时候选择想要保留的原子不就好了

谢谢您，我现在也是这么想的，但由于后面要从top文件中减少这部分被删除的水分子，用select选择要保留的原子会让top文件的自动更新计算变麻烦，不过应该能行，我试一下

作者
Author: pal 时间: 2025-12-13 11:00

Eianghuan 发表于 2025-12-13 10:32
谢谢您，我现在也是这么想的，但由于后面要从top文件中减少这部分被删除的水分子，用select选择要保留的 ...

这也不是啥问题，顶多多做一次select，两个select都做一次就好了

作者
Author: Eianghuan 时间: 2025-12-13 12:07

pal 发表于 2025-12-13 11:00
这也不是啥问题，顶多多做一次select，两个select都做一次就好了

好嘞，很有用，谢谢您的帮助

作者
Author: Eianghuan 时间: 2025-12-16 17:08
感谢公社的各位老师，模拟已经能够正常运行了，在这里分享一下我的脚本，后续有做相同工作的同学可以参考一下，希望各位老师如果发现有什么错误可以不吝指出：

———————————————————————————————————————分隔线——————————————————————————————————————
#!/bin/bash

# 初始结构（请修改为自己的初始输入文件的名称）
initial_structure="nvt.gro"
# 初始拓扑文件（请修改为自己的top文件名称）
initial_topology="topol.top"
# 循环次数（设置一个合适的循环次数）
num_cycles=8000
# 液相区阈值 (nm) - 高于此值的水分子去除，低于此值的所有分子保留（请按照自己的盒子大小设置一个合适的值，建议高于液相界面1~2 nm）
liquid_threshold=7

# GPU设置（按自己的硬件合理配置）
GPU_OPTIONS="-pin on -pinoffset 0 -ntmpi 1 -ntomp 24 -nb gpu -gpu_id 0"

# 检测从哪个循环开始（选择1则从第一个循环开始跑，若中间运行中断，可从最后一个循环重新开始续跑）
start_cycle=1
current_structure="$initial_structure"

# 检查是否存在已完成的循环
for ((i=1; i<=num_cycles; i++))
do
if [[ -f "evap_${i}_modified.gro" ]]; then
      echo "Found completed cycle $i"
      start_cycle=$((i + 1))
      current_structure="evap_${i}_modified.gro"
else
      break
fi
done

echo "Starting from cycle $start_cycle with structure: $current_structure"

# 如果起始循环大于总循环数，说明已经完成
if [[ $start_cycle -gt $num_cycles ]]; then
echo "All cycles already completed!"
exit 0
fi

# 复制初始拓扑文件
cp "$initial_topology" "topol_${start_cycle}.top"

for ((i=start_cycle; i<=num_cycles; i++))
do
echo "Starting evaporation cycle $i"

# 检查当前结构文件是否存在
if [[ ! -f "$current_structure" ]]; then
      echo "ERROR: Structure file $current_structure not found!"
      exit 1
fi

# 检查当前拓扑文件是否存在
current_topology="topol_${i}.top"
if [[ ! -f "$current_topology" ]]; then
      echo "ERROR: Topology file $current_topology not found!"
      exit 1
fi

# 1. 短时间MD模拟 (NVT) - 使用GPU加速
echo "Running MD simulation for cycle $i..."
gmx grompp -f md_evap.mdp -c "$current_structure" -p "$current_topology" -o "evap_${i}.tpr"
gmx mdrun -v -deffnm "evap_${i}" $GPU_OPTIONS

# 检查MD模拟是否成功完成
if [[ ! -f "evap_${i}.gro" ]]; then
      echo "ERROR: MD simulation failed to produce evap_${i}.gro"
      exit 1
fi

# 2. 选择要保留的分子（水和溶质DMAPS）
echo "Creating selection for molecules to keep (resname SOL or MOL in Z-range)..."
gmx select -f "evap_${i}.gro" -s "evap_${i}.tpr" -on "index_keepMOL_${i}.ndx" \
      -select "whole_mol_com z < $liquid_threshold"

# 3. 保留所需分子，输出新的gro坐标文件
echo "Generating modified structure with only selected molecules..."
gmx trjconv -f "evap_${i}.gro" -s "evap_${i}.tpr" -n "index_keepMOL_${i}.ndx" \
      -o "evap_${i}_modified.gro"

# 4. 统计新结构中的水分子数量
new_sol_count=$(grep "OW" "evap_${i}_modified.gro" | wc -l)
echo "New water molecules after deletion: $new_sol_count"

# 5. 为下一轮循环准备拓扑文件
next_topology="topol_$((i+1)).top"
cp "$current_topology" "$next_topology"

# 6. 更新拓扑文件中的SOL数量
sed -i "s/^SOL.*/SOL $new_sol_count/" "$next_topology"
echo "Updated $next_topology with new SOL count: $new_sol_count"

# 7. 为下一轮循环更新当前结构
current_structure="evap_${i}_modified.gro"

# 清理临时文件
rm -f "index_keepMOL_${i}.ndx"

echo "Completed cycle $i successfully"

# 保存检查点
echo $i > last_completed_cycle.txt
echo $current_structure > current_structure.txt

done

echo "All evaporation cycles completed successfully!"
#准备存储输出文件的文件夹
mkdir -p output
mv evap_* ./output
mv topol_* ./output
echo "Moved output files to directory: output"
——————————————————————————————————分隔线——————————————————————————————————————————

脚本用到的mdp文件如下：
define =
integrator = md
dt       = 0.002 ; ps
nsteps    = 100000 ; 0.2ns
comm-grps  = system
energygrps =
;
nstxout = 0
nstvout = 0
nstfout = 0
nstlog  = 1000
nstenergy = 1000
nstxout-compressed = 1000
compressed-x-grps  = system
;
pbc = xyz
cutoff-scheme = Verlet
coulombtype = PME
rcoulomb    = 1.0
vdwtype    = cut-off
rvdw       = 1.0
DispCorr    = no
;
Tcoupl  = V-rescale
tau_t = 0.2
tc_grps = system
ref_t = 298.15
;
Pcoupl    = no
pcoupltype = isotropic
tau_p = 2.0
ref_p = 1.0
compressibility = 4.5e-5
;
gen_vel  = no
gen_temp = 298.15
gen_seed = -1
;
freezegrps  =
freezedim =
constraints = hbonds
—————————————————————————————————————分隔线———————————————————————————————————————

经过我的个人测试，每个循环的运行时间对于蒸发走的水分子数量影响不大，50循环/每循环0.1ns 一共蒸发走20个水分子。50循环/每循环1ns一共蒸发走18个水分子。如果想蒸发走更多的水分子，应该增加循环次数，或者推测提高模拟温度也可以，但因为我的体系是室温蒸发条件，因此没有测试，有需要的同学可以测试一下试试。

最后需要指出的是，室温下的蒸发是一个比较缓慢的过程，因此这个脚本模拟下来占用的磁盘空间较大，2w个原子的盒子模拟8000次0.2ns的循环大约产生了300G的数据，使用前请确保电脑有足够的空间。

作者
Author: Eianghuan 时间: 2025-12-16 17:12
标题: 【已解决】Gromacs模拟水溶液蒸发的脚本
感谢各位老师的帮助

欢迎光临计算化学公社 (http://bbs.keinsci.com/)