CamCASP/Bugs/8

CamCASP => Bugs => Energy Scan

Error in Energy Scan

Spotted by Anthony with the adenine-thymine system. Data in

/home/am592/DistProp/systems/adenine-thymine/scan

Anthony had calculated all energy components *before* calling ENERGY-SCAN. One geometry in ENERGY-SCAN was the original geometry, and for this case, the energy components seem correct. That is, they agree with those calculated before ENERGY-SCAN was called.

What this seems to suggest is that not all integrals are being correctly re-calculated when geometries are changed.

Now, tests done on He2 and (H2O)2 do not show these problems. These calculations are in

/home/am592/CamCASP/tests/energy_scan/

But these are small systems. So perhaps this is a large-system problem. I've tried reducing the memory given to the (H2O)2 job to try to force this error, but it is still correct. But this test is not conclusive as the way objects are released from memory is hard to control.

Pyridine also shows this problem. Here, <math>E^{(1)}_{\rm exch}</math> comes out correct, but not <math>E^{(1)}_{\rm elst}</math>. Data in

/home/am592/DistProp/systems/pyridine/scans/overlap/JK-tzvpp/scan3-3000pts

For the larger adenine-thymine system, not even <math>E^{(1)}_{\rm exch}</math> is correct.

The scan fails even when the monomers are only translated.

Pyridine tests in

/home/am592/CamCASP/tests/energy_scan/pyridine2/

• Translations seem correct. Tested using single point calculations within CamCASP (i.e., translation defined in MOLECULE block)
• Rotations definitely wrong. Compared with earlier calculations in MC+ basis. Elst should be very similar. it is is vastly different. See

   /home/am592/CamCASP/tests/energy_scan/pyridine2/OUT_3pts_old_scan

  Compare with data in

  /home/am592/DistProp/systems/pyridine/overlap_models/old-fits

  file

  energy_sadlejMC+_scan2_L4_C6iso_scale.dat

• Did single point calculation. Rotated one molecule using CamCASP MOLECULE specification. Results differ from CamCASP calculation with ENERGY-SCAN. Both E1elst and E1exch different. See outputs in

  /home/am592/CamCASP/tests/energy_scan/pyridine2/single_pts/OUT_geom1_old_set

  . This means that the error lies in the way rotations are handled. Why don't we see this with smaller molecules???
• Another single point calculation: this time translations only. Agree with ENERGY-SCAN module.

Need a clear example without ambiguities:

• Complete MCBS calculation with MOs obtained with DALTON (i.e. not using rotations in CamCASP).
• Then compare with rotations done in CamCASP.
• Print out integrals for E1elst. There are only three of these.
• Use pyridine.
• Use geometry: MOLECULE pyridine2 at -5.870040E-01 3.957260E+00 9.016280E+00 rotated by 1.227930E+02 about -9.097780E-01 -1.008600E-01 4.026560E-01
• BUG in CamCASP: Global units need to be used in this line in the MOLECULE block. At present, no conversion is made, i.e. assuming BOHR and DEGREES.

Pyridine dimer

• Reference calculations: DALTON+CamCASP: single point calculations at three geometries. results in

  /home/am592/CamCASP/tests/energy_scan/pyridine2/ref

• Comparisons with same geometries but calculations done in CamCASP with ENERGY-SCAN module: <math>E^{(1)}_{\rm elst}</math> definitely rubbish. Small errors (3rd/4th place) in <math>E^{(1)}_{\rm exch}</math>. I'm not convinced it is completely correct. But the errors could be due to inevitable differences in the MOs obtained by rotation (in CamCASP) and by calculation using DALTON.
• The original reference calculations were performed using the MC main basis and DC aux basis (MC/DC type). This is not right as the energy-scan used MC/MC. Redoing the CamCASP calculations in the MC/MC type.

Basis!!!

<math>E^{(1)}_{\rm elst}</math> is very sensitive to the auxiliary basis used (something I've come across before but have forgotten! Here's some data. All calculations in tati:/home/am592/CamCASP/tests/energy_scan/pyridine2/overlap-scans/

Geometry 2 (Bohr Degree):

-5.199820E-01  3.505430E+00  7.986820E+00  1.227930E+02 -9.097780E-01 -1.008600E-01  4.026560E-01

Reference (sadlej/aTZ MC/MC single point calculation):

 E^{1}_{elst} =      -826.795351 CM-1
 Electron--electron energy =  43394789.727237 CM-1
 Electron--nuclear (A) energy = -43272011.950414 CM-1
 Electron--nuclear (B) energy = -43327497.537571 CM-1
 Nuclear --Nuclear  energy =  43203892.965397 CM-1

   E1 =       845.426620 CM-1
   E2a=     -1119.375223 CM-1
   E2b=     -1165.959718 CM-1
   E3a=       -23.684813 CM-1
   E3b=        -9.622773 CM-1
   E4 =       129.948984 CM-1
 E^{1}_{exch}(S2) =      2686.533842 CM-1

Total density-overlap S[rhoA,rhoB] =   322.944665 CM-1

ENERGY-SCAN: JK-tzvpp aux basis:

 E^{1}_{elst} =     -3045.759516 CM-1
 Electron--electron energy =  43392594.787698 CM-1
 Electron--nuclear (A) energy = -43272011.960473 CM-1
 Electron--nuclear (B) energy = -43327521.561670 CM-1
 Nuclear --Nuclear  energy =  43203892.974930 CM-1

   E1 =       845.035154 CM-1
   E2a=     -1119.653037 CM-1
   E2b=     -1166.014566 CM-1
   E3a=       -20.753604 CM-1
   E3b=        -8.242699 CM-1
   E4 =       129.813805 CM-1
E^{1}_{exch}(S2) = 2679.63 cm-1
Total density-overlap = 292.2 cm-1

ENERGY-SCAN: aTZ aux basis:

 E^{1}_{elst} =      -819.355723 CM-1
 Electron--electron energy =  43394821.191491 CM-1
 Electron--nuclear (A) energy = -43272011.960473 CM-1
 Electron--nuclear (B) energy = -43327521.561670 CM-1
 Nuclear --Nuclear  energy =  43203892.974930 CM-1

   E1 =       845.505864 CM-1
   E2a=     -1119.581265 CM-1
   E2b=     -1166.007448 CM-1
   E3a=       -23.685578 CM-1
   E3b=        -9.638240 CM-1
   E4 =       129.964225 CM-1
E^{1}_{exch}(S2) = 2686.885
Total density-overlap = 323.022 cm-1

Summary

                 Ref(sadlej/aTZ)   Eng-Scan(sad/JK)   End-Scan(sad/aTZ)
========================================================================
E^{1}_{elst}      -826.79          -3045.75  (-268%)   -819.35 (0.89%)
E^{1}_{exch}(S2)  2686.53           2679.63  (0.25%)   2686.88 (0.01%)
Sovr               322.94            292.2   (9.5%)     323.02 (-0.02%)
========================================================================

We know about the large errors (9.5%) in Sovr with the JK-tzvpp basis, but there is little we can do about this at present. The exchange-repulsion is not bad in the JK basis, but is essentially exact in the aTZ aux basis. It should be exact, but errors creep in from small differences in the molecular geometries obtained using the two routes (single-point calculation with DALTON+CamCASP versus ENERGY-SCAN and molecule rotations). Though perhaps it is convergence errors in DALTON that are the main culprit here. This should be looked into.

What's very disturbing is the large error in the electrostatic energy calculated with the JK-tzvpp auxiliary basis. It is awful. The error arises from the e-e repulsion term (see energy break-up shown above). Is this typical?

Well, if this is the only problem, there is a simple fix: the electrostatic energy can be calculated using the larger aTZ aux basis. Needs checking. In principle, only this component need be calculated in this way. The exchange-repulsion is already good enough.

Recalculating <math>E^{(1)}_{\rm elst}</math>

All reference energies calculated as single-points with DALTON-CamCASP using the Sadlej/aTZ MC/MC basis. All energies in cm-1.

          Reference              redo-DF                no-redo-DF(rotate)
      Sadlej/aTZ     -----------------------------   ---------------------
          MC/MC        SphAux      SphAux   CartAux    SphAux    
         Sph/Cart       sw=0        sw=1      sw=1      sw=0       sw=1
===========================================================================
Geom1    -50.15       12.98        -348.98    -45.21     12.98       NA*
Geom2   -826.79     -870.64       -1024.76   -819.35   -870.64       NA
Geom3   -547.57     -576.60        -780.54   -551.18   -576.60       NA
===========================================================================
*: Doesn't work with integral switch=1 (exact overlap and nuclear integrals)
 as there is an inconsistency in the rotations and MOs. 

Clearly, the best approach is to use the same basis (main & aux) as the reference calculation. In this case, use a spherical main basis and a Cartesian auxiliary basis. Agreement is very good for all geometries. Residual errors have been (briefly) discussed above.

So, method of choice: REDO DF, Cartesian auxiliary basis, integral switch = 1 (exact 1-e integrals)

Adenine-Thymine

• ENERGY-SCAN using translations only.
• Calculated only density overlap. With and without redoing DF. Data in

/home/am592/DistProp/systems/adenine-thymine/scan/OUT_AJM_5pts_ovr_sw1_redoDF

and same dir name with 'NOredoDF'.

• No change in Total density overlap with increasing separation. Therefore OVRL_XaXb not updated although it should have been.

Friday 18th Nov

Back to this problem. Running on mmp12.

• Overlaps only
• Geometry of adenine is not being updated. So Ovr never changes.