【已解决】Gromacs模拟水溶液蒸发的脚本

Eianghuan

请教各位老师，我在用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中的原子?
谢谢各位老师了

Eianghuan

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

———————————————————————————————————————分隔线——————————————————————————————————————
#!/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的数据，使用前请确保电脑有足够的空间。

Eianghuan

感谢各位老师的帮助

Eianghuan

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

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

pal

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

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

Eianghuan

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

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

Eianghuan

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

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

pal

在select的时候选择想要保留的原子不就好了

SharkYYX2025

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