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