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 /* Deck zfsdrv */ 19 SUBROUTINE ZFSDRV(WORK,LWORK) 20C 21C CALCULATE AVERAGE VALUE OF PROPERTIES 22C 23C 24#include "implicit.h" 25 DIMENSION WORK(LWORK) 26#include "dummy.h" 27C 28 PARAMETER ( D1 = 1.0D0 , D2 = 2.0D0 , D8 = 8.0D0 ) 29 PARAMETER (DP5=0.5D0, D1P5=1.5D0) 30 PARAMETER ( D1INF = 0.99999D0 , CKMXPR = 1.0D-6 ) 31 DIMENSION ISYMDM(2), IFCTYP(2) 32 CHARACTER SPD(7) 33 DATA SPD/'S','P','D','F','G','H','I'/ 34C 35#include "maxorb.h" 36#include "maxash.h" 37#include "inforb.h" 38#include "infrsp.h" 39#include "infind.h" 40#include "wrkrsp.h" 41#include "rspprp.h" 42#include "infave.h" 43#include "infpri.h" 44#include "inftap.h" 45#include "iratdef.h" 46#include "aovec.h" 47#include "codata.h" 48#include "gfac.h" 49#include "priunit.h" 50#include "blocks.h" 51#include "zfs.h" 52#include "infesr.h" 53#include "infinp.h" 54C 55 LOGICAL ANTI, PANTI, NODPTR, NODV, NOPV, NOCONT, TTIME, 56 & RETUR, NOBLK, DEBUG, DIA2SO, ZFS2EL 57 CHARACTER*5 STRING 58 IF (.NOT.ZFSCAL) RETURN 59 60 CALL QENTER('ZFSDRV') 61C 62 STRING=" " 63 IPRINT = IPRESR 64 NODPTR = IPRINT.GT.10 65 NOBLK = .FALSE. 66 KFREE = 1 67 LFREE = LWORK 68 CALL HEADER("Output from ZFSDRV",0) 69 70C 71C ******************************************************* 72C ***** Set up COMMON /BLOCKS/ for PSORT and TWOINT ***** 73C ******************************************************* 74C 75 CALL MEMGET('INTE',KJSTRS,MXSHEL*MXAOVC*2,WORK,KFREE,LFREE) 76 CALL MEMGET('INTE',KNPRIM,MXSHEL*MXAOVC*2,WORK,KFREE,LFREE) 77 CALL MEMGET('INTE',KNCONT,MXSHEL*MXAOVC*2,WORK,KFREE,LFREE) 78 CALL MEMGET('INTE',KIORBS,MXSHEL*MXAOVC ,WORK,KFREE,LFREE) 79 CALL MEMGET('INTE',KJORBS,MXSHEL*MXAOVC ,WORK,KFREE,LFREE) 80 CALL MEMGET('INTE',KKORBS,MXSHEL*MXAOVC ,WORK,KFREE,LFREE) 81C 82 IPRPAO=0 83 CALL PAOVEC(WORK(KJSTRS),WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 84 & WORK(KJORBS),WORK(KKORBS),0,NOBLK,IPRPAO) 85C 86 CALL MEMREL('AVEDIA:PAOVEC',WORK,KJORBS,KJORBS,KFREE,LFREE) 87C 88 IF (HSROHF) THEN 89 NODV = .FALSE. 90 NOPV = .TRUE. 91 NDMAT = 2 92 ELSE 93 NODV = .TRUE. 94 NOPV = NASHT .LE. 1 95 NDMAT = 0 96 END IF 97C 98 IF (.NOT.NOPV) THEN 99C 100C Transform two-electron density matrix to SO basis 101C 102 ANTI = .FALSE. 103 PANTI = .FALSE. 104 DIA2SO = .FALSE. 105 ZFS2EL = .TRUE. 106 KDTAO = KFREE 107 JPRINT = IPRINT 108 CALL GPOPEN(LUPAO,'ZFSPAO','NEW',' ',' ',IDUMMY,.FALSE.) 109 CALL PTRAN(NODPTR,WORK(KFREE),LFREE,JPRINT,ANTI,PANTI,DIA2SO, 110 & ZFS2EL) 111 CALL PSORG(WORK(KFREE),WORK(KFREE),LFREE,WORK(KNCONT),JPRINT, 112 & ANTI,PANTI) 113C 114 ELSE IF (HSROHF) THEN 115 CALL MEMGET('REAL',KDMAT,NNASHX,WORK,KFREE,LFREE) 116 CALL MEMGET('REAL',KDTAO,2*N2BASX,WORK,KFREE,LFREE) 117 KDVAO=KDTAO+N2BASX 118 DO I=1,NASHT 119 II= KDMAT + I*(I+1)/2 - 1 120 WORK(II) = D1 121 END DO 122 CALL DZERO(WORK(KDTAO),N2BASX) 123 CALL GETDMT(WORK(KDTAO),NDMAT,WORK(KFREE),LFREE,.TRUE.,.FALSE., 124 & .TRUE.,IPRINT) 125 ELSE 126 CALL QUIT( 127 & 'Zero-field splitting requires at least two open shells') 128 END IF 129C 130C Call HERMIT to evaluate expectation value. A lot of these variables 131C may be of interest to control through a input routine. 132C 133 ISYMDM(1) = 1 134 ISYMDM(2) = 1 135 IFCTYP(1) = 14 136 IFCTYP(2) = 14 137 ITYPE = 11 138 MAXDIF = 2 139 JATOM = 0 140 NOCONT = .FALSE. 141 TTIME = .FALSE. 142 RETUR = .FALSE. 143 IPRNTA = 0 144 IPRNTB = 0 145 IPRNTC = 0 146 IPRNTD = 0 147 I2TYP = 0 148 CALL TWOINT(WORK(KFREE),LFREE,WORK(KFREE),WORK(KFREE),WORK(KDTAO), 149 & NDMAT,IREPDM,IFCTYP,DUMMY,IDUMMY,IDUMMY,1,ITYPE, 150 & MAXDIF,JATOM,NODV,NOPV,NOCONT,TTIME,IPRINT,IPRNTA, 151 & IPRNTB,IPRNTC,IPRNTD,RETUR,IDUMMY,I2TYP,WORK(KJSTRS), 152 & WORK(KNPRIM),WORK(KNCONT),WORK(KIORBS), 153 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY, 154 & DUMMY,.FALSE.,.false.) 155C 156 IF (.NOT. NOPV) CALL GPCLOSE(LUPAO,'DELETE') 157 CALL MEMREL('ZFSDRV',WORK,1,1,KFREE,LFREE) 158C 159C Analyze and print 160C 161 CALL ZFSANA(WORK,LWORK,IPRINT) 162 CALL QEXIT('ZFSDRV') 163 RETURN 164 END 165C /* Deck zfsana */ 166 SUBROUTINE ZFSANA(WORK,LWORK,IPRINT) 167#include "implicit.h" 168 DIMENSION WORK(LWORK) 169#include "maxorb.h" 170#include "infinp.h" 171 INTEGER S 172 CALL QENTER('ZFSANA') 173C 174C Assuming determinant basis and high spin projection 175C 176 MULT=ISPIN 177 MULT2=MULT*MULT 178 179 180 KFREE=1 181 LFREE=LWORK 182 CALL MEMGET('COMP',KHZFS,MULT2,WORK,KFREE,LFREE) 183 CALL MEMGET('REAL',KRZFS,MULT2,WORK,KFREE,LFREE) 184 CALL MEMGET('REAL',KCG,MULT2,WORK,KFREE,LFREE) 185 CALL MEMGET('REAL',KHEIG,MULT,WORK,KFREE,LFREE) 186 CALL MEMGET('REAL',KRWORK,3*MULT-2,WORK,KFREE,LFREE) 187 CALL ZFSAN1(MULT,WORK(KHZFS),WORK(KRZFS),WORK(KCG),WORK(KHEIG), 188 & WORK(KRWORK),WORK,KFREE,LFREE,IPRINT) 189 CALL MEMREL('ZFSANA',WORK,1,1,KFREE,LFREE) 190 191 CALL QEXIT('ZFSANA') 192 END 193C /* Deck zfsan1 */ 194 SUBROUTINE ZFSAN1(MULT,HZFS,RZFS,CG,HEIG,RWORK,WORK,KFREE,LFREE, 195 & IPRINT) 196#include "implicit.h" 197 DOUBLE COMPLEX HZFS(MULT,MULT) 198 DIMENSION RZFS(MULT,MULT) 199 DIMENSION CG(MULT,MULT) 200 DIMENSION HEIG(MULT) 201 DIMENSION WORK(*), RWORK(*) 202#include "priunit.h" 203#include "zfs.h" 204#include "codata.h" 205#include "gfac.h" 206 DOUBLE COMPLEX U(-2:2) 207 DIMENSION RU(-2:2) 208 DIMENSION ZFSEIG(3) 209 DOUBLE COMPLEX I 210 PARAMETER ( D1 = 1.0D0 , D2 = 2.0D0 , D6=6.0D0, D8 = 8.0D0 ) 211 PARAMETER (DP5=0.5D0, D1P5=1.5D0, D0=0.0D0) 212 LOGICAL TESTZ 213 DATA TESTZ /.FALSE./ 214 215 CALL QENTER('ZFSAN1') 216 217CBS I=(D0,D1) 218 I=(0D0,1D0) 219 S=DBLE(MULT-1)/2 220 221 IF (IPRINT.GT.5) THEN 222 CALL HEADER('2-electron field gradient tensor',0) 223 CALL OUTPUT(ZFS,1,3,1,3,3,3,1,LUPRI) 224 END IF 225C 226C Quintet density needs proper scaling 227C 228 DENFAC=D1/(3*S*S - S*(S+1)) 229C 230C Traceless part 231C 232 TR3=(ZFS(1,1)+ZFS(2,2)+ZFS(3,3))/3 233 DO K=1,3 234 ZFS(K,K)=ZFS(K,K)-TR3 235 END DO 236C 237C All in cm-1 238C 239 ZFSFAC=-DENFAC*XTKAYS*(GFAC/2)**2*ALPHAC**2 240 CALL DSCAL(9,ZFSFAC,ZFS,1) 241 CALL HEADER('Trace-less zero field splitting tensor (cm-1)',0) 242 CALL OUTPUT(ZFS,1,3,1,3,3,3,1,LUPRI) 243C 244C We now have a pure second rank tensor, get the spherical representation 245C 246 U2R=(ZFS(1,1)-ZFS(2,2))/2 247 U2I=(ZFS(1,2)+ZFS(2,1))/2 248 U1R=(ZFS(1,3)+ZFS(3,1))/2 249 U1I=(ZFS(2,3)+ZFS(3,2))/2 250 U( 2) = U2R + I*U2I 251 U(-2) = DCONJG(U(2)) 252 U( 1) =-U1R - I*U1I 253 U(-1) =-DCONJG(U(1)) 254 U( 0) = (2*ZFS(3,3)-ZFS(1,1)-ZFS(2,2))/SQRT(D6) 255 IF (IPRINT.GT.5) THEN 256 CALL HEADER('Spherical zero field splitting tensor (cm-1)',0) 257 CALL COUTPUT(U,1,5,1,1,5,1,1,LUPRI) 258 END IF 259C 260C Tabulate Clebsh Gordan CG_S2S(M,N) = <S2NM-N|SM> 261C 262 263 CALL DZERO(CG,MULT*MULT) 264 DO IM=1,MULT 265 RM = -S-1+IM 266 DO IN=1,MULT 267 RN = -S-1+IN 268 MN=(IM-IN) 269 IF (MN.EQ.2) THEN 270 CG(IM,IN)=SQRT( 271 & (3*(S+RM-1)*(S+RM)*(S-RM+1)*(S-RM+2))/ 272 & ((2*S-1)*2*S*(S+1)*(2*S+3)) 273 & ) 274 ELSE IF (MN.EQ.1) THEN 275 CG(IM,IN)=(1-2*RM) 276 & * SQRT( 277 & DBLE(3*(S-RM+1)*(S+RM)) / 278 & ((2*S-1)*S*(2*S+2)*(2*S+3)) 279 & ) 280 ELSE IF (MN .EQ. 0) THEN 281 CG(IM,IN)=DBLE(3*RM*RM-S*(S+1))/ 282 & SQRT( DBLE((2*S-1)*S*(S+1)*(2*S+3)) ) 283 ELSE IF (MN.EQ.-1) THEN 284 CG(IM,IN)=(2*RM+1) 285 & * SQRT( 286 & DBLE(3*(S-RM)*(S+RM+1))/ 287 & ((2*S-1)*S*(2*S+2)*(2*S+3)) 288 & ) 289 ELSE IF (MN.EQ.-2) THEN 290 CG(IM,IN)=SQRT( 291 & DBLE(3*(S-RM-1)*(S-RM)*(S+RM+1)*(S+RM+2))/ 292 & ((2*S-1)*S*(2*S+2)*(2*S+3)) 293 & ) 294 END IF 295 END DO 296 END DO 297 IF (IPRINT.GT.5) THEN 298 CALL HEADER('Clebsh Gordan coefficients',0) 299 CALL OUTPUT(CG,1,MULT,1,MULT,MULT,MULT,1,LUPRI) 300 END IF 301C 302C Reduced matrix element of rank 2 evaluated as 303C 304C TSS=<SS|T(2,0)|SS>/<S2S0|SS> 305C 306 TSS=SQRT( (2*S-1) * S*(S+1) * (2*S+3)/6 ) 307C 308C Build hamiltonian matrix 309C 310C H(M,N)=<SM|T|SN>U(M-N)* = <S||T||S><S2NM-N|SM>U(M-N)* 311 312 DO IM=1,MULT 313 DO IN=1,MULT 314 HZFS(IM,IN) = TSS*CG(IM,IN)*DCONJG(U(IM-IN)) 315 END DO 316 END DO 317 IF (IPRINT.GT.5) THEN 318 CALL HEADER('ZFS Hamiltonian (spherical basis)',0) 319 CALL COUTPUT(HZFS,1,MULT,1,MULT,MULT,MULT,1,LUPRI) 320 END IF 321C 322C Diagonalize HZFS for energy splittings 323C 324 CALL ZHEEV('N','U',MULT,HZFS,MULT,HEIG,WORK(KFREE),LFREE, 325 & RWORK,INFO) 326 IF (INFO.LT.0) THEN 327 WRITE(LUERR)'ZFSDRV:CSYEV:INFO=',INFO 328 CALL QUIT('ZFSDRV:LAPACK DIAGONALIZATION CSYEV FAILED') 329 END IF 330 CALL HEADER('ZFS energy eigenvalues (cm-1)',0) 331 CALL OUTPUT(HEIG,1,MULT,1,1,MULT,1,1,LUPRI) 332C 333C Conventions for triplet states 334C 335 IF (MULT.EQ.3) THEN 336C 337C Z the level of greatest splitting, X mid level 338C 339 DZ1=ABS(HEIG(1)) 340 DZ2=ABS(HEIG(2)) 341 DZ3=ABS(HEIG(3)) 342 IF (DZ1.GE.DZ3) THEN 343 IZ=1 344 IX=2 345 IY=3 346 ELSE 347 IZ=3 348 IX=2 349 IY=1 350 END IF 351 DX=HEIG(IX) 352 DY=HEIG(IY) 353 DZ=HEIG(IZ) 354C 355C D and E values for triplet reference states 356C 357 D=-D1P5*DZ 358 IF (D.GT.0) THEN 359 E=ABS(DX-DY)/2 360 ELSE 361 E=-ABS(DX-DY)/2 362 END IF 363 CALL HEADER('Zero field splitting paramters',0) 364 WRITE(LUPRI,'(A,F10.6,A,F10.2,A)')'@ZFS parameter D = ', 365 & D, ' cm-1 = ',D*XTHZ*1D-6/XTKAYS,' MHz' 366 WRITE(LUPRI,'(A,F10.6,A,F10.2,A)')'@ZFS parameter E = ', 367 & E, ' cm-1 = ',E*XTHZ*1D-6/XTKAYS,' MHz' 368C 369C Check convention (ref Langhoff...) 370C 371 IF (ABS(D).LT.3*ABS(E) .OR. D*E.LT.0) THEN 372 WRITE(LUPRI,*) 'WARNING:ZFS Principle axis test failed' 373 NWARN=NWARN+1 374 END IF 375C 376C Get principal axes by diagonalizing D 377C 378 CALL DSYEV('V','U',MULT,ZFS,3,ZFSEIG,WORK(KFREE),LFREE, 379 & INFO) 380 IF (INFO.LT.0) THEN 381 WRITE(LUERR)'ZFSDRV:DSYEV:INFO=',INFO 382 CALL QUIT('ZFSDRV:LAPACK DIAGONALIZATION DSYEV FAILED') 383 END IF 384 CALL HEADER('ZFS tensor eigenvalues and cosines',0) 385 CM2MHZ=XTHZ*1D-6/XTKAYS 386 WRITE(LUPRI,'(2A16,A24)')'cm-1','MHz',' direction cosines' 387 DO K=1,3 388 WRITE(LUPRI,'(2F16.6,3F8.4)') ZFSEIG(K),ZFSEIG(K)*CM2MHZ, 389 & (ZFS(J,K),J=1,3) 390 END DO 391 IF (TESTZ) THEN 392C 393C Diagonalize D first -> U, H will be real 394C 395 396 RU(2)=(ZFSEIG(1)-ZFSEIG(2))/2 397 RU(-2)=RU(2) 398 RU(1)=D0 399 RU(-1)=D0 400 RU( 0) = (2*ZFSEIG(3)-ZFSEIG(1)-ZFSEIG(2))/SQRT(D6) 401C 402C Build hamiltonian matrix 403C 404C H(M,N)=<SM|T|SN>U(M-N)* = <S||T||S><S2NM-N|SM>U(M-N)* 405 406 DO IM=1,MULT 407 DO IN=1,MULT 408 RZFS(IM,IN) = TSS*CG(IM,IN)*RU(IM-IN) 409 END DO 410 END DO 411 CALL HEADER('RZFS Hamiltonian (spherical basis)',0) 412 CALL OUTPUT(RZFS,1,MULT,1,MULT,MULT,MULT,1,LUPRI) 413 CALL DSYEV('N','U',MULT,RZFS,MULT,ZFSEIG,WORK(KFREE),LFREE, 414 & INFO) 415 IF (INFO.LT.0) THEN 416 WRITE(LUERR)'ZFSDRV:DSYEV:INFO=',INFO 417 CALL QUIT('ZFSDRV:LAPACK DIAGONALIZATION DSYEV FAILED') 418 END IF 419 CALL HEADER('ZFS energy eigenvalues (cm-1)',0) 420 CALL OUTPUT(ZFSEIG,1,MULT,1,1,MULT,1,1,LUPRI) 421 END IF 422 END IF 423C 424C Clean up and exit 425C 426 CALL QEXIT('ZFSAN1') 427 END 428! --- end of DALTON/rsp/rspzfs.F --- 429