← Back to branch summary

~ubuntu-branches/ubuntu/saucy/nwchem/saucy

« back to all changes in this revision

Viewing changes to src/nwdft/scf_dft_cg/dft_uks_energy.F

  • Committer: Package Import Robot
  • Author(s): Michael Banck, Michael Banck, Daniel Leidert
  • Date: 2012-02-09 20:02:41 UTC
  • mfrom: (1.1.1)
  • Revision ID: package-import@ubuntu.com-20120209200241-jgk03qfsphal4ug2
Tags: 6.1-1
http://bugs.debian.org/632611
https://launchpad.net/bugs/791308
* New upstream release.

[ Michael Banck ]
* debian/patches/02_makefile_flags.patch: Updated.
* debian/patches/02_makefile_flags.patch: Use internal blas and lapack code.
* debian/patches/02_makefile_flags.patch: Define GCC4 for LINUX and LINUX64
  (Closes: #632611 and LP: #791308).
* debian/control (Build-Depends): Added openssh-client.
* debian/rules (USE_SCALAPACK, SCALAPACK): Removed variables (Closes:
  #654658).
* debian/rules (LIBDIR, USE_MPIF4, ARMCI_NETWORK): New variables.
* debian/TODO: New file.
* debian/control (Build-Depends): Removed libblas-dev, liblapack-dev and
  libscalapack-mpi-dev.
* debian/patches/04_show_testsuite_diff_output.patch: New patch, shows the
  diff output for failed tests.
* debian/patches/series: Adjusted.
* debian/testsuite: Optionally run all tests if "all" is passed as option.
* debian/rules: Run debian/testsuite with "all" if DEB_BUILD_OPTIONS
  contains "checkall".

[ Daniel Leidert ]
* debian/control: Used wrap-and-sort. Added Vcs-Svn and Vcs-Browser fields.
  (Priority): Moved to extra according to policy section 2.5.
  (Standards-Version): Bumped to 3.9.2.
  (Description): Fixed a typo.
* debian/watch: Added.
* debian/patches/03_hurd-i386_define_path_max.patch: Added.
  - Define MAX_PATH if not defines to fix FTBFS on hurd.
* debian/patches/series: Adjusted.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/nwdft/scf_dft_cg/dft_uks_energy.F
 
1
      subroutine dft_uks_energy( g_vecs, eone, etwo, exc, enrep, energy,
 
2
     $                           g_grad, nexc )
 
3
      implicit none
 
4
#include "errquit.fh"
 
5
#include "mafdecls.fh"
 
6
#include "global.fh"
 
7
#include "geom.fh"
 
8
#include "cuhf.fh"
 
9
#include "cscf.fh"
 
10
#include "util.fh"
 
11
#include "cscfps.fh"
 
12
#include "rtdb.fh"
 
13
#include "bgj.fh"
 
14
#include "case.fh"
 
15
c     
 
16
c     $Id: dft_uks_energy.F 21176 2011-10-10 06:35:49Z d3y133 $
 
17
c
 
18
      integer g_vecs(2)
 
19
      double precision energy
 
20
      integer g_grad
 
21
c     
 
22
      double precision eone, etwo, enrep, Exc(2)
 
23
      integer nExc
 
24
      integer gtype, grow, gcol
 
25
      integer d(4), f(8), nfock
 
26
      integer g_a_dens, g_a_coul, g_a_exch, g_a_xc
 
27
      integer g_b_dens, g_b_coul, g_b_exch, g_b_xc
 
28
      integer g_hcore
 
29
      integer g_tmp(6),ifock
 
30
      double precision jfac(4), kfac(4), one, zero, mone
 
31
      parameter (one=1.0d0, zero=0.0d0, mone=-1.0d0)
 
32
      double precision e_a_coul, e_a_exch, e_b_coul, e_b_exch,
 
33
     &     e_a_xc, e_b_xc
 
34
      double precision errmaxa, errmaxb
 
35
      integer ga_create_atom_blocked
 
36
      external ga_create_atom_blocked
 
37
      logical odebug
 
38
      logical cphf_uhf
 
39
      logical xc_gotxc
 
40
      external xc_gotxc
 
41
      integer dims(3),chunk(3)
 
42
c
 
43
      double precision xc_hfexch
 
44
      external xc_hfexch
 
45
c
 
46
      data f/8*0/
 
47
c     
 
48
c     Check
 
49
c     
 
50
      odebug = util_print('uks_debug', print_debug)
 
51
      call uhf_jkfac(jfac,kfac)
 
52
      if (.not.cuhf_init_flag)
 
53
     $     call errquit('dft_uks_energy: UKS internal block invalid',0,
 
54
     &       UNKNOWN_ERR)
 
55
      call ga_inquire(g_grad, gtype, grow, gcol)
 
56
      if ((grow.ne.cuhf_vlen).or.(gcol.ne.1))
 
57
     $     call errquit('dft_uks_energy: invalid vector length',0,
 
58
     $                  GA_ERR)
 
59
      cphf_uhf = .false.
 
60
      if (.not. rtdb_get(bgj_get_rtdb_handle(), 
 
61
     &     'cphf_solve:cphf_uhf', mt_log, 1, cphf_uhf)) then
 
62
         cphf_uhf = .false.
 
63
      endif
 
64
c     
 
65
c     Arrays for AO density, coulomb and exchange matrices
 
66
c
 
67
      g_a_coul = ga_create_atom_blocked(geom, basis, 'uks:a coul')
 
68
      g_b_coul = ga_create_atom_blocked(geom, basis, 'uks:b coul')
 
69
      g_a_exch = ga_create_atom_blocked(geom, basis, 'uks:a exch')
 
70
      g_b_exch = ga_create_atom_blocked(geom, basis, 'uks:b exch')
 
71
      if(cphf_uhf .or. xc_gotxc())then
 
72
         g_a_xc   = ga_create_atom_blocked(geom, basis, 'uks:a xc')
 
73
         g_b_xc   = ga_create_atom_blocked(geom, basis, 'uks:b xc')
 
74
      endif
 
75
      g_a_dens = ga_create_atom_blocked(geom, basis, 'uks:a dens')
 
76
      g_b_dens = ga_create_atom_blocked(geom, basis, 'uks:b dens')
 
77
      call ga_zero(g_a_dens)
 
78
      call ga_zero(g_b_dens)
 
79
c
 
80
c     Make the densites and build the fock matrices
 
81
c
 
82
      call ga_dgemm('n', 't', nbf, nbf, nalpha, one, g_vecs(1),
 
83
     $     g_vecs(1), zero, g_a_dens)
 
84
      if (nbeta .gt. 0) then
 
85
         call ga_dgemm('n', 't', nbf, nbf, nbeta, one, g_vecs(2),
 
86
     $        g_vecs(2), zero, g_b_dens)
 
87
      else
 
88
         call ga_zero(g_b_dens)
 
89
      endif
 
90
c
 
91
c     Since UHF can break spatial symmetry by localizing the orbitals
 
92
c     the densities may not be totally symmetric, but since the
 
93
c     Hamiltonian is symmetric contraction with the integrals projects
 
94
c     out the totally symmetric component ... hence we can symmetrize
 
95
c     the densities and exploit symmetry. Compute the max change in any 
 
96
c     element due to symmetrizing and print a warning if it is big.
 
97
c
 
98
c     !! If this is the case then where does the 'force' for symmetry
 
99
c     breaking come from?  Must be missing something?
 
100
c     !! Yes, the symmetrization does not take to irrep of the wave-
 
101
c     function into account. Instead it generates some totally symmetric
 
102
c     density, which is wrong if the wavefunction has a different
 
103
c     irrep.
 
104
c
 
105
      call ga_copy(g_a_dens,g_a_coul)
 
106
      call ga_copy(g_b_dens,g_b_coul)
 
107
      if (oskel) then
 
108
         if (oscfps) call pstat_on(ps_sym_sym)
 
109
         call sym_symmetrize(geom, basis, .true., g_a_dens)
 
110
         if (oscfps) call pstat_off(ps_sym_sym)
 
111
         if (oscfps) call pstat_on(ps_sym_sym)
 
112
         call sym_symmetrize(geom, basis, .true., g_b_dens)
 
113
         if (oscfps) call pstat_off(ps_sym_sym)
 
114
      endif
 
115
      call ga_dadd(one, g_a_dens, mone, g_a_coul, g_a_coul)
 
116
      call ga_dadd(one, g_b_dens, mone, g_b_coul, g_b_coul)
 
117
      call ga_maxelt(g_a_coul, errmaxa)
 
118
      call ga_maxelt(g_b_coul, errmaxb)
 
119
      if (max(errmaxa,errmaxb).gt.1d-4) then
 
120
         if (ga_nodeid().eq.0) then
 
121
            write(6,77) errmaxa,errmaxb
 
122
 77         format(' Warning: spatial symmetry breaking in UKS: ',
 
123
     $           1p,2d9.2)
 
124
            call util_flush(6)
 
125
         endif
 
126
      endif
 
127
c
 
128
      if (odebug) then
 
129
         call ga_print(g_vecs(1))
 
130
         call ga_print(g_vecs(2))
 
131
         call ga_print(g_a_dens)
 
132
         call ga_print(g_b_dens)
 
133
      endif
 
134
 
 
135
      call ga_zero(g_a_coul)
 
136
      call ga_zero(g_b_coul)
 
137
      call ga_zero(g_a_exch)
 
138
      call ga_zero(g_b_exch)
 
139
      if(cphf_uhf .or. xc_gotxc())then
 
140
         call ga_zero(g_a_xc)
 
141
         call ga_zero(g_b_xc)
 
142
      endif
 
143
      d(1) = g_a_dens
 
144
      d(2) = g_a_dens
 
145
      d(3) = g_b_dens
 
146
      d(4) = g_b_dens
 
147
      f(1) = g_a_coul
 
148
      f(2) = g_a_exch
 
149
      f(3) = g_b_coul
 
150
      f(4) = g_b_exch
 
151
      if(cphf_uhf .or. xc_gotxc())then
 
152
         f(5) = g_a_xc
 
153
         f(6) = g_b_xc
 
154
      endif
 
155
c
 
156
c     two extra ga's are passed to fock_2e to get the xc matrix
 
157
c
 
158
      if(cphf_uhf)then
 
159
         nfock = 6
 
160
      else
 
161
         nfock = 4
 
162
      endif
 
163
      call do_riscf (.false.)
 
164
      if (.not.cam_exch) then
 
165
        kfac(1)=0d0
 
166
        kfac(2)=xc_hfexch()
 
167
        kfac(3)=0d0
 
168
        kfac(4)=xc_hfexch()
 
169
        call fock_2e(geom, basis, nfock, jfac, kfac, tol2e,
 
170
     $     oskel, d, f, .false.)
 
171
      else ! for attenuated calculations
 
172
c
 
173
c       calculate the CAM exchange
 
174
c
 
175
        do ifock = 1,nfock
 
176
           g_tmp(ifock) = ga_create_atom_blocked(geom, basis, 'tmp')
 
177
        end do
 
178
        call ga_zero(g_tmp)
 
179
c
 
180
        call case_setflags(.true.)
 
181
        jfac(1)=0d0
 
182
        jfac(2)=0d0
 
183
        jfac(3)=0d0
 
184
        jfac(4)=0d0
 
185
        kfac(1)=0d0
 
186
        kfac(2)=xc_hfexch()
 
187
        kfac(3)=0d0
 
188
        kfac(4)=xc_hfexch()
 
189
        call fock_2e_cam(geom, basis, nfock, jfac, kfac, tol2e, 
 
190
     &     oskel, d, g_tmp, .false., .false.)
 
191
        call ga_dadd(1d0,f,1d0,g_tmp,f)
 
192
c
 
193
c       calculate the full Coulomb
 
194
c
 
195
        call case_setflags(.false.)
 
196
        jfac(1)=1d0
 
197
        jfac(2)=0d0
 
198
        jfac(3)=1d0
 
199
        jfac(4)=0d0
 
200
        kfac(1)=0d0
 
201
        kfac(2)=0d0
 
202
        kfac(3)=0d0
 
203
        kfac(4)=0d0
 
204
        call fock_2e_cam(geom, basis, nfock, jfac, kfac, tol2e, 
 
205
     &     oskel, d, g_tmp, .false., .true.)
 
206
        call ga_dadd(1d0,f,1d0,g_tmp,f)
 
207
c
 
208
c       destroy work space
 
209
        if (.not. ga_destroy(g_tmp)) 
 
210
     &   call errquit('uhf: ga corrupt?',0, GA_ERR)
 
211
      end if
 
212
      call do_riscf (.true.)
 
213
c
 
214
c     do DFT stuff
 
215
c
 
216
      Exc(1) = 0.0d0
 
217
      Exc(2) = 0.0d0
 
218
      if (xc_gotxc()) then
 
219
        call fock_xc(geom, nbf, basis, 4*2, d, f(5), Exc, nExc, .false.)
 
220
      endif
 
221
c
 
222
c     enddo DFT stuff
 
223
c
 
224
      e_a_coul = 0.5d0*
 
225
     $     (ga_ddot(g_a_dens,g_a_coul) + ga_ddot(g_a_dens,g_b_coul))
 
226
      e_b_coul = 0.5d0*
 
227
     $     (ga_ddot(g_b_dens,g_a_coul) + ga_ddot(g_b_dens,g_b_coul))
 
228
      e_a_exch = 0.5d0*ga_ddot(g_a_dens,g_a_exch)
 
229
      e_b_exch = 0.5d0*ga_ddot(g_b_dens,g_b_exch)
 
230
      if (xc_gotxc()) then
 
231
        etwo = e_a_coul + e_b_coul 
 
232
        exc(1) = exc(1) - e_a_exch - e_b_exch
 
233
      else
 
234
        etwo = e_a_coul + e_b_coul - e_a_exch - e_b_exch 
 
235
      endif
 
236
c take this out? it seems that the xc_gotxc part does essentially the
 
237
c same thing apart from the fact that the energy is added on properly.
 
238
      if(cphf_uhf)then
 
239
         e_a_xc = ga_ddot(g_a_dens,g_a_xc)
 
240
         e_b_xc = ga_ddot(g_b_dens,g_b_xc)
 
241
         etwo = etwo + e_a_xc + e_b_xc        
 
242
      endif
 
243
c
 
244
      if (odebug .and. ga_nodeid().eq.0) then
 
245
         write(6,*) ' coulomb energies', e_a_coul, e_b_coul
 
246
         write(6,*) ' exchang energies', e_a_exch, e_b_exch
 
247
         call util_flush(6)
 
248
      endif
 
249
      if (odebug) then
 
250
         call ga_print(g_a_coul)
 
251
         call ga_print(g_a_exch)
 
252
      endif
 
253
c
 
254
c     Form energies and AO fock matrices
 
255
c
 
256
c     Fa (in g_a_coul) = h + J(a) + J(b) - K(a)
 
257
c     Fb (in g_b_coul) = h + J(a) + J(b) - K(b)
 
258
c
 
259
c     E = ((Da + Db)*h + Da*Fa + Db*Fb) / 2
 
260
c     Eone = h * (Da + Db)
 
261
c
 
262
c     2e denotes 2-electron components only
 
263
c
 
264
      call ga_dadd(one, g_a_coul, one, g_b_coul, g_a_coul)
 
265
      call ga_copy(g_a_coul, g_b_coul)
 
266
      call ga_dadd(one, g_a_coul, mone, g_a_exch, g_a_coul)
 
267
      call ga_dadd(one, g_b_coul, mone, g_b_exch, g_b_coul)
 
268
      if(cphf_uhf.or.xc_gotxc())then
 
269
         call ga_dadd(one, g_a_coul, one, g_a_xc, g_a_coul)
 
270
         call ga_dadd(one, g_b_coul, one, g_b_xc, g_b_coul)
 
271
      endif
 
272
c
 
273
c     reuse g_a_exch to hold the 1-e integrals
 
274
c
 
275
      g_hcore = g_a_exch
 
276
      call ga_zero(g_hcore)
 
277
      call int_1e_ga(basis, basis, g_hcore, 'kinetic', oskel)
 
278
      call int_1e_ga(basis, basis, g_hcore, 'potential', oskel)
 
279
      eone = 
 
280
     $     (ga_ddot(g_a_dens,g_hcore) + ga_ddot(g_b_dens,g_hcore))
 
281
      call ga_dadd(one, g_hcore, one, g_a_coul, g_a_coul)
 
282
      call ga_dadd(one, g_hcore, one, g_b_coul, g_b_coul)
 
283
c
 
284
      if (oskel) then
 
285
         if (oscfps) call pstat_on(ps_sym_sym)
 
286
         call sym_symmetrize(geom, basis, .false., g_a_coul)
 
287
         if (oscfps) call pstat_off(ps_sym_sym)
 
288
         if (oscfps) call pstat_on(ps_sym_sym)
 
289
         call sym_symmetrize(geom, basis, .false., g_b_coul)
 
290
         if (oscfps) call pstat_off(ps_sym_sym)
 
291
      endif
 
292
c
 
293
      if (odebug) then
 
294
         call ga_print(g_a_coul)
 
295
         call ga_print(g_b_coul)
 
296
      endif
 
297
c
 
298
c     Transform the Fock matrices to the MO basis using g_a_dens
 
299
c     for scratch
 
300
c
 
301
      call two_index_transf(g_a_coul, g_vecs(1), g_vecs(1),
 
302
     $     g_a_dens, cuhf_g_falpha)
 
303
      call two_index_transf(g_b_coul, g_vecs(2), g_vecs(2),
 
304
     $     g_a_dens, cuhf_g_fbeta)
 
305
c
 
306
      if (odebug) then
 
307
         call ga_print(cuhf_g_falpha)
 
308
         call ga_print(cuhf_g_fbeta)
 
309
      endif
 
310
c
 
311
c     Free up dead global arrays
 
312
c
 
313
      if (.not. ga_destroy(g_a_dens)) call errquit('uks_e: destroy',0,
 
314
     &       GA_ERR)
 
315
      if (.not. ga_destroy(g_b_dens)) call errquit('uks_e: destroy',0,
 
316
     &       GA_ERR)
 
317
      if(cphf_uhf .or. xc_gotxc())then
 
318
         if (.not. ga_destroy(g_a_xc)) call errquit('uks_e: destroy',0,
 
319
     &       GA_ERR)
 
320
         if (.not. ga_destroy(g_b_xc)) call errquit('uks_e: destroy',0,
 
321
     &       GA_ERR)
 
322
      endif
 
323
      if (.not. ga_destroy(g_a_exch)) call errquit('uks_e: destroy',0,
 
324
     &       GA_ERR)
 
325
      if (.not. ga_destroy(g_b_exch)) call errquit('uks_e: destroy',0,
 
326
     &       GA_ERR)
 
327
      if (.not. ga_destroy(g_a_coul)) call errquit('uks_e: destroy',0,
 
328
     &       GA_ERR)
 
329
      if (.not. ga_destroy(g_b_coul)) call errquit('uks_e: destroy',0,
 
330
     &       GA_ERR)
 
331
c
 
332
c     extract the gradient
 
333
c
 
334
      call uhf_get_grad(g_grad)
 
335
c
 
336
      if (odebug) call ga_print(g_grad)
 
337
c
 
338
      if (.not. geom_nuc_rep_energy(geom, enrep))
 
339
     $     call errquit('dft_uks_energy: no repulsion energy?', 0,
 
340
     $                  GEOM_ERR)
 
341
      energy = eone + etwo + enrep
 
342
      if (xc_gotxc()) then
 
343
        energy = energy + Exc(1) + Exc(2)
 
344
      endif
 
345
c
 
346
      if (odebug .and. ga_nodeid().eq.0) then
 
347
         write(6,*) ' eone, etwo, enrep, energy ',
 
348
     $        eone, etwo, enrep, energy
 
349
      endif
 
350
c
 
351
      end