1      subroutine rohf_canon(oaufbau, oprint)
2C$Id$
3      implicit none
4#include "errquit.fh"
5#include "mafdecls.fh"
6#include "tcgmsg.fh"
7#include "global.fh"
8#include "cscfps.fh"
9#include "cscf.fh"
10#include "crohf.fh"
11      logical oaufbau
12      logical oprint
13c
14      integer g_u, g_fock, g_tmp
15      double precision one, zero
16      data one, zero/1.d0, 0.d0/
17c
18c     This routine assumes that rohf_energy/rohf_fock have been called
19c     so that the contents of /crohf/ are current.
20c
21c     Diagonalize the ROHF 'Fock' matrix
22c
23c     If (oaufbau)
24c        diagonalize the whole thing and allow mixing of closed-open-virt
25c     else
26c        diagonalize separately the cloed-closed, open-open, and
27c        virt-virt parts
28c
29c     Transform Fock matrices and MO coefficients into the new canonical basis
30c
31      if (oscfps) call pstat_on(ps_diag)
32c
33      if (.not. ga_create(MT_DBL, nmo, nmo, 'rohf_canon: u',
34     $     32, 32, g_u)) call errquit('rohf_canon: ga failed for u', 0,
35     &       GA_ERR)
36      call ga_zero(g_u)         ! ESSENTIAL FOR SYMM BLOCKED SUBSPACE CANON
37      if (.not. ga_create(MT_DBL, nmo, nmo, 'rohf_canon: fock',
38     $     32, 32, g_fock)) call errquit
39     $     ('rohf_canon: ga failed for fock', 0, GA_ERR)
40c
41      call rohf_get_fock(g_fock)
42c
43      if (oaufbau) then
44#if defined(PARALLEL_DIAG)
45#ifdef SCALAPACK
46         call ga_pdsyev(g_fock, g_u, dbl_mb(k_eval), 0)
47#else
48         call ga_diag_std(g_fock, g_u, dbl_mb(k_eval))
49#endif
50#else
51         call ga_diag_std_seq(g_fock, g_u, dbl_mb(k_eval))
52#endif
53      else
54c
55c     closed-closed piece
56c
57         if (nclosed .gt. 0) call rohf_canon_subspace
58     $        (g_fock, g_u, dbl_mb(k_eval), 1, nclosed,int_mb(k_irs))
59c
60c     open-open piece
61c
62         if (nopen .gt. 0) call rohf_canon_subspace
63     $        (g_fock, g_u, dbl_mb(k_eval+nclosed),
64     $        nclosed+1, nclosed+nopen, int_mb(k_irs))
65c
66c     virt-virt piece
67c
68         if (nmo-nclosed-nopen .gt. 0)
69     $        call rohf_canon_subspace(g_fock, g_u,
70     $        dbl_mb(k_eval+nclosed+nopen),
71     $        nclosed+nopen+1, nmo, int_mb(k_irs))
72      endif
73c
74      if (oscfps) call pstat_off(ps_diag)
75c
76      call movecs_fix_phase(g_u)
77c
78c     Apply rotation to orbitals and fock matrix
79c
80*ga:1:0
81      if (.not. ga_duplicate(g_movecs, g_tmp, 'rohf_canon: tmp'))
82     $     call errquit
83     $     ('rohf_canon: ga_dup for tmp', 0, GA_ERR)
84c
85      call ga_copy(g_movecs, g_tmp)
86      call ga_dgemm('n', 'n', nbf, nmo, nmo, one, g_tmp, g_u,
87     $     zero, g_movecs)
88c
89      if (nbf .ne. nmo) then
90         if (.not. ga_destroy(g_tmp))
91     $        call errquit('rohf_canon: ga_destroy?',0, GA_ERR)
92         if (.not. ga_create(MT_DBL, nmo, nmo, 'rohf_canon: tmp',
93     $        32, 32, g_tmp)) call errquit
94     $        ('rohf_canon: ga failed for tmp', 0, GA_ERR)
95      endif
96c
97      call ga_dgemm('n', 'n', nmo, nmo, nmo, one, crohf_g_fcv, g_u,
98     $     zero, g_tmp)
99      call ga_dgemm('t', 'n', nmo, nmo, nmo, one, g_u, g_tmp,
100     $     zero, crohf_g_fcv)
101c
102      if (nopen .gt. 0) then
103         call ga_dgemm('n', 'n', nmo, nmo, nmo, one, crohf_g_fcp, g_u,
104     $        zero, g_tmp)
105         call ga_dgemm('t', 'n', nmo, nmo, nmo, one, g_u, g_tmp,
106     $        zero, crohf_g_fcp)
107         call ga_dgemm('n', 'n', nmo, nmo, nmo, one, crohf_g_fpv, g_u,
108     $        zero, g_tmp)
109         call ga_dgemm('t', 'n', nmo, nmo, nmo, one, g_u, g_tmp,
110     $        zero, crohf_g_fpv)
111      endif
112c
113      if (oprint .and. ga_nodeid().eq.0) then
114         write(6,*)
115         write(6,*)
116         call util_print_centered(6, 'Eigenvalues', 20, .true.)
117         call output(dbl_mb(k_eval), 1, min(nclosed+nopen+5,nmo),
118     $        1, 1, nmo, 1, 1)
119         call util_flush(6)
120      endif
121c
122      if (.not. ga_destroy(g_u))
123     $     call errquit('rohf_canon: destroy', 0, GA_ERR)
124      if (.not. ga_destroy(g_fock))
125     $     call errquit('rohf_canon: destroy', 0, GA_ERR)
126      if (.not. ga_destroy(g_tmp))
127     $     call errquit('rohf_canon: destroy', 0, GA_ERR)
128c
129      end
130      subroutine rohf_canon_subspace(g_fock, g_u, evals, ilo, ihi,
131     $     irs)
132      implicit none
133#include "errquit.fh"
134#include "global.fh"
135#include "mafdecls.fh"
136      integer g_fock, g_u, ilo, ihi, irs(*)
137      double precision evals(*)
138c
139      integer k_map,l_map
140      integer k_work,l_work,k_tmp,l_tmp
141      integer sizew
142c
143c     map dim ihi-ilo+1
144c
145!     sizew=ihi-ilo+1
146      sizew=ihi
147      if (.not.ma_push_get(MT_int,sizew,'map',l_map,k_map))
148     &   call errquit('rohf_can: cannot allocate map',0, MA_ERR)
149      if (.not.ma_push_get(MT_Dbl,sizew,'work',l_work,k_work))
150     &   call errquit('rohf_can: cannot allocate work',0, MA_ERR)
151      if (.not.ma_push_get(MT_Dbl,sizew,'tmp',l_tmp,k_tmp))
152     &   call errquit('rohf_can: cannot allocate tmp',0, MA_ERR)
153      call rohf_canon_sub0(g_fock, g_u, evals, ilo, ihi,
154     $     irs, int_mb(k_map),
155     ,     dbl_mb(k_work), dbl_mb(k_tmp))
156      if (.not.ma_chop_stack(l_map))
157     &   call errquit('rohf_can: cannot pop stack',1, MA_ERR)
158      return
159      end
160c
161      subroutine rohf_canon_sub0(g_fock, g_u, evals, ilo, ihi,
162     $     irs, map, work, tmp)
163      implicit none
164#include "errquit.fh"
165#include "global.fh"
166#include "mafdecls.fh"
167      integer g_fock, g_u, ilo, ihi, irs(*)
168      double precision evals(*)
169c
170      integer map(*)        ! Should be dynamically allocated
171      double precision work(*)
172      double precision tmp(*)
173c
174      integer n, i, j, nirs, irrep
175      integer g_v
176      integer g_tmp
177c
178      nirs = 0
179      do i = ilo, ihi
180         nirs = max(irs(i),nirs)
181      enddo
182c
183      do irrep = 1, nirs
184         n = 0
185         do i = ilo, ihi
186            if (irs(i) .eq. irrep) then
187               n = n + 1
188               map(n) = i
189            endif
190         enddo
191         if (n .gt. 0) then
192c
193*ga:1:0
194            if (.not. ga_create(MT_DBL, n, n, 'rohf_canon: tmp',
195     $           0, 0, g_tmp)) call errquit
196     $           ('rohf_canon: ga failed for tmp', 0, GA_ERR)
197*ga:1:0
198            if (.not. ga_create(MT_DBL, n, n, 'rohf_canon: v',
199     $           n, 0, g_v)) call errquit
200     $           ('rohf_canon: ga failed for v', 0, GA_ERR)
201c
202            call ga_sync
203            do i = 1+ga_nodeid(),n,ga_nnodes()
204               call ga_get(g_fock,ilo,ihi,map(i),map(i),work(ilo),
205     I              ihi-ilo+11)
206               do j = 1, n
207                  tmp(j) = work(map(j))
208               enddo
209               call ga_put(g_tmp,1, n, i, i, tmp, n)
210            enddo
211            call ga_sync
212c
213#if defined(PARALLEL_DIAG)
214#ifdef SCALAPACK
215      call ga_pdsyev(g_tmp, g_v, tmp, 0)
216#else
217      call ga_diag_std(g_tmp, g_v, tmp)
218#endif
219#else
220      call ga_diag_std_seq(g_tmp, g_v, tmp)
221#endif
222            call ga_sync
223            do i = 1, n
224               evals(map(i)-ilo+1) = tmp(i)
225            enddo
226            do i = 1+ga_nodeid(),n,ga_nnodes()
227               call ga_get(g_v,1, n, i, i, tmp, n)
228               call dfill((ihi-ilo+1), 0.0d0, work(ilo), 1)
229               do j = 1, n
230                  work(map(j)) = tmp(j)
231               enddo
232               call ga_acc(g_u,ilo,ihi,map(i),map(i),work(ilo),
233     I              ihi-ilo+1,1.0d0)
234            enddo
235            call ga_sync
236            if (.not. (ga_destroy(g_tmp) .and. ga_destroy(g_v)))
237     $           call errquit('rohf_canon: ga_destroy ?', 0, GA_ERR)
238         endif
239      enddo
240c
241      if (nirs .eq. 1) return
242c
243c     Now we have diagonalized each symmetry block we need to reorder
244c     so that have Auf Bau ordering within the entire block
245c
246      do i = ilo, ihi
247	 map(i) = i
248	 do j = ilo, i-1
249	    if (evals(map(i)-ilo+1).lt.evals(map(j)-ilo+1)) then
250	       n = map(i)
251	       map(i) = map(j)
252	       map(j) = n
253            endif
254         enddo
255      enddo
256      do i = ilo, ihi
257	 if (map(i) .ne. i) goto 100
258      enddo
259      return
260c
261100   continue
262*      do i = ilo, ihi
263*         if (map(i) .ne. i) write(6,*) ' FLIP ', i, map(i)
264*      enddo
265      n = ihi - ilo + 1
266*ga:1:0
267      if (.not. ga_create(MT_DBL, n, n, 'rohf_canon: tmp',
268     $        n, 0, g_tmp)) call errquit
269     $        ('rohf_canon: ga failed for tmp', 0, GA_ERR)
270      do i = ilo+ga_nodeid(),ihi,ga_nnodes()
271         call ga_get(g_u,ilo,ihi, map(i), map(i), tmp(ilo), ihi-ilo+1)
272         call ga_put(g_tmp,1,n,i-ilo+1,i-ilo+1,tmp(ilo),n)
273      enddo
274      call ga_sync
275      do i = ilo+ga_nodeid(),ihi,ga_nnodes()
276         call ga_get(g_tmp,1,n,i-ilo+1,i-ilo+1,work(ilo),n)
277         call ga_put(g_u,ilo,ihi, i, i, work(ilo), ihi-ilo+1)
278      enddo
279*      write(6,*) ' input '
280*      call output(evals,ilo,ihi,1,1,ihi,1,1)
281      do i = ilo, ihi
282	 tmp(i) = evals(map(i)-ilo+1)
283	 work(i)= irs(map(i))
284      enddo
285      do i = ilo, ihi
286	 evals(i-ilo+1) = tmp(i)
287	 irs(i)   = work(i)
288      enddo
289*      write(6,*) ' output '
290*      call output(evals,ilo,ihi,1,1,ihi,1,1)
291c
292      if (.not. ga_destroy(g_tmp)) call errquit('rohf_canon: ga?',0,
293     &       GA_ERR)
294c
295      end
296