How to start a quantum-refinement job

The general procedure to set up a quantum refinement job consists of the following steps:

1. phenix.fetch_pdb <pdb_code> --mtz
   

2. Build a model.

   • The model in the QM regions needs to make chemical sense. This  means that the QM regions should be protonated and complete (but the model outside of the QM region can be incomplete, e.g. with no H atoms).
   • phenix.ready_set add_h_to_water=True optimise_final_geometry_of_hydrogens=False <pdb_file> can be useful for this purpose.
   • Instructions how to build a proper QM region can be found in standard QM/MM reviews, e.g.  10.1016/bs.mie.2016.05.014, 10.1002/anie.200802019 and 10.1021/cr5004419.

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.
   • Cut it out from the full file
   • phenix.elbow --final_geometry=<pdb_file>
   • The output if the .cif file
     

4. Prepare syst1 files that define the QM regions.

   • If there is only one QM region, the file should be named syst1. It there are multiple QM regions, the names should be syst11, syst12, etc.
   • Atoms included in the QM regions are 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 the script sort_pdb.py is also provided with QRef; this script takes as input the name of a PDB file and outputs a new PDB file with the suffix _sorted.pdb where the atoms are sorted in the same order as that  Phenix uses internally. You should  use this file as the input model for refinement, as well as the reference when defining the syst1 files.
       sort_pdb.py <model>.pdb
       Gives a file <model>_sorted.pdb
       
   • The syst1 files allows for multiple atoms or intervals of atoms to be specified on a single line, where ,or blank works as delimiters and- is used to indicate an interval.
   • # and ! can be used to include comments in the syst1 files.
   • A 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.
   • Sample syst1 file.

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  bonds for the link-atoms. If another junctfactor file is to be used this can be specified with the -j or --junctfactor option. Our standard file is available on GitHub.
     

   • The theory used for the ideal QM distance can be changed with the -l or --ltype option. Default is type 12. Current 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.
     

   • Parametrisation is done by simply optimising a small model of the QM system with the QM method and then insert the optimised bond length in the juncfactor file:
     
     CYM CB   CA   H1  # Cym CB junction
      7   1.1137d0     # MeS-, Cs, BP/RI, 6-31G*, 29/10-02
      9   1.1307d0     # MeS-, Cs, BP/RI, def2-SV(P), UR 24/6-08
     10   1.1306d0     # MeS-, Cs, PBE/RI, def2-SV(P), H in plane, UR 22/6-11
     11   1.1201d0     # MeS-, Cs, B3LYP/RI, def2-SV(P), UR 18/11-09
     12   1.1238d0     # MeS-, Cs, TPSS/RI, def2-SV(P), UR 29/12-09
     
     The first line contains the residue name (CYM), the QM atom bound to the link atom (CB) and the link atom (CA), i.e. the atom that is converted to a H atom in the QM calculation. The fourth entry is the desired atom type of the link atom in the MM force filed. It is not used in this implementation.
     The coming lines give the theory (ltype) and the optiised bond length in Å. Characters after "#" are comments and are ignored. The entry ends with an empty line.
     

   • 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 2 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 (1 + 2 n_syst1) files as well as selection strings:

     • 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).
     • qm_i_c.pdb, which is the model used to calculate the MM1 energy.
     • qm_i_h.pdb, which is the model used to calculate the QM1 energy.
       

     The output PDB files can 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 to select the model to refine when crafting the input to either phenix.refine or phenix.real_space_refine, see point 6 below.

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

6. Check that you have got the correct path to orca in file qref.dat
   Note that if you run Orca in parallel, it needs to be compiled with static linking
   

7. The next step is to prepare input files for Orca.
   

   • 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 (note that this overrules any Opt keyword, so that only a gradient setp is done)
! 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.
   • Sample input file. More examples can be found in the examples folder on GitHub.
   • Tutorial how to set up QM calculations.
     
8. Enable Phenix and Orca, e.g. by (on Cosmos)
   source /lunarc/nobackup/projects/teobio/Bio/Qref/phenix-1.20.1-4487-qref/phenix_env.sh
   

9. Now, you can start the quantum refinement. However, it is recommended to first create an empty file named qm.lock (this disables QR; touch qm.lock), then start a refinement so that an .eff file is obtained:
   phenix.refine <model>_sorted.pdb <model>.mtz
   The .eff file contains all the Phenix refinement settings for the current experimental data and model). Rename this file to phenix.params and edit it to make sure the following options are set:

   • For reciprocal space refinement (phenix.refine):
     • refinement.pdb_interpretation.restraints_library.cdl = False
(search for cdl)

     • refinement.pdb_interpretation.restraints_library.mcl = False
     • refinement.pdb_interpretation.restraints_library.cis_pro_eh99 = False
(usually already disabled; do the master key)

     • 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>
(from the qref_prep.py script)

     • refinement.hydrogens.refine = *individual riding Auto
     • refinement.hydrogens.real_space_optimize_x_h_orientation = False
     • refinement.main.nqh_flips = False
   • For neutron data, also
         scattering_table = wk1995 it1992 n_gaussian electron *neutron
     
   • 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.mdl = False
     • pdb_interpretation.restraints_library.cis_pro_eh99 = 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
   • Other options can be set as needed.

10. To run the quantum refinement job, delete the qm.lock file and 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.
    phenix.refine phenix.params --overwrite
    or
    phenix_real_space.refine phenix.params --overwrite
    If the job is run in parallel, you need to load modules for openMPI
    
    If you need to do a restart, do
    rm *.lock

1. In the directory where you have run the QRef refinement, run
   prep_geo_opt_qm_constrained.py <pdb_file_qm_refined>
   
   
2. It creeates the folder
   geo_opt_qm_constrained_1
   with files needed for the Orca optimisation
3. Copy or construct a sbatch file to run the orca optimisation and submit it
   

Calculating RSCC

1. With sort_atoms=False Phenix does not translate the molecule back to the original fram (i.e. in agreement with the map). Therefore:
   At the end of the slurm file there is a command
   phenix.pdbtools translate="x y z"
   
   
2. Run this command, with the final pdb file appended xxx_sorted_realspace_refined_001.pdb.
   
   
3. Run
   phenix.map_model_cc translated_pdb_file map resolution=xxx

Sample syst1 file

Sample ORCA input file

! TPSS D4 DEF2-SV(P)
! ENGRAD
%PAL NPROCS 4 END
*pdbfile +1 1 qm_1_h.pdb

Input files

• syst1 defines the QM system with atom numbers corresponding to the input pdb file xxxx_sorted.pdb
• qref.dat defines the junctions
• qm_1.inp is the Orca input file, showing the charge, multiplicity and the QM method. If there is more than one QM region, the corresponding files are called qm_2.inp, qm_3.inp, etc.
  

Where to find final results

• qr.log shows the QM energy in the various iterations. The last line is the final QM energy
  
• qm_1.out shows that last QM calculation. It includes charges and spin densities.
  
• grep 'FINAL SINGLE' qm_1.out     
  gives the final QM energy (of the first QM region).
• for x in A B ; do
  cd $x/geo_opt_constrained_qm_1
  grep 'FINAL SINGLE' qm_1.out | tail -n1
  cd ../..
  done
  To get strain energies
  
• qm_1_h.pdb is the final input coordinates of the QM system to ORCA. It is not necessarily the final coordinates, because the refinement might discard some updates of the coordinates. qm_1_c.pdb contains the corresponding coordinates for the MM1 calculation (with C-link atoms, instead of H-link atoms). If there are several QM regions, there are more files, like qm_2_h.pdb and qm_2_c.pdb.
  
• xxxx_refine_001.pdb contains the final coordinates. It contains also crystallographic information, e.g. final R-factors
  
• xxxx_refine_001.log is the log file from Phenix. It contains the largest residuals for the force field, the scale factors used. It shows the final R-factors at the end. grep "wxc = " *.log
  shows the weight factors (reciprocal) or (real-space):
  grep "overall best weight" *.log
  
• All residuals for the MM terms are given in xxxx_refine_001.geo.
  This is by default for the starting model. If you want it for the final model, set
  write_final_geo_file True
  in the Phenix parameter file.
  
• slurm*.out contains essentially the information in the previous file, but also some information about the minimisation:
• grep line slurm*.out
• grep 'HIS A  26' edstats.out
  The RSZD score is in the 6th last column (ZDa; before the negative one); the RSR result is in the Ra column (first of 5 results with 3 decimals at the end of the output); the RSCC score is in the CCSs column (second last of the middle group of results with 3 decimals; CCSm for water and metals).
• Molden often refuses to show the ORCA pdb files.
  Changepdb can fix this:
  rt (remove title)
  rn (set residue name to ???; not needed)
  f 1 2 (set the atom format to 2; not needed)
  w (overwrite the file)
  in fact, it is only the first 6 lines (ORIGIN and SCALE) that need to be deleted
  
• This file is not available in the GitHub implementation: file lbfgs.log contains some information from the LBFGS minimisations, but it is unclear to what use. It is typically >800 minimisations and >10000 steps.

Additional comments

Use the same altloc for a specific occupancy (also, qref_prep.py does give a warning if a syst1 definition encompasses different altlocs

Occupancies:
Kristoffer 16/8-24: I have now added change_occ_pdb.py as part of the QRef repository (https://github.com/krlun/QRef/blob/main/scripts/change_occ_pdb.py), which might be useful in terms of adjusting occupancies in a model. It now also supports using a list (using columns 17-26 in a PDB file).

As part of the rewrite now one of the options --syst1or --residues has to be specified (and they are mutually exclusive), i.e. to use a syst1 file the syntax is

$ change_occ_pdb.py pdb occupancy --syst1 syst1_file

or using a residue list

$ change_occ_pdb.py pdb occupancy --residues residue_file

One might also have to manually update the shebang at line 1 so that it points towards an existing Python interpreter.

Do note that this rudimentary implementation is kind of a “all or only syst1” version; I recommend visually inspecting the starting model to find where different altlocs merge, then change the occupancies for the altlocs for the entire part of the subsystem to a single value as to not introduce weird artefacts from the experimental gradients (i.e. construct a complementary syst1 file (which does not need to contain a specification of the junction atoms) covering the entire altloc part of the model and use this new file to change the occupancies).

Also, it is possible to change occupancies using phenix.pdbtools, see the documentation at https://phenix-online.org/documentation/reference/pdbtools.html.
It is also probably doable through Coot, using the “Measures” -> “Residue Info…” tool.

Installation

A somewhat more hands-on description of the installation than described on GitHub.

1. Download Phenix or make a copy of a existing installation of Phenix.
   If moving an already installed version, change the path in the phenix_env.sh file to the new location.
   
2. Download QRef from GitHub (Code - Download zip).
3. Move it to a proper directiory and
   unzip QRef-main.zip
4. Alternatively get the link from GitHub (Code - Right-click Download zip - Copy Link) and then
   git clone <link>
5. source path_to_phenix/phenix_env.sh
   
6. Go to the QRef-main/modules directory
7. cp -pr qref $PHENIX/modules
8. cp cctbx/geometry_restraints/energies.py $PHENIX/modules//cctbx_project/cctbx/geometry_restraints
9. cp cctbx_project/mmtbx/model/model.py $PHENIX/modules/cctbx_project/mmtbx/model
10. Uncomment the line
    # if not os.path.exists('qm.lock') and (os.path.exists('xyz_reciprocal.lock') or os.path.exists('xyz.lock')):
    and comment out the next line (if not os.path.exists('qm.lock'):)
    in the file $PHENIX/modules/cctbx_project/cctbx/geometry_restraints/energies.py
11. Add the line
    import os
    in the beginning of the file
    $PHENIX/modules/phenix/phenix/refinement/xyz_reciprocal_space.py
    Also add the two lines
          with open('xyz_reciprocal.lock', 'w'):
            pass
    before the line containing mmtbx.refinement.minimization.lbfgs
    and this line after it:
          os.remove('xyz_reciprocal.lock')
    
12. Do similar changes in the file $PHENIX/modules/phenix/phenix/refinement/macro_cycle_real_space.py:
    Add the line
    import os
    in the beginning of the file.
    Add the two lines
        with open('xyz.lock', 'w'):
          pass
    before the lines containing either self.minimization_no_ncs() or self.minimization_ncs(), i.e. two places
    and this line after them:
        os.remove('xyz.lock')
13. which cctbx.python
14. Copy the answer of this command to the first line of the files
    path_to_Qref-main/scripts/qref_prep.py and sort_pdb.py
15. Ensure that the path_to_Qref-main/scripts directory is in you PATH variable.
    export PATH=$PATH:path_to_Qref-main/scripts