Dear all
Could this be the symmetry-breaking problem that
was lurking in the code some years ago?
there were situations in which (x,y) .ne. z, e.g. in the
forces (very slightly and not worrying quantitatively) in
situations where they should be equal (atom or cubic
system in cubic box with well chosen position w.r.t.
grid points). I don't think it was ever solved (just not
re-tested carefully enough).
Emilio
------------------------------------------------------------------------
Subject:
Re: [SIESTA-L] About dS12/dR12 calculated by matel
From:
Andrei Postnikov <[EMAIL PROTECTED]>
Date:
Thu, 26 Jan 2006 02:16:44 +0100
To:
SIESTA-L@LISTSERV.UAM.ES
Dear Siesteros,
On 2 Jan 2006, Hongjun Xiang reported the following problem:
| xianghjun wrote:
| | >Dear all,
| >
| >I want to get the values of dS_{\alpha,\beta}/dR, where
dS_{\alpha,\beta} is
| >the overlap matrix between atomic orbital |\alpha> centered at atom R_i
| >and |\beta>
| >centered at atom R_j. When \alpha =\beta, and i=j, I found that the
| >third component
| >of grSij from matel is not zero,
| >
| >For example, using the following setting,
| >#############################################################
| >and a call like this:
| >call matel( 'S', is, js, ioa, joa, xij(1,jn),Sij, grSij )
| >
| >I get Sij and grSij:
| >io,jo= 1 1 Sij 0.999999978 grSij 0.000000000 0.000000000 -0.040022905
| >io,jo= 2 2 Sij 0.999999990 grSij 0.000000000 0.000000000 -0.125447740
| >io,jo= 3 3 Sij 0.999999990 grSij 0.000000000 0.000000000 -0.125447740
| >io,jo= 4 4 Sij 0.999999990 grSij 0.000000000 0.000000000 -0.125447740
| >io,jo= 5 5 Sij 1.000000011 grSij 0.000000000 0.000000000 -0.023958303
| >io,jo= 6 6 Sij 1.000000011 grSij 0.000000000 0.000000000 -0.023958303
| >
| >In my opinion, grSij(1:3) should be zero when io == jo, so
| >I'm a little confused by the value of grSij calculated by matel.
| (...)
| I also test the results using the finite difference method and find that
| the values of grSij
| when io==jo from the finite difference method are different from those
| from matel. It seems
| that there might be a bug in matel.
| | Best regards,
| Hongjun Xiang
|
I tried to look into the issue, and here are my related notes.
First of all, when the problem is solved, the on-site gradient of S
is indeed zero, and the numerical derivative on moving the two coinciding
centers apart agrees with the analytical one, here are
my tests:
#
# Test of calculated gradients in matel
# System: two Al atoms (s + 3*p basis) in a box,
# separated along (001) by z, i.e.
# %block Atomic_Coordinates_and_Atomic_Species
# 0.000 0.000 -0.100 1
# 0.000 0.000 0.000 1
# %endblock Atomic_Coordinates_and_Atomic_Species
# Columns: z, property, analytical differentiation from matel;
# numerical differentiation as [F(i+1)-F(i-1)]/0.2
#
# ia,io,ioa,ja,jo,joa: 1 1 1 2 5 1 (S-S overlap)
#
#z(Bohr) S d(S)/dz num. 0.0000 1.00000 0.00000
0.1000 0.99892 -0.01303 -0.01415 0.2000 0.99717 -0.02297 -0.02400
0.3000 0.99412 -0.03849 -0.03820 0.4000 0.98953 -0.05304 -0.05270
0.5000 0.98358 -0.06546 #
#z(Bohr) T d(T)/dz num. 0.0000 0.40039 0.00000
0.1000 0.39874 -0.01874 -0.01985 0.2000 0.39642 -0.02872 -0.02975
0.3000 0.39279 -0.04419 -0.04380 0.4000 0.38766 -0.05780 -0.05725
0.5000 0.38134 -0.06825 #
#z(Bohr) Z d(Z)/dz num.
0.0000 0.00000 -0.50000
0.1000 -0.04995 -0.49881 -0.49860
0.2000 -0.09972 -0.49629 -0.49585
0.3000 -0.14912 -0.49129 -0.49095
0.4000 -0.19791 -0.48416 -0.48390
0.5000 -0.24590 -0.47543
#
#
# ia,io,ioa,ja,jo,joa: 1 1 1 2 6 2 (S-Py overlap)
z(Bohr) Y d(Y)/dz num.
0.0000 -1.43087 0.00000 0.1000 -1.42975 0.01756 0.01815
0.2000 -1.42724 0.03338 0.03395 0.3000 -1.42296 0.05234 0.05220
0.4000 -1.41680 0.07071 0.07050 0.5000 -1.40886 0.08779
Further come the notes to the problem, and how to fix it.
I suggest to fix this quasi-bug in matel, and I think it can be done
relatively easy, but I'd suggest the authors of matel choose
the most elegant way (by somehow employing the symmeries of the sherical
harmonics, as I suggest).
Problem:
--------
erroneous values of some gradients as delivered by matel.
I.e., the analytical gradient of the overlap of each basis function
with itself is not zero, contrary to expectation. An example for Al
s-function:
ia,io,ioa,ja,jo,joa: 1 1 1 1 1 1 xij= 0.0000 0.0000 0.0000
S = 1.00000E+00 grad S = 0.00000E+00 0.00000E+00 -9.71745E-03
T = 4.00390E-01 grad T = 0.00000E+00 0.00000E+00 -1.54155E-02
The values in the last column must be zero.
Explanation:
------------
This behavior can be traced to the construction in subr. matel
which avoids R12=0 (vector separating the centers of orbitals),
as follows:
...
X12(1) = R12(1)
X12(2) = R12(2)
X12(3) = R12(3)
R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) )
FOUND = .FALSE.
IF (R .GT. RCUT(IS1,IO1)+RCUT(IS2,IO2) ) THEN
FOUND = .TRUE.
ELSEIF (R .LT. TINY) THEN
X12(3) = TINY
R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) )
ENDIF
i.e. R12=0 is arbitrarily shifted to (0, 0, TINY)
and used in the evaluation of sperical harmonics.
The marix element and its gradient are recovered at the very end of matel,
in the summation over (l,m):
...
CALL SPLINT( RMAX/NR, FFR(0,1,JFFR),
. FFR(0,2,JFFR), NR+1, R, SR, DSRDR )
...
S12 = S12 + SR * FFY(IFFY) * Y(JLM)
DO IX = 1,3
DSDR(IX) = DSDR(IX) +
. DSRDR * FFY(IFFY) * Y(JLM) * X12(IX) / R +
. SR * FFY(IFFY) * DYDR(IX,JLM)
ENDDO
...
Subr. SPLINT makes interpolation of the radial function S_{l,m}(R)/R**L,
where R can be zero without problem, because an extrapolation is properly
done; FFY() are Gaunt coefficients, and Y(JLM) - spherical harmonics times
R**L.
The accumulation of S12 is along the formulae (25-28) of the basic
Siesta paper, J.Phys.Cond.Matter 14, 2745 (2002).
The construction of the gradient is straightforward, as
grad(S12) = grad(SR) * Y(JLM) + SR * grad( Y(JLM) ).
As SR depend on the radius only,
grad(SR) = d/dr(SR) * [x,y,z]/r ==> DSRDR * X12(IX) / R;
grad( Y(JLM) ) is returned from
CALL RLYLM( LMAX, X12, Y, DYDR )
as analytical gradient of Y_lm * R**L.
The avoided R12=0 affects only the calculation of spherical harmonics
and their gradients. In place of (0,0,0), they are evaluated
at (0,0,TINY). This gives a negligible error in the spherical harmonic
value, but sometimes noticeble one for the gradient.
Solution:
---------
I propose the evaluation of spherical harmonics at two split points,
(0,0,TINY) and (0,0,-TINY), and taking the average over them,
for both spherical harmonic value and it gradient. This will reduce
the error to the 3d order in TINY. A "fast and dirty" patch is
...
if (.not.RTINY) then
S12 = S12 + SR * FFY(IFFY) * Y(JLM)
DO IX = 1,3
DSDR(IX) = DSDR(IX) +
. DSRDR * FFY(IFFY) * Y(JLM) * X12(IX) / R +
. SR * FFY(IFFY) * DYDR(IX,JLM)
ENDDO
else S12 = S12 + SR * FFY(IFFY) *
(Y(JLM)+Y2(JLM)) / 2.
DO IX = 1,3
DSDR(IX) = DSDR(IX) +
. DSRDR * FFY(IFFY) * 0.5d0 * .
(Y(JLM) * X12(IX) + Y2(JLM) * X21(IX)) / R +
. SR * FFY(IFFY) * .
(DYDR(IX,JLM) + DYDR2(IX,JLM)) / 2.
ENDDO
endif
...
and the switch RTINY is introduced earlier in combination with
the "mirror vector" X21:
...
C Find if orbitals are out of range and avoid R12=0
RTINY = .FALSE.
X12(1) = R12(1)
X12(2) = R12(2) X12(3) = R12(3)
R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) )
FOUND = .FALSE.
IF (R .GT. RCUT(IS1,IO1)+RCUT(IS2,IO2) ) THEN
FOUND = .TRUE.
ELSEIF (R .LT. TINY) THEN
RTINY = .TRUE.
X12(3) = TINY
X21(1) = R12(1)
X21(2) = R12(2)
X21(3) = -TINY
R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) )
ENDIF
...
Moreover one has to calculate spherical harmonics for the "mirror vector"
X21:
...
CALL RLYLM( LMAX, X12, Y, DYDR )
if (RTINY) then
CALL RLYLM( LMAX, X21, Y2, DYDR2 )
endif
...
and allocate space for Y2 and DYDR2:
CALL RE_ALLOC( Y, 1,NILM, MYNAME//'Y' )
CALL RE_ALLOC( Y2, 1,NILM, MYNAME//'Y2' )
CALL RE_ALLOC( DYDR, 1,3, 1,NILM, MYNAME//'DYDR' )
CALL RE_ALLOC( DYDR2,1,3, 1,NILM, MYNAME//'DYDR2')
corespondingly in CALL DE_ALLOC(...) and declaring all new variables.
This doubled allocations certainly an overkill, and a better idea could be
to use symmetry relations of the spherical harmonics (and their gradients)
with respect to R -> -R, i.e. to write a function which immediately
returns Y(JLM) and DYDR(1..3,JLM) for -X12.
Best regards,
Andrei Postnikov
+-- Dr. Andrei Postnikov ---- Tel. +33-387315873 ----- mobile +33-666784053
---+
| Paul Verlaine University - Institute de Physique Electronique et Chimie,
|
| Laboratoire de Physique des Milieux Denses, 1 Bd Arago, F-57078 Metz,
France |
+-- [EMAIL PROTECTED] ------------
http://www.home.uni-osnabrueck.de/apostnik/ --+
