ComQum-Pol-EXAFS

• QM+pEXAFS, i.e. QM suplemented by polEXAFS (or an automatic EXAFS-refinement protocol, where QM is used to suplement the EXAFS data to get a full set of coordinates, like MM is used in crystallographic refinement).
• QM/MM+pEXAFS, i.e. QM/MM suplemented by EXAFS inside a full protein.
  
• ComQum-X +pEXAFS, i.e. quantum-refinement suplemented by polEXAFS.

• U. Ryde, X. Li, P. E. M. Siegbahn, E. Sproviero, V. Batista (2010), to be submitted.

How to set up a ComQum-EXAFS calculation

1. Currently, the program is designed for exerimental data, provided only in R space.
   
   
2. Start from Turbomole files of the desired QM system.
   Run the program turbo2pfeff to get 3 feff[123].inp files.
   The program reads the file transform to perform some coordinate transformations; This is partly hardcoded now (some of the transformations are double transformations), so it will probably not fully work for other proteins.
   Currently, the parameters of feff are hard coded, so you need to change the code to modify the parameters.
   Consider to change the keywords TITLE, RPATH, NLEG, CRITERIA, and HOLE, (cf. the FEFF page).
   It is not possible calculate Debye-Waller factors. Instead, we use fixed factors.
   
   
3. Run cqpefixenergy to get all three polarised EXAFS calculations.
   This program will run feff83 three times. This is the only version of FEFF that seems to work for polarised EXAFS.
   The only output of feff83 needed are the xmu.dat file, which then is postprocessed by other programs (xmu2chi and ifeffit) to finally obtain a Chi^2 estimate and a R-factor.
   The file all.dat contains the three experimental and calculated spectra, ready to be plotted by gnuplot
   Template commands for gnuplot are:
   plot "all.dat" using 1:2 with lines,"all.dat" using 1:5 with lines
   plot "all.dat" using 1:3 with lines,"all.dat" using 1:6 with lines
   plot "all.dat" using 1:4 with lines,"all.dat" using 1:7 with lines
   
   To run this command, you need a comqum.dat file with the following lines:
   $weight_exafs
    1.0
   $weight_qm
    1.0
   $r_range
     1.0 4.3
   $e0
    6543.3
   $expfile
   pexafs_exp.txt
   
   $r_range gives rmin and rmax, used in the calculation of chi^2.
   $e0 is E0
   $expfile is the name of the file with the experimental data
   
   This is all needed to do single-point energy calculations
   
   
4. To run a geometry optimisation, you need to add some additional data to the file comqum.dat file (sample file).
   You should check the atoms for which you want to calculate the gradients in the keyword $exafs_gradients (see below for an explanation); by default internal gradients are calculated for all paths shorter than 3.0 Å.
   You may change $weight_exafs and perhaps also $weight_qm.
   
   You man also turn on the double-sided gradients (increases the time by a factor 2 in the EXAFS calculations):
   $double_sided_gradient on
   Finally, you may change the step size for the numerical EXAFS gradients (although the default is normally OK). However, note that you should not use values smaller than ~1E-4, because there are only five decimals in the xmu.dat file, so that you will loose precision and may get zero gradients.
   $g_step
    1.0D-4
   
   
5. Set up the Turbomole ComQum, or ComQum-X calculations as usual, and add the needed keywords to the comqum.dat file (explained above).
   
   In addition, you need the file with the exafs raw data, the transform file, a file process.iff with the following content (no changes should be needed to this file - it only converts the calculated data from k to R-space):
   
   # This is the file process.iff, set up for ComQum-PExafs calculations
   # It converts the calculated data from k-space to R-space.
   # Ulf Ryde, 15/3-10, based on templates from V. Batista and E. Sproviero.
   #
   read_data(file = fe_chi.DAT, type = chi, group = fe_chi)
   fek_chi.k = sqrt(0.2624683*fe_chi.k)
   newfe.k = range(0,ceil(fek_chi.k),0.05)
   newfe.chi = qinterp(fek_chi.k, fe_chi.chi, newfe.k)
   fftf(newfe.chi, k=newfe.k, kmin=floor(fek_chi.k), kmax=ceil(fek_chi.k), dk=1,kweight=3, kwindow='hanning')
   write_data(file=fe_FT.out,newfe.r,newfe.chir_mag)
   exit
   
   
6. It is normally wise to first run (in a separate directory)
   ridft>logd
   rdgrad>logg
   cqpefixenergy>loge
   cqpefixgrad>logi
   And then move to the original comqum.dat file
   $norm_scale
    0.67531351E-09 (or whatever)
   Otherwise, the first energy will be different from the following ones.
   
   
7. Start ComQum-E by
   for QM+pEXAFS
   jobexpe -backup -ri -c 200
   or for QM/MM+pEXAFS:
   comqumpe -backup -ri -c 200
   or for ComQum-X+pEXAFS:
   comqumxpe -backup -ri -c 200
   
   
8. When finished, the results are in the files energy (the first energy is the ComQum-EXAFS energy, the third is the Turbomole energy, and the fourth is chi_squared EXAFS energy) and all.out.
   $energy      SCF               SCFKIN            SCFPOT
        2 -4043.63646887853  4024.47271869400 -4043.63702859200 22325.04453929200
   
   
9.  If you have problem with convergence, you may:

• reduce dqmax to 0.03 (or smaller; instead of 0.3) to improve the convergence, especially if w_Exafs is high.
  $coordinateupdate
     dqmax=0.03
     interpolate  on
     statistics    5
• Use steepest descent
  $forceupdate
     none
     allow=0.0 threig=0.005  reseig=0.005  thrbig=3.0  scale=1.00  damping=0.0
• Restart the forceapprox file
  $forceinit on
     diag=default
  
  
  

1. Finally, you can remove unnecessary files by typing
   purtur
   

• jobexpe or comqumpe or comqumxpe
  
• cqefixenergy
• cqefixgrad
• turbo2pfeff
• xmu2chi