← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_xm11.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/nwdft/xc/xc_xm11.F
 
1
      Subroutine xc_xm11(tol_rho, fac,lfac,nlfac, rho, delrho, 
 
2
     &                     Amat, Cmat, nq, ipol, Ex, 
 
3
     &                     qwght, ldew, func, tau, Mmat, ijzy)
 
4
 
 
5
 
 
6
c   
 
7
c$Id: xc_xm11.F 23096 2012-11-14 01:09:04Z edo $
 
8
c
 
9
c
 
10
c
 
11
c**********************************************************************c
 
12
c                                                                      c
 
13
c  xc_xm11 evaluates the exchange part of the M08 and M11 suite of     c
 
14
c  functionals on the grid.                                            c
 
15
c  !!! Second derivatives are not available yet.                       c
 
16
c                                                                      c
 
17
c  Ref: (a) Zhao, Y.  and Truhlar, D. G. JCTC, 2008, 4 , 1849          c
 
18
c       (b) Peverati, R. and Truhlar, D. G. J.P.C.Lett. 2011, 2, 2810  c
 
19
c       (c) Peverati, R. and Truhlar, D. G. J.P.C.Lett. 2011, 3, 117   c
 
20
c                                                                      c
 
21
c       ijzy - 1 M08-HX (a)                                            c
 
22
c       ijzy - 2 M08-SO (a)                                            c
 
23
c       ijzy - 3 M11 (b)                                               c
 
24
c       ijzy - 4 M11-L (c)                                             c
 
25
c                                                                      c
 
26
c Coded by Roberto Peverati (12/11)                                    c
 
27
c                                                                      c
 
28
c**********************************************************************c    
 
29
c
 
30
      implicit none
 
31
c      
 
32
      double precision fac, Ex
 
33
      integer nq, ipol
 
34
      logical lfac, nlfac,ldew,   uselc
 
35
      double precision func(*)  ! value of the functional [output]
 
36
c
 
37
c     Charge Density & Its Cube Root
 
38
c
 
39
      double precision rho(nq,ipol*(ipol+1)/2)
 
40
c
 
41
c     Charge Density Gradient
 
42
c
 
43
      double precision delrho(nq,3,ipol)
 
44
c
 
45
c     Quadrature Weights
 
46
c
 
47
      double precision qwght(nq)
 
48
c
 
49
c     Sampling Matrices for the XC Potential & Energy
 
50
c
 
51
      double precision Amat(nq,ipol), Cmat(nq,*), Mmat(nq,*)
 
52
      double precision tol_rho, pi
 
53
c
 
54
c     kinetic energy density   or  tau
 
55
c
 
56
      double precision tau(nq,ipol)
 
57
      double precision tauN,tauu,DTol
 
58
c
 
59
c      functional derivatives
 
60
c
 
61
      double precision dWdT, dTdR, dTdTau
 
62
c
 
63
c     Intermediate derivative results, etc.
 
64
c
 
65
      integer n, ijzy
 
66
c
 
67
      double precision Ax, s, s2
 
68
      double precision kapa,kapas,mu,mus
 
69
c
 
70
      double precision f0,f1,f2,f3,f4,f5,f6,f7,f8,f9,f10,f11
 
71
      double precision F1o3,F2o3,F3o5,F4o3,F5o3,F48,F81 
 
72
      double precision Fsig1, Fsig2, Fsig3, Fsig4, Fx1, Fx2
 
73
      double precision ElSR, ElLR
 
74
      double precision PDUM
 
75
      double precision GGA1, GGA2, GGA3, GGA4
 
76
      double precision Emu, X, deno
 
77
      double precision ds2drho, ds2dg, dfx1ds2
 
78
      double precision dfx2ds2, df1dw
 
79
      double precision dfx1drho,dfx1dg,dfx2drho,dfx2dg,df2dw
 
80
      double precision df3dw, df4dw, delsrdr, dellrdr
 
81
      double precision dgga1dr, dgga2dr, dgga3dr, dgga4dr
 
82
      double precision df1dr, df1dtau, df2dr, df2dtau
 
83
      double precision df3dr, df3dtau, df4dr, df4dtau
 
84
      double precision dgga1dg, dgga2dg, dgga3dg, dgga4dg
 
85
c
 
86
      double precision at00,at01,at02,at03,at04,at05,at06
 
87
      double precision at07,at08,at09,at10,at11
 
88
      double precision bt00,bt01,bt02,bt03,bt04,bt05,bt06
 
89
      double precision bt07,bt08,bt09,bt10,bt11
 
90
      double precision ct00,ct01,ct02,ct03,ct04,ct05,ct06
 
91
      double precision ct07,ct08,ct09,ct10,ct11
 
92
      double precision dt00,dt01,dt02,dt03,dt04,dt05,dt06
 
93
      double precision dt07,dt08,dt09,dt10,dt11
 
94
      double precision rho43, rho13, rhoo, rho53
 
95
      double precision Gamma2, Gamma
 
96
      double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6
 
97
      double precision W7, W8, W9, W10, W11
 
98
c
 
99
      parameter( F0=0.0D+00,  F1=1.0D+00,  F2=2.0D+00,
 
100
     $           F3=3.0D+00,  F4=4.0D+00,  F5=5.0D+00,
 
101
     $           F6=6.0D+00,  F7=7.0D+00,  F8=8.0D+00,
 
102
     $           F9=9.0D+00,  F10=10.0D+00,F11=11.0D+00)
 
103
c
 
104
        pi=acos(-1d0)      
 
105
c
 
106
        ct00= 0D+00
 
107
        ct01= 0D+00
 
108
        ct02= 0D+00
 
109
        ct03= 0D+00
 
110
        ct04= 0D+00
 
111
        ct05= 0D+00
 
112
        ct06= 0D+00
 
113
        ct07= 0D+00
 
114
        ct08= 0D+00
 
115
        ct09= 0D+00
 
116
        ct10= 0D+00
 
117
        ct11= 0D+00
 
118
C
 
119
        dt00= 0D+00
 
120
        dt01= 0D+00
 
121
        dt02= 0D+00
 
122
        dt03= 0D+00
 
123
        dt04= 0D+00
 
124
        dt05= 0D+00
 
125
        dt06= 0D+00
 
126
        dt07= 0D+00
 
127
        dt08= 0D+00
 
128
        dt09= 0D+00
 
129
        dt10= 0D+00
 
130
        dt11= 0D+00
 
131
        at00=0.000000D+00
 
132
        at01=0.000000D+00
 
133
        at02=0.000000D+00
 
134
        at03=0.000000D+00
 
135
        at04=0.000000D+00
 
136
        at05=0.000000D+00
 
137
        at06=0.000000D+00
 
138
        at07=0.000000D+00
 
139
        at08=0.000000D+00
 
140
        at09=0.000000D+00
 
141
        at10=0.000000D+00 
 
142
        at11=0.000000D+00 
 
143
        bt00=0.000000D+00
 
144
        bt01=0.000000D+00
 
145
        bt02=0.000000D+00
 
146
        bt03=0.000000D+00
 
147
        bt04=0.000000D+00
 
148
        bt05=0.000000D+00
 
149
        bt06=0.000000D+00
 
150
        bt07=0.000000D+00
 
151
        bt08=0.000000D+00
 
152
        bt09=0.000000D+00
 
153
        bt10=0.000000D+00 
 
154
        bt11=0.000000D+00 
 
155
        UseLC=.False.
 
156
C
 
157
      if (ijzy.eq.1) then
 
158
C     Parameters for M08-HX
 
159
        at00=  1.3340172D+00
 
160
        at01= -9.4751087D+00
 
161
        at02= -1.2541893D+01
 
162
        at03=  9.1369974D+00
 
163
        at04=  3.4717204D+01
 
164
        at05=  5.8831807D+01
 
165
        at06=  7.1369574D+01
 
166
        at07=  2.3312961D+01
 
167
        at08=  4.8314679D+00
 
168
        at09= -6.5044167D+00
 
169
        at10= -1.4058265D+01
 
170
        at11=  1.2880570D+01
 
171
 
 
172
        bt00= -8.5631823D-01
 
173
        bt01=  9.2810354D+00
 
174
        bt02=  1.2260749D+01
 
175
        bt03= -5.5189665D+00
 
176
        bt04= -3.5534989D+01
 
177
        bt05= -8.2049996D+01
 
178
        bt06= -6.8586558D+01
 
179
        bt07=  3.6085694D+01
 
180
        bt08= -9.3740983D+00
 
181
        bt09= -5.9731688D+01
 
182
        bt10=  1.6587868D+01
 
183
        bt11=  1.3993203D+01
 
184
C
 
185
        UseLC=.False.
 
186
C
 
187
       elseif (ijzy.eq.2) then
 
188
C     Parameters for M08-SO
 
189
        at00= -3.4888428D-01
 
190
        at01= -5.8157416D+00
 
191
        at02=  3.7550810D+01
 
192
        at03=  6.3727406D+01
 
193
        at04= -5.3742313D+01
 
194
        at05= -9.8595529D+01
 
195
        at06=  1.6282216D+01
 
196
        at07=  1.7513468D+01
 
197
        at08= -6.7627553D+00
 
198
        at09=  1.1106658D+01
 
199
        at10=  1.5663545D+00
 
200
        at11=  8.7603470D+00
 
201
 
 
202
        bt00=  7.8098428D-01
 
203
        bt01=  5.4538178D+00
 
204
        bt02= -3.7853348D+01
 
205
        bt03= -6.2295080D+01
 
206
        bt04=  4.6713254D+01
 
207
        bt05=  8.7321376D+01
 
208
        bt06=  1.6053446D+01
 
209
        bt07=  2.0126920D+01
 
210
        bt08= -4.0343695D+01
 
211
        bt09= -5.8577565D+01
 
212
        bt10=  2.0890272D+01
 
213
        bt11=  1.0946903D+01
 
214
C
 
215
        UseLC=.False.
 
216
C
 
217
      elseif (ijzy.eq.3) then
 
218
C     Parameters for M11
 
219
        at00= -0.18399900D+00
 
220
        at01= -1.39046703D+01
 
221
        at02=  1.18206837D+01
 
222
        at03=  3.10098465D+01
 
223
        at04= -5.19625696D+01
 
224
        at05=  1.55750312D+01
 
225
        at06= -6.94775730D+00
 
226
        at07= -1.58465014D+02
 
227
        at08= -1.48447565D+00
 
228
        at09=  5.51042124D+01
 
229
        at10= -1.34714184D+01
 
230
        at11=  0.00000000D+00
 
231
 
 
232
        bt00=  0.75599900D+00
 
233
        bt01=  1.37137944D+01
 
234
        bt02= -1.27998304D+01
 
235
        bt03= -2.93428814D+01
 
236
        bt04=  5.91075674D+01
 
237
        bt05= -2.27604866D+01
 
238
        bt06= -1.02769340D+01
 
239
        bt07=  1.64752731D+02
 
240
        bt08=  1.85349258D+01
 
241
        bt09= -5.56825639D+01
 
242
        bt10=  7.47980859D+00
 
243
        bt11=  0.00000000D+00
 
244
C
 
245
        UseLC=.True.
 
246
        Emu =0.25D+00
 
247
C
 
248
      elseif (ijzy.eq.4) then
 
249
C     Parameters for M11-L
 
250
        at00=  8.121131D-01
 
251
        at01=  1.738124D+01
 
252
        at02=  1.154007D+00
 
253
        at03=  6.869556D+01
 
254
        at04=  1.016864D+02
 
255
        at05= -5.887467D+00
 
256
        at06=  4.517409D+01
 
257
        at07= -2.773149D+00
 
258
        at08= -2.617211D+01
 
259
        at09=  0.000000D+00
 
260
        at10=  0.000000D+00 
 
261
        at11=  0.000000D+00
 
262
C
 
263
        bt00=  1.878869D-01
 
264
        bt01= -1.653877D+01
 
265
        bt02=  6.755753D-01
 
266
        bt03= -7.567572D+01
 
267
        bt04= -1.040272D+02
 
268
        bt05=  1.831853D+01
 
269
        bt06= -5.573352D+01
 
270
        bt07= -3.520210D+00
 
271
        bt08=  3.724276D+01
 
272
        bt09=  0.000000D+00
 
273
        bt10=  0.000000D+00
 
274
        bt11=  0.000000D+00
 
275
C
 
276
        ct00= -4.386615D-01
 
277
        ct01= -1.214016D+02
 
278
        ct02= -1.393573D+02
 
279
        ct03= -2.046649D+00
 
280
        ct04=  2.804098D+01
 
281
        ct05= -1.312258D+01
 
282
        ct06= -6.361819D+00
 
283
        ct07= -8.055758D-01
 
284
        ct08=  3.736551D+00
 
285
        ct09=  0.000000D+00
 
286
        ct10=  0.000000D+00
 
287
        ct11=  0.000000D+00
 
288
C
 
289
        dt00=  1.438662D+00
 
290
        dt01=  1.209465D+02
 
291
        dt02=  1.328252D+02
 
292
        dt03=  1.296355D+01
 
293
        dt04=  5.854866D+00
 
294
        dt05= -3.378162D+00
 
295
        dt06= -4.423393D+01
 
296
        dt07=  6.844475D+00
 
297
        dt08=  1.949541D+01
 
298
        dt09=  0.000000D+00
 
299
        dt10=  0.000000D+00
 
300
        dt11=  0.000000D+00
 
301
C
 
302
        UseLC=.True.
 
303
        Emu =0.25D+00
 
304
C
 
305
      endif
 
306
      DTol=tol_rho
 
307
      F1o3 = F1/F3 
 
308
      F2o3 = F2/F3
 
309
      F3o5 = F3/F5
 
310
      F4o3 = F4/F3 
 
311
      F5o3 = F5/F3
 
312
      F48 = 48.0d0
 
313
      F81 = 81.0d0
 
314
      Ax = -(F3/F2) * (F4o3*Pi)**(-F1o3) 
 
315
C     RPBE parameters
 
316
      Mus = F10/F81
 
317
      kapas = 0.552d0
 
318
C     PBE parameters 
 
319
      Mu = 0.21951d0
 
320
      kapa = 0.804d0
 
321
c
 
322
      if (ipol.eq.1 )then
 
323
c
 
324
c        ======> SPIN-RESTRICTED <======
 
325
c                     or
 
326
c                SPIN-UNPOLARIZED
 
327
c
 
328
c
 
329
         do 10 n = 1, nq
 
330
            if (rho(n,1).lt.DTol) goto 10
 
331
            rhoo = rho(n,1)/F2
 
332
            rho43 = rhoo**F4o3  
 
333
            rho13 = rho43/rhoo
 
334
            rho53 = rhoo**F5o3
 
335
c            
 
336
            tauN = tau(n,1)
 
337
         if(taun.lt.dtol) goto 10
 
338
            tauu=tauN
 
339
            TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53
 
340
            Tsig =TauUEG/tauu
 
341
            Wsig =(Tsig - F1)/(Tsig + F1)
 
342
            W1=Wsig 
 
343
            W2=Wsig*W1
 
344
            W3=Wsig*W2
 
345
            W4=Wsig*W3
 
346
            W5=Wsig*W4
 
347
            W6=Wsig*W5
 
348
            W7=Wsig*W6
 
349
            W8=Wsig*W7
 
350
            W9=Wsig*W8
 
351
            W10=Wsig*W9
 
352
            W11=Wsig*W10
 
353
            Fsig1 =(at00    + at01*W1 + at02*W2 + at03*W3
 
354
     $            + at04*W4 + at05*W5 + at06*W6 + at07*W7
 
355
     $            + at08*W8 + at09*W9 + at10*W10+ at11*W11)
 
356
            Fsig2 =(bt00    + bt01*W1 + bt02*W2 + bt03*W3
 
357
     $            + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7
 
358
     $            + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11)
 
359
            Fsig3 =(ct00    + ct01*W1 + ct02*W2 + ct03*W3
 
360
     $            + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7
 
361
     $            + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11)
 
362
            Fsig4 =(dt00    + dt01*W1 + dt02*W2 + dt03*W3
 
363
     $            + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7
 
364
     $            + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11)
 
365
 
 
366
            Gamma2 =(delrho(n,1,1)*delrho(n,1,1) +
 
367
     &              delrho(n,2,1)*delrho(n,2,1) +
 
368
     &              delrho(n,3,1)*delrho(n,3,1))/F4
 
369
            Gamma = dsqrt(Gamma2)
 
370
         if(gamma.lt.dtol) goto 10
 
371
         
 
372
         X = GAMMA/RHO43
 
373
         S = X/(F48*PI*PI)**F1o3
 
374
         s2     = s*s
 
375
         Deno = (F1 + Mu*s2/kapa)
 
376
         fx1=F1+kapa*(F1-F1/Deno)
 
377
         fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas))
 
378
         If(UseLC) then
 
379
           CALL LRCLSDA(EMU,RHOO,ElSR,PDUM)
 
380
           ElLR = Ax*Rho43-ElSR
 
381
         else
 
382
           ElSR = Ax*Rho43
 
383
           ElLR = F0
 
384
         endIf
 
385
         GGA1 = ElSR*fx1
 
386
         GGA2 = ElSR*fx2
 
387
         GGA3 = ElLR*fx1
 
388
         GGA4 = ElLR*fx2
 
389
C
 
390
          Ex = Ex +F2*(GGA1*Fsig1 + GGA2*Fsig2
 
391
     $            +    GGA3*Fsig3 + GGA4*Fsig4)*qwght(n)
 
392
          if(ldew) func(n)=func(n)+F2*(GGA1*Fsig1+GGA2*Fsig2
 
393
     $                            +    GGA3*Fsig3+GGA4*Fsig4)
 
394
 
 
395
c
 
396
c     functional derivatives 
 
397
c
 
398
            ds2dRho = -(F8/F3) * s2/rhoo
 
399
            ds2dG = s2/Gamma2
 
400
C
 
401
            dfx1ds2 = Mu*(F1/(Deno*Deno)) 
 
402
            dfx1dRho = dfx1ds2*ds2dRho
 
403
            dfx1dG = dfx1ds2*ds2dG
 
404
C
 
405
            dfx2ds2 = Mus*DExp(-Mus*s2/kapas)
 
406
            dfx2dRho = dfx2ds2*ds2dRho
 
407
            dfx2dG = dfx2ds2*ds2dG
 
408
c
 
409
            dF1dW = (at01 + F2*at02*W1 + F3*at03*W2
 
410
     $                    + F4*at04*W3 + F5*at05*W4
 
411
     $                    + F6*at06*W5 + F7*at07*W6
 
412
     $                    + F8*at08*W7 + F9*at09*W8
 
413
     $                    + F10*at10*W9+F11*at11*W10)
 
414
            dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2
 
415
     $                    + F4*bt04*W3 + F5*bt05*W4
 
416
     $                    + F6*bt06*W5 + F7*bt07*W6
 
417
     $                    + F8*bt08*W7 + F9*bt09*W8
 
418
     $                    + F10*Bt10*W9+F11*Bt11*W10)
 
419
            dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2
 
420
     $                    + F4*ct04*W3 + F5*ct05*W4
 
421
     $                    + F6*ct06*W5 + F7*ct07*W6
 
422
     $                    + F8*ct08*W7 + F9*ct09*W8
 
423
     $                    + F10*ct10*W9+F11*ct11*W10)
 
424
            dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2
 
425
     $                    + F4*dt04*W3 + F5*dt05*W4
 
426
     $                    + F6*dt06*W5 + F7*dt07*W6
 
427
     $                    + F8*dt08*W7 + F9*dt09*W8
 
428
     $                    + F10*dt10*W9+F11*dt11*W10)
 
429
c
 
430
            dWdT = F2/((F1 + Tsig)**F2)
 
431
            dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN
 
432
            dTdTau = -TauUEG/tauN**F2
 
433
C
 
434
           If(UseLC) then
 
435
             dElSRdR = PDUM
 
436
             dElLRdR = Ax*F4o3*Rho13-PDUM
 
437
           else
 
438
             dElSRdR=Ax*F4o3*Rho13
 
439
             dElLRdR=F0
 
440
           endIf  
 
441
           dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho
 
442
           dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 
 
443
           dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho
 
444
           dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 
 
445
c
 
446
           dF1dR = dF1dW*dWdT*dTdR
 
447
           dF1dTau=dF1dW*dWdT*dTdTau
 
448
           dF2dR = dF2dW*dWdT*dTdR
 
449
           dF2dTau=dF2dW*dWdT*dTdTau
 
450
           dF3dR = dF3dW*dWdT*dTdR
 
451
           dF3dTau=dF3dW*dWdT*dTdTau
 
452
           dF4dR = dF4dW*dWdT*dTdR
 
453
           dF4dTau=dF4dW*dWdT*dTdTau
 
454
c
 
455
           dGGA1dG = ElSR*dfx1dG
 
456
           dGGA2dG = ElSR*dfx2dG
 
457
           dGGA3dG = ElLR*dfx1dG
 
458
           dGGA4dG = ElLR*dfx2dG
 
459
c
 
460
           Amat(n,1) = Amat(n,1)   +dGGA1dR*Fsig1 + GGA1*dF1dR
 
461
     $                             +dGGA2dR*Fsig2 + GGA2*dF2dR
 
462
     $                             +dGGA3dR*Fsig3 + GGA3*dF3dR
 
463
     $                             +dGGA4dR*Fsig4 + GGA4*dF4dR
 
464
           Cmat(n,1)=  Cmat(n,1)  +dGGA1dG*Fsig1 + dGGA2dG*Fsig2
 
465
     $                            +dGGA3dG*Fsig3 + dGGA4dG*Fsig4
 
466
           Mmat(n,1)=  Mmat(n,1)   +GGA1*dF1dTau + GGA2*dF2dTau
 
467
     $                             +GGA3*dF3dTau + GGA4*dF4dTau
 
468
c    
 
469
10      continue
 
470
c
 
471
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
472
      else  ! ipol=2
 
473
c
 
474
c        ======> SPIN-UNRESTRICTED <======
 
475
 
 
476
c
 
477
c  use spin density functional theory ie n-->2n
 
478
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
479
c
 
480
c     Alpha            ALPHA               ALPHA
 
481
c
 
482
         do 20 n = 1, nq
 
483
           if (rho(n,1).lt.DTol) goto 20
 
484
           if (rho(n,2).lt.DTol) goto 25           
 
485
            rhoo  = rho(n,2)
 
486
            rho43 = rhoo**F4o3  
 
487
            rho13 = rho43/rhoo
 
488
            rho53 = rhoo**F5o3
 
489
c       
 
490
            tauN = tau(n,1)*F2
 
491
     
 
492
         if(taun.lt.dtol) goto 25
 
493
            tauu=tauN
 
494
            TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53
 
495
            Tsig =TauUEG/tauu
 
496
            Wsig =(Tsig - F1)/(Tsig + F1)
 
497
            W1=Wsig 
 
498
            W2=Wsig*W1
 
499
            W3=Wsig*W2
 
500
            W4=Wsig*W3
 
501
            W5=Wsig*W4
 
502
            W6=Wsig*W5
 
503
            W7=Wsig*W6
 
504
            W8=Wsig*W7
 
505
            W9=Wsig*W8
 
506
            W10=Wsig*W9
 
507
            W11=Wsig*W10
 
508
            Fsig1 =(at00    + at01*W1 + at02*W2 + at03*W3
 
509
     $            + at04*W4 + at05*W5 + at06*W6 + at07*W7
 
510
     $            + at08*W8 + at09*W9 + at10*W10+ at11*W11)
 
511
            Fsig2 =(bt00    + bt01*W1 + bt02*W2 + bt03*W3
 
512
     $            + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7
 
513
     $            + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11)
 
514
            Fsig3 =(ct00    + ct01*W1 + ct02*W2 + ct03*W3
 
515
     $            + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7
 
516
     $            + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11)
 
517
            Fsig4 =(dt00    + dt01*W1 + dt02*W2 + dt03*W3
 
518
     $            + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7
 
519
     $            + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11)
 
520
 
 
521
            Gamma2 =(delrho(n,1,1)*delrho(n,1,1) +
 
522
     &              delrho(n,2,1)*delrho(n,2,1) +
 
523
     &              delrho(n,3,1)*delrho(n,3,1))
 
524
            Gamma = dsqrt(Gamma2)
 
525
         if(gamma.lt.dtol) goto 25
 
526
         
 
527
         X = GAMMA/RHO43
 
528
         S = X/(F48*PI*PI)**F1o3
 
529
         s2     = s*s
 
530
         Deno = (F1 + Mu*s2/kapa)
 
531
         fx1=F1+kapa*(F1-F1/Deno)
 
532
         fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas))
 
533
         If(UseLC) then
 
534
           CALL LRCLSDA(EMU,RHOO,ElSR,PDUM)
 
535
           ElLR = Ax*Rho43-ElSR
 
536
         else
 
537
           ElSR = Ax*Rho43
 
538
           ElLR = F0
 
539
         endIf
 
540
         GGA1 = ElSR*fx1
 
541
         GGA2 = ElSR*fx2
 
542
         GGA3 = ElLR*fx1
 
543
         GGA4 = ElLR*fx2
 
544
C
 
545
          Ex = Ex +   (GGA1*Fsig1 + GGA2*Fsig2
 
546
     $            +    GGA3*Fsig3 + GGA4*Fsig4)*qwght(n)
 
547
          if(ldew) func(n)=func(n)+   (GGA1*Fsig1+GGA2*Fsig2
 
548
     $                            +    GGA3*Fsig3+GGA4*Fsig4)
 
549
c
 
550
c     functional derivatives 
 
551
c
 
552
            ds2dRho = -(F8/F3) * s2/rhoo
 
553
            ds2dG = s2/Gamma2
 
554
C
 
555
            dfx1ds2 = Mu*(F1/(Deno*Deno)) 
 
556
            dfx1dRho = dfx1ds2*ds2dRho
 
557
            dfx1dG = dfx1ds2*ds2dG
 
558
C
 
559
            dfx2ds2 = Mus*DExp(-Mus*s2/kapas)
 
560
            dfx2dRho = dfx2ds2*ds2dRho
 
561
            dfx2dG = dfx2ds2*ds2dG
 
562
c
 
563
            dF1dW = (at01 + F2*at02*W1 + F3*at03*W2
 
564
     $                    + F4*at04*W3 + F5*at05*W4
 
565
     $                    + F6*at06*W5 + F7*at07*W6
 
566
     $                    + F8*at08*W7 + F9*at09*W8
 
567
     $                    + F10*at10*W9+F11*at11*W10)
 
568
            dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2
 
569
     $                    + F4*bt04*W3 + F5*bt05*W4
 
570
     $                    + F6*bt06*W5 + F7*bt07*W6
 
571
     $                    + F8*bt08*W7 + F9*bt09*W8
 
572
     $                    + F10*Bt10*W9+F11*Bt11*W10)
 
573
            dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2
 
574
     $                    + F4*ct04*W3 + F5*ct05*W4
 
575
     $                    + F6*ct06*W5 + F7*ct07*W6
 
576
     $                    + F8*ct08*W7 + F9*ct09*W8
 
577
     $                    + F10*ct10*W9+F11*ct11*W10)
 
578
            dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2
 
579
     $                    + F4*dt04*W3 + F5*dt05*W4
 
580
     $                    + F6*dt06*W5 + F7*dt07*W6
 
581
     $                    + F8*dt08*W7 + F9*dt09*W8
 
582
     $                    + F10*dt10*W9+F11*dt11*W10)
 
583
 
 
584
            dWdT = F2/((F1 + Tsig)**F2)
 
585
            dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN
 
586
            dTdTau = -TauUEG/tauN**F2
 
587
C
 
588
           If(UseLC) then
 
589
             dElSRdR = PDUM
 
590
             dElLRdR = Ax*F4o3*Rho13-PDUM
 
591
           else
 
592
             dElSRdR=Ax*F4o3*Rho13
 
593
             dElLRdR=F0
 
594
           endIf  
 
595
           dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho
 
596
           dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 
 
597
           dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho
 
598
           dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 
 
599
c
 
600
           dF1dR = dF1dW*dWdT*dTdR
 
601
           dF1dTau=dF1dW*dWdT*dTdTau
 
602
           dF2dR = dF2dW*dWdT*dTdR
 
603
           dF2dTau=dF2dW*dWdT*dTdTau
 
604
           dF3dR = dF3dW*dWdT*dTdR
 
605
           dF3dTau=dF3dW*dWdT*dTdTau
 
606
           dF4dR = dF4dW*dWdT*dTdR
 
607
           dF4dTau=dF4dW*dWdT*dTdTau
 
608
c
 
609
           dGGA1dG = ElSR*dfx1dG
 
610
           dGGA2dG = ElSR*dfx2dG
 
611
           dGGA3dG = ElLR*dfx1dG
 
612
           dGGA4dG = ElLR*dfx2dG
 
613
c
 
614
           Amat(n,1) = Amat(n,1)   +dGGA1dR*Fsig1 + GGA1*dF1dR
 
615
     $                             +dGGA2dR*Fsig2 + GGA2*dF2dR
 
616
     $                             +dGGA3dR*Fsig3 + GGA3*dF3dR
 
617
     $                             +dGGA4dR*Fsig4 + GGA4*dF4dR
 
618
           Cmat(n,1)=  Cmat(n,1)   +dGGA1dG*Fsig1 + dGGA2dG*Fsig2
 
619
     $                             +dGGA3dG*Fsig3 + dGGA4dG*Fsig4
 
620
           Mmat(n,1)=  Mmat(n,1)   +GGA1*dF1dTau  + GGA2*dF2dTau
 
621
     $                             +GGA3*dF3dTau  + GGA4*dF4dTau
 
622
c
 
623
25         continue
 
624
c
 
625
c     Beta               BETA           BETA
 
626
c
 
627
            if (rho(n,3).lt.DTol) goto 20
 
628
            rhoo  = rho(n,3)
 
629
            rho43 = rhoo**F4o3  
 
630
            rho13 = rho43/rhoo
 
631
            rho53 = rhoo**F5o3
 
632
c
 
633
             
 
634
            tauN = tau(n,2)*F2
 
635
     
 
636
         if(taun.lt.dtol) goto 20
 
637
            tauu=tauN
 
638
            TAUUEG=F3O5*((F6*PI*PI)**F2O3)*RHO53
 
639
            Tsig =TauUEG/tauu
 
640
            Wsig =(Tsig - F1)/(Tsig + F1)
 
641
            W1=Wsig 
 
642
            W2=Wsig*W1
 
643
            W3=Wsig*W2
 
644
            W4=Wsig*W3
 
645
            W5=Wsig*W4
 
646
            W6=Wsig*W5
 
647
            W7=Wsig*W6
 
648
            W8=Wsig*W7
 
649
            W9=Wsig*W8
 
650
            W10=Wsig*W9
 
651
            W11=Wsig*W10
 
652
            Fsig1 =(at00    + at01*W1 + at02*W2 + at03*W3
 
653
     $            + at04*W4 + at05*W5 + at06*W6 + at07*W7
 
654
     $            + at08*W8 + at09*W9 + at10*W10+ at11*W11)
 
655
            Fsig2 =(bt00    + bt01*W1 + bt02*W2 + bt03*W3
 
656
     $            + bt04*W4 + bt05*W5 + bt06*W6 + bt07*W7
 
657
     $            + bt08*W8 + bt09*W9 + bt10*W10+ bt11*W11)
 
658
            Fsig3 =(ct00    + ct01*W1 + ct02*W2 + ct03*W3
 
659
     $            + ct04*W4 + ct05*W5 + ct06*W6 + ct07*W7
 
660
     $            + ct08*W8 + ct09*W9 + ct10*W10+ ct11*W11)
 
661
            Fsig4 =(dt00    + dt01*W1 + dt02*W2 + dt03*W3
 
662
     $            + dt04*W4 + dt05*W5 + dt06*W6 + dt07*W7
 
663
     $            + dt08*W8 + dt09*W9 + dt10*W10+ dt11*W11)
 
664
 
 
665
            Gamma2 =(delrho(n,1,2)*delrho(n,1,2) +
 
666
     &              delrho(n,2,2)*delrho(n,2,2) +
 
667
     &              delrho(n,3,2)*delrho(n,3,2))
 
668
            Gamma = dsqrt(Gamma2)
 
669
         if(gamma.lt.dtol) goto 20
 
670
         
 
671
         X = GAMMA/RHO43
 
672
         S = X/(F48*PI*PI)**F1o3
 
673
         s2     = s*s
 
674
         Deno = (F1 + Mu*s2/kapa)
 
675
         fx1=F1+kapa*(F1-F1/Deno)
 
676
         fx2=F1+kapas*(F1-DExp(-Mus*s2/kapas))
 
677
         If(UseLC) then
 
678
           CALL LRCLSDA(EMU,RHOO,ElSR,PDUM)
 
679
           ElLR = Ax*Rho43-ElSR
 
680
         else
 
681
           ElSR = Ax*Rho43
 
682
           ElLR = F0
 
683
         endIf
 
684
         GGA1 = ElSR*fx1
 
685
         GGA2 = ElSR*fx2
 
686
         GGA3 = ElLR*fx1
 
687
         GGA4 = ElLR*fx2
 
688
C
 
689
          Ex = Ex +   (GGA1*Fsig1 + GGA2*Fsig2
 
690
     $            +    GGA3*Fsig3 + GGA4*Fsig4)*qwght(n)
 
691
          if(ldew) func(n)=func(n)+   (GGA1*Fsig1+GGA2*Fsig2
 
692
     $                            +    GGA3*Fsig3+GGA4*Fsig4)
 
693
 
 
694
c
 
695
c     functional derivatives 
 
696
c
 
697
            ds2dRho = -(F8/F3) * s2/rhoo
 
698
            ds2dG = s2/Gamma2
 
699
C
 
700
            dfx1ds2 = Mu*(F1/(Deno*Deno)) 
 
701
            dfx1dRho = dfx1ds2*ds2dRho
 
702
            dfx1dG = dfx1ds2*ds2dG
 
703
C
 
704
            dfx2ds2 = Mus*DExp(-Mus*s2/kapas)
 
705
            dfx2dRho = dfx2ds2*ds2dRho
 
706
            dfx2dG = dfx2ds2*ds2dG
 
707
c
 
708
            dF1dW = (at01 + F2*at02*W1 + F3*at03*W2
 
709
     $                    + F4*at04*W3 + F5*at05*W4
 
710
     $                    + F6*at06*W5 + F7*at07*W6
 
711
     $                    + F8*at08*W7 + F9*at09*W8
 
712
     $                    + F10*at10*W9+F11*at11*W10)
 
713
            dF2dW = (bt01 + F2*bt02*W1 + F3*bt03*W2
 
714
     $                    + F4*bt04*W3 + F5*bt05*W4
 
715
     $                    + F6*bt06*W5 + F7*bt07*W6
 
716
     $                    + F8*bt08*W7 + F9*bt09*W8
 
717
     $                    + F10*Bt10*W9+F11*Bt11*W10)
 
718
            dF3dW = (ct01 + F2*ct02*W1 + F3*ct03*W2
 
719
     $                    + F4*ct04*W3 + F5*ct05*W4
 
720
     $                    + F6*ct06*W5 + F7*ct07*W6
 
721
     $                    + F8*ct08*W7 + F9*ct09*W8
 
722
     $                    + F10*ct10*W9+F11*ct11*W10)
 
723
            dF4dW = (dt01 + F2*dt02*W1 + F3*dt03*W2
 
724
     $                    + F4*dt04*W3 + F5*dt05*W4
 
725
     $                    + F6*dt06*W5 + F7*dt07*W6
 
726
     $                    + F8*dt08*W7 + F9*dt09*W8
 
727
     $                    + F10*dt10*W9+F11*dt11*W10)
 
728
 
 
729
            dWdT = F2/((F1 + Tsig)**F2)
 
730
            dTdR = ((F6*PI*PI)**F2o3)*(rhoo**F2o3)/tauN
 
731
            dTdTau = -TauUEG/tauN**F2
 
732
C
 
733
           If(UseLC) then
 
734
             dElSRdR = PDUM
 
735
             dElLRdR = Ax*F4o3*Rho13-PDUM
 
736
           else
 
737
             dElSRdR=Ax*F4o3*Rho13
 
738
             dElLRdR=F0
 
739
           endIf  
 
740
           dGGA1dR = dElSRdR*fx1 + ElSR*dfx1dRho
 
741
           dGGA2dR = dElSRdR*fx2 + ElSR*dfx2dRho 
 
742
           dGGA3dR = dElLRdR*fx1 + ElLR*dfx1dRho
 
743
           dGGA4dR = dElLRdR*fx2 + ElLR*dfx2dRho 
 
744
c
 
745
           dF1dR = dF1dW*dWdT*dTdR
 
746
           dF1dTau=dF1dW*dWdT*dTdTau
 
747
           dF2dR = dF2dW*dWdT*dTdR
 
748
           dF2dTau=dF2dW*dWdT*dTdTau
 
749
           dF3dR = dF3dW*dWdT*dTdR
 
750
           dF3dTau=dF3dW*dWdT*dTdTau
 
751
           dF4dR = dF4dW*dWdT*dTdR
 
752
           dF4dTau=dF4dW*dWdT*dTdTau
 
753
c
 
754
           dGGA1dG = ElSR*dfx1dG
 
755
           dGGA2dG = ElSR*dfx2dG
 
756
           dGGA3dG = ElLR*dfx1dG
 
757
           dGGA4dG = ElLR*dfx2dG
 
758
c
 
759
           Amat(n,2) = Amat(n,2)   +dGGA1dR*Fsig1 + GGA1*dF1dR
 
760
     $                             +dGGA2dR*Fsig2 + GGA2*dF2dR
 
761
     $                             +dGGA3dR*Fsig3 + GGA3*dF3dR
 
762
     $                             +dGGA4dR*Fsig4 + GGA4*dF4dR
 
763
           Cmat(n,3)=  Cmat(n,3)   +dGGA1dG*Fsig1 + dGGA2dG*Fsig2
 
764
     $                             +dGGA3dG*Fsig3 + dGGA4dG*Fsig4
 
765
           Mmat(n,2)=  Mmat(n,2)   +GGA1*dF1dTau  + GGA2*dF2dTau
 
766
     $                             +GGA3*dF3dTau  + GGA4*dF4dTau
 
767
c
 
768
20      continue
 
769
      endif
 
770
      return
 
771
      end
 
772
c
 
773
      Subroutine xc_xm11_d2()
 
774
      call errquit(' not coded ',0,0)
 
775
      return
 
776
      end
 
777
c
 
778
      SUBROUTINE LRCLSDA(Emu,Rho,F,D1F)
 
779
c
 
780
c***********************************************
 
781
c                                               
 
782
c   INPUT:                                      
 
783
c      Emu - Value of mu (or omega)
 
784
c      Rho - Spin density                 
 
785
c                                               
 
786
c   OUTPUT:                                     
 
787
c      F      - Functional value               
 
788
c      D1F    - First derivative               
 
789
c                                               
 
790
c***********************************************
 
791
c
 
792
      IMPLICIT REAL*8 (a-h,o-z)
 
793
      Save F1, F2, F3, F4, F5, F6, F7, F8, F9
 
794
      DATA F1/1.0D+00/,F2/2.0D+00/,F3/3.0D+00/,F4/4.0D+00/,F5/5.0D+00/,
 
795
     $     F6/6.0D+00/,F7/7.0D+00/,F8/8.0D+00/,F9/9.0D+00/
 
796
C
 
797
      PARAMETER( PI = 3.1415926535897932384626433832795D+00 )
 
798
C
 
799
      F1o2 = F1 / F2
 
800
      F1o3 = F1 / F3
 
801
      F1o4 = F1 / F4
 
802
      F4o3 = F4 / F3
 
803
      F8o3 = F8 / F3
 
804
      PI12 = SQRT(Pi)
 
805
C
 
806
      AX   = -(F3/F2) * (F4o3*PI)**(-F1o3)
 
807
      Cmu  = (F6*Pi**F2)**F1o3   
 
808
C
 
809
      Rho13 = Rho**F1o3
 
810
      Rho43 = Rho**F4o3
 
811
c
 
812
      tmu  = Emu/(F2*Cmu*Rho13)
 
813
      tmu2 = tmu*tmu
 
814
      tmu3 = tmu*tmu2
 
815
c
 
816
      W    = DExp(-F1o4/tmu2)
 
817
      ERFV = DErf( F1o2/tmu)
 
818
      dtmudR = -F1o3*tmu / Rho
 
819
c
 
820
      Fsr = F1-F4o3*tmu*(-F6*tmu+F8*tmu3+W*
 
821
     $        (F4*tmu-F8*tmu3)+F2*PI12*ERFV)
 
822
      dFsrdtmu = F8o3*(F2*tmu*(F3-F8*tmu2+W*
 
823
     $          (-F1+F8*tmu2))-PI12*ERFV)
 
824
c
 
825
      F = Ax*Rho43*Fsr
 
826
      D1F = Ax*F4o3*Rho13*Fsr + Ax*Rho43*(dFsrdtmu*dtmudR)
 
827
c
 
828
      RETURN
 
829
      END
 
830