1c
2c     Determine the HOMO-LUMO gap
3c
4      subroutine calc_homolumogap(k_eval,nelec,rlshift,homo,lumo,gap)
5c
6      implicit none
7c
8#include "errquit.fh"
9#include "stdio.fh"
10#include "mafdecls.fh"
11c
12      integer k_eval(2)      ! Pointers to the eigenvalue arrays. Only (1) is used
13      integer nelec
14      double precision rlshift
15      double precision homo
16      double precision lumo
17      double precision gap
18c
19c     Determine the HOMO, LUMO and gap
20c
21      homo = Dbl_MB(k_eval(1)+nelec-1)   ! Extract the contents of the eigenvalue array
22      lumo = Dbl_MB(k_eval(1)+nelec)
23      gap = min(gap, (lumo-homo-rlshift))
24c
25      return
26      end
27c
28c     Transform Fock To Orthonormal
29c
30      subroutine trans_fock_to_ortho(g_tmp2,nbf_mo,g_sm12,g_fockso)
31c
32      implicit none
33c
34#include "errquit.fh"
35#include "global.fh"
36#include "consts.fh"
37c
38      integer g_tmp2          ! temp scratch array
39      integer nbf_mo          ! number of molecular orbitals
40      integer g_sm12          ! S^(-1/2)
41      integer g_fockso(2)     ! 1: real, 2: imag
42c
43c     Real part of the Fock matrix
44c
45      call ga_zero(g_tmp2)
46      call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one,
47     &                 g_sm12, g_fockso(1), zero, g_tmp2)
48      call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
49     &                 g_tmp2, g_sm12, zero, g_fockso(1))
50c
51c     Imag part of the Fock matrix
52c
53      call ga_zero(g_tmp2)
54      call ga_dgemm('T', 'N', nbf_mo, nbf_mo, nbf_mo, one,
55     &                 g_sm12, g_fockso(2), zero, g_tmp2)
56      call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
57     &                 g_tmp2, g_sm12, zero, g_fockso(2))
58c
59      return
60      end
61c
62c     Diagonalize complex Fock Matrix
63c
64      subroutine diag_fock(nbf_mo,ia,g_fockso,ibuff,g_moso,iwork,irwork,
65     &     k_eval,trace,llwork,info)
66c
67      implicit none
68c
69#include "errquit.fh"
70#include "global.fh"
71#include "mafdecls.fh"
72#include "stdio.fh"
73#include "util.fh"
74c
75      integer nbf_mo
76      integer ia
77      integer g_fockso(2)
78      integer ibuff
79      integer g_moso(2)
80      integer iwork
81      integer irwork
82      integer k_eval(2)
83      double precision trace
84      integer llwork
85      integer info
86c
87      integer i,j,i1
88      double precision ddot
89      external ddot
90c
91c     Prepare arrays for diagonalization
92c
93      trace = 0.d0
94      do i = 1, nbf_mo
95         do j = 1, nbf_mo
96            DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=dcmplx(0.0, 0.0)
97         enddo
98      enddo
99      do i = 1, nbf_mo
100         call ga_get(g_fockso(1), 1,i, i,i, dbl_mb(ibuff),1)
101         do j=1,i
102            DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=
103     =           dcmplx(dbl_mb(ibuff+j-1),0d0)
104         enddo
105         call ga_get(g_fockso(2), 1,i, i,i, dbl_mb(ibuff),1)
106         do j=1,i
107            DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))=
108     $               DCpl_mb(ia+(nbf_mo)*(i-1)+(j-1))
109     $           +dcmplx(0d0,dbl_mb(ibuff+j-1))
110         enddo
111      enddo
112      call ga_zero(g_moso(1))
113      call ga_zero(g_moso(2))
114c
115c     Call the diagonalizer (complex diagonalizer)
116c
117      call zheev( 'V', 'U', nbf_mo, DCpl_mb(ia), nbf_mo,
118     $            Dbl_mb(k_eval(1)),
119     $            DCpl_mb(iwork), llwork, Dbl_mb(irwork), info )
120      do i = 1, nbf_mo
121         do j = 1, nbf_mo
122            dbl_mb(ibuff+j-1)=0.0d0
123            dbl_mb(ibuff+j-1)=dble(DCpl_mb(ia+nbf_mo*(i-1)+(j-1)))
124         enddo
125         i1=i
126         call ga_put(g_moso(1),1,nbf_mo,i1,i1,dbl_mb(ibuff),1)
127         trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1)
128         do j = 1, nbf_mo
129            dbl_mb(ibuff+j-1)=0.0d0
130            dbl_mb(ibuff+j-1)=
131     $             dimag(dcmplx(DCpl_mb(ia+nbf_mo*(i-1)+(j-1))))
132         enddo
133         i1=i
134         call ga_put(g_moso(2),1,nbf_mo,i1,i1,dbl_mb(ibuff),1)
135         trace = ddot(nbf_mo,dbl_mb(ibuff),1,dbl_mb(ibuff),1)
136      enddo
137c
138      return
139      end
140c
141c     Back-transform eigenvectors with S^-1/2.
142c
143      subroutine trans_vec_to_ao(nbf_mo,g_sm12,g_moso,g_fockso)
144c
145      implicit none
146c
147#include "errquit.fh"
148#include "global.fh"
149#include "stdio.fh"
150#include "consts.fh"
151c
152      integer nbf_mo
153      integer g_sm12      ! contains S^-1/2
154      integer g_moso(2)   ! MO vecs 1: real, 2: imag
155      integer g_fockso(2) ! being used as scratch
156c
157c     Back-transform eigenvectors with S^-1/2.
158c
159      call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
160     &     g_sm12, g_moso(1), zero, g_fockso(1))
161      call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
162     &     g_sm12, g_moso(2), zero, g_fockso(2))
163c
164c     Transfer into MO vec arrays
165c
166      call ga_zero(g_moso(1))
167      call ga_zero(g_moso(2))
168      call ga_copy(g_fockso(1), g_moso(1))
169      call ga_copy(g_fockso(2), g_moso(2))
170c
171      return
172      end
173c
174c     Calculate S powers: S,S^(-1/2),S^(1/2),S^(-1)
175c
176      subroutine calc_s_powers(g_scr,g_tmp,nbf_ao,toll_s,svals,g_svecs,
177     &   g_sp1,g_sm12,g_sp12,g_sm1)
178c
179      implicit none
180c
181#include "errquit.fh"
182#include "global.fh"
183#include "stdio.fh"
184#include "consts.fh"
185c
186      integer iw
187      integer g_scr
188      integer g_tmp
189      integer nbf_ao
190      double precision toll_s
191      double precision svals(*)
192      integer g_svecs
193      integer g_sp1
194      integer g_sm12
195      integer g_sp12
196      integer g_sm1
197c
198c     Calculate S^1: g_sp1
199c
200      iw = 1
201      call ga_zero(g_scr)
202      call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
203     &     g_tmp, nbf_ao, iw)
204      call ga_zero(g_sp1)
205      call ga_fock_sf(g_scr, g_sp1, nbf_ao) ! Map onto larger array
206c
207c     Calculate S^-1/2: g_sm12
208c
209      iw = 2
210      call ga_zero(g_scr)
211      call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
212     &     g_tmp, nbf_ao, iw)
213      call ga_zero(g_sm12)
214      call ga_fock_sf(g_scr, g_sm12, nbf_ao) ! Map onto larger array
215c
216c     Calculate S^+1/2: g_sp12
217c
218      iw = 3
219      call ga_zero(g_scr)
220      call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
221     &     g_tmp, nbf_ao, iw)
222      call ga_zero(g_sp12)
223      call ga_fock_sf(g_scr, g_sp12, nbf_ao) ! Map onto larger array
224c
225c     Calculate S^-1: g_sm1
226c
227      iw = 4
228      call ga_zero(g_scr)
229      call diis_bld12_so(toll_s, svals, g_svecs, g_scr,
230     &     g_tmp, nbf_ao, iw)
231      call ga_zero(g_sm1)
232      call ga_fock_sf(g_scr, g_sm1, nbf_ao) ! Map onto larger array
233c
234      return
235      end
236c
237c     Print energies
238c
239      subroutine print_energies(etnew,enuc,ecore,ecoul,exc,
240     &   nexc,rho_n,dft_time)
241c
242      implicit none
243c
244#include "global.fh"
245#include "errquit.fh"
246#include "stdio.fh"
247c
248      double precision etnew
249      double precision enuc
250      double precision ecore
251      double precision ecoul
252      double precision exc(2)
253      double precision nexc
254      double precision rho_n
255      double precision dft_time
256c
257c     Print
258c
259      if (nexc.le.1) then
260       write(luout,222)etnew+enuc,ecore,ecoul,exc(1),enuc
261      else
262       write(luout,223)etnew+enuc,ecore,ecoul,exc(1),exc(2),enuc
263      end if
264      write(luout,2222) rho_n
265      write(luout,2223) dft_time
266c
267 222  format(//
268     &     '      Total SO-DFT energy =', f22.12/
269     &     '      One electron energy =', f22.12/
270     &     '           Coulomb energy =', f22.12/
271     &     '    Exchange-Corr. energy =', f22.12/
272     &     ' Nuclear repulsion energy =', f22.12/)
273c
274 223  format(//
275     &     '      Total SO-DFT energy =', f22.12/
276     &     '      One electron energy =', f22.12/
277     &     '           Coulomb energy =', f22.12/
278     &     '          Exchange energy =', f22.12/
279     &     '       Correlation energy =', f22.12/
280     &     ' Nuclear repulsion energy =', f22.12/)
281c
282 2222 format(' Numeric. integr. density =', f22.12/)
283 2223 format('     Total iterative time =', f9.1,'s'//)
284c
285      return
286      end
287c
288c     Level shifting is implemented here (similarity
289c     transformation before standard eigensolver).  Note,
290c     levelshifting is appropriate once a transformation
291c     is available which makes the resulting Fock matrix
292c     diagonally dominant, e.g., in an approximate MO basis.
293c     Also note, there are many matrix multiplies with S^+-1/2
294c     which are redundant if one is sure that the former basis
295c     is orthonormal.
296c
297      subroutine levelshift_fock(nbf_mo,ntotocc,g_tmp2,g_sp12,g_moso,
298     &   g_scr2,g_fockso,rlshift)
299c
300      implicit none
301c
302#include "global.fh"
303#include "errquit.fh"
304#include "stdio.fh"
305#include "consts.fh"
306c
307         integer nbf_mo
308         integer ntotocc
309         integer g_tmp2
310         integer g_sp12
311         integer g_moso(2)
312         integer g_scr2
313         integer g_fockso(2)
314         double precision rlshift
315c
316         integer j
317         integer me
318         integer nproc
319c
320c        Preliminaries
321c
322         me = ga_nodeid()
323         nproc = ga_nnodes()
324c
325c        Build a matrix which is diagonal in the "MO" rep,
326c        back-transform, and shift the current Fock matrix
327c
328c        Use S^+1/2 (g_sp12)  * old movecs (as a transform).
329c
330c        Real part
331         call ga_zero(g_tmp2)
332         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
333     &                 g_sp12, g_moso(1), zero, g_tmp2)
334         call ga_copy(g_tmp2,  g_moso(1))
335
336c        Imag part
337         call ga_zero(g_tmp2)
338         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
339     &                 g_sp12, g_moso(2), zero, g_tmp2)
340         call ga_copy(g_tmp2,  g_moso(2))
341c
342c        Build diagonal matrix for the shift
343         call ga_zero(g_tmp2)
344         do j = nTotOcc+1+me, nbf_mo, nproc
345            call ga_put(g_tmp2, j, j, j, j, rlshift, 1)
346         enddo
347c
348c        Transform this into "AO" basis and add to current
349c        Fock matrix
350c
351c        Real part
352         call ga_zero(g_scr2)  ! used as a work area
353         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
354     &                 g_moso(1), g_tmp2, zero, g_scr2)
355         call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
356     &                 g_scr2, g_moso(1), one, g_fockso(1))
357         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
358     &                 g_moso(2), g_tmp2, zero, g_scr2)
359         call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
360     &                 g_scr2, g_moso(2), one, g_fockso(1))
361c
362c        Imag part
363         call ga_zero(g_scr2)  ! used as a work area
364         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
365     &                 g_moso(1), g_tmp2, zero, g_scr2)
366         call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, mone,
367     &                 g_scr2, g_moso(2), one, g_fockso(2))
368         call ga_dgemm('N', 'N', nbf_mo, nbf_mo, nbf_mo, one,
369     &                 g_moso(2), g_tmp2, zero, g_scr2)
370         call ga_dgemm('N', 'T', nbf_mo, nbf_mo, nbf_mo, one,
371     &                 g_scr2, g_moso(1), one, g_fockso(2))
372c
373      return
374      end
375c $Id$
376