QRef

QRef is a plugin for the macromolecular software suite Phenix enabling what in the literature commonly is referred to as "quantum refinement" (QR) in both real and reciprocal space, utilising the (for academic users) free software Orca as the quantum chemistry engine. This first version of QRef using Phenix was implemented under the paradigm "correctness and completeness first, performance and adherence to coding standards later".

Theory

Refinement of macromolecules in both real and and reciprocal space relies on previous knowledge (i.e. a Bayesian prior) for the structure. This is usually encoded as a (pseudo-energy) penalty term, $E_{restraints}(\mathbf{R})$, giving rise to a target function for the refinement with the general appearance

$$E_{total}\left(\mathbf{R}\right) = E_{exp}\left(\mathbf{R}\right) + wE_{restraints}\left(\mathbf{R}\right),$$

where $\mathbf{R}$ is the coordinate set for the current model and $w$ is a weight factor. $E_{restraints}(\mathbf{R})$ can in turn be broken down into its components:

$$E_{restraints}(\mathbf{R}) = E_{chem}(\mathbf{R}) + E_{SS}(\mathbf{R}) + E_{NCS}(\mathbf{R}) + \dots$$

In this QR implementation $E_{chem}(\mathbf{R})$ (which traditionally is a molecular mechanics force field) is replaced using a subtractive QM/MM scheme (using hydrogen link-atoms) according to

$$E_{chem}(\mathbf{R}) = \sum_{i} \left(w_{QM}E_{QM1, i}(\mathbf{R_{syst1, i}}) - E_{MM1, i}(\mathbf{R})\right) + E_{MM}(\mathbf{R}),$$

where index 1 in turn indicates small, but interesting, parts of the structure. The QM regions must not have QM atoms in common: a restraint between atoms that belong to several regions is included only once in $E_{MM}(\mathbf{R})$, but is subtracted once in each of their $E_{MM1, i}(\mathbf{R})$. Additionally another scaling factor, $w_{QM}$, is needed due to the fact that crystallographic MM force fields are of a statistical nature, whereas $E_{QM1, i}$ represents physical energies. $\mathbf{R_{syst1, i}}$ in turn is the coordinate set for the i:th region of QM atoms. Placing the hydrogen link-atoms at

$$\overline{r_{H_L}} = \overline{r_X} + g_{bond}\left(\overline{r_{C_L}} - \overline{r_X}\right)$$

implies that $\mathbf{R_{syst1, i}} = \mathbf{R_{syst1, i}}(\mathbf{R})$, thus the gradient for the chemical restraints is then obtained as

$$\nabla E_{restraints}(\mathbf{R}) = \sum_{i} \left( w_{QM} \nabla E_{QM1, i}(\mathbf{R_{syst1, i}(R)}) \cdot J(\mathbf{R_{syst1,i}}; \mathbf{R}) - \nabla E_{MM1, i}(\mathbf{R}) \right) + \nabla E_{MM}(\mathbf{R})$$

where $J(\mathbf{R_{syst1,i}}; \mathbf{R})$ is the Jacobian between $\mathbf{R_{syst1,i}}$ and $\mathbf{R}$. While this obviosuly is a matrix of size $3N_{syst1, i} \times 3N$, for a single junction in one QM system this becomes a 6x6 matrix with the general shape

$$\begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 \\\ 0 & 1 & 0 & 0 & 0 & 0 \\\ 0 & 0 & 1 & 0 & 0 & 0 \\\ (1 - g_{bond}) & 0 & 0 & g_{bond} & 0 & 0 \\\ 0 & (1 - g_{bond}) & 0 & 0 & g_{bond} & 0 \\\ 0 & 0 & (1 - g_{bond}) & 0 & 0 & g_{bond} \end{pmatrix}.$$

Installation

Installer

install.py installs QRef into the installation of Phenix given by $PHENIX:

source /path/to/phenix/phenix_env.sh
./install.py

The qref module is copied into that installation, and energies.py in cctbx and model.py in mmtbx are modified so that QRef is called during refinement. Each modified file is backed up as <file>.qref-orig.

--dry-run lists the changes without making them.

--uninstall reverses the installation.

--phenix DIR installs into DIR instead of $PHENIX, which is useful when more than one version of Phenix is present.

--source DIR installs from the QRef directory DIR instead of the one install.py itself is in.

--link makes the installation import qref from the QRef directory rather than copying it, so that several installations can share one copy of QRef. It also installs the scripts (see below) as wrappers that run them from the QRef directory, so a git pull updates both.

A file that is already patched is left alone, so running install.py again reports what it skipped and changes nothing.

--force patches an already patched file again, starting over from the .qref-orig backup.

QRef has been verified to work with Phenix versions 1.20.1-4487, 1.21-5207, 1.21.1-5286, 1.21.2-5419, 2.0-5936, 2.1-6084 and 2.2.1-6174.

Guard clause

QRef is only needed while coordinates are being refined. To restrict it to that, install.py also modifies

$PHENIX/.../phenix/refinement/xyz_reciprocal_space.py
$PHENIX/.../phenix/refinement/macro_cycle_real_space.py

so that a lock file exists while coordinates are refined, and QRef runs only when it does. Use --phenix-edits report to print these modifications and make them by hand instead, or --guard always to leave the two files untouched, at the cost of running QM calculations outside coordinate refinement as well.

Scripts

./install.py --scripts ~/bin --templates ~/qref

installs the scripts into a directory of your choice, which should be on $PATH. --templates copies junctfactor and qm_1.inp.

Per default, the supporting scripts run under whichever installation of Phenix is sourced when they are called. A pinned version of cctbx.python can be specified through --pin-shebang, so that they run without phenix_env.sh having been sourced.

Shifted models

phenix.real_space_refine does not shift the model back to the original coordinate frame when pdb_interpretation.sort_atoms = False. --shift-hint makes phenix.real_space_refine print the phenix.pdbtools translate= command at the end of refinement (can be found in the .log file), which can then be used to correct the output model.

Orca

Orca can be found at orcaforum.kofo.mpg.de - follow their guide for installation. QRef has been verified to work with Orca versions 5.0.4, 6.0.0, 6.1.0 and 6.1.1.

Usage

The general procedure to set up a quantum refinement job consists of

1. Select QM region(s). Best practises for selecting proper QM region(s) can be found at for example:

   • doi.org/10.1002/anie.200802019
   • doi.org/10.1016/bs.mie.2016.05.014
   • doi.org/10.1021/cr5004419

2. Build a model.

   • The model in the QM regions needs to make chemical sense. This for example means that the QM regions should be protonated as well as being complete.
   • The model outside of the QM region (as well as the protonation of the carbon link atom) can be incomplete.
   • phenix.ready_set add_h_to_water=True can be useful for this purpose.

3. Prepare restraint files for unknown residues and ligands. The script qref_prep.py will tell you if there are any missing restraint files.

   • This can be achieved using phenix.ready_set and phenix.elbow.

4. Prepare syst1 files; these files define the QM regions.

   • If there is only one QM region the default is to look for a file named syst1 by the software. For multiple QM regions the recommended, and default, naming scheme is syst11, syst12, etc.
   • QM regions must not share atoms other than link atoms; qref_prep.py exits if any other atom occurs in more than one syst1 file. For a QM region with alternative conformations, duplicate the QM region into separate altlocs in the model and define one QM region per altloc, e.g. syst11 for altloc A and syst12 for altloc B.
   • Which atoms to include in the QM regions is defined using the serial number from the PDB file describing the entire model.
     • While setting sort_atoms = False in the input to Phenix should ensure that the ordering in the input model is preserved, we have encountered instances where this is not adhered to. Thus it is recommended to use iotbx.pdb.sort_atoms (supplied with Phenix) which will give you a new PDB file with the suffix _sorted.pdb where the atoms are sorted in the same order as that which Phenix uses internally. It is recommended to use the _sorted.pdb file as the input model for refinement, as well as the reference when defining the syst1 files.
   • The syst1 files allows for multiple atoms or intervals of atoms to be specified on a single line, where , or blank works as delimiters; - is used to indicate an interval.
   • # and ! can be used to include comments in the syst1 files.
   • The second occurence of an atom in the syst1 files will indicate that this is a link atom, i.e. it will be replaced by a hydrogen at the appropriate position in the QM calculation.
   • Examples are included in the examples folder.

5. Run qref_prep.py <model>_sorted.pdb to generate qref.dat (a file containing settings for the QR interface), as well as PDB files describing the QM regions.

   • The junctfactor file needs to be present in the same directory as where qref_prep.py is run. The junctfactor file contain ideal QM distances for the $C_L - H_L$ bonds for the link-atoms. If another junctfactor file is to be used this can be specified with the -j or --junctfactor option.

   • The theory used for the ideal $C_L - H_L$ QM distance can be changed with the -l or --ltype option. Default is type 12. The options are:

     • 6: B3LYP/6-31G*
     • 7: BP(RI)/6-31G*
     • 8: BP(RI)/SVP
     • 9: BP(RI)/def2-SV(P)
     • 10: PBE(RI)/def2-SVP
     • 11: B3LYP(RI)/def2-SV(P)
     • 12: TPSS(RI)/def2-SV(P)
     • 13: B97-D(RI)/def2-SV(P)

     There needs to be a parametrisation in the junctfactor file for the bond one intends to cleave; it is recommended that the user inspects the junctfactor file to verify that there is support to cleave the intended bond type. In the case parametrisation is lacking another selection for the QM system (and in particular where the link between QM and MM occurs) needs to be made or appropriate parametrisation added to the junctfactor file.

   • Ideally only the input model is needed as an argument for qref_prep.py. If there was a need to prepare restraint files for novel residues or ligands in point 3 above, qref_prep.py needs to be made aware of these. This can be achieved with the -c or --cif option.

   • The output from qref_prep.py should be $\left(1+2 n_{syst1}\right)$ files as well as recommended selection strings for both real and reciprocal space. Additionally, a warning is given if the syst1 file covers more than one conformation. In the case that all atoms in the syst1 definition belong to the same conformation, this will be indicated by an altloc specifier.

     • qref.dat, which contains the settings for the QR interface. This file can be changed manually and it is a good idea to inspect that the value for orca_binary is the correct path for the actual Orca binary file (qref_prep.py tries to locate this file automatically but may sometimes fail).
     • mm_i_c.pdb, which is the model used to calculate $E_{MM1, i}$.
     • qm_i_h.pdb, which is the model used to calculate $E_{QM1, i}$.

     The output PDB files can, and probably should, be used to inspect that the QM selection is proper.

     • Two selection strings are printed on the screen, one for reciprocal space and one for real space. They are intended to be used in regards to which selection of the model to refine when crafting the input to either phenix.refine or phenix.real_space_refine, see point 7 below.

   • Harmonic (bond) distance restraints can be added through the -rd or --restraint_distance option, using the syntax i atom1_serial atom2_serial desired_distance_in_Å force_constant. The force constant corresponds to $1/\sigma^2$ of a restraint in Phenix, with $\sigma$ in Å. Experience has shown that it needs to be $\geq$ 2500 ($\sigma \leq$ 0.02 Å) to achieve adherence to the restraint.

   • Harmonic (bond) angle restraints can be added through the -ra or --restraint_angle option, using the syntax i atom1_serial atom2_serial atom3_serial desired_angle_in_degrees force_constant, where atom2_serial defines the angle tip. The force constant corresponds to $1/\sigma^2$ of a restraint in Phenix, with $\sigma$ in degrees. A force constant of 1 ($\sigma$ = 1°) is a reasonable starting point.

   • Symmetry interactions are handled through the -t or --transform option, using the syntax i "<atoms>" R11 R12 R13 R21 R22 R23 R31 R32 R33 t1 t2 t3 for each of the desired transforms, where <atoms> follow the same syntax as for syst1 files (do note the usage of quotation marks, i.e. <atoms> should be given as a string). R is the rotation matrix (in Cartesian coordinates) in row-wise order and t is the translation vector (in Cartesian coordinates). To obtain $R$ and $t$, find for example the fractional rotation matrix $R_{frac}$ and the fractional translation vector $t_{frac}$ for the symmetry operator of interest (which can be found through using for example Coot), together with the desired fractional unit cell translation vector, $t_{u}$. $R$ is then calculated as $R = S^{-1}R_{frac}S$ and $t = S^{-1}(t_{frac} + t_{u})$, where $S$ is the Cartesian to fractional conversion matrix for the space group the crystal belongs to (which can be found in the SCALEn records in a PDB file).

   • All available options for qref_prep.py can be seen through -h or --help.

6. The next step is to prepare the input files for Orca. Examples can be found in the examples folder.

   • The input files should be named qm_i.inp, where i refers to the i:th QM region; counting starts at 1.
   • The level of theory should match the junction type; using the default (type 12) the corresponding input to Orca then becomes ! TPSS D4 DEF2-SV(P).
   • Energy and gradient needs to be written to disk. This is achieved through the keyword ! ENGRAD.
   • To read coordinates from a PDB file (qm_i_h.pdb) use *pdbfile <charge> <multiplicity> qm_i_h.pdb, where again i refers to the i:th QM region.
   • Custom settings can also be supplied (see the examples folder).

7. At this point it is possible to start a quantum refinement. It is however recommended to first create an empty file named qm.lock (this disables QRef - additionally, if qref.dat is not present in the folder QRef will not run, i.e. this step can be run before starting to set up the QR job) through for example the touch command, then start a refinement so that an .eff file is obtained (this file contains all the Phenix refinement settings for the current experimental data and model). A good idea is to rename this file to phenix.params or similar, then edit this file and make sure the following options are set:

   • For reciprocal space refinement (phenix.refine):

     • refinement.pdb_interpretation.restraints_library.cdl = False

     • refinement.pdb_interpretation.restraints_library.mcl = False

     • refinement.pdb_interpretation.restraints_library.cis_pro_eh99 = False

     • refinement.pdb_interpretation.secondary_structure.enabled = False

     • refinement.pdb_interpretation.sort_atoms = False

     • refinement.pdb_interpretation.flip_symmetric_amino_acids = False

     • refinement.refine.strategy = *individual_sites individual_sites_real_space rigid_body *individual_adp group_adp tls occupancies group_anomalous

     • refinement.refine.sites.individual = <reciprocal selection string>

     • refinement.hydrogens.refine = *individual riding Auto

     • refinement.hydrogens.real_space_optimize_x_h_orientation = False

     • refinement.main.nqh_flips = False

     • If the QM regions have been duplicated into separate altlocs (see point 4) and occupancies are refined, the occupancies of the whole QM region should be constrained as one group. Using the reciprocal selection string printed by qref_prep.py:

       refinement.refine.occupancies {
         constrained_group {
           selection = altloc A and (<reciprocal selection string>)
           selection = altloc B and (<reciprocal selection string>)
         }
       }
       

   • For real space refinement (phenix.real_space_refine):
     • refinement.run = *minimization_global rigid_body local_grid_search morphing simulated_annealing adp occupancy nqh_flips
     • pdb_interpretation.restraints_library.cdl = False
     • pdb_interpretation.restraints_library.mcl = False
     • pdb_interpretation.restraints_library.cis_pro_eh99 = False
     • pdb_interpretation.flip_symmetric_amino_acids = False
     • pdb_interpretation.sort_atoms = False
     • pdb_interpretation.secondary_structure = False
     • pdb_interpretation.reference_coordinate_restraints.enabled = True
     • pdb_interpretation.reference_coordinate_restraints.selection = <real space selection string>
     • pdb_interpretation.reference_coordinate_restraints.sigma = 0.01
     • pdb_interpretation.ramachandran_plot_restraints.enabled = False
   • Other options can be set as needed.

8. To run the quantum refinement job, make sure that the qm.lock file has been deleted, then execute either phenix.refine phenix.params or phenix.real_space_refine phenix.params. If there is a need to restart the job with different settings for Phenix, make sure to delete the file settings.pickle.

Notes

For COSMOS@LUNARC users

When using ORCA/6.0.0 either manually run (before you submit your job) or add to the beginning of your submit script the exports below:

export OMPI_MCA_btl='^uct,ofi'
export OMPI_MCA_pml='ucx'
export OMPI_MCA_mtl='^ofi'

Todo

• Add symmetry support for the QM calculations. Done.
• Add support for distance restraints. Done.
• Add support for angle restraints. Done.
• Refactor code to be OOP.
• Turn QRef into a proper restraint_manager class.

Citation

Lundgren, K. J. M., Caldararu, O., Oksanen, E., & Ryde, U. (2024). "Quantum refinement in real and reciprocal space using the Phenix and ORCA software", IUCrJ, 11(6), 921-937.
doi.org/10.1107/S2052252524008406