← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_camxpbe96.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_camxpbe96.F
1
 
c
2
 
c     Coulomb attenuated PBE exchange functional 
3
 
c
4
 
c     References:
5
 
c     [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996).
6
 
c     [b] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).;
7
 
c                                               40, 3399 (1989) (E).
8
 
c     Hammer, Hansen and Norskov, PRB 59, 7413 (1999) [RPBE]
9
 
c     Zhang and Yang, PRL 80, 890 (1998) [RevPBE]
10
 
c
11
 
#ifndef SECOND_DERIV
12
 
      Subroutine xc_camxpbe96(whichf,
13
 
     W     tol_rho, fac, lfac, nlfac, rho, delrho, 
14
 
     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
15
 
#else
16
 
      Subroutine xc_camxpbe96_d2(whichf,
17
 
     W     tol_rho, fac, lfac, nlfac, rho, delrho, 
18
 
     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
19
 
     &                        qwght,ldew,func)
20
 
#endif
21
 
c
22
 
c$Id: xc_camxpbe96.F,v 1.3 2008/12/09 01:42:28 niri Exp $
23
 
c
24
 
      implicit none
25
 
c
26
 
#include "dft2drv.fh"
27
 
c      
28
 
      character*4 whichf
29
 
      double precision fac, Ex
30
 
      integer nq, ipol
31
 
      logical lfac, nlfac,ldew
32
 
      double precision func(*)  ! value of the functional [output]
33
 
c
34
 
c     Charge Density & Its Cube Root
35
 
c
36
 
      double precision rho(nq,ipol*(ipol+1)/2)
37
 
c
38
 
c     Charge Density Gradient
39
 
c
40
 
      double precision delrho(nq,3,ipol)
41
 
c
42
 
c     Quadrature Weights
43
 
c
44
 
      double precision qwght(nq)
45
 
c
46
 
c     Sampling Matrices for the XC Potential & Energy
47
 
c
48
 
      double precision Amat(nq,ipol), Cmat(nq,*)
49
 
      double precision Atmp, Ctmp, Etmp
50
 
c
51
 
#ifdef SECOND_DERIV
52
 
      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
53
 
      double precision A2tmp, C2tmp, C3tmp
54
 
#endif
55
 
c
56
 
      double precision tol_rho, pi, um, uk, umk,ukrev,umkrev
57
 
      double precision C, Cs
58
 
      double precision F43, F13
59
 
#ifdef SECOND_DERIV
60
 
      double precision F73
61
 
#endif
62
 
      parameter(um=0.2195149727645171d0, uk=0.8040d0, umk=um/uk)
63
 
      parameter(ukrev=1.245d0, umkrev=um/ukrev)
64
 
      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
65
 
#ifdef SECOND_DERIV
66
 
      parameter (F73=7.d0/3.d0)
67
 
#endif
68
 
c
69
 
      integer n
70
 
      double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2),
71
 
     &      d, g, gp, d1g(2)
72
 
#ifdef SECOND_DERIV
73
 
      double precision rhom23, d2s(3), gpp, d2g(3)
74
 
#endif
75
 
      double precision gpbe0,gpbe1,gpbe2
76
 
      double precision grpbe0,grpbe1,grpbe2
77
 
      double precision grevpbe0,grevpbe1,grevpbe2
78
 
      gpbe0(s)= uk*(1d0 - 1d0/(1d0+umk*s*s))
79
 
      gpbe1(s)= 2d0*um*s/(1d0+umk*s*s)**2
80
 
      gpbe2(s)= 2d0*um*(1d0-4d0*umk*s*s/(1d0+umk*s*s))/(1d0+umk*s*s)**2
81
 
      grevpbe0(s)= ukrev*(1d0 - 1d0/(1d0+umkrev*s*s))
82
 
      grevpbe1(s)= 2d0*um*s/(1d0+umkrev*s*s)**2
83
 
      grevpbe2(s)= 2d0*um*(1d0-4d0*umkrev*s*s/(1d0+umkrev*s*s))/
84
 
     /     (1d0+umkrev*s*s)**2
85
 
      grpbe0(s)= uk*(1d0 - exp(-umk*s*s))
86
 
      grpbe1(s)= 2d0*um*s*exp(-umk*s*s)
87
 
      grpbe2(s)= 2d0*um*exp(-umk*s*s)*(1d0-2d0*umk*s*s)
88
 
c
89
 
      pi = acos(-1.d0)
90
 
      C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13
91
 
      Cs = 0.5d0/(3d0*pi*pi)**F13
92
 
      Cs = Cs * C               ! account for including C in rho43
93
 
c
94
 
      if (ipol.eq.1 )then
95
 
c
96
 
c        ======> SPIN-RESTRICTED <======
97
 
c
98
 
#ifdef IFCV81
99
 
CDEC$ NOSWP
100
 
#endif
101
 
         do 10 n = 1, nq
102
 
            if (rho(n,1).lt.tol_rho) goto 10
103
 
            rho43 = C*rho(n,1)**F43
104
 
            rrho = 1d0/rho(n,1)
105
 
            rho13 = F43*rho43*rrho
106
 
#ifdef SECOND_DERIV
107
 
            rhom23 = F13*rho13*rrho
108
 
#endif
109
 
            Etmp = 0.d0
110
 
            Atmp = 0.d0
111
 
            Ctmp = 0.d0
112
 
            if (lfac) then
113
 
               Etmp = rho43*fac
114
 
               Atmp = rho13*fac
115
 
#ifdef SECOND_DERIV
116
 
               A2tmp =  2d0*rhom23*fac
117
 
#endif
118
 
            endif
119
 
c
120
 
            gamma = delrho(n,1,1)*delrho(n,1,1) +
121
 
     &              delrho(n,2,1)*delrho(n,2,1) +
122
 
     &              delrho(n,3,1)*delrho(n,3,1)
123
 
            gam12 = dsqrt(gamma)
124
 
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10
125
 
c
126
 
            s = Cs*gam12/rho43
127
 
            d1s(1) = -F43*s*rrho
128
 
            d1s(2) = 0.5d0*s/gamma
129
 
c
130
 
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
131
 
c
132
 
            if(whichf.eq.'revp') then
133
 
               g=grevpbe0(s)
134
 
               gp=grevpbe1(s)
135
 
            elseif(whichf.eq.'rpbe') then
136
 
               g=grpbe0(s)
137
 
               gp=grpbe1(s)
138
 
            else
139
 
               g=gpbe0(s)
140
 
               gp=gpbe1(s)
141
 
            endif
142
 
c
143
 
            d1g(1) = gp*d1s(1)
144
 
            d1g(2) = gp*d1s(2)
145
 
            Etmp = Etmp + rho43*g*fac
146
 
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
147
 
            Ctmp = 2d0*rho43*d1g(2)*fac
148
 
#ifdef SECOND_DERIV
149
 
            d2s(1) = -F73*d1s(1)*rrho
150
 
            d2s(2) = -F43*d1s(2)*rrho
151
 
            d2s(3) = -0.5d0*d1s(2)/gamma
152
 
            if(whichf.eq.'revp') then
153
 
               gpp=grevpbe2(s)
154
 
            elseif(whichf.eq.'rpbe') then
155
 
               gpp=grpbe2(s)
156
 
            else
157
 
               gpp=gpbe2(s)
158
 
            endif
159
 
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
160
 
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
161
 
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
162
 
            A2tmp = A2tmp
163
 
     &           +(rhom23*g + 2.d0*rho13*d1g(1) + rho43*d2g(1))*fac*2d0
164
 
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac*4d0
165
 
            C3tmp = rho43*d2g(3)*fac*8d0
166
 
c
167
 
            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
168
 
     &           C2tmp,C3tmp)
169
 
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
170
 
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
171
 
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
172
 
#else
173
 
            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
174
 
#endif
175
 
            Ex = Ex + qwght(n)*Etmp
176
 
            if (ldew) func(n) = func(n) + Etmp
177
 
            Amat(n,1) = Amat(n,1) + Atmp
178
 
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
179
 
 10      continue
180
 
c
181
 
      else
182
 
c
183
 
c        ======> SPIN-UNRESTRICTED <======
184
 
c
185
 
#ifdef IFCV81
186
 
CDEC$ NOSWP
187
 
#endif
188
 
         do 20 n = 1, nq
189
 
            if (rho(n,1).lt.tol_rho) goto 20
190
 
c
191
 
c     Alpha
192
 
c
193
 
            if (rho(n,2).lt.tol_rho) goto 25
194
 
            rho43 = C*(2d0*rho(n,2))**F43
195
 
            rrho = 0.5d0/rho(n,2)
196
 
            rho13 = F43*rho43*rrho
197
 
#ifdef SECOND_DERIV
198
 
            rhom23 = F13*rho13*rrho
199
 
#endif
200
 
            Etmp = 0.d0
201
 
            Atmp = 0.d0
202
 
            Ctmp = 0.d0
203
 
            if (lfac) then
204
 
               Etmp = rho43*fac*0.5d0
205
 
               Atmp = rho13*fac
206
 
#ifdef SECOND_DERIV
207
 
               A2tmp = 2d0*rhom23*fac
208
 
#endif
209
 
            endif
210
 
c
211
 
            gamma = delrho(n,1,1)*delrho(n,1,1) +
212
 
     &              delrho(n,2,1)*delrho(n,2,1) +
213
 
     &              delrho(n,3,1)*delrho(n,3,1)
214
 
            gam12 = 2d0*dsqrt(gamma)
215
 
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25
216
 
c
217
 
            s = Cs*gam12/rho43
218
 
            d1s(1) = -F43*s*rrho
219
 
            d1s(2) = 0.5d0*s/gamma
220
 
c
221
 
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
222
 
c
223
 
            if(whichf.eq.'revp') then
224
 
               g=grevpbe0(s)
225
 
               gp=grevpbe1(s)
226
 
            elseif(whichf.eq.'rpbe') then
227
 
               g=grpbe0(s)
228
 
               gp=grpbe1(s)
229
 
            else
230
 
               g=gpbe0(s)
231
 
               gp=gpbe1(s)
232
 
            endif
233
 
c
234
 
            d1g(1) = gp*d1s(1)
235
 
            d1g(2) = gp*d1s(2)
236
 
            Etmp = Etmp + rho43*g*fac*0.5d0
237
 
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
238
 
            Ctmp = 0.5d0*rho43*d1g(2)*fac
239
 
#ifdef SECOND_DERIV
240
 
            d2s(1) = -F73*d1s(1)*rrho
241
 
            d2s(2) = -F43*d1s(2)*rrho
242
 
            d2s(3) = -0.5d0*d1s(2)/gamma
243
 
            if(whichf.eq.'revp') then
244
 
               gpp=grevpbe2(s)
245
 
            elseif(whichf.eq.'rpbe') then
246
 
               gpp=grpbe2(s)
247
 
            else
248
 
               gpp=gpbe2(s)
249
 
            endif
250
 
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
251
 
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
252
 
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
253
 
            A2tmp = A2tmp + (rhom23*g + 2.d0*rho13*d1g(1)
254
 
     &           + rho43*d2g(1))*fac*2d0
255
 
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac
256
 
            C3tmp = rho43*d2g(3)*fac*0.5d0
257
 
c
258
 
            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
259
 
     &           C2tmp,C3tmp)
260
 
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
261
 
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
262
 
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
263
 
#else
264
 
            call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
265
 
#endif
266
 
            Ex = Ex + qwght(n)*Etmp
267
 
            if (ldew) func(n) = func(n) + Etmp
268
 
            Amat(n,1) = Amat(n,1) + Atmp
269
 
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
270
 
c
271
 
c     Beta
272
 
c
273
 
 25         continue
274
 
            if (rho(n,3).lt.tol_rho) goto 20
275
 
            rho43 = C*(2d0*rho(n,3))**F43
276
 
            rrho = 0.5d0/rho(n,3)
277
 
            rho13 = F43*rho43*rrho
278
 
#ifdef SECOND_DERIV
279
 
            rhom23 = F13*rho13*rrho
280
 
#endif
281
 
            Etmp = 0.d0
282
 
            Atmp = 0.d0
283
 
            Ctmp = 0.d0
284
 
            if (lfac) then
285
 
               Etmp = rho43*fac*0.5d0
286
 
               Atmp = rho13*fac
287
 
#ifdef SECOND_DERIV
288
 
               A2tmp= 2.d0*rhom23*fac
289
 
#endif
290
 
            endif
291
 
c
292
 
            gamma = delrho(n,1,2)*delrho(n,1,2) +
293
 
     &              delrho(n,2,2)*delrho(n,2,2) +
294
 
     &              delrho(n,3,2)*delrho(n,3,2)
295
 
            gam12 = 2d0*dsqrt(gamma)
296
 
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20
297
 
c
298
 
            s = Cs*gam12/rho43
299
 
            d1s(1) = -F43*s*rrho
300
 
            d1s(2) = 0.5d0*s/gamma
301
 
c
302
 
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
303
 
c
304
 
            if(whichf.eq.'revp') then
305
 
               g=grevpbe0(s)
306
 
               gp=grevpbe1(s)
307
 
            elseif(whichf.eq.'rpbe') then
308
 
               g=grpbe0(s)
309
 
               gp=grpbe1(s)
310
 
            else
311
 
               g=gpbe0(s)
312
 
               gp=gpbe1(s)
313
 
            endif
314
 
c
315
 
            d1g(1) = gp*d1s(1)
316
 
            d1g(2) = gp*d1s(2)
317
 
            Etmp = Etmp + rho43*g*fac*0.5d0
318
 
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
319
 
            Ctmp = 0.5d0*rho43*d1g(2)*fac
320
 
#ifdef SECOND_DERIV
321
 
            d2s(1) = -F73*d1s(1)*rrho
322
 
            d2s(2) = -F43*d1s(2)*rrho
323
 
            d2s(3) = -0.5d0*d1s(2)/gamma
324
 
            if(whichf.eq.'revp') then
325
 
               gpp=grevpbe2(s)
326
 
            elseif(whichf.eq.'rpbe') then
327
 
               gpp=grpbe2(s)
328
 
            else
329
 
               gpp=gpbe2(s)
330
 
            endif
331
 
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
332
 
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
333
 
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
334
 
c
335
 
            A2tmp = A2tmp + (rhom23*g + 2.d0*rho13*d1g(1)
336
 
     &           + rho43*d2g(1))*fac*2d0
337
 
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac
338
 
            C3tmp =  rho43*d2g(3)*fac*0.5d0
339
 
c
340
 
            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
341
 
     &           C2tmp,C3tmp)
342
 
            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
343
 
            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
344
 
            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
345
 
#else
346
 
            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
347
 
#endif
348
 
            Ex = Ex + qwght(n)*Etmp
349
 
            if (ldew) func(n) = func(n) + Etmp
350
 
            Amat(n,2) = Amat(n,2) + Atmp
351
 
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
352
 
c
353
 
 20      continue
354
 
      endif
355
 
c
356
 
      return
357
 
      end
358
 
#ifndef SECOND_DERIV
359
 
#define SECOND_DERIV
360
 
c
361
 
c     Compile source again for the 2nd derivative case
362
 
c
363
 
#include "xc_camxpbe96.F"
364
 
#endif
 
1
c
 
2
c     Coulomb attenuated PBE exchange functional 
 
3
c
 
4
c     References:
 
5
c     [a] J.P. Perdew, K. Burke, and M. Ernzerhof, PRL 77, 3865 (1996).
 
6
c     [b] J.P. Perdew and Y. Wang, Phys. Rev. B 33, 8800 (1986).;
 
7
c                                               40, 3399 (1989) (E).
 
8
c     Hammer, Hansen and Norskov, PRB 59, 7413 (1999) [RPBE]
 
9
c     Zhang and Yang, PRL 80, 890 (1998) [RevPBE]
 
10
c
 
11
#ifndef SECOND_DERIV
 
12
      Subroutine xc_camxpbe96(whichf,
 
13
     W     tol_rho, fac, lfac, nlfac, rho, delrho, 
 
14
     &                     Amat, Cmat, nq, ipol, Ex, qwght,ldew,func)
 
15
#else
 
16
      Subroutine xc_camxpbe96_d2(whichf,
 
17
     W     tol_rho, fac, lfac, nlfac, rho, delrho, 
 
18
     &                        Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
 
19
     &                        qwght,ldew,func)
 
20
#endif
 
21
c
 
22
c$Id: xc_camxpbe96.F 20247 2011-04-28 18:58:49Z d3y133 $
 
23
c
 
24
      implicit none
 
25
c
 
26
#include "dft2drv.fh"
 
27
c      
 
28
      character*4 whichf
 
29
      double precision fac, Ex
 
30
      integer nq, ipol
 
31
      logical lfac, nlfac,ldew
 
32
      double precision func(*)  ! value of the functional [output]
 
33
c
 
34
c     Charge Density & Its Cube Root
 
35
c
 
36
      double precision rho(nq,ipol*(ipol+1)/2)
 
37
c
 
38
c     Charge Density Gradient
 
39
c
 
40
      double precision delrho(nq,3,ipol)
 
41
c
 
42
c     Quadrature Weights
 
43
c
 
44
      double precision qwght(nq)
 
45
c
 
46
c     Sampling Matrices for the XC Potential & Energy
 
47
c
 
48
      double precision Amat(nq,ipol), Cmat(nq,*)
 
49
      double precision Atmp, Ctmp, Etmp
 
50
c
 
51
#ifdef SECOND_DERIV
 
52
      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
 
53
      double precision A2tmp, C2tmp, C3tmp
 
54
#endif
 
55
c
 
56
      double precision tol_rho, pi, um, uk, umk,ukrev,umkrev
 
57
      double precision C, Cs
 
58
      double precision F43, F13
 
59
#ifdef SECOND_DERIV
 
60
      double precision F73
 
61
#endif
 
62
      parameter(um=0.2195149727645171d0, uk=0.8040d0, umk=um/uk)
 
63
      parameter(ukrev=1.245d0, umkrev=um/ukrev)
 
64
      parameter (F43=4.d0/3.d0, F13=1.d0/3.d0)
 
65
#ifdef SECOND_DERIV
 
66
      parameter (F73=7.d0/3.d0)
 
67
#endif
 
68
c
 
69
      integer n
 
70
      double precision rrho, rho43, rho13, gamma, gam12, s, d1s(2),
 
71
     &      d, g, gp, d1g(2)
 
72
#ifdef SECOND_DERIV
 
73
      double precision rhom23, d2s(3), gpp, d2g(3)
 
74
#endif
 
75
      double precision gpbe0,gpbe1,gpbe2
 
76
      double precision grpbe0,grpbe1,grpbe2
 
77
      double precision grevpbe0,grevpbe1,grevpbe2
 
78
      gpbe0(s)= uk*(1d0 - 1d0/(1d0+umk*s*s))
 
79
      gpbe1(s)= 2d0*um*s/(1d0+umk*s*s)**2
 
80
      gpbe2(s)= 2d0*um*(1d0-4d0*umk*s*s/(1d0+umk*s*s))/(1d0+umk*s*s)**2
 
81
      grevpbe0(s)= ukrev*(1d0 - 1d0/(1d0+umkrev*s*s))
 
82
      grevpbe1(s)= 2d0*um*s/(1d0+umkrev*s*s)**2
 
83
      grevpbe2(s)= 2d0*um*(1d0-4d0*umkrev*s*s/(1d0+umkrev*s*s))/
 
84
     /     (1d0+umkrev*s*s)**2
 
85
      grpbe0(s)= uk*(1d0 - exp(-umk*s*s))
 
86
      grpbe1(s)= 2d0*um*s*exp(-umk*s*s)
 
87
      grpbe2(s)= 2d0*um*exp(-umk*s*s)*(1d0-2d0*umk*s*s)
 
88
c
 
89
      pi = acos(-1.d0)
 
90
      C = -3d0/(4d0*pi)*(3d0*pi*pi)**F13
 
91
      Cs = 0.5d0/(3d0*pi*pi)**F13
 
92
      Cs = Cs * C               ! account for including C in rho43
 
93
c
 
94
      if (ipol.eq.1 )then
 
95
c
 
96
c        ======> SPIN-RESTRICTED <======
 
97
c
 
98
#ifdef IFCV81
 
99
CDEC$ NOSWP
 
100
#endif
 
101
         do 10 n = 1, nq
 
102
            if (rho(n,1).lt.tol_rho) goto 10
 
103
            rho43 = C*rho(n,1)**F43
 
104
            rrho = 1d0/rho(n,1)
 
105
            rho13 = F43*rho43*rrho
 
106
#ifdef SECOND_DERIV
 
107
            rhom23 = F13*rho13*rrho
 
108
#endif
 
109
            Etmp = 0.d0
 
110
            Atmp = 0.d0
 
111
            Ctmp = 0.d0
 
112
            if (lfac) then
 
113
               Etmp = rho43*fac
 
114
               Atmp = rho13*fac
 
115
#ifdef SECOND_DERIV
 
116
               A2tmp =  2d0*rhom23*fac
 
117
#endif
 
118
            endif
 
119
c
 
120
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
121
     &              delrho(n,2,1)*delrho(n,2,1) +
 
122
     &              delrho(n,3,1)*delrho(n,3,1)
 
123
            gam12 = dsqrt(gamma)
 
124
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 10
 
125
c
 
126
            s = Cs*gam12/rho43
 
127
            d1s(1) = -F43*s*rrho
 
128
            d1s(2) = 0.5d0*s/gamma
 
129
c
 
130
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
 
131
c
 
132
            if(whichf.eq.'revp') then
 
133
               g=grevpbe0(s)
 
134
               gp=grevpbe1(s)
 
135
            elseif(whichf.eq.'rpbe') then
 
136
               g=grpbe0(s)
 
137
               gp=grpbe1(s)
 
138
            else
 
139
               g=gpbe0(s)
 
140
               gp=gpbe1(s)
 
141
            endif
 
142
c
 
143
            d1g(1) = gp*d1s(1)
 
144
            d1g(2) = gp*d1s(2)
 
145
            Etmp = Etmp + rho43*g*fac
 
146
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
 
147
            Ctmp = 2d0*rho43*d1g(2)*fac
 
148
#ifdef SECOND_DERIV
 
149
            d2s(1) = -F73*d1s(1)*rrho
 
150
            d2s(2) = -F43*d1s(2)*rrho
 
151
            d2s(3) = -0.5d0*d1s(2)/gamma
 
152
            if(whichf.eq.'revp') then
 
153
               gpp=grevpbe2(s)
 
154
            elseif(whichf.eq.'rpbe') then
 
155
               gpp=grpbe2(s)
 
156
            else
 
157
               gpp=gpbe2(s)
 
158
            endif
 
159
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
 
160
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
 
161
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
 
162
            A2tmp = A2tmp
 
163
     &           +(rhom23*g + 2.d0*rho13*d1g(1) + rho43*d2g(1))*fac*2d0
 
164
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac*4d0
 
165
            C3tmp = rho43*d2g(3)*fac*8d0
 
166
c
 
167
            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
168
     &           C2tmp,C3tmp)
 
169
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
170
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
171
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
172
#else
 
173
            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
 
174
#endif
 
175
            Ex = Ex + qwght(n)*Etmp
 
176
            if (ldew) func(n) = func(n) + Etmp
 
177
            Amat(n,1) = Amat(n,1) + Atmp
 
178
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
179
 10      continue
 
180
c
 
181
      else
 
182
c
 
183
c        ======> SPIN-UNRESTRICTED <======
 
184
c
 
185
#ifdef IFCV81
 
186
CDEC$ NOSWP
 
187
#endif
 
188
         do 20 n = 1, nq
 
189
            if (rho(n,1).lt.tol_rho) goto 20
 
190
c
 
191
c     Alpha
 
192
c
 
193
            if (rho(n,2).lt.tol_rho) goto 25
 
194
            rho43 = C*(2d0*rho(n,2))**F43
 
195
            rrho = 0.5d0/rho(n,2)
 
196
            rho13 = F43*rho43*rrho
 
197
#ifdef SECOND_DERIV
 
198
            rhom23 = F13*rho13*rrho
 
199
#endif
 
200
            Etmp = 0.d0
 
201
            Atmp = 0.d0
 
202
            Ctmp = 0.d0
 
203
            if (lfac) then
 
204
               Etmp = rho43*fac*0.5d0
 
205
               Atmp = rho13*fac
 
206
#ifdef SECOND_DERIV
 
207
               A2tmp = 2d0*rhom23*fac
 
208
#endif
 
209
            endif
 
210
c
 
211
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
212
     &              delrho(n,2,1)*delrho(n,2,1) +
 
213
     &              delrho(n,3,1)*delrho(n,3,1)
 
214
            gam12 = 2d0*dsqrt(gamma)
 
215
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 25
 
216
c
 
217
            s = Cs*gam12/rho43
 
218
            d1s(1) = -F43*s*rrho
 
219
            d1s(2) = 0.5d0*s/gamma
 
220
c
 
221
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
 
222
c
 
223
            if(whichf.eq.'revp') then
 
224
               g=grevpbe0(s)
 
225
               gp=grevpbe1(s)
 
226
            elseif(whichf.eq.'rpbe') then
 
227
               g=grpbe0(s)
 
228
               gp=grpbe1(s)
 
229
            else
 
230
               g=gpbe0(s)
 
231
               gp=gpbe1(s)
 
232
            endif
 
233
c
 
234
            d1g(1) = gp*d1s(1)
 
235
            d1g(2) = gp*d1s(2)
 
236
            Etmp = Etmp + rho43*g*fac*0.5d0
 
237
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
 
238
            Ctmp = 0.5d0*rho43*d1g(2)*fac
 
239
#ifdef SECOND_DERIV
 
240
            d2s(1) = -F73*d1s(1)*rrho
 
241
            d2s(2) = -F43*d1s(2)*rrho
 
242
            d2s(3) = -0.5d0*d1s(2)/gamma
 
243
            if(whichf.eq.'revp') then
 
244
               gpp=grevpbe2(s)
 
245
            elseif(whichf.eq.'rpbe') then
 
246
               gpp=grpbe2(s)
 
247
            else
 
248
               gpp=gpbe2(s)
 
249
            endif
 
250
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
 
251
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
 
252
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
 
253
            A2tmp = A2tmp + (rhom23*g + 2.d0*rho13*d1g(1)
 
254
     &           + rho43*d2g(1))*fac*2d0
 
255
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac
 
256
            C3tmp = rho43*d2g(3)*fac*0.5d0
 
257
c
 
258
            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
259
     &           C2tmp,C3tmp)
 
260
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
261
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
262
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
263
#else
 
264
            call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
 
265
#endif
 
266
            Ex = Ex + qwght(n)*Etmp
 
267
            if (ldew) func(n) = func(n) + Etmp
 
268
            Amat(n,1) = Amat(n,1) + Atmp
 
269
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
270
c
 
271
c     Beta
 
272
c
 
273
 25         continue
 
274
            if (rho(n,3).lt.tol_rho) goto 20
 
275
            rho43 = C*(2d0*rho(n,3))**F43
 
276
            rrho = 0.5d0/rho(n,3)
 
277
            rho13 = F43*rho43*rrho
 
278
#ifdef SECOND_DERIV
 
279
            rhom23 = F13*rho13*rrho
 
280
#endif
 
281
            Etmp = 0.d0
 
282
            Atmp = 0.d0
 
283
            Ctmp = 0.d0
 
284
            if (lfac) then
 
285
               Etmp = rho43*fac*0.5d0
 
286
               Atmp = rho13*fac
 
287
#ifdef SECOND_DERIV
 
288
               A2tmp= 2.d0*rhom23*fac
 
289
#endif
 
290
            endif
 
291
c
 
292
            gamma = delrho(n,1,2)*delrho(n,1,2) +
 
293
     &              delrho(n,2,2)*delrho(n,2,2) +
 
294
     &              delrho(n,3,2)*delrho(n,3,2)
 
295
            gam12 = 2d0*dsqrt(gamma)
 
296
            if (.not.(nlfac.and.gam12.gt.tol_rho**2)) goto 20
 
297
c
 
298
            s = Cs*gam12/rho43
 
299
            d1s(1) = -F43*s*rrho
 
300
            d1s(2) = 0.5d0*s/gamma
 
301
c
 
302
c     Evaluate the GC part of F(s), i.e. g(s) = F(s) - 1
 
303
c
 
304
            if(whichf.eq.'revp') then
 
305
               g=grevpbe0(s)
 
306
               gp=grevpbe1(s)
 
307
            elseif(whichf.eq.'rpbe') then
 
308
               g=grpbe0(s)
 
309
               gp=grpbe1(s)
 
310
            else
 
311
               g=gpbe0(s)
 
312
               gp=gpbe1(s)
 
313
            endif
 
314
c
 
315
            d1g(1) = gp*d1s(1)
 
316
            d1g(2) = gp*d1s(2)
 
317
            Etmp = Etmp + rho43*g*fac*0.5d0
 
318
            Atmp = Atmp + (rho13*g+rho43*d1g(1))*fac
 
319
            Ctmp = 0.5d0*rho43*d1g(2)*fac
 
320
#ifdef SECOND_DERIV
 
321
            d2s(1) = -F73*d1s(1)*rrho
 
322
            d2s(2) = -F43*d1s(2)*rrho
 
323
            d2s(3) = -0.5d0*d1s(2)/gamma
 
324
            if(whichf.eq.'revp') then
 
325
               gpp=grevpbe2(s)
 
326
            elseif(whichf.eq.'rpbe') then
 
327
               gpp=grpbe2(s)
 
328
            else
 
329
               gpp=gpbe2(s)
 
330
            endif
 
331
            d2g(1) = gp*d2s(1) + gpp*d1s(1)*d1s(1)
 
332
            d2g(2) = gp*d2s(2) + gpp*d1s(1)*d1s(2)
 
333
            d2g(3) = gp*d2s(3) + gpp*d1s(2)*d1s(2)
 
334
c
 
335
            A2tmp = A2tmp + (rhom23*g + 2.d0*rho13*d1g(1)
 
336
     &           + rho43*d2g(1))*fac*2d0
 
337
            C2tmp = (rho13*d1g(2) + rho43*d2g(2))*fac
 
338
            C3tmp =  rho43*d2g(3)*fac*0.5d0
 
339
c
 
340
            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
341
     &           C2tmp,C3tmp)
 
342
            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
 
343
            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
 
344
            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
 
345
#else
 
346
            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
 
347
#endif
 
348
            Ex = Ex + qwght(n)*Etmp
 
349
            if (ldew) func(n) = func(n) + Etmp
 
350
            Amat(n,2) = Amat(n,2) + Atmp
 
351
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
 
352
c
 
353
 20      continue
 
354
      endif
 
355
c
 
356
      return
 
357
      end
 
358
#ifndef SECOND_DERIV
 
359
#define SECOND_DERIV
 
360
c
 
361
c     Compile source again for the 2nd derivative case
 
362
c
 
363
#include "xc_camxpbe96.F"
 
364
#endif