1************************************************************************
2      subroutine gucc_fock(irefspc,itrefspc,
3     &                     ccvec1,ccvec2,ccvec3,ccvec4,
4     &                     civec1,civec2,c2vec,
5     &                     n_cc_typ,i_cc_typ,
6     &                     namp_cc_typ,ioff_cc_typ,
7     &                     iopsym,mxb_ci,
8     &                     luc,luamp,luomg,ludia,
9     &                     luec,luhc)
10************************************************************************
11*
12* purpose : driver for tests of fock-space GUCC
13*
14*  ak, early 2004
15*
16************************************************************************
17*
18* units:
19*   luc   = definition of reference function
20*   luamp = amplitude vectors (also output for most recent vector)
21*   luampold = scratch containing old vectors from previous iterations
22*           (on input it may also be a first trial vector)
23*   luomg = error vectors
24*   ludia = diagonal preconditioner
25*   luec  = scratch for exp(G)|ref>
26*   luhc  = scratch for H exp(G)|ref>
27
28* diverse inludes with commons and paramters
29c      include 'implicit.inc'
30c      include 'mxpdim.inc'
31      include 'wrkspc.inc'
32      include 'crun.inc'
33      include 'cstate.inc'
34      include 'cgas.inc'
35      include 'ctcc.inc'
36      include 'gasstr.inc'
37      include 'strinp.inc'
38      include 'orbinp.inc'
39      include 'lucinp.inc'
40      include 'cprnt.inc'
41      include 'corbex.inc'
42      include 'csm.inc'
43      include 'cecore.inc'
44      include 'gtbce.inc'
45      include 'cintfo.inc'
46      include 'glbbas.inc'
47      include 'opti.inc'
48* constants
49      integer, parameter ::
50     &     ntest = 5
51
52* arrays
53      integer ::
54     &     ioff_cc_typ(n_cc_typ), namp_cc_typ(n_cc_typ),
55     &     i_cc_typ(4*ngas,n_cc_typ)
56      real*8 ::
57     &     ccvec1(n_cc_amp), ccvec2(n_cc_amp), ccvec3(n_cc_amp)
58* local
59      logical ::
60     &     calc_Omg, calc_gradE, tstgrad, tst_hss, comm_ops,
61     &     do_eag, do_foo, do_hss
62      character*8 cctype
63      integer ::
64     &     ictp(n_cc_typ)
65* external functions
66      real*8, external :: inprod, inprdd
67
68* ============================================================
69* initialize : restart, set coefs to zero
70* ============================================================
71
72      call atim(cpu0,wall0)
73
74      nprint = 1
75      lblk = -1
76
77      if (ntest.ge.5) then
78        write(6,*) '======================='
79        write(6,*) 'entered gucc_fock with:'
80        write(6,*) '======================='
81        write(6,*) ' iopsym = ',iopsym
82      end if
83
84* unit inits
85      lusc1 = iopen_nus('GTBSCR1')
86      lusc2 = iopen_nus('GTBSCR2')
87      lusc3 = iopen_nus('GTBSCR3')
88      lusc4 = iopen_nus('GTBSCR4')
89      lusc5 = iopen_nus('GTBSCR5')
90c      lusc6 = iopen_nus('GTBSCR6')
91c      lusc7 = iopen_nus('GTBSCR7')
92c      lusc8 = iopen_nus('GTBSCR8')
93c      lusc9 = iopen_nus('GTBSCR9')
94
95      call memman(idum,idum,'MARK  ',idum,'FOCKSP')
96
97! be sure to have the correct h1 integrals:
98      call copvec(work(kint1o),work(kint1),nint1)
99
100      npairs = 4*ntoob**2
101      nheff = npairs*(npairs+1)/2
102      call memman(kheff,nheff,'ADDL  ',2,'H_EFF ')
103      call memman(kheff2,npairs*npairs,'ADDL  ',2,'H_EFF ')
104      nwmat = npairs*npairs
105      call memman(kwmat,nwmat,'ADDL  ',2,'W_MAT ')
106      call memman(keig,npairs,'ADDL  ',2,'EIGEN ')
107
108      npairs = (2*ntoob-1)*(2*ntoob)/2
109      nel = nactel
110
111      print *,' NEL = ', nel
112
113      call diag_effH(work(kheff),work(kheff2),
114     &     work(kwmat),work(keig),nel)
115
116      call memchk2('a effH')
117
118c     call logwmat(work(kwmat))
119
120      call cleanvec(work(kwmat),npairs**2)
121
122      call rearr_g(work(kwmat),ccvec1,
123     &     n_cc_typ,i_cc_typ,ioff_cc_typ)
124
125      igtbmod = 0
126      iopsym = 0
127
128      ! important:
129      call scalve(ccvec1,-1d0,n_cc_amp)
130      call gtbce_E(igtbmod,elen1,variance,ovl,
131     &               ecore,
132     &               ccvec1,iopsym,ccvec4,
133     &               civec1,civec2,c2vec,
134     &               n_cc_amp,mxb_ci,
135     &               luc,luec,luhc,lusc1,lusc2)
136      h1 = elen1+ecore
137
138      ! effective Hamiltonian:
139      call rearr_g(work(kheff2),ccvec3,
140     &     n_cc_typ,i_cc_typ,ioff_cc_typ)
141
142      call memchk2('a r g')
143
144      ! recalc energy with effective hamilt.
145      isigden = 1
146      call sigden_cc(civec1,civec2,luc,lusc1,ccvec3,isigden)
147      hf1 = ecore - inprdd(civec1,civec2,luc,lusc1,1,-1)
148      print *,'TEST: E(HF) = ',hf1
149
150
151      ! create a 2-fold excited determinant
152      ccvec2(ioff_cc_typ(25)+2) = 1d0
153      isigden = 1
154      call sigden_cc(civec1,civec2,luc,lusc1,ccvec2,isigden)
155
156      call gtbce_E(igtbmod,elen2,variance,ovl,
157     &               ecore,
158     &               ccvec1,iopsym,ccvec4,
159     &               civec1,civec2,c2vec,
160     &               n_cc_amp,mxb_ci,
161     &               lusc1,lusc2,lusc3,lusc4,lusc5)
162                    !ref   eGref HeGref
163
164      h2 = elen2+ecore
165      ! <ex|e(-G)He(G)|0>:
166      h12 = inprdd(civec1,civec2,lusc2,luhc,1,-1)
167      h21 = inprdd(civec1,civec2,lusc3,luec,1,-1)
168
169      print *,' :: ',h1 , h21
170      print *,' :: ',h12, h2
171
172      ! Heff|expG 1> on lusc4
173      isigden = 1
174      call sigden_cc(civec1,civec2,luec,lusc4,ccvec3,isigden)
175      h1t = ecore - inprdd(civec1,civec2,luec,lusc4,1,-1)
176      ! Heff|expG 2> on lusc5
177      isigden = 1
178      call sigden_cc(civec1,civec2,lusc2,lusc5,ccvec3,isigden)
179      h2t = ecore - inprdd(civec1,civec2,lusc2,lusc5,1,-1)
180
181      h12t = -inprdd(civec1,civec2,lusc2,lusc4,1,-1)
182      h21t = -inprdd(civec1,civec2,luec,lusc5,1,-1)
183
184      print *,'TEST:'
185      print *,' :: ',h1t , h21t
186      print *,' :: ',h12t, h2t
187
188
189      call relunit(lusc1,'delete')
190      call relunit(lusc2,'delete')
191      call relunit(lusc3,'delete')
192      call relunit(lusc4,'delete')
193      call relunit(lusc5,'delete')
194c      call relunit(lusc6,'delete')
195c      call relunit(lusc7,'delete')
196c      call relunit(lusc8,'delete')
197c      call relunit(lusc9,'delete')
198
199      call memman(idum,idum,'FLUSM ',idum,'FOCKSP')
200
201      return
202
203      end
204
205************************************************************************
206
207      subroutine diag_effH(h_eff,h_eff_out,wmat,eig,nel)
208************************************************************************
209*
210*     set up
211*           H(pq,rs) = (pq|rs) + 1/(N-1)(h(pq)dlt(rs)+h(rs)dlt(pq))
212*
213*     and diagonalize in the sense:
214*
215*       W(tu,pq)H(pq,rs)W(rs,vw) = e(tu)dlt(tu,vw)
216*
217************************************************************************
218
219      include 'wrkspc.inc'
220      include 'crun.inc'
221      include 'cstate.inc'
222      include 'cgas.inc'
223      include 'ctcc.inc'
224      include 'gasstr.inc'
225      include 'strinp.inc'
226      include 'orbinp.inc'
227      include 'cprnt.inc'
228      include 'corbex.inc'
229      include 'csm.inc'
230      include 'cecore.inc'
231      include 'gtbce.inc'
232      include 'opti.inc'
233      include 'cintfo.inc'
234      include 'glbbas.inc'
235      include 'oper.inc'
236
237      dimension h_eff(*), wmat(*), eig(*), h_eff_out(*)
238
239      ! set up h_eff
240
241      npairs = 2*ntoob*(2*ntoob-1)/2
242      h_eff(1:npairs*(npairs+1)/2) = 0d0
243      xel = dble(nel)
244
245      ipr = 0
246
247      i_unrorb = 0
248      ispcas = 0
249      i_use_simtrh = 0
250
251      do mp = 1, 2
252c old:
253c      do ip = 1, ntoob
254c        do mr = 1, mp
255c
256c new:
257      do mr = 1, mp
258        do ip = 1, ntoob
259c
260          irmx = ntoob
261          if (mr.eq.mp) irmx = ip-1
262        do ir = 1, irmx
263          ipr = ipr + 1
264!          ips = (imp-1)*2*ntoob + ims
265
266          iqs = 0
267
268          do mq = 1, 2
269c old:
270c          do iq = 1, ntoob
271c            do ms = 1, mq
272c
273c new:
274          do ms = 1, mq
275            do iq = 1, ntoob
276c
277              ismx = ntoob
278              if (ms.eq.mq) ismx = iq-1
279            do is = 1, ismx
280              iqs = iqs + 1
281!              iqr = (imq-1)*2*ntoob + imr
282
283              if (ipr.lt.iqs) cycle ! hermitian operator
284
285              idxpqrs = ipr*(ipr-1)/2 + iqs
286
287c              write(6,*) 'ip,iq,ir,is: ',ip,iq,ir,is
288c              write(6,*) 'mp,mq,mr,ms: ',mp,mq,mr,ms
289cc
290c              write(6,*) 'ipr,iqs,idxpqrs: ', ipr,iqs,idxpqrs
291cc
292              if (mp.eq.mq.and.mr.eq.ms) then
293c
294c                write(6,*) ' before    : ', h_eff(idxpqrs)
295c
296                h_eff(idxpqrs) = h_eff(idxpqrs) + gtijkl(ip,iq,ir,is)
297c
298c                write(6,*) ' +2-el int : ', h_eff(idxpqrs)
299c
300                if (ip.eq.iq)
301     &               h_eff(idxpqrs) = h_eff(idxpqrs) +
302     &                                  1d0/(xel-1d0)*geth1i_2(ir,is)
303
304                if (ir.eq.is)
305     &               h_eff(idxpqrs) = h_eff(idxpqrs) +
306     &                                  1d0/(xel-1d0)*geth1i_2(ip,iq)
307
308c                write(6,*) ' +1-el int : ', h_eff(idxpqrs)
309              end if
310              if (mp.eq.ms.and.mq.eq.mr) then
311c
312c                write(6,*) ' before    : ', h_eff(idxpqrs)
313c
314                h_eff(idxpqrs) = h_eff(idxpqrs) - gtijkl(ip,is,ir,iq)
315
316c                write(6,*) ' +2-el int : ', h_eff(idxpqrs)
317
318                if (ip.eq.is)
319     &               h_eff(idxpqrs) = h_eff(idxpqrs) -
320     &                                  1d0/(xel-1d0)*geth1i_2(ir,iq)
321
322                if (ir.eq.iq)
323     &               h_eff(idxpqrs) = h_eff(idxpqrs) -
324     &                                  1d0/(xel-1d0)*geth1i_2(ip,is)
325
326c                write(6,*) ' +1-el int : ', h_eff(idxpqrs)
327
328              end if
329
330              idx2pqrs = (ipr-1)*npairs+iqs
331              idx2qpsr = (iqs-1)*npairs+ipr
332
333              h_eff_out(idx2pqrs) = h_eff(idxpqrs)
334              h_eff_out(idx2qpsr) = h_eff(idxpqrs)
335
336            end do
337            end do
338          end do
339          end do
340        end do
341        end do
342      end do
343      end do
344
345      write(6,*) 'the hamiltonian:'
346c      write(6,*) ' norm^2 = ',ddot(npairs**2,h_eff,1,h_eff,1)
347
348      call prtrlt(h_eff,npairs)
349
350      call copdia(h_eff,eig,npairs,1)
351
352      write(6,*) 'the diagonal of the hamiltonian:'
353
354      call wrtmat(eig,npairs,1,npairs,1)
355
356      ! call diagonalizer
357
358c      call eigen(h_eff,wmat,npairs,0,1)
359      call eigen(h_eff,wmat,npairs,0,0)
360
361      call copdia(h_eff,eig,npairs,1)
362
363      write(6,*) 'the eigenvalues of the hamiltonian:'
364
365      call wrtmat(eig,npairs,1,npairs,1)
366
367      call memman(keigr,npairs,'ADDL  ',2,'EIG RE')
368      call memman(keigi,npairs,'ADDL  ',2,'EIG IM')
369
370      call dumsort(eig,npairs,npairs,work(keigr),work(keigi))
371      call reo_vec(work(keigr),eig,npairs,work(keigi),1)
372      call copvec(work(keigi),eig,npairs)
373
374      write(6,*) 'the eigenvalues of the hamiltonian (sorted):'
375
376      call wrtmat(eig,npairs,1,npairs,1)
377
378      write(6,*) 'E_nuc = ',ecore
379
380      ener = ecore
381      do ii = 1, nel*(nel-1)/2
382        ener = ener + eig(ii)
383        print *,ii,eig(ii),ener
384      end do
385
386      write(6,*) 'the W-matrix: ',sqrt(ddot(npairs**2,wmat,1,wmat,1))
387
388      call wrtmat2(wmat,npairs,npairs,npairs,npairs)
389
390c      return
391
392c      call memman(keigr,npairs,'ADDL  ',2,'EIG RE')
393c      call memman(keigi,npairs,'ADDL  ',2,'EIG IM')
394      call memman(klneig,npairs,'ADDL  ',2,'EIG IM')
395      call memman(keigv,npairs**2,'ADDL  ',2,'EIGVEC')
396      call memman(keigvr,npairs**2,'ADDL  ',2,'EIGVEC')
397      call memman(keigvi,npairs**2,'ADDL  ',2,'EIGVEC')
398      call memman(keigvr2,npairs**2,'ADDL  ',2,'EIGVEC')
399      call memman(keigvi2,npairs**2,'ADDL  ',2,'EIGVEC')
400      call memman(kgmatr,npairs**2,'ADDL  ',2,'GMAT R')
401      call memman(kgmati,npairs**2,'ADDL  ',2,'GMAT I')
402      call memman(kiscr,npairs,'ADDL  ',1,'INTSCR')
403      call memman(kscr,npairs,'ADDL  ',2,'RELSCR')
404
405      call logumat(npairs,work(kgmatr),wmat,
406     &     work(keigv),work(keigvr),work(keigvi))
407c
408c      call copvec(wmat,work(keigvr),npairs**2)
409c
410c      call rg(npairs,npairs,work(keigvr),
411c     &     work(keigr),work(keigi),1,work(keigv),
412c     &     work(kiscr),work(kscr),ierr)
413c      call nrmvec(npairs,work(keigv),work(keigi))
414c
415c      print *,'eigenvalues of W-matrix:'
416c
417c      do ii = 1, npairs
418c        eigr = work(keigr-1+ii)
419c        eigi = work(keigi-1+ii)
420c        rad  = eigr*eigr+eigi*eigi
421c        costau = eigr
422c        tau = acos(eigr)*sign(1d0,eigi)
423c        work(klneig-1+ii) = tau
424c        print '(i4,2(2x,e20.10),2(2x,f15.10))',ii,eigr,eigi,rad,
425c     &       tau/3.1415926535897932d0*180d0
426c      end do
427c
428cc      print *,'the eigenvectors'
429cc      call wrtmat2(work(keigv),npairs,npairs,npairs,npairs)
430c
431c      ! sort components into real and imaginary matrices
432c
433c      ii = 0
434c      do while(ii.lt.npairs)
435c        ii = ii+1
436c        eigr = work(keigr-1+ii)
437c        eigi = work(keigi-1+ii)
438c        if (eigi.eq.0d0) then
439c          call copvec(work(keigv+(ii-1)*npairs),
440c     &                work(keigvr+(ii-1)*npairs),npairs)
441c          call setvec(work(keigvi+(ii-1)*npairs),0d0,npairs)
442c        else
443c          call copvec(work(keigv+(ii-1)*npairs),
444c     &                work(keigvr+(ii-1)*npairs),npairs)
445c          call copvec(work(keigv+(ii-1)*npairs),
446c     &                work(keigvr+(ii)*npairs),npairs)
447c          call copvec(work(keigv+(ii)*npairs),
448c     &                work(keigvi+(ii-1)*npairs),npairs)
449c          call copvec(work(keigv+(ii)*npairs),
450c     &                work(keigvi+(ii)*npairs),npairs)
451c          call scalve(work(keigvi+(ii)*npairs),-1d0,npairs)
452c          ii = ii+1  ! additional increment
453c        end if
454c
455c      end do
456c
457cc      print *,'the eigenvectors (re):'
458cc      call wrtmat2(work(keigvr),npairs,npairs,npairs,npairs)
459cc      print *,'the eigenvectors (im):'
460cc      call wrtmat2(work(keigvi),npairs,npairs,npairs,npairs)
461c
462cc      ! test the trafo matrices:
463cc
464cc      ! A^T A - B^T B:
465cc      call dgemm('t','n',npairs,npairs,npairs,
466cc     &     1d0,work(keigvr),npairs,
467cc     &         work(keigvr),npairs,
468cc     &     0d0,work(keigvr2),npairs)
469cc      call dgemm('t','n',npairs,npairs,npairs,
470cc     &     1d0,work(keigvi),npairs,
471cc     &         work(keigvi),npairs,
472cc     &     1d0,work(keigvr2),npairs)
473cc
474cc      ! A^T B + B^T A:
475cc      call dgemm('t','n',npairs,npairs,npairs,
476cc     &     1d0,work(keigvr),npairs,
477cc     &         work(keigvi),npairs,
478cc     &     0d0,work(keigvi2),npairs)
479cc      call dgemm('t','n',npairs,npairs,npairs,
480cc     &    -1d0,work(keigvi),npairs,
481cc     &         work(keigvr),npairs,
482cc     &     1d0,work(keigvi2),npairs)
483cc
484cc      print *,'U^T U, real part:'
485cc      call wrtmat2(work(keigvr2),npairs,npairs,npairs,npairs)
486cc
487cc      print *,'U^T U, imaginary part:'
488cc      call wrtmat2(work(keigvi2),npairs,npairs,npairs,npairs)
489cc
490cc
491cc      ! W A
492cc      call dgemm('n','n',npairs,npairs,npairs,
493cc     &           1d0,wmat,npairs,
494cc     &               work(keigvr),npairs,
495cc     &           0d0,work(keigvr2),npairs)
496cc
497cc      ! W B
498cc      call dgemm('n','n',npairs,npairs,npairs,
499cc     &           1d0,wmat,npairs,
500cc     &               work(keigvi),npairs,
501cc     &           0d0,work(keigvi2),npairs)
502cc
503cc      ! Re D = A^T W A - B^T W B
504cc      call dgemm('t','n',npairs,npairs,npairs,
505cc     &           1d0,work(keigvi),npairs,
506cc     &               work(keigvi2),npairs,
507cc     &           0d0,work(kgmatr),npairs)
508cc      call dgemm('t','n',npairs,npairs,npairs,
509cc     &           1d0,work(keigvr),npairs,
510cc     &               work(keigvr2),npairs,
511cc     &           1d0,work(kgmatr),npairs)
512cc
513cc      ! Im D = A^T W B + B^T W A
514cc      call dgemm('t','n',npairs,npairs,npairs,
515cc     &           1d0,work(keigvr),npairs,
516cc     &               work(keigvi2),npairs,
517cc     &           0d0,work(kgmati),npairs)
518cc      call dgemm('t','n',npairs,npairs,npairs,
519cc     &          -1d0,work(keigvi),npairs,
520cc     &               work(keigvr2),npairs,
521cc     &           1d0,work(kgmati),npairs)
522cc
523cc      print *,'U^T W U, real part:'
524cc      xnrm = sqrt(ddot(npairs**2,work(kgmatr),1,work(kgmatr),1 ))
525cc      print *,'norm = ',xnrm
526cc
527cc      call wrtmat2(work(kgmatr),npairs,npairs,npairs,npairs)
528cc
529cc      print *,'U^T W U, imaginary part:'
530cc      call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs)
531c
532c      ! D A^+
533c      do ii = 1, npairs
534c        tau = work(klneig-1+ii)
535c        do jj = 1, npairs
536c          ijadr = (jj-1)*npairs+ii-1  ! uahhh
537c          jiadr = (ii-1)*npairs+jj-1  ! uahhh
538c          work(keigvr2+ijadr) = tau * work(keigvr+jiadr)
539c        end do
540c      end do
541c
542c      ! D B^+
543c      do ii = 1, npairs
544c        tau = work(klneig-1+ii)
545c        do jj = 1, npairs
546c          ijadr = (jj-1)*npairs+ii-1
547c          jiadr = (ii-1)*npairs+jj-1
548c          work(keigvi2+ijadr) = tau * work(keigvi+jiadr)
549c        end do
550c      end do
551c
552c      ! Re G = - B D A^+ + A D B^+
553c      call dgemm('n','n',npairs,npairs,npairs,
554c     &           -1d0,work(keigvi),npairs,
555c     &                work(keigvr2),npairs,
556c     &            0d0,work(kgmatr),npairs)
557c      call dgemm('n','n',npairs,npairs,npairs,
558c     &           +1d0,work(keigvr),npairs,
559c     &                work(keigvi2),npairs,
560c     &            1d0,work(kgmatr),npairs)
561c
562c      ! Im G = A D A^+ + B D B^+
563c      call dgemm('n','n',npairs,npairs,npairs,
564c     &            1d0,work(keigvr),npairs,
565c     &                work(keigvr2),npairs,
566c     &            0d0,work(kgmati),npairs)
567c      call dgemm('n','n',npairs,npairs,npairs,
568c     &            1d0,work(keigvi),npairs,
569c     &                work(keigvi2),npairs,
570c     &            1d0,work(kgmati),npairs)
571c
572c      print *,'G, real part:'
573c
574c      call wrtmat2(work(kgmatr),npairs,npairs,npairs,npairs)
575c
576c      print *,'G, imaginary part:'
577c      xnrm = sqrt(ddot(npairs**2,work(kgmatr),1,work(kgmatr),1 ))
578c      print *,'norm = ',xnrm
579c
580c      if (xnrm.gt.1d-10)
581c     &     call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs)
582
583      call expmat(work(kgmati),work(kgmatr),work(keigv),work(keigvr),
584     &     npairs,1d-20)
585
586      print *,'W from exp(G):'
587      call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs)
588**
589
590      call copvec(work(kgmatr),wmat,npairs**2)
591
592c      stop 'end of test'
593
594      end
595
596      subroutine expmat(expx,xmat,xscr1,xscr2,ndim,thrsh)
597
598*     calculate exp(X), returned on expx, of (ndim,ndim)-matrix X,
599*     input on xmat, by Taylor expansion (threshold thrsh)
600*     xscr is a scratch matrix of the same dimensions as xmat, expx
601
602      implicit none
603
604      integer, parameter :: ntest = 100, maxn = 100
605
606      integer, intent(in) ::
607     &     ndim
608      real(8), intent(in) ::
609     &     thrsh
610      real(8), intent(inout) ::
611     &     expx(ndim,ndim), xmat(ndim,ndim),
612     &     xscr1(ndim,ndim),  xscr2(ndim,ndim)
613
614      logical ::
615     &     conv
616      integer ::
617     &     n, ndim2, ii
618      real(8) ::
619     &     xnrm, fac
620
621      real(8), external ::
622     &     inprod
623
624      expx(1:ndim,1:ndim) = xmat(1:ndim,1:ndim)
625      xscr2(1:ndim,1:ndim) = xmat(1:ndim,1:ndim)
626
627      do ii = 1, ndim
628        expx(ii,ii) = expx(ii,ii) + 1d0
629      end do
630
631      ndim2 = ndim*ndim
632      n = 1
633      conv = .false.
634
635      do while (.not.conv)
636        n = n+1
637        if (n.gt.maxn) exit
638
639        fac = 1d0/dble(n)
640
641        ! Xscr = 1/N Xscr * X
642        call matml7(xscr1,xscr2,xmat,
643     &              ndim,ndim,
644     &              ndim,ndim,
645     &              ndim,ndim,0d0,fac,0)
646
647        xnrm = sqrt(inprod(xscr1,xscr1,ndim2))
648        if (xnrm.lt.thrsh) conv = .true.
649
650        if (ntest.ge.10)
651     &       write(6,*) ' N = ',n,'  |1/N! X^N| = ',xnrm
652
653        call vecsum(expx,expx,xscr1,1d0,1d0,ndim2)
654
655        call copvec(xscr1,xscr2,ndim2)
656
657      end do
658
659      if (.not.conv) then
660        write(6,*) ' Taylor expansion of exp(X) did not converge!'
661        stop 'expmat'
662      end if
663
664      return
665      end
666c
667      subroutine rearr_g(gmat_in,gmat_out,ntss_tp,itss_tp,ibtss_tp)
668
669      include 'implicit.inc'
670      include 'mxpdim.inc'
671      include 'cgas.inc'
672      include 'multd2h.inc'
673      include 'orbinp.inc'
674      include 'csm.inc'
675      include 'ctcc.inc'
676      include 'cc_exc.inc'
677
678      integer, parameter ::
679     &     ntest = 1000
680
681* input
682      real(8), intent(inout) ::
683     &     gmat_in(2*ntoob*(2*ntoob-1)/2,2*ntoob*(2*ntoob-1)/2),
684     &     gmat_out(*)
685
686c      input needed: itss_tp <-- work(klsobex), ntss_tp <-- nspobex_tp
687      integer, intent(in) ::
688     &     ntss_tp,
689     &     itss_tp(ngas,4,ntss_tp),
690     &     ibtss_tp(ntss_tp)
691
692* local
693      integer ::
694     &     igrp_ca(mxpngas), igrp_cb(mxpngas),
695     &     igrp_aa(mxpngas), igrp_ab(mxpngas),
696     &     iocc_ca(mx_st_tsoso_blk_mx),
697     &     iocc_cb(mx_st_tsoso_blk_mx),
698     &     iocc_aa(mx_st_tsoso_blk_mx),
699     &     iocc_ab(mx_st_tsoso_blk_mx),
700     &     idx_c(4), idx_s(4)
701
702
703      npairs = 2*ntoob*(2*ntoob-1)/2
704
705      if (ntest.ge.1000) then
706        write(6,*) ' input amplitudes: '
707        call wrtmat2(gmat_in,npairs,npairs,npairs,npairs)
708        write(6,*) 'ibtss_tp:'
709        call iwrtma(ibtss_tp,1,ntss_tp,1,ntss_tp)
710      end if
711
712
713        ! loop over types
714        do itss = 1, ntss_tp
715          idx = ibtss_tp(itss) - 1
716
717          ! identify two-particle excitations:
718          nel_ca = ielsum(itss_tp(1,1,itss),ngas)
719          nel_cb = ielsum(itss_tp(1,2,itss),ngas)
720          nel_aa = ielsum(itss_tp(1,3,itss),ngas)
721          nel_ab = ielsum(itss_tp(1,4,itss),ngas)
722          nc = nel_ca + nel_cb
723          na = nel_aa + nel_ab
724          if (na.ne.2) cycle
725          ! transform occupations to groups
726          call occ_to_grp(itss_tp(1,1,itss),igrp_ca,1)
727          call occ_to_grp(itss_tp(1,2,itss),igrp_cb,1)
728          call occ_to_grp(itss_tp(1,3,itss),igrp_aa,1)
729          call occ_to_grp(itss_tp(1,4,itss),igrp_ab,1)
730
731          if (mscomb_cc.ne.0) then
732            call diag_exc_cc(itss_tp(1,1,itss),itss_tp(1,2,itss),
733     &           itss_tp(1,3,itss),itss_tp(1,4,itss),
734     &           ngas,idiag)
735          else
736            idiag = 0
737          end if
738
739          ! loop over symmetry blocks
740          ism = 1 ! totally symmetric operators
741          do ism_c = 1, nsmst
742            ism_a = multd2h(ism,ism_c)
743            do ism_ca = 1, nsmst
744            ism_cb = multd2h(ism_c,ism_ca)
745            do ism_aa = 1, nsmst
746              ism_ab = multd2h(ism_a,ism_aa)
747              ! get alpha and beta symmetry index
748              ism_alp = (ism_aa-1)*nsmst+ism_ca  ! = (sym Ca,sym Aa)
749              ism_bet = (ism_ab-1)*nsmst+ism_cb  ! = (sym Cb,sym Ab)
750
751              ! restrict to (sym Ca,sym Aa) >= (sym Cb,sym Ab)
752              if (idiag.eq.1.and.ism_bet.gt.ism_alp) cycle
753              if (idiag.eq.0.or.ism_alp.gt.ism_bet) then
754                irestr = 0
755              else
756                irestr = 1
757              end if
758
759              ! get the strings
760              call getstr2_totsm_spgp(igrp_ca,ngas,ism_ca,nel_ca,
761     &             lca,iocc_ca,norb,0,idum,idum)
762              call getstr2_totsm_spgp(igrp_cb,ngas,ism_cb,nel_cb,
763     &             lcb,iocc_cb,norb,0,idum,idum)
764              call getstr2_totsm_spgp(igrp_aa,ngas,ism_aa,nel_aa,
765     &             laa,iocc_aa,norb,0,idum,idum)
766              call getstr2_totsm_spgp(igrp_ab,ngas,ism_ab,nel_ab,
767     &             lab,iocc_ab,norb,0,idum,idum)
768
769              ! length of strings in this symmetry block
770              if (lca*lcb*laa*lab.eq.0) cycle
771
772              do iab = 1, lab
773                if (irestr.eq.1) then
774                  iaa_min = iab
775                else
776                  iaa_min = 1
777                end if
778                do iaa = iaa_min, laa
779                  do icb = 1, lcb
780                    if (irestr.eq.1.and.iaa.eq.iab) then
781                      ica_min = icb
782                    else
783                      ica_min = 1
784                    end if
785                    do ica = ica_min, lca
786                      idx = idx + 1
787                      ! translate into canonical index quadrupel
788                      ii = 0
789                      do iel = 1, nel_ca
790                        ii = ii + 1
791                        idx_c(ii) = iocc_ca((ica-1)*nel_ca+iel)
792                        idx_s(ii) = 1
793                      end do
794                      do iel = 1, nel_cb
795                        ii = ii + 1
796                        idx_c(ii) = iocc_cb((icb-1)*nel_cb+iel)
797                        idx_s(ii) = 2
798                      end do
799                      do iel = 1, nel_aa
800                        ii = ii + 1
801                        idx_c(ii) = iocc_aa((iaa-1)*nel_aa+iel)
802                        idx_s(ii) = 1
803                        idx_s(ii) = 1
804                      end do
805                      do iel = 1, nel_ab
806                        ii = ii + 1
807                        idx_c(ii) = iocc_ab((iab-1)*nel_ab+iel)
808                        idx_s(ii) = 2
809                      end do
810
811                      ! idx_c contains:
812                      ! p1, p2, h1, h2
813                      ! idx_s contains
814                      ! mp1, mp2, mh1, mh2
815
816                      print *,'p, r: ',idx_c(1),idx_c(2)
817                      print *,'      ',idx_s(1),idx_s(2)
818                      print *,'q, s: ',idx_c(3),idx_c(4)
819                      print *,'      ',idx_s(3),idx_s(4)
820
821                      idxpr = igtidx(idx_c(1),idx_c(2),
822     &                              idx_s(1),idx_s(2),ntoob)
823                      idxqs = igtidx(idx_c(3),idx_c(4),
824     &                              idx_s(3),idx_s(4),ntoob)
825
826                      print *,' >>> ',idx,' ---> ',idxpr,idxqs
827
828                      gmat_out(idx) = gmat_in(idxpr,idxqs)
829
830c testing coverage:
831                      if (gmat_in(idxpr,idxqs).ne.1d100) then
832                        gmat_in(idxpr,idxqs) = 1d100
833                      else
834                        gmat_in(idxpr,idxqs) =
835     &                       gmat_in(idxpr,idxqs)+1d100
836                      end if
837
838                    end do ! ica
839                  end do ! icb
840                end do ! iaa
841              end do ! iab
842
843            end do ! ism_aa
844            end do ! ism_ca
845          end do ! ism_c
846
847        end do ! itss
848
849      if (ntest.ge.1000) then
850        write(6,*) 'coverage:'
851        call wrtmat2(gmat_in,npairs,npairs,npairs,npairs)
852        write(6,*) ' output amplitudes: '
853        call wrt_cc_vec2(gmat_out,6,'GEN_CC')
854      end if
855
856      return
857      end
858
859      integer function igtidx(ip,ir,mp,mr,ndim)
860
861      implicit none
862
863      integer, intent(in) ::
864     &     ip, ir,
865     &     mp, mr, ndim
866
867      integer ::
868     &     ioff, idx1, idx2
869
870      ! get ipr:
871      if (mp.eq.1.and.mr.eq.1) then
872        ioff = 0
873        idx1 = max(ip,ir)
874        idx2 = min(ip,ir)
875      else if (mp.eq.1.and.mr.eq.2) then
876        ioff = ndim*(ndim-1)/2 ! NB: no diagonal in a/a block
877        idx1 = ir
878        idx2 = ip
879      else if (mp.eq.2.and.mr.eq.1) then
880        ioff = ndim*(ndim-1)/2
881        idx1 = ip
882        idx2 = ir
883      else if (mp.eq.2.and.mr.eq.2) then
884        ioff = ndim*(ndim-1)/2 + ndim*ndim
885        idx1 = max(ip,ir)
886        idx2 = min(ip,ir)
887      else
888        stop 'mp, mr ?'
889      end if
890
891      print *,'ioff,idx1,idx2: ',ioff,idx1,idx2
892
893      if (mp.eq.mr) then
894        if (ip.eq.ir) stop 'ip, ir ?'
895        if (idx1.eq.1) stop 'max(ip,ir) ?'
896        igtidx = ioff + (idx1-1)*(idx1-2)/2+idx2
897      else
898        igtidx = ioff + (idx1-1)*ndim+idx2
899      end if
900
901      print *,' >> gtidx = ',igtidx
902
903      return
904
905      end
906
907      subroutine cleanvec(vec,ndim)
908
909      implicit none
910
911      integer, intent(in) ::
912     &     ndim
913      real(8), intent(inout) ::
914     &     vec(ndim)
915
916      integer ::
917     &     ii
918      real(8) ::
919     &     thr
920
921      thr = 1000d0*epsilon(1d0)
922
923      print *,'thr = ',thr
924      print *,'ndim = ',ndim
925
926      do ii = 1, ndim
927        if (abs(vec(ii)).lt.thr) vec(ii)=0d0
928      end do
929
930      return
931
932      end
933c $Id$
934