1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck cclr_dia */ 20 SUBROUTINE CCLR_DIA(DIAA,ISYMTR,TRIPLET,WORK,LWORK) 21C 22C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 23C 24C Written by Ove Christiansen 4-9-1994 25C Triplet extensions by Christof Haettig & Kasper Hald, May 1999 26C R12 extension by Christof Haettig, Jun 2003 27C 28C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 29C 30 IMPLICIT NONE 31#include "dummy.h" 32#include "priunit.h" 33#include "inftap.h" 34#include "ccorb.h" 35#include "ccsdsym.h" 36#include "ccsdinp.h" 37Cholesky 38#include "maxorb.h" 39#include "ccdeco.h" 40Cholesky 41 42 LOGICAL TRIPLET, LOCDBG 43 PARAMETER (LOCDBG = .FALSE.) 44 INTEGER LWORK 45 46#if defined (SYS_CRAY) 47 REAL WORK(*), DIAA(*) 48 REAL ONE, HALF, TMPKH 49#else 50 DOUBLE PRECISION WORK(*), DIAA(*) 51 DOUBLE PRECISION ONE, HALF, TMPKH 52#endif 53 PARAMETER (ONE = 1.0D0, HALF = 0.5D00, TMPKH= 1.0D08) 54C 55 INTEGER KOFF1, KOFF2, KOFF2P, KOFF2M, KOFF12, KFOCKD, KEND1, 56 & LWRK1, ISYMA, ISYMI, MI, MA, KAI, ISYMBJ, ISYMAI, 57 & ISYMB, ISYMJ, MJ, MB, NAI, NBJ, NAIBJ, NAIAI, ISYMTR, 58 & INDEX 59C 60 INDEX(I,J) = MAX(I,J)*(MAX(I,J) -3 )/2 + I + J 61C 62 CALL QENTER('CCLR_DIA') 63 64C ---------------------------------------------------- 65C set start addresses for the individual blocks parts: 66C (singles, doubles, R12 doubles) 67C ---------------------------------------------------- 68 KOFF1 = 1 69 KOFF2 = KOFF1 + NT1AM(ISYMTR) 70 IF (TRIPLET) THEN 71 KOFF2P = KOFF2 72 KOFF2M = KOFF2P + NT2AM(ISYMTR) 73 KOFF12 = KOFF2M + NT2AMA(ISYMTR) 74 ELSE 75 KOFF12 = KOFF2 + NT2AM(ISYMTR) 76 END IF 77C 78C----------------------------- 79C Allocation of workspace. 80C----------------------------- 81C 82 KFOCKD = 1 83 KEND1 = KFOCKD + NORBTS 84 LWRK1 = LWORK - KEND1 85C 86 IF ( LWRK1 .LT. 0 ) THEN 87 CALL QUIT('Insufficient space in CCLR_DIA ') 88 ENDIF 89C 90C------------------------------------- 91C Read canonical orbital energies. 92C------------------------------------- 93 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY, 94 & .FALSE.) 95 REWIND LUSIFC 96 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 97 READ (LUSIFC) 98 READ (LUSIFC) (WORK(KFOCKD+I-1), I=1,NORBTS) 99 CALL GPCLOSE(LUSIFC,'KEEP') 100C 101 IF (IPRINT.GT.100) THEN 102 CALL AROUND('FOCK DIAGONAL ') 103 CALL OUTPUT(WORK(KFOCKD),1,NORBTS,1,1,NORBTS,1,1,LUPRI) 104 ENDIF 105C 106 IF (FROIMP .OR. FROEXP) 107 * CALL CCSD_DELFRO(WORK(KFOCKD),WORK(KEND1),LWRK1) 108 109C------------------------------------------------ 110C FOCK based approximative A(1,1) diagonal. 111C------------------------------------------------ 112 DO ISYMA = 1,NSYM 113 ISYMI = MULD2H(ISYMA,ISYMTR) 114 DO I = 1,NRHF(ISYMI) 115 MI = IORB(ISYMI) + I 116 DO A = 1,NVIR(ISYMA) 117 MA = IORB(ISYMA) + NRHF(ISYMA) + A 118 KAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 119 DIAA(KAI) = WORK(MA) - WORK(MI) 120 END DO 121 END DO 122 END DO 123C 124 IF (IPRINT.GT.20 .OR. DEBUG .OR. LOCDBG) THEN 125 CALL AROUND('FOCK BASED DIAGONAL single exc. ') 126 CALL OUTPUT(DIAA,1,NT1AM(ISYMTR),1,1,NT1AM(ISYMTR),1,1, 127 & LUPRI) 128 ENDIF 129C 130 IF (CCS .OR. CIS) GOTO 9999 131Cholesky 132 IF (CHOINT .AND. CC2) THEN 133 CALL QEXIT('CCLR_DIA') 134 RETURN 135 END IF 136Cholesky 137C 138C------------------------------------------------ 139C FOCK based approximative A(2,2) diagonal. 140C------------------------------------------------ 141 DO 200 ISYMBJ = 1,NSYM 142C 143 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 144C 145 DO 210 ISYMJ = 1,NSYM 146C 147 ISYMB = MULD2H(ISYMJ,ISYMBJ) 148C 149 DO 220 ISYMI = 1,NSYM 150C 151 ISYMA = MULD2H(ISYMI,ISYMAI) 152C 153 DO 230 J = 1,NRHF(ISYMJ) 154C 155 MJ = IORB(ISYMJ) + J 156C 157 DO 240 B = 1,NVIR(ISYMB) 158C 159 NBJ = IT1AM(ISYMB,ISYMJ) 160 * + NVIR(ISYMB)*(J - 1) + B 161C 162 MB = IORB(ISYMB) + NRHF(ISYMB) + B 163C 164 DO 250 I = 1,NRHF(ISYMI) 165C 166 MI = IORB(ISYMI) + I 167C 168 DO 260 A = 1,NVIR(ISYMA) 169C 170 NAI = IT1AM(ISYMA,ISYMI) 171 * + NVIR(ISYMA)*(I - 1) + A 172C 173 MA = IORB(ISYMA) + NRHF(ISYMA) + A 174C 175 IF (ISYMAI.EQ.ISYMBJ) THEN 176 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 177 * + INDEX(NAI,NBJ) 178 ELSE IF (ISYMAI.LT.ISYMBJ) THEN 179 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 180 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 181 ELSE IF (ISYMAI.GT.ISYMBJ) THEN 182 NAIBJ = IT2AM(ISYMBJ,ISYMAI) 183 * + NT1AM(ISYMBJ)*(NAI-1) + NBJ 184 ENDIF 185C 186 DIAA(NAIBJ+NT1AM(ISYMTR)) = 187 * WORK(MA)+WORK(MB)-WORK(MI)-WORK(MJ) 188C 189 260 CONTINUE 190 250 CONTINUE 191 240 CONTINUE 192 230 CONTINUE 193 220 CONTINUE 194 210 CONTINUE 195 200 CONTINUE 196C 197 IF (TRIPLET) THEN 198C -------------------------------------------------------- 199C 2- block as for singlet vector, but with the ai,ai 200C diagonal set to a huge value to ensure that the 201C corresponding elements in the trial vectors are zero. 202C -------------------------------------------------------- 203 CALL DCOPY(NT2AM(ISYMTR),DIAA(KOFF2P),1,DIAA(KOFF2M),1) 204 IF (ISYMTR.EQ.1) THEN 205 DO ISYMAI = 1,NSYM 206 DO ISYMI = 1,NSYM 207 ISYMA = MULD2H(ISYMI,ISYMAI) 208 DO I = 1,NRHF(ISYMI) 209 MI = IORB(ISYMI) + I 210 DO A = 1,NVIR(ISYMA) 211 MA = IORB(ISYMA) + NRHF(ISYMA) + A 212 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1)+A 213 NAIAI = IT2AM(ISYMAI,ISYMAI) + INDEX(NAI,NAI) 214 DIAA(KOFF2M-1+NAIAI) = 1.0D08 215 END DO 216 END DO 217 END DO 218 END DO 219 END IF 220 221C ----------------------------------------------------------- 222C 2+ block as for singlet vector, but times a factor of 0.5 223C and with the elements a<=b or i<=j set to huge value 224C to ensure that the corresponding elements in the trial 225C vectors are zero 226C ----------------------------------------------------------- 227 CALL DCOPY(NT2AM(ISYMTR),TMPKH,0,DIAA(KOFF2P),1) 228 DO ISYMBJ = 1,NSYM 229 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 230C 231 IF (ISYMAI .LE. ISYMBJ) THEN 232! 233 DO ISYMI = 1,NSYM 234 ISYMA = MULD2H(ISYMI,ISYMAI) 235 DO ISYMJ = 1,NSYM 236 ISYMB = MULD2H(ISYMJ,ISYMBJ) 237C 238 DO I = 1,NRHF(ISYMI) 239 DO J = 1,NRHF(ISYMJ) 240C 241 MI = IORB(ISYMI) + I 242 MJ = IORB(ISYMJ) + J 243 IF (MI .NE. MJ) THEN 244C 245 DO A = 1,NVIR(ISYMA) 246 DO B = 1,NVIR(ISYMB) 247C 248 MA = IORB(ISYMA) + NRHF(ISYMA) + A 249 MB = IORB(ISYMB) + NRHF(ISYMB) + B 250 IF (MA .NE. MB) THEN 251C 252 NBJ = IT1AM(ISYMB,ISYMJ)+NVIR(ISYMB)*(J-1)+B 253 NAI = IT1AM(ISYMA,ISYMI)+NVIR(ISYMA)*(I-1)+A 254C 255 IF (ISYMAI.EQ.ISYMBJ) THEN 256C 257 IF ((NAI .LT. NBJ) .AND. 258 * (MA .LT. MB)) THEN 259C 260 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 261 * + INDEX(NAI,NBJ) 262 DIAA(NT1AM(ISYMTR)+NAIBJ) = 263 * DIAA(KOFF2M-1+NAIBJ) 264! 265 ENDIF 266 ELSE IF (ISYMAI.LT.ISYMBJ) THEN 267! 268 IF (((MA .LT. MB) .AND. 269 * (MI .LT. MJ)) .OR. 270 * ((MA .GT. MB) .AND. 271 * (MI .GT. MJ))) THEN 272! 273 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 274 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 275 DIAA(NT1AM(ISYMTR)+NAIBJ) = 276 * DIAA(KOFF2M-1+NAIBJ) 277! 278 ENDIF 279 ENDIF 280 281 ENDIF 282! 283 END DO 284 END DO 285! 286 ENDIF 287C 288 END DO 289 END DO 290C 291 END DO 292 END DO 293C 294! 295 ENDIF 296! 297 END DO 298C 299 ENDIF 300C 301 IF (CCR12) THEN 302C CALL CCLR_DIAR12(DIAA(KOFF12),WORK,LWORK) 303C CALL QUIT('Missing subroutine in CCLR_DIA...') 304c Write(LUPRI,*)'Warning Missing subroutine in CCLR_DIA...' 305c Write(LUPRI,*)'Set R12 Jacobian diagonal to 1.0D01...' 306 CALL DCOPY(NTR12AM(ISYMTR),1.0D01,0,DIAA(KOFF12),1) 307 END IF 308C 309 IF (IPRINT.GT. 20 .OR. DEBUG .OR. LOCDBG) THEN 310 CALL AROUND('APPROXIMATIVE DIAGONAL 22 BLOCK ') 311 IF (TRIPLET) THEN 312 CALL OUTPUT(DIAA(1+NT1AM(ISYMTR)),1,2*NT2AM(ISYMTR), 313 * 1,1,2*NT2AM(ISYMTR),1,1,LUPRI) 314 ELSE 315 CALL OUTPUT(DIAA(1+NT1AM(ISYMTR)),1,1*NT2AM(ISYMTR), 316 * 1,1,1*NT2AM(ISYMTR),1,1,LUPRI) 317 END IF 318 IF (CCR12) THEN 319 WRITE(LUPRI,*) 'R12 doubles part:' 320 CALL OUTPUT(DIAA(KOFF12),1,NTR12AM(ISYMTR), 321 * 1,1,NTR12AM(ISYMTR),1,1,LUPRI) 322 END IF 323 ENDIF 324C 325 9999 CONTINUE 326C 327 CALL QEXIT('CCLR_DIA') 328 RETURN 329 END 330C /* Deck cclr_diascl */ 331 SUBROUTINE CCLR_DIASCL(OMEGA2,X,ISYMTR) 332C 333C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 334C 335C Written by Ove Christiansen 6-5-1994 336C symmetry 16-2-1995. 337C 338C Purpose: Scale diagonal of omega2 with x 339C 340C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 341C 342#include "implicit.h" 343 DIMENSION OMEGA2(*) 344#include "ccorb.h" 345#include "ccsdsym.h" 346C 347 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 348C 349 IF (ISYMTR.EQ.1) THEN 350C 351 DO 50 ISYMBJ = 1, NSYM 352C 353 ISYMAI = MULD2H(ISYMBJ,ISYMTR) 354C 355 DO 100 NBJ = 1,NT1AM(ISYMBJ) 356C 357 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NBJ,NBJ) 358 OMEGA2(NAIBJ) = OMEGA2(NAIBJ)*X 359C 360 100 CONTINUE 361C 362 50 CONTINUE 363 364 ENDIF 365C 366 RETURN 367 END 368C /* Deck cclr_1c1f */ 369 SUBROUTINE CCLR_1C1F(RHO1,FOCKC,ISYMTR) 370C 371C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 372C 373C Ove Christiansen 19-Jan-1994 374C 375C Purpose: Add Sum-k-Laikk-bar to rho1(ai) 376C 377C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 378C 379#include "implicit.h" 380 PARAMETER (ONE = 1.0D0, ZERO = 0.0D0) 381 DIMENSION RHO1(*),FOCKC(*) 382#include "ccorb.h" 383#include "ccsdsym.h" 384C 385C---------------------- 386C Add contribution. 387C---------------------- 388C 389 ISYMAI = ISYMTR 390C 391 KOFF1 = 1 392C 393 DO 110 ISYMI = 1,NSYM 394C 395 ISYMA = MULD2H(ISYMI,ISYMAI) 396C 397 DO 120 I = 1, NRHF(ISYMI) 398C 399 K1 = KOFF1 + NORB(ISYMA)*(I-1) + NRHF(ISYMA) 400 K2 = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + 1 401 CALL DAXPY(NVIR(ISYMA),ONE,FOCKC(K1),1,RHO1(K2),1) 402C 403 120 CONTINUE 404C 405 KOFF1 = KOFF1 + NORB(ISYMA)*NRHF(ISYMI) 406C 407 110 CONTINUE 408C 409 END 410C /* Deck cclr_e1c1 */ 411 SUBROUTINE CCLR_E1C1(RHO1,C1AM,EMAT1,EMAT2,WORK,LWORK, 412 * ISYMV,ISYMIM,TRANS) 413C 414C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 415C Written by Ove Christiansen 31-Jan-1995 416C 417C input: Ematrices transformed in ccrhs_E 418C Purpose: Calculate E-terms in 1C1 part of linear transformation. 419C 420C Calculates rho(ai) = Sum(c)C1AM(c,i)*E1(a,c) (If trans then E1(c,a) 421C - Sum(k)C1AM(a,k)*E2(k,i) (If trans then E2(i,k) 422C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 423C 424#include "implicit.h" 425 PARAMETER (ONE = 1.0D0, XMONE=-1.0D0,TWO = 2.0D0) 426 DIMENSION EMAT1(*), EMAT2(*),C1AM(*) 427 DIMENSION WORK(LWORK),RHO1(*) 428 CHARACTER*1 TRANS 429#include "ccorb.h" 430#include "ccsdsym.h" 431#include "ccsdinp.h" 432C 433C---------------------------------------- 434C Calculate 1C1 contribution from E1. 435C---------------------------------------- 436C 437 ISYMAI = MULD2H(ISYMIM,ISYMV) 438C 439 DO 100 ISYMA = 1, NSYM 440C 441 ISYMI = MULD2H(ISYMAI,ISYMA) 442 ISYMC = MULD2H(ISYMV ,ISYMI) 443C 444 NVIRA = MAX(NVIR(ISYMA),1) 445 NVIRC = MAX(NVIR(ISYMC),1) 446 KOFF2 = IT1AM(ISYMC,ISYMI) + 1 447C 448 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 449C 450 IF (TRANS.EQ.'N') THEN 451C 452 KOFF1 = IMATAB(ISYMA,ISYMC) + 1 453C 454 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC), 455 * ONE,EMAT1(KOFF1),NVIRA,C1AM(KOFF2),NVIRC,ONE, 456 * RHO1(KOFF3),NVIRA) 457C 458 ELSE IF (TRANS.EQ.'T') THEN 459C 460 KOFF1 = IMATAB(ISYMC,ISYMA) + 1 461C 462 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NVIR(ISYMC), 463 * ONE,EMAT1(KOFF1),NVIRC,C1AM(KOFF2),NVIRC,ONE, 464 * RHO1(KOFF3),NVIRA) 465 ENDIF 466 100 CONTINUE 467C 468C---------------------------------------- 469C Calculate 1C1 contribution from E2. 470C---------------------------------------- 471C 472 DO 200 ISYMA = 1, NSYM 473C 474 ISYMI = MULD2H(ISYMAI,ISYMA) 475 ISYMK = MULD2H(ISYMV ,ISYMA) 476C 477 KOFF1 = IT1AM(ISYMA,ISYMK) + 1 478 KOFF3 = IT1AM(ISYMA,ISYMI) + 1 479C 480 NVIRA = MAX(NVIR(ISYMA),1) 481 NRHFK = MAX(NRHF(ISYMK),1) 482 NRHFI = MAX(NRHF(ISYMI),1) 483C 484 IF (TRANS.EQ.'N') THEN 485C 486 KOFF2 = IMATIJ(ISYMK,ISYMI) + 1 487C 488 CALL DGEMM('N','N',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 489 * XMONE,C1AM(KOFF1),NVIRA,EMAT2(KOFF2),NRHFK,ONE, 490 * RHO1(KOFF3),NVIRA) 491C 492 ELSE IF (TRANS.EQ.'T') THEN 493C 494 KOFF2 = IMATIJ(ISYMI,ISYMK) + 1 495C 496 CALL DGEMM('N','T',NVIR(ISYMA),NRHF(ISYMI),NRHF(ISYMK), 497 * XMONE,C1AM(KOFF1),NVIRA,EMAT2(KOFF2),NRHFI,ONE, 498 * RHO1(KOFF3),NVIRA) 499C 500 ENDIF 501C 502 200 CONTINUE 503C 504 RETURN 505C 506 END 507c*DECK CCLR_JAC 508 SUBROUTINE CCLR_JAC(ISIDE,WORK,LWORK,APROXR12,TRIPLET) 509C 510C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 511C 512C Written by Ove Christiansen May/November 1994, 513C to calculate the Jacobian by finite difference and by 514C construction with transformation from the right on the CC jacobian. 515C 516C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 517C 518#include "implicit.h" 519#include "maxorb.h" 520#include "iratdef.h" 521#include "ccorb.h" 522#include "aovec.h" 523C 524 LOGICAL RSPIM2 525 DIMENSION WORK(LWORK) 526 CHARACTER*24 CBLOCK 527 CHARACTER*(*) APROXR12 528 LOGICAL TRIPLET 529C 530#include "ccsdinp.h" 531C Note cclr.h includes ccexcinf.h 532#include "cclr.h" 533#include "ccsdsym.h" 534#include "ccsdio.h" 535#include "leinf.h" 536#include "priunit.h" 537C 538 IF (IPRINT .GT. 10) THEN 539 CALL AROUND(' START OF CCLR_JAC ') 540 ENDIF 541C 542 IF ((NSYM.GT.1).AND.(JACTST.OR.FDJAC)) THEN 543 WRITE(LUPRI,*) 544 * 'Finite difference Joacobian does only work without symmetry' 545 CALL QUIT( 546 * 'Finite difference Joacobian does only work without symmetry') 547 ENDIF 548C 549C-------------------- 550C Initialization. 551C-------------------- 552C 553 DO ISYMTR = 1, NSYM 554 IPRLE = IPRINT 555 IF (JACTST.OR.JACEXP.OR.FDJAC) THEN 556 NC1VEC = NT1AM(ISYMTR) 557 NC2VEC = NT2AM(ISYMTR) 558 IF (TRIPLET) NC2VEC = 2*NC2VEC 559c IF (NT1AMX.GT.3) THEN 560c NC1VEC = 3 561c ENDIF 562c IF (NT2AMX.GT.5) THEN 563c NC2VEC = 5 564c ENDIF 565 ENDIF 566 IF (CCR12) THEN 567 NTEMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) + NTR12AM(ISYMTR) 568 NTEMP2 = NTEMP*(NC1VEC + NC2VEC + NTR12AM(ISYMTR)) 569 ELSE 570 NTEMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 571 IF (TRIPLET) NTEMP = NTEMP + NT2AM(ISYMTR) 572 NTEMP2 = NTEMP*(NC1VEC + NC2VEC) 573 END IF 574C 575 KEND1 = 1 576 LWRK1 = LWORK 577C 578C------------------------------------------------------------ 579C Calculate Jacobian by Transformation with unit vectors. 580C First elements in workspace contains the jacobian. 581C------------------------------------------------------------ 582C 583 IF (IPRINT.GT.10) WRITE(LUPRI,*) 'CCLR_JAC Workspace:',LWORK 584C 585 IF (JACTST.OR.JACEXP) THEN 586 IF(ISYMTR.EQ.1) THEN 587 CALL AROUND( 'Calculation of analytical CC Jacobian') 588 END IF 589 KJACAN = KEND1 590 KEND1 = KJACAN + NTEMP2 591 LWRK1 = LWORK - KEND1 592 IF ( LWRK1 .LE. 0 ) THEN 593 WRITE(LUPRI,*) 'out of memory:' 594 WRITE(LUPRI,*) 'need:',KEND1 595 WRITE(LUPRI,*) 'have:',LWORK 596 CALL QUIT('Too little workspace in CCLR_JAC') 597 END IF 598 CALL CCLR_JACAN(NC1VEC,NC2VEC,ISIDE,WORK(KJACAN),LWORK, 599 & APROXR12,TRIPLET) 600 ENDIF 601 602 END DO ! ISYMTR = 1, NSYM 603 604 ISYMTR = 1 605C 606C------------------------------------------------------- 607C Calculate Jacobian by finite difference. 608C both for jactst and jacexp. 609C First elements in workspace contains the jacobian. 610C------------------------------------------------------- 611C 612 IF (JACTST.OR.FDJAC) THEN 613 KFDJAC = KEND1 614 KEND2 = KFDJAC + NTEMP2 615 LWRK2 = LWORK - KEND2 616 IF ( LWRK2 .LE. 0 ) THEN 617 WRITE(LUPRI,*) 'out of memory:' 618 WRITE(LUPRI,*) 'need:',KEND1 619 WRITE(LUPRI,*) 'have:',LWORK 620 CALL QUIT('Too little workspace in CCLR_JAC-2') 621 END IF 622 CALL AROUND('Calculation of finite difference CC Jacobian') 623 RSPIM2 = RSPIM 624 RSPIM = .FALSE. 625 CALL CCLR_FDJAC(NC1VEC,NC2VEC,WORK(KFDJAC),LWRK1,APROXR12) 626 RSPIM = RSPIM2 627 END IF 628C 629C-------------------------------------------------------- 630C calculate difference between the two jacobiants. 631C-------------------------------------------------------- 632C 633 IF (JACTST) THEN 634 CALL DAXPY(NTEMP2,-1.0D00,WORK(KFDJAC),1,WORK(KJACAN),1) 635 IDXMX = IDAMAX(NTEMP2,WORK(KJACAN),1) 636 IDXCOL = 1 + (IDXMX-1)/NTEMP 637 IDXROW = IDXMX - ((IDXCOL-1)*NTEMP) 638 IF (IDXCOL.LE.NC1VEC) THEN 639 CBLOCK = 'singles block:' 640 ELSE IF (IDXCOL.LE.(NC1VEC+NC2VEC)) THEN 641 CBLOCK = 'doubles block:' 642 IDXCOL = IDXCOL - NC1VEC 643 ELSE IF (IDXCOL.LE.NCCVAR) THEN 644 CBLOCK = 'R12 doubles block:' 645 IDXCOL = IDXCOL - NC1VEC - NC2VEC 646 ELSE 647 CALL QUIT('Error in CCLR_JAC') 648 END IF 649 IF ( IPRINT .GE. 5) THEN 650 CALL AROUND( 'DIFFERENCE. CC JACOBIAN - singles PART ' ) 651 CALL OUTPUT(WORK(KJACAN),1,NTEMP,1,NC1VEC,NTEMP,NC1VEC,1, 652 & LUPRI) 653 CALL AROUND( 'DIFFERENCE. CC JACOBIAN - doubles PART ' ) 654 CALL OUTPUT(WORK(KJACAN+NTEMP*NC1VEC),1,NTEMP,1,NC2VEC, 655 * NTEMP,NC2VEC,1,LUPRI) 656 IF (CCR12) THEN 657 CALL AROUND('DIFFERENCE. CC JACOBIAN - R12 doubles PART') 658 CALL OUTPUT(WORK(KJACAN+NTEMP*(NC1VEC+NC2VEC)), 659 * 1,NTEMP,1,NTR12AM(1),NTEMP,NTR12AM(1),1,LUPRI) 660 END IF 661 ENDIF 662 DIFNOR = DDOT(NTEMP2,WORK(KJACAN),1,WORK(KJACAN),1) 663 WRITE(LUPRI,'(/A,E20.10/)') ' The norm of the difference'// 664 * ' between fd and analytical jacobian is ',SQRT(DIFNOR) 665 WRITE(LUPRI,'(/2A/,3I5,E20.10/)') 666 * ' Largest element in difference is in ', 667 * CBLOCK, IDXMX, IDXCOL, IDXROW, WORK(KJACAN-1+IDXMX) 668 WRITE(LUPRI,'(A/)') ' END OF JACTST TEST' 669 ENDIF 670C 671 IF (IPRINT .GT. 10) THEN 672 CALL AROUND(' END OF CCLR_JAC ') 673 ENDIF 674C 675 RETURN 676 END 677c*DECK CCLR_JACAN 678 SUBROUTINE CCLR_JACAN(NC1VEC,NC2VEC,ISIDE,WORK,LWORK, 679 & APROXR12,TRIPLET) 680C 681C------------------------------------------------------------------------- 682C Test routine for calculating the CC Jacobian by Transformation of 683C unit vectors. 684C calculates the first nc1vec rows of the C1 blocks 685C and nc2vec of the rows of the c2 blocks. 686C Written by Ove Christiansen 26-5-1994 687C------------------------------------------------------------------------- 688C 689#include "implicit.h" 690#include "maxorb.h" 691#include "iratdef.h" 692#include "ccorb.h" 693#include "aovec.h" 694#include "ccsdinp.h" 695#include "cclr.h" 696#include "ccsdsym.h" 697#include "ccsdio.h" 698#include "leinf.h" 699#include "priunit.h" 700C 701 DIMENSION WORK(LWORK) 702 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-08) 703 CHARACTER*(*) APROXR12 704 LOGICAL TRIPLET 705C 706 INTEGER :: NT2 707C 708 IF (IPRINT .GT. 10) THEN 709 CALL AROUND(' START OF CCLR_JACAN') 710 ENDIF 711C 712C---------------------------- 713C Work space allocations. 714C---------------------------- 715C 716 !ISYMTR = 1 717 IF (CCR12) THEN 718 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) + NTR12AM(ISYMTR) 719 NTAMP2 = NTAMP*(NC1VEC + NC2VEC + NTR12AM(ISYMTR)) 720 ELSE 721 NT2 = NT2AM(ISYMTR) 722 IF (TRIPLET) NT2 = 2*NT2 723 NTAMP = NT1AM(ISYMTR) + NT2 724 NTAMP2 = NTAMP*(NC1VEC + NC2VEC ) 725 END IF 726 727 KJACOBI = 1 728 KEND1 = 1 + NTAMP2 729 IF (CCR12) THEN 730 KMETRIC = KEND1 731 KEND1 = KMETRIC + NTAMP2 732 END IF 733 LWRK1 = LWORK - KEND1 734 IF (LWRK1 .LT. NTAMP) 735 * CALL QUIT('INSUFFICIENT WORK SPACE IN CCLR_JACAN') 736C 737C-------------------------------------------------------------- 738C Put up nc1vec C1 unit vectors and nc2vec C2 unit vectors. 739C-------------------------------------------------------------- 740C 741 KC1 = KJACOBI 742 KC2 = KC1 + NTAMP*NC1VEC 743 KC12 = KC2 + NTAMP*NC2VEC 744 CALL CCLR_UNIT(NTAMP,NC1VEC,NC2VEC, 745 & WORK(KC1),WORK(KC2),WORK(KC12)) 746 747 IF (CCR12) THEN 748 ! for R12 initialize Metric with the same unit vectors: 749 CALL DCOPY(NTAMP2,WORK(KJACOBI),1,WORK(KMETRIC),1) 750 END IF 751C 752C--------------------------------------------------------- 753C Make the transformation of the unit vectors. 754C Transformed vectors are returned as first 755C elements in work. 756C Work follows immediately after vectors as it should. 757C--------------------------------------------------------- 758C 759 NLOOP = NC1VEC + NC2VEC 760 IF (CCR12) NLOOP = NLOOP + NTR12AM(1) 761 762 DO I = 1, NLOOP 763 CALL DCOPY(NTAMP,WORK(KC1+NTAMP*(I-1)),1,WORK(KEND1),1) 764 CALL CC_ATRR(0.0D0,ISYMTR,ISIDE,WORK(KEND1),LWRK1, 765 & .FALSE.,DUMMY,APROXR12,TRIPLET) 766 CALL DCOPY(NTAMP,WORK(KEND1),1,WORK(KC1+NTAMP*(I-1)),1) 767 768 IF (CCR12) THEN 769 KS12AM = KEND1 + NTAMP 770 KM12AM = KMETRIC + NT1AM(ISYMTR)+NT2AM(ISYMTR)+NTAMP*(I-1) 771 CALL DCOPY(NTR12AM(ISYMTR),WORK(KS12AM),1,WORK(KM12AM),1) 772 END IF 773 END DO 774C 775C---------------------------------------------------------------- 776C Calculate norms to compare with finite difference jacobian. 777C---------------------------------------------------------------- 778C 779 XNJ = DDOT(NTAMP2,WORK(KJACOBI),1,WORK(KJACOBI),1) 780C 781 WRITE(LUPRI,*) ' ' 782 WRITE(LUPRI,*) ' NORM OF UNIT VEC. TRA. JAC.', SQRT(XNJ) 783 WRITE(LUPRI,*) ' ' 784C 785 CALL CCLR_NORMS(X11,X21,XR1,X12,X22,XR2,X1R,X2R,XRR, 786 & WORK(KJACOBI),NTAMP,NC1VEC,NC2VEC) 787C 788 WRITE(LUPRI,*) 'NORM OF 11 PART OF UNIT VECT. JAC. ', SQRT(X11) 789 WRITE(LUPRI,*) 'NORM OF 21 PART OF UNIT VECT. JAC. ', SQRT(X21) 790 WRITE(LUPRI,*) 'NORM OF R1 PART OF UNIT VECT. JAC. ', SQRT(XR1) 791 WRITE(LUPRI,*) 'NORM OF 12 PART OF UNIT VECT. JAC. ', SQRT(X12) 792 WRITE(LUPRI,*) 'NORM OF 22 PART OF UNIT VECT. JAC. ', SQRT(X22) 793 WRITE(LUPRI,*) 'NORM OF R2 PART OF UNIT VECT. JAC. ', SQRT(XR2) 794 WRITE(LUPRI,*) 'NORM OF 1R PART OF UNIT VECT. JAC. ', SQRT(X1R) 795 WRITE(LUPRI,*) 'NORM OF 2R PART OF UNIT VECT. JAC. ', SQRT(X2R) 796 WRITE(LUPRI,*) 'NORM OF RR PART OF UNIT VECT. JAC. ', SQRT(XRR) 797C 798C------------------------------------------------ 799C Print the columns of the jacobian obtained. 800C------------------------------------------------ 801C 802 IF (IPRINT .GE. 5) THEN 803 CALL AROUND( 'CC JACOBIANT I PRESUME - A*C1 PART ' ) 804 CALL OUTPUT(WORK(KJACOBI),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1, 805 & LUPRI) 806 CALL AROUND( 'CC JACOBIANT I PRESUME - A*C2 PART ' ) 807 CALL OUTPUT(WORK(KJACOBI+NTAMP*NC1VEC),1,NTAMP,1,NC2VEC, 808 * NTAMP,NC2VEC,1,LUPRI) 809 IF ( CCR12) THEN 810 CALL AROUND( 'CC JACOBIANT I PRESUME - A*C12 PART ' ) 811 CALL OUTPUT(WORK(KJACOBI+NTAMP*(NC1VEC+NC2VEC)),1,NTAMP, 812 * 1,NTR12AM(1),NTAMP,NTR12AM(1),1,LUPRI) 813 CALL AROUND( 'CC-R12 METRIC') 814 CALL OUTPUT(WORK(KMETRIC),1,NTAMP,1,NC1VEC+NC2VEC+NTR12AM(1), 815 * NTAMP,NC1VEC+NC2VEC+NTR12AM(1),1,LUPRI) 816 END IF 817 END IF 818 819 IF (IPRINT .GE. 2) THEN 820 CALL AROUND( 'FULL ANALYTICAL CC JACOBIANT' ) 821 IF (ISIDE.EQ.-1) THEN 822 CALL CCLR_TRANSSQ(WORK(KJACOBI),NTAMP) 823 ENDIF 824 if (.not. CCR12) then 825 call output(WORK(KJACOBI),1,NTAMP,1,NC1VEC+NC2VEC,NTAMP, 826 & NC1VEC+NC2VEC,1,LUPRI) 827 else 828 call output(WORK(KJACOBI),1,NTAMP,1,NC1VEC+NC2VEC+NTR12AM(1), 829 & NTAMP,NC1VEC+NC2VEC+NTR12AM(1),1,LUPRI) 830 end if 831 832 CALL AROUND(' END OF CCLR_JACAN ') 833 ENDIF 834C 835 RETURN 836 END 837C 838c*DECK CCLR_UNIT 839 SUBROUTINE CCLR_UNIT(NTAMP,NC1VEC,NC2VEC,C1AM,C2AM,C12AM) 840C 841C Put up the first NC1VEC and NC2VEC unit vectors in the single and double 842C excitation space respectively. 843C 844#include "implicit.h" 845#include "ccsdsym.h" 846#include "ccsdinp.h" 847#include "cclr.h" 848C 849 DIMENSION C1AM(NTAMP,NC1VEC),C2AM(NTAMP,NC2VEC) 850 DIMENSION C12AM(NTAMP,*) 851C 852C-------------------------------------- 853C Put up the required unit vectors. 854C-------------------------------------- 855C 856 DO J = 1, NC1VEC 857 CALL DZERO(C1AM(1,J),NTAMP) 858 C1AM(J,J) = 1.0D00 859 END DO 860 861 DO J = 1, NC2VEC 862 CALL DZERO(C2AM(1,J),NTAMP) 863 C2AM(J + NT1AM(ISYMTR),J) = 1.0D00 864 END DO 865 866 IF (CCR12) THEN 867 DO J = 1, NTR12AM(1) 868 CALL DZERO(C12AM(1,J),NTAMP) 869 C12AM(J + NT1AM(ISYMTR)+NT2AM(ISYMTR),J) = 1.0D00 870 END DO 871 END IF 872C 873 RETURN 874 END 875c*DECK CCLR_NORMS 876 SUBROUTINE CCLR_NORMS(X11,X21,XR1,X12,X22,XR2,X1R,X2R,XRR, 877 & XJAC,NTAMP,NC1VEC,NC2VEC) 878#include "implicit.h" 879#include "ccsdsym.h" 880#include "ccsdinp.h" 881#include "cclr.h" 882C 883 DIMENSION XJAC(NTAMP,*) 884C 885 X11 = 0.0D00 886 X12 = 0.0D00 887 X1R = 0.0D00 888 X21 = 0.0D00 889 X22 = 0.0D00 890 X2R = 0.0D00 891 XR1 = 0.0D00 892 XR2 = 0.0D00 893 XRR = 0.0D00 894C 895 DO 100 I = 1, NC1VEC 896C 897 X11 = X11 + DDOT(NT1AMX,XJAC(1,I),1,XJAC(1,I),1) 898 X21 = X21 + DDOT(NT2AMX,XJAC(1 + NT1AMX,I),1, 899 * XJAC(1 + NT1AMX,I),1) 900 IF (CCR12) THEN 901 XR1 = XR1 + DDOT(NTR12AM(1),XJAC(1+NT1AMX+NT2AMX,I),1, 902 * XJAC(1+NT1AMX+NT2AMX,I),1) 903 END IF 904C 905 100 CONTINUE 906C 907 DO 200 I = NC1VEC+1, NC1VEC + NC2VEC 908C 909 X12 = X12 + DDOT(NT1AMX,XJAC(1,I),1,XJAC(1,I),1) 910 X22 = X22 + DDOT(NT2AMX,XJAC(1 + NT1AMX,I),1, 911 * XJAC(1 + NT1AMX,I),1) 912 IF (CCR12) THEN 913 XR2 = XR2 + DDOT(NTR12AM(1),XJAC(1+NT1AMX+NT2AMX,I),1, 914 * XJAC(1+NT1AMX+NT2AMX,I),1) 915 END IF 916C 917 200 CONTINUE 918 919 IF (CCR12) THEN 920 DO I = NC1VEC+NC2VEC+1,NC1VEC+NC2VEC+NTR12AM(1) 921 X1R = X1R + DDOT(NT1AMX,XJAC(1,I),1,XJAC(1,I),1) 922 X2R = X2R + DDOT(NT2AMX,XJAC(1 + NT1AMX,I),1, 923 * XJAC(1 + NT1AMX,I),1) 924 XRR = XRR + DDOT(NTR12AM(1),XJAC(1+NT1AMX+NT2AMX,I),1, 925 * XJAC(1+NT1AMX+NT2AMX,I),1) 926 END DO 927 END IF 928C 929 RETURN 930 END 931c*DECK CCLR_DX 932 SUBROUTINE CCLR_DX(VEC,X,N) 933C 934C Written by Ove Christiansen 26-5-1994 935C 936#include "implicit.h" 937C 938 DIMENSION VEC(*) 939 DO 100 I = 1, N 940 VEC(I) = X 941100 CONTINUE 942 RETURN 943 END 944C----------------------------------------------------------------------- 945c*DECK CCLR_FDJAC 946 SUBROUTINE CCLR_FDJAC(NC1VEC,NC2VEC,WORK,LWORK,APROXR12) 947C----------------------------------------------------------------------- 948C Test routine for calculating the CC Jacobian by finite 949C difference on the CC vector function. 950C IF JACTST then first elements in work contains the jacobian. 951C calculates the first nc1vec rows of the C1 blocks 952C and nc2vec of the rows of the c2 blocks. 953C Written by Ove Christiansen 26-5-1994 954C Changes for R12 by Christof Haettig 21-5-2003 955C----------------------------------------------------------------------- 956 IMPLICIT NONE 957#include "maxorb.h" 958#include "ccorb.h" 959#include "ccsdsym.h" 960#include "priunit.h" 961#include "dummy.h" 962#include "second.h" 963#include "ccsdinp.h" 964#include "cclr.h" 965#include "r12int.h" 966#include "ccr12int.h" 967 968 INTEGER NC1VEC, NC2VEC, LWORK 969 970#if defined (SYS_CRAY) 971 REAL WORK(LWORK) 972 REAL XHALF, XMTWO, DELTA 973 REAL TI,X11,X12,X21,X22,XNJ,X1R,XR1,X2R,XR2,XRR,DDOT 974#else 975 DOUBLE PRECISION WORK(LWORK) 976 DOUBLE PRECISION XHALF, XMTWO, DELTA, DELTAI 977 DOUBLE PRECISION TI,X11,X12,X21,X22,XNJ,X1R,XR1,X2R,XR2,XRR,DDOT 978#endif 979 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, 980C & DELTA = 1.0D-05, DELTAI = 1.0D00/DELTA) 981 & DELTA = 1.0D-07, DELTAI = 1.0D00/DELTA) 982 CHARACTER*10 MODEL,MODELR12 983 CHARACTER*8 FC12AM,FS12AM,FC2AM,FS2AM 984 CHARACTER*(*) APROXR12 985 PARAMETER (FC12AM='CCR_C12M',FS12AM='CCR_S12M') 986 PARAMETER (FC2AM='CCR_C2AM',FS2AM='CCR_S2DM') 987 988 INTEGER NTAMP,NTAMP2,KJACOBI,KJACCO,KRHO1,KRHO2,KRHO1D,KRHO2D, 989 & KC1AM,KC2AM,KEND1,LWRK1,KJACOBI2,IOPT,LUJAC, KJACOBIR, 990 & NAI, NBJ, NKI, NLJ, KCR12AM, KRHOR12, KRHOR12D, 991 & LUNIT, LUFC12, LUFS12, LUMET, 992 & LUFS2,LUFC2,LUMET2,LUMET3,KC2AMPCK,KSCR1 993C 994 INTEGER INDEX 995 INDEX(I,J) = MAX(I,J)*(MAX(I,J)-3)/2 + I + J 996C 997 CALL QENTER('FDJAC') 998 999 LUFC2 = -1 1000 LUFS2 = -1 1001 LUFS12 = -1 1002 LUFC12 = -1 1003 1004 IF (OMEGSQ) CALL QUIT('CCLR_FDJAC incompatible with OMEGSQ.') 1005 IF (LCOR.OR.LSEC) 1006 & CALL QUIT('CCLR_FDJAC incompatible with FCORE and FSECON.') 1007 IF (NSYM.NE.1) CALL QUIT('CCLR_FDJAC must run in C1 symmetry.') 1008C 1009 IF (IPRINT.GT.5) THEN 1010 CALL AROUND( 'IN CCLR_FDJAC: MAKING FINITE DIFF. CC JACOBIANT') 1011 WRITE(LUPRI,*) 'Number of single excitations:',NC1VEC 1012 WRITE(LUPRI,*) 'Number of double excitations:',NC2VEC 1013 IF (CCR12) THEN 1014 WRITE(LUPRI,*) 'This is a CC-R12 calculation...' 1015 WRITE(LUPRI,*) 'Number of R12 double excitations:',NTR12AM(1) 1016 END IF 1017 ENDIF 1018C 1019C-------------------------------------------------- 1020C Write the jacobian to disk if FDJAC is done. 1021C-------------------------------------------------- 1022C 1023 IF (FDJAC ) THEN 1024 LUJAC = -1 1025 CALL GPOPEN(LUJAC,'CCLR_JAC','UNKNOWN',' ','UNFORMATTED', 1026 * IDUMMY,.FALSE.) 1027 IF (IPRINT.GT.5) WRITE(LUPRI,'(3X,A,//)') 1028 & 'FINITE DIFF. JACOBIAN WILL BE WRITTEN TO DISK' 1029 IF (CCR12) THEN 1030 LUMET = -1 1031 CALL GPOPEN(LUMET,'CCLR_MET','UNKNOWN',' ','UNFORMATTED', 1032 * IDUMMY,.FALSE.) 1033 IF (IPRINT.GT.5) WRITE(LUPRI,'(3X,A,//)') 1034 & 'METRIC MATRIX WILL BE WRITTEN TO DISK' 1035 CALL WOPEN2(LUFC12,FC12AM,64,0) 1036 CALL WOPEN2(LUFS12,FS12AM,64,0) 1037 IF (IANR12.EQ.2) THEN 1038 LUMET2 = -1 1039 CALL GPOPEN(LUMET2,'CCLR_MET2','UNKNOWN',' ','UNFORMATTED', 1040 * IDUMMY,.FALSE.) 1041 LUMET3 = -1 1042 CALL GPOPEN(LUMET3,'CCLR_MET3','UNKNOWN',' ','UNFORMATTED', 1043 * IDUMMY,.FALSE.) 1044 CALL WOPEN2(LUFS2,FS2AM,64,0) 1045 CALL WOPEN2(LUFC2,FC2AM,64,0) 1046 END IF 1047 END IF 1048 ENDIF 1049C 1050C---------------------------- 1051C Work space allocations. 1052C---------------------------- 1053C 1054 ISYMTR = 1 1055 IF (CCR12) THEN 1056 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) + NTR12AM(1) 1057 NTAMP2 = NTAMP*(NC1VEC + NC2VEC + NTR12AM(1)) 1058 ELSE 1059 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 1060 NTAMP2 = NTAMP*(NC1VEC + NC2VEC) 1061 END IF 1062 1063 KEND1 = 1 1064 1065 IF (JACTST) THEN 1066 KJACOBI = KEND1 1067 KEND1 = KJACOBI + NTAMP2 1068 ELSE 1069 KJACOBI = -999999 1070 ENDIF 1071 1072 IF (FDJAC ) THEN 1073 KJACCO = KEND1 1074 KEND1 = KJACCO + NTAMP 1075 ELSE 1076 KJACCO = -999999 1077 ENDIF 1078 1079 KRHO1 = KEND1 1080 KRHO2 = KRHO1 + NT1AMX 1081 KC1AM = KRHO2 + MAX(NT2AMX,NT2AO(ISYMOP),2*NT2ORT(ISYMOP)) 1082 KRHO1D = KC1AM + NT1AMX 1083 KRHO2D = KRHO1D + NT1AM(ISYMOP) 1084 KC2AM = KRHO2D + MAX(NT2AMX,NT2AO(ISYMOP),2*NT2ORT(ISYMOP)) 1085 KEND1 = KC2AM + MAX(NT2SQ(ISYMOP),NT2R12(1)) 1086 1087 IF (CCR12) THEN 1088 KCR12AM = KEND1 1089 KRHOR12 = KCR12AM + NTR12AM(1) 1090 KRHOR12D = KRHOR12 + NTR12AM(1) 1091 KEND1 = KRHOR12D + NTR12AM(1) 1092 IF (IANR12.EQ.2) THEN 1093 KC2AMPCK = KEND1 1094 KSCR1 = KC2AMPCK + NT2AM(ISYMTR) 1095 KEND1 = KSCR1 + MAX(NT2AM(ISYMTR),NTR12AM(ISYMTR)) 1096 END IF 1097 ELSE 1098 KCR12AM = -999999 1099 KRHOR12 = -999999 1100 KRHOR12D = -999999 1101 ENDIF 1102 1103 LWRK1 = LWORK - KEND1 1104 1105 IF (LWRK1.LT.0) THEN 1106 WRITE (LUPRI,*) 'Too little work space in cc_jacfd ' 1107 WRITE (LUPRI,*) 'AVAILABLE: LWORK = ',LWORK 1108 WRITE (LUPRI,*) 'NEEDED (AT LEAST) = ',KEND1 1109 CALL QUIT('TOO LITTLE WORKSPACE IN CCLR_FDJAC ') 1110 ENDIF 1111 1112 IF (IPRINT .GT. 15) THEN 1113 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KJACOBI = ',KJACOBI 1114 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KJACCO = ',KJACCO 1115 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHO1 = ',KRHO1 1116 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHO2 = ',KRHO2 1117 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KC1AM = ',KC1AM 1118 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHO1D = ',KRHO1D 1119 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHO2D = ',KRHO2D 1120 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KC2AM = ',KC2AM 1121 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KCR12AM = ',KCR12AM 1122 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHOR12 = ',KRHOR12 1123 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KRHOR12D= ',KRHOR12D 1124 WRITE (LUPRI,*) ' IN CCLR_FDJAC: KEND1 = ',KEND1 1125 ENDIF 1126 1127 KJACOBI2 = KJACOBI + NC1VEC*NTAMP 1128 KJACOBIR = KJACOBI2 + NC2VEC*NTAMP 1129C 1130C--------------------- 1131C Initializations. 1132C--------------------- 1133C 1134 IF (JACTST) CALL DZERO(WORK(KJACOBI),NTAMP2) 1135 IF (FDJAC ) CALL DZERO(WORK(KJACCO),NTAMP) 1136 1137 X11 = 0.0D00 1138 X12 = 0.0D00 1139 X1R = 0.0D00 1140 X21 = 0.0D00 1141 X22 = 0.0D00 1142 X2R = 0.0D00 1143 XR1 = 0.0D00 1144 XR2 = 0.0D00 1145 XRR = 0.0D00 1146C 1147C------------------------------------------------ 1148C Read the CC reference amplitudes From disk. 1149C------------------------------------------------ 1150C 1151 IOPT = 3 1152 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM)) 1153C 1154C-------------------------------------- 1155C Calculate reference omega vector. 1156C-------------------------------------- 1157C 1158 RSPIM = .FALSE. 1159C 1160 IF (.FALSE.) THEN 1161 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 1162 & WORK(KC2AM), 1163 * WORK(KEND1),LWRK1,APROXR12) 1164C 1165 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1),1) 1166 CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2),1) 1167 IF (CCR12) THEN 1168 LUNIT = -1 1169 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1170 * IDUMMY,.FALSE.) 1171 READ(LUNIT) (WORK(KRHOR12+N), N = 0, NTR12AM(1)-1 ) 1172 CALL GPCLOSE(LUNIT,'KEEP') 1173 END IF 1174C 1175 IF (IPRINT.GT.125) THEN 1176 CALL AROUND( 'RHO1' ) 1177 CALL OUTPUT(WORK(KRHO1),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 1178 CALL AROUND( 'RHO2' ) 1179 CALL OUTPAK(WORK(KRHO2),NT1AMX,1,LUPRI) 1180 IF (CCR12) THEN 1181 CALL AROUND( 'RHO-R12' ) 1182 CALL OUTPAK(WORK(KRHOR12),NMATKI(1),1,LUPRI) 1183 END IF 1184 ENDIF 1185 END IF 1186C 1187C--------------------------------------------- 1188C Calculate Jacobiant by finite difference. 1189C--------------------------------------------- 1190C 1191 DO I = 1, NC1VEC 1192 TI = SECOND() 1193 1194C Compute Omega(t - 0.5 delta) 1195 WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) - 0.5D0*DELTA 1196 1197 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2),WORK(KC1AM), 1198 & WORK(KC2AM), 1199 * WORK(KEND1),LWRK1,APROXR12) 1200 IF (CCR12) THEN 1201 LUNIT = -1 1202 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1203 & IDUMMY,.FALSE.) 1204 READ(LUNIT) (WORK(KRHOR12+N), N = 0, NTR12AM(1)-1 ) 1205 CALL GPCLOSE(LUNIT,'KEEP') 1206 END IF 1207 1208C Compute Omega(t + 0.5 delta) 1209 WORK(KC1AM +I -1) = WORK(KC1AM +I -1 ) + DELTA 1210 1211 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D), 1212 & WORK(KC1AM),WORK(KC2AM), 1213 * WORK(KEND1),LWRK1,APROXR12) 1214 1215C Restore t amplitudes 1216 WORK(KC1AM +I -1) = WORK(KC1AM +I -1) - 0.5D0*DELTA 1217 1218 IF (IPRINT.GT.125) THEN 1219 CALL AROUND( 'RHO1D ' ) 1220 CALL OUTPUT(WORK(KRHO1D),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1,LUPRI) 1221 CALL AROUND( 'RHO2D ' ) 1222 CALL OUTPAK(WORK(KRHO2D),NT1AMX,1,LUPRI) 1223 ENDIF 1224 1225C Calculate [Omega(t+0.5delta)-(omega(t-0.5delta)]/delta 1226 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 1227 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 1228 CALL CCLR_DIASCL(WORK(KRHO2D),XHALF,ISYMTR) 1229 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 1230 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 1231 1232 IF (CCR12) THEN 1233 LUNIT = -1 1234 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1235 & IDUMMY,.FALSE.) 1236 READ(LUNIT) (WORK(KRHOR12D+N), N = 0, NTR12AM(1)-1 ) 1237 CALL GPCLOSE(LUNIT,'KEEP') 1238 IF (IPRINT.GT.125) THEN 1239 CALL AROUND( 'RHO-R12' ) 1240 CALL OUTPAK(WORK(KRHO2),NMATKI(1),1,LUPRI) 1241 END IF 1242 CALL DAXPY(NTR12AM(1),-1.0D0,WORK(KRHOR12),1, 1243 & WORK(KRHOR12D),1) 1244 CALL CCLR_DIASCLR12(WORK(KRHOR12D),BRASCL,1) 1245 CALL DSCAL(NTR12AM(1),DELTAI,WORK(KRHOR12D),1) 1246 END IF 1247 1248 IF (JACTST) THEN 1249 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 1250 * WORK(KJACOBI+NTAMP*(I-1)),1) 1251 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1252 * WORK(KJACOBI+NTAMP*(I-1)+NT1AMX),1) 1253 IF (CCR12) CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1254 * WORK(KJACOBI+NTAMP*(I-1)+NT1AMX+NT2AMX),1) 1255 ENDIF 1256 IF (FDJAC ) THEN 1257 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KJACCO),1) 1258 CALL DCOPY(NT2AMX,WORK(KRHO2D),1,WORK(KJACCO+NT1AMX),1) 1259 IF (CCR12) CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1260 * WORK(KJACCO+NT1AMX+NT2AMX),1) 1261 WRITE(LUJAC) (WORK(KJACCO +J-1), J = 1, NTAMP) 1262 ENDIF 1263 X11 = X11 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1264 X21 = X21 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1265 IF (CCR12) 1266 * XR1 = XR1 + DDOT(NTR12AM(1),WORK(KRHOR12D),1,WORK(KRHOR12D),1) 1267C 1268 TI = SECOND() - TI 1269 IF (IPRINT.GT.5 ) THEN 1270 WRITE(LUPRI,*) ' ' 1271 WRITE(LUPRI,*) 'FDJAC ROW NR. ',I,' DONE IN ',TI,' SEC.' 1272 ENDIF 1273C 1274 END DO 1275C 1276C---------------------------------------------------------------- 1277C Loop over T2 amplitudes. Take care of diagonal t2 elements 1278C is in a different convention in the energy code. 1279C Factor 1/2 from right , and factor 2 from left. 1280C---------------------------------------------------------------- 1281C 1282 DO NAI = 1, NT1AMX 1283 DO NBJ = 1, NAI 1284 I = INDEX(NAI,NBJ) 1285C 1286 IF (I.LE.NC2VEC) THEN 1287 TI = SECOND() 1288 1289C Compute Omega(t - 0.5 delta) 1290 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) - 0.5D0*DELTA 1291 1292 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2),WORK(KC1AM), 1293 * WORK(KC2AM),WORK(KEND1),LWRK1,APROXR12) 1294 IF (CCR12) THEN 1295 LUNIT = -1 1296 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1297 & IDUMMY,.FALSE.) 1298 READ(LUNIT) (WORK(KRHOR12+N), N = 0, NTR12AM(1)-1 ) 1299 CALL GPCLOSE(LUNIT,'KEEP') 1300 END IF 1301 1302C Compute Omega(t + 0.5 delta) 1303 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) + DELTA 1304 1305 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D), 1306 & WORK(KC1AM), 1307 * WORK(KC2AM),WORK(KEND1),LWRK1,APROXR12) 1308 1309C Restore t amplitudes 1310 WORK(KC2AM + I -1) = WORK(KC2AM+I -1) - 0.5D0*DELTA 1311C 1312 IF (IPRINT.GT.125) THEN 1313 CALL AROUND( 'RHO1D ' ) 1314 CALL OUTPUT(WORK(KRHO1D),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1, 1315 & LUPRI) 1316 CALL AROUND( 'RHO2D ' ) 1317 CALL OUTPAK(WORK(KRHO2D),NT1AMX,1,LUPRI) 1318 ENDIF 1319C 1320C Calculate [Omega(t+0.5delta)-(omegat-0.5delta)]/delta 1321 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 1322 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 1323 CALL CCLR_DIASCL(WORK(KRHO2D),XHALF,ISYMTR) 1324 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 1325 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 1326 1327 IF (CCR12) THEN 1328 LUNIT = -1 1329 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ', 1330 & 'UNFORMATTED',IDUMMY,.FALSE.) 1331 READ(LUNIT) (WORK(KRHOR12D+N), N = 0, NTR12AM(1)-1 ) 1332 CALL GPCLOSE(LUNIT,'KEEP') 1333 IF (IPRINT.GT.125) THEN 1334 CALL AROUND( 'RHO-R12' ) 1335 CALL OUTPAK(WORK(KRHO2),NMATKI(1),1,LUPRI) 1336 END IF 1337 CALL DAXPY(NTR12AM(1),-1.0D0,WORK(KRHOR12),1, 1338 * WORK(KRHOR12D),1) 1339 CALL CCLR_DIASCLR12(WORK(KRHOR12D),BRASCL,1) 1340 CALL DSCAL(NTR12AM(1),DELTAI,WORK(KRHOR12D),1) 1341 END IF 1342C 1343 IF (NAI.EQ.NBJ) THEN 1344 CALL DSCAL(NT1AMX,2.0D00,WORK(KRHO1D),1) 1345 CALL DSCAL(NT2AMX,2.0D00,WORK(KRHO2D),1) 1346 IF (CCR12) CALL DSCAL(NTR12AM(1),2.0D00,WORK(KRHOR12D),1) 1347 ENDIF 1348C 1349 IF (JACTST) THEN 1350 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 1351 * WORK(KJACOBI2+NTAMP*(I-1)),1) 1352 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1353 * WORK(KJACOBI2+NTAMP*(I-1)+NT1AMX),1) 1354 IF (CCR12) CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1355 * WORK(KJACOBI2+NTAMP*(I-1)+NT1AMX+NT2AMX),1) 1356 ENDIF 1357C 1358 IF (FDJAC ) THEN 1359 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KJACCO ),1) 1360 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1361 * WORK(KJACCO+NT1AMX),1) 1362 IF (CCR12) CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1363 * WORK(KJACCO+NT1AMX+NT2AMX),1) 1364 WRITE(LUJAC) (WORK(KJACCO +J-1), J = 1, NTAMP) 1365 IF (CCR12.AND.IANR12.EQ.2) THEN 1366c make unit vector for this element (0 1 0) and 1367c compute lower nondiagonal part for metric and 1368c store it on file CCLR_MET3 1369 CALL DZERO(WORK(KC2AMPCK),NT2AM(1)) 1370 WORK(KC2AMPCK+I-1) = 1.0d0 1371 CALL CC_WVEC(LUFC2,FC2AM,NT2AM(1),NT2AM(1),1, 1372 & WORK(KC2AMPCK)) 1373 CALL DZERO(WORK(KSCR1),NTR12AM(1)) 1374 CALL CC_WVEC(LUFC12,FC12AM,NTR12AM(1),NTR12AM(1),1, 1375 & WORK(KSCR1)) 1376 CALL CC_R12METRIC(1,BRASCL,KETSCL,WORK(KEND1),LWRK1, 1377 & FC2AM,LUFC2,FC12AM,LUFC12,FS12AM, 1378 & LUFS12,FS2AM,LUFS2,1,.FALSE.,DUMMY) 1379 CALL CC_RVEC(LUFS12,FS12AM,NTR12AM(1),NTR12AM(1),1, 1380 & WORK(KSCR1)) 1381 WRITE(LUMET3) (WORK(KSCR1 +J-1), J = 1, NTR12AM(1)) 1382 END IF 1383 ENDIF 1384C 1385 X12 = X12 + DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1386 X22 = X22 + DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1387 IF (CCR12) 1388 * XR2=XR2 + DDOT(NTR12AM(1),WORK(KRHOR12D),1,WORK(KRHOR12D),1) 1389 1390 TI = SECOND() - TI 1391 IF (IPRINT.GT.5 ) THEN 1392 WRITE(LUPRI,*) ' ' 1393 WRITE(LUPRI,*) 'FDJAC ROW NR. ',I+NT1AMX, 1394 * ' DONE IN ',TI,' SEC.' 1395 ENDIF 1396C 1397 ENDIF 1398C 1399 END DO 1400 END DO 1401C 1402C---------------------------------------------------------------- 1403C Loop over R12 doubles amplitudes. 1404C Take care of the diagonal elements, which are in a 1405C different convention than in the energy code: 1406C factor 1/2 from right , and factor 2 from left. 1407C (Analogous to the conventional doubles amplitudes t^ab_ij) 1408C---------------------------------------------------------------- 1409C 1410 IF (CCR12) THEN 1411 ! read R12 doubles amplitudes from file 1412 IOPT = 32 1413 CALL CC_RDRSP('R0 ',0,1,IOPT,MODELR12,DUMMY,WORK(KCR12AM)) 1414 1415 ! make copy from which correct amplitudes can be 1416 ! restored after finite difference calculation 1417 LUNIT = -1 1418 CALL GPOPEN(LUNIT,'CCR12_TEMP','UNKNOWN',' ','UNFORMATTED', 1419 & IDUMMY,.FALSE.) 1420 WRITE(LUNIT) (WORK(KCR12AM+N), N = 0, NTR12AM(1)-1 ) 1421 CALL GPCLOSE(LUNIT,'KEEP') 1422 1423 DO NKI = 1, NMATKI(1) 1424 DO NLJ = 1, NKI 1425 I = INDEX(NKI,NLJ) 1426C 1427 TI = SECOND() 1428 1429C Compute Omega(t - 0.5 delta) 1430 WORK(KCR12AM + I -1) = WORK(KCR12AM+I -1) - 0.5D0*DELTA 1431 IOPT = 32 1432 CALL CC_WRRSP('R0 ',0,1,IOPT,MODELR12,DUMMY,DUMMY, 1433 * WORK(KCR12AM),WORK(KEND1),LWRK1) 1434 1435 CALL CCRHSN(WORK(KRHO1),WORK(KRHO2),WORK(KC1AM), 1436 * WORK(KC2AM),WORK(KEND1),LWRK1,APROXR12) 1437 LUNIT = -1 1438 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1439 & IDUMMY,.FALSE.) 1440 READ(LUNIT) (WORK(KRHOR12+N), N = 0, NTR12AM(1)-1 ) 1441 CALL GPCLOSE(LUNIT,'KEEP') 1442 1443C Compute Omega(t + 0.5 delta) 1444 WORK(KCR12AM + I -1) = WORK(KCR12AM+I -1) + DELTA 1445 IOPT = 32 1446 CALL CC_WRRSP('R0 ',0,1,IOPT,MODELR12,DUMMY,DUMMY, 1447 * WORK(KCR12AM),WORK(KEND1),LWRK1) 1448 1449 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D), 1450 & WORK(KC1AM), 1451 * WORK(KC2AM),WORK(KEND1),LWRK1,APROXR12) 1452 1453C Restore t amplitudes 1454 WORK(KCR12AM + I -1) = WORK(KCR12AM+I -1) - 0.5D0*DELTA 1455C 1456 IF (IPRINT.GT.125) THEN 1457 CALL AROUND( 'RHO1D ' ) 1458 CALL OUTPUT(WORK(KRHO1D),1,NVIRT,1,NRHFT,NVIRT,NRHFT,1, 1459 & LUPRI) 1460 CALL AROUND( 'RHO2D ' ) 1461 CALL OUTPAK(WORK(KRHO2D),NT1AMX,1,LUPRI) 1462 ENDIF 1463 1464C Calculate [Omega(t+0.5delta)-(omegat-0.5delta)]/delta 1465 CALL DAXPY(NT1AMX,-1.0D00,WORK(KRHO1),1,WORK(KRHO1D),1) 1466 CALL DAXPY(NT2AMX,-1.0D00,WORK(KRHO2),1,WORK(KRHO2D),1) 1467 CALL CCLR_DIASCL(WORK(KRHO2D),XHALF,ISYMTR) 1468 CALL DSCAL(NT1AMX,DELTAI,WORK(KRHO1D),1) 1469 CALL DSCAL(NT2AMX,DELTAI,WORK(KRHO2D),1) 1470 1471 LUNIT = -1 1472 CALL GPOPEN(LUNIT,'CC_OMEGAR12','UNKNOWN',' ','UNFORMATTED', 1473 & IDUMMY,.FALSE.) 1474 READ(LUNIT) (WORK(KRHOR12D+N), N = 0, NTR12AM(1)-1 ) 1475 CALL GPCLOSE(LUNIT,'KEEP') 1476 IF (IPRINT.GT.125) THEN 1477 CALL AROUND( 'RHO-R12' ) 1478 CALL OUTPAK(WORK(KRHO2),NMATKI(1),1,LUPRI) 1479 END IF 1480 CALL DAXPY(NTR12AM(1),-1.0D0,WORK(KRHOR12),1, 1481 * WORK(KRHOR12D),1) 1482 CALL CCLR_DIASCLR12(WORK(KRHOR12D),BRASCL,1) 1483 CALL DSCAL(NTR12AM(1),DELTAI,WORK(KRHOR12D),1) 1484C 1485 IF (NKI.EQ.NLJ) THEN 1486 CALL DSCAL(NT1AMX,KETSCL,WORK(KRHO1D),1) 1487 CALL DSCAL(NT2AMX,KETSCL,WORK(KRHO2D),1) 1488 IF (CCR12) CALL DSCAL(NTR12AM(1),KETSCL,WORK(KRHOR12D),1) 1489 ENDIF 1490C 1491 IF (JACTST) THEN 1492 CALL DCOPY(NT1AMX,WORK(KRHO1D),1, 1493 * WORK(KJACOBIR+NTAMP*(I-1)),1) 1494 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1495 * WORK(KJACOBIR+NTAMP*(I-1)+NT1AMX),1) 1496 IF (CCR12) CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1497 * WORK(KJACOBIR+NTAMP*(I-1)+NT1AMX+NT2AMX),1) 1498 ENDIF 1499C 1500 IF (FDJAC ) THEN 1501 CALL DCOPY(NT1AMX,WORK(KRHO1D),1,WORK(KJACCO ),1) 1502 CALL DCOPY(NT2AMX,WORK(KRHO2D),1, 1503 * WORK(KJACCO+NT1AMX),1) 1504 CALL DCOPY(NTR12AM(1),WORK(KRHOR12D),1, 1505 * WORK(KJACCO+NT1AMX+NT2AMX),1) 1506 WRITE(LUJAC) (WORK(KJACCO +J-1), J = 1, NTAMP) 1507 1508 ! make unit vector for this element and compute 1509 ! S * R_12 and store on file CCLR_MET: 1510 CALL DZERO(WORK(KRHOR12D),NTR12AM(1)) 1511 WORK(KRHOR12D + I -1) = 1.0D0 1512 CALL CC_WVEC(LUFC12,FC12AM,NTR12AM(1),NTR12AM(1),1, 1513 & WORK(KRHOR12D)) 1514 IF (IANR12.EQ.2) THEN 1515 CALL DZERO(WORK(KC2AMPCK),NT2AM(1)) 1516 CALL CC_WVEC(LUFC2,FC2AM,NT2AM(1),NT2AM(1),1, 1517 & WORK(KC2AMPCK)) 1518 END IF 1519 CALL CC_R12METRIC(1,BRASCL,KETSCL,WORK(KEND1),LWRK1, 1520 & FC2AM,LUFC2,FC12AM,LUFC12,FS12AM, 1521 & LUFS12,FS2AM,LUFS2,1,.FALSE.,DUMMY) 1522 CALL CC_RVEC(LUFS12,FS12AM,NTR12AM(1),NTR12AM(1),1, 1523 & WORK(KRHOR12D)) 1524 IF (IANR12.EQ.2) THEN 1525 CALL CC_RVEC(LUFS2,FS2AM,NT2AM(1),NT2AM(1),1, 1526 & WORK(KC2AMPCK)) 1527c write mixed S matrix on file 1528 WRITE(LUMET2) (WORK(KC2AMPCK +J-1), J = 1, NT2AM(1)) 1529 END IF 1530 1531 WRITE(LUMET) (WORK(KRHOR12D +J-1), J = 1, NTR12AM(1)) 1532 ENDIF 1533C 1534 X1R = X1R +DDOT(NT1AMX,WORK(KRHO1D),1,WORK(KRHO1D),1) 1535 X2R = X2R +DDOT(NT2AMX,WORK(KRHO2D),1,WORK(KRHO2D),1) 1536 XRR = XRR +DDOT(NTR12AM(1),WORK(KRHOR12D),1,WORK(KRHOR12D),1) 1537 1538 TI = SECOND() - TI 1539 IF (IPRINT.GT.5 ) THEN 1540 WRITE(LUPRI,*) ' ' 1541 WRITE(LUPRI,*) 'FDJAC ROW NR. ',I+NT1AMX+NT2AMX, 1542 * ' DONE IN ',TI,' SEC.' 1543 ENDIF 1544C 1545 END DO 1546 END DO 1547 1548 ! restore R12 amplitudes 1549 LUNIT = -1 1550 CALL GPOPEN(LUNIT,'CCR12_TEMP','UNKNOWN',' ','UNFORMATTED', 1551 & IDUMMY,.FALSE.) 1552 READ(LUNIT) (WORK(KCR12AM+N), N = 0, NTR12AM(1)-1 ) 1553 CALL GPCLOSE(LUNIT,'DELETE') 1554 IOPT = 32 1555 CALL CC_WRRSP('R0 ',0,1,IOPT,MODELR12,DUMMY,DUMMY, 1556 * WORK(KCR12AM),WORK(KEND1),LWRK1) 1557 END IF 1558C 1559 WRITE(LUPRI,*) ' ' 1560 WRITE(LUPRI,*) '** FINITE DIFF WITH DELTA ',DELTA, '**' 1561 WRITE(LUPRI,*) ' ' 1562 IF (FDJAC) CALL GPCLOSE(LUJAC,'KEEP') 1563 IF (FDJAC .AND. CCR12) THEN 1564 CALL GPCLOSE(LUMET,'KEEP') 1565 CALL WCLOSE2(LUFC12,FC12AM,'DELETE') 1566 CALL WCLOSE2(LUFS12,FS12AM,'DELETE') 1567 IF (IANR12.EQ.2) THEN 1568 CALL WCLOSE2(LUFS2,FS2AM,'DELETE') 1569 CALL WCLOSE2(LUFC2,FC2AM,'DELETE') 1570 CALL GPCLOSE(LUMET2,'KEEP') 1571 CALL GPCLOSE(LUMET3,'KEEP') 1572 END IF 1573 END IF 1574 IF (IPRINT.GE.5 .AND. JACTST) THEN 1575 CALL AROUND( 'FINITE DIFF. CC JACOBIANT - singles PART ' ) 1576 CALL OUTPUT(WORK(KJACOBI),1,NTAMP,1,NC1VEC,NTAMP,NC1VEC,1, 1577 & LUPRI) 1578 CALL AROUND( 'FINITE DIFF. CC JACOBIANT - doubles PART ' ) 1579 CALL OUTPUT(WORK(KJACOBI2),1,NTAMP,1,NC2VEC, 1580 * NTAMP,NC2VEC,1,LUPRI) 1581 IF (CCR12) THEN 1582 CALL AROUND( 'FINITE DIFF. CC JACOBIANT - R12 doubles PART') 1583 CALL OUTPUT(WORK(KJACOBIR),1,NTAMP,1,NTR12AM(1), 1584 * NTAMP,NTR12AM(1),1,LUPRI) 1585 END IF 1586 ENDIF 1587 IF (JACTST) THEN 1588 XNJ = X11 + X12 + X1R + X21 + X22 + X2R + XR1 + XR2 + XRR 1589 WRITE(LUPRI,*) ' ' 1590 WRITE(LUPRI,*) ' NORM OF FIN. DIFF. JAC.', SQRT(XNJ) 1591 WRITE(LUPRI,*) ' ' 1592 WRITE(LUPRI,*) 'NORM OF 11 PART OF FIN. DIFF. JAC. ', SQRT(X11) 1593 WRITE(LUPRI,*) 'NORM OF 21 PART OF FIN. DIFF. JAC. ', SQRT(X21) 1594 WRITE(LUPRI,*) 'NORM OF R1 PART OF FIN. DIFF. JAC. ', SQRT(XR1) 1595 WRITE(LUPRI,*) 'NORM OF 12 PART OF FIN. DIFF. JAC. ', SQRT(X12) 1596 WRITE(LUPRI,*) 'NORM OF 22 PART OF FIN. DIFF. JAC. ', SQRT(X22) 1597 WRITE(LUPRI,*) 'NORM OF R2 PART OF FIN. DIFF. JAC. ', SQRT(XR2) 1598 WRITE(LUPRI,*) 'NORM OF 1R PART OF FIN. DIFF. JAC. ', SQRT(X1R) 1599 WRITE(LUPRI,*) 'NORM OF 2R PART OF FIN. DIFF. JAC. ', SQRT(X2R) 1600 WRITE(LUPRI,*) 'NORM OF RR PART OF FIN. DIFF. JAC. ', SQRT(XRR) 1601 ENDIF 1602C 1603C----------------------------------- 1604C restore intermediates on disk: 1605C----------------------------------- 1606C 1607 IOPT = 3 1608 RSPIM = .TRUE. 1609 CALL CC_RDRSP('R0',0,1,IOPT,MODEL,WORK(KC1AM),WORK(KC2AM)) 1610 CALL CCRHSN(WORK(KRHO1D),WORK(KRHO2D),WORK(KC1AM), 1611 & WORK(KC2AM), 1612 * WORK(KEND1),LWRK1,APROXR12) 1613C 1614 IF (IPRINT .GT. 10) THEN 1615 CALL AROUND(' END OF CCLR_FDJAC ') 1616 ENDIF 1617C 1618 CALL QEXIT('FDJAC') 1619C 1620 RETURN 1621 END 1622C----------------------------------------------------------------------- 1623c*DECK CCLR_DUMTRR 1624 SUBROUTINE CCLR_DUMTRR(CVEC1,RHO1,RHO2,RHO12,S12AM,S2AM, 1625 & WORK,LWORK) 1626C------------------------------------------------------------------- 1627C Written by Ove Christiansen 21-11-1994 1628C 1629C Makes the transformation from vectors onto the jacobian 1630C by using a finite difference jacobian read in from disk. 1631C 1632C------------------------------------------------------------------- 1633C 1634#include "implicit.h" 1635#include "dummy.h" 1636C 1637#include "priunit.h" 1638#include "ccorb.h" 1639#include "maxorb.h" 1640#include "ccsdsym.h" 1641#include "ccsdinp.h" 1642#include "cclr.h" 1643#include "ccsdio.h" 1644#include "aovec.h" 1645#include "leinf.h" 1646#include "r12int.h" 1647 PARAMETER ( TWO = 2.0D00,XHALF=0.5D00 ) 1648 CHARACTER*1 TRANS 1649 DIMENSION CVEC1(*),WORK(LWORK) 1650 DIMENSION RHO1(*),RHO2(*), RHO12(*), S12AM(*), S2AM(*) 1651C 1652 IF (IPRINT .GT. 10) THEN 1653 CALL AROUND(' START OF CCLR_DUMTRR ') 1654 ENDIF 1655C 1656C---------------------------------- 1657C Print transformed vectors. 1658C---------------------------------- 1659C 1660 IF (IPRINT.GT.10) THEN 1661 CALL AROUND(' CCLR_DUMTRR: vector to be transformed ') 1662 CALL CC_PRP(CVEC1,CVEC1(1+NT1AM(ISYMTR)),ISYMTR,1,1) 1663 CALL AROUND('R12 double excitation part of vector:') 1664 CALL OUTPUT(CVEC1(1+NT1AM(ISYMTR)+NT2AM(ISYMTR)),1, 1665 * NTR12AM(ISYMTR),1,1,NTR12AM(ISYMTR),1,1,LUPRI) 1666 ENDIF 1667C 1668 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 1669 IF (CCR12) NTAMP = NTAMP + NTR12AM(ISYMTR) 1670C 1671C------------------------ 1672C Open Jacobian file. 1673C------------------------ 1674C 1675 LUJAC = -1 1676 CALL GPOPEN(LUJAC,'CCLR_JAC','UNKNOWN',' ','UNFORMATTED',IDUMMY, 1677 & .FALSE.) 1678 REWIND(LUJAC) 1679C 1680C--------------------------- 1681C Work space allocation. 1682C--------------------------- 1683C 1684 KJACCR= 1 1685 KEND1 = KJACCR+ NTAMP 1686 LWRK1 = LWORK - KEND1 1687 IF (LWRK1 .LT. 0) 1688 & CALL QUIT(' INSUFFICIENT WORK SPACE IN CCLR_DUMTRR ') 1689C 1690C--------------------------------------------------- 1691C Transform vectors by one column/row at a time. 1692C--------------------------------------------------- 1693C 1694 IF (IPRINT.GT.110) THEN 1695 CALL AROUND( 'PRINTOUT FROM CCLR_DUMTRR') 1696 ENDIF 1697 DO 100 I = 1, NT1AM(ISYMTR) 1698C 1699 IF (NSIDE.EQ.-1) THEN 1700C 1701 READ(LUJAC) (WORK(KJACCR+K-1),K=1,NTAMP) 1702 IF (IPRINT.GT.110) THEN 1703 WRITE (LUPRI,*) 1704 * 'THE ',I,'-th coulumn of the jacobian is ' 1705 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1706 ENDIF 1707C 1708 ELSE IF (NSIDE.EQ.1) THEN 1709C 1710 CALL CCLR_JACCR(WORK(KJACCR),I,LUJAC,WORK(KEND1),LWRK1) 1711 IF (IPRINT.GT.110) THEN 1712 WRITE (LUPRI,*) 'THE ',I,'-th row of the jacobian is ' 1713 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1714 ENDIF 1715C 1716 ENDIF 1717C 1718 RHO1(I) = DDOT(NTAMP,WORK(KJACCR),1,CVEC1,1) 1719C 1720 100 CONTINUE 1721C 1722 DO 200 I = 1, NT2AM(ISYMTR) 1723C 1724 IF (NSIDE.EQ.-1) THEN 1725C 1726 READ(LUJAC) (WORK(KJACCR+K-1),K=1,NTAMP) 1727 IF (IPRINT.GT.110) THEN 1728 WRITE (LUPRI,*) 'THE',I+NT1AM(ISYMTR), 1729 * '-th coulumn of the jacobian ' 1730 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1731 ENDIF 1732C 1733 ELSE IF (NSIDE.EQ.1) THEN 1734C 1735 CALL CCLR_JACCR(WORK(KJACCR),I+NT1AMX,LUJAC,WORK(KEND1), 1736 * LWRK1) 1737 IF (IPRINT.GT.110) THEN 1738 WRITE (LUPRI,*) 'THE ',I+NT1AM(ISYMTR), 1739 * '-th row of the jacobian ' 1740 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1741 ENDIF 1742C 1743 ENDIF 1744C 1745 RHO2(I) = DDOT(NTAMP,WORK(KJACCR),1,CVEC1,1) 1746C 1747 200 CONTINUE 1748 1749 IF (CCR12) THEN 1750 DO I = 1, NTR12AM(ISYMTR) 1751 IF (NSIDE.EQ.-1) THEN 1752 READ(LUJAC) (WORK(KJACCR+K-1),K=1,NTAMP) 1753 IF (IPRINT.GT.110) THEN 1754 WRITE (LUPRI,*) 'THE',I+NT1AM(ISYMTR)+NT2AM(ISYMTR), 1755 * '-th coulumn of the jacobian ' 1756 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1757 ENDIF 1758 ELSE IF (NSIDE.EQ.1) THEN 1759 CALL CCLR_JACCR(WORK(KJACCR),I+NT1AMX+NT2AMX,LUJAC, 1760 * WORK(KEND1),LWRK1) 1761 IF (IPRINT.GT.110) THEN 1762 WRITE (LUPRI,*) 'THE ',I+NT1AM(ISYMTR)+NT2AM(ISYMTR), 1763 * '-th row of the jacobian ' 1764 CALL OUTPUT(WORK(KJACCR),1,NTAMP,1,1,NTAMP,1,1,LUPRI) 1765 ENDIF 1766 ENDIF 1767 RHO12(I) = DDOT(NTAMP,WORK(KJACCR),1,CVEC1,1) 1768 END DO 1769 1770 KSMAT = 1 1771 KEND1 = KSMAT + NTR12AM(ISYMTR)*NTR12AM(ISYMTR) 1772 IF (IANR12.EQ.2) THEN 1773 KS2MAT = KEND1 1774 KEND1 = KS2MAT + NT2AM(ISYMTR)*NT2AM(ISYMTR) 1775 END IF 1776 LWRK1 = LWORK - KEND1 1777 IF (LWRK1 .LT. 0) 1778 & CALL QUIT(' INSUFFICIENT WORK SPACE IN CCLR_DUMTRR ') 1779 1780 LUMET = -1 1781 CALL GPOPEN(LUMET,'CCLR_MET','UNKNOWN',' ','UNFORMATTED', 1782 * IDUMMY,.FALSE.) 1783 CALL DZERO(WORK(KSMAT),NTR12AM(ISYMTR)*NTR12AM(ISYMTR)) 1784 DO I = 1, NTR12AM(1) 1785 IOFFS = KSMAT-1+NTR12AM(ISYMTR)*(I-1) 1786 READ(LUMET) (WORK(IOFFS+J),J=1,NTR12AM(ISYMTR)) 1787 END DO 1788 CALL GPCLOSE(LUMET,'KEEP') 1789 1790 IF (IANR12.EQ.2) THEN 1791 LUMET2 = -1 1792 CALL GPOPEN(LUMET2,'CCLR_MET2','UNKNOWN',' ','UNFORMATTED', 1793 * IDUMMY,.FALSE.) 1794 CALL DZERO(WORK(KS2MAT),NT2AM(ISYMTR)*NT2AM(ISYMTR)) 1795 DO I = 1, NT2AM(ISYMTR) 1796 IOFFS2 = KS2MAT -1+NT2AM(ISYMTR)*(I-1) 1797 READ(LUMET2)(WORK(IOFFS2+J),J=1,NT2AM(ISYMTR)) 1798 END DO 1799 CALL GPCLOSE(LUMET2,'KEEP') 1800 END IF 1801 1802 IF (NSIDE.EQ.1) THEN 1803 TRANS = 'N' 1804 ELSE IF (NSIDE.EQ.-1) THEN 1805 TRANS = 'T' 1806 ELSE 1807 CALL QUIT('Illegal NSIDE in CCLR_DUMTRR.') 1808 END IF 1809 KR12 = NT1AM(ISYMTR)+NT2AM(ISYMTR)+1 1810 CALL DGEMV(TRANS,NTR12AM(ISYMTR),NTR12AM(ISYMTR),1.0D0, 1811 * WORK(KSMAT),NTR12AM(ISYMTR),CVEC1(KR12),1,0.0D0, 1812 * S12AM,1) 1813 IF (IANR12.EQ.2) THEN 1814 KR12 = NT1AM(ISYMTR)+1 1815 CALL DGEMV(TRANS,NT2AM(ISYMTR),NT2AM(ISYMTR),1.0D0, 1816 & WORK(KS2MAT),NT2AM(ISYMTR),CVEC1(KR12),1,0.0D0, 1817 & S2AM,1) 1818 1819 END IF 1820 END IF 1821C 1822C---------------------------------- 1823C Print transformed vectors. 1824C---------------------------------- 1825C 1826 IF (IPRINT .GT. 10) THEN 1827 CALL AROUND(' CCLR_DUMTRR: transformed vectors ') 1828 CALL CC_PRP(RHO1,RHO2,ISYMTR,1,1) 1829 CALL AROUND('R12 double excitation part of vector:') 1830 CALL OUTPUT(RHO12,1,NTR12AM(ISYMTR),1,1, 1831 * NTR12AM(ISYMTR),1,1,LUPRI) 1832 CALL AROUND('R12 double excitation part of metric:') 1833 CALL OUTPUT(S12AM,1,NTR12AM(ISYMTR),1,1, 1834 * NTR12AM(ISYMTR),1,1,LUPRI) 1835 ENDIF 1836C 1837 IF (IPRINT .GT. 10) THEN 1838 CALL AROUND(' END OF CCLR_DUMTRR ') 1839 ENDIF 1840C 1841 CALL GPCLOSE(LUJAC,'KEEP') 1842C 1843 RETURN 1844 END 1845c*DECK CCLR_JACCR 1846 SUBROUTINE CCLR_JACCR(XJACCR,I,LU,WORK,LWORK) 1847C 1848C--------------------------------------------------------------- 1849C Reads the I-th row of the jacobian stored in the file 1850C specified by unit number LU. 1851C---------------------------------------------------------------- 1852C 1853#include "implicit.h" 1854C 1855#include "priunit.h" 1856#include "ccorb.h" 1857#include "maxorb.h" 1858#include "ccsdsym.h" 1859#include "ccsdinp.h" 1860#include "cclr.h" 1861#include "ccsdio.h" 1862#include "aovec.h" 1863C 1864 DIMENSION XJACCR(*),WORK(LWORK) 1865C 1866 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 1867 IF (CCR12) NTAMP = NTAMP + NTR12AM(ISYMTR) 1868C 1869C--------------------------- 1870C Work space allocation. 1871C--------------------------- 1872C 1873 KROW = 1 1874 KEND1 = KROW + NTAMP 1875 LWRK1 = LWORK - KEND1 1876 IF (LWRK1 .LT. 0) 1877 & CALL QUIT('INSUFFICIENT WORK SPACE IN CCLR_JACCR') 1878C 1879C--------------------------------------------------- 1880C Read in vectors by one column at a time. 1881C--------------------------------------------------- 1882C 1883 REWIND(LU) 1884C 1885 DO 100 J = 1, NTAMP 1886C 1887 READ(LU) (WORK(KROW +K -1),K= 1, NTAMP) 1888C 1889 XJACCR(J) = WORK(KROW + I -1) 1890C 1891 100 CONTINUE 1892C 1893 RETURN 1894 END 1895C /* Deck cclr_lamtra */ 1896 SUBROUTINE CCLR_LAMTRA(XLAMDP,XLAMPC,XLAMDH,XLAMHC,C1AM,ISYMTR) 1897C 1898C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1899C Written by Ove Christiansen 10-2-1995 1900C 1901C PURPOSE: 1902C transform lamda matrices wtih C1AM 1903C debugged 10-2-1995 1904C 1905C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1906C 1907#include "implicit.h" 1908#include "iratdef.h" 1909 PARAMETER(ZERO = 0.0D0, ONE = 1.0D0, XMONE= -1.0D00) 1910 DIMENSION C1AM(*),XLAMDP(*),XLAMDH(*),XLAMPC(*),XLAMHC(*) 1911#include "ccorb.h" 1912#include "ccsdsym.h" 1913#include "ccsdinp.h" 1914C 1915 IF (IPRINT .GT.25) THEN 1916 CALL AROUND('IN CCLR_LAMTRA ') 1917 ENDIF 1918C 1919C----------------------------------------- 1920C Transform lamda particle matrix. 1921C LaP~(al,a) = -sum(k)[ La(al,k)*C(a,k)] 1922C NB!! notice the minus sign. 1923C----------------------------------------- 1924C 1925 CALL DZERO(XLAMPC,NGLMDT(ISYMTR)) 1926C 1927 DO 100 ISYMA = 1,NSYM 1928C 1929 ISYMK = MULD2H(ISYMTR,ISYMA) 1930 ISYMAL = ISYMK 1931C 1932 NBASAL = MAX(NBAS(ISYMAL),1) 1933 NBASA = MAX(NVIR(ISYMA),1) 1934C 1935 KOFF1 = IGLMRH(ISYMAL,ISYMK) + 1 1936 KOFF2 = IT1AM(ISYMA,ISYMK) + 1 1937 KOFF3 = IGLMVI(ISYMAL,ISYMA) + 1 1938C 1939 CALL DGEMM('N','T',NBAS(ISYMAL),NVIR(ISYMA),NRHF(ISYMK), 1940 * XMONE,XLAMDP(KOFF1),NBASAL,C1AM(KOFF2),NBASA, 1941 * ZERO,XLAMPC(KOFF3),NBASAL) 1942C 1943 100 CONTINUE 1944C 1945C----------------------------------------- 1946C Transform lamda hole matrix. 1947C LaH~(al,i) = sum(c)[ La(al,c)*C(c,i)] 1948C----------------------------------------- 1949C 1950 CALL DZERO(XLAMHC,NGLMDT(ISYMTR)) 1951C 1952 DO 200 ISYMI = 1,NSYM 1953C 1954 ISYMC = MULD2H(ISYMTR,ISYMI) 1955 ISYMAL = ISYMC 1956C 1957 NBASAL = MAX(NBAS(ISYMAL),1) 1958 NBASC = MAX(NVIR(ISYMC),1) 1959C 1960 KOFF1 = IGLMVI(ISYMAL,ISYMC) + 1 1961 KOFF2 = IT1AM(ISYMC,ISYMI) + 1 1962 KOFF3 = IGLMRH(ISYMAL,ISYMI) + 1 1963C 1964 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI),NVIR(ISYMC), 1965 * ONE,XLAMDH(KOFF1),NBASAL,C1AM(KOFF2),NBASC, 1966 * ZERO,XLAMHC(KOFF3),NBASAL) 1967C 1968 200 CONTINUE 1969C 1970 RETURN 1971 END 1972C /* Deck cc_prtspv */ 1973 SUBROUTINE CC_PRTSPV(T1AM,T2AM) 1974C 1975C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1976C Ove Christiansen 24-1-1994 1977C PRint Total Symmetric Packed Vector (PRTSPV) (T1,T2) 1978C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 1979C 1980#include "implicit.h" 1981#include "ccorb.h" 1982#include "ccsdsym.h" 1983#include "ccsdinp.h" 1984#include "priunit.h" 1985 DIMENSION T1AM(NT1AMX),T2AM(NT2AMX) 1986C 1987C------------------------------------ 1988C Write single excitation vector. 1989C------------------------------------ 1990C 1991 CALL AROUND('single excitation part of vector ') 1992 DO 300 ISYM = 1,NSYM 1993 WRITE(LUPRI,*) 'Symmetry block number : ',ISYM 1994 KOFF = IT1AM(ISYM,ISYM) + 1 1995 CALL OUTPUT(T1AM(KOFF),1,NVIR(ISYM),1,NRHF(ISYM), 1996 * NVIR(ISYM),NRHF(ISYM),1,LUPRI) 1997 WRITE(LUPRI,*) ' ' 1998 300 CONTINUE 1999C 2000C------------------------------------ 2001C Write double excitation vector. 2002C------------------------------------ 2003C 2004 CALL AROUND('double excitation part of vector ') 2005 DO 310 ISYM = 1,NSYM 2006 WRITE(LUPRI,*) 'Symmetry block number : ',ISYM 2007 KOFF = IT2AM(ISYM,ISYM) + 1 2008 CALL OUTPAK(T2AM(KOFF),NT1AM(ISYM),1,LUPRI) 2009 WRITE(LUPRI,*) ' ' 2010 310 CONTINUE 2011C 2012 RETURN 2013 END 2014C /* Deck cc_prp */ 2015 SUBROUTINE CC_PRP(T1AM,T2AM,ISYM,NT1,NT2) 2016C 2017C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2018C Ove Christiansen 24-1-1994 2019C PRint Packed Vector - PRP (in general non. tot. sym.) (t1,t2). 2020C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2021C 2022#include "implicit.h" 2023#include "ccorb.h" 2024#include "ccsdsym.h" 2025#include "ccsdinp.h" 2026#include "priunit.h" 2027C 2028 DIMENSION T1AM(*),T2AM(*) 2029C 2030C------------------------------------ 2031C set dimensions 2032C------------------------------------ 2033C 2034 LT1AM = NT1AM(ISYM) 2035 LT2AM = NT2AM(ISYM) 2036 2037C 2038C------------------------------------ 2039C Write single excitation vector. 2040C------------------------------------ 2041C 2042 DO 10 I = 1, NT1 2043C 2044 CALL AROUND('single excitation part of vector ') 2045 IF (NT1.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2046C 2047 DO 100 ISYMI = 1,NSYM 2048C 2049 ISYMA = MULD2H(ISYMI,ISYM) 2050C 2051 WRITE(LUPRI,*) 'Symmetry block number(a,i): ',ISYMA,ISYMI 2052 KOFF = LT1AM*(I-1) + IT1AM(ISYMA,ISYMI) + 1 2053 CALL OUTPUT(T1AM(KOFF),1,NVIR(ISYMA),1,NRHF(ISYMI), 2054 * NVIR(ISYMA),NRHF(ISYMI),1,LUPRI) 2055 WRITE(LUPRI,*) ' ' 2056C 2057 100 CONTINUE 2058C 2059 10 CONTINUE 2060C 2061C------------------------------------ 2062C Write double excitation vector. 2063C------------------------------------ 2064C 2065 DO 20 I = 1,NT2 2066C 2067 CALL AROUND('double excitation part of vector ') 2068 IF (NT2.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2069C 2070 DO 200 ISYMBJ = 1,NSYM 2071C 2072 ISYMAI = MULD2H(ISYMBJ,ISYM) 2073C 2074 KOFF = LT2AM*(I-1) + IT2AM(ISYMAI,ISYMBJ) + 1 2075 2076 IF (ISYMAI.EQ.ISYMBJ) THEN 2077C 2078 WRITE(LUPRI,*) 'Symmetry block number(ai,bj): ', 2079 & ISYMAI,ISYMBJ 2080 CALL OUTPAK(T2AM(KOFF),NT1AM(ISYMAI),1,LUPRI) 2081 WRITE(LUPRI,*) ' ' 2082C 2083 ELSE IF (ISYMBJ.GT.ISYMAI) THEN 2084C 2085 WRITE(LUPRI,*) 'Symmetry block number(ai,bj): ', 2086 & ISYMAI,ISYMBJ 2087 CALL OUTPUT(T2AM(KOFF),1,NT1AM(ISYMAI),1,NT1AM(ISYMBJ), 2088 * NT1AM(ISYMAI),NT1AM(ISYMBJ),1,LUPRI) 2089 WRITE(LUPRI,*) ' ' 2090C 2091 ENDIF 2092C 2093 200 CONTINUE 2094C 2095 20 CONTINUE 2096C 2097 RETURN 2098 END 2099C /* Deck cc_prpao */ 2100 SUBROUTINE CC_PRPAO(T1AO,T2AO,ISYM,NT1,NT2) 2101C 2102C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2103C Ove Christiansen 24-1-1994 2104C PRint Packed Vector - PRP (in general non. tot. sym.) (t1,t2). 2105C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2106C 2107#include "implicit.h" 2108#include "ccorb.h" 2109#include "ccsdsym.h" 2110#include "ccsdinp.h" 2111#include "priunit.h" 2112C 2113 DIMENSION T1AO(*),T2AO(*) 2114C 2115C------------------------------------ 2116C set dimensions 2117C------------------------------------ 2118C 2119 LT1AO = NT1AO(ISYM) 2120 LT2AO = NT2AO(ISYM) 2121 2122C 2123C------------------------------------ 2124C Write single excitation vector. 2125C------------------------------------ 2126C 2127 DO 10 I = 1, NT1 2128C 2129 CALL AROUND('single excitation part of vector ') 2130 IF (NT1.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2131C 2132 DO 100 ISYMI = 1,NSYM 2133C 2134 ISYMA = MULD2H(ISYMI,ISYM) 2135C 2136 WRITE(LUPRI,*) 'Symmetry block number(a,i): ',ISYMA,ISYMI 2137 KOFF = LT1AO*(I-1) + IT1AO(ISYMA,ISYMI) + 1 2138 CALL OUTPUT(T1AO(KOFF),1,NBAS(ISYMA),1,NRHF(ISYMI), 2139 * NBAS(ISYMA),NRHF(ISYMI),1,LUPRI) 2140 WRITE(LUPRI,*) ' ' 2141C 2142 100 CONTINUE 2143C 2144 10 CONTINUE 2145C 2146C------------------------------------ 2147C Write double excitation vector. 2148C------------------------------------ 2149C 2150 DO 20 I = 1,NT2 2151C 2152 CALL AROUND('double excitation part of vector ') 2153 IF (NT2.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2154C 2155 DO 200 ISYMBJ = 1,NSYM 2156C 2157 ISYMAI = MULD2H(ISYMBJ,ISYM) 2158C 2159 KOFF = LT2AO*(I-1) + IT2AO(ISYMAI,ISYMBJ) + 1 2160 2161 IF (ISYMAI.EQ.ISYMBJ) THEN 2162C 2163 WRITE(LUPRI,*) 'Symmetry block number(ai,bj): ', 2164 & ISYMAI,ISYMBJ 2165 CALL OUTPAK(T2AO(KOFF),NT1AO(ISYMAI),1,LUPRI) 2166 WRITE(LUPRI,*) ' ' 2167C 2168 ELSE IF (ISYMBJ.GT.ISYMAI) THEN 2169C 2170 WRITE(LUPRI,*) 'Symmetry block number(ai,bj): ', 2171 & ISYMAI,ISYMBJ 2172 CALL OUTPUT(T2AO(KOFF),1,NT1AO(ISYMAI),1,NT1AO(ISYMBJ), 2173 * NT1AO(ISYMAI),NT1AO(ISYMBJ),1,LUPRI) 2174 WRITE(LUPRI,*) ' ' 2175C 2176 ENDIF 2177C 2178 200 CONTINUE 2179C 2180 20 CONTINUE 2181C 2182 RETURN 2183 END 2184C /* Deck cc_prsq */ 2185 SUBROUTINE CC_PRSQ(T1AM,T2AM,ISYM,NT1,NT2) 2186C 2187C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2188C Ove Christiansen 24-1-1994 2189C PRint SQuared vector. general non. tot. sym. (t1,t2). 2190C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2191C 2192#include "implicit.h" 2193#include "ccorb.h" 2194#include "ccsdsym.h" 2195#include "ccsdinp.h" 2196#include "priunit.h" 2197C 2198 DIMENSION T1AM(*),T2AM(*) 2199 2200C 2201C------------------------------------ 2202C Set dimensions 2203C------------------------------------ 2204C 2205 LT1AM = NT1AM(ISYM) 2206 LT2AM = NT2SQ(ISYM) 2207C 2208C------------------------------------ 2209C Write single excitation vector. 2210C------------------------------------ 2211C 2212 DO 10 I = 1, NT1 2213C 2214 CALL AROUND('single excitation part of vector ') 2215 IF (NT1.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2216C 2217 DO 100 ISYMI = 1,NSYM 2218C 2219 ISYMA = MULD2H(ISYMI,ISYM) 2220C 2221 WRITE(LUPRI,*) 'Symmetry block number(a,i): ',ISYMA,ISYMI 2222 2223 KOFF = LT1AM*(I-1) + IT1AM(ISYMA,ISYMI) + 1 2224 CALL OUTPUT(T1AM(KOFF),1,NVIR(ISYMA),1,NRHF(ISYMI), 2225 * NVIR(ISYMA),NRHF(ISYMI),1,LUPRI) 2226 WRITE(LUPRI,*) ' ' 2227C 2228 100 CONTINUE 2229C 2230 10 CONTINUE 2231C 2232C------------------------------------ 2233C Write double excitation vector. 2234C------------------------------------ 2235C 2236 DO 20 I = 1, NT2 2237C 2238 CALL AROUND('double excitation part of vector ') 2239 IF (NT2.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 2240C 2241 DO 200 ISYMBJ = 1,NSYM 2242C 2243 ISYMAI = MULD2H(ISYMBJ,ISYM) 2244C 2245 WRITE(LUPRI,*) 'Symmetry block number(ai,bj): ', 2246 & ISYMAI,ISYMBJ 2247C 2248 KOFF = LT2AM*(I-1) + IT2SQ(ISYMAI,ISYMBJ) + 1 2249 CALL OUTPUT(T2AM(KOFF),1,NT1AM(ISYMAI),1,NT1AM(ISYMBJ), 2250 * NT1AM(ISYMAI),NT1AM(ISYMBJ),1,LUPRI) 2251C 2252 WRITE(LUPRI,*) ' ' 2253C 2254 200 CONTINUE 2255C 2256 20 CONTINUE 2257C 2258 RETURN 2259 END 2260c*DECK CCLR_FDEXCI 2261 SUBROUTINE CCLR_FDEXCI(WORK,LWORK) 2262C 2263C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2264C 2265C Purpose: Calculation of Excitation energies from 2266C explicit jacobian constructed from finite 2267C difference. 2268C 2269C Written by Ove Christiansen November 1994 2270C 2271C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2272C 2273#include "implicit.h" 2274#include "maxorb.h" 2275#include "ccorb.h" 2276#include "iratdef.h" 2277#include "codata.h" 2278#include "priunit.h" 2279#include "ccsdinp.h" 2280#include "cclr.h" 2281#include "ccsdsym.h" 2282#include "ccsdio.h" 2283c 2284#include "r12int.h" 2285C 2286 LOGICAL OPNSTA, LOCDBG 2287 DIMENSION WORK(LWORK) 2288 PARAMETER (THRZER = 1.0D-99) 2289C 2290#include "leinf.h" 2291C 2292C--------------------------------------------- 2293C Header of Excitation Energy calculation. 2294C--------------------------------------------- 2295C 2296 CALL QENTER('FDEXCI') 2297c 2298 CALL TIMER('START ',TIMEIN,TIMOUT) 2299 2300 IF (IPRINT .GT. 0) THEN 2301 WRITE (LUPRI,'(A,/)') ' ' 2302 WRITE (LUPRI,'(A,/)') 2303 * '1 <<<<<<<<<< Output from CCLR_FDEXCI >>>>>>>>> ' 2304 IF (CCR12) WRITE (LUPRI,'(A,/)') 'This is a R12 run.' 2305 WRITE (LUPRI,'(A,/)') ' ' 2306 ENDIF 2307C 2308C--------------------------------- 2309C Initialization of variables. 2310C--------------------------------- 2311C 2312 MATZ = 1 2313 NREDH = NT1AMX + NT2AMX 2314 NTAMP = NT1AMX + NT2AMX 2315 IF (CCR12) THEN 2316 NREDH = NREDH + NTR12AM(1) 2317 NTAMP = NTAMP + NTR12AM(1) 2318 END IF 2319C 2320C---------------------------- 2321C Alloction of Workspace. 2322C---------------------------- 2323C 2324 IF (NTAMP.GT.MAXRED .OR. NREDH.GT.MAXRED) THEN 2325 WRITE (LUPRI,*)'MAXRED, NTAMP, NREDH', MAXRED, NTAMP, NREDH 2326 CALL QUIT('MAXRED too small in CCLR_FDEXCI.') 2327 END IF 2328C 2329 KWI = 1 2330 KAMAT = KWI + MAXRED 2331 KIV1 = KAMAT+ MAXRED*MAXRED 2332 KFV1 = KIV1 + MAXRED 2333 KEND = KFV1 + MAXRED 2334 LEND = LWORK - KEND 2335 KEIVAL = KEND 2336 KSOLEQ = KEIVAL + MAXRED 2337 KAMAT2 = KSOLEQ + MAXRED*MAXRED 2338 KWRK1 = KAMAT2 + MAXRED*MAXRED 2339 IF (CCR12) THEN 2340 KDENOM = KWRK1 2341 KSMAT = KDENOM + MAXRED 2342 KWRK1 = KSMAT + MAXRED*MAXRED 2343 IF (IANR12.EQ.2) THEN 2344 KSEIVAL = KWRK1 2345 KSWI = KSEIVAL + MAXRED 2346 KSOL = KSWI + MAXRED 2347 KSOLEQS = KSOL + MAXRED*MAXRED 2348 KSCOPY = KSOLEQS + MAXRED*MAXRED 2349 KWRK1 = KSCOPY + MAXRED*MAXRED 2350 END IF 2351 END IF 2352 LWRK1 = LWORK - KWRK1 2353 IF (LWRK1.LT. 0) 2354 & CALL QUIT('TOO LITTLE WORK SPACE IN CCLR_FDEXCI ') 2355C 2356C--------------------- 2357C Readin Jacobian. 2358C--------------------- 2359C 2360 CALL CCLR_JACIN(WORK(KSOLEQ),WORK(KSMAT)) 2361C 2362 IF (IPRINT .GE. 10) THEN 2363 CALL AROUND( 'FIN. DIF. CC JACOBIANT - A*C1 PART ' ) 2364 CALL OUTPUT(WORK(KSOLEQ),1,NTAMP,1,NT1AMX,MAXRED,NTAMP,1,LUPRI) 2365 CALL AROUND( 'FIN. DIF. CC JACOBIANT - A*C2 PART ' ) 2366 CALL OUTPUT(WORK(KSOLEQ +MAXRED*NT1AMX),1,NTAMP,1,NT2AMX, 2367 * MAXRED,NT2AMX,1,LUPRI) 2368 IF (CCR12) THEN 2369 CALL AROUND( 'FIN. DIF. CC JACOBIANT - A*CR12 PART ' ) 2370 CALL OUTPUT(WORK(KSOLEQ +MAXRED*(NT1AMX+NT2AMX)),1,NTAMP, 2371 * 1,NTR12AM(1),MAXRED,NTR12AM(1),1,LUPRI) 2372 CALL AROUND( 'CC METRIC' ) 2373 CALL OUTPUT(WORK(KSMAT),1,NTAMP,1,NTAMP,MAXRED,NTAMP,1,LUPRI) 2374 END IF 2375 ENDIF 2376C 2377C--------------------------- 2378C Solve for eigenvalues. 2379C--------------------------- 2380C 2381 IF (CCR12) THEN 2382 CALL DZERO(WORK(KAMAT),MAXRED*MAXRED) 2383 CALL DCOPY(MAXRED*MAXRED,WORK(KSOLEQ),1,WORK(KAMAT),1) 2384 WRITE(LUPRI,'(/A)')'Jacobi matrix:' 2385 CALL OUTPUT(WORK(KAMAT),1,MAXRED,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2386 2387C lujacobi = 0 2388C call gpopen(lujacobi,'JACOBIMAT','UNKNOWN',' ','FORMATTED', 2389C & IDUMMY,.FALSE.) 2390 2391C do i = 1, nrhf(1) 2392C do a = 1, nvir(1) 2393C nai = nvir(1)*(i-1) + a 2394C do j = 1, nrhf(1) 2395C do b = 1, nvir(1) 2396C nbj = nvir(1)*(j-1) + b 2397C if (nbj.le.nai) then 2398C naibj = nai*(nai-1)/2 + nbj 2399C idxaibj = naibj + nvir(1)*nrhf(1) 2400 2401C do k = 1, nrhf(1) 2402CCN do m = 1, nrhf(1) 2403C do m = 1, nrhfb(1) 2404CCN nmk = nrhf(1)*(k-1) + m 2405C nmk = nrhfb(1)*(k-1) + m 2406C do l = 1, nrhf(1) 2407CCN do n = 1, nrhf(1) 2408C do n = 1, nrhfb(1) 2409CCN nnl = nrhf(1)*(l-1) + n 2410C nnl = nrhfb(1)*(l-1) + n 2411C if (nnl.le.nmk) then 2412C nmknl = nmk*(nmk-1)/2 + nnl 2413C idxmknl = nmknl + nvir(1)*nrhf(1) + 2414C & (nrhf(1)*nvir(1))*(nrhf(1)*nvir(1)+1)/2 2415 2416C idx23 = maxred*(idxaibj-1) + idxmknl 2417C idx32 = maxred*(idxmknl-1) + idxaibj 2418 2419C if ((dabs(work(kamat-1+idx23))+ 2420C & dabs(work(ksmat-1+idx23))+ 2421C & dabs(work(kamat-1+idx32))+ 2422C & dabs(work(ksmat-1+idx32))).gt.1.0d-8) then 2423 2424C write(lujacobi,'(8i2,4e20.10)') 2425C & i,a,j,b,k,m,l,n, 2426C & work(kamat-1+idx23), 2427C & work(ksmat-1+idx23), 2428C & work(kamat-1+idx32), 2429C & work(ksmat-1+idx32) 2430C end if 2431 2432C end if 2433C end do 2434C end do 2435C end do 2436C end do 2437C end if 2438C end do 2439C end do 2440C end do 2441C end do 2442 2443C call gpclose(lujacobi,'KEEP') 2444 2445c get idea of stability matrix 2446 call cc_r12stabmat(work(kamat),maxred,nt1am(1),nt2am(1), 2447 & ntr12am(1),work(kwrk1),lwrk1) 2448 2449 write(lupri,*)'A(2,3) Block:' 2450 call output(work(kamat),nt1am(1)+1, 2451 & nt1am(1)+nt2am(1),nt1am(1)+nt2am(1)+1, 2452 & nt1am(1)+nt2am(1)+ntr12am(1),maxred,maxred,1,lupri) 2453 write(lupri,*)'A(3,2) Block:' 2454 call output(work(kamat),nt1am(1)+nt2am(1)+1, 2455 & nt1am(1)+nt2am(1)+ntr12am(1),nt1am(1)+1, 2456 & nt1am(1)+nt2am(1),maxred,maxred,1,lupri) 2457c 2458 IF (IANR12.EQ.2) THEN 2459 WRITE(LUPRI,*)'S matrix:' 2460 CALL OUTPUT(WORK(KSMAT),1,NTAMP,1,NTAMP,MAXRED,NTAMP,1,LUPRI) 2461 CALL DCOPY(MAXRED*MAXRED,WORK(KSMAT),1,WORK(KSCOPY),1) 2462chf 2463 write(lupri,*)'S(2,3) Block:' 2464 call output(work(ksmat),nt1am(1)+1, 2465 & nt1am(1)+nt2am(1),nt1am(1)+nt2am(1)+1, 2466 & nt1am(1)+nt2am(1)+ntr12am(1),maxred,maxred,1,lupri) 2467 write(lupri,*)'S(3,2) Block:' 2468 call output(work(ksmat),nt1am(1)+nt2am(1)+1, 2469 & nt1am(1)+nt2am(1)+ntr12am(1),nt1am(1)+1, 2470 & nt1am(1)+nt2am(1),maxred,maxred,1,lupri) 2471c get eigenvalues of S matrix to check if it is positive definite 2472c caution rg routine destroys old S matrix 2473 CALL RG(MAXRED,NREDH,WORK(KSCOPY),WORK(KSEIVAL),WORK(KSWI), 2474 & MATZ,WORK(KSOLEQS),WORK(KIV1),WORK(KFV1),IERR) 2475 WRITE(LUPRI,'(/A)')'Eigenvalues of S matrix:' 2476 CALL OUTPUT(WORK(KSEIVAL),1,NREDH,1,1,NREDH,MAXRED,1,LUPRI) 2477 END IF 2478c 2479 CALL RGG(MAXRED,NREDH,WORK(KAMAT),WORK(KSMAT),WORK(KEIVAL), 2480 * WORK(KWI),WORK(KDENOM),MATZ,WORK(KSOLEQ),IERR) 2481 DO I = 1, NREDH 2482 IF (ABS(WORK(KDENOM-1+I)).GT.THRZER) THEN 2483 WORK(KEIVAL-1+I) = WORK(KEIVAL-1+I)/WORK(KDENOM-1+I) 2484 WORK(KWI-1+I) = WORK(KWI-1+I) /WORK(KDENOM-1+I) 2485 ELSE 2486 WORK(KEIVAL-1+I) = 1.0D0/THRZER 2487 WORK(KWI-1+I) = WORK(KWI-1+I)/THRZER 2488 END IF 2489 END DO 2490 ELSE 2491 CALL DCOPY(MAXRED*MAXRED,WORK(KSOLEQ),1,WORK(KAMAT),1) 2492 CALL DCOPY(MAXRED*MAXRED,WORK(KSOLEQ),1,WORK(KAMAT2),1) 2493 CALL RG(MAXRED,NREDH,WORK(KAMAT),WORK(KEIVAL),WORK(KWI),MATZ, 2494 * WORK(KSOLEQ),WORK(KIV1),WORK(KFV1),IERR) 2495 END IF 2496C 2497 IF (IPRINT .GE. 70 .OR. IERR .NE. 0) THEN 2498 WRITE (LUPRI,'(/A)') ' EIGENVALUES real part:' 2499 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2500 CALL OUTPUT(WORK(KEIVAL),1,NREDH,1,1,NREDH,MAXRED,1,LUPRI) 2501 WRITE (LUPRI,'(/A)') 2502 * ' EIGENVALUES imaginary part:' 2503 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2504 CALL OUTPUT(WORK(KWI),1,NREDH,1,1,NREDH,MAXRED,1,LUPRI) 2505 WRITE (LUPRI,'(/A)') ' EIGENVECTORS :' 2506 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2507 CALL OUTPUT(WORK(KSOLEQ),1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2508 END IF 2509C 2510 CALL RGORD(MAXRED,NREDH,WORK(KEIVAL),WORK(KWI),WORK(KSOLEQ), 2511 * .FALSE.) 2512C 2513 IF (IPRINT .GE. 70 ) THEN 2514 WRITE (LUPRI,'(/A)') ' EIGENVECTORS :' 2515 WRITE (LUPRI,'(A)') ' After sort of eigenvalues' 2516 CALL OUTPUT(WORK(KSOLEQ),1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2517 END IF 2518 IF ( IERR.NE.0 ) THEN 2519 WRITE(LUPRI,'(/A,I5)') 2520 * ' EIGENVALUE PROBLEM NOT CONVERGED IN RG, IERR =',IERR 2521 CALL QUIT(' CCLR_EXPJAC: EIGENVALUE EQUATION NOT CONVERGED ') 2522 END IF 2523C 2524c CALL RGTST(LUPRI,IPRINT,MAXRED,NREDH,WORK(KAMAT2),WORK(KSOLEQ), 2525c * WORK(KWRK1),LWRK1,INFO) 2526C 2527C---------------------------------- 2528C Write out Excitation energies. 2529C---------------------------------- 2530C 2531 WRITE (LUPRI,'(//A/A/A//A/A/)') 2532 *' CCSD excitation energies :', 2533 *' ============================', 2534 *' (conversion factor used: 1 au = 27.2113957 eV)', 2535 *' Excitation no. Hartree eV', 2536 *' -------------- ------- --' 2537 DO 400 I = 1,NREDH 2538 WRITE (LUPRI,'(I10,2F20.6)') 2539 * I,WORK(KEIVAL-1+I),WORK(KEIVAL-1+I)*XTEV 2540 400 CONTINUE 2541C 2542C--------------------------------- 2543C Analysis of solution vectors. 2544C--------------------------------- 2545C 2546 THRESH = 0.05 2547 MAXLIN = 100 2548C 2549 DO 500 I = 1,NREDH 2550 WRITE(LUPRI,'(//5X,A,I2)') 2551 *'Analysis of the Coupled Cluster Excitation Vector Number : ',I 2552 WRITE(LUPRI,'(5X,A)') 2553 *'-------------------------------------------------------------' 2554 WRITE(LUPRI,'(/15X,A,F10.5,A)') 2555 * 'Excitation Energy : ',WORK(KEIVAL+I-1)*XTEV, 2556 * ' eV' 2557 CALL CC_PRAM(WORK(KSOLEQ + (I-1)*MAXRED),PT1,1,.FALSE.) 2558 500 CONTINUE 2559C 2560 CALL QEXIT('FDEXCI') 2561c 2562 RETURN 2563 END 2564c*DECK CCLR_JACIN 2565 SUBROUTINE CCLR_JACIN(FDJACO,SMATRIX) 2566C 2567C------------------------------------------------------------------- 2568C Written by Ove Christiansen 21-11-1994 2569C 2570C Reads the jacobian from disk and put it into JAC. 2571C 2572C------------------------------------------------------------------- 2573C 2574#include "implicit.h" 2575#include "priunit.h" 2576#include "ccorb.h" 2577#include "maxorb.h" 2578#include "ccsdsym.h" 2579#include "ccsdinp.h" 2580#include "cclr.h" 2581#include "ccsdio.h" 2582#include "aovec.h" 2583#include "cbirea.h" 2584#include "r12int.h" 2585C 2586 DIMENSION FDJACO(MAXRED,MAXRED) 2587 DIMENSION SMATRIX(MAXRED,MAXRED) 2588C 2589 NTAMP = NT1AMX + NT2AMX 2590 IF (LMULBS) NTAMP = NTAMP + NTR12AM(1) 2591C 2592C------------------------ 2593C Open Jacobian file. 2594C------------------------ 2595C 2596 LUJAC = -1 2597 CALL GPOPEN(LUJAC,'CCLR_JAC','UNKNOWN',' ','UNFORMATTED',IDUMMY, 2598 & .FALSE.) 2599 REWIND(LUJAC) 2600C 2601 DO I = 1, NTAMP 2602 READ(LUJAC) (FDJACO(J,I),J=1,NTAMP) 2603 END DO 2604C 2605 CALL GPCLOSE(LUJAC,'KEEP') 2606 2607 2608 IF (CCR12) THEN 2609 LUMET = -1 2610 CALL GPOPEN(LUMET,'CCLR_MET','UNKNOWN',' ','UNFORMATTED', 2611 * IDUMMY,.FALSE.) 2612 CALL DZERO(SMATRIX,MAXRED*MAXRED) 2613 DO I = 1, NT1AMX+NT2AMX 2614 SMATRIX(I,I) = 1.0D0 2615 END DO 2616chf 2617 IF (IANR12.EQ.2) THEN 2618 LUMET2 = -1 2619 CALL GPOPEN(LUMET2,'CCLR_MET2','UNKNOWN',' ','UNFORMATTED', 2620 * IDUMMY,.FALSE.) 2621 LUMET3 = -1 2622 CALL GPOPEN(LUMET3,'CCLR_MET3','UNKNOWN',' ','UNFORMATTED', 2623 * IDUMMY,.FALSE.) 2624 2625 DO I = NT1AM(1)+NT2AM(1)+1, NT1AM(1)+NT2AM(1)+NTR12AM(1) 2626 READ(LUMET2)(SMATRIX(NT1AM(1)+J,I),J=1,NT2AM(1)) 2627 END DO 2628c 2629 DO I = NT1AM(1)+1, NT1AM(1)+NT2AM(1) 2630 READ(LUMET3)(SMATRIX(NT1AM(1)+NT2AM(1)+J,I),J=1,NTR12AM(1)) 2631 END DO 2632 CALL GPCLOSE(LUMET2,'KEEP') 2633 CALL GPCLOSE(LUMET3,'KEEP') 2634 END IF 2635c 2636 DO I = NT1AMX+NT2AMX+1, NT1AMX+NT2AMX+NTR12AM(1) 2637 READ(LUMET) (SMATRIX(NT1AMX+NT2AMX+J,I),J=1,NTR12AM(1)) 2638 END DO 2639 CALL GPCLOSE(LUMET,'KEEP') 2640c WRITE(LUPRI,*)'S-Matrix out of cclr_jacin' 2641c CALL OUTPUT(SMATRIX,1,NT1AMX+NT2AMX+NTR12AM(1),1, 2642c & NT1AMX+NT2AMX+NTR12AM(1),MAXRED,MAXRED,1,LUPRI) 2643 END IF 2644C 2645 RETURN 2646C 2647 END 2648C /* Deck cc_core */ 2649 SUBROUTINE CC_CORE(T1AM,T2AM,ISYMTR) 2650C 2651C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2652C 30-5-1995 Ove Christiansen 2653C 2654C Purpose: To zero core and secondary part of 2655C CC vector. 2656C ISYMTR is symmetry of vector. 2657C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2658C 2659#include "implicit.h" 2660#include "ccorb.h" 2661#include "ccsdsym.h" 2662#include "ccsdinp.h" 2663C 2664 PARAMETER ( ZERO= 0.0D0) 2665C 2666 DIMENSION T1AM(*),T2AM(*) 2667C 2668 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2669C 2670 IF (IPRINT .GT. 45) THEN 2671 CALL AROUND('START OF CC_CORE: (T1,T2) vector packed ') 2672 CALL CC_PRP(T1AM,T2AM,ISYMTR,1,1) 2673 ENDIF 2674C 2675 ISYMAI = MULD2H(ISYMTR,ISYMOP) 2676C 2677C------------------------------------------- 2678C Loop through single excitation vector. 2679C------------------------------------------- 2680C 2681 DO 100 ISYMA = 1,NSYM 2682C 2683 ISYMI = MULD2H(ISYMAI,ISYMA) 2684C 2685 IF (LCOR) THEN 2686C 2687 DO 110 I = 1,NRHF(ISYMI) 2688C 2689 IF (I .LE. ICOR(ISYMI)) THEN 2690 KOFF = IT1AM(ISYMA,ISYMI) 2691 * + NVIR(ISYMA)*(I-1) + 1 2692 CALL DZERO(T1AM(KOFF),NVIR(ISYMA)) 2693 ENDIF 2694C 2695 110 CONTINUE 2696C 2697 ENDIF 2698 IF (LSEC) THEN 2699C 2700 DO 120 A=1,NVIR(ISYMA) 2701C 2702 IA = NVIR(ISYMA) - A + 1 2703C 2704 IF (IA .LE. ISEC(ISYMA)) THEN 2705 DO 130 I = 1, NRHF(ISYMI) 2706 NAI = IT1AM(ISYMA,ISYMI) 2707 * + NVIR(ISYMA)*(I-1) + A 2708 T1AM(NAI) = ZERO 2709 130 CONTINUE 2710 ENDIF 2711C 2712 120 CONTINUE 2713C 2714 ENDIF 2715C 2716 100 CONTINUE 2717C 2718C-------------------------------------------- 2719C Loop through Doublee excitation vector. 2720C If not ccs or ccp2 2721C-------------------------------------------- 2722C 2723 IF ( CCS ) RETURN 2724C 2725 ISAIBJ = MULD2H(ISYMTR,ISYMOP) 2726C 2727 DO 200 ISYMAI = 1,NSYM 2728C 2729 ISYMBJ = MULD2H(ISYMAI,ISAIBJ) 2730C 2731 DO 210 ISYMJ = 1,NSYM 2732C 2733 ISYMB = MULD2H(ISYMJ,ISYMBJ) 2734C 2735 DO 220 ISYMI = 1,NSYM 2736C 2737 ISYMA = MULD2H(ISYMI,ISYMAI) 2738C 2739 DO 230 J = 1,NRHF(ISYMJ) 2740C 2741 DO 240 B = 1,NVIR(ISYMB) 2742C 2743 NBJ = IT1AM(ISYMB,ISYMJ) 2744 * + NVIR(ISYMB)*(J - 1) + B 2745C 2746 DO 250 I = 1,NRHF(ISYMI) 2747C 2748 DO 260 A = 1,NVIR(ISYMA) 2749C 2750 NAI = IT1AM(ISYMA,ISYMI) 2751 * + NVIR(ISYMA)*(I - 1) + A 2752C 2753 IF (((ISYMAI.EQ.ISYMBJ).AND. 2754 * (NAI .GT. NBJ)).OR.(ISYMAI.GT.ISYMBJ)) 2755 * GOTO 260 2756C 2757 IF (ISYMAI.EQ.ISYMBJ) THEN 2758 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2759 * + INDEX(NAI,NBJ) 2760 ELSE 2761 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 2762 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 2763 ENDIF 2764C 2765 IF (LCOR) THEN 2766C 2767 IF (I .LE. ICOR(ISYMI)) THEN 2768C 2769 T2AM(NAIBJ) = ZERO 2770C 2771 ELSE IF (J .LE. ICOR(ISYMJ)) THEN 2772C 2773 T2AM(NAIBJ) = ZERO 2774C 2775 ENDIF 2776C 2777 ENDIF 2778C 2779 IF (LSEC) THEN 2780C 2781 IA = NVIR(ISYMA) - A + 1 2782 IB = NVIR(ISYMB) - B + 1 2783 IF (IA .LE. ISEC(ISYMA)) THEN 2784C 2785 T2AM(NAIBJ) = ZERO 2786C 2787 ELSE IF (IB .LE. ISEC(ISYMB)) THEN 2788C 2789 T2AM(NAIBJ) = ZERO 2790C 2791 ENDIF 2792C 2793 ENDIF 2794C 2795 260 CONTINUE 2796 250 CONTINUE 2797 240 CONTINUE 2798 230 CONTINUE 2799 220 CONTINUE 2800 210 CONTINUE 2801 200 CONTINUE 2802C 2803 IF (IPRINT .GT. 45) THEN 2804 CALL AROUND('END OF CC_CORE: (T1,T2) vector packed ') 2805 CALL CC_PRP(T1AM,T2AM,ISYMTR,1,1) 2806 ENDIF 2807C 2808 RETURN 2809 END 2810C /* Deck cc_pram */ 2811 SUBROUTINE CC_PRAM(CAM,PT1,ISYMTR,LGRS) 2812C 2813C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2814C 30-5-1995 Ove Christiansen 2815C 2816C Purpose: Writes out vector: 2817C %T1 and %T2 and ||T1||/||T2|| 2818C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 2819C 2820#include "implicit.h" 2821C 2822 PARAMETER (TWO = 2.0D0, THPRT = 1.0D-9, HUNDRED = 100.0D0) 2823C 2824#include "ccorb.h" 2825#include "ccsdsym.h" 2826#include "ccsdinp.h" 2827#include "priunit.h" 2828Cholesky 2829#include "maxorb.h" 2830#include "ccdeco.h" 2831Cmlcc 2832#include "peract.h" 2833#include "mxcent.h" 2834#include "center.h" 2835Cmlcc 2836C 2837 LOGICAL CCSEFF, first1, first2 2838Cholesky 2839C 2840C 2841 LOGICAL LGRS 2842 DIMENSION CAM(*) 2843C 2844 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 2845C 2846Cholesky 2847 CCSEFF = CCS .OR. (CHOINT.AND.CC2) 2848Cholesky 2849C 2850C------------------------ 2851C Add up the vectors. 2852C------------------------ 2853C 2854 C1NOSQ = 0.0D0 2855 C2NOSQ = 0.0D0 2856 KC1 = 1 2857 KC2 = KC1 + NT1AM(ISYMTR) 2858 C1NOSQ = DDOT(NT1AM(ISYMTR),CAM(KC1),1,CAM(KC1),1) 2859Chol IF (.NOT. CCS) C2NOSQ = DDOT(NT2AM(ISYMTR),CAM(KC2),1,CAM(KC2),1) 2860 IF (.NOT. CCSEFF) 2861 & C2NOSQ = DDOT(NT2AM(ISYMTR),CAM(KC2),1,CAM(KC2),1) 2862 CNOSQ = C1NOSQ + C2NOSQ 2863C 2864 IF (.NOT. (CCSEFF.OR.CCP2) .AND. CNOSQ.NE.0.0D0) THEN 2865C 2866 WRITE(LUPRI,'(//10X,A)') 2867 * 'CC_PRAM:Overall Contribution of the Different Components' 2868 WRITE(LUPRI,'(10X,A//)') 2869 * '--------------------------------------------------------' 2870 WRITE(LUPRI,'(/10X,A,10X,F10.4,A)') 2871 * 'Single Excitation Contribution : ', 2872 * (C1NOSQ/CNOSQ)*HUNDRED,' %' 2873 WRITE(LUPRI,'(/10X,A,10X,F10.4,A)') 2874 * 'Double Excitation Contribution : ', 2875 * (C2NOSQ/CNOSQ)*HUNDRED,' %' 2876 WRITE(LUPRI,'(/10X,A,10X,F10.4)') 2877 * '||T1||/||T2|| : ', 2878 * SQRT(C1NOSQ/C2NOSQ) 2879 IF (LGRS) WRITE(LUPRI,'(/10X,A,10X,F10.4)') 2880 * 'Tau1 diagnostic : ', 2881 * SQRT(C1NOSQ/(TWO*DFLOAT(NRHFT))) 2882 PT1 = (C1NOSQ/CNOSQ)*HUNDRED 2883 ELSE 2884 PT1 = HUNDRED 2885 ENDIF 2886 WRITE(LUPRI,'(/10X,A,10X,F10.4)') 2887 * 'Norm of Total Amplitude Vector : ',SQRT(CNOSQ) 2888C 2889 CALL FLSHFO(LUPRI) 2890C 2891C---------------------------------------------- 2892C Initialize threshold etc from Printlevel. 2893C---------------------------------------------- 2894C 2895 NL = MAX(1,2*IPRINT) 2896C 2897 CNOSQ = MAX(CNOSQ,THPRT) 2898C 2899 THR1 = SQRT(C1NOSQ/CNOSQ)/NL 2900 THR2 = SQRT(C2NOSQ/CNOSQ)/NL 2901 THR1 = MAX(THR1,1.0D-02) 2902 THR2 = MAX(THR2,1.0D-02) 2903 SUMOFP = 0.0D00 2904 first1 = .true. 2905 first2 = .true. 2906 sumoftt = 0.0D00 2907 sumofts = 0.0D00 2908 sumofst = 0.0D00 2909 sumofss = 0.0D00 2910 sumofrr = 0.0D00 2911 sumsing = 0.0D00 2912 sum2int = 0.0D00 2913 sum2se = 0.0D00 2914 sum2ext = 0.0D00 2915C 2916 IF (DEBUG) THR1 = 0.0D0 2917C 2918C------------------------------------- 2919C Loop until a few is Printed out. 2920C------------------------------------- 2921C 2922C 2923C--------------------------------------- 2924C Loop through One excitation part. 2925C--------------------------------------- 2926C 2927 WRITE(LUPRI,'(//A)') 2928 * ' +==============================================' 2929 * //'===============================+' 2930 WRITE(LUPRI,'(1X,A)') 2931 * '| symmetry| orbital index | Excitation Numbers' 2932 * //' | Amplitude |' 2933 WRITE(LUPRI,'(1X,A)') 2934 * '| Index | a b i j | NAI NBJ |' 2935 * //' NAIBJ | |' 2936 WRITE(LUPRI,'(A)') 2937 * ' +==============================================' 2938 * //'===============================+' 2939C 2940 ISYMAI = MULD2H(ISYMTR,ISYMOP) 2941C 2942 1 CONTINUE 2943 N1 = 0 2944C 2945 DO 100 ISYMA = 1,NSYM 2946C 2947 ISYMI = MULD2H(ISYMAI,ISYMA) 2948C 2949 DO 110 I = 1,NRHF(ISYMI) 2950C 2951 MI = IORB(ISYMI) + I 2952C 2953 DO 120 A=1,NVIR(ISYMA) 2954C 2955 NAI = IT1AM(ISYMA,ISYMI) + NVIR(ISYMA)*(I-1) + A 2956C 2957 MA = IORB(ISYMA) + NRHF(ISYMA) + A 2958C 2959 IF (ABS(CAM(NAI)) .GE. THR1 ) THEN 2960C 2961 WRITE(LUPRI,9990) ISYMA,ISYMI,A,I,NAI,CAM(NAI) 2962C 2963 N1 = N1 + 1 2964 SUMOFP = SUMOFP + CAM(NAI)*CAM(NAI) 2965C 2966 ENDIF 2967C 2968 if (first1) then 2969 if ((mlcc3 .and. pertcc) 2970 * .and. nspace .eq. 1) then 2971 if (actorb(i,isymi) 2972 * .and. actorb(a+nrhf(isyma),isyma)) then 2973C 2974 sumoftt = sumoftt + cam(nai)*cam(nai) 2975C 2976 else if (actorb(i,isymi) .and. 2977 * .not. actorb(a+nrhf(isyma),isyma)) then 2978C 2979 sumofts = sumofts + cam(nai)*cam(nai) 2980C 2981 else if (.not. actorb(i,isymi) 2982 * .and. actorb(a+nrhf(isyma),isyma)) then 2983C 2984 sumofst = sumofst + cam(nai)*cam(nai) 2985C 2986 else if (.not. actorb(i,isymi) 2987 * .and. .not. actorb(a+nrhf(isyma),isyma)) then 2988C 2989 sumofss = sumofss + cam(nai)*cam(nai) 2990C 2991 endif 2992 endif 2993C 2994C 2995 if ((mlcc3 .and. pertcc) 2996 * .and. nspace .eq. 2) then 2997 if (iacorb(i,isymi) .eq. -5 2998 * .or. iacorb(a+nrhf(isyma),isyma) .eq. -5) then 2999 3000 sumofrr = sumofrr + cam(nai)*cam(nai) 3001 3002 else if (iacorb(i,isymi) .eq. 1 3003 * .and. iacorb(a+nrhf(isyma),isyma) .eq. 1) then 3004 3005 sumoftt = sumoftt + cam(nai)*cam(nai) 3006 3007 else if (iacorb(i,isymi) .eq. 1 .and. 3008 * .not. iacorb(a+nrhf(isyma),isyma) .eq. 1) then 3009 3010 sumofts = sumofts + cam(nai)*cam(nai) 3011 3012 else if (.not. iacorb(i,isymi) .eq. 1 3013 * .and. iacorb(a+nrhf(isyma),isyma) .eq. 1) then 3014 3015 sumofst = sumofst + cam(nai)*cam(nai) 3016 3017 else if (.not. iacorb(i,isymi) .eq. 1 .and. 3018 * .not. iacorb(a+nrhf(isyma),isyma) .eq. 1) then 3019 3020 sumofss = sumofss + cam(nai)*cam(nai) 3021 3022 endif 3023 endif 3024 endif 3025C 3026C 3027 120 CONTINUE 3028 110 CONTINUE 3029 100 CONTINUE 3030C 3031 sumsing = sumoftt + sumofst + sumofts + sumofss 3032C 3033 IF ((N1.LT.1).AND.(SQRT(C1NOSQ/CNOSQ).GT.1.0D-3)) THEN 3034 THR1 = THR1/5.0D0 3035 first1 = .false. 3036 GOTO 1 3037 ENDIF 3038C 3039 CALL FLSHFO(LUPRI) 3040C 3041C-------------------------------------------- 3042C Loop through Double excitation vector. 3043C If not ccs or ccp2 3044C-------------------------------------------- 3045C 3046 IF (.NOT. ( CCSEFF .OR. CCP2 )) THEN 3047C 3048 WRITE(LUPRI,'(A)') 3049 * ' +----------------------------------------------' 3050 * //'-------------------------------+' 3051C 3052 3053 2 CONTINUE 3054 N2 = 0 3055C 3056 DO 200 ISYMAI = 1,NSYM 3057C 3058 ISYMBJ = MULD2H(ISYMAI,ISYMTR) 3059C 3060 DO 210 ISYMJ = 1,NSYM 3061C 3062 ISYMB = MULD2H(ISYMJ,ISYMBJ) 3063C 3064 DO 220 ISYMI = 1,NSYM 3065C 3066 ISYMA = MULD2H(ISYMI,ISYMAI) 3067C 3068 DO 230 J = 1,NRHF(ISYMJ) 3069C 3070 MJ = IORB(ISYMJ) + J 3071C 3072 DO 240 B = 1,NVIR(ISYMB) 3073C 3074 NBJ = IT1AM(ISYMB,ISYMJ) 3075 * + NVIR(ISYMB)*(J - 1) + B 3076C 3077 MB = IORB(ISYMB) + NRHF(ISYMB) + B 3078C 3079 DO 250 I = 1,NRHF(ISYMI) 3080C 3081 MI = IORB(ISYMI) + I 3082C 3083 DO 260 A = 1,NVIR(ISYMA) 3084C 3085 NAI = IT1AM(ISYMA,ISYMI) 3086 * + NVIR(ISYMA)*(I - 1) + A 3087C 3088 MA = IORB(ISYMA) + NRHF(ISYMA) + A 3089C 3090 IF (((ISYMAI.EQ.ISYMBJ).AND. 3091 * (NAI .LT. NBJ)).OR.(ISYMAI.LT.ISYMBJ)) 3092 * GOTO 260 3093C 3094 IF (ISYMAI.EQ.ISYMBJ) THEN 3095 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 3096 * + INDEX(NAI,NBJ) 3097 ELSE 3098 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 3099 * + NT1AM(ISYMAI)*(NBJ-1) + NAI 3100 ENDIF 3101C 3102 KAIBJ = NAIBJ + NT1AM(ISYMTR) 3103 IF (ABS(CAM(KAIBJ)) .GT. THR2 ) THEN 3104C 3105 WRITE(LUPRI,9991) ISYMA,ISYMB,ISYMI,ISYMJ, 3106 * A,B,I,J,NAI,NBJ,NAIBJ, 3107 * CAM(KAIBJ) 3108 N2 = N2 + 1 3109C 3110 SUMOFP = SUMOFP + CAM(KAIBJ)*CAM(KAIBJ) 3111C 3112 ENDIF 3113C 3114 if (first2) then 3115 if ((mlcc3 .and. pertcc) 3116 * .and. nspace .eq. 1) then 3117 if(actorb(i,isymi) .and. actorb(j,isymj) 3118 * .and. actorb(a+nrhf(isyma),isyma) 3119 * .and. actorb(b+nrhf(isymb),isymb)) then 3120 3121 sum2int = sum2int+cam(kaibj)*cam(kaibj) 3122 3123 else if (.not. actorb(i,isymi) 3124 * .and. .not. actorb(j,isymj) 3125 * .and. .not. actorb(a+nrhf(isyma),isyma) 3126 * .and. .not. actorb(b+nrhf(isymb),isymb)) 3127 * then 3128 3129 sum2ext = sum2ext+cam(kaibj)*cam(kaibj) 3130 3131 else 3132 3133 sum2se = sum2se+cam(kaibj)*cam(kaibj) 3134 3135 endif 3136 endif 3137C 3138C 3139 if ((mlcc3 .and. pertcc) 3140 * .and. nspace .eq. 2) then 3141C 3142 itest = iacorb(i,isymi) + iacorb(j,isymj) 3143 * + iacorb(a+nrhf(isyma),isyma) 3144 * + iacorb(b+nrhf(isymb),isymb) 3145C 3146 if(itest .eq. 4) then 3147 3148 sum2int = sum2int+cam(kaibj)*cam(kaibj) 3149 3150 else if(itest .eq. 0) then 3151 3152 sum2ext = sum2ext+cam(kaibj)*cam(kaibj) 3153 3154 else 3155 3156 sum2se = sum2se+cam(kaibj)*cam(kaibj) 3157 3158 endif 3159C 3160 endif 3161 endif 3162C 3163C 3164 260 CONTINUE 3165 250 CONTINUE 3166 240 CONTINUE 3167 230 CONTINUE 3168 220 CONTINUE 3169 210 CONTINUE 3170 200 CONTINUE 3171C 3172 IF ((N2.LT.1).AND.(SQRT(C2NOSQ/CNOSQ).GT.1.0D-3)) THEN 3173 THR2 = THR2/5D00 3174 first2 = .false. 3175 GOTO 2 3176 ENDIF 3177C 3178 ENDIF 3179C 3180 WRITE(LUPRI,'(A)') 3181 * ' +==============================================' 3182 * //'===============================+' 3183C 3184 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3185 * 'Norm of Printed Amplitude Vector : ',SQRT(SUMOFP) 3186 WRITE(LUPRI,'(/10X,A43,F9.6)') 3187 * 'Printed all single excitations greater than',THR1 3188 IF (.NOT. (CCSEFF.OR.CCP2)) THEN 3189 WRITE(LUPRI,'(/10X,A43,F9.6)') 3190 * 'Printed all double excitations greater than',THR2 3191 ENDIF 3192C 3193 if ((mlcc3 .and. pertcc) .and. nspace .le. 2) then 3194 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3195 * 'Norm of T to T single amplitudes : ',sumoftt 3196 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3197 * 'Norm of T to S single amplitudes : ',sumofts 3198 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3199 * 'Norm of S to T single amplitudes : ',sumofst 3200 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3201 * 'Norm of S to S single amplitudes : ',sumofss 3202 if (nspace .eq. 2) then 3203 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3204 * 'Norm of R single amplitudes : ',sumofrr 3205 endif 3206 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3207 * 'Norm of T to T double amplitudes : ',sum2int 3208 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3209 * 'Norm of X to Y double amplitudes : ',sum2se 3210 WRITE(LUPRI,'(//10X,A,8X,F10.4)') 3211 * 'Norm of S to S double amplitudes : ',sum2ext 3212 endif 3213C 3214C 3215 9990 FORMAT(1X,'| ',I1,3X,I1,2X,' | ',I3,5X,I3,4X,' | ',I8,9x, 3216 * ' | ',12x,' | ',1x, F10.6,' |') 3217 9991 FORMAT(1X,'| ',I1,1X,I1,1X,I1,1X,I1,' | ', 3218 * I3,1X,I3,1X,I3,1X,I3,' | ', 3219 * I8,1x,I8,' | ',I12,' | ',1x,F10.6,' |') 3220C 3221 RETURN 3222 END 3223C /* Deck cc_prtm */ 3224 SUBROUTINE CC_PRTM(TM,ISYMD,ISYMTR) 3225C 3226#include "implicit.h" 3227#include "iratdef.h" 3228 DIMENSION TM(*) 3229#include "ccorb.h" 3230#include "ccsdsym.h" 3231#include "priunit.h" 3232C 3233 DO 100 ISYMJ = 1,NSYM 3234C 3235 ISYMDJ = MULD2H(ISYMJ,ISYMD) 3236 ISYMCI = MULD2H(ISYMTR,ISYMDJ) 3237C 3238 NTOTCI = MAX(NT1AM(ISYMCI),1) 3239C 3240 WRITE(LUPRI,*) 'Transformede double exc. ampitudes Tm ' 3241 WRITE(LUPRI,*) ' ' 3242 WRITE(LUPRI,*) 'ISYMCI,ISYMJ = ',ISYMCI,ISYMJ 3243 DO 110 J = 1,NRHF(ISYMJ) 3244C 3245 KOFF3 = IT2BCD(ISYMCI,ISYMJ) 3246 * + NT1AM(ISYMCI)*(J-1) + 1 3247C 3248 WRITE(LUPRI,*) 'FOR J= ',J 3249 CALL OUTPUT(TM(KOFF3),1,NT1AM(ISYMCI),1,1, 3250 * NT1AM(ISYMCI),1,1,LUPRI) 3251C 3252 110 CONTINUE 3253 100 CONTINUE 3254C 3255 END 3256C /* Deck cc_prlam */ 3257 SUBROUTINE CC_PRLAM(XLAMDP,XLAMDH,ISYML) 3258C 3259C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3260C 3261C Ove Christiansen 13-7-1995 based on lammat by Henrik Koch. 3262C 3263C PURPOSE: 3264C Printout lambda matrices 3265C for usual lambda in CC opt ISYML = 1 3266C for c1 transformed ISYML = ISYMTR 3267C 3268C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3269C 3270#include "implicit.h" 3271 DIMENSION XLAMDH(*),XLAMDP(*) 3272#include "ccorb.h" 3273#include "ccsdsym.h" 3274#include "priunit.h" 3275C 3276 WRITE(LUPRI,*) 3277 WRITE(LUPRI,*) 'In cc_prlam, symmetry of lambda matrices is :' 3278 * ,ISYML 3279 CALL AROUND('Lambda Particle matrix ') 3280 DO 200 ISYMAO = 1,NSYM 3281 ISYMMO = MULD2H(ISYML,ISYMAO) 3282 WRITE(LUPRI,1) ISYMAO,ISYMMO 3283 WRITE(LUPRI,2) 3284 WRITE(LUPRI,3) 3285 IF (NRHF(ISYMMO) .EQ. 0) THEN 3286 WRITE(LUPRI,4) 3287 GOTO 210 3288 ENDIF 3289 KOFF1 = IGLMRH(ISYMAO,ISYMMO) + 1 3290 CALL OUTPUT(XLAMDP(KOFF1),1,NBAS(ISYMAO),1,NRHF(ISYMMO), 3291 * NBAS(ISYMAO),NRHF(ISYMMO),1,LUPRI) 3292 210 WRITE(LUPRI,5) 3293 WRITE(LUPRI,6) 3294 IF (NVIR(ISYMMO) .EQ. 0) THEN 3295 WRITE(LUPRI,4) 3296 GOTO 220 3297 ENDIF 3298 KOFF2 = IGLMVI(ISYMAO,ISYMMO) + 1 3299 CALL OUTPUT(XLAMDP(KOFF2),1,NBAS(ISYMAO),1,NVIR(ISYMMO), 3300 * NBAS(ISYMAO),NVIR(ISYMMO),1,LUPRI) 3301C 3302 220 CONTINUE 3303 200 CONTINUE 3304C 3305 CALL AROUND('Lambda Hole matrix ') 3306 DO 300 ISYMAO = 1,NSYM 3307 ISYMMO = MULD2H(ISYML,ISYMAO) 3308 WRITE(LUPRI,1) ISYMAO,ISYMMO 3309 WRITE(LUPRI,7) 3310 WRITE(LUPRI,8) 3311 IF (NRHF(ISYMMO) .EQ. 0) THEN 3312 WRITE(LUPRI,4) 3313 GOTO 310 3314 ENDIF 3315 KOFF1 = IGLMRH(ISYMAO,ISYMMO) + 1 3316 CALL OUTPUT(XLAMDH(KOFF1),1,NBAS(ISYMAO),1,NRHF(ISYMMO), 3317 * NBAS(ISYMAO),NRHF(ISYMMO),1,LUPRI) 3318 310 WRITE(LUPRI,9) 3319 WRITE(LUPRI,10) 3320 IF (NVIR(ISYMMO) .EQ. 0) THEN 3321 WRITE(LUPRI,4) 3322 GOTO 320 3323 ENDIF 3324 KOFF2 = IGLMVI(ISYMAO,ISYMMO) + 1 3325 CALL OUTPUT(XLAMDH(KOFF2),1,NBAS(ISYMAO),1,NVIR(ISYMMO), 3326 * NBAS(ISYMAO),NVIR(ISYMMO),1,LUPRI) 3327C 3328 320 CONTINUE 3329 300 CONTINUE 3330C 3331 RETURN 3332C 3333 1 FORMAT(/,/,7X,'Symmetry numbers ao,mo:',I5,I5) 3334 2 FORMAT(/,/,7X,'Lambda particle occupied part') 3335 3 FORMAT(7X,'-----------------------------') 3336 4 FORMAT(/,/,7X,'This symmetry is empty') 3337 5 FORMAT(/,/,7X,'Lambda particle virtual part') 3338 6 FORMAT(7X,'----------------------------') 3339 7 FORMAT(/,/,7X,'Lambda hole occupied part') 3340 8 FORMAT(7X,'-------------------------') 3341 9 FORMAT(/,/,7X,'Lambda hole virtual part') 3342 10 FORMAT(7X,'------------------------') 3343C 3344 END 3345C /* Deck cc_print */ 3346 SUBROUTINE CC_PRINT(XINT,DSRHF,ISYDIS) 3347C 3348C Written by Ove Christiansen 24-7-1995 3349C 3350C Purpose: Write out integrals. 3351C 3352#include "implicit.h" 3353C 3354 DIMENSION XINT(*),DSRHF(*) 3355C 3356#include "ccorb.h" 3357#include "ccsdsym.h" 3358#include "priunit.h" 3359C 3360 KOFF1 = 1 3361 KOFF2 = 1 3362 DO 100 ISYMJ = 1,NSYM 3363C 3364 ISYMG = ISYMJ 3365 ISYMAB = MULD2H(ISYMG,ISYDIS) 3366C 3367 NNBSAB = MAX(NNBST(ISYMAB),1) 3368 NBASG = MAX(NBAS(ISYMG),1) 3369C 3370 CALL AROUND( ' AO - INTEGRALS ') 3371 WRITE(LUPRI,*) 'ISYMG: ',ISYMG 3372 CALL OUTPUT(XINT(KOFF1),1,NNBST(ISYMAB),1,NBAS(ISYMG), 3373 * NNBSAB,NBASG,1,LUPRI) 3374C 3375 CALL AROUND( 'DSRHF ') 3376 WRITE(LUPRI,*) 'ISYMJ: ',ISYMJ 3377 CALL OUTPUT(DSRHF(KOFF2),1,NNBST(ISYMAB),1,NRHF(ISYMJ), 3378 * NNBSAB,NBASG,1,LUPRI) 3379C 3380 KOFF1 = KOFF1 + NNBST(ISYMAB)*NBAS(ISYMG) 3381 KOFF2 = KOFF2 + NNBST(ISYMAB)*NRHF(ISYMJ) 3382C 3383 100 CONTINUE 3384C 3385 RETURN 3386 END 3387C /* Deck cc_prei */ 3388 SUBROUTINE CC_PREI(EI1,EI2,ISYMEI,LEI1MO) 3389C 3390#include "implicit.h" 3391#include "iratdef.h" 3392 DIMENSION EI1(*),EI2(*) 3393#include "ccorb.h" 3394#include "ccsdsym.h" 3395#include "priunit.h" 3396C 3397 CALL AROUND( 'EI1 -intermediate ') 3398C 3399 DO 100 ISYMB = 1,NSYM 3400C 3401 ISYMA = MULD2H(ISYMB,ISYMEI) 3402C 3403 WRITE(LUPRI,*) ' ' 3404 WRITE(LUPRI,*) 'ISYMA,ISYMB = ',ISYMA,ISYMB 3405 WRITE(LUPRI,*) ' ' 3406C 3407 IF (LEI1MO.EQ.1) THEN 3408 KOFF1 = IMATAB(ISYMA,ISYMB) + 1 3409 LB = NVIR(ISYMB) 3410 ELSE 3411 KOFF1 = IEMAT1(ISYMA,ISYMB) + 1 3412 LB = NBAS(ISYMB) 3413 ENDIF 3414 CALL OUTPUT(EI1(KOFF1),1,NVIR(ISYMA),1,LB, 3415 * NVIR(ISYMA),LB,1,LUPRI) 3416C 3417 100 CONTINUE 3418C 3419 CALL AROUND( 'EI2 -intermediate ') 3420C 3421 DO 200 ISYMJ = 1,NSYM 3422C 3423 ISYMI = MULD2H(ISYMJ,ISYMEI) 3424C 3425 WRITE(LUPRI,*) ' ' 3426 WRITE(LUPRI,*) 'ISYMI,ISYMJ = ',ISYMI,ISYMJ 3427 WRITE(LUPRI,*) ' ' 3428C 3429 KOFF1 = 1 + IMATIJ(ISYMI,ISYMJ) 3430 CALL OUTPUT(EI2(KOFF1),1,NRHF(ISYMI),1,NRHF(ISYMJ), 3431 * NRHF(ISYMI),NRHF(ISYMJ),1,LUPRI) 3432C 3433 200 CONTINUE 3434C 3435 END 3436C /* Deck cc_praoden */ 3437 SUBROUTINE CC_PRAODEN(DENS,ISYDEN) 3438C 3439C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3440C 3441C Ove Christiansen 13-7-1995 3442C 3443C Purpose: Print density matrices 3444C calculated from ordinary lambda matrices 3445C and from C1 transformed lambda hole matrix. 3446C 3447C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3448C 3449#include "implicit.h" 3450 DIMENSION DENS(*) 3451#include "ccorb.h" 3452#include "ccsdsym.h" 3453#include "priunit.h" 3454C 3455 CALL AROUND('Lamda density matrix') 3456 KOFF1 = 1 3457 DO 110 ISYMB = 1,NSYM 3458C 3459 ISYMA = MULD2H(ISYMB,ISYDEN) 3460 WRITE(LUPRI,*) 'Symmetry number alfa,beta: ',ISYMA,ISYMB 3461 NBASA = NBAS(ISYMA) 3462 NBASB = NBAS(ISYMB) 3463 CALL OUTPUT(DENS(KOFF1),1,NBASA,1,NBASB,NBASA,NBASB,1,LUPRI) 3464 KOFF1 = KOFF1 + NBASA*NBASB 3465C 3466 110 CONTINUE 3467C 3468 END 3469C /* Deck cc_prfckao */ 3470 SUBROUTINE CC_PRFCKAO(FOCK,ISYFCK) 3471C 3472C Ove Christiansen 14-7-1994 3473C 3474C Purpose: Print AO Fock matrix. 3475C 3476#include "implicit.h" 3477 DIMENSION FOCK(*) 3478#include "ccorb.h" 3479#include "ccsdsym.h" 3480#include "ccsdinp.h" 3481#include "priunit.h" 3482C 3483 CALL AROUND('The AO Fock matrix') 3484 KOFF1 = 1 3485 DO 50 ISYMB = 1,NSYM 3486 ISYMA = MULD2H(ISYMB,ISYFCK) 3487 WRITE(LUPRI,*) 'Symmetry of A,B : ',ISYMA,ISYMB 3488 NBASA = NBAS(ISYMA) 3489 NBASB = NBAS(ISYMB) 3490 CALL OUTPUT(FOCK(KOFF1),1,NBASA,1,NBASB,NBASA,NBASB,1,LUPRI) 3491 KOFF1 = KOFF1 + NBAS(ISYMA)*NBAS(ISYMB) 3492 50 CONTINUE 3493C 3494 END 3495C /* Deck cc_prfckmo */ 3496 SUBROUTINE CC_PRFCKMO(FOCK,ISYFCK) 3497C 3498C Ove Christiansen 14-7-1994 3499C 3500C Purpose: Print MO Fock matrix. 3501C 3502#include "implicit.h" 3503 DIMENSION FOCK(*) 3504#include "ccorb.h" 3505#include "ccsdsym.h" 3506#include "ccsdinp.h" 3507#include "priunit.h" 3508C 3509 CALL AROUND('The MO Fock matrix ') 3510C 3511 KOFF1 = 1 3512 KOFF2 = NLRHFR(ISYFCK) + 1 3513 DO 300 ISYMQ = 1,NSYM 3514C 3515 ISYMP = MULD2H(ISYMQ,ISYFCK) 3516C 3517 WRITE(LUPRI,1) ISYMP,ISYMQ 3518 WRITE(LUPRI,2) 3519 WRITE(LUPRI,3) 3520 IF (NRHF(ISYMQ) .EQ. 0) THEN 3521 WRITE(LUPRI,4) 3522 GOTO 310 3523 ENDIF 3524 CALL OUTPUT(FOCK(KOFF1),1,NORB(ISYMP),1,NRHF(ISYMQ), 3525 * NORB(ISYMP),NRHF(ISYMQ),1,LUPRI) 3526 310 CONTINUE 3527C 3528 WRITE(LUPRI,5) 3529 WRITE(LUPRI,6) 3530 IF (NVIR(ISYMQ) .EQ. 0) THEN 3531 WRITE(LUPRI,4) 3532 GOTO 320 3533 ENDIF 3534 CALL OUTPUT(FOCK(KOFF2),1,NORB(ISYMP),1,NVIR(ISYMQ), 3535 * NORB(ISYMP),NVIR(ISYMQ),1,LUPRI) 3536C 3537 320 CONTINUE 3538C 3539 KOFF1 = KOFF1 + NORB(ISYMP)*NRHF(ISYMQ) 3540 KOFF2 = KOFF2 + NORB(ISYMP)*NVIR(ISYMQ) 3541C 3542 300 CONTINUE 3543C 3544 1 FORMAT(/,/,7X,'Symmetry of P,Q :',I5,I5) 3545 2 FORMAT(/,/,7X,'Occupied part') 3546 3 FORMAT(7X,'-------------') 3547 4 FORMAT(/,/,7X,'This symmetry is empty') 3548 5 FORMAT(/,/,7X,'Virtual part') 3549 6 FORMAT(7X,'------------') 3550C 3551 END 3552C /* Deck cc_rvec */ 3553 SUBROUTINE CC_RVEC(LU,FIL,LLEN,LEN,NR,VEC) 3554C 3555C Ove Christiansen 22-1-1996: 3556C 3557C Read vector NR on file FIL with unit number LU and 3558C put into VEC. The position is calculated relative to 3559C LLEN which is the length of the vectors according to 3560C which the vectors was put there in the first place. 3561C 3562C tbp: Use fortran I/O for Cholesky.... 3563C 3564C 3565#include "implicit.h" 3566Cholesky 3567#include "maxorb.h" 3568#include "ccdeco.h" 3569#include "priunit.h" 3570Cholesky 3571 DIMENSION VEC(*) 3572 CHARACTER FIL*(*) 3573C 3574 IF (CHOINT) THEN 3575 3576C Fortran I/O section. 3577C -------------------- 3578 3579 IF (LLEN .NE. LEN) THEN 3580 WRITE(LUPRI,'(//,5X,A,/,5X,A,I10,/,5X,A,I10,/)') 3581 & 'Fatal error in CC_RVEC: LLEN .ne. LEN:', 3582 & 'LLEN = ',LLEN, 3583 & 'LEN = ',LEN 3584 CALL QUIT('Fatal [FORTRAN] I/O error in CC_RVEC') 3585 ENDIF 3586 3587 CALL CHO_MOREAD(VEC,LEN,1,NR,LU) 3588 3589 ELSE 3590 3591C crayio section. 3592C --------------- 3593C 3594 IOFF = 1 + LLEN*(NR-1) 3595C 3596 IF (LEN .GT. 0) THEN 3597 CALL GETWA2(LU,FIL,VEC,IOFF,LEN) 3598C 3599 ENDIF 3600C 3601 END IF 3602C 3603 RETURN 3604 END 3605C /* Deck cc_wvec */ 3606 SUBROUTINE CC_WVEC(LU,FIL,LLEN,LEN,NR,VEC) 3607C 3608C Ove Christiansen 22-1-1996: 3609C 3610C Write vector NR given in VEC on file FIL with unit number LU. 3611C The position is calculated relative to 3612C LLEN which is the length of the vectors according to 3613C which the vectors was put there in the first place. 3614C 3615#include "implicit.h" 3616Cholesky 3617#include "maxorb.h" 3618#include "ccdeco.h" 3619#include "priunit.h" 3620Cholesky 3621 DIMENSION VEC(*) 3622 CHARACTER FIL*(*) 3623C 3624 IF (CHOINT) THEN 3625 3626C Fortran I/O section. 3627C -------------------- 3628 3629 IF (LLEN .NE. LEN) THEN 3630 WRITE(LUPRI,'(//,5X,A,/,5X,A,I10,/,5X,A,I10,/)') 3631 & 'Fatal error in CC_WVEC: LLEN .ne. LEN:', 3632 & 'LLEN = ',LLEN, 3633 & 'LEN = ',LEN 3634 CALL QUIT('Fatal [FORTRAN] I/O error in CC_WVEC') 3635 ENDIF 3636 3637 CALL CHO_MOWRITE(VEC,LEN,1,NR,LU) 3638 3639 ELSE 3640 3641C crayio section. 3642C --------------- 3643 3644 IOFF = 1 + LLEN*(NR-1) 3645C 3646 IF (LEN .GT. 0) THEN 3647 CALL PUTWA2(LU,FIL,VEC,IOFF,LEN) 3648 ENDIF 3649C 3650 ENDIF 3651C 3652 RETURN 3653 END 3654c*DECK CC_JACEXP 3655 SUBROUTINE CC_JACEXP(WORK,LWORK) 3656C 3657C---------------------------------------------------------------------- 3658C Calculates the CC Jacobian by Transformation of 3659C unit vectors and writes it to disk in CCLR_JAC. 3660C 3661C Written by Ove Christiansen 13-10-1995 3662C 3663C Christof Haettig, October 1998: 3664C changed to a dummy routine, because call of CCLR_RHTR and 3665C CCLHST via CCLR_TRR was switched off since quite some time... 3666C after change of F matrix implementation CCLR_TRR was 3667C completly dummy and removed... 3668C 3669C---------------------------------------------------------------------- 3670C 3671#include "implicit.h" 3672#include "dummy.h" 3673#include "maxorb.h" 3674#include "iratdef.h" 3675#include "ccorb.h" 3676#include "aovec.h" 3677#include "ccsdinp.h" 3678#include "cclr.h" 3679#include "ccsdsym.h" 3680#include "ccsdio.h" 3681#include "leinf.h" 3682#include "priunit.h" 3683C 3684 DIMENSION WORK(LWORK) 3685 PARAMETER (XHALF = 0.5D00,XMTWO = -2.0D00, DELTA = 1.0D-08) 3686C 3687CCH IF (IPRINT .GT. 10) THEN 3688 CALL AROUND(' START OF CC_JACEXP') 3689 CALL AROUND(' The routine is outdated... Leaving CC_JACEXP.') 3690 RETURN 3691CCH ENDIF 3692#ifdef OUTDATED_ROUTINE 3693C 3694C------------------- 3695C Open jac file. 3696C------------------- 3697C 3698 LUJAC = -1 3699 CALL GPOPEN(LUJAC,'CCLR_JAC','UNKNOWN',' ','UNFORMATTED',IDUMMY, 3700 * .FALSE.) 3701 REWIND(LUJAC) 3702C 3703C---------------------------- 3704C Work space allocations. 3705C---------------------------- 3706C 3707 NTAMP = NT1AM(ISYMTR) + NT2AM(ISYMTR) 3708C 3709 KV = 1 3710 KEND1 = 1 + NTAMP 3711 LWRK1 = LWORK - KEND1 3712 IF (LWRK1 .LT. 0 ) 3713 & CALL QUIT('INSUFFICIENT WORK SPACE IN CC_JACEXP') 3714C 3715C--------------------------------------------------------------------- 3716C Loop over vectors and save on disk the coloumns of the jacobian. 3717C--------------------------------------------------------------------- 3718C 3719 CALL DZERO(WORK(KV),NTAMP) 3720C 3721 JACEXP = .FALSE. 3722C 3723 DO 100 I = 1,NTAMP 3724C 3725 WORK(KV+I-1) = 1.0D00 3726CCH CALL CCLR_TRR(1,0,WORK(KV),DUMMY, 3727CCH * XLINLE,WORK(KEND1),LWRK1) 3728 WRITE(LUJAC) (WORK(KEND1+J-1),J=1,NTAMP) 3729 WORK(KV+I-1) = 0.0D00 3730C 3731 100 CONTINUE 3732C 3733 JACEXP = .TRUE. 3734C 3735 IF (IPRINT .GT. 10) THEN 3736 CALL AROUND(' END OF CC_JACEXP ') 3737 ENDIF 3738C 3739 CALL GPCLOSE(LUJAC,'KEEP') 3740C 3741 RETURN 3742#endif 3743 END 3744C 3745C /* Deck cc_pronelao */ 3746 SUBROUTINE CC_PRONELAO(FOCK,ISYFCK) 3747C 3748C Purpose: Print AO one-electron matrix. 3749C 3750#include "implicit.h" 3751 DIMENSION FOCK(*) 3752#include "ccorb.h" 3753#include "ccsdsym.h" 3754#include "ccsdinp.h" 3755#include "priunit.h" 3756C 3757 KOFF1 = 1 3758 DO ISYMB = 1,NSYM 3759 ISYMA = MULD2H(ISYMB,ISYFCK) 3760 WRITE(LUPRI,*) 'Symmetry of A,B : ',ISYMA,ISYMB 3761 NBASA = NBAS(ISYMA) 3762 NBASB = NBAS(ISYMB) 3763 CALL OUTPUT(FOCK(KOFF1),1,NBASA,1,NBASB,NBASA,NBASB,1,LUPRI) 3764 KOFF1 = KOFF1 + NBAS(ISYMA)*NBAS(ISYMB) 3765 END DO 3766C 3767 END 3768C 3769C /* Deck cc_prpr12 */ 3770 SUBROUTINE CC_PRPR12(TR12AM,ISYM,NTR12,LHEAD) 3771C 3772C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3773C Christian Neiss, Feb. 2005 3774C PRint Packed R12-vector - PRPR12 (in general non. tot. sym.) 3775C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3776C 3777#include "implicit.h" 3778#include "ccorb.h" 3779#include "ccsdsym.h" 3780#include "ccsdinp.h" 3781#include "priunit.h" 3782C 3783 LOGICAL LHEAD 3784 DIMENSION TR12AM(*) 3785 3786C 3787C------------------------------------ 3788C Write R12 vector. 3789C------------------------------------ 3790C 3791 DO 20 I = 1,NTR12 3792C 3793 IF (LHEAD) CALL AROUND('R12 part of vector ') 3794 IF (NTR12.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 3795C 3796 DO 200 ISYMLJ = 1,NSYM 3797C 3798 ISYMKI = MULD2H(ISYMLJ,ISYM) 3799C 3800 KOFF = ITR12AM(ISYMKI,ISYMLJ) + 1 3801 3802 IF (ISYMKI.EQ.ISYMLJ) THEN 3803C 3804 WRITE(LUPRI,*) 'Symmetry block number(ki,lj): ', 3805 & ISYMKI,ISYMLJ 3806 IF (NMATKI(ISYMKI).EQ.0) THEN 3807 WRITE(LUPRI,*) 'This symmetry is empty' 3808 ELSE 3809 CALL OUTPAK(TR12AM(KOFF),NMATKI(ISYMKI),1,LUPRI) 3810 END IF 3811 WRITE(LUPRI,*) ' ' 3812C 3813 ELSE IF (ISYMLJ.GT.ISYMKI) THEN 3814C 3815 WRITE(LUPRI,*) 'Symmetry block number(ki,lj): ', 3816 & ISYMKI,ISYMLJ 3817 IF (NMATKI(ISYMKI).EQ.0 .OR. NMATKI(ISYMLJ).EQ.0) THEN 3818 WRITE(LUPRI,*) 'This symmetry is empty' 3819 ELSE 3820 CALL OUTPUT(TR12AM(KOFF),1,NMATKI(ISYMKI), 3821 & 1,NMATKI(ISYMLJ), 3822 * NMATKI(ISYMKI),NMATKI(ISYMLJ),1,LUPRI) 3823 END IF 3824 WRITE(LUPRI,*) ' ' 3825C 3826 ENDIF 3827C 3828 200 CONTINUE 3829C 3830 20 CONTINUE 3831C 3832 RETURN 3833 END 3834C 3835C /* Deck cc_prsqr12 */ 3836 SUBROUTINE CC_PRSQR12(TR12AM,ISYM,TRANS,NTR12,LHEAD) 3837C 3838C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3839C Christian Neiss Feb. 2005 3840C PRint SQuared R12 vector. general non. tot. sym. 3841C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3842C 3843#include "implicit.h" 3844#include "ccorb.h" 3845#include "ccsdsym.h" 3846#include "ccsdinp.h" 3847#include "priunit.h" 3848C 3849 LOGICAL LHEAD 3850 CHARACTER*1 TRANS 3851 DIMENSION TR12AM(*) 3852C 3853 IF ((TRANS.NE.'T').AND.(TRANS.NE.'N')) 3854 & CALL QUIT('Illegal value for "TRANS" in CC_PRSQR12') 3855C 3856C 3857C------------------------------------ 3858C Write R12 vector. 3859C------------------------------------ 3860C 3861 DO 20 I = 1, NTR12 3862C 3863 IF (LHEAD) CALL AROUND('R12 part of vector ') 3864 3865 IF (NTR12.GT.1) WRITE(LUPRI,*) 'Vector Number:',I 3866C 3867 DO 200 ISYMIJ = 1,NSYM 3868C 3869 ISYMKL = MULD2H(ISYMIJ,ISYM) 3870C 3871 IF (TRANS.EQ.'T') THEN 3872 WRITE(LUPRI,*) 'Symmetry block (ij,kl):', 3873 & ISYMIJ,ISYMKL 3874 ELSE IF (TRANS.EQ.'N') THEN 3875 WRITE(LUPRI,*) 'Symmetry block (kl,ij):', 3876 & ISYMKL,ISYMIJ 3877 END IF 3878 IF (NMATKL(ISYMKL).EQ.0 .OR. NMATIJ(ISYMIJ).EQ.0) THEN 3879 WRITE(LUPRI,*) 'This symmetry is empty' 3880 ELSE 3881 IF (TRANS.EQ.'T') THEN 3882 KOFF = ITR12SQT(ISYMIJ,ISYMKL) + 1 3883 CALL OUTPUT(TR12AM(KOFF),1,NMATIJ(ISYMIJ),1, 3884 * NMATKL(ISYMKL),NMATIJ(ISYMIJ),NMATKL(ISYMKL), 3885 * 1,LUPRI) 3886 ELSE IF (TRANS.EQ.'N') THEN 3887 KOFF = ITR12SQ(ISYMKL,ISYMIJ) + 1 3888 CALL OUTPUT(TR12AM(KOFF),1,NMATKL(ISYMKL),1, 3889 * NMATIJ(ISYMIJ),NMATKL(ISYMKL),NMATIJ(ISYMIJ), 3890 * 1,LUPRI) 3891 END IF 3892 END IF 3893C 3894 WRITE(LUPRI,*) ' ' 3895C 3896 200 CONTINUE 3897C 3898 20 CONTINUE 3899C 3900 RETURN 3901 END 3902C 3903C /* Deck cclr_transsq */ 3904 SUBROUTINE CCLR_TRANSSQ(MATRIX,NROWS) 3905C 3906C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3907C Christian Neiss Mar. 2005 3908C Transpose a square matrix 'in place' 3909C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3910C 3911 implicit none 3912#include "priunit.h" 3913C 3914 INTEGER IDXIJ,IDXJI,NROWS,I,J 3915#if defined (SYS_CRAY) 3916 REAL X1,X2,matrix(*) 3917#else 3918 DOUBLE PRECISION X1,X2,matrix(*) 3919#endif 3920C 3921 do i = 1, nrows 3922 do j = 1, i 3923 idxij = nrows*(j-1) + i 3924 idxji = nrows*(i-1) + j 3925 x1 = matrix(idxij) 3926 x2 = matrix(idxji) 3927 matrix(idxij) = x2 3928 matrix(idxji) = x1 3929 end do 3930 end do 3931 3932 RETURN 3933 END 3934C 3935C /* Deck cclr_trsqr12 */ 3936 SUBROUTINE CCLR_TRSQR12(MATRIX,ISYM) 3937C 3938C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3939C Christian Neiss April 2005 3940C Transpose a symmetry packed square matrix X_kl^ij: 3941C X(kl,ij) --> X(ij,kl) 3942C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 3943C 3944 implicit none 3945#include "priunit.h" 3946#include "ccsdsym.h" 3947#include "ccorb.h" 3948C 3949 INTEGER IDXKLIJ,IDXIJKL,IDXKL,IDXIJ,ISYM,ISYMKL,ISYMIJ 3950#if defined (SYS_CRAY) 3951 REAL X1,X2,matrix(*) 3952#else 3953 DOUBLE PRECISION X1,X2,matrix(*) 3954#endif 3955 LOGICAL LOCDBG 3956 PARAMETER (LOCDBG = .FALSE.) 3957C 3958 if (locdbg) then 3959 do I = 1, NSYM 3960 write(lupri,*) 'ITR12SQ = ', (ITR12SQ(I,J), J=1,NSYM) 3961 write(lupri,*) 'ITR12SQT= ', (ITR12SQT(I,J),J=1,NSYM) 3962 end do 3963 call around('Input Matrix in CCLR_TRSQR12') 3964 call cc_prsqr12(matrix,isym,'N',1,.false.) 3965 end if 3966C 3967 do isymkl = 1, nsym 3968 isymij = muld2h(isym,isymkl) 3969 do idxkl = 1, nmatkl(isymkl) 3970 do idxij = 1, nmatij(isymij) 3971 idxklij = itr12sq(isymkl,isymij) + 3972 & nmatkl(isymkl)*(idxij-1) + idxkl 3973 idxijkl = itr12sqt(isymij,isymkl) + 3974 & nmatij(isymij)*(idxkl-1) + idxij 3975 if (idxklij.lt.idxijkl) then 3976 x1 = matrix(idxklij) 3977 x2 = matrix(idxijkl) 3978 matrix(idxklij) = x2 3979 matrix(idxijkl) = x1 3980 end if 3981 end do 3982 end do 3983 end do 3984C 3985 if (locdbg) then 3986 call around('Transposed Matrix in CCLR_TRSQR12') 3987 call cc_prsqr12(matrix,isym,'T',1,.false.) 3988 end if 3989C 3990 RETURN 3991 END 3992C 3993