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