1 SUBROUTINE MV10(CCIN,CCOUT) 2* 3* Initial outer routine for cc Jacobian times vector 4* 5* Jeppe Olsen, June 2000 6* 7c INCLUDE 'implicit.inc' 8c INCLUDE 'mxpdim.inc' 9 INCLUDE 'wrkspc.inc' 10 INCLUDE 'glbbas.inc' 11 INCLUDE 'ctcc.inc' 12 INCLUDE 'lorr.inc' 13 INCLUDE 'crun.inc' 14C COMMON/LORR/L_OR_R 15* 16C JAC_T_VEC(L_OR_R,CC_AMP,JAC_VEC,TVEC,VEC1,VEC2, 17C & CCVEC) 18 CALL JAC_T_VEC(L_OR_R,WORK(KCC1),CCOUT,CCIN, 19 & WORK(KVEC1P),WORK(KVEC2P),WORK(KCC5)) 20* 21 NTEST = 00 22 IF(NTEST.GE.100) THEN 23 WRITE(6,*) 24 WRITE(6,*) ' =============================== ' 25 WRITE(6,*) ' Input and output from JAC_T_VEC ' 26 WRITE(6,*) ' =============================== ' 27 WRITE(6,*) 28 WRITE(6,*) ' L_OR_R: ', L_OR_R 29 CALL WRTMAT(CCIN,1,N_CC_AMP,1,LEN_T_VEC_MX) 30 WRITE(6,*) 31 CALL WRTMAT(CCOUT,1,N_CC_AMP,1,LEN_T_VEC_MX) 32 END IF 33* 34 RETURN 35 END 36 SUBROUTINE CC_EXC_E(CCVEC1,CCVEC2,VEC1,VEC2, 37 & LUCCAMP,IEXC_RESTRT) 38* 39* Master routine for coupled cluster linear response calculation 40* of excitation energies 41* 42* Output : CC-excitation operators are stored on LU_CCEXC_OP 43* 44* Jeppe Olsen, May 2000 45* 46c INCLUDE 'implicit.inc' 47c INCLUDE 'mxpdim.inc' 48 INCLUDE 'wrkspc.inc' 49 INCLUDE 'crun.inc' 50 INCLUDE 'csm.inc' 51 INCLUDE 'ctcc.inc' 52 INCLUDE 'clunit.inc' 53 INCLUDE 'cprnt.inc' 54 INCLUDE 'glbbas.inc' 55 INCLUDE 'cc_exc.inc' 56 INCLUDE 'lorr.inc' 57 INCLUDE 'opti.inc' 58*. Scratch vectors ( To show Jesper that I am on the right track ) 59 DIMENSION VEC1(*),VEC2(*),CCVEC1(*),CCVEC2(*) 60*. Local scratch, assuming atmost 1000 roots per sym 61 INTEGER IREO(1000) 62 LOGICAL DID_RHS 63* 64 IDUM = 0 65 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'CC_EXC ') 66* 67 WRITE(6,*) 68 WRITE(6,*) ' *******************************************' 69 WRITE(6,*) ' * *' 70 WRITE(6,*) ' * CCLR calculation of excitation energies *' 71 WRITE(6,*) ' * *' 72 WRITE(6,*) ' * Version of May 2000 *' 73 WRITE(6,*) ' * Revised June 2004 *' 74 WRITE(6,*) ' * *' 75 WRITE(6,*) ' *******************************************' 76 WRITE(6,*) 77 IF(IEXC_RESTRT.NE.0) THEN 78 WRITE(6,*) ' Restarted calculation ' 79 END IF 80 IEXC_RESTRT_L = IEXC_RESTRT 81* 82* Save - if required - excitation vectors from previous run 83* 84 LU_INIAMP = IOPEN_NUS('CC_INIAMP') 85 LU_SCR0 = IOPEN_NUS('CC_EXCSCR0') 86 LU_SCR1 = IOPEN_NUS('CC_EXCSCR1') 87 LU_SCR2 = IOPEN_NUS('CC_EXCSCR2') 88 89 DID_RHS = .FALSE. 90* Loop over right and left equations 91 DO ISIDE = 2, 1, -1 92 93*. should this side be done? 94 IF (IAND(ICCEX_SLEQ,ISIDE).NE.ISIDE) CYCLE 95 96 WRITE(6,'(X,"+",44("-"),"+")') 97 IF (ISIDE.EQ.2) 98 & WRITE(6,*)'| Solving right-hand-side eigenvalue problem |' 99 IF (ISIDE.EQ.1) 100 & WRITE(6,*)'| Solving left-hand-side eigenvalue problem |' 101 WRITE(6,'(X,"+",44("-"),"+")') 102 103 IF (ISIDE.EQ.2) THEN 104 LU_CCEXC_OP_R = IOPEN_NUS('CC_EXCAMP_R') 105 LU_CCEXC_OP = LU_CCEXC_OP_R 106 LU_CCEXC_OP_REST = LU_CCEXC_OP_R 107 END IF 108 IF (ISIDE.EQ.1) THEN 109 LU_CCEXC_OP_L = IOPEN_NUS('CC_EXCAMP_L') 110 LU_CCEXC_OP = LU_CCEXC_OP_L 111 LU_CCEXC_OP_REST = LU_CCEXC_OP_L 112 ! restart LHS EVP from RHS solutions 113 IF (IEXC_RESTRT.EQ.0.AND.DID_RHS) 114 & LU_CCEXC_OP_REST = LU_CCEXC_OP_R 115 END IF 116 117 IF(IEXC_RESTRT.NE.0.OR.DID_RHS) THEN 118*. Copy initial sets of excitation operators from file 119*. LU_CCEXC_OP to LU_INIAMP 120 IEXC_RESTRT_L = 1 121 122 IF (IEXC_RESTRT.EQ.0.AND.DID_RHS) 123 & WRITE(6,*) ' Initializing with right-hand-side amplitudes ...' 124 125 CALL REWINO(LU_CCEXC_OP_REST) 126 CALL REWINO(LU_INIAMP) 127* 128 DO ISM = 1, NSMST 129 IF(NEXC_PER_SYM(ISM).NE.0) THEN 130C IMSCOMB_CC = 0 131 CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,ISM, 132 & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 133 & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VEC, 134 & MSCOMB_CC,MX_SBSTR, 135 & WORK(KISPX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 136 & NTCONF,IPRCC) 137C FRMDSC(ARRAY,NDIM,MBLOCK,IFILE,IMZERO,I_AM_PACKED) 138 NROOT_CC = NEXC_PER_SYM(ISM) 139 CALL IFRMDS(LEN_T_VEC_PREV,1,-1,LU_CCEXC_OP_REST) 140 DO IROOT_CC = 1, NROOT_CC 141 ZERO = 0.0D0 142 CALL SETVEC(CCVEC1,ZERO,LEN_T_VEC) 143 CALL FRMDSC(CCVEC1,LEN_T_VEC_PREV,-1,LU_CCEXC_OP_REST, 144 & IMZERO,I_AM_PACKED) 145 CALL TODSC(CCVEC1,LEN_T_VEC,-1,LU_INIAMP) 146 END DO 147 END IF 148 END DO 149 END IF 150* ^ End of this run is a restart 151 152 CALL REWINO(LU_CCEXC_OP) 153 IF(IEXC_RESTRT_L.NE.0) CALL REWINO(LU_INIAMP) 154 DO ISM = 1, NSMST 155 IF(NEXC_PER_SYM(ISM).NE.0) THEN 156 WRITE(6,*) 157 WRITE(6,*) ' =============================================' 158 WRITE(6,*) ' Information for excitations of symmetry',ISM 159 WRITE(6,*) ' =============================================' 160 WRITE(6,*) 161 WRITE(6,*) ' Number of roots required : ', NEXC_PER_SYM(ISM) 162 NROOT_CC = NEXC_PER_SYM(ISM) 163 ITEX_SM = ISM 164* 165*. Number of Amplitudes for this symmetry 166* 167C IMSCOMB_CC = 0 168 CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,ISM, 169 & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 170 & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VEC, 171 & MSCOMB_CC,MX_SBSTR, 172 & WORK(KISPX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 173 & NTCONF,IPRCC) 174 WRITE(6,*) ' Number of CC amplitudes ', LEN_T_VEC 175 IF (LEN_T_VEC.GT.N_CC_AMP) STOP 'DIMENSIONS' 176* 177*. Set up Diagonal for this symmetry 178* 179 IMOD = 2 ! use alpha/beta Fock-matrix 180 CALL GENCC_F_DIAG_M(IMOD,WORK(KLSOBEX),NSPOBEX_TP,CCVEC1,ISM, 181 & XDUM,IDUM,IDUM,0, 182 & VEC1,VEC2,MX_ST_TSOSO,MX_ST_TSOSO_BLK) 183*. Save on LUDIA for future use 184 CALL REWINO(LUDIA) 185 CALL TODSC(CCVEC1,LEN_T_VEC,-1,LUDIA) 186* 187*. Initialization : Restart or lowest diagonal elements 188* 189* 190 IF(IEXC_RESTRT_L.EQ.0) THEN 191*. Lowest elements of diagonal 192 CALL DUMSORT(CCVEC1,LEN_T_VEC,NROOT_CC,IREO,CCVEC2) 193*. Save corresponding initial guesses 194 CALL REWINO(LU_SCR0) 195 DO IROOT = 1, NROOT_CC 196 ZERO = 0.0D0 197 CALL SETVEC(CCVEC1,ZERO,LEN_T_VEC) 198 CCVEC1(IREO(IROOT)) = 1.0D0 199 CALL TODSC(CCVEC1,LEN_T_VEC,-1,LU_SCR0) 200C? WRITE(6,*) ' Initial vector for IROOT = ', IROOT 201C? CALL WRTMAT(CCVEC1,1,LEN_T_VEC,1,LEN_T_VEC) 202 END DO 203 ELSE IF(IEXC_RESTRT_L.NE.0) THEN 204 CALL REWINO(LU_SCR0) 205*. Read initial vectors from LU_INIAMP 206 DO IROOT =1, NROOT_CC 207 CALL FRMDSC(CCVEC1,LEN_T_VEC,-1,LU_INIAMP, 208 & IMZERO,I_AM_PACKED) 209 CALL TODSC(CCVEC1,LEN_T_VEC,-1,LU_SCR0) 210 END DO 211 END IF 212* ^ End of shift between different forms of initialization 213*. Allocate scratch space for the diagonalization 214*. Length of APROJ : MAXVEC**2 215*. Length of AVEC : MAXVEC**2 216* Length of WORK : 4*MAXVEC**2 + 7*MAXVEC 217* Length of H0SCR : 2*(NP1+NP2) ** 2 + 4 * (NP1+NP2+NQ) 218 MAXVEC = MXCIV*NROOT_CC 219 CALL MEMMAN(KLAPROJ,MAXVEC**2,'ADDL ',2,'APROJ ') 220 CALL MEMMAN(KLAVEC ,MAXVEC**2,'ADDL ',2,'AVEC ') 221 LWORK = 4*MAXVEC**2 + 7*MAXVEC 222 CALL MEMMAN(KLWORK,LWORK,'ADDL ',2,'LWORK ') 223 CALL MEMMAN(KLSCR ,0,'ADDL ',2,'LSCR ') 224 LEN = NROOT_CC*MAXIT 225 CALL MEMMAN(KLRNRM,LEN,'ADDL ',2,'RNRM ') 226 CALL MEMMAN(KLEIG ,LEN,'ADDL ',2,'EIG ') 227 CALL MEMMAN(KLFINEIG,NROOT_CC,'ADDL ',2,'FINEIG') 228*. And then to the work 229C GENMAT_DIAG(VEC1,VEC2,LU1,LU2,RNRM,EIG,FINEIG,MAXIT,NVAR, 230C & LU3,LUDIA,NROOT,MAXVEC,NINVEC, 231C & APROJ,AVEC,WORK,IPRT, 232C & NPRDIM,H0,IPNTR,NP1,NP2,NQ,H0SCR,EIGSHF, 233C & IOLSEN,IPICO) 234*. We will obtain right eigenvectors. Tell program 235*. to perform Jacobian times right vectors 236 L_OR_R = ISIDE 237 IPRDIA_L = 10 238 CONV = CCEX_CONV 239 240 CALL GENMAT_DIAG(L_OR_R, 241 & CCVEC1,CCVEC2,VEC1,VEC2, 242 & LU_SCR0,LU_SCR1,LUCCAMP, 243 & WORK(KLRNRM),WORK(KLEIG),WORK(KLFINEIG),MAXIT,LEN_T_VEC, 244 & N_CC_AMP, 245 & LU_SCR2,LUDIA,NROOT_CC,MAXVEC,NROOT_CC, 246 & WORK(KLAPROJ),WORK(KLAVEC),WORK(KLWORK),IPRDIA_L, 247 & 0,0.0D0,0,0,0,0,0.0D0,0.0D0, 248 & 0,0,CONV) 249* 250* Analyze excitation vectors - and copy to LU_CCEXC 251* 252 CALL REWINO(LU_SCR0) 253 CALL ITODS(LEN_T_VEC,1,-1,LU_CCEXC_OP) 254 DO IROOT = 1, NROOT_CC 255 WRITE(6,'(/,X,"+",75("="),"+")') 256 WRITE(6,'(X,"|",2X,A,I2,A,I4,39X,"|")') 257 & 'Analysis for symmetry ',ISM,' root ',IROOT 258 WRITE(6,'(X,"|",2X,A,F20.12,34X,"|")') 259 & 'Excitation energy: ',WORK(KLFINEIG+IROOT-1) 260 WRITE(6,'(X,"|",2X,A,F20.12,34X,"|")') 261 & 'Residual norm: ',WORK(KLRNRM+IROOT-1) 262 WRITE(6,'(X,"+",75("="),"+")') 263 CALL FRMDSC(CCVEC1,LEN_T_VEC,-1,LU_SCR0,IMZERO,IAMPACK) 264 CALL ANA_GENCC(CCVEC1,ISM) 265 CALL TODSC(CCVEC1,LEN_T_VEC,-1,LU_CCEXC_OP) 266 END DO 267* 268 END IF 269 END DO 270* 271 IF (ISIDE.EQ.2) DID_RHS = .TRUE. 272 273 END DO ! ISIDE 274* 275 IF (IAND(1,ICCEX_SLEQ).EQ.1) CALL RELUNIT(LU_CCEXC_OP_L,'KEEP') 276 IF (IAND(2,ICCEX_SLEQ).EQ.2) CALL RELUNIT(LU_CCEXC_OP_R,'KEEP') 277 CALL RELUNIT(LU_INIAMP,'DELETE') 278 CALL RELUNIT(LU_SCR0,'DELETE') 279 CALL RELUNIT(LU_SCR1,'DELETE') 280 CALL RELUNIT(LU_SCR2,'DELETE') 281* 282 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'CC_EXC ') 283 RETURN 284 END 285 286 SUBROUTINE GET_NONSYM_SUBMAT(LUC,LUHC,NVEC,NVECP,SUBMAT,NDIM, 287 & VEC1,VEC2,SUBMATP) 288* 289* Obtain subspace matrix from trialvectors and matrix times 290* trial vectors. No assumption about symmetry of matrix 291* 292* Jeppe Olsen, May 2000 293* 294* Version assuming the ability to hold two vectors in core 295* 296 INCLUDE 'implicit.inc' 297 REAL*8 INPROD 298*. Input : Previous subspace matrix 299 DIMENSION SUBMATP(*) 300*. Output : New subspace matrix 301 DIMENSION SUBMAT(*) 302*. Scratch 303 DIMENSION VEC1(*),VEC2(*) 304* 305 LBLK = -1 306 NTEST = 00 307 IF(NTEST.GE.100) THEN 308* 309 WRITE(6,*) ' Input vectors to GET_NONSYM_SUBMAT ' 310 WRITE(6,*) ' ===================================' 311* 312 WRITE(6,*) ' C vectors ' 313 CALL REWINO(LUC) 314 DO IVEC = 1, NVEC 315 CALL FRMDSC(VEC1,NDIM,LBLK,LUC,IMZERO,IAMPACK) 316 WRITE(6,*) ' Input C vector number = ', IVEC 317 CALL WRTMAT(VEC1,1,NDIM,1,NDIM) 318 END DO 319* 320 WRITE(6,*) ' HC vectors ' 321 CALL REWINO(LUHC) 322 DO IVEC = 1, NVEC 323 CALL FRMDSC(VEC1,NDIM,LBLK,LUHC,IMZERO,IAMPACK) 324 WRITE(6,*) ' Input HC vector number = ', IVEC 325 CALL WRTMAT(VEC1,1,NDIM,1,NDIM) 326 END DO 327* 328 END IF 329* 330 IMZERO = 0 331 IAMPACK = 0 332*.Reform previous matrix from (NVECP,NVECP) to (NVEC,NVEC) format 333 IF(NVECP.GT.0) THEN 334C COPMAT(AIN,AOUT,NIN,NOUT) 335 CALL COPMAT(SUBMATP,SUBMAT,NVECP,NVEC) 336 END IF 337* 338*. Add new columns of subspace matrix 339* 340C? WRITE(6,*) ' NVEC, NVECP = ', NVEC,NVECP 341 CALL REWINO(LUHC) 342 DO IVEC = 1, NVECP 343 CALL FRMDSC(VEC1,NDIM,LBLK,LUHC,IMZERO,IAMPACK) 344 END DO 345 DO J = 1, NVEC-NVECP 346*. Read in Hc(nvecp+j) 347C? WRITE(6,*) ' J = ', J 348 CALL FRMDSC(VEC1,NDIM,LBLK,LUHC,IMZERO,IAMPACK) 349 CALL REWINO(LUC) 350 DO I = 1, NVEC 351C? WRITE(6,*) ' I = ', I 352*. Read in c(i) 353 CALL FRMDSC(VEC2,NDIM,LBLK,LUC,IMZERO,IAMPACK) 354* Submat(i,nvecp+j) = c(i)t Hc(nvecp+j) 355 SUBMAT((NVECP+J-1)*NVEC+I) = INPROD(VEC2,VEC1,NDIM) 356 END DO 357 END DO 358* 359*. Add new rows of subspace matrix 360* 361C? WRITE(6,*) ' Part 2 ' 362 CALL REWINO(LUC) 363 DO I = 1, NVECP 364 CALL FRMDSC(VEC2,NDIM,LBLK,LUC,IMZERO,IAMPACK) 365 END DO 366C CALL SKPRC3(NVECP,LUC) 367 DO I = 1, NVEC-NVECP 368C? WRITE(6,*) ' I = ', I 369*. Read in c(nvecp+i) 370 CALL FRMDSC(VEC2,NDIM,LBLK,LUC,IMZERO,IAMPACK) 371 CALL REWINO(LUHC) 372 DO J = 1, NVECP 373C? WRITE(6,*) ' J = ', J 374*. Read in Hc(j) 375 CALL FRMDSC(VEC1,NDIM,LBLK,LUHC,IMZERO,IAMPACK) 376* Submat(i,nvecp+j) = c(i)t Hc(nvecp+j) 377 SUBMAT((J-1)*NVEC+I+NVECP) = INPROD(VEC2,VEC1,NDIM) 378 END DO 379 END DO 380* 381 IF(NTEST.GE.100) THEN 382 WRITE(6,*) ' Updated subspace matrix ' 383 CALL WRTMAT(SUBMAT,NVEC,NVEC,NVEC,NVEC) 384 END IF 385* 386 RETURN 387 END 388 SUBROUTINE GENMAT_DIAG(L_OR_R, 389 & VEC1,VEC2,SCRVEC1,SCRVEC2, 390 & LU1,LU2,LUCCAMP,!LUINTM1,LUINTM2,LUINTM3, 391 & RNRM,EIG,FINEIG,MAXIT,NVAR,NVAR0, 392 & LU3,LUDIA,NROOT,MAXVEC,NINVEC, 393 & APROJ,AVEC,WORK,IPRT, 394 & NPRDIM,H0,IPNTR,NP1,NP2,NQ,H0SCR,EIGSHF, 395 & IOLSEN,IPICO,TEST) 396* 397* p.t. matrix vector routine is hardwired 398* 399* Iterative solution of eigenvalue problem for general matrix 400* 401* Final and intermediate eigenvalues and eigenvectors are assumed real 402* 403* MIN version requiring two vectors in core 404* 405* 406* Jeppe Olsen May 2000, from MINDV4 407* 408* Input : 409* ======= 410* LU1 : Initial set of vectors 411* VEC1,VEC2 : Two vectors, each must be dimensioned to hold 412* complete vector 413* LU2,LU3 : Scatch files 414* LUDIA : File containing diagonal of matrix 415* NROOT : Number of eigenvectors to be obtained 416* MAXVEC : Largest allowed number of vectors 417* must atleast be 2 * NROOT 418* NINVEC : Number of initial vectors ( atleast NROOT ) 419* NPRDIM : Dimension of subspace with 420* nondiagonal preconditioning 421* (NPRDIM = 0 indicates no such subspace ) 422* For NPRDIM .gt. 0: 423* PEIGVC : EIGENVECTORS OF MATRIX IN PRIMAR SPACE 424* Holds preconditioner matrices 425* PHP,PHQ,QHQ in this order !! 426* PEIGVL : EIGENVALUES OF MATRIX IN PRIMAR SPACE 427* IPNTR : IPNTR(I) IS ORIGINAL ADRESS OF SUBSPACE ELEMENT I 428* NP1,NP2,NQ : Dimension of the three subspaces 429* 430* H0SCR : Scratch space for handling H0, at least 2*(NP1+NP2) ** 2 + 431* 4 (NP1+NP2+NQ) 432* On input LU1 is supposed to hold initial guess to eigenvectors 433* 434* IOLSEN : Use inverse iteration modified Davidson 435* 436 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 437 LOGICAL CONVER 438 DIMENSION VEC1(*),VEC2(*) 439 REAL * 8 INPROD 440 DIMENSION RNRM(MAXIT,NROOT),EIG(MAXIT,NROOT) 441 DIMENSION FINEIG(1) 442 DIMENSION H0(*),IPNTR(*) 443* 444*. Scratch through argument list 445* 446*. Length of APROJ : MAXVEC**2 447*. Length of AVEC : MAXVEC**2 448* Length of WORK : 4*MAXVEC**2 + 7*MAXVEC 449* Length of H0SCR : 2*(NP1+NP2) ** 2 + 4 * (NP1+NP2+NQ) 450 DIMENSION APROJ(*),AVEC(*),WORK(*) 451 DIMENSION H0SCR(*) 452*. Local scratch - atmost 1000 roots 453 LOGICAL RTCNV(1000) 454*. Ad hoc arrays for root hooming 455 DIMENSION XJEP(3000),IXJEP(3000) 456* 457* 458 CALL ATIM(CPU0,WALL0) 459* 460 IOLSTM = 0 461 IF(IPRT.GT.1.AND.(IOLSEN.NE.0.AND.IPICO.EQ.0)) 462 & WRITE(6,*) ' Inverse iteration modified Davidson, Variational' 463 IF(IPRT.GT.1.AND.(IOLSEN.NE.0.AND.IPICO.NE.0)) 464 & WRITE(6,*) ' Inverse iteration modified Davidson, Perturbational' 465 IF(IPRT.GT.1.AND.(IOLSEN.EQ.0.AND.IPICO.EQ.0)) 466 & WRITE(6,*) ' Normal Davidson, Variational ' 467 IF(IPRT.GT.1.AND.(IOLSEN.EQ.0.AND.IPICO.NE.0)) 468 & WRITE(6,*) ' Normal Davidson, Perturbational' 469C IF( MAXVEC .LT. 2 * NROOT ) THEN 470C WRITE(6,*) ' SORRY MINDV2 WOUNDED , MAXVEC .LT. 2*NROOT ' 471C STOP ' ENFORCED STOP IN MINDV2' 472C END IF 473 WRITE(6,*) ' Convergence threshold for residual : ', TEST 474* 475 I_DO_ROOTHOOMING = 0 476 IF(I_DO_ROOTHOOMING.EQ.1) THEN 477 WRITE(6,*) ' Root hooming active ' 478 END IF 479* 480*. Scratch files 481 LUINTM1 = IOPEN_NUS('CCGMATSCR1') 482 LUINTM2 = IOPEN_NUS('CCGMATSCR2') 483 LUINTM3 = IOPEN_NUS('CCGMATSCR3') 484* 485 IF(IPICO.NE.0) THEN 486 MAXVEC = 2*NROOT 487 END IF 488*. Division of scratch memory 489 KAPROJ = 1 490 KFREE = KAPROJ + MAXVEC**2 491* 492 KARVAL = KFREE 493 KFREE = KARVAL + MAXVEC 494* 495 KAIVAL = KFREE 496 KFREE = KAIVAL + MAXVEC 497* 498 KARVEC = KFREE 499 KFREE = KARVEC + MAXVEC**2 500* 501 KAIVEC = KFREE 502 KFREE = KAIVEC + MAXVEC**2 503* 504 KZ = KFREE 505 KFREE = KZ + MAXVEC**2 506* 507 KW = KFREE 508 KFREE = KW + MAXVEC 509* 510 KSCR1 = KFREE 511 KFREE = KSCR1 + 4*MAXVEC 512* 513C TEST = 1.0D-10 514 CONVER = .FALSE. 515 DO 1234 MACRO = 1,1 516* 517*. INITAL ITERATION 518* 519* =========================================================== 520*. The initial vectors does not neccessarily constitute an 521*. orthonormal basis. Start by orthogonalizing 522* =========================================================== 523* 524 DO IROOT = 1, NINVEC 525*. Read vector IROOT in 526 CALL REWINO(LU1) 527 DO ISKP = 1, IROOT 528 CALL FRMDSC(VEC1,NVAR,-1,LU1,IMZERO,IAMPACK) 529 END DO 530*. Diagonal element 531 WORK(KAPROJ-1+(IROOT-1)*NINVEC+IROOT) = 532 & INPROD(VEC1,VEC1,NVAR) 533 DO JROOT = IROOT+1, NINVEC 534 CALL FRMDSC(VEC2,NVAR,-1,LU1,IMZERO,IAMPACK) 535 WORK(KAPROJ-1+(IROOT-1)*NINVEC+JROOT) = 536 & INPROD(VEC1,VEC2,NVAR) 537 WORK(KAPROJ-1+(JROOT-1)*NINVEC+IROOT) = 538 & WORK(KAPROJ-1+(IROOT-1)*NINVEC+JROOT) 539 END DO 540 END DO 541* 542 IF(IPRT.GE.1) THEN 543 WRITE(6,*) ' Overlap matrix of initial basis ' 544 CALL WRTMAT(WORK(KAPROJ),NROOT,NROOT,NROOT,NROOT) 545 END IF 546*. Orthonormal basis for the NROOT lowest roots 547 CALL MGS3(WORK(KARVEC),WORK(KAPROJ),NINVEC,WORK(KAIVEC)) 548*. Orthogonalize, save orthogonlized vectors on LU3 549 CALL REWINO( LU3) 550 DO IROOT = 1, NINVEC 551 CALL REWINO( LU1) 552 CALL SETVEC(VEC1,0.0D0,NVAR) 553 DO IVEC = 1, NINVEC 554 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 555 FACTOR = WORK(KARVEC-1+(IROOT-1)*NINVEC+IVEC) 556 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 557 END DO 558 CALL TODSC(VEC1,NVAR,-1 ,LU3) 559 END DO 560* 561 CALL REWINO( LU1) 562 CALL REWINO( LU3) 563 DO IVEC = 1,NINVEC 564 CALL FRMDSC(VEC1,NVAR,-1 ,LU3,IMZERO,IAMPACK) 565 CALL TODSC (VEC1,NVAR,-1, LU1) 566 END DO 567* 568 ITER = 1 569 CALL REWINO( LU1 ) 570 CALL REWINO( LU2 ) 571* Mat * initial vectors 572 DO JVEC = 1,NINVEC 573 IADD_RHS = 0 574 IF (JVEC.EQ.1) THEN 575 LUINTM1_L = -LUINTM1 576 LUINTM2_L = -LUINTM2 577 ELSE 578 LUINTM1_L = LUINTM1 579 LUINTM2_L = LUINTM2 580 END IF 581 IWRMOD=1 582 XDUM = -1 583 LUDUM = 0 584 CALL JAC_T_VEC2(L_OR_R,0,IWRMOD,0,0, 585 & VEC1,VEC2,SCRVEC1,SCRVEC2, 586 & NVAR,NVAR0, 587 & XDUM,XNORM,RNORM, 588 & LUCCAMP,LUDUM,LU1,LU2,LUDUM, 589 & LUINTM1_L,LUINTM2_L,LUINTM3) 590c CALL FRMDSC(VEC1,NVAR,-1 ,LU1,IMZERO,IAMPACK) 591c CALL MV10(VEC1,VEC2) 592c CALL TODSC(VEC2,NVAR,-1 ,LU2) 593 END DO 594*. Projected matrix 595 CALL GET_NONSYM_SUBMAT(LU1,LU2,NINVEC,0,APROJ,NVAR, 596 & VEC1,VEC2,0.0D0) 597 IF( IPRT .GE.10 ) THEN 598 WRITE(6,*) ' INITIAL PROJECTED MATRIX ' 599 CALL WRTMAT(APROJ,NINVEC,NINVEC,NINVEC,NINVEC) 600 END IF 601* Diagonalize initial subspace matrix 602 CALL COPVEC(APROJ,WORK(KAPROJ),NINVEC*NINVEC) 603 CALL EIGGMT3(WORK(KAPROJ),NINVEC,WORK(KARVAL),WORK(KAIVAL), 604 & WORK(KARVEC),WORK(KAIVEC),WORK(KZ),WORK(KW), 605 & WORK(KSCR1),1,1) 606C EIGGMT3(AMAT,NDIM,ARVAL,AIVAL,ARVEC,AIVEC, 607C & Z,W,SCR,IORD,IEIGVC) 608 DO 20 IROOT = 1, NROOT 609 EIG(1,IROOT) = WORK(KARVAL-1+IROOT) 610 20 CONTINUE 611 CALL COPVEC(WORK(KARVEC),AVEC,NINVEC**2) 612* 613 IF( IPRT .GE. 3 ) THEN 614 WRITE(6,'(A,I4)') ' Initial set of eigenvalues ' 615 WRITE(6,'(5F18.13)') 616 & ( (EIG(ITER,IROOT)+EIGSHF),IROOT=1,NROOT) 617 END IF 618 NVEC = NINVEC 619 IF (MAXIT .EQ. 1 ) GOTO 901 620* 621 WRITE(6,'(">>>",A)') 622 & ' Iter. Root exc. energy res. norm conv.' 623 WRITE(6,'(">>>",A)') 624 & '-------------------------------------------------' 625* 626** LOOP OVER ITERATIONS 627* 628 1000 CONTINUE 629 IF(IPRT .GE. 10 ) THEN 630 WRITE(6,*) ' INFO FORM ITERATION .... ', ITER 631 END IF 632 633 634 ITER = ITER + 1 635* 636** 1 NEW DIRECTION TO BE INCLUDED 637* 638* 1.1 : R = H*X - EIGAPR*X 639 IADD = 0 640 CONVER = .TRUE. 641 RNRMMX = 0d0 642 DO 100 IROOT = 1, NROOT 643 CALL SETVEC(VEC1,0.0D0,NVAR) 644* 645 CALL REWINO( LU2) 646 DO 60 IVEC = 1, NVEC 647 CALL FRMDSC(VEC2,NVAR,-1 ,LU2,IMZERO,IAMPACK) 648 FACTOR = AVEC((IROOT-1)*NVEC+IVEC) 649 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 650 60 CONTINUE 651* 652 EIGAPR = EIG(ITER-1,IROOT) 653 CALL REWINO( LU1) 654 DO 50 IVEC = 1, NVEC 655 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 656 FACTOR = (-EIGAPR)*AVEC((IROOT-1)*NVEC+ IVEC) 657 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 658 50 CONTINUE 659 IF ( IPRT .GE.600 ) THEN 660 WRITE(6,*) ' ( HX - EX ) ' 661 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 662 END IF 663* STRANGE PLACE TO TEST CONVERGENCE , BUT .... 664 RNORM = SQRT( INPROD(VEC1,VEC1,NVAR) ) 665 RNRM(ITER-1,IROOT) = RNORM 666 RNRMMX = MAX(RNRMMX,RNORM) 667 IF(RNORM.LT. TEST ) THEN 668 RTCNV(IROOT) = .TRUE. 669 ELSE 670 RTCNV(IROOT) = .FALSE. 671 CONVER = .FALSE. 672 END IF 673* 674 WRITE(6,'(">>>",5X,I5,F20.12,2X,E10.4,4X,L)') 675 & IROOT,EIGAPR,RNORM,RTCNV(IROOT) 676* 677 IF( ITER .GT. MAXIT) GOTO 100 678*. 1.2 : MULTIPLY WITH INVERSE HESSIAN APROXIMATION TO GET NEW DIRECTIO 679 IF( .NOT. RTCNV(IROOT) ) THEN 680 IADD = IADD + 1 681 CALL REWINO( LUDIA) 682 CALL FRMDSC(VEC2,NVAR,-1 ,LUDIA,IMZERO,IAMPACK) 683 CALL H0M1TV(VEC2,VEC1,VEC1,NVAR,NPRDIM,IPNTR, 684 & H0,-EIGAPR,H0SCR,XDUMMY,NP1,NP2,NQ, 685 & IPRT) 686 IF ( IPRT .GE. 600) THEN 687 WRITE(6,*) ' (D-E)-1 *( HX - EX ) ' 688 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 689 END IF 690* 691 IF(IOLSTM .NE. 0 ) THEN 692* add Olsen correction if neccessary 693 CALL REWINO(LU3) 694 CALL TODSC(VEC1,NVAR,-1,LU3) 695* Current eigen vector 696 CALL REWINO( LU1) 697 CALL SETVEC(VEC1,0.0D0,NVAR) 698 DO 59 IVEC = 1, NVEC 699 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 700 FACTOR = AVEC((IROOT-1)*NVEC+ IVEC) 701 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 702 59 CONTINUE 703 IF ( IPRT .GE. 600 ) THEN 704 WRITE(6,*) ' And X ' 705 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 706 END IF 707 CALL TODSC(VEC1,NVAR,-1,LU3) 708* (H0 - E )-1 * X 709 CALL REWINO( LUDIA) 710 CALL FRMDSC(VEC2,NVAR,-1 ,LUDIA,IMZERO,IAMPACK) 711 CALL H0M1TV(VEC2,VEC1,VEC2,NVAR,NPRDIM,IPNTR, 712 & H0,-EIGAPR,H0SCR,XDUMMY,NP1,NP2,NQ, 713 & IPRT) 714 CALL TODSC(VEC2,NVAR,-1,LU3) 715* Gamma = X(T) * (H0 - E) ** -1 * X 716 GAMMA = INPROD(VEC2,VEC1,NVAR) 717* is X an eigen vector for (H0 - 1 ) - 1 718 CALL VECSUM(VEC2,VEC1,VEC2,GAMMA,-1.0D0,NVAR) 719 VNORM = SQRT(MAX(0.0D0,INPROD(VEC2,VEC2,NVAR))) 720 IF(VNORM .GT. 1.0D-7 ) THEN 721 IOLSAC = 1 722 ELSE 723 IOLSAC = 0 724 END IF 725 IF(IOLSAC .EQ. 1 ) THEN 726 IF(IPRT.GE.5) WRITE(6,*) ' Olsen Correction active ' 727 CALL REWINO(LU3) 728 CALL FRMDSC(VEC2,NVAR,-1,LU3,IMZERO,IAMPACK) 729 DELTA = INPROD(VEC1,VEC2,NVAR) 730 CALL FRMDSC(VEC1,NVAR,-1,LU3,IMZERO,IAMPACK) 731 CALL FRMDSC(VEC1,NVAR,-1,LU3,IMZERO,IAMPACK) 732 FACTOR = (-DELTA)/GAMMA 733 IF(IPRT.GE.5) WRITE(6,*) ' DELTA,GAMMA,FACTOR' 734 IF(IPRT.GE.5) WRITE(6,*) DELTA,GAMMA,FACTOR 735 CALL VECSUM(VEC1,VEC1,VEC2,FACTOR,1.0D0,NVAR) 736 IF(IPRT.GE.600) THEN 737 WRITE(6,*) ' Modified new trial vector ' 738 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 739 END IF 740 ELSE 741 IF(IPRT.GT.0) WRITE(6,*) 742 & ' Inverse correction switched of' 743 CALL REWINO(LU3) 744 CALL FRMDSC(VEC1,NVAR,-1,LU3,IMZERO,IAMPACK) 745 END IF 746 END IF 747*. 1.3 ORTHOGONALIZE TO ALL PREVIOUS VECTORS 748 XNRMI = INPROD(VEC1,VEC1,NVAR) 749 CALL REWINO( LU1 ) 750 751 DO 80 IVEC = 1,NVEC+IADD-1 752 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 753 OVLAP = INPROD(VEC1,VEC2,NVAR) 754 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,-OVLAP,NVAR) 755 80 CONTINUE 756*. 1.4 Normalize vector and check for linear dependency 757 SCALE = INPROD(VEC1,VEC1,NVAR) 758 IF(ABS(SCALE)/XNRMI .LT. 1.0D-10) THEN 759*. Linear dependency 760 IADD = IADD - 1 761 IF ( IPRT .GE. 10 ) THEN 762 WRITE(6,*) ' Trial vector linear dependent so OUT !!! ' 763 END IF 764 ELSE 765 C1NRM = SQRT(SCALE) 766 FACTOR = 1.0D0/SQRT(SCALE) 767 CALL SCALVE(VEC1,FACTOR,NVAR) 768* 769 CALL TODSC(VEC1,NVAR,-1 ,LU1) 770 IF ( IPRT .GE.600 ) THEN 771 WRITE(6,*) 'ORTHONORMALIZED (D-E)-1 *( HX - EX ) ' 772 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 773 END IF 774 END IF 775* 776 END IF 777 100 CONTINUE 778 WRITE(6,'(">>>",I5,25X,(2X,E10.4))') 779 & ITER-1,RNRMMX 780 IF( CONVER ) GOTO 901 781 IF( ITER.GT. MAXIT) THEN 782 ITER = MAXIT 783 GOTO 1001 784 END IF 785* 786** 2 : OPTIMAL COMBINATION OF NEW AND OLD DIRECTION 787* 788* 2.1: MULTIPLY NEW DIRECTION WITH MATRIX 789 CALL REWINO( LU1) 790 CALL REWINO( LU2) 791 DO 110 IVEC = 1, NVEC 792 CALL FRMDSC(VEC1,NVAR,-1,LU1,IMZERO,IAMPACK) 793 CALL FRMDSC(VEC1,NVAR,-1,LU2,IMZERO,IAMPACK) 794 110 CONTINUE 795* 796 DO 150 IVEC = 1, IADD 797 IWRMOD = 1 798 XDUM = -1 799 LUDUM = 0 800 CALL JAC_T_VEC2(L_OR_R,0,IWRMOD,0,0, 801 & VEC1,VEC2,SCRVEC1,SCRVEC2, 802 & NVAR,NVAR0, 803 & XDUM,XNORM,RNORM, 804 & LUCCAMP,LUDUM,LU1,LU2,LUDUM, 805 & LUINTM1,LUINTM2,LUINTM3) 806c CALL FRMDSC(VEC1,NVAR,-1 ,LU1,IMZERO,IAMPACK) 807c CALL MV10(VEC1,VEC2) 808c CALL TODSC(VEC2,NVAR,-1 ,LU2) 809 150 CONTINUE 810*.Augment projected matrix 811C GET_NONSYM_SUBMAT(LUC,LUHC,NVEC,NVECP,SUBMAT,NDIM, 812C & VEC1,VEC2,SUBMATP) 813 CALL COPVEC(APROJ,WORK(KAPROJ),NVEC**2) 814 CALL GET_NONSYM_SUBMAT(LU1,LU2,NVEC+IADD,NVEC,APROJ,NVAR, 815 & VEC1,VEC2,WORK(KAPROJ)) 816* DIAGONALIZE PROJECTED MATRIX 817 NVEC = NVEC + IADD 818 CALL COPVEC(APROJ,WORK(KAPROJ),NVEC*NVEC) 819 CALL EIGGMT3(WORK(KAPROJ),NVEC,WORK(KARVAL),WORK(KAIVAL), 820 & WORK(KARVEC),WORK(KAIVEC),WORK(KZ),WORK(KW), 821 & WORK(KSCR1),1,1) 822 823 DO IROOT = 1, NROOT 824 EIG(ITER,IROOT) = WORK(KARVAL-1+IROOT) 825 END DO 826 CALL COPVEC(WORK(KARVEC),AVEC,NROOT*NVEC) 827 IF(I_DO_ROOTHOOMING.EQ.1) THEN 828* 829*. Reorder roots so the NROOT with the largest overlap with 830* the original roots become the first 831* 832*. Norm of wavefunction in previous space 833 DO IVEC = 1, NVEC 834 XJEP(IVEC) = INPROD(AVEC(1+(IVEC-1)*NROOT), 835 & AVEC(1+(IVEC-1)*NROOT),NROOT) 836 END DO 837 WRITE(6,*) 838 & ' Norm of projections to previous vector space ' 839 CALL WRTMAT(XJEP,1,NVEC,1,NVEC) 840*. My sorter arranges in increasing order, multiply with minus 1 841* so the eigenvectors with largest overlap comes out first 842 ONEM = -1.0D0 843 CALL SCALVE(XJEP,ONEM,NVEC) 844 CALL SORLOW(XJEP,XJEP(1+NVEC),IXJEP,NVEC,NVEC,NSORT,IPRT) 845 IF(NSORT.LT.NVEC) THEN 846 WRITE(6,*) ' Warning : Some elements lost in sorting ' 847 WRITE(6,*) ' NVEC,NSORT = ', NSORT,NVEC 848 END IF 849 IF(IPRT.GE.0) THEN 850 WRITE(6,*) ' New roots choosen as vectors ' 851 CALL IWRTMA(IXJEP,1,NROOT,1,NROOT) 852 END IF 853*. Reorder 854 DO INEW = 1, NVEC 855 IOLD = IXJEP(INEW) 856 CALL COPVEC(AVEC(1+(IOLD-1)*NVEC),XJEP(1+(INEW-1)*NVEC),NVEC) 857 END DO 858 CALL COPVEC(XJEP,AVEC,NVEC*NVEC) 859 DO INEW = 1, NVEC 860 IOLD = IXJEP(INEW) 861 XJEP(INEW) = WORK(KARVAL-1+IOLD) 862 END DO 863 DO INEW = 1, NROOT 864 EIG(ITER,INEW)=XJEP(INEW) 865 END DO 866* 867 IF(IPRT.GE.3) THEN 868 WRITE(6,*) ' Reordered AVEC arrays ' 869 CALL WRTMAT(AVEC,NVEC,NVEC,NVEC,NVEC) 870 END IF 871* 872 END IF 873* ^ End of root homing procedure 874 IF(IPRT .GE. 3 ) THEN 875 WRITE(6,'(A,I4)') ' Eigenvalues of iteration ..', ITER 876 WRITE(6,'(5F18.13)') 877 & ( (EIG(ITER,IROOT)+EIGSHF) ,IROOT=1,NROOT) 878 WRITE(6,*) ' Norms of residuals ' 879 WRITE(6,'(10E13.5)') (RNRM(ITER-1,KROOT),KROOT=1,NROOT) 880 END IF 881* 882 IF( IPRT .GE. 5 ) THEN 883 WRITE(6,*) ' PROJECTED MATRIX AND EIGEN PAIRS ' 884 CALL WRTMAT(AVEC,NVEC,NROOT,NVEC,NROOT) 885 WRITE(6,'(2X,E13.7)') (EIG(ITER,IROOT),IROOT = 1, NROOT) 886 END IF 887* 888** PERHAPS RESET OR ASSEMBLE CONVERGED EIGENVECTORS 889* 890 901 CONTINUE 891* 892 IPULAY = 0 893 IF(IPULAY.EQ.1 .AND. MAXVEC.EQ.3 .AND.NVEC.GE.2. 894 & .AND. .NOT.CONVER) THEN 895* Save trial vectors : 1 -- current trial vector 896* 2 -- previous trial vector orthogonalized 897 CALL REWINO( LU3) 898 CALL REWINO( LU1) 899*. Current trial vector 900 CALL SETVEC(VEC1,0.0D0,NVAR) 901 DO 2200 IVEC = 1, NVEC 902 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 903 FACTOR = AVEC(IVEC) 904 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 905 2200 CONTINUE 906 SCALE = INPROD(VEC1,VEC1,NVAR) 907 SCALE = 1.0D0/SQRT(SCALE) 908 CALL SCALVE(VEC1,SCALE,NVAR) 909 CALL TODSC(VEC1,NVAR,-1 ,LU3) 910* Previous trial vector orthonormalized 911 CALL REWINO(LU1) 912 CALL FRMDSC(VEC2,NVAR,-1,LU1,IMZERO,IAMPACK) 913 OVLAP = INPROD(VEC1,VEC2,NVAR) 914 CALL VECSUM(VEC2,VEC2,VEC1,1.0D0,-OVLAP,NVAR) 915 SCALE2 = INPROD(VEC2,VEC2,NVAR) 916 SCALE2 = 1.0D0/SQRT(SCALE2) 917 CALL SCALVE(VEC2,SCALE2,NVAR) 918 CALL TODSC(VEC2,NVAR,-1,LU3) 919* 920 CALL REWINO( LU1) 921 CALL REWINO( LU3) 922 DO 2411 IVEC = 1,2 923 CALL FRMDSC(VEC1,NVAR,-1 ,LU3,IMZERO,IAMPACK) 924 CALL TODSC (VEC1,NVAR,-1, LU1) 925 2411 CONTINUE 926*. Corresponding sigma vectors 927 CALL REWINO ( LU3) 928 CALL REWINO( LU2) 929 CALL SETVEC(VEC1,0.0D0,NVAR) 930 DO 2250 IVEC = 1, NVEC 931 CALL FRMDSC(VEC2,NVAR,-1 ,LU2,IMZERO,IAMPACK) 932 FACTOR = AVEC(IVEC) 933 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 934 2250 CONTINUE 935* 936 CALL SCALVE(VEC1,SCALE,NVAR) 937 CALL TODSC(VEC1,NVAR,-1, LU3) 938* Sigma vector corresponding to second vector on LU1 939 CALL REWINO(LU2) 940 CALL FRMDSC(VEC2,NVAR,-1,LU2,IMZERO,IAMPACK) 941 CALL VECSUM(VEC2,VEC2,VEC1,1.0D0,-OVLAP,NVAR) 942 CALL SCALVE(VEC2,SCALE2,NVAR) 943 CALL TODSC(VEC2,NVAR,-1,LU3) 944* 945 CALL REWINO( LU2) 946 CALL REWINO( LU3) 947 DO 2400 IVEC = 1,2 948 CALL FRMDSC(VEC2,NVAR,-1 ,LU3,IMZERO,IAMPACK) 949 CALL TODSC (VEC2,NVAR,-1 ,LU2) 950 2400 CONTINUE 951 NVEC = 2 952* 953 CALL SETVEC(AVEC,0.0D0,NVEC**2) 954 DO 2410 IROOT = 1,NVEC 955 AVEC((IROOT-1)*NVEC+IROOT) = 1.0D0 956 2410 CONTINUE 957*.Projected hamiltonian 958 CALL REWINO( LU1 ) 959 DO 2010 IVEC = 1,NVEC 960 CALL FRMDSC(VEC1,NVAR,-1 ,LU1,IMZERO,IAMPACK) 961 CALL REWINO( LU2) 962 DO 2008 JVEC = 1, IVEC 963 CALL FRMDSC(VEC2,NVAR,-1 ,LU2,IMZERO,IAMPACK) 964 IJ = IVEC*(IVEC-1)/2 + JVEC 965 APROJ(IJ) = INPROD(VEC1,VEC2,NVAR) 966 2008 CONTINUE 967 2010 CONTINUE 968 END IF 969 IF(NVEC+NROOT.GT.MAXVEC .OR. CONVER .OR. MAXIT .EQ.ITER)THEN 970*. Select space spanning lowest NROOT as new subspace 971*. Note that subspace is required to be orthogonal although 972*. the eigenvectors in general not are orthogonal 973* 974*. Overlap matrix of lowest NROOT eigenvectors 975 DO IROOT = 1, NROOT 976 DO JROOT = 1, NROOT 977 WORK(KAPROJ-1+(IROOT-1)*NROOT+JROOT) = 978 & INPROD(AVEC((IROOT-1)*NVEC+1), 979 & AVEC((JROOT-1)*NVEC+1),NVEC) 980 END DO 981 END DO 982* 983 IF(IPRT.GE.1) THEN 984 WRITE(6,*) ' Overlap matrix in new basis ' 985 CALL WRTMAT(WORK(KAPROJ),NROOT,NROOT,NROOT,NROOT) 986 END IF 987*. Orthonormal basis for the NROOT lowest roots 988C MGS3(X,S,NDIM,SCR1) 989 CALL MGS3(WORK(KARVEC),WORK(KAPROJ),NROOT,WORK(KAIVEC)) 990*. Orthogonalization matrix is now in WORK(KARVEC) 991* Obtain new basis - if iteration procedure continues 992 IF(ITER.LT.MAXIT .AND. (.NOT.CONVER)) THEN 993 CALL MATML4(WORK(KAIVEC),AVEC,WORK(KARVEC), 994 & NVEC,NROOT,NVEC,NROOT,NROOT,NROOT,0) 995 CALL COPVEC(WORK(KAIVEC),AVEC,NROOT*NVEC) 996 END IF 997* 998 CALL REWINO( LU3) 999 DO 320 IROOT = 1, NROOT 1000 CALL REWINO( LU1) 1001 CALL SETVEC(VEC1,0.0D0,NVAR) 1002 DO 200 IVEC = 1, NVEC 1003 CALL FRMDSC(VEC2,NVAR,-1 ,LU1,IMZERO,IAMPACK) 1004 FACTOR = AVEC((IROOT-1)*NVEC+IVEC) 1005 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 1006 200 CONTINUE 1007* 1008 SCALE = INPROD(VEC1,VEC1,NVAR) 1009 SCALE = 1.0D0/SQRT(SCALE) 1010 CALL SCALVE(VEC1,SCALE,NVAR) 1011 CALL TODSC(VEC1,NVAR,-1 ,LU3) 1012 320 CONTINUE 1013 CALL REWINO( LU1) 1014 CALL REWINO( LU3) 1015 DO 411 IVEC = 1,NROOT 1016 CALL FRMDSC(VEC1,NVAR,-1 ,LU3,IMZERO,IAMPACK) 1017 CALL TODSC (VEC1,NVAR,-1, LU1) 1018 411 CONTINUE 1019* CORRESPONDING SIGMA VECTOR 1020 CALL REWINO ( LU3) 1021 DO 329 IROOT = 1, NROOT 1022 CALL REWINO( LU2) 1023 CALL SETVEC(VEC1,0.0D0,NVAR) 1024 DO 250 IVEC = 1, NVEC 1025 CALL FRMDSC(VEC2,NVAR,-1 ,LU2,IMZERO,IAMPACK) 1026 FACTOR = AVEC((IROOT-1)*NVEC+IVEC) 1027 CALL VECSUM(VEC1,VEC1,VEC2,1.0D0,FACTOR,NVAR) 1028 250 CONTINUE 1029* 1030 CALL SCALVE(VEC1,SCALE,NVAR) 1031 CALL TODSC(VEC1,NVAR,-1, LU3) 1032 329 CONTINUE 1033* PLACE C IN LU1 AND HC IN LU2 1034 CALL REWINO( LU2) 1035 CALL REWINO( LU3) 1036 DO 400 IVEC = 1,NROOT 1037 CALL FRMDSC(VEC2,NVAR,-1 ,LU3,IMZERO,IAMPACK) 1038 CALL TODSC (VEC2,NVAR,-1 ,LU2) 1039 400 CONTINUE 1040 NVEC = NROOT 1041* 1042 CALL SETVEC(AVEC,0.0D0,NVEC**2) 1043 DO 410 IROOT = 1,NROOT 1044 AVEC((IROOT-1)*NROOT+IROOT) = 1.0D0 1045 410 CONTINUE 1046* 1047 CALL GET_NONSYM_SUBMAT(LU1,LU2,NROOT,0,APROJ,NVAR, 1048 & VEC1,VEC2,WORK(KAPROJ)) 1049C 1050 END IF 1051C 1052C IF( ITER .LT. MAXIT .AND. .NOT. CONVER) GOTO 1000 1053 IF( ITER .LE. MAXIT .AND. .NOT. CONVER) GOTO 1000 1054 1001 CONTINUE 1055*. Place first eigenvector in vec1 1056 CALL REWINO(LU1) 1057 CALL FRMDSC(VEC1,NVAR,-1 ,LU1,IMZERO,IAMPACK) 1058 1059 1060 1061* ( End of loop over iterations ) 1062* 1063* 1064* 1065 IF( .NOT. CONVER ) THEN 1066* CONVERGENCE WAS NOT OBTAINED 1067 IF(IPRT .GE. 2 ) 1068 & WRITE(6,1170) MAXIT 1069 1170 FORMAT('0 Convergence was not obtained in ',I3,' iterations') 1070 ELSE 1071* CONVERGENCE WAS OBTAINED 1072 ITER = ITER - 1 1073 IF (IPRT .GE. 2 ) 1074 & WRITE(6,1180) ITER 1075 1180 FORMAT(1H0,' Convergence was obtained in ',I3,' iterations') 1076 END IF 1077*. Final eigenvalues 1078 DO 1601 IROOT = 1, NROOT 1079 FINEIG(IROOT) = EIG(ITER,IROOT)+EIGSHF 1080 1601 CONTINUE 1081* 1082 IF ( IPRT .GT. 1 ) THEN 1083 DO 1600 IROOT = 1, NROOT 1084 WRITE(6,*) 1085 WRITE(6,'(A,I3)') 1086 & ' Information about convergence for root... ' ,IROOT 1087 WRITE(6,*) 1088 & '============================================' 1089 WRITE(6,*) 1090 WRITE(6,1190) FINEIG(IROOT) 1091 1190 FORMAT(' The final approximation to eigenvalue ',F18.10) 1092 IF(IPRT.GE.400) THEN 1093 WRITE(6,1200) 1094 1200 FORMAT(1H0,'The final approximation to eigenvector') 1095 CALL REWINO( LU1) 1096 CALL FRMDSC(VEC1,NVAR,-1 ,LU1,IMZERO,IAMPACK) 1097 CALL WRTMAT(VEC1,1,NVAR,1,NVAR) 1098 END IF 1099 WRITE(6,1300) 1100 1300 FORMAT(1H0,' Summary of iterations ',/,1H 1101 + ,' ----------------------') 1102 WRITE(6,1310) 1103 1310 FORMAT 1104 & (1H0,' Iteration point Eigenvalue Residual ') 1105 DO 1330 I=1,ITER 1106 1330 WRITE(6,1340) I,EIG(I,IROOT)+EIGSHF,RNRM(I,IROOT) 1107 1340 FORMAT(1H ,6X,I4,8X,F20.13,2X,E12.5) 1108 1600 CONTINUE 1109 END IF 1110* 1111 IF(IPRT .EQ. 1 ) THEN 1112 DO 1607 IROOT = 1, NROOT 1113 WRITE(6,'(A,2I3,E13.6,2E10.3)') 1114 & ' >>> CI-OPT Iter Root E g-norm g-red', 1115 & ITER,IROOT,FINEIG(IROOT), 1116 & RNRM(ITER,IROOT), 1117 & RNRM(1,IROOT)/RNRM(ITER,IROOT) 1118 1607 CONTINUE 1119 END IF 1120 1234 CONTINUE 1121C 1122* 1123*. Return final residual norm of roots: 1124 DO IROOT = 1, NROOT 1125 RNRM(IROOT,1) = RNRM(ITER,IROOT) 1126 END DO 1127* 1128 CALL RELUNIT(LUINTM1,'DELETE') 1129 CALL RELUNIT(LUINTM2,'DELETE') 1130 CALL RELUNIT(LUINTM3,'DELETE') 1131 1132 CALL ATIM(CPU,WALL) 1133 CALL PRTIM(6,'time in GEVP solver',CPU-CPU0,WALL-WALL0) 1134* 1135 RETURN 1136 1030 FORMAT(1H0,2X,7F15.8,/,(1H ,2X,7F15.8)) 1137 1120 FORMAT(1H0,2X,I3,7F15.8,/,(1H ,5X,7F15.8)) 1138 END 1139 SUBROUTINE DUMSORT(VEC,NDIM,NELMNT,IREO,ISCR) 1140* 1141* Extremely stupid routine for finding the ordering of 1142* elements in VEC. Only the lowest NELMNT elements are obtained 1143* 1144* On output : IREO: ordered => unordered index 1145* 1146* Author : Prefers to be anonymous, May 2000 1147* 1148* 1149 INCLUDE 'implicit.inc' 1150*. Input 1151 DIMENSION VEC(NDIM) 1152*. Output 1153 INTEGER IREO(NELMNT) 1154*. Scratch through parameter list 1155 INTEGER ISCR(NDIM) 1156* 1157 IZERO = 0 1158 CALL ISETVC(ISCR,IZERO,NDIM) 1159 XMAX = FNDMNX(VEC,NDIM,2) 1160C? WRITE(6,*) ' XMAX = ', XMAX 1161* 1162 DO I = 1, NELMNT 1163*. Find the I'th lowest element 1164 IELMNT = 0 1165 XVAL = XMAX 1166 DO J = 1, NDIM 1167C? WRITE(6,*) ' I,J,XVAL,VEC(J) = ', I,J,XVAL,VEC(J) 1168 IF(VEC(J).LE.XVAL.AND.ISCR(J).EQ.0) THEN 1169C? WRITE(6,*) ' Update of lowest element I and J = ', I,J 1170 XVAL = VEC(J) 1171 IELMNT = J 1172 END IF 1173 END DO 1174 IREO(I) = IELMNT 1175 ISCR(IELMNT) = I 1176 END DO 1177* 1178 NTEST = 10 1179 IF(NTEST.GE.100) THEN 1180 WRITE(6,*) ' Elements to be sorted : ' 1181 CALL WRTMAT(VEC,1,NDIM,1,NDIM) 1182 END IF 1183 IF(NTEST.GE.10) THEN 1184 WRITE(6,*) ' Reorder array : new => old order ' 1185 CALL IWRTMA(IREO,1,NELMNT,1,NELMNT) 1186 END IF 1187* 1188 RETURN 1189 END 1190 SUBROUTINE REO_VEC(IREO,VECIN,NDIM,VECOUT,IWAY) 1191* 1192* Reorder vector 1193* 1194* IWAY = 1 : VECOUT(I) = VECIN(IREO(I)) 1195* IWAY = 2 : VECOUT(IREO(I)) = VECIN(I) 1196* 1197* Jeppe Olsen, May 2000 1198* 1199 INCLUDE 'implicit.inc' 1200*. Input 1201 DIMENSION VECIN(NDIM) 1202 INTEGER IREO(NDIM) 1203*. Output 1204 DIMENSION VECOUT(NDIM) 1205* 1206 IF(IWAY.EQ.1) THEN 1207 DO I = 1, NDIM 1208 VECOUT(I) = VECIN(IREO(I)) 1209 END DO 1210 ELSE 1211 DO I = 1, NDIM 1212 VECOUT(IREO(I)) = VECIN(I) 1213 END DO 1214 END IF 1215* 1216 RETURN 1217 END 1218 SUBROUTINE REO_COL_MAT(IREO,AIN,AOUT,NR,NC,IWAY) 1219* 1220* Reorder columns of matrix 1221* 1222* IWAY = 1 : AOUT(I,J) = AIN(I,IREO(J)) 1223* IWAY = 2 : AOUT(I,IREO(J)) = AIN(I,J) 1224* 1225* Jeppe Olsen, May of 2000 1226* 1227 INCLUDE 'implicit.inc' 1228*. Input 1229 DIMENSION AIN(NR,NC) 1230 INTEGER IREO(*) 1231*. Output 1232 DIMENSION AOUT(NR,NC) 1233* 1234 IF(IWAY.EQ.1) THEN 1235 DO J = 1, NC 1236 CALL COPVEC(AIN(1,IREO(J)),AOUT(1,J),NR) 1237 END DO 1238 ELSE 1239 DO J = 1, NC 1240 CALL COPVEC(AIN(1,J),AOUT(1,IREO(J)),NR) 1241 END DO 1242 END IF 1243* 1244 RETURN 1245 END 1246 SUBROUTINE EIGGMT3(AMAT,NDIM,ARVAL,AIVAL,ARVEC,AIVEC, 1247 & Z,W,SCR,IORD,IEIGVC) 1248* 1249* IF IEIGVC = 1 : Calculate eigenvalues and eigenvalues 1250* = 0 : Calculate only eigenvalues 1251* 1252* Outer routine for calculating eigenvectors and eigenvalues 1253* of a general real matrix 1254* 1255* Version employing EISPACK path RG 1256* 1257* Current implementation is rather wastefull with respect to 1258* memory but at allows one to work with real arithmetic 1259* outside this routine 1260* 1261* If IORD.EQ.1, the eigenvalues are oredered according to the 1262* size of the real part of the eigenvalues 1263* 1264 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1265 REAL * 8 INPROD 1266 DIMENSION AMAT(NDIM,NDIM),SCR(2*NDIM) 1267 DIMENSION ARVAL(NDIM),AIVAL(NDIM) 1268 DIMENSION ARVEC(NDIM,NDIM),AIVEC(NDIM,NDIM) 1269 DIMENSION Z(NDIM,NDIM),W(NDIM) 1270* 1271* Diagonalize 1272* 1273 NSCR = 2*NDIM 1274 IF(IEIGVC.EQ.1) THEN 1275*. Eigenvalues and eigenvectors 1276 CALL RG(NDIM,NDIM,AMAT,ARVAL,AIVAL,1,Z,SCR(1),SCR(1+NDIM),IERR) 1277 ELSE 1278* only eigenvalues : 1279 CALL RG(NDIM,NDIM,AMAT,ARVAL,AIVAL,0,DUMMY,SCR(1),SCR(1+NDIM), 1280 & IERR) 1281 END IF 1282 IF( IERR.NE.0) THEN 1283 WRITE(6,*) ' Problem in EIGGMTN, no convergence ' 1284 WRITE(6,*) ' I have to stop ' 1285 STOP ' No convergence in EIGGMTN ' 1286 END IF 1287* 1288* Extract real and imaginary parts according to Eispack manual p.89 1289* 1290 IF(IEIGVC.EQ.1) THEN 1291 DO 150 K = 1, NDIM 1292* 1293 IF(AIVAL(K).NE.0.0D0) GOTO 110 1294 CALL COPVEC(Z(1,K),ARVEC(1,K),NDIM) 1295 CALL SETVEC(AIVEC(1,K),0.0D0,NDIM) 1296 GOTO 150 1297* 1298 110 CONTINUE 1299 IF(AIVAL(K).LT.0.0D0) GOTO 130 1300 CALL COPVEC(Z(1,K),ARVEC(1,K),NDIM) 1301 CALL COPVEC(Z(1,K+1),AIVEC(1,K),NDIM) 1302 GOTO 150 1303* 1304 130 CONTINUE 1305 CALL COPVEC(ARVEC(1,K-1),ARVEC(1,K),NDIM) 1306 CALL VECSUM(AIVEC(1,K),AIVEC(1,K),AIVEC(1,K-1), 1307 & 0.0D0,-1.0D0,NDIM) 1308* 1309 150 CONTINUE 1310 1311* 1312* explicit orthogonalization of eigenvectors with 1313* (degenerate eigenvalues are not orthogonalized by DGEEV) 1314* 1315 GOTO 201 1316 TEST = 1.0D-11 1317 DO 200 IVEC = 1, NDIM 1318 RNORM = INPROD(ARVEC(1,IVEC),ARVEC(1,IVEC),NDIM) 1319 & + INPROD(AIVEC(1,IVEC),AIVEC(1,IVEC),NDIM) 1320 FACTOR = 1.0d0/SQRT(RNORM) 1321 CALL SCALVE(ARVEC(1,IVEC),FACTOR,NDIM) 1322 CALL SCALVE(AIVEC(1,IVEC),FACTOR,NDIM) 1323 DO 190 JVEC = IVEC+1,NDIM 1324 IF(ARVAL(IVEC)-ARVAL(JVEC).LE.TEST) THEN 1325* orthogonalize jvec to ivec 1326 OVLAPR = INPROD(ARVEC(1,IVEC),ARVEC(1,JVEC),NDIM) 1327 & + INPROD(AIVEC(1,JVEC),AIVEC(1,IVEC),NDIM) 1328 OVLAPI = INPROD(ARVEC(1,IVEC),AIVEC(1,JVEC),NDIM) 1329 & - INPROD(AIVEC(1,IVEC),ARVEC(1,JVEC),NDIM) 1330 CALL VECSUM(ARVEC(1,JVEC),ARVEC(1,JVEC),ARVEC(1,IVEC), 1331 & 1.0D0,-OVLAPR,NDIM ) 1332 CALL VECSUM(AIVEC(1,JVEC),AIVEC(1,JVEC),AIVEC(1,IVEC), 1333 & 1.0D0,-OVLAPR,NDIM ) 1334 CALL VECSUM(ARVEC(1,JVEC),ARVEC(1,JVEC),AIVEC(1,IVEC), 1335 & 1.0D0,OVLAPI,NDIM ) 1336 CALL VECSUM(AIVEC(1,JVEC),AIVEC(1,JVEC),ARVEC(1,IVEC), 1337 & 1.0D0,-OVLAPI,NDIM ) 1338 END IF 1339 190 CONTINUE 1340 200 CONTINUE 1341 201 CONTINUE 1342 1343* 1344* Normalize eigenvectors 1345* 1346 DO 300 L = 1, NDIM 1347 XNORM = INPROD(ARVEC(1,L),ARVEC(1,L),NDIM) 1348 & + INPROD(AIVEC(1,L),AIVEC(1,L),NDIM) 1349 FACTOR = 1.0D0/SQRT(XNORM) 1350 CALL SCALVE(ARVEC(1,L),FACTOR,NDIM) 1351 CALL SCALVE(AIVEC(1,L),FACTOR,NDIM) 1352 300 CONTINUE 1353 END IF 1354* 1355* Order eigensolutions after size of real part of A 1356* 1357 IF(IORD.EQ.1) THEN 1358*. Get reorder array 1359 CALL DUMSORT(ARVAL,NDIM,NDIM,W,SCR) 1360C DUMSORT(VEC,NDIM,NELMNT,IREO,ISCR) 1361*. Reorder eigenvalues 1362C REO_VEC(IREO,VECIN,NDIM,VECOUT,IWAY) 1363 CALL REO_VEC(W,ARVAL,NDIM,SCR,1) 1364 CALL COPVEC(SCR,ARVAL,NDIM) 1365 CALL REO_VEC(W,AIVAL,NDIM,SCR,1) 1366 CALL COPVEC(SCR,AIVAL,NDIM) 1367*. Reorder eigenvectors 1368C REO_COL_MAT(IREO,AIN,AOUT,NR,NC,IWAY) 1369 IF(IEIGVC.EQ.1) THEN 1370 CALL REO_COL_MAT(W,ARVEC,Z,NDIM,NDIM,1) 1371 CALL COPVEC(Z,ARVEC,NDIM**2) 1372 CALL REO_COL_MAT(W,AIVEC,Z,NDIM,NDIM,1) 1373 CALL COPVEC(Z,AIVEC,NDIM**2) 1374 END IF 1375* 1376 END IF 1377* 1378 NTEST = 0 1379 IF(NTEST .GE. 1 ) THEN 1380 WRITE(6,*) ' Output from EIGGMT ' 1381 WRITE(6,*) ' ================== ' 1382 WRITE(6,*) ' Real and imaginary parts of eigenvalues ' 1383 CALL WRTMAT_EP(ARVAL,1,NDIM,1,NDIM) 1384 CALL WRTMAT_EP(AIVAL,1,NDIM,1,NDIM) 1385 END IF 1386* 1387 IF(NTEST.GE.10.AND.IEIGVC.EQ.1) THEN 1388 WRITE(6,*) ' real part of eigenvectors ' 1389 CALL WRTMAT(ARVEC,NDIM,NDIM,NDIM,NDIM) 1390 WRITE(6,*) ' imaginary part of eigenvectors ' 1391 CALL WRTMAT(AIVEC,NDIM,NDIM,NDIM,NDIM) 1392 END IF 1393* 1394 RETURN 1395 END 1396 SUBROUTINE COPMAT(AIN,AOUT,NIN,NOUT) 1397C 1398C COPY MATRIX AIN OF DIMENSION NIN,NIN INTO 1399C MATRIX AOUT OF DIMENSAION NOUT,NOUT 1400C 1401 IMPLICIT DOUBLE PRECISION (A-H,O-Z) 1402 DIMENSION AIN(NIN,NIN) 1403 DIMENSION AOUT(NOUT,NOUT) 1404C 1405 DO 100 J = 1, NOUT 1406 CALL COPVEC(AIN(1,J),AOUT(1,J),NOUT) 1407 100 CONTINUE 1408C 1409 RETURN 1410 END 1411 1412c $Id$ 1413