← Back to branch summary

~ubuntu-branches/ubuntu/trusty/nwchem/trusty-proposed

« back to all changes in this revision

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