MolarVerse / PQAnalysis

PQAnalysis is a API/CLI python package for the analysis of MD simulations
https://molarverse.github.io/PQAnalysis/
MIT License
4 stars 2 forks source link

shake_generator comment function #71

Closed Lucky-Armin closed 1 month ago

Lucky-Armin commented 1 month ago

Issue: When creating a topology file using the write_topology function I would like to have a comment appended for the equivalent bonds that I have averaged. In the code Snippet below i have a molecule with two equivalent hydrogen bond lengths which I labelled "inner_hydrogen" and "outer_hydrogen". It would be nice to have a custom comment appended in the shake.top output file based on averaged bonds, in this example "# inner H-C" and "# outer H-C". I have attached two files to show how it currently looks and how I want it to look (the .top file format is not supported by GitHub, so i appended a .txt).

Code Snippet:

from PQAnalysis.topology.shake_topology import ShakeTopologyGenerator
from PQAnalysis.io import TrajectoryReader
from PQAnalysis.topology import Selection
import numpy as np

file = ["molecule.xyz"]
reader = TrajectoryReader(file)
traj = reader.read()

selection = Selection("H")
h_indices = selection.select(traj.topology)
inner_hydrogen = np.array(h_indices)[[0, 5, 9, 10]]
outer_hydrogen = np.setdiff1d(h_indices, inner_hydrogen)

shake_generator = ShakeTopologyGenerator("H")
shake_generator.generate_topology(traj)
shake_generator.average_equivalents([inner_hydrogen, outer_hydrogen])
shake_generator.write_topology("shake.top", mode="o")

current_shake.top.txt desired_shake.top.txt

97gamjak commented 1 month ago

I have two possible suggestions - probably implementing both would be the way to go.

1) i add a add_line_comment(indices, comment) function

2) I add a function that takes in general a list of comments as input add_line_comments(list_of_comments)

The major difference is that in the second solution the user has to provide a list that is exactly as long as the number of shake lines generated - however, for the first solution the user has to now the exact indices of the shake atoms in the system (should not be a problem - they should be known in most of the cases)

Lucky-Armin commented 1 month ago

If I understand correctly the comments are then always based on the number of indices provided. Would it also be possible to provide a number of comments corresponding to the number of averaged bonds? So in this example 2: one for the inner and one for the outer H-C bonds.