← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_camb88.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/xc/xc_camb88.F
1
 
#if defined(FUJITSU_VPP)
2
 
!ocl scalar
3
 
#endif
4
 
#ifndef SECOND_DERIV
5
 
      Subroutine xc_camb88(tol_rho, fac, lfac, nlfac, rho, delrho, 
6
 
     &                      Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
7
 
#else
8
 
      Subroutine xc_camb88_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 
9
 
     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
10
 
     &                         qwght,ldew,func)
11
 
#endif
12
 
c
13
 
C$Id: xc_camb88.F,v 1.4 2008/12/06 23:08:58 niri Exp $
14
 
c
15
 
c     Coulomb attenuated Becke88 functional
16
 
c
17
 
      implicit none
18
 
c      
19
 
#include "dft2drv.fh"
20
 
c
21
 
      double precision tol_rho, fac, Ex
22
 
      integer nq, ipol
23
 
      logical lfac, nlfac,ldew
24
 
      double precision func(*)  ! value of the functional [output]
25
 
c
26
 
c     Charge Density
27
 
c
28
 
      double precision rho(nq,ipol*(ipol+1)/2)
29
 
c
30
 
c     Charge Density Gradient
31
 
c
32
 
      double precision delrho(nq,3,ipol)
33
 
c
34
 
c     Quadrature Weights
35
 
c
36
 
      double precision qwght(nq)
37
 
c
38
 
c     Sampling Matrices for the XC Potential
39
 
c
40
 
      double precision Amat(nq,ipol), Cmat(nq,*)
41
 
c
42
 
#ifdef SECOND_DERIV
43
 
c
44
 
c     Second Derivatives of the Exchange Energy Functional
45
 
c
46
 
      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
47
 
      double precision A2tmp, C2tmp, C3tmp
48
 
#endif
49
 
c
50
 
      double precision BETA
51
 
      Parameter (BETA = 0.0042D0)
52
 
c
53
 
c References:
54
 
c
55
 
c    Becke, Phys. Rev. A 38, 3098 (1988)
56
 
c    Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993)
57
 
c
58
 
c***************************************************************************
59
 
c
60
 
      integer n
61
 
      double precision arcsinh, darcsinh
62
 
      double precision C, rho13, rho43, gamma, x, g, gdenom, dg,
63
 
     &     dgdenom, t, Etmp, Atmp, Ctmp
64
 
#ifdef SECOND_DERIV
65
 
      double precision rhom23, d2g, d2gdenom
66
 
#endif
67
 
c
68
 
      arcsinh(x)=log(x+dsqrt(1d0+x*x))
69
 
      darcsinh(x)=1d0/dsqrt(1d0+x*x)
70
 
c
71
 
c     Uniform electron gas constant
72
 
c
73
 
      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
74
 
c
75
 
      if (ipol.eq.1) then
76
 
c
77
 
c        ======> SPIN-RESTRICTED <======
78
 
c
79
 
         do 10 n = 1, nq
80
 
            if (rho(n,1).lt.tol_rho) goto 10
81
 
c
82
 
c           Spin alpha:
83
 
c
84
 
            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
85
 
            rho43 = rho13**4 
86
 
            gamma = delrho(n,1,1)*delrho(n,1,1) +
87
 
     &              delrho(n,2,1)*delrho(n,2,1) +
88
 
     &              delrho(n,3,1)*delrho(n,3,1)
89
 
            if (dsqrt(gamma).gt.tol_rho)then
90
 
               gamma = 0.25d0 * gamma
91
 
               x = sqrt(gamma) / rho43
92
 
            else
93
 
               x = 0d0
94
 
            endif
95
 
c
96
 
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
97
 
            g = -BETA*x*x / gdenom
98
 
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
99
 
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
100
 
c
101
 
            Etmp = 0.d0
102
 
            Atmp = 0.d0
103
 
            Ctmp = 0.d0
104
 
            if (lfac) then
105
 
               Etmp = 2d0*rho43*C*fac
106
 
               Atmp = (4d0/3d0)*rho13*C*fac
107
 
            endif
108
 
c
109
 
            if (nlfac) then
110
 
               Etmp = Etmp + 2d0*rho43*g*fac
111
 
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
112
 
            endif
113
 
c
114
 
            if (x.gt.tol_rho) then
115
 
               Ctmp = 0.5d0 * dg / sqrt(gamma) * fac
116
 
            endif
117
 
c
118
 
#ifdef SECOND_DERIV
119
 
            A2tmp = 0d0
120
 
            C2tmp = 0d0
121
 
            C3tmp = 0d0
122
 
            if(lfac) g = g + C           ! Add local contribution back to g
123
 
            rhom23 = rho13 / (0.5d0*rho(n,1))
124
 
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
125
 
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
126
 
     &           + BETA*x*x*d2gdenom/gdenom**2
127
 
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
128
 
c
129
 
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
130
 
            C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac
131
 
            if (x.gt.tol_rho) then
132
 
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
133
 
            endif
134
 
#endif
135
 
c
136
 
#ifdef SECOND_DERIV
137
 
            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
138
 
     &           C2tmp,C3tmp)
139
 
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
140
 
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
141
 
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
142
 
#else
143
 
            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
144
 
#endif
145
 
            Ex = Ex + qwght(n)*Etmp
146
 
            if(ldew) func(n) = func(n) + Etmp
147
 
            Amat(n,1) = Amat(n,1) + Atmp
148
 
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
149
 
 10      continue
150
 
c
151
 
      else
152
 
c
153
 
c        ======> SPIN-UNRESTRICTED <======
154
 
c
155
 
         do 20 n = 1, nq
156
 
            if (dabs(rho(n,1)).lt.tol_rho) goto 20
157
 
            if (dabs(rho(n,2)).lt.tol_rho) goto 25
158
 
c
159
 
c           Spin alpha:
160
 
c
161
 
            rho13 = abs(rho(n,2))**(1.d0/3.d0)*sign(1d0,rho(n,2))
162
 
            rho43 = rho13**4 
163
 
            gamma = delrho(n,1,1)*delrho(n,1,1) +
164
 
     &              delrho(n,2,1)*delrho(n,2,1) +
165
 
     &              delrho(n,3,1)*delrho(n,3,1)
166
 
            if (dsqrt(gamma).gt.tol_rho)then
167
 
               x = sqrt(gamma) / rho43
168
 
            else
169
 
               x = 0d0
170
 
            endif
171
 
c
172
 
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
173
 
            g = -BETA*x*x / gdenom
174
 
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
175
 
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
176
 
c
177
 
            Etmp = 0.d0
178
 
            Atmp = 0.d0
179
 
            Ctmp = 0.d0
180
 
            if (lfac) then
181
 
               Etmp = rho43*C*fac
182
 
               Atmp = (4d0/3d0)*rho13*C*fac
183
 
            endif
184
 
c
185
 
            if (nlfac) then
186
 
               Etmp = Etmp + rho43*g*fac
187
 
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
188
 
            endif
189
 
c
190
 
            if (x.gt.tol_rho) then
191
 
               Ctmp = 0.5d0*dg / sqrt(gamma) * fac
192
 
            endif
193
 
c
194
 
#ifdef SECOND_DERIV
195
 
            if (lfac) g = g + C           ! Add local contribution back to g
196
 
            rhom23 = rho13 / rho(n,2)
197
 
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
198
 
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
199
 
     &           + BETA*x*x*d2gdenom/gdenom**2
200
 
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
201
 
c
202
 
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
203
 
            C2tmp =  (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac
204
 
            if (x.gt.tol_rho) then
205
 
               C3tmp = -0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
206
 
            endif
207
 
#endif
208
 
#ifdef SECOND_DERIV
209
 
            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
210
 
     &           C3tmp)
211
 
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
212
 
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
213
 
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
214
 
#else
215
 
           call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
216
 
#endif
217
 
            Ex = Ex + qwght(n)*Etmp
218
 
            if(ldew) func(n) = func(n) + Etmp
219
 
            Amat(n,1) = Amat(n,1) + Atmp
220
 
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
221
 
c
222
 
 25         continue
223
 
c
224
 
c           Spin beta:
225
 
c
226
 
            if (dabs(rho(n,3)).lt.tol_rho) goto 20
227
 
c
228
 
            rho13 = abs(rho(n,3))**(1.d0/3.d0)*sign(1d0,rho(n,3))
229
 
            rho43 = rho13**4 
230
 
            gamma = delrho(n,1,2)*delrho(n,1,2) +
231
 
     &              delrho(n,2,2)*delrho(n,2,2) +
232
 
     &              delrho(n,3,2)*delrho(n,3,2)
233
 
            if (dsqrt(gamma).gt.tol_rho)then
234
 
               x = sqrt(gamma) / rho43
235
 
            else
236
 
               x = 0d0
237
 
            endif
238
 
c
239
 
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
240
 
            g = -BETA*x*x / gdenom
241
 
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
242
 
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
243
 
c
244
 
            Etmp = 0.d0
245
 
            Atmp = 0.d0
246
 
            Ctmp = 0.d0
247
 
            if (lfac) then
248
 
               Etmp = rho43*C*fac
249
 
               Atmp = (4d0/3d0)*rho13*C*fac
250
 
            endif
251
 
c
252
 
            if (nlfac) then
253
 
               Etmp = Etmp + rho43*g*fac
254
 
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
255
 
            endif
256
 
c
257
 
            if (x.gt.tol_rho) then
258
 
               Ctmp = 0.5d0*dg / sqrt(gamma) * fac
259
 
            endif
260
 
c
261
 
#ifdef SECOND_DERIV
262
 
            A2tmp = 0d0
263
 
            C2tmp = 0d0
264
 
            C3tmp = 0d0
265
 
            if(lfac) g = g + C           ! Add local contribution back to g
266
 
            rhom23 = rho13 / rho(n,3)
267
 
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
268
 
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
269
 
     &           + BETA*x*x*d2gdenom/gdenom**2
270
 
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
271
 
c
272
 
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
273
 
            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac
274
 
            if (x.gt.tol_rho) then
275
 
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
276
 
            endif
277
 
#endif
278
 
c
279
 
#ifdef SECOND_DERIV
280
 
            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
281
 
     &           C3tmp)
282
 
            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
283
 
            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
284
 
            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
285
 
#else
286
 
            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
287
 
#endif
288
 
            Ex = Ex + qwght(n)*Etmp
289
 
            if(ldew) func(n) = func(n) + Etmp
290
 
            Amat(n,2) = Amat(n,2) + Atmp
291
 
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
292
 
 20      continue
293
 
c
294
 
      endif
295
 
c
296
 
      return
297
 
      end
298
 
#ifndef SECOND_DERIV
299
 
#define SECOND_DERIV
300
 
c
301
 
c     Compile source again for the 2nd derivative case
302
 
c
303
 
#include "xc_camb88.F"
304
 
#endif
 
1
#if defined(FUJITSU_VPP)
 
2
!ocl scalar
 
3
#endif
 
4
#ifndef SECOND_DERIV
 
5
      Subroutine xc_camb88(tol_rho, fac, lfac, nlfac, rho, delrho, 
 
6
     &                      Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
 
7
#else
 
8
      Subroutine xc_camb88_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 
 
9
     &                         Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
 
10
     &                         qwght,ldew,func)
 
11
#endif
 
12
c
 
13
C$Id: xc_camb88.F 20247 2011-04-28 18:58:49Z d3y133 $
 
14
c
 
15
c     Coulomb attenuated Becke88 functional
 
16
c
 
17
      implicit none
 
18
c      
 
19
#include "dft2drv.fh"
 
20
c
 
21
      double precision tol_rho, fac, Ex
 
22
      integer nq, ipol
 
23
      logical lfac, nlfac,ldew
 
24
      double precision func(*)  ! value of the functional [output]
 
25
c
 
26
c     Charge Density
 
27
c
 
28
      double precision rho(nq,ipol*(ipol+1)/2)
 
29
c
 
30
c     Charge Density Gradient
 
31
c
 
32
      double precision delrho(nq,3,ipol)
 
33
c
 
34
c     Quadrature Weights
 
35
c
 
36
      double precision qwght(nq)
 
37
c
 
38
c     Sampling Matrices for the XC Potential
 
39
c
 
40
      double precision Amat(nq,ipol), Cmat(nq,*)
 
41
c
 
42
#ifdef SECOND_DERIV
 
43
c
 
44
c     Second Derivatives of the Exchange Energy Functional
 
45
c
 
46
      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
 
47
      double precision A2tmp, C2tmp, C3tmp
 
48
#endif
 
49
c
 
50
      double precision BETA
 
51
      Parameter (BETA = 0.0042D0)
 
52
c
 
53
c References:
 
54
c
 
55
c    Becke, Phys. Rev. A 38, 3098 (1988)
 
56
c    Johnson, Gill & Pople, J. Chem. Phys. 98, 5612 (1993)
 
57
c
 
58
c***************************************************************************
 
59
c
 
60
      integer n
 
61
      double precision arcsinh, darcsinh
 
62
      double precision C, rho13, rho43, gamma, x, g, gdenom, dg,
 
63
     &     dgdenom, t, Etmp, Atmp, Ctmp
 
64
#ifdef SECOND_DERIV
 
65
      double precision rhom23, d2g, d2gdenom
 
66
#endif
 
67
c
 
68
      arcsinh(x)=log(x+dsqrt(1d0+x*x))
 
69
      darcsinh(x)=1d0/dsqrt(1d0+x*x)
 
70
c
 
71
c     Uniform electron gas constant
 
72
c
 
73
      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
 
74
c
 
75
      if (ipol.eq.1) then
 
76
c
 
77
c        ======> SPIN-RESTRICTED <======
 
78
c
 
79
         do 10 n = 1, nq
 
80
            if (rho(n,1).lt.tol_rho) goto 10
 
81
c
 
82
c           Spin alpha:
 
83
c
 
84
            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
 
85
            rho43 = rho13**4 
 
86
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
87
     &              delrho(n,2,1)*delrho(n,2,1) +
 
88
     &              delrho(n,3,1)*delrho(n,3,1)
 
89
            if (dsqrt(gamma).gt.tol_rho)then
 
90
               gamma = 0.25d0 * gamma
 
91
               x = dsqrt(gamma) / rho43
 
92
            else
 
93
               x = 0d0
 
94
            endif
 
95
c
 
96
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
 
97
            g = -BETA*x*x / gdenom
 
98
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
 
99
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
 
100
c
 
101
            Etmp = 0.d0
 
102
            Atmp = 0.d0
 
103
            Ctmp = 0.d0
 
104
            if (lfac) then
 
105
               Etmp = 2d0*rho43*C*fac
 
106
               Atmp = (4d0/3d0)*rho13*C*fac
 
107
            endif
 
108
c
 
109
            if (nlfac) then
 
110
               Etmp = Etmp + 2d0*rho43*g*fac
 
111
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
112
            endif
 
113
c
 
114
            if (x.gt.tol_rho) then
 
115
               Ctmp = 0.5d0 * dg / sqrt(gamma) * fac
 
116
            endif
 
117
c
 
118
#ifdef SECOND_DERIV
 
119
            A2tmp = 0d0
 
120
            C2tmp = 0d0
 
121
            C3tmp = 0d0
 
122
            if(lfac) g = g + C           ! Add local contribution back to g
 
123
            rhom23 = rho13 / (0.5d0*rho(n,1))
 
124
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
 
125
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
 
126
     &           + BETA*x*x*d2gdenom/gdenom**2
 
127
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
 
128
c
 
129
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
130
            C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac
 
131
            if (x.gt.tol_rho) then
 
132
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
133
            endif
 
134
#endif
 
135
c
 
136
#ifdef SECOND_DERIV
 
137
            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
138
     &           C2tmp,C3tmp)
 
139
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
140
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
141
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
142
#else
 
143
            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
 
144
#endif
 
145
            Ex = Ex + qwght(n)*Etmp
 
146
            if(ldew) func(n) = func(n) + Etmp
 
147
            Amat(n,1) = Amat(n,1) + Atmp
 
148
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
149
 10      continue
 
150
c
 
151
      else
 
152
c
 
153
c        ======> SPIN-UNRESTRICTED <======
 
154
c
 
155
         do 20 n = 1, nq
 
156
            if (rho(n,1).lt.tol_rho) goto 20
 
157
            if (rho(n,2).lt.tol_rho) goto 25
 
158
c
 
159
c           Spin alpha:
 
160
c
 
161
            rho13 = rho(n,2)**(1.d0/3.d0)
 
162
            rho43 = rho13*rho(n,2)
 
163
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
164
     &              delrho(n,2,1)*delrho(n,2,1) +
 
165
     &              delrho(n,3,1)*delrho(n,3,1)
 
166
            if (dsqrt(gamma).gt.tol_rho)then
 
167
               x = dsqrt(gamma) / rho43
 
168
            else
 
169
               x = 0d0
 
170
            endif
 
171
c
 
172
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
 
173
            g = -BETA*x*x / gdenom
 
174
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
 
175
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
 
176
c
 
177
            Etmp = 0.d0
 
178
            Atmp = 0.d0
 
179
            Ctmp = 0.d0
 
180
            if (lfac) then
 
181
               Etmp = rho43*C*fac
 
182
               Atmp = (4d0/3d0)*rho13*C*fac
 
183
            endif
 
184
c
 
185
            if (nlfac) then
 
186
               Etmp = Etmp + rho43*g*fac
 
187
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
188
            endif
 
189
c
 
190
            if (x.gt.tol_rho) then
 
191
               Ctmp = 0.5d0*dg / sqrt(gamma) * fac
 
192
            endif
 
193
c
 
194
#ifdef SECOND_DERIV
 
195
            if (lfac) g = g + C           ! Add local contribution back to g
 
196
            rhom23 = rho13 / rho(n,2)
 
197
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
 
198
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
 
199
     &           + BETA*x*x*d2gdenom/gdenom**2
 
200
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
 
201
c
 
202
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
203
            C2tmp =  (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac
 
204
            if (x.gt.tol_rho) then
 
205
               C3tmp = -0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
206
            endif
 
207
#endif
 
208
#ifdef SECOND_DERIV
 
209
            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
 
210
     &           C3tmp)
 
211
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
212
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
213
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
214
#else
 
215
           call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
 
216
#endif
 
217
            Ex = Ex + qwght(n)*Etmp
 
218
            if(ldew) func(n) = func(n) + Etmp
 
219
            Amat(n,1) = Amat(n,1) + Atmp
 
220
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
221
c
 
222
 25         continue
 
223
c
 
224
c           Spin beta:
 
225
c
 
226
            if (rho(n,3).lt.tol_rho) goto 20
 
227
c
 
228
            rho13 = rho(n,3)**(1.d0/3.d0)
 
229
            rho43 = rho13*rho(n,3)
 
230
            gamma = delrho(n,1,2)*delrho(n,1,2) +
 
231
     &              delrho(n,2,2)*delrho(n,2,2) +
 
232
     &              delrho(n,3,2)*delrho(n,3,2)
 
233
            if (dsqrt(gamma).gt.tol_rho)then
 
234
               x = dsqrt(gamma) / rho43
 
235
            else
 
236
               x = 0d0
 
237
            endif
 
238
c
 
239
            gdenom = 1d0 + 6d0*BETA*x*arcsinh(x)
 
240
            g = -BETA*x*x / gdenom
 
241
            dgdenom = 6d0*BETA*(arcsinh(x) + x*darcsinh(x))
 
242
            dg = BETA*x*(x*dgdenom - 2d0*gdenom) / gdenom**2
 
243
c
 
244
            Etmp = 0.d0
 
245
            Atmp = 0.d0
 
246
            Ctmp = 0.d0
 
247
            if (lfac) then
 
248
               Etmp = rho43*C*fac
 
249
               Atmp = (4d0/3d0)*rho13*C*fac
 
250
            endif
 
251
c
 
252
            if (nlfac) then
 
253
               Etmp = Etmp + rho43*g*fac
 
254
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
255
            endif
 
256
c
 
257
            if (x.gt.tol_rho) then
 
258
               Ctmp = 0.5d0*dg / sqrt(gamma) * fac
 
259
            endif
 
260
c
 
261
#ifdef SECOND_DERIV
 
262
            A2tmp = 0d0
 
263
            C2tmp = 0d0
 
264
            C3tmp = 0d0
 
265
            if(lfac) g = g + C           ! Add local contribution back to g
 
266
            rhom23 = rho13 / rho(n,3)
 
267
            d2gdenom = 6d0*BETA*darcsinh(x)*(2d0 - x*x/(x*x+1d0))
 
268
            d2g = -2d0*BETA/gdenom + 4d0*BETA*x*dgdenom/gdenom**2
 
269
     &           + BETA*x*x*d2gdenom/gdenom**2
 
270
     &           - 2d0*BETA*x*x*(dgdenom)**2/gdenom**3
 
271
c
 
272
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
273
            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac
 
274
            if (x.gt.tol_rho) then
 
275
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
276
            endif
 
277
#endif
 
278
c
 
279
#ifdef SECOND_DERIV
 
280
            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,C2tmp,
 
281
     &           C3tmp)
 
282
            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
 
283
            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
 
284
            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
 
285
#else
 
286
            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
 
287
#endif
 
288
            Ex = Ex + qwght(n)*Etmp
 
289
            if(ldew) func(n) = func(n) + Etmp
 
290
            Amat(n,2) = Amat(n,2) + Atmp
 
291
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
 
292
 20      continue
 
293
c
 
294
      endif
 
295
c
 
296
      return
 
297
      end
 
298
#ifndef SECOND_DERIV
 
299
#define SECOND_DERIV
 
300
c
 
301
c     Compile source again for the 2nd derivative case
 
302
c
 
303
#include "xc_camb88.F"
 
304
#endif