1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8! 9 module m_new_matel 10 public :: new_matel 11 CONTAINS 12 SUBROUTINE new_MATEL( OPERAT, IG1, IG2, R12, S12, DSDR ) 13C ******************************************************************* 14C Finds two-center matrix elements between 'atomic orbitals' 15C with finite radial and angular momentum cutoffs. 16C Written by J.M.Soler. April 1995. 17C Matrix elements of the position operator by DSP. June, 1999 18C Electrostatic interaction added by JMS. July 2002. 19C Introduction of unified 'global' indexes by AG. July 2011. 20C ************************* INPUT *********************************** 21C CHARACTER OPERAT : Operator to be used. The valid options are: 22C 'S' => Unity (overlap). Uppercase required for all values. 23C 'T' => -Laplacian 24C 'U' => 1/|r'-r| (with evaluate returning charge distributions) 25C 'X' => x, returning <phi1(r-R12)|x|phi2(r)> (origin on second atom) 26C 'Y' => y, returning <phi1(r-R12)|y|phi2(r)> 27C 'Z' => z, returning <phi1(r-R12)|z|phi2(r)> 28C INTEGER IG1 : Global index of 1st function (must be positive) 29C INTEGER IG2 : Gloabal index of 2nd function 30C Indexes IG1, IG2 are used only to call 31C routines LCUT, RCUT and EVALUATE (see below), and 32C may have other meanings within those routines 33C REAL*8 R12(3) : Vector from first to second atom 34C ************************* OUTPUT ********************************** 35C REAL*8 S12 : Matrix element between orbitals. 36C REAL*8 DSDR(3) : Derivative (gradient) of S12 with respect to R12. 37C ************************* ROUTINES CALLED ************************* 38C The following functions must exist: 39C 40C INTEGER FUNCTION LCUT(IG) 41C Returns the maximum angular momentum of the functions 42C Input: 43C INTEGER IG : Global function index 44C 45C REAL*8 FUNCTION RCUT(IG) 46C Returns cutoff radius of the function 47C Input: 48C INTEGER IG : Global function index 49C 50C SUBROUTINE EVALUATE(IG,R,PHI,GRPHI) 51C Returns the value of the functions to be integrated 52C Input: 53C INTEGER IG : Global function index 54C REAL*8 R(3) : Position with respect to atom 55C Output: 56C REAL*8 PHI : Value of orbital at point R 57C REAL*8 GRPHI(3) : Gradient of PHI at point R 58C ************************* UNITS *********************************** 59C Length units are arbitrary, but must be consistent in MATEL, RCUT 60C and EVALUATE. The laplacian unit is (length unit)**(-2). 61C ************************* BEHAVIOUR ******************************* 62C 1) If |R12| > RCUT(IS1,IO1) + RCUT(IS2,IO2), returns U(Rmax)*Rmax/R 63C for OPERAT='U', and exactly zero in all other cases. 64C 2) If (IS1.LE.0 .OR. IS2.LE.0) all internal tables are erased for 65C reinitialization and nothing is calculated. 66C ******************************************************************* 67C 6 10 20 30 40 50 60 7072 68 69C Modules ----------------------------------------------------------- 70 use precision, only : dp 71 use parallel, only: node 72 USE ALLOC 73 USE m_matel_registry, ONLY: EVALUATE, LCUT, RCUT 74 USE m_matel_registry, ONLY: EVALUATE_X, EVALUATE_Y, EVALUATE_Z 75 use interpolation, only: spline, splint 76 use m_errorf, only: derf 77 use spher_harm, only: rlylm, ylmexp, ylmylm, lofilm 78 use spher_harm, only: reset_spher_harm 79 use sys, only: die 80 use m_radfft 81 use siesta_options, only: matel_NRTAB 82C ------------------------------------------------------------------- 83 84C Argument types and dimensions ------------------------------------- 85 IMPLICIT NONE 86 CHARACTER OPERAT 87 INTEGER IG1, IG2 88 real(dp) DSDR(3), R12(3), S12 89C ------------------------------------------------------------------- 90 91C Internal precision parameters ------------------------------------ 92C NRTAB is the number of radial points for matrix-element tables. 93C NQ is the number of radial points in reciprocal space. 94C EXPAND is a factor to expand some array sizes 95C Q2CUT is the required planewave cutoff for orbital expansion 96C (in Ry if lengths are in Bohr). 97C GWBYDR is the width of a gaussian used to neutralize the charge 98C distributions for operat='U', in units of the radial interval 99C FFTOL is the tolerance for considering equal the radial part of 100C two orbitals. 101C TINY is a small number to avoid a zero denominator 102 INTEGER, save :: NRTAB = 1024 103 INTEGER, PARAMETER :: NQ = 1024 104 real(dp), PARAMETER :: EXPAND = 1.20_dp 105 real(dp), PARAMETER :: Q2CUT = 2.50e3_dp 106 real(dp), PARAMETER :: GWBYDR = 1.5_dp 107 real(dp), PARAMETER :: FFTOL = 1.e-8_dp 108 real(dp), PARAMETER :: TINY = 1.e-12_dp 109 CHARACTER(LEN=*), PARAMETER :: MYNAME = 'MATEL' 110C ------------------------------------------------------------------- 111 112C Internal variable types and dimensions ---------------------------- 113 INTEGER :: 114 & I, IF1, IF2, IFF, IFFY, IFLM1, IFLM2, 115 & IG, IOPER, IQ, IR, IX, 116 & JF1, JF2, JFF, JFFR, JFFY, JFLM1, JFLM2, JLM, 117 & JG1, JG2, JR, 118 & L, L1, L2, L3, LMAX, 119 & NILM, NILM1, NILM2, NJLM1, NJLM2 120 INTEGER, SAVE :: 121 & MF=0, MFF=0, MFFR=0, MFFY=0, 122 & NF=0, NFF=0, NFFR=0, NFFY=0, NR=NQ 123 124 INTEGER, POINTER, SAVE :: 125 & IFFR(:), ILM(:), ILMFF(:), INDF(:,:), INDFF(:,:,:), 126 & INDFFR(:), INDFFY(:), NLM(:,:) 127 128 real(dp) :: 129 & C, CH(2), CPROP, DFFR0, DFFRMX, DSRDR, ERRF, 130 & FFL(0:NQ), FFQ(0:NQ), FR(0:NQ,2), FQ(0:NQ,2), GAUSS, 131 & Q, R, SR, VR(0:NQ,2), VQ(0:NQ,2), X12(3) 132 133 real(dp), SAVE :: 134 & DQ, DR, DRTAB, GWIDTH, PI, QMAX, RMAX 135 136 real(dp), POINTER, SAVE :: 137 & CFFR(:), DYDR(:,:), F(:,:), FFR(:,:,:), FFY(:,:), Y(:) 138 139 LOGICAL :: 140 & FAR, FOUND, PROPOR 141 142 LOGICAL, SAVE :: 143 & NULLIFIED=.FALSE. 144 145 TYPE(allocDefaults) :: 146 & OLDEFS 147 148 EXTERNAL PROPOR, TIMER 149C ------------------------------------------------------------------- 150C Start time counter 151* CALL TIMER( MYNAME, 1 ) 152 153C Set allocation defaults 154 CALL ALLOC_DEFAULT( OLD=OLDEFS, COPY=.TRUE., SHRINK=.FALSE. ) 155 156C Nullify pointers 157 IF (.NOT.NULLIFIED) THEN 158 NULLIFY( IFFR, ILM, ILMFF, INDF, INDFF, INDFFR, INDFFY, NLM, 159 & CFFR, DYDR, F, FFR, FFY, Y ) 160 CALL RE_ALLOC( INDF, 1, 1, 1, 1, 'INDF', MYNAME ) 161 CALL RE_ALLOC( INDFFY, 0,MFF, 'INDFFY', MYNAME ) 162 INDFFY(0) = 0 163 NULLIFIED = .TRUE. 164 165! Set only at specific times. If matel_nrtab is changed the new 166! value will not take effect unless there is a re-initialization, 167! to avoid corruption of the tables. 168 NRTAB = matel_NRTAB 169 170 ENDIF 171 172C Check if tables must be re-initialized 173 IF ( IG1.LE.0 .OR. IG2.LE.0 ) THEN 174 CALL RESET_SPHER_HARM( ) 175 CALL RESET_RADFFT( ) 176 CALL DE_ALLOC( IFFR, 'IFFR', MYNAME ) 177 CALL DE_ALLOC( ILM, 'ILM', MYNAME ) 178 CALL DE_ALLOC( ILMFF, 'ILMFF', MYNAME ) 179 CALL DE_ALLOC( INDF, 'INDF', MYNAME ) 180 CALL DE_ALLOC( INDFF, 'INDFF', MYNAME ) 181 CALL DE_ALLOC( INDFFR, 'INDFFR', MYNAME ) 182 CALL DE_ALLOC( INDFFY, 'INDFFY', MYNAME ) 183 CALL DE_ALLOC( NLM, 'NLM', MYNAME ) 184 CALL DE_ALLOC( CFFR, 'CFFR', MYNAME ) 185 CALL DE_ALLOC( DYDR, 'DYDR', MYNAME ) 186 CALL DE_ALLOC( F, 'F', MYNAME ) 187 CALL DE_ALLOC( FFR, 'FFR', MYNAME ) 188 CALL DE_ALLOC( FFY, 'FFY', MYNAME ) 189 CALL DE_ALLOC( Y, 'Y', MYNAME ) 190 MF = 0 191 MFF = 0 192 MFFR = 0 193 MFFY = 0 194 NF = 0 195 NFF = 0 196 NFFR = 0 197 NFFY = 0 198 NULLIFIED = .FALSE. 199 GOTO 900 200 ENDIF 201 202C Check argument OPERAT 203 IF ( OPERAT .EQ. 'S' ) THEN 204 IOPER = 1 205 ELSEIF ( OPERAT .EQ. 'T' ) THEN 206 IOPER = 2 207 ELSEIF ( OPERAT .EQ. 'U' ) THEN 208 IOPER = 3 209 ELSEIF ( OPERAT .EQ. 'X' ) THEN 210 IOPER = 4 211 ELSEIF ( OPERAT .EQ. 'Y' ) THEN 212 IOPER = 5 213 ELSEIF ( OPERAT .EQ. 'Z' ) THEN 214 IOPER = 6 215 ELSE 216 call die('MATEL: Invalid value of argument OPERAT') 217 ENDIF 218 219C Check size of orbital index table 220 IF ( MAX(IG1,IG2).GT.UBOUND(INDF,1) .OR. 221 & MAX(IOPER,3).GT.UBOUND(INDF,2) ) THEN 222 CALL RE_ALLOC( INDF, 1, MAX(UBOUND(INDF,1),IG1,IG2), 223 & 1,MAX(UBOUND(INDF,2),IOPER,3), 224 & 'INDF', MYNAME) 225 CALL RE_ALLOC( NLM, 1,MAX(UBOUND(INDF,1),IG1,IG2), 226 & 1,MAX(UBOUND(INDF,2),IOPER,3), 227 & 'NLM', MYNAME ) 228 ENDIF 229 230C Find radial expansion of each orbital ----------------------------- 231 DO I = 1,2 232 IF (I .EQ. 1) THEN 233 IG = IG1 234 FOUND = (INDF(IG,1) .NE. 0) 235 ELSE 236 IG = IG2 237 FOUND = (INDF(IG,IOPER) .NE. 0) 238 ENDIF 239 IF (.NOT.FOUND) THEN 240* CALL TIMER( 'MATEL1', 1 ) 241 PI = 4._dp * ATAN(1._dp) 242C Factor 2 below is because we will expand the product of 243C two orbitals 244 QMAX = 2._dp * SQRT( Q2CUT ) 245 DQ = QMAX / NQ 246 DR = PI / QMAX 247 NR = NQ 248 RMAX = NR * DR 249 DRTAB = RMAX / NRTAB 250 IF ( RCUT(IG) .GT. RMAX ) THEN 251 call die('MATEL: NQ too small for required cutoff.') 252 ENDIF 253 254C Reallocate arrays 255 L = LCUT( IG ) 256 NILM = (L+2)**2 257 IF (NF+NILM .GT. MF) MF = EXPAND * (NF+NILM) 258 CALL RE_ALLOC( F, 0,NQ, 1,MF, 'F', MYNAME ) 259 CALL RE_ALLOC( ILM, 1,MF, 'ILM', MYNAME ) 260 CALL RE_ALLOC( INDFF, 1,MF, 1,MF, 1, MAX(IOPER,3), 261 & 'INDFF', MYNAME ) 262 263C Expand orbital in spherical harmonics 264 IF ((I.EQ.1) .OR. (IOPER.LE.3)) THEN 265 CALL YLMEXP( L, RLYLM, EVALUATE, IG, 0, NQ, RMAX, 266 & NILM, ILM(NF+1:), F(:,NF+1:) ) 267 INDF(IG,1) = NF+1 268 INDF(IG,2) = NF+1 269 INDF(IG,3) = NF+1 270 NLM(IG,1) = NILM 271 NLM(IG,2) = NILM 272 NLM(IG,3) = NILM 273 ELSE 274 IF(IOPER.EQ.4) THEN 275 CALL YLMEXP( L+1, RLYLM, EVALUATE_X, IG, 0, NQ, RMAX, 276 & NILM, ILM(NF+1:), F(:,NF+1:) ) 277 ELSEIF(IOPER.EQ.5) THEN 278 CALL YLMEXP( L+1, RLYLM, EVALUATE_Y, IG, 0, NQ, RMAX, 279 & NILM, ILM(NF+1:), F(:,NF+1:) ) 280 ELSE 281 CALL YLMEXP( L+1, RLYLM, EVALUATE_Z, IG, 0, NQ, RMAX, 282 & NILM, ILM(NF+1:), F(:,NF+1:) ) 283 ENDIF 284 INDF(IG,IOPER) = NF+1 285 NLM(IG,IOPER) = NILM 286 ENDIF 287 288C Store orbital in k-space 289 DO JLM = 1,NILM 290 NF = NF + 1 291 L = LOFILM( ILM(NF) ) 292 CALL RADFFT( L, NQ, RMAX, F(0:NQ,NF), F(0:NQ,NF) ) 293* F(NQ,NF) = 0._dp 294 ENDDO 295 296* CALL TIMER( 'MATEL1', 2 ) 297 ENDIF 298 ENDDO 299C ------------------------------------------------------------------- 300 301C Find radial expansion of overlap ---------------------------------- 302 IF1 = INDF(IG1,1) 303 IF2 = INDF(IG2,IOPER) 304 305 IF ( INDFF(IF1,IF2,IOPER) .EQ. 0 ) THEN 306* CALL TIMER( 'MATEL2', 1 ) 307 308 NILM1 = NLM(IG1,1) 309 NILM2 = NLM(IG2,IOPER) 310 DO IFLM1 = IF1,IF1+NILM1-1 311 DO IFLM2 = IF2,IF2+NILM2-1 312 313C Check interaction range 314 IF (RCUT(IG1)+RCUT(IG2) .GT. RMAX) THEN 315 call die('MATEL: NQ too small for required cutoff.') 316 ENDIF 317 318C Special case for operat='U' 319 IF (OPERAT.EQ.'U') THEN 320 321C Check that charge distributions are spherical 322 IF (ILM(IFLM1).NE.1 .OR. ILM(IFLM2).NE.1) THEN 323 CALL DIE('matel: ERROR: operat=U for L=0 only') 324 ENDIF 325 326C Find L=0, Q=0 components of charge distributions 327C The actual charges are these times Y00=1/sqrt(4*pi) 328 CH(1) = (2._dp*PI)**1.5_dp * F(0,IFLM1) 329 CH(2) = (2._dp*PI)**1.5_dp * F(0,IFLM2) 330 331C Add gaussian neutralizing charge and find potential 332 GWIDTH = GWBYDR * DR 333 DO IQ = 1,NQ 334 Q = IQ * DQ 335 GAUSS = EXP( -0.5_dp*(Q*GWIDTH)**2 ) 336 FQ(IQ,1) = F(IQ,IFLM1) - F(0,IFLM1) * GAUSS 337 FQ(IQ,2) = F(IQ,IFLM2) - F(0,IFLM2) * GAUSS 338 VQ(IQ,1) = FQ(IQ,1) * 4._dp*PI/(Q*Q) 339 VQ(IQ,2) = FQ(IQ,2) * 4._dp*PI/(Q*Q) 340 ENDDO 341 FQ(0,1) = 0._dp 342 FQ(0,2) = 0._dp 343 VQ(0,1) = 0._dp 344 VQ(0,2) = 0._dp 345 346C Return to real space 347 L = 0 348 CALL RADFFT( L, NQ, NQ*PI/RMAX, FQ(0:NQ,1), FR(0:NQ,1) ) 349 CALL RADFFT( L, NQ, NQ*PI/RMAX, FQ(0:NQ,2), FR(0:NQ,2) ) 350 CALL RADFFT( L, NQ, NQ*PI/RMAX, VQ(0:NQ,1), VR(0:NQ,1) ) 351 CALL RADFFT( L, NQ, NQ*PI/RMAX, VQ(0:NQ,2), VR(0:NQ,2) ) 352 353C Subtract neutralizing charge 354 DO IR = 0,NR 355 R = IR * DR 356 GAUSS = EXP(-0.5_dp*(R/GWIDTH)**2) / 357 & (2._dp*PI*GWIDTH**2)**1.5_dp 358 IF (IR.EQ.0) THEN 359 ERRF = SQRT(2._dp/PI) / GWIDTH 360 ELSE 361 ERRF = DERF(SQRT(0.5_dp)*R/GWIDTH) / R 362 ENDIF 363 FR(IR,1) = FR(IR,1) + CH(1) * GAUSS 364 FR(IR,2) = FR(IR,2) + CH(2) * GAUSS 365 VR(IR,1) = VR(IR,1) + CH(1) * ERRF 366 VR(IR,2) = VR(IR,2) + CH(2) * ERRF 367 ENDDO 368 369C Select NRTAB out of NQ points 370 DO IR = 0,NRTAB 371 JR = IR * NR / NRTAB 372 FR(IR,1) = FR(JR,1) 373 FR(IR,2) = FR(JR,2) 374 VR(IR,1) = VR(JR,1) 375 VR(IR,2) = VR(JR,2) 376 ENDDO 377 378 ELSE 379 DO IQ = 0,NQ 380 FQ(IQ,1) = F(IQ,IFLM1) 381 FQ(IQ,2) = F(IQ,IFLM2) 382 ENDDO 383 ENDIF 384 385C Find orbitals convolution by multiplication in k-space 386 C = ( 2.0_dp * PI )**1.5_dp 387 DO IQ = 0,NQ 388 Q = IQ * DQ 389 FFQ(IQ) = C * FQ(IQ,1) * FQ(IQ,2) 390 IF ( OPERAT .EQ. 'T' ) THEN 391 FFQ(IQ) = FFQ(IQ) * Q*Q 392 ELSEIF (OPERAT.EQ.'U' .AND. IQ.GT.0) THEN 393 FFQ(IQ) = FFQ(IQ) * 4._dp*PI/(Q*Q) 394 ENDIF 395 ENDDO 396 397C Loop on possible values of l quantum number of product 398 L1 = LOFILM( ILM(IFLM1) ) 399 L2 = LOFILM( ILM(IFLM2) ) 400 CALL RE_ALLOC( IFFR, 0,L1+L2, 'IFFR', MYNAME ) 401 CALL RE_ALLOC( CFFR, 0,L1+L2, 'CFFR', MYNAME ) 402 DO L3 = ABS(L1-L2), L1+L2, 2 403 404C Return to real space 405 CALL RADFFT( L3, NQ, NQ*PI/RMAX, FFQ, FFL ) 406* FFL(NQ) = 0._dp 407 IF (MOD(ABS(L1-L2-L3)/2,2) .NE. 0) THEN 408 DO IR = 0,NR 409 FFL(IR) = - FFL(IR) 410 ENDDO 411 ENDIF 412 413C Divide by R**L 414 IF (L3 .NE. 0) THEN 415 DO IR = 1,NR 416 R = IR * DR 417 FFL(IR) = FFL(IR) / R**L3 418 ENDDO 419C Parabolic extrapolation to R=0 420 FFL(0) = ( 4.0_dp * FFL(1) - FFL(2) ) / 3.0_dp 421 ENDIF 422 423C Select NRTAB out of NR points 424 IF (MOD(NR,NRTAB) .NE. 0) 425 & CALL DIE('matel ERROR: NQ must be multiple of NRTAB') 426 DO IR = 0,NRTAB 427 JR = IR * NR / NRTAB 428 FFL(IR) = FFL(JR) 429 ENDDO 430 431C Find if radial function is already in table 432 FOUND = .FALSE. 433 SEARCH: DO JG1 = 1,SIZE(INDF,1) 434 DO JG2 = 1,SIZE(INDF,1) 435 JF1 = INDF(JG1,1) 436 JF2 = INDF(JG2,IOPER) 437 IF (JF1.NE.0 .AND. JF2.NE.0) THEN 438 NJLM1 = NLM(JG1,1) 439 NJLM2 = NLM(JG2,IOPER) 440 DO JFLM1 = JF1,JF1+NJLM1-1 441 DO JFLM2 = JF2,JF2+NJLM2-1 442 JFF = INDFF(JFLM1,JFLM2,IOPER) 443 IF (JFF .NE. 0) THEN 444 DO JFFY = INDFFY(JFF-1)+1, INDFFY(JFF) 445 JFFR = INDFFR(JFFY) 446 IF ( PROPOR(NRTAB,FFL(1),FFR(1,1,JFFR), 447 & FFTOL,CPROP) ) THEN 448 FOUND = .TRUE. 449 IFFR(L3) = JFFR 450 CFFR(L3) = CPROP 451 EXIT SEARCH 452 ENDIF 453 ENDDO 454 ENDIF 455 ENDDO 456 ENDDO 457 ENDIF 458 ENDDO 459 ENDDO SEARCH 460 461C Store new radial function 462 IF (.NOT.FOUND) THEN 463 NFFR = NFFR + 1 464 IF (NFFR .GT. MFFR) THEN 465 MFFR = EXPAND * NFFR 466 CALL RE_ALLOC( FFR, 0,NRTAB, 1, 2, 1, MFFR, 'FFR', 467 & MYNAME ) 468 ENDIF 469 IFFR(L3) = NFFR 470 CFFR(L3) = 1._dp 471 DO IR = 0,NRTAB 472 FFR(IR,1,NFFR) = FFL(IR) 473 ENDDO 474 475C Add neutralizing-charge corrections 476 IF (OPERAT .EQ. 'U') THEN 477 DO IR = 0,NRTAB 478 IF (IR.EQ.0) THEN 479 ERRF = 1._dp / SQRT(PI) / GWIDTH 480 ELSE 481 R = IR * DRTAB 482 ERRF = DERF(0.5_dp*R/GWIDTH) / R 483 ENDIF 484 FFR(IR,1,NFFR) = 485 & FFR(IR,1,NFFR) 486 & - CH(1) * CH(2) * ERRF 487 & + CH(1) * VR(IR,2) + CH(2) * VR(IR,1) 488 & - 2._dp*PI*GWIDTH**2 * CH(1) * FR(IR,2) 489 & - 2._dp*PI*GWIDTH**2 * CH(2) * FR(IR,1) 490 ENDDO 491 ENDIF 492 493C Setup spline interpolation 494!! Force derivative, rather than second derivative, to zero 495!! DFFR0 = HUGE(1.0_dp) 496 DFFR0 = 0.0_dp 497 DFFRMX = 0.0_dp 498 CALL SPLINE( RMAX/NRTAB, FFR(0:NRTAB,1,NFFR), NRTAB+1, 499 & DFFR0, DFFRMX, FFR(0:NRTAB,2,NFFR) ) 500 ENDIF 501 ENDDO 502 503C Reallocate some arrays 504 NILM = (L1+L2+1)**2 505 IF (NFFY+NILM .GT. MFFY) THEN 506 MFFY = EXPAND * (NFFY+NILM) 507 CALL RE_ALLOC( FFY, 1, 1, 1,MFFY, 'FFY', MYNAME ) 508 CALL RE_ALLOC( ILMFF, 1,MFFY, 'ILMFF', MYNAME ) 509 CALL RE_ALLOC( INDFFR, 1,MFFY, 'INDFFR', MYNAME ) 510 ENDIF 511 CALL RE_ALLOC( Y, 1, NILM, 'Y', MYNAME, .FALSE. ) 512 CALL RE_ALLOC( DYDR, 1, 3, 1, NILM, 'DYDR', MYNAME, .FALSE. ) 513 514C Expand the product of two spherical harmonics (SH) also in SH 515 CALL YLMEXP( L1+L2, RLYLM, YLMYLM, ILM(IFLM1), ILM(IFLM2), 516 & 1, 1, 1.0_dp, NILM, ILMFF(NFFY+1:), 517 & FFY(:,NFFY+1:)) 518 519C Loop on possible lm values of orbital product 520 DO I = 1,NILM 521 NFFY = NFFY + 1 522 JLM = ILMFF(NFFY) 523 L3 = LOFILM( JLM ) 524 INDFFR(NFFY) = IFFR(L3) 525 FFY(1,NFFY) = FFY(1,NFFY) * CFFR(L3) 526 ENDDO 527 528 NFF = NFF + 1 529 IF (NFF .GT. MFF) THEN 530 MFF = EXPAND * NFF 531 CALL RE_ALLOC( INDFFY, 0,MFF, 'INDFFY', MYNAME ) 532 ENDIF 533 INDFFY(NFF) = NFFY 534 INDFF(IFLM1,IFLM2,IOPER) = NFF 535 ENDDO 536 ENDDO 537 538* CALL TIMER( 'MATEL2', 2 ) 539 ENDIF 540C End of initialization section ------------------------------------- 541 542C Find value of matrix element and its gradient --------------------- 543* CALL TIMER( 'MATEL3', 1 ) 544 545C Initialize output 546 S12 = 0.0_dp 547 DSDR(1) = 0.0_dp 548 DSDR(2) = 0.0_dp 549 DSDR(3) = 0.0_dp 550 551C Avoid R12=0 552 X12(1) = R12(1) 553 X12(2) = R12(2) 554 X12(3) = R12(3) 555 R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) ) 556 IF (R .LT. TINY) THEN 557 X12(3) = TINY 558 R = SQRT( X12(1)*X12(1) + X12(2)*X12(2) + X12(3)*X12(3) ) 559 ENDIF 560 561C Find if orbitals are far (out of range) 562 IF (R .GT. RCUT(IG1)+RCUT(IG2) ) THEN 563 FAR = .TRUE. 564 IF (OPERAT.EQ.'U') THEN 565 IF1 = INDF(IG1,1) 566 IF2 = INDF(IG2,IOPER) 567 IFF = INDFF(IF1,IF2,IOPER) 568 IFFY = INDFFY(IFF) 569 JFFR = INDFFR(IFFY) 570 S12 = FFR(NRTAB,1,JFFR) * RMAX / R 571 DO IX = 1,3 572 DSDR(IX) = - S12 * R12(IX) / R**2 573 ENDDO 574 ENDIF 575 ELSE 576 FAR = .FALSE. 577 ENDIF 578 579C Find spherical harmonics times R**L 580 IF (.NOT.FAR) THEN 581 IF1 = INDF(IG1,1) 582 IF2 = INDF(IG2,IOPER) 583 NILM1 = NLM(IG1,1) 584 NILM2 = NLM(IG2,IOPER) 585 DO IFLM1 = IF1,IF1+NILM1-1 586 DO IFLM2 = IF2,IF2+NILM2-1 587 IFF = INDFF(IFLM1,IFLM2,IOPER) 588 LMAX = 0 589 DO IFFY = INDFFY(IFF-1)+1, INDFFY(IFF) 590 JLM = ILMFF(IFFY) 591 LMAX = MAX( LMAX, LOFILM(JLM) ) 592 ENDDO 593 CALL RLYLM( LMAX, X12, Y, DYDR ) 594 595C Interpolate radial functions and obtain SH expansion 596 DO IFFY = INDFFY(IFF-1)+1, INDFFY(IFF) 597 JFFR = INDFFR(IFFY) 598 CALL SPLINT( RMAX/NRTAB, FFR(0:NRTAB,1,JFFR), 599 & FFR(0:NRTAB,2,JFFR), NRTAB+1, R, SR, DSRDR ) 600 JLM = ILMFF(IFFY) 601 S12 = S12 + SR * FFY(1,IFFY) * Y(JLM) 602 DO IX = 1,3 603 DSDR(IX) = DSDR(IX) + 604 & DSRDR * FFY(1,IFFY) * Y(JLM) * X12(IX) / R + 605 & SR * FFY(1,IFFY) * DYDR(IX,JLM) 606 ENDDO 607 ENDDO 608 ENDDO 609 ENDDO 610 ENDIF 611 612* CALL TIMER( 'MATEL3', 2 ) 613C ------------------------------------------------------------------- 614 615C Restore allocation defaults 616 900 CONTINUE 617 CALL ALLOC_DEFAULT( RESTORE=OLDEFS ) 618 619C Stop time counter 620* CALL TIMER( MYNAME, 2 ) 621 622!------------------------ Internal procedures 623 624 END SUBROUTINE new_MATEL 625 end module m_new_matel 626