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