← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_cams12x.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_cams12x.F
 
1
#if defined(FUJITSU_VPP)
 
2
!ocl scalar
 
3
#endif
 
4
#ifndef SECOND_DERIV
 
5
      Subroutine xc_cams12x(tol_rho, fac, lfac, nlfac, rho, delrho, 
 
6
     &              Amat, Cmat, nq, ipol, Ex, qwght,ldew,func,is12x)
 
7
#else
 
8
      Subroutine xc_cams12x_d2(tol_rho, fac, lfac, nlfac, rho, delrho, 
 
9
     &               Amat, Amat2, Cmat, Cmat2, nq, ipol, Ex,
 
10
     &                         qwght,ldew,func,is12x)
 
11
#endif
 
12
c
 
13
C$Id: xc_cams12x.F 20247 2011-04-28 18:58:49Z d3y133 $
 
14
c
 
15
      implicit none
 
16
c      
 
17
#include "dft2drv.fh"
 
18
c
 
19
      double precision tol_rho, fac, Ex
 
20
      integer nq, ipol, is12x
 
21
      logical lfac, nlfac,ldew
 
22
      double precision func(*)  ! value of the functional [output]
 
23
c
 
24
c     Charge Density
 
25
c
 
26
      double precision rho(nq,ipol*(ipol+1)/2)
 
27
c
 
28
c     Charge Density Gradient
 
29
c
 
30
      double precision delrho(nq,3,ipol)
 
31
c
 
32
c     Quadrature Weights
 
33
c
 
34
      double precision qwght(nq)
 
35
c
 
36
c     Sampling Matrices for the XC Potential
 
37
c
 
38
      double precision Amat(nq,ipol), Cmat(nq,*)
 
39
      double precision Atmp, Ctmp, Etmp
 
40
      double precision A2tmp, C2tmp, C3tmp
 
41
c
 
42
#ifdef SECOND_DERIV
 
43
c
 
44
c     Second Derivatives of the Exchange Energy Functional
 
45
c
 
46
      double precision Amat2(nq,NCOL_AMAT2), Cmat2(nq,NCOL_CMAT2)
 
47
#endif
 
48
c
 
49
      double precision rB, rC, rD, rA, rK, rE, rH, rG2, rH2
 
50
      double precision ums, vms, dudx, dvdx, d2udx2, d2vdx2
 
51
c
 
52
c References:
 
53
c
 
54
c    Swart, J. Chem. Phys. (2012/2013)
 
55
c
 
56
c***************************************************************************
 
57
c
 
58
      integer n
 
59
      double precision C, rho13, rho43, gamma, x, g, gdenom, dg
 
60
      double precision dgdenom, t, x2, x3, x4, g1, g2
 
61
      double precision g1h1, g2h1, g1h2, g2h2, g1h3, g3h1
 
62
      double precision hdenom, dhdenom, d2hdenom, PI, rM
 
63
      double precision gc4, dgc4, d2gc4
 
64
      parameter (rM=60.770665d0)
 
65
      parameter (PI = 3.1415926535897932385d0)
 
66
#ifdef SECOND_DERIV
 
67
      double precision rhom23, d2g, d2gdenom, d2g1, d2g2, d2g3
 
68
#endif
 
69
c
 
70
      if (is12x.eq.1) then
 
71
c
 
72
cswar1      1.03323556     0.00417251     0.00115216     0.75700000     0.00000000
 
73
cswar2      0.00706184     1.20250451     0.86124355     0.00000000     0.34485046
 
74
cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     1.52420731
 
75
c
 
76
        rA = 1.03323556d0
 
77
        rK = 0.757d0
 
78
        rC = 0.00417251d0
 
79
        rD = 0.00115216d0
 
80
        rE = 0.00706184d0
 
81
      elseif (is12x.eq.2) then
 
82
c
 
83
cswar1      1.02149642     0.00825905     0.00235804     0.75700000     0.25000000
 
84
cswar2      0.00654977     1.08034183     0.37999939     0.00000000     0.35897845
 
85
cswar3      1.00000000     0.00000000     0.00000000     0.00000000     1.00000000     0.48516891
 
86
c
 
87
        rA = 1.02149642d0
 
88
        rK = 0.757d0
 
89
        rC = 0.00825905d0
 
90
        rD = 0.00235804d0
 
91
        rE = 0.00654977d0
 
92
      else
 
93
        stop 'error in xc_cams12x.F'
 
94
      endif
 
95
      rB = 1d0 + rK - rA
 
96
c
 
97
c     Uniform electron gas constant
 
98
c
 
99
      C = -(1.5d0)*(0.75d0/acos(-1d0))**(1d0/3d0)
 
100
c
 
101
      if (ipol.eq.1) then
 
102
c
 
103
c        ======> SPIN-RESTRICTED <======
 
104
c
 
105
         do 10 n = 1, nq
 
106
            if (rho(n,1).lt.tol_rho) goto 10
 
107
c
 
108
c           Spin alpha:
 
109
c
 
110
            rho13 = (0.5d0*rho(n,1))**(1.d0/3.d0)
 
111
            rho43 = rho13**4 
 
112
c
 
113
            Etmp = 0.d0
 
114
            Atmp = 0.d0
 
115
            Ctmp = 0.d0
 
116
            if (lfac) then
 
117
               Etmp = rA * 2d0*rho43*C*fac
 
118
               Atmp = rA * (4d0/3d0)*rho13*C*fac
 
119
            endif
 
120
c
 
121
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
122
     &              delrho(n,2,1)*delrho(n,2,1) +
 
123
     &              delrho(n,3,1)*delrho(n,3,1)
 
124
            if (dsqrt(gamma).gt.tol_rho)then
 
125
               gamma = 0.25d0 * gamma
 
126
               x = dsqrt(gamma) / rho43
 
127
               x2 = x*x
 
128
            else
 
129
               x = 0d0
 
130
               x2 = 0d0
 
131
            endif
 
132
c
 
133
            gdenom = 1d0 + rC*x2 + rD*x2*x2
 
134
            hdenom = 1d0 + rE*x2
 
135
            ums = 1d0 - 1d0 / gdenom
 
136
            vms = 1d0 - 1d0 / hdenom
 
137
            g = C*rB*ums*vms
 
138
c
 
139
            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
 
140
            dvdx = 2d0*rE*x/(hdenom**2)
 
141
            dg = C*rB*(dudx*vms + ums*dvdx)
 
142
c
 
143
            if (nlfac) then
 
144
               Etmp = Etmp + 2d0*rho43*g*fac
 
145
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
146
            endif
 
147
c
 
148
            if (x.gt.tol_rho) then
 
149
               Ctmp = 0.5d0 * dg / sqrt(gamma) * fac
 
150
            endif
 
151
c
 
152
#ifdef SECOND_DERIV
 
153
c
 
154
c           Add local contribution back to g
 
155
c
 
156
            if(lfac) g = g + rA * C
 
157
c
 
158
            rhom23 = rho13 / (0.5d0*rho(n,1))
 
159
 
 
160
            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
 
161
     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
 
162
            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
 
163
            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
 
164
c
 
165
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
166
            C2tmp = - (4d0/3d0)*(rhom23**2/rho(n,1))*d2g*fac
 
167
            if (x.gt.tol_rho) then
 
168
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
169
            else
 
170
               C3tmp = 0d0
 
171
            endif
 
172
c
 
173
            call xc_att_xc_d2(rho(n,1),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
174
     &           C2tmp,C3tmp)
 
175
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
176
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
177
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
178
#else
 
179
            call xc_att_xc(rho(n,1),ipol,Etmp,Atmp,Ctmp)
 
180
#endif
 
181
            Ex = Ex + qwght(n)*Etmp
 
182
            if (ldew) func(n) = func(n) + Etmp
 
183
            Amat(n,1) = Amat(n,1) + Atmp
 
184
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
185
c
 
186
 10      continue
 
187
c
 
188
      else
 
189
c
 
190
c        ======> SPIN-UNRESTRICTED <======
 
191
c
 
192
         do 20 n = 1, nq
 
193
            if (rho(n,1).lt.tol_rho) goto 20
 
194
            if (rho(n,2).lt.tol_rho) goto 25
 
195
c
 
196
c           Spin alpha:
 
197
c
 
198
            rho13 = rho(n,2)**(1.d0/3.d0)
 
199
            rho43 = rho13*rho(n,2)
 
200
c
 
201
            Etmp = 0.d0
 
202
            Atmp = 0.d0
 
203
            Ctmp = 0.d0
 
204
            if (lfac) then
 
205
               Etmp = rA * rho43*C*fac
 
206
               Atmp = rA * (4d0/3d0)*rho13*C*fac
 
207
            endif
 
208
c
 
209
            gamma = delrho(n,1,1)*delrho(n,1,1) +
 
210
     &              delrho(n,2,1)*delrho(n,2,1) +
 
211
     &              delrho(n,3,1)*delrho(n,3,1)
 
212
            if (dsqrt(gamma).gt.tol_rho)then
 
213
               x = dsqrt(gamma) / rho43
 
214
               x2 = x*x
 
215
            else
 
216
               x = 0d0
 
217
               x2 = 0d0
 
218
            endif
 
219
c
 
220
            gdenom = 1d0 + rC*x2 + rD*x2*x2
 
221
            hdenom = 1d0 + rE*x2
 
222
            ums = 1d0 - 1d0 / gdenom
 
223
            vms = 1d0 - 1d0 / hdenom
 
224
            g = C*rB*ums*vms
 
225
c
 
226
            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
 
227
            dvdx = 2d0*rE*x/(hdenom**2)
 
228
            dg = C*rB*(dudx*vms + ums*dvdx)
 
229
c
 
230
            if (nlfac) then
 
231
               Etmp = Etmp + rho43*g*fac
 
232
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
233
            endif
 
234
c
 
235
            if (x.gt.tol_rho) then
 
236
               t = dg / sqrt(gamma) * fac
 
237
               Ctmp = t * 0.5d0
 
238
            endif
 
239
c
 
240
#ifdef SECOND_DERIV
 
241
c
 
242
c           Add local contribution back to g
 
243
c
 
244
            if (lfac) g = g + rA * C
 
245
c
 
246
            rhom23 = rho13 / rho(n,2)
 
247
c
 
248
            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
 
249
     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
 
250
            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
 
251
            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
 
252
c
 
253
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
254
            C2tmp = - (2d0/3d0)*(rhom23**2/rho(n,2))*d2g*fac
 
255
            if (x.gt.tol_rho) then
 
256
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
257
            else
 
258
               C3tmp = 0d0
 
259
            endif
 
260
 
 
261
            call xc_att_xc_d2(rho(n,2),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
262
     &           C2tmp,C3tmp)
 
263
            Amat2(n,D2_RA_RA) = Amat2(n,D2_RA_RA) + A2tmp
 
264
            Cmat2(n,D2_RA_GAA) = Cmat2(n,D2_RA_GAA) + C2tmp
 
265
            Cmat2(n,D2_GAA_GAA) = Cmat2(n,D2_GAA_GAA) + C3tmp
 
266
#else
 
267
            call xc_att_xc(rho(n,2),ipol,Etmp,Atmp,Ctmp)
 
268
#endif
 
269
            Ex = Ex + qwght(n)*Etmp
 
270
            if (ldew) func(n) = func(n) + Etmp
 
271
            Amat(n,1) = Amat(n,1) + Atmp
 
272
            Cmat(n,D1_GAA) = Cmat(n,D1_GAA) + Ctmp
 
273
c
 
274
 25         continue
 
275
c
 
276
c           Spin beta:
 
277
c
 
278
            if (rho(n,3).lt.tol_rho) goto 20
 
279
c
 
280
            rho13 = rho(n,3)**(1.d0/3.d0)
 
281
            rho43 = rho13*rho(n,3)
 
282
c
 
283
            Etmp = 0.d0
 
284
            Atmp = 0.d0
 
285
            Ctmp = 0.d0
 
286
            if (lfac) then
 
287
               Etmp = rA * rho43*C*fac
 
288
               Atmp = rA * (4d0/3d0)*rho13*C*fac
 
289
            endif
 
290
c
 
291
            gamma = delrho(n,1,2)*delrho(n,1,2) +
 
292
     &              delrho(n,2,2)*delrho(n,2,2) +
 
293
     &              delrho(n,3,2)*delrho(n,3,2)
 
294
            if (dsqrt(gamma).gt.tol_rho)then
 
295
               x = dsqrt(gamma) / rho43
 
296
               x2 = x*x
 
297
            else
 
298
               x = 0d0
 
299
               x2 = 0d0
 
300
            endif
 
301
c
 
302
            gdenom = 1d0 + rC*x2 + rD*x2*x2
 
303
            hdenom = 1d0 + rE*x2
 
304
            ums = 1d0 - 1d0 / gdenom
 
305
            vms = 1d0 - 1d0 / hdenom
 
306
            g = C*rB*ums*vms
 
307
c
 
308
            dudx = (2d0*rC*x + 4d0*rD*x2*x)/(gdenom**2)
 
309
            dvdx = 2d0*rE*x/(hdenom**2)
 
310
            dg = C*rB*(dudx*vms + ums*dvdx)
 
311
c
 
312
            if (nlfac) then
 
313
               Etmp = Etmp + rho43*g*fac
 
314
               Atmp = Atmp + (4d0/3d0)*rho13*(g-x*dg)*fac
 
315
            endif
 
316
c
 
317
            if (x.gt.tol_rho) then
 
318
               t = dg / sqrt(gamma) * fac
 
319
               Ctmp = t * 0.5d0
 
320
            endif
 
321
c
 
322
#ifdef SECOND_DERIV
 
323
c
 
324
c           Add local contribution back to g
 
325
c
 
326
            if (lfac) g = g + rA * C
 
327
c
 
328
            rhom23 = rho13 / rho(n,3)
 
329
c
 
330
            d2udx2 = (2d0*rC-6d0*rC*rC*x2+12d0*rD*x2-18d0*rC*rD*x2*x2
 
331
     &                -20d0*rD*rD*x2*x2*x2)/(gdenom**3)
 
332
            d2vdx2 = (2d0*rE - 6d0*rE*rE*x2)/(hdenom**3)
 
333
            d2g = C*rB*(d2udx2*vms + 2d0*dudx*dvdx + ums*d2vdx2)
 
334
 
 
335
c
 
336
            A2tmp = (4d0/9d0)*rhom23*(g-x*dg+4d0*x*x*d2g)*fac
 
337
            C2tmp = -(2d0/3d0)*(rhom23**2/rho(n,3))*d2g*fac
 
338
            if (x.gt.tol_rho) then
 
339
               C3tmp = - 0.25d0*gamma**(-1.5d0)*(dg-x*d2g)*fac
 
340
            else
 
341
               C3tmp = 0d0
 
342
            endif
 
343
c
 
344
            call xc_att_xc_d2(rho(n,3),ipol,Etmp,Atmp,Ctmp,A2tmp,
 
345
     &           C2tmp,C3tmp)
 
346
            Amat2(n,D2_RB_RB) = Amat2(n,D2_RB_RB) + A2tmp
 
347
            Cmat2(n,D2_RB_GBB) = Cmat2(n,D2_RB_GBB) + C2tmp
 
348
            Cmat2(n,D2_GBB_GBB) = Cmat2(n,D2_GBB_GBB) + C3tmp
 
349
#else
 
350
            call xc_att_xc(rho(n,3),ipol,Etmp,Atmp,Ctmp)
 
351
#endif
 
352
            Ex = Ex + qwght(n)*Etmp
 
353
            if (ldew) func(n) = func(n) + Etmp
 
354
            Amat(n,2) = Amat(n,2) + Atmp
 
355
            Cmat(n,D1_GBB) = Cmat(n,D1_GBB) + Ctmp
 
356
c
 
357
 20      continue
 
358
c
 
359
      endif
 
360
c
 
361
      return
 
362
      end
 
363
#ifndef SECOND_DERIV
 
364
#define SECOND_DERIV
 
365
c
 
366
c     Compile source again for the 2nd derivative case
 
367
c
 
368
#include "xc_cams12x.F"
 
369
#endif