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 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.

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
       
   • 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 
! 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
     • refinement.pdb_interpretation.restraints_library.mcl = False
     • refinement.pdb_interpretation.restraints_library.cis_pro_eh99 = False
     • refinement.pdb_interpretation.sort_atoms = 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.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
    

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, compared 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).
• 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).
  
• 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.

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