1 SUBROUTINE lambda_ccsd_t(d_t1,k_t1_offset,d_t2,k_t2_offset, 2 1 d_y1,k_y1_offset,d_y2,k_y2_offset, 3 2 d_f1,k_f1_offset,d_v2,k_v2_offset, 4 3 energy1,energy2) 5C 6C $Id$ 7C 8 IMPLICIT NONE 9#include "global.fh" 10#include "mafdecls.fh" 11#include "util.fh" 12#include "errquit.fh" 13#include "tce.fh" 14#include "tce_main.fh" 15 integer d_t1 16 integer k_t1_offset 17 integer d_t2 18 integer k_t2_offset 19 integer d_y1 20 integer k_y1_offset 21 integer d_y2 22 integer k_y2_offset 23 integer d_f1 24 integer k_f1_offset 25 integer d_v2 26 integer k_v2_offset 27 integer t_h1b, t_h1 28 integer t_h2b, t_h2 29 integer t_h3b, t_h3 30 integer t_p4b, t_p4 31 integer t_p5b, t_p5 32 integer t_p6b, t_p6 33 integer k_tdoubles,l_tdoubles 34 integer k_ysingles,l_ysingles 35 integer k_ydoublesA,l_ydoublesA ! L3 form, before sorting 36 integer k_ydoubles,l_ydoubles ! T3 form, after sorting 37 integer size,i 38 integer g_energy 39 integer nxtask 40 integer next 41 integer nprocs 42 integer count 43 double precision energy1,energy2 44 double precision factor 45 external nxtask 46C 47#ifdef DEBUG_PRINT 48 print*,"entering lambda_ccsd_t" 49#endif 50 if (.not.ga_create(mt_dbl,1,1,'perturbative',1,1,g_energy)) 51 1 call errquit('lambda_ccsd_t: GA problem',0,GA_ERR) 52 nprocs = GA_NNODES() 53 count = 0 54 next = nxtask(nprocs,1) 55 energy1=0.0d0 56 energy2=0.0d0 57 do t_p4b = noab+1,noab+nvab 58 do t_p5b = t_p4b,noab+nvab 59 do t_p6b = t_p5b,noab+nvab 60 do t_h1b = 1,noab 61 do t_h2b = t_h1b,noab 62 do t_h3b = t_h2b,noab 63 if (int_mb(k_spin+t_p4b-1) 64 1 +int_mb(k_spin+t_p5b-1) 65 2 +int_mb(k_spin+t_p6b-1) 66 3 .eq.int_mb(k_spin+t_h1b-1) 67 4 +int_mb(k_spin+t_h2b-1) 68 5 +int_mb(k_spin+t_h3b-1)) then 69 if ((.not.restricted).or. 70 1 (int_mb(k_spin+t_p4b-1) 71 1 +int_mb(k_spin+t_p5b-1) 72 2 +int_mb(k_spin+t_p6b-1) 73 3 +int_mb(k_spin+t_h1b-1) 74 4 +int_mb(k_spin+t_h2b-1) 75 5 +int_mb(k_spin+t_h3b-1).le.8)) then 76 if (ieor(int_mb(k_sym+t_p4b-1), 77 1 ieor(int_mb(k_sym+t_p5b-1), 78 2 ieor(int_mb(k_sym+t_p6b-1), 79 3 ieor(int_mb(k_sym+t_h1b-1), 80 4 ieor(int_mb(k_sym+t_h2b-1), 81 5 int_mb(k_sym+t_h3b-1)))))).eq.0) then 82c 83 if (next.eq.count) then 84c 85 size = int_mb(k_range+t_p4b-1) 86 1 * int_mb(k_range+t_p5b-1) 87 2 * int_mb(k_range+t_p6b-1) 88 3 * int_mb(k_range+t_h1b-1) 89 4 * int_mb(k_range+t_h2b-1) 90 5 * int_mb(k_range+t_h3b-1) 91c 92 if (.not.MA_PUSH_GET(mt_dbl,size,'(T) right doubles', 93 1 l_tdoubles,k_tdoubles)) 94 2 call errquit('lambda_ccsd_t: MA error',301,MA_ERR) 95 if (.not.MA_PUSH_GET(mt_dbl,size,'(T) left singles', 96 1 l_ysingles,k_ysingles)) 97 2 call errquit('lambda_ccsd_t: MA error',302,MA_ERR) 98 if (.not.MA_PUSH_GET(mt_dbl,size,'(T) left doubles', 99 1 l_ydoubles,k_ydoubles)) 100 2 call errquit('lambda_ccsd_t: MA error',303,MA_ERR) 101c 102 do i = 1, size 103 dbl_mb(k_tdoubles+i-1) = 0.0d0 104 dbl_mb(k_ysingles+i-1) = 0.0d0 105 dbl_mb(k_ydoubles+i-1) = 0.0d0 106 enddo 107c 108#ifdef DEBUG_PRINT 109 print*,"before ccsd_t_doubles" 110#endif 111c 112c <T|{W*T2}|0> : i0(p4 p5 p6 h1 h2 h3)_vt += -P(9)*Sum(h7)*t(p4 p5 h1 h7)_t*v(h7 p6 h2 h3)_v 113c : i0(p4 p5 p6 h1 h2 h3)_vt += -P(9)*Sum(p7)*t(p4 p7 h1 h2)_t*v(p5 p6 h3 p7)_v 114c 115 call ccsd_t_doubles(dbl_mb(k_tdoubles),d_t2,d_v2, 116 1 k_t2_offset,k_v2_offset, 117 2 t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2) 118c 119#ifdef DEBUG_PRINT 120 print*,"before lambda_ccsd_t_left" 121#endif 122c 123c <0|L(F+W)|T> : i0(h4 h5 h6 p1 p2 p3)_yv += P(9)*y(h4 p1)_y*v(h5 h6 p2 p3)_v 124c 125 call lambda_ccsd_t_left(dbl_mb(k_ysingles),d_f1,d_v2, 126 1 d_y1,d_y2,k_f1_offset,k_v2_offset, 127 2 k_y1_offset,k_y2_offset, 128 3 t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 129c 130c <0|L(F+W)|T> : i0(h4 h5 h6 p1 p2 p3)_yf += P(9)*y(h4 h5 p1 p2)_y*f(h6 p3)_f 131c : i0(h4 h5 h6 p1 p2 p3)_yv += -P(9)*Sum(h7)*y(h4 h7 p1 p2)_y*v(h5 h6 h7 p3)_v 132c : i0(h4 h5 h6 p1 p2 p3)_yv += -P(9)*Sum(p7)*y(h4 h5 p1 p7)_y*v(h6 p7 p2 p3)_v 133c 134 call lambda_ccsd_t_left(dbl_mb(k_ydoubles),d_f1,d_v2, 135 1 d_y1,d_y2,k_f1_offset,k_v2_offset, 136 2 k_y1_offset,k_y2_offset, 137 3 t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2) 138c 139 if (restricted) then 140 factor = 2.0d0 141 else 142 factor = 1.0d0 143 endif 144 if ((t_p4b.eq.t_p5b).and.(t_p5b.eq.t_p6b)) then 145 factor = factor / 6.0d0 146 else if ((t_p4b.eq.t_p5b).or.(t_p5b.eq.t_p6b)) then 147 factor = factor / 2.0d0 148 endif 149 if ((t_h1b.eq.t_h2b).and.(t_h2b.eq.t_h3b)) then 150 factor = factor / 6.0d0 151 else if ((t_h1b.eq.t_h2b).or.(t_h2b.eq.t_h3b)) then 152 factor = factor / 2.0d0 153 endif 154 i = 0 155 do t_p4 = 1, int_mb(k_range+t_p4b-1) 156 do t_p5 = 1, int_mb(k_range+t_p5b-1) 157 do t_p6 = 1, int_mb(k_range+t_p6b-1) 158 do t_h1 = 1, int_mb(k_range+t_h1b-1) 159 do t_h2 = 1, int_mb(k_range+t_h2b-1) 160 do t_h3 = 1, int_mb(k_range+t_h3b-1) 161 i = i + 1 162 energy1 = energy1 + factor * dbl_mb(k_tdoubles+i-1) 163 1 * dbl_mb(k_ydoubles+i-1) 164 2 / (-dbl_mb(k_evl_sorted+int_mb(k_offset+t_p4b-1)+t_p4-1) 165 3 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p5b-1)+t_p5-1) 166 4 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p6b-1)+t_p6-1) 167 5 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h1b-1)+t_h1-1) 168 6 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h2b-1)+t_h2-1) 169 7 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h3b-1)+t_h3-1)) 170 energy2 = energy2 + factor * dbl_mb(k_tdoubles+i-1) 171 1 * (dbl_mb(k_ysingles+i-1) + dbl_mb(k_ydoubles+i-1)) 172 2 / (-dbl_mb(k_evl_sorted+int_mb(k_offset+t_p4b-1)+t_p4-1) 173 3 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p5b-1)+t_p5-1) 174 4 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p6b-1)+t_p6-1) 175 5 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h1b-1)+t_h1-1) 176 6 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h2b-1)+t_h2-1) 177 7 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h3b-1)+t_h3-1)) 178 enddo 179 enddo 180 enddo 181 enddo 182 enddo 183 enddo 184 if (.not.MA_POP_STACK(l_ydoubles)) 185 1 call errquit('lambda_ccsd_t',314,MA_ERR) 186 if (.not.MA_POP_STACK(l_ysingles)) 187 1 call errquit('lambda_ccsd_t',312,MA_ERR) 188 if (.not.MA_POP_STACK(l_tdoubles)) 189 1 call errquit('lambda_ccsd_t',311,MA_ERR) 190c 191 next = nxtask(nprocs,1) 192 endif 193 count = count + 1 194c 195 endif 196 endif 197 endif 198 enddo 199 enddo 200 enddo 201 enddo 202 enddo 203 enddo 204 next = nxtask(-nprocs,1) 205 call ga_zero(g_energy) 206 call ga_acc(g_energy,1,1,1,1,energy1,1,1.0d0) 207 call ga_sync() 208 call ga_get(g_energy,1,1,1,1,energy1,1) 209 call ga_zero(g_energy) 210 call ga_acc(g_energy,1,1,1,1,energy2,1,1.0d0) 211 call ga_sync() 212 call ga_get(g_energy,1,1,1,1,energy2,1) 213 if (.not.ga_destroy(g_energy)) 214 1 call errquit('lambda_ccsd_t: GA problem',1,GA_ERR) 215 return 216 end 217