1!#define DEBUG_RSPOLI 2#undef DEBUG_RSPOLI 3! 4! Dalton, a molecular electronic structure program 5! Copyright (C) by the authors of Dalton. 6! 7! This program is free software; you can redistribute it and/or 8! modify it under the terms of the GNU Lesser General Public 9! License version 2.1 as published by the Free Software Foundation. 10! 11! This program is distributed in the hope that it will be useful, 12! but WITHOUT ANY WARRANTY; without even the implied warranty of 13! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 14! Lesser General Public License for more details. 15! 16! If a copy of the GNU LGPL v2.1 was not distributed with this 17! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 18! 19! 20C 21C FILE: rsp/rspoli.F ("Orbital LInear transformation") 22C 23!======================================================================== 24!990427-hjaaj: use if DXAO is symmetric or antisymmetric in IFCTYP 25!961228-kr Transfered SOPPA particle-hole matrices in paramater list to RSPOLI 26!950513-kr Added comdeck INFVAR and and changed to ISYMDM()=JWOPSY 27!940708-hjaaj: SOPPA changes 28!940602-hjaaj: do not call QONEDI and QTD for nasht.eq.1 29!931124-martin packer: added call to HRPA which adds A(2) and B(2) matrices to FCX 30!920721-hinne hettema: implemented RSPSUP (super symmetry averaging) 31!======================================================================== 32 33C /* Deck rspoli */ 34 SUBROUTINE RSPOLI(NCSIM,NOSIM,UDV,ZYCVEC,LZYCVEC,FC,FV,PVX,ZYMAT, 35 * FCX,FVX,QAX,QBX, FVTD,QATD,QBTD, EVECS, 36 * XINDX,CMO,WRK,LWRK) 37 38C 39C Copyright 13 FEB 1986 Poul Joergensen and Hans Joergen Aa. Jensen 40C 41C 21-1-1992: Changed for Direct RPA (KEYWORD: DIRFCK), HA. 42C 21-7-1992: Supersymmetry averaging added. Hinne Hettema 43C 28-2-1995: Changed for Direct RPA with AO integrals read 44C from disk, NOITRA=.TRUE. P&D /is not included?/hjaaj 45C Oct-2007-hjaaj: added LZYCVEC because in lagrang.F RSPOLI is called 46C with ZYCVEC=CREF, and for some RHF cases KZCONF=0 and NCREF=1 47C 48C DFT modifications T. Helgaker 49C 50C PURPOSE: 51C CALCULATE ORBITAL PART OF THE LINEAR TRANSFORMATIONS E[2]*X 52C 53C FLOW: 54C 1) READ IN MULLIKEN INTEGRAL DISTRIBUTION AND ADD CONTRIBUTIONS 55C TO FOCK CORE , FOCK VALENCE , QA AND QB MATRICES 56C 2) READ IN DIRAC INTEGRAL DISTRIBUTION AND ADD CONTRIBUTIONS 57C TO QAX AND QBX MATRICES 58C 3) DISTRIBUTE FOCK CORE , FOCK VALENCE , QA AND QB MATRICES 59C INTO ORBITAL PART OF LINEAR TRANSFORMED VECTORS 60C 61C Naming conventions: 62C FCX, FVX, ... are FC, FV, ... with one-index transformed integrals 63C FVTD, ... are FV calculated with Transition dens.mat. 64C 65#include "implicit.h" 66C 67 DIMENSION UDV(NASHDI,*),ZYCVEC(*),FC(*),FV(*),PVX(*) 68 DIMENSION ZYMAT(N2ORBX,*),FCX(NORBT,NORBT,*) 69 DIMENSION FVX(NORBT,NORBT,*),QAX(NORBT,NASHDI,*) 70 DIMENSION QBX(NORBT,NASHDI,*),FVTD(NORBT,NORBT,*) 71 DIMENSION QATD(NORBT,NASHDI,*),QBTD(NORBT,NASHDI,*) 72 DIMENSION EVECS(*) 73 DIMENSION XINDX(*),WRK(*),CMO(*) 74C 75#include "thrzer.h" 76#include "dummy.h" 77 78#include "cbihr2.h" 79C 80C INFDIM : NASHDI 81C INFINP : DIRFCK, HSROHF, DODFT, ? 82C INFRSP : ? 83C INFSOP : A2EXIST,KABSAD,KABTAD,... 84C 85#include "gnrinf.h" 86#include "maxorb.h" 87#include "maxash.h" 88#include "priunit.h" 89#include "inforb.h" 90#include "infind.h" 91#include "infinp.h" 92#include "infdim.h" 93#include "infpri.h" 94#include "infrsp.h" 95#include "infvar.h" 96#include "inftra.h" 97#include "orbtypdef.h" 98#include "wrkrsp.h" 99#include "infsop.h" 100#include "dftcom.h" 101 102CSONIA SONIA SONIA control direct for kappabar of CCSD(T) 103#include "grdccpt.h" 104CSONIA SONIA SONIA end control direct for kappabar of CCSD(T) 105 106C 107C -- local constants 108C 109 PARAMETER ( D1 = 1D0, DP5 = 0.5D0, D2 = 2.0D0, DM1 = -1.0D0 ) 110 LOGICAL DFTADX 111 LOGICAL DOFXC, DOFXV, DOFTC, DOFTV 112 integer, allocatable :: ISYMDM(:), IFCTYP(:) 113 DIMENSION NEEDMU(-4:6),NEEDDI(-4:6) 114 115 CALL QENTER('RSPOLI') 116C 117 IF (IPRRSP .GT. 5) THEN 118 write(lupri,'(//A,2I10/A,5L10/A,L10,I10)') 119 & 'Output from "RSPOLI". NCSIM,NOSIM=',NCSIM,NOSIM, 120 & 'TRPLET,HSROHF,DOHSRDFT,DOMCSRDFT,DFTADD', 121 & TRPLET,HSROHF,DOHFSRDFT,DOMCSRDFT,DFTADD, 122 & 'TDHF,NASHT',TDHF,NASHT 123 END IF 124C 125 NEEDMU(-4:6) = 0 126 NEEDDI(-4:6) = 0 127 IF (SOPPA) THEN 128 NEEDMU(1:4) = 1 129C exclude virt-virt contributions to B(2) matrix 130C for .OPTORB iterations from ORPCTL 131 IF (NOSIM .GT. 0 .AND. KZCONF .GT. 0) NEEDMU(6) = 1 132 NEEDDI(1) = 1 133 ELSE IF (DIRFCK) THEN ! only act-act needed for FQ if MCSCF/MC-srDFT (KZCONF .gt. 0) 134 NEEDMU(3) = 1 135 NEEDDI(3) = 1 136 ELSE 137 NEEDMU(1:3) = 1 138 NEEDDI(1:3) = 1 139 END IF 140 141 DOFXC = NISHT .GT. 0 .OR. DODFT .OR. DOHFSRDFT .OR. DOMCSRDFT 142 DOFXV = NASHT .GT. 0 143 DOFTC = DOMCSRDFT 144 DOFTV = NASHT .GT. 0 ! i.e. .false. for SOPPA 145 146 147C 148C ALLOCATE WORK SPACE FOR DTV and PTVD and SOPPA 149C 150 KFRSAV = 1 151 KFREE = KFRSAV 152 LFREE = LWRK 153 CALL MEMGET2('REAL','DTV',KDTV ,NCSIM*N2ASHX,WRK,KFREE,LFREE) 154 CALL MEMGET2('REAL','PTVD',KPTVD,NCSIM*N2ASHX*N2ASHX, 155 & WRK,KFREE,LFREE) 156 IF (KZCONF.GT.0) THEN 157 LH2XAC = NOSIM*N2ASHX*NNASHX 158 ELSE 159 LH2XAC = 0 160 END IF 161 CALL MEMGET2('REAL','H2XAC',KH2XAC,LH2XAC,WRK,KFREE,LFREE) 162 IF (SOPPA) THEN 163 CALL MEMGET2('REAL','H2XP',KH2XP ,N2ORBX,WRK,KFREE,LFREE) 164 CALL MEMGET2('REAL','COEUN',KCOEUN,N2ORBX,WRK,KFREE,LFREE) 165 CALL MEMGET2('REAL','TZYMT',KTZYMT,NOSIM*N2ORBX, 166 & WRK,KFREE,LFREE) 167 ELSE 168 CALL MEMGET2('REAL','H2XP',KH2XP ,0,WRK,KFREE,LFREE) 169 CALL MEMGET2('REAL','COEUN',KCOEUN,0,WRK,KFREE,LFREE) 170 CALL MEMGET2('REAL','TZYMT',KTZYMT,0,WRK,KFREE,LFREE) 171 ENDIF 172 KWRK0 = KFREE 173C 174C CALCULATE TRANSITION DENSITY MATRIX 175C 176 IF ( (NCSIM.GT.0) .AND. (.NOT.RSPCI) .AND. (.NOT.SOPPA) ) THEN 177 CALL MEMGET2('REAL','CREF',KCREF,NCREF,WRK,KFREE,LFREE) 178 IF (DOMCSRDFT) KWRK0 = KFREE ! must keep CREF for later 179 CALL GETREF(WRK(KCREF),NCREF) 180 IF (TRPLET) THEN 181 ISPIN1 = 1 182 ISPIN2 = 0 183 ELSE 184 ISPIN1 = 0 185 ISPIN2 = 0 186 END IF 187 CALL RSPTDM(NCSIM,IREFSY,KSYMST,NCREF,LZYCVEC,WRK(KCREF), 188 * ZYCVEC,WRK(KDTV),WRK(KPTVD), 189 * ISPIN1,ISPIN2,.TRUE.,.FALSE., 190 * XINDX,WRK,KFREE,LFREE) 191C CALL RSPTDM(NCSIM,ILRESY,IRSYM,NCLREF,NCRDIM,CLREF, 192C * CR, RHO1,RHO2, ISPIN1,ISPIN2,TDM,NORHO2, 193C * XNDXCI,WORK,KFREE,LFREE) 194 CALL MEMREL('RSPOLI.tdm',WRK,KFRSAV,KWRK0,KFREE,LFREE) 195 END IF 196 197C 198C INITIALIZE ONE INDEX TRANSFORMED (X) AND TRANSITION DENSITY (TD) 199C Q, FC AND FV MATRICES 200C 201 IF (NOSIM.GT.0) THEN 202 IF (SOPPA) THEN 203C ... A(2) matrix is saved permanently in XINDX 204C Initialisation of XINDX is done in SET2SOPPA now. 205C IF (.NOT.A2EXIST) CALL DZERO(XINDX,N2ORBX) 206 DO I = 1,NOSIM 207 JTZYMT = KTZYMT + (I-1)*N2ORBX 208 CALL MTRSP(NORBT,NORBT,ZYMAT(1,I),NORBT,WRK(JTZYMT),NORBT) 209C CALL MTRSP(NROWA,NCOLA,A,NRDIMA,B,NRDIMB) 210 END DO 211 END IF 212 CALL DZERO(FCX,N2ORBX*NOSIM) 213 IF (NASHT.GT.0) THEN 214 CALL DZERO(FVX,N2ORBX*NOSIM) 215 NTOT = NORBT*NASHT*NOSIM 216 CALL DZERO(QAX,NTOT) 217 CALL DZERO(QBX,NTOT) 218 IF (KZCONF.GT.0 ) THEN 219 NTOT = NOSIM*N2ASHX*NNASHX 220 CALL DZERO(WRK(KH2XAC),NTOT) 221 END IF 222 END IF 223 END IF 224 IF (NCSIM.GT.0) THEN 225 NTOT = N2ORBX * NCSIM 226 IF (DOMCSRDFT) NTOT = 2*NTOT 227C ... we also need FCTD for MC-srDFT 228 CALL DZERO(FVTD,NTOT) 229 NTOT = NORBT*NASHT*NCSIM 230 IF (NTOT .GT. 0) THEN 231C ... is automatically zero for SOPPA (as NASHT.eq.0 then) 232 CALL DZERO(QATD,NTOT) 233 CALL DZERO(QBTD,NTOT) 234 END IF 235 END IF 236C 237C If DIRFCK compute Fock matrices in AO basis 238C If also single determinant, then nothing to do in MO part 239C 240 !IF (DIRFCK .AND. (TDHF .OR. DODFT)) GO TO 1000 241 !IF (DIRFCK .AND. KZCONF .EQ. 0) GO TO 1000 242 IF (HSROHF .OR. DODFT) GO TO 1000 243 ! ... KS-DFT only works with direct code, even if DIRFCK is false /hjaaj 244 IF (DIRFCK .AND. TDHF .AND. NASHT.EQ.1) GO TO 1000 ! .OPEN SHELL HF 245 IF (DIRFCK .AND. KZCONF.EQ.0 .AND. NASHT.EQ.0) GO TO 1000 246 ! new test Dec. 2011, also valid for SOPPA .and. DIRFCK 247 ! where TDHF true (KZCONF > 0 and NASHT == 0 for SOPPA) /hjaaj 248 249CSONIA SONIA SONIA 250CSONIA SONIA SONIA control direct for kappabar of CCSD(T) 251CSONIA SONIA SONIA 252 IF (LGRDCCPT) THEN 253 write(lupri,*)'Warning RSPOLI: LGRDCCPT = ', LGRDCCPT 254 GO TO 1000 255 END IF 256CSONIA SONIA SONIA 257CSONIA SONIA SONIA end control direct for kappabar of CCSD(T) 258CSONIA SONIA SONIA 259C 260C If not DIRFCK, i.e. 261C calculate Fock matrices etc. in MO basis, then 262C ALLOCATE WORK SPACE FOR MULLIKEN and DIRAC DISTRIBUTIONs 263C to Fock matrices 264C 265 IF (.NOT.DIRFCK) THEN 266 CALL MEMGET2('REAL','DENA',KDENA ,NOSIM*NORBT*NASHT, 267 & WRK,KFREE,LFREE) 268 CALL MEMGET2('REAL','DENB',KDENB ,NOSIM*NORBT*NASHT, 269 & WRK,KFREE,LFREE) 270 CALL MEMGET2('REAL','FCOCO',KFCOCO,NOSIM*N2ORBX, 271 & WRK,KFREE,LFREE) 272 CALL DZERO(WRK(KFCOCO),N2ORBX*NOSIM) 273 IF (NASHT.GT.0) THEN 274 CALL MEMGET2('REAL','FVOCO',KFVOCO,NOSIM*N2ORBX, 275 & WRK,KFREE,LFREE) 276 CALL DZERO(WRK(KFVOCO),N2ORBX*NOSIM) 277 ELSE 278 CALL MEMGET2('REAL','FVOCO',KFVOCO,0,WRK,KFREE,LFREE) 279 END IF 280 IF ((NASHT .GT. 0) .AND. (NOSIM .GT.0)) THEN 281 CALL RSPTR1(NOSIM,UDV,ZYMAT,WRK(KDENA),WRK(KDENB)) 282 END IF 283 ELSE 284 CALL MEMGET2('REAL','DENA' ,KDENA ,0,WRK,KFREE,LFREE) 285 CALL MEMGET2('REAL','DENB' ,KDENB ,0,WRK,KFREE,LFREE) 286 CALL MEMGET2('REAL','FCOCO',KFCOCO,0,WRK,KFREE,LFREE) 287 CALL MEMGET2('REAL','FVOCO',KFVOCO,0,WRK,KFREE,LFREE) 288 END IF 289 CALL MEMGET2('REAL','H2' ,KH2 ,N2ORBX,WRK,KFREE,LFREE) 290 CALL MEMGET2('REAL','H2X',KH2X,N2ORBX,WRK,KFREE,LFREE) 291 KWRK1 = KFREE 292C 293C setup for reading Mulliken MO integrals ... 294C NEEDED DISTRIBUTIONS DEFINED IN NEEDMU(-4:6) 295C 296 CALL MEMGET2('REAL','PVCD',KPVCD , N2ASHX,WRK,KFREE,LFREE) 297 IDIST = 0 298 90 CALL NXTH2M(IC,ID,WRK(KH2),NEEDMU,WRK,KFREE,LFREE,IDIST) 299 IF (IDIST .LT. 0) GO TO 95 300C ... IF IDIST.LT.0 NO MORE DISTRIBUTIONS 301 KWRK2 = KFREE 302 LWRK2 = LFREE 303C (KFREE,LFREE are updated inside NXTH2M) 304 ICDSYM= MULD2H(ISMO(IC),ISMO(ID)) 305 ICW = ISW(IC) 306 IDW = ISW(ID) 307 IDI = ID 308 ICI = IC 309 ICIW = ICW 310 IDIW = IDW 311C 312C find distribution type ITYPCD = 313C 1:inactive-inactive 2:active-inactive 3:active-active 314C 4:secondary-inactive 5:secondary-active 6:secondary-secondary 315C We do not need type 6 (except for SOPPA). 316C 317 ITYPC = IOBTYP(IC) 318 ITYPD = IOBTYP(ID) 319 ITYPCD = IDBTYP(ITYPC,ITYPD) 320C 321C ORDER INDICES OF C AND D SUCH THAT C < D 322C 323 IF ( ICW.GT.IDW ) THEN 324 ICIW = IDW 325 IDIW = ICW 326 ICI = ID 327 IDI = IC 328 ISWAP = ITYPC 329 ITYPC = ITYPD 330 ITYPD = ISWAP 331 ENDIF 332 NCIW = ICIW - NISHT 333 NDIW = IDIW - NISHT 334C 335 336 IF (NOSIM.GT.0) THEN 337C 338C CONTRIBUTIONS TO FCX AND FVX 339C 340 IF (.NOT.DIRFCK .AND. ITYPCD .LE. 3) 341 & CALL FONEMU(NOSIM,ICI,IDI,WRK(KH2), 342 & FCX,FVX,ZYMAT,WRK(KDENA),WRK(KDENB), 343 & WRK(KWRK2),LWRK2) 344C CALL FONEMU(NSIM,ICI1,IDI1,H2, 345C * FCOEX,FVOEX,ZYMAT,DENA,DENB,WRK,LWRK) 346C 347Cmjp CONTRIBTUIONS FORM THE HRPA TERMS 348 IF (SOPPA) THEN 349 CALL HRPA(FCX,UDV,PVX,WRK(KH2),ZYMAT,WRK(KTZYMT), 350 * WRK(KH2X),WRK(KH2XP),WRK(KCOEUN),ICI,IDI,ICDSYM, 351 * ITYPCD,NOSIM,XINDX(KIADR1),WRK(KWRK2),LWRK2) 352 IF (ITYPCD .EQ. 4) THEN 353 IF (.NOT.HIRPA .AND. KZCONF .GT. 0) THEN 354C ... skip SOPH2X for .OPTORB iterations 355 CALL SOPH2X(EVECS(1+NCSIM*KZYVAR),ZYMAT,WRK(KTZYMT), 356 * WRK(KH2),ICI,IDI,ICDSYM,NOSIM,XINDX(KABSAD), 357 & XINDX(KABTAD),XINDX(KIJSAD),XINDX(KIJTAD), 358 & XINDX(KIJ1AD),XINDX(KIJ2AD),XINDX(KIJ3AD), 359 * WRK(KWRK2),LWRK2) 360 END IF 361C Construct the A(2) matrix explicitly (but only once) 362C for the symmetrization of A(2)b 363 IF (.NOT.A2EXIST) THEN 364 CALL HRPAA2(PVX,WRK(KH2),XINDX(KAB2),WRK(KCOEUN), 365 & XINDX(KIADR1),ICI,IDI,ICDSYM) 366 END IF 367 END IF 368 END IF 369 IF ((NASHT.GT.1).AND.(ITYPCD.EQ.3)) THEN 370C 371C CONTRIBUTIONS TO QAX AND QBX 372C 373 CALL QONEMU(NOSIM,NCIW,NDIW,ICDSYM, 374 * QAX,QBX,ZYMAT,WRK(KH2),WRK(KH2X),WRK(KH2XAC), 375 * PVX,WRK(KPVCD),WRK(KWRK2),LWRK2) 376C CALL QONEMU(NOSIM,NCIW,NDIW,ICDSYM,QAX,QBX, 377C * ZYMAT,H2,H2X,H2XAC,PVX,PVCD,WRK,LWRK) 378 END IF 379 END IF 380 IF (NCSIM.GT.0) THEN 381 IF (SOPPA) THEN 382 IF (.NOT.HIRPA .AND. ITYPCD.EQ.4) THEN 383 CALL SOPPAF(FVTD,WRK(KH2),ZYCVEC,ICI,IDI,ICDSYM,NCSIM, 384 * XINDX,WRK(KWRK2),LWRK2) 385 END IF 386 ELSE 387 IF (.NOT. DIRFCK .AND. 388 & (ITYPCD.EQ.2 .OR. ITYPCD.EQ.3 .OR. ITYPCD.EQ.5) ) THEN 389C 390C CONTRIBUTIONS TO FVTD 391C 392 CALL FTDMU(NCSIM,ICI,IDI,FVTD,WRK(KDTV),WRK(KH2)) 393C CALL FTDMU(NSIM,ICI1,IDI1,FVTD,DTV,WRK(KH2)) 394 END IF 395 IF (ITYPCD.EQ.3) THEN 396C 397C CONTRIBUTIONS TO QATD AND QBTD 398C 399 CALL QTD(NCSIM,NCIW,NDIW,ICDSYM,QATD,QBTD,WRK(KH2), 400 * WRK(KPTVD),WRK(KPVCD),WRK(KWRK2),LWRK2) 401C CALL QTD(NCSIM,NCIW,NDIW,ICDSYM,QATD,QBTD,H2,PTVD,PVDEN, 402C * WRK,LWRK) 403 END IF 404 END IF 405 END IF 406 GO TO 90 407C 408C ALL MULLIKEN DISTRIBUTIONS HAVE BEEN READ IN 409C 410 95 CONTINUE 411 IF ((IPRRSP.GT.50).AND.(NOSIM.GT.0).AND.(NASHT.GT.1)) THEN 412 DO 1001 IOSIM = 1,NOSIM 413 WRITE(LUPRI,'(/A)') ' MULLIKEN DISTRIBUTION CONTRIBUTION' 414 WRITE(LUPRI,'(/A,I5,A)') 415 * ' QAX FOR',IOSIM,' ORBITAL TRIAL VECTOR' 416 CALL OUTPUT(QAX(1,1,IOSIM), 417 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 418 WRITE(LUPRI,'(/A,I5,A)') 419 * ' QBX FOR',IOSIM,' ORBITAL TRIAL VECTOR' 420 CALL OUTPUT(QBX(1,1,IOSIM), 421 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 422 1001 CONTINUE 423 END IF 424 IF ((IPRRSP.GT.50).AND.(NOSIM.GT.0).AND. .NOT.DIRFCK) THEN 425 DO 1002 IOSIM = 1,NOSIM 426 WRITE(LUPRI,'(/2A,I5)')' FCOEX MULLIKEN DISTRIBUTION CONTRB' 427 * ,' FOR VECTOR ',IOSIM 428 CALL OUTPUT(FCX(1,1,IOSIM), 429 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 430 IF (NASHT .GT. 0) THEN 431 WRITE(LUPRI,'(/2A,I5)')' FVOEX MULLIKEN DISTRIBUTION CONTRB' 432 * ,' FOR VECTOR ',IOSIM 433 CALL OUTPUT(FVX(1,1,IOSIM), 434 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 435 END IF 436 1002 CONTINUE 437 END IF 438 IF ((IPRRSP.GT.50).AND.(NCSIM.GT.0).AND.(NASHT.GT.1)) THEN 439 DO 1003 ICSIM = 1,NCSIM 440 WRITE(LUPRI,'(/A)') ' MULLIKEN DISTRIBUTION CONTRIBUTION' 441 WRITE(LUPRI,'(/A,I5,A)') 442 * ' QATD FOR',ICSIM,' CONFIGURATION TRIAL VECTOR' 443 CALL OUTPUT(QATD(1,1,ICSIM), 444 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 445 WRITE(LUPRI,'(/A,I5,A)') 446 * ' QBTD FOR',ICSIM,' CONFIGURATION TRIAL VECTOR' 447 CALL OUTPUT(QBTD(1,1,ICSIM), 448 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 449 1003 CONTINUE 450 END IF 451 IF ((IPRRSP.GT.50).AND.(NCSIM.GT.0).AND. .NOT.DIRFCK) THEN 452 DO 1004 ICSIM = 1,NCSIM 453 WRITE(LUPRI,'(/2A,I5)')' FVTD MULLIKEN DISTRIBUTION CONTRB' 454 * ,' FOR VECTOR ',ICSIM 455 CALL OUTPUT(FVTD(1,1,ICSIM), 456 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 457 1004 CONTINUE 458 END IF 459 CALL MEMREL('RSPOLI.h2m',WRK,KFRSAV,KWRK1,KFREE,LFREE) 460C 461C ********************************************************** 462C 463C ADD CONTRIBUTIONS FROM DIRAC INTEGRAL DISTRIBUTIONS 464C 465C ********************************************************* 466C 467 LUINTD = -1 ! file MODRCINT will be opened in NXTH2D, if needed 468 IF (SOPPA .AND. DIRFCK) GO TO 1095 469C 470C ALLOCATE WORK SPACE FOR DIRAC DISTRIBUTION 471C 472 CALL MEMGET2('REAL','PVCD2',KPVCD2,N2ASHX,WRK,KFREE,LFREE) 473 CALL MEMGET2('REAL','PVCD3',KPVCD3,N2ASHX,WRK,KFREE,LFREE) 474 CALL MEMGET2('REAL','PVDC2',KPVDC2,N2ASHX,WRK,KFREE,LFREE) 475 CALL MEMGET2('REAL','PVDC3',KPVDC3,N2ASHX,WRK,KFREE,LFREE) 476C 477 IDIST = 0 478 1090 CALL NXTH2D(IC,ID,WRK(KH2),NEEDDI,LUINTD,WRK,KFREE,LFREE,IDIST) 479 IF (IDIST .LT. 0) GO TO 1095 480C IF IDIST.LT.0 NO MORE DISTRIBUTIONS 481 ICDSYM= MULD2H(ISMO(IC),ISMO(ID)) 482 ICW = ISW(IC) 483 IDW = ISW(ID) 484 IDI = ID 485 ICI = IC 486 ICIW = ICW 487 IDIW = IDW 488C 489C ORDER INDICES OF C AND D SUCH THAT C =>D 490C 491 IF ( IDW.GT.ICW ) THEN 492 ICIW = IDW 493 IDIW = ICW 494 ICI = ID 495 IDI = IC 496 CALL DGETRN(WRK(KH2),NORBT,NORBT) 497 ENDIF 498 NCIW = ICIW - NISHT 499 NDIW = IDIW - NISHT 500C 501 IF ( IPRRSP.GT.150 ) THEN 502 WRITE(LUPRI,'(/A,2I8)')' DIRAC DISTRIBUTION : IC,ID ',IC,ID 503 CALL OUTPUT(WRK(KH2),1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 504 write(lupri,*) 'NOSIM=', NOSIM 505 ENDIF 506 IF (NOSIM.GT.0) THEN 507 IF (.NOT.DIRFCK) CALL FONEDR(NOSIM,ICI,IDI,WRK(KH2), 508 & WRK(KFCOCO),WRK(KFVOCO),FCX,FVX,ZYMAT, 509 & WRK(KDENA),WRK(KDENB),WRK(KFREE),LFREE) 510C CALL FONEDR(NSIM,ICI1,IDI1,H2D,FCOCO,FVOCO, 511C & FCOEX,FVOEX,ZYMAT,DENA,DENB,WRK,LWRK) 512C 513C CONTRIBUTIONS TO QAX AND QBX 514C 515 IF ( (IOBTYP(ICI).EQ.JTACT) .AND. (IOBTYP(IDI).EQ.JTACT)) 516 & CALL QONEDI(NOSIM,NCIW,NDIW,ICDSYM,QAX,QBX,ZYMAT, 517 & WRK(KH2),WRK(KH2X),PVX,WRK(KPVCD2),WRK(KPVCD3), 518 & WRK(KPVDC2),WRK(KPVDC3),WRK(KFREE),LFREE) 519C CALL QONEDI(NOSIM,NCIW,NDIW,ICDSYM, 520C & QAX,QBX,ZYMAT,H2,H2X, 521C & PVX,PVCD2,PVCD3,PVDC2,PVDC3,WRK,LWRK) 522 END IF 523#ifdef DEBUG_RSPOLI 524 525 IF (NOSIM.GT.0) THEN 526 DO IOSIM = 1,NOSIM 527 WRITE(LUPRI,'(/A)') ' DIRAC DISTRIBUTION CONTRIBUTION added' 528 WRITE(LUPRI,'(/A,I5,A)') 529 * ' QAX FOR',IOSIM,' ORBITAL TRIAL VECTOR' 530 CALL OUTPUT(QAX(1,1,IOSIM), 531 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 532 WRITE(LUPRI,'(/A,I5,A)') 533 * ' QBX FOR',IOSIM,' ORBITAL TRIAL VECTOR' 534 CALL OUTPUT(QBX(1,1,IOSIM), 535 * 1,NORBT,1,NASHT,NORBT,NASHT,-1,LUPRI) 536 END DO 537 END IF 538 539#endif 540 IF (NCSIM.GT.0 .AND. .NOT.SOPPA .AND. .NOT.DIRFCK) THEN 541 CALL FTDDR(NCSIM,ICI,IDI,FVTD,WRK(KDTV),WRK(KH2)) 542C CALL FTDDR(NSIM,ICI1,IDI1,FVTD,DTV,H2D) 543 END IF 544C 545 GO TO 1090 546C 547C ALL DIRAC DISTRIBUTIONS HAVE BEEN READ IN 548C 549 1095 CONTINUE 550 IF (LUINTD .GT. 0) CALL GPCLOSE(LUINTD,'KEEP') 551C 552 1096 CONTINUE 553 IF (DIRFCK) THEN 554 CALL MEMREL('RSPOLI.h2m+d',WRK,KFRSAV,KWRK0,KFREE,LFREE) 555 GO TO 1000 556C we now go to the direct part of rspoli in order to build 557C the Fock operators 558 ENDIF 559C 560C ADD CONTRIBUTION TO FCX AND FVX FROM ONE INDEX TRANSFORMED 561C TOTAL SYMMETRIC FOCK MATRICES 562C 563 IF ( NOSIM .GT. 0 ) THEN 564C 565C SUM UP COULOMB CONTRIBUTIONS WHICH ARE OBTAINED USING THE 566C TERM IS SYMMETRIC 567C 568 IF (IPRRSP.GT.50) THEN 569 DO, IOSIM = 1,NOSIM 570 WRITE(LUPRI,'(/A,I5)') 571 & ' FCOEX contribution for vector',IOSIM 572 CALL OUTPUT(FCX(1,1,IOSIM), 573 & 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 574 WRITE(LUPRI,'(/A,I5)') 575 & ' FCOCO contribution for vector',IOSIM 576 CALL OUTPUT(WRK(KFCOCO+ (IOSIM-1)*NORBT*NORBT), 577 & 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 578 IF (NASHT.GT.0) THEN 579 WRITE(LUPRI,'(/A,I5)') 580 & ' FVOEX contribution for vector',IOSIM 581 CALL OUTPUT(FVX(1,1,IOSIM), 582 & 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 583 WRITE(LUPRI,'(/A,I5)') 584 & ' FVOCO contribution for vector',IOSIM 585 CALL OUTPUT(WRK(KFVOCO+ (IOSIM-1)*NORBT*NORBT), 586 & 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 587 END IF 588 END DO 589 END IF 590 CALL FONEDN(NOSIM,WRK(KFCOCO),WRK(KFVOCO),FCX,FVX) 591 END IF 592C 593 CALL MEMREL('RSPOLI.h2m+d x',WRK,KFRSAV,KWRK0,KFREE,LFREE) 594 GO TO 2000 595 1000 CONTINUE 596C 597C end up here if DIRFCK or HSROHF or DODFT or NASHT.eq.1 598C 599C Allocate work memory: 600C FXCAO,FXVAO,FTCAO,FTVAO must be stored consecutively 601C DXCAO,DXVAO,DTCAO,DTVAO must be stored consecutively 602C (requirement for SIRFCK call) 603C 604 NFOMAT = 0 605 IF (DOFXC) NFOMAT = NFOMAT + NOSIM 606C ... for FXCAO (+ VXxcAO when DFT or srDFT) 607 IF (DOFXV) NFOMAT = NFOMAT + NOSIM 608C ... for FXVAO 609 NFCMAT = 0 610 IF (DOFTC) NFCMAT = NFCMAT + NCSIM 611C ... for FTCAO (with V[2c]xcAO when MC-srDFT) 612 IF (DOFTV) NFCMAT = NFCMAT + NCSIM 613C ... for FTVAO 614 NFMAT = NFOMAT + NFCMAT 615 616 IF (NFMAT .EQ. 0) GOTO 2000 ! may happen for SOPPA 617C 618 NDMAT = NFMAT 619 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 620 NDMAT = NDMAT + 2 621C ... for DtotAO, DXtotAO for SRDFT:SRDFTLTR 622 IF (NFOMAT .GT. 0 .AND. NFCMAT .gt. 0) THEN 623 CALL QUIT( 624 & 'NCSIM .gt. 0 .and. NOSIM .gt. 0 not implemented for MCSRDFT') 625 END IF 626 END IF 627C 628 CALL MEMGET2('REAL','FXCAO',KFXCAO,NFMAT*N2BASX,WRK,KFREE,LFREE) 629 CALL MEMGET2('REAL','DXCAO',KDXCAO,NDMAT*N2BASX,WRK,KFREE,LFREE) 630 IF (DOFXC) THEN 631 KFXVAO = KFXCAO + NOSIM*N2BASX 632 KDXVAO = KDXCAO + NOSIM*N2BASX 633 ELSE 634 KFXVAO = KFXCAO 635 KDXVAO = KDXCAO 636 END IF 637 IF (DOFXV) THEN 638 KFTCAO = KFXVAO + NOSIM*N2BASX 639 KDTCAO = KDXVAO + NOSIM*N2BASX 640 ELSE 641 KFTCAO = KFXVAO 642 KDTCAO = KDXVAO 643 END IF 644 IF (DOFTC) THEN 645 KFTVAO = KFTCAO + NCSIM*N2BASX 646 KDTVAO = KDTCAO + NCSIM*N2BASX 647 ELSE 648 KFTVAO = KFTCAO 649 KDTVAO = KDTCAO 650 END IF 651 IF (DOFTV) THEN 652 KDTOTAO = KDTVAO + NCSIM*N2BASX 653 ELSE 654 KDTOTAO = KDTVAO 655 END IF 656 657#ifdef DEBUG_RSPOLI 658 write(lupri,*) 'rspoli: ncsim,nosim',ncsim,nosim 659 write(lupri,*) 'rspoli: nfmat,ndmat',nfmat,ndmat 660 write(lupri,*) 'rspoli: dohfsrdft,domcsrdft,soppa', 661 & dohfsrdft,domcsrdft,soppa 662 write(lupri,*) 'rspoli: dofxc,dofxv,doftc,doftv', 663 & dofxc,dofxv,doftc,doftv 664 write(lupri,*) 665 & 'rspoli: KDXCAO,KDXVAO,KDTCAO,KDTVAO,KDTOTAO,kfree', 666 & KDXCAO,KDXVAO,KDTCAO,KDTVAO,KDTOTAO,KFREE 667#endif 668 669#ifdef MOD_SRDFT 670 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 671 CALL SRDFT_DENS(WRK(KDTOTAO),UDV,PVX,CMO,WRK,KFREE,LFREE) 672C Get DtotAO, allocation for DXtotAO in SIRFCK.SRDFT can here be 673C used temporarily for DVAO if NASHT.gt.0 674C 675 END IF 676#endif 677C 678C 679C Manu about MC-srDFT: Note that for ncsim >0 and doftc 680C the core transition DM (DTCAO) is set to zero 681C --> the corresponding IFCTYP will therefore be zero. 682C 683C 684C Transform any transition DMs to the AO basis 685 IF (NCSIM.GT.0 .AND. .NOT.SOPPA) THEN 686 IF (DOFTC) CALL DZERO(WRK(KDTCAO),NCSIM*N2BASX) 687 JAO = KDTVAO 688 JMO = KDTV 689 DO ICSIM = 1, NCSIM 690 CALL FCKDEN2(.FALSE.,.TRUE.,DUMMY,WRK(JAO),CMO,WRK(JMO), 691 & WRK(KFREE),LFREE) 692 JAO = JAO + N2BASX 693 JMO = JMO + N2ASHX 694 END DO 695 696 IF (IPRRSP .GT. 50) THEN 697 IF (DOFTC) THEN 698 WRITE(LUPRI,*) 'RSPOLI: first DTCAO matrix of',NCSIM 699 CALL OUTPUT(WRK(KDTCAO),1,NBAST,1,NBAST, 700 & NBAST,NBAST,-1,LUPRI) 701 END IF 702 IF (DOFTV) THEN 703 WRITE(LUPRI,*) 'RSPOLI: first DTVAO matrix of',NCSIM 704 CALL OUTPUT(WRK(KDTVAO),1,NBAST,1,NBAST, 705 & NBAST,NBAST,-1,LUPRI) 706 END IF 707 END IF 708 END IF 709 710 IF (IPRRSP .GT. 50) THEN 711 IF (DOHFSRDFT .OR. DOMCSRDFT) THEN 712 WRITE(LUPRI,'(/A)') 713 & ' RSPOLI SRDFT: total density matrix AO basis DTOTAO' 714 CALL OUTPUT(WRK(KDTOTAO),1,NBAST,1,NBAST, 715 & NBAST,NBAST,-1,LUPRI) 716 END IF 717 END IF 718C 719C Construct NOSIM Inactive and Active density matrices according to eq.27 720C in Chem Phys 119, 297 (1988). 721C 722C They should be consecutive in core. We transpose the inactive 723C and active D-matrices and multiply the inactive D-matrix by 2.0 724C to get things to work for RPA !! 725C 726 IF (DOFXC .AND. NISHT .EQ. 0) THEN ! FXC needed for DFT/srDFT, also when NISHT.eq.0 727 CALL DZERO(WRK(KDXCAO),NOSIM*N2BASX) 728 END IF 729 JOFFAO = 0 730 DO IOSIM = 1,NOSIM 731 CALL DEQ27(CMO,ZYMAT(1,IOSIM),UDV,WRK(KDXCAO+JOFFAO), 732 & WRK(KDXVAO+JOFFAO),WRK(KFREE),LFREE) 733 JOFFAO = JOFFAO + N2BASX 734 END DO 735C 736C Set ISYMDM and IFCTYP information 737C 738 739 allocate(isymdm(nfmat)) 740 allocate(ifctyp(nfmat)) 741 742 IF (HSROHF) THEN 743 CALL DAXPY(NOSIM*N2BASX,D1,WRK(KDXVAO),1,WRK(KDXCAO),1) 744 CALL DSCAL(NOSIM*N2BASX,-D1,WRK(KDXVAO),1) 745 END IF 746 747 CALL DZERO(WRK(KFXCAO),NFMAT*N2BASX) 748 749 JDXCAO = KDXCAO 750 DO I = 1,NFMAT 751C 752C Changed ISYMDM() from 1 to JWOPSY by hint from hjj, kr-may-95 753C 754 ISYMDM(I) = JWOPSY 755C 756C IFCTYP = XY 757C X indicates symmetry about diagonal 758C X = 0 No symmetry 759C X = 1 Symmetric 760C X = 2 Anti-symmetric 761C Y indicates contributions 762C Y = 0 no contribution ! 763C Y = 1 Coulomb 764C Y = 2 Exchange 765C Y = 3 Coulomb + Exchange 766C 767C Check if density matrix is unsymmetric (IX=0), 768C symmetric (IX=10), antisymmetric (IX=20), or zero matrix (IX=30) 769C to threshold THRZER 770C 771 IX = 10 * MATSYM(NBAST,NBAST,WRK(JDXCAO),THRZER) 772C INTEGER FUNCTION MATSYM(N,NDIM,AMAT,THRZER) 773C 774 IF (IX .EQ. 30) THEN 775C zero density matrix, do nothing ! 776 IFCTYP(I) = 0 777 ELSE IF (IX .EQ. 20) THEN 778C Only exchange if antisymmetric density matrix 779 IFCTYP(I) = IX + 2 780 ELSE IF (TRPLET) THEN 781C only exchange 782 IFCTYP(I) = IX + 2 783 ELSE 784C Coulomb+exchange 785 IFCTYP(I) = IX + 3 786 END IF 787 JDXCAO = JDXCAO + N2BASX 788C 789 IF (IPRRSP .GT. 10) write(lupri,*) 790 & 'Fock matrix I, IX, IFCTYP(I) = ', I, IX ,IFCTYP(I) 791 792 END DO ! I = 1,NFMAT 793C 794C 795 IF (HSROHF) THEN 796C Change IFCTYP for all Fa matrices 797 DO I=1,NOSIM 798 IF (TRPLET) THEN 799 IFCTYP(I)=10*(IFCTYP(I)/10) + 2 800 IFCTYP(NOSIM+I)=10*(IFCTYP(NOSIM+I)/10) + 3 801 ELSE 802 IFCTYP(I)=10*(IFCTYP(I)/10) + 3 803 IFCTYP(NOSIM+I)=10*(IFCTYP(NOSIM+I)/10) + 2 804 END IF 805 END DO 806 END IF 807 808 IF (TRPLET) THEN 809 DO I = 1,NFMAT 810 IFCTYP(I) = -IFCTYP(I) ! tell SIRFCK this is triplet 811 ! (needed for MC-srDFT) 812 END DO 813 END IF 814 815 IF (IPRRSP .GT. 60) THEN 816 WRITE(LUPRI,*) 'RSPOLI:SIRFCK NFMAT DIRFCK',NFMAT,DIRFCK 817 WRITE(LUPRI,*) 'ISYMDM ',(ISYMDM(I),I=1,NFMAT) 818 WRITE(LUPRI,*) 'IFCTYP ',(IFCTYP(I),I=1,NFMAT) 819 DO I=1,NOSIM 820 JDXCAO = KDXCAO + (I-1)*N2BASX 821 WRITE(LUPRI,*) 'RSPOLI:SIRFCK DXCAO no.',I 822 CALL OUTPUT(WRK(JDXCAO),1,NBAST,1,NBAST, 823 & NBAST,NBAST,-1,LUPRI) 824 END DO 825 END IF 826#ifdef DEBUG_RSPOLI 827 write(lupri,*) 'RSPOLI: before calling sirfck' 828 IF (DOFXC.AND.NOSIM.gt.0) THEN 829 write(lupri,*) 'WRK(KFXCAO)=' 830 CALL OUTPUT(WRK(KFXCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 831 END IF 832 IF (DOFXV.and.NOSIM.gt.0) THEN 833 write(lupri,*) 'WRK(KFXVAO)=' 834 CALL OUTPUT(WRK(KFXVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 835 END IF 836 IF (DOFTC.and.NCSIM.gt.0) THEN 837 write(lupri,*) 'WRK(KFTCAO)=' 838 CALL OUTPUT(WRK(KFTCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 839 END IF 840 IF (DOFTV.and.NCSIM.gt.0) THEN 841 write(lupri,*) 'WRK(KFTVAO)=' 842 CALL OUTPUT(WRK(KFTVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 843 END IF 844!!! 845!!! print also density matrices 846!!! 847 IF (DOFXC.AND.NOSIM.gt.0) THEN 848 write(lupri,*) 'WRK(KDXCAO)=' 849 CALL OUTPUT(WRK(KDXCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 850 END IF 851 IF (DOFXV.and.NOSIM.gt.0) THEN 852 write(lupri,*) 'WRK(KDXVAO)=' 853 CALL OUTPUT(WRK(KDXVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 854 END IF 855 IF (DOFTC.and.NCSIM.gt.0) THEN 856 write(lupri,*) 'WRK(KDTCAO)=' 857 CALL OUTPUT(WRK(KDTCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 858 END IF 859 IF (DOFTV.and.NCSIM.gt.0) THEN 860 write(lupri,*) 'WRK(KDTVAO)=' 861 CALL OUTPUT(WRK(KDTVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 862 END IF 863#endif 864 865! Manu test !!!!!!!!!!! 866 DFTADX = DFTADD 867 DFTADD = .FALSE. 868 CALL SIRFCK(WRK(KFXCAO),WRK(KDXCAO),NFMAT,ISYMDM,IFCTYP, 869 & DIRFCK,WRK(KFREE),LFREE) 870 DFTADD = DFTADX 871 872#ifdef DEBUG_RSPOLI 873 write(lupri,*) 'just after call sirfck in rspoli' 874 IF (DOFXC.AND.NOSIM.gt.0) THEN 875 write(lupri,*) 'WRK(KFXCAO)=' 876 CALL OUTPUT(WRK(KFXCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 877 END IF 878 IF (DOFXV.and.NOSIM.gt.0) THEN 879 write(lupri,*) 'WRK(KFXVAO)=' 880 CALL OUTPUT(WRK(KFXVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 881 END IF 882 IF (DOFTC.and.NCSIM.gt.0) THEN 883 write(lupri,*) 'WRK(KFTCAO)=' 884 CALL OUTPUT(WRK(KFTCAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 885 END IF 886 IF (NCSIM.gt.0) THEN 887 write(lupri,*) 'WRK(KFTVAO)=' 888 CALL OUTPUT(WRK(KFTVAO),1,NBAST,1,NBAST,NBAST,NBAST,-1,LUPRI) 889 END IF 890#endif 891 892 IF (HSROHF) THEN 893C 894C Calculated matrices are (FC+FV,-FV(exch)=FV-Q) 895C Input fock matrices (from SIRIFC) are of the form 896C (FC+Q,FV-Q), where Q=FV(coul) + 2*FV(exch) 897C adapt to this form (which works in RSPORB) 898C 899C FC+Q 900 CALL DAXPY(NOSIM*N2BASX,-D1,WRK(KFXVAO),1,WRK(KFXCAO),1) 901 END IF 902C ***** Transform FXCAO and FXVAO to MO basis and add to FCX 903C ***** and FVX, resp. (the one-index transformations of FC 904C ***** and FV are added after this routine). 905C CALL AUTPV(ISYM,JSYM,U,V, 906C & PRPAO,NBAS,NBAST,PRPMO,NORB,NORBT,WRK,LWRK) 907C 908C 909 IF (NOSIM .GT. 0) THEN 910 JOFFAO = 0 911 DO 5900 IOSIM = 1,NOSIM 912 DO 300 ISYM=1,NSYM 913 JSYM = MULD2H(ISYM,KSYMOP) 914 NORBI = NORB(ISYM) 915 NORBJ = NORB(JSYM) 916 IF (NORBI.EQ.0 .OR. NORBJ.EQ.0) GO TO 300 917 IF (DOFXC) THEN 918 CALL AUTPV(ISYM,JSYM,CMO(ICMO(ISYM)+1),CMO(ICMO(JSYM)+1), 919 & WRK(KFXCAO+JOFFAO),NBAS,NBAST,FCX(1,1,IOSIM),NORB, 920 & NORBT,WRK(KFREE),LFREE) 921 END IF 922 IF (DOFXV) THEN 923 CALL AUTPV(ISYM,JSYM,CMO(ICMO(ISYM)+1),CMO(ICMO(JSYM)+1), 924 & WRK(KFXVAO+JOFFAO),NBAS,NBAST,FVX(1,1,IOSIM),NORB, 925 & NORBT,WRK(KFREE),LFREE) 926 END IF 927 300 CONTINUE 928C 929 JOFFAO = JOFFAO + N2BASX 930 5900 CONTINUE 931 END IF 932C 933C Manu+hjaaj Aug 09: 934C transform FTCAO and FTVAO to FTC(in FTV) and FTV 935C (for DOMCSRDFT and maybe soon also DOMC (standard MCSCF)) 936C 937 IF (NCSIM .GT. 0) THEN 938 JOFFAO = 0 939 DO ICSIM = 1,NCSIM 940 DO 320 ISYM=1,NSYM 941 JSYM = MULD2H(ISYM,KSYMOP) 942 NORBI = NORB(ISYM) 943 NORBJ = NORB(JSYM) 944 IF (NORBI.EQ.0 .OR. NORBJ.EQ.0) GO TO 320 945 IF (DOFTC) THEN 946C FCT exists with srDFT contribution for srDFT!!! 947C NB! saved in FVTD(*,*,NCSIM+ICSIM) after FVT 948 CALL AUTPV(ISYM,JSYM,CMO(ICMO(ISYM)+1),CMO(ICMO(JSYM)+1), 949 & WRK(KFTCAO+JOFFAO),NBAS,NBAST,FVTD(1,1,NCSIM+ICSIM), 950 & NORB,NORBT,WRK(KFREE),LFREE) 951 END IF 952 IF (DOFTV) THEN 953 CALL AUTPV(ISYM,JSYM,CMO(ICMO(ISYM)+1),CMO(ICMO(JSYM)+1), 954 & WRK(KFTVAO+JOFFAO),NBAS,NBAST,FVTD(1,1,ICSIM), 955 & NORB,NORBT,WRK(KFREE),LFREE) 956 END IF 957 320 CONTINUE 958C 959 JOFFAO = JOFFAO + N2BASX 960 END DO ! ICSIM = 1,NCSIM 961 END IF ! NCSIM .gt. 0 962C 963 CALL MEMREL('RSPOLI.sirfck',WRK,KFRSAV,KWRK0,KFREE,LFREE) 964C 965 IF (IPRRSP.GT.30) THEN ! Manu debug1 966 DO, ICSIM = 1,NCSIM 967 WRITE(LUPRI,'(/2A,I5)')' FVTD before aotomo' 968 * ,' FOR VECTOR ',ICSIM 969 CALL OUTPUT(FVTD(1,1,ICSIM), 970 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 971 IF (DOMCSRDFT) THEN 972 WRITE(LUPRI,'(/2A,I5)')'V^[2c]sr_Hxc before aotomo' 973 * ,' FOR VECTOR ',ICSIM 974 CALL OUTPUT(FVTD(1,1,NCSIM+ICSIM), 975 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 976 END IF 977 END DO 978 END IF 979C Manu debug2 980C 981C The following would comply for a call from LAGRAN in HF 982C 983 if (IPRRSP .gt. 10) then 984 write(lupri,*) 'LAGRAN? NCSIM,TDHF,TRPLET,HSROHF,DFTADD:', 985 & NCSIM,TDHF,TRPLET,HSROHF,DFTADD 986 end if 987 IF (NCSIM.EQ.1 .AND. TDHF .AND. TRPLET) THEN 988C 989C For high spin - dft we just unpack FC and FV 990C 991 IF (HSROHF.OR.DFTADD) THEN 992 if (IPRRSP .gt. 10) then 993 write(lupri,*) 'HSROHF .or. DFTADD, FV=' 994 CALL OUTPKB(FV,NORB,NSYM,-1,LUPRI) 995 end if 996 CALL MEMGET2('REAL','FC',KFC,N2ORBX,WRK,KFREE,LFREE) 997 CALL MEMGET2('REAL','TMP',KTMP,NNORBX,WRK,KFREE,LFREE) 998 CALL PKSYM1(WRK(KTMP),FV,NORB,NSYM,-1) 999 CALL DSPTSI(NORBT,WRK(KTMP),FVTD) 1000 CALL PKSYM1(WRK(KTMP),FC,NORB,NSYM,-1) 1001 CALL DSPTSI(NORBT,WRK(KTMP),WRK(KFC)) 1002 ELSE 1003 CALL MEMGET2('REAL','DA',KDA,N2BASX,WRK,KFREE,LFREE) 1004 CALL MEMGET2('REAL','FA',KFA,N2BASX,WRK,KFREE,LFREE) 1005 CALL FCKDEN( 1006 & .FALSE.,.TRUE.,DUMMY,WRK(KDA),CMO,WRK(KDTV), 1007 & WRK(KFREE),LFREE 1008 & ) 1009 CALL DZERO(WRK(KFA),N2BASX) 1010 ISYMDM(1)=1 1011 IFCTYP(1)=12 1012 CALL SIRFCK( 1013 & WRK(KFA),WRK(KDA),1,ISYMDM,IFCTYP,DIRFCK,WRK(KFREE),LFREE) 1014 CALL AOTOMO(WRK(KFA),FVTD,CMO,1,WRK(KFREE),LFREE) 1015C CALL AOTOMO(XAO,XMO,CMO,XSYM,WRK,LWRK) 1016 END IF 1017 END IF 1018 1019 deallocate(isymdm) 1020 deallocate(ifctyp) 1021 1022 2000 CONTINUE 1023C 1024C From here the code is as before DIRFCK 1025C 1026 IF (NOSIM .GT. 0) THEN 1027C 1028C 1) add DFT contributions 1029C 1030 IF (DODFT) THEN 1031#ifdef DEBUG_RSPOLI 1032 write(lupri,*) 'DODFT: before call dft_lin_resp*' 1033 IF (DOFXC) THEN 1034 write(lupri,*) 'FCX before dft_lin_resp*=' 1035 CALL OUTPUT(FCX,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1036 END IF 1037 IF (DOFXV) THEN 1038 write(lupri,*) 'FCV before dft_lin_resp*=' 1039 CALL OUTPUT(FCV,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1040 END IF 1041#endif 1042C Special case of alpha density only (old code) 1043 IF((NISHT.EQ.0) .AND. (NASHT.GT.0)) THEN 1044 DO IOSIM = 1, NOSIM 1045 call dft_lin_respab(fcx(1,1,iosim),fvx(1,1,iosim), 1046 & cmo,zymat(1,iosim),trplet,ksymop,wrk(kfree), 1047 & lfree,iprrsp) 1048 END DO 1049 ENDIF 1050 1051 1052 ! Proceed normal way otherwise 1053 IF ((NASHT.GT.0) .AND. (NISHT.GT.0)) THEN 1054 call dft_lin_respab_b(NOSIM,fcx,fvx,cmo,zymat, 1055 & trplet, ksymop,wrk(kfree),lfree,iprdft) 1056 ELSE 1057 call dft_lin_respf(NOSIM,fcx,cmo,zymat, 1058 & trplet, ksymop,wrk(kfree),lfree,iprdft) 1059 ENDIF 1060#ifdef DEBUG_RSPOLI 1061 write(lupri,*) 'DODFT: after call dft_lin_resp*' 1062 IF (DOFXC.AND.NOSIM.gt.0) THEN 1063 write(lupri,*) 'FCX after dft_lin_resp*=' 1064 CALL OUTPUT(FCX,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1065 END IF 1066 IF (DOFXV.and.NOSIM.gt.0) THEN 1067 write(lupri,*) 'FCV after dft_lin_resp*=' 1068 CALL OUTPUT(FCV,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1069 END IF 1070#endif 1071 END IF ! IF (DODFT) THEN 1072C 1073C ADD CONTRIBUTION FROM THE A(0)S(2) TERM IF HRPA 1074C 1075 IF (SOPPA) THEN 1076 CALL A0S2(FCX,FC,UDV,ZYMAT,WRK(KTZYMT),WRK(KH2X),WRK(KH2XP), 1077 * WRK(KCOEUN),NOSIM, WRK,LWRK) 1078C 1079C Cancel out the nonsymmetric parts of A(2) matrix. Use the explicitly 1080C constructed A(2) matrix in XINDX. 1081C 1082 IF (.NOT. A2EXIST) THEN 1083 CALL DSCAL(N2ORBX,DP5,XINDX(KAB2),1) 1084 IF (IPRRSP .GT. 40) THEN 1085 WRITE(LUPRI,'(/A)') 1086 & ' SOPPA A(2)(ai,bj) packed as .5*A(2)(i,j) and .5*A(2)(a,b)' 1087 CALL OUTPUT(XINDX,1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1088 END IF 1089 A2EXIST = .TRUE. 1090 END IF 1091 CALL HRPAHM(FCX,XINDX(KAB2),ZYMAT,NOSIM) 1092 END IF 1093C 1094C ADD CONTRIBUTION TO FCX AND FVX FROM ONE-INDEX TRANSFORMED 1095C TOTAL SYMMETRIC FOCK MATRICES 1096C 1097 CALL FCKOIN(NOSIM,FC,FV,ZYMAT,FCX,FVX) 1098C 1099C 1100 IF ( (KZCONF .GT. 0) .AND. .NOT.SOPPA) THEN 1101C 1102C HALF TRANSFORMED INTEGRALS ARE NOW IN H2XAC AND CONFIGURATION 1103C PART OF LINEAR TRANSFORMATION CAN BE CARRIED OUT FOR EACH ORBITAL 1104C TRIAL VECTOR 1105C 1106 CALL H2XSIG(NOSIM,FCX,ZYMAT,WRK(KH2XAC), 1107 * EVECS(1+NCSIM*KZYVAR), 1108 * XINDX,WRK(KFREE),LFREE) 1109C 1110 END IF 1111C 1112C 1113C DISTRIBUTE FOCK AND Q MATRICES IN LINEAR TRANSFORMED VECTORS 1114C 1115 IF ((IPRRSP.GT.30).AND.(NASHT.GT.1)) THEN 1116 DO 700 IOSIM = 1,NOSIM 1117 WRITE(LUPRI,'(/A,I5,A)') ' QAX FOR ',IOSIM, 1118 * ' ORBITAL TRIAL VECTOR' 1119 CALL OUTPUT(QAX(1,1,IOSIM),1,NORBT,1,NASHT, 1120 * NORBT,NASHT,-1,LUPRI) 1121 WRITE(LUPRI,'(/A,I5,A)') ' QBX FOR',IOSIM, 1122 * ' ORBITAL TRIAL VECTOR' 1123 CALL OUTPUT(QBX(1,1,IOSIM),1,NORBT,1,NASHT, 1124 * NORBT,NASHT,-1,LUPRI) 1125 700 CONTINUE 1126 END IF 1127 IF (IPRRSP.GT.50) THEN 1128 DO 1012 IOSIM = 1,NOSIM 1129 WRITE(LUPRI,'(/A,I5)') 1130 * ' FCX total contribution for vector',IOSIM 1131 CALL OUTPUT(FCX(1,1,IOSIM), 1132 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1133 IF (DOFXV) THEN 1134 WRITE(LUPRI,'(/A,I5)') 1135 * ' FVX total contribution for vector',IOSIM 1136 CALL OUTPUT(FVX(1,1,IOSIM), 1137 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1138 END IF 1139 1012 CONTINUE 1140 END IF 1141#ifdef DEBUG_RSPOLI 1142 write(lupri,*) 'I (RSPOLI) call RSPORB with ONEIND true' 1143#endif 1144 CALL RSPORB(.TRUE.,NOSIM,FC,FCX,FVX,QAX,QBX,UDV, 1145 * EVECS(NCSIM*KZYVAR+1)) 1146C CALL RSPORB(ONEIND,NSIM,FC,FCX,FVX,QAX,QBX,UDV,EVECS) 1147 1148 END IF ! IF (NOSIM.GT.0) THEN 1149C ------------------------------------------------------------------ 1150 IF (NCSIM.GT.0) THEN 1151 IF (DOMCSRDFT) THEN 1152C NOTE, the TDM passed to SIRFCK is really 1153C minus the TDM, thus we need to multiply 1154C the FTC matrices in FVTD(1,1,NCSIM+ICSIM) by minus 1 1155C to get the right contribution from the DFT 1156C module in srDFT. 1157C 1158 CALL DSCAL(NCSIM*N2ORBX,DM1,FVTD(1,1,NCSIM+1),1) 1159 END IF 1160C 1161 IF (IPRRSP.GT.30) THEN 1162 DO ICSIM = 1,NCSIM 1163 WRITE(LUPRI,'(/A,I5)') 1164 * ' FVTD total contribution FOR VECTOR ',ICSIM 1165 CALL OUTPUT(FVTD(1,1,ICSIM), 1166 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1167 IF (.NOT.SOPPA) THEN 1168 WRITE(LUPRI,'(/A,I5,A)') 1169 * ' QATD FOR ',ICSIM,' CONFIGURATION TRIAL VECTOR' 1170 CALL OUTPUT(QATD(1,1,ICSIM),1,NORBT,1,NASHT, 1171 * NORBT,NASHT,-1,LUPRI) 1172 WRITE(LUPRI,'(/A,I5,A)') 1173 * ' QBTD FOR',ICSIM,' CONFIGURATION TRIAL VECTOR' 1174 CALL OUTPUT(QBTD(1,1,ICSIM),1,NORBT,1,NASHT, 1175 * NORBT,NASHT,-1,LUPRI) 1176 END IF 1177 IF (DOMCSRDFT) THEN 1178 WRITE(LUPRI,'(/A,I5)') 1179 * 'V^[2c]sr_Hxc total contribution for conf vector',ICSIM 1180 CALL OUTPUT(FVTD(1,1,NCSIM+ICSIM), 1181 * 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1182 END IF 1183 END DO 1184 END IF 1185 1186 IF (HSROHF.OR.DODFT) THEN 1187 CALL RSPORB(.TRUE.,NCSIM,FC,WRK(KFC),FVTD,QATD,QBTD,UDV, 1188 & EVECS) 1189C 1190C Note: This is to cancel the change of sign in RSPORB 1191C In ordinary ROHF/MCSCF, RSPORB is called with negative 1192C matrices (due to negative active densities produced by 1193C the call to RSPTDM above when the CI vector is the reference state) 1194C 1195C 1196 CALL DSCAL(KZYVAR,DM1,EVECS,1) 1197 ELSE 1198C (FVTD(1,1,NCSIM+1:NCSIM+NCSIM) contains the V[2c]xcsr for srDFT, 1199C it is NOT used when normal MCSCF) 1200 CALL RSPORB(.FALSE.,NCSIM,FC,FVTD(1,1,NCSIM+1),FVTD, 1201 * QATD,QBTD,WRK(KDTV),EVECS) 1202 END IF 1203#ifdef MOD_SRDFT 1204 IF (DOMCSRDFT) THEN 1205C ... special MCSRDFT contribution to csf sigma-vectors, 1206C stored in FTV(1,NCSIM+ICSIM) in SIRTR1 after the 1207C standard FTV matrices. 1208C In paper: V^([2c]xc-SR) matrix. 1209 1210 IF (TRPLET) THEN 1211 JSPIN1 = 1 1212 JSPIN2 = 0 1213 ELSE 1214 JSPIN1 = 0 1215 JSPIN2 = 0 1216 END IF 1217C 1218C Aug. 09: structure inspired from rspsol:SLVSC routine /hjaaj 1219C 1220 CALL MEMGET2('REAL','SRCVEC',KSRCVEC,KZCONF,WRK,KFREE,LFREE) 1221 CALL MEMGET2('REAL','SRW' ,KSRW ,N2ASHX,WRK,KFREE,LFREE) 1222 JEVECS = 1 1223 DO ICSIM = 1,NCSIM 1224 CALL GETAC1(FVTD(1,1,NCSIM+ICSIM),WRK(KSRW)) 1225 IF (IPRRSP .GT. 40) THEN 1226 WRITE(LUPRI,'(/A,I3)') 1227 & 'MCSRDFT "FTC"="V^([2c],xc-SR)" no.',ICSIM 1228 CALL OUTPUT(FVTD(1,1,NCSIM+ICSIM), 1229 & 1,NORBT,1,NORBT,NORBT,NORBT,-1,LUPRI) 1230 WRITE(LUPRI,'(/A,I3)') 'MCSRDFT SRLTRAC no.',ICSIM 1231 CALL OUTPUT(WRK(KSRW), 1232 & 1,NASHT,1,NASHT,NASHT,NASHT,-1,LUPRI) 1233 END IF 1234C 1235#ifdef DEBUG_RSPOLI 1236 IF (IPRRSP.GT.101) THEN 1237 WRITE(LUPRI,*) 1238 * ' Z CONF PART OF LINEAR TRANSFORMED', 1239 * ' before Vsr[2c] contr.', 1240 * 'CONF TRIAL VECTOR' 1241 write(lupri,*) 'kzvar = ', kzvar 1242 write(lupri,*) 'kzconf = ', kzconf 1243 CALL OUTPUT(EVECS(JEVECS),1,KZVAR,1,2, 1244 * KZVAR,2,-1,LUPRI) 1245 END IF 1246#endif 1247C 1248C Z part of linear transformation (operator: V[2c]) 1249C (The Y part configuration contribution is zero) ---- Manu: no, it is 1250C NOT ! 1251C 1252 CALL CISIGD(IREFSY,KSYMST,NCREF,KZCONF,WRK(KCREF), 1253 & EVECS(JEVECS),WRK(KSRW),DUMMY, 1254 & .TRUE.,.FALSE.,XINDX,JSPIN1,JSPIN2, 1255 & WRK(KFREE),LFREE) 1256#ifdef DEBUG_RSPOLI 1257 write(lupri,*) 'Z part ... V[2c]' 1258 write(lupri,*) 'IREFSY = ', IREFSY 1259 write(lupri,*) 'KSYMST = ', KSYMST 1260#endif 1261 1262 IF (IREFSY .EQ. KSYMST .AND. .NOT.TRPLET) THEN 1263 FAC = DDOT(KZCONF,EVECS(JEVECS),1,WRK(KCREF),1) 1264 CALL DAXPY(KZCONF,(-FAC),WRK(KCREF),1,EVECS(JEVECS),1) 1265 END IF 1266C 1267 IF (IPRRSP.GT.101) THEN 1268 WRITE(LUPRI,*) 1269 * ' srDFT Z CONF PART OF LINEAR TRANSFORMED ', 1270 * 'CONF TRIAL VECTOR' 1271 WRITE(LUPRI,*)' REFERENCE COMPONENT PROJECTED OUT', 1272 & ', factor = ',-FAC 1273 CALL OUTPUT(EVECS(JEVECS),1,KZVAR,1,2, 1274 * KZVAR,2,-1,LUPRI) 1275 END IF 1276C 1277! Manu: build the Y part of the sigma vector (due to the Vsr[2c] operator) 1278C 1279C 1280 IF (.NOT. (TDA .OR. CISRPA)) THEN 1281 JYEVECS=JEVECS+KZVAR 1282 CALL CISIGD(IREFSY,KSYMST,NCREF,KZCONF,WRK(KCREF), 1283 & EVECS(JYEVECS),WRK(KSRW),DUMMY, 1284 & .TRUE.,.FALSE.,XINDX,JSPIN1,JSPIN2, 1285 & WRK(KFREE),LFREE) 1286 IF (IREFSY .EQ. KSYMST .AND. .NOT.TRPLET) THEN 1287 FAC = DDOT(KZCONF,EVECS(JYEVECS),1,WRK(KCREF),1) 1288 CALL DAXPY(KZCONF,(-FAC),WRK(KCREF),1, 1289 & EVECS(JYEVECS),1) 1290 END IF 1291! 1292! Manu: For the Y conf. part, the transition density matrix changes sign. 1293! Therefore we have to multiply the Vsr[2c] contribution to 1294! the Y part of the sigma vector by -1 1295! 1296 CALL DSCAL(KZCONF,DM1,EVECS(JYEVECS),1) 1297! 1298 IF (IPRRSP.GT.101) THEN 1299 WRITE(LUPRI,*) 1300 * ' srDFT Z and Y CONF PART OF LINEAR TRANSFORMED ', 1301 * 'CONF TRIAL VECTOR' 1302 WRITE(LUPRI,*)' REFERENCE COMPONENT PROJECTED OUT', 1303 & ', factor = ',-FAC 1304 CALL OUTPUT(EVECS(JEVECS),1,KZVAR,1,2, 1305 * KZVAR,2,-1,LUPRI) 1306 END IF 1307 END IF ! not (TDA or CISRPA) 1308C 1309C 1310C Special MCSRDFT correction, corresponding to the effective operator 1311C V^([2c],xc-SR) 1312C is saved in FVTD(*,*,NCSIM+ICSIM) in RSPOLI. 1313C 1314 CALL RSP_SRDFTSO(UDV,FVTD(1,1,NCSIM+ICSIM),EVECS(JEVECS)) 1315C 1316 JEVECS = JEVECS + KZYVAR 1317 END DO ! ICSIM = 1,NCSIM 1318 END IF ! DOMCSRDFT 1319#endif 1320 END IF ! NCSIM .gt. 0 1321 CALL MEMREL('RSPOLI.bottom',WRK,KFRSAV,KFRSAV,KFREE,LFREE) 1322C 1323C *** Perform supersymmetry averaging 1324C 1325 IF ( RSPSUP .AND. (KZWOPT.GT. 0) ) THEN 1326 CALL RFANTI(NOSIM,EVECS(NCSIM*KZYVAR+1),ZYMAT,WRK,LWRK) 1327 END IF 1328 IF ( RSPSUP .AND. (KSYMOP. EQ. 1) .AND. (KZWOPT .GT. 0)) THEN 1329 CALL RSPAVE(EVECS(KZCONF+1),KZVAR,2*(NCSIM+NOSIM)) 1330 END IF 1331C 1332C END OF RSPOLI 1333C 1334 CALL QEXIT('RSPOLI') 1335 RETURN 1336 END 1337C --- end of rsp/rspoli.F --- 1338