CamCASP/Bugs/8 - Revision history

import>Am592: /* Friday 18th Nov */

2011-11-21T10:25:10Z

Friday 18th Nov

← Older revision		Revision as of 10:25, 21 November 2011
Line 151:		Line 151:
	* Geometry of adenine is not being updated. So Ovr never changes.		* 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.		* 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:
			<pre>
			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
			</pre>

			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:
			<pre>
			!
			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
			!
			</pre>
			This works for all cases.

import>Am592: /* Adenine-Thymine */

2011-11-18T14:58:27Z

Adenine-Thymine

← Older revision		Revision as of 14:58, 18 November 2011
Line 150:		Line 150:
	* Overlaps only		* Overlaps only
	* Geometry of adenine is not being updated. So Ovr never changes.		* 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.

import>Am592: /* Adenine-Thymine */

2011-11-18T12:18:03Z

Adenine-Thymine

← Older revision		Revision as of 12:18, 18 November 2011
Line 145:		Line 145:
	<pre>/home/am592/DistProp/systems/adenine-thymine/scan/OUT_AJM_5pts_ovr_sw1_redoDF</pre> and same dir name with 'NOredoDF'.		<pre>/home/am592/DistProp/systems/adenine-thymine/scan/OUT_AJM_5pts_ovr_sw1_redoDF</pre> 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.		* 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.

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

2011-07-19T18:19:01Z

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

← Older revision		Revision as of 18:19, 19 July 2011
Line 125:		Line 125:
	Reference redo-DF no-redo-DF(rotate)		Reference redo-DF no-redo-DF(rotate)
	Sadlej/aTZ ----------------------------- ---------------------		Sadlej/aTZ ----------------------------- ---------------------
	MC/MC ~~CartAux~~ ~~CartAux~~ ~~SphAux~~ ~~CartAux~~		MC/MC SphAux SphAux CartAux SphAux
	Sph/Cart sw=0 sw=1 sw=1 sw=0 sw=1		Sph/Cart sw=0 sw=1 sw=1 sw=0 sw=1
	===========================================================================		===========================================================================
Line 137:		Line 137:

	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.		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==		==Adenine-Thymine==

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

2011-07-19T18:17:09Z

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

← Older revision		Revision as of 18:17, 19 July 2011
Line 123:		Line 123:
	All reference energies calculated as single-points with DALTON-CamCASP using the Sadlej/aTZ MC/MC basis. All energies in cm-1.		All reference energies calculated as single-points with DALTON-CamCASP using the Sadlej/aTZ MC/MC basis. All energies in cm-1.
	<pre>		<pre>
	Reference redo-DF no-redo-DF(rotate)		Reference redo-DF no-redo-DF(rotate)
	--------------------- ---------------------		Sadlej/aTZ ----------------------------- ---------------------
	~~sw=0~~ ~~sw=1~~ ~~sw=0 sw=1~~		MC/MC CartAux CartAux SphAux CartAux
			Sph/Cart sw=0 sw=1 sw=1 sw=0 sw=1
	====================================================================		===========================================================================
⚫	~~Geom1~~ -50.15 12.98 -~~348~~.98 12.98 NA*
	~~Geom2~~ -~~826~~.79 ~~-870~~.64 -~~1024~~.76 -~~870~~.64 NA		Geom1 -50.15 12.98 -348.98 -45.21 12.98 NA*
	~~Geom3~~ -~~547~~.57 -~~576~~.60 -~~780~~.54 -~~576~~.60 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)		*: Doesn't work with integral switch=1 (exact overlap and nuclear integrals)
	as there is an inconsistency in the rotations and MOs.		as there is an inconsistency in the rotations and MOs.
	</pre>		</pre>

			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.

	==Adenine-Thymine==		==Adenine-Thymine==

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

2011-07-19T18:11:20Z

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

← Older revision		Revision as of 18:11, 19 July 2011
Line 127:		Line 127:
	sw=0 sw=1 sw=0 sw=1		sw=0 sw=1 sw=0 sw=1
	====================================================================		====================================================================
	Geom1 -50.15 12.98 -348.98 12.98		Geom1 -50.15 12.98 -348.98 12.98 NA*
	Geom2 -826.79 -870.64 -1024.76 -870.64		Geom2 -826.79 -870.64 -1024.76 -870.64 NA
	Geom3 -547.57 -576.60 -780.54 -576.60		Geom3 -547.57 -576.60 -780.54 -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.
	</pre>		</pre>

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

2011-07-19T18:05:35Z

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

← Older revision		Revision as of 18:05, 19 July 2011
Line 121:		Line 121:

	===Recalculating <math>E^{(1)}_{\rm elst}</math>===		===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 reference energies calculated as single-points with DALTON-CamCASP using the Sadlej/aTZ MC/MC basis. All energies in cm-1.
	<pre>		<pre>
	Reference redo-DF no-redo-DF(rotate)		Reference redo-DF no-redo-DF(rotate)
	--------------------- ---------------------		--------------------- ---------------------
	sw=0 sw=1 sw=0 sw=1		sw=0 sw=1 sw=0 sw=1
	====================================================================		====================================================================
			Geom1 -50.15 12.98 -348.98 12.98
	Geom1
			Geom2 -826.79 -870.64 -1024.76 -870.64
	Geom2
			Geom3 -547.57 -576.60 -780.54 -576.60
	Geom3
	====================================================================		====================================================================
	</pre>		</pre>

import>Am592: /* Basis!!! */

2011-07-19T17:59:15Z

Basis!!!

← Older revision		Revision as of 17:59, 19 July 2011
Line 119:		Line 119:

	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.		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.
			<pre>
			Reference redo-DF no-redo-DF(rotate)
			--------------------- ---------------------
			sw=0 sw=1 sw=0 sw=1
			====================================================================
			Geom1
			Geom2
			Geom3
			====================================================================
			</pre>

	==Adenine-Thymine==		==Adenine-Thymine==

import>Am592: /* Basis!!! */

2011-07-19T15:57:19Z

Basis!!!

← Older revision		Revision as of 15:57, 19 July 2011
Line 112:		Line 112:
	========================================================================		========================================================================
	</pre>		</pre>

			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.

	==Adenine-Thymine==		==Adenine-Thymine==

import>Am592: /* Basis!!! */

2011-07-19T15:48:31Z

Basis!!!

@@ Line 101: / Line 101: @@
 E^{1}_{exch}(S2) = 2686.885
 Total density-overlap = 323.022 cm-1
 </pre>

CamCASP/Bugs/8 - Revision history

import>Am592: /* Friday 18th Nov */

import>Am592: /* Adenine-Thymine */

import>Am592: /* Adenine-Thymine */

import>Am592: /* Recalculating E^{(1)}_{\rm elst} */

import>Am592: /* Recalculating E^{(1)}_{\rm elst} */

import>Am592: /* Recalculating E^{(1)}_{\rm elst} */

import>Am592: /* Recalculating E^{(1)}_{\rm elst} */

import>Am592: /* Basis!!! */

import>Am592: /* Basis!!! */

import>Am592: /* Basis!!! */

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */

import>Am592: /* Recalculating $E^{(1)}_{\rm elst}$ */