How to start a quantum-refinement job

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

1. Ensure that you have proper paths to Phenix and Orca:
   export PHENIX=/lunarc/nobackup/projects/teobio/Bio/Qref/phenix-1.21.2-5419
   export PHENIX_VERSION=1.21.2-5419
   . $PHENIX/build/setpaths.sh
   module add GCC/12.3.0  OpenMPI/4.1.5 ORCA/6.0.1
   
   or on Tetralith:
   export PHENIX=/proj/teobio/Gayathri/softwares/phenix/phenix-installer-1.21.1-5286-intel-linux-2.6-x86_64-centos6/phenix-1.21.1-5286
   export PHENIX_VERSION=1.21.1-5286
   . $PHENIX/build/setpaths.sh
   module load ORCA/6.1.0-bdist
   
   
2. Get the deposited pdb file and structure factor files (the latter converted to mtz format):
   phenix.fetch_pdb <pdb_code> --mtz >fetch.out
   If you do not get the mtz file, download the structure factor file from PDB and then
   phenix.cif_as_mtz <structure_facto_file>
   The latter often gives problem. Try to solve them with the instructions given from Phenix.
   However, Phenix can run also on the cif file.
   

3. Protonate the pdb file and construct restraints files for all non-standard residues.
phenix.ready_set add_h_to_water=True optimise_final_geometry_of_hydrogens=False <pdb_file> >ready_set.out
   The file <pdb_code>.ligands.cif contains all restraints needed.
   Move only protons of the QM system to the original PDB file.
   

4. Already done in previous step but may be needed if the step fails.
   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 out the non-standard ligand from the full pdb file and put it into a separate pdb file.
   phenix.elbow --final_geometry=<pdb_file> >elbow.out
   The output is in the xxx.cif file
   

5. Sort the pdb file (so that the order does not change when running Phenix.
   sort_pdb.py <model>.pdb
   This gives a file <model>_sorted.pdb.
   

6. Prepare the syst1 file that defines the QM region(s).
   The QM region needs to make chemical sense. This means that the QM region should be protonated and complete (but the model outside of the QM region can be incomplete, e.g. with no H atoms).
   Some instructions how to select the QM regions are given below.
   

   • 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.
     
   • 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 file  (after a blank line) indicates 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.

7. 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. Note that the pdb file name must come first, before the options.
   
   qref_prep.py6qi1_sorted.pdb -s syst11 syst12 -j /lunarc/nobackup/projects/teobio/Bio/Data/Comqum/junctfactor -c cif >prep.out


   • If you have several syst1 files, their names need to be specified with the -s option:
     -s syst11 syst12
     

   • The junctfactor file needs to be present in the same directory as where qref_prep.py is run. The junctfactor file contains 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. For example on gefion or Cosmos:
     -j /temp4/bio/Data/Comqum//junctfactor
     -j /lunarc/nobackup/projects/teobio/Bio/Data/Comqum/junctfactor
     Our standard file is available on GitHub.
     Some instructions how to construct new entries are found below.

   • The theory used for the ideal QM distance can be changed with the -l or --ltype option.
     The default is 12, TPSS/def2-SV(P). Other options are listed below.

   • If you prepared restraint files for novel residues in point 3 above, qref_prep.py needs to be made aware of these. This can be achieved with the -c or --cif option, e.g.
     qref_prep.py 5b6w_sorted.pdb -c 5b6w.ligands.cif
     

   • 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.
       Note that this file keeps the atom names of the CL atoms (for technical reasons; uses cbx routines) and changes only the element column to H. This may confuse some viewers. It also contains anisotropic B-factors, not really used by QM software.
       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.

   • • 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. Experience has shown that the force constant needs to be $\ge$ 2500 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. Experience has shown that the force constant needs to be $\ge$ 10 (?) to achieve adherence to the restraint.

     • 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 the -h or --help option.

8. 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 (only Orca 5.2):
   /lunarc/nobackup/projects/teobio/Kristoffer/orca/orca_5_0_4-static/orca
   sed -i s#/usr/bin/orca#/lunarc/nobackup/projects/teobio/Kristoffer/orca/orca_5_0_4-static/orca# qref.dat
   

9. 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.
     ! TPSS D4 DEF2-SV(P)
     ! ENGRAD
     %PAL NPROCS 4 END
     *pdbfile +1 1 qm_1_h.pdb
   • Tutorial how to set up QM calculations.
     

10. Run a standard (non-QM) Phenix refinement
    touch qm.lock
    phenix.refine --overwrite <model>_sorted.pdb <model>.mtz (cif_fil)
    Find the wxc factor used:
    grep "wxc = " *.log
    | wxc = 23.293294             wxc_scale = 0.500            wc = 1.000         |
    | wxc = 3.484164              wxc_scale = 0.500            wc = 1.000         |
    | wxc = 3.656784              wxc_scale = 0.500            wc = 1.000         |
    and calculate the average of the latter two. Note that it should be divided by two, because wxc_scale=0.5.
    
    

11. Construct a  Phenix parameter file.
    This is the .eff file from the Phenix refinement in the previous point. It contains all the Phenix refinement settings for the current experimental data and model).
    Rename the .eff file to phenix.params
    cp 5b6w_sorted_refine_001.eff phenix.params
    and edit it to make sure the following options are set:
    Change pdb and mtz file names if needed
    

    • For reciprocal space refinement (phenix.refine):
      sed -i "s/\*individual_sites \*individual_sites_real_space/\*individual_sites individual_sites_real_space/" phenix.params
      sed -i "s/\*individual_adp group_adp tls \*occupancies/\*individual_adp group_adp tls occupancies/" phenix.params
      sed -i "s/nqh_flips = True/nqh_flips = False/" phenix.params
      sed -i "s/individual \*riding Auto/\*individual riding Auto/" phenix.params
      sed -i "s/real_space_optimize_x_h_orientation = True/real_space_optimize_x_h_orientation = False/" phenix.params
      sed -i "s/cdl = True/cdl = False/" phenix.params
      sed -i "s/mcl = True/mcl = False/" phenix.params
      sed -i "s/cis_pro_eh99 = True/cis_pro_eh99 = False/" phenix.params
      sed -i "s/sort_atoms = True/sort_atoms = False/" phenix.params
      sed -i "s/flip_symmetric_amino_acids = True/flip_symmetric_amino_acids = False/" phenix.params
      sed -i "s/wxc_scale = 0.5/wxc_scale = 1.0/" phenix.params
      
      # Insert proper values here
      sed -i "s#6qi1.mtz#/lunarc/nobackup/projects/teobio/Ulf/Time_resolv/Facd/12/6qi1.mtz#" phenix.params
      sed -i "s/6qi1.pdb/6qi1_sorted.pdb/" phenix.params
      sed -i "s/6qi1_refine/6qi1_sorted_refine/" phenix.params
      sed -i "s/fix_wxc = None/fix_wxc = 3.91/" phenix.params
      grep 'enabled = False' phenix.params
      
      # Only this needs to be done by hand individual = None (insert string from prep command)
      
      
      • Lines 242–245:
refine {
    strategy = *individual_sites individual_sites_real_space rigid_body \
               *individual_adp group_adp tls occupancies group_anomalous den
    sites {
      individual = <reciprocal selection string> (string from the qref_prep.py script)
      torsion_angles = None
e.g.
      individual = "((resseq 216 or resseq 300) and chain A)"
      • Line 282:
refinement.main.nqh_flips = False
  main {
    show_process_info = False
    stop_for_unknowns = True
    bulk_solvent_and_scale = True
    apply_overall_isotropic_scale_to_adp = True
    flip_peptides = False
    nqh_flips = False
      • Lines 438–443:
hydrogens {
    refine = *individual riding Auto
    force_riding_adp = None
    optimize_scattering_contribution = True
    contribute_to_f_calc = True
    high_resolution_limit_to_include_scattering_from_h = 1.6
    real_space_optimize_x_h_orientation = False
      • Lines 585–600:
pdb_interpretation {
    restraints_library {
      cdl = False
      mcl = False
      cis_pro_eh99 = False

      • ~12 lines further down
    }
    sort_atoms = False
    flip_symmetric_amino_acids = False
      • ~79 lines further down:
      max_representative_chains = 100
      enabled = False   (normally already set to False)
    }
    c_beta_restraints = True
        
      • With quantum refinement, you should always run with a fixed wxc factor. We normally use the factor suggested by Phenix (obtained in the previous point).
Lines 1095–1100:
  target_weights {
    optimize_xyz_weight = False
    optimize_adp_weight = False
    wxc_scale = 1
    wxu_scale = 1
    wc = 1
    wu = 1
    fix_wxc = 1.80
    fix_wxu = None

    • 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="all and not <reisdues with any atom in the QM system>"

      • pdb_interpretation.reference_coordinate_restraints.sigma = 0.01
      • refinement.refine.sites.individual = <reciprocal selection string>
(from the qref_prep.py script)
Note that it is opposite to what is used for reciprocal-space refinement, something like:
        individual = "all and not (((resseq 133 or resseq 232 or resseq 233 or resseq 237 or resseq 282 or resseq 283 or resseq 376 or resseq 384 or resseq 449) and chain A) or ((resseq 601 or resseq 603 or resseq 614) and chain A) or ((resseq 783) and chain A))"
      • To run with a fixed weight for restraints (should normally always be done):
        weight = 10
    • Other options can be set as needed.
    • Check what non-default parameters you have in your parameter file:
      phenix.refine --diff_params "phenix.params"

12. 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 (those needed for ORCA). Phenix does not run in parallel.
    
    If you need to do a restart, do
    rm *.lock

13. Set up the strain-energies with the command
    prep_geo_opt_qm_constrained.py <pdb_file_qm_refined>
    Go to the directories geo_opt_qm_constrained_1, copy a sbatch script and run them
    
14. Run edstats to get the RSZD scores:
    First run
    phenix.maps
    without any arguments to get the file maps.params.
    
    Insert the final pdb file and the reflection file into the file maps.params:
       input {
         pdb_file_name = "5b6w_sorted_refine_001.pdb"
         reflection_data {
           file_name = "5b6w_sorted_refine_001.mtz"
    It works with the original xxxx.mtz file also (they normally give the same results), but it is more safe with the one from Phenix (ensures that all fields are the same as used for the refinement).
    
    Then run
    phenix.maps maps.params
    
    Ensure that you can run ccp4 programs:
    . /lunarc/nobackup/projects/bio/CCP4/ccp4-7.0/bin/ccp4.setup-sh
    Construct a file edstats.sh and run it
    ./edstats.sh
    This is the content of the edstats.sh file (file names in bold face may need to be changed):
    cad hklin1 pdb_sorted_refine_001_fmodel.mtz hklin2 pdb_sorted_refine_001_map_coeffs.mtz hklout hkl_ed.mtz<<EOF
    LABIN FILE 1 E1=FOBS E2=SIGFOBS E3=FCALC E4=PHIFCALC E5=FOM E6=R_FREE_FLAGS
    LABOUT FILE 1 E1=FOBS E2=SIGFOBS E3=FC E4=PHIC E5=FOM E6=FREE
    LABIN FILE 2 ALL
    EOF
    edstats.pl -xyzin pdb_sorted_refine_001.pdb -hklin hkl_ed.mtz -flabel FOBS -noerror
    
    This gives the file edstats.out
    It contains results from all residues.
    You normally want to grep out only the interesting results, e.g.
    grep "CYS C   62" edstats.out
    grep "GLY C   87" edstats.out
    ...
    
    There are three sections in the output: one for the mainchain (m), one for the sidechain (s) and one for all atoms in the residue (a).
    The interesting results are in the ZD-, ZD+ and ZD columns (maximum value for the negative or positive difference density, as well as the maximum absolute value of these two.
    
    If you run edstats.py without any options, the program will tell you the content of all columns.
    
    
15. Remove unnecessary files for a restart calculation:
    \rm *_sorted.pdb *_sorted_refine_001.cif *_sorted_refine_001.eff *_sorted_refine_001.geo *_sorted_refine_001.log *_sorted_refine_001.mtz *_sorted_refine_001_2mFo-DFc_map.ccp4 *_sorted_refine_001_fmodel.mtz *_sorted_refine_001_map_coeffs.mtz *_sorted_refine_002.def edstats.out hkl_ed.mtz qm_?.bibtex qm_?.densities qm_?.densitiesinfo qm_?.engrad qm_?.gbw qm_?.ges qm_?.out qm_?.pdb qm_?.property.txt qref_?.log settings.pickle slurm-*.out 
    You also need to rename the *_sorted_refine_001.pdb file to *_sorted.pdb
16. This is a script to set occupancies in a restart calculation:
    https://github.com/krlun/QRef/blob/main/scripts/change_occ_pdb.py
    
    
17. QRef is completely reproducible, meaning that when started with the same files, the same results are obtained.
    

• Harmonic 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. Experience has shown that the force constant needs to be $\ge$ 2500 to achieve adherence to the restraint.
  or added in the qref.dat file:
          "restraints_distance": [
              [
                  1,
                  20912,
                  3.12,
                  2500.0
              ]
          ],
  

• Harmonic 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. Experience has shown that the force constant needs to be $\ge$ 10 (?) to achieve adherence to the restraint.

• Groups in symmetry-related copies can be included in the QM system by adding the symmetry operations needed in the qref.dat file:
          "transforms": [
              {
                  "R": [
                      [
                          0.5,
                          -0.86601334,
                          0.0
                      ],
                      [
                          -0.86603747,
                          -0.5,
                          0.0
                      ],
                      [
                          0.0,
                          0.0,
                          -1.0
                      ]
                  ],
                  "atoms": "2534-2541, 2546-2550",
                  "t": [
                      40.66363045,
                      70.43245528,
                      202.4127601
                  ]
              }
        ]
  
  

1. Set up the syst1 file(s)
2. Protonate the file, if necessary.
3. Run qref_prep.py:
qref_prep.py<file>.pdb -s syst11 syst12 -j /lunarc/nobackup/projects/teobio/Bio/Data/Comqum/junctfactor -c cif >prep.out
4. mkdir Qm1 Qm2
   cp qm_1_h.pdb Qm1
   cp qm_2_h.pdb Qm2
5. Set up the following qm_1.inp file (and similar for Qm2):
   ! TPSS D4 DEF2-SV(P) Opt pal8
   *pdbfile +1 1 qm_1_h.pdb
   %geom optimizehydrogens true
   end
6. Start the calculations with something like:
   #!/bin/sh
   #SBATCH -n 8
   #SBATCH -t 48:00:00
   #SBATCH -A LU2025-2-123
   
   module add GCC/12.3.0  OpenMPI/4.1.5
   alias orc=/lunarc/nobackup/projects/teobio/Bio/ORCA/orca_6_1_0_linux_x86-64_shared_openmpi418/orca
   
   orc qm_1.inp >qm_1.out
   cd ../Qm2
   orc qm_2.inp >qm_2.out
7. Results are found by
   grep FINAL *.out
   
   

Sometimes when dealing with QR structures (especially for time-resolved structures), it can be desirable to modify the occupancies within the QM region. Then, we recommend the following procedure:

1. Start with the best QR structure
2. Perform a normal Phenix (non-QM) refinement of occupancies:
3. Copy to a separate directory
4. Remove all QR files:
   \rm mm* qm* qref* syst* slurm*
5. mv 5b6w_sorted_refine_001.pdb 5b6w_sorted.pdb
6. \rm 5b6w_sorted_refine_001* 5b6w_sorted_refine_002* edstats.out settings.pickle hkl_ed.mtz
   
7. In phenix.param, do the following changes:
8.   refine {
       strategy = individual_sites individual_sites_real_space rigid_body \
                  individual_adp group_adp tls *occupancies group_anomalous den
   i.e. turn off coordinate and adp refinement and on occupancy refinement.
9. Next lines:
      sites {
         individual = None
         torsion_angles = None
         rigid_body = None
       }
     i.e. delete any string after individual (probably strictly not needed) 
   
10.     occupancies {
          individual = None
          remove_selection = None
          constrained_group {
            selection = "(resseq 110) and (chain A) and (altloc A)"
            selection = "(resseq 110) and (chain A) and (altloc B)"
            selection = "(resseq 110 or resseq 401) and (chain A) and (altloc C)"
            selection = "(resseq 110 or resseq 401) and (chain A) and (altloc D)"
          }
        }
    i.e. define the various altloc groups, one for each QM region
11. If you have H atoms in the structure (outside QM region) you also need change from individual to riding (otherwise all of them will get non-unity occupancies):
      hydrogens {
        refine = individual *riding Auto
    
12. Modify the sbatch file so that only Phenix is run
    phenix.refine phenix.params --overwrite
    and then submit it.
    
    
    

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 frame (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

Restart a QR calculation from an old QR result

1. Copy files
   cp phenix.params syst1 qref.dat qm_1.inp qm_1_h.pdb mm_1_c.pdb *refined_000.pdb q ../New
   
1. phenix.params
2. start pdb file
3. syst1
   
4. qref.dat
5. qm_1.inp
6. qm_1_h.pdb
7. mm_1_c.pdb
8. slurm file
9. experimental data (normally specified in phenix.param and no need to copy)
10. Possibly, xxx.cif files
    
2. Change in phenix.param the name of the pdb file in two places
3. Change output_serial = 1 in phenix.param
4. Or alternatively, change the name of the pdb file to the original input-file name (xxx_sorted.pdb); then there is no need to change phenix.param, but it is a bit harder to know that it is a restart and where it was started from.
   
5. Change in qm_1.inp (charge and multiplicity)

Calculate single-point Phenix on deposited structure
phenix.model_vs_data pdb mtz
Calculates statistics of the current model, e.g. Rfactors

Sample syst1 file

Sample ORCA input file

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

or for metal-containing systems:

Input files

1. phenix.param (parameter file for phenix)
2. The start model (a pdb or mmcif file)
   
3. syst1 defines the QM system with atom numbers corresponding to the input pdb file.
   
4. qref.dat defines the junctions (it also tells phenix to run QRef)
   
5. 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.
6. qm_1_h.pdb
7. mm_1_c.odb
8. A slurm script file
9. Experimental data (mtz or map file). The name is specified in the phenix.param file and it might be outside the current director (since it does not change, that is normally and preferably the case).
10. Possibly, additional restrants files for phenix (xxx.cif)

• .eff = the parameter file actually used (similar to phenix.param, but not exactly)
• .pdb = the final pdb file
• .cif = the final cif file (same info as pdb file)
• _initial.geo = the restraint (MM) file
• .log = the logt file
• settings.pickle - a byte-dump file of all settings to phenix.refine. Needed by QRef MM1 calculation. Can be deleted.
  

Orca gives 6 output files (per QM system) with names qm_1.*

• .ges = old gbw file
  
• .densities = unformatted file with densities
  
• .gbw = unformatted file will all Orca info, including wavefunction
  
• .engrad = formatted file with coordinates, energy and gradients
  
• _property.text = formatted file with various properties: energy, DFT components, Mayer population and bond order, DFT-D correction, dipole moment and coordinates
  
• .out = Orca single-point output file
• qm_1.pdb (coordinates of the QM system, but without much info)
  

QRef gives one additional file and updates two files

• qref_1.log: QM, MM12, MM1 and QM/MM energies in each step
• qm_1_h.pdb: QM system with HL atoms (but the name is still that of CL); input to Orca; Contains pdb entries and ANISOU
  
• mm_1_c.pdb: same file as qm_1_h.pdb, but with CL atoms (positions and element)

Where to find final results

• qr_1.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.
  However, geo_opt_qm_constrained/qm_1_h.pdb always contains the final coordinates (if Kristoffer's strain script has been used).
• xxxx_refine_001.pdb contains the final coordinates. It contains also crystallographic information, e.g. final R-factors
  grep 'R VALUE            (WORKING SET)' *_refine_001.pdb
  grep 'FREE R VALUE  ' *_refine_001.pdb
• 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.
  Start R-work = 0.1609, R-free = 0.1795
  Final R-work = 0.1561, R-free = 0.1781
  
• 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).
• Spin populations are found in qm_1.out.
  Search for "MULLIKEN ATOMIC CHARGES AND SPIN POPULATIONS"
  
• 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.

Instructions how to build the QM regions

• 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.
  
• In general:
• Include all atoms of chemical interest.
• Never truncate within aromatic systems.
• Include all ligands of metal sites.
• Preferably, cleave C–C single bonds.
• Include a methyl group to truncate side chains.
• Many residues can be truncated at the CA atom (converted to H). Larger residues may be cleaved further out in the side chain, e.g. CB for Glu and Gln, and CD for Lys.
• Backbone is best modelled by CH3CONH-CH2-CONH2CH3 units.
• Include residues hydrogen bonding to the system of interest, if possible.

Instructions regarding the junctfactor file

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

• A sample entry from the 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.

• 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)
  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:
  

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 (https://github.com/krlun/QRef.; 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. Overwrite the file $PHENIX/modules/cctbx_project/cctbx/geometry_restraints/energy.py with the proper file from cctbx_project/cctbx/geometry_restraints
   For example (for Phenix version 1.21.2-5419):
   cp cctbx_project/cctbx/geometry_restraints/energies_1.21.2-5419.py $PHENIX/modules/cctbx_project/cctbx/geometry_restraints/energies.py
9. Do the same for the model.py file.
   For example:
   cp cctbx_project/mmtbx/model/model_1.21.2-5419.py $PHENIX/modules/cctbx_project/mmtbx/model/model.py
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
          minimized = mmtbx.refinement.minimization.lbfgs(
    and this line after it:
          os.remove('xyz_reciprocal.lock')
    See below
    
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 both the line containing self.minimization_no_ncs() and the line containing self.minimization_ncs(), i.e. two places
    and this line after them:
        os.remove('xyz.lock')
    See below
    
13. Copy the answer of this command
    which cctbx.python
    to the first line of the files
    path_to_Qref-main/scripts/qref_prep.py, prep_geo_opt_qm_constrained.py and sort_pdb.py
    
    
14. Ensure that the path_to_Qref-main/scripts directory is in you PATH variable.
    export PATH=$PATH:path_to_Qref-main/scripts

Thus, point 11 should look like this:
      if(self.all_params.main.scattering_table == "neutron"):
        is_neutron_scat_table = True
      with open('xyz_reciprocal.lock', 'w'):
        pass
      minimized = mmtbx.refinement.minimization.lbfgs(
        restraints_manager       = self.model.restraints_manager,
        fmodels                  = self.fmodels,
        model                    = self.model,
        refine_xyz               = True,
        lbfgs_termination_params = lbfgs_termination_params,
        target_weights           = self.target_weights,
        verbose                  = 0,
        correct_special_position_tolerance =
          self.all_params.main.correct_special_position_tolerance,
        is_neutron_scat_table    = is_neutron_scat_table,
        qblib_params             = getattr(self.all_params, "qbio", None),
        macro_cycle              = self.macro_cycle,
        h_params                 = self.h_params,
        log                      = self.log)
      os.remove('xyz_reciprocal.lock')
    elif(self.all_params.ncs.constraints.apply_to_coordinates):

Likewise, point 12 should look like this:
    if(ncs_groups is None or len(ncs_groups)==0):
      with open('xyz.lock', 'w'):
        pass
      self.minimization_no_ncs()
      os.remove('xyz.lock')
    else:
      with open('xyz.lock', 'w'):
        pass
      self.minimization_ncs()
      os.remove('xyz.lock')

  def nqh_flips(self):

Template qref.dat file:
{
    "cif": [
        "HXG.cif"
    ],
    "ltype": 12,
    "n_atoms": 25259,
    "orca_binary": "/proj/teobio/users/x_gayyu/softwares/orca/orca_5_0_4-static/orca",
    "restart": null,
    "syst1": {
        "g": {
            "3": 0.7259016393442623,
            "27": 0.7308552631578947,
            "872": 0.7303273213092851,
            "892": 0.7303273213092851
        },
        "link_pairs": {
            "3": 2,
            "27": 28,
            "872": 873,
            "892": 893
        },
        "restraints_angle": [],
        "restraints_distance": [
            [
                1,
                20912,
                3.12,
                2500.0
            ]
        ],
        "transforms": []
    },
    "syst1_files": [
        "syst1"
    ],
    "w_qm": 7.5
}