How to get rid of part of molecules in one snapshot becuase of PBC

Dear Lammps users,

I want to minimize a box of water and then after minimization, I want to load the last snapshot of dump file in VMD and save that as a pdb file. Then I want to use this pdb to build another system. I want to combine this minimized water box with other components.
Now the problem is that when I load the last frame in VMD, there is H and O-H in surface because of PBC. They make problem for my system when I create data file where there is extra wrong bonds and angles in the data file.
So I tried to build a larger box of water (larger than the box that I need) and after minimization, I used delete_atom command to remove the surface in all the 3 directions.
In this way I could get rid of the H and O-H in surface but another problem appeared.
The box now with the below command does not have time to be adapted with the new size and the density of water is wrong. So I need to correct the size manually to have a correct density. This reduces the size in for example z direction (if x and y be constant) as much as 10 A. But with this size in z direction, Lammps gives me an error: atom counts is inconsistent … and in my case it is because of smaller size. I need to increase size in z direction little by little to remove that error. If I increase the size more than a number, then I will have bubble in water since its density becomes less that 0.997 g/cm^3.
Also another problem is that I want that the pressure be 1 bar. This is OK before using delete-atom command since after that the pressure increases until 800 bar.
I appreciate any comment to solve the problem from the first step. I mean how can I get rid of those H and O-H parts in one snapshot when dump file is loaded in VMD to be saved as pdb, if I do not use delete_atom command?

Also I see WARNING: Ignoring ‘compress yes’ for molecular system (src/delete_atoms.cpp:125). Does it make any problem in my results?
Thank you.

input file:


fix 1 all nvt temp 298.15 298.15 0.1
dump mydump all xyz 500 dump.xyz
dump_modify mydump element H O
run 80000

region extra block -63.1245 63.1245 -4.4856 4.4856 -20.000 20.000 side out
delete_atoms region extra mol yes bond yes
dump mydump1 all xyz 1 dump1.xyz
dump_modify mydump1 element H O
run 0

log file:


Step Temp Density Press Volume
70000 296.29871 0.98423631 -695.348 59907.364
71000 295.5356 0.98423631 -660.56891 59907.364
72000 299.87829 0.98423631 403.43043 59907.364
73000 300.86958 0.98423631 348.74989 59907.364
74000 300.31389 0.98423631 -40.800278 59907.364
75000 296.33583 0.98423631 61.091599 59907.364
76000 293.50274 0.98423631 -152.90919 59907.364
77000 292.23952 0.98423631 -61.899885 59907.364
78000 287.06727 0.98423631 -150.65725 59907.364
79000 296.18986 0.98423631 -155.34794 59907.364
80000 295.18662 0.98423631 286.77278 59907.364


delete_atoms region extra mol yes bond yes
WARNING: Ignoring ‘compress yes’ for molecular system (src/delete_atoms.cpp:125)
Deleted 2025 atoms, new total = 3888
Deleted 1350 bonds, new total = 2592
Deleted 675 angles, new total = 1296
thermo_style custom step temp density press vol
thermo 1
dump mydumpr all xyz 1 dumpwaterr.xyz
dump_modify mydumpr element H O
run 0


Step Temp Density Press Volume
80000 298.2442 0.64716908 800.61142 59907.364

The main issue at hand is that you are writing out “wrapped” coordinate data, if you would be writing unwrapped coordinates there would be no issue.
However, that is in part due to using unsuitable file formats like .xyz and .pdb.
.xyz files store no information about the simulation box, topology or much else, just coordinates and atom types. .pdb files are not much better, but they have a major deficiency in writing out coordinates in extremely low precision (not a problem for their original purpose of storing crystallographic data, which has very low precision).

So for post-processing system data to set up future simulations it would be much better to preserve as much information about the system as you can. That could be done, for example, by writing out a data file with the write_data command and then loading that into VMD with the TopoTools readlammpsdata subcommand. That will read and parse not just the principal cell coordinates but also the image flags and thus reconstruct the unwrapped coordinates. Topology data is read and stored as well. See the TopoTools documentation for more info.

Everything else you describe and ask about are just consequences of your bad initial choices and thus not worth discussing.

Thank you for your comments.
I have VMD 1.9.4, on Windows, and when I typed package require topotools in TKConsole, it returned 1.8 while in this website Homepage of Axel Kohlmeyer - TopoTools I see only 1.7 as the last version. Does 1.7 work on all new versions of VMD after 1.9.3?

Right now I found that there is already topotools1.8 in VMD 1.9.4 installation.

The website is outdated. Just use the version bundled with VMD.