1 SUBROUTINE cr_eomccsd_t_6dts(d_x1,k_x1_offset,d_x2,k_x2_offset, 2 1 d_t1,k_t1_offset,d_t2,k_t2_offset, 3 2 d_f1,k_f1_offset,d_v2,k_v2_offset, 4 3 d_e,k_e_offset, 5 4 d_ex1,k_ex1_offset,size_ex1, 6 5 d_ex2,k_ex2_offset,size_ex2, 7 6 d_c1,k_c1_offset,size_c1, 8 7 d_c2,k_c2_offset,size_c2, 9 8 excit,energy1,energy2,size_t1,size_x1, 10 8 r0xx,xmem) 11C 12C $Id$ 13C 14c from this point on we assume that corresponding one- and two- 15c body components of R(d_x1,d_x2) and T(d_t1,d_t2) operators, 16c and corresponding excitation energy are available. 17c objects c1 and c2 are created in tce_energy right before 18c calling this procedure. The same structure of x1 and c1 and 19c x2 and c2 is assumed (irrep_c = irrep_x) 20c 21c 1. calculate R0 22c 2. calculate c1 and c2 vectors (c1excit,c2excit) 23c 3. calculate d 24c 25c local representations go only to "triple" procedures 26c see size_t1 size_x1 in the header of this procedure 27c 28 IMPLICIT NONE 29#include "global.fh" 30#include "mafdecls.fh" 31#include "util.fh" 32#include "errquit.fh" 33#include "tce.fh" 34#include "tce_main.fh" 35 integer d_t1 36 integer k_t1_offset 37 integer d_t2 38 integer k_t2_offset 39c --- sizes --- 40 integer size_c1,size_c2 41 integer size_ex1,size_ex2 42 integer size_d2 43c ------------- 44 integer d_x1 45 integer k_x1_offset 46 integer d_x2 47 integer k_x2_offset 48 double precision r0xx ! r0 49 double precision dr0xx ! r0*r0 50 double precision d1xx ! <(singles)+singles> 51 double precision d2xx ! <(doubles)+doubles> 52 double precision d12xx ! <(singles+doubles)+(singles+doubles)> 53 double precision d1xxt,d2xxt 54 double precision d1xxr,d2xxr 55 double precision d1xxtr,d2xxtr 56 double precision excit ! eomsd excitation energy 57 logical lr0 ! (true) r0*M3-calculated 58 integer d_ex1 59 integer k_ex1_offset 60 integer d_ex2 61 integer k_ex2_offset 62 integer d_c1 63 integer k_c1_offset 64 integer d_c2 65 integer k_c2_offset 66 integer d_d2 67 integer l_d2_offset 68 integer k_d2_offset 69c -- 6DTS -- 70 integer dp4,dp5,ii,jj,istart,istop,jstart,jstop 71 integer mdp4,mdp5 72 integer maxp4,maxp5 73 integer slice_dp4,slice_dp5,qp4,qp5 74ccx double precision xlocal 75 integer xlocal 76 integer range_p4,range_p5,range_p6 77 integer range_h1,range_h2,range_h3 78c ---- 79c === jaguar === 80 integer xmem 81c ============== 82 integer d_f1 83 integer k_f1_offset 84 integer d_v2 85 integer k_v2_offset 86 integer d_e 87 integer k_e_offset 88 integer t_h1b, t_h1 89 integer t_h2b, t_h2 90 integer t_h3b, t_h3 91 integer t_p4b, t_p4 92 integer t_p5b, t_p5 93 integer t_p6b, t_p6 94 integer k_den,l_den 95 integer k_right,l_right !r0*M3(T1,T2) 96c --- cr_ccsd_t_E --- 97 integer k_left,l_left 98c ------------------- 99c - T1/X1 LOCALIZATION ------------------- 100 integer l_t1_local,k_t1_local 101 integer l_x1_local,k_x1_local 102 integer size_t1,size_x1 103c --------------------------------------- 104 integer size,i 105 integer g_energy 106 integer nxtask 107 integer next 108 integer nprocs 109 integer count 110c --- new intermediates --- 111c cr_ccsd_t_N1 or cr_ccsd_t_N 112 integer d_i1_1,d_i1_2,d_i1_3 113 integer k_i1_offset_1,k_i1_offset_2,k_i1_offset_3 114 integer l_i1_offset_1,l_i1_offset_2,l_i1_offset_3 115c cr_ccsd_t_E 116 integer d_i1_4,k_i1_offset_4,l_i1_offset_4 117c cr_ccsd_t_N2 118 integer d_i2_1,d_i2_2,d_i2_3,d_i2_4,d_i2_5,d_i2_6 119 integer k_i2_offset_1,k_i2_offset_2,k_i2_offset_3 120 integer k_i2_offset_4,k_i2_offset_5,k_i2_offset_6 121 integer l_i2_offset_1,l_i2_offset_2,l_i2_offset_3 122 integer l_i2_offset_4,l_i2_offset_5,l_i2_offset_6 123c q3rexpt 124 integer d_i3_1,d_i3_2 125 integer k_i3_offset_1,k_i3_offset_2 126 integer l_i3_offset_1,l_i3_offset_2 127c ---------------------------- 128 double precision energy1,energy2 129 double precision factor 130 double precision den1,num1 131 double precision den2,num2 132 double precision denex 133 double precision cpu,wall 134 character*255 filename 135 logical nodezero 136 external nxtask 137c 138c === jaguar ==== 139c xmem max local memory in MB 140c =============== 141c 142c 143c 144c 145 nodezero=(ga_nodeid().eq.0) 146c 147c 148 if(nodezero) then 149 write(6,*)'I AM IN CR_EOMCCSD_T_6DTS' 150 call util_flush(6) 151 end if 152c 153c Getting R0 154c 155c === jaguar === 156 if(.not.read_in3) then ! --- read_in3 --- 157 r0xx=0.0d0 158 call nr0(d_f1,d_e,d_t1,d_v2,d_x1,d_x2,k_f1_offset, 159 &k_e_offset,k_t1_offset,k_v2_offset,k_x1_offset,k_x2_offset) 160 call reconcilefile(d_e,1) 161 call get_block(d_e,r0xx,1,0) 162 if(dabs(excit).gt.1.0d-7) then 163 r0xx = r0xx/excit 164 else 165 write(6,1000) 166 end if 167 end if ! --- read_in3 168 dr0xx = r0xx*r0xx 169 lr0 = .true. 170 if(dabs(r0xx).lt.1.0d-7) lr0 = .false. 171c === jaguar === 172 if(nodezero) write(6,*)'R0: ',r0xx 173 call util_flush(6) 174c ============== 175c 176c Now on ga with handle d_e we store corresponding R0 value 177c 178c 179c 180c Calculating one- and two-body overlaps 181c 182 if(lr0) then !symmetry of the reference 183 call tce_zero(d_ex1,size_ex1) 184 call tce_zero(d_ex2,size_ex2) 185 call c1_c1(d_t1,d_ex1,k_t1_offset,k_ex1_offset) 186 call reconcilefile(d_ex1,1) 187 d1xxt = 0.0d0 188 call get_block(d_ex1,d1xxt,1,0) 189 call t2t12(d_c2,d_t1,d_t2,k_c2_offset,k_t1_offset,k_t2_offset) 190 call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset) 191 call reconcilefile(d_ex2,1) 192 d2xxt = 0.0d0 193 call get_block(d_ex2,d2xxt,1,0) 194c forming vector R1T1+R2, d_d2's irrep corresponds to irrep_x 195c (in this case this is fully symmetric situation) 196c on d_c2 we have (T2+1/2T1*T1)|Phi\rangle 197 irrep_d=irrep_c 198 call tce_zero(d_ex1,size_ex1) 199 call tce_x2_offset(l_d2_offset,k_d2_offset,size_d2) 200 call tce_filename('d2',filename) 201 call createfile(filename,d_d2,size_d2) 202 call tce_zero(d_d2,size_d2) 203 call c2excit2(d_d2,d_t1,d_x1,d_x2,k_d2_offset,k_t1_offset,k_ 204 & x1_offset,k_x2_offset) 205c <C2+ D2> here 206 call c2_d2(d_c2,d_d2,d_ex1,k_c2_offset,k_d2_offset,k_ex1_offset) 207 call reconcilefile(d_ex1,1) 208 d2xxtr = 0.0d0 209 call get_block(d_ex1,d2xxtr,1,0) 210 call deletefile(d_d2) 211 if (.not.ma_pop_stack(l_d2_offset)) 212 1 call errquit("tce_energy: MA problem",36,MA_ERR) 213c 214 call tce_zero(d_ex1,size_ex1) 215 call tce_zero(d_ex2,size_ex2) 216 call x1_t1(d_ex1,d_t1,d_x1,k_ex1_offset,k_t1_offset, 217 & k_x1_offset) 218 call reconcilefile(d_ex1,1) 219 d1xxtr=0.0d0 220 call get_block(d_ex1,d1xxtr,1,0) 221c 222 call tce_zero(d_ex1,size_ex1) 223 call tce_zero(d_ex2,size_ex2) 224 call tce_zero(d_c1,size_c1) 225 call tce_zero(d_c2,size_c2) 226 call c1excit2(d_c1,d_x1,k_c1_offset,k_x1_offset) 227 call c2excit2(d_c2,d_t1,d_x1,d_x2,k_c2_offset,k_t1_offset,k_ 228 & x1_offset,k_x2_offset) 229 call c1_c1(d_c1,d_ex1,k_c1_offset,k_ex1_offset) 230 call reconcilefile(d_ex1,1) 231 d1xxr = 0.0d0 232 call get_block(d_ex1,d1xxr,1,0) 233 call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset) 234 call reconcilefile(d_ex2,1) 235 d2xxr = 0.0d0 236 call get_block(d_ex2,d2xxr,1,0) 237 d12xx = d1xxt*r0xx*r0xx+d2xxt*r0xx*r0xx+d1xxr+d2xxr+ 238 & 2.0d0*d1xxtr*r0xx+2.0d0*d2xxtr*r0xx 239 else !symmetry different form the symmetry of the reference 240c <(R1+R0*T1)^{\dagger} (R1+R0*T1)> 241 call c1excit2(d_c1,d_x1,k_c1_offset,k_x1_offset) 242 call c2excit2(d_c2,d_t1,d_x1,d_x2,k_c2_offset,k_t1_offset,k_ 243 & x1_offset,k_x2_offset) 244 call c1_c1(d_c1,d_ex1,k_c1_offset,k_ex1_offset) 245 call reconcilefile(d_ex1,1) 246 d1xx = 0.0d0 247 call get_block(d_ex1,d1xx,1,0) 248c <(R2+R1T1)^{\dagger} (R2+R1T1)> 249 call c2_c2(d_c2,d_ex2,k_c2_offset,k_ex2_offset) 250 call reconcilefile(d_ex2,1) 251 d2xx = 0.0d0 252 call get_block(d_ex2,d2xx,1,0) 253 d12xx = d1xx+d2xx 254 end if 255c - T1/X1 LOCALIZATION ---------- 256c opening l_t1_local and l_x1_local 257 if (.not.MA_PUSH_GET(mt_dbl,size_t1,'t1_local', 258 1 l_t1_local,k_t1_local)) 259 1 call errquit('t1_local',1,MA_ERR) 260 if (.not.MA_PUSH_GET(mt_dbl,size_x1,'x1_local', 261 1 l_x1_local,k_x1_local)) 262 1 call errquit('x1_local',1,MA_ERR) 263 call ma_zero(dbl_mb(k_t1_local),size_t1) 264 call ma_zero(dbl_mb(k_x1_local),size_x1) 265c copy d_t1 ==> l_t1_local 266c copy x1(ivec) ==> l_x1_local 267 call ga_get(d_t1,1,size_t1,1,1,dbl_mb(k_t1_local),1) 268 call ga_get(d_x1,1,size_x1,1,1,dbl_mb(k_x1_local),1) 269c ------------------------------- 270c 271c Caution! k_right & k_den are not even allocated yet 272c but they will not be used. 273c --------------- initialization ----------------------------------- 274 cpu=-util_cpusec() 275 wall=-util_wallsec() 276 if(lr0) then 277 call cr_ccsd_t_N(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2, 278 1 k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2, 279 2 k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1, 280 3 l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 281 end if 282c 283 call creomsd_t_n2_mem(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2, 284 &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2,k_x1_local, 285 &d_x2,k_f1_offset,k_i2_offset_1, 286 &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset, 287 &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2, 288 &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 289c 290c 291 if(lr0) then 292 call cr_ccsd_t_E(dbl_mb(k_left),d_i1_4, 293 1 k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset, 294 2 l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 295 call q3rexpt2(dbl_mb(k_left),d_i3_1, 296 &k_t1_local,d_t2,k_x1_local,d_x2, 297 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 298 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 299 else 300 call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2, 301 &k_x1_local,d_x2, 302 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 303 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,1) 304 end if 305 cpu=cpu+util_cpusec() 306 wall=wall+util_wallsec() 307 if(nodezero) then 308 write(6,*)' ' 309 write(6,111) cpu,wall 310 if(lr0) then 311 write(6,*)'cr_N i1_1 nr of boxes ',int_mb(k_i1_offset_1) 312 write(6,*)'cr_N i1_2 nr of boxes ',int_mb(k_i1_offset_2) 313 end if 314 write(6,*)'creom_N i2_1 nr of boxes ',int_mb(k_i2_offset_1) 315 write(6,*)'creom_N i2_2 nr of boxes ',int_mb(k_i2_offset_2) 316 write(6,*)'creom_N i2_3 nr of boxes ',int_mb(k_i2_offset_3) 317 write(6,*)'creom_N i2_4 nr of boxes ',int_mb(k_i2_offset_4) 318 write(6,*)' ' 319 call util_flush(6) 320 end if 321 111 format('CR-EOMCCSD(T) Intermediates cpu wall time ',2f15.1) 322c ------------------------------------------------------------------ 323c 324c Get the numerator 325c 326 num1 = 0.0d0 327 den1 = 0.0d0 328c 329 num2 = 0.0d0 330 den2 = 0.0d0 331c 332 cpu=-util_cpusec() 333 wall=-util_wallsec() 334c 335 if (.not.ga_create(mt_dbl,1,1,'perturbative',1,1,g_energy)) 336 1 call errquit('ccsd_t: GA problem',0,GA_ERR) 337 nprocs = GA_NNODES() 338 count = 0 339 next = nxtask(nprocs,1) 340 do t_p4b = noab+1,noab+nvab 341 do t_p5b = t_p4b,noab+nvab 342 do t_p6b = t_p5b,noab+nvab 343 do t_h1b = 1,noab 344 do t_h2b = t_h1b,noab 345 do t_h3b = t_h2b,noab 346 if (int_mb(k_spin+t_p4b-1) 347 1 +int_mb(k_spin+t_p5b-1) 348 2 +int_mb(k_spin+t_p6b-1) 349 3 .eq.int_mb(k_spin+t_h1b-1) 350 4 +int_mb(k_spin+t_h2b-1) 351 5 +int_mb(k_spin+t_h3b-1)) then 352 if ((.not.restricted).or. 353 1 (int_mb(k_spin+t_p4b-1) 354 1 +int_mb(k_spin+t_p5b-1) 355 2 +int_mb(k_spin+t_p6b-1) 356 3 +int_mb(k_spin+t_h1b-1) 357 4 +int_mb(k_spin+t_h2b-1) 358 5 +int_mb(k_spin+t_h3b-1).le.8)) then 359 if (ieor(int_mb(k_sym+t_p4b-1), 360 1 ieor(int_mb(k_sym+t_p5b-1), 361 2 ieor(int_mb(k_sym+t_p6b-1), 362 3 ieor(int_mb(k_sym+t_h1b-1), 363 4 ieor(int_mb(k_sym+t_h2b-1), 364 5 int_mb(k_sym+t_h3b-1)))))).eq.irrep_x) then 365c 366 if (next.eq.count) then 367c 368c Symmetry control (above) 369c -- 6DTS -- 370 range_p4 = int_mb(k_range+t_p4b-1) 371 range_p5 = int_mb(k_range+t_p5b-1) 372 range_p6 = int_mb(k_range+t_p6b-1) 373 range_h1 = int_mb(k_range+t_h1b-1) 374 range_h2 = int_mb(k_range+t_h2b-1) 375 range_h3 = int_mb(k_range+t_h3b-1) 376c 377 dp4=1 378 dp5=1 379 300 continue 380 xlocal=((int_mb(k_range+t_p4b-1)/dp4) 381 1 * (int_mb(k_range+t_p5b-1)/dp5)) 382 2 * int_mb(k_range+t_p6b-1) 383 3 * int_mb(k_range+t_h1b-1) 384 4 * int_mb(k_range+t_h2b-1) 385 5 * int_mb(k_range+t_h3b-1) 386c 387 slice_dp4=range_p4/dp4 388 qp4=range_p4/slice_dp4 389 slice_dp5=range_p5/dp5 390 qp5=range_p5/slice_dp5 391c 392 if((xlocal)/((1024*1024)/8).gt.(xmem)) then 393 if(dp4.lt.int_mb(k_range+t_p4b-1)) then 394 dp4=dp4+1 395 else 396 dp5=dp5+1 397 end if 398 go to 300 399 end if 400c 401 if(slice_dp4*qp4.ne.range_p4) then 402 mdp4=qp4+1 403 else 404 mdp4=qp4 405 end if 406 if(slice_dp5*qp5.ne.range_p5) then 407 mdp5=qp5+1 408 else 409 mdp5=qp5 410 end if 411c 412 size = slice_dp4 413 6 * slice_dp5 414 2 * range_p6 415 3 * range_h1 416 4 * range_h2 417 5 * range_h3 418c 419 maxp4=slice_dp4 420 maxp5=slice_dp5 421c 422 if (.not.MA_PUSH_GET(mt_dbl,size,'moment 2,3', 423 1 l_right,k_right)) call errquit('eomccsd_t',3,MA_ERR) 424 if (.not.MA_PUSH_GET(mt_dbl,size,'denominator', 425 1 l_left,k_left)) call errquit('ccsd_t',3,MA_ERR) 426c 427c 428 do ii=1,mdp4 !do ii --- 429 do jj=1,mdp5 !do jj --- 430 if(ii.ne.mdp4) then 431 istart=(ii-1)*(slice_dp4)+1 432 istop=ii*(slice_dp4) 433 else 434 istart=(ii-1)*(slice_dp4)+1 435 istop=int_mb(k_range+t_p4b-1) 436 end if 437 if(jj.ne.mdp5) then 438 jstart=(jj-1)*(slice_dp5)+1 439 jstop=jj*(slice_dp5) 440 else 441 jstart=(jj-1)*(slice_dp5)+1 442 jstop=int_mb(k_range+t_p5b-1) 443 end if 444c zeroing --- 445 call dfill(size, 0.0d0, dbl_mb(k_right), 1) 446 call dfill(size, 0.0d0, dbl_mb(k_left), 1) 447c ----------- 448c 449c Moments are calculated here 450c 451c 452 if(lr0) then 453 call cr_ccsd_t_N_6dts(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2, 454 1 k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2, 455 2 k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1, 456 3 l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2, 457 4 istart,istop,jstart,jstop,maxp4,maxp5) 458c 459cc do i = 1, size 460cc dbl_mb(k_right+i-1) = r0xx*dbl_mb(k_right+i-1) 461cc enddo 462 call dscal(size,r0xx,dbl_mb(k_right),1) 463c 464 end if 465c 466 call creomsd_t_n2_mem_6dts(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2, 467 &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2, 468 &k_x1_local,d_x2,k_f1_offset,k_i2_offset_1, 469 &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset, 470 &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2, 471 &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2, 472 &istart,istop,jstart,jstop,maxp4,maxp5) 473c 474c 475c Q3(R0+R1+R2)exp(T1+T2)|Ref> calculated here 476c 477 if(lr0) then 478 call cr_ccsd_t_E_6dts(dbl_mb(k_left),d_i1_4, 479 1 k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset, 480 2 l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2, 481 3 istart,istop,jstart,jstop,maxp4,maxp5) 482c 483cc do i = 1, size 484cc dbl_mb(k_left+i-1) = r0xx*dbl_mb(k_left+i-1) 485cc enddo 486 call dscal(size,r0xx,dbl_mb(k_left),1) 487c 488 call q3rexpt2_6dts(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2, 489 &k_x1_local,d_x2, 490 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 491 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2, 492 &istart,istop,jstart,jstop,maxp4,maxp5) 493 else 494 call q3rexpt2_6dts(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2, 495 &k_x1_local,d_x2, 496 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 497 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,2, 498 &istart,istop,jstart,jstop,maxp4,maxp5) 499 end if 500c 501 if (restricted) then 502 factor = 2.0d0 503 else 504 factor = 1.0d0 505 endif 506 if ((t_p4b.eq.t_p5b).and.(t_p5b.eq.t_p6b)) then 507 factor = factor / 6.0d0 508 else if ((t_p4b.eq.t_p5b).or.(t_p5b.eq.t_p6b)) then 509 factor = factor / 2.0d0 510 endif 511 if ((t_h1b.eq.t_h2b).and.(t_h2b.eq.t_h3b)) then 512 factor = factor / 6.0d0 513 else if ((t_h1b.eq.t_h2b).or.(t_h2b.eq.t_h3b)) then 514 factor = factor / 2.0d0 515 endif 516c 517c 518c 519 i = 0 520ccx do t_p4 = 1, int_mb(k_range+t_p4b-1) 521ccx do t_p5 = 1, int_mb(k_range+t_p5b-1) 522 do t_p4 = istart,istop 523 do t_p5 = jstart,jstop 524 do t_p6 = 1, int_mb(k_range+t_p6b-1) 525 do t_h1 = 1, int_mb(k_range+t_h1b-1) 526 do t_h2 = 1, int_mb(k_range+t_h2b-1) 527 do t_h3 = 1, int_mb(k_range+t_h3b-1) 528 i = i + 1 529 denex=-dbl_mb(k_evl_sorted+int_mb(k_offset+t_p4b-1)+t_p4-1) 530 4 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p5b-1)+t_p5-1) 531 5 -dbl_mb(k_evl_sorted+int_mb(k_offset+t_p6b-1)+t_p6-1) 532 6 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h1b-1)+t_h1-1) 533 7 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h2b-1)+t_h2-1) 534 8 +dbl_mb(k_evl_sorted+int_mb(k_offset+t_h3b-1)+t_h3-1) 535 9 +excit 536c numerator 537 num1 = num1 + factor * 538 1 (dbl_mb(k_right+i-1)*dbl_mb(k_right+i-1))/denex 539 num1 = num1 + factor * 540 1 dbl_mb(k_left+i-1)*dbl_mb(k_right+i-1) 541c denominator 542c 543 den1 = den1 + factor * 544 1 (dbl_mb(k_left+i-1)*dbl_mb(k_right+i-1))/denex 545 den1 = den1 + factor * 546 1 (dbl_mb(k_left+i-1)*dbl_mb(k_left+i-1)) 547c 548 enddo 549 enddo 550 enddo 551 enddo 552 enddo 553 enddo 554c 555 enddo !do ii --- 556 enddo !do jj --- 557c 558 if (.not.MA_POP_STACK(l_left)) 559 1 call errquit('eomccsd_t',6,MA_ERR) 560 if (.not.MA_POP_STACK(l_right)) 561 1 call errquit('eomccsd_t',6,MA_ERR) 562c 563 next = nxtask(nprocs,1) 564 endif 565 count = count + 1 566c 567 endif 568 endif 569 endif 570 enddo 571 enddo 572 enddo 573 enddo 574 enddo 575 enddo 576 next = nxtask(-nprocs,1) 577c 578 cpu=cpu+util_cpusec() 579 wall=wall+util_wallsec() 580 if(nodezero) then 581 write(6,*)' ' 582 write(6,112) cpu,wall 583 write(6,*)' ' 584 call util_flush(6) 585 end if 586 112 format('CR-EOMCCSD(T) triples loop cpu wall ',2f15.1) 587c --- toggle = 3 --- 588 if(lr0) then 589 call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2, 590 &k_x1_local,d_x2, 591 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 592 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3) 593 call cr_ccsd_t_E(dbl_mb(k_left),d_i1_4, 594 1 k_t1_local,d_t2,k_i1_offset_4,k_t1_offset,k_t2_offset, 595 2 l_i1_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3) 596 else 597 call q3rexpt2(dbl_mb(k_left),d_i3_1,k_t1_local,d_t2, 598 &k_x1_local,d_x2, 599 &k_i3_offset_1,k_t1_offset,k_t2_offset,k_x1_offset,k_x2_offset, 600 &l_i3_offset_1,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3) 601 end if 602c 603 call creomsd_t_n2_mem(dbl_mb(k_right),d_f1,d_i2_1,d_i2_2, 604 &d_i2_3,d_i2_4,k_t1_local,d_t2,d_v2, 605 &k_x1_local,d_x2,k_f1_offset,k_i2_offset_1, 606 &k_i2_offset_2,k_i2_offset_3,k_i2_offset_4,k_t1_offset,k_t2_offset, 607 &k_v2_offset,k_x1_offset,k_x2_offset,l_i2_offset_1,l_i2_offset_2, 608 &l_i2_offset_3,l_i2_offset_4,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3) 609c 610c 611 if(lr0) then 612 call cr_ccsd_t_N(dbl_mb(k_right),d_f1,d_i1_1,d_i1_2, 613 1 k_t1_local,d_t2,d_v2,k_f1_offset,k_i1_offset_1,k_i1_offset_2, 614 2 k_t1_offset,k_t2_offset,k_v2_offset,l_i1_offset_1, 615 3 l_i1_offset_2,t_h1b,t_h2b,t_h3b,t_p4b,t_p5b,t_p6b,3) 616 end if 617c ------------------ 618c 619 call ga_zero(g_energy) 620 call ga_acc(g_energy,1,1,1,1,num1,1,1.0d0) 621 call ga_sync() 622 call ga_get(g_energy,1,1,1,1,num1,1) 623c 624 call ga_zero(g_energy) 625 call ga_acc(g_energy,1,1,1,1,den1,1,1.0d0) 626 call ga_sync() 627 call ga_get(g_energy,1,1,1,1,den1,1) 628c 629 call ga_zero(g_energy) 630 call ga_acc(g_energy,1,1,1,1,num2,1,1.0d0) 631 call ga_sync() 632 call ga_get(g_energy,1,1,1,1,num2,1) 633c 634 call ga_zero(g_energy) 635 call ga_acc(g_energy,1,1,1,1,den2,1,1.0d0) 636 call ga_sync() 637 call ga_get(g_energy,1,1,1,1,den2,1) 638c 639 if (.not.ga_destroy(g_energy)) 640 1 call errquit('creom_t: GA problem',1,GA_ERR) 641 energy1 = num1/(dr0xx+d12xx+den1) 642 call util_flush(6) 643c - T1/X1 LOCALIZATION ------ 644 if(.not.MA_POP_STACK(l_x1_local)) 645 & call errquit('l_x1_local',4,MA_ERR) 646 if(.not.MA_POP_STACK(l_t1_local)) 647 & call errquit('l_t1_local',4,MA_ERR) 648c --------------------------- 649 1000 format('corresponding excitation energy = 0') 650c 651 return 652 end 653c 654 655