1C /* Deck cc_choecc2 */ 2 SUBROUTINE CC_CHOECC2(WORK,LWORK,ESCF,ECC2,RSPIM2) 3C 4C Thomas Bondo Pedersen, August 2002. 5C 6C Purpose: 7C Optimization of the CC2 wave function using Cholesky 8C decomposed two-electron integrals. The doubles amplitudes 9C are never stored in memory (unless memory actually allows 10C this). 11C 12C Overall structure is based on CCSD_ENERGY by Henrik Koch. 13C 14C IF (RSPIM2): Calculate the E intermediates for CCLR. 15C 16C WARNING: NOCCIT has not been implemented: some essential intermediates 17C are assumed to be generated in the course of the calculation, 18C such as transformed Cholesky vectors and F(ia). These will be 19C needed if this is a response calculation. Thus, at least one 20C iteration must be carried out. (Unless, of course, you opt for 21C a rewriting of the code....) 22C 23#include "implicit.h" 24 DIMENSION WORK(LWORK) 25 LOGICAL RSPIM2 26 COMMON /LUDIIS/ LUTDIS, LUSDIS 27#include "maxorb.h" 28#include "ccorb.h" 29#include "ccsdsym.h" 30#include "ccsdinp.h" 31#include "ccdeco.h" 32#include "priunit.h" 33#include "chocc2.h" 34#include "cclr.h" 35#include "dummy.h" 36#include "ccfro.h" 37 38 LOGICAL LCONVG,LXDIS 39 40 CHARACTER*10 MODEL 41 42 PARAMETER (LDIIST = 28, LDIISS = 29) 43 44 CHARACTER*10 SECNAM 45 PARAMETER (SECNAM = 'CC_CHOECC2') 46 47C Start timing. 48C ------------- 49 50 CALL QENTER(SECNAM) 51 CALL GETTIM(CSTR,WSTR) 52 TIMTOT = SECOND() 53 54C Check that this is a Cholesky run. 55C ---------------------------------- 56 57 IF (.NOT. CHOINT) THEN 58 WRITE(LUPRI,'(//,5X,A,A,A,/,5X,A,/)') 59 & 'FATAL ERROR IN ',SECNAM,':', 60 & '- integrals must be Cholesky decomposed!' 61 CALL QUIT('Fatal error in '//SECNAM) 62 ENDIF 63 64C Set MODEL (for CC_WRRSP and CC_RDRSP). 65C -------------------------------------- 66 67 MODEL = 'CC2 ' 68 69C Open DIIS file (CRAY only). 70C --------------------------- 71 72#if defined (SYS_CRAY) 73 CALL WOPEN('CC_DIIS',64,0,IERR) 74 IF (IERR .NE. 0) THEN 75 WRITE(LUPRI,'(//,5X,A,A,A,/,5X,A,I10,/)') 76 & 'Error opening file CC_DIIS in ',SECNAM,':', 77 & 'Subroutine WOPEN returned IERR = ',IERR 78 CALL QUIT('Error opening CC_DIIS in '//SECNAM) 79 ENDIF 80#endif 81 82C Print header. 83C ------------- 84 85 CALL HEADER('Output from '//SECNAM,-1) 86 WRITE(LUPRI,'(5X,A)') 87 & 'Cholesky based optimization of the CC2 wave function:' 88 IF (CHOT2) THEN 89 WRITE(LUPRI,'(5X,A)') 90 & 'The doubles amplitudes will be separately decomposed.' 91 ELSE IF (CHOMO2) THEN 92 WRITE(LUPRI,'(5X,A)') 93 & 'The (ai|bj) integral matrix will be separately decomposed.' 94 ENDIF 95 WRITE(LUPRI,'(5X,A,1P,D15.6,/,5X,A,1P,D15.6)') 96 & 'Threshold for energy : ',THRENR, 97 & 'Threshold for vector function: ',THRVEC 98 WRITE(LUPRI,'(5X,A,5X,I10,/)') 99 & 'Maximum number of iterations : ',MAXITE 100 101 CALL FLSHFO(LUPRI) 102 103C Calculate T1-independent MO transformations. 104C -------------------------------------------- 105 106 CALL CC2_TRF0(WORK,LWORK) 107 108C Allocation. 109C ----------- 110 111 KFOCKD = 1 112 KT1AM = KFOCKD + NORBTS 113 KOMEG1 = KT1AM + NT1AMX 114 KEND1 = KOMEG1 + NT1AMX 115 LWRK1 = LWORK - KEND1 + 1 116 117 IF (LWRK1 .LT. 0) THEN 118 WRITE(LUPRI,'(//,5X,A,A,A)') 119 & 'Insufficient memory in ',SECNAM,' - allocation: T1' 120 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 121 & 'Need (more than): ',KEND1-1, 122 & 'Available : ',LWORK 123 CALL QUIT('Insufficient memory in '//SECNAM) 124 ENDIF 125 126C Get SCF and orbital energies. 127C ----------------------------- 128 129 CALL CC2_GETSI(WORK(KFOCKD),WORK(KEND1),LWRK1,ESCF,POTNUC) 130 131C Set up initial amplitudes. 132C -------------------------- 133 134 CALL CC2_INITAM(WORK(KFOCKD),WORK(KT1AM),WORK(KOMEG1), 135 & WORK(KEND1),LWRK1,ECC2,IRST) 136 137 138C Start of iterative loop. 139C ------------------------ 140 141 ITER = 0 142 EN1 = ECC2 143 144 IF (IRST .EQ. 1) THEN 145 WRITE(LUPRI,'(1X,A,I3,A,F23.16)') 146 & 'Iter.',ITER,': Coupled cluster RSTAR energy : ',ECC2 147 CALL FLSHFO(LUPRI) 148 ENDIF 149 150 100 CONTINUE 151 152C Update iteration counter. 153C ------------------------- 154 155 ITER = ITER + 1 156 157 IF (ITER .GT. MAXITE) THEN 158 WRITE(LUPRI,'(//,5X,A,A,I4,A,/)') 159 & SECNAM,': CC2 wave function not converged in ',MAXITE, 160 & ' iterations.' 161 CALL QUIT(SECNAM//': CC2 wave function not converged') 162 ENDIF 163 164 IF (IPRINT .GT. 2) THEN 165 WRITE(LUPRI,'(/,7X,A,I3)') 'Iteration no.:',ITER 166 WRITE(LUPRI,'(7X,A,/)') '-----------------' 167 ENDIF 168 169C Calculate the CC2 vector function and energy. 170C --------------------------------------------- 171 172c lwsav = lwrk1 173c lwrk1 = 10000 174c write(lupri,*) 'calling rhs with lwork : ',LWRK1 175 176 CALL CC2_CHORHS(WORK(KFOCKD),WORK(KT1AM),WORK(KOMEG1), 177 & WORK(KEND1),LWRK1,T2NM,ECC2,ESCF) 178 EN2 = ECC2 179 180c lwrk1 = lwsav 181 182 IF (IPRINT .GE. 5) THEN 183 T1NM = DDOT(NT1AMX,WORK(KT1AM),1,WORK(KT1AM),1) 184 WRITE(LUPRI,'(7X,A,1P,D25.10)') 185 & 'Norm of t1am after ccvec: ',DSQRT(T1NM) 186 WRITE(LUPRI,'(7X,A,1P,D25.10)') 187 & 'Norm of t2am after ccvec: ',DSQRT(T2NM) 188 WRITE(LUPRI,'(7X,A,1P,D25.10)') 189 & 'Total norm of amplitudes : ',DSQRT(T1NM+T2NM) 190 ENDIF 191 192 O1NM = DSQRT(DDOT(NT1AMX,WORK(KOMEG1),1,WORK(KOMEG1),1)) 193 IF (IPRINT .GE. 3) THEN 194 WRITE(LUPRI,'(7X,A,1P,D25.10)') 195 & 'Norm of omega1 after ccvec: ',O1NM 196 IF (IPRINT .GE. 25) THEN 197 CALL HEADER('CC2 Vector Function',-1) 198 DO ISYM = 1,NSYM 199 WRITE(LUPRI,'(10X,A,I1,A,/)') 200 & 'Symmetry block ',ISYM,':' 201 NA = NVIR(ISYM) 202 NI = NRHF(ISYM) 203 IF ((NA.GT.0) .AND. (NI.GT.0)) THEN 204 KOFF = KOMEG1 + IT1AM(ISYM,ISYM) 205 CALL OUTPUT(WORK(KOFF),1,NA,1,NI,NA,NI,1,LUPRI) 206 WRITE(LUPRI,*) 207 ELSE 208 WRITE(LUPRI,'(10X,A,/)') 209 & ' - empty block' 210 ENDIF 211 ENDDO 212 ENDIF 213 ENDIF 214 215 IF ((ITER.EQ.1) .AND. (IRST.EQ.0)) THEN 216 WRITE(LUPRI,'(1X,A,I3,A,F23.16)') 217 & 'Iter.',ITER,': Coupled cluster MP2 energy : ',ECC2 218 ELSE 219 WRITE(LUPRI,'(1X,A,I3,A,F23.16)') 220 & 'Iter.',ITER,': Coupled cluster CC2 energy : ',ECC2 221 ENDIF 222 IF (IPRINT .GE. 2) THEN 223 WRITE(LUPRI,'(33X,A,9X,1P,D14.6)') 224 & 'DeltaE : ',EN2-EN1 225 END IF 226 227 CALL FLSHFO(LUPRI) 228 229C Save information for this iteration. 230C ------------------------------------ 231 232 CALL CC2_SAVTAM(WORK(KT1AM),WORK(KOMEG1),ECC2) 233 234 235C Check convergence. 236C ------------------ 237 238 DELE = DABS(EN2-EN1) 239 LCONVG = (DELE.LE.THRENR) .AND. (O1NM.LE.THRVEC) 240 IF (.NOT. LCONVG) THEN 241 CALL CCSD_NXTAM(WORK(KT1AM),DUM1,DUM12,WORK(KOMEG1), 242 & DUM2,DUM22,WORK(KFOCKD),.FALSE.,1,0.0D0) 243 CALL CCSD_DIIS(WORK(KOMEG1),WORK(KT1AM),NT1AMX,ITER) 244 EN1 = EN2 245 GOTO 100 246 ENDIF 247 248C CC2 wave function has converged. 249C -------------------------------- 250 251 WRITE(LUPRI,'(/,1X,A,D11.4,A,F16.9)') 252 & 'CC2 energy converged to within ',THRENR,' is ',ECC2 253 WRITE(LUPRI,'(1X,A,F16.9)') 254 & '- the CC2 correlation energy contribution is ',ECC2-ESCF 255 WRITE(LUPRI,'(1X,A,6X,1P,D14.4)') 256 & 'Final norm of the CC2 vector function is ',O1NM 257 TIMTOT = SECOND() - TIMTOT 258 WRITE(LUPRI,'(1X,A,7X,F14.2,A,/)') 259 & 'Total time used for CC2 optimization is ',TIMTOT,' seconds' 260 261 CALL FLSHFO(LUPRI) 262 263C Calculate and write Fock matrix 264C ------------------------------- 265 266 IF (RSPIM2) THEN 267 CALL CHO_FCKH(WORK(KT1AM),WORK(KEND1),LWRK1) 268 END IF 269 270C Close DIIS files. 271C ----------------- 272 273#if defined (SYS_CRAY) 274 CALL WCLOSE('CC_DIIS',IERR) 275 INFO = ISHELL('rm CC_DIIS') 276#endif 277#if !defined (SYS_CRAY) 278 IF (LUTDIS .GT. 0) THEN 279 INQUIRE(UNIT=LUTDIS,OPENED=LXDIS,ERR=121) 280 IF (LXDIS) CALL GPCLOSE(LUTDIS,'DELETE') 281 ENDIF 282 IF (LUSDIS .GT. 0) THEN 283 INQUIRE(UNIT=LUSDIS,OPENED=LXDIS,ERR=121) 284 CALL GPCLOSE(LUSDIS,'DELETE') 285 ENDIF 286 121 CONTINUE 287#endif 288 289C Print largest components of 0th order wave function. 290C ---------------------------------------------------- 291 292 IF (IPRINT .GT. 2) THEN 293 CALL AROUND('Largest amplitudes in converged solution') 294 CALL CC_PRAM(WORK(KT1AM),PT1,1) 295 ENDIF 296 297C Save a copy on rsp-file system. 298C ------------------------------- 299 300 KT0AM = KT1AM + NT1AMX 301 KEND1 = KT0AM + 2*NALLAI(1) 302 LWRK1 = LWORK - KEND1 + 1 303 304 IF (LWRK1 .LT. 0) THEN 305 WRITE(LUPRI,'(//,5X,A,A)') 306 & 'Insufficient memory for amplitude saving in ',SECNAM 307 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 308 & 'Need : ',KEND1-1, 309 & 'Available: ',LWORK 310 CALL QUIT('Insufficient memory in '//SECNAM) 311 ENDIF 312 313 CALL DZERO(WORK(KT0AM),2*NALLAI(1)) 314 315 IOPT = 7 316 CALL CC_WRRSP('R0',0,1,IOPT,MODEL,WORK(KT0AM),WORK(KT1AM), 317 & DUMMY,WORK(KEND1),LWRK1) 318 319C Calculate global E intermediates for response. 320C ---------------------------------------------- 321 322 IF (RSPIM2.AND.(.NOT.IMSKIP)) THEN 323 324 RSPIM = RSPIM2 325 326 WRITE(LUPRI,'(/)') 327 CALL AROUND('Calculating Intermediates for CCLR') 328 WRITE(LUPRI,'(/)') 329 330 KEND1 = KT1AM + NT1AMX 331 LWRK1 = LWORK - KEND1 + 1 332 CALL CC_CHOEIM(WORK(KFOCKD),WORK(KT1AM),WORK(KEND1),LWRK1, 333 & .TRUE.) 334 335 IF (IPRINT .GT. 1) WRITE(LUPRI,'(/)') 336 WRITE(LUPRI,'(12X,A)') 'E-intermediates calculated' 337 338 ELSE IF (RSPIM2.AND.IMSKIP) THEN 339 340 RSPIM = RSPIM2 341 WRITE(LUPRI,'(12X,A)') 342 & 'Intermediates assumed to be restart IM. ' 343 344 ENDIF 345 346C Delete files with Y intermediates. 347C (Might have been done already in cc_choeim). 348C -------------------------------------------- 349 350 DO ISYCHO = 1,NSYM 351 CALL CC_CYIOP(-1,ISYCHO,0) 352 CALL CC_CYIOP(0,ISYCHO,0) 353 ENDDO 354 355 CALL GETTIM(CEND,WEND) 356 CTOT = CEND - CSTR 357 WTOT = WEND - WSTR 358 CALL TIMTXT(' Total CPU time used in '//SECNAM//':',CTOT, 359 & LUPRI) 360 CALL TIMTXT(' Total wall time used in '//SECNAM//':',WTOT, 361 & LUPRI) 362 CALL QEXIT(SECNAM) 363 RETURN 364 END 365C /* Deck cc2_getsi */ 366 SUBROUTINE CC2_GETSI(FOCKD,WORK,LWORK,ESCF,POTNUC) 367C 368C Thomas Bondo Pedersen, July 2002. 369C 370C Purpose: Read Fock diagonal and SCF energy from SIRIUS 371C interface file. 372C 373#include "implicit.h" 374 DIMENSION FOCKD(*), WORK(LWORK) 375#include "dummy.h" 376 377 CALL CHO_RDSIR(POTNUC,ESCF,FOCKD,DUMMY,WORK,LWORK,.TRUE.,.TRUE., 378 & .FALSE.) 379 380 RETURN 381 END 382C /* Deck cc2_initam */ 383 SUBROUTINE CC2_INITAM(FOCKD,T1AM,OMEGA1,WORK,LWORK,ECC2,IRST) 384C 385C Thomas Bondo Pedersen, August 2002. 386C 387C Purpose: Set up initial Cholesky CC2 amplitudes. 388C Restart if requested and if possible. 389C 390C On exit: 391C 392C IRST = 0: no restart [if requested, it failed] 393C IRST = 1: "full" restart: T1AM, OMEGA1 and ECC2 read from disk. 394C IRST = 2: restart only with T1AM [as in conventional code]. 395C 396C The reason for trying to read not only T1AM but also OMEGA1 and ECC2 397C is that with conventional amplitude restart, the Cholesky CC2 energy 398C solver will spend 2 iterations before realizing that the amplitudes 399C are converged: the first iteration to get energy and vector function, 400C and convergence will be found after the second (if the initial ones 401C were in fact converged, obviously). 402C 403#include "implicit.h" 404 DIMENSION FOCKD(*), T1AM(*), OMEGA1(*), WORK(LWORK) 405#include "ccorb.h" 406#include "ccsdsym.h" 407#include "ccsdinp.h" 408#include "priunit.h" 409#include "dummy.h" 410 411 CHARACTER*10 SECNAM 412 PARAMETER (SECNAM = 'CC2_INITAM') 413 414 LOGICAL RSTART 415 416 CHARACTER*10 FILRST, MODEL 417 PARAMETER (FILRST = 'CHOCC2.RST') 418 419 PARAMETER (ZERO = 0.0D0) 420 421C Set model for reading in CC_RDRSP. 422C ---------------------------------- 423 424 IF (CC2) THEN 425 MODEL = 'CC2 ' 426 ELSE 427 WRITE(LUPRI,'(//,5X,A,A)') 428 & SECNAM,': Error: only CC2 model recognized!' 429 CALL QUIT('Error in '//SECNAM) 430 ENDIF 431 432C Default is MP2 guess. 433C --------------------- 434 435 IRST = 0 436 437C If requested, try to restart. 438C ----------------------------- 439 440 IF (CCRSTR) THEN 441 442C First try full restart. 443C ----------------------- 444 445 INQUIRE(FILE=FILRST,EXIST=RSTART,ERR=200) 446 IF (RSTART) THEN 447 IRST = 1 448 ELSE 449 IRST = 2 450 ENDIF 451 GO TO 100 452 453 200 CONTINUE 454 IF (IPRINT .GE. 1) THEN 455 WRITE(LUPRI,'(5X,A,A)') 456 & SECNAM,': an error ocurred inquering file ',FILRST 457 WRITE(LUPRI,'(5X,A,A)') 458 & SECNAM,': trying to restart from amplitudes only.' 459 ENDIF 460 IRST = 2 461 462 100 CONTINUE 463 464 IF (IRST .EQ. 1) THEN 465 466C FILRST exists, try to read data while checking that the 467C amount of data [NT1AMX] is consistent with current run. 468C If not, reset IRST to 2. 469C ------------------------------------------------------- 470 471 IFAIL = 0 472 CALL CHO_RDRST(ECC2,T1AM,OMEGA1,.TRUE.,.TRUE.,.TRUE.,IFAIL) 473 IF (IFAIL .EQ. 0) THEN 474 IF (IPRINT .GE. 1) THEN 475 WRITE(LUPRI,'(5X,A,A,A,A)') 476 & SECNAM,': file ',FILRST, 477 & ' detected and accepted for full restart.' 478 ENDIF 479 ELSE 480 IF (IPRINT .GE. 1) THEN 481 WRITE(LUPRI,'(5X,A,A,A,A)') 482 & SECNAM,': file ',FILRST, 483 & ' detected but rejected for restart.' 484 ENDIF 485 IRST = 2 486 ENDIF 487 488 ENDIF 489 490 ENDIF 491 492C Action taken according to restart (and which) or not. 493C ===================================================== 494 495 IF (IRST .EQ. 1) THEN 496 497C T1AM and OMEGA1 (and ECC2) were successfully read from FILRST. 498C Generate new amplitudes by perturbation theory. 499C -------------------------------------------------------------- 500 501 ISYALF = 1 502 CALL CCSD_NXTAM(T1AM,DUMMY,DUM12,OMEGA1,DUMMY,DUM22,FOCKD, 503 & .FALSE.,ISYALF,0.0D0) 504 CALL DCOPY(NT1AMX,OMEGA1,1,T1AM,1) 505 506 ELSE IF (IRST .EQ. 2) THEN 507 508C For whatever reason, restart from FILRST failed. Now try 509C to get amplitudes from 'R0' file or, if that fails, the 510C original ("old restart") CCSD_TAM file. 511C -------------------------------------------------------- 512 513 ECC2 = ZERO 514 515 ILST = 1 516 ISYM = 1 517 IOPT = 33 518 CALL CC_RDRSP('R0 ',ILST,ISYM,IOPT,MODEL,T1AM,DUMMY) 519 520 IF (IOPT .EQ. 33) THEN 521 INQUIRE(FILE='CCSD_TAM',EXIST=RSTART,ERR=1000) 522 GO TO 1100 523 1000 RSTART = .FALSE. 524 1100 CONTINUE 525 IF (RSTART) THEN 526 LUTAM = -1 527 CALL GPOPEN(LUTAM,'CCSD_TAM','UNKNOWN',' ','UNFORMATTED', 528 * IDUMMY,.FALSE.) 529 REWIND(LUTAM) 530 READ(LUTAM) (T1AM(I), I = 1,NT1AMX) 531 CALL GPCLOSE(LUTAM,'KEEP') 532 IF (IPRINT .GE. 1) THEN 533 WRITE(LUPRI,'(5X,A,A)') 534 & SECNAM,': "old" amplitude restart succeeded.' 535 ENDIF 536 ELSE 537 IF (IPRINT .GE. 1) THEN 538 WRITE(LUPRI,'(5X,A,A)') 539 & SECNAM,': amplitude restart failed.' 540 WRITE(LUPRI,'(5X,A,A)') 541 & SECNAM,': setting up MP2 initial guess instead.' 542 ENDIF 543 CALL DZERO(T1AM,NT1AMX) 544 IRST = 0 545 ENDIF 546 ELSE 547 IF (IPRINT .GE. 1) THEN 548 WRITE(LUPRI,'(5X,A,A)') 549 & SECNAM,': amplitude restart succeeded.' 550 ENDIF 551 ENDIF 552 553 ELSE IF (IRST .EQ. 0) THEN 554 555C MP2 initial guess. 556C ------------------ 557 558 IF (IPRINT .GT. 10) THEN 559 WRITE(LUPRI,'(5X,A,A)') 560 & SECNAM,': Setting up MP2 initial guess.' 561 ENDIF 562 563 ECC2 = ZERO 564 CALL DZERO(T1AM,NT1AMX) 565 566 ELSE 567 568C Logical error. 569C -------------- 570 571 WRITE(LUPRI,'(//,5X,A,A)') 572 & 'Logical error ocurred in ',SECNAM 573 WRITE(LUPRI,'(5X,A,I10,/)') 574 & 'IRST out of bounds: ',IRST 575 CALL QUIT('Logical error in '//SECNAM) 576 577 ENDIF 578 579C Print. 580C ------ 581 582 IF (CCRSTR .AND. (IRST.EQ.0)) THEN 583 WRITE(LUPRI,'(5X,A,A)') 584 & SECNAM,': Restart failed. Initial amplitudes are MP2 guess' 585 ENDIF 586 587 IF ((IPRINT.GT.100) .OR. DEBUG) THEN 588 T1NRM = DSQRT(DDOT(NT1AMX,T1AM,1,T1AM,1)) 589 CALL AROUND('Initial CC2 amplitudes in '//SECNAM) 590 WRITE(LUPRI,'(10X,A,I2)') 591 & 'IRST =',IRST 592 WRITE(LUPRI,'(10X,A,1P,D22.15,/)') 593 & 'Norm of amplitudes: ',T1NRM 594 CALL NOCC_PRT(T1AM,1,'AI ') 595 ENDIF 596 597 RETURN 598 END 599C /* Deck cc2_savtam */ 600 SUBROUTINE CC2_SAVTAM(T1AM,OMEGA1,ECC2) 601C 602C Thomas Bondo Pedersen, August 2002. 603C 604C Purpose: Save CC2 amplitudes, vector function, and energy 605C on disk. 606C 607#include "implicit.h" 608 DIMENSION T1AM(*), OMEGA1(*) 609 610 CALL CHO_WRRST(ECC2,T1AM,OMEGA1) 611 612 RETURN 613 END 614C /* Deck cc2_trf0 */ 615 SUBROUTINE CC2_TRF0(WORK,LWORK) 616C 617C Thomas Bondo Pedersen, July 2002. 618C 619C Purpose: T1-independent MO transformations of Cholesky 620C vectors; results are written to disk. 621C 622#include "implicit.h" 623 DIMENSION WORK(LWORK) 624#include "ccorb.h" 625#include "ccsdsym.h" 626#include "priunit.h" 627 628 CHARACTER*8 SECNAM 629 PARAMETER (SECNAM = 'CC2_TRF0') 630 631C Allocation. 632C ----------- 633 634 KCMO = 1 635 KEND1 = KCMO + NLAMDS 636 LWRK1 = LWORK - KEND1 + 1 637 638 IF (LWRK1 .LT. 0) THEN 639 WRITE(LUPRI,'(//,5X,A,A,A)') 640 & 'Insufficient memory in ',SECNAM,' - allocation: 1' 641 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 642 & 'Need (more than): ',KEND1-1, 643 & 'Available : ',LWORK 644 CALL QUIT('Insufficient memory in '//SECNAM) 645 ENDIF 646 647C Read the MO coefficient matrix. Change from SIRIUS to CC ordering. 648C ------------------------------------------------------------------ 649 650 CALL CHO_RDSIR(DUM1,DUM2,DUM3,WORK(KCMO),WORK(KEND1),LWRK1, 651 & .FALSE.,.FALSE.,.TRUE.) 652 653C Transform Cholesky vectors: L(ia,J). 654C Result vectors are stored on disk. 655C ------------------------------------ 656 657 CALL CHO_TRFAI(WORK(KCMO),WORK(KEND1),LWRK1) 658 659C Transform one-electron integrals: h(ia). 660C Result is stored on disk. 661C ---------------------------------------- 662Ctbp: obsolete in new version 663Ctbp CALL CC2_ONETRF(WORK(KCMO),WORK(KEND1),LWRK1) 664Ctbp 665 RETURN 666 END 667C /* Deck cc2_onetrf */ 668 SUBROUTINE CC2_ONETRF(CMO,WORK,LWORK) 669C 670C Thomas Bondo Pedersen, July 2002. 671C 672C Purpose: Partially transform AO one-electron integrals 673C to MO basis. Result array h(ia) stored on disk. 674C 675#include "implicit.h" 676 DIMENSION CMO(*), WORK(LWORK) 677#include "ccorb.h" 678#include "ccsdsym.h" 679#include "priunit.h" 680 681 CHARACTER*10 SECNAM 682 PARAMETER (SECNAM = 'CC2_ONETRF') 683 684 PARAMETER (ITYP = 1, IOPEN = -1, IKEEP = 1) 685 PARAMETER (ZERO = 0.00D0, ONE = 1.00D0) 686 687C Allocation. 688C ----------- 689 690 MAXALI = 0 691 DO ISYM = 1,NSYM 692 MAXALI = MAX(MAXALI,NBAS(ISYM)*NRHF(ISYM)) 693 ENDDO 694 695 KONEL = 1 696 KHIA = KONEL + N2BST(1) 697 KSCR1 = KHIA + NT1AM(1) 698 KEND1 = KSCR1 + MAXALI 699 LWRK1 = LWORK - KEND1 + 1 700 701 IF (LWRK1 .LT. 0) THEN 702 WRITE(LUPRI,'(//,5X,A,A,A)') 703 & 'Insufficient memory in ',SECNAM,' - allocation: ONEL' 704 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 705 & 'Need : ',KEND1-1, 706 & 'Available: ',LWORK 707 CALL QUIT('Insufficient memory in '//SECNAM) 708 ENDIF 709 710C Read AO integrals. 711C ------------------ 712 713 CALL CCRHS_ONEAO(WORK(KONEL),WORK(KEND1),LWRK1) 714 715C Transform. 716C ---------- 717 718 DO ISYMI = 1,NSYM 719 720 ISYMG = ISYMI 721 ISYMD = ISYMG 722 ISYMA = ISYMD 723 724 NTOTG = MAX(NBAS(ISYMG),1) 725 NTOTD = MAX(NBAS(ISYMD),1) 726 NTOTA = MAX(NVIR(ISYMA),1) 727 728 KOFF1 = KONEL + IAODIS(ISYMG,ISYMD) 729 KOFF2 = IGLMRH(ISYMG,ISYMI) + 1 730 KOFF3 = IGLMVI(ISYMD,ISYMA) + 1 731 KOFF4 = KHIA + IT1AM(ISYMA,ISYMI) 732 733 CALL DGEMM('T','N',NBAS(ISYMD),NRHF(ISYMI),NBAS(ISYMG), 734 & ONE,WORK(KOFF1),NTOTG,CMO(KOFF2),NTOTG, 735 & ZERO,WORK(KSCR1),NTOTD) 736 737 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMD), 738 & ONE,CMO(KOFF3),NTOTD,WORK(KSCR1),NTOTD, 739 & ZERO,WORK(KOFF4),NTOTA) 740 741 ENDDO 742 743C Write h(ia) to disk. 744C -------------------- 745 746 CALL ONEL_OP(IOPEN,ITYP,LUHIA) 747 CALL CHO_MOWRITE(WORK(KHIA),NT1AM(1),1,1,LUHIA) 748 CALL ONEL_OP(IKEEP,ITYP,LUHIA) 749 750 RETURN 751 END 752C /* Deck cc2_chorhs */ 753 SUBROUTINE CC2_CHORHS(FOCKD,T1AM,OMEGA1,WORK,LWORK,T2NRM,ECC2, 754 & ESCF) 755C 756C Thomas Bondo Pedersen, August 2002. 757C 758C Purpose: Calculate the CC2 vector function from Cholesky 759C decomposed two-electron integrals. The CC2 energy 760C ECC2 is calculated as a byproduct; note that it 761C contains the SCF energy. 762C 763C Notes: 764C (1) The T1-independent MO transformations are assumed 765C available on disk, i.e. L(ia,J). 766C (2) The T2 amplitudes may be treated in 3 different ways: 767C - by straightforward calculation from (ai|bj) integrals 768C represented by the transformed AO Cholesky vectors. 769C - by calculation from separately Cholesky decomposed 770C (ai|bj) integrals. 771C - by decomposition of the amplitudes themselves. 772C 773#include "implicit.h" 774 DIMENSION FOCKD(*), T1AM(*), OMEGA1(*), WORK(LWORK) 775#include "maxorb.h" 776#include "ccdeco.h" 777#include "ccorb.h" 778#include "ccsdsym.h" 779#include "ccsdinp.h" 780#include "priunit.h" 781 782 CHARACTER*10 SECNAM 783 PARAMETER (SECNAM = 'CC2_CHORHS') 784 785 LOGICAL LOCDBG 786 PARAMETER (LOCDBG = .FALSE.) 787 788 PARAMETER (INFTOT = 2, INFO = 3) 789 PARAMETER (ZERO = 0.0D0) 790 791C Start timing this routine. 792C -------------------------- 793 794 TIMTOT = SECOND() 795 796 IF (LOCDBG) THEN 797 OESUM = ZERO 798 DO ISYMI = 1,NSYM 799 DO I = 1,NRHF(ISYMI) 800 KOFFI = IRHF(ISYMI) + I 801 OESUM = OESUM + FOCKD(KOFFI) 802 ENDDO 803 ENDDO 804 WRITE(LUPRI,'(5X,A,A,F16.9)') 805 & SECNAM,': debug: Sum of occupied orbital energies: ',OESUM 806 ENDIF 807 808C Initialization. 809C --------------- 810 811 ECC2 = ESCF 812 T2NRM = ZERO 813 814C Allocation. 815C ----------- 816 817 KTRACE = 1 818 KLAMDP = KTRACE + NUMCHO(1) 819 KLAMDH = KLAMDP + NGLMDT(1) 820 KEND1 = KLAMDH + NGLMDT(1) 821 LWRK1 = LWORK - KEND1 + 1 822 823 IF (LWRK1 .LT. 0) THEN 824 WRITE(LUPRI,'(//,5X,A,A,A)') 825 & 'Insufficient memory in ',SECNAM,' - allocation: Lambda' 826 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 827 & 'Need (more than): ',KEND1-1, 828 & 'Available : ',LWORK 829 CALL QUIT('Insufficient memory in '//SECNAM) 830 ENDIF 831 832C Calculate Lambda matrices. 833C -------------------------- 834 835 TIMLAM = SECOND() 836 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),T1AM,WORK(KEND1),LWRK1) 837 TIMLAM = SECOND() - TIMLAM 838 839C Calculate T1-dependent MO transformations. 840C Calculate T1-dependent energy contribution. 841C ------------------------------------------- 842 843 TIMTRL = SECOND() 844 ENERGY = ZERO 845 CALL CC2_TRFL(T1AM,WORK(KLAMDP),WORK(KLAMDH),WORK(KTRACE), 846 & WORK(KEND1),LWRK1,ENERGY) 847 ECC2 = ECC2 + ENERGY 848 TIMTRL = SECOND() - TIMTRL 849 850 IF (LOCDBG) THEN 851 WRITE(LUPRI,'(5X,A,A,F16.9,/,5X,A,A,F16.9)') 852 & SECNAM,': debug: T1 contribution to energy: ',ENERGY, 853 & SECNAM,': debug: Accumulated CC2 energy : ',ECC2 854 ENDIF 855 856C Initialize result vector. 857C ------------------------- 858 859 CALL CC2_INIOME(FOCKD,T1AM,OMEGA1,1) 860 861C Calculate Y-intermediates and the I-term. 862C The Y-intermediates are stored on disk, 863C and the I-term is added into OMEGA1. 864C Calculate T2-dependent energy contribution. 865C ------------------------------------------- 866 867 TIMYI = SECOND() 868 ENERGY = ZERO 869 CALL CHOCC2_YI(FOCKD,OMEGA1,WORK(KEND1),LWRK1,T2NRM,ENERGY) 870 ECC2 = ECC2 + ENERGY 871 TIMYI = SECOND() - TIMYI 872 873 IF (LOCDBG) THEN 874 OMNRM = DSQRT(DDOT(NT1AMX,OMEGA1,1,OMEGA1,1)) 875 WRITE(LUPRI,'(5X,A,A,F16.9,/,5X,A,A,F16.9)') 876 & SECNAM,': debug: T2 contribution to energy: ',ENERGY, 877 & SECNAM,': debug: Accumulated CC2 energy : ',ECC2 878 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 879 & SECNAM,': debug: Norm of Omega1 after CHOCC2_YI: ',OMNRM 880 ENDIF 881 882C Calculate the JG-term. 883C N.B.: TRACE is scaled by 2 on exit from CHOCC2_JG. 884C -------------------------------------------------- 885 886 TIMJG = SECOND() 887 CALL CHOCC2_JG(T1AM,OMEGA1,WORK(KLAMDP),WORK(KLAMDH), 888 & WORK(KTRACE),WORK(KEND1),LWRK1) 889 TIMJG = SECOND() - TIMJG 890 891 IF (LOCDBG) THEN 892 OMNRM = DSQRT(DDOT(NT1AMX,OMEGA1,1,OMEGA1,1)) 893 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 894 & SECNAM,': debug: Norm of Omega1 after CHOCC2_JG: ',OMNRM 895 ENDIF 896 897C Calculate the H term. 898C Pass full memory: trace and Lambda matrices no 899C longer needed. 900C ---------------------------------------------- 901 902 TIMJH = SECOND() 903 CALL CHOCC2_H(OMEGA1,WORK,LWORK) 904 TIMJH = SECOND() - TIMJH 905 906 IF (LOCDBG) THEN 907 OMNRM = DSQRT(DDOT(NT1AMX,OMEGA1,1,OMEGA1,1)) 908 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 909 & SECNAM,': debug: Norm of Omega1 after CHOCC2_H: ',OMNRM 910 ENDIF 911 912C Print section. 913C -------------- 914 915 IF (IPRINT .GT. INFTOT) THEN 916 TIMTOT = SECOND() - TIMTOT 917 WRITE(LUPRI,'(7X,A,F17.2,A)') 918 & 'Time used in RHS - TOTAL : ',TIMTOT,' seconds' 919 IF (IPRINT .GE. INFO) THEN 920 WRITE(LUPRI,'(7X,A,F17.2,A)') 921 & 'Time used in LAMMAT : ',TIMLAM,' seconds' 922 WRITE(LUPRI,'(7X,A,F17.2,A)') 923 & 'Time used in CC2_TRFL : ',TIMTRL,' seconds' 924 WRITE(LUPRI,'(7X,A,F17.2,A)') 925 & 'Time used in CHOCC2_YI : ',TIMYI,' seconds' 926 WRITE(LUPRI,'(7X,A,F17.2,A)') 927 & 'Time used in CHOCC2_JG : ',TIMJG,' seconds' 928 WRITE(LUPRI,'(7X,A,F17.2,A)') 929 & 'Time used in CHOCC2_H : ',TIMJH,' seconds' 930 ENDIF 931 ENDIF 932 933 RETURN 934 END 935C /* Deck cc2_trfl */ 936 SUBROUTINE CC2_TRFL(T1AM,XLAMDP,XLAMDH,TRACE,WORK,LWORK,ECC2) 937C 938C Thomas Bondo Pedersen, August 2002. 939C 940C Purpose: Calculate T1-dependent MO transformations for 941C Cholesky CC2. 942C The results are written to disk: 943C L(ij,J), L(ai,J), F(ia). 944C Calculate T1-dependent energy contribution. 945C On exit, the TRACE array contains the vector 946C TRACE(J) = sum(ai) L(ia,J) * T1AM(ai). 947C 948#include "implicit.h" 949 DIMENSION T1AM(*), XLAMDP(*), XLAMDH(*), TRACE(*), WORK(LWORK) 950#include "maxorb.h" 951#include "ccdeco.h" 952#include "ccorb.h" 953#include "ccsdsym.h" 954#include "priunit.h" 955 956 CHARACTER*8 SECNAM 957 PARAMETER (SECNAM = 'CC2_TRFL') 958 959 LOGICAL LOCDBG 960 PARAMETER (LOCDBG = .FALSE.) 961 962 PARAMETER (IOPEN = -1, IKEEP = 1) 963 PARAMETER (ZERO = 0.00D0, ONE = 1.00D0, TWO = 2.00D0) 964 965C Calculate F(ia), L(ij,J), and L(ai,J). 966C Calculate T1-dependent energy contributions. 967C -------------------------------------------- 968 969 KFOCK = 1 970 KEND1 = KFOCK + NT1AMX 971 LWRK1 = LWORK - KEND1 + 1 972 973 IF (LWRK1 .LT. 0) THEN 974 WRITE(LUPRI,'(//,5X,A,A,A)') 975 & 'Insufficient memory in ',SECNAM,' - allocation: Fock' 976 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 977 & 'Need : ',KEND1-1, 978 & 'Available: ',LWORK 979 CALL QUIT('Insufficient memory in '//SECNAM) 980 ENDIF 981 982 CALL DZERO(WORK(KFOCK),NT1AMX) 983 984 CALL CC2_TRF1(T1AM,XLAMDP,XLAMDH,TRACE,WORK(KFOCK), 985 & WORK(KEND1),LWRK1,ECC2) 986 987 IF (LOCDBG) THEN 988 FNORM = DSQRT(DDOT(NT1AMX,WORK(KFOCK),1,WORK(KFOCK),1)) 989 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 990 & SECNAM,': debug: Norm of F(ia) after CC2_TRF1: ',FNORM 991 ENDIF 992 993 CALL ONEL_OP(IOPEN,3,LUFOCK) 994 CALL CHO_MOWRITE(WORK(KFOCK),NT1AMX,1,1,LUFOCK) 995 CALL ONEL_OP(IKEEP,3,LUFOCK) 996 997 RETURN 998 END 999C /* Deck cc2_trf1 */ 1000 SUBROUTINE CC2_TRF1(T1AM,XLAMDP,XLAMDH,TRACE,FOCK,WORK,LWORK,ECC2) 1001C 1002C Thomas Bondo Pedersen, August 2002. 1003C 1004C Purpose: Calculate T1-dependent MO transformations 1005C F(ia), L(ji,J), and L(ai,J). 1006C The Fock matrix must be initialized. 1007C Calculate energy contribution from singles: 1008C 1009C ECC2 = ECC2 + 2 * sum(J) [sum(ai) T1AM(ai) * L(ia,J)]**2 1010C 1011C - sum(jiJ) [sum(a) T1AM(aj) * L(ia,J)] 1012C 1013C * [sum(b) T1AM(bi) * L(jb,J)] 1014C 1015C On exit, the TRACE array contains the vector 1016C 1017C TRACE(J) = sum(ai) L(ia,J) * T1AM(ai). 1018C 1019#include "implicit.h" 1020 DIMENSION T1AM(*), XLAMDP(*), XLAMDH(*), TRACE(*), FOCK(*) 1021 DIMENSION WORK(LWORK) 1022#include "maxorb.h" 1023#include "ccdeco.h" 1024#include "ccorb.h" 1025#include "ccsdsym.h" 1026#include "chocc2.h" 1027#include "priunit.h" 1028 1029 CHARACTER*8 SECNAM 1030 PARAMETER (SECNAM = 'CC2_TRF1') 1031 1032 LOGICAL LOCDBG 1033 PARAMETER (LOCDBG = .FALSE.) 1034 1035 PARAMETER (IOPEN = -1, IKEEP = 1, IOPTR = 2) 1036 PARAMETER (ITYIA = 1, ITYAI = 3, ITYJI = 4) 1037 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 1038 1039C Initialize energy contributions. 1040C -------------------------------- 1041 1042 ENRC = ZERO 1043 ENRX = ZERO 1044 1045C Read reduce index array. 1046C ------------------------ 1047 1048 KIND1 = 1 1049 CALL CC_GETIND1(WORK(KIND1),LWORK,LIND1) 1050 KEND0 = KIND1 + LIND1 1051 LWRK0 = LWORK - KEND0 + 1 1052 1053 IF (LWRK0 .LT. 0) THEN 1054 WRITE(LUPRI,'(//,5X,A,A,A)') 1055 & 'Insufficient memory in ',SECNAM,' - allocation: index' 1056 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1057 & 'Need (more than): ',KEND0-1, 1058 & 'Available : ',LWORK 1059 CALL QUIT('Insufficient memory in '//SECNAM) 1060 ENDIF 1061 1062C Start Cholesky symmetry loop. 1063C ----------------------------- 1064 1065 DO ISYCHO = 1,NSYM 1066 1067 IF (NUMCHO(ISYCHO) .LE. 0) GO TO 999 1068 IF (N2BST(ISYCHO) .LE. 0) GO TO 999 1069 1070C Open MO Cholesky files for this symmetry. 1071C ----------------------------------------- 1072 1073 CALL CHO_MOP(IOPEN,ITYJI,ISYCHO,LUJI,1,1) 1074 CALL CHO_MOP(IOPEN,ITYAI,ISYCHO,LUAI,1,1) 1075 CALL CHO_MOP(IOPEN,ITYIA,ISYCHO,LUIA,1,1) 1076 1077C Calculate size of scratch space needed for read 1078C and first half-transformation [L(gamma i)]. 1079C ----------------------------------------------- 1080 1081 MAXGI = 0 1082 DO ISYMI = 1,NSYM 1083 ISYMG = MULD2H(ISYMI,ISYCHO) 1084 MAXGI = MAX(MAXGI,NBAS(ISYMG)*NRHF(ISYMI)) 1085 ENDDO 1086 LREAD = MEMRD(1,ISYCHO,IOPTR) 1087 LSCR1 = MAX(LREAD,MAXGI) 1088 1089C Allocation. 1090C ----------- 1091 1092 KCHOL = KEND0 1093 KSCR1 = KCHOL + N2BST(ISYCHO) 1094 KEND1 = KSCR1 + LSCR1 1095 LWRK1 = LWORK - KEND1 + 1 1096 1097 IF (LWRK1 .LE. 0) THEN 1098 WRITE(LUPRI,'(//,5X,A,A,A)') 1099 & 'Insufficient memory in ',SECNAM,' - allocation: Cholesky' 1100 MREQ = KEND1 + NT1AM(ISYCHO) + NMATIJ(ISYCHO) - 1 1101 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1102 & 'Need (minimum): ',MREQ, 1103 & 'Available : ',LWORK 1104 CALL QUIT('Insufficient memory in '//SECNAM) 1105 ENDIF 1106 1107C Set up batch. 1108C ------------- 1109 1110 MINMEM = NT1AM(ISYCHO) + NMATIJ(ISYCHO) 1111 NVEC = MIN(LWRK1/MINMEM,NUMCHO(ISYCHO)) 1112 1113 IF (NVEC .LE. 0) THEN 1114 WRITE(LUPRI,'(//,5X,A,A)') 1115 & 'Insufficient memory for batch in ',SECNAM 1116 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10,/,5X,A,I10,/)') 1117 & 'Minimum memory required : ',MINMEM, 1118 & 'Memory available for batch: ',LWRK1, 1119 & 'Memory available in total : ',LWORK, 1120 & 'Number of Cholesky vectors: ',NUMCHO(ISYCHO) 1121 CALL QUIT('Insufficient memory for batch in '//SECNAM) 1122 ENDIF 1123 1124 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 1125 1126C Batch loop. 1127C ----------- 1128 1129 DO IBATCH = 1,NBATCH 1130 1131 NUMV = NVEC 1132 IF (IBATCH .EQ. NBATCH) THEN 1133 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 1134 ENDIF 1135 1136 JVEC1 = NVEC*(IBATCH - 1) + 1 1137 1138C Complete allocation. 1139C -------------------- 1140 1141 KCHJI = KEND1 1142 KCHAI = KCHJI + NMATIJ(ISYCHO)*NUMV 1143 1144 IF (LOCDBG) THEN 1145 KEND2 = KCHAI + NT1AM(ISYCHO)*NUMV 1146 LWRK2 = LWORK - KEND2 + 1 1147 IF (LWRK2 .LT. 0) THEN 1148 WRITE(LUPRI,'(//,5X,A,A,A)') 1149 & 'Error in batching in ',SECNAM,':' 1150 WRITE(LUPRI,*) 1151 & 'ISYCHO, NUMCHO: ',ISYCHO,NUMCHO(ISYCHO) 1152 WRITE(LUPRI,*) 1153 & 'NVEC, NUMV : ',NVEC,NUMV 1154 WRITE(LUPRI,*) 1155 & 'IBATCH, NBATCH: ',IBATCH,NBATCH 1156 WRITE(LUPRI,*) 1157 & 'LWRK1, MINMEM : ',LWRK1,MINMEM 1158 WRITE(LUPRI,*) 1159 & 'LWORK, KEND1-1: ',LWORK,KEND1-1 1160 CALL QUIT('Error in '//SECNAM) 1161 ENDIF 1162 ENDIF 1163 1164C Loop over vectors in this batch. 1165C -------------------------------- 1166 1167 DO LVEC = 1,NUMV 1168 1169C Read current vector. 1170C -------------------- 1171 1172 JVEC = JVEC1 + LVEC - 1 1173 CALL CHO_READN(WORK(KCHOL),JVEC,1,WORK(KIND1),IDUM2, 1174 & ISYCHO,IOPTR,WORK(KSCR1),LSCR1) 1175 1176C Transform. 1177C ---------- 1178 1179 DO ISYMI = 1,NSYM 1180 1181 NI = NRHF(ISYMI) 1182 IF (NI .LE. 0) GOTO 998 1183 1184 ISYMA = MULD2H(ISYMI,ISYCHO) 1185 ISYMJ = ISYMA 1186 ISYMG = ISYMA 1187 ISYMD = ISYMI 1188 1189 NG = NBAS(ISYMG) 1190 NTOTG = MAX(NG,1) 1191 ND = NBAS(ISYMD) 1192 NTOTD = MAX(ND,1) 1193 NA = NVIR(ISYMA) 1194 NTOTA = MAX(NA,1) 1195 NJ = NRHF(ISYMJ) 1196 NTOTJ = MAX(NJ,1) 1197 1198 KOFF1 = KCHOL + IAODIS(ISYMG,ISYMD) 1199 KOFF2 = IGLMRH(ISYMD,ISYMI) + 1 1200 1201 CALL DGEMM('N','N',NG,NI,ND, 1202 & ONE,WORK(KOFF1),NTOTG,XLAMDH(KOFF2),NTOTD, 1203 & ZERO,WORK(KSCR1),NTOTG) 1204 1205 KOFF1 = IGLMRH(ISYMG,ISYMJ) + 1 1206 KOFF2 = KCHJI + NMATIJ(ISYCHO)*(LVEC - 1) 1207 & + IMATIJ(ISYMJ,ISYMI) 1208 1209 CALL DGEMM('T','N',NJ,NI,NG, 1210 & ONE,XLAMDP(KOFF1),NTOTG,WORK(KSCR1),NTOTG, 1211 & ZERO,WORK(KOFF2),NTOTJ) 1212 1213 KOFF1 = IGLMVI(ISYMG,ISYMA) + 1 1214 KOFF2 = KCHAI + NT1AM(ISYCHO)*(LVEC - 1) 1215 & + IT1AM(ISYMA,ISYMI) 1216 1217 CALL DGEMM('T','N',NA,NI,NG, 1218 & ONE,XLAMDP(KOFF1),NTOTG,WORK(KSCR1),NTOTG, 1219 & ZERO,WORK(KOFF2),NTOTA) 1220 1221 998 CONTINUE 1222 1223 ENDDO 1224 1225 ENDDO 1226 1227C Write this batch to disk. 1228C ------------------------- 1229 1230 CALL CHO_MOWRITE(WORK(KCHJI),NMATIJ(ISYCHO),NUMV,JVEC1,LUJI) 1231 CALL CHO_MOWRITE(WORK(KCHAI),NT1AM(ISYCHO),NUMV,JVEC1,LUAI) 1232 1233C Read L(ia,J) into the space of L(ai,J). 1234C --------------------------------------- 1235 1236 CALL CHO_MOREAD(WORK(KCHAI),NT1AM(ISYCHO),NUMV,JVEC1,LUIA) 1237 1238 IF (ISYCHO .EQ. 1) THEN 1239 1240C Calculate TRACE vector. 1241C ----------------------- 1242 1243 NT1T = MAX(NT1AMX,1) 1244 1245 CALL DGEMV('T',NT1AMX,NUMV, 1246 & ONE,WORK(KCHAI),NT1T,T1AM,1, 1247 & ZERO,TRACE(JVEC1),1) 1248 1249C Calculate Coulomb contribution to F(ia). 1250C ---------------------------------------- 1251 1252 CALL DGEMV('N',NT1AMX,NUMV, 1253 & TWO,WORK(KCHAI),NT1T,TRACE(JVEC1),1, 1254 & ONE,FOCK,1) 1255 1256 ENDIF 1257 1258C Calculate exchange contribution to F(ia) and energy. 1259C First L(ji,J) vector is overwritten here!!! 1260C ---------------------------------------------------- 1261 1262 KUMAT = KCHJI 1263 1264 DO LVEC = 1,NUMV 1265 1266C Calculate U(ji) = sum(a) T1AM(aj) * L(ia,J). 1267C -------------------------------------------- 1268 1269 DO ISYMI = 1,NSYM 1270 1271 ISYMJ = MULD2H(ISYMI,ISYCHO) 1272 ISYMA = ISYMJ 1273 1274 KOFF1 = IT1AM(ISYMA,ISYMJ) + 1 1275 KOFF2 = KCHAI + NT1AM(ISYCHO)*(LVEC - 1) 1276 & + IT1AM(ISYMA,ISYMI) 1277 KOFF3 = KUMAT + IMATIJ(ISYMJ,ISYMI) 1278 1279 NA = NVIR(ISYMA) 1280 NTOTA = MAX(NA,1) 1281 NJ = NRHF(ISYMJ) 1282 NTOTJ = MAX(NJ,1) 1283 NI = NRHF(ISYMI) 1284 1285 CALL DGEMM('T','N',NJ,NI,NA, 1286 & ONE,T1AM(KOFF1),NTOTA,WORK(KOFF2),NTOTA, 1287 & ZERO,WORK(KOFF3),NTOTJ) 1288 1289 ENDDO 1290 1291C Calculate exchange contribution to F(ia). 1292C ----------------------------------------- 1293 1294 DO ISYMK = 1,NSYM 1295 1296 ISYMI = MULD2H(ISYMK,ISYCHO) 1297 ISYMA = ISYMI 1298 1299 KOFF1 = KCHAI + NT1AM(ISYCHO)*(LVEC - 1) 1300 & + IT1AM(ISYMA,ISYMK) 1301 KOFF2 = KUMAT + IMATIJ(ISYMK,ISYMI) 1302 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 1303 1304 NA = NVIR(ISYMA) 1305 NTOTA = MAX(NA,1) 1306 NI = NRHF(ISYMI) 1307 NK = NRHF(ISYMK) 1308 NTOTK = MAX(NK,1) 1309 1310 CALL DGEMM('N','N',NA,NI,NK, 1311 & XMONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOTK, 1312 & ONE,FOCK(KOFF3),NTOTA) 1313 1314 ENDDO 1315 1316C Calculate "exchange" energy contribution. 1317C ----------------------------------------- 1318 1319 DO ISYMI = 1,NSYM 1320 1321 ISYMJ = MULD2H(ISYMI,ISYCHO) 1322 1323 IF (NRHF(ISYMJ) .GT. 0) THEN 1324 1325 DO I = 1,NRHF(ISYMI) 1326 1327 KOFF1 = KUMAT + IMATIJ(ISYMJ,ISYMI) 1328 & + NRHF(ISYMJ)*(I - 1) 1329 KOFF2 = KUMAT + IMATIJ(ISYMI,ISYMJ) 1330 & + I - 1 1331 1332 ENRX = ENRX 1333 & + DDOT(NRHF(ISYMJ),WORK(KOFF1),1, 1334 & WORK(KOFF2),NRHF(ISYMI)) 1335 1336 ENDDO 1337 1338 ENDIF 1339 1340 ENDDO 1341 1342 ENDDO 1343 1344 ENDDO 1345 1346C Close MO Cholesky files for this symmetry. 1347C ------------------------------------------ 1348 1349 CALL CHO_MOP(IKEEP,ITYIA,ISYCHO,LUIA,1,1) 1350 CALL CHO_MOP(IKEEP,ITYAI,ISYCHO,LUAI,1,1) 1351 CALL CHO_MOP(IKEEP,ITYJI,ISYCHO,LUJI,1,1) 1352 1353 999 CONTINUE 1354 1355 ENDDO 1356 1357C Calculate Coulomb energy contribution and sum. 1358C ---------------------------------------------- 1359 1360 ENRC = TWO*DDOT(NUMCHO(1),TRACE,1,TRACE,1) 1361 ECC2 = ECC2 + ENRC - ENRX 1362 1363 IF (LOCDBG) THEN 1364 FNORM = DSQRT(DDOT(NT1AMX,FOCK,1,FOCK,1)) 1365 TNORM = DSQRT(DDOT(NUMCHO(1),TRACE,1,TRACE,1)) 1366 WRITE(LUPRI,'(5X,A,A,F16.9,/,5X,A,A,F16.9,/,5X,A,A,F16.9)') 1367 & SECNAM,': debug: Coulomb contribution to energy: ',ENRC, 1368 & SECNAM,': debug: Exch. contribution to energy: ',-ENRX, 1369 & SECNAM,': debug: Accumulated CC2 energy : ',ECC2 1370 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 1371 & SECNAM,': debug: Norm of F(ia): ',FNORM 1372 WRITE(LUPRI,'(5X,A,A,1P,D15.6)') 1373 & SECNAM,': debug: Norm of TRACE: ',TNORM 1374 ENDIF 1375 1376 RETURN 1377 END 1378C /* Deck cc2_choamd */ 1379 SUBROUTINE CC2_CHOAMD(DIAG,FOCKD,WORK,LWORK,ISYM,LUCHMO) 1380C 1381C Thomas Bondo Pedersen, August 2002. 1382C 1383C Purpose: Calculate (minus) the diagonal of 1384C the CC2 amplitudes of symmetry ISYM. The calculation 1385C uses the (ai|bj) integral representation 1386C from file LUCHMO; this file is assumed 1387C to be open and to contain the transformed 1388C AO Cholesky vectors L(ai,J). 1389C 1390#include "implicit.h" 1391 DIMENSION DIAG(*), FOCKD(*), WORK(LWORK) 1392#include "priunit.h" 1393#include "ccorb.h" 1394#include "ccsdsym.h" 1395#include "ccsdinp.h" 1396 1397 INTEGER AI 1398 1399 CHARACTER*10 SECNAM 1400 PARAMETER (SECNAM = 'CC2_CHOAMD') 1401 1402 PARAMETER (INFO = 20) 1403 PARAMETER (TWO = 2.0D0) 1404 1405C Start timing. 1406C ------------- 1407 1408 TIMTOT = SECOND() 1409 1410C Calculate the (ai|ai) integral diagonal. 1411C ---------------------------------------- 1412 1413 TIMDIA = SECOND() 1414 CALL CHO_AIBJD(DIAG,WORK,LWORK,ISYM,LUCHMO) 1415 TIMDIA = SECOND() - TIMDIA 1416 1417C Divide by the orbital energy denominators. 1418C t(ai,ai) = (ai|ai)/[2*e(a)-2*e(i)]. 1419C ------------------------------------------ 1420 1421 TIMDIV = SECOND() 1422 DO ISYMI = 1,NSYM 1423 1424 ISYMA = MULD2H(ISYMI,ISYM) 1425 1426 DO I = 1,NRHF(ISYMI) 1427 1428 KOFFI = IRHF(ISYMI) + I 1429 1430 DO A = 1,NVIR(ISYMA) 1431 1432 KOFFA = IVIR(ISYMA) + A 1433 1434 AI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 1435 1436 DENOM = TWO*(FOCKD(KOFFA)-FOCKD(KOFFI)) 1437 DIAG(AI) = DIAG(AI)/DENOM 1438 1439 ENDDO 1440 1441 ENDDO 1442 1443 ENDDO 1444 TIMDIV = SECOND() - TIMDIV 1445 1446C Print. 1447C ------ 1448 1449 IF (IPRINT .GE. INFO) THEN 1450 TIMTOT = SECOND() - TIMTOT 1451 XNRM = DSQRT(DDOT(NT1AM(ISYM),DIAG,1,DIAG,1)) 1452 CALL HEADER('Information from '//SECNAM,-1) 1453 WRITE(LUPRI,'(5X,A,I1)') 1454 & 'CC2 amplitude diagonal calculated, symmetry block:',ISYM 1455 WRITE(LUPRI,'(5X,A,1P,D30.15,/)') 1456 & 'Norm of calculated diagonal: ',XNRM 1457 WRITE(LUPRI,'(5X,A,F10.2,A)') 1458 & 'Time for calculating (ai|ai) diagonal: ',TIMDIA,' seconds' 1459 WRITE(LUPRI,'(5X,A,F10.2,A)') 1460 & 'Time for orbital energy denominators : ',TIMDIV,' seconds' 1461 WRITE(LUPRI,'(5X,A)') 1462 & '---------------------------------------------------------' 1463 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 1464 & 'Total time : ',TIMTOT,' seconds' 1465 ENDIF 1466 1467 RETURN 1468 END 1469C /* Deck cc2_choamd_i */ 1470 SUBROUTINE CC2_CHOAMD_I(DIAG,FOCKD,CHAI,ISYM) 1471C 1472C Thomas Bondo Pedersen, August 2002. 1473C 1474C Purpose: Calculate (minus) the diagonal of 1475C the CC2 amplitudes of symmetry ISYM. 1476C CHAI must contain the L(ai,J) vectors. 1477C 1478#include "implicit.h" 1479 DIMENSION DIAG(*), FOCKD(*), CHAI(*) 1480#include "priunit.h" 1481#include "ccorb.h" 1482#include "ccsdsym.h" 1483#include "ccsdinp.h" 1484 1485 INTEGER AI 1486 1487 CHARACTER*12 SECNAM 1488 PARAMETER (SECNAM = 'CC2_CHOAMD_I') 1489 1490 PARAMETER (INFO = 20) 1491 PARAMETER (TWO = 2.0D0) 1492 1493C Start timing. 1494C ------------- 1495 1496 TIMTOT = SECOND() 1497 1498C Calculate the (ai|ai) integral diagonal. 1499C ---------------------------------------- 1500 1501 TIMDIA = SECOND() 1502 CALL CHO_AIBJD_I(DIAG,CHAI,ISYM) 1503 TIMDIA = SECOND() - TIMDIA 1504 1505C Divide by the orbital energy denominators. 1506C t(ai,ai) = (ai|ai)/[2*e(a)-2*e(i)]. 1507C ------------------------------------------ 1508 1509 TIMDIV = SECOND() 1510 DO ISYMI = 1,NSYM 1511 1512 ISYMA = MULD2H(ISYMI,ISYM) 1513 1514 DO I = 1,NRHF(ISYMI) 1515 1516 KOFFI = IRHF(ISYMI) + I 1517 1518 DO A = 1,NVIR(ISYMA) 1519 1520 KOFFA = IVIR(ISYMA) + A 1521 1522 AI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 1523 1524 DENOM = TWO*(FOCKD(KOFFA)-FOCKD(KOFFI)) 1525 DIAG(AI) = DIAG(AI)/DENOM 1526 1527 ENDDO 1528 1529 ENDDO 1530 1531 ENDDO 1532 TIMDIV = SECOND() - TIMDIV 1533 1534C Print. 1535C ------ 1536 1537 IF (IPRINT .GE. INFO) THEN 1538 TIMTOT = SECOND() - TIMTOT 1539 XNRM = DSQRT(DDOT(NT1AM(ISYM),DIAG,1,DIAG,1)) 1540 CALL HEADER('Information from '//SECNAM,-1) 1541 WRITE(LUPRI,'(5X,A,I1)') 1542 & 'CC2 amplitude diagonal calculated, symmetry block:',ISYM 1543 WRITE(LUPRI,'(5X,A,1P,D30.15,/)') 1544 & 'Norm of calculated diagonal: ',XNRM 1545 WRITE(LUPRI,'(5X,A,F10.2,A)') 1546 & 'Time for calculating (ai|ai) diagonal: ',TIMDIA,' seconds' 1547 WRITE(LUPRI,'(5X,A,F10.2,A)') 1548 & 'Time for orbital energy denominators : ',TIMDIV,' seconds' 1549 WRITE(LUPRI,'(5X,A)') 1550 & '---------------------------------------------------------' 1551 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 1552 & 'Total time : ',TIMTOT,' seconds' 1553 ENDIF 1554 1555 RETURN 1556 END 1557C /* Deck cc2_choam */ 1558 SUBROUTINE CC2_CHOAM(T2AM,FOCKD,WORK,LWORK,LISTBJ,NUMBJ,ISYMBJ, 1559 & LUCHMO) 1560C 1561C Thomas Bondo Pedersen, August 2002. 1562C 1563C Purpose: 1564C Calculate minus CC2 amplitudes for a list, LISTBJ, of NUMBJ column 1565C (bj) indices of symmetry ISYMBJ. 1566C 1567C Notes: 1568C (a) MO Cholesky vectors, L(ai,J), are assumed available on disk on 1569C opened unit LUCHMO (in subroutine CHO_AIBJ). 1570C 1571#include "implicit.h" 1572 DIMENSION T2AM(*), FOCKD(*), WORK(LWORK) 1573 INTEGER LISTBJ(*) 1574#include "ccorb.h" 1575#include "ccsdsym.h" 1576#include "ccsdinp.h" 1577#include "priunit.h" 1578 1579 INTEGER AI,BJ,AIBJ 1580 1581 CHARACTER*9 SECNAM 1582 PARAMETER (SECNAM = 'CC2_CHOAM') 1583 1584 PARAMETER (INFO = 20) 1585 PARAMETER (ZERO = 0.0D0) 1586 1587C Start timing. 1588C ------------- 1589 1590 TIMTOT = SECOND() 1591 TIMINV = ZERO 1592 1593C Calculate (ai|#[bj]) integrals. 1594C ------------------------------- 1595 1596 TIMINT = SECOND() 1597 CALL CHO_AIBJ(T2AM,WORK,LWORK,LISTBJ,NUMBJ,ISYMBJ,LUCHMO) 1598 TIMINT = SECOND() - TIMINT 1599 1600C Divide by the orbital energy denominators. 1601C t(ai,bj) = (ai|bj)/[e(a)-e(i)+e(b)-e(j)]. 1602C ------------------------------------------ 1603 1604 TIMDIV = SECOND() 1605 ISYMAI = ISYMBJ 1606 DO IBJ = 1,NUMBJ 1607 1608 DTIME = SECOND() 1609 BJ = LISTBJ(IBJ) 1610 CALL CC2_INVBJ(BJ,ISYMBJ,B,ISYMB,J,ISYMJ) 1611 DTIME = SECOND() - DTIME 1612 TIMINV = TIMINV + DTIME 1613 1614 KOFFJ = IRHF(ISYMJ) + J 1615 KOFFB = IVIR(ISYMB) + B 1616 DENBJ = FOCKD(KOFFB) - FOCKD(KOFFJ) 1617 1618 DO ISYMI = 1,NSYM 1619 1620 ISYMA = MULD2H(ISYMI,ISYMAI) 1621 1622 DO I = 1,NRHF(ISYMI) 1623 1624 KOFFI = IRHF(ISYMI) + I 1625 1626 DO A = 1,NVIR(ISYMA) 1627 1628 KOFFA = IVIR(ISYMA) + A 1629 1630 AI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 1631 AIBJ = NT1AM(ISYMAI)*(IBJ - 1) + AI 1632 1633 DENOM = DENBJ + FOCKD(KOFFA) - FOCKD(KOFFI) 1634 T2AM(AIBJ) = T2AM(AIBJ)/DENOM 1635 1636 ENDDO 1637 1638 ENDDO 1639 1640 ENDDO 1641 1642 ENDDO 1643 TIMDIV = SECOND() - TIMDIV 1644 1645C Print. 1646C ------ 1647 1648 IF (IPRINT .GE. INFO) THEN 1649 TIMTOT = SECOND() - TIMTOT 1650 LEN = NT1AM(ISYMBJ)*NUMBJ 1651 XNRM = DSQRT(DDOT(LEN,T2AM,1,T2AM,1)) 1652 CALL HEADER('Information from '//SECNAM,-1) 1653 WRITE(LUPRI,'(5X,A,I1)') 1654 & 'CC2 amplitudes calculated, symmetry block:',ISYMBJ 1655 WRITE(LUPRI,'(5X,A,I8,A)') 1656 & 'CC2 amplitudes calculated for the following ',NUMBJ, 1657 & ' columns:' 1658 WRITE(LUPRI,'(10I8)') (LISTBJ(I),I=1,NUMBJ) 1659 WRITE(LUPRI,'(5X,A,1P,D30.15,/)') 1660 & 'Norm of calculated amplitudes: ',XNRM 1661 WRITE(LUPRI,'(5X,A,F10.2,A)') 1662 & 'Time for calculating (ai|bj) integrals: ',TIMINT,' seconds' 1663 WRITE(LUPRI,'(5X,A,F10.2,A)') 1664 & 'Time for orbital energy denominators : ',TIMDIV,' seconds' 1665 WRITE(LUPRI,'(5X,A,F10.2,A)') 1666 & ' - of which CC2_INVBJ required : ',TIMINV,' seconds' 1667 WRITE(LUPRI,'(5X,A)') 1668 & '----------------------------------------------------------' 1669 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 1670 & 'Total time : ',TIMTOT,' seconds' 1671 ENDIF 1672 1673 RETURN 1674 END 1675C /* Deck cc2_choam_i */ 1676 SUBROUTINE CC2_CHOAM_I(T2AM,FOCKD,CHAI,CHBJ,LISTBJ,NUMBJ,ISYMBJ) 1677C 1678C Thomas Bondo Pedersen, September 2002. 1679C 1680C Purpose: 1681C Calculate minus CC2 amplitudes for a list, LISTBJ, of NUMBJ column 1682C (bj) indices of symmetry ISYMBJ. 1683C 1684C Notes: 1685C (a) MO Cholesky vectors, L(ai,J), are assumed stored in CHAI. 1686C 1687#include "implicit.h" 1688 DIMENSION T2AM(*), FOCKD(*), CHAI(*), CHBJ(*) 1689 INTEGER LISTBJ(*) 1690#include "ccorb.h" 1691#include "ccsdsym.h" 1692#include "ccsdinp.h" 1693#include "priunit.h" 1694 1695 INTEGER AI,BJ,AIBJ 1696 1697 CHARACTER*11 SECNAM 1698 PARAMETER (SECNAM = 'CC2_CHOAM_I') 1699 1700 PARAMETER (INFO = 20) 1701 PARAMETER (ZERO = 0.0D0) 1702 1703C Start timing. 1704C ------------- 1705 1706 TIMTOT = SECOND() 1707 TIMINV = ZERO 1708 1709C Calculate (ai|#[bj]) integrals. 1710C ------------------------------- 1711 1712 TIMINT = SECOND() 1713 CALL CHO_AIBJ_I(T2AM,CHAI,CHBJ,LISTBJ,NUMBJ,ISYMBJ) 1714 TIMINT = SECOND() - TIMINT 1715 1716C Divide by the orbital energy denominators. 1717C t(ai,bj) = (ai|bj)/[e(a)-e(i)+e(b)-e(j)]. 1718C ------------------------------------------ 1719 1720 TIMDIV = SECOND() 1721 ISYMAI = ISYMBJ 1722 DO IBJ = 1,NUMBJ 1723 1724 DTIME = SECOND() 1725 BJ = LISTBJ(IBJ) 1726 CALL CC2_INVBJ(BJ,ISYMBJ,B,ISYMB,J,ISYMJ) 1727 DTIME = SECOND() - DTIME 1728 TIMINV = TIMINV + DTIME 1729 1730 KOFFJ = IRHF(ISYMJ) + J 1731 KOFFB = IVIR(ISYMB) + B 1732 DENBJ = FOCKD(KOFFB) - FOCKD(KOFFJ) 1733 1734 DO ISYMI = 1,NSYM 1735 1736 ISYMA = MULD2H(ISYMI,ISYMAI) 1737 1738 DO I = 1,NRHF(ISYMI) 1739 1740 KOFFI = IRHF(ISYMI) + I 1741 1742 DO A = 1,NVIR(ISYMA) 1743 1744 KOFFA = IVIR(ISYMA) + A 1745 1746 AI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 1747 AIBJ = NT1AM(ISYMAI)*(IBJ - 1) + AI 1748 1749 DENOM = DENBJ + FOCKD(KOFFA) - FOCKD(KOFFI) 1750 T2AM(AIBJ) = T2AM(AIBJ)/DENOM 1751 1752 ENDDO 1753 1754 ENDDO 1755 1756 ENDDO 1757 1758 ENDDO 1759 TIMDIV = SECOND() - TIMDIV 1760 1761C Print. 1762C ------ 1763 1764 IF (IPRINT .GE. INFO) THEN 1765 TIMTOT = SECOND() - TIMTOT 1766 LEN = NT1AM(ISYMBJ)*NUMBJ 1767 XNRM = DSQRT(DDOT(LEN,T2AM,1,T2AM,1)) 1768 CALL HEADER('Information from '//SECNAM,-1) 1769 WRITE(LUPRI,'(5X,A,I1)') 1770 & 'CC2 amplitudes calculated, symmetry block:',ISYMBJ 1771 WRITE(LUPRI,'(5X,A,I8,A)') 1772 & 'CC2 amplitudes calculated for the following ',NUMBJ, 1773 & ' columns:' 1774 WRITE(LUPRI,'(10I8)') (LISTBJ(I),I=1,NUMBJ) 1775 WRITE(LUPRI,'(5X,A,1P,D30.15,/)') 1776 & 'Norm of calculated amplitudes: ',XNRM 1777 WRITE(LUPRI,'(5X,A,F10.2,A)') 1778 & 'Time for calculating (ai|bj) integrals: ',TIMINT,' seconds' 1779 WRITE(LUPRI,'(5X,A,F10.2,A)') 1780 & 'Time for orbital energy denominators : ',TIMDIV,' seconds' 1781 WRITE(LUPRI,'(5X,A,F10.2,A)') 1782 & ' - of which CC2_INVBJ required : ',TIMINV,' seconds' 1783 WRITE(LUPRI,'(5X,A)') 1784 & '----------------------------------------------------------' 1785 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 1786 & 'Total time : ',TIMTOT,' seconds' 1787 ENDIF 1788 1789 RETURN 1790 END 1791C /* Deck cc2_invbj */ 1792 SUBROUTINE CC2_INVBJ(BJ,ISYMBJ,B,ISYMB,J,ISYMJ) 1793C 1794C Thomas Bondo Pedersen, August 2002. 1795C 1796C Purpose: "Invert" the compund index BJ of symmetry ISYMBJ 1797C to get B,ISYMB and J,ISYMJ. 1798C 1799#include "implicit.h" 1800 INTEGER BJ 1801#include "ccorb.h" 1802#include "ccsdsym.h" 1803#include "priunit.h" 1804 1805 LOGICAL LOCDBG 1806 PARAMETER (LOCDBG = .FALSE.) 1807 1808 CHARACTER*9 SECNAM 1809 PARAMETER (SECNAM = 'CC2_INVBJ') 1810 1811C Check for errors. 1812C ----------------- 1813 1814 IF ((BJ.LE.0) .OR. (BJ.GT.NT1AM(ISYMBJ))) THEN 1815 WRITE(LUPRI,'(//,5X,A,A,A)') 1816 & 'BJ out of bounds in ',SECNAM,':' 1817 WRITE(LUPRI,'(5X,A,I10,A,I1)') 1818 & 'BJ = ',BJ,', ISYMBJ = ',ISYMBJ 1819 IF ((ISYMBJ.LE.0) .OR. (ISYMBJ.GT.NSYM)) THEN 1820 WRITE(LUPRI,'(5X,A)') 1821 & '(The symmetry is out of bounds)' 1822 ENDIF 1823 WRITE(LUPRI,*) 1824 CALL QUIT('Error in '//SECNAM) 1825 ENDIF 1826 1827C Invert. 1828C Detailed search only in relevant symmetry block. 1829C ------------------------------------------------ 1830 1831 DO ISYMO = 1,NSYM 1832 1833 ISYMV = MULD2H(ISYMO,ISYMBJ) 1834 1835 IBJF = IT1AM(ISYMV,ISYMO) + 1 1836 IBJL = IBJF + NVIR(ISYMV)*NRHF(ISYMO) - 1 1837 1838 IF ((BJ.GE.IBJF) .AND. (BJ.LE.IBJL)) THEN 1839 1840 ICOUNT = IBJF - 1 1841 1842 DO IJ = 1,NRHF(ISYMO) 1843 DO IB = 1,NVIR(ISYMV) 1844 1845 ICOUNT = ICOUNT + 1 1846 1847 IF (ICOUNT .EQ. BJ) THEN 1848 J = IJ 1849 ISYMJ = ISYMO 1850 B = IB 1851 ISYMB = ISYMV 1852 GOTO 100 1853 ENDIF 1854 1855 ENDDO 1856 ENDDO 1857 1858 ENDIF 1859 1860 ENDDO 1861 1862C Debug: print. 1863C ------------- 1864 1865 100 IF (LOCDBG) THEN 1866 WRITE(LUPRI,'(/,5X,A,A)') 1867 & SECNAM,': Debug flag set. Testing bj -> b,j:' 1868 WRITE(LUPRI,'(5X,A,A,I10,1X,I1)') 1869 & SECNAM,': Input : BJ,ISYMBJ: ',BJ,ISYMBJ 1870 WRITE(LUPRI,'(5X,A,A,I10,1X,I1)') 1871 & SECNAM,': Output: B,ISYMB : ',B,ISYMB 1872 WRITE(LUPRI,'(5X,A,A,I10,1X,I1)') 1873 & SECNAM,': Output: J,ISYMJ : ',J,ISYMJ 1874 IERRS = 0 1875 IERRC = 0 1876 ISYMVO = MULD2H(ISYMB,ISYMJ) 1877 IF (ISYMVO .NE. ISYMBJ) IERRS = 1 1878 IBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 1879 IF (IBJ .NE. BJ) IERRC = 1 1880 IERR = IERRS + IERRC 1881 IF (IERR .NE. 0) THEN 1882 WRITE(LUPRI,'(5X,A,A,I1,A)') 1883 & SECNAM,': ',IERR,' error(s) detected.' 1884 IF (IERRS .NE. 0) THEN 1885 WRITE(LUPRI,'(5X,A,A)') 1886 & SECNAM,': Symmetry error detected.' 1887 ENDIF 1888 IF (IERRC .NE. 0) THEN 1889 WRITE(LUPRI,'(5X,A,A)') 1890 & SECNAM,': Index error detected.' 1891 ENDIF 1892 WRITE(LUPRI,*) 1893 CALL QUIT('Error detected in '//SECNAM) 1894 ELSE 1895 WRITE(LUPRI,'(5X,A,A,/)') 1896 & SECNAM,': No errors detected.' 1897 ENDIF 1898 ENDIF 1899 1900 RETURN 1901 END 1902C /* Deck chocc2_yi */ 1903 SUBROUTINE CHOCC2_YI(FOCKD,OMEGA1,WORK,LWORK,T2NRM,ECC2) 1904C 1905C Thomas Bondo Pedersen, August 2002. 1906C 1907C Purpose: Calculate Y-intermediates and the I-term: 1908C 1909C Y(ai,J) = sum(bj) s(ai,bj) * L(jb,J) 1910C 1911C OMEGA1(ai) = OMEGA1(ai) 1912C + sum(bj) s(ai,bj) * F(jb) 1913C 1914C where s(ai,bj) = 2*t(ai,bj) - t(aj,bi), and 1915C t(ai,bj) = -(ai|bj)/[e(a)-e(i)+e(b)-e(j)]. 1916C 1917C Calculate contribution to energy from doubles 1918C amplitudes: 1919C 1920C ECC2 = ECC2 + sum(aiJ) Y(ai,J) * L(ia,J) 1921C 1922C Notes: 1923C (1) The doubles amplitudes are never stored in full, neither on 1924C disk nor in core. There are two main routes, controlled by 1925C flag CHOT2: 1926C (a) calculation from (ai|bj) integrals (CHOT2=.FALSE.): 1927C - represented by the transformed AO Cholesky vectors, OR 1928C - represented by vectors from separately Cholesky 1929C decomposed (ai|bj) integrals. 1930C (b) decomposition of the amplitudes themselves (CHOT2=.TRUE.). 1931C (2) The norm of the packed amplitudes is returned in T2NRM, 1932C partly for debugging and partly for mimicking the output 1933C of the "standard" code. 1934C 1935#include "implicit.h" 1936 DIMENSION FOCKD(*), OMEGA1(*), WORK(LWORK) 1937#include "maxorb.h" 1938#include "ccdeco.h" 1939#include "ccorb.h" 1940#include "ccsdsym.h" 1941#include "ccsdinp.h" 1942#include "chocc2.h" 1943#include "priunit.h" 1944#include "iratdef.h" 1945#include "dummy.h" 1946 1947 CHARACTER*9 SECNAM 1948 PARAMETER (SECNAM = 'CHOCC2_YI') 1949 1950 PARAMETER (ZERO = 0.0D0) 1951 1952C Decomposition as requested from input. 1953C -------------------------------------- 1954 1955 IF (CHOT2 .OR. CHOMO2) THEN 1956 1957C Set threshold and span factor. 1958C ------------------------------ 1959 1960 IF (THRCC2 .LT. ZERO) THRCC2 = THRCOM 1961 IF (SPACC2 .LT. ZERO) SPACC2 = SPAN 1962 1963C Set requested decomposition. 1964C Amplitude deco. precedes integral deco. 1965C --------------------------------------- 1966 1967 IF (CHOT2) THEN 1968 IMAT = 3 1969 ELSE IF (CHOMO2) THEN 1970 IMAT = 2 1971 ENDIF 1972 1973C Scratch allocations. 1974C -------------------- 1975 1976 KDIANL = 1 1977 KCOLUM = KDIANL + 3*NSYM 1978 KINDIA = KCOLUM + (NSYM - 1)/IRAT + 1 1979 KEND1 = KINDIA + (MXDEC2 - 1)/IRAT + 1 1980 LWRK1 = LWORK - KEND1 1981 1982 IF (LWRK1 .LT. 0) THEN 1983 WRITE(LUPRI,'(//,5X,A,A,A)') 1984 & 'Insufficient memory in ',SECNAM,' - allocation: decomp.' 1985 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 1986 & 'Need (more than): ',KEND1, 1987 & 'Available : ',LWORK 1988 CALL QUIT('Insufficient memory in '//SECNAM) 1989 ENDIF 1990 1991C Decompose. 1992C ---------- 1993 1994 CALL CC_DECMO(WORK(KDIANL),NCHMO2,WORK(KCOLUM),WORK(KINDIA), 1995 & THRCC2,SPACC2, 1996 & FCC2S,.FALSE.,MXDEC2,NCHRD2,NT1AM,IPRINT,IMAT, 1997 & NSYM,THZCC2,FOCKD,WORK(KEND1),LWRK1) 1998 1999 ENDIF 2000 2001C Calculate Y-intermediates, I-term and energy contribution. 2002C ---------------------------------------------------------- 2003 2004 IF (IALCC2 .EQ. 1) THEN 2005 2006 CALL CC2_CHOYI1(FOCKD,OMEGA1,WORK,LWORK,T2NRM,ECC2) 2007 2008 ELSE IF (IALCC2 .EQ. 2) THEN 2009 2010 KFCKIA = 1 2011 KDUM1 = KFCKIA + NT1AM(1) 2012 KDUM2 = KDUM1 + 1 2013 KEND1 = KDUM2 + 1 2014 LWRK1 = LWORK - KEND1 + 1 2015 2016 IF (LWRK1 .LT. 0) THEN 2017 WRITE(LUPRI,'(//,5X,A,A,A)') 2018 & 'Insufficient memory in ',SECNAM, 2019 & ' - allocation: F(ia)' 2020 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2021 & 'Need (more than): ',KEND1-1, 2022 & 'Available : ',LWORK 2023 CALL QUIT('Insufficient memory in '//SECNAM) 2024 ENDIF 2025 2026 CALL ONEL_OP(-1,3,LUFIA) 2027 CALL CHO_MOREAD(WORK(KFCKIA),NT1AM(1),1,1,LUFIA) 2028 CALL ONEL_OP(1,3,LUFIA) 2029 2030 CALL CC_CYIV0(FOCKD,OMEGA1,WORK(KFCKIA),WORK(KEND1),LWRK1, 2031 & T2NRM,T2CNM,.FALSE.,.FALSE.) 2032 CALL CC_CYIENR(WORK,LWORK,ECC2) 2033 2034 ELSE 2035 2036 WRITE(LUPRI,'(//,5X,A,A)') 2037 & 'Error in ',SECNAM 2038 WRITE(LUPRI,'(5X,A,I10,/)') 2039 & 'Algorithm option not recognized: IALCC2 = ',IALCC2 2040 CALL QUIT('Error in '//SECNAM) 2041 2042 ENDIF 2043 2044 RETURN 2045 END 2046C /* Deck cc_cyiv0 */ 2047 SUBROUTINE CC_CYIV0(FOCKD,OMEGA1,FIA,WORK,LWORK,T2NRM,T2CNM, 2048 & DELP1,DELP2) 2049C 2050C Thomas Bondo Pedersen, February 2003. 2051C 2052C Purpose: Calculate Y intermediates and the I term for the CC2 2053C vector function. 2054C 2055#include "implicit.h" 2056 DIMENSION FOCKD(*), OMEGA1(*), FIA(*), WORK(LWORK) 2057 LOGICAL DELP1, DELP2 2058#include "maxorb.h" 2059#include "ccdeco.h" 2060#include "ccorb.h" 2061#include "ccsdsym.h" 2062#include "ccsdinp.h" 2063#include "priunit.h" 2064#include "chocc2.h" 2065 2066 INTEGER NCHP12(8) 2067 2068C Set Cholesky vector type for amplitude assembly. 2069C ------------------------------------------------ 2070 2071 FAC1 = 1.0D0 2072 FAC2 = 0.0D0 2073 NUMP12 = 1 2074 NUMUV = 0 2075 ISYMP1 = 1 2076 ISYMP2 = ISYMP1 2077 SCD = 1.0D0 2078 SCDG = 1.0D0 2079 IOPTDN = 1 2080 FREQ = 0.0D0 2081 IOPTCE = 1 2082 IF (CHOT2) THEN 2083 JTYP1 = 6 2084 IOPTDN = 0 2085 DO ISYCHO = 1,NSYM 2086 NCHP12(ISYCHO) = NCHMO2(ISYCHO) 2087 ENDDO 2088 ELSE IF (CHOMO2) THEN 2089 JTYP1 = 5 2090 DO ISYCHO = 1,NSYM 2091 NCHP12(ISYCHO) = NCHMO2(ISYCHO) 2092 ENDDO 2093 ELSE 2094 JTYP1 = 3 2095 DO ISYCHO = 1,NSYM 2096 NCHP12(ISYCHO) = NUMCHO(ISYCHO) 2097 ENDDO 2098 ENDIF 2099 JTYP2 = JTYP1 2100 2101C Set info for Y intermediates. 2102C ----------------------------- 2103 2104 NUMZ = 1 2105 JTYPZ = 1 2106 ISYMZ = 1 2107 JTYPY = 0 2108 2109C Calculate Y intermediates and I-term. 2110C ------------------------------------- 2111 2112 CALL CC_CYI(JTYP1,ISYMP1,JTYP2,ISYMP2,NCHP12,FAC1,NUMP12, 2113 & DUMMY,IDUMMY,DUMMY,IDUMMY,FAC2,NUMUV, 2114 & SCD,SCDG,IOPTDN,FOCKD,FREQ,IOPTCE, 2115 & JTYPZ,ISYMZ,NUMCHO,NUMZ,JTYPY, 2116 & FIA,1,1,OMEGA1, 2117 & WORK,LWORK,T2NRM,T2CNM,DELP1,DELP2) 2118 2119 RETURN 2120 END 2121C /* Deck cc2_choyi1 */ 2122 SUBROUTINE CC2_CHOYI1(FOCKD,OMEGA1,WORK,LWORK,T2NRM,ECC2) 2123C 2124C Thomas Bondo Pedersen, August 2002. 2125C 2126C Purpose: Calculate Y-intermediates and the I-term. 2127C Calculate energy contribution from doubles 2128C amplitudes. 2129C 2130#include "implicit.h" 2131 DIMENSION FOCKD(*), OMEGA1(*), WORK(LWORK) 2132#include "maxorb.h" 2133#include "ccdeco.h" 2134#include "ccisvi.h" 2135#include "priunit.h" 2136#include "ccorb.h" 2137#include "ccsdsym.h" 2138#include "ccsdinp.h" 2139#include "chocc2.h" 2140 2141 INTEGER ABEG, AEND 2142 INTEGER LVIR(8), IOFA1(8), IX1AM(8,8), NX1AM(8), IX2SQ(8) 2143 2144 CHARACTER*10 SECNAM 2145 PARAMETER (SECNAM = 'CC2_CHOYI1') 2146 2147 PARAMETER (INFO = 3) 2148 PARAMETER (IOPEN = -1, IKEEP = 1) 2149 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 2150 PARAMETER (WTOMB = 7.62939453125D-6, D100 = 100.00D0) 2151 2152C Initilize timings. 2153C ------------------ 2154 2155 TIMT = SECOND() 2156 TIMG = ZERO 2157 TIM2 = ZERO 2158 TIMI = ZERO 2159 TIMR = ZERO 2160 TIMY = ZERO 2161 2162C Initial check of memory. 2163C We must be able to hold at least 1 virtual in each symmetry, 2164C as well as at least 1 Cholesky vector in each symmetry. 2165C Note: the minimal requirements for amplitude generation, 2166C for the Y-intermediates, and for the I-term are identical. 2167C ------------------------------------------------------------ 2168 2169 NEED = 0 2170 DO ISYMA = 1,NSYM 2171 NEEDG = 0 2172 DO ISYCHO = 1,NSYM 2173 ISYMI = MULD2H(ISYMA,ISYCHO) 2174 NTAMG = NT1AM(ISYCHO) + NRHF(ISYMI) 2175 NEEDG = MAX(NEEDG,NTAMG) 2176 ENDDO 2177 NEEDA = NCKI(ISYMA) + NEEDG 2178 NEED = MAX(NEED,NEEDG) 2179 ENDDO 2180 IF (LWORK .LT. NEED) THEN 2181 WRITE(LUPRI,'(//,5X,A,A)') 2182 & 'Insufficient memory in ',SECNAM 2183 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2184 & 'Minimum memory required: ',NEED, 2185 & 'Available memory : ',LWORK 2186 CALL QUIT('Insufficient memory in '//SECNAM) 2187 ENDIF 2188 2189C Initializations. 2190C ---------------- 2191 2192 MXMEMT = 0 2193 MXMEML = 0 2194 XLWORK = ONE*LWORK 2195 2196C Determine memory split. 2197C ----------------------- 2198 2199 CALL MP2_MEMSPL2(LWORK,NUMCHO,LWRK) 2200 2201C Print. 2202C ------ 2203 2204 IF (IPRINT .GE. INFO) THEN 2205 CALL AROUND('Output from '//SECNAM) 2206 XMB = XLWORK*WTOMB 2207 LEFT = LWORK - LWRK 2208 WRITE(LUPRI,'(12X,A,/,12X,A,I10,A,F7.1,A)') 2209 & 'CC2 amplitudes will be calculated on the fly.', 2210 & 'Available memory: ',LWORK,' words (',XMB,' Mb).' 2211 WRITE(LUPRI,'(12X,A,1X,I10,A,/,12X,A,1X,I10,A)') 2212 & 'Maximum memory allowed for CC2 amplitudes: ',LWRK,' words.', 2213 & 'Maximum memory allowed for Cholesky vecs.: ',LEFT,' words.' 2214 WRITE(LUPRI,'(A)') ' ' 2215 WRITE(LUPRI,'(5X,A,A,/,5X,A,A)') 2216 & ' #a First Last Mem. Usage', 2217 & ' Time ', 2218 & ' (a, sym. a) (a, sym. a) (%) ', 2219 & ' (seconds)' 2220 WRITE(LUPRI,'(5X,A,A)') 2221 & '----------------------------------------------------', 2222 & '---------------' 2223 ENDIF 2224 2225C Start batched loop over a. 2226C -------------------------- 2227 2228 ABEG = 1 2229 NUMA = 0 2230 NPASS = 0 2231 100 CONTINUE 2232 2233 ABEG = ABEG + NUMA 2234 IF (ABEG .GT. NVIRT) GOTO 200 2235 2236C Time this batch. 2237C ---------------- 2238 2239 TABAT = SECOND() 2240 2241C Find out how many a's can be treated. 2242C ------------------------------------- 2243 2244 CALL GET_BTCH(LVIR,ABEG,NUMA,MEM,LWRK) 2245 2246C Check for errors. 2247C ----------------- 2248 2249 AEND = ABEG + NUMA - 1 2250 IF (AEND .GT. NVIRT) THEN 2251 WRITE(LUPRI,'(//,5X,A,A,A)') 2252 & 'Batch error in ',SECNAM,' !' 2253 WRITE(LUPRI,'(5X,A)') 2254 & '- dumping presumably most relevant info:' 2255 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10)') 2256 & 'First A : ',ABEG, 2257 & 'Last A : ',AEND, 2258 & 'Number of virtuals: ',NVIRT 2259 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2260 & 'Available memory : ',LWRK, 2261 & 'Needed memory : ',MEM 2262 CALL QUIT('Batch error in '//SECNAM) 2263 ENDIF 2264 2265 IF (NUMA .LE. 0) THEN 2266 WRITE(LUPRI,'(//,5X,A,A)') 2267 & 'Insufficient memory for a-batch in ',SECNAM 2268 WRITE(LUPRI,'(5X,A,I10,/)') 2269 & 'Available memory: ',LWRK 2270 CALL QUIT('Insufficient memory in '//SECNAM) 2271 ENDIF 2272 2273C Get symmetries of first and last a in batch. 2274C -------------------------------------------- 2275 2276 ISABEG = ISVI(ABEG) 2277 ISAEND = ISVI(AEND) 2278 2279C Complete allocation for amplitudes. 2280C ----------------------------------- 2281 2282 KT2AM = 1 2283 KEND1 = KT2AM + MEM 2284 LWRK1 = LWORK - KEND1 + 1 2285 2286 KLAST = KEND1 - 1 2287 MXMEML = MAX(MXMEML,KLAST) 2288 IF (MXMEML .GT. LWORK) THEN 2289 ID = 1 2290 GOTO 2000 2291 ENDIF 2292 2293C Set up local index arrays. 2294C -------------------------- 2295 2296 CALL IZERO(IOFA1,NSYM) 2297 2298 IA1 = ABEG 2299 DO ISYMA = ISABEG,ISAEND 2300 IOFA1(ISYMA) = IA1 + NRHFT - IVIR(ISYMA) 2301 IA1 = IA1 + LVIR(ISYMA) 2302 ENDDO 2303 2304 ICOUN2 = 0 2305 DO ISYMAI = 1,NSYM 2306 2307 ICOUN1 = 0 2308 DO ISYMI = 1,NSYM 2309 ISYMA = MULD2H(ISYMI,ISYMAI) 2310 IX1AM(ISYMA,ISYMI) = ICOUN1 2311 ICOUN1 = ICOUN1 + LVIR(ISYMA)*NRHF(ISYMI) 2312 ENDDO 2313 NX1AM(ISYMAI) = ICOUN1 2314 2315 IX2SQ(ISYMAI) = ICOUN2 2316 ICOUN2 = ICOUN2 + NT1AM(ISYMAI)*NX1AM(ISYMAI) 2317 2318 ENDDO 2319 NX2SQ = ICOUN2 2320 2321 IF (NX2SQ .NE. MEM) THEN 2322 WRITE(LUPRI,'(//,5X,A,A)') 2323 & 'Error detected in ',SECNAM 2324 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,/)') 2325 & 'NX2SQ is calculated to be ',NX2SQ, 2326 & 'MEM determines alloc. as ',MEM, 2327 & 'These numbers must be equal....' 2328 CALL QUIT('Error in '//SECNAM) 2329 ENDIF 2330 2331C Calculate the CC2 amplitudes. 2332C Scale by -1 to get the right sign. 2333C ---------------------------------- 2334 2335 DTIME = SECOND() 2336 CALL CC2_CHOAMG(WORK(KT2AM),FOCKD,WORK(KEND1),LWRK1, 2337 & LVIR,IOFA1,IX1AM,NX1AM,IX2SQ,NX2SQ, 2338 & T2NRM,.TRUE.,MEMLEN) 2339 CALL DSCAL(NX2SQ,XMONE,WORK(KT2AM),1) 2340 DTIME = SECOND() - DTIME 2341 TIMG = TIMG + DTIME 2342 2343 KTEST = KLAST + MEMLEN 2344 MXMEML = MAX(MXMEML,KTEST) 2345 IF (MXMEML .GT. LWORK) THEN 2346 ID = 2 2347 GOTO 2000 2348 ENDIF 2349 2350C Set up 2 Coulomb minus exchange (almost) in place. 2351C -------------------------------------------------- 2352 2353 DTIME = SECOND() 2354 CALL CC2_TCME(WORK(KT2AM),WORK(KEND1),LWRK1,LVIR,IOFA1,IX1AM, 2355 & NX1AM,IX2SQ,NX2SQ,MEMLEN) 2356 DTIME = SECOND() - DTIME 2357 TIM2 = TIM2 + DTIME 2358 2359 KTEST = KLAST + MEMLEN 2360 MXMEML = MAX(MXMEML,KTEST) 2361 IF (MXMEML .GT. LWORK) THEN 2362 ID = 3 2363 GOTO 2000 2364 ENDIF 2365 2366C Calculate the I-term. 2367C --------------------- 2368 2369 DTIME = SECOND() 2370 CALL CC2_ITERM(OMEGA1,WORK(KT2AM),WORK(KEND1),LWRK1,LVIR,IOFA1, 2371 & IX1AM,NX1AM,IX2SQ,MEMLEN) 2372 DTIME = SECOND() - DTIME 2373 TIMI = TIMI + DTIME 2374 2375 KTEST = KLAST + MEMLEN 2376 MXMEML = MAX(MXMEML,KTEST) 2377 IF (MXMEML .GT. LWORK) THEN 2378 ID = 4 2379 GOTO 2000 2380 ENDIF 2381 2382C Calculate the Y-intermediates and energy contribution 2383C from doubles amplitudes. 2384C ----------------------------------------------------- 2385 2386 DTIME = SECOND() 2387 CALL CC2_Y(WORK(KT2AM),WORK(KEND1),LWRK1,LVIR,IOFA1,IX1AM, 2388 & NX1AM,IX2SQ,MEMLEN,ECC2) 2389 DTIME = SECOND() - DTIME 2390 TIMY = TIMY + DTIME 2391 2392 KTEST = KLAST + MEMLEN 2393 MXMEML = MAX(MXMEML,KTEST) 2394 IF (MXMEML .GT. LWORK) THEN 2395 ID = 5 2396 GOTO 2000 2397 ENDIF 2398 2399C Print information from this a-batch. 2400C ------------------------------------ 2401 2402 IF (IPRINT .GE. INFO) THEN 2403 TABAT = SECOND() - TABAT 2404 XMUSE = (D100*MXMEML)/XLWORK 2405 LAFST = IOFA1(ISABEG) 2406 LALST = IOFA1(ISAEND) + LVIR(ISAEND) - 1 2407 WRITE(LUPRI,'(5X,I6,4X,I6,1X,I1,8X,I6,1X,I1,9X,F6.2,8X,F10.2)') 2408 & NUMA,LAFST,ISABEG,LALST,ISAEND,XMUSE,TABAT 2409 ENDIF 2410 2411C Update memory info. 2412C ------------------- 2413 2414 MXMEMT = MAX(MXMEMT,MXMEML) 2415 MXMEML = 0 2416 2417C Update NPASS and go to next a-batch. 2418C (Which will exit immediately if no more a's.) 2419C --------------------------------------------- 2420 2421 NPASS = NPASS + 1 2422 GO TO 100 2423 2424C Exit a-batch: calculation done. 2425C ------------------------------- 2426 2427 200 CONTINUE 2428 2429C Print. 2430C ------ 2431 2432 IF (IPRINT .GE. INFO) THEN 2433 WRITE(LUPRI,'(5X,A,A,/)') 2434 & '----------------------------------------------------', 2435 & '---------------' 2436 TIMT = SECOND() - TIMT 2437 CALL HEADER('Timing of '//SECNAM,-1) 2438 WRITE(LUPRI,'(5X,A,I10,/)') 2439 & 'Passes through Cholesky file(s) : ',NPASS 2440 WRITE(LUPRI,'(5X,A,F10.2,A)') 2441 & 'Time used for generating amplitudes : ',TIMG,' seconds' 2442 WRITE(LUPRI,'(5X,A,F10.2,A)') 2443 & 'Time used for two Coulomb - exchange: ',TIM2,' seconds' 2444 WRITE(LUPRI,'(5X,A,F10.2,A)') 2445 & 'Time used for the I-term : ',TIMI,' seconds' 2446 WRITE(LUPRI,'(5X,A,F10.2,A)') 2447 & 'Time used for the Y-intermediates : ',TIMY,' seconds' 2448 WRITE(LUPRI,'(5X,A,F10.2,A)') 2449 & '--------------------------------------------------------' 2450 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 2451 & 'Total time : ',TIMT,' seconds' 2452 ENDIF 2453 2454C Normal exit. 2455C ------------ 2456 2457 RETURN 2458 2459C Memory error. 2460C ------------- 2461 2462 2000 WRITE(LUPRI,'(//,5X,A,A,A)') 2463 & 'Allocation error detected in ',SECNAM,':' 2464 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10)') 2465 & 'Last memory location used: ',MXMEML, 2466 & 'Memory available : ',LWORK 2467 WRITE(LUPRI,'(5X,A,I3,/)') 2468 & 'Identifier: ',ID 2469 CALL QUIT('Allocation error in '//SECNAM) 2470 2471 END 2472C /* Deck cc2_y */ 2473 SUBROUTINE CC2_Y(T2AM,WORK,LWORK,LVIR,IOFA1,IX1AM,NX1AM,IX2SQ, 2474 & MEMLEN,ECC2) 2475C 2476C Thomas Bondo Pedersen, August 2002. 2477C 2478C Purpose: Calculate the Y-intermediates and write them to disk. 2479C Calculate contribution to energy from doubles. 2480C 2481C Note: The Y-intermediate I/O section is probably crappy.... 2482C 2483#include "implicit.h" 2484 DIMENSION T2AM(*), WORK(LWORK) 2485 INTEGER LVIR(8), IOFA1(8), IX1AM(8,8), NX1AM(8), IX2SQ(8) 2486#include "maxorb.h" 2487#include "ccdeco.h" 2488#include "ccorb.h" 2489#include "ccsdsym.h" 2490#include "chocc2.h" 2491#include "ccsdinp.h" 2492#include "priunit.h" 2493 2494 INTEGER LROW, ICOL, ADDR 2495 2496 CHARACTER*5 SECNAM 2497 PARAMETER (SECNAM = 'CC2_Y') 2498 2499 LOGICAL LOCDBG 2500 PARAMETER (LOCDBG = .FALSE.) 2501 2502 PARAMETER (IOPEN = -1, IKEEP = 1, ITYP = 1, INFO = 5) 2503 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 2504 2505C Initialize timings. 2506C ------------------- 2507 2508 TIMTOT = SECOND() 2509 TIMRMO = ZERO 2510 TIMYCL = ZERO 2511 TIMENR = ZERO 2512 TIMWRY = ZERO 2513 2514C Initialize MEMLEN. 2515C ------------------ 2516 2517 MEMLEN = 0 2518 2519C Start Cholesky symmetry loop. 2520C ----------------------------- 2521 2522 DO ISYCHO = 1,NSYM 2523 2524 IF (NUMCHO(ISYCHO) .LE. 0) GOTO 1000 2525 IF (NT1AM(ISYCHO) .LE. 0) GOTO 1000 2526 IF (NX1AM(ISYCHO) .LE. 0) GOTO 1000 2527 2528C Open MO file for this symmetry [L(jb,J)]. 2529C ----------------------------------------- 2530 2531 DTIME = SECOND() 2532 CALL CHO_MOP(IOPEN,ITYP,ISYCHO,LUCHMO,1,1) 2533 DTIME = SECOND() - DTIME 2534 TIMRMO = TIMRMO + DTIME 2535 2536C Open Y-intermediate file for this symmetry. 2537C ------------------------------------------- 2538 2539 DTIME = SECOND() 2540 CALL WOPEN2(LCC2YM,FCC2YM(ISYCHO),64,0) 2541 DTIME = SECOND() - DTIME 2542 TIMWRY = TIMWRY + DTIME 2543 2544C Set up batch. 2545C ------------- 2546 2547 MINMEM = NT1AM(ISYCHO) + NX1AM(ISYCHO) 2548 NVEC = MIN(LWORK/MINMEM,NUMCHO(ISYCHO)) 2549 2550 IF (NVEC .LE. 0) THEN 2551 WRITE(LUPRI,'(//,5X,A,A)') 2552 & 'Insufficient memory for batch in ',SECNAM 2553 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2554 & 'Need (pref. multiple) : ',MINMEM, 2555 & 'Total memory available: ',LWORK 2556 CALL QUIT('Insufficient memory in '//SECNAM) 2557 ENDIF 2558 2559 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 2560 2561C Batch loop. 2562C ----------- 2563 2564 DO IBATCH = 1,NBATCH 2565 2566 NUMV = NVEC 2567 IF (IBATCH .EQ. NBATCH) THEN 2568 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 2569 ENDIF 2570 2571 JVEC1 = NVEC*(IBATCH - 1) + 1 2572 2573C Complete allocation. 2574C -------------------- 2575 2576 KRESY = 1 2577 KCHOL = KRESY + NX1AM(ISYCHO)*NUMV 2578 KLAST = KCHOL + NT1AM(ISYCHO)*NUMV - 1 2579 2580 MEMLEN = MAX(MEMLEN,KLAST) 2581 2582C Read vectors. 2583C ------------- 2584 2585 DTIME = SECOND() 2586 CALL CHO_MOREAD(WORK(KCHOL),NT1AM(ISYCHO),NUMV,JVEC1,LUCHMO) 2587 DTIME = SECOND() - DTIME 2588 TIMRMO = TIMRMO + DTIME 2589 2590C Calculate Y-intermediates. 2591C -------------------------- 2592 2593 DTIME = SECOND() 2594 2595 KOFFT = IX2SQ(ISYCHO) + 1 2596 2597 NBJ = NT1AM(ISYCHO) 2598 NAI = NX1AM(ISYCHO) 2599 2600 CALL DGEMM('T','N',NAI,NUMV,NBJ, 2601 & ONE,T2AM(KOFFT),NBJ,WORK(KCHOL),NBJ, 2602 & ZERO,WORK(KRESY),NAI) 2603 2604 DTIME = SECOND() - DTIME 2605 TIMYCL = TIMYCL + DTIME 2606 2607C Calculate energy contribution from doubles. 2608C ------------------------------------------- 2609 2610 DTIME = SECOND() 2611 IF (NX1AM(ISYCHO) .EQ. NT1AM(ISYCHO)) THEN 2612 LTOT = NT1AM(ISYCHO)*NUMV 2613 ECC2 = ECC2 + DDOT(LTOT,WORK(KRESY),1,WORK(KCHOL),1) 2614 ELSE 2615 DO IVEC = 1,NUMV 2616 DO ISYMI = 1,NSYM 2617 ISYMA = MULD2H(ISYMI,ISYCHO) 2618 IF (LVIR(ISYMA) .GT. 0) THEN 2619 IF (LVIR(ISYMA) .EQ. NVIR(ISYMA)) THEN 2620 KOFF1 = KRESY + NX1AM(ISYCHO)*(IVEC - 1) 2621 & + IX1AM(ISYMA,ISYMI) 2622 KOFF2 = KCHOL + NT1AM(ISYCHO)*(IVEC - 1) 2623 & + IT1AM(ISYMA,ISYMI) 2624 LENAI = NVIR(ISYMA)*NRHF(ISYMI) 2625 ECC2 = ECC2 2626 & + DDOT(LENAI,WORK(KOFF1),1, 2627 & WORK(KOFF2),1) 2628 ELSE 2629 DO I = 1,NRHF(ISYMI) 2630 KOFF1 = KRESY + NX1AM(ISYCHO)*(IVEC - 1) 2631 & + IX1AM(ISYMA,ISYMI) 2632 & + LVIR(ISYMA)*(I - 1) 2633 KOFF2 = KCHOL + NT1AM(ISYCHO)*(IVEC - 1) 2634 & + IT1AM(ISYMA,ISYMI) 2635 & + NVIR(ISYMA)*(I - 1) 2636 & + IOFA1(ISYMA) - 1 2637 ECC2 = ECC2 2638 & + DDOT(LVIR(ISYMA),WORK(KOFF1),1, 2639 & WORK(KOFF2),1) 2640 ENDDO 2641 ENDIF 2642 ENDIF 2643 ENDDO 2644 ENDDO 2645 ENDIF 2646 DTIME = SECOND() - DTIME 2647 TIMENR = TIMENR + DTIME 2648 2649C Save Y-intermediates on disk. 2650C ----------------------------- 2651 2652 DTIME = SECOND() 2653 IF (NX1AM(ISYCHO) .EQ. NT1AM(ISYCHO)) THEN 2654 2655C Full vector; simply write. 2656C -------------------------- 2657 2658 ICOL = JVEC1 - 1 2659 LROW = NT1AM(ISYCHO) 2660 ADDR = LROW*ICOL + 1 2661 LEN = NT1AM(ISYCHO)*NUMV 2662 2663 IF (LOCDBG) THEN 2664 WRITE(LUPRI,'(5X,A,A)') 2665 & SECNAM,': debug: calling PUTWA2:' 2666 WRITE(LUPRI,'(7X,A,I1)') 2667 & '- full amplitude array of sym. ',ISYCHO 2668 ENDIF 2669 2670 CALL PUTWA2(LCC2YM,FCC2YM(ISYCHO),WORK(KRESY),ADDR,LEN) 2671 2672 ELSE IF (NX1AM(ISYCHO) .LT. NT1AM(ISYCHO)) THEN 2673 2674C Part of vector; write to disk in small chunks. 2675C ---------------------------------------------- 2676 2677 DO IVEC = 1,NUMV 2678 2679 JVEC = JVEC1 + IVEC - 1 2680 ICOL = JVEC - 1 2681 LROW = NT1AM(ISYCHO) 2682 2683 DO ISYMI = 1,NSYM 2684 2685 ISYMA = MULD2H(ISYMI,ISYCHO) 2686 2687 IF (LVIR(ISYMA) .GT. 0) THEN 2688 2689 IF (LVIR(ISYMA) .EQ. NVIR(ISYMA)) THEN 2690 2691 KOFFY = KRESY 2692 & + NX1AM(ISYCHO)*(IVEC - 1) 2693 & + IX1AM(ISYMA,ISYMI) 2694 ADDR = LROW*ICOL 2695 & + IT1AM(ISYMA,ISYMI) + 1 2696 2697 IF (LOCDBG) THEN 2698 WRITE(LUPRI,'(5X,A,A)') 2699 & SECNAM,': debug: calling PUTWA2:' 2700 WRITE(LUPRI,'(7X,A,I1,1X,I1)') 2701 & '- full symmetry block a,i: ',ISYMA,ISYMI 2702 ENDIF 2703 2704c write(lupri,*) 'LROW, ICOL: ',LROW,ICOL 2705c write(lupri,*) 'LROW*ICOL : ',LROW*ICOL 2706c write(lupri,*) 'IT1AM : ',IT1AM(ISYMA,ISYMI) 2707c write(lupri,*) 'ADDR expr.: ',LROW*ICOL+IT1AM(ISYMA,ISYMI) 2708c write(lupri,*) 'ADDR pass.: ',ADDR 2709 2710 CALL PUTWA2(LCC2YM,FCC2YM(ISYCHO), 2711 & WORK(KOFFY),ADDR, 2712 & NVIR(ISYMA)*NRHF(ISYMI)) 2713 2714 ELSE 2715 2716 DO I = 1,NRHF(ISYMI) 2717 2718 KOFFY = KRESY 2719 & + NX1AM(ISYCHO)*(IVEC - 1) 2720 & + IX1AM(ISYMA,ISYMI) 2721 & + LVIR(ISYMA)*(I - 1) 2722 ADDR = LROW*ICOL 2723 & + IT1AM(ISYMA,ISYMI) 2724 & + NVIR(ISYMA)*(I - 1) 2725 & + IOFA1(ISYMA) 2726 2727 IF (LOCDBG) THEN 2728 WRITE(LUPRI,'(5X,A,A)') 2729 & SECNAM,': debug: calling PUTWA2:' 2730 WRITE(LUPRI,'(7X,A,I1,1X,I1)') 2731 & '- symmetry block a,i: ',ISYMA,ISYMI 2732 WRITE(LUPRI,'(7X,A,I10,A,I10)') 2733 & '- writing ',LVIR(ISYMA),' for i = ',I 2734 ENDIF 2735 2736 CALL PUTWA2(LCC2YM,FCC2YM(ISYCHO), 2737 & WORK(KOFFY),ADDR, 2738 & LVIR(ISYMA)) 2739 2740 ENDDO 2741 2742 ENDIF 2743 2744 ENDIF 2745 2746 ENDDO 2747 2748 ENDDO 2749 2750 ELSE 2751 2752C Error. 2753C ------ 2754 2755 WRITE(LUPRI,'(//,5X,A,A,A)') 2756 & 'Error detected in ',SECNAM,':' 2757 WRITE(LUPRI,'(5X,A,I10,A,I1,A)') 2758 & 'Row dimension of Y-intermediate is ',NX1AM(ISYCHO), 2759 & ' (symmetry ',ISYCHO,')' 2760 WRITE(LUPRI,'(5X,A,I10)') 2761 & 'Maximum possible dimension is ',NT1AM(ISYCHO) 2762 WRITE(LUPRI,'(5X,A,/)') 2763 & 'The error is probably caused by calling routine.' 2764 CALL QUIT('Error detected in '//SECNAM) 2765 2766 ENDIF 2767 DTIME = SECOND() - DTIME 2768 TIMWRY = TIMWRY + DTIME 2769 2770 ENDDO 2771 2772C Close Y-intermediate file for this symmetry. 2773C -------------------------------------------- 2774 2775 DTIME = SECOND() 2776 CALL WCLOSE2(LCC2YM,FCC2YM(ISYCHO),'KEEP') 2777 DTIME = SECOND() - DTIME 2778 TIMWRY = TIMWRY + DTIME 2779 2780C Close MO file for this symmetry [L(jb,J)]. 2781C ------------------------------------------ 2782 2783 DTIME = SECOND() 2784 CALL CHO_MOP(IKEEP,ITYP,ISYCHO,LUCHMO,1,1) 2785 DTIME = SECOND() - DTIME 2786 TIMRMO = TIMRMO + DTIME 2787 2788 1000 CONTINUE 2789 2790 ENDDO 2791 2792C Print. 2793C ------ 2794 2795 IF (IPRINT .GE. INFO) THEN 2796 TIMTOT = SECOND() - TIMTOT 2797 CALL HEADER('Timing of '//SECNAM,-1) 2798 WRITE(LUPRI,'(5X,A,F10.2,A)') 2799 & 'Time used for I/O of MO Cholesky vectors : ',TIMRMO,' seconds' 2800 WRITE(LUPRI,'(5X,A,F10.2,A)') 2801 & 'Time used for I/O of the Y-intermediates: ',TIMWRY,' seconds' 2802 WRITE(LUPRI,'(5X,A,F10.2,A)') 2803 & 'Time used for calculating Y-intermediates: ',TIMYCL,' seconds' 2804 WRITE(LUPRI,'(5X,A,F10.2,A)') 2805 & 'Time used for calculating energy contr. : ',TIMENR,' seconds' 2806 WRITE(LUPRI,'(5X,A)') 2807 & '-------------------------------------------------------------' 2808 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 2809 & 'Total time : ',TIMTOT,' seconds' 2810 ENDIF 2811 2812 RETURN 2813 END 2814C /* Deck cc2_iterm */ 2815 SUBROUTINE CC2_ITERM(OMEGA1,T2AM,WORK,LWORK,LVIR,IOFA1,IX1AM, 2816 & NX1AM,IX2SQ,MEMLEN) 2817C 2818C Thomas Bondo Pedersen, August 2002. 2819C 2820C Purpose: Calculate the I-term from a batch of doubles amplitudes. 2821C 2822#include "implicit.h" 2823 DIMENSION OMEGA1(*), T2AM(*), WORK(LWORK) 2824 INTEGER LVIR(8), IOFA1(8), IX1AM(8,8), NX1AM(8), IX2SQ(8) 2825#include "ccorb.h" 2826#include "ccsdsym.h" 2827#include "priunit.h" 2828 2829 CHARACTER*9 SECNAM 2830 PARAMETER (SECNAM = 'CC2_ITERM') 2831 2832 PARAMETER (IOPEN = -1, IKEEP = 1, ITYP = 3) 2833 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 2834 2835C Initialize MEMLEN. 2836C ------------------ 2837 2838 MEMLEN = 0 2839 2840C Test if anything to do. 2841C ----------------------- 2842 2843 IF ((NX1AM(1).LE.0) .OR. (NT1AMX.LE.0)) RETURN 2844 2845C Allocation. 2846C ----------- 2847 2848 KFOCK = 1 2849 KSCR1 = KFOCK + NT1AMX 2850 KEND1 = KSCR1 + NX1AM(1) 2851 LWRK1 = LWORK - KEND1 + 1 2852 2853 KLAST = KEND1 - 1 2854 MEMLEN = MAX(MEMLEN,KLAST) 2855 2856 IF (LWRK1 .LT. 0) THEN 2857 WRITE(LUPRI,'(//,5X,A,A)') 2858 & 'Insufficient memory in ',SECNAM 2859 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2860 & 'Need : ',KEND1-1, 2861 & 'Available: ',LWORK 2862 CALL QUIT('Insufficient memory in '//SECNAM) 2863 ENDIF 2864 2865C Read Fock matrix block F(jb). 2866C ----------------------------- 2867 2868 CALL ONEL_OP(IOPEN,ITYP,LUFOCK) 2869 CALL CHO_MOREAD(WORK(KFOCK),NT1AMX,1,1,LUFOCK) 2870 CALL ONEL_OP(IKEEP,ITYP,LUFOCK) 2871 2872C Calculate contribution to I-term in scratch space. 2873C -------------------------------------------------- 2874 2875 KOFFT = IX2SQ(1) + 1 2876 CALL DGEMV('T',NT1AMX,NX1AM(1),ONE,T2AM(KOFFT),NT1AMX, 2877 & WORK(KFOCK),1,ZERO,WORK(KSCR1),1) 2878 2879C Add into OMEGA1. 2880C ---------------- 2881 2882 DO ISYMI = 1,NSYM 2883 2884 ISYMA = ISYMI 2885 2886 IF (LVIR(ISYMA) .LE. 0) GOTO 100 2887 2888 DO I = 1,NRHF(ISYMI) 2889 2890 KOFF1 = KSCR1 + IX1AM(ISYMA,ISYMI) + LVIR(ISYMA)*(I - 1) 2891 KOFF2 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) 2892 & + IOFA1(ISYMA) 2893 2894 CALL DAXPY(LVIR(ISYMA),ONE,WORK(KOFF1),1,OMEGA1(KOFF2),1) 2895 2896 ENDDO 2897 2898 100 CONTINUE 2899 2900 ENDDO 2901 2902 RETURN 2903 END 2904C /* Deck cc2_tcme */ 2905 SUBROUTINE CC2_TCME(T2AM,WORK,LWORK,LVIR,IOFB1,IX1AM,NX1AM,IX2SQ, 2906 & NX2SQ,MEMLEN) 2907C 2908C Thomas Bondo Pedersen, August 2002. 2909C 2910C Purpose: Set up two Coulomb minus exchange in place 2911C for a batch of CC2 amplitudes, t(ai,#bj). 2912C LWORK must be appr. 2*V. 2913C 2914C The overall structure is taken from the standard "exchange" 2915C CCSD_T2TP routine by A. Sanchez and H. Koch. 2916C 2917#include "implicit.h" 2918 DIMENSION T2AM(*), WORK(LWORK) 2919 INTEGER LVIR(8), IOFB1(8), IX1AM(8,8), NX1AM(8), IX2SQ(8) 2920#include "ccorb.h" 2921#include "ccsdsym.h" 2922#include "priunit.h" 2923 2924 CHARACTER*8 SECNAM 2925 PARAMETER (SECNAM = 'CC2_TCME') 2926 2927 LOGICAL LOCDBG 2928 PARAMETER (LOCDBG = .FALSE.) 2929 2930 PARAMETER (XMONE = -1.0D0, TWO = 2.0D0) 2931 PARAMETER (TOL = 1.0D-15) 2932 2933C Initialize MEMLEN. 2934C ------------------ 2935 2936 MEMLEN = 0 2937 2938C Debug: set up tcme in a straightforward fashion 2939C using (potentially a lot of) work space. 2940C ----------------------------------------------- 2941 2942 IF (LOCDBG) THEN 2943 MAXVIR = 0 2944 DO ISYM = 1,NSYM 2945 MAXVIR = MAX(MAXVIR,NVIR(ISYM)) 2946 ENDDO 2947 KSCRT = 2*MAXVIR + 1 2948 KEND0 = KSCRT + NX2SQ 2949 LWRK0 = LWORK - KEND0 + 1 2950 KLAST = KEND0 - 1 2951 MEMLEN = MAX(MEMLEN,KLAST) 2952 IF (LWRK0 .LT. 0) THEN 2953 WRITE(LUPRI,'(//,5X,A,A)') 2954 & 'Insufficient memory for debug in ',SECNAM 2955 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 2956 & 'Need : ',KEND0-1, 2957 & 'Available: ',LWORK 2958 CALL QUIT('Insufficient memory for debug in '//SECNAM) 2959 ENDIF 2960 DO ISYMBJ = 1,NSYM 2961 ISYMAI = ISYMBJ 2962 DO ISYMJ = 1,NSYM 2963 ISYMB = MULD2H(ISYMJ,ISYMBJ) 2964 DO J = 1,NRHF(ISYMJ) 2965 DO LB = 1,LVIR(ISYMB) 2966 LBJ = IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) 2967 & + LB 2968 DO ISYMI = 1,NSYM 2969 ISYMA = MULD2H(ISYMI,ISYMAI) 2970 ISYMAJ = MULD2H(ISYMA,ISYMJ) 2971 ISYMBI = MULD2H(ISYMB,ISYMI) 2972 DO I = 1,NRHF(ISYMI) 2973 LBI = IX1AM(ISYMB,ISYMI) 2974 & + LVIR(ISYMB)*(I - 1) + LB 2975 DO A = 1,NVIR(ISYMA) 2976 NAI = IT1AM(ISYMA,ISYMI) 2977 & + NVIR(ISYMA)*(I - 1) + A 2978 NAJ = IT1AM(ISYMA,ISYMJ) 2979 & + NVIR(ISYMA)*(J - 1) + A 2980 NAIBJ = IX2SQ(ISYMBJ) 2981 & + NT1AM(ISYMAI)*(LBJ - 1) + NAI 2982 NAJBI = IX2SQ(ISYMBI) 2983 & + NT1AM(ISYMAJ)*(LBI - 1) + NAJ 2984 WORK(KSCRT+NAIBJ-1) = TWO*T2AM(NAIBJ) 2985 & - T2AM(NAJBI) 2986 ENDDO 2987 ENDDO 2988 ENDDO 2989 ENDDO 2990 ENDDO 2991 ENDDO 2992 ENDDO 2993 ENDIF 2994 2995C Set up TCME. 2996C ------------ 2997 2998 DO ISYMJ = 1,NSYM 2999 DO J = 1,NRHF(ISYMJ) 3000 DO ISYMB = 1,NSYM 3001 3002 ISYMBJ = MULD2H(ISYMB,ISYMJ) 3003 ISYMAI = ISYMBJ 3004 3005 DO LB = 1,LVIR(ISYMB) 3006 3007 LBJ = IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) + LB 3008 3009 DO ISYMI = 1,ISYMJ 3010 3011 ISYMA = MULD2H(ISYMI,ISYMAI) 3012 ISYMAJ = MULD2H(ISYMA,ISYMJ) 3013 ISYMBI = MULD2H(ISYMB,ISYMI) 3014 3015 KAIBJ = 1 3016 KAJBI = KAIBJ + NVIR(ISYMA) 3017 KEND1 = KAJBI + NVIR(ISYMA) 3018 LWRK1 = LWORK - KEND1 + 1 3019 3020 KLAST = KEND1 - 1 3021 MEMLEN = MAX(MEMLEN,KLAST) 3022 3023 IF (LWRK1 .LT. 0) THEN 3024 WRITE(LUPRI,'(//,5X,A,A,/,5X,A,I10,/,5X,A,I10,/)') 3025 & 'Insufficient memory in ',SECNAM, 3026 & 'Need : ',KEND1-1, 3027 & 'Available: ',LWORK 3028 CALL QUIT('Insufficient memory in '//SECNAM) 3029 ENDIF 3030 3031 IF (ISYMI .EQ. ISYMJ) THEN 3032 NRHFI = J - 1 3033 ELSE 3034 NRHFI = NRHF(ISYMI) 3035 ENDIF 3036 3037 DO I = 1,NRHFI 3038 3039 LBI = IX1AM(ISYMB,ISYMI) + LVIR(ISYMB)*(I - 1) 3040 & + LB 3041 3042 NAIBJ = IX2SQ(ISYMAI) + NT1AM(ISYMAI)*(LBJ - 1) 3043 & + IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) 3044 & + 1 3045 NAJBI = IX2SQ(ISYMAJ) + NT1AM(ISYMAJ)*(LBI - 1) 3046 & + IT1AM(ISYMA,ISYMJ) + NVIR(ISYMA)*(J - 1) 3047 & + 1 3048 3049 CALL DCOPY(NVIR(ISYMA),T2AM(NAIBJ),1, 3050 & WORK(KAIBJ),1) 3051 CALL DCOPY(NVIR(ISYMA),T2AM(NAJBI),1, 3052 & WORK(KAJBI),1) 3053 CALL DSCAL(NVIR(ISYMA),TWO,T2AM(NAIBJ),1) 3054 CALL DSCAL(NVIR(ISYMA),TWO,T2AM(NAJBI),1) 3055 CALL DAXPY(NVIR(ISYMA),XMONE,WORK(KAJBI),1, 3056 & T2AM(NAIBJ),1) 3057 CALL DAXPY(NVIR(ISYMA),XMONE,WORK(KAIBJ),1, 3058 & T2AM(NAJBI),1) 3059 3060 ENDDO 3061 3062 ENDDO 3063 3064 ENDDO 3065 3066 ENDDO 3067 ENDDO 3068 ENDDO 3069 3070C Debug: compare results. 3071C ----------------------- 3072 3073 IF (LOCDBG) THEN 3074 CALL AROUND('Diff. TCME Amplitudes in '//SECNAM) 3075 CALL DAXPY(NX2SQ,XMONE,T2AM,1,WORK(KSCRT),1) 3076 DNRM = DDOT(NX2SQ,WORK(KSCRT),1,WORK(KSCRT),1) 3077 XNX2 = 1.0D0*NX2SQ 3078 RMS = DSQRT(DNRM/XNX2) 3079 WRITE(LUPRI,'(10X,A,I10,1X,1P,D14.6,/)') 3080 & 'Number of amplitudes and RMS error: ',NX2SQ,RMS 3081 IF (RMS .LE. TOL) GOTO 100 3082 WRITE(LUPRI,'(10X,A,1P,D14.6,A,/)') 3083 & 'Printing all differences larger than ',TOL,':' 3084 WRITE(LUPRI,'(1X,A)') 3085 & ' A ISYMA I ISYMI B ISYMB J ISYMJ Amplitude' 3086 WRITE(LUPRI,'(1X,A)') 3087 & '-------------------------------------------------------------' 3088 DO ISYMBJ = 1,NSYM 3089 ISYMAI = ISYMBJ 3090 DO ISYMJ = 1,NSYM 3091 ISYMB = MULD2H(ISYMJ,ISYMBJ) 3092 DO J = 1,NRHF(ISYMJ) 3093 DO LB = 1,LVIR(ISYMB) 3094 LBJ = IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) + LB 3095 B = IOFB1(ISYMB) + LB - 1 3096 DO ISYMI = 1,NSYM 3097 ISYMA = MULD2H(ISYMI,ISYMAI) 3098 DO I = 1,NRHF(ISYMI) 3099 DO A = 1,NVIR(ISYMA) 3100 NAI = IT1AM(ISYMA,ISYMI) 3101 & + NVIR(ISYMA)*(I - 1) + A 3102 NAIBJ = KSCRT + IX2SQ(ISYMBJ) 3103 & + NT1AM(ISYMAI)*(LBJ - 1) + NAI - 1 3104 IF (DABS(WORK(NAIBJ)) .GT. TOL) THEN 3105 WRITE(LUPRI,99) 3106 & A,ISYMA,I,ISYMI,B,ISYMB,J,ISYMJ, 3107 & WORK(NAIBJ) 3108 ENDIF 3109 ENDDO 3110 ENDDO 3111 ENDDO 3112 ENDDO 3113 ENDDO 3114 ENDDO 3115 ENDDO 3116 WRITE(LUPRI,'(1X,A)') 3117 & '-------------------------------------------------------------' 3118 100 CONTINUE 3119 ENDIF 3120 3121 RETURN 3122 99 FORMAT(1X,I5,3X,I1,3X,I5,3X,I1,3X,I5,3X,I1,3X,I5,3X,I1,2X,1P, 3123 & D14.6) 3124 END 3125C /* Deck cc2_choamg */ 3126 SUBROUTINE CC2_CHOAMG(T2AM,FOCKD,WORK,LWORK,LVIR,IOFB1,IX1AM, 3127 & NX1AM,IX2SQ,NX2SQ,T2NRM,CLNRM,MEMLN) 3128C 3129C Thomas Bondo Pedersen, August 2002. 3130C 3131C Purpose: Generate a batch of (minus) CC2 amplitudes t(ai,#bj). 3132C 3133C Input: 3134C 3135C FOCKD -- Orbital energies in CC ordering. 3136C Dummy if amplitudes are generated directly 3137C from T2 decomposition. (CHOT2 = .TRUE. in chocc2.h). 3138C 3139C WORK -- Work space, dimension: LWORK. 3140C 3141C LVIR -- Number of b's in each symmetry. 3142C I.e., a "local" version of NVIR(*). 3143C 3144C IOFB1 -- Offset to first b in each symmetry. 3145C I.e., IOFB1(ISYMB) returns the global index 3146C of the first b (in LVIR) of symmetry ISYMB. 3147C 3148C IX1AM -- Offset to symmetry blocks of bj. 3149C I.e., a "local" version of IT1AM(*,*). 3150C 3151C NX1AM -- Number of bj's in each symmetry. 3152C I.e., a "local" version of NT1AM(*). 3153C 3154C IX2SQ -- Offset to symmetry blocks ai,bj in T2AM. 3155C I.e., a "local" version of IT2SQ(*,*) 3156C except that only 1 symmetry argument is 3157C used, as T2AM is tot. symmetric. 3158C 3159C NX2SQ -- The total dimension of T2AM. 3160C I.e., a "local" version of NT2SQ(1). 3161C 3162C CLNRM -- .TRUE.: calculate T2NRM. 3163C 3164C Output: 3165C 3166C T2AM -- The amplitudes sorted according to IX2SQ. 3167C 3168C T2NRM -- Norm of the PACKED T2AM array. 3169C (Must be initialized on input, 3170C only calculated if CLNRM = .TRUE.) 3171C 3172C MEMLN -- Length of work space used (for statistics). 3173C 3174C Notes: 3175C 3176C There are currently 3 possible generations controlled 3177C by input in chocc2.h: 3178C 3179C (1) Generation from L(ai,J) vectors, i.e. transformed 3180C AO Cholesky vectors (CHOMO2 = .FALSE. and CHOT2 = .FALSE.). 3181C 3182C (2) Generation from separate decomposition of (ai|bj) 3183C integrals (CHOMO2 = .TRUE. and CHOT2 = .FALSE.). 3184C 3185C (3) Generation from separately decomposed CC2 amplitudes 3186C (CHOT2 = .TRUE.) 3187C 3188#include "implicit.h" 3189 DIMENSION T2AM(*), FOCKD(*), WORK(LWORK) 3190 INTEGER LVIR(8), IOFB1(8), IX1AM(8,8), NX1AM(8), IX2SQ(8) 3191 LOGICAL CLNRM 3192#include "maxorb.h" 3193#include "ccdeco.h" 3194#include "chocc2.h" 3195#include "ccorb.h" 3196#include "ccsdsym.h" 3197#include "priunit.h" 3198#include "ccsdinp.h" 3199 3200 INTEGER AI,BJ,AIBJ 3201 INTEGER NVECTOT(8) 3202 3203 LOGICAL LOCDBG 3204 PARAMETER (LOCDBG = .FALSE.) 3205 3206 CHARACTER*10 SECNAM 3207 PARAMETER (SECNAM = 'CC2_CHOAMG') 3208 3209 PARAMETER (IOPEN = -1, IKEEP = 1, INFO = 3) 3210 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0) 3211 3212C Initialize timings. 3213C ------------------- 3214 3215 TIMT = SECOND() 3216 TIMR = ZERO 3217 TIMS = ZERO 3218 TIMC = ZERO 3219 TIMD = ZERO 3220 TIMN = ZERO 3221 3222C Initialize MEMLN. 3223C ----------------- 3224 3225 MEMLN = 0 3226 3227C Debug: print input. 3228C ------------------- 3229 3230 IF (LOCDBG) THEN 3231 CALL HEADER('Input Variables Entering '//SECNAM,-1) 3232 WRITE(LUPRI,'(5X,A)') 3233 & 'LVIR - local NVIR array:' 3234 WRITE(LUPRI,'(8I10)') (LVIR(ISYM),ISYM=1,NSYM) 3235 WRITE(LUPRI,'(5X,A)') 3236 & 'IOFB1 - offsets to virtuals:' 3237 WRITE(LUPRI,'(8I10)') (IOFB1(ISYM),ISYM=1,NSYM) 3238 WRITE(LUPRI,'(5X,A)') 3239 & 'IX1AM - local IT1AM array:' 3240 DO JSYM = 1,NSYM 3241 WRITE(LUPRI,'(8I10)') (IX1AM(ISYM,JSYM),ISYM=1,NSYM) 3242 ENDDO 3243 WRITE(LUPRI,'(5X,A)') 3244 & 'NX1AM - local NT1AM array:' 3245 WRITE(LUPRI,'(8I10)') (NX1AM(ISYM),ISYM=1,NSYM) 3246 WRITE(LUPRI,'(5X,A)') 3247 & 'IX2SQ - local IT2SQ array:' 3248 WRITE(LUPRI,'(8I10)') (IX2SQ(ISYM),ISYM=1,NSYM) 3249 WRITE(LUPRI,'(5X,A,I10,/)') 3250 & 'NX2SQ - local dimension of T2AM: ',NX2SQ 3251 ENDIF 3252 3253C If nothing to do, return. 3254C ------------------------- 3255 3256 IF (NX2SQ .LE. 0) THEN 3257 IF (LOCDBG) THEN 3258 WRITE(LUPRI,'(5X,A,A,I10,/)') 3259 & SECNAM,' returns immediately: NX2SQ = ',NX2SQ 3260 ENDIF 3261 RETURN 3262 ENDIF 3263 3264C Initialize T2AM. 3265C ---------------- 3266 3267 CALL DZERO(T2AM,NX2SQ) 3268 3269C Set number of vectors and the type of vectors in CHO_MOP. 3270C --------------------------------------------------------- 3271 3272 IF (CHOT2) THEN 3273 ITYP = 6 3274 DO ISYM = 1,NSYM 3275 NVECTOT(ISYM) = NCHMO2(ISYM) 3276 ENDDO 3277 ELSE IF (CHOMO2) THEN 3278 ITYP = 5 3279 DO ISYM = 1,NSYM 3280 NVECTOT(ISYM) = NCHMO2(ISYM) 3281 ENDDO 3282 ELSE 3283 ITYP = 3 3284 DO ISYM = 1,NSYM 3285 NVECTOT(ISYM) = NUMCHO(ISYM) 3286 ENDDO 3287 ENDIF 3288 3289C Start Cholesky symmetry loop. 3290C ----------------------------- 3291 3292 DO ISYCHO = 1,NSYM 3293 3294 ISYMAI = ISYCHO 3295 ISYMBJ = ISYMAI 3296 3297C If nothing to do, skip this symmetry. 3298C ------------------------------------- 3299 3300 IF (NVECTOT(ISYCHO) .LE. 0) GOTO 1000 3301 IF (NT1AM(ISYMAI) .LE. 0) GOTO 1000 3302 IF (NX1AM(ISYMBJ) .LE. 0) GOTO 1000 3303 3304C Open file for this symmetry. 3305C ---------------------------- 3306 3307 DTIME = SECOND() 3308 CALL CHO_MOP(IOPEN,ITYP,ISYCHO,LUCHMO,1,1) 3309 DTIME = SECOND() - DTIME 3310 TIMR = TIMR + DTIME 3311 3312C Set up batch. 3313C ------------- 3314 3315 MINMEM = NT1AM(ISYMAI) + NX1AM(ISYMBJ) 3316 NVEC = MIN(LWORK/MINMEM,NVECTOT(ISYCHO)) 3317 3318 IF (NVEC .LE. 0) THEN 3319 WRITE(LUPRI,'(//,5X,A,A)') 3320 & 'Insufficient memory for batch in ',SECNAM 3321 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 3322 & 'Minimum memory required: ',MINMEM, 3323 & 'Total memory available : ',LWORK 3324 CALL QUIT('Insufficient memory in '//SECNAM) 3325 ENDIF 3326 3327 NBATCH = (NVECTOT(ISYCHO) - 1)/NVEC + 1 3328 3329C Batch loop. 3330C ----------- 3331 3332 DO IBATCH = 1,NBATCH 3333 3334 NUMV = NVEC 3335 IF (IBATCH .EQ. NBATCH) THEN 3336 NUMV = NVECTOT(ISYCHO) - NVEC*(NBATCH - 1) 3337 ENDIF 3338 3339 JVEC1 = NVEC*(IBATCH - 1) + 1 3340 3341C Complete allocation. 3342C -------------------- 3343 3344 KCHAI = 1 3345 KCHBJ = KCHAI + NT1AM(ISYMAI)*NUMV 3346 KLAST = KCHBJ + NX1AM(ISYMBJ)*NUMV - 1 3347 3348 MEMLN = MAX(MEMLN,KLAST) 3349 3350C Read vectors. 3351C ------------- 3352 3353 DTIME = SECOND() 3354 CALL CHO_MOREAD(WORK(KCHAI),NT1AM(ISYMAI),NUMV,JVEC1,LUCHMO) 3355 DTIME = SECOND() - DTIME 3356 TIMR = TIMR + DTIME 3357 3358C Copy out the #bj block. 3359C ----------------------- 3360 3361 DTIME = SECOND() 3362 DO IVEC = 1,NUMV 3363 DO ISYMJ = 1,NSYM 3364 3365 ISYMB = MULD2H(ISYMJ,ISYMBJ) 3366 IF (LVIR(ISYMB) .LE. 0) GOTO 999 3367 3368 DO J = 1,NRHF(ISYMJ) 3369 3370 KOFF1 = KCHAI + NT1AM(ISYCHO)*(IVEC - 1) 3371 & + IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) 3372 & + IOFB1(ISYMB) - 1 3373 KOFF2 = KCHBJ + NX1AM(ISYCHO)*(IVEC - 1) 3374 & + IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) 3375 3376 CALL DCOPY(LVIR(ISYMB),WORK(KOFF1),1, 3377 & WORK(KOFF2),1) 3378 3379 ENDDO 3380 3381 999 CONTINUE 3382 3383 ENDDO 3384 ENDDO 3385 DTIME = SECOND() - DTIME 3386 TIMS = TIMS + DTIME 3387 3388C Calculate the product. 3389C ---------------------- 3390 3391 DTIME = SECOND() 3392 3393 NTOAI = MAX(NT1AM(ISYMAI),1) 3394 NTOBJ = MAX(NX1AM(ISYMBJ),1) 3395 3396 KOFFT = IX2SQ(ISYCHO) + 1 3397 3398 CALL DGEMM('N','T',NT1AM(ISYMAI),NX1AM(ISYMBJ),NUMV, 3399 & ONE,WORK(KCHAI),NTOAI,WORK(KCHBJ),NTOBJ, 3400 & ONE,T2AM(KOFFT),NTOAI) 3401 3402 DTIME = SECOND() - DTIME 3403 TIMC = TIMC + DTIME 3404 3405 ENDDO 3406 3407C Close file for this symmetry. 3408C ----------------------------- 3409 3410 DTIME = SECOND() 3411 CALL CHO_MOP(IKEEP,ITYP,ISYCHO,LUCHMO,1,1) 3412 DTIME = SECOND() - DTIME 3413 TIMR = TIMR + DTIME 3414 3415 1000 CONTINUE 3416 3417 ENDDO 3418 3419C If necessary, divide by orbital energy denominators. 3420C ---------------------------------------------------- 3421 3422 IF (.NOT. CHOT2) THEN 3423 3424 TIMD = SECOND() 3425 3426 DO ISYMBJ = 1,NSYM 3427 3428 IF (NX1AM(ISYMBJ) .LE. 0) GOTO 1002 3429 3430 ISYMAI = ISYMBJ 3431 3432 DO ISYMJ = 1,NSYM 3433 3434 ISYMB = MULD2H(ISYMJ,ISYMBJ) 3435 3436 IF (LVIR(ISYMB) .LE. 0) GOTO 1001 3437 3438 DO J = 1,NRHF(ISYMJ) 3439 3440 KOFFJ = IRHF(ISYMJ) + J 3441 3442 DO LB = 1,LVIR(ISYMB) 3443 3444 B = IOFB1(ISYMB) + LB - 1 3445 KOFFB = IVIR(ISYMB) + B 3446 3447 LBJ = IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) + LB 3448 3449 DENBJ = FOCKD(KOFFB) - FOCKD(KOFFJ) 3450 3451 DO ISYMI = 1,NSYM 3452 3453 ISYMA = MULD2H(ISYMI,ISYMAI) 3454 3455 DO I = 1,NRHF(ISYMI) 3456 3457 KOFFI = IRHF(ISYMI) + I 3458 3459 DNBJI = DENBJ - FOCKD(KOFFI) 3460 3461 DO A = 1,NVIR(ISYMA) 3462 3463 KOFFA = IVIR(ISYMA) + A 3464 3465 AI = IT1AM(ISYMA,ISYMI) 3466 & + NVIR(ISYMA)*(I - 1) + A 3467 3468 AIBJ = IX2SQ(ISYMBJ) 3469 & + NT1AM(ISYMAI)*(LBJ - 1) + AI 3470 3471 DENOM = DNBJI + FOCKD(KOFFA) 3472 3473 T2AM(AIBJ) = T2AM(AIBJ)/DENOM 3474 3475 ENDDO 3476 3477 ENDDO 3478 3479 ENDDO 3480 3481 ENDDO 3482 3483 ENDDO 3484 3485 1001 CONTINUE 3486 3487 ENDDO 3488 3489 1002 CONTINUE 3490 3491 ENDDO 3492 3493 TIMD = SECOND() - TIMD 3494 3495 ENDIF 3496 3497C If requested, calculate T2NRM. 3498C ------------------------------ 3499 3500 IF (CLNRM) THEN 3501 3502 TIMN = SECOND() 3503 3504 DO ISYMBJ = 1,NSYM 3505 ISYMAI = ISYMBJ 3506 DO ISYMJ = 1,NSYM 3507 ISYMB = MULD2H(ISYMJ,ISYMBJ) 3508 IF (LVIR(ISYMB) .GT. 0) THEN 3509 DO J = 1,NRHF(ISYMJ) 3510 DO LB = 1,LVIR(ISYMB) 3511 LBJ = IX1AM(ISYMB,ISYMJ) + LVIR(ISYMB)*(J - 1) 3512 & + LB 3513 B = IOFB1(ISYMB) + LB - 1 3514 BJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) 3515 & + B 3516 DO AI = 1,BJ 3517 AIBJ = IX2SQ(ISYMBJ) 3518 & + NT1AM(ISYMAI)*(LBJ - 1) + AI 3519 T2NRM = T2NRM + T2AM(AIBJ)*T2AM(AIBJ) 3520 ENDDO 3521 ENDDO 3522 ENDDO 3523 ENDIF 3524 ENDDO 3525 ENDDO 3526 3527 TIMN = SECOND() - TIMN 3528 3529 ENDIF 3530 3531C Print. 3532C ------ 3533 3534 IF (IPRINT .GE. INFO) THEN 3535 TIMT = SECOND() - TIMT 3536 CALL HEADER('Information from '//SECNAM,-1) 3537 IF (CHOT2) THEN 3538 WRITE(LUPRI,'(5X,A,/)') 3539 & 'CC2 doubles amplitudes separately decomposed' 3540 ELSE IF (CHOMO2) THEN 3541 WRITE(LUPRI,'(5X,A,/)') 3542 & 'CC2 doubles amplitudes from separately decomposed (ai|bj)' 3543 ELSE 3544 WRITE(LUPRI,'(5X,A,/)') 3545 & 'CC2 doubles amplitudes calculated from (ai|bj)' 3546 ENDIF 3547 IF (CLNRM) THEN 3548 WRITE(LUPRI,'(5X,A,1P,D30.15)') 3549 & 'Accumulated amplitude norm: ',DSQRT(T2NRM) 3550 ENDIF 3551 WRITE(LUPRI,'(5X,A,F10.2,A)') 3552 & 'Time used for I/O of Cholesky vectors: ',TIMR,' seconds' 3553 WRITE(LUPRI,'(5X,A,F10.2,A)') 3554 & 'Time used for sort of MO vectors : ',TIMS,' seconds' 3555 WRITE(LUPRI,'(5X,A,F10.2,A)') 3556 & 'Time used for multiplying vectors : ',TIMC,' seconds' 3557 IF (.NOT. CHOT2) THEN 3558 WRITE(LUPRI,'(5X,A,F10.2,A)') 3559 & 'Time used for orbital denominators : ',TIMD,' seconds' 3560 ENDIF 3561 IF (CLNRM) THEN 3562 WRITE(LUPRI,'(5X,A,F10.2,A)') 3563 & 'Time used for calculating T2 norm : ',TIMN,' seconds' 3564 ENDIF 3565 WRITE(LUPRI,'(5X,A)') 3566 & '---------------------------------------------------------' 3567 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 3568 & 'Time used in total : ',TIMT,' seconds' 3569 ENDIF 3570 3571 RETURN 3572 END 3573C /* Deck chocc2_jg */ 3574 SUBROUTINE CHOCC2_JG(T1AM,OMEGA1,XLAMDP,XLAMDH,TRACE,WORK,LWORK) 3575C 3576C Thomas Bondo Pedersen, August 2002. 3577C 3578C Purpose: Calculate the CC2 JG-term. 3579C 3580C Notice: TRACE is scaled by 2 on exit !!! 3581C 3582#include "implicit.h" 3583 DIMENSION T1AM(*), OMEGA1(*), XLAMDP(*), XLAMDH(*), TRACE(*) 3584 DIMENSION WORK(LWORK) 3585#include "priunit.h" 3586#include "maxorb.h" 3587#include "ccdeco.h" 3588#include "ccorb.h" 3589#include "ccsdsym.h" 3590#include "ccsdinp.h" 3591#include "chocc2.h" 3592 3593 INTEGER IOFFC(8), IOFFZ(8), IOFFY(8) 3594 3595 INTEGER LROW, ICOL, ADDR 3596 3597 CHARACTER*9 SECNAM 3598 PARAMETER (SECNAM = 'CHOCC2_JG') 3599 3600 PARAMETER (IOPEN = -1, IKEEP = 1, ITYKI = 4, IOPTR = 2) 3601 PARAMETER (INFO = 5) 3602 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 3603 3604C Initialize timings. 3605C ------------------- 3606 3607 TIMTOT = SECOND() 3608 TIMYIO = ZERO 3609 TIMCIO = ZERO 3610 TIMMIO = ZERO 3611 TIMYAD = ZERO 3612 TIMYBK = ZERO 3613 TIMYSO = ZERO 3614 TIMCSO = ZERO 3615 TIMOMA = ZERO 3616 TIMOMM = ZERO 3617 3618C Scale trace by 2. 3619C ----------------- 3620 3621 CALL DSCAL(NUMCHO(1),TWO,TRACE,1) 3622 3623C Read reduce index array. 3624C ------------------------ 3625 3626 DTIME = SECOND() 3627 3628 KIND1 = 1 3629 CALL CC_GETIND1(WORK(KIND1),LWORK,LIND1) 3630 KEND0 = KIND1 + LIND1 3631 LWRK0 = LWORK - KEND0 + 1 3632 3633 IF (LWRK0 .LT. 0) THEN 3634 WRITE(LUPRI,'(//,5X,A,A,A)') 3635 & 'Insufficient memory in ',SECNAM,' - allocation: index' 3636 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 3637 & 'Need (more than): ',KEND0-1, 3638 & 'Available : ',LWORK 3639 CALL QUIT('Insufficient memory in '//SECNAM) 3640 ENDIF 3641 3642 DTIME = SECOND() - DTIME 3643 TIMCIO = TIMCIO + DTIME 3644 3645C Allocation. 3646C ----------- 3647 3648 KOMAO = KEND0 3649 KEND1 = KOMAO + NT1AO(1) 3650 LWRK1 = LWORK - KEND1 + 1 3651 3652 IF (LWRK1 .LT. 0) THEN 3653 WRITE(LUPRI,'(//,5X,A,A,A)') 3654 & 'Insufficient memory in ',SECNAM,' - allocation: Omega' 3655 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 3656 & 'Need (more than): ',KEND1-1, 3657 & 'Available : ',LWORK 3658 CALL QUIT('Insufficient memory in '//SECNAM) 3659 ENDIF 3660 3661C Initialize AO result vector. 3662C ---------------------------- 3663 3664 CALL DZERO(WORK(KOMAO),NT1AO(1)) 3665 3666C Start Cholesky symmetry loop. 3667C ----------------------------- 3668 3669 DO ISYCHO = 1,NSYM 3670 3671 IF (NUMCHO(ISYCHO) .LE. 0) GOTO 1000 3672 IF (N2BST(ISYCHO) .LE. 0) GOTO 1000 3673 IF (NT1AM(ISYCHO) .LE. 0) GOTO 1000 3674 IF (NT1AO(ISYCHO) .LE. 0) GOTO 1000 3675 3676C Open file with Y-intermediates. 3677C ------------------------------- 3678 3679 DTIME = SECOND() 3680 CALL WOPEN2(LCC2YM,FCC2YM(ISYCHO),64,0) 3681 DTIME = SECOND() - DTIME 3682 TIMYIO = TIMYIO + DTIME 3683 3684C Open file with L(ki,J). 3685C ----------------------- 3686 3687 DTIME = SECOND() 3688 CALL CHO_MOP(IOPEN,ITYKI,ISYCHO,LUCHKI,1,1) 3689 DTIME = SECOND() - DTIME 3690 TIMMIO = TIMMIO + DTIME 3691 3692C Allocation. 3693C ----------- 3694 3695 LREAD = MEMRD(1,ISYCHO,IOPTR) 3696 3697 KSCRC = KEND1 3698 KREAD = KSCRC + N2BST(ISYCHO) 3699 KEND2 = KREAD + LREAD 3700 LWRK2 = LWORK - KEND2 + 1 3701 3702 IF (LWRK2 .LE. 0) THEN 3703 WRITE(LUPRI,'(//,5X,A,A,A)') 3704 & 'Insufficient memory in ',SECNAM,' - allocation: Cholesky' 3705 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 3706 & 'Need (more than): ',KEND2-1, 3707 & 'Available : ',LWORK 3708 CALL QUIT('Insufficient memory in '//SECNAM) 3709 ENDIF 3710 3711C Set up batch. 3712C ------------- 3713 3714 LENC = MAX(N2BST(ISYCHO),NT1AM(ISYCHO),NMATIJ(ISYCHO)) 3715 LENZ = MAX(NT1AO(ISYCHO),NT1AM(ISYCHO)) 3716 MINMEM = LENC + LENZ 3717 NVEC = MIN(LWRK2/MINMEM,NUMCHO(ISYCHO)) 3718 3719 IF (NVEC .LE. 0) THEN 3720 WRITE(LUPRI,'(//,5X,A,A)') 3721 & 'Insufficient memory for batch in ',SECNAM 3722 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10,/)') 3723 & 'Minimum memory needed: ',MINMEM, 3724 & 'Available for batch : ',LWRK2, 3725 & 'Available in total : ',LWORK 3726 CALL QUIT('Insufficient memory in '//SECNAM) 3727 ENDIF 3728 3729 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 3730 3731C Start batch loop. 3732C ----------------- 3733 3734 DO IBATCH = 1,NBATCH 3735 3736 NUMV = NVEC 3737 IF (IBATCH .EQ. NBATCH) THEN 3738 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 3739 ENDIF 3740 3741C Set first vector in this batch. 3742C ------------------------------- 3743 3744 JVEC1 = NVEC*(IBATCH - 1) + 1 3745 3746C Set up index arrays IOFFC, IOFFY, and IOFFZ. 3747C -------------------------------------------- 3748 3749 ICOUN1 = 0 3750 ICOUN2 = 0 3751 ICOUN3 = 0 3752 3753 DO ISYMB = 1,NSYM 3754 3755 ISYMA = MULD2H(ISYMB,ISYCHO) 3756 ISYMI = ISYMA 3757 ISYMC = ISYMB 3758 3759 IOFFC(ISYMB) = ICOUN1 3760 IOFFZ(ISYMB) = ICOUN2 3761 IOFFY(ISYMC) = ICOUN3 3762 3763 LENBJ = NBAS(ISYMB)*NUMV 3764 LENCJ = NVIR(ISYMC)*NUMV 3765 ICOUN1 = ICOUN1 + NBAS(ISYMA)*LENBJ 3766 ICOUN2 = ICOUN2 + LENBJ*NRHF(ISYMI) 3767 ICOUN3 = ICOUN3 + LENCJ*NRHF(ISYMI) 3768 3769 ENDDO 3770 3771C Complete allocation. 3772C -------------------- 3773 3774 KCHOL = KEND2 3775 KZMAT = KCHOL + LENC*NUMV 3776 3777C Read Y(ci,#J) in place of Z. 3778C ---------------------------- 3779 3780 DTIME = SECOND() 3781 3782 ICOL = JVEC1 - 1 3783 LROW = NT1AM(ISYCHO) 3784 ADDR = LROW*ICOL + 1 3785 LEN = NT1AM(ISYCHO)*NUMV 3786 CALL GETWA2(LCC2YM,FCC2YM(ISYCHO),WORK(KZMAT),ADDR,LEN) 3787 3788 DTIME = SECOND() - DTIME 3789 TIMYIO = TIMYIO + DTIME 3790 3791C Read L(ki,#J) in place of L(al,be,#J). 3792C -------------------------------------- 3793 3794 DTIME = SECOND() 3795 CALL CHO_MOREAD(WORK(KCHOL),NMATIJ(ISYCHO),NUMV,JVEC1, 3796 & LUCHKI) 3797 DTIME = SECOND() - DTIME 3798 TIMMIO = TIMMIO + DTIME 3799 3800C Subtract sum(k) T1AM(ck) * L(ki,#J) from Y(ci,#J). 3801C -------------------------------------------------- 3802 3803 DTIME = SECOND() 3804 3805 DO IVEC = 1,NUMV 3806 DO ISYMI = 1,NSYM 3807 3808 ISYMK = MULD2H(ISYMI,ISYCHO) 3809 ISYMC = ISYMK 3810 3811 KOFF1 = IT1AM(ISYMC,ISYMK) + 1 3812 KOFF2 = KCHOL + NMATIJ(ISYCHO)*(IVEC - 1) 3813 & + IMATIJ(ISYMK,ISYMI) 3814 KOFF3 = KZMAT + NT1AM(ISYCHO)*(IVEC - 1) 3815 & + IT1AM(ISYMC,ISYMI) 3816 3817 NC = NVIR(ISYMC) 3818 NTOTC = MAX(NC,1) 3819 NI = NRHF(ISYMI) 3820 NK = NRHF(ISYMK) 3821 NTOTK = MAX(NK,1) 3822 3823 CALL DGEMM('N','N',NC,NI,NK, 3824 & XMONE,T1AM(KOFF1),NTOTC,WORK(KOFF2),NTOTK, 3825 & ONE,WORK(KOFF3),NTOTC) 3826 3827 ENDDO 3828 ENDDO 3829 3830 DTIME = SECOND() - DTIME 3831 TIMYAD = TIMYAD + DTIME 3832 3833C Reorder to obtain Y(c#J,i) in place of L(al,be,#J). 3834C --------------------------------------------------- 3835 3836 DTIME = SECOND() 3837 3838 DO IVEC = 1,NUMV 3839 DO ISYMI = 1,NSYM 3840 3841 ISYMC = MULD2H(ISYMI,ISYCHO) 3842 3843 LENCJ = NVIR(ISYMC)*NUMV 3844 3845 DO I = 1,NRHF(ISYMI) 3846 3847 KOFF1 = KZMAT + NT1AM(ISYCHO)*(IVEC - 1) 3848 & + IT1AM(ISYMC,ISYMI) + NVIR(ISYMC)*(I - 1) 3849 KOFF2 = KCHOL + IOFFY(ISYMC) + LENCJ*(I - 1) 3850 & + NVIR(ISYMC)*(IVEC - 1) 3851 3852 CALL DCOPY(NVIR(ISYMC),WORK(KOFF1),1,WORK(KOFF2),1) 3853 3854 ENDDO 3855 3856 ENDDO 3857 ENDDO 3858 3859 DTIME = SECOND() - DTIME 3860 TIMYSO = TIMYSO + DTIME 3861 3862C Calculate Z(be#J,i) by backtransforming virtual index of Y. 3863C ----------------------------------------------------------- 3864 3865 DTIME = SECOND() 3866 DO ISYMI = 1,NSYM 3867 3868 ISYMC = MULD2H(ISYMI,ISYCHO) 3869 ISYMB = ISYMC 3870 3871 KOFF1 = IGLMVI(ISYMB,ISYMC) + 1 3872 KOFF2 = KCHOL + IOFFY(ISYMC) 3873 KOFF3 = KZMAT + IOFFZ(ISYMB) 3874 3875 NTOTB = MAX(NBAS(ISYMB),1) 3876 NTOTC = MAX(NVIR(ISYMC),1) 3877 3878 CALL DGEMM('N','N', 3879 & NBAS(ISYMB),NUMV*NRHF(ISYMI),NVIR(ISYMC), 3880 & ONE,XLAMDH(KOFF1),NTOTB,WORK(KOFF2),NTOTC, 3881 & ZERO,WORK(KOFF3),NTOTB) 3882 3883 ENDDO 3884 DTIME = SECOND() - DTIME 3885 TIMYBK = TIMYBK + DTIME 3886 3887C Include Coulomb part from J-term. 3888C N.B.: TRACE has already been scaled by 2 above. 3889C ----------------------------------------------- 3890 3891 IF (ISYCHO .EQ. 1) THEN 3892 3893 DTIME = SECOND() 3894 DO ISYMB = 1,NSYM 3895 3896 IF (NBAS(ISYMB) .GT. 0) THEN 3897 3898 ISYMI = ISYMB 3899 3900 LENBJ = NBAS(ISYMB)*NUMV 3901 3902 DO I = 1,NRHF(ISYMI) 3903 3904 KOFF2 = IGLMRH(ISYMB,ISYMI) 3905 & + NBAS(ISYMB)*(I - 1) + 1 3906 3907 DO IVEC = 1,NUMV 3908 3909 KOFF1 = JVEC1 + IVEC - 1 3910 KOFF3 = KZMAT + IOFFZ(ISYMB) 3911 & + LENBJ*(I - 1) 3912 & + NBAS(ISYMB)*(IVEC - 1) 3913 3914 CALL DAXPY(NBAS(ISYMB),TRACE(KOFF1), 3915 & XLAMDH(KOFF2),1,WORK(KOFF3),1) 3916 3917 ENDDO 3918 3919 ENDDO 3920 3921 ENDIF 3922 3923 ENDDO 3924 DTIME = SECOND() - DTIME 3925 TIMYAD = TIMYAD + DTIME 3926 3927 ENDIF 3928 3929C Set up L(al,be#J). 3930C ------------------ 3931 3932 DO IVEC = 1,NUMV 3933 3934 DTIME = SECOND() 3935 JVEC = JVEC1 + IVEC - 1 3936 CALL CHO_READN(WORK(KSCRC),JVEC,1,WORK(KIND1),IDUM2, 3937 & ISYCHO,IOPTR,WORK(KREAD),LREAD) 3938 DTIME = SECOND() - DTIME 3939 TIMCIO = TIMCIO + DTIME 3940 3941 DTIME = SECOND() 3942 DO ISYMB = 1,NSYM 3943 3944 ISYMA = MULD2H(ISYMB,ISYCHO) 3945 3946 LENAB = NBAS(ISYMA)*NBAS(ISYMB) 3947 3948 KOFF1 = KSCRC + IAODIS(ISYMA,ISYMB) 3949 KOFF2 = KCHOL + IOFFC(ISYMB) + LENAB*(IVEC - 1) 3950 3951 CALL DCOPY(LENAB,WORK(KOFF1),1,WORK(KOFF2),1) 3952 3953 ENDDO 3954 DTIME = SECOND() - DTIME 3955 TIMCSO = TIMCSO + DTIME 3956 3957 ENDDO 3958 3959C Calculate contribution in AO basis. 3960C ----------------------------------- 3961 3962 DTIME = SECOND() 3963 DO ISYMB = 1,NSYM 3964 3965 ISYMI = MULD2H(ISYMB,ISYCHO) 3966 ISYMA = ISYMI 3967 3968 KOFF1 = KCHOL + IOFFC(ISYMB) 3969 KOFF2 = KZMAT + IOFFZ(ISYMB) 3970 KOFF3 = KOMAO + IT1AO(ISYMA,ISYMI) 3971 3972 LENBJ = NBAS(ISYMB)*NUMV 3973 3974 NTOTA = MAX(NBAS(ISYMA),1) 3975 NTOBJ = MAX(LENBJ,1) 3976 3977 CALL DGEMM('N','N',NBAS(ISYMA),NRHF(ISYMI),LENBJ, 3978 & ONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOBJ, 3979 & ONE,WORK(KOFF3),NTOTA) 3980 3981 ENDDO 3982 DTIME = SECOND() - DTIME 3983 TIMOMA = TIMOMA + DTIME 3984 3985 ENDDO 3986 3987C Close file with L(ki,J). 3988C ------------------------ 3989 3990 DTIME = SECOND() 3991 CALL CHO_MOP(IKEEP,ITYKI,ISYCHO,LUCHKI,1,1) 3992 DTIME = SECOND() - DTIME 3993 TIMMIO = TIMMIO + DTIME 3994 3995C Close file with Y-intermediates. 3996C -------------------------------- 3997 3998 DTIME = SECOND() 3999 CALL WCLOSE2(LCC2YM,FCC2YM(ISYCHO),'KEEP') 4000 DTIME = SECOND() - DTIME 4001 TIMYIO = TIMYIO + DTIME 4002 4003 1000 CONTINUE 4004 4005 ENDDO 4006 4007C Transform AO index to virtual. 4008C ------------------------------ 4009 4010 DTIME = SECOND() 4011 DO ISYMI = 1,NSYM 4012 4013 ISYMA = ISYMI 4014 ISYMG = ISYMA 4015 4016 NTOTG = MAX(NBAS(ISYMG),1) 4017 NTOTA = MAX(NVIR(ISYMA),1) 4018 4019 KOFF1 = IGLMVI(ISYMG,ISYMA) + 1 4020 KOFF2 = KOMAO + IT1AO(ISYMG,ISYMI) 4021 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 4022 4023 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMG), 4024 & ONE,XLAMDP(KOFF1),NTOTG,WORK(KOFF2),NTOTG, 4025 & ONE,OMEGA1(KOFF3),NTOTA) 4026 4027 ENDDO 4028 DTIME = SECOND() - DTIME 4029 TIMOMM = TIMOMM + DTIME 4030 4031C Print. 4032C ------ 4033 4034 IF (IPRINT .GE. INFO) THEN 4035 TIMTOT = SECOND() - TIMTOT 4036 CALL HEADER('Timing of '//SECNAM,-1) 4037 WRITE(LUPRI,'(5X,A,F10.2,A)') 4038 & 'Time used for I/O of Y-intermediates : ',TIMYIO,' seconds' 4039 WRITE(LUPRI,'(5X,A,F10.2,A)') 4040 & 'Time used for J-contributions to Z : ',TIMYAD,' seconds' 4041 WRITE(LUPRI,'(5X,A,F10.2,A)') 4042 & 'Time used for reorderering Y : ',TIMYSO,' seconds' 4043 WRITE(LUPRI,'(5X,A,F10.2,A)') 4044 & 'Time used for backtransformation of Y : ',TIMYBK,' seconds' 4045 WRITE(LUPRI,'(5X,A,F10.2,A)') 4046 & 'Time used for I/O of MO Cholesky vecs.: ',TIMMIO,' seconds' 4047 WRITE(LUPRI,'(5X,A,F10.2,A)') 4048 & 'Time used for I/O of AO Cholesky vecs.: ',TIMCIO,' seconds' 4049 WRITE(LUPRI,'(5X,A,F10.2,A)') 4050 & 'Time used for reorder. AO Chol. vecs.: ',TIMCSO,' seconds' 4051 WRITE(LUPRI,'(5X,A,F10.2,A)') 4052 & 'Time used for calculating AO Omega1 : ',TIMOMA,' seconds' 4053 WRITE(LUPRI,'(5X,A,F10.2,A)') 4054 & 'Time used for transforming AO Omega1 : ',TIMOMM,' seconds' 4055 WRITE(LUPRI,'(5X,A)') 4056 & '----------------------------------------------------------' 4057 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 4058 & 'Total time : ',TIMTOT,' seconds' 4059 ENDIF 4060 4061 RETURN 4062 END 4063C /* Deck chocc2_h */ 4064 SUBROUTINE CHOCC2_H(OMEGA1,WORK,LWORK) 4065C 4066C Thomas Bondo Pedersen, August 2002. 4067C 4068C Purpose: Calculate the CC2 H-term. 4069C 4070#include "implicit.h" 4071 DIMENSION OMEGA1(*), WORK(LWORK) 4072#include "maxorb.h" 4073#include "ccdeco.h" 4074#include "ccorb.h" 4075#include "ccsdsym.h" 4076#include "ccsdinp.h" 4077#include "chocc2.h" 4078#include "priunit.h" 4079 4080 INTEGER IOFFX(8), IOFFC(8) 4081 4082 INTEGER LROW, ICOL, ADDR 4083 4084 CHARACTER*8 SECNAM 4085 PARAMETER (SECNAM = 'CHOCC2_H') 4086 4087 PARAMETER (IOPEN = -1, IKEEP = 1, INFO = 5) 4088 PARAMETER (ITYKI = 4) 4089 PARAMETER (XMONE = -1.0D0, ZERO = 0.0D0, ONE = 1.0D0) 4090 4091C Initialize timings. 4092C ------------------- 4093 4094 TIMTOT = SECOND() 4095 TIMCIO = ZERO 4096 TIMYIO = ZERO 4097 TIMXIM = ZERO 4098 TIMCAL = ZERO 4099 4100C Start Cholesky symmetry loop. 4101C ----------------------------- 4102 4103 DO ISYCHO = 1,NSYM 4104 4105 IF (NUMCHO(ISYCHO) .LE. 0) GOTO 1000 4106 IF (NT1AM(ISYCHO) .LE. 0) GOTO 1000 4107 IF (NMATIJ(ISYCHO) .LE. 0) GOTO 1000 4108 4109C Open files for MO Cholesky vectors. 4110C ----------------------------------- 4111 4112 DTIME = SECOND() 4113 CALL CHO_MOP(IOPEN,ITYKI,ISYCHO,LUCHKI,1,1) 4114 DTIME = SECOND() - DTIME 4115 TIMCIO = TIMCIO + DTIME 4116 4117C Open file for Y-intermediates. 4118C ------------------------------ 4119 4120 DTIME = SECOND() 4121 CALL WOPEN2(LCC2YM,FCC2YM(ISYCHO),64,0) 4122 DTIME = SECOND() - DTIME 4123 TIMYIO = TIMYIO + DTIME 4124 4125C Allocation: scratch space for read. 4126C ----------------------------------- 4127 4128 KSCR1 = 1 4129 KEND1 = KSCR1 + MAX(NT1AM(ISYCHO),NMATIJ(ISYCHO)) 4130 LWRK1 = LWORK - KEND1 + 1 4131 4132 IF (LWRK1 .LE. 0) THEN 4133 WRITE(LUPRI,'(//,5X,A,A,A)') 4134 & 'Insufficient memory in ',SECNAM,' - allocation: scratch' 4135 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/)') 4136 & 'Need (more than): ',KEND1-1, 4137 & 'Available : ',LWORK 4138 CALL QUIT('Insufficient memory in '//SECNAM) 4139 ENDIF 4140 4141C Set up batch. 4142C ------------- 4143 4144 MINMEM = NT1AM(ISYCHO) + NMATIJ(ISYCHO) 4145 NVEC = MIN(LWRK1/MINMEM,NUMCHO(ISYCHO)) 4146 4147 IF (NVEC .LE. 0) THEN 4148 WRITE(LUPRI,'(//,5X,A,A)') 4149 & 'Insufficient memory for batch in ',SECNAM 4150 WRITE(LUPRI,'(5X,A,I10,/,5X,A,I10,/,5X,A,I10,/)') 4151 & 'Minimum memory required for batch: ',MINMEM, 4152 & 'Available for batch : ',LWRK1, 4153 & 'Available in total : ',LWORK 4154 CALL QUIT('Insufficient memory for batch in '//SECNAM) 4155 ENDIF 4156 4157 NBATCH = (NUMCHO(ISYCHO) - 1)/NVEC + 1 4158 4159C Start batch loop. 4160C ----------------- 4161 4162 DO IBATCH = 1,NBATCH 4163 4164 NUMV = NVEC 4165 IF (IBATCH .EQ. NBATCH) THEN 4166 NUMV = NUMCHO(ISYCHO) - NVEC*(NBATCH - 1) 4167 ENDIF 4168 4169 JVEC1 = NVEC*(IBATCH - 1) + 1 4170 4171C Complete allocation. 4172C -------------------- 4173 4174 KYMAT = KEND1 4175 KCHKI = KYMAT + NT1AM(ISYCHO)*NUMV 4176 4177C Set up index arrays IOFFX and IOFFC. 4178C ------------------------------------ 4179 4180 ICOUN1 = 0 4181 ICOUN2 = 0 4182 4183 DO ISYMK = 1,NSYM 4184 4185 ISYMA = MULD2H(ISYMK,ISYCHO) 4186 ISYMI = ISYMA 4187 4188 IOFFX(ISYMK) = ICOUN1 4189 IOFFC(ISYMK) = ICOUN2 4190 4191 LENKJ = NRHF(ISYMK)*NUMV 4192 ICOUN1 = ICOUN1 + NVIR(ISYMA)*LENKJ 4193 ICOUN2 = ICOUN2 + LENKJ*NRHF(ISYMI) 4194 4195 ENDDO 4196 4197 DO IVEC = 1,NUMV 4198 4199C Set current vector. 4200C ------------------- 4201 4202 JVEC = JVEC1 + IVEC - 1 4203 4204C Read L(ki,J). 4205C ------------- 4206 4207 DTIME = SECOND() 4208 CALL CHO_MOREAD(WORK(KSCR1),NMATIJ(ISYCHO),1,JVEC,LUCHKI) 4209 DTIME = SECOND() - DTIME 4210 TIMCIO = TIMCIO + DTIME 4211 4212C Reorder: L(k#J,i). 4213C ------------------ 4214 4215 DTIME = SECOND() 4216 DO ISYMI = 1,NSYM 4217 4218 ISYMK = MULD2H(ISYMI,ISYCHO) 4219 4220 IF (NRHF(ISYMK) .GT. 0) THEN 4221 LENKJ = NRHF(ISYMK)*NUMV 4222 DO I = 1,NRHF(ISYMI) 4223 4224 KOFF1 = KSCR1 + IMATIJ(ISYMK,ISYMI) 4225 & + NRHF(ISYMK)*(I - 1) 4226 KOFF2 = KCHKI + IOFFC(ISYMK) 4227 & + LENKJ*(I - 1) + NRHF(ISYMK)*(IVEC - 1) 4228 4229 CALL DCOPY(NRHF(ISYMK),WORK(KOFF1),1, 4230 & WORK(KOFF2),1) 4231 4232 ENDDO 4233 ENDIF 4234 4235 ENDDO 4236 DTIME = SECOND() - DTIME 4237 TIMXIM = TIMXIM + DTIME 4238 4239C Read Y(ak,J). 4240C ------------- 4241 4242 DTIME = SECOND() 4243 ICOL = JVEC - 1 4244 LROW = NT1AM(ISYCHO) 4245 ADDR = LROW*ICOL + 1 4246 CALL GETWA2(LCC2YM,FCC2YM(ISYCHO),WORK(KSCR1),ADDR, 4247 & NT1AM(ISYCHO)) 4248 DTIME = SECOND() - DTIME 4249 TIMYIO = TIMYIO + DTIME 4250 4251C Reorder: Y(a,k#J). 4252C ------------------ 4253 4254 DTIME = SECOND() 4255 DO ISYMK = 1,NSYM 4256 4257 ISYMA = MULD2H(ISYMK,ISYCHO) 4258 4259 LENAK = NVIR(ISYMA)*NRHF(ISYMK) 4260 4261 KOFF1 = KSCR1 + IT1AM(ISYMA,ISYMK) 4262 KOFF2 = KYMAT + IOFFX(ISYMK) 4263 & + LENAK*(IVEC - 1) 4264 4265 CALL DCOPY(LENAK,WORK(KOFF1),1,WORK(KOFF2),1) 4266 4267 ENDDO 4268 DTIME = SECOND() - DTIME 4269 TIMXIM = TIMXIM + DTIME 4270 4271 ENDDO 4272 4273C Calculate contribution. 4274C ----------------------- 4275 4276 DTIME = SECOND() 4277 DO ISYMK = 1,NSYM 4278 4279 ISYMA = MULD2H(ISYMK,ISYCHO) 4280 ISYMI = ISYMA 4281 4282 KOFF1 = KYMAT + IOFFX(ISYMK) 4283 KOFF2 = KCHKI + IOFFC(ISYMK) 4284 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 4285 4286 LENKJ = NRHF(ISYMK)*NUMV 4287 NTOTA = MAX(NVIR(ISYMA),1) 4288 NTOKJ = MAX(LENKJ,1) 4289 4290 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),LENKJ, 4291 & XMONE,WORK(KOFF1),NTOTA,WORK(KOFF2),NTOKJ, 4292 & ONE,OMEGA1(KOFF3),NTOTA) 4293 4294 ENDDO 4295 DTIME = SECOND() - DTIME 4296 TIMCAL = TIMCAL + DTIME 4297 4298 ENDDO 4299 4300C Close file for Y-intermediates. 4301C ------------------------------- 4302 4303 DTIME = SECOND() 4304 CALL WCLOSE2(LCC2YM,FCC2YM(ISYCHO),'KEEP') 4305 DTIME = SECOND() - DTIME 4306 TIMYIO = TIMYIO + DTIME 4307 4308C Close file for MO Cholesky vectors. 4309C ----------------------------------- 4310 4311 DTIME = SECOND() 4312 CALL CHO_MOP(IKEEP,ITYKI,ISYCHO,LUCHKI,1,1) 4313 DTIME = SECOND() - DTIME 4314 TIMCIO = TIMCIO + DTIME 4315 4316 1000 CONTINUE 4317 4318 ENDDO 4319 4320C Print. 4321C ------ 4322 4323 IF (IPRINT .GE. INFO) THEN 4324 TIMTOT = SECOND() - TIMTOT 4325 CALL HEADER('Timing of '//SECNAM,-1) 4326 WRITE(LUPRI,'(5X,A,F10.2,A)') 4327 & 'Time used for I/O of Y-intermediates : ',TIMYIO,' seconds' 4328 WRITE(LUPRI,'(5X,A,F10.2,A)') 4329 & 'Time used for I/O of MO Chol. vectors : ',TIMCIO,' seconds' 4330 WRITE(LUPRI,'(5X,A,F10.2,A)') 4331 & 'Time used for reordering factor arrays : ',TIMXIM,' seconds' 4332 WRITE(LUPRI,'(5X,A,F10.2,A)') 4333 & 'Time used for calculating contribution : ',TIMCAL,' seconds' 4334 WRITE(LUPRI,'(5X,A)') 4335 & '-----------------------------------------------------------' 4336 WRITE(LUPRI,'(5X,A,F10.2,A,/)') 4337 & 'Total time : ',TIMTOT,' seconds' 4338 ENDIF 4339 4340 RETURN 4341 END 4342C /* Deck cc2_iniome */ 4343 SUBROUTINE CC2_INIOME(FOCKD,T1AM,OMEGA1,ISYM) 4344C 4345C Thomas Bondo Pedersen, August 2002. 4346C 4347C Purpose: Initialize the CC2 vector function: 4348C 4349C OMEGA1(ai) = [e(a) - e(i)] * T1AM(ai) 4350C 4351C 4352#include "implicit.h" 4353 DIMENSION FOCKD(*), T1AM(*), OMEGA1(*) 4354#include "ccorb.h" 4355#include "ccsdsym.h" 4356 4357 INTEGER AI 4358 4359 DO ISYMI = 1,NSYM 4360 4361 ISYMA = MULD2H(ISYMI,ISYM) 4362 4363 DO I = 1,NRHF(ISYMI) 4364 4365 KOFFI = IRHF(ISYMI) + I 4366 4367 DO A = 1,NVIR(ISYMA) 4368 4369 KOFFA = IVIR(ISYMA) + A 4370 AI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I - 1) + A 4371 4372 OMEGA1(AI) = (FOCKD(KOFFA) - FOCKD(KOFFI))*T1AM(AI) 4373 4374 ENDDO 4375 4376 ENDDO 4377 4378 ENDDO 4379 4380 RETURN 4381 END 4382C 4383C /* Deck cho_fckh */ 4384 SUBROUTINE CHO_FCKH(T1AM,WORK,LWORK) 4385C 4386C Write Fock matrix in CC_FCKH needed in CC property calculation 4387C Use Cholesky vectors 4388C 4389C asm September 2005 4390C 4391#include "implicit.h" 4392#include "maxorb.h" 4393#include "priunit.h" 4394#include "dummy.h" 4395#include "inftap.h" 4396#include "ccorb.h" 4397#include "ccdeco.h" 4398#include "chocc2.h" 4399#include "ccsdsym.h" 4400#include "ccsdinp.h" 4401#include "ccfield.h" 4402C 4403 LOGICAL TINFO, IOPTPV 4404C 4405 DIMENSION T1AM(*) 4406 DIMENSION WORK(LWORK) 4407C 4408C 4409C--------------- 4410C Initialize 4411C--------------- 4412C 4413 TINFO = IPRINT .GE. 5 4414 ISYDEN = 1 4415C 4416C----------------- 4417C Allocation.1 4418C----------------- 4419C 4420 KDENSI = 1 4421 KLAMDP = KDENSI + N2BST(ISYDEN) 4422 KLAMDH = KLAMDP + NLAMDT 4423 KEND0 = KLAMDH + NLAMDT 4424 LWRK0 = LWORK - KEND0 + 1 4425C 4426 IF (LWRK0 .LT. 0) CALL QUIT('Not enough space in CHO_FCKH.0') 4427C 4428C------------------------------ 4429C Construct density matrix 4430C------------------------------ 4431C 4432c IOPT = 1 4433c CALL CC_RDRSP('R0 ',1,1,IOPT,MODEL,WORK(KT1AM),DUMMY) 4434C 4435 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),T1AM,WORK(KEND0),LWRK0) 4436C 4437 IC = 1 4438 CALL CC_AODENS(WORK(KLAMDP),WORK(KLAMDH),WORK(KDENSI),ISYDEN, 4439 & IC,WORK(KEND0),LWRK0) 4440C 4441 IF (IPRINT .GT. 50) THEN 4442 CALL AROUND( 'Density matrix' ) 4443 CALL CC_PRFCKAO(WORK(KDENSI),ISYDEN) 4444 END IF 4445C 4446C----------------- 4447C Allocation.2 4448C----------------- 4449C 4450 KFOCK = KLAMDP 4451 KEND1 = KFOCK + N2BST(ISYDEN) 4452 LWRK1 = LWORK - KEND1 + 1 4453C 4454 IF (LWRK1 .LT. 0) CALL QUIT('Not enough space in CHO_FCKH.1') 4455C 4456C-------------------------- 4457C Construct Fock matrix 4458C-------------------------- 4459C 4460 CALL DZERO(WORK(KFOCK),N2BST(ISYDEN)) 4461C 4462 CALL CCRHS_ONEAO(WORK(KFOCK),WORK(KEND1),LWRK1) 4463C 4464 DO I = 1, NFIELD 4465 FIELD = EFIELD(I) 4466 CALL CC_ONEP(WORK(KFOCK),WORK(KEND1),LWRK1,FIELD,1,LFIELD(I)) 4467 END DO 4468C 4469 IOPTPV = .TRUE. 4470 CALL SCF_CHOFCK3(WORK(KDENSI),WORK(KFOCK),WORK(KEND1),LWRK1, 4471 & ISYDEN,TINFO,IOPTPV) 4472C 4473C------------------ 4474C Write do disk 4475C------------------ 4476C 4477 LUFCK = -1 4478 CALL GPOPEN(LUFCK,'CC_FCKH','UNKNOWN',' ','UNFORMATTED',IDUMMY, 4479 & .FALSE.) 4480 REWIND(LUFCK) 4481 WRITE(LUFCK)(WORK(KFOCK+I-1), I = 1,N2BST(ISYDEN)) 4482 CALL GPCLOSE(LUFCK,'KEEP' ) 4483C 4484 IF (IPRINT .GT. 50) THEN 4485 CALL AROUND( 'Fock AO matrix written to disk' ) 4486 CALL CC_PRFCKAO(WORK(KFOCK),ISYDEN) 4487 END IF 4488C 4489 RETURN 4490 END 4491 4492