1      subroutine tce_mrcc_jacobi_t2(d_r2,d_t2,k_t2_offset,iter,iref)
2c
3      implicit none
4#include "global.fh"
5#include "mafdecls.fh"
6#include "sym.fh"
7#include "util.fh"
8#include "stdio.fh"
9#include "errquit.fh"
10#include "tce.fh"
11#include "tce_main.fh"
12#include "tce_mrcc.fh"
13      integer d_r2
14      integer d_t2
15      integer p1b
16      integer p2b
17      integer h3b
18      integer h4b
19      integer p1
20      integer p2
21      integer h3
22      integer h4
23      integer k_t2_offset
24      integer size
25      integer l_r2,k_r2
26      integer l_t2,k_t2
27      integer i
28      integer nprocs
29      integer count
30      integer next
31cc      integer nxtask
32c
33      integer iter
34      integer iref
35      double precision denom,cit
36      double precision d1,d2,d3,d4
37      integer orbspin1,orbspin2,orbspin3,orbspin4
38      integer orbindex1,orbindex2,orbindex3,orbindex4
39      integer blck1,blck2,blck3,blck4
40      integer off1,off2,off3,off4
41c
42cc      external nxtask
43      INTEGER NXTASK
44      EXTERNAL NXTASK
45      logical nodezero
46      logical noloadbalance
47c *** shift ***
48      double precision shift
49c *************
50      if(iter.le.100) then
51       shift=-(2.0d0)*zlshift
52      else
53       shift=0.0d0
54      end if
55c
56c     ================
57c     Loop over blocks
58c     ================
59c
60      nodezero = (ga_nodeid().eq.0)
61      noloadbalance = ((ioalg.eq.4).or.
62     1                ((ioalg.eq.6).and.(.not.fileisga(d_r2))))
63      nprocs = ga_nnodes()
64      count = 0
65cc      next = nxtask(nprocs,1)
66      call ga_sync()
67      next = NXTASK(nprocs, 1)
68      do p1b = noab+1,noab+nvab
69        do p2b = p1b,noab+nvab
70          do h3b = 1,noab
71            do h4b = h3b,noab
72              if (noloadbalance.or.(next.eq.count)) then
73                if (int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1)
74     1            .eq. int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1)) then
75                if ((.not.restricted).or.
76     1            (int_mb(k_spin+p1b-1)+int_mb(k_spin+p2b-1)+
77     2            int_mb(k_spin+h3b-1)+int_mb(k_spin+h4b-1).ne.8)) then
78                if (ieor(int_mb(k_sym+p1b-1),ieor(int_mb(k_sym+p2b-1),
79     1            ieor(int_mb(k_sym+h3b-1),int_mb(k_sym+h4b-1))))
80     2            .eq. 0) then
81                  size = int_mb(k_range+p1b-1) * int_mb(k_range+p2b-1)
82     1                 * int_mb(k_range+h3b-1) * int_mb(k_range+h4b-1)
83
84                  if (.not.ma_push_get(mt_dbl,size,'r2',l_r2,k_r2))
85     1              call errquit('tce_jacobi_t2: MA problem',0,MA_ERR)
86                  call get_hash_block(d_r2,dbl_mb(k_r2),size,
87     1              int_mb(k_t2_offset),((((p1b-noab-1)*nvab+p2b-noab-1)
88     2              *noab+h3b-1)*noab+h4b-1))
89c                if(lsubterm) then
90                  if (.not.ma_push_get(mt_dbl,size,'t2',l_t2,k_t2))
91     1              call errquit('tce_jacobi_t2: MA problem',0,MA_ERR)
92                  call get_hash_block(d_t2,dbl_mb(k_t2),size,
93     1              int_mb(k_t2_offset),((((p1b-noab-1)*nvab+p2b-noab-1)
94     2              *noab+h3b-1)*noab+h4b-1))
95c                endif
96
97                  i = 0
98c *** shift added ***
99ccc               shift=0.50d0
100c *******************
101                  do p1 = 1,int_mb(k_range+p1b-1)
102                    do p2 = 1,int_mb(k_range+p2b-1)
103                      do h3 = 1,int_mb(k_range+h3b-1)
104                        do h4 = 1,int_mb(k_range+h4b-1)
105                          i = i + 1
106
107        denom = (-dbl_mb(k_evl_sorted+int_mb(k_offset+p1b-1)+p1-1)
108     2  -dbl_mb(k_evl_sorted+int_mb(k_offset+p2b-1)+p2-1)
109     3  +dbl_mb(k_evl_sorted+int_mb(k_offset+h3b-1)+h3-1)
110     4  +dbl_mb(k_evl_sorted+int_mb(k_offset+h4b-1)+h4-1)
111     5  +mrccshift)
112
113c       if((abs(denom).lt.0.01d0).and.
114c     1 (abs(dbl_mb(k_r2+i-1)/denom).gt.1.0d0))
115c     1  then
116c       write(6,"('2DENOM CLOSE TO ZERO: ',F16.10,F16.10,F16.10)")
117c     1 denom,dbl_mb(k_t2+i-1),dbl_mb(k_r2+i-1)
118c       endif
119
120       if(lsubterm) then
121        cit = dbl_mb(k_r2+i-1)-(mrccshift*dbl_mb(k_t2+i-1))
122       else
123        cit = dbl_mb(k_r2+i-1)
124       endif
125
126       if(lusesamefock_it) then
127         orbspin1 = int_mb(k_spin+p1b-1)-1
128         orbspin2 = int_mb(k_spin+p2b-1)-1
129         orbspin3 = int_mb(k_spin+h3b-1)-1
130         orbspin4 = int_mb(k_spin+h4b-1)-1
131
132         orbindex1 = (1-orbspin1+int_mb(k_mo_indexm(iref)+
133     1 int_mb(k_offset+p1b-1)+p1-1))/2
134         orbindex2 = (1-orbspin2+int_mb(k_mo_indexm(iref)+
135     1 int_mb(k_offset+p2b-1)+p2-1))/2
136         orbindex3 = (1-orbspin3+int_mb(k_mo_indexm(iref)+
137     1 int_mb(k_offset+h3b-1)+h3-1))/2
138         orbindex4 = (1-orbspin4+int_mb(k_mo_indexm(iref)+
139     1 int_mb(k_offset+h4b-1)+h4-1))/2
140
141       blck1 = orbinblck(orbindex1,orbspin1+1,1)
142       blck2 = orbinblck(orbindex2,orbspin2+1,1)
143       blck3 = orbinblck(orbindex3,orbspin3+1,1)
144       blck4 = orbinblck(orbindex4,orbspin4+1,1)
145
146       off1 = offsetinblck(orbindex1,orbspin1+1,1)
147       off2 = offsetinblck(orbindex2,orbspin2+1,1)
148       off3 = offsetinblck(orbindex3,orbspin3+1,1)
149       off4 = offsetinblck(orbindex4,orbspin4+1,1)
150
151      d1 = -dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)+blck1-1)+off1)
152      d2 = -dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)+blck2-1)+off2)
153      d3 = dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)+blck3-1)+off3)
154      d4 = dbl_mb(k_evl_sortedm(1)+int_mb(k_offsetm(1)+blck4-1)+off4)
155
156       else
157        d1 = -dbl_mb(k_evl_sorted+int_mb(k_offset+p1b-1)+p1-1)
158        d2 = -dbl_mb(k_evl_sorted+int_mb(k_offset+p2b-1)+p2-1)
159        d3 = dbl_mb(k_evl_sorted+int_mb(k_offset+h3b-1)+h3-1)
160        d4 = dbl_mb(k_evl_sorted+int_mb(k_offset+h4b-1)+h4-1)
161       endif
162
163       dbl_mb(k_r2+i-1) = cit/(d1+d2+d3+d4+shift+mrccshift)
164
165c       if(abs(dbl_mb(k_r2+i-1)).gt.20.0d0) then
166c           if(dbl_mb(k_r2+i-1).lt.0.0d0) then
167c               dbl_mb(k_r2+i-1)=-20.0d0
168c           else
169c               dbl_mb(k_r2+i-1)=20.0d0
170c           endif
171c           write(6,"('2RESIDUE HAS BEEN CUTED')")
172c       endif
173      if(iter .lt. 4) then
174       if(.not. lreadt) then
175        if(abs(dbl_mb(k_r2+i-1)).gt.0.1d0) then
176         if(dbl_mb(k_r2+i-1).lt.0.0d0) then
177               dbl_mb(k_r2+i-1)=-0.01d0
178         else
179               dbl_mb(k_r2+i-1)=0.01d0
180         endif
181c        if(nodezero)
182c          write(luout,*)"t2 residue has been modified to 0.01"
183c        if(nodezero) call util_flush(LuOut)
184        endif
185       endif
186      endif
187
188                        enddo
189                      enddo
190                    enddo
191                  enddo
192                  call add_hash_block(d_t2,dbl_mb(k_r2),size,
193     1              int_mb(k_t2_offset),((((p1b-noab-1)*nvab+p2b-noab-1)
194     2              *noab+h3b-1)*noab+h4b-1))
195c update of the res.-double vector to the form of increment used in DIIS proc.
196c                  call put_hash_block(d_r2,dbl_mb(k_r2),size,
197c     1              int_mb(k_t2_offset),((((p1b-noab-1)*nvab+p2b-noab-1)
198c     2              *noab+h3b-1)*noab+h4b-1))
199c ----------------------------------------------------------------------------
200                  if (nodezero.and.util_print('t2',print_debug))
201     1              then
202                    call get_hash_block(d_t2,dbl_mb(k_r2),size,
203     1              int_mb(k_t2_offset),((((p1b-noab-1)*nvab+p2b-noab-1)
204     2              *noab+h3b-1)*noab+h4b-1))
205                    call ma_print_compact
206     2              (dbl_mb(k_r2),size,1,'t2')
207                  endif
208c                 if(lsubterm) then
209                  if (.not.ma_pop_stack(l_t2))
210     1              call errquit('tce_jacobi_t2: MA problem',1,MA_ERR)
211c                 endif
212                  if (.not.ma_pop_stack(l_r2))
213     1              call errquit('tce_jacobi_t2: MA problem',1,MA_ERR)
214                endif
215                endif
216                endif
217cc                next = nxtask(nprocs,1)
218      next = NXTASK(nprocs, 1)
219              endif
220              count = count + 1
221            enddo
222          enddo
223        enddo
224      enddo
225cc      next = nxtask(-nprocs,1)
226      next = NXTASK(-nprocs, 1)
227      call ga_sync()
228      call util_flush(LuOut)
229      return
230      end
231c
232c
233c
234
235c $Id$
236