1 subroutine movecs_anal_so(basis, ilo, ihi, thresh, 2 $ g_vecs, title, 3 $ oevals, evals, oirs, irs, oocc, occ) 4* 5* $Id$ 6* 7 implicit none 8#include "errquit.fh" 9#include "mafdecls.fh" 10#include "bas.fh" 11#include "global.fh" 12#include "inp.fh" 13#include "cscfps.fh" 14c 15 integer basis 16 integer ilo, ihi ! [input] Range of vectors to print 17 double precision thresh ! [input] Print coeffs with absval >= thresh 18 double precision norm, normn 19 integer g_vecs(2) 20 character*(*) title 21 logical oevals ! [input] If true print eigenvalues 22 double precision evals(*) 23 logical oirs ! [input] If true print irreps 24 integer irs(*) 25 logical oocc ! [input] If true print occupations 26 double precision occ(*) 27 logical or2 ! [input] If true print orbital center and r^2 28c 29c Print a summary of the MO vectors in the specified range. 30c 31 integer l_vecsre, k_vecsre, i, j, k_tags, l_tags, k_list, l_list 32 integer l_vecsim, k_vecsim 33 integer n, k, m, klo, khi, ibuf 34 character buf*80 35 double precision cur_thresh 36c 37 integer type, nbf, nmo 38 integer len_r_r2, l_r_r2, k_r_r2 39 integer maxop, maxireps, geom 40 parameter (maxop = 120, maxireps=20) 41 integer nop, nir, class_dim(maxireps) 42 character*8 zir(maxireps), zclass(maxireps) 43 character*20 zname 44 double precision chars(maxireps*maxireps) 45 logical sym_char_table_so 46 external sym_char_table_so 47c 48 or2 = .true. 49 if (.not. bas_cando_mpoles(basis)) or2 = .false. 50c 51 call ga_sync 52c call ga_summarize(.true.) 53 if (oscfps) call pstat_on(ps_moanal) 54 call ga_inquire(g_vecs, type, nbf, nmo) 55 if (or2) then 56c 57c local array for x, y, z, x^2, y^2, z^2 and r^2 for each MO 58c 59 len_r_r2 = 7*nmo 60 if (.not. ma_push_get(mt_dbl, len_r_r2, 'nmo:r_and_r2', 61 & l_r_r2, k_r_r2))call errquit 62 & ('movecs_anal_so: cannot allocate r_r2', len_r_r2, MA_ERR) 63 call so_r_and_r2(basis, g_vecs, dbl_mb(k_r_r2)) 64 endif 65c 66 if (ga_nodeid() .eq. 0) then 67 if (.not. bas_numbf(basis, nbf/2)) call errquit 68 $ ('movecs_anal: basis bad?',basis, BASIS_ERR) 69 if (.not. ma_push_get(mt_dbl,nbf,'movecs_anal1',l_vecsre, 70 & k_vecsre)) 71 $ call errquit('movecs_anal: ma 1 failed', nbf, MA_ERR) 72 if (.not. ma_push_get(mt_dbl,nbf,'movecs_anal1',l_vecsim, 73 & k_vecsim)) 74 $ call errquit('movecs_anal: ma 1.1 failed', nbf, MA_ERR) 75 if (.not. ma_push_get(mt_int,nbf,'movecs_anal2',l_list,k_list)) 76 $ call errquit('movecs_anal: ma 2 failed', nbf, MA_ERR) 77 if (.not. ma_push_get(mt_byte,nbf*16,'movecs_anal3', 78 $ l_tags,k_tags)) 79 $ call errquit('movecs_anal: ma 3 failed', nbf*16, MA_ERR) 80c 81 if (oirs) then 82 if (.not. bas_geom(basis, geom)) call errquit 83 $ ('movecs_anal: bas geom', basis, BASIS_ERR) 84 call sym_group_name(geom, zname) 85 if (.not. sym_char_table_so(zname, nop, nir, class_dim, 86 $ zir, zclass, chars)) call errquit 87 $ ('movecs_anal: no char table available ',geom, 88 & GEOM_ERR) 89 endif 90c 91 call bas_vec_info(basis, byte_mb(k_tags)) 92 call bas_vec_info(basis, byte_mb(k_tags+16*nbf/2)) 93c 94 write(6,*) 95 call util_print_centered(6,title, 40, .true.) 96 write(6,*) 97c 98 99 do i = ilo, ihi 100c 101 call ga_get(g_vecs(1), 1, nbf, i, i, dbl_mb(k_vecsre), 1) 102 call ga_get(g_vecs(2), 1, nbf, i, i, dbl_mb(k_vecsim), 1) 103c 104c Identify significant coefficients and sort by size 105c 106 n = 0 107 cur_thresh = thresh 108 111 do j = 0, nbf-1 109 norm=sqrt(dbl_mb(k_vecsre+j)**2+dbl_mb(k_vecsim+j)**2) 110 if (norm.ge.cur_thresh) then 111 int_mb(k_list + n) = j 112 n = n + 1 113 endif 114 enddo 115 if (n.eq.0 .and. cur_thresh.le.64*thresh) then 116 cur_thresh = cur_thresh*2 117 goto 111 ! Go back if found nothing to print 118 endif 119 do j = 0, n-1 120 do k = 0, j 121 norm = sqrt(dbl_mb(k_vecsre+int_mb(k_list+k))**2 122 & + dbl_mb(k_vecsim+int_mb(k_list+k))**2) 123 normn = sqrt(dbl_mb(k_vecsre+int_mb(k_list+j))**2 124 & + dbl_mb(k_vecsim+int_mb(k_list+j))**2) 125 if (norm.lt.normn) then 126 m = int_mb(k_list+j) 127 int_mb(k_list+j) = int_mb(k_list+k) 128 int_mb(k_list+k) = m 129 endif 130 enddo 131 enddo 132c 133c Construct optional output line 134c 135 ibuf = 1 136 buf = ' ' 137 if (oocc) then 138 write(buf(ibuf:),'(''Occ='',1p,d12.6)') occ(i) 139 ibuf = ibuf + 18 140 endif 141 if (oevals) then 142 write(buf(ibuf:),'(''E='',1p,d13.6)') evals(i) 143 ibuf = ibuf + 17 144 endif 145 if (oirs) then 146 write(buf(ibuf:),'(''Symmetry='',a4)') zir(irs(i)) 147 ibuf = ibuf + 18 148 endif 149 write(6,1) i, buf(1:max(inp_strlen(buf),1)) 150 1 format(' Vector',i5,2x,a) 151 if (or2) then 152c 153c Construct optional 2nd output line 154c 155 ibuf = 1 156 buf = ' ' 157 write(buf(ibuf:),'(''MO Center='',1p,3(1x,d8.1,'',''))') 158 & dbl_mb(k_r_r2+(i-1)*7), 159 & dbl_mb(k_r_r2+(i-1)*7+1), 160 & dbl_mb(k_r_r2+(i-1)*7+2) 161 ibuf = ibuf + 41 162 write(buf(ibuf:),'(''r^2='',1p,d8.1)') 163 & dbl_mb(k_r_r2+(i-1)*7+6) 164 ibuf = ibuf + 14 165 write(6,3) buf(1:max(inp_strlen(buf),1)) 166 3 format(' ',7x,a) 167 endif 168c 169c Output the analysis 170c 171 write(6,22) 172 22 format(1x,2(' Bfn. Coefficient Function ', 173 $ 4x,4x)) 174 write(6,23) 175 23 format(1x,2(' ---- ------------------- ------------', 176 $ 4x,4x)) 177 do klo = 0, min(n-1,9), 2 178 khi = min(klo+1,n-1) 179 write(6,2) ( 180 $ int_mb(k_list+k)+1, 181 $ dbl_mb(k_vecsre+int_mb(k_list+k)), 182 $ dbl_mb(k_vecsim+int_mb(k_list+k)), 183 $ (byte_mb(k_tags+int_mb(k_list+k)*16+m),m=0,15), 184 $ k = klo,khi) 185 2 format(1x,2(i5,2x,f11.6,f11.6,2x,16a1,4x)) 186 enddo 187 write(6,*) 188 enddo 189 call util_flush(6) 190 if (.not. ma_chop_stack(l_vecsim)) call errquit 191 $ ('bas_vec_info: ma pop?:l_vecsim', 0, MA_ERR) 192 if (.not. ma_chop_stack(l_vecsre)) call errquit 193 $ ('bas_vec_info: ma pop?l_vecsre', 0, MA_ERR) 194 endif 195 if (or2) then 196 if (.not. ma_pop_stack(l_r_r2)) call errquit('movecs_anal?',0, 197 & MA_ERR) 198 endif 199c 200 if (oscfps) call pstat_off(ps_moanal) 201 call ga_sync() 202c 203 end 204 subroutine so_r_and_r2 (basis, g_movecs, r_and_r2) 205c 206 implicit none 207#include "errquit.fh" 208c 209#include "global.fh" 210#include "mafdecls.fh" 211c 212 integer basis ! [input] basis 213 double precision r_and_r2(7,*) ! [output] x,y,z,x^2,y^2,z^2,r^2 for each mo 214 integer g_movecs(2) ! [input] GA mo vectors 215c 216 double precision center(3) 217 integer i, j 218 integer type, nbf, nmo 219 integer g_xlm 220 data center/3*0.0d0/ 221c 222c compute r and r^2 for each MO at the origin 223c 224 call ga_inquire(g_movecs, type, nbf, nmo) 225c 226c create a global array to store x, y, z and x^2, y^2, z^2 for each AO 227c 228 if (.not. ga_create(mt_dbl, 6*nbf, nbf, 'GXLM', 229 $ 32,32,g_xlm)) call errquit('so_r_and_r2 : g_xlm',6*nbf*nbf, 230 & GA_ERR) 231c 232c compute dipoles and quadrupole components for each AO 233c 234 call xlm_so_ao_poles(basis, center, g_xlm) 235c 236c compute dipoles and quadrupole components for each MO 237c 238 call xlm_ao_to_mo_so(g_movecs, g_xlm) 239c 240c put global data into local array on node 0 (node 0 241c will be doing the printing of the desired info) 242c 243 call dfill(nmo*7, 0.0d0, r_and_r2, 1) 244 if (ga_nodeid() .eq. 0) then 245 do i = 1, nmo 246 call ga_get(g_xlm,(i-1)*6+1,(i-1)*6+6,i,i, 247 $ r_and_r2(1,i), 1) 248 enddo 249c 250c write(6,*)' x, y, z, x^2, y^2, z^2 at (0,0,0)' 251c call output(r_and_r2, 1, 7, 1, nmo, 7, nmo, 1) 252c 253c shift and sum 254c 255 do i = 1, nmo 256 do j = 1, 3 257 r_and_r2(j+3,i) = r_and_r2(j+3,i) - r_and_r2(j,i)**2 258 enddo 259 r_and_r2(7,i) = r_and_r2(4,i) + r_and_r2(5,i) + r_and_r2(6,i) 260 enddo 261c write(6,*)' x, y, z, x^2, y^2, z^2, r^2 shifted' 262c call output(r_and_r2, 1, 7, 1, nmo, 7, nmo, 1) 263 endif 264c 265 if (.not. ga_destroy(g_xlm)) call errquit('so_r_and_r2 : ga?',0, 266 & GA_ERR) 267c 268 return 269 end 270 subroutine xlm_so_ao_poles(basis, center, g_xlm) 271c 272 implicit none 273#include "errquit.fh" 274c 275#include "mafdecls.fh" 276#include "global.fh" 277#include "bas.fh" 278#include "geom.fh" 279#include "util_params.fh" 280c 281 integer basis ! [input] basis 282 double precision center(3) ! [input] the expansion center 283 integer g_xlm ! [input] GA that will return the mpoles 284c 285 double precision one, two 286 parameter (one=1.d0, two=2.d0) 287 double precision autoang2 288 parameter (autoang2 = cau2ang*cau2ang) 289c 290 integer mcart, l_xlm, k_xlm 291 integer geom 292 integer nshell, noperators, maxscr, me, nproc 293 integer nbf_max, lmpmax, ishell, jshell, ijshell 294 integer ilo, ihi, jlo, jhi, idim, jdim, ind, i, j, l, ioff 295 integer l_scr, k_scr, l_mp, k_mp 296 integer lmax ! Maximum value of L = 2 297c 298 if (.not. bas_geom(basis, geom)) call errquit 299 $ ('multiplole: bad basis', 0, BASIS_ERR) 300 if (.not. bas_numcont(basis, nshell)) call errquit 301 $ ('xlm_pole: bas_numcont failed for basis', basis, BASIS_ERR) 302 if (.not. bas_nbf_cn_max(basis,nbf_max)) call errquit( 303 & 'xlm_pole: bas_nbf_cn_max failed',20, BASIS_ERR) 304c 305c 306c note lmax is hardwired to 2 307c 308 lmax = 2 309c 310c length of int_mpole integral output for full square list 311c includes l_pole = 0,...,lmax, where l_pole = 0 is simply 312c the 2-c overlap matrix. (cartesian or sphericalcomponents). 313c 314 noperators = (lmax+1)*(lmax+2)*(lmax+3)/6 315 call int_mem_dipole(lmpmax,maxscr,basis,basis,lmax) 316 maxscr = max(100000,maxscr) 317c 318c allocate necessary local temporary arrays on the stack 319c 320 if(.not. ma_push_get(mt_dbl, lmpmax, 'mult:mp', l_mp, k_mp)) 321 & call errquit('xlm_pole: cannot allocate mp', lmpmax, MA_ERR) 322 if(.not. ma_push_get(mt_dbl, lmpmax, 'mult:xlm', l_xlm, k_xlm)) 323 & call errquit('xlm_pole: cannot allocate xlm', lmpmax, MA_ERR) 324 if(.not. ma_push_get(mt_dbl, maxscr, 'mult:scr', l_scr, k_scr)) 325 & call errquit('xlm_pole: cannot allocate scratch', maxscr, 326 & MA_ERR) 327c 328 call ga_zero(g_xlm) 329c 330 ijshell = -1 331 me = ga_nodeid() 332 nproc = ga_nnodes() 333 do ishell = 1, nshell 334 if (.not. bas_cn2bfr(basis, ishell, ilo, ihi)) call errquit 335 & ('xlm_pole: bas_cn2bfr failed for basis', basis, 336 & BASIS_ERR) 337 idim = ihi - ilo + 1 338 339 do jshell = 1, nshell 340 ijshell = ijshell + 1 341 if (mod(ijshell,nproc) .eq. me) then 342 if (.not. bas_cn2bfr(basis, jshell, jlo, jhi)) 343 & call errquit('xlm_pole: bas_cn2bfr', basis, 344 & BASIS_ERR) 345 jdim = jhi - jlo + 1 346c 347 call int_mpole(basis, ishell, basis, jshell, 348 & lmax, center, maxscr, dbl_mb(k_scr), 349 & lmpmax, dbl_mb(k_mp)) 350c 351c output from int_mpole is: overlap, dipole, q-pole, ... 352c within a multipole block, the order is <i|m|j> j fastest, 353c then m, then i ... we must put m first 354c 355 call dfill(6*idim*jdim,0.0d0,dbl_mb(k_xlm),1) 356c 357c write(6,*)' ishell, jshell = ', ishell, jshell 358 ind = k_mp 359 do l = 0, lmax 360 do i = 1, idim 361 do mcart = 1, ((l+1)*(l+2))/2 362 do j = 1, jdim 363 ioff = k_xlm + 6*(j-1 + jdim*(i-1)) 364c write(6,*)' l, i, mcart, j, ind-k_mp, dbl_mb(ind) ', 365c & l, i, mcart, j, ind-k_mp, dbl_mb(ind) 366 if (l.eq.1.and.mcart.eq.1)then 367 dbl_mb(ioff) = dbl_mb(ind)*cau2ang 368c write(6,*)' ioff-k_xlm ', ioff-k_xlm 369 endif 370 if (l.eq.1.and.mcart.eq.2)then 371 ioff = ioff + 1 372 dbl_mb(ioff) = dbl_mb(ind)*cau2ang 373c write(6,*)' ioff-k_xlm ', ioff-k_xlm 374 endif 375 if (l.eq.1.and.mcart.eq.3)then 376 ioff = ioff + 2 377 dbl_mb(ioff) = dbl_mb(ind)*cau2ang 378c write(6,*)' ioff-k_xlm ', ioff-k_xlm 379 endif 380 if (l.eq.2.and.mcart.eq.1)then 381 ioff = ioff + 3 382 dbl_mb(ioff) = dbl_mb(ind)*autoang2 383c write(6,*)' ioff-k_xlm ', ioff-k_xlm 384 endif 385 if (l.eq.2.and.mcart.eq.4)then 386 ioff = ioff + 4 387 dbl_mb(ioff) = dbl_mb(ind)*autoang2 388c write(6,*)' ioff-k_xlm ', ioff-k_xlm 389 endif 390 if (l.eq.2.and.mcart.eq.6)then 391 ioff = ioff + 5 392 dbl_mb(ioff) = dbl_mb(ind)*autoang2 393c write(6,*)' ioff-k_xlm ', ioff-k_xlm 394 endif 395 ind = ind + 1 396 end do 397 end do 398 end do 399 end do 400c 401 call ga_put(g_xlm,1+(jlo-1)*6,jhi*6,ilo,ihi, 402 $ dbl_mb(k_xlm), 6*jdim) 403c 404 end if 405 end do 406 end do 407c 408 call ga_sync 409c 410c write(6,*) ' THE AO MPOLES ' 411c call ga_print(g_xlm) 412c 413c clean up stack 414c 415 if (.not. ma_pop_stack(l_scr)) call errquit('xlm_pole: ma?',0, 416 & MA_ERR) 417 if (.not. ma_pop_stack(l_xlm)) call errquit('xlm_pole: ma?',0, 418 & MA_ERR) 419 if (.not. ma_pop_stack(l_mp)) call errquit('xlm_pole: ma?',0, 420 & MA_ERR) 421c 422 end 423 subroutine xlm_ao_to_mo_so(g_vecs, g_xlm) 424 implicit none 425#include "errquit.fh" 426#include "global.fh" 427#include "mafdecls.fh" 428 integer lmax, g_vecs(2), g_xlm 429c 430c Transform the multipoles from AO to MO basis overwriting 431c the input AO set. 432c 433 integer type, nbf, nmo, g_tmp1, g_tmp2, k_tmp, l_tmp, 434 $ i, lm 435c 436c note lmax is hardwired to 2 437c 438 lmax = 2 439c 440 call ga_inquire(g_vecs, type, nbf, nmo) 441 if (.not. ga_create(mt_dbl, nbf, nbf, 'aomotmp', 442 $ 32,32,g_tmp1)) call errquit('xlm_ao_mo: tmp1',nbf*nbf, 443 & GA_ERR) 444 if (.not. ga_create(mt_dbl, nbf, nmo, 'aomotmp2', 445 $ 32,32,g_tmp2)) call errquit('xlm_ao_mo: tmp2',nmo*nbf, 446 & GA_ERR) 447 if (.not. ma_push_get(mt_dbl,nbf,'xlmtpm',l_tmp, k_tmp)) 448 $ call errquit('xlm_ao_mo: tmp', nbf, MA_ERR) 449c 450c Must transform the LHS index one mpole at a time so 451c might as well do both at the same time since this will 452c use less memory. 453c 454 do lm = 1, 6 455 call ga_sync 456 do i = 1+ga_nodeid(), nbf, ga_nnodes() 457 call ga_get(g_xlm,lm+(i-1)*6,lm+(i-1)*6,1,nbf, 458 $ dbl_mb(k_tmp), 1) 459 call ga_put(g_tmp1, i, i, 1, nbf, dbl_mb(k_tmp), 1) 460 end do 461 call ga_sync() 462 call ga_dgemm('n','n',nbf,nmo,nbf,1.0d0,g_tmp1,g_vecs, 463 $ 0.0d0,g_tmp2) 464 call ga_dgemm('t','n',nmo,nmo,nbf,1.0d0,g_vecs,g_tmp2, 465 $ 0.0d0,g_tmp1) 466 call ga_sync 467 do i = 1+ga_nodeid(), nmo, ga_nnodes() 468 call ga_get(g_tmp1, i, i, 1, nmo, dbl_mb(k_tmp), 1) 469 call ga_put(g_xlm,lm+(i-1)*6,lm+(i-1)*6,1,nmo, 470 $ dbl_mb(k_tmp), 1) 471 end do 472 call ga_sync() 473 end do 474c 475c write(6,*) ' THE MO MPOLES' 476c call ga_print(g_xlm) 477c 478 if (.not. ga_destroy(g_tmp1)) call errquit('xlm_ao_mo?',1, GA_ERR) 479 if (.not. ga_destroy(g_tmp2)) call errquit('xlm_ao_mo?',2, GA_ERR) 480 if (.not. ma_pop_stack(l_tmp)) call errquit('xlm_ao_mo?',3, 481 & MA_ERR) 482c 483 end 484