Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
psf = open(filename)
line = psf.readline()
i_line = 1
while line:
line = line.strip()
if line.strip().endswith('!NATOM'):
n_atoms = int(line.split('!')[0])
break
line = psf.readline()
i_line += 1
if name is None:
name = os.path.splitext(os.path.split(filename)[1])[0]
else:
name = str(name)
if ag is None:
ag = prody.AtomGroup(name)
else:
if n_atoms != ag.getNumOfAtoms():
raise ValueError('ag and PSF file must have same number of atoms')
serials = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['serial'].dtype)
segnames = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['segment'].dtype)
resnums = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['resnum'].dtype)
resnames = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['resname'].dtype)
atomnames = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['name'].dtype)
atomtypes = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['type'].dtype)
charges = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['charge'].dtype)
masses = np.zeros(n_atoms, ATOMIC_DATA_FIELDS['mass'].dtype)
lines = psf.readlines(71 * (n_atoms + 5))
index = True
split = False
def convert_chimera_molecule_to_prody(molecule):
"""
Function that transforms a chimera molecule into a prody atom group
Parameters
----------
molecule : chimera.Molecule
Returns
-------
prody_molecule : prody.AtomGroup()
chimera2prody : dict
dictionary: chimera2prody[chimera_atom.coordIndex] = i-thm element prody getCoords() array
"""
prody_molecule = prody.AtomGroup()
try:
coords, elements, names, resnums, chids, betas, masses = [], [], [], [], [], [], []
chimera2prody = {}
offset_chimera_residue = min(r.id.position for r in molecule.residues)
for i, atm in enumerate(molecule.atoms):
chimera2prody[atm.coordIndex] = i
coords.append(tuple(atm.coord())) # array documentation to improve
elements.append(atm.element.name)
names.append(atm.name)
resnums.append(atm.residue.id.position - offset_chimera_residue)
chids.append(atm.residue.id.chainId)
masses.append(atm.element.mass)
betas.append(atm.bfactor)
prody_molecule.setCoords(coords)
ag = opt.psf or opt.pdb
if ag:
if os.path.splitext(ag)[1].lower() == '.psf':
ag = prody.parsePSF(ag)
else:
ag = prody.parsePDB(ag)
dcd = prody.DCDFile(opt.coords)
if len(dcd) < 2:
opt.subparser("DCD file must contain multiple frames.")
if ag:
dcd.setAtomGroup(ag)
select = dcd.select(selstr)
LOGGER.info('{0:d} atoms are selected for calculations.'
.format(len(select)))
else:
select = prody.AtomGroup()
select.setCoords(dcd.getCoords())
pca = prody.PCA(dcd.getTitle())
if len(dcd) > 1000:
pca.buildCovariance(dcd)
pca.calcModes(dcd)
else:
pca.performSVD(dcd[:])
else:
pdb = prody.parsePDB(opt.coords)
if pdb.numCoordsets() < 2:
opt.subparser("PDB file must contain multiple models.")
if prefix == '_pca':
prefix = pdb.getTitle() + '_pca'
select = pdb.select(selstr)
LOGGER.info('{0:d} atoms are selected for calculations.'
.format(len(select)))
def convert_chimera_molecule_to_prody(molecule):
"""
Function that transforms a chimera molecule into a prody atom group
Parameters
----------
molecule : chimera.Molecule
Returns
-------
prody_molecule : prody.AtomGroup()
chimera2prody : dict
dictionary: chimera2prody[chimera_atom.coordIndex] = i-thm element prody getCoords() array
"""
prody_molecule = prody.AtomGroup()
try:
coords, elements, serials = [], [], []
names, resnums, resnames = [], [], []
chids, betas, masses = [], [], []
chimera2prody = {}
offset_chimera_residue = min(r.id.position for r in molecule.residues)
for i, atm in enumerate(molecule.atoms):
chimera2prody[atm.serialNumber] = i
coords.append(tuple(atm.coord())) # array documentation to improve
elements.append(atm.element.name)
serials.append(atm.serialNumber)
names.append(atm.name)
resnums.append(atm.residue.id.position - offset_chimera_residue)
resnames.append(atm.residue.type)
chids.append(atm.residue.id.chainId)
:class:`~prody.atomic.AtomGroup` instance *ag*. STRIDE output file must
be in the new format used from July 1995 and onwards. When *stride* file
is parsed, following attributes are added to *ag*:
* *stride_resnum*: STRIDE's sequential residue number, starting at the
first residue actually in the data set.
* *stride_phi*, *stride_psi*: peptide backbone torsion angles phi and psi
* *stride_area*: residue solvent accessible area
.. versionadded:: 0.8"""
if not os.path.isfile(stride):
raise IOError('{0:s} is not a valid file path'.format(stride))
if not isinstance(ag, prody.AtomGroup):
raise TypeError('ag argument must be an AtomGroup instance')
stride = open(stride)
n_atoms = ag.getNumOfAtoms()
NUMBER = np.zeros(n_atoms, int)
AREA = np.zeros(n_atoms, float)
PHI = np.zeros(n_atoms, float)
PSI = np.zeros(n_atoms, float)
ag.setSecondaryStrs(np.zeros(n_atoms))
for line in stride:
if not line.startswith('ASG '):
continue
res = ag[(line[9], int(line[10:15]), line[15].strip())]
if res is None:
elif subset.lower() not in _PDBSubsets:
raise ValueError('"{0:s}" is not a valid subset'.format(subset))
name_suffix = '_' + _PDBSubsets[subset]
if chain is not None:
if not isinstance(chain, str):
raise TypeError('chain must be a string')
elif len(chain) == 0:
raise ValueError('chain must not be an empty string')
name_suffix = '_' + chain + name_suffix
if 'ag' in kwargs:
ag = kwargs['ag']
if not isinstance(ag, prody.AtomGroup):
raise TypeError('ag must be an AtomGroup instance')
n_csets = ag.getNumOfCoordsets()
else:
ag = prody.AtomGroup(str(kwargs.get('name', 'unknown')) + name_suffix)
n_csets = 0
biomol = kwargs.get('biomol', False)
secondary = kwargs.get('secondary', False)
split = 0
hd = None
if model != 0:
lines = stream.readlines()
if header or biomol or secondary:
hd, split = _getHeaderDict(lines)
start = time.time()
_parsePDBLines(ag, lines, split, model, chain, subset, altloc)
if ag.getNumOfAtoms() > 0:
LOGGER.info('{0:d} atoms and {1:d} coordinate sets were '
'parsed in {2:.2f}s.'.format(ag.getNumOfAtoms(),
ag.getNumOfCoordsets() - n_csets, time.time()-start))
select = pdb.select(selstr)
dcd.setAtoms(select)
LOGGER.info('{0} atoms are selected for calculations.'
.format(len(select)))
else:
select = pdb.select(selstr)
if select.numAtoms() != dcd.numAtoms():
raise ValueError('number of selected atoms ({0}) does '
'not match number of atoms in the DCD '
'file ({1})'.format(select.numAtoms(),
dcd.numAtoms()))
if pdb.numCoordsets():
dcd.setCoords(select.getCoords())
else:
select = prody.AtomGroup()
select.setCoords(dcd.getCoords())
pca = prody.PCA(dcd.getTitle())
if len(dcd) > 1000:
pca.buildCovariance(dcd, aligned=kwargs.get('aligned'))
pca.calcModes(nmodes)
ensemble = dcd
else:
ensemble = dcd[:]
if not kwargs.get('aligned'):
ensemble.iterpose()
pca.performSVD(ensemble)
else:
pdb = prody.parsePDB(coords)
if pdb.numCoordsets() < 2:
raise ValueError('PDB file must contain multiple models')