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
fileenergy_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.
- Found some changes made to df_integral_utilities.F90 on mmp12 that were not committed. I seem to have found a fix to the problem but hadn't committed the changes. Testing these now.
Monday 21st Nov
Found the source of the problem on Saturday. Anthony independently found a work-around. He had defined Adenine (= MonA) and Thymine (=MonB). But contrary to the usual convention, was scanning Thymine with Adenine, i.e., MonB with MonA. This should have been allowed, but it failed.
The problem lay in part with subroutine compare_geometries from module molecule_operations and in part with a choice of default flags I had set in subroutine check_if_integrals_are_done.
Here's the BUGGY version of subroutine compare_geometries:
subroutine compare_geometries(A,B,same,translated,info)
!Compare the molecule geometries specified by two coordinate matrices.
!Determine if they are the same. If not, check if they are related by a
!translation. If not related by a translation, we will assume they are
!related by a rotation.
!
!info = 0 Normal exit
! = -1 Problems with the definition of A and/or B.
! = -2 The geometries supplied have different numbers of sites.
!
use parameters, only : epsilon_significant
implicit none
real(dp), dimension(:,:), intent(in) :: A, B
logical, intent(out) :: same, translated
integer, intent(out) :: info
!------------
character(len=*), parameter :: this_routine='compare_geometries'
integer :: i, nsites
real(dp) :: R, Rref
real(dp), external :: distance
!------------
info = 0
!
same = .true.
translated = .true.
!
...
...
!
R = 0.0_dp
Rref = -1.0_dp
nsites = size(A,2)
!
do i = 1, nsites
R = distance(A(:,i),B(:,i))
if (Rref.lt.0.0_dp) then
!Initialize Rref and if Rref <> 0, assume we have a translation
Rref = R
if (Rref.gt.epsilon_significant) translated = .true.
endif
if (R.gt.epsilon_significant) then
!These are not the same geometries...
same = .false.
!...but they could be related by a translation...
if (abs(Rref-R).gt.epsilon_significant) then
!...no.
translated = .false.
endif
endif
enddo
!
return
end subroutine compare_geometries
A rather simple routine. If the coordinated in A are different from those in B set same=F and if they are not related by a translation set translation=F. The problem lay with the latter. The routine decides that B is a translation of A if all coordinate differ by the same amount Rref which is set to be the first nonzero difference. So far so good. But the test that sets translation=F is performed only if dist(A,B)>0. This is wrong and fails for the case encountered in Anthony's example when testing to see if the dimer (constructed from MonA and MonB in that order) has altered. For here, we have, say, Rref=2.0 (since MonA is the one translated and the coordinates of MonA appear first in the dimer. But when we come to the coordinates of MonB, R=0.0 as MonB is at its original location. And the test for translations is not performed. So we get the flags: same=F,translated=T. This is wrong. It needs to set both to false!
Here's the corrected fragment of the code:
!
same = .true.
translated = .true.
!
...
...
R = 0.0_dp
Rref = -1.0_dp
nsites = size(A,2)
!
do i = 1, nsites
R = distance(A(:,i),B(:,i))
if (Rref<0.0_dp) then
!Initialize Rref and if Rref <> 0, assume we have a translation
Rref = R
if (Rref>epsilon_significant) translated = .true.
endif
if (R>epsilon_significant) then
!These are not the same geometries...
same = .false.
endif
!...but they could be related by a translation...
if (abs(Rref-R)>epsilon_significant) then
!...no.
translated = .false.
endif
enddo
!
This works for all cases.