← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwxc/nwxc_x_m05.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/nwxc/nwxc_x_m05.F
 
1
C> \ingroup nwxc
 
2
C> @{
 
3
C>
 
4
C> \file nwxc_c_m05.F
 
5
C> Implementation of the M05 correlation functional
 
6
C>
 
7
C> @}
 
8
C>
 
9
C> \ingroup nwxc_priv
 
10
C> @{
 
11
C>
 
12
C> \brief The M05 correlation functional
 
13
C>
 
14
C> The M05 functional [1,2] is a meta-GGA of which this evaluates
 
15
C> the correlation component.
 
16
C>
 
17
C> ### References ###
 
18
C>
 
19
C> [1] Y Zhao, NE Schultz, DG Truhlar,
 
20
C>     "Exchange-correlation functional with broad accuracy for 
 
21
C>     metallic and nonmetallic compounds, kinetics, and 
 
22
C>     noncovalent interactions",
 
23
C>     J.Chem.Phys. <b>123</b>, 161103-161106 (2005), DOI:
 
24
C>     <a href="http://dx.doi.org/10.1063/1.2126975">
 
25
C>     10.1063/1.2126975</a>.
 
26
C>
 
27
C> [2] Y Zhao, NE Schultz, DG Truhlar,
 
28
C>     "Design of density functionals by combining the method of
 
29
C>     constraint satisfaction parametrization for thermochemistry,
 
30
C>     thermochemical kinetics, and noncovalent interactions",
 
31
C>     J.Chem.Theory.Comput. <b>2</b>, 364-382 (2006), DOI:
 
32
C>     <a href="http://dx.doi.org/10.1021/ct0502763">
 
33
C>     10.1021/ct0502763</a>.
 
34
C>
 
35
 
 
36
c   M05 or M05-2X  exchange functional  
 
37
c           META GGA
 
38
C         utilizes ingredients:
 
39
c                              rho   -  density
 
40
c                              delrho - gradient of density
 
41
c                              tau - K.S kinetic energy density
 
42
c                              tauU - uniform-gas KE density
 
43
c                              ijzy - 1  M05
 
44
c                              ijzy - 2  M05-2X  
 
45
c     References: 
 
46
c     [a]       Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Phys. 2005, 123, 161103;
 
47
c      Note that in this communication we interchanged cCab,i and cCss,i in Table 1. 
 
48
c     [b]       Zhao, Y.; Schultz, N. E.; Truhlar, D. G. J. Chem. Theory Comput. 2006, in press.    
 
49
 
 
50
 
 
51
      Subroutine nwxc_x_m05(param, tol_rho, ipol, nq, wght, rho, rgamma,
 
52
     &                      tau, func, Amat, Cmat, Mmat)
 
53
c   
 
54
c$Id: xc_xm05.F 22535 2012-05-31 01:22:18Z edo $
 
55
c
 
56
      implicit none
 
57
c
 
58
#include "nwxc_param.fh"
 
59
c
 
60
c
 
61
c     Input and other parameters
 
62
c
 
63
      double precision param(*) !< [Input] Parameters of functional
 
64
                                !< - param(1): \f$ a_{1} \f$
 
65
                                !< - param(2): \f$ a_{2} \f$
 
66
                                !< - param(3): \f$ a_{3} \f$
 
67
                                !< - param(4): \f$ a_{4} \f$
 
68
                                !< - param(5): \f$ a_{5} \f$
 
69
                                !< - param(6): \f$ a_{6} \f$
 
70
                                !< - param(7): \f$ a_{7} \f$
 
71
                                !< - param(8): \f$ a_{8} \f$
 
72
                                !< - param(9): \f$ a_{9} \f$
 
73
                                !< - param(10): \f$ a_{10} \f$
 
74
                                !< - param(11): \f$ a_{11} \f$
 
75
      double precision tol_rho !< [Input] The lower limit on the density
 
76
      integer nq               !< [Input] The number of points
 
77
      integer ipol             !< [Input] The number of spin channels
 
78
      double precision wght    !< [Input] The weight of the functional
 
79
c
 
80
c     Charge Density 
 
81
c
 
82
      double precision rho(nq,*) !< [Input] The density
 
83
c
 
84
c     Charge Density Gradient Norm
 
85
c
 
86
      double precision rgamma(nq,*) !< [Input] The density gradient norm
 
87
c
 
88
c     Kinetic Energy Density
 
89
c
 
90
      double precision tau(nq,*) !< [Input] The kinetic energy density
 
91
c
 
92
c     Functional values
 
93
c
 
94
      double precision func(*) !< [Output] The functional value
 
95
c
 
96
c     Sampling Matrices for the XC Potential
 
97
c
 
98
      double precision Amat(nq,*) !< [Output] Derivative wrt density
 
99
      double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma
 
100
      double precision Mmat(nq,*) !< [Output] Derivative wrt tau
 
101
c      
 
102
      double precision pi
 
103
c
 
104
      integer n
 
105
      double precision at1, at2, at3, at4, at5, at6, at7, at8, at9
 
106
      double precision at, at10, at11, C1, C2, fL, fNL
 
107
      double precision rrho, rho43, rho13, rhoo, rho53
 
108
      double precision Gamma2, Gamma
 
109
      double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6
 
110
      double precision W7, W8, W9, W10, W11, Fsig
 
111
c
 
112
c     kinetic energy density   or  tau
 
113
c
 
114
      double precision tauN,tauu,DTol
 
115
      double precision F83, F23, F53, F43, F13, F1o3
 
116
      double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5
 
117
      double precision One, Two, Three, Four, Five, Six, Seven, Eight
 
118
      double precision Nine, F10, F11
 
119
      double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd
 
120
 
 
121
c      functional derivatives below FFFFFFFFFFFF
 
122
 
 
123
       double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR
 
124
       double precision dFdTau, dGGAdG,tauW
 
125
 
 
126
c      functional derivatives above FFFFFFFFFFFF
 
127
 
 
128
 
 
129
       parameter( pi = 3.1415926535897932384626433832795d0 )
 
130
         
 
131
        
 
132
        parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 
 
133
     &             F3o2=3.d0/2.d0,F13=1.d0/3.d0)
 
134
        parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0)
 
135
        parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0)
 
136
        parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 
 
137
     &             Five=5.0d0,Six=6.0d0, Seven=7.0d0,
 
138
     &             Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0)
 
139
 
 
140
 
 
141
      at1=  param(1)
 
142
      at2=  param(2)
 
143
      at3=  param(3)
 
144
      at4=  param(4)
 
145
      at5=  param(5)
 
146
      at6=  param(6)
 
147
      at7=  param(7)
 
148
      at8=  param(8)
 
149
      at9=  param(9)
 
150
      at10= param(10)
 
151
      at11= param(11)
 
152
c     if (ijzy.eq.1) then
 
153
c     Parameters for M05 
 
154
c       at1=    0.08151d0
 
155
c       at2=    -0.43956d0
 
156
c       at3=    -3.22422d0
 
157
c       at4=    2.01819d0
 
158
c       at5=    8.79431d0
 
159
c       at6=    -0.00295d0
 
160
c       at7=    9.82029d0
 
161
c       at8=    -4.82351d0
 
162
c       at9=    -48.17574d0
 
163
c       at10=   3.64802d0
 
164
c       at11=   34.02248d0
 
165
c     elseif (ijzy.eq.2) then
 
166
c     Parameters for M05-2X
 
167
c       at1=    -0.56833d0
 
168
c       at2=    -1.30057d0
 
169
c       at3=    5.50070d0
 
170
c       at4=    9.06402d0
 
171
c       at5=    -32.21075d0
 
172
c       at6=    -23.73298d0
 
173
c       at7=    70.22996d0
 
174
c       at8=    29.88614d0
 
175
c       at9=    -60.25778d0
 
176
c       at10=   -13.22205d0
 
177
c       at11=   15.23694d0
 
178
c     endif
 
179
      
 
180
      at=1.0d0
 
181
      C1 = 0.2195149727645171d0
 
182
      C2 = C1/0.804d0 
 
183
      DTol=tol_rho
 
184
C
 
185
C     Scale factors for local and non-local contributions.
 
186
C
 
187
      fL  =  wght
 
188
      fNL =  wght
 
189
      Cs = 0.5d0/(3.0d0*pi*pi)**F13
 
190
      P32 = (3.d0*pi**2)**F23
 
191
         
 
192
c     
 
193
       Ax = (-0.75d0)*(3.0d0/pi)**F13
 
194
 
 
195
 
 
196
c
 
197
      if (ipol.eq.1 )then
 
198
c
 
199
c        ======> SPIN-RESTRICTED <======
 
200
c                     or
 
201
c                SPIN-UNPOLARIZED
 
202
c
 
203
c
 
204
         do 10 n = 1, nq
 
205
           
 
206
            if (rho(n,R_T).lt.DTol) goto 10
 
207
            rhoo = rho(n,R_T)
 
208
            rho43 = rhoo**F4o3  
 
209
            rrho = 1d0/rhoo       ! reciprocal of rho
 
210
            rho13 = rho43*rrho
 
211
            rho53 = rhoo**F53
 
212
 
 
213
c
 
214
             
 
215
            tauN = tau(n,T_T)
 
216
            tauu=tauN 
 
217
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
218
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
219
c    &              delrho(n,3,1)*delrho(n,3,1)
 
220
            Gamma2 = rgamma(n,G_TT)
 
221
            Gamma = dsqrt(Gamma2)
 
222
            if (gamma.lt.DTol) goto 10
 
223
            TauUEG=0.3d0*P32*rho53
 
224
            Tsig =TauUEG/tauu
 
225
            Wsig =(Tsig-One)/(Tsig+One)
 
226
            W1=Wsig 
 
227
            W2=Wsig*W1
 
228
            W3=Wsig*W2
 
229
            W4=Wsig*W3
 
230
            W5=Wsig*W4
 
231
            W6=Wsig*W5
 
232
            W7=Wsig*W6
 
233
            W8=Wsig*W7
 
234
            W9=Wsig*W8
 
235
            W10=Wsig*W9
 
236
            W11=Wsig*W10
 
237
            Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
238
     &          + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
239
     &          + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
240
 
 
241
            s      = Cs*Gamma/rho43
 
242
            s2     = s*s
 
243
            En     = C1*s2
 
244
            Ed     = One + C2*s2
 
245
            E      = En/Ed
 
246
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n)
 
247
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)
 
248
c
 
249
c     functional derivatives 
 
250
c
 
251
            dEn   = Two*C1*s
 
252
            dEd   = Two*C2*s
 
253
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
254
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
255
     &             + Four*at4*W3 + Five*at5*W4
 
256
     &             + Six*at6*W5 + Seven*at7*W6
 
257
     &             + Eight*at8*W7 + Nine*at9*W8
 
258
     &             + F10*at10*W9 + F11*at11*W10)
 
259
            dWdT = Two/((One + Tsig)**2)
 
260
            dTdR = (0.5d0*P32*rho13*rho13)/tauu
 
261
            dTdTau = -TauUEG/tauu**2
 
262
            dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
263
            dFdR = dFdW*dWdT*dTdR
 
264
            dFdTau=dFdW*dWdT*dTdTau
 
265
            dGGAdG =(fNL*dE*s/(Two*Gamma2))
 
266
            Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*(One+Fsig)
 
267
     &        + (fL+fNL*E)*Ax*rho43*dFdR
 
268
            Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + 
 
269
     &                    Two*dGGAdG*Ax*rho43*(One+Fsig) 
 
270
            Mmat(n,D1_TA)=  Mmat(n,D1_TA)
 
271
     &                   + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
272
 
 
273
 
 
274
 
 
275
10      continue
 
276
 
 
277
 
 
278
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
279
      else  ! ipol=2
 
280
c
 
281
c        ======> SPIN-UNRESTRICTED <======
 
282
 
 
283
c
 
284
c  use spin density functional theory ie n-->2n
 
285
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
286
 
 
287
         do 20 n = 1, nq
 
288
           if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
 
289
c
 
290
c     Alpha            ALPHA               ALPHA
 
291
c
 
292
            if (rho(n,R_A).lt.0.5d0*DTol) goto 25
 
293
             rhoo = Two*rho(n,R_A)
 
294
             rho43 = rhoo**F4o3
 
295
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
296
             rho13 = rho43*rrho
 
297
             rho53 = rhoo**F53
 
298
 
 
299
c
 
300
 
 
301
             tauN = tau(n,T_A)
 
302
             tauu = Two*tauN
 
303
             TauUEG=0.3d0*P32*rho53
 
304
             Tsig =TauUEG/tauu
 
305
             Wsig =(Tsig-One)/(Tsig+One)
 
306
             W1=Wsig
 
307
             W2=Wsig*W1
 
308
             W3=Wsig*W2
 
309
             W4=Wsig*W3
 
310
             W5=Wsig*W4
 
311
             W6=Wsig*W5
 
312
             W7=Wsig*W6
 
313
             W8=Wsig*W7
 
314
             W9=Wsig*W8
 
315
             W10=Wsig*W9
 
316
             W11=Wsig*W10
 
317
             Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
318
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
319
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
320
 
 
321
 
 
322
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
323
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
324
c    &              delrho(n,3,1)*delrho(n,3,1)
 
325
            Gamma2 = rgamma(n,G_AA)
 
326
            Gamma2 = Four*Gamma2
 
327
            Gamma = dsqrt(Gamma2)
 
328
            if (gamma.lt.DTol) goto 25
 
329
 
 
330
            s      = Cs*Gamma/rho43
 
331
            s2     = s*s
 
332
            En     = C1*s2
 
333
            Ed     = One + C2*s2
 
334
            E      = En/Ed
 
335
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n)
 
336
            func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0
 
337
c
 
338
c     functional derivatives
 
339
c
 
340
            dEn   = Two*C1*s
 
341
            dEd   = Two*C2*s
 
342
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
343
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
344
     &             + Four*at4*W3 + Five*at5*W4
 
345
     &             + Six*at6*W5 + Seven*at7*W6
 
346
     &             + Eight*at8*W7 + Nine*at9*W8
 
347
     &             + F10*at10*W9 + F11*at11*W10)
 
348
            dWdT = Two/((One + Tsig)**2)
 
349
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
350
            dTdTau = -Two*TauUEG/tauu**2
 
351
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
352
            dFdR = dFdW*dWdT*dTdR
 
353
            dFdTau=dFdW*dWdT*dTdTau
 
354
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
355
 
 
356
            Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(One+Fsig)
 
357
     &        + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
358
            Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + 
 
359
     &                      dGGAdG*Ax*rho43*(One+Fsig)*0.5d0
 
360
            Mmat(n,D1_TA)=  Mmat(n,D1_TA) + 
 
361
     &                  0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
362
 
 
363
c             write (*,*) "Ex,Amat(n,1),Cmat(n,1)",
 
364
c     &        Ex,Amat(n,1),Cmat(n,1)
 
365
 
 
366
c
 
367
c     Beta               BETA           BETA
 
368
c
 
369
 
 
370
25         continue
 
371
 
 
372
c
 
373
c     Beta
 
374
c
 
375
            if (rho(n,R_B).lt.0.5d0*DTol) goto 20
 
376
             rhoo = Two*rho(n,R_B)
 
377
             rho43 = rhoo**F4o3
 
378
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
379
             rho13 = rho43*rrho
 
380
             rho53 = rhoo**F53
 
381
 
 
382
c
 
383
 
 
384
             tauN = tau(n,T_B)
 
385
             tauu = Two*tauN
 
386
             TauUEG=0.3d0*P32*rho53
 
387
             Tsig =TauUEG/tauu
 
388
             Wsig =(Tsig-One)/(Tsig+One)
 
389
             W1=Wsig
 
390
             W2=Wsig*W1
 
391
             W3=Wsig*W2
 
392
             W4=Wsig*W3
 
393
             W5=Wsig*W4
 
394
             W6=Wsig*W5
 
395
             W7=Wsig*W6
 
396
             W8=Wsig*W7
 
397
             W9=Wsig*W8
 
398
             W10=Wsig*W9
 
399
             W11=Wsig*W10
 
400
             Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
401
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
402
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
403
 
 
404
 
 
405
c           Gamma2 = delrho(n,1,2)*delrho(n,1,2) +
 
406
c    &              delrho(n,2,2)*delrho(n,2,2) +
 
407
c    &              delrho(n,3,2)*delrho(n,3,2)
 
408
            Gamma2 = rgamma(n,G_BB)
 
409
            Gamma2 = Four*Gamma2
 
410
            Gamma = dsqrt(Gamma2)
 
411
            if (gamma.lt.DTol) goto 20
 
412
            s      = Cs*Gamma/rho43
 
413
            s2     = s*s
 
414
            En     = C1*s2
 
415
            Ed     = One + C2*s2
 
416
            E      = En/Ed
 
417
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n)
 
418
            func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0
 
419
c
 
420
c     functional derivatives
 
421
c
 
422
            dEn   = Two*C1*s
 
423
            dEd   = Two*C2*s
 
424
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
425
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
426
     &             + Four*at4*W3 + Five*at5*W4
 
427
     &             + Six*at6*W5 + Seven*at7*W6
 
428
     &             + Eight*at8*W7 + Nine*at9*W8
 
429
     &             + F10*at10*W9 + F11*at11*W10)
 
430
            dWdT = Two/((One + Tsig)**2)
 
431
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
432
            dTdTau = -Two*TauUEG/tauu**2
 
433
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
434
            dFdR = dFdW*dWdT*dTdR
 
435
            dFdTau=dFdW*dWdT*dTdTau
 
436
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
437
 
 
438
            Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(One+Fsig)
 
439
     &        + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
440
            Cmat(n,D1_GBB)=  Cmat(n,D1_GBB) + 
 
441
     &                   dGGAdG*Ax*rho43*(One+Fsig)*0.5d0
 
442
            Mmat(n,D1_TB)=  Mmat(n,D1_TB) +
 
443
     &                  0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
444
     
 
445
 
 
446
c
 
447
20      continue
 
448
      endif
 
449
c
 
450
      return
 
451
      end
 
452
 
 
453
 
 
454
 
 
455
 
 
456
      Subroutine nwxc_x_m05_d2()
 
457
      call errquit(' xm05: d2 not coded ',0,0)
 
458
      return
 
459
      end
 
460
 
 
461
 
 
462
c----------------------------------------------------------------------
 
463
C> \brief Calculate the dlDF correlation functional
 
464
C>
 
465
C> Calculate the dlDF correlation functional [1].
 
466
C>
 
467
C> ### References ###
 
468
C>
 
469
C> [1] K Pernal, R Podeszwa, K Patkowski, K Szalewicz,
 
470
C> "Dispersionless density functional theory",
 
471
C> Phys.Rev.Lett. <b>103</b>, 263201-263204 (2009), DOI:
 
472
C> <a href="http://dx.doi.org/10.1103/PhysRevLett.103.263201">
 
473
C> 10.1103/PhysRevLett.103.263201</a>.
 
474
C>
 
475
c   dlDF  exchange functional  
 
476
c           META GGA
 
477
C         utilizes ingredients:
 
478
c                              rho   -  density
 
479
c                              delrho - gradient of density
 
480
c                              tau - K.S kinetic energy density
 
481
c                              tauU - uniform-gas KE density
 
482
c                              
 
483
c     References: 
 
484
c     [a]       Pernal,Podeszwa,Patkowski,Szalewicz, PRL 103 263201 (2009)
 
485
 
 
486
      Subroutine nwxc_x_dldf(tol_rho, ipol, nq, wght, rho, rgamma, tau,
 
487
     &                       func, Amat, Cmat, Mmat)
 
488
c
 
489
      implicit none
 
490
c
 
491
#include "nwxc_param.fh"
 
492
c      
 
493
c
 
494
c     Input and other parameters
 
495
c
 
496
      double precision tol_rho !< [Input] The lower limit on the density
 
497
      integer nq               !< [Input] The number of points
 
498
      integer ipol             !< [Input] The number of spin channels
 
499
      double precision wght    !< [Input] The weight of the functional
 
500
c
 
501
c     Charge Density 
 
502
c
 
503
      double precision rho(nq,*) !< [Input] The density
 
504
c
 
505
c     Charge Density Gradient Norm
 
506
c
 
507
      double precision rgamma(nq,*) !< [Input] The density gradient norm
 
508
c
 
509
c     Kinetic Energy Density
 
510
c
 
511
      double precision tau(nq,*) !< [Input] The kinetic energy density
 
512
c
 
513
c     Functional values
 
514
c
 
515
      double precision func(*) !< [Output] The functional value
 
516
c
 
517
c     Sampling Matrices for the XC Potential
 
518
c
 
519
      double precision Amat(nq,*) !< [Output] Derivative wrt density
 
520
      double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma
 
521
      double precision Mmat(nq,*) !< [Output] Derivative wrt tau
 
522
c
 
523
      double precision pi
 
524
c
 
525
      integer n
 
526
      double precision at1, at2, at3, at4, at5, at6, at7, at8, at9
 
527
      double precision at, at10, at11, C1, C2, fL, fNL
 
528
      double precision rrho, rho43, rho13, rhoo, rho53
 
529
      double precision Gamma2, Gamma
 
530
      double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6
 
531
      double precision W7, W8, W9, W10, W11, Fsig
 
532
c
 
533
c     kinetic energy density   or  tau
 
534
c
 
535
      double precision tauN,tauu,DTol
 
536
      double precision F83, F23, F53, F43, F13, F1o3
 
537
      double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5
 
538
      double precision One, Two, Three, Four, Five, Six, Seven, Eight
 
539
      double precision Nine, F10, F11
 
540
      double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd
 
541
 
 
542
c      functional derivatives below FFFFFFFFFFFF
 
543
 
 
544
       double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR
 
545
       double precision dFdTau, dGGAdG,tauW
 
546
 
 
547
c      functional derivatives above FFFFFFFFFFFF
 
548
 
 
549
 
 
550
       parameter( pi = 3.1415926535897932384626433832795d0 )
 
551
         
 
552
        
 
553
        parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 
 
554
     &             F3o2=3.d0/2.d0,F13=1.d0/3.d0)
 
555
        parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0)
 
556
        parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0)
 
557
        parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 
 
558
     &             Five=5.0d0,Six=6.0d0, Seven=7.0d0,
 
559
     &             Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0)
 
560
   
 
561
c     Parameters for dlDF
 
562
        at1=    -1.637571d-01
 
563
        at2=    -1.880028d-01
 
564
        at3=    -4.490609d-01
 
565
        at4=    -8.2359d-03
 
566
        at5=     0.0d0
 
567
        at6=     0.0d0
 
568
        at7=     0.0d0
 
569
        at8=     0.0d0
 
570
        at9=     0.0d0
 
571
        at10=    0.0d0
 
572
        at11=    0.0d0
 
573
      
 
574
      
 
575
      at=1.0d0
 
576
      C1 = 0.3511128d0
 
577
      C2 = C1/4.8827323d0 
 
578
      DTol=tol_rho
 
579
C
 
580
C     Scale factors for local and non-local contributions.
 
581
C
 
582
      fL  =  wght
 
583
      fNL =  wght
 
584
      Cs = 0.5d0/(3.0d0*pi*pi)**F13
 
585
      P32 = (3.d0*pi**2)**F23
 
586
         
 
587
c     
 
588
       Ax = (-0.75d0)*(3.0d0/pi)**F13
 
589
 
 
590
 
 
591
c
 
592
      if (ipol.eq.1 )then
 
593
c
 
594
c        ======> SPIN-RESTRICTED <======
 
595
c                     or
 
596
c                SPIN-UNPOLARIZED
 
597
c
 
598
c
 
599
         do 10 n = 1, nq
 
600
           
 
601
            if (rho(n,R_T).lt.DTol) goto 10
 
602
            rhoo = rho(n,R_T)
 
603
            rho43 = rhoo**F4o3  
 
604
            rrho = 1d0/rhoo       ! reciprocal of rho
 
605
            rho13 = rho43*rrho
 
606
            rho53 = rhoo**F53
 
607
 
 
608
c
 
609
             
 
610
            tauN = tau(n,T_T)
 
611
            tauu=tauN 
 
612
cedo            if (taun.lt.sqrt(DTol)) goto 10
 
613
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
614
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
615
c    &              delrho(n,3,1)*delrho(n,3,1)
 
616
            Gamma2 = rgamma(n,G_TT)
 
617
            Gamma = dsqrt(Gamma2)
 
618
            if (gamma.lt.DTol) goto 10
 
619
            TauUEG=0.3d0*P32*rho53
 
620
            Tsig =TauUEG/tauu
 
621
            Wsig =(Tsig-One)/(Tsig+One)
 
622
            W1=Wsig 
 
623
            W2=Wsig*W1
 
624
            W3=Wsig*W2
 
625
            W4=Wsig*W3
 
626
            W5=Wsig*W4
 
627
            W6=Wsig*W5
 
628
            W7=Wsig*W6
 
629
            W8=Wsig*W7
 
630
            W9=Wsig*W8
 
631
            W10=Wsig*W9
 
632
            W11=Wsig*W10
 
633
            Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
634
     &          + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
635
     &          + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
636
 
 
637
            s      = Cs*Gamma/rho43
 
638
            s2     = s*s
 
639
            En     = C1*s2
 
640
            Ed     = One + C2*s2
 
641
            E      = En/Ed
 
642
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*qwght(n)
 
643
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)
 
644
c
 
645
c     functional derivatives 
 
646
c
 
647
            dEn   = Two*C1*s
 
648
            dEd   = Two*C2*s
 
649
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
650
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
651
     &             + Four*at4*W3 + Five*at5*W4
 
652
     &             + Six*at6*W5 + Seven*at7*W6
 
653
     &             + Eight*at8*W7 + Nine*at9*W8
 
654
     &             + F10*at10*W9 + F11*at11*W10)
 
655
            dWdT = Two/((One + Tsig)**2)
 
656
            dTdR = (0.5d0*P32*rho13*rho13)/tauu
 
657
            dTdTau = -TauUEG/tauu**2
 
658
            dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
659
            dFdR = dFdW*dWdT*dTdR
 
660
            dFdTau=dFdW*dWdT*dTdTau
 
661
            dGGAdG =(fNL*dE*s/(Two*Gamma2))
 
662
            Amat(n,D1_RA) = Amat(n,D1_RA) + dGGAdR*(One+Fsig)
 
663
     &        + (fL+fNL*E)*Ax*rho43*dFdR
 
664
            Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + 
 
665
     &                    Two*dGGAdG*Ax*rho43*(One+Fsig) 
 
666
            Mmat(n,D1_TA)=  Mmat(n,D1_TA)
 
667
     &                   + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
668
 
 
669
 
 
670
 
 
671
10      continue
 
672
 
 
673
 
 
674
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
675
      else  ! ipol=2
 
676
c
 
677
c        ======> SPIN-UNRESTRICTED <======
 
678
 
 
679
c
 
680
c  use spin density functional theory ie n-->2n
 
681
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
682
 
 
683
         do 20 n = 1, nq
 
684
           if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
 
685
c
 
686
c     Alpha            ALPHA               ALPHA
 
687
c
 
688
            if (rho(n,R_A).lt.0.5d0*DTol) goto 25
 
689
             rhoo = Two*rho(n,R_A)
 
690
             rho43 = rhoo**F4o3
 
691
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
692
             rho13 = rho43*rrho
 
693
             rho53 = rhoo**F53
 
694
 
 
695
c
 
696
 
 
697
             tauN = tau(n,T_A)
 
698
             tauu = Two*tauN
 
699
             TauUEG=0.3d0*P32*rho53
 
700
             Tsig =TauUEG/tauu
 
701
             Wsig =(Tsig-One)/(Tsig+One)
 
702
             W1=Wsig
 
703
             W2=Wsig*W1
 
704
             W3=Wsig*W2
 
705
             W4=Wsig*W3
 
706
             W5=Wsig*W4
 
707
             W6=Wsig*W5
 
708
             W7=Wsig*W6
 
709
             W8=Wsig*W7
 
710
             W9=Wsig*W8
 
711
             W10=Wsig*W9
 
712
             W11=Wsig*W10
 
713
             Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
714
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
715
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
716
 
 
717
 
 
718
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
719
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
720
c    &              delrho(n,3,1)*delrho(n,3,1)
 
721
            Gamma2 = rgamma(n,G_AA)
 
722
            Gamma2 = Four*Gamma2
 
723
            Gamma = dsqrt(Gamma2)
 
724
            if (gamma.lt.DTol) goto 25
 
725
 
 
726
            s      = Cs*Gamma/rho43
 
727
            s2     = s*s
 
728
            En     = C1*s2
 
729
            Ed     = One + C2*s2
 
730
            E      = En/Ed
 
731
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n)
 
732
            func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0
 
733
c
 
734
c     functional derivatives
 
735
c
 
736
            dEn   = Two*C1*s
 
737
            dEd   = Two*C2*s
 
738
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
739
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
740
     &             + Four*at4*W3 + Five*at5*W4
 
741
     &             + Six*at6*W5 + Seven*at7*W6
 
742
     &             + Eight*at8*W7 + Nine*at9*W8
 
743
     &             + F10*at10*W9 + F11*at11*W10)
 
744
            dWdT = Two/((One + Tsig)**2)
 
745
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
746
            dTdTau = -Two*TauUEG/tauu**2
 
747
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
748
            dFdR = dFdW*dWdT*dTdR
 
749
            dFdTau=dFdW*dWdT*dTdTau
 
750
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
751
 
 
752
            Amat(n,D1_RA) = Amat(n,D1_RA) + (dGGAdR*(One+Fsig)
 
753
     &        + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
754
            Cmat(n,D1_GAA)=  Cmat(n,D1_GAA) + 
 
755
     &                      dGGAdG*Ax*rho43*(One+Fsig)*0.5d0
 
756
            Mmat(n,D1_TA)=  Mmat(n,D1_TA) + 
 
757
     &                  0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
758
 
 
759
c             write (*,*) "Ex,Amat(n,1),Cmat(n,1)",
 
760
c     &        Ex,Amat(n,1),Cmat(n,1)
 
761
 
 
762
c
 
763
c     Beta               BETA           BETA
 
764
c
 
765
 
 
766
25         continue
 
767
 
 
768
c
 
769
c     Beta
 
770
c
 
771
            if (rho(n,R_B).lt.0.5d0*DTol) goto 20
 
772
             rhoo = Two*rho(n,R_B)
 
773
             rho43 = rhoo**F4o3
 
774
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
775
             rho13 = rho43*rrho
 
776
             rho53 = rhoo**F53
 
777
 
 
778
c
 
779
 
 
780
             tauN = tau(n,T_B)
 
781
             tauu = Two*tauN
 
782
             TauUEG=0.3d0*P32*rho53
 
783
             Tsig =TauUEG/tauu
 
784
             Wsig =(Tsig-One)/(Tsig+One)
 
785
             W1=Wsig
 
786
             W2=Wsig*W1
 
787
             W3=Wsig*W2
 
788
             W4=Wsig*W3
 
789
             W5=Wsig*W4
 
790
             W6=Wsig*W5
 
791
             W7=Wsig*W6
 
792
             W8=Wsig*W7
 
793
             W9=Wsig*W8
 
794
             W10=Wsig*W9
 
795
             W11=Wsig*W10
 
796
             Fsig =at*(at1*W1+ at2*W2 + at3*W3
 
797
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
798
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
799
 
 
800
 
 
801
c           Gamma2 = delrho(n,1,2)*delrho(n,1,2) +
 
802
c    &              delrho(n,2,2)*delrho(n,2,2) +
 
803
c    &              delrho(n,3,2)*delrho(n,3,2)
 
804
            Gamma2 = rgamma(n,G_BB)
 
805
            Gamma2 = Four*Gamma2
 
806
            Gamma = dsqrt(Gamma2)
 
807
            if (gamma.lt.DTol) goto 20
 
808
            s      = Cs*Gamma/rho43
 
809
            s2     = s*s
 
810
            En     = C1*s2
 
811
            Ed     = One + C2*s2
 
812
            E      = En/Ed
 
813
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(One+Fsig)*0.5d0*qwght(n)
 
814
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(One+Fsig)*.5d0
 
815
c
 
816
c     functional derivatives
 
817
c
 
818
            dEn   = Two*C1*s
 
819
            dEd   = Two*C2*s
 
820
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
821
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
822
     &             + Four*at4*W3 + Five*at5*W4
 
823
     &             + Six*at6*W5 + Seven*at7*W6
 
824
     &             + Eight*at8*W7 + Nine*at9*W8
 
825
     &             + F10*at10*W9 + F11*at11*W10)
 
826
            dWdT = Two/((One + Tsig)**2)
 
827
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
828
            dTdTau = -Two*TauUEG/tauu**2
 
829
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
830
            dFdR = dFdW*dWdT*dTdR
 
831
            dFdTau=dFdW*dWdT*dTdTau
 
832
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
833
 
 
834
            Amat(n,D1_RB) = Amat(n,D1_RB) + (dGGAdR*(One+Fsig)
 
835
     &        + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
836
            Cmat(n,D1_GBB)=  Cmat(n,D1_GBB) + 
 
837
     &                   dGGAdG*Ax*rho43*(One+Fsig)*0.5d0
 
838
            Mmat(n,D1_TB)=  Mmat(n,D1_TB) +
 
839
     &                  0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
840
     
 
841
 
 
842
c
 
843
20      continue
 
844
      endif
 
845
c
 
846
      return
 
847
      end
 
848
 
 
849
 
 
850
 
 
851
      Subroutine nwxc_x_dldf_d2()
 
852
      call errquit('xdldf: d2 not coded',0,0)
 
853
      return
 
854
      end
 
855
C> @}
 
856
 
 
857