← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_bnl.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_bnl.F
 
1
c
 
2
c     -----------------------------------------------------------------------
 
3
c     Uniform electron gas exchange functional for the erfc(r)/r interaction
 
4
c     as implemented in the following paper:
 
5
c     "A well-tempered density functional theory of electrons in molecules"
 
6
c     Ester Livshits & Roi Baer, Phys. Chem. Chem. Phys., 9, 2932 (2007)
 
7
c     The other relevant publication is: 
 
8
c     R. Baer, D. Neuhauser, Phys. Rev. Lett., 94, 043002 (2005)
 
9
c     -----------------------------------------------------------------------
 
10
c
 
11
#ifndef SECOND_DERIV
 
12
      subroutine xc_bnl(tol_rho, fac, lfac, nlfac, rho, Amat, nq,
 
13
     &                    ipol, Ex, qwght, ldew, func)
 
14
#else
 
15
c     For locations of 2nd derivatives of functionals in array
 
16
#include "dft2drv.fh"
 
17
      subroutine xc_bnl_d2(tol_rho, fac, lfac, nlfac, rho, Amat,
 
18
     &                       Amat2, nq, ipol, Ex, qwght, ldew, func)
 
19
#endif
 
20
c
 
21
      implicit none
 
22
c
 
23
#include "errquit.fh"
 
24
#include "stdio.fh"
 
25
c
 
26
ccase...start
 
27
#include "case.fh"
 
28
ccase...end
 
29
c
 
30
      integer nq, ipol, n
 
31
      double precision fac, Ex, total
 
32
      logical ldew, lfac, nlfac
 
33
      double precision func(*) ! value of the functional [output]
 
34
      double precision tol_rho
 
35
      double precision rho(nq,(ipol*(ipol+1))/2) ! charge density
 
36
      double precision qwght(nq)                 ! quadrature weights
 
37
      double precision Amat(nq,ipol)             ! partial first derivatives  
 
38
      double precision F(nq),RA(nq),RB(nq)
 
39
      double precision rhoA, rhoB, rhoTotal, rhoA1, rhoB1
 
40
      double precision gamma
 
41
      double precision fA, fB, fpA, fpB, fppA, fppB
 
42
      double precision EpsX
 
43
      double precision EpsXprime
 
44
      double precision EpsTwoXprime
 
45
c
 
46
#ifdef SECOND_DERIV
 
47
      double precision Amat2(nq,NCOL_AMAT2)               ! partial second derivatives
 
48
#endif
 
49
c
 
50
c     -----------------------------------------------------------------------
 
51
c     Preliminaries
 
52
c     -----------------------------------------------------------------------
 
53
c
 
54
      gamma = cam_omega
 
55
c
 
56
      do n = 1,nq
 
57
         if (ipol.eq.1) then   ! spin-restricted
 
58
            rA(n) = rho(n,1)
 
59
            rB(n) = 0.d0
 
60
         else                  ! spin-unrestricted
 
61
            rA(n) = rho(n,2)
 
62
            rB(n) = rho(n,3)
 
63
         end if
 
64
      end do
 
65
c
 
66
c     -----------------------------------------------------------------------
 
67
c     Calculate the first and second derivatives
 
68
c     -----------------------------------------------------------------------
 
69
c
 
70
      total = 0.d0
 
71
      do n = 1,nq
 
72
         rhoA = rA(n)
 
73
         rhoB = rB(n)
 
74
         rhoTotal  = rhoA + rhoB   ! total density at point
 
75
         if (rhoTotal.gt.tol_rho) then
 
76
 
 
77
            if (ipol.eq.1) then    ! spin-restricted
 
78
              rhoA1 = rhoA
 
79
              rhoB1 = rhoB
 
80
            else                   ! spin-unrestricted
 
81
              rhoA1 = rhoA*2.d0
 
82
              rhoB1 = rhoB*2.d0
 
83
            end if
 
84
 
 
85
            fA   = EpsX(rhoA1,gamma)
 
86
            fB   = EpsX(rhoB1,gamma)
 
87
            fpA  = EpsXprime(rhoA1,gamma)
 
88
            fpB  = EpsXprime(rhoB1,gamma)
 
89
 
 
90
            f(n) = fA * rhoA + fB * rhoB
 
91
            Amat(n,1) = Amat(n,1) + (fpA*rhoA1+fA)*fac
 
92
            if (ipol.gt.1) then
 
93
              Amat(n,2) = Amat(n,2) + (fpB*rhoB1+fB)*fac
 
94
            end if
 
95
 
 
96
#ifdef SECOND_DERIV
 
97
            fppA = EpsTwoXprime(rhoA1,gamma)
 
98
            fppB = EpsTwoXprime(rhoB1,gamma)
 
99
c
 
100
            if (ipol.eq.1) then
 
101
             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
 
102
     &         (fppA*rhoA+2.0d0*fpA)*fac*2.0d0
 
103
            else
 
104
             Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) +
 
105
     &         (fppA*rhoA+fpA)*fac*4.0d0
 
106
c            Guard against case of no beta electrons, e.g. H atom
 
107
             if (rho(n,3).gt.tol_rho) then
 
108
                 Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) +
 
109
     &               (fppB*rhoB+fpB)*fac*4.0d0
 
110
             end if
 
111
            end if
 
112
#endif
 
113
            if (ldew) func(n) = func(n) + f(n)*fac 
 
114
            total = total + f(n)*qwght(n)
 
115
         end if
 
116
      end do
 
117
 
 
118
      Ex = Ex + total*fac
 
119
 
 
120
      return
 
121
      end
 
122
c
 
123
#ifndef SECOND_DERIV
 
124
#define SECOND_DERIV
 
125
c
 
126
#include "xc_bnl.F"
 
127
c
 
128
c     ---------------------------------------------------------------------------------------
 
129
c     Utility functions
 
130
c     ---------------------------------------------------------------------------------------
 
131
c
 
132
c     ---------------------------------------------------------------------------------------
 
133
c     Return the value of pi
 
134
c     ---------------------------------------------------------------------------------------
 
135
c
 
136
      double precision function ValueOfPi()
 
137
c 
 
138
      implicit none
 
139
#include "xc_params.fh"
 
140
c      
 
141
      ValueOfPi = pi          
 
142
 
 
143
      return
 
144
      end
 
145
c
 
146
c     ---------------------------------------------------------------------------------------
 
147
c     Evaluates the actual function
 
148
c     ---------------------------------------------------------------------------------------
 
149
c
 
150
      double precision function HqBNL(q)
 
151
 
 
152
      implicit none
 
153
#include "xc_params.fh"
 
154
 
 
155
      double precision q,TwoSqrtPi,OneOverQ,q2,DERF
 
156
 
 
157
      OneOverQ = 1.0d0/q
 
158
      TwoSqrtPi = 2.0d0*dsqrt(pi) 
 
159
      q2 = q**2.0d0
 
160
 
 
161
      if (q .lt. 1D-15) then
 
162
         HqBNL=1.d0
 
163
         return
 
164
      end if
 
165
 
 
166
      if (q .lt. 0.1d0) then
 
167
         HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi-q+q*(q2-2.0d0))
 
168
         return
 
169
      end if
 
170
 
 
171
      HqBNL=1.0d0-q*2.0d0/3.0d0*(TwoSqrtPi*DERF(OneOverQ)-q+
 
172
     $     q*(q2-2.0d0)*(1.0d0-exp(-OneOverQ*OneOverQ)))
 
173
 
 
174
      return
 
175
      end
 
176
c
 
177
c     ---------------------------------------------------------------------------------------
 
178
c     Calculate the local Fermi vector for the provided density
 
179
c     ---------------------------------------------------------------------------------------
 
180
c
 
181
      double precision function FermiK(den)
 
182
 
 
183
      implicit none
 
184
#include "xc_params.fh"
 
185
 
 
186
      double precision F13, den
 
187
 
 
188
      F13 = 1.0D0 / 3.0D0
 
189
      FermiK = (3.d0*pi*pi*den)**F13
 
190
 
 
191
      return
 
192
      end
 
193
c
 
194
c     ---------------------------------------------------------------------------------------
 
195
c     Calculate the function EpsX at the given density value and gamma
 
196
c     ---------------------------------------------------------------------------------------
 
197
c
 
198
      double precision function EpsX(Rho,gamma)
 
199
 
 
200
      implicit none
 
201
#include "xc_params.fh"
 
202
 
 
203
      double precision  kF,RHO,gamma,Cs
 
204
      double precision HqBNL
 
205
      double precision FermiK
 
206
 
 
207
      if (RHO.le.0D0) then
 
208
         EpsX = 0.0D0
 
209
         return
 
210
      end if
 
211
 
 
212
      kF = FermiK(Rho)
 
213
      Cs = -3.0D0/(4.0d0*pi)
 
214
      EpsX = Cs * kF * HqBNL(gamma/kF)
 
215
 
 
216
      return
 
217
      end      
 
218
c
 
219
c     ---------------------------------------------------------------------------------------
 
220
c     Calculate the first derivative of the function
 
221
c     ---------------------------------------------------------------------------------------
 
222
c
 
223
      double precision function HqBNLPrime(q)
 
224
 
 
225
      implicit none
 
226
#include "xc_params.fh"
 
227
 
 
228
      double precision q,OneOverQ,q2,q3,DERF
 
229
 
 
230
      OneOverQ = 1.0d0/q
 
231
      q2 = q**2.0d0
 
232
      q3 = q**3.0d0
 
233
 
 
234
      if (q .lt. 0.1d0) then
 
235
        HqBNLPrime = -4.0d0/3.0d0*(dsqrt(Pi)+2.0d0*q3-3.0d0*q)
 
236
        return
 
237
      end if
 
238
 
 
239
      HqBNLPrime = 4.0d0/3.0d0*(q*(exp(-OneOverQ*OneOverQ)*(2.0d0*q2
 
240
     $     -1.0d0)+(3.0d0-2.0d0*q2))-dsqrt(Pi)*DERF(OneOverQ))
 
241
 
 
242
      return
 
243
      end
 
244
c
 
245
c     ---------------------------------------------------------------------------------------
 
246
c     Calculate the first derivative of the local Fermi vector (it depends on the density)
 
247
c     ---------------------------------------------------------------------------------------
 
248
c
 
249
      double precision function FermiKPrime(den)
 
250
 
 
251
      implicit none
 
252
#include "xc_params.fh"
 
253
   
 
254
      double precision F23, den
 
255
 
 
256
      F23 = 2.0D0 / 3.0D0
 
257
      FermiKPrime = (Pi/(3.0d0*den))**F23
 
258
 
 
259
      return
 
260
      end
 
261
c
 
262
c     ---------------------------------------------------------------------------------------
 
263
c     Calculate the first derivative of q (q=gamma/kf) (it implicitly depends on the density)
 
264
c     ---------------------------------------------------------------------------------------
 
265
c
 
266
      double precision function QPrime(gamma,kF)
 
267
 
 
268
      implicit none
 
269
 
 
270
      double precision  kF, FermiK2, gamma
 
271
 
 
272
      FermiK2 = kF**2.0d0
 
273
      QPrime = -gamma/FermiK2
 
274
 
 
275
      return
 
276
      end
 
277
c
 
278
c     ---------------------------------------------------------------------------------------
 
279
c     Calculate the first derivative of EpsX
 
280
c     ---------------------------------------------------------------------------------------
 
281
c
 
282
      double precision function EpsXprime(Rho,gamma)
 
283
 
 
284
      implicit none
 
285
#include "xc_params.fh"
 
286
 
 
287
      double precision Rho,gamma
 
288
      double precision Cs,kF,CsPrime
 
289
 
 
290
      double precision HqBNL
 
291
      double precision HqBNLPrime
 
292
      double precision QPrime
 
293
      double precision FermiK
 
294
      double precision FermiKPrime
 
295
 
 
296
      kF = FermiK(Rho)
 
297
      CsPrime = -3.0D0/(4.0d0*Pi)
 
298
      Cs = CsPrime*kF
 
299
 
 
300
      if (Rho.le.0d0) then
 
301
         EpsXprime = 0.0d0
 
302
         return
 
303
      end if
 
304
 
 
305
      EpsXprime = FermiKPrime(Rho)*(CsPrime*HqBNL(gamma/kF)+
 
306
     $     QPrime(gamma,kF)*HqBNLPrime(gamma/kF)*Cs)
 
307
 
 
308
      return
 
309
      end
 
310
c
 
311
c     ---------------------------------------------------------------------------------------
 
312
c     Calculate the second derivative of the main function that consititutes the functional
 
313
c     ---------------------------------------------------------------------------------------
 
314
c
 
315
      double precision function HqBNLTwoPrime(q)
 
316
 
 
317
      implicit none
 
318
#include "xc_params.fh"
 
319
 
 
320
      double precision q,OneOverQ,q2
 
321
 
 
322
      OneOverQ = 1.0d0/q
 
323
      q2 = q**2.0d0
 
324
 
 
325
      if (q .lt. 0.1d0) then
 
326
         HqBNLTwoPrime = 4.0d0-8.0d0*q2
 
327
         return
 
328
      end if
 
329
 
 
330
      HqBNLTwoPrime = exp(-OneOverQ*OneOverQ)*(4.0d0+8.0d0*q2)
 
331
     $     -8.0d0*q2+4.0d0
 
332
 
 
333
      return
 
334
      end
 
335
c
 
336
c     ---------------------------------------------------------------------------------------
 
337
c     Calculate the second derivative of the local Fermi vector
 
338
c     ---------------------------------------------------------------------------------------
 
339
c
 
340
      double precision function FermiKTwoPrime(den)
 
341
 
 
342
      implicit none
 
343
#include "xc_params.fh"
 
344
 
 
345
      double precision F13, den
 
346
 
 
347
      F13 = 1.0D0/3.0D0
 
348
      FermiKTwoPrime =  -(8.0d0*Pi**2.0d0/(243.0d0*den**5.0d0))**F13
 
349
 
 
350
      return
 
351
      end
 
352
c
 
353
c     ---------------------------------------------------------------------------------------
 
354
c     Calculate the second derivative of q    
 
355
c     ---------------------------------------------------------------------------------------
 
356
c
 
357
      double precision function QTwoPrime(gamma,kF)
 
358
 
 
359
      implicit none
 
360
 
 
361
      double precision gamma, kF, FermiK3
 
362
 
 
363
      FermiK3 = kF**3.0d0
 
364
      QTwoPrime = (2.0d0*gamma)/FermiK3
 
365
 
 
366
      return
 
367
      end
 
368
c
 
369
c     ---------------------------------------------------------------------------------------
 
370
c     Calculate the second derivative of EpsX
 
371
c     ---------------------------------------------------------------------------------------
 
372
c
 
373
      double precision function EpsTwoXprime(Rho,gamma)
 
374
 
 
375
      implicit none
 
376
#include "xc_params.fh"
 
377
 
 
378
      double precision Rho,gamma
 
379
      double precision kF,kFPrim,kFprim2,kF2prim
 
380
      double precision q,qprim,qprim2,q2prim
 
381
      double precision g,gprim,g2prim
 
382
      double precision Cs,CsPrim
 
383
 
 
384
      double precision FermiK
 
385
      double precision FermiKPrime
 
386
      double precision FermiKTwoPrime
 
387
      double precision QPrime
 
388
      double precision QTwoPrime
 
389
      double precision HqBNL
 
390
      double precision HqBNLPrime
 
391
      double precision HqBNLTwoPrime
 
392
 
 
393
      if (Rho.le.0d0) then
 
394
         EpsTwoXprime = 0.0d0
 
395
         return
 
396
      end if
 
397
 
 
398
      kF = FermiK(Rho)
 
399
      kFPrim = FermiKPrime(Rho)
 
400
      kFPrim2=kFPrim**2.0d0
 
401
      kF2prim = FermiKTwoPrime(Rho)
 
402
      CsPrim = -3.0d0/(4.0d0*Pi)
 
403
      Cs = CsPrim * kF
 
404
      q = gamma / kF
 
405
      qprim = QPrime(gamma,kF)
 
406
      Qprim2=qprim**2.0d0
 
407
      q2prim = QTwoPrime(gamma,kF)
 
408
      g = HqBNL(q)
 
409
      gprim = HqBNLPrime(q)
 
410
      g2prim = HqBNLTwoPrime(q)
 
411
 
 
412
      EpsTwoXprime = 
 
413
     $     kFPrim2*(2.0d0*CsPrim*gprim*qprim
 
414
     $     +Cs*(QPrim2*g2prim+gprim*Q2Prim))
 
415
     $     +kF2Prim*(g*CsPrim+Cs*gprim*qprim)
 
416
 
 
417
      return
 
418
      end
 
419
c
 
420
#endif
 
421
c $Id: xc_bnl.F 24192 2013-05-04 20:14:33Z niri $