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 19#ifdef REVLOG 20=========================================================================== 21Revision 1.2 2000/05/24 19:09:06 hjj 22inserted Dalton header 23some changes for triplet response with CSF 24=========================================================================== 25#endif 26C============================================================================= 27C /* Deck E3IEF */ 28C============================================================================= 29 SUBROUTINE E3IEF(VECA, VEC1, VEC2,ETRS,XINDX,ZYM1,ZYM2, 30 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 31 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP) 32C 33C Purpose: 34C Outer driver routine for IEF-PCM solvent contribution 35C to E[3] times two vectors. 36C Derived from its "older brother", rspsol1.F 37C 38C Main comments: 39C 1) POPGET contains two cycles over tesserae because the terms 40C are no longer diagonal: the charge on each tessera 41C does not depend only on the potential on that tessera 42C 2) For the time beeing I have not implemented the C3IEF 43C for the memory efficent SCF calculations 44C 45#include "implicit.h" 46C 47#include "maxorb.h" 48#include "priunit.h" 49#include "inforb.h" 50#include "infdim.h" 51#include "infinp.h" 52#include "infvar.h" 53#include "infrsp.h" 54#include "infpri.h" 55#include "rspprp.h" 56#include "infcr.h" 57C 58 DIMENSION ETRS(KZYVR),XINDX(*) 59 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 60 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 61 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 62 DIMENSION MJWOP(2,MAXWOP,8) 63 INTEGER DBGFLG(10) 64C 65 NSIM = 1 66 KFREE = 1 67 CALL MEMGET('REAL',KTA ,N2ORBX,WRK,KFREE,LFREE) 68 CALL MEMGET('REAL',KTB ,N2ORBX,WRK,KFREE,LFREE) 69 CALL MEMGET('REAL',KTB1,N2ORBX,WRK,KFREE,LFREE) 70 CALL MEMGET('REAL',KTB2,N2ORBX,WRK,KFREE,LFREE) 71 CALL MEMGET('REAL',KTC1,N2ORBX,WRK,KFREE,LFREE) 72 CALL MEMGET('REAL',KTC2,N2ORBX,WRK,KFREE,LFREE) 73 CALL MEMGET('REAL',KTD1,N2ORBX,WRK,KFREE,LFREE) 74 CALL MEMGET('REAL',KTD2,N2ORBX,WRK,KFREE,LFREE) 75 CALL MEMGET('REAL',KTE ,N2ORBX,WRK,KFREE,LFREE) 76C 77 CALL DZERO(WRK(KTA), N2ORBX) 78 CALL DZERO(WRK(KTB), N2ORBX) 79 CALL DZERO(WRK(KTB1),N2ORBX) 80 CALL DZERO(WRK(KTB2),N2ORBX) 81 CALL DZERO(WRK(KTC1),N2ORBX) 82 CALL DZERO(WRK(KTC2),N2ORBX) 83 CALL DZERO(WRK(KTD1),N2ORBX) 84 CALL DZERO(WRK(KTD2),N2ORBX) 85 CALL DZERO(WRK(KTE), N2ORBX) 86C 87 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYM1,MJWOP) 88 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,ZYM2,MJWOP) 89C 90C DBGFLG initialization 91C 1 2 3 4 5 6 7 8 9 10 92C DATA DBGFLG/-1,-2,-3,-4,-5,-6,-7,-8,-9,-10/ 93 DATA DBGFLG/ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10/ 94C 95C VECA is only available if E3TEST is set 96C 97 VAL = DDOT(KZYVR,VECA,1,ETRS,1) 98C 99 CALL PCASE1(VECA, VEC1, VEC2,WRK(KTA),WRK(KTB), 100 * ETRS,XINDX,ZYM1,ZYM2, 101 * DEN1,UDV,WRK(KFREE),LFREE, 102 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 103 * DBGFLG) 104C 105 CALL PCASE2(VECA,VEC1,VEC2,WRK(KTC1), 106 * ETRS,XINDX,ZYM1,ZYM2, 107 * DEN1,UDV,WRK(KFREE),LFREE, 108 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 109 * DBGFLG) 110C 111 CALL PCASE2(VECA,VEC2,VEC1,WRK(KTC2), 112 * ETRS,XINDX,ZYM2,ZYM1, 113 * DEN1,UDV,WRK(KFREE),LFREE, 114 * KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1,CMO,MJWOP, 115 * DBGFLG) 116C 117 CALL PCASE3(VECA,VEC1,VEC2,WRK(KTB),WRK(KTB1),WRK(KTC1),WRK(KTD1), 118 * ETRS,XINDX,ZYM1,ZYM2, 119 * DEN1,UDV,WRK(KFREE),LFREE, 120 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP) 121C 122 CALL PCASE3(VECA,VEC2,VEC1,WRK(KTB),WRK(KTB2),WRK(KTC2),WRK(KTD2), 123 * ETRS,XINDX,ZYM2,ZYM1, 124 * DEN1,UDV,WRK(KFREE),LFREE, 125 * KZYVR,KZYV2,KZYV1,IGRSYM,ISYMV2,ISYMV1,CMO,MJWOP) 126C 127 CALL PCASE4(VECA, VEC1, VEC2,WRK(KTA),WRK(KTB1), 128 * WRK(KTB2),WRK(KTC1),WRK(KTC2),WRK(KTE), 129 * ETRS,XINDX,ZYM1,ZYM2, 130 * DEN1,UDV,WRK(KFREE),LFREE, 131 * KZYVR,KZYV1,KZYV2,IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 132 * DBGFLG) 133C 134 IF (CRCAL .OR. E3TEST) THEN 135 WRITE (LUPRI,'(A,F20.12)') 'Total contribution from E3IEF:', 136 * DDOT(KZYVR,VECA,1,ETRS,1) - VAL 137 END IF 138C 139 RETURN 140 END 141C============================================================================= 142C /* Deck POPGET */ 143C============================================================================= 144 SUBROUTINE POPGET(ZYM1,ZYM2,ZYM3,IKLVL,OVLAP,ISYMDN, 145 * ISYMV1,ISYMV2,ISYMV3,TOP,DEN1,NUCFLG, 146 * CMO,WRK,LWRK,IPRLVL,PRSTR,DFCTR, 147 * ADDFLG) 148C 149C Output: 150C 151C TOP = sum(l,m) f(l) <L| T(l,m)(k1,k2,..) |R> TE(l,m) 152C 153C POP = sum_i <L|q_i(k1,k2,..)|R> V_i 154C 155C Input: 156C 157C OVLAP = <L|R>, DEN1 = <L|...|R> 158C NUCFLG indicates if T(l,m) = TN(l,m) or T(l,m) = TE(l,m) 159C IKLVL is the number of times T(l,m) is to be one-index tranformed 160C 161C PCM Notes 162C 1) The two T (TE and TN) are now replaced by charges q and 163C potentials V 164C 2) Potentials are always the electronic potentials derived from 165C the MOs (J1 integrals) 166C 3) Charges can be - qtot = qed + qei +qn (total charges) 167C - qe = qed + qei (total electronic charges) 168C - qed (dynamic electronic charges) 169C 4) All transformations are carried out on the charge part and 170C never on the potential part but since the final term is always Vq 171C this is not relevant 172C 5) For the response calculation starting from a nonequilibrium 173C state some coding will be needed in order to get all the charges 174C properly: qei in particular and the total charges 175 176C 177C 178#include "implicit.h" 179#include "dummy.h" 180C 181#include "maxorb.h" 182#include "mxcent.h" 183#include "pcmdef.h" 184#include "priunit.h" 185#include "infdim.h" 186#include "orgcom.h" 187#include "inforb.h" 188#include "infpri.h" 189#include "infpar.h" 190#include "inftap.h" 191#include "infinp.h" 192#include "infrsp.h" 193#include "pcm.h" 194#include "pcmlog.h" 195C 196Clf PCM NUCFLG 197C NUCFLG = 0 MULTIPLY BY ELECTRONIC GHARGES 198C NUCFLG = 1 MULTIPLY BY NUCLEAR GHARGES 199C NUCFLG = 2 MULTIPLY BY ELECTRONIC + NUCLEAR GHARGES 200 CHARACTER*(*)PRSTR 201 LOGICAL FNDLAB, CLCCHG, EXP1VL, TOFILE, TRIMAT 202Clf 203 INTEGER ADDFLG 204C 205 DIMENSION DEN1(NASHDI,NASHDI), INTREP(9*MXCENT), INTADR(9*MXCENT) 206 DIMENSION WRK(*), CMO(*) 207 DIMENSION ZYM1(NORBT,NORBT),ZYM2(NORBT,NORBT),ZYM3(NORBT,NORBT) 208 DIMENSION TOP(N2ORBX) 209 CHARACTER*8 LABINT(9*MXCENT) 210C 211 NSIM = 1 212 213C When IKLVL is passed as -1 we initialize CLCCHG as false 214 IF(IKLVL.LT.0) THEN 215 CLCCHG = .FALSE. 216 ELSE 217 CLCCHG = .TRUE. 218 END IF 219 IF (IKLVL.LE.0) ISYM = 1 220 IF (IKLVL.EQ.1) ISYM = ISYMV1 221 IF (IKLVL.EQ.2) ISYM = MULD2H(ISYMV1,ISYMV2) 222 IF (IKLVL.EQ.3) ISYM = MULD2H(MULD2H(ISYMV1,ISYMV2),ISYMV3) 223 ISYMT = MULD2H(ISYM,ISYMDN) 224C 225 KFREE = 1 226 CALL MEMGET('REAL',KCHGAO,NSYM*NNBASX,WRK,KFREE,LWRK) 227 CALL MEMGET('REAL',KUCMO ,NORBT*NBAST,WRK,KFREE,LWRK) 228 CALL MEMGET('REAL',KTELM ,N2ORBX,WRK,KFREE,LWRK) 229 CALL MEMGET('REAL',KTLMA ,N2ORBX,WRK,KFREE,LWRK) 230 CALL MEMGET('REAL',KTLMB ,N2ORBX,WRK,KFREE,LWRK) 231 CALL MEMGET('REAL',KTLMC ,N2ORBX,WRK,KFREE,LWRK) 232 CALL MEMGET('REAL',KTCHG ,NTS,WRK,KFREE,LWRK) 233 CALL MEMGET('REAL',KQCHG ,NTS,WRK,KFREE,LWRK) 234 CALL MEMGET('REAL',KMULT ,NTS,WRK,KFREE,LWRK) 235C 236 CALL DZERO(WRK(KTLMC),N2ORBX) 237Clf 238C 239C Unpack symmetry blocked CMO 240C 241 CALL UPKCMO(CMO,WRK(KUCMO)) 242C 243C Loop over tesserae 244C 245ckr REWIND (LUPROP) 246C 247C One-index transform charges IKLVL times. 248C The result will be in WRK(KTLMA) and of symmetry ISYM. 249C (ISYM should equal ISYMDN.) 250C 251C 252C We calculate charges only if necessary 253C Otherwise we use QSE or QSN 254C 255 IF (CLCCHG) THEN 256 CALL DZERO(WRK(KTCHG),NTS) 257 XI = DIPORG(1) 258 YI = DIPORG(2) 259 ZI = DIPORG(3) 260#if defined (VAR_MPI) 261 IF (NODTOT .GE. 1) THEN 262 CALL J2XP(IKLVL,ISYMV1,ISYMV2,ISYMV3,ISYMDN,ZYM1,ZYM2,ZYM3, 263 & WRK(KTCHG),OVLAP,WRK(KUCMO),DEN1,WRK(KFREE),LWRK) 264 ELSE 265#endif 266 DO ITS=1,NTSIRR 267 CALL DZERO(WRK(KCHGAO),NNBASX) 268 L = 1 269 NCOMP = NSYM 270 DIPORG(1) = XTSCOR(ITS) 271 DIPORG(2) = YTSCOR(ITS) 272 DIPORG(3) = ZTSCOR(ITS) 273 EXP1VL = .FALSE. 274 TOFILE = .FALSE. 275 KPATOM = 0 276 TRIMAT = .TRUE. 277 CALL GET1IN(WRK(KCHGAO),'NPETES ',NCOMP,WRK(KFREE),LWRK, 278 & LABINT,INTREP,INTADR,L,TOFILE,KPATOM,TRIMAT, 279 & DUMMY,EXP1VL,DUMMY,IPRRSP) 280 JCHGAO = KCHGAO 281 DO ILOP = 1, NSYM 282 ISYM = ILOP 283 JTS = (ILOP - 1)*NTSIRR + ITS 284 CALL UTHU(WRK(JCHGAO),WRK(KTLMA),WRK(KUCMO),WRK(KFREE), 285 $ NBAST,NORBT) 286 CALL DSPTSI(NORBT,WRK(KTLMA),WRK(KTELM)) 287C with no transformation we copy the array from KTELM to KTLMA 288 IF (IKLVL.EQ.0) THEN 289 CALL DZERO(WRK(KTLMA),N2ORBX) 290 CALL DCOPY(N2ORBX,WRK(KTELM),1,WRK(KTLMA),1) 291 END IF 292C First transformation of charges: q^e_{ab} --> q^e_{ab}({}^1\kappa) 293 IF (IKLVL.GE.1) THEN 294 CALL DZERO(WRK(KTLMA),N2ORBX) 295 CALL OITH1(ISYMV1,ZYM1,WRK(KTELM),WRK(KTLMA),ISYM) 296 ISYM = MULD2H(ISYM,ISYMV1) 297 END IF 298C Second transformation of charges: q^e_{ab}({}^1\kappa) --> q^e_{ab}({}^1\kappa {}^2\kappa) 299 IF (IKLVL.GE.2) THEN 300 CALL DZERO(WRK(KTLMB),N2ORBX) 301 CALL OITH1(ISYMV2,ZYM2,WRK(KTLMA),WRK(KTLMB),ISYM) 302 CALL DCOPY(N2ORBX,WRK(KTLMB),1,WRK(KTLMA),1) 303 ISYM = MULD2H(ISYM,ISYMV2) 304 END IF 305C Third transformation of charges: hope you can figure out the formula..... 306 IF (IKLVL.GE.3) THEN 307 CALL DZERO(WRK(KTLMA),N2ORBX) 308 CALL OITH1(ISYMV3,ZYM3,WRK(KTLMB),WRK(KTLMA),ISYM) 309 ISYM = MULD2H(ISYM,ISYMV3) 310 END IF 311C Contract transformed charges with the density 312 CALL MELONE(WRK(KTLMA),ISYM,DEN1,OVLAP,FACT,200,'POPGET') 313 WRK(KTCHG + JTS - 1) = FACT 314 JCHGAO = JCHGAO + NNBASX 315 END DO 316 END DO 317#if defined (VAR_MPI) 318 END IF 319#endif 320 CALL V2Q(WRK(KFREE),WRK(KTCHG),WRK(KQCHG),QTEXS,NEQRSP) 321 END IF 322 IF(LUPCMD .GT. 0) THEN 323 CALL GPCLOSE(LUPCMD,'KEEP') 324 END IF 325C 326C Make charges of the n-index transformed potentials 327C 328C 329C Construction of the operator (J, X, or whatever....) 330C 331 IF (CLCCHG) THEN 332 CALL DCOPY(NTSIRR,WRK(KQCHG+(ISYMT - 1)*NTSIRR),1,WRK(KMULT),1) 333 ELSE IF (NUCFLG .EQ. 0) THEN 334 IF (NEQRSP.AND.(ABS(ADDFLG).EQ.1)) THEN 335 CALL DCOPY(NTS,QSENEQ,1,WRK(KMULT),1) 336 ELSE 337 CALL DCOPY(NTS,QSE,1,WRK(KMULT),1) 338 END IF 339 CALL DSCAL(NTS,-1.0D0,WRK(KMULT),1) 340 ELSE IF (NUCFLG .EQ. 1) THEN 341 CALL DCOPY(NTS,QSN,1,WRK(KMULT),1) 342 CALL DSCAL(NTS,-1.0D0,WRK(KMULT),1) 343 ELSE IF (NUCFLG .EQ. 2) THEN 344 CALL DCOPY(NTS,QSN,1,WRK(KMULT),1) 345 CALL DAXPY(NTS,1.0D0,QSE,1,WRK(KMULT),1) 346 CALL DSCAL(NTS,-1.0D0,WRK(KMULT),1) 347 ELSE 348 CALL QUIT('CASE NOT DEFINED IN POPGET') 349 END IF 350 NOSIM = 1 351 CALL DZERO(WRK(KCHGAO),NSYM*NNBASX) 352 CALL J1INT(WRK(KMULT),.FALSE.,WRK(KCHGAO),NOSIM,.FALSE.,'NPETES ', 353 & ISYMT,WRK(KFREE),LWRK) 354 CALL UTHU(WRK(KCHGAO),WRK(KTLMA),WRK(KUCMO),WRK(KFREE), 355 & NBAST,NORBT) 356 CALL DSPTSI(NORBT,WRK(KTLMA),WRK(KTLMC)) 357 IF (ADDFLG .GT. 0) CALL DAXPY(N2ORBX,DFCTR,WRK(KTLMC),1,TOP,1) 358C 359 DIPORG(1) = XI 360 DIPORG(2) = YI 361 DIPORG(3) = ZI 362C 363 IF (IPRRSP.GE.IPRLVL) THEN 364 WRITE(LUPRI,'(/3A,2D22.14)') 'Norm of TOPGET in ', PRSTR, 365 * ' : ',DNRM2(N2ORBX,wrk(ktlmc),1),DNRM2(N2ORBX,TOP,1) 366 END IF 367C 368 RETURN 369C 370 END 371C============================================================================= 372C /* Deck PCASE1 */ 373C============================================================================= 374 SUBROUTINE PCASE1(VECA, VEC1, VEC2,TA,TB, 375 * ETRS,XINDX,ZYM1,ZYM2, 376 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 377 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP, 378 * DBGFLG) 379#include "implicit.h" 380#include "dummy.h" 381C 382 PARAMETER ( D1 = 1.0D0, DM1 = -1.0D0, D2 = 2.0D0, DM2 = -2.0D0, 383 $ DP5= 0.5D0) 384C 385#include "maxorb.h" 386C#include "pcmdef.h" 387#include "infdim.h" 388#include "inforb.h" 389#include "wrkrsp.h" 390#include "infrsp.h" 391#include "infpri.h" 392#include "infvar.h" 393#include "qrinf.h" 394#include "infspi.h" 395#include "infden.h" 396#include "infinp.h" 397C 398 DIMENSION ETRS(KZYVR),XINDX(*) 399 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 400 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 401 DIMENSION TA(NORBT,NORBT),TB(NORBT,NORBT) 402 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 403 DIMENSION MJWOP(2,MAXWOP,8) 404 INTEGER DBGFLG(10) 405C 406 LOGICAL TDM, LREF, NORHO2 407C 408C Initialise variables 409C 410 JSPIN = 0 411 TDM = .TRUE. 412 KFREE = 1 413 NORHO2 = .TRUE. 414 NSIM = 1 415C 416C TA = 2*sum(l,m) f(l) <0| TE(l,m) |0> TE(l,m) 417C PCM sum_its QSE(its)*V(its) 418C 419Clf IF ((ISYMV1.EQ.ISYMV2) .AND. (MZCONF(ISYMV1).GT.0)) THEN 420 CALL DZERO(TA,N2ORBX) 421 IKLVL = -1 422Clf Call n.1 423Clf 424 CALL POPGET(DUMMY,DUMMY,DUMMY,IKLVL,D1,1, 425 * IDUMMY,IDUMMY,IDUMMY,TA,UDV,0, 426 * CMO,WRK(KFREE),LFREE,100,'TA',D1, 427 * DBGFLG(1)) 428C 429C Scaling not needed for PCM 430Clf CALL DSCAL(N2ORBX,D2,TA,1) 431Clf END IF 432C 433C PCM TB = 2*sum(its) [qse(its)+qsn(its)]*TE(its) (TE is the potential operator). 434C 435C Onsager:TB = 2*sum(l,m) f(l)*( <0|TE(l,m)|0> - Tn(l,m) )*TE(l,m) 436C 437Clf IF (NEQRSP) THEN 438Clf EPSTOP = EPSINF 439Clf ELSE 440Clf EPSTOP = EPS 441Clf END IF 442 CALL DZERO(TB,N2ORBX) 443Clf Call n.2 444Clf 445 IKLVL = -1 446 CALL POPGET(DUMMY,DUMMY,DUMMY,IKLVL,DUMMY,1, 447 * IDUMMY,IDUMMY,IDUMMY,TB,DUMMY,1, 448 * CMO,WRK(KFREE),LFREE,100,'TB',D1, 449 * DBGFLG(2)) 450Clf CALL DSCAL(N2ORBX,DM2,TB,1) 451Clf CALL DSCAL(N2ORBX,DM1,TB,1) 452Clf Call n.3 453Clf 454 IKLVL = -1 455 CALL POPGET(DUMMY,DUMMY,DUMMY,IKLVL,D1,1, 456 * IDUMMY,IDUMMY,IDUMMY,TB,UDV,0, 457 * CMO,WRK(KFREE),LFREE,100,'TB',D1, 458 * DBGFLG(3)) 459C 460 IF (MZCONF(ISYMV1) .EQ. 0 .OR. MZCONF(ISYMV2) .EQ. 0) RETURN 461C 462C / <01L| [qj,TB] |02R> + <02L| [qj,TB] |01R> \ 463C | 0 | 464C | <01L| [qj+,TB] |02R> + <02L| [qj+,TB] |01R> | 465C \ 0 / 466C 467C Construct <01L|..|02R> + <02L|..|01R> density 468C 469 ILSYM = MULD2H(IREFSY,ISYMV1) 470 IRSYM = MULD2H(IREFSY,ISYMV2) 471 NCL = MZCONF(ISYMV1) 472 NCR = MZCONF(ISYMV2) 473 KZVARL = MZYVAR(ISYMV1) 474 KZVARR = MZYVAR(ISYMV2) 475 LREF = .FALSE. 476 ISYMDN = MULD2H(ILSYM,IRSYM) 477 CALL DZERO(DEN1,NASHT*NASHT) 478 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 479 * VEC1,VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM,NORHO2, 480 * XINDX,WRK,KFREE,LFREE,LREF) 481C 482C Make the gradient 483C 484 IF ( MZWOPT(IGRSYM) .GT. 0 ) THEN 485 CALL ORBSX(NSIM,IGRSYM,KZYVR,ETRS,TB,OVLAP,ISYMDN, 486 * DEN1,MJWOP,WRK(KFREE),LFREE) 487 END IF 488C 489 CALL PRIRES(ETRS,VECA,IGRSYM,'PCASE1') 490C 491 RETURN 492 END 493C============================================================================= 494C /* Deck PCASE2 */ 495C============================================================================= 496 SUBROUTINE PCASE2(VECA, VEC1, VEC2,TC1, 497 * ETRS,XINDX,ZYM1,ZYM2, 498 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 499 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP,DBGFLG) 500#include "implicit.h" 501#include "dummy.h" 502C 503 PARAMETER ( D1 = 1.0D0, D2 = 2.0D0 , DM1 = -1.0D0 ) 504C 505#include "maxorb.h" 506#include "infdim.h" 507#include "inforb.h" 508#include "wrkrsp.h" 509#include "infrsp.h" 510#include "infpri.h" 511#include "infvar.h" 512#include "qrinf.h" 513#include "infspi.h" 514#include "infden.h" 515#include "infinp.h" 516C 517 DIMENSION ETRS(KZYVR),XINDX(*) 518 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 519 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 520 DIMENSION TC1(NORBT,NORBT) 521 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 522 DIMENSION MJWOP(2,MAXWOP,8) 523 INTEGER DBGFLG(10) 524C 525 LOGICAL TDM, LREF, NORHO2 526C 527C Initialise variables 528C 529 JSPIN = 0 530 TDM = .TRUE. 531 IPRONE = 200 532 KFREE = 1 533 NORHO2 = .TRUE. 534 NSIM = 1 535C 536 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 537 CALL GETREF(WRK(KCREF),MZCONF(1)) 538C 539C Onsager: TC1 = 2*sum(l,m) f(l) <0| TE(l,m) (k1) |0> TE(l,m) + ... 540C PCM: TC1 = 2*sum(its) <0| q_e(its)(k1) |0> V(its) 541C 542 CALL DZERO(TC1,N2ORBX) 543Clf CALL OUTPUT(zym1,1,NORBT,1,NORBT,NORBT,NORBT,1, 544Clf & 2) 545 IF (MZWOPT(ISYMV1).GT.0) THEN 546Clf Call n.4 547Clf 548 IKLVL = 1 549 CALL POPGET(ZYM1,DUMMY,DUMMY,IKLVL,D1,1, 550 * ISYMV1,IDUMMY,IDUMMY,TC1,UDV,IDUMMY, 551 * CMO,WRK(KFREE),LFREE,100,'TC1 cont1',DM1, 552 * DBGFLG(4)) 553 END IF 554C 555C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m) |0> + <0| TE(l,m) |01R> ) TE(l,m) 556C 557 IF (MZCONF(ISYMV1).GT.0) THEN 558C 559C Construct the density matrix <01L|..|0> + <0|..|01R> 560C 561 ILSYM = IREFSY 562 IRSYM = MULD2H(IREFSY,ISYMV1) 563 NCL = MZCONF(1) 564 NCR = MZCONF(ISYMV1) 565 KZVARL = MZCONF(1) 566 KZVARR = MZYVAR(ISYMV1) 567 LREF = .TRUE. 568 ISYMDN = MULD2H(ILSYM,IRSYM) 569 CALL DZERO(DEN1,NASHT*NASHT) 570 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 571 * WRK(KCREF),VEC1,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 572 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 573C 574Clf Call n.5 575Clf 576 IKLVL = 0 577 CALL POPGET(DUMMY,DUMMY,DUMMY,IKLVL,OVLAP,ISYMDN, 578 * IDUMMY,IDUMMY,IDUMMY,TC1,DEN1,0, 579 * CMO,WRK(KFREE),LFREE,100,'TC1 cont2',DM1, 580 * DBGFLG(5)) 581 END IF 582C 583Clf CALL DSCAL(N2ORBX,D2,TC1,1) 584Clf CALL DSCAL(N2ORBX,DM1,TC1,1) 585C 586 IF (MZCONF(ISYMV2).LE.0) RETURN 587C 588C 589C / 0 \ 590C | Sj(2) | * <0| TC1 |0> 591C | 0 | 592C \ Sj(2)' / 593C 594 IF (IGRSYM.EQ.ISYMV2) THEN 595 OVLAP = D1 596 CALL MELONE(TC1,1,UDV,OVLAP,FACT,IPRONE,'FACT in PCASE2') 597 NZCONF = MZCONF(IGRSYM) 598 NZVAR = MZVAR(IGRSYM) 599 CALL DAXPY(NZCONF,FACT,VEC2,1,ETRS,1) 600 CALL DAXPY(NZCONF,FACT,VEC2(NZVAR+1),1,ETRS(NZVAR+1),1) 601 END IF 602C 603 CALL PRIRES(ETRS,VECA,IGRSYM,'PCASE2') 604C 605 RETURN 606 END 607C============================================================================= 608C /* Deck PCASE3 */ 609C============================================================================= 610 SUBROUTINE PCASE3(VECA, VEC1, VEC2,TB,TB1,TC1,TD1, 611 * ETRS,XINDX,ZYM1,ZYM2, 612 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 613 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP) 614#include "implicit.h" 615#include "dummy.h" 616C 617 PARAMETER ( D1 = 1.0D0 ) 618C 619#include "maxorb.h" 620#include "infdim.h" 621#include "inforb.h" 622#include "wrkrsp.h" 623#include "infrsp.h" 624#include "infpri.h" 625#include "infvar.h" 626#include "qrinf.h" 627#include "infspi.h" 628#include "infden.h" 629#include "infinp.h" 630C 631 DIMENSION ETRS(KZYVR),XINDX(*) 632 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 633 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 634 DIMENSION TB(NORBT,NORBT),TB1(NORBT,NORBT) 635 DIMENSION TC1(NORBT,NORBT),TD1(NORBT,NORBT) 636 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 637 DIMENSION MJWOP(2,MAXWOP,8) 638C 639 LOGICAL LCON, LORB 640 LOGICAL TDM, LREF, NORHO2 641C 642C Initialise variables 643C 644 JSPIN = 0 645 TDM = .TRUE. 646 KFREE = 1 647 IPRONE = -1 648 NORHO2 = .TRUE. 649 NSIM = 1 650C 651C TD1 = TB1 + TC1, TB1 = TB(k1) 652C 653 CALL DZERO(TB1,N2ORBX) 654 CALL DZERO(TD1,N2ORBX) 655 CALL OITH1(ISYMV1,ZYM1,TB,TB1,1) 656C 657Clf Moved the check before the two DAXPY: TD1 is used only in this subroutine) 658C 659 IF (MZCONF(ISYMV2).LE.0) RETURN 660C 661 CALL DAXPY(N2ORBX,D1,TB1,1,TD1,1) 662 CALL DAXPY(N2ORBX,D1,TC1,1,TD1,1) 663C the check was here 664 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 665 CALL GETREF(WRK(KCREF),MZCONF(1)) 666C 667C MCSCF to be checked!!! 668C 669C / <0| [qj,TD1] |02R> + <02L| [qj,TD1] |0> \ 670C | <j| TD1 |02R> | 671C | <0| [qj+,TD1] |02R> + <02L| [qj+,TD1] |0> | 672C \ -<02L| TD1 |j> / 673C 674C Construct the density matrix <02L|..|0> + <0|..|02R> 675C 676 ILSYM = IREFSY 677 IRSYM = MULD2H(IREFSY,ISYMV2) 678 NCL = MZCONF(1) 679 NCR = MZCONF(ISYMV2) 680 KZVARL = MZCONF(1) 681 KZVARR = MZYVAR(ISYMV2) 682 LREF = .TRUE. 683 ISYMDN = MULD2H(ILSYM,IRSYM) 684 CALL DZERO(DEN1,NASHT*NASHT) 685 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 686 * WRK(KCREF),VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 687 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 688C 689C Make the gradient 690C 691 ISYMST = MULD2H(IGRSYM,IREFSY) 692 IF ( ISYMST .EQ. IREFSY ) THEN 693 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 694 ELSE 695 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 696 END IF 697 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 698 LREF = .FALSE. 699 NZYVEC = MZYVAR(ISYMV2) 700 NZCVEC = MZCONF(ISYMV2) 701 CALL RSP1GR(NSIM,KZYVR,IDUMMY,JSPIN,IGRSYM,JSPIN,ISYMV2,ETRS, 702 * VEC2,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,TD1, 703 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 704C 705 CALL PRIRES(ETRS,VECA,IGRSYM,'PCASE3') 706C 707 RETURN 708 END 709C============================================================================= 710C /* Deck PCASE4 */ 711C============================================================================= 712 SUBROUTINE PCASE4(VECA, VEC1, VEC2,TA,TB1,TB2,TC1,TC2,TE, 713 * ETRS,XINDX,ZYM1,ZYM2, 714 * DEN1,UDV,WRK,LFREE,KZYVR,KZYV1,KZYV2, 715 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP,DBGFLG) 716#include "implicit.h" 717#include "dummy.h" 718C 719 PARAMETER ( D1 = 1.0D0, D2 = 2.0D0, DH = 0.5D0 , DMP5 = -0.5D0, 720 $ DM1 = -1.0D0) 721C 722#include "maxorb.h" 723#include "infdim.h" 724#include "inforb.h" 725#include "wrkrsp.h" 726#include "infrsp.h" 727#include "infpri.h" 728#include "infvar.h" 729#include "qrinf.h" 730#include "infspi.h" 731#include "infden.h" 732#include "infinp.h" 733#include "priunit.h" 734C 735 DIMENSION ETRS(KZYVR),XINDX(*) 736 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 737 DIMENSION ZYM1(*),ZYM2(*),WRK(*),CMO(*) 738 DIMENSION TA(NORBT,NORBT),TE(NORBT,NORBT) 739 DIMENSION TB1(NORBT,NORBT),TB2(NORBT,NORBT) 740 DIMENSION TC1(NORBT,NORBT),TC2(NORBT,NORBT) 741 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 742 DIMENSION MJWOP(2,MAXWOP,8) 743 INTEGER DBGFLG(10) 744C 745 LOGICAL LCON, LORB 746 LOGICAL TDM, LREF, NORHO2 747C 748C Initialise variables 749C 750 JSPIN = 0 751 TDM = .TRUE. 752 KFREE = 1 753 IPRONE = 100 754 NORHO2 = .TRUE. 755 NSIM = 1 756C 757 CALL MEMGET('REAL',KCREF,MZCONF(1),WRK,KFREE,LFREE) 758 CALL GETREF(WRK(KCREF),MZCONF(1)) 759C 760C TE = 1/2 * TB1(k2) + 1/2 * TB2(k1) + TC1(k2) + TC2(k1) + ... 761C 762 CALL DZERO(TE,N2ORBX) 763 CALL OITH1(ISYMV2,ZYM2,TB1,TE,ISYMV1) 764 CALL OITH1(ISYMV1,ZYM1,TB2,TE,ISYMV2) 765 CALL DSCAL(N2ORBX,DH,TE,1) 766 CALL OITH1(ISYMV2,ZYM2,TC1,TE,ISYMV1) 767 CALL OITH1(ISYMV1,ZYM1,TC2,TE,ISYMV2) 768C 769C ... + ( S(1)S(2)' + S(2)S(1)' ) * TA + ... 770C 771 IF ((ISYMV1.EQ.ISYMV2) .AND. (MZCONF(ISYMV1).GT.0)) THEN 772 NZCONF = MZCONF(ISYMV1) 773 NZVAR = MZVAR(ISYMV1) 774 FACT = DDOT(NZCONF,VEC1,1,VEC2(NZVAR+1),1) + 775 * DDOT(NZCONF,VEC2,1,VEC1(NZVAR+1),1) 776 CALL DAXPY(N2ORBX,FACT,TA,1,TE,1) 777 END IF 778C 779C ... + sum(l,m) f(l) <0| TE(l,m)(k1,k2) |0> TE(l,m) 780C + sum(l,m) f(l) <0| TE(l,m)(k2,k1) |0> TE(l,m) + ... 781C 782 IF (MZWOPT(ISYMV1).GT.0 .AND. MZWOPT(ISYMV2).GT.0) THEN 783Clf Call n.6 784Clf 785 IKLVL = 2 786 CALL POPGET(ZYM1,ZYM2,DUMMY,IKLVL,D1,1, 787 * ISYMV1,ISYMV2,IDUMMY,TE,UDV,IDUMMY, 788 * CMO,WRK(KFREE),LFREE,100,'TE cont1a',DMP5, 789 * DBGFLG(6)) 790Clf Call n.7 791Clf 792 IKLVL = 2 793 CALL POPGET(ZYM2,ZYM1,DUMMY,IKLVL,D1,1, 794 * ISYMV2,ISYMV1,IDUMMY,TE,UDV,IDUMMY, 795 * CMO,WRK(KFREE),LFREE,100,'TE cont1b',DMP5, 796 * DBGFLG(7)) 797 END IF 798C 799C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m)(k2) |0> + 800C <0| TE(l,m)(k2) |01R> ) TE(l,m) + ... 801C 802CLF WE DO NOT NEED THIS FOR PCM 803C Put the factor two into one of the vectors. 804C CALL DSCAL(KZYV1,D2,VEC1,1) 805C CALL DSCAL(NORBT*NORBT,D2,ZYM1,1) 806C 807C 808 IF (MZCONF(ISYMV1).GT.0 .AND. MZWOPT(ISYMV2).GT.0) THEN 809C 810C Construct the density matrix <01L|..|0> + <0|..|01R> 811C 812 ILSYM = IREFSY 813 IRSYM = MULD2H(IREFSY,ISYMV1) 814 NCL = MZCONF(1) 815 NCR = MZCONF(ISYMV1) 816 KZVARL = MZCONF(1) 817 KZVARR = MZYVAR(ISYMV1) 818 LREF = .TRUE. 819 ISYMDN = MULD2H(ILSYM,IRSYM) 820 CALL DZERO(DEN1,NASHT*NASHT) 821 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 822 * WRK(KCREF),VEC1,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 823 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 824C 825 IKLVL = 1 826Clf Call n.8 827Clf 828 CALL POPGET(ZYM2,DUMMY,DUMMY,IKLVL,OVLAP,ISYMDN, 829 * ISYMV2,IDUMMY,IDUMMY,TE,DEN1,IDUMMY, 830 * CMO,WRK(KFREE),LFREE,100,'TE cont2a',DM1, 831 * DBGFLG(8)) 832 END IF 833C 834C ... + 2*sum(l,m) f(l) ( <02L| TE(l,m)(k1) |0> + 835C <0| TE(l,m)(k1) |02R> ) TE(l,m) + ... 836C 837C The factor two is already included in one of the vectors. 838C 839 IF (MZCONF(ISYMV2).GT.0 .AND. MZWOPT(ISYMV1).GT.0) THEN 840C 841C Construct the density matrix <02L|..|0> + <0|..|02R> 842C 843 ILSYM = IREFSY 844 IRSYM = MULD2H(IREFSY,ISYMV2) 845 NCL = MZCONF(1) 846 NCR = MZCONF(ISYMV2) 847 KZVARL = MZCONF(1) 848 KZVARR = MZYVAR(ISYMV2) 849 LREF = .TRUE. 850 ISYMDN = MULD2H(ILSYM,IRSYM) 851 CALL DZERO(DEN1,NASHT*NASHT) 852 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 853 * WRK(KCREF),VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM, 854 * NORHO2,XINDX,WRK,KFREE,LFREE,LREF) 855C 856Clf Call n.9 857Clf 858 IKLVL = 1 859 CALL POPGET(ZYM1,DUMMY,DUMMY,IKLVL,OVLAP,ISYMDN, 860 * ISYMV1,IDUMMY,IDUMMY,TE,DEN1,IDUMMY, 861 * CMO,WRK(KFREE),LFREE,100,'TE cont2b',DM1, 862 * DBGFLG(9)) 863 END IF 864C 865C ... + 2*sum(l,m) f(l) ( <01L| TE(l,m) |02R> + 866C <02L| TE(l,m) |01R> ) TE(l,m) + ... 867C 868C The factor two is already included in one of the vectors. 869C 870 IF (MZCONF(ISYMV1) .GT. 0 .AND. MZCONF(ISYMV2) .GT. 0) THEN 871C 872C Construct <01L|..|02R> + <02L|..|01R> density 873C 874 ILSYM = MULD2H(IREFSY,ISYMV1) 875 IRSYM = MULD2H(IREFSY,ISYMV2) 876 NCL = MZCONF(ISYMV1) 877 NCR = MZCONF(ISYMV2) 878 KZVARL = MZYVAR(ISYMV1) 879 KZVARR = MZYVAR(ISYMV2) 880 LREF = .FALSE. 881 ISYMDN = MULD2H(ILSYM,IRSYM) 882 CALL DZERO(DEN1,NASHT*NASHT) 883 CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 884 * VEC1,VEC2,OVLAP,DEN1,DUMMY,JSPIN,JSPIN,TDM,NORHO2, 885 * XINDX,WRK,KFREE,LFREE,LREF) 886C 887Clf Call n.10 888Clf 889 IKLVL = 0 890 CALL POPGET(DUMMY,DUMMY,DUMMY,IKLVL,OVLAP,ISYMDN, 891 * IDUMMY,IDUMMY,IDUMMY,TE,DEN1,0, 892 * CMO,WRK(KFREE),LFREE,100,'TE cont3',DM1, 893 * DBGFLG(10)) 894 END IF 895C 896C / <0| [qj ,TE] |0> \ 897C | <j| TE |0> | 898C | <0| [qj+,TE] |0> | 899C \ -<0| TE |j> / 900C 901 ISYMDN = 1 902 OVLAP = D1 903 ISYMV = IREFSY 904 ISYMST = MULD2H(IGRSYM,IREFSY) 905 IF ( ISYMST .EQ. IREFSY ) THEN 906 LCON = ( MZCONF(IGRSYM) .GT. 1 ) 907 ELSE 908 LCON = ( MZCONF(IGRSYM) .GT. 0 ) 909 END IF 910 LORB = ( MZWOPT(IGRSYM) .GT. 0 ) 911 LREF = .TRUE. 912 NZYVEC = MZCONF(1) 913 NZCVEC = MZCONF(1) 914 CALL RSP1GR(NSIM,KZYVR,IDUMMY,JSPIN,IGRSYM,JSPIN,ISYMV,ETRS, 915 * WRK(KCREF),NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,TE, 916 * XINDX,MJWOP,WRK(KFREE),LFREE,LORB,LCON,LREF) 917C 918Clf and consequently we do not need this..... 919CC Restore the vector 920C 921C CALL DSCAL(KZYV1,DH,VEC1,1) 922CC 923Clf 924C 925 CALL PRIRES(ETRS,VECA,IGRSYM,'PCASE4') 926C 927 RETURN 928 END 929 930