CamCASP/Programming/0

From CUC3
Jump to navigation Jump to search

CamCASP => Programming => Status

ENERGY-SCAN

What's working:

  • REDO_DF_ON_ROTATION = .TRUE.
    • E1elst, E1exch, E2ind(UC),E2exind(UC),E2exdisp(UC),E2disp(UC), Overlap
  • REDO_DF_ON_ROTATION = .FALSE.
    • E1elst, E1exch

What's not working and why:

  • REDO_DF_ON_ROTATION = .TRUE.
    • E2ind & E2disp: The DF-FDDS is created correctly for the first dimer configuration, but as it is not updated, it is incorrect for all subsequent configurations.
    • E2exind & E2exdisp: Since these are obtained by scaling E2exind(UC) & E2exdisp(UC) using E2ind & E2disp (and their UC counterparts), these are wrong too.

What's going on with the Hessians?

Here's what happens when the make_prop subroutine is called:

  1. make_prop
  2. densfit_prop
    1. init_prop
      1. df_monomer : Done each time.
      2. OVOV AAAA, VVOOAAAA, OVOVBBBB, VVOOBBBB : Done each time.
      3. make_h1_h2 : Done only once!
      4. 4 H2 D : Calculated only once.
  3. Solution of LR-DFT equations.
  4. DF-PROP done


  • So, if the molecule rotates then DF is repeated irrespective of REDO_DO_ON_ROTATE. This is an error.
  • Then, even if REDO_DF_ON_ROTATE=.true. and DF is correctly repeated, (4 H2 D) is calculated only once, so this will not be consistent with the new DF. This is an error.

If we want to get the Cartesian/REDO_DF_ON_ROTATE=.true. route working, then we have to repeat the calculation of (4 H2 D) every time the DF is repeated.

Actually, I don't think this is correct. Here's what's happenning:

  • The MOs and Hessians are calculated with the molecule in its reference geometry.
  • IF the density-fitting and DF-FDDS are calculated in this reference geometry, all will be well.
  • However, IF the molecule is rotated, and the MOs rotated and the DF is done in the rotated geometry, then, while the DF is correct, the rotated MOs are no longer consistent with the MOs used to calculate the Hessian integrals obtained from DALTON.
  • Therefore, we must conclude that the DF-FDDS can be calculated *only* in the REFERENCE molecular geometry.

How then are these DF-FDDSs to be used when the molecule rotates?

Water dimer

Sadlej/MC PBE0/AC

Work dir: /home/am592/DistProp/systems/water/scans/sadlej_mc

Units: kJ/mol

Water geometry: 

MOLECULE water1
  Units Bohr
  ! Vibrationally averaged geom.
  O1          8.0    0.00000000     0.00000000     0.00000000
  H1          1.0   -1.45365196     0.00000000    -1.12168732
  H2          1.0    1.45365196     0.00000000    -1.12168732
END

Dimer geometries used: 

     Rx      Ry      Rz      alpha   Nx      Ny      Nz       
    1   7.5489   0.0000   0.0000 180.0000   1.0000   0.0000   0.0000 <--water2 unchanged
    2   6.5478   0.0000   0.0000 180.0000   1.0000   0.0000   0.0000 <--water2 unchanged
    3  -4.8676   0.0000   2.8103 180.0000  -0.7071   0.0000  -0.7071

Since the rotation stays the same for the first two, the Cartesian scan which uses REDO_DF_ON_ROTATION=.TRUE. should get all energies correct for these two geometries. But not for the third. Let's see.

Reference energies are calculated using single point calculations (no energy-scan) with Cartesian GTOs in the aux basis. Why not reference energies with spherical GTOs in the aux basis??? See section on Spherical and Cartesian auxiliary bases. 

Energies -----------Geom 1---------------    ------------Geom2-------------  --------------Geom3------------
         reference         JK-tzvpp          reference          JK-tzvpp      reference          JK-tzvpp
         JK     aTZ     Spherical Cartesian  JK     aTZ    Spherical Cart     JK     aTZ    Spherical  Cart
------------------------------------------------------------------------------------------------------------
REDO-DF                    F         T                      F         T                        F        T
------------------------------------------------------------------------------------------------------------
E1elst   -3.736 -3.610  -3.527    -3.736   -6.897  -6.574  -6.718    -6.897  8.037  8.250  8.087  8.037
E1exch    0.137  0.135   0.129     0.137    1.718   1.722   1.624     1.718  5.369  5.412  1.008  5.368
E2ind    -0.181 -0.203 -37.08     -0.181   -0.556  -0.664 -97.861    -0.556 -0.727 -0.856 -6494. -11.62 
E2exind   0.001  0.002   0.101     0.001    0.026   0.051   3.616     0.026  0.022  0.193  263.1  1.300
E2disp   -0.919 -0.985  -1.263    -0.919   -2.351  -2.490  -3.384    -2.351 -3.334 -3.551 -8.835 -62.95
E2exdisp  0.007  0.011   0.003     0.007    0.054   0.072   0.030     0.054  0.204  0.247  0.744  3.859
============================================================================================================
  • Cartesian/REDO-DF:
    • E1elst & E1exch are always correct.
    • Second-order energies wrong once molecule is rotated. This was expected. See argument above.
  • Spherical/NO-REDO-DF:
    • E1elst & E1exch are always correct.
    • Second-order energies are always wrong. Why?

Are the JK-tzvpp Aux bases any good?

  • <math>\pi</math>-systems?
  • H-bonded systems?

CamCASP and truncated MO space

When molecules are too long and/or basis sets are too diffuse, DALTON will often need to truncate the MO space to enable the SCF cycle to converge. So the effective number of MOs will be less than the size of the basis used. Since CamCASP assumes these two are equal (though it doesn't need to), this results in errors.

Where do changes need to be made to fix this?

  • Reading in MOs.
  • Constructing the density-matrix.
  • Rotating MOs.
  • DF
  • Transformation code.
  • mol%ndim = mol%main%size : both fields need to be modified.

User-defined MOs

These can be read in ASCII form and this could be ideal for testing the code as we do not need to mess around with the DALTON interface to get this working.

DALTON interface

This needs to be modified to handle MO-space truncation. I have copied it to src/interfaces/ in the CamCASP directory tree. We cannot distribute it, but I will modify this version a lot. It writes integrals we do not need and which just occupy disk space.

When reading the basis data in subroutine readnbas, the basis size and orbital size is printed. I'm not sure what the difference is, but it could be useful if one is the actual size of the MO space and the other is the size of the un-truncated basis.

Integral Rotation

Still not complete. Things to do:

  • Hessians.
  • subroutine set_forcerotate_flag in df_int_operations.F90


Spherical versus Cartesian auxiliary bases

We definitely have a significant loss in accuracy when spherical GTOs are used in the auxiliary basis. The additional flexibilty of the Cartesian basis sets seems very important. Robust integrals may help reduce or even remove the problem, but can a larger, more diffuse shperical aux basis help?

Some thoughts:

  • The Cartesian d-functions include an s-function. But looking at the basis exponents it is not immediately obvious why this is an advantage as the most diffuse s- and d-functions have similar exponents. Can it be that the extra s-function in the Cartesian d-function, because it is fof the form <math> r^2 e^{-\alpha r^2}</math> is effectively even more diffuse than a standard s-function which is of the form <math> e^{-\alpha r^2}</math>?

If so, including even more diffuse s-functions may help. It's worth a try.

Retrieved from "http://wikis.ch.cam.ac.uk/cuc3/wiki/index.php?title=CamCASP/Programming/0&oldid=3695"