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! 18C 19C /* Deck rspdm */ 20 SUBROUTINE RSPDM(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR, RHO1,RHO2, 21 * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WORK, 22 * KFREE,LFREE) 23C 24C CONSTRUCT ONE (RHO1) AND TWO (RHO2) PARTICLE DENSITY 25C MATRICES TO BE USED IN RESPONSE CALCULATION. 26C 27C OUTPUT: 28C 29C RHO1(KL) = (L/ E(KL) /R) 30C 31C RHO2 : TWO-BODY TRANSITION DENSITY MATRIX FOR STATE 32C (L/ AND /R) PACKED WITH SQUARE DISTRIBUTIONS 33C 34C RHO2(IJKL) = RHO2[NASHT,NASHT,NASHT,NASHT] 35C 36C RHO2(IJKL) = (L/ E(IJ,KL) - DELTA(JK) E(IL) /R ) 37C 38#include "implicit.h" 39C 40 DIMENSION CL(*), CR(*),RHO1(NASHT,*) 41 DIMENSION RHO2(NASHT,NASHT,NASHT,*) 42 DIMENSION XINDX(*),WORK(*) 43C 44 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0 , SMALL = 1.0D-8) 45C 46C Used from common blocks: 47C INFINP : LSYM,FLAG(*) 48C INFORB : NASHT,N2ASHX,NNASHX 49C INFVAR : ?? not used ?? /920920-hjaaj 50C wrkrsp.h : DETERM 51C 52#include "maxorb.h" 53#include "priunit.h" 54#include "infinp.h" 55#include "inforb.h" 56#include "infvar.h" 57#include "infrsp.h" 58#include "wrkrsp.h" 59C 60C 61 INTEGER RHOTYP 62 LOGICAL USEPSM, CSFEXP, TDM, NORHO2 63 CALL QENTER('RSPDM') 64C 65C Treat NASHT .eq. 1 as a special case 66C 67 IF (NASHT .EQ. 0) GO TO 9999 68 IF (NASHT .EQ. 1) THEN 69 IF ( (ILSYM.EQ.IRSYM) .AND. (ABS(CL(1)-D1).LT.SMALL) 70 * .AND. (ABS(CR(1)-D1).LT.SMALL) ) THEN 71 RHO1(1,1) = D1 72C RHO1 = 1.0 both for ISPIN1 = 0 and ISPIN1 = 1 73 ELSE 74 RHO1(1,1) = D0 75 END IF 76 IF (.NOT. NORHO2) 77 * RHO2(1,1,1,1) = D0 78 GO TO 9999 79 ELSE IF (HSROHF) THEN 80 IF ( (ILSYM.EQ.IRSYM) .AND. (ABS(CL(1)-D1).LT.SMALL) 81 * .AND. (ABS(CR(1)-D1).LT.SMALL) ) THEN 82 CALL DUNIT(RHO1,NASHT) 83 ELSE 84 CALL DZERO(RHO1,N2ASHX) 85 END IF 86 IF (.NOT. NORHO2) THEN 87 DO I=1,NASHT 88 DO J=1,NASHT 89 RHO2(I,I,J,J)=D1 90 RHO2(I,J,J,I)=RHO2(I,J,J,I)-D1 91 END DO 92 END DO 93 END IF 94 GO TO 9999 95 END IF 96C 97C Set RHOTYP and USEPSM 98C RHOTYP = 2 for non-symmetrized density matrix 99C (PV(ij,kl) usually .ne. PV(ij,lk)) 100C USEPSM tells densid if it is allowed to use permutation symmetry 101C for Ms = 0. 102C 103 RHOTYP = 2 104 IF (FLAG(66)) THEN 105 USEPSM = .FALSE. 106 ELSE 107 USEPSM = .TRUE. 108 END IF 109C 110C 111 CSFEXP = .NOT.DETERM 112C 113 IF ( IPRRSP.GT.65 ) THEN 114 WRITE(LUPRI,'(/A)')' ***RSPDM BEFORE CALLING DENSID' 115 WRITE(LUPRI,'(/A,/2I6,2I4,4I8)') 116 * ' ILSYM IRSYM ISPIN1/2 NCLDIM NCRDIM KFREE LFREE:', 117 * ILSYM,IRSYM,ISPIN1,ISPIN2,NCLDIM,NCRDIM,KFREE,LFREE 118 END IF 119 IF ( IPRRSP.GT.120 ) THEN 120 WRITE(LUPRI,'(/A)')' *RSPDM* LEFT REFERENCE VECTOR' 121 CALL OUTPUT(CL,1,NCLDIM,1,1,NCLDIM,1,1,LUPRI) 122 WRITE(LUPRI,'(/A)')' *RSPDM* RIGHT REFERENCE VECTOR' 123 CALL OUTPUT(CR,1,NCRDIM,1,1,NCRDIM,1,1,LUPRI) 124 END IF 125 IF (.NOT. NORHO2) 126 *CALL SETVEC(RHO2,D0,N2ASHX*N2ASHX) 127 CALL SETVEC(RHO1,D0,N2ASHX) 128 CALL DENSID(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR,RHO1,RHO2, 129 * RHOTYP,CSFEXP,USEPSM,NORHO2,ISPIN1,ISPIN2,TDM, 130 * XINDX,WORK,KFREE,LFREE) 131C CALL DENSID(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR,RHO1,RHO2, 132C * RHOTYP,CSFEXP,USEPSM,NORHO2,ISPIN1,ISPIN2,TDM, 133C * XNDXCI,WORK,KFREE,LFREE) 134C 135 if(ci_program .eq. 'LUCITA ')then 136 KRHO1 = KFREE 137 call dzero(work(krho1),n2ashx) 138 CALL DSPTSI(NASHT,RHO1,work(krho1)) 139 call dcopy(n2ashx,work(krho1),1,rho1,1) 140C ... unpack RHO1 using CALL DSPTSI(N,ASP,ASI) 141 end if 142 143 IF ( IPRRSP.GT.90 ) THEN 144 WRITE(LUPRI,'(/A)')' *** LEAVING RSPDM' 145 WRITE(LUPRI,'(/A,/2I6,2I4,4I8)') 146 * ' ILSYM IRSYM ISPIN1/2 NCLDIM NCRDIM KFREE LFREE:', 147 * ILSYM,IRSYM,ISPIN1,ISPIN2,NCLDIM,NCRDIM,KFREE,LFREE 148 END IF 149 IF ( IPRRSP.GT.90 ) THEN 150 WRITE (LUPRI,1100) 151 CALL OUTPUT(RHO1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 152 ENDIF 153 1100 FORMAT(/' RSPDM: One-el. density matrix, active part, MO-basis') 154C 155 IF ( (.NOT. NORHO2 ) .AND. IPRRSP .GE. 120 ) THEN 156 WRITE(LUPRI,'(/A)') 157 * ' RSPDM: Two-body density matrix RHO2 (non-zero elements):' 158 DO 240 L = 1, NASHT 159 DO 240 K = 1, NASHT 160 WRITE(LUPRI,'(/A,2I5)') ' RHO2 distribution K,L:',K,L 161 CALL OUTPUT(RHO2(1,1,K,L),1,NASHT,1,NASHT, 162 * NASHT,NASHT,1,LUPRI) 163 240 CONTINUE 164 END IF 165C 166 9999 CONTINUE 167 CALL QEXIT('RSPDM') 168C ... end of RSPDM. 169 END 170C /* Deck rsptdm */ 171 SUBROUTINE RSPTDM(NCSIM,ILRESY,IRSYM,NCLREF,NCRDIM,CLREF, 172 * CR, RHO1,RHO2, ISPIN1,ISPIN2,TDM,NORHO2, 173 * XINDX,WORK,KFREE,LFREE) 174C 175C CONSTRUCT ONE (RHO1) AND TWO (RHO2) PARTICLE TRANSITION DENSITY 176C MATRICES TO BE USED IN RESPONSE CALCULATION. 177C 178C OUTPUT: 179C 180C RHO1(KL) = (CLREF/ E(KL) /-CR) 181C 182C RHO2 : TWO-BODY TRANSITION DENSITY MATRIX FOR STATE 183C (CLREF/ AND /CR) PACKED WITH SQUARE DISTRIBUTIONS 184C 185C RHO2(IJKL) = RHO2[NASHT,NASHT,NASHT,NASHT] 186C 187C RHO2(IJKL) = (CLREF/ E(IJ,KL) - DELTA(JK) E(IL) /-CR ) 188C 189#include "implicit.h" 190C 191 DIMENSION CLREF(*), CR(KZCONF,*),RHO1(NASHDI,NASHDI,*) 192 DIMENSION RHO2(NASHDI*NASHDI,NASHDI*NASHDI,*) 193 DIMENSION XINDX(*),WORK(*) 194C 195#include "priunit.h" 196#include "wrkrsp.h" 197#include "infdim.h" 198#include "infpri.h" 199#include "infrsp.h" 200#include "inforb.h" 201C 202 PARAMETER ( DM1 = -1.0D0 ) 203 LOGICAL TDM, NORHO2 204C 205 CALL QENTER('RSPTDM') 206C IF (SOPPA) THEN 207C WRITE (LUPRI,*) 'RSPTDM FATAL ERROR: called in SOPPA calc.' 208C CALL QUIT('RSPTDM FATAL ERROR: called in SOPPA calc.') 209C END IF 210 DO 100 ICSIM = 1,NCSIM 211 IF ( IPRRSP.GT.65 ) THEN 212 WRITE(LUPRI,*) '*** RSPTDM: calling RSPDM for ICSIM=',ICSIM 213 END IF 214 CALL RSPDM(ILRESY,IRSYM,NCLREF,NCRDIM,CLREF,CR(1,ICSIM), 215 * RHO1(1,1,ICSIM),RHO2(1,1,ICSIM), 216 * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WORK, 217 * KFREE,LFREE) 218C CALL RSPDM(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR, RHO1,RHO2, 219C * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WORK, 220C * KFREE,LFREE) 221C 222C TAKE CARE OF MINUS SIGN IN CR 223C 224 IF( .NOT.NORHO2 ) 225 * CALL DSCAL(N2ASHX*N2ASHX,DM1,RHO2(1,1,ICSIM),1) 226 CALL DSCAL(N2ASHX, DM1,RHO1(1,1,ICSIM), 1) 227 100 CONTINUE 228C 229 CALL QEXIT('RSPTDM') 230 RETURN 231C ... end of RSPTDM. 232 END 233C /* Deck rspgdm */ 234 SUBROUTINE RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 235 * CL,CR,OVLAP,RHO1,RHO2,ISPIN1,ISPIN2,TDM,NORHO2, 236 * XINDX,WRK,KFREE,LFREE,LREF) 237C 238C CONSTRUCT ONE (RHO1) AND TWO (RHO2) PARTICLE TRANSITION DENSITY 239C MATRICES TO BE USED IN RESPONSE CALCULATION. 240C 241C Adapted from RSPTDM 242C 243C Instead of RSPTDM, this routine can be called with a whole 244C vector, whereas RSPTDM needs to be called with the configura- 245C tion part only. Any state can be put left and right, the 246C logical LREFR determining whether we have the reference state 247C on the right hand side, LREFL determines whether we have it on the 248C right hand side. 249C 250C OUTPUT: <0(L) /../0R> + <0L/../0(L)> 251C L,R between brackets may be reference state 252C 253C RHO1(KL) = (CL/ E(KL) /-CR) + (CR / E(KL) / -CL) 254C 255C RHO2 : TWO-BODY TRANSITION DENSITY MATRIX FOR STATE 256C (CL/ AND /CR) PACKED WITH SQUARE DISTRIBUTIONS 257C 258C RHO2(IJKL) = RHO2[NASHT,NASHT,NASHT,NASHT] 259C 260C RHO2(IJKL) = (CL/ E(IJ,KL) - DELTA(JK) E(IL) /-CR ) 261C 262#include "implicit.h" 263C 264 DIMENSION CL(KZVARL,*), CR(KZVARR,*) 265 DIMENSION RHO1(NASHDI,NASHDI,*) 266 DIMENSION RHO2(NASHDI*NASHDI,NASHDI*NASHDI,*) 267 DIMENSION XINDX(*),WRK(*) 268C 269#include "maxorb.h" 270#include "priunit.h" 271#include "infvar.h" 272#include "inforb.h" 273#include "infrsp.h" 274#include "wrkrsp.h" 275#include "infdim.h" 276#include "qrinf.h" 277C 278 PARAMETER ( D0 = 0.0D0, D1 = 1.0D0, DM1 = -1.0D0 ) 279 LOGICAL TDM, NORHO2, LREF 280C 281 IF (SOPPA) THEN 282 WRITE (LUPRI,*) 'RSPGDM FATAL ERROR: called in SOPPA calc.' 283 CALL QUIT('RSPGDM FATAL ERROR: called in SOPPA calc.') 284 END IF 285 IF ( IPRRSP .GT. 200) THEN 286 WRITE(LUPRI,'(//A)') ' Left hand side vector in RSPGDM' 287 IF ( LREF ) WRITE(LUPRI,'(A)') ' (Reference state)' 288 CALL OUTPUT(CL,1,KZVARL,1,NSIM,KZVARR,NSIM,1,LUPRI) 289 WRITE(LUPRI,'(//A)') ' Right hand side vector in RSPGDM' 290 CALL OUTPUT(CR,1,KZVARR,1,NSIM,KZVARR,NSIM,1,LUPRI) 291 ENDIF 292C 293 IF (.NOT.TDM ) THEN 294 WRITE(LUPRI,'(A)') ' FATAL ERROR: Illegal call of RSPGDM;' 295 WRITE(LUPRI,'(A)') ' TDM is FALSE, must be TRUE' 296 CALL QUIT(' Illegal call of RSPGDM; TDM = FALSE') 297 END IF 298C 299 OVLAP = D0 300C 301C Treat the case NASHT <= 1 302C 303 IF ( NASHT .EQ. 0 ) RETURN 304C 305 IF (NCL .EQ. 0 .OR. NCR .EQ. 0 .OR. NASHT .EQ. 1) THEN 306 CALL DZERO(RHO1,NSIM*N2ASHX) 307 IF (.NOT. NORHO2) CALL DZERO(RHO2,NSIM*N2ASHX*N2ASHX) 308 RETURN 309 END IF 310C 311C Do <0(L) .. 0R> 312C 313 IF( LREF ) THEN 314 IOFF = 1 315 NDREF = MZCONF(1) 316C ... may be NCDETS instead of NCONF inf triplet 317C therefore not NCREF which always is NCONF 318 IF (KZVARL .NE. NDREF ) THEN 319 CALL QUIT(' Illegal call of RSPGDM: KZVARL ne NDREF') 320 ENDIF 321 IF (NCL .NE. NDREF ) THEN 322 CALL QUIT(' Illegal call of RSPGDM: NCL ne NDREF') 323 ENDIF 324 ELSE 325 IOFF = MZVAR(MULD2H(IREFSY,ILSYM)) + 1 326 ENDIF 327C 328 DO 100 ISIM = 1,NSIM 329 IF ( IPRRSP.GT.65 ) THEN 330 WRITE(LUPRI,*) '*** RSPGDM: 1. call of RSPDM for ISIM=',ISIM 331 END IF 332 CALL RSPDM(ILSYM,IRSYM,NCL,NCR,CL(IOFF,ISIM),CR(1,ISIM), 333 * RHO1(1,1,ISIM),RHO2(1,1,ISIM), 334 * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WRK, 335 * KFREE,LFREE) 336C CALL RSPDM(ILSYM,IRSYM,NCLDIM,NCRDIM,CL,CR, RHO1,RHO2, 337C * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WRK, 338C * KFREE,LFREE) 339C 340 FACT = DM1 341 CALL DSCAL(N2ASHX,FACT,RHO1(1,1,ISIM),1) 342 IF (.NOT.NORHO2) THEN 343 CALL DSCAL(N2ASHX*N2ASHX,FACT,RHO2(1,1,ISIM),1) 344 ENDIF 345 IF (ILSYM .EQ. IRSYM) THEN 346 OVLAP = OVLAP - DDOT(NCL,CL(IOFF,ISIM),1,CR(1,ISIM),1) 347 END IF 348C 349C Do <0L .. 0(R)> 350C 351 IOFF = MZVAR(MULD2H(IREFSY,IRSYM)) + 1 352 KRHO1 = KFREE 353 KRHO2 = KRHO1 + N2ASHX 354 IF (NORHO2) THEN 355 KFREE1 = KRHO2 356 ELSE 357 KFREE1 = KRHO2 + N2ASHX*N2ASHX 358 END IF 359 LFREE1 = LFREE - KFREE1 360 IF (LFREE1.LT.0) CALL ERRWRK('RSPGDM',KFREE1-1,LFREE) 361C 362 IF ( IPRRSP.GT.65 ) THEN 363 WRITE(LUPRI,*) '*** RSPGDM: 2. call of RSPDM for ISIM=',ISIM 364 END IF 365 CALL RSPDM(IRSYM,ILSYM,NCR,NCL,CR(IOFF,ISIM),CL(1,ISIM), 366 * WRK(KRHO1),WRK(KRHO2), 367 * ISPIN1,ISPIN2,TDM,NORHO2,XINDX,WRK, 368 * KFREE1,LFREE1) 369C 370C Take care of minus sign in CR 371C ( in case of no LREF both terms have a CR, scale the result) 372C 373 IF ( LREF ) THEN 374 FACT = D1 375 ELSE 376 FACT = DM1 377 END IF 378 CALL DAXPY(N2ASHX,FACT,WRK(KRHO1),1,RHO1(1,1,ISIM),1) 379 IF (.NOT.NORHO2) THEN 380 CALL DAXPY(N2ASHX*N2ASHX,FACT,WRK(KRHO2),1,RHO2(1,1,ISIM),1) 381 ENDIF 382 IF (ILSYM .EQ. IRSYM) THEN 383 OVLAP = OVLAP + FACT*DDOT(NCR,CR(IOFF,ISIM),1,CL(1,ISIM),1) 384 END IF 385 100 CONTINUE 386C 387 IF (IPRRSP .GT. 60) THEN 388 WRITE (LUPRI,*) 'RSPGDM: Overlap factor =',OVLAP 389 END IF 390C 391 RETURN 392C ... end of RSPGDM. 393 END 394C /* Deck pvxdis */ 395 SUBROUTINE PVXDIS(K,L,PVDEN,PVX,IPVDIS) 396C 397C GET 2-ELECTRON DENSITY DISTRIBUTIONS OF VARIOUS TYPE FROM PVX 398C 399C IPVDIS = 1 [**,KL] + { 1 - DELTA(K,L) } [**,LK] 400C IN PVDEN[*,*] 401C SQARE PACKED DISTRIBUTIONS IN PVX 402C 403C IPVDIS = 2 [*K,*L] IN PVDEN[*,*] 404C SQARE PACKED DISTRIBUTIONS IN PVX 405C 406C IPVDIS = 3 [*K,L*] IN PVDEN[*,*] 407C SQARE PACKED DISTRIBUTIONS IN PVX 408C 409C IPVDIS = 4 [**,KL] IN PVDEN[*,*] 410C SQARE PACKED DISTRIBUTIONS IN PVX 411C 412#include "implicit.h" 413C 414 DIMENSION PVDEN(NASHDI,*),PVX(NASHDI,NASHDI,NASHDI,*) 415C 416#include "priunit.h" 417#include "infdim.h" 418#include "inforb.h" 419#include "infrsp.h" 420C 421 PARAMETER ( D1 = 1.0D0 ) 422 IF ( IPRRSP.GT.2000 ) THEN 423 WRITE(LUPRI,'(/A)')' ********* PVXDIS **********' 424 DO 2000 I=1,NASHT 425 DO 2100 J=1,NASHT 426 WRITE(LUPRI,'(/A,2I6)') 427 * ' TWO-BODY DENSITY DISTRIBUTION: I,J ',I,J 428 CALL OUTPUT(PVX(1,1,J,I),1,NASHT,1,NASHT, 429 * NASHT,NASHT,1,LUPRI) 430 2100 CONTINUE 431 2000 CONTINUE 432 END IF 433C 434 GO TO (1,2,3,4),IPVDIS 435 WRITE(LUERR,'(/A,I5)') 436 * ' PVXDIS :INCORRECT SPECIFICATION OF IPVDIS ,IPVDIS:',IPVDIS 437 CALL QUIT('PVXDIS INCORRECT SPECIFICATION OF IPVDIS') 438 1 CONTINUE 439 CALL DCOPY(N2ASHX,PVX(1,1,K,L),1,PVDEN(1,1),1) 440 IF (K.NE.L) CALL DAXPY(N2ASHX,D1,PVX(1,1,L,K),1,PVDEN(1,1),1) 441 GO TO 100 442 2 CONTINUE 443 DO 200 I=1,NASHT 444 CALL DCOPY(NASHT,PVX(1,K,I,L),1,PVDEN(1,I),1) 445 200 CONTINUE 446 GO TO 100 447 3 CONTINUE 448 DO 300 I=1,NASHT 449 CALL DCOPY(NASHT,PVX(1,K,L,I),1,PVDEN(1,I),1) 450 300 CONTINUE 451 GO TO 100 452 4 CONTINUE 453 CALL DCOPY(N2ASHX,PVX(1,1,K,L),1,PVDEN(1,1),1) 454 GO TO 100 455 100 CONTINUE 456 IF (IPRRSP.GT.150) THEN 457 WRITE(LUPRI,'(/A,I5,A,2I5)') ' PVXDIS DISTRIBUTION TYPE', 458 * IPVDIS,' DISTRIBUTION: K,L',K,L 459 CALL OUTPUT(PVDEN,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 460 END IF 461 RETURN 462 END 463