Difference between revisions of "CamCASP/Programming/0"

Latest revision as of 15:20, 8 December 2011

CamCASP => Programming => Status

Propagator

There have been major changes to the propagator modules. See Propagator for details. --alston 18:42, 5 May 2010 (BST)

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

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:

make_prop
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.
Solution of LR-DFT equations.
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?

The DF-FDDS is a purely monomer property: <math>C_{kl} = \langle k | \hat{C} | l \rangle</math>. On rotation by rotation operator <math>\hat{R}</math>, it transforms as

<math> \langle \hat{R} k | \hat{C} | \hat{R} l \rangle = {R^{A'A}_{k'k}}^{-1} \langle k' | \hat{C} | l' \rangle {R^{A'A}_{l'l}}^{-1} </math>

or

<math> R^{AA'}_{kk'} \langle k' | \hat{C} | l' \rangle R^{AA'}_{ll'} </math>

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 is always correct (apart from issues because of the basis). Why is E1exch wrong for Geom3?????
- Second-order energies are always wrong. Why?

I've found an even simpler instance of a bug

Water dimer. Geom3. Single point calculation without Energy-Scan. The REDO-DF-ON-ROTATION flag should have no effect on single point calculations but it does. If set to true all is well. But if false, integrals are rotated even though no rotations are needed as MOs and molecular geometries are consistent. Here are the commands that lead to the error:

SET DF
  Eta = 0.0
  Lambda = 0.0
  REDO-DF-ON-ROTATION  False
END

Edit water2
  ! This is geom 3:
  Units Bohr Degrees
  Rotate by 180.0000 Degrees about -0.7071   0.0000  -0.7071
  Place at XYZ -4.8676   0.0000   2.8103
  Update Geom & MOs
End

Begin DF
  Dimer
  Type OV
End
BEGIN DF
  Molecule water1
  Type FULL
  Print Nothing
END
BEGIN DF
  Molecule water2
  Type FULL
  Print Nothing
END

Begin df-ints
  Int     Type NUCA Desc __XB Print
  Int     Type NUCA Desc OOBB Print
End

Begin E1elst
End

Finish

The NUCA,__XB integrals will be incorrect if REDO-DF-ON-ROTATION is set False. This is because the EDIT water2 command rotates water2 and updates both the geometry and molecular orbitals. IN fact, even if no update commands were used here, both geometry and orbitals would be automatically updated in the monomer DF routines. So there would be no need for the integrals to be rotated (rotations are done in df_integrals.F90).

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.

The problem probably arises because the basis includes p functions such as <math> x e^{-\alpha r^2}</math>, so the density includes substantial contributions from <math> x^2 e^{-2\alpha r^2}</math> etc. and needs <math> (x^2+y^2+z^2) e^{-2\alpha r^2}</math> to describe these properly. Just adding more diffuse s functions won't do it.

It occurred to me that we could pinpoint the reason for the loss in accuracy by using the water dimer in a configuration where the penetration part of the electrostatic energy is negligible, but the electrostatic energy is still large because of the permanent moments of water. In this case, both the spherical and Cartesian bases should result in the same energy particularly if we impose the charge constraint.

Test: water dimer

We need a good reference, but this is now quite hard because the old scripts that could run SAPT are no longer supported. These could and probably should be re-activated.

But we do have a water scan at a single orientation and this should be used for the present. In addition, I will do a scan with about 50 random geometries.

Some questions that need answering:

Is MC+ always more accurate than MC?
- Sadlej/MC versus Sadlej/MC+
What aux basis do we use for the mid-bonds in the MC+ type?
- Do variations in this basis have a significant impact on energies?
Which functions should we add to the MC bases to improve accuracy?
How sensitive is the electrostatic self-energy to the auxiliary basis? This quantity is a monomer quantity and should be less sensitive.