# -*- encoding: utf-8 -*- # möpß # # ---------------------------------------------------------------------------- # "THE BEER-WARE LICENSE" (Revision 42): # wrote this file. As long as you retain # this notice you can do whatever you want with this stuff. If we meet some day, # and you think this stuff is worth it, you can buy me a beer in return. # Daniel Kratzert # ---------------------------------------------------------------------------- # __doc__ = """ This is a full implementation of the SHELXL file syntax. Additionally it is able to edit SHELX properties with Python. The implementation is Python3-only and supports SHELXL after 2017 (You should not use old versions anyway). The parser is quiet about the most errors unless you enable debug or verbose during initialization. The parser will try to read the SHELX file even if it has syntax errors, but if for example, the SFAC and UNIT instruction is not consistent it will fail. """ import re import sys from contextlib import suppress from pathlib import Path from typing import Union, List, Optional from shelxfile.atoms.atom import Atom from shelxfile.atoms.atoms import Atoms from shelxfile.cif.cif_write import CifFile from shelxfile.misc.dsrmath import Array from shelxfile.misc.elements import weight_from_symbol # noinspection PyUnresolvedReferences from shelxfile.misc.misc import ParseOrderError, ParseNumError, ParseUnknownParam, \ multiline_test, dsr_regex, wrap_line, ParseSyntaxError from shelxfile.refine.refine import ShelxlRefine from shelxfile.shelx.cards import ACTA, FVAR, FVARs, REM, BOND, Restraints, DEFS, NCSY, ISOR, FLAT, \ BUMP, DFIX, DANG, SADI, SAME, RIGU, SIMU, DELU, CHIV, EADP, EXYZ, DAMP, HFIX, HKLF, SUMP, SYMM, LSCycles, \ SFACTable, UNIT, BASF, TWIN, WGHT, BLOC, SymmCards, CONN, CONF, BIND, DISP, GRID, HTAB, MERG, FRAG, FREE, FMAP, \ MOVE, PLAN, PRIG, RTAB, SHEL, SIZE, SPEC, STIR, TWST, WIGL, WPDB, XNPD, ZERR, CELL, LATT, MORE, MPLA, AFIX, PART, \ RESI, ABIN, ANIS, Residues, SWAT, Command, Restraint from shelxfile.shelx.sdm import SDM from shelxfile.version import VERSION __version__ = VERSION """ TODO: - Handle BEDE & LONE plus their results - Rotate ellipsoids with kabsch - killall.Q, killall.C - Q-peak printing is wrong: Q1 1 0.9828 1.1159 0.3148 11.00000 0.04 0.00 - Write out parts and afix in grow mode - Delete atoms (H) in AFIX -> delete entire afix group - Is backup file and reccovery from failed refinement working? ------------------------------------------------------------ - add remove_hydrogen_atoms(atom) method. - shx.remove_all_H([list of atoms], or all) - bond list - implement an add_afix(afixm, afixn, atoms, frag_fend=False, position=None, afix_options=None) default position is directly behind FVAR or FRAG/FEND if enabled - will read in lst file after refinement to fill shx.lst_file properties. - shx.move_in_part([list of atoms]) - shx.move_in_resi([list of atoms]) - shx.sort.residue(3) (low priority) - shx.lst_file.residuals -> density, r-values, goof, movements, bad reflections bad restraints, shelx *** errors - shx.unused_atom_name('C') -> get a carbon atom with unused number - shx.sort -> sort file (low priority) """ SHX_CARDS = ('TITL', 'CELL', 'ZERR', 'LATT', 'SYMM', 'SFAC', 'UNIT', 'LIST', 'L.S.', 'CGLS', 'BOND', 'FMAP', 'PLAN', 'TEMP', 'ACTA', 'CONF', 'SIMU', 'RIGU', 'WGHT', 'FVAR', 'DELU', 'SAME', 'DISP', 'LAUE', 'REM ', 'MORE', 'TIME', 'END ', 'HKLF', 'OMIT', 'SHEL', 'BASF', 'TWIN', 'EXTI', 'SWAT', 'HOPE', 'MERG', 'SPEC', 'RESI', 'MOVE', 'ANIS', 'AFIX', 'HFIX', 'FRAG', 'FEND', 'EXYZ', 'EADP', 'EQIV', 'CONN', 'BIND', 'FREE', 'DFIX', 'BUMP', 'SADI', 'CHIV', 'FLAT', 'DEFS', 'ISOR', 'NCSY', 'SUMP', 'BLOC', 'DAMP', 'STIR', 'MPLA', 'RTAB', 'HTAB', 'SIZE', 'WPDB', 'GRID', 'MOLE', 'XNPD', 'REST', 'CHAN', 'FLAP', 'RNUM', 'SOCC', 'PRIG', 'WIGL', 'RANG', 'TANG', 'ADDA', 'STAG', 'NEUT', 'ABIN', 'ANSC', 'ANSR', 'NOTR', 'TWST', 'PART', 'DANG', 'BEDE', 'LONE', 'REM', 'END') class Shelxfile(): """ Class for data from a SHELXL res file. Includes Atoms, cards and unit cell. """ _r1_regex = re.compile(r'^REM\s+R1\s+=', re.IGNORECASE) _wr2_regex = re.compile(r'^REM\s+wR2\s+=', re.IGNORECASE) _parameters_regex = re.compile(r'^REM\s+\d+\s+parameters\s+refined', re.IGNORECASE) _diff_peak_regex = re.compile(r'^REM\sHighest\sdifference', re.IGNORECASE) _goof_regex = re.compile(r'^REM\swR2\s=\s.*,\sGooF', re.IGNORECASE) _spgrp_regex = re.compile(r'^REM\s+\S+\s+in\s+\S+', re.IGNORECASE) def __init__(self, verbose: bool = False, debug: bool = False) -> None: if debug and verbose: raise ValueError("Either 'verbose' or 'debug' allowed, not both.") self.debug = False self.verbose = False if debug: self.debug = True if verbose: self.verbose = True # print(f'DEBUG: {self.debug}, VERBOSE: {self.verbose}') self.restraint_errors: List[str] = [] self.temp_in_kelvin: float = 0.0 self.shelx_max_line_length: int = 79 # maximum character lenth per line in SHELXL self.cell: Optional[CELL] = None self.ansc: List[float] = [] self.abin: Optional[ABIN] = None self.acta: Optional[ACTA] = None self.fmap: Optional[FMAP] = None self.xnpd: Optional[XNPD] = None self.wpdb: Optional[WPDB] = None self.wigl: Optional[WIGL] = None self.temp: Union[int, float] = 20 self.swat: Optional[SWAT] = None self.stir: Optional[STIR] = None self.spec: Optional[SPEC] = None self.twst: Optional[TWST] = None self.plan: Optional[PLAN] = None self.prig: Optional[PRIG] = None self.merg: Optional[MERG] = None self.more: Optional[MORE] = None self.move: Optional[MOVE] = None self.defs: Optional[DEFS] = None self.zerr: Optional[ZERR] = None self.wght: Optional[WGHT] = None self.frag: Optional[FRAG] = None self.twin: Optional[TWIN] = None self.basf: Optional[BASF] = None self.latt: Optional[LATT] = None self.anis: Optional[ANIS] = None self.damp: Optional[DAMP] = None self.unit: Optional[UNIT] = None self.size: Optional[SIZE] = None self.htab: Optional[HTAB] = None self.shel: Optional[SHEL] = None self.mpla: Optional[MPLA] = None self.hklf: Optional[HKLF] = None self.grid: Optional[GRID] = None self.conn: Optional[CONN] = None self.conf: Optional[CONF] = None self.afix: Optional[AFIX] = None self.rtab: List[RTAB] = [] self.omit: List[str] = [] self.free: List[FREE] = [] self.eqiv: List[str] = [] self.bonds: List[BOND] = [] self.disp: List[BOND] = [] self.bind: List[BIND] = [] self.bloc: List[BLOC] = [] self.part: PART = PART(self, ['PART', '0']) self.resi: RESI = RESI(self, ['RESI', '0']) self.residues = Residues(self) self.dsrlines: List[str] = [] self.dsrline_nums: List[int] = [] self.symmcards: SymmCards = SymmCards(self) self.hfixes: List[HFIX] = [] self.sump: List[SUMP] = [] self.wght_suggested: Optional[WGHT] = None self.Z: int = 1 self.titl: str = "" self.exti: float = 0.0 self.ansr: float = 0.001 self.rem: List[REM] = [] self.atoms: Atoms = Atoms(self) self.fvars: FVARs = FVARs(self) self.restraints: Restraints = Restraints() self.sfac_table: SFACTable = SFACTable(self) self.cycles: Optional[LSCycles] = None self.R1: Optional[float] = None self.wr2: Optional[float] = None self.goof: Optional[float] = None self.rgoof: Optional[float] = None self.space_group: Optional[str] = None self.data: Optional[int] = None self.parameters: Optional[int] = None self.dat_to_param: Optional[float] = None self.num_restraints: Optional[int] = None self.highest_peak: Optional[float] = None self.deepest_hole: Optional[float] = None self.formula_weight: Optional[float] = None self.end: bool = False self.maxsof: float = 1.0 self.delete_on_write: set = set() self.wavelength: float = 0.0 self.global_sadi: Optional[int] = None self.list: int = 0 self.theta_full: float = 0.0 self.error_line_num: int = -1 # Only used to tell the line number during an exception. self.resfile: Optional[Path] = None self._reslist: List[Union[str, Command, SFACTable, FVARs, Atom, SYMM]] = [] def write_shelx_file(self, filename: Union[str, Path, None] = None) -> None: if not self._reslist: print('*** No file was loaded for writing. ***') return None if not filename: filename = self.resfile if isinstance(filename, str): filename = Path(filename) with open(filename, 'w') as f: for num, line in enumerate(self._reslist): if num in self.delete_on_write: if self.debug: # print('Deleted line {}'.format(num + 1)) pass continue if line == '': # and self._reslist[num + 1] == '': continue # Prevent wrapping long lines with \n breaks by splitting first: line = "\n".join([wrap_line(x) for x in str(line).split("\n")]) f.write(str(line) + '\n') if self.verbose or self.debug: print(f'*** File successfully written to {filename.resolve()} ***') def read_file(self, resfile: Union[Path, str]) -> None: """ Read input from a file path. """ self.__init__(debug=self.debug, verbose=self.verbose) if isinstance(resfile, str): resfile = Path(resfile) self.resfile = resfile.resolve() if self.debug: print(f'Resfile is: {resfile}') try: self._reslist: List = resfile.read_text().splitlines(keepends=False) self._test_if_file_is_valid(resfile) except UnicodeDecodeError: if self.debug or self.verbose: print(f'*** Unable to read file {resfile} ***') return self._find_included_files() self.parse_cards() def read_string(self, resfile_string: str): """ Read input as string. This will not read files included with "+filename" syntax! """ self.__init__(debug=self.debug, verbose=self.verbose) self._reslist = resfile_string.splitlines(keepends=False) self.parse_cards() def parse_cards(self) -> None: try: self._parse_cards() except Exception as e: if self.debug or self.verbose: self.show_line_where_error_occured(e) if self.debug: raise else: return self.restraint_errors = self._assign_atoms_to_restraints() def _assign_atoms_to_restraints(self) -> List[str]: warnings = [] for restraint in self.restraints: bad_atoms = [] for restraint_atom in restraint.atoms: if restraint_atom in ('>', '<', '='): continue if (restraint.residue_class or sum(restraint.residue_number) > 0) and '_' not in restraint_atom: for num in restraint.residue_number: self.does_atom_exist(f'{restraint_atom}_{num}', bad_atoms, f'{restraint_atom}_{num}') elif '_' in restraint_atom: self.does_atom_exist(f'{restraint_atom}', bad_atoms, restraint_atom) else: self.does_atom_exist(f'{restraint_atom}_{0}', bad_atoms, restraint_atom) if bad_atoms: sorted_atoms = list(set(bad_atoms)) sorted_atoms.sort() warnings.append(f'*** Unknown atom{"s" if len(bad_atoms) > 0 else ""} in restraint: {restraint}, ' f'line {restraint.index + 1} ***') warnings.append(f'*** Atom list has no --> {", ".join(sorted_atoms)} ***') bad_atoms.clear() if restraint.residue_class and sum(restraint.residue_number) == 0: warnings.append(f"*** Restraint '{restraint}', line {restraint.index + 1}, " f"has a residue class, but no residues are defined. ***") if self.debug or self.verbose: print('\n'.join(warnings)) return warnings def does_atom_exist(self, atom_name: str, bad_atoms: List[str], restraint_atom: str): residue_number_is_wildcard = '_' in atom_name and atom_name.split('_')[-1] == '*' if atom_name.startswith('$'): return None elif residue_number_is_wildcard: for num in self.residues.residue_numbers.keys(): residue_atom = f"{atom_name.split('_')[0]}_{num}" if not self.atoms.get_atom_by_name(residue_atom): bad_atoms.append(residue_atom) else: if not self.atoms.get_atom_by_name(atom_name): bad_atoms.append(restraint_atom) def _test_if_file_is_valid(self, resfile: Path) -> None: if len(self._reslist) < 20 and (self.debug or self.verbose): print('*** Not a SHELXL file: {} ***'.format(resfile)) if self.debug: sys.exit() def show_line_where_error_occured(self, e): try: print(f'Error near:\n {self._reslist[self.error_line_num]}') except IndexError: pass print(e) print(f"*** Syntax error found in file {self.resfile}, line {self.error_line_num + 1} ***") def _find_included_files(self) -> None: # Tracks the file names of included files in order to find recursive inclusion: includefiles = [] for line_num, line in enumerate(self._reslist): if line.startswith('+'): try: file_included_in_includefile = self._read_included_file(includefiles, line) if file_included_in_includefile: for line_num_includefile, l in enumerate(file_included_in_includefile): reslist_position = line_num + 1 + line_num_includefile # '+filename' include files are not copied to res file, # so I have to delete these lines on write. # '++filename' copies them to the .res file where appropriate # I leave this out, because I am not SHELXL: # if l.startswith('+') and l[:2] != '++': # self.delete_on_write.update([lnum]) self._reslist.insert(reslist_position, l) continue except IndexError: if self.debug or self.verbose: print(f'*** CANNOT READ INCLUDE FILE {line} ***') # Not sure if this is a good idea: del reslist[n] def _read_included_file(self, includefiles: List[str], line: str): include_filename: Path = self.resfile.resolve().parent.joinpath(line[1:]) # Detect recursive file inclusion: if include_filename.name in includefiles: raise ValueError('*** Recoursive include files detected! ***') includefiles.append(include_filename.name) try: newfile = include_filename.read_text().splitlines(keepends=False) except IOError as e: if self.debug or self.verbose: print(e) print(f'*** CANNOT OPEN NESTED INPUT FILE {include_filename} ***') return [] return newfile def reload(self) -> None: """ Reloads the shelx file and parses it again. """ if self.debug or self.verbose: print(f'*** reloading file: {self.resfile} ***') self.read_file(self.resfile.resolve()) def _parse_cards(self) -> None: last_nonhydrogen_atom: Optional[Atom] = None lastcard = '' fvarnum = 1 for line_num, line in enumerate(self._reslist): self.error_line_num = line_num # For exception during parsing. list_of_lines = [line_num] # list of lines where a card appears, e.g. for atoms with two lines if line.startswith(' ') or line == '': continue wrapindex = 0 # This while loop makes wrapped lines look like they are not wrapped. The following lines are then # beginning with a space character and thus are ignored. The 'lines' list holds the line nnumbers where # 'line' is located ([line_num]) plus the wrapped lines. if multiline_test(self._reslist[line_num]): multiline = True else: multiline = False while multiline: # Glue together the two lines wrapped with "=": wrapindex += 1 line = line.rpartition('=')[0] + self._reslist[line_num + wrapindex] # self.delete_on_write.update([line_num + wrapindex]) list_of_lines.append(line_num + wrapindex) # list containing the lines of a multiline command # Do not activate this, otherwise, the unwrapping stops after two lines. if multiline_test(self._reslist[line_num + wrapindex]): multiline = True else: multiline = False self._reslist[line_num + wrapindex] = '' # The current line split: spline: list = line.split('!')[0].split() # Ignore comments with "!", see how this performes # The current line as string: line = line.upper().split('!')[0] # Ignore comments with "!", see how this performes word = line[:4] # get RESI: if line.startswith(('END', 'HKLF')) and self.resi: self.resi.num = 0 if self.debug or self.verbose: print('RESI in line {} was not closed'.format(line_num + 1)) # Do not continue here, otherwise HKLF is not parsed # continue if word == 'RESI': self.resi = self._assign_card(RESI(self, spline), line_num) if self.resi.residue_number > 0: self.residues.append(self.resi) continue # Now collect the PART: if line.startswith(('END', 'HKLF')) and self.part: self.part.n = 0 if self.debug or self.verbose: print('PART in line {} was not closed'.format(line_num + 1)) # Do not continue here, otherwise HKLF is not parsed # continue if word == 'PART': self.part = self._assign_card(PART(self, spline), line_num) continue # collect AFIX: if line.startswith(('END', 'HKLF')) and self.afix: self.afix.mn = 0 if self.debug or self.verbose: print('AFIX in line {} was not closed'.format(line_num + 1)) elif word == 'AFIX': self.afix = self._assign_card(AFIX(self, spline), line_num) elif self.is_atom(line): # A SHELXL atom: # F9 4 0.395366 0.177026 0.601546 21.00000 0.03231 ( 0.03248 = # 0.03649 -0.00522 -0.01212 0.00157 ) a = Atom(self) if last_nonhydrogen_atom: a.pivot = last_nonhydrogen_atom a.parse_line(spline, list_of_lines, part=self.part, afix=self.afix, resi=self.resi) if not a.is_hydrogen: last_nonhydrogen_atom = a a.pivot = None self._append_card(self.atoms, a, line_num) elif word == 'SADI': # SADI s[0.02] pairs of atoms # or SADI if len(spline) == 1: self.global_sadi = line_num self._append_card(self.restraints, SADI(self, spline), line_num) elif word == 'DFIX': # DFIX d s[0.02] atom pairs self._append_card(self.restraints, DFIX(self, spline), line_num) elif word == 'SIMU': # SIMU s[0.04] st[0.08] dmax[2.0] atomnames self._append_card(self.restraints, SIMU(self, spline), line_num) elif word == 'DELU': # DELU s1[0.01] s2[0.01] atomnames self._append_card(self.restraints, DELU(self, spline), line_num) elif word == 'RIGU': # RIGU s1[0.004] s2[0.004] atomnames self._append_card(self.restraints, RIGU(self, spline), line_num) elif word == 'BASF': # BASF scale factors self._assign_card(BASF(self, spline), line_num) elif word == 'HFIX': # HFIX mn U[#] d[#] atomnames self._append_card(self.hfixes, HFIX(self, spline), line_num) elif word == 'DANG': # DANG d s[0.04] atom pairs self._append_card(self.restraints, DANG(self, spline), line_num) elif word == 'EADP': self._append_card(self.restraints, EADP(self, spline), line_num) elif line.startswith('REM'): if dsr_regex.match(line): self.dsrlines.append(" ".join(spline)) self.dsrline_nums.extend(list_of_lines) self._append_card(self.rem, REM(self, spline), line_num) self._get_residuals(spline, line) elif word == 'CELL': # CELL λ a b c α β γ if lastcard != 'TITL' and (self.debug or self.verbose): print('*** TITL is missing. ***') self.cell = self._assign_card(CELL(self, spline), line_num) self.orthogonal_matrix = self.cell.o self.wavelength = self.cell.wavelen lastcard = 'CELL' elif word == "ZERR": # ZERR Z esd(a) esd(b) esd(c) esd(α) esd(β) esd(γ) if lastcard != 'CELL': if self.debug or self.verbose: print('*** Invalid SHELX file: CELL must occur before ZERR. ***') if self.debug: raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) if not self.cell: raise ParseOrderError('*** Cell parameters missing! ***', debug=shx.debug, verbose=shx.verbose) if len(spline) >= 8: self.zerr: ZERR = self._assign_card(ZERR(self, spline), line_num) self.Z = self.zerr.Z if self.Z < 1: self.Z = 1 if self.verbose or self.debug: print('*** Warning: Z value is zero. ***') lastcard = 'ZERR' elif word == "LATT": # LATT N[1] # 1=P, 2=I, 3=rhombohedral obverse on hexagonal axes, 4=F, 5=A, 6=B, 7=C. # negative is non-centrosymmetric self.latt = self._assign_card(LATT(self, spline), line_num) if lastcard != 'ZERR' and (self.verbose or self.debug): print('*** ZERR instruction is missing! ***') if self.latt.centric: self.symmcards.set_centric(True) elif word == "SYMM": # SYMM symmetry operation # Being more greedy, because many files do this wrong: # if not lastcard == 'ZERR': # raise ParseOrderError # if not self.zerr: # raise ParseOrderError s = SYMM(self, spline) if not self.latt and (self.debug or self.verbose): print("*** LATT instruction is missing! ***") if self.debug: raise ParseSyntaxError(debug=self.debug, verbose=self.verbose) # Have to do this after parsing, because P-1 has no SYMM! # if self.latt.centric: # self.symmcards.set_centric(True) self.symmcards.append(s.symmcard) if s not in self._reslist: self._reslist[line_num] = s else: self.delete_on_write.update([line_num]) self._reslist[line_num] = ' ' lastcard = 'SYMM' elif word == 'SFAC': # SFAC elements or # SFAC E a1 b1 a2 b2 a3 b3 a4 b4 c f' f" mu r wt # Being less strict to be able to parse files without cell errors: # if not (lastcard == 'LATT' or lastcard == 'ZERR'): # raise ParseOrderError # if not self.symmcards: # raise ParseOrderError if len(spline) <= 1: continue self.sfac_table.parse_element_line(spline) if self.sfac_table not in self._reslist: self._reslist[line_num] = self.sfac_table else: self.delete_on_write.update([line_num]) self._reslist[line_num] = ' ' lastcard = 'SFAC' elif word == 'UNIT': # UNIT n1 n2 ... # Number of atoms of each type in the unit-cell, in SFAC order. if lastcard != 'SFAC': raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) if self.sfac_table: try: self.unit = self._assign_card(UNIT(self, spline), line_num) except ValueError: if self.debug or self.verbose: print('*** Non-numeric value in SFAC instruction! ***') if self.debug: raise else: raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) if len(self.unit.values) != len(self.sfac_table.elements_list) and (self.debug or self.verbose): print('*** Number of UNIT and SFAC values differ! ***') if self.debug: raise ParseNumError(debug=self.shx.debug, verbose=self.shx.verbose) lastcard = 'UNIT' elif word in ['L.S.', 'CGLS']: # CGLS nls[0] nrf[0] nextra[0] # L.S. nls[0] nrf[0] nextra[0] self.cycles = self._assign_card(LSCycles(self, spline), line_num) elif word == "LIST": # LIST m[#] mult[1] (mult is for list 4 only) self.list = int(spline[1]) elif word == "FVAR": # FVAR osf[1] free variables for fvvalue in spline[1:]: fvarnum += 1 self._append_card(self.fvars, FVAR(fvarnum, float(fvvalue)), line_num) if self.fvars not in self._reslist: self._reslist[line_num] = self.fvars else: self.delete_on_write.update([line_num]) elif word == 'ANIS': # ANIS n or ANIS names # Must be before Atom(), to know which atom is anis. self.anis = self._assign_card(ANIS(self, spline), line_num) elif word == 'WGHT': # WGHT a[0.1] b[0] c[0] d[0] e[0] f[.33333] if self.end: self.wght_suggested = self._assign_card(WGHT(self, spline), line_num) continue self.wght = self._assign_card(WGHT(self, spline), line_num) elif word == 'ACTA': # ACTA 2θfull[#] -> optional parameter NOHKL self.acta = self._assign_card(ACTA(self, spline), line_num) elif word == 'DAMP': # DAMP damp[0.7] limse[15] self.damp = self._assign_card(DAMP(self, spline), line_num) elif word == 'ABIN': # ABIN n1 n2 self.abin = self._assign_card(ABIN(self, spline), line_num) elif word == 'ANSC': # ANSC six coefficients if len(spline) == 7: self.ansc = [float(x) for x in spline[:1]] elif word == 'ANSR': # ANSR anres[0.001] if len(spline) == 2: self.ansr = float(spline[1]) elif word == 'BIND': # BIND atom1 atom2 if len(spline) == 3: self._append_card(self.bind, BIND(self, spline), line_num) elif word == 'BLOC': # BLOC n1 n2 atomnames self._append_card(self.bloc, BLOC(self, spline), line_num) elif word == 'BOND': # BOND atomnames self._append_card(self.bonds, BOND(self, spline), line_num) elif word == 'BUMP': # BUMP s [0.02] self._assign_card(BUMP(self, spline), line_num) elif word == 'CHIV': # CHIV V[0] s[0.1] atomnames self._append_card(self.restraints, CHIV(self, spline), line_num) elif word == 'CONF': # CONF atomnames max_d[1.9] max_a[170] self.conf = self._assign_card(CONF(self, spline), line_num) elif word == 'CONN': # CONN bmax[12] r[#] atomnames or CONN bmax[12] # bonded are d < (r1 + r2 + 0.5) Å self.conn = self._assign_card(CONN(self, spline), line_num) elif word == 'DEFS': # DEFS sd[0.02] sf[0.1] su[0.01] ss[0.04] maxsof[1] self.defs = self._assign_card(DEFS(self, spline), line_num) elif word == 'DISP': # DISP E f' f"[#] mu[#] if lastcard != 'SFAC': raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) self._append_card(self.disp, DISP(self, spline), line_num) elif word == 'EQIV': # EQIV $n symmetry operation # TODO: implement EQUIV class if len(spline) > 1 and spline[1].startswith('$'): self.eqiv.append(spline[1:]) elif word == 'EXTI': # EXTI x[0] self.exti = float(spline[1]) elif word == 'EXYZ': # EXYZ atomnames self._append_card(self.restraints, EXYZ(self, spline), line_num) elif word == 'FRAG': # FRAG code[17] a[1] b[1] c[1] α[90] β[90] γ[90] if len(spline) == 8: self.frag = self._assign_card(FRAG(self, spline), line_num) elif word == 'FEND': # FEND (must follow FRAG) if not self.frag: raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) self.frag = None # Turns frag mode off. elif word == 'FLAT': # FLAT s[0.1] four or more atoms self._append_card(self.restraints, FLAT(self, spline), line_num) elif word == 'FREE': # FREE atom1 atom2 self._append_card(self.free, FREE(self, spline), line_num) elif word == 'GRID': # GRID sl[#] sa[#] sd[#] dl[#] da[#] dd[#] self.grid = self._assign_card(GRID(self, spline), line_num) elif word == 'HKLF': # HKLF N[0] S[1] r11...r33[1 0 0 0 1 0 0 0 1] sm[1] m[0] self.hklf = self._assign_card(HKLF(self, spline), line_num) elif line.startswith('END'): # END (after HKLF or ends an include file) self.end = True elif word == 'HTAB': # HTAB dh[2.0] or HTAB donor-atom acceptor-atom self.htab = self._assign_card(HTAB(self, spline), line_num) elif word == 'ISOR': # ISOR s[0.1] st[0.2] atomnames self._append_card(self.restraints, ISOR(self, spline), line_num) elif word == 'LAUE': # LAUE E # I completely do not understand the LAUE instruction description in the manual! continue elif word == 'MERG': # MERG n[2] self.merg = self._assign_card(MERG(self, spline), line_num) elif word == 'MORE': # MORE m[1] self.more = self._assign_card(MORE(self, spline), line_num) elif word == 'FMAP': # FMAP code[2] axis[#] nl[53] self.fmap = self._assign_card(FMAP(self, spline), line_num) elif word == 'MOVE': # MOVE dx[0] dy[0] dz[0] sign[1] self.move = self._assign_card(MOVE(self, spline), line_num) elif word == 'MPLA': # MPLA na atomnames self.mpla = self._assign_card(MPLA(self, spline), line_num) elif word == 'NCSY': # NCSY DN sd[0.1] su[0.05] atoms self._append_card(self.restraints, NCSY(self, spline), line_num) elif word == 'NEUT': # NEUT # TODO: Implement NEUT class if lastcard != 'SYMM': raise ParseOrderError(debug=shx.debug, verbose=shx.verbose) elif word == 'OMIT': # OMIT atomnames or OMIT s[-2] 2θ(lim)[180] or OMIT h k l # TODO: Implement OMIT class self.omit.append(spline[1:]) elif word == 'PLAN': # PLAN npeaks[20] d1[#] d2[#] self.plan = self._assign_card(PLAN(self, spline), line_num) elif word == 'PRIG': # PRIG p[#] self.prig = self._assign_card(PRIG(self, spline), line_num) elif word == 'RTAB': # RTAB codename atomnames --> codename: e.g. 'omeg' gets tabualted in the lst self._append_card(self.rtab, RTAB(self, spline), line_num) elif word == 'SAME': # SAME s1[0.02] s2[0.04] atomnames self._append_card(self.restraints, SAME(self, spline), line_num) elif word == 'SHEL': # SHEL lowres[infinite] highres[0] self.shel = self._assign_card(SHEL(self, spline), line_num) elif word == 'SIZE': # SIZE dx dy dz self.size = self._assign_card(SIZE(self, spline), line_num) elif word == 'SPEC': # SPEC del[0.2] if len(spline) > 1: self.spec = self._assign_card(SPEC(self, spline), line_num) elif word == 'STIR': # STIR sres step[0.01] -> stepwise improvement in the resolution sres self.stir = self._assign_card(STIR(self, spline), line_num) elif word == 'SUMP': # SUMP c sigma c1 m1 c2 m2 ... self._append_card(self.sump, SUMP(self, spline), line_num) elif word == 'SWAT': # SWAT g[0] U[2] self.swat = self._assign_card(SWAT(self, spline), line_num) elif word == 'TEMP': # TEMP T[20] -> in Celsius self.temp = float(spline[1].split('(')[0]) self.temp_in_kelvin = self.temp + 273.15 elif word == 'TWIN': # TWIN 3x3 matrix [-1 0 0 0 -1 0 0 0 -1] N[2] self.twin = self._assign_card(TWIN(self, spline), line_num) elif word == 'TWST': # TWST N[0] (N[1] after SHELXL-2018/3) if len(spline) > 1: self.twst = self._assign_card(TWST(self, spline), line_num) elif word == 'WIGL': # WIGL del[0.2] dU[0.2] self.wigl = self._assign_card(WIGL(self, spline), line_num) elif word == 'WPDB': # WPDB n[1] self.wpdb = self._assign_card(WPDB(self, spline), line_num) elif word == 'XNPD': # XNPD Umin[-0.001] self.xnpd = self._assign_card(XNPD(self, spline), line_num) elif word == 'BEDE': # Later... continue elif word == 'LONE': # Later... continue elif word == 'MOLE': # print('*** MOLE is deprecated! Do not use it! ***') pass elif word == 'HOPE': # print('*** HOPE is deprecated! Do not use it! ***') pass elif line.startswith('+'): pass elif word == 'TITL': self.titl = line[5:76] lastcard = 'TITL' else: if not line.strip(): continue if self.debug or self.verbose: print("Error in line: {} -> {}".format(line_num + 1, line)) if self.debug: raise ParseUnknownParam(debug=self.debug, verbose=self.verbose) def add_atom(self, name: str = None, coordinates: list = None, element='C', uvals: list = None, part: int = 0, sof: float = 11.0): """ Adds an atom to the ShelxFile.atoms list. If no element is given, carbon atoms are assumed. """ if uvals is None: uvals = [0.04] part = PART(self, f'PART {part}'.split()) afix = AFIX(self, 'AFIX 0'.split()) resi = RESI(self, 'RESI 0'.split()) a = Atom(self) sfac_num = self.elem2sfac(element) self.atoms._atomsdict.clear() a.set_atom_parameters(name=name, sfac_num=sfac_num, coords=coordinates, part=part, afix=afix, resi=resi, site_occupation=sof, uvals=uvals) self._append_card(self.atoms, a, 0) self.atoms._atomsdict.clear() def frac_to_cart(self, coordinates: list) -> Array: """ fractional to cartesian coordinates by applying the orthogonal matrix. """ return self.orthogonal_matrix * Array(coordinates) def __repr__(self): """ Represents the shelxl object. """ resl = [] for num, line in enumerate(self._reslist): if num in self.delete_on_write: continue try: if line == '' and self._reslist[num + 1] == '': continue except IndexError: pass # Prevent wrapping long lines with \n breaks by splitting first: line = "".join([wrap_line(x) for x in str(line).split("\n")]) resl.append(line) return "\n".join(resl) def grow(self, with_qpeaks: bool = False): """ Returns a list of atoms that represent the complete molecules of the structure. """ sdm = SDM(self) needsymm = sdm.calc_sdm() packed_atoms = sdm.packer(sdm, needsymm, with_qpeaks=with_qpeaks) return packed_atoms def refine(self, cycles: Union[int, None] = None, backup_before: bool = True) -> bool: if self.resfile: filen = self.resfile.stem if cycles is not None: self.cycles.number = cycles ref = ShelxlRefine(self, self.resfile) ref.remove_acta_card(self.acta) self.write_shelx_file(filen + '.ins') ref.run_shelxl(backup_before=backup_before) self.reload() ref.restore_acta_card() # self.write_shelx_file(filen + '.res') return True return False def refine_weight_convergence(self, stop_after: int = 10): """ Tries to refine weigting sheme from SHELXL until it converged (self.weight_difference() is zero) or stopt_after cycles are reached. """ for _ in range(stop_after): difference = self.wght.difference() print("Weighting difference = {} {}".format(*difference)) if self._weight_converged(difference): return True else: self.update_weight() self.refine(9) print("Maximum number of refinement cycles reached, but no WGHT convergence.") return False def _weight_converged(self, diff: List[float]) -> bool: return diff == [0.0, 0.0] def _append_card(self, obj, card, line_num: int) -> Command: """ Appends SHELX card to an object list, e.g. self.restraints and assigns the line_num in reslist with the card instance. """ obj.append(card) self._reslist[line_num] = card return card def _assign_card(self, card, line_num: int): self._reslist[line_num] = card return card @staticmethod def is_atom(atomline: str) -> bool: """ Returns True is line contains an atom. """ # no empty line, not in cards and not space at start: if atomline[:4].upper() not in SHX_CARDS: # exclude all non-atom cards spline: List[str] = atomline.split() # Too few parameter for an atom: if len(spline) < 5: return False # means sfac number is missing: if '.' in spline[1]: return False if Shelxfile._coordinates_are_unrealistic(spline): return False # Exclude lone pairs: if len(spline) > 5 and spline[5] == '!': return False return True else: return False @staticmethod def _coordinates_are_unrealistic(spline: List[str]) -> bool: return any(float(y) > 4.0 for y in spline[2:5]) def to_cif(self, filename: str = None, template: Optional[str] = None) -> None: """ Writes a CIF file from the ShelxFile object. """ if not filename: filename = self.resfile.stem + '.cif' CifFile(self, template).write_cif(Path(filename)) def elem2sfac(self, atom_type: str) -> int: """ returns an sfac-number for the element given in "atom_type" """ for num, element in enumerate(self.sfac_table, 1): if atom_type.capitalize() == element.capitalize(): return num # return sfac number # Element was not found in sfac table return 0 def sfac2elem(self, sfacnum: int) -> str: """ returns an element and needs an sfac-number :param sfacnum: string like '2' """ try: elem = self.sfac_table[int(sfacnum)] except IndexError: return '' return elem def add_line(self, linenum: int, line: str) -> None: """ Adds a new SHELX card to the reslist after linenum. e.g. shx.add_line(shx.unit.position, 'ANIS') """ self._reslist.insert(linenum + 1, line) def replace_line(self, obj, new_line: str) -> None: """ Replaces a single line in the res file with new_line. """ self._reslist[self.index_of(obj)] = new_line def index_of(self, obj: Union[Atom, Restraint, Command]) -> int: return self._reslist.index(obj) @property def sum_formula(self) -> str: """ The sum formula of the structure in regard to the UNIT instruction. """ formstring = '' formula_weight = 0.0 try: val = self.unit.values eli = self.sfac_table.elements_list except AttributeError: return '' if len(val) == len(eli): for el, num in zip(self.sfac_table.elements_list, self.unit.values): try: elcount = num / self.Z formula_weight += elcount * float(weight_from_symbol(el.capitalize())) formstring += f"{el}{elcount :,g} " except ZeroDivisionError: return '' self.formula_weight = round(formula_weight, 3) return formstring.strip() def update_weight(self) -> None: try: self.wght.a = self.wght_suggested.a self.wght.b = self.wght_suggested.b self.wght.c = self.wght_suggested.c self.wght.d = self.wght_suggested.d self.wght.e = self.wght_suggested.e self.wght.f = self.wght_suggested.f except AttributeError: return def insert_anis(self, atoms: str = '', residue: str = ''): """ Inserts ANIS into a results file for refinement of anisotropic displacement parameters. TODO: implement ANIS n @param atoms: Specify secific atoms or wildcards of atoms. @param residue: Specify a residue like ANIS_ABC with residue='ABC' or ANIS_* (residue='*') """ if atoms: self.add_line(self.unit.position, f'ANIS{"_" if residue else ""}{residue} {atoms}') else: self.add_line(self.unit.position, 'ANIS') @property def sum_formula_exact(self) -> str: """ The sum formula of the structure with all atom occupancies summed together as string. """ formstring = '' sumdict = self.sum_formula_exact_as_dict() for el in sumdict: formstring += f"{el.capitalize()}{round(sumdict[el], 2):,g} " return formstring.strip() def sum_formula_exact_as_dict(self) -> dict: """ The sum formula of the structure with all atom occupancies summed together as dictionary. """ sumdict = {} for el in self.sfac_table.elements_list: for atom in self.atoms: if atom.element.upper() == el.upper() and not atom.qpeak: if el in sumdict: sumdict[el] += atom.occupancy else: sumdict[el] = atom.occupancy if el not in sumdict: sumdict[el] = 0.0 return sumdict def insert_frag_fend_entry(self, dbatoms: list, cell: list): """ Inserts the FRAG ... FEND entry in the res file. :param dbatoms: list of atoms in the database entry :param cell: string with "FRAG 17 cell" from the database entry """ dblist = [] for line in dbatoms: dblist.append("{:4} {:<4} {:>8} {:>8} {:>8}".format(*line)) dblines = ' The following is from DSR:\n' dblines = dblines + 'FRAG 17 {} {} {} {} {} {}'.format(*cell) + '\n' dblines = dblines + '\n'.join(dblist) dblines = dblines + '\nFEND\n' # insert the db entry right after FVAR self.add_line(self.fvars.position, dblines) def _get_residuals(self, spline: List[str], line: str) -> None: if Shelxfile._r1_regex.match(line): self._get_r1(spline) if Shelxfile._wr2_regex.match(line): self._get_wr2(spline) if Shelxfile._parameters_regex.match(line): self._get_params_and_restraints(spline) if Shelxfile._diff_peak_regex.match(line): self._get_peak_hole(spline) if Shelxfile._goof_regex.match(line): self._get_goof(spline) if Shelxfile._spgrp_regex.match(line): self._get_space_group(spline) def _get_space_group(self, spline: List[str]) -> None: try: self.space_group = spline[3] except(IndexError, ValueError): pass def _get_goof(self, spline: List[str]) -> None: with suppress(IndexError, ValueError): self.goof = float(spline[8].split(',')[0]) self.rgoof = float(spline[12].split(',')[0]) def _get_peak_hole(self, spline: List[str]) -> None: # REM Highest difference peak 0.407, deepest hole -0.691, 1-sigma level 0.073 with suppress(IndexError, ValueError): self.highest_peak = float(spline[4].split(",")[0]) self.deepest_hole = float(spline[7].split(",")[0]) def _get_params_and_restraints(self, spline: List[str]) -> None: with suppress(IndexError): self.parameters = int(spline[1]) if self.data and self.parameters: self.dat_to_param = float(self.data) / float(self.parameters) with suppress(IndexError, ValueError): self.num_restraints = int(spline[-2]) def _get_wr2(self, spline: List[str]) -> None: with suppress(IndexError, ValueError): self.wr2 = float(spline[3].split(",")[0]) def _get_r1(self, spline: List[str]) -> None: with suppress(IndexError, ValueError): self.R1 = float(spline[3]) with suppress(IndexError, ValueError): self.data = int(spline[-2]) if __name__ == "__main__": print(Path('.').resolve()) # file = r'../shelxfile/tests/resources/p21c.res' file = r'tests/resources/p-31c.res' shx = Shelxfile() shx.read_file(file) print(shx.atoms) print(shx.sum_formula_exact) print(shx.sum_formula) print(shx.sum_formula_exact_as_dict()) print(shx.restraints) print(shx.atoms.nameslist) sys.exit() # noinspection PyUnreachableCode """ #To get all available SHELX commands: def get_shelx_commands(): url = "http://shelx.uni-goettingen.de/shelxl_html.php" response = urlopen('{}/version.txt'.format(url)) html = response.read().decode('UTF-8') #res = BeautifulSoup(html, "html5lib") tags = res.findAll("p", {"class": 'instr'}) for l in tags: if l: print(str(l).split(">")[1].split("<")[0]) """