1 SUBROUTINE GET_SUPSYM_INFO 2* 3* Obtain supersymmetry info 4* 5*. Jeppe Olsen, May 2012 6* Revised July 2012 7* Last revision; Jeppe Olsen; March 5, 2013; call to GET_SUPSYM_FOR_BASIS changed to include 8* irrep info 9* 10* Last modified July 8, 2012 (Jeppe) 11* 12 INCLUDE 'implicit.inc' 13 INCLUDE 'mxpdim.inc' 14 INCLUDE 'wrkspc-static.inc' 15 INCLUDE 'glbbas.inc' 16 INCLUDE 'orbinp.inc' 17 INCLUDE 'lucinp.inc' 18 INCLUDE 'crun.inc' 19 INCLUDE 'cgas.inc' 20#include "errquit.fh" 21#include "mafdecls.fh" 22#include "global.fh" 23* 24 IDUM = 0 25 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'GTSPSM') 26*. Parity for standard symmetry 27 IF(CSUPSYM.EQ.'LINEAR'.AND.INVCNT.EQ.1) THEN 28 CALL SET_PARITY_FOR_STASYM(IPA_FOR_STASYM,NIRREP) 29 ELSE 30 IZERO = 0 31 CALL ISETVC(IPA_FOR_STASYM,IZERO,NIRREP) 32 END IF 33* Labels for the basis functions 34 CALL GET_SUPSYM_LABELS_FOR_ORBITALS 35 &(int_mb(KCSUPSYM_FOR_ORB),int_mb(KLVAL_FOR_ORB), 36 & int_mb(KMLVAL_FOR_ORB),int_mb(KPA_FOR_ORB)) 37*. Relation between supersymmetry and irreps 38 CALL SYM_AND_IRREP_FOR_SUPSYM( 39 & int_mb(KL_FOR_SUPSYM),int_mb(KML_FOR_SUPSYM), 40 & int_mb(KPA_FOR_SUPSYM),int_mb(KIRREP_FOR_SUPSYM), 41 & int_mb(KNSUPSYM_FOR_IRREP), int_mb(KIBSUPSYM_FOR_IRREP), 42 & int_mb(KISUPSYM_FOR_IRREP) ) 43C SYM_AND_IRREP_FOR_SUPSYM 44C & (L_FOR_SUPSYM,ML_FOR_SUPSYM,IPA_FOR_SUPSYM,IRREP_FOR_SUPSYM, 45C & NSUPSYM_FOR_IRREP, IBSUPSYM_FOR_IRREP,ISUPSYM_FOR_IRREP) 46* 47*. Info on the supersymmetry of the basis set 48* 49 CALL GET_SUPSYM_FOR_BASIS(int_mb(KISUPSYM_FOR_BAS), 50 & int_mb(KNBAS_FOR_SUP_STA_SYM),int_mb(KIBBAS_FOR_SUP_STA_SYM), 51 & int_mb(KIBAS_FOR_SUP_STA_SYM),int_mb(KISHELL_FOR_BAS), 52 & int_mb(KNBAS_FOR_SHELL),int_mb(KIBBAS_FOR_SHELL), 53 & int_mb(KIBAS_FOR_SHELL), 54 & int_mb(KNSUPSYM_FOR_IRREP),int_mb(KIBSUPSYM_FOR_IRREP), 55 & int_mb(KISUPSYM_FOR_IRREP)) 56* 57* Obtain mappings from standard to super symmetry components 58* simply by reading info in NBAS_FOR_SUP_STA_SYM 59 CALL SUP_TO_STASYM(int_mb(KNBAS_FOR_SUP_STA_SYM)) 60*. Number, offsets and actual numbers of orbitals of 61*. given symmetry 62* 63* The GENSMOB arrays allowing a symmetric treatment of standard 64* and super symmetry 65*. Just let the supersymmetry components be general symmetry 66 NGENSMOB= N_SUPSYM 67 CALL ICOPVE(NBAS_SUPSYM,NBAS_GENSMOB,NGENSMOB) 68 CALL ICOPVE(IBBAS_SUPSYM,IBBAS_GENSMOB,NGENSMOB) 69 CALL ICOPVE(IBAS_SUPSYM,ISTA_TO_GENSM_REO,NTOOB) 70 CALL ICOPVE(ISTASM_FOR_SUPSYM,ISTASM_FOR_GENSM,NGENSMOB) 71* 72 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'GTSPSM') 73 RETURN 74 END 75 SUBROUTINE GET_MAX_SUPSYM_IRREP 76* 77* Obtain Max supersymmetry Irrep from read of labels of orbitals 78*. (saved in ORBINP) 79* 80* Jeppe Olsen, May 2012 81* 82 INCLUDE 'implicit.inc' 83 INCLUDE 'mxpdim.inc' 84 INCLUDE 'glbbas.inc' 85 INCLUDE 'orbinp.inc' 86 INCLUDE 'crun.inc' 87* 88*. General input 89 CHARACTER*4 AO_CENT, AO_TYPE, CHAR4 90 COMMON/AOLABELS/AO_CENT(MXPORB),AO_TYPE(MXPORB) 91*. A label is in general of the form nlm, where l is the 92*. standard spectroscopic term 93* 94 NTEST = 100 95 IF(NTEST.GE.100) THEN 96 WRITE(6,*) 97 WRITE(6,*) ' ==============================' 98 WRITE(6,*) ' Info from GET_MAX_SUPSYM_IRREP' 99 WRITE(6,*) ' ==============================' 100 WRITE(6,*) 101 END IF 102*. Obtained largest L-value 103 LMAX = 0 104 DO IORB = 1, NTOOB 105 CHAR4 = AO_TYPE(IORB) 106 IF(NTEST.GE.1000) WRITE(6,'(A,A)') ' AO_TYPE = ', AO_TYPE(IORB) 107 IF(CHAR4(2:2).EQ.'s') THEN 108 LMAX = MAX(0,LMAX) 109 ELSE IF(CHAR4(2:2).EQ.'p') THEN 110 LMAX = MAX(1,LMAX) 111 ELSE IF(CHAR4(2:2).EQ.'d') THEN 112 LMAX = MAX(2,LMAX) 113 ELSE IF(CHAR4(2:2).EQ.'f') THEN 114 LMAX = MAX(3,LMAX) 115 ELSE IF(CHAR4(2:2).EQ.'g') THEN 116 LMAX = MAX(4,LMAX) 117 ELSE IF(CHAR4(2:2).EQ.'h') THEN 118 LMAX = MAX(5,LMAX) 119 ELSE IF(CHAR4(2:2).EQ.'i') THEN 120 LMAX = MAX(6,LMAX) 121 ELSE 122 WRITE(6,*) ' Unknown form of orbital ', CHAR4 123 STOP ' Unknown form of orbital ' 124 END IF 125 END DO 126* 127C? IF(NTEST.GE.100) WRITE(6,*) ' Lmax = ', LMAX 128* 129 IF(CSUPSYM.EQ.'ATOMIC') THEN 130 N_SUPSYM = LMAX + 1 + LMAX*(LMAX+1) 131 N_SUPSYM_IRREP = LMAX +1 132 ELSE IF (CSUPSYM.EQ.'LINEAR') THEN 133 N_SUPSYM_IRREP = (INVCNT+1)*(LMAX + 1) 134 N_SUPSYM = (INVCNT+1)*(2*LMAX + 1) 135 END IF 136 LMAX_ORB = LMAX 137* 138 IF(NTEST.GE.1) THEN 139 WRITE(6,*) ' N_SUPSYM_IRREP = ', N_SUPSYM_IRREP 140 WRITE(6,*) ' N_SUPSYM = ', N_SUPSYM 141 WRITE(6,*) ' Lmax = ', LMAX 142 WRITE(6,*) ' INVCNT = ', INVCNT 143 END IF 144* 145 RETURN 146 END 147 SUBROUTINE GET_SUPSYM_LABELS_FOR_ORBITALS 148 & (CSUPSYM_FOR_ORB,LVAL_FOR_ORB,MLVAL_FOR_ORB, 149 & IPA_FOR_ORB) 150* 151* Obtain labels for symmetry for the various basis functions 152* 153*. Jeppe Olsen, early morning May 23, 2012 (when I should do other things...) 154* 155 INCLUDE 'implicit.inc' 156 INCLUDE 'mxpdim.inc' 157 INCLUDE 'crun.inc' 158 INCLUDE 'orbinp.inc' 159*. Output 160 CHARACTER*3 CSUPSYM_FOR_ORB(NTOOB) 161 INTEGER LVAL_FOR_ORB(NTOOB) 162 INTEGER MLVAL_FOR_ORB(NTOOB) 163 INTEGER IPA_FOR_ORB(NTOOB) 164* 165* 166 CHARACTER*4 AO_CENT, AO_TYPE, CHAR4, CHAR4B 167 COMMON/AOLABELS/AO_CENT(MXPORB),AO_TYPE(MXPORB) 168* 169 NTEST = 00 170 IF(NTEST.GE.100) THEN 171 WRITE(6,*) 172 WRITE(6,*) ' ========================================' 173 WRITE(6,*) ' Info from GET_SUPSYM_LABELS_FOR_ORBITALS' 174 WRITE(6,*) ' ========================================' 175 WRITE(6,*) 176 END IF 177* 178 DO IORB = 1, NTOOB 179 CHAR4 = AO_TYPE(IORB) 180 IF(CSUPSYM.EQ.'ATOMIC') THEN 181 CSUPSYM_FOR_ORB(IORB) = CHAR4(1:2) 182 ELSE IF (CSUPSYM.EQ.'LINEAR') THEN 183 IF(CHAR4(1:2).EQ.'1s') THEN 184 CSUPSYM_FOR_ORB(IORB) = '0 ' 185 CHAR4B = '0 ' 186 ELSE 187 CHAR4B = ' ' 188 CHAR4B(1:2) = CHAR4(3:4) 189 CSUPSYM_FOR_ORB(IORB) = CHAR4B(1:2) 190 END IF 191*. And parity labels 192 IF(INVCNT.EQ.1) THEN 193 IF(IPA_FOR_STASYM(ISMFSO(IORB)).EQ.1.) THEN 194 CHAR4B(3:3) = 'g' 195 IPA_FOR_ORB(IORB) = 1 196 ELSE 197 CHAR4B(3:3) = 'u' 198 IPA_FOR_ORB(IORB) =-1 199 ENDIF 200 CSUPSYM_FOR_ORB(IORB) = CHAR4B(1:3) 201 END IF !parity is present 202 END IF 203 CALL ORB_LABEL_TO_LML 204 & (AO_TYPE(IORB),LVAL_FOR_ORB(IORB),MLVAL_FOR_ORB(IORB)) 205 END DO 206* 207 IF(NTEST.GE.100) THEN 208 WRITE(6,*) ' Labels for basis functions ' 209 WRITE(6,*) ' ===========================' 210 DO IORB = 1, NTOOB 211 WRITE(6,'(2X,I4,2X,A3)') IORB, CSUPSYM_FOR_ORB(IORB) 212 END DO 213* 214 IF(INVCNT.EQ.0) THEN 215 WRITE(6,*) ' L and ML-values for orbitals ' 216 WRITE(6,*) ' =============================' 217 WRITE(6,*) 218 WRITE(6,*) ' Orbital L Ml ' 219 WRITE(6,*) ' ===============' 220 DO IORB = 1, NTOOB 221 WRITE(6,'(2X,I3,4X,I2,2X,I2)') 222 & IORB,LVAL_FOR_ORB(IORB),MLVAL_FOR_ORB(IORB) 223 END DO 224 ELSE 225 WRITE(6,*) ' L, ML, and parity for orbitals ' 226 WRITE(6,*) ' =============================' 227 WRITE(6,*) 228 WRITE(6,*) ' Orbital L Ml Parity' 229 WRITE(6,*) ' =====================' 230 DO IORB = 1, NTOOB 231 WRITE(6,'(2X,I3,4X,I2,2X,I2,3X,I3)') 232 & IORB,LVAL_FOR_ORB(IORB),MLVAL_FOR_ORB(IORB), 233 & IPA_FOR_ORB(IORB) 234 END DO 235 ENDIF! explicit center of inversion 236* 237 END IF 238* 239 RETURN 240 END 241 SUBROUTINE ORB_LABEL_TO_LML(ORB_LABEL,LVAL,MLVAL) 242* 243* A char*4 lavel ORB_LABEL is given, obtain 244* corresponding values of L and ML 245* 246* Very unelegant routine, but gets the work done (I hope) 247* 248* Jeppe Olsen, May 2012 249* 250*. Input 251 CHARACTER*4 ORB_LABEL 252* 253 NTEST = 00 254* 255*. Lvalue 256* 257 IF(ORB_LABEL(2:2).EQ.'s') THEN 258 LVAL = 0 259 ELSE IF( ORB_LABEL(2:2).EQ.'p') THEN 260 LVAL = 1 261 ELSE IF( ORB_LABEL(2:2).EQ.'d') THEN 262 LVAL = 2 263 ELSE IF(ORB_LABEL(2:2).EQ.'f') THEN 264 LVAL = 3 265 ELSE IF(ORB_LABEL(2:2).EQ.'g') THEN 266 LVAL = 4 267 ELSE IF(ORB_LABEL(2:2).EQ.'h') THEN 268 LVAL = 5 269 ELSE IF(ORB_LABEL(2:2).EQ.'i') THEN 270 LVAL = 6 271 ELSE 272 WRITE(6,*) ' Unknown ORB_LABEL: ', ORB_LABEL 273 STOP ' Unknown ORB_LABEL ' 274 END IF 275* 276*. Ml-value 277* 278 IF(ORB_LABEL(2:2).EQ.'s') THEN 279 MLVAL = 0 280 ELSE IF(ORB_LABEL(2:3).EQ.'px') THEN 281 MLVAL = 1 282 ELSE IF(ORB_LABEL(2:3).EQ.'py') THEN 283 MLVAL = -1 284 ELSE IF(ORB_LABEL(2:3).EQ.'pz') THEN 285 MLVAL = 0 286 ELSE IF(ORB_LABEL(3:3).EQ.'0') THEN 287 MLVAL = 0 288 ELSE IF(ORB_LABEL(3:4).EQ.'1+') THEN 289 MLVAL = 1 290 ELSE IF(ORB_LABEL(3:4).EQ.'1-') THEN 291 MLVAL = -1 292 ELSE IF(ORB_LABEL(3:4).EQ.'2+') THEN 293 MLVAL = 2 294 ELSE IF(ORB_LABEL(3:4).EQ.'2-') THEN 295 MLVAL = -2 296 ELSE IF(ORB_LABEL(3:4).EQ.'3+') THEN 297 MLVAL = 3 298 ELSE IF(ORB_LABEL(3:4).EQ.'3-') THEN 299 MLVAL = -3 300 ELSE IF(ORB_LABEL(3:4).EQ.'4+') THEN 301 MLVAL = 4 302 ELSE IF(ORB_LABEL(3:4).EQ.'4-') THEN 303 MLVAL = -4 304 ELSE IF(ORB_LABEL(3:4).EQ.'5+') THEN 305 MLVAL = 5 306 ELSE IF(ORB_LABEL(3:4).EQ.'5-') THEN 307 MLVAL = -5 308 ELSE IF(ORB_LABEL(3:4).EQ.'6+') THEN 309 MLVAL = 6 310 ELSE IF(ORB_LABEL(3:4).EQ.'6-') THEN 311 MLVAL = -6 312 ELSE 313 WRITE(6,*) ' Unknown ORB_LABEL: ', ORB_LABEL 314 STOP ' Unknown ORB_LABEL: ' 315 END IF 316* 317 IF(NTEST.GE.100) THEN 318 WRITE(6,'(A,A,4X,I2,2X,I2)') 319 & ' ORB_LABEL, LVAL, MLVAL ', ORB_LABEL,LVAL,MLVAL 320 END IF 321* 322 RETURN 323 END 324 SUBROUTINE SYM_AND_IRREP_FOR_SUPSYM( 325 & L_FOR_SUPSYM,ML_FOR_SUPSYM, 326 & IPA_FOR_SUPSYM,IRREP_FOR_SUPSYM, 327 & NSUPSYM_FOR_IRREP, IBSUPSYM_FOR_IRREP, 328 & ISUPSYM_FOR_IRREP) 329*. Obtain tables L_FOR_SUPSYM,ML_FOR_SUPSYM,IPA_FOR_SUPSYM 330* giving L, ML and parity for symmetry 331*. for linear molecules, L is set to zero 332* Obtain IRREP_FOR_SUPSYM giving irrep for symmetry 333* Obtain NSUPSYM_FOR_IRREP, IBSUPSYM_FOR_IRREP,ISUPSYM_FOR_IRREP giving symmetry for irrep 334* 335* Irrep: well irrep, say s,p,d for ATOMIC case, !Lambda! for LINEAR, 336* |Lambda| g for linear with inversion 337* symmetry: given component of irrep, say p+1 for atomic, Pi+1 for linear 338* 339*. Jeppe Olsen, May 23, 2012, inversion added June11 340* 341 INCLUDE 'implicit.inc' 342 INCLUDE 'mxpdim.inc' 343 INCLUDE 'crun.inc' 344 INCLUDE 'orbinp.inc' 345*. Output 346 INTEGER L_FOR_SUPSYM(*), ML_FOR_SUPSYM(*),IPA_FOR_SUPSYM(*) 347 INTEGER IRREP_FOR_SUPSYM(*) 348 INTEGER NSUPSYM_FOR_IRREP(*), IBSUPSYM_FOR_IRREP(*) 349 INTEGER ISUPSYM_FOR_IRREP(*) 350* 351 NTEST = 100 352 IF(NTEST.GE.100) THEN 353 WRITE(6,*) 354 WRITE(6,*) ' ===================================' 355 WRITE(6,*) ' Info from SYM_AND_IRREP_FOR_SUPSYM ' 356 WRITE(6,*) ' ===================================' 357 WRITE(6,*) 358 WRITE(6,*) ' LMAX_ORB = ', LMAX_ORB 359 END IF 360* 361 IF(CSUPSYM.EQ.'ATOMIC') THEN 362 ISYM = 0 363 IRREP = 0 364 DO L = 0, LMAX_ORB 365 IRREP = IRREP + 1 366 NDEG = 2*L + 1 367 NSUPSYM_FOR_IRREP(IRREP) = NDEG 368 IBSUPSYM_FOR_IRREP(IRREP) = ISYM + 1 369 DO ML = -L,L 370 ISYM = ISYM + 1 371 L_FOR_SUPSYM(ISYM) = L 372 ML_FOR_SUPSYM(ISYM) = ML 373 ISUPSYM_FOR_IRREP(ISYM) = ISYM 374 IRREP_FOR_SUPSYM(ISYM) = IRREP 375 END DO 376 END DO 377 ELSE IF (CSUPSYM.EQ.'LINEAR') THEN 378 ISYM = 0 379 IRREP = 0 380 IF(INVCNT.EQ.0) THEN 381 NPARITY = 1 382 ELSE 383 NPARITY = 2 384 END IF 385 DO LAMBDA = 0, LMAX_ORB 386 DO IPARITY = 1, NPARITY 387 IRREP = IRREP + 1 388 IF(LAMBDA.EQ.0) THEN 389 NDEG = 1 390 ELSE 391 NDEG = 2 392 END IF 393 NSUPSYM_FOR_IRREP(IRREP) = NDEG 394 IBSUPSYM_FOR_IRREP(IRREP) = ISYM + 1 395 DO ICOMP = 1, NDEG 396 ISYM = ISYM + 1 397 IRREP_FOR_SUPSYM(ISYM) = IRREP 398 ISUPSYM_FOR_IRREP(ISYM) = ISYM 399 L_FOR_SUPSYM(ISYM) = 0 400 IF(NDEG.EQ.1) THEN 401 ML_FOR_SUPSYM(ISYM) = 0 402 ELSE 403 IF(ICOMP.EQ.1) THEN 404 ML_FOR_SUPSYM(ISYM) = -LAMBDA 405 ELSE 406 ML_FOR_SUPSYM(ISYM) = LAMBDA 407 END IF 408 END IF 409 IF(INVCNT.EQ.0) THEN 410 IPA_FOR_SUPSYM(ISYM) = 0 411 ELSE 412 IF(IPARITY.EQ.1) THEN 413 IPA_FOR_SUPSYM(ISYM) = 1 414 ELSE 415 IPA_FOR_SUPSYM(ISYM) =-1 416 END IF 417 END IF! Parity present 418 END DO ! loop over ICOMP 419 END DO! Loop over parity 420 END DO !loop over lambda 421 END IF! Switch of SUPSYM 422* 423 IF(NTEST.GE.100) THEN 424 IF(CSUPSYM.EQ.'ATOMIC') THEN 425 WRITE(6,*) ' L and ML values for symmetries ' 426 WRITE(6,*) ' ============================== ' 427 WRITE(6,*) 428 WRITE(6,*) ' Symmetry L Ml ' 429 WRITE(6,*) ' =====================' 430 DO ISYM = 1, N_SUPSYM 431 WRITE(6,'(4X,I2,6X,I2,4X,I2)') 432 & ISYM, L_FOR_SUPSYM(ISYM), ML_FOR_SUPSYM(ISYM) 433 END DO 434 ELSE IF(CSUPSYM.EQ.'LINEAR') THEN 435 IF(INVCNT.EQ.0) THEN 436 WRITE(6,*) ' ML values for symmetries ' 437 WRITE(6,*) ' ============================== ' 438 WRITE(6,*) 439 WRITE(6,*) ' Symmetry Ml ' 440 WRITE(6,*) ' =====================' 441 DO ISYM = 1, N_SUPSYM 442 WRITE(6,'(3X,I2,7X,I2)') 443 & ISYM, ML_FOR_SUPSYM(ISYM) 444 END DO 445 ELSE 446 WRITE(6,*) ' Ml and parity for symmetries ' 447 WRITE(6,*) ' ============================== ' 448 WRITE(6,*) 449 WRITE(6,*) ' Symmetry Ml, parity ' 450 WRITE(6,*) ' =====================' 451 DO ISYM = 1, N_SUPSYM 452 WRITE(6,'(3X,I2,7X,I2,5X,I2)') 453 & ISYM, ML_FOR_SUPSYM(ISYM),IPA_FOR_SUPSYM(ISYM) 454 END DO 455 END IF ! parity is considered 456 END IF ! form of supersymmetry 457* 458 WRITE(6,*) ' Irrep for the various symmetries ' 459 WRITE(6,*) ' ================================ ' 460 WRITE(6,*) 461 WRITE(6,*) ' Symmetry Irrep ' 462 WRITE(6,*) ' ================= ' 463 DO ISYM = 1, N_SUPSYM 464 WRITE(6,'(5X,I2,7X,I2)') 465 & ISYM, IRREP_FOR_SUPSYM(ISYM) 466 END DO 467* 468 WRITE(6,*) ' Symmetry for the various irreps ' 469 WRITE(6,*) ' ================================' 470 WRITE(6,*) 471 DO IRREP = 1, N_SUPSYM_IRREP 472 WRITE(6,*) ' Supersymmetries for IRREP = ', IRREP 473 IB = IBSUPSYM_FOR_IRREP(IRREP) 474 NDEG = NSUPSYM_FOR_IRREP(IRREP) 475 CALL IWRTMA(ISUPSYM_FOR_IRREP(IB),1,NDEG,1,NDEG) 476 END DO 477 END IF 478* 479 RETURN 480 END 481 SUBROUTINE SYM_AND_IRREP_FOR_LMLPA(L,ML,IPA,ISYM,IRREP) 482* 483* An orbital has label L and ML, and parity IPA 484* obtain corresponding symmetry and irrep number 485* 486*. Jeppe Olsen, May 2012 487* June 2012: parity added 488* 489 INCLUDE 'implicit.inc' 490 INCLUDE 'mxpdim.inc' 491 INCLUDE 'crun.inc' 492 INCLUDE 'orbinp.inc' 493 INCLUDE 'glbbas.inc' 494 INCLUDE 'wrkspc-static.inc' 495#include "errquit.fh" 496#include "mafdecls.fh" 497#include "global.fh" 498* 499 NTEST = 00 500* 501 CALL SYM_AND_IRREP_FOR_LMLPA_S(L,ML,IPA,ISYM,IRREP, 502 & CSUPSYM,int_mb(KL_FOR_SUPSYM),int_mb(KML_FOR_SUPSYM), 503 & int_mb(KPA_FOR_SUPSYM), 504 & int_mb(KIRREP_FOR_SUPSYM),N_SUPSYM) 505* 506 IF(NTEST.GE.100) THEN 507 WRITE(6,*) ' Output from SYM_AND_IRREP_FOR_LML ' 508 WRITE(6,*) ' Input: L and ML ', L, ML 509 WRITE(6,*) ' Output: ISYM, IRREP ', ISYM, IRREP 510 END IF 511* 512 RETURN 513 END 514C SYM_AND_IRREP_FOR_LMLPA_S(L,ML,IPA,ISYM,IRREP, 515C & CSUPSYM,WORK(KL_FOR_SUPSYM),WORK(KML_FOR_SUPSYM), 516C & WORK(KIRREP_FOR_SUPSYM)) 517 SUBROUTINE SYM_AND_IRREP_FOR_LMLPA_S(L,ML,IPA,ISYM,IRREP, 518 & CSUPSYM,L_FOR_SUPSYM,ML_FOR_SUPSYM, 519 & IPA_FOR_SUPSYM, 520 & IRREP_FOR_SUPSYM,N_SUPSYM) 521* 522* Obtain supersymmetry and irrep for orbital with given L and ML, 523* and parity 524* 525* Jeppe Olsen, May 23, 2012 526* 527 INCLUDE 'implicit.inc' 528 CHARACTER*6 CSUPSYM 529*. General input 530 INTEGER L_FOR_SUPSYM(N_SUPSYM),ML_FOR_SUPSYM(N_SUPSYM) 531 INTEGER IPA_FOR_SUPSYM(N_SUPSYM) 532 INTEGER IRREP_FOR_SUPSYM(N_SUPSYM) 533* 534 NTEST = 00 535 ISYM = 0 536 IF(CSUPSYM.EQ.'ATOMIC') THEN 537 DO ISUPSYM = 1, N_SUPSYM 538 IF(L_FOR_SUPSYM(ISUPSYM).EQ.L.AND. 539 & ML_FOR_SUPSYM(ISUPSYM).EQ.ML ) THEN 540 ISYM = ISUPSYM 541 END IF 542 END DO 543 ELSE IF(CSUPSYM.EQ.'LINEAR') THEN 544 DO ISUPSYM = 1, N_SUPSYM 545 IF(ML_FOR_SUPSYM(ISUPSYM).EQ.ML.AND. 546 & IPA_FOR_SUPSYM(ISUPSYM).EQ.IPA) THEN 547 ISYM = ISUPSYM 548 END IF 549 END DO 550 ELSE 551 WRITE(6,*) ' Unknown type of supersymmetry ', CSUPSYM 552 STOP ' Unknown type of supersymmetry ' 553 END IF 554* 555 IF(ISYM.EQ.0) THEN 556 WRITE(6,*) 557 & ' Supersymmetry was not found for L, ML, IPA = ', 558 & L, ML, IPA 559 STOP ' Supersymmetry was not found for L, ML, IPA ' 560 END IF 561*. And irrep 562 IRREP = IRREP_FOR_SUPSYM(ISYM) 563* 564 IF(NTEST.GE.100) THEN 565 WRITE(6,'(A,5I3)') ' L, ML, IPA (in) ISYM, IRREP(out) = ', 566 & L, ML, IPA, ISYM, IRREP 567 END IF 568* 569 RETURN 570 END 571C GET_SUPSYM_FOR_BASIS(WORK(KISUPSYM_FOR_BAS), 572C & WORK(KNBAS_FOR_SUP_STA_SYM),WORK(KIBBAS_FOR_SUP_STA_SYM), 573C & WORK(KIBAS_FOR_SUP_STA_SYM),WORK(KIRREP_FOR_BAS), 574C & WORK(KNBAS_FOR_IRREP),WORK(KIBBAS_FOR_IRREP), 575C & WORK(KIBAS_FOR_IRREP)) 576 SUBROUTINE GET_SUPSYM_FOR_BASIS( 577 & ISUPSYM_FOR_BAS, NBAS_FOR_SUP_STA_SYM, 578 & IBBAS_FOR_SUP_STA_SYM, IBAS_FOR_SUP_STA_SYM, 579 & ISHELL_FOR_BAS,NBAS_FOR_SHELL,IBBAS_FOR_SHELL, 580 & IBAS_FOR_SHELL, NSUPSYM_FOR_IRREP,IBSUPSYM_FOR_IRREP, 581 & ISUPSYM_FOR_IRREP) 582* 583* Obtain supersymmetry info for basis set from labels 584* of basis functions 585* 586*. Jeppe Olsen, May 23, 2012 587* Last modification; Jeppe Olsen; March 6, 2013; Info on irrep <=> basis added 588* 589 INCLUDE 'implicit.inc' 590 INCLUDE 'mxpdim.inc' 591 INCLUDE 'orbinp.inc' 592 INCLUDE 'lucinp.inc' 593 INCLUDE 'glbbas.inc' 594 INCLUDE 'wrkspc-static.inc' 595#include "errquit.fh" 596#include "mafdecls.fh" 597#include "global.fh" 598 599*. General input 600 CHARACTER*4 AO_CENT, AO_TYPE, CHAR4 601 COMMON/AOLABELS/AO_CENT(MXPORB),AO_TYPE(MXPORB) 602 INTEGER NSUPSYM_FOR_IRREP(*), IBSUPSYM_FOR_IRREP(*) 603 INTEGER ISUPSYM_FOR_IRREP(*) 604*. Output 605 INTEGER ISUPSYM_FOR_BAS(*) 606 INTEGER NBAS_FOR_SUP_STA_SYM(N_SUPSYM,NSMOB), 607 & IBBAS_FOR_SUP_STA_SYM(N_SUPSYM,NSMOB), 608 & IBAS_FOR_SUP_STA_SYM(*) 609 INTEGER ISHELL_FOR_BAS(NTOOB),NBAS_FOR_SHELL(NTOOB) 610 INTEGER IBBAS_FOR_SHELL(NTOOB),IBAS_FOR_SHELL(NTOOB) 611*. Output is also the arrays NBAS_SUPSYM, IBBAS_SUPSYM, IBAS_SUPSYM in /ORBINP/ 612* 613 NTEST = 100 614* 615 CALL GET_SUPSYM_FOR_BASIS_S(ISUPSYM_FOR_BAS, 616 & int_mb(KLVAL_FOR_ORB),int_mb(KMLVAL_FOR_ORB), 617 & int_mb(KPA_FOR_ORB),NTOOB) 618* 619* 620 IZERO = 0 621 CALL ISETVC(NBAS_FOR_SUP_STA_SYM,IZERO,NSMOB*N_SUPSYM) 622 CALL ISETVC(NBAS_SUPSYM,IZERO,N_SUPSYM) 623 IBAS= 0 624 DO ISTASYM = 1, NSMOB 625 NBASS = NTOOBS(ISTASYM) 626 DO IIBAS = 1, NBASS 627 IBAS = IBAS + 1 628 ISUPSYM = ISUPSYM_FOR_BAS(IBAS) 629 NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) = 630 & NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) + 1 631 NBAS_SUPSYM(ISUPSYM) = NBAS_SUPSYM(ISUPSYM) + 1 632 END DO 633 END DO 634* 635 IF(NTEST.GE.10) THEN 636 WRITE(6,*) 637 & ' Number of basis functions per super- and standard-sym' 638 WRITE(6,*) 639 & ' =====================================================' 640 CALL IWRTMA3(NBAS_FOR_SUP_STA_SYM,N_SUPSYM,NSMOB,N_SUPSYM,NSMOB) 641 WRITE(6,*) 642 WRITE(6,*) 643 & ' Number of basis functions per super-sym' 644 WRITE(6,*) 645 & ' =======================================' 646 CALL IWRTMA3(NBAS_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 647 END IF 648* 649*. And then the basisfunctions of a given super and standard symmetry 650* 651 IB = 1 652 DO ISTASYM = 1, NSMOB 653 DO ISUPSYM = 1, N_SUPSYM 654 IBBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) = IB 655 NBAS = NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) 656 IB = IB + NBAS 657 END DO 658 END DO 659* 660 CALL ISETVC(NBAS_FOR_SUP_STA_SYM,IZERO,NSMOB*N_SUPSYM) 661 IBAS = 0 662 DO ISTASYM = 1, NSMOB 663 NBAS = NTOOBS(ISTASYM) 664 DO IIBAS = 1, NBAS 665 IBAS = IBAS + 1 666 ISUPSYM = ISUPSYM_FOR_BAS(IBAS) 667C? WRITE(6,*) ' TEST: ISTASYM, ISUPSYM = ', 668C? & ISTASYM, ISUPSYM 669 NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) = 670 & NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) + 1 671 N = NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) 672 IB = IBBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) 673 IBAS_FOR_SUP_STA_SYM(IB-1+N) = IBAS 674 END DO 675 END DO 676* 677* and the orbitals of a given supersymmetry 678* 679 IB = 1 680 DO ISUPSYM = 1, N_SUPSYM 681 IBBAS_SUPSYM(ISUPSYM) = IB 682 NBAS = NBAS_SUPSYM(ISUPSYM) 683 IB = IB + NBAS 684 END DO 685* 686 CALL ISETVC(NBAS_SUPSYM,IZERO,N_SUPSYM) 687 DO IBAS = 1, NTOOB 688 ISUPSYM = ISUPSYM_FOR_BAS(IBAS) 689 NBAS_SUPSYM(ISUPSYM) = 690 & NBAS_SUPSYM(ISUPSYM) + 1 691 N = NBAS_SUPSYM(ISUPSYM) 692 IB = IBBAS_SUPSYM(ISUPSYM) 693 IBAS_SUPSYM(IB-1+N) = IBAS 694 END DO 695* 696 IF(NTEST.GE.100) THEN 697 WRITE(6,*) ' Basis functions of given super and standard sym ' 698 WRITE(6,*) ' ================================================' 699 DO ISTASYM = 1, NSMOB 700 DO ISUPSYM = 1, N_SUPSYM 701 N = NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) 702 IB = IBBAS_FOR_SUP_STA_SYM(ISUPSYM,ISTASYM) 703 IF(N.NE.0) THEN 704 WRITE(6,'(A,2(1X,I3))') 705 & ' Standard and supersymmetry ', ISTASYM, ISUPSYM 706 CALL IWRTMA3(IBAS_FOR_SUP_STA_SYM(IB),1,N,1,N) 707 END IF 708 END DO 709 END DO 710* 711 WRITE(6,*) ' Basis functions of given supersym ' 712 WRITE(6,*) ' ==================================' 713 DO ISUPSYM = 1, N_SUPSYM 714 N = NBAS_SUPSYM(ISUPSYM) 715 IB = IBBAS_SUPSYM(ISUPSYM) 716 IF(N.NE.0) THEN 717 WRITE(6,'(A,1(1X,I3))') 718 & ' Supersymmetry ', ISUPSYM 719 CALL IWRTMA3(IBAS_SUPSYM(IB),1,N,1,N) 720 END IF 721 END DO 722 END IF 723* 724*. Relation between basis functions /MO's in symmetry order and 725*. Shells in symmetry order 726* 727* 728*. a) The number of the shell corresponding to each 729* basis functions 730*. A bit of inefficient coding (but I do not want to create a scratch array today...) 731* 732 IB_SHELL = 1 733COLD DO IRREP = 1, NIRREP 734 DO IRREP = 1, N_SUPSYM_IRREP 735*. Supersymmmetries of this irrep 736 IB = IBSUPSYM_FOR_IRREP(IRREP) 737 N = NSUPSYM_FOR_IRREP(IRREP) 738 IF(NTEST.GE.10) WRITE(6,*) ' IRREP, IB, N = ', IRREP, IB, N 739 DO IISUPSYM = IB, IB-1+N 740 ISUPSYM = ISUPSYM_FOR_IRREP(IISUPSYM) 741 IBO = IBBAS_SUPSYM(ISUPSYM) 742 NO = NBAS_SUPSYM(ISUPSYM) 743 IBS = IB_SHELL 744 IF(NTEST.GE.10) THEN 745 WRITE(6,*) ' ISUPSYM, IBO, NO, IBS = ', 746 & ISUPSYM, IBO, NO, IBS 747 END IF 748 DO IIORB = IBO, IBO-1+NO 749 IORB = IBAS_SUPSYM(IIORB) 750 ISHELL_FOR_BAS(IORB) = IBS-1+IIORB-IBO+1 751 IF(NTEST.GE.10) THEN 752 WRITE(6,*) ' IIORB, IORB, RHS = ', 753 & IIORB, IORB, IBS-1+IIORB-IBO+1 754 END IF 755 END DO 756 END DO 757 IB_SHELL = IB_SHELL + NO 758 END DO ! loop over IRREPS 759 NTOSH = IB_SHELL - 1 760 WRITE(6,*) ' Total number of shells = ', NTOSH 761* 762 IF(NTEST.GE.0) THEN 763 WRITE(6,*) ' Shell number for orbitals in sym.order' 764 WRITE(6,*) ' ======================================' 765 CALL IWRTMA3(ISHELL_FOR_BAS,1,NTOOB,1,NTOOB) 766 END IF 767* 768 RETURN 769 END 770 SUBROUTINE GET_SUPSYM_FOR_BASIS_S(ISUPSYM_FOR_BAS, 771 & LVAL_FOR_ORB,MLVAL_FOR_ORB,IPA_FOR_ORB, 772 & NTOOB) 773* 774 INCLUDE 'implicit.inc' 775*. Input 776 INTEGER LVAL_FOR_ORB(NTOOB),MLVAL_FOR_ORB(NTOOB) 777 INTEGER IPA_FOR_ORB(NTOOB) 778*. Output 779 INTEGER ISUPSYM_FOR_BAS(NTOOB) 780* 781 NTEST = 00 782 IF(NTEST.GE.10) THEN 783 WRITE(6,*) 784 WRITE(6,*) ' Info for GET_SUPSYM_FOR_BASIS_S' 785 WRITE(6,*) ' ===============================' 786 WRITE(6,*) 787 END IF 788 789* 790 DO IORB = 1, NTOOB 791 L = LVAL_FOR_ORB(IORB) 792 ML = MLVAL_FOR_ORB(IORB) 793 IPA = IPA_FOR_ORB(IORB) 794 CALL SYM_AND_IRREP_FOR_LMLPA(L,ML,IPA,ISYM,IRREP) 795C SYM_AND_IRREP_FOR_LMLPA(L,ML,IPA,ISYM,IRREP) 796 ISUPSYM_FOR_BAS(IORB) = ISYM 797 END DO 798* 799 IF(NTEST.GE.100) THEN 800 WRITE(6,*) ' Basis functions and their supersymmetry ' 801 WRITE(6,*) ' ========================================' 802 WRITE(6,*) 803 WRITE(6,*) ' Orbital Supersymmetry ' 804 WRITE(6,*) ' ======================= ' 805 DO IORB = 1, NTOOB 806 WRITE(6,'(2X,I3,5X,I2)') IORB, ISUPSYM_FOR_BAS(IORB) 807 END DO 808 END IF 809* 810 RETURN 811 END 812 SUBROUTINE GEN_GENSMOB(NGENSMOB,NBAS_GENSMOB,IBBAS_GENSMOB, 813 & ISTA_TO_GENSM_REO, 814 & NBAS_FOR_SUP_STA_SYM,IBBAS_FOR_SUP_STA_SYM, 815 & IBAS_FOR_SUP_STA_SYM,N_SUPSYM,NSMOB,ICHECK) 816* 817* Set the GENSMOB arrays for supersymmetry case 818* 819*. Jeppe Olsen, May 23, 2012 820* 821 INCLUDE 'implicit.inc' 822*. Input 823 INTEGER NBAS_FOR_SUP_STA_SYM(N_SUPSYM,NSMOB), 824 & IBBAS_FOR_SUP_STA_SYM(N_SUPSYM,NSMOB), 825 & IBAS_FOR_SUP_STA_SYM(*) 826*. Output 827 INTEGER NBAS_GENSMOB(N_SUPSYM*NSMOB), 828 & IBBAS_GENSMOB(N_SUPSYM*NSMOB), 829 & ISTA_TO_GENSM_REO(*) 830* 831 NTEST = 100 832 IF(NTEST.GE.100) THEN 833 WRITE(6,*) ' GEN_GENSMOB reporting' 834 WRITE(6,*) ' =====================' 835 WRITE(6,*) 836 WRITE(6,*) ' ICHECK = ', ICHECK 837 WRITE(6,*) ' N_SUPSYM, NSMOB = ', N_SUPSYM, NSMOB 838 END IF 839* 840 NGENSMOB = N_SUPSYM*NSMOB 841 IGENSMOB = 0 842 IB_OUT = 1 843* 844 DO ISMOB = 1, NSMOB 845 DO ISUPSYM = 1, N_SUPSYM 846 IF(NTEST.GE.1000) WRITE(6,'(A,2(2X,I2))') 847 & ' ISMOB, ISUPSUM = ', ISMOB, ISUPSUM 848 IGENSMOB = IGENSMOB + 1 849 N = NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISMOB) 850 IB_IN = IBBAS_FOR_SUP_STA_SYM(ISUPSYM,ISMOB) 851 NBAS_GENSMOB(IGENSMOB) = N 852 IBBAS_GENSMOB(IGENSMOB) = IB_OUT 853 WRITE(6,*) ' TEST: IB, IB_OUT, N = ', 854 & IB, IB_OUT, N 855 CALL ICOPVE( 856 & IBAS_FOR_SUP_STA_SYM(IB),ISTA_TO_GENSM_REO(IB_OUT),N) 857 IB_OUT = IB_OUT + N 858 END DO 859 END DO 860 NORB = IB_OUT - 1 861* 862 IF(NTEST.GE.100) THEN 863 WRITE(6,*) 864 WRITE(6,*) ' Information on GENSMOB arrays ' 865 WRITE(6,*) ' ==============================' 866 WRITE(6,*) 867 WRITE(6,*) ' Number of general orbital symmetries = ', NGENSMOB 868 WRITE(6,*) ' Number of orbitals per general orbital symmetries' 869 CALL IWRTMA3(NBAS_GENSMOB,1,NGENSMOB,1,NGENSMOB) 870 WRITE(6,*) ' Offsets of orbitals for general orbital symmetry' 871 CALL IWRTMA3(IBBAS_GENSMOB,1,NGENSMOB,1,NGENSMOB) 872 WRITE(6,*) ' Reorder array, standard to general symmetry ' 873 CALL IWRTMA3(ISTA_TO_GENSM_REO,1,NORB,1,NSMOB) 874 END IF 875* 876 RETURN 877 END 878 SUBROUTINE REFORM_CMO_STA_GEN(CMO_STA,CMO_GEN, 879 & IDO_REORDER,IREO,IWAY) 880* 881* Reform CMO matrix between standard and general symmetry forms 882* 883* IWAY = 1: Standard => general symmetry 884* IWAY = 2: general symmetry => standard 885* 886* IF IDO_REORDER = 0, then the orbitals are ordered as 887* the basis-functions 888* = 1, the reorder array IREO from the above order 889* is used 890*. Jeppe Olsen, May 24, 2012 891* 892 INCLUDE 'implicit.inc' 893 INCLUDE 'mxpdim.inc' 894 INCLUDE 'orbinp.inc' 895 INCLUDE 'lucinp.inc' 896* 897*. Input and output 898 DIMENSION CMO_STA(*), CMO_GEN(*) 899 INTEGER IREO(*) 900* 901 NTEST = 00 902 IF(NTEST.GE.100) THEN 903 WRITE(6,*) ' Info from REFORM_CMO_STA_GEN' 904 WRITE(6,*) ' ============================' 905 END IF 906* 907 IOFF_C_GEN = 1 908 DO IGENSM = 1, NGENSMOB 909 NGEN = NBAS_GENSMOB(IGENSM) 910 IBGEN = IBBAS_GENSMOB(IGENSM) 911 ISYM = ISTASM_FOR_GENSM(IGENSM) 912 NSTA = NTOOBS(ISYM) 913 IF(NTEST.GE.1000) WRITE(6,'(A,4I4)') 914 & 'IGENSM, ISYM, NGEN, NSTA = ', IGENSM, ISYM, NGEN, NSTA 915*. Start of symmetry block in C of standard expansion 916 IBSTA_C = 1 917 DO JSYM = 1, ISYM-1 918 IBSTA_C = IBSTA_C + NTOOBS(JSYM)**2 919 END DO 920*. Start of orbitals of symmetry ISYM 921 IBSTA = ITOOBS(ISYM) 922* 923 DO IGEN = 1, NGEN 924 ISTA = ISTA_TO_GENSM_REO(IBGEN-1+IGEN) 925 IF(IDO_REORDER.EQ.1) THEN 926 ISTA = IREO(ISTA) 927 END IF 928 IOFF_C_STA = IBSTA_C-1 + (ISTA-IBSTA)*NSTA + 1 929C REFORM_SINGLE_CMO_STA_GEN( CMO_STA,CMO_GEN,ISYM,IGENSM,IWAY) 930 IF(NTEST.GE.1000) WRITE(6,'(A,2I6)') 931 & ' IOFF_C_STA, IOFF_C_GEN = ', IOFF_C_STA, IOFF_C_GEN 932 CALL REFORM_SINGLE_CMO_STA_GEN( 933 & CMO_STA(IOFF_C_STA),CMO_GEN(IOFF_C_GEN),ISYM,IGENSM,IWAY) 934 IOFF_C_GEN = IOFF_C_GEN + NGEN 935 END DO !loop over orbitals of given general symmetry 936 END DO! Loop over general symmetries 937* 938 IF(NTEST.GE.100) THEN 939 WRITE(6,*) ' CMO coefficients in general form ' 940 CALL APRBLM2(CMO_GEN,NBAS_GENSMOB,NBAS_GENSMOB,NGENSMOB,0) 941 WRITE(6,*) ' CMO coefficients in standard form ' 942 CALL APRBLM2(CMO_STA,NTOOBS,NTOOBS,NSMOB,0) 943 END IF 944* 945 RETURN 946 END 947 SUBROUTINE REFORM_SINGLE_CMO_STA_GEN 948 & (CMO_STA,CMO_GEN,ISYM,IGENSM,IWAY) 949* 950* Reform between standard and general symmetry form of expansions of 951* a single orbital of symmetry ISYM and general symmetry IGENSM 952* 953* IWAY = 1: Standard => general symmetry 954* IWAY = 2: general symmetry => standard 955* 956*. Jeppe Olsen, May 23 - FCN has just become Danish Champions in soccer- 957* FCK is down to silver... 958* 959 INCLUDE 'implicit.inc' 960*. Input and output: Coefs of a single MO 961 DIMENSION CMO_STA(*), CMO_GEN(*) 962* 963 INCLUDE 'mxpdim.inc' 964 INCLUDE 'orbinp.inc' 965* 966 NTEST = 000 967 IF(NTEST.GE.100) WRITE(6,*) ' Info from REFORM_CMO... ' 968* 969 NGEN = NBAS_GENSMOB(IGENSM) 970 IBGEN = IBBAS_GENSMOB(IGENSM) 971 NSTA = NTOOBS(ISYM) 972 IF(NTEST.GE.1000) WRITE(6,'(A,4I4)') 973 &' IGENSM, ISYM, NGEN, NSTA = ', IGENSM, ISYM, NGEN, NSTA 974*. Start for given symmetry 975 IB_SYM = ITOOBS(ISYM) 976* 977 IF(IWAY.EQ.1) THEN 978*. Standard => supersymmetry 979 DO IGEN = 1, NGEN 980 CMO_GEN(IGEN) = 981 & CMO_STA(ISTA_TO_GENSM_REO(IBGEN-1+IGEN)-IB_SYM+1) 982 END DO 983 ELSE 984*. Supersymmetry to standard 985 ZERO = 0.0D0 986 CALL SETVEC(CMO_STA,ZERO,NSTA) 987 DO IGEN = 1, NGEN 988 CMO_STA(ISTA_TO_GENSM_REO(IBGEN-1+IGEN)-IB_SYM+1) = 989 & CMO_GEN(IGEN) 990 END DO 991 END IF 992* 993 IF(NTEST.GE.100) THEN 994 IF(IWAY.EQ.1) THEN 995 WRITE(6,*) ' Standard => supersymmetry ' 996 ELSE 997 WRITE(6,*) ' Supersymmetry => standard ' 998 END IF 999 WRITE(6,*) ' Single MO in standard form ' 1000 CALL WRTMAT(CMO_STA,1,NSTA,1,NSTA) 1001 WRITE(6,*) ' Single MO in general symmetry form ' 1002 CALL WRTMAT(CMO_GEN,1,NGEN,1,NGEN) 1003 END IF 1004* 1005 RETURN 1006 END 1007 SUBROUTINE REFORM_MAT_STA_SUP(ASTA,ASUP,IPACK,IWAY) 1008* 1009* Reform a matrix between standard and supersymmetry order 1010* 1011* IWAY = 1: Standard => supersymmetry 1012* IWAY = 2: Supersymmetry => standard 1013* 1014* IPACK = 0: Full blocked matrix 1015* = 1: Lower half blocked matrix 1016* 1017*. Jeppe Olsen, May 23 1018* 1019 INCLUDE 'implicit.inc' 1020 INCLUDE 'mxpdim.inc' 1021 INCLUDE 'orbinp.inc' 1022 INCLUDE 'lucinp.inc' 1023*. Input and output 1024 DIMENSION ASTA(*),ASUP(*) 1025* 1026 NTEST = 00 1027 IF(NTEST.GE.100) THEN 1028 WRITE(6,*) 1029 WRITE(6,*) ' Output from REFORM_MAT_STA_SUP ' 1030 WRITE(6,*) ' ===============================' 1031 WRITE(6,*) 1032 WRITE(6,*) ' IWAY, IPACK = ', IWAY, IPACK 1033 END IF 1034* 1035 IB_STA = 1 1036 IB_SUP = 1 1037 IJ_SUP = 0 1038 IJB_SYM = 1 1039 DO ISYM = 1, NSMOB 1040 IF(NTEST.GE.1000) THEN 1041 WRITE(6,*) 1042 WRITE(6,*) ' Info for ISYM = ', ISYM 1043 WRITE(6,*) ' ======================= ' 1044 WRITE(6,*) 1045 END IF 1046 NSTA = NTOOBS(ISYM) 1047 IBSTA = ITOOBS(ISYM) 1048 IF(IWAY.EQ.2) THEN 1049*. Zero symmetry block 1050 IF(IPACK.EQ.0) THEN 1051 NSMBLK = NSTA*NSTA 1052 ELSE 1053 NSMBLK = NSTA*(NSTA+1)/2 1054 END IF 1055 ZERO = 0.0D0 1056C? WRITE(6,*) ' IJB_SYM, NSMBLK = ', IJB_SYM, NSMBLK 1057 CALL SETVEC(ASTA(IJB_SYM),ZERO,NSMBLK) 1058 END IF 1059 NSUPSYM = NSUP_FOR_STA_SYM(ISYM) 1060 IBSUPSYM = IBSUP_FOR_STA_SYM(ISYM) 1061 DO IISUPSYM = IBSUPSYM, IBSUPSYM - 1 + NSUPSYM 1062 ISUPSYM = ISUP_FOR_STA_SYM(IISUPSYM) 1063 IF(NTEST.GE.1000) 1064 & WRITE(6,*) ' IISUPSYM, ISUPSYM= ', IISUPSYM, ISUPSYM 1065 NSUP = NBAS_GENSMOB(ISUPSYM) 1066 IBSUP = IBBAS_GENSMOB(ISUPSYM) 1067*. Start of general symmetry block ISUPSYM 1068 IB_SUP_MAT = 1 1069 DO JSUPSYM = 1, ISUPSYM-1 1070 NS = NBAS_GENSMOB(JSUPSYM) 1071 IF(IPACK.EQ.0) THEN 1072 IB_SUP_MAT = IB_SUP_MAT + NS**2 1073 ELSE 1074 IB_SUP_MAT = IB_SUP_MAT + NS*(NS+1)/2 1075 END IF 1076 END DO 1077 IJ_SUP = IB_SUP_MAT - 1 1078* 1079 DO ISUP = 1, NSUP 1080 IF(IPACK.EQ.0) THEN 1081 JSUP_MAX = NSUP 1082 ELSE 1083 JSUP_MAX = ISUP 1084 END IF 1085 DO JSUP = 1, JSUP_MAX 1086 IJ_SUP = IJ_SUP + 1 1087*. Indeces of orbitals and matrix in standard form 1088 I_STA = ISTA_TO_GENSM_REO(IBSUP-1+ISUP)-IBSTA+1 1089 J_STA = ISTA_TO_GENSM_REO(IBSUP-1+JSUP)-IBSTA+1 1090 IF(NTEST.GE.1000) 1091 & WRITE(6,*) ' I_STA, J_STA =', I_STA, J_STA 1092 IF(IPACK.EQ.0) THEN 1093 IJ_STA = (J_STA-1)*NSTA+I_STA 1094 ELSE 1095 IJ_STA = I_STA*(I_STA-1)/2 + J_STA 1096 END IF 1097 IF(NTEST.GE.1000) WRITE(6,*) 1098 & ' ISUP, JSUP, IJ_STA = ', ISUP, JSUP, IJ_STA 1099 IF(NTEST.GE.1000) WRITE(6,*) 1100 & 'IJ_SUP, IJB_SYM-1+IJ_STA ', IJ_SUP, IJB_SYM-1+IJ_STA 1101 IF(IWAY.EQ.1) THEN 1102 ASUP(IJ_SUP) = ASTA(IJB_SYM-1+IJ_STA) 1103 ELSE 1104 ASTA(IJB_SYM-1+IJ_STA) = ASUP(IJ_SUP) 1105 END IF 1106 END DO! Loop over JSUP 1107 END DO ! Loop over ISUP 1108 END DO! Loop over super symmetry 1109*. Update pointer to start of standard 1110 IF(IPACK.EQ.0) THEN 1111 IJB_SYM = IJB_SYM + NSTA**2 1112 ELSE 1113 IJB_SYM = IJB_SYM + NSTA*(NSTA+1)/2 1114 END IF 1115C? WRITE(6,*) ' NSTA and updated IJB_SYM = ', NSTA, IJB_SYM 1116 END DO ! End of loop over Symmetries 1117* 1118 IF(NTEST.GE.100) THEN 1119 IF(IWAY.EQ.1) THEN 1120 WRITE(6,*) ' Standard => supersymmetry ' 1121 ELSE 1122 WRITE(6,*) ' Supersymmetry => Standard' 1123 END IF 1124 WRITE(6,*) ' Matrix in standard symmetry form ' 1125 CALL APRBLM2(ASTA,NTOOBS,NTOOBS,NSMOB,IPACK) 1126 WRITE(6,*) ' Matrix in Super symmetry form ' 1127 CALL APRBLM2(ASUP,NBAS_GENSMOB, NBAS_GENSMOB,NGENSMOB,IPACK) 1128 END IF 1129* 1130 RETURN 1131 END 1132 SUBROUTINE N_SUPSYM_IRREP_TO_SUPSYM(N_PER_IRREP,N_PER_SUPSYM) 1133* 1134* Obtain number of orbitals per super symmetry from 1135* number of shells per supersymmetry irrep 1136* 1137*. Jeppe Olsen, May 23, 2012 1138* 1139 INCLUDE 'implicit.inc' 1140 INCLUDE 'mxpdim.inc' 1141 INCLUDE 'orbinp.inc' 1142 INCLUDE 'glbbas.inc' 1143 INCLUDE 'wrkspc-static.inc' 1144#include "errquit.fh" 1145#include "mafdecls.fh" 1146#include "global.fh" 1147*. Input 1148 INTEGER N_PER_IRREP(*) 1149* Output 1150 INTEGER N_PER_SUPSYM(*) 1151* 1152 NTEST = 000 1153 IF(NTEST.GE.100) THEN 1154 WRITE(6,*) ' Info from N_SUPSYM_IRREP_TO_SUPSYM ' 1155 WRITE(6,*) ' ================================== ' 1156 WRITE(6,*) 1157 END IF 1158* 1159 CALL N_SUPSYM_IRREP_TO_SUPSYM_S( 1160 & N_PER_IRREP,N_PER_SUPSYM, 1161 & N_SUPSYM_IRREP, N_SUPSYM, 1162 & int_mb(KNSUPSYM_FOR_IRREP), 1163 & int_mb(KIBSUPSYM_FOR_IRREP), 1164 & int_mb(KISUPSYM_FOR_IRREP) ) 1165* 1166 IF(NTEST.GE.100) THEN 1167 WRITE(6,*) ' Number of shells per supsym irrep ' 1168 CALL IWRTMA3(N_PER_IRREP,1,N_SUPSYM_IRREP,1,N_SUPSYM_IRREP) 1169 WRITE(6,*) ' Number of orbitals per supersymmetry ' 1170 CALL IWRTMA3(N_PER_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 1171 END IF 1172* 1173 RETURN 1174 END 1175 SUBROUTINE N_SUPSYM_IRREP_TO_SUPSYM_S( 1176 & N_PER_IRREP,N_PER_SUPSYM, 1177 & N_SUPSYM_IRREP, N_SUPSYM, 1178 & NSUPSYM_FOR_IRREP, 1179 & IBSUPSYM_FOR_IRREP, 1180 & ISUPSYM_FOR_IRREP) 1181* 1182* Obtain number of orbitals per super symmetry from 1183* number of shells per supersymmetry irrep - slave routine 1184* 1185*. Jeppe Olsen, May 23, 2012 1186* 1187 INCLUDE 'implicit.inc' 1188*. Input 1189 INTEGER N_PER_IRREP(N_SUPSYM_IRREP) 1190*. Output 1191 INTEGER N_PER_SUPSYM(N_SUPSYM) 1192*. General info 1193 INTEGER NSUPSYM_FOR_IRREP(*), 1194 & IBSUPSYM_FOR_IRREP(*), 1195 & ISUPSYM_FOR_IRREP(*) 1196* 1197 NTEST = 00 1198 IF(NTEST.GE.100) WRITE(6,*) 'From N_SUPSYM_IRREP_TO_SUPSYM_S' 1199 IZERO = 0 1200 CALL ISETVC(N_PER_SUPSYM,IZERO,N_SUPSYM) 1201* 1202 DO IRREP = 1, N_SUPSYM_IRREP 1203 NSHL = N_PER_IRREP(IRREP) 1204 NCOMP = NSUPSYM_FOR_IRREP(IRREP) 1205 IBCOMP = IBSUPSYM_FOR_IRREP(IRREP) 1206 IF(NTEST.GE.1000) 1207 & WRITE(6,*) ' IRREP, NCOMP, IBCOMP = ', IRREP, NCOMP, IBCOMP 1208 DO IICOMP = IBCOMP, IBCOMP-1+NCOMP 1209 ICOMP = ISUPSYM_FOR_IRREP(IICOMP) 1210 IF(NTEST.GE.1000) 1211 & WRITE(6,*) ' IICOMP, ICOMP = ', IICOMP, ICOMP 1212 N_PER_SUPSYM(ICOMP) = NSHL 1213 END DO 1214 END DO 1215* 1216 RETURN 1217 END 1218 SUBROUTINE SUP_TO_STASYM(NBAS_FOR_SUP_STA_SYM) 1219* 1220* Super-symmetry components for each standard symmetry 1221* 1222* Output is in /ORBIBP/: NSUP_FOR_STA_SYM,IBSUP_FOR_STA_SYM,ISUP_FOR_STA_SYM, 1223* ISTASM_FOR_SUPSYM 1224*. Jeppe Olsen, May 24, 2012 1225* 1226* 1227 INCLUDE 'implicit.inc' 1228 INCLUDE 'mxpdim.inc' 1229 INCLUDE 'orbinp.inc' 1230 INCLUDE 'lucinp.inc' 1231*. Specific input 1232 INTEGER NBAS_FOR_SUP_STA_SYM(N_SUPSYM,NSMOB) 1233* 1234 DO ISYM = 1, NSMOB 1235 NCOMP = 0 1236 DO ISUPSYM = 1, N_SUPSYM 1237 IF(NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISYM).NE.0) NCOMP = NCOMP + 1 1238 END DO 1239 NSUP_FOR_STA_SYM(ISYM) = NCOMP 1240 END DO 1241* 1242 IB = 1 1243 DO ISYM = 1, NSMOB 1244 IBSUP_FOR_STA_SYM(ISYM) = IB 1245 IB = IB + NSUP_FOR_STA_SYM(ISYM) 1246 END DO 1247* 1248 DO ISYM = 1, NSMOB 1249 IB = IBSUP_FOR_STA_SYM(ISYM) 1250 NCOMP = 0 1251 DO ISUPSYM = 1, N_SUPSYM 1252 IF(NBAS_FOR_SUP_STA_SYM(ISUPSYM,ISYM).NE.0) THEN 1253 NCOMP = NCOMP + 1 1254 ISUP_FOR_STA_SYM(IB-1+NCOMP) = ISUPSYM 1255 ISTASM_FOR_SUPSYM(ISUPSYM) = ISYM 1256 END IF 1257 END DO 1258 END DO 1259* 1260 NTEST = 100 1261 IF(NTEST.GE.1000) THEN 1262 WRITE(6,*) ' Supersymmetries for standard symmetry ' 1263 WRITE(6,*) ' ======================================' 1264 WRITE(6,*) 1265 DO ISYM = 1, NSMOB 1266 WRITE(6,*) 1267 WRITE(6,*) ' Supersymmetries for standard symmetry', ISYM 1268 IB = IBSUP_FOR_STA_SYM(ISYM) 1269 N = NSUP_FOR_STA_SYM(ISYM) 1270 CALL IWRTMA3(ISUP_FOR_STA_SYM(IB),1,N,1,N) 1271 END DO 1272 END IF 1273* 1274 IF(NTEST.GE.100) THEN 1275 WRITE(6,*) ' Standard symmetry of supersymmetries ' 1276 WRITE(6,*) ' ==================================== ' 1277 WRITE(6,*) 1278 WRITE(6,*) ' Super Standard ' 1279 WRITE(6,*) ' ==================' 1280 DO ISUP = 1, N_SUPSYM 1281 WRITE(6,'(2X,I2,8X,I2)') ISUP, ISTASM_FOR_SUPSYM(ISUP) 1282 END DO 1283 END IF 1284* 1285 RETURN 1286 END 1287 SUBROUTINE REFORM_MAT_STA_GEN(ASTA,AGEN,IPACK,IWAY) 1288* 1289* Reform a matrix between standard and general symmetry order 1290* as defined by the *GENSM* arrays - may be super or standard 1291* symmetry 1292* 1293* IWAY = 1: Standard => supersymmetry 1294* IWAY = 2: Supersymmetry => standard 1295* 1296* IPACK = 0: Full blocked matrix 1297* = 1: Lower half blocked matrix 1298* 1299*. Jeppe Olsen, May 23, rewritten May 24 1300*. Last modification: Oct. 1, 2012: Jeppe Olsen: Removed bug for nonsymmetric matrices 1301* 1302 INCLUDE 'implicit.inc' 1303 INCLUDE 'mxpdim.inc' 1304 INCLUDE 'orbinp.inc' 1305 INCLUDE 'lucinp.inc' 1306*. Input and output 1307 DIMENSION ASTA(*),AGEN(*) 1308* 1309 NTEST = 00 1310 IF(NTEST.GE.100) THEN 1311 WRITE(6,*) 1312 WRITE(6,*) ' Output from REFORM_MAT_STA_GEN ' 1313 WRITE(6,*) ' ===============================' 1314 WRITE(6,*) 1315 WRITE(6,*) ' ISTA_TO_GENSM_REO ' 1316 CALL IWRTMA3(ISTA_TO_GENSM_REO,1,NTOOB,1,NTOOB) 1317 END IF 1318* 1319 IF(IWAY.EQ.2) THEN 1320*. Zero matrix 1321 LEN_STA = LEN_BLMAT(NSMOB,NTOOBS,NTOOBS,IPACK) 1322 ZERO = 0.0D0 1323 CALL SETVEC(ASTA,ZERO,LEN_STA) 1324 END IF 1325*. Loop over blocks of general symmetry 1326 IB_GEN_MAT = 1 1327 IJ_GEN = 0 1328 DO IGENSM = 1, NGENSMOB 1329*. Standard symmetry 1330 ISTASM = ISTASM_FOR_GENSM(IGENSM) 1331 IB_STA_ORB = ITOOBS(ISTASM) 1332 NSTA = NTOOBS(ISTASM) 1333 IF(NTEST.GE.1000) THEN 1334 WRITE(6,*) ' IGENSM, ISTASM = ', IGENSM, ISTASM 1335 END IF 1336 IB_STA_MAT = 1 1337 DO JSTASM = 1, ISTASM-1 1338 N = NTOOBS(JSTASM) 1339 IF(IPACK.EQ.0) THEN 1340 IB_STA_MAT = IB_STA_MAT + N**2 1341 ELSE 1342 IB_STA_MAT = IB_STA_MAT + N*(N+1)/2 1343 END IF 1344 END DO 1345 IF(NTEST.GE.1000) WRITE(6,*) ' IB_STA_MAT = ', IB_STA_MAT 1346*. Loop over pairs of orbitals in general symmetry block 1347 NGEN = NBAS_GENSMOB(IGENSM) 1348 IBGEN = IBBAS_GENSMOB(IGENSM) 1349 IJ_GEN0 = IJ_GEN 1350 DO IBASGN = 1, NGEN 1351 IF(IPACK.EQ.0) THEN 1352 JBASGN_MAX = NGEN 1353 ELSE 1354 JBASGN_MAX = IBASGN 1355 END IF 1356 DO JBASGN = 1, JBASGN_MAX 1357*. The corresponding orbitals in standard order 1358 IBAS = ISTA_TO_GENSM_REO(IBGEN-1+IBASGN)-IB_STA_ORB + 1 1359 JBAS = ISTA_TO_GENSM_REO(IBGEN-1+JBASGN)-IB_STA_ORB + 1 1360 IF(IPACK.EQ.0) THEN 1361 IJ_STA = IB_STA_MAT - 1 + (JBAS-1)*NSTA + IBAS 1362 ELSE 1363 IJ_STA = IB_STA_MAT - 1 + IBAS*(IBAS-1)/2 + JBAS 1364 END IF 1365 IF(NTEST.GE.1000) 1366 & WRITE(6,'(A,4I3)') ' IBASGN, IBAS, JBASGN, JBAS = ', 1367 & IBASGN, IBAS, JBASGN, JBAS 1368* 1369 IF(IPACK.EQ.0) THEN 1370 IJ_GEN = IJ_GEN0 + (JBASGN-1)*NGEN + IBASGN 1371 ELSE 1372 IJ_GEN = IJ_GEN + 1 1373 END IF 1374 IF(NTEST.GE.1000) 1375 & WRITE(6,'(A,2I5)') ' IJ_GEN, IJ_STA = ', IJ_GEN, IJ_STA 1376 IF(IWAY.EQ.1) THEN 1377 AGEN(IJ_GEN) = ASTA(IJ_STA) 1378 ELSE 1379 ASTA(IJ_STA) = AGEN(IJ_GEN) 1380 END IF 1381 END DO! loop over IBASGN 1382 END DO! loop over JBASGN 1383 END DO! loop over IGENSM 1384* 1385 IF(NTEST.GE.100) THEN 1386 IF(IWAY.EQ.1) THEN 1387 WRITE(6,*) ' Standard => supersymmetry ' 1388 ELSE 1389 WRITE(6,*) ' Supersymmetry => Standard' 1390 END IF 1391 WRITE(6,*) ' Matrix in standard symmetry form ' 1392 CALL APRBLM2(ASTA,NTOOBS,NTOOBS,NSMOB,IPACK) 1393 WRITE(6,*) ' Matrix in general symmetry form ' 1394 CALL APRBLM2(AGEN,NBAS_GENSMOB, NBAS_GENSMOB,NGENSMOB,IPACK) 1395 END IF 1396* 1397 RETURN 1398 END 1399 SUBROUTINE N_SUPSYM_TO_STASYM(N_PER_SUPSYM,N_PER_STASYM) 1400* A list of integers for supersymmetries are given 1401* Reform to list of integers over standard symmetries 1402* 1403*. Jeppe Olsen, May 24, 2012 1404* 1405 INCLUDE 'implicit.inc' 1406 INCLUDE 'mxpdim.inc' 1407 INCLUDE 'orbinp.inc' 1408 INCLUDE 'lucinp.inc' 1409*. Input 1410 INTEGER N_PER_SUPSYM(*) 1411*. Output 1412 INTEGER N_PER_STASYM(*) 1413* 1414 NTEST = 00 1415* 1416 DO ISTASYM = 1, NSMOB 1417 NSUP = NSUP_FOR_STA_SYM(ISTASYM) 1418 IBSUP = IBSUP_FOR_STA_SYM(ISTASYM) 1419 N = 0 1420 DO IISUPSYM = IBSUP, IBSUP + NSUP - 1 1421 ISUPSYM = ISUP_FOR_STA_SYM(IISUPSYM) 1422 IF(NTEST.GE.1000) 1423 & WRITE(6,'(A,3I4)') ' ISTASYM, IISUPSYM, ISUPSYM = ', 1424 & ISTASYM, IISUPSYM, ISUPSYM 1425 N = N + N_PER_SUPSYM(ISUPSYM) 1426 END DO 1427 N_PER_STASYM(ISTASYM) = N 1428 END DO 1429* 1430 IF(NTEST.GE.100) THEN 1431 WRITE(6,*) ' Output from N_SUPSYM_TO_STASYM:' 1432 WRITE(6,*) ' ===============================' 1433 WRITE(6,*) ' Integer list over supersymmetries (input) ' 1434 CALL IWRTMA3(N_PER_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 1435 WRITE(6,*) ' Integer list over standard symmetries (output) ' 1436 CALL IWRTMA3(N_PER_STASYM,1,NSMOB,1,NSMOB) 1437 END IF 1438* 1439 RETURN 1440 END 1441 SUBROUTINE ORDER_SUPSYM_ORBITALS(NSPC,ISPC,MO_SUPSYM,IREO, 1442 & ISUPSYM_FOR_BAS) 1443* 1444* A set of orbital spaces ISPC defined in terms of 1445* supersymmetries are given. Order the orbitals accordingly. 1446* 1447*. Well, the deal is that the basis functions within a 1448*. given standard symmetry is defined from input, whereas the 1449*. ordering of the molecular orbitals in a given standard 1450*. symmetry is not specified. LUCIA does as a standard 1451*. use a canonical order where the MO's are ordered in the 1452*. same way as the basis functions. 1453*. However, we typically define say the occupied orbitals 1454* or the CAS orbitals of a given symmetry as the lowest in 1455* a given stadard symmetry. To accomplish this, it is 1456* useful to introduce a reordering of the orbitals 1457* in a given symmetry. THis is what this routine is about... 1458* 1459* Jeppe Olsen, May 24, 2012 1460* 1461 INCLUDE 'implicit.inc' 1462 INCLUDE 'mxpdim.inc' 1463 INCLUDE 'orbinp.inc' 1464 INCLUDE 'lucinp.inc' 1465*. Input 1466 INTEGER ISPC(MXP_NSUPSYM,NSPC), ISUPSYM_FOR_BAS(*) 1467*. Output: IREO: New order index from standard order index 1468 INTEGER MO_SUPSYM(NTOOB),IREO(*) 1469*. Local scratch: dim of number of irreps 1470 INTEGER ISCR1(1000),ISCR2(1000),ISCR3(1000) 1471* 1472 NTEST = 100 1473 IF(NTEST.GE.100) THEN 1474 WRITE(6,*) ' Output from ORDER_SUPSYM_ORBITALS ' 1475 WRITE(6,*) ' The requested division of orbitals' 1476 CALL IWRTMA3(ISPC,N_SUPSYM,NSPC,MXP_NSUPSYM,NSPC) 1477 END IF 1478* 1479* Obtain first MO_SUPSYM 1480* 1481 IORB = 0 1482 DO ISYM = 1, NSMOB 1483 NSUP = NSUP_FOR_STA_SYM(ISYM) 1484 IBSUP = IBSUP_FOR_STA_SYM(ISYM) 1485*. Loop over the spac division for this symmetry 1486 DO JSPC = 1, NSPC 1487 DO IISUPSYM = IBSUP, IBSUP + NSUP -1 1488 ISUPSYM = ISUP_FOR_STA_SYM(IISUPSYM) 1489 NSUPSPC = ISPC(ISUPSYM,JSPC) 1490 NORB_SUPSPC = ISPC(ISUPSYM,JSPC) 1491 DO IIORB = 1, NORB_SUPSPC 1492 IORB = IORB + 1 1493 MO_SUPSYM(IORB) = ISUPSYM 1494 IF(NTEST.GE.10000) 1495 & WRITE(6,*) ' IORB, ISUPSYM ', IORB, ISUPSYM 1496 END DO! loop over orbital IORB 1497 END DO! Loop over supersymmetries 1498 END DO ! Loop over orbital spaces 1499 END DO ! Loop over standard symmetries 1500* 1501 IF(NTEST.GE.100) THEN 1502 WRITE(6,*) ' Supersymmetry of occ-ordered orbitals' 1503 CALL IWRTMA3(MO_SUPSYM,1,NTOOB,1,NTOOB) 1504 END IF 1505* 1506*. Counting index in ISCR2 for orbital with given supersymmetry 1507* 1508 IZERO = 0 1509 CALL ISETVC(ISCR1,IZERO,N_SUPSYM) 1510 DO IORB = 1, NTOOB 1511 ISUPSYM = MO_SUPSYM(IORB) 1512 ICOUNT = ISCR1(ISUPSYM)+1 1513 ISCR2(IORB) = ICOUNT 1514 ISCR1(ISUPSYM) = ISCR1(ISUPSYM) + 1 1515 END DO 1516 IF(NTEST.GE.1000) THEN 1517 WRITE(6,*) ' Counting index of reordered orbitals' 1518 CALL IWRTMA3(ISCR2,1,NTOOB,1,NTOOB) 1519 END IF 1520* 1521*. Loop over orbitals in standard order and find the corresponding new 1522*. order. Not elegant(quadratic scaling!!), but I'm to tired for elegance 1523* 1524 IZERO = 0 1525 CALL ISETVC(ISCR1,IZERO, N_SUPSYM) 1526 MONE = -1 1527 CALL ISETVC(IREO,MONE,NTOOB) 1528 DO ISYM = 1, NSMOB 1529 NORB = NTOOBS(ISYM) 1530 IBORB = ITOOBS(ISYM) 1531 DO IORB_STA = IBORB, IBORB + NORB -1 1532 ISUPSYM = ISUPSYM_FOR_BAS(IORB_STA) 1533C? WRITE(6,*) ' ISUPSYM_FOR_BAS(1) = ', ISUPSYM_FOR_BAS(1) 1534 IF(NTEST.GE.1000) 1535 & WRITE(6,*) ' IORB_STA, ISUPSYM = ', IORB_STA, ISUPSYM 1536 ISCR1(ISUPSYM) = ISCR1(ISUPSYM) + 1 1537 ICOUNT = ISCR1(ISUPSYM) 1538 IF(NTEST.GE.1000) 1539 & WRITE(6,'(A,3(1X,I2))') ' IORB_STA, ISUPSUM, ICOUNT = ', 1540 & IORB_STA, ISUPSYM, ICOUNT 1541*. Find orbital with this supersymmetry and count number 1542 DO IORB = 1, NTOOB 1543 IF(ISUPSYM.EQ.MO_SUPSYM(IORB).AND. 1544 & ICOUNT.EQ.ISCR2(IORB)) THEN 1545*. Match 1546 IREO(IORB_STA) = IORB 1547 END IF 1548 END DO! Loop over orbitals in new order 1549 END DO 1550 END DO 1551* 1552 IF(NTEST.GE.100) THEN 1553 WRITE(6,*) ' Reordering of orbitals (standard => occ/gas)' 1554 CALL IWRTMA3(IREO,1,NTOOB,1,NTOOB) 1555 END IF 1556* 1557 RETURN 1558 END 1559 SUBROUTINE REO_CMOAO(CIN,COUT,IREO,ICOPY,IWAY) 1560* 1561* Reorder coefficients of MOAO matrix 1562* 1563* IWAY = 1: COUT(IREO(I)) = CIN(I) 1564* IWAY = 2: COUT(I) = CIN(IREO(I) 1565* 1566*. If ICOPY = 1, COUT is copied over CIN 1567* 1568* Jeppe Olsen, May 25, 2012 1569* 1570 INCLUDE 'implicit.inc' 1571 INCLUDE 'mxpdim.inc' 1572 INCLUDE 'orbinp.inc' 1573 INCLUDE 'lucinp.inc' 1574*. Input (and perhaps output) 1575 DIMENSION CIN(*), IREO(NTOOB) 1576*. Output 1577 DIMENSION COUT(*) 1578* 1579 NTEST = 00 1580 IF(NTEST.GE.100) THEN 1581 WRITE(6,*) ' Info from REO_CMOAO' 1582 WRITE(6,*) ' ===================' 1583 WRITE(6,*) 1584 WRITE(6,*) ' Input matrix ' 1585C CALL APRBLM2(CIN,NTOOBS,NTOOBS,NSMOB,0) 1586 CALL APRBLM_F7(CIN,NTOOBS,NTOOBS,NSMOB,0) 1587 WRITE(6,*) ' Reorder array ' 1588 CALL IWRTMA3(IREO,1,NTOOB,1,NTOOB) 1589 END IF 1590* 1591 IB_SYM_MAT = 1 1592 DO ISYM = 1, NSMOB 1593 NORB = NTOOBS(ISYM) 1594 IBORB = ITOOBS(ISYM) 1595 DO IORB_IN = IBORB, IBORB + NORB - 1 1596 IORB_REO = IREO(IORB_IN) 1597 IOFF_IN = IB_SYM_MAT + (IORB_IN-IBORB)*NORB 1598 IOFF_REO = IB_SYM_MAT + (IORB_REO-IBORB)*NORB 1599 IF(NTEST.GE.1000) WRITE(6,'(A,4(1X,I3))') 1600 & 'IORB_IN, IORB_REO, IOFF_IN,IOFF_REO = ', 1601 & IORB_IN, IORB_REO, IOFF_IN,IOFF_REO 1602 IF(IWAY.EQ.1) THEN 1603 CALL COPVEC(CIN(IOFF_IN),COUT(IOFF_REO),NORB) 1604 ELSE 1605 CALL COPVEC(CIN(IOFF_REO),COUT(IOFF_IN),NORB) 1606 END IF 1607 END DO ! Loop over IORB_IN 1608 IB_SYM_MAT = IB_SYM_MAT + NORB*NORB 1609 END DO! Loop over ISYM 1610* 1611 LEN_C = IB_SYM_MAT 1612 IF(ICOPY.EQ.1) CALL COPVEC(COUT,CIN,LEN_C) 1613* 1614 IF(NTEST.GE.100) THEN 1615 WRITE(6,*) 1616 WRITE(6,*) ' Reordered MO matrix ' 1617 WRITE(6,*) ' ====================' 1618 WRITE(6,*) 1619 IF(ICOPY.EQ.1) THEN 1620C CALL APRBLM2(CIN,NTOOBS,NTOOBS,NSMOB,0) 1621 CALL APRBLM_F7(CIN,NTOOBS,NTOOBS,NSMOB,0) 1622 ELSE 1623C CALL APRBLM2(COUT,NTOOBS,NTOOBS,NSMOB,0) 1624 CALL APRBLM_F7(COUT,NTOOBS,NTOOBS,NSMOB,0) 1625 END IF 1626 END IF 1627* 1628 RETURN 1629 END 1630 SUBROUTINE SET_HF_DIST_SUPSYM 1631* 1632* Define the matrix HF_DSV_SUPSYM,HF_DSV_STASYM,HF_DSV_GNSYM 1633* giving number of orbitals per supsym, stasym, gnsym 1634* in Doubly, Singly and virtual spaces 1635* 1636* Input and output are all in ORBINP 1637* 1638* Jeppe Olsen, May 26, 2012 1639* 1640 INCLUDE 'implicit.inc' 1641 INCLUDE 'mxpdim.inc' 1642 INCLUDE 'orbinp.inc' 1643 INCLUDE 'lucinp.inc' 1644* 1645 NTEST = 100 1646 IF(NTEST.GE.100) THEN 1647 WRITE(6,*) 1648 WRITE(6,*) ' Info from SET_HF_DIST_SUPSYM ' 1649 WRITE(6,*) ' =============================' 1650 WRITE(6,*) 1651 WRITE(6,*) ' NHFD_SUPSYM, NHFS_SUPSYM, NBAS_SUPSYM( input ) ' 1652 CALL IWRTMA3(NHFD_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 1653 CALL IWRTMA3(NHFS_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 1654 CALL IWRTMA3(NBAS_SUPSYM,1,N_SUPSYM,1,N_SUPSYM) 1655* 1656 WRITE(6,*) ' NHFD_STASYM, NHFS_STASYM, ( input ) ' 1657 CALL IWRTMA3(NHFD_STASYM,1,NSMOB,1,MXPOBS) 1658 CALL IWRTMA3(NHFS_STASYM,1,NSMOB,1,MXPOBS) 1659 END IF 1660* 1661*. Just copy the info on doubly and singly occupied orbitals 1662* 1663 CALL ICOPVE(NHFD_SUPSYM,NHF_DSV_SUPSYM(1,1),N_SUPSYM) 1664 CALL ICOPVE(NHFS_SUPSYM,NHF_DSV_SUPSYM(1,2),N_SUPSYM) 1665* 1666 CALL ICOPVE(NHFD_STASYM,NHF_DSV_STASYM(1,1),N_SUPSYM) 1667 CALL ICOPVE(NHFS_STASYM,NHF_DSV_STASYM(1,2),N_SUPSYM) 1668*. And the virtual spaces as the NBAS - NDOUBLE - NSINGLE 1669 IONE = 1 1670 IMONE = -1 1671C IVCSUM(IA,IB,IC,IFACB,IFACC,NDIM) 1672 CALL IVCSUM(NHF_DSV_SUPSYM(1,3),NBAS_SUPSYM,NHF_DSV_SUPSYM(1,1), 1673 & IONE,IMONE,N_SUPSYM) 1674 CALL IVCSUM(NHF_DSV_SUPSYM(1,3),NHF_DSV_SUPSYM(1,3), 1675 & NHF_DSV_SUPSYM(1,2),IONE,IMONE,N_SUPSYM) 1676* 1677 CALL IVCSUM(NHF_DSV_STASYM(1,3),NTOOBS,NHF_DSV_STASYM(1,1), 1678 & IONE,IMONE,NSMOB) 1679 CALL IVCSUM(NHF_DSV_STASYM(1,3),NHF_DSV_STASYM(1,3), 1680 & NHF_DSV_STASYM(1,2),IONE,IMONE,NSMOB) 1681* 1682 CALL ICOPVE(NHF_DSV_SUPSYM(1,1),NHF_DSV_GNSYM(1,1),N_SUPSYM) 1683 CALL ICOPVE(NHF_DSV_SUPSYM(1,2),NHF_DSV_GNSYM(1,2),N_SUPSYM) 1684 CALL ICOPVE(NHF_DSV_SUPSYM(1,3),NHF_DSV_GNSYM(1,3),N_SUPSYM) 1685* 1686 IF(NTEST.GE.100) THEN 1687 WRITE(6,*) 1688 WRITE(6,*) 1689 & ' Doubly, singly, and unoccupied orbitals per supersymmetry' 1690 WRITE(6,*) 1691 & ' =========================================================' 1692 CALL IWRTMA3(NHF_DSV_SUPSYM,N_SUPSYM,3,MXP_NSUPSYM,3) 1693 WRITE(6,*) 1694 WRITE(6,*) 1695 & ' Doubly, singly, and unoccupied orbitals per standard symmetry' 1696 WRITE(6,*) 1697 & ' =============================================================' 1698 CALL IWRTMA3(NHF_DSV_STASYM,NSMOB,3,MXPOBS,3) 1699 WRITE(6,*) 1700 WRITE(6,*) 1701 & ' Doubly, singly, and unoccupied orbitals per gensymmetry' 1702 WRITE(6,*) 1703 & ' =========================================================' 1704 CALL IWRTMA3(NHF_DSV_GNSYM,N_SUPSYM,3,MXP_NSUPSYM,3) 1705 END IF 1706* 1707 RETURN 1708 END 1709 SUBROUTINE SET_HF_DIST_STASYM 1710* 1711* Define the matrix HF_DSV_STASYM,HF_DSV_GNSYM 1712* giving number of orbitals per supsym, stasym, gnsym 1713* in Doubly, Singly and virtual spaces 1714* 1715* Input and output are all in ORBINP. 1716* The number of doubly and singly occupied orbitals were given in 1717* NHFD_IRREP_SUPSYM, NHFS_IRREP_SUPSYM - not logical, I admit 1718* 1719* Jeppe Olsen, May 26, 2012 1720* 1721 INCLUDE 'implicit.inc' 1722 INCLUDE 'mxpdim.inc' 1723 INCLUDE 'orbinp.inc' 1724 INCLUDE 'lucinp.inc' 1725* 1726 NTEST = 100 1727 IF(NTEST.GE.100) THEN 1728 WRITE(6,*) 1729 WRITE(6,*) ' Info from SET_HF_DIST_STASYM ' 1730 WRITE(6,*) ' =============================' 1731 WRITE(6,*) 1732 END IF 1733* 1734*. Just copy the info on doubly and singly occupied orbitals 1735* 1736 CALL ICOPVE(NHFD_IRREP_SUPSYM, NHFD_SUPSYM,NSMOB) 1737 CALL ICOPVE(NHFS_IRREP_SUPSYM, NHFS_SUPSYM,NSMOB) 1738* 1739 CALL ICOPVE(NHFD_SUPSYM,NHF_DSV_STASYM(1,1),NSMOB) 1740 WRITE(6,*) ' NHF_DSV_STASYM(*,1): ' 1741 CALL IWRTMA3(NHF_DSV_STASYM(1,1),1,NSMOB,1,NSMOB) 1742 CALL ICOPVE(NHFS_SUPSYM,NHF_DSV_STASYM(1,2),NSMOB) 1743*. And the virtual spaces as the NBAS - NDOUBLE - NSINGLE 1744 IONE = 1 1745 IMONE = -1 1746C IVCSUM(IA,IB,IC,IFACB,IFACC,NDIM) 1747* 1748 CALL IVCSUM(NHF_DSV_STASYM(1,3),NTOOBS,NHF_DSV_STASYM(1,1), 1749 & IONE,IMONE,NSMOB) 1750 CALL IVCSUM(NHF_DSV_STASYM(1,3),NHF_DSV_STASYM(1,3), 1751 & NHF_DSV_STASYM(1,2),IONE,IMONE,NSMOB) 1752* 1753 CALL ICOPVE(NHF_DSV_STASYM(1,1),NHF_DSV_GNSYM(1,1),NSMOB) 1754 CALL ICOPVE(NHF_DSV_STASYM(1,2),NHF_DSV_GNSYM(1,2),NSMOB) 1755 CALL ICOPVE(NHF_DSV_STASYM(1,3),NHF_DSV_GNSYM(1,3),NSMOB) 1756* 1757 IF(NTEST.GE.100) THEN 1758 WRITE(6,*) 1759 WRITE(6,*) 1760 & ' Doubly, singly, and unoccupied orbitals per standard symmetry' 1761 WRITE(6,*) 1762 & ' =============================================================' 1763 CALL IWRTMA3(NHF_DSV_STASYM,NSMOB,3,MXPOBS,3) 1764 WRITE(6,*) 1765 WRITE(6,*) 1766 & ' Doubly, singly, and unoccupied orbitals per gensymmetry' 1767 WRITE(6,*) 1768 & ' =========================================================' 1769 CALL IWRTMA3(NHF_DSV_GNSYM,NSMOB,3,MXP_NSUPSYM,3) 1770 END IF 1771* 1772 RETURN 1773 END 1774 SUBROUTINE ORDER_GAS_SUPSYM_ORBITALS 1775* 1776* Obtain the order of the orbitals according to the 1777* symmetry specified by NGAS_IRREP_SUPSYM 1778* 1779* Jeppe Olsen, May 26 1780* 1781 INCLUDE 'implicit.inc' 1782 INCLUDE 'mxpdim.inc' 1783 INCLUDE 'orbinp.inc' 1784 INCLUDE 'lucinp.inc' 1785 INCLUDE 'cgas.inc' 1786 INCLUDE 'glbbas.inc' 1787 INCLUDE 'wrkspc-static.inc' 1788#include "errquit.fh" 1789#include "mafdecls.fh" 1790#include "global.fh" 1791 1792* 1793 NTEST = 100 1794 IF(NTEST.GE.100) THEN 1795 WRITE(6,*) 1796 WRITE(6,*) ' Output from ORDER_GAS_SUPSYM_ORBITALS ' 1797 WRITE(6,*) ' ===================================== ' 1798 WRITE(6,*) 1799 END IF 1800*. Reform from irrep to super-symmetry and standard symmetry 1801 DO IGAS = 0, NGAS + 1 1802 CALL N_SUPSYM_IRREP_TO_SUPSYM( 1803 & NGAS_IRREP_SUPSYM(1,IGAS),NGAS_SUPSYM(1,IGAS)) 1804 CALL N_SUPSYM_TO_STASYM( 1805 & NGAS_SUPSYM(1,IGAS),NGAS_STASYM(1,IGAS)) 1806 END DO 1807* 1808 IF(NTEST.GE.100) THEN 1809 WRITE(6,*) ' GA-spaces over supersymmetries ' 1810 WRITE(6,*) ' ===============================' 1811 WRITE(6,*) 1812 CALL IWRTMA3(NGAS_SUPSYM,N_SUPSYM,NGAS+2,MXP_NSUPSYM,NGAS+2) 1813 WRITE(6,*) ' GA-spaces over standard symmetries ' 1814 WRITE(6,*) ' ===================================' 1815 WRITE(6,*) 1816 CALL IWRTMA3(NGAS_STASYM,NSMOB,NGAS+2,MXPOBS,NGAS+2) 1817 END IF 1818* 1819* Obtain the reordering for this set of orbitals 1820* 1821 CALL ORDER_SUPSYM_ORBITALS(NGAS+2,NGAS_SUPSYM, 1822 & int_mb(KMO_SUPSYM),int_mb(KMO_STA_TO_ACT_REO), 1823 & int_mb(KISUPSYM_FOR_BAS) ) 1824* 1825 IF(NTEST.GE.10) THEN 1826 WRITE(6,*) ' Super-symmetry reordering of orbitals, GAS ' 1827 CALL IWRTMA3(int_mb(KMO_STA_TO_ACT_REO),1,NTOOB,1,NTOOB) 1828 END IF 1829* 1830 RETURN 1831 END 1832 SUBROUTINE SET_PARITY_FOR_STASYM(IPA_FOR_STASYM,NIRREP) 1833* 1834* Define parity (+1/-1) for the various standard symmetries 1835* Currently assuming D2H in the standard MOLECULE order 1836* 1837* Jeppe Olsen, June 2012 1838* 1839 INCLUDE 'implicit.inc' 1840*. Output 1841 INTEGER IPA_FOR_STASYM(*) 1842* 1843 NTEST = 100 1844* 1845 IF(NIRREP.EQ.8) THEN 1846*. the gerade symmetries 1847 IPA_FOR_STASYM(1) = 1 1848 IPA_FOR_STASYM(4) = 1 1849 IPA_FOR_STASYM(6) = 1 1850 IPA_FOR_STASYM(7) = 1 1851*. and the ungerade symmetries 1852 IPA_FOR_STASYM(2) =-1 1853 IPA_FOR_STASYM(3) =-1 1854 IPA_FOR_STASYM(5) =-1 1855 IPA_FOR_STASYM(8) =-1 1856 ELSE 1857 WRITE(6,*) ' SET_PARITY_FOR_STASYM not set for NIRREP = ',NIRREP 1858 STOP ' SET_PARITY_FOR_STASYM not set for NIRREP ' 1859 END IF 1860* 1861 IF(NTEST.GE.100) THEN 1862 WRITE(6,*) ' Parity of standard symmetries ' 1863 CALL IWRTMA3(IPA_FOR_STASYM,1,NIRREP,1,NIRREP) 1864 END IF 1865* 1866 RETURN 1867 END 1868 SUBROUTINE SUPSYM_FROM_CMOAO(CMOAO,ISUPSYM_FOR_BAS, 1869 & ISUPSYM_FOR_MOS) 1870* 1871* A CMOAO matrix is given in CMOAO and the supersymmetry of the basis functions is given in 1872* ISUPSYM_FOR_BAS. Obtain the supersymmetry of the MO's 1873* 1874*. Jeppe Olsen, July 3, 2012 1875*. Last modification; Oct. 3, 2012; Jeppe Olsen; Removed bug 1876* 1877 INCLUDE 'implicit.inc' 1878 INCLUDE 'mxpdim.inc' 1879 INCLUDE 'orbinp.inc' 1880 INCLUDE 'lucinp.inc' 1881*. Input 1882 INTEGER ISUPSYM_FOR_BAS(*) 1883 DIMENSION CMOAO(*) 1884*. Output 1885 INTEGER ISUPSYM_FOR_MOS(NTOOB) 1886* 1887 NTEST = 100 1888 IF(NTEST.GE.100) THEN 1889 WRITE(6,*) ' Info from SUPSYM_FROM_CMOAO ' 1890 WRITE(6,*) ' ============================' 1891 END IF 1892 IF(NTEST.GE.1000) THEN 1893 WRITE(6,*) 1894 WRITE(6,*) ' Input CMO ' 1895 WRITE(6,*) 1896 CALL PRINT_CMOAO(CMOAO) 1897 END IF 1898* 1899 IJOFF = 1 1900 IORB = 0 1901 DO ISM = 1, NSMOB 1902 NI = NTOOBS(ISM) 1903 DO I = 1, NI 1904 IORB = IORB + 1 1905 ICOFF = IJOFF-1+(I-1)*NI + 1 1906 ISUPSYM = ISUPSYM_FOR_MO(CMOAO(ICOFF),ISM,ISUPSYM_FOR_BAS) 1907 ISUPSYM_FOR_MOS(IORB) = ISUPSYM 1908 END DO 1909 IJOFF = IJOFF + NI**2 1910 END DO! loop over ISM 1911*. Check that all orbitals were assigned supesymmetry 1912 NZERO = 0 1913 DO IOB = 1, NTOOB 1914 IF(ISUPSYM_FOR_MOS(IOB).EQ.0) NZERO = NZERO + 1 1915 END DO 1916* 1917 IF(NTEST.GE.100.OR.NZERO.NE.0) THEN 1918 WRITE(6,*) ' Supersymmetries determined from CMOAO matrix ' 1919 CALL IWRTMA3(ISUPSYM_FOR_MOS,1,NTOOB,1,NTOOB) 1920 END IF 1921* 1922 IF(NZERO.NE.0) THEN 1923 WRITE(6,*) ' Not all orbitals have well-defined super-symmetry' 1924 WRITE(6,*) ' Number of orbitals with unassigned super-sym ', 1925 & NZERO 1926 STOP ' Not all orbitals have well-defined super-symmetry' 1927 END IF 1928* 1929 RETURN 1930 END 1931 FUNCTION ISUPSYM_FOR_MO(CMO,ISM, ISUPSYM_FOR_BAS) 1932* 1933* The MO-coefficients of an MO of standard symmetry ISM is given in CMO. 1934*. Obtain the super symmetry of this orbital. If the orbital does not have 1935* a well-defined super-symmetry, a zero is returned. 1936* 1937*. Jeppe Olsen, July 3, 2012 1938* 1939 INCLUDE 'implicit.inc' 1940 INCLUDE 'mxpdim.inc' 1941 INCLUDE 'orbinp.inc' 1942*. Input 1943 INTEGER ISUPSYM_FOR_BAS(*) 1944 DIMENSION CMO(*) 1945* 1946 NTEST = 00 1947 IF(NTEST.GE.100) THEN 1948 WRITE(6,*) ' Info from ISUPSYM_FOR_MO ' 1949 WRITE(6,*) ' ==========================' 1950 END IF 1951 IF(NTEST.GE.1000) THEN 1952 WRITE(6,*) ' Standard symmetry of CMO ', ISM 1953 WRITE(6,*) ' Expansion CMO: ' 1954 NI = NTOOBS(ISM) 1955 CALL WRTMAT_F7(CMO,1,NI,1,NI) 1956 END IF 1957* 1958*. Offset to basis functions of standard symmetry ISM 1959 I_OFF = 1 1960 DO IISM = 1, ISM-1 1961 I_OFF = I_OFF + NTOOBS(IISM) 1962 END DO 1963 NI = NTOOBS(ISM) 1964*. Find supersymmetry of basis function with largest expansion coefficient 1965 ISUPSYM_MAX = 0 1966 IMAX_BAS = 0 1967 C_MAX = 0.0D0 1968 DO IBAS = 1, NI 1969 IF(ABS(CMO(IBAS)).GT.ABS(C_MAX))THEN 1970 C_MAX = CMO(IBAS) 1971 IMAX_BAS = IBAS 1972 END IF 1973 END DO 1974* 1975 IF(C_MAX.EQ.0.0D0) THEN 1976*. vanishing expansion, supersymmetry is set to zero 1977 WRITE(6,*) ' Warning: vanishing CMO in ISUPSYM_FOR_MO' 1978 ISUPSYM = 0 1979 ELSE 1980 ISUPSYM = ISUPSYM_FOR_BAS(I_OFF-1+IMAX_BAS) 1981 IF(NTEST.GE.1000) 1982 & WRITE(6,*) ' Super-sym of bf with max coef ', ISUPSYM 1983*. Make sure that all non-vanishing coefficients have the same supersymmetry 1984*.(could be relaxed by using a non-vanishing threshold) 1985*. Note, I am pt using a rather high threshold, so some contamination 1986*. is allowed. Done in connection with some MCVB calculations 1987 THRES = 1.0D-6 1988 NMISS = 0 1989 CSUM = 0.0D0 1990 DO IBAS = 1, NI 1991 IF(ABS(CMO(IBAS)).GT.THRES) THEN 1992 JSUPSYM = ISUPSYM_FOR_BAS(I_OFF-1+IBAS) 1993 IF(NTEST.GE.1000) 1994 & WRITE(6,*) ' IBAS, JSUPSYM = ', IBAS,JSUPSYM 1995 IF(JSUPSYM.NE.ISUPSYM) THEN 1996 NMISS = NMISS + 1 1997 CSUM = CSUM + ABS(CMO(IBAS)) 1998 ILAST = IBAS 1999 END IF 2000 END IF 2001 END DO 2002 IF(NMISS.NE.0) THEN 2003*. Not well defined supersymmetry, set to zero 2004 WRITE(6,*) ' Problem, ISM, NI, ISUPSYM = ', ISM, NI, ISUPSYM 2005 WRITE(6,*) 2006 & ' Number of coefficients with deviating supersymmetry ', NMISS 2007 WRITE(6,*) 2008 & ' Sum of coefficients with deviating supersymmetry ', CSUM 2009 WRITE(6,*) 2010 & ' Last included basis function with wrong supersymmetry ', 2011 & ILAST 2012 WRITE(6,*) ' Expansion CMO: ' 2013 CALL WRTMAT_F7(CMO,1,NI,1,NI) 2014 ISUPSYM = 0 2015 END IF 2016 END IF ! CMO was nonvanishing 2017* 2018 ISUPSYM_FOR_MO = ISUPSYM 2019* 2020 IF(NTEST.GE.100) THEN 2021 WRITE(6,*) ' Supersymmetry of MO: ', ISUPSYM 2022 END IF 2023* 2024 RETURN 2025 END 2026 SUBROUTINE EXTR_CP_GASBLKS_FROM_GENSYM_MAT 2027 & (AS,ASG,IEORC,IGAS_F,IGAS_L,IPAK) 2028* 2029* Reform between two diagonal blockings of orbital matrices 2030* S: Blocks over general symmetry, all blocks between IGAS_F and 2031* IGAS_L 2032* SG: Blocks over general symmetry and GAS 2033* 2034*. Jeppe Olsen, July 2012 2035*. Last modification, Sept. 24 2012, Jeppe Olsen, Debugged.. 2036* 2037 INCLUDE 'implicit.inc' 2038 INCLUDE 'mxpdim.inc' 2039 INCLUDE 'orbinp.inc' 2040*. Input and output 2041 DIMENSION AS(*),ASG(*) 2042* 2043 NTEST = 00 2044 IF(NTEST.GE.100) THEN 2045 WRITE(6,*) ' Info from EXTR_CP_GASBLKS_FROM_GENSYM_MAT' 2046 WRITE(6,*) ' =========================================' 2047 END IF 2048* 2049 IJOFF_S = 1 2050 IJOFF_SG = 1 2051 IOFF_S = 1 2052 IOFF_SG = 1 2053* 2054 DO IGENSM = 1, NGENSMOB 2055 IF(NTEST.GE.1000) WRITE(6,*) ' IGENSM = ', IGENSM 2056*. Number of orbitals in AFULL for this supersymmetry 2057 NORB_S = 0 2058 DO IGAS = IGAS_F, IGAS_L 2059 NORB_S = NORB_S + NGAS_GNSYM(IGENSM,IGAS) 2060 END DO 2061 IOFF_S = 1 2062 IOFF_SG = 1 2063 DO IGAS = IGAS_F, IGAS_L 2064 IF(NTEST.GE.1000) WRITE(6,*) ' IGAS = ', IGAS 2065 NORB_SG = NGAS_GNSYM(IGENSM,IGAS) 2066*. Loop over pairs of orbitals in this symmetry-gas block 2067 DO IORB_SG = 1,NORB_SG 2068 IORB_S = IOFF_SG - 1 + IORB_SG 2069 IF(NTEST.GE.1000) 2070 & WRITE(6,*) ' IORB_S, IORB_SG = ', IORB_S, IORB_SG 2071 IF(IPAK.EQ.1) THEN 2072 JORB_MX = IORB_SG 2073 ELSE 2074 JORB_MX = NORB_SG 2075 END IF 2076 DO JORB_SG = 1, JORB_MX 2077 JORB_S = IOFF_SG - 1 + JORB_SG 2078 IF(NTEST.GE.1000) 2079 & WRITE(6,*) ' JORB_S, JORB_SG = ', JORB_S, JORB_SG 2080 IF(IPAK.EQ.0) THEN 2081 IJ_SG = (JORB_SG-1)*NORB_SG + IORB_SG 2082 IJ_S = (JORB_S-1)*NORB_S + IORB_S 2083 ELSE 2084 IJ_SG = IORB_SG*(IORB_SG-1)/2 + JORB_SG 2085 IJ_S = IORB_S*(IORB_S-1)/2 + JORB_S 2086 END IF 2087 IF(NTEST.GE.1000) 2088 & WRITE(6,*) ' IJ_SG, IJ_S = ', IJ_SG, IJ_S 2089 IJ_SG_ABS = IJOFF_SG - 1 + IJ_SG 2090 IJ_S_ABS = IJOFF_S - 1 + IJ_S 2091 IF(NTEST.GE.1000) 2092 & WRITE(6,*) ' IJ_SG_ABS, IJ_S_ABS = ', IJ_SG_ABS, IJ_S_ABS 2093* 2094 IF(IEORC.EQ.1) THEN 2095 AS(IJ_S_ABS) = ASG(IJ_SG_ABS) 2096 ELSE 2097 ASG(IJ_SG_ABS) = AS(IJ_S_ABS) 2098 END IF 2099 END DO ! Loop over JORB_SG 2100 END DO ! Loop over IORB_SG 2101* 2102 IOFF_SG = IOFF_SG + NORB_SG 2103 IF(IPAK.EQ.0) THEN 2104 IJOFF_SG = IJOFF_SG + NORB_SG**2 2105 ELSE 2106 IJOFF_SG = IJOFF_SG + NORB_SG*(NORB_SG+1)/2 2107 END IF 2108 IF(NTEST.GE.1000) 2109 & WRITE(6,*) ' Updated IJOFF_SG = ', IJOFF_SG 2110 END DO ! loop over IGAS 2111* 2112 IF(IPAK.EQ.0) THEN 2113 IJOFF_S = IJOFF_S + NORB_S**2 2114 ELSE 2115 IJOFF_S = IJOFF_S + NORB_S*(NORB_S+1)/2 2116 END IF 2117 IF(NTEST.GE.1000) 2118 & WRITE(6,*) ' Updated IJOFF_S = ', IJOFF_S 2119 END DO !loop over IGENSM 2120* 2121 IF(NTEST.GE.100) THEN 2122 WRITE(6,*) ' Matrix over general symmetries ' 2123 CALL WRT_SG_MAT(AS,1,IGAS_F,IGAS_L,IPAK,1) 2124 WRITE(6,*) ' Matrix over general symmetries and Gasblocks' 2125 CALL WRT_SG_MAT(ASG,2,IGAS_F,IGAS_L,IPAK,1) 2126 END IF 2127* 2128 RETURN 2129 END 2130 SUBROUTINE WRT_SG_MAT(A,IS_OR_SG,IGAS_F,IGAS_L,IPAK,IEXT) 2131* 2132* Matrix A consists of diagonal general symmetry-blocks (IS_OR_SG = 1) 2133* or general symmetry-gas blocks(IS_OR_SG=2), both for gas-paces between 2134* IGAS_F and IGAS_L. Print this 2135* IEXT = 0 => Compact (F7) 2136* IEXT = 1 => Standard output 2137* 2138*. Jeppe Olsen, July 2012 2139* 2140 INCLUDE 'implicit.inc' 2141 INCLUDE 'mxpdim.inc' 2142 INCLUDE 'orbinp.inc' 2143*. Input 2144 DIMENSION A(*) 2145* 2146 IOFF = 1 2147 DO IGENSM = 1, NGENSMOB 2148 NI_S = 0 2149 DO IGAS = IGAS_F, IGAS_L 2150 NI_S = NI_S + NGAS_GNSYM(IGENSM, IGAS) 2151 END DO 2152 WRITE(6,'(A, I2)') ' Symmetry block ', IGENSM 2153 IF(IS_OR_SG.EQ.1) THEN 2154 IF(IPAK.EQ.0) THEN 2155 IF(IEXT.EQ.0) THEN 2156 CALL WRTMAT_F7(A(IOFF),NI_S,NI_S,NI_S,NI_S) 2157 ELSE 2158 CALL WRTMAT(A(IOFF),NI_S,NI_S,NI_S,NI_S) 2159 END IF 2160 ELSE 2161 IF(IEXT.EQ.0) THEN 2162 CALL PRSYM_F7(A(IOFF),NI_S) 2163 ELSE 2164 CALL PRSYM(A(IOFF),NI_S) 2165 END IF 2166 END IF 2167 ELSE IF(IS_OR_SG .EQ. 2) THEN 2168 DO IGAS = IGAS_F, IGAS_L 2169 WRITE(6,*) ' Diagonal block with GAS = ', IGAS 2170 NI_SG = NGAS_GNSYM(IGENSM,IGAS) 2171C? WRITE(6,*) ' IGENSM, IGAS, NI_SG = ', 2172C? & IGENSM, IGAS, NI_SG 2173 IF(IPAK.EQ.0) THEN 2174 IF(IEXT.EQ.0) THEN 2175 CALL WRTMAT_F7(A(IOFF),NI_SG,NI_SG,NI_SG,NI_SG) 2176 ELSE 2177 CALL WRTMAT(A(IOFF),NI_SG,NI_SG,NI_SG,NI_SG) 2178 END IF 2179 IOFF = IOFF + NI_SG**2 2180 ELSE 2181 IF(IEXT.EQ.0) THEN 2182 CALL PRSYM_F7(A(IOFF),NI_SG) 2183 ELSE 2184 CALL PRSYM(A(IOFF),NI_SG) 2185 END IF 2186 IOFF = IOFF + NI_SG*(NI_SG+1)/2 2187 END IF 2188 END DO 2189 END IF !IS_OR_SG switch 2190* 2191 IF(IS_OR_SG.EQ.1) THEN 2192 IF(IPAK.EQ.0) THEN 2193 IOFF = IOFF + NI_S**2 2194 ELSE 2195 IOFF = IOFF + NI_S*(NI_S+1)/2 2196 END IF 2197 END IF 2198* 2199 END DO ! Loop over general symmetry blocks 2200* 2201 RETURN 2202 END 2203 SUBROUTINE REFORM_RHO1_TO_GNSM( 2204 & RHO1_ST,RHO1_GNSM_ST,IWAY,IREO_GNSYM_TO_TS) 2205* 2206* Reform between standard and general symmetry-order 2207* of total symmetric density matrix 2208* 2209*. Jeppe Olsen, July 2012 2210* 2211* Last modified, July 8, 2012 (Jeppe) 2212* 2213* IWAY = 1 => Standard form to general symmetry blocked 2214* IWAY = 2 => General symmetry blockes to standard from 2215* 2216* 2217 INCLUDE 'implicit.inc' 2218 INCLUDE 'mxpdim.inc' 2219 INCLUDE 'orbinp.inc' 2220 INCLUDE 'cgas.inc' 2221*. Input 2222 INTEGER IREO_GNSYM_TO_TS(*) 2223*. Input and output 2224 DIMENSION RHO1_ST(*), RHO1_GNSM_ST(*) 2225* 2226 NTEST = 00 2227 IF(NTEST.GE.100) THEN 2228 WRITE(6,*) ' Info from REFORM_RHO1_TO_GNSM ' 2229 WRITE(6,*) ' IWAY = ', IWAY 2230 WRITE(6,*) ' IREO_GNSYM_TO_TS: ' 2231 CALL IWRTMA(IREO_GNSYM_TO_TS,1,NACOB,1,NACOB) 2232 END IF 2233* 2234 IJOFF_S = 1 2235 IOFF_S = 1 2236 DO IGENSM = 1, NGENSMOB 2237 NACT_S = 0 2238 DO IGAS = 1, NGAS 2239 NACT_S = NACT_S + NGAS_GNSYM(IGENSM,IGAS) 2240 END DO 2241 DO IORB_S = 1, NACT_S 2242 DO JORB_S = 1, NACT_S 2243* 2244 IORB = IREO_GNSYM_TO_TS(IOFF_S-1+IORB_S) 2245 JORB = IREO_GNSYM_TO_TS(IOFF_S-1+JORB_S) 2246C? WRITE(6,*) ' IGENSM, IORB_S, JORB_S, IORB, JORB = ', 2247C? & IGENSM, IORB_S, JORB_S, IORB, JORB 2248* 2249 IF(IWAY.EQ.1) THEN 2250 RHO1_GNSM_ST(IJOFF_S-1+(JORB_S-1)*NACT_S+IORB_S) = 2251 & RHO1_ST((JORB-1)*NACOB + IORB) 2252 ELSE 2253 RHO1_ST((JORB-1)*NACOB + IORB) = 2254 & RHO1_GNSM_ST(IJOFF_S-1+(JORB_S-1)*NACT_S+IORB_S) 2255 END IF ! IWAY switch 2256 END DO 2257 END DO! End of loops over orbitals IORB_S, JORB_S 2258 IOFF_S = IOFF_S + NACT_S 2259 IJOFF_S = IJOFF_S + NACT_S**2 2260 END DO ! Loop over symmetries 2261* 2262 IF(NTEST.GE.100) THEN 2263 WRITE(6,*) ' 1-body density as NACOB x NACOB matrix' 2264 CALL WRTMAT(RHO1_ST,NACOB,NACOB,NACOB,NACOB) 2265 WRITE(6,*) 2266 WRITE(6,*) ' 1-body density as blocks over general symmetry' 2267C WRT_SG_MAT(A,IS_OR_SG,IGAS_F,IGAS_L,IPAK,IEXT) 2268 CALL WRT_SG_MAT(RHO1_GNSM_ST,1,1,NGAS,0,1) 2269 END IF 2270* 2271 RETURN 2272 END 2273 SUBROUTINE REO_ACT_ORB_TO_GNSM(IMO_GNSYM, IREO_GNSYM_TO_TS) 2274* 2275* Reorder array for active orbitals from General symmetry to 2276* standard order for active orbitals. 2277* 2278*. Jeppe Olsen, July 2012 2279* 2280 INCLUDE 'implicit.inc' 2281 INCLUDE 'mxpdim.inc' 2282 INCLUDE 'orbinp.inc' 2283 INCLUDE 'cgas.inc' 2284*. Input 2285 INTEGER IMO_GNSYM(NTOOB) 2286*. Output 2287 INTEGER IREO_GNSYM_TO_TS(NACOB) 2288* 2289 NTEST = 00 2290* 2291 IF(NTEST.GE.100) THEN 2292 WRITE(6,*) ' Output from REO_ACT_ORB_TO_GNSM' 2293 WRITE(6,*) ' Input: IMO_GNSYM ' 2294 CALL IWRTMA3(IMO_GNSYM,1,NTOOB,1,NTOOB) 2295 END IF 2296* 2297 IOFF_GNSYM = 1 2298 IACOB = 0 2299 DO IGNSYM = 1, NGENSMOB 2300C? WRITE(6,*) ' IGNSYM = ', IGNSYM 2301 NINOB_S = NGAS_GNSYM(IGNSYM,0) 2302 NACOB_S = 0 2303 DO IGAS = 1, NGAS 2304 NACOB_S = NACOB_S + NGAS_GNSYM(IGNSYM,IGAS) 2305 END DO 2306C? WRITE(6,*) ' NINOB_S, NACOB_S = ', NINOB_S, NACOB_S 2307 IOB_S = 0 2308*. Find orbitals NINOB_S+1, NINOB_S + NACOB_S of sym IGNSYM 2309 DO IORB = 1, NTOOB 2310 IF(NTEST.GE.1000) 2311 & WRITE(6,*) ' TESTJ: IORB, GNSYM = ', IORB, IMO_GNSYM(IORB) 2312 IF(IMO_GNSYM(IORB).EQ.IGNSYM) THEN 2313 IOB_S = IOB_S + 1 2314 IF(NINOB_S.LT.IOB_S.AND.IOB_S.LE.NINOB_S+NACOB_S) THEN 2315* Orbital IORB is active orbital IACOB in general sym order 2316 IACOB = IACOB + 1 2317*. Orbital IORB address in standard type order 2318 IACOB_STA = IREOST(IORB)-NINOB 2319C? WRITE(6,*) ' IACOB, IACOB_STA ' 2320 IREO_GNSYM_TO_TS(IACOB) = IACOB_STA 2321 END IF 2322 END IF 2323 END DO! loop over IORB 2324 END DO ! Loop over IGNSYM 2325* 2326 IF(NTEST.GE.100) THEN 2327 WRITE(6,*) 2328 & ' Reorder array for active orbs: General => Type order ' 2329 CALL IWRTMA3(IREO_GNSYM_TO_TS,1,NACOB,1,NACOB) 2330 END IF 2331* 2332 RETURN 2333 END 2334 SUBROUTINE LEN_GAS_GS_BLOCKS(LEN_GAS_GS,N_GAS_GS,IGAS_F,IGAS_L) 2335* 2336* Obtain in LEN_GAS_GS the dimension of each symmetry-gas block 2337* for general symmetry 2338* 2339*. Jeppe Olsen, July 2012 2340* 2341*. Last revision: July 8, 2012 2342* 2343 INCLUDE 'implicit.inc' 2344 INCLUDE 'mxpdim.inc' 2345 INCLUDE 'orbinp.inc' 2346 INCLUDE 'cgas.inc' 2347*. Output 2348 INTEGER LEN_GAS_GS((NGAS+2)*NGENSMOB) 2349* 2350 NTEST = 000 2351* 2352 NBLK = 0 2353 DO IGENSM = 1, NGENSMOB 2354 DO IGAS = IGAS_F, IGAS_L 2355 NBLK = NBLK + 1 2356 LEN_GAS_GS(NBLK) = NGAS_GNSYM(IGENSM,IGAS) 2357 END DO 2358 END DO 2359 N_GAS_GS = NBLK 2360* 2361 IF(NTEST.GE.100) THEN 2362 WRITE(6,*) ' Number of General symmetry gas-blocks ', NBLK 2363 WRITE(6,*) ' Length of general-symmetry gas-block ' 2364 CALL IWRTMA3(LEN_GAS_GS,1,NBLK,1,NLBLK) 2365 END IF 2366* 2367 RETURN 2368 END 2369 SUBROUTINE INVERT_REO(IREO,IREO_INV, NDIM) 2370* Obtain inverse mapping of reordering IREO(I) 2371* If K = IREO(I) then IREO_INV(K) = I 2372* 2373* IREO_INV(IREO(I)) = I 2374* 2375*. Jeppe Olsen, Oct. 2, 2012 2376*. Last modification; Oct. 2, 2012; Jeppe Olsen, original version 2377* 2378 IMPLICIT REAL*8(A-H,O-Z) 2379*. Input 2380 INTEGER IREO(NDIM) 2381*. Output 2382 INTEGER IREO_INV(NDIM) 2383* 2384 NTEST = 100 2385 IF(NTEST.GE.100) THEN 2386 WRITE(6,*) ' Output from INVERT_REO ' 2387 END IF 2388* 2389 DO I = 1, NDIM 2390 K = IREO(I) 2391 IREO_INV(K) = I 2392 END DO 2393* 2394 IF(NTEST.GE.100) THEN 2395 WRITE(6,*) ' IREO and IREO_INV ' 2396 CALL IWRTMA3(IREO,1,NDIM,1,NDIM) 2397 WRITE(6,*) 2398 CALL IWRTMA3(IREO_INV,1,NDIM,1,NDIM) 2399 END IF 2400* 2401 RETURN 2402 END 2403 SUBROUTINE COMB_TWO_REO(IREO3,IREO2,IREO1,NDIM) 2404* 2405* Two reorder arrays IREO1, IREO2 are given, combine these 2406* Obtain IREO3(I) = IREO2(IREO1(I)) 2407* 2408*. Jeppe Olsen, Oct. 2, 2012 2409*. Last modification; Oct. 2, 2012; Jeppe Olsen, original version 2410* 2411 INCLUDE 'implicit.inc' 2412*. Input 2413 INTEGER IREO1(NDIM), IREO2(NDIM) 2414*. Output 2415 INTEGER IREO3(NDIM) 2416* 2417 NTEST = 00 2418 IF(NTEST.GE.100) THEN 2419 WRITE(6,*) ' Output from COMB_TWO_REO ' 2420 END IF 2421 2422 DO I = 1, NDIM 2423 IREO3(I) = IREO2(IREO1(I)) 2424 END DO 2425* 2426 IF(NTEST.GE.100) THEN 2427 WRITE(6,*) ' Input: IREO1, IREO2 ' 2428 CALL IWRTMA3(IREO1,1,NDIM,1,NDIM) 2429 WRITE(6,*) 2430 CALL IWRTMA3(IREO2,1,NDIM,1,NDIM) 2431 WRITE(6,*) ' Output: IREO3 ' 2432 CALL IWRTMA3(IREO3,1,NDIM,1,NDIM) 2433 END IF 2434* 2435 RETURN 2436 END 2437 SUBROUTINE GET_IACT_TO_GENSM_REO(IACT_TO_GENSM_REO, 2438 & ISTA_TO_GENSM_REO, MO_STA_TO_ACT_REO, NTOOB) 2439* 2440* Obtain reorder array going from orbitals in super-symmetry order to 2441* actual supersymmetry order 2442* 2443*. Jeppe Olsen, Oct.2, 2012 2444*. Last revision; Oct. 3, 2012; Jeppe Olsen; debugged 2445* 2446 INCLUDE 'implicit.inc' 2447COLD INCLUDE 'mxpdim.inc' 2448*. Input 2449 INTEGER ISTA_TO_GENSM_REO(NTOOB),MO_STA_TO_ACT_REO(NTOOB) 2450*. Output 2451 INTEGER IACT_TO_GENSM_REO(NTOOB) 2452*. local scratch 2453COLD INTEGER IREO(MXPORB) 2454*. Array from actual order to standard order 2455COLD INVERT_REO(IREO,IREO_INV, NDIM) 2456COLD CALL INVERT_REO(MO_STA_TO_ACT_REO,IREO,NTOOB) 2457*. Combine IACT_TO_GENSM_REO(I) = MO_STA_TO_ACT_REO(ISTA_TO_GENSM_REO((I)) 2458C COMB_TWO_REO(IREO3,IREO2,IREO1,NDIM) 2459 CALL COMB_TWO_REO(IACT_TO_GENSM_REO, 2460 & MO_STA_TO_ACT_REO, ISTA_TO_GENSM_REO,NTOOB) 2461* 2462 NTEST = 100 2463 IF(NTEST.GE.100) THEN 2464 WRITE(6,*) ' IACT_TO_GENSM_REO array ' 2465 WRITE(6,*) ' ======================= ' 2466 CALL IWRTMA3(IACT_TO_GENSM_REO,1,NTOOB,1,NTOOB) 2467 END IF 2468* 2469 RETURN 2470 END 2471 SUBROUTINE REFORM_CMO_SUP_SHL(CMO_SUP,CMO_SHL,IWAY) 2472* 2473* Reform CMO matrix between supersymmetry order and shell order. 2474*. Outer routine 2475* 2476*. Jeppe Olsen, March 5, 2013 2477* 2478 INCLUDE 'implicit.inc' 2479 INCLUDE 'mxpdim.inc' 2480 INCLUDE 'wrkspc-static.inc' 2481 INCLUDE 'glbbas.inc' 2482#include "errquit.fh" 2483#include "mafdecls.fh" 2484#include "global.fh" 2485* 2486* IWAY = 1: Super-symmetry to shell order 2487* IWAY = 2: Shell-order to super-symmetry order 2488* 2489 CALL REFORM_CMO_SUP_SHL_IN(CMO_SUP,CMO_SHL,IWAY, 2490 & int_mb(KNSUPSYM_FOR_IRREP),int_mb(KIBSUPSYM_FOR_IRREP), 2491 & int_mb(KISUPSYM_FOR_IRREP)) 2492* 2493 RETURN 2494 END 2495 SUBROUTINE REFORM_CMO_SUP_SHL_IN(CMO_SUP,CMO_SHL,IWAY, 2496 & NSUPSYM_FOR_IRREP,IBSUPSYM_FOR_IRREP, 2497 & ISUPSYM_FOR_IRREP) 2498 2499* 2500* In shell-order the orbitals are arranged as 2501* 2502* Loop over irreducible representations 2503* Loop over shells of this supersymmetry 2504* Loop over subshells for a given shell 2505* End of loop over subshell of a given shell 2506* End of loop over shells 2507* End of loop over irreps 2508* 2509*. Jeppe Olsen, March 5, 2013 2510* 2511* IWAY = 1: Super-symmetry to shell order 2512* IWAY = 2: Shell-order to super-symmetry order 2513* 2514 INCLUDE 'implicit.inc' 2515 INCLUDE 'mxpdim.inc' 2516 INCLUDE 'orbinp.inc' 2517 INCLUDE 'lucinp.inc' 2518* 2519 INTEGER IBSUPSYM_FOR_IRREP(*),NSUPSYM_FOR_IRREP(*) 2520 INTEGER ISUPSYM_FOR_IRREP(*) 2521* 2522*. Input and output 2523 DIMENSION CMO_SUP(*), CMO_SHL(*) 2524* 2525 NTEST = 00 2526 IF(NTEST.GE.10) THEN 2527 WRITE(6,*) ' Info from REFORM_CMO_SUP_SHL' 2528 WRITE(6,*) ' ============================' 2529 WRITE(6,*) 2530 WRITE(6,*) ' IWAY = ', IWAY 2531C? WRITE(6,*) ' First 10 elements of CMO_SUP ' 2532C? CALL WRTMAT(CMO_SUP,1,10,1,10) 2533 END IF 2534* 2535 IB_SHL = 1 2536 DO IRREP = 1, N_SUPSYM_IRREP 2537*. Number of shells for this irrep 2538 NSHELL = NBAS_SUPSYM(IBSUPSYM_FOR_IRREP(IRREP)) 2539* 2540 IB = IBSUPSYM_FOR_IRREP(IRREP) 2541 NDEG = NSUPSYM_FOR_IRREP(IRREP) 2542 DO ISHELL = 1, NSHELL 2543 DO IISUPSYM = IB, IB + NDEG-1 2544 ISUPSYM = ISUPSYM_FOR_IRREP(IISUPSYM) 2545*. Offset to shell ISHELL of supersymmetry ISUPSYM 2546 IB_SUP = IB_CMOSUP_ORB(ISUPSYM,ISHELL) 2547 IF(IWAY.EQ.1) THEN 2548 CALL COPVEC(CMO_SUP(IB_SUP),CMO_SHL(IB_SHL),NSHELL) 2549C? WRITE(6,*) ' Elements copied ' 2550C? CALL WRTMAT(CMO_SHL(IB_SHL),1,NSHELL,1,NSHELL) 2551 ELSE 2552 CALL COPVEC(CMO_SHL(IB_SHL),CMO_SUP(IB_SUP),NSHELL) 2553 END IF 2554C? WRITE(6,'(A,3I4)') ' ISHELL IISUPSYM, ISHELL = ', 2555C? & ISHELL IISUPSYM, ISHELL 2556C? WRITE(6,'(A,2I4)') ' IB_SHL, IB_SUP = ', IB_SHL, IB_SUP 2557 IB_SHL = IB_SHL + NSHELL 2558 END DO ! loop over supersymmetries 2559 END DO ! End of loop over shells 2560 END DO ! End of loop over irreps 2561* 2562 IF(NTEST.GE.100) THEN 2563 WRITE(6,*) ' C arranged according to shells ' 2564 WRITE(6,*) ' ================================' 2565 WRITE(6,*) 2566 CALL PRINT_CSHELL(CMO_SHL) 2567 END IF 2568* 2569 RETURN 2570 END 2571 FUNCTION IB_CMOSUP_ORB(ISUPSYM,ISHELL) 2572* 2573* Determine offset of orbital with given supersymmetry and shell number 2574* in supersymmetry-ordered matrix 2575* 2576*. Jeppe Olsen, March 5, 2013 2577* 2578 INCLUDE 'implicit.inc' 2579 INCLUDE 'mxpdim.inc' 2580 INCLUDE 'orbinp.inc' 2581* 2582 NTEST = 00 2583*. Offset to supersymmetry block 2584 IOFF = 1 2585 DO JSUPSYM = 1, ISUPSYM-1 2586 NSHELL = NBAS_SUPSYM(JSUPSYM) 2587 IOFF = IOFF + NSHELL**2 2588 END DO 2589*. And to the given shell in the block 2590 IOFF = IOFF + (ISHELL-1)*NBAS_SUPSYM(ISUPSYM) 2591* 2592 IB_CMOSUP_ORB = IOFF 2593* 2594 IF(NTEST.GE.100) THEN 2595 WRITE(6,*) ' Info from IB_CMOSUP_ORB: ' 2596 WRITE(6,*) ' ISUPSYM, ISHELL, IB_CMOSUP_ORB = ', 2597 & ISUPSYM, ISHELL, IB_CMOSUP_ORB 2598 END IF 2599* 2600 RETURN 2601 END 2602 SUBROUTINE PRINT_CSHELL(CSHELL) 2603* 2604* A MO-AO expansion CSHELL is given in SHELL ordered form. 2605* Print it! 2606*. Outer part 2607*. Jeppe Olsen, March 5, 2013 2608* 2609 INCLUDE 'implicit.inc' 2610 INCLUDE 'mxpdim.inc' 2611 INCLUDE 'orbinp.inc' 2612 INCLUDE 'glbbas.inc' 2613 INCLUDE 'wrkspc-static.inc' 2614#include "errquit.fh" 2615#include "mafdecls.fh" 2616#include "global.fh" 2617*. Input 2618 DIMENSION CSHELL(*) 2619* 2620 CALL PRINT_CSHELLIN(CSHELL,int_mb(KNSUPSYM_FOR_IRREP), 2621 & int_mb(KIBSUPSYM_FOR_IRREP)) 2622* 2623 RETURN 2624 END 2625 SUBROUTINE PRINT_CSHELLIN(CSHELL,NSUPSYM_FOR_IRREP, 2626 & IBSUPSYM_FOR_IRREP) 2627* 2628* A MO-AO expansion CSHELL is given in SHELL ordered form. 2629* Print it! 2630*. Outer part 2631*. Jeppe Olsen, March 5, 2013 2632* 2633 INCLUDE 'implicit.inc' 2634 INCLUDE 'mxpdim.inc' 2635 INCLUDE 'orbinp.inc' 2636 INCLUDE 'glbbas.inc' 2637 INCLUDE 'wrkspc-static.inc' 2638* 2639 DIMENSION NSUPSYM_FOR_IRREP(*) 2640 DIMENSION IBSUPSYM_FOR_IRREP(*) 2641* 2642*. Jeppe Olsen, March 5, 2013 2643* 2644*. CMO to be printed 2645 DIMENSION CSHELL(*) 2646* 2647 IB = 1 2648 DO IRREP = 1, N_SUPSYM_IRREP 2649 WRITE(6,*) ' Irrep number ', IRREP 2650 WRITE(6,*) ' =====================' 2651 NDEG = NSUPSYM_FOR_IRREP(IRREP) 2652 NSHELL = NBAS_SUPSYM(IBSUPSYM_FOR_IRREP(IRREP)) 2653 DO ISHELL = 1, NSHELL 2654 WRITE(6,*) ' Subshells for shell ', ISHELL 2655 CALL WRTMAT(CSHELL(IB),NSHELL,NDEG,NSHELL,NDEG) 2656 IB = IB + NDEG*NSHELL 2657C? WRITE(6,*) ' TESTY, NDEG, NSHELL, IB = ', 2658C? & NDEG, NSHELL, IB 2659 END DO 2660 END DO 2661* 2662 RETURN 2663 END 2664 SUBROUTINE EXC_OO_TO_SS(NOOEX,IOOEX,NSSEX,ISSEX, 2665 & NOOFSSX,IBOOFSSX,IOOFSSX) 2666* 2667* A set of NOOEXC orbital excitations are given by IOOEXC. Obtain the 2668* corresponding shell excitations and obtain the mappings between the 2669* Orbital and shell excitations 2670* 2671*. Jeppe Olsen, March 5/6 2013 2672* 2673 INCLUDE 'implicit.inc' 2674 INCLUDE 'mxpdim.inc' 2675 INCLUDE 'glbbas.inc' 2676 INCLUDE 'wrkspc-static.inc' 2677*. Input 2678 RETURN 2679 END 2680 SUBROUTINE NONRED_SS_EXC(NOOEX,IOOEXC,NSSEX) 2681* 2682* A set of orbital excitation is given by IOOEXC 2683* Obtain the corresponding set of shell-excitations 2684* 2685*. Input: 2686* NOOEX: Number of orbital excitations 2687* IOOEX: Orbital excitations 2688*.Output: (mainly as pointers) 2689* NSSEX: Number of shell excitations 2690* KISSEXC: THe shell excitations in compact form 2691* KNIOOFSS: Number of orbital excitations for a gvien shell excitation 2692* KIBOOFSS: The offset for orbitial excitations 2693* KIOOFSS: The actual orbital excitations for a given shell excitation 2694* 2695* The information returned is pointers defined in the routine 2696* 2697* 2698*. Jeppe Olsen, March 6, 2013 2699* 2700 INCLUDE 'implicit.inc' 2701 INCLUDE 'mxpdim.inc' 2702 INCLUDE 'glbbas.inc' 2703 INCLUDE 'orbinp.inc' 2704 INCLUDE 'wrkspc-static.inc' 2705#include "errquit.fh" 2706#include "mafdecls.fh" 2707#include "global.fh" 2708*. Input 2709 INTEGER IOOEXC(2,NOOEX) 2710* 2711 IDUM = 0 2712*. No mark as allocated memory should be conserved 2713* 2714 CALL MEMMAN(KLSSEXE,NTOSH**2,'ADDL ',1,'SS_EXE') 2715 CALL MEMMAN(KLACT_TO_STA,NTOOB,'ADDL ',1,'REACST') 2716*. Obtain mapping from actual to standard mapping of orbitals from inverse 2717C INVERT_MAP(MAP,MAPINV,NELMNT) 2718 CALL INVERT_MAP(int_mb(KMO_STA_TO_ACT_REO),int_mb(KLACT_TO_STA), 2719 & NTOOB) 2720 2721* 2722*. Obtain the number of shell-shell excitations 2723* 2724C GET_SSEXC(NOOEX,IOOEXC,NSSEX,ISSEXE,ISSEXC, 2725C & NOOFSS,IBOOFSS,IOOFSS,IFLAG,IACT_TO_STA,ISHL_FOR_STA) 2726 CALL GET_SSEXC(NOOEX,IOOEXC,NSSEX,int_mb(KLSSEXE), 2727 & IDUM,IDUM,IDUM,IDUM,1, 2728 & int_mb(KLACT_TO_STA),int_mb(KISHELL_FOR_BAS)) 2729*. Allocate space for the shell excitations 2730 CALL MEMMAN(KISSEXC,2*NSSEX,'ADDL ',1,'ISSEXC') 2731 CALL MEMMAN(KNOOFSS,NSSEX,'ADDL ',1,'NOOFSS') 2732 CALL MEMMAN(KIBOOFSS,NSSEX,'ADDL ',1,'BOOFSS') 2733 CALL MEMMAN(KIOOFSS,NOOEX,'ADDL ',1,'IOOFSS') 2734* 2735*. And the actual shell-shell excitations 2736* 2737C GET_SSEXC(NOOEX,IOOEXC,NSSEX,ISSEXE,ISSEXC, 2738C & NOOFSS,IBOOFSS,IOOFSS,IFLAG,IACT_TO_STA,ISHL_FOR_STA) 2739 CALL GET_SSEXC(NOOEX,IOOEXC,NSSEX,int_mb(KLSSEXE), 2740 & int_mb(KISSEXC),int_mb(KNOOFSS),int_mb(KIBOOFSS), 2741 & int_mb(KIOOFSS),2, 2742 & int_mb(KLACT_TO_STA),int_mb(KISHELL_FOR_BAS)) 2743* 2744 RETURN 2745 END 2746 SUBROUTINE GET_SSEXC(NOOEX,IOOEXC,NSSEX,ISSEXE,ISSEXC, 2747 & NOOFSS,IBOOFSS,IOOFSS,IFLAG,IACT_TO_STA,ISHL_FOR_STA) 2748* 2749* IFLAG = 1: Obtain the allowed shell-shell excitation in expanded form in ISSEXE 2750* IFLAG = 2: Use the ISSEXE array to obtain the actual components of the 2751* shell-excitations in compact form 2752* 2753*. Jeppe Olsen, March 6, 2013 2754* 2755 INCLUDE 'implicit.inc' 2756 INCLUDE 'mxpdim.inc' 2757 INCLUDE 'orbinp.inc' 2758 INCLUDE 'wrkspc-static.inc' 2759*. Input 2760 INTEGER IOOEXC(2,NOOEX) 2761*. Mapping from actual(symmetry-type) order to standard order 2762 INTEGER IACT_TO_STA(*) 2763*. Shell number of given basis function/standard ordered MO 2764 INTEGER ISHL_FOR_STA(*) 2765*. Output 2766 INTEGER ISSEXE(NTOSH,NTOSH) 2767 INTEGER ISSEXC(2,*),NOOFSS(*),IBOOFSS(*),IOOFSS(*) 2768* 2769 NTEST = 00 2770 IF(NTEST.GE.100) THEN 2771 WRITE(6,*) ' Output from GET_SSEXC ' 2772 WRITE(6,*) ' ======================' 2773 WRITE(6,*) 2774 WRITE(6,*) ' IFLAG = ', IFLAG 2775 WRITE(6,*) ' NTOSH = ', NTOSH 2776 WRITE(6,*) ' NOOEX = ', NOOEX 2777 END IF 2778* 2779*. Set up array of shell-shell excitations 2780* 2781 IF(IFLAG.EQ.1) THEN 2782* 2783* Obtain the array ISSEXE(ISHL,JSHL) giving the number of 2784* orbital excitations between these shells. 2785* 2786*. INITIALIZE 2787 IZERO = 0 2788 CALL ISETVC(ISSEXE,IZERO,NTOSH**2) 2789 DO JOOEX = 1, NOOEX 2790*. Orbitals in symmetry order 2791 IIORB = IREOTS(IOOEXC(1,JOOEX)) 2792 JJORB = IREOTS(IOOEXC(2,JOOEX)) 2793C? WRITE(6,*) ' JOOEX, IIORB, JJORB = ', JOOEX, IIORB, JJORB 2794*. The same numbers in standard order 2795 IORB = IACT_TO_STA(IIORB) 2796 JORB = IACT_TO_STA(JJORB) 2797*. Shell numbers of these orbitals 2798 ISHL = ISHL_FOR_STA(IORB) 2799C? WRITE(6,*) ' IORB, ISHL = ', IORB, ISHL 2800 JSHL = ISHL_FOR_STA(JORB) 2801C? WRITE(6,*) ' JORB, JSHL = ', JORB, JSHL 2802*. Enroll 2803 ISSEXE(ISHL,JSHL) = ISSEXE(ISHL,JSHL) + 1 2804 END DO 2805*. Total number of shell-shell excitations 2806 NSSEX = 0 2807 DO ISHL = 1, NTOSH 2808 DO JSHL = 1, NTOSH 2809 IF(ISSEXE(ISHL,JSHL).NE.0) THEN 2810 NSSEX = NSSEX + 1 2811 END IF 2812 END DO 2813 END DO 2814* 2815 IF(NTEST.GE.1000) THEN 2816 WRITE(6,*) ' Shell-Shell excitation array ' 2817 CALL IWRTMA3(ISSEXE,NTOSH,NTOSH,NTOSH,NTOSH) 2818 END IF 2819* 2820 END IF ! IFLAG = 1 2821* 2822 IF(IFLAG.EQ.2) THEN 2823* 2824 IF(NTEST.GE.1000) THEN 2825 WRITE(6,*) ' Shell-Shell excitation array(input now) ' 2826 CALL IWRTMA3(ISSEXE,NTOSH,NTOSH,NTOSH,NTOSH) 2827 END IF 2828*. Obtain the shell-shell-excitations in compact form and pointer to start of OO for given SS 2829 ISSEX = 0 2830 DO ISHL = 1, NTOSH 2831 DO JSHL = 1, NTOSH 2832 IF(ISSEXE(ISHL,JSHL).NE.0) THEN 2833 ISSEX = ISSEX + 1 2834C? WRITE(6,*) ' ISHL, JSHL, ISSEX = ', ISHL, JSHL, ISSEX 2835 ISSEXC(1,ISSEX) = ISHL 2836 ISSEXC(2,ISSEX) = JSHL 2837 NOOFSS(ISSEX) = ISSEXE(ISHL,JSHL) 2838 END IF 2839 END DO 2840 END DO 2841C? WRITE(6,*) ' NOOFSS: ' 2842C? CALL IWRTMA3(NOOFSS,1,ISSEX,1,ISSEX) 2843*. Pointers to start of Orbital excitations for given shell excitation 2844 IBS = 1 2845 DO ISSEX = 1, NSSEX 2846 IBOOFSS(ISSEX) = IBS 2847 IBS = IBS + NOOFSS(ISSEX) 2848*. And clear for later use 2849 NOOFSS(ISSEX) = 0 2850 END DO 2851C? WRITE(6,*) ' IBOOFSS: ' 2852C? CALL IWRTMA3(IBOOFSS,1,ISSEX,1,ISSEX) 2853*. Change ISSEXE to give the number of a given shell excitation 2854 ISSEX = 0 2855 DO ISHL = 1, NTOSH 2856 DO JSHL = 1, NTOSH 2857 IF(ISSEXE(ISHL,JSHL).NE.0) THEN 2858 ISSEX = ISSEX + 1 2859 ISSEXE(ISHL,JSHL) = ISSEX 2860 END IF 2861 END DO 2862 END DO 2863*. And the orbital excitations of a given shell excitation 2864 DO IOOEX = 1, NOOEX 2865*. Orbitals in symmetry-type order 2866 IIORB = IREOTS(IOOEXC(1,IOOEX)) 2867 JJORB = IREOTS(IOOEXC(2,IOOEX)) 2868*. The same numbers in standard order 2869 IORB = IACT_TO_STA(IIORB) 2870 JORB = IACT_TO_STA(JJORB) 2871*. Shell numbers of these orbitals 2872 ISHL = ISHL_FOR_STA(IORB) 2873 JSHL = ISHL_FOR_STA(JORB) 2874 IJEXC = ISSEXE(ISHL,JSHL) 2875 NOOFSS(IJEXC) = NOOFSS(IJEXC) + 1 2876 IB = IBOOFSS(IJEXC) 2877 IOOFSS(IB + NOOFSS(IJEXC)-1) = IOOEX 2878 END DO 2879* 2880 IF(NTEST.GE.100) THEN 2881 WRITE(6,*) ' The orbital excitations of shell excitations' 2882 WRITE(6,*) ' ============================================' 2883 DO ISSEX = 1, NSSEX 2884 WRITE(6,*) 2885 WRITE(6,*) ' Shell excitation ', ISSEX 2886 IB = IBOOFSS(ISSEX) 2887 N = NOOFSS(ISSEX) 2888 CALL IWRTMA3(IOOFSS(IB),1,N,1,N) 2889 END DO 2890 END IF ! NTEST large 2891 END IF ! IFLAG = 2 2892* 2893 2894* 2895 RETURN 2896 END 2897 SUBROUTINE INVERT_MAP(MAP,MAPINV,NELMNT) 2898* 2899* A mappping MAP(I) is given. Obtain inverse mapping 2900* 2901*. Jeppe Olsen, March 6, 2013 2902 INCLUDE 'implicit.inc' 2903*. Input 2904 INTEGER MAP(NELMNT) 2905*. Output 2906 INTEGER MAPINV(NELMNT) 2907* 2908 DO I = 1, NELMNT 2909 MAPINV(MAP(I)) = I 2910 END DO 2911* 2912 NTEST = 00 2913 IF(NTEST.GE.100) THEN 2914 WRITE(6,*) ' Map (input) and inverted map(output)' 2915 WRITE(6,*) ' =====================================' 2916 WRITE(6,*) 2917 CALL IWRTMA3(MAP,1,NELMNT,1,NELMNT) 2918 CALL IWRTMA3(MAPINV,1,NELMNT,1,NELMNT) 2919 END IF 2920* 2921 RETURN 2922 END 2923 SUBROUTINE SHELL_AVERAGE_ORBEXC(VECIN,NSSEX,NOOFSS,IBOOFSS, 2924 & IOOFSS,VECUT,NOOEX,ICOPY) 2925* 2926* A vector V over orbital excitation is given. 2927* Average over orbital excitations belonging to identical shell excitation. 2928* 2929* If ICOPY = 1, the output vector is copied to the input vector 2930* 2931*. Jeppe Olsen, March 7, 2013 2932* 2933 INCLUDE 'implicit.inc' 2934*. General input 2935 INTEGER NOOFSS(NSSEX),IBOOFSS(NSSEX),IOOFSS(*) 2936*. Specific input 2937 DIMENSION VECIN(*) 2938*. Output 2939 DIMENSION VECUT(*) 2940* 2941 NTEST = 00 2942* 2943 DO ISSEX = 1, NSSEX 2944*. Obtain average value for this shell-excitation 2945 IB = IBOOFSS(ISSEX) 2946 N = NOOFSS(ISSEX) 2947 AVE = 0.0D0 2948 DO I = 1, N 2949 AVE = AVE + VECIN(IOOFSS(IB-1+I)) 2950 END DO 2951 AVE = AVE/FLOAT(N) 2952*. And spread out 2953 DO I = 1, N 2954 VECUT(IOOFSS(IB-1+I)) = AVE 2955 END DO 2956 END DO 2957* 2958 IF(NTEST.GE.100) THEN 2959 WRITE(6,*) 2960 & ' Input and output from averaging over shell-components ' 2961 WRITE(6,*) 2962 & ' ======================================================' 2963 CALL WRT_2VEC(VECIN,VECUT,NOOEX) 2964 END IF 2965* 2966 IF(ICOPY.NE.0) THEN 2967 CALL COPVEC(VECUT,VECIN,NOOEX) 2968 END IF 2969* 2970 RETURN 2971 END 2972 SUBROUTINE AVE_DENS_OVER_SUBSHELLS(RHO1,RHO1AVE) 2973* 2974* Average the one-particle density over subshells belonging to a given shell 2975* Outer routine 2976* 2977*. Jeppe Olsen, March 8, 2013 2978* 2979 INCLUDE 'implicit.inc' 2980 INCLUDE 'mxpdim.inc' 2981 INCLUDE 'wrkspc-static.inc' 2982 INCLUDE 'orbinp.inc' 2983 INCLUDE 'glbbas.inc' 2984 INCLUDE 'cgas.inc' 2985#include "errquit.fh" 2986#include "mafdecls.fh" 2987#include "global.fh" 2988*. Input 2989 DIMENSION RHO1(NACOB,NACOB) 2990*. Output 2991 DIMENSION RHO1AVE(NACOB,NACOB) 2992* 2993 CALL AVE_DENS_OVER_SUBSHELLS_IN(RHO1,RHO1VE,NSHTO,NACOB, 2994 & int_mb(KNBAS_FOR_SHELL),int_mb(KIBBAS_FOR_SHELL), 2995 & int_mb(KIBAS_FOR_SHELL),ITPFSO,NGAS) 2996* 2997 RETURN 2998 END 2999 SUBROUTINE AVE_DENS_OVER_SUBSHELLS_IN(RHO1,RHO1AVE,NSHTO,NACOB, 3000 & NBAS_FOR_SHELL,IBBAS_FOR_SHELL,IBAS_FOR_SHELL, 3001 & ITPFSO,NGAS) 3002* 3003* Average the one-particle density over subshells belonging to a given shell 3004* 3005*. We do not have a list of active shells. What we instead will do is to 3006* connect a shell with the first subshell of this shell, and then we can check 3007* the gas space of this subshell 3008* 3009* 3010 INCLUDE 'implicit.inc' 3011 INTEGER NBAS_FOR_SHELL(*), IBBAS_FOR_SHELL(*) 3012 INTEGER IBAS_FOR_SHELL(*), ITPFSO(*) 3013*. Input 3014 DIMENSION RHO1(NACOB,NACOB) 3015*. Output 3016 DIMENSION RHO1AVE(NACOB,NACOB) 3017* 3018 NTEST = 10 3019 IF(NTEST.GE.10) THEN 3020 WRITE(6,*) ' Info from AVE_DENS_OVER_SUBSHELLS_IN ' 3021 END IF 3022 WRITE(6,*) ' THIS ROUTINE HAS NEVER BEEN TESTED OR USED ' 3023* 3024 DO ISHELL = 1, NSHTO 3025 IB = IBBAS_FOR_SHELL(ISHELL) 3026 IORB = IBAS_FOR_SHELL(IB) 3027 IF(0.LT.ITPFSO(IORB).AND.ITPFSO(IORB).LE.NGAS) THEN 3028*. Shell is active 3029 DO JSHELL = 1, NSHTO 3030 JB = IBBAS_FOR_SHELL(JSHELL) 3031 JORB = IBAS_FOR_SHELL(JB) 3032 IF(0.LT.ITPFSO(JORB).AND.ITPFSO(JORB).LE.NGAS) THEN 3033* 3034* Obtain average value for these shells, non-diagonal and diagonal 3035* 3036 AVE = 0.0D0 3037 AVED = 0.0D0 3038* 3039 NI = NBAS_FOR_SHELL(ISHELL) 3040 NJ = NBAS_FOR_SHELL(JSHELL) 3041 DO IIORB = IB, IB - 1 + NI 3042 DO JJORB = JB, JB - 1 + NI 3043 IORB = IBAS_FOR_SHELL(IIORB) 3044 JORB = IBAS_FOR_SHELL(JJORB) 3045* 3046 IF(IORB.NE.JORB) THEN 3047 AVE = AVE + RHO1(IORB,JORB) 3048 ELSE 3049 AVED = AVED + RHO1(IORB,JORB) 3050 END IF 3051 END DO 3052 END DO 3053* 3054 IF(ISHELL.EQ.JSHELL) THEN 3055 AVED = AVED/FLOAT(NI) 3056 AVE = AVE/(FLOAT(NI)*(FLOAT(NI)-1)) 3057 ELSE 3058 AVE = AVE/(FLOAT(NI)*(FLOAT(NI))) 3059 END IF 3060* 3061* Scatter average values out 3062* 3063 DO IIORB = IB, IB - 1 + NI 3064 DO JJORB = JB, JB - 1 + NI 3065 IORB = IBAS_FOR_SHELL(IIORB) 3066 JORB = IBAS_FOR_SHELL(JJORB) 3067 RHO1AVE(IORB,JORB) = AVE 3068 END DO 3069 END DO 3070* 3071 IF(ISHELL.EQ.JSHELL) THEN 3072 DO IIORB = IB, IB-1+NI 3073 IORB = IBAS_FOR_SHELL(IIORB) 3074 RHO1AVE(IORB,IORB) = AVED 3075 END DO 3076 END IF 3077* 3078 END IF 3079 END DO ! For jshell 3080 END IF 3081 END DO ! For Ishell 3082* 3083 IF(NTEST.GE.10) THEN 3084 WRITE(6,*) ' Original and averaged density matrix ' 3085 WRITE(6,*) ' =====================================' 3086 WRITE(6,*) 3087 CALL WRTMAT(RHO1,NACOB,NACOB,NACOB,NACOB) 3088 WRITE(6,*) 3089 CALL WRTMAT(RHO1AVE,NACOB,NACOB,NACOB,NACOB) 3090 END IF 3091* 3092 RETURN 3093 END 3094 SUBROUTINE AVE_SUPSYM_MAT(ASUP,NOBPSPSM,IPACK) 3095* 3096* A matrix ASUP over blocks of supersymmetry is given 3097* Average over blocks belonging to the same irrep. 3098* 3099*. Outer part 3100* 3101*. Jeppe Olsen, March 9, 2013 3102* 3103 INCLUDE 'implicit.inc' 3104 INCLUDE 'mxpdim.inc' 3105 INCLUDE 'orbinp.inc' 3106 INCLUDE 'glbbas.inc' 3107 INCLUDE 'wrkspc-static.inc' 3108#include "errquit.fh" 3109#include "mafdecls.fh" 3110#include "global.fh" 3111* 3112*. Input 3113 DIMENSION NOBPSPSM(*) 3114*. Input and output 3115 DIMENSION ASUP(*) 3116* 3117 CALL AVE_SUPSYM_MAT_IN(ASUP,NOBPSPSM,IPACK, 3118 & N_SUPSYM_IRREP,N_SUPSYM,int_mb(KNSUPSYM_FOR_IRREP), 3119 & int_mb(KIBSUPSYM_FOR_IRREP),int_mb(KISUPSYM_FOR_IRREP)) 3120* 3121 RETURN 3122 END 3123 SUBROUTINE AVE_SUPSYM_MAT_IN(ASUP,NOBPSPSM,IPACK, 3124 & N_SUPSYM_IRREP,N_SUPSYM,NSUPSYM_FOR_IRREP, 3125 & IBSUPSYM_FOR_IRREP,ISUPSYM_FOR_IRREP) 3126* 3127* A matrix ASUP over supersymmmetries is given with NOBPSPSM orbitals per supersymmetry 3128* Average over supersymmetries belonging to the same super-symmetry irrep 3129* 3130*. Jeppe Olsen, March 9, 2013 3131* 3132*. General input 3133 INCLUDE 'implicit.inc' 3134 INTEGER NSUPSYM_FOR_IRREP(*), IBSUPSYM_FOR_IRREP(*) 3135 INTEGER ISUPSYM_FOR_IRREP(*) 3136*. 3137 INTEGER NOBPSPSM(*) 3138*. Input and output 3139 DIMENSION ASUP(*) 3140* 3141 NTEST = 100 3142 ONE = 1.0D0 3143* 3144 IF(NTEST.GE.100) THEN 3145 WRITE(6,*) ' Info from AVE_SUPSYM_MAT_IN ' 3146 WRITE(6,*) ' ============================' 3147 WRITE(6,*) 3148 WRITE(6,*) ' N_SUPSYM = ', N_SUPSYM 3149 WRITE(6,*) ' NOBPSPSM: ' 3150 CALL IWRTMA(NOBPSPSM,1,N_SUPSYM,1,N_SUPSYM) 3151 WRITE(6,*) ' Input matrix over supersymmetries ' 3152 WRITE(6,*) 3153 CALL APRBLM2(ASUP,NOBPSPSM,NOBPSPSM,N_SUPSYM,IPACK) 3154 END IF 3155 3156* 3157 DO IRREP = 1, N_SUPSYM_IRREP 3158 NSPSM = NSUPSYM_FOR_IRREP(IRREP) 3159 IF(NTEST.GE.1000) WRITE(6,*) ' Info for IRREP = ', IRREP 3160 IF(NSPSM.GT.1) THEN 3161* 3162*. Average over the various components in the first supersymmetry of the given irrep 3163* 3164 IB = IBSUPSYM_FOR_IRREP(IRREP) 3165 N = NSUPSYM_FOR_IRREP(IRREP) 3166 ISUPSYM1 = ISUPSYM_FOR_IRREP(IB) 3167 L1 = NOBPSPSM(ISUPSYM1) 3168 IF(IPACK.EQ.0) THEN 3169 LBLK = L1*L1 3170 ELSE 3171 LBLK = L1*(L1+1)/2 3172 END IF 3173 IOFF1 = IOFF_BLCK(ISUPSYM1,NOBPSPSM,NOBPSPSM,IPAK) 3174 WRITE(6,*) ' ISUPSYM1, IOFF1 = ', ISUPSYM1, IOFF1 3175*. And terms from the remaining supersymmetries 3176 DO IISUPSYM = 2, N 3177 ISUPSYM = ISUPSYM_FOR_IRREP(IB-1+IISUPSYM) 3178 IOFF = IOFF_BLCK(ISUPSYM,NOBPSPSM,NOBPSPSM,IPAK) 3179 CALL VECSUM(ASUP(IOFF1),ASUP(IOFF1),ASUP(IOFF),ONE,ONE,LBLK) 3180 END DO 3181 FACTOR = 1.0D0/FLOAT(N) 3182C? WRITE(6,*) ' FACTOR = ', FACTOR 3183 CALL SCALVE(ASUP(IOFF1),FACTOR,LBLK) 3184 IF(NTEST.GE.1000) THEN 3185 WRITE(6,*) ' Averaged block ' 3186 CALL APRBLM2(ASUP(IOFF1),L1,L1,1,IPACK) 3187 END IF 3188* 3189*. And copy the average to the remaining blocks 3190* 3191 DO IISUPSYM = 2, N 3192 ISUPSYM = ISUPSYM_FOR_IRREP(IB-1+IISUPSYM) 3193 IOFF = IOFF_BLCK(ISUPSYM,NOBPSPSM,NOBPSPSM,IPAK) 3194 CALL COPVEC(ASUP(IOFF1),ASUP(IOFF),LBLK) 3195 END DO 3196* 3197 END IF ! Irrep was degenerate 3198 END DO ! Loop over irreps 3199* 3200 IF(NTEST.GE.100) THEN 3201 WRITE(6,*) ' Averaged matrix ' 3202 CALL APRBLM2(ASUP,NOBPSPSM,NOBPSPSM,N_SUPSYM,IPACK) 3203 END IF 3204* 3205 RETURN 3206 END 3207 FUNCTION IOFF_BLCK(IBLK,LR,LC,IPAK) 3208* 3209* Offset to block IBLK in matrix with LR/LC row/colomn elements per block 3210* 3211* Jeppe Olsen, March 9, (did it really take me about 25 years to write this function) 3212* 3213 INCLUDE 'implicit.inc' 3214* 3215 INTEGER LR(*), LC(*) 3216* 3217 IOFF = 1 3218 DO JBLK = 1, IBLK - 1 3219 IF(IPAK.EQ.0) THEN 3220 LBLK = LR(JBLK)*LC(JBLK) 3221 ELSE 3222 LBLK = LR(JBLK)*(LR(JBLK)+1)/2 3223 END IF 3224 IOFF = IOFF + LBLK 3225 END DO 3226* 3227 IOFF_BLCK = IOFF 3228* 3229 RETURN 3230 END 3231 SUBROUTINE PRINT_CMO_AS_SHELLS(CMO,IFORM) 3232* 3233* A CMOA matrix CMO is given in form defined by CMO. 3234* Print as subshells of a given shell 3235* 3236* IFORM = 1: Input CMO is in standard order 3237* IFORM = 2: Input CMO is in actual(gas) order 3238* IFORM = 3: Input CMO is in supersymmetry-order 3239* IFORM = 4: Input CMO is in Shell-order 3240* 3241*. Jeppe Olsen, March 2013 3242* 3243 INCLUDE 'implicit.inc' 3244 INCLUDE 'mxpdim.inc' 3245 INCLUDE 'wrkspc-static.inc' 3246 INCLUDE 'glbbas.inc' 3247 INCLUDE 'orbinp.inc' 3248 INCLUDE 'lucinp.inc' 3249#include "errquit.fh" 3250#include "mafdecls.fh" 3251#include "global.fh" 3252*. Input 3253 DIMENSION CMO(*) 3254* 3255 IDUM = 0 3256 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'PR_SHL') 3257* 3258 LCMO = LEN_BLMAT(NSMOB,NTOOBS,NTOOBS,0) 3259 CALL MEMMAN(KLMO1,LCMO,'ADDL ',2,'LMO1 ') 3260* 3261*. Reform to shell format 3262* 3263 CALL REFORM_CMO(CMO,IFORM,dbl_mb(KLMO1),4) 3264* 3265*. And print 3266* 3267 CALL PRINT_CSHELL(dbl_mb(KLMO1)) 3268* 3269 3270 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'PR_SHL') 3271* 3272 RETURN 3273 END 3274 SUBROUTINE REFORM_CMO(C_IN,IFORM_IN, C_OUT, IFORM_OUT) 3275* 3276* Reform matrix from FORM IFORM_IN to FORM IFORM_OUT 3277* 3278* IFORM = 1 => standard form 3279* IFORM = 2 => gas-ordered form 3280* IFORM = 3 => super-symmetry form 3281* IFORM = 4 => shell form 3282* 3283*. Jeppe Olsen, March 2013 3284* 3285 INCLUDE 'implicit.inc' 3286 INCLUDE 'mxpdim.inc' 3287 INCLUDE 'lucinp.inc' 3288 INCLUDE 'orbinp.inc' 3289 INCLUDE 'wrkspc-static.inc' 3290#include "errquit.fh" 3291#include "mafdecls.fh" 3292#include "global.fh" 3293*. Input 3294 DIMENSION C_IN(*) 3295*. Output 3296 DIMENSION C_OUT(*) 3297* 3298 NTEST = 000 3299* 3300 IF(NTEST.GE.10) THEN 3301 WRITE(6,*) ' Info from REFORM_CMO ' 3302 WRITE(6,*) ' ==================== ' 3303 WRITE(6,*) 3304 WRITE(6,*) ' IFORM_IN,IFORM_OUT = ', IFORM_IN, IFORM_OUT 3305 END IF 3306* 3307 IDUM = 0 3308 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'REFCMO') 3309 LCMO = LEN_BLMAT(NSMOB,NTOOBS,NTOOBS,0) 3310 CALL MEMMAN(KLMO1,LCMO,'ADDL ',2,'LMO1 ') 3311* 3312*. Reform from input form to standard form 3313* 3314C REFORM_CMO_TO_STANDARD(CMO,CMOST,IFORM,IWAY) 3315 CALL REFORM_CMO_TO_STANDARD(C_IN,dbl_mb(KLMO1),IFORM_IN,1) 3316* 3317 IF(NTEST.GE.1000) THEN 3318 WRITE(6,*) ' Intermediate CMO in standard form ' 3319 CALL APRBLM2(dbl_mb(KLMO1),NTOOBS,NTOOBS,NSMOB,0) 3320 END IF 3321* 3322*. Reform from standard to output form 3323* 3324 CALL REFORM_CMO_TO_STANDARD(C_OUT,dbl_mb(KLMO1),IFORM_OUT,2) 3325* 3326 IF(NTEST.GE.100) THEN 3327 WRITE(6,*) ' Input CMO to REFORM_CMO: ' 3328 CALL PRINT_CMO_ARBFORM(C_IN,IFORM_IN) 3329 WRITE(6,*) ' Output CMO from REFORM_CMO: ' 3330 CALL PRINT_CMO_ARBFORM(C_OUT,IFORM_OUT) 3331 END IF 3332* 3333 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'REFCMO') 3334* 3335 RETURN 3336 END 3337 SUBROUTINE REFORM_CMO_TO_STANDARD(CMO,CMOST,IFORM,IWAY) 3338* 3339* Reform between various super-symmetry forms of a CMO matrix. 3340* 3341* A CMO matrix, CMO is given in a from specified by IFORM. Reform this to 3342* standard form and save in CMOST 3343* 3344* IFORM = 1 => CMO is in standard form 3345* IFORM = 2 => CMO is in gas-ordered form 3346* IFORM = 3 => CMO is in super-symmetry form 3347* IFORM = 4 => CMO is in shell form 3348* 3349* IWAY = 1: From the general form in CMO to CMOST 3350* IWAY = 2: From the standard form in CMOST to CMO 3351* 3352*. Jeppe Olsen, March 9, 2013 3353* 3354 INCLUDE 'implicit.inc' 3355 INCLUDE 'mxpdim.inc' 3356 INCLUDE 'glbbas.inc' 3357 INCLUDE 'orbinp.inc' 3358 INCLUDE 'wrkspc-static.inc' 3359 INCLUDE 'lucinp.inc' 3360#include "errquit.fh" 3361#include "mafdecls.fh" 3362#include "global.fh" 3363*. Input or output 3364 DIMENSION CMO(*), CMOST(*) 3365* 3366 NTEST = 00 3367 IF(NTEST.GE.100) THEN 3368 WRITE(6,*) 3369 WRITE(6,*) ' Info from REFORM_CMO_TO_STANDARD' 3370 WRITE(6,*) ' ================================' 3371 WRITE(6,*) 3372 WRITE(6,*) ' Iway and Iform = ', IWAY, IFORM 3373 WRITE(6,*) ' Input matrix: ' 3374 IF(IWAY.EQ.1) THEN 3375C PRINT_CMO_ARBFORM(CMO,IFORM) 3376 CALL PRINT_CMO_ARBFORM(CMO,IFORM) 3377 ELSE 3378 CALL PRINT_CMO_ARBFORM(CMOST,1) 3379 END IF 3380* 3381 END IF !Ntest is large 3382* 3383 IDUM = 0 3384 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'RFSPCM') 3385 LCMO = LEN_BLMAT(NSMOB,NTOOBS,NTOOBS,0) 3386 CALL MEMMAN(KLMO1,LCMO,'ADDL ',2,'LMO1 ') 3387* 3388 IF(IWAY.EQ.1) THEN 3389* 3390* CMO in some form => CMOST in standard form 3391* 3392* 3393 IF(IFORM.EQ.4) THEN 3394*. Shell in CMO => supersymmetry form in KLMO1 3395C REFORM_CMO_SUP_SHL(CMO_SUP,CMO_SHL,IWAY) 3396 CALL REFORM_CMO_SUP_SHL(dbl_mb(KLMO1),CMO,2) 3397 END IF 3398* 3399 IF (IFORM.GE.3) THEN 3400 IF(IFORM.EQ.3) CALL COPVEC(CMO,dbl_mb(KLMO1),LCMO) 3401*. Supersymmetry in KLMO1 => standard in CMOST 3402C REFORM_CMO_STA_GEN(CMO_STA,CMO_GEN,IDO_REORDER,IREO,IWAY) 3403 CALL REFORM_CMO_STA_GEN(CMOST,dbl_mb(KLMO1),0,0,2) 3404 ELSE IF (IFORM.EQ.2) THEN 3405*. Actual => standard 3406C REO_CMOAO(CIN,COUT,IREO,ICOPY,IWAY) 3407 CALL REO_CMOAO(CMO,CMOST, 3408 & int_mb(KMO_STA_TO_ACT_REO),0,2) 3409 ELSE IF(IFORM.EQ.1) THEN 3410*. Easy living, just copy 3411 CALL COPVEC(CMO,CMOST,LCMO) 3412 END IF 3413 ELSE IF (IWAY.EQ.2) THEN 3414* 3415* CMOST in standard form => CMO in some form 3416* 3417 IF(IFORM.GE.3) THEN 3418* Standard in CMOST to supersymmetry in CMO 3419 CALL REFORM_CMO_STA_GEN(CMOST,CMO,0,0,1) 3420 END IF 3421 IF(IFORM.EQ.4) THEN 3422 CALL COPVEC(CMO,dbl_mb(KLMO1),LCMO) 3423*. Supersymmetry to shell 3424C REFORM_CMO_SUP_SHL(CMO_SUP,CMO_SHL,IWAY) 3425 CALL REFORM_CMO_SUP_SHL(dbl_mb(KLMO1),CMO,1) 3426 END IF 3427* 3428 IF(IFORM.EQ.2) THEN 3429*. Standard => Actual/gas ordered 3430C REO_CMOAO(CIN,COUT,IREO,ICOPY,IWAY) 3431 CALL REO_CMOAO(CMOST,CMO, 3432 & int_mb(KMO_STA_TO_ACT_REO),0,1) 3433 ELSE IF(IFORM.EQ.1) THEN 3434*. Standard => standard 3435 CALL COPVEC(CMOST,CMO,LCMO) 3436 END IF 3437 END IF ! switch IWAY 3438* 3439 IF(NTEST.GE.100) THEN 3440 WRITE(6,*) ' Output matrix from REFORM_CMO_TO_STANDARD: ' 3441 IF(IWAY.EQ.2) THEN 3442 WRITE(6,*) ' Output matrix is of general type' 3443C PRINT_CMO_ARBFORM(CMO,IFORM) 3444 CALL PRINT_CMO_ARBFORM(CMO,IFORM) 3445 ELSE 3446 WRITE(6,*) ' Output matrix is STANDARD type ' 3447 CALL PRINT_CMO_ARBFORM(CMOST,1) 3448 END IF 3449 END IF 3450* 3451 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'RFSPCM') 3452* 3453 RETURN 3454 END 3455 SUBROUTINE PRINT_CMO_ARBFORM(CMO,IFORM) 3456* 3457* Print a CMO matrix with form defined by IFORM 3458* 3459*. IFORM = 1: Standard form 3460*. IFORM = 2: Actual/gas ordered 3461*. IFORM = 3: Supersymmetry ordered 3462*. IFORM = 4: Shell ordered 3463* 3464*. Jeppe Olsen, March 10, 2013 3465* 3466 INCLUDE 'implicit.inc' 3467 INCLUDE 'mxpdim.inc' 3468 INCLUDE 'orbinp.inc' 3469 INCLUDE 'lucinp.inc' 3470*. Input 3471 DIMENSION CMO(*) 3472* 3473 IF(IFORM.EQ.1.OR.IFORM.EQ.2) THEN 3474 CALL PRINT_CMOAO(CMO) 3475 ELSE IF(IFORM.EQ.3) THEN 3476 CALL APRBLM2(CMO,NBAS_SUPSYM,NBAS_SUPSYM,N_SUPSYM) 3477 ELSE IF(IFORM.EQ.4) THEN 3478 CALL PRINT_CSHELL(CMO) 3479 ELSE 3480 WRITE(6,*) ' PRINT_CMO_ARBFORM: Unknown IFORM = ', IFORM 3481 STOP ' PRINT_CMO_ARBFORM: Unknown IFORM ' 3482 END IF 3483* 3484 RETURN 3485 END 3486 SUBROUTINE REO_2SUPSYM_ORDERS(ISUPSYM1,ISUPSYM2,IREO12) 3487* Two sequences of orbitals defined by their supersymmetries are given 3488* 3489* Obtain reordering array IREO12(I1) = I2 3490* 3491* Jeppe Olsen, March 2013 3492* 3493 INCLUDE 'implicit.inc' 3494 INCLUDE 'mxpdim.inc' 3495 INCLUDE 'orbinp.inc' 3496 INCLUDE 'lucinp.inc' 3497*. Input 3498 INTEGER ISUPSYM1(NTOOB),ISUPSYM2(NTOOB) 3499*. Output 3500 INTEGER IREO12(NTOOB) 3501* 3502 NTEST = 1000 3503 IF(NTEST.GE.100) THEN 3504 WRITE(6,*) ' Output REO_2SUPSYM_ORDERS' 3505 WRITE(6,*) ' Supersymmetry lists 1 and 2 ' 3506 CALL IWRTMA3(ISUPSYM1,1,NTOOB,1,NTOOB) 3507 CALL IWRTMA3(ISUPSYM2,1,NTOOB,1,NTOOB) 3508 END IF 3509* 3510 DO ISUPSYM = 1, N_SUPSYM 3511 NORB1 = 0 3512 DO IORB1 = 1, NTOOB 3513 IF(ISUPSYM1(IORB1).EQ.ISUPSYM) THEN 3514 NORB1 = NORB1 + 1 3515*. Find orbital NORB1 of symmetry ISUPSYM in ISUPSYM2 3516 NORB2 = 0 3517 DO IIORB2 = 1, NTOOB 3518 IF(ISUPSYM2(IIORB2).EQ.ISUPSYM) THEN 3519 NORB2 = NORB2 + 1 3520 IF(NORB2.EQ.NORB1) THEN 3521 IORB2 = IIORB2 3522 END IF 3523 END IF 3524 END DO 3525 IREO12(IORB1) = IORB2 3526 END IF 3527 END DO 3528 END DO 3529* 3530 IF(NTEST.GE.100) THEN 3531 WRITE(6,*) ' Obtained redorder array ' 3532 CALL IWRTMA3(IREO12,1,NTOOB,1,NTOOB) 3533 END IF 3534* 3535 RETURN 3536 END 3537 SUBROUTINE GET_OCC_ORDER_SUPSYM(IMO_OCCORD_SUPSYM) 3538* 3539* Obtain the super-symmetry for the required order of orbitals 3540* in ocupation/gas order 3541* 3542*. Jeppe Olsen, March 2013 3543* 3544* 3545*. 3546 INCLUDE 'implicit.inc' 3547 INCLUDE 'mxpdim.inc' 3548 INCLUDE 'cgas.inc' 3549 INCLUDE 'lucinp.inc' 3550 INCLUDE 'orbinp.inc' 3551*. Output 3552 INTEGER IMO_OCCORD_SUPSYM(NTOOB) 3553* 3554 NTEST = 100 3555 IF(NTEST.GE.100) THEN 3556 WRITE(6,*) ' Output from GET_OCC_ORDER_SUPSYM ' 3557 WRITE(6,*) ' =================================' 3558 END IF 3559 IF(NTEST.GE.1000) THEN 3560 WRITE(6,*) ' NGAS_SUPSYM: ' 3561 CALL IWRTMA3(NGAS_SUPSYM,N_SUPSYM,NGAS+2,MXP_NSUPSYM,NGAS+2) 3562 END IF 3563* 3564 IZERO = 0 3565 CALL ISETVC(IMO_OCCORD_SUPSYM,IZERO,NTOOB) 3566* 3567 IORB = 0 3568 DO ISYM = 1, NSMOB 3569 NSUP = NSUP_FOR_STA_SYM(ISYM) 3570 IBSUP = IBSUP_FOR_STA_SYM(ISYM) 3571 IF(NTEST.GE.1000) WRITE(6,*) ' ISYM, NSUP = ', ISYM, NSUP 3572*. Loop over the spac division for this symmetry 3573 DO JGAS= 0, NGAS + 1 3574 DO IISUPSYM = IBSUP, IBSUP + NSUP -1 3575 ISUPSYM = ISUP_FOR_STA_SYM(IISUPSYM) 3576 NORB = NGAS_SUPSYM(ISUPSYM,JGAS) 3577 IF(NTEST.GE.1000) THEN 3578 WRITE(6,*) ' ISUPSYM, JGAS, NORB = ', 3579 & ISUPSYM, JGAS, NORB 3580 END IF 3581 DO IIORB = 1, NORB 3582 IORB = IORB + 1 3583 IMO_OCCORD_SUPSYM(IORB) = ISUPSYM 3584 IF(NTEST.GE.10000) 3585 & WRITE(6,*) ' IORB, ISUPSYM ', IORB, ISUPSYM 3586 END DO! loop over orbital IIORB 3587 END DO! Loop over supersymmetries for given standard 3588 END DO ! Loop over orbital spaces 3589 END DO ! Loop over standard symmetries 3590* 3591 IF(NTEST.GE.100) THEN 3592 WRITE(6,*) 3593 & ' GET_OCC_ORDER: Supersymmetry of occ-ordered orbitals' 3594 CALL IWRTMA3(IMO_OCCORD_SUPSYM,1,NTOOB,1,NTOOB) 3595 END IF 3596* 3597 RETURN 3598 END 3599 SUBROUTINE ANA_SUBSHELLS_CMO(CMO,IFORM,XMAX,MAXIRR,MAXSHL, 3600 & IALIGN) 3601* 3602* A CMOA matrix CMO is given in form defined by CMO. 3603* Analyize the matrix for differences between shells belonginh 3604* to the same shell, and align if required. 3605* 3606* IFORM = 1: Input CMO is in standard order 3607* IFORM = 2: Input CMO is in actual(gas) order 3608* IFORM = 3: Input CMO is in supersymmetry-order 3609* IFORM = 4: Input CMO is in Shell-order 3610* 3611*. Jeppe Olsen, March 2013 3612* 3613 INCLUDE 'implicit.inc' 3614 INCLUDE 'mxpdim.inc' 3615 INCLUDE 'wrkspc-static.inc' 3616 INCLUDE 'glbbas.inc' 3617 INCLUDE 'orbinp.inc' 3618 INCLUDE 'lucinp.inc' 3619#include "errquit.fh" 3620#include "mafdecls.fh" 3621#include "global.fh" 3622*. Input 3623 DIMENSION CMO(*) 3624* 3625 IDUM = 0 3626 CALL MEMMAN(IDUM,IDUM,'MARK ',IDUM,'AN_SHL') 3627* 3628 LCMO = LEN_BLMAT(NSMOB,NTOOBS,NTOOBS,0) 3629 CALL MEMMAN(KLMO1,LCMO,'ADDL ',2,'LMO1 ') 3630* 3631*. Reform to shell format 3632* 3633 CALL REFORM_CMO(CMO,IFORM,dbl_mb(KLMO1),4) 3634* 3635*. And analyze /align 3636* 3637 CALL ANA_SUBSHELLS_CMO_IN(dbl_mb(KLMO1), 3638 & int_mb(KNSUPSYM_FOR_IRREP),int_mb(KIBSUPSYM_FOR_IRREP), 3639 & XMAX,MAXIRR,MAXSHL,IALIGN) 3640* 3641 IF(IALIGN.EQ.1) THEN 3642*. Pump aligned MOs back in CMO 3643 CALL REFORM_CMO(dbl_mb(KLMO1),4,CMO,IFORM) 3644 END IF 3645 3646 CALL MEMMAN(IDUM,IDUM,'FLUSM ',IDUM,'AN_SHL') 3647* 3648 RETURN 3649 END 3650 SUBROUTINE ANA_SUBSHELLS_CMO_IN(CSHELL,NSUPSYM_FOR_IRREP, 3651 & IBSUPSYM_FOR_IRREP,XMAX,MAXIRR,MAXSHL,IALIGN) 3652* 3653* A MO-AO expansion CSHELL is given in SHELL ordered form. 3654* 3655* Find largest deviation between equivalent subshells 3656* align the subshells if IALIGN = 1 3657* 3658*. Jeppe Olsen, March 5, 2013 3659* 3660 INCLUDE 'implicit.inc' 3661 INCLUDE 'mxpdim.inc' 3662 INCLUDE 'orbinp.inc' 3663 INCLUDE 'glbbas.inc' 3664 INCLUDE 'wrkspc-static.inc' 3665* 3666 INTEGER NSUPSYM_FOR_IRREP(*) 3667 INTEGER IBSUPSYM_FOR_IRREP(*) 3668* 3669*. Jeppe Olsen, March, 2013 3670* 3671*. CMO to be analyzed 3672 DIMENSION CSHELL(*) 3673* 3674 NTEST = 100 3675* 3676 IB = 1 3677 XMAX = -1.0D0 3678 DO IRREP = 1, N_SUPSYM_IRREP 3679 NDEG = NSUPSYM_FOR_IRREP(IRREP) 3680 NSHELL = NBAS_SUPSYM(IBSUPSYM_FOR_IRREP(IRREP)) 3681C? WRITE(6,*) ' IRREP, NDEG, NSHELL = ', IRREP, NSHELL, NSHELL 3682 DO ISHELL = 1, NSHELL 3683 IBS = IB + (ISHELL-1)*NSHELL*NDEG 3684*. Check differences to first sub shell 3685 DO ISUB = 2, NDEG 3686 DO IBAS = 1, NSHELL 3687 DIF = CSHELL(IBS-1+(ISUB-1)*NSHELL+IBAS) - 3688 & CSHELL(IBS-1+(1 -1)*NSHELL+IBAS) 3689C? WRITE(6,*) ' IR,ISH,IB, DIF = ', 3690C? & IRREP,ISHELL,IBAS,DIF 3691C? WRITE(6,*) 'CSHELL1, CSHELLI =', 3692C? & CSHELL(IBS-1+(1-1)*NSHELL+IBAS), 3693C? & CSHELL(IBS-1+(ISUB-1)*NSHELL+IBAS) 3694 IF(ABS(DIF).GT.XMAX) THEN 3695 XMAX = ABS(DIF) 3696 MAXIRR = IRREP 3697 MAXSHL= ISHELL 3698 MAXBAS = IBAS 3699 END IF 3700 END DO 3701 END DO 3702 END DO 3703* 3704 IF(IALIGN.EQ.1) THEN 3705*. Align: Copy the first subshell to the remaining 3706 DO ISHELL = 1, NSHELL 3707 IBS = IB + (ISHELL-1)*NSHELL*NDEG 3708 IB1 = IBS 3709 DO ISUB = 2, NDEG 3710 IBI =IBS+(ISUB-1)*NSHELL 3711 CALL COPVEC(CSHELL(IB1),CSHELL(IBI),NSHELL) 3712 END DO 3713 END DO 3714 END IF 3715* 3716 IB = IB + NDEG*NSHELL*NSHELL 3717C? WRITE(6,*) ' TESTY, NDEG, NSHELL, IB = ', 3718C? & NDEG, NSHELL, IB 3719 END DO 3720* 3721 IF(NTEST.GE.100) THEN 3722 WRITE(6,*) ' Info from ANA_SUBSHELLS_ ...' 3723 WRITE(6,*) ' ============================' 3724 WRITE(6,*) 3725 WRITE(6,*) ' Largest difference of subshells: ',XMAX 3726 WRITE(6,*) ' Occurs for, IRREP, shell, basis func ', 3727 & MAXIRR,MAXSHL, MAXBAS 3728 END IF 3729* 3730 RETURN 3731 END 3732 3733 3734c $Id$ 3735