1* 2* Note on the transformation 3* For a given type of external state and a given symmetry of 4* internal states, the transformation from the elementary EI 5* operators to the orthonormal basis goes as 6* 7* 1: a orthogonal transformation X1 that diagonalizes metric is 8* set up: X1(T) S X1 = Sigma, X1(T) X(1) = 1, dimension: LINT*LINT 9* 2: The nonsingular part of X1 is scaled so it gives eigenvectors 10* of unit norm: X1s = X1 * sigma(-1/2), X1s(T) S X1S = 1, dimension: 11* LINT*LORTN 12* 3: The zeroorder Hamiltonian is set up in the X1s basis and 13* diagonalized by a orthogonal matrix: 14* H0(tilde) = X1s(T) H0 X1s, X2(T) H0(tilde) X2 = epsilon 15* 16* The transformation form the elementary basis to the orthonormal 17* basis reads therefore X = X1 sigma(-1/2) X2, and its inverse is 18* X(-1) = X2(T) sigma(1/2) X1(T). To obtain inverse X(-1) readily 19* the matrices X1, X2 and the vector sigma are stored. 20* 21* Missing: Calculate diagonal 22* <0!0(INT_J)O(EXT(J)! H0!O(EXT_I)O(INT_I)|0> 23*. Assuming that H0 = H0(INT) + H0(EXT) gives 24* <0!0(INT_J)O(EXT(J)! H0!O(EXT_I)O(INT_I)|0> = 25* 26 SUBROUTINE PROJ_TO_NONRED(VECIN,VECOUT,ITSYM,VECSCR) 27* 28*. Project vector to nonredundant space 29* 30*. VECOUT = X1 Sigma^-1 X1^T S VECIN 31* 32*. Jeppe Olsen, Oct 17, 2009 33* 34 INCLUDE 'wrkspc.inc' 35 INCLUDE 'cei.inc' 36 INCLUDE 'ctcc.inc' 37 INCLUDE 'lucinp.inc' 38* 39 REAL*8 40 &INPROD 41*. Input 42 DIMENSION VECIN(*) 43*. Output 44 DIMENSION VECOUT(*) 45*. Scratch 46 DIMENSION VECSCR(*) 47* 48 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'PROJNR') 49 CALL LUCIAQENTER('PRJNR') 50* 51 NTEST = 100 52 53 CALL PROJ_TO_NONRED_SLAVE(VECIN,VECOUT,ITSYM,VECSCR, 54 &WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_SG_INT_EI_FOR_SE), 55 &WORK(KL_S_INT_EI_FOR_SE), 56 &WORK(KL_NDIM_EX_ST),WORK(KL_NDIM_IN_SE), 57 &WORK(KL_N_ORTN_FOR_SE), 58 &WORK(KL_IREO_EI_ST), 59 &N_EXTOP_TP,NSMOB,N_ZERO_EI,NDIM_EI,I_INCLUDE_UNI) 60* 61 XNRM_IN = SQRT(INPROD(VECIN,VECIN,NDIM_EI-1)) 62 XNRM_OUT = SQRT(INPROD(VECOUT,VECOUT,NDIM_EI-1)) 63 XNRM_DIFF = XNORM_DIFF(VECIN,VECOUT,NDIM_EI-1) 64* 65* 66 NTEST = 0 67 IF(NTEST.GE.100) THEN 68 WRITE(6,*) ' Input and output vectors from PROJ_TO_NONRED' 69 CALL WRT_2VEC(VECIN,VECOUT,NDIM_EI) 70 END IF 71 IF(NTEST.GE.1) THEN 72 WRITE(6,*) ' PROJ_TO_NONRED: XNRM_IN, XNRM_OUT,XNRM_DIFF =', 73 & XNRM_IN, XNRM_OUT,XNRM_DIFF 74 END IF 75* 76 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'PROJNR') 77 CALL LUCIAQEXIT('PRJNR') 78 79* 80 RETURN 81 END 82 FUNCTION XNORM_DIFF(VEC1,VEC2,NDIM) 83* 84* Norm of difference of two vectors 85* 86*. Jeppe Olsen, Oct. 18, 2009 87* 88 INCLUDE 'implicit.inc' 89 DIMENSION VEC1(NDIM),VEC2(*) 90* 91 XNORM = 0.0D0 92 DO I = 1, NDIM 93 XNORM = XNORM + (VEC1(I)-VEC2(I))**2 94 END DO 95 XNORM = SQRT(XNORM) 96* 97 XNORM_DIFF = XNORM 98* 99 RETURN 100 END 101* 102 SUBROUTINE PROJ_TO_NONRED_SLAVE(VECIN,VECOUT,ITSYM,VECSCR, 103 &X1,SG,S, 104 &NDIM_EX_ST,NDIM_IN_ST,NDIM_ORT_ST, 105 &IREO_EI_ST, 106 &N_EXTP,NSMOB,N_ZERO_EI,NDIM_EI,I_INCLUDE_UNI) 107* 108*. Project a vector to nonredundant basis 109* 110* Vecout = X1 Sigma^-1 X1^T S Vecin 111* 112* Input and output vectors in CAAB order 113* 114*. Jeppe Olsen, Oct18, 2009 115* 116 INCLUDE 'implicit.inc' 117 INCLUDE 'multd2h.inc' 118*.General input 119* 120 DIMENSION X1(*),SG(*),S(*) 121 INTEGER IREO_EI_ST(*) 122* 123 DIMENSION NDIM_EX_ST(NSMOB,N_EXTP), 124 & NDIM_ORT_ST(NSMOB,N_EXTP), 125 & NDIM_IN_ST(NSMOB,N_EXTP) 126*. Input 127 DIMENSION VECIN(*) 128*. Output 129 DIMENSION VECOUT(*) 130*. Scratch 131 DIMENSION VECSCR(*) 132* 133 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'PROJNS') 134*.1: Reorder to EI order: CAAB(in VECIN) => EI-order(in VECSCR) 135 CALL SGATVEC(VECSCR,VECIN,IREO_EI_ST,NDIM_EI) 136 IF(I_INCLUDE_UNI.EQ.1) THEN 137 VECSCR(NDIM_EI+1) = VECIN(NDIM_EI+1) 138 END IF 139*.2 Perform transformation, save in VECSCR 140 IOFF_SG=1 141 IOFF_S=1 142 IOFF_X1=1 143 IOFF_ORT=1 144 IOFF_VEC = 1 145* 146 DO I_EXTP = 1, N_EXTP 147 DO I_EXSM = 1, NSMOB 148 I_INSM = MULTD2H(I_EXSM,ITSYM) 149* 150 N_EX = NDIM_EX_ST(I_EXSM,I_EXTP) 151 N_IN = NDIM_IN_ST(I_INSM,I_EXTP) 152 N_ORT = NDIM_ORT_ST(I_INSM,I_EXTP) 153* 154 CALL PROJ_EI_BL(X1(IOFF_X1),S(IOFF_S),SG(IOFF_SG), 155 & VECSCR(IOFF_VEC),VECOUT(IOFF_VEC),N_EX,N_IN,N_ORT) 156 CALL COPVEC(VECOUT(IOFF_VEC),VECSCR(IOFF_VEC),N_EX*N_IN) 157* 158 IOFF_SG = IOFF_SG + N_ORT 159 IOFF_S = IOFF_S + N_IN*(N_IN+1)/2 160 IOFF_X1 = IOFF_X1 + N_ORT*N_IN 161 IOFF_VEC = IOFF_VEC + N_EX*N_IN 162* 163 END DO 164 END DO 165*. Reorder to CAAB order 166 CALL SSCAVEC(VECOUT,VECSCR,IREO_EI_ST,NDIM_EI) 167 IF(I_INCLUDE_UNI.EQ.1) THEN 168 VECOUT(NDIM_EI+1) = VECSCR(NDIM_EI+1) 169 END IF 170* 171 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'PROJNS') 172 RETURN 173 END 174 SUBROUTINE PROJ_EI_BL(X1,S,SG,VECIN,VECOUT,N_EX,N_IN,N_ORT) 175* 176* Project a block of a vector to non-redundant basis 177* Input vector is destroyed in the process 178* 179* Vecout = X1 Sigma^-1 X1^T S Vecin 180* 181*. Jeppe Olsen, oct. 18, 2009 182* 183 INCLUDE 'implicit.inc' 184*. General input 185 DIMENSION X1(N_ORT,N_IN),S(N_IN*(N_IN+1)/2),SG(N_ORT) 186*. input and scratch 187C DIMENSION VECIN(N_IN,N_EX) 188 DIMENSION VECIN(N_IN*N_EX) 189*. output 190 DIMENSION VECOUT(N_IN,N_EX) 191* 192 NTEST = 00 193 IF(NTEST.GE.100) THEN 194 WRITE(6,*) ' Input vector to PROJ_EI_BL' 195 CALL WRTMAT(VECIN,N_IN,N_EX,N_IN,N_EX) 196 WRITE(6,*) 'N_IN, N_ORT, N_EX = ', N_IN, N_ORT, N_EX 197 END IF 198*1: S Vecin in Vecout 199 DO I_EX = 1, N_EX 200C CALL SYMTVC(S,VECIN(1,I_EX),VECOUT(1,I_EX),N_IN) 201 CALL SYMTVC(S,VECIN((I_EX-1)*N_IN+1),VECOUT(1,I_EX),N_IN) 202C SYMTVC(A,VECIN,VECOUT,NDIM) 203 END DO 204 205*2: X1^T S Vecin in Vecin 206 FACTOR_AB = 1.0D0 207 FACTOR_C = 0.0D0 208* 209 CALL MATML7(VECIN,X1,VECOUT,N_ORT,N_EX,N_IN,N_ORT,N_IN,N_EX, 210 & FACTOR_C,FACTOR_AB,1) 211*3: Sigma^-1 X1^T S Vecin in Vecin 212 DO J = 1, N_EX 213 DO I = 1, N_ORT 214C VECIN(I,J) = VECIN(I,J)/SG(I) 215 VECIN((J-1)*N_ORT+I) = VECIN((J-1)*N_ORT+I)/SG(I) 216 END DO 217 END DO 218*4: X1 Sigma^-1 X1^T S Vecin in Vecout 219 CALL MATML7(VECOUT,X1,VECIN,N_IN,N_EX,N_IN,N_ORT,N_ORT,N_EX, 220 & FACTOR_C,FACTOR_AB,0) 221* 222 IF(NTEST.GE.100) THEN 223 WRITE(6,*) ' Output block from PROJ_EI_BL' 224 CALL WRTMAT(VECOUT,N_IN,N_EX,N_IN,N_EX) 225 END IF 226* 227 RETURN 228 END 229 SUBROUTINE SYMMAT_TV(SYMMAT,VECIN,VECOUT,NDIM) 230* 231* A symmetric matrix is stored rowwise, lowerhalf 232* multiply with vector VECIN 233* 234*. NO cache-considerations have been invoked 235* 236*. Jeppe Olsen, oct. 18, 2009 237* 238 INCLUDE 'implicit.inc' 239*.input 240 DIMENSION SYMMAT(*) 241 DIMENSION VECIN(*) 242*.Output 243 DIMENSION VECOUT(*) 244* 245 ZERO = 0.0D0 246 CALL SETVEC(VECOUT,ZERO,NDIM) 247 IJ = 0 248 DO I = 1, NDIM 249 DO J = 1, I 250 IJ = IJ + 1 251 VECOUT(J) = VECOUT(J) + SYMMAT(IJ)*VECIN(I) 252 VECOUT(I) = VECOUT(I) + SYMMAT(IJ)*VECIN(J) 253 END DO 254 END DO 255* 256 NTEST = 100 257 IF(NTEST.GE.100) THEN 258 WRITE(6,*) ' Output vector from SYMMAT_TV' 259 CALL WRTMAT(VECOUT,1,NDIM,1,NDIM) 260 END IF 261* 262 RETURN 263 END 264 SUBROUTINE MODDIAG(H0DIAG,NDIM,XMIN) 265* 266*. Jacobian Diagonal of H0 is given. 267*. Check this and replace all values smaller than XMIN 268* by XMIN 269* 270*. Jeppe Olsen 271* 272 INCLUDE 'implicit.inc' 273 DIMENSION H0DIAG(NDIM) 274* 275 XMIN_FOUND = ABS(H0DIAG(1)) 276 NSHIFT = 0 277* 278 DO I = 1, NDIM 279 IF(ABS(H0DIAG(I)).LT.XMIN_FOUND) THEN 280 XMIN_FOUND = H0DIAG(I) 281 END IF 282 IF(ABS(H0DIAG(I)).LT.XMIN) THEN 283 H0DIAG(I) = XMIN 284 NSHIFT = NSHIFT + 1 285 END IF 286 END DO 287* 288 NTEST = 10 289 IF(NTEST.GE.1) THEN 290 WRITE(6,*) ' Check of J(H0DIAG)' 291 WRITE(6,*) ' Imposed lower value =', XMIN 292 WRITE(6,*) ' Number of elements shifted =', NSHIFT 293 WRITE(6,*) ' Lowest value found =', XMIN_FOUND 294 END IF 295* 296 RETURN 297 END 298 SUBROUTINE TEST_GENMAT(MSTV,NVAR,NVAR_INT,I_DO_DIAG) 299* 300* A method using a matrix vector routine MSTV is 301* tested by constructing complete matrix and metric and maybe 302* diagonalizing (If I_DO_DIAG = 1) 303* 304*. NVAT_INT: Sometimes a different number of variables 305*. are used internally in matrix vector. The vectors 306*. should have this dimension 307* 308* Jeppe Olsen, Hotel room in Dusseldorf, 2009 309* 310 INCLUDE 'wrkspc.inc' 311 EXTERNAL MSTV 312* 313 IDUM = 0 314 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'TST_GN') 315* 316 WRITE(6,*) ' NVAR, NVAR_INT =',NVAR,NVAR_INT 317 CALL MEMMAN(KLMATH,NVAR**2 ,'ADDL ',2,'MATH ') 318 CALL MEMMAN(KLMATS,NVAR**2 ,'ADDL ',2,'MATS ') 319 CALL MEMMAN(KLVEC1,NVAR_INT,'ADDL ',2,'VEC1 ') 320 CALL MEMMAN(KLVEC2,NVAR_INT,'ADDL ',2,'VEC2 ') 321 CALL MEMMAN(KLVEC3,NVAR_INT,'ADDL ',2,'VEC3 ') 322* 323*. Space for diagonalization(I use a simple routine for this) 324 IF(I_DO_DIAG.EQ.1) THEN 325 CALL MEMMAN(KLSCRM1,NVAR**2,'ADDL ',2,'MAT1 ') 326 CALL MEMMAN(KLSCRM2,NVAR**2,'ADDL ',2,'MAT2 ') 327 ELSE 328 KLSCRM1 = 1 329 KLSCRM2 = 1 330 END IF 331* 332 CALL TEST_GENMAT_INNER(MSTV,NVAR,I_DO_DIAG, 333 & WORK(KLMATH),WORK(KLMATS), 334 & WORK(KLVEC1),WORK(KLVEC2),WORK(KLVEC3), 335 & WORK(KLSCRM1),WORK(KLSCRM2)) 336* 337 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'TST_GN') 338* 339 RETURN 340 END 341 SUBROUTINE TEST_GENMAT_INNER(MSTV,NVAR,I_DO_DIAG, 342 & HMAT,SMAT,VEC1,VEC2,VEC3,SCRM1,SCRM2) 343* 344* A method using a matrix vector routine MSTV is 345* tested by constructing complete matrix and metric and maybe 346* diagonalizing - Inner routine ( Well, sounds better than slave 347* routine) 348* 349* Jeppe Olsen, Hotel room in Dusseldorf, 2009 350* 351 INCLUDE 'wrkspc.inc' 352 EXTERNAL MSTV 353 DIMENSION HMAT(NVAR,NVAR),SMAT(NVAR,NVAR) 354 DIMENSION VEC1(NVAR),VEC2(NVAR),VEC3(NVAR) 355 DIMENSION SCRM1(*),SCRM2(*) 356* 357 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'TS_GNI') 358* 359 NTEST = 1000 360 IF(NTEST.GE.100) THEN 361 WRITE(6,*) ' TEST_GENMAT_INNER reporting ' 362 END IF 363*. Test: just parameter 136 364C? ZERO = 0.0D0 365C? CALL SETVEC(VEC1,ZERO,NVAR) 366C? VEC1(136) = 1.0D0 367C? CALL MSTV(VEC1,VEC2,VEC3,1,1) 368C? STOP 'Enforced stop after call to test MSTV' 369* 370 DO I = 1, NVAR 371 ZERO = 0.0D0 372 CALL SETVEC(VEC1,ZERO,NVAR) 373 ONE = 1.0D0 374 VEC1(I) = ONE 375 WRITE(6,*) ' Constructing HV, SV for element = ', I 376 WRITE(6,*) ' --------------------------------------' 377 CALL MSTV(VEC1,VEC2,VEC3,1,1) 378 CALL COPVEC(VEC2,HMAT(1,I),NVAR) 379 CALL COPVEC(VEC3,SMAT(1,I),NVAR) 380 END DO 381*. Compare to unit matrix 382 CALL COMPARE_TO_UNI(SMAT,NVAR) 383C SUBROUTINE COMPARE_TO_UNI(A,NDIM) 384 IF(NTEST.GE.100) THEN 385 WRITE(6,*) ' H-matrix: ' 386 CALL WRTMAT(HMAT,NVAR,NVAR,NVAR,NVAR) 387 WRITE(6,*) ' S-matrix: ' 388 CALL WRTMAT(SMAT,NVAR,NVAR,NVAR,NVAR) 389 END IF 390* 391 IF (I_DO_DIAG.EQ.1) THEN 392*. Diagonalize 393C GENDIA(HIN,SIN,VOUT,EIGENV,PVEC,NDIM,ISORT) 394 CALL GENDIA(HMAT,SMAT,SCRM1,VEC1,SCRM2,NVAR,1) 395 END IF 396* 397 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'TS_GNI') 398 RETURN 399 END 400* 401 SUBROUTINE PRINT_ZERO_TRMAT 402* 403* Print transformation matrices for obtaining orthogonal 404* internal zero-order states 405* 406*. Jeppe Olsen, the train to dusseldorf, aug14, 2009 407* 408 INCLUDE 'wrkspc.inc' 409 INCLUDE 'lucinp.inc' 410 INCLUDE 'cei.inc' 411* 412 WRITE(6,*) ' ----------------------------------------------' 413 WRITE(6,*) ' Transformation matrices for orthonormal states' 414 WRITE(6,*) ' ----------------------------------------------' 415 WRITE(6,*) ' PRINT_ZERO.(a): NSMOB, N_EXTOP_TP=',NSMOB,N_EXTOP_TP 416 417 CALL PRINT_ZERO_TRMAT_SLAVE(NSMOB,N_EXTOP_TP, 418 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 419 & WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_X2_INT_EI_FOR_SE), 420 & WORK(KL_SG_INT_EI_FOR_SE), 421 & WORK(KL_IBX1_INT_EI_FOR_SE),WORK(KL_IBX2_INT_EI_FOR_SE), 422 & WORK(KL_IBSG_INT_EI_FOR_SE)) 423* 424 RETURN 425 END 426 SUBROUTINE PRINT_ZERO_TRMAT_SLAVE(NSMOB,N_EXTP, 427 & N_ORTN_FOR_SE,N_INT_FOR_SE, 428 & X1,X2,SG, 429 & IBX1_INT_EI_FOR_SE,IBX2_INT_EI_FOR_SE, 430 & IBSG_INT_EI_FOR_SE) 431*. Print matrices relevant for generation of orthonormal 432*. zero-order stated 433* 434*. Jeppe Olsen, Aug.14, 2009 435* 436 INCLUDE 'implicit.inc' 437 DIMENSION X1(*),X2(*),SG(*) 438* 439 INTEGER N_ORTN_FOR_SE(NSMOB,N_EXTP),N_INT_FOR_SE(NSMOB,N_EXTP) 440* 441 INTEGER IBX1_INT_EI_FOR_SE(NSMOB,N_EXTP) 442 INTEGER IBX2_INT_EI_FOR_SE(NSMOB,N_EXTP) 443 INTEGER IBSG_INT_EI_FOR_SE(NSMOB,N_EXTP) 444* 445 NTEST = 0 446 WRITE(6,*) ' PRINT_ZERO.. : N_EXTP, NSMOB = ', N_EXTP, NSMOB 447 DO I_EXTP = 1, N_EXTP 448 DO I_EXSM = 1, NSMOB 449 WRITE(6,*) ' Info for external type and sym =', I_EXTP,I_EXSM 450 N_INT = N_INT_FOR_SE(I_EXSM,I_EXTP) 451 N_ORTN= N_ORTN_FOR_SE(I_EXSM,I_EXTP) 452* 453 IOFF_X1 = IBX1_INT_EI_FOR_SE(I_EXSM,I_EXTP) 454 IOFF_X2 = IBX2_INT_EI_FOR_SE(I_EXSM,I_EXTP) 455 IOFF_SG = IBSG_INT_EI_FOR_SE(I_EXSM,I_EXTP) 456*. Test output for printing routines - I must be getting old 457 IF(NTEST.GE.100) THEN 458 WRITE(6,*) ' N_INT, N_ORTN = ', N_INT, N_ORTN 459 WRITE(6,*) ' IOFF_X1, IOFF_X2 ,IOFF_SG =', 460 & IOFF_X1, IOFF_X2 ,IOFF_SG 461 END IF 462* 463 WRITE(6,*) ' Block of X1, X2, SG ' 464 CALL WRTMAT(X1(IOFF_X1),N_INT,N_ORTN,N_INT,N_ORTN) 465 WRITE(6,*) 466 CALL WRTMAT(X2(IOFF_X2),N_ORTN,N_ORTN,N_ORTN,N_ORTN) 467 WRITE(6,*) 468 CALL WRTMAT(SG(IOFF_SG),1,N_ORTN,1,N_ORTN) 469 WRITE(6,*) 470 END DO 471 END DO 472* 473 RETURN 474 END 475 SUBROUTINE TRANS_CAAB_ORTN(T_CAAB,T_ORTN,ITSYM,ICO,ILR,SCR, 476 & ICOCON) 477* 478*. Transform a vector between standard CAAB order and EI order 479*. with orthornormal internal states 480* 481* ICO = 1: CAAB => Ortn 482* ICO = 2: Ortn => CAAB 483* 484* ICOCON = 1 => Covariant transformation 485* ICOCON = 2 => Contravariant transformation 486* 487*. Jeppe Olsen, July 29, 2009, sidder p� R�nnevej 12 med Jette 488* g�ttende kryds og tv�rs... 489* 490*. NOTE: at the moment no signs are used. When transforming tp 491* spin-adapted operators, sign changes going from CAAB tp EI form 492* must be considered. 493* 494 INCLUDE 'wrkspc.inc' 495 INCLUDE 'cei.inc' 496*. Input/output 497 DIMENSION T_ORTN(*),T_CAAB(*) 498*. Scratch: Dimension of full CAAB expansion 499 DIMENSION SCR(*) 500* 501 NTEST = 000 502 IF(NTEST.GE.100) WRITE(6,*) ' Starting with transformation' 503* 504* 505 IF(ICO.EQ.1) THEN 506*. CAAB => Ortn is done in two steps 507* CAAB(in T_CAAB) => EI-order(in SCR) => Ortn(in T_ORTN) 508 CALL SGATVEC(SCR,T_CAAB,WORK(KL_IREO_EI_ST),NDIM_EI) 509 IF(I_INCLUDE_UNI.EQ.1) THEN 510 SCR(NDIM_EI+1) = T_CAAB(NDIM_EI+1) 511 END IF 512 CALL TRANS_EI_ORTN(SCR,T_ORTN,ITSYM,1,ILR,ICOCON) 513 ELSE 514*. Ortn => CAAB is done in two steps 515*. Ortn(in T_ORTN) => EI-order(in SCR) => CAAB (in T_CAAB) 516 CALL TRANS_EI_ORTN(SCR,T_ORTN,ITSYM,2,ILR,ICOCON) 517 CALL SSCAVEC(T_CAAB,SCR,WORK(KL_IREO_EI_ST),NDIM_EI) 518* 519 IF(I_INCLUDE_UNI.EQ.1) THEN 520 T_CAAB(NDIM_EI+1) = SCR(NDIM_EI+1) 521 END IF 522 END IF 523* 524 IF(NTEST.GE.100) THEN 525 IF(ICO.EQ.1) THEN 526 WRITE(6,*) ' CAAB => ORTN transformation ' 527 ELSE 528 WRITE(6,*) ' ORTN => CAAB transformation ' 529 END IF 530 IF(ICOCON.EQ.1) THEN 531 WRITE(6,*) ' Covariant transformation' 532 ELSE 533 WRITE(6,*) ' Contravariant transformation' 534 END IF 535* 536 WRITE(6,*) ' The T vector in zero-order EI basis: ' 537 CALL PRINT_T_EI(T_ORTN,2,ITSYM) 538COLD IF(I_INCLUDE_UNI.EQ.1) THEN 539COLD WRITE(6,*) ' Element corresponding to unit operator', 540COLD & T_ORTN(N_ZERO_EI+1) 541COLD END IF 542 WRITE(6,*) ' The T vector in CAAB form' 543 CALL WRTMAT(T_CAAB,1,NDIM_EI,1,NDIM_EI) 544 IF(I_INCLUDE_UNI.EQ.1) THEN 545 WRITE(6,*) ' Element corresponding to unit operator', 546 & T_CAAB(NDIM_EI+1) 547 END IF 548* 549 END IF 550 IF(NTEST.GE.100) WRITE(6,*) ' Finished with transformation' 551* 552 RETURN 553 END 554 SUBROUTINE PRINT_T_EI(T,IEO,ITSYM) 555* 556* Print a T(I,E) vector given in elementary(IEO=1) or 557* orthonormal(IEO=2) form 558* 559*. Jeppe Olsen, July 29, 2009 560* 561 INCLUDE 'wrkspc.inc' 562 INCLUDE 'cei.inc' 563 INCLUDE 'lucinp.inc' 564*. Specific input 565 DIMENSION T(*) 566* 567 WRITE(6,*) 568 WRITE(6,*) ' ------------------------------------' 569 WRITE(6,*) ' T(I,E) vector with symmetry ', ITSYM 570 WRITE(6,*) ' ------------------------------------' 571 WRITE(6,*) 572 IF(IEO.EQ.1) THEN 573 CALL PRINT_T_EI_SLAVE(T,ITSYM,N_EXTOP_TP, 574 & WORK(KL_NDIM_EX_ST),WORK(KL_NDIM_IN_SE),NSMOB) 575 IF(I_INCLUDE_UNI.EQ.1) THEN 576 WRITE(6,*) ' Element corresponding to unit-operator', 577 & T(NDIM_EI+1) 578 END IF 579 ELSE 580 CALL PRINT_T_EI_SLAVE(T,ITSYM,N_EXTOP_TP, 581 & WORK(KL_NDIM_EX_ST),WORK(KL_N_ORTN_FOR_SE),NSMOB) 582 IF(I_INCLUDE_UNI.EQ.1) THEN 583 WRITE(6,*) ' Element corresponding to unit-operator', 584 & T(N_ZERO_EI+1) 585 END IF 586 END IF 587* 588 RETURN 589 END 590 SUBROUTINE PRINT_T_EI_SLAVE(T,ITSYM,N_EXTP, 591 & NDIM_EX_ST,NDIM_IN_ST,NSMOB) 592* 593 INCLUDE 'implicit.inc' 594 INCLUDE 'multd2h.inc' 595*. Specific input 596 DIMENSION T(*) 597 DIMENSION NDIM_EX_ST(NSMOB,N_EXTP),NDIM_IN_ST(NSMOB,N_EXTP) 598* 599 IOFF = 1 600 DO I_EXTP = 1, N_EXTP 601 DO I_EXSM = 1, NSMOB 602 I_INSM = MULTD2H(I_EXSM,ITSYM) 603 WRITE(6,*) ' Block with external type and sym :', I_EXTP, I_EXSM 604 N_EX = NDIM_EX_ST(I_EXSM,I_EXTP) 605 N_IN = NDIM_IN_ST(I_INSM,I_EXTP) 606C? WRITE(6,*) ' N_EX, N_IN = ', N_EX, N_IN 607 CALL WRTMAT(T(IOFF),N_IN,N_EX,N_IN,N_EX) 608 IOFF = IOFF + N_EX*N_IN 609 END DO 610 END DO 611* 612 RETURN 613 END 614 SUBROUTINE GET_DIAG_H0_EI(DIAG) 615* 616* Construct diagonal of H0 over orthonormal states 617* 618*. Jeppe Olsen, March 2009 619* 620 INCLUDE 'wrkspc.inc' 621 INCLUDE 'orbinp.inc' 622 INCLUDE 'cgas.inc' 623 INCLUDE 'cei.inc' 624 INCLUDE 'ctcc.inc' 625 INCLUDE 'glbbas.inc' 626 INCLUDE 'cintfo.inc' 627 INCLUDE 'lucinp.inc' 628*. Output 629 DIMENSION DIAG(*) 630* 631 NTEST = 100 632 IDUM = 0 633 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GT_DIA') 634C? WRITE(6,*) ' Entering GET_DIAG.., WORK(KINT1) = ', WORK(KINT1) 635*. Generate H0 as FI + FA - starting from scatch to be on the safe side.. 636 CALL COPVEC(WORK(KH),WORK(KHINA),NINT1) 637*. Inactive Fock matric 638 CALL FISM(WORK(KHINA),ECC) 639*. Active Fock matrix 640 CALL FAM(WORK(KFIFA)) 641*. and add 642 ONE = 1.0D0 643 CALL VECSUM(WORK(KFIFA),WORK(KFIFA),WORK(KHINA), 644 & ONE,ONE,NINT1) 645* 646 IF(NTEST.GE.100) THEN 647 WRITE(6,*) ' FIFA from GET_DIAG_H0_EI: ' 648 CALL APRBLM2(WORK(KFIFA),NTOOBS,NTOOBS,NSMOB,1) 649 END IF 650*. External part of zero-order energy: External part is assumed to 651*. be doubly occupied hole space and unoccupied particle space 652 E0REF_EXT = 0.0D0 653*. Obtain diagonal of H0 654 CALL MEMMAN(KLH0DIAS,NTOOB,'ADDL ',2,'H0DIAS') 655 CALL MEMMAN(KLH0DIA,NTOOB,'ADDL ',2,'H0DIA ') 656 CALL GET_DIAG_BLOC_MAT(WORK(KFIFA),WORK(KLH0DIAS),NSMOB,NTOOBS,1) 657C GET_DIAG_BLOC_MAT(A,ADIAG,NBLOCK,LBLOCK,ISYM) 658 IF(NTEST.GE.100) THEN 659 WRITE(6,*)' Diagonal of FIFA in sym-order' 660 CALL WRTMAT(WORK(KLH0DIAS),1,NTOOB,1,NTOOB) 661 END IF 662*. Was obtained in symmetry ordered basis, type used below so 663 DO I= 1,NTOOB 664 WORK(KLH0DIA-1+IREOST(I)) = WORK(KLH0DIAS-1+I) 665 END DO 666 667*. 4 Blocks for occupation of external strings 668 CALL MEMMAN(KL_IST_EX_CA,MX_ST_TSOSO_BLK_MX,'ADDL ',1,'STE_CA') 669 CALL MEMMAN(KL_IST_EX_CB,MX_ST_TSOSO_BLK_MX,'ADDL ',1,'STE_CB') 670 CALL MEMMAN(KL_IST_EX_AA,MX_ST_TSOSO_BLK_MX,'ADDL ',1,'STE_AA') 671 CALL MEMMAN(KL_IST_EX_AB,MX_ST_TSOSO_BLK_MX,'ADDL ',1,'STE_AB') 672* 673C GET_DIAG_H0_EI_SLAVE(DIAG, 674C & ISPOBEX_TP,N_EXTP,E0REF_EXT, 675C & N_ORTN_FOR_SE, 676C & I_IN_TP,IB_INTP,IB_EXTP,H0DIAG, 677C & IST_EX_CA, IST_EX_CB, IST_EX_AA,IST_EX_AB) 678C? WRITE(6,*) ' I_INT_OFF before callto GET_DI..SLAVE',I_INT_OFF 679 CALL GET_DIAG_H0_EI_SLAVE(DIAG,WORK(KLSOBEX),N_EXTOP_TP,E0REF_EXT, 680 & WORK(KL_N_ORTN_FOR_SE),I_IN_TP, 681 & I_INT_OFF,I_EXT_OFF,WORK(KLH0DIA), 682 & WORK(KL_IST_EX_CA),WORK(KL_IST_EX_CB), 683 & WORK(KL_IST_EX_AA),WORK(KL_IST_EX_AB)) 684*. Reinstall one-electron integrals 685 CALL COPVEC(WORK(KINT1O),WORK(KFIFA),NINT1) 686C? WRITE(6,*) ' Leaving GET_DIAG.., WORK(KINT1) = ', WORK(KINT1) 687* 688 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'GT_DIA') 689* 690 RETURN 691 END 692 SUBROUTINE GET_DIAG_H0_EI_SLAVE(DIAG, 693 & ISPOBEX_TP,N_EXTP,E0REF_EXT, 694 & N_ORTN_FOR_SE, 695 & I_IN_TP,IB_INTP,IB_EXTP,H0DIAG, 696 & IST_EX_CA, IST_EX_CB, IST_EX_AA,IST_EX_AB) 697* 698* Generate diagonal of H0 over zero-order states using EI approach 1 699* 700* Jeppe Olsen, March 2009, trying to do a bit of science 701* 702 INCLUDE 'wrkspc.inc' 703 REAL*8 INPROD 704* 705 INCLUDE 'cgas.inc' 706 INCLUDE 'gasstr.inc' 707 INCLUDE 'lucinp.inc' 708 INCLUDE 'ctcc.inc' 709 INCLUDE 'cstate.inc' 710 INCLUDE 'cands.inc' 711 INCLUDE 'multd2h.inc' 712 INCLUDE 'glbbas.inc' 713 INCLUDE 'clunit.inc' 714 INCLUDE 'oper.inc' 715 INCLUDE 'cintfo.inc' 716*. E0REF_EXT is external parts of E0 of reference state 717*. Input: 718 INTEGER ISPOBEX_TP(NGAS,4,*) 719*. Number of orthonormal internal states per symmetry and external type 720 INTEGER N_ORTN_FOR_SE(NSMOB,N_EXTP) 721*. Diagonal part of H0 in orbital basis 722 DIMENSION H0DIAG(*) 723*. Output 724 DIMENSION DIAG(*) 725*. Scratch 726*. Blocks for holding group of external strings with given sym 727 DIMENSION IST_EX_CA(*),IST_EX_CB(*), IST_EX_AA(*),IST_EX_AB(*) 728*. Local scratch 729 INTEGER IGRP(MXPNGAS) 730* 731 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GT_DIA') 732* 733 NTEST = 10 734 IF(NTEST.GE.10) THEN 735 WRITE(6,*) ' -------------- ' 736 WRITE(6,*) ' GET_DIAG_H0_EI ' 737 WRITE(6,*) ' -------------- ' 738 WRITE(6,*) 739C? WRITE(6,*) ' E0REF_EXT at start of GET_DIAG.. ', E0REF_EXT 740 END IF 741*. Symmetry of O+(e)O+(i) 742 ISYM_EI = 1 743*. Largest number of internal zero-order states for given symmetry 744C IMNXVC(IVEC,NDIM,MXMN,IVAL,IPLACE) 745 CALL IMNXVC(N_ORTN_FOR_SE,NSMOB*N_EXTP,1,NMAX_ORTN_FOR_SE,IPLACE) 746 IF(NTEST.GE.1000) THEN 747 WRITE(6,*) ' NMAX_ORTN_FOR_SE = ', NMAX_ORTN_FOR_SE 748 CALL IWRTMA(N_ORTN_FOR_SE,NSMOB,N_EXTP,NSMOB,N_EXTP) 749 END IF 750*. an array for internal state energies 751 CALL MEMMAN(KLE0INT,NMAX_ORTN_FOR_SE,'ADDL ',2,'E0_INT') 752* 753 IEI_ORTN = 0 754 DO J_EXTP = 1, N_EXTP 755 J_EXTP_ABS = IB_EXTP-1+J_EXTP 756* 757 NEL_EX_CA = IELSUM(ISPOBEX_TP(1,1,J_EXTP_ABS),NGAS) 758 NEL_EX_CB = IELSUM(ISPOBEX_TP(1,2,J_EXTP_ABS),NGAS) 759 NEL_EX_AA = IELSUM(ISPOBEX_TP(1,3,J_EXTP_ABS),NGAS) 760 NEL_EX_AB = IELSUM(ISPOBEX_TP(1,4,J_EXTP_ABS),NGAS) 761* 762*. Loop over strings in EI order 763 DO I_EX_SM = 1, NSMOB 764 I_IN_SM = MULTD2H(I_EX_SM,ISYM_EI) 765 IF(NTEST.GE.20) WRITE(6,*) ' J_EXTP, I_EX_SM,I_IN_SM =', 766 & J_EXTP, I_EX_SM,I_IN_SM 767*. Obtain diagonal of internal energy contributions 768C GET_BLOCK_OF_HS_IN_INTERNAL_SPACE( 769C & IEXTP,IINTSM,I_HS,HSBLCK,I_INT_TP,I_ONLY_DIA) 770 CALL GET_BLOCK_OF_HS_IN_INTERNAL_SPACE( 771 & J_EXTP,I_IN_SM,1,WORK(KLE0INT),I_INT_TP,1) 772 IF(NTEST.GE.20) THEN 773 L_INTOP = N_ORTN_FOR_SE(I_IN_SM,J_EXTP) 774 WRITE(6,*) ' H0 energies of internal states' 775 CALL WRTMAT(WORK(KLE0INT),1,L_INTOP,1,L_INTOP) 776 END IF 777*. Loop over symmetries of external operators 778 DO I_EX_C_SM = 1, NSMOB 779 I_EX_A_SM = MULTD2H(I_EX_C_SM,I_EX_SM) 780 DO I_EX_CA_SM = 1, NSMOB 781 I_EX_CB_SM = MULTD2H(I_EX_CA_SM,I_EX_C_SM) 782 DO I_EX_AA_SM = 1, NSMOB 783 I_EX_AB_SM = MULTD2H(I_EX_AA_SM,I_EX_A_SM) 784* 785 IF(NTEST.GE.1000) THEN 786 WRITE(6,'(A,4I4)') 787 & 'I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM = ', 788 & I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM 789 END IF 790*. Occupation of external strings 791*. CA in IST_EX_CA 792 CALL OCC_TO_GRP(ISPOBEX_TP(1,1,J_EXTP_ABS),IGRP,1) 793 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 794 & I_EX_CA_SM,NEL_EX_CA,NSTR_EX_CA,IST_EX_CA,NTOOB,0, 795 & IDUM,IDUM) 796*. CB in IST_EX_CB 797 CALL OCC_TO_GRP(ISPOBEX_TP(1,2,J_EXTP_ABS),IGRP,1) 798 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 799 & I_EX_CB_SM,NEL_EX_CB,NSTR_EX_CB,IST_EX_CB,NTOOB,0, 800 & IDUM,IDUM) 801*. AA in IST_EX_AA 802 CALL OCC_TO_GRP(ISPOBEX_TP(1,3,J_EXTP_ABS),IGRP,1) 803 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 804 & I_EX_AA_SM,NEL_EX_AA,NSTR_EX_AA,IST_EX_AA,NTOOB,0, 805 & IDUM,IDUM) 806*. AB in IST_EX_AB 807 CALL OCC_TO_GRP(ISPOBEX_TP(1,4,J_EXTP_ABS),IGRP,1) 808 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 809 & I_EX_AB_SM,NEL_EX_AB,NSTR_EX_AB,IST_EX_AB,NTOOB,0, 810 & IDUM,IDUM) 811C? WRITE(6,*) ' End of external ' 812*. Loop over External strings (CA, CB, AA, AB) 813 DO I_EX_AB = 1, NSTR_EX_AB 814*. Contribution from I_EX_AB to diagonal 815C E1_FOR_STRING(HDIAG,ISTRING,NEL) 816 D_EX_AB = 817 & E1_FOR_STRING(H0DIAG,IST_EX_AB(1+(NEL_EX_AB)*(I_EX_AB-1)), 818 & NEL_EX_AB) 819 DO I_EX_AA = 1, NSTR_EX_AA 820 D_EX_AA = 821 & E1_FOR_STRING(H0DIAG,IST_EX_AA(1+(NEL_EX_AA)*(I_EX_AA-1)), 822 & NEL_EX_AA) 823 DO I_EX_CB = 1, NSTR_EX_CB 824 D_EX_CB = 825 & E1_FOR_STRING(H0DIAG,IST_EX_CB(1+(NEL_EX_CB)*(I_EX_CB-1)), 826 & NEL_EX_CB) 827 DO I_EX_CA = 1, NSTR_EX_CA 828 D_EX_CA = 829 & E1_FOR_STRING(H0DIAG,IST_EX_CA(1+(NEL_EX_CA)*(I_EX_CA-1)), 830 & NEL_EX_CA) 831 E_EXT = E0REF_EXT + D_EX_CB +D_EX_CA -D_EX_AB - D_EX_AA 832 IF(NTEST.GE.1000) THEN 833 WRITE(6,*) ' D_EX_CB,D_EX_CA,D_EX_AB,D_EX_AA =', 834 & D_EX_CB,D_EX_CA,D_EX_AB,D_EX_AA 835 END IF 836CM*. Obtain diagonal of internal energy contributions 837C? CALL GET_BLOCK_OF_HS_IN_INTERNAL_SPACE( 838C? & J_EXTP,I_IN_SM,1,WORK(KLE0INT),I_INT_TP,1) 839* 840 L_INTOP = N_ORTN_FOR_SE(I_IN_SM,J_EXTP) 841 DO I_INTOP = 1, L_INTOP 842 IEI_ORTN = IEI_ORTN + 1 843 IF(NTEST.GE.1000) THEN 844 WRITE(6,*) ' I, Internal, external term =', 845 & WORK(KLE0INT-1+I_INTOP), E_EXT 846 END IF 847 DIAG(IEI_ORTN) = WORK(KLE0INT-1+I_INTOP) + E_EXT 848 849 END DO 850*. ^ End of loop over internal states 851 END DO 852 END DO 853 END DO 854 END DO 855* ^ End of loop over external CA,CB,AA,AB strings of given sym 856 END DO 857 END DO 858 END DO 859*. ^ End of loop over symmetry of external CA, CB, AA, AB strings 860 END DO 861*. ^ End of loop over symmetry of external operators 862 END DO 863*. ^ End of loop over external types 864* 865 IF(NTEST.GE.1000) THEN 866 WRITE(6,*) ' Diagonal of zero-order energy of ortn states' 867 CALL WRTMAT(DIAG,1,IEI_ORTN,1,IEI_ORTN) 868 END IF 869* 870 CALL MEMMAN(IDUM,IDUM,'FLUSM',1,'GT_DIA') 871* 872 RETURN 873 END 874 SUBROUTINE LUCIA_IC_EI 875 & (IREFSPCE,ITREFSPC,ICTYP,EREF,I_DO_CUMULANTS, 876 & EFINAL,CONVER,VNFINAL) 877* 878* 879* Master routine for internally contracted calculation 880* 881* Specialized for reference states that allows division into 882* internal and external parts 883* 884* Jeppe Olsen, September 06 885* 886* The spaces are assumed organized as 887* 1 : Reference space on which excitations are performed (IREFSPC) 888* 2 : Space that defines excitations (ITREFSPC) 889* 3 : Space where calculation is performed (ITREFSPC) 890* 891* Deviations from this will cause trouble 892* 893 INCLUDE 'wrkspc.inc' 894 REAL*8 895 &INPROD 896 INCLUDE 'crun.inc' 897 INCLUDE 'cstate.inc' 898 INCLUDE 'cgas.inc' 899 INCLUDE 'ctcc.inc' 900 INCLUDE 'gasstr.inc' 901 INCLUDE 'strinp.inc' 902 INCLUDE 'orbinp.inc' 903 INCLUDE 'cprnt.inc' 904 INCLUDE 'corbex.inc' 905 INCLUDE 'csm.inc' 906 INCLUDE 'cicisp.inc' 907 INCLUDE 'cecore.inc' 908 INCLUDE 'glbbas.inc' 909 INCLUDE 'clunit.inc' 910 INCLUDE 'lucinp.inc' 911 INCLUDE 'cc_exc.inc' 912 INCLUDE 'cintfo.inc' 913 INCLUDE 'cei.inc' 914 INCLUDE 'oper.inc' 915 INCLUDE 'newccp.inc' 916*. Transfer common block for communicating with H_EFF * vector routines 917 COMMON/COM_H_S_EFF_ICCI_TV/ 918 & C_0X,KLTOPX,NREFX,IREFSPCX,ITREFSPCX,NCAABX, 919 & IUNIOPX,NSPAX,IPROJSPCX 920*. A bit of local scratch 921 DIMENSION ICASCR(MXPNGAS) 922 CHARACTER*6 ICTYP 923* 924 EXTERNAL MTV_FUSK, STV_FUSK 925 EXTERNAL H_S_EFF_ICCI_TV,H_S_EXT_ICCI_TV 926 EXTERNAL HOME_SD_INV_T_ICCI 927 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'ICCI ') 928* 929 I_SPIN_ADAPT = 0 930 IREFSPC = IREFSPCE 931 MSCOMB_CC = 0 932* 933C? WRITE(6,*) ' I_INC_AA = ', I_INC_AA 934 IF(I_INC_AA.EQ.0) THEN 935 NOAAEX = 1 936 ELSE 937 NOAAEX = 0 938 END IF 939*. Form of internal states 940* 1 => diagonalize metric 941* 2 => Diagonalize zero-order Hamiltonian matrix 942* 3 => Diagonalize zero-order Jacobian => gives different left and right 943* zero-order stated 944 I_IN_TP = 2 945*. This is using splitting into internals and externals so 946 I_DO_EI = 1 947* Internal inactive -> inactive, sec -> sec will be eliminated 948 IF(ICEXC_INT.NE.0) THEN 949 WRITE(6,*) 950 & ' inactive -> inactive, sec -> sec excitations turned off' 951 ICEXC_INT = 0 952 END IF 953* 954 NTEST = 2 955 IPRNT = NTEST 956 IF(NTEST.GE.1) THEN 957 WRITE(6,*) 958 WRITE(6,*) ' Internal contracted section entered ' 959 WRITE(6,*) ' ==================================== ' 960 WRITE(6,*) 961 WRITE(6,*) ' Version exploiting external/internal division' 962 WRITE(6,*) 963 WRITE(6,*) 964 WRITE(6,'(A,A)') ' Form of calculation ', ICTYP 965 WRITE(6,*) 966 WRITE(6,*) ' Symmetri of reference vector ' , IREFSM 967 WRITE(6,*) ' Space of Reference vector ', IREFSPC 968 WRITE(6,*) ' Space of Internal contracted vector ', ITREFSPC-1 969 WRITE(6,*) 970 WRITE(6,*) ' Parameters defining operator manifold:' 971 WRITE(6,*) ' Max operator rank ', ICOP_RANK_MAX 972 WRITE(6,*) ' Max number of active indeces ', 973 & ICEXC_MAX_ACT 974 WRITE(6,*) ' Max number of external indeces ', 975 & ICEXC_MAX_EXT 976 IF(ICEXC_INT.EQ.1) THEN 977 WRITE(6,*) ' ', 978 & 'Internal (ina->ina, sec->sec) excitations allowed' 979 ELSE 980 WRITE(6,*) ' ', 981 & 'Internal (ina->ina, sec->sec) excitations not allowed' 982 END IF 983 IF(NOAAEX.EQ.1) THEN 984 WRITE(6,*) ' ', 985 & 'Pure active excitations are not included' 986 ELSE 987 WRITE(6,*) ' ', 988 & ' Pure active excitations are included' 989 END IF 990 WRITE(6,*) 991 WRITE(6,*) 992 & ' Threshold for nonsingular eigenvalues of metric', 993 & THRES_SINGU 994* 995 IF(IRESTRT_IC.EQ.1) THEN 996 WRITE(6,*) ' Restarted calculation : ' 997 WRITE(6,*) ' IC coefficients read from LUSC54' 998 WRITE(6,*) ' CI for reference read from LUEXC ' 999 END IF 1000 WRITE(6,*) ' Zero-order states obtained by:' 1001 IF(I_IN_TP.EQ.1) THEN 1002 WRITE(6,*) ' diagonalizing metric ' 1003 ELSE IF(I_IN_TP.EQ.2) THEN 1004 WRITE(6,*) ' Diagonalizing zero-order Hamiltonian matrix' 1005 ELSE IF (I_IN_TP.EQ.3) THEN 1006 WRITE(6,*) ' Diagonalizing zero-order Jacobian matrix' 1007 END IF 1008 END IF 1009* 1010*. Divide orbital spaces into inactive, active, secondary using 1011*. space 1 1012 CALL CC_AC_SPACES(1,IREFTYP) 1013*. Obtain the orbital excitations 1014* (copied more or less from LUCIA_GENCC) 1015 IATP = 1 1016 IBTP = 2 1017* 1018 NAEL = NELEC(IATP) 1019 NBEL = NELEC(IBTP) 1020 NEL = NAEL + NBEL 1021* 1022COLD ICSPC = ICISPC 1023COLD ISSPC = ICISPC 1024* 1025C? WRITE(6,*) ' Zero-order Hamiltonian with zero-order density ' 1026C. IPHGAS1 should be used to divide into H,P,V, but IPHGAS is used, so swap 1027 CALL ISWPVE(IPHGAS(1),IPHGAS1(1),NGAS) 1028 CALL COPVEC(WORK(KINT1O),WORK(KFIFA),NINT1) 1029 CALL FIFAM(WORK(KFIFA)) 1030 IF(NTEST.GE.10) THEN 1031 WRITE(6,*) ' FIFA: ' 1032 CALL APRBLM2(WORK(KFIFA),NTOOBS,NTOOBS,NSMOB,1) 1033 END IF 1034*. External part of zero-order energy: External part is assumed to 1035*. be doubly occupied hole space and unoccupied particle space 1036C GET_E0REF_EXT(FI,IPHGASX) 1037C E0_REF_EXT = GET_E0REF_EXT(WORK(KFIFA),IPHGAS) 1038*. And clean up 1039 CALL ISWPVE(IPHGAS,IPHGAS1,NGAS) 1040 1041* 1042* 1043* ======================== 1044* info for reference space 1045* ======================== 1046* 1047*. Make sure that there is just a single occupation space 1048 CALL OCCLSE(1,NOCCLS_REF,IOCCLS,NEL,IREFSPC,0,0,NOBPT) 1049 IF(NOCCLS_REF.NE.1) THEN 1050 WRITE(6,*) ' Problem in LUCIA_IC_NEW : ' 1051 WRITE(6,*) 1052 & ' Reference space is not a single occupation space' 1053 STOP 1054 & ' Reference space is not a single occupation space' 1055 END IF 1056*. and the reference occupation space 1057 CALL MEMMAN(KLOCCLS_REF,NGAS,'ADDL ',1,'OCC_RF') 1058 CALL OCCLSE(2,NOCCLS_REF,WORK(KLOCCLS_REF),NEL,IREFSPC,0,0,NOBPT) 1059* 1060* ==================================== 1061* Info for space defining excitations 1062* ==================================== 1063* 1064*. Number 1065C IT2REFSPC = ITREFSPC - 1 1066 IT2REFSPC = ITREFSPC 1067 CALL OCCLSE(1,NOCCLS,IOCCLS,NEL,IT2REFSPC,0,0,NOBPT) 1068*. And the occupation classes 1069 CALL MEMMAN(KLOCCLS,NOCCLS*NGAS,'ADDL ',1,'OCCLS ') 1070 CALL OCCLSE(2,NOCCLS,WORK(KLOCCLS),NEL,IT2REFSPC,0,0,NOBPT) 1071*. Number of occupation classes for T-operators 1072 NTOCCLS = NOCCLS 1073*. It could be an idea to check that reference space is included 1074* ======================== 1075* Orbital excitation types 1076* ======================== 1077* 1078*. Number of excitation types 1079 IFLAG = 1 1080 IDUM = 1 1081* 1082 MX_NCREA = ICOP_RANK_MAX 1083 MX_NANNI = ICOP_RANK_MAX 1084 MX_AAEXC = ICEXC_MAX_ACT 1085 I_OOCC = 0 1086*. Pure AA excitations (without external part?) 1087 IF(I_INC_AA.EQ.0) THEN 1088 NOAAEX = 1 1089 ELSE 1090 NOAAEX = 0 1091 END IF 1092C? WRITE(6,*) ' NOAAEX =', NOAAEX 1093 IFLAG = 1 1094 CALL TP_OBEX3(NOCCLS,NEL,NGAS,WORK(IDUM), 1095 & WORK(IDUM),WORK(IDUM), 1096 & WORK(KLOCCLS),WORK(KLOCCLS_REF),MX_NCREA,MX_NANNI, 1097 & MX_EXC_LEVEL,WORK(IDUM),MX_AAEXC,IFLAG, 1098 & I_OOCC,NOBEX_TP,NOAAEX,ICEXC_MAX_EXT,IPRCC) 1099 IF(IPRNT.GE.5) 1100 &WRITE(6,*) ' NOBEX_TP,MX_EXC_LEVEL = ', NOBEX_TP,MX_EXC_LEVEL 1101*. And the actual orbital excitations 1102*. An orbital excition operator is defined by 1103* 1 : Number of creation operators 1104* 2 : Number of annihilation operators 1105* 3 : The actual creation and annihilation operators 1106*. The number of orbital excitations is increased by one to include 1107*. excitations within the reference space 1108 NOBEX_TPE = NOBEX_TP+1 1109 CALL MEMMAN(KLCOBEX_TP,NOBEX_TPE,'ADDL ',1,'LCOBEX') 1110 CALL MEMMAN(KLAOBEX_TP,NOBEX_TPE,'ADDL ',1,'LAOBEX') 1111 CALL MEMMAN(KOBEX_TP ,NOBEX_TPE*2*NGAS,'ADDL ',1,'IOBE_X') 1112*. Excitation type => Original occupation class 1113 CALL MEMMAN(KEX_TO_OC,NOBEX_TPE,'ADDL ',1,'EX__OC') 1114 IFLAG = 0 1115*. (Unit operator is created even if only NOBEX_TP is transferred 1116 CALL TP_OBEX3(NOCCLS,NEL,NGAS,WORK(KOBEX_TP), 1117 & WORK(KLCOBEX_TP),WORK(KLAOBEX_TP), 1118 & WORK(KLOCCLS),WORK(KLOCCLS_REF),MX_NCREA,MX_NANNI, 1119 & MX_EXC_LEVEL,WORK(KEX_TO_OC),MX_AAEXC,IFLAG, 1120 & I_OOCC,NOBEX_TP,NOAAEX,ICEXC_MAX_EXT,IPRCC) 1121 1122* 1123* ======================= 1124* Spinorbital excitations 1125* ======================= 1126* 1127*. Spin combinations of CC excitations : Currently we assume that 1128*. The T-operator is a singlet, can 'easily' be changed 1129* 1130*. Notice : The first time in OBEX_TO_SPOBEX we always use MSCOMB_CC = 0. 1131*. This may lead to the allocation of too much space for 1132*. spinorbital excitations, but MSCOMB_CC = 1, requires access 1133*. to WORK(KLSOBEX) which has not been defined 1134* 1135* 1136* Combined external-internal excitations 1137*. Largest spin-orbital excitation level 1138 IF(MXSPOX.NE.0) THEN 1139 MXSPOX_L = MXSPOX 1140 ELSE 1141 MXSPOX_L = MX_EXC_LEVEL 1142 END IF 1143 IF(NTEST.GE.10) 1144 &WRITE(6,*) ' MXSPOX, MXSPOX_L, MX_EXC_LEVEL = ', 1145 & MXSPOX, MXSPOX_L, MX_EXC_LEVEL 1146 IZERO = 0 1147 CALL OBEX_TO_SPOBEX(1,WORK(KOBEX_TP),WORK(KLCOBEX_TP), 1148 & WORK(KLAOBEX_TP),NOBEX_TP,IDUMMY,NSPOBEX_TP,NGAS, 1149 & NOBPT,0,IZERO ,IAAEXC_TYP,IACT_SPC,IPRCC,IDUMMY, 1150 & MXSPOX_L,WORK(KNSOX_FOR_OX), 1151 & WORK(KIBSOX_FOR_OX),WORK(KISOX_FOR_OX),NAEL,NBEL,IREFSPC) 1152*. Extended number of spin-orbital excitations : Include 1153*. unit operator as last spinorbital excitation operator 1154 NSPOBEX_TPE = NSPOBEX_TP + 1 1155 IF(IPRNT.GE.1) THEN 1156 WRITE(6,*) ' Number of spinorbital excitations(with unit)', 1157 & NSPOBEX_TPE 1158 END IF 1159 IF(IPRNT.GE.10) WRITE(6,*) ' NSPOBEX_TP = ', NSPOBEX_TP 1160*. Allocate space for EI, E, I (external-internal, external, internal ) 1161*. spinorbital excitations. As the number of E and I are not known 1162*. we set their dimension to NSPOBEX_TP ( an upper limit) 1163*. And the actual spinorbital excitation operators 1164 CALL MEMMAN(KLSOBEX,4*NGAS*3*NSPOBEX_TPE,'ADDL ',1,'SPOBEX') 1165*. Map spin-orbital exc type => orbital exc type 1166 CALL MEMMAN(KLSOX_TO_OX,3*NSPOBEX_TPE,'ADDL ',1,'SPOBEX') 1167*. First SOX of given OX ( including zero operator ) 1168 CALL MEMMAN(KIBSOX_FOR_OX,NOBEX_TP+1,'ADDL ',1,'IBSOXF') 1169*. Number of SOX's for given OX 1170 CALL MEMMAN(KNSOX_FOR_OX,NOBEX_TP+1,'ADDL ',1,'IBSOXF') 1171*. SOX for given OX 1172 CALL MEMMAN(KISOX_FOR_OX,NSPOBEX_TP+1,'ADDL ',1,'IBSOXF') 1173 CALL OBEX_TO_SPOBEX(2,WORK(KOBEX_TP),WORK(KLCOBEX_TP), 1174 & WORK(KLAOBEX_TP),NOBEX_TP,WORK(KLSOBEX),NSPOBEX_TP,NGAS, 1175 & NOBPT,0,0,IAAEXC_TYP,IACT_SPC, 1176 & IPRCC,WORK(KLSOX_TO_OX),MXSPOX_L,WORK(KNSOX_FOR_OX), 1177 & WORK(KIBSOX_FOR_OX),WORK(KISOX_FOR_OX),NAEL,NBEL,IREFSPC) 1178 NSPOBEX_TPE = NSPOBEX_TP + 1 1179*. Add unit-operator as last spinorbital excitation 1180 IZERO = 0 1181 CALL ISTVC3(WORK(KLSOBEX),(NSPOBEX_TPE-1)*4*NGAS+1,IZERO,4*NGAS) 1182 IF(IPRNT.GE.5) THEN 1183 WRITE(6,*) ' Extended list of spin-orbital excitations : ' 1184 CALL WRT_SPOX_TP(WORK(KLSOBEX),NSPOBEX_TPE) 1185 END IF 1186* 1187*. Construct the various external and internal operators of the 1188*. various operatortpys and set up mappings from IE operator 1189*. to I,E operators 1190*. Question Jeppe : Should the zero-operator also be splitted. 1191*. Yes, I am doing this from today -aug20 1192 I_EXT_OFF = NSPOBEX_TPE + 1 1193*. Offset in KLSOBEX to internal part is obtained in SPLIT_IE_SPOBEXTP 1194 I_INT_OFF = 0 1195*. The various internal operators for the same external operators 1196*. will be collected. Mappings for this 1197*. Number of internal types per external type 1198 CALL MEMMAN(KL_N_INT_FOR_EXT,NSPOBEX_TPE,'ADDL ',1,'N_INT_') 1199*. Offsets for internals for given external 1200 CALL MEMMAN(KL_IB_INT_FOR_EXT,NSPOBEX_TPE,'ADDL ',1,'IB_INT') 1201*. And the actual internals for each external 1202 CALL MEMMAN(KL_I_INT_FOR_EXT,NSPOBEX_TPE,'ADDL ',1,'I_INT_') 1203* 1204 CALL SPLIT_IE_SPOBEXTP(WORK(KLSOBEX),NSPOBEX_TPE, N_EXTOP_TP, 1205 & N_INTOP_TP,I_EXT_OFF, I_INT_OFF,WORK(KL_N_INT_FOR_EXT), 1206 & WORK(KL_IB_INT_FOR_EXT), WORK(KL_I_INT_FOR_EXT),NGAS, 1207 & IHPVGAS ) 1208*. Obtain reorder array EI-order => standard order 1209 CALL MEMMAN(KL_I_EI_TP_REO,NSPOBEX_TPE,'ADDL ',1,'EITPRE') 1210C? WRITE(6,*) ' NSPOBEX_TPE before EITP.. =', NSPOBEX_TPE 1211 CALL EITP_TO_SPOXTP_REO(WORK(KL_I_EI_TP_REO), 1212 & WORK(KLSOBEX),NSPOBEX_TPE, 1213 & I_EXT_OFF,I_INT_OFF,NGAS,N_EXTOP_TP,N_INTOP_TP, 1214 & WORK(KL_N_INT_FOR_EXT), WORK(KL_IB_INT_FOR_EXT), 1215 & WORK(KL_I_INT_FOR_EXT) ) 1216C? WRITE(6,*) ' I_INT_OFF after sec call to EITP_TO.. ', I_INT_OFF 1217C EITP_TO_SPOXTP_REO(I_EI_TP_REO,ISPOBEX_TP,NSPOBEX_TP, 1218C & IB_EXTP, IB_INTP,NGAS,N_EXTP,N_INTP, 1219C & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX) 1220*. We have now stored information about the external types 1221*. in WORK(KLSOBEX) from I_EXT_OFF and for external types in 1222*. I_INT_OFF. We do however not increase NSPOBEX_TP, 1223*. so in the following there are hopefully invisible 1224*. 1225*. Obtain info on the dimension of the various internal and external types 1226* 1227C DIMENSION_EI_EXP(ISPOBEX_TP,IB_EXTP,IB_INTP,N_EXTP, 1228C & N_INTP,N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX, 1229C & NDIM_EX_ST,NDIM_IN_ST,NDIM_EI,NDIM_IN_SE, NSMOB) 1230 CALL MEMMAN(KL_NDIM_EX_ST,NSMOB*N_EXTOP_TP,'ADDL ',1,'NDIM_E') 1231 CALL MEMMAN(KL_NDIM_IN_ST,NSMOB*N_INTOP_TP,'ADDL ',1,'NDIM_I') 1232 CALL MEMMAN(KL_NDIM_IN_SE,NSMOB*N_EXTOP_TP,'ADDL ',1,'NDIMSE') 1233 CALL DIMENSION_EI_EXP(WORK(KLSOBEX),I_EXT_OFF,I_INT_OFF, 1234 & N_EXTOP_TP,N_INTOP_TP,WORK(KL_N_INT_FOR_EXT), 1235 & WORK(KL_IB_INT_FOR_EXT), WORK(KL_I_INT_FOR_EXT), 1236 & WORK(KL_NDIM_EX_ST),WORK(KL_NDIM_IN_ST),NDIM_EI, 1237 & WORK(KL_NDIM_IN_SE),NSMOB) 1238* 1239 CALL ISTVC3(WORK(KLSOX_TO_OX),NSPOBEX_TPE,NOBEX_TP+1,1) 1240*. Mapping spinorbital excitations => occupation classes 1241 CALL MEMMAN(KIBSOX_FOR_OCCLS,NOCCLS,'ADDL ',1,'IBSXOC') 1242 CALL MEMMAN(KNSOX_FOR_OCCLS,NOCCLS,'ADDL ',1,' NSXOC') 1243 CALL MEMMAN(KISOX_FOR_OCCLS,NSPOBEX_TPE,'ADDL ',1,' ISXOC') 1244C SPOBEX_FOR_OCCLS( 1245C & IEXTP_TO_OCCLS,NOCCLS,ISOX_TO_OX,NSOX, 1246C & NSOX_FOR_OCCLS,ISOX_FOR_OCCLS,IBSOX_FOR_OCCLS) 1247 CALL SPOBEX_FOR_OCCLS(WORK(KEX_TO_OC),NOCCLS,WORK(KLSOX_TO_OX), 1248 & NSPOBEX_TPE,WORK(KNSOX_FOR_OCCLS),WORK(KISOX_FOR_OCCLS), 1249 & WORK(KIBSOX_FOR_OCCLS)) 1250* 1251*. Frozen spin-orbital excitation types 1252 CALL MEMMAN(KLSPOBEX_FRZ, NSPOBEX_TPE,'ADDL ',1,'SPOBFR') 1253 CALL FRZ_SPOBEX(WORK(KLSPOBEX_FRZ),WORK(KLCOBEX_TP),NSPOBEX_TP, 1254 & WORK(KLSOX_TO_OX),IFRZ_CC_AR,NFRZ_CC) 1255 IZERO = 0 1256 CALL ISTVC3(WORK(KLSPOBEX_FRZ),NSPOBEX_TPE,IZERO,1) 1257*. Spin-orbital excitation types related by spin-flip 1258 CALL MEMMAN(KLSPOBEX_AB,NSPOBEX_TPE,'ADDL ',1,'SPOBAB') 1259 CALL SPOBEXTP_PAIRS(NSPOBEX_TPE,WORK(KLSOBEX),NGAS, 1260 & WORK(KLSPOBEX_AB)) 1261C SPOBEXTP_PAIRS(NSPOBEX_TP,ISPOBEX,NGAS,ISPOBEX_PAIRS) 1262C SELECT_AB_TYPES(NSPOBEX_TP,ISPOBEX_TP, 1263C & ISPOBEX_PAIRS,NGAS) 1264 CALL SELECT_AB_TYPES(NSPOBEX_TPE,WORK(KLSOBEX), 1265 & WORK(KLSPOBEX_AB),NGAS) 1266*. Alpha- and beta-excitations constituting the spinorbital excitations 1267*. Number 1268 CALL SPOBEX_TO_ABOBEX(WORK(KLSOBEX),NSPOBEX_TP,NGAS, 1269 & 1,NAOBEX_TP,NBOBEX_TP,IDUMMY,IDUMMY) 1270*. And the alpha-and beta-excitations 1271 LENA = 2*NGAS*NAOBEX_TP 1272 LENB = 2*NGAS*NBOBEX_TP 1273 CALL MEMMAN(KLAOBEX,LENA,'ADDL ',2,'IAOBEX') 1274 CALL MEMMAN(KLBOBEX,LENB,'ADDL ',2,'IAOBEX') 1275 CALL SPOBEX_TO_ABOBEX(WORK(KLSOBEX),NSPOBEX_TP,NGAS, 1276 & 0,NAOBEX_TP,NBOBEX_TP,WORK(KLAOBEX),WORK(KLBOBEX)) 1277*. Max dimensions of CCOP !KSTR> = !ISTR> maps 1278*. For alpha excitations 1279 IATP = 1 1280 IOCTPA = IBSPGPFTP(IATP) 1281 NOCTPA = NSPGPFTP(IATP) 1282 CALL LEN_GENOP_STR_MAP( 1283 & NAOBEX_TP,WORK(KLAOBEX),NOCTPA,NELFSPGP(1,IOCTPA), 1284 & NOBPT,NGAS,MAXLENA) 1285 IBTP = 2 1286 IOCTPB = IBSPGPFTP(IBTP) 1287 NOCTPB = NSPGPFTP(IBTP) 1288 CALL LEN_GENOP_STR_MAP( 1289 & NBOBEX_TP,WORK(KLBOBEX),NOCTPB,NELFSPGP(1,IOCTPB), 1290 & NOBPT,NGAS,MAXLENB) 1291 MAXLEN_I1 = MAX(MAXLENA,MAXLENB) 1292 IF(IPRNT.GE.10) WRITE(6,*) ' MAXLEN_I1 = ', MAXLEN_I1 1293* 1294* Max Dimension of spinorbital excitation operators 1295* 1296 CALL MEMMAN(KLLSOBEX,NSPOBEX_TPE,'ADDL ',1,'LSPOBX') 1297 CALL MEMMAN(KLIBSOBEX,NSPOBEX_TPE,'ADDL ',1,'LSPOBX') 1298 CALL MEMMAN(KLSPOBEX_AC,NSPOBEX_TPE,'ADDL ',1,'SPOBAC') 1299*. ALl spinorbital excitations are initially active 1300 IONE = 1 1301 CALL ISETVC(WORK(KLSPOBEX_AC),IONE,NSPOBEX_TPE) 1302* 1303 MX_ST_TSOSO_MX = 0 1304 MX_ST_TSOSO_BLK_MX = 0 1305 MX_TBLK_MX = 0 1306 MX_TBLK_AS_MX = 0 1307 LEN_T_VEC_MX = 0 1308 DO ICCAMP_SM = 1, NSMST 1309*. Dimension without zero-particle operator 1310 CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,ICCAMP_SM, 1311 & MX_ST_TSOSOL,MX_ST_TSOSO_BLKL,MX_TBLKL, 1312 & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VECL, 1313 & MSCOMB_CC,MX_TBLK_AS, 1314 & WORK(KISOX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 1315 & NTCONF,IPRCC) 1316* 1317 MX_ST_TSOSO_MX = MAX(MX_ST_TSOSO_MX,MX_ST_TSOSOL) 1318 MX_ST_TSOSO_BLK_MX = MAX(MX_ST_TSOSO_BLK_MX,MX_ST_TSOSO_BLKL) 1319 MX_TBLK_MX = MAX(MX_TBLK_MX,MX_TBLKL) 1320 MX_TBLK_AS_MX = MAX(MX_TBLK_AS_MX,MX_TBLK_AS) 1321 LEN_T_VEC_MX = MAX(LEN_T_VEC_MX, LEN_T_VECL) 1322* 1323 END DO 1324 IF(IPRNT.GE.10) WRITE(6,*) ' MX_TBLK_AS_MX = ', MX_TBLK_AS_MX 1325 IF(IPRNT.GE.10) WRITE(6,*) ' LEN_T_VEC_MX = ', LEN_T_VEC_MX 1326*. And dimensions for symmetry 1 1327 ITOP_SM = 1 1328 CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,ITOP_SM, 1329 & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 1330 & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VEC, 1331 & MSCOMB_CC,MX_SBSTR, 1332 & WORK(KISOX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 1333 & NTCONF,IPRCC) 1334*. NCAAB and N_CC_AMP does not include zero-particle operator 1335 N_CC_AMP = LEN_T_VEC 1336 NCAAB = N_CC_AMP 1337 IF(IPRNT.GE.5) 1338 &WRITE(6,*) ' LUCIA_GENCC : N_CC_AMP = ', N_CC_AMP 1339 IF(IPRNT.GE.5) THEN 1340 WRITE(6,*) ' Number of amplitudes per operator type: ' 1341 CALL IWRTMA(WORK(KLLSOBEX),NSPOBEX_TP,1,NSPOBEX_TP,1) 1342 END IF 1343*. Hard wire info for unit operator stored as last spinorbital excitation 1344C ISTVC2(IVEC,IBASE,IFACT,NDIM) 1345 IONE = 1 1346 CALL ISTVC3(WORK(KLLSOBEX),NSPOBEX_TPE,IONE,1) 1347 N_CC_AMPP1 = N_CC_AMP + 1 1348 CALL ISTVC3(WORK(KLIBSOBEX),NSPOBEX_TPE,N_CC_AMPP1,1) 1349*. Obtain mapping between EI order and standard order 1350 CALL MEMMAN(KL_IREO_EI_ST,N_CC_AMP+1,'ADDL ',1,'IREOST') 1351C? WRITE(6,*) ' I_EXT_OFF, I_INT_OFF = ', I_EXT_OFF,I_INT_OFF 1352 CALL EI_REORDER_ARRAYS(WORK(KL_IREO_EI_ST),WORK(KLSOBEX), 1353 & I_EXT_OFF,I_INT_OFF,WORK(KLLSOBEX),WORK(KLIBSOBEX), 1354 & NSPOBEX_TPE,N_EXTOP_TP,N_INTOP_TP, 1355 & WORK(KL_N_INT_FOR_EXT),WORK(KL_IB_INT_FOR_EXT), 1356 & WORK(KL_I_INT_FOR_EXT),1) 1357C EI_REORDER_ARRAYS(IREO_EI_ST, 1358C & ISPOBEX_TP,IB_EXTP,IB_INTP, 1359C & L_SPOBEX_TP,IB_SPOBEX_TP,NSPOBEX_TP, 1360C & N_EXTP,N_INTP, 1361C & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX,ISYM) 1362*. Determine dimensions of the various internal and external 1363*. types 1364* ============= 1365* Scratch space 1366* ============= 1367* 1368 CALL GET_3BLKS_GCC(KVEC1,KVEC2,KVEC3,MXCJ) 1369 KVEC1P = KVEC1 1370 KVEC2P = KVEC2 1371* 1372 IF(NTEST.GE.100) THEN 1373 WRITE(6,*) ' Standard list of CAAB operators' 1374 CALL PRINT_CAAB_LIST(1) 1375 END IF 1376* 1377* 1378* -------------------------- 1379* Obtain the internal states 1380* -------------------------- 1381* 1382 CALL MEMMAN(KL_N_ORTN_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1,'ORT_SE') 1383 CALL MEMMAN(KL_N_INT_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1,'INT_SE') 1384*. Dimension of matrices with given symmetry and external types 1385C? WRITE(6,*) ' NDIM_IN_SE array : ' 1386C? CALL IWRTMA(WORK(KL_NDIM_IN_SE), 1387C? & NSMOB,N_EXTOP_TP,NSMOB,N_EXTOP_TP) 1388* 1389 I_SE_DIM2 = ISQELSUM(WORK(KL_NDIM_IN_SE),NSMOB*N_EXTOP_TP,0) 1390 WRITE(6,*) ' Dimension of transformation matrices X1, X2', 1391 & I_SE_DIM2 1392 I_SE_DIM2S= ISQELSUM(WORK(KL_NDIM_IN_SE),NSMOB*N_EXTOP_TP,1) 1393 WRITE(6,*) ' Dimension of internal overlap matrix ', 1394 & I_SE_DIM2S 1395 I_SE_DIM = IELSUM(WORK(KL_NDIM_IN_SE),NSMOB*N_EXTOP_TP) 1396*. Space for transformation matrices: For each block matrices 1397*. X1, X2 and a vector sigma 1398 CALL MEMMAN(KL_X1_INT_EI_FOR_SE,I_SE_DIM2,'ADDL ',2,'X1INEI') 1399 CALL MEMMAN(KL_X2_INT_EI_FOR_SE,I_SE_DIM2,'ADDL ',2,'X2INEI') 1400 CALL MEMMAN(KL_SG_INT_EI_FOR_SE,I_SE_DIM ,'ADDL ',2,'SGINEI') 1401 CALL MEMMAN(KL_S_INT_EI_FOR_SE,I_SE_DIM2S,'ADDL ',2,'SINEI ') 1402 CALL MEMMAN(KL_IBX1_INT_EI_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1, 1403 & 'IBX1IN') 1404 CALL MEMMAN(KL_IBX2_INT_EI_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1, 1405 & 'IBX2IN') 1406 CALL MEMMAN(KL_IBSG_INT_EI_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1, 1407 & 'IBSGIN') 1408 CALL MEMMAN(KL_IBS_INT_EI_FOR_SE,NSMOB*N_EXTOP_TP,'ADDL ',1, 1409 & 'IBSIN ') 1410*. For left hand side transformation to internal state basis 1411*. if Jacobian is diagonalized 1412*. matrix diagonalizing metric is common for L and R as is sigma, 1413*. so only an extra matrix for the last diagonalization is needed 1414 IF(I_IN_TP.LE.2) THEN 1415 KL_X2L_INT_EI_FOR_SE = KL_X2_INT_EI_FOR_SE 1416 ELSE 1417 CALL MEMMAN(KL_X2L_INT_EI_FOR_SE,I_SE_DIM2,'ADDL ',2,'XRINEI') 1418 END IF 1419* 1420*.--------------------------------------------- 1421*. Obtain internal states by diagonalizing H0 1422*.--------------------------------------------- 1423C? WRITE(6,*) ' Call to IDIM_TCC just before zero-order states' 1424C? CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,1, 1425C? & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 1426C? & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VEC, 1427C? & MSCOMB_CC,MX_SBSTR, 1428C? & WORK(KISOX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 1429C? & NTCONF,IPRCC) 1430* 1431C? WRITE(6,*) ' Before GET_INT, I_INT_OFF,I_EXT_OFF ', 1432C? & I_INT_OFF,I_EXT_OFF 1433 CALL GET_INTERNAL_STATES(N_EXTOP_TP,N_INTOP_TP, 1434 & WORK(KLSOBEX),WORK(KL_N_INT_FOR_EXT),WORK(KL_IB_INT_FOR_EXT), 1435 & WORK(KL_I_INT_FOR_EXT),WORK(KL_NDIM_IN_SE), 1436 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 1437 & WORK(KL_X1_INT_EI_FOR_SE), WORK(KL_X2_INT_EI_FOR_SE), 1438 & WORK(KL_SG_INT_EI_FOR_SE),WORK(KL_S_INT_EI_FOR_SE), 1439 & WORK(KL_IBX1_INT_EI_FOR_SE), WORK(KL_IBX2_INT_EI_FOR_SE), 1440 & WORK(KL_IBSG_INT_EI_FOR_SE),WORK(KL_IBS_INT_EI_FOR_SE), 1441 & WORK(KL_X2L_INT_EI_FOR_SE), 1442 & I_IN_TP,I_INT_OFF,I_EXT_OFF) 1443* 1444 IF(NTEST.GE.10) THEN 1445 WRITE(6,*) 1446 & ' Number of internal states per sym and ext type' 1447 CALL IWRTMA(WORK(KL_N_INT_FOR_SE), 1448 & NSMOB,N_EXTOP_TP, NSMOB,N_EXTOP_TP) 1449 WRITE(6,*) 1450 & ' Number of orthn internal states per sym and ext type' 1451 CALL IWRTMA(WORK(KL_N_ORTN_FOR_SE), 1452 & NSMOB,N_EXTOP_TP, NSMOB,N_EXTOP_TP) 1453 END IF 1454*. Largest number of internal states for given sym and external type 1455C IMNNMX(IVEC,NDIM,MINMAX) 1456 N_INT_MAX = IMNMX(WORK(KL_N_INT_FOR_SE),N_EXTOP_TP*NSMOB,2) 1457*. Largest number of zero-order states of given sym and external type 1458 N_ORTN_MAX = IMNMX(WORK(KL_N_ORTN_FOR_SE),N_EXTOP_TP*NSMOB,2) 1459 WRITE(6,*) ' N_INT_MAX, N_ORTN_MAX = ', N_INT_MAX, N_ORTN_MAX 1460*. Largest transformation block 1461 N_XEO_MAX = N_INT_MAX*N_ORTN_MAX 1462 IF(NTEST.GE.5) WRITE(6,*) ' Largest (EL,ORTN) block = ', N_XEO_MAX 1463* 1464*. Number of zero-order states - does now include the unit-operator 1465 N_ZERO_EI = N_ZERO_ORDER_STATES(WORK(KL_N_ORTN_FOR_SE), 1466 & WORK(KL_NDIM_EX_ST),N_EXTOP_TP,1) 1467 WRITE(6,*) ' Number of zero-order states with sym 1 = ', N_ZERO_EI 1468C N_ZERO_ORDER_STATES(NORTN_FOR_SE,NDIM_EX_ST,N_EXTP,ITOTSYM) 1469C? WRITE(6,*) ' Call to IDIM_TCC just after zero-order states' 1470C? CALL IDIM_TCC(WORK(KLSOBEX),NSPOBEX_TP,1, 1471C? & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 1472C? & WORK(KLLSOBEX),WORK(KLIBSOBEX),LEN_T_VEC, 1473C? & MSCOMB_CC,MX_SBSTR, 1474C? & WORK(KISOX_FOR_OCCLS),NOCCLS,WORK(KIBSOX_FOR_OCCLS), 1475C? & NTCONF,IPRCC) 1476* ============= 1477* Scratch space 1478* ============= 1479* 1480*. Scratch space for CI - behind the curtain 1481COLD CALL GET_3BLKS_GCC(KVEC1,KVEC2,KVEC3,MXCJ) 1482*. Pointers to KVEC1 and KVEC2, transferred through GLBBAS 1483COLD KVEC1P = KVEC1 1484COLD KVEC2P = KVEC2 1485C? WRITE(6,*) ' KVEC3 after GET_3BLKS.. ', KVEC3 1486*. and two CC vectors , extra element for unexcited SD at end of vectors 1487 N_SD_INT = 1 1488 LENNY = LEN_T_VEC_MX + N_SD_INT 1489 if (i_obcc.eq.1.or.i_oocc.eq.1.or.i_bcc.eq.1) then 1490C! lenny = max(lenny,nooexc(1)+nooexc(2)) 1491 STOP ' Jeppe copied out (nooexc not defined)' 1492 end if 1493 CALL MEMMAN(KCC1,LENNY,'ADDL ',2,'CC1_VE') 1494 CALL MEMMAN(KCC2,LENNY,'ADDL ',2,'CC2_VE') 1495*. And the CC diagonal 1496 CALL MEMMAN(KDIA,LENNY,'ADDL ',2,'CC_DIA') 1497* 1498*. Max dimensions of CCOP !KSTR> = !ISTR> maps 1499*. For alpha excitations 1500 IATP = 1 1501 IOCTPA = IBSPGPFTP(IATP) 1502 NOCTPA = NSPGPFTP(IATP) 1503 CALL LEN_GENOP_STR_MAP( 1504 & NAOBEX_TP,WORK(KLAOBEX),NOCTPA,NELFSPGP(1,IOCTPA), 1505 & NOBPT,NGAS,MAXLENA) 1506 IBTP = 2 1507 IOCTPB = IBSPGPFTP(IBTP) 1508 NOCTPB = NSPGPFTP(IBTP) 1509 CALL LEN_GENOP_STR_MAP( 1510 & NBOBEX_TP,WORK(KLBOBEX),NOCTPB,NELFSPGP(1,IOCTPB), 1511 & NOBPT,NGAS,MAXLENB) 1512 MAXLEN_I1 = MAX(MAXLENA,MAXLENB) 1513 IF(NTEST.GE.5) WRITE(6,*) ' MAXLEN_I1 = ', MAXLEN_I1 1514* 1515 IF(NTEST.GE.100) CALL PRINT_ZERO_TRMAT 1516* 1517 IF(ICTYP(1:4).EQ.'ICCI') THEN 1518* 1519* ============================== 1520* Internal contracted CI section 1521* ============================== 1522* 1523 CALL LUCIA_ICCI(IREFSPC,ITREFSPC,ICTYP,EREF, 1524 & EFINAL,CONVER,VNFINAL) 1525* 1526 ELSE IF(ICTYP(1:4).EQ.'ICPT') THEN 1527* 1528* ========================================== 1529* Internal contracted Perturbation expansion 1530* ========================================== 1531* 1532 CALL LUCIA_ICPT(IREFSPC,ITREFSPC,ICTYP,EREF, 1533 & EFINAL,CONVER,VNFINAL) 1534* 1535 ELSE IF(ICTYP(1:4).EQ.'ICCC') THEN 1536* 1537* ====================================== 1538* Internal contracted Coupled Cluster 1539* ======================================= 1540* 1541 CALL LUCIA_ICCC(IREFSPC,ITREFSPC,ICTYP,EREF, 1542 & EFINAL,CONVER,VNFINAL) 1543 END IF 1544* 1545*. 1546 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'ICCI ') 1547* 1548 RETURN 1549 END 1550 SUBROUTINE SPLIT_IE_SPOBEXTP(ISOBEX_TP, NSPOBEX_TP, 1551 & N_EXTOP_TP, N_INTOP_TP,IB_EXTOP_TP,IB_INTOP_TP, 1552 & N_INT_FOR_EXT, IB_INT_FOR_EXT, I_INT_FOR_EXT, 1553 & NGAS, IHPVGAS) 1554*. A set of spinorbital excitations is given. 1555*. Split these into internal and external parts 1556*. and store these in ISOBEX_TP 1557*. IHPVGAS is used to identify Internal/external parts 1558* 1559*. Jeppe Olsen, October 2006 1560* 1561*. Input : 1562* ISOBEX_TP : The spinorbital types in CAAB form 1563* NSPOBEX_TP : Number of spinorbitaltypes 1564* IB_EXTOP_TP : offset for storing external operators in ISOBEX_TP 1565*. output 1566* N_EXTOP_TP : number of external spinorbitalexcitation types 1567* N_INTOP_TP : number of internal spinorbitalexcitation types 1568* N_INT_FOR_EXT(IEXT) : number of internal types for external type IEXT 1569* I_INT_FOR_EXT() : gives the internal types for each external 1570* IB_INT_FOR_EXT(IEXT) : Offset for given external types in I_INT_FOR_EXT 1571* 1572* Note the assymmetry : IB_EXTOP_TP is input whereas IB_INTOP_TP is output 1573* 1574* 1575 INCLUDE 'wrkspc.inc' 1576*. Input (output added) 1577 INTEGER ISOBEX_TP(NGAS,4,*), IHPVGAS(NGAS) 1578*. Output 1579 INTEGER N_INT_FOR_EXT(*), IB_INT_FOR_EXT(*), I_INT_FOR_EXT(*) 1580*. Scratch : Internal and external parts 1581 INTEGER I_INT_TP(4*MXPNGAS) 1582 INTEGER I_EXT_TP(4*MXPNGAS) 1583* 1584 NTEST = 100 1585* 1586*. Obtain the various external parts 1587* 1588 IZERO = 0 1589 CALL ISETVC(N_INT_FOR_EXT,IZERO,NSPOBEX_TP) 1590 N_EXTOP_TP = 0 1591 DO I_EI_OP = 1, NSPOBEX_TP 1592*. Split operator into ext and int parts 1593 CALL SPLIT_EIOP(ISOBEX_TP(1,1,I_EI_OP),I_EXT_TP, I_INT_TP, 1594 & NGAS,IHPVGAS) 1595*. Is this a new external operator ? 1596 NEW = 1 1597 IDENT = 0 1598 DO J_EXTOP_TP = 1, N_EXTOP_TP 1599C COMPARE_TWO_INTARRAYS(IA,IB,NAB,IDENT) 1600 CALL COMPARE_TWO_INTARRAYS(I_EXT_TP, 1601 & ISOBEX_TP(1,1,IB_EXTOP_TP-1+J_EXTOP_TP),4*NGAS,IDENT) 1602 IF(IDENT.EQ.1) THEN 1603*. Match with previous type 1604 N_INT_FOR_EXT(J_EXTOP_TP) = N_INT_FOR_EXT(J_EXTOP_TP) + 1 1605 NEW = 0 1606 END IF 1607 END DO 1608 IF(NEW.EQ.1) THEN 1609 N_EXTOP_TP = N_EXTOP_TP + 1 1610 N_INT_FOR_EXT(N_EXTOP_TP) = 1 1611 CALL ICOPVE(I_EXT_TP, 1612 & ISOBEX_TP(1,1,IB_EXTOP_TP-1+N_EXTOP_TP),4*NGAS) 1613 END IF 1614 END DO 1615* 1616 IF(NTEST.GE.100) THEN 1617 WRITE(6,*) ' Number of external spinorbitalexcitation types ', 1618 & N_EXTOP_TP 1619 WRITE(6,*) ' And the external spinorbitalexcitation types : ' 1620 CALL WRT_SPOX_TP(ISOBEX_TP(1,1,IB_EXTOP_TP),N_EXTOP_TP) 1621 END IF 1622*. Offsets for the various internal parts for each external part 1623 IB_INT_FOR_EXT(1) = 1 1624 DO J_EXTOP_TP = 2, N_EXTOP_TP 1625 IB_INT_FOR_EXT(J_EXTOP_TP) = IB_INT_FOR_EXT(J_EXTOP_TP-1) 1626 & + N_INT_FOR_EXT(J_EXTOP_TP-1) 1627 END DO 1628C? WRITE(6,*) ' N_INT_FOR_EXT, IB_INT_FOR_EXT ' 1629C? CALL IWRTMA(N_INT_FOR_EXT,1,N_EXTOP_TP,1,N_EXTOP_TP) 1630C? CALL IWRTMA(IB_INT_FOR_EXT,1,N_EXTOP_TP,1,N_EXTOP_TP) 1631*. Offsets for internal parts 1632 IB_INTOP_TP = IB_EXTOP_TP+N_EXTOP_TP 1633C? WRITE(6,*) ' IB_INTOP_TP, IB_EXTOP_TP, N_EXTOP_TP = ', 1634C? & IB_INTOP_TP, IB_EXTOP_TP, N_EXTOP_TP 1635*. And generate the various internal parts 1636 N_INTOP_TP = 0 1637 CALL ISETVC(N_INT_FOR_EXT,IZERO,NSPOBEX_TP) 1638 DO I_EI_OP = 1, NSPOBEX_TP 1639C? WRITE(6,*) ' Info for I_EI_OP = ', I_EI_OP 1640*. Split operator into ext and int parts 1641 CALL SPLIT_EIOP(ISOBEX_TP(1,1,I_EI_OP),I_EXT_TP, I_INT_TP, 1642 & NGAS,IHPVGAS) 1643*. Type of this external operator 1644 JJ_EXTOP_TP = -1 1645 DO J_EXTOP_TP = 1, N_EXTOP_TP 1646 CALL COMPARE_TWO_INTARRAYS(I_EXT_TP, 1647 & ISOBEX_TP(1,1,IB_EXTOP_TP-1+J_EXTOP_TP),4*NGAS,IDENT) 1648 IF(IDENT.EQ.1) JJ_EXTOP_TP = J_EXTOP_TP 1649 END DO 1650C? WRITE(6,*) 'JJ_EXTOP_TP = ', JJ_EXTOP_TP 1651*. Is internal type new ? 1652 NEW = 1 1653 IDENT = 0 1654 DO J_INTOP_TP = 1, N_INTOP_TP 1655 CALL COMPARE_TWO_INTARRAYS(I_INT_TP, 1656 & ISOBEX_TP(1,1,IB_INTOP_TP-1+J_INTOP_TP),4*NGAS,IDENT) 1657 IF(IDENT.EQ.1) THEN 1658*. Match with previous type 1659 NEW = 0 1660 JJ_INTOP_TP = J_INTOP_TP 1661 END IF 1662C? WRITE(6,*) ' IDENT = ', IDENT 1663 END DO 1664C? WRITE(6,*) ' NEW = ', NEW 1665 IF(NEW.EQ.1) THEN 1666 N_INTOP_TP = N_INTOP_TP + 1 1667 JJ_INTOP_TP = N_INTOP_TP 1668 CALL ICOPVE(I_INT_TP,ISOBEX_TP(1,1,IB_INTOP_TP-1+N_INTOP_TP), 1669 & 4*NGAS) 1670 END IF 1671 N_INT_FOR_EXT(JJ_EXTOP_TP) = N_INT_FOR_EXT(JJ_EXTOP_TP) + 1 1672 I_INT_FOR_EXT(IB_INT_FOR_EXT(JJ_EXTOP_TP)-1 1673 & +N_INT_FOR_EXT(JJ_EXTOP_TP)) = JJ_INTOP_TP 1674C? WRITE(6,*) ' IB_INT_FOR_EXT(JJ_EXTOP_TP)-1 + ... ', 1675C? & IB_INT_FOR_EXT(JJ_EXTOP_TP)-1 + N_INT_FOR_EXT(JJ_EXTOP_TP) 1676C? WRITE(6,*) ' JJ_EXTOP_TP = ', JJ_EXTOP_TP 1677 END DO 1678* 1679 IF(NTEST.GE.100) THEN 1680 WRITE(6,*) ' Number of internal excitationtypes ', N_INTOP_TP 1681 WRITE(6,*) ' And the internal excitationtypes : ' 1682 CALL WRT_SPOX_TP(ISOBEX_TP(1,1,IB_INTOP_TP),N_INTOP_TP) 1683 WRITE(6,*) ' And the various EI types for each E-type' 1684 DO J_EXTOP_TP = 1, N_EXTOP_TP 1685 WRITE(6,*) ' External operator type ', J_EXTOP_TP 1686 WRITE(6,*) ' Number of internal types for this ', 1687 & N_INT_FOR_EXT(J_EXTOP_TP) 1688 WRITE(6,*) ' and the internal operator types ' 1689 LEN = N_INT_FOR_EXT(J_EXTOP_TP) 1690 CALL IWRTMA(I_INT_FOR_EXT(IB_INT_FOR_EXT(J_EXTOP_TP)), 1691 & 1,LEN,1,LEN) 1692 END DO 1693 END IF 1694* 1695 RETURN 1696 END 1697 SUBROUTINE SPLIT_EIOP(I_EIOP,I_EOP,I_IOP,NGAS,IHPVGAS) 1698* 1699* Split operator IEOP into external and internal parts 1700* according to IHPVGAS 1701* 1702*. Jeppe Olsen, October 2006 1703* 1704 INCLUDE 'wrkspc.inc' 1705*. Input 1706 INTEGER I_EIOP(NGAS,4), IHPVGAS(NGAS) 1707*. Output 1708 INTEGER I_IOP(NGAS,4),I_EOP(NGAS,4) 1709* 1710 IZERO = 0 1711 CALL ISETVC(I_EOP,IZERO,4*NGAS) 1712 CALL ISETVC(I_IOP,IZERO,4*NGAS) 1713* 1714 DO IGAS = 1, NGAS 1715 IF(IHPVGAS(IGAS).LE.2) THEN 1716*. External : Secondary or inactive 1717 DO ICAAB = 1, 4 1718 I_EOP(IGAS,ICAAB) = I_EIOP(IGAS,ICAAB) 1719 END DO 1720 ELSE 1721 DO ICAAB = 1, 4 1722 I_IOP(IGAS,ICAAB) = I_EIOP(IGAS,ICAAB) 1723 END DO 1724 END IF 1725 END DO 1726* 1727 NTEST = 00 1728 IF(NTEST.GE.100) THEN 1729 WRITE(6,*) ' Splitting EI operator : ' 1730 WRITE(6,*) ' Input EI and output E, I ' 1731 CALL WRT_SPOX_TP(I_EIOP,1) 1732 CALL WRT_SPOX_TP(I_EOP,1) 1733 CALL WRT_SPOX_TP(I_IOP,1) 1734 END IF 1735* 1736 RETURN 1737 END 1738 SUBROUTINE DIMENSION_EI_EXP(ISPOBEX_TP,IB_EXTP,IB_INTP,N_EXTP, 1739 & N_INTP,N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX, 1740 / NDIM_EX_ST,NDIM_IN_ST,NDIM_EI,NDIM_IN_SE,NSYM) 1741* 1742* NOTE: Symmetry of EI is assumed to be 1 (totsym) 1743* 1744* Obtain dimension for EI expansion 1745* ISPOBEX_TP contains all spinorbital excitations with external types 1746* starting at IB_EXTP and internal types starting at IB_INTP. 1747* The combinations of internal and external types are specified 1748* by N_IN_FOR_EX, IB_IN_FOR_EX, I_IN_FOR_EX. 1749* 1750*. On output 1751* NDIM_IN_ST(ISYM,INTP) : Number of internal strings of sym ISYM and type INTP 1752* NDIM_EX_ST(ISYM,INTP) : Number of external strings of sym ISYM and type INTP 1753* NDIM_IN_SE(ISYM,IEXTP) : Number of internal strings per symmetry and extenal type 1754* N_DIM_EI : Dimension of complete expansion 1755* 1756*. Jeppe Olsen, October 2006 1757* 1758 INCLUDE 'wrkspc.inc' 1759 INCLUDE 'cgas.inc' 1760 INCLUDE 'cprnt.inc' 1761*. Input 1762 INTEGER ISPOBEX_TP(NGAS,4,*) 1763 INTEGER N_IN_FOR_EX(N_EXTP),IB_IN_FOR_EX(N_EXTP), 1764 & I_IN_FOR_EX(N_EXTP) 1765*. Local scratch 1766 INTEGER ISCR1(MXPSTT),ISCR2(MXPSTT) 1767*. Output 1768 INTEGER NDIM_EX_ST(NSYM,N_EXTP), NDIM_IN_ST(NSYM,N_INTP) 1769 INTEGER NDIM_IN_SE(NSYM,N_EXTP) 1770* 1771 NTEST = 100 1772 NTEST = MAX(NTEST,IPRCC) 1773* 1774*. External parts 1775 MX_ST_TSOSO_EX = 0 1776 MX_ST_TSOSO_BLK_EX = 0 1777 MX_TBLK_EX = 0 1778 MX_TBLK_AS_EX = 0 1779* 1780 DO ISYM = 1, NSYM 1781C IDIM_TCC(ITSOSO_TP,NTSOSO_TP,ISYM, 1782C & MX_ST_TSOSO,MX_ST_TSOSO_BLK,MX_TBLK, 1783C & LTSOSO_TP,IBTSOSO_TP,IDIM_T,MSCOMB_CC, 1784C & MX_TBLK_AS,ISPOX_FOR_OCCLS,NOCCLS,IBSPOX_FOR_OCCLS, 1785C & NTCONF,IPRCC) 1786*. External parts 1787 CALL IDIM_TCC(ISPOBEX_TP(1,1,IB_EXTP),N_EXTP,ISYM, 1788 & MX_ST_TSOSO_EXL,MX_ST_TSOSO_BLK_EXL,MX_TBLK_EXL, 1789 & ISCR1, ISCR2, IDIM_T_EXL, 0, 1790 & MX_TBLK_AS_EXL,IDUM,0,IDIM,IDUM,IPRCC) 1791 MX_ST_TSOSO_EX = MAX(MX_ST_TSOSO_EX,MX_ST_TSOSO_EXL) 1792 DO I_EXTP = 1, N_EXTP 1793 NDIM_EX_ST(ISYM,I_EXTP) = ISCR1(I_EXTP) 1794 END DO 1795 MX_ST_TSOSO_BLK_EX = MAX(MX_ST_TSOSO_BLK_EX,MX_ST_TSOSO_BLK_EXL) 1796 MX_TBLK_EX = MAX(MX_TBLK_EX,MX_TBLK_EXL) 1797 MX_TBLK_AS_EX = MAX(MX_TBLK_AS_EX,MX_TBLK_AS_EXL) 1798 END DO 1799*. Internal parts 1800 MX_ST_TSOSO_IN = 0 1801 MX_ST_TSOSO_BLK_IN = 0 1802 MX_TBLK_IN = 0 1803 MX_TBLK_AS_IN = 0 1804* 1805 DO ISYM = 1, NSYM 1806 CALL IDIM_TCC(ISPOBEX_TP(1,1,IB_INTP),N_INTP,ISYM, 1807 & MX_ST_TSOSO_INL,MX_ST_TSOSO_BLK_INL,MX_TBLK_INL, 1808 & ISCR1, ISCR2, IDIM_T_INL, 0, 1809 & MX_TBLK_AS_INL,IDUM,0,IDIM,IDUM,IPRCC) 1810 MX_ST_TSOSO_IN = MAX(MX_ST_TSOSO_IN,MX_ST_TSOSO_INL) 1811 DO I_INTP = 1, N_INTP 1812 NDIM_IN_ST(ISYM,I_INTP) = ISCR1(I_INTP) 1813 END DO 1814 MX_ST_TSOSO_BLK_IN = MAX(MX_ST_TSOSO_BLK_IN,MX_ST_TSOSO_BLK_INL) 1815 MX_TBLK_IN = MAX(MX_TBLK_IN,MX_TBLK_INL) 1816 MX_TBLK_AS_IN = MAX(MX_TBLK_AS_IN,MX_TBLK_AS_INL) 1817 END DO 1818*. Obtain number of internals per symmetry and external. 1819 DO I_EXTP = 1, N_EXTP 1820 IB = IB_IN_FOR_EX(I_EXTP) 1821 N_IN = N_IN_FOR_EX(I_EXTP) 1822 DO ISYM = 1, NSYM 1823 NDIM_IN_SE(ISYM,I_EXTP) = 0 1824 DO II_INTP = 1, N_IN 1825 I_INTP = I_IN_FOR_EX(IB-1+II_INTP) 1826 NDIM_IN_SE(ISYM,I_EXTP) = 1827 & NDIM_IN_SE(ISYM,I_EXTP) + NDIM_IN_ST(ISYM,I_INTP) 1828 END DO 1829 END DO 1830 END DO 1831* 1832 IF(NTEST.GE.100) THEN 1833 WRITE(6,*) ' Dimension of externals: ' 1834 WRITE(6,*) ' ======================== ' 1835 DO I_EXTP = 1, N_EXTP 1836 WRITE(6,*) ' External type : ', I_EXTP 1837 CALL IWRTMA(NDIM_EX_ST(1,I_EXTP),1,NSYM,1,NSYM) 1838 END DO 1839 WRITE(6,*) ' Dimension of internals: ' 1840 WRITE(6,*) ' ======================== ' 1841 DO I_INTP = 1, N_INTP 1842 WRITE(6,*) ' Internal type : ', I_INTP 1843 CALL IWRTMA(NDIM_IN_ST(1,I_INTP),1,NSYM,1,NSYM) 1844 END DO 1845 WRITE(6,*) ' Number of internals per symmetry and external:' 1846 WRITE(6,*) ' ============================================== ' 1847 CALL IWRTMA(NDIM_IN_SE,NSYM,N_EXTP,NSYM,N_EXTP) 1848 END IF 1849*. And then for the various combinations of internals and externals for sym1 1850*. (for comparison with standard test) 1851 NDIM_EI = 0 1852 DO I_EXTP = 1, N_EXTP 1853 DO II_INTP = 1, N_IN_FOR_EX(I_EXTP) 1854 I_INTP = I_IN_FOR_EX(IB_IN_FOR_EX(I_EXTP)-1+II_INTP) 1855 DO ISYM = 1, NSYM 1856 NDIM_EI = NDIM_EI 1857 & + NDIM_EX_ST(ISYM,I_EXTP)*NDIM_IN_ST(ISYM,I_INTP) 1858 END DO 1859 END DO 1860 END DO 1861* 1862 IF(NTEST.GE.100) THEN 1863 WRITE(6,*) ' Number of EI combinations with sym 1 = ', NDIM_EI 1864 END IF 1865* 1866 RETURN 1867 END 1868 SUBROUTINE EI_REORDER_ARRAYS(IREO_EI_ST, 1869 & ISPOBEX_TP,IB_EXTP,IB_INTP, 1870 & L_SPOBEX_TP,IB_SPOBEX_TP,NSPOBEX_TPL, 1871 & N_EXTP,N_INTP, 1872 & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX,ITSYM) 1873* 1874* Generate reorder array for going between EI and standard CAAB type 1875* orders 1876* Jeppe Olsen, October 2006, modified March 2009: Changed into 1877* T(Internal, external ordering), with all internal and external 1878* strings corresponding to given type of external string and 1879* symmetry of internal string(sic) are a single matrix block 1880*. Total symmetry of operator is ITSYM 1881* 1882 INCLUDE 'wrkspc.inc' 1883 INCLUDE 'cgas.inc' 1884 INCLUDE 'lucinp.inc' 1885 INCLUDE 'orbinp.inc' 1886 INCLUDE 'strinp.inc' 1887 INCLUDE 'ctcc.inc' 1888*. Some of the dimensions. 1889 INTEGER N_IN_FOR_EX(N_EXTP) 1890 1891* 1892 IDUM = 0 1893 IATP = 1 1894 IBTP = 2 1895* 1896 NAEL = NELEC(IATP) 1897 NBEL = NELEC(IBTP) 1898 NELMAX = MAX(NAEL,NBEL) 1899*. largest number of external types for given internal type 1900 MX_INTP = IMNMX(N_IN_FOR_EX,N_EXTP,2) 1901C? WRITE(6,*) ' MX_INTP = ', MX_INTP 1902 IF(MX_INTP.GT.MXP_NINTP_FOR_EX) THEN 1903 WRITE(6,*) ' Problem with fixed dimension ' 1904 WRITE(6,*) ' Observed in routine EI_REORDER_ARRAYS' 1905 WRITE(6,*) ' MXP_NINTP_FOR_EX: actual and required value', 1906 & MX_INTP,MXP_NINTP_FOR_EX 1907 WRITE(6,*) ' Increase value of MXP_NINTP_FOR_EX' 1908 STOP ' Increase value of MXP_NINTP_FOR_EX' 1909 END IF 1910* 1911 CALL MEMMAN(IDUM,IDUM,'MARK ',1,'EI_REO') 1912*. 1913 CALL MEMMAN(KL_IST_SCR,MX_ST_TSOSO_BLK_MX,'ADDL ',1,'I_STR_') 1914 LEN_Z = NELMAX*NTOOB 1915 CALL MEMMAN(KL_IZ_CA,LEN_Z*MX_INTP,'ADDL ',1,'IZ_CA ') 1916 CALL MEMMAN(KL_IZ_CB,LEN_Z*MX_INTP,'ADDL ',1,'IZ_CB ') 1917 CALL MEMMAN(KL_IZ_AA,LEN_Z*MX_INTP,'ADDL ',1,'IZ_AA ') 1918 CALL MEMMAN(KL_IZ_AB,LEN_Z*MX_INTP,'ADDL ',1,'IZ_AB ') 1919 L_IZ_SCR = (NELMAX+1)*(NTOOB+1) + 2 * NTOOB 1920 CALL MEMMAN(KL_IZ_SCR,L_IZ_SCR,'ADDL ',1,'IZ_SCR') 1921* 1922 LEN_REO= NSMOB*MX_ST_TSOSO_MX 1923 CALL MEMMAN(KL_IREO_CA,MX_INTP*LEN_REO,'ADDL ',1,'IRE_CA') 1924 CALL MEMMAN(KL_IREO_CB,MX_INTP*LEN_REO,'ADDL ',1,'IRE_CB') 1925 CALL MEMMAN(KL_IREO_AA,MX_INTP*LEN_REO,'ADDL ',1,'IRE_AA') 1926 CALL MEMMAN(KL_IREO_AB,MX_INTP*LEN_REO,'ADDL ',1,'IRE_AB') 1927*. 1928 LEN_STBLK = NELMAX*MX_ST_TSOSO_BLK_MX 1929 CALL MEMMAN(KL_IST_EX_CA,LEN_STBLK*NSMOB,'ADDL ',1,'STE_CA') 1930 CALL MEMMAN(KL_IST_EX_CB,LEN_STBLK*NSMOB,'ADDL ',1,'STE_CB') 1931 CALL MEMMAN(KL_IST_EX_AA,LEN_STBLK*NSMOB,'ADDL ',1,'STE_AA') 1932 CALL MEMMAN(KL_IST_EX_AB,LEN_STBLK*NSMOB,'ADDL ',1,'STE_AB') 1933* 1934 CALL MEMMAN(KL_IST_IN_CA,LEN_STBLK*NSMOB*MX_INTP,'ADDL ',1, 1935 & 'STI_CA') 1936 CALL MEMMAN(KL_IST_IN_CB,LEN_STBLK*NSMOB*MX_INTP,'ADDL ',1, 1937 & 'STI_CB') 1938 CALL MEMMAN(KL_IST_IN_AA,LEN_STBLK*NSMOB*MX_INTP,'ADDL ',1, 1939 & 'STI_AA') 1940 CALL MEMMAN(KL_IST_IN_AB,LEN_STBLK*NSMOB*MX_INTP,'ADDL ',1, 1941 & 'STI_AB') 1942*. Obtain array for reordering types from EI to standard 1943C? WRITE(6,*) ' IB_EXTP, IB_INTP =', IB_EXTP, IB_INTP 1944 CALL MEMMAN(KL_EI_TP_REO,NSPOBEX_TP+1,'ADDL ',1,'EI_TPR') 1945 CALL EITP_TO_SPOXTP_REO(WORK(KL_EI_TP_REO),ISPOBEX_TP, 1946 & NSPOBEX_TPL,IB_EXTP,IB_INTP,NGAS,N_EXTP,N_INTP, 1947 / N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX) 1948C EITP_TO_SPOXTP_REO(I_EI_TP_REO,ISPOBEX_TP,NSPOBEX_TP, 1949C & IB_EXTP, IB_INTP,NGAS,N_EXTP,N_INTP, 1950C & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX) 1951 CALL MEMCHK2('AFT_TP') 1952* 1953 LCHECK = 123456789 1954 CALL EI_REORDER_ARRAYS_S(IREO_EI_ST, 1955 & ISPOBEX_TP,IB_EXTP,IB_INTP,NGAS, 1956 & L_SPOBEX_TP,IB_SPOBEX_TP,NOBPT,NTOOB, 1957 & NSPOBEX_TPL,N_EXTP,N_INTP, 1958 & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX, 1959 & WORK(KL_EI_TP_REO),ITSYM,NSMOB, 1960 & WORK(KL_IST_SCR),WORK(KL_IZ_CA),WORK(KL_IZ_CB), 1961 / WORK(KL_IZ_AA),WORK(KL_IZ_AB),WORK(KL_IZ_SCR), 1962 / WORK(KL_IREO_CA),WORK(KL_IREO_CB), 1963 & WORK(KL_IREO_AA),WORK(KL_IREO_AB), 1964 & WORK(KL_IST_EX_CA),WORK(KL_IST_EX_CB), 1965 & WORK(KL_IST_EX_AA),WORK(KL_IST_EX_AB), 1966 & WORK(KL_IST_IN_CA),WORK(KL_IST_IN_CB), 1967 & WORK(KL_IST_IN_AA),WORK(KL_IST_IN_AB),LEN_Z,LEN_REO, 1968 & LEN_STBLK,LCHECK) 1969 CALL MEMMAN(IDUM,IDUM,'FLUSM ',1,'EI_REO') 1970 RETURN 1971 END 1972 SUBROUTINE EI_REORDER_ARRAYS_S(IREO_EI_ST, 1973 & ISPOBEX_TP,IB_EXTP, IB_INTP,NGAS, 1974 & L_SPOBEX_TP,IB_SPOBEX_TP,NOBPT,NTOOB, 1975 & NSPOBEX_TP,N_EXTP,N_INTP, 1976 & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX, 1977 & I_EI_TP_REO,ITSYM,NSMOB, 1978 & IST_SCR,IZ_CA,IZ_CB,IZ_AA,IZ_AB,IZ_SCR, 1979 & IREO_CA,IREO_CB,IREO_AA,IREO_AB, 1980 & IST_EX_CA,IST_EX_CB,IST_EX_AA,IST_EX_AB, 1981 & IST_IN_CA,IST_IN_CB,IST_IN_AA,IST_IN_AB,LEN_Z,LEN_REO, 1982 & LEN_STBLK,LCHECK) 1983*. Obtain ordering : EI-ordering => standard-ordering 1984* 1985*. Jeppe Olsen, October 2006, modified march 2009 1986* 1987* Obtain reorder array going from standard CAAB ordering into 1988* EI ordering. The EI order is T(internal, external) 1989* 1990 INCLUDE 'wrkspc.inc' 1991 INCLUDE 'multd2h.inc' 1992*. Input 1993 INTEGER ISPOBEX_TP(NGAS,4,*) 1994 INTEGER NOBPT(NGAS) 1995*. Length and offset for each spinorbitalexcitation 1996 INTEGER L_SPOBEX_TP(NSPOBEX_TP),IB_SPOBEX_TP(NSPOBEX_TP) 1997 INTEGER N_IN_FOR_EX(N_EXTP),IB_IN_FOR_EX(N_EXTP),I_IN_FOR_EX(*) 1998 INTEGER I_EI_TP_REO(NSPOBEX_TP) 1999*. Output : reorder array with signs included 2000 INTEGER IREO_EI_ST(*) 2001*. Scratch through argument list 2002*. Space for occupations of strings 2003 INTEGER IST_SCR(*) 2004*. Z arrays for lexical ordering and scratch array for constructing then 2005 INTEGER IZ_CA(LEN_Z,*),IZ_CB(LEN_Z,*), 2006 & IZ_AA(LEN_Z,*),IZ_AB(LEN_Z,*),IZ_SCR(*) 2007 2008*. Reorder arrays 2009 INTEGER IREO_CA(LEN_REO,*), IREO_CB(LEN_REO,*), 2010 & IREO_AA(LEN_REO,*), IREO_AB(LEN_REO,*) 2011*. Occupation of external and internal strings 2012 INTEGER IST_EX_CA(LEN_STBLK,NSMOB), IST_EX_CB(LEN_STBLK,NSMOB) 2013 INTEGER IST_EX_AA(LEN_STBLK,NSMOB), IST_EX_AB(LEN_STBLK,NSMOB) 2014 INTEGER IST_IN_CA(LEN_STBLK,NSMOB,*), IST_IN_CB(LEN_STBLK,NSMOB,*) 2015 INTEGER IST_IN_AA(LEN_STBLK,NSMOB,*), IST_IN_AB(LEN_STBLK,NSMOB,*) 2016*. Local scratch ( DIMENSIONS ???? ) 2017 INTEGER I_EI_CA(MXPORB),I_EI_CB(MXPORB),I_EI_AA(MXPORB), 2018 & I_EI_AB(MXPORB) 2019 INTEGER IGRP(MXPNGAS) 2020* 2021 INTEGER NEL_IN_CA(MXP_NINTP_FOR_EX),NEL_IN_CB(MXP_NINTP_FOR_EX) 2022 INTEGER NEL_IN_AA(MXP_NINTP_FOR_EX),NEL_IN_AB(MXP_NINTP_FOR_EX) 2023 INTEGER NEL_CA(MXP_NINTP_FOR_EX),NEL_CB(MXP_NINTP_FOR_EX) 2024 INTEGER NEL_AA(MXP_NINTP_FOR_EX),NEL_AB(MXP_NINTP_FOR_EX) 2025 INTEGER I_ST_TP(MXP_NINTP_FOR_EX) 2026 INTEGER IB_ST_TP(MXP_NINTP_FOR_EX) 2027 INTEGER LEN_ST(MXP_NINTP_FOR_EX) 2028 INTEGER IB_ST_SSS(8,8,8,MXP_NINTP_FOR_EX) 2029* 2030 INTEGER NST_CA(8,MXP_NINTP_FOR_EX), 2031 & NST_CB(8,MXP_NINTP_FOR_EX), 2032 & NST_AA(8,MXP_NINTP_FOR_EX), 2033 & NST_AB(8,MXP_NINTP_FOR_EX) 2034* 2035 INTEGER NST_EX_CA(8),NST_EX_CB(8),NST_EX_AA(8),NST_EX_AB(8) 2036* 2037 INTEGER NST_IN_CA(8,MXP_NINTP_FOR_EX), 2038 & NST_IN_CB(8,MXP_NINTP_FOR_EX), 2039 & NST_IN_AA(8,MXP_NINTP_FOR_EX), 2040 & NST_IN_AB(8,MXP_NINTP_FOR_EX) 2041 2042* 2043 NTEST = 10 2044 IF(NTEST.GE.10) THEN 2045 WRITE(6,*) 2046 WRITE(6,*) ' ------------------------------------------------' 2047 WRITE(6,*) ' Information from subroutine EI_REORDER_ARRAYS_S ' 2048 WRITE(6,*) ' ------------------------------------------------' 2049 WRITE(6,*) 2050 END IF 2051C? WRITE(6,*) 'LEN_STBLK = ', LEN_STBLK 2052C? WRITE(6,*) ' I_EI_TP_REO at start of EI_REORDER' 2053C? CALL IWRTMA(I_EI_TP_REO,1,NSPOBEX_TP,1,NSPOBEX_TP) 2054* 2055 I_EI_TP = 0 2056 I_EI_OP = 0 2057 DO J_EXTP = 1, N_EXTP 2058 J_EXTP_ABS = IB_EXTP-1+J_EXTP 2059* 2060 NEL_EX_CA = IELSUM(ISPOBEX_TP(1,1,J_EXTP_ABS),NGAS) 2061 NEL_EX_CB = IELSUM(ISPOBEX_TP(1,2,J_EXTP_ABS),NGAS) 2062 NEL_EX_AA = IELSUM(ISPOBEX_TP(1,3,J_EXTP_ABS),NGAS) 2063 NEL_EX_AB = IELSUM(ISPOBEX_TP(1,4,J_EXTP_ABS),NGAS) 2064* 2065 N_INTP_J = N_IN_FOR_EX(J_EXTP) 2066 IF(NTEST.GE.100) THEN 2067 WRITE(6,*) 2068 WRITE(6,*) ' Information for external type = ', J_EXTP 2069 WRITE(6,*) ' Number of internal types =', N_INTP_J 2070 WRITE(6,*) ' I_EI_OP at start ', I_EI_OP 2071 WRITE(6,*) 2072 END IF 2073 IB_INTP_J = IB_IN_FOR_EX(J_EXTP) 2074* 2075*. Construct information for the various internal and combined ei types 2076*. that will be used for this external type 2077 DO JJ_INTP = 1, N_INTP_J 2078 I_EI_TP = I_EI_TP + 1 2079 J_INTP = I_IN_FOR_EX(IB_INTP_J-1+JJ_INTP) 2080 J_INTP_ABS = IB_INTP-1+J_INTP 2081*. Corresponding type in standard order 2082 I_ST_TP(JJ_INTP) = I_EI_TP_REO(I_EI_TP) 2083 IF(NTEST.GE.10000) 2084 & WRITE(6,*) ' JJ_INTP, I_EI_TP, I_ST, I_EI ', 2085 & JJ_INTP,I_EI_TP,I_EI_TP_REO(I_EI_TP),I_ST_TP(JJ_INTP) 2086 IB_ST_TP(JJ_INTP) = IB_SPOBEX_TP(I_ST_TP(JJ_INTP)) 2087C? WRITE(6,*) ' IB_ST_TP = ', IB_ST_TP 2088* 2089 IF(NTEST.GE.10000) THEN 2090 WRITE(6,'(A,4I4)') ' J_INTP, IB_INTP, J_INTP_ABS =', 2091 & J_INTP, IB_INTP, J_INTP_ABS 2092 END IF 2093*. Number of operators in internal part 2094 NEL_IN_CA(JJ_INTP) = IELSUM(ISPOBEX_TP(1,1,J_INTP_ABS),NGAS) 2095 NEL_IN_CB(JJ_INTP) = IELSUM(ISPOBEX_TP(1,2,J_INTP_ABS),NGAS) 2096 NEL_IN_AA(JJ_INTP) = IELSUM(ISPOBEX_TP(1,3,J_INTP_ABS),NGAS) 2097 NEL_IN_AB(JJ_INTP) = IELSUM(ISPOBEX_TP(1,4,J_INTP_ABS),NGAS) 2098*. Number of operators in combined external/internal part 2099 NEL_CA(JJ_INTP) =IELSUM(ISPOBEX_TP(1,1,I_ST_TP(JJ_INTP)),NGAS) 2100 NEL_CB(JJ_INTP) =IELSUM(ISPOBEX_TP(1,2,I_ST_TP(JJ_INTP)),NGAS) 2101 NEL_AA(JJ_INTP) =IELSUM(ISPOBEX_TP(1,3,I_ST_TP(JJ_INTP)),NGAS) 2102 NEL_AB(JJ_INTP) =IELSUM(ISPOBEX_TP(1,4,I_ST_TP(JJ_INTP)),NGAS) 2103* 2104 IF(NTEST.GE.10000) THEN 2105 WRITE(6,'(A,4I4)') ' << E, I, EI, ST types >> ', 2106 & J_EXTP, J_INTP, I_EI_TP, I_ST_TP(JJ_INTP) 2107 END IF 2108 IF(NTEST.GE.10000) THEN 2109 WRITE(6,'(A,4I4)') 'NEL_CA, NEL_CB, NEL_AA, NEL_AB = ', 2110 & NEL_CA(JJ_INTP), NEL_CB(JJ_INTP), 2111 & NEL_AA(JJ_INTP), NEL_AB(JJ_INTP) 2112 END IF 2113*. Construct reorder arrays for the combined CA, CB, AA, AB strings 2114C WEIGHT_SPGP(Z,NORBTP,NELFTP,NORBFTP,ISCR,NTEST) 2115 CALL WEIGHT_SPGP(IZ_CA(1,JJ_INTP),NGAS, 2116 & ISPOBEX_TP(1,1,I_ST_TP(JJ_INTP)),NOBPT,IZ_SCR,0) 2117 CALL WEIGHT_SPGP(IZ_CB(1,JJ_INTP),NGAS, 2118 & ISPOBEX_TP(1,2,I_ST_TP(JJ_INTP)),NOBPT,IZ_SCR,0) 2119 CALL WEIGHT_SPGP(IZ_AA(1,JJ_INTP),NGAS, 2120 & ISPOBEX_TP(1,3,I_ST_TP(JJ_INTP)),NOBPT,IZ_SCR,0) 2121 CALL WEIGHT_SPGP(IZ_AB(1,JJ_INTP),NGAS, 2122 & ISPOBEX_TP(1,4,I_ST_TP(JJ_INTP)),NOBPT,IZ_SCR,0) 2123 DO ISYM = 1, NSMOB 2124*. CA 2125 CALL OCC_TO_GRP(ISPOBEX_TP(1,1,I_ST_TP(JJ_INTP)),IGRP,1) 2126C GETSTR2_TOTSM_SPGP(IGRP,NIGRP,ISPGRPSM,NEL,NSTR,ISTR, 2127C & NORBT,IDOREO,IZ,IREO) 2128 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS,ISYM,NEL_CA(JJ_INTP), 2129 & NST_CA_L,IST_SCR,NTOOB,1,IZ_CA(1,JJ_INTP), 2130 & IREO_CA(1,JJ_INTP)) 2131 NST_CA(ISYM,JJ_INTP) = NST_CA_L 2132*. CB 2133 CALL OCC_TO_GRP(ISPOBEX_TP(1,2,I_ST_TP(JJ_INTP)),IGRP,1) 2134 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS,ISYM,NEL_CB(JJ_INTP), 2135 & NST_CB_L,IST_SCR,NTOOB,1,IZ_CB(1,JJ_INTP), 2136 & IREO_CB(1,JJ_INTP)) 2137 NST_CB(ISYM,JJ_INTP) = NST_CB_L 2138 2139*. AA 2140 CALL OCC_TO_GRP(ISPOBEX_TP(1,3,I_ST_TP(JJ_INTP)),IGRP,1) 2141 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS,ISYM,NEL_AA(JJ_INTP), 2142 & NST_AA_L,IST_SCR,NTOOB,1,IZ_AA(1,JJ_INTP), 2143 & IREO_AA(1,JJ_INTP)) 2144 NST_AA(ISYM,JJ_INTP) = NST_AA_L 2145*. AB 2146 CALL OCC_TO_GRP(ISPOBEX_TP(1,4,I_ST_TP(JJ_INTP)),IGRP,1) 2147 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS,ISYM,NEL_AB(JJ_INTP), 2148 & NST_AB_L,IST_SCR,NTOOB,1,IZ_AB(1,JJ_INTP), 2149 & IREO_AB(1,JJ_INTP)) 2150 NST_AB(ISYM,JJ_INTP) = NST_AB_L 2151* 2152*. And the internal strings 2153*. CA 2154 CALL OCC_TO_GRP(ISPOBEX_TP(1,1,J_INTP_ABS),IGRP,1) 2155 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2156 & ISYM,NEL_IN_CA(JJ_INTP),NST_IN_CA(ISYM,JJ_INTP), 2157 & IST_IN_CA(1,ISYM,JJ_INTP),NTOOB,0, 2158 & IDUM,IDUM) 2159*. CB 2160 CALL OCC_TO_GRP(ISPOBEX_TP(1,2,J_INTP_ABS),IGRP,1) 2161 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2162 & ISYM,NEL_IN_CB(JJ_INTP),NST_IN_CB(ISYM,JJ_INTP), 2163 & IST_IN_CB(1,ISYM,JJ_INTP),NTOOB,0, 2164 & IDUM,IDUM) 2165*. AA 2166 CALL OCC_TO_GRP(ISPOBEX_TP(1,3,J_INTP_ABS),IGRP,1) 2167 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2168 & ISYM,NEL_IN_AA(JJ_INTP),NST_IN_AA(ISYM,JJ_INTP), 2169 & IST_IN_AA(1,ISYM,JJ_INTP),NTOOB,0, 2170 & IDUM,IDUM) 2171C? WRITE(6,*) ' The internal AA strings right after delivery' 2172C? CALL IWRTMA(IST_IN_AA(1,ISYM,JJ_INTP),NEL_IN_AA(JJ_INTP), 2173C? & NST_IN_AA(1,JJ_INTP) ,NEL_IN_AA(JJ_INTP), 2174C? & NST_IN_AA(1,J_INTPJ)) 2175*. AB 2176 CALL OCC_TO_GRP(ISPOBEX_TP(1,4,J_INTP_ABS),IGRP,1) 2177 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2178 & ISYM,NEL_IN_AB(JJ_INTP),NST_IN_AB(ISYM,JJ_INTP), 2179 & IST_IN_AB(1,ISYM,JJ_INTP),NTOOB,0, 2180 & IDUM,IDUM) 2181* 2182 2183 END DO 2184* ^ End of loops over symmetries 2185*. Array for accessing CA,CB,AA,AB blocks of given sym 2186 IF(NTEST.GE.10000) THEN 2187 WRITE(6,*) ' NST_CA, NST_CB, NST_AA, NST_AB:' 2188 CALL IWRTMA(NST_CA(1,JJ_INTP),1,NSMOB,1,NSMOB) 2189 CALL IWRTMA(NST_CB(1,JJ_INTP),1,NSMOB,1,NSMOB) 2190 CALL IWRTMA(NST_AA(1,JJ_INTP),1,NSMOB,1,NSMOB) 2191 CALL IWRTMA(NST_AB(1,JJ_INTP),1,NSMOB,1,NSMOB) 2192 END IF 2193* 2194 CALL Z_TCC_OFF2(IB_ST_SSS(1,1,1,JJ_INTP),LEN_ST(JJ_INTP), 2195 & NST_CA(1,JJ_INTP),NST_CB(1,JJ_INTP), 2196 & NST_AA(1,JJ_INTP),NST_AB(1,JJ_INTP),ITSYM,NSMOB) 2197 END DO 2198* ^ End of loop over internal types for given external type 2199C? WRITE(6,*) ' The NEL_IN_AA array in a test' 2200C? CALL IWRTMA(NEL_IN_AA,1, N_INTP_J,1, N_INTP_J) 2201C? WRITE(6,*) ' The NST_IN_AA-array in a test' 2202C? CALL IWRTMA(NST_IN_AA,1, N_INTP_J,8, N_INTP_J) 2203C? WRITE(6,*) ' The IN_AA strings in sym 1 ' 2204C? DO JX = 1, N_INTP_J 2205C? CALL IWRTMA(IST_IN_AA(1,1,JX),NEL_IN_AA(JX), 2206C? & NST_IN_AA(1,JX) ,NEL_IN_AA(JX), 2207C? & NST_IN_AA(1,JX)) 2208C? END DO 2209*. Occupation of external strings 2210*. CA 2211 DO ISYM = 1, NSMOB 2212 CALL OCC_TO_GRP(ISPOBEX_TP(1,1,J_EXTP_ABS),IGRP,1) 2213 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2214 & ISYM,NEL_EX_CA,NST_EX_CA(ISYM),IST_EX_CA(1,ISYM), 2215 & NTOOB,0,IDUM,IDUM) 2216*. CB 2217 CALL OCC_TO_GRP(ISPOBEX_TP(1,2,J_EXTP_ABS),IGRP,1) 2218 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2219 & ISYM,NEL_EX_CB,NST_EX_CB(ISYM),IST_EX_CB(1,ISYM), 2220 & NTOOB,0,IDUM,IDUM) 2221*. AA 2222 CALL OCC_TO_GRP(ISPOBEX_TP(1,3,J_EXTP_ABS),IGRP,1) 2223 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2224 & ISYM,NEL_EX_AA,NST_EX_AA(ISYM),IST_EX_AA(1,ISYM), 2225 & NTOOB,0,IDUM,IDUM) 2226*. AB 2227 CALL OCC_TO_GRP(ISPOBEX_TP(1,4,J_EXTP_ABS),IGRP,1) 2228 CALL GETSTR2_TOTSM_SPGP(IGRP,NGAS, 2229 & ISYM,NEL_EX_AB,NST_EX_AB(ISYM),IST_EX_AB(1,ISYM), 2230 & NTOOB,0,IDUM,IDUM) 2231 END DO 2232C? WRITE(6,*) ' End of external ' 2233C? WRITE(6,*) ' End of combined ' 2234*. Loop over the individual strings of this External type 2235* and obtain strings in EI-order 2236*. Loop over strings in EI order 2237 DO I_EX_SM = 1, NSMOB 2238 I_IN_SM = MULTD2H(I_EX_SM,ITSYM) 2239*. Loop over symmetries of external operators 2240 DO I_EX_C_SM = 1, NSMOB 2241 I_EX_A_SM = MULTD2H(I_EX_C_SM,I_EX_SM) 2242 DO I_EX_CA_SM = 1, NSMOB 2243 I_EX_CB_SM = MULTD2H(I_EX_CA_SM,I_EX_C_SM) 2244 DO I_EX_AA_SM = 1, NSMOB 2245 I_EX_AB_SM = MULTD2H(I_EX_AA_SM,I_EX_A_SM) 2246*. Loop over External strings (CA, CB, AA, AB) 2247 DO I_EX_AB = 1, NST_EX_AB(I_EX_AB_SM) 2248 DO I_EX_AA = 1, NST_EX_AA(I_EX_AA_SM) 2249 DO I_EX_CB = 1, NST_EX_CB(I_EX_CB_SM) 2250 DO I_EX_CA = 1, NST_EX_CA(I_EX_CA_SM) 2251 IF(NTEST.GE.10000) THEN 2252 WRITE(6,*) ' I_EX_AB, I_EX_AA, I_EX_CB, I_EX_CA = ', 2253 & I_EX_AB, I_EX_AA, I_EX_CB, I_EX_CA 2254 END IF 2255*. Loop over the internal types for this external types 2256 DO JJ_INTP = 1, N_INTP_J 2257 J_INTP = I_IN_FOR_EX(IB_INTP_J-1+JJ_INTP) 2258 J_INTP_ABS = IB_INTP-1+J_INTP 2259*. Loop over symmetries of internal operators 2260 DO I_IN_C_SM = 1, NSMOB 2261 I_IN_A_SM = MULTD2H(I_IN_C_SM,I_IN_SM) 2262 DO I_IN_CA_SM = 1, NSMOB 2263 I_IN_CB_SM = MULTD2H(I_IN_C_SM,I_IN_CA_SM) 2264 DO I_IN_AA_SM = 1, NSMOB 2265 I_IN_AB_SM = MULTD2H(I_IN_A_SM,I_IN_AA_SM) 2266 IF(NTEST.GE.10000) THEN 2267 WRITE(6,'(A,4I4)') 2268 & 'I_IN_A_SM, I_IN_CB_SM,I_IN_AB_SM,I_IN_AB_SM = ', 2269 & I_IN_A_SM, I_IN_CB_SM,I_IN_AB_SM,I_IN_AB_SM 2270 END IF 2271C? WRITE(6,*) ' I_IN_C_SM, I_IN_A_SM = ', 2272C? & I_IN_C_SM, I_IN_A_SM 2273 IF(NTEST.GE.10000) THEN 2274 WRITE(6,'(A,4I4)') 2275 & 'I_IN_CA_SM, I_IN_CB_SM, I_IN_AA_SM, I_IN_AB_SM = ', 2276 & I_IN_CA_SM, I_IN_CB_SM, I_IN_AA_SM, I_IN_AB_SM 2277 WRITE(6,'(A,4I4)') 2278 & 'I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM = ', 2279 & I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM 2280 WRITE(6,'(A,4I4)') 2281 & 'I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM = ', 2282 & I_EX_CA_SM, I_EX_CB_SM, I_EX_AA_SM, I_EX_AB_SM 2283 END IF 2284*. Symmetry of complete CA,CB,AA and AB strings 2285*. Symmetry of complete CA,CB,AA and AB strings 2286 I_CA_SM = MULTD2H(I_EX_CA_SM,I_IN_CA_SM) 2287 I_CB_SM = MULTD2H(I_EX_CB_SM,I_IN_CB_SM) 2288 I_AA_SM = MULTD2H(I_EX_AA_SM,I_IN_AA_SM) 2289 I_AB_SM = MULTD2H(I_EX_AB_SM,I_IN_AB_SM) 2290* Sign to bring ECA ECB EAA EAB ICA ICB IAA IAB into 2291* ECA ICA ECB ICB EAA IAB EAB IAB 2292 NPERMG = NEL_EX_AB* 2293 & (NEL_IN_CA(JJ_INTP)+NEL_IN_CB(JJ_INTP)+NEL_IN_AA(JJ_INTP)) 2294 & + NEL_EX_AA* 2295 & (NEL_IN_CA(JJ_INTP)+NEL_IN_CB(JJ_INTP)) 2296 & + NEL_EX_CB* 2297 & NEL_IN_CA(JJ_INTP) 2298 IF(MOD(NPERMG,2).EQ.0) THEN 2299 ISIGNG = 1 2300 ELSE 2301 ISIGNG = -1 2302 END IF 2303* 2304 IF(NTEST.GE.10000) THEN 2305 WRITE(6,'(A,4I4)') 2306 & 'I_CA_SM, I_CB_SM, I_AA_SM, I_AB_SM = ', 2307 & I_CA_SM, I_CB_SM, I_AA_SM, I_AB_SM 2308 END IF 2309*. And internal strings 2310 DO I_IN_AB = 1, NST_IN_AB(I_IN_AB_SM,JJ_INTP) 2311 DO I_IN_AA = 1, NST_IN_AA(I_IN_AA_SM,JJ_INTP) 2312 DO I_IN_CB = 1, NST_IN_CB(I_IN_CB_SM,JJ_INTP) 2313 DO I_IN_CA = 1, NST_IN_CA(I_IN_CA_SM,JJ_INTP) 2314 IF(NTEST.GE.10000) THEN 2315 WRITE(6,*) ' I_IN_AB, I_IN_AA, I_IN_CB, I_IN_CA = ', 2316 & I_IN_AB, I_IN_AA, I_IN_CB, I_IN_CA 2317 END IF 2318 I_EI_OP = I_EI_OP + 1 2319*. Merge strings to obtain complete CA,CB,AA,AB strings 2320*. CA 2321 IOFF_IN_CA = 1+(I_IN_CA-1)*NEL_IN_CA(JJ_INTP) 2322 IOFF_EX_CA = 1+(I_EX_CA-1)*NEL_EX_CA 2323 CALL PROD_TWO_STRINGS(I_EI_CA, 2324 & IST_EX_CA(IOFF_EX_CA,I_EX_CA_SM), 2325 & IST_IN_CA(IOFF_IN_CA,I_IN_CA_SM,JJ_INTP), 2326 & NEL_EX_CA,NEL_IN_CA(JJ_INTP),IS_CA) 2327*. CB 2328 IOFF_IN_CB = 1+(I_IN_CB-1)*NEL_IN_CB(JJ_INTP) 2329 IOFF_EX_CB = 1+(I_EX_CB-1)*NEL_EX_CB 2330 CALL PROD_TWO_STRINGS(I_EI_CB, 2331 & IST_EX_CB(IOFF_EX_CB,I_EX_CB_SM), 2332 & IST_IN_CB(IOFF_IN_CB,I_IN_CB_SM,JJ_INTP), 2333 & NEL_EX_CB,NEL_IN_CB(JJ_INTP),IS_CB) 2334*. AA 2335 IOFF_IN_AA = 1+(I_IN_AA-1)*NEL_IN_AA(JJ_INTP) 2336 IOFF_EX_AA = 1+(I_EX_AA-1)*NEL_EX_AA 2337 CALL PROD_TWO_STRINGS(I_EI_AA, 2338 & IST_EX_AA(IOFF_EX_AA,I_EX_AA_SM), 2339 & IST_IN_AA(IOFF_IN_AA,I_IN_AA_SM,JJ_INTP), 2340 & NEL_EX_AA,NEL_IN_AA(JJ_INTP),IS_AA) 2341*. AB 2342 IOFF_IN_AB = 1+(I_IN_AB-1)*NEL_IN_AB(JJ_INTP) 2343 IOFF_EX_AB = 1+(I_EX_AB-1)*NEL_EX_AB 2344 CALL PROD_TWO_STRINGS(I_EI_AB, 2345 & IST_EX_AB(IOFF_EX_AB,I_EX_AB_SM), 2346 & IST_IN_AB(IOFF_IN_AB,I_IN_AB_SM,JJ_INTP), 2347 & NEL_EX_AB,NEL_IN_AB(JJ_INTP),IS_AB) 2348*. And Adresses of combined strings 2349C ISTRNM(IOCC,NORB,NEL,Z,NEWORD,IREORD) 2350*. CA 2351 I_CA_ADR = ISTRNM(I_EI_CA,NTOOB,NEL_CA(JJ_INTP), 2352 & IZ_CA(1,JJ_INTP),IREO_CA(1,JJ_INTP),1) 2353*. CB 2354 I_CB_ADR = ISTRNM(I_EI_CB,NTOOB,NEL_CB(JJ_INTP), 2355 & IZ_CB(1,JJ_INTP),IREO_CB(1,JJ_INTP),1) 2356*. AA 2357 I_AA_ADR = ISTRNM(I_EI_AA,NTOOB,NEL_AA(JJ_INTP), 2358 & IZ_AA(1,JJ_INTP),IREO_AA(1,JJ_INTP),1) 2359*. AB 2360 I_AB_ADR = ISTRNM(I_EI_AB,NTOOB,NEL_AB(JJ_INTP), 2361 & IZ_AB(1,JJ_INTP),IREO_AB(1,JJ_INTP),1) 2362* 2363 IF(NTEST.GE.10000) THEN 2364 WRITE(6,*) ' CA, CB, AA, AB strings : ' 2365 CALL IWRTMA(I_EI_CA,1,NEL_CA(JJ_INTP),1,NEL_CA(JJ_INTP)) 2366 CALL IWRTMA(I_EI_CB,1,NEL_CB(JJ_INTP),1,NEL_CB(JJ_INTP)) 2367 CALL IWRTMA(I_EI_AA,1,NEL_AA(JJ_INTP),1,NEL_AA(JJ_INTP)) 2368 CALL IWRTMA(I_EI_AB,1,NEL_AB(JJ_INTP),1,NEL_AB(JJ_INTP)) 2369 WRITE(6,'(A,4I5)') ' Adresses of CA, CB, AA, AB', 2370 & I_CA_ADR, I_CB_ADR, I_AA_ADR, I_AB_ADR 2371 WRITE(6,*) ' Offset of symmetry-blocks in ST', 2372 & IB_ST_SSS(I_CA_SM,I_CB_SM,I_AA_SM,JJ_INTP) 2373 END IF 2374* CA, CB, AA, AB adress 2375 IST_ADR = 2376 & (I_AB_ADR-1)*NST_AA(I_AA_SM,JJ_INTP)*NST_CB(I_CB_SM,JJ_INTP) 2377 & *NST_CA(I_CA_SM,JJ_INTP) 2378 & +(I_AA_ADR-1) *NST_CB(I_CB_SM,JJ_INTP) 2379 & *NST_CA(I_CA_SM,JJ_INTP) 2380 & +(I_CB_ADR-1)*NST_CA(I_CA_SM,JJ_INTP) 2381 & + I_CA_ADR + IB_ST_TP(JJ_INTP) - 1 2382 & + IB_ST_SSS(I_CA_SM,I_CB_SM,I_AA_SM,JJ_INTP)-1 2383* 2384*. Signs 2385*. The sign is composed of two parts 2386* 1: Sign to bring ECA ECB EAA EAB ICA ICB IAA IAB into 2387* ECA ICA ECB ICB EAA IAB EAB IAB 2388* 2. Sign to bring each of ECA ICA, ECB ICB, EAA IAB, EAB IAB 2389* into standard order 2390 IST_SIGN = IS_CA*IS_CB*IS_AA*IS_AB*ISIGNG 2391 IF(NTEST.GE.10000) 2392 & WRITE(6,*) ' I_EI_OP, IST_ADR, ISIGNG; IST_SIGN', 2393 & I_EI_OP, IST_ADR, ISIGNG, IST_SIGN 2394*. And now : what we all have been waiting for a few hundred lines.. 2395 IREO_EI_ST(I_EI_OP) = IST_SIGN*IST_ADR 2396* 2397 END DO 2398 END DO 2399 END DO 2400 END DO 2401*. ^ End of loop over internal CA,CB,AA,AB strings of given sym 2402 END DO 2403 END DO 2404 END DO 2405*. ^ End of loop over symmetry of internal CA, CB, AA, AB strings 2406 END DO 2407*. ^ End of loop over internal types 2408 END DO 2409 END DO 2410 END DO 2411 END DO 2412* ^ End of loop over external CA, CB,AA,AB strings of given sym 2413 END DO 2414 END DO 2415 END DO 2416*. ^ End of loop over symmetry of external CA, CB, AA, AB strings 2417 END DO 2418*. ^ End of loop over symmetry of external operators 2419 END DO 2420*. ^ End of loop over external types 2421 N_EI_OP = I_EI_OP 2422* 2423*. Sum check 2424 I_DO_SUMCHECK = 1 2425 IF(I_DO_SUMCHECK.EQ.1.AND.N_EI_OP.LT.30000) THEN 2426*. Check that sum SUM_I IREO_EI_ST(I) = N_EI_OP*(N_EI_OP+1)/2 2427 ISUM = 0 2428 DO I = 1, N_EI_OP 2429 ISUM = ISUM + ABS(IREO_EI_ST(I)) 2430 END DO 2431 IF(ISUM.EQ.N_EI_OP*(N_EI_OP+1)/2) THEN 2432 WRITE(6,*) ' Sumcheck passed ' 2433 ELSE 2434 WRITE(6,*) ' Sumcheck failed in REO_EI..' 2435 WRITE(6,'(A,2I7)') 2436 & ' Observed and expected sum', ISUM,N_EI_OP*(N_EI_OP+1)/2 2437*. Determine elements that were not obtained exactly one time- 2438*, N-squared algorithm in use, as I do not want to allocate another 2439*. array.. 2440 DO ITARGET = 1, N_EI_OP 2441 NFOUND = 0 2442 DO I = 1,N_EI_OP 2443 IF(IREO_EI_ST(I).EQ.ITARGET) NFOUND = NFOUND + 1 2444 END DO 2445 IF(NFOUND.NE.1) THEN 2446 WRITE(6,*) 2447 & ' Element ', ITARGET ,' Found ', NFOUND ,' times' 2448 END IF 2449 END DO 2450*. Print reorder array before stop 2451 IF(NTEST.GE.1000) THEN 2452 WRITE(6,*) ' The reorder array EI => standard order ' 2453 CALL IWRTMA(IREO_EI_ST,1,N_EI_OP,1,N_EI_OP,1) 2454 END IF 2455C STOP ' sumcheck failed ' 2456 END IF 2457 END IF 2458 IF(NTEST.GE.1000) THEN 2459 WRITE(6,*) ' The reorder array EI => standard order ' 2460 CALL IWRTMA(IREO_EI_ST,1,N_EI_OP,1,N_EI_OP,1) 2461 END IF 2462* 2463 RETURN 2464 END 2465 SUBROUTINE PROD_TWO_STRINGS(ISTR_AB,ISTR_A,ISTR_B,NEL_A,NEL_B,IS) 2466* 2467* Two strings are given in standard ascending order 2468* Obtain product of string in ascending order, and sign 2469* required for conversion 2470* 2471*. Jeppe Olsen, October 2006 2472*. Input 2473 INTEGER ISTR_A(NEL_A),ISTR_B(NEL_B) 2474*. Output 2475 INTEGER ISTR_AB(NEL_A+NEL_B) 2476* 2477 IEL_A = 1 2478 IEL_B = 1 2479 NPERM = 0 2480 DO IEL_AB = 1, NEL_A+NEL_B 2481 IF(IEL_B.GT.NEL_B) THEN 2482*. No more B-electrons so next electron is from A 2483 ISTR_AB(IEL_AB) = ISTR_A(IEL_A) 2484 IEL_A = IEL_A + 1 2485 ELSE IF (IEL_A.GT.NEL_A) THEN 2486*. No more A-electrons so next electron is from B 2487 ISTR_AB(IEL_AB) = ISTR_B(IEL_B) 2488 IEL_B = IEL_B + 1 2489 ELSE 2490*. There are both A and B-electrons left so compare 2491 IF(ISTR_A(IEL_A).LE.ISTR_B(IEL_B)) THEN 2492*. Next electron is from IEL_A 2493 ISTR_AB(IEL_AB) = ISTR_A(IEL_A) 2494 IEL_A = IEL_A + 1 2495 ELSE 2496*. Next electron is from IEL_B 2497 ISTR_AB(IEL_AB) = ISTR_B(IEL_B) 2498 IEL_B = IEL_B + 1 2499*. Number of permutations over A-electrons to get it in place 2500 NPERM = NPERM + NEL_A-IEL_A+1 2501 END IF 2502 END IF 2503 END DO 2504*. Well, I am not sure that the above count of permutations is correct 2505*. so here is another try-starting with the smallest b, and checking how 2506*. many a's it should be moved passed 2507 NPERM = 0 2508 DO IEL_B = 1, NEL_B 2509 IORB_B = ISTR_B(IEL_B) 2510*. Find largest index in A that is smaller than IORB_B 2511 I_SMALL = 0 2512 DO IEL_A = 1, NEL_A 2513 IF(ISTR_A(IEL_A).LT.IORB_B) I_SMALL = IEL_A 2514 END DO 2515 NPERM = NPERM + NEL_A-I_SMALL 2516 END DO 2517* 2518 IF(MOD(NPERM,2).EQ.0) THEN 2519 IS = 1 2520 ELSE 2521 IS = -1 2522 END IF 2523* 2524 NTEST = 00 2525 IF(NTEST.GE.100) THEN 2526 NEL_AB = NEL_A + NEL_B 2527 WRITE(6,*) 'Merging two strings : AB, A, B ' 2528 WRITE(6,*) 'Number of operators in A,B', NEL_A, NEL_B 2529 CALL IWRTMA(ISTR_AB,1,NEL_AB,1,NEL_AB) 2530 CALL IWRTMA(ISTR_A,1,NEL_A,1,NEL_A) 2531 CALL IWRTMA(ISTR_B,1,NEL_B,1,NEL_B) 2532 END IF 2533* 2534 RETURN 2535 END 2536 SUBROUTINE EITP_TO_SPOXTP_REO(I_EI_TP_REO,ISPOBEX_TP,NSPOBEX_TP, 2537 & IB_EXTP, IB_INTP,NGAS,N_EXTP,N_INTP, 2538 & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX) 2539* 2540* Obtain reorder array of spinorbital excitation types from EI 2541* to standard order 2542* 2543*. Jeppe Olsen 2544* 2545*. October 2006 2546* 2547 INCLUDE 'wrkspc.inc' 2548*. Input : CAAB for all spinorbitalexcitations 2549 INTEGER ISPOBEX_TP(4,NGAS,*) 2550 INTEGER N_IN_FOR_EX(N_EXTP),IB_IN_FOR_EX(N_EXTP), I_IN_FOR_EX(*) 2551*. Output 2552 INTEGER I_EI_TP_REO(NSPOBEX_TP) 2553*. Local scratch 2554 INTEGER I_EI_MERGE(4*MXPNGAS) 2555* 2556 NTEST = 100 2557 IF(NTEST.GE.100) THEN 2558 WRITE(6,*) ' Info from EITP_TO_SPOXTP_REO ' 2559 WRITE(6,*) ' =============================' 2560 END IF 2561 IZERO = 0 2562 CALL ISETVC(I_EI_MERGE,IZERO,4*NGAS) 2563 I_EI_TP = 0 2564* 2565C? WRITE(6,*) ' NSPOBEX_TP, N_EXTP = ', NSPOBEX_TP, N_EXTP 2566 DO J_EXTP = 1, N_EXTP 2567C? WRITE(6,*) ' J_EXTP = ', J_EXTP 2568 L_INTP = N_IN_FOR_EX(J_EXTP) 2569 IB_JEX = IB_IN_FOR_EX(J_EXTP) 2570 DO JJ_INTP = 1, L_INTP 2571 I_EI_TP = I_EI_TP + 1 2572 J_INTP = I_IN_FOR_EX(IB_JEX-1+JJ_INTP) 2573 IF(NTEST.GE.100) THEN 2574 WRITE(6,*) ' J_EXTP, J_INTP, I_EI_TP = ', 2575 & J_EXTP, J_INTP, I_EI_TP 2576 END IF 2577 I_IN_ADR = IB_INTP - 1 + J_INTP 2578 I_EX_ADR = IB_EXTP - 1 + J_EXTP 2579C IVCSUM(IA,IB,IC,IFACB,IFACC,NDIM) 2580*. Merge Internal and external types 2581 IONE = 1 2582 IF(NTEST.GE.100) THEN 2583 WRITE(6,*) ' J_EXTP, J_INTP: ' 2584 CALL WRT_SPOX_TP(ISPOBEX_TP(1,1,I_EX_ADR),1) 2585 CALL WRT_SPOX_TP(ISPOBEX_TP(1,1,I_IN_ADR),1) 2586 END IF 2587* 2588 CALL IVCSUM(I_EI_MERGE, 2589 & ISPOBEX_TP(1,1,I_IN_ADR),ISPOBEX_TP(1,1,I_EX_ADR), 2590 & 1,1,4*NGAS) 2591 IF(NTEST.GE.100) THEN 2592 WRITE(6,*) ' I_EI_MERGE: ' 2593 CALL WRT_SPOX_TP(I_EI_MERGE,1) 2594 END IF 2595*. Find this type in SPOBEX 2596 IFOUND = 0 2597 DO JSPOBEX_TP = 1, NSPOBEX_TP 2598C COMPARE_TWO_INTARRAYS(IA,IB,NAB,IDENT) 2599 CALL COMPARE_TWO_INTARRAYS( 2600 & I_EI_MERGE,ISPOBEX_TP(1,1,JSPOBEX_TP),4*NGAS,IDENT) 2601 IF(IDENT.NE.0) IFOUND = JSPOBEX_TP 2602 END DO 2603 I_EI_TP_REO(I_EI_TP) = IFOUND 2604 END DO 2605 END DO 2606* 2607 IF(NTEST.GE.10) THEN 2608 WRITE(6,*) ' Reorder array for types: EI => standard order ' 2609 CALL IWRTMA(I_EI_TP_REO,1,NSPOBEX_TP,1,NSPOBEX_TP) 2610 END IF 2611* 2612 RETURN 2613 END 2614C CALL GET_INTERNAL_STATES(N_EXTOP_TP,N_INTOP_TP, 2615C & WORK(KLSOBEX),WORK(KL_N_INT_FOR_EXT),WORK(KL_IB_INT_FOR_EXT), 2616C & WORK(KL_I_INT_FOR_EXT),WORK(KL_NDIM_IN_SE), 2617C & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 2618C & WORK(KL_X1_INT_EI_FOR_SE), WORK(KL_X2_INT_EI_FOR_SE), 2619C & WORK(KL_SG_INT_EI_FOR_SE) 2620C & WORK(KL_IBX1_INT_EI_FOR_SE), WORK(KL_IBX2_INT_EI_FOR_SE), 2621C & WORK(KL_IBSG_INT_EI_FOR_SE) 2622C & WORK(KL_X2L_INT_EI_FOR_SE), 2623C & I_IN_TP,I_INT_OFF,I_EXT_OFF) 2624 SUBROUTINE GET_INTERNAL_STATES(N_EXTP,N_INTP,ISPOBEX, 2625 & N_IN_FOR_EX,IB_IN_FOR_EX,I_IN_FOR_EX, 2626 & N_IEL_FOR_SE,N_ORTN_FOR_SE,N_INT_FOR_SE, 2627 & X1_INT_EI_FOR_SE,X2_INT_EI_FOR_SE,SG_INT_EI_FOR_SE, 2628 & S_INT_EI_FOR_SE, 2629 & IBX1_INT_EI_FOR_SE,IBX2_INT_EI_FOR_SE, 2630 & IBSG_INT_EI_FOR_SE, IBS_INT_EI_FOR_SE, 2631 & X2L_INT_EI_FOR_SE, 2632 & I_INT_TP,IB_INTP,IB_EXTP) 2633* 2634*. Obtain the orthonormal internal states for each EI class. 2635*. If I_INT_TP = 1, then the internal states are the 2636* the nonorthogonal states of the metric 2637*. IF I_INT_TP = 2, then the internal states are 2638*. obtained by diagonalizing the internal Hamiltonian in the 2639*. space of orthonormal internal states 2640*. If I_INT_TP = 3, then the Jacobian in the internal space 2641*. is constructed and diagonalized and the left and right hand 2642*. side eigenvectors are obtained. 2643*. The transformation matrix between actual orthonormal and elementary internal states 2644* 2645* Currently a CASSCF state is assumed 2646* 2647*. October 28, 2007 in Crete, being at conference at Stines 28 years birthday.. 2648* 2649*. Continued March 17, 2008, on the train back from Warsaw 2650* 2651* Finalized(!) in Pisa Oct 6 (2008 (no comments, please...)) 2652*. Nope, modified march 12, 2009- and silence is still appreciated 2653* 2654*. Metric over internal states added Oct. 2009 2655* 2656*. Improving, Oct. 2012 2657* 2658* Last Modification; Oct. 18, 2012; Jeppe Olsen, Madrid; Batching to disc, etc 2659*.(Debugged in the air from Madrid to Amsterdam....) 2660*. 2661* 2662*. Form of internal Hamiltonian is defined by I_INT_HAM from CRUN 2663 INCLUDE 'wrkspc.inc' 2664 REAL*8 INPROD 2665* 2666 INCLUDE 'cgas.inc' 2667 INCLUDE 'gasstr.inc' 2668 INCLUDE 'lucinp.inc' 2669 INCLUDE 'ctcc.inc' 2670 INCLUDE 'cstate.inc' 2671 INCLUDE 'cands.inc' 2672 INCLUDE 'multd2h.inc' 2673 INCLUDE 'glbbas.inc' 2674 INCLUDE 'clunit.inc' 2675 INCLUDE 'oper.inc' 2676 INCLUDE 'cintfo.inc' 2677 INCLUDE 'crun.inc' 2678*. Input 2679 INTEGER N_IN_FOR_EX(N_EXTP),IB_IN_FOR_EX(N_EXTP),I_IN_FOR_EX(*) 2680*. Number of elementary internal operators of given sym 2681*. for given external type 2682 INTEGER N_IEL_FOR_SE(NSMOB,N_EXTP) 2683 INTEGER ISPOBEX(NGAS,4,*) 2684*. Output 2685*. Number of orthonormal internal states per symmetry and external type 2686 INTEGER N_ORTN_FOR_SE(NSMOB,N_EXTP) 2687*. Number of included internal states per symmetry and external types 2688 INTEGER N_INT_FOR_SE(NSMOB,N_EXTP) 2689*. Transformation matrix X1(EI,INT) Diagonalizing metric 2690*. various internal states 2691 DIMENSION X1_INT_EI_FOR_SE(*) 2692*. Nonvanishing eigenvalues sigma of metric 2693 DIMENSION SG_INT_EI_FOR_SE(*) 2694*. Overlap of internal states 2695 DIMENSION S_INT_EI_FOR_SE(*) 2696*. Tranformation matrix diagonalizing zero order Hamiltonian in 2697* orthonormal basis 2698 DIMENSION X2_INT_EI_FOR_SE(*) 2699*. lefthand side transformation for X2 if I_IN_TP = 3 2700 DIMENSION X2L_INT_EI_FOR_SE(*) 2701*. and offsets for X1 2702 INTEGER IBX1_INT_EI_FOR_SE(NSMOB,N_EXTP) 2703*. and offsets for X2 2704 INTEGER IBX2_INT_EI_FOR_SE(NSMOB,N_EXTP) 2705*. and offsets for Sigma 2706 INTEGER IBSG_INT_EI_FOR_SE(NSMOB,N_EXTP) 2707*. and offsets for S 2708 INTEGER IBS_INT_EI_FOR_SE(NSMOB,N_EXTP) 2709* 2710 NTEST = 10 2711 IF(NTEST.GE.10) THEN 2712 WRITE(6,*) ' ---------------------------- ' 2713 WRITE(6,*) ' GET_INTERNAL_STATES Speaking ' 2714 WRITE(6,*) ' ---------------------------- ' 2715 WRITE(6,*) 2716 WRITE(6,*) ' I_INT_TP, I_INT_HAM = ', I_INT_TP,I_INT_HAM 2717 END IF 2718* 2719 I12_SAVE = I12 2720 PSSIGN_SAVE = PSSIGN 2721 IDC_SAVE = IDC 2722* 2723 PSSIGN = 0.0D0 2724 IDC = 1 2725*. Incore or out-of core construction of matrices 2726 I_IN_OR_OUT = 2 2727 IF(I_IN_OR_OUT.EQ.1) THEN 2728 WRITE(6,*) ' In core construction of matrices' 2729 LUSCR_INT = -1 2730 LUSCR_INT2 = -1 2731 ELSE 2732 WRITE(6,*) ' Out of core construction of matrices' 2733C FILEMAN_MINI(IFILE,ITASK) 2734 CALL FILEMAN_MINI(LUSCR_INT,'ASSIGN') 2735 CALL FILEMAN_MINI(LUSCR_INT2,'ASSIGN') 2736 CALL REWINO(LUSCR_INT) 2737 CALL REWINO(LUSCR_INT2) 2738 END IF 2739* 2740 IB_X1 = 1 2741 IB_X2 = 1 2742 IB_SG = 1 2743 IB_S = 1 2744* 2745 CALL MEMMAN(IDUM,IDUM,'MARK ',1,'GT_INS') 2746 CALL LUCIAQENTER('GT_INS') 2747*. Scratch space : Two matrices S(INT1,INT2), where INT1, INT2 runs over 2748* internal states with given symmetry 2749* belonging to a given type of external state 2750* 2751*. Obtain occupations of alpha- and beta-strings of reference space - 2752*. is currently assumed to be a single space 2753C GET_REF_ALBE_OCC(IREFSPC,IREF_AL,IREF_BEC 2754 CALL MEMMAN(KL_REF_AL,NGAS,'ADDL ',2,'REF_AL') 2755 CALL MEMMAN(KL_REF_BE,NGAS,'ADDL ',2,'REF_BE') 2756 CALL MEMMAN(KL_IREF_AL,NGAS,'ADDL ',2,'IEF_AL') 2757 CALL MEMMAN(KL_IREF_BE,NGAS,'ADDL ',2,'IEF_BE') 2758*. NOTE : reference space is assumed to be space 1 2759 IREFSPC = 1 2760 CALL GET_REF_ALBE_OCC(IREFSPC,WORK(KL_REF_AL),WORK(KL_REF_BE)) 2761 IF(NTEST.GE.1000) THEN 2762 WRITE(6,*) ' Occupation of alpha- and beta-strings in reference' 2763 CALL IWRTMA(WORK(KL_REF_AL),1,NGAS,1,NGAS) 2764 CALL IWRTMA(WORK(KL_REF_BE),1,NGAS,1,NGAS) 2765 END IF 2766*Two matrices over internal states for given external type 2767 NDIM_II = IMXMN(1,N_IEL_FOR_SE,NSMOB*N_EXTP) 2768 IF(NTEST.GE.10) WRITE(6,*) 2769 &' Largest number of internal state for given sym and E-type', 2770 & NDIM_II 2771 CALL MEMMAN(KL_MATII1,NDIM_II**2,'ADDL ',2,'MATII1') 2772 CALL MEMMAN(KL_MATII2,NDIM_II**2,'ADDL ',2,'MATII2') 2773*. and two vectors over internal states 2774 CALL MEMMAN(KL_VECII1,NDIM_II,'ADDL ',2,'VECII1') 2775 CALL MEMMAN(KL_VECII2,NDIM_II,'ADDL ',2,'VECII2') 2776 2777*. In the following we are going to use/misuse the standard 2778*. spinorbital types stored as the first elements in the spinorbitalexcitation arrays. 2779*. Save the information usually stored there 2780*. Arrays 2781 CALL MEMMAN(KLSOBEX_SAVE,(NSPOBEX_TP+1)*NGAS*4,'ADDL ',1,'SPOXSV') 2782 CALL ICOPVE(ISPOBEX,WORK(KLSOBEX_SAVE),(NSPOBEX_TP+1)*NGAS*4) 2783 NSPOBEX_TP_SAVE = NSPOBEX_TP 2784*. Loop over the various EI spaces 2785 DO I_EXTP = 1, N_EXTP 2786 IF(NTEST.GE.10) THEN 2787 WRITE(6,*) 'Output for external type = ', I_EXTP 2788 WRITE(6,*) '===============================' 2789 END IF 2790 L_INTP = N_IN_FOR_EX(I_EXTP) 2791 IOFF = IB_IN_FOR_EX(I_EXTP) 2792C? WRITE(6,*) ' L_INTP, IOFF = ', L_INTP, IOFF 2793*. Prepare the arrays defining calculations for this space 2794 NSPOBEX_TP = L_INTP 2795 DO II_INTP = 1, L_INTP 2796 I_INTP = I_IN_FOR_EX(IOFF-1+II_INTP) 2797 IF(NTEST.GE.1000) THEN 2798 WRITE(6,*) ' II_INTP, I_INTP ', II_INTP, I_INTP 2799 WRITE(6,*) ' IB_INTP,NGAS =', IB_INTP,NGAS 2800 END IF 2801 CALL ICOPVE(ISPOBEX(1,1,IB_INTP-1+I_INTP), 2802 & ISPOBEX(1,1,II_INTP),4*NGAS) 2803 END DO 2804* 2805*---------------------------------- 2806*. Obtain metric in internal space 2807* --------------------------------- 2808* 2809* 2810* Type of C-coefficients will be type of reference 2811 ICATP = 1 2812 ICBTP = 2 2813* 2814 N_ALEL_C = NELFTP(ICATP) 2815 N_BEEL_C = NELFTP(ICBTP) 2816*. Find how action of internal operator changes occupation, 2817*. as all operators makes identical changes of occupations, 2818*. it is sufficient to look on the first operator 2819COLD I_INTP1 = I_IN_FOR_EX(IB_INTP-1+1) 2820 IALDEL = IELSUM(ISPOBEX(1,1,1),NGAS)-IELSUM(ISPOBEX(1,3,1),NGAS) 2821 IBEDEL = IELSUM(ISPOBEX(1,2,1),NGAS)-IELSUM(ISPOBEX(1,4,1),NGAS) 2822C? WRITE(6,*) ' IALDEL, IBEDEL = ', IALDEL, IBEDEL 2823*. Occupation of internal operator times reference space 2824C CAAB_T_ABSTR(ICAAB,IAOC_IN,IBOC_IN,IAOC_OUT,IBOC_OUT,NGAS ) 2825 CALL CAAB_T_ABSTR(ISPOBEX(1,1,1),WORK(KL_REF_AL),WORK(KL_REF_BE), 2826 & WORK(KL_IREF_AL),WORK(KL_IREF_BE),NGAS) 2827*. We now have occupation of resulting strings, find string numbers 2828C FIND_SPGRP_FROM_OCC(IOCC,ISPGRP_NUM,ITP) 2829*. Types of alpha and beta strings 2830 N_ALEL_S = N_ALEL_C + IALDEL 2831 N_BEEL_S = N_BEEL_C + IBEDEL 2832 IF(NTEST.GE.1000) THEN 2833 WRITE(6,*) ' Number of alpha-electrons in IOP*|ref> ',N_ALEL_S 2834 WRITE(6,*) ' Number of beta- electrons in IOP*|ref> ',N_BEEL_S 2835 END IF 2836*. Find supergroup types with these number of electrons 2837 ISATP = 0 2838 ISBTP = 0 2839 DO ISPGP_TP = 1, NSTTP 2840 IF(NELFTP(ISPGP_TP).EQ.N_ALEL_S) ISATP = ISPGP_TP 2841 IF(NELFTP(ISPGP_TP).EQ.N_BEEL_S) ISBTP = ISPGP_TP 2842 END DO 2843 IF(NTEST.GE.1000) THEN 2844 WRITE(6,*) ' a-,b-Types of internal operator times ref', 2845 & ISATP,ISBTP 2846 END IF 2847*. Below is not neccesary 2848 CALL FIND_SPGRP_FROM_OCC(WORK(KL_IREF_AL),IAL_SPGP_S,ISATP) 2849 CALL FIND_SPGRP_FROM_OCC(WORK(KL_IREF_BE),IBE_SPGP_S,ISBTP) 2850 IF(NTEST.GE.1000) THEN 2851 WRITE(6,*) ' a-b-supergroups of internal operator times ref', 2852 & ISATP,ISBTP 2853 END IF 2854* 2855 IF(ISATP.EQ.0.OR.ISBTP.EQ.0) THEN 2856 WRITE(6,*) 2857 & ' String type with modified electroncount does not exist' 2858 WRITE(6,*) ' NAEL, NBEL, ISATP, ISBTP = ', 2859 & N_ALEL_S, N_BEEL_S, ISATP, ISBTP 2860 STOP 2861 & ' String type with modified electroncount does not exist' 2862 END IF 2863* 2864*. Occupations of internal operator times reference space 2865* 2866* Generate a space with this occupation. This space is 2867* stored as the last CI space and combination space. If the numbers 2868* of these equals MXPICI I am trouble 2869* 2870 IF(NCISPC.EQ.MXPICI.OR.NCMBSPC.EQ.MXPICI) THEN 2871 WRITE(6,*) ' Not space for an extra CI space ' 2872 WRITE(6,*) ' Increase parameter MXPICI' 2873 WRITE(6,*) ' NCISPC, MXPICI ', NCISPC, MXPICI 2874 WRITE(6,*) ' NCMBSPC, MXPICI ', NCMBSPC, MXPICI 2875 STOP ' Increase parameter MXPICI' 2876 END IF 2877 IMODREF_CISPC = NCISPC + 1 2878 IMODREF_CMBSPC = NCMBSPC + 1 2879 LCMBSPC(IMODREF_CMBSPC) = 1 2880 ICMBSPC(1,IMODREF_CMBSPC) = IMODREF_CISPC 2881* 2882 IALTP_FOR_GAS(IMODREF_CMBSPC) = ISATP 2883 IBETP_FOR_GAS(IMODREF_CMBSPC) = ISBTP 2884*. Occupation constraints for modified reference : electrons are 2885*. added or removed from the active space. At the moment a single 2886* active space is assumed. 2887*. Copy reference space occupation 2888 IF(NTEST.GE.1000) THEN 2889 WRITE(6,*) ' Min and max for reference space ' 2890 CALL IWRTMA(IGSOCCX(1,1,IREFSPC),1,NGAS,1,NGAS) 2891 CALL IWRTMA(IGSOCCX(1,2,IREFSPC),1,NGAS,1,NGAS) 2892 END IF 2893* 2894 CALL ICOPVE(IGSOCCX(1,1,IREFSPC),IGSOCCX(1,1,IMODREF_CISPC), 2895 & NGAS) 2896 CALL ICOPVE(IGSOCCX(1,2,IREFSPC),IGSOCCX(1,2,IMODREF_CISPC), 2897 & NGAS) 2898*. The active space 2899 IACT_GAS = 0 2900 DO IGAS = 1, NGAS 2901 IF(IHPVGAS(IGAS).EQ.3) IACT_GAS = IGAS 2902 END DO 2903 IF(NTEST.GE.1000) WRITE(6,*) ' The active GASpace ', IACT_GAS 2904*. Change number of electrons in active space 2905 DO IGAS = 1, NGAS 2906 IF(IGAS.GE.IACT_GAS) THEN 2907 IGSOCCX(IGAS,1,IMODREF_CISPC) = 2908 & IGSOCCX(IGAS,1,IMODREF_CISPC) + IALDEL + IBEDEL 2909 IGSOCCX(IGAS,2,IMODREF_CISPC) = 2910 & IGSOCCX(IGAS,2,IMODREF_CISPC) + IALDEL + IBEDEL 2911 END IF 2912 END DO 2913* 2914 IF(NTEST.GE.1000) THEN 2915 WRITE(6,*) ' Min and max for modified reference space ' 2916 CALL IWRTMA(IGSOCCX(1,1,IMODREF_CISPC),1,NGAS,1,NGAS) 2917 CALL IWRTMA(IGSOCCX(1,2,IMODREF_CISPC),1,NGAS,1,NGAS) 2918 END IF 2919*. Loop over symmetry of internal operator 2920 DO IOPSM = 1, NSMOB 2921 IF(NTEST.GE.10) WRITE(6,'(A,2I5)') 2922 & ' Info for I_EXTP, IOPSM = ', I_EXTP, IOPSM 2923 IBX1_INT_EI_FOR_SE(IOPSM,I_EXTP) = IB_X1 2924 IBX2_INT_EI_FOR_SE(IOPSM,I_EXTP) = IB_X2 2925 IBSG_INT_EI_FOR_SE(IOPSM,I_EXTP) = IB_SG 2926 IBS_INT_EI_FOR_SE(IOPSM,I_EXTP) = IB_S 2927*.(IB_X1, IB_X2, IB_SG will be updated at end of loop over IOPSM) 2928*. Symmetry of operator times reference state 2929 IOPREFSM = MULTD2H(IOPSM,IREFSM) 2930 IF(NTEST.GE.1000) THEN 2931 WRITE(6,*) ' IOPSM, IREFSM, IOPREFSM = ', 2932 & IOPSM, IREFSM, IOPREFSM 2933 END IF 2934*. Number of determinants in modified internal space with symmetry ISYM 2935 LDET = LEN_CISPC(IMODREF_CMBSPC,IOPREFSM,NTEST) 2936*. Number of internal operators of this sym 2937 LINT = N_IEL_FOR_SE(IOPSM,I_EXTP) 2938*. Copy to N_INT_FOR_SE 2939 N_INT_FOR_SE(IOPSM,I_EXTP) = LINT 2940 2941 IF(NTEST.GE.10) 2942 & WRITE(6,*) ' IOPSM, LDET, LINT = ', IOPSM,LDET,LINT 2943 IF(LDET.EQ.0.OR.LINT.EQ.0) THEN 2944*. No novanishing states of this symmetry, so 2945 N_ORTN_FOR_SE(IOPSM,I_EXTP) = 0 2946 ELSE 2947*. Not trivial zero, so look further ( for some hundred of lines..) 2948*. Space for expansion Int |ref> in SD for all Internal states 2949 IDUM = 0 2950 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'INTMAT') 2951 IF(I_IN_OR_OUT.EQ.1) THEN 2952 LDIM = LINT*LDET 2953 ELSE 2954 LDIM = MAX(LINT,LDET) 2955 END IF 2956 CALL MEMMAN(KL_INTMAT,LDIM,'ADDL ',2,'INTMAT') 2957 CALL MEMMAN(KL_INTMAT2,LDIM,'ADDL ',2,'INTMT2') 2958 CALL MEMMAN(KL_INTOP,LINT**2,'ADDL ',2,'INTOP ') 2959 CALL MEMMAN(KL_INTOP2,LINT**2,'ADDL ',2,'INTOP2') 2960 CALL MEMMAN(KL_INTOPV,LINT,'ADDL ',2,'INTOP ') 2961* 2962 LWORK = 4*LINT 2963 CALL MEMMAN(KL_WORK,LWORK,'ADDL ',2,'INTOP ') 2964 CALL MEMMAN(KL_IWORK,LWORK,'ADDL ',2,'INTOP ') 2965* 2966 ICSPC = IREFSPC 2967 ISSPC = IMODREF_CISPC 2968 ICSM = IREFSM 2969 ISSM = IOPREFSM 2970 IF(I_IN_OR_OUT.EQ.2) CALL REWINO(LUSCR_INT) 2971 DO INTOP = 1, LINT 2972 ZERO = 0.0D0 2973 CALL SETVEC(WORK(KL_INTOP),ZERO,LINT) 2974 WORK(KL_INTOP-1+INTOP) = 1.0D0 2975 ICSPC = IREFSPC 2976 ISSPC = IMODREF_CISPC 2977 ICSM = IREFSM 2978 ISSM = IOPREFSM 2979 CALL REWINO(LUC) 2980 CALL REWINO(LUHC) 2981 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 2982 & WORK(KL_INTOP),1) 2983 CALL REWINO(LUHC) 2984 IF(NTEST.GE.10000) THEN 2985 WRITE(6,*) ' Result of SIGDEN_CC on LUHC ' 2986 CALL WRTVCD(WORK(KVEC1P),LUHC,1,-1) 2987 END IF 2988 CALL REWINO(LUHC) 2989 IF(I_IN_OR_OUT.EQ.1) THEN 2990 CALL FRMDSCN(WORK(KL_INTMAT+(INTOP-1)*LDET),-1,-1,LUHC) 2991 IF(NTEST.GE.10000) THEN 2992 WRITE(6,*) ' INTOP * |ref> for INTOP = ', INTOP 2993 CALL WRTMAT(WORK(KL_INTMAT+(INTOP-1)*LDET),1,LDET,1,LDET) 2994 END IF 2995 ELSE 2996 CALL COPVCD(LUHC,LUSCR_INT,WORK(KL_INTMAT),0,-1) 2997 END IF 2998C? STOP ' Enforced stop after first call to SIGDEN_CC' 2999 END DO ! Loop over internal states 3000 IF(NTEST.GE.1000) THEN 3001 WRITE(6,*) ' The matrix X(IDET,IELOP) ' 3002 IF(I_IN_OR_OUT.EQ.1) THEN 3003 CALL WRTMAT(WORK(KL_INTMAT),LDET,LINT,LDET,LINT) 3004 ELSE 3005 CALL REWINO(LUSCR_INT) 3006 DO INTOP = 1, LINT 3007 CALL WRTVCD(WORK(KL_INTMAT),LUSCR_INT,0,-1) 3008 END DO 3009 END IF 3010 END IF !NTEST is sufficiently large 3011*. Set up the overlap matrix S(I,J) = <ref!op(i)+ op(j)!ref> 3012 IF(I_IN_OR_OUT.EQ.1) THEN 3013*. In house construction 3014 IJ = 0 3015 DO I = 1, LINT 3016 DO J = 1, I 3017 IJ = IJ + 1 3018 WORK(KL_MATII1-1+I*(I-1)/2+J) = 3019 & INPROD(WORK(KL_INTMAT+(I-1)*LDET), 3020 & WORK(KL_INTMAT+(J-1)*LDET),LDET) 3021 END DO 3022 END DO 3023 ELSE 3024*. Disc based construction 3025*. Determine amount of memory that can be used for storing expansion 3026*. of external states 3027C MEMMAN(KBASE,KADD,TASK,IR,IDENT) 3028 KBASE = 0 3029 CALL MEMMAN(KBASE,KFREE,'FREE ',IDUM,'CHKFRE') 3030 KMAX = 100000000 3031 KDISC = MIN(KMAX,KFREE,LDET*LINT) 3032 WRITE(6,*) 3033 & ' Amount of memory to be used for batches of int. states ', 3034 & KDISC 3035*. Number of states per batch 3036 NSTA_BAT = KDISC/LDET 3037 IF(NSTA_BAT.EQ.0) THEN 3038 WRITE(6,*) 3039 & ' Problem: Batch cannot contain a single internal state' 3040 WRITE(6,*) ' KDISC, LDET = ', KDISC, LDET 3041 STOP 3042 & ' Problem: Batch cannot contain a single internal state' 3043 END IF 3044 NBAT = LINT/NSTA_BAT 3045 IF(NBAT*NSTA_BAT.LT.LINT) NBAT = NBAT + 1 3046 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'BATINT') 3047 CALL MEMMAN(KLCBAT,KDISC,'ADDL ',2,'CBAT ') 3048 DO JBAT = 1, NBAT 3049 J_INI = (JBAT-1)*NSTA_BAT + 1 3050 J_END = MIN(JBAT*NSTA_BAT, LINT) 3051 LBAT = J_END + 1 - J_INI 3052*. Read J-states in 3053 CALL SKPVCD(LUSCR_INT,J_INI-1,WORK(KL_INTMAT),1,-1) 3054 DO JSTA = 1, LBAT 3055 CALL FRMDSCN(WORK(KLCBAT+(JSTA-1)*LDET),-1,-1,LUSCR_INT) 3056 END DO 3057*. Loop over I-states and generate overlap 3058 CALL REWINO(LUSCR_INT) 3059 DO I = 1, LINT 3060 CALL FRMDSCN(WORK(KL_INTMAT),-1,-1,LUSCR_INT) 3061 DO J = J_INI, I 3062 IJ = I*(I-1)/2 + J 3063 WORK(KL_MATII1-1+I*(I-1)/2+J) = 3064 & INPROD(WORK(KL_INTMAT), 3065 & WORK(KLCBAT+(J-J_INI)*LDET),LDET) 3066 END DO 3067 END DO ! End of loop over pair of states 3068 END DO ! End of loop over batches of states 3069 CALL MEMMAN(IDUM,IDUN,'FLUSM ',IDUM,'BATINT') 3070 END IF ! in house or disc version 3071* 3072 IF(NTEST.GE.100) THEN 3073 WRITE(6,*) ' The overlap <ref!op(i)+ op(j)!ref>' 3074 CALL PRSYM(WORK(KL_MATII1),LINT) 3075 END IF 3076 CALL 3077 & COPVEC(WORK(KL_MATII1),S_INT_EI_FOR_SE(IB_S),LINT*(LINT+1)/2) 3078* 3079* ------------------- 3080*. Diagonalize metric 3081* ------------------- 3082* 3083*. Expand to complete form 3084C TRIPAK(AUTPAK,APAK,IWAY,MATDIM,NDIM) 3085 CALL TRIPAK(WORK(KL_MATII2),WORK(KL_MATII1),2,LINT,LINT) 3086*. Obtain nonsingular orthogonal internal states 3087C E CHK_S_FOR_SING(S,NDIM,NSING,X,SCR,SCR2) 3088 CALL CHK_S_FOR_SING2(WORK(KL_MATII2),LINT,NSING, 3089 & WORK(KL_MATII2),WORK(KL_VECII1),WORK(KL_VECII2), 3090 & THRES_SINGU) 3091 NNONSING = LINT - NSING 3092 IF(NTEST.GE.10) 3093 & WRITE(6,*) ' Number of singular and nonsingular states ', 3094 & NSING,NNONSING 3095 N_ORTN_FOR_SE(IOPSM,I_EXTP) = NNONSING 3096*. The nonsingular basis are the last NNONSING vectors in WORK(KLMATII2) 3097*. copy these to first X1 3098 DO IORTN = 1, NNONSING 3099 CALL COPVEC(WORK(KL_MATII2+(NSING+IORTN-1)*LINT), 3100 & X1_INT_EI_FOR_SE(IB_X1+(IORTN-1)*LINT),LINT) 3101 SG_INT_EI_FOR_SE(IB_SG-1+IORTN) = 3102 & WORK(KL_VECII1-1+NSING+IORTN) 3103 END DO 3104* 3105 IF(NTEST.GE.1000) THEN 3106 WRITE(6,*) ' internal states diagonalizing metric' 3107 CALL WRTMAT(X1_INT_EI_FOR_SE(IB_X1), 3108 & LINT,NNONSING,LINT,NNONSING) 3109 END IF 3110 IF(NTEST.GE.20) THEN 3111 WRITE(6,*) ' Eigenvalues of metric' 3112 CALL WRTMAT(SG_INT_EI_FOR_SE(IB_SG),1,NNONSING,1,NNONSING) 3113 END IF 3114 IF(NTEST.GE.20) THEN 3115 WRITE(6,*) 3116 & ' Info for construction of matrices over orth. states' 3117 END IF 3118* 3119 IF(I_INT_TP.GE.2) THEN 3120*. Obtain - if required - zero-order Hamiltonian in internal 3121*. space and diagonalize this 3122C. H0 Intop |ref> for all orthonormal states 3123 IF(I_IN_OR_OUT.EQ.2) CALL REWINO(LUSCR_INT2) 3124 IF(I_IN_OR_OUT.EQ.2) CALL REWINO(LUSCR_INT) 3125 DO IORTN = 1, NNONSING 3126 IF(NTEST.GE.1000) WRITE(6,*) ' IORTN = ', IORTN 3127*. Scale to obtain orthonormal state 3128 FACTOR = 1.0D0/SQRT(SG_INT_EI_FOR_SE(IB_SG-1+IORTN)) 3129 CALL COPVEC(X1_INT_EI_FOR_SE(IB_X1+(IORTN-1)*LINT), 3130 & WORK(KL_VECII1),LINT) 3131 CALL SCALVE(WORK(KL_VECII1),FACTOR,LINT) 3132* 3133 ICSPC = IREFSPC 3134 ISSPC = IMODREF_CMBSPC 3135 ICSM = IREFSM 3136 ISSM = IOPREFSM 3137C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 3138 CALL REWINO(LUHC) 3139 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 3140 & WORK(KL_VECII1),1) 3141 IF(NTEST.GE.10000) THEN 3142 WRITE(6,*) ' After sigden ' 3143 WRITE(6,*) ' Resulting vector : ' 3144 CALL WRTVCD(WORK(KVEC1P),LUHC,1,-1) 3145 END IF 3146 IF(I_IN_OR_OUT.EQ.1) THEN 3147 CALL REWINO(LUHC) 3148 CALL FRMDSCN(WORK(KL_INTMAT+(IORTN-1)*LDET),-1,-1,LUHC) 3149 IF(NTEST.GE.10000) THEN 3150 WRITE(6,*) ' Intop (IORTN) |ref>' 3151 CALL WRTMAT(WORK(KL_INTMAT+(IORTN-1)*LDET),1,LDET,1,LDET) 3152 END IF 3153 ELSE 3154 CALL REWINO(LUHC) 3155 CALL COPVCD(LUHC,LUSCR_INT,WORK(KL_INTMAT),0,-1) 3156 END IF 3157*. zero-order Hamiltonian times Int(iortn)|ref> 3158 ICSPC = IMODREF_CMBSPC 3159 ISSPC = IMODREF_CMBSPC 3160 ICSM = IOPREFSM 3161 ISSM = IOPREFSM 3162 IF(I_INT_HAM.EQ.1) THEN 3163C? WRITE(6,*) 3164C? & ' One-body H0 used for internal zero-order states' 3165 I12 = 1 3166 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 3167 ELSE 3168C? WRITE(6,*) 3169C? & ' Two-body H used for internal zero-order states' 3170 END IF 3171* 3172 CALL REWINO(LUHC) 3173 CALL REWINO(LUSC51) 3174 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51,0,0) 3175 CALL MEMCHK2('AFTMV7') 3176 IF(NTEST.GE.10000) THEN 3177 WRITE(6,*) ' IMODREF_CMBSPC =', IMODREF_CMBSPC 3178 WRITE(6,*) ' After mv7 ' 3179 WRITE(6,*) ' Resulting vector : ' 3180 CALL WRTVCD(WORK(KVEC1P),LUSC51,1,-1) 3181 END IF 3182 IF(I_INT_HAM.EQ.1) 3183 & CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 3184 CALL REWINO(LUSC51) 3185 IF(I_IN_OR_OUT.EQ.1) THEN 3186 CALL FRMDSCN(WORK(KL_INTMAT2+(IORTN-1)*LDET),-1,-1,LUSC51) 3187 ELSE 3188 CALL COPVCD(LUSC51,LUSCR_INT2,WORK(KL_INTMAT),0,-1) 3189 END IF 3190 END DO 3191* 3192 IF(NTEST.GE.1000) THEN 3193 WRITE(6,*) ' Oint!0> and H0 Oint!0> expansions:' 3194 IF(I_IN_OR_OUT.EQ.1) THEN 3195 CALL WRTMAT(WORK(KL_INTMAT),LDET,NNONSING, 3196 & LDET,NNONSING) 3197 WRITE(6,*) 3198 CALL WRTMAT(WORK(KL_INTMAT2),LDET,NNONSING, 3199 & LDET,NNONSING) 3200 WRITE(6,*) 3201 ELSE 3202 CALL REWINO(LUSCR_INT) 3203 DO I = 1, NNONSING 3204 CALL WRTVCD(WORK(KVEC1P),LUSCR_INT,0,-1) 3205 END DO 3206 WRITE(6,*) 3207 CALL REWINO(LUSCR_INT2) 3208 DO I = 1, NNONSING 3209 CALL WRTVCD(WORK(KVEC1P),LUSCR_INT2,0,-1) 3210 END DO 3211 END IF 3212 END IF 3213 3214*. <ref|intop+ H+ intop|ref>: for test complete matrix and metric are 3215*. calculated 3216 IF(I_IN_OR_OUT.EQ.1) THEN 3217 DO IORTN = 1, NNONSING 3218 DO JORTN = 1, NNONSING 3219 IJ = (JORTN-1)*NNONSING + IORTN 3220*. H0 3221 WORK(KL_INTOP + IJ - 1) = 3222 & INPROD(WORK(KL_INTMAT2+(JORTN-1)*LDET), 3223 & WORK(KL_INTMAT+(IORTN-1)*LDET),LDET) 3224*. S 3225 WORK(KL_INTOP2+ IJ - 1) = 3226 & INPROD(WORK(KL_INTMAT+(JORTN-1)*LDET), 3227 & WORK(KL_INTMAT+(IORTN-1)*LDET),LDET) 3228 END DO 3229 END DO 3230 ELSE 3231*. Disk version 3232*. Determine amount of memory that can be used for storing expansion 3233*. of external states 3234C MEMMAN(KBASE,KADD,TASK,IR,IDENT) 3235 KBASE = 0 3236 CALL MEMMAN(KBASE,KFREE,'FREE ',IDUM,'CHKFRE') 3237 KMAX = 100000000 3238 KDISC = MIN(KMAX,KFREE,LDET*LINT) 3239 WRITE(6,*) 3240 & ' Amount of memory to be used for batches of int. states ', 3241 & KDISC 3242*. Number of states per batch 3243 NSTA_BAT = KDISC/LDET 3244 IF(NSTA_BAT.EQ.0) THEN 3245 WRITE(6,*) 3246 & ' Problem: Batch cannot contain a single internal state' 3247 WRITE(6,*) ' KDISC, LDET = ', KDISC, LDET 3248 STOP 3249 & ' Problem: Batch cannot contain a single internal state' 3250 END IF 3251 NBAT = NNONSING/NSTA_BAT 3252 IF(NBAT*NSTA_BAT.LT.NNONSING) NBAT = NBAT + 1 3253 IF(NTEST.GE.1000) THEN 3254 WRITE(6,*) ' NSTA_BAT, NBAT = ', 3255 & NSTA_BAT, NBAT 3256 END IF 3257 CALL MEMMAN(IDUM,IDUN,'MARK ',IDUM,'BATINT') 3258 CALL MEMMAN(KLCBAT,KDISC,'ADDL ',2,'CBAT ') 3259 DO JBAT = 1, NBAT 3260 J_INI = (JBAT-1)*NSTA_BAT + 1 3261 J_END = MIN(JBAT*NSTA_BAT, NNONSING) 3262 LBAT = J_END + 1 - J_INI 3263 IF(NTEST.GE.1000) THEN 3264 WRITE(6,*) ' JBAT, J_INI, J_END = ', 3265 & JBAT, J_INI, J_END 3266 END IF 3267*. Read states in 3268 CALL SKPVCD(LUSCR_INT,J_INI-1,WORK(KL_INTMAT),1,-1) 3269 IF(NTEST.GE.1000) WRITE(6,*) ' SKPVCD passed ' 3270 DO JSTA = 1, LBAT 3271 CALL FRMDSCN(WORK(KLCBAT+(JSTA-1)*LDET),-1,-1,LUSCR_INT) 3272 IF(NTEST.GE.1000) 3273 & WRITE(6,*) ' FRMDSCN passed for JSTA =', JSTA 3274 END DO 3275* <I!HJ> 3276 CALL REWINO(LUSCR_INT2) 3277 DO I = 1, NNONSING 3278 CALL FRMDSCN(WORK(KL_INTMAT),-1,-1,LUSCR_INT2) 3279 IF(NTEST.GE.1000) WRITE(6,*) 3280 & ' FRMDSCN, LUSCR_INT2 passed for I = ', I 3281 DO J = J_INI, J_END 3282 IJ = (J-1)*NNONSING + I 3283 WORK(KL_INTOP-1+IJ) = 3284 & INPROD(WORK(KL_INTMAT),WORK(KLCBAT+(J-J_INI)*LDET), 3285 & LDET) 3286 END DO 3287 END DO 3288*.<I!J> 3289 CALL REWINO(LUSCR_INT) 3290 DO I = 1, NNONSING 3291 CALL FRMDSCN(WORK(KL_INTMAT),-1,-1,LUSCR_INT) 3292 DO J = J_INI, J_END 3293 IJ = (J-1)*NNONSING+ I 3294 WORK(KL_INTOP2-1+IJ) = 3295 & INPROD(WORK(KL_INTMAT),WORK(KLCBAT+(J-J_INI)*LDET), 3296 & LDET) 3297 END DO 3298 END DO ! End of loop over pair of states 3299 END DO ! End of loop over batches of states 3300 CALL MEMMAN(IDUM,IDUN,'FLUSM ',IDUM,'BATINT') 3301 END IF ! disk version 3302* 3303 IF(NTEST.GE.100) THEN 3304 WRITE(6,*) ' Metric over orthonormal states ' 3305 CALL WRTMAT(WORK(KL_INTOP2),NNONSING,NNONSING, 3306 & NNONSING,NNONSING) 3307 WRITE(6,*) ' H0 over orthonormal states ' 3308 CALL WRTMAT(WORK(KL_INTOP),NNONSING,NNONSING, 3309 & NNONSING,NNONSING) 3310 END IF 3311 IF(I_INT_TP.EQ.2.AND.NNONSING.NE.0) THEN 3312*. Diagonalize zero-Hamiltonian 3313C DIAG_SYMMAT_EISPACK(A,EIGVAL,SCRVEC,NDIM,IRETURN) 3314 CALL DIAG_SYMMAT_EISPACK(WORK(KL_INTOP),WORK(KL_VECII1), 3315 & WORK(KL_VECII2),NNONSING,IRETURN) 3316*. In output eigenvectors are in WORK(KL_INTOP) and eigenvalues are in 3317*. WORK(KL_VECII1) 3318 CALL COPVEC(WORK(KL_INTOP),X2_INT_EI_FOR_SE(IB_X2), 3319 & NNONSING*NNONSING) 3320 IF(NTEST.GE.20) THEN 3321 WRITE(6,*) ' Energies of internal states ' 3322 CALL WRTMAT(WORK(KL_VECII1),1,NNONSING,1,NNONSING) 3323 ELSE IF(NTEST.GE.10.AND.NNONSING.GE.1) THEN 3324 WRITE(6,*) ' Energy of lowest internal state ' 3325 CALL WRTMAT(WORK(KL_VECII1),1,1,1,1) 3326 END IF 3327 IF(NTEST.GE.1000) THEN 3328 WRITE(6,*) 3329 & ' Transformation from zero-order to orthonormal states' 3330 CALL WRTMAT(X2_INT_EI_FOR_SE(IB_X2), 3331 & NNONSING,NNONSING,NNONSING,NNONSING) 3332 END IF 3333 ELSE IF(I_INT_TP.EQ.3) THEN 3334 IF(I_IN_OR_OUT.EQ.2) THEN 3335 WRITE(6,*) ' Disc version for I_INT_TP = 3 not programmed' 3336 STOP ' Disc version for I_INT_TP = 3 not programmed' 3337 END IF 3338*. <0!O+I [H0,OJ]|0> should be obtained and diagonalized 3339*. On input WORK(KL_INTOP) contains <0!O+I H0 OJ|0>. 3340*. Obtain <0!O+I OJ H0 |0> 3341*. H0 |0> on LUHC 3342 ICSPC = IREFSPC 3343 ISSPC = IREFSPC 3344 ICSM = IREFSM 3345 ISSM = IREFSM 3346 IF(I_INT_HAM.EQ.1) THEN 3347C? WRITE(6,*) 3348C? & ' One-body H0 used for internal zero-order states' 3349 I12 = 1 3350 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 3351 ELSE 3352C? WRITE(6,*) 3353C? & ' Two-body H used for internal zero-order states' 3354 END IF 3355 CALL REWINO(LUC) 3356 CALL REWINO(LUHC) 3357 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC,0,0) 3358 IF(I_INT_HAM.EQ.1) 3359 & CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 3360*. Generate OI H0 |ref> 3361 DO IORTN = 1, NNONSING 3362 IF(NTEST.GE.1000) WRITE(6,*) ' IORTN = ', IORTN 3363*. Scale to obtain orthonormal state 3364 FACTOR = 1.0D0/SQRT(SG_INT_EI_FOR_SE(IB_SG-1+IORTN)) 3365 CALL COPVEC(X1_INT_EI_FOR_SE(IB_X1+(IORTN-1)*LINT), 3366 & WORK(KL_VECII1),LINT) 3367 CALL SCALVE(WORK(KL_VECII1),FACTOR,LINT) 3368 ICSPC = IREFSPC 3369 ISSPC = IMODREF_CMBSPC 3370C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 3371 CALL REWINO(LUSC51) 3372 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51, 3373 & WORK(KL_VECII1),1) 3374 IF(NTEST.GE.10000) THEN 3375 WRITE(6,*) ' After sigden ' 3376 WRITE(6,*) ' Resulting vector : ' 3377 CALL WRTVCD(WORK(KVEC1P),LUHC,1,-1) 3378 END IF 3379 CALL FRMDSCN(WORK(KL_INTMAT2+(IORTN-1)*LDET),-1,-1,LUSC51) 3380 END DO 3381*. Generate OI |ref> 3382 DO JORTN = 1, NNONSING 3383 IF(NTEST.GE.10000) WRITE(6,*) ' JORTN = ', JORTN 3384*. Scale to obtain orthonormal state 3385 FACTOR = 1.0D0/SQRT(SG_INT_EI_FOR_SE(IB_SG-1+IORTN)) 3386 CALL COPVEC(X1_INT_EI_FOR_SE(IB_X1+(JORTN-1)*LINT), 3387 & WORK(KL_VECII1),LINT) 3388 CALL SCALVE(WORK(KL_VECII1),FACTOR,LINT) 3389 ICSPC = IREFSPC 3390 ISSPC = IMODREF_CMBSPC 3391C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 3392 CALL REWINO(LUHC) 3393 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 3394 & WORK(KL_VECII1),1) 3395 CALL REWINO(LUSC51) 3396 CALL FRMDSCN(WORK(KL_INTMAT+(JORTN-1)*LDET),-1,-1,LUSC51) 3397 END DO 3398* Obtain <ref!o+i [H0,o j]|ref> = <ref!o+i H0 o j|ref> - 3399* <ref!o+i o j H0|ref> 3400 DO IORTN = 1, NNONSING 3401 DO JORTN = 1, NNONSING 3402 IJ = (JORTN-1)*NNONSING + IORTN 3403*. H0 3404 WORK(KL_INTOP + IJ - 1) = 3405 & WORK(KL_INTOP + IJ - 1) - 3406 & INPROD(WORK(KL_INTMAT2+(JORTN-1)*LDET), 3407 & WORK(KL_INTMAT+(IORTN-1)*LDET),LDET) 3408 END DO 3409 END DO 3410* 3411 IF(NTEST.GE.100) THEN 3412 WRITE(6,*) 3413 & ' <ref!o+i [H0,o j]|ref> over orthonormal states' 3414 CALL WRTMAT(WORK(KL_INTOP),NNONSING,NNONSING, 3415 & NNONSING,NNONSING) 3416 END IF 3417*. Diagonalize: 3418*. Obtain eigenvalues in WORK(KL_VECII1),lefteigenvectors in 3419*. X2L,righteigenvectors in X2.., work(KL_MATII1), 3420*- ONLY PROGRAMMED FOR REAL EIGENVALUES... 3421 CALL GET_LR_EIGVEC_GENMAT(WORK(KL_INTOP),NNONSING, 3422 & X2_INT_EI_FOR_SE(IB_X2),X2L_INT_EI_FOR_SE(IB_X1), 3423 & WORK(KL_VECII1), 3424 & WORK(KL_VECII2),WORK(KL_WORK),WORK(KL_IWORK), 3425 & LWORK,WORK(KL_MATII1)) 3426C GET_LR_EIGVEC_GENMAT(A,NDIM,VCR,VCL,VLR, 3427C & VLI,WORK,IWORK,LWORK,SCRMAT) 3428 END IF 3429*. ^ End if I_INT_TP = 3 3430 END IF 3431* ^ End if I_INT_TP.GE.2 3432 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'INTMAT') 3433 IB_X1 = IB_X1 + LINT*NNONSING 3434 IB_X2 = IB_X2 + NNONSING*NNONSING 3435 IB_SG = IB_SG + NNONSING 3436 IB_S = IB_S + LINT*(LINT+1)/2 3437 END IF 3438* ^ End if there was a nonvanishing number of EI states 3439 END DO 3440* ^ End of loop over symmetries 3441 END DO 3442* ^ End of loop over external states 3443*- Clean-up time: Move Spin-orbital excitation types back 3444 NSPOBEX_TP = NSPOBEX_TP_SAVE 3445 CALL ICOPVE(WORK(KLSOBEX_SAVE),ISPOBEX,(NSPOBEX_TP+1)*NGAS*4) 3446 I12 = I12_SAVE 3447 PSSIGN = PSSIGN_SAVE 3448 IDC = IDC_SAVE 3449* 3450 CALL MEMMAN(IDUM,IDUM,'FLUSM',1,'GT_INS') 3451* 3452 IF(I_IN_OR_OUT.EQ.2) THEN 3453C FILEMAN_MINI(IFILE,ITASK) 3454 CALL FILEMAN_MINI(LUSCR_INT,'FREE ') 3455 CALL FILEMAN_MINI(LUSCR_INT2,'FREE ') 3456 END IF 3457* 3458 CALL LUCIAQEXIT('GT_INS') 3459 RETURN 3460 END 3461 FUNCTION IMXMN(MAX_OR_MIN,IVEC,NELMNT) 3462* 3463*. Find largest (MAX_OR_MIN = 1) or smallest (MAX_OR_MIN = 2) element 3464*. in integer vector IVEC 3465* 3466 INCLUDE 'implicit.inc' 3467 INTEGER IVEC(NELMNT) 3468* 3469 IVAL = IVEC(1) 3470 IF(MAX_OR_MIN.EQ.1) THEN 3471 DO I = 2, NELMNT 3472 IVAL = MAX(IVAL,IVEC(I)) 3473 END DO 3474 ELSE 3475 DO I = 2, NELMNT 3476 IVAL = MIN(IVAL,IVEC(I)) 3477 END DO 3478 END IF 3479* 3480 IMXMN = IVAL 3481* 3482 NTEST = 100 3483 IF(NTEST.GE.100) THEN 3484 IF(MAX_OR_MIN.EQ.1) THEN 3485 WRITE(6,*) ' Largest element found by IMXMN ', IMXMN 3486 ELSE 3487 WRITE(6,*) ' Smallest element found by IMXMN ', IMXMN 3488 END IF 3489 END IF 3490* 3491 RETURN 3492 END 3493 SUBROUTINE CAAB_T_ABSTR(ICAAB,IAOC_IN,IBOC_IN, 3494 & IAOC_OUT,IBOC_OUT,NGAS ) 3495* 3496* A CAAB operator ICAAB and occupation of alpha and betastrings 3497* IAOC_IN, IBOC_IN are given 3498* 3499*. Find occupation of resulting strings 3500* 3501* Jeppe Olsen, March 17, 2007, trying once again to get momentum to MRCC 3502* 3503 INCLUDE 'implicit.inc' 3504*. Input 3505 INTEGER ICAAB(NGAS,4), IAOC_IN(NGAS),IBOC_IN(NGAS) 3506*. Output 3507 INTEGER IAOC_OUT(NGAS),IBOC_OUT(NGAS) 3508* 3509 DO IGAS = 1, NGAS 3510 IAOC_OUT(IGAS) = IAOC_IN(IGAS) + ICAAB(IGAS,1)-ICAAB(IGAS,3) 3511 IBOC_OUT(IGAS) = IBOC_IN(IGAS) + ICAAB(IGAS,2)-ICAAB(IGAS,4) 3512 END DO 3513* 3514 NTEST = 00 3515 IF(NTEST.GE.100) THEN 3516 WRITE(6,*) ' Output from CAAB_T_ABSTR ' 3517 WRITE(6,*) ' Input CAAB type ' 3518 CALL WRT_SPOX_TP(ICAAB,1) 3519 WRITE(6,*) ' Input alpha- and beta-types' 3520 CALL IWRTMA(IAOC_IN,1,NGAS,1,NGAS) 3521 CALL IWRTMA(IBOC_IN,1,NGAS,1,NGAS) 3522 END IF 3523* 3524 RETURN 3525 END 3526 FUNCTION ISQELSUM(IVEC,NVEC,ISYM) 3527* 3528* ISYM = 0: SUM_I IVEC(I)*IVEC(I) 3529* ISYM = 1: SUM_I IVEC(I)*(IVEC(I)+1)/2 3530* 3531*. Jeppe Olsen, March 2007 3532* 3533 INCLUDE 'implicit.inc' 3534*. Input 3535 INTEGER IVEC(NVEC) 3536* 3537 ISUM = 0 3538 IF(ISYM.EQ.0) THEN 3539 DO I = 1, NVEC 3540 ISUM = ISUM + IVEC(I)**2 3541 END DO 3542 ELSE 3543 DO I = 1, NVEC 3544 ISUM = ISUM + IVEC(I)*(IVEC(I)+1)/2 3545 END DO 3546 END IF 3547 ISQELSUM = ISUM 3548* 3549 NTEST = 00 3550 IF(NTEST.GE.100) THEN 3551 WRITE(6,*) ' ISQELSUM : sum of squared elements = ', ISUM 3552 END IF 3553* 3554 RETURN 3555 END 3556 FUNCTION LEN_CISPC(JCMBSPC,ISYM,IPRNT) 3557* 3558* Number of dets and combinations for given sym and combination space 3559* 3560* Jeppe Olsen, obtained from LCISPC 3561* 3562* =================== 3563*.Input common blocks 3564* =================== 3565* 3566 INCLUDE 'wrkspc.inc' 3567 INCLUDE 'lucinp.inc' 3568 INCLUDE 'cstate.inc' 3569 INCLUDE 'strinp.inc' 3570 INCLUDE 'strbas.inc' 3571 INCLUDE 'csm.inc' 3572 INCLUDE 'stinf.inc' 3573 INCLUDE 'cgas.inc' 3574 INCLUDE 'gasstr.inc' 3575 INCLUDE 'cicisp.inc' 3576* 3577 IDUM = 0 3578 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'LNCISP') 3579 CALL LUCIAQENTER('LCISP') 3580*. Obtain types 3581 ICISPC1 = ICMBSPC(1,JCMBSPC) 3582*. Type of alpha- and beta strings 3583 IF(ICISPC1.LE.NCISPC) THEN 3584 IATP = 1 3585 IBTP = 2 3586 ELSE 3587 IATP = IALTP_FOR_GAS(ICISPC1) 3588 IBTP = IBETP_FOR_GAS(ICISPC1) 3589 END IF 3590* 3591 NOCTPA = NOCTYP(IATP) 3592 NOCTPB = NOCTYP(IBTP) 3593* 3594 IOCTPA = IBSPGPFTP(IATP) 3595 IOCTPB = IBSPGPFTP(IBTP) 3596*.Local memory 3597 CALL MEMMAN(KLBLTP,NSMST,'ADDL ',2,'KLBLTP') 3598 IF(IDC.EQ.3 .OR. IDC .EQ. 4 ) 3599 &CALL MEMMAN(KLCVST,NSMST,'ADDL ',2,'KLCVST') 3600 CALL MEMMAN(KLIOIO,NOCTPA*NOCTPB, 'ADDL ',2,'KLIOIO') 3601*. Obtain array giving symmetry of sigma v reflection times string 3602*. symmetry. 3603 IF(IDC.EQ.3.OR.IDC.EQ.4) 3604 &CALL SIGVST(WORK(KLCVST),NSMST) 3605* 3606*. Note : size of max blocks not recalculated 3607*. Array defining symmetry combinations of internal strings 3608 CALL SMOST(NSMST,NSMCI,MXPCSM,ISMOST) 3609*. allowed combination of types 3610 CALL IAIBCM(JCMBSPC,WORK(KLIOIO)) 3611 CALL ZBLTP(ISMOST(1,ISYM),NSMST,IDC,WORK(KLBLTP),WORK(KLCVST)) 3612 CALL NGASDT(IGSOCCX(1,1,1),IGSOCCX(1,2,1),NGAS,ISYM, 3613 & NSMST,NOCTPA,NOCTPB,WORK(KNSTSO(IATP)), 3614 & WORK(KNSTSO(IBTP)), 3615 & ISPGPFTP(1,IBSPGPFTP(IATP)), 3616 & ISPGPFTP(1,IBSPGPFTP(IBTP)),MXPNGAS,NELFGP, 3617 & NCOMB,XNCOMB,MXS,MXSOO,WORK(KLBLTP),NTTSBL, 3618 & LCOL,WORK(KLIOIO),MXSOO_AS,XMXSOO,XMXSOO_AS) 3619* 3620 LEN_CISPC = NCOMB 3621* 3622 NTEST = 0 3623 NTEST = MAX(NTEST,IPRNT) 3624 IF(NTEST.GE.100) THEN 3625 WRITE(6,*) ' CMB space and symmetry ', JCMBSPC,ISYM 3626 WRITE(6,*) ' Number of determinants ', LEN_CISPC 3627 END IF 3628* 3629 CALL LUCIAQEXIT('LCISP') 3630 CALL MEMMAN(IDUM,IDUM,'FLUSM',IDUM,'LNCISP') 3631* 3632 RETURN 3633 END 3634 SUBROUTINE FIND_TYPSTR_WITH_TOTOCC(NEL,ITYP) 3635* 3636* Find type with NEL electrons. The program finds 3637* the last set type with given number of electrons 3638* 3639*. Jeppe Olsen, Billund airport, March 22, 2006 3640* 3641 INCLUDE 'wrkspc.inc' 3642 INCLUDE 'gasstr.inc' 3643#include "global.fh" 3644* 3645 ITYP = 0 3646 DO ISTTP = 1, NSTTP 3647 IF(NELFTP(ISTTP).EQ.NEL) ITYP = ISTTP 3648 END DO 3649* 3650 IF(ISTTP.EQ.0) THEN 3651 WRITE(6,*) 3652 & ' Error, stringtype with given number of elecs not found' 3653 WRITE(6,*) ' Required number of electrons = ', NEL 3654 STOP 3655 & ' Error, stringtype with given number of elecs not found' 3656 END IF 3657* 3658 NTEST = 00 3659 IF(NTEST.GE.100.and.ga_nodeid().eq.0) THEN 3660 WRITE(6,*) ' Number of electrons and stringtype ', 3661 & NEL, ITYP 3662 END IF 3663* 3664 RETURN 3665 END 3666 SUBROUTINE MATRIX_TIMES_SPARSEVEC 3667 &(A,VECELM,IVECIND,NDIM,NVEC,VECOUT) 3668* 3669*. VECOUT(I) = sum_(k=1,NVEC) A(I,(IVECIND(K))*VECELM(K) 3670* 3671 INCLUDE 'implicit.inc' 3672*. Input 3673 DIMENSION A(NDIM,*) 3674 DIMENSION VECELM(NVEC),IVECIND(NVEC) 3675*. Output 3676 DIMENSION VECOUT(NDIM) 3677* 3678 ZERO = 0.0D0 3679 CALL SETVEC(VECOUT,ZERO,NDIM) 3680* 3681 DO K = 1, NVEC 3682 KCOL = IVECIND(K) 3683 COEF = VECELM(K) 3684 DO I = 1, NDIM 3685 VECOUT(I) = VECOUT(I) + COEF*A(I,KCOL) 3686 END DO 3687 END DO 3688* 3689 RETURN 3690 END 3691 SUBROUTINE GET_LR_EIGVEC_GENMAT(A,NDIM,VCR,VCL,VLR, 3692 & VLI,WORK,IWORK,LWORK,SCRMAT) 3693* 3694*. Obtain eigenvalues and left and right eigenvectors of real 3695*. general matrix - outer routine for DGEEV 3696* 3697* Has not been generalized to complex eigenvalues 3698* and eigenvectors 3699* 3700*. Jeppe Olsen, September 07 3701* 3702 INCLUDE 'implicit.inc' 3703*. Input 3704 DIMENSION A(NDIM,NDIM) 3705*. Output 3706 DIMENSION VCR(NDIM,NDIM),VCL(NDIM,NDIM) 3707 DIMENSION VLR(NDIM),VLI(NDIM) 3708*. Scratch: WORK: Atleast 4*NDIM, preferably longer, IWORK 3709 DIMENSION WORK(LWORK), IWORK(LWORK) 3710*. And a scratch matrix - used for constructing overlap of eigenvectors 3711 DIMENSION SCRMAT(NDIM*NDIM) 3712* 3713 REAL*8 INPROD 3714 INFO = 0 3715 CALL DGEEV('V','V',NDIM,A,NDIM,VLR,VLI,VCL,NDIM,VCR,NDIM, 3716 & WORK,LWORK,INFO) 3717 IF(INFO.NE.0) THEN 3718 WRITE(6,*) ' Error in GET_LT_EIGVEC_GENMAT' 3719 WRITE(6,*) ' INFO=',INFO 3720 STOP ' Error in GET_LT_EIGVEC_GENMAT' 3721 END IF 3722* 3723 NTEST = 100 3724 IF(NTEST.GE.100) THEN 3725 WRITE(6,*) ' Output from DGEEV ' 3726 WRITE(6,*) ' List of eigenvalues, real and imaginary parts' 3727 CALL WRT_2VEC(VLR,VLI,NDIM) 3728 END IF 3729*. test if there are imaginary parts of eigenvalues 3730 THRES = 1.0D-10 3731 XMAX_IMAG = 0.0D0 3732 N_IMAG = 0 3733 DO IVAL = 1, NDIM 3734 IF(ABS(VLI(IVAL)).GT.THRES) N_IMAG = N_IMAG + 1 3735 XMAX_IMAG = MAX(XMAX_IMAG,ABS(VLI(IVAL))) 3736 END DO 3737* 3738 IF(N_IMAG.NE.0) THEN 3739 WRITE(6,*) 3740 & ' GET_LR_EIGVEC_GENMAT Complex eigenvalues encountered ' 3741 WRITE(6,*) ' Number of imaginary eigenvalues ', N_IMAG 3742 WRITE(6,*) ' Largest imaginary component ', XMAX_IMAG 3743 STOP 3744 & ' GET_LR_EIGVEC_GENMAT Complex eigenvalues encountered ' 3745 END IF 3746 K1 = 1 3747 K2 = K1 + NDIM 3748 K3 = K2 + NDIM 3749* 3750 DO I = 1, NDIM 3751 IWORK(K1-1+I) = 0 3752 END DO 3753*. Make sure that <L(I)!R(J)> = delta(I,J) also holds 3754*. degenerate eigenvectors 3755* 3756*. Loop over eigenvalues and collect degenerate sets 3757* 3758 DO I = 1, NDIM 3759*. Ensure that this eigenvalue has not been studied before 3760 IF(IWORK(K1-1+I).EQ.0) THEN 3761*. Check if eigenvalue I is degenerate with other 3762 NDEG = 0 3763 DO J = 1, NDIM 3764 IF(IWORK(K1-1+J).EQ.0.AND.ABS(VLR(I)-VLR(J)).LE.THRES) THEN 3765 NDEG = NDEG + 1 3766 IWORK(K2-1+NDEG) = J 3767 END IF 3768 END DO 3769*. The NDEG eigenvalues IWORK(K2-1+J),J=1,NDEG are degenerate, set 3770*. up overlap matrix in this space 3771 IF(NDEG.EQ.1) THEN 3772*. No problem- just scale L(I) so <L(I)!R(I)> = 1 3773 IWORK(K1-1+I) = 1 3774 XLR = INPROD(VCL(1,I),VCR(1,I),NDIM) 3775 SCALE = 1.0D0/SCALE 3776 CALL SCALVE(VCL(1,I),SCALE,NDIM) 3777 ELSE 3778*. Mark that these vectors have been accessed 3779 DO K = 1, NDEG 3780 IWORK(K1-1+IWORK(K2-1+K)) = 1 3781 END DO 3782*. construct overlap matrix of degenerate vectors 3783 DO II = 1, NDEG 3784 DO J = 1, NDEG 3785 SCRMAT((J-1)*NDEG + II) = 3786 & INPROD(VCL(1,II),VCR(1,J),NDIM) 3787 END DO 3788 END DO 3789 IF(NTEST.GE.100) THEN 3790 WRITE(6,*) 3791 & ' Overlap matrix <L(I)!R(J)> for degenerate eigvectors' 3792 WRITE(6,*) ' Eigenvectors: ', (IWORK(K2-1+II),II=1,NDEG) 3793 CALL WRTMAT(SCRMAT,NDEG,NDEG,NDEG) 3794 END IF 3795*. Find inverse of overlap matrix- use original matrix as scratch space, 3796*. inverse is returned in SCRMAT 3797 CALL INVMAT(SCRMAT,A,NDIM,NDIM,ISING) 3798C INVMAT(A,B,MATDIM,NDIM,ISING) 3799 IF(ISING.NE.0) THEN 3800 WRITE(6,*) ' Problems with matrix inversion' 3801 WRITE(6,*) ' I was programmed by an optimistic person' 3802 WRITE(6,*) ' so I continue' 3803 END IF 3804*.Transform the left eigenvectors, L'(i) = sum_k L_k S^-1 _ik, save in A 3805 DO II = 1, NDEG 3806 DO K = 1, NDEG 3807 WORK(K3-1+K) = SCRMAT((K-1)*NDEG+II) 3808 END DO 3809C MATRIX_TIMES_SPARSEVEC(A,VECELM,IVECIND,NDIM,NVEC,VECOUT) 3810 CALL MATRIX_TIMES_SPARSEVEC(VCL,WORK(K3),IWORK(K2),NDIM, 3811 & NDEG,A(1,II)) 3812 END DO 3813*. And copy back to the matrix of lefteigenvectors 3814 DO II = 1, NDEG 3815 ICOL = IWORK(K2-1+II) 3816 CALL COPVEC(A(1,II),VCL(1,ICOL),NDIM) 3817 END DO 3818 END IF 3819* ^ End if degenerate space contained more than one element 3820 END IF 3821* ^ End if I had not been used before 3822 END DO 3823* ^ End of loop over I 3824* 3825 IF(NTEST.GE.100) THEN 3826 WRITE(6,*) ' Eigenvalues: real and imaginary parts' 3827 CALL WRT_2VEC(VLR,VLI,NDIM) 3828* WRT_2VEC(VEC1,VEC2,NDIM) 3829 WRITE(6,*) ' Space of bioorthonormal L and R eigenvectors' 3830 WRITE(6,*) ' (Real eigenvalues assumed...)' 3831 CALL WRTMAT(VCL,NDIM,NDIM,NDIM,NDIM) 3832 CALL WRTMAT(VCR,NDIM,NDIM,NDIM,NDIM) 3833 END IF 3834* 3835 RETURN 3836 END 3837 SUBROUTINE GET_BLOCK_OF_HS_IN_INTERNAL_SPACE( 3838 & IEXTP,IINTSM,I_HS,HSBLCK,I_INT_TP,I_ONLY_DIA) 3839* 3840*. Obtain block of H0(I_HS=1) or S(I_HS=2) over internal states 3841*. for external type IEXTP and internal symmetry IINTSM. 3842* IF_I_ONLY_DIA = 1, then only the diagonal is calculated 3843* 3844* Jeppe Olsen, March 11, 2009 3845* 3846 INCLUDE 'wrkspc.inc' 3847 INCLUDE 'lucinp.inc' 3848 INCLUDE 'cei.inc' 3849 INCLUDE 'ctcc.inc' 3850*. Output 3851 DIMENSION HSBLCK(*) 3852* 3853 CALL GET_BLOCK_OF_HS_IN_INTERNAL_SPACE_SLAVE( 3854 & IEXTP,IINTSM,I_HS,HSBLCK, 3855 & N_INTOP_TP,WORK(KLSOBEX),N_EXTOP_TP, 3856 & WORK(KL_N_INT_FOR_EXT),WORK(KL_IB_INT_FOR_EXT), 3857 & WORK(KL_I_INT_FOR_EXT), 3858 & WORK(KL_N_INT_FOR_SE),WORK(KL_N_ORTN_FOR_SE), 3859 & I_IN_TP,I_INT_OFF,I_EXT_OFF,I_ONLY_DIA) 3860* 3861 RETURN 3862 END 3863 SUBROUTINE GET_BLOCK_OF_HS_IN_INTERNAL_SPACE_SLAVE( 3864 % IEXTP,IINTSM,I_HS,HSBLCK, 3865 & N_INTP,ISPOBEX,N_EXTP, 3866 & N_IN_FOR_EX,IB_IN_FOR_EX, 3867 & I_IN_FOR_EX, 3868 & N_INT_FOR_SE,N_ORTN_FOR_SE, 3869 & I_IN_TP,IB_INTP,IB_EXTP,I_ONLY_DIA) 3870* 3871*. Obtain block of H (I_HS=1) or S (I_HS=2) for external 3872* type IEXTP and Internal symmetry IINTSM 3873* 3874* If I_ONLY_DIA = 1, then only the diagonal is constructed and stored in 3875* HSBLCK 3876* 3877* X_INT_EI_FOR_SE contains righthandside zero-order states 3878* whereas XL_INT_EI_FOR_SE contains lefthand side zero-order states 3879* 3880*. Jeppe Olsen, October 6, 2008 3881* 3882*. Note current version assumes that expansion of all internal 3883*. states over determinants may be kept in memory.. 3884 INCLUDE 'wrkspc.inc' 3885 REAL*8 INPROD 3886* 3887 INCLUDE 'cgas.inc' 3888 INCLUDE 'gasstr.inc' 3889 INCLUDE 'lucinp.inc' 3890 INCLUDE 'ctcc.inc' 3891 INCLUDE 'cstate.inc' 3892 INCLUDE 'cands.inc' 3893 INCLUDE 'multd2h.inc' 3894 INCLUDE 'glbbas.inc' 3895 INCLUDE 'clunit.inc' 3896 INCLUDE 'oper.inc' 3897 INCLUDE 'cintfo.inc' 3898 INCLUDE 'crun.inc' 3899*. Input: Complete matrices, i.e. matrices for all external types 3900*. and internal symmetries are used. 3901*. Number of internal types for given external type, base for internal 3902*- types for given external types and the actual internal types for 3903* given external type 3904 INTEGER N_IN_FOR_EX(*),IB_IN_FOR_EX(*),I_IN_FOR_EX(*) 3905*. Number of elementary internal operators of given sym 3906*. for given external type 3907 INTEGER ISPOBEX(NGAS,4,*) 3908*. Number of orthonormal internal states per symmetry and external type 3909 INTEGER N_ORTN_FOR_SE(NSMOB,*) 3910*. Number of internal states per symmetry and external types 3911 INTEGER N_INT_FOR_SE(NSMOB,*) 3912*. Output: complete matrix or diagonal 3913 DIMENSION HSBLCK(*) 3914* 3915 NTEST = 0 3916 IF(NTEST.GE.5) THEN 3917 WRITE(6,*) ' ---------------------------------- ' 3918 WRITE(6,*) ' GET_BLOCK_OF_HS_IN_INTERNAL_SPACE ' 3919 WRITE(6,*) ' ---------------------------------- ' 3920* 3921 IF(I_HS.EQ.1) THEN 3922 WRITE(6,*) ' H-block will be constructed' 3923 ELSE 3924 WRITE(6,*) ' S-block wil be constructed' 3925 END IF 3926 WRITE(6,*) ' External type of block: ', IEXTP 3927 WRITE(6,*) ' Symmetry of internal block: ', IINTSM 3928 IF(I_ONLY_DIA.EQ.1) 3929 & WRITE(6,*) ' Only diagonal terms calculated' 3930 WRITE(6,*) ' I_IN_TP, I_INT_HAM = ', I_IN_TP, I_INT_HAM 3931C? WRITE(6,*) ' WORK(KINT1) = ', WORK(KINT1) 3932 END IF 3933* 3934 CALL MEMMAN(IDUM,IDUM,'MARK ',1,'GT_HS ') 3935* 3936 ICSM_SAVE = ICSM 3937 ISSM_SAVE = ISSM 3938 PSSIGN_SAVE = PSSIGN 3939 IDC_SAVE = IDC 3940 PSSIGN = 0.0D0 3941 IDC = 1 3942* 3943* --------------------------------------------------------------- 3944*. Offsets and general info for given external type and internal 3945* symmetry 3946* --------------------------------------------------------------- 3947* 3948*. NOTE : reference space is assumed to be space 1 3949 IREFSPC = 1 3950*. Number of elementary internal states 3951 LINT = N_INT_FOR_SE(IINTSM,IEXTP) 3952*. Number of orthonormal internal states 3953 LINT_ORTN = N_ORTN_FOR_SE(IINTSM,IEXTP) 3954C? WRITE(6,*) ' LINT,LINT_ORTN, ', LINT,LINT_ORTN 3955*. Number of internal types for given external type 3956 L_INTP = N_IN_FOR_EX(IEXTP) 3957*. Offset in I_IN_FOR_EX for given external type 3958 IOFF = IB_IN_FOR_EX(IEXTP) 3959C? WRITE(6,*) ' L_INTP, IOFF = ', L_INTP, IOFF 3960*. The active space 3961 IACT_GAS = 0 3962 DO IGAS = 1, NGAS 3963 IF(IHPVGAS(IGAS).EQ.3) IACT_GAS = IGAS 3964 END DO 3965C? WRITE(6,*) ' The active GASpace ', IACT_GAS 3966*. Symmetry of internal operator times reference state 3967 IINTREFSM = MULTD2H(IINTSM,IREFSM) 3968C? WRITE(6,*) ' IINTSM, IREFSM, IINTREFSM = ', 3969C? & IINTSM, IREFSM, IINTREFSM 3970*. Two vectors of length = number of elementary internal operators 3971 CALL MEMMAN(KLVEC1,LINT,'ADDL ',2,'LVEC1 ') 3972 CALL MEMMAN(KLVEC2,LINT,'ADDL ',2,'LVEC1 ') 3973*. Obtain occupations of alpha- and beta-strings of reference space - 3974*. is currently assumed to be a single space 3975C GET_REF_ALBE_OCC(IREFSPC,IREF_AL,IREF_BEC 3976 CALL MEMMAN(KL_REF_AL,NGAS,'ADDL ',2,'REF_AL') 3977 CALL MEMMAN(KL_REF_BE,NGAS,'ADDL ',2,'REF_BE') 3978 CALL MEMMAN(KL_IREF_AL,NGAS,'ADDL ',2,'IEF_AL') 3979 CALL MEMMAN(KL_IREF_BE,NGAS,'ADDL ',2,'IEF_BE') 3980 CALL GET_REF_ALBE_OCC(IREFSPC,WORK(KL_REF_AL),WORK(KL_REF_BE)) 3981 IF(NTEST.GE.100) THEN 3982 WRITE(6,*) ' Occupation of alpha- and beta-strings in reference' 3983 CALL IWRTMA(WORK(KL_REF_AL),1,NGAS,1,NGAS) 3984 CALL IWRTMA(WORK(KL_REF_BE),1,NGAS,1,NGAS) 3985 END IF 3986 3987*. In the following we are going to use/misuse the standard 3988*. spinorbital types stored as the first elements in the spinorbitalexcitation arrays. 3989*. Save the information usually stored there 3990*. Arrays 3991 CALL MEMMAN(KLSOBEX_SAVE,(NSPOBEX_TP+1)*NGAS*4,'ADDL ',1,'SPOXSV') 3992 CALL ICOPVE(ISPOBEX,WORK(KLSOBEX_SAVE),(NSPOBEX_TP+1)*NGAS*4) 3993 NSPOBEX_TP_SAVE = NSPOBEX_TP 3994*. Prepare the arrays defining calculations for this space 3995 NSPOBEX_TP = L_INTP 3996 DO II_INTP = 1, L_INTP 3997 I_INTP = I_IN_FOR_EX(IOFF-1+II_INTP) 3998C? WRITE(6,*) ' II_INTP, I_INTP ', II_INTP, I_INTP 3999C? WRITE(6,*) ' IB_INTP = ', IB_INTP 4000 CALL ICOPVE(ISPOBEX(1,1,IB_INTP-1+I_INTP), 4001 & ISPOBEX(1,1,II_INTP),4*NGAS) 4002 END DO 4003* Type of C-coefficients will be type of reference 4004 ICATP = 1 4005 ICBTP = 2 4006* 4007 N_ALEL_C = NELFTP(ICATP) 4008 N_BEEL_C = NELFTP(ICBTP) 4009*. Find how action of internal operator changes occupation, 4010*. as all operators makes identical changes of occupations, 4011*. it is sufficient to look on the first operator 4012 IALDEL = IELSUM(ISPOBEX(1,1,1),NGAS)-IELSUM(ISPOBEX(1,3,1),NGAS) 4013 IBEDEL = IELSUM(ISPOBEX(1,2,1),NGAS)-IELSUM(ISPOBEX(1,4,1),NGAS) 4014C? WRITE(6,*) ' IALDEL, IBEDEL = ', IALDEL, IBEDEL 4015*. Occupation of internal operator times reference space 4016 CALL CAAB_T_ABSTR(ISPOBEX(1,1,1),WORK(KL_REF_AL),WORK(KL_REF_BE), 4017 & WORK(KL_IREF_AL),WORK(KL_IREF_BE),NGAS) 4018*. We now have occupation of resulting strings, find corresponding 4019*. supergroups 4020 N_ALEL_S = N_ALEL_C + IALDEL 4021 N_BEEL_S = N_BEEL_C + IBEDEL 4022 IF(NTEST.GE.100) THEN 4023 WRITE(6,*) ' Number of alpha-electrons in IOP*|ref> ',N_ALEL_S 4024 WRITE(6,*) ' Number of beta- electrons in IOP*|ref> ',N_BEEL_S 4025 END IF 4026*. Find supergroup types with these number of electrons 4027 ISATP = 0 4028 ISBTP = 0 4029 DO ISPGP_TP = 1, NSTTP 4030 IF(NELFTP(ISPGP_TP).EQ.N_ALEL_S) ISATP = ISPGP_TP 4031 IF(NELFTP(ISPGP_TP).EQ.N_BEEL_S) ISBTP = ISPGP_TP 4032 END DO 4033 IF(NTEST.GE.100) 4034 &WRITE(6,*) ' a-,b-supergroups of internal operator times ref(I)', 4035 &ISATP,ISBTP 4036*. Below is not neccesary 4037 CALL FIND_SPGRP_FROM_OCC(WORK(KL_IREF_AL),IAL_SPGP_S,ISATP) 4038 CALL FIND_SPGRP_FROM_OCC(WORK(KL_IREF_BE),IBE_SPGP_S,ISBTP) 4039 IF(NTEST.GE.100) THEN 4040 WRITE(6,*) 4041 & ' a-,b-supergroups of internal operator times ref(II)', 4042 & ISATP,ISBTP 4043 END IF 4044*. End of not neccesary 4045* 4046* ------------------------------------------------------------------- 4047* Generate the space with this occupation. This space is 4048* stored as the last CI space and combination space. If the numbers 4049* of these equals MXPICI I am trouble 4050* ------------------------------------------------------------------- 4051* 4052 IF(NCISPC.EQ.MXPICI.OR.NCMBSPC.EQ.MXPICI) THEN 4053 WRITE(6,*) ' Not space for an extra CI space ' 4054 WRITE(6,*) ' Increase parameter MXPICI' 4055 WRITE(6,*) ' NCISPC, MXPICI ', NCISPC, MXPICI 4056 WRITE(6,*) ' NCMBSPC, MXPICI ', NCMBSPC, MXPICI 4057 STOP ' Increase parameter MXPICI' 4058 END IF 4059 IMODREF_CISPC = NCISPC + 1 4060 IMODREF_CMBSPC = NCMBSPC + 1 4061C? WRITE(6,*) ' IMODREF_CISPC, IMODREF_CMBSPC =', 4062C? & IMODREF_CISPC, IMODREF_CMBSPC 4063 LCMBSPC(IMODREF_CMBSPC) = 1 4064 ICMBSPC(1,IMODREF_CMBSPC) = IMODREF_CISPC 4065* 4066 IALTP_FOR_GAS(IMODREF_CMBSPC) = ISATP 4067 IBETP_FOR_GAS(IMODREF_CMBSPC) = ISBTP 4068*. Occupation constraints for modified reference: electrons are 4069*. added or removed from the active space. At the moment a single 4070* active space is assumed. 4071*. Copy reference space occupation 4072 CALL ICOPVE(IGSOCCX(1,1,IREFSPC),IGSOCCX(1,1,IMODREF_CISPC), 4073 & NGAS) 4074 CALL ICOPVE(IGSOCCX(1,2,IREFSPC),IGSOCCX(1,2,IMODREF_CISPC), 4075 & NGAS) 4076*. Change number of electrons in active space 4077 DO IGAS = 1, NGAS 4078 IF(IGAS.GE.IACT_GAS) THEN 4079 IGSOCCX(IGAS,1,IMODREF_CISPC) = 4080 & IGSOCCX(IGAS,1,IMODREF_CISPC) + IALDEL + IBEDEL 4081 IGSOCCX(IGAS,2,IMODREF_CISPC) = 4082 & IGSOCCX(IGAS,2,IMODREF_CISPC) + IALDEL + IBEDEL 4083 END IF 4084 END DO 4085* 4086 IF(NTEST.GE.1000) THEN 4087 WRITE(6,*) ' Min and max for modified reference space ' 4088 CALL IWRTMA(IGSOCCX(1,1,IMODREF_CISPC),1,NGAS,1,NGAS) 4089 CALL IWRTMA(IGSOCCX(1,2,IMODREF_CISPC),1,NGAS,1,NGAS) 4090 END IF 4091*. Number of determinants in modified reference space 4092 LDET = LEN_CISPC(IMODREF_CMBSPC,IINTREFSM,NTEST) 4093C WRITE(6,*) ' LDET = ', LDET 4094*. Space for expansion Int |ref> in SD for all Internal states 4095 IDUM = 0 4096 IF(I_ONLY_DIA.EQ.0) THEN 4097 NTERMS = LINT_ORTN 4098 ELSE 4099 NTERMS = LDET 4100 END IF 4101 LDIM = LINT_ORTN*NTERMS 4102 CALL MEMMAN(KL_INTSDMAT, LDIM,'ADDL ',2,'INTMAT') 4103 CALL MEMMAN(KL_INTSDMAT2,LDIM,'ADDL ',2,'INTMAT') 4104* 4105 ICSPC = IREFSPC 4106 ISSPC = IMODREF_CISPC 4107 ICSM = IREFSM 4108 ISSM = IINTREFSM 4109 IF(I_HS.EQ.2) THEN 4110 IF(I_ONLY_DIA.EQ.0) THEN 4111*. Obtain complete block 4112*. O+(I,right)|ref> in KL_INTSDMAT 4113 DO INTOP = 1, LINT_ORTN 4114 CALL REWINO(LUC) 4115 CALL REWINO(LUHC) 4116*. obtain internal operator INTOP 4117C GET_ZERO_ORDER_STATE(IEXTP,INTSM,IORTN,X,ILR) 4118 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4119 ICSM = IREFSM 4120 ISSM = IINTREFSM 4121 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4122 & WORK(KLVEC1),1) 4123 CALL REWINO(LUHC) 4124 CALL FRMDSCN(WORK(KL_INTSDMAT+(INTOP-1)*LDET),-1,-1,LUHC) 4125 IF(NTEST.GE.1000) THEN 4126 WRITE(6,*) ' O+(I,right) |ref> in SD-basis for I = ', INTOP 4127 CALL WRTMAT(WORK(KL_INTSDMAT+(INTOP-1)*LDET),1,LDET,1,LDET) 4128 END IF 4129 END DO 4130*. O+(J.left)|ref> in KL_INTSDMAT2 4131 DO INTOP = 1, LINT_ORTN 4132 CALL REWINO(LUC) 4133 CALL REWINO(LUHC) 4134*. obtain internal operator INTOP 4135C GET_ZERO_ORDER_STATE(IEXTP,INTSM,IORTN,X,ILR) 4136 ICSM = IREFSM 4137 ISSM = IINTREFSM 4138 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4139 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4140 & WORK(KLVEC1),1) 4141 CALL REWINO(LUHC) 4142 CALL FRMDSCN(WORK(KL_INTSDMAT2+(INTOP-1)*LDET),-1,-1,LUHC) 4143 IF(NTEST.GE.1000) THEN 4144 WRITE(6,*) ' O+(I,left) |ref> in SD-basis for I = ', INTOP 4145 CALL WRTMAT(WORK(KL_INTSDMAT2+(INTOP-1)*LDET),1,LDET,1,LDET) 4146 END IF 4147 END DO 4148*. S(I,J) = <ref!op(i,left) op(j,right)!ref> 4149 DO I = 1, LINT_ORTN 4150 DO J = 1, LINT_ORTN 4151 HSBLCK((J-1)*LINT_ORTN + I) = 4152 & INPROD(WORK(KL_INTSDMAT2+(I-1)*LDET), 4153 & WORK(KL_INTSDMAT+(J-1)*LDET),LDET) 4154 END DO 4155 END DO 4156* 4157 IF(NTEST.GE.1000) THEN 4158 WRITE(6,*) ' The overlap <ref!op(i)+ op(j)!ref>' 4159 CALL PRSYM(HSBLCK,LINT) 4160 END IF 4161 ELSE IF(I_ONLY_DIA.EQ.1) THEN 4162*. Obtain only diagonal block 4163 DO INTOP = 1, LINT_ORTN 4164 CALL REWINO(LUC) 4165 CALL REWINO(LUHC) 4166 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4167 ICSM = IREFSM 4168 ISSM = IINTREFSM 4169 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4170 & WORK(KLVEC1),1) 4171 CALL REWINO(LUHC) 4172 CALL FRMDSCN(WORK(KL_INTSDMAT),-1,-1,LUHC) 4173 IF(NTEST.GE.1000) THEN 4174 WRITE(6,*) ' O+(I,right) |ref> in SD-basis for I = ', INTOP 4175 CALL WRTMAT(WORK(KL_INTSDMAT),1,LDET,1,LDET) 4176 END IF 4177*. O+(I.left)|ref> in KL_INTSDMAT2 4178 CALL REWINO(LUC) 4179 CALL REWINO(LUHC) 4180 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4181 ICSM = IREFSM 4182 ISSM = IINTREFSM 4183 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4184 & WORK(KLVEC1),1) 4185 CALL REWINO(LUHC) 4186 CALL FRMDSCN(WORK(KL_INTSDMAT2),-1,-1,LUHC) 4187 IF(NTEST.GE.1000) THEN 4188 WRITE(6,*) ' O+(I,left) |ref> in SD-basis for I = ', INTOP 4189 CALL WRTMAT(WORK(KL_INTSDMAT2),1,LDET,1,LDET) 4190 END IF 4191*. S(I,I) = <ref!op(i,left) op(i,right)!ref> 4192 HSBLCK(INTOP) = 4193 & INPROD(WORK(KL_INTSDMAT2),WORK(KL_INTSDMAT),LDET) 4194 END DO 4195* 4196 IF(NTEST.GE.1000) THEN 4197 WRITE(6,*) ' The overlapdiagonal <ref!op(i)+ op(i)!ref>' 4198 CALL WRTMAT(HSBLCK,LINT_ORTN) 4199 END IF 4200 END IF 4201* ^ End if I_ONLY_DIA = 1 4202 ELSE IF(I_HS.EQ.1) THEN 4203 IF(I_ONLY_DIA.EQ.0) THEN 4204 IF(I_IN_TP.GE.2) THEN 4205*. Obtain zero-order Hamiltonian in internal space 4206C. H0 Intop(right) |ref> for all orthonormal states 4207 DO INTOP = 1, LINT_ORTN 4208 CALL REWINO(LUHC) 4209 ICSPC = IREFSPC 4210 ISSPC = IMODREF_CMBSPC 4211 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4212 ICSM = IREFSM 4213 ISSM = IINTREFSM 4214 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4215 & WORK(KLVEC1),1) 4216*. zero-order Hamiltonian times Int(intop,right)|ref> in INTSDMAT 4217 ICSPC = IMODREF_CMBSPC 4218 ISSPC = IMODREF_CMBSPC 4219 ISSM = IINTREFSM 4220 ICSM = IINTREFSM 4221 IF(I_INT_HAM.EQ.1) THEN 4222C? WRITE(6,*) 4223C? & ' One-body H0 used for internal zero-order states' 4224 I12 = 1 4225 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4226 ELSE 4227C? WRITE(6,*) 4228C? & ' Two-body H used for internal zero-order states' 4229 END IF 4230 CALL REWINO(LUHC) 4231 CALL REWINO(LUSC51) 4232 ICSM = IINTREFSM 4233 ISSM = IINTREFSM 4234 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51,0,0) 4235 IF(NTEST.GE.1000) THEN 4236 WRITE(6,*) ' After mv7 ' 4237 WRITE(6,*) ' Resulting vector : ' 4238 CALL WRTVCD(WORK(KVEC1P),LUSC51,1,-1) 4239 END IF 4240 IF(I_INT_HAM.EQ.1) 4241 & CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4242* 4243 CALL REWINO(LUSC51) 4244 CALL FRMDSCN(WORK(KL_INTSDMAT+(INTOP-1)*LDET),-1,-1, 4245 & LUSC51) 4246 END DO 4247C. Intop(left) |ref> for all orthonormal states IN INTSDMAT2 4248 DO INTOP = 1, LINT_ORTN 4249 ICSPC = IREFSPC 4250 ISSPC = IMODREF_CMBSPC 4251 ICSM = IREFSM 4252 ISSM = IINTREFSM 4253C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4254 CALL REWINO(LUC) 4255 CALL REWINO(LUHC) 4256 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4257 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4258 & WORK(KLVEC1),1) 4259 CALL REWINO(LUHC) 4260 CALL FRMDSCN(WORK(KL_INTSDMAT2+(INTOP-1)*LDET),-1,-1,LUHC) 4261 END DO 4262*. <ref|intop(left) H0 intop(right)|ref> 4263 DO IORTN = 1, LINT_ORTN 4264 DO JORTN = 1, LINT_ORTN 4265 IJ = (JORTN-1)*LINT_ORTN + IORTN 4266*. H0 4267 HSBLCK(IJ) = 4268 & INPROD(WORK(KL_INTSDMAT+(JORTN-1)*LDET), 4269 & WORK(KL_INTSDMAT2+(IORTN-1)*LDET),LDET) 4270 END DO 4271 END DO 4272* 4273 IF(NTEST.GE.1000) THEN 4274 WRITE(6,*) ' H0 over orthonormal states ' 4275 CALL WRTMAT(HSBLCK,LINT_ORTN,LINT_ORTN,LINT_ORTN,LINT_ORTN) 4276 END IF 4277 END IF 4278 IF(I_IN_TP.EQ.3) THEN 4279*. <0!O+I [H0,OJ]|0> should be obtained 4280*. On input HSBLCK contains <0!O+I H0 OJ|0>. 4281*. Obtain <0!O+I OJ H0 |0> 4282*. H0 |0> on LUHC 4283 ICSPC = IREFSPC 4284 ISSPC = IREFSPC 4285 ICSM = IREFSM 4286 ISSM = IREFSM 4287 IF(I_INT_HAM.EQ.1) THEN 4288C? WRITE(6,*) 4289C? & ' One-body H0 used for internal zero-order states' 4290 I12 = 1 4291 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4292 ELSE 4293C? WRITE(6,*) 4294C? & ' Two-body H used for internal zero-order states' 4295 END IF 4296 CALL REWINO(LUC) 4297 CALL REWINO(LUHC) 4298 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC,0,0) 4299 IF(I_INT_HAM.EQ.1) CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4300*. Generate OI H0 |ref> and save in KL_INTSDMAT 4301 DO INTOP= 1, LINT_ORTN 4302 ICSPC = IREFSPC 4303 ISSPC = IMODREF_CMBSPC 4304C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4305 CALL REWINO(LUSC51) 4306 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4307 ICSM = IINTREFSM 4308 ISSM = IREFSM 4309 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51, 4310 & WORK(KLVEC1),1) 4311 CALL FRMDSCN(WORK(KL_INTSDMAT+(INTOP-1)*LDET),-1,-1,LUSC51) 4312 END DO 4313*. Generate OI |ref> in KL_INTSDMAT2 4314 DO INTOP = 1, LINT_ORTN 4315 ICSPC = IREFSPC 4316 ISSPC = IMODREF_CMBSPC 4317C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4318 CALL REWINO(LUHC) 4319 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4320 ICSM = IREFSM 4321 ISSM = IINTREFSM 4322 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4323 & WORK(KLVEC1),1) 4324 CALL REWINO(LUSC51) 4325 CALL FRMDSCN(WORK(KL_INTSDMAT2+(INTOP-1)*LDET),-1,-1, 4326 & LUSC51) 4327 END DO 4328* Obtain <ref!o+i [H0,o j]|ref> = <ref!o+i H0 o j|ref> - 4329* <ref!o+i o j H0|ref> 4330 DO IORTN = 1, LINT_ORTN 4331 DO JORTN = 1, LINT_ORTN 4332 IJ = (JORTN-1)*LINT_ORTN + IORTN 4333*. H0 4334 HSBLCK(IJ) = HSBLCK(IJ) - 4335 & INPROD(WORK(KL_INTSDMAT2+(JORTN-1)*LDET), 4336 & WORK(KL_INTSDMAT+(IORTN-1)*LDET),LDET) 4337 END DO 4338 END DO 4339* 4340 IF(NTEST.GE.1000) THEN 4341 WRITE(6,*) 4342 & ' <ref!o+i [H0,o j]|ref> over orthonormal states' 4343 CALL WRTMAT(HSBLCK,LINT_ORTN,LINT_ORTN,LINT_ORTN,LINT_ORTN) 4344 END IF 4345* 4346 END IF 4347*. ^ End if I_IN_TP = 3 4348 ELSE IF(I_ONLY_DIA.EQ.1) THEN 4349 IF(I_IN_TP.GE.2) THEN 4350*. Obtain diagonal of zero-order Hamiltonian in internal space 4351 DO INTOP = 1, LINT_ORTN 4352C. H0 Intop(right) |ref> 4353 CALL REWINO(LUHC) 4354 ICSPC = IREFSPC 4355 ISSPC = IMODREF_CMBSPC 4356 ICSM = IREFSM 4357 ISSM = IINTREFSM 4358 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4359 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4360 & WORK(KLVEC1),1) 4361*. zero-order Hamiltonian times Int(intop,right)|ref> in INTSDMAT 4362 IF(I_INT_HAM.EQ.1) THEN 4363C? WRITE(6,*) 4364C? & ' One-body H0 used for internal zero-order states' 4365 I12 = 1 4366 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4367 ELSE 4368C? WRITE(6,*) 4369C? & ' Two-body H used for internal zero-order states' 4370 END IF 4371 ICSPC = IMODREF_CMBSPC 4372 ISSPC = IMODREF_CMBSPC 4373 CALL REWINO(LUHC) 4374 CALL REWINO(LUSC51) 4375 ICSM = IINTREFSM 4376 ISSM = IINTREFSM 4377 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51,0,0) 4378 IF(I_INT_HAM.EQ.1) 4379 & CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4380 IF(NTEST.GE.1000) THEN 4381 WRITE(6,*) ' After mv7 ' 4382 WRITE(6,*) ' Resulting vector : ' 4383 CALL WRTVCD(WORK(KVEC1P),LUSC51,1,-1) 4384 END IF 4385 CALL REWINO(LUSC51) 4386 CALL FRMDSCN(WORK(KL_INTSDMAT),-1,-1, 4387 & LUSC51) 4388C. Intop(left) |ref> INTSDMAT2 4389 ICSPC = IREFSPC 4390 ISSPC = IMODREF_CMBSPC 4391 ICSM = IREFSM 4392 ISSM = IINTREFSM 4393 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4394C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4395 CALL REWINO(LUC) 4396 CALL REWINO(LUHC) 4397 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4398 & WORK(KLVEC1),1) 4399 CALL REWINO(LUHC) 4400 CALL FRMDSCN(WORK(KL_INTSDMAT2),-1,-1,LUHC) 4401*. <ref|intop(left) H0 intop(right)|ref> 4402 HSBLCK(INTOP) = 4403 & INPROD(WORK(KL_INTSDMAT),WORK(KL_INTSDMAT2),LDET) 4404 END DO 4405* 4406 IF(NTEST.GE.1000) THEN 4407 WRITE(6,*) ' Diagonal of H0 over orthonormal states ' 4408 CALL WRTMAT(HSBLCK,LINT_ORTN,1,LINT_ORTN,1) 4409 END IF 4410 END IF 4411 IF(I_IN_TP.EQ.3) THEN 4412*. Not prepared for I_INT_HAM = 2 !!! 4413*. <0!O+I [H0,OI]|0> should be obtained 4414*. On input HSBLCK contains <0!O+I H0 OI|0>. 4415*. Obtain <0!O+I OI H0 |0> 4416*. H0 |0> on LUHC 4417 ICSPC = IREFSPC 4418 ISSPC = IREFSPC 4419 ICSM = IREFSM 4420 ISSM = IREFSM 4421 I12 = 1 4422 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4423 CALL REWINO(LUC) 4424 CALL REWINO(LUHC) 4425 CALL MV7(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC,0,0) 4426 CALL SWAPVE(WORK(KINT1),WORK(KFIFA),NINT1) 4427*. Generate OI H0 |ref> and save in KL_INTSDMAT 4428 DO INTOP= 1, LINT_ORTN 4429C? WRITE(6,*) ' INTOP= ', INTOP 4430 ICSPC = IREFSPC 4431 ISSPC = IMODREF_CMBSPC 4432 ICSM = IREFSM 4433 ISSM = IINTREFSM 4434C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4435 CALL REWINO(LUSC51) 4436 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),2) 4437 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUHC,LUSC51, 4438 & WORK(KLVEC1),1) 4439 CALL FRMDSCN(WORK(KL_INTSDMAT),-1,-1,LUSC51) 4440*. Generate OI |ref> in KL_INTSDMAT2 4441 ICSPC = IREFSPC 4442 ISSPC = IMODREF_CMBSPC 4443C SIGDEN_CC(C,HC,LUC,LUHC,T,ISIGDEN) 4444 CALL REWINO(LUHC) 4445 CALL GET_ZERO_ORDER_STATE(IEXTP,IINTSM,INTOP,WORK(KLVEC1),1) 4446 CALL SIGDEN_CC(WORK(KVEC1P),WORK(KVEC2P),LUC,LUHC, 4447 & WORK(KLVEC1),1) 4448 CALL REWINO(LUSC51) 4449 CALL FRMDSCN(WORK(KL_INTSDMAT2),-1,-1, 4450 & LUSC51) 4451* Obtain <ref!o+i [H0,o j]|ref> = <ref!o+i H0 o j|ref> - 4452 HSBLCK(INTOP) = HSBLCK(INTOP) - 4453 & INPROD(WORK(KL_INTSDMAT2),WORK(KL_INTSDMAT),LDET) 4454 END DO 4455* 4456 IF(NTEST.GE.1000) THEN 4457 WRITE(6,*) 4458 & ' <ref!o+i [H0,o i]|ref> over orthonormal states' 4459 CALL WRTMAT(HSBLCK,1,LINT_ORTN,1,LINT_ORTN) 4460 END IF 4461* 4462 END IF 4463*. ^ End if I_IN_TP = 3 4464 END IF 4465*. ^ End of I_ONLY_DIA switch 4466 END IF 4467* ^ 4468* End if I_HS = 1 4469*. Clean up: restore spinorbitaltypes and symmetries 4470 NSPOBEX_TP = NSPOBEX_TP_SAVE 4471 CALL ICOPVE(WORK(KLSOBEX_SAVE),ISPOBEX,(NSPOBEX_TP+1)*NGAS*4) 4472* 4473 ICSM = ICSM_SAVE 4474 ISSM = ISSM_SAVE 4475 PSSIGN = PSSIGN_SAVE 4476 IDC = IDC_SAVE 4477* 4478 CALL MEMMAN(IDUM,IDUM,'FLUSM',1,'GT_HS ') 4479* 4480 RETURN 4481 END 4482 FUNCTION E1_FOR_STRING(HDIAG,ISTRING,NEL) 4483* 4484* Obtain one-electron contribution to energy from given string 4485* 4486*. Jeppe Olsen, March 2009 4487* 4488 INCLUDE 'implicit.inc' 4489* 4490*. Input: 4491 DIMENSION ISTRING(NEL) 4492 DIMENSION HDIAG(*) 4493* 4494 E1 = 0.0D0 4495 DO IEL = 1, NEL 4496 E1 = E1 + HDIAG(ISTRING(IEL)) 4497 END DO 4498* 4499 E1_FOR_STRING = E1 4500* 4501 RETURN 4502 END 4503 FUNCTION GET_E0REF_EXT(FI,IPHGASX) 4504* 4505*. Obtain external part of zero-order energy. 4506*. State is assumed to contain a double occupied hole part 4507*. and an unoccupied particle part. P/H is flagged by IPHGASX) 4508* 4509*. Jeppe Olsen, March 11, 2009 4510* 4511 INCLUDE 'wrkspc.inc' 4512 INCLUDE 'cgas.inc' 4513 INCLUDE 'orbinp.inc' 4514*. Specific input 4515 DIMENSION FI(*) 4516 DIMENSION IPHGASX(*) 4517* 4518 NTEST = 100 4519 WRITE(6,*) ' I am not working' 4520 STOP ' GET_E0REF_EXT not working' 4521* 4522 E0REF_EXT = 0.0D0 4523 DO IGAS = 1, NGAS 4524 IF(IGAS.EQ.1) THEN 4525 IB = 1 4526 ELSE 4527 IB = IB + NOBPT(IGAS) 4528 END IF 4529 IF(IPHGASX(IGAS).EQ.2) THEN 4530 DO I = 1, NOBPT(IGAS) 4531*. Absolute number of orbital in ST ordering 4532 I_ABS = IREOTS(IB-1+I) 4533 E0REF_EXT = E0REF_EXT + 2.0D0*FI((I_ABS+1)*I_ABS/2) 4534 END DO 4535 END IF 4536 END DO 4537* 4538 GET_E0REF_EXT = E0REF_EXT 4539* 4540 IF(NTEST.GE.100) THEN 4541 WRITE(6,*) ' External part of zero-order energy ', E0REF_EXT 4542 END IF 4543* 4544 RETURN 4545 END 4546 FUNCTION N_ZERO_ORDER_STATES(NORTN_FOR_SE,NDIM_EX_ST,N_EXTP, 4547 & ITOTSM) 4548* 4549*. Number of zero-order states of given symmetry 4550* 4551*. Jeppe Olsen, March 2009 4552* 4553 INCLUDE 'wrkspc.inc' 4554*. General input 4555 INCLUDE 'multd2h.inc' 4556 INCLUDE 'lucinp.inc' 4557*. Specific input 4558 INTEGER NORTN_FOR_SE(NSMOB,N_EXTP),NDIM_EX_ST(NSMOB,N_EXTP) 4559* 4560 NTEST = 1000 4561*. 4562 N = 0 4563 DO I_EXTP = 1, N_EXTP 4564 DO I_EXSM = 1, NSMOB 4565 I_INSM = MULTD2H(I_EXSM,ITOTSM) 4566 N = N + NORTN_FOR_SE(I_INSM,I_EXTP)*NDIM_EX_ST(I_EXSM,I_EXTP) 4567C? WRITE(6,*) ' I_EXTP, I_EXSM, I_INSM = ',I_EXTP,I_EXSM,I_INSM 4568C? WRITE(6,*) ' NORTN_FOR_SE, NDIM_EX_ST = ', 4569C? & NORTN_FOR_SE(I_INSM,I_EXTP),NDIM_EX_ST(I_EXSM,I_EXTP) 4570 END DO 4571 END DO 4572* 4573 N_ZERO_ORDER_STATES = N 4574* 4575 IF(NTEST.GE.100) THEN 4576 WRITE(6,*) ' Number of zero-order states = ', 4577 & N_ZERO_ORDER_STATES 4578 END IF 4579* 4580 RETURN 4581 END 4582 SUBROUTINE GET_ZERO_ORDER_STATE(IEXTP,INTSM,IORTN,X,ILR) 4583* 4584* Obtain internal state IORTN of symmetri INTSM and corresponding to 4585* external type IEXTP 4586* 4587*. Master code 4588* 4589 INCLUDE 'wrkspc.inc' 4590 INCLUDE 'lucinp.inc' 4591 INCLUDE 'cei.inc' 4592*.Output 4593 DIMENSION X(*) 4594* 4595 IDUM = 0 4596 CALL MEMMAN(IDUM,IDUM,'MARK ', IDUM,'GT_ZOST') 4597 CALL MEMMAN(KLSCR,N_INT_MAX,'ADDL ',2,'SCRINT') 4598* 4599 IF(ILR.EQ.2) THEN 4600 CALL GET_ZERO_ORDER_STATE_SLAVE(IEXTP,INTSM,IORTN,X, 4601 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 4602 & WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_X2_INT_EI_FOR_SE), 4603 & WORK(KL_SG_INT_EI_FOR_SE), 4604 & WORK(KL_IBX1_INT_EI_FOR_SE),WORK(KL_IBX2_INT_EI_FOR_SE), 4605 & WORK(KL_IBSG_INT_EI_FOR_SE),WORK(KLSCR), 4606 & NSMOB) 4607 ELSE 4608 CALL GET_ZERO_ORDER_STATE_SLAVE(IEXTP,INTSM,IORTN,X, 4609 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 4610 & WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_X2L_INT_EI_FOR_SE), 4611 & WORK(KL_SG_INT_EI_FOR_SE), 4612 & WORK(KL_IBX1_INT_EI_FOR_SE),WORK(KL_IBX2_INT_EI_FOR_SE), 4613 & WORK(KL_IBSG_INT_EI_FOR_SE),WORK(KLSCR), 4614 & NSMOB) 4615 END IF 4616* 4617 CALL MEMMAN(IDUM,IDUM,'FLUSM', IDUM,'GT_ZOST') 4618* 4619 RETURN 4620 END 4621 SUBROUTINE GET_ZERO_ORDER_STATE_SLAVE(IEXTP,INTSM,IORTN,X, 4622 & N_ORTN_FOR_SE,N_INT_FOR_SE, 4623 & X1,X2,SG,IBX1,IBX2,IBSG,XSCR,NSMOB) 4624* 4625*. Obtain zero-order state IORTN of symmetry INTSM corresponding 4626* to external type IEXTP 4627* 4628*. Jeppe Olsen, March 12, 2009 4629* 4630*. General Input 4631 include 'implicit.inc' 4632* 4633 INTEGER N_ORTN_FOR_SE(NSMOB,*),N_INT_FOR_SE(NSMOB,*) 4634 DIMENSION X1(*),X2(*),SG(*) 4635 INTEGER IBX1(NSMOB,*),IBX2(NSMOB,*),IBSG(NSMOB,*) 4636*. Scratch 4637 DIMENSION XSCR(*) 4638*. Output 4639 DIMENSION X(*) 4640* 4641 NTEST = 00 4642 IF(NTEST.GE.100) THEN 4643 WRITE(6,*) ' GET_ZERO_ORDER_STATE_SLAVE in action ' 4644 WRITE(6,*) ' IEXTP,INTSM = ', IEXTP,INTSM 4645 END IF 4646* 4647* Zero-order state X(IORTN)(J) =(X1 SG(-1/2( X2)) K,IORTN 4648* 4649 IB_X1 = IBX1(INTSM,IEXTP) 4650 IB_X2 = IBX2(INTSM,IEXTP) 4651 IB_SG = IBSG(INTSM,IEXTP) 4652* 4653 LINT = N_INT_FOR_SE(INTSM,IEXTP) 4654 LORTN = N_ORTN_FOR_SE(INTSM,IEXTP) 4655 IF(NTEST.GE.1000) THEN 4656 WRITE(6,*) ' LINT, LORTN = ',LINT,LORTN 4657 WRITE(6,*) ' IB_X1, IB_X2, IB_SG = ', IB_X1,IB_X2,IB_SG 4658 WRITE(6,*) ' Blocks of Sigma, X1, X2:' 4659 CALL WRTMAT(SG(IB_SG),1,LORTN,1,LORTN) 4660 CALL WRTMAT(X1(IB_X1),LINT,LORTN,LINT,LORTN) 4661 CALL WRTMAT(X2(IB_X2),LORTN,LORTN,LORTN,LORTN) 4662 END IF 4663* SG(-1/2) X(2) K, IORTN 4664 DO K = 1, LORTN 4665 XSCR(K) = 4666 & 1.0D0/SQRT(SG(IB_SG-1+K))*X2(IB_X2-1+(IORTN-1)*LORTN + K) 4667 END DO 4668* (X1 SG(-1/2( X2)) K,IORTN 4669C MATVCC(A,VIN,VOUT,NROW,NCOL,ITRNS) 4670 CALL MATVCC(X1(IB_X1),XSCR,X,LINT,LORTN,0) 4671* 4672 IF(NTEST.GE.100) THEN 4673 WRITE(6,*) ' Internal orthonormal zero-order state ' 4674 WRITE(6,*) ' IEXTP, INTSM, IORTN =', IEXTP,INTSM,IORTN 4675 WRITE(6,*) ' Expansion in elementary operators ' 4676 CALL WRTMAT(X,LINT,1,LINT,1) 4677 END IF 4678* 4679 RETURN 4680 END 4681 SUBROUTINE TRANS_EI_ORTN(T_EI,T_ORTN,ITSYM,IEO,ILR,ICOCON) 4682* 4683* Transform between elementary ei and orthonormal form of 4684* vector 4685* 4686* IEO = 1 => elementary ei to orthonormal order 4687* IEO = 2 => orthonormal to elementary ei order 4688* 4689* ICOCON = 1 => Covariant transformation (transform as eigenvectors) 4690* ICOCON = 2 => Contravariant transformation (transform to ensure 4691* invariance of scalar) 4692* 4693* Zero-order state is explicitly transferred as last elements in 4694* T_EI and T_ORTN, respectively (if I_INCLUDE_UNI=1) 4695* 4696* Jeppe Olsen, July 2009, finished on the train to Dusseldorf, aug14 4697* 4698 INCLUDE 'wrkspc.inc' 4699 INCLUDE 'cei.inc' 4700 INCLUDE 'lucinp.inc' 4701* 4702*. Explicit input and output 4703* 4704 DIMENSION T_EI(*), T_ORTN(*) 4705* 4706 IF(IEO.EQ.2.AND.ICOCON.EQ.1) THEN 4707 WRITE(6,*) 4708 & ' Covariant backtransformation to elementary basis not defined' 4709 STOP 4710 & ' Covariant backtransformation to elementary basis not defined' 4711 END IF 4712* 4713 CALL TRANS_EI_ORTN_SLAVE(T_EI,T_ORTN,ITSYM,IEO,ILR,N_EXTOP_TP, 4714 & WORK(KL_NDIM_EX_ST),WORK(KL_NDIM_IN_SE), 4715 & WORK(KL_N_ORTN_FOR_SE),NSMOB,ICOCON,I_INCLUDE_UNI) 4716* 4717 RETURN 4718 END 4719 SUBROUTINE TRANS_EI_ORTN_SLAVE(T_EI,T_ORTN,ITSYM,IEO,ILR,N_EXTP, 4720 & NDIM_EX_ST,NDIM_IN_ST,NDIM_ORT_ST,NSMOB,ICOCON, 4721 & I_INCLUDE_UNI) 4722 INCLUDE 'implicit.inc' 4723*. General input 4724 INCLUDE 'multd2h.inc' 4725 DIMENSION NDIM_EX_ST(NSMOB,N_EXTP),NDIM_IN_ST(NSMOB,N_EXTP), 4726 & NDIM_ORT_ST(NSMOB,N_EXTP) 4727*. Specific input and output 4728 DIMENSION T_EI(*), T_ORTN(*) 4729* 4730 NTEST = 00 4731* 4732 IOFF_EI=1 4733 IOFF_ORT=1 4734 DO I_EXTP = 1, N_EXTP 4735 DO I_EXSM = 1, NSMOB 4736 I_INSM = MULTD2H(I_EXSM,ITSYM) 4737* 4738 N_EX = NDIM_EX_ST(I_EXSM,I_EXTP) 4739 N_IN = NDIM_IN_ST(I_INSM,I_EXTP) 4740 N_ORT = NDIM_ORT_ST(I_INSM,I_EXTP) 4741* 4742 CALL TRANS_EI_ORTN_BL(T_EI(IOFF_EI),T_ORTN(IOFF_ORT), 4743 & I_EXTP,I_INSM,N_EX,IEO,ILR,ICOCON) 4744 IOFF_EI = IOFF_EI + N_EX*N_IN 4745 IOFF_ORT = IOFF_ORT + N_EX*N_ORT 4746 END DO 4747 END DO 4748* 4749 IF(I_INCLUDE_UNI.EQ.1) THEN 4750 IF(IEO.EQ.1) THEN 4751 T_ORTN(IOFF_ORT) = T_EI(IOFF_EI) 4752 ELSE 4753 T_EI(IOFF_EI) = T_ORTN(IOFF_ORT) 4754 END IF 4755 END IF 4756* 4757 IF(NTEST.GE.100) THEN 4758* 4759 IF(IEO.EQ.1) THEN 4760 WRITE(6,*) ' Transformation: T(I,E) => T(ORT,E)' 4761 ELSE 4762 WRITE(6,*) ' Transformation: T(ORT,E) => T(I,E)' 4763 END IF 4764* 4765 IF(ICOCON.EQ.1) THEN 4766 WRITE(6,*) ' Covariant transformation ' 4767 ELSE 4768 WRITE(6,*) ' Contravariant transformation' 4769 END IF 4770* 4771 WRITE(6,*) ' Vector T(I,E):' 4772 CALL PRINT_T_EI(T_EI,1,ITSYM) 4773 WRITE(6,*) ' Vector T(Ort,E):' 4774 CALL PRINT_T_EI(T_ORTN,2,ITSYM) 4775* 4776COLD IF(I_INCLUDE_UNI.EQ.1) THEN 4777COLD WRITE(6,*) ' And the coefficient for unitop:' 4778COLD IF(IEO.EQ.1) THEN 4779COLD WRITE(6,*) T_EI(IOFF_EI) 4780COLD ELSE 4781COLD WRITE(6,*) T_ORTN(IOFF_ORT) 4782COLD END IF 4783COLD END IF 4784 END IF 4785* 4786 RETURN 4787 END 4788 SUBROUTINE TRANS_EI_ORTN_BL(T_EI,T_ORTN,IEXTP,INTSM,NVEC,IEO,ILR, 4789 & ICOCON) 4790* 4791*. A block of coefficients given by external type IEXTP and 4792*. internal symmetry INTSM 4793*. Transform between orthonormal form and elementary ei form of 4794*. operators. 4795* IEO = 1 => elementary ei to orthonormal order 4796* IEO = 2 => orthonormal to elementaty ei order 4797* 4798* ICOCON = 1 => Covariant transformation (transform as eigenvectors) 4799* ICOCON = 2 => Contravariant transformation (transform to ensure 4800* invariance) 4801* 4802* The elements are supposed to be in EI order 4803* 4804*. Jeppe Olsen, March 12, 2009, on the train to Koeln, 4805* July 29, 2009, ICOCON added 4806* 4807* Last modification; Oct. 18, 2012; Jeppe Olsen, reduced allocation 4808* 4809 INCLUDE 'wrkspc.inc' 4810 INCLUDE 'lucinp.inc' 4811 INCLUDE 'cei.inc' 4812*. Specific input and output 4813 DIMENSION T_EI(*), T_ORTN(*) 4814* 4815 CALL MEMMAN(IDUM,IDUM,'MARK ', IDUM,'TREIOR') 4816 CALL MEMMAN(KLSCR,N_INT_MAX*NVEC,'ADDL ',2,'SCRINT') 4817* 4818 IF(ILR.EQ.2) THEN 4819 CALL TRANS_EI_ORTN_EL_SLAVE 4820 & (T_EI,T_ORTN,IEXTP,INTSM,NVEC,IEO, 4821 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 4822 & WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_X2_INT_EI_FOR_SE), 4823 & WORK(KL_SG_INT_EI_FOR_SE), 4824 & WORK(KL_IBX1_INT_EI_FOR_SE),WORK(KL_IBX2_INT_EI_FOR_SE), 4825 & WORK(KL_IBSG_INT_EI_FOR_SE),WORK(KLSCR), 4826 & NSMOB,ICOCON) 4827 ELSE 4828 CALL TRANS_EI_ORTN_EL_SLAVE 4829 & (T_EI,T_ORTN,IEXTP,INTSM,NVEC,IEO, 4830 & WORK(KL_N_ORTN_FOR_SE),WORK(KL_N_INT_FOR_SE), 4831 & WORK(KL_X1_INT_EI_FOR_SE),WORK(KL_X2L_INT_EI_FOR_SE), 4832 & WORK(KL_SG_INT_EI_FOR_SE), 4833 & WORK(KL_IBX1_INT_EI_FOR_SE),WORK(KL_IBX2_INT_EI_FOR_SE), 4834 & WORK(KL_IBSG_INT_EI_FOR_SE),WORK(KLSCR), 4835 & NSMOB,ICOCON) 4836 END IF 4837* 4838 CALL MEMMAN(IDUM,IDUM,'FLUSM ', IDUM,'TREIOR') 4839* 4840 RETURN 4841 END 4842 SUBROUTINE TRANS_EI_ORTN_EL_SLAVE( 4843 & T_EI,T_ORTN,IEXTP,INTSM,NVEC,IEO, 4844 & N_ORTN_FOR_SE,N_INT_FOR_SE, 4845 & X1,X2,SG,IBX1,IBX2,IBSG,XSCR,NSMOB,ICOCON) 4846* 4847* Transform between elementary and orthonormal form of 4848* a block of coefficients given by external type IEXTP and 4849* internal symmetry INTSM 4850* 4851* 4852* IEO = 1 => elementary to orthonormal order 4853* IEO = 2 => orthonormal to elementaty order 4854* 4855* ICOCON = 1 => Covariant transformation (transform as eigenvectors) 4856* ICOCON = 2 => Contravariant transformation (transform to ensure 4857* invariance) 4858* 4859*. Jeppe Olsen, March 12, 2009 4860* 4861*. Contravariant transformation: 4862* 4863* T_ORTN = X(2)T sigma(1/2) X(1)T T_EI 4864* T_EI = X(1) Sigma(-1/2) X(2) T_ORTN 4865* 4866*. Covariant transformation: 4867* 4868* V_ORTN = X(2)T sigma(-1/2)X1(T) V_EI 4869* Transformation from V_ORTN to V_EI is not defined or needed (I hope) 4870 4871 4872 4873 4874 4875*. Generel Input 4876 include 'implicit.inc' 4877* 4878 INTEGER N_ORTN_FOR_SE(NSMOB,*),N_INT_FOR_SE(NSMOB,*) 4879 DIMENSION X1(*),X2(*),SG(*) 4880 INTEGER IBX1(NSMOB,*),IBX2(NSMOB,*),IBSG(NSMOB,*) 4881*. Scratch- should hold a block of coefficients 4882 DIMENSION XSCR(*) 4883*. Input and Output (Depending on IEO) 4884 DIMENSION T_EI(*), T_ORTN(*) 4885* 4886 NTEST = 00 4887 IF(NTEST.GE.100) THEN 4888 WRITE(6,*) ' TRANS_EI_ORTN_EL_SLAVE in action' 4889 WRITE(6,*) ' IEXTP,INTSM,IEO,ICOCON =', IEXTP,INTSM,IEO,ICOCON 4890 END IF 4891* 4892 IF(IEO.EQ.2.AND.ICOCON.EQ.1) THEN 4893 WRITE(6,*) 4894 & ' Covariant backtransformation to elementary basis not defined' 4895 STOP 4896 & ' Covariant backtransformation to elementary basis not defined' 4897 END IF 4898* 4899 IB_X1 = IBX1(INTSM,IEXTP) 4900 IB_X2 = IBX2(INTSM,IEXTP) 4901 IB_SG = IBSG(INTSM,IEXTP) 4902* 4903 LINT = N_INT_FOR_SE(INTSM,IEXTP) 4904 LORTN = N_ORTN_FOR_SE(INTSM,IEXTP) 4905 IF(NTEST.GE.1000) THEN 4906 WRITE(6,*) ' LINT, LORTN = ',LINT,LORTN 4907 WRITE(6,*) ' IB_X1, IB_X2, IB_SG = ', IB_X1,IB_X2,IB_SG 4908 END IF 4909* 4910 FACTORC = 0.0D0 4911 FACTORAB = 1.0D0 4912* 4913 IF(IEO.EQ.1) THEN 4914*. Elementary => orthonormal transformation 4915* T_ORTN = X(2)T sigma(1/2) X(1)T T_EI 4916*. X(1)T T_EI 4917 CALL MATML7(XSCR,X1(IB_X1),T_EI,LORTN,NVEC,LINT,LORTN,LINT,NVEC, 4918 & FACTORC, FACTORAB,1) 4919* Sigma(+/- 1/2) (X(1)T T_EI) - done here columnwise to reduce number of 4920* sqrts - probable set up an array of sigma(-1/2) and proceed rowwise 4921 DO IROW = 1, LORTN 4922 IF(ICOCON.EQ.1) THEN 4923 FACTOR = 1.0D0/SQRT(SG(IB_SG-1+IROW)) 4924 ELSE 4925 FACTOR = SQRT(SG(IB_SG-1+IROW)) 4926 END IF 4927 DO ICOL = 1, NVEC 4928 XSCR((ICOL-1)*LORTN+IROW) = FACTOR*XSCR((ICOL-1)*LORTN+IROW) 4929 END DO 4930 END DO 4931* X(2)T (sigma(+/- 1/2) X(1)T T_EI) = X(2)T XSCR 4932 CALL MATML7(T_ORTN,X2(IB_X2),XSCR, 4933 & LORTN,NVEC,LORTN,LORTN,LORTN,NVEC, 4934 & FACTORC,FACTORAB,1) 4935 ELSE IF(IEO.EQ.2) THEN 4936*. Orthonormal to elementary order 4937* T_EI = X(1) Sigma(-1/2) X(2) T_ORTN 4938* 4939*. X(2) T_ORTN in XSCR 4940 CALL MATML7(XSCR,X2(IB_X2),T_ORTN, 4941 & LORTN,NVEC,LORTN,LORTN,LORTN,NVEC, 4942 & FACTORC,FACTORAB,0) 4943* Sigma(-1/2) X(2) T_ORTN = Sigma(-1/2) XSCR 4944 DO IROW = 1, LORTN 4945 FACTOR = 1.0D0/SQRT(SG(IB_SG-1+IROW)) 4946 DO ICOL = 1, NVEC 4947 XSCR((ICOL-1)*LORTN+IROW) = FACTOR*XSCR((ICOL-1)*LORTN+IROW) 4948 END DO 4949 END DO 4950* X(1) Sigma(-1/2) X(2) T_ORTN = X(1) XSCR in T_EI 4951 CALL MATML7(T_EI,X1(IB_X1),XSCR, 4952 & LINT,NVEC,LINT,LORTN,LORTN,NVEC, 4953 & FACTORC, FACTORAB) 4954 END IF 4955* ^End of IEO switch 4956* 4957 IF(NTEST.GE.100) THEN 4958 WRITE(6,*) 4959 & ' Transformation between orthonormal and internal expansion' 4960 WRITE(6,*) ' IEXTP, INTSM =', IEXTP,INTSM 4961 IF(IEO.EQ.1) THEN 4962 WRITE(6,*) ' Elementary => orthonormal transformation' 4963 ELSE 4964 WRITE(6,*) ' Orthonormal=> elementary transformation' 4965 END IF 4966 IF(ICOCON.EQ.1) THEN 4967 WRITE(6,*) ' Convariant transformation ' 4968 ELSE 4969 WRITE(6,*) ' Contravariant transformation ' 4970 END IF 4971 WRITE(6,*) ' Coefficients in elementary basis' 4972 CALL WRTMAT(T_EI,LINT,NVEC,LINT,NVEC) 4973 WRITE(6,*) ' Coefficients in orthonormal basis' 4974 CALL WRTMAT(T_ORTN,LORTN,NVEC,LORTN,NVEC) 4975 END IF 4976* 4977 RETURN 4978 END 4979 FUNCTION LARGEST_BLOCK_IN_MAT(NBLK,LR,LC) 4980* 4981* A matrix is given with NBLK blocks with row dim LR anc column dim LC 4982* Find largest block 4983* 4984* Jeppe Olsen, March 12, 2009 4985 INCLUDE 'implicit.inc' 4986*. Input 4987 INTEGER LR(*), LC(*) 4988* 4989 LMAX = 0 4990 DO I = 1, NBLK 4991 LMAX = MAX(LMAX,LR(I)*LC(I)) 4992 END DO 4993* 4994 NTEST = 100 4995 IF(NTEST.GE.100) THEN 4996 WRITE(6,*) ' Largest block = ', LMAX 4997 END IF 4998* 4999 LARGEST_BLOCK_IN_MAT = LMAX 5000* 5001 RETURN 5002 END 5003 SUBROUTINE COMPARE_TO_UNI(A,NDIM) 5004* 5005* A matrix A is given. Find largest deviation from unit matrix 5006* Jeppe Olsen, Aug. 2009, Red Roof Hotel, Washington 5007* 5008 INCLUDE 'implicit.inc' 5009 DIMENSION A(NDIM,NDIM) 5010* 5011 IOFF_R = 0 5012 IOFF_C = 0 5013 I_DIAG = 0 5014 XMAX_OFF = 0.0D0 5015 XMAX_DIAG = 0.0D0 5016 DO I = 1, NDIM 5017 DO J = 1,NDIM 5018 IF(I.NE.J) THEN 5019 IF(ABS(A(I,J)).GT.XMAX_OFF) THEN 5020 XMAX_OFF = ABS(A(I,J)) 5021 IOFF_R = I 5022 IOFF_C = J 5023 END IF 5024 ELSE 5025 IF(ABS(A(I,I)-1.0D0).GT.XMAX_DIAG) THEN 5026 XMAX_DIAG = ABS(A(I,I)-1) 5027 I_DIAG = I 5028 END IF 5029 END IF 5030*. ^ End of diagonal/off diagonal switch 5031 END DO 5032 END DO 5033* 5034 WRITE(6,*) ' Comparison of matrix to unit matrix ' 5035 WRITE(6,*) ' Largest offdiagonal element : value, row, column =', 5036 & XMAX_OFF,IOFF_R, IOFF_C 5037 WRITE(6,*) 5038 &' Largest deviation of unit element from 1: value and row ', 5039 & XMAX_DIAG,I_DIAG 5040* 5041 RETURN 5042 END 5043 SUBROUTINE NORM_T_EI(T,IEO,ITSYM,XNORM_EI,IPRT) 5044* 5045* Norm of the various blocks of a T(I,E) vector given in elementary(IEO=1) or 5046* orthonormal(IEO=2) form 5047* 5048*. Jeppe Olsen, Nov. 12, 2009 5049* 5050 INCLUDE 'wrkspc.inc' 5051 INCLUDE 'cei.inc' 5052 INCLUDE 'lucinp.inc' 5053*. Specific input 5054 DIMENSION T(*) 5055*. Output 5056 DIMENSION XNORM_EI(*) 5057* 5058 IF(IEO.EQ.1) THEN 5059 CALL NORM_T_EI_SLAVE(T,ITSYM,N_EXTOP_TP, 5060 & WORK(KL_NDIM_EX_ST),WORK(KL_NDIM_IN_SE),NSMOB,XNORM_EI) 5061 ELSE 5062 CALL NORM_T_EI_SLAVE(T,ITSYM,N_EXTOP_TP, 5063 & WORK(KL_NDIM_EX_ST),WORK(KL_N_ORTN_FOR_SE),NSMOB,XNORM_EI) 5064 END IF 5065* 5066 IF(IPRT.NE.0) THEN 5067 WRITE(6,*) ' Norm of T-EI vector for various E-types' 5068 CALL WRTMAT(XNORM_EI,1,N_EXTOP_TP,1,N_EXTOP_TP) 5069 END IF 5070* 5071 RETURN 5072 END 5073 SUBROUTINE NORM_T_EI_SLAVE(T,ITSYM,N_EXTP, 5074 & NDIM_EX_ST,NDIM_IN_ST,NSMOB,XNORM_EI) 5075* 5076 INCLUDE 'implicit.inc' 5077 INCLUDE 'multd2h.inc' 5078* 5079 REAL*8 5080 &INPROD 5081*. Specific input 5082 DIMENSION T(*) 5083 DIMENSION NDIM_EX_ST(NSMOB,N_EXTP),NDIM_IN_ST(NSMOB,N_EXTP) 5084*. Output 5085 DIMENSION XNORM_EI(*) 5086* 5087 IOFF = 1 5088 DO I_EXTP = 1, N_EXTP 5089 X = 0.0D0 5090 DO I_EXSM = 1, NSMOB 5091 I_INSM = MULTD2H(I_EXSM,ITSYM) 5092 N_EX = NDIM_EX_ST(I_EXSM,I_EXTP) 5093 N_IN = NDIM_IN_ST(I_INSM,I_EXTP) 5094C? WRITE(6,*) ' N_EX, N_IN = ', N_EX, N_IN 5095 X = X +INPROD(T(IOFF),T(IOFF),N_EX*N_IN) 5096 IOFF = IOFF + N_EX*N_IN 5097 END DO 5098 XNORM_EI(I_EXTP) = SQRT(X) 5099 END DO 5100* 5101 RETURN 5102 END 5103 SUBROUTINE TEST_E 5104* 5105* Test calc of E 5106* 5107 INCLUDE 'wrkspc.inc' 5108 INCLUDE 'clunit.inc' 5109 INCLUDE 'cecore.inc' 5110 INCLUDE 'glbbas.inc' 5111 INCLUDE 'cgas.inc' 5112 REAL*8 INPRDD 5113* 5114 DIMENSION XJ1(10000),XJ2(10000) 5115* 5116 WRITE(6,*) ' First element of INT1 = ', WORK(KINT1) 5117 WRITE(6,*) ' IPHGAS: ' 5118 CALL IWRTMA(IPHGAS,1,NGAS,1,NGAS) 5119 CALL MV7(XJ1,XJ2,LUC,LUHC,0,0) 5120* 5121 EREF = INPRDD(XJ1,XJ2,LUC,LUHC,1,-1) 5122 WRITE(6,*) ' EREF, ECORE calc in ....', EREF+ECORE, ECORE 5123* 5124 RETURN 5125 END 5126 5127 5128 5129c $Id$ 5130