1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18! FILE: rspsoppa.F 19C 20C /* Deck setsoppa */ 21 SUBROUTINE SETSOPPA 22C 23C Set up SOPPA-variables 24C 25C Written 5/4-94 by Erik K. Dalskov 26C 27C Modified by K.Ruud in Dec.-96 in order to avoid static allocation of 28C large particle-hole matrices 29C 30#include "implicit.h" 31C 32C Used from common blocks: 33C INFORB : NSYM,MULD2H(8,8),NRHF(),NVIR(),... 34C INFMP2 : ...,NFRMP2(),IFRMP2(NORBT),... 35C 36#include "maxorb.h" 37#include "maxash.h" 38#include "infpri.h" 39#include "inforb.h" 40#include "infmp2.h" 41#include "infsop.h" 42#include "infind.h" 43#include "infdim.h" 44#include "iratdef.h" 45C 46 CALL QENTER('SETSOPPA') 47C 48C Find offsets in XINDX for AB and IJ address arrays and for 49C the A(2) matrix. They will be placed in infsop.h. 50C 51 KABSAD = 1 52 KABTAD = KABSAD + (N2ORBX + 1)/IRAT 53 KIJSAD = KABTAD + (N2ORBX + 1)/IRAT 54 KIJTAD = KIJSAD + (N2ISHX * NSYM + 1)/IRAT 55 KIJ1AD = KIJTAD + (N2ISHX * NSYM + 1)/IRAT 56 KIJ2AD = KIJ1AD + (N2ISHX * NSYM + 1)/IRAT 57 KIJ3AD = KIJ2AD + (N2ISHX * NSYM + 1)/IRAT 58 KIADR1 = KIJ3AD + (N2ISHX * NSYM + 1)/IRAT 59 KAB2 = KIADR1 + (NH * NP + 1)/IRAT 60C 61C A(2) matrix does not exist in XINDX yet 62C 63 A2EXIST = .FALSE. 64C 65C Redefine some variables in infdim: 66C 67 LCINDX = KAB2 + N2ORBX 68C ...no PHP for SOPPA 69 MAXPHP = 0 70C 71C Find no. of vir/vir-pairs with a>=b (NSVV(SYM)) and a>b (NTVV(SYM)) 72C and no. of occ/occ-pairs with i>=j (NSOO(SYM)) and i>j (NTOO(SYM)) 73C 74 DO I=1,NSYM 75 NSOO(I) = 0 76 NTOO(I) = 0 77 NSVV(I) = 0 78 NTVV(I) = 0 79 ENDDO 80C 81 DO IJSYM=2,NSYM 82 DO ISYM = 1,NSYM 83 JSYM = MULD2H(ISYM,IJSYM) 84 IF (ISYM .LT. JSYM) GOTO 100 85 NHOLEI = NRHF(ISYM) - NFRMP2(ISYM) 86 NHOLEJ = NRHF(JSYM) - NFRMP2(JSYM) 87 NSOO(IJSYM) = NSOO(IJSYM) + NHOLEI*NHOLEJ 88 NTOO(IJSYM) = NTOO(IJSYM) + NHOLEI*NHOLEJ 89 NSVV(IJSYM) = NSVV(IJSYM) + NVIR(ISYM)*NVIR(JSYM) 90 NTVV(IJSYM) = NTVV(IJSYM) + NVIR(ISYM)*NVIR(JSYM) 91 100 CONTINUE 92 ENDDO 93 ENDDO 94C 95 DO ISYM=1,NSYM 96 NHOLEI = NRHF(ISYM) - NFRMP2(ISYM) 97 NSOO(1) = NSOO(1) + NHOLEI*(NHOLEI + 1) / 2 98 NTOO(1) = NTOO(1) + NHOLEI*(NHOLEI - 1) / 2 99 NSVV(1) = NSVV(1) + NVIR(ISYM)*(NVIR(ISYM) + 1) / 2 100 NTVV(1) = NTVV(1) + NVIR(ISYM)*(NVIR(ISYM) - 1) / 2 101 ENDDO 102C 103C Find total no. of 2p-2h variables in each symmetry 104C (singlet and triplet case) 105C 106 DO I=1,NSYM 107 NS2P2H(I) = 0 108 NT2P2H(I) = 0 109 N12P2H(I) = 0 110 N22P2H(I) = 0 111 N32P2H(I) = 0 112 ENDDO 113 114C 115 DO ISYM=1,NSYM 116 DO JSYM=1,NSYM 117 IJSYM = MULD2H(ISYM,JSYM) 118 NS2P2H(IJSYM) = NS2P2H(IJSYM) + NSVV(ISYM) * NSOO(JSYM) 119 NT2P2H(IJSYM) = NT2P2H(IJSYM) + NTVV(ISYM) * NTOO(JSYM) 120 N12P2H(IJSYM) = N12P2H(IJSYM) + NTVV(ISYM) * NTOO(JSYM) 121 N22P2H(IJSYM) = N22P2H(IJSYM) + NSVV(ISYM) * NTOO(JSYM) 122 N32P2H(IJSYM) = N32P2H(IJSYM) + NTVV(ISYM) * NSOO(JSYM) 123 ENDDO 124 ENDDO 125C 126 DO I=1,NSYM 127 N2P2HS(I) = NS2P2H(I) + NT2P2H(I) 128 N2P2HT(I) = N12P2H(I) + N22P2H(I) + N32P2H(I) 129 ENDDO 130 CALL QEXIT('SETSOPPA') 131 RETURN 132 END 133C 134C /* Deck set2soppa */ 135 SUBROUTINE SET2SOPPA(ABSADR,ABTADR,IJSADR,IJTADR, 136 & IJ1ADR,IJ2ADR,IJ3ADR,IADR1) 137#include "implicit.h" 138 INTEGER A,AA,B,BEND,ASYM,BSYM,ABSYM,ABSNDX,ABTNDX 139 PARAMETER (IBIG = -100 000 000) 140#include "maxorb.h" 141#include "maxash.h" 142#include "infpri.h" 143#include "inforb.h" 144#include "infsop.h" 145#include "infmp2.h" 146#include "infind.h" 147C 148 INTEGER ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 149 & IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM), 150 & IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM), 151 & IJ3ADR(NISHT,NISHT,NSYM), IADR1(NP,NH) 152C 153 CALL QENTER('SET2SOPPA') 154C 155C Find offset-arrays for the symmetry packed 2p-2h vectors in SOPPA 156C (singlet properties) and new IJ offsets for triplet properties 157C 158C Moved out of SETSOPPA by K.Ruud, Dec.-96 159C 160 ABSADR(:,:) = 0 161 ABTADR(:,:) = 0 162 DO ABSYM=1,NSYM 163 ABSNDX = 0 164 ABTNDX = 0 165 DO A=1,NP 166 AA = IPHORD(NH+A) 167 ASYM = ISMO(AA) 168 BSYM = MULD2H(ASYM,ABSYM) 169 IF (BSYM .LE. ASYM) THEN 170 BEND = MIN(ISSH(BSYM)+NVIR(BSYM),A) 171 DO B=ISSH(BSYM)+1,BEND 172 IF (A .EQ. B) THEN 173 ABSNDX = ABSNDX + 1 174 ABSADR(A,A) = ABSNDX 175 ABTADR(A,A) = IBIG 176 ELSE 177 ABSNDX = ABSNDX + 1 178 ABTNDX = ABTNDX + 1 179 ABSADR(A,B) = ABSNDX 180 ABTADR(A,B) = ABTNDX 181 ABSADR(B,A) = ABSNDX 182 ABTADR(B,A) = ABTNDX 183 ENDIF 184 ENDDO 185 ENDIF 186 ENDDO 187 ENDDO 188C 189C 190 DO LSYMOP=1,NSYM 191 IJSNDX = 0 192 IJTNDX = NS2P2H(LSYMOP) 193 IJ1NDX = 0 194 IJ2NDX = N12P2H(LSYMOP) 195 IJ3NDX = N22P2H(LSYMOP) + N12P2H(LSYMOP) 196 DO I=1,NH 197 II = IPHORD(I) 198 ISYM = ISMO(II) 199 DO J=1,I 200 JJ = IPHORD(J) 201 JSYM = ISMO(JJ) 202 IJSYM = MULD2H(ISYM,JSYM) 203 ABSYM = MULD2H(IJSYM,LSYMOP) 204 IF (I .EQ. J) THEN 205 IF (IFRMP2(II) .EQ. 0) THEN 206 IJSADR(I,I,LSYMOP) = IJSNDX 207 IJTADR(I,I,LSYMOP) = IBIG 208 IJ1ADR(I,I,LSYMOP) = IBIG 209 IJ2ADR(I,I,LSYMOP) = IBIG 210 IJ3ADR(I,I,LSYMOP) = IJ3NDX 211 IJSNDX = IJSNDX + NSVV(ABSYM) 212 IJ3NDX = IJ3NDX + NTVV(ABSYM) 213 ELSE 214 IJSADR(I,I,LSYMOP) = IBIG 215 IJTADR(I,I,LSYMOP) = IBIG 216 IJ1ADR(I,I,LSYMOP) = IBIG 217 IJ2ADR(I,I,LSYMOP) = IBIG 218 IJ3ADR(I,I,LSYMOP) = IBIG 219 END IF 220 ELSE 221 IF (IFRMP2(II) .EQ. 0 .AND. IFRMP2(JJ) .EQ. 0) THEN 222 IJSADR(I,J,LSYMOP) = IJSNDX 223 IJSADR(J,I,LSYMOP) = IJSNDX 224 IJTADR(I,J,LSYMOP) = IJTNDX 225 IJTADR(J,I,LSYMOP) = IJTNDX 226 IJ1ADR(I,J,LSYMOP) = IJ1NDX 227 IJ1ADR(J,I,LSYMOP) = IJ1NDX 228 IJ2ADR(I,J,LSYMOP) = IJ2NDX 229 IJ2ADR(J,I,LSYMOP) = IJ2NDX 230 IJ3ADR(I,J,LSYMOP) = IJ3NDX 231 IJ3ADR(J,I,LSYMOP) = IJ3NDX 232 IJSNDX = IJSNDX + NSVV(ABSYM) 233 IJTNDX = IJTNDX + NTVV(ABSYM) 234 IJ1NDX = IJ1NDX + NTVV(ABSYM) 235 IJ2NDX = IJ2NDX + NSVV(ABSYM) 236 IJ3NDX = IJ3NDX + NTVV(ABSYM) 237 ELSE 238 IJSADR(I,J,LSYMOP) = IBIG 239 IJSADR(J,I,LSYMOP) = IBIG 240 IJTADR(I,J,LSYMOP) = IBIG 241 IJTADR(J,I,LSYMOP) = IBIG 242 IJ1ADR(I,J,LSYMOP) = IBIG 243 IJ1ADR(J,I,LSYMOP) = IBIG 244 IJ2ADR(I,J,LSYMOP) = IBIG 245 IJ2ADR(J,I,LSYMOP) = IBIG 246 IJ3ADR(I,J,LSYMOP) = IBIG 247 IJ3ADR(J,I,LSYMOP) = IBIG 248 ENDIF 249 ENDIF 250 ENDDO 251 ENDDO 252 ENDDO 253C 254C Find also IADR1 255C 256 CALL PHADR1(IADR1,1) 257C 258C End of SET2SOPPA 259C 260 CALL QEXIT('SET2SOPPA') 261 RETURN 262 END 263C /* Deck a0s2 */ 264 SUBROUTINE A0S2(FCONE,FC,DONEPT,ZYMAT,TZYMAT,STWO, 265 * ORBE,EPSIL,NSIM,WRK,LWRK) 266C 267C Copyright by Martin Packer 251193. In formulas below 268C i,j,k.. are occupied and a,b,c... are virtual orbitals. 269C We must form FCONE = 0.5*(A(0)S(2)ZYMAT + S(2)A(0)ZYMAT) 270C IN ORDER TO ACCOUNT FOR THE NON-HERMITICITY. 271C 272C 940726-hjaaj: corrected code for NSIM .gt. 1 273C 274#include "implicit.h" 275C 276C 277#include "maxorb.h" 278#include "maxash.h" 279#include "inforb.h" 280#include "infrsp.h" 281#include "wrkrsp.h" 282C 283 DIMENSION FCONE(N2ORBX,*), ORBE(NORBT,NORBT), FC(*) 284 DIMENSION WRK(*), ZYMAT(NORBT,NORBT,*), TZYMAT(NORBT,NORBT,*) 285 DIMENSION DONEPT(NORBT,NORBT), STWO(NORBT,NORBT), EPSIL(NORBT) 286C 287 PARAMETER( DMP25 = -0.25D0 ) 288C 289C 290C First PUT THE ORBITAL ENERGIES FROM FC INTO THE ARRAY ORBE 291C 292 IS = 0 293 DO IM = 1,NSYM 294 IORBM = IORB(IM) 295 NORBM = NORB(IM) 296 IIORBM = IIORB(IM) 297 IF(NORBM.NE.0) THEN 298 DO IJ = 1,NORBM 299 IS = IS + 1 300 IPL = IIORBM + (IJ*(IJ+1))/2 301 EPSIL(IS) = FC(IPL) 302 ENDDO 303 ENDIF 304 ENDDO 305C 306C FORM S(2).ZYMAT MATRIX AND PLACE IN STWO 307C 308 DO I = 1,NSIM 309 CALL DZERO(STWO,N2ORBX) 310 DO IASYM = 1,NSYM 311 IBSYM = MULD2H(IASYM,KSYMOP) 312 NOCCA = NOCC(IASYM) 313 NSSHB = NSSH(IBSYM) 314 IF ((NOCCA.NE.0).AND.(NSSHB.NE.0)) THEN 315C S(AI) TERMS 316 IAST = IORB(IASYM) + 1 317 IBST = IORB(IBSYM) + 1 + NOCC(IBSYM) 318 CALL DGEMM('N','N',NSSHB,NOCCA,NOCCA,1.D0, 319 & ZYMAT(IBST,IAST,I),NORBT, 320 & DONEPT(IAST,IAST),NORBT,1.D0, 321 & STWO(IBST,IAST),NORBT) 322 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHB,-1.D0, 323 & DONEPT(IBST,IBST),NORBT, 324 & ZYMAT(IBST,IAST,I),NORBT,1.D0, 325 & STWO(IBST,IAST),NORBT) 326C S(IA) TERMS 327 CALL DGEMM('T','N',NOCCA,NSSHB,NSSHB,1.D0, 328 & TZYMAT(IBST,IAST,I),NORBT, 329 & DONEPT(IBST,IBST),NORBT,1.D0, 330 & STWO(IAST,IBST),NORBT) 331 CALL DGEMM('N','T',NOCCA,NSSHB,NOCCA,-1.D0, 332 & DONEPT(IAST,IAST),NORBT, 333 & TZYMAT(IBST,IAST,I),NORBT,1.D0, 334 & STWO(IAST,IBST),NORBT) 335 ENDIF 336 ENDDO 337C 338C NOW FORM A(0)S(2).ZYMAT = G, AND STORE IN STWO: 339C G = SUM(JB) (EPSILON(A)-EPSILON(I))*S(AI,JB).ZYMAT(BJ) 340C and FORM A(0).ZYMAT = H, AND STORE IN ORBE: 341C H(J,B) = (EPSILON(B)-EPSILON(J))ZYMAT(BJ) 342C 343C CALL DZERO(ORBE,N2ORBX) 344 DO ICSYM = 1,NSYM 345 IDSYM = MULD2H(ICSYM,KSYMOP) 346 IORBC = IORB(ICSYM) 347 IORBD = IORB(IDSYM) 348 NORBD = NORB(IDSYM) 349 NOCCC = NOCC(ICSYM) 350 NOCCD = NOCC(IDSYM) 351 DO IK = (IORBC + 1),(IORBC + NOCCC) 352 DO IL = (IORBD + 1 + NOCCD),(IORBD + NORBD) 353 STWO(IK,IL) = 354 * STWO(IK,IL)*(EPSIL(IL)-EPSIL(IK)) 355 STWO(IL,IK) = 356 * STWO(IL,IK)*(EPSIL(IL)-EPSIL(IK)) 357 ORBE(IK,IL) = 358 * ZYMAT(IK,IL,I)*(EPSIL(IL)-EPSIL(IK)) 359 ORBE(IL,IK) = 360 * ZYMAT(IL,IK,I)*(EPSIL(IL)-EPSIL(IK)) 361 ENDDO 362 ENDDO 363 ENDDO 364C 365C FORM S(2)*ORBE, WHICH IS THE CONJUGATE TERM TO H 366 DO IASYM = 1,NSYM 367 IBSYM = MULD2H(IASYM,KSYMOP) 368 NOCCA = NOCC(IASYM) 369 NSSHB = NSSH(IBSYM) 370 IF ((NOCCA.NE.0).AND.(NSSHB.NE.0)) THEN 371C S(AI) TERMS 372 IAST = IORB(IASYM) + 1 373 IBST = IORB(IBSYM) + 1 + NOCC(IBSYM) 374 CALL DGEMM('N','N',NSSHB,NOCCA,NOCCA,1.D0, 375 & ORBE(IBST,IAST),NORBT, 376 & DONEPT(IAST,IAST),NORBT,1.D0, 377 & STWO(IBST,IAST),NORBT) 378 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHB,-1.D0, 379 & DONEPT(IBST,IBST),NORBT, 380 & ORBE(IBST,IAST),NORBT,1.D0, 381 & STWO(IBST,IAST),NORBT) 382C S(IA) TERMS 383 CALL DGEMM('N','N',NOCCA,NSSHB,NSSHB,1.D0, 384 & ORBE(IAST,IBST),NORBT, 385 & DONEPT(IBST,IBST),NORBT,1.D0, 386 & STWO(IAST,IBST),NORBT) 387 CALL DGEMM('N','N',NOCCA,NSSHB,NOCCA,-1.D0, 388 & DONEPT(IAST,IAST),NORBT, 389 & ORBE(IAST,IBST),NORBT,1.D0, 390 & STWO(IAST,IBST),NORBT) 391 ENDIF 392 ENDDO 393C PLACE STWO INTO FCONE, WITH THE CORRECT SCALE FACTOR 394 CALL DAXPY(N2ORBX,DMP25, STWO,1, FCONE(1,I),1) 395 ENDDO 396 RETURN 397 END 398C /* Deck sopudv */ 399 SUBROUTINE SOPUDV(DONEPT) 400C 401C Copyright by Martin Packer 251193. In formulas below 402C i,j,k.. are occupied and a,b,c... are virtual orbitals. 403C 404C 405#include "implicit.h" 406C 407C 408#include "maxorb.h" 409#include "maxash.h" 410#include "inforb.h" 411#include "infind.h" 412#include "infrsp.h" 413C 414 DIMENSION DONEPT(NORBT,NORBT) 415C 416 PARAMETER(D2 = 2.0D0) 417C 418C REMOVE THE DIAGONAL S(0) CONTRIBUTION FROM DONEPT (ADDED IN MP2FAC) 419 DO ISY = 1,NSYM 420 IORBSY = IORB(ISY) 421 NOCCSY = NOCC(ISY) 422 NORBSY = NORB(ISY) 423 IF (NORBSY.NE.0) THEN 424 DO IA = (IORBSY+1),(IORBSY+NOCCSY) 425 DONEPT(IA,IA) = DONEPT(IA,IA) - D2 426 ENDDO 427 ENDIF 428 ENDDO 429 RETURN 430 END 431C /* Deck hrpa */ 432 SUBROUTINE HRPA (FCONE,DONEPT,COEMP2,H2,ZYMAT,TZYMAT,H2X,H2XP, 433 * COEUNP,ICI,IDI,ICDSYM,ITYPCD,NSIM,IADR1,WRK,LWRK) 434C 435C Copyright by Martin Packer 191193. In formulas below 436C i,j,k.. are occupied and a,b,c... are virtual orbitals. 437C 438C 439#include "implicit.h" 440C 441C 442#include "maxorb.h" 443#include "maxash.h" 444#include "infmp2.h" 445 DIMENSION IADR1(NP,NH) 446#include "inforb.h" 447#include "orbtypdef.h" 448#include "infrsp.h" 449#include "wrkrsp.h" 450C 451 DIMENSION H2(NORBT,NORBT), FCONE(N2ORBX,*), COEMP2(*) 452 DIMENSION DONEPT(NORBT,NORBT) 453 DIMENSION WRK(*), ZYMAT(NORBT,NORBT,*), TZYMAT(NORBT,NORBT,*) 454 DIMENSION H2X(NORBT,NORBT), H2XP(NORBT,NORBT), COEUNP(*) 455C 456 CALL QENTER('HRPA ') 457C 458C 459C use distribution type ITYPCD = 460C 1:inactive-inactive 2:active-inactive 3:active-active 461C 4:secondary-inactive 5:secondary-active 6:secondary-secondary 462C We only need types 1, 4 and 6. 463C 464 IF (ITYPCD .EQ. 1) THEN 465 IF (IFRMP2(ICI) .NE. 0 .AND. IFRMP2(IDI) .NE. 0) GO TO 9999 466 ELSE IF (ITYPCD .EQ. 4) THEN 467 IF (IFRMP2(ICI) .NE. 0) GO TO 9999 468 ELSE IF (ITYPCD .NE. 6) THEN 469 GO TO 9999 470 END IF 471C 472C One index transform the two-e integrals 473 DO I = 1,NSIM 474 CALL DZERO(H2X ,N2ORBX) 475 IF (ITYPCD .NE. 4) CALL DZERO(H2XP,N2ORBX) 476 DO IPSYM = 1,NSYM 477 IBSYM = MULD2H(IPSYM,ICDSYM) 478 IASYM = MULD2H(IPSYM,KSYMOP) 479 NOCCP = NOCC(IPSYM) 480 NOCCB = NOCC(IBSYM) 481 NOCCA = NOCC(IASYM) 482 NSSHP = NSSH(IPSYM) 483 NSSHB = NSSH(IBSYM) 484 NSSHA = NSSH(IASYM) 485cmjp branch to code for each distribution type 486 IF (ITYPCD .EQ. 6) GOTO 600 487 IF (ITYPCD .EQ. 4) GOTO 400 488 IF (ITYPCD .NE. 1) GOTO 9999 489C 490cmjp occupied-occupied distributions 491C 492C This is the transform of (CD|**) with C,D occupied 493C H2X(K,C) = SUM(M) H2(K,M)ZYMAT(M,C) - SUM(E) ZYMAT(K,E)H2(E,C): F(AI) 494C H2X(C,L) = SUM(E)H2(C,E)ZYMAT(E,L) -SUM(M)ZYMAY(C,M)H2(M,L): F(IA) 495C 496C TRANSFORM OCC-OCC PART OF H2 497C 498 IF (NOCCB .GT. 0 .AND. NOCCP .GT. 0) THEN 499 IAST = IORB(IASYM) + 1 + NOCCA 500 IBST = IORB(IBSYM) + 1 501 IPST = IORB(IPSYM) + 1 502C Now transform for those parts which enter F(ai) 503 CALL DGEMM('N','N',NSSHA,NOCCB,NOCCP,1.D0, 504 & TZYMAT(IAST,IPST,I),NORBT, 505 & H2(IPST,IBST),NORBT,1.D0, 506 & H2XP(IAST,IBST),NORBT) 507C CALL AMPAB(H2 (IBST,IPST),NOCCB,NOCCP,NORBT,NORBT, 508C * ZYMAT (IPST,IAST,I),NOCCP,NSSHA,NORBT,NORBT, 509C * H2X (IBST,IAST),NORBT,NORBT ) 510C Now transform for those parts which enter F(ia) 511 CALL DGEMM('N','N',NSSHA,NOCCB,NOCCP,-1.D0, 512 & ZYMAT(IAST,IPST,I),NORBT, 513 & H2(IPST,IBST),NORBT,1.D0, 514 & H2X(IAST,IBST),NORBT) 515 END IF 516C TRANSFORM VIRT-VIRT PART OF H2 517 IF (NOCCA .GT. 0) THEN 518 IAST = IORB(IASYM) + 1 519 IBST = IORB(IBSYM) + 1 + NOCCB 520 IPST = IORB(IPSYM) + 1 + NOCCP 521C Now transform for those parts which enter F(ia) 522 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHP,1.D0, 523 & H2(IBST,IPST),NORBT, 524 & ZYMAT(IPST,IAST,I),NORBT,1.D0, 525 & H2X(IBST,IAST),NORBT) 526C Now transform for those parts which enter F(ai) 527 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHP,-1.D0, 528 & H2(IBST,IPST),NORBT, 529 & TZYMAT(IPST,IAST,I),NORBT,1.D0, 530 & H2XP(IBST,IAST),NORBT) 531C CALL SMPAB( ZYMAT(IAST,IPST,I),NOCCA,NSSHP,NORBT,NORBT, 532C * H2 (IPST,IBST),NSSHP,NSSHB,NORBT,NORBT, 533C * H2X (IAST,IBST),NORBT,NORBT ) 534 END IF 535C 536 GOTO 500 537C 538cmjp occupied-virtual distributions 539C Perform two transforms here, which enter different parts 540C of h2x, where h2 is an occ-virt distribution. 541C H2X(A,D)=SUM(M){ZYMAT(A,M)H2(M,D)-H2(A,M)ZYMAT(M,D)} 542C H2X(L,I)=SUM(E){H2(L,E)ZYMAT(E,I)-ZYMAT(L,E)H2(E,I)} 543C 544C 545C Transform occ-virt part of H2 546C 547 400 CONTINUE 548 IF (NOCCP .GT. 0) THEN 549 IAST = IORB(IASYM) + 1 + NOCCA 550 IBST = IORB(IBSYM) + 1 + NOCCB 551 IPST = IORB(IPSYM) + 1 552 CALL DGEMM('N','N',NSSHA,NSSHB,NOCCP,1.D0, 553 & ZYMAT(IAST,IPST,I),NORBT, 554 & H2(IPST,IBST),NORBT,1.D0, 555 & H2X(IAST,IBST),NORBT) 556C Transform virt-occ part of H2 557 CALL DGEMM('N','N',NSSHB,NSSHA,NOCCP,-1.D0, 558 & H2(IBST,IPST),NORBT, 559 & ZYMAT(IPST,IAST,I),NORBT,1.D0, 560 & H2X(IBST,IAST),NORBT) 561 END IF 562C 563C Now do transform which goes into H2X(occ,occ) 564C Transform occ-virt part of H2 565C 566 IF (NOCCB .GT. 0 .AND. NOCCA .GT. 0) THEN 567 IAST = IORB(IASYM) + 1 568 IBST = IORB(IBSYM) + 1 569 IPST = IORB(IPSYM) + 1 + NOCCP 570 CALL DGEMM('T','N',NOCCB,NOCCA,NSSHP,1.D0, 571 & H2(IPST,IBST),NORBT, 572 & ZYMAT(IPST,IAST,I),NORBT,1.D0, 573 & H2X(IBST,IAST),NORBT) 574C Transform virt-occ part of H2 575 CALL DGEMM('T','N',NOCCA,NOCCB,NSSHP,-1.D0, 576 & TZYMAT(IPST,IAST,I),NORBT, 577 & H2(IPST,IBST),NORBT,1.D0, 578 & H2X(IAST,IBST),NORBT) 579 END IF 580 GOTO 500 581C 582cmjp virtual-virtual distributions 583C transform H2X(L,C)=SUM(E)ZYMAT(L,E)H2(E,C) - SUM(M)H2(L,M)ZYMAT(M,C) 584C where h2 is a virt-virt distribution. This contributes to F(ai). 585C 586C transform H2X(D,L)=SUM(M)ZYMAT(D,M)H2(M,L) - SUM(E)H2(D,E)ZYMAT(E,L) 587C where h2 is a virt-virt distribution. This contributes to F(ia). 588C 589C 590C TRANSFORM VIRT-VIRT PART OF H2 591C 592 600 CONTINUE 593 IF (NOCCA .GT. 0) THEN 594 IAST = IORB(IASYM) + 1 595 IBST = IORB(IBSYM) + 1 + NOCCB 596 IPST = IORB(IPSYM) + 1 + NOCCP 597C NOW DO F(AI) CONTRIBUTION 598 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHP,1.D0, 599 & H2(IBST,IPST),NORBT, 600 & TZYMAT(IPST,IAST,I),NORBT,1.D0, 601 & H2XP(IBST,IAST),NORBT) 602C CALL AMPAB( ZYMAT (IAST,IPST,I),NOCCA,NSSHP,NORBT,NORBT, 603C * H2 (IPST,IBST),NSSHP,NSSHB,NORBT,NORBT, 604C * H2X (IAST,IBST),NORBT,NORBT ) 605C NOW DO F(IA) CONTRIBUTION 606 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHP,-1.D0, 607 & H2(IBST,IPST),NORBT, 608 & ZYMAT(IPST,IAST,I),NORBT,1.D0, 609 & H2X(IBST,IAST),NORBT) 610 END IF 611C TRANSFORM OCC-OCC PART OF H2 612 IF (NOCCP .GT. 0 .AND. NOCCB .GT. 0) THEN 613 IAST = IORB(IASYM) + 1 + NOCCA 614 IBST = IORB(IBSYM) + 1 615 IPST = IORB(IPSYM) + 1 616C NOW DO F(IA) CONTRIBUTION 617 CALL DGEMM('N','N',NSSHA,NOCCB,NOCCP,1.D0, 618 & ZYMAT(IAST,IPST,I),NORBT, 619 & H2(IPST,IBST),NORBT,1.D0, 620 & H2X(IAST,IBST),NORBT) 621C NOW DO F(AI) CONTRIBUTION 622 CALL DGEMM('N','N',NSSHA,NOCCB,NOCCP,-1.D0, 623 & TZYMAT(IAST,IPST,I),NORBT, 624 & H2(IPST,IBST),NORBT,1.D0, 625 & H2XP(IAST,IBST),NORBT) 626C CALL SMPAB( H2 (IBST,IPST),NOCCB,NOCCP,NORBT,NORBT, 627C * ZYMAT (IPST,IAST,I),NOCCP,NSSHA,NORBT,NORBT, 628C * H2X (IBST,IAST),NORBT,NORBT ) 629C 630 END IF 631 500 CONTINUE 632 ENDDO 633C Now call HRPAF to add contributions to the FCONE matrix 634C HRPAF assumes C .ge. D, we know C .le. D therefore 635C we switch ICI and IDI in call of HRPAF 636C 637 CALL HRPAF(FCONE(1,I),DONEPT,COEMP2,H2X,H2XP,COEUNP,IDI,ICI, 638 * ICDSYM,ITYPCD,IADR1,WRK,LWRK) 639 ENDDO 640C 641C End of HRPA 642C 643 9999 CALL QEXIT('HRPA ') 644 RETURN 645 END 646C /* Deck hrpaf */ 647 SUBROUTINE HRPAF (FCONE,DONEPT,COEMP2,H2X,H2XP,COEUNP, 648 * ICI,IDI,ICDSYM,ISEL,IADR1,WRK,LWRK) 649C 650C Copyright by Martin Packer 191193. In the formulas below 651C i,j,k.. are occupied and a,b,c... are virtual orbitals. 652C Construct the HRPA parts of the transformed Fock matrix. 653C ISEL=1 OCC-OCC DISTRIBUTION 654C ISEL=4 OCC-VIRT DISTRIBUTION 655C ISEL=6 VIRT-VIRT DISTRIBUTION 656C 657#include "implicit.h" 658C 659#include "maxorb.h" 660#include "maxash.h" 661#include "infmp2.h" 662 DIMENSION IADR1(NP,NH) 663#include "inforb.h" 664#include "infind.h" 665#include "infrsp.h" 666#include "wrkrsp.h" 667C 668 DIMENSION FCONE(NORBT,NORBT), WRK(*), H2XP(NPHMAX) 669 DIMENSION COEMP2(*), H2X(NORBT,NORBT) 670 DIMENSION COEUNP(NORBT,NORBT), DONEPT(NORBT,NORBT) 671C 672 IMSYM = MULD2H(ICDSYM,KSYMOP) 673 IF (TRPLET) THEN 674 MPOFF = LPVMAT 675 ELSE 676 MPOFF = 0 677 END IF 678 IF (ISEL .EQ. 1) THEN 679C Occ-occ distribution 680C Use INDXPH to indicate which occ or virt orbital ICI and IDI 681C correspond to. 682 ICP = INDXPH(ICI) 683 IDP = INDXPH(IDI) 684C F(ai) contribution from H2X(occ,virt) stored in H2XP(virt,occ) 685C which we pack into COEUNP 686 CALL PHPACK(H2XP,COEUNP,IMSYM,0,-1) 687C F(ia) contribution from H2X(virt,occ) packed into H2XP 688 CALL PHPACK(H2X ,H2XP, IMSYM,0,-1) 689 DO IASYM = 1,NSYM 690 IBSYM = MULD2H(IASYM,ISMO(ICI)) 691 ICSYM = MULD2H(IASYM,ISMO(IDI)) 692 NOCCA = NOCC(IASYM) 693 NSSHA = NSSH(IASYM) 694 IAST = IORB(IASYM) + 1 + NOCCA 695 IASTM = -INDXPH(IAST) 696 IF((KSYMOP.NE.IBSYM).OR.(ICSYM.NE.IMSYM)) GOTO 20 697 IF (IFRMP2(IDI).EQ.0) THEN 698 DO J = 0,NSSHA-1 699 FCONE(IAST+J,ICI) = FCONE(IAST+J,ICI) + 700 * DDOT(NPHSYM(IMSYM),COEUNP,1, 701 * COEMP2(MPOFF+IADR1(IASTM+J,IDP)+1),1) 702 FCONE(ICI,IAST+J) = FCONE(ICI,IAST+J) + 703 * DDOT(NPHSYM(IMSYM),H2XP ,1, 704 * COEMP2(MPOFF+IADR1(IASTM+J,IDP)+1),1) 705 ENDDO 706 ENDIF 707 20 IF((KSYMOP.NE.ICSYM).OR.(IBSYM.NE.IMSYM)) GOTO 30 708 IF (IFRMP2(ICI).EQ.0) THEN 709 IF(ICI .NE. IDI) THEN 710 DO J = 0,NSSHA-1 711 FCONE(IAST+J,IDI) = FCONE(IAST+J,IDI) + 712 * DDOT(NPHSYM(IMSYM),COEUNP,1, 713 * COEMP2(MPOFF+IADR1(IASTM+J,ICP)+1),1) 714 FCONE(IDI,IAST+J) = FCONE(IDI,IAST+J) + 715 * DDOT(NPHSYM(IMSYM),H2XP ,1, 716 * COEMP2(MPOFF+IADR1(IASTM+J,ICP)+1),1) 717 ENDDO 718 ENDIF 719 ENDIF 720 30 CONTINUE 721 ENDDO 722 ENDIF 723C NOW FOR VIRT-VIRT DISTRIBUTION 724 IF (ISEL .EQ. 6) THEN 725 ICP = -INDXPH(ICI) 726 IDP = -INDXPH(IDI) 727C F(ai) contribution from H2X(occ,virt) stored in H2XP(virt,occ) 728 CALL PHPACK(H2XP,COEUNP,IMSYM,0,-1) 729C F(ia) contribution from H2X(virt,occ) 730 CALL PHPACK(H2X ,H2XP ,IMSYM,0,-1) 731 DO IASYM = 1,NSYM 732 IBSYM = MULD2H(IASYM,ISMO(ICI)) 733 ICSYM = MULD2H(IASYM,ISMO(IDI)) 734 NOCCA = NOCC(IASYM) 735 NSSHA = NSSH(IASYM) 736 IAST = IORB(IASYM) + 1 737 IASTM = INDXPH(IAST) 738 IF((KSYMOP.NE.IBSYM).OR.(ICSYM.NE.IMSYM)) GOTO 60 739 DO J = NFRMP2(IASYM),NOCCA-1 740 FCONE(ICI,IAST+J) = FCONE(ICI,IAST+J) + 741 * DDOT(NPHSYM(IMSYM),COEUNP,1, 742 * COEMP2(MPOFF+IADR1(IDP,IASTM+J)+1),1) 743 FCONE(IAST+J,ICI) = FCONE(IAST+J,ICI) + 744 * DDOT(NPHSYM(IMSYM),H2XP ,1, 745 * COEMP2(MPOFF+IADR1(IDP,IASTM+J)+1),1) 746 ENDDO 747 60 IF((KSYMOP.NE.ICSYM).OR.(IBSYM.NE.IMSYM)) GOTO 70 748 IF (ICI .NE. IDI) THEN 749 DO J = NFRMP2(IASYM),NOCCA-1 750 FCONE(IDI,IAST+J) = FCONE(IDI,IAST+J) + 751 * DDOT(NPHSYM(IMSYM),COEUNP,1, 752 * COEMP2(MPOFF+IADR1(ICP,IASTM+J)+1),1) 753 FCONE(IAST+J,IDI) = FCONE(IAST+J,IDI) + 754 * DDOT(NPHSYM(IMSYM),H2XP ,1, 755 * COEMP2(MPOFF+IADR1(ICP,IASTM+J)+1),1) 756 ENDDO 757 ENDIF 758 70 CONTINUE 759 ENDDO 760 ENDIF 761C NOW DO OCCUPIED-VIRTUAL DISTRIBUTIONS 762 IF (ISEL.EQ.4) THEN 763C skip if IDI is frozen in MP2: 764 IF (IFRMP2(IDI).NE.0) GO TO 9999 765 ICP = -INDXPH(ICI) 766 IDP = INDXPH(IDI) 767C unpack ph block in COEUNP (i.e. COEUNP(i,A) not defined) 768 CALL PHPACK(COEUNP,COEMP2(IADR1(ICP,IDP)+1),ICDSYM,-1,-1) 769 DO IASYM = 1, NSYM 770 IBSYM = MULD2H(IASYM,ICDSYM) 771 ICSYM = MULD2H(IASYM,IMSYM) 772C CONSTRUCT F(AI) TERMS FROM VIRT-VIRT PART OF H2X 773 NOCCB = NOCC(IBSYM) 774 IF (NOCCB .GT. 0) THEN 775 NSSHA = NSSH(IASYM) 776 NSSHC = NSSH(ICSYM) 777 IAST = IORB(IASYM) + 1 + NOCC(IASYM) 778 IBST = IORB(IBSYM) + 1 779 ICST = IORB(ICSYM) + 1 + NOCC(ICSYM) 780 CALL DGEMM('N','N',NSSHC,NOCCB,NSSHA,1.D0, 781 & H2X(ICST,IAST),NORBT, 782 & COEUNP(IAST,IBST),NORBT,1.D0, 783 & FCONE(ICST,IBST),NORBT) 784C CONSTRUCT F(IA) TERMS FROM VIRT-VIRT PART OF H2X 785 CALL DGEMM('T','N',NOCCB,NSSHC,NSSHA,1.D0, 786 & COEUNP(IAST,IBST),NORBT, 787 & H2X(IAST,ICST),NORBT,1.D0, 788 & FCONE(IBST,ICST),NORBT) 789 ENDIF 790C 791C CONSTRUCT F(AI) TERMS FROM OCC-OCC PART OF H2X 792C 793 NOCCA = NOCC(IASYM) 794 NOCCC = NOCC(ICSYM) 795 IF (NOCCA .GT. 0 .AND. NOCCC .GT. 0) THEN 796 NSSHB = NSSH(IBSYM) 797 IAST = IORB(IASYM) + 1 798 IBST = IORB(IBSYM) + 1 + NOCCB 799 ICST = IORB(ICSYM) + 1 800 CALL DGEMM('N','N',NSSHB,NOCCC,NOCCA,1.D0, 801 & COEUNP(IBST,IAST),NORBT, 802 & H2X(IAST,ICST),NORBT,1.D0, 803 & FCONE(IBST,ICST),NORBT) 804C CONSTRUCT F(IA) TERMS FROM OCC-OCC PART OF H2X 805 CALL DGEMM('N','T',NOCCC,NSSHB,NOCCA,1.D0, 806 & H2X(ICST,IAST),NORBT, 807 & COEUNP(IBST,IAST),NORBT,1.D0, 808 & FCONE(ICST,IBST),NORBT) 809 ENDIF 810 ENDDO 811 ENDIF 812 9999 RETURN 813C 814C End of HRPAF 815C 816 END 817C /* Deck rsps2m */ 818 SUBROUTINE RSPS2M(NSIM,DONEPT,SVEC,ZYMAT,STWO) 819C 820C Copyright by Martin Packer 251193. In formulas below 821C i,j,k.. are occupied and a,b,c... are virtual orbitals. 822C 823C 824#include "implicit.h" 825C 826C 827#include "maxorb.h" 828#include "maxash.h" 829#include "infmp2.h" 830#include "inforb.h" 831#include "infrsp.h" 832#include "infvar.h" 833#include "wrkrsp.h" 834C 835 DIMENSION ZYMAT(NORBT,NORBT,*), DONEPT(NORBT,*) 836 DIMENSION STWO(NORBT,NORBT,*), SVEC(KZYVAR,*) 837C 838 PARAMETER(D2 = 2.0D0) 839C 840 KYCONF = KZCONF + KZVAR 841C 842 CALL DZERO(STWO,NORBT*NORBT*NSIM) 843C FORM S(2).ZYMAT MATRIX AND PLACE IN STWO 844 DO I = 1,NSIM 845 DO IASYM = 1,NSYM 846 IBSYM = MULD2H(IASYM,KSYMOP) 847 NOCCA = NOCC(IASYM) 848 NSSHB = NSSH(IBSYM) 849 IF ((NOCCA.NE.0).AND.(NSSHB.NE.0)) THEN 850C S(AI) TERMS 851 IAST = IORB(IASYM) + 1 852 IBST = IORB(IBSYM) + 1 + NOCC(IBSYM) 853 CALL DGEMM('N','N',NSSHB,NOCCA,NOCCA,1.D0, 854 & ZYMAT(IBST,IAST,I),NORBT, 855 & DONEPT(IAST,IAST),NORBT,1.D0, 856 & STWO(IBST,IAST,I),NORBT) 857 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHB,-1.D0, 858 & DONEPT(IBST,IBST),NORBT, 859 & ZYMAT(IBST,IAST,I),NORBT,1.D0, 860 & STWO(IBST,IAST,I),NORBT) 861C S(IA) TERMS 862 CALL DGEMM('N','N',NOCCA,NSSHB,NSSHB,1.D0, 863 & ZYMAT(IAST,IBST,I),NORBT, 864 & DONEPT(IBST,IBST),NORBT,1.D0, 865 & STWO(IAST,IBST,I),NORBT) 866 CALL DGEMM('N','N',NOCCA,NSSHB,NOCCA,-1.D0, 867 & DONEPT(IAST,IAST),NORBT, 868 & ZYMAT(IAST,IBST,I),NORBT,1.D0, 869 & STWO(IAST,IBST,I),NORBT) 870 ENDIF 871 ENDDO 872 ENDDO 873C NOW PACK STWO INTO SVEC USING CODE ADAPTED FROM RSPSOR 874 DO IG = 1,KZWOPT 875 K = JWOP(1,IG) 876 L = JWOP(2,IG) 877 DO ISIM = 1 ,NSIM 878 SVEC(KZCONF+IG,ISIM) = SVEC(KZCONF+IG,ISIM) 879 * + STWO(L,K,ISIM) 880 SVEC(KYCONF+IG,ISIM) = SVEC(KYCONF+IG,ISIM) 881 * + STWO(K,L,ISIM) 882 SVEC(KZCONF+IG,ISIM) = SVEC(KZCONF+IG,ISIM) 883 * + D2 * ZYMAT(L,K,ISIM) 884 SVEC(KYCONF+IG,ISIM) = SVEC(KYCONF+IG,ISIM) 885 * - D2 * ZYMAT(K,L,ISIM) 886 ENDDO 887 ENDDO 888 RETURN 889C 890C End of RSPS2M 891C 892 END 893C /* Deck hrpaa2 */ 894 SUBROUTINE HRPAA2(COEMP2,H2,A2TEMP,COEUNP,IADR1,ICI,IDI,ICDSYM) 895C 896C Copyright by Martin Packer 110194. In the formulas below 897C i,j,k.. are occupied and a,b,c... are virtual orbitals. 898C Add the contribution which symmetrises the A matrix. 899C 900C Changed to construction of A(2) only. 940916-ekd 901C A new routine is then added to symmetrise A. 902C 903#include "implicit.h" 904C 905#include "maxorb.h" 906#include "maxash.h" 907#include "infmp2.h" 908 DIMENSION IADR1(NP,NH) 909#include "inforb.h" 910#include "infind.h" 911#include "infrsp.h" 912#include "wrkrsp.h" 913C 914 DIMENSION COEMP2(*),H2(NORBT,NORBT),A2TEMP(NORBT,NORBT) 915 DIMENSION COEUNP(NORBT,NORBT) 916C 917C C occupied, D virtual 918C 919 IF (IFRMP2(ICI) .NE. 0) RETURN 920 ICP = INDXPH(ICI) 921 IDP = -INDXPH(IDI) 922C unpack ph block in COEUNP (i.e. COEUNP(i,A) not defined) 923 CALL PHPACK(COEUNP,COEMP2(IADR1(IDP,ICP)+1),ICDSYM,-1,-1) 924 DO IASYM = 1, NSYM 925 IBSYM = MULD2H(IASYM,ICDSYM) 926 NSSHA = NSSH(IASYM) 927 NOCCB = NOCC(IBSYM) 928 IF ((NSSHA.NE.0).AND.(NOCCB.NE.0)) THEN 929C CONSTRUCT A2TEMP(IJ) TERMS 930 IAST = IORB(IASYM) + 1 + NOCC(IASYM) 931 IBST = IORB(IBSYM) + 1 932C we use H2 symmetric to get transpose H2 933 CALL DGEMM('T','N',NOCCB,NOCCB,NSSHA,1.D0, 934 & H2(IAST,IBST),NORBT, 935 & COEUNP(IAST,IBST),NORBT,1.D0, 936 & A2TEMP(IBST,IBST),NORBT) 937C CONSTRUCT A2TEMP(AB) TERMS 938 CALL DGEMM('N','N',NSSHA,NSSHA,NOCCB,1.D0, 939 & COEUNP(IAST,IBST),NORBT, 940 & H2(IBST,IAST),NORBT,1.D0, 941 & A2TEMP(IAST,IAST),NORBT) 942 ENDIF 943 ENDDO 944C 945C END OF HRPAA2 946C 947 RETURN 948 END 949C /* Deck hrpahm */ 950 SUBROUTINE HRPAHM(FCONE,A2TEMP,ZYMAT,NOSIM) 951C 952C This routine cancels the nonsymmetric parts of the A matrix 953C in SOPPA by using the explicitly constructed A(2) matrix 954C which is in A2TEMP (outside called XINDX). 955C Adapted by Erik K. Dalskov from HERM_A 956C written by Martin J. Packer. 957C 958#include "implicit.h" 959C 960#include "maxorb.h" 961#include "maxash.h" 962#include "infmp2.h" 963#include "inforb.h" 964#include "infind.h" 965#include "infrsp.h" 966#include "wrkrsp.h" 967C 968 DIMENSION FCONE(NORBT,NORBT,*) 969 DIMENSION A2TEMP(NORBT,NORBT), ZYMAT(NORBT,NORBT,*) 970C Now cancel the nonsymmetric parts of A(2) 971 DO I = 1,NOSIM 972 DO IASYM = 1,NSYM 973 IBSYM = MULD2H(IASYM,KSYMOP) 974 NOCCA = NOCC(IASYM) 975 NSSHB = NSSH(IBSYM) 976 IF ((NOCCA.NE.0).AND.(NSSHB.NE.0)) THEN 977C F(AI) TERMS 978 IAST = IORB(IASYM) + 1 979 IBST = IORB(IBSYM) + 1 + NOCC(IBSYM) 980C 981C sum(j)(-b(bj)*A(ai,bj) 982 CALL DGEMM('N','N',NSSHB,NOCCA,NOCCA,-1.D0, 983 & ZYMAT(IBST,IAST,I),NORBT, 984 & A2TEMP(IAST,IAST),NORBT,1.D0, 985 & FCONE(IBST,IAST,I),NORBT) 986C +b(bj)*A(bj,ai)) 987 CALL DGEMM('N','T',NSSHB,NOCCA,NOCCA,1.D0, 988 & ZYMAT(IBST,IAST,I),NORBT, 989 & A2TEMP(IAST,IAST),NORBT,1.D0, 990 & FCONE(IBST,IAST,I),NORBT) 991C sum(b)(-b(bj)*A(ai,bj) 992 CALL DGEMM('N','N',NSSHB,NOCCA,NSSHB,-1.D0, 993 & A2TEMP(IBST,IBST),NORBT, 994 & ZYMAT(IBST,IAST,I),NORBT,1.D0, 995 & FCONE(IBST,IAST,I),NORBT) 996C +b(bj)*A(bj,ai)) 997 CALL DGEMM('T','N',NSSHB,NOCCA,NSSHB,1.D0, 998 & A2TEMP(IBST,IBST),NORBT, 999 & ZYMAT(IBST,IAST,I),NORBT,1.D0, 1000 & FCONE(IBST,IAST,I),NORBT) 1001C F(IA) TERMS 1002C 1003C sum(b)(-b(jb)*A(jb,ia) 1004 CALL DGEMM('N','N',NOCCA,NSSHB,NSSHB,-1.D0, 1005 & ZYMAT(IAST,IBST,I),NORBT, 1006 & A2TEMP(IBST,IBST),NORBT,1.D0, 1007 & FCONE(IAST,IBST,I),NORBT) 1008C +b(jb)*A(ia,jb)) 1009 CALL DGEMM('N','T',NOCCA,NSSHB,NSSHB,1.D0, 1010 & ZYMAT(IAST,IBST,I),NORBT, 1011 & A2TEMP(IBST,IBST),NORBT,1.D0, 1012 & FCONE(IAST,IBST,I),NORBT) 1013C sum(j)(-b(jb)*A(jb,ia) 1014 CALL DGEMM('N','N',NOCCA,NSSHB,NOCCA,-1.D0, 1015 & A2TEMP(IAST,IAST),NORBT, 1016 & ZYMAT(IAST,IBST,I),NORBT,1.D0, 1017 & FCONE(IAST,IBST,I),NORBT) 1018C +b(jb)*A(ia,jb)) 1019 CALL DGEMM('T','N',NOCCA,NSSHB,NOCCA,1.D0, 1020 & A2TEMP(IAST,IAST),NORBT, 1021 & ZYMAT(IAST,IBST,I),NORBT,1.D0, 1022 & FCONE(IAST,IBST,I),NORBT) 1023 ENDIF 1024 ENDDO 1025 ENDDO 1026C 1027C END OF HRPAHM 1028C 1029 RETURN 1030 END 1031C /* Deck sopdiag */ 1032 SUBROUTINE SOPDIAG(KSYMOP,FC,DIAG,ABSADR,ABTADR,IJSADR,IJTADR, 1033 & IJ1ADR,IJ2ADR,IJ3ADR,WRK,LWRK,IPRINT) 1034C 1035C Copyright 11/4-1994 by Erik K. Dalskov 1036C 1037C Purpose: Construct D-matrix in SOPPA 1038C 1039#include "implicit.h" 1040#include "priunit.h" 1041C 1042C 1043#include "maxorb.h" 1044#include "maxash.h" 1045#include "infpri.h" 1046#include "inforb.h" 1047#include "infsop.h" 1048#include "infmp2.h" 1049#include "infrsp.h" 1050#include "infind.h" 1051C 1052 INTEGER ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 1053 & IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM), 1054 & IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM), 1055 & IJ3ADR(NISHT,NISHT,NSYM) 1056 REAL*8 DIAG(*),FC(*),WRK(*) 1057 INTEGER A, AA, B, BB, ABSYM, BSYM 1058C 1059 KFRSAV = 1 1060 KFREE = KFRSAV 1061 LFREE = LWRK 1062 CALL MEMGET2('REAL','EPSIL',KEPSIL,NORBT,WRK,KFREE,LFREE) 1063 JEPSIL = KEPSIL - 1 1064 IS = 0 1065 DO IM=1,NSYM 1066 IORBM = IORB(IM) 1067 NORBM = NORB(IM) 1068 IIORBM = IIORB(IM) 1069 IF (NORBM .NE. 0) THEN 1070 DO IJ=1,NORBM 1071 IL = IIORBM + (IJ*(IJ+1))/2 1072 WRK((JEPSIL+IS)+IJ) = FC(IL) 1073 ENDDO 1074 IS = IS + NORBM 1075 ENDIF 1076 ENDDO 1077 IF (IPRINT .GT. 100) THEN 1078 write (lupri,'(/A)') 'Orbital energies in SOPDIAG:' 1079 write (lupri,'(10F13.6)') (WRK(JEPSIL+I),I=1,NORBT) 1080 write (lupri,'(/A)') 'Orbitals frozen in MP2 in SOPDIAG:' 1081 write (lupri,'(4(3X,5I3))') (IFRMP2(I),I=1,NORBT) 1082 write (lupri,'(/A,I4)') '2p-2h E[2] diagonal, symmetry',KSYMOP 1083 IF (TRPLET) THEN 1084 ELSE 1085 write (lupri,'(/A)') 'IJSADR :' 1086 call ioutput(IJSADR(1,1,KSYMOP),1,NISHT,1,NISHT, 1087 & NISHT,NISHT,-1,LUPRI) 1088 write (lupri,'(/A)') 'IJTADR :' 1089 call ioutput(IJTADR(1,1,KSYMOP),1,NISHT,1,NISHT, 1090 & NISHT,NISHT,-1,LUPRI) 1091 write (lupri,'(/A)') 'ABSADR :' 1092 call ioutput(ABSADR(1,1),1,NORBT,1,NORBT, 1093 & NORBT,NORBT,-1,LUPRI) 1094 write (lupri,'(/A)') 'ABTADR :' 1095 call ioutput(ABTADR(1,1),1,NORBT,1,NORBT, 1096 & NORBT,NORBT,-1,LUPRI) 1097 END IF 1098 END IF 1099C 1100 DO II=1,NH 1101 I = IPHORD(II) 1102 IF (IFRMP2(I) .NE. 0) GO TO 100 1103 DO JJ=1,II 1104 J = IPHORD(JJ) 1105 IF (IFRMP2(J) .NE. 0) GO TO 200 1106 IJSYM = MULD2H(ISMO(I),ISMO(J)) 1107 ABSYM = MULD2H(IJSYM,KSYMOP) 1108 ENIJ = WRK(JEPSIL+I) + WRK(JEPSIL+J) 1109 IF (TRPLET) THEN 1110 IJ1OFF = IJ1ADR(II,JJ,KSYMOP) 1111 IJ2OFF = IJ2ADR(II,JJ,KSYMOP) 1112 IJ3OFF = IJ3ADR(II,JJ,KSYMOP) 1113 ELSE 1114 IJSOFF = IJSADR(II,JJ,KSYMOP) 1115 IJTOFF = IJTADR(II,JJ,KSYMOP) 1116 IF (IPRINT .GT. 100) THEN 1117 write (lupri,'(A,5I5,2I20)') 1118 & 'II,I,JJ,J,IJSYM,IJSOFF,IJTOFF', 1119 & II,I,JJ,J,IJSYM,IJSOFF,IJTOFF 1120 END IF 1121 END IF 1122 DO AA=1,NP 1123 A = IPHORD(NH+AA) 1124 ENAIJ = WRK(JEPSIL+A) - ENIJ 1125 BSYM = MULD2H(ISMO(A),ABSYM) 1126 IF (BSYM .GT. ISMO(A)) GO TO 300 1127 DO BB=1,AA 1128 B = IPHORD(NH+BB) 1129 IF (ISMO(B) .EQ. BSYM) THEN 1130 ENABIJ = ENAIJ + WRK(JEPSIL+B) 1131 IF (TRPLET) THEN 1132 NABT = ABTADR(AA,BB) 1133 IF (I .NE. J) THEN 1134 IF (A .NE. B) THEN 1135 DIAG(IJ1OFF+NABT)=ENABIJ 1136 DIAG(IJ3OFF+NABT)=ENABIJ 1137 ENDIF 1138 NABS = ABSADR(AA,BB) 1139 DIAG(IJ2OFF+NABS)=ENABIJ 1140 ELSE 1141 IF (A .NE. B) THEN 1142 DIAG(IJ3OFF+NABT)=ENABIJ 1143 END IF 1144 END IF 1145 ELSE 1146 NABS = ABSADR(AA,BB) 1147 DIAG(IJSOFF+NABS) = ENABIJ 1148 IF ((I .NE. J) .AND. (A .NE. B)) THEN 1149 NABT = ABTADR(AA,BB) 1150 DIAG(IJTOFF+NABT) = ENABIJ 1151 ELSE 1152 NABT = -1 1153 ENDIF 1154 IF (IPRINT .GT. 100) THEN 1155 write (lupri,'(A,5I5,2I20,F15.8)') 1156 & 'AA,A,BB,B,ABSYM,NABS,NABT,ENABIJ', 1157 & AA,A,BB,B,ABSYM,NABS,NABT,ENABIJ 1158 END IF 1159 ENDIF 1160 ENDIF 1161 ENDDO 1162 300 CONTINUE 1163 ENDDO 1164 200 CONTINUE 1165 ENDDO 1166 100 CONTINUE 1167 ENDDO 1168 CALL MEMREL('SOPDIAG',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1169C 1170C End of SOPDIAG 1171C 1172 RETURN 1173 END 1174C /* Deck soppaf */ 1175 SUBROUTINE SOPPAF(FVTD,H2,ZYCVEC,IJ,IB,JBSYM,NCSIM,XINDX, 1176 * WRK,LWRK) 1177C 1178C Purpose: This routine calculates the product of a 1179C 2p-2h trial vector with a C-matrix for SOPPA. 1180C 1181C Copyright 11/7-1994 by Erik K. Dalskov 1182C 1183#include "implicit.h" 1184C 1185C 1186#include "maxorb.h" 1187#include "maxash.h" 1188#include "wrkrsp.h" 1189#include "infpri.h" 1190#include "inforb.h" 1191#include "infmp2.h" 1192#include "infsop.h" 1193#include "infind.h" 1194#include "scbrhf.h" 1195C 1196 DIMENSION H2(NORBT,*),FVTD(NORBT,NORBT,*),ZYCVEC(KZCONF,*),WRK(*) 1197 DIMENSION XINDX(*) 1198 INTEGER CSYM, ASYM 1199C 1200C skip this distribution if J is frozen in MP2 1201 IF (IFRMP2(IJ) .NE. 0) GO TO 9999 1202C 1203 KFRSAV = 1 1204 KFREE = KFRSAV 1205 LFREE = LWRK 1206 CALL MEMGET('REAL',KUNP,NP*NH,WRK,KFREE,LFREE) 1207C 1208 DO ISIM = 1,NCSIM 1209 CALL GET2PH(WRK(KUNP),ZYCVEC(1,ISIM),XINDX(KABSAD), 1210 & XINDX(KABTAD),XINDX(KIJSAD),XINDX(KIJTAD), 1211 & XINDX(KIJ1AD),XINDX(KIJ2AD),XINDX(KIJ3AD), 1212 & IJ,IB,JBSYM) 1213 DO KSYM=1,NSYM 1214 CSYM = MULD2H(KSYM,KSYMOP) 1215 NOCCK = NOCC(KSYM) - NFRRHF(KSYM) 1216 NSSHC = NSSH(CSYM) 1217 IF ((NOCCK.NE.0).AND.(NSSHC.NE.0)) THEN 1218 IKST = IORB(KSYM) + NFRRHF(KSYM) + 1 1219 ICST = IORB(CSYM) + NOCC(CSYM) + 1 1220C 1221C F(c,k) = - SUM(i)[S(c,i)*H2(i,k)] 1222C NB! in this summation k may be frozen in MP2 1223C 1224 ISYM = MULD2H(KSYM,JBSYM) 1225 NOCCI = NOCC(ISYM) - NFRMP2(ISYM) 1226 IF (NOCCI .GT. 0) THEN 1227 IIST = IORB(ISYM) + NFRMP2(ISYM) + 1 1228 ICIST = KUNP + ( IOCC(ISYM) + NFRMP2(ISYM) ) * NP 1229 CALL DGEMM('N','N',NSSHC,NOCCK,NOCCI,-1.D0, 1230 & WRK(ICIST),NP, 1231 & H2(IIST,IKST),NORBT,1.D0, 1232 & FVTD(ICST,IKST,ISIM),NORBT) 1233 END IF 1234C 1235C + SUM(a)[H2(c,a)*S(a,k)] 1236C NB! in this summation k must not be frozen in MP2 1237C 1238 ASYM = MULD2H(CSYM,JBSYM) 1239 NSSHA = NSSH(ASYM) 1240 NOCCK = NOCC(KSYM) - NFRMP2(KSYM) 1241 IF (NSSHA .GT. 0 .AND. NOCCK.GT.0) THEN 1242 IAST = IORB(ASYM) + NOCC(ASYM) + 1 1243 IAKST = KUNP + ( IOCC(KSYM) + NFRMP2(KSYM) ) * NP 1244 IKST = IORB(KSYM) + NFRMP2(KSYM) + 1 1245 CALL DGEMM('N','N',NSSHC,NOCCK,NSSHA,1.D0, 1246 & H2(ICST,IAST),NORBT, 1247 & WRK(IAKST),NP,1.D0, 1248 & FVTD(ICST,IKST,ISIM),NORBT) 1249 END IF 1250 ENDIF 1251 ENDDO 1252 ENDDO 1253C 1254C END OF SOPPAF 1255C 1256 CALL MEMREL('SOPPAF',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1257 9999 RETURN 1258 END 1259C /* Deck get2ph */ 1260 SUBROUTINE GET2PH(UNPS,ZYCVEC,ABSADR,ABTADR,IJSADR,IJTADR,IJ1ADR, 1261 & IJ2ADR,IJ3ADR,JJ,BB,JBSYM) 1262C 1263C This routine finds the 2p-2h trial vector elements to be used in 1264C routine SOPPAF. If a .ne. b and i .ne. j then we have that an 1265C unpacked trial matrix are S(1) + ROOT3 * S(2), where the 1 and 2 1266C denotes the singlet and triplet coupled elements, respectively. 1267C 1268C Copyright 11/7-1994 by Erik K. Dalskov 1269C 1270C Triplet implemented 22 Feb. 1995 by Thomas Enevoldsen (tec) 1271C 1272#include "implicit.h" 1273C 1274C 1275#include "maxorb.h" 1276#include "maxash.h" 1277#include "wrkrsp.h" 1278#include "infpri.h" 1279#include "inforb.h" 1280#include "infsop.h" 1281#include "infmp2.h" 1282#include "infind.h" 1283#include "infrsp.h" 1284C 1285 DIMENSION UNPS(NP,*),ZYCVEC(*) 1286 DIMENSION ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT) 1287 DIMENSION IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM) 1288 DIMENSION IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM) 1289 DIMENSION IJ3ADR(NISHT,NISHT,NSYM) 1290 INTEGER A,AA,B,BB,ASYM,AOFF,ABSADR,ABTADR 1291C 1292 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, D3 = 3.0D0 , DP5 = 0.5D0 ) 1293C 1294 LOGICAL FIRST 1295 SAVE FIRST,ROOT3H,ROOT2I 1296 DATA FIRST /.TRUE./ 1297C 1298 IF (FIRST) THEN 1299 FIRST = .FALSE. 1300 ROOT3H = SQRT(D3) * DP5 1301 ROOT2I = SQRT(DP5) 1302 END IF 1303C 1304 J = INDXPH(JJ) 1305 B = -INDXPH(BB) 1306 IASYM = MULD2H(KSYMOP,JBSYM) 1307C Test for triplet or singlet 1308 IF (.NOT. TRPLET) THEN 1309 DO ISYM=1,NSYM 1310 ASYM = MULD2H(ISYM,IASYM) 1311 IOFF = IORB(ISYM) 1312 AOFF = IORB(ASYM) + NOCC(ASYM) 1313 DO II=1,NOCC(ISYM) 1314 IF (IFRMP2(IOFF + II) .NE. 0) GO TO 100 1315 I = INDXPH(IOFF + II) 1316 IJSOFF = IJSADR(I,J,KSYMOP) 1317 IJTOFF = IJTADR(I,J,KSYMOP) 1318 IF (I .GE. J) THEN 1319 FACIJ = ROOT3H 1320 ELSE 1321 FACIJ = -ROOT3H 1322 ENDIF 1323 DO AA=1,NSSH(ASYM) 1324 A = -INDXPH(AOFF + AA) 1325 NIJABS = IJSOFF + ABSADR(A,B) 1326 IF (I .NE. J) THEN 1327 IF (A .NE. B) THEN 1328 IF (A .GE. B) THEN 1329 FAC = FACIJ 1330 ELSE 1331 FAC = -FACIJ 1332 ENDIF 1333 NIJABT = IJTOFF + ABTADR(A,B) 1334 UNPS(AA,I) = ZYCVEC(NIJABS)*DP5 1335 * + ZYCVEC(NIJABT)*FAC 1336 ELSE 1337 UNPS(AA,I) = ZYCVEC(NIJABS)*ROOT2I 1338 ENDIF 1339 ELSE 1340 IF (A .NE. B) THEN 1341 UNPS(AA,I) = ZYCVEC(NIJABS)*ROOT2I 1342 ELSE 1343 UNPS(AA,I) = ZYCVEC(NIJABS) 1344 ENDIF 1345 ENDIF 1346 ENDDO 1347 100 CONTINUE 1348 ENDDO 1349 ENDDO 1350 ELSE 1351C If TRPLET 1352 DO ISYM=1,NSYM 1353 ASYM = MULD2H(ISYM,IASYM) 1354 IOFF = IORB(ISYM) 1355 AOFF = IORB(ASYM) + NOCC(ASYM) 1356 DO II=1,NOCC(ISYM) 1357 IF (IFRMP2(IOFF + II) .NE. 0) GO TO 110 1358 I = INDXPH(IOFF + II) 1359 IJ1OFF = IJ1ADR(I,J,KSYMOP) 1360 IJ2OFF = IJ2ADR(I,J,KSYMOP) 1361 IJ3OFF = IJ3ADR(I,J,KSYMOP) 1362 IF (I .GE. J) THEN 1363 TIJSGN = D1 1364 ELSE 1365 TIJSGN = -D1 1366 ENDIF 1367 IF (I .EQ. J) THEN 1368 C3FAC = ROOT2I 1369 ELSE 1370 C3FAC = DP5 1371 ENDIF 1372 DO AA=1,NSSH(ASYM) 1373 A = -INDXPH(AOFF + AA) 1374 NABS = ABSADR(A,B) 1375 NABT = ABTADR(A,B) 1376 IF (A .GE. B) THEN 1377 TABSGN = D1 1378 ELSE 1379 TABSGN = -D1 1380 ENDIF 1381 IF (A .EQ. B) THEN 1382 C2FAC = ROOT2I 1383 ELSE 1384 C2FAC = DP5 1385 ENDIF 1386C C(1) CONTRIBUTION 1387 IF ((I .NE. J) .AND. (A .NE. B)) THEN 1388 UNPS(AA,I) = TABSGN * TIJSGN 1389 * * ROOT2I * ZYCVEC(IJ1OFF + NABT) 1390 ELSE 1391 UNPS(AA,I) = D0 1392 ENDIF 1393C C(2) CONTRIBUTION 1394 IF (I .NE. J) THEN 1395 UNPS(AA,I) = UNPS(AA,I) 1396 * - C2FAC * TIJSGN * ZYCVEC(IJ2OFF + NABS) 1397 ENDIF 1398C C(3) CONTRIBUTION 1399 IF (A .NE. B) THEN 1400 UNPS(AA,I) = UNPS(AA,I) 1401 * - C3FAC * TABSGN * ZYCVEC(IJ3OFF + NABT) 1402 ENDIF 1403 ENDDO 1404 110 CONTINUE 1405 ENDDO 1406 ENDDO 1407 ENDIF 1408C 1409C END OF GET2PH 1410C 1411 RETURN 1412 END 1413C /* Deck soph2x */ 1414 SUBROUTINE SOPH2X(EVECS,ZYMAT,TZYMAT,H2,IJ,IB,JBSYM,NOSIM, 1415 & ABSADR,ABTADR,IJSADR,IJTADR,IJ1ADR,IJ2ADR, 1416 & IJ3ADR,WRK,LWRK) 1417C 1418C This routine calculates ph trial vector times a C-matrix 1419C in SOPPA. 1420C 1421C Copyright 11/7-1994 by Erik K. Dalskov 1422C Revision 940711/940726 hjaaj 1423C Triplet case implementet 21 feb. 1995 by Thomas Enevoldsen (tec) 1424#include "implicit.h" 1425C 1426C 1427#include "maxorb.h" 1428#include "maxash.h" 1429#include "wrkrsp.h" 1430#include "infpri.h" 1431#include "inforb.h" 1432#include "infsop.h" 1433#include "infmp2.h" 1434#include "infind.h" 1435#include "infrsp.h" 1436C 1437C 1438 DIMENSION EVECS(KZYVAR,*),H2(NORBT,*),WRK(*) 1439 DIMENSION ZYMAT(NORBT,NORBT,*),TZYMAT(NORBT,NORBT,*) 1440 DIMENSION ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT) 1441 DIMENSION IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM) 1442 DIMENSION IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM) 1443 DIMENSION IJ3ADR(NISHT,NISHT,NSYM) 1444 INTEGER A,AA,B,ASYM,CSYM, ABSADR, ABTADR 1445C 1446 PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, D3 = 3.0D0 ) 1447C 1448 LOGICAL FIRST 1449 SAVE FIRST,ROOT2,ROOT3 1450 DATA FIRST /.TRUE./ 1451C 1452 IF (FIRST) THEN 1453 FIRST = .FALSE. 1454 ROOT2 = SQRT(D2) 1455 ROOT3 = SQRT(D3) 1456 ENDIF 1457C 1458C skip this distribution if J is frozen in MP2 1459 IF (IFRMP2(IJ) .NE. 0) GO TO 9999 1460C 1461 KFRSAV = 1 1462 KFREE = KFRSAV 1463 LFREE = LWRK 1464 CALL MEMGET('REAL',KB1,NP*NH,WRK,KFREE,LFREE) 1465 CALL MEMGET('REAL',KB2,NP*NH,WRK,KFREE,LFREE) 1466 LB1 = KB1 - 1 1467 LB2 = KB2 - 1 1468C 1469C 1470 J = INDXPH(IJ) 1471 B = -INDXPH(IB) 1472 IASYM = MULD2H(KSYMOP,JBSYM) 1473 DO ISYM=1,NSYM 1474 ASYM = MULD2H(ISYM,IASYM) 1475 NOCCI = NOCC(ISYM) - NFRMP2(ISYM) 1476 NSSHA = NSSH(ASYM) 1477 IF (NOCCI .GT. 0 .AND. NSSHA .GT. 0) THEN 1478 IIST = IORB(ISYM) + NFRMP2(ISYM) + 1 1479 IAST = IORB(ASYM) + NOCC(ASYM) + 1 1480C 1481 KSYM = MULD2H(ISYM,JBSYM) 1482 NOCCK = NOCC(KSYM) 1483 CSYM = MULD2H(ASYM,JBSYM) 1484 NSSHC = NSSH(CSYM) 1485 IF (NOCCK .GT. 0 .OR. NSSHC .GT. 0) THEN 1486 DO ISIM = 1,NOSIM 1487 IF (NOCCK .EQ. 0) THEN 1488 CALL DZERO(WRK(KB1),NSSHA*NOCCI) 1489 CALL DZERO(WRK(KB2),NSSHA*NOCCI) 1490 ELSE 1491 IKST = IORB(KSYM) + 1 1492 CALL DGEMM('N','N',NSSHA,NOCCI,NOCCK,1.D0, 1493 & ZYMAT(IAST,IKST,ISIM),NORBT, 1494 & H2(IKST,IIST),NORBT,0.D0, 1495 & WRK(KB1),NSSHA) 1496C 1497 CALL DGEMM('N','N',NSSHA,NOCCI,NOCCK,1.D0, 1498 & TZYMAT(IAST,IKST,ISIM),NORBT, 1499 & H2(IKST,IIST),NORBT,0.D0, 1500 & WRK(KB2),NSSHA) 1501 ENDIF 1502 IF (NSSHC .GT. 0) THEN 1503 ICST = IORB(CSYM) + NOCC(CSYM) + 1 1504 CALL DGEMM('N','N',NSSHA,NOCCI,NSSHC,-1.D0, 1505 & H2(IAST,ICST),NORBT, 1506 & ZYMAT(ICST,IIST,ISIM),NORBT,1.D0, 1507 & WRK(KB1),NSSHA) 1508C 1509 CALL DGEMM('N','N',NSSHA,NOCCI,NSSHC,-1.D0, 1510 & H2(IAST,ICST),NORBT, 1511 & TZYMAT(ICST,IIST,ISIM),NORBT,1.D0, 1512 & WRK(KB2),NSSHA) 1513 ENDIF 1514 IF (.NOT. TRPLET) THEN 1515C SINGLET CASE 1516 DO II=1,NOCCI 1517 I = INDXPH(IIST + II - 1) 1518 IJSOFF = IJSADR(I,J,KSYMOP) 1519 IJTOFF = IJTADR(I,J,KSYMOP) 1520 IF (I .GE. J) THEN 1521 FACIJ = ROOT3 1522 ELSE 1523 FACIJ = -ROOT3 1524 ENDIF 1525 DO AA=1,NSSHA 1526 A = -INDXPH(IAST + AA - 1) 1527 IF (A .GE. B) THEN 1528 FAC = FACIJ 1529 ELSE 1530 FAC = - FACIJ 1531 ENDIF 1532 NABS = ABSADR(A,B) 1533 NABT = ABTADR(A,B) 1534 IF (I .NE. J) THEN 1535 IF (A .NE. B) THEN 1536 EVECS(IJSOFF+NABS,ISIM) = 1537 * EVECS(IJSOFF+NABS,ISIM) + 1538 * WRK(LB1+(II-1)*NSSHA+AA) 1539 EVECS(KZVAR+IJSOFF+NABS,ISIM) = 1540 * EVECS(KZVAR+IJSOFF+NABS,ISIM) + 1541 * WRK(LB2+(II-1)*NSSHA+AA) 1542 EVECS(IJTOFF+NABT,ISIM) = 1543 * EVECS(IJTOFF+NABT,ISIM) + 1544 * WRK(LB1+(II-1)*NSSHA+AA) * FAC 1545 EVECS(KZVAR+IJTOFF+NABT,ISIM) = 1546 * EVECS(KZVAR+IJTOFF+NABT,ISIM) + 1547 * WRK(LB2+(II-1)*NSSHA+AA) * FAC 1548 ELSE 1549 EVECS(IJSOFF+NABS,ISIM) = 1550 * EVECS(IJSOFF+NABS,ISIM) + 1551 * WRK(LB1+(II-1)*NSSHA+AA) * ROOT2 1552 EVECS(KZVAR+IJSOFF+NABS,ISIM) = 1553 * EVECS(KZVAR+IJSOFF+NABS,ISIM) + 1554 * WRK(LB2+(II-1)*NSSHA+AA) * ROOT2 1555 ENDIF 1556 ELSE 1557 IF (A .NE. B) THEN 1558 EVECS(IJSOFF+NABS,ISIM) = 1559 * EVECS(IJSOFF+NABS,ISIM) + 1560 * WRK(LB1+(II-1)*NSSHA+AA) * ROOT2 1561 EVECS(KZVAR+IJSOFF+NABS,ISIM) = 1562 * EVECS(KZVAR+IJSOFF+NABS,ISIM) + 1563 * WRK(LB2+(II-1)*NSSHA+AA) * ROOT2 1564 ELSE 1565 EVECS(IJSOFF+NABS,ISIM) = 1566 * EVECS(IJSOFF+NABS,ISIM) + 1567 * WRK(LB1+(II-1)*NSSHA+AA) * D2 1568 EVECS(KZVAR+IJSOFF+NABS,ISIM) = 1569 * EVECS(KZVAR+IJSOFF+NABS,ISIM) + 1570 * WRK(LB2+(II-1)*NSSHA+AA) * D2 1571 ENDIF 1572 ENDIF 1573 ENDDO 1574 ENDDO 1575 ELSE 1576C Triplet case 1577 DO II=1,NOCCI 1578 I = INDXPH(IIST + II - 1) 1579 IJ1OFF = IJ1ADR(I,J,KSYMOP) 1580 IJ2OFF = IJ2ADR(I,J,KSYMOP) 1581 IJ3OFF = IJ3ADR(I,J,KSYMOP) 1582 IF (I .GE. J) THEN 1583 TIJSGN = D1 1584 ELSE 1585 TIJSGN = -D1 1586 ENDIF 1587 IF (I .EQ. J) THEN 1588 C3FAC = ROOT2 1589 ELSE 1590 C3FAC = D1 1591 ENDIF 1592 DO AA=1,NSSHA 1593 A = -INDXPH(IAST + AA - 1) 1594 IF (A .GE. B) THEN 1595 TABSGN = D1 1596 ELSE 1597 TABSGN = -D1 1598 ENDIF 1599 IF (A .EQ. B) THEN 1600 C2FAC = ROOT2 1601 ELSE 1602 C2FAC = D1 1603 ENDIF 1604 NABS = ABSADR(A,B) 1605 NABT = ABTADR(A,B) 1606C C(1) PART OF VECTOR 1607 IF ((I .NE. J) .AND. (A .NE. B)) THEN 1608 EVECS(IJ1OFF+NABT,ISIM) = 1609 * EVECS(IJ1OFF+NABT,ISIM) + 1610 * WRK(LB1+(II-1)*NSSHA+AA) * TIJSGN 1611 * * TABSGN * ROOT2 1612 EVECS(KZVAR+IJ1OFF+NABT,ISIM) = 1613 * EVECS(KZVAR+IJ1OFF+NABT,ISIM) + 1614 * WRK(LB2+(II-1)*NSSHA+AA) * TIJSGN 1615 * * TABSGN * ROOT2 1616 ENDIF 1617C C(2) PART OF VECTOR 1618 IF (I .NE. J) THEN 1619 EVECS(IJ2OFF+NABS,ISIM) = 1620 * EVECS(IJ2OFF+NABS,ISIM) - 1621 * WRK(LB1+(II-1)*NSSHA+AA)*TIJSGN*C2FAC 1622 EVECS(KZVAR+IJ2OFF+NABS,ISIM) = 1623 * EVECS(KZVAR+IJ2OFF+NABS,ISIM) - 1624 * WRK(LB2+(II-1)*NSSHA+AA)*TIJSGN*C2FAC 1625 ENDIF 1626C C(3) PART OF VECTOR 1627 IF (A .NE. B) THEN 1628 EVECS(IJ3OFF+NABT,ISIM) = 1629 * EVECS(IJ3OFF+NABT,ISIM) - 1630 * WRK(LB1+(II-1)*NSSHA+AA)*TABSGN*C3FAC 1631 EVECS(KZVAR+IJ3OFF+NABT,ISIM) = 1632 * EVECS(KZVAR+IJ3OFF+NABT,ISIM) - 1633 * WRK(LB2+(II-1)*NSSHA+AA)*TABSGN*C3FAC 1634 ENDIF 1635 ENDDO 1636 ENDDO 1637 ENDIF 1638 END DO 1639 ENDIF 1640 ENDIF 1641 ENDDO 1642C 1643 CALL MEMREL('SOPH2X',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1644C 1645C END OF SOPH2X 1646C 1647 9999 RETURN 1648 END 1649C /* Deck sopcli */ 1650 SUBROUTINE SOPCLI(NCSIM,ZYCVEC,EVECS,SOPCL,LWRK) 1651C 1652C This routine calculates D(0)*S, where D(0) is the D-matrix in 1653C SOPPA and S is a 2p-2h trial vector 1654C 1655C Copyright 11/7-1994 by Erik K. Dalskov 1656C 1657C 1658#include "implicit.h" 1659C 1660C 1661#include "wrkrsp.h" 1662#include "inftap.h" 1663#include "infrsp.h" 1664C 1665 DIMENSION EVECS(KZYVAR,*),ZYCVEC(KZCONF,*),SOPCL(*) 1666C 1667 IF (KZCONF .GT. LWRK) CALL ERRWRK('SOPCLI',KZCONF,LWRK) 1668 REWIND (LURSP4) 1669 CALL READT(LURSP4,KZCONF,SOPCL) 1670C 1671C 1672C 1673 DO ISIM=1,NCSIM 1674 DO I=1,KZCONF 1675 EVECS(I,ISIM) = EVECS(I,ISIM) + ZYCVEC(I,ISIM)*SOPCL(I) 1676 ENDDO 1677 ENDDO 1678 RETURN 1679 END 1680C /* Deck sopijab */ 1681 SUBROUTINE SOPIJAB(IG,ISYMV,II,ISYM,JJ,JSYM,AA,ASYM,BB, 1682 & BSYM,STTYP,ABSADR,ABTADR,IJSADR,IJTADR, 1683 & IJ1ADR,IJ2ADR,IJ3ADR) 1684C 1685C This routine finds i,j,a,b-indices corresponding to a specific 1686C element IG in a 2p-2h vector 1687C 1688C Copyright 11/7-94 by Erik K. Dalskov 1689C 1690#include "implicit.h" 1691#include "dummy.h" 1692C 1693C 1694#include "maxorb.h" 1695#include "maxash.h" 1696#include "priunit.h" 1697C 1698 DIMENSION ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 1699 & IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM), 1700 & IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM), 1701 & IJ3ADR(NISHT,NISHT,NSYM) 1702#include "infpri.h" 1703#include "inforb.h" 1704#include "infsop.h" 1705#include "infmp2.h" 1706#include "infind.h" 1707#include "infrsp.h" 1708C 1709C 1710 INTEGER A,AA,B,BB,C,CC,D,DD,CSYM,DSYM,ABSYM,DIFF,DISC,ASYM,BSYM, 1711 & ABSADR,ABTADR 1712 CHARACTER*1 STTYP 1713C 1714 A = IDUMMY 1715 B = IDUMMY 1716 IF (.NOT. TRPLET) THEN 1717 IF (IG .GT. NS2P2H(ISYMV)) THEN 1718 STTYP = 'T' 1719 KOLD = 0 1720 LOLD = 0 1721 DO K=2,NH 1722 IF (IFRMP2(IPHORD(K)) .EQ. 0) THEN 1723 DO L=1,K-1 1724 IF (IFRMP2(IPHORD(L)) .EQ. 0) THEN 1725 DIFF = IG - IJTADR(L,K,ISYMV) 1726 IF (DIFF .LE. 0) THEN 1727 IF (L .GT. LOLD) THEN 1728 I = K 1729 J = LOLD 1730 ELSE 1731 I = KOLD 1732 J = LOLD 1733 ENDIF 1734 GO TO 120 1735 ENDIF 1736 LOLD = L 1737 END IF 1738 ENDDO 1739 KOLD = K 1740 END IF 1741 ENDDO 1742 I = KOLD 1743 J = LOLD 1744C 1745 120 II = IPHORD(I) 1746 JJ = IPHORD(J) 1747 IJSYM = MULD2H(ISMO(II),ISMO(JJ)) 1748 ABSYM = MULD2H(IJSYM,ISYMV) 1749 DISC = IG - IJTADR(I,J,ISYMV) 1750 DO C =2,NP 1751 CC = IPHORD(NH+C) 1752 CSYM = ISMO(CC) 1753 DO D=1,C-1 1754 DD = IPHORD(NH+D) 1755 DSYM = ISMO(DD) 1756 IF (ABSYM .EQ. MULD2H(DSYM,CSYM)) THEN 1757 IF (DISC .EQ. ABTADR(D,C)) THEN 1758 A = C 1759 B = D 1760 GO TO 140 1761 ENDIF 1762 ENDIF 1763 ENDDO 1764 ENDDO 1765 WRITE (LUPRI,*) 'SOPIJAB singlet type 2 ERROR CODE 140' 1766 WRITE (LUPRI,*) 'I,J,II,JJ,IJSYM,ABSYM',I,J,II,JJ,IJSYM, 1767 & ABSYM 1768 WRITE (LUPRI,*) 'IG,DISC',IG,DISC 1769 WRITE (LUPRI,*) 'ABTADR array:' 1770 CALL IWRTMA(ABTADR,NP,NP,NORBT,NORBT) 1771 WRITE (LUPRI,*) 'IJTADR(*,*,ISYMV) array:' 1772 CALL IWRTMA(IJTADR(1,1,ISYMV),NH,NH,MAXOCC,MAXOCC) 1773C CALL IWRTMA(IMAT,NROW,NCOL,MAXROW,MAXCOL) 1774 CALL QUIT('SOPIJAB ERROR CODE 140') 1775 140 CONTINUE 1776 ELSE 1777 STTYP = 'S' 1778 KOLD = 0 1779 LOLD = 0 1780 DO K=1,NH 1781 IF (IFRMP2(IPHORD(K)) .NE. 0) GO TO 333 1782 DO L=1,K 1783 IF (IFRMP2(IPHORD(L)) .NE. 0) GO TO 444 1784 DIFF = IG - IJSADR(L,K,ISYMV) 1785 IF (DIFF .LE. 0) THEN 1786 IF (L .GT. LOLD) THEN 1787 I = K 1788 J = LOLD 1789 ELSE 1790 I = KOLD 1791 J = LOLD 1792 ENDIF 1793 GO TO 110 1794 ENDIF 1795 LOLD = L 1796 444 CONTINUE 1797 ENDDO 1798 KOLD = K 1799 333 CONTINUE 1800 ENDDO 1801 I = KOLD 1802 J = LOLD 1803C 1804 110 II = IPHORD(I) 1805 JJ = IPHORD(J) 1806 IJSYM = MULD2H(ISMO(II),ISMO(JJ)) 1807 ABSYM = MULD2H(IJSYM,ISYMV) 1808 DISC = IG - IJSADR(I,J,ISYMV) 1809 DO C =1,NP 1810 CC = IPHORD(NH+C) 1811 CSYM = ISMO(CC) 1812 DO D=1,C 1813 DD = IPHORD(NH+D) 1814 DSYM = ISMO(DD) 1815 IF (ABSYM .EQ. MULD2H(DSYM,CSYM)) THEN 1816 IF (DISC .EQ. ABSADR(D,C)) THEN 1817 A = C 1818 B = D 1819 GO TO 130 1820 ENDIF 1821 ENDIF 1822 ENDDO 1823 ENDDO 1824 WRITE (LUPRI,*) 'SOPIJAB singlet type 1 ERROR CODE 130' 1825 WRITE (LUPRI,*) 'I,J,II,JJ',I,J,II,JJ 1826 WRITE (LUPRI,*) 'A,B,IJSYM,ABSYM',A,B,IJSYM,ABSYM 1827 WRITE (LUPRI,*) 'IG,DISC',IG,DISC 1828 WRITE (LUPRI,*) 'ABSADR array:' 1829 CALL IWRTMA(ABSADR,NP,NP,NORBT,NORBT) 1830 WRITE (LUPRI,*) 'IJSADR(*,*,ISYMV) array:' 1831 CALL IWRTMA(IJSADR(1,1,ISYMV),NH,NH,MAXOCC,MAXOCC) 1832C CALL IWRTMA(IMAT,NROW,NCOL,MAXROW,MAXCOL) 1833 CALL QUIT('SOPIJAB ERROR CODE 130') 1834 130 CONTINUE 1835 END IF 1836 ELSE 1837C TRIPLET CASE 1838 IF (IG .GT. (N12P2H(ISYMV)+N22P2H(ISYMV))) THEN 1839 STTYP = '3' 1840 KOLD = 0 1841 LOLD = 0 1842 DO K=1,NH 1843 IF (IFRMP2(IPHORD(K)) .NE. 0) GO TO 300 1844 DO L=1,K 1845 IF (IFRMP2(IPHORD(L)) .NE. 0) GO TO 400 1846 DIFF = IG - IJ3ADR(L,K,ISYMV) 1847 IF (DIFF .LE. 0) THEN 1848 IF (L.GT.LOLD) THEN 1849 I = K 1850 J = LOLD 1851 ELSE 1852 I = KOLD 1853 J = LOLD 1854 ENDIF 1855 GO TO 200 1856 ENDIF 1857 LOLD = L 1858 400 CONTINUE 1859 ENDDO 1860 KOLD = K 1861 300 CONTINUE 1862 ENDDO 1863 I = KOLD 1864 J = LOLD 1865C 1866 200 II = IPHORD(I) 1867 JJ = IPHORD(J) 1868 IJSYM = MULD2H(ISMO(II),ISMO(JJ)) 1869 ABSYM = MULD2H(IJSYM,ISYMV) 1870 DISC = IG - IJ3ADR(I,J,ISYMV) 1871 DO C =2,NP 1872 CC = IPHORD(NH+C) 1873 CSYM = ISMO(CC) 1874 DO D=1,C-1 1875 DD = IPHORD(NH+D) 1876 DSYM = ISMO(DD) 1877 IF (ABSYM .EQ. MULD2H(DSYM,CSYM)) THEN 1878 IF (DISC .EQ. ABTADR(D,C)) THEN 1879 A = C 1880 B = D 1881 GO TO 100 1882 ENDIF 1883 ENDIF 1884 ENDDO 1885 ENDDO 1886 WRITE (LUPRI,*) 'SOPIJAB triplet type 3 ERROR CODE 100' 1887 WRITE (LUPRI,*) 'I,J,II,JJ,IJSYM,ABSYM',I,J,II,JJ,IJSYM, 1888 & ABSYM 1889 WRITE (LUPRI,*) 'IG,DISC',IG,DISC 1890 WRITE (LUPRI,*) 'ABTADR array:' 1891 CALL IWRTMA(ABTADR,NP,NP,NORBT,NORBT) 1892 WRITE (LUPRI,*) 'IJ3ADR(*,*,ISYMV) array:' 1893 CALL IWRTMA(IJ3ADR(1,1,ISYMV),NH,NH,MAXOCC,MAXOCC) 1894C CALL IWRTMA(IMAT,NROW,NCOL,MAXROW,MAXCOL) 1895 CALL QUIT('SOPIJAB ERROR CODE 100') 1896 100 CONTINUE 1897 ELSE IF (IG .GT. N12P2H(ISYMV)) THEN 1898 STTYP = '2' 1899 KOLD = 0 1900 LOLD = 0 1901 DO K=2,NH 1902 IF (IFRMP2(IPHORD(K)) .EQ. 0) THEN 1903 DO L=1,K-1 1904 IF (IFRMP2(IPHORD(L)) .EQ. 0) THEN 1905 DIFF = IG - IJ2ADR(L,K,ISYMV) 1906 IF (DIFF .LE. 0) THEN 1907 IF (L .GT. LOLD) THEN 1908 I = K 1909 J = LOLD 1910 ELSE 1911 I = KOLD 1912 J = LOLD 1913 ENDIF 1914 GO TO 201 1915 ENDIF 1916 LOLD = L 1917 END IF 1918 ENDDO 1919 KOLD = K 1920 END IF 1921 ENDDO 1922 I = KOLD 1923 J = LOLD 1924C 1925C 1926 201 II = IPHORD(I) 1927 JJ = IPHORD(J) 1928 IJSYM = MULD2H(ISMO(II),ISMO(JJ)) 1929 ABSYM = MULD2H(IJSYM,ISYMV) 1930 DISC = IG - IJ2ADR(I,J,ISYMV) 1931 DO C =1,NP 1932 CC = IPHORD(NH+C) 1933 CSYM = ISMO(CC) 1934 DO D=1,C 1935 DD = IPHORD(NH+D) 1936 DSYM = ISMO(DD) 1937 IF (ABSYM .EQ. MULD2H(DSYM,CSYM)) THEN 1938 IF (DISC .EQ. ABSADR(D,C)) THEN 1939 A = C 1940 B = D 1941 GO TO 101 1942 ENDIF 1943 ENDIF 1944 ENDDO 1945 ENDDO 1946 WRITE (LUPRI,*) 'SOPIJAB triplet type 2 ERROR CODE 101' 1947 WRITE (LUPRI,*) 'I,J,II,JJ',I,J,II,JJ 1948 WRITE (LUPRI,*) 'A,B,IJSYM,ABSYM',A,B,IJSYM,ABSYM 1949 WRITE (LUPRI,*) 'IG,DISC',IG,DISC 1950 WRITE (LUPRI,*) 'ABSADR array:' 1951 CALL IWRTMA(ABSADR,NP,NP,NORBT,NORBT) 1952 WRITE (LUPRI,*) 'IJ2ADR(*,*,ISYMV) array:' 1953 CALL IWRTMA(IJ2ADR(1,1,ISYMV),NH,NH,MAXOCC,MAXOCC) 1954C CALL IWRTMA(IMAT,NROW,NCOL,MAXROW,MAXCOL) 1955 CALL QUIT('SOPIJAB ERROR CODE 101') 1956 101 CONTINUE 1957 ELSE 1958 STTYP = '1' 1959 KOLD = 0 1960 LOLD = 0 1961 DO K=2,NH 1962 IF (IFRMP2(IPHORD(K)) .EQ. 0) THEN 1963 DO L=1,K-1 1964 IF (IFRMP2(IPHORD(L)) .EQ. 0) THEN 1965 DIFF = IG - IJ1ADR(L,K,ISYMV) 1966 IF (DIFF .LE. 0) THEN 1967 IF (L .GT. LOLD) THEN 1968 I = K 1969 J = LOLD 1970 ELSE 1971 I = KOLD 1972 J = LOLD 1973 ENDIF 1974 GO TO 202 1975 ENDIF 1976 LOLD = L 1977 END IF 1978 ENDDO 1979 KOLD = K 1980 END IF 1981 ENDDO 1982 I = KOLD 1983 J = LOLD 1984C 1985 202 II = IPHORD(I) 1986 JJ = IPHORD(J) 1987 IJSYM = MULD2H(ISMO(II),ISMO(JJ)) 1988 ABSYM = MULD2H(IJSYM,ISYMV) 1989 DISC = IG - IJ1ADR(I,J,ISYMV) 1990 DO C =2,NP 1991 CC = IPHORD(NH+C) 1992 CSYM = ISMO(CC) 1993 DO D=1,C-1 1994 DD = IPHORD(NH+D) 1995 DSYM = ISMO(DD) 1996 IF (ABSYM .EQ. MULD2H(DSYM,CSYM)) THEN 1997 IF (DISC .EQ. ABTADR(D,C)) THEN 1998 A = C 1999 B = D 2000 GO TO 102 2001 ENDIF 2002 ENDIF 2003 ENDDO 2004 ENDDO 2005 WRITE (LUPRI,*) 'SOPIJAB triplet type 1 ERROR CODE 102' 2006 WRITE (LUPRI,*) 'I,J,II,JJ,IJSYM,ABSYM',I,J,II,JJ,IJSYM, 2007 & ABSYM 2008 WRITE (LUPRI,*) 'IG,DISC',IG,DISC 2009 WRITE (LUPRI,*) 'ABTADR array:' 2010 CALL IWRTMA(ABTADR,NP,NP,NORBT,NORBT) 2011 WRITE (LUPRI,*) 'IJ1ADR(*,*,ISYMV) array:' 2012 CALL IWRTMA(IJ1ADR(1,1,ISYMV),NH,NH,MAXOCC,MAXOCC) 2013C CALL IWRTMA(IMAT,NROW,NCOL,MAXROW,MAXCOL) 2014 CALL QUIT('SOPIJAB ERROR CODE 102') 2015 102 CONTINUE 2016 ENDIF 2017 ENDIF 2018 AA = IPHORD(NH+A) 2019 ASYM = ISMO(AA) 2020 BB = IPHORD(NH+B) 2021 BSYM = ISMO(BB) 2022 ISYM = ISMO(II) 2023 JSYM = ISMO(JJ) 2024C 2025C END OF SOPIJAB 2026C 2027 RETURN 2028 END 2029C /* Deck trpkap */ 2030 SUBROUTINE TRPKAP(COEMP2,TCOEMP,IADR1,IADR2) 2031C 2032C This routine makes the combination: 2033C Kt(ai,bj) = 1/3 ( K(ai,bj) + 2 K(aj,bi) ) 2034C and puts the result in PV(LPVMAT+1)=TCOEMP 2035C for use in SOPPA for triplet properties. 2036C 2037C Copyright 941122 Erik K. Dalskov 2038C 2039#include "implicit.h" 2040C 2041#include "maxorb.h" 2042#include "maxash.h" 2043#include "infmp2.h" 2044 DIMENSION IADR1(NP,NH) 2045#include "inforb.h" 2046#include "infind.h" 2047#include "infrsp.h" 2048#include "wrkrsp.h" 2049C 2050 INTEGER ASYM,BSYM,AST,AEND,BB1,B1,B1IOFF,AA,A,AIOFF1,BB 2051C 2052 DIMENSION COEMP2(*), TCOEMP(*), IADR2(NP,NH) 2053 PARAMETER ( D3I = 1.0D0/3.0D0 ) 2054C 2055 CALL DCOPY(LPVMAT,COEMP2(1),1,TCOEMP(1),1) 2056C 2057C call PHADR2 to get the index array IADR2 2058C 2059 CALL PHADR2(IADR2,NP,NH) 2060C 2061 DO IASYM=1,NSYM 2062 DO JASYM=1,NSYM 2063 DO I=1,NH 2064 II = IPHORD(I) 2065 ISYM = ISMO(II) 2066 ASYM = MULD2H(IASYM,ISYM) 2067 JSYM = MULD2H(JASYM,ASYM) 2068 BSYM = MULD2H(IASYM,JSYM) 2069 AST = IORB(ASYM) + NOCC(ASYM) + 1 2070 AEND = IORB(ASYM) + NORB(ASYM) 2071 JST = IORB(JSYM) + NFRMP2(JSYM) + 1 2072 JEND = IORB(JSYM) + NOCC(JSYM) 2073 BB1 = IORB(BSYM) + NOCC(BSYM) + 1 2074 B1 = -INDXPH(BB1) 2075 B1IOFF = IADR2(B1,I) - 1 2076 NSSHB= NSSH(BSYM) 2077 DO AA=AST,AEND 2078 A = - INDXPH(AA) 2079 AIOFF1 = IADR1(A,I) - 1 2080 DO JJ=JST,JEND 2081 J = INDXPH(JJ) 2082 IAJOFF = IADR1(A,J) + B1IOFF 2083 JAIOFF = AIOFF1 + IADR2(B1,J) 2084 DO BB=1,NSSHB 2085 TCOEMP(JAIOFF + BB) = 2086 & TCOEMP(JAIOFF + BB) + 2 * COEMP2(IAJOFF + BB) 2087 END DO 2088 END DO 2089 END DO 2090 END DO 2091 END DO 2092 END DO 2093 CALL DSCAL(LPVMAT,-D3I,TCOEMP,1) 2094C 2095C END OF TRPKAP 2096C 2097 RETURN 2098 END 2099C /* Deck sopdw4 */ 2100 SUBROUTINE SOPDW4(GP,DINVGP,DIAG,EFREQ) 2101C 2102C 10-Mar-1995 Hans Joergen Aa. Jensen 2103C 2104C This routine calculates [D(0) +/- EFREQ]**(-1) * GP 2105C where D(0) is the D-matrix in SOPPA and 2106C GP is the 2p-2h property gradient vector 2107C 2108C 2109#include "implicit.h" 2110C 2111C 2112#include "wrkrsp.h" 2113#include "inftap.h" 2114#include "infrsp.h" 2115C 2116 DIMENSION GP(KZYVAR), DINVGP(KZYVAR), DIAG(*) 2117C 2118 REWIND (LURSP4) 2119 CALL READT(LURSP4,KZCONF,DIAG) 2120C 2121C Zero ph part of DINVGP 2122C 2123 CALL DZERO(DINVGP(KZCONF+1),KZWOPT) 2124 CALL DZERO(DINVGP(KZVAR+KZCONF+1),KZWOPT) 2125C 2126C Calculate 2p2h part of DINVGP 2127C 2128 DO 400 I = 1,KZCONF 2129 DINVGP(I) = GP(I) / (DIAG(I) - EFREQ) 2130 DINVGP(KZVAR+I) = GP(KZVAR + I) / (DIAG(I) + EFREQ) 2131 400 CONTINUE 2132 RETURN 2133 END 2134C /* Deck onemp2 */ 2135 SUBROUTINE ONEMP2(PRPMO,DONEPT,ONEACT) 2136C 2137C CALCULATES THE SECOND-ORDER CORRECTION TO AVERAGE VALUE 2138C 2139C Programmed 14/12-1993 by Erik K. Dalskov 2140C 2141#include "implicit.h" 2142C 2143 DIMENSION PRPMO(NORBT,NORBT), DONEPT(NORBT,NORBT) 2144C 2145#include "maxorb.h" 2146#include "inforb.h" 2147#include "infdim.h" 2148#include "infmp2.h" 2149C 2150C 2151 DO 50 IASYM = 1,NSYM 2152 NORBA = NORB(IASYM) 2153 IORBA = IORB(IASYM) 2154 DO 40 I = 1,NORBA 2155 DO 30 J = 1, NORBA 2156 ONEACT = ONEACT + DONEPT(IORBA+J,IORBA+I) * 2157 * PRPMO(IORBA+J,IORBA+I) 215830 CONTINUE 215940 CONTINUE 216050 CONTINUE 2161C 2162C END OF ONEMP2 2163C 2164 RETURN 2165 END 2166C /* Deck prpomp */ 2167 SUBROUTINE PRPOMP(PRPMO,PRPSE,DONEPT) 2168C 2169C Written by Martin Packer 030394. Calculate the second order 2170C corrections to the ph transition moments. 2171C (q|A)_ai = sum_b A_ab D_bi - sum_j D_aj A_ji 2172C + sum_j A_aj D_ji - sum_b D_ab A_bi 2173C (q|A)_ia = sum_j A_ij D_ja - sum_b D_ib P_ba 2174C + sum_b A_ib D_ba - sum_j D_ij A_ja 2175C i,j,k.. are occupied and a,b,c... are virtual orbitals. 2176C 2177C 2178#include "implicit.h" 2179C 2180C 2181#include "maxorb.h" 2182#include "maxash.h" 2183#include "infmp2.h" 2184#include "infpri.h" 2185#include "inforb.h" 2186#include "infrsp.h" 2187#include "infvar.h" 2188#include "wrkrsp.h" 2189C 2190 DIMENSION DONEPT(NORBT,NORBT), PRPMO(NORBT,NORBT) 2191 DIMENSION PRPSE(NORBT,NORBT) 2192C 2193 CALL DZERO(PRPSE,NORBT*NORBT) 2194 DO IBSYM = 1,NSYM 2195 IASYM = MULD2H(IBSYM,KSYMOP) 2196 NORBA = NORB(IASYM) 2197 NORBB = NORB(IBSYM) 2198 NOCCA = NOCC(IASYM) 2199 NSSHA = NSSH(IASYM) 2200 NOCCB = NOCC(IBSYM) 2201 NSSHB = NSSH(IBSYM) 2202 IF ((NORBA.NE.0).AND.(NORBB.NE.0)) THEN 2203C PRPSE(AI) a terms of jo.cpr2 2204 IAST = IORB(IASYM) + 1 + NOCCA 2205 IBST = IORB(IBSYM) + 1 + NOCCB 2206 IBST1= IORB(IBSYM) + 1 2207 CALL DGEMM('N','N',NSSHA,NOCCB,NSSHB,1.D0, 2208 & PRPMO(IAST,IBST),NORBT, 2209 & DONEPT(IBST,IBST1),NORBT,1.D0, 2210 & PRPSE(IAST,IBST1),NORBT) 2211C PRPSE(IA) a terms of jo.cpr2 2212 CALL DGEMM('N','N',NOCCB,NSSHA,NSSHB,-1.D0, 2213 & DONEPT(IBST1,IBST),NORBT, 2214 & PRPMO(IBST,IAST),NORBT,1.D0, 2215 & PRPSE(IBST1,IAST),NORBT) 2216C PRPSE(AI) a terms of jo.cpr2 2217 IAST = IORB(IASYM) + 1 2218 IBST = IORB(IBSYM) + 1 + NOCCB 2219 IBST1= IORB(IBSYM) + 1 2220 CALL DGEMM('N','N',NSSHB,NOCCA,NOCCB,-1.D0, 2221 & DONEPT(IBST,IBST1),NORBT, 2222 & PRPMO(IBST1,IAST),NORBT,1.D0, 2223 & PRPSE(IBST,IAST),NORBT) 2224C PRPSE(IA) a terms of jo.cpr2 2225 CALL DGEMM('N','N',NOCCA,NSSHB,NOCCB,1.D0, 2226 & PRPMO(IAST,IBST1),NORBT, 2227 & DONEPT(IBST1,IBST),NORBT,1.D0, 2228 & PRPSE(IAST,IBST),NORBT) 2229C PRPSE(AI) b terms of jo.cpr2 2230 IAST = IORB(IASYM) + 1 + NOCCA 2231 IBST = IORB(IBSYM) + 1 2232 CALL DGEMM('N','N',NSSHA,NOCCB,NOCCB,1.D0, 2233 & PRPMO(IAST,IBST),NORBT, 2234 & DONEPT(IBST,IBST),NORBT,1.D0, 2235 & PRPSE(IAST,IBST),NORBT) 2236 CALL DGEMM('N','N',NSSHA,NOCCB,NSSHA,-1.D0, 2237 & DONEPT(IAST,IAST),NORBT, 2238 & PRPMO(IAST,IBST),NORBT,1.D0, 2239 & PRPSE(IAST,IBST),NORBT) 2240C PRPSE(IA) b terms of jo.cpr2 2241 CALL DGEMM('N','N',NOCCB,NSSHA,NSSHA,1.D0, 2242 & PRPMO(IBST,IAST),NORBT, 2243 & DONEPT(IAST,IAST),NORBT,1.D0, 2244 & PRPSE(IBST,IAST),NORBT) 2245 CALL DGEMM('N','N',NOCCB,NSSHA,NOCCB,-1.D0, 2246 & DONEPT(IBST,IBST),NORBT, 2247 & PRPMO(IBST,IAST),NORBT,1.D0, 2248 & PRPSE(IBST,IAST),NORBT) 2249 ENDIF 2250 ENDDO 2251 RETURN 2252C 2253C End of PRPOMP 2254C 2255 END 2256C /* Deck prpors */ 2257 SUBROUTINE PRPORS(PRPMO,PRPSE,GP) 2258C 2259C WRITTEN 14-FEB 1986 2260C 2261C PURPOSE: 2262C DISTRIBUTE PROPERTY MO INTEGRALS INTO GP VECTORS IN SOPPA 2263C 2264C List of updates 2265C 21-Jul-1992 Hinne Hettema Average orbital part if necessary. 2266C 030394 - mjp written for SOPPA case 2267C where PRPSE are the second order corrections 2268C 2269#include "implicit.h" 2270C 2271 DIMENSION PRPMO(NORBT,NORBT),PRPSE(NORBT,NORBT),GP(KZYVAR) 2272C 2273C Used from common blocks: 2274C INFORB: NORBT 2275C INFVAR: JWOP() 2276C INFRSP: IPRRSP 2277C WRKRSP: KZYVAR 2278C 2279#include "maxorb.h" 2280#include "priunit.h" 2281#include "inforb.h" 2282#include "infvar.h" 2283#include "infrsp.h" 2284#include "wrkrsp.h" 2285C 2286C -- local constants 2287C 2288 PARAMETER ( D2 = 2.0D0 , DP5 = 0.5D0 ) 2289C 2290C DISTRIBUTE PROPERTY MATRIX IN GP 2291C 2292 ROOT2I = SQRT(DP5) 2293 DO 1400 IG = 1,KZWOPT 2294 K = JWOP(1,IG) 2295 L = JWOP(2,IG) 2296 GP(KZCONF+IG) = GP(KZCONF+IG) 2297 * + D2 * PRPMO(L,K) + PRPSE(L,K) 2298 GP(KZVAR+KZCONF+IG) = GP(KZVAR+KZCONF+IG) 2299 * - D2 * PRPMO(K,L) + PRPSE(K,L) 2300 1400 CONTINUE 2301C 2302C CALL DSCAL(KZWOPT,0.0D0,GP(KZCONF+1),1) 2303C CALL DSCAL(KZWOPT,0.0D0,GP(KZVAR+KZCONF+1),1) 2304C *** Perform averaging 2305C 2306 IF (RSPSUP .AND. (KSYMOP .EQ. 1)) THEN 2307 CALL RSPAVE(GP(KZCONF+1),KZVAR,2) 2308 END IF 2309C 2310 IF (IPRRSP.GT.20) THEN 2311 WRITE(LUPRI,'(/A)') 2312 & ' (PRPORS) Z and Y property GP vector in SOPPA approximation' 2313 CALL OUTPUT(GP,1,KZVAR,1,2,KZVAR,2,1,LUPRI) 2314 END IF 2315 RETURN 2316 END 2317C /* Deck prpcmp */ 2318 SUBROUTINE PRPCMP(PRPMO,COEMP2,GP,ABSADR,ABTADR,IJSADR,IJTADR, 2319 & IADR1,WRK,KFREE,LFREE) 2320C 2321C Copyright 8/4-94 by Erik K. Dalskov 2322C 2323C Purpose: Construct 2p-2h corrections to transition moments 2324C in SOPPA 2325C 2326#include "implicit.h" 2327C 2328C Used from common blocks: 2329C INFRSP : IPRRSP 2330C WRKRSP : KSYMOP,KZYVAR 2331C 2332#include "maxorb.h" 2333#include "maxash.h" 2334#include "priunit.h" 2335#include "infpri.h" 2336#include "inforb.h" 2337#include "infsop.h" 2338#include "infmp2.h" 2339 DIMENSION IADR1(NP,NH) 2340#include "infind.h" 2341#include "infrsp.h" 2342#include "wrkrsp.h" 2343C 2344 INTEGER BSYM,CSYM,A,AA,B,BB,ABSADR,ABTADR 2345 PARAMETER ( D1 = 1.0D0 , D3 = 3.0D0 , D4 = 4.0D0 ) 2346 PARAMETER ( D2 = 2.0D0 , D8 = 8.0D0 ) 2347C 2348 DIMENSION PRPMO(NORBT,*), COEMP2(*), GP(KZYVAR) 2349 DIMENSION ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 2350 & IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM) 2351 DIMENSION WRK(*) 2352C 2353 ROOT3I = D1 / SQRT(D3) 2354 ROOT2 = SQRT(D2) 2355 ROOT8 = SQRT(D8) 2356C 2357C Work space allocation 2358C 2359 KFRSAV = KFREE 2360 CALL MEMGET('REAL',KMPUNP,NORBT*NORBT,WRK,KFREE,LFREE) 2361 CALL MEMGET('REAL',KRESZ,NP*NH,WRK,KFREE,LFREE) 2362 CALL MEMGET('REAL',KRESY,NP*NH,WRK,KFREE,LFREE) 2363 LMPUNP = KMPUNP - 1 2364 LRESZ = KRESZ - 1 2365 LRESY = KRESY - 1 2366C 2367 DO II=1,NH 2368 I = IPHORD(II) 2369 IF (IFRMP2(I) .NE. 0) GO TO 9999 2370 DO AA=1,NP 2371 A = IPHORD(NH+AA) 2372 IASYM = MULD2H(ISMO(A),ISMO(I)) 2373 JBSYM = MULD2H(IASYM,KSYMOP) 2374 CALL PHPACK(WRK(KMPUNP),COEMP2(IADR1(AA,II)+1),IASYM,-1,-1) 2375 DO JSYM=1,NSYM 2376 BSYM = MULD2H(JSYM,JBSYM) 2377 CSYM = MULD2H(JSYM,IASYM) 2378 KSYM = MULD2H(BSYM,IASYM) 2379 NOCCJ = NOCC(JSYM) - NFRMP2(JSYM) 2380 NOCCK = NOCC(KSYM) 2381 NOCCB = NOCC(BSYM) 2382 NOCCC = NOCC(CSYM) 2383 NSSHB = NSSH(BSYM) 2384 NSSHC = NSSH(CSYM) 2385 IF ((NOCCJ .NE. 0) .AND. (NSSHB .NE. 0)) THEN 2386 IBST = IORB(BSYM) + NOCCB + 1 2387 ICST = IORB(CSYM) + NOCCC + 1 2388 IJST = IORB(JSYM) + NFRMP2(JSYM) + 1 2389 IKST = IORB(KSYM) + 1 2390 JCST = ( IORB(JSYM) + NFRMP2(JSYM) ) * NORBT + ICST 2391 KBST = IORB(KSYM) * NORBT + IBST 2392 IF (NSSHC .EQ. 0) THEN 2393 CALL DZERO(WRK(KRESZ),NSSHB*NOCCJ) 2394 ELSE 2395 CALL DGEMM('N','N',NSSHB,NOCCJ,NSSHC,1.D0, 2396 & PRPMO(IBST,ICST),NORBT, 2397 & WRK(LMPUNP+JCST),NORBT,0.D0, 2398 & WRK(KRESZ),NSSHB) 2399 END IF 2400 CALL DGEMM('N','N',NSSHB,NOCCJ,NOCCK,-1.D0, 2401 & WRK(LMPUNP+KBST),NORBT, 2402 & PRPMO(IKST,IJST),NORBT,1.D0, 2403 & WRK(KRESZ),NSSHB) 2404 IF (NOCCK .EQ. 0) THEN 2405 CALL DZERO(WRK(KRESY),NSSHB*NOCCJ) 2406 ELSE 2407 CALL DGEMM('N','T',NSSHB,NOCCJ,NOCCK,1.D0, 2408 & WRK(LMPUNP+KBST),NORBT, 2409 & PRPMO(IJST,IKST),NORBT,0.D0, 2410 & WRK(KRESY),NSSHB) 2411 END IF 2412 CALL DGEMM('T','N',NSSHB,NOCCJ,NSSHC,-1.D0, 2413 & PRPMO(ICST,IBST),NORBT, 2414 & WRK(LMPUNP+JCST),NORBT,1.D0, 2415 & WRK(KRESY),NSSHB) 2416 DO J=1,NOCCJ 2417 JJ = INDXPH(IJST-1+J) 2418 IJSOFF = IJSADR(II,JJ,KSYMOP) 2419 IJTOFF = IJTADR(II,JJ,KSYMOP) 2420 IF (II .GE. JJ) THEN 2421 TIJSGN = ROOT3I 2422 ELSE 2423 TIJSGN = -ROOT3I 2424 ENDIF 2425 DO B = 1,NSSHB 2426 BB = - INDXPH(IBST-1+B) 2427 NABS = ABSADR(AA,BB) 2428 NABT = ABTADR(AA,BB) 2429 IF (II .EQ. JJ) THEN 2430 IF (AA .EQ. BB) THEN 2431 FAC = D2 2432 ELSE 2433 FAC = ROOT2 2434 ENDIF 2435 ELSE 2436 IF (AA .EQ. BB) THEN 2437 FAC = ROOT2 2438 ELSE 2439 FAC = D1 2440 ENDIF 2441 ENDIF 2442 GP(IJSOFF + NABS) = 2443 * GP(IJSOFF + NABS) 2444 * - WRK(LRESZ+(J-1)*NSSHB+B) * FAC 2445 GP(KZVAR + IJSOFF + NABS) = 2446 * GP(KZVAR + IJSOFF + NABS) 2447 * - WRK(LRESY+(J-1)*NSSHB+B) * FAC 2448 IF ((II .NE. JJ) .AND. (AA .NE. BB)) THEN 2449 IF (AA .GT .BB) THEN 2450 TSIGN = -TIJSGN 2451 ELSE 2452 TSIGN = TIJSGN 2453 ENDIF 2454 GP(IJTOFF + NABT) = GP(IJTOFF + NABT) 2455 * + TSIGN * WRK(LRESZ+(J-1)*NSSHB+B) 2456 GP(KZVAR + IJTOFF + NABT) = 2457 * GP(KZVAR + IJTOFF + NABT) 2458 * + TSIGN * WRK(LRESY+(J-1)*NSSHB+B) 2459 ENDIF 2460 ENDDO 2461 ENDDO 2462 ENDIF 2463 ENDDO 2464 ENDDO 2465 9999 CONTINUE 2466 ENDDO 2467 IF ( IPRRSP.GT.110) THEN 2468 WRITE(LUPRI,'(/A)') 2469 * ' 2p-2h transition moments for SOPPA (Z-part and Y-part)' 2470 CALL OUTPUT(GP,1,KZCONF,1,2,KZVAR,2,1,LUPRI) 2471 END IF 2472C 2473C End of PRPCMP 2474C 2475 CALL MEMREL('PRPCMP',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 2476 RETURN 2477 END 2478C /* Deck prpcmpt */ 2479 SUBROUTINE PRPTCMP(PRPMO,COEMP2,GP,ABSADR,ABTADR,IJ1ADR,IJ2ADR, 2480 & IJ3ADR,IADR1,WRK,KFREE,LFREE) 2481C 2482C Written 14-Feb-1995 by Thomas Enevoldsen (tec) 2483C 2484C Purpose: Construct 2p-2h contributions to transition moments 2485C in triplet SOPPA 2486C 2487C 2488C Notice COEMP2 = Kappa / 2 2489C 2490C 2491#include "implicit.h" 2492C 2493C Used from common blocks: 2494C INFRSP : IPRRSP 2495C WRKRSP : KSYMOP,KZYVAR 2496C 2497#include "maxorb.h" 2498#include "maxash.h" 2499#include "priunit.h" 2500#include "infpri.h" 2501#include "inforb.h" 2502#include "infsop.h" 2503#include "infmp2.h" 2504 DIMENSION IADR1(NP,NH) 2505#include "infind.h" 2506#include "infrsp.h" 2507#include "wrkrsp.h" 2508C 2509 INTEGER BSYM,CSYM,A,AA,B,BB,ABSADR,ABTADR 2510 PARAMETER ( D1 = 1.0D0 , D2 = 2.0D0 , D3 = 3.0D0 , D6 = 6.0D0 ) 2511C 2512 DIMENSION PRPMO(NORBT,*), COEMP2(*), GP(KZYVAR) 2513 DIMENSION ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 2514 & IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM), 2515 & IJ3ADR(NISHT,NISHT,NSYM) 2516 DIMENSION WRK(*) 2517C 2518 ROOT2 = SQRT(D2) 2519 D3I = D1 / D3 2520 T1FAC = ROOT2 * D3I 2521C 2522C Work space allocation 2523C 2524 KFRSAV = KFREE 2525 CALL MEMGET('REAL',KMPUNP,NORBT*NORBT,WRK,KFREE,LFREE) 2526 CALL MEMGET('REAL',KRESZ1,NP*NH,WRK,KFREE,LFREE) 2527 CALL MEMGET('REAL',KRESZ2,NP*NH,WRK,KFREE,LFREE) 2528 CALL MEMGET('REAL',KRESY1,NP*NH,WRK,KFREE,LFREE) 2529 CALL MEMGET('REAL',KRESY2,NP*NH,WRK,KFREE,LFREE) 2530 LMPUNP = KMPUNP - 1 2531 LRESZ1 = KRESZ1 - 1 2532 LRESY1 = KRESY1 - 1 2533 LRESZ2 = KRESZ2 - 1 2534 LRESY2 = KRESY2 - 1 2535C 2536 DO II=1,NH 2537 I = IPHORD(II) 2538 IF (IFRMP2(I) .NE. 0) GO TO 9999 2539 DO AA=1,NP 2540 A = IPHORD(NH+AA) 2541 IASYM = MULD2H(ISMO(A),ISMO(I)) 2542 JBSYM = MULD2H(IASYM,KSYMOP) 2543 CALL PHPACK(WRK(KMPUNP),COEMP2(IADR1(AA,II)+1),IASYM,-1,-1) 2544 DO JSYM=1,NSYM 2545 BSYM = MULD2H(JSYM,JBSYM) 2546 CSYM = MULD2H(JSYM,IASYM) 2547 KSYM = MULD2H(BSYM,IASYM) 2548 NOCCJ = NOCC(JSYM) - NFRMP2(JSYM) 2549 NOCCK = NOCC(KSYM) 2550 NOCCB = NOCC(BSYM) 2551 NOCCC = NOCC(CSYM) 2552 NSSHB = NSSH(BSYM) 2553 NSSHC = NSSH(CSYM) 2554 IF ((NOCCJ .NE. 0) .AND. (NSSHB .NE. 0)) THEN 2555 IBST = IORB(BSYM) + NOCCB + 1 2556 ICST = IORB(CSYM) + NOCCC + 1 2557 IJST = IORB(JSYM) + NFRMP2(JSYM) + 1 2558 IKST = IORB(KSYM) + 1 2559 JCST = ( IORB(JSYM) + NFRMP2(JSYM) ) * NORBT + ICST 2560 KBST = IORB(KSYM) * NORBT + IBST 2561 IF (NSSHC .EQ. 0) THEN 2562 CALL DZERO(WRK(KRESZ1),NSSHB*NOCCJ) 2563 CALL DZERO(WRK(KRESY2),NSSHB*NOCCJ) 2564 ELSE 2565 CALL DGEMM('N','N',NSSHB,NOCCJ,NSSHC,1.D0, 2566 & PRPMO(IBST,ICST),NORBT, 2567 & WRK(LMPUNP+JCST),NORBT,0.D0, 2568 & WRK(KRESZ1),NSSHB) 2569 CALL DGEMM('T','N',NSSHB,NOCCJ,NSSHC,1.D0, 2570 & PRPMO(ICST,IBST),NORBT, 2571 & WRK(LMPUNP+JCST),NORBT,0.D0, 2572 & WRK(KRESY2),NSSHB) 2573 END IF 2574 IF (NOCCK .EQ. 0) THEN 2575 CALL DZERO(WRK(KRESZ2),NSSHB*NOCCJ) 2576 CALL DZERO(WRK(KRESY1),NSSHB*NOCCJ) 2577 ELSE 2578 CALL DGEMM('N','N',NSSHB,NOCCJ,NOCCK,1.D0, 2579 & WRK(LMPUNP+KBST),NORBT, 2580 & PRPMO(IKST,IJST),NORBT,0.D0, 2581 & WRK(KRESZ2),NSSHB) 2582 CALL DGEMM('N','T',NSSHB,NOCCJ,NOCCK,1.D0, 2583 & WRK(LMPUNP+KBST),NORBT, 2584 & PRPMO(IJST,IKST),NORBT,0.D0, 2585 & WRK(KRESY1),NSSHB) 2586 ENDIF 2587 DO J=1,NOCCJ 2588 JJ = INDXPH(IJST-1+J) 2589 IJ1OFF = IJ1ADR(II,JJ,KSYMOP) 2590 IJ2OFF = IJ2ADR(II,JJ,KSYMOP) 2591 IJ3OFF = IJ3ADR(II,JJ,KSYMOP) 2592 IF (II .GT. JJ) THEN 2593 TIJSGN = D1 2594 ELSE 2595 TIJSGN = -D1 2596 ENDIF 2597 IF (II .EQ. JJ) THEN 2598 T3FAC = ROOT2 2599 ELSE 2600 T3FAC = D1 2601 ENDIF 2602 DO B = 1,NSSHB 2603 BB = - INDXPH(IBST-1+B) 2604 NABS = ABSADR(AA,BB) 2605 NABT = ABTADR(AA,BB) 2606 IF (AA .GT. BB) THEN 2607 TABSGN = D1 2608 ELSE 2609 TABSGN = -D1 2610 ENDIF 2611 IF (AA .EQ. BB) THEN 2612 T2FAC = ROOT2 2613 ELSE 2614 T2FAC = D1 2615 ENDIF 2616C T(1) part of GP-vector 2617 IF ((II .NE. JJ) .AND. (AA .NE. BB)) THEN 2618 GP(IJ1OFF + NABT) = GP(IJ1OFF + NABT) 2619 * - TABSGN * TIJSGN * T1FAC * 2620 * (WRK(LRESZ1+(J-1)*NSSHB+B)-WRK(LRESZ2+(J-1) 2621 * *NSSHB+B)) 2622 GP(KZVAR + IJ1OFF + NABT) = 2623 * GP(KZVAR + IJ1OFF + NABT) 2624 * - TABSGN * TIJSGN * T1FAC * 2625 * (WRK(LRESY1+(J-1)*NSSHB+B)-WRK(LRESY2+(J-1) 2626 * *NSSHB+B)) 2627 ENDIF 2628C T(2) part of GP-vector 2629 IF (II .NE. JJ) THEN 2630 GP(IJ2OFF + NABS) = GP(IJ2OFF + NABS) 2631 * + TIJSGN * T2FAC * 2632 * (WRK(LRESZ2+(J-1)*NSSHB+B) - WRK(LRESZ1 2633 * +(J-1)*NSSHB+B) * D3I) 2634 GP(KZVAR + IJ2OFF + NABS) = GP(KZVAR + IJ2OFF 2635 * + NABS) + TIJSGN * T2FAC * 2636 * (WRK(LRESY2+(J-1)*NSSHB+B) * D3I - WRK(LRESY1 2637 * +(J-1)*NSSHB+B) ) 2638 ENDIF 2639C T(3) part of GP-vector 2640 IF (AA .NE. BB) THEN 2641 GP(IJ3OFF + NABT) = GP(IJ3OFF + NABT) 2642 * + TABSGN * T3FAC * 2643 * (WRK(LRESZ2+(J-1)*NSSHB+B) * D3I - WRK(LRESZ1 2644 * +(J-1)*NSSHB+B) ) 2645 GP(KZVAR + IJ3OFF + NABT) = GP(KZVAR + IJ3OFF 2646 * + NABT) + TABSGN * T3FAC * 2647 * (WRK(LRESY2+(J-1)*NSSHB+B) - WRK(LRESY1 2648 * +(J-1)*NSSHB+B) * D3I ) 2649 ENDIF 2650 ENDDO 2651 ENDDO 2652 ENDIF 2653 ENDDO 2654 ENDDO 2655 9999 CONTINUE 2656 ENDDO 2657 IF ( IPRRSP.GT.110) THEN 2658 WRITE(LUPRI,'(/A)') 2659 * ' 2p-2h transition moments for SOPPA (Z-part and Y-part)' 2660 CALL OUTPUT(GP(1),1,KZCONF,1,2,KZVAR,2,1,LUPRI) 2661 END IF 2662C 2663C End of PRPTCMP 2664C 2665 CALL MEMREL('PRPTCMP',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 2666 RETURN 2667 END 2668C /* Deck rspact */ 2669 SUBROUTINE RSPACT 2670C 2671C 16-Nov-1993 mjp+hjaaj 2672C 2673#include "implicit.h" 2674C 2675C Used from common blocks: 2676C INFRSP: SOPPA,LPVMAT,NACTT,NACT(8),IACT(8),JACT(8) 2677C INFORB: NSYM,*ASH*,*ORB*,NVIRT,NRHFT 2678C INFMP2: NPHTOT(1) 2679C 2680#include "maxorb.h" 2681#include "infrsp.h" 2682#include "inforb.h" 2683#include "infmp2.h" 2684C 2685 IF (SOPPA) THEN 2686 LPVMAT = NPHTOT(1) 2687 NACTT = NORBT 2688 DO ISYM = 1,NSYM 2689 NACT(ISYM) = NORB(ISYM) 2690 IACT(ISYM) = IORB(ISYM) 2691 JACT(ISYM) = IORB(ISYM) 2692 END DO 2693 ELSE 2694 LPVMAT = N2ASHX*N2ASHX 2695 NACTT = NASHT 2696 DO ISYM = 1,NSYM 2697 NACT(ISYM) = NASH(ISYM) 2698 IACT(ISYM) = IASH(ISYM) 2699 JACT(ISYM) = IORB(ISYM) + NISH(ISYM) 2700 END DO 2701 END IF 2702 RETURN 2703 END 2704C /* Deck sopij */ 2705 SUBROUTINE SOPIJ(ISYMV,IOFFY,IUNIT,RSPVEC,ABSADR,ABTADR,IJSADR, 2706 & IJTADR,IJ1ADR,IJ2ADR,IJ3ADR,THPSOP) 2707C 2708C This routine sums up all excitation contributions to a specific 2709C I,J pair in the 2p2h solution vector. 2710C 2711C Copyright 23/2-96 by Erik K. Dalskov and Thomas Enevoldsen 2712C 2713#include "implicit.h" 2714C 2715C 2716#include "maxorb.h" 2717#include "maxash.h" 2718#include "infpri.h" 2719#include "inforb.h" 2720#include "infsop.h" 2721#include "infmp2.h" 2722#include "infind.h" 2723#include "infrsp.h" 2724C 2725 INTEGER ABSYM,BASYM,ASYM,BSYM,A,B,ABSOFF,ABTOFF,ABSADR,ABTADR 2726 DIMENSION RSPVEC(*), ABSADR(NORBT,NORBT), ABTADR(NORBT,NORBT), 2727 & IJSADR(NISHT,NISHT,NSYM), IJTADR(NISHT,NISHT,NSYM), 2728 & IJ1ADR(NISHT,NISHT,NSYM), IJ2ADR(NISHT,NISHT,NSYM), 2729 & IJ3ADR(NISHT,NISHT,NSYM) 2730 PARAMETER ( ZERO = 0.0D0 ) 2731C 2732 IKKEPR = 0 2733 EXCITOT = ZERO 2734 THPR = THPSOP * THPSOP 2735 IF (.NOT.TRPLET) THEN 2736 WRITE(IUNIT,1000) 2737 DO I=1,NH 2738 II = IPHORD(I) 2739 IF (IFRMP2(II) .EQ. 0) THEN 2740 ISYM=ISMO(II) 2741 DO J=1,I 2742 JJ = IPHORD(J) 2743 IF (IFRMP2(JJ).EQ.0) THEN 2744 JSYM=ISMO(JJ) 2745 IJSYM=MULD2H(ISYM,JSYM) 2746 ABSYM=MULD2H(ISYMV,IJSYM) 2747 IJSOFF=IJSADR(I,J,ISYMV) 2748 IJTOFF=IJTADR(I,J,ISYMV) 2749 EXCISZ = ZERO 2750 EXCITZ = ZERO 2751 EXCISY = ZERO 2752 EXCITY = ZERO 2753 DO A=1,NP 2754 ASYM=ISMO(IPHORD(NH+A)) 2755 DO B=1,A 2756 BSYM=ISMO(IPHORD(NH+B)) 2757 BASYM=MULD2H(BSYM,ASYM) 2758 IF (ABSYM .EQ. BASYM) THEN 2759 ABSOFF = ABSADR(A,B) 2760 ABTOFF = ABTADR(A,B) 2761 EXCISZ = EXCISZ + RSPVEC(IJSOFF+ABSOFF) 2762 & * RSPVEC(IJSOFF+ABSOFF) 2763 EXCISY = EXCISY 2764 & + RSPVEC(IOFFY+IJSOFF+ABSOFF) 2765 & * RSPVEC(IOFFY+IJSOFF+ABSOFF) 2766 IF (I.NE.J .AND. A.NE.B) THEN 2767 EXCITZ = EXCITZ + RSPVEC(IJTOFF+ABTOFF) 2768 & * RSPVEC(IJTOFF+ABTOFF) 2769 EXCITY = EXCITY 2770 & + RSPVEC(IOFFY+IJTOFF+ABTOFF) 2771 & * RSPVEC(IOFFY+IJTOFF+ABTOFF) 2772 ENDIF 2773 END IF 2774 ENDDO 2775 ENDDO 2776 EXCIMAX = MAX(EXCISZ,EXCISY,EXCITZ,EXCITY) 2777 IF (EXCIMAX .GT. THPR) THEN 2778 IKKEPR = IKKEPR + 1 2779 EXCISZ = SQRT(EXCISZ) 2780 EXCISY = SQRT(EXCISY) 2781 EXCITZ = SQRT(EXCITZ) 2782 EXCITY = SQRT(EXCITY) 2783 WRITE(IUNIT,1001) II,ISYM,JJ,JSYM, 2784 & EXCISZ,EXCISY,EXCITZ,EXCITY 2785 ELSE 2786 EXCITOT = EXCITOT + 2787 & EXCISZ + EXCISY + EXCITZ + EXCITY 2788 ENDIF 2789 ENDIF 2790 ENDDO 2791 ENDIF 2792 ENDDO 2793 ELSE 2794 WRITE(IUNIT,1100) 2795 DO I=1,NH 2796 II = IPHORD(I) 2797 IF (IFRMP2(II) .EQ. 0) THEN 2798 ISYM=ISMO(II) 2799 DO J=1,I 2800 JJ = IPHORD(J) 2801 IF (IFRMP2(JJ) .EQ. 0) THEN 2802 JSYM=ISMO(JJ) 2803 IJSYM=MULD2H(ISYM,JSYM) 2804 ABSYM=MULD2H(ISYMV,IJSYM) 2805 IJ1OFF=IJ1ADR(I,J,ISYMV) 2806 IJ2OFF=IJ2ADR(I,J,ISYMV) 2807 IJ3OFF=IJ3ADR(I,J,ISYMV) 2808 EXCI1Z = ZERO 2809 EXCI2Z = ZERO 2810 EXCI3Z = ZERO 2811 EXCI1Y = ZERO 2812 EXCI2Y = ZERO 2813 EXCI3Y = ZERO 2814 DO A=1,NP 2815 ASYM=ISMO(IPHORD(NH+A)) 2816 DO B=1,A 2817 BSYM=ISMO(IPHORD(NH+B)) 2818 BASYM=MULD2H(BSYM,ASYM) 2819 IF (ABSYM .EQ. BASYM) THEN 2820 ABSOFF = ABSADR(A,B) 2821 ABTOFF = ABTADR(A,B) 2822 IF (I.NE.J) THEN 2823 IF(A.NE.B) THEN 2824 EXCI1Z = EXCI1Z + 2825 & RSPVEC(IJ1OFF+ABTOFF) * 2826 & RSPVEC(IJ1OFF+ABTOFF) 2827 EXCI2Z = EXCI2Z + 2828 & RSPVEC(IJ2OFF+ABSOFF) * 2829 & RSPVEC(IJ2OFF+ABSOFF) 2830 EXCI3Z = EXCI3Z + 2831 & RSPVEC(IJ3OFF+ABTOFF) * 2832 & RSPVEC(IJ3OFF+ABTOFF) 2833 EXCI1Y = EXCI1Y + 2834 & RSPVEC(IOFFY+IJ1OFF+ABTOFF) * 2835 & RSPVEC(IOFFY+IJ1OFF+ABTOFF) 2836 EXCI2Y = EXCI2Y + 2837 & RSPVEC(IOFFY+IJ2OFF+ABSOFF) * 2838 & RSPVEC(IOFFY+IJ2OFF+ABSOFF) 2839 EXCI3Y = EXCI3Y + 2840 & RSPVEC(IOFFY+IJ3OFF+ABTOFF) * 2841 & RSPVEC(IOFFY+IJ3OFF+ABTOFF) 2842 ELSE 2843 EXCI2Z = EXCI2Z + 2844 & RSPVEC(IJ2OFF+ABSOFF) * 2845 & RSPVEC(IJ2OFF+ABSOFF) 2846 EXCI2Y = EXCI2Y + 2847 & RSPVEC(IOFFY+IJ2OFF+ABSOFF) * 2848 & RSPVEC(IOFFY+IJ2OFF+ABSOFF) 2849 ENDIF 2850 ELSE 2851 IF (A.NE.B) THEN 2852 EXCI3Z = EXCI3Z + 2853 & RSPVEC(IJ3OFF+ABTOFF) * 2854 & RSPVEC(IJ3OFF+ABTOFF) 2855 EXCI3Y = EXCI3Y + 2856 & RSPVEC(IOFFY+IJ3OFF+ABTOFF) * 2857 & RSPVEC(IOFFY+IJ3OFF+ABTOFF) 2858 ENDIF 2859 ENDIF 2860 END IF 2861 ENDDO 2862 ENDDO 2863 EXCIMAX = MAX(EXCI1Z,EXCI1Y,EXCI2Z, 2864 & EXCI2Y,EXCI3Z,EXCI3Y) 2865 IF (EXCIMAX .GT. THPR) THEN 2866 IKKEPR = IKKEPR + 1 2867 EXCI1Z = SQRT(EXCI1Z) 2868 EXCI2Z = SQRT(EXCI2Z) 2869 EXCI3Z = SQRT(EXCI3Z) 2870 EXCI1Y = SQRT(EXCI1Y) 2871 EXCI2Y = SQRT(EXCI2Y) 2872 EXCI3Y = SQRT(EXCI3Y) 2873 WRITE(IUNIT,1101) II,ISYM,JJ,JSYM, 2874 & EXCI1Z,EXCI1Y,EXCI2Z,EXCI2Y,EXCI3Z,EXCI3Y 2875 ELSE 2876 EXCITOT = EXCITOT + EXCI1Z + EXCI2Z + EXCI3Z 2877 & + EXCI1Y + EXCI2Y + EXCI3Y 2878 ENDIF 2879 ENDIF 2880 ENDDO 2881 ENDIF 2882 ENDDO 2883 ENDIF 2884 EXCITOT = SQRT(EXCITOT) 2885 IKKEPR = NH*(NH+1)/2 - IKKEPR 2886 WRITE(IUNIT,1500) IKKEPR,THPSOP,EXCITOT 2887C 2888C End of SOPIJ 2889C 2890 RETURN 2891 1000 FORMAT(/' i j S Z oper. Y oper.', 2892 & ' T Z oper. Y oper.', 2893 & /' --- --- --- -------- --------', 2894 & ' --- -------- --------') 2895 1001 FORMAT(I4,'(',I1,')',I4,'(',I1,')',3X,2F10.6,3X,2F10.6) 2896 1100 FORMAT(/' i j 1 Z oper. Y oper.', 2897 & ' 2 Z oper. Y oper.', 2898 & ' 3 Z oper. Y oper.', 2899 & /' --- --- --- ------- -------', 2900 & ' --- ------- -------', 2901 & ' --- ------- -------') 2902 1101 FORMAT(I4,'(',I1,')',I4,'(',I1,')',3X,2F9.5,3X,2F9.5,3X,2F9.5) 2903 1500 FORMAT(/5X,I5,' elements with norm less than',1P,D9.2, 2904 & ' not printed', 2905 & //' The total norm of the excitation out of these', 2906 & ' elements is',1P,D11.4) 2907 END 2908! -- end of rspsoppa.F -- 2909