如何正确筛除amber拓扑和轨迹文件中的水分子？

yorksy

各位老师好，我在用amber分析一次md的结果。因为这次模拟的水分子很多，我想用strip去除所有在我的配体7埃外的水分子。我在parmed和cpptraj都用了strip，reference是模拟过程中的第一个restart文件。输出的拓扑（.parm）和轨迹文件(.mdcrd)能成功地在vmd里面可视化出来。但是我看到所有水分子都是散落在模拟盒子各处的。我检查了一部分水分子跟配体的距离，发现在某一个时间点它们跟配体很近，大概三到四埃的样子。但是其他时间离配体都很远。请问这个问题我该怎么解决？
以下是我cpptraj的代码：

parm /home/data/chongyao/6pru_hip523_4/trajectory/0_initial.parm7
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0000.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0100.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0200.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0300.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0400.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0500.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0600.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0700.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0800.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_0900.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_1000.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_1100.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_1200.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_1300.nc 1 last 100
trajin /home/data/chongyao/6pru_hip523_4/trajectory/3_prod_npt_1400.nc 1 last 100
autoimagei
strip !(:PVB<:7.0)&:WAT
trajout output_water7A.mdcrd mdcrd

Serious

1）我觉得一帧一帧单独trajout比较好；
2）假定reference是3_prot_npt_0000.rst，strip是依照reference去找保留的分子。比如第100位水分子在reference中是保留的，符合strip距离的设定，那么trajin的所有帧都保留第100位水分子，所以会有“发现在某一个时间点它们跟配体很近，大概三到四埃的样子。但是其他时间离配体都很远”的情况；总的来说，strip的分子是按照reference定的；

yorksy

Serious 发表于 2024-2-12 23:57
1）我觉得一帧一帧单独trajout比较好；
2）假定reference是3_prot_npt_0000.rst，strip是依照reference去 ...

我也试过每一帧单独trajout，但这样就导致我不知道用什么拓扑结构。我的拓扑文件里面把所有水分子的编号都写上了。如果不能导出对应的拓扑文件的话，就会导致我的拓扑和结构合不到一起去，VMD里面显示一团乱。

Serious

yorksy 发表于 2024-2-14 17:11
我也试过每一帧单独trajout，但这样就导致我不知道用什么拓扑结构。我的拓扑文件里面把所有水分子的编号 ...

1）如果要导出相应top文件，可以试试cpptraj里parmwrite 和 parmstrip；
2）不过这样子挺麻烦的。。如果只是看结构，导pdb挺好；VMD没用过，不知论坛里的同志有没有办法。。