1# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- 2# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 3# 4# MDAnalysis --- https://www.mdanalysis.org 5# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors 6# (see the file AUTHORS for the full list of names) 7# 8# Released under the GNU Public Licence, v2 or any higher version 9# 10# Please cite your use of MDAnalysis in published work: 11# 12# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, 13# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. 14# MDAnalysis: A Python package for the rapid analysis of molecular dynamics 15# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th 16# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. 17# 18# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. 19# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. 20# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 21# 22 23""" 24LAMMPSParser 25============ 26 27Parses data_ or dump_ files produced by LAMMPS_. 28 29.. _LAMMPS: http://lammps.sandia.gov/ 30.. _data: DATA file format: :http://lammps.sandia.gov/doc/2001/data_format.html 31.. _dump: http://lammps.sandia.gov/doc/dump.html 32 33 34.. _atom_style_kwarg: 35 36Atom styles 37----------- 38 39By default parsers and readers for Lammps data files expect either an 40*atomic* or *full* `atom_style`_. This can be customised by passing 41the `atom_style` keyword argument. This should be a space separated 42string indicating the position of the `id`, `type`, `resid`, `charge`, 43`x`, `y` and `z` fields. The `resid` and `charge` fields are optional 44and any other specified field will be ignored. 45 46For example to read a file with the following format, where there is no resid:: 47 48 Atoms # atomic 49 50 1 1 3.7151744275286681e+01 1.8684434743140471e+01 1.9285127961842125e+01 0 0 0 51 52 53The following code could be used:: 54 55 >>> import MDAnalysis as mda 56 >>> 57 >>> u = mda.Universe('myfile.data', atom_style='id type x y z') 58 59 60.. _`atom_style`: http://lammps.sandia.gov/doc/atom_style.html 61 62Classes 63------- 64 65.. autoclass:: DATAParser 66 :members: 67 :inherited-members: 68 69.. autoclass:: LammpsDumpParser 70 :members: 71 72Deprecated classes 73------------------ 74 75.. autoclass:: LAMMPSDataConverter 76 :members: 77 78""" 79from __future__ import absolute_import, print_function 80 81from six.moves import range 82 83import numpy as np 84import logging 85import string 86import functools 87import warnings 88 89from . import guessers 90from ..lib.util import openany, conv_float 91from ..lib.mdamath import triclinic_box 92from .base import TopologyReaderBase, squash_by 93from ..core.topology import Topology 94from ..core.topologyattrs import ( 95 Atomtypes, 96 Atomids, 97 Angles, 98 Bonds, 99 Charges, 100 Dihedrals, 101 Impropers, 102 Masses, 103 Resids, 104 Resnums, 105 Segids, 106) 107 108logger = logging.getLogger("MDAnalysis.topology.LAMMPS") 109 110 111# Sections will all start with one of these words 112# and run until the next section title 113SECTIONS = set([ 114 'Atoms', # Molecular topology sections 115 'Velocities', 116 'Masses', 117 'Ellipsoids', 118 'Lines', 119 'Triangles', 120 'Bodies', 121 'Bonds', # Forcefield sections 122 'Angles', 123 'Dihedrals', 124 'Impropers', 125 'Pair', 126 'Pair LJCoeffs', 127 'Bond Coeffs', 128 'Angle Coeffs', 129 'Dihedral Coeffs', 130 'Improper Coeffs', 131 'BondBond Coeffs', # Class 2 FF sections 132 'BondAngle Coeffs', 133 'MiddleBondTorsion Coeffs', 134 'EndBondTorsion Coeffs', 135 'AngleTorsion Coeffs', 136 'AngleAngleTorsion Coeffs', 137 'BondBond13 Coeffs', 138 'AngleAngle Coeffs', 139]) 140# We usually check by splitting around whitespace, so check 141# if any SECTION keywords will trip up on this 142# and add them 143for val in list(SECTIONS): 144 if len(val.split()) > 1: 145 SECTIONS.add(val.split()[0]) 146 147 148HEADERS = set([ 149 'atoms', 150 'bonds', 151 'angles', 152 'dihedrals', 153 'impropers', 154 'atom types', 155 'bond types', 156 'angle types', 157 'dihedral types', 158 'improper types', 159 'extra bond per atom', 160 'extra angle per atom', 161 'extra dihedral per atom', 162 'extra improper per atom', 163 'extra special per atom', 164 'ellipsoids', 165 'lines', 166 'triangles', 167 'bodies', 168 'xlo xhi', 169 'ylo yhi', 170 'zlo zhi', 171 'xy xz yz', 172]) 173 174 175class DATAParser(TopologyReaderBase): 176 """Parse a LAMMPS DATA file for topology and coordinates. 177 178 Note that LAMMPS DATA files can be used standalone. 179 180 Both topology and coordinate parsing functionality is kept in this 181 class as the topology and coordinate reader share many common 182 functions 183 184 By default the parser expects either *atomic* or *full* `atom_style` 185 however this can be by passing an `atom_style` keyword argument, 186 see :ref:`atom_style_kwarg`. 187 188 .. versionadded:: 0.9.0 189 """ 190 format = 'DATA' 191 192 def iterdata(self): 193 with openany(self.filename) as f: 194 for line in f: 195 line = line.partition('#')[0].strip() 196 if line: 197 yield line 198 199 def grab_datafile(self): 200 """Split a data file into dict of header and sections 201 202 Returns 203 ------- 204 header - dict of header section: value 205 sections - dict of section name: content 206 """ 207 f = list(self.iterdata()) 208 209 starts = [i for i, line in enumerate(f) 210 if line.split()[0] in SECTIONS] 211 starts += [None] 212 213 header = {} 214 for line in f[:starts[0]]: 215 for token in HEADERS: 216 if line.endswith(token): 217 header[token] = line.split(token)[0] 218 continue 219 220 sects = {f[l]:f[l+1:starts[i+1]] 221 for i, l in enumerate(starts[:-1])} 222 223 return header, sects 224 225 @staticmethod 226 def _interpret_atom_style(atom_style): 227 """Transform a string description of atom style into a dict 228 229 Required fields: id, type, x, y, z 230 Optional fields: resid, charge 231 232 eg: "id resid type charge x y z" 233 {'id': 0, 234 'resid': 1, 235 'type': 2, 236 'charge': 3, 237 'x': 4, 238 'y': 5, 239 'z': 6, 240 } 241 """ 242 style_dict = {} 243 244 atom_style = atom_style.split() 245 246 for attr in ['id', 'type', 'resid', 'charge', 'x', 'y', 'z']: 247 try: 248 location = atom_style.index(attr) 249 except ValueError: 250 pass 251 else: 252 style_dict[attr] = location 253 254 reqd_attrs = ['id', 'type', 'x', 'y', 'z'] 255 missing_attrs = [attr for attr in reqd_attrs if attr not in style_dict] 256 if missing_attrs: 257 raise ValueError("atom_style string missing required field(s): {}" 258 "".format(', '.join(missing_attrs))) 259 260 return style_dict 261 262 def parse(self, **kwargs): 263 """Parses a LAMMPS_ DATA file. 264 265 Returns 266 ------- 267 MDAnalysis Topology object. 268 """ 269 # Can pass atom_style to help parsing 270 try: 271 self.style_dict = self._interpret_atom_style(kwargs['atom_style']) 272 except KeyError: 273 self.style_dict = None 274 275 head, sects = self.grab_datafile() 276 277 try: 278 masses = self._parse_masses(sects['Masses']) 279 except KeyError: 280 masses = None 281 282 if 'Atoms' not in sects: 283 raise ValueError("Data file was missing Atoms section") 284 285 try: 286 top = self._parse_atoms(sects['Atoms'], masses) 287 except: 288 raise ValueError( 289 "Failed to parse atoms section. You can supply a description " 290 "of the atom_style as a keyword argument, " 291 "eg mda.Universe(..., atom_style='id resid x y z')" 292 ) 293 294 # create mapping of id to index (ie atom id 10 might be the 0th atom) 295 mapping = {atom_id: i for i, atom_id in enumerate(top.ids.values)} 296 297 for attr, L, nentries in [ 298 (Bonds, 'Bonds', 2), 299 (Angles, 'Angles', 3), 300 (Dihedrals, 'Dihedrals', 4), 301 (Impropers, 'Impropers', 4) 302 ]: 303 try: 304 type, sect = self._parse_bond_section(sects[L], nentries, mapping) 305 except KeyError: 306 pass 307 else: 308 top.add_TopologyAttr(attr(sect, type)) 309 310 return top 311 312 def read_DATA_timestep(self, n_atoms, TS_class, TS_kwargs, 313 atom_style=None): 314 """Read a DATA file and try and extract x, v, box. 315 316 - positions 317 - velocities (optional) 318 - box information 319 320 Fills this into the Timestep object and returns it 321 322 .. versionadded:: 0.9.0 323 .. versionchanged:: 0.18.0 324 Added atom_style kwarg 325 """ 326 if atom_style is None: 327 self.style_dict = None 328 else: 329 self.style_dict = self._interpret_atom_style(atom_style) 330 331 header, sects = self.grab_datafile() 332 333 unitcell = self._parse_box(header) 334 335 try: 336 positions, ordering = self._parse_pos(sects['Atoms']) 337 except KeyError as err: 338 raise IOError("Position information not found: {}".format(err)) 339 340 if 'Velocities' in sects: 341 velocities = self._parse_vel(sects['Velocities'], ordering) 342 else: 343 velocities = None 344 345 ts = TS_class.from_coordinates(positions, 346 velocities=velocities, 347 **TS_kwargs) 348 ts.dimensions = unitcell 349 350 return ts 351 352 def _parse_pos(self, datalines): 353 """Strip coordinate info into np array""" 354 pos = np.zeros((len(datalines), 3), dtype=np.float32) 355 # TODO: could maybe store this from topology parsing? 356 # Or try to reach into Universe? 357 # but ugly because assumes lots of things, and Reader should be standalone 358 ids = np.zeros(len(pos), dtype=np.int32) 359 360 if self.style_dict is None: 361 if len(datalines[0].split()) in (7, 10): 362 style_dict = {'id': 0, 'x': 4, 'y': 5, 'z': 6} 363 else: 364 style_dict = {'id': 0, 'x': 3, 'y': 4, 'z': 5} 365 else: 366 style_dict = self.style_dict 367 368 for i, line in enumerate(datalines): 369 line = line.split() 370 371 ids[i] = line[style_dict['id']] 372 373 pos[i, :] = [line[style_dict['x']], 374 line[style_dict['y']], 375 line[style_dict['z']]] 376 377 order = np.argsort(ids) 378 pos = pos[order] 379 380 # return order for velocities 381 return pos, order 382 383 def _parse_vel(self, datalines, order): 384 """Strip velocity info into np array 385 386 Parameters 387 ---------- 388 datalines : list 389 list of strings from file 390 order : np.array 391 array which rearranges the velocities into correct order 392 (from argsort on atom ids) 393 394 Returns 395 ------- 396 velocities : np.ndarray 397 """ 398 vel = np.zeros((len(datalines), 3), dtype=np.float32) 399 400 for i, line in enumerate(datalines): 401 line = line.split() 402 vel[i] = line[1:4] 403 404 vel = vel[order] 405 406 return vel 407 408 def _parse_bond_section(self, datalines, nentries, mapping): 409 """Read lines and strip information 410 411 Arguments 412 --------- 413 datalines : list 414 the raw lines from the data file 415 nentries : int 416 number of integers per line 417 mapping : dict 418 converts atom_ids to index within topology 419 """ 420 section = [] 421 type = [] 422 for line in datalines: 423 line = line.split() 424 # map to 0 based int 425 section.append(tuple([mapping[int(x)] for x in line[2:2 + nentries]])) 426 type.append(line[1]) 427 return tuple(type), tuple(section) 428 429 def _parse_atoms(self, datalines, massdict=None): 430 """Creates a Topology object 431 432 Adds the following attributes 433 - resid 434 - type 435 - masses (optional) 436 - charge (optional) 437 438 Lammps atoms can have lots of different formats, 439 and even custom formats 440 441 http://lammps.sandia.gov/doc/atom_style.html 442 443 Treated here are 444 - atoms with 7 fields (with charge) "full" 445 - atoms with 6 fields (no charge) "molecular" 446 447 Arguments 448 --------- 449 datalines - the relevent lines from the data file 450 massdict - dictionary relating type to mass 451 452 Returns 453 ------- 454 top - Topology object 455 """ 456 logger.info("Doing Atoms section") 457 458 n_atoms = len(datalines) 459 460 if self.style_dict is None: 461 sd = {'id': 0, 462 'resid': 1, 463 'type': 2} 464 # Fields per line 465 n = len(datalines[0].split()) 466 if n in (7, 10): 467 sd['charge'] = 3 468 else: 469 sd = self.style_dict 470 471 has_charge = 'charge' in sd 472 has_resid = 'resid' in sd 473 474 # atom ids aren't necessarily sequential 475 atom_ids = np.zeros(n_atoms, dtype=np.int32) 476 types = np.zeros(n_atoms, dtype=object) 477 if has_resid: 478 resids = np.zeros(n_atoms, dtype=np.int32) 479 else: 480 resids = np.ones(n_atoms, dtype=np.int32) 481 if has_charge: 482 charges = np.zeros(n_atoms, dtype=np.float32) 483 484 for i, line in enumerate(datalines): 485 line = line.split() 486 487 # these numpy array are already typed correctly, 488 # so just pass the raw strings 489 # and let numpy handle the conversion 490 atom_ids[i] = line[sd['id']] 491 if has_resid: 492 resids[i] = line[sd['resid']] 493 types[i] = line[sd['type']] 494 if has_charge: 495 charges[i] = line[sd['charge']] 496 497 # at this point, we've read the atoms section, 498 # but it's still (potentially) unordered 499 # TODO: Maybe we can optimise by checking if we need to sort 500 # ie `if np.any(np.diff(atom_ids) > 1)` but we want to search 501 # in a generatorish way, np.any() would check everything at once 502 order = np.argsort(atom_ids) 503 atom_ids = atom_ids[order] 504 types = types[order] 505 if has_resid: 506 resids = resids[order] 507 if has_charge: 508 charges = charges[order] 509 510 attrs = [] 511 attrs.append(Atomtypes(types)) 512 if has_charge: 513 attrs.append(Charges(charges)) 514 if massdict is not None: 515 masses = np.zeros(n_atoms, dtype=np.float64) 516 for i, at in enumerate(types): 517 masses[i] = massdict[at] 518 attrs.append(Masses(masses)) 519 else: 520 # Guess them 521 masses = guessers.guess_masses(types) 522 attrs.append(Masses(masses, guessed=True)) 523 524 residx, resids = squash_by(resids)[:2] 525 n_residues = len(resids) 526 527 attrs.append(Atomids(atom_ids)) 528 attrs.append(Resids(resids)) 529 attrs.append(Resnums(resids.copy())) 530 attrs.append(Segids(np.array(['SYSTEM'], dtype=object))) 531 532 top = Topology(n_atoms, n_residues, 1, 533 attrs=attrs, 534 atom_resindex=residx) 535 536 return top 537 538 def _parse_masses(self, datalines): 539 """Lammps defines mass on a per atom type basis. 540 541 This reads mass for each type and stores in dict 542 """ 543 logger.info("Doing Masses section") 544 545 masses = {} 546 for line in datalines: 547 line = line.split() 548 masses[line[0]] = float(line[1]) 549 550 return masses 551 552 def _parse_box(self, header): 553 x1, x2 = np.float32(header['xlo xhi'].split()) 554 x = x2 - x1 555 y1, y2 = np.float32(header['ylo yhi'].split()) 556 y = y2 - y1 557 z1, z2 = np.float32(header['zlo zhi'].split()) 558 z = z2 - z1 559 560 if 'xy xz yz' in header: 561 # Triclinic 562 unitcell = np.zeros((3, 3), dtype=np.float32) 563 564 xy, xz, yz = np.float32(header['xy xz yz'].split()) 565 566 unitcell[0][0] = x 567 unitcell[1][0] = xy 568 unitcell[1][1] = y 569 unitcell[2][0] = xz 570 unitcell[2][1] = yz 571 unitcell[2][2] = z 572 573 unitcell = triclinic_box(*unitcell) 574 else: 575 # Orthogonal 576 unitcell = np.zeros(6, dtype=np.float32) 577 unitcell[:3] = x, y, z 578 unitcell[3:] = 90., 90., 90. 579 580 return unitcell 581 582 583class LammpsDumpParser(TopologyReaderBase): 584 """Parses Lammps ascii dump files in 'atom' format 585 586 Only reads atom ids. Sets all masses to 1.0. 587 588 .. versionadded:: 0.19.0 589 """ 590 format = 'LAMMPSDUMP' 591 592 def parse(self, **kwargs): 593 with openany(self.filename) as fin: 594 fin.readline() # ITEM TIMESTEP 595 fin.readline() # 0 596 597 fin.readline() # ITEM NUMBER OF ATOMS 598 natoms = int(fin.readline()) 599 600 fin.readline() # ITEM BOX 601 fin.readline() # x 602 fin.readline() # y 603 fin.readline() # z 604 605 indices = np.zeros(natoms, dtype=int) 606 types = np.zeros(natoms, dtype=object) 607 608 fin.readline() # ITEM ATOMS 609 for i in range(natoms): 610 idx, atype, _, _, _ = fin.readline().split() 611 612 indices[i] = idx 613 types[i] = atype 614 615 order = np.argsort(indices) 616 indices = indices[order] 617 types = types[order] 618 619 attrs = [] 620 attrs.append(Atomids(indices)) 621 attrs.append(Atomtypes(types)) 622 attrs.append(Masses(np.ones(natoms, dtype=np.float64), guessed=True)) 623 warnings.warn('Guessed all Masses to 1.0') 624 attrs.append(Resids(np.array([1], dtype=int))) 625 attrs.append(Resnums(np.array([1], dtype=int))) 626 attrs.append(Segids(np.array(['SYSTEM'], dtype=object))) 627 628 return Topology(natoms, 1, 1, attrs=attrs) 629 630 631@functools.total_ordering 632class LAMMPSAtom(object): # pragma: no cover 633 __slots__ = ("index", "name", "type", "chainid", "charge", "mass", "_positions") 634 635 def __init__(self, index, name, type, chain_id, charge=0, mass=1): 636 self.index = index 637 self.name = repr(type) 638 self.type = type 639 self.chainid = chain_id 640 self.charge = charge 641 self.mass = mass 642 643 def __repr__(self): 644 return "<LAMMPSAtom " + repr(self.index + 1) + ": name " + repr(self.type) + " of chain " + repr( 645 self.chainid) + ">" 646 647 def __lt__(self, other): 648 return self.index < other.index 649 650 def __eq__(self, other): 651 return self.index == other.index 652 653 def __hash__(self): 654 return hash(self.index) 655 656 def __getattr__(self, attr): 657 if attr == 'pos': 658 return self._positions[self.index] 659 else: 660 super(LAMMPSAtom, self).__getattribute__(attr) 661 662 def __iter__(self): 663 pos = self.pos 664 return iter((self.index + 1, self.chainid, self.type, self.charge, self.mass, pos[0], pos[1], pos[2])) 665 666 667class LAMMPSDataConverter(object): # pragma: no cover 668 """Class to parse a LAMMPS_ DATA file and convert it to PSF/PDB. 669 670 The DATA file contains both topology and coordinate information. 671 672 The :class:`LAMMPSDataConverter` class can extract topology information and 673 coordinates from a LAMMPS_ data file. For instance, in order to 674 produce a PSF file of the topology and a PDB file of the coordinates 675 from a data file "lammps.data" you can use:: 676 677 from MDAnalysis.topology.LAMMPSParser import LAMPPSData 678 d = LAMMPSDataConverter("lammps.data") 679 d.writePSF("lammps.psf") 680 d.writePDB("lammps.pdb") 681 682 You can then read a trajectory (e.g. a LAMMPS DCD, see 683 :class:`MDAnalysis.coordinates.LAMMPS.DCDReader`) with :: 684 685 u = MDAnalysis.Unverse("lammps.psf", "lammps.dcd", format="LAMMPS") 686 687 .. deprecated:: 0.9.0 688 689 .. versionchanged:: 0.9.0 690 Renamed from ``LAMMPSData`` to ``LAMMPSDataConverter``. 691 """ 692 header_keywords = [ 693 "atoms", "bonds", "angles", "dihedrals", "impropers", 694 "atom types", "bond types", "angle types", 695 "dihedral types", "improper types", 696 "xlo xhi", "ylo yhi", "zlo zhi"] 697 698 connections = dict([ 699 ["Bonds", ("bonds", 3)], 700 ["Angles", ("angles", 3)], 701 ["Dihedrals", ("dihedrals", 4)], 702 ["Impropers", ("impropers", 2)]]) 703 704 coeff = dict([ 705 ["Masses", ("atom types", 1)], 706 ["Velocities", ("atoms", 3)], 707 ["Pair Coeffs", ("atom types", 4)], 708 ["Bond Coeffs", ("bond types", 2)], 709 ["Angle Coeffs", ("angle types", 4)], 710 ["Dihedral Coeffs", ("dihedral types", 3)], 711 ["Improper Coeffs", ("improper types", 2)]]) 712 713 def __init__(self, filename=None): 714 self.names = {} 715 self.headers = {} 716 self.sections = {} 717 if filename is None: 718 self.title = "LAMMPS data file" 719 else: 720 # Open and check validity 721 with openany(filename) as file: 722 file_iter = file.xreadlines() 723 self.title = file_iter.next() 724 # Parse headers 725 headers = self.headers 726 for l in file_iter: 727 line = l.strip() 728 if len(line) == 0: 729 continue 730 found = False 731 for keyword in self.header_keywords: 732 if line.find(keyword) >= 0: 733 found = True 734 values = line.split() 735 if keyword in ("xlo xhi", "ylo yhi", "zlo zhi"): 736 headers[keyword] = (float(values[0]), float(values[1])) 737 else: 738 headers[keyword] = int(values[0]) 739 if found is False: 740 break 741 742 # Parse sections 743 # XXX This is a crappy way to do it 744 with openany(filename) as file: 745 file_iter = file.xreadlines() 746 # Create coordinate array 747 positions = np.zeros((headers['atoms'], 3), np.float64) 748 sections = self.sections 749 for l in file_iter: 750 line = l.strip() 751 if len(line) == 0: 752 continue 753 if line in self.coeff: 754 h, numcoeff = self.coeff[line] 755 # skip line 756 file_iter.next() 757 data = [] 758 for i in range(headers[h]): 759 fields = file_iter.next().strip().split() 760 data.append(tuple([conv_float(el) for el in fields[1:]])) 761 sections[line] = data 762 elif line in self.connections: 763 h, numfields = self.connections[line] 764 # skip line 765 file_iter.next() 766 data = [] 767 for i in range(headers[h]): 768 fields = file_iter.next().strip().split() 769 data.append(tuple(np.int64(fields[1:]))) 770 sections[line] = data 771 elif line == "Atoms": 772 file_iter.next() 773 data = [] 774 for i in range(headers["atoms"]): 775 fields = file_iter.next().strip().split() 776 index = int(fields[0]) - 1 777 a = LAMMPSAtom(index=index, name=fields[2], type=int(fields[2]), chain_id=int(fields[1]), 778 charge=float(fields[3])) 779 a._positions = positions 780 data.append(a) 781 positions[index] = np.array([float(fields[4]), float(fields[5]), float(fields[6])]) 782 sections[line] = data 783 elif line == "Masses": 784 file_iter.next() 785 data = [] 786 for i in range(headers["atom type"]): 787 fields = file_iter.next().strip().split() 788 print("help") 789 self.positions = positions 790 791 def writePSF(self, filename, names=None): 792 """Export topology information to a simple PSF file.""" 793 # Naveen formatted -- works with MDAnalysis verison 52 794 #psf_atom_format = " %5d %-4s %-4d %-4s %-4s %-4s %10.6f %7.4f %1d\n" 795 # Liz formatted -- works with MDAnalysis verison 59 796 #psf_atom_format = "%8d %4.4s %-4.4s %-4.4s %-4.4s %-4.4s %16.8e %1s %-7.4f %7.7s %s\n" 797 # Oli formatted -- works with MDAnalysis verison 81 798 psf_atom_format = "%8d %4s %-4s %4s %-4s% 4s %-14.6f%-14.6f%8s\n" 799 with openany(filename, 'w') as file: 800 file.write("PSF\n\n") 801 file.write(string.rjust('0', 8) + ' !NTITLE\n\n') 802 file.write(string.rjust(str(len(self.sections["Atoms"])), 8) + ' !NATOM\n') 803 #print self.sections["Masses"] 804 for i, atom in enumerate(self.sections["Atoms"]): 805 if names is not None: 806 resname, atomname = names[i] 807 else: 808 resname, atomname = 'TEMP', 'XXXX' 809 for j, liz in enumerate(self.sections["Masses"]): 810 liz = liz[0] 811 #print j+1, atom.type, liz 812 if j + 1 == atom.type: 813 line = [ 814 i + 1, 'TEMP', 815 str(atom.chainid), resname, atomname, str(atom.type + 1), atom.charge, 816 float(liz), 0.] 817 else: 818 continue 819 #print line 820 file.write(psf_atom_format % tuple(line)) 821 822 file.write("\n") 823 num_bonds = len(self.sections["Bonds"]) 824 bond_list = self.sections["Bonds"] 825 file.write(string.rjust(str(num_bonds), 8) + ' !NBOND\n') 826 for index in range(0, num_bonds, 4): 827 try: 828 bonds = bond_list[index:index + 4] 829 except IndexError: 830 bonds = bond_list[index:-1] 831 bond_line = [string.rjust(str(bond[1]), 8) + string.rjust(str(bond[2]), 8) for bond in bonds] 832 file.write(''.join(bond_line) + '\n') 833 834 def writePDB(self, filename): 835 """Export coordinates to a simple PDB file.""" 836 atom_format = "%6s%.5s %4s %4s %.4s %8.3f%8.3f%8.3f%6.2f%6.2f %2s \n" 837 p = self.positions 838 with openany(filename, 'w') as file: 839 for i, atom in enumerate(self.sections["Atoms"]): 840 line = [ 841 "ATOM ", str(i + 1), 'XXXX', 'TEMP', str(atom.type + 1), p[i, 0], p[i, 1], p[i, 2], 0.0, 0.0, 842 str(atom.type)] 843 file.write(atom_format % tuple(line)) 844