1C 2C /* Deck so_tmltr */ 3 SUBROUTINE SO_TMLTR(T2AM,SCAL,ISYOPE) 4C 5C This routine is part of the atomic integral direct SOPPA program. 6C 7C Andrea Ligabue, January 2004 8C 9C Developed starting frm CCSD_TCMEPKX 10C Henrik Koch and Alfredo Sanchez. Dec 1994 11C Made workable for non-symmetric T2AM, Keld Bak, Dec 1996 12C 13C Purpose: calculate the left T_aibj starting from the right ones 14C (I need to call that subroutine with SCAL = HALF) 15C 16#include "implicit.h" 17#include "priunit.h" 18C 19 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0, THREE = 3.0D0) 20 PARAMETER (FOUR = 4.0D0) 21C 22 DIMENSION T2AM(*) 23C 24#include "ccorb.h" 25#include "ccsdinp.h" 26#include "ccsdsym.h" 27C 28 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 29C 30C------------------ 31C Add to trace. 32C------------------ 33C 34 CALL QENTER('SO_TMLTR') 35C 36 FAC = TWO/THREE 37C 38 DO 100 ISYMIJ = 1,NSYM 39C 40 ISYMAB = MULD2H(ISYMIJ,ISYOPE) 41C 42 DO 110 ISYMJ = 1,NSYM 43C 44 ISYMI = MULD2H(ISYMJ,ISYMIJ) 45C 46 IF (ISYMI .GT. ISYMJ) GOTO 110 47C 48 DO 120 ISYMB = 1,NSYM 49C 50 ISYMA = MULD2H(ISYMB,ISYMAB) 51C 52 IF (ISYMA .GT. ISYMB) GOTO 120 53C 54 ISYMAI = MULD2H(ISYMA,ISYMI) 55 ISYMBJ = MULD2H(ISYMB,ISYMJ) 56 ISYMBI = MULD2H(ISYMB,ISYMI) 57 ISYMAJ = MULD2H(ISYMA,ISYMJ) 58C 59 DO 130 J = 1,NRHF(ISYMJ) 60C 61 IF (ISYMI .EQ. ISYMJ) THEN 62 NRHFI = J 63 ELSE 64 NRHFI = NRHF(ISYMI) 65 ENDIF 66C 67 IF ( ISYMAI .EQ. ISYMBJ ) THEN 68C 69 DO 140 I = 1,NRHFI 70C 71 DO 150 B = 1,NVIR(ISYMB) 72C 73 IF (ISYMB .EQ. ISYMA) THEN 74 NVIRA = B 75 ELSE 76 NVIRA = NVIR(ISYMA) 77 ENDIF 78C 79 NBI = IT1AM(ISYMB,ISYMI) 80 * + NVIR(ISYMB)*(I - 1) + B 81 NBJ = IT1AM(ISYMB,ISYMJ) 82 * + NVIR(ISYMB)*(J - 1) + B 83C 84 DO 160 A = 1,NVIRA 85C 86 NAI = IT1AM(ISYMA,ISYMI) 87 * + NVIR(ISYMA)*(I - 1) + A 88 NAJ = IT1AM(ISYMA,ISYMJ) 89 * + NVIR(ISYMA)*(J - 1) + A 90C 91 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 92 * + INDEX(NAI,NBJ) 93C 94 NAJBI = IT2AM(ISYMAJ,ISYMBI) 95 * + INDEX(NAJ,NBI) 96C 97 XAIBJ = FAC*(TWO*T2AM(NAIBJ)+T2AM(NAJBI)) 98 XAJBI = FAC*(TWO*T2AM(NAJBI)+T2AM(NAIBJ)) 99C 100 T2AM(NAIBJ) = XAIBJ 101 T2AM(NAJBI) = XAJBI 102C 103 160 CONTINUE 104 150 CONTINUE 105 140 CONTINUE 106C 107 ELSE IF ((ISYMAI.LT.ISYMBJ).AND.(ISYMAJ.LT.ISYMBI)) THEN 108C 109 DO 240 I = 1,NRHFI 110C 111 DO 250 B = 1,NVIR(ISYMB) 112C 113 IF (ISYMB .EQ. ISYMA) THEN 114 NVIRA = B 115 ELSE 116 NVIRA = NVIR(ISYMA) 117 ENDIF 118C 119 NBI = IT1AM(ISYMB,ISYMI) 120 * + NVIR(ISYMB)*(I - 1) + B 121 NBJ = IT1AM(ISYMB,ISYMJ) 122 * + NVIR(ISYMB)*(J - 1) + B 123C 124 DO 260 A = 1,NVIRA 125C 126 NAI = IT1AM(ISYMA,ISYMI) 127 * + NVIR(ISYMA)*(I - 1) + A 128 NAJ = IT1AM(ISYMA,ISYMJ) 129 * + NVIR(ISYMA)*(J - 1) + A 130C 131 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 132 * + NT1AM(ISYMAI) * (NBJ - 1) + NAI 133C 134 NAJBI = IT2AM(ISYMAJ,ISYMBI) 135 * + NT1AM(ISYMAJ) * (NBI - 1) + NAJ 136C 137 XAIBJ = FAC*(TWO*T2AM(NAIBJ)+T2AM(NAJBI)) 138 XAJBI = FAC*(TWO*T2AM(NAJBI)+T2AM(NAIBJ)) 139C 140 T2AM(NAIBJ) = XAIBJ 141 T2AM(NAJBI) = XAJBI 142C 143 260 CONTINUE 144 250 CONTINUE 145 240 CONTINUE 146C 147 ELSE IF ((ISYMAI.GT.ISYMBJ).AND.(ISYMAJ.GT.ISYMBI)) THEN 148C 149 DO 340 I = 1,NRHFI 150C 151 DO 350 B = 1,NVIR(ISYMB) 152C 153 IF (ISYMB .EQ. ISYMA) THEN 154 NVIRA = B 155 ELSE 156 NVIRA = NVIR(ISYMA) 157 ENDIF 158C 159 NBI = IT1AM(ISYMB,ISYMI) 160 * + NVIR(ISYMB)*(I - 1) + B 161 NBJ = IT1AM(ISYMB,ISYMJ) 162 * + NVIR(ISYMB)*(J - 1) + B 163C 164 DO 360 A = 1,NVIRA 165C 166 NAI = IT1AM(ISYMA,ISYMI) 167 * + NVIR(ISYMA)*(I - 1) + A 168 NAJ = IT1AM(ISYMA,ISYMJ) 169 * + NVIR(ISYMA)*(J - 1) + A 170C 171 NAIBJ = IT2AM(ISYMBJ,ISYMAI) 172 * + NT1AM(ISYMBJ) * (NAI - 1) + NBJ 173C 174 NAJBI = IT2AM(ISYMBI,ISYMAJ) 175 * + NT1AM(ISYMBI) * (NAJ - 1) + NBI 176C 177 XAIBJ = FAC*(TWO*T2AM(NAIBJ) + T2AM(NAJBI)) 178 XAJBI = FAC*(TWO*T2AM(NAJBI) + T2AM(NAIBJ)) 179C 180 T2AM(NAIBJ) = XAIBJ 181 T2AM(NAJBI) = XAJBI 182C 183 360 CONTINUE 184 350 CONTINUE 185 340 CONTINUE 186C 187 ELSE IF ((ISYMAI.LT.ISYMBJ).AND.(ISYMAJ.GT.ISYMBI)) THEN 188C 189 DO 440 I = 1,NRHFI 190C 191 DO 450 B = 1,NVIR(ISYMB) 192C 193 IF (ISYMB .EQ. ISYMA) THEN 194 NVIRA = B 195 ELSE 196 NVIRA = NVIR(ISYMA) 197 ENDIF 198C 199 NBI = IT1AM(ISYMB,ISYMI) 200 * + NVIR(ISYMB)*(I - 1) + B 201 NBJ = IT1AM(ISYMB,ISYMJ) 202 * + NVIR(ISYMB)*(J - 1) + B 203C 204 DO 460 A = 1,NVIRA 205C 206 NAI = IT1AM(ISYMA,ISYMI) 207 * + NVIR(ISYMA)*(I - 1) + A 208 NAJ = IT1AM(ISYMA,ISYMJ) 209 * + NVIR(ISYMA)*(J - 1) + A 210C 211 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 212 * + NT1AM(ISYMAI) * (NBJ - 1) + NAI 213C 214 NAJBI = IT2AM(ISYMBI,ISYMAJ) 215 * + NT1AM(ISYMBI) * (NAJ - 1) + NBI 216C 217 XAIBJ = FAC*(TWO*T2AM(NAIBJ) + T2AM(NAJBI)) 218 XAJBI = FAC*(TWO*T2AM(NAJBI) + T2AM(NAIBJ)) 219C 220 T2AM(NAIBJ) = XAIBJ 221 T2AM(NAJBI) = XAJBI 222C 223 460 CONTINUE 224 450 CONTINUE 225 440 CONTINUE 226C 227 ELSE IF ((ISYMAI.GT.ISYMBJ).AND.(ISYMAJ.LT.ISYMBI)) THEN 228C 229 DO 540 I = 1,NRHFI 230C 231 DO 550 B = 1,NVIR(ISYMB) 232C 233 IF (ISYMB .EQ. ISYMA) THEN 234 NVIRA = B 235 ELSE 236 NVIRA = NVIR(ISYMA) 237 ENDIF 238C 239 NBI = IT1AM(ISYMB,ISYMI) 240 * + NVIR(ISYMB)*(I - 1) + B 241 NBJ = IT1AM(ISYMB,ISYMJ) 242 * + NVIR(ISYMB)*(J - 1) + B 243C 244 DO 560 A = 1,NVIRA 245C 246 NAI = IT1AM(ISYMA,ISYMI) 247 * + NVIR(ISYMA)*(I - 1) + A 248 NAJ = IT1AM(ISYMA,ISYMJ) 249 * + NVIR(ISYMA)*(J - 1) + A 250C 251 NAIBJ = IT2AM(ISYMBJ,ISYMAI) 252 * + NT1AM(ISYMBJ) * (NAI - 1) + NBJ 253C 254 NAJBI = IT2AM(ISYMAJ,ISYMBI) 255 * + NT1AM(ISYMAJ) * (NBI - 1) + NAJ 256C 257 XAIBJ = FAC*(TWO*T2AM(NAIBJ) + T2AM(NAJBI)) 258 XAJBI = FAC*(TWO*T2AM(NAJBI) + T2AM(NAIBJ)) 259C 260 T2AM(NAIBJ) = XAIBJ 261 T2AM(NAJBI) = XAJBI 262C 263 560 CONTINUE 264 550 CONTINUE 265 540 CONTINUE 266C 267 END IF 268C 269 130 CONTINUE 270 120 CONTINUE 271 110 CONTINUE 272 100 CONTINUE 273C 274C--------------------------------------- 275C Scale diagonal elements of result. 276C--------------------------------------- 277C 278 IF ((ISYOPE .NE. 1).OR.(SCAL.EQ.1.0D0)) GOTO 1000 279C 280 DO 600 ISYMAI = 1,NSYM 281 DO 610 NAI = 1,NT1AM(ISYMAI) 282 NAIAI = IT2AM(ISYMAI,ISYMAI) + INDEX(NAI,NAI) 283 T2AM(NAIAI) = SCAL*T2AM(NAIAI) 284 610 CONTINUE 285 600 CONTINUE 286C 287C----------------------- 288C Remove from trace. 289C----------------------- 290C 291 1000 CALL QEXIT('SO_TMLTR') 292C 293 RETURN 294 END 295C 296 297 SUBROUTINE SO_TRANTRIP(X2SQ,X2PACK,ISYMTR) 298 299 use so_info, only: sop_dp 300 implicit none 301 302#include "ccorb.h" 303#include "ccsdsym.h" 304#include "soppinf.h" 305 306 real(sop_dp),intent(out) :: X2SQ(*) 307 real(sop_dp),intent(in) :: x2pack(*) 308 integer, intent(in) :: isymtr 309 310 character(len=*), parameter :: myname = 'SO_TRANTRIP' 311 integer :: isymbj, isymai 312 integer :: KOUT 313 314 CALL QENTER(myname) 315 316 KOUT = 1 317 do isymbj = 1, nsym 318 isymai = muld2h(isymbj,isymtr) 319C 320C Calculate the intermediate 321C ~ _ 322C x(aibj) = - /2 x1(abij) - x2(abij) - x3(abij) 323C 324C Note that x2 changes sign when permuting i and j, 325C x3 when permuting a and b, and x1 change sign on both 326C these permutations. 327C 328 CALL SQUAREXT(X2SQ(KOUT),X2PACK,ISYMBJ,ISYMTR) 329 KOUT = KOUT + NT1AM(ISYMAI)*NT1AM(ISYMBJ) 330 331 end do 332 333 IF ( KOUT-1.NE. NT2SQ(ISYMTR) ) THEN 334 PRINT *, 'WARNING', KOUT, NT2SQ(ISYMTR) 335 END IF 336 337 CALL QEXIT(myname) 338 339 contains 340 341 SUBROUTINE SQUAREXT(X2SQ,X2PACK,ISYMR,ISYMT) 342 real(sop_dp),intent(out) :: X2SQ(*) 343 real(sop_dp),intent(in) :: x2pack(*) 344 integer, intent(in) :: isymr, isymt 345 346 real(sop_dp),parameter :: factd = -SQRT(2.0D0), 347 & fact1 = -SQRT(2.0D0), 348 & fact23 = -1.D0, 349 & zero = 0.0d0, 350 & one = 1.0D0, 351 & onem = -1.0D0 352 353 integer :: isymbj, isymai, isymaj, isymbi, isymij, isymab, 354 & isyma, isymb, isymi, isymj 355 integer :: sab1, sab2, nvira, nrhfi, nvirb, nrhfj 356 integer :: nai, naj, nbi, nbj, naibj, najbi, nbiaj, nbjai, 357 & nbibj, nbjbi, najbj, nbjaj, nbjbj, 358 & nabij1, nabij2, nabij3, nbbij2, nabjj3 359 integer :: i, j, b, a 360 integer :: ioff1, ioffj1, ioffbj1 361 integer :: ioff2, ioffj2, ioffbj2 362 integer :: ioff3, ioffj3, ioffbj3 363 real(sop_dp) :: f, fij, fab 364 365 isymbj = isymr 366 isymai = muld2h(isymr,isymt) 367C 368C Calculate the intermediate 369C ~ _ 370C x(aibj) = - /2 x1(abij) - x2(abij) - x3(abij) 371C 372C Note that x2 changes sign when permuting i and j, 373C x3 when permuting a and b, and x1 change sign on both 374C these permutations. 375C 376 377 if (isymt .eq. 1 ) then 378 do isymj = 1, nsym 379 isymb = muld2h(isymj,isymbj) 380 381 ! Ensure isymi <= isymj 382 do isymi = 1, isymj 383 isyma = muld2h(isymi,isymai) 384 ! and isuma <= isymb 385C if (isyma.gt.isymb) cycle 386C 387 isymaj = muld2h(isyma,isymj) 388 isymbi = muld2h(isymb,isymi) 389 isymab = muld2h(isyma,isymb) 390 isymij = muld2h(isymi,isymj) 391C 392 sab1 = ntvv(isymab) 393 sab2 = nsvv(isymab) 394 ioff1 = it2amt1(isymij,isymab) + 395 & sab1*itoo(isymj,isymi) + itvv(isymb,isyma) 396 ioff2 = it2amt2(isymij,isymab) + 397 & sab2*itoo(isymj,isymi) + isvv(isymb,isyma) 398 ioff3 = it2amt3(isymij,isymab) + 399 & sab1*isoo(isymj,isymi) + itvv(isymb,isyma) 400C 401 if (isymi.eq.isymj) then ! isyma == isymb 402 do j = 1, nrhf(isymj) 403 nrhfi = j - 1 404 ioffj1 = ioff1 + ((j-1)*(j-2)/2)*sab1 405 ioffj2 = ioff2 + ((j-1)*(j-2)/2)*sab2 406 ioffj3 = ioff3 + (j*(j-1)/2)*sab1 407 do b = 1, nvir(isymb) 408 nbj = pair_pos(isymb,b,isymj,j) 409 nvira = b - 1 410 ioffbj1 = ioffj1 + (b-1)*(b-2)/2 411 ioffbj2 = ioffj2 + b*(b-1)/2 412 ioffbj3 = ioffj3 + (b-1)*(b-2)/2 413 414 do i = 1, nrhfi 415 nbi = pair_pos(isymb,b,isymi,i) 416 do a = 1, nvira 417 nai = pair_pos(isyma,a,isymi,i) 418 naj = pair_pos(isyma,a,isymj,j) 419 nabij1 = ioffbj1+(i-1)*sab1+a 420 nabij2 = ioffbj2+(i-1)*sab2+a 421 nabij3 = ioffbj3+(i-1)*sab1+a 422 naibj = quad_pos(isymai,nai,isymbj,nbj) 423 najbi = quad_pos(isymaj,naj,isymbi,nbi) 424 nbiaj = quad_pos(isymbi,nbi,isymaj,naj) 425 nbjai = quad_pos(isymbj,nbj,isymai,nai) 426 x2sq(naibj) = fact1*x2pack(nabij1) 427 & + fact23*( x2pack(nabij2) + x2pack(nabij3) ) 428 x2sq(najbi) = -fact1*x2pack(nabij1) 429 & - fact23*( x2pack(nabij2) - x2pack(nabij3) ) 430 x2sq(nbiaj) = -fact1*x2pack(nabij1) 431 & + fact23*( x2pack(nabij2) - x2pack(nabij3) ) 432 x2sq(nbjai) = fact1*x2pack(nabij1) 433 & - fact23*( x2pack(nabij2) + x2pack(nabij3) ) 434 end do ! a 435C Handle a == b, i<j 436 nbbij2 = ioffbj2+(i-1)*sab2+b 437 nbjbi = quad_pos(isymbj,nbj,isymbi,nbi) 438 nbibj = quad_pos(isymbi,nbi,isymbj,nbj) 439 x2sq(nbibj) = factd*x2pack(nbbij2) 440 x2sq(nbjbi) = - factd*x2pack(nbbij2) 441 end do ! i 442C Handle i = j, a < b 443 do a = 1, nvira 444 naj = pair_pos(isyma,a,isymj,j) 445 nabjj3 = ioffbj3+(j-1)*sab1+a 446 najbj = quad_pos(isymaj,naj,isymbj,nbj) 447 nbjaj = quad_pos(isymbj,nbj,isymaj,naj) 448 x2sq(najbj) = factd*x2pack(nabjj3) 449 x2sq(nbjaj) = -factd*x2pack(nabjj3) 450 end do 451 ! i == j, a == b 452C if (isyma.eq.isymb) then 453 nbjbj = quad_pos(isymbj,nbj,isymbj,nbj) 454 x2sq(nbjbj) = zero 455C end if 456 457 end do ! b 458 end do ! j 459 else if (isyma .lt. isymb ) then 460 do j = 1, nrhf(isymj) 461 nrhfi = nrhf(isymi) 462 ioffj1 = ioff1 + (j-1)*nrhfi*sab1 463 ioffj2 = ioff2 + (j-1)*nrhfi*sab2 464 ioffj3 = ioff3 + (j-1)*nrhfi*sab1 465 do b = 1, nvir(isymb) 466 nbj = pair_pos(isymb,b,isymj,j) 467 nvira = nvir(isyma) 468 ioffbj1 = ioffj1 + (b-1)*nvira 469 ioffbj2 = ioffj2 + (b-1)*nvira 470 ioffbj3 = ioffj3 + (b-1)*nvira 471 do i = 1, nrhfi 472 nbi = pair_pos(isymb,b,isymi,i) 473 do a = 1, nvira 474 nai = pair_pos(isyma,a,isymi,i) 475 naj = pair_pos(isyma,a,isymj,j) 476 nabij1 = ioffbj1+(i-1)*sab1+a 477 nabij2 = ioffbj2+(i-1)*sab2+a 478 nabij3 = ioffbj3+(i-1)*sab1+a 479 naibj = quad_pos(isymai,nai,isymbj,nbj) 480 nbjai = quad_pos(isymbj,nbj,isymai,nai) 481 x2sq(naibj) = fact1*x2pack(nabij1) 482 & + fact23*( x2pack(nabij2) + x2pack(nabij3) ) 483 x2sq(nbjai) = fact1*x2pack(nabij1) 484 & - fact23*( x2pack(nabij2) + x2pack(nabij3) ) 485 end do ! a 486 end do ! i 487 end do ! b 488 end do ! j 489 else ! isyma > isymb -> a > b ! 490 ioff1 = it2amt1(isymij,isymab) + 491 & sab1*itoo(isymj,isymi) + itvv(isyma,isymb) 492 ioff2 = it2amt2(isymij,isymab) + 493 & sab2*itoo(isymj,isymi) + isvv(isyma,isymb) 494 ioff3 = it2amt3(isymij,isymab) + 495 & sab1*isoo(isymj,isymi) + itvv(isyma,isymb) 496 nvirb = nvir(isymb) 497 do j = 1, nrhf(isymj) 498 nrhfi = nrhf(isymi) 499 ioffj1 = ioff1 + (j-1)*nrhfi*sab1 500 ioffj2 = ioff2 + (j-1)*nrhfi*sab2 501 ioffj3 = ioff3 + (j-1)*nrhfi*sab1 502 do b = 1, nvir(isymb) 503 nbj = pair_pos(isymb,b,isymj,j) 504 nvira = nvir(isyma) 505 ioffbj1 = ioffj1 + b 506 ioffbj2 = ioffj2 + b 507 ioffbj3 = ioffj3 + b 508 do i = 1, nrhfi 509 nbi = pair_pos(isymb,b,isymi,i) 510 do a = 1, nvira 511 nai = pair_pos(isyma,a,isymi,i) 512 naj = pair_pos(isyma,a,isymj,j) 513 nabij1 = ioffbj1+(i-1)*sab1+(a-1)*nvirb 514 nabij2 = ioffbj2+(i-1)*sab2+(a-1)*nvirb 515 nabij3 = ioffbj3+(i-1)*sab1+(a-1)*nvirb 516 naibj = quad_pos(isymai,nai,isymbj,nbj) 517 nbjai = quad_pos(isymbj,nbj,isymai,nai) 518 x2sq(naibj) = -fact1*x2pack(nabij1) 519 & + fact23*( x2pack(nabij2) - x2pack(nabij3) ) 520 x2sq(nbjai) = -fact1*x2pack(nabij1) 521 & - fact23*( x2pack(nabij2) - x2pack(nabij3) ) 522 end do ! a1 523 end do ! i 524 end do ! b 525 end do ! j 526 end if 527 528 end do ! loop isymi 529 end do ! loop isymj 530C 531 else ! isymai != isymbj 532C 533 do isymj = 1, nsym 534 isymb = muld2h(isymj,isymbj) 535 536 do isymi = 1, nsym 537 isyma = muld2h(isymi,isymai) 538 ! and isuma <= isymb 539C if (isyma.gt.isymb) cycle 540C 541 isymaj = muld2h(isyma,isymj) 542 isymbi = muld2h(isymb,isymi) 543 isymab = muld2h(isyma,isymb) 544 isymij = muld2h(isymi,isymj) 545C 546 sab1 = ntvv(isymab) 547 sab2 = nsvv(isymab) 548 ioff1 = it2amt1(isymij,isymab) + 549 & sab1*itoo(isymj,isymi) + itvv(isymb,isyma) 550 ioff2 = it2amt2(isymij,isymab) + 551 & sab2*itoo(isymj,isymi) + isvv(isymb,isyma) 552 ioff3 = it2amt3(isymij,isymab) + 553 & sab1*isoo(isymj,isymi) + itvv(isymb,isyma) 554 555 if ((isymi .eq. isymj)) then ! isyma != isymb 556 if (isyma .lt. isymb ) then 557 nvira = nvir(isyma) 558 nvirb = 1 559 f = one 560 else 561 nvira = 1 562 nvirb = nvir(isymb) 563 f = onem 564 end if 565 do j = 1, nrhf(isymj) 566 nrhfi = j - 1 567 ioffj1 = ioff1 + ((j-1)*(j-2)/2)*sab1 568 ioffj2 = ioff2 + ((j-1)*(j-2)/2)*sab2 569 ioffj3 = ioff3 + (j*(j-1)/2)*sab1 570 do b = 1, nvir(isymb) 571 nbj = pair_pos(isymb,b,isymj,j) 572 ioffbj1 = ioffj1 + (b-1)*nvira 573 ioffbj2 = ioffj2 + (b-1)*nvira 574 ioffbj3 = ioffj3 + (b-1)*nvira 575 do i = 1, nrhfi 576 nbi = pair_pos(isymb,b,isymi,i) 577 do a = 1, nvir(isyma) 578 nai = pair_pos(isyma,a,isymi,i) 579 naj = pair_pos(isyma,a,isymj,j) 580 nabij1 = ioffbj1 + (i-1)*sab1+ 581 & (a-1)*nvirb + 1 582 nabij2 = ioffbj2 + (i-1)*sab2+ 583 & (a-1)*nvirb + 1 584 nabij3 = ioffbj3 + (i-1)*sab1+ 585 & (a-1)*nvirb + 1 586 naibj = quad_pos(isymai,nai,isymbj,nbj) 587 najbi = quad_pos(isymaj,naj,isymbi,nbi) 588 x2sq(naibj) = f*fact1*x2pack(nabij1) 589 & + fact23*( x2pack(nabij2) + f*x2pack(nabij3) ) 590 x2sq(najbi) = -f*fact1*x2pack(nabij1) 591 & - fact23*( x2pack(nabij2) - f*x2pack(nabij3) ) 592 end do ! a 593 end do ! i 594C Handle i = j, a != b 595 do a = 1, nvir(isyma) 596 naj = pair_pos(isyma,a,isymj,j) 597 nabjj3 = ioffbj3+(j-1)*sab1+(a-1)*nvirb+1 598 najbj = quad_pos(isymaj,naj,isymbj,nbj) 599 x2sq(najbj) = f*factd*x2pack(nabjj3) 600 end do 601 end do ! b 602 end do ! j 603 else if (isyma .eq. isymb ) then ! isymi != isymj 604 if (isymi .lt. isymj ) then 605 nrhfi = nrhf(isymi) 606 nrhfj = 1 607 fij = one 608 else 609 nrhfi = 1 610 nrhfj = nrhf(isymj) 611 fij = onem 612 end if 613 do j = 1, nrhf(isymj) 614 ioffj1 = ioff1 + (j-1)*nrhfi*sab1 615 ioffj2 = ioff2 + (j-1)*nrhfi*sab2 616 ioffj3 = ioff3 + (j-1)*nrhfi*sab1 617 do b = 1, nvir(isymb) 618 nbj = pair_pos(isymb,b,isymj,j) 619 nvira = b - 1 620 ioffbj1 = ioffj1 + (b-1)*(b-2)/2 621 ioffbj2 = ioffj2 + b*(b-1)/2 622 ioffbj3 = ioffj3 + (b-1)*(b-2)/2 623 624 do i = 1, nrhf(isymi) 625 nbi = pair_pos(isymb,b,isymi,i) 626 do a = 1, nvira 627 nai = pair_pos(isyma,a,isymi,i) 628 naj = pair_pos(isyma,a,isymj,j) 629 nabij1 = ioffbj1+(i-1)*nrhfj*sab1+a 630 nabij2 = ioffbj2+(i-1)*nrhfj*sab2+a 631 nabij3 = ioffbj3+(i-1)*nrhfj*sab1+a 632 naibj = quad_pos(isymai,nai,isymbj,nbj) 633 nbiaj = quad_pos(isymbi,nbi,isymaj,naj) 634 x2sq(naibj) = 635 & fij*fact1*x2pack(nabij1) 636 & + fact23*( fij*x2pack(nabij2) + 637 & x2pack(nabij3) ) 638 x2sq(nbiaj) = 639 & -fij*fact1*x2pack(nabij1) 640 & + fact23*( fij*x2pack(nabij2) - 641 & x2pack(nabij3) ) 642 end do ! a 643C Handle a == b, i<j 644 nbbij2 = ioffbj2+(i-1)*nrhfj*sab2+b 645 nbibj = quad_pos(isymbi,nbi,isymbj,nbj) 646 x2sq(nbibj) = fij*factd*x2pack(nbbij2) 647 end do ! i 648 end do ! b 649 end do ! j 650 651 else ! a != b, i != j 652 if (isyma .lt. isymb ) then 653 nvira = nvir(isyma) 654 nvirb = 1 655 fab = one 656 else 657 nvira = 1 658 nvirb = nvir(isymb) 659 fab = onem 660 end if 661 if (isymi .lt. isymj ) then 662 nrhfi = nrhf(isymi) 663 nrhfj = 1 664 fij = one 665 else 666 nrhfi = 1 667 nrhfj = nrhf(isymj) 668 fij = onem 669 end if 670 do j = 1, nrhf(isymj) 671 ioffj1 = ioff1 + (j-1)*nrhfi*sab1 672 ioffj2 = ioff2 + (j-1)*nrhfi*sab2 673 ioffj3 = ioff3 + (j-1)*nrhfi*sab1 674 do b = 1, nvir(isymb) 675 nbj = pair_pos(isymb,b,isymj,j) 676 ioffbj1 = ioffj1 + (b-1)*nvira 677 ioffbj2 = ioffj2 + (b-1)*nvira 678 ioffbj3 = ioffj3 + (b-1)*nvira 679 do i = 1, nrhf(isymi) 680 do a = 1, nvir(isyma) 681 nai = pair_pos(isyma,a,isymi,i) 682 nabij1 = ioffbj1 + 683 & (i-1)*nrhfj*sab1+ 684 & (a-1)*nvirb + 1 685 nabij2 = ioffbj2 + 686 & (i-1)*nrhfj*sab2+ 687 & (a-1)*nvirb + 1 688 nabij3 = ioffbj3 + 689 & (i-1)*nrhfj*sab1+ 690 & (a-1)*nvirb + 1 691 naibj = quad_pos(isymai,nai,isymbj,nbj) 692 x2sq(naibj) = 693 & fab*fij*fact1*x2pack(nabij1) 694 & + fact23*( fij*x2pack(nabij2) + 695 & fab*x2pack(nabij3) ) 696 end do ! a 697 end do ! i 698 end do ! b 699 end do ! j 700 701 end if 702 end do ! isymi 703 end do ! isymj 704 end if 705 END SUBROUTINE 706 707 pure function pair_pos(isyma,na,isymi,ni) 708 integer :: pair_pos 709 integer, intent(in) :: na, ni, isyma, isymi 710 711 pair_pos = it1am(isyma,isymi) + nvir(isyma)*(ni-1) + na 712 return 713 end function 714 715 pure function quad_pos(isymai,nai,isymbj,nbj) 716 integer :: quad_pos 717 integer, intent(in) :: isymai, isymbj, nai, nbj 718 quad_pos = nt1am(isymai)*(nbj-1) + nai 719 return 720 end function 721 722 END SUBROUTINE 723 724 725