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#define DFTQRC dft_qr_respons 19c#define DFTQRC dftqrcf 20C 21#ifdef REVLOG 22C========================================================================= 23CRevision 1.2 2000/05/24 19:06:44 hjj 24Cnew getref calls with appropriate NDREF instead of NCREF 25C(fixing error for triplet with CSF) 26C941223-hjaaj: new QFOCK call 27C931015-hjaaj: GTZYMT, TRZYMT: stop if NSIM .ne. 1 28C=========================================================================== 29#endif 30C /* Deck t3drv */ 31 SUBROUTINE T3DRV(NSIM,ISYMA,ISYMB,ISYMC,VECB,VECC,ATEST,VECA, 32 * OMEGAB,OMEGAC,XINDX,UDV,PV,MJWOP,WRK,LWRK, 33 * CMO,FC,FV) 34C 35C Purpose: 36C Compute E[3] times two vectors 37C Compute S[3] times two vectors, multiply with omega, and add to result 38C 39C OUTPUT is put in WRK: 40C First NSIM * KZYVA elements are result vector 41C 42#include "implicit.h" 43C 44 PARAMETER ( D0 = 0.0D0 , DM1 = -1.0D0 , DEM10 = 1.0D-10 ) 45C 46#include "priunit.h" 47#include "thrzer.h" 48#include "infvar.h" 49#include "qrinf.h" 50#include "rspprp.h" 51#include "infhyp.h" 52#include "infrsp.h" 53#include "wrkrsp.h" 54#include "infpri.h" 55#include "infspi.h" 56#include "gnrinf.h" 57C 58 LOGICAL ATEST, TEST 59 DATA TEST /.FALSE./ 60C 61 DIMENSION WRK(*) 62 DIMENSION XINDX(*) 63 DIMENSION UDV(*) 64 DIMENSION PV(*), MJWOP(2,MAXWOP,8) 65 DIMENSION VECB(*), VECC(*), VECA(*) 66 DIMENSION CMO(*),FC(*) 67C 68 CALL QENTER('T3DRV') 69C 70C Compute the length of storage needed 71C KZYVA is the length of the resulting vector 72C 73 KZYVA = MZYVAR(ISYMA) 74 KZYVB = MZYVAR(ISYMB) 75 KZYVC = MZYVAR(ISYMC) 76C 77 IBEQC = 0 78 IF ( HYPCAL .AND. ISYMB .EQ. ISYMC .AND. ISPINB.EQ.ISPINC .AND. 79 & ABS(OMEGAB) .EQ. ABS(OMEGAC) .AND. 80 & .NOT. (RSPCI .OR. TDA .OR. CISRPA) ) THEN 81 IF (OMEGAB .EQ. OMEGAC) THEN 82 BMINC = D0 83 DO 110 I = 1,KZYVB 84 BMINC = BMINC + ABS(VECB(I) - VECC(I)) 85 110 CONTINUE 86 IF (ATEST) WRITE(LUPRI,*) '1norm(B - C) = BMINC =',BMINC 87 IF (BMINC .LE. THRZER) IBEQC = 1 88 ELSE 89 BMINC = D0 90 KZVB = MZVAR(ISYMB) 91 DO 120 I = 1,KZVB 92 BMINC = MAX(BMINC,ABS(VECB(I) + VECC(KZVB+I) 93 * + VECC(I) + VECB(KZVB+I)) ) 94 120 CONTINUE 95 IF (ATEST) WRITE(LUPRI,*) 'max(B + Ct) = BMINC =',BMINC 96 IF (BMINC .LE. DEM10) IBEQC = -1 97C MAERKE: note that this will not catch dipole velocity type 98C where VECB(I) = VECC(KZVB+I) /910905-pj+hjaaj 99 END IF 100 END IF 101C 102C Assumed order of spin operators in triplet quadratic response 103C 104C 105 IF( IPRRSP .GT. 50) THEN 106 WRITE(LUPRI, '(//A)') 107 * ' Characteristics in T3DRV routine' 108 WRITE(LUPRI, '( //A,2I10 )') 109 * ' Symmetry of vectors b and c : ', ISYMB,ISYMC 110 WRITE(LUPRI, '( A,2I10 )') 111 * ' Spin of vectors b and c : ', ISPINB,ISPINC 112 WRITE(LUPRI, '( A,2I10 )') 113 * ' Length of orbital part : ', 114 * MZWOPT(ISYMB),MZWOPT(ISYMC) 115 WRITE(LUPRI, '( A,2I10 )') 116 * ' Length of configuration part : ', MZCONF(ISYMB),MZCONF(ISYMC) 117 WRITE(LUPRI, '( A,2I10 )') 118 * ' Length of vectors : ', KZYVB,KZYVC 119 WRITE(LUPRI, '( A,2F10.7 )') 120 * ' Frequencies of vectors : ', OMEGAB,OMEGAC 121 WRITE(LUPRI, '( A,I10 )') 122 * ' vector b equal to vector c ? : ', IBEQC 123 END IF 124C 125 KE3 = 1 126 KS3 = KE3 + KZYVA 127 KWRK = KS3 + KZYVA 128 LWRKE = LWRK - KS3 129 LWRKS = LWRK - KWRK 130 IF (LWRKS.LT.0) CALL ERRWRK('T3DRV',KWRK-1,LWRK) 131C 132C (Keep WRK as the true workspace and split off the computed arrays) 133C 134 IF( (NSIM .GT. 1) ) THEN 135 WRITE(LUERR,'(//A,/A,I5,/A)') 136 * ' --> Error : Routine T3DRV not implemented with NSIM > 1', 137 * ' NSIM = ', NSIM, 138 * ' Halting execution.' 139 CALL QTRACE(LUERR) 140 CALL QUIT(' T3DRV ERROR: NSIM GREATER THAN 1 IN T3DRV') 141 END IF 142C 143C Initialise the result vector 144C 145 CALL DZERO(WRK(KE3),KZYVA) 146C 147C Return if we have a CI response calculation 148C 149 IF ( RSPCI .OR. TDA .OR. CISRPA) RETURN 150C 151 IF ( TEST ) THEN 152 CALL DCOPY(MZVAR(ISYMB),VECB,1,VECB(1+MZVAR(ISYMB)),1) 153 CALL DSCAL(MZVAR(ISYMB),DM1,VECB(1+MZVAR(ISYMB)),1) 154 CALL DCOPY(MZVAR(ISYMC),VECC,1,VECC(1+MZVAR(ISYMC)),1) 155 CALL DSCAL(MZVAR(ISYMC),DM1,VECC(1+MZVAR(ISYMC)),1) 156 END IF 157C 158C Compute ( E3(jkl) + E3(jlk) ) N(l) N(k) 159C 160 CALL E3INIT(VECB,VECC,VECA,ATEST,IBEQC, 161 * WRK(KE3),XINDX,UDV,PV,WRK(KS3),LWRKE, 162 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC, 163 * ISPINA,ISPINB,ISPINC,CMO,FC,FV,MJWOP) 164C 165C Compute S3(jlk) N(l) N(k) 166C 167 IF (OMEGAB .NE. D0 ) THEN 168 CALL S3INIT(KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC, 169 * ISPINA,ISPINB,ISPINC, 170 * VECB,VECC,VECA,ATEST, 171 * WRK(KS3),XINDX,UDV,MJWOP,WRK(KWRK),LWRKS) 172C 173C Multiply with omega-b 174C 175 CALL DAXPY(KZYVA,OMEGAB,WRK(KS3),1,WRK(KE3),1) 176 END IF 177C 178C Compute S3(jkl) N(l) N(k) 179C 180 IF (OMEGAC .NE. D0 ) THEN 181 CALL S3INIT(KZYVA,KZYVC,KZYVB,ISYMA,ISYMC,ISYMB, 182 * ISPINA,ISPINC,ISPINB, 183 * VECC,VECB,VECA,ATEST, 184 * WRK(KS3),XINDX,UDV,MJWOP,WRK(KWRK),LWRKS) 185C 186C Multiply with omega-c 187C 188 CALL DAXPY(KZYVA,OMEGAC,WRK(KS3),1,WRK(KE3),1) 189 END IF 190C 191C Result is now in WRK(KE3) 192C 193 IF (IPRRSP .GT. 100 ) THEN 194 KZVA = MZVAR(ISYMA) 195 WRITE(LUPRI,'(//A/A)') 196 & ' Final result in T3DRV (Col 1 = Z, Col 2 = Y)', 197 & ' ======================' 198 CALL OUTPUT(WRK(KE3),1,KZVA,1,2,KZVA,2,1,LUPRI) 199 END IF 200C 201 CALL QEXIT('T3DRV') 202 RETURN 203 END 204C /* Deck dft3drv */ 205 SUBROUTINE DFT3DRV(VECB, VECC,FI,FA,ZYMB,ZYMC, 206 & KZYVA,KZYVB,KZYVC, 207 & ISYMA,ISYMB,ISYMC, 208 & ISPINA,ISPINB,ISPINC, 209 & CMO,MJWOP,WRK,LFREE,ADDFOCK) 210#include "implicit.h" 211#include "infvar.h" 212#include "inforb.h" 213#include "maxorb.h" 214#include "infinp.h" 215#include "dftcom.h" 216 DIMENSION VECB(*),VECC(*),ZYMB(*),ZYMC(*),WRK(*),CMO(*) 217 DIMENSION FI(NORBT,NORBT),FA(NORBT,NORBT),MJWOP(2,MAXWOP,8) 218 LOGICAL ADDFOCK 219c unpack response vectors. 220 NSIM = 1 221 CALL GTZYMT(NSIM,VECB,KZYVB,ISYMB,ZYMB,MJWOP) 222 CALL GTZYMT(NSIM,VECC,KZYVC,ISYMC,ZYMC,MJWOP) 223c compute dft contribution 224C ZR testing dft_qr_ab 225 IF (NASHT.GT.0) THEN 226 CALL dft_qr_ab(FI,FA,CMO,ZYMB,ISYMB,ISPINB,ZYMC,ISYMC,ISPINC, 227 & ADDFOCK,WRK,LFREE,IPRDFT) 228 ELSE 229 CALL DFTQRC(FI,CMO,ZYMB,ISYMB,ISPINB,ZYMC,ISYMC,ISPINC, 230 & ADDFOCK,WRK,LFREE,IPRDFT) 231 END IF 232c pack stuff back. 233 END 234C /* Deck e3init */ 235 SUBROUTINE E3INIT(VECB,VECC,VECA,ATEST,IBEQC, 236 * ETRS,XINDX,UDV,PV,WRK,LWRK,KZYVA,KZYVB,KZYVC, 237 * ISYMA,ISYMB,ISYMC,ISPINA,ISPINB,ISPINC,CMO,FC, 238 * FV,MJWOP) 239 240C Compute ( E3(jkl) + E3(jlk) ) N(l) N(k) 241C 242C VECB = N(l) (input) 243C VECC = N(k) (input) 244C VECA = N(j) (input, only used if ATEST true) 245C ATEST (input, true to print individiual E3 contributions based on VECA) 246C IBEQC (input, true if VECB = VECC) 247C ETRS = ( E3(jkl) + E3(jlk) ) N(l) N(k) (output) 248 249 use qmcmm_response, only: qmcmm_qr 250 use pelib_interface, only: use_pelib, pelib_ifc_qro, 251 & pelib_ifc_gspol 252C 253#include "implicit.h" 254#include "priunit.h" 255C 256 PARAMETER ( D1 = 1.0D0 ) 257C 258C infinp.h : HSROHF, DODFT, ? 259C 260#include "infdim.h" 261#include "maxorb.h" 262#include "maxash.h" 263#include "infinp.h" 264#include "inforb.h" 265#include "infvar.h" 266#include "infrsp.h" 267#include "rspprp.h" 268#include "infcr.h" 269#include "inftpa.h" 270#include "inftmo.h" 271#include "pcmlog.h" 272#include "infpp.h" 273#include "infsmo.h" 274#include "infhyp.h" 275#include "esg.h" 276#include "dftcom.h" 277#include "absorp.h" 278#include "abslrs.h" 279#include "abscrs.h" 280#include "gnrinf.h" 281#include "infpar.h" 282C 283 LOGICAL ATEST, TRPLETSAVE, DFTSAV, ADDFOCK 284 REAL*8 VECB(*), VECC(*), VECA(*), ETRS(KZYVA) 285 REAL*8 XINDX(*), UDV(*), PV(*), WRK(*) 286 REAL*8 CMO(*), FC(*), FV(*) 287 INTEGER MJWOP(2,MAXWOP,8) 288C 289C Initialise variables and layout some workspace 290C 291 CALL QENTER('E3INIT') 292 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 293 CALL QUIT('srDFT not implemented in E3INIT') 294 END IF 295 IF (CISRPA .OR. TDA) THEN 296 ! E[3] is zero for CI, i.e. also for CIS = TDA 297 IF (EMBEDDING) THEN 298 ! embedding may give an E[3] contribution ... 299 CALL QUIT('QR embedding not implemented for CISRPA or TDA') 300 END IF 301 ETRS(1:KZYVA) = 0.0D0 302 GO TO 9000 303 END IF 304 KFREE = 1 305 LFREE = LWRK 306 CALL MEMGET('REAL',KZYMB,NORBT*NORBT,WRK,KFREE,LFREE) 307 CALL MEMGET('REAL',KZYMC,NORBT*NORBT,WRK,KFREE,LFREE) 308 CALL MEMGET('REAL',KFC,NORBT*NORBT,WRK,KFREE,LFREE) 309 IF (NASHT.GT.0) THEN 310 CALL MEMGET('REAL',KFO,NORBT*NORBT,WRK,KFREE,LFREE) 311 ELSE 312 CALL MEMGET('REAL',KFO,0,WRK,KFREE,LFREE) 313 END IF 314C 315C Zero Q matrix for ORBEX (and FO for closed shell) 316C 317 CALL MEMGET('REAL',KDUM,N2ORBX,WRK,KFREE,LFREE) 318 CALL DZERO(WRK(KDUM),N2ORBX) 319 IF (NASHT.EQ.0) KFO=KDUM 320C 321C Fock matrices 322C 323 IF (TDHF.AND.(NASHT.EQ.0.OR.HSROHF)) THEN ! excluding single open shell RHF with .OPEN SHELL 324 IF (CRCAL.OR.TOMOM.OR.TPAMP.OR.ESG) THEN 325 CALL C3FOCK(IBEQC,WRK(KFC),VECB,VECC, 326 & WRK(KZYMB),WRK(KZYMC), 327 & KZYVB,KZYVC, 328 & ISYMA,ISYMB,ISYMC, 329 & CMO,FC,MJWOP,WRK(KFREE),LFREE) 330 ELSE IF (HYPCAL.OR.SOMOM.OR.EXMOM.OR.ABS_BETA 331 & .OR.ABS_MCD.OR.ABSLRS_MCD.OR.ABSLRS_MCHD 332 & .OR.ABSLRS_NSCD.OR.ABSLRS_GAMMA.OR.ABSLRS_IDRI) THEN 333 334 DFTSAV=DODFT 335 DODFT=.FALSE. 336 CALL Q3FOCK( VECB,VECC, 337 & KZYVB,KZYVC, ISYMB,ISYMC, ISPINB,ISPINC, 338 & CMO, MJWOP, FC, FV, 339 & WRK(KFC),WRK(KFO), 340 & WRK(KFREE),LFREE ) 341 DODFT=DFTSAV 342 ELSE 343 CALL QUIT('ERROR in E3INIT for HF or DFT') 344 END IF 345C 346 TRPLETSAVE=TRPLET 347 TRPLET=MOD(ISPINA+ISPINB+ISPINC,2).NE.0 348 CALL ORBEX( 349 & ISYMA,1,KZYVA,WRK(KFC),WRK(KFO),WRK(KDUM),WRK(KDUM), 350 & ETRS,D1,UDV,MJWOP 351 & ) 352 TRPLET=TRPLETSAVE 353C 354C Special consideration, ingoing FC does contain the DFT contribution 355C which is transformed in Q3FOCK as [k,k,FC] so the corresponding 356C contribution must be skipped in DFT3DRV 357C 358 ADDFOCK=.FALSE. 359C 360 ELSE 361 ADDFOCK=.NOT.DIRFCK 362 CALL MEMGET('REAL',KDEN1,NASHT*NASHT,WRK,KFREE,LFREE) 363 IF (TRPFLG) THEN 364 CALL MEMGET('REAL',KDEN2,2*N2ASHX*N2ASHX,WRK,KFREE,LFREE) 365 ELSE 366 CALL MEMGET('REAL',KDEN2,N2ASHX*N2ASHX,WRK,KFREE,LFREE) 367 END IF 368 CALL E3DRV(VECB, VECC, VECA, ATEST, IBEQC, 369 * ETRS,XINDX,WRK(KFC),WRK(KZYMB),WRK(KZYMC), 370 * WRK(KDEN1),WRK(KDEN2),UDV,PV,WRK(KFREE),LFREE, 371 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,FC,MJWOP) 372 END IF 373C 374C DFT contribution 375C 376 IF(DODFT) THEN 377 CALL DZERO(WRK(KFC),N2ORBX) 378 IF (NASHT.GT.0) CALL DZERO(WRK(KFO),N2ORBX) 379 CALL DFT3DRV(VECB,VECC,WRK(KFC),WRK(KFO),WRK(KZYMB),WRK(KZYMC), 380 * KZYVA,KZYVB,KZYVC, 381 * ISYMA,ISYMB,ISYMC, 382 * ISPINA,ISPINB,ISPINC, 383 * CMO,MJWOP,WRK(KFREE),LFREE,ADDFOCK) 384 IF (E3TEST) THEN 385 E3 = -DDOT(KZYVA,ETRS,1,VECA,1) 386 END IF 387 TRPLETSAVE=TRPLET 388 TRPLET=MOD(ISPINA+ISPINB+ISPINC,2).NE.0 389 CALL ORBEX(ISYMA,1,KZYVA,WRK(KFC),WRK(KFO),WRK(KDUM), 390 & WRK(KDUM),ETRS,D1,UDV,MJWOP) 391 TRPLET=TRPLETSAVE 392 IF (E3TEST) THEN 393 E3TOT = -DDOT(KZYVA,ETRS,1,VECA,1) 394 E3DFT = E3TOT - E3 395 WRITE(LUPRI,'(/A,2F20.8)') 396 * ' COU CONTRIBUTION TO HYPVAL',E3,E3 397 WRITE(LUPRI,'(/A,2F20.8)') 398 * ' DFT CONTRIBUTION TO HYPVAL',E3DFT,E3TOT 399 END IF 400 END IF 401 402 IF (USE_PELIB()) THEN 403 IF (.NOT. TDHF .AND. .NOT. TRPLET) THEN 404 CALL QUIT('ERROR: PE-MCSCF quadratic response not'// 405 & ' implemented using the PE library.') 406 ELSE IF (NASHT > 0) THEN 407 CALL QUIT('ERROR: quadratic response not implemented for'// 408 & ' open-shell systems using the PE library.') 409 ELSE IF (TRPLET) THEN 410 CALL QUIT('ERROR: triplets not implemented for quadratic'// 411 & ' response using the PE library.') 412 ENDIF 413 IF (.NOT. PELIB_IFC_GSPOL()) THEN 414 CALL PELIB_IFC_QRO(VECB, VECC, ETRS, XINDX, WRK(KZYMB), 415 & WRK(KZYMC), UDV, WRK(KFREE), LFREE, 416 & KZYVA, KZYVB, KZYVC, ISYMA, ISYMB, 417 & ISYMC, CMO, MJWOP) 418 END IF 419 END IF 420 421 IF (QM3) THEN 422 IF (TDHF) THEN 423#if defined(VAR_MPI) 424 IF (NODTOT.GE.1) THEN 425 CALL QM3QRO_P(VECB,VECC, 426 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 427 * UDV,WRK(KFREE),LFREE, 428 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP) 429 CALL QM3QRO_P(VECC,VECB, 430 * ETRS,XINDX,WRK(KZYMC),WRK(KZYMB), 431 * UDV,WRK(KFREE),LFREE, 432 * KZYVA,KZYVC,KZYVB,ISYMA,ISYMC,ISYMB,CMO,MJWOP) 433 ELSE 434#endif 435 CALL QM3QRO(VECB,VECC, 436 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 437 * UDV,WRK(KFREE),LFREE, 438 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP) 439 CALL QM3QRO(VECC,VECB, 440 * ETRS,XINDX,WRK(KZYMC),WRK(KZYMB), 441 * UDV,WRK(KFREE),LFREE, 442 * KZYVA,KZYVC,KZYVB,ISYMA,ISYMC,ISYMB,CMO,MJWOP) 443#if defined(VAR_MPI) 444 ENDIF 445#endif 446 END IF 447 END IF 448 449 IF (QMMM) THEN 450 IF (TDHF) THEN 451 CALL QMMMQRO(VECB,VECC, 452 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 453 * UDV,WRK(KFREE),LFREE, 454 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP, 455 * ISPINA,ISPINB,ISPINC) 456 CALL QMMMQRO(VECC,VECB, 457 * ETRS,XINDX,WRK(KZYMC),WRK(KZYMB), 458 * UDV,WRK(KFREE),LFREE, 459 * KZYVA,KZYVC,KZYVB,ISYMA,ISYMC,ISYMB,CMO,MJWOP, 460 * ISPINA,ISPINC,ISPINB) 461 END IF 462 END IF 463C 464 IF (QMNPMM) THEN 465 IF (TDHF) THEN 466 CALL qmcmm_qr(VECB,VECC,ETRS,XINDX,WRK(KZYMB), 467 * WRK(KZYMC),UDV,WRK(KFREE),LFREE, KZYVA, 468 * KZYVB,KZYVC, ISYMA,ISYMB,ISYMC,CMO,MJWOP, 469 * ISPINA,ISPINB,ISPINC) 470 CALL qmcmm_qr(VECC,VECB,ETRS,XINDX,WRK(KZYMC), 471 * WRK(KZYMB), UDV,WRK(KFREE),LFREE,KZYVA, 472 * KZYVC,KZYVB,ISYMA,ISYMC,ISYMB,CMO,MJWOP, 473 * ISPINA,ISPINC,ISPINB) 474 END IF 475 END IF 476C 477 IF (FLAG(16)) THEN 478 CALL MEMGET('INTE',KSYRLM,(2*LSOLMX+1),WRK,KFREE,LFREE) 479 IF (TDHF) THEN 480 CALL C3SOL(VECA, VECB, VECC, 481 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 482 * UDV,WRK(KFREE),LFREE, 483 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP, 484 * WRK(KSYRLM)) 485 CALL C3SOL(VECA, VECC, VECB, 486 * ETRS,XINDX,WRK(KZYMC),WRK(KZYMB), 487 * UDV,WRK(KFREE),LFREE, 488 * KZYVA,KZYVC,KZYVB,ISYMA,ISYMC,ISYMB,CMO,MJWOP, 489 * WRK(KSYRLM)) 490 ELSE 491 CALL E3SOL(VECA, VECB, VECC, 492 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 493 * WRK(KDEN1),UDV,WRK(KFREE),LFREE, 494 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP, 495 * WRK(KSYRLM)) 496 END IF 497C 498C 499Clf 10/06/2002 500C quadratic PCM response (for the time beeing without memory efficent SCF) 501C 502 ELSE IF (PCM) THEN 503 IF (.NOT.NEWQR) THEN 504 CALL E3IEF(VECA, VECB, VECC, 505 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 506 * WRK(KDEN1),UDV,WRK(KFREE),LFREE, 507 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP) 508 ELSE 509 CALL E3IEF2(VECA, VECB, VECC, 510 * ETRS,XINDX,WRK(KZYMB),WRK(KZYMC), 511 * WRK(KDEN1),UDV,WRK(KFREE),LFREE, 512 * KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC,CMO,MJWOP) 513 END IF 514 END IF 515C 516 9000 CALL QEXIT('E3INIT') 517 RETURN 518 END 519C /* Deck e3drv */ 520 SUBROUTINE E3DRV(VECB, VECC, VECA, ATEST, IBEQC, 521 * ETRS,XINDX,FI,ZYMB,ZYMC,DEN1,DEN2, 522 * UDV,PV,WRK,LWRK,KZYVA,KZYVB,KZYVC, 523 * ISYMA,ISYMB,ISYMC,CMO,FC,MJWOP) 524C 525C Purpose: 526C Outer driver routine for E[3] times two vectors. This subroutine 527C calls the setup of the integral transformation, constructs a den- 528C sity matrix if necessary and calls RSP2GR to compute the gradient. 529C 530#include "implicit.h" 531C 532 PARAMETER ( D0 = 0.0D0, DH = 0.5D0 , D1 = 1.0D0 , D2 = 2.0D0 ) 533 PARAMETER ( DM1=-1.0D0 ) 534C 535#include "priunit.h" 536#include "infdim.h" 537#include "inforb.h" 538#include "wrkrsp.h" 539#include "infrsp.h" 540#include "infpri.h" 541#include "infvar.h" 542#include "qrinf.h" 543#include "infspi.h" 544#include "infden.h" 545C 546 DIMENSION ETRS(KZYVA), MJWOP(2,MAXWOP,8) 547 DIMENSION XINDX(*) 548 DIMENSION UDV(NASHDI,NASHDI) 549 DIMENSION PV(N2ASHX*N2ASHX,*) 550 DIMENSION DEN1(NASHDI,NASHDI), DEN2(N2ASHX*N2ASHX,*) 551 DIMENSION FI(NORBT,NORBT) 552 DIMENSION ZYMB(*), ZYMC(*) 553 DIMENSION WRK(*) 554 DIMENSION CMO(*),FC(*) 555 DIMENSION VECB(KZYVB), VECC(KZYVC), VECA(KZYVA) 556C 557 LOGICAL ATEST, LCON, LORB, LREFST 558 LOGICAL TDM, NORHO2, LREF 559 LOGICAL LORBB, LORBC, LCONB, LCONC 560C 561C Initialise variables and layout some workspace 562C 563 CALL QENTER('E3DRV') 564 IGRSYM = ISYMA 565 IGRSPI = ISPINA 566 TDM = .TRUE. 567 NORHO2 = .FALSE. 568 IDAX0 = 0 569C 570C for E3TEST 571C 572 E3TMZC = D0 573 E3TMZO = D0 574 E3TMYC = D0 575 E3TMYO = D0 576 E3TERM = D0 577 KZCA = MZCONF(ISYMA) 578 KZOA = MZWOPT(ISYMA) 579 KAZC = 1 580 KAZO = KAZC + KZCA 581 KAYC = KAZO + KZOA 582 KAYO = KAYC + KZCA 583C 584 KH1 = 1 585 KH1X = KH1 + N2ORBX 586 KWOIT = KH1X + N2ORBX 587 LWOIT = LWRK - KWOIT 588 IF (LWOIT.LT.0) CALL ERRWRK('E3DRV',KWOIT-1,LWRK) 589 KFROIT = 1 590 KCREF = KWOIT 591 KFREE = KCREF + MZCONF(1) 592 LFREE = LWRK - KFREE 593 IF (LFREE.LT.0) CALL ERRWRK('E3DRV',KFREE-1,LWRK) 594C 595 IPROIT = MAX(IPRRSP/5,2) 596C 597 CALL RSPH1(WRK(KH1),CMO,WRK(KWOIT),LWOIT) 598 IF (IPRRSP. GT. 200 ) THEN 599 WRITE(LUPRI,'(//A)') ' One electron matrix in E3DRV' 600 WRITE(LUPRI,'(A)') ' ============================' 601 CALL OUTPUT(WRK(KH1),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 602 END IF 603C 604C Gradient flags 605C 606 607 ISYMST = MULD2H(IGRSYM,IREFSY) 608 IF ( IGRSYM .EQ. 1 ) THEN 609 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 610 ELSE 611 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 612 END IF 613 IF ( ISYMB .EQ. 1 ) THEN 614 LCONB = ( MZCONF(ISYMB) .GT. 1 ) 615 ELSE 616 LCONB = ( MZCONF(ISYMB) .GT. 0 ) 617 END IF 618 IF ( ISYMC .EQ. 1 ) THEN 619 LCONC = ( MZCONF(ISYMC) .GT. 1 ) 620 ELSE 621 LCONC = ( MZCONF(ISYMC) .GT. 0 ) 622 END IF 623 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 624 LORBB = ( MZWOPT(ISYMB) .GT. 0 ) 625 LORBC = ( MZWOPT(ISYMC) .GT. 0 ) 626C 627C Case 1: 628C ======= 629C IF ( MZWOPT(ISYMB) .EQ. 0 ) GO TO 3000 630 IF (.NOT. LORBB) GO TO 3000 631C Transform the integrals with kappa(1) 632C 633 NSIM = 1 634 CALL GTZYMT(NSIM,VECB,KZYVB,ISYMB,ZYMB,MJWOP) 635 JTRLVL= 2 636 IF (.NOT.DIROIT) THEN 637 CALL RSPOIT(IDAX0,IDAX1,JTRLVL,ISYMB,ZYMB, 638 * WRK(KWOIT),KFROIT,LWOIT,IPROIT) 639 END IF 640 CALL DZERO(WRK(KH1X),N2ORBX) 641 CALL OITH1(ISYMB,ZYMB,WRK(KH1),WRK(KH1X),1) 642C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 643 CALL DCOPY(N2ORBX,WRK(KH1X),1,FI,1) 644C 645C IF (MZCONF(ISYMC) .EQ. 0) GO TO 2000 646 IF (.NOT.LCONC) GO TO 2000 647C 648C Construct the density matrix <02L|..|0> + <0|..|02R> 649C 650 CALL GETREF(WRK(KCREF),MZCONF(1)) 651C 652 ILSYM = IREFSY 653 IRSYM = MULD2H(IREFSY,ISYMC) 654 NCL = MZCONF(1) 655 NCR = MZCONF(ISYMC) 656 KZVARL = MZCONF(1) 657 KZVARR = MZYVAR(ISYMC) 658 LREF = .TRUE. 659 ISYMDN = MULD2H(ILSYM,IRSYM) 660 NSIM = 1 661 ISPIN1 = MULSP(IGRSPI,ISPINB) 662 ISPIN2 = 0 663 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 664 * WRK(KCREF),VECC,OVLAP,DEN1,DEN2,ISPIN1,ISPIN2,TDM, 665 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 666 IF (IGRSPI.EQ.1) THEN 667 ISPIN1 = ISPINB 668 ISPIN2 = IGRSPI 669 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 670 * WRK(KCREF),VECC,OVLAP,DEN1,DEN2(1,2),ISPIN1,ISPIN2,TDM, 671 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 672 END IF 673 CALL IPSET(MULSP(IGRSPI,ISPINB),0,ISPINB,IGRSPI) 674C 675C Make the gradient 676C 677 ISYMV = ISYMC 678 INTSYM = ISYMB 679 LREFST = .FALSE. 680 NZYVEC = MZYVAR(ISYMC) 681 NZCVEC = MZCONF(ISYMC) 682 IKLVL = 1 683 IF (DIROIT) THEN 684 ISPIN1 = ISPINB 685 ISPIN2 = 0 686 IDAX = IDAX0 687 KINTSY = 1 688 ELSE 689 IDAX = IDAX1 690 KINTSY = INTSYM 691 END IF 692 CALL RSP2CR(IGRSYM,IGRSPI,ETRS,VECC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 693 * XINDX,OVLAP,DEN1,DEN2,FI,WRK(KWOIT),LWOIT,KZYVA, 694 * LCON,LORB,LREFST,IDAX,KINTSY,ISPIN1,ISPIN2,ISPIN3, 695 * IKLVL,ZYMB,ISYMB,DUM,IDUM,DUM,IDUM,CMO,FC,MJWOP) 696C 697 IF (ATEST) THEN 698 E3NWZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 699 E3NWZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 700 E3NWYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 701 E3NWYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 702 E3NEW = E3NWZC + E3NWZO + E3NWYC + E3NWYO 703 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 1 :',E3NEW-E3TERM 704 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 1 ZC :',E3NWZC-E3TMZC 705 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 1 ZO :',E3NWZO-E3TMZO 706 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 1 YC :',E3NWYC-E3TMYC 707 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 1 YO :',E3NWYO-E3TMYO 708 E3TMZC = E3NWZC 709 E3TMZO = E3NWZO 710 E3TMYC = E3NWYC 711 E3TMYO = E3NWYO 712 E3TERM = E3NEW 713 END IF 714 IF (IPRRSP .GT. 100) THEN 715 WRITE(LUPRI,'(A)') ' Case 1 for E3 term' 716 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 717 END IF 718C 719C Case 2: 720C ======= 721C2000 IF ((MZWOPT(ISYMC) .EQ. 0).OR.(MZWOPT(ISYMB) .EQ. 0)) THEN 722 2000 IF (.NOT.LORBC .OR. .NOT.LORBB) THEN 723 IF (.NOT.DIROIT) CALL OITCLO(IDAX1,'DELETE') 724 GO TO 3000 725 END IF 726C 727C Transform the integrals with 2k,1k 728C 729C Put the factor of one half present in this term into the 730C ZY matrix used for transforming the integrals 731C 732 NSIM = 1 733 CALL GTZYMT(NSIM,VECC,KZYVC,ISYMC,ZYMC,MJWOP) 734C 735C This half would dissappear were we to use the commutator 736C [1k,[2k,H]] = [2k,[1k,H]] + [[2k,1k],H] 737C 738 CALL DSCAL(NORBT*NORBT,DH,ZYMC,1) 739 JTRLVL= 1 740 IF (.NOT.DIROIT) THEN 741 CALL RSPOIT(IDAX1,IDAX21,JTRLVL,ISYMC,ZYMC, 742 * WRK(KWOIT),KFROIT,LWOIT,IPROIT) 743 CALL OITCLO(IDAX1,'DELETE') 744 END IF 745 CALL DZERO(FI,N2ORBX) 746 CALL OITH1(ISYMC,ZYMC,WRK(KH1X),FI,ISYMB) 747C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 748C 749C We have the density matrices over the reference state already 750C 751 OVLAP = D1 752 ISYMDN = 1 753C 754C Make the gradient 755C 756 ISYMV = IREFSY 757 INTSYM = MULD2H(ISYMB,ISYMC) 758 LREFST = .TRUE. 759 NZYVEC = MZCONF(1) 760 NZCVEC = MZCONF(1) 761 VECDUM = D0 762C 763C Vector may be dummy if LREFST = .TRUE. 764C 765C (~~| )e(SbSc,+) + (~|~)e(Sb,Sc) 766C 767 IKLVL=2 768 IF (DIROIT) THEN 769 IDAX = IDAX0 770 KINTSY = 1 771 ELSE 772 IDAX = IDAX21 773 KINTSY = INTSYM 774 END IF 775 CALL IPSET(0,0,1,1) 776 CALL RSP2CR(IGRSYM,IGRSPI,ETRS,VECDUM,NZYVEC,NZCVEC,ISYMV,ISYMDN, 777 * XINDX,OVLAP,UDV,PV,FI,WRK(KWOIT),LWOIT,KZYVA, 778 * LCON,LORB,LREFST,IDAX,KINTSY,ISPINB,ISPINC,ISPIN3, 779 * IKLVL,ZYMB,ISYMB,ZYMC,ISYMC,DUM,IDUM,CMO,FC,MJWOP) 780C 781 IF (.NOT.DIROIT) CALL OITCLO(IDAX21,'DELETE') 782C 783 IF (ATEST) THEN 784 E3NWZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 785 E3NWZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 786 E3NWYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 787 E3NWYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 788 E3NEW = E3NWZC + E3NWZO + E3NWYC + E3NWYO 789 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 2 :',E3NEW-E3TERM 790 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 2 ZC :',E3NWZC-E3TMZC 791 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 2 ZO :',E3NWZO-E3TMZO 792 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 2 YC :',E3NWYC-E3TMYC 793 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 2 YO :',E3NWYO-E3TMYO 794 E3TMZC = E3NWZC 795 E3TMZO = E3NWZO 796 E3TMYC = E3NWYC 797 E3TMYO = E3NWYO 798 E3TERM = E3NEW 799 END IF 800 IF (IPRRSP .GT. 100 ) THEN 801 WRITE(LUPRI,'(A)') ' Case 2 for E3 term' 802 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 803 END IF 804C 805C Case 3: 806C ======= 807C3000 IF ( MZWOPT(ISYMC) .EQ. 0 ) GO TO 5000 8083000 IF ( .NOT. LORBC ) GO TO 5000 809 IF ( IBEQC .EQ. 1 .AND. .NOT.E3TEST) THEN 810 CALL DSCAL(KZYVA,D2,ETRS,1) 811 IF (ATEST .OR. IPRRSP .GT. 100) THEN 812 WRITE(LUPRI,*) 813 & ' E3TEST Case 3 = Case 1 because vecb .eq. vecc' 814 WRITE(LUPRI,*) 815 & ' E3TEST Case 4 = Case 2 because vecb .eq. vecc' 816 E3TMZC = D2*E3TMZC 817 E3TMZO = D2*E3TMZO 818 E3TMYC = D2*E3TMYC 819 E3TMYO = D2*E3TMYO 820 E3TERM = D2*E3TERM 821 END IF 822 GO TO 5000 823 ELSE IF ( IBEQC .EQ. -1 .AND. .NOT.E3TEST) THEN 824 KZVA = MZVAR(ISYMA) 825 CALL DAXPY(KZVA,DM1,ETRS(KZVA+1),1,ETRS,1) 826 CALL DCOPY(KZVA,ETRS,1,ETRS(KZVA+1),1) 827 CALL DSCAL(KZVA,DM1,ETRS(KZVA+1),1) 828 IF (ATEST .OR. IPRRSP .GT. 100) THEN 829 WRITE(LUPRI,*) 830 * ' E3TEST Case 3 = (Case 1)T because vecb .eq. -(vecc)T' 831 WRITE(LUPRI,*) 832 * ' E3TEST Case 4 = (Case 2)T because vecb .eq. -(vecc)T' 833 E3TMZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 834 E3TMZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 835 E3TMYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 836 E3TMYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 837 E3TERM = E3TMZC + E3TMZO + E3TMYC + E3TMYO 838 END IF 839 GO TO 5000 840 END IF 841C 842C Transform the integrals with kappa(2) 843C 844 NSIM = 1 845 CALL GTZYMT(NSIM,VECC,KZYVC,ISYMC,ZYMC,MJWOP) 846 JTRLVL= 2 847 IF (.NOT.DIROIT) THEN 848 CALL RSPOIT(IDAX0,IDAX2,JTRLVL,ISYMC,ZYMC, 849 * WRK(KWOIT),KFROIT,LWOIT,IPROIT) 850 END IF 851 CALL DZERO(WRK(KH1X),N2ORBX) 852 CALL OITH1(ISYMC,ZYMC,WRK(KH1),WRK(KH1X),1) 853C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 854 CALL DCOPY(N2ORBX,WRK(KH1X),1,FI,1) 855C 856C IF (MZCONF(ISYMB) .EQ. 0) GO TO 4000 857 IF (.NOT.LCONB) GO TO 4000 858C 859C Construct the density matrix <01L|..|0> + <0|..|01R> 860C 861 CALL GETREF(WRK(KCREF),MZCONF(1)) 862C 863 ILSYM = IREFSY 864 IRSYM = MULD2H(ISYMB,IREFSY) 865 NCL = MZCONF(1) 866 NCR = MZCONF(ISYMB) 867 KZVARL = MZCONF(1) 868 KZVARR = MZYVAR(ISYMB) 869 LREF = .TRUE. 870 ISYMDN = MULD2H(ILSYM,IRSYM) 871 NSIM = 1 872 ISPIN1 = MULSP(IGRSPI,ISPINC) 873 ISPIN2 = 0 874 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 875 * WRK(KCREF),VECB,OVLAP,DEN1,DEN2,ISPIN1,ISPIN2,TDM, 876 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 877 IF (IGRSPI.EQ.1) THEN 878 ISPIN1 = ISPINC 879 ISPIN2 = IGRSPI 880 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 881 * WRK(KCREF),VECB,OVLAP,DEN1,DEN2(1,2),ISPIN1,ISPIN2,TDM, 882 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 883 END IF 884 CALL IPSET(MULSP(IGRSPI,ISPINC),0,ISPINC,IGRSPI) 885C 886C Make the gradient 887C 888 ISYMV = ISYMB 889 INTSYM = ISYMC 890 LREFST = .FALSE. 891 NZYVEC = MZYVAR(ISYMB) 892 NZCVEC = MZCONF(ISYMB) 893 IKLVL = 1 894 IF (DIROIT) THEN 895 ISPIN1 = ISPINC 896 ISPIN2 = 0 897 IDAX = IDAX0 898 KINTSY = 1 899 ELSE 900 IDAX = IDAX2 901 KINTSY = INTSYM 902 END IF 903 CALL RSP2CR(IGRSYM,IGRSPI,ETRS,VECB,NZYVEC,NZCVEC,ISYMV,ISYMDN, 904 * XINDX,OVLAP,DEN1,DEN2,FI,WRK(KWOIT),LWOIT,KZYVA, 905 * LCON,LORB,LREFST,IDAX,KINTSY,ISPIN1,ISPIN2,ISPIN3, 906 * IKLVL,ZYMC,ISYMC,DUM,IDUM,DUM,IDUM,CMO,FC,MJWOP) 907C 908 IF (ATEST) THEN 909 E3NWZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 910 E3NWZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 911 E3NWYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 912 E3NWYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 913 E3NEW = E3NWZC + E3NWZO + E3NWYC + E3NWYO 914 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 3 :',E3NEW-E3TERM 915 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 3 ZC :',E3NWZC-E3TMZC 916 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 3 ZO :',E3NWZO-E3TMZO 917 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 3 YC :',E3NWYC-E3TMYC 918 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 3 YO :',E3NWYO-E3TMYO 919 E3TMZC = E3NWZC 920 E3TMZO = E3NWZO 921 E3TMYC = E3NWYC 922 E3TMYO = E3NWYO 923 E3TERM = E3NEW 924 END IF 925 IF (IPRRSP .GT. 100 ) THEN 926 WRITE(LUPRI,'(A)') ' Case 3 for E3 term' 927 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 928 END IF 929C 930C Case 4: 931C ======= 932C4000 IF ((MZWOPT(ISYMB) .EQ. 0).OR.(MZWOPT(ISYMC) .EQ. 0 )) THEN 933 4000 IF (.NOT.LORBB .OR. .NOT.LORBC) THEN 934 IF (.NOT.DIROIT) CALL OITCLO(IDAX2,'DELETE') 935 GO TO 5000 936 END IF 937C 938C 939C Put the factor of one half present in this term into the 940C ZY matrix used for transforming the integrals 941C 942 NSIM = 1 943 CALL GTZYMT(NSIM,VECB,KZYVB,ISYMB,ZYMB,MJWOP) 944 CALL DSCAL(NORBT*NORBT,DH,ZYMB,1) 945 JTRLVL= 1 946 IF (.NOT.DIROIT) THEN 947 CALL RSPOIT(IDAX2,IDAX12,JTRLVL,ISYMB,ZYMB, 948 * WRK(KWOIT),KFROIT,LWOIT,IPROIT) 949 CALL OITCLO(IDAX2,'DELETE') 950 END IF 951 CALL DZERO(FI,N2ORBX) 952 CALL OITH1(ISYMB,ZYMB,WRK(KH1X),FI,ISYMC) 953C CALL OITH1(LSYOIT,OITMAT,H1,H1X,IH1SYM) 954C 955C We have the density matrices over the reference state already 956C 957 OVLAP = D1 958 ISYMDN = 1 959C 960C Make the gradient 961C 962 ISYMV = IREFSY 963 INTSYM = MULD2H(ISYMB,ISYMC) 964 LREFST = .TRUE. 965 NZYVEC = MZCONF(1) 966 NZCVEC = MZCONF(1) 967 VECDUM = D0 968C 969C (~~| )e(ScSb,+) + (~|~)e(Sc,Sb) 970C 971 IKLVL=2 972 IF (DIROIT) THEN 973 IDAX = IDAX0 974 KINTSY = 1 975 ELSE 976 IDAX = IDAX12 977 KINTSY = INTSYM 978 END IF 979 CALL IPSET(0,0,1,1) 980 CALL RSP2CR(IGRSYM,IGRSPI,ETRS,VECDUM,NZYVEC,NZCVEC,ISYMV,ISYMDN, 981 * XINDX,OVLAP,UDV,PV,FI,WRK(KWOIT),LWOIT,KZYVA, 982 * LCON,LORB,LREFST,IDAX,KINTSY,ISPINC,ISPINB,ISPIN3, 983 * IKLVL,ZYMC,ISYMC,ZYMB,ISYMB,DUM,IDUM,CMO,FC,MJWOP) 984C 985 IF (.NOT.DIROIT) CALL OITCLO(IDAX12,'DELETE') 986C 987 IF (ATEST) THEN 988 E3NWZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 989 E3NWZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 990 E3NWYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 991 E3NWYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 992 E3NEW = E3NWZC + E3NWZO + E3NWYC + E3NWYO 993 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 4 :',E3NEW-E3TERM 994 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 4 ZC :',E3NWZC-E3TMZC 995 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 4 ZO :',E3NWZO-E3TMZO 996 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 4 YC :',E3NWYC-E3TMYC 997 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 4 YO :',E3NWYO-E3TMYO 998 E3TMZC = E3NWZC 999 E3TMZO = E3NWZO 1000 E3TMYC = E3NWYC 1001 E3TMYO = E3NWYO 1002 E3TERM = E3NEW 1003 END IF 1004 IF (IPRRSP .GT. 100 ) THEN 1005 WRITE(LUPRI,'(A)') ' Case 4 for E3 term' 1006 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 1007 END IF 1008C 1009C Case 5: 1010C ======= 1011C5000 IF (MZCONF(ISYMB) .EQ. 0 .OR. MZCONF(ISYMC) .EQ. 0) GO TO 6000 10125000 IF (.NOT.LCONB .OR. .NOT.LCONC) GO TO 6000 1013C 1014 CALL DCOPY(N2ORBX,WRK(KH1),1,FI,1) 1015C 1016C Construct <01L|..|02R> + <02L|..|01R> 1017C 1018 ILSYM = MULD2H(IREFSY,ISYMB) 1019 IRSYM = MULD2H(IREFSY,ISYMC) 1020 NCL = MZCONF(ISYMB) 1021 NCR = MZCONF(ISYMC) 1022 KZVARL = MZYVAR(ISYMB) 1023 KZVARR = MZYVAR(ISYMC) 1024 LREF = .FALSE. 1025 ISYMDN = MULD2H(ILSYM,IRSYM) 1026 NSIM = 1 1027 ISPIN1 = IGRSPI 1028 ISPIN2 = 0 1029 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 1030 * VECB,VECC,OVLAP,DEN1,DEN2,ISPIN1,ISPIN2,TDM,NORHO2, 1031 * XINDX,WRK,KFREE,LFREE,LREF) 1032 IF (IGRSPI.EQ.1) THEN 1033 ISPIN1 = 0 1034 ISPIN2 = IGRSPI 1035 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 1036 * VECB,VECC,OVLAP,DEN1,DEN2(1,2),ISPIN1,ISPIN2,TDM,NORHO2, 1037 * XINDX,WRK,KFREE,LFREE,LREF) 1038 END IF 1039 CALL IPSET(IGRSPI,0,0,IGRSPI) 1040C 1041C Make the gradient 1042C 1043 KFREE = 1 1044 LFREE = LWRK - KFREE 1045C 1046 LCON = .FALSE. 1047 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 1048 LREFST = .FALSE. 1049 IDUM = 1 1050 VECDUM = D0 1051C Calls for configuration part of lt are dummies: 1052 INTSYM = 1 1053 IKLVL = 0 1054 ISPIN1 = 0 1055 ISPIN2 = 0 1056 CALL RSP2CR(IGRSYM,IGRSPI,ETRS,VECDUM,IDUM,IDUM,IDUM,ISYMDN, 1057 * XINDX,OVLAP,DEN1,DEN2,FI,WRK(KFREE),LFREE,KZYVA, 1058 * LCON,LORB,LREFST,IDAX0,INTSYM,ISPIN1,ISPIN2,ISPIN3, 1059 * IKLVL,DUM,IDUM,DUM,IDUM,DUM,IDUM,CMO,FC,MJWOP) 1060C 1061 IF (ATEST) THEN 1062 E3NWZC = -DDOT(KZCA,VECA(KAZC),1,ETRS(KAZC),1) 1063 E3NWZO = -DDOT(KZOA,VECA(KAZO),1,ETRS(KAZO),1) 1064 E3NWYC = -DDOT(KZCA,VECA(KAYC),1,ETRS(KAYC),1) 1065 E3NWYO = -DDOT(KZOA,VECA(KAYO),1,ETRS(KAYO),1) 1066 E3NEW = E3NWZC + E3NWZO + E3NWYC + E3NWYO 1067 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 5 :',E3NEW-E3TERM 1068 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 5 ZC :',E3NWZC-E3TMZC 1069 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 5 ZO :',E3NWZO-E3TMZO 1070 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 5 YC :',E3NWYC-E3TMYC 1071 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Case 5 YO :',E3NWYO-E3TMYO 1072 E3TMZC = E3NWZC 1073 E3TMZO = E3NWZO 1074 E3TMYC = E3NWYC 1075 E3TMYO = E3NWYO 1076 E3TERM = E3NEW 1077 END IF 1078 IF (IPRRSP .GT. 100 ) THEN 1079 WRITE(LUPRI,'(A)') ' Case 5 for E3 term added' 1080 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 1081 END IF 1082C 1083C End of transformation 1084C 1085 6000 CONTINUE 1086 IF (ATEST) THEN 1087 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Final result :',E3TERM 1088 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Final result ZC :',E3TMZC 1089 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Final result ZO :',E3TMZO 1090 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Final result YC :',E3TMYC 1091 WRITE(LUPRI,'(A,F24.8)') ' E3TEST Final result YO :',E3TMYO 1092 END IF 1093C 1094 CALL QEXIT('E3DRV') 1095 RETURN 1096 END 1097C /* Deck rsp2cr */ 1098 SUBROUTINE RSP2CR(IGRSYM,IGRSPI, 1099 * ETRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 1100 * XINDX,OVLAP,DEN1,DEN2,FI,WRK,LWRK,KZYVA, 1101 * LCON,LORB,LREFST,IDAX,INTSYM, 1102 * ISPIN1,ISPIN2,ISPIN3, 1103 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3, 1104 * CMO,FC,MJWOP) 1105C 1106C Layout the core for RSP2GR 1107C 1108#include "implicit.h" 1109C 1110#include "infdim.h" 1111#include "inforb.h" 1112#include "infvar.h" 1113#include "infrsp.h" 1114#include "wrkrsp.h" 1115#include "qrinf.h" 1116#include "infpri.h" 1117#include "infspi.h" 1118C 1119 DIMENSION WRK(*) 1120 DIMENSION ETRS(KZYVA) 1121 DIMENSION XINDX(*) 1122 DIMENSION DEN1(*), DEN2(N2ASHX*N2ASHX,*) 1123 DIMENSION FI(*) 1124 DIMENSION VEC(NZYVEC),CMO(*),FC(*) 1125 DIMENSION ZYMAT1(*), ZYMAT2(*), ZYMAT3(*) 1126C 1127 LOGICAL LCON, LORB, LREFST, HSORUN 1128C 1129 CALL QENTER('RSP2CR') 1130C 1131C Layout of WRK: 1132C We keep the H2AX array at the beginning, in order to 1133C release extra workspace in the gradient routine after processing 1134C the orbital part of the linear transformation. 1135C 1136 KH2AX = 1 1137 IF (LCON) THEN 1138 KFA = KH2AX + N2ASHX * N2ASHX 1139 IF (DIROIT) THEN 1140 IF (IKLVL.EQ.2) KFA = KH2AX + N2ASHX * N2ASHX * 2 1141 IF (IKLVL.EQ.3) KFA = KH2AX + N2ASHX * N2ASHX * 4 1142 END IF 1143 ELSE 1144 KFA = KH2AX 1145 END IF 1146 KQA = KFA + NORBT * NORBT 1147 KQB = KQA + NORBT * NASHDI 1148 KPVD = KQB + NORBT * NASHDI 1149 KH2 = KPVD + N2ASHX * N2ASHX 1150 KH2X = KH2 + N2ORBX 1151 KWRKO = KH2X + N2ORBX 1152C 1153C Get free workspace for orbital part of linear transformation 1154C 1155 LWRKO = LWRK - KWRKO 1156 IF (LWRKO.LT.0) CALL ERRWRK('RSP2CR 1',KWRKO-1,LWRK) 1157C 1158C Get space that can be used in configuration part 1159C We need the arrays H2XAC and FI, so keep them intact 1160C 1161 KWRKC = KFA 1162 LWRKC = LWRK - KWRKC 1163 IF (LWRKC.LT.0) CALL ERRWRK('RSP2CR 2',KWRKC-1,LWRK) 1164 HSORUN = .FALSE. 1165 CALL RSP2GR (IGRSYM,IGRSPI,ETRS,VEC,NZYVEC,NZCVEC,ISYMV,ISYMDN, 1166 * FI,WRK(KFA),WRK(KQA),WRK(KQB),WRK(KH2),WRK(KH2X), 1167 * OVLAP,DEN1,DEN2,WRK(KPVD),WRK(KH2AX),XINDX, 1168 * WRK(KWRKO),LWRKO,WRK(KWRKC),LWRKC,KZYVA, 1169 * LCON,LORB,LREFST,IDAX,INTSYM,ISPIN1,ISPIN2,ISPIN3, 1170 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3, 1171 * HSORUN,CMO,FC,MJWOP) 1172C 1173 CALL QEXIT('RSP2CR') 1174 RETURN 1175 END 1176C /* Deck rsp2gr */ 1177 SUBROUTINE RSP2GR(IGRSYM,IGRSPI,ETRS,VEC,NZYVEC,NZCVEC,ISYMV, 1178 * ISYMDN, 1179 * FI,FA,QA,QB,H2,H2X, 1180 * OVLAP,DEN1,DEN2,PVDEN,H2AX,XINDX, 1181 * WRKO,LWRKO,WRKC,LWRKC,KZYVA, 1182 * LCON,LORB,LREFST,IDAX,INTSYM, 1183 * ISPIN1,ISPIN2,ISPIN3, 1184 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3, 1185 * HSORUN,CMO,FC,MJWOP) 1186C 1187C This routine computes the gradient, given a density matrix, 1188C the two vectors and a set of (transformed) integrals. In the deriva- 1189C vation of the implemented formulas it is assumed that the integrals 1190C have particle permutation symmetry. 1191C 1192C The result is returned added to the values present in ETRS upon 1193C calling. 1194C 1195C Memory layout: 1196C 1. WRKO is the workarray for the orbital part of the lt. 1197C 2. WRKC is the workarray for the CI part of the lt. Since we do 1198C not need the active fock matrix and the Q matrices here, 1199C some space can be released. Organisation should be 1200C done by the calling routine. (See E3NBNC for example) 1201C 1202#include "implicit.h" 1203C 1204 PARAMETER ( D0 = 0.0D0, DH = 0.50D0 ) 1205C 1206C INFINP : DIRFCK 1207C 1208#include "dftcom.h" 1209#include "maxorb.h" 1210#include "maxash.h" 1211#include "priunit.h" 1212#include "wrkrsp.h" 1213#include "infrsp.h" 1214#include "infvar.h" 1215#include "qrinf.h" 1216#include "inforb.h" 1217#include "infind.h" 1218#include "infdim.h" 1219#include "infpri.h" 1220#include "infspi.h" 1221#include "infinp.h" 1222C 1223 DIMENSION ETRS(KZYVA), MJWOP(2,MAXWOP,8) 1224 DIMENSION FI(NORBT,NORBT),FA(NORBT,NORBT) 1225 DIMENSION QA(NORBT,NASHDI),QB(NORBT,NASHDI) 1226 DIMENSION H2 (NORBT,NORBT) 1227 DIMENSION H2X (NORBT,NORBT) 1228 DIMENSION H2AX(NASHDI,NASHDI,NASHDI,NASHDI,*) 1229 DIMENSION DEN1(NASHDI,NASHDI), DEN2(NASHDI,NASHDI,NASHDI,NASHDI,*) 1230 DIMENSION PVDEN(NASHDI,NASHDI) 1231 DIMENSION VEC(NZYVEC),CMO(*),FC(*) 1232 DIMENSION XINDX(*) 1233 DIMENSION WRKO(*), WRKC(*) 1234C 1235 LOGICAL LCON, LORB, LREFST, HSORUN, DFTSAV, ADDFOCK 1236C 1237 CALL QENTER('RSP2GR') 1238C 1239 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 1240 CALL QUIT('srDFT not implemented in RSP2GR') 1241 END IF 1242C 1243 CALL DZERO(FA,NORBT*NORBT) 1244 IF (NASHT .GT. 0) THEN 1245 CALL DZERO(QA,NASHT*NORBT) 1246 CALL DZERO(QB,NASHT*NORBT) 1247 IF (LCON) THEN 1248 LH2AX = N2ASHX*N2ASHX 1249 IF (DIROIT) THEN 1250 IF (IKLVL.EQ.2) LH2AX = LH2AX*2 1251 IF (IKLVL.EQ.3) LH2AX = LH2AX*4 1252 END IF 1253 CALL DZERO(H2AX,LH2AX) 1254 END IF 1255 END IF 1256C 1257C 1258 ONEFAC = D0 1259 IF (ISPIN1.EQ.ISPIN2) THEN 1260 DO 105 IW = 1, NISHT 1261 IX = ISX(IW) 1262 ONEFAC = ONEFAC + FI(IX,IX) 1263105 CONTINUE 1264 END IF 1265C 1266 IF( IPRRSP .GT. 200 ) THEN 1267 WRITE(LUPRI,'(/A,F15.8)') ' One electron factor : ',ONEFAC 1268 WRITE(LUPRI,'(/A)') ' One electron matrix in RSP2GR' 1269 WRITE(LUPRI,'(A)') ' =============================' 1270 CALL OUTPUT(FI,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1271 END IF 1272C 1273C Construct FI, FA, QA, QB, H2AX matrices 1274C 1275 KFREE = 1 1276 LFREE = LWRKO 1277 IF (DIROIT) THEN 1278 IF (HSORUN) THEN 1279 CALL HSOFXD (FI,FA,QA,QB,H2AX, 1280 * IDAX,INTSYM, ISYMDN,DEN1,DEN2, 1281 * PVDEN,H2,H2X,WRKO,KFREE,LFREE, 1282 * LCON,LORB,IPRRSP,IGRSPI,ISPIN1,ISPIN2, 1283 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2) 1284 ELSE IF (DIRFCK .AND. NASHT.EQ.0) THEN 1285 DFTSAV = DFTADD 1286 DFTADD = .FALSE. 1287c 1288c note that QFOCK returns complete FOCK/Kohn-Sham matrix, including 1289c DFT contribution, as opposed to RSPFXD/RSPFX which return electronic 1290c part only (w/o DFT contribution). If this is the case (i.e. .NOT.DIRFCK) 1291c we add the DFT contribution to KS matrix later in DFTQRC. 1292c 1293 IF (NASHT.GT.0) THEN 1294 !This point in code is reached for the special case 1295 !single open shell direct HF 1296 WRITE(LUPRI,*)'Input combination not implemented' 1297 WRITE(LUPRI,*) 1298 & 'Remove .DIRECT or reformulate problem as DFT GGAKey HF=1' 1299 CALL QUIT('QFOCK not implemented for open shells') 1300 END IF 1301 CALL QFOCK(1,ISYM1,ISYM2,IGRSPI,ISPIN1,ISPIN2, 1302 & ZYMAT1,ZYMAT2,FC,CMO,FI,WRKO,KFREE,LFREE) 1303 DFTADD = DFTSAV 1304 ELSE 1305 CALL RSPFXD (FI,FA,QA,QB,H2AX, 1306 * IDAX,INTSYM, ISYMDN,DEN1,DEN2, 1307 * PVDEN,H2,H2X,WRKO,KFREE,LFREE, 1308 * LCON,LORB,IPRRSP,IGRSPI,ISPIN1,ISPIN2,ISPIN3, 1309 * IKLVL,ZYMAT1,ISYM1,ZYMAT2,ISYM2,ZYMAT3,ISYM3) 1310 END IF 1311 ELSE 1312 CALL RSPFX (FI,FA,QA,QB,H2AX, 1313 * IDAX,INTSYM, ISYMDN,DEN1,DEN2, 1314 * PVDEN,H2X,WRKO,KFREE,LFREE, 1315 * LCON,LORB,IPRRSP) 1316 END IF 1317C 1318 IF( IPRRSP .GT. 200 ) THEN 1319 WRITE(LUPRI,'(/A)') ' COMPLETED MATRICES IN RSP2GR' 1320 WRITE(LUPRI,'(A)') ' =============================' 1321 WRITE(LUPRI,'(///A)') ' Inactive fock matrix' 1322 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1323 WRITE(LUPRI,'(///A)') ' Active fock matrix' 1324 CALL OUTPUT(FA,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1325 WRITE(LUPRI,'(///A)') ' Qa matrix' 1326 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1327 WRITE(LUPRI,'(///A)') ' Qb matrix' 1328 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1329 WRITE(LUPRI,'(/A)') ' Active two-electron integrals' 1330 CALL PRIAC2(H2AX,NASHT,LUPRI) 1331 END IF 1332C 1333C Now distribute the computed matrices over the orbital part 1334C of the gradient: 1335C 1336 IF ( LORB ) THEN 1337 TRPLET = IGRSPI.NE.MULSP(ISPIN1,ISPIN2) 1338 CALL ORBEX(IGRSYM,ISYMDN,KZYVA,FI,FA,QA,QB,ETRS,OVLAP,DEN1, 1339 & MJWOP) 1340 END IF 1341C 1342C and over the configuration part: 1343C 1344 IF ( LCON ) THEN 1345 ISGN1 = (-1)**ISPIN1 1346 ISGN2 = (-1)**ISPIN2 1347 ISYMJ = MULD2H( IGRSYM, IREFSY ) 1348 NZCSTJ = MZCONF(IGRSYM) 1349 IF ( LREFST ) THEN 1350 LFREE = LWRKC - MZCONF(1) 1351 IF (LFREE.LT.0) CALL ERRWRK('RSP2GR',MZCONF(1),LWRKC) 1352 CALL GETREF(WRKC,MZCONF(1)) 1353 CALL CONEX(KZYVA,IGRSYM,FI,ONEFAC,H2AX,WRKC,MZCONF(1), 1354 * MZCONF(1),IREFSY,NZCSTJ,ISYMJ,LREFST,ETRS, 1355 * XINDX,ISPIN1,ISPIN2,WRKC(NCREF+1),LFREE) 1356 ELSE 1357 CALL CONEX(KZYVA,IGRSYM,FI,ONEFAC,H2AX,VEC, 1358 * NZYVEC,NZCVEC, 1359 * ISYMV, NZCSTJ,ISYMJ,LREFST,ETRS,XINDX, 1360 * ISPIN1,ISPIN2,WRKC,LWRKC) 1361 END IF 1362 END IF 1363C 1364C and that's it. 1365C 1366 CALL QEXIT('RSP2GR') 1367 RETURN 1368 END 1369C /* Deck rspfx */ 1370 SUBROUTINE RSPFX (FI,FA,QA,QB,H2AX, 1371 * IDAX,INTSYM, ISYMDN,DEN1,DEN2, 1372 * PVDEN,H2X,WRK,KFREE,LFREE, 1373 * LCON,LORB,IPRFX) 1374C 1375C 1376C 15-Nov-1991 hjaaj (pulled out of RSP2GR) 1377C 1378C This routine computes the FX matrices, that is, the Fock 1379C type matrices FI, FA, QA, and QB with transformed ("X") 1380C integrals. In addition, If LCON true then the H2AX matrix 1381C with transformed integrals is extracted for the CI routines. 1382C 1383#include "implicit.h" 1384C 1385 PARAMETER ( D0 = 0.0D0, DH = 0.50D0 ) 1386C 1387C Used from common blocks: 1388C ? 1389C 1390#include "maxorb.h" 1391#include "maxash.h" 1392#include "priunit.h" 1393#include "inforb.h" 1394#include "infind.h" 1395#include "infdim.h" 1396#include "infpri.h" 1397C 1398C 1399 DIMENSION FI(NORBT,NORBT),FA(NORBT,NORBT) 1400 DIMENSION QA(NORBT,NASHDI),QB(NORBT,NASHDI) 1401 DIMENSION H2X (NORBT,NORBT) 1402 DIMENSION H2AX(NASHDI,NASHDI,NASHDI,NASHDI) 1403 DIMENSION DEN1(NASHDI,NASHDI), DEN2(NASHDI,NASHDI,NASHDI,NASHDI) 1404 DIMENSION PVDEN(NASHDI,NASHDI) 1405 DIMENSION WRK(*) 1406 DIMENSION NEEDED(-4:6) 1407C 1408 LOGICAL LCON, LORB 1409C 1410 CALL QENTER('RSPFX') 1411C 1412C Put up the structure for reading the (transformed) integrals: 1413C IEND < 0 flags the completeness of the distributions 1414C 1415 NEEDED(-4:6) = 0 1416 NEEDED( 1:3) = 1 1417 IFCKSY = MULD2H(INTSYM,ISYMDN) 1418 IDIST = 0 1419 NEEDED(-4:6) = 0 1420 NEEDED( 1:3) = 1 14211000 CALL NXTH2X(IP,IQ,H2X,IDAX,NEEDED,WRK,KFREE,LFREE,IDIST) 1422 IF ( IDIST .LT. 0) GO TO 2000 1423C 1424 IPTYP = IOBTYP(IP) 1425 IQTYP = IOBTYP(IQ) 1426C 1427 IF (IPRFX .GT. 200) THEN 1428 WRITE(LUPRI,'(//A,/A,/A,I5,/A,I5,/A)') 1429 * ' Getting from NXTH2X: ', 1430 * ' ==================== ', 1431 * ' IP = ',IP, 1432 * ' IQ = ',IQ, 1433 * ' ====================' 1434 END IF 1435C 1436C Determine type of distribution 1437C 1438C 1 : inactive - inactive 1439C 2 : active - inactive 1440C 3 : inactive - active 1441C 4 : active - active 1442C 1443 IPQTYP = 0 1444 IF (( IPTYP .EQ. JTINAC) .AND. ( IQTYP .EQ. JTINAC)) IPQTYP = 1 1445 IF (( IPTYP .EQ. JTACT ) .AND. ( IQTYP .EQ. JTINAC)) IPQTYP = 2 1446 IF (( IPTYP .EQ. JTINAC) .AND. ( IQTYP .EQ. JTACT )) IPQTYP = 3 1447 IF (( IPTYP .EQ. JTACT ) .AND. ( IQTYP .EQ. JTACT )) IPQTYP = 4 1448C 1449 IF (IPQTYP .EQ. 0) THEN 1450 WRITE(LUPRI,'(//A,/A,I5,A,I5,/A)') 1451 * ' --> Warning : Unidentifyable distribution found ; ', 1452 * ' IPTYP = ', IPTYP,' IQTYP = ',IQTYP, 1453 * ' Continuing execution by skipping this load' 1454 GO TO 1000 1455 END IF 1456C 1457C Generally, an integral is labeled by (pq|rs). We know that p and q 1458C are occupied. When going into the specific routines, we substitute 1459C p,q,r,s for the special cases. We get symmetry characteristics: 1460C 1461 IPSYM = ISMO( IP ) 1462 IQSYM = ISMO( IQ ) 1463 IPQSYM = MULD2H(IPSYM,IQSYM) 1464C 1465 IF ( IPRFX .GT. 200) THEN 1466 WRITE(LUPRI,'(//A)') 'Characteristics of distribution:' 1467 WRITE(LUPRI,'(A)') '================================' 1468 WRITE(LUPRI,'(A,I5)') 'IPTYP = ', IPTYP 1469 WRITE(LUPRI,'(A,I5)') 'IQTYP = ', IQTYP 1470 WRITE(LUPRI,'(A,I5)') 'IPQTYP = ', IPQTYP 1471 WRITE(LUPRI,'(A,I5)') 'IPSYM = ', IPSYM 1472 WRITE(LUPRI,'(A,I5)') 'IQSYM = ', IQSYM 1473 WRITE(LUPRI,'(A,I5)') 'IPQSYM = ', IPQSYM 1474 WRITE(LUPRI,'(A)') '================================' 1475C 1476 WRITE(LUPRI,'(//A/A/)') 1477 * ' Two-electron integrals from this distribution', 1478 * ' =============================================' 1479 CALL OUTPUT(H2X,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1480 END IF 1481C 1482 IF (IPQTYP .LE. 3) THEN 1483 CALL FIXADD(INTSYM,IP,IQ,IPQTYP,IPQSYM,IPSYM,IQSYM, 1484 * FI,H2X ) 1485 END IF 1486 IF ((IPQTYP .NE. 1) .AND. LORB) THEN 1487 CALL FAXADD(IFCKSY,INTSYM,IP,IQ,IPQTYP,IPQSYM,IPSYM, 1488 * IQSYM,FA,H2X,DEN1 ) 1489 END IF 1490 IF ( (IPQTYP .EQ. 4) ) THEN 1491 ISWIP = ISW( IP ) - NISHT 1492 ISWIQ = ISW( IQ ) - NISHT 1493 IF (LORB) THEN 1494 CALL QXADD(IFCKSY,INTSYM,ISWIP,ISWIQ,IP,IQ,IPQTYP, 1495 * IPQSYM,IPSYM,IQSYM, 1496 * QA,QB,H2X,DEN2,PVDEN) 1497 END IF 1498C 1499C Store (xy|uv) in H2ACX for use in configuration part of lt. 1500C 1501 IF ( LCON ) THEN 1502 IF ( IPRFX .GT. 200) THEN 1503 WRITE(LUPRI,'(//A,/A,2(/A,I5),/A)' ) 1504 * 'Storing in H2AX :', 1505 * '=================', 1506 * 'ISWIP = ',ISWIP, 1507 * 'ISWIQ = ',ISWIQ, 1508 * '=================' 1509 END IF 1510 IUVSYM = MULD2H(IPQSYM, INTSYM) 1511 CALL DSH2AX(H2AX(1,1,ISWIP,ISWIQ),H2X,IUVSYM) 1512 END IF 1513C 1514 IF( IPRFX .GT. 200 ) THEN 1515 WRITE(LUPRI,'(/A)') ' PARTIAL MATRICES IN RSP2GR' 1516 WRITE(LUPRI,'(A)') ' ===========================' 1517 WRITE(LUPRI,'(//A)') ' Inactive fock matrix' 1518 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1519 WRITE(LUPRI,'(//A)') ' Active fock matrix' 1520 CALL OUTPUT(FA,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1521 WRITE(LUPRI,'(//A)') ' Qa matrix' 1522 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1523 WRITE(LUPRI,'(//A)') ' Qb matrix' 1524 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1525 END IF 1526C 1527 END IF 1528C 1529C We have now completed processing this load, so get the next 1530C 1531 GO TO 1000 1532C 15332000 CONTINUE 1534C 1535C Account for the fact that we calculated 2 times FA 1536 CALL DSCAL(NORBT*NORBT,DH,FA,1) 1537C 1538 CALL QEXIT('RSPFX') 1539 RETURN 1540 END 1541C /* Deck fixadd */ 1542 SUBROUTINE FIXADD(INTSYM,IP,IQ,IPQTYP,IPQSYM,IPSYM,IQSYM, 1543 * FI,H2X ) 1544C 1545C Add contribution presently in H2X to the inactive Fock matrix. 1546C Generic indices are i1,i2, for the rest we adhere to convention. 1547C Expression: 1548C 1549C iF(i1,i2) = iF(i1,i2) + sum(j) [ 2 (jj|i1 i2) - (i1 j|j i2) ] 1550C 1551C iF(i1,i2) = iF(i1,i2) + sum(j) [ 2 (jj|i1 i2) - (j i2|i1 j) ] 1552C 1553C Procedure: 1554C ---------- 1555C 1) Check if for the inactive - inactive distributions we have 1556C i = j. If so, add the direct contribution to the inactive Fock 1557C matrix. 1558C 1559C 2) For occupied - inactive distributions we have to subtract 1560C (oj|js) from iF(o,i2) for all s in the distribution 1561C 1562C 3) For secondary - inactive distributions we have to subtract 1563C (jo|aj) from iF(i1,o) 1564C 1565C Needed inactive Fock matrices and corresponding contributions: 1566C (First column exchange = case 2, second = case 3) 1567C---------------------------------------------------------------------- 1568C Label direct exchange Label direct exchange 1569C---------------------------------------------------------------------- 1570C ia (jj|ia) (ij|ja) ai (jj|ai) (aj|ji) = (ji|aj) 1571C ti (jj|ti) (tj|ji) at (jj|at) (aj|jt) = (jt|aj) 1572C it (jj|it) (ij|jt) 1573C ta (jj|ta) (tj|ja) 1574C tu (jj|tu) (tj|ju) 1575C---------------------------------------------------------------------- 1576C 1577#include "implicit.h" 1578C 1579 PARAMETER ( D2= 2.0D0, DM1 = -1.0D0 ) 1580C 1581#include "priunit.h" 1582#include "wrkrsp.h" 1583#include "rspprp.h" 1584#include "infhyp.h" 1585#include "infdim.h" 1586#include "inforb.h" 1587#include "infpri.h" 1588C 1589 DIMENSION FI(NORBT,NORBT),H2X (NORBT,NORBT) 1590C 1591 IF (IPQTYP .EQ. 4) THEN 1592C (We have no contributing distributions, but print a warning) 1593 WRITE(LUPRI,'(/A)') 1594 * ' --> Warning: invalid call of FIXADD, IPQTYP = 4' 1595 RETURN 1596 END IF 1597 IF (( IP .EQ. IQ) .AND. (IPQTYP .EQ. 1)) THEN 1598C Process case 1) inactive - inactive distributions coulomb part 1599C 1600C iF(i1,i2) = iF(i1,i2) + 2.0d0 (jj|i1 i2) 1601C 1602 CALL DAXPY(NORBT*NORBT,D2,H2X,1,FI,1) 1603 END IF 1604C 1605C Now process all exchange contributions 1606C 1607 IF (( IPQTYP .EQ. 1) .OR. (IPQTYP .EQ. 2) ) THEN 1608C Process case 2) : occupied - inactive distribution 1609C 1610C iF(o,i2) = iF(o,i2) - (oj|j i2) 1611C 1612 I2SYM = MULD2H(IPSYM,INTSYM) 1613 NORBI2 = NORB(I2SYM) 1614 IOFFI2 = IORB(I2SYM) + 1 1615 IF( NORBI2 .GT. 0) THEN 1616 CALL DAXPY(NORBI2,DM1,H2X(IQ,IOFFI2),NORBT, 1617 * FI(IP,IOFFI2),NORBT) 1618 END IF 1619 END IF 1620 IF (( IPQTYP .EQ. 1) .OR. (IPQTYP .EQ. 3) ) THEN 1621C Process case 3) : inactive - occupied distribution 1622C 1623C iF(a,o) = iF(a,o) - (jo|aj) 1624C 1625 IASYM = MULD2H(IQSYM,INTSYM) 1626 NSSHA = NSSH(IASYM) 1627 IOFFA = IORB(IASYM) + 1 + NOCC(IASYM) 1628 IF(NSSHA .GT. 0) THEN 1629 CALL DAXPY(NSSHA,DM1,H2X (IOFFA,IP),1, 1630 * FI(IOFFA,IQ),1) 1631 END IF 1632 END IF 1633 RETURN 1634 END 1635C /* Deck faxadd */ 1636 SUBROUTINE FAXADD(IFCKSY,INTSYM,IP,IQ,IPQTYP, 1637 * IPQSYM,IPSYM,IQSYM,FA,H2X,DEN1) 1638C 1639C Add contribution presently in H2X to the active Fock matrix. 1640C Note: twice the active Fock matrix is computed. 1641C 1642C aF(i1,i2) = aF(i1,i2) + [ (xy|i1 i2) - (i1 y|x i2) ] D(x,y) 1643C 1644C aF(i1,i2) = aF(i1,i2) + [ (xy|i1 i2) - (x i2|i1 y) ] D(x,y) 1645C 1646C Procedure: 1647C ---------- 1648C 1) Check for the active - active distributions and add the 1649C direct contribution to the inactive Fock matrix. 1650C 1651C 2) For occupied - active distributions we have to subtract 1652C (oy|xs) D(x,y) from aF(o,i2) for all s in the distribution 1653C 1654C 3) For active - occupied distributions we have to subtract 1655C (xo|ay) D(x,y) from aF(i1,o) 1656C 1657C Needed active Fock matrices and corresponding contributions: 1658C (First column exchange = case 2, second = case 3) 1659C---------------------------------------------------------------------- 1660C Label direct exchange Label direct exchange 1661C---------------------------------------------------------------------- 1662C ia (xy|ia) (iy|xa) ai (xy|ai) (ay|xi) = (xi|ay) 1663C ti (xy|ti) (ty|xi) 1664C it (xy|it) (iy|xt) 1665C---------------------------------------------------------------------- 1666C 1667#include "implicit.h" 1668C 1669 PARAMETER ( D2= 2.0D0 ) 1670C 1671#include "wrkrsp.h" 1672#include "rspprp.h" 1673#include "infhyp.h" 1674#include "maxorb.h" 1675#include "maxash.h" 1676#include "priunit.h" 1677#include "infdim.h" 1678#include "infind.h" 1679#include "inforb.h" 1680#include "infpri.h" 1681C 1682 DIMENSION FA(NORBT,NORBT),H2X (NORBT,NORBT) 1683 DIMENSION DEN1(NASHDI,NASHDI) 1684C 1685 IF (IPQTYP .EQ. 1) THEN 1686C (We have no contributing distributions, but print a warning) 1687 WRITE(LUPRI,'(/A)') 1688 * ' --> Warning: invalid call of FAXADD, IPQTYP = 1' 1689 RETURN 1690 END IF 1691 IF ( IPQTYP .EQ. 4 ) THEN 1692C Process case 1) active - active distributions coulomb part 1693C 1694C aF(i1,i2) = aF(i1,i2) + (xy|i1 i2) D(x,y) 1695C 1696 IXDEN = ISW(IP) - NISHT 1697 IYDEN = ISW(IQ) - NISHT 1698 FACT = D2 * DEN1(IXDEN,IYDEN) 1699 CALL DAXPY(NORBT*NORBT,FACT,H2X,1,FA,1) 1700 END IF 1701C 1702C Now process all exchange contributions 1703C 1704 IF (( IPQTYP .EQ. 3) .OR. (IPQTYP .EQ. 4) ) THEN 1705C Process case 2) : occupied - active distribution 1706C 1707C aF(o,i2) = aF(o,i2) - sum(x) D(x,y) (oy|x i2) 1708C 1709 I2SYM = MULD2H(IPSYM,IFCKSY) 1710 IXSYM = MULD2H( MULD2H( IPQSYM, I2SYM), INTSYM) 1711 NASHX = NASH(IXSYM) 1712 NORBI2 = NORB(I2SYM) 1713 IF ( (NORBI2.GT.0) .AND. (NASHX.GT.0) ) THEN 1714 IOFFI2 = IORB(I2SYM) + 1 1715 IOFFX = IORB(IXSYM) + 1 + NISH(IXSYM) 1716 IDENX = ISW( IOFFX ) - NISHT 1717 IDENY = ISW( IQ ) - NISHT 1718 CALL DGEMM('T','N',1,NORBI2,NASHX,-1.D0, 1719 & DEN1(IDENX,IDENY),NASHT, 1720 & H2X(IOFFX,IOFFI2),NORBT,1.D0, 1721 & FA(IP,IOFFI2),NORBT) 1722 END IF 1723 END IF 1724 IF ( ( IPQTYP .EQ. 2) .OR. ( IPQTYP .EQ. 4) ) THEN 1725C Process case 3) : active - occupied distribution 1726C 1727C aF(a,o) = aF(a,o) - sum(y) (xo|ay) D(x,y) 1728C 1729 IASYM = MULD2H(IQSYM,IFCKSY) 1730 IYSYM = MULD2H( MULD2H( IPQSYM, IASYM), INTSYM) 1731 NASHY = NASH(IYSYM) 1732 NSSHA = NSSH(IASYM) 1733 IF ( (NSSHA.GT.0) .AND. (NASHY.GT.0) ) THEN 1734 IOFFA = IORB(IASYM) + 1 + NOCC(IASYM) 1735 IOFFY = IORB(IYSYM) + 1 + NISH(IYSYM) 1736 IDENX = ISW( IP ) - NISHT 1737 IDENY = ISW( IOFFY ) - NISHT 1738 CALL DGEMM('N','T',NSSHA,1,NASHY,-1.D0, 1739 & H2X(IOFFA,IOFFY),NORBT, 1740 & DEN1(IDENX,IDENY),NASHT,1.D0, 1741 & FA(IOFFA,IQ),NORBT) 1742 END IF 1743 END IF 1744 RETURN 1745 END 1746C /* Deck qxadd */ 1747 SUBROUTINE QXADD(IFCKSY,INTSYM,ISWIP,ISWIQ,IP,IQ,IPQTYP, 1748 * IPQSYM,IPSYM,IQSYM,QA,QB, 1749 * H2X,DEN2,PVDEN) 1750C 1751C Process contributions from the Q matrices. Formulas: 1752C 1753C aQ(i1,y) = aQ(i1,y) + sum(w) (xv|w i1) d(xv;wy) 1754C 1755C bQ(i1,y) = bQ(i1,y) + sum(w) (xv|i1 w) d(xv;yw) 1756C 1757C The sum over x and v is taken by reading in all active contribu- 1758C tions. 1759C 1760#include "implicit.h" 1761#include "maxorb.h" 1762#include "maxash.h" 1763#include "priunit.h" 1764#include "inforb.h" 1765#include "infdim.h" 1766#include "infind.h" 1767#include "wrkrsp.h" 1768#include "infpri.h" 1769C 1770 DIMENSION QA(NORBT,NASHDI),QB(NORBT,NASHDI) 1771 DIMENSION DEN2(N2ASHX,N2ASHX), PVDEN(NASHDI,NASHDI) 1772 DIMENSION H2X(NORBT,NORBT) 1773C 1774 IF (IPQTYP .NE. 4) THEN 1775C (We have no contributing distributions, but print a warning) 1776 WRITE(LUPRI,'(/A)') 1777 * ' --> Warning: invalid call of QXADD, IPQTYP is not 4' 1778 RETURN 1779 END IF 1780C 1781C 1782C First get d(xv;**) 1783C 1784 CALL PVXDIS(ISWIP,ISWIQ,PVDEN,DEN2,4) 1785C 1786 IXSYM = IPSYM 1787 IVSYM = IQSYM 1788 DO 200 IWSYM = 1, NSYM 1789 II1SYM = MULD2H( MULD2H( IPQSYM, IWSYM) , INTSYM) 1790 IYSYM = MULD2H(II1SYM,IFCKSY) 1791 NORBI1 = NORB(II1SYM) 1792 NASHW = NASH(IWSYM) 1793 NASHY = NASH(IYSYM) 1794 IF((NORBI1.GT.0) .AND. (NASHW .GT. 0) .AND. 1795 * (NASHY .GT. 0) ) THEN 1796 IOFFI1 = IORB(II1SYM) + 1 1797 IOFFW = IORB(IWSYM) + NISH(IWSYM) + 1 1798 IOFFY = IORB(IYSYM) + NISH(IYSYM) + 1 1799 IDENW = ISW( IOFFW ) - NISHT 1800 IDENY = ISW( IOFFY ) - NISHT 1801C 1802C aQ(i1,y) = aQ(i1,y) + sum(w) (xv|w i1) d(xv;wy) 1803C 1804 CALL DGEMM('T','N',NORBI1,NASHY,NASHW,1.D0, 1805 & H2X(IOFFW,IOFFI1),NORBT, 1806 & PVDEN(IDENW,IDENY),NASHT,1.D0, 1807 & QA(IOFFI1,IDENY),NORBT) 1808C 1809C bQ(i1,y) = bQ(i1,y) + sum(w) (xv|i1 w) d(xv;yw) 1810C 1811 CALL DGEMM('N','T',NORBI1,NASHY,NASHW,1.D0, 1812 & H2X(IOFFI1,IOFFW),NORBT, 1813 & PVDEN(IDENY,IDENW),NASHT,1.D0, 1814 & QB(IOFFI1,IDENY),NORBT) 1815C 1816 END IF 1817200 CONTINUE 1818 RETURN 1819 END 1820C /* Deck orbex */ 1821 SUBROUTINE ORBEX(IGRSYM,ISYMDN,KZYVA, 1822 * FI,FA,QA,QB,ETRS,OVLAP,DEN1,MJWOP) 1823C 1824C Pupose: 1825C Get the gradient vector from the active ,inactive, Qa and Qb ma- 1826C trices. We have the expression: 1827C 1828C F(l,k) = <0| [ E(k,l), K ] |0> 1829C 1830C We have the expressions: 1831C 1832C 1) F(a,i) = 2 * OVLAP * iF(a,i) + 2 aF(a,i) 1833C 1834C 2) F(i,a) = - 2 * OVLAP * iF(i,a) - 2 aF(i,a) 1835C 1836C 3) F(t,i) = 2 * OVLAP * iF(t,i) + 2 aF(t,i) - sum(x) iF(x,i) D(t,i) 1837C - aQ(i,t) 1838C 1839C 4) F(i,t) = sum(x) iF(t,x) + bQ(i,t) - 2 * OVLAP * iF(i,t) 1840C - 2 aF(i,t) 1841C 1842C 5) F(t,a) = -sum(x) iF(x,a) D(x,t) - aQ(a,t) 1843C 1844C 6) F(a,t) = sum(x) iF(a,x) D(t,x) + bQ(a,t) 1845C 1846C 7) F(u,t) = sum(x) iF(u,x) D(t,x) + bQ(u,t) 1847C - sum(x) iF(x,t) D(x,u) - aQ(t,u) 1848C 1849#include "implicit.h" 1850C 1851 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0, DM1 = -1.0D0 ) 1852C 1853#include "maxorb.h" 1854#include "maxash.h" 1855#include "priunit.h" 1856#include "infvar.h" 1857#include "inforb.h" 1858#include "infind.h" 1859#include "infdim.h" 1860#include "infpri.h" 1861#include "infrsp.h" 1862#include "wrkrsp.h" 1863#include "qrinf.h" 1864#include "infspi.h" 1865C 1866 DIMENSION ETRS(KZYVA) 1867 DIMENSION FI(NORBT,NORBT),FA(NORBT,NORBT) 1868 DIMENSION QA(NORBT,NASHDI),QB(NORBT,NASHDI) 1869 DIMENSION DEN1(NASHDI,NASHDI), MJWOP(2,MAXWOP,8) 1870C 1871C 1872 IF( IPRRSP .GT. 200 ) THEN 1873 WRITE(LUPRI,'(/A)') ' Matrices in ORBEX' 1874 WRITE(LUPRI,'(A)') ' =================' 1875 WRITE(LUPRI,'(/A)') ' Density matrix in ORBEX' 1876 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,-1,LUPRI) 1877 WRITE(LUPRI,'(/A)') ' Inactive fock matrix in ORBEX' 1878 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1879 WRITE(LUPRI,'(/A)') ' Active fock matrix in ORBEX' 1880 CALL OUTPUT(FA,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1881 WRITE(LUPRI,'(/A)') ' Qa matrix in ORBEX' 1882 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 1883 WRITE(LUPRI,'(/A)') ' Qb matrix in ORBEX' 1884 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 1885 END IF 1886C 1887C 1888C Compute offset for the F(k,l) elements 1889C 1890 NZCONF = MZCONF(IGRSYM) 1891 NYCONF = MZCONF(IGRSYM) + MZVAR(IGRSYM) 1892C 1893C distribute fock matrices in the gradient: 1894C 1895 KSYM1 = 0 1896 DO 1300 IG = 1,MZWOPT(IGRSYM) 1897 K = MJWOP(1,IG,IGRSYM) 1898 L = MJWOP(2,IG,IGRSYM) 1899 KSYM = ISMO(K) 1900 LSYM = ISMO(L) 1901 IF( KSYM.NE.KSYM1 ) THEN 1902 KSYM1 = KSYM 1903 NORBK = NORB(KSYM) 1904 IORBK = IORB(KSYM) 1905 NASHK = NASH(KSYM) 1906 NISHK = NISH(KSYM) 1907 IASHK = IASH(KSYM) 1908 KXSYM = MULD2H(KSYM,ISYMDN) 1909 NORBKX= NORB(KXSYM) 1910 IORBKX= IORB(KXSYM) 1911 NASHKX= NASH(KXSYM) 1912 NISHKX= NISH(KXSYM) 1913 IASHKX= IASH(KXSYM) 1914 IORBL = IORB(LSYM) 1915 NASHL = NASH(LSYM) 1916 NISHL = NISH(LSYM) 1917 NORBL = NORB(LSYM) 1918 IASHL = IASH(LSYM) 1919 LXSYM = MULD2H(LSYM,ISYMDN) 1920 NORBLX= NORB(LXSYM) 1921 IORBLX= IORB(LXSYM) 1922 NASHLX= NASH(LXSYM) 1923 NISHLX= NISH(LXSYM) 1924 IASHLX= IASH(LXSYM) 1925 END IF 1926 NK = K - IORBK 1927 NL = L - IORBL 1928 ITYPK = IOBTYP(K) 1929 ITYPL = IOBTYP(L) 1930 IF ( ITYPK.EQ.JTINAC )THEN 1931C 1932C Do F(a,i) and F(i,a) and a part of F(t,i) and F(i,t) 1933C 1934 IF (.NOT.TRPLET) THEN 1935 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 1936 * + D2 * OVLAP * FI(L,K) 1937 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 1938 * - D2 * OVLAP * FI(K,L) 1939 END IF 1940 IF ( NASHT . GT . 0 ) THEN 1941 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) + D2 * FA(L,K) 1942 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) - D2 * FA(K,L) 1943 END IF 1944 IF ( ITYPL.EQ.JTACT ) THEN 1945C 1946C Complete F(t,i) and F(i,t) 1947C 1948 NWL = ISW(L) - NISHT 1949 IF (NASHLX .GT. 0) THEN 1950 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 1951 * - DDOT(NASHLX,FI(IORBLX+NISHLX+1,K),1, 1952 * DEN1(IASHLX+1,NWL),1) 1953 DO IX = 1,NASHLX 1954 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 1955 * + FI(K,IORBLX+NISHLX+IX)*DEN1(NWL,IASHLX+IX) 1956 END DO 1957 END IF 1958 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) - QA(K,NWL) 1959 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) + QB(K,NWL) 1960 END IF 1961 ELSE 1962C 1963C We now know that K labels an active orbital 1964C 1965 IF (ITYPL.EQ.JTACT) THEN 1966C 1967C Do F(t,u) 1968C 1969 NWL = ISW(L) - NISHT 1970 NWK = ISW(K) - NISHT 1971 IF (NASHLX .GT. 0) THEN 1972 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 1973 * -DDOT(NASHLX,FI(IORBLX+NISHLX+1,K),1, 1974 * DEN1(IASHLX+1,NWL),1) 1975 DO IX = 1,NASHLX 1976 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 1977 * + FI(K,IORBLX+NISHLX+IX)*DEN1(NWL,IASHLX+IX) 1978 END DO 1979 END IF 1980 IF (NASHKX .GT. 0) THEN 1981 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 1982 * -DDOT(NASHKX,FI(IORBKX+NISHKX+1,L),1, 1983 * DEN1(IASHKX+1,NWK),1) 1984 DO IX = 1,NASHKX 1985 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 1986 * + FI(L,IORBKX+NISHKX+IX)*DEN1(NWK,IASHKX+IX) 1987 END DO 1988 END IF 1989 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 1990 * +QB(L,NWK) - QA(K,NWL) 1991 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 1992 * +QB(K,NWL) - QA(L,NWK) 1993 ELSE 1994C 1995C Do F(a,t) and F(t,a) 1996C 1997 NWK = ISW(K) - NISHT 1998 IF (NASHKX .GT. 0) THEN 1999 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) 2000 * - DDOT(NASHKX,FI(IORBKX+NISHKX+1,L),1, 2001 * DEN1(IASHKX+1,NWK),1) 2002 DO IX = 1,NASHKX 2003 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) 2004 * + FI(L,IORBKX+NISHKX+IX)*DEN1(NWK,IASHKX+IX) 2005 END DO 2006 END IF 2007 ETRS(NYCONF+IG) = ETRS(NYCONF+IG) - QA(L,NWK) 2008 ETRS(NZCONF+IG) = ETRS(NZCONF+IG) + QB(L,NWK) 2009 ENDIF 2010 ENDIF 2011 1300 CONTINUE 2012C 2013 IF (IPRRSP .GT. 200 ) THEN 2014 WRITE(LUPRI,'(//A/A)') 2015 & ' ORBEX: Orbital part of linear transformation', 2016 & ' ============================================' 2017 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 2018 END IF 2019 RETURN 2020 END 2021C /* Deck dsh2ax */ 2022 SUBROUTINE DSH2AX(H2AX,H2X,IUVSYM) 2023C 2024C 891013-hjaaj 2025C Modified for no symmetry and name changed 901028 hh 2026C 2027C Add contributions from (xy) distribution to 2028C H2AC(uv,xy) from H2XY(uv); (uv) distribution has symmetry IUVSYM 2029C 2030C Subroutine can also be used for contributions to H2XAC. 2031C 2032#include "implicit.h" 2033C 2034C Used from common blocks: 2035C INFORB : NSYM, MULD2H, NORBT, NASHT, IASH, NASH 2036C INFIND : ISX, NSM, IROW 2037C 2038#include "maxash.h" 2039#include "maxorb.h" 2040#include "priunit.h" 2041#include "inforb.h" 2042#include "infind.h" 2043#include "infrsp.h" 2044#include "wrkrsp.h" 2045#include "infpri.h" 2046C 2047 DIMENSION H2AX(NASHT,NASHT),H2X(NORBT,NORBT) 2048C 2049 DO 200 NVW = 1,NASHT 2050 IVSYM = NSM(NVW) 2051 IUSYM = MULD2H(IVSYM, IUVSYM) 2052 IV = ISX(NISHT+NVW) 2053 NUWST = IASH(IUSYM) + 1 2054 NUWEND= IASH(IUSYM) + NASH(IUSYM) 2055 DO 100 NUW = NUWST,NUWEND 2056 IU = ISX(NISHT+NUW) 2057 H2AX(NUW,NVW) = H2AX(NUW,NVW) + H2X(IU,IV) 2058 100 CONTINUE 2059 200 CONTINUE 2060C 2061 IF( IPRRSP .GT. 1000 ) THEN 2062 WRITE(LUPRI,'(//A)') ' Two electron integrals in H2ACXY' 2063 CALL OUTPUT(H2X,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 2064 WRITE(LUPRI,'(//A)') ' Active distribution in H2AX ' 2065 CALL OUTPUT(H2AX,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 2066 END IF 2067C 2068 RETURN 2069 END 2070C /* Deck conex */ 2071 SUBROUTINE CONEX(KZYVA,IGRSYM,FI,ONEFAC,H2AX,VEC,NZYVEC, 2072 * NZCVEC,ISYMV,NZCSTJ,ISYMJ,LREFST,ETRS,XINDX, 2073 * ISPIN1,ISPIN2,WRK,LWRK) 2074C 2075C Purpose: 2076C Perform the configuration part of the linear transformation. 2077C Expressions: 2078C 2079C < j | K | 0(R) > 2080C - < 0(L) | K | j > 2081C 2082C 2083C The expression of the operator K between two arbitrary CSF's, in 2084C terms of iF where appropriate, is: 2085C 2086C <CSF1 | K | CSF2 > = E(inactive) + 2087C 2088C sum(xy) iF(x,y) D(x,y) + 0.5 sum(xyvw) (xy|vw) d(xy;vw) 2089C 2090#include "implicit.h" 2091C 2092 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0, D2 = 2.0D0 ) 2093C 2094#include "maxorb.h" 2095#include "maxash.h" 2096#include "priunit.h" 2097#include "wrkrsp.h" 2098#include "infrsp.h" 2099#include "infopt.h" 2100#include "inforb.h" 2101#include "infind.h" 2102#include "infpri.h" 2103#include "rspprp.h" 2104#include "infhyp.h" 2105#include "infvar.h" 2106#include "qrinf.h" 2107#include "cbgetdis.h" 2108#include "infspi.h" 2109#include "infdis.h" 2110C 2111 DIMENSION FI(NORBT,NORBT) 2112 DIMENSION VEC(NZYVEC) 2113 DIMENSION WRK(*) 2114 DIMENSION H2AX(N2ASHX,N2ASHX) 2115 DIMENSION ETRS(KZYVA) 2116 DIMENSION XINDX(*) 2117C 2118 LOGICAL NOH2, IH8SM, LREFST,ITEST 2119 DATA ITEST /.FALSE./ 2120C 2121C Allocate work space for copying inactive Fock matrix and keeping 2122C the CI part of the vector, initialise if necessary 2123C 2124 KFIAC = 1 2125 KTRS = KFIAC + N2ASHX 2126 KFREE = KTRS + MZCONF(IGRSYM) 2127 LFREE = LWRK - KFREE 2128 IF (LFREE.LT.0) CALL ERRWRK('CONEX',KFREE-1,LWRK) 2129C 2130 IADH2 = -1 2131C ... IADH2 .lt. 0: H2AC in core 2132 NOH2 = .FALSE. 2133 IH8SM = .FALSE. 2134 IF ( LREFST ) THEN 2135 ISYMCI = IREFSY 2136 ELSE 2137 ISYMCI = MULD2H(IREFSY,ISYMV) 2138 END IF 2139C 2140C 2141 IF( IPRRSP .GT. 200) THEN 2142 WRITE(LUPRI,'(/A/)') ' Inactive Fock matrix in CONEX' 2143 CALL OUTPUT(FI,1,NORBT,1,NORBT,NORBT,NORBT, 1,LUPRI) 2144 WRITE(LUPRI,'(/A/)') ' Vector in CONEX' 2145 CALL OUTPUT(VEC,1,NZYVEC,1,1,NZYVEC,1, 1,LUPRI) 2146 END IF 2147C 2148C Do < j | K | 0(R) > 2149C ========================== 2150C 2151C Extract inactive contribution from FI in EINFI 2152C Copy transpose active-active blocks of inactive Fock matrix to WRK 2153C 2154 CALL DZERO(WRK(KTRS),MZCONF(IGRSYM)) 2155 EINFI = ONEFAC 2156 IF (ISPIN1.EQ.ISPIN2) THEN 2157 DO 105 IW = 1, NISHT 2158 IX = ISX(IW) 2159 EINFI = EINFI + FI(IX,IX) 2160 105 CONTINUE 2161 END IF 2162C 2163 DO 110 IW = 1, NASHT 2164 IX = ISX(NISHT + IW) 2165 DO 111 JW = 1, NASHT 2166 JX = ISX(NISHT + JW) 2167 IND = (IW-1) * NASHT + JW + KFIAC - 1 2168 WRK(IND) = FI(IX,JX) 2169111 CONTINUE 2170110 CONTINUE 2171C 2172 IF( IPRRSP .GT. 200) THEN 2173 WRITE(LUPRI,'(/A,F20.10)') 2174 * ' Inactive contribution from FI in CONEX',EINFI 2175 WRITE(LUPRI,'(/A/)') ' Active parts of FI in CONEX' 2176 CALL OUTPUT(WRK(KFIAC),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 2177 END IF 2178C 2179C Tell CISIGD to use the transposed H2AX array 2180 IF (DIROIT) THEN 2181 DISTYP = INFDIS(1) 2182 ELSE 2183 DISTYP = 9 2184 END IF 2185C 2186 CALL CISIGD(ISYMCI,ISYMJ,NZCVEC,NZCSTJ,VEC, 2187 * WRK(KTRS),WRK(KFIAC),H2AX, 2188 * NOH2,IH8SM,XINDX,ISPIN1,ISPIN2,WRK(KFREE), 2189 * LFREE) 2190C 2191C Add inactive contribution 2192C 2193 IF (ISYMCI .EQ. ISYMJ .AND. ISPIN1.EQ.ISPIN2) THEN 2194 IF ( NISHT .GT. 0 ) 2195 * CALL DAXPY(NZCVEC,EINFI,VEC,1,WRK(KTRS),1) 2196 END IF 2197C 2198C 2199C If we are testing, control for the CREF component 2200C 2201 IF( ITEST .AND. ISPIN1 .EQ. 0 .AND. ISPIN2 .EQ. 0 ) THEN 2202 IF (LREFST .AND. ISYMCI.EQ.ISYMJ) THEN 2203 T1 = DDOT(NZCVEC,VEC,1,WRK(KTRS),1) 2204 IF( IPRRSP .GT. 200) THEN 2205 WRITE(LUPRI,*) ' CREF factor in CONEX, T1 =',T1 2206 WRITE(LUPRI,'(/A/)') 2207 * ' Z result in CONEX before removing CREF component' 2208 CALL OUTPUT(WRK(KTRS),1,NZCVEC,1,1,NZCVEC,1,1,LUPRI) 2209 END IF 2210 CALL DAXPY(NZCVEC,(-T1),VEC,1,WRK(KTRS),1) 2211 END IF 2212 END IF 2213C 2214 IF( IPRRSP .GT. 200) THEN 2215 WRITE(LUPRI,'(/A/)') 'Final Z result in CONEX' 2216 CALL OUTPUT(WRK(KTRS),1,NZCVEC,1,1,NZCVEC,1,1,LUPRI) 2217 END IF 2218C 2219C Take care of minus sign in CR if we have anything but the 2220C reference state 2221C 2222 IF (LREFST) THEN 2223 FACT = D1 2224 ELSE 2225 FACT = DM1 2226 END IF 2227 CALL DAXPY(MZCONF(IGRSYM),FACT,WRK(KTRS),1,ETRS,1) 2228C 2229 IF (IPRRSP .GT. 200 ) THEN 2230 WRITE(LUPRI,'(/A,/A,/)') 2231 * ' Y component of configuration part of linear', 2232 * ' transformation added. Gradient is now:' 2233 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 2234 END IF 2235C 2236C Do - < 0(L) | K | j > 2237C ========================== 2238C 2239 CALL DZERO(WRK(KTRS),MZCONF(IGRSYM)) 2240C Tell CISIGD to use the H2AX array directly 2241 IF (DIROIT) THEN 2242 DISTYP = INFDIS(2) 2243 ELSE 2244 DISTYP = 10 2245 END IF 2246C 2247 IF( LREFST ) THEN 2248 ICOFF = 1 2249 ELSE 2250 ICOFF = MZVAR(ISYMV) + 1 2251 END IF 2252 IEOFF = MZVAR(IGRSYM) + 1 2253C 2254C Copy active-active blocks of inactive Fock matrix to WRK 2255C 2256 DO 210 IW = 1, NASHT 2257 IX = ISX(NISHT + IW) 2258 DO 211 JW = 1, NASHT 2259 JX = ISX(NISHT + JW) 2260 IND = (IW-1) * NASHT + JW + KFIAC - 1 2261 WRK(IND) = FI(JX,IX) 2262211 CONTINUE 2263210 CONTINUE 2264C 2265 IF( IPRRSP .GT. 200) THEN 2266 WRITE(LUPRI,'(/A,F20.10)') 2267 * ' Inactive contribution from FI in CONEX',EINFI 2268 WRITE(LUPRI,'(/A/)') 2269 * ' Transposed active parts of FI in CONEX' 2270 CALL OUTPUT(WRK(KFIAC),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 2271 END IF 2272C 2273 CALL CISIGD(ISYMCI,ISYMJ,NZCVEC,NZCSTJ,VEC(ICOFF), 2274 * WRK(KTRS),WRK(KFIAC),H2AX, 2275 * NOH2,IH8SM,XINDX,ISPIN1,ISPIN2,WRK(KFREE), 2276 * LFREE) 2277C 2278C 2279C Add inactive contribution 2280C 2281 IF (ISYMCI .EQ. ISYMJ .AND. ISPIN1.EQ.ISPIN2) THEN 2282 IF ( NISHT .GT. 0 ) 2283 * CALL DAXPY(NZCVEC,EINFI,VEC(ICOFF), 2284 * 1,WRK(KTRS),1) 2285 END IF 2286C 2287 IF( ITEST .AND. ISPIN1 .EQ. 0 .AND. ISPIN2 .EQ. 0 ) THEN 2288C 2289 IF (LREFST .AND. ISYMCI.EQ.ISYMJ) THEN 2290 T1 = DDOT(NZCVEC,VEC(ICOFF),1,WRK(KTRS),1) 2291 IF( IPRRSP .GT. 200) THEN 2292 WRITE(LUPRI,*) ' CREF factor in CONEX, T1 =',T1 2293 WRITE(LUPRI,'(/A/)') 2294 * ' Y result in CONEX before removing CREF component' 2295 CALL OUTPUT(WRK(KTRS),1,NZCVEC,1,1,NZCVEC,1,1,LUPRI) 2296 END IF 2297 CALL DAXPY(NZCVEC,(-T1),VEC(ICOFF),1,WRK(KTRS),1) 2298 END IF 2299 END IF 2300C 2301 IF( IPRRSP .GT. 200) THEN 2302 WRITE(LUPRI,'(/A/)') ' Final result in CONEX' 2303 CALL OUTPUT(WRK(KTRS),1,NZCVEC,1,1,NZCVEC,1,1,LUPRI) 2304 END IF 2305C 2306C Do the minus sign of the whole term 2307C 2308 CALL DAXPY(MZCONF(IGRSYM),DM1,WRK(KTRS),1,ETRS(IEOFF),1) 2309C 2310 IF (IPRRSP .GT. 200 ) THEN 2311 WRITE(LUPRI,'(/A,/A,/)') 2312 * ' Z component of configuration part of linear', 2313 * ' transformation added. Gradient is now:' 2314 CALL OUTPUT(ETRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 2315 END IF 2316C 2317 RETURN 2318 END 2319C /* Deck s3init */ 2320 SUBROUTINE S3INIT(KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC, 2321 * IGRSPI,ISPIN1,ISPIN2, 2322 * VECB,VECC,VECA,ATEST, 2323 * S3TRS,XINDX,UDV,MJWOP,WRK,LWRK) 2324C 2325C Layout the core for the calculation of S3 times two vectors 2326C 2327#include "implicit.h" 2328#include "infdim.h" 2329#include "inforb.h" 2330#include "maxash.h" 2331#include "infvar.h" 2332#include "infrsp.h" 2333#include "wrkrsp.h" 2334#include "rspprp.h" 2335#include "infhyp.h" 2336#include "qrinf.h" 2337#include "infpri.h" 2338C 2339 LOGICAL ATEST 2340 DIMENSION WRK(*) 2341 DIMENSION S3TRS(KZYVA) 2342 DIMENSION VECB(KZYVB) 2343 DIMENSION VECC(KZYVC) 2344 DIMENSION VECA(KZYVA) 2345 DIMENSION XINDX(*) 2346 DIMENSION UDV(NASHDI,NASHDI), MJWOP(2,MAXWOP,8) 2347C 2348C Initialise the gradient to zero. 2349C 2350 CALL QENTER('S3INIT') 2351 CALL DZERO(S3TRS,KZYVA) 2352C 2353 KZYMAT = 1 2354 KDEN1 = KZYMAT + NORBT * NORBT 2355 KFREE = KDEN1 + NASHT * NASHT 2356 LWRKF = LWRK - KFREE 2357 IF (LWRKF.LT.0) CALL ERRWRK('S3INIT',KFREE-1,LWRK) 2358 2359 CALL S3DRV(KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC, 2360 * IGRSPI,ISPIN1,ISPIN2, 2361 * S3TRS,VECB,VECC,VECA,ATEST, 2362 * UDV,WRK(KZYMAT),WRK(KDEN1),XINDX,MJWOP,WRK(KFREE), 2363 * LWRKF) 2364C 2365 CALL QEXIT('S3INIT') 2366 RETURN 2367 END 2368C /* Deck s3drv */ 2369 SUBROUTINE S3DRV(KZYVA,KZYVB,KZYVC,ISYMA,ISYMB,ISYMC, 2370 * IGRSPI,ISPIN1,ISPIN2, 2371 * S3TRS,VECB,VECC,VECA,ATEST, 2372 * UDV,ZYMAT,DEN1,XINDX,MJWOP,WRK,LWRK) 2373C 2374C Drive the computation of S[3] times two vectors 2375C 2376#include "implicit.h" 2377#include "priunit.h" 2378#include "infdim.h" 2379#include "inforb.h" 2380#include "maxorb.h" 2381#include "maxash.h" 2382#include "infrsp.h" 2383#include "wrkrsp.h" 2384#include "rspprp.h" 2385#include "infhyp.h" 2386#include "infvar.h" 2387#include "infind.h" 2388#include "qrinf.h" 2389#include "infpri.h" 2390#include "infspi.h" 2391C 2392 PARAMETER(DM1=-1.0D0, DMH=-0.5D0, D0=0.0D0, D1=1.0D0, D2=2.0D0 ) 2393C 2394 DIMENSION WRK(*) 2395 DIMENSION S3TRS(KZYVA) 2396 DIMENSION VECB(KZYVB) 2397 DIMENSION VECC(KZYVC) 2398 DIMENSION VECA(KZYVA) 2399 DIMENSION DEN1(NASHDI,NASHDI) 2400 DIMENSION ZYMAT(NORBT,NORBT) 2401 DIMENSION XINDX(*) 2402 DIMENSION UDV(NASHDI,NASHDI), MJWOP(2,MAXWOP,8) 2403C 2404 LOGICAL ATEST, NORHO2, TDM, LORB, LCON, LREFST, LREF 2405C 2406 CALL QENTER('S3DRV') 2407 CALL DZERO(S3TRS, KZYVA) 2408C 2409C Layout some workspace 2410C 2411 KCREF = 1 2412 KFREE = KCREF + MZCONF(1) 2413 LFREE = LWRK - KFREE + 1 2414 IF (LFREE.LT.0) CALL ERRWRK('S3DRV 1',KFREE-1,LWRK) 2415C 2416 NSIM = 1 2417 IGRSYM = ISYMA 2418 TDM = .TRUE. 2419 NORHO2 = .TRUE. 2420 S3TMZC = D0 2421 S3TMZO = D0 2422 S3TMYC = D0 2423 S3TMYO = D0 2424 S3TERM = D0 2425 KZCA = MZCONF(ISYMA) 2426 KZOA = MZWOPT(ISYMA) 2427 KAZC = 1 2428 KAZO = KAZC + KZCA 2429 KAYC = KAZO + KZOA 2430 KAYO = KAYC + KZCA 2431C 2432C Case 1 2433C ====== 2434 IF ((MZWOPT(ISYMB) .EQ. 0).OR.(MZWOPT(ISYMC).EQ. 0)) GO TO 2000 2435C Transform the kappa's 2436C 2437 CALL TRZYMT(NSIM,VECC,VECB,KZYVC,KZYVB,ISYMC,ISYMB,ZYMAT,MJWOP, 2438 * WRK,LWRK) 2439C 2440C Make copies of the MC density matrix in DEN1 2441C 2442 CALL DCOPY(NASHT*NASHT,UDV,1,DEN1,1) 2443 OVLAP = D1 2444 ISYMDN = 1 2445C 2446C Put the factor of minus one half into the operator matrix 2447C 2448 CALL DSCAL(NORBT*NORBT,DMH,ZYMAT,1) 2449C 2450C Make the gradient 2451C 2452 CALL GETREF(WRK(KCREF),MZCONF(1)) 2453C 2454 INTSYM = MULD2H(ISYMB,ISYMC) 2455 INTSPI = MULSP(ISPIN1,ISPIN2) 2456 IF ( INTSYM .NE. IGRSYM ) THEN 2457 WRITE(LUPRI,'(/A,I5)') ' INTSYM = ', INTSYM 2458 WRITE(LUPRI,'(/A,I5)') ' IGRSYM = ', IGRSYM 2459 CALL QUIT(' INTSYM is not equal to IGRSYM in S3DRV') 2460 ENDIF 2461 IF ( INTSPI .NE. IGRSPI ) THEN 2462 WRITE(LUPRI,'(/A,I5)') ' INTSPI = ', INTSPI 2463 WRITE(LUPRI,'(/A,I5)') ' IGRSPI = ', IGRSPI 2464 CALL QUIT(' INTSPI is not equal to IGRSPI in S3DRV') 2465 END IF 2466 ISYMV = IREFSY 2467 ISYMST = MULD2H(IGRSYM,IREFSY) 2468 IF ( ISYMST .EQ. IREFSY ) THEN 2469 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 2470 ELSE 2471 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 2472 END IF 2473 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 2474 LREFST = .TRUE. 2475 NZYVEC = MZCONF(1) 2476 NZCVEC = MZCONF(1) 2477 CALL RSP1GR(NSIM,KZYVA,INTSYM,INTSPI,IGRSYM,IGRSPI,ISYMV,S3TRS, 2478 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,ZYMAT, 2479 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREFST) 2480C 2481 IF (ATEST) THEN 2482 S3NWZC = DDOT(KZCA,VECA(KAZC),1,S3TRS(KAZC),1) 2483 S3NWZO = DDOT(KZOA,VECA(KAZO),1,S3TRS(KAZO),1) 2484 S3NWYC = DDOT(KZCA,VECA(KAYC),1,S3TRS(KAYC),1) 2485 S3NWYO = DDOT(KZOA,VECA(KAYO),1,S3TRS(KAYO),1) 2486 S3NEW = S3NWZC + S3NWZO + S3NWYC + S3NWYO 2487 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 1 :',S3NEW-S3TERM 2488 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 1 ZC :',S3NWZC-S3TMZC 2489 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 1 ZO :',S3NWZO-S3TMZO 2490 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 1 YC :',S3NWYC-S3TMYC 2491 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 1 YO :',S3NWYO-S3TMYO 2492 S3TMZC = S3NWZC 2493 S3TMZO = S3NWZO 2494 S3TMYC = S3NWYC 2495 S3TMYO = S3NWYO 2496 S3TERM = S3NEW 2497 END IF 2498C 2499 IF (IPRRSP .GT. 100 ) THEN 2500 WRITE(LUPRI,'(A)') ' Case 1 for S3 term' 2501 CALL OUTPUT(S3TRS,1,KZYVA,1,1,KZYVA,1,1,LUPRI) 2502 END IF 2503C 2504C Case 2 2505C ====== 2506 2000 IF((MZWOPT(ISYMB).EQ.0).OR.(MZCONF(ISYMC).EQ.0)) GO TO 3000 2507C 2508C Construct the density matrix <02L|..|0> + <0|..|02R> 2509C 2510 OVLAP = D0 2511 CALL GETREF(WRK(KCREF),MZCONF(1)) 2512C 2513 CALL DZERO(DEN1,NASHT*NASHT) 2514C 2515 ILSYM = IREFSY 2516 IRSYM = MULD2H(ISYMC,IREFSY) 2517 NCL = MZCONF(1) 2518 NCR = MZCONF(ISYMC) 2519 KZVARL = MZCONF(1) 2520 KZVARR = MZYVAR(ISYMC) 2521 LREF = .TRUE. 2522 ISYMDN = MULD2H(ILSYM,IRSYM) 2523 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 2524 * WRK(KCREF),VECC,OVLAP,DEN1,DUMMY,IGRSPI,ISPIN1,TDM, 2525 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 2526C 2527C Construct the one electron matrix 2528C 2529 CALL GTZYMT(NSIM,VECB,KZYVB,ISYMB,ZYMAT,MJWOP) 2530C 2531C Put the factor of minus one into the operator matrix 2532C 2533 CALL DSCAL(NORBT*NORBT,DM1,ZYMAT,1) 2534C 2535C Make the gradient 2536C 2537 INTSYM = ISYMB 2538 ISYMV = ISYMC 2539 ISYMST = MULD2H(IGRSYM,IREFSY) 2540 IF ( ISYMST .EQ. IREFSY ) THEN 2541 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 2542 ELSE 2543 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 2544 END IF 2545 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 2546 LREFST = .FALSE. 2547 NZCVEC = MZCONF(ISYMC) 2548 NZYVEC = MZYVAR(ISYMC) 2549 CALL RSP1GR(NSIM,KZYVA,INTSYM,ISPIN1,IGRSYM,IGRSPI,ISYMV,S3TRS, 2550 * VECC,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,ZYMAT, 2551 * XINDX,MJWOP,WRK,LWRK,LORB,LCON,LREFST) 2552C 2553C 2554 IF (ATEST) THEN 2555 S3NWZC = DDOT(KZCA,VECA(KAZC),1,S3TRS(KAZC),1) 2556 S3NWZO = DDOT(KZOA,VECA(KAZO),1,S3TRS(KAZO),1) 2557 S3NWYC = DDOT(KZCA,VECA(KAYC),1,S3TRS(KAYC),1) 2558 S3NWYO = DDOT(KZOA,VECA(KAYO),1,S3TRS(KAYO),1) 2559 S3NEW = S3NWZC + S3NWZO + S3NWYC + S3NWYO 2560 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 2 :',S3NEW-S3TERM 2561 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 2 ZC :',S3NWZC-S3TMZC 2562 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 2 ZO :',S3NWZO-S3TMZO 2563 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 2 YC :',S3NWYC-S3TMYC 2564 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 2 YO :',S3NWYO-S3TMYO 2565 S3TMZC = S3NWZC 2566 S3TMZO = S3NWZO 2567 S3TMYC = S3NWYC 2568 S3TMYO = S3NWYO 2569 S3TERM = S3NEW 2570 END IF 2571 IF (IPRRSP .GT. 100 ) THEN 2572 WRITE(LUPRI,'(A)') ' Case 2 for S3 term' 2573 CALL OUTPUT(S3TRS,1,KZYVA,1,NSIM,KZYVA,NSIM,1,LUPRI) 2574 END IF 2575C 2576C Case 3 2577C ====== 2578 3000 IF((MZCONF(ISYMB).EQ.0).OR.(MZCONF(ISYMC).EQ.0)) GO TO 4000 2579C 2580C Construct <01L|..|02R> + <02L|..|01R> 2581C 2582 OVLAP = D0 2583 CALL DZERO(DEN1,NASHT*NASHT) 2584C 2585 ILSYM = MULD2H(IREFSY,ISYMB) 2586 IRSYM = MULD2H(IREFSY,ISYMC) 2587 NCL = MZCONF(ISYMB) 2588 NCR = MZCONF(ISYMC) 2589 KZVARL = MZYVAR(ISYMB) 2590 KZVARR = MZYVAR(ISYMC) 2591 LREF = .FALSE. 2592 ISYMDN = MULD2H(ILSYM,IRSYM) 2593 IF (IGRSYM .EQ. ISYMDN) THEN 2594 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 2595 * VECB,VECC,OVLAP,DEN1,DUMMY,IGRSPI,0,TDM, 2596 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 2597C 2598C Make the gradient on the fly: 2599C 2600 NZCONF = MZCONF(IGRSYM) 2601 NYCONF = MZCONF(IGRSYM) + MZVAR(IGRSYM) 2602C 2603 DO 500 IC = 1, MZWOPT(IGRSYM) 2604 K = MJWOP(1,IC,IGRSYM) 2605 L = MJWOP(2,IC,IGRSYM) 2606 ITYPK = IOBTYP(K) 2607 ITYPL = IOBTYP(L) 2608 IF(ITYPK.EQ.JTACT .AND. ITYPL. EQ. JTACT) THEN 2609 NWK = ISW(K) - NISHT 2610 NWL = ISW(L) - NISHT 2611 S3TRS(NZCONF+IC) = S3TRS(NZCONF+IC) 2612 * + DEN1(NWK,NWL) 2613 S3TRS(NYCONF+IC) = S3TRS(NYCONF+IC) 2614 * + DEN1(NWL,NWK) 2615 END IF 2616500 CONTINUE 2617 END IF 2618C 2619 IF (ATEST) THEN 2620 S3NWZC = DDOT(KZCA,VECA(KAZC),1,S3TRS(KAZC),1) 2621 S3NWZO = DDOT(KZOA,VECA(KAZO),1,S3TRS(KAZO),1) 2622 S3NWYC = DDOT(KZCA,VECA(KAYC),1,S3TRS(KAYC),1) 2623 S3NWYO = DDOT(KZOA,VECA(KAYO),1,S3TRS(KAYO),1) 2624 S3NEW = S3NWZC + S3NWZO + S3NWYC + S3NWYO 2625 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 3 :',S3NEW-S3TERM 2626 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 3 ZC :',S3NWZC-S3TMZC 2627 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 3 ZO :',S3NWZO-S3TMZO 2628 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 3 YC :',S3NWYC-S3TMYC 2629 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 3 YO :',S3NWYO-S3TMYO 2630 S3TMZC = S3NWZC 2631 S3TMZO = S3NWZO 2632 S3TMYC = S3NWYC 2633 S3TMYO = S3NWYO 2634 S3TERM = S3NEW 2635 END IF 2636C 2637 IF (IPRRSP .GT. 100 ) THEN 2638 WRITE(LUPRI,'(A)') ' Case 3 for S3 term' 2639 CALL OUTPUT(S3TRS,1,KZYVA,1,NSIM,KZYVA,NSIM,1,LUPRI) 2640 END IF 2641C 2642C Case 4 2643C ====== 2644 4000 IF((MZWOPT(ISYMB).EQ.0).OR.(MZCONF(ISYMC).EQ.0)) GO TO 5000 2645C Do Sj(2) <0 |-kappa(1)| 0> 2646C 2647 IF (ISYMB .EQ. 1) THEN 2648C 2649C Compute the expectation value 2650C 2651 KSYMOP = ISYMB 2652 CALL GTZYMT(NSIM,VECB,KZYVB,ISYMB,ZYMAT,MJWOP) 2653C 2654 IPRONE = 75 2655 TRPLET = ISPIN1.EQ.1 2656 CALL PRPONE(ZYMAT,UDV,FACT,IPRONE,'<0 |-kappa(1)| 0>') 2657C 2658 IF (FACT .NE. D0) THEN 2659 NZCONF = MZCONF(IGRSYM) 2660 NZVAR = MZVAR(IGRSYM) 2661 CALL DAXPY(NZCONF,-FACT,VECC,1,S3TRS,1) 2662 CALL DAXPY(NZCONF,-FACT,VECC(NZVAR+1),1, 2663 * S3TRS(NZVAR+1),1) 2664 END IF 2665 END IF 2666C 2667 IF (ATEST) THEN 2668 S3NWZC = DDOT(KZCA,VECA(KAZC),1,S3TRS(KAZC),1) 2669 S3NWZO = DDOT(KZOA,VECA(KAZO),1,S3TRS(KAZO),1) 2670 S3NWYC = DDOT(KZCA,VECA(KAYC),1,S3TRS(KAYC),1) 2671 S3NWYO = DDOT(KZOA,VECA(KAYO),1,S3TRS(KAYO),1) 2672 S3NEW = S3NWZC + S3NWZO + S3NWYC + S3NWYO 2673 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 4 :',S3NEW-S3TERM 2674 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 4 ZC :',S3NWZC-S3TMZC 2675 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 4 ZO :',S3NWZO-S3TMZO 2676 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 4 YC :',S3NWYC-S3TMYC 2677 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 4 YO :',S3NWYO-S3TMYO 2678 S3TMZC = S3NWZC 2679 S3TMZO = S3NWZO 2680 S3TMYC = S3NWYC 2681 S3TMYO = S3NWYO 2682 S3TERM = S3NEW 2683 END IF 2684C 2685 IF (IPRRSP .GT. 100 ) THEN 2686 WRITE(LUPRI,'(A)') ' Case 4 for S3 term' 2687 CALL OUTPUT(S3TRS,1,KZYVA,1,NSIM,KZYVA,NSIM,1,LUPRI) 2688 END IF 2689C 2690C Case 5: 2691C ======= 2692 5000 IF((MZCONF(ISYMB).EQ.0).OR.(MZCONF(ISYMC).EQ.0)) GO TO 6000 2693C 2694 IF ( (ISYMB .EQ. ISYMC) .AND. (NASHT .GT. 0) ) THEN 2695 FACT = D0 2696 ICOFF = MZVAR(ISYMC) 2697 DO 730 IC = 1, MZCONF(ISYMC) 2698 FACT = FACT + VECB(IC) * VECC(ICOFF + IC) 2699 FACT = FACT + VECB(ICOFF + IC) * VECC(IC) 2700 730 CONTINUE 2701 IF (IPRRSP .GT. 100) THEN 2702 WRITE(LUPRI,'(/A,F20.10)') ' Factor in case 5 = ',FACT 2703 END IF 2704C 2705 NZCONF = MZCONF(IGRSYM) 2706 NYCONF = MZCONF(IGRSYM) + MZVAR(IGRSYM) 2707C 2708 IF (IGRSPI.EQ.1) THEN 2709 KCREF = 1 2710 KFREE = KCREF + MZCONF(1) 2711 LFREE = LWRK - KFREE + 1 2712 IF (LFREE.LT.0) CALL ERRWRK('S3DRV 2',KFREE-1,LWRK) 2713 CALL GETREF(WRK(KCREF),MZCONF(1)) 2714 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 2715 * WRK(KCREF),WRK(KCREF), DEN1,RHO2, 2716 * IGRSPI,0,.FALSE.,.TRUE.,XINDX,WRK, KFREE,LFREE) 2717 ELSE 2718 CALL DCOPY(N2ASHX,UDV,1,DEN1,1) 2719 END IF 2720C 2721 DO 710 IC = 1, MZWOPT(IGRSYM) 2722 K = MJWOP(1,IC,IGRSYM) 2723 L = MJWOP(2,IC,IGRSYM) 2724 ITYPK = IOBTYP(K) 2725 ITYPL = IOBTYP(L) 2726 IF(ITYPK.EQ.JTACT .AND. ITYPL. EQ. JTACT) THEN 2727 NWK = ISW(K) - NISHT 2728 NWL = ISW(L) - NISHT 2729 S3TRS(NZCONF+IC) = S3TRS(NZCONF+IC) 2730 * + FACT * DEN1(NWK,NWL) 2731 S3TRS(NYCONF+IC) = S3TRS(NYCONF+IC) 2732 * + FACT * DEN1(NWL,NWK) 2733 END IF 2734710 CONTINUE 2735 END IF 2736C 2737 IF (ATEST) THEN 2738 S3NWZC = DDOT(KZCA,VECA(KAZC),1,S3TRS(KAZC),1) 2739 S3NWZO = DDOT(KZOA,VECA(KAZO),1,S3TRS(KAZO),1) 2740 S3NWYC = DDOT(KZCA,VECA(KAYC),1,S3TRS(KAYC),1) 2741 S3NWYO = DDOT(KZOA,VECA(KAYO),1,S3TRS(KAYO),1) 2742 S3NEW = S3NWZC + S3NWZO + S3NWYC + S3NWYO 2743 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 5 :',S3NEW-S3TERM 2744 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 5 ZC :',S3NWZC-S3TMZC 2745 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 5 ZO :',S3NWZO-S3TMZO 2746 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 5 YC :',S3NWYC-S3TMYC 2747 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Case 5 YO :',S3NWYO-S3TMYO 2748 S3TMZC = S3NWZC 2749 S3TMZO = S3NWZO 2750 S3TMYC = S3NWYC 2751 S3TMYO = S3NWYO 2752 S3TERM = S3NEW 2753 END IF 2754C 2755 IF (IPRRSP .GT. 100 ) THEN 2756 WRITE(LUPRI,'(A)') ' Case 5 for S3 term' 2757 CALL OUTPUT(S3TRS,1,KZYVA,1,NSIM,KZYVA,NSIM,1,LUPRI) 2758 END IF 2759C 2760C 27616000 CONTINUE 2762C 2763C End of transformation; print final result if required 2764C 2765 IF (ATEST) THEN 2766 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Final result :',S3TERM 2767 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Final result ZC :',S3TMZC 2768 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Final result ZO :',S3TMZO 2769 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Final result YC :',S3TMYC 2770 WRITE(LUPRI,'(A,F24.8)') ' S3TEST Final result YO :',S3TMYO 2771 END IF 2772 CALL QEXIT('S3DRV') 2773 RETURN 2774 END 2775C /* Deck gtzymt */ 2776 SUBROUTINE GTZYMT(NSIM,VEC,KZYV,ISYMV,ZYMAT,MJWOP) 2777C 2778C This subroutine unpacks the ZY matrix from the vector. It uses 2779C the Z and the Y part of the vector. 2780C 2781#include "implicit.h" 2782#include "maxorb.h" 2783#include "maxash.h" 2784#include "priunit.h" 2785#include "infvar.h" 2786#include "inforb.h" 2787#include "infind.h" 2788#include "infrsp.h" 2789#include "wrkrsp.h" 2790#include "qrinf.h" 2791#include "infpri.h" 2792C 2793 DIMENSION VEC(KZYV) 2794 DIMENSION ZYMAT(NORBT,NORBT), MJWOP(2,MAXWOP,8) 2795C 2796 IF (NSIM .NE. 1) THEN 2797 CALL QUIT('GTZYMT ERROR: NSIM .ne. 1 not implemented') 2798 END IF 2799 IF( IPRRSP .GT. 200 ) THEN 2800 WRITE(LUPRI,'(/A)') ' GTZYMT called with the vector' 2801 WRITE(LUPRI,'(A)') ' =============================' 2802 CALL OUTPUT(VEC,1,KZYV,1,NSIM,KZYV,NSIM,1,LUPRI) 2803 END IF 2804C 2805 IOB = MZCONF(ISYMV) 2806 IOBT = MZVAR(ISYMV) + MZCONF(ISYMV) 2807 CALL DZERO(ZYMAT,NORBT*NORBT*NSIM) 2808 DO 100 ISIM = 1, NSIM 2809 DO 200 IG = 1, MZWOPT(ISYMV) 2810 I = MJWOP(1,IG,ISYMV) 2811 J = MJWOP(2,IG,ISYMV) 2812 ZYMAT(J,I) = VEC(IOB + IG) 2813 ZYMAT(I,J) = VEC(IOBT + IG) 2814200 CONTINUE 2815C 2816 IF( IPRRSP .GT. 200 ) THEN 2817 WRITE(LUPRI,'(/A)') ' Vector unpacked in matrix' 2818 WRITE(LUPRI,'(A)') ' =========================' 2819 CALL OUTPUT(ZYMAT,1,NORBT,1,NORBT,NORBT,NORBT, 2820 * 1,LUPRI) 2821 END IF 2822C 2823100 CONTINUE 2824C 2825 RETURN 2826 END 2827C /* Deck trzymt */ 2828 SUBROUTINE TRZYMT(NSIM,VEC1,VEC2,KZYV1,KZYV2,ISYM1,ISYM2,ZYMAT, 2829 * MJWOP,WRK,LWRK) 2830C 2831C This subroutine unpacks the ZY matrices from the two vectors and 2832C does the transformation 2833C 2834#include "implicit.h" 2835#include "maxorb.h" 2836#include "maxash.h" 2837#include "priunit.h" 2838#include "infvar.h" 2839#include "inforb.h" 2840#include "infind.h" 2841#include "infrsp.h" 2842#include "wrkrsp.h" 2843#include "qrinf.h" 2844#include "infpri.h" 2845C 2846 DIMENSION VEC1(KZYV1), VEC2(KZYV2) 2847 DIMENSION ZYMAT(NORBT,NORBT), MJWOP(2,MAXWOP,8) 2848 DIMENSION WRK(*) 2849C 2850 IF (NSIM .NE. 1) THEN 2851 CALL QUIT('TRZYMT ERROR: NSIM .ne. 1 not implemented') 2852 END IF 2853 KZY1 = 1 2854 KZY2 = KZY1 + NSIM * NORBT * NORBT 2855 KFREE = KZY2 + NSIM * NORBT * NORBT 2856 LFREE = LWRK - KFREE 2857 IF (LFREE.LT.0) CALL ERRWRK('TRZYMT',KFREE-1,LWRK) 2858C 2859 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYM1,WRK(KZY1),MJWOP ) 2860 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYM2,WRK(KZY2),MJWOP ) 2861C 2862 DO 100 ISIM = 1, NSIM 2863 IZY1 = (ISIM - 1) * NORBT * NORBT + KZY1 2864 IZY2 = (ISIM - 1) * NORBT * NORBT + KZY2 2865C 2866 CALL DGEMM('N','N',NORBT,NORBT,NORBT,1.D0, 2867 & WRK(IZY1),NORBT, 2868 & WRK(IZY2),NORBT,0.D0, 2869 & ZYMAT,NORBT) 2870C 2871 CALL DGEMM('N','N',NORBT,NORBT,NORBT,-1.D0, 2872 & WRK(IZY2),NORBT, 2873 & WRK(IZY1),NORBT,1.D0, 2874 & ZYMAT,NORBT) 2875C 2876 IF( IPRRSP .GT. 200 ) THEN 2877 WRITE(LUPRI,'(/A)') ' Final result in TRZYMT' 2878 WRITE(LUPRI,'(A)') ' ======================' 2879 CALL OUTPUT(ZYMAT,1,NORBT,1,NORBT,NORBT,NORBT, 2880 * 1,LUPRI) 2881 END IF 2882C 2883100 CONTINUE 2884C 2885 RETURN 2886 END 2887C /* Deck x2init */ 2888 SUBROUTINE X2INIT(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 2889 * GRVEC,VEC,X2TRS,XINDX,UDV,PV,OPLBL,IOPSYM,IOPSPI, 2890 * CMO,MJWOP,WRK,LWRK) 2891C 2892C Layout the core for the calculation of X2 times one vector 2893C 2894#include "implicit.h" 2895#include "infdim.h" 2896#include "inforb.h" 2897#include "maxash.h" 2898#include "priunit.h" 2899#include "infvar.h" 2900#include "infrsp.h" 2901#include "wrkrsp.h" 2902#include "rspprp.h" 2903#include "infhyp.h" 2904#include "qrinf.h" 2905#include "infpri.h" 2906#include "infspi.h" 2907C 2908 CHARACTER*8 OPLBL 2909C 2910 DIMENSION WRK(*) 2911 DIMENSION X2TRS(KZYVR), MJWOP(2,MAXWOP,8) 2912 DIMENSION GRVEC(KZYVR), VEC(KZYV) 2913 DIMENSION XINDX(*) 2914 DIMENSION CMO(*) 2915 DIMENSION UDV(NASHDI,NASHDI), PV(*) 2916C 2917C Any variables in this symmetry? 2918C 2919 IF (KZYVR .EQ. 0) RETURN 2920C 2921 IF( (NSIM .GT. 1) ) THEN 2922 WRITE(LUERR,'(//A,/A,I5,/A)') 2923 * ' --> Error : Routine X2INIT not implemented with NSIM > 1', 2924 * ' NSIM = ', NSIM, 2925 * ' Halting execution.' 2926 CALL QTRACE(LUERR) 2927 CALL QUIT(' RSPQR ERROR: NSIM GREATER THAN 1 IN X2INIT') 2928 END IF 2929C 2930C Initialise the gradient to zero. 2931C 2932 CALL DZERO(X2TRS,KZYVR) 2933C 2934 KOPMAT = 1 2935 KDEN1 = KOPMAT + NORBT * NORBT 2936 KDEN2 = KDEN1 + NASHT * NASHT 2937 IF (OPLBL(3:8).EQ.'SPNORB' .OR. OPLBL(3:8).EQ.'MNF-SO' 2938 & .OR. OPLBL(3:8).EQ.'SPNSCA') THEN 2939 KFREE = KDEN2 + N2ASHX*N2ASHX*2 2940 ELSE 2941 KFREE = KDEN2 2942 END IF 2943 LWRKF = LWRK - KFREE + 1 2944 IF (LWRKF.LT.0) CALL ERRWRK('X2INIT',KFREE-1,LWRK) 2945C 2946 IF (IPRRSP .GT. 50) THEN 2947 WRITE(LUPRI,'(//A)') ' Characteristics of X2 gradient' 2948 WRITE(LUPRI,'(A)') ' ==============================' 2949 WRITE(LUPRI,'(A,I8)') ' Gradient symmetry :',IGRSYM 2950 WRITE(LUPRI,'(A,I8)') ' Gradient spin :',IGRSPI 2951 WRITE(LUPRI,'(A,I8)') ' Gradient length :',KZYVR 2952 WRITE(LUPRI,'(A,I8)') ' Vector symmetry :',ISYMV 2953 WRITE(LUPRI,'(A,I8)') ' Vector spin :',ISPINV 2954 WRITE(LUPRI,'(A,I8)') ' Vector length :',KZYV 2955 WRITE(LUPRI,'(A,I8)') ' Operator symmetry :',IOPSYM 2956 WRITE(LUPRI,'(A,I8)') ' Operator spin :',IOPSPI 2957 WRITE(LUPRI,'(A,A8)') ' Operator label :',OPLBL 2958 WRITE(LUPRI,'(A//)') ' ==============================' 2959 END IF 2960C 2961 IF (OPLBL(3:8).EQ.'SPNORB' .OR. OPLBL(3:8).EQ.'MNF-SO' 2962 & .OR. OPLBL(3:8).EQ.'SPNSCA') THEN 2963 CALL X2HSO(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 2964 * OPLBL,IOPSYM,IOPSPI,GRVEC,X2TRS, 2965 * VEC,UDV,PV,WRK(KOPMAT),WRK(KDEN1),WRK(KDEN2), 2966 * XINDX,CMO,MJWOP,WRK(KFREE),LWRKF) 2967 ELSE 2968 CALL X2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 2969 * OPLBL,IOPSYM,IOPSPI,GRVEC,X2TRS, 2970 * VEC,UDV,WRK(KOPMAT),WRK(KDEN1), 2971 * XINDX,CMO,MJWOP,WRK(KFREE),LWRKF) 2972 END IF 2973C 2974 RETURN 2975 END 2976C /* Deck x2drv */ 2977 SUBROUTINE X2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 2978 * OPLBL,IOPSYM,IOPSPI, 2979 * GRVEC,X2TRS,VEC,UDV,OPMAT,DEN1,XINDX,CMO, 2980 * MJWOP,WRK,LWRK) 2981C 2982C Drive the computation of X[2] times one vector 2983C 2984#include "implicit.h" 2985C 2986 PARAMETER (D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0 ) 2987C 2988#include "maxorb.h" 2989#include "maxash.h" 2990#include "priunit.h" 2991#include "inforb.h" 2992#include "infind.h" 2993#include "infvar.h" 2994#include "infdim.h" 2995#include "infrsp.h" 2996#include "wrkrsp.h" 2997#include "rspprp.h" 2998#include "infhyp.h" 2999#include "qrinf.h" 3000#include "infpri.h" 3001#include "infspi.h" 3002C 3003 CHARACTER*8 OPLBL 3004C 3005 DIMENSION WRK(*) 3006 DIMENSION X2TRS(KZYVR), MJWOP(2,MAXWOP,8) 3007 DIMENSION GRVEC(KZYVR), VEC(KZYV) 3008 DIMENSION DEN1(NASHDI,NASHDI) 3009 DIMENSION OPMAT(NORBT,NORBT) 3010 DIMENSION XINDX(*) 3011 DIMENSION CMO(*) 3012 DIMENSION UDV(NASHDI,NASHDI) 3013C 3014 LOGICAL NORHO2, TDM, LORB, LCON, LREFST, LREF 3015C 3016C Layout some workspace 3017C 3018 KCREF = 1 3019 KZYM = KCREF + MZCONF(1) 3020 KX2X = KZYM + N2ORBX 3021 KFREE = KX2X + N2ORBX 3022 LFREE = LWRK - KFREE 3023 IF (LFREE.LT.0) CALL ERRWRK('X2DRV 1',KFREE-1,LWRK) 3024C 3025 IF (X2TEST) THEN 3026 X2TMZC = D0 3027 X2TMZO = D0 3028 X2TMYC = D0 3029 X2TMYO = D0 3030 X2TM = D0 3031 KZX2O = MZWOPT(IGRSYM) 3032 KZX2C = MZCONF(IGRSYM) 3033 KZC = 1 3034 KZO = KZC + KZX2C 3035 KYC = KZO + KZX2O 3036 KYO = KYC + KZX2C 3037 END IF 3038C 3039 KSYMOP = IOPSYM 3040 TDM = .TRUE. 3041 NORHO2 = .TRUE. 3042C 3043C Get the operator matrix 3044C Put the minus sign in the whole term into the operator 3045C 3046 KSYMP = -1 3047 CALL PRPGET(OPLBL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRRSP) 3048 IF (KSYMP .NE. KSYMOP) CALL QUIT( 3049 & 'ERROR: Unexpected symmetry of property matrix') 3050 CALL DSCAL(NORBT*NORBT,DM1,OPMAT,1) 3051C 3052C Case 1 3053C ====== 3054 IF (MZWOPT(ISYMV) .EQ. 0 ) GOTO 2000 3055C 3056 ISPIN = MULSP(ISPINV,IOPSPI) 3057 IF (NASHT.EQ.0 .AND. ISPIN.NE.IGRSPI) THEN 3058 WRITE(LUPRI,*) "FOUND: ", ISPIN, " Expected:", IGRSPI 3059 WRITE(LUPRI,*) "parameters: ", ispinv, iopspi 3060 CALL QUIT('X2DRV: SPIN ERROR') 3061 ENDIF 3062C 3063C Transform the operator 3064C 3065 CALL GTZYMT(NSIM,VEC,KZYV,ISYMV,WRK(KZYM),MJWOP) 3066 CALL DZERO(WRK(KX2X),N2ORBX) 3067 CALL OITH1(ISYMV,WRK(KZYM),OPMAT,WRK(KX2X),IOPSYM) 3068C 3069C Make copies of the MC density matrix in DEN1 3070C 3071 IF (ISPIN.EQ.IGRSPI) THEN 3072 CALL DCOPY(N2ASHX,UDV,1,DEN1,1) 3073 ELSE 3074 CALL DCOPY(N2ASHX,UDV(1+N2ASHX,1),1,DEN1,1) 3075 END IF 3076 OVLAP = D1 3077 ISYMDN = 1 3078C 3079 CALL GETREF(WRK(KCREF),MZCONF(1)) 3080C 3081C Make the gradient 3082C 3083 ISYMR = IREFSY 3084 INTSYM = MULD2H(KSYMOP,ISYMV) 3085 ISYMST = MULD2H(IGRSYM,IREFSY) 3086 IF ( ISYMST .EQ. IREFSY ) THEN 3087 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 3088 ELSE 3089 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 3090 END IF 3091 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 3092 NZYVEC = MZCONF(1) 3093 NZCVEC = MZCONF(1) 3094 LREFST = .TRUE. 3095 CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,X2TRS, 3096 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,WRK(KX2X), 3097 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREFST) 3098C 3099 IF (X2TEST) THEN 3100 X2ZO = DDOT(KZX2O,GRVEC(KZO),1,X2TRS(KZO),1) 3101 X2ZC = DDOT(KZX2C,GRVEC(KZC),1,X2TRS(KZC),1) 3102 X2YO = DDOT(KZX2O,GRVEC(KYO),1,X2TRS(KYO),1) 3103 X2YC = DDOT(KZX2C,GRVEC(KYC),1,X2TRS(KYC),1) 3104 X2TOT = X2ZO+X2ZC+X2YO+X2YC 3105 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 ZO :',X2ZO - X2TMZO 3106 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 ZC :',X2ZC - X2TMZC 3107 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 YO :',X2YO - X2TMYO 3108 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 YC :',X2YC - X2TMYC 3109 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 1 :',X2TOT - X2TM 3110 X2TMZO = X2ZO 3111 X2TMZC = X2ZC 3112 X2TMYO = X2YO 3113 X2TMYC = X2YC 3114 X2TM = X2TOT 3115 END IF 3116 IF ( IPRRSP .GT. 100 ) THEN 3117 WRITE(LUPRI,'(A)') ' Case 1 for X2 term' 3118 CALL OUTPUT(X2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3119 END IF 3120C 3121C Case 2 3122C ====== 3123 2000 IF (MZCONF(ISYMV) .EQ. 0 ) GOTO 3000 3124C Construct the density matrix <02L|..|0> + <0|..|02R> 3125C 3126C 3127 CALL GETREF(WRK(KCREF),MZCONF(1)) 3128C 3129 CALL DZERO(DEN1,NSIM*NASHT*NASHT) 3130 OVLAP = D0 3131C 3132 ISPIN1 = IGRSPI 3133 ISPIN2 = IOPSPI 3134 ILSYM = IREFSY 3135 IRSYM = MULD2H(IREFSY,ISYMV) 3136 NCL = MZCONF(1) 3137 NCR = MZCONF(ISYMV) 3138 KZVARL = MZCONF(1) 3139 KZVARR = MZYVAR(ISYMV) 3140 LREF = .TRUE. 3141 ISYMDN = MULD2H(ILSYM,IRSYM) 3142 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 3143 * WRK(KCREF),VEC,OVLAP,DEN1,DUMMY,ISPIN1,ISPIN2, 3144 * TDM,NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 3145C 3146C Make the gradient 3147C 3148 ISYMR = ISYMV 3149 INTSYM = KSYMOP 3150 ISPIN = IOPSPI 3151 ISYMST = MULD2H(IGRSYM,IREFSY) 3152 IF ( ISYMST .EQ. IREFSY ) THEN 3153 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 3154 ELSE 3155 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 3156 END IF 3157 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 3158 NZYVEC = MZYVAR(ISYMV) 3159 NZCVEC = MZCONF(ISYMV) 3160 LREFST = .FALSE. 3161 CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,X2TRS, 3162 * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,OPMAT, 3163 * XINDX,MJWOP,WRK,LWRK,LORB,LCON,LREFST) 3164C 3165 IF (X2TEST) THEN 3166 X2ZO = DDOT(KZX2O,GRVEC(KZO),1,X2TRS(KZO),1) 3167 X2ZC = DDOT(KZX2C,GRVEC(KZC),1,X2TRS(KZC),1) 3168 X2YO = DDOT(KZX2O,GRVEC(KYO),1,X2TRS(KYO),1) 3169 X2YC = DDOT(KZX2C,GRVEC(KYC),1,X2TRS(KYC),1) 3170 X2TOT = X2ZO+X2ZC+X2YO+X2YC 3171 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 ZO :',X2ZO - X2TMZO 3172 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 ZC :',X2ZC - X2TMZC 3173 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 YO :',X2YO - X2TMYO 3174 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 YC :',X2YC - X2TMYC 3175 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 2 :',X2TOT - X2TM 3176 X2TMZO = X2ZO 3177 X2TMZC = X2ZC 3178 X2TMYO = X2YO 3179 X2TMYC = X2YC 3180 X2TM = X2TOT 3181 END IF 3182 IF ( IPRRSP .GT. 100 ) THEN 3183 WRITE(LUPRI,'(A)') ' Case 2 for X2 term' 3184 CALL OUTPUT(X2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3185 END IF 3186C 3187C Case 3 3188C ====== 3189C Do Sj(2) <0 |op| 0> 3190C 3191C Compute the expectation value 3192C 3193 IF (KSYMOP .EQ. 1) THEN 3194C 3195 ISPIN1 = IOPSPI 3196 ISPIN2 = 0 3197 TRPLET = IOPSPI.EQ.1 3198 IF (RSPCI.OR.TRPLET) THEN 3199 KCREF = 1 3200 KDEN1 = KCREF + MZCONF(1) 3201 KFREE = KDEN1 + N2ASHX 3202 LFREE = LWRK - KFREE 3203 IF (KFREE.LT.0) CALL ERRWRK('X2DRV 2',KFREE-1,LWRK) 3204 CALL GETREF(WRK(KCREF),MZCONF(1)) 3205 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 3206 * WRK(KCREF),WRK(KCREF),WRK(KDEN1),RHO2, 3207 * ISPIN1,ISPIN2,.FALSE.,.TRUE.,XINDX,WRK, 3208 * KFREE,LFREE) 3209 IPRONE = 75 3210 CALL PRPONE(OPMAT,WRK(KDEN1),FACT,IPRONE,' Sj(2) <0|op|0>') 3211 ELSE 3212 IPRONE = 75 3213 CALL PRPONE(OPMAT,UDV,FACT,IPRONE,' Sj(2) <0|op|0>') 3214 END IF 3215C 3216 IF (FACT .NE. D0) THEN 3217 NZCONF = MZCONF(IGRSYM) 3218 NZVAR = MZVAR(IGRSYM) 3219 CALL DAXPY(NZCONF,FACT,VEC ,1,X2TRS ,1) 3220 CALL DAXPY(NZCONF,FACT,VEC(NZVAR+1),1,X2TRS(NZVAR+1),1) 3221 END IF 3222 END IF 3223C 3224 IF (X2TEST) THEN 3225 X2ZO = DDOT(KZX2O,GRVEC(KZO),1,X2TRS(KZO),1) 3226 X2ZC = DDOT(KZX2C,GRVEC(KZC),1,X2TRS(KZC),1) 3227 X2YO = DDOT(KZX2O,GRVEC(KYO),1,X2TRS(KYO),1) 3228 X2YC = DDOT(KZX2C,GRVEC(KYC),1,X2TRS(KYC),1) 3229 X2TOT = X2ZO+X2ZC+X2YO+X2YC 3230 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 3 ZO :',X2ZO - X2TMZO 3231 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 3 ZC :',X2ZC - X2TMZC 3232 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 3 YO :',X2YO - X2TMYO 3233 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 3 YC :',X2YC - X2TMYC 3234 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Case 3 :',X2TOT - X2TM 3235 X2TMZO = X2ZO 3236 X2TMZC = X2ZC 3237 X2TMYO = X2YO 3238 X2TMYC = X2YC 3239 X2TM = X2TOT 3240 END IF 3241 IF ( IPRRSP .GT. 100 ) THEN 3242 WRITE(LUPRI,'(A)') ' Case 3 for X2 term' 3243 CALL OUTPUT(X2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3244 END IF 3245C 3246 3000 CONTINUE 3247 IF (X2TEST) THEN 3248 WRITE(LUPRI,'(/A,F24.8)')' X2TEST Final result: ZO', X2ZO 3249 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: ZC', X2ZC 3250 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: YO', X2YO 3251 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: YC', X2YC 3252 WRITE(LUPRI,'(A,F24.8)') ' X2TEST Final result: ', X2TOT 3253 END IF 3254 RETURN 3255 END 3256C /* Deck a2init */ 3257 SUBROUTINE A2INIT(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 3258 * GRVEC,VEC,A2TRS,XINDX,UDV,PV, 3259 * OPLBL,IOPSYM,IOPSPI, 3260 * CMO,MJWOP,WRK,LWRK) 3261C 3262C Layout the core for the calculation of A[2] times one vector 3263C 3264#include "implicit.h" 3265#include "maxash.h" 3266#include "priunit.h" 3267#include "infdim.h" 3268#include "inforb.h" 3269#include "infvar.h" 3270#include "infrsp.h" 3271#include "wrkrsp.h" 3272#include "qrinf.h" 3273#include "infpri.h" 3274C 3275 CHARACTER*8 OPLBL 3276C 3277 DIMENSION WRK(*) 3278 DIMENSION A2TRS(KZYVR), MJWOP(2,MAXWOP,8) 3279 DIMENSION GRVEC(KZYVR), VEC(KZYV) 3280 DIMENSION XINDX(*) 3281 DIMENSION UDV(NASHDI,NASHDI) 3282 DIMENSION CMO(*) 3283C 3284C Any variables in this symmetry? 3285C 3286 IF (KZYVR .EQ. 0) RETURN 3287C 3288 IF( (NSIM .GT. 1) ) THEN 3289 WRITE(LUERR,'(//A,/A,I5,/A)') 3290 * ' --> Error : Routine A2INIT not implemented with NSIM > 1', 3291 * ' NSIM = ', NSIM, 3292 * ' Halting execution.' 3293 CALL QTRACE(LUERR) 3294 CALL QUIT(' RSPQR ERROR: NSIM GREATER THAN 1 IN A2INIT') 3295 END IF 3296C 3297C Initialise the gradient to zero. 3298C 3299 CALL DZERO(A2TRS,KZYVR) 3300C 3301 KOPMAT = 1 3302 KFREE = KOPMAT + NORBT * NORBT 3303 LWRKF = LWRK - KFREE + 1 3304 IF (LWRKF.LT.0) CALL ERRWRK('A2INIT',KFREE-1,LWRK) 3305 CALL DZERO(WRK,KFREE-1) 3306C 3307 IF (IPRRSP .GT. 50 ) THEN 3308 WRITE(LUPRI,'(//A)') ' Characteristics of A2 gradient' 3309 WRITE(LUPRI,'(A)') ' ==============================' 3310 WRITE(LUPRI,'(A,I8)') ' Gradient symmetry :',IGRSYM 3311 WRITE(LUPRI,'(A,I8)') ' Gradient spin :',IGRSPI 3312 WRITE(LUPRI,'(A,I8)') ' Gradient length :',KZYVR 3313 WRITE(LUPRI,'(A,I8)') ' Vector symmetry :',ISYMV 3314 WRITE(LUPRI,'(A,I8)') ' Vector spin :',ISPINV 3315 WRITE(LUPRI,'(A,I8)') ' Vector length :',KZYV 3316 WRITE(LUPRI,'(A,I8)') ' Operator symmetry :',IOPSYM 3317 WRITE(LUPRI,'(A,I8)') ' Operator spin :',IOPSPI 3318 WRITE(LUPRI,'(A,A8)') ' Operator label :',OPLBL 3319 WRITE(LUPRI,'(A//)') ' ==============================' 3320 END IF 3321C 3322 IF (OPLBL(3:8).EQ.'SPNORB' .OR. OPLBL(3:8).EQ.'MNF-SO' 3323 & .OR. OPLBL(3:8).EQ.'SPNSCA') THEN 3324 CALL A2HSO(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 3325 * OPLBL,IOPSYM,IOPSPI, 3326 * GRVEC,A2TRS,VEC,UDV,PV,WRK(KOPMAT),XINDX,CMO, 3327 * MJWOP,WRK(KFREE),LWRKF) 3328 ELSE 3329 CALL A2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 3330 * OPLBL,IOPSYM,IOPSPI, 3331 * GRVEC,A2TRS,VEC,UDV,WRK(KOPMAT),XINDX,CMO, 3332 * MJWOP,WRK(KFREE),LWRKF) 3333C 3334 END IF 3335 RETURN 3336 END 3337C /* Deck a2drv */ 3338 SUBROUTINE A2DRV(NSIM,KZYVR,KZYV,IGRSYM,IGRSPI,ISYMV,ISPINV, 3339 * OPLBL,IOPSYM,IOPSPI, 3340 * GRVEC,A2TRS,VEC,UDV,OPMAT,XINDX,CMO, 3341 * MJWOP,WRK,LWRK) 3342C 3343C Drive the computation of A[2] times one vector 3344C 3345#include "implicit.h" 3346C 3347 PARAMETER( D0 = 0.0D0, D1 = 1.0D0, DH = 0.5D0 ) 3348C 3349#include "maxorb.h" 3350#include "maxash.h" 3351#include "priunit.h" 3352#include "infdim.h" 3353#include "inforb.h" 3354#include "infrsp.h" 3355#include "wrkrsp.h" 3356#include "infvar.h" 3357#include "infind.h" 3358#include "qrinf.h" 3359#include "infpri.h" 3360#include "infspi.h" 3361C 3362 CHARACTER*8 OPLBL 3363C 3364 DIMENSION WRK(*) 3365 DIMENSION A2TRS(KZYVR), MJWOP(2,MAXWOP,8) 3366 DIMENSION GRVEC(KZYVR), VEC(KZYV) 3367 DIMENSION OPMAT(NORBT,NORBT) 3368 DIMENSION XINDX(*) 3369 DIMENSION UDV(NASHDI,NASHDI) 3370 DIMENSION CMO(*) 3371C 3372 LOGICAL NORHO2, TDM, LORB, LCON, LREFST, LREF 3373C 3374C Layout some workspace 3375C 3376 KCREF = 1 3377 KZYM = KCREF + MZCONF(1) 3378 KA2X = KZYM + N2ORBX 3379 KFREE = KA2X + N2ORBX 3380 LFREE = LWRK - KFREE 3381 IF (LFREE.LT.0) CALL ERRWRK('A2DRV 1',KFREE-1,LWRK) 3382C 3383 TDM = .TRUE. 3384 NORHO2 = .TRUE. 3385C 3386 IF (A2TEST) THEN 3387 A2TMZC = D0 3388 A2TMZO = D0 3389 A2TMYC = D0 3390 A2TMYO = D0 3391 A2TM = D0 3392 KZA2O = MZWOPT(IGRSYM) 3393 KZA2C = MZCONF(IGRSYM) 3394 KZC = 1 3395 KZO = KZC + KZA2C 3396 KYC = KZO + KZA2O 3397 KYO = KYC + KZA2C 3398 END IF 3399C Get the operator matrix 3400C 910408-hjaaj: transpose the operator matrix to 3401C get right sign for imaginary operators 3402C 3403 KSYMOP = IOPSYM 3404 KSYMP = -1 3405 CALL PRPGET(OPLBL,CMO,OPMAT,KSYMP,ANTSYM,WRK,LWRK,IPRRSP) 3406 IF (KSYMP .NE. KSYMOP) CALL QUIT( 3407 & 'ERROR: Unexpected symmetry of property matrix') 3408 CALL DGETRN(OPMAT,NORBT,NORBT) 3409C 3410C Case 1: 3411C ======= 3412C / 1/2 <0| [qj+,A(K)] |0> \ 3413C | - <0| A(K) |j> | 3414C | 1/2 <0| [qj ,A(K)] |0> | 3415C \ <j| A(K) |0> / 3416C 3417 IF ( MZWOPT(ISYMV) .EQ. 0 ) GOTO 2000 3418C 3419 ISPIN = MULSP(ISPINV,IOPSPI) 3420 IF (NASHT.EQ.0 .AND. ISPIN.NE.IGRSPI) THEN 3421 WRITE(LUPRI,*) "FOUND: ", ISPIN, " Expected:", IGRSPI 3422 WRITE(LUPRI,*) "parameters: ", ispinv, iopspi 3423 CALL QUIT('A2DRV: SPIN ERROR') 3424 END IF 3425C 3426C Transform the operator 3427C 3428 INTSYM = MULD2H(ISYMV,IOPSYM) 3429 CALL GTZYMT(NSIM,VEC,KZYV,ISYMV,WRK(KZYM),MJWOP) 3430 CALL DZERO(WRK(KA2X),N2ORBX) 3431 CALL OITH1(ISYMV,WRK(KZYM),OPMAT,WRK(KA2X),IOPSYM) 3432C 3433 CALL GETREF(WRK(KCREF),MZCONF(1)) 3434C 3435C Make the gradient 3436C 3437 OVLAP = D1 3438 ISYMDN = 1 3439 ISYMR = IREFSY 3440 ISYMST = MULD2H(IGRSYM,IREFSY) 3441 IF(ISYMST .EQ. IREFSY ) THEN 3442 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 3443 ELSE 3444 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 3445 END IF 3446 LORB = (MZWOPT(IGRSYM) .GT. 0 ) 3447 NZYVEC = MZCONF(1) 3448 NZCVEC = MZCONF(1) 3449 LREFST = .TRUE. 3450 IF (LCON) THEN 3451 CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 3452 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV, 3453 * WRK(KA2X),XINDX,MJWOP,WRK(KFREE),LFREE, 3454 * .FALSE.,LCON,LREFST) 3455 END IF 3456 IF (LORB) THEN 3457C multiply A(K) with 1/2 for orbital part 3458 CALL DSCAL(N2ORBX,DH,WRK(KA2X),1) 3459 CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 3460 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV, 3461 * WRK(KA2X),XINDX,MJWOP,WRK(KFREE),LFREE, 3462 * LORB,.FALSE.,LREFST) 3463 END IF 3464C 3465 IF (A2TEST) THEN 3466 A2ZO = DDOT(KZA2O,GRVEC(KZO),1,A2TRS(KYO),1) 3467 A2ZC = DDOT(KZA2C,GRVEC(KZC),1,A2TRS(KYC),1) 3468 A2YO = DDOT(KZA2O,GRVEC(KYO),1,A2TRS(KZO),1) 3469 A2YC = DDOT(KZA2C,GRVEC(KYC),1,A2TRS(KZC),1) 3470 A2TOT = A2ZO+A2ZC+A2YO+A2YC 3471 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 ZO :',A2ZO - A2TMZO 3472 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 ZC :',A2ZC - A2TMZC 3473 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 YO :',A2YO - A2TMYO 3474 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 YC :',A2YC - A2TMYC 3475 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 1 :',A2TOT - A2TM 3476 A2TMZO = A2ZO 3477 A2TMZC = A2ZC 3478 A2TMYO = A2YO 3479 A2TMYC = A2YC 3480 A2TM = A2TOT 3481 END IF 3482 IF ( IPRRSP .GT. 30 ) THEN 3483 WRITE(LUPRI,'(A)') ' Case 1 for A2 term' 3484 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3485 END IF 3486C 3487C Case 2: 3488C ======= 34892000 IF (MZCONF(ISYMV) .EQ. 0 ) GO TO 4000 3490C 3491C Multiply the factor of one half into the operator matrix 3492C for evaluation of Case 2 and 3. 3493C 3494 CALL DSCAL(NSIM*NORBT*NORBT,DH,OPMAT,1) 3495C 3496 ISPIN = IOPSPI 3497C 3498C Make the gradient 3499C 3500 OVLAP = D0 3501 ISYMDN = 1 3502 INTSYM = IOPSYM 3503 ISYMR = ISYMV 3504 NZYVEC = MZYVAR(ISYMV) 3505 NZCVEC = MZCONF(ISYMV) 3506 LORB = .FALSE. 3507 ISYMST = MULD2H(IGRSYM,IREFSY) 3508 IF(ISYMST .EQ. IREFSY ) THEN 3509 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 3510 ELSE 3511 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 3512 END IF 3513 LREFST = .FALSE. 3514 CALL RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI,ISYMR,A2TRS, 3515 * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,OPMAT, 3516 * XINDX,MJWOP,WRK,LWRK,LORB,LCON,LREFST) 3517C 3518 IF (A2TEST) THEN 3519 A2ZO = DDOT(KZA2O,GRVEC(KZO),1,A2TRS(KYO),1) 3520 A2ZC = DDOT(KZA2C,GRVEC(KZC),1,A2TRS(KYC),1) 3521 A2YO = DDOT(KZA2O,GRVEC(KYO),1,A2TRS(KZO),1) 3522 A2YC = DDOT(KZA2C,GRVEC(KYC),1,A2TRS(KZC),1) 3523 A2TOT = A2ZO+A2ZC+A2YO+A2YC 3524 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 ZO :',A2ZO - A2TMZO 3525 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 ZC :',A2ZC - A2TMZC 3526 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 YO :',A2YO - A2TMYO 3527 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 YC :',A2YC - A2TMYC 3528 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 2 :',A2TOT - A2TM 3529 A2TMZO = A2ZO 3530 A2TMZC = A2ZC 3531 A2TMYO = A2YO 3532 A2TMYC = A2YC 3533 A2TM = A2TOT 3534 END IF 3535 IF ( IPRRSP .GT. 100 ) THEN 3536 WRITE(LUPRI,'(A)') ' Case 2 for A2 term' 3537 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3538 END IF 3539C 3540C Case 3: 3541C ======= 3542C Do Sj(2) <0 |op| 0> 3543C 3544C Compute the expectation value (factor of 0.5 multiplied into OPMAT) 3545C 3546 IF ( IOPSYM .EQ. 1 ) THEN 3547C 3548 TRPLET = IOPSPI.EQ.1 3549 ISPIN1 = IOPSPI 3550 ISPIN2 = 0 3551 IF (RSPCI.OR.TRPLET) THEN 3552 KCREF = 1 3553 KDEN1 = KCREF + MZCONF(1) 3554 KFREE = KDEN1 + N2ASHX 3555 LFREE = LWRK - KFREE 3556 IF (LFREE.LT.0) CALL ERRWRK('A2DRV 2',KFREE-1,LWRK) 3557 CALL GETREF(WRK(KCREF),MZCONF(1)) 3558 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 3559 * WRK(KCREF),WRK(KCREF),WRK(KDEN1),RHO2, 3560 * ISPIN1,ISPIN2,.FALSE.,.TRUE.,XINDX,WRK, 3561 * KFREE,LFREE) 3562 IPRONE = 75 3563 CALL PRPONE(OPMAT,WRK(KDEN1),FACT,IPRONE, 3564 * ' A(2) TERM <0|op|0>') 3565 ELSE 3566 IPRONE = 75 3567 CALL PRPONE(OPMAT,UDV,FACT,IPRONE, 3568 * ' A(2) TERM <0|op|0>') 3569 END IF 3570C 3571 IF ( FACT .NE. D0 ) THEN 3572 NZCONF = MZCONF(IGRSYM) 3573 NZVAR = MZVAR(IGRSYM) 3574 CALL DAXPY(NZCONF,FACT,VEC,1,A2TRS,1) 3575 CALL DAXPY(NZCONF,FACT,VEC(NZVAR+1),1,A2TRS(NZVAR+1),1) 3576 END IF 3577 END IF 3578C 3579 IF (A2TEST) THEN 3580 A2ZO = DDOT(KZA2O,GRVEC(KZO),1,A2TRS(KYO),1) 3581 A2ZC = DDOT(KZA2C,GRVEC(KZC),1,A2TRS(KYC),1) 3582 A2YO = DDOT(KZA2O,GRVEC(KYO),1,A2TRS(KZO),1) 3583 A2YC = DDOT(KZA2C,GRVEC(KYC),1,A2TRS(KZC),1) 3584 A2TOT = A2ZO+A2ZC+A2YO+A2YC 3585 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 3 ZO :',A2ZO - A2TMZO 3586 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 3 ZC :',A2ZC - A2TMZC 3587 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 3 YO :',A2YO - A2TMYO 3588 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 3 YC :',A2YC - A2TMYC 3589 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Case 3 Tot:',A2TOT - A2TM 3590 A2TMZO = A2ZO 3591 A2TMZC = A2ZC 3592 A2TMYO = A2YO 3593 A2TMYC = A2YC 3594 A2TM = A2TOT 3595 END IF 3596 IF ( IPRRSP .GT. 100 ) THEN 3597 WRITE(LUPRI,'(A)') ' Case 3 for A2 term' 3598 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3599 END IF 3600C 36014000 CONTINUE 3602 IF (A2TEST) THEN 3603 WRITE(LUPRI,'(/A,F24.8)')' A2TEST Final result: ZO', A2ZO 3604 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: ZC', A2ZC 3605 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: YO', A2YO 3606 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: YC', A2YC 3607 WRITE(LUPRI,'(A,F24.8)') ' A2TEST Final result: ', A2TOT 3608 END IF 3609 IF ( IPRRSP .GT. 150 ) THEN 3610 WRITE(LUPRI,'(//A)') ' Gradient before swapping in A2DRV' 3611 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3612 END IF 3613C 3614C Swap the final result to conform with the notation for A[2] 3615C 3616 CALL DSWAP(MZVAR(IGRSYM), A2TRS,1, A2TRS(MZVAR(IGRSYM)+1),1) 3617C 3618 IF ( IPRRSP .GT. 100 ) THEN 3619 WRITE(LUPRI,'(//A)') ' Gradient after swapping in A2DRV' 3620 CALL OUTPUT(A2TRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3621 END IF 3622C 3623 RETURN 3624 END 3625C /* Deck rsp1gr */ 3626 SUBROUTINE RSP1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI, 3627 * ISYMV,OTRS, 3628 * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,OPMAT, 3629 * XINDX,MJWOP,WRK,LWRK,LORB,LCON,LREFST) 3630C 3631C Compute the gradient resulting from multiplying the S matrix 3632C with one ore more vectors. 3633C 3634#include "implicit.h" 3635#include "infdim.h" 3636#include "inforb.h" 3637#include "priunit.h" 3638#include "infrsp.h" 3639#include "wrkrsp.h" 3640#include "rspprp.h" 3641#include "infhyp.h" 3642#include "infvar.h" 3643#include "qrinf.h" 3644#include "infpri.h" 3645C 3646 LOGICAL LORB, LCON, LREFST 3647C 3648 DIMENSION OTRS (KZYVR), MJWOP(2,MAXWOP,8) 3649 DIMENSION VEC(NZYVEC) 3650 DIMENSION DEN1(NASHDI,NASHDI) 3651 DIMENSION OPMAT(NORBT,NORBT) 3652 DIMENSION XINDX(*) 3653 DIMENSION WRK(*) 3654C 3655 IF ( IPRRSP .GT. 150 ) THEN 3656 WRITE(LUPRI,'(A)') ' Vector in RSP1GR' 3657 IF (LREFST) WRITE(LUPRI,'(A)') ' (Reference state)' 3658 CALL OUTPUT(VEC,1,NZYVEC,1,NSIM,NZYVEC,NSIM,1,LUPRI) 3659 IF ( LORB ) THEN 3660 WRITE(LUPRI,'(//A)') ' Density matrix in RSP1GR' 3661 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 3662 END IF 3663 WRITE(LUPRI,'(//A)') ' One electron matrix in RSP1GR' 3664 CALL OUTPUT(OPMAT,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 3665 END IF 3666C 3667 IF ( LORB ) THEN 3668 TRPLET = IGRSPI.NE.ISPIN 3669 CALL ORBSX(NSIM,IGRSYM,KZYVR,OTRS,OPMAT,OVLAP, 3670 * ISYMDN,DEN1,MJWOP,WRK,LWRK) 3671 END IF 3672C 3673 IF ( LCON ) THEN 3674 ISYMJ = MULD2H( IGRSYM, IREFSY ) 3675 NZCSTJ = MZCONF( IGRSYM ) 3676C 3677 TRPLET = ISPIN.EQ.1 3678 CALL CONSX(NSIM,KZYVR,IGRSYM,OPMAT,VEC,NZYVEC, 3679 * NZCVEC,ISYMV,NZCSTJ,ISYMJ,LREFST,OTRS,XINDX, 3680 * WRK,LWRK) 3681 END IF 3682C 3683 RETURN 3684 END 3685C /* Deck orbsx */ 3686 SUBROUTINE ORBSX(NSIM,IGRSYM,KZYVR,OTRS,OPMAT,OVLAP,ISYMDN,DEN1, 3687 * MJWOP,WRK,LWRK) 3688C 3689C Pupose: 3690C Get the gradient vector from the zymat one -electron matrix. 3691C We have the expression: 3692C 3693C k(l,k) = <0| [ E(k,l), zymat ] |0> 3694C 3695C We have the expressions: 3696C 3697C 1) k(a,i) = 2 * OVLAP * zymat(a,i) 3698C 3699C 2) k(i,a) = - 2 * OVLAP * zymat(i,a) 3700C 3701C 3) k(t,i) = 2 * OVLAP * zymat(t,i)- sum(x) zymat(x,i) D(x,t) 3702C 3703C 4) k(i,t) = sum(x) zymat(i,x) D(t,x) - 2 * OVLAP * zymat(i,t) 3704C 3705C 5) k(t,a) = - sum(x) zymat(x,a) D(t,x) 3706C 3707C 6) k(a,t) = sum(x) zymat(a,x) D(t,x) 3708C 3709C 7) k(u,t) = sum(x) zymat(u,x) D(t,x) - zymat(x,t) D(x,u) 3710C 3711#include "implicit.h" 3712#include "maxorb.h" 3713#include "maxash.h" 3714#include "priunit.h" 3715#include "infvar.h" 3716#include "inforb.h" 3717#include "infind.h" 3718#include "infdim.h" 3719#include "infpri.h" 3720#include "infrsp.h" 3721#include "wrkrsp.h" 3722#include "qrinf.h" 3723C 3724 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0, DM1 = -1.0D0 ) 3725C 3726 DIMENSION OTRS(KZYVR) 3727 DIMENSION DEN1(NASHDI,NASHDI) 3728 DIMENSION OPMAT(NORBT,NORBT), MJWOP(2,MAXWOP,8) 3729 DIMENSION WRK(*) 3730C 3731 IF (IPRRSP.GT.100) THEN 3732 WRITE(LUPRI,'(//A)') ' One electron matrix in ORBSX' 3733 CALL OUTPUT(OPMAT,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 3734 WRITE(LUPRI,'(//A)') 3735 * ' Sigma vector before one electron operator' 3736 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3737 END IF 3738C 3739C Compute offset for the k(k,l) elements 3740C 3741 NZCONF = MZCONF(IGRSYM) 3742 NYCONF = MZCONF(IGRSYM) + MZVAR(IGRSYM) 3743C 3744C distribute one electron matrix in the gradient: 3745C 3746 KSYM1 = 0 3747 DO 1300 IG = 1,MZWOPT(IGRSYM) 3748 K = MJWOP(1,IG,IGRSYM) 3749 L = MJWOP(2,IG,IGRSYM) 3750 KSYM = ISMO(K) 3751 LSYM = ISMO(L) 3752 IF( KSYM.NE.KSYM1 ) THEN 3753 KSYM1 = KSYM 3754 IORBK = IORB(KSYM) 3755 NASHK = NASH(KSYM) 3756 NISHK = NISH(KSYM) 3757 IASHK = IASH(KSYM) 3758 IIORBK= IIORB(KSYM) 3759 KXSYM = MULD2H(KSYM,ISYMDN) 3760 NORBKX= NORB(KXSYM) 3761 IORBKX= IORB(KXSYM) 3762 NASHKX= NASH(KXSYM) 3763 NISHKX= NISH(KXSYM) 3764 IASHKX= IASH(KXSYM) 3765 IORBL = IORB(LSYM) 3766 NASHL = NASH(LSYM) 3767 NISHL = NISH(LSYM) 3768 IASHL = IASH(LSYM) 3769 IIORBL= IIORB(LSYM) 3770 LXSYM = MULD2H(LSYM,ISYMDN) 3771 NORBLX= NORB(LXSYM) 3772 IORBLX= IORB(LXSYM) 3773 NASHLX= NASH(LXSYM) 3774 NISHLX= NISH(LXSYM) 3775 IASHLX= IASH(LXSYM) 3776 END IF 3777 ITYPK = IOBTYP(K) 3778 ITYPL = IOBTYP(L) 3779 IF ( ITYPK.EQ.JTINAC )THEN 3780C 3781C Do k(a,i) and k(i,a) 3782C 3783 IF (.NOT.TRPLET) THEN 3784 OTRS(NZCONF+IG) = OTRS(NZCONF+IG) 3785 * + D2 * OVLAP * OPMAT(L,K) 3786 OTRS(NYCONF+IG) = OTRS(NYCONF+IG) 3787 * - D2 * OVLAP * OPMAT(K,L) 3788 END IF 3789 IF ( ITYPL.EQ.JTACT .AND. NASHLX .GT. 0) THEN 3790C 3791C Do k(t,i) and k(i,t) 3792C 3793 NWL = ISW(L) - NISHT 3794 OTRS(NZCONF+IG) = OTRS(NZCONF+IG) 3795 * - DDOT(NASHLX,OPMAT(IORBLX+NISHLX+1,K),1, 3796 * DEN1(IASHLX+1,NWL),1) 3797 DO 730 IX = 1,NASHLX 3798 OTRS(NYCONF+IG) = OTRS(NYCONF+IG) 3799 * + OPMAT(K,IORBLX+NISHLX+IX)* 3800 * DEN1(NWL,IASHLX+IX) 3801 730 CONTINUE 3802 END IF 3803 ELSE 3804C 3805C We now know that K labels an active orbital 3806C 3807 IF (ITYPL.EQ.JTACT) THEN 3808C 3809C Do k(t,u) 3810C 3811 NWL = ISW(L) - NISHT 3812 NWK = ISW(K) - NISHT 3813 IF (NASHLX .GT. 0) THEN 3814 OTRS(NZCONF+IG) = OTRS(NZCONF+IG) 3815 * -DDOT(NASHLX,OPMAT(IORBLX+NISHLX+1,K),1, 3816 * DEN1(IASHLX+1,NWL),1) 3817 DO 740 IX = 1,NASHLX 3818 OTRS(NYCONF+IG) = OTRS(NYCONF+IG) 3819 * +OPMAT(K,IORBLX+NISHLX+IX)*DEN1(NWL,IASHLX+IX) 3820 740 CONTINUE 3821 END IF 3822 IF (NASHKX .GT. 0) THEN 3823 OTRS(NYCONF+IG) = OTRS(NYCONF+IG) 3824 * -DDOT(NASHKX,OPMAT(IORBKX+NISHKX+1,L),1, 3825 * DEN1(IASHKX+1,NWK),1) 3826 DO 750 IX = 1,NASHKX 3827 OTRS(NZCONF+IG) = OTRS(NZCONF+IG) 3828 * +OPMAT(L,IORBKX+NISHKX+IX)*DEN1(NWK,IASHKX+IX) 3829 750 CONTINUE 3830 END IF 3831 ELSE IF (NASHKX .GT. 0) THEN 3832C 3833C Do k(a,t) and k(t,a) 3834C 3835 NWK = ISW(K) - NISHT 3836 OTRS(NYCONF+IG) = OTRS(NYCONF+IG) 3837 * -DDOT(NASHKX,OPMAT(IORBKX+NISHKX+1,L),1, 3838 * DEN1(IASHKX+1,NWK),1) 3839 DO 760 IX = 1,NASHKX 3840 OTRS(NZCONF+IG) = OTRS(NZCONF+IG) 3841 * +OPMAT(L,IORBKX+NISHKX+IX)*DEN1(NWK,IASHKX+IX) 3842 760 CONTINUE 3843 ENDIF 3844 ENDIF 3845 1300 CONTINUE 3846C 3847 IF (IPRRSP.GT.100) THEN 3848 WRITE(LUPRI,'(//A)') ' Sigma vector from one electron operator' 3849 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3850 END IF 3851C 3852 RETURN 3853 END 3854C /* Deck consx */ 3855 SUBROUTINE CONSX(NSIM,KZYVR,IGRSYM,OPMAT,VEC,NZYVEC,NZCVEC, 3856 * ISYMV,NZCSTJ,ISYMJ,LREFST,OTRS,XINDX,WRK,LWRK) 3857C 3858C Purpose: 3859C Perform the configuration part of the linear transformation. 3860C Expressions: 3861C 3862C < j | zymat | 0(R) > 3863C - < 0(L) | zymat | j > 3864C 3865C 3866C The expression of the operator zymat between two arbitrary 3867C CSF's, in terms of zymat, is: 3868C 3869C <CSF1 | K | CSF2 > = sum(xy) zymat(x,y) D(x,y) 3870C 3871#include "implicit.h" 3872C 3873 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0, D2 = 2.0D0 ) 3874C 3875#include "maxorb.h" 3876#include "maxash.h" 3877#include "priunit.h" 3878#include "wrkrsp.h" 3879#include "infrsp.h" 3880#include "infopt.h" 3881#include "inforb.h" 3882#include "infind.h" 3883#include "infpri.h" 3884#include "rspprp.h" 3885#include "infhyp.h" 3886#include "infvar.h" 3887#include "qrinf.h" 3888#include "cbgetdis.h" 3889C 3890 DIMENSION OPMAT(NORBT,NORBT), VEC(NZYVEC) 3891 DIMENSION WRK(*), OTRS(KZYVR), XINDX(*) 3892C 3893 LOGICAL NOH2,IH8SM,LREFST, ITEST 3894 DATA ITEST /.FALSE./ 3895C 3896C Allocate work space for copying inactive part of one 3897C electron matrix 3898C 3899 KZYAC = 1 3900 KFREE = KZYAC + N2ASHX 3901 LFREE = LWRK - KFREE 3902 IF (LFREE.LT.0) CALL ERRWRK('CONSX',KFREE-1,LWRK) 3903C 3904 NOH2 = .TRUE. 3905 IH8SM = .TRUE. 3906 IF (TRPLET) THEN 3907 ISPIN1 = 1 3908 ISPIN2 = 0 3909 ELSE 3910 ISPIN1 = 0 3911 ISPIN2 = 0 3912 END IF 3913 IF ( LREFST ) THEN 3914 ISYMCI = IREFSY 3915 ELSE 3916 ISYMCI = MULD2H(IREFSY,ISYMV) 3917 END IF 3918C 3919 IF( IPRRSP .GT. 200) THEN 3920 WRITE(LUPRI,'(//A/)') ' Vector in CONSX' 3921 CALL OUTPUT(VEC,1,NZYVEC,1,NSIM,NZYVEC,NSIM,1,LUPRI) 3922 END IF 3923C 3924C 3925 IF( IPRRSP .GT. 200) THEN 3926 WRITE(LUPRI,'(//A/)') ' One electron matrix in CONSX' 3927 CALL OUTPUT(OPMAT,1,NORBT,1,NORBT,NORBT,NORBT, 3928 * 1,LUPRI) 3929 END IF 3930C 3931C Do <j | A | CR> 3932C =============== 3933C 3934C Copy active-active blocks of one electron matrix to WRK 3935C Put the sign of CR into the one electron matrix if .NOT. LREFST 3936C 3937 IF (. NOT. LREFST ) THEN 3938 SIGN = DM1 3939 ELSE 3940 SIGN = D1 3941 END IF 3942C 3943 OPINAC = D0 3944 IF (.NOT.TRPLET) THEN 3945 DO 105 IW = 1,NISHT 3946 IX = ISX(IW) 3947 OPINAC = OPINAC + OPMAT(IX,IX) 3948 105 CONTINUE 3949 OPINAC = SIGN * D2 * OPINAC 3950 END IF 3951C 3952C We need transposed matrix because 3953C CISIGD calculates <0|A|j>, and <j|A|0> = <0|At|j> 3954C 3955 DO 110 IW = 1, NASHT 3956 IX = ISX(NISHT + IW) 3957 DO 111 JW = 1, NASHT 3958 JX = ISX(NISHT + JW) 3959 IND = (IW-1) * NASHT + JW + KZYAC - 1 3960 WRK(IND) = SIGN * OPMAT(IX,JX) 3961111 CONTINUE 3962110 CONTINUE 3963C 3964 IF( IPRRSP .GT. 200) THEN 3965 WRITE(LUPRI,'(//A/)') 3966 * ' Active part of one electron matrix' 3967 CALL OUTPUT(WRK(KZYAC),1,NASHT,1,NASHT,NASHT,NASHT, 3968 * 1,LUPRI) 3969 END IF 3970C 3971 CALL CISIGD(ISYMCI,ISYMJ,NZCVEC,NZCSTJ,VEC, 3972 * OTRS,WRK(KZYAC),DUMMY, 3973 * NOH2,IH8SM,XINDX,ISPIN1,ISPIN2,WRK(KFREE), 3974 * LFREE) 3975C 3976 IF (ISYMCI .EQ. ISYMJ .AND..NOT.TRPLET) THEN 3977 IF ( NISHT .GT. 0 ) 3978 * CALL DAXPY(NZCVEC,OPINAC,VEC,1, 3979 * OTRS,1) 3980 END IF 3981C 3982 IF( ITEST .AND. .NOT.TRPLET ) THEN 3983 IF ( LREFST .AND. (ISYMCI .EQ .ISYMJ) ) THEN 3984 WRITE(LUPRI,'(//A,/A)') 3985 * ' Sigma vector after Z part of linear transformation', 3986 * ' Reference state included' 3987 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3988 T1 = DDOT(NZCVEC,VEC,1,OTRS,1) 3989 WRITE(LUPRI,'(//A)') ' Reference state projected out' 3990 WRITE(LUPRI,'(A,F14.7)') ' Average value T1 =', T1 3991 CALL DAXPY(NZCVEC,(-T1),VEC,1,OTRS,1) 3992 END IF 3993 END IF 3994C 3995 IF ( IPRRSP .GT. 200 ) THEN 3996 WRITE(LUPRI,'(//A)') 3997 * ' Sigma vector after Z part of linear transformation' 3998 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 3999 END IF 4000C 4001C Do -<CL | A | j> 4002C =============== 4003C 4004C Copy active-active blocks of one electron matrix to WRK 4005C Put the minus sign into the operator 4006C 4007 IF ( LREFST ) OPINAC = -OPINAC 4008 DO 210 IW = 1, NASHT 4009 IX = ISX(NISHT + IW) 4010 DO 211 JW = 1, NASHT 4011 JX = ISX(NISHT + JW) 4012 IND = (IW-1) * NASHT + JW + KZYAC - 1 4013 WRK(IND) = DM1 * OPMAT(JX,IX) 4014211 CONTINUE 4015210 CONTINUE 4016C 4017 IF( IPRRSP .GT. 200) THEN 4018 WRITE(LUPRI,'(//A/)') 4019 * ' Transposed active part of one electron matrix' 4020 CALL OUTPUT(WRK(KZYAC),1,NASHT,1,NASHT,NASHT,NASHT, 4021 * 1,LUPRI) 4022 END IF 4023C 4024 IF( LREFST ) THEN 4025 ICOFF = 1 4026 ELSE 4027 ICOFF = MZCONF(ISYMV) + MZWOPT(ISYMV) + 1 4028 ENDIF 4029 ISOFF = MZCONF(IGRSYM) + MZWOPT(IGRSYM) + 1 4030C 4031 CALL CISIGD(ISYMCI,ISYMJ,NZCVEC,NZCSTJ,VEC(ICOFF), 4032 * OTRS(ISOFF),WRK(KZYAC),DUMMY, 4033 * NOH2,IH8SM,XINDX,ISPIN1,ISPIN2,WRK(KFREE), 4034 * LFREE) 4035C 4036 IF (ISYMCI .EQ. ISYMJ .AND. .NOT.TRPLET) THEN 4037 IF ( NISHT .GT. 0 ) 4038 * CALL DAXPY(NZCVEC,OPINAC,VEC(ICOFF),1, 4039 * OTRS(ISOFF),1) 4040 END IF 4041C 4042 IF( ITEST .AND. .NOT.TRPLET ) THEN 4043 IF ((IREFSY .EQ .ISYMJ).AND. LREFST ) THEN 4044 WRITE(LUPRI,'(//A,/A)') 4045 * ' Sigma vector after Y part of linear transformation', 4046 * ' Reference state included' 4047 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 4048 T1 = DDOT(NZCVEC,VEC(ICOFF),1, 4049 * OTRS(ISOFF),1) 4050 WRITE(LUPRI,'(//A)') ' Reference state projected out' 4051 WRITE(LUPRI,'(A,F14.7)') ' Average value T1 =', T1 4052 CALL DAXPY(NZCVEC,(-T1),VEC(ICOFF),1, 4053 * OTRS(ISOFF),1) 4054 END IF 4055 END IF 4056C 4057 IF ( IPRRSP .GT. 200 ) THEN 4058 WRITE(LUPRI,'(//A)') 4059 * ' Sigma vector after Y part of linear transformation' 4060 CALL OUTPUT(OTRS,1,KZYVR,1,NSIM,KZYVR,NSIM,1,LUPRI) 4061 END IF 4062C 4063 RETURN 4064 END 4065C /* Deck rsph1 */ 4066 SUBROUTINE RSPH1(H1,CMO,WRK,LWRK) 4067C 4068C Get the one-index transformed one electron integrals and put 4069C them in H1 4070C 4071#include "implicit.h" 4072 DIMENSION H1(NORBT,NORBT), CMO(NCMOT), WRK(LWRK) 4073C 4074C Used from common blocks: 4075C INFORB : NORBT,NBAST,... 4076C INFRSP : IPRRSP 4077C 4078#include "priunit.h" 4079#include "inforb.h" 4080#include "infrsp.h" 4081C 4082 LNEED = N2BASX + NNBAST 4083 IF (LNEED .GT. LWRK) THEN 4084 CALL QUIT('Insufficient WRK space in RSPH1') 4085 END IF 4086 KH1AO = 1 4087 KWRK1 = KH1AO + N2BASX 4088 LWRK1 = LWRK + 1 - KWRK1 4089 IF (LWRK1.LT.0) CALL ERRWRK('RSPH1',KWRK1-1,LWRK) 4090C use SIRH1, not RDONEL, in order get field dependent terms added 4091C CALL RDONEL('ONEHAMIL',.TRUE.,WRK,NNBAST) 4092 CALL SIRH1(WRK(KH1AO),WRK(KWRK1),LWRK1) 4093 CALL PKSYM1(WRK(KWRK1),WRK,NBAS,NSYM,-1) 4094 CALL DSPTSI(NBAST,WRK(KWRK1),WRK(KH1AO)) 4095C 4096 CALL DZERO(H1,N2ORBX) 4097 DO 600 ISYM=1,NSYM 4098 IF(NORB(ISYM).LE.0)GO TO 600 4099 CALL UTHV(CMO(ICMO(ISYM)+1),WRK(KH1AO),CMO(ICMO(ISYM)+1), 4100 & ISYM,ISYM,NBAS(ISYM),NBAS(ISYM), 4101 & H1,WRK(KWRK1)) 4102C CALL UTHV(U,PRPAO,V,ISYM,JSYM,NBASI,NBASJ,PRPMO,WRK) 4103 600 CONTINUE 4104 IF (IPRRSP.GT.1000) THEN 4105 WRITE(LUPRI,'(/A)')' One-electron Hamiltonian in MO basis' 4106 CALL OUTPUT(H1,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 4107 END IF 4108C 4109 RETURN 4110 END 4111 SUBROUTINE QRGP(WORD,GP,CMO,XINDX,ANTSYM_GP,WORK,LWORK) 4112#include "implicit.h" 4113 CHARACTER*8 WORD 4114 DIMENSION GP(*), CMO(*), XINDX(*), WORK(LWORK) 4115#include "inforb.h" 4116#include "infrsp.h" 4117#include "infvar.h" 4118#include "rspprp.h" 4119#include "inflr.h" 4120#include "wrkrsp.h" 4121#include "priunit.h" 4122#include "qrinf.h" 4123#include "inflin.h" 4124#include "infrank.h" 4125 LOGICAL LREFST, LORB, LCON, TRPSAVE 4126 PARAMETER (D1 = 1.0D0) 4127C 4128 CALL QENTER('QRGP') 4129C 4130 KWORK = 1 4131 KFREE = KWORK 4132 LFREE = LWORK 4133C 4134C Operator 4135C 4136 CALL MEMGET('REAL',KOP,N2ORBX,WORK,KFREE,LFREE) 4137 KSYMP=KSYMOP 4138 IPR=IPRRSP 4139 CALL PRPGET(WORD,CMO,WORK(KOP),KSYMP,ANTSYM,WORK(KFREE),LFREE,IPR) 4140 IOPSYM=KSYMOP 4141 IOPSPI=OPRANK(INDPRP(WORD)) 4142C 4143C Whereas ANTSYM from PRPGET refers to matrix symmetry, ANTSYM_GP will 4144C from here on refer to the antisymmetry of the response property vector, 4145C thus we change sign on it. 4146C 4147 ANTSYM_GP = -ANTSYM 4148C 4149C Density 4150C 4151 IGRSYM=KSYMOP 4152 IF (TRPLET) THEN 4153 IGRSPI=1 4154 ELSE 4155 IGRSPI=0 4156 END IF 4157 CALL MEMGET('REAL',KD,N2ASHX,WORK,KFREE,LFREE) 4158 IF (TDHF) THEN 4159 CALL DUNIT(WORK(KD),NASHT) 4160 CALL MEMGET('REAL',KCREF,0,WORK,KFREE,LFREE) 4161 ELSE 4162 CALL MEMGET('REAL',KCREF,MZCONF(1),WORK,KFREE,LFREE) 4163 CALL MEMGET('REAL',KP,1,WORK,KFREE,LFREE) 4164 CALL GETREF(WORK(KCREF),MZCONF(1)) 4165 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 4166 & WORK(KCREF),WORK(KCREF),WORK(KD),WORK(KP), 4167 & IGRSPI,IOPSPI,.FALSE.,.TRUE.,XINDX,WORK,KFREE,LFREE) 4168 END IF 4169C 4170C Gradient 4171C 4172 LGP = KZYVAR 4173 IGPSYM = IGRSYM 4174 IDSYM = 1 4175 LREFST = .TRUE. 4176 LORB = KZWOPT.GT.0 4177 IF (IGRSYM.EQ.1) THEN 4178 LCON = KZCONF.GT.1 4179 ELSE 4180 LCON = KZCONF.GT.0 4181 END IF 4182 CALL MEMGET('INTE',KMJWOP,2*MAXWOP*8,WORK,KFREE,LFREE) 4183 CALL SETZY(WORK(KMJWOP)) 4184 CALL DZERO(GP,LGP) 4185 TRPSAVE=TRPLET 4186 CALL RSP1GR(1,LGP,IOPSYM,IOPSPI,IGRSYM,IGRSPI, 4187 & IGPSYM,GP, 4188 & WORK(KCREF),MZCONF(1),MZCONF(1),D1,IDSYM,WORK(KD),WORK(KOP), 4189 & XINDX,WORK(KMJWOP),WORK(KFREE),LFREE,LORB,LCON,.TRUE.) 4190 TRPLET=TRPSAVE 4191#ifdef DEBUG 4192 CALL HEADER('QRGP:Operator '//WORD,0) 4193 CALL OUTPUT(WORK(KOP),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 4194 CALL HEADER('QRGP:One electron density',0) 4195 CALL OUTPUT(WORK(KD),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 4196 CALL HEADER('QRGP:Property vector',0) 4197 CALL OUTPUT(GP,1,LGP/2,1,2,LGP/2,2,1,LUPRI) 4198#endif 4199 CALL MEMREL('QRGP',WORK,KWORK,KWORK,KFREE,LFREE) 4200 CALL QEXIT('QRGP') 4201 END 4202