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