1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18 SUBROUTINE RADLMG(RADNDE,RADWT,NR,RADERR,NRADPT,NUCORB,AA,IPRINT) 19#include "implicit.h" 20 PARAMETER(D0 = 0D0, D1 = 1D0,D2 = 2D0, D3 = 3D0) 21#include "dummy.h" 22#include "maxaqn.h" 23#include "mxcent.h" 24#include "ccom.h" 25#include "nuclei.h" 26#include "priunit.h" 27 DIMENSION NUCORB(NHTYP,2),AA(2,NHTYP,2) 28 DIMENSION RADWT(NRADPT), RADNDE(NRADPT) 29 CHARACTER SPDCAR*1 30c 31C Grid spacing to H and inner grid point to AH 32 IF(IPRINT.GE.3)WRITE(LUPRI,'(A)') '* Grid spacing' 33 NR = 0 34 H = DUMMY 35 AH = 0D0 36 DO LL = 1,NHTYP 37 L = LL-1 38 NBAS=NUCORB(LL,1)+NUCORB(LL,2) 39 IF(NBAS.GT.0) THEN 40 HTMP = GRID_DISERR(RADERR,L) 41 H = MIN(H,HTMP) 42 IF(IPRINT.GE.3) THEN 43 WRITE(LUPRI,'(3X,A1,A,F6.3)') 44 & SPDCAR(L),'-orbitals --> ',HTMP 45 ENDIF 46 ENDIF 47c AH = MAX(AA(1,LL,1),AA(1,LL,2)) 48 AH = MAX(AH,AA(1,LL,1)) 49 ENDDO 50 IF(AH .EQ. 0d0) RETURN 51 EPH = EXP(H) 52 IF(IPRINT.GE.3)WRITE(LUPRI,'(A,F12.3)') ' Value chosen:',H 53C... Inner grid point AA->R transformation. 54 AH = D2*AH 55 IF(IPRINT.GE.3)WRITE(LUPRI,*) 'AH = ',AH 56 RL = ((1.9D0+LOG(RADERR))/D3)-(LOG(AH)/D2) 57 RL = EXP(RL) 58 IF(IPRINT.GE.3) 59 & WRITE(LUPRI,'(A,1P,E12.5)') '* Inner grid point:',RL 60C... Outer point 61 IF(IPRINT.GE.3) WRITE(LUPRI,'(A)') '* Outer point:' 62 RH = D0 63 DO LL = 1,NHTYP 64 L = LL-1 65 AL=DUMMY 66 IF(NUCORB(LL,1).GT.0) AL=AA(2,LL,1) 67 IF(NUCORB(LL,2).GT.0) AL=MIN(AL,AA(2,LL,2)) 68 IF(AL.LT.DUMMY) THEN 69 AL = AL+AL 70 RHTMP = GRID_OUTERR(AL,L,RADERR) 71 RH=MAX(RH,RHTMP) 72 IF(IPRINT.GE.3) THEN 73 WRITE(LUPRI,'(3X,A1,A,F6.3)') 74 & SPDCAR(L),'-orbitals --> ',RHTMP 75 ENDIF 76 ENDIF 77 ENDDO 78 IF(IPRINT.GE.3)WRITE(LUPRI,'(A,F12.3)') ' Value chosen:',RH 79 GRDC = RL/(EPH-D1) 80 IF(IPRINT.GE.3)WRITE(LUPRI,'(A,1P,E12.5)') ' Constant c: ',GRDC 81 NR = NINT(LOG(D1+(RH/GRDC))/H) 82 IF(IPRINT.GE.3)WRITE(LUPRI,'(A,I9)') ' Number of points:',NR 83 IF(NR.GT.NRADPT) CALL QUIT('Too many radial points.') 84 RADNDE(NR) = RL 85 RADWT(NR) = (RL+GRDC)*RL*RL*H 86 DO IR = NR-1,1,-1 87 RADNDE(IR) = (RADNDE(IR+1)+GRDC)*EPH-GRDC 88 RADWT(IR) = (RADNDE(IR)+GRDC)*RADNDE(IR)*RADNDE(IR)*H 89 ENDDO 90 END 91CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 92C C 93 FUNCTION GRID_DISERR(RD,L) 94C C 95C Provide grid spacing h for given angular momentum L C 96C and discretization error RD C 97C C 98C Based on eqs. (17) and (18) of C 99C R. Lindh, P.-Aa. Malmqvist and L. Gagliardi C 100C "Molecular integrals by numerical quadrature", C 101C Theor. Chem. Acc. 106 (2001) 178-187 C 102C C 103C The array CF(4,L) contains coefficients of a 3rd order C 104C polynomial fit to provide start values for the C 105C determination of H by a Newton-Raphson search. C 106C C 107C Written by T. Saue July 2002 C 108C C 109CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 110 111 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 112#include "priunit.h" 113 PARAMETER(ACC=1.0D-5) 114 PARAMETER(DP5=0.5D0, D2=2.0D0) 115 PARAMETER (PI = 3.14159 26535 89793 D00) 116 PARAMETER(MXIT=20) 117 DIMENSION CF(4,0:4) 118 DATA CF/0.91570D0,0.78806D-1,0.28056D-2,3.4197D-05, 119 & 0.74912D0,0.61502D-1,0.21558D-2,2.6100D-05, 120 & 0.65449D0,0.52322D-1,0.18217D-2,2.2004D-05, 121 & 0.59321D0,0.46769D-1,0.16261D-2,1.9649D-05, 122 & 0.55125D0,0.43269D-1,0.15084D-2,1.8270D-05/ 123C 124C Initialization 125C 126chj start 127c FAC = SQRT(D2)*D2*D2 128 IFAC = 1 129 DO I = 1,L 130c FAC = FAC*D2 131 IFAC = IFAC*(2*I+1) 132 ENDDO 133c FAC = FAC/IFAC 134 FAC = IFAC 135 FAC = SQRT(D2) * D2**(L+2) / IFAC 136chj end 137 LM = MIN(L,4) 138 RDLOG = LOG(RD) 139 GRID_DISERR = POLVAL(3,CF(1,LM),RDLOG) 140 HTLOG = LOG(GRID_DISERR) 141C Newton-Raphson search 142 DO IT = 1,MXIT 143 PIH = PI/GRID_DISERR 144 PIHL = PIH 145 PIEX = PI*PIH*DP5 146 DO I = 1,L 147 PIHL = PIHL*PIH 148 ENDDO 149 U0 = FAC*PIHL*EXP(-PIEX) 150 U1 = U0*((PIEX/GRID_DISERR)-(L+1)/PIH) 151 F0 = LOG(U0)-RDLOG 152 F1 = GRID_DISERR*U1/U0 153 DX = F0/F1 154 HTLOG = HTLOG - DX 155 GRID_DISERR = EXP(HTLOG) 156 IF (ABS(DX).LT.ACC) RETURN 157 ENDDO 158 159 WRITE (LUPRI,*) "Error in GRID_DISERR in dft_grid.F" 160 WRITE (LUPRI,*) "RD, L =", RD,L 161 CALL QUIT('Error in GRID_DISERR in dft_grid.F') 162 163 END 164 165CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 166C C 167 FUNCTION GRID_OUTERR(AL,L,RD) 168C C 169C Provide outer grid point for given angular momentum L C 170C outer exponent AL and discretization error RD C 171C C 172C Based on eq. (19) of C 173C R. Lindh, P.-Aa. Malmqvist and L. Gagliardi C 174C "Molecular integrals by numerical quadrature", C 175C Theor. Chem. Acc. 106 (2001) 178-187 C 176C C 177C The variable U = AL*R*R is found by a Newton-Raphson search. C 178C C 179C Written by T. Saue July 2002 C 180C C 181CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 182 183 IMPLICIT DOUBLE PRECISION(A-H,O-Z) 184#include "priunit.h" 185 PARAMETER(ACC=1.0D-6) 186 PARAMETER(D0=0.0D0,D1=1.0D0,D2=2.0D0,TI=1.0D1) 187 PARAMETER (PI = 3.14159 26535 89793 D00) 188 PARAMETER(MXIT=8) 189C 190C Initialization 191C 192 TOLEN = D2 193 FAC = D1 194 DO I = 1,L 195 TOLEN = TOLEN*D2 196 FAC = FAC*(2*I+1) 197 ENDDO 198 EXPL = (2*L+1)/D2 199 A = SQRT(PI)*FAC/TOLEN 200 ALN = LOG(A) 201 RLN = LOG(RD) 202 U = 35.0D0 203C Newton-Raphson search 204 DO IT = 1,MXIT 205 F0HLN = ALN+EXPL*LOG(U)-U-RLN 206 F1HLN = EXPL/U-D1 207 DX = F0HLN/F1HLN 208 U = U - DX 209 IF(ABS(DX).LT.ACC) THEN 210 GRID_OUTERR = SQRT(U/AL) 211 RETURN 212 ENDIF 213 ENDDO 214 215 WRITE (LUPRI,*) "Error in GRID_OUTERR in dft_grid.F" 216 WRITE (LUPRI,*) "RD, AL, L =", RD,AL,L 217 CALL QUIT('Error in GRID_OUTERR in dft_grid.F') 218 219 RETURN 220 END 221 FUNCTION GETMAXL() 222#include "implicit.h" 223#include "maxaqn.h" 224#include "ccom.h" 225 INTEGER GETMAXL 226 GETMAXL = NHTYP 227 END 228 SUBROUTINE NUCBAS(NUCORB,AA,IPRINT) 229C*********************************************************************** 230C 231C Extract basis information for all centers 232C 233C Written by T.Saue March 12 2001 234C 235C*********************************************************************** 236#include "implicit.h" 237#include "priunit.h" 238 PARAMETER(D0=0.0D0) 239C 240#include "maxaqn.h" 241#include "mxcent.h" 242#include "maxorb.h" 243#include "aovec.h" 244#include "dummy.h" 245C 246#include "nuclei.h" 247#include "ccom.h" 248#include "shells.h" 249#include "primit.h" 250 CHARACTER SPDCAR*1 251 DIMENSION NUCORB(NHTYP,2,NUCIND),AA(2,NHTYP,2,NUCIND) 252C 253C Initialize 254C 255 NDIM = 2*(NHTYP)*NUCIND 256 CALL IZERO(NUCORB,NDIM) 257C 258 JCENT = 0 259 JPRIM = -1 260 JC = -1 261 JLVAL = -1 262 DO ISHELL = 1,KMAX 263 ICENT = NCENT(ISHELL) 264 IF(ICENT.NE.JCENT) THEN 265 JCENT = ICENT 266 JLVAL = 0 267 ENDIF 268 IC = LCLASS(ISHELL) 269 IF(IC.NE.JC) THEN 270 JC = IC 271 JLVAL = 0 272 ENDIF 273 ILVAL = NHKT(ISHELL) 274 IF(ILVAL.NE.JLVAL) THEN 275 JLVAL = ILVAL 276 NUCORB(ILVAL,IC,ICENT) = 0 277 AA(1,ILVAL,IC,ICENT)=D0 278 AA(2,ILVAL,IC,ICENT)=DUMMY 279 ENDIF 280 NUCORB(ILVAL,IC,ICENT)=NUCORB(ILVAL,IC,ICENT)+1 281 IPRIM = JSTRT(ISHELL) 282 IF(IPRIM.NE.JPRIM) THEN 283 JPRIM = IPRIM 284 NPRIM = NUCO(ISHELL) 285 DO IEXP = 1,NPRIM 286 A=PRIEXP(IPRIM+IEXP) 287 AA(1,ILVAL,IC,ICENT)=MAX(AA(1,ILVAL,IC,ICENT),A) 288 AA(2,ILVAL,IC,ICENT)=MIN(AA(2,ILVAL,IC,ICENT),A) 289 ENDDO 290 ENDIF 291 ENDDO 292C 293 IF(IPRINT.GE.2) THEN 294 CALL HEADER('NUCORB:Basis set information:',-1) 295 DO I = 1,NUCIND 296 WRITE(LUPRI,'(/A,A4,A/)') '*** Center: ',NAMN(I),' ***' 297 WRITE(LUPRI,'(2X,A)') '* Large components:' 298 IC = 1 299 DO LL = 1,NHTYP 300 IF(NUCORB(LL,IC,I).GT.0) THEN 301 L=LL-1 302 WRITE(LUPRI,'(3X,A1,A,I6,2(3X,A,E12.5))') 303 & SPDCAR(L),'-orbitals:',NUCORB(LL,IC,I), 304 & 'Alpha_H :',AA(1,LL,IC,I), 305 & 'Alpha_L :',AA(2,LL,IC,I) 306 ENDIF 307 ENDDO 308 WRITE(LUPRI,'(2X,A)') '* Small components:' 309 IC = 2 310 DO LL = 1,NHTYP 311 IF(NUCORB(LL,IC,I).GT.0) THEN 312 L=LL-1 313 WRITE(LUPRI,'(3X,A1,A,I6,2(3X,A,E12.5))') 314 & SPDCAR(L),'-orbitals:',NUCORB(LL,IC,I), 315 & 'Alpha_H: ',AA(1,LL,IC,I), 316 & 'Alpha_L: ',AA(2,LL,IC,I) 317 ENDIF 318 ENDDO 319 ENDDO 320 ENDIF 321C 322 RETURN 323 END 324c Trond: polynomial evaluation. 325 FUNCTION POLVAL(NORDER,B,XVAL) 326#include "implicit.h" 327 DIMENSION B(NORDER+1) 328 POLVAL = B(1) 329 XBUF = 1 330 DO 10 I = 2,(NORDER+1) 331 XBUF = XBUF*XVAL 332 POLVAL = POLVAL + XBUF*B(I) 333 10 CONTINUE 334 END 335