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