← Back to branch summary

~ubuntu-branches/ubuntu/utopic/nwchem/utopic

« back to all changes in this revision

Viewing changes to src/nwpw/nwpwlib/utilities/paw_utilities/nwpw_gamma_function.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Daniel Leidert, Andreas Tille, Michael Banck
  • Date: 2013-07-04 12:14:55 UTC
  • mfrom: (1.1.2)
  • Revision ID: package-import@ubuntu.com-20130704121455-5tvsx2qabor3nrui
Tags: 6.3-1
http://bugs.debian.org/696361
* New upstream release.
* Fixes anisotropic properties (Closes: #696361).
* New features include:
  + Multi-reference coupled cluster (MRCC) approaches
  + Hybrid DFT calculations with short-range HF 
  + New density-functionals including Minnesota (M08, M11) and HSE hybrid
    functionals
  + X-ray absorption spectroscopy (XAS) with TDDFT
  + Analytical gradients for the COSMO solvation model
  + Transition densities from TDDFT 
  + DFT+U and Electron-Transfer (ET) methods for plane wave calculations
  + Exploitation of space group symmetry in plane wave geometry optimizations
  + Local density of states (LDOS) collective variable added to Metadynamics
  + Various new XC functionals added for plane wave calculations, including
    hybrid and range-corrected ones
  + Electric field gradients with relativistic corrections 
  + Nudged Elastic Band optimization method
  + Updated basis sets and ECPs 

[ Daniel Leidert ]
* debian/watch: Fixed.

[ Andreas Tille ]
* debian/upstream: References

[ Michael Banck ]
* debian/upstream (Name): New field.
* debian/patches/02_makefile_flags.patch: Refreshed.
* debian/patches/06_statfs_kfreebsd.patch: Likewise.
* debian/patches/07_ga_target_force_linux.patch: Likewise.
* debian/patches/05_avoid_inline_assembler.patch: Removed, no longer needed.
* debian/patches/09_backported_6.1.1_fixes.patch: Likewise.
* debian/control (Build-Depends): Added gfortran-4.7 and gcc-4.7.
* debian/patches/10_force_gcc-4.7.patch: New patch, explicitly sets
  gfortran-4.7 and gcc-4.7, fixes test suite hang with gcc-4.8 (Closes:
  #701328, #713262).
* debian/testsuite: Added tests for COSMO analytical gradients and MRCC.
* debian/rules (MRCC_METHODS): New variable, required to enable MRCC methods.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwpw/nwpwlib/utilities/paw_utilities/nwpw_gamma_function.F
1
 
c $Id: nwpw_gamma_function.F 21411 2011-11-05 06:41:27Z d3y133 $
2
 
 
3
 
!  *************************************************************************
4
 
!  
5
 
!        This is a  modified version of Log(Gamma) function
6
 
!        program from Num. Rec.
7
 
!       Serves as backup if intrinic Gamma function
8
 
!       is not available
9
 
!  ************************************************************************
10
 
 
11
 
      DOUBLE PRECISION FUNCTION nwpw_LN_GAMMA (XX)
12
 
      implicit none
13
 
      DOUBLE PRECISION XX
14
 
      INTEGER J
15
 
      DOUBLE PRECISION SER
16
 
      DOUBLE PRECISION STP
17
 
      DOUBLE PRECISION TMP
18
 
      DOUBLE PRECISION X
 
1
*     ******************************************************
 
2
*     *                                                    *
 
3
*     *             nwpw_SpecialKummer                     *
 
4
*     *                                                    *
 
5
*     ******************************************************
 
6
*
 
7
*     Calculates a special case of the Kummer confluent hypergeometric 
 
8
*     function, M(n+1/2,l+3/2,z) for z .LE. 0
 
9
*
 
10
*     This function was created by  Marat Valiev, and  modified by Eric Bylaska.
 
11
*     See Abramowitz and Stegun for the formulas used in this function.
 
12
*
 
13
      real*8 function nwpw_SpecialKummer(n,l,z)
 
14
      implicit none
 
15
      integer n,l
 
16
      real*8  z
 
17
 
 
18
*     *** local variables ***
 
19
      real*8 eps
 
20
      parameter (eps=1.0d-16)
 
21
 
 
22
      integer i
 
23
      real*8 a,b,m1,m3,s
 
24
 
 
25
*     **** external functions ****
 
26
      real*8   nwpw_gamma,nwpw_gammp
 
27
      external nwpw_gamma,nwpw_gammp
 
28
 
 
29
      nwpw_SpecialKummer = 0.0d0
 
30
 
 
31
*    *** cannot handle positive z ***
 
32
      if (z.gt.0.0d0) then
 
33
         call errquit('nwpw_SpecialKummer:invalid parameter, z>0',0,1)
 
34
      end if
 
35
 
 
36
 
 
37
*    *** solution for z==0 ***
 
38
      if (z.eq.0.0d0) then
 
39
         nwpw_SpecialKummer = 1.0d0
 
40
         return 
 
41
      end if
 
42
 
 
43
*     ***** M(a,a+1,z) = a * (-z)**(-a) * igamma(a,-z) = a * (-z)**(-a) * P(a,-z) *Gamma(a)  where z is real and a = (n+0.5)  ****
 
44
      if (n.eq.l) then
 
45
         nwpw_SpecialKummer = nwpw_gammp(n+0.5d0,(-z))
 
46
     >                       *(n+0.5d0) 
 
47
     >                       *((-z)**((-n)- 0.5d0))
 
48
     >                       *nwpw_gamma(n+0.5d0)
 
49
         return 
 
50
 
 
51
*     ***** M(a,a,z) = exp(z)  where a = (n+0.5)  ****
 
52
      else if (n.eq.(l+1)) then
 
53
         nwpw_SpecialKummer = dexp(z)
 
54
         return 
 
55
      end if
 
56
 
 
57
!     *** do inifinite series for small z
 
58
      if (dabs(z).le.1.0d0) then
 
59
 
 
60
         nwpw_SpecialKummer = 1.0d0
 
61
         s = 1.0d0
 
62
         a = n + 0.5d0
 
63
         b = l + 1.5d0
 
64
         do i=1,10000
 
65
            s = s*(a+i-1)*z/((b+i-1)*i)
 
66
            nwpw_SpecialKummer = nwpw_SpecialKummer + s
 
67
            if (dabs(s).lt.eps) return 
 
68
         end do
 
69
         call errquit("nwpw_SpecialKummer:cannot converge",0,1)
 
70
         return 
 
71
      end if
 
72
 
 
73
      if (n.lt.l) then
 
74
 
 
75
      !*** starting point n=l or b=a+1***
 
76
         a = n + 0.5d0
 
77
         b = n + 1.5d0
 
78
 
 
79
      !*** m1 = M(a,b-1) ***
 
80
      !*** m2 = M(a,b,z) ***
 
81
         m1 = dexp(z)
 
82
         nwpw_SpecialKummer = nwpw_gammp(a,(-z))*a/(-z)**a*nwpw_gamma(a)
 
83
 
 
84
      !**********************************************
 
85
      ! using recursion formula
 
86
      ! z(a-b)M(a,b+1,z)=b(b-1)M(a,b-1,z)+b(1-b-z)M(a,b,z)
 
87
      ! obtain M(1/2,3/2+l  ,z) --> m2
 
88
      !        M(1/2,3/2+l-1,z) --> m2
 
89
      !**********************************************
 
90
         do i=1,l-n
 
91
            m3=(b*(b-1.0d0)*m1+b*(1.0d0-b-z)*nwpw_SpecialKummer)
 
92
     >         /(z*(a-b))
 
93
            b = b + 1
 
94
            m1 = nwpw_SpecialKummer
 
95
            nwpw_SpecialKummer = m3
 
96
         end do
 
97
 
 
98
      else if (n.gt.(l+1)) then
 
99
 
 
100
      !*** starting point n=l+1 or b=a ***
 
101
         a = l + 1.5d0
 
102
         b = l + 1.5d0
 
103
 
 
104
      !*** m1 = M(a-1,b) ***
 
105
      !*** m2 = M(a,a,z) ***
 
106
         m1 = nwpw_gammp(a-1.0d0,(-z))*(a-1.0d0)/(-z)**(a-1.0d0)*
 
107
     >      nwpw_gamma(a-1.0d0)
 
108
         nwpw_SpecialKummer = dexp(z)
 
109
 
 
110
      !**********************************************
 
111
      ! using recursion formula
 
112
      ! aM(a+1,b,z)=(b-a)M(a-1,b,z)+(2a-b+z)M(a,b,z)
 
113
      ! obtain M(n+1/2-1,3/2,z)   --> m1
 
114
      !        M(n+1/2  ,3/2,z)   --> m2
 
115
      !**********************************************
 
116
         do i=1,n-l-1
 
117
            m3 = ((b-a)*m1+(2*a-b+z)*nwpw_SpecialKummer)/a
 
118
            m1 = nwpw_SpecialKummer
 
119
            nwpw_SpecialKummer = m3
 
120
            a = a + 1
 
121
         end do
 
122
      end if
 
123
 
 
124
      return
 
125
      end 
 
126
 
 
127
*     ***************************************   
 
128
*     *                                     *
 
129
*     *          nwpw_ln_gamma              *
 
130
*     *                                     *
 
131
*     ***************************************   
 
132
 
 
133
*        This function computes the Log(Gamma) function
 
134
*       Serves as backup if intrinic Gamma function is not available
 
135
 
 
136
      real*8 function nwpw_ln_gamma(XX)
 
137
      implicit none
 
138
      real*8 XX
 
139
      integer J
 
140
      real*8 SER
 
141
      real*8 STP
 
142
      real*8 TMP
 
143
      real*8 X
19
144
      DOUBLE PRECISION Y, COF(6)
20
145
      SAVE STP, COF
21
 
      DATA COF, STP/ 76.18009172947146D0, -86.50532032941677D0, 
22
 
     1   24.01409824083091D0, -1.231739572450155D0, 
23
 
     2   0.1208650973866179D-2, -.5395239384953D-5, 2.5066282746310005D0
 
146
      DATA COF, STP/ 76.18009172947146d0, -86.50532032941677d0, 
 
147
     1   24.01409824083091d0, -1.231739572450155d0, 
 
148
     2   0.1208650973866179d-2, -.5395239384953d-5, 2.5066282746310005d0
24
149
     3   / 
25
150
      X = XX
26
151
      Y = X
27
 
      TMP = X + 5.5D0
28
 
      TMP = (X + 0.5D0)*DLOG(TMP) - TMP
29
 
      SER = 1.000000000190015D0
 
152
      TMP = X + 5.5d0
 
153
      TMP = (X + 0.5d0)*dlog(TMP) - TMP
 
154
      SER = 1.000000000190015d0
30
155
      DO J = 1, 6
31
 
         Y = Y + 1.0D0
 
156
         Y = Y + 1.0d0
32
157
         SER = SER + COF(J)/Y
33
158
      END DO
34
 
 
35
 
      nwpw_LN_GAMMA = TMP + DLOG(STP*SER/X)
36
 
 
37
 
      RETURN 
38
 
      END 
39
 
 
40
 
!  *************************************************
41
 
!  
42
 
!     Name    :
43
 
!  
44
 
!  
45
 
!     Purpose :
46
 
!  
47
 
!  
48
 
!     Created :
49
 
!  
50
 
!  *************************************************
51
 
      DOUBLE PRECISION FUNCTION nwpw_GAMMA(X)
 
159
      nwpw_ln_gamma = TMP + dlog(STP*SER/X)
 
160
      return 
 
161
      end 
 
162
 
 
163
*     **********************************************
 
164
*     *                                            *
 
165
*     *            nwpw_gamma                      *
 
166
*     *                                            *
 
167
*     **********************************************
 
168
*
 
169
*     This function computes the Gamma function
 
170
*
 
171
      real*8 function nwpw_gamma(X)
52
172
      implicit none
53
 
      DOUBLE PRECISION X
 
173
      real*8 X
54
174
 
55
 
      DOUBLE PRECISION XX
56
 
      DOUBLE PRECISION nwpw_LN_GAMMA
 
175
      real*8 XX
 
176
      real*8 nwpw_ln_gamma
 
177
      external nwpw_ln_gamma
57
178
 
58
179
      XX = X
59
 
      nwpw_GAMMA = DEXP(nwpw_LN_GAMMA(XX))
60
 
 
61
 
      return
62
 
      END 
63
 
 
64
 
 
65
 
 
66
 
!  *************************************************
67
 
!  
68
 
!     Name    :
69
 
!  
70
 
!  
71
 
!     Purpose :
72
 
!  
73
 
!  
74
 
!     Created :
75
 
!  
76
 
!  *************************************************
77
 
      DOUBLE PRECISION FUNCTION nwpw_GAMMP (A, X)
78
 
      implicit none
79
 
      DOUBLE PRECISION A, X
80
 
 
81
 
      DOUBLE PRECISION GAMMCF, GAMSER, GLN
82
 
 
83
 
      IF (X .LT. A+1.0D0) THEN
84
 
 
85
 
         CALL nwpw_GSER(GAMSER, A, X, GLN)
86
 
         nwpw_GAMMP = GAMSER
87
 
 
88
 
      ELSE
89
 
 
90
 
         CALL nwpw_GCF(GAMMCF, A, X, GLN)
91
 
         nwpw_GAMMP = 1.0D0 - GAMMCF
92
 
 
93
 
      ENDIF
94
 
 
95
 
      return
96
 
      END 
97
 
 
98
 
!  *************************************************
99
 
!  
100
 
!     Name    :
101
 
!  
102
 
!  
103
 
!     Purpose :
104
 
!  
105
 
!  
106
 
!     Created :
107
 
!  
108
 
!  *************************************************
109
 
      SUBROUTINE nwpw_GCF(GAMMCF, A, X, GLN)
110
 
      implicit none
111
 
      INTEGER ITMAX
112
 
      DOUBLE PRECISION A, GAMMCF, GLN, X, EPS, FPMIN
113
 
      PARAMETER (ITMAX = 100, EPS = 3.D-16, FPMIN = 1.D-30)
114
 
      DOUBLE PRECISION AN, B, C, D, DEL, H
115
 
      INTEGER I
 
180
      nwpw_gamma = dexp(nwpw_ln_gamma(XX))
 
181
      return
 
182
      end 
 
183
 
 
184
*     **********************************************
 
185
*     *                                            *
 
186
*     *           nwpw_gammp                       *
 
187
*     *                                            *
 
188
*     **********************************************
 
189
*
 
190
*     This function computes the incomplete gamma function P(a,x) = igamma(a,x)/Gamma(a)
 
191
*  
 
192
*                             /x
 
193
*       P(a,x) = 1/Gamma(a) * | exp(-t) * t**(a-1) dt
 
194
*                             /0
 
195
*
 
196
      real*8 function nwpw_gammp(A, X)
 
197
      implicit none
 
198
      real*8 a,x
 
199
      real*8 gammcf,gamser,gln
 
200
 
 
201
      if (x.lt. (a+1.0d0)) then
 
202
         call nwpw_gser(gamser, a, x, gln)
 
203
         nwpw_gammp = gamser
 
204
      else
 
205
         call nwpw_gcf(gammcf,a,x,gln)
 
206
         nwpw_gammp = 1.0d0 - gammcf
 
207
      end if
 
208
      return
 
209
      end 
 
210
 
 
211
*     **********************************************
 
212
*     *                                            *
 
213
*     *               nwpw_gcf                     *
 
214
*     *                                            *
 
215
*     **********************************************
 
216
      subroutine nwpw_gcf(GAMMCF, A, X, GLN)
 
217
      implicit none
 
218
      integer ITMAX
 
219
      real*8 A,GAMMCF,GLN,X,EPS,FPMIN
 
220
      parameter (ITMAX = 100, EPS = 3.0d-16, FPMIN = 1.0d-30)
 
221
      real*8 AN,B,C,D,DEL,H
 
222
      integer I
116
223
 
117
 
      DOUBLE PRECISION nwpw_LN_GAMMA
118
 
      EXTERNAL         nwpw_LN_GAMMA
 
224
      real*8   nwpw_ln_gamma
 
225
      external nwpw_ln_gamma
119
226
 
120
 
      GLN = nwpw_LN_GAMMA(A)
121
 
      B = X + 1.0D0 - A
122
 
      C = 1.0D0/1.D-30
123
 
      D = 1.0D0/B
 
227
      GLN = nwpw_ln_gamma(A)
 
228
      B = X + 1.0d0 - A
 
229
      C = 1.0d0/1.0d-30
 
230
      D = 1.0d0/B
124
231
      H = D
125
232
      DO I = 1, 100
126
233
         AN = -I*(I - A)
127
 
         B = B + 2.0D0
 
234
         B = B + 2.0d0
128
235
         D = AN*D + B
129
 
         IF (DABS(D) .LT. 1.D-30) D = 1.D-30
 
236
         IF (dabs(D) .LT. 1.0d-30) D = 1.0d-30
130
237
         C = B + AN/C
131
 
         IF (DABS(C) .LT. 1.D-30) C = 1.D-30
132
 
         D = 1.0D0/D
 
238
         IF (dabs(C) .LT. 1.0d-30) C = 1.0d-30
 
239
         D = 1.0d0/D
133
240
         DEL = D*C
134
241
         H = H*DEL
135
 
         IF (DABS(DEL - 1.0D0) .LT. 3.D-16) GO TO 1
 
242
         IF (dabs(DEL-1.0d0) .LT. 3.0d-16) GO TO 1
136
243
      END DO
137
 
      PAUSE 'a too large, ITMAX too small in gcf'
138
 
    1 CONTINUE
139
 
      GAMMCF = DEXP((-X) + A*DLOG(X) - GLN)*H
140
 
 
 
244
      call errquit('a too large, ITMAX too small in nwpw_gcf',0,0)
 
245
 
 
246
    1 continue
 
247
      GAMMCF = dexp((-X) + A*dlog(X) - GLN)*H
141
248
      return
142
 
      END 
 
249
      end 
143
250
 
144
 
!  *************************************************
145
 
!  
146
 
!     Name    :
147
 
!  
148
 
!  
149
 
!     Purpose :
150
 
!  
151
 
!  
152
 
!     Created :
153
 
!  
154
 
!  *************************************************
155
 
      SUBROUTINE nwpw_GSER(GAMSER, A, X, GLN)
 
251
*     **********************************************
 
252
*     *                                            *
 
253
*     *             nwpw_gser                      *
 
254
*     *                                            *
 
255
*     **********************************************
 
256
      subroutine nwpw_gser(GAMSER, A, X, GLN)
156
257
      implicit none
157
 
      DOUBLE PRECISION A, X
158
 
      DOUBLE PRECISION GAMSER, GLN
 
258
      real*8 A,X
 
259
      real*8 GAMSER,GLN
159
260
 
160
261
!    *** local variables ***
161
 
      INTEGER ITMAX
162
 
      PARAMETER (ITMAX = 100)
163
 
      DOUBLE PRECISION EPS
164
 
      PARAMETER (EPS = 3.0D-16)
165
 
      INTEGER N
166
 
      DOUBLE PRECISION AP, DEL, SUM
 
262
      integer ITMAX
 
263
      parameter (ITMAX = 100)
 
264
      real*8 EPS
 
265
      parameter (EPS = 3.0d-16)
 
266
      integer N
 
267
      real*8 AP,DEL,SUM
167
268
 
168
 
      DOUBLE PRECISION nwpw_LN_GAMMA
169
 
      EXTERNAL         nwpw_LN_GAMMA
170
 
 
171
 
 
172
 
      GLN = nwpw_LN_GAMMA(A)
173
 
 
174
 
      IF (X .LE. 0.0D0) THEN
175
 
      IF(X.lt.0.0d0) CALL errquit("x < 0 in nwpw_GSER",0,1)
176
 
         GAMSER = 0.0D0
177
 
         RETURN 
178
 
      ENDIF
 
269
      real*8   nwpw_ln_gamma
 
270
      external nwpw_ln_gamma
 
271
 
 
272
      GLN = nwpw_ln_gamma(A)
 
273
 
 
274
      if (X .le. 0.0d0) then
 
275
         if (X.lt.0.0d0) call errquit("x < 0 in nwpw_gser",0,1)
 
276
         GAMSER = 0.0d0
 
277
         return
 
278
      end if
179
279
 
180
280
      AP = A
181
 
      SUM = 1.0D0/A
 
281
      SUM = 1.0d0/A
182
282
      DEL = SUM
183
 
      DO N = 1, 100
184
 
         AP = AP + 1.0D0
 
283
      do N = 1, 100
 
284
         AP = AP + 1.0d0
185
285
         DEL = DEL*X/AP
186
286
         SUM = SUM + DEL
187
 
         IF (DABS(DEL) .LT. DABS(SUM)*3.0D-16) GO TO 1
188
 
 
189
 
      END DO
190
 
 
191
 
      CALL errquit
192
 
     >     ("a too large,ITMAX too small in nwpw_GSER",0,1)
193
 
 
194
 
    1 CONTINUE
195
 
      GAMSER = SUM*DEXP((-X) + A*DLOG(X) - GLN)
 
287
         if (dabs(DEL) .lt. dabs(SUM)*3.0d-16) GO TO 1
 
288
 
 
289
      end do
 
290
 
 
291
      call errquit
 
292
     >     ("a too large,ITMAX too small in nwpw_gser",0,1)
 
293
 
 
294
    1 continue
 
295
      GAMSER = SUM*dexp((-X) + A*dlog(X) - GLN)
196
296
 
197
297
      return 
198
 
      END 
199
 
 
200
 
 
 
298
      end 
201
299