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