1 subroutine uccsdt_makex(g_omega,g_x_big,spini,spinj) 2csh subroutine uccsdt_makex(g_omega,d_x,spini,spinj) 3c 4c x(k,l,i,j) = Sum(u,v) t(u,v,i,j)C(u,k,spini)C(v,l,spinj) 5c where u & v are SO's (no spin), k,b,i,j are MO's with spin labels 6c 7 implicit none 8#include "mafdecls.fh" 9#include "cuccsdtP.fh" 10#include "amplitudes.fh" 11#include "global.fh" 12 integer g_x_big ! [output] DRA handle for X(k,l,i,j) 13 integer g_omega ! [input] GA handle for Omega(u,v,i,j) 14 integer spini, spinj ! [input] Spins of i and j 15 integer spink, spinl ! Spins of k and l 16 integer g_tmp ! Temporary GA x(u,l,i,j) & x(u,v,i,j) 17 integer g_x ! Temporary GA x(k,l,i,j) 18 integer k_buf,l_buf ! Temporary MA for tranformation 19 integer k_tmp,l_tmp ! Temporary MA for tranformation 20 integer k_c,l_c ! Temporary MA for MO coeffs 21 integer lenij,lenul,lenuv,lenu,lenv,lenl 22 integer lenik,lenjl,leni,lenk,leniu 23 integer maxlenuv 24 integer symij,symj,symi,symul,syml 25 integer symuv,symu,symv,symjl,symik,symk 26 integer joff(nbf,0:7) 27 integer loff(nbf,0:7) 28 integer symvoff(0:7,0:7) 29 integer symloff(0:7,0:7) 30 integer ij,jl,ik 31 integer ijlo,ijhi,jllo,jlhi 32 integer i,j,k,l 33 integer llo,lhi,vlo,klo,ilo,ihi,ulo,uhi 34 integer dummy 35csh 36 logical odebug 37 odebug = .false. 38csh 39c 40csh 41 if (odebug) then 42 write(*,*) 'no = ',no(1),no(2) 43 write(*,*) 'nv = ',nv(1),nv(2) 44 write(*,*) 'no_sym = ',(no_sym(i,1),i=0,7),(no_sym(i,2),i=0,7) 45 write(*,*) 'nv_sym = ',(nv_sym(i,1),i=0,7),(nv_sym(i,2),i=0,7) 46 write(*,*) 'o_sym(a) = ',(o_sym(1,i,1),i=0,7),(o_sym(2,i,1),i=0,7) 47 write(*,*) 'o_sym(b) = ',(o_sym(1,i,2),i=0,7),(o_sym(2,i,2),i=0,7) 48 write(*,*) 'v_sym(a) = ',(v_sym(1,i,1),i=0,7),(v_sym(2,i,1),i=0,7) 49 write(*,*) 'v_sym(b) = ',(v_sym(1,i,2),i=0,7),(v_sym(2,i,2),i=0,7) 50 write(*,*) 'bf_per_ir = ',(bf_per_ir(i),i=0,7) 51 endif 52csh 53 spink = spini 54 spinl = spinj 55csh ***** We have to create an empty g_b_big 56 dummy = 0 57 do symjl = 0,7 58 do syml = 0,7 59 symj = ieor(symjl,syml) 60 do l = 1,no_sym(syml,spinl) 61 do j = 1,no_sym(symj,spinj) 62 do symk = 0,7 63 symi = ieor(symjl,symk) 64 do k = 1,no_sym(symk,spink) 65 do i = 1,no_sym(symi,spini) 66 dummy = dummy + 1 67 enddo 68 enddo 69 enddo 70 enddo 71 enddo 72 enddo 73 enddo 74 if (odebug) write(*,*) 'Size of g_x_big = ',dummy 75 if (.not.ga_create(mt_dbl,dummy,1,'x_big', 76 $ -1,1,g_x_big)) call errquit 77 $ ('uccsdt_makex: room for x_big?',dummy) 78csh 79c 80c Addressing & sizes 81c 82 lenij = 0 83 do symij = 0,7 84 do symj = 0,7 85 symi = ieor(symij,symj) 86 do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj) 87 joff(j,symij) = lenij 88 lenij = lenij + no_sym(symi,spini) 89 end do 90 end do 91 end do 92 if (odebug) write(*,*) 'lenij = ',lenij 93c 94c Addressing & sizes 95c 96 do symul = 0,7 97 lenul = 0 98 do syml = 0,7 99 symu = ieor(symul,syml) 100 do l = o_sym(1,syml,spinl),o_sym(2,syml,spinl) 101 loff(l,symul) = lenul 102 lenul = lenul + bf_per_ir(symu) 103 end do 104 end do 105 end do 106c 107c Addressing & sizes 108c 109 do symul = 0,7 110 lenul = 0 111 do syml = 0,7 112 symu = ieor(symul,syml) 113 symloff(syml,symul) = lenul 114 lenul = bf_per_ir(symu)*no_sym(syml,spinl) 115 end do 116 end do 117c 118c Addressing & sizes 119c 120 do symuv = 0,7 121 lenuv = 0 122 do symv = 0,7 123 symu = ieor(symuv,symv) 124 symvoff(symv,symuv) = lenuv 125 lenuv = lenuv + bf_per_ir(symu)*bf_per_ir(symv) 126 end do 127 if (lenuv .gt. maxlenuv) maxlenuv = lenuv 128 end do 129 if (odebug) write(*,*) 'maxlenuv = ',maxlenuv 130c 131c Transformed MO coefficients 132c 133 if (.not. ma_push_get(mt_dbl, nbf*no(spinl), 'c', 134 $ l_c, k_c)) 135 $ call errquit('ma? nbf*no',nbf*no(spinl),0) 136 call ga_get(g_part(spinl),1,nbf,1,no(spinl), 137 $ dbl_mb(k_c),nbf) 138c 139c Allocate temporary GA 140c 141 if (.not.ga_create(mt_dbl,maxlenuv,lenij,'X tmp', 142 $ maxlenuv,-1,g_tmp)) call errquit 143 $ ('uccsdt_makex: room for tmp?',maxlenuv*lenij) 144c 145c Data parallel transform Omega(u,v,i,j)C(v,l) = tmp(u,l,i,j) 146c (tmp was allocated big enough to do this). 147c 148 call ga_distribution(g_omega,ga_nodeid(),ijlo,ijhi,dummy,dummy) 149 ij = 0 150 do symij = 0,7 151 symuv = symij 152 symul = symij 153 do symj = 0,7 154 symi = ieor(symij,symj) 155 do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj) 156 do i = o_sym(1,symi,spini),o_sym(2,symi,spini) 157 ij = ij + 1 158 if (ij.ge.ijlo .and. ij.le.ijhi) then 159 do syml = 0,7 160 symu = ieor(symij,syml) 161 symv = syml 162 lenl = no_sym(syml,spinl) 163 lenu = bf_per_ir(symu) 164 lenv = bf_per_ir(symv) 165 lenul = lenu * lenl 166 lenuv = lenu * lenv 167 llo = o_sym(1,syml,spinl) 168 vlo = bf_per_ir_cum(symv) + 1 169csh 170 if (odebug) write(*,*) 'symij,symj,j,i,syml = ', 171 $ symij,symj,j,i,syml 172csh 173 if (lenuv.gt.0 .and. lenl.gt.0) then 174csh 175 if (odebug) write(*,*) 'lenul,lenuv = ',lenul,lenuv 176csh 177 if (.not. ma_push_get(mt_dbl, lenuv, 'buf', 178 $ l_buf, k_buf)) 179 $ call errquit('ma? lenuv',lenuv,0) 180 if (.not. ma_push_get(mt_dbl, lenul, 'tmp', 181 $ l_tmp, k_tmp)) 182 $ call errquit('ma? lenul',lenul,0) 183c 184c (uv) -> (ul) for l (occupied) & v (SO) of same irrep 185c 186 call ga_get(g_omega,1,lenuv,ij,ij, 187 $ dbl_mb(k_buf),lenu) 188csh 189 if (odebug) call ma_print(dbl_mb(k_buf),lenu,lenv,'UV') 190csh 191 call dgemm('n','n',lenu,lenl,lenv, 192 $ 1.0d0,dbl_mb(k_buf),lenu, 193 $ dbl_mb(k_c),nbf, 194 $ 0.0d0,dbl_mb(k_tmp),lenu) 195csh 196 if (odebug) call ma_print(dbl_mb(k_tmp),lenu,lenl,'UL') 197 if (odebug) write(*,*) 'ga_put in the range ... ', 198 $ symloff(syml,symul)+1,'~',symloff(syml,symul)+lenul, 199 $ ij,'~',ij 200csh 201 call ga_put(g_tmp, 202 $ symloff(syml,symul)+1,symloff(syml,symul)+lenul, 203 $ ij,ij,dbl_mb(k_tmp),lenu) 204 if (.not. ma_pop_stack(l_tmp)) 205 $ call errquit('ma_pop?',l_tmp,0) 206 if (.not. ma_pop_stack(l_buf)) 207 $ call errquit('ma_pop?',l_buf,0) 208 end if 209 end do 210 end if 211 end do 212 end do 213 end do 214 end do 215csh 216 if (odebug) write(*,*) '**************************' 217 if (odebug) write(*,*) 'transformation v -> l done' 218 if (odebug) write(*,*) '**************************' 219csh 220c 221c Perform quarter transformation from tmp(u,l,i,j) to t(k,l,i,j) 222c t(i,k,symk,j,l,syml,symjl) 223c 224 if (.not. ma_pop_stack(l_c)) 225 $ call errquit('ma_pop?',l_c,0) 226 if (.not. ma_push_get(mt_dbl, nbf*no(spink), 'c', 227 $ l_c, k_c)) 228 $ call errquit('ma? nbf*no',nbf*no(spink),0) 229 call ga_get(g_part(spink),1,nbf,1,no(spink), 230 $ dbl_mb(k_c),nbf) 231 do symjl = 0,7 232 symik = symjl 233 do syml = 0,7 234 symj = ieor(symjl,syml) 235c 236c Allocate X block 237c 238 llo = o_sym(1,syml,spinl) 239 lhi = o_sym(2,syml,spinl) 240 lenik = 0 241 do symk = 0,7 242 symi = ieor(symik,symk) 243 lenik = lenik + no_sym(symi,spini) * no_sym(symk,spink) 244 end do 245 lenjl = (lhi-llo+1)*no_sym(symj,spinj) 246 if ((lenik.gt.0).and.(lenjl.gt.0)) then 247csh 248 if (odebug) write(*,*) 'lenik,lenjl = ',lenik,lenjl 249csh 250 if (.not.ga_create(mt_dbl,lenik,lenjl,'l', 251 $ lenik,-1,g_x)) call errquit 252 $ ('uccsdt_makex: room for x?',lenik*lenjl) 253c 254c Data-parallel transformation of k 255c 256 call ga_distribution(g_x,ga_nodeid(),jllo,jlhi,dummy,dummy) 257 jl = 0 258 do l = llo, lhi 259 do j = o_sym(1,symj,spinj),o_sym(2,symj,spinj) 260 jl = jl + 1 261 if (jl.ge.jllo .and. jl.le.jlhi) then 262 ik = 0 263c Compare the next line with the comments in Amplitude.F 264c I think the comments are not accurate 265 do symk = 0, 7 266 symi = ieor(symik,symk) 267 symij = ieor(symi,symj) 268 symu = symk 269 symul = ieor(symu,syml) 270 lenk = no_sym(symk,spink) 271 leni = no_sym(symi,spini) 272 lenu = bf_per_ir(symu) 273 lenik = leni*lenk 274 leniu = leni*lenu 275 klo = o_sym(1,symk,spink) 276 ilo = o_sym(1,symi,spini) 277 ihi = o_sym(2,symi,spini) 278 ulo = bf_per_ir_cum(symu) + 1 279 uhi = ulo + lenu - 1 280 if (lenu.gt.0) then 281 if (.not. ma_push_get(mt_dbl, leniu, 'buf', 282 $ l_buf, k_buf)) 283 $ call errquit('ma? leniu',leniu,0) 284 if (.not. ma_push_get(mt_dbl, lenik, 'tmp', 285 $ l_tmp, k_tmp)) 286 $ call errquit('ma? lenik',lenik,0) 287c 288c (iu) -> (ik) for k (occupied) & u (SO) of same irrep 289c 290 ij = joff(j,symij) 291 do i = ilo,ihi 292 ij = ij + 1 293csh 294 if (odebug) write(*,*) 'i = ',i 295 if (odebug) write(*,*) 'get from range ... ', 296 $ loff(l,symul)+1,'~',loff(l,symul)+lenu, 297 $ ij,'~',ij, 298 $ ' to k_buf+',lenu*(i-ilo) 299csh 300 call ga_get(g_tmp, 301 $ loff(l,symul)+1,loff(l,symul)+lenu, 302 $ ij,ij,dbl_mb(k_buf+lenu*(i-ilo)),lenu) 303 end do 304csh 305 if (odebug) call ma_print(dbl_mb(k_buf),lenu,leni,'X(UI)') 306csh 307 call dgemm('t','n',leni,lenk,lenu, 308 $ 1.0d0,dbl_mb(k_buf),lenu, 309 $ dbl_mb(k_c),nbf, 310 $ 0.0d0,dbl_mb(k_tmp),leni) 311csh 312 if (odebug) call ma_print(dbl_mb(k_tmp),leni,lenk,'X(IK)') 313 if (odebug) write(*,*) 'writing to range ... ', 314 $ ik+1,'~',ik+lenik,jl,'~',jl 315csh 316 call ga_put(g_x,ik+1,ik+lenik,jl,jl, 317 $ dbl_mb(k_tmp),leni) 318 ik = ik + lenik 319 if (.not. ma_pop_stack(l_tmp)) 320 $ call errquit('ma_pop?',l_tmp,0) 321 if (.not. ma_pop_stack(l_buf)) 322 $ call errquit('ma_pop?',l_buf,0) 323 end if 324 end do 325 end if 326 end do 327 end do 328c 329c Write one symmetry block of X to disk 330c 331 if (spini.eq.spinj) call ga_scale(g_x,2.0d0) 332 if (odebug) call ga_print(g_x) 333 call sh_put_x(g_x,syml,symjl,spini,spinj,g_x_big) 334csh if (.not. uccsdt_ampfile_write_t2(d_x, 335csh $ spini, spink, spinj, spinl, symik, 336csh $ llo, lhi, g_x)) 337csh $ call errquit('write_t2 failed',0) 338 if (.not. ga_destroy(g_x)) call errquit('GA 999',0,0) 339 end if 340 end do 341 end do 342 if (.not. ga_destroy(g_tmp)) call errquit('GA 999',0,0) 343 if (.not. ma_pop_stack(l_c)) 344 $ call errquit('ma_pop?',l_c,0) 345c 346 return 347 end 348c $Id$ 349