Source code for prody.trajectory.psffile

# -*- coding: utf-8 -*-
"""This module defines a function for parsing protein structure files in
`PSF format`_.

.. _PSF format:
   http://www.ks.uiuc.edu/Training/Tutorials/namd/
   namd-tutorial-unix-html/node21.html"""

import os.path

from numpy import fromstring, zeros, ones, array, add

from prody import PY2K
from prody.atomic import ATOMIC_FIELDS, AtomGroup
from prody.utilities import openFile

if PY2K:
    range = xrange

__all__ = ['parsePSF', 'writePSF']

[docs]def parsePSF(filename, title=None, ag=None): """Return an :class:`.AtomGroup` instance storing data parsed from X-PLOR format PSF file *filename*. Atom and bond information is parsed from the file. If *title* is not given, *filename* will be set as the title of the :class:`.AtomGroup` instance. An :class:`.AtomGroup` instance may be provided as *ag* argument. When provided, *ag* must have the same number of atoms in the same order as the file. Data from PSF file will be added to the *ag*. This may overwrite present data if it overlaps with PSF file content. Note that this function does not evaluate angles, dihedrals, and impropers sections.""" if ag is not None: if not isinstance(ag, AtomGroup): raise TypeError('ag must be an AtomGroup instance') psf = openFile(filename, 'rb') line = psf.readline() i_line = 1 while line: line = line.strip() if line.endswith('!NATOM'): n_atoms = int(line.split('!')[0]) break line = psf.readline() i_line += 1 if title is None: title = os.path.splitext(os.path.split(filename)[1])[0] else: title = str(title) if ag is None: ag = AtomGroup(title) else: if n_atoms != ag.numAtoms(): raise ValueError('ag and PSF file must have same number of atoms') serials = zeros(n_atoms, ATOMIC_FIELDS['serial'].dtype) segnames = zeros(n_atoms, ATOMIC_FIELDS['segment'].dtype) resnums = zeros(n_atoms, ATOMIC_FIELDS['resnum'].dtype) resnames = zeros(n_atoms, ATOMIC_FIELDS['resname'].dtype) atomnames = zeros(n_atoms, ATOMIC_FIELDS['name'].dtype) atomtypes = zeros(n_atoms, ATOMIC_FIELDS['type'].dtype) charges = zeros(n_atoms, ATOMIC_FIELDS['charge'].dtype) masses = zeros(n_atoms, ATOMIC_FIELDS['mass'].dtype) lines = psf.readlines(71 * (n_atoms + 5)) if len(lines) < n_atoms: raise IOError('number of lines in PSF is less than the number of ' 'atoms') for i, line in enumerate(lines): if i == n_atoms: break i_line += 1 if len(line) <= 71: serials[i] = line[:8] segnames[i] = line[9:13].strip() resnums[i] = line[14:19] resnames[i] = line[19:23].strip() atomnames[i] = line[24:28].strip() atomtypes[i] = line[29:35].strip() charges[i] = line[35:44] masses[i] = line[50:60] else: items = line.split() serials[i] = items[0] segnames[i] = items[1] resnums[i] = items[2] resnames[i] = items[3] atomnames[i] = items[4] atomtypes[i] = items[5] charges[i] = items[6] masses[i] = items[7] i = n_atoms while 1: line = lines[i].split() if len(line) >= 2 and line[1] == '!NBOND:': n_bonds = int(line[0]) break i += 1 lines = ''.join(lines[i+1:]) + psf.read(n_bonds/4 * 71) array = fromstring(lines, count=n_bonds*2, dtype=int, sep=' ') if len(array) != n_bonds*2: raise IOError('number of bonds expected and parsed do not match') psf.close() ag.setSerials(serials) ag.setSegnames(segnames) ag.setResnums(resnums) ag.setResnames(resnames) ag.setNames(atomnames) ag.setTypes(atomtypes) ag.setCharges(charges) ag.setMasses(masses) array = add(array, -1, array) ag.setBonds(array.reshape((n_bonds, 2))) return ag
PSFLINE = ('%8d %-4s %-4d %-4s %-4s %-4s %10.6f %13.4f %11d\n')
[docs]def writePSF(filename, atoms): """Write atoms in X-PLOR format PSF file with name *filename* and return *filename*. This function will write available atom and bond information only.""" try: n_atoms, segments, rnums, rnames, names, types, charges, masses = ( atoms.numAtoms(), atoms._getSegnames(), atoms._getResnums(), atoms._getResnames(), atoms._getNames(), atoms._getTypes(), atoms._getCharges(), atoms._getMasses()) except AttributeError: raise TypeError('atoms must be an Atomic instance') if segments is None: segments = atoms._getChids() if segments is None: segments = ['UNK'] * n_atoms if rnums is None: rnums = ones(n_atoms, int) if rnames is None: rnames = ['UNK'] * n_atoms if names is None: raise ValueError('atom names are not set') if types is None: atomtypes = zeros(n_atoms, array(['a']).dtype.char + '1') long_fields = array([len(tp) for tp in types]).max() > 4 out = openFile(filename, 'w') write = out.write write('PSF{0}\n'.format( ' NAMD' if long_fields else '')) write('\n') write('{0:8d} !NTITLE\n'.format(1)) write(' REMARKS {0}\n'.format(str(atoms))) write('\n') write('{0:8d} !NATOM\n'.format(n_atoms)) for i in range(n_atoms): write(PSFLINE % (i + 1, segments[i], rnums[i], rnames[i], names[i], types[i], charges[i], masses[i], 0)) bonds = list(atoms._iterBonds()) if bonds: bonds = array(bonds, int) + 1 write('\n') write('{0:8d} !NBOND: bonds\n'.format(len(bonds))) for i, bond in enumerate(bonds): write('%8s%8s' % (bond[0], bond[1])) if i % 4 == 3: write('\n') if i % 4 != 3: write('\n') out.close() return filename
Read the Docs v: v1.5
Versions
latest
v1.5
Downloads
On Read the Docs
Project Home
Builds

Free document hosting provided by Read the Docs.