HOWTO for qmdff
===============

-------------------------------------------------------
---------------------- TURBOMOLE ----------------------
-------------------------------------------------------

1. Prepare a coord file.

2. Optimize.

3. Calculate the hessian (aforce, snf, driver [e.g. dftb or mopac])
   the name should be <hessian> (<hessian_driver> for driver runs)      

4. Do an ORCA (V.>2.9)  single point with for example
! pbe cosmo
! sv(p)
! grid2
! loosescf
! nofinalgrid
%output
       print[P_Hirshfeld] 1
       print[P_BondOrder_M] 1
       print[P_Mayer] 1
 end 

to get the charges and WBO (any other reasonable functional/basis will work as well). Large
charged (bio)systems can be problematic in the SCF and in these cases pbe/cosmo can be
replaced by e.g. pbe0. The following script does the single-point calculation and makes 
the Hirshfeld charges:

#!/bin/bash
orca inp > orca.out
cat orca.out | sed -n '/HIRSHFELD/,/TOTAL[0-9. ]*/ p' | grep '^ \+[0-9]\+ [a-zA-Z]\+' | gawk '{print $3}' > charges

The output must be on <orca.out> (read by qmdff). If for some reason ORCA does not work you
can take the charges from another source. The format of the file <charges> is just
#atoms lines with the charge of the atom in each line.

--------------------------------------------------
---------------------- ORCA ----------------------
--------------------------------------------------

1. Add to the above ORCA input the keyword
! freq
to get additionally the Hessian. The frequency and charges/WBO generation steps of course
can be separated and done at different levels. The code always tries to read <orca.out> 
and <my_orca_name.hess> (i.e. if you run ORCA with orca <my_orca_name>).

------------------------------------------------------
---------------------- Gaussian ----------------------
------------------------------------------------------

1. Run Gaussian with the following single-point(!) input/settings (of course
on the optimized geometry):

... freq pop(hirshfeld) density=current iop(7/33=1) iop(6/80=1) ...

-------------------------------------------------------------
------------------- running the FF generator ----------------
-------------------------------------------------------------

0. IMPORTANT IMPORTANT IMPORTANT IMPORTANT IMPORTANT IMPORTANT
 Do not forget to
 set ulimit -s unlimited 
in your shell! 

Use setenv OMP_NUM_THREADS <nproc> for parallel runs.
System charges must be specified in a file <.CHRG> with one line
and the appropriate system charge (used for checking the input
charges only).


1. Run it:

qmdff <QM_input_file_name> [-options] > out

Here 
QM_input_file_name=coord (TURBOMOLE) 
                  =my_orca_file_name.hess (ORCA) 
                  =my_Gaussian_name.out.
The names with extensions (coord, .hess, .out) are mandatory for automatic file detection.
The ORCA output should be in <orca.out>.
The final FF is in readable format in the file <solvent>.


2. Test the FF by doing a 50 ps MD:

qmdff <QM_input_file_name> -rd -md > out

(followed by molden qmdff.trj)


3. Optionally: test the FF by doing an optimization with the FF:

qmdff <QM_input_file_name> -rd -opt > out

(followed by molden qmdff.opt)
The optimizer is rather primitive at the moment and not intended for routine use.

-------------------------------------------------------------
------------------- A (complicated) example  ----------------
-------------------------------------------------------------

1. Generate the following orca input file (<inp> in this case):

! pbe SV(P) grid4 tightscf opt freq
* xyz 0 1 
   F         -1.95690068      0.20534125      1.19397795
   B         -0.60384932      0.17936152      1.14699629
   F          0.02030072      1.35150674      1.39207153
   F         -0.02279995     -0.91855264      1.69643953
   C          1.10757233     -0.13087716     -0.89274184
   O          1.67763054      0.74371187     -1.45632138
   N         -0.29829481     -0.08424579     -0.61283845
   H          1.57716236     -1.04280700     -0.50181492
   H         -0.75590393     -0.96831284     -0.82642105
   H         -0.74491725      0.66487405     -1.13934766
*
%output
  Print[ P_Hirshfeld ] 1 
end 


2. run ORCA: 
orca inp > orca.out

3. create the charge file:
cat orca.out | sed -n '/HIRSHFELD/,/TOTAL[0-9. ]*/ p' | grep '^ \+[0-9]\+ [a-zA-Z]\+' | gawk '{print $3}' > charges

4. run the generator:
qmdff inp.hess > qmdff.out

5. You will find a relatively large imaginary mode in the QMDFF Hessian. If you would try to
optimize the system for testing purposes with:
qmdff inp.hess -opt > qmdff.out
you will find that the BF3 group dissociates from the N-atom in contradiction to the DFT
input geometry. The reason can be found by careful inspection of the output. The B-N 
WBO is only about 0.4 which is a borderline case for covalent bonding. The qmdff program
internally has a threshold when a bond is considered to be covalent which is 0.5 by
default. Hence, in this example the amide and BF3 are considered as being only
non-covalently interacting and therefore the generated FF corresponds to something
different from the input. Such cases can not be handled fully automatically
and require chemical knowledge of the user. The problem is solved by calling the
program with a different threshold e.g.:

6. qmdff inp.hess -wbothr 0.4 > qmdff.out
will produce a reaonable FF with the weak but covalent B-N stretch term.
Similar problems normally do not appear for TM complexes because metal
atoms are internally handled somewhat differently.