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#ifdef VAR_DEBUG 19#define HSODEBUG .TRUE. 20#else 21#define HSODEBUG .FALSE. 22#endif 23C /* Deck hsoctl */ 24 SUBROUTINE HSOCTL (WORD,GP,CMO,UDV,PV,XINDX,ANTSYM,WRK,LWRK) 25C 26C Olav Vahtras 27C Jun 20, 1990 28C 29C Driver routine for the construction of spin-orbit property vector 30C 31C Input: Spin-orbit component (as read from LUPROP), WORD 32C MO coefficients,CMO 33C First order reduced density matrix, UDV 34C Second order reduced density matrix, PV 35C Index vector, XINDX 36C 37C Output: Spin-orbit property vector returned in first elements of WRK 38C 39#include "implicit.h" 40 DIMENSION GP(*),CMO(*),UDV(*),PV(*),XINDX(*),WRK(LWRK) 41 CHARACTER*8 WORD 42C 43 PARAMETER (D0=0.D0, D1=1.D0, DM1=-1.D0, D2=2.D0) 44 LOGICAL SO2TRA, FILE_EXISTS, FOPEN, DO1,DO2, LORB, LCON, TRPDEN 45 CHARACTER*8 LABEL 46 LOGICAL FNDLAB 47 EXTERNAL FNDLAB 48C 49C Used from common blocks: 50C MULD2H,ISMO,ISW,IOBTYP,IDBTYP,NISHT,NORBT,NASHT,NNASHX 51C INFORB: MULD2H,NISHT,NASHT,NORBT,NNASHX 52C INFIND: ISW,IOBTYP 53C TRHSO : ILXYZ, KSYMSO, OLDTRA 54C INFHSO: IPRHSO, TESTZY, DOSO1, DOSO2 55C INFHYP: HYPCAL 56C INFSMO: SOMOM 57C INFPP : EXMOM 58C INFRSP: IPRRSP,SOPPA,??? 59C 60C-- common blocks: 61#include "maxash.h" 62#include "maxorb.h" 63#include "aovec.h" 64#include "priunit.h" 65#include "dummy.h" 66#include "inforb.h" 67#include "inftap.h" 68#include "infind.h" 69#include "infvar.h" 70#include "infrsp.h" 71#include "wrkrsp.h" 72#include "infpri.h" 73#include "trhso.h" 74#include "rspprp.h" 75#include "infhso.h" 76#include "infhyp.h" 77#include "infsmo.h" 78#include "infpp.h" 79#include "iratdef.h" 80#include "codata.h" 81#include "inflr.h" 82#include "qrinf.h" 83#include "cbihr2.h" 84#include "infesr.h" 85#include "dftcom.h" 86#include "infrank.h" 87C 88 CALL QENTER('HSOCTL') 89C 90 IF (SOPPA) CALL QUIT('HSOCTL: SOPPA not implemented yet!') 91 IPRHSO = MAX(IPRHSO,IPRRSP,IPRPP,IPRLR,IPRESR) 92 KSYMSO=KSYMOP 93 CALL DZERO(GP,KZYVAR) 94 HSOFAC=ALPHAC**2/4 95 IF (WORD(2:2).EQ.'1') THEN 96 OPRANK(INDPRP(WORD)) =1 97 CALL QRGP(WORD,GP,CMO,XINDX,ANTSYM,WRK,LWRK) 98 CALL DSCAL(KZYVAR,HSOFAC,GP,1) 99 GO TO 9999 100 END IF 101 ANTSYM = 1.0D0 102 IF (WORD(2:2).EQ.'2') THEN 103 DO1 = .FALSE. 104 DO2 = .TRUE. 105 END IF 106 IF (WORD(2:2).EQ.' ') THEN 107 DO1 = DOSO1 108 DO2 = DOSO2 109 END IF 110C 111 IF (IPRHSO.GT.2) THEN 112 CALL TIMER('START ',HSOSTA,HSOTIM) 113 CALL HEADER('Output from HSOCTL',-1) 114 WRITE(LUPRI,'(/2A,3X,A)') 115 * ' Spin-orbit property vector calculation', 116 * ' component = ',WORD 117 WRITE(LUPRI,'(/A,I5)')' Print level in HSOCTL: ',IPRHSO 118 IF (TESTZY) WRITE(LUPRI,'(/A)') 119 *' Z and Y parts of configurational property vector explicitly' 120 IF (.NOT.DO1) WRITE(LUPRI,'(/A)') 121 *' Skip one-electron spin-orbit contributions' 122 IF (.NOT.DO2) WRITE(LUPRI,'(/A)') 123 *' Skip two-electron spin-orbit contributions' 124 END IF 125 LORB = KZWOPT.GT.0 126 IF (KSYMOP.EQ.1) THEN 127 LCON = KZCONF.GT.1 128 ELSE 129 LCON = KZCONF.GT.0 130 END IF 131C 132C Check if gradient is on file 133C 134 INQUIRE(FILE='HSOGRAD',EXIST=FILE_EXISTS) 135 LUHSO = -1 136 CALL GPOPEN(LUHSO,'HSOGRAD','UNKNOWN',' ','UNFORMATTED',IDUMMY, 137 & .FALSE.) 138 REWIND LUHSO 139 IF (FILE_EXISTS) THEN 140C ... aug07-hjaaj: gfortran doesn't like BACKSPACE(LUHSO) in FNDLAB on empty file! 141 IF (FNDLAB(WORD,LUHSO)) THEN 142 CALL READT(LUHSO,KZYVAR,GP) 143 CALL GPCLOSE(LUHSO,'KEEP') 144 GO TO 9999 145 END IF 146 END IF 147C 148C ALLOCATE WORK SPACE 149C 150 KWRK1 = 1 151 LWRK1 = LWRK 152 CALL MEMGET2('REAL','FC', KFC , NORBT*NORBT,WRK,KWRK1,LWRK1) 153 IF (LORB) THEN 154 CALL MEMGET2('REAL','FV', KFV , NORBT*NORBT,WRK,KWRK1,LWRK1) 155 CALL MEMGET2('REAL','QA', KQA , NORBT*NASHT,WRK,KWRK1,LWRK1) 156 CALL MEMGET2('REAL','QB', KQB , NORBT*NASHT,WRK,KWRK1,LWRK1) 157 END IF 158 CALL MEMGET2('REAL','H1', KH1 , NORBT*NORBT,WRK,KWRK1,LWRK1) 159 CALL MEMGET2('REAL','H2', KH2 , NORBT*NORBT,WRK,KWRK1,LWRK1) 160 IF (LCON) THEN 161 CALL MEMGET2('REAL','H2AC',KH2AC,NNASHX*NNASHX,WRK,KWRK1,LWRK1) 162 ELSE 163 CALL MEMGET2('REAL','H2AC',KH2AC,0 ,WRK,KWRK1,LWRK1) 164 END IF 165C 166 IF (IPRHSO.GT.20 .AND. LORB .AND. NASHT.GT.0) THEN 167 WRITE(LUPRI,'(/A)') 168 * ' First order density matrix:' 169 CALL OUTPUT(UDV,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 170 WRITE(LUPRI,'(/A)') 171 * ' Second order density matrix:' 172 CALL PRIAC2(PV,NASHT,LUPRI) 173 END IF 174C 175C 176C Read AO one-electron property integrals and transform to MO basis. 177C 178 IF (DO1) THEN 179 LABEL = WORD 180 LABEL(2:2) = '1' 181 KSYMP = -1 182 CALL PRPGET (LABEL,CMO,WRK(KH1),KSYMP,ANTSYM,WRK(KWRK1),LWRK1, 183 & IPRHSO) 184 IF (KSYMP.NE.KSYMOP) THEN 185 WRITE (LUPRI,'(/A/2A/A,2I5/A,F10.2)') 186 & 'FATAL ERROR: KSYMOP .ne. KSYMP from PRPGET', 187 & ' Property label : ',WORD, 188 & ' KSYMOP and KSYMP:',KSYMOP,KSYMP, 189 & ' ANTSYM :',ANTSYM 190 CALL QUIT('KSYMOP .ne. KSYMP from PRPGET') 191 END IF 192 END IF 193C 194C Print atomic and molecular property integrals, if desired 195C 196 IF (IPRHSO.GT.20 .AND. DO1) THEN 197 WRITE (LUPRI,'(/2A)')' Atomic property integrals:', LABEL 198 CALL OUTPUT(WRK(KWRK1),1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 199 WRITE (LUPRI,'(/2A)')' Molecular property integrals:', LABEL 200 CALL OUTPUT(WRK(KH1),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 201 END IF 202C IF (.NOT.DO2) GOTO 95 203 IF (.NOT.TDHF) CALL TRANSFORM_HSO(WORD,CMO,WRK,KWRK1,LWRK1) 204C 205 IF (TRPLET) THEN 206 TRPDEN=.FALSE. 207 CALL HSOFCK(WORD,WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB), 208 & WRK(KH2AC),UDV,PV,CMO,LORB,LCON,WRK,KWRK1,LWRK1) 209 ELSE 210 TRPDEN=.TRUE. 211 CALL MEMGET2('REAL','D',KD,N2ASHX,WRK,KWRK1,LWRK1) 212 CALL MEMGET2('REAL','P',KP,2*N2ASHX*N2ASHX,WRK,KWRK1,LWRK1) 213 CALL DZERO(WRK(KD),N2ASHX) 214 CALL DZERO(WRK(KP),2*N2ASHX*N2ASHX) 215 IF (TDHF) THEN 216 IF (NASHT.GE.1) THEN 217 CALL DUNIT(WRK(KD),NASHT) 218 END IF 219 ELSE 220C 221C Need new density matrices 222C 223 CALL MEMGET2('REAL','CREF',KCREF,MZCONF(1),WRK,KWRK1,LWRK1) 224 CALL GETREF(WRK(KCREF),MZCONF(1)) 225 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 226 & WRK(KCREF),WRK(KCREF), 227 & WRK(KD),WRK(KP),1,0,.FALSE.,.FALSE.,XINDX, 228 & WRK,KWRK1,LWRK1) 229 KP1=KP 230 KP2=KP+N2ASHX*N2ASHX 231 CALL MTRSP(N2ASHX,N2ASHX,WRK(KP1),N2ASHX,WRK(KP2),N2ASHX) 232 CALL MEMREL('HSOCTL<-RSPDM',WRK,KD,KCREF,KWRK1,LWRK1) 233 END IF 234 CALL HSOFCK(WORD,WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB), 235 & WRK(KH2AC),WRK(KD),WRK(KP),CMO,LORB,LCON,WRK,KWRK1,LWRK1) 236 CALL MEMREL('HSOCTL<-HSOFCK',WRK,KD,KP,KWRK1,LWRK1) 237 END IF 238 95 CONTINUE 239C 240C 241C Add the one-electron and two electron parts of the inactive Fock matrix 242C 243 IF (DO1) THEN 244 CALL DAXPY(N2ORBX,D1,WRK(KH1),1,WRK(KFC),1) 245 IF (TDHF.AND.NASHT.GT.0) THEN 246 CALL DAXPY(N2ORBX,DM1,WRK(KH1),1,WRK(KFV),1) 247 END IF 248 END IF 249C 250C 251C 252C Construct orbital property vector 253C 254 IF (LORB) THEN 255 IF (IPRHSO.GT.2) CALL TIMER('START ',ORBSTA,ORBTIM) 256 IF (TRPLET) THEN 257 CALL HSOORB(.TRUE.,1,DUMMY1, 258 * WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB), 259 * UDV ,GP,TRPDEN) 260 ELSE 261C 262C Quick fix for high spin matrices which do not fit neatly into HSOORB 263C 264 IF (TDHF.AND.NASHT.GT.0) TRPDEN=.FALSE. 265C 266 CALL HSOORB(.TRUE.,1,DUMMY1, 267 * WRK(KFC),WRK(KFV),WRK(KQA),WRK(KQB), 268 * WRK(KD),GP,TRPDEN) 269 END IF 270C CALL HSOORB(ONEIND,NSIM,FC,FCX,FVX,QAX,QBX,UDV, 271C * EVECS,TRPDEN) 272C 273C 274 IF (IPRHSO.GT.2) CALL TIMER('HSOORB',ORBSTA,ORBTIM) 275 END IF 276C 277C Compensate for the sign in HSOORB 278C 279 CALL DSCAL(KZWOPT,DM1,GP(1+KZCONF),1) 280 CALL DSCAL(KZWOPT,DM1,GP(1+KZCONF+KZVAR),1) 281C 282C Construct configuration property vector 283C 284 IF (LCON) THEN 285 IF (IPRHSO.GT.2) CALL TIMER('START ',ORBSTA,ORBTIM) 286 CALL HSOSIG(WRK(KFC),WRK(KH2AC), 287 * GP,XINDX,WRK(KWRK1),LWRK1) 288C 289C CALL HSOSIG(FC,H2AC, GP,XINDX,WRK,LWRK) 290C 291C 292 IF (IPRHSO.GT.2) CALL TIMER('HSOSIG',ORBSTA,ORBTIM) 293 END IF 294 IF (IPRHSO.GT.10) THEN 295 WRITE(LUPRI,'(//A//A)') 296 * ' Configuration property vector: ', 297 * ' Z part Y part' 298 CALL OUTPUT(GP,1,KZCONF,1,2,KZVAR,2,1,LUPRI) 299 WRITE(LUPRI,'(//A//A)') 300 * ' Orbital property vector: ', 301 * ' Z part Y part' 302 CALL OUTPUT(GP(1+KZCONF),1,KZWOPT,1,2,KZVAR,2,1,LUPRI) 303 ELSE IF (IPRHSO.GT.2) THEN 304 IF (LORB) CALL RSPPRO (GP(1+KZCONF),KZVAR,UDV,LUPRI) 305 IF (LCON) CALL RSPPRC (GP,KZCONF,KZVAR,LUPRI) 306 END IF 307C 308 IF (IPRHSO.GT.2) CALL TIMER('HSOCTL',HSOSTA,HSOTIM) 309 IF (X2GRAD) THEN 310 CALL MEMREL('X2GRAD',WRK,KFC,KFC,KWRK1,LWRK1) 311 CALL MEMGET2('INTE','MJWOP',KMJWOP,16*MAXWOP,WRK,KWRK1,LWRK1) 312 CALL HEADER('X2GRAD TEST FOR SPIN-ORBIT GRADIENT ELEMENTS',3) 313 CALL MEMGET2('REAL','GP2',KGP2,KZYVAR,WRK,KWRK1,LWRK1) 314 CALL DZERO(WRK(KGP2),KZYVAR) 315 CALL SETZY(WRK(KMJWOP)) 316 CALL HSOAL2 (WORD,WRK(KGP2),CMO,UDV,PV,XINDX,WRK(KMJWOP), 317 & WRK(KWRK1),LWRK1) 318 WRITE(LUPRI,'(/A)')' GP1 GP2' 319 WRITE(LUPRI,'(/A)')' ---------------------------------' 320 DMAXGP = D0 321 DO 77, K = 0,KZYVAR-1 322 WRITE(LUPRI,'(10X,2F14.8)') GP(1+K), WRK(KGP2+K) 323 DMAXGP = MAX(DMAXGP,ABS(GP(1+K)-WRK(KGP2+K))) 32477 CONTINUE 325 WRITE(LUPRI,'(/A)')' ---------------------------------' 326 WRITE(LUPRI,'(/A,E20.8)') 327 * 'LARGEST DIFFERENCE OF SPIN-ORBIT GRADIENT ELEMENTS' , 328 * DMAXGP 329 CALL MEMREL('X2GRAD',WRK,KGP2,KGP2,KWRK1,LWRK1) 330 END IF 331 CALL DSCAL(KZYVAR,HSOFAC,GP,1) 332C 333C Save on file 334C 335 REWIND LUHSO 336 CALL NEWLAB(WORD,LUHSO,LUERR) 337 CALL WRITT(LUHSO,KZYVAR,GP) 338 CALL MEMREL('HSOCTL',WRK,1,1,KWRK1,LWRK1) 339 CALL GPCLOSE(LUHSO,'KEEP') 340C 341 9999 CALL QEXIT('HSOCTL') 342 RETURN 343 END 344C /* Deck hsoinp */ 345 SUBROUTINE HSOINP(WORD) 346C 347#include "implicit.h" 348C 349#include "priunit.h" 350#include "infrsp.h" 351#include "wrkrsp.h" 352#include "infs0.h" 353#include "infpri.h" 354#include "infdim.h" 355#include "trhso.h" 356#include "infhso.h" 357C 358 LOGICAL NEWDEF 359 PARAMETER ( NTABLE = 8 ) 360 CHARACTER PROMPT*1, WORD*7, TABLE(NTABLE)*7, WORD1*7 361C 362 DATA TABLE /'.TESTZY', '.SO1ONL', '.SO2ONL', '.PRINT', 363 * '.OLDTRA', '.X2MAT ', '.A2MAT ', '.X2GRAD'/ 364C 365 NEWDEF = (WORD .EQ. '*SPIN-O') 366 ICHANG = 0 367 IF (NEWDEF) THEN 368 WORD1 = WORD 369 100 CONTINUE 370 READ (LUCMD, '(A7)') WORD 371 CALL UPCASE(WORD) 372 PROMPT = WORD(1:1) 373 IF (PROMPT .EQ. '!' .OR. PROMPT .EQ. '#') GO TO 100 374 IF (PROMPT .EQ. '.') THEN 375 ICHANG = ICHANG + 1 376 DO I=1, NTABLE 377 IF (TABLE(I) .EQ. WORD) THEN 378 GO TO (1,2,3,4,5,6,7,8), I 379 END IF 380 END DO 381 IF (WORD .EQ. '.OPTION') THEN 382 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 383 GO TO 100 384 END IF 385 WRITE (LUPRI,'(/,3A,/)') ' KEYWORD "',WORD, 386 * '" NOT RECOGNIZED IN HSOINP.' 387 CALL PRTAB(NTABLE,TABLE,WORD1//' input keywords',LUPRI) 388 CALL QUIT(' ILLEGAL KEYWORD IN HSOINP ') 389 1 CONTINUE 390 TESTZY = .TRUE. 391 GO TO 100 392 2 CONTINUE 393 DOSO2 = .FALSE. 394 GO TO 100 395 3 CONTINUE 396 DOSO1 = .FALSE. 397 GO TO 100 398 4 CONTINUE 399 READ(LUCMD,*)IPRHSO 400 GO TO 100 401 5 CONTINUE 402 OLDTRA = .TRUE. 403 GO TO 100 404 6 CONTINUE 405 X2MAT = .TRUE. 406 GO TO 100 407 7 CONTINUE 408 A2MAT = .TRUE. 409 GO TO 100 410 8 CONTINUE 411 X2GRAD = .TRUE. 412 GO TO 100 413 ELSE IF (PROMPT.EQ.'*') THEN 414 GO TO 300 415 ELSE 416 WRITE (LUPRI,'(/,3A,/)') ' PROMPT "',WORD, 417 * '" NOT RECOGNIZED IN RSPINP.' 418 CALL QUIT(' ILLEGAL PROMPT IN HSOINP ') 419 END IF 420 GO TO 100 421 END IF 422 300 CONTINUE 423 IF (ICHANG .GT. 0) THEN 424 CALL HEADER('CHANGES OF DEFAULTS FOR HSOINP:',0) 425 IF (TESTZY) WRITE(LUPRI,'(/A/A,L1)') 426 * ' Both parts of configuration property vector explicitely', 427 * ' TESTZY = ',TESTZY 428 IF (IPRHSO.NE.2) WRITE(LUPRI,'(/A/A,I5)') 429 * ' Print level in spin-orbit property vector calculation', 430 * ' IPRHSO = ',IPRHSO 431 IF (OLDTRA) WRITE(LUPRI,'(/A/A,L1)') 432 * ' Use existing two-electron spin-orbit integral file', 433 * ' OLDTRA = ',OLDTRA 434 IF (.NOT.DOSO1) WRITE(LUPRI,'(/A)') 435 * 'Skip one-particle part in HSOCTL' 436 IF (.NOT.DOSO2) WRITE(LUPRI,'(/A)') 437 * 'Skip two-particle part in HSOCTL' 438 IF (X2MAT) WRITE(LUPRI,'(/A)') 439 * 'X2MAT: Calculate full X2 matrix (quadratic response)' 440 IF (A2MAT) WRITE(LUPRI,'(/A)') 441 * 'A2MAT: not implemented' 442 END IF 443C 444C *** END OF HSOINP 445C 446 RETURN 447 END 448C /* Deck fsomu */ 449 SUBROUTINE FSOMU(ICI1,IDI1,H2,FCSO,FVSO,UDV,LORB,WRK,LWRK) 450C 451C Olav Vahtras 452C Apr 11, 1990 453C 454C CALCULATE ALL CONTRIBUTIONS TO INACTIVE AND ACTIVE SPIN-ORBIT 455C FOCK MATRICES FROM MULLIKEN DISTRIBUTIONS 456C 457C FCSO(P,Q) = SUM (K) 2*(KK|P^Q) - SUM (K) 3*(PK|K^Q) 458C - SUM (K) 3*(KQ|P^K) 459C 460C FVSO(P,Q) = SUM (X,Y) (XY|P^Q) D(XY) - SUM (K) 3/2*(PX|Y^Q) D(XY) 461C - SUM (K) 3/2*(XQ|P^Y) D(XY) 462C 463C 464#include "implicit.h" 465C 466C Used from common blocks: 467C 468C INFORB : NORBT,NISHT,NASHT,MULD2H 469C INFIND : ISMO,IOBTYP,ICH,ISMO,IASH,IORB,NISH,NASH,NOCC,NORB 470C INFHSO : 471C WRKRSP : KSYMOP 472C 473#include "maxash.h" 474#include "maxorb.h" 475#include "priunit.h" 476#include "inforb.h" 477#include "infind.h" 478#include "wrkrsp.h" 479#include "infrsp.h" 480#include "infpri.h" 481#include "orbtypdef.h" 482#include "trhso.h" 483#include "infhso.h" 484C 485 DIMENSION H2(NORBT,NORBT) 486 DIMENSION FCSO(NORBT,NORBT),FVSO(NORBT,NORBT) 487 DIMENSION UDV(NASHT,NASHT),WRK(*) 488 LOGICAL LORB 489C 490 PARAMETER (D1=1.0D0, D1P5=1.5D0, D2=2.0D0) 491C 492 CALL QENTER('FSOMU') 493C 494C 495C Local print level 496C 497 IPRFSO = 20 498C 499C Order (C,D) index such that C .ge. D 500C in inactive-active-secondary order (using ISW) 501C 502 IF (ISW(ICI1) .GE. ISW(IDI1)) THEN 503 ICI = ICI1 504 IDI = IDI1 505 DISFAC=D1 506 ELSE 507 ICI = IDI1 508 IDI = ICI1 509 DISFAC=-D1 510 END IF 511 IF (TRPLET) THEN 512 COUFAC=D1 513 ELSE 514 COUFAC=D2 515 END IF 516C 517C Find distribution type ITYPCD 518C 519 ITYPC = IOBTYP(ICI) 520 ITYPD = IOBTYP(IDI) 521 ITYPCD = IDBTYP(ITYPC,ITYPD) 522C 523C We only need secondary-occupied distributions 524C 525 IF (ITYPCD .EQ. JTSESE) GO TO 9999 526C 527 IF (ITYPC .EQ. JTACT) NCIW = ICH(ICI) 528 IF (ITYPD .EQ. JTACT) NDIW = ICH(IDI) 529C 530 ICSYM = ISMO(ICI) 531 IDSYM = ISMO(IDI) 532 ICDSYM = MULD2H(ICSYM,IDSYM) 533 KCDSYM = MULD2H(KSYMOP,ICDSYM) 534C 535 IF (IPRHSO.GT.IPRFSO) THEN 536 WRITE(LUPRI,'(/A)')' ------ Output from FSOMU ------' 537 WRITE(LUPRI,'(/A,2I5,5X,A,2X,A)')' Distribution CD',ICI1,IDI1, 538 * COBTYP(ITYPC),COBTYP(ITYPD) 539 WRITE(LUPRI,'(A,2I5)')' Reordered ',ICI,IDI 540 WRITE(LUPRI,'(A,2I5)')' Symmetry ',ICSYM,IDSYM 541 ENDIF 542C 543C 544C Inactive Fock matrix 545C 546C Direct terms 547C 548C FCSO(P,Q) = FCSO(P,Q) + SUM(K) 2*(KK|P^Q) 549C 550C here: 551C SUM(K) 2*(KK|C^D) -> FCSO(C,D) 552C - SUM(K) 2*(KK|C^D) -> FCSO(D,C) 553C 554 IF ( KSYMOP.EQ.ICDSYM ) THEN 555 DO 10 I=1,NISHT 556 IX=ISX(I) 557 WRK(I)=H2(IX,IX) 55810 CONTINUE 559 IF (IPRHSO.GT.IPRFSO) THEN 560 WRITE(LUPRI,'(/A)')' Inactive direct terms' 561 WRITE(LUPRI,'(A)')' (KK|CD) diagonal' 562 CALL OUTPUT(WRK(1),1,NISHT,1,1,NISHT,1,1,LUPRI) 563 END IF 564 FAC=2*DISFAC*DSUM(NISHT,WRK(1),1) 565 FCSO(ICI,IDI)=FCSO(ICI,IDI) + FAC 566 FCSO(IDI,ICI)=FCSO(IDI,ICI) - FAC 567 ENDIF 568C 569C Exchange terms: rearranged with D position inactive 570C 571C FCSO(P,Q) = FCSO(P,Q) + SUM(K) 3*(PK|Q^K) - 3*(QK|P^K) 572C 573C 574 IF (IPRHSO.GT.IPRFSO) THEN 575 WRITE(LUPRI,'(/A)')' Loop over symmetries' 576 END IF 577 DO 200 ISYM = 1,NSYM 578 IPSYM = ISYM 579 IOFFP = IORB(IPSYM) 580 NASHP = NASH(IPSYM) 581 NOCCP = NOCC(IPSYM) 582 NORBP = NORB(IPSYM) 583 IASHP=IASH(IPSYM) 584 IOFFPA=IOFFP+NISH(IPSYM) 585 IF (IPRHSO.GT.IPRFSO) THEN 586 WRITE(LUPRI,'(/A,I5)')' IPSYM',IPSYM 587 ENDIF 588 IF ( (NORBP.EQ.0) ) GO TO 200 589C 590 ICPSYM = MULD2H(ICSYM,IPSYM) 591 IDPSYM = MULD2H(IDSYM,IPSYM) 592C 593C For the case D inactive 594C 595 IF ((ITYPCD.EQ.JTININ).OR. 596 * (ITYPCD.EQ.JTACIN).OR. 597 * (ITYPCD.EQ.JTSEIN)) THEN 598C 599C here: 600C + 3*(PD|C^D) -> FCSO(P,C) 601C - 3*(PD|C^D) -> FCSO(C,P) 602C 603C such that (PC) is at most secondary-active 604C 605 IF (IPRHSO.GT.IPRFSO) THEN 606 IF (ICPSYM.EQ.KSYMOP .AND. 607 * (ITYPCD.EQ.JTSEIN .AND. NOCCP.GT.0 608 * .OR. 609 * ITYPCD.NE.JTSEIN ) 610 * .OR. 611 * IDPSYM.EQ.KSYMOP .AND. ITYPCD.EQ.JTININ) THEN 612 WRITE(LUPRI,'(/A)')' Inactive exchange terms' 613 END IF 614 END IF 615 FAC = 3*DISFAC 616 IF ( ICPSYM.EQ.KSYMOP ) THEN 617 IF (ITYPCD.EQ.JTSEIN) THEN 618 NDIMP = NOCCP 619 ELSE 620 NDIMP = NORBP 621 ENDIF 622 IF (NDIMP.GT.0) THEN 623 CALL DAXPY(NDIMP,FAC,H2(IOFFP+1,IDI),1, 624 * FCSO(IOFFP+1,ICI),1) 625 CALL DAXPY(NDIMP,-FAC,H2(IOFFP+1,IDI),1, 626 * FCSO(ICI,IOFFP+1),NORBT) 627 IF (IPRHSO.GT.IPRFSO) THEN 628 WRITE(LUPRI,'(A,I3,A,I3,I5)')' PC contribution ', 629 * IOFFP+1,':',IOFFP+NDIMP,ICI 630 CALL OUTPUT(H2(IOFFP+1,IDI),1,NDIMP,1,1, 631 * NORBT,NORBT,1,LUPRI) 632 ENDIF 633 END IF 634 END IF 635C 636C if both C and D are inactive we also have 637C 638C - 3*(PC|C^D) -> FCSO(P,D) 639C + 3*(PC|C^D) -> FCSO(D,P) 640C 641 IF ( IDPSYM.EQ.KSYMOP .AND. ITYPCD.EQ.JTININ ) THEN 642 IF (IPRHSO.GT.IPRFSO) THEN 643 WRITE(LUPRI,'(A,I3,A,I3,I5)')' PD contribution ', 644 * IOFFP+1,':',IOFFP+NDIMP,IDI 645 CALL OUTPUT(H2(IOFFP+1,ICI),1,NORBP,1,1, 646 * NORBT,NORBT,1,LUPRI) 647 ENDIF 648 CALL DAXPY(NORBP,-FAC,H2(IOFFP+1,ICI),1, 649 * FCSO(IOFFP+1,IDI),1) 650 CALL DAXPY(NORBP,FAC,H2(IOFFP+1,ICI),1, 651 * FCSO(IDI,IOFFP+1),NORBT) 652 ENDIF 653 ENDIF 654C 655 IF (LORB) THEN 656C 657C Active Fock matrix 658C 659C Direct terms: 660C 661C FVSO(P,Q) = FVSO(P,Q)+ SUM(X,Y) (XY|P^Q)*D(XY) 662C 663C 664C here: 665C + SUM(X,Y) (XY|C^D)*D(XY) -> FVSO(C,D) 666C - SUM(X,Y) (XY|C^D)*D(XY) -> FVSO(D,C) 667C 668C where the sum is taken over diagonal symmetry blocks (X,Y) 669C 670 FAC = DISFAC*COUFAC 671 IXSYM = ISYM 672 IASHX = IASH(IXSYM) 673 NASHX = NASH(IXSYM) 674 IOFFXA = IORB(IXSYM) + NISH(IXSYM) 675 IF (KSYMOP.EQ.ICDSYM) THEN 676 IF (IPRHSO.GT.IPRFSO) THEN 677 WRITE(LUPRI,'(/A)')' Active direct terms' 678 END IF 679 DO 20 IX=1,NASHX 680 WRK(IX) = DDOT(NASHX,H2(IOFFXA+1,IOFFXA+IX),1, 681 * UDV(IASHX+1,IASHX+IX),1) 682 20 CONTINUE 683 XYSUM = DSUM(NASHX,WRK(1),1) 684 FVSO(ICI,IDI) = FVSO(ICI,IDI) + FAC*XYSUM 685 FVSO(IDI,ICI) = FVSO(IDI,ICI) - FAC*XYSUM 686 END IF 687C 688C Exchange terms: rearranged with C position active 689C 690C FVSO(P,Q) = FVSO(P,Q) - SUM(X,Y) 3/2*(PX|Y^Q)*D(X,Y) 691C + SUM(X,Y) 3/2*(QX|Y^P)*D(X,Y) 692C 693C 694C here: 695C - SUM(X) 3/2*(PX|C^D)*D(X,C) -> FVSO(P,D) 696C + SUM(X) 3/2*(PX|C^D)*D(X,C) -> FVSO(D,P) 697C 698C for active-active and active-inactive distributions 699C 700 FAC=D1P5*DISFAC 701 IXSYM = MULD2H(IPSYM,KCDSYM) 702 IF (IPRHSO.GT.IPRFSO) THEN 703 IF (IXSYM.EQ.ICSYM .AND. 704 * (ITYPCD.EQ.JTACIN .OR. ITYPCD.EQ.JTACAC) 705 * .OR. 706 * IXSYM.EQ.IDSYM .AND. 707 * (ITYPCD.EQ.JTACAC .OR. 708 * NOCCP.GT.0 .AND. ITYPCD.EQ.JTSEAC)) THEN 709 WRITE(LUPRI,'(/A)')' Active exchange terms' 710 END IF 711 END IF 712 IF (IXSYM.EQ.ICSYM) THEN 713 IOFFXA = IORB(IXSYM) + NISH(IXSYM) 714 NASHX = NASH(IXSYM) 715 IASHX = IASH(IXSYM) 716 IF (ITYPCD.EQ.JTACIN .OR. ITYPCD.EQ.JTACAC) THEN 717 CALL DGEMM('N','N',NORBP,1,NASHX,1.0D0, 718 & H2(IOFFP+1,IOFFXA+1),NORBT, 719 & UDV(IASHX+1,NCIW),NASHT,0.0D0,WRK(1),NORBP) 720 CALL DAXPY(NORBP,-FAC,WRK(1),1, 721 * FVSO(IOFFP+1,IDI),1) 722 CALL DAXPY(NORBP,FAC,WRK(1),1, 723 * FVSO(IDI,IOFFP+1),NORBT) 724 IF (IPRHSO.GT.IPRFSO) THEN 725 WRITE(LUPRI,'(A,I3,A,I3,I5)')' PD contribution ', 726 * IOFFP+1,':',IOFFP+NORBP,IDI 727 CALL OUTPUT(WRK(1),1,NORBP,1,1,NORBP,1,1,LUPRI) 728 ENDIF 729 END IF 730 END IF 731C 732C Exchange terms: rearranged with D position active 733C 734C FVSO(P,Q) = FVSO(P,Q) + SUM(X,Y) 3/2*(PX|Q^Y)*D(X,Y) 735C - SUM(X,Y) 3/2*(QX|P^Y)*D(X,Y) 736C 737C 738C here: 739C SUM(X) 3/2*(PX|C^D)*D(X,D) -> FVSO(P,C) 740C - SUM(X) 3/2*(PX|C^D)*D(X,D) -> FVSO(C,P) 741C 742C such that (CP) is at most secondary-active 743C 744C for active-active and secondary-active distributions 745C 746 IF (IXSYM.EQ.IDSYM) THEN 747 IOFFXA = IORB(IXSYM) + NISH(IXSYM) 748 NASHX = NASH(IXSYM) 749 IASHX = IASH(IXSYM) 750 IF (ITYPCD.EQ.JTACAC .OR. ITYPCD.EQ.JTSEAC) THEN 751 IF (ITYPCD.EQ.JTSEAC) THEN 752 NDIMP = NOCCP 753 ELSE 754 NDIMP = NORBP 755 END IF 756 IF (NDIMP .GT. 0) THEN 757 CALL DGEMM('N','N',NDIMP,1,NASHX,1.0D0, 758 & H2(IOFFP+1,IOFFXA+1),NORBT, 759 & UDV(IASHX+1,NDIW),NASHT,0.0D0, 760 & WRK(1),NDIMP) 761 CALL DAXPY(NDIMP,FAC,WRK(1),1, 762 * FVSO(IOFFP+1,ICI),1) 763 CALL DAXPY(NDIMP,-FAC,WRK(1),1, 764 * FVSO(ICI,IOFFP+1),NORBT) 765 IF (IPRHSO.GT.IPRFSO) THEN 766 WRITE(LUPRI,'(A,I3,A,I3,I5)')' PC contribution ', 767 * IOFFP+1,':',IOFFP+NDIMP,ICI 768 CALL OUTPUT(WRK(1),NDIMP,1,1,1,NDIMP,1,1,LUPRI) 769 ENDIF 770 END IF 771 END IF 772 END IF 773C 774 END IF 775C (RSPCI) 776200 CONTINUE 7779999 CALL QEXIT('FSOMU') 778 RETURN 779 END 780C /* Deck qsomu */ 781 SUBROUTINE QSOMU(ICI,IDI,QASO,QBSO, 782 * H2,PVX,PV12,PV21, 783 * WRK,LWRK) 784C 785C Olav Vahtras 786C Apr 11, 1990 787C 788C Purpose: 789C Calculate all contributions to QA and QB spin-orbit matrices 790C from Mulliken (**|C^D) integral distributions 791C 792C In general: 793C 794C QBSO(P,Q) = SUM(X,Y,W) (PW|X^Y)*( 2*PV(++)(Q,W,X,Y) + PV(--)(Q,W,X,Y) ) 795C + SUM(X,Y,W) (XY|P^W)*( 2*PV(--)(X,Y,Q,W) + PV(++)(X,Y,Q,W) ) 796C 797C QASO(P,Q) = SUM(X,Y,W) (WP|X^Y)*( 2*PV(++)(W,Q,X,Y) + PV(--)(W,Q,X,Y) ) 798C + SUM(X,Y,W) (XY|W^P)*( 2*PV(--)(X,Y,W,Q) + PV(++)(X,Y,W,Q) ) 799C 800C where PV(++)(P,Q,R,S) = <0| e(+,+)(pqrs) |0> 801C and PV(--)(P,Q,R,S) = <0| e(-,-)(pqrs) |0> 802C 803#include "implicit.h" 804C 805#include "maxash.h" 806#include "maxorb.h" 807#include "priunit.h" 808#include "inforb.h" 809#include "infind.h" 810#include "wrkrsp.h" 811#include "infrsp.h" 812#include "infpri.h" 813#include "orbtypdef.h" 814#include "trhso.h" 815#include "infhso.h" 816C 817C 818 DIMENSION QASO(NORBT,NASHT),QBSO(NORBT,NASHT) 819 DIMENSION PV12(NASHT,NASHT),PV21(NASHT,NASHT) 820 DIMENSION H2(NORBT,NORBT),PVX(*),WRK(*) 821C 822 PARAMETER(D2=2.0D0) 823C 824 CALL QENTER('QSOMU') 825C 826C Local print level 827C 828 IPRQSO = 20 829C 830 IF (LWRK.LT.N2ASHX) CALL ERRWRK('QSOMU',N2ASHX,LWRK) 831C 832C 833C Symmetry to type order 834C 835 ICIW = ISW(ICI) 836 IDIW = ISW(IDI) 837C 838C Orbital type, at least one has to be active to contribute 839C 840 ITYPC = IOBTYP(ICI) 841 ITYPD = IOBTYP(IDI) 842 IF (ITYPC.NE.JTACT .AND. ITYPD.NE.JTACT) GO TO 9999 843C 844C Distribution type and symmetry 845C 846 ICSYM = ISMO(ICI) 847 IDSYM = ISMO(IDI) 848 ICDSYM = MULD2H(ICSYM,IDSYM) 849 KCDSYM = MULD2H(KSYMOP,ICDSYM) 850 ITYPCD=IDBTYP(ITYPC,ITYPD) 851C 852C Order within actives for the case that C or D are active 853C 854 IF (ITYPC .EQ. JTACT) NCIW = ICIW - NISHT 855 IF (ITYPD .EQ. JTACT) NDIW = IDIW - NISHT 856C 857 IF (IPRHSO.GT.IPRQSO) THEN 858 WRITE(LUPRI,'(/A)') ' ------ Output from QSOMU ------' 859 WRITE(LUPRI,'(/A,2I5,5X,2A)')' Distribution CD',ICI,IDI, 860 * COBTYP(ITYPC),COBTYP(ITYPD) 861 WRITE(LUPRI,'(A,2I5)') ' Symmetry ',ICSYM,IDSYM 862 ENDIF 863C 864C 865 IPP = 1 866 IMM = N2ASHX*N2ASHX+1 867C 868C Both C and D are active: 869C 870 IF (ITYPCD.EQ.JTACAC) THEN 871C 872C Get (C,D) density distributions in the form 2PV(++)+PV(--) 873C 874 NCDOFF = (NCIW-1 + (NDIW-1)*NASHT)*N2ASHX 875 CALL DCOPY(N2ASHX,PVX(NCDOFF+IMM),1,PV12,1) 876 CALL DAXPY(N2ASHX,D2,PVX(NCDOFF+IPP),1,PV12,1) 877 IF (IPRHSO.GT.IPRQSO) THEN 878 WRITE(LUPRI,'(/A)')' Active diagonal terms' 879 WRITE(LUPRI,'(/A,2I5)') 880 * ' Total CD density distribution' ,NCIW,NDIW 881 CALL OUTPUT(PV12,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 882 END IF 883C 884C Get (D,C) density distributions 885C 886 NDCOFF = (NDIW-1 + (NCIW-1)*NASHT)*N2ASHX 887 CALL DCOPY(N2ASHX,PVX(NDCOFF+IMM),1,PV21,1) 888 CALL DAXPY(N2ASHX,D2,PVX(NDCOFF+IPP),1,PV21,1) 889 IF (IPRHSO.GT.IPRQSO) THEN 890 WRITE(LUPRI,'(/A,2I5)') 891 * ' Total DC density distribution' ,NDIW,NCIW 892 CALL OUTPUT(PV21,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 893 END IF 894C 895C Add contibutions to QASO and QBSO from (C,D) and (D,C) density 896C distributions 897C 898C 899 IF (IPRHSO.GT.IPRQSO) THEN 900 WRITE(LUPRI,'(/A)') ' Loop over symmetry blocks ' 901 END IF 902 DO 100 IPSYM = 1,NSYM 903 IWSYM = MULD2H(IPSYM,KCDSYM) 904 IQSYM = MULD2H(IWSYM,ICDSYM) 905 NORBP = NORB(IPSYM) 906 NASHP = NASH(IPSYM) 907 NASHQ = NASH(IQSYM) 908 NASHW = NASH(IWSYM) 909 IF (IPRHSO.GT.IPRQSO) THEN 910 WRITE(LUPRI,'(/A,I5)')' IPSYM',IPSYM 911 END IF 912 IF ( NORBP.GT.0 .AND. NASHW.GT.0 .AND. NASHQ.GT.0 ) 913 * THEN 914 IOFFP = IORB(IPSYM) 915 IOFFW = IORB(IWSYM) 916 IOFFPA = IOFFP + NISH(IPSYM) 917 IOFFWA = IOFFW + NISH(IWSYM) 918 IASHP = IASH(IPSYM) 919 IASHQ = IASH(IQSYM) 920 IASHW = IASH(IWSYM) 921 IF (IPRHSO.GT.IPRQSO) THEN 922 WRITE(LUPRI,'(/A,2I5)') 923 * ' Integral symmetry block PW',IPSYM,IWSYM 924 CALL OUTPUT(H2,IOFFP+1,IOFFP+NORBP, 925 * IOFFWA+1,IOFFWA+NASHW, 926 * NORBT,NORBT,1,LUPRI) 927 WRITE(LUPRI,'(/A,2I5)') 928 * ' Density symmetry block QW',IQSYM,IWSYM 929 CALL OUTPUT(PV12,IASHQ+1,IASHQ+NASHQ, 930 * IASHW+1,IASHW+NASHW, 931 * NASHT,NASHT,1,LUPRI) 932 END IF 933C 934C QBSO(P,Q) += SUM(W) (PW|C^D)*(2*PV(++)(Q,W,C,D) + PV(--)(Q,W,C,D)) 935C 936C 937 CALL DGEMM('N','T',NORBP,NASHQ,NASHW,1.D0, 938 & H2(IOFFP+1,IOFFWA+1),NORBT, 939 & PV12(IASHQ+1,IASHW+1),NASHT,1.D0, 940 & QBSO(IOFFP+1,IASHQ+1),NORBT) 941C 942C QBSO(P,Q) += SUM(W) (PW|D^C)*(2*PV(++)(Q,W,D,C) + PV(--)(Q,W,D,C)) 943C ( - ... C^D ) 944C 945 CALL DGEMM('N','T',NORBP,NASHQ,NASHW,-1.D0, 946 & H2(IOFFP+1,IOFFWA+1),NORBT, 947 & PV21(IASHQ+1,IASHW+1),NASHT,1.D0, 948 & QBSO(IOFFP+1,IASHQ+1),NORBT) 949C 950C QASO(P,Q) += SUM(W) (WP|C^D)*(2*PV(++)(W,Q,C,D) + PV(--)(W,Q,C,D)) 951C 952 CALL DGEMM('T','N',NORBP,NASHQ,NASHW,1.D0, 953 & H2(IOFFWA+1,IOFFP+1),NORBT, 954 & PV12(IASHW+1,IASHQ+1),NASHT,1.D0, 955 & QASO(IOFFP+1,IASHQ+1),NORBT) 956C 957C QASO(P,Q) += SUM(W) (WP|D^C)*(2*PV(++)(W,Q,D,C) + PV(--)(W,Q,D,C)) 958C ( - ... C^D ) 959C 960 CALL DGEMM('T','N',NORBP,NASHQ,NASHW,-1.D0, 961 & H2(IOFFWA+1,IOFFP+1),NORBT, 962 & PV21(IASHW+1,IASHQ+1),NASHT,1.D0, 963 & QASO(IOFFP+1,IASHQ+1),NORBT) 964C 965 END IF 966 100 CONTINUE 967 END IF 968C 969C Following integral distributions contribute if either 970C C or D is active. 971C 972C 973C QBSO(C,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,Q,D)+PV(++)(X,Y,Q,D)) 974C 975 IF (ITYPD.EQ.JTACT) THEN 976 IQSYM = MULD2H(KSYMOP,ICSYM) 977 NASHQ = NASH(IQSYM) 978 IASHQ = IASH(IQSYM) 979 IF (NASHQ.GT.0) THEN 980 IF (IPRHSO.GT.IPRQSO) THEN 981 WRITE(LUPRI,'(/A,I5)') 982 * ' Loop over actives in symmetry',IQSYM 983 END IF 984 DO 200 IQ=1,NASHQ 985C 986C For each Q read a new (Q,D) density distribution in the form 987C 2PV(--) + PV(++) 988C 989 IF (IPRHSO.GT.IPRQSO) THEN 990 WRITE(LUPRI,'(/A,2I5)')' IQ',IASHQ+IQ 991 END IF 992 NQDOFF = (IASHQ+IQ-1 + (NDIW-1)*NASHT)*N2ASHX 993 CALL DCOPY(N2ASHX,PVX(NQDOFF+IPP),1,PV12,1) 994 CALL DAXPY(N2ASHX,D2,PVX(NQDOFF+IMM),1,PV12,1) 995 IF (IPRHSO.GT.IPRQSO) THEN 996 WRITE(LUPRI,'(/A,2I5)') 997 * ' Total QD density distribution', IASHQ+IQ,NDIW 998 CALL OUTPUT(PV12,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 999 END IF 1000C 1001C 1002 DO 20 IYSYM = 1,NSYM 1003 NASHY = NASH(IYSYM) 1004 IASHY = IASH(IYSYM) 1005 IOFFYA = IORB(IYSYM) + NISH(IYSYM) 1006 IXSYM = MULD2H(IYSYM,KCDSYM) 1007 NASHX = NASH(IXSYM) 1008 IASHX = IASH(IXSYM) 1009 IOFFXA = IORB(IXSYM) + NISH(IXSYM) 1010 IF (NASHX.GT.0 .AND .NASHY .GT.0) THEN 1011C 1012C QBSO(C,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,Q,D)+PV(++)(X,Y,Q,D)) 1013C 1014C 1015C QASO(C,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,D,Q)+PV(++)(X,Y,D,Q)) 1016C ( - ... C^D ) (Q,D) (Q,D) 1017C 1018 IF (IPRHSO.GT.IPRQSO) THEN 1019 WRITE(LUPRI,'(/A,2I5)') 1020 * ' Density symmetry block XY',IXSYM,IYSYM 1021 CALL OUTPUT(PV12,IASHX+1,IASHX+NASHX, 1022 * IASHY+1,IASHY+NASHY, 1023 * NASHT,NASHT,1,LUPRI) 1024 END IF 1025 DO 21 IY=1,NASHY 1026 WRK(IY) = DDOT(NASHX,H2(IOFFXA+1,IOFFYA+IY),1, 1027 * PV12(IASHX+1,IASHY+IY),1) 1028 21 CONTINUE 1029 XYSUM = DSUM(NASHY,WRK,1) 1030 QBSO(ICI,IASHQ+IQ) 1031 * = QBSO(ICI,IASHQ+IQ) + XYSUM 1032 QASO(ICI,IASHQ+IQ) 1033 * = QASO(ICI,IASHQ+IQ) - XYSUM 1034 END IF 1035 20 CONTINUE 1036 200 CONTINUE 1037 END IF 1038 ENDIF 1039C 1040C QBSO(D,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,Q,C)+PV(++)(X,Y,Q,C)) 1041C 1042 IF (ITYPC.EQ.JTACT) THEN 1043 IQSYM = MULD2H(KSYMOP,IDSYM) 1044 NASHQ = NASH(IQSYM) 1045 IASHQ = IASH(IQSYM) 1046 IF (NASHQ.GT.0) THEN 1047 IF (IPRHSO.GT.IPRQSO) THEN 1048 WRITE(LUPRI,'(/A,I5)') 1049 * ' Loop over actives in symmetry',IQSYM 1050 END IF 1051 DO 300 IQ = 1,NASHQ 1052 IF (IPRHSO.GT.IPRQSO) THEN 1053 WRITE(LUPRI,'(/A,2I5)')' IQ',IASHQ+IQ 1054 END IF 1055C 1056C For each Q read a new (Q,C) density distribution in the form 1057C 2PV(--) + PV(++) 1058C 1059 NQCOFF = (IASHQ+IQ-1 + (NCIW-1)*NASHT)*N2ASHX 1060 CALL DCOPY(N2ASHX,PVX(NQCOFF+IPP),1,PV21,1) 1061 CALL DAXPY(N2ASHX,D2,PVX(NQCOFF+IMM),1,PV21,1) 1062 IF (IPRHSO.GT.IPRQSO) THEN 1063 WRITE(LUPRI,'(/A,2I5)') 1064 * ' Total QC density distribution', IASHQ+IQ,NCIW 1065 CALL OUTPUT(PV21,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 1066 END IF 1067C 1068C 1069 DO 30 IYSYM = 1,NSYM 1070 NASHY = NASH(IYSYM) 1071 IASHY = IASH(IYSYM) 1072 IOFFYA = IORB(IYSYM) + NISH(IYSYM) 1073 IXSYM = MULD2H(IYSYM,KCDSYM) 1074 NASHX = NASH(IXSYM) 1075 IASHX = IASH(IXSYM) 1076 IOFFXA = IORB(IXSYM) + NISH(IXSYM) 1077 IF (NASHX.GT.0 .AND. NASHY.GT.0) THEN 1078C 1079C QBSO(D,Q) += SUM(X,Y) (XY|D^C)*(2*PV(--)(X,Y,Q,C)+PV(++)(X,Y,Q,C)) 1080C ( - ... C^D ) 1081C 1082C 1083C QASO(D,Q) += SUM(X,Y) (XY|C^D)*(2*PV(--)(X,Y,C,Q)+PV(++)(X,Y,C,Q)) 1084C (Q,C) (Q,C) 1085 IF (IPRHSO.GT.IPRQSO) THEN 1086 WRITE(LUPRI,'(/A,2I5)') 1087 * ' Density symmetry block XY',IXSYM,IYSYM 1088 CALL OUTPUT(PV21,IASHX+1,IASHX+NASHX, 1089 * IASHY+1,IASHY+NASHY, 1090 * NASHT,NASHT,1,LUPRI) 1091 END IF 1092 DO 31 IY=1,NASHY 1093 WRK(IY) = DDOT(NASHX,H2(IOFFXA+1,IOFFYA+IY),1, 1094 * PV21(IASHX+1,IASHY+IY),1) 1095 31 CONTINUE 1096 XYSUM = DSUM(NASHY,WRK,1) 1097 QBSO(IDI,IASHQ+IQ) 1098 * = QBSO(IDI,IASHQ+IQ) - XYSUM 1099 QASO(IDI,IASHQ+IQ) 1100 * = QASO(IDI,IASHQ+IQ) + XYSUM 1101 END IF 1102 30 CONTINUE 1103C 1104 300 CONTINUE 1105 END IF 1106 END IF 1107C 1108 9999 CALL QEXIT('QSOMU') 1109 RETURN 1110 END 1111C /* Deck hsosig */ 1112 SUBROUTINE HSOSIG(FC,H2AC, 1113 * GP,XINDX,WRK,LWRK) 1114C 1115C Calculate the configuration part of the spin-orbit property vector 1116C 1117C ( <J,HSO,0> ) 1118C GP(J) = ( ) 1119C (-<0,HSO,J> ) 1120C 1121C 1122C H2AC contain active two-electron spin-orbit integrals in doubly 1123C triangularly packed form 1124C 1125C Olav Vahtras 1126C Sep 14, 1990 1127#include "implicit.h" 1128C 1129 DIMENSION FC(NORBT,NORBT) 1130 DIMENSION H2AC(NNASHX,NNASHX) 1131 DIMENSION GP(KZYVAR),XINDX(*),WRK(*) 1132C 1133#include "maxorb.h" 1134#include "maxash.h" 1135#include "priunit.h" 1136#include "wrkrsp.h" 1137#include "infrsp.h" 1138#include "infopt.h" 1139#include "inforb.h" 1140#include "infind.h" 1141#include "infpri.h" 1142#include "infdim.h" 1143#include "cbgetdis.h" 1144#include "trhso.h" 1145#include "infhso.h" 1146C 1147C 1148 PARAMETER ( DM1 = -1.0D0 ) 1149C 1150 CALL QENTER('HSOSIG') 1151C 1152C Local print level 1153C 1154 IPRLOC = 20 1155 IF (IPRHSO.GT.IPRLOC) THEN 1156 WRITE(LUPRI,'(/A)') ' ------ Output from HSOSIG ------' 1157 END IF 1158C 1159C ALLOCATE WORK SPACE 1160C 1161 IF (IREFSY .EQ. KSYMST) THEN 1162 NDREF = KZCONF 1163C ... if KZCONF .ne. NCREF, we need CREF in determinants 1164 ELSE 1165 NDREF = NCREF 1166 END IF 1167 KFREE = 1 1168 LFREE = LWRK 1169 CALL MEMGET2('REAL','FCAC' ,KFCAC,N2ASHX,WRK,KFREE,LFREE) 1170 CALL MEMGET2('REAL','REFCO',KREFCO,NDREF,WRK,KFREE,LFREE) 1171 CALL MEMGET2('REAL','HSQSQ',KHSQSQ,N2ASHX*N2ASHX,WRK,KFREE,LFREE) 1172 CALL DZERO(WRK(KFCAC),N2ASHX) 1173 CALL DZERO(WRK(KREFCO),NDREF) 1174 CALL DZERO(WRK(KHSQSQ),N2ASHX*N2ASHX) 1175 IF (TESTZY) THEN 1176 CALL MEMGET2('REAL','HSQTR',KHSQTR,N2ASHX*N2ASHX, 1177 & WRK,KFREE,LFREE) 1178 CALL DZERO(WRK(KHSQTR),N2ASHX*N2ASHX) 1179 END IF 1180 IADINT = -1 1181C 1182C 1183C Obtain square packed combined two-electron integrals 1184C 1185 CALL H2ACSO(H2AC,WRK(KHSQSQ)) 1186C 1187 ISPIN1 = 0 1188 ISPIN2 = 1 1189C 1190 CALL GETREF(WRK(KREFCO),NDREF) 1191C 1192C CREATE Z PART <0,HSO,J> OF SPIN-ORBIT GP VECTOR 1193C 1194C Get FCAC matrix for Z sigma vector 1195C (note: CISIGD requires UFCAC(I,J) = FCXAC(J,I)) 1196C 1197 DO 220 IW = 1,NASHT 1198 IX = ISX(NISHT+IW) 1199 DO 230 JW = 1,NASHT 1200 JX = ISX(NISHT+JW) 1201 IJW = (IW-1) * NASHT + JW 1202 WRK(KFCAC-1+IJW) = FC(IX,JX) 1203 230 CONTINUE 1204 220 CONTINUE 1205C 1206 IF (IPRHSO.GT.IPRLOC) THEN 1207 WRITE(LUPRI,'(/A)')' Inactive Fock matrix' 1208 CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1209 WRITE(LUPRI,'(/A)')' Active part of inactive Fock matrix' 1210 CALL OUTPUT(WRK(KFCAC),1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 1211 END IF 1212 DISTYP = 6 1213C 1214 KSYMST = MULD2H(IREFSY,KSYMSO) 1215 CALL CISIGD(IREFSY, KSYMST, NDREF, KZCONF, WRK(KREFCO),GP, 1216 * WRK(KFCAC), WRK(KHSQSQ), .FALSE., .FALSE., 1217 * XINDX, ISPIN1, ISPIN2, WRK(KFREE), LFREE) 1218C CALL RSPSIG(ICSYM,IHCSYM,NCDET,NHCDET,C,HC,UFCAC,H2AC,IFLAG, 1219C * NOH2,WORK,KFREE,LFREE) 1220C 1221C IF ((NDREF.NE.NCDET).OR.(KZCONF.NE.NHCDET)) THEN 1222C WRITE(LUPRI,'(/2(A,I5),/3(A,I5))') 1223C * ' NUMBER OF REFERENCE DETERMINANTS ,NDREF:',NDREF, 1224C * ' CALCULATED NUMBER ,NCDET:',NCDET, 1225C * ' NUMBER OF DETERMINANTS FOR SYMMETRY',KSYMOP, 1226C * ' IS:',KZCONF,' CALCULATED NUMBER, NHCDET:',NHCDET 1227C CALL QUIT(' H2XSIG, INCORRECT CALCULATION OF DETERMINANTS') 1228C END IF 1229 IF (IPRHSO.GT.IPRLOC) THEN 1230 WRITE(LUPRI,'(/A)') 1231 * ' Configuration part of spin-orbit property vector: Z part' 1232 CALL OUTPUT(GP,1,KZCONF,1,1,KZCONF,1,1,LUPRI) 1233 END IF 1234C 1235C CREATE Y PART <0,HSO,J> OF SPIN-ORBIT GP VECTOR 1236C 1237C Normal procedure: copy the Z part 1238C Test procedure: do it explicitely 1239C 1240C Get transposed FCXAC matrix for Y sigma vector 1241C (note: CISIGD requires UFCAC(I,J) = FCXAC(J,I)) 1242C 1243 IF (TESTZY) THEN 1244C 1245C ISPIN1 = 0 1246C ISPIN2 = 1 1247C 1248 DO 120 IW = 1,NASHT 1249 IX = ISX(NISHT+IW) 1250 DO 130 JW = 1,NASHT 1251 JX = ISX(NISHT+JW) 1252 IJW = (IW-1) * NASHT + JW 1253 WRK(KFCAC-1+IJW) = FC(JX,IX) 1254 130 CONTINUE 1255 120 CONTINUE 1256 IF (IPRHSO.GT.IPRLOC) THEN 1257 WRITE(LUPRI,'(/A)')' Inactive Fock matrix' 1258 CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1259 WRITE(LUPRI,'(/A)')' Active part of inactive Fock matrix' 1260 CALL OUTPUT(WRK(KFCAC),1,NASHT,1,NASHT,NASHT,NASHT,1, 1261 & LUPRI) 1262 END IF 1263C 1264C Get transposed integrals 1265C 1266C CALL MTRSP(N2ASHX,N2ASHX,WRK(KHSQSQ),N2ASHX, 1267C * WRK(KHSQTR),N2ASHX) 1268 CALL DAXPY(N2ASHX*N2ASHX,DM1,WRK(KHSQSQ),1,WRK(KHSQTR),1) 1269 DISTYP = 6 1270 CALL CISIGD(IREFSY, KSYMST, NDREF, KZCONF, WRK(KREFCO), 1271 * GP(KZVAR+1), 1272 * WRK(KFCAC), WRK(KHSQTR), .FALSE., .FALSE., 1273 * XINDX, ISPIN1, ISPIN2, WRK(KFREE), LFREE) 1274 CALL DSCAL(KZCONF,DM1,GP(1+KZVAR),1) 1275 IF (IPRHSO.GT.IPRLOC) THEN 1276 WRITE(LUPRI,'(/A)') 1277 & ' Configuration part of spin-orbit property vector: Y part' 1278 CALL OUTPUT(GP(1+KZVAR),1,KZCONF,1,1,KZCONF,1,1,LUPRI) 1279 END IF 1280C 1281 ELSE 1282 CALL DCOPY(KZCONF,GP,1,GP(1+KZVAR),1) 1283 END IF 1284 100 CONTINUE 1285 CALL MEMREL('HSOSIG',WRK,1,1,KFREE,LFREE) 1286 CALL QEXIT('HSOSIG') 1287 RETURN 1288 END 1289C /* Deck h2acso */ 1290 SUBROUTINE H2ACSO(H2ACPK,H2SOAC) 1291C 1292C 20-Jun-1990 hjaaj 1293C 1294#include "implicit.h" 1295 DIMENSION H2ACPK(NNASHX,NNASHX), H2SOAC(N2ASHX,NASHT,NASHT) 1296C 1297 PARAMETER ( DM1 = -1.0D0 ) 1298C 1299C Used from common blocks: 1300C INFORB : NASHT,NNASHX,N2ASHX 1301C INFIND : IROW() 1302C 1303#include "maxash.h" 1304#include "maxorb.h" 1305#include "inforb.h" 1306#include "infind.h" 1307C 1308 CALL QENTER('H2ACSO') 1309C 1310C 1. Unpack H2ACPK into H2SOAC 1311C 1312 DO 140 K = 1,NASHT 1313 DO 130 L = 1,K-1 1314 KL = IROW(K) + L 1315 CALL DSPTSI(NASHT,H2ACPK(1,KL),H2SOAC(1,K,L)) 1316 CALL DAXPY(N2ASHX,DM1,H2SOAC(1,K,L),1,H2SOAC(1,L,K),1) 1317 130 CONTINUE 1318 CALL DZERO(H2SOAC(1,K,K),N2ASHX) 1319 140 CONTINUE 1320 CALL QEXIT('H2ACSO') 1321 RETURN 1322 END 1323 SUBROUTINE TRANSFORM_HSO(WORD,CMO,WRK,KWRK1,LWRK1) 1324#include "implicit.h" 1325 CHARACTER*8 WORD 1326 DIMENSION WRK(*), CMO(*) 1327#include "dummy.h" 1328#include "rspprp.h" 1329#include "infhso.h" 1330#include "infhyp.h" 1331#include "infpp.h" 1332#include "infsmo.h" 1333#include "inftap.h" 1334#include "priunit.h" 1335#include "trhso.h" 1336#include "wrkrsp.h" 1337 LOGICAL FOPEN 1338 LOGICAL MOEXIS, SO2TRA 1339C 1340C Set spin-orbit component ILXYZ for TRHSO 1341C 1342 CALL QENTER('TRANSFORM_HSO') 1343 IF (WORD(1:1).EQ.'X') THEN 1344 ILXYZ = 1 1345 ELSE IF (WORD(1:1).EQ.'Y') THEN 1346 ILXYZ = 2 1347 ELSE IF (WORD(1:1).EQ.'Z') THEN 1348 ILXYZ = 3 1349 ELSE 1350 CALL QUIT('Wrong property in TRANSFORM_HSO, WORD = '//WORD) 1351 END IF 1352C 1353C Transform spin-orbit two-electron integrals, we need sec-occ in linear 1354C response and sec-sec in quadratic response 1355C 1356 IF (HYPCAL.OR.SOMOM.OR.EXMOM) THEN 1357 ITRLSO = 2 1358 ELSE 1359 ITRLSO = 1 1360 END IF 1361 KSYMSO = KSYMOP 1362C defined in RSPSET (920922-ov) 1363 SO2TRA = .TRUE. 1364C 1365C If .OLDTRA has been specified check that MO2SOINT exists 1366C If it does not exist transform AO2SOINT as usual 1367C 1368 IF (OLDTRA) THEN 1369 INQUIRE(FILE='MOHSOINT',EXIST=MOEXIS) 1370 IF (MOEXIS) THEN 1371 SO2TRA = .FALSE. 1372 ELSE 1373 WRITE(LUPRI,'(/3A/A)') ' WARNING: Expected transformed ', 1374 * ' two-electron spin-orbit integrals not', 1375 * ' found.' , 1376 * ' - file MOHSOINT does not exist' 1377 NWARN = NWARN + 1 1378 END IF 1379 END IF 1380 IF (SO2TRA) THEN 1381 IF (IPRHSO.GT.0) THEN 1382 CALL TIMER('START ',TRASTA,TRATIM) 1383 END IF 1384 CALL GPOPEN(LUAHSO,'AO2SOINT','OLD',' ','UNFORMATTED',IDUMMY, 1385 & .FALSE.) 1386 CALL TRAHSO(ITRLSO,CMO,WRK(KWRK1),LWRK1) 1387 CALL GPCLOSE(LUAHSO,'KEEP') 1388 IF (IPRHSO.GT.0) CALL TIMER('TRAHSO',TRASTA,TRATIM) 1389 ELSE 1390 INQUIRE(FILE='MOHSOINT',OPENED=FOPEN) 1391 IF (.NOT.FOPEN) THEN 1392 CALL DAOPEN(LUMHSO,'MOHSOINT') 1393 END IF 1394 END IF 1395 CALL QEXIT('TRANSFORM_HSO') 1396 END 1397 SUBROUTINE HSOFCK(WORD,FC,FV,QA,QB,H2AC,UDV,PV,CMO,LORB,LCON, 1398 & WRK,KWRK1,LWRK1) 1399#include "implicit.h" 1400 CHARACTER*8 WORD 1401 DIMENSION FC(*), FV(*), QA(*), QB(*), H2AC(*), UDV(*), PV(*), 1402 & CMO(*), WRK(*) 1403 LOGICAL LORB, LCON 1404#include "infrsp.h" 1405#include "inforb.h" 1406 CALL QENTER('HSOFCK') 1407 IF (TDHF) THEN 1408 IF (NASHT.GT.0) THEN 1409 CALL DZERO(QA,NORBT*NASHT) 1410 CALL DZERO(QB,NORBT*NASHT) 1411 END IF 1412 CALL HSOFCKAO(WORD,CMO,UDV,FC,FV,WRK,KWRK1,LWRK1) 1413 ELSE 1414 CALL HSOFCKMO(FC,FV,QA,QB,H2AC,UDV,PV,LORB,LCON, 1415 & WRK,KWRK1,LWRK1) 1416 END IF 1417 CALL QEXIT('HSOFCK') 1418 END 1419C 1420C /* Deck hsofckao */ 1421C 1422 SUBROUTINE HSOFCKAO(WORD,CMO,UDV,FC,FV,WRK,KFREE,LFREE) 1423 IMPLICIT NONE 1424 CHARACTER*(*) WORD 1425 REAL*8 CMO(*), UDV(*), FC(*), FV(*), WRK(*) 1426 INTEGER KFREE,LFREE 1427C 1428C Build HSO inactive fockmatrix from AO integrals 1429C If DIRFCK=.TRUE. call TWOINT 1430C else read from AO2SOINT 1431C 1432#include "priunit.h" 1433#include "maxaqn.h" 1434#include "mxcent.h" 1435#include "maxorb.h" 1436#include "dummy.h" 1437#include "infinp.h" 1438#include "wrkrsp.h" 1439#include "inforb.h" 1440#include "aovec.h" 1441#include "infrsp.h" 1442#include "dftcom.h" 1443#include "symmet.h" 1444 REAL*8 D1 1445 PARAMETER (D1=1.0D0) 1446 INTEGER NDMAT,KDMAO,KFMAO,KDMSO,I 1447 INTEGER KDC,KDV,KFC,KFV,KFMO,KFAO 1448 INTEGER KJSTRS , KNPRIM , KNCONT , KIORBS , KJORBS , NPAO , KKORBS 1449 INTEGER I2TYP 1450 INTEGER IPRTWO,IPRNTA,IPRNTB,IPRNTC,IPRNTD 1451 INTEGER IRNTYP , NUMDIS, IFCTYP(2), ISYMDM(2) 1452 REAL*8 TIMEND,TIMSTR,SO1WAL,SO1CPU 1453 REAL*8 SO2WAL,SO2CPU,SOCPU,SOWAL 1454 REAL*8 HFXSAV 1455 LOGICAL TKTIME,RTNTWO 1456 LOGICAL TRIPLET(2) 1457 CALL QENTER('HSOFCKAO') 1458 1459 NDMAT=1 1460 IF (NASHT.GT.0) NDMAT = 2 1461 CALL MEMGET2('REAL','FC',KFC,NDMAT*N2BASX,WRK,KFREE,LFREE) 1462 CALL MEMGET2('REAL','DC',KDC,NDMAT*N2BASX,WRK,KFREE,LFREE) 1463 KFV = KFC + N2BASX 1464 KDV = KDC + N2BASX 1465C 1466C Full density in 1, spin (=active) density in 2, swap for singlet 1467C 1468 CALL GTDMSO(UDV,CMO,WRK(KDC),WRK(KDV),WRK(KFREE)) 1469 IF (NASHT.GT.0) THEN 1470 CALL DAXPY(N2BASX,D1,WRK(KDV),1,WRK(KDC),1) 1471 IF (.NOT.TRPLET) THEN 1472 CALL DSWAP(N2BASX,WRK(KDC),1,WRK(KDV),1) 1473 END IF 1474 END IF 1475 IF (HSODEBUG) THEN 1476 CALL MAOPRI(WRK(KDC),'HSOFCKAO:Input D1') 1477 IF (NASHT.GT.0) THEN 1478 CALL MAOPRI(WRK(KDC+N2BASX),'HSOFCKAO:Input D2') 1479 END IF 1480 END IF 1481 IF (DIRFCK) THEN 1482 IF (WORD(1:1).EQ.'X') IFCTYP(1)=1 1483 IF (WORD(1:1).EQ.'Y') IFCTYP(1)=2 1484 IF (WORD(1:1).EQ.'Z') IFCTYP(1)=3 1485 IF (TRPLET) IFCTYP(1) = -IFCTYP(1) 1486 IFCTYP(2) = - IFCTYP(1) 1487C 1488C Transform density to AO basis 1489C 1490 CALL DSOTAO(WRK(KDC),WRK(KFC),NBAST,0,IPRRSP) 1491 IF (NASHT.GT.0) THEN 1492 CALL DSOTAO(WRK(KDV),WRK(KFV),NBAST,0,IPRRSP) 1493 END IF 1494 CALL DCOPY(NDMAT*N2BASX,WRK(KFC),1,WRK(KDC),1) 1495C 1496C Setup for TWOINT 1497C 1498 NPAO = MXSHEL*MXAOVC 1499 CALL MEMGET('INTE',KJSTRS,NPAO*2,WRK,KFREE,LFREE) 1500 CALL MEMGET('INTE',KNPRIM,NPAO*2,WRK,KFREE,LFREE) 1501 CALL MEMGET('INTE',KNCONT,NPAO*2,WRK,KFREE,LFREE) 1502 CALL MEMGET('INTE',KIORBS,NPAO ,WRK,KFREE,LFREE) 1503 CALL MEMGET('INTE',KJORBS,NPAO ,WRK,KFREE,LFREE) 1504 CALL MEMGET('INTE',KKORBS,NPAO ,WRK,KFREE,LFREE) 1505 CALL PAOVEC(WRK(KJSTRS),WRK(KNPRIM),WRK(KNCONT), 1506 & WRK(KIORBS),WRK(KJORBS),WRK(KKORBS),0, 1507 & .FALSE.,IPRRSP) 1508 CALL MEMREL('HERINT.PAOVEC',WRK,1,KJORBS,KFREE,LFREE) 1509 CALL TIMER('START ',TIMSTR,TIMEND) 1510 CALL GETTIM(SO1CPU,SO1WAL) 1511 I2TYP = 0 1512 IRNTYP = -20 1513 IF (HSODEBUG) THEN 1514 IPRTWO = 3 1515 ELSE 1516 IPRTWO = 0 1517 END IF 1518 TKTIME = .FALSE. 1519 CALL DZERO(WRK(KFC),NDMAT*N2BASX) 1520C Always full exchange for spin-orbit integrals: 1521 HFXSAV=HFXFAC 1522 HFXFAC=D1 1523 CALL TWOINT(WRK(KFREE),LFREE,WRK(KFREE), 1524 & WRK(KFC),WRK(KDC),NDMAT, 1525 & IDUMMY, IFCTYP, 1526 & DUMMY,IDUMMY,NUMDIS,1,IRNTYP,0,0,.TRUE.,.TRUE., 1527 & .FALSE.,TKTIME,IPRTWO,IPRNTA,IPRNTB,IPRNTC, 1528 & IPRNTD,RTNTWO,IDUMMY,I2TYP,WRK(KJSTRS), 1529 & WRK(KNPRIM),WRK(KNCONT),WRK(KIORBS), 1530 & IDUMMY,IDUMMY,DUMMY,DUMMY,DUMMY,DUMMY, 1531 & .FALSE.,.false.) 1532 HFXFAC=HFXSAV 1533 CALL MEMREL('HSOFCKAO.TWOINT',WRK,KFC,KJSTRS,KFREE,LFREE) 1534 IF (HSODEBUG) THEN 1535 CALL MAOPRI(WRK(KFC),':TWOINT:Calculated F1') 1536 IF (NASHT.GT.0) THEN 1537 CALL MAOPRI(WRK(KFC+N2BASX),':TWOINT:Calculated F2') 1538 END IF 1539 END IF 1540C 1541C Transform to SO 1542C 1543 ISYMDM(1)=0 1544 ISYMDM(2)=0 1545 CALL SKLFCK(WRK(KFC),DUMMY,WRK(KFREE),LFREE, 1546 & IPRTWO,.FALSE., 1547 & .FALSE.,.FALSE.,.FALSE.,.TRUE.,IDUMMY,.TRUE.,NDMAT, 1548 & ISYMDM,IFCTYP,IDUMMY,.TRUE.) 1549C 1550 IF (HSODEBUG) THEN 1551 CALL MAOPRI(WRK(KFC),':SKLFCK:Calculated F1') 1552 IF (NASHT.GT.0) THEN 1553 CALL MAOPRI(WRK(KFC+N2BASX),':SKLFCK:Calculated F2') 1554 END IF 1555 END IF 1556 CALL GETTIM(SO2CPU,SO2WAL) 1557 SOCPU=SO2CPU-SO1CPU 1558 SOWAL=SO2WAL-SO1WAL 1559 CALL TIMER('TWOINT',TIMSTR,TIMEND) 1560 CALL FLSHFO(LUPRI) 1561 WRITE(LUPRI,'(2(/A),2(/A,F5.0,A))') 1562 & ' Two-electron spin-orbit integrals', 1563 & ' =================================', 1564 & ' Spin-orbit 2-electron CPU time ',SOCPU,' seconds', 1565 & ' Spin-orbit 2-electron wall time ',SOWAL,' seconds' 1566 CALL MEMCHK('HSOFCKAO.SKLFCK',WRK,KFC) 1567 ELSE 1568C 1569C Read AO integrals from disk (AO2SOINT) 1570C 1571 TRIPLET(1) = TRPLET 1572 TRIPLET(2) = .NOT.TRPLET 1573 CALL DZERO(WRK(KFC),NDMAT*N2BASX) 1574 CALL GET_FSO_AO(WORD,TRIPLET,WRK(KFC),WRK(KDC),NDMAT) 1575 END IF 1576 IF (HSODEBUG) THEN 1577 CALL MAOPRI(WRK(KFC),'HSOFCKAO:Calculated F1') 1578 IF (NASHT.GT.0) THEN 1579 CALL MAOPRI(WRK(KFC+N2BASX),'HSOFCKAO:Calculated F2') 1580 END IF 1581 END IF 1582C 1583C Reorder matrices to match gradient routine 1584C (fa+fb,fa-fb)/2 = (fc,fo-fc) -> (fo,fc-fo) 1585C 1586 IF (NASHT.GT.0) THEN 1587 CALL DAXPY(N2BASX,D1,WRK(KFV),1,WRK(KFC),1) !now (fo,fo-fc) 1588C CALL DSWAP(N2BASX,WRK(KFC),1,WRK(KFV),1) 1589 CALL DSCAL(N2BASX,-D1,WRK(KFV),1) !now (fo,fc-fo) 1590 END IF 1591C 1592 CALL DZERO(FC,N2ORBX) 1593 CALL AOTOMO(WRK(KFC),FC,CMO,KSYMOP,WRK(KFREE),LFREE) 1594 IF (NASHT.GT.0) THEN 1595 CALL DZERO(FV,N2ORBX) 1596 CALL AOTOMO(WRK(KFV),FV,CMO,KSYMOP,WRK(KFREE),LFREE) 1597 END IF 1598 IF (HSODEBUG) THEN 1599 CALL MAOPRI(WRK(KFC),'HSOFCKAO:FC(AO)') 1600 IF (NASHT.GT.0) CALL MAOPRI(WRK(KFV),'HSOFCKAO:FC-FO(AO)') 1601 CALL MMOPRI(FC,'HSOFCKAO:FO(MO)') 1602 IF (NASHT.GT.0) CALL MMOPRI(FV,'HSOFCKAO:FC-FO(MO)') 1603 END IF 1604 CALL MEMREL('HSOFCKAO',WRK,KFC,KFC,KFREE,LFREE) 1605 CALL QEXIT('HSOFCKAO') 1606 END 1607C /* Deck hsofckmo */ 1608 SUBROUTINE HSOFCKMO(FC,FV,QA,QB,H2AC,UDV,PV,LORB,LCON, 1609 & WRK,KWRK1,LWRK1) 1610#include "implicit.h" 1611 DIMENSION FC(*), FV(*), QA(*), QB(*), H2AC(*), UDV(*), PV(*), 1612 & WRK(*) 1613 LOGICAL LORB, LCON 1614#include "infhso.h" 1615#include "maxorb.h" 1616#include "maxash.h" 1617#include "infind.h" 1618#include "inforb.h" 1619#include "priunit.h" 1620#include "wrkrsp.h" 1621 INTEGER NEEDSO(-4:6) 1622 CALL QENTER('HSOFCKMO') 1623 CALL DZERO(FC,N2ORBX) 1624 IF (LORB) THEN 1625 CALL DZERO(FV,N2ORBX) 1626 CALL DZERO(QA,NORBT*NASHT) 1627 CALL DZERO(QB,NORBT*NASHT) 1628 END IF 1629 IF (LCON) THEN 1630 CALL DZERO(H2AC,NNASHX*NNASHX) 1631 END IF 1632 IF (.NOT.DOSO2) THEN 1633 CALL QEXIT('HSOFCK') 1634 RETURN 1635 END IF 1636 CALL MEMGET2('REAL','H2' ,KH2 ,N2ORBX,WRK,KWRK1,LWRK1) 1637 CALL MEMGET2('REAL','PV12',KPV12,N2ASHX,WRK,KWRK1,LWRK1) 1638 CALL MEMGET2('REAL','PV21',KPV21,N2ASHX,WRK,KWRK1,LWRK1) 1639C 1640C Read distributions 1641C 1642 NEEDSO(-4:6) = 0 1643 NEEDSO( 1:5) = 1 1644 IDIST = 0 1645 KFREE = 1 1646 LFREE = LWRK1 1647 90 CALL NXTHSO(IC,ID,WRK(KH2),NEEDSO,WRK(KWRK1),KFREE,LFREE,IDIST) 1648 IF (IDIST.LT.0) GOTO 95 1649 IF (IC.EQ.ID) GOTO 90 1650 IF (IPRHSO.GT.20) THEN 1651 WRITE(LUPRI,'(//A,2I5)')' Integral distribution ',IC,ID 1652 CALL OUTPUT(WRK(KH2),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1653 END IF 1654C 1655C Construct inactive and active spin-orbit Fock matrices 1656C 1657 CALL FSOMU(IC,ID,WRK(KH2),FC,FV,UDV,LORB, 1658 * WRK(KWRK1-1+KFREE),LFREE) 1659C CALL FSOMU(ICI1,IDI1,H2,FCSO,FVSO,UDV,LORB,WRK,LWRK) 1660C 1661C 1662C Construct spin-orbit Q matrices 1663C 1664 IF (NASHT.GT.0 .AND. LORB) 1665 * CALL QSOMU(IC,ID,QA,QB,WRK(KH2), 1666 * PV,WRK(KPV12),WRK(KPV21), 1667 * WRK(KWRK1-1+KFREE),LFREE) 1668C 1669C Add active-active distributions to H2AC(uv,xy) 1670C CALL ADH2AC(H2ACXY,H2XY,IUVSYM) 1671C 1672 IF (IOBTYP(IC).EQ.JTACT .AND. IOBTYP(ID).EQ.JTACT .AND. 1673 * LCON) THEN 1674 ICSYM = ISMO(IC) 1675 IDSYM = ISMO(ID) 1676 ICDSYM = MULD2H(ICSYM,IDSYM) 1677 KCDSYM = MULD2H(KSYMOP,ICDSYM) 1678 NCW = ICH(IC) 1679 NDW = ICH(ID) 1680 IF (NCW.GT.NDW) THEN 1681 NCDW = IROW(NCW) + NDW 1682 KH2XY = 1 + (NCDW-1)*NNASHX 1683 CALL ADH2AC(H2AC(KH2XY),WRK(KH2),KCDSYM) 1684 ELSE 1685 NCDW = IROW(NDW) + NCW 1686 KH2XY = 1 + (NCDW-1)*NNASHX 1687 CALL DAXPY(N2ORBX,-1.0D0,WRK(KH2),1,WRK(KWRK1-1+KFREE),1) 1688 CALL ADH2AC(H2AC(KH2XY),WRK(KWRK1-1+KFREE),KCDSYM) 1689 END IF 1690 END IF 1691 GOTO 90 1692 95 CONTINUE 1693C 1694 IF (IPRHSO.GT.10) THEN 1695 WRITE(LUPRI,'(//A)') ' Final Inactive Fock matrix' 1696 CALL OUTPUT(FC,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1697C 1698 IF (LORB .AND. NASHT.GT.0) THEN 1699 WRITE(LUPRI,'(//A)') ' Final Active Fock matrix' 1700 CALL OUTPUT(FV,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1701C 1702 WRITE(LUPRI,'(//A)') ' Spin-orbit QA matrix' 1703 CALL OUTPUT(QA,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1704 WRITE(LUPRI,'(//A)') ' Spin-orbit QB matrix' 1705 CALL OUTPUT(QB,1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1706 END IF 1707C 1708 IF (LCON) THEN 1709 DO 71 K=1,NASHT 1710 DO 72 L=1,K 1711 KL = IROW(K)+L 1712 KLOFF = (KL-1)*NNASHX 1713 WRITE(LUPRI,'(/A,2I3,A)') 1714 * ' Active integral distribution (**', K,L,')' 1715 CALL OUTPAK(H2AC(1+KLOFF),NASHT,1,LUPRI) 1716 72 CONTINUE 1717 71 CONTINUE 1718 END IF 1719 END IF 1720 CALL MEMREL('HSOFCKMO',WRK,KH2,KH2,KWRK1,LWRK1) 1721 CALL QEXIT('HSOFCKMO') 1722 END 1723C /* Deck hsoorb */ 1724 SUBROUTINE HSOORB(ONEIND,NSIM,FC,FCX,FVX,QAX,QBX,UDV, 1725 * EVECS,TRPDEN) 1726C 1727C WRITTEN 14-FEB 1986 1728C adapted 5-Apr 2001 1729C 1730C PURPOSE: 1731C 1)DISTRIBUTE INACTIVE FCX AND ACTIVE FVX FOCK MATRICES AND QAX AND 1732C QBX MATRICES INTO ORBITAL PART OF LINEAR TRANSFORMED VECTOR 1733C 1734C ( [ E(K,L) , H ] ) K<L 1735C E[2]*N = - ( ) 1736C ( [ E(L,K) , H ] ) K>L 1737C X MAY REFER TO EITHER ONE INDEX TRANSFORMED MATRICES (N IS A 1738C ORBITAL TRIAL VECTOR ) OR TRANSITION DENSITY MATRICES (N IS A 1739C CONFIGURATION TRIAL VECTOR). FOR CONFIGURATION TRIAL VECTORS 1740C FCX IS IDENTICALLY ZERO. FURTHER OVERLAP IS ZERO BECAUSE TRIAL 1741C VECTORS ARE CHOSEN ORTHOGONAL TO REFERENCE STATE. 1742C 1743C 2)CREATE LINEAR TRANSFORMED VECTOR S[2]*N FOR N EQUAL TO EITHER A 1744C ORBITAL OR A CONFIGURATION TRIAL VECTOR 1745C 1746C ONEIND = .TRUE. FOR A ORBITAL TRIAL VECTOR 1747C 1748C ****************************** 1749C 1750C [ E(P,Q) , H ] = 1751C 1752C ACTIVE-INACTIVE (P,Q) = (T,I) 1753C L,K 1754C FCX(I,X)*DV(T,X)-2*FCX(I,T)-2FVX(I,T)+QBX(I,T) 1755C K L K,L K,L K,L 1756C 1757C INACTIVE-ACTIVE (I,U) 1758C K,L 1759C -FCX(X,I)*DV(X,U)+2*FCX(U,I)+2*FVX(U,I)-QAX(I,U) 1760C K L L,K L,K K,L 1761C 1762C ACTIVE-ACTIVE (T,U) 1763C K,L 1764C L,K 1765C FCX(U,X)*DV(T,X)-FCX(X,T)*DV(X,U)+QBX(U,T)-QAX(T,U) 1766C L K K L L,K K,L 1767C K L L K K,L L,K 1768C 1769C ACTIVE-SECUNDARY (T,A) 1770C K,L 1771C FCX(A,X)*DV(T,X)+QBX(A,T) 1772C L K L,K 1773C 1774C SECUNDARY-ACTIVE (A,U) 1775C L,K 1776C -FCX(X,A)*DV(X,U)-QAX(A,U) 1777C L K L,K 1778C 1779C INACTIVE-SECUNDARY (I,A) 1780C K,L 1781C 2*FCX(A,I)+2FVX(A,I) 1782C L,K L,K 1783C 1784C SECUNDARY-INACTIVE (A,J) 1785C L,K 1786C -2*FCX(J,A)-2*FVX(J,A) 1787C K,L K,L 1788C 1789C 1790#include "implicit.h" 1791#include "priunit.h" 1792C 1793 LOGICAL ONEIND, TRPDEN 1794C 1795 DIMENSION FC(*),FCX(NORBT,NORBT,*),FVX(NORBT,NORBT,*) 1796 DIMENSION QAX(NORBT,NASHDI,*),QBX(NORBT,NASHDI,*) 1797 DIMENSION UDV(NASHDI,NASHDI,*), EVECS(KZYVAR,*) 1798C 1799C INFDIM : NASHDI 1800C 1801#include "maxorb.h" 1802#include "maxash.h" 1803#include "infvar.h" 1804#include "inforb.h" 1805#include "infind.h" 1806#include "infdim.h" 1807#include "infpri.h" 1808#include "infrsp.h" 1809#include "wrkrsp.h" 1810C 1811C -- local constants 1812C 1813 PARAMETER ( D0 = 0.0D0 , D2 =2.0D0 , DM1 = -1.0D0 ) 1814C 1815C 1816 KYCONF = KZCONF + KZVAR 1817 IF (TRPDEN) THEN 1818 DC = D0 1819 ELSE 1820 DC = D2 1821 END IF 1822C 1823C DISTRIBUTE FOCK MATRICES IN EVECS 1824C 1825 KSYM1 = 0 1826 DO 1300 IG = 1,KZWOPT 1827 K = JWOP(1,IG) 1828 L = JWOP(2,IG) 1829 KSYM = ISMO(K) 1830 LSYM = ISMO(L) 1831 IF( KSYM.NE.KSYM1 ) THEN 1832 KSYM1 = KSYM 1833 NORBK = NORB(KSYM) 1834 IORBK = IORB(KSYM) 1835 NASHK = NASH(KSYM) 1836 NISHK = NISH(KSYM) 1837 IASHK = IASH(KSYM) 1838 IIORBK= IIORB(KSYM) 1839 IORBL = IORB(LSYM) 1840 NASHL = NASH(LSYM) 1841 NISHL = NISH(LSYM) 1842 NORBL = NORB(LSYM) 1843 IASHL = IASH(LSYM) 1844 IIORBL= IIORB(LSYM) 1845 END IF 1846 NK = K - IORBK 1847 NL = L - IORBL 1848 ITYPK = IOBTYP(K) 1849 ITYPL = IOBTYP(L) 1850 IF ( ITYPK.EQ.JTINAC )THEN 1851 DO 2000 ISIM = 1 ,NSIM 1852 IF (ONEIND) THEN 1853 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1854 * + DC * FCX(L,K,ISIM) 1855 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1856 * - DC * FCX(K,L,ISIM) 1857 IF ( NASHT . GT . 0 ) THEN 1858 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1859 * + D2 * FVX(L,K,ISIM) 1860 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1861 * - D2 * FVX(K,L,ISIM) 1862 END IF 1863 ELSE 1864 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1865 * + D2 * FVX(L,K,ISIM) 1866 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1867 * - D2 * FVX(K,L,ISIM) 1868 END IF 1869 IF ( ITYPL.EQ.JTACT ) THEN 1870 NWL = ISW(L) - NISHT 1871 IF (ONEIND) THEN 1872 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1873 * - DDOT(NASHL,FCX(IORBL+NISHL+1,K,ISIM),1, 1874 * UDV(IASHL+1,NWL,1),1) 1875 DO 730 IX = 1,NASHL 1876 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1877 * + FCX(K,IORBL+NISHL+IX,ISIM)* 1878 * UDV(NWL,IASHL+IX,1) 1879 730 CONTINUE 1880 ELSE 1881 TEMP1 = D0 1882 TEMP2 = D0 1883 NX = NISHK 1884 DO 825 NWX = IASHK+1,IASHK+NASHK 1885 NX = NX + 1 1886 IF (NX.LE.NK) THEN 1887 FCXK = FC(IIORBK+IROW(NK)+NX) 1888 ELSE 1889 FCXK = FC(IIORBK+IROW(NX)+NK) 1890 END IF 1891 TEMP1 = TEMP1 + FCXK * UDV(NWX,NWL,ISIM) 1892 TEMP2 = TEMP2 + FCXK * UDV(NWL,NWX,ISIM) 1893 825 CONTINUE 1894 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1895 * - TEMP1 1896 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1897 * + TEMP2 1898 END IF 1899 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1900 * - QAX(K,NWL,ISIM) 1901 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1902 * + QBX(K,NWL,ISIM) 1903 END IF 1904 2000 CONTINUE 1905 ELSE 1906 IF (ITYPL.EQ.JTACT) THEN 1907 NWL = ISW(L) - NISHT 1908 NWK = ISW(K) - NISHT 1909 DO 2020 ISIM=1,NSIM 1910 IF (ONEIND) THEN 1911 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1912 * -DDOT(NASHL,FCX(IORBL+NISHL+1,K,ISIM),1, 1913 * UDV(IASHL+1,NWL,1),1) 1914 DO 740 IX = 1,NASHL 1915 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1916 * +FCX(K,IORBL+NISHL+IX,ISIM)* 1917 * UDV(NWL,IASHL+IX,1) 1918 740 CONTINUE 1919 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1920 * -DDOT(NASHK,FCX(IORBK+NISHK+1,L,ISIM),1, 1921 * UDV(IASHK+1,NWK,1),1) 1922 DO 750 IX = 1,NASHK 1923 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1924 * +FCX(L,IORBK+NISHK+IX,ISIM)* 1925 * UDV(NWK,IASHK+IX,1) 1926 750 CONTINUE 1927 ELSE 1928 TEMP1 = D0 1929 TEMP2 = D0 1930 NX = NISHL 1931 DO 835 NWX = IASHL+1,IASHL+NASHL 1932 NX = NX + 1 1933 IF (NX.LE.NL) THEN 1934 FCXL = FC(IIORBL+IROW(NL)+NX) 1935 ELSE 1936 FCXL = FC(IIORBL+IROW(NX)+NL) 1937 END IF 1938 TEMP1 = TEMP1 + FCXL * UDV(NWX,NWK,ISIM) 1939 TEMP2 = TEMP2 + FCXL * UDV(NWK,NWX,ISIM) 1940 835 CONTINUE 1941 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1942 * + TEMP2 1943 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1944 * - TEMP1 1945 TEMP1 = D0 1946 TEMP2 = D0 1947 NX = NISHK 1948 DO 845 NWX = IASHK+1,IASHK+NASHK 1949 NX = NX + 1 1950 IF (NX.LE.NK) THEN 1951 FCXK = FC(IIORBK+IROW(NK)+NX) 1952 ELSE 1953 FCXK = FC(IIORBK+IROW(NX)+NK) 1954 END IF 1955 TEMP1 = TEMP1 + FCXK * UDV(NWX,NWL,ISIM) 1956 TEMP2 = TEMP2 + FCXK * UDV(NWL,NWX,ISIM) 1957 845 CONTINUE 1958 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1959 * - TEMP1 1960 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1961 * + TEMP2 1962 END IF 1963 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1964 * +QBX(L,NWK,ISIM) -QAX(K,NWL,ISIM) 1965 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1966 * +QBX(K,NWL,ISIM) - QAX(L,NWK,ISIM) 1967 2020 CONTINUE 1968 ELSE 1969 NWK = ISW(K) - NISHT 1970 DO 2030 ISIM=1,NSIM 1971 IF (ONEIND) THEN 1972 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1973 * -DDOT(NASHK,FCX(IORBK+NISHK+1,L,ISIM),1, 1974 * UDV(IASHK+1,NWK,1),1) 1975 DO 760 IX = 1,NASHK 1976 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1977 * +FCX(L,IORBK+NISHK+IX,ISIM)* 1978 * UDV(NWK,IASHK+IX,1) 1979 760 CONTINUE 1980 ELSE 1981 TEMP1 = D0 1982 TEMP2 = D0 1983 NX = NISHL 1984 DO 860 NWX = IASHL+1,IASHL+NASHL 1985 NX = NX + 1 1986 IF (NX.LE.NL) THEN 1987 FCXL = FC(IIORBL+IROW(NL)+NX) 1988 ELSE 1989 FCXL = FC(IIORBL+IROW(NX)+NL) 1990 END IF 1991 TEMP2 = TEMP2 + FCXL * UDV(NWX,NWK,ISIM) 1992 TEMP1 = TEMP1 + FCXL * UDV(NWK,NWX,ISIM) 1993 860 CONTINUE 1994 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 1995 * + TEMP1 1996 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 1997 * - TEMP2 1998 END IF 1999 EVECS(KYCONF+IG,ISIM) = EVECS(KYCONF+IG,ISIM) 2000 * -QAX(L,NWK,ISIM) 2001 EVECS(KZCONF+IG,ISIM) = EVECS(KZCONF+IG,ISIM) 2002 * +QBX(L,NWK,ISIM) 2003 2030 CONTINUE 2004 ENDIF 2005 ENDIF 2006 1300 CONTINUE 2007C 2008C CHANGE SIGN ON ORBITAL PART OF LINEAR TRANSFORMATION E[2]*N 2009C 2010 DO 3000 ISIM = 1,NSIM 2011 CALL DSCAL(KZWOPT,DM1,EVECS(KZCONF+1,ISIM),1) 2012 CALL DSCAL(KZWOPT,DM1,EVECS(KYCONF+1,ISIM),1) 2013 3000 CONTINUE 2014C 2015 IF (IPRRSP.GT.110) THEN 2016 WRITE(LUPRI,*)' HSOORB: LINEAR TRANSFORMATION WITH E(2)' 2017 CALL OUTPUT(EVECS,1,KZYVAR,1,NSIM,KZYVAR,NSIM,1,LUPRI) 2018 END IF 2019C 2020C END OF HSOORB 2021C 2022 RETURN 2023 END 2024C /* Deck hsoal2 */ 2025 SUBROUTINE HSOAL2 (WORD,GP,CMO,UDV,PV,XINDX,MJWOP,WRK,LWRK) 2026C 2027C Alternative routine for constructing spin-orbit gradient 2028C by way of X2HSO Case 2 2029C 2030#include "implicit.h" 2031 DIMENSION GP(*), CMO(*),UDV(*),PV(*),XINDX(*),WRK(*) 2032 CHARACTER*8 WORD 2033C 2034 PARAMETER (D1=1.0D0, DM1=-1.0D0, DH=0.5D0) 2035 LOGICAL DO1,DO2,LORB,LCON,LREFST,TRPSAV 2036 CHARACTER*8 LABEL 2037C 2038C Used from common blocks 2039C 2040C INFINP: FLAG() 2041C INFORB: NORBT... 2042C INFVAR: NCONF 2043C WRKRSP: KZYVAR...,KSYMOP 2044C INFRSP: NCREF,IREFSY,DIROIT,SOPPA 2045C INFHSO: DOSO1 2046C INFTAP: LUINDX, LUMHSO 2047C 2048#include "priunit.h" 2049#include "maxorb.h" 2050#include "infinp.h" 2051#include "inforb.h" 2052#include "wrkrsp.h" 2053#include "infrsp.h" 2054#include "infvar.h" 2055#include "qrinf.h" 2056 DIMENSION MJWOP(2,MAXWOP,8) 2057#include "infhso.h" 2058#include "trhso.h" 2059#include "inftap.h" 2060C 2061 CALL QENTER('HSOAL2') 2062 CALL HEADER('HSOAL2: TEST OF SPIN-ORBIT GRADIENT',1) 2063 IF (SOPPA) CALL QUIT('HSOAL2: SOPPA not implemented yet!') 2064 IF (.NOT.TDHF .AND. .NOT.FLAG(27)) 2065 & CALL QUIT('HSOAL2 is only implemented for .DETERMINANTS') 2066 IF (.NOT.TDHF .AND. NCONF.NE.NCREF) 2067 & CALL QUIT('HSOAL2 inconsistency : NCONF.ne.NCREF') 2068C 2069C Initialise MZYVAR... for symmetries KSYMOP and 1 (from INFVAR) 2070C 2071C Symmetry KSYMOP: 2072C 2073 CALL SETZY(MJWOP) 2074C 2075C Symmetry 1 2076C 2077 LUINDX = -1 2078 CALL GPOPEN(LUINDX,'LUINDF','UNKNOWN',' ','UNFORMATTED',IDUMMY, 2079 & .FALSE.) 2080 REWIND LUINDX 2081 CALL MOLLAB('EXOPSYM1',LUINDX,LUERR) 2082 READ (LUINDX) IWOPT,IWOP,IWOPH 2083 CALL GPCLOSE(LUINDX,'KEEP') 2084 IVAR = IWOPT + NCONF 2085 MZVAR(1) = IVAR 2086 MZCONF(1) = NCONF 2087 MZWOPT(1) = IWOPT 2088 MZYVAR(1) = 2*IVAR 2089 MZYCON(1) = 2*NCONF 2090 MZYWOP(1) = 2*IWOPT 2091 WRITE(LUPRI,'(/A)') ' Number of variables symmetry 1' 2092 WRITE(LUPRI,'( A,I5)') ' Rotations: ',MZWOPT(1) 2093 WRITE(LUPRI,'( A,I5)') ' Configurations:',MZCONF(1) 2094 WRITE(LUPRI,'(/A,I2)') ' Number of variables symmetry',KSYMOP 2095 WRITE(LUPRI,'( A,I5)') ' Rotations: ',MZWOPT(KSYMOP) 2096 WRITE(LUPRI,'( A,I5)') ' Configurations:',MZCONF(KSYMOP) 2097 WRITE(LUPRI,'(/A)') ' Rotational indices K L' 2098 DO 30 I=1,MZWOPT(KSYMOP) 209930 WRITE(LUPRI,'(3I5)') I,MJWOP(1,I,KSYMOP),MJWOP(2,I,KSYMOP) 2100C 2101C Allocate work space 2102C 2103 KFREE = 1 2104 LFREE = LWRK 2105 CALL MEMGET('REAL',KVEC,2*IVAR,WRK,KFREE,LFREE) 2106 CALL MEMGET('REAL',KFI,N2ORBX,WRK,KFREE,LFREE) 2107C KFI = KVEC + 2*IVAR 2108C KFREE = KFI + N2ORBX 2109C LFREE = LWRK - KFREE + 1 2110C IF (LFREE.LT.0) CALL ERRWRK('HSOAL2',KFREE-1,LWRK) 2111C 2112C Convention: X1SPNORB: one-electron parts only 2113C X2SPNORB: two-electron parts only 2114C X SPNORB: both (default) 2115C overridden by SO1ONLY or SO2ONLY in HSOINP 2116C 2117 IF (WORD(2:2).EQ.'1' .OR. WORD(2:2).EQ.' '.AND.DOSO1) THEN 2118 LABEL = WORD 2119 LABEL(2:2) = '1' 2120 KSYMP = -1 2121 CALL PRPGET(LABEL,CMO,WRK(KFI),KSYMP,ANTSYM, 2122 & WRK(KFREE),LFREE,IPRHSO) 2123 IF (KSYMP.NE.KSYMOP) THEN 2124 WRITE (LUPRI,'(/A/2A/A,2I5/A,F10.2)') 2125 & 'FATAL ERROR: KSYMOP .ne. KSYMP from PRPGET', 2126 & ' Property label : ',LABEL, 2127 & ' KSYMOP and KSYMP:',KSYMOP,KSYMP, 2128 & ' ANTSYM :',ANTSYM 2129 CALL QUIT('KSYMOP .ne. KSYMP from PRPGET') 2130 END IF 2131 ELSE 2132 CALL DZERO(WRK(KFI),N2ORBX) 2133 END IF 2134C 2135C 2136C Set up the vector to be N_o = ( 0 ) , N_c = ( - Ref ) 2137C ( 0 ) ( Ref ) 2138C 2139 CALL DZERO(WRK(KVEC),2*IVAR) 2140 CALL GETREF(WRK(KVEC + IVAR),NCONF) 2141 CALL DAXPY(NCONF,DM1,WRK(KVEC + IVAR),1,WRK(KVEC),1) 2142 WRITE(LUPRI,'(/A)') 'REFERENCE CI VECTOR' 2143 CALL OUTPUT(WRK(KVEC + IVAR),1,NCONF,1,1,NCONF,1,1,LUPRI) 2144 WRITE(LUPRI,'(/A)') 'CALLING HSO2CR WITH VECTOR' 2145 CALL RSPPRC(WRK(KVEC),NCONF,IVAR,LUPRI) 2146 CALL RSPPRW(WRK(KVEC+NCONF),MJWOP,IWOPT,1,IVAR,LUPRI) 2147 ISYMV = 1 2148 ISPINV = 0 2149 IKLVL = 0 2150 DIROIT = .TRUE. 2151C 2152C Gradient setup 2153C 2154 IGRSYM = KSYMOP 2155 IF (TRPLET) THEN 2156 IGRSPI = 1 2157 ELSE 2158 IGRSPI = 0 2159 END IF 2160C 2161C Operator setup 2162C 2163 IOPSYM = KSYMOP 2164 IOPSPI = 1 2165C 2166C Orbital gradient 2167C 2168 LORB = KZWOPT.GT.0 2169 LPV = 2*N2ASHX*N2ASHX 2170 CALL MEMGET('REAL',KD,N2ASHX,WRK,KFREE,LFREE) 2171 CALL MEMGET('REAL',KP,LPV,WRK,KFREE,LFREE) 2172 KP1=KP 2173 KP2=KP+N2ASHX*N2ASHX 2174 IF (TRPLET) THEN 2175 CALL DCOPY(N2ASHX,UDV,1,WRK(KD),1) 2176 CALL DCOPY(LPV,PV,1,WRK(KP),1) 2177 CALL IPSET(0,0,1,1) 2178 ELSE 2179 IF (TDHF) THEN 2180 IF (NASHT.EQ.1) THEN 2181 WRK(KD)=1.0D0 2182 END IF 2183 ELSE 2184 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 2185 CALL GETREF(WRK(KCREF),MZCONF(1)) 2186 CALL RSPDM(IREFSY,IREFSY,MZCONF(1),MZCONF(1), 2187 & WRK(KCREF),WRK(KCREF), 2188 & WRK(KD),WRK(KP),0,1,.FALSE.,.FALSE.,XINDX, 2189 & WRK,KFREE,LFREE) 2190 CALL MEMREL('HSOCTL<-RSPDM',WRK,KD,KCREF,KFREE,LFREE) 2191 CALL MTRSP(N2ASHX,N2ASHX,WRK(KP1),N2ASHX,WRK(KP2),N2ASHX) 2192 END IF 2193 CALL IPSET(0,1,1,0) 2194 END IF 2195 OVLAP = D1 2196 ISYMDN = 1 2197C 2198C Configurational gradient 2199C 2200 IF ( IGRSYM.EQ.1) THEN 2201 LCON = KZCONF.GT.1 2202 ELSE 2203 LCON = KZCONF.GT.0 2204 END IF 2205 LREFST = .FALSE. 2206C 2207C Create the gradient 2208C 2209 TRPSAV=TRPLET 2210 CALL HSO2CR(IGRSYM,IGRSPI,GP,WRK(KVEC), 2211 * 2*IVAR,NCONF,ISYMV,ISYMDN, 2212 * XINDX,OVLAP,WRK(KD),WRK(KP),WRK(KFI), 2213 * WRK(KFREE),LFREE, 2214 * KZYVAR,LCON,LORB,LREFST,LUMHSO,KSYMOP, 2215 * IOPSPI,ISPINV,IKLVL,DUM,IDUM,DUM,IDUM,MJWOP) 2216 CALL HEADER('SPIN-ORBIT GRADIENT FROM HSOAL2',3) 2217 TRPLET=TRPSAV 2218C 2219 IF (IPRRSP.GT.100) THEN 2220 WRITE(LUPRI,'(/A)') ' Orbital part' 2221 CALL OUTPUT(GP(1+KZCONF),1,KZWOPT,1,2,KZVAR,2,1,LUPRI) 2222 WRITE(LUPRI,'(/A)') ' Configuration part' 2223 CALL OUTPUT(GP,1,KZCONF,1,2,KZVAR,2,1,LUPRI) 2224 ELSE 2225 CALL RSPPRW(GP(1+KZCONF),MJWOP,KZWOPT,KSYMOP,KZVAR,LUPRI) 2226 CALL RSPPRC(GP,KZCONF,KZVAR,LUPRI) 2227 END IF 2228 IF (KSYMOP.EQ.1 .AND. LCON) THEN 2229 WRITE(LUPRI,'(/A)') 'Project out reference component' 2230 KCREF = KVEC + IVAR 2231 OVLAP = DDOT(KZCONF,WRK(KCREF),1,GP,1) 2232 CALL DAXPY(KZCONF,-OVLAP,WRK(KCREF),1,GP,1) 2233 OVLAP = DDOT(KZCONF,WRK(KCREF),1,GP(1+KZVAR),1) 2234 CALL DAXPY(KZCONF,-OVLAP,WRK(KCREF),1,GP(1+KZVAR),1) 2235 CALL RSPPRW(GP(1+KZCONF),MJWOP,KZWOPT,KSYMOP,KZVAR,LUPRI) 2236 CALL RSPPRC(GP,KZCONF,KZVAR,LUPRI) 2237 END IF 2238 CALL MEMREL('HSOAL2',WRK,1,1,KFREE,LFREE) 2239 CALL QEXIT('HSOAL2') 2240 RETURN 2241 END 2242C 2243C /* Deck get_FSO_AO */ 2244C 2245 SUBROUTINE GET_FSO_AO(WORD,TRIPLET,F,D,ND) 2246C 2247C Generalize GETFMAT to handle non-symmetric densities (quadratic 2248C response spin-orbit) 2249C 2250#include "implicit.h" 2251 CHARACTER*8 WORD 2252 INTEGER ND 2253 LOGICAL TRIPLET(ND) 2254#include "inforb.h" 2255 DIMENSION F(NBAST,NBAST,ND), D(NBAST,NBAST,ND) 2256 2257#include "infhso.h" 2258#include "priunit.h" 2259#include "iratdef.h" 2260#include "eribuf.h" 2261C 2262C LOCAL 2263C 2264 PARAMETER (LBUF_alloc = 600) 2265 REAL*8 BUF(2*LBUF_alloc+1) 2266 INTEGER IINDX4(4,LBUF_alloc), LENGTH 2267 INTEGER IB, I,J,K,L,ICOOR, LUHSOAO, COMPONENT,ID 2268 REAL*8 LEFT(2), RIGHT(2) 2269 REAL*8 GIJKL,XIJKL,DIJ,DKL, DP5, D1P5, D1, D2 2270 PARAMETER(DP5=0.5D0, D1P5=1.5D0, D1=1.0D0, D2=2.0D0) 2271C 2272 CALL QENTER('GET_FSO_AO') 2273 IF (ND .GT. 2) THEN 2274 CALL QUIT('ND .gt. fixed dimension of 2') 2275 END IF 2276C 2277C READ AND DISTRIBUTE SPIN-ORBIT AO INTEGRALS 2278C 2279 DO ID=1,ND 2280 IF (TRIPLET(ID)) THEN 2281 LEFT(ID)= D1 2282 RIGHT(ID)=D2 2283 ELSE 2284 LEFT(ID) =D2 2285 RIGHT(ID)=D1 2286 END IF 2287 END DO 2288 2289 LUHSOAO=-1 2290 CALL GPOPEN( 2291 & LUHSOAO,'AO2SOINT','OLD','SEQUENTIAL','UNFORMATTED',-1,.FALSE. 2292 & ) 2293 REWIND (LUHSOAO) 2294 2295 CALL ERIBUF_INI ! set NIBUF, NBITS, IBIT1, IBIT2 2296 LBUF = LBUF_alloc 2297 2298 LEN_BUF4 = 2*LBUF + NIBUF*LBUF + 1 ! BUF(LBUF), IBUF4(NIBUF,LBUF), LENGTH4 in integer*4 units 2299 2300 CALL MOLLAB('AO2SOINT',LUHSOAO,LUPRI) 2301 IF (WORD(1:1).EQ.'X') COMPONENT = 1 2302 IF (WORD(1:1).EQ.'Y') COMPONENT = 2 2303 IF (WORD(1:1).EQ.'Z') COMPONENT = 3 2304 2305 IF (IPRHSO .GT. 100 .OR. HSODEBUG) THEN 2306 WRITE(LUPRI,*) 'Hello from GET_FSO_AO. WORD: ',WORD,COMPONENT 2307 WRITE(LUPRI,*) ' TRIPLET? ',TRIPLET(1:ND) 2308 WRITE(LUPRI,*) ' LBUF, NIBUF, NBITS, IBIT1, IBIT2, LEN_BUF4', 2309 & LBUF, NIBUF, NBITS, IBIT1, IBIT2, LEN_BUF4 2310 END IF 2311 2312 NRECS = 0 2313cF90 READFILE: DO 2314 100 CONTINUE 2315 CALL READI4(LUHSOAO, LEN_BUF4, BUF) 2316 NRECS = NRECS + 1 2317 CALL AOLAB4(BUF(LBUF+1),LBUF,NIBUF,NBITS,IINDX4,LENGTH) 2318 IF (HSODEBUG) THEN 2319 IF (NRECS .LE. 2) THEN 2320 write (lupri,*) 'AO2SOINT record no., length:',NRECS, length 2321 do ib = 1,length 2322 write(lupri,*) iindx4(1:4,ib), buf(ib) 2323 end do 2324 END IF 2325 END IF 2326cF90 IF (LENGTH.LE.0) EXIT 2327 IF (LENGTH.LE.0) GO TO 300 2328cF90 BUFFER: DO IB=1,LENGTH 2329 DO IB=1,LENGTH 2330 K = IINDX4(3,IB) 2331 L = IINDX4(4,IB) 2332 IF (L.EQ.0) THEN 2333 ICOOR = K 2334cF90 CYCLE BUFFER 2335 GO TO 200 2336 END IF 2337cF90 IF (ICOOR.NE.COMPONENT) CYCLE BUFFER 2338 IF (ICOOR.NE.COMPONENT) GO TO 200 2339 I = IINDX4(1,IB) 2340 J = IINDX4(2,IB) 2341 GIJKL=BUF(IB) 2342 IF (I.EQ.J) GIJKL = DP5*GIJKL 2343 XIJKL=D1P5*GIJKL 2344 DO ID=1,ND 2345 DIJ = LEFT(ID)*(D(I,J,ID) + D(J,I,ID)) 2346 DKL = RIGHT(ID)*(D(L,K,ID) - D(K,L,ID)) 2347 F(I,J,ID) = F(I,J,ID) + GIJKL*DKL 2348 F(J,I,ID) = F(J,I,ID) + GIJKL*DKL 2349 F(K,L,ID) = F(K,L,ID) + GIJKL*DIJ 2350 F(L,K,ID) = F(L,K,ID) - GIJKL*DIJ 2351 F(I,L,ID) = F(I,L,ID) - XIJKL*D(J,K,ID) 2352 F(L,I,ID) = F(L,I,ID) + XIJKL*D(K,J,ID) 2353 F(I,K,ID) = F(I,K,ID) + XIJKL*D(J,L,ID) 2354 F(K,I,ID) = F(K,I,ID) - XIJKL*D(L,J,ID) 2355 F(J,L,ID) = F(J,L,ID) - XIJKL*D(I,K,ID) 2356 F(L,J,ID) = F(L,J,ID) + XIJKL*D(K,I,ID) 2357 F(J,K,ID) = F(J,K,ID) + XIJKL*D(I,L,ID) 2358 F(K,J,ID) = F(K,J,ID) - XIJKL*D(L,I,ID) 2359 END DO 2360cF90 END DO BUFFER 2361 200 CONTINUE 2362 END DO 2363cF90 END DO READFILE 2364 GO TO 100 2365 300 CONTINUE 2366C 2367 CALL GPCLOSE(LUHSOAO,'KEEP') 2368 CALL QEXIT('GET_FSO_AO') 2369cF90 END SUBROUTINE GET_FSO_AO 2370 END 2371