Editing velocity files

[Saleh-OM4R]

Dear support team,

Is there a way to edit output.vel files to remove a ligand for example? I tried to load the pdb file then read the output.vel so I could edit it but it did not work and I got the following error. p.s. I ran the simulation using CHARMM36m forcefield,acemd3 MDengine, and htmd version 2.2.0. I would really appreciate your help on this issue.

v = Molecule('structure.pdb')
v.read('../output.vel', type='coor')
v.filter('not segid P1 P3')
v.write('x_out_test.vel', type='coor')

2023-09-12 16:44:36,378 - moleculekit.readers - WARNING - Non-integer values were read from the PDB "serial" field. Dropping PDB values and assigning new ones.
2023-09-12 16:44:36,452 - moleculekit.readers - WARNING - Reading PDB file with more than 99999 atoms. Bond information can be wrong.
---------------------------------------------------------------------------
MemoryError                               Traceback (most recent call last)
Cell In [35], line 3
      1 v = Molecule('structure.pdb')
      2 v.read('../output.vel', type='coor')
----> 3 v.filter('not segid P1 P3')
      4 v.write('x_out_test.vel', type='coor')

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/molecule.py:871, in Molecule.filter(self, sel, _logger)
    858 def filter(self, sel, _logger=True):
    859     """Removes all atoms not included in the selection
    860 
    861     Parameters
   (...)
    869     >>> mol.filter('protein')
    870     """
--> 871     s = self.atomselect(sel)
    872     if np.all(s):  # If all are selected do nothing
    873         return np.array([], dtype=np.int32)

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/molecule.py:825, in Molecule.atomselect(self, sel, indexes, strict, fileBonds, guessBonds, _debug)
    820     s = np.ones(self.numAtoms, dtype=bool)
    821 elif isinstance(sel, str):
    822     s = atomselect(
    823         self,
    824         sel,
--> 825         bonds=self._getBonds(fileBonds, guessBonds),
    826         _return_ast=_debug,
    827         _debug=_debug,
    828     )
    829     if _debug:
    830         s, ast = s

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/molecule.py:778, in Molecule._getBonds(self, fileBonds, guessBonds)
    776     bonds = np.vstack((bonds, self.bonds))
    777 if guessBonds:
--> 778     bonds = np.vstack((bonds, self._guessBonds()))
    779 return bonds.astype(np.uint32)

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/molecule.py:940, in Molecule._guessBonds(self)
    934 """Tries to guess the bonds in the Molecule
    935 
    936 Can fail badly when non-bonded atoms are very close together. Use with extreme caution.
    937 """
    938 from moleculekit.bondguesser import guess_bonds
--> 940 return guess_bonds(self)

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/bondguesser.py:145, in guess_bonds(mol)
    142 # Setting the grid box distance to 1.2 times the largest VdW radius
    143 grid_cutoff = np.max(radii) * 1.2
--> 145 bonds = bond_grid_search(coords, grid_cutoff, is_hydrogen, radii)
    147 # bonds = guess_bonds(mol.coords, radii, is_hydrogen, grid_cutoff)
    148 return bonds

File ~/miniconda3/envs/htmd/lib/python3.10/site-packages/moleculekit/bondguesser.py:208, in bond_grid_search(coords, grid_cutoff, is_hydrogen, radii, max_boxes, cutoff_incr)
    206 result_list = []
    207 for boxidx in list(atoms_in_box.keys()):
--> 208     pairs = grid_bonds(
    209         coords,
    210         radii,
    211         is_hydrogen,
    212         pairdist,
    213         boxidx,
    214         atoms_in_box_full,
    215         gridlist,
    216     )
    217     if len(pairs) != 0:
    218         pairs = np.array(pairs, dtype=np.float32).reshape(-1, 2)

File moleculekit/bondguesser_utils/bondguesser_utils.pyx:160, in moleculekit.bondguesser_utils.grid_bonds()

MemoryError: std::bad_alloc

[stefdoerr]

v = Molecule('structure.pdb')
sel = v.atomselect('not segid P1 P3')

v.read('../output.vel', type='coor')
v.filter(sel)
v.write('x_out_test.vel', type='coor')

Try this, it might work. Doing an atomselection after loading the vel file will do rubbish because the velocities are loaded into the coords field of the object and as you see in the error it tries to guess bonds from atom coordinates (which are velocities here).

This is bad design from my side that we load the velocities into the coordinates field but handling the vel files were such an edge case that I didn't bother creating a new class attribute to store them.

[Saleh-OM4R]

Thank you @stefdoerr I've tried your code and it worked now,and thank you also for the explanation.