1!  Dalton, a molecular electronic structure program
2!  Copyright (C) by the authors of Dalton.
3!
4!  This program is free software; you can redistribute it and/or
5!  modify it under the terms of the GNU Lesser General Public
6!  License version 2.1 as published by the Free Software Foundation.
7!
8!  This program is distributed in the hope that it will be useful,
9!  but WITHOUT ANY WARRANTY; without even the implied warranty of
10!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
11!  Lesser General Public License for more details.
12!
13!  If a copy of the GNU LGPL v2.1 was not distributed with this
14!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
15C
16C /* Deck erf */
17C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&
18      FUNCTION ERF(X)
19C*****************************************************************************
20C
21C     Written by Jesper Kielberg Pedersen, Mar. 2003
22C
23C     Purpose : Calculate value of error-function in point X
24C               Based on code from numerical recipies.
25C
26C*****************************************************************************
27#include "implicit.h"
28      PARAMETER (D0 = 0.0D0, DP5 = 0.5D0)
29C
30      IF(X.LT.D0)THEN
31        ERF = -GAMMP2(DP5,X**2)
32      ELSE
33        ERF =  GAMMP2(DP5,X**2)
34      ENDIF
35      RETURN
36      END
37C
38C*****************************************************************************
39C
40      FUNCTION GAMMP2(A,X)
41#include "implicit.h"
42      PARAMETER (D0 = 0.0D0, D1 = 1.0D0)
43      IF (X .LT. D0 .OR. A.LE. D0) CALL QUIT('BAD ARGUMENTS IN GAMMP2')
44      IF (X.LT.(A+D1)) THEN
45        CALL GSER2(GAMSER,A,X)
46        GAMMP2=GAMSER
47      ELSE
48        CALL GCF2(GAMMCF,A,X)
49        GAMMP2=D1-GAMMCF
50      ENDIF
51      RETURN
52      END
53C
54C*****************************************************************************
55C
56      SUBROUTINE GSER2(GAMSER,A,X)
57#include "implicit.h"
58#include "priunit.h"
59      PARAMETER (ITMAX=100, EPS=3.0D-7, D0=0.0D0, D1 = 1.0D0)
60      GLN=GAMMLN2(A)
61      IF(X.LE.D0)THEN
62        IF(X.LT.D0)
63     >    WRITE(LUPRI,'(//3X//,A)') 'WARNING IN ERF : X < 0 IN GSER'
64        GAMSER=D0
65        RETURN
66      ENDIF
67      AP=A
68      SUM=D1/A
69      DEL=SUM
70      DO 11 N=1,ITMAX
71        AP=AP+D1
72        DEL=DEL*X/AP
73        SUM=SUM+DEL
74        IF(ABS(DEL).LT.ABS(SUM)*EPS)GOTO 1
7511    CONTINUE
76      IF(X.LT.D0)
77     >  WRITE(LUPRI,'(//3X//,A)') 'A TOO LARGE, ITMAX TOO SMALL IN GSER'
781     GAMSER=SUM*EXP(-X+A*LOG(X)-GLN)
79      RETURN
80      END
81C
82C*****************************************************************************
83C
84      SUBROUTINE GCF2(GAMMCF,A,X)
85#include "implicit.h"
86#include "priunit.h"
87      PARAMETER (ITMAX=100,EPS=3.0D-7,FPMIN=1.0D-30, D1 = 1.0D0,
88     >           D2 = 2.0D0)
89      GLN=GAMMLN(A)
90      B=X+D1-A
91      C=D1/FPMIN
92      D=D1/B
93      H=D
94      DO 11 I=1,ITMAX
95        AN=-I*(I-A)
96        B=B+D2
97        D=AN*D+B
98        IF(ABS(D).LT.FPMIN)D=FPMIN
99        C=B+AN/C
100        IF(ABS(C).LT.FPMIN)C=FPMIN
101        D=D1/D
102        DEL=D*C
103        H=H*DEL
104        IF(ABS(DEL-D1).LT.EPS)GOTO 1
10511    CONTINUE
106      WRITE(LUPRI,'(//3X//,A)') 'A TOO LARGE, ITMAX TOO SMALL IN GCF2'
1071     GAMMCF=EXP(-X+A*LOG(X)-GLN)*H
108      RETURN
109      END
110C
111C*****************************************************************************
112C
113      FUNCTION GAMMLN2(XX)
114#include "implicit.h"
115      DOUBLE PRECISION STP,COF(6)
116      SAVE COF,STP
117      DATA COF,STP/76.18009172947146D0,-86.50532032941677D0,
118     *24.01409824083091D0,-1.231739572450155D0,.1208650973866179D-2,
119     *-.5395239384953D-5,2.5066282746310005D0/
120      X=XX
121      Y=X
122      TMP=X+5.5D0
123      TMP=(X+0.5D0)*LOG(TMP)-TMP
124      SER=1.000000000190015D0
125      DO 11 J=1,6
126        Y=Y+1.0D0
127        SER=SER+COF(J)/Y
12811    CONTINUE
129      GAMMLN2=TMP+LOG(STP*SER/X)
130      RETURN
131      END
132C*****************************************************************************
133      FUNCTION DERF(X)
134C*****************************************************************************
135C
136C     Written by Jesper Kielberg Pedersen, Mar. 2003
137C
138C     Purpose : Calculate value of error-function in point X
139C
140C               Based on f90-code from :
141C               Naval Surface Warfare Center Mathematical Library
142C               (http://www.math.iastate.edu/burkardt/f_src/nswc/nswc.html)
143C
144C*****************************************************************************
145#include "implicit.h"
146      PARAMETER (C = .564189583547756D0, ONE = 1.0D0, HALF = 0.5D0,
147     >           ZERO = 0.0D0)
148      DOUBLE PRECISION A(5), B(3), P(8), Q(8), R(5), S(4)
149      SAVE A,B,P,Q,R,S
150      DATA A / .771058495001320D-04, -.133733772997339D-02,
151     >         .323076579225834D-01,  .479137145607681D-01,
152     >         .128379167095513D+00 /
153      DATA B / .301048631703895D-02,  .538971687740286D-01,
154     >         .375795757275549D+00 /
155      DATA P / -1.36864857382717D-07, 5.64195517478974D-01,
156     >          7.21175825088309D+00, 4.31622272220567D+01,
157     >          1.52989285046940D+02, 3.39320816734344D+02,
158     >          4.51918953711873D+02, 3.00459261020162D+02 /
159      DATA Q /  1.00000000000000D+00, 1.27827273196294D+01,
160     >          7.70001529352295D+01, 2.77585444743988D+02,
161     >          6.38980264465631D+02, 9.31354094850610D+02,
162     >          7.90950925327898D+02, 3.00459260956983D+02 /
163      DATA R /  2.10144126479064D+00, 2.62370141675169D+01,
164     >          2.13688200555087D+01, 4.65807828718470D+00,
165     >          2.82094791773523D-01 /
166      DATA S /  9.41537750555460D+01, 1.87114811799590D+02,
167     >          9.90191814623914D+01, 1.80124575948747D+01 /
168C
169C
170      AX = ABS(X)
171C
172      IF (AX .LE. HALF) THEN
173        T = X*X
174        TOP = ((((A(1)*T + A(2))*T + A(3))*T + A(4))*T + A(5)) + ONE
175        BOT = ((B(1)*T + B(2))*T + B(3))*T + ONE
176        FN_VAL = X*(TOP/BOT)
177      ELSE IF (AX .LE. 4.0D0) THEN
178        TOP = ((((((P(1)*AX + P(2))*AX + P(3))*AX + P(4))*AX + P(5))*AX
179     >        + P(6))*AX + P(7))*AX + P(8)
180        BOT = ((((((Q(1)*AX + Q(2))*AX + Q(3))*AX + Q(4))*AX + Q(5))*AX
181     >        + Q(6))*AX + Q(7))*AX + Q(8)
182        FN_VAL = HALF + (HALF - DEXP(-X*X)*TOP/BOT)
183        IF (X .LT. ZERO) FN_VAL = -FN_VAL
184      ELSE IF (AX .LT. 5.8D0) THEN
185        X2 = X*X
186        T = ONE / X2
187        TOP = (((R(1)*T + R(2))*T + R(3))*T + R(4))*T + R(5)
188        BOT = (((S(1)*T + S(2))*T + S(3))*T + S(4))*T + ONE
189        FN_VAL = (C - TOP/(X2*BOT)) / AX
190        FN_VAL = HALF + (HALF - DEXP(-X2)*FN_VAL)
191        IF (X .LT. ZERO) FN_VAL = -FN_VAL
192      ELSE
193        FN_VAL = SIGN(ONE, X)
194      END IF
195C
196      DERF = FN_VAL
197      RETURN
198      END
199C*****************************************************************************
200      FUNCTION DERFC(X)
201C*****************************************************************************
202C
203C     Written by Jesper Kielberg Pedersen, Mar. 2003
204C
205C     Purpose : Calculate value of complimentary error-function in point X
206C
207C               Based on f90-code from :
208C               Naval Surface Warfare Center Mathematical Library
209C               (http://www.math.iastate.edu/burkardt/f_src/nswc/nswc.html)
210C
211C*****************************************************************************
212#include "implicit.h"
213      DERFC = 1.0D0 - DERF(X)
214      RETURN
215      END
216C*****************************************************************************
217