Jose, Emilio: with respect to the matel issue discussed yesterday, I suggest now a more elegant substitution of my yesterday's fast and dirty fix. It explicitly uses the symmetry, that is that the spherical harmonics change with R -> -R as (-1)**l whereas gradients as (-1)**(l-1). The changes in matel amount to two variables introduced and the final block re-written; they are all marked as C --- AVP begin . LVAL, C --- AVP end etc. in the attached matel.f I could reproduce my numerical tests of yesterday (and even better, with some noice of the order of ~E-29 now removed). However, my test system included only (s,p) functions so far. Please advise. I'll put this new fix in the mailing list as well, for the case someone would like to do more tests.
Best regards, Andrei +-- 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/ --+
! ! Copyright (c) Fundacion General Universidad Autonoma de Madrid: ! E.Artacho, J.Gale, A.Garcia, J.Junquera, P.Ordejon, D.Sanchez-Portal ! and J.M.Soler, 1996-2003 ! ! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS ! "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR ! A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, ! SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, ! DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY ! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ! (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE ! OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! SUBROUTINE MATEL( OPERAT, IS1, IS2, IO1, IO2, R12, S12, DSDR ) C ******************************************************************* C Finds the overlap or laplacian matrix elements between two atomic C orbitals. Written by J.M.Soler. April 1995. C Matrix elements of the position operator by DSP. June, 1999 C Last modifications by JMS. January 2001. C ************************* INPUT *********************************** C CHARACTER OPERAT : Operator to be used: 'S' => Unity (overlap) C 'T' => -Laplacian C 'X' => x C 'Y' => y C 'Z' => z C What is actually computed is the matrix elements: C 'X' => < phi1(r-R12)|x|phi2(r)> C 'Y' => < phi1(r-R12)|y|phi2(r)> C 'Z' => < phi1(r-R12)|z|phi2(r)> C i.e the origin is taken at the position of the second atom. C INTEGER IS1 : Species index of 1st orbital atom C INTEGER IS2 : Species index of 2nd orbital atom C Allowed range of IS1, IS2: (1, MAXS), where C MAXS and MAXL (below) are internal parameters. C INTEGER IO1 : Orbital index of 1st orbital (within atom) C INTEGER IO2 : Orbital index of 2nd orbital (within atom) C Allowed range of IO1, IO2: (-MAXL**2, MAXL**2) C Indexes IS1,IS2,IO1,IO2 are used to call routines C RCUT and PHIATM (see below). C REAL*8 R12(3) : Vector from first to second atom C ************************* OUTPUT ********************************** C REAL*8 S12 : Matrix element between orbitals. C REAL*8 DSDR(3) : Derivative (gradient) of S12 with respect to R12. C ************************* ROUTINES CALLED ************************* C The following functions must exist: C C INTEGER FUNCTION LOFIO(IS,IO) C Returns the total angular momentum of orbitals and KB proyectors. C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C C REAL*8 FUNCTION RCUT(IS,IO) C Returns cutoff radius of orbitals and KB projectors. C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C C SUBROUTINE PHIATM(IS,IO,R,PHI,GRPHI) C Finds values and gradients of: C a) basis orbitals (IO > 0) C b) KB proyectors (IO < 0) C c) Local part of screened pseudopotential (IO = 0) ( b) and c) are C not required if MATEL is called only with IO > 0 ) C Input: C INTEGER IS : Species index C INTEGER IO : Orbital index C REAL*8 R(3) : Position with respect to atom C Output: C REAL*8 PHI : Value of orbital or KB projector at point R C REAL*8 GRPHI(3) : Gradient of PHI at point R C C SUBROUTINE XPHIATM(IS,IO,R,PHI,GRPHI) C The same as PHIATM but multiply by x. C C SUBROUTINE YPHIATM(IS,IO,R,PHI,GRPHI) C The same as PHIATM but multiply by y. C C SUBROUTINE ZPHIATM(IS,IO,R,PHI,GRPHI) C The same as PHIATM but multiply by z. C C ************************* UNITS *********************************** C Length units are arbitrary, but must be consistent in MATEL, RCUT C and PHIATM. The laplacian unit is (length unit)**(-2). C ************************* BEHAVIOUR ******************************* C 1) Returns exactly zero if |R12| > RCUT(IS1,IO1) + RCUT(IS2,IO2) C 2) If (IS1.LE.0 .OR. IS2.LE.0) all internal tables are erased for C reinitialization. C ******************************************************************* C 6 10 20 30 40 50 60 7072 C Modules ----------------------------------------------------------- USE ALLOC USE ATMFUNCS, ONLY: . LOFIO, PHIATM, RCUT, XPHIATM, YPHIATM, ZPHIATM use spher_harm use m_radfft C ------------------------------------------------------------------- C Argument types and dimensions ------------------------------------- IMPLICIT NONE CHARACTER OPERAT INTEGER IO1, IO2, IS1, IS2 DOUBLE PRECISION DSDR(3), R12(3), S12 C ------------------------------------------------------------------- C Internal precision parameters ------------------------------------ C NR is the number of radial points for matrix-element tables. C NQ is the number of radial points in reciprocal space. C EXPAND is a factor to expand some array sizes C Q2CUT is the required planewave cutoff for orbital expansion C (in Ry if lengths are in Bohr). C FFTOL is the tolerance for considering equal the radial part of C two orbitals. C TINY is a small number to add to a zero denominator INTEGER, PARAMETER :: NR = 128 INTEGER, PARAMETER :: NQ = 1024 DOUBLE PRECISION, PARAMETER :: EXPAND = 1.20D0 DOUBLE PRECISION, PARAMETER :: Q2CUT = 2.50D3 DOUBLE PRECISION, PARAMETER :: FFTOL = 1.D-8 DOUBLE PRECISION, PARAMETER :: TINY = 1.D-12 CHARACTER(LEN=*), PARAMETER :: MYNAME = 'MATEL ' C ------------------------------------------------------------------- C Internal variable types and dimensions ---------------------------- INTEGER :: . I, IF1, IF2, IFF, IFFY, IFLM1, IFLM2, . IO, IOPER, IQ, IR, IS, IX, . JF1, JF2, JFF, JFFR, JFFY, JFLM1, JFLM2, JLM, . JO, JO1, JO2, JR, . L, L1, L2, L3, LMAX, C --- AVP begin . LVAL, C --- AVP end . N, NILM, NILM1, NILM2, NJLM1, NJLM2 INTEGER, SAVE :: . MF=0, MFF=0, MFFR=0, MFFY=0, . NF=0, NFF=0, NFFR=0, NFFY=0 INTEGER, POINTER, SAVE :: . IFFR(:), ILM(:), ILMFF(:), INDF(:,:,:), INDFF(:,:,:), . INDFFR(:), INDFFY(:), NLM(:,:,:) DOUBLE PRECISION :: . C, CPROP, DFFR0, DFFRMX, DSRDR, FFL(0:NQ), FFQ(0:NQ), . Q, R, SR, X12(3) DOUBLE PRECISION, SAVE :: . PI, QMAX, RMAX DOUBLE PRECISION, POINTER, SAVE :: . CFFR(:), DYDR(:,:), F(:,:), FFR(:,:,:), FFY(:), Y(:) LOGICAL :: . FOUND, PROPOR C --- AVP begin . , RTINY C --- AVP end LOGICAL, SAVE :: . NULLIFIED=.FALSE. ! logical, save :: opened=.false. ! JMS debug TYPE(allocDefaults) :: . OLDEFS EXTERNAL . PROPOR, SPLINE, SPLINT, TIMER C ------------------------------------------------------------------- C Start time counter * CALL TIMER( MYNAME, 1 ) C Nullify pointers IF (.NOT.NULLIFIED) THEN NULLIFY( IFFR, ILM, ILMFF, INDF, INDFF, INDFFR, INDFFY, NLM, . CFFR, DYDR, F, FFR, FFY, Y ) ALLOCATE( INDF(0,0,0) ) NULLIFIED = .TRUE. ENDIF C Set allocation defaults CALL ALLOC_DEFAULT( OLD=OLDEFS, ROUTINE=MYNAME, . COPY=.TRUE., SHRINK=.FALSE. ) C Check if tables must be re-initialized IF ( IS1.LE.0 .OR. IS2.LE.0 ) THEN CALL DE_ALLOC( IFFR, MYNAME//'IFFR' ) CALL DE_ALLOC( ILM, MYNAME//'ILM' ) CALL DE_ALLOC( ILMFF, MYNAME//'ILMFF' ) CALL DE_ALLOC( INDF, MYNAME//'INDF' ) CALL DE_ALLOC( INDFF, MYNAME//'INDFF' ) CALL DE_ALLOC( INDFFR, MYNAME//'INDFFR' ) CALL DE_ALLOC( INDFFY, MYNAME//'INDFFY' ) CALL DE_ALLOC( NLM, MYNAME//'NLM' ) CALL DE_ALLOC( CFFR, MYNAME//'CFFR' ) CALL DE_ALLOC( DYDR, MYNAME//'DYDR' ) CALL DE_ALLOC( F, MYNAME//'F' ) CALL DE_ALLOC( FFR, MYNAME//'FFR' ) CALL DE_ALLOC( FFY, MYNAME//'FFY' ) CALL DE_ALLOC( Y, MYNAME//'Y' ) MF = 0 MFF = 0 MFFR = 0 MFFY = 0 NF = 0 NFF = 0 NFFR = 0 NFFY = 0 ENDIF C Check argument OPERAT IF ( OPERAT .EQ. 'S' ) THEN IOPER = 1 ELSEIF ( OPERAT .EQ. 'T' ) THEN IOPER = 2 ELSEIF ( OPERAT .EQ. 'X' ) THEN IOPER = 3 ELSEIF ( OPERAT .EQ. 'Y' ) THEN IOPER = 4 ELSEIF ( OPERAT .EQ. 'Z' ) THEN IOPER = 5 ELSE call die('MATEL: Invalid value of argument OPERAT') ENDIF C Check size of orbital index table IF ( MAX(IS1,IS2).GT.SIZE(INDF,1) .OR. . MIN(IO1,IO2).LT.LBOUND(INDF,2) .OR. . MAX(IO1,IO2).GT.UBOUND(INDF,2) .OR. . MAX(IOPER,2).GT.SIZE(INDF,3) ) THEN CALL RE_ALLOC( INDF, 1,MAX(SIZE(INDF,1),IS1,IS2), . MIN(LBOUND(INDF,2),IO1,IO2),MAX(UBOUND(INDF,2), . IO1,IO2),1,MAX(SIZE(INDF,3),IOPER,2), . MYNAME//'INDF' ) CALL RE_ALLOC( NLM, 1,MAX(SIZE(INDF,1),IS1,IS2), . MIN(LBOUND(INDF,2),IO1,IO2),MAX(UBOUND(INDF,2), . IO1,IO2),1,MAX(SIZE(INDF,3),IOPER,2), . MYNAME//'NLM' ) ENDIF C JMS debug c if (.not.opened) then c open(41,file='matel.out',form='formatted',status='unknown') c opened = .true. c endif C Find radial expansion of each orbital ----------------------------- DO I = 1,2 IF (I .EQ. 1) THEN IS = IS1 IO = IO1 FOUND = (INDF(IS,IO,1) .NE. 0) ELSE IS = IS2 IO = IO2 FOUND = (INDF(IS,IO,IOPER) .NE. 0) ENDIF IF (.NOT.FOUND) THEN * CALL TIMER( 'MATEL1', 1 ) PI = 4.D0 * ATAN(1.D0) C Factor 2 below is because we will expand the product of C two orbitals QMAX = 2.D0 * SQRT( Q2CUT ) RMAX = PI * NQ / QMAX IF ( RCUT(IS,IO) .GT. RMAX ) THEN call die('MATEL: NQ too small for required cutoff.') ENDIF C Reallocate arrays L = LOFIO( IS, IO ) NILM = (L+2)**2 IF (NF+NILM .GT. MF) MF = EXPAND * (NF+NILM) CALL RE_ALLOC( F, 0,NQ, 1,MF, MYNAME//'F' ) CALL RE_ALLOC( ILM, 1,MF, MYNAME//'ILM' ) CALL RE_ALLOC( INDFF, 1,MF, 1,MF, 1,MAX(IOPER,2), . MYNAME//'INDFF' ) C Expand orbital in spherical harmonics IF ((I.EQ.1) .OR. (IOPER.LE.2)) THEN CALL YLMEXP( L, RLYLM, PHIATM, IS, IO, 0, NQ, RMAX, . NILM, ILM(NF+1:), F(0:,NF+1:) ) INDF(IS,IO,1) = NF+1 INDF(IS,IO,2) = NF+1 NLM(IS,IO,1) = NILM NLM(IS,IO,2) = NILM ELSE IF(IOPER.EQ.3) THEN CALL YLMEXP( L+1, RLYLM, XPHIATM, IS, IO, 0, NQ, RMAX, . NILM, ILM(NF+1:), F(0:,NF+1:) ) ELSEIF(IOPER.EQ.4) THEN CALL YLMEXP( L+1, RLYLM, YPHIATM, IS, IO, 0, NQ, RMAX, . NILM, ILM(NF+1:), F(0:,NF+1:) ) ELSE CALL YLMEXP( L+1, RLYLM, ZPHIATM, IS, IO, 0, NQ, RMAX, . NILM, ILM(NF+1:), F(0:,NF+1:) ) ENDIF INDF(IS,IO,IOPER) = NF+1 NLM(IS,IO,IOPER) = NILM ENDIF C Store orbital in k-space DO JLM = 1,NILM NF = NF + 1 L = LOFILM( ILM(NF) ) CALL RADFFT( L, NQ, RMAX, F(0:NQ,NF), F(0:NQ,NF) ) * F(NQ,NF) = 0.D0 ENDDO C JMS debug c write(41,'(/,a,4i6)') 'is,io,l,nf =', is, io, l, nf c write(41,'(i6,e25.12)') (iq,f(iq,nf),iq=0,nq) * CALL TIMER( 'MATEL1', 2 ) ENDIF ENDDO C ------------------------------------------------------------------- C Find radial expansion of overlap ---------------------------------- IF1 = INDF(IS1,IO1,1) IF2 = INDF(IS2,IO2,IOPER) IF ( INDFF(IF1,IF2,IOPER) .EQ. 0 ) THEN * CALL TIMER( 'MATEL2', 1 ) NILM1 = NLM(IS1,IO1,1) NILM2 = NLM(IS2,IO2,IOPER) DO IFLM1 = IF1,IF1+NILM1-1 DO IFLM2 = IF2,IF2+NILM2-1 C Check interaction range IF (RCUT(IS1,IO1)+RCUT(IS2,IO2) .GT. RMAX) THEN call die('MATEL: NQ too small for required cutoff.') ENDIF C Find orbitals convolution by multiplication in k-space C = ( 2.D0 * PI )**1.5D0 DO IQ = 0,NQ FFQ(IQ) = C * F(IQ,IFLM1) * F(IQ,IFLM2) IF ( OPERAT .EQ. 'T' ) THEN Q = IQ * QMAX / NQ FFQ(IQ) = FFQ(IQ) * Q * Q ENDIF ENDDO C Loop on possible values of l quantum number of product L1 = LOFILM( ILM(IFLM1) ) L2 = LOFILM( ILM(IFLM2) ) CALL RE_ALLOC( IFFR, 0,L1+L2, MYNAME//'IFFR' ) CALL RE_ALLOC( CFFR, 0,L1+L2, MYNAME//'CFFR' ) DO L3 = ABS(L1-L2), L1+L2, 2 C Return to real space CALL RADFFT( L3, NQ, NQ*PI/RMAX, FFQ, FFL ) * FFL(NQ) = 0.D0 IF (MOD(ABS(L1-L2-L3)/2,2) .NE. 0) THEN DO IR = 0,NQ FFL(IR) = - FFL(IR) ENDDO ENDIF C Divide by R**L IF (L3 .NE. 0) THEN DO IR = 1,NQ R = IR * RMAX / NQ FFL(IR) = FFL(IR) / R**L3 ENDDO C Parabolic extrapolation to R=0 FFL(0) = ( 4.D0 * FFL(1) - FFL(2) ) / 3.D0 ENDIF C Select NR out of NQ points C Copy NQ to a variable, to fool the compiler N = NQ IF (MOD(N,NR) .NE. 0) . call die('MATEL: NQ must be multiple of NR') DO IR = 0,NR JR = IR * NQ / NR FFL(IR) = FFL(JR) ENDDO C Find if radial function is already in table FOUND = .FALSE. SEARCH: DO JO1 = LBOUND(INDF,2),UBOUND(INDF,2) DO JO2 = LBOUND(INDF,2),UBOUND(INDF,2) JF1 = INDF(IS1,JO1,1) JF2 = INDF(IS2,JO2,IOPER) IF (JF1.NE.0 .AND. JF2.NE.0) THEN NJLM1 = NLM(IS1,JO1,1) NJLM2 = NLM(IS2,JO2,IOPER) DO JFLM1 = JF1,JF1+NJLM1-1 DO JFLM2 = JF2,JF2+NJLM2-1 JFF = INDFF(JFLM1,JFLM2,IOPER) IF (JFF .NE. 0) THEN DO JFFY = INDFFY(JFF-1)+1, INDFFY(JFF) JFFR = INDFFR(JFFY) IF ( PROPOR(NR,FFL(1),FFR(1,1,JFFR), . FFTOL,CPROP) ) THEN FOUND = .TRUE. IFFR(L3) = JFFR CFFR(L3) = CPROP EXIT SEARCH ENDIF ENDDO ENDIF ENDDO ENDDO ENDIF ENDDO ENDDO SEARCH C Store new radial function and setup spline interpolation IF (.NOT.FOUND) THEN NFFR = NFFR + 1 IF (NFFR .GT. MFFR) THEN MFFR = EXPAND * NFFR CALL RE_ALLOC( FFR, 0,NR, 1,2, 1,MFFR, MYNAME//'FFR' ) ENDIF IFFR(L3) = NFFR CFFR(L3) = 1.D0 DO IR = 0,NR FFR(IR,1,NFFR) = FFL(IR) ENDDO DFFR0 = HUGE(1.D0) DFFRMX = 0.D0 CALL SPLINE( RMAX/NR, FFR(0,1,NFFR), NR+1, DFFR0, DFFRMX, . FFR(0,2,NFFR) ) C JMS debug c write(41,'(/,a,4i6)') 'if1,if2,nffr =', if1,if2,nffr c write(41,'(i6,2e25.12)') c . (ir,ffr(ir,1,nffr),ffr(ir,2,nffr),ir=0,nr) ENDIF ENDDO C Reallocate some arrays NILM = (L1+L2+1)**2 IF (NFFY+NILM .GT. MFFY) THEN MFFY = EXPAND * (NFFY+NILM) CALL RE_ALLOC( FFY, 1,MFFY, MYNAME//'FFY' ) CALL RE_ALLOC( ILMFF, 1,MFFY, MYNAME//'ILMFF' ) CALL RE_ALLOC( INDFFR, 1,MFFY, MYNAME//'INDFFR' ) ENDIF CALL RE_ALLOC( Y, 1,NILM, MYNAME//'Y' ) CALL RE_ALLOC( DYDR, 1,3, 1,NILM, MYNAME//'DYDR' ) C Expand the product of two spherical harmonics (SH) also in SH CALL YLMEXP( L1+L2, RLYLM, YLMYLM, ILM(IFLM1), ILM(IFLM2), . 1, 1, 1.D0, NILM, ILMFF(NFFY+1:), FFY(NFFY+1:) ) C Loop on possible lm values of orbital product DO I = 1,NILM NFFY = NFFY + 1 JLM = ILMFF(NFFY) L3 = LOFILM( JLM ) INDFFR(NFFY) = IFFR(L3) FFY(NFFY) = FFY(NFFY) * CFFR(L3) ENDDO NFF = NFF + 1 IF (NFF .GT. MFF) THEN MFF = EXPAND * NFF CALL RE_ALLOC( INDFFY, 0,MFF, MYNAME//'INDFFY' ) ENDIF INDFFY(NFF) = NFFY INDFF(IFLM1,IFLM2,IOPER) = NFF ENDDO ENDDO * CALL TIMER( 'MATEL2', 2 ) ENDIF C End of initialization section ------------------------------------- C Find value of matrix element and its gradient --------------------- * CALL TIMER( 'MATEL3', 1 ) C Initialize output S12 = 0.D0 DSDR(1) = 0.D0 DSDR(2) = 0.D0 DSDR(3) = 0.D0 C Find if orbitals are out of range and avoid R12=0 C --- AVP begin RTINY = .FALSE. C --- AVP end 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 C --- AVP begin RTINY = .TRUE. C --- AVP end X12(3) = TINY R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) ) ENDIF C Find spherical harmonics times R**L IF (.NOT.FOUND) THEN IF1 = INDF(IS1,IO1,1) IF2 = INDF(IS2,IO2,IOPER) NILM1 = NLM(IS1,IO1,1) NILM2 = NLM(IS2,IO2,IOPER) DO IFLM1 = IF1,IF1+NILM1-1 DO IFLM2 = IF2,IF2+NILM2-1 IFF = INDFF(IFLM1,IFLM2,IOPER) LMAX = 0 DO IFFY = INDFFY(IFF-1)+1, INDFFY(IFF) JLM = ILMFF(IFFY) LMAX = MAX( LMAX, LOFILM(JLM) ) ENDDO CALL RLYLM( LMAX, X12, Y, DYDR ) C Interpolate radial functions and obtain SH expansion DO IFFY = INDFFY(IFF-1)+1, INDFFY(IFF) JFFR = INDFFR(IFFY) CALL SPLINT( RMAX/NR, FFR(0,1,JFFR), . FFR(0,2,JFFR), NR+1, R, SR, DSRDR ) JLM = ILMFF(IFFY) C --- AVP begin if (.not.RTINY) then C straightforward calculation of function and gradient 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 C special treatment near R=0: take an average of values C at (0,0,TINY) and (0,0,-TINY) for both function and its gradient. C The averaging depends on the symmetry properties C of spherical harmonics and hence on the L value. C Extract L-value from the global index: LVAL = SQRT(real(JLM-1)) if (mod(LVAL,2).eq.0) then C Y(JLM) modifies with R -> -R as (-1)^LVAL, therefore C when averaged over (R,-R) contriutes only for LVAL even, C otherwise stricly zero: S12 = S12 + SR * FFY(IFFY) * Y(JLM) else C In the gradient, both Y(JLM)*X12(IX) and DYDR change C with R -> -R as (-1)^[LVAL+1], therefore averaging C over (R, -R) yields non-zero only for LVAL odd: DO IX = 1,3 DSDR(IX) = DSDR(IX) + . DSRDR * FFY(IFFY) * Y(JLM) * X12(IX) / R + . SR * FFY(IFFY) * DYDR(IX,JLM) ENDDO endif endif ! end of special treatment R=0 C --- AVP end ENDDO ENDDO ENDDO ENDIF * CALL TIMER( 'MATEL3', 2 ) C ------------------------------------------------------------------- C Restore allocation defaults CALL ALLOC_DEFAULT( RESTORE=OLDEFS ) C Stop time counter * CALL TIMER( MYNAME, 2 ) END SUBROUTINE MATEL