1C /* Deck so_t2m1 */ 2 SUBROUTINE SO_T2M1(T2M1,LT2M1,T2MP,LT2MP,CMO,LCMO,IDEL, 3 & ISYMD,ISYDIS,WORK,LWORK) 4C 5C This routine is part of the atomic integral direct SOPPA program. 6C 7C Keld Bak and Henrik Koch, September 1995 8C Stephan P. A. Sauer: 10.11.2003: merge with Dalton 2.0 9C Rasmus Faber 2016: Rewrite for better memory access 10C 11C PURPOSE: Calculate MP2 T2-amplitudes with one back-transformed 12C index. 13C 14#include "implicit.h" 15#include "priunit.h" 16 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 17C 18 DIMENSION T2M1(LT2M1), T2MP(LT2MP), CMO(LCMO) 19 DIMENSION WORK(LWORK) 20C 21#include "ccorb.h" 22#include "ccsdinp.h" 23#include "ccsdsym.h" 24C 25C------------------ 26C Add to trace. 27C------------------ 28C 29 CALL QENTER('SO_T2M1') 30C 31 ISYMB = ISYMD 32C 33 LCDB = NVIR(ISYMB) 34CPi If no virtual orbitals in this symmetry, skip this transformation 35 IF (LCDB .EQ. 0) THEN 36 CALL QEXIT('SO_T2M1') 37 RETURN 38 END IF 39C 40 KCDB = 1 41 KEND1 = KCDB + LCDB 42 LWORK1 = LWORK - KEND1 43C 44 CALL SO_MEMMAX ('SO_T2M1.1',LWORK1) 45 IF (LWORK1 .LT. 0) CALL STOPIT('SO_T2M1.1',' ',KEND1,LWORK) 46C 47 KOFF1 = ILMVIR(ISYMD) + IDEL - IBAS(ISYMD) 48C 49C-------------------------------------------------- 50C Copy delta MO-coefficients to the vector CDB. 51C-------------------------------------------------- 52C 53 CALL DCOPY(NVIR(ISYMB),CMO(KOFF1),NBAS(ISYMD),WORK(KCDB),1) 54C 55C------------------------------------------------------------------ 56C Loop over symmetry-blocks, do reductions 57C------------------------------------------------------------------ 58C 59C$OMP PARALLEL PRIVATE(ISYMBJ,ISYMJ,ISYMAI,NAIBJ1,NAIBJ,FACT,NBJ,NBJ1, 60C$OMP& J,B,NAI,NAIBJ2,ISYMI,ISYMA,I,A,KOFFBJ,TMP,ISTART) 61 DO 100 ISYMJ = 1,NSYM 62C 63 ISYMBJ = MULD2H(ISYMJ,ISYMB) 64C 65 ISYMAI = ISYMBJ 66C 67 NAIBJ1 = IT2AM(ISYMAI,ISYMBJ) 68C$OMP DO 69 DO 120 J = 1,NRHF(ISYMJ) 70C 71 KOFF = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J-1) 72 NBJ1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) 73C 74C Zero output array 75C (only the part accessed in ai<bj loop) 76 DO NAI = 1, NBJ1+NVIR(ISYMB) 77 T2M1(KOFF+NAI) = ZERO 78 END DO 79C 80C Do reduction for this symmetry-block 81C T2 (ai<=bj) * C(b) -> T2M1 (ai,j) 82C 83 DO 130 B = 1,NVIR(ISYMB) 84C 85 NBJ = NBJ1 + B 86 FACT = WORK(KCDB-1+B) 87 NAIBJ2 = NAIBJ1 + NBJ*(NBJ-1)/2 88C 89C Stride 1 loop when ai <= bj 90 DO NAI = 1, NBJ 91 NAIBJ = NAIBJ2 + NAI 92 T2M1(KOFF+NAI) = FACT*T2MP(NAIBJ)+T2M1(KOFF+NAI) 93 END DO 94C 95 130 CONTINUE ! LOOP B 96C 97C Do the same for ai>bi : 98C T2(bj<ai)*C(b) -> T2M1 (ai,j) 99 DO ISYMI = ISYMJ, NSYM 100 ISYMA = MULD2H(ISYMAI,ISYMI) 101 IF (ISYMI .EQ. ISYMJ) THEN 102 ! Same symmetries: two cases I>J and I=J 103 ! DO I = J Seperately 104 I = J 105 DO A = 2, NVIR(ISYMA) 106 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 107 TMP = ZERO 108 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 109 ! In this case B<A, The B>=A has allready been 110 ! handled 111 DO B = 1, A-1 112 NAIBJ = KOFFBJ + B 113 TMP = TMP + T2MP(NAIBJ)*WORK(KCDB-1+B) 114 END DO 115 T2M1(KOFF+NAI) = T2M1(KOFF+NAI) + TMP 116 END DO 117 ! I > J taken care of in the following 118 ISTART = J + 1 119 ELSE 120 ! ISYMI > ISYMJ, so I>J is implicit, start loop from 121 ! one 122 ISTART = 1 123 END IF 124 125 DO I = ISTART, NRHF(ISYMI) 126 ! Since I>J, we have no restrictions on A,B 127 DO A = 1, NVIR(ISYMA) 128 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 129 TMP = ZERO 130 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 131 ! This could be a DDOT 132 DO B = 1, NVIR(ISYMB) 133 NAIBJ = KOFFBJ + B 134 TMP = TMP + T2MP(NAIBJ)*WORK(KCDB-1+B) 135 END DO 136 T2M1(KOFF+NAI) = TMP 137 END DO 138 END DO 139 END DO 140C 141 120 CONTINUE ! LOOP J 142C$OMP END DO NOWAIT 143 100 CONTINUE ! LOOP ISYMJ 144C$OMP END PARALLEL 145C 146C----------------------- 147C Remove from trace. 148C----------------------- 149C 1504321 CALL QEXIT('SO_T2M1') 151C 152 RETURN 153 END 154C /* Deck so_x2m1 */ 155 SUBROUTINE SO_X2M1(X2M1,LX2M1,X2MP,LX2MP,CMO,LCMO,IDEL, 156 & ISYMD,ISYMTR,WORK,LWORK) 157C 158C This routine is part of the atomic integral direct SOPPA program. 159C 160C Rasmus Faber 2016: Based one SO_T2M1, handles non-totally 161C symmetric vectors 162C ~ 163C PURPOSE: Calculate X2-trial vectors with one back-transformed 164C index. 165C 166C 167#include "implicit.h" 168#include "priunit.h" 169 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 170 DOUBLE PRECISION, PARAMETER :: INVSQ2 = SQRT(HALF) 171C 172 DIMENSION X2M1(LX2M1), X2MP(LX2MP), CMO(LCMO) 173 DIMENSION WORK(LWORK) 174C 175#include "ccorb.h" 176#include "ccsdinp.h" 177#include "ccsdsym.h" 178C 179C------------------ 180C Add to trace. 181C------------------ 182C 183 CALL QENTER('SO_X2M1') 184C 185 ISYMB = ISYMD 186C 187 LCDB = NVIR(ISYMB) 188C 189 KCDB = 1 190 KEND1 = KCDB + LCDB 191 LWORK1 = LWORK - KEND1 192C 193 CALL SO_MEMMAX ('SO_X2M1.1',LWORK1) 194 IF (LWORK1 .LT. 0) CALL STOPIT('SO_X2M1.1',' ',KEND1,LWORK) 195C 196 KOFF1 = ILMVIR(ISYMD) + IDEL - IBAS(ISYMD) 197C 198C-------------------------------------------------- 199C Copy delta MO-coefficients to the vector CDB. 200C-------------------------------------------------- 201C 202 CALL DCOPY(NVIR(ISYMB),CMO(KOFF1),NBAS(ISYMD),WORK(KCDB),1) 203C 204C Incorporate factor 1/sqrt(2) 205 CALL DSCAL(NVIR(ISYMB),INVSQ2,WORK(KCDB),1) 206C 207C------------------------------------------------------------------ 208C Loop over symmetry-blocks, do reductions 209C------------------------------------------------------------------ 210C 211 IF (ISYMTR .EQ. 1) THEN 212C$OMP PARALLEL PRIVATE(ISYMBJ,ISYMJ,ISYMAI,ISYMI,ISYMA, 213C$OMP& NAIBJ1,NAIBJ,NAIBJ2, 214C$OMP& FACT,KOFF,TMP,ISTART, 215C$OMP& NBJ,NBJ1,NAI, 216C$OMP& J,B,I,A) 217 DO 100 ISYMJ = 1,NSYM 218C 219 ISYMBJ = MULD2H(ISYMJ,ISYMB) 220C 221 ISYMAI = ISYMBJ 222C 223 NAIBJ1 = IT2AM(ISYMAI,ISYMBJ) 224C$OMP DO 225 DO 120 J = 1,NRHF(ISYMJ) 226C 227 KOFF = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J-1) 228 NBJ1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) 229C 230C Zero output array 231C (only the part accessed in ai<bj loop) 232 DO NAI = 1, NBJ1+NVIR(ISYMB) 233 X2M1(KOFF+NAI) = ZERO 234 END DO 235C 236C Do reduction for this symmetry-block 237C X2 (ai<=bj) * C(b) -> X2M1 (ai,j) 238C 239 DO 130 B = 1,NVIR(ISYMB) 240C 241 NBJ = NBJ1 + B 242 FACT = WORK(KCDB-1+B) 243 NAIBJ2 = NAIBJ1 + NBJ*(NBJ-1)/2 244C 245C Stride 1 loop when ai <= bj 246 DO NAI = 1, NBJ 247 NAIBJ = NAIBJ2 + NAI 248 X2M1(KOFF+NAI) = FACT*X2MP(NAIBJ)+X2M1(KOFF+NAI) 249 END DO 250C 251 130 CONTINUE ! LOOP B 252C 253C Do the same for ai>bi : 254C X2(bj<ai)*C(b) -> X2M1 (ai,j) 255 DO ISYMI = ISYMJ, NSYM 256 ISYMA = MULD2H(ISYMAI,ISYMI) 257 IF (ISYMI .EQ. ISYMJ) THEN 258 ! Same symmetries: two cases I>J and I=J 259 ! DO I = J Seperately 260 I = J 261 DO A = 2, NVIR(ISYMA) 262 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 263 TMP = ZERO 264 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 265 ! In this case B<A, The B>=A has allready been 266 ! handled 267 DO B = 1, A-1 268 NAIBJ = KOFFBJ + B 269 TMP = TMP + X2MP(NAIBJ)*WORK(KCDB-1+B) 270 END DO 271 X2M1(KOFF+NAI) = X2M1(KOFF+NAI) + TMP 272 END DO 273 ! I > J taken care of in the following 274 ISTART = J + 1 275 ELSE 276 ! ISYMI > ISYMJ, so I>J is implicit, start loop from 277 ! one 278 ISTART = 1 279 END IF 280 281 DO I = ISTART, NRHF(ISYMI) 282 ! Since I>J, we have no restrictions on A,B 283 DO A = 1, NVIR(ISYMA) 284 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 285 TMP = ZERO 286 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 287 ! This could be a DDOT 288 DO B = 1, NVIR(ISYMB) 289 NAIBJ = KOFFBJ + B 290 TMP = TMP + X2MP(NAIBJ)*WORK(KCDB-1+B) 291 END DO 292 X2M1(KOFF+NAI) = TMP 293 END DO 294 END DO 295 END DO 296C 297 120 CONTINUE ! LOOP J 298C$OMP END DO NOWAIT 299 100 CONTINUE ! LOOP ISYMJ 300C$OMP END PARALLEL 301 ELSE 302C$OMP PARALLEL PRIVATE(ISYMBJ,ISYMJ,ISYMAI, 303C$OMP& FACT,DSUM, 304C$OMP& KOFF,IOFFAIBJ,KOFFBJ,KOUT,KINP,KOFFINP, 305C$OMP& NBJ,NBJ1,NAI, 306C$OMP& J,B) 307 KOFF = 0 308 DO ISYMJ = 1, NSYM 309 ISYMBJ = MULD2H(ISYMJ,ISYMB) 310 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 311 IOFFAIBJ = IT2AM(ISYMAI,ISYMBJ) 312 KOFFBJ = IT1AM(ISYMB,ISYMJ) 313 314 IF (ISYMAI.LT.ISYMBJ) THEN 315 316 KOFFINP = IOFFAIBJ + NT1AM(ISYMAI)*KOFFBJ 317C$OMP DO 318 DO J = 1, NRHF(ISYMJ) 319 KOUT = KOFF + NT1AM(ISYMAI)*(J-1)+1 320 KINP = KOFFINP + NT1AM(ISYMAI)*NVIR(ISYMB)*(J-1)+1 321 CALL DGEMV('N',NT1AM(ISYMAI),NVIR(ISYMB),ONE, 322 & X2MP(KINP),MAX(1,NT1AM(ISYMAI)), 323 & WORK(KCDB),1,ZERO,X2M1(KOUT),1) 324 325 END DO 326C$OMP END DO NOWAIT 327 328 ELSE ! ISYMAI > ISYMBJ 329 330C$OMP DO 331 DO J = 1,NRHF(ISYMJ) 332 KOUT = KOFF + NT1AM(ISYMAI)*(J-1) 333 KOFFINP = IOFFAIBJ + KOFFBJ + NVIR(ISYMB)*(J-1) 334 DO NAI = 1,NT1AM(ISYMAI) 335 DSUM = 0.0D0 336 DO B = 1, NVIR(ISYMB) 337 DSUM = DSUM + X2MP(KOFFINP+B)*WORK(KCDB-1+B) 338 END DO 339 X2M1(KOUT+NAI) = DSUM 340 KOFFINP = KOFFINP + NT1AM(ISYMBJ) 341 END DO 342 END DO 343C$OMP END DO NOWAIT 344 END IF 345 346 KOFF = KOFF + NT1AM(ISYMAI)*NRHF(ISYMJ) 347 348 END DO 349C$OMP END PARALLEL 350 ENDIF 351C 352C----------------------- 353C Remove from trace. 354C----------------------- 355C 356 CALL QEXIT('SO_X2M1') 357C 358 RETURN 359 END 360C /* Deck so_x2sm1 */ 361 SUBROUTINE SO_X2SM1(X2M1,LX2M1,X2SQ,LX2SQ,CMO,LCMO,IDEL, 362 & ISYMD,ISYMTR,WORK,LWORK) 363C 364C This routine is part of the atomic integral direct SOPPA program. 365C 366C Rasmus Faber 2017: Version of SO_X2M1, that handles squared, non-totally 367C symmetric vectors 368C ~ 369C PURPOSE: Calculate X2-trial vectors with one back-transformed 370C index. 371C 372C 373#include "implicit.h" 374#include "priunit.h" 375 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 376 DOUBLE PRECISION, PARAMETER :: INVSQ2 = SQRT(HALF) 377 character(len=*), parameter :: myname = 'SO_X2SM1' 378C 379 DIMENSION X2M1(LX2M1), X2SQ(LX2SQ), CMO(LCMO) 380 DIMENSION WORK(LWORK) 381C 382#include "ccorb.h" 383#include "ccsdinp.h" 384#include "ccsdsym.h" 385C 386C------------------ 387C Add to trace. 388C------------------ 389C 390 CALL QENTER(myname) 391C 392 ISYMB = ISYMD 393C 394 LCDB = NVIR(ISYMB) 395C 396 KCDB = 1 397 KEND1 = KCDB + LCDB 398 LWORK1 = LWORK - KEND1 399C 400 CALL SO_MEMMAX (myname//'.1',LWORK1) 401 IF (LWORK1 .LT. 0) CALL STOPIT(myname//'.1',' ',KEND1,LWORK) 402C 403 KOFF1 = ILMVIR(ISYMD) + IDEL - IBAS(ISYMD) 404C 405C-------------------------------------------------- 406C Copy delta MO-coefficients to the vector CDB. 407C-------------------------------------------------- 408C 409 CALL DCOPY(NVIR(ISYMB),CMO(KOFF1),NBAS(ISYMD),WORK(KCDB),1) 410C 411C Incorporate factor 1/sqrt(2) 412 CALL DSCAL(NVIR(ISYMB),INVSQ2,WORK(KCDB),1) 413C 414C------------------------------------------------------------------ 415C Loop over symmetry-blocks, do reductions 416C------------------------------------------------------------------ 417C 418 KOFF = 0 419 DO ISYMJ = 1, NSYM 420 ISYMBJ = MULD2H(ISYMJ,ISYMB) 421 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 422 IOFFAIBJ = IT2SQ(ISYMAI,ISYMBJ) 423 KOFFBJ = IT1AM(ISYMB,ISYMJ) 424 425 426 KOFFINP = IOFFAIBJ + NT1AM(ISYMAI)*KOFFBJ 427 DO J = 1, NRHF(ISYMJ) 428 KOUT = KOFF + NT1AM(ISYMAI)*(J-1)+1 429 KINP = KOFFINP + NT1AM(ISYMAI)*NVIR(ISYMB)*(J-1)+1 430 CALL DGEMV('N',NT1AM(ISYMAI),NVIR(ISYMB),ONE, 431 & X2SQ(KINP),MAX(1,NT1AM(ISYMAI)), 432 & WORK(KCDB),1,ONE,X2M1(KOUT),1) 433 END DO 434 435 436 KOFF = KOFF + NT1AM(ISYMAI)*NRHF(ISYMJ) 437 438 END DO 439C 440C----------------------- 441C Remove from trace. 442C----------------------- 443C 444 CALL QEXIT(myname) 445C 446 RETURN 447 END 448 449C /* Deck so_t2x1 */ 450 SUBROUTINE SO_T2X1(X2M1,LX2M1,T2MP,LT2MP,TJ1,LTJ1,IDEL, 451 & ISYMD,ISYMTR,WORK,LWORK) 452C 453C Rasmus Faber, 2016: Bases on SO_T2M1 454C 455C PURPOSE: Calculates a T2(ai, bj) * x1( b, delta) intermediate and 456C and adds it to the back-transformed x2 trial-vectors. 457C This allows term (1) and (5) of the B-matrix to be 458C calculated using SO_RES_TCB 459C 460#include "implicit.h" 461#include "priunit.h" 462 PARAMETER (ZERO = 0.0D0, HALF = 0.5D0, ONE = 1.0D0, TWO = 2.0D0) 463C 464 DIMENSION X2M1(LX2M1), T2MP(LT2MP), TJ1(LTJ1) 465 DIMENSION WORK(LWORK) 466C 467#include "ccorb.h" 468#include "ccsdinp.h" 469#include "ccsdsym.h" 470C 471C------------------ 472C Add to trace. 473C------------------ 474C 475 CALL QENTER('SO_T2X1') 476C 477 ISYMB = MULD2H(ISYMD,ISYMTR) 478C 479 LCDB = NVIR(ISYMB) 480CRF If there's nothing to do here, skip straight to exit 481 IF ( LCDB .LE. 0 ) GOTO 4321 482C 483 KCDB = 1 484 KEND1 = KCDB + LCDB 485 LWORK1 = LWORK - KEND1 486C 487 CALL SO_MEMMAX ('SO_T2X1.1',LWORK1) 488 IF (LWORK1 .LT. 0) CALL STOPIT('SO_T2X1.1',' ',KEND1,LWORK) 489C 490 KOFF1 = IMATAV(ISYMD,ISYMB) + IDEL - IBAS(ISYMD) 491C 492C-------------------------------------------------- 493C Copy delta X1-coefficients to the vector CDB. 494C-------------------------------------------------- 495C 496 CALL DCOPY(NVIR(ISYMB),TJ1(KOFF1),NBAS(ISYMD),WORK(KCDB),1) 497C 498C------------------------------------------------------------------ 499C Loop over symmetry-blocks, do reductions 500C------------------------------------------------------------------ 501C 502C$OMP PARALLEL PRIVATE(ISYMBJ,ISYMJ,ISYMAI,NAIBJ1,NAIBJ,FACT,NBJ,NBJ1, 503C$OMP& J,B,NAI,NAIBJ2,ISYMI,ISYMA,I,A,KOFFBJ,TMP,ISTART) 504 DO 100 ISYMJ = 1,NSYM 505C 506 ISYMBJ = MULD2H(ISYMJ,ISYMB) 507C 508 ISYMAI = ISYMBJ 509C 510 NAIBJ1 = IT2AM(ISYMAI,ISYMBJ) 511C$OMP DO 512 DO 120 J = 1,NRHF(ISYMJ) 513C 514 KOFF = IT2BCD(ISYMAI,ISYMJ) + NT1AM(ISYMAI)*(J-1) 515 NBJ1 = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J-1) 516C 517C T2 (ai<=bj) * C(b) -> T2M1 (ai,j) 518C 519 DO 130 B = 1,NVIR(ISYMB) 520C 521 NBJ = NBJ1 + B 522 FACT = WORK(KCDB-1+B) 523 NAIBJ2 = NAIBJ1 + NBJ*(NBJ-1)/2 524C 525C Stride 1 loop when ai <= bj 526 DO NAI = 1, NBJ 527 NAIBJ = NAIBJ2 + NAI 528 X2M1(KOFF+NAI) = FACT*T2MP(NAIBJ)+X2M1(KOFF+NAI) 529 END DO 530C 531 130 CONTINUE ! LOOP B 532C 533C Do the same for ai>bi : 534C T2(bj<ai)*C(b) -> T2M1 (ai,j) 535 DO ISYMI = ISYMJ, NSYM 536 ISYMA = MULD2H(ISYMAI,ISYMI) 537 IF (ISYMI .EQ. ISYMJ) THEN 538 ! Same symmetries: two cases I>J and I=J 539 ! DO I = J Seperately 540 I = J 541 DO A = 2, NVIR(ISYMA) 542 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 543 TMP = ZERO 544 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 545 ! In this case B<A, The B>=A has allready been 546 ! handled 547 DO B = 1, A-1 548 NAIBJ = KOFFBJ + B 549 TMP = TMP + T2MP(NAIBJ)*WORK(KCDB-1+B) 550 END DO 551 X2M1(KOFF+NAI) = X2M1(KOFF+NAI) + TMP 552 END DO 553 ! I > J taken care of in the following 554 ISTART = J + 1 555 ELSE 556 ! ISYMI > ISYMJ, so I>J is implicit, start loop from 557 ! one 558 ISTART = 1 559 END IF 560C 561 DO I = ISTART, NRHF(ISYMI) 562 ! Since I>J, we have no restrictions on A,B 563 DO A = 1, NVIR(ISYMA) 564 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 565 TMP = ZERO 566 KOFFBJ = NAIBJ1 + NAI*(NAI-1)/2 + NBJ1 567 ! This could be a DDOT 568 DO B = 1, NVIR(ISYMB) 569 NAIBJ = KOFFBJ + B 570 TMP = TMP + T2MP(NAIBJ)*WORK(KCDB-1+B) 571 END DO 572 X2M1(KOFF+NAI) = X2M1(KOFF+NAI) + TMP 573 END DO 574 END DO 575 END DO 576C 577 120 CONTINUE ! LOOP J 578C$OMP END DO NOWAIT 579 100 CONTINUE ! LOOP ISYMJ 580C$OMP END PARALLEL 581C 582C----------------------- 583C Remove from trace. 584C----------------------- 585C 5864321 CALL QEXIT('SO_T2X1') 587C 588 RETURN 589 END 590 591 SUBROUTINE SO_M1SHUF(T2M1,LT2M1,ISYMD,ISYMTR) 592C 593C Calculate 594C v ~ ~ 595C t(ai,j) = -1/3 ( t(ai,j) + 2 t(aj,i) ) 596C in place. 597 use so_info, only: sop_dp 598 implicit none 599 600#include "ccorb.h" 601#include "ccsdinp.h" 602#include "ccsdsym.h" 603 real(sop_dp), intent(inout) :: T2M1(LT2M1) 604 integer, intent(in) :: LT2M1, ISYMD, ISYMTR 605 606 character(len=*), parameter :: myname = 'SO_M1SHUF' 607 real(sop_dp), parameter :: inv3m = -1.0D0/3.0D0, two = 2.0D0 608 609 integer :: isymi, isymj, isymai, isyma, isymjd, isymaj 610 integer :: ioffj, ioffi, ioffaij, ioffaji 611 integer :: nrhfi 612 real(sop_dp) :: xaij, xaji 613 614 CALL QENTER(MYNAME) 615 616 DO ISYMJ = 1, NSYM 617 isymjd = muld2h(isymj,isymd) 618 do ISYMI = 1, isymj ! Ensure i < j 619 isymai = muld2h(isymjd,isymtr) 620 isyma = muld2h(isymai,isymi) 621 isymaj = muld2h(isymj,isyma) 622 623 do j = 1, nrhf(isymj) 624 ioffj = it2bcd(isymai,isymj) + nt1am(isymai)*(j-1) 625 if (isymi.eq.isymj) then 626 nrhfi = j 627 else 628 nrhfi = nrhf(isymi) 629 end if 630 do i = 1, nrhfi 631 ioffi = it2bcd(isymaj,isymi)+nt1am(isymaj)*(i-1) 632 ioffaij = ioffj + it1am(isyma,isymi) 633 & + nvir(isyma)*(i-1) 634 ioffaji = ioffi + it1am(isyma,isymj) 635 & + nvir(isyma)*(j-1) 636 do a = 1, nvir(isyma) 637 xaij = t2m1(ioffaij+a) 638 xaji = t2m1(ioffaji+a) 639 t2m1(ioffaij+a) = inv3m * (xaij+two*xaji) 640 t2m1(ioffaji+a) = inv3m * (xaji+two*xaij) 641 end do ! a 642 end do ! i 643 end do ! j 644 end do ! isymi 645 END DO ! isymj 646 647 CALL QEXIT(MYNAME) 648 END SUBROUTINE 649