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