1 Subroutine dft_mxspin_ovlp(nao,nmo,basis,g_alpha,g_beta,g_tmp) 2 3C$Id$ 4 Implicit none 5#include "errquit.fh" 6 integer basis 7 integer g_s ! overlap 8 integer g_alpha ! alpha eigenvecs [input] 9 integer g_beta ! beta eigenvecs [input] 10 integer g_tmp ! scratch space 11 integer nao ! # of basis functions 12 integer nmo ! # of molecular orbitals 13 14#include "bas.fh" 15#include "cdft.fh" 16#include "mafdecls.fh" 17#include "global.fh" 18#include "tcgmsg.fh" 19#include "msgids.fh" 20#include "stdio.fh" 21#include "util.fh" 22c 23c local 24c 25 integer me,nproc 26c 27 integer i,j,jbig,n,ichunks,nbe,nend 28 integer k_tmpr1, l_tmpr1, k_tmpr2, l_tmpr2, k_tmpi1, l_tmpi1 29 integer g_ss, g_vt, g_u, g_t, k_vals, l_vals, g_alphaT 30 integer nalp, g_ualpha, k_unp, l_unp, l_part, k_part 31 integer nct, l_non, k_non, ind 32c 33 integer ga_create_atom_blocked 34 external ga_create_atom_blocked 35c 36 double precision prodbig, prodtmp 37 double precision alp_thresh, eval, ovlmax, ovl, absovl 38c 39 me=ga_nodeid() 40 nproc=ga_nnodes() 41c 42 g_s = ga_create_atom_blocked(geom, basis, 'AO ovl') 43c 44 if(.not.MA_Push_Get(MT_Dbl,nao,'real_tmp1',l_tmpr1, k_tmpr1)) 45 & call errquit('dft_mxspin_ovlp: cannot allocate real_tmp1',0, 46 & MA_ERR) 47 if(.not.MA_Push_Get(MT_Dbl,nao,'real_tmp2',l_tmpr2, k_tmpr2)) 48 & call errquit('dft_mxspin_ovlp: cannot allocate real_tmp2',0, 49 & MA_ERR) 50 if(.not.MA_Push_Get(MT_Int,nao,'int_tmp1',l_tmpi1, k_tmpi1)) 51 & call errquit('dft_mxspin_ovlp: cannot allocate int_tmp1',0, 52 & MA_ERR) 53c 54 call ga_sync() 55 call ga_zero(g_s) 56 call int_1e_ga(basis,basis,g_s,'overlap',.false.) 57c 58c Compute matrix mult (C_alpha)T * S * C_beta = S` 59c 60 call ga_dgemm('T','N',nmo,nao,nao,1.d0,g_alpha,g_s,0.d0,g_tmp) 61 call ga_dgemm('N','N',nmo,nmo,nao,1.d0,g_tmp,g_beta,0.d0,g_s) 62c 63c 64 if(me.eq.0) then 65 write(LuOut, 9996) 66c call ga_print(g_s) 67 endif 68 if(me.eq.0) then 69 jbig = 1 ! take care of compiler warnings 70 do i = 1, nmo 71c 72c get row of g_s 73c 74 call ga_get(g_s,i,i,1,nmo,DBL_MB(k_tmpr1),1) 75 prodbig=0.0d0 76 do j = 1, nmo 77 prodtmp = abs(dbl_mb(k_tmpr1+j-1)) 78 if(prodtmp.gt.prodbig) then 79 prodbig = prodtmp 80 jbig = j 81 endif 82 enddo 83 dbl_mb(k_tmpr2+i-1) = prodbig 84 int_mb(k_tmpi1+i-1) = jbig 85 enddo 86 ichunks = nmo/10 87 if(ichunks*10.lt.nbf)ichunks = ichunks +1 88 do i = 1, ichunks 89 nbe = 10*(i-1) + 1 90 nend = nbe + 9 91 if(nend.gt.nmo)nend = nmo 92 write(LuOut,9997)(n,n=nbe,nend) 93 write(LuOut,9998)(int_mb(k_tmpi1+n-1),n=nbe,nend) 94 write(LuOut,9999)(dbl_mb(k_tmpr2+n-1),n=nbe,nend) 95 enddo 96 endif 97 9996 format(/,1x,' alpha - beta orbital overlaps ',/, 98 & 1x,' ----------------------------- ',/) 99 9997 format(/,1x,' alpha ',10(1x,i5,1x)) 100 9998 format( 1x,' beta ',10(1x,i5,1x)) 101 9999 format( 1x,'overlap ',10(f7.3),/) 102c 103c 104c ------------------------------------------------------------- 105c Print some useful information 106c about the un-paired alpha orbitals. 107c 108 if(.not.util_print('alpha partner info', print_high)) goto 2001 109 if(noc(2).eq.0) goto 2001 110c 111c Find alpha orbital which best overlaps with each beta 112c orbital, then the remaining alpha orbitals are the un-partnered 113c ones: 114c 115 if(.not.MA_Push_Get(MT_int,noc(2),'partners',l_part, k_part)) 116 & call errquit('dft_mxspin_ovlp: cannot allocate values',0, 117 & MA_ERR) 118c 119c find alpha partners for the beta orbitals: 120c 121 do j = 1,noc(2) 122 ovlmax = 0.0d0 123 do i = 1,noc(1) 124 call ga_get(g_s,i,i,j,j,ovl,1) 125 absovl = dabs(ovl) 126 if (absovl.gt.ovlmax) then 127 ovlmax = absovl 128 int_mb(k_part+j-1) = i 129 endif 130 enddo !i 131 enddo !j 132c 133 if(me.eq.0) then 134 write(luout,*) 135 write(luout,*)'ALPHA/BETA PARTNERS' 136 endif 137c 138 call ga_sync() 139 if(me.eq.0) then 140 ichunks = noc(2)/10 141 if(ichunks*10.lt.noc(2))ichunks = ichunks +1 142 do i = 1, ichunks 143 nbe = 10*(i-1) + 1 144 nend = nbe + 9 145 if(nend.gt.noc(2))nend = noc(2) 146 write(LuOut,9997)(int_mb(k_part+n-1),n=nbe,nend) 147 write(LuOut,9998)(n,n=nbe,nend) 148 enddo 149 write(luout,*) 150 endif 151c 152c determine which alpha orbitals are left over: 153c 154 call ga_sync() 155 if(.not.MA_Push_Get(MT_int,noc(1),'non',l_non, k_non)) 156 & call errquit('dft_mxspin_ovlp: cannot allocate values',0, 157 & MA_ERR) 158c 159 nct=1 160 do 160 i = 1,noc(1) 161 do 150 j = 1,noc(2) 162 if(int_mb(k_part+j-1).EQ.i) goto 160 163c if(nct.le.noc(1)-noc(2)) then 164 int_mb(k_non+nct-1) = i 165c endif 166 150 continue 167 nct = nct + 1 168 160 continue 169c 170 nct=nct-1 171c 172 if(me.eq.0) then 173 if (nct.Ge.1) then 174 write(luout,9995) nct, (int_mb(k_non+n-1),n=1,nct) 175 endif 176 endif 177c 178 9995 format(/,1X,'THERE ARE ',i3,' UN-PARTNERED ALPHA ORBITALS :',20I4) 179c 180c print the un-partnered alpha orbitals 181c 182 if (nct.GE.1) then !only do the rest if there are un-partnered orbitals 183c 184 call ga_sync() 185c 186 call movecs_print_anal(basis,int_mb(k_non),int_mb(k_non) 187 & ,0.15d0,g_alpha,'Alpha Orbitals without Beta Partners', 188 & .false., 0.0 ,.false., 0 , .false., 0 ) 189c 190 if (nct.GE.2) then 191 do i = 2,nct 192 ind = int_mb(k_non+i-1) 193 call movecs_print_anal(basis,ind,ind 194 & ,0.15d0,g_alpha,' ', 195 & .false., 0.0 ,.false., 0 , .false., 0 ) 196 enddo 197 endif 198c 199c Overlap Diagonalization 200c 201c SVD diagonalizes Sab*SabT (Sab is the S' calculated above). 202c Following diagonalization, it is clearer which 203c alpha orbitals do not have beta partners. For alpha orbitals 204c with partners, the overlap is very nearly 1.0. For 205c un-partnered orbitals, the overlap is zero. 206c See J. Chem. Phys. (1967) 47, 1936. 207c 208c 209c if the SVD vectors and overlap aren't desired in the output, 210c don't bother calculating them. 211c 212c calculate Sab*SabT 213c 214 call ga_sync() 215 if(.not.ga_create(mt_dbl,noc(1),noc(1),'SS',0,0,g_ss)) 216 $ call errquit('ga_create failed', g_ss, GA_ERR) 217c 218 call ga_dgemm('N','T',noc(1),noc(1),noc(2),1.d0,g_s,g_s,0.d0,g_ss) 219c 220c call ga_print(g_ss) 221c create arrays needed for SVD 222c 223 if(.not.ga_create(mt_dbl,noc(1),noc(1),'u',0,0,g_u)) 224 $ call errquit('ga_create failed', g_u, GA_ERR) 225c 226 if(.not.ga_create(mt_dbl,noc(1),noc(1),'vt',0,0,g_vt)) 227 $ call errquit('ga_create failed', g_vt, GA_ERR) 228c 229 if(.not.MA_Push_Get(MT_dbl,noc(1),'values',l_vals, k_vals)) 230 & call errquit('dft_mxspin_ovlp: cannot allocate values',0, 231 & MA_ERR) 232 call ga_sync() 233c 234c perform SVD on Sab*SabT to determine unpaired alpha MO's 235c 236 call ga_svd_seq(g_ss,g_u,g_vt,dbl_mb(k_vals)) 237c 238c call ga_print(g_u) 239c call ga_print(g_vt) 240c 241 if (.not. ga_destroy(g_vt)) call errquit 242 & ('dft_mxspin_ovlp: could not destroy g_vt', 0, GA_ERR) 243 if (.not. ga_destroy(g_ss)) call errquit 244 & ('dft_mxspin_ovlp: could not destroy g_ss', 0, GA_ERR) 245c 246 call ga_sync() 247c write(luout,*)'SVD eigenvalues' 248c write(luout,*) (dbl_mb(k_vals+i),i=0,noc(1)-1) 249c call ga_sync() 250c 251c calculate transformed alpha vectors, alphaT 252c 253 if(.not.ga_create(mt_dbl,noc(1),noc(1),'t',0,0,g_t)) 254 $ call errquit('ga_create failed', g_t, GA_ERR) 255c 256 call ga_zero(g_t) 257c 258 call ga_sync() 259 do i = 1,noc(1) 260 call ga_put(g_t,i,i,i,i,dbl_mb(k_vals+i-1),k_vals) 261 enddo 262 call ga_sync() 263c 264c call ga_print(g_t) 265c 266 if(.not.ga_create(mt_dbl,nao,noc(1),'alphaT',0,0,g_alphaT)) 267 $ call errquit('ga_create failed', g_alphaT, GA_ERR) 268c 269 call ga_dgemm('N','N',nao,noc(1),noc(1),1.d0 270 & ,g_alpha,g_u,0.d0,g_alphaT) 271c 272c call ga_print(g_alpha) 273 if (.not. ga_destroy(g_u)) call errquit 274 & ('dft_mxspin_ovlp: could not destroy g_u', 0, GA_ERR) 275c 276c call ga_print(g_alphaT) 277c 278 if (.not. ga_destroy(g_t)) call errquit 279 & ('dft_mxspin_ovlp: could not destroy g_t', 0, GA_ERR) 280c 281c create array containing only alpha MO's which don't 282c have beta partners 283c 284 call ga_sync() 285 nalp = 0 286 alp_thresh = 1.0d-10 287c 288 if(.not.MA_Push_Get(MT_int,noc(1),'unpaired',l_unp, k_unp)) 289 & call errquit('dft_mxspin_ovlp: cannot allocate values',0, 290 & MA_ERR) 291c 292 do i = 1,noc(1) 293 eval = dbl_mb(k_vals+i-1) 294 if (dabs(eval).LT.alp_thresh) then 295 nalp = nalp + 1 296 int_mb(k_unp+i-1) = 1 297 else 298 int_mb(k_unp+i-1) = 0 299 endif 300 enddo 301c 302c 303c write(luout,*) 'paired/unpaired alpha orbitals' 304c write(luout,*) (int_mb(k_unp+i-1),i=1,noc(1)) 305c 306 call ga_sync() 307 if(.not.ga_create(mt_dbl,nao,max(nalp,1),'unp alphaT', 308 $ 0,0,g_ualpha)) 309 $ call errquit('ga_create failed', g_ualpha, GA_ERR) 310c 311 nalp=0 312 do i = 1,noc(1) 313 if(int_mb(k_unp+i-1).EQ.1) then 314 nalp = nalp + 1 315 call ga_copy_patch('N',g_alphaT,1,nbf,i,i 316 & ,g_ualpha,1,nao,nalp,nalp) 317 endif 318 enddo 319 call ga_sync() 320 321c call ga_print(g_ualpha) 322c 323c associate SVD transformed alpha orbitals with the original 324c alpha orbitals: 325c 326 if(me.eq.0) then 327 write(luout,9989) 328 & '==================================================' 329 write(luout,9989) 330 & ' Performing Singular Value Decomposition (SVD) ' 331 write(luout,9989) 332 & 'to diagonalize and maximize the alpha/beta overlap' 333 write(luout,9989) 334 & ' See J. Chem. Phys. (1967) 47, 1936. ' 335 write(luout,9989) 336 write(luout,9989) 337 & 'NOTE: The vector numbering of the SVD transformed ' 338 write(luout,9989) 339 & 'orbitals is different from the original orbitals. ' 340 write(luout,9989) 341 & '==================================================' 342 write(luout,9989) 343 endif 344 9989 format(13x,a51) 345c 346c if(me.eq.0) then 347c if (nalp.GT.1) then 348c write(luout,9990) nalp 349c endif 350c endif 351c 9990 format(/,18x,'THERE ARE',i3,1x,'UN-PARTNERED ALPHA ORBITALS') 352c 353 call movecs_print_anal(basis, 1, nalp, 0.15d0, g_ualpha, 354 & 'Alpha Orb. w/o Beta Partners (after maxim. alpha/beta overlap)', 355 & .false., 0.0 ,.false., 0 , .false., 0 ) 356c 357c print the SVD eigenvalues 358c 359 call ga_sync() 360 if(me.eq.0) then 361 write(LuOut,9994) 362 ichunks = noc(1)/10 363 if(ichunks*10.lt.noc(1))ichunks = ichunks +1 364 do i = 1, ichunks 365 nbe = 10*(i-1) + 1 366 nend = nbe + 9 367 if(nend.gt.noc(1)) nend = noc(1) 368 write(LuOut,9997)(n,n=nbe,nend) 369 write(LuOut,9998)(n,n=nbe,nend) 370 write(LuOut,9999)(dbl_mb(k_vals+n-1),n=nbe,nend) 371 enddo 372 write(luout,*) 373 endif 374c 375 9994 format(/,1x,' SVD maximized alpha - beta orbital overlaps ',/, 376 & 1x,' ------------------------------------------- ',/) 377c 378c 379 if (.not. ga_destroy(g_ualpha)) call errquit 380 & ('dft_mxspin_ovlp: could not destroy g_ualpha', 0, GA_ERR) 381c 382 if (.not. ga_destroy(g_alphaT)) call errquit 383 & ('dft_mxspin_ovlp: could not destroy g_alphaT', 0, GA_ERR) 384c 385 if(.not.MA_Pop_Stack(l_unp)) 386 & call errquit('dft_mxspin_ovlp: cannot pop stack',0, MA_ERR) 387c 388 if(.not.MA_Pop_Stack(l_vals)) 389 & call errquit('dft_mxspin_ovlp: cannot pop stack',0, MA_ERR) 390c 391 endif !if na>nb 392c 393 if(.not.MA_Pop_Stack(l_non)) 394 & call errquit('dft_mxspin_ovlp: cannot pop stack',0, MA_ERR) 395c 396 if(.not.MA_Pop_Stack(l_part)) 397 & call errquit('dft_mxspin_ovlp: cannot pop stack',0, MA_ERR) 398c 399c--------------------------------------------------------------- 4002001 continue 401 if (.not. ga_destroy(g_s)) call errquit 402 & ('dft_mxspin_ovlp: could not destroy g_s', 0, GA_ERR) 403c 404 if(.not.MA_chop_Stack(l_tmpr1)) 405 & call errquit('dft_mxspin_ovlp: cannot pop stack',0, MA_ERR) 406c 407 return 408 end 409