1C /* Deck cc_choatr */ 2 SUBROUTINE CC_CHOATR(SIGMA1,X1AM,FREQ,WORK,LWORK,ISYMX,NUMX, 3 & ISIDE) 4C 5C Thomas Bondo Pedersen, February 2003. 6C 7C Purpose: Calculate NUMX linear transformations with the effective CC2 8C Jacobian from left (ISIDE=-1) or right (ISIDE=1). 9C 10C ISIDE = -1: SIGMA1 = X1AM * A(FREQ) [NOTE THE SIGN OF FREQ] 11C 12C ISIDE = 1: SIGMA1 = A(FREQ) * X1AM 13C 14C The (global) E intermediates must be available on disk (files opened 15C through subroutine CHO_IMOP). 16C 17#include "implicit.h" 18 DIMENSION SIGMA1(*), X1AM(*), FREQ(NUMX), WORK(LWORK) 19#include "ccorb.h" 20#include "ccsdsym.h" 21#include "priunit.h" 22 23 LOGICAL LOCDBG 24 PARAMETER (LOCDBG = .FALSE.) 25 26 CHARACTER*9 SECNAM 27 PARAMETER (SECNAM = 'CC_CHOATR') 28 29C Debug: print. 30C ------------- 31 32 IF (LOCDBG) THEN 33 WRITE(LUPRI,'(/,5X,A,A,I10)') 34 & SECNAM,': ISIDE = ',ISIDE 35 WRITE(LUPRI,'(5X,A,A,I10)') 36 & SECNAM,': ISYMX = ',ISYMX 37 WRITE(LUPRI,'(5X,A,A,I10)') 38 & SECNAM,': NUMX = ',NUMX 39 WRITE(LUPRI,'(5X,A,A)') 40 & SECNAM,': Frequencies:' 41 WRITE(LUPRI,'(1P,D20.12)') (FREQ(I), I = 1,NUMX) 42 WRITE(LUPRI,*) 43 ENDIF 44 45C Calculate contribution from global E intermediates 46C (initialization of result vectors in SIGMA1). 47C --------------------------------------------------- 48 49 CALL CC_CHOATR0(SIGMA1,X1AM,WORK,LWORK,ISYMX,NUMX,ISIDE) 50 51 IF (LOCDBG) THEN 52 DO I = 1,NUMX 53 KOFF = NT1AM(ISYMX)*(I - 1) + 1 54 SNRM = DSQRT(DDOT(NT1AM(ISYMX),SIGMA1(KOFF),1, 55 & SIGMA1(KOFF),1)) 56 WRITE(LUPRI,'(5X,A,I3,A,1P,D22.15)') 57 & 'Norm of transformed vector',I,' after CC_CHOATR0: ',SNRM 58 ENDDO 59 ENDIF 60 61C Read the canonical orbital energies and calculate Lambda matrices. 62C ------------------------------------------------------------------ 63 64 KFOCKD = 1 65 KLAMDP = KFOCKD + NORBTS 66 KLAMDH = KLAMDP + NGLMDT(1) 67 KT1AM = KLAMDH + NGLMDT(1) 68 KEND1 = KT1AM + NT1AM(1) 69 LWRK1 = LWORK - KEND1 + 1 70 71 IF (LWRK1 .LT. 0) THEN 72 WRITE(LUPRI,'(//,5X,A,A,A)') 73 & 'Insufficient memory in ',SECNAM,' - allocation: orb. en.' 74 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 75 & 'Need (more than): ',KEND1-1, 76 & 'Available : ',LWORK 77 CALL QUIT('Insufficient memory in '//SECNAM) 78 ENDIF 79 80 CALL CHO_RDSIR(DUM1,DUM2,WORK(KFOCKD),DUM3,WORK(KEND1),LWRK1, 81 & .FALSE.,.TRUE.,.FALSE.) 82 IFAIL = -1 83 CALL CHO_RDRST(DUM1,WORK(KT1AM),DUM2,.FALSE.,.TRUE.,.FALSE.,IFAIL) 84 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM),WORK(KEND1), 85 & LWRK1) 86 87C Read unperturbed F(ia) into T1AM (no longer needed). 88C ---------------------------------------------------- 89 90 KFIA = KT1AM 91 CALL ONEL_OP(-1,3,LUFIA) 92 CALL CHO_MOREAD(WORK(KFIA),NT1AM(1),1,1,LUFIA) 93 CALL ONEL_OP(1,3,LUFIA) 94 95C Calculate C intermediates for ISIDE = -1, or 96C calculate the I2J contribution for ISIDE = 1. 97C --------------------------------------------- 98 99 IF (ISIDE .EQ. -1) THEN 100 CALL CC_CHOCIM(WORK(KFOCKD),X1AM,WORK(KEND1),LWRK1,ISYMX,NUMX) 101 ELSE 102 CALL CC_CHORTRI2J(WORK(KFOCKD),WORK(KLAMDP),WORK(KLAMDH), 103 & SIGMA1,X1AM,WORK(KEND1),LWRK1,ISYMX,NUMX) 104 IF (LOCDBG) THEN 105 DO I = 1,NUMX 106 KOFF = NT1AM(ISYMX)*(I - 1) + 1 107 SNRM = DSQRT(DDOT(NT1AM(ISYMX),SIGMA1(KOFF),1, 108 & SIGMA1(KOFF),1)) 109 WRITE(LUPRI,'(5X,A,I3,A,1P,D22.15)') 110 & 'Norm of transformed vector',I,' after CC_CHORTRI2J: ', 111 & SNRM 112 ENDDO 113 ENDIF 114 ENDIF 115 116C Loop over trial vectors for calculating doubles-dependent terms. 117C ================================================================ 118 119 DO IX = 1,NUMX 120 121C Offset to vectors in SIGMA1 and X1AM. 122C ------------------------------------- 123 124 KOFIX = NT1AM(ISYMX)*(IX - 1) + 1 125 126C Calculate perturbed Cholesky vectors, store on disk. 127C ---------------------------------------------------- 128 129 CALL CC_CHOTG(WORK(KLAMDP),WORK(KLAMDH),X1AM(KOFIX), 130 & WORK(KEND1),LWRK1,1,ISYMX,ISIDE) 131 132C Calculate Y intermediates. For ISIDE = 1, also 133C calculate the I1 contribution. Perturbed Cholesky 134C vector files are deleted at the end of CC_CYI. 135C ------------------------------------------------- 136 137c CALL CC_CYI(WORK(KFOCKD),SIGMA1(KOFIX),X1AM(KOFIX),WORK(KFIA), 138c & WORK(KEND1),LWRK1,ISIDE,ISYMX,1,FREQ(IX),X2NRM, 139c & X2CNM,.TRUE.) 140 141 IF (ISIDE .EQ. -1) THEN 142 CALL CC_CYILTR(WORK(KFOCKD),X1AM(KOFIX),WORK(KFIA), 143 & WORK(KEND1),LWRK1,ISYMX,FREQ(IX),X2NRM, 144 & X2CNM,.FALSE.,.TRUE.) 145 ELSE IF (ISIDE .EQ. 1) THEN 146 CALL CC_CYIRTR(WORK(KFOCKD),SIGMA1(KOFIX),WORK(KFIA), 147 & WORK(KEND1),LWRK1,ISYMX,FREQ(IX),X2NRM, 148 & X2CNM,.FALSE.,.TRUE.) 149 ELSE 150 CALL QUIT('ISIDE error in '//SECNAM) 151 ENDIF 152 153 IF (LOCDBG) THEN 154 SNRM = DSQRT(DDOT(NT1AM(ISYMX),SIGMA1(KOFIX),1, 155 & SIGMA1(KOFIX),1)) 156 WRITE(LUPRI,'(5X,A,I3,A,1P,D22.15)') 157 & 'Norm of transformed vector',IX,' after CC_CYI: ',SNRM 158 ENDIF 159 160C For ISIDE=-1: Calculate GIJ and H terms. 161C For ISIDE=+1: Calculate G and H terms. 162C Delete Y intermediate files after use. 163C ----------------------------------------- 164 165 CALL CC_CHOATR1(WORK(KLAMDP),WORK(KLAMDH),SIGMA1(KOFIX), 166 & X1AM(KOFIX),WORK(KEND1),LWRK1,ISYMX,ISIDE, 167 & .TRUE.,IX) 168 169 IF (LOCDBG) THEN 170 SNRM = DSQRT(DDOT(NT1AM(ISYMX),SIGMA1(KOFIX),1, 171 & SIGMA1(KOFIX),1)) 172 WRITE(LUPRI,'(5X,A,I3,A,1P,D22.15)') 173 & 'Norm of transformed vector',IX,' after CC_CHOATR1: ',SNRM 174 ENDIF 175 176 ENDDO 177 178 RETURN 179 END 180C /* Deck cc_cyirtr */ 181 SUBROUTINE CC_CYIRTR(FOCKD,SIGMA1,FIA,WORK,LWORK,ISYMX,FREQ,X2NRM, 182 & X2CNM,DELP1,DELP2) 183C 184C Thomas Bondo Pedersen, February 2003. 185C 186C Purpose: Calculate Y intermediates for right-hand Jacobian 187C transformation. 188C 189#include "implicit.h" 190 DIMENSION FOCKD(*), SIGMA1(*), FIA(*), WORK(LWORK) 191 LOGICAL DELP1, DELP2 192#include "maxorb.h" 193#include "ccdeco.h" 194#include "ccorb.h" 195#include "ccsdsym.h" 196#include "priunit.h" 197#include "dummy.h" 198 199C Set Cholesky vector type for amplitude assembly. 200C ------------------------------------------------ 201 202 FAC1 = 1.0D0 203 FAC2 = 0.0D0 204 SCD = 1.0D0 205 SCDG = 1.0D0 206 207 NUMP12 = 1 208 JTYP1 = 3 209 JTYP2 = 9 210 ISYMP1 = 1 211 ISYMP2 = ISYMX 212 213 IOPTDN = 1 214 IOPTCE = 1 215 216C Set info for Y intermediates. 217C ----------------------------- 218 219 NUMZ = 1 220 JTYPZ = 1 221 ISYMZ = 1 222 JTYPY = 1 223 224C Calculate Y intermediates. 225C -------------------------- 226 227 CALL CC_CYI(JTYP1,ISYMP1,JTYP2,ISYMP2,NUMCHO,FAC1,NUMP12, 228 & DUMMY,IDUMMY,DUMMY,IDUMMY,FAC2,0, 229 & SCD,SCDG,IOPTDN,FOCKD,FREQ,IOPTCE, 230 & JTYPZ,ISYMZ,NUMCHO,NUMZ,JTYPY, 231 & FIA,1,1,SIGMA1, 232 & WORK,LWORK,X2NRM,X2CNM,DELP1,DELP2) 233 234 RETURN 235 END 236C /* Deck cc_cyiltr */ 237 SUBROUTINE CC_CYILTR(FOCKD,X1AM,FIA,WORK,LWORK,ISYMX,FREQ,X2NRM, 238 & X2CNM,DELP1,DELP2) 239C 240C Thomas Bondo Pedersen, February 2003. 241C 242C Purpose: Calculate Y intermediates for left-hand Jacobian 243C transformation. 244C 245#include "implicit.h" 246 DIMENSION FOCKD(*), X1AM(*), FIA(*), WORK(LWORK) 247 LOGICAL DELP1, DELP2 248#include "maxorb.h" 249#include "ccdeco.h" 250#include "ccorb.h" 251#include "ccsdsym.h" 252#include "priunit.h" 253#include "dummy.h" 254 255C Set Cholesky vector type for amplitude assembly. 256C ------------------------------------------------ 257 258 FAC1 = 1.0D0 259 FAC2 = 1.0D0 260 SCD = 1.0D0 261 SCDG = 1.0D0 262 263 NUMP12 = 1 264 JTYP1 = 1 265 JTYP2 = 8 266 ISYMP1 = 1 267 ISYMP2 = ISYMX 268 269 IOPTDN = 1 270 IOPTCE = 1 271 272C Set info for Y intermediates. 273C ----------------------------- 274 275 NUMZ = 1 276 JTYPZ = 3 277 ISYMZ = 1 278 JTYPY = -1 279 280C Calculate Y intermediates. 281C -------------------------- 282 283 CALL CC_CYI(JTYP1,ISYMP1,JTYP2,ISYMP2,NUMCHO,FAC1,NUMP12, 284 & X1AM,ISYMX,FIA,1,FAC2,1, 285 & SCD,SCDG,IOPTDN,FOCKD,FREQ,IOPTCE, 286 & JTYPZ,ISYMZ,NUMCHO,NUMZ,JTYPY, 287 & DUMMY,IDUMMY,0,DUMMY, 288 & WORK,LWORK,X2NRM,X2CNM,DELP1,DELP2) 289 290 RETURN 291 END 292C /* Deck cc_choatr1 */ 293 SUBROUTINE CC_CHOATR1(XLAMDP,XLAMDH,SIGMA1,X1AM,WORK,LWORK,ISYMX, 294 & ISIDE,DELY,IX) 295C 296C Thomas Bondo Pedersen, February 2003. 297C 298C Purpose: Calculate the GIJ and H terms for ISIDE=-1, or 299C the G and H terms for ISIDE=1. 300C 301#include "implicit.h" 302 DIMENSION XLAMDP(*), XLAMDH(*), SIGMA1(*), X1AM(*), WORK(LWORK) 303 LOGICAL DELY 304#include "ccorb.h" 305#include "ccsdsym.h" 306#include "priunit.h" 307 308 CHARACTER*10 SECNAM 309 PARAMETER (SECNAM = 'CC_CHOATR1') 310 311C Call appropriate routine. 312C ------------------------- 313 314 IF (ISIDE .EQ. -1) THEN 315 316 KCIM = 1 317 KEND1 = KCIM + NT1AM(ISYMX) 318 LWRK1 = LWORK - KEND1 + 1 319 320 IF (LWRK1 .LT. 0) THEN 321 WRITE(LUPRI,'(//,5X,A,A)') 322 & 'Insufficient memory in ',SECNAM 323 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 324 & 'Need (more than): ',KEND1-1, 325 & 'Available : ',LWORK 326 CALL QUIT('Insufficient memory in '//SECNAM) 327 ENDIF 328 329 CALL CHO_IMOP(-1,3,LUCIM,ISYMX) 330 CALL CHO_MOREAD(WORK(KCIM),NT1AM(ISYMX),1,IX,LUCIM) 331 CALL CHO_IMOP(1,3,LUCIM,ISYMX) 332 333 CALL CC_CHOLTRGIJH(XLAMDP,XLAMDH,SIGMA1,X1AM,WORK(KCIM), 334 & WORK(KEND1),LWRK1,ISYMX) 335 336 ELSE IF (ISIDE .EQ. 1) THEN 337 338 CALL CC_CHORTRGH(XLAMDP,XLAMDH,SIGMA1,WORK,LWORK,ISYMX) 339 340 ELSE 341 342 WRITE(LUPRI,'(//,5X,A,A)') 343 & 'ISIDE must be -1 or +1 in ',SECNAM 344 WRITE(LUPRI,'(5X,A,I10,/)') 345 & 'ISIDE = ',ISIDE 346 CALL QUIT('Error in '//SECNAM) 347 348 ENDIF 349 350C Delete Y intermediate files if requested. 351C ----------------------------------------- 352 353 IF (DELY) THEN 354 DO ISYCHO = 1,NSYM 355 CALL CC_CYIOP(-1,ISYCHO,ISIDE) 356 CALL CC_CYIOP(0,ISYCHO,ISIDE) 357 ENDDO 358 ENDIF 359 360 RETURN 361 END 362C /* Deck cc_chortrgh */ 363 SUBROUTINE CC_CHORTRGH(XLAMDP,XLAMDH,SIGMA1,WORK,LWORK,ISYMY) 364C 365C Thomas Bondo Pedersen, February 2003. 366C 367C Purpose: Calculate the G and H terms for right-hand Jacobian 368C transformations. 369C 370C Formula: 371C 372C SIGMA1(ai) = SIGMA1(ai) 373C + sum(g) LamP(g,a) sum(Jd) L(gd,J) sum(b) LamH(d,b) * Y(bi,J) 374C - sum(Jj) Y(aj,J) * L(ji,J) 375C 376C where g,d refer to AO basis. 377C 378#include "implicit.h" 379 DIMENSION XLAMDP(*), XLAMDH(*), SIGMA1(*), WORK(LWORK) 380#include "maxorb.h" 381#include "ccdeco.h" 382#include "ccorb.h" 383#include "ccsdsym.h" 384#include "priunit.h" 385 386 INTEGER IOFFL(8), IOFFZ(8), IOFFY(8) 387 388 CHARACTER*11 SECNAM 389 PARAMETER (SECNAM = 'CC_CHORTRGH') 390 391 PARAMETER (IOPTR = 2) 392 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 393 394C Read reduce index array. 395C ------------------------ 396 397 KIND1 = 1 398 CALL CC_GETIND1(WORK(KIND1),LWORK,LIND1) 399 KEND0 = KIND1 + LIND1 400 LWRK0 = LWORK - KEND0 + 1 401 402 IF (LWRK0 .LT. 0) THEN 403 WRITE(LUPRI,'(//,5X,A,A,A)') 404 & 'Insufficient memory in ',SECNAM,' - allocation: index' 405 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 406 & 'Need (more than): ',KEND0-1, 407 & 'Available : ',LWORK 408 CALL QUIT('Insufficient memory in '//SECNAM) 409 ENDIF 410 411C Allocation. 412C ----------- 413 414 KSIGMA = KEND0 415 KEND1 = KSIGMA + NT1AO(ISYMY) 416 LWRK1 = LWORK - KEND1 + 1 417 418 IF (LWRK1 .LT. 0) THEN 419 WRITE(LUPRI,'(//,5X,A,A)') 420 & 'Insufficient memory in ',SECNAM 421 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 422 & 'Need (more than): ',KEND1-1, 423 & 'Available : ',LWORK 424 CALL QUIT('Insufficient memory in '//SECNAM) 425 ENDIF 426 427C Initialize AO G term. 428C --------------------- 429 430 CALL DZERO(WORK(KSIGMA),NT1AO(ISYMY)) 431 432C Start Cholesky symmetry loop. 433C ----------------------------- 434 435 DO ISYCHO = 1,NSYM 436 437 IF (NUMCHO(ISYCHO) .LE. 0) GO TO 1000 438 439 ISYMBI = MULD2H(ISYCHO,ISYMY) 440 ISYMAJ = ISYMBI 441 ISYMDI = ISYMBI 442 443C Allocation for one AO Cholesky vector. 444C -------------------------------------- 445 446 LREAD = MEMRD(1,ISYCHO,IOPTR) 447 448 KCHOL = KEND1 449 KREAD = KCHOL + N2BST(ISYCHO) 450 KEND2 = KREAD + LREAD 451 LWRK2 = LWORK - KEND2 + 1 452 453 IF (LWRK2 .LE. 0) THEN 454 WRITE(LUPRI,'(//,5X,A,A)') 455 & 'Insufficient memory in ',SECNAM 456 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 457 & 'Need (more than): ',KEND2-1, 458 & 'Available : ',LWORK 459 CALL QUIT('Insufficient memory in '//SECNAM) 460 ENDIF 461 462C Set up batch. 463C SCRL: L(ji,#J),Y(b#J,i),L(g,d#J). 464C SCRZ: Y(bi,#J),Z(d#J,i). 465C --------------------------------- 466 467 LENL = MAX(NMATIJ(ISYCHO),NT1AM(ISYMBI),N2BST(ISYCHO)) 468 LENZ = MAX(NT1AM(ISYMBI),NT1AO(ISYMDI)) 469 470 MINMEM = LENL + LENZ 471 IF (MINMEM .GT. 0) THEN 472 NVEC = MIN(LWRK2/MINMEM,NUMCHO(ISYCHO)) 473 ELSE IF (MINMEM .EQ. 0) THEN 474 GO TO 1000 475 ELSE 476 WRITE(LUPRI,'(//,5X,A,A,A,I10,/)') 477 & 'MINMEM is negative in ',SECNAM,': ',MINMEM 478 CALL QUIT('Error in '//SECNAM) 479 ENDIF 480 481 IF (NVEC .LE. 0) THEN 482 WRITE(LUPRI,'(//,5X,A,A)') 483 & 'Insufficient memory for Cholesky batch in ',SECNAM 484 WRITE(LUPRI,'(5X,A,I2,A)') 485 & '(Cholesky symmetry:',ISYCHO,')' 486 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10,/)') 487 & 'Minimum memory required for batch: ',MINMEM, 488 & 'Memory available for batch : ',LWRK2, 489 & 'Total memory available : ',LWORK 490 CALL QUIT('Insufficient memory in '//SECNAM) 491 ENDIF 492 493 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 494 495C Open Y intermediate file. 496C ------------------------- 497 498 CALL CC_CYIOP(-1,ISYCHO,1) 499 500C Open occupied Cholesky vector file. 501C ----------------------------------- 502 503 CALL CHO_MOP(-1,4,ISYCHO,LUCHJI,1,1) 504 505C Start batch loop. 506C ----------------- 507 508 DO IBATCH = 1,NBATCH 509 510 NUMV = NVEC 511 IF (IBATCH .EQ. NBATCH) THEN 512 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 513 ENDIF 514 JVEC1 = NVEC*(IBATCH - 1) + 1 515 516C Set up local index arrays. 517C -------------------------- 518 519 ICOUNL = 0 520 ICOUNZ = 0 521 ICOUNY = 0 522 DO ISYMD = 1,NSYM 523 ISYMI = MULD2H(ISYMD,ISYMDI) 524 ISYMG = MULD2H(ISYMD,ISYCHO) 525 ISYMB = ISYMD 526 IOFFL(ISYMD) = ICOUNL 527 IOFFZ(ISYMD) = ICOUNZ 528 IOFFY(ISYMB) = ICOUNY 529 LENDJ = NBAS(ISYMD)*NUMV 530 LENBJ = NVIR(ISYMB)*NUMV 531 ICOUNL = ICOUNL + NBAS(ISYMG)*LENDJ 532 ICOUNZ = ICOUNZ + LENDJ*NRHF(ISYMI) 533 ICOUNY = ICOUNY + LENBJ*NRHF(ISYMI) 534 ENDDO 535 536C Complete allocation. 537C -------------------- 538 539 KSCRL = KEND2 540 KSCRZ = KSCRL + LENL*NUMV 541 542C Read Y(bi,#J) in place of Z. 543C ---------------------------- 544 545 CALL CC_CYIRDF(WORK(KSCRZ),NUMV,JVEC1,ISYCHO,ISYMY,1) 546 547C Read L(ji,#J) in place of L(g,d#J). 548C ----------------------------------- 549 550 CALL CHO_MOREAD(WORK(KSCRL),NMATIJ(ISYCHO),NUMV,JVEC1, 551 & LUCHJI) 552 553C Calculate H term. 554C ----------------- 555 556 DO IVEC = 1,NUMV 557 DO ISYMJ = 1,NSYM 558 559 ISYMI = MULD2H(ISYMJ,ISYCHO) 560 ISYMA = MULD2H(ISYMJ,ISYMAJ) 561 562 KOFF1 = KSCRZ + NT1AM(ISYMAJ)*(IVEC - 1) 563 & + IT1AM(ISYMA,ISYMJ) 564 KOFF2 = KSCRL + NMATIJ(ISYCHO)*(IVEC - 1) 565 & + IMATIJ(ISYMJ,ISYMI) 566 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 567 568 NTOTA = MAX(NVIR(ISYMA),1) 569 NTOTJ = MAX(NRHF(ISYMJ),1) 570 571 CALL DGEMM('N','N', 572 & NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 573 & XMONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOTJ, 574 & ONE,SIGMA1(KOFF3),NTOTA) 575 576 ENDDO 577 ENDDO 578 579C Resort Y intermediate as Y(b#J,i), store in place of L(g,d#J). 580C -------------------------------------------------------------- 581 582 DO IVEC = 1,NUMV 583 DO ISYMI = 1,NSYM 584 585 ISYMB = MULD2H(ISYMI,ISYMBI) 586 LENBJ = NVIR(ISYMB)*NUMV 587 588 DO I = 1,NRHF(ISYMI) 589 590 KOFF1 = KSCRZ + NT1AM(ISYMBI)*(IVEC - 1) 591 & + IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) 592 KOFF2 = KSCRL + IOFFY(ISYMB) + LENBJ*(I - 1) 593 & + NVIR(ISYMB)*(IVEC - 1) 594 595 CALL DCOPY(NVIR(ISYMB),WORK(KOFF1),1,WORK(KOFF2),1) 596 597 ENDDO 598 599 ENDDO 600 ENDDO 601 602C Backtransform the Y intermediates to get Z(d#J,i). 603C -------------------------------------------------- 604 605 DO ISYMB = 1,NSYM 606 607 ISYMI = MULD2H(ISYMB,ISYMBI) 608 ISYMD = ISYMB 609 610 KOFF1 = IGLMVI(ISYMD,ISYMB) + 1 611 KOFF2 = KSCRL + IOFFY(ISYMB) 612 KOFF3 = KSCRZ + IOFFZ(ISYMD) 613 614 NTOTD = MAX(NBAS(ISYMD),1) 615 NTOTB = MAX(NVIR(ISYMB),1) 616 617 CALL DGEMM('N','N', 618 & NBAS(ISYMD),NUMV*NRHF(ISYMI),NVIR(ISYMB), 619 & ONE,XLAMDH(KOFF1),NTOTD,WORK(KOFF2),NTOTB, 620 & ZERO,WORK(KOFF3),NTOTD) 621 622 ENDDO 623 624C Read AO Cholesky vectors and store as L(g,d#J). 625C ----------------------------------------------- 626 627 DO IVEC = 1,NUMV 628 629 JVEC = JVEC1 + IVEC - 1 630 CALL CHO_READN(WORK(KCHOL),JVEC,1,WORK(KIND1),IDUM2, 631 & ISYCHO,IOPTR,WORK(KREAD),LREAD) 632 633 DO ISYMD = 1,NSYM 634 635 ISYMG = MULD2H(ISYMD,ISYCHO) 636 637 LENGD = NBAS(ISYMG)*NBAS(ISYMD) 638 639 KOFF1 = KCHOL + IAODIS(ISYMG,ISYMD) 640 KOFF2 = KSCRL + IOFFL(ISYMD) + LENGD*(IVEC - 1) 641 642 CALL DCOPY(LENGD,WORK(KOFF1),1,WORK(KOFF2),1) 643 644 ENDDO 645 646 ENDDO 647 648C Calculate G term in AO basis. 649C ----------------------------- 650 651 DO ISYMD = 1,NSYM 652 653 ISYMG = MULD2H(ISYMD,ISYCHO) 654 ISYMI = MULD2H(ISYMD,ISYMDI) 655 656 KOFF1 = KSCRL + IOFFL(ISYMD) 657 KOFF2 = KSCRZ + IOFFZ(ISYMD) 658 KOFF3 = KSIGMA + IT1AO(ISYMG,ISYMI) 659 660 LENDJ = NBAS(ISYMD)*NUMV 661 662 NTOTG = MAX(NBAS(ISYMG),1) 663 NTODJ = MAX(LENDJ,1) 664 665 CALL DGEMM('N','N',NBAS(ISYMG),NRHF(ISYMI),LENDJ, 666 & ONE,WORK(KOFF1),NTOTG,WORK(KOFF2),NTODJ, 667 & ONE,WORK(KOFF3),NTOTG) 668 669 ENDDO 670 671 ENDDO 672 673C Close occupied Cholesky vector file. 674C ------------------------------------ 675 676 CALL CHO_MOP(1,4,ISYCHO,LUCHJI,1,1) 677 678C Close Y intermediate file. 679C -------------------------- 680 681 CALL CC_CYIOP(1,ISYCHO,1) 682 683C Escape point for empty symmetry. 684C -------------------------------- 685 686 1000 CONTINUE 687 688 ENDDO 689 690C Transform the G term to MO basis. 691C --------------------------------- 692 693 DO ISYMI = 1,NSYM 694 695 ISYMA = MULD2H(ISYMI,ISYMY) 696 ISYMG = ISYMA 697 698 NTOTG = MAX(NBAS(ISYMG),1) 699 NTOTA = MAX(NVIR(ISYMA),1) 700 701 KOFF1 = IGLMVI(ISYMG,ISYMA) + 1 702 KOFF2 = KSIGMA + IT1AO(ISYMG,ISYMI) 703 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 704 705 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 706 & ONE,XLAMDP(KOFF1),NTOTG,WORK(KOFF2),NTOTG, 707 & ONE,SIGMA1(KOFF3),NTOTA) 708 709 ENDDO 710 711 RETURN 712 END 713C /* Deck cc_choltrgijh */ 714 SUBROUTINE CC_CHOLTRGIJH(XLAMDP,XLAMDH,SIGMA1,X1AM,CIM,WORK,LWORK, 715 & ISYMX) 716C 717C Thomas Bondo Pedersen, February 2003. 718C 719C Purpose: Calculate the GIJ and H terms for left-hand Jacobian 720C transformations. 721C 722C Formula: 723C 724C SIGMA1(ai) = SIGMA1(ai) 725C - sum(Jj) Y(aj,J) * L(ij,J) 726C + sum(g) LamH(g,a) 727C * sum(Jd) L(gd,J) 728C * [ sum(b) LamP(d,b) * {Y(bi,J) - sum(j) X(bj)*L(ij,J)} 729C - sum(j) LamP(d,j) * sum(b) L(ib,J)*C(bj) 730C + 2 LamP(d,i) * sum(bj) {L(jb,J)*C(bj) + L(bj,J)*X(bj)} ] 731C 732C where g,d refer to AO basis, and L(gd,J) = L(dg,J) is used. 733C 734#include "implicit.h" 735 DIMENSION XLAMDP(*), XLAMDH(*), SIGMA1(*), X1AM(*), CIM(*) 736 DIMENSION WORK(LWORK) 737#include "maxorb.h" 738#include "ccdeco.h" 739#include "ccorb.h" 740#include "ccsdsym.h" 741#include "priunit.h" 742 743 INTEGER IOFFL(8), IOFFZ(8), IOFFY(8) 744 745 CHARACTER*13 SECNAM 746 PARAMETER (SECNAM = 'CC_CHOLTRGIJH') 747 748 PARAMETER (IOPTR = 2) 749 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 750 751c snrm = dsqrt(ddot(nt1am(isymx),sigma1,1,sigma1,1)) 752c xnrm = dsqrt(ddot(nt1am(isymx),x1am,1,x1am,1)) 753c cnrm = dsqrt(ddot(nt1am(isymx),cim,1,cim,1)) 754c write(lupri,'(A,A,1P,D22.15)') 755c & SECNAM,': norm of sigma1: ',snrm 756c write(lupri,'(A,A,1P,D22.15)') 757c & SECNAM,': norm of x1am : ',xnrm 758c write(lupri,'(A,A,1P,D22.15)') 759c & SECNAM,': norm of cim : ',cnrm 760 761C Set symmetries. 762C --------------- 763 764 ISYCIM = ISYMX 765 ISYMY = ISYMX 766 767C Read reduce index array. 768C ------------------------ 769 770 KIND1 = 1 771 CALL CC_GETIND1(WORK(KIND1),LWORK,LIND1) 772 KEND0 = KIND1 + LIND1 773 LWRK0 = LWORK - KEND0 + 1 774 775 IF (LWRK0 .LT. 0) THEN 776 WRITE(LUPRI,'(//,5X,A,A,A)') 777 & 'Insufficient memory in ',SECNAM,' - allocation: index' 778 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 779 & 'Need (more than): ',KEND0-1, 780 & 'Available : ',LWORK 781 CALL QUIT('Insufficient memory in '//SECNAM) 782 ENDIF 783 784C Allocation. 785C ----------- 786 787 KSIGMA = KEND0 788 KEND1 = KSIGMA + NT1AO(ISYMY) 789 LWRK1 = LWORK - KEND1 + 1 790 791 IF (LWRK1 .LT. 0) THEN 792 WRITE(LUPRI,'(//,5X,A,A)') 793 & 'Insufficient memory in ',SECNAM 794 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 795 & 'Need (more than): ',KEND1-1, 796 & 'Available : ',LWORK 797 CALL QUIT('Insufficient memory in '//SECNAM) 798 ENDIF 799 800C Initialize AO GIJ term. 801C ----------------------- 802 803 CALL DZERO(WORK(KSIGMA),NT1AO(ISYMY)) 804 805C Start Cholesky symmetry loop. 806C ----------------------------- 807 808 DO ISYCHO = 1,NSYM 809 810 IF (NUMCHO(ISYCHO) .LE. 0) GO TO 1000 811 812 ISYMBI = MULD2H(ISYCHO,ISYMY) 813 ISYMAJ = ISYMBI 814 ISYMDI = ISYMBI 815 ISYMKI = ISYMBI 816 817C Allocation for one AO Cholesky vector. 818C -------------------------------------- 819 820 LREAD = MEMRD(1,ISYCHO,IOPTR) 821 MAXKI = 0 822 DO ISYMI = 1,NSYM 823 ISYMK = MULD2H(ISYMI,ISYMKI) 824 NTEST = NRHF(ISYMK)*NRHF(ISYMI) 825 MAXKI = MAX(MAXKI,NTEST) 826 ENDDO 827 LSCR = MAX(LREAD,MAXKI) 828 829 KCHOL = KEND1 830 KREAD = KCHOL + N2BST(ISYCHO) 831 KEND2 = KREAD + LSCR 832 LWRK2 = LWORK - KEND2 + 1 833 834 IF (LWRK2 .LE. 0) THEN 835 WRITE(LUPRI,'(//,5X,A,A)') 836 & 'Insufficient memory in ',SECNAM 837 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 838 & 'Need (more than): ',KEND2-1, 839 & 'Available : ',LWORK 840 CALL QUIT('Insufficient memory in '//SECNAM) 841 ENDIF 842 843C Set up batch. 844C SCRL: L(ij,#J),Y(b#J,i),L(ck,#J),L(ic,#J),L(g,d#J). 845C SCRZ: Y(bi,#J),Z(d#J,i). 846C --------------------------------------------------- 847 848 LENL = MAX(NMATIJ(ISYCHO),NT1AM(ISYMBI),NT1AM(ISYCHO), 849 & N2BST(ISYCHO)) 850 LENZ = MAX(NT1AM(ISYMBI),NT1AO(ISYMDI)) 851 852 MINMEM = LENL + LENZ 853 IF (ISYCHO .EQ. ISYMX) MINMEM = MINMEM + 1 854 IF (MINMEM .GT. 0) THEN 855 NVEC = MIN(LWRK2/MINMEM,NUMCHO(ISYCHO)) 856 ELSE IF (MINMEM .EQ. 0) THEN 857 GO TO 1000 858 ELSE 859 WRITE(LUPRI,'(//,5X,A,A,A,I10,/)') 860 & 'MINMEM is negative in ',SECNAM,': ',MINMEM 861 CALL QUIT('Error in '//SECNAM) 862 ENDIF 863 864 IF (NVEC .LE. 0) THEN 865 WRITE(LUPRI,'(//,5X,A,A)') 866 & 'Insufficient memory for Cholesky batch in ',SECNAM 867 WRITE(LUPRI,'(5X,A,I2,A)') 868 & '(Cholesky symmetry:',ISYCHO,')' 869 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10,/)') 870 & 'Minimum memory required for batch: ',MINMEM, 871 & 'Memory available for batch : ',LWRK2, 872 & 'Total memory available : ',LWORK 873 CALL QUIT('Insufficient memory in '//SECNAM) 874 ENDIF 875 876 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 877 878C Open Y intermediate file. 879C ------------------------- 880 881 CALL CC_CYIOP(-1,ISYCHO,-1) 882 883C Open MO Cholesky vector files. 884C ------------------------------ 885 886 CALL CHO_MOP(-1,1,ISYCHO,LUCHIA,1,1) 887 IF (ISYCHO .EQ. ISYMX) CALL CHO_MOP(-1,3,ISYCHO,LUCHAI,1,1) 888 CALL CHO_MOP(-1,4,ISYCHO,LUCHIJ,1,1) 889 890C Start batch loop. 891C ----------------- 892 893 DO IBATCH = 1,NBATCH 894 895 NUMV = NVEC 896 IF (IBATCH .EQ. NBATCH) THEN 897 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 898 ENDIF 899 JVEC1 = NVEC*(IBATCH - 1) + 1 900 901C Set up local index arrays. 902C -------------------------- 903 904 ICOUNL = 0 905 ICOUNZ = 0 906 ICOUNY = 0 907 DO ISYMD = 1,NSYM 908 ISYMI = MULD2H(ISYMD,ISYMDI) 909 ISYMG = MULD2H(ISYMD,ISYCHO) 910 ISYMB = ISYMD 911 IOFFL(ISYMD) = ICOUNL 912 IOFFZ(ISYMD) = ICOUNZ 913 IOFFY(ISYMB) = ICOUNY 914 LENDJ = NBAS(ISYMD)*NUMV 915 LENBJ = NVIR(ISYMB)*NUMV 916 ICOUNL = ICOUNL + NBAS(ISYMG)*LENDJ 917 ICOUNZ = ICOUNZ + LENDJ*NRHF(ISYMI) 918 ICOUNY = ICOUNY + LENBJ*NRHF(ISYMI) 919 ENDDO 920 921C Complete allocation. 922C -------------------- 923 924 KSCRL = KEND2 925 KSCRZ = KSCRL + LENL*NUMV 926 KUVEC = KSCRZ + LENZ*NUMV 927 928C Read Y(bi,#J) in place of Z. 929C ---------------------------- 930 931 CALL CC_CYIRDF(WORK(KSCRZ),NUMV,JVEC1,ISYCHO,ISYMY,-1) 932 933c ynrm1 = dsqrt(ddot(nt1am(isymbi)*numv,work(kscrz),1, 934c & work(kscrz),1)) 935c write(lupri,*) '...ynrm at read : ',ynrm1 936 937C Read L(ij,#J) in place of L(g,d#J). 938C ----------------------------------- 939 940 CALL CHO_MOREAD(WORK(KSCRL),NMATIJ(ISYCHO),NUMV,JVEC1, 941 & LUCHIJ) 942 943C Calculate H term. 944C ----------------- 945 946 DO IVEC = 1,NUMV 947 DO ISYMJ = 1,NSYM 948 949 ISYMI = MULD2H(ISYMJ,ISYCHO) 950 ISYMA = MULD2H(ISYMJ,ISYMAJ) 951 952 KOFF1 = KSCRZ + NT1AM(ISYMAJ)*(IVEC - 1) 953 & + IT1AM(ISYMA,ISYMJ) 954 KOFF2 = KSCRL + NMATIJ(ISYCHO)*(IVEC - 1) 955 & + IMATIJ(ISYMI,ISYMJ) 956 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 957 958 NTOTA = MAX(NVIR(ISYMA),1) 959 NTOTI = MAX(NRHF(ISYMI),1) 960 961 CALL DGEMM('N','T', 962 & NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 963 & XMONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOTI, 964 & ONE,SIGMA1(KOFF3),NTOTA) 965 966 ENDDO 967 ENDDO 968 969C Subtract sum(j) X(bj) * L(ij,#J) from Y(bi,#J). 970C ----------------------------------------------- 971 972 DO IVEC = 1,NUMV 973 DO ISYMJ = 1,NSYM 974 975 ISYMI = MULD2H(ISYMJ,ISYCHO) 976 ISYMB = MULD2H(ISYMJ,ISYMX) 977 978 KOFF1 = IT1AM(ISYMB,ISYMJ) + 1 979 KOFF2 = KSCRL + NMATIJ(ISYCHO)*(IVEC - 1) 980 & + IMATIJ(ISYMI,ISYMJ) 981 KOFF3 = KSCRZ + NT1AM(ISYMBI)*(IVEC - 1) 982 & + IT1AM(ISYMB,ISYMI) 983 984 NTOTB = MAX(NVIR(ISYMB),1) 985 NTOTI = MAX(NRHF(ISYMI),1) 986 987 CALL DGEMM('N','T', 988 & NVIR(ISYMB),NRHF(ISYMI),NRHF(ISYMJ), 989 & XMONE,X1AM(KOFF1),NTOTB,WORK(KOFF2),NTOTI, 990 & ONE,WORK(KOFF3),NTOTB) 991 992 ENDDO 993 ENDDO 994 995c ynrm2 = dsqrt(ddot(nt1am(isymbi)*numv,work(kscrz),1, 996c & work(kscrz),1)) 997c write(lupri,*) '...ynrm after X-subtraction: ',ynrm2 998 999C Resort Y intermediate as Y(b#J,i), store in place of L(g,d#J). 1000C -------------------------------------------------------------- 1001 1002 DO IVEC = 1,NUMV 1003 DO ISYMI = 1,NSYM 1004 1005 ISYMB = MULD2H(ISYMI,ISYMBI) 1006 LENBJ = NVIR(ISYMB)*NUMV 1007 1008 DO I = 1,NRHF(ISYMI) 1009 1010 KOFF1 = KSCRZ + NT1AM(ISYMBI)*(IVEC - 1) 1011 & + IT1AM(ISYMB,ISYMI) + NVIR(ISYMB)*(I - 1) 1012 KOFF2 = KSCRL + IOFFY(ISYMB) + LENBJ*(I - 1) 1013 & + NVIR(ISYMB)*(IVEC - 1) 1014 1015 CALL DCOPY(NVIR(ISYMB),WORK(KOFF1),1,WORK(KOFF2),1) 1016 1017 ENDDO 1018 1019 ENDDO 1020 ENDDO 1021 1022c ynrm3 = dsqrt(ddot(nt1am(isymbi)*numv,work(kscrl),1, 1023c & work(kscrl),1)) 1024c write(lupri,*) '...ynrm after resort : ',ynrm3 1025 1026C Backtransform the Y intermediates to get Z(d#J,i). 1027C -------------------------------------------------- 1028 1029 DO ISYMB = 1,NSYM 1030 1031 ISYMI = MULD2H(ISYMB,ISYMBI) 1032 ISYMD = ISYMB 1033 1034 KOFF1 = IGLMVI(ISYMD,ISYMB) + 1 1035 KOFF2 = KSCRL + IOFFY(ISYMB) 1036 KOFF3 = KSCRZ + IOFFZ(ISYMD) 1037 1038 NTOTD = MAX(NBAS(ISYMD),1) 1039 NTOTB = MAX(NVIR(ISYMB),1) 1040 1041 CALL DGEMM('N','N', 1042 & NBAS(ISYMD),NUMV*NRHF(ISYMI),NVIR(ISYMB), 1043 & ONE,XLAMDP(KOFF1),NTOTD,WORK(KOFF2),NTOTB, 1044 & ZERO,WORK(KOFF3),NTOTD) 1045 1046 ENDDO 1047 1048 IF (ISYCHO .EQ. ISYMX) THEN 1049 1050C Read L(ck,#J) in place of L(g,d#J). 1051C ----------------------------------- 1052 1053 CALL CHO_MOREAD(WORK(KSCRL),NT1AM(ISYCHO),NUMV,JVEC1, 1054 & LUCHAI) 1055 1056C Calculate U(#J) = 2 * sum(ck) L(ck,#J) * X(ck). 1057C ----------------------------------------------- 1058 1059 NTOCK = MAX(NT1AM(ISYCHO),1) 1060 CALL DGEMV('T',NT1AM(ISYCHO),NUMV, 1061 & TWO,WORK(KSCRL),NTOCK,X1AM,1, 1062 & ZERO,WORK(KUVEC),1) 1063 1064 ENDIF 1065 1066C Read L(ic,#J) in place of L(g,d#J). 1067C ----------------------------------- 1068 1069 CALL CHO_MOREAD(WORK(KSCRL),NT1AM(ISYCHO),NUMV,JVEC1,LUCHIA) 1070 1071 IF (ISYCHO .EQ. ISYCIM) THEN 1072 1073C Calculate U(#J) = U(#J) + 2 * sum(ck) L(kc,#J) * CIM(ck). 1074C --------------------------------------------------------- 1075 1076 NTOCK = MAX(NT1AM(ISYCHO),1) 1077 CALL DGEMV('T',NT1AM(ISYCHO),NUMV, 1078 & TWO,WORK(KSCRL),NTOCK,CIM,1, 1079 & ONE,WORK(KUVEC),1) 1080 1081C Set up Z contribution: 1082C Z(d#J,i) = Z(d#J,i) + LamP(d,i) * U(#J) 1083C --------------------------------------- 1084 1085 DO ISYMD = 1,NSYM 1086 ISYMI = ISYMD 1087 DO I = 1,NRHF(ISYMI) 1088 KOFF1 = IGLMRH(ISYMD,ISYMI) 1089 & + NBAS(ISYMD)*(I - 1) + 1 1090 DO IVEC = 1,NUMV 1091 FAC = WORK(KUVEC+IVEC-1) 1092 KOFF2 = KSCRZ + IOFFZ(ISYMD) 1093 & + NBAS(ISYMD)*NUMV*(I - 1) 1094 & + NBAS(ISYMD)*(IVEC - 1) 1095 CALL DAXPY(NBAS(ISYMD),FAC,XLAMDP(KOFF1),1, 1096 & WORK(KOFF2),1) 1097 ENDDO 1098 ENDDO 1099 ENDDO 1100 1101 ENDIF 1102 1103C Transform: -sum(k) LamP(d,k) {sum(c) C(ck) * L(ic,#J)} 1104C and add result into Z(d#J,i). Store intermediate result 1105C in the scratch space for reading AO Cholesky vectors. 1106C ------------------------------------------------------- 1107 1108 DO IVEC = 1,NUMV 1109 DO ISYMI = 1,NSYM 1110 1111 ISYMC = MULD2H(ISYMI,ISYCHO) 1112 ISYMK = MULD2H(ISYMC,ISYCIM) 1113 1114 KOFF1 = IT1AM(ISYMC,ISYMK) + 1 1115 KOFF2 = KSCRL + NT1AM(ISYCHO)*(IVEC - 1) 1116 & + IT1AM(ISYMC,ISYMI) 1117 1118 NTOTC = MAX(NVIR(ISYMC),1) 1119 NTOTK = MAX(NRHF(ISYMK),1) 1120 1121 CALL DGEMM('T','N', 1122 & NRHF(ISYMK),NRHF(ISYMI),NVIR(ISYMC), 1123 & ONE,CIM(KOFF1),NTOTC,WORK(KOFF2),NTOTC, 1124 & ZERO,WORK(KREAD),NTOTK) 1125 1126 ISYMD = ISYMK 1127 KOFF1 = IGLMRH(ISYMD,ISYMK) + 1 1128 NTOTD = MAX(NBAS(ISYMD),1) 1129 1130 DO I = 1,NRHF(ISYMI) 1131 1132 KOFF2 = KREAD + NRHF(ISYMK)*(I - 1) 1133 KOFF3 = KSCRZ + IOFFZ(ISYMD) 1134 & + NBAS(ISYMD)*NUMV*(I - 1) 1135 & + NBAS(ISYMD)*(IVEC - 1) 1136 1137 CALL DGEMV('N',NBAS(ISYMD),NRHF(ISYMK), 1138 & XMONE,XLAMDP(KOFF1),NTOTD,WORK(KOFF2),1, 1139 & ONE,WORK(KOFF3),1) 1140 1141 ENDDO 1142 1143 ENDDO 1144 ENDDO 1145 1146C Read AO Cholesky vectors and store as L(g,d#J). 1147C ----------------------------------------------- 1148 1149c cnrm1 = 0.0d0 1150 1151 DO IVEC = 1,NUMV 1152 1153 JVEC = JVEC1 + IVEC - 1 1154 CALL CHO_READN(WORK(KCHOL),JVEC,1,WORK(KIND1),IDUM2, 1155 & ISYCHO,IOPTR,WORK(KREAD),LREAD) 1156 1157c cnrm1 = cnrm1 + ddot(n2bst(isycho),work(kchol),1,work(kchol),1) 1158 1159 DO ISYMD = 1,NSYM 1160 1161 ISYMG = MULD2H(ISYMD,ISYCHO) 1162 1163 LENGD = NBAS(ISYMG)*NBAS(ISYMD) 1164 1165 KOFF1 = KCHOL + IAODIS(ISYMG,ISYMD) 1166 KOFF2 = KSCRL + IOFFL(ISYMD) + LENGD*(IVEC - 1) 1167 1168 CALL DCOPY(LENGD,WORK(KOFF1),1,WORK(KOFF2),1) 1169 1170 ENDDO 1171 1172 ENDDO 1173 1174c cnrm1 = dsqrt(cnrm1) 1175c cnrm2 = dsqrt(ddot(n2bst(isycho)*numv,work(kscrl),1, 1176c & work(kscrl),1)) 1177c write(lupri,*) '...Norm of Cholesky, read: ',cnrm1 1178c write(lupri,*) '...Norm of Cholesky, sort: ',cnrm2 1179 1180C Calculate GIJ term in AO basis. 1181C ------------------------------- 1182 1183 DO ISYMD = 1,NSYM 1184 1185 ISYMG = MULD2H(ISYMD,ISYCHO) 1186 ISYMI = MULD2H(ISYMD,ISYMDI) 1187 1188 KOFF1 = KSCRL + IOFFL(ISYMD) 1189 KOFF2 = KSCRZ + IOFFZ(ISYMD) 1190 KOFF3 = KSIGMA + IT1AO(ISYMG,ISYMI) 1191 1192 LENDJ = NBAS(ISYMD)*NUMV 1193 1194 NTOTG = MAX(NBAS(ISYMG),1) 1195 NTODJ = MAX(LENDJ,1) 1196 1197 CALL DGEMM('N','N',NBAS(ISYMG),NRHF(ISYMI),LENDJ, 1198 & ONE,WORK(KOFF1),NTOTG,WORK(KOFF2),NTODJ, 1199 & ONE,WORK(KOFF3),NTOTG) 1200 1201 ENDDO 1202 1203 ENDDO 1204 1205C Close MO Cholesky vector files. 1206C ------------------------------- 1207 1208 CALL CHO_MOP(1,4,ISYCHO,LUCHIJ,1,1) 1209 IF (ISYCHO .EQ. ISYMX) CALL CHO_MOP(1,3,ISYCHO,LUCHAI,1,1) 1210 CALL CHO_MOP(1,1,ISYCHO,LUCHIA,1,1) 1211 1212C Close Y intermediate file. 1213C -------------------------- 1214 1215 CALL CC_CYIOP(1,ISYCHO,-1) 1216 1217C Escape point for empty symmetry. 1218C -------------------------------- 1219 1220 1000 CONTINUE 1221 1222 ENDDO 1223 1224C Transform the GIJ term to MO basis. 1225C ----------------------------------- 1226 1227 DO ISYMI = 1,NSYM 1228 1229 ISYMA = MULD2H(ISYMI,ISYMY) 1230 ISYMG = ISYMA 1231 1232 NTOTG = MAX(NBAS(ISYMG),1) 1233 NTOTA = MAX(NVIR(ISYMA),1) 1234 1235 KOFF1 = IGLMVI(ISYMG,ISYMA) + 1 1236 KOFF2 = KSIGMA + IT1AO(ISYMG,ISYMI) 1237 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1238 1239 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 1240 & ONE,XLAMDH(KOFF1),NTOTG,WORK(KOFF2),NTOTG, 1241 & ONE,SIGMA1(KOFF3),NTOTA) 1242 1243 ENDDO 1244 1245 RETURN 1246 END 1247C /* Deck cc_choatr0 */ 1248 SUBROUTINE CC_CHOATR0(SIGMA1,X1AM,WORK,LWORK,ISYMX,NUMX,ISIDE) 1249C 1250C Thomas Bondo Pedersen, February 2003. 1251C 1252C Purpose: Calculate contributions from global (tot. sym.) E intermediates 1253C to left-hand (ISIDE=-1) or right-hand (ISIDE=1) Jacobian 1254C tranformations. 1255C 1256C ISIDE = -1: 1257C 1258C SIGMA1(ai) = sum(b) E(ba) * X(bi) - sum(j) X(aj) * E(ij) 1259C 1260C ISIDE = 1: 1261C 1262C SIGMA1(ai) = sum(b) E(ab) * X(bi) - sum(j) X(aj) * E(ji) 1263C 1264C NOTICE: the content of SIGMA1 is overwritten! 1265C 1266#include "implicit.h" 1267 DIMENSION SIGMA1(*), X1AM(*), WORK(LWORK) 1268#include "ccorb.h" 1269#include "ccsdsym.h" 1270#include "priunit.h" 1271 1272 LOGICAL LOCDBG 1273 PARAMETER (LOCDBG = .FALSE.) 1274 1275 CHARACTER*10 SECNAM 1276 PARAMETER (SECNAM = 'CC_CHOATR0') 1277 1278 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 1279 1280C Read the E intermediates. 1281C ------------------------- 1282 1283 KEIJ = 1 1284 KEAB = KEIJ + NMATIJ(1) 1285 KEND = KEAB + NMATAB(1) 1286 LWRK = LWORK - KEND + 1 1287 1288 IF (LWRK .LT. 0) THEN 1289 WRITE(LUPRI,'(//,5X,A,A)') 1290 & 'Insufficient memory in ',SECNAM 1291 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1292 & 'Need : ',KEND-1, 1293 & 'Available: ',LWORK 1294 CALL QUIT('Insufficient memory in '//SECNAM) 1295 ENDIF 1296 1297 CALL CHO_IMOP(-1,1,LUE,1) 1298 CALL CHO_MOREAD(WORK(KEIJ),NMATIJ(1),1,1,LUE) 1299 CALL CHO_IMOP(1,1,LUE,1) 1300 1301 CALL CHO_IMOP(-1,2,LUE,1) 1302 CALL CHO_MOREAD(WORK(KEAB),NMATAB(1),1,1,LUE) 1303 CALL CHO_IMOP(1,2,LUE,1) 1304 1305 IF (LOCDBG) THEN 1306 EIJNRM = DSQRT(DDOT(NMATIJ(1),WORK(KEIJ),1,WORK(KEIJ),1)) 1307 EABNRM = DSQRT(DDOT(NMATAB(1),WORK(KEAB),1,WORK(KEAB),1)) 1308 WRITE(LUPRI,'(A,A,1P,D22.15)') 1309 & SECNAM,': norm of EIJ: ',EIJNRM 1310 WRITE(LUPRI,'(A,A,1P,D22.15)') 1311 & SECNAM,': norm of EAB: ',EABNRM 1312 ENDIF 1313 1314C Calculate contributions. 1315C ------------------------ 1316 1317 IF (ISIDE .EQ. -1) THEN 1318 CALL CC_CHOATR0M1(SIGMA1,X1AM,WORK(KEIJ),WORK(KEAB),ISYMX,1, 1319 & NUMX) 1320 ELSE IF (ISIDE .EQ. 1) THEN 1321 CALL CC_CHOATR0P1(SIGMA1,X1AM,WORK(KEIJ),WORK(KEAB),ISYMX,1, 1322 & NUMX) 1323 ELSE 1324 WRITE(LUPRI,'(//,5X,A,A,A,I10)') 1325 & 'ISIDE must be -1 or +1 in ',SECNAM,': ISIDE = ',ISIDE 1326 CALL QUIT('Error in '//SECNAM) 1327 ENDIF 1328 1329 RETURN 1330 END 1331C /* Deck cc_choatr0m1 */ 1332 SUBROUTINE CC_CHOATR0M1(SIGMA1,X1AM,EIJ,EAB,ISYMX,ISYME,NUMX) 1333C 1334C Thomas Bondo Pedersen, February 2003. 1335C 1336C Purpose: Calculate E-intermediate contributions, left-hand side style: 1337C 1338C SIGMA1(ai) = sum(b) E(ba) * X(bi) - sum(j) X(aj) * E(ij) 1339C 1340C NOTICE: the content of SIGMA1 is overwritten! 1341C 1342#include "implicit.h" 1343 DIMENSION SIGMA1(*), X1AM(*), EIJ(*), EAB(*) 1344#include "ccorb.h" 1345#include "ccsdsym.h" 1346 1347 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 1348 1349C Get result symmetry. 1350C -------------------- 1351 1352 ISYRES = MULD2H(ISYMX,ISYME) 1353 1354 DO IX = 1,NUMX 1355 1356 KOFX1 = NT1AM(ISYMX)*(IX - 1) + 1 1357 KOFS1 = NT1AM(ISYRES)*(IX - 1) + 1 1358 1359C Virtual part. 1360C ------------- 1361 1362 DO ISYMI = 1,NSYM 1363 1364 ISYMB = MULD2H(ISYMI,ISYMX) 1365 ISYMA = MULD2H(ISYMB,ISYME) 1366 1367 NTOTA = MAX(NVIR(ISYMA),1) 1368 NTOTB = MAX(NVIR(ISYMB),1) 1369 1370 KOFFE = IMATAB(ISYMB,ISYMA) + 1 1371 KOFFX = KOFX1 + IT1AM(ISYMB,ISYMI) 1372 KOFFS = KOFS1 + IT1AM(ISYMA,ISYMI) 1373 1374 CALL DGEMM('T','N', 1375 & NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1376 & ONE,EAB(KOFFE),NTOTB,X1AM(KOFFX),NTOTB, 1377 & ZERO,SIGMA1(KOFFS),NTOTA) 1378 1379 ENDDO 1380 1381C Occupied part. 1382C -------------- 1383 1384 DO ISYMJ = 1,NSYM 1385 1386 ISYMA = MULD2H(ISYMJ,ISYMX) 1387 ISYMI = MULD2H(ISYMJ,ISYME) 1388 1389 NTOTA = MAX(NVIR(ISYMA),1) 1390 NTOTI = MAX(NRHF(ISYMI),1) 1391 1392 KOFFX = KOFX1 + IT1AM(ISYMA,ISYMJ) 1393 KOFFE = IMATIJ(ISYMI,ISYMJ) + 1 1394 KOFFS = KOFS1 + IT1AM(ISYMA,ISYMI) 1395 1396 CALL DGEMM('N','T', 1397 & NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1398 & XMONE,X1AM(KOFFX),NTOTA,EIJ(KOFFE),NTOTI, 1399 & ONE,SIGMA1(KOFFS),NTOTA) 1400 1401 ENDDO 1402 1403 ENDDO 1404 1405 RETURN 1406 END 1407C /* Deck cc_choatr0p1 */ 1408 SUBROUTINE CC_CHOATR0P1(SIGMA1,X1AM,EIJ,EAB,ISYMX,ISYME,NUMX) 1409C 1410C Thomas Bondo Pedersen, February 2003. 1411C 1412C Purpose: Calculate E-intermediate contributions, right-hand side style: 1413C 1414C SIGMA1(ai) = sum(b) E(ab) * X(bi) - sum(j) X(aj) * E(ji) 1415C 1416C NOTICE: the content of SIGMA1 is overwritten! 1417C 1418#include "implicit.h" 1419 DIMENSION SIGMA1(*), X1AM(*), EIJ(*), EAB(*) 1420#include "ccorb.h" 1421#include "ccsdsym.h" 1422 1423 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 1424 1425C Get result symmetry. 1426C -------------------- 1427 1428 ISYRES = MULD2H(ISYMX,ISYME) 1429 1430 DO IX = 1,NUMX 1431 1432 KOFX1 = NT1AM(ISYMX)*(IX - 1) + 1 1433 KOFS1 = NT1AM(ISYRES)*(IX - 1) + 1 1434 1435C Virtual part. 1436C ------------- 1437 1438 DO ISYMI = 1,NSYM 1439 1440 ISYMB = MULD2H(ISYMI,ISYMX) 1441 ISYMA = MULD2H(ISYMB,ISYME) 1442 1443 NTOTA = MAX(NVIR(ISYMA),1) 1444 NTOTB = MAX(NVIR(ISYMB),1) 1445 1446 KOFFE = IMATAB(ISYMA,ISYMB) + 1 1447 KOFFX = KOFX1 + IT1AM(ISYMB,ISYMI) 1448 KOFFS = KOFS1 + IT1AM(ISYMA,ISYMI) 1449 1450 CALL DGEMM('N','N', 1451 & NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMB), 1452 & ONE,EAB(KOFFE),NTOTA,X1AM(KOFFX),NTOTB, 1453 & ZERO,SIGMA1(KOFFS),NTOTA) 1454 1455 ENDDO 1456 1457C Occupied part. 1458C -------------- 1459 1460 DO ISYMI = 1,NSYM 1461 1462 ISYMJ = MULD2H(ISYMI,ISYME) 1463 ISYMA = MULD2H(ISYMJ,ISYMX) 1464 1465 NTOTA = MAX(NVIR(ISYMA),1) 1466 NTOTJ = MAX(NRHF(ISYMJ),1) 1467 1468 KOFFX = KOFX1 + IT1AM(ISYMA,ISYMJ) 1469 KOFFE = IMATIJ(ISYMJ,ISYMI) + 1 1470 KOFFS = KOFS1 + IT1AM(ISYMA,ISYMI) 1471 1472 CALL DGEMM('N','N', 1473 & NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMJ), 1474 & XMONE,X1AM(KOFFX),NTOTA,EIJ(KOFFE),NTOTJ, 1475 & ONE,SIGMA1(KOFFS),NTOTA) 1476 1477 ENDDO 1478 1479 ENDDO 1480 1481 RETURN 1482 END 1483C /* Deck cc_chortri2j */ 1484 SUBROUTINE CC_CHORTRI2J(FOCKD,XLAMDP,XLAMDH,SIGMA1,X1AM, 1485 & WORK,LWORK,ISYMX,NUMX) 1486C 1487C Thomas Bondo Pedersen, February 2003. 1488C 1489C Purpose: Calculate the I2 and J terms for right-hand Jacobian 1490C transformation. 1491C 1492C The J term: 1493C 1494C SIGMA1(ai) = SIGMA1(ai) + sum(ck) [2(ia|kc)-(ac|ka)] * X1AM(ck) 1495C 1496C The I2 term: 1497C 1498C SIGMA1(ai) = SIGMA1(ai) 1499C + sum(bj) [2t(ai,bj)-t(aj,bi)] 1500C * sum(ck) [2(ai|kc)-(ac|ki)] * X1AM(ck) 1501C 1502C The calculation is done through construction of the Fock-matrix type 1503C intermediate 1504C 1505C W(al,be) = sum(ck) [2(al be|kc)-(al c|k be)] * X1AM(ck) 1506C 1507C and, for the I2 term, through Cholesky decomposed 0th order doubles 1508C amplitudes (as C intermediates in left-hand transformation). The 1509C canonical orbital energies (FOCKD) are needed for this decomposition, 1510C if it has not already been done. 1511C 1512#include "implicit.h" 1513 DIMENSION FOCKD(*), XLAMDP(*), XLAMDH(*), SIGMA1(*), X1AM(*) 1514 DIMENSION WORK(LWORK) 1515#include "priunit.h" 1516#include "ccorb.h" 1517#include "ccsdsym.h" 1518 1519 CHARACTER*12 SECNAM 1520 PARAMETER (SECNAM = 'CC_CHORTRI2J') 1521 1522C Return if nothing to do. 1523C ------------------------ 1524 1525 IF (NUMX .LE. 0) RETURN 1526 1527C Allocation. 1528C ----------- 1529 1530 KLAMD2 = 1 1531 KWMAT = KLAMD2 + MAX(NGLMRH(ISYMX),NT1AM(ISYMX))*NUMX 1532 KEND1 = KWMAT + N2BST(ISYMX)*NUMX 1533 LWRK1 = LWORK - KEND1 + 1 1534 1535 KFJBP = KLAMD2 1536 1537 IF (LWRK1 .LT. 0) THEN 1538 WRITE(LUPRI,'(//,5X,A,A,A)') 1539 & 'Insufficient memory in ',SECNAM,' - allocation: W' 1540 WRITE(LUPRI,'(5X,A,I10,A,I1,A)') 1541 & 'Number of trial vectors: ',NUMX,' (symmetry ',ISYMX,')' 1542 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1543 & 'Need : ',KEND1-1, 1544 & 'Available: ',LWORK 1545 CALL QUIT('Insufficient memory in '//SECNAM) 1546 ENDIF 1547 1548C Backtransform X1AM to obtain perturbed Lambda-hole matrices. 1549C ------------------------------------------------------------ 1550 1551 CALL CC_CHOBCKTR(XLAMDH,X1AM,WORK(KLAMD2),ISYMX,NUMX) 1552 1553C Calculate perturbed Fock matrices (W intermediates). 1554C ---------------------------------------------------- 1555 1556 CALL DZERO(WORK(KWMAT),N2BST(ISYMX)*NUMX) 1557 CALL CC_CHOFCKP(XLAMDP,WORK(KLAMD2),WORK(KWMAT),WORK(KEND1),LWRK1, 1558 & 1,ISYMX,NUMX) 1559 1560C Calculate J terms. 1561C ------------------ 1562 1563 CALL CC_CHORTRJ(XLAMDP,XLAMDH,SIGMA1,WORK(KWMAT),WORK(KEND1), 1564 & LWRK1,1,1,ISYMX,NUMX) 1565 1566C Calculate perturbed F(jb). 1567C -------------------------- 1568 1569 CALL DZERO(WORK(KFJBP),NT1AM(ISYMX)*NUMX) 1570 CALL CC_CHORTRFJB(XLAMDP,XLAMDH,WORK(KFJBP),WORK(KWMAT), 1571 & WORK(KEND1),LWRK1,1,1,ISYMX,NUMX) 1572 1573C Calculate I2 terms. 1574C ------------------- 1575 1576 KEND1 = KFJBP + NT1AM(ISYMX)*NUMX 1577 LWRK1 = LWORK - KEND1 + 1 1578 CALL CC_CHOCIM1(FOCKD,WORK(KFJBP),SIGMA1,WORK(KEND1),LWRK1,ISYMX, 1579 & NUMX) 1580 1581 RETURN 1582 END 1583C /* Deck cc_chortrfjb */ 1584 SUBROUTINE CC_CHORTRFJB(XLAMDP,XLAMDH,SIGMA1,WMAT,WORK,LWORK, 1585 & ISYMP,ISYMH,ISYMW,NUMW) 1586C 1587C Thomas Bondo Pedersen, February 2003. 1588C 1589C Purpose: Transform NUMW AO matrices to MO basis as 1590C 1591C SIGMA1(ia) = SIGMA1(ia) 1592C + sum(g,d) XLAMDH(g,a) * WMAT(d,g) * XLAMDP(d,i) 1593C 1594C [with SIGMA1(ia) stored as ai] 1595C 1596#include "implicit.h" 1597 DIMENSION XLAMDP(*), XLAMDH(*), SIGMA1(*), WMAT(*), WORK(LWORK) 1598#include "ccorb.h" 1599#include "ccsdsym.h" 1600#include "priunit.h" 1601 1602 CHARACTER*12 SECNAM 1603 PARAMETER (SECNAM = 'CC_CHORTRFJB') 1604 1605 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1606 1607C Allocation. 1608C ----------- 1609 1610 NEED = 0 1611 DO ISYMG = 1,NSYM 1612 ISYMD = MULD2H(ISYMG,ISYMW) 1613 ISYMI = MULD2H(ISYMD,ISYMP) 1614 NEEDS = NBAS(ISYMG)*NRHF(ISYMI) 1615 NEED = MAX(NEED,NEEDS) 1616 ENDDO 1617 1618 IF (NEED .GT. LWORK) THEN 1619 WRITE(LUPRI,'(//,5X,A,A)') 1620 & 'Insufficient memory in ',SECNAM 1621 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1622 & 'Need : ',NEED, 1623 & 'Available: ',LWORK 1624 CALL QUIT('Insufficient memory in '//SECNAM) 1625 ENDIF 1626 1627C Calculate J terms. 1628C ------------------ 1629 1630 ISYMIM = MULD2H(ISYMP,ISYMW) 1631 ISYMS = MULD2H(ISYMH,ISYMIM) 1632 1633 DO IW = 1,NUMW 1634 DO ISYMG = 1,NSYM 1635 1636 ISYMD = MULD2H(ISYMG,ISYMW) 1637 ISYMI = MULD2H(ISYMD,ISYMP) 1638 ISYMA = MULD2H(ISYMG,ISYMH) 1639 1640 NTOTA = MAX(NVIR(ISYMA),1) 1641 NTOTG = MAX(NBAS(ISYMG),1) 1642 NTOTD = MAX(NBAS(ISYMD),1) 1643 1644 KOFFW = N2BST(ISYMW)*(IW - 1) + IAODIS(ISYMD,ISYMG) + 1 1645 KOFFP = IGLMRH(ISYMD,ISYMI) + 1 1646 KOFFH = IGLMVI(ISYMG,ISYMA) + 1 1647 KOFFS = NT1AM(ISYMS)*(IW - 1) + IT1AM(ISYMA,ISYMI) + 1 1648 1649 CALL DGEMM('T','N',NBAS(ISYMG),NRHF(ISYMI),NBAS(ISYMD), 1650 & ONE,WMAT(KOFFW),NTOTD,XLAMDP(KOFFP),NTOTD, 1651 & ZERO,WORK,NTOTG) 1652 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 1653 & ONE,XLAMDH(KOFFH),NTOTG,WORK,NTOTG, 1654 & ONE,SIGMA1(KOFFS),NTOTA) 1655 1656 ENDDO 1657 ENDDO 1658 1659 RETURN 1660 END 1661C /* Deck cc_chortrj */ 1662 SUBROUTINE CC_CHORTRJ(XLAMDP,XLAMDH,SIGMA1,WMAT,WORK,LWORK, 1663 & ISYMP,ISYMH,ISYMW,NUMW) 1664C 1665C Thomas Bondo Pedersen, February 2003. 1666C 1667C Purpose: Calculate J term for right-hand Jacobian transformation 1668C from the "perturbed" AO Fock matrix WMAT, 1669C 1670C SIGMA1(ai) = SIGMA1(ai) 1671C + sum(g,d) XLAMDP(g,a) * WMAT(g,d) * XLAMDH(d,i) 1672C 1673#include "implicit.h" 1674 DIMENSION XLAMDP(*), XLAMDH(*), SIGMA1(*), WMAT(*), WORK(LWORK) 1675#include "ccorb.h" 1676#include "ccsdsym.h" 1677#include "priunit.h" 1678 1679 CHARACTER*10 SECNAM 1680 PARAMETER (SECNAM = 'CC_CHORTRJ') 1681 1682 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 1683 1684C Allocation. 1685C ----------- 1686 1687 NEED = 0 1688 DO ISYMD = 1,NSYM 1689 ISYMG = MULD2H(ISYMD,ISYMW) 1690 ISYMI = MULD2H(ISYMD,ISYMH) 1691 NEEDS = NBAS(ISYMG)*NRHF(ISYMI) 1692 NEED = MAX(NEED,NEEDS) 1693 ENDDO 1694 1695 IF (NEED .GT. LWORK) THEN 1696 WRITE(LUPRI,'(//,5X,A,A)') 1697 & 'Insufficient memory in ',SECNAM 1698 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1699 & 'Need : ',NEED, 1700 & 'Available: ',LWORK 1701 CALL QUIT('Insufficient memory in '//SECNAM) 1702 ENDIF 1703 1704C Calculate J terms. 1705C ------------------ 1706 1707 ISYMIM = MULD2H(ISYMH,ISYMW) 1708 ISYMS = MULD2H(ISYMP,ISYMIM) 1709 1710 DO IW = 1,NUMW 1711 DO ISYMD = 1,NSYM 1712 1713 ISYMG = MULD2H(ISYMD,ISYMW) 1714 ISYMI = MULD2H(ISYMD,ISYMH) 1715 ISYMA = MULD2H(ISYMG,ISYMP) 1716 1717 NTOTA = MAX(NVIR(ISYMA),1) 1718 NTOTG = MAX(NBAS(ISYMG),1) 1719 NTOTD = MAX(NBAS(ISYMD),1) 1720 1721 KOFFW = N2BST(ISYMW)*(IW - 1) + IAODIS(ISYMG,ISYMD) + 1 1722 KOFFH = IGLMRH(ISYMD,ISYMI) + 1 1723 KOFFP = IGLMVI(ISYMG,ISYMA) + 1 1724 KOFFS = NT1AM(ISYMS)*(IW - 1) + IT1AM(ISYMA,ISYMI) + 1 1725 1726 CALL DGEMM('N','N',NBAS(ISYMG),NRHF(ISYMI),NBAS(ISYMD), 1727 & ONE,WMAT(KOFFW),NTOTG,XLAMDH(KOFFH),NTOTD, 1728 & ZERO,WORK,NTOTG) 1729 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 1730 & ONE,XLAMDP(KOFFP),NTOTG,WORK,NTOTG, 1731 & ONE,SIGMA1(KOFFS),NTOTA) 1732 1733 ENDDO 1734 ENDDO 1735 1736 RETURN 1737 END 1738C /* Deck cc_choaeffd */ 1739 SUBROUTINE CC_CHOAEFFD(WORK,LWORK,MULSOL,TBAR) 1740C 1741C Thomas Bondo Pedersen, February 2003. 1742C 1743C Purpose: Test RH and LH transformations by calculating 1744C the effective CC2 Jacobian via transformations with 1745C unit vectors. 1746C 1747C NOTE: this routine requires a *lot* of memory! 1748C 1749#include "implicit.h" 1750 DIMENSION WORK(LWORK), TBAR(*) 1751 LOGICAL MULSOL 1752#include "ccorb.h" 1753#include "ccsdsym.h" 1754#include "dccsdsym.h" 1755#include "ccsdinp.h" 1756#include "priunit.h" 1757 1758 CHARACTER*11 SECNAM 1759 PARAMETER (SECNAM = 'CC_CHOAEFFD') 1760 1761 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1762 1763 CALL AROUND(SECNAM//': Test of Jacobian tranformations') 1764 1765C Test memory. 1766C ------------ 1767 1768 XLWORK = ONE*LWORK 1769 XNEED = TWO*XT2SQ(1) 1770 IF (XNEED .GT. XLWORK) THEN 1771 WRITE(LUPRI,'(5X,A)') 1772 & 'Insufficient memory for test:' 1773 WRITE(LUPRI,'(5X,A,F15.1,A,/,5X,A,F15.1,A)') 1774 & 'Effective Jacobians alone require: ',XNEED,' words', 1775 & 'Total memory available is : ',XLWORK,' words' 1776 WRITE(LUPRI,'(5X,A)') 1777 & 'Test is skipped!' 1778 GO TO 1000 1779 ENDIF 1780 1781C Allocation for effective Jacobians. 1782C ----------------------------------- 1783 1784 KLEFT = 1 1785 KRIGHT = KLEFT + NT2SQ(1) 1786 KEND1 = KRIGHT + NT2SQ(1) 1787 LWRK1 = LWORK - KEND1 + 1 1788 1789 IF (LWRK1 .LT. 0) THEN 1790 WRITE(LUPRI,'(//,5X,A,A)') 1791 & 'Insufficient memory in ',SECNAM,' - allocation: Aeff' 1792 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1793 & 'Need (more than): ',KEND1-1, 1794 & 'Available : ',LWORK 1795 CALL QUIT('Insufficient memory in '//SECNAM) 1796 ENDIF 1797 1798C Reduce printing. 1799C ---------------- 1800 1801 IPRSAV = IPRINT 1802 IPRINT = -10 1803 1804C Left hand generation (row-wise). 1805C -------------------------------- 1806 1807 ISIDE = -1 1808 FREQ = ZERO 1809 CALL CC_CHOAEFFD2(WORK(KLEFT),WORK(KEND1),LWRK1,FREQ,ISIDE) 1810 1811c CALL AROUND('Aeff from left-hand transformation (row-wise)') 1812c CALL NOCC_PRT(WORK(KLEFT),1,'AIBJ') 1813 1814C Solve for multipliers. 1815C ---------------------- 1816 1817 IF (MULSOL) THEN 1818 1819 KETA = KEND1 1820 KSCR = KETA + NT1AM(1) 1821 KTBAR = KSCR + NT1AM(1) 1822 KEND = KTBAR + NT1AM(1) 1823 LWRK = LWORK - KEND + 1 1824 1825 IF (LWRK .LT. 0) THEN 1826 WRITE(LUPRI,'(//,5X,A,A)') 1827 & 'Insufficient memory for MULSOL in ',SECNAM 1828 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1829 & 'Need : ',KEND-1, 1830 & 'Available: ',LWORK 1831 CALL QUIT('Insufficient memory in '//SECNAM) 1832 ENDIF 1833 1834C Get right-hand side. 1835C -------------------- 1836 1837 CALL CC_CHOETA(WORK(KETA),WORK(KEND),LWRK) 1838 CALL DCOPY(NT1AM(1),WORK(KETA),1,WORK(KTBAR),1) 1839 CALL DSCAL(NT1AM(1),XMONE,WORK(KTBAR),1) 1840 1841C Copy out sym. 1 block. 1842C ---------------------- 1843 1844 KOFF1 = KLEFT + IT2SQ(1,1) 1845 LAIBJ = NT1AM(1)*NT1AM(1) 1846 CALL DCOPY(LAIBJ,WORK(KOFF1),1,WORK(KRIGHT),1) 1847 1848C Solve [Aeff]^T TBAR = -ETA. 1849C --------------------------- 1850 1851 INFO = -1 1852 CALL DGESOL(1,NT1AM(1),NT1AM(1),NT1AM(1), 1853 & WORK(KRIGHT),WORK(KTBAR),WORK(KSCR),INFO) 1854 1855 IF (INFO .NE. 0) THEN 1856 WRITE(LUPRI,'(//,5X,A,A,I10,/)') 1857 & SECNAM,': Error: DGESOL returned INFO = ',INFO 1858 CALL QUIT('DGESOL error in '//SECNAM) 1859 ENDIF 1860 1861C Check solution. 1862C --------------- 1863 1864 ETANRM = DSQRT(DDOT(NT1AM(1),WORK(KETA),1,WORK(KETA),1)) 1865 TBRNRM = DSQRT(DDOT(NT1AM(1),WORK(KTBAR),1,WORK(KTBAR),1)) 1866 1867 WRITE(LUPRI,'(//,5X,A,A,1P,D22.15)') 1868 & SECNAM,': norm of ETA : ',ETANRM 1869 WRITE(LUPRI,'(5X,A,A,1P,D22.15,/,5X,A)') 1870 & SECNAM,': norm of TBAR: ',TBRNRM, 1871 & '- going to calculate residual by trf. of solution vector...' 1872 1873 CALL CC_CHORSDTST(WORK(KTBAR),WORK(KEND),LWRK) 1874 1875 WRITE(LUPRI,'(5X,A)') 1876 & '- going to calculate trf. vector by DGEMV using Aeff...' 1877 1878 KOFF1 = KLEFT + IT2SQ(1,1) 1879 LAIBJ = NT1AM(1)*NT1AM(1) 1880 CALL DCOPY(LAIBJ,WORK(KOFF1),1,WORK(KRIGHT),1) 1881 CALL DGEMV('N',NT1AM(1),NT1AM(1),1.0D0,WORK(KRIGHT),NT1AM(1), 1882 & WORK(KTBAR),1,0.0D0,WORK(KSCR),1) 1883 TRFNRM = DSQRT(DDOT(NT1AM(1),WORK(KSCR),1,WORK(KSCR),1)) 1884 WRITE(LUPRI,'(5X,A,A,1P,D22.15,/)') 1885 & SECNAM,': norm of trf. vector from DGEMV calc.: ',TRFNRM 1886 1887 CALL DBGDIFAI(TBAR,WORK(KTBAR),TRGNRM,TSTNRM,ERRMIN,ERRMAX,1) 1888 WRITE(LUPRI,'(5X,A,A,1P,D22.15)') 1889 & SECNAM,': norm of TBAR from conv. calc.: ',TRGNRM 1890 WRITE(LUPRI,'(5X,A,A,1P,D22.15)') 1891 & SECNAM,': norm of TBAR from Chol. calc.: ',TBRNRM 1892 WRITE(LUPRI,'(5X,A,A,1P,D22.15)') 1893 & SECNAM,': Min. difference : ',ERRMIN 1894 WRITE(LUPRI,'(5X,A,A,1P,D22.15)') 1895 & SECNAM,': Max. difference : ',ERRMAX 1896 1897C Calculate density and properties. 1898C --------------------------------- 1899 1900 IF (NT2SQ(1) .GT. N2BST(1)) THEN 1901 1902 IPRINT = IPRSAV 1903 1904 KDEN = KRIGHT 1905 1906 WRITE(LUPRI,'(/,5X,A,A)') 1907 & SECNAM,': Calculating density matrix...' 1908 CALL CC_CHODEN(WORK(KTBAR),WORK(KDEN),WORK(KEND),LWRK) 1909 1910 WRITE(LUPRI,'(5X,A,A)') 1911 & SECNAM,': Calculating FOPs...' 1912 CALL CC_CHOFOP1(WORK(KDEN),WORK(KEND),LWRK,'CC2 ') 1913 1914 IPRSAV = IPRINT 1915 IPRINT = -10 1916 1917 ENDIF 1918 1919 ENDIF 1920 1921C Transpose left-hand result. 1922C --------------------------- 1923 1924 CALL DCOPY(NT2SQ(1),WORK(KLEFT),1,WORK(KRIGHT),1) 1925 DO ISYMBJ = 1,NSYM 1926 ISYMAI = ISYMBJ 1927 DO NBJ = 1,NT1AM(ISYMBJ) 1928 KOFF1 = KRIGHT + IT2SQ(ISYMBJ,ISYMAI) + NBJ - 1 1929 KOFF2 = KLEFT + IT2SQ(ISYMAI,ISYMBJ) 1930 & + NT1AM(ISYMAI)*(NBJ - 1) 1931 CALL DCOPY(NT1AM(ISYMAI),WORK(KOFF1),NT1AM(ISYMBJ), 1932 & WORK(KOFF2),1) 1933 ENDDO 1934 ENDDO 1935 1936c CALL AROUND('Aeff from left-hand transformation (column-wise)') 1937c CALL NOCC_PRT(WORK(KLEFT),1,'AIBJ') 1938 1939C Right hand generation (column-wise). 1940C ------------------------------------ 1941 1942 ISIDE = 1 1943 FREQ = ZERO 1944 CALL CC_CHOAEFFD2(WORK(KRIGHT),WORK(KEND1),LWRK1,FREQ,ISIDE) 1945 1946c CALL AROUND('Aeff from right-hand transformation (column-wise)') 1947c CALL NOCC_PRT(WORK(KRIGHT),1,'AIBJ') 1948 1949C Restore print level. 1950C -------------------- 1951 1952 IPRINT = IPRSAV 1953 1954C Compare matrices obtained from LH and RH transformations. 1955C --------------------------------------------------------- 1956 1957 CALL DBGDIFAIBJ(WORK(KLEFT),WORK(KRIGHT),ALNRM,ARNRM,ERRMIN, 1958 & ERRMAX,1,ISMNAI,ISMNBJ,MNAI,MNBJ,ISMXAI,ISMXBJ, 1959 & MXAI,MXBJ,TARMIN,TARMAX) 1960 1961 CALL DAXPY(NT2SQ(1),XMONE,WORK(KRIGHT),1,WORK(KLEFT),1) 1962 DIFNM2 = DDOT(NT2SQ(1),WORK(KLEFT),1,WORK(KLEFT),1) 1963 DIFNRM = DSQRT(DIFNM2) 1964 DIFRMS = DSQRT(DIFNM2/NT2SQ(1)) 1965 1966 CALL HEADER('Testing LH and RH eff. Jacobians',-1) 1967 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1968 & 'Norm of left-hand eff. Jacobian: ',ALNRM 1969 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1970 & 'Norm of right-hand eff. Jacobian: ',ARNRM 1971 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1972 & 'Difference of norms : ',ALNRM-ARNRM 1973 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1974 & 'Norm of difference : ',DIFNRM 1975 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1976 & 'Minimum absolute error : ',ERRMIN 1977 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1978 & 'Maximum absolute error : ',ERRMAX 1979 WRITE(LUPRI,'(5X,A,1P,D22.15,/)') 1980 & 'RMS error : ',DIFRMS 1981 WRITE(LUPRI,'(5X,A)') 1982 & 'Location, min. abs. err.:' 1983 WRITE(LUPRI,'(5X,A,I10,A,I1)') 1984 & ' BJ = ',MNBJ,' of sym. ',ISMNBJ 1985 WRITE(LUPRI,'(5X,A,I10,A,I1)') 1986 & ' AI = ',MNAI,' of sym. ',ISMNAI 1987 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1988 & ' Value of LH Aeff: ',TARMIN 1989 WRITE(LUPRI,'(5X,A)') 1990 & 'Location, max. abs. err.:' 1991 WRITE(LUPRI,'(5X,A,I10,A,I1)') 1992 & ' BJ = ',MXBJ,' of sym. ',ISMXBJ 1993 WRITE(LUPRI,'(5X,A,I10,A,I1)') 1994 & ' AI = ',MXAI,' of sym. ',ISMXAI 1995 WRITE(LUPRI,'(5X,A,1P,D22.15)') 1996 & ' Value of LH Aeff: ',TARMAX 1997 1998c CALL AROUND('L-R Aeff difference matrix') 1999c CALL NOCC_PRT(WORK(KLEFT),1,'AIBJ') 2000 2001 1000 CALL AROUND('End of '//SECNAM) 2002 2003 RETURN 2004 END 2005C /* Deck dbgdifaibj */ 2006 SUBROUTINE DBGDIFAIBJ(TARGET,TEST,TRGNRM,TSTNRM,ERRMIN,ERRMAX, 2007 & ISYM,ISMNAI,ISMNBJ,MNAI,MNBJ,ISMXAI,ISMXBJ, 2008 & MXAI,MXBJ,TARMIN,TARMAX) 2009C 2010C Thomas Bondo Pedersen, February 2003. 2011C 2012#include "implicit.h" 2013 DIMENSION TARGET(*), TEST(*) 2014#include "ccorb.h" 2015#include "ccsdsym.h" 2016 2017 TRGNRM = DSQRT(DDOT(NT2SQ(1),TARGET,1,TARGET,1)) 2018 TSTNRM = DSQRT(DDOT(NT2SQ(1),TEST,1,TEST,1)) 2019 ERRMIN = 1.0D10 2020 ERRMAX = -1.0D10 2021 2022 DO ISYMBJ = 1,NSYM 2023 ISYMAI = MULD2H(ISYMBJ,ISYM) 2024 DO NBJ = 1,NT1AM(ISYMBJ) 2025 DO NAI = 1,NT1AM(ISYMAI) 2026 NAIBJ = IT2SQ(ISYMAI,ISYMBJ) + NT1AM(ISYMAI)*(NBJ - 1) 2027 & + NAI 2028 DIFF = DABS(TARGET(NAIBJ) - TEST(NAIBJ)) 2029 IF (DIFF .LT. ERRMIN) THEN 2030 ISMNAI = ISYMAI 2031 ISMNBJ = ISYMBJ 2032 MNAI = NAI 2033 MNBJ = NBJ 2034 ERRMIN = DIFF 2035 TARMIN = TARGET(NAIBJ) 2036 ENDIF 2037 IF (DIFF .GT. ERRMAX) THEN 2038 ISMXAI = ISYMAI 2039 ISMXBJ = ISYMBJ 2040 MXAI = NAI 2041 MXBJ = NBJ 2042 ERRMAX = DIFF 2043 TARMAX = TARGET(NAIBJ) 2044 ENDIF 2045 ENDDO 2046 ENDDO 2047 ENDDO 2048 2049 RETURN 2050 END 2051C /* Deck cc_choaeffd2 */ 2052 SUBROUTINE CC_CHOAEFFD2(AEFF,WORK,LWORK,FREQ,ISIDE) 2053C 2054C Thomas Bondo Pedersen, February 2003. 2055C 2056C Purpose: Set up effective Jacobian by transformations 2057C with unit vectors. For ISIDE = -1, AEFF is stored 2058C transposed (i.e. row-wise). 2059C 2060#include "implicit.h" 2061 DIMENSION AEFF(*), WORK(LWORK) 2062#include "ccorb.h" 2063#include "ccsdsym.h" 2064#include "priunit.h" 2065 2066 CHARACTER*12 SECNAM 2067 PARAMETER (SECNAM = 'CC_CHOAEFFD2') 2068 2069 CHARACTER*5 ISITXT 2070 2071 PARAMETER (ONE = 1.0D0) 2072 2073 IF (ISIDE .EQ. -1) THEN 2074 ISITXT = 'left ' 2075 ELSE IF (ISIDE .EQ. 1) THEN 2076 ISITXT = 'right' 2077 ELSE 2078 WRITE(LUPRI,'(//,5X,A,A,A,/,5X,A,I10,/)') 2079 & 'ISIDE must be -1 or +1 in ',SECNAM,':', 2080 & 'ISIDE = ',ISIDE 2081 CALL QUIT('ISIDE error in '//SECNAM) 2082 ENDIF 2083 2084 WRITE(LUPRI,'(/,A,A,A,A,1P,D14.7,/)') 2085 & SECNAM,': eff. Jacobian transformations from ',ISITXT, 2086 & ' with FREQ=',FREQ 2087 2088 DO ISYMTR = 1,NSYM 2089 2090 KTRIAL = 1 2091 KEND1 = KTRIAL + NT1AM(ISYMTR) 2092 LWRK1 = LWORK - KEND1 + 1 2093 2094 IF (LWRK1 .LT. 0) THEN 2095 WRITE(LUPRI,'(//,5X,A,A)') 2096 & 'Insufficient memory in ',SECNAM,' - allocation: Trial' 2097 WRITE(LUPRI,'(5X,A,I1)') 2098 & 'Trial vector symmetry: ',ISYMTR 2099 WRITE(LUPRI,'(5X,A,I10)') 2100 & 'ISIDE = ',ISIDE 2101 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2102 & 'Need (more than): ',KEND1-1, 2103 & 'Available : ',LWORK 2104 CALL QUIT('Insufficient memory in '//SECNAM) 2105 ENDIF 2106 2107 DO IVEC = 1,NT1AM(ISYMTR) 2108 2109 WRITE(LUPRI,'(A,A,I10,A,I1,A,A)') 2110 & SECNAM,': transforming unit vector number ',IVEC, 2111 & ' of sym. ',ISYMTR,' from ',ISITXT 2112 2113 CALL DZERO(WORK(KTRIAL),NT1AM(ISYMTR)) 2114 WORK(KTRIAL+IVEC-1) = ONE 2115 2116 KOFFA = IT2SQ(ISYMTR,ISYMTR) + NT1AM(ISYMTR)*(IVEC - 1) + 1 2117 CALL CC_CHOATR(AEFF(KOFFA),WORK(KTRIAL),FREQ,WORK(KEND1), 2118 & LWRK1,ISYMTR,1,ISIDE) 2119 2120 ENDDO 2121 2122 ENDDO 2123 2124 RETURN 2125 END 2126