Becksteinlab / GromacsWrapper

GromacsWrapper wraps system calls to GROMACS tools into thin Python classes (GROMACS 4.6.5 - 2024 supported).
https://gromacswrapper.readthedocs.org
GNU General Public License v3.0
169 stars 53 forks source link

"neutral" option in gromacs.genion #241

Closed marcodigennaro closed 1 year ago

marcodigennaro commented 1 year ago

Hello, I am trying to reproduce the tutorial from http://www.mdtutorials.com/gmx/lysozyme/index.html with GromacsWrapper. I am stuck on the ionic addtion for charge neutrality, specifically:

gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

I could not find the equivalent of "-neutral" in GromacsWrapper gromacs.genion There is a comment here: https://github.com/Becksteinlab/GromacsWrapper/blob/11f0daf0f054bff56e8c53888b13d7271c7bb471/gromacs/setup.py#L476 but it's one function in gromacs.setup, not in the gromacs.tools <class 'gromacs.tools.Genion'>

Can you advise on this? Thanks

marcodigennaro commented 1 year ago

You have to pass neutral=True to the command: gromacs.genion(s='ions.tpr', o='1AKI_solv_ions.gro', p=topol, pname="NA", nname="CL", neutral=True) and the expected behaviour is met. Going forward, how do I pass the index of molecules to be removed?

In the example, I get:

Command line:
  gmx genion -s ions.tpr -o 1AKI_solv_ions.gro -p topol.top -pname NA -nname CL -neutral

Reading file ions.tpr, VERSION 2023.1 (single precision)
Reading file ions.tpr, VERSION 2023.1 (single precision)
Will try to add 0 NA ions and 8 CL ions.
Select a continuous group of solvent molecules
Group     0 (         System) has 35968 elements
Group     1 (        Protein) has  1960 elements
Group     2 (      Protein-H) has  1001 elements
Group     3 (        C-alpha) has   129 elements
Group     4 (       Backbone) has   387 elements
Group     5 (      MainChain) has   517 elements
Group     6 (   MainChain+Cb) has   634 elements
Group     7 (    MainChain+H) has   646 elements
Group     8 (      SideChain) has  1314 elements
Group     9 (    SideChain-H) has   484 elements
Group    10 (    Prot-Masses) has  1960 elements
Group    11 (    non-Protein) has 34008 elements
Group    12 (          Water) has 34008 elements
Group    13 (            SOL) has 34008 elements
Group    14 (      non-Water) has  1960 elements

and I have to input 13 from keyboard. How do I pass it in the python command?

marcodigennaro commented 1 year ago

The answer is here:

gromacs.genion(s='ions.tpr', o='1AKI_solv_ions.gro', p=topol, pname="NA", nname="CL", neutral=True, input=("13"))

pass the desired value as a tuple to the input field of gromacs.genion

I guess I just needed to do some googling at this stage on my own to get the right answers

marcodigennaro commented 1 year ago

where can I find the energy module? i.e. the equivalent of gmx energy -f em.edr -o potential.xvg?