← Back to branch summary

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

« back to all changes in this revision

Viewing changes to src/nwdft/xc/xc_xsogga.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_xsogga.F
 
1
      Subroutine xc_xsogga(tol_rho, fac,lfac,nlfac, rho, delrho, 
 
2
     &                     Amat, Cmat, nq, ipol, Ex, 
 
3
     &                     qwght, ldew, func, ijzy)
 
4
c   
 
5
c$Id: xc_xsogga.F 22503 2012-05-20 06:58:57Z d3y133 $
 
6
c
 
7
c
 
8
c
 
9
c**********************************************************************c
 
10
c                                                                      c
 
11
c  SOGGA11X evaluates the exchange part of the SOGGA, SOGGA11          c
 
12
c  and SOGGA11-X functionals on the grid.                              c
 
13
c                                                                      c
 
14
c     a) Zhao and Truhlar, J.Chem.Phys., 128, 184109 (2008)            c
 
15
c     b) Peverati, Zhao and Truhlar, J.Phys.Chem.Lett, 2, 1991 (2011)  c
 
16
c     c) Peverati and Truhlar, J.Chem.Phys, 135, 191102 (2011)         c
 
17
c                                                                      c
 
18
c      ijzy = 1 - SOGGA functional (a) - it requres PBE correlation    c
 
19
c      ijzy = 2 - SOGGA11 functional (b)                               c
 
20
c      ijzy = 3 - SOGGA11-X functional (c)                             c
 
21
c                                                                      c
 
22
c Coded by Roberto Peverati (12/11)                                    c
 
23
c                                                                      c
 
24
c**********************************************************************c    
 
25
c
 
26
      implicit none
 
27
c      
 
28
      double precision fac, Ex
 
29
      integer nq, ipol
 
30
      logical lfac, nlfac,ldew
 
31
      double precision func(*)  ! value of the functional [output]
 
32
c
 
33
c     Charge Density & Its Cube Root
 
34
c
 
35
      double precision rho(nq,ipol*(ipol+1)/2)
 
36
c
 
37
c     Charge Density Gradient
 
38
c
 
39
      double precision delrho(nq,3,ipol)
 
40
c
 
41
c     Quadrature Weights
 
42
c
 
43
      double precision qwght(nq)
 
44
c
 
45
c     Sampling Matrices for the XC Potential & Energy
 
46
c
 
47
      double precision Amat(nq,ipol), Cmat(nq,*)
 
48
      double precision tol_rho, pi
 
49
c
 
50
c     Intermediate derivative results, etc.
 
51
c
 
52
      integer n, ijzy
 
53
c
 
54
      double precision CxA0,CxA1,CxA2,CxA3,CxA4,CxA5
 
55
      double precision CxB0,CxB1,CxB2,CxB3,CxB4,CxB5
 
56
      double precision rho43, rho13, rhoo
 
57
c
 
58
      double precision AS, ASO, AX, DELOCDR
 
59
      double precision DFA1DG, DFA1DR, DFA1DY
 
60
      double precision DFA2DG, DFA2DR, DFA2DY
 
61
      double precision DFA3DG, DFA3DR, DFA3DY
 
62
      double precision DFA4DG, DFA4DR, DFA4DY
 
63
      double precision DFA5DG, DFA5DR, DFA5DY
 
64
      double precision DFB1DG, DFB1DR, DFB1DY
 
65
      double precision DFB2DG, DFB2DR, DFB2DY
 
66
      double precision DFB3DG, DFB3DR, DFB3DY
 
67
      double precision DFB4DG, DFB4DR, DFB4DY
 
68
      double precision DFB5DG, DFB5DR, DFB5DY
 
69
      double precision DFEXPDPON, DFFRACDPON, DFGGAXDG, DFGGAXDR
 
70
      double precision DYDG, DYDR, DTOL, ELOC
 
71
      double precision FA0, FA1, FA2, FA3, FA4, FA5
 
72
      double precision FB0, FB1, FB2, FB3, FB4, FB5
 
73
      double precision FEXP, FFRAC, FGGAX, MU, PON, S, X, Y       
 
74
      double precision Gam12, Gam      
 
75
c      
 
76
      double precision f1,f2,f3,f4,f5,f8
 
77
      double precision F1o3,F4o3,F48
 
78
      parameter( F1=1.0D+00,  F2=2.0D+00,  F3=3.0D+00,  
 
79
     $           F4=4.0D+00,  F5=5.0D+00,  F8=8.0D+00,
 
80
     $           F48=48.0D+00)
 
81
c
 
82
        pi=acos(-1d0)      
 
83
c
 
84
      if (ijzy.eq.1) then
 
85
c SOGGA11
 
86
       CxA0 = 0.5d0
 
87
       CxA1 = 0.276d0
 
88
       CxA2 = 0d0
 
89
       CxA3 = 0d0
 
90
       CxA4 = 0d0
 
91
       CxA5 = 0d0
 
92
       CxB0 = 0.5d0
 
93
       CxB1 = 0.276d0
 
94
       CxB2 = 0d0
 
95
       CxB3 = 0d0
 
96
       CxB4 = 0d0
 
97
       CxB5 = 0d0
 
98
      elseif (ijzy.eq.2) then
 
99
c SOGGA11
 
100
       CxA0 =  0.50000d0
 
101
       CxA1 = -2.95535d0
 
102
       CxA2 =  15.7974d0
 
103
       CxA3 = -91.1804d0
 
104
       CxA4 =  96.2030d0
 
105
       CxA5 =  0.18683d0
 
106
       CxB0 =  0.50000d0
 
107
       CxB1 =  3.50743d0
 
108
       CxB2 = -12.9523d0
 
109
       CxB3 =  49.7870d0
 
110
       CxB4 = -33.2545d0
 
111
       CxB5 = -11.1396d0
 
112
      elseif (ijzy.eq.3) then
 
113
c SOGGA11-X
 
114
       CxA0 =  2.99250d-01
 
115
       CxA1 =  3.21638d+00
 
116
       CxA2 = -3.55605d+00
 
117
       CxA3 =  7.65852d+00
 
118
       CxA4 = -1.12830d+01
 
119
       CxA5 =  5.25813d+00
 
120
       CxB0 =  2.99250d-01
 
121
       CxB1 = -2.88595d+00
 
122
       CxB2 =  3.23617d+00
 
123
       CxB3 = -2.45393d+00
 
124
       CxB4 = -3.75495d+00
 
125
       CxB5 =  3.96613d+00
 
126
      endif
 
127
c
 
128
      DTol = tol_rho
 
129
      F1o3 = F1/F3 
 
130
      F4o3 = F4/F3
 
131
      Pi   = ACos(-F1)
 
132
      AsO  = (F48*PI*PI)**F1o3
 
133
      As   = F1/AsO
 
134
      Ax   = -(F3/F2) * (F4o3*Pi)**(-F1o3)
 
135
      mu = 0.2236536053d0
 
136
c
 
137
      if (ipol.eq.1 )then
 
138
c
 
139
c        ======> SPIN-RESTRICTED <======
 
140
c                     or
 
141
c                SPIN-UNPOLARIZED
 
142
c
 
143
c
 
144
         do 10 n = 1, nq
 
145
            if (rho(n,1).lt.DTol) goto 10
 
146
            rhoo = rho(n,1)/F2
 
147
            rho43 = rhoo**F4o3  
 
148
            rho13 = rho43/rhoo
 
149
            Gam =(delrho(n,1,1)*delrho(n,1,1) +
 
150
     &              delrho(n,2,1)*delrho(n,2,1) +
 
151
     &              delrho(n,3,1)*delrho(n,3,1))/F4
 
152
            Gam12 = dsqrt(Gam)
 
153
            if(gam12.lt.dtol) goto 10
 
154
c
 
155
          Eloc = Ax*Rho43
 
156
          x = Gam12/Rho43
 
157
          s = As*x
 
158
          y = s*s
 
159
          PON = mu*y
 
160
          Ffrac = F1-F1/(F1+PON)
 
161
          Fexp  = F1-exp(-PON)
 
162
          fa0 = CxA0
 
163
          fa1 = CxA1 *Ffrac
 
164
          fa2 = CxA2 *Ffrac**F2
 
165
          fa3 = CxA3 *Ffrac**F3
 
166
          fa4 = CxA4 *Ffrac**F4
 
167
          fa5 = CxA5 *Ffrac**F5
 
168
          fb0 = CxB0
 
169
          fb1 = CxB1 *Fexp
 
170
          fb2 = CxB2 *Fexp**F2
 
171
          fb3 = CxB3 *Fexp**F3
 
172
          fb4 = CxB4 *Fexp**F4
 
173
          fb5 = CxB5 *Fexp**F5
 
174
c
 
175
          Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
 
176
     $            fb0+fb1+fb2+fb3+fb4+fb5
 
177
C
 
178
C     1st derivatives.
 
179
C
 
180
 
 
181
          dElocdR=Ax*F4o3*Rho13
 
182
          dydR = -(F8/F3)*y/Rhoo
 
183
          dydG   = y/Gam
 
184
          dFfracdPON = F1/((F1+PON)**F2)
 
185
          dFexpdPON  = exp(-PON)
 
186
          dfa1dy = CxA1 *mu*dFfracdPON
 
187
          dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
 
188
          dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
 
189
          dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
 
190
          dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
 
191
          dfa1dR = dfa1dy *dydR
 
192
          dfa2dR = dfa2dy *dydR
 
193
          dfa3dR = dfa3dy *dydR
 
194
          dfa4dR = dfa4dy *dydR
 
195
          dfa5dR = dfa5dy *dydR
 
196
          dfa1dG = dfa1dy *dydG
 
197
          dfa2dG = dfa2dy *dydG
 
198
          dfa3dG = dfa3dy *dydG
 
199
          dfa4dG = dfa4dy *dydG
 
200
          dfa5dG = dfa5dy *dydG
 
201
          dfb1dy = CxB1 *mu*dFexpdPON
 
202
          dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
 
203
          dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
 
204
          dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
 
205
          dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
 
206
          dfb1dR = dfb1dy *dydR
 
207
          dfb2dR = dfb2dy *dydR
 
208
          dfb3dR = dfb3dy *dydR
 
209
          dfb4dR = dfb4dy *dydR
 
210
          dfb5dR = dfb5dy *dydR
 
211
          dfb1dG = dfb1dy *dydG
 
212
          dfb2dG = dfb2dy *dydG
 
213
          dfb3dG = dfb3dy *dydG
 
214
          dfb4dG = dfb4dy *dydG
 
215
          dfb5dG = dfb5dy *dydG
 
216
c
 
217
          dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
 
218
     $               dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
 
219
c                     
 
220
          dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
 
221
     $               dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
 
222
c
 
223
           Ex = Ex +F2*(Eloc*Fggax)*qwght(n)
 
224
           if(ldew) func(n)=func(n)+F2*(Eloc*Fggax)
 
225
           Amat(n,1) = Amat(n,1)   +dElocdR*Fggax + Eloc*dFggaxdR        
 
226
           Cmat(n,1)=  Cmat(n,1)  + Eloc*dFggaxdG
 
227
           
 
228
10      continue
 
229
c
 
230
c UUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUUnrestricted
 
231
      else  ! ipol=2
 
232
c
 
233
c        ======> SPIN-UNRESTRICTED <======
 
234
 
 
235
c
 
236
c  use spin density functional theory ie n-->2n
 
237
c                               Ex=(1/2)Ex[2*alpha] + (1/2)Ex[2*beta]
 
238
c
 
239
c     Alpha            ALPHA               ALPHA
 
240
c
 
241
         do 20 n = 1, nq
 
242
           if (rho(n,1).lt.DTol) goto 20
 
243
           if (rho(n,2).lt.DTol) goto 25           
 
244
            rhoo  = rho(n,2)
 
245
            rho43 = rhoo**F4o3  
 
246
            rho13 = rho43/rhoo
 
247
            Gam =(delrho(n,1,1)*delrho(n,1,1) +
 
248
     &              delrho(n,2,1)*delrho(n,2,1) +
 
249
     &              delrho(n,3,1)*delrho(n,3,1))
 
250
            Gam12 = dsqrt(Gam)
 
251
         if(gam12.lt.dtol) goto 25
 
252
c
 
253
          Eloc = Ax*Rho43
 
254
          x = Gam12/Rho43
 
255
          s = As*x
 
256
          y = s*s
 
257
          PON = mu*y
 
258
          Ffrac = F1-F1/(F1+PON)
 
259
          Fexp  = F1-exp(-PON)
 
260
          fa0 = CxA0
 
261
          fa1 = CxA1 *Ffrac
 
262
          fa2 = CxA2 *Ffrac**F2
 
263
          fa3 = CxA3 *Ffrac**F3
 
264
          fa4 = CxA4 *Ffrac**F4
 
265
          fa5 = CxA5 *Ffrac**F5
 
266
          fb0 = CxB0
 
267
          fb1 = CxB1 *Fexp
 
268
          fb2 = CxB2 *Fexp**F2
 
269
          fb3 = CxB3 *Fexp**F3
 
270
          fb4 = CxB4 *Fexp**F4
 
271
          fb5 = CxB5 *Fexp**F5
 
272
c         
 
273
          Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
 
274
     $            fb0+fb1+fb2+fb3+fb4+fb5
 
275
C
 
276
C     1st derivatives.
 
277
C
 
278
           dElocdR=Ax*F4o3*Rho13
 
279
           dydR = -(F8/F3)*y/Rhoo
 
280
           dydG   = y/Gam
 
281
           dFfracdPON = F1/((F1+PON)**F2)
 
282
           dFexpdPON  = exp(-PON)
 
283
           dfa1dy = CxA1 *mu*dFfracdPON
 
284
           dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
 
285
           dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
 
286
           dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
 
287
           dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
 
288
           dfa1dR = dfa1dy *dydR
 
289
           dfa2dR = dfa2dy *dydR
 
290
           dfa3dR = dfa3dy *dydR
 
291
           dfa4dR = dfa4dy *dydR
 
292
           dfa5dR = dfa5dy *dydR
 
293
           dfa1dG = dfa1dy *dydG
 
294
           dfa2dG = dfa2dy *dydG
 
295
           dfa3dG = dfa3dy *dydG
 
296
           dfa4dG = dfa4dy *dydG
 
297
           dfa5dG = dfa5dy *dydG
 
298
           dfb1dy = CxB1 *mu*dFexpdPON
 
299
           dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
 
300
           dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
 
301
           dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
 
302
           dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
 
303
           dfb1dR = dfb1dy *dydR
 
304
           dfb2dR = dfb2dy *dydR
 
305
           dfb3dR = dfb3dy *dydR
 
306
           dfb4dR = dfb4dy *dydR
 
307
           dfb5dR = dfb5dy *dydR
 
308
           dfb1dG = dfb1dy *dydG
 
309
           dfb2dG = dfb2dy *dydG
 
310
           dfb3dG = dfb3dy *dydG
 
311
           dfb4dG = dfb4dy *dydG
 
312
           dfb5dG = dfb5dy *dydG
 
313
c
 
314
           dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
 
315
     $                dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
 
316
c                      
 
317
           dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
 
318
     $                dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
 
319
c
 
320
           Ex = Ex + (Eloc*Fggax)*qwght(n)
 
321
           if(ldew) func(n)=func(n)+ Eloc*Fggax
 
322
           Amat(n,1) = Amat(n,1)   + dElocdR*Fggax + Eloc*dFggaxdR
 
323
           Cmat(n,1)=  Cmat(n,1)   + Eloc*dFggaxdG
 
324
c
 
325
25         continue
 
326
c
 
327
c     Beta               BETA           BETA
 
328
c
 
329
            if (rho(n,3).lt.DTol) goto 20
 
330
            rhoo  = rho(n,3)
 
331
            rho43 = rhoo**F4o3  
 
332
            rho13 = rho43/rhoo
 
333
c
 
334
            Gam =(delrho(n,1,2)*delrho(n,1,2) +
 
335
     &            delrho(n,2,2)*delrho(n,2,2) +
 
336
     &            delrho(n,3,2)*delrho(n,3,2))
 
337
            Gam12 = dsqrt(Gam)
 
338
         if(gam12.lt.dtol) goto 20
 
339
c
 
340
          Eloc = Ax*Rho43
 
341
          x = Gam12/Rho43
 
342
          s = As*x
 
343
          y = s*s
 
344
          PON = mu*y
 
345
          Ffrac = F1-F1/(F1+PON)
 
346
          Fexp  = F1-exp(-PON)
 
347
          fa0 = CxA0
 
348
          fa1 = CxA1 *Ffrac
 
349
          fa2 = CxA2 *Ffrac**F2
 
350
          fa3 = CxA3 *Ffrac**F3
 
351
          fa4 = CxA4 *Ffrac**F4
 
352
          fa5 = CxA5 *Ffrac**F5
 
353
          fb0 = CxB0
 
354
          fb1 = CxB1 *Fexp
 
355
          fb2 = CxB2 *Fexp**F2
 
356
          fb3 = CxB3 *Fexp**F3
 
357
          fb4 = CxB4 *Fexp**F4
 
358
          fb5 = CxB5 *Fexp**F5
 
359
c
 
360
          Fggax = fa0+fa1+fa2+fa3+fa4+fa5 +
 
361
     $            fb0+fb1+fb2+fb3+fb4+fb5
 
362
C
 
363
C     1st derivatives.
 
364
C
 
365
 
 
366
          dElocdR=Ax*F4o3*Rho13
 
367
          dydR = -(F8/F3)*y/Rhoo
 
368
          dydG   = y/Gam
 
369
          dFfracdPON = F1/((F1+PON)**F2)
 
370
          dFexpdPON  = exp(-PON)
 
371
          dfa1dy = CxA1 *mu*dFfracdPON
 
372
          dfa2dy = CxA2 *mu*F2*Ffrac*dFfracdPON
 
373
          dfa3dy = CxA3 *mu*(F3 *Ffrac**F2)*dFfracdPON
 
374
          dfa4dy = CxA4 *mu*(F4 *Ffrac**F3)*dFfracdPON
 
375
          dfa5dy = CxA5 *mu*(F5 *Ffrac**F4)*dFfracdPON
 
376
          dfa1dR = dfa1dy *dydR
 
377
          dfa2dR = dfa2dy *dydR
 
378
          dfa3dR = dfa3dy *dydR
 
379
          dfa4dR = dfa4dy *dydR
 
380
          dfa5dR = dfa5dy *dydR
 
381
          dfa1dG = dfa1dy *dydG
 
382
          dfa2dG = dfa2dy *dydG
 
383
          dfa3dG = dfa3dy *dydG
 
384
          dfa4dG = dfa4dy *dydG
 
385
          dfa5dG = dfa5dy *dydG
 
386
          dfb1dy = CxB1 *mu*dFexpdPON
 
387
          dfb2dy = CxB2 *mu*F2*Fexp*dFexpdPON
 
388
          dfb3dy = CxB3 *mu*(F3 *Fexp**F2)*dFexpdPON
 
389
          dfb4dy = CxB4 *mu*(F4 *Fexp**F3)*dFexpdPON
 
390
          dfb5dy = CxB5 *mu*(F5 *Fexp**F4)*dFexpdPON
 
391
          dfb1dR = dfb1dy *dydR
 
392
          dfb2dR = dfb2dy *dydR
 
393
          dfb3dR = dfb3dy *dydR
 
394
          dfb4dR = dfb4dy *dydR
 
395
          dfb5dR = dfb5dy *dydR
 
396
          dfb1dG = dfb1dy *dydG
 
397
          dfb2dG = dfb2dy *dydG
 
398
          dfb3dG = dfb3dy *dydG
 
399
          dfb4dG = dfb4dy *dydG
 
400
          dfb5dG = dfb5dy *dydG
 
401
c
 
402
          dFggaxdR = dfa1dR+dfa2dR+dfa3dR+dfa4dR+dfa5dR +
 
403
     $               dfb1dR+dfb2dR+dfb3dR+dfb4dR+dfb5dR
 
404
c                     
 
405
          dFggaxdG = dfa1dG+dfa2dG+dfa3dG+dfa4dG+dfa5dG +
 
406
     $               dfb1dG+dfb2dG+dfb3dG+dfb4dG+dfb5dG
 
407
c
 
408
          Ex = Ex + (Eloc*Fggax)*qwght(n)
 
409
          if(ldew) func(n)=func(n)+ Eloc*Fggax
 
410
 
 
411
 
 
412
           Amat(n,2) = Amat(n,2)  +dElocdR*Fggax + Eloc*dFggaxdR
 
413
                      
 
414
           Cmat(n,3)=  Cmat(n,3)  + Eloc*dFggaxdG
 
415
    
 
416
c
 
417
20      continue
 
418
      endif
 
419
      return
 
420
      end
 
421
c
 
422
      Subroutine xc_xsogga_d2()
 
423
      call errquit(' not coded ',0,0)
 
424
      return
 
425
      end
 
426