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 cceq_sol */ 20 SUBROUTINE CCEQ_SOL(LIST,LPROJECT,ISTATPRJ,ECURR, 21 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 22 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 23 * LREDS,REDS,FS12AM,LUFS12,FS2AM,LUFS2, 24 * LINQCC,TRIPLET,ISIDE,ISTRVE,NVEC,NUPVEC, 25 * NREDH,REDH,EIVAL,SOLEQ, 26 * WRK,LWRK,APROXR12) 27C 28C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 29C 30C input 31C 32C FRHO1,FRHO2,FRHO12,FC1AM,FC2AM,FC12AM are file names for files where 33C transformed vectors (rho1,rho2,rhor12) and trial vectors 34C (c1am,c2am,cr12am) are stored. 35C LUFR1,LUFR2,LUFR12,LUFC1,LUFC2,LUFC12 are the corresponding unit 36C numbers. 37C 38C TRIPLET.EQ.T solve triplet equations 39C .EQ.F solve singlet equations 40C 41C ISIDE.EQ.1 SOLVE EQUATIONS FROM RIGHT 42C ISIDE.EQ.-1 SOLVE EQUATIONS FROM LEFT 43C 44C LPROJECT = LOGICAL FOR PROJECTION, (sonia) 45C ISTATPRJ = STATE INDEX FOR PROJECTION (sonia) 46C 47C LINQCC.EQ.T ,SOLVE SET OF NVEC LINEAR EQUATIONS 48C ISTRVE = START NUMBER IN THE LIST FOR SOLVING LINEAR EQUATIONS 49C NVEC = NUMBER OF EQUATIONS THAT MUST BE SOLVED 50C NUPVEC = NUMBER OF CURRENT LINEAR INDEPENDENT NEW VECTORS. 51C 52C EIVAL CONTAINS NVEC FREQUENCIES 53C EIVAL( ISTRVE ... (ISTRVE+NVEC-1) ) 54C 55C LINQCC.EQ.F ,SOLVE EIGENVALUE EQUATIONS FOR NVEC LOWEST ROOTS 56C 57C IF (LINQCC) THEN 58C 59C THIS SUBROUTINE DIRECT THE SOLUTIONS OF NVECS SIMULTANEOUS 60C SET OF NEWTON-RAPHSON EQUATIONS 61C 62C IF (ISIDE.EQ.1) THEN 63C 64C (A-EIVAL*I)*X(J) + Fright(J) = 0 65C 66C ELSE IF ( ISIDE.EQ.-1) THEN 67C 68C [ (A-EIVAL*I) ] T *X(J) + Fleft(J) = 0 69C ( the transposed reduced matrix is constructed, this is 70C done without change in code because the linear 71C transformations are from the left, the reduced space is 72C set up as 73C REDH(I,J) = B(I) * S(J) ) 74C 75C END IF 76C ELSE 77C 78C SOLVE for NVECS lowest eigenvalues 79C 80C IF (ISIDE.EQ.1) THEN 81C 82C (A-EIVAL*I)*X(J) = 0 83C 84C ELSE IF ( ISIDE.EQ.-1) THEN 85C 86C [ (A-EIVAL*I) ] T *X(J) = 0 87C ( the transposed reduced matrix is constructed, this is 88C done without change in code because the linear 89C transformations are from the left) 90C 91C END IF 92C 93C END IF 94C 95C 96C---------------------------------------------------------------- 97C CC_TRDRV carry out linear transformations on new trial vectors 98C on LUFC1,LUFC2,LUFC12 and return linear transformed 99C trial vectors on LUFR1,LUFR2,LUFR12 100C CCRED finds the solution in reduced space. 101C CCNEX finds residual and the next trial vectors (saved on disk) 102C 103C Work space flow: 104C The whole work space is released after each routine is called 105C all communication takes place via file with trial vectors 106C and linear transformed vectors on LUFC1,LUFC2,LUFC12 and 107C LUFR1,LUFR2,LUFR12 respectively 108C 109C MAXLE: MAXIMUM NUMBER OF MICROITERATIONS 110C from common LEINF 111C MAXRED MAXIMUM DIMENSION OF REDUCED SPACE from common CCLR 112C---------------------------------------------------------------- 113C 114C---------------------------------------------------------------- 115C 116C JK&OC 080802: 117C Dirty hack to allow small number of iterations for CCSLV/MM calc. 118C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 119C 120 USE PELIB_INTERFACE, ONLY: USE_PELIB 121#include "implicit.h" 122#include "priunit.h" 123#include "ccsdinp.h" 124#include "ccsections.h" 125#include "cclr.h" 126#include "leinf.h" 127#include "ccslvinf.h" 128Cholesky 129#include "maxorb.h" 130#include "ccdeco.h" 131C 132C#include "ccsdsym.h" 133CSonia: 134! Consider moving the CVSEPARA in the argument list for better control 135! when not needed 136#include "ccexcicvs.h" 137 PARAMETER (TOL = 1.0D-14) 138Cholesky 139 PARAMETER (D0 = 0.0D0 ) 140C 141 CHARACTER*(*) FRHO1,FRHO2,FC1AM,FC2AM,FRHO12,FC12AM,FS12AM,FS2AM 142 CHARACTER*(*) LIST 143 CHARACTER*3 APROXR12 144 DIMENSION REDH(*), SOLEQ(*), REDS(*) 145 DIMENSION EIVAL(*) 146 DIMENSION WRK(*) 147 LOGICAL LINQCC, LPROJECT, TRIPLET, LREDS 148C 149 CALL QENTER('CCEQ_SOL') 150 IF (IPRLE .GT. 3 ) THEN 151 WRITE(LUPRI,*) 'CCEQ_SOL: IPRLE = ',IPRLE 152 WRITE(LUPRI,*) 'CCEQ_SOL: TRIPLET = ',TRIPLET 153 WRITE(LUPRI,*) 'CCEQ_SOL: LINQCC = ',LINQCC 154 WRITE(LUPRI,*) 'CCEQ_SOL: NREDH = ',NREDH 155 WRITE(LUPRI,*) 'CCEQ_SOL: MAXRED = ',MAXRED 156 WRITE(LUPRI,*) 'CCEQ_SOL: MAXLE = ',MAXLE 157 WRITE(LUPRI,*) 'CCEQ_SOL: ISIDE = ',ISIDE 158 WRITE(LUPRI,*) 'CCEQ_SOL: ISTRVE = ',ISTRVE 159 WRITE(LUPRI,*) 'CCEQ_SOL: NVEC = ',NVEC 160 WRITE(LUPRI,*) 'CCEQ_SOL: NUPVEC = ',NUPVEC 161 WRITE(LUPRI,*) 'CCEQ_SOL: NCCVAR = ',NCCVAR 162 WRITE(LUPRI,*) 'CCEQ_SOL: CVSEPARA= ',LCVSEXCI 163 WRITE(LUPRI,*) 'CCEQ_SOL: LIONISAT= ',LIONIZEXCI 164 WRITE(LUPRI,*) 'CCEQ_SOL: REMOVE_CORE = ',LRMCORE 165 ENDIF 166C 167 CALL GETTIM(CSTR,WSTR) 168C 169 TIMMIC = SECOND() 170 TIMLIN = D0 171 TIMRED = D0 172 TIMNEX = D0 173 ITLE = 0 174C 175Cholesky 176C 177C Check that all frequencies in EIVAL are identical for Cholesky CC2. 178C 179 IF (CHOINT .AND. CC2.AND. (NVEC.GT.0)) THEN 180 IF (.NOT. CHOEXC) THEN 181 ICOUN = 0 182 REFRQ = EIVAL(1) 183 DO I = 2,NVEC 184 DIFF = DABS(EIVAL(I) - REFRQ) 185 IF (DIFF .GT. TOL) ICOUN = ICOUN + 1 186 ENDDO 187 IF (ICOUN .GT. 0) THEN 188 WRITE(LUPRI,'(//,A)') 189 & 'CCEQ_SOL: frequencies MUST be identical for all', 190 & ' perturbations for Cholesky CCC2.' 191 WRITE(LUPRI,'(A)') 192 & 'Frequencies passed to CCEQ_SOL in EIVAL array:' 193 WRITE(LUPRI,'(4D20.12)') (EIVAL(I), I = 1,NVEC) 194 CALL QUIT('CCEQ_SOL: frequency error for Cholesky CC2') 195 ENDIF 196 ELSE 197 DO I = 1,NVEC 198 EIVAL(I) = ECURR 199 END DO 200 END IF 201 ENDIF 202C 203Cholesky 204C 205C print banner for solver: 206C 207 IF (IPRLE.GT.0) THEN 208 WRITE (LUPRI,'(///A/)') 209 & ' ----- COUPLED CLUSTER RESPONSE SOLVER -----' 210Chol 211 IF (CHOINT .AND. CC2.AND. (NVEC.GT.0)) 212 & WRITE(LUPRI,'(A,F12.5,/)') 213 & ' Frequnecy :',EIVAL(1) 214Chol 215 WRITE (LUPRI,'(3X,A)') 216 & ' Iter #Vectors time (min) residual' 217 WRITE (LUPRI,'(3X,A)') 218 & ' --------------------------------------' 219 END IF 220C 221C 222C Loop over micro iterations. 223C 224C 225 100 CONTINUE 226 ITLE = ITLE + 1 227C 228 IF (IPRLE .GE. 2) TIMIT = SECOND() 229C 230 TIM = SECOND() 231C 232C-------------------------------- 233C LINEAR TRANSFORMATIONS ON NUPVEC TRIAL VECTORS ON LUB 234C THE LINEAR TRANSFORMED VECTORS ARE RETURNED ON LUS 235C-------------------------------- 236C 237 IF (IPRLE .GT. 5) THEN 238 WRITE (LUPRI,*) 239 * ' --- Call CC_TRDRV with TRIPLET and ECURR set to:', 240 * TRIPLET, ECURR 241 ENDIF 242C 243Cholesky Pass EIVAL (before it was DUMMY) to CC_TRDRV 244Cholesky to be used by CC_CHOATR 245C 246 IST = NREDH - NUPVEC + 1 247 CALL CC_TRDRV(ECURR,FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 248 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 249 * TRIPLET,.FALSE.,IDUMMY,EIVAL,FS12AM,LUFS12, 250 * FS2AM,LUFS2, 251 * ISIDE,IST,NUPVEC,WRK,LWRK,APROXR12) 252C 253 CALL FLSHFO(LUPRI) 254 TIM = SECOND() - TIM 255 TIMLIN = TIMLIN + TIM 256 NUPVECSAVE = NUPVEC 257C 258C 259C-------------------------------- 260C CALL CCRED(), INCREASE DIMENSION OF REDUCED SPACE, 261C AND FIND SOLUTIONS 262C-------------------------------- 263C 264 TIMRED = TIMRED - SECOND() 265 266 CALL CCRED(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 267 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 268 * CCR12,FS12AM,LUFS12,FS2AM,LUFS2, 269 * LINQCC,TRIPLET,ISIDE,LIST, 270 * LPROJECT,ISTATPRJ,LREDS,REDS, 271 * REDH,NREDH,ISTRVE,NUPVEC,NVEC, 272 * EIVAL,SOLEQ,WRK,LWRK,DEBUG) 273 CALL FLSHFO(LUPRI) 274 TIMRED = TIMRED + SECOND() 275C 276C CREATE NVEC NEW LINEAR INDEPENDENT TRIAL VECTORS IN CCNEX() 277C FROM THE SOLUTIONS OF THE REDUCED EQUATIONS 278C REDUCE THE NUMBER FOR CONVERGED VECTORS AND IF LINEAR 279C DEPENDENCE IS ENCOUNTERED. 280C NUPVEC IS NVEC REDUCED FOR CONVERGENCE AND 281C LINEAR DEPENDENCE 282C 283 IF (ITLE .LT. MAXLE) THEN 284 JCONV = 0 285 ELSE 286 JCONV = -1 287 END IF 288C 289C-------------------------------- 290C Find the next trial vector. 291C-------------------------------- 292C 293 TIMNEX = TIMNEX - SECOND() 294 CALL CCNEX(LIST,ITLE,LPROJECT,ISTATPRJ, 295 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 296 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 297 * FS12AM,LUFS12,FS2AM,LUFS2, 298 * LINQCC,TRIPLET,ISIDE,NREDH,ISTRVE, 299 * NVEC,NUPVEC,JCONV,EIVAL,SOLEQ,WRK,LWRK,RMXNORM) 300C 301C----------------------------------------------------- 302C Print timing & convergence statistic 303C----------------------------------------------------- 304C 305 IF (IPRLE.GT.0) THEN 306 WRITE (LUPRI,'(2X,I5,3X,I5,1X,F12.2,1X,E12.2)') 307 * ITLE,NUPVECSAVE,TIM/60.0D0,RMXNORM 308 END IF 309 CALL FLSHFO(LUPRI) 310C 311C----------------------------------------------------- 312C Print out the lowest eigenvalues, Ove 24-9-1995. 313C----------------------------------------------------- 314C 315 IF (IPRLE .GT. 2 ) THEN 316C 317 IF (.NOT. LINQCC) THEN 318 WRITE (LUPRI,'(/A)') 319 * ' Lowest eigenvalues in this iteration: ' 320 CALL OUTPUT(EIVAL,1,NVEC,1,1,NVEC,1,1,LUPRI) 321 ENDIF 322C 323 WRITE (LUPRI,'(A,I5)') 324 * ' Dimension of the reduced space: ',NREDH 325C 326 ENDIF 327C 328 CALL FLSHFO(LUPRI) 329 TIMNEX = TIMNEX + SECOND() 330C 331 IF (IPRLE .GE. 2) THEN 332 TIMIT = SECOND() - TIMIT 333 WRITE (LUPRI,'(A,F12.2,/)') 334 * ' --- TIME USED IN THIS MICRO ITERATION',TIMIT 335 END IF 336C 337 IF (JCONV.LT.0) THEN 338C ( LINEAR DEPENDENCE BETWEEN ALL NEW TRIAL VECTORS ) 339 WRITE (LUPRI,'(/A/A)') 340 * ' *** CCEQ_SOL -MICROITERATIONS STOPPED ', 341 * ' LINEAR DEPENDENCE BETWEEN ALL NEW TRIAL VECTORS' 342CCN IF (RMXNORM.GT.THRLE) THEN 343CCN WRITE(LUPRI,*) 344CCN * ' BUT NOT ALL SOLUTION VECTORS ARE'// 345CCN * ' CONVERGED!' 346CCN CALL QUIT(' SOLUTION VECTORS NOT CONVERGED ') 347CCN END IF 348 ELSE IF (JCONV.GT.0) THEN 349C (CONVERGED) 350 IF (IPRLE .GT. 10) THEN 351 WRITE(LUPRI,'(/,(A,I4,A,A,I4,A))') 352 * ' MICROITERATIONS CONVERGED FOR',NVEC, 353 * ' SOLUTION VECTORS', 354 * ' IN',ITLE,' ITERATIONS.' 355 END IF 356 IF (CCSLV.OR.USE_PELIB()) NSLVINIT = NSLVINIT + NREDH 357 IF (CCSLV.OR.USE_PELIB()) LRSPFUL = .TRUE. 358 ELSE IF (NREDH + NUPVEC .GT. MAXRED) THEN 359C (NOT CONVERGED) 360 WRITE(LUPRI,'(/A,I5,A)') 361 * ' *** CCEQ_SOL-MAXIMUM DIMENSION OF REDUCED EQUATIONS', 362 * MAXRED,' REACHED.' 363 WRITE(LUPRI,*)' ALL SOLUTION VECTORS NOT CONVERGED ' 364 CALL QUIT(' ALL SOLUTION VECTORS NOT CONVERGED ') 365 ELSE IF (ITLE .LT. MAXLE) THEN 366 CALL FLSHFO(LUPRI) 367 IF (LINQCC.AND.(CCSLV.OR.USE_PELIB()).AND. 368 & (ITLE.GE.MXLINIT).AND. 369 & (LIST(1:1) .EQ. 'L')) THEN 370C additions to reduced space should be ignored. 371C 372 NREDH = NREDH-NUPVEC 373 NSLVINIT = NSLVINIT + NREDH 374 LRSPFUL = .FALSE. 375 WRITE(LUPRI,*) 376 * 'CCSLV/PE: We stop for now though not fully converged yet' 377 WRITE(LUPRI,*) 378 * ' Accumulated inner iterations are:', NSLVINIT 379C WRITE(LUPRI,*) 'last - NREDH,NUPVEC:',NREDH,NUPVEC 380 ELSE 381 GO TO 100 382 END IF 383 ELSE 384 WRITE(LUPRI,'(/A,I5,A)') 385 * ' *** CCEQ_SOL-MAXIMUM NUMBER OF MICROITERATIONS', 386 * ITLE,' REACHED.' 387 CALL QUIT(' *** CCEQ_SOL-MAX. MICROITERATIONS REACHED') 388 END IF 389C 390C===================== 391C End of 100 Loop. 392C===================== 393C 394 TIMMIC = SECOND() - TIMMIC 395 CALL GETTIM(CEND,WEND) 396 CTOT = CEND - CSTR 397 WTOT = WEND - WSTR 398C 399 IF (IPRLE .GT. 0) THEN 400 WRITE (LUPRI,'(3X,A)') 401 & ' --------------------------------------' 402 WRITE (LUPRI,'(3X,A,I3,A)') ' converged in',ITLE,' iterations' 403 WRITE (LUPRI,'(3X,A,E12.2)') ' threshold:',THRLE 404 WRITE (LUPRI,'(//T10,A)') 'Routine Time (min)' 405 WRITE (LUPRI,'(T10,A)') '---------------------------' 406 WRITE (LUPRI,'(T10,A,F15.2)') 'CC_TRDRV ',TIMLIN/60.0D0 407 WRITE (LUPRI,'(T10,A,F15.2)') 'CCRED ',TIMRED/60.0D0 408 WRITE (LUPRI,'(T10,A,F15.2)') 'CCNEX ',TIMNEX/60.0D0 409 WRITE (LUPRI,'(T10,A)') '---------------------------' 410 WRITE (LUPRI,'(T10,A,F14.2,//)') 'Total time',TIMMIC/60.0D0 411 CALL TIMTXT(' Total CPU time used in CCEQ_SOLV:',CTOT, 412 & LUPRI) 413 CALL TIMTXT(' Total wall time used in CCEQ_SOLV:',WTOT, 414 & LUPRI) 415 WRITE (LUPRI,'(//A/)') 416 & ' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx' 417 END IF 418C 419 IF (IPRLE .GT. 0 .AND. .NOT. LINQCC) THEN 420 WRITE (LUPRI,'(//A/A/A,1P,D10.2,A/)') 421 * ' Final eigenvalues :', 422 * ' =====================', 423 * ' (convergence threshold:',THRLE,')' 424 DO 920 I = 1,NVEC 425 WRITE (LUPRI,'(I5,F15.10)') I,EIVAL(I) 426 920 CONTINUE 427 END IF 428 IF (IPRLE .GE. 10) THEN 429 IF (LINQCC) THEN 430 WRITE (LUPRI,'(//A)') 431 * ' FINAL LINEAR SOLUTION VECTORS IN REDUCED BASIS FOR :' 432 ELSE 433 WRITE (LUPRI,'(//A)') 434 * ' FINAL EIGENVECTORS IN REDUCED BASIS FOR :' 435 END IF 436 IOFF = 0 437 DO 930 I = 1,NVEC 438 IF (LINQCC) THEN 439 WRITE (LUPRI,'(/A,I4/)') ' - GRADIENT VECTOR NO.',I 440 ELSE 441 WRITE (LUPRI,'(/A,I4,A,1P,D15.6/A)') 442 * ' - EIGENVALUE NO.',I, ' - EIGENVALUE',EIVAL(I), 443 * ' - EIGENVECTOR :' 444 END IF 445 WRITE (LUPRI,'(5F15.8)') (SOLEQ(IOFF+J),J=1,NREDH) 446 IOFF = IOFF + MAXRED 447 930 CONTINUE 448 END IF 449 CALL FLSHFO(LUPRI) 450C 451 CALL QEXIT('CCEQ_SOL') 452 RETURN 453 END 454C /* Deck cceq_str */ 455 SUBROUTINE CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 456 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 457 * TRIPLET,ISIDE,ISTRVE,NVEC,NUPVEC, 458 * NREDH,EIVAL,IPLACE,WORK,LWORK,LIST) 459C 460C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 461C 462C Based on original EVALST routine. 463C 464C Modified by Ove Christiansen to calculate diagonal when needed. 465C Included projection, Sonia Coriani 1997 466C Included CVS, Sonia Coriani 2015 467C 468C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 469C 470#include "implicit.h" 471 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 , 472 * THRLDP = 1.0D-20 , THRDIA = 1.0D-6 ) 473#include "priunit.h" 474#include "ibndxdef.h" 475#include "ccorb.h" 476#include "cclr.h" 477#include "ccinf.h" 478#include "ccsdinp.h" 479#include "ccsdsym.h" 480#include "leinf.h" 481#include "ccinf.h" 482#include "ccexci.h" 483#include "dummy.h" 484#include "r12int.h" 485C 486#include "ccexcicvs.h" 487 488 LOGICAL LOCDBG, LNOREAD 489 PARAMETER (LOCDBG = .FALSE.) 490 CHARACTER*(*) FC1AM,FC2AM,FC12AM,FS12AM,FS2AM 491 CHARACTER*(*) LIST 492 CHARACTER*10 MODEL 493 DIMENSION IPLACE(*) , WORK(*),EIVAL(*) 494 LOGICAL NEWSTR, LPROJECT, TRIPLET 495C 496 CALL QENTER('CCEQ_STR') 497 IF (IPRLE. GT. 10 ) THEN 498 CALL AROUND( ' START OF CCEQ_STR ' ) 499 WRITE(LUPRI,*) 'CCEQ_STR: LIST =',LIST(1:3) 500 WRITE(LUPRI,*) 'CCEQ_STR: TRIPLET =',TRIPLET 501 WRITE(LUPRI,*) 'CCEQ_STR: LPROJECT =',LPROJECT 502 WRITE(LUPRI,*) 'NCCVAR =',NCCVAR 503 WRITE(LUPRI,*) 'CCEQ_STR: NCCVAR =',NCCVAR 504 WRITE(LUPRI,*) 'CCEQ_STR: CVSEPARA =',LCVSEXCI 505 WRITE(LUPRI,*) 'CCEQ_STR: LIONISAT =',LIONIZEXCI 506 WRITE(LUPRI,*) 'CCEQ_STR: REMOVE_CORE = ',LRMCORE 507 WRITE(LUPRI,*) 'CCEQ_STR: LCOR in GS = ',LCOR 508 CALL FLSHFO(LUPRI) 509 ENDIF 510 CALL FLSHFO(LUPRI) 511C 512 MODEL = 'CCSD ' 513 IF (CCS) MODEL = 'CCS ' 514 IF (CIS) MODEL = 'CIS ' 515 IF (CC2) MODEL = 'CC2 ' 516 IF (CC3) MODEL = 'CC3 ' 517 IF (CC1A) MODEL = 'CC1A ' 518C 519 IF (TRIPLET) THEN 520 IMULT = 3 521 ELSE 522 IMULT = 1 523 ENDIF 524C 525C ---------------------------- 526C Initialize NUPVEC and NREDH: 527C ---------------------------- 528C 529 NUPVEC = 0 530 NREDH = 0 531C 532 KV = 1 533 KEND1 = KV + NCCVAR 534 LEND1 = LWORK - KEND1 535C 536 KDIA = KEND1 537 KEND2 = KDIA + NCCVAR 538 LEND2 = LWORK - KEND2 539 IF (LEND2 .LT. 0) CALL QUIT('Too little workspace in CCEQ_STR-1') 540C 541C ------------------------------------- 542C Get an approximate diagonal Jacobian: 543C ------------------------------------- 544 IF (LIST(2:2) .EQ.'E') THEN 545 CALL CCLR_DIA(WORK(KDIA),ISYMTR,TRIPLET,WORK(KEND2),LEND2) 546 IF (LOCDBG) THEN 547 CALL AROUND(' CCEQ_STR: approximate diagonal Jacobian') 548 CALL OUTPUT(WORK(KDIA),1,NCCVAR,1,1,NCCVAR,1,1,LUPRI) 549 ENDIF 550 IF (LCVSEXCI.or.LIONIZEXCI) THEN 551 write(lupri,*)'CCEQ_STR: freeze VALENCE OR VIRT' 552 if (triplet) then 553 KCAM1 = KDIA 554 KCAMP = KCAM1 + NT1AM(ISYMTR) 555 KCAMM = KCAMP + NT2AM(ISYMTR) 556 CALL CC_FREEZE_TRIPLETstart(WORK(KCAM1), 557 & WORK(KCAMP),WORK(KCAMM),ISYMTR, 558 & MAXCORE,MAXION, 559 & NRHFCORE,IRHFCORE, 560 & NVIRION,IVIRION,LBOTHEXCI) 561 else 562 CALL CC_FREEZE_start(WORK(KDIA),ISYMTR, 563 & MAXCORE,MAXION, 564 & NRHFCORE,IRHFCORE, 565 & NVIRION,IVIRION,LBOTHEXCI) 566 end if 567 END IF 568 ENDIF 569C 570C --------------------------------------------------------------- 571C For eigenvalue problems select elements for unit start vectors: 572C --------------------------------------------------------------- 573 IF (STVEC.AND.(LIST(2:2) .EQ.'E')) THEN 574 WRITE (LUPRI,'(/A/A/A)') 575 * ' CC elements selected for start vectors:', 576 * ' =======================================', 577 * ' Start vector no. CC element no.' 578 DO I = 1,NVEC 579 IPLACE(I) = ISTVEC(I,ISYMTR) 580 WRITE (LUPRI,'(I10,I20)') I,IPLACE(I) 581 END DO 582 ELSE IF (LIST(2:2) .EQ.'E') THEN 583 MXELMN = MIN(NVEC+3,NCCVAR,MAXRED) 584 CALL FNDMN3(WORK(KDIA),NCCVAR,MXELMN,IPLACE, 585 * NELMN,IPRLE,THRDIA) 586 DO I = 1, NVEC 587 IF (WORK(KDIA-1+IPLACE(I)) .GT. 1.0D6) THEN 588 WRITE (LUPRI,*) 'Problem in CCEQ_STR:' 589 WRITE (LUPRI,*) 'not enough allowed elements in vector:' 590 WRITE (LUPRI,*) 'requested:',NVEC 591 WRITE (LUPRI,*) 'found :',I-1 592 WRITE (LUPRI,*) 'reset NVEC variable to ',I-1 593 NELMN = MIN(NELMN,I-1) 594 NVEC = MIN(NVEC,I-1) 595 END IF 596 END DO 597 END IF 598C 599C-------------------------------------------------------- 600C Add one vector at the time to list of startvectors. 601C-------------------------------------------------------- 602C 603 DO 100 I = 1,NVEC 604 NEWSTR = .TRUE. 605C 606C ---------------------------------------------- 607C try to restart from vectors available on file: 608C ---------------------------------------------- 609 IF (CCRSTR) THEN 610 IOPT = 33 611 ILSTNR = ISTRVE + I - 1 612 ISYM = ILSTSYM(LIST,ILSTNR) 613 CALL DZERO(WORK(KV),NCCVAR) 614 CALL CC_RDRSP(LIST,ILSTNR,ISYM,IOPT,MODEL, 615 & WORK(KV),WORK(KV+NT1AM(ISYM))) 616 IF (CCR12 .AND. IOPT.NE.33) THEN 617 IOPTR12 = -32 618 CALL CC_RDRSP(LIST,ILSTNR,ISYM,IOPTR12,MODEL, 619 & DUMMY,WORK(KV+NT1AM(ISYM)+NT2AM(ISYM))) 620 END IF 621 ! check if readin of vector was succesfull - else calc. new. 622 LEN_LIST = LEN(LIST) 623 IF (IOPT .EQ. 33) THEN 624 NEWSTR = .TRUE. 625 WRITE(LUPRI,'(A,I3,3A,I3,A)') ' Vector nr.',ILSTNR, 626 & ' of type ',LIST(1:LEN_LIST),' and symmetry',ISYM, 627 & ' was not found on file.' 628 ELSE 629 NEWSTR = .FALSE. 630 WRITE(LUPRI,'(A,I3,A,I3,A)') ' Vector nr.',ILSTNR, 631 & ' of symmetry',ISYM,' found on file - RESTART SUCCESS' 632 WRITE(LUPRI,'(4A)') ' Start vector is a ', 633 & MODEL,LIST(1:3),' vector' 634Casm & MODEL,LIST(1:LEN_LIST),' vector' 635 ENDIF 636 ENDIF 637C 638C ----------------------------------------------------------- 639C for left eigenvectors try to start from right eigenvectors: 640C ----------------------------------------------------------- 641 IF (NEWSTR .AND. (LIST .EQ. 'LE')) THEN 642 IOPT = 33 643 ILSTNR = ISTRVE + I - 1 644 ISYM = ILSTSYM(LIST,ILSTNR) 645 WRITE(LUPRI,'(2A)') 'Try to use right eigenvector as start', 646 & ' for a left eigenvector.' 647 CALL DZERO(WORK(KV),NCCVAR) 648 CALL CC_RDRSP('RE ',ILSTNR,ISYM,IOPT,MODEL, 649 & WORK(KV),WORK(KV+NT1AM(ISYM))) 650 IF (CCR12 .AND. IOPT.NE.33) THEN 651 IOPTR12 = -32 652 CALL CC_RDRSP('RE ',ILSTNR,ISYM,IOPTR12,MODEL, 653 & DUMMY,WORK(KV+NT1AM(ISYM)+NT2AM(ISYM))) 654 END IF 655 IF (IOPT .EQ. 33) THEN 656 NEWSTR = .TRUE. 657 WRITE(LUPRI,'(1X,A,I3,A,I3,A)') 'Vector nr.',ILSTNR, 658 * ' of type RE and symmetry',ISYM,' was not found on file.' 659 ELSE 660 NEWSTR = .FALSE. 661 WRITE(LUPRI,'(1X,A,I3,A,I3,A)') 'Vector nr.',ILSTNR, 662 * ' of symmetry',ISYM,' found on file - RESTART SUCCESS' 663 WRITE(LUPRI,'(1X,A,A13,A)') 'Start vector is a ', 664 * MODEL//'RE',' vector' 665 ENDIF 666 ENDIF 667C 668C ------------------------------------------------------------ 669C if a restart was not possible generate a trial vector from 670C the gradient vector (for lin.eq.) or by choosing a unit vec. 671C ------------------------------------------------------------ 672 IF (NEWSTR) THEN 673 674 IF (LIST(2:2) .EQ.'E') THEN 675C --------------------------------------------------------- 676C for excited state pick unit vector according to diagonal: 677C --------------------------------------------------------- 678 WRITE(LUPRI,'(1X,A)')'Start vector guessed from diagonal' 679 WRITE(LUPRI,'(1X,A,I3)')'... selected element no.', 680 & IPLACE(I) 681 CALL DZERO(WORK(KV),NCCVAR) 682 WORK(KV-1+IPLACE(I)) = D1 683 ELSE 684C -------------------- 685C get gradient vector: 686C -------------------- 687 ILSTNR = ISTRVE + I - 1 688 CALL CC_GETGD(WORK(KV),ILSTNR,ISIDE,LIST, 689 * WORK(KEND1),LEND1) 690 IF (IPRLE.GT.0) THEN 691 IF (LIST(1:2).EQ.'L0') THEN 692 WRITE(LUPRI,'(5X,A)') 693 * 'Start vector generated from gradient' 694 ELSE 695 WRITE(LUPRI,'(5X,2A,I3,A,I3,A)') LIST, 696 * ' start vector nr.',ILSTNR, 697 * ' of symmetry',ISYMTR,' generated from gradient' 698 END IF 699 END IF 700C ------------------------------------------ 701C Scale with diagonal.(including frequency). 702C (not for Cauchy vectors. CH.H. 20-11-97) 703C ------------------------------------------ 704 IF ( LIST(1:2).NE.'RC' .AND. LIST(1:2).NE.'LC' .AND. 705 * LIST(1:3).NE.'CRC'.AND. LIST(1:3).NE.'CLC' ) THEN 706 KOMEG1 = KV 707 KOMEG2 = KV + NT1AM(ISYMTR) 708 CALL CC_VSCAL(WORK(KOMEG1),WORK(KOMEG2),EIVAL(I), 709 * WORK(KEND1),LEND1,ISYMTR) 710 ENDIF 711 ENDIF 712 713 ENDIF 714C 715 IF ( DEBUG .OR. LOCDBG ) THEN 716 RHO = DDOT(NCCVAR,WORK(KV),1,WORK(KV),1) 717 WRITE(LUPRI,*) 'Norm of Start vector before CCORT: ',RHO 718 ENDIF 719C 720 IF (LOCDBG) THEN 721 CALL AROUND(' CCEQ_STR:VECTORS BEFORE CCORT/PRJDRV: ') 722 WRITE(LUPRI,*) 'this is vector number:',I 723 CALL OUTPUT(WORK(KV),1,NCCVAR,1,1,NCCVAR,1,1,LUPRI) 724 ENDIF 725 726 NFIN = 1 727C 728C------------------------------------------------------------- 729C if requested, project excited reference state out: 730C ------------------------------------------------------------- 731C 732 IF (LPROJECT) THEN 733 IF (TRIPLET) CALL QUIT('LPROJECT & TRIPLET not tested!') 734 IF (DEBUG) THEN 735 WRITE(LUPRI,*)' ---- CCEQ_STR: PROJECTION REQUIRED ----' 736 END IF 737 KPRJC = KEND2 738 ISYMPR = ILSTSYM(LIST,ILSTNR) 739 CALL PRJDRV(ISIDE,ISTATPRJ,ISYMPR,NFIN,NCCVAR, 740 * WORK(KV),WORK(KPRJC)) 741 ENDIF 742C 743C------------------------------------------------------------- 744C if requested, project out valence excitation elements 745C or core excitation elements or virtuals 746C ------------------------------------------------------------- 747C 748 IF (LCVSEXCI.or.LIONIZEXCI) THEN 749 !write(lupri,*)'CCEQ_STR: freeze VALENCE OR VIRT' 750 if (triplet) then 751 KCAM1 = KV 752 KCAMP = KCAM1 + NT1AM(ISYMTR) 753 KCAMM = KCAMP + NT2AM(ISYMTR) 754 CALL CC_FREEZE_TRIPLETEXCI(WORK(KCAM1), 755 & WORK(KCAMP),WORK(KCAMM),ISYMTR, 756 & MAXCORE,MAXION, 757 & NRHFCORE,IRHFCORE, 758 & NVIRION,IVIRION,LBOTHEXCI) 759 else 760 CALL CC_FREEZE_EXCI(WORK(KV),ISYMTR, 761 & MAXCORE,MAXION, 762 & NRHFCORE,IRHFCORE, 763 & NVIRION,IVIRION,LBOTHEXCI) 764 end if 765 END IF 766 IF (LRMCORE) THEN 767 !write(lupri,*)'CCEQ_STR: freeze CORE ' 768 if (TRIPLET) then 769 KCAM1 = KV 770 KCAMP = KCAM1 + NT1AM(ISYMTR) 771 KCAMM = KCAMP + NT2AM(ISYMTR) 772 CALL CC_FREEZE_TRIPLETCORE(WORK(KCAM1), WORK(KCAMP), 773 & WORK(KCAMM),ISYMTR, 774 & MAXCORE,MAXION, 775 & NRHFCORE,IRHFCORE, 776 & NVIRION,IVIRION,LBOTHEXCI) 777 else 778 CALL CC_FREEZE_CORE(WORK(KV),ISYMTR, 779 & MAXCORE, 780 & NRHFCORE,IRHFCORE) 781 end if 782 END IF 783 !SONIA FIXME CHECK IF IT GIVES PROBLEMS WITH CCS 784 IF ((LIST.eq.'L0').and.(LCOR.or.LSEC)) then 785 !write(lupri,*) 'Project out core from start vector' 786 call cc_core(WORK(KV),WORK(KV+NT1AM(ISYMTR)),ISYMTR) 787 end if 788!CN 789! 790! calculate metric for new b-vectors: 791! 792 IF (CCR12) THEN 793 LNOREAD = .TRUE. 794 IF (ISIDE .EQ. 1) THEN 795 CALL CC_R12METRIC(ISYMTR,BRASCL,KETSCL, 796 * WORK(KEND1),LEND1,DUMMY,DUMMY, 797 * DUMMY,DUMMY,FS12AM,LUFS12, 798 * FS2AM,LUFS2,NREDH+1,LNOREAD,WORK(KV)) 799 ELSE IF (ISIDE .EQ. -1) THEN 800 CALL CC_R12METRIC(ISYMTR,0.5D0*KETSCL,2.0D0*BRASCL, 801 * WORK(KEND1),LEND1,DUMMY,DUMMY, 802 * DUMMY,DUMMY,FS12AM,LUFS12, 803 * FS2AM,LUFS2,NREDH+1,LNOREAD,WORK(KV)) 804 END IF 805 END IF 806CCN 807C 808 CALL CCORT(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 809 * LUFS2,FS2AM,LUFS12,FS12AM, 810 * TRIPLET,NFIN,NREDH,NCCVAR, 811 * WORK(KV),ISYMTR,THRLDP,WORK(KEND1),LEND1,IPRLE) 812 NUPVEC = NUPVEC + NFIN 813 IF ( DEBUG .OR. LOCDBG ) THEN 814 RHO = DDOT(NCCVAR,WORK(KV),1,WORK(KV),1) 815 WRITE(LUPRI,*) 'Norm of Start vectors after CCORT:',RHO 816 ENDIF 817C 818 100 CONTINUE 819 IF (IPRLE .GT. 10) 820 * WRITE (LUPRI,'(/A,I5)') ' CCEQ_STR: NREDH =',NREDH 821C 822 IF (IPRLE. GT. 10 ) THEN 823 CALL AROUND( ' END OF CCEQ_STR ' ) 824 ENDIF 825C 826 CALL QEXIT('CCEQ_STR') 827 RETURN 828 END 829C /* Deck ccconv */ 830 SUBROUTINE CCCONV(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 831 * TRIPLET,CCR12,NREDH,LST,LED,SOLEQ,RD,WRK) 832C 833C PURPOSE: 834C Calculate converged solutions to the LED sets of linear 835C equations in RD starting with equation number LEDST 836C 837#include "implicit.h" 838#include "priunit.h" 839#include "cclr.h" 840#include "leinf.h" 841C 842 CHARACTER*(*) C1FIL,C2FIL,C12FIL 843 LOGICAL TRIPLET,CCR12 844C 845 DIMENSION SOLEQ(MAXRED,*),RD(NCCVAR,*),WRK(*) 846C 847 CALL QENTER('CCCONV') 848 CALL DZERO(RD,LED*NCCVAR) 849 DO I = 1, NREDH 850 CALL CC_GETVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 851 * TRIPLET,I,WRK) 852 DO K = 1, LED 853 KVENU = LST -1 + K 854 CALL DAXPY(NCCVAR,SOLEQ(I,KVENU),WRK,1,RD(1,K),1) 855 END DO 856 END DO 857 CALL QEXIT('CCCONV') 858 RETURN 859 END 860C /* Deck ccnex */ 861 SUBROUTINE CCNEX(LIST,ITLE,LPROJECT,ISTATPRJ, 862 * FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 863 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 864 * FS12AM,LUFS12,FS2AM,LUFS2, 865 * LINQCC,TRIPLET,ISIDE,NREDH,ISTRVE, 866 * NVEC,NUPVEC,JCONV,EIVAL,SOLEQ,WRK,LWRK,RMXNORM) 867C 868C PURPOSE: 1) Construct residual (A)*X(I) - EIVAL(I)*X(I) + F 869C for solution X(I) of reduced equations 870C ( F is zero for eigenvalue equation ) 871C 2) Test for convergence of solutions 872C convergence criterium: 873C ||((A)*X(I) - EIVAL(I)*X(I) + F|| / ||X|| .LE. THRLE 874C 3) Use generalized conjugate gradient algorithm 875C or davidson (for eigenvalue equation) 876C to determine next guess of trial vectors 877C 878C JCONV input: if JCONV .lt. 0 do not calculate new trial vectors. 879C output: = 1 converged 880C = 0 not converged 881C = -1 not converged, linear dependency among all 882C trial vectors. 883C 884C RMXNORM : maximum norm of residuals 885C 886#include "implicit.h" 887#include "priunit.h" 888#include "ccsdinp.h" 889#include "ccexci.h" 890#include "ccsdsym.h" 891#include "cclr.h" 892#include "ccfro.h" 893#include "leinf.h" 894#include "r12int.h" 895Cholesky 896#include "maxorb.h" 897#include "ccdeco.h" 898Cholesky 899#include "ccsections.h" 900CVS Sonia 901#include "ccexcicvs.h" 902C 903 LOGICAL LOCDBG,LNOREAD 904 PARAMETER ( LOCDBG = .FALSE. ) 905 PARAMETER ( THRLDP = 1.D-20 ) 906 PARAMETER ( DTEST = 1.0D-04 ) 907 PARAMETER ( DM1=-1.0D0,D1 =1.0D0, D0=0.0D0 ) 908C 909 CHARACTER*(*) FRHO1,FRHO2,FRHO12,FC1AM,FC2AM,FC12AM,FS12AM,FS2AM 910 CHARACTER*(*) LIST 911 INTEGER :: NREDHOLD 912 LOGICAL LINQCC,CCRSTRS,LPROJECT,TRIPLET 913 DIMENSION EIVAL(*), SOLEQ(MAXRED,*),WRK(*) 914C 915C Space for CCNEX: 916C 917C MAXNEX: Maximum number of simultaneous vectors in CCNEX 918C 919 IF (DEBUG) CALL AROUND(' Start of CCNEX ') 920 921 MAXNEX = (LWRK-3*NCCVAR)/NCCVAR 922 NUPVEC = 0 923 NOTCON = 0 924 !RF : Note that NREDH is modified in the following loop 925 ! in CCORT, but we need to know the previous size! 926 NREDHOLD = NREDH 927 928 DO 4000 ISIMC = 1,NVEC,MAXNEX 929 NBX = MIN(MAXNEX,(NVEC+1-ISIMC)) 930c NBX = MIN(MAXNEX,(NVEC+1-ISIMC),(NREDH+1-ISIMC)) 931C 932C CONSTRUCT RESIDUAL IN WRK(KRES) 933C RESIDUAL: (A)*X(I)-EIVAL(I)*X(I) +F 934C 935 KRES = 1 936 KTST = KRES + NBX*NCCVAR 937 KWRK1 = KTST + NBX 938 LWRK1 = LWRK - KWRK1 939 940 IF (LINQCC) THEN 941 LGD = KRES 942 DO JR = 1,NBX 943 944C --------------------- 945C Get gradient vectors: 946C --------------------- 947 IVENU = ISIMC-1+JR 948 ILSTNR = IVENU + ISTRVE - 1 949 950 CALL CC_GETGD(WRK(LGD),ILSTNR,ISIDE,LIST, 951 * WRK(KWRK1),LWRK1) 952C 953C ---------------------- 954C Project from gradient: 955C ---------------------- 956 IF (LPROJECT) THEN 957 ISYPR = ILSTSYM(LIST,ILSTNR) 958 IF (DEBUG) WRITE(LUPRI,*) 959 & '--- CCNEX projection from gradient --- ' 960 CALL PRJDRV(ISIDE,ISTATPRJ,ISYPR,1,NCCVAR, 961 & WRK(LGD),WRK(KWRK1)) 962 END IF 963 IF (LCVSEXCI.or.LIONIZEXCI) THEN 964 if (triplet) then 965 KCAM1 = LGD 966 KCAMP = KCAM1 + NT1AM(ISYMTR) 967 KCAMM = KCAMP + NT2AM(ISYMTR) 968 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1), 969 & WRK(KCAMP), WRK(KCAMM),ISYMTR, 970 & MAXCORE,MAXION, 971 & NRHFCORE,IRHFCORE, 972 & NVIRION,IVIRION,LBOTHEXCI) 973 else 974 !write(lupri,*)'CCNEX:CVS or IONISATION (1)' 975 CALL CC_FREEZE_EXCI(WRK(LGD),ISYMTR, 976 & MAXCORE,MAXION, 977 & NRHFCORE,IRHFCORE, 978 & NVIRION,IVIRION,LBOTHEXCI) 979 end if 980 END IF 981 IF (LRMCORE) THEN 982 !write(lupri,*)'CCNEX: freeze the core' 983 IF (TRIPLET) THEN 984! Eirik & Sonia 985 KCAM1 = LGD 986 KCAMP = KCAM1 + NT1AM(ISYMTR) 987 KCAMM = KCAMP + NT2AM(ISYMTR) 988 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 989 & WRK(KCAMM),ISYMTR, 990 & MAXCORE,MAXION, 991 & NRHFCORE,IRHFCORE, 992 & NVIRION,IVIRION,.false.) 993 ELSE 994 CALL CC_FREEZE_CORE(WRK(LGD),ISYMTR, 995 & MAXCORE, 996 & NRHFCORE,IRHFCORE) 997 END IF 998 END IF 999!Sonia: fixme check for CCS 1000 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1001 write(lupri,*)'REMOVE CORE FROM L0-RHS (CCNEX)' 1002 CALL CC_CORE(WRK(LGD),WRK(LGD+NT1AM(ISYMTR)),isymtr) 1003 ENDIF 1004 ! 1005 LGD = LGD + NCCVAR 1006 END DO 1007 ELSE 1008 CALL DZERO(WRK(KRES),NBX*NCCVAR) 1009 END IF 1010 1011 IF (IPRLE.GT.110 .OR. LOCDBG) THEN 1012 WRITE(LUPRI,*) 'CCNEX> RESIDUAL AFTER GD CONTRIBUTION' 1013 CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NBX,NCCVAR,NBX,1,LUPRI) 1014 END IF 1015C 1016C ----------------------------------------------------------- 1017C add -EIVAL(I)*X(I), where X(I) is the i'th solution vector: 1018C ----------------------------------------------------------- 1019 ISYM = ILSTSYM(LIST,ISTRVE) 1020 NVARPT = NCCVAR + 2*NALLAI(ISYM) 1021 MAXBVE = ( LWRK1 - NVARPT )/NCCVAR 1022 DO 60 JR = 1,NBX,MAXBVE 1023 NBVEC = MIN(MAXBVE,NBX+1-JR) 1024 KB = KWRK1 1025 KWRK2 = KB + NCCVAR*NBVEC 1026 LWRK2 = LWRK - KWRK2 1027 IF (LWRK2 .LT. 0 ) THEN 1028 CALL QUIT('Too little work in CCNEX xx') 1029 END IF 1030C ---------------------------------------------------------- 1031C Construct the solution vectors in full space: 1032C Save norm of solution vector in WRK(KTST) and the 1033C vectors itself for restart on the files CC//LIST//{ivec} 1034C ---------------------------------------------------------- 1035 IVECNU = ISIMC -1 + JR 1036 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM,TRIPLET, 1037 * CCR12,NREDHOLD,IVECNU,NBVEC,SOLEQ,WRK(KB), 1038 & WRK(KWRK2)) 1039 1040 DO IBVEC = 1, NBVEC 1041 WRK(KTST+JR-1+IBVEC-1) = 1042 * DNRM2(NCCVAR,WRK(KB+(IBVEC-1)*NCCVAR),1) 1043 IVEC = ISIMC - 1 + JR -1 + IBVEC 1044 CALL CC_SAVE(WRK(KB+(IBVEC-1)*NCCVAR),IVEC+ISTRVE-1, 1045 * LIST,WRK(KWRK2),LWRK2) 1046 END DO 1047 1048 IF (IPRLE.GT.105 .OR. LOCDBG) THEN 1049 WRITE (LUPRI,'(/A)')' CCNEX: solution vectors' 1050 CALL OUTPUT(WRK(KB),1,NCCVAR,1,NBVEC,NCCVAR,NBVEC,1,LUPRI) 1051 CALL PROVLP(WRK(KB),WRK(KB),NCCVAR,NBVEC,WRK(KWRK2),LUPRI) 1052 END IF 1053 1054C ---------------------------------------------------------- 1055C Construct S * X in full space and add - Eval * S * X to 1056C the residual vectors: 1057C Note for conventional CC the metric S is a unit matrix 1058C and thus S * X are the solution vectors already 1059C constructed above. Only for CCR12 we need to construct 1060C S * X from the transformations with the metric matrix 1061C ---------------------------------------------------------- 1062 IF (CCR12) THEN 1063 IVECNU = ISIMC -1 + JR 1064 IF (IANR12.EQ.1 .OR. IANR12.EQ.3) THEN 1065 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFS12,FS12AM, 1066 * TRIPLET,CCR12,NREDHOLD,IVECNU,NBVEC,SOLEQ, 1067 * WRK(KB),WRK(KWRK2)) 1068 ELSE IF (IANR12.EQ.2) THEN 1069 CALL CCCONV(LUFC1,FC1AM,LUFS2,FS2AM,LUFS12,FS12AM, 1070 * TRIPLET,CCR12,NREDHOLD,IVECNU,NBVEC,SOLEQ, 1071 * WRK(KB),WRK(KWRK2)) 1072 ELSE 1073 CALL QUIT('This R12 ansatz is not implemented yet.') 1074 END IF 1075 END IF 1076 1077 DO IBVEC = 1,NBVEC 1078 IVEC = ISIMC - 1 + JR -1 + IBVEC 1079 IVEOFF = JR -1 + IBVEC -1 1080 CALL DAXPY(NCCVAR,-EIVAL(IVEC),WRK(KB+(IBVEC-1)*NCCVAR), 1081 * 1,WRK(KRES+IVEOFF*NCCVAR),1) 1082 END DO 1083C 1084 60 CONTINUE 1085 1086 IF (IPRLE.GT.110 .OR. LOCDBG) THEN 1087 WRITE(LUPRI,*)' RESIDUAL AFTER - w * S * X CONTRIBUTION' 1088 CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NBX,NCCVAR,NBX,1,LUPRI) 1089 END IF 1090 1091C ----------------------------------------------- 1092C Add A * X(I) where X(I) is the solution to the 1093C I'th set of Newton-Raphson equations 1094C ----------------------------------------------- 1095 DO K = 1,NREDHOLD 1096 CALL CC_GETVEC(LUFR1,FRHO1,LUFR2,FRHO2,LUFR12,FRHO12, 1097 * TRIPLET,K,WRK(KWRK1)) 1098 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1099 if (triplet) then 1100 KCAM1 = KWRK1 1101 KCAMP = KCAM1 + NT1AM(ISYMTR) 1102 KCAMM = KCAMP + NT2AM(ISYMTR) 1103 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1), 1104 & WRK(KCAMP),WRK(KCAMM),ISYMTR, 1105 & MAXCORE,MAXION, 1106 & NRHFCORE,IRHFCORE, 1107 & NVIRION,IVIRION,LBOTHEXCI) 1108 else 1109 CALL CC_FREEZE_EXCI(WRK(KWRK1),ISYMTR, 1110 & MAXCORE,MAXION, 1111 & NRHFCORE,IRHFCORE, 1112 & NVIRION,IVIRION,LBOTHEXCI) 1113 end if 1114 END IF 1115 IF (LRMCORE) THEN 1116 IF (TRIPLET) THEN 1117C Eirik & Sonia 1118 KCAM1 = KWRK1 1119 KCAMP = KCAM1 + NT1AM(ISYMTR) 1120 KCAMM = KCAMP + NT2AM(ISYMTR) 1121 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 1122 & WRK(KCAMM),ISYMTR, 1123 & MAXCORE,MAXION, 1124 & NRHFCORE,IRHFCORE, 1125 & NVIRION,IVIRION,.false.) 1126 ELSE 1127!write(lupri,*)'CCNEX: freeze the core' 1128 CALL CC_FREEZE_CORE(WRK(KWRK1),ISYMTR, 1129 & MAXCORE, 1130 & NRHFCORE,IRHFCORE) 1131 END IF 1132 END IF 1133 !Sonia: fixme check for CCS 1134 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1135 write(lupri,*)'REMOVE CORE FROM L0 in CCNEX' 1136 CALL CC_CORE(WRK(KWRK1), 1137 & WRK(KWRK1+NT1AM(ISYMTR)),isymtr) 1138 ENDIF 1139! 1140 DO JR = 1,NBX 1141 JROOTJ = ISIMC - 1 + JR 1142 EVAL1 = SOLEQ(K,JROOTJ) 1143 CALL DAXPY(NCCVAR,EVAL1,WRK(KWRK1),1, 1144 * WRK(KRES+(JR-1)*NCCVAR),1) 1145 END DO 1146 END DO 1147 1148C ------------------------------- 1149C Residual is now in WRK(KRES) 1150C ------------------------------- 1151 IF (IPRLE.GT.45 .OR. LOCDBG) THEN 1152 WRITE (LUPRI,'(/A)') ' CCNEX: residual vectors ' 1153 CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NBX,NCCVAR,NBX,1,LUPRI) 1154 CALL PROVLP(WRK(KRES),WRK(KRES),NCCVAR,NBX,WRK(KWRK1),LUPRI) 1155 END IF 1156 1157C ----------------------------------------- 1158C Projection of eigenvectors from residual: 1159C ----------------------------------------- 1160 IF (DEBUG) WRITE (LUPRI,*) 'CCNEX> LPROJECT = ', LPROJECT 1161 IF (LPROJECT) THEN 1162 KPRJ = KWRK2 1163 ISYMPR = ILSTSYM(LIST,ISTRVE) 1164 IF (DEBUG) THEN 1165 WRITE(LUPRI,*)'-- CCNEX projection from residual -- ' 1166 WRITE(LUPRI,*)' NBX IN CCNEX',NBX 1167 END IF 1168 CALL PRJDRV(ISIDE,ISTATPRJ,ISYMPR,NBX,NCCVAR, 1169 * WRK(KRES),WRK(KWRK2)) 1170 ENDIF 1171C ------------------------------------------------------- 1172C Remove specific excitations from residual if requested: 1173C ------------------------------------------------------- 1174 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1175 DO JR = 1,NBX 1176 !WRITE(LUPRI,*)'CCNEX CVS or IONISATION (resid)' 1177 KOFF = KRES+(JR-1)*NCCVAR 1178 if (triplet) then 1179 KCAM1 = Koff 1180 KCAMP = KCAM1 + NT1AM(ISYMTR) 1181 KCAMM = KCAMP + NT2AM(ISYMTR) 1182 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1), WRK(KCAMP), 1183 & WRK(KCAMM),ISYMTR, 1184 & MAXCORE,MAXION, 1185 & NRHFCORE,IRHFCORE, 1186 & NVIRION,IVIRION,LBOTHEXCI) 1187 else 1188 CALL CC_FREEZE_EXCI(WRK(Koff),ISYMTR, 1189 & MAXCORE,MAXION, 1190 & NRHFCORE,IRHFCORE, 1191 & NVIRION,IVIRION,LBOTHEXCI) 1192 end if 1193 END DO 1194 END IF 1195 IF (LRMCORE) THEN 1196 DO JR = 1,NBX 1197 !WRITE(LUPRI,*)'freeze core in resid' 1198 KOFF = KRES+(JR-1)*NCCVAR 1199 IF (TRIPLET) THEN 1200C Eirik 1201 KCAM1 = KOFF 1202 KCAMP = KCAM1 + NT1AM(ISYMTR) 1203 KCAMM = KCAMP + NT2AM(ISYMTR) 1204 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 1205 & WRK(KCAMM),ISYMTR, 1206 & MAXCORE,MAXION, 1207 & NRHFCORE,IRHFCORE, 1208 & NVIRION,IVIRION,.FALSE.) 1209 ELSE 1210 CALL CC_FREEZE_CORE(WRK(Koff),ISYMTR, 1211 & MAXCORE, 1212 & NRHFCORE,IRHFCORE) 1213 END IF 1214 END DO 1215 END IF 1216 !Sonia, fixme CCS!!!! 1217 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1218 DO JR = 1,NBX 1219 KOFF = KRES+(JR-1)*NCCVAR 1220 write(lupri,*)'SONIA: REMOVE CORE FROM L0 in CCNEX' 1221 CALL CC_CORE(WRK(Koff), 1222 & WRK(Koff+NT1AM(ISYMTR)),isymtr) 1223 end do 1224 ENDIF 1225C ------------------------------------------ 1226C Test for convergence of the N-R solutions: 1227C ------------------------------------------ 1228 RMXNORM = 0.0D0 1229 NFIN = 0 1230 IF (IPRLE .GT. 1) WRITE(LUPRI,*) ' ' 1231 DO 2000 JR = 1,NBX 1232 JROOTJ = ISIMC - 1 + JR 1233 QNORM = DNRM2(NCCVAR,WRK(KRES+(JR-1)*NCCVAR),1) 1234 * /MAX(WRK(KTST+JR-1),1.0D-14) 1235 RMXNORM = MAX(RMXNORM,QNORM) 1236C 1237 IF (IPRLE .GT. 1 .OR. DEBUG) THEN 1238 WRITE(LUPRI,'(A,I3,1P,1(A,D12.5))') 1239 * ' ROOT:',JROOTJ,' RESIDUAL TOT:',QNORM 1240 END IF 1241 1242 IF (QNORM.GT.THRLE) THEN 1243 NFIN = NFIN + 1 1244 IF ( NFIN.NE.JR) THEN 1245 CALL DCOPY(NCCVAR,WRK(KRES+(JR-1)*NCCVAR), 1246 * 1,WRK(KRES+(NFIN-1)*NCCVAR),1) 1247 END IF 1248 ! Overwrite KTST with frequency of this vector 1249 WRK(KTST+NFIN-1) = EIVAL(JROOTJ) 1250 ELSE 1251 IF (IPRLE .GT. 1) THEN 1252 WRITE(LUPRI,'(A,I3,1P,1(A,D12.5))') 1253 * ' ROOT:',JROOTJ,' HAS CONVERGED. RESIDUAL TOT:',QNORM 1254 END IF 1255 END IF 1256 2000 CONTINUE 1257 NOTCON = NOTCON + NFIN 1258C 1259 IF (NFIN.EQ.0 .OR. JCONV.LT.0) GO TO 3999 1260C 1261C -------------------------------------------- 1262C Use generalized conjugate gradient algorithm 1263C to form new trial vectors 1264C -------------------------------------------- 1265Cholesky 1266C 1267C Cholesky CC2 did not converge using precond, so 1268C uncomment next lines and go back to old safe method 1269C 1270 IF (CHOINT.OR.CCSLV) THEN 1271 KDIA = KWRK1 1272 KWRK2 = KDIA + NCCVAR 1273 LWRK2 = LWRK - KWRK2 1274 IF (LWRK2 .LT. 0 ) THEN 1275 CALL QUIT('Too little work in CCNEX') 1276 END IF 1277 CALL CCLR_DIA(WRK(KDIA),ISYMTR,TRIPLET,WRK(KWRK2),LWRK2) 1278 CALL LEDIAG(NCCVAR,WRK(KDIA),DTEST) 1279C 1280 DO JR = 1,NFIN 1281 IOFF = (JR-1)*NCCVAR 1282 DO I=1,NCCVAR 1283 WRK(IOFF+I) = WRK(IOFF+I)*WRK(KDIA-1+I) 1284 END DO 1285 END DO 1286C 1287 ELSE 1288CCH 1289C WRITE(LUPRI,'(/I4,1X,A)') NFIN, 1290C * 'new trial vectors before PRECOND:' 1291C WRITE (LUPRI,*) 'NORM: ',DDOT(NCCVAR*NFIN,WRK(KRES),1, 1292C * WRK(KRES),1) 1293C CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NFIN,NCCVAR,NFIN,1,LUPRI) 1294CCH 1295 CALL CC_PRECOND(WRK(KRES),WRK(KTST), 1296 & NFIN,NCCVAR,ISYMTR,TRIPLET, 1297 & WRK(KWRK1),LWRK1) 1298 END IF 1299C 1300Cholesky 1301C 1302 IF ((IPRLE.GT.105 .OR. LOCDBG).AND. NFIN .GT. 0) THEN 1303 WRITE(LUPRI,'(/I4,1X,A)') NFIN, 1304 * 'new trial vectors before CCORT/PRJDRV:' 1305 CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NFIN,NCCVAR,NFIN,1,LUPRI) 1306 CALL PROVLP(WRK(KRES),WRK(KRES),NCCVAR,NFIN,WRK(KWRK1),LUPRI) 1307 CALL FLSHFO(LUPRI) 1308 END IF 1309 1310C ------------------------------ 1311C Projection from updated trials 1312C ------------------------------ 1313 IF (LPROJECT) THEN 1314 KPRJ = KWRK2 1315 ISYMPR = ILSTSYM(LIST,ISTRVE) 1316 IF (DEBUG) THEN 1317 WRITE(LUPRI,*)' NFIN IN CCNEX',NFIN 1318 WRITE(LUPRI,*) 1319 & ' ---- CCNEX Before projection on trials ----' 1320 END IF 1321 CALL PRJDRV(ISIDE,ISTATPRJ,ISYMPR,NFIN,NCCVAR, 1322 * WRK(KRES),WRK(KPRJ)) 1323 ENDIF 1324C ----------------------------------------------------------- 1325C Projection specific orbital excitations from updated trials 1326C ----------------------------------------------------------- 1327 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1328 DO JR = 1,NFIN 1329 !write(lupri,*)'CCNEX CVS or IONISATION upd trials' 1330 koff = KRES+(JR-1)*NCCVAR 1331 if (triplet) then 1332 KCAM1 = Koff 1333 KCAMP = KCAM1 + NT1AM(ISYMTR) 1334 KCAMM = KCAMP + NT2AM(ISYMTR) 1335 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1),WRK(KCAMP), 1336 & WRK(KCAMM),ISYMTR, 1337 & MAXCORE,MAXION, 1338 & NRHFCORE,IRHFCORE, 1339 & NVIRION,IVIRION,LBOTHEXCI) 1340 else 1341 CALL CC_FREEZE_EXCI(WRK(Koff),ISYMTR, 1342 & MAXCORE,MAXION, 1343 & NRHFCORE,IRHFCORE, 1344 & NVIRION,IVIRION,LBOTHEXCI) 1345 end if 1346 END DO 1347 END IF 1348 IF (LRMCORE) THEN 1349 DO JR = 1,NFIN 1350 !write(lupri,*)'CCNEX CVS or IONISATION upd trials' 1351 koff = KRES+(JR-1)*NCCVAR 1352 IF (TRIPLET) THEN 1353C Eirik 1354 KCAM1 = KOFF 1355 KCAMP = KCAM1 + NT1AM(ISYMTR) 1356 KCAMM = KCAMP + NT2AM(ISYMTR) 1357 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 1358 & WRK(KCAMM),ISYMTR, 1359 & MAXCORE,MAXION, 1360 & NRHFCORE,IRHFCORE, 1361 & NVIRION,IVIRION,.false.) 1362 ELSE 1363 CALL CC_FREEZE_CORE(WRK(Koff),ISYMTR, 1364 & MAXCORE,NRHFCORE,IRHFCORE) 1365 END IF 1366 END DO 1367 END IF 1368 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1369 do JR = 1,NFIN 1370 koff = KRES+(JR-1)*NCCVAR 1371 write(lupri,*)'SONIA REMOVE CORE FROM L0 in CCNEX' 1372 !Sonia FIXME for CCS 1373 CALL CC_CORE(WRK(Koff),WRK(Koff+NT1AM(ISYMTR)),isymtr) 1374 end do 1375 ENDIF 1376CCN 1377C ----------------------------------- 1378C calculate metric for trial vectors: 1379C ----------------------------------- 1380 IF (CCR12) THEN 1381 DO I = 1, NFIN 1382 LNOREAD = .TRUE. 1383 IF (ISIDE .EQ. 1) THEN 1384 CALL CC_R12METRIC(ISYMTR,BRASCL,KETSCL, 1385 * WRK(KWRK1),LWRK1,DUMMY,DUMMY, 1386 * DUMMY,DUMMY,FS12AM,LUFS12, 1387 * FS2AM,LUFS2,NREDH+I,LNOREAD, 1388 * WRK(KRES+NCCVAR*(I-1))) 1389 ELSE IF (ISIDE .EQ. -1) THEN 1390 CALL CC_R12METRIC(ISYMTR,0.5D0*KETSCL,2.0D0*BRASCL, 1391 * WRK(KWRK1),LWRK1,DUMMY,DUMMY, 1392 * DUMMY,DUMMY,FS12AM,LUFS12, 1393 * FS2AM,LUFS2,NREDH+I,LNOREAD, 1394 * WRK(KRES+NCCVAR*(I-1))) 1395 END IF 1396 END DO 1397 END IF 1398CCN 1399C -------------------------------------------------------------- 1400C Orthogonalize trial vectors and examine for linear dependence: 1401C -------------------------------------------------------------- 1402 CALL CCORT(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 1403 * LUFS2,FS2AM,LUFS12,FS12AM, 1404 * TRIPLET,NFIN,NREDH,NCCVAR, 1405 * WRK(KRES),ISYMTR,THRLDP,WRK(KWRK1),LWRK1,IPRLE) 1406C 1407 IF ((IPRLE.GT.105.OR.LOCDBG) .AND. NFIN.GT.0) THEN 1408 WRITE (LUPRI,'(/I4,1X,A)') NFIN, 1409 * 'new trial vectors after CCORT/PRJDRV:' 1410 CALL OUTPUT(WRK(KRES),1,NCCVAR,1,NFIN,NCCVAR,NFIN,1,LUPRI) 1411 CALL PROVLP(WRK(KRES),WRK(KRES),NCCVAR,NFIN,WRK(KWRK1),LUPRI) 1412 END IF 1413 1414 3999 CONTINUE 1415 NUPVEC = NUPVEC + NFIN 1416 4000 CONTINUE 1417 1418C we have now: 1419C NUPVEC : number of new trial vectors 1420C NOTCON : number of equations not yet converged 1421C 1422 IF (NOTCON .EQ. 0 .AND. NUPVEC.EQ.0) THEN 1423 ! ALL EQUATIONS HAVE CONVERGED 1424 IF (IPRLE.GT.10) WRITE(LUPRI,'(A)') ' *** EQUATIONS CONVERGED' 1425 JCONV = 1 1426 RETURN 1427 END IF 1428C 1429C Do NOT calculate new trial vectors 1430C 1431 IF (JCONV .LT. 0) THEN 1432 JCONV = 0 1433 RETURN 1434 END IF 1435C 1436C Linear dependence between trial vectors 1437C 1438 IF (NUPVEC .EQ. 0 .AND. NOTCON.GT.0) THEN 1439 WRITE(LUPRI,5010) 1440 WRITE(LUPRI,*)(NVEC-NOTCON),' trial vectors converged' 1441 WRITE(LUPRI,*)' LINEAR DEPENDENCE BETWEEN',NOTCON, 1442 * ' NON CONVERGED TRIAL VECTORS' 1443 JCONV = -1 1444 RETURN 1445 ENDIF 1446 5010 FORMAT(/' *** MICROITERATIONS STOPPED DUE TO LINEAR ', 1447 * 'DEPENDENCE BETWEEN ALL NEW TRIAL VECTORS') 1448C 1449C----------------------------------------------------------- 1450C Ove Christiansen 2-4-1996: 1451C Quick and dirty way of skipping all vectors up to this 1452C point and restart from the optimal solution so far. 1453C----------------------------------------------------------- 1454C 1455 IF ( NREDH .GE. MXLRV ) THEN 1456C IF ( IPRLE .GT. 20 ) THEN 1457 WRITE(LUPRI,*) ' NREDH : ',NREDH 1458 WRITE(LUPRI,*) ' MXLRV : ',MXLRV 1459C ENDIF 1460C 1461 WRITE(LUPRI,'(/,1X,A)') 1462 * 'ATTENTION!!!! ' 1463 WRITE(LUPRI,'(/,1X,A,/,A)') 1464 * 'All trial and transformed vectors are skipped '// 1465 * 'and the iterative procedure ', 1466 * ' is continued from the current optimal solution ' 1467 WRITE(LUPRI,'(/1X,A,/)') 1468 * 'Notice: You asked for it by setting a low MXLRV.' 1469C 1470 IF(ITLE.LE.1) THEN 1471 WRITE(LUPRI,*)' MXLRV:',MXLRV, 1472 & ' ALLOWS ONE ITERATION ,ITLE:',ITLE 1473 WRITE(LUPRI,*)' NEW INFORMATION NOT CONTAINED IN' 1474 WRITE(LUPRI,*)' NEW REDUCED SPACE, SPECIFY NEW MXLRV' 1475 CALL QUIT('CCNEX SPECIFY NEW MXLRV') 1476 END IF 1477 CCRSTRS = CCRSTR 1478 CCRSTR = .TRUE. 1479 KIPLAC = 1 1480 KWRK1 = KIPLAC + MAXRED 1481 LWRK1 = LWRK - KWRK1 1482 CALL CCEQ_STR(FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 1483 * FS12AM,LUFS12,FS2AM,LUFS2,LPROJECT,ISTATPRJ, 1484 * TRIPLET,ISIDE,ISTRVE,NVEC,NUPVEC, 1485 * NREDH,EIVAL,WRK(KIPLAC), 1486 * WRK(KWRK1),LWRK1,LIST) 1487 CCRSTR = CCRSTRS 1488 IF ( IPRLE .GT. 20 ) THEN 1489 WRITE(LUPRI,*) ' NREDH : ', NREDH 1490 WRITE(LUPRI,*) ' NUPVEC: ', NUPVEC 1491 WRITE(LUPRI,*) ' NVEC: ', NVEC 1492 WRITE(LUPRI,*) ' MAXRED: ', MAXRED 1493 WRITE(LUPRI,*) ' MXLRV : ', MXLRV 1494 ENDIF 1495 CALL FLSHFO(LUPRI) 1496C 1497 ENDIF 1498C 1499 IF (DEBUG) THEN 1500 CALL AROUND(' End of CCNEX ') 1501 CALL FLSHFO(LUPRI) 1502 ENDIF 1503C 1504C End of CCNEX 1505C 1506 RETURN 1507 END 1508C /* Deck ccred */ 1509 SUBROUTINE CCRED(FRHO1,LUFR1,FRHO2,LUFR2,FRHO12,LUFR12, 1510 * FC1AM,LUFC1,FC2AM,LUFC2,FC12AM,LUFC12, 1511 * CCR12,FS12AM,LUFS12,FS2AM,LUFS2, 1512 * LINQCC,TRIPLET,ISIDE,LIST, 1513 * LPROJECT,ISTATPRJ,LREDS,REDS, 1514 * REDH,NREDH,ISTRVE,NUPVEC,NVEC, 1515 * EIVAL,SOLEQ,WRK,LWRK,DEBUG) 1516C 1517C 1518C Input: 1519C ISIDE, 1 for right transformation ,-1 for left 1520C NREDH, dimension of the (updated) Hessian/Jacobian matrix in the 1521C reduced space (NREDH is not changed in CCRED) 1522C ISTRVE is startnumber on list of equations to be solved. 1523C NUPVEC, number of new b-vectors for which REDH is extended 1524C NVEC, number of linear equations or roots to be solved for 1525C in the reduced space 1526C 1527C LINQCC.EQ.T : solve nvec set of linear equations 1528C EIVAL contains NVEC frequencies 1529C EIVAL( ISTRVE ... (ISTRVE+NVEC-1) ) 1530C 1531C LINQCC.EQ.F : solve for NVEC lowest eigenvalues 1532C 1533C (The reduced PROJECTED HESSIAN matrix is the projection 1534C on the basis of b-vectors. 1535C REDH(I,J) = B(I) * S(J) ) 1536C 1537C Output: 1538C REDH, Jacobian in reduced space (dimension: NREDH) 1539C REDS, overlap in reduced space (dimension: NREDH) 1540C SOLEQ, solutions to the NVEC set of NEWTON equations 1541C or eigenvalue equations 1542C EIVAL, eigenvalues or frequencies 1543C 1544C Flow: 1545C 1. if NUPVEC gt 0 update reduced jacobian (and metric) with nupvec 1546C vectors 1547C 2. determine NVEC solution vectors, returned in SOLEQ for eigenvalue 1548C equation. The reduced eigenvalues are returned in EIVAL 1549C 1550#include "implicit.h" 1551#include "priunit.h" 1552#include "cclr.h" 1553#include "ccfro.h" 1554#include "leinf.h" 1555#include "ccrc1rsp.h" 1556#include "r12int.h" 1557#include "ccsdsym.h" 1558CSonia: orbital excitations projection 1559#include "ccexcicvs.h" 1560C 1561 PARAMETER (D1 = 1.0D0 , D0 = 0.0D0, THRZER = 1.0D-99 ) 1562C 1563 CHARACTER*(*) FRHO1,FRHO2,FC1AM,FC2AM,FRHO12,FC12AM,FS12AM 1564 CHARACTER*(*) LIST,FS2AM 1565 LOGICAL LINQCC, LPROJECT, TRIPLET, DEBUG, CCR12, LREDS, LOCDBG 1566 PARAMETER (LOCDBG = .FALSE.) 1567 DIMENSION REDH(*),SOLEQ(MAXRED,*),REDS(*) 1568 DIMENSION EIVAL(*),WRK(*) 1569C 1570C ************************************************************** 1571C Section 1: 1572C extend reduced PROJECTED HESSIAN matrix with NUPVEC new b-vectors 1573C ************************************************************** 1574C 1575 IF (IPRLE.GT.5) THEN 1576 WRITE(LUPRI,*)' CCRED ' 1577 WRITE(LUPRI,*)' ISIDE ',ISIDE 1578 WRITE(LUPRI,*)' LREDS ',LREDS 1579 WRITE(LUPRI,*)' NREDH ',NREDH 1580 WRITE(LUPRI,*)' ISTRVE ',ISTRVE 1581 WRITE(LUPRI,*)' NUPVEC ',NUPVEC 1582 WRITE(LUPRI,*)' NVEC ',NVEC 1583 WRITE(LUPRI,*)' LWRK ',LWRK 1584 WRITE(LUPRI,*)' NCCVAR ', NCCVAR 1585 WRITE(LUPRI,*)' CVSEPARA', LCVSEXCI 1586 WRITE(LUPRI,*)' IONISATI', LIONIZEXCI 1587 WRITE(LUPRI,*)' REMOVE_CORE', LRMCORE 1588 WRITE(LUPRI,*)' LCOR,LSEC in GS', LCOR, LSEC 1589 END IF 1590C 1591 IF (NUPVEC.GT.0) THEN 1592 IF (LREDS) THEN 1593 MAXVEC = (LWRK-NCCVAR)/(2*NCCVAR) 1594 ELSE 1595 MAXVEC = (LWRK-NCCVAR)/NCCVAR 1596 END IF 1597 1598 IREDH = NREDH - NUPVEC 1599 DO 100 IVEC = 1,NUPVEC,MAXVEC 1600 NSIM = MIN(MAXVEC,NUPVEC+1-IVEC) 1601C 1602C work space allocation 1603C 1604 KSVEC = 1 1605 KBVEC = KSVEC + NSIM*NCCVAR 1606 KSBVEC = KBVEC + NCCVAR 1607C 1608C read in sigma vectors in batches 1609C 1610 LSVEC = KSVEC 1611 DO INUM = 1,NSIM 1612 IVENU = IREDH + IVEC-1 + INUM 1613 CALL CC_GETVEC(LUFR1,FRHO1,LUFR2,FRHO2,LUFR12,FRHO12, 1614 * TRIPLET,IVENU,WRK(LSVEC)) 1615!sonia 1616 !freeze valence or virtual 1617 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1618 if (triplet) then 1619 KCAM1 = LSVEC 1620 KCAMP = KCAM1 + NT1AM(ISYMTR) 1621 KCAMM = KCAMP + NT2AM(ISYMTR) 1622 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1), 1623 & WRK(KCAMP),WRK(KCAMM),ISYMTR, 1624 & MAXCORE,MAXION, 1625 & NRHFCORE,IRHFCORE, 1626 & NVIRION,IVIRION,LBOTHEXCI) 1627 else 1628 write(lupri,*)'CCRED CVS/IONIS sigma_',ivenu, 1629 & 'project out valence occupied or virtual' 1630 CALL CC_FREEZE_EXCI(WRK(LSVEC),ISYMTR, 1631 & MAXCORE,MAXION, 1632 & NRHFCORE,IRHFCORE, 1633 & NVIRION,IVIRION,LBOTHEXCI) 1634 end if 1635 END IF 1636 !freeze core in exci calculation 1637 IF (LRMCORE) THEN 1638 write(lupri,*)'CCRED CVS/IONIS sigma_',ivenu, 1639 & 'project out core or virtual' 1640 IF (TRIPLET) THEN 1641C Eirik 1642 KCAM1 = LSVEC 1643 KCAMP = KCAM1 + NT1AM(ISYMTR) 1644 KCAMM = KCAMP + NT2AM(ISYMTR) 1645 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 1646 & WRK(KCAMM),ISYMTR, 1647 & MAXCORE,MAXION, 1648 & NRHFCORE,IRHFCORE, 1649 & NVIRION,IVIRION,.false.) 1650 ELSE 1651 CALL CC_FREEZE_CORE(WRK(LSVEC),ISYMTR, 1652 & MAXCORE, 1653 & NRHFCORE,IRHFCORE) 1654 END IF 1655 END IF 1656!Sonia: FIXME for CCS 1657 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1658 write(lupri,*)'SONIA REMOVE CORE FROM L0 in CCRED' 1659 CALL CC_CORE(WRK(LSVEC),WRK(LSVEC+NT1AM(ISYMTR)), 1660 & isymtr) 1661 ENDIF 1662!sonia 1663 LSVEC = LSVEC + NCCVAR 1664 END DO 1665 IF (LOCDBG) THEN 1666 WRITE(LUPRI,'(a,i3,a,i3)') 1667 * 'A * C for ',IVEC,' - ',IVEC+NSIM-1 1668 CALL OUTPUT(WRK(KSVEC),1,NCCVAR,1,NSIM, 1669 * NCCVAR,NSIM,1,LUPRI) 1670 END IF 1671 1672 IF (LREDS) THEN 1673 ! read in S * b vectors in batches, assumes that the 1674 ! conventional CC amplitudes are in the orthogonal basis 1675 LSVEC = KSBVEC 1676 DO INUM = 1,NSIM 1677 IVENU = IREDH + IVEC-1 + INUM 1678 IF (IANR12.EQ.1 .OR. IANR12.EQ.3) THEN 1679 CALL CC_GETVEC(LUFC1,FC1AM,LUFC2,FC2AM,LUFS12,FS12AM, 1680 * TRIPLET,IVENU,WRK(LSVEC)) 1681 ELSE IF (IANR12.EQ.2) THEN 1682 CALL CC_GETVEC(LUFC1,FC1AM,LUFS2,FS2AM,LUFS12,FS12AM, 1683 * TRIPLET,IVENU,WRK(LSVEC)) 1684 ELSE 1685 CALL QUIT('This R12 ansatz is not implemented yet') 1686 END IF 1687 LSVEC = LSVEC + NCCVAR 1688 END DO 1689 IF (LOCDBG) THEN 1690 WRITE(LUPRI,'(a,i3,a,i3)') 1691 * 'S * C for ',IVEC,' - ',IVEC+NSIM-1 1692 CALL OUTPUT(WRK(KSBVEC),1,NCCVAR,1,NSIM, 1693 * NCCVAR,NSIM,1,LUPRI) 1694 END IF 1695 END IF 1696C 1697C read in b vectors and extend columns of the Jacobian 1698C (and the metric) in the reduced space 1699C 1700 DO I = 1,NREDH 1701 CALL CC_GETVEC(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 1702 * TRIPLET,I,WRK(KBVEC)) 1703!sonia 1704! IF (CVSEPARA) THEN 1705! write(lupri,*)'SONIAREDCVS 2: Project out i=/ core' 1706! CALL CC_FREEZE_EXCI(WRK(KBVEC),ISYMTR, 1707! & MAXCORE,MAXION, 1708! & NRHFCORE,IRHFCORE, 1709! & NVIRION,IVIRION,LBOTHEXCI) 1710! END IF 1711!sonia 1712 IF (LOCDBG) THEN 1713 WRITE(LUPRI,'(a,i3)') 'C for ',I 1714 CALL OUTPUT(WRK(KBVEC),1,NCCVAR,1,1,NCCVAR,1,1,LUPRI) 1715 END IF 1716 DO JNUM = 1,NSIM 1717 J = IREDH + IVEC -1 + JNUM 1718 REDH( I + (J-1)*MAXRED ) = 1719 * DDOT(NCCVAR,WRK(KBVEC),1,WRK(KSVEC+(JNUM-1)*NCCVAR),1) 1720 IF (LREDS) THEN 1721 REDS( I + (J-1)*MAXRED ) = DDOT(NCCVAR, 1722 * WRK(KBVEC),1,WRK(KSBVEC+(JNUM-1)*NCCVAR),1) 1723 END IF 1724 END DO 1725 END DO 1726 1727 100 CONTINUE 1728C 1729C extend row vectors in reduced spaces matrices if IREDH gt 0 1730C 1731 IF (IREDH.GT.0) THEN 1732 IF (LREDS) THEN 1733 MAXVEC = (LWRK-2*NCCVAR)/NCCVAR 1734 ELSE 1735 MAXVEC = (LWRK-NCCVAR)/NCCVAR 1736 END IF 1737 1738 DO 600 IVEC = 1,NUPVEC,MAXVEC 1739 NSIM = MIN(MAXVEC,NUPVEC+1-IVEC) 1740C 1741 KBVEC = 1 1742 KSVEC = KBVEC + NSIM*NCCVAR 1743 KSBVEC = KSVEC + NCCVAR 1744C 1745 LBVEC = KBVEC 1746 DO INUM = 1,NSIM 1747 IVENU = IREDH + IVEC -1 + INUM 1748 CALL CC_GETVEC(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 1749 * TRIPLET,IVENU,WRK(LBVEC)) 1750 LBVEC = LBVEC + NCCVAR 1751 END DO 1752 IF (LOCDBG) THEN 1753 WRITE(LUPRI,'(a,i3,a,i3)')'C for ',IVEC,'-',IVEC+NSIM-1 1754 CALL OUTPUT(WRK(KBVEC),1,NCCVAR,1,NSIM, 1755 * NCCVAR,NSIM,1,LUPRI) 1756 END IF 1757C 1758C read in s vectors and extend rows of projected matrices 1759!Sonia memo: sigma vectors were not saved on file as projected, so 1760!we need to project again 1761C 1762 DO J = 1,IREDH 1763 CALL CC_GETVEC(LUFR1,FRHO1,LUFR2,FRHO2,LUFR12,FRHO12, 1764 * TRIPLET,J,WRK(KSVEC)) 1765 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1766 if (triplet) then 1767 KCAM1 = KSVEC 1768 KCAMP = KCAM1 + NT1AM(ISYMTR) 1769 KCAMM = KCAMP + NT2AM(ISYMTR) 1770 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1), 1771 & WRK(KCAMP),WRK(KCAMM),ISYMTR, 1772 & MAXCORE,MAXION, 1773 & NRHFCORE,IRHFCORE, 1774 & NVIRION,IVIRION,LBOTHEXCI) 1775 else 1776 !write(lupri,*)'CCRED CVS or IONISATION sigma_',J 1777 CALL CC_FREEZE_EXCI(WRK(KSVEC),ISYMTR, 1778 & MAXCORE,MAXION, 1779 & NRHFCORE,IRHFCORE, 1780 & NVIRION,IVIRION,LBOTHEXCI) 1781 end if 1782 END IF 1783 IF (LRMCORE) THEN 1784 !write(lupri,*)'CCRED remove core orbitals',J 1785 IF (TRIPLET) THEN 1786C Eirik 1787 KCAM1 = KSVEC 1788 KCAMP = KCAM1 + NT1AM(ISYMTR) 1789 KCAMM = KCAMP + NT2AM(ISYMTR) 1790 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), 1791 & WRK(KCAMP),WRK(KCAMM),ISYMTR, 1792 & MAXCORE,MAXION, 1793 & NRHFCORE,IRHFCORE, 1794 & NVIRION,IVIRION,.false.) 1795 ELSE 1796 CALL CC_FREEZE_CORE(WRK(KSVEC),ISYMTR, 1797 & MAXCORE, 1798 & NRHFCORE,IRHFCORE) 1799 END IF 1800 END IF 1801 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1802 write(lupri,*)'SONIA REMOVE CORE FROM L0 in CCRED' 1803 !sonia fixme CCS case 1804 CALL CC_CORE(WRK(KSVEC), 1805 & WRK(KSVEC+NT1AM(ISYMTR)),isymtr) 1806 ENDIF 1807 1808 IF (LREDS) THEN 1809 IF (IANR12.EQ.1 .OR. IANR12.EQ.3) THEN 1810 CALL CC_GETVEC(LUFC1,FC1AM,LUFC2,FC2AM,LUFS12, 1811 * FS12AM,TRIPLET,J,WRK(KSBVEC)) 1812 ELSE IF (IANR12.EQ.2) THEN 1813 CALL CC_GETVEC(LUFC1,FC1AM,LUFS2,FS2AM,LUFS12, 1814 * FS12AM,TRIPLET,J,WRK(KSBVEC)) 1815 ELSE 1816 CALL QUIT('This R12 ansatz is not implemented yet') 1817 END IF 1818 END IF 1819 DO INUM = 1,NSIM 1820 I = IREDH + IVEC -1 + INUM 1821 REDH( I + (J-1)*MAXRED ) = DDOT(NCCVAR, 1822 * WRK(KSVEC),1, WRK(KBVEC+(INUM-1)*NCCVAR),1) 1823 IF (LREDS) THEN 1824 REDS( I + (J-1)*MAXRED ) = DDOT(NCCVAR, 1825 * WRK(KSBVEC),1,WRK(KBVEC+(INUM-1)*NCCVAR),1) 1826 END IF 1827 END DO 1828 END DO 1829 1830 600 CONTINUE 1831 END IF 1832C 1833 IF (LINQCC) THEN 1834 NCAU = 0 1835 IF (NEWCAU) THEN 1836 IF (LIST(1:2).EQ.'RC') THEN 1837 DO IVEC = 1, NVEC 1838 ILSTNR = ISTRVE-1 + IVEC 1839 NCAU = MAX(NCAU,ILRCAU(ILSTNR)) 1840 END DO 1841 IF (IPRLE.GE.1) THEN 1842 WRITE(LUPRI,'(3X,A,I5)') 1843 & 'Nb. of Cauchy iterations in reduced space:', NCAU 1844 END IF 1845 END IF 1846 END IF 1847 END IF 1848C 1849 END IF ! (NUPVEC.GT.0) THEN 1850C 1851 1852 1111 CONTINUE 1853 1854C 1855 IF (NUPVEC.GT.0) THEN 1856 IF (LINQCC) THEN 1857C 1858C set up reduced gradient 1859C 1860C determine maximum number of simultaneous gradients 1861C 1862 MAXVEC = (LWRK-3*NCCVAR)/NCCVAR 1863 DO 1000 IVEC = 1,NVEC,MAXVEC 1864 NSIM = MIN(MAXVEC,NVEC+1-IVEC) 1865C 1866 KGDVEC = 1 1867 KWRK1 = KGDVEC + NSIM*NCCVAR 1868 LWRK1 = LWRK - KWRK1 1869 IF (LWRK1 .LT. 0 ) CALL QUIT('Too little work in CCRED') 1870C 1871 LGDVEC = KGDVEC 1872 DO 1100 INUM = 1,NSIM 1873C 1874C------------------------------- 1875C Get gradient vector. 1876C------------------------------- 1877C 1878C GETGD bruger GD + 3 extra. NCCVAR - skal vaere foerst i memory. 1879 ILSTNR = ISTRVE + IVEC -1 + INUM - 1 1880 CALL CC_GETGD(WRK(LGDVEC),ILSTNR,ISIDE,LIST, 1881 * WRK(KWRK1),LWRK1) 1882C 1883cs ------------------------------- 1884C Project from gradient 1885cs ------------------------------- 1886 IF (LPROJECT) THEN 1887 ISYPR = ILSTSYM(LIST,ILSTNR) 1888 IF (DEBUG) THEN 1889 WRITE(LUPRI,*)'CCRED projection from gradient' 1890 END IF 1891 CALL PRJDRV(ISIDE,ISTATPRJ,ISYPR,1,NCCVAR, 1892 * WRK(LGDVEC),WRK(KWRK1)) 1893 END IF 1894cs ------------------ 1895 IF (LCVSEXCI.or.LIONIZEXCI) THEN 1896 1897 IF (TRIPLET) THEN 1898C Eirik 1899 KCAM1 = LGDVEC 1900 KCAMP = KCAM1 + NT1AM(ISYMTR) 1901 KCAMM = KCAMP + NT2AM(ISYMTR) 1902 CALL CC_FREEZE_TRIPLETEXCI(WRK(KCAM1),WRK(KCAMP), 1903 & WRK(KCAMM),ISYMTR,MAXCORE,MAXION,NRHFCORE, 1904 & IRHFCORE,NVIRION,IVIRION,LBOTHEXCI) 1905 ELSE 1906 CALL CC_FREEZE_EXCI(WRK(LGDVEC),ISYMTR, 1907 & MAXCORE,MAXION, 1908 & NRHFCORE,IRHFCORE, 1909 & NVIRION,IVIRION,LBOTHEXCI) 1910 END IF 1911 END IF 1912 IF (LRMCORE) THEN 1913 write(lupri,*)'CCRED: freeze core in grad' 1914 IF (TRIPLET) THEN 1915C Eirik 1916 KCAM1 = LGDVEC 1917 KCAMP = KCAM1 + NT1AM(ISYMTR) 1918 KCAMM = KCAMP + NT2AM(ISYMTR) 1919 CALL CC_FREEZE_TRIPLETCORE(WRK(KCAM1), WRK(KCAMP), 1920 & WRK(KCAMM),ISYMTR, 1921 & MAXCORE,MAXION, 1922 & NRHFCORE,IRHFCORE, 1923 & NVIRION,IVIRION,.false.) 1924 ELSE 1925 CALL CC_FREEZE_CORE(WRK(LGDVEC),ISYMTR, 1926 & MAXCORE, 1927 & NRHFCORE,IRHFCORE) 1928 END IF 1929 END IF 1930 IF ((LIST.EQ.'L0').and.(LCOR .OR. LSEC )) THEN 1931 write(lupri,*)'SONIA REMOVE CORE in RHS of L0 CCRED' 1932 !Sonia FIXME for CCS 1933 CALL CC_CORE(WRK(LGDVEC), 1934 & WRK(LGDVEC+NT1AM(ISYMTR)),isymtr) 1935 ENDIF 1936C 1937 LGDVEC = LGDVEC + NCCVAR 1938 1100 CONTINUE 1939C 1940 KBVEC = KWRK1 1941 KWRK2 = KBVEC + NCCVAR 1942 LWRK2 = LWRK - KWRK2 1943 IF (LWRK2 .LT. 0 ) CALL QUIT('Too little work in CCRED') 1944C 1945C read in b vectors and set up reduced gradient in SOLEQ 1946C 1947 DO 1200 I = 1,NREDH 1948 CALL CC_GETVEC(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM, 1949 * TRIPLET,I,WRK(KBVEC)) 1950 DO 1300 JNUM = 1,NSIM 1951 J = IVEC -1 + JNUM 1952 SOLEQ(I,J ) = 1953 * -DDOT(NCCVAR,WRK(KBVEC),1, 1954 * WRK(KGDVEC+(JNUM-1)*NCCVAR),1) 1955 1300 CONTINUE 1956 1200 CONTINUE 1957 IF (IPRLE.GT.105) THEN 1958 WRITE(LUPRI,'(/A)')' ** REDUCED GRADIENT i constr **' 1959 CALL OUTPUT(SOLEQ(1,1),1,NREDH,1,NSIM, 1960 * MAXRED,NSIM,1,LUPRI) 1961 END IF 1962 1000 CONTINUE 1963C 1964 END IF 1965 1966 IF (IPRLE.GE.15 .OR. LOCDBG) THEN 1967 IF (LREDS) THEN 1968 WRITE(LUPRI,'(/A)')' ** REDUCED OVERLAP MATRIX **' 1969 CALL OUTPUT(REDS,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 1970 END IF 1971 WRITE(LUPRI,'(/A)') ' ** REDUCED HESSIAN MATRIX **' 1972 CALL OUTPUT(REDH,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 1973 END IF 1974C 1975 END IF ! (NUPVEC.GT.0) THEN 1976C 1977C ************************************************************** 1978C Section 2 1979C find solution vectors in reduced space 1980C ************************************************************** 1981C 1982 IF ( LINQCC ) THEN 1983C 1984C 1985C ************************************************************** 1986C Section 2: FIND SOLUTIONS TO THE NVEC SET OF NEWTON-RAPHSON EQ. 1987C ************************************************************** 1988C 1989C SOLVE NVEC SET OF N-R EQ. 1990C Copy REDH to SLEVEC 1991C 1992C SOLEQ contain reduced gradients 1993C 1994 KMA = MAXRED*MAXRED 1995 DO 2100 IVEC = 1,NVEC 1996 CALL DCOPY(KMA,REDH, 1,WRK,1) 1997 IF (CCR12 .AND. DABS(EIVAL(IVEC)).GT.1.0D-15) THEN 1998 CALL DAXPY(NREDH*NREDH,-EIVAL(IVEC),REDS,1,WRK,1) 1999 ELSE 2000 DO 2200 I = 1,NREDH 2001 IADR = (I-1)*MAXRED + I 2002 WRK(IADR) = WRK(IADR) - EIVAL(IVEC) 2003 2200 CONTINUE 2004 END IF 2005 IF (IPRLE.GE.15 .OR. LOCDBG) THEN 2006 IF (LREDS) THEN 2007 WRITE(LUPRI,'(/A)')' ** REDUCED OVERLAP MATRIX **' 2008 CALL OUTPUT(REDS,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2009 END IF 2010 WRITE(LUPRI,'(/A)')' ** REDUCED HESSIAN - w1 MATRIX **' 2011 WRITE(LUPRI,'(A,F10.6)') ' FREQUENCY: ',EIVAL(IVEC) 2012 CALL OUTPUT(WRK,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2013 WRITE(LUPRI,'(/A)')' ** REDUCED GRADIENT **' 2014 CALL OUTPUT(SOLEQ(1,IVEC),1,NREDH,1,1,MAXRED,1,1,LUPRI) 2015 END IF 2016 NGD = 1 2017 CALL DGESOL(NGD,NREDH,MAXRED,MAXRED,WRK, 2018 * SOLEQ(1,IVEC),WRK(1+KMA),INFO) 2019C CALL DGESOL(NSIM,N,LDA,LDB,A,B,KPVT,INFO) 2020 IF (INFO.NE.0) WRITE(LUPRI,7000)IVEC,INFO 2021 7000 FORMAT(/' ***CCRED*** SOLUTION NOT OBTAINED TO LINEAR EQUATIONS' 2022 * ,I5,/' CHECK IF HESSIAN MATRIX IS SINGULAR,', 2023 * ' INFO from DGESOL/DSPSOL =',I7) 2024 IF (IPRLE.GE.15 .OR. LOCDBG) THEN 2025 WRITE(LUPRI,'(/A)')' ** REDUCED SOLUTION VECTOR **' 2026 CALL OUTPUT(SOLEQ(1,IVEC),1,NREDH,1,1,MAXRED,1,1,LUPRI) 2027 END IF 2028 2100 CONTINUE 2029 IF (IPRLE.GT.15 .OR. LOCDBG) THEN 2030 WRITE(LUPRI,*)' * SOLUTION TO',NVEC,' NEWTON EQUATIONS *' 2031 CALL OUTPUT(SOLEQ(1,1),1,NREDH,1,NVEC,MAXRED,NVEC, 2032 * 1,LUPRI) 2033 END IF 2034C 2035 END IF 2036C 2037 IF (LINQCC.AND.NEWCAU.AND.NCAU.GT.1) THEN 2038 NCAU = NCAU - 1 2039 ISYM = ILSTSYM(LIST,ISTRVE) 2040 NVARPT = NCCVAR + 2*NALLAI(ISYM) 2041 MAXVEC = (LWRK-NVARPT)/NCCVAR 2042 DO IVEC = 1, NVEC, MAXVEC 2043 NBVEC = MIN(MAXVEC,NVEC+1-IVEC) 2044 KWRK1 = 1 + NCCVAR*NBVEC 2045 LWRK1 = LWRK - KWRK1 2046 IF (LWRK1 .LT. 0) THEN 2047 CALL QUIT('Out of memory in CCRED') 2048 END IF 2049 CALL CCCONV(LUFC1,FC1AM,LUFC2,FC2AM,LUFC12,FC12AM,TRIPLET, 2050 * CCR12,NREDH,IVEC,NBVEC,SOLEQ,WRK(1),WRK(KWRK1)) 2051 DO I = 1, NBVEC 2052 IDX = IVEC-1+I 2053 CALL CC_SAVE(WRK(1+(I-1)*NCCVAR),ISTRVE-1+IDX, 2054 * LIST,WRK(KWRK1),LWRK1) 2055 END DO 2056 END DO 2057 GOTO 1111 2058 END IF 2059C 2060 IF (.NOT. LINQCC ) THEN 2061C 2062C 2063C ************************************************************** 2064C Section 2: Solve reduced generalized eigenvalue problem in subspace 2065C ************************************************************** 2066C 2067 MATZ = 1 2068 KWI = 1 2069 KAMAT= KWI + MAXRED 2070 KIV1 = KAMAT+ MAXRED*MAXRED 2071 KFV1 = KIV1 + MAXRED 2072 KEND = KFV1 + MAXRED 2073 IF (LREDS) THEN 2074 KSMAT = KEND 2075 KEND = KSMAT + MAXRED*MAXRED 2076 KDENOM = KFV1 2077 END IF 2078 2079 LEND = LWRK - KEND 2080 IF (KEND .GT. LWRK) CALL ERRWRK('LERED.RG',KEND,LWRK) 2081 2082 IF (LREDS) THEN 2083 ! solve generalized eigenvalue problem for a real general 2084 ! jacobian and a nontrivial metric 2085 CALL DCOPY(MAXRED*MAXRED,REDH,1,WRK(KAMAT),1) 2086 CALL DCOPY(MAXRED*MAXRED,REDS,1,WRK(KSMAT),1) 2087c 2088 CALL RGG(MAXRED,NREDH,WRK(KAMAT),WRK(KSMAT),EIVAL, 2089 * WRK(KWI),WRK(KDENOM),MATZ,SOLEQ,IERR) 2090 DO I = 1, NREDH 2091 IF (ABS(WRK(KDENOM-1+I)).GT.THRZER) THEN 2092 EIVAL(I) = EIVAL(I)/WRK(KDENOM-1+I) 2093 WRK(KWI-1+i) = WRK(KWI-1+I)/WRK(KDENOM-1+I) 2094 ELSE 2095 EIVAL(I) = D1/THRZER 2096 WRK(KWI-1+i) = WRK(KWI-1+I)/THRZER 2097 END IF 2098 END DO 2099 ELSE 2100 ! solve eigenvalue problem for a real general jacobian 2101 CALL DCOPY(MAXRED*MAXRED,REDH,1,WRK(KAMAT),1) 2102 CALL RG(MAXRED,NREDH,WRK(KAMAT),EIVAL,WRK(KWI),MATZ,SOLEQ, 2103 * WRK(KIV1),WRK(KFV1),IERR) 2104 END IF 2105 2106 IF (IPRLE .GE. 70 .OR. IERR .NE. 0) THEN 2107 WRITE (LUPRI,'(/A)') ' REDUCED EIGENVALUES real part:' 2108 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2109 CALL OUTPUT(EIVAL,1,NREDH,1,1,NREDH,MAXRED,1,LUPRI) 2110 WRITE (LUPRI,'(/A)') 2111 * ' REDUCED EIGENVALUES imaginary part:' 2112 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2113 CALL OUTPUT(WRK(KWI),1,NREDH,1,1,NREDH,MAXRED,1,LUPRI) 2114 WRITE (LUPRI,'(/A)') ' REDUCED EIGENVECTORS :' 2115 WRITE (LUPRI,'(A)') ' before sort of eigenvalues' 2116 CALL OUTPUT(SOLEQ,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2117 END IF 2118 IF ( IERR.NE.0 ) THEN 2119 WRITE(LUPRI,'(/A,I5)') 2120 * ' EIGENVALUE PROBLEM NOT CONVERGED IN RG, IERR =',IERR 2121 CALL QUIT(' CCRED: EIGENVALUE EQUATION NOT CONVERGED ') 2122 END IF 2123C 2124 CALL RGORD(MAXRED,NREDH,EIVAL,WRK(KWI),SOLEQ,.FALSE.) 2125C 2126 ICPLX = 0 2127 DO I=1,NVEC 2128 IF (WRK(I) .NE. D0) THEN 2129 ICPLX = ICPLX + 1 2130 WRITE(LUPRI,'(I10,1P,2D15.8,A/)') I,EIVAL(I),WRK(I), 2131 * ' *** CCRED WARNING **** COMPLEX VALUE.' 2132 END IF 2133 END DO 2134C 2135 IF (IPRLE .GE.1 .OR. ICPLX .GT. 0) THEN 2136 WRITE (LUPRI,'(/A)') ' REDUCED EIGENVALUES real part:' 2137 CALL OUTPUT(EIVAL,1,NVEC,1,1,MAXRED,1,1,LUPRI) 2138 ENDIF 2139C 2140 IF (IPRLE .GE.11 .OR. ICPLX .GT. 0) THEN 2141 WRITE (LUPRI,'(/A)') ' REDUCED EIGENVALUES imaginary part:' 2142 CALL OUTPUT(WRK(KWI),1,NVEC,1,1,MAXRED,1,1,LUPRI) 2143 END IF 2144C 2145 IF (IPRLE .GE. 15) THEN 2146 WRITE (LUPRI,'(/A)') ' REDUCED EIGENVECTORS :' 2147 CALL OUTPUT(SOLEQ,1,NREDH,1,NVEC,MAXRED,MAXRED,1,LUPRI) 2148 END IF 2149C 2150 IF (ICPLX .GT. 0) 2151 * WRITE(LUPRI,*) '**WARNING CCRED: COMPLEX EIGENVALUES.' 2152C 2153 END IF 2154C 2155 IF (IPRLE.GE.15 ) THEN 2156 WRITE(LUPRI,'(/A)')' REDUCED HESSIAN MATRIX:' 2157 CALL OUTPUT(REDH,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2158 IF (LREDS) THEN 2159 WRITE(LUPRI,'(/A)')' REDUCED METRIC MATRIX:' 2160 CALL OUTPUT(REDS,1,NREDH,1,NREDH,MAXRED,MAXRED,1,LUPRI) 2161 END IF 2162 END IF 2163C 2164C *** End of subroutine CCRED 2165C 2166 RETURN 2167 END 2168C /* Deck rgord */ 2169 SUBROUTINE RGORD(NM,N,WR,WI,Z,LSORT) 2170C 2171C 15-Aug-1989 Hans Joergen Aa. Jensen 2172C modified by Ove Christiansen 20-dec 1994 2173C to normalize correct! 2174C 2175C Reorder and normalize eigenpairs from RG. 2176C 2177C LSORT = .FALSE. --> ascending order 2178C LSORT = .TRUE. --> descending order 2179#include "implicit.h" 2180 LOGICAL LSORT 2181 DIMENSION WR(N), WI(N), Z(NM,N) 2182 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 ) 2183 DO 100 I = 1,N-1 2184 ILOW = I 2185 XLOW = WR(I) 2186 DO 50 J = I+1,N 2187 IF (.NOT.LSORT) THEN 2188 IF (WR(J) .LT. XLOW) THEN 2189 ILOW = J 2190 XLOW = WR(J) 2191 END IF 2192 ELSE 2193 IF (WR(J) .GT. XLOW) THEN 2194 ILOW = J 2195 XLOW = WR(J) 2196 END IF 2197 END IF 2198 50 CONTINUE 2199 IF (ILOW .NE. I) THEN 2200 WR(ILOW) = WR(I) 2201 WR(I) = XLOW 2202 XLOW = WI(ILOW) 2203 WI(ILOW) = WI(I) 2204 WI(I) = XLOW 2205 CALL DSWAP(N,Z(1,I),1,Z(1,ILOW),1) 2206 END IF 2207 100 CONTINUE 2208C 2209C Normalize eigenvectors 2210C 2211 DO 200 I = 1,N 2212 ZNRM = DDOT(N,Z(1,I),1,Z(1,I),1) 2213 ZNRM = D1 / SQRT(ZNRM) 2214 IMAX = IDAMAX(N,Z(1,I),1) 2215 IF (Z(IMAX,I) .LT. D0) ZNRM = -ZNRM 2216 CALL DSCAL(N,ZNRM,Z(1,I),1) 2217 200 CONTINUE 2218 RETURN 2219 END 2220C /* Deck ccort */ 2221 SUBROUTINE CCORT(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 2222 * LUFS2,S2FIL,LUFS12,S12FIL, 2223 * TRIPLET,NBADD,NBPREV,NCCVAR, 2224 * BVECS,ISYMTR,THRLDP,WRK,LWRK,IPRLE) 2225C 2226C Modified version of ORTBVLE of Hans Jorgen Aa Jensen 2227C specialized for CC 2228C 2229C Purpose: 2230C Orthogonalize new b-vectors against all previous b-vectors 2231C and among themselves, and renormalize. 2232C (Orthogonalization is performed twice if round-off is large, 2233C if larger than THRRND). 2234C 2235C Input: 2236C LUFC1 - file with singles parts of trial vectors with name C1FIL 2237C LUFC2 - file with doubles parts of trial vectors with name C2FIL 2238C LUFC12 - file with R12 doubles parts of trial vec. with name C12FIL 2239C BVECS - non-orthogonal b-vectors 2240C NBADD - number of BVECS 2241C NBPREV - number of trial vectors on LUB 2242C THRLDP - threshold for linear dependence 2243C 2244C 2245C Output: 2246C NBADD trial vectors added to previous 2247C NBPREV updated to contain NBADD new trial vectors 2248C 2249C Scratch: 2250C require WRK of length NCCVAR 2251C 2252#include "implicit.h" 2253#include "priunit.h" 2254#include "ccsdinp.h" 2255#include "ccsdsym.h" 2256#include "r12int.h" 2257C 2258 CHARACTER C1FIL*(*),C2FIL*(*),C12FIL*(*),S2FIL*(*),S12FIL*(*) 2259 DIMENSION BVECS(NCCVAR,*),WRK(NCCVAR) 2260 LOGICAL TRIPLET, LOCDBG 2261 PARAMETER (THRRND=1.D-4, THRTT=1.D-6, THRZER=1.0D-35) 2262 PARAMETER (THRTST=1.0D-12) 2263 PARAMETER (D1 = 1.0D0) 2264 PARAMETER (LOCDBG = .FALSE.) 2265 INTEGER SVEC, RVEC, LVEC 2266C 2267C 2268C 2269 NSAVE = NBADD 2270C 2271 IF ((IPRLE.GT.120) .OR. LOCDBG) THEN 2272 WRITE (LUPRI,'(//A/)') ' --- Debug output from CCORT ---' 2273 END IF 2274C 2275C work space allocation 2276C 2277 KEND1 = 1 2278 IF (CCR12) THEN 2279 KSVECS = KEND1 2280 KEND1 = KSVECS+ NCCVAR*NBADD 2281 END IF 2282 LWRK1 = LWRK - KEND1 2283 IF (LWRK1.LT.0) THEN 2284 WRITE(LUPRI,*) 'LWRK, KEND1: ',LWRK, KEND1 2285 CALL QUIT('Insufficient memory in CCORT') 2286 END IF 2287C 2288C Normalize input BVECS and remove 2289C vectors with norm .le. THRZER 2290C 2291 NNUM = 0 2292 DO I = 1, NBADD 2293 2294 IF (CCR12) THEN 2295 CALL DCOPY(NCCVAR,BVECS(1,I),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2296 CALL CC_RVEC(LUFS12,S12FIL,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2297 * NBPREV+I,WRK(KSVECS+NCCVAR*(I-1)+NT1AM(ISYMTR)+ 2298 * NT2AM(ISYMTR))) 2299 IF (IANR12.EQ.2) THEN 2300 CALL CC_RVEC(LUFS2,S2FIL,NT2AM(ISYMTR),NT2AM(ISYMTR), 2301 * NBPREV+I,WRK(KSVECS+NCCVAR*(I-1)+ 2302 * NT1AM(ISYMTR))) 2303 END IF 2304 END IF 2305 2306 IF (CCR12) THEN 2307C Check norm of r12 part: 2308 TST = DDOT(NTR12AM(ISYMTR),BVECS((1+NT1AM(ISYMTR)+ 2309 * NT2AM(ISYMTR)),I),1,WRK(KSVECS+NCCVAR*(I-1)+ 2310 * NT1AM(ISYMTR)+NT2AM(ISYMTR)),1) 2311 IF (TST .LT. -THRTST) THEN 2312 WRITE(LUPRI,*) '(CCORT) Norm^2 (R12-part) = ',TST 2313 WRITE(LUPRI,'(A,I4,A)') 'WARNING: R12 part of input '// 2314 * 'vector no. ',I,' in CCORT has norm < 0' 2315C WRITE(LUPRI,*) '------- Check your '// 2316C * 'auxiliary basis and/or approximations!' 2317C CALL QUIT('R12 input vector in CCORT has norm < 0') 2318 END IF 2319 TT = DDOT(NCCVAR,BVECS(1,I),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2320 ELSE 2321 TT = DDOT(NCCVAR,BVECS(1,I),1,BVECS(1,I),1) 2322 END IF 2323 IF ((IPRLE.GT.20) .OR. LOCDBG) THEN 2324 WRITE (LUPRI,*) ' I, initial norm of BVECS(1,I) :',I,TT 2325 END IF 2326 IF (TT.LT.THRZER) THEN 2327 IF (IPRLE.GT.2) WRITE (LUPRI,3240) I,TT 2328 3240 FORMAT(/' (CCORT) new b-vector no.',I3, 2329 * /' is removed at input because of too small norm;' 2330 * /' norm of input vector is ',1P,D12.5) 2331 ELSE 2332 NNUM = NNUM + 1 2333 TT = SQRT(TT) 2334 IF (TT.LT.THRTT) THEN 2335 CALL DSCAL(NCCVAR,(D1/THRTT),BVECS(1,I),1) 2336 IF (CCR12) CALL 2337 & DSCAL(NCCVAR,(D1/THRTT),WRK(KSVECS+NCCVAR*(I-1)),1) 2338 IF (CCR12) THEN 2339 TT=DDOT(NCCVAR,BVECS(1,I),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2340 ELSE 2341 TT = DDOT(NCCVAR,BVECS(1,I),1,BVECS(1,I),1) 2342 END IF 2343 TT = SQRT(TT) 2344 END IF 2345 CALL DSCAL (NCCVAR,(D1/TT),BVECS(1,I),1) 2346 IF (CCR12) CALL 2347 & DSCAL(NCCVAR,(D1/TT),WRK(KSVECS+NCCVAR*(I-1)),1) 2348 IF ( I.GT.NNUM) THEN 2349 CALL DCOPY(NCCVAR,BVECS(1,I),1,BVECS(1,NNUM),1) 2350 IF (CCR12) CALL DCOPY(NCCVAR,WRK(KSVECS+NCCVAR*(I-1)),1, 2351 * WRK(KSVECS+NCCVAR*(NNUM-1)),1) 2352 END IF 2353 END IF 2354 END DO 2355C 2356 IF (NBADD.NE.NNUM) THEN 2357 WRITE(LUPRI,*)NNUM-NBADD,' vectors with norm le:', 2358 * THRZER,' REMOVED' 2359 NBADD = NNUM 2360 END IF 2361 IF (NBADD .EQ. 0) THEN 2362 WRITE (LUPRI,'(/1X,2A,I5)') 2363 * 'CCORT,**WARNING** number vectors added is zero', 2364 * ' after check for vector norm , NBADD :',NBADD 2365 END IF 2366C 2367C work space allocation 2368C 2369 KBPRE = KEND1 2370 KEND1 = KBPRE + NCCVAR 2371 IF (CCR12) THEN 2372 KSPRE = KEND1 2373 KEND1 = KSPRE + NCCVAR 2374 END IF 2375 IROUND = 0 2376 ITURN = 0 2377 2378 1500 CONTINUE 2379 ITURN = ITURN + 1 2380C 2381C Orthogonalize new b-vectors against previous b-vectors 2382C 2383 DO K = 1, NBPREV 2384 IF (NBADD .GT. 0) THEN 2385 CALL CC_GETVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 2386 * TRIPLET,K,WRK(KBPRE)) 2387 IF (CCR12) THEN 2388C put S * b into WRK(KSPRE) 2389 CALL DCOPY(NCCVAR,WRK(KBPRE),1,WRK(KSPRE),1) 2390 IF (IANR12.EQ.1 .OR. IANR12.EQ.3) THEN 2391 CALL CC_GETVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFS12,S12FIL, 2392 * TRIPLET,K,WRK(KSPRE)) 2393 ELSE IF (IANR12.EQ.2) THEN 2394 CALL CC_GETVEC(LUFC1,C1FIL,LUFS2,S2FIL,LUFS12,S12FIL, 2395 * TRIPLET,K,WRK(KSPRE)) 2396 ELSE 2397 CALL QUIT('Unknown IANR12 in CCORT') 2398 END IF 2399 END IF 2400 DO J = 1, NBADD 2401 IF (CCR12) THEN 2402 TT=-DDOT(NCCVAR,WRK(KBPRE),1,WRK(KSVECS+NCCVAR*(J-1)),1) 2403 ELSE 2404 TT = -DDOT(NCCVAR,WRK(KBPRE),1,BVECS(1,J),1) 2405 END IF 2406 IF ((IPRLE.GT.20) .OR. LOCDBG) THEN 2407 WRITE (LUPRI,*) 2408 * ' K (prev.vec.), J (new type a), overlap:', 2409 * K,J,TT 2410 END IF 2411 CALL DAXPY(NCCVAR,TT,WRK(KBPRE),1,BVECS(1,J),1) 2412 IF (CCR12) THEN 2413 CALL DAXPY(NCCVAR,TT,WRK(KSPRE),1,WRK(KSVECS+ 2414 * NCCVAR*(J-1)),1) 2415 END IF 2416 END DO 2417 ELSE 2418 WRITE (LUPRI,'(/1X,A,2I5)') 2419 * 'CCORT, number vectors added is zero , NBADD :',NBADD 2420 END IF 2421 END DO 2422CCN 2423C IF (LOCDBG) THEN 2424C IF (CCR12) THEN 2425C KOLP = KEND1 2426C KEND2 = KOLP + NBADD*NBADD 2427C IF (KEND2.LT.0) 2428C * CALL QUIT('Insufficient memory in CCORT') 2429C DO I = 1, NBADD 2430C DO J = 1, NBADD 2431C WRK(KOLP-1+NBADD*(J-1)+I) = DDOT(NTR12AM(ISYMTR), 2432C * BVECS((1+NT1AM(ISYMTR)+NT2AM(ISYMTR)),I),1, 2433C * WRK(KSVECS+NCCVAR*(J-1)+ 2434C * NT1AM(ISYMTR)+NT2AM(ISYMTR)),1) 2435C END DO 2436C END DO 2437C WRITE(LUPRI,*)' R12 PART OF OVERLAP:' 2438C CALL OUTPUT(WRK(KOLP),1,NBADD,1,NBADD,NBADD,NBADD,1,LUPRI) 2439C WRITE(LUPRI,*)' R12 PART OF B vectors:' 2440C CALL OUTPUT(BVECS,1+NT1AM(ISYMTR)+NT2AM(ISYMTR),NCCVAR, 2441C * 1,NBADD,NCCVAR,NBADD,1,LUPRI) 2442C WRITE(LUPRI,*)' R12 PART OF S x B vectors:' 2443C CALL OUTPUT(WRK(KSVECS),1+NT1AM(ISYMTR)+NT2AM(ISYMTR),NCCVAR, 2444C * 1,NBADD,NCCVAR,NBADD,1,LUPRI) 2445C END IF 2446C END IF 2447CCN 2448C 2449C Normalize and orthogonalize new vectors against each other 2450C and remove vectors that are linear dependent 2451C 2452 NNUM = 0 2453 DO I = 1, NBADD 2454 DO J = 1, NNUM 2455 IF (CCR12) THEN 2456 T1 = DDOT(NCCVAR,BVECS(1,J),1,WRK(KSVECS+NCCVAR*(J-1)),1) 2457 TT = -DDOT(NCCVAR,BVECS(1,J),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2458 ELSE 2459 T1 = DDOT(NCCVAR,BVECS(1,J),1,BVECS(1,J),1) 2460 TT = -DDOT(NCCVAR,BVECS(1,J),1,BVECS(1,I),1) 2461 END IF 2462 IF ((IPRLE.GT.20) .OR. LOCDBG) THEN 2463 WRITE (LUPRI,*) 2464 * ' I (new type a), J (new type a), norm bJ, overlap:', 2465 * I,J,T1,TT 2466 END IF 2467 TT = TT / T1 2468 CALL DAXPY(NCCVAR,TT,BVECS(1,J),1,BVECS(1,I),1) 2469 IF (CCR12) THEN 2470 CALL DAXPY(NCCVAR,TT,WRK(KSVECS+NCCVAR*(J-1)),1, 2471 * WRK(KSVECS+NCCVAR*(I-1)),1) 2472 END IF 2473 END DO 2474 IF (CCR12) THEN 2475C Check norm of r12 part: 2476 TST = DDOT(NTR12AM(ISYMTR),BVECS((1+NT1AM(ISYMTR)+ 2477 * NT2AM(ISYMTR)),I),1,WRK(KSVECS+NCCVAR*(I-1)+ 2478 * NT1AM(ISYMTR)+NT2AM(ISYMTR)),1) 2479 IF (TST .LT. -THRTST) THEN 2480 WRITE(LUPRI,*) '(CCORT) Norm^2 (R12-part) = ',TST 2481 WRITE(LUPRI,'(A,I4,A)') 'WARNING: R12 part of result '// 2482 * 'vector no. ',I,' in CCORT has norm < 0' 2483C WRITE(LUPRI,*) '------- Check your '// 2484C * 'auxiliary basis and/or approximations!' 2485C CALL QUIT('R12 result vector in CCORT has norm < 0') 2486 END IF 2487 TT = DDOT(NCCVAR,BVECS(1,I),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2488 ELSE 2489 TT = DDOT(NCCVAR,BVECS(1,I),1,BVECS(1,I),1) 2490 END IF 2491 IF ((IPRLE.GT.20) .OR. LOCDBG) THEN 2492 WRITE (LUPRI,*) 2493 * ' I (new type a), norm before renormalization:', 2494 * I,SQRT(TT) 2495 END IF 2496 IF (TT .LE. THRLDP) THEN 2497 IF (IPRLE.GT.2) WRITE (LUPRI,3250) I,TT 2498 3250 FORMAT(/' (CCORT) new b-vector no.',I3, 2499 * /' is removed because of linear dependence;' 2500 * /' norm of vector after Gram-Schmidt''s orthogonalization ' 2501 * ,1P,D12.5) 2502 ELSE 2503 NNUM = NNUM + 1 2504 IF (TT .LT. THRRND) IROUND = IROUND+1 2505 TT = SQRT(TT) 2506 IF (TT.LT.THRTT) THEN 2507 CALL DSCAL(NCCVAR,(D1/THRTT),BVECS(1,I),1) 2508 IF (CCR12) CALL 2509 & DSCAL(NCCVAR,(D1/THRTT),WRK(KSVECS+NCCVAR*(I-1)),1) 2510 IF (CCR12) THEN 2511 TT=DDOT(NCCVAR,BVECS(1,I),1,WRK(KSVECS+NCCVAR*(I-1)),1) 2512 ELSE 2513 TT = DDOT(NCCVAR,BVECS(1,I),1,BVECS(1,I),1) 2514 END IF 2515 TT = SQRT(TT) 2516 END IF 2517 CALL DSCAL(NCCVAR,(D1/TT),BVECS(1,I),1) 2518 IF (CCR12) CALL 2519 & DSCAL(NCCVAR,(D1/TT),WRK(KSVECS+NCCVAR*(I-1)),1) 2520 IF ( I.NE.NNUM) THEN 2521 CALL DCOPY(NCCVAR,BVECS(1,I),1,BVECS(1,NNUM),1) 2522 IF (CCR12) CALL DCOPY(NCCVAR,WRK(KSVECS+NCCVAR*(I-1)),1, 2523 * WRK(KSVECS+NCCVAR*(NNUM-1)),1) 2524 END IF 2525 END IF 2526 END DO 2527C 2528 IF (NBADD.NE.NNUM) THEN 2529 IF (IPRLE.GT.2) THEN 2530 WRITE(LUPRI,*)' After Gram Schmidt orthonormalization' 2531 WRITE(LUPRI,*) NNUM-NBADD,' vectors with norm le:' 2532 * ,THRZER,' REMOVED' 2533 END IF 2534 NBADD = NNUM 2535 END IF 2536 IF (IROUND.GT.0 .AND. ITURN.EQ.1 ) GO TO 1500 2537C 2538C Save new vector together with old ones on LUS 2539C 2540 NLAST = NBPREV + NBADD 2541 DO I = 1, NBADD 2542 IVENU = NBPREV + I 2543 CALL CC_PUTVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 2544 * TRIPLET,IVENU,BVECS(1,I)) 2545 IF (CCR12) THEN 2546 CALL CC_WVEC(LUFS12,S12FIL,NTR12AM(ISYMTR),NTR12AM(ISYMTR), 2547 * IVENU,WRK(KSVECS+NT1AM(ISYMTR)+NT2AM(ISYMTR)+ 2548 * NCCVAR*(I-1))) 2549 IF (IANR12.EQ.2) THEN 2550 CALL CC_WVEC(LUFS2,S2FIL,NT2AM(ISYMTR),NT2AM(ISYMTR), 2551 * IVENU,WRK(KSVECS+NT1AM(ISYMTR)+NCCVAR*(I-1))) 2552 END IF 2553 END IF 2554 END DO 2555C 2556C update NBPREV 2557C 2558 NBPREV = NBPREV + NBADD 2559C 2560C 2561 IF ( NSAVE .NE. NBADD ) THEN 2562 IF (IPRLE.GT.2) WRITE(LUPRI,4310)NSAVE,NBADD 2563 4310 FORMAT(/' (CCORT) trial vectors reduced from',I3,' to',I3) 2564 END IF 2565C 2566C test if trial vectors on disk are orthonormal 2567C 2568 IF (( DEBUG .AND. IPRLE .GT. 20) .OR. LOCDBG) THEN 2569 KVEC = 1 2570 KOLP = KVEC + NBPREV*NCCVAR 2571 KWRK1 = KOLP + NBPREV*NBPREV 2572 IF (CCR12) THEN 2573 SVEC = KWRK1 2574 KWRK1 = SVEC + NBPREV*NCCVAR 2575 END IF 2576 LWRK1 = LWRK - KWRK1 2577 IF ( LWRK1 .LE.0) CALL QUIT('CCORT TOO LITTLE WRK') 2578 LVEC = KVEC 2579 RVEC = SVEC 2580 DO I = 1, NBPREV 2581 CALL CC_GETVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFC12,C12FIL, 2582 * TRIPLET,I,WRK(LVEC)) 2583 IF (CCR12) THEN 2584 IF (IANR12.EQ.1 .OR. IANR12.EQ.3) THEN 2585 CALL CC_GETVEC(LUFC1,C1FIL,LUFC2,C2FIL,LUFS12,S12FIL, 2586 * TRIPLET,I,WRK(RVEC)) 2587 ELSE IF (IANR12.EQ.2) THEN 2588 CALL CC_GETVEC(LUFC1,C1FIL,LUFS2,S2FIL,LUFS12,S12FIL, 2589 * TRIPLET,I,WRK(RVEC)) 2590 ELSE 2591 CALL QUIT('Unknown IANR12 in CCORT') 2592 END IF 2593 END IF 2594 LVEC = LVEC + NCCVAR 2595 RVEC = RVEC + NCCVAR 2596 END DO 2597 WRITE(LUPRI,*)'*** CCORT TEST ****' 2598 WRITE(LUPRI,*)' OVERLAP BETWEEN',NBPREV,' TRIAL VECTORS' 2599 IF (CCR12) THEN 2600 CALL PROVLP(WRK(KVEC),WRK(SVEC),NCCVAR,NBPREV,WRK(KOLP),LUPRI) 2601 ELSE 2602 CALL PROVLP(WRK(KVEC),WRK(KVEC),NCCVAR,NBPREV,WRK(KOLP),LUPRI) 2603 END IF 2604 IF (IPRLE.GT.50) THEN 2605 CALL OUTPUT(WRK(KVEC),1,NCCVAR,1,NBPREV,NCCVAR,NBPREV, 2606 * 1,LUPRI) 2607 END IF 2608 CALL FLSHFO(LUPRI) 2609 END IF 2610C 2611C *** End of subroutine CCORT 2612C 2613 IF (IPRLE .GT. 20) THEN 2614 WRITE (LUPRI,'(/A//)')' --- End of debug output from CCORT ---' 2615 END IF 2616 2617 RETURN 2618 END 2619C /* Deck provlp */ 2620 SUBROUTINE PROVLP(VECS,SVECS,LVEC,NVECS,OVLP,LUPRI) 2621C 2622C 22-Aug-1989 Hans Joergen Aa. Jensen 2623C 2624C Print overlap matrix between the NVECS input vectors. 2625C 2626#include "implicit.h" 2627 DIMENSION VECS(LVEC,NVECS), OVLP(*), SVECS(LVEC,NVECS) 2628C 2629 IJ = 0 2630 DO 200 I = 1,NVECS 2631 DO 100 J = 1,NVECS 2632 IJ = IJ + 1 2633 OVLP(IJ) = DDOT(LVEC,VECS(1,I),1,SVECS(1,J),1) 2634 100 CONTINUE 2635 200 CONTINUE 2636 WRITE (LUPRI,'(/A)') ' Overlap matrix :' 2637 CALL OUTPUT(OVLP,1,NVECS,1,NVECS,NVECS,NVECS,1,LUPRI) 2638 RETURN 2639 END 2640C /* Deck lediag */ 2641 SUBROUTINE LEDIAG(NDIAG,DIAG,DTEST) 2642C 2643C 5-Nov-1988 Hans Joergen Aa. Jensen 2644C 2645C Calculate inverse diagonal, 2646C for generalized conjugate gradient algorithm. 2647C 2648#include "implicit.h" 2649#include "priunit.h" 2650 DIMENSION DIAG(NDIAG) 2651 PARAMETER ( D1 = 1.0D0 ) 2652C 2653 DO 100 I = 1,NDIAG 2654 IF (ABS(DIAG(I)).LE.DTEST) THEN 2655 DIAG(I) = D1 / SIGN(DTEST,DIAG(I)) 2656 ELSE 2657 DIAG(I) = D1 / DIAG(I) 2658 END IF 2659 100 CONTINUE 2660 RETURN 2661 END 2662C /* Deck cc_save */ 2663 SUBROUTINE CC_SAVE(VEC,ILSTNR,LIST,WORK,LWORK) 2664C 2665C Ove Christiansen 14-11-1996 2666C Write vec to file based on list and list number. 2667C 2668#include "implicit.h" 2669#include "priunit.h" 2670#include "cclr.h" 2671#include "ccorb.h" 2672#include "leinf.h" 2673#include "ccsdsym.h" 2674#include "ccsdinp.h" 2675Cholesky 2676#include "maxorb.h" 2677#include "ccdeco.h" 2678Cholesky 2679C 2680 CHARACTER*(*) LIST 2681 CHARACTER*14 MODEL, LABR12 2682C 2683 DIMENSION VEC(*), WORK(LWORK) 2684C 2685 INTEGER ILSTSYM 2686C 2687 MODEL = ' ' 2688 IF (CCR12 .AND. (.NOT.(CCS.OR.CIS))) THEN 2689 ILEN = 4 2690 LABR12 = '-R12' 2691 ELSE 2692 ILEN = 1 2693 LABR12 = ' ' 2694 END IF 2695 2696 MODEL = 'CCSD'//LABR12(1:ILEN) 2697 IF (CCS) MODEL = 'CCS' 2698 IF (CC2) MODEL = 'CC2'//LABR12(1:ILEN) 2699 IF (CC3) MODEL = 'CC3'//LABR12(1:ILEN) 2700 IF (CC1A) MODEL = 'CCSDT-1A'//LABR12(1:ILEN) 2701 IF (CC1B) MODEL = 'CCSDT-1B'//LABR12(1:ILEN) 2702C 2703 ISYM = ILSTSYM(LIST,ILSTNR) 2704C 2705 KT1 = 1 2706 KT2 = 1 + NT1AM(ISYM) 2707C 2708 IOPT = 3 2709Cholesky 2710 IF (CHOINT .AND. CC2) IOPT = 1 2711Cholesky 2712 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY,VEC(KT1), 2713 * VEC(KT2),WORK,LWORK) 2714 IF (CCR12 .AND. (.NOT.(CCS.OR.CIS))) THEN 2715 IOPT = 32 2716 KT12 = KT2 + NT2AM(ISYM) 2717 CALL CC_WRRSP(LIST,ILSTNR,ISYM,IOPT,MODEL,DUMMY,DUMMY, 2718 * VEC(KT12),WORK,LWORK) 2719 END IF 2720C 2721 END 2722*=====================================================================* 2723 SUBROUTINE PRJDRV(ISIDE,ISTATPRJ,ISYMPRJ,NVECS,NDIM,VECT,WORK) 2724* 2725* driver subroutine for projection. 2726* ISIDE = 1 : solve right response equations (A-w)R = X 2727* ISIDE =-1 : solve left response equations L(A+w) = X 2728* ISTATPRJ : index on list of the excited state used for proj. 2729* ISYMPRJ : symmetry of BOTH the vector and the eigevector. 2730* Please note: The vector we project from and the eigenvector 2731* we project away MUST have the same symmetry 2732* 2733* Sonia Coriani sept 1997. Revised march 2000 2734*=====================================================================* 2735#include "implicit.h" 2736#include "priunit.h" 2737#include "ccsdinp.h" 2738#include "leinf.h" 2739#include "ccexci.h" 2740#include "ccsdsym.h" 2741#include "cbirea.h" 2742* 2743 DIMENSION VECT(*), WORK(*) 2744 CHARACTER*10 MODFIL 2745 LOGICAL LOCDBG 2746 PARAMETER (LOCDBG = .FALSE.) 2747* 2748 IF (LOCDBG) THEN 2749 WRITE(LUPRI,*) 'PRJDRV: IPRLE = ',IPRLE 2750 WRITE(LUPRI,*) 'PRJDRV: ISIDE = ',ISIDE 2751 WRITE(LUPRI,*) 'PRJDRV: ISTATPRJ = ',ISTATPRJ 2752 WRITE(LUPRI,*) 'PRJDRV: ISYMPRJ = ',ISYMPRJ 2753 WRITE(LUPRI,*) 'PRJDRV: NCCVAR = ',NDIM 2754 WRITE(LUPRI,*) 'PRJDRV: NVECS = ',NVECS 2755 ENDIF 2756* 2757 IF (LMULBS) CALL QUIT('PRJDRV not yet adapted for CC-R12') 2758 2759 IF (CCS) THEN 2760 IOPT = 1 2761 ELSE 2762 IOPT = 3 2763 END IF 2764* 2765* set dimension of work space 2766* 2767 NCCVAR1 = NT1AM(ISYMPRJ) 2768 NCCVAR2 = NT2AM(ISYMPRJ) 2769* 2770* allocate WORK space for RE and LE 2771* 2772 KREVT1 = 1 2773 KREVT2 = KREVT1 + NCCVAR1 2774 KLEVT1 = KREVT2 + NCCVAR2 2775 KLEVT2 = KLEVT1 + NCCVAR1 2776* 2777* puntatore (offset) 2778* 2779 IOFFSTP = ISTATPRJ 2780* 2781* get right eigenvector RE 2782* 2783 CALL CC_RDRSP('RE',IOFFSTP,ISYMPRJ,IOPT,MODFIL, 2784 * WORK(KREVT1),WORK(KREVT2)) 2785* 2786* get left eigenvector LE 2787* 2788 CALL CC_RDRSP('LE',IOFFSTP,ISYMPRJ,IOPT,MODFIL, 2789 * WORK(KLEVT1),WORK(KLEVT2)) 2790* 2791* check dimension matching 2792* 2793 IF (CCS) THEN 2794 NDIMV = NCCVAR1 2795 ELSE 2796 NDIMV = NCCVAR1+NCCVAR2 2797 ENDIF 2798 IF (NDIMV.NE.NDIM) THEN 2799 WRITE(LUPRI,*)' Warning PRJDRV: vects dimension inconsistent' 2800 CALL QUIT('PRJDRV: NDIMV.ne.NDIM') 2801 ENDIF 2802c 2803 DO IVC = 1, NVECS 2804 IF (ISIDE.EQ.1) THEN 2805* 2806* use LE for constant, RE for subtraction 2807* 2808 KOFFV = (IVC-1)*NDIMV+1 !offset vector to be projected 2809 CALL PRJECT(NDIMV,VECT(KOFFV),WORK(KLEVT1),WORK(KREVT1)) 2810* 2811 CHECKL = DDOT(NDIMV,VECT(KOFFV),1,WORK(KLEVT1),1) 2812 CHECKR = DDOT(NDIMV,VECT(KOFFV),1,WORK(KREVT1),1) 2813 IF (LOCDBG) THEN 2814 WRITE(LUPRI,*) 'After PRJECT from right: <LE|VP> = ', CHECKL 2815 WRITE(LUPRI,*) 'After PRJECT from right: <RET|VP> = ',CHECKR 2816 END IF 2817* 2818 ELSE 2819* 2820* use RE for the constant, LE for subtraction 2821* 2822 KOFFV = (IVC-1)*NDIMV+1 !offset vector to be projected 2823 CALL PRJECT(NDIMV,VECT(KOFFV),WORK(KREVT1),WORK(KLEVT1)) 2824C 2825 CHECKR = DDOT(NDIMV,VECT(KOFFV),1,WORK(KREVT1),1) 2826 CHECKL = DDOT(NDIMV,VECT(KOFFV),1,WORK(KLEVT1),1) 2827ccs CALL PRJECT(NDIMV,VECT,WORK(KREVT1),WORK(KLEVT1)) 2828ccs CHECKR = DDOT(NDIMV,VECT,1,WORK(KREVT1),1) 2829ccs CHECKL = DDOT(NDIMV,VECT,1,WORK(KLEVT1),1) 2830 IF (LOCDBG) THEN 2831 WRITE(LUPRI,*) 'After PRJECT from left: <VP|RE> = ', CHECKR 2832 WRITE(LUPRI,*) 'After PRJECT from left: <VP|LET> = ',CHECKL 2833 END IF 2834 ENDIF 2835 END DO 2836 2837 RETURN 2838 END 2839*================================================================* 2840 SUBROUTINE PRJECT(NDIMVPR,VECT,CVECT,SVECT) 2841*================================================================* 2842* Sonia Coriani sept. 1999 2843* subroutine for projecting out eigenstates 2844* NDIMVPR = dimension of the vectors 2845* VECT = input vector to project/output projected vector 2846* CVECT = left/right vector for the constant 2847* SVECT = right/left vector to be subtracted 2848*================================================================* 2849#include "implicit.h" 2850#include "priunit.h" 2851#include "ccsdinp.h" 2852* 2853 DIMENSION VECT(NDIMVPR),CVECT(NDIMVPR),SVECT(NDIMVPR) 2854 LOGICAL LOCDBG 2855 PARAMETER (LOCDBG = .FALSE.) 2856* 2857* compute the constant for projection 2858* const = <left|vect> or const = <vect|right> 2859* 2860 CONST = - DDOT(NDIMVPR,VECT,1,CVECT,1) 2861 IF (LOCDBG) THEN 2862 WRITE(LUPRI,*) '=== We are inside PRJECT === ' 2863 WRITE(LUPRI,*) 'constant for projection = ', CONST 2864 END IF 2865* 2866* get projected vector 2867* |vectp> = |vect> - |right>*<left|vect> 2868* or 2869* <vectp| = <vect| - <vect|right>*<left| 2870* 2871 CALL DAXPY(NDIMVPR,CONST,SVECT,1,VECT,1) 2872* 2873 RETURN 2874 END 2875*======================================================================* 2876 SUBROUTINE CC_PRECOND(VECTORS,EIVAL,NVEC,NCCVAR,ISYMTR,TRIPLET, 2877 & WORK,LWORK) 2878 IMPLICIT NONE 2879#include "priunit.h" 2880#include "ccorb.h" 2881#include "ccsdsym.h" 2882#include "ccsdinp.h" 2883#include "dummy.h" 2884#include "r12int.h" 2885 2886 LOGICAL TRIPLET, LCCEQ 2887 INTEGER NVEC, NCCVAR, ISYMTR, LWORK 2888 2889#if defined (SYS_CRAY) 2890 REAL VECTORS(NCCVAR,NVEC), WORK(*), DTEST, 2891 & EIVAL(NVEC) 2892#else 2893 DOUBLE PRECISION VECTORS(NCCVAR,NVEC), WORK(*), DTEST, 2894 & EIVAL(NVEC) 2895#endif 2896 PARAMETER ( DTEST = 1.0D-04 ) 2897 2898 INTEGER KDIA, KEND1, KSCR, LWRK1, NCCVAR0, JVEC 2899 2900 KDIA = 1 2901 KEND1 = KDIA + NCCVAR 2902 2903 IF (CCR12) THEN 2904 KSCR = KEND1 2905 KEND1= KSCR + NTR12AM(ISYMTR) 2906 END IF 2907 2908 LWRK1 = LWORK - KEND1 2909 IF (LWRK1.LT.0) THEN 2910 CALL QUIT('Too little work space in CC_PRECOND') 2911 END IF 2912 2913 CALL CCLR_DIA(WORK(KDIA),ISYMTR,TRIPLET,WORK(KEND1),LWRK1) 2914C CALL LEDIAG(NCCVAR,WORK(KDIA),DTEST) 2915 2916 IF (CCS) THEN 2917 NCCVAR0 = NT1AM(ISYMTR) 2918 ELSE 2919 NCCVAR0 = NT1AM(ISYMTR) + NT2AM(ISYMTR) 2920 IF (TRIPLET) NCCVAR0 = NCCVAR0 + NT2AMA(ISYMTR) 2921 END IF 2922 2923 DO JVEC = 1,NVEC 2924 ! for conventional part divide by orbital energy differences 2925 DO I=1,NCCVAR0 2926 VECTORS(I,JVEC) = VECTORS(I,JVEC) 2927 & * SAFE_DENOM(WORK(KDIA-1+I)-EIVAL(JVEC)) 2928 END DO 2929 2930 IF (CCR12) THEN 2931 ! for R12 part divide by B matrix 2932 LCCEQ = .FALSE. 2933 CALL CC_R12NXTAM(VECTORS(NCCVAR0+1,JVEC),ISYMTR, 2934 & WORK(KSCR),LCCEQ, 2935 & DUMMY,WORK(KEND1),LWRK1) 2936 CALL DSCAL(NTR12AM(ISYMTR),-1.0D0,VECTORS(NCCVAR0+1,JVEC),1) 2937 END IF 2938 END DO 2939 2940 RETURN 2941 2942 CONTAINS 2943 2944 PURE FUNCTION SAFE_DENOM(A) 2945 DOUBLE PRECISION, PARAMETER :: ONE = 1.0D0 2946 DOUBLE PRECISION, INTENT(IN) :: A 2947 DOUBLE PRECISION :: SAFE_DENOM 2948 2949 SAFE_DENOM = SIGN(ONE,A)/MAX(ABS(A),DTEST) 2950 RETURN 2951 END FUNCTION 2952 2953 END 2954*======================================================================* 2955