RecursionError when read/write psf file #1175
Description
Hi,
I got a RecursionError when I load the attached psf file to Parmed (v3.4.1) using Python 3.7.
My psf file has the topology of a coarse-grained protein, including only C-alpha atoms and the bonds between two adjacent C-alpha atoms. In total, there are 1024 atoms and 1023 bonds in this structure. This is a small system but when I load it using the following command, I got the RecursionError.
import sys
import parmed as pmd
a=load_file('1jz7_model_clean_ca.psf')~/software/miniconda3/lib/python3.7/site-packages/parmed/formats/registry.py in load_file(filename, *args, **kwargs)
195 _prune_argument(cls.parse, kwargs, 'hasbox')
196 _prune_argument(cls.parse, kwargs, 'skip_bonds')
--> 197 return cls.parse(filename, *args, **kwargs)
198 elif hasattr(cls, 'open_old'):
199 _prune_argument(cls.open_old, kwargs, 'structure')
~/software/miniconda3/lib/python3.7/site-packages/parmed/formats/psf.py in parse(filename)
57 The PSF file instance with all information loaded
58 """
---> 59 return CharmmPsfFile(filename)
60
61 #===================================================
~/software/miniconda3/lib/python3.7/site-packages/parmed/charmm/psf.py in newfunc(*args, **kwargs)
31 """ Catch the index error """
32 try:
---> 33 return func(*args, **kwargs)
34 except IndexError as e:
35 raise CharmmError('Array is too short: %s' % e)
~/software/miniconda3/lib/python3.7/site-packages/parmed/charmm/psf.py in __init__(self, psf_name)
312 # Assign all of the atoms to molecules recursively
313 tmp = psfsections['MOLNT'][1]
--> 314 set_molecules(self.atoms)
315 molecule_list = [a.marked for a in self.atoms]
316 if len(tmp) == len(self.atoms):
~/software/miniconda3/lib/python3.7/site-packages/parmed/charmm/psf.py in set_molecules(atoms)
695 if not atoms[i].marked:
696 tmp = [i]
--> 697 _set_owner(atoms, tmp, i, molecule_number)
698 # Make sure the atom indexes are sorted
699 tmp.sort()
~/software/miniconda3/lib/python3.7/site-packages/parmed/charmm/psf.py in _set_owner(atoms, owner_array, atm, mol_id)
710 if not partner.marked:
711 owner_array.append(partner.idx)
--> 712 _set_owner(atoms, owner_array, partner.idx, mol_id)
713 assert partner.marked == mol_id, 'Atom in multiple molecules!'
714
... last 1 frames repeated, from the frame below ...
~/software/miniconda3/lib/python3.7/site-packages/parmed/charmm/psf.py in _set_owner(atoms, owner_array, atm, mol_id)
710 if not partner.marked:
711 owner_array.append(partner.idx)
--> 712 _set_owner(atoms, owner_array, partner.idx, mol_id)
713 assert partner.marked == mol_id, 'Atom in multiple molecules!'
714
RecursionError: maximum recursion depth exceeded in comparison
I have checked the source code in parmed/charmm/psf.py. I found you have already set the recursion limits to prohibit this issue but I still can get the RecursionError. I tried to modify the set_molecules and _set_owner functions in the source code to monitor how many times _set_owner is called:
def set_molecules(atoms):
"""
Correctly sets the molecularity of the system based on connectivity.
"""
from sys import setrecursionlimit, getrecursionlimit
# Since we use a recursive function here, we make sure that the recursion
# limit is large enough to handle the maximum possible recursion depth we'll
# need (NATOM). We don't want to shrink it, though, since we use list
# comprehensions in list constructors in some places that have an implicit
# (shallow) recursion, therefore, reducing the recursion limit too much here
# could raise a recursion depth exceeded exception during a _Type/Atom/XList
# creation. Therefore, set the recursion limit to the greater of the current
# limit or the number of atoms
setrecursionlimit(max(len(atoms), getrecursionlimit()))
print(getrecursionlimit())
# Unmark all atoms so we can track which molecule each goes into
atoms.unmark()
# The molecule "ownership" list
owner = []
# The way I do this is via a recursive algorithm, in which
# the "set_owner" method is called for each bonded partner an atom
# has, which in turn calls set_owner for each of its partners and
# so on until everything has been assigned.
molecule_number = 1 # which molecule number we are on
for i in range(len(atoms)):
# If this atom has not yet been "owned", make it the next molecule
# However, we only increment which molecule number we're on if
# we actually assigned a new molecule (obviously)
if not atoms[i].marked:
tmp = [i]
_set_owner.count = 0
_set_owner(atoms, tmp, i, molecule_number)
print(_set_owner.count)
# Make sure the atom indexes are sorted
tmp.sort()
owner.append(tmp)
molecule_number += 1
return owner
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
def _set_owner(atoms, owner_array, atm, mol_id):
_set_owner.count += 1
""" Recursively sets ownership of given atom and all bonded partners """
atoms[atm].marked = mol_id
for partner in atoms[atm].bond_partners:
if not partner.marked:
owner_array.append(partner.idx)
_set_owner(atoms, owner_array, partner.idx, mol_id)
assert partner.marked == mol_id, 'Atom in multiple molecules!'
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++When I set the limit to 3000, I got the output
3000
1024
which means that the script set the limit to 3000 and _set_owner was called 1024 times. This makes sense because I have 1024 atoms in this molecule.
However, when I change the limit back to 1000 and then load the psf file, I found the code set the limit to 1024, which is expected, but after _set_owner was called 997 times, I got the RecursionError.
Does anyone know how could this happen? Is this a bug in parmed or some intrinsic feature of python?
Many thanks!
Yang