rsref module

Real space refinement; model-map fit.

usage:

• RSRef [options] [map-file] [input-pdb-file]

• rsref.py [options] [map-file] [input-pdb-file]

Module author: Michael S. Chapman <chapmanms@missouri.edu>

Authors:

Michael S. Chapman <chapmanms@missouri.edu>

Oregon Health & Science University & University of Missouri

Version:

1, Nov 26, 2024

Changed in version 12/10/10: Started

Changed in version 07/03/17: Rotamer command being added in

Changed in version 1.0.0: 09/26/20 Python 2.7 –> 3.6 cmd2: optparse –> argparse: changed calls and decorators

Todo

Extend Map.calculate (and evaluate) to accept aggregate of selections should the user only want to include refining atoms & neighbors in statistics. Expect up to 10-fold increase in speed for local refinement or high symmetry cases. The refining atoms would be in AtomicParameterization.changing, but not sure about symmetry neighbors. Care needed, because: (1) should not be applying for image refinement; (2) scaling would vary by selection.

Todo

Support optional input of the grad convergence criterion in fractional units relative to the starting objective function value.

Todo

copy impact from superpose - might require heavy duty consolidation.

class rsref.Arguments(imports=[], main=None, *args, **kwargs)

Bases: Arguments

Manager for command-line arguments from main and imported modules.

This is a template to be copied into modules, for the handling of command- line options.

methods export() and domestic() should be overridden by the module subclass with declarations of the command-line arguments needed by the module.

Parameters:

• imports ((list of) ArgumentParser(s)) – (list of) module(s) from which to obtain “parent” ArgumentParser objects for inclusion.

• main (bool|NoneType) – Add information and arguments appropriate for a main program (or not); if None, will determine by whether this class is defined in __main__.

domestic()

Defines options used only when main program and not when imported.

export()

Defines options used in both stand-alone and imported modes.

finalize()

Completion of the parser after program sub-class initialized.

class rsref.OldArguments(imports=[], main=True, *args, **kwargs)

Bases: ArgumentParser

Subclass, adding program-specific command-line options.

Parameters:

• imports ((list of) ArgumentParser method(s)) – method(s) adding arguments to ArgumentParser instance. (Alternative to ArgumentParser parents w/ better maintained groups.)

• main (bool) – Parser for main program (i.e. not listed within an ArgumentParser parents argument). Adds default program information.

static export(self)

Adding options used in both stand-alone and imported modes.

rsref.RSRef

alias of Tasks

class rsref.Tasks

Bases: Tasks, Tasks, Tasks

Program-specific tasks, supporting just-in-time execution.

Methods with names ending _prereq define tasks whose prerequisites will be met, if possible, prior to execution.

Superclass for atoms-related common tasks.

Warning

If this class defines pre-requisites on which other tasks depend, then this superclass must be instantiated before the other tasks are defined. If this class depends on pre-requisites defined elsewhere, then it must be instantiated afterwards. It may not be possible to satisfy both requirements!

analyze(itemize=True)

Analyze the gradients and shifts of prior refinement.

analyze_prereq(**kwargs)

Pass-through, trying to meet the prerequisites.

property atoms_being_evaluated

property atoms_refined

calculate_map(quiet=False, verbose=False)

Calculate atomic density for points in map.

calculate_map_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

property default_group

diagnostics()

Cumulative statistics/breakdown on calculation of derivatives.

evaluate(quiet=False, verbose=False, new_file=None)

Statistics comparing model & map (correlation, R, residual).

evaluate_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

image_refine(optimize_magnification=False, optimize_overall_b=False, optimize_envelope=False, optimize_filter=False, max_cycles=10, min_improvement=0.0, min_grad=0.0, finite_difference=False, defaults=None, verbosity=0, confirm=False)

Optimization of imaging (map) parameters for agreement w/ model.

Warning

Convergence is often poor when refining both magnification and and (overall B or envelope) - try refining B/envelope first.

See also

More discussion of convergence in density.ModelMap.image_refine().

Parameters:

• optimize_magnification (bool) – optimize magnification.

• optimize_overall_b (bool) – optimize overall_b.

• optimize_envelope (bool) – optimize envelope.

• optimize_filter (bool) – optimize the resolution of a low-pass filter to be applied to the model.

• max_cycles (int) – maximum number of max_cycles.

• min_improvement (float) – stopped when relative change in function < min_improvement, i.e. 1.e-2 would stop when change < 1% of current residual value. Rocognized by some optimizers (eg. bound-constrained l_bfgs_b), but ignored by others (eg. bfgs).

• min_grad (float) – refinement stopped when the largest derivative magnitude falls below this value (internal parameter units).

• defaults (object or NoneType) – optional object with attributes units_per_magnification, units_per_envelope & units_per_resolution, with which local defaults are over-written if all present.

• confirm (bool) – request user confirmation before accepting results.

image_refine_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

implied_parameterization(refining_xyz=True, refining_b=False, refining_occ=False)

Creates AtomicParameterization from Atoms.evaluate.

This is applicable (in embedded implementations) when there is only a single flag designating whether an atom is to be refined, and there is not the fine-grain control native to RSRef. Even though parameterization (rigid group, torsion, etc.) might be controlled by the embedding application, RSRef needs an AtomicParameterization instance to control which partial derivatives are to be calculated.

Note

requires self.atoms.evaluate to have been populated (see import_coord).

Parameters:

• refining_xyz (bool) – turn on the calculation of partial derivatives needed for any refinement of atomic positions.

• refining_b (bool) – turn on the calculation of partial derivatives needed for any refinement of atomic B-factors.

• refining_occ (bool) – turn on the calculation of partial derivatives needed for any refinement of occupancies.

Note

don’t like this implementation, because all partials for all atoms will be calculated.

import_coord(attribute, mylist)

Imports attribute into Atoms object from mylist array.

Parameters:

• attribute (str) – name of atoms attribute or x, y, or z for one of the xyz coordinates.

• mylist (list) – value for each atom

neighbors_deprecated(target_atoms=None, distance=None, verbose=True)

Identify self.atom neighbors of target (or self.atoms) + symm equivs.

Pass through routine for Atoms.neighbors or Symmetry.neighbors.

Parameters:

distance (float) – Angstrom within which considered close. Default calculated from Q.option.density_radius + Q.option.use_pts.

Todo

Consider adding to distance to allow for coordinate movement. May not be important as density_radius typically generous.

neighbors_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

parameterization_init()

Initialize parameterization - directives on what/how to refine.

parameterization_init_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

parameterize()

Pass-back after checking prerequisites for parameterize sub-comds.

parameterize_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

perturb(individual=False, displacement=1.0, direction=None, steps=1, repeats=1, verbose=True, exception=False, seed=None)

Shifts coordinates in steps, calculating comparison stats & gradients.

Purpose: to evaluate the convergence radius with the current options.

Perturbation may be rigid group or individual-atom, usually in random directions, and repeated to generate mean values & standard deviations. For a single rigid perturbation, a direction may be set. For rigid-group perturbation, displacement is the actual value, whereas for individual-atom, it is the mean of a normal distribution.

Parameters:

• individual (bool) – shift atoms independently, else as a rigid body.

• displacement (float) – overall Angstrom shift (rigid body) or rms shift (individual atom).

• direction (nd.array-like, shape(3,1) or Nonetype) – unit vector for rigid displacement or None if to be chosen randomly.

• steps (int) – number of steps into which the displacement is to be divided.

• repeats (int) – Number of random-direction shifts over-which statistics are to be averaged.

• verbose (bool)

• exception (bool) – terminate on error, else return None

• seed (int) – for random number generator

Returns:

(magnitude(s), correlation coefficient(s), rms, R-factor(s)).

Return type:

tuple of 4 {float or nd.array.size(steps)}

Attention

over-writes atomic map, but as the last loop is with zero-error, this should be benign.

perturb_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

profile(b=10.0, atom_type='C', max_dist=4.0, points=17)

Print calculated density profile of an isolated atom.

Parameters:

• b (float) – thermal factor

• atom_type (str)

• max_dist (float) – tabulate up to this distance from atom center

• points (int) – number of columns in table

randomize(**arguments)

randomize_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

reconcile_types()

Reconcile atom types of coordinate file with those of form factor data base.

reconcile_types_prereq()

Pass-through, trying to meet the prerequisites.

refine(continuation=False, restart=False, method='l_bfgs_b', max_cycles=0, min_improvement=0.0, accuracy=0.0, min_grad=-1.0, reference=None, verbosity=0)

Refine atomic model (currently free-atom w/o stereochemical restraints).

Warning

Joint refinement of occupancy and B-factors for individual atoms requires very high resolution data.

Note

Convergence can be improved by setting stringent convergence criteria, eg. min_improvement = 1e.-9.

Todo

arguments for link, unlink (brought from input).

See also

Until joint refinement implemented here, see alternative image optimization in image_refine() etc.

Parameters:

• continuation (bool) – If not restart, treat as a continuation of a prior batch if no change in parameterization and atoms since last refinement. False - force a new batch.

• restart (bool) – from start of prior refinement.

• method (str) – tested are ‘l_bfgs_b’, ‘bfgs’, but ‘Powell’, ‘CG’, ‘Newton_CG’ & ‘TNC’ are also available (see scipy.optimize.minimize documentation)

Note

l_bfgs is usually recommended as the memory usage of bfgs increases with the square of the number of parameters (atoms). l_bfgs_b also supports bound-constrained refinement.

Note

The convergence criteria below are used by the methods listed in parentheses (see scipy documentation), and are otherwise ignored.

Parameters:

• min_improvement (float) – refinement (l_bfgs_b, TNC) is stopped when relative change in function < min_improvement, i.e. precision of 1.e-2 stops when change < 1% of current residual value.

• accuracy (float) – criteria for convergence if estimated fractional error for all parameters is less than accuracy (Powell, Newton_CG).

• max_cycles (int) – maximum number of iterations (Powell, bfgs, Newton-CG, l_bfgs_b, TNC).

• min_grad (float) – refinement stopped when the largest derivative magnitude falls below this value (internal parameter units as reported by analyze) (CG, bfgs, l_bfgs_b, TNC).

• reference (Atoms or NoneType) – if torsion angle parameter changes are to be restrained (option.torsion_weight), an optional reference state of the moving structure (same atoms, order etc.) can be provided (eg. input structure); else any restraint will apply from the start of each refinement batch.

• verbosity (int) – -1 for silent, 0 for terse, 1 for verbose, > 1 for (debugging) information on each function or gradient evaluation.

Publications describing work using l_bfgs_b, or commercial products using it, should cite one of the Nocedal references [Byrd-1995] [Zhu-1997].

[Byrd-1995]

R. H. Byrd, P. Lu and J. Nocedal. A Limited Memory Algorithm for Bound Constrained Optimization, (1995), SIAM Journal on Scientific and Statistical Computing , 16, 5, pp. 1190-1208.

[Zhu-1997]

C. Zhu, R. H. Byrd and J. Nocedal. L-BFGS-B: Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization (1997), ACM Transactions on Mathematical Software, Vol 23, Num. 4, pp. 550 - 560.

Todo

Run extract at end with refined.x to save a copy of final parameters for later analysis if required.

refine_magnification()

Refine (only) magnification.

Deprecated since version (mostly): as functionality in self.image_refine, so have not defined prerequisites.

refine_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

property residual

Weighted local (rho_o - rho:sub:’c`)², per current options.

Return type:

float

property residual_db

Weighted partials of local residual re B-factors.

d((rho_o - rho:sub:’c`)²)/dB.

Return type:

1D ndarray length=#atoms

property residual_do

Weighted partials of residual re occupancies.

d((rho_o - rho:sub:’c`)²)/dO.

Return type:

1D ndarray length=#atoms

property residual_dx

Weighted partials of residual re x atomic coordinates.

d((rho_o - rho:sub:’c`)²)/dx.

Return type:

1D ndarray length=#atoms

property residual_dxyz

Weighted partials of residual re atomic positions.

d((rho_o - rho:sub:’c`)²)/dxyz.

Return type:

ndarray shape(3,#atoms)

property residual_dy

Weighted partials of residual re y atomic coordinates.

d((rho_o - rho:sub:’c`)²)/dy.

Return type:

1D ndarray length=#atoms

property residual_dz

Weighted partials of residual re z atomic coordinates: d((rho_o - rho:sub:’c`)²)/dz.

Return type:

1D ndarray length=#atoms

restrain_to_reference = True

rotamer(selection_arg, **kwargs)

Find best fitting standard rotamers for each side chain.

rotamer_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

run_test(name, arguments=None)

Program development testing.

property symmetry_current

True if current coordinates were same as last-expanded.

symmetry_expand(target_atoms=None, center=None, distance=None)

Identify self.atom neighbors of target (or self.atoms) + symm equivs.

Pass through for Symmetry.neighbors.

Parameters:

• center ([float, float, float]) – atoms w/in distance of this point instead of other atoms.

• distance (float) – Angstrom within which considered close.

symmetry_expand_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

symmetry_init(distance=None, force=False)

Initialize symmetry calculation & identify neighbors.

By default, initialization is conditional upon presence of space_group or local_symmetry, so lattice symmetry will only be applied if the space group is declared as at least P1. Unit_cell is not sufficient to invoke lattice symmetry, because it is required to interpret input maps and could be either a real lattice repeat or an arbitrary volume enclosing a single particle (in EM).

Parameters:

• distance (float) – neighbor cut-off; default = density radius + use_pts.

• force (bool) – initialize even if no space_group or local_symmetry.

symmetry_init_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

verbose = False

write_map(filename, observed=False, calculated=False, difference=False, scaled=False, title='')

Write map file

Parameters:

• filename (str) – output map

• observed (bool) – output experimental map

• scaled (bool) – output experimental map, scaled to atomic model

• calculated (bool) – output map calculated from atomic model

• difference (bool) – output difference map, observed minus model

• title (str) – put in the map header

write_map_prereq(*args, **kwargs)

Pass-through, trying to meet the prerequisites.

class rsref.TopCommands(tasks, initial_commands=None, py_locals={}, parser=None)

Bases: Commands, Commands, Commands, TopCmd, ProgCmd

Top-level commands for Program Prog.

Parameters:

• tasks (Contol) – class instance that defines methods to be called.

• initial_commands (str or list of str) – file(s) of commands to run before interactive loop at the top level, and before any command-line commands.

• py_locals (dict) – objects added to the local namespace in the python interpreter invoked by ‘py’ commands.

• parser (ArgumentParser (argparse') or :py:class:`OptionParser (optparse)) – to which -t / –test option will be added for cmd2 regression testing. Should only be used once / program.

do_evaluate(args)

Calculate statistics comparing current model to map.

do_image_refine(args)

Refine image (map) parameters to maximize agreement with atomic model.

do_map(args)

Write density map(s) on the grid of the input map.

do_neighbors(args)

Identify neighbors within distance of (selected) atoms.

do_parameterize(args)

Deisgnate how chosen atomic parameter(s) are to be refined.

do_perturb(args)

Perturb model & calculate statistics.

do_profile(args)

Density of single atom vs. distance; a help in setting cut-off radii

do_randomize(args)

Randomize coordinates (with normal error distributions).

do_refine(args)

Refine atomic model, per previously defined parameterization.

do_rotamer(args)

Compare to, and optionally reset amino acids to standard rotamers.

do_test(rest)

Run a program developer’s test.

evaluate_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

help_test()

image_refine_group = <argparse._MutuallyExclusiveGroup object>

image_refine_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

map_parser = ArgumentParser(prog='sphinx-build', usage=None, description='Write density map(s) on the grid of the input map.', formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

neighbors_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

parameterize_parser = ArgumentParser(prog='sphinx-build', usage=None, description='Deisgnate how chosen atomic parameter(s) are to be refined.', formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

perturb_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

profile_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

randomize_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

refine_convergence = <argparse._ArgumentGroup object>

refine_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

rotamer_parser = ArgumentParser(prog='sphinx-build', usage=None, description=None, formatter_class=<class 'argparse.HelpFormatter'>, conflict_handler='error', add_help=True)

rsref.main()

rsref.startup()

Program initialization, reading options etc..

Returns:

parser

Return type:

Arguments.ArgumentParser