Solvation_analysis subpackage multi-atom solute giving trouble

[aksharmaChem]

Hello,

I have a simulation trajectory from DFTB+ for 1 NaNO3 in 32 H2O molecules. The trajectory is in xyz format and I converted a frame of xyz to pdb using openbabel to serve as topology input.

I am following the instructions from the solvation analysis package documentation. However, I keep getting an assertion error

All residues must have the same length

I don't know how to get past this. My codes are pasted below from my Jupyter notebook.

`import MDAnalysis as mda
u=mda.Universe("./single-frame-1NaNO3-32-H2O.pdb", "./nano3-1-water-32-some-frames.xyz")
import solvation_analysis
nano3Atoms=u.atoms.select_atoms("id 1 or id 2 or id 3 or id 4 or id 5")
oAtomsWater=u.atoms.select_atoms("type O and resname HOH")
hAtomsWater=u.atoms.select_atoms("type H and resname HOH")
solventAtoms=u.atoms - nano3Atoms

from solvation_analysis.solute import Solute
soluteObject=Solute.from_atoms(nano3Atoms,{"Wat":solventAtoms},solute_name="salt")`

All the intermediate steps seem to present reasonable outputs with nano3Atoms showing all atoms that correspond to NaNO3. The last step when soluteObject is defined is when the assertion error is thrown. Any help would be greatly appreciated.

A pdb file and some frames from the simulation output are presented here.
Thank you,
frames.zip

[hmacdope]

Tagging original author @orionarcher :)

[orionarcher]

Thank you for the detailed bug report @aksharmaChem. Nothing stands out to me as to why this shouldn't work. Would you quickly check two things for me:

len(nano3Atoms) == 5
len(solventAtoms) == 96

If both of those evaluate to True, I'll look at this later this week and see if I can identify the issue.

[aksharmaChem]

Yes, both of those expressions evaluate to True.
Thank you so much for looking into this issue.

[aksharmaChem]

@orionarcher: Pinging you to send a reminder so that this request rises again. Thank you so much for looking into this.

[orionarcher]

Thanks for the ping, sorry this fell under my radar, will take a look later this week.

[orionarcher]

Hi @aksharmaChem

The problem is that the nano3Atoms group is actually two residues, one composed of NO3 and one composed of Na.

print([res.atoms.elements for res in nano3Atoms.residues])
[array(['N', 'O', 'O', 'O']), array(['Na'])]

The solute is intended to be only a single type of species, so it could be many Na residues or many NO3 residues, but not both. If you are interested in the solvation around both molecules, I'd suggest creating separate Solute's for each.

It's possible there is a way to forcibly merge the residues together but you'd have to check the MDAnalysis documentation.

[orbeckst]

@aksharmaChem does this answer your question? (If yes, please mark as answer, if no please say what's not working.)