← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwxc/nwxc_x_m06.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_m06.F
 
1
C> \ingroup nwxc
 
2
C> @{
 
3
C>
 
4
C> \file nwxc_x_m06.F
 
5
C> Implementation of the M06 exchange functional
 
6
C>
 
7
C> @}
 
8
C>
 
9
C> \ingroup nwxc_priv
 
10
C> @{
 
11
C>
 
12
C> \brief The M06 exchange functional
 
13
C>
 
14
C> The M06 functional [1,2] is a meta-GGA of which this evaluates
 
15
C> the exchange component.
 
16
C>
 
17
C> ### References ###
 
18
C>
 
19
C> [1] Y Zhao, DG Truhlar,
 
20
C> "A new local density functional for main-group thermochemistry,
 
21
C> transition metal bonding, thermochemical kinetics, and noncovalent
 
22
C> interactions",
 
23
C> J. Chem. Phys. <b>125</b>, 194101 (2006), DOI:
 
24
C> <a href="http://dx.doi.org/10.1063/1.2370993">
 
25
C> 10.1063/1.2370993</a>.
 
26
C>
 
27
C> [2] Y Zhao, DG Truhlar,
 
28
C> "Density functional for spectroscopy: No long-range self-interaction
 
29
C> error, good performance for Rydberg and charge-transfer states,
 
30
C> and better performance on average than B3LYP for ground states",
 
31
C> J. Phys. Chem. A <b>110</b>, 13126-13130 (2006), DOI:
 
32
C> <a href="http://dx.doi.org/10.1021/jp066479k">
 
33
C> 10.1021/jp066479k</a>.
 
34
C>
 
35
c   M06 suite  exchange functional  
 
36
c           META GGA
 
37
C         utilizes ingredients:
 
38
c                              rho   -  density
 
39
c                              delrho - gradient of density
 
40
c                              tau - K.S kinetic energy density
 
41
c                              tauU - uniform-gas KE density
 
42
c                              ijzy - 1  M06-L  
 
43
c                              ijzy - 2  M06-HF
 
44
c                              ijzy - 3  M06
 
45
c                              ijzy - 4  M06-2X
 
46
c     References: 
 
47
c     [a]       Zhao, Y. and  Truhlar, D. G. J. Chem. Phys. 2006, 125, 194101;
 
48
c     [b]       Zhao, Y. and  Truhlar, D. G. J. Phys. Chem. A (2006),110(49),13126-13130.    
 
49
 
 
50
 
 
51
      Subroutine nwxc_x_m06(param, tol_rho, ipol, nq, wght, rho, rgamma,
 
52
     &                      tau, func, Amat, Cmat, Mmat)
 
53
c   
 
54
c$Id: xc_xm06.F 21740 2012-01-11 00:25:15Z edo $
 
55
c
 
56
      implicit none
 
57
c
 
58
#include "nwxc_param.fh"
 
59
c
 
60
      double precision param(*) !< [Input] Parameters of the functional
 
61
                                !< - param( 1): \f$ d_0 \f$
 
62
                                !< - param( 2): \f$ d_1 \f$
 
63
                                !< - param( 3): \f$ d_2 \f$
 
64
                                !< - param( 4): \f$ d_3 \f$
 
65
                                !< - param( 5): \f$ d_4 \f$
 
66
                                !< - param( 6): \f$ d_5 \f$
 
67
                                !< - param( 7): \f$ a_0 \f$
 
68
                                !< - param( 8): \f$ a_1 \f$
 
69
                                !< - param( 9): \f$ a_2 \f$
 
70
                                !< - param(10): \f$ a_3 \f$
 
71
                                !< - param(11): \f$ a_4 \f$
 
72
                                !< - param(12): \f$ a_5 \f$
 
73
                                !< - param(13): \f$ a_6 \f$
 
74
                                !< - param(14): \f$ a_7 \f$
 
75
                                !< - param(15): \f$ a_8 \f$
 
76
                                !< - param(16): \f$ a_9 \f$
 
77
                                !< - param(17): \f$ a_10 \f$
 
78
                                !< - param(18): \f$ a_11 \f$
 
79
      double precision tol_rho !< [Input] The lower limit on the density
 
80
      integer nq               !< [Input] The number of points
 
81
      integer ipol             !< [Input] The number of spin channels
 
82
      double precision wght    !< [Input] The weight of the functional
 
83
c
 
84
c     Charge Density
 
85
c
 
86
      double precision rho(nq,*) !< [Input] The density
 
87
c
 
88
c     Charge Density Gradient Norm
 
89
c
 
90
      double precision rgamma(nq,*) !< [Input] The density gradient norm
 
91
c
 
92
c     Kinetic Energy Density
 
93
c
 
94
      double precision tau(nq,*) !< [Input] The kinetic energy density
 
95
c
 
96
c     Functional values
 
97
c
 
98
      double precision func(*) !< [Output] The functional value
 
99
c
 
100
c     Sampling Matrices for the XC Potential
 
101
c
 
102
      double precision Amat(nq,*) !< [Output] Derivative wrt density
 
103
      double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma
 
104
      double precision Mmat(nq,*) !< [Output] Derivative wrt tau
 
105
c      
 
106
      double precision pi
 
107
c
 
108
      integer n
 
109
      double precision at1, at2, at3, at4, at5, at6, at7, at8, at9
 
110
      double precision at, at10, at11, at0, C1, C2, fL, fNL
 
111
      double precision rrho, rho43, rho13, rhoo, rho53
 
112
      double precision Gamma2, Gamma
 
113
      double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6
 
114
      double precision W7, W8, W9, W10, W11, Fsig
 
115
c
 
116
      double precision tauN,tauu,DTol
 
117
      double precision F83, F23, F53, F43, F13, F1o3
 
118
      double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5
 
119
      double precision One, Two, Three, Four, Five, Six, Seven, Eight
 
120
      double precision Nine, F10, F11
 
121
      double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd
 
122
 
 
123
c      functional derivatives below FFFFFFFFFFFF
 
124
 
 
125
      double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR
 
126
      double precision dFdTau, dGGAdG,tauW
 
127
 
 
128
c     functional derivatives above FFFFFFFFFFFF
 
129
 
 
130
 
 
131
cedo       parameter( pi = 3.1415926535897932384626433832795d0 )
 
132
         
 
133
      
 
134
      parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 
 
135
     &     F3o2=3.d0/2.d0,F13=1.d0/3.d0)
 
136
      parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0)
 
137
      parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0)
 
138
      parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 
 
139
     &     Five=5.0d0,Six=6.0d0, Seven=7.0d0,
 
140
     &     Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0)
 
141
      pi=acos(-1d0)      
 
142
 
 
143
      at0  = param(7)
 
144
      at1  = param(8)
 
145
      at2  = param(9)
 
146
      at3  = param(10)
 
147
      at4  = param(11)
 
148
      at5  = param(12)
 
149
      at6  = param(13)
 
150
      at7  = param(14)
 
151
      at8  = param(15)
 
152
      at9  = param(16)
 
153
      at10 = param(17)
 
154
      at11 = param(18)
 
155
c     if (ijzy.eq.1) then
 
156
c       at0=    3.987756D-01
 
157
c       at1=    2.548219D-01
 
158
c       at2=    3.923994D-01
 
159
c       at3=    -2.103655D+00
 
160
c       at4=    -6.302147D+00
 
161
c       at5=    1.097615D+01
 
162
c       at6=    3.097273D+01
 
163
c       at7=    -2.318489D+01
 
164
c       at8=    -5.673480D+01
 
165
c       at9=    2.160364D+01
 
166
c       at10=   3.421814D+01
 
167
c       at11=   -9.049762D+00
 
168
c      elseif (ijzy.eq.2) then
 
169
c     Parameters for M06-HF
 
170
c       at0=    1.179732D-01
 
171
c       at1=    -1.066708D+00
 
172
c       at2=    -1.462405D-01
 
173
c       at3=    7.481848D+00
 
174
c       at4=    3.776679D+00
 
175
c       at5=    -4.436118D+01
 
176
c       at6=    -1.830962D+01
 
177
c       at7=    1.003903D+02
 
178
c       at8=    3.864360D+01
 
179
c       at9=    -9.806018D+01
 
180
c       at10=   -2.557716D+01
 
181
c       at11=   3.590404D+01
 
182
c      elseif (ijzy.eq.3) then
 
183
c     Parameters for M06
 
184
c       at0=    5.877943D-01
 
185
c       at1=    -1.371776D-01
 
186
c       at2=    2.682367D-01
 
187
c       at3=    -2.515898D+00
 
188
c       at4=    -2.978892D+00
 
189
c       at5=    8.710679D+00
 
190
c       at6=    1.688195D+01
 
191
c       at7=    -4.489724D+00
 
192
c       at8=    -3.299983D+01
 
193
c       at9=    -1.449050D+01
 
194
c       at10=   2.043747D+01
 
195
c       at11=   1.256504D+01
 
196
c      elseif (ijzy.eq.4) then
 
197
c     Parameters for M06-2X
 
198
c       at0=    4.600000D-01
 
199
c       at1=    -2.206052D-01
 
200
c       at2=    -9.431788D-02
 
201
c       at3=    2.164494D+00
 
202
c       at4=    -2.556466D+00
 
203
c       at5=    -1.422133D+01
 
204
c       at6=    1.555044D+01
 
205
c       at7=    3.598078D+01
 
206
c       at8=    -2.722754D+01
 
207
c       at9=    -3.924093D+01
 
208
c       at10=   1.522808D+01
 
209
c       at11=   1.522227D+01
 
210
c     endif
 
211
 
 
212
      at=1.0d0
 
213
      call nwxc_x_vs98(param,tol_rho, ipol, nq, wght, rho, rgamma, tau,
 
214
     &                 func, Amat, Cmat, Mmat)
 
215
 
 
216
 
 
217
      C1 = 0.2195149727645171d0
 
218
      C2 = C1/0.804d0 
 
219
cedo      DTol=1.0D-8
 
220
      DTol=tol_rho
 
221
C
 
222
C     Scale factors for local and non-local contributions.
 
223
C
 
224
      fL  =  wght
 
225
      fNL =  wght
 
226
      Cs = 0.5d0/(3.0d0*pi*pi)**F13
 
227
      P32 = (3.d0*pi**2)**F23
 
228
         
 
229
c     
 
230
       Ax = (-0.75d0)*(3.0d0/pi)**F13
 
231
 
 
232
 
 
233
c
 
234
      if (ipol.eq.1 )then
 
235
c
 
236
c        ======> SPIN-RESTRICTED <======
 
237
c                     or
 
238
c                SPIN-UNPOLARIZED
 
239
c
 
240
c
 
241
         do 10 n = 1, nq
 
242
            if (rho(n,R_T).lt.DTol) goto 10
 
243
           
 
244
            rhoo = rho(n,R_T)
 
245
            rho43 = rhoo**F4o3  
 
246
            rrho = 1d0/rhoo       ! reciprocal of rho
 
247
            rho13 = rho43*rrho
 
248
            rho53 = rhoo**F53
 
249
 
 
250
c
 
251
             
 
252
            tauN = tau(n,T_T)
 
253
         if(taun.lt.dtol) goto 10
 
254
            tauu=tauN 
 
255
            TauUEG=0.3d0*P32*rho53
 
256
            Tsig =TauUEG/tauu
 
257
            Wsig =(Tsig-One)/(Tsig+One)
 
258
            W1=Wsig 
 
259
            W2=Wsig*W1
 
260
            W3=Wsig*W2
 
261
            W4=Wsig*W3
 
262
            W5=Wsig*W4
 
263
            W6=Wsig*W5
 
264
            W7=Wsig*W6
 
265
            W8=Wsig*W7
 
266
            W9=Wsig*W8
 
267
            W10=Wsig*W9
 
268
            W11=Wsig*W10
 
269
            Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3
 
270
     &          + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
271
     &          + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
272
 
 
273
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
274
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
275
c    &              delrho(n,3,1)*delrho(n,3,1)
 
276
            Gamma2 = rgamma(n,G_TT)
 
277
            Gamma = dsqrt(Gamma2)
 
278
         if(gamma.lt.dtol) goto 10
 
279
            s      = Cs*Gamma/rho43
 
280
            s2     = s*s
 
281
            En     = C1*s2
 
282
            Ed     = One + C2*s2
 
283
            E      = En/Ed
 
284
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*Fsig*qwght(n)
 
285
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*Fsig
 
286
c
 
287
c     functional derivatives 
 
288
c
 
289
            dEn   = Two*C1*s
 
290
            dEd   = Two*C2*s
 
291
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
292
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
293
     &             + Four*at4*W3 + Five*at5*W4
 
294
     &             + Six*at6*W5 + Seven*at7*W6
 
295
     &             + Eight*at8*W7 + Nine*at9*W8
 
296
     &             + F10*at10*W9 + F11*at11*W10)
 
297
            dWdT = Two/((One + Tsig)**2)
 
298
            dTdR = (0.5d0*P32*rho13*rho13)/tauu
 
299
            dTdTau = -TauUEG/tauu**2
 
300
            dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
301
            dFdR = dFdW*dWdT*dTdR
 
302
            dFdTau=dFdW*dWdT*dTdTau
 
303
            dGGAdG =(fNL*dE*s/(Two*Gamma2))
 
304
            Amat(n,D1_RA)  = Amat(n,D1_RA) + dGGAdR*Fsig
 
305
     &                     + (fL+fNL*E)*Ax*rho43*dFdR
 
306
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA)
 
307
     &                     + Two*dGGAdG*Ax*rho43*Fsig 
 
308
            Mmat(n,D1_TA)  = Mmat(n,D1_TA)
 
309
     &                     + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
310
 
 
311
10      continue
 
312
 
 
313
 
 
314
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
315
      else  ! ipol=2
 
316
c
 
317
c        ======> SPIN-UNRESTRICTED <======
 
318
 
 
319
c
 
320
c  use spin density functional theory ie n-->2n
 
321
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
322
 
 
323
         do 20 n = 1, nq
 
324
           if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
 
325
c
 
326
c     Alpha            ALPHA               ALPHA
 
327
c
 
328
            if (rho(n,R_A).lt.0.5d0*DTol) goto 25
 
329
             rhoo = Two*rho(n,R_A)
 
330
             rho43 = rhoo**F4o3
 
331
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
332
             rho13 = rho43*rrho
 
333
             rho53 = rhoo**F53
 
334
c
 
335
             tauN = tau(n,T_A)
 
336
             tauu = Two*tauN
 
337
             TauUEG=0.3d0*P32*rho53
 
338
             Tsig =TauUEG/tauu
 
339
             Wsig =(Tsig-One)/(Tsig+One)
 
340
             W1=Wsig
 
341
             W2=Wsig*W1
 
342
             W3=Wsig*W2
 
343
             W4=Wsig*W3
 
344
             W5=Wsig*W4
 
345
             W6=Wsig*W5
 
346
             W7=Wsig*W6
 
347
             W8=Wsig*W7
 
348
             W9=Wsig*W8
 
349
             W10=Wsig*W9
 
350
             W11=Wsig*W10
 
351
             Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3
 
352
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
353
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
354
 
 
355
 
 
356
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
357
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
358
c    &              delrho(n,3,1)*delrho(n,3,1)
 
359
            Gamma2 = rgamma(n,G_AA)
 
360
            Gamma2 = Four*Gamma2
 
361
            Gamma = dsqrt(Gamma2)
 
362
 
 
363
            s      = Cs*Gamma/rho43
 
364
            s2     = s*s
 
365
            En     = C1*s2
 
366
            Ed     = One + C2*s2
 
367
            E      = En/Ed
 
368
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n)
 
369
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0
 
370
c
 
371
c     functional derivatives
 
372
c
 
373
            dEn   = Two*C1*s
 
374
            dEd   = Two*C2*s
 
375
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
376
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
377
     &             + Four*at4*W3 + Five*at5*W4
 
378
     &             + Six*at6*W5 + Seven*at7*W6
 
379
     &             + Eight*at8*W7 + Nine*at9*W8
 
380
     &             + F10*at10*W9 + F11*at11*W10)
 
381
            dWdT = Two/((One + Tsig)**2)
 
382
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
383
            dTdTau = -Two*TauUEG/tauu**2
 
384
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
385
            dFdR = dFdW*dWdT*dTdR
 
386
            dFdTau=dFdW*dWdT*dTdTau
 
387
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
388
 
 
389
            Amat(n,D1_RA)  = Amat(n,D1_RA) + (dGGAdR*(Fsig)
 
390
     &                     + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
391
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA)
 
392
     &                     + dGGAdG*Ax*rho43*(Fsig)*0.5d0
 
393
            Mmat(n,D1_TA)  = Mmat(n,D1_TA)
 
394
     &                     + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
395
 
 
396
c             write (*,*) "Ex,Amat(n,1),Cmat(n,1)",
 
397
c     &        Ex,Amat(n,1),Cmat(n,1)
 
398
 
 
399
c
 
400
c     Beta               BETA           BETA
 
401
c
 
402
 
 
403
25         continue
 
404
 
 
405
c
 
406
c     Beta
 
407
c
 
408
            if (rho(n,R_B).lt.0.5d0*DTol) goto 20
 
409
             rhoo = Two*rho(n,R_B)
 
410
             rho43 = rhoo**F4o3
 
411
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
412
             rho13 = rho43*rrho
 
413
             rho53 = rhoo**F53
 
414
c
 
415
             tauN = tau(n,T_B)
 
416
             tauu = Two*tauN
 
417
             TauUEG=0.3d0*P32*rho53
 
418
             Tsig =TauUEG/tauu
 
419
             Wsig =(Tsig-One)/(Tsig+One)
 
420
             W1=Wsig
 
421
             W2=Wsig*W1
 
422
             W3=Wsig*W2
 
423
             W4=Wsig*W3
 
424
             W5=Wsig*W4
 
425
             W6=Wsig*W5
 
426
             W7=Wsig*W6
 
427
             W8=Wsig*W7
 
428
             W9=Wsig*W8
 
429
             W10=Wsig*W9
 
430
             W11=Wsig*W10
 
431
             Fsig =at*(at0+at1*W1+ at2*W2 + at3*W3
 
432
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
433
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
434
 
 
435
 
 
436
c           Gamma2 = delrho(n,1,2)*delrho(n,1,2) +
 
437
c    &              delrho(n,2,2)*delrho(n,2,2) +
 
438
c    &              delrho(n,3,2)*delrho(n,3,2)
 
439
            Gamma2 = rgamma(n,G_BB)
 
440
            Gamma2 = Four*Gamma2
 
441
            Gamma = dsqrt(Gamma2)
 
442
            s      = Cs*Gamma/rho43
 
443
            s2     = s*s
 
444
            En     = C1*s2
 
445
            Ed     = One + C2*s2
 
446
            E      = En/Ed
 
447
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n)
 
448
            func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0
 
449
c
 
450
c     functional derivatives
 
451
c
 
452
            dEn   = Two*C1*s
 
453
            dEd   = Two*C2*s
 
454
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
455
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
456
     &             + Four*at4*W3 + Five*at5*W4
 
457
     &             + Six*at6*W5 + Seven*at7*W6
 
458
     &             + Eight*at8*W7 + Nine*at9*W8
 
459
     &             + F10*at10*W9 + F11*at11*W10)
 
460
            dWdT = Two/((One + Tsig)**2)
 
461
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
462
            dTdTau = -Two*TauUEG/tauu**2
 
463
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
464
            dFdR = dFdW*dWdT*dTdR
 
465
            dFdTau=dFdW*dWdT*dTdTau
 
466
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
467
 
 
468
            Amat(n,D1_RB)  = Amat(n,D1_RB) + (dGGAdR*(Fsig)
 
469
     &                     + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
470
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB)
 
471
     &                     + dGGAdG*Ax*rho43*(Fsig)*0.5d0
 
472
            Mmat(n,D1_TB)  = Mmat(n,D1_TB)
 
473
     &                     + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
474
c
 
475
20      continue
 
476
      endif
 
477
c
 
478
      return
 
479
      end
 
480
C>
 
481
C> \brief Evaluate the M06-2X exchange functional
 
482
C>
 
483
C> This routine evaluates the M06-2X exchange functional [1,2]. This functional
 
484
C> is closely related to the M06, M06-HF and M06-L exchange functionals
 
485
C> except that it does not use the VS98 exchange functional, hence the
 
486
C> parameter list is defined differently.
 
487
C>
 
488
C> ### References ###
 
489
C>
 
490
C> [1] Y Zhao, DG Truhlar,
 
491
C> "A new local density functional for main-group thermochemistry,
 
492
C> transition metal bonding, thermochemical kinetics, and noncovalent
 
493
C> interactions",
 
494
C> J. Chem. Phys. <b>125</b>, 194101 (2006), DOI:
 
495
C> <a href="http://dx.doi.org/10.1063/1.2370993">
 
496
C> 10.1063/1.2370993</a>.
 
497
C>
 
498
C> [2] Y Zhao, DG Truhlar,
 
499
C> "Density functional for spectroscopy: No long-range self-interaction
 
500
C> error, good performance for Rydberg and charge-transfer states,
 
501
C> and better performance on average than B3LYP for ground states",
 
502
C> J. Phys. Chem. A <b>110</b>, 13126-13130 (2006), DOI:
 
503
C> <a href="http://dx.doi.org/10.1021/jp066479k">
 
504
C> 10.1021/jp066479k</a>.
 
505
C>
 
506
 
 
507
      Subroutine nwxc_x_m06_2x(param, tol_rho, ipol, nq, wght, 
 
508
     &                         rho, rgamma, tau, func, Amat, Cmat, Mmat)
 
509
c   
 
510
c$Id: xc_xm06.F 21740 2012-01-11 00:25:15Z edo $
 
511
c
 
512
      implicit none
 
513
c
 
514
#include "nwxc_param.fh"
 
515
c
 
516
      double precision param(*) !< [Input] Parameters of the functional
 
517
                                !< - param( 1): \f$ a_0 \f$
 
518
                                !< - param( 2): \f$ a_1 \f$
 
519
                                !< - param( 3): \f$ a_2 \f$
 
520
                                !< - param( 4): \f$ a_3 \f$
 
521
                                !< - param( 5): \f$ a_4 \f$
 
522
                                !< - param( 6): \f$ a_5 \f$
 
523
                                !< - param( 7): \f$ a_6 \f$
 
524
                                !< - param( 8): \f$ a_7 \f$
 
525
                                !< - param( 9): \f$ a_8 \f$
 
526
                                !< - param(10): \f$ a_9 \f$
 
527
                                !< - param(11): \f$ a_10 \f$
 
528
                                !< - param(12): \f$ a_11 \f$
 
529
      double precision tol_rho !< [Input] The lower limit on the density
 
530
      integer nq               !< [Input] The number of points
 
531
      integer ipol             !< [Input] The number of spin channels
 
532
      double precision wght    !< [Input] The weight of the functional
 
533
c
 
534
c     Charge Density
 
535
c
 
536
      double precision rho(nq,*) !< [Input] The density
 
537
c
 
538
c     Charge Density Gradient Norm
 
539
c
 
540
      double precision rgamma(nq,*) !< [Input] The density gradient norm
 
541
c
 
542
c     Kinetic Energy Density
 
543
c
 
544
      double precision tau(nq,*) !< [Input] The kinetic energy density
 
545
c
 
546
c     Functional values
 
547
c
 
548
      double precision func(*) !< [Output] The functional value
 
549
c
 
550
c     Sampling Matrices for the XC Potential
 
551
c
 
552
      double precision Amat(nq,*) !< [Output] Derivative wrt density
 
553
      double precision Cmat(nq,*) !< [Output] Derivative wrt rgamma
 
554
      double precision Mmat(nq,*) !< [Output] Derivative wrt tau
 
555
c      
 
556
      double precision pi
 
557
c
 
558
      integer n
 
559
      double precision at1, at2, at3, at4, at5, at6, at7, at8, at9
 
560
      double precision at, at10, at11, at0, C1, C2, fL, fNL
 
561
      double precision rrho, rho43, rho13, rhoo, rho53
 
562
      double precision Gamma2, Gamma
 
563
      double precision TauUEG, Tsig, Wsig, W1, W2, W3, W4, W5, W6
 
564
      double precision W7, W8, W9, W10, W11, Fsig
 
565
c
 
566
      double precision tauN,tauu,DTol
 
567
      double precision F83, F23, F53, F43, F13, F1o3
 
568
      double precision F1o4, F2o3, F3o2, F4o3, F4o9, F3o5
 
569
      double precision One, Two, Three, Four, Five, Six, Seven, Eight
 
570
      double precision Nine, F10, F11
 
571
      double precision Cs, Ax, P32, s, s2, En, Ed, E, dE, dEn, dEd
 
572
 
 
573
c      functional derivatives below FFFFFFFFFFFF
 
574
 
 
575
      double precision dFdW, dWdT, dTdR, dTdTau, dGGAdR, dFdR
 
576
      double precision dFdTau, dGGAdG,tauW
 
577
 
 
578
c     functional derivatives above FFFFFFFFFFFF
 
579
 
 
580
 
 
581
cedo       parameter( pi = 3.1415926535897932384626433832795d0 )
 
582
         
 
583
      
 
584
      parameter (F1o3=1.d0/3.d0, F1o4=1.d0/4.d0, F2o3=2.d0/3.d0, 
 
585
     &     F3o2=3.d0/2.d0,F13=1.d0/3.d0)
 
586
      parameter (F4o3=4.d0/3.d0, F4o9=4.d0/9.d0, F3o5=3.d0/5.d0)
 
587
      parameter (F83=8.d0/3.0d0, F23=2.0d0/3.d0, F53=5.d0/3.d0)
 
588
      parameter (One=1.0d0, Two=2.0d0, Three=3.0d0, Four=4.0d0, 
 
589
     &     Five=5.0d0,Six=6.0d0, Seven=7.0d0,
 
590
     &     Eight=8.0d0, Nine=9.0d0,F10=10.d0, F11=11.d0)
 
591
      pi=acos(-1d0)      
 
592
 
 
593
      at0  = param(1)
 
594
      at1  = param(2)
 
595
      at2  = param(3)
 
596
      at3  = param(4)
 
597
      at4  = param(5)
 
598
      at5  = param(6)
 
599
      at6  = param(7)
 
600
      at7  = param(8)
 
601
      at8  = param(9)
 
602
      at9  = param(10)
 
603
      at10 = param(11)
 
604
      at11 = param(12)
 
605
c     if (ijzy.eq.1) then
 
606
c       at0=    3.987756D-01
 
607
c       at1=    2.548219D-01
 
608
c       at2=    3.923994D-01
 
609
c       at3=    -2.103655D+00
 
610
c       at4=    -6.302147D+00
 
611
c       at5=    1.097615D+01
 
612
c       at6=    3.097273D+01
 
613
c       at7=    -2.318489D+01
 
614
c       at8=    -5.673480D+01
 
615
c       at9=    2.160364D+01
 
616
c       at10=   3.421814D+01
 
617
c       at11=   -9.049762D+00
 
618
c      elseif (ijzy.eq.2) then
 
619
c     Parameters for M06-HF
 
620
c       at0=    1.179732D-01
 
621
c       at1=    -1.066708D+00
 
622
c       at2=    -1.462405D-01
 
623
c       at3=    7.481848D+00
 
624
c       at4=    3.776679D+00
 
625
c       at5=    -4.436118D+01
 
626
c       at6=    -1.830962D+01
 
627
c       at7=    1.003903D+02
 
628
c       at8=    3.864360D+01
 
629
c       at9=    -9.806018D+01
 
630
c       at10=   -2.557716D+01
 
631
c       at11=   3.590404D+01
 
632
c      elseif (ijzy.eq.3) then
 
633
c     Parameters for M06
 
634
c       at0=    5.877943D-01
 
635
c       at1=    -1.371776D-01
 
636
c       at2=    2.682367D-01
 
637
c       at3=    -2.515898D+00
 
638
c       at4=    -2.978892D+00
 
639
c       at5=    8.710679D+00
 
640
c       at6=    1.688195D+01
 
641
c       at7=    -4.489724D+00
 
642
c       at8=    -3.299983D+01
 
643
c       at9=    -1.449050D+01
 
644
c       at10=   2.043747D+01
 
645
c       at11=   1.256504D+01
 
646
c      elseif (ijzy.eq.4) then
 
647
c     Parameters for M06-2X
 
648
c       at0=    4.600000D-01
 
649
c       at1=    -2.206052D-01
 
650
c       at2=    -9.431788D-02
 
651
c       at3=    2.164494D+00
 
652
c       at4=    -2.556466D+00
 
653
c       at5=    -1.422133D+01
 
654
c       at6=    1.555044D+01
 
655
c       at7=    3.598078D+01
 
656
c       at8=    -2.722754D+01
 
657
c       at9=    -3.924093D+01
 
658
c       at10=   1.522808D+01
 
659
c       at11=   1.522227D+01
 
660
c     endif
 
661
 
 
662
      at=1.0d0
 
663
 
 
664
      C1 = 0.2195149727645171d0
 
665
      C2 = C1/0.804d0 
 
666
cedo      DTol=1.0D-8
 
667
      DTol=tol_rho
 
668
C
 
669
C     Scale factors for local and non-local contributions.
 
670
C
 
671
      fL  =  wght
 
672
      fNL =  wght
 
673
      Cs = 0.5d0/(3.0d0*pi*pi)**F13
 
674
      P32 = (3.d0*pi**2)**F23
 
675
         
 
676
c     
 
677
       Ax = (-0.75d0)*(3.0d0/pi)**F13
 
678
 
 
679
 
 
680
c
 
681
      if (ipol.eq.1 )then
 
682
c
 
683
c        ======> SPIN-RESTRICTED <======
 
684
c                     or
 
685
c                SPIN-UNPOLARIZED
 
686
c
 
687
c
 
688
         do 10 n = 1, nq
 
689
            if (rho(n,R_T).lt.DTol) goto 10
 
690
           
 
691
            rhoo = rho(n,R_T)
 
692
            rho43 = rhoo**F4o3  
 
693
            rrho = 1d0/rhoo       ! reciprocal of rho
 
694
            rho13 = rho43*rrho
 
695
            rho53 = rhoo**F53
 
696
 
 
697
c
 
698
             
 
699
            tauN = tau(n,T_T)
 
700
         if(taun.lt.dtol) goto 10
 
701
            tauu=tauN 
 
702
            TauUEG=0.3d0*P32*rho53
 
703
            Tsig =TauUEG/tauu
 
704
            Wsig =(Tsig-One)/(Tsig+One)
 
705
            W1=Wsig 
 
706
            W2=Wsig*W1
 
707
            W3=Wsig*W2
 
708
            W4=Wsig*W3
 
709
            W5=Wsig*W4
 
710
            W6=Wsig*W5
 
711
            W7=Wsig*W6
 
712
            W8=Wsig*W7
 
713
            W9=Wsig*W8
 
714
            W10=Wsig*W9
 
715
            W11=Wsig*W10
 
716
            Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3
 
717
     &          + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
718
     &          + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
719
 
 
720
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
721
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
722
c    &              delrho(n,3,1)*delrho(n,3,1)
 
723
            Gamma2 = rgamma(n,G_TT)
 
724
            Gamma = dsqrt(Gamma2)
 
725
         if(gamma.lt.dtol) goto 10
 
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)*Fsig*qwght(n)
 
732
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*Fsig
 
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 = (0.5d0*P32*rho13*rho13)/tauu
 
746
            dTdTau = -TauUEG/tauu**2
 
747
            dGGAdR = F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
748
            dFdR = dFdW*dWdT*dTdR
 
749
            dFdTau=dFdW*dWdT*dTdTau
 
750
            dGGAdG =(fNL*dE*s/(Two*Gamma2))
 
751
            Amat(n,D1_RA)  = Amat(n,D1_RA) + dGGAdR*Fsig
 
752
     &                     + (fL+fNL*E)*Ax*rho43*dFdR
 
753
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA)
 
754
     &                     + Two*dGGAdG*Ax*rho43*Fsig 
 
755
            Mmat(n,D1_TA)  = Mmat(n,D1_TA)
 
756
     &                     + 0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
757
 
 
758
10      continue
 
759
 
 
760
 
 
761
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
762
      else  ! ipol=2
 
763
c
 
764
c        ======> SPIN-UNRESTRICTED <======
 
765
 
 
766
c
 
767
c  use spin density functional theory ie n-->2n
 
768
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
769
 
 
770
         do 20 n = 1, nq
 
771
           if (rho(n,R_A)+rho(n,R_B).lt.DTol) goto 20
 
772
c
 
773
c     Alpha            ALPHA               ALPHA
 
774
c
 
775
            if (rho(n,R_A).lt.0.5d0*DTol) goto 25
 
776
             rhoo = Two*rho(n,R_A)
 
777
             rho43 = rhoo**F4o3
 
778
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
779
             rho13 = rho43*rrho
 
780
             rho53 = rhoo**F53
 
781
c
 
782
             tauN = tau(n,T_A)
 
783
             tauu = Two*tauN
 
784
             TauUEG=0.3d0*P32*rho53
 
785
             Tsig =TauUEG/tauu
 
786
             Wsig =(Tsig-One)/(Tsig+One)
 
787
             W1=Wsig
 
788
             W2=Wsig*W1
 
789
             W3=Wsig*W2
 
790
             W4=Wsig*W3
 
791
             W5=Wsig*W4
 
792
             W6=Wsig*W5
 
793
             W7=Wsig*W6
 
794
             W8=Wsig*W7
 
795
             W9=Wsig*W8
 
796
             W10=Wsig*W9
 
797
             W11=Wsig*W10
 
798
             Fsig =at*(at0 + at1*W1+ at2*W2 + at3*W3
 
799
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
800
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
801
 
 
802
 
 
803
c           Gamma2 = delrho(n,1,1)*delrho(n,1,1) +
 
804
c    &              delrho(n,2,1)*delrho(n,2,1) +
 
805
c    &              delrho(n,3,1)*delrho(n,3,1)
 
806
            Gamma2 = rgamma(n,G_AA)
 
807
            Gamma2 = Four*Gamma2
 
808
            Gamma = dsqrt(Gamma2)
 
809
 
 
810
            s      = Cs*Gamma/rho43
 
811
            s2     = s*s
 
812
            En     = C1*s2
 
813
            Ed     = One + C2*s2
 
814
            E      = En/Ed
 
815
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n)
 
816
            func(n)=func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0
 
817
c
 
818
c     functional derivatives
 
819
c
 
820
            dEn   = Two*C1*s
 
821
            dEd   = Two*C2*s
 
822
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
823
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
824
     &             + Four*at4*W3 + Five*at5*W4
 
825
     &             + Six*at6*W5 + Seven*at7*W6
 
826
     &             + Eight*at8*W7 + Nine*at9*W8
 
827
     &             + F10*at10*W9 + F11*at11*W10)
 
828
            dWdT = Two/((One + Tsig)**2)
 
829
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
830
            dTdTau = -Two*TauUEG/tauu**2
 
831
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
832
            dFdR = dFdW*dWdT*dTdR
 
833
            dFdTau=dFdW*dWdT*dTdTau
 
834
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
835
 
 
836
            Amat(n,D1_RA)  = Amat(n,D1_RA) + (dGGAdR*(Fsig)
 
837
     &                     + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
838
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA)
 
839
     &                     + dGGAdG*Ax*rho43*(Fsig)*0.5d0
 
840
            Mmat(n,D1_TA)  = Mmat(n,D1_TA)
 
841
     &                     + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
842
 
 
843
c             write (*,*) "Ex,Amat(n,1),Cmat(n,1)",
 
844
c     &        Ex,Amat(n,1),Cmat(n,1)
 
845
 
 
846
c
 
847
c     Beta               BETA           BETA
 
848
c
 
849
 
 
850
25         continue
 
851
 
 
852
c
 
853
c     Beta
 
854
c
 
855
            if (rho(n,R_B).lt.0.5d0*DTol) goto 20
 
856
             rhoo = Two*rho(n,R_B)
 
857
             rho43 = rhoo**F4o3
 
858
             rrho = 1.0d0/rhoo       ! reciprocal of rho
 
859
             rho13 = rho43*rrho
 
860
             rho53 = rhoo**F53
 
861
c
 
862
             tauN = tau(n,T_B)
 
863
             tauu = Two*tauN
 
864
             TauUEG=0.3d0*P32*rho53
 
865
             Tsig =TauUEG/tauu
 
866
             Wsig =(Tsig-One)/(Tsig+One)
 
867
             W1=Wsig
 
868
             W2=Wsig*W1
 
869
             W3=Wsig*W2
 
870
             W4=Wsig*W3
 
871
             W5=Wsig*W4
 
872
             W6=Wsig*W5
 
873
             W7=Wsig*W6
 
874
             W8=Wsig*W7
 
875
             W9=Wsig*W8
 
876
             W10=Wsig*W9
 
877
             W11=Wsig*W10
 
878
             Fsig =at*(at0+at1*W1+ at2*W2 + at3*W3
 
879
     &           + at4*W4 + at5*W5 + at6*W6 + at7*W7
 
880
     &           + at8*W8 + at9*W9 + at10*W10 + at11*W11)
 
881
 
 
882
 
 
883
c           Gamma2 = delrho(n,1,2)*delrho(n,1,2) +
 
884
c    &              delrho(n,2,2)*delrho(n,2,2) +
 
885
c    &              delrho(n,3,2)*delrho(n,3,2)
 
886
            Gamma2 = rgamma(n,G_BB)
 
887
            Gamma2 = Four*Gamma2
 
888
            Gamma = dsqrt(Gamma2)
 
889
            s      = Cs*Gamma/rho43
 
890
            s2     = s*s
 
891
            En     = C1*s2
 
892
            Ed     = One + C2*s2
 
893
            E      = En/Ed
 
894
c           Ex = Ex + rho43*Ax*(fL+fNL*E)*(Fsig)*0.5d0*qwght(n)
 
895
            func(n)= func(n)+rho43*Ax*(fL+fNL*E)*(Fsig)*.5d0
 
896
c
 
897
c     functional derivatives
 
898
c
 
899
            dEn   = Two*C1*s
 
900
            dEd   = Two*C2*s
 
901
            dE    = (dEn*Ed-En*dEd)/(Ed*Ed)
 
902
            dFdW = at*(at1 + Two*at2*W1 + Three*at3*W2
 
903
     &             + Four*at4*W3 + Five*at5*W4
 
904
     &             + Six*at6*W5 + Seven*at7*W6
 
905
     &             + Eight*at8*W7 + Nine*at9*W8
 
906
     &             + F10*at10*W9 + F11*at11*W10)
 
907
            dWdT = Two/((One + Tsig)**2)
 
908
            dTdR = Two*(0.5d0*P32*rho13*rho13)/tauu
 
909
            dTdTau = -Two*TauUEG/tauu**2
 
910
            dGGAdR = Two*F4o3*rho13*Ax*(fL+fNL*(E-s*dE))
 
911
            dFdR = dFdW*dWdT*dTdR
 
912
            dFdTau=dFdW*dWdT*dTdTau
 
913
            dGGAdG =Four*(fNL*dE*s/(Two*Gamma2))
 
914
 
 
915
            Amat(n,D1_RB)  = Amat(n,D1_RB) + (dGGAdR*(Fsig)
 
916
     &                     + (fL+fNL*E)*Ax*rho43*dFdR)*0.5d0
 
917
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB)
 
918
     &                     + dGGAdG*Ax*rho43*(Fsig)*0.5d0
 
919
            Mmat(n,D1_TB)  = Mmat(n,D1_TB)
 
920
     &                     + 0.5d0*0.5d0*rho43*Ax*(fL+fNL*E)*dFdTau
 
921
c
 
922
20      continue
 
923
      endif
 
924
c
 
925
      return
 
926
      end
 
927
 
 
928
 
 
929
 
 
930
 
 
931
      Subroutine nwxc_x_m06_d2()
 
932
      call errquit(' xm06: d2 not coded ',0,0)
 
933
      return
 
934
      end
 
935
C> @}