1*----------------------------------------------------------------------*
2      subroutine occ_orbgrad(i_hybrid,ivcc,
3     &                       ccvec1,ccvec2,vec1,vec2,
4     &                       ittact,iooexcc,iooexc,
5     &                       nooexc,namp,nspin,
6     &                       xkapnrm,xkresnrm,xlresnrm,
7     &                       lu_ccamp,lu_lamp,lu_sig,
8     &                       lukappa,lu_ogrd,lu_odia,
9     &                       lu_1den,lu_2den,
10     &                       lu1int_o,lu2int_o,
11     &                       luintm1,luintm2,luref)
12*----------------------------------------------------------------------*
13*
14*     purpose: driver routine to calculate orbital gradient
15*
16*----------------------------------------------------------------------*
17      include "wrkspc.inc"
18      include "orbinp.inc"
19      include "glbbas.inc"
20      include "cgas.inc"
21      include "cands.inc"
22      include "cc_exc.inc"
23      include "oper.inc"
24      include "cstate.inc"
25      include "ctcc.inc"
26
27      integer, parameter ::
28     &     ntest = 50
29
30      integer, intent(in) ::
31     &     namp,
32     &     i_hybrid,ivcc,
33     &     lu_ccamp, lu_lamp, lu_sig,
34     &     lukappa, luintm1, luintm2, luref,
35     &     lu1int_o,lu2int_o,
36     &     iooexcc(*), iooexc(*), nooexc(nspin), ittact(*)
37      real(8), intent(inout) ::
38     &     ccvec1(*), ccvec2(*), vec1(*), vec2(*)
39
40      character(8) :: cctype
41
42      logical ::
43     &     l_ahap
44      real(8), external ::
45     &     inprod, inprdd
46
47      call atim(cpu0,wall0)
48
49      if (ntest.ge.10) then
50        write(6,*) 'Welcome to occ_orbgrad'
51        write(6,*) '======================'
52        write(6,*) ' namp = ',namp
53        write(6,*) ' i_hybrid = ',i_hybrid
54        write(6,*) ' lukappa, lu_ogrd = ',lukappa, lu_ogrd
55        write(6,*) ' lu_ccamp, lu_lamp = ',lu_ccamp, lu_lamp
56        write(6,*) ' luintm1, luintm2, luref = ',luintm1, luintm2, luref
57      end if
58
59      lblk = -1
60      cctype(1:6) = 'GEN_CC'
61      idum = 0
62      icspc = ietspc
63      icsm  = irefsm
64      isspc = ietspc
65      issm  = irefsm
66      mx_term = 100
67      icc_exc = 1
68      xconv = 1.0d-20
69
70      nooexc_tot = nooexc(1)
71      if (nspin.eq.2) nooexc_tot = nooexc_tot + nooexc(2)
72
73      call memman(idum,idum,'MARK  ',idum,'ORBGRD')
74
75      lu_lambda = iopen_nus('CCLAMBDA')
76      if (luintm1.lt.0.or.luintm2.lt.0) then
77        lusc1 = iopen_nus('CCDENSCR1')
78        lusc2 = iopen_nus('CCDENSCR2')
79        lusc3 = iopen_nus('CCDENSCR3')
80        lusc4 = iopen_nus('CCDENSCR4')
81      else
82        lusc3 = iopen_nus('CCDENSCR3')
83      end if
84
85* calculate 1 and 2 density
86      ! get right state
87      if (luintm1.gt.0) then
88        lurst = luintm1
89      else
90        lurst = -luintm1
91        call vec_from_disc(ccvec1,namp,1,lblk,lu_ccamp)
92        call expt_ref2(luref,lurst,lusc1,lusc2,lusc3,xconv,mx_term,
93     &               ccvec1,dum,vec1,vec2,namp,cctype,0)
94      end if
95
96      if (ivcc.eq.1) then
97        ! normalize state
98        xnorm = sqrt(inprdd(vec1,vec1,lurst,lurst,1,lblk))
99        call sclvcd(lurst,lusc3,1d0/xnorm,vec1,1,lblk)
100        lurst = lusc3
101      end if
102
103      if (ivcc.eq.0) then
104        ! get left state, if necessary
105        ! get refstate in ITSPC:
106        call expciv(irefsm,ietspc,luref,itspc,lusc3,
107     &            lblk,lusc4,1,0,idc,0)
108        if (luintm2.gt.0) then
109          ! add |ref> to obtain lambda
110          lulst = lu_lambda
111          call vecsmd(vec1,vec2,1d0,1d0,lusc3,luintm2,lu_lambda,1,lblk)
112        else
113          lulst = lu_lambda
114          ! sum(mu) L(mu) tau(mu) |HF>
115          call vec_from_disc(ccvec2,namp,1,lblk,lu_lamp)
116          icspc = itspc
117          isspc = itspc
118          call sig_gcc(vec1,vec2,lusc3,lusc2,ccvec2)
119          ! add |hf>
120          call vecsmd(vec1,vec2,1d0,1d0,lusc3,lusc2,lusc1,1,lblk)
121          ! and multiply with exp(-t\dag)
122          if (luintm1.gt.0)
123     &        call vec_from_disc(ccvec1,namp,1,lblk,lu_ccamp)
124          call scalve(ccvec1,-1d0,namp)
125          ! conjugate amplitudes (reorder)  and operators
126          call conj_ccamp(ccvec1,1,ccvec2)
127          call conj_t
128          ! exp(-t\dag), result on lulst
129          call expt_ref2(lusc1,lulst,lusc2,lusc3,lusc4,xconv,mx_term,
130     &       ccvec2,dum,vec1,vec2,namp,cctype,0)
131          call conj_t
132        end if
133      else
134        ! just get a copy of the right state
135        lulst = lu_lambda
136        call copvcd(lurst,lulst,vec1,1,lblk)
137      end if
138
139      ntbsq = ntoob*ntoob
140      lrho1 = nspin*ntbsq
141      lrho2 = nspin*ntbsq*(ntbsq+1)/2 + (nspin-1)*ntbsq**2
142
143      iden = 2
144      isepab = 0
145      if (nspin.eq.2) isepab = 1
146c      call densi2(iden,work(krho1),work(krho2),
147c     &     vec1,vec2,lulst,lurst,exps2,
148c     &     0,work(ksrho1))
149      if (ivcc.eq.0) then
150        icspc = ietspc
151        isspc = itspc
152      else
153        ! well, no games with VCC; it always needs the full (FCI) space
154        icspc = ietspc
155        isspc = ietspc
156      end if
157c dbg
158      print *,'call to densi2_ab'
159      call flush(6)
160c dbg
161      call densi2_ab(iden,work(krho1),work(krho2),
162     &     vec1,vec2,lulst,lurst,exps2,
163     &     0,work(ksrho1),isepab)
164c dbg
165      print *,'after densi2_ab'
166      call flush(6)
167c dbg
168
169      ! symmetrize one-body density matrix
170      do ispin = 1, nspin
171        ioff = (ispin-1)*ntbsq
172        call sym_blmat(work(krho1+ioff),1,ntoob)
173      end do
174      ! symmetrize two-body density matrix
175      do ispc = 1, nspin*2-1
176        ioff = (ispc-1)*ntbsq*(ntbsq+1)/2
177        imod = 0
178        if (ispc.eq.3) imod = 1
179
180        call sym_2dens(work(krho2+ioff),ntoob,imod,0.01d0)
181      end do
182
183      if (lu_1den.gt.0)
184     &     call vec_to_disc(work(krho1),lrho1,1,lblk,lu_1den)
185      if (lu_2den.gt.0) then
186        call vec_to_disc(work(krho2),lrho2,1,lblk,lu_2den)
187      end if
188
189* calculate eff. fock matrix
190      icc_exc = 0
191      call memman(klfoo,nspin*ntoob**2,'ADDL  ',2,'FEFF  ')
192      call fock_mat_ab(work(klfoo),2,nspin)
193c      call fock_mat(work(klfoo),2)
194* get orbital gradient
195      call memman(kle1,nooexc_tot,'ADDL  ',2,'E1    ')
196
197      ! Hybrid: the p/h gradient is already on disc, read it
198      if (i_hybrid.eq.1)
199     &     call vec_from_disc(work(kle1),nooexc_tot,1,lblk,lu_ogrd)
200
201* orbital gradient types:
202* 2 :  general orbital gradient at kappa != 0
203* 1 :  CEP (current expansion point) orbital gradient
204* 3 :  numerical general orbital gradient
205      iogtyp = 1
206
207      if (iogtyp.eq.1) then
208        do ispin = 1, nspin
209          ifoo = (ispin-1)*ntoob**2
210          ie1  = (ispin-1)*nooexc(1)
211          imode = 0
212          if (i_hybrid.eq.1) imode=-1
213          call f_to_e1(work(klfoo+ifoo),work(kle1+ie1),
214     &         1,iooexcc(2*ie1+1),
215     &         nooexc(ispin),imode)
216          if (ntest.ge.100) then
217            if (i_hybrid.eq.1)
218     &           write(6,*) 'After adding inactive/active rotations:'
219            if (nspin.eq.1)
220     &           write(6,*) 'The orbital gradient:'
221            if (nspin.eq.2.and.ispin.eq.1)
222     &           write(6,*) 'The orbital gradient (alpha):'
223            if (nspin.eq.2.and.ispin.eq.2)
224     &           write(6,*) 'The orbital gradient (beta):'
225            call wrt_excvec(work(kle1+ie1),
226     &           iooexcc(2*ie1+1),nooexc(ispin))
227          end if
228        end do
229
230        if (nspin.eq.2) then
231          iexc2 = nooexc(1)
232          do iexc = 1, nooexc(1)
233            iorb = iooexcc((iexc-1)*2+1)
234            jorb = iooexcc((iexc-1)*2+2)
235            itp = itpfto(iorb)
236            jtp = itpfto(jorb)
237            l_ahap = i_iadx(itp).eq.2.and.i_iadx(jtp).eq.2.and.
238     &           ihpvgas(itp).ne.ihpvgas(jtp)
239            if (.not.l_ahap) then
240              do
241                iexc2 = iexc2+1
242                if (iexc2.gt.nooexc_tot) stop 'oha!'
243                iorb2 = iooexcc((iexc2-1)*2+1)
244                jorb2 = iooexcc((iexc2-1)*2+2)
245                if (iorb.eq.iorb2.and.jorb.eq.jorb2) then
246                  avg = 0.5d0*(work(kle1-1+iexc)+work(kle1-1+iexc2))
247                  work(kle1-1+iexc) = avg
248                  work(kle1-1+iexc2)= avg
249                  exit
250                end if
251              end do
252            end if
253          end do
254          if (ntest.ge.100) then
255            if (i_hybrid.eq.1)
256     &           write(6,*) 'After averaging non-ap/ah rotations:'
257            do ispin = 1, nspin
258              ie1  = (ispin-1)*nooexc(1)
259              if (nspin.eq.1)
260     &             write(6,*) 'The orbital gradient:'
261              if (nspin.eq.2.and.ispin.eq.1)
262     &             write(6,*) 'The orbital gradient (alpha):'
263              if (nspin.eq.2.and.ispin.eq.2)
264     &             write(6,*) 'The orbital gradient (beta):'
265              call wrt_excvec(work(kle1+ie1),
266     &             iooexcc(2*ie1+1),nooexc(ispin))
267            end do
268          end if
269
270        end if
271
272      else if (iogtyp.eq.2) then
273        stop 'do your work'
274      else if (iogtyp.eq.3) then
275cedo disabled since routine is missing
276        write(6,*) 'orbgrd_num2 missing '
277        stop 'orbgrd_num2'
278c        call orbgrd_num2(0,work(kle1),lukappa,
279c     &       lu1int_o,lu2int_o,iooexcc,nooexc)
280      else
281        write(6,*) 'unknown iogtyp (',iogtyp,')'
282        stop 'occ_orbgrad'
283      end if
284      xkresnrm = sqrt(inprod(work(kle1),work(kle1),nooexc_tot))
285
286      if (ntest.ge.10) then
287        do ispin = 1, nspin
288          ie1  = (ispin-1)*nooexc(1)
289          if (nspin.eq.1)
290     &         write(6,*) 'The orbital gradient:'
291          if (nspin.eq.2.and.ispin.eq.1)
292     &         write(6,*) 'The orbital gradient (alpha):'
293          if (nspin.eq.2.and.ispin.eq.2)
294     &         write(6,*) 'The orbital gradient (beta):'
295          call wrt_excvec(work(kle1+ie1),
296     &         iooexcc(2*ie1+1),nooexc(ispin))
297        end do
298      end if
299
300      call vec_to_disc(work(kle1),nooexc_tot,1,lblk,lu_ogrd)
301
302      call relunit(lu_lambda,'delete')
303
304      if (luintm1.lt.0.or.luintm2.lt.0) then
305        call relunit(lusc1,'delete')
306        call relunit(lusc2,'delete')
307        call relunit(lusc3,'delete')
308        call relunit(lusc4,'delete')
309      else
310        call relunit(lusc3,'delete')
311      end if
312
313      ! add orbital-gradient contributions to residual of
314      ! left-hand vector
315      if (i_hybrid.eq.1) then
316        call memman(ko1c,ntoob,'ADDL  ',1,'OGRD C')
317        call memman(ko1a,ntoob,'ADDL  ',1,'OGRD A')
318
319        call vec_from_disc(ccvec1,namp,1,lblk,lu_sig)
320
321        do ispin = 1, nspin
322          call compress_t1s(ccvec1,work(klfoo),1,work(ko1c),work(ko1a),
323     &     work(klsobex),work(klsox_to_ox),
324     &     work(klcobex_tp),work(klibsobex),
325     &     nspobex_tp,ispin,-2)
326        end do
327
328        xlresnrm = sqrt(inprod(ccvec1,ccvec1,namp))
329
330        call vec_to_disc(ccvec1,namp,1,lblk,lu_sig)
331
332      end if
333
334      iaprhess=1
335c      iaprhess=0
336      if (lu_odia.le.0) iaprhess=0
337      if (iaprhess.eq.1) then
338CCCCC
339c          call memman(kle3,nooexc_tot,'ADDL  ',2,'NT KAP')
340c          call memman(kle4,nooexc_tot,'ADDL  ',2,'NT KAP')
341CCCCC
342        imode = 0
343        do ispin = 1, nspin
344          ifoo = (ispin-1)*ntoob**2
345          ie1  = (ispin-1)*nooexc(1)
346          itt  = (ispin-1)*ngas**2
347
348          if (nspin.eq.1) ispc = 0
349          if (nspin.eq.2) ispc = ispin
350          call diag_orbhes(work(kle1+ie1),work(klfoo+ifoo),
351     &         iooexc(1+ifoo),nooexc(ispin),1,ittact(1+itt),ispc)
352CCCCC
353c          xinc = 0.00001d0
354c          luoinc = iopen_uus()
355c          lumodu = iopen_uus()
356c
357c          idx = 0
358c          if (ispin.eq.2) idx = nooexc(1)
359c          do iexc = 1, nooexc(ispin)
360c            idx = idx+1
361c            ! increment kappa
362c            work(kle3:kle3-1+nooexc_tot) = 0d0
363c            work(kle3-1+idx) = xinc
364c            call vec_to_disc(work(kle3),nooexc_tot,1,-1,luoinc)
365c            ! get transf.-mat. from current kappa
366c            call kap2u(3,-1000,lukappa,lumodu,iooexcc,nooexc,nspin)
367c            ! update transf.mat.
368c            call kap2u(-3,luoinc,-1000,lumodu,iooexcc,nooexc,nspin)
369c            ! get new integrals with mod. transf.mat.
370c            call tra_kappa(-1,lumodu,iooexcc,nooexc,nspin,
371c     &                      1,lu1int_o,lu2int_o)
372c            ! get Fock
373c            call fock_mat_ab(work(klfoo),2,nspin)
374c            ! get E1(+)
375c            do jspin = 1, nspin
376c              jfoo = (jspin-1)*ntoob**2
377c              je1  = (jspin-1)*nooexc(1)
378c              imode = 0
379c              call f_to_e1(work(klfoo+jfoo),work(kle3+je1),
380c     &             1,iooexcc(2*ie1+1),
381c     &             nooexc(jspin),imode)
382c            end do
383c
384c            ! increment kappa -
385c            work(kle4:kle4-1+nooexc_tot) = 0d0
386c            work(kle4-1+idx) = -xinc
387c            call vec_to_disc(work(kle4),nooexc_tot,1,-1,luoinc)
388c            ! get transf.-mat. from current kappa
389c            call kap2u(3,-1000,lukappa,lumodu,iooexcc,nooexc,nspin)
390c            ! update transf.mat.
391c            call kap2u(-3,luoinc,-1000,lumodu,iooexcc,nooexc,nspin)
392c            ! get new integrals with mod. transf.mat.
393c            call tra_kappa(-1,lumodu,iooexcc,nooexc,nspin,
394c     &                      1,lu1int_o,lu2int_o)
395c            ! get Fock
396c            call fock_mat_ab(work(klfoo),2,nspin)
397c            ! get E1(-)
398c            do jspin = 1, nspin
399c              jfoo = (jspin-1)*ntoob**2
400c              je1  = (jspin-1)*nooexc(1)
401c              imode = 0
402c              call f_to_e1(work(klfoo+jfoo),work(kle4+je1),
403c     &             1,iooexcc(2*je1+1),
404c     &             nooexc(jspin),imode)
405c            end do
406c
407c            iorb = iooexcc((idx-1)*2+1)
408c            jorb = iooexcc((idx-1)*2+2)
409c
410c            xdia = (work(kle3-1+idx)-work(kle4-1+idx))/(2d0*xinc)
411c            print *,' i,j: ',iorb,jorb
412c            print *,' numerical  ',xdia
413c            print *,' analytical ',work(kle1-1+idx)
414c
415c            work(kle1-1+idx) = xdia
416c          end do
417c
418c          ! restore law and order:
419c          call kap2u(3,-1000,lukappa,lumodu,iooexcc,nooexc,nspin)
420c          ! update transf.mat.
421c          call tra_kappa(-1,lumodu,iooexcc,nooexc,nspin,
422c     &                      1,lu1int_o,lu2int_o)
423c
424c          call relunit(luoinc,'delete')
425c          call relunit(lumodu,'delete')
426CCCCC
427          if (ntest.ge.100) then
428            write(6,*) 'The diagonal orbital Hessian (unmodified):'
429            if (ispin.eq.1.and.nspin.eq.2) write(6,*) '   alpha part'
430            if (ispin.eq.2.and.nspin.eq.2) write(6,*) '   beta part'
431            call wrt_excvec(work(kle1+ie1),
432     &                      iooexcc(1+2*ie1),nooexc(ispin))
433          end if
434        end do
435
436        if (nspin.eq.2) then
437          iexc2 = nooexc(1)
438          do iexc = 1, nooexc(1)
439            iorb = iooexcc((iexc-1)*2+1)
440            jorb = iooexcc((iexc-1)*2+2)
441            itp = itpfto(iorb)
442            jtp = itpfto(jorb)
443            l_ahap = i_iadx(itp).eq.2.and.i_iadx(jtp).eq.2.and.
444     &           ihpvgas(itp).ne.ihpvgas(jtp)
445            if (.not.l_ahap) then
446              do
447                iexc2 = iexc2+1
448                if (iexc2.gt.nooexc_tot) stop 'oha2!'
449                iorb2 = iooexcc((iexc2-1)*2+1)
450                jorb2 = iooexcc((iexc2-1)*2+2)
451                if (iorb.eq.iorb2.and.jorb.eq.jorb2) then
452                  avg = 0.5d0*(work(kle1-1+iexc)+work(kle1-1+iexc2))
453                  work(kle1-1+iexc) = avg
454                  work(kle1-1+iexc2)= avg
455                  exit
456                end if
457              end do
458            end if
459          end do
460        end if
461
462        xmin = 1d15
463        do iexc = 1, nooexc_tot
464c          work(kle1-1+iexc) = work(kle1-1+iexc)
465c          if (abs(work(kle1-1+iexc)).le.1d-4) work(kle1-1+iexc)=1d+13
466          if (work(kle1-1+iexc).lt.-0.2d0)
467     &         work(kle1-1+iexc)=-work(kle1-1+iexc)
468          xmin = min(xmin,work(kle1-1+iexc))
469        end do
470        xsh = 0d0
471        xshm = 2d-2
472        if (xmin.lt.xshm) then
473          xsh = xshm - xmin!2d-1 - xmin
474          work(kle1:kle1-1+nooexc_tot) =
475     &         work(kle1:kle1-1+nooexc_tot)+xsh
476        end if
477        if (ntest.ge.50) then
478          write(6,*) 'The diagonal orbital Hessian (modified):'
479          write(6,*) ' shift was: ',xsh
480          do ispin = 1, nspin
481            ie1  = (ispin-1)*nooexc(1)
482            if (ispin.eq.1.and.nspin.eq.2) write(6,*) '   alpha part'
483            if (ispin.eq.2.and.nspin.eq.2) write(6,*) '   beta part'
484            call wrt_excvec(work(kle1+ie1),
485     &                      iooexcc(1+2*ie1),nooexc(ispin))
486          end do
487        end if
488        call vec_to_disc(work(kle1),nooexc_tot,1,lblk,lu_odia)
489      end if
490
491      call memman(idum,idum,'FLUSM  ',idum,'ORBGRD')
492
493      call atim(cpu,wall)
494
495      call prtim(6,'time for orbital gradient',cpu-cpu0,wall-wall0)
496
497      return
498
499      end
500*----------------------------------------------------------------------*
501*----------------------------------------------------------------------*
502      subroutine mdfy_hss(xhss,ittact,iooexcc,itpfto,
503     &                    ihpvgas_ab,ld_hpv,
504     &                    nooexc,ngas,nspin)
505*----------------------------------------------------------------------*
506*     a diagonal hessian approximated as
507*
508*      H(ij,ij) = F_ii - F_jj
509*
510*     is input. Some empirical rules are used to modify i.pt.
511*     inactive/active rotations
512*
513*----------------------------------------------------------------------*
514      implicit none
515
516      integer, intent(in) ::
517     &     nooexc(2),ngas,nspin,ld_hpv,
518     &     ittact(ngas,ngas,nspin),iooexcc(2,*),itpfto(*),
519     &     ihpvgas_ab(ld_hpv,nspin)
520      real(8), intent(inout) ::
521     &     xhss(*)
522
523      integer, parameter ::
524     &     ntest = 00
525
526      integer ::
527     &     nooexc_tot,
528     &     ispin, iexc, idx, iorb, jorb, itp, jtp
529
530      nooexc_tot = nooexc(1)
531      if (nspin.eq.2) nooexc_tot = nooexc_tot + nooexc(2)
532
533      idx = 0
534      do ispin = 1, nspin
535c        iooff = (ispin-1)*nooexc(1)
536        do iexc = 1, nooexc(ispin)
537          idx = idx + 1
538          iorb = iooexcc(1,idx)
539          jorb = iooexcc(2,idx)
540          itp = itpfto(iorb)
541          jtp = itpfto(jorb)
542          if (ihpvgas_ab(itp,ispin).eq.
543     &        ihpvgas_ab(jtp,ispin)    ) then
544            xhss(idx) = 0.025d0 * xhss(idx)
545          end if
546          if (xhss(idx).lt.0d0) xhss(idx)=0.1d0
547
548        end do
549      end do
550
551      if (ntest.ge.100) then
552        write(6,*) 'Modified approx. orbital hessian:'
553        idx = 0
554        do ispin = 1, nspin
555          if (nspin.eq.2.and.ispin.eq.1) write(6,*) 'alpha:'
556          if (nspin.eq.2.and.ispin.eq.2) write(6,*) 'beta:'
557          do iexc = 1, nooexc(ispin)
558            idx = idx+1
559            write(6,'(x,2i5,g20.10)')
560     &           iooexcc(1,idx),
561     &           iooexcc(2,idx),
562     &           xhss(idx)
563          end do
564        end do
565      end if
566
567      return
568
569      end
570*----------------------------------------------------------------------*
571*----------------------------------------------------------------------*
572      subroutine new_fock(fock,rho1,diag_h,nooexc,iooexcc)
573*----------------------------------------------------------------------*
574
575      include "wrkspc.inc"
576      include "glbbas.inc"
577      include "cintfo.inc"
578      include "orbinp.inc"
579      include "lucinp.inc"
580
581      integer, parameter ::
582     &     ntest = 00
583
584      real(8), intent(inout) ::
585     &     fock(*), rho1(*), diag_h(*)
586      integer, intent(in) ::
587     &     iooexcc(2,*), nooexc
588
589      integer ::
590     &     ioff(nsmob)
591
592      call swapve(work(krho1),rho1,ntoob*ntoob)
593      call copvec(work(kint1),fock,nint1)
594      call fifam(fock)
595
596      if (ntest.ge.100) then
597        write(6,*) 'new fock matrix:'
598        call aprblm2(fock,ntoobs,ntoobs,nsmob,1)
599      end if
600
601      ! precalculate offsets of symmetry blocks
602      idx = 0
603      do ism = 1, nsmob
604        ioff(ism) = idx
605        idx = idx + (ntoobs(ism)+1)*ntoobs(ism)/2
606      end do
607
608      do iexc = 1, nooexc
609        idx = iooexcc(1,iexc)
610        jdx = iooexcc(2,iexc)
611        ! convert type to symmetry ordering:
612        ism = ismfto(idx)
613        jsm = ismfto(jdx)
614        ii = ireots(idx) - ibso(ism) + 1
615        jj = ireots(jdx) - ibso(jsm) + 1
616        ! diagonal element addresses in *upper* triangles
617        iid = ioff(ism) + (ii+1)*ii/2
618        jjd = ioff(jsm) + (jj+1)*jj/2
619
620          print *,'iexc, idx, jdx, ii, jj: ',iexc, idx, jdx, ii, jj
621          print *,'   ism, jsm, ndimi, ndimj: ', ism, jsm, ibso(ism)
622          print *,'   iid, jjd: ', iid, jjd
623          print *,'             ',fock(iid) , fock(jjd)
624
625        diag_h(iexc) =  4d0 * (fock(iid) - fock(jjd))
626      end do
627      if (ntest.ge.50) then
628        write(6,*) 'diagonal Orbital Hessian:'
629        call wrtmat(diag_h,nooexc,1,nooexc,1)
630      end if
631
632      call swapve(work(krho1),rho1,ntoob*ntoob)
633
634      return
635      end
636*----------------------------------------------------------------------*
637*----------------------------------------------------------------------*
638      subroutine omg2ogrd(luomg,lu_ogrd,ccvec1,
639     &             iooexcc,nooexc,n_cc_amp,nspin,
640     &             xkresnrm,xomgnrm)
641*----------------------------------------------------------------------*
642*
643*     rearrage single-excitation part of Omega as orbital gradient
644*     for Brueckner CC
645*
646*----------------------------------------------------------------------*
647      include "wrkspc.inc"
648      include "orbinp.inc"
649      include "lucinp.inc"
650      include "ctcc.inc"
651
652      integer, parameter ::
653     &     ntest = 00
654
655      integer, intent(in) ::
656     &     luomg, lu_ogrd,
657     &     nooexc(2), n_cc_amp, iooexcc(*), nspin
658      real(8), intent(inout) ::
659     &     ccvec1(n_cc_amp)
660      real(8), intent(out) ::
661     &     xkresnrm,xomgnrm
662
663      real(8), external ::
664     &     inprod
665
666      if (ntest.ge.10) then
667        write(6,*) 'OMG2OGRD'
668        write(6,*) '========'
669        write(6,*) ' luomg, lu_ogrd: ',luomg, lu_ogrd
670        write(6,*) ' nooexc, n_cc_amp: ',nooexc, n_cc_amp, nspin
671      end if
672
673      lblk = -1
674      idum = 0
675      call memman(idum,idum,'MARK  ',idum,'OM2OGR')
676
677      nooexc_tot = nooexc(1)
678      if (nspin.eq.2) nooexc_tot = nooexc_tot + nooexc(2)
679
680      ! get workspace
681      logrd = 0
682      do ism = 1, nsmob
683        logrd = logrd + ntoobs(ism)*ntoobs(ism)
684      end do
685
686      call memman(kogrd,logrd,'ADDL  ',2,'OGRD  ')
687      call memman(kogrdc,nooexc_tot,'ADDL  ',2,'OGRDC ')
688      call memman(ko1c,ntoob,'ADDL  ',1,'OGRD C')
689      call memman(ko1a,ntoob,'ADDL  ',1,'OGRD A')
690
691      ! read from disc
692      call vec_from_disc(ccvec1,n_cc_amp,1,lblk,luomg)
693
694      ! and sort singles into new array
695      do ispin = 1, nspin
696c      iab = 1 ! prelim
697
698        call expand_t1s_new(ccvec1,work(kogrd),1,work(ko1c),work(ko1a),
699     &     work(klsobex),work(klsox_to_ox),
700     &     work(klcobex_tp),work(klibsobex),
701     &     nspobex_tp,ispin)
702
703        ! compress orbital gradient
704        ioff = (ispin-1)*nooexc(1)*2
705        koff = (ispin-1)*nooexc(1)
706        call comprs_kap(work(kogrd),work(kogrdc+koff),
707     &       iooexcc(ioff+1),nooexc(ispin))
708
709      end do
710
711      ! save modified omega to disc
712      call zero_t1(ccvec1)
713      xomgnrm = sqrt(inprod(ccvec1,ccvec1,n_cc_amp))
714
715      call vec_to_disc(ccvec1,n_cc_amp,1,lblk,luomg)
716
717      ! scale it with factor
718      if (nspin.eq.1) fac = 4d0
719      if (nspin.eq.2) fac = 2d0
720      call scalve(work(kogrdc),fac,nooexc_tot)
721
722      xkresnrm = sqrt(inprod(work(kogrdc),work(kogrdc),nooexc_tot))
723
724      ! save orbital gradient to disc
725      call vec_to_disc(work(kogrdc),nooexc_tot,1,lblk,lu_ogrd)
726
727      if (ntest.ge.100) then
728        write(6,*) ' Brueckner orbital gradient'
729        write(6,*) ' =========================='
730        do ispin = 1, nspin
731          ioff = (ispin-1)*nooexc(1)*2
732          koff = (ispin-1)*nooexc(1)
733          if (ispin.eq.1.and.nspin.eq.2) write(6,*) '  alpha'
734          if (ispin.eq.2.and.nspin.eq.2) write(6,*) '  beta'
735          call wrt_excvec(work(kogrdc+koff),iooexcc(ioff+1),
736     &                    nooexc(ispin))
737        end do
738      end if
739
740      idum = 0
741      call memman(idum,idum,'FLUSM ',idum,'OM2OGR')
742
743      return
744      end
745
746      subroutine occ_scnd_num(i_obcc,
747     &           ccvec1,ccvec2,vec1,vec2,
748     &           ittact,iooexcc,iooexc,
749     &           nooexc,n_cc_amp,nspin,
750     &           lu1into,lu2into,luref,
751     &           luccamp,lu_lamp,lukappa,
752     &           lutrvec,lutrv_l,lutrv_o,
753     &           lursig,lusig_l,lusig_o,
754     &           iactt,iactl,iacto)
755
756      include 'implicit.inc'
757
758      integer, intent(in) ::
759     &     nooexc(nspin)
760      real(8), intent(inout) ::
761     &     ccvec1(n_cc_amp), ccvec2(n_cc_amp)
762
763      real(8) ::
764     &     xinc(2)
765      integer ::
766     &     luomg(2), lugrd(2), luogrd(2)
767
768      real(8), external ::
769     &     inprdd
770
771      ninc = 2
772      delta = 1d-5
773      xinc(1) = delta
774      xinc(2) = -delta
775
776      lumodt = iopen_nus('NUMMODT')
777      lumodl = iopen_nus('NUMMODL')
778      lumodo = iopen_nus('NUMMODO')
779
780      luomg(1) = iopen_nus('OMGINC1')
781      luomg(2) = iopen_nus('OMGINC2')
782      lugrd(1) = iopen_nus('GRDINC1')
783      lugrd(2) = iopen_nus('GRDINC2')
784      luogrd(1) = iopen_nus('OGRDINC1')
785      luogrd(2) = iopen_nus('OGRDINC2')
786
787      luintm1 = iopen_nus('NUMCCINTM1')
788      luintm2 = iopen_nus('NUMCCINTM2')
789      luintm3 = iopen_nus('NUMCCINTM3')
790
791      lurhs = iopen_nus('NUM_RHS')
792
793      write(6,*) ' NUMERICAL HESSIAN FOR OCC:'
794
795      xnt2 = inprdd(ccvec1,ccvec1,lutrvec,lutrvec,1,-1)
796      xnl2 = inprdd(ccvec1,ccvec1,lutrv_l,lutrv_l,1,-1)
797      xno2 = inprdd(ccvec1,ccvec1,lutrv_o,lutrv_o,1,-1)
798      write(6,*) ' NORM OF INPUT VECTOR: ',sqrt(xnt2+xnl2+xno2)
799      write(6,*) ' NORM OF COMPONENTS:   ',
800     &     sqrt(xnt2),sqrt(xnl2),sqrt(xno2)
801
802      do iinc = 1, ninc
803        write(6,*) ' OUTPUT FOR XINC = ',XINC(IINC),' FOLLOWS:'
804        ! modify parameters by increment
805        call vecsmd(ccvec1,ccvec2,1d0,xinc(iinc),
806     &       luccamp,lutrvec,lumodt,1,-1)
807        call vecsmd(ccvec1,ccvec2,1d0,xinc(iinc),
808     &       lu_lamp,lutrv_l,lumodl,1,-1)
809c V A:
810c        call vecsmd(ccvec1,ccvec2,1d0,xinc(iinc),
811c     &       lukappa,lutrv_o,lumodo,1,-1)
812c
813c V B:
814        ! get transf.-mat. from current kappa
815        call kap2u(3,-1000,lukappa,lumodo,iooexcc,nooexc,nspin)
816        call sclvcd(lutrv_o,luintm1,xinc(iinc),ccvec1,1,-1)
817        ! update transf.mat.
818        call kap2u(-3,luintm1,-1000,lumodo,iooexcc,nooexc,nspin)
819c
820
821        ! get new orbitals with mod. transf.mat.
822        call tra_kappa(-1,lumodo,iooexcc,nooexc,nspin,
823     &                 1,lu1into,lu2into)
824
825        ! call vector function
826        call cc_vec_fnc2(ccvec1,ccvec2,ecc,eccl,
827     &                   xampnrm,xomgnrm,xdum,
828     &                   vec1,vec2,1,'GEN_CC',
829     &                   xdum,
830     &                   lumodt,luomg(iinc),lumodl,
831     &                   luintm1,luintm2,luintm3)
832
833        ! call rhs and jacobian lh trafo
834        call zero_ord_rhs(ccvec1,vec1,vec2,lumodt,lurhs,luintm1)
835        lr_switch = 1
836        iadd_rhs = 1
837        itex_sm = 1
838        call jac_t_vec2(lr_switch,iadd_rhs,0,1,1,
839     &           ccvec1,ccvec2,vec1,vec2,
840     &           n_cc_amp,n_cc_amp,
841     &           ecc1,xlampnrm,xlresnrm,
842     &           lumodt,luomg(iinc),lumodl,lugrd(iinc),lurhs,
843     &           luintm1,luintm2,luintm3)
844
845        ! call orbital gradient
846        if (iacto.eq.0) then
847          namp = sum(nooexc(1:nspin))
848          ccvec1(1:namp) = 0d0
849          call vec_to_disc(ccvec1,namp,1,-1,luogrd(iinc))
850        else
851          if (i_obcc.eq.1) then
852            call omg2ogrd(luomg(iinc),luogrd(iinc),ccvec1,
853     &           iooexcc,nooexc,n_cc_amp,nspin,
854     &           xkresnrm,xomgnrm)
855          end if
856          ivcc = 0
857          call occ_orbgrad(i_obcc,ivcc,
858     &         ccvec1,ccvec2,vec1,vec2,
859     &         ittact,iooexcc,iooexc,
860     &         nooexc,n_cc_amp,nspin,
861     &         xkapnrm,xkresnrm,xlresnrm,
862     &         lumodt,lumodl,lugrd(iinc),
863     &         ludum,luogrd(iinc),-1,
864     &         -1,-1,
865     &         lu1into,lu2into,
866     &         luintm1,luintm3,luref)
867        end if
868
869      end do
870
871      ! calculate numerical sigma-vectors
872      fac = 1d0/(2d0*delta)
873      call vecsmd(ccvec1,ccvec2,fac,-fac,luomg(1),luomg(2),lursig,1,-1)
874      call vecsmd(ccvec1,ccvec2,fac,-fac,lugrd(1),lugrd(2),lusig_l,1,-1)
875c TEST --- no t/l contribution
876c          do ii = 1, 10
877c            print *,'T,L set to 0D0 !!!'
878c          end do
879c          ccvec1(1:n_cc_amp)=0d0
880c          call vec_to_disc(ccvec1,n_cc_amp,1,-1,lursig)
881c          call vec_to_disc(ccvec1,n_cc_amp,1,-1,lusig_l)
882c TEST
883      call vecsmd(ccvec1,ccvec2,fac,-fac,
884     &     luogrd(1),luogrd(2),lusig_o,1,-1)
885c TEST --- no kappa contribution
886c          if (nspin.eq.2) stop 'not possible'
887c          do ii = 1, 10
888c            print *,'kappa set to 0D0 !!!'
889c          end do
890c          ccvec1(1:nooexc(1))=0d0
891c          call vec_to_disc(ccvec1,nooexc(1),1,-1,lusig_o)
892c TEST
893
894      xnt2 = inprdd(ccvec1,ccvec1,lursig,lursig,1,-1)
895      xnl2 = inprdd(ccvec1,ccvec1,lusig_l,lusig_l,1,-1)
896      xno2 = inprdd(ccvec1,ccvec1,lusig_o,lusig_o,1,-1)
897      write(6,*) ' NORM OF OUTPUT VECTOR: ',sqrt(xnt2+xnl2+xno2)
898      write(6,*) ' NORM OF COMPONENTS:   ',
899     &     sqrt(xnt2),sqrt(xnl2),sqrt(xno2)
900
901      ! delete files
902      call relunit(lumodt,'delete')
903      call relunit(lumodl,'delete')
904      call relunit(lumodo,'delete')
905
906      call relunit(luomg(1),'delete')
907      call relunit(luomg(2),'delete')
908      call relunit(lugrd(1),'delete')
909      call relunit(lugrd(2),'delete')
910      call relunit(luogrd(1),'delete')
911      call relunit(luogrd(2),'delete')
912
913      call relunit(luintm1,'delete')
914      call relunit(luintm2,'delete')
915      call relunit(luintm3,'delete')
916
917      call relunit(lurhs,'delete')
918
919      ! restore the old orbitals
920      call tra_kappa(lukappa,-1,iooexcc,nooexc,nspin,
921     &               1,lu1into,lu2into)
922
923      end
924c $Id$
925