1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C============================================================================= 19C /* Deck E3IEF */ 20C============================================================================= 21 SUBROUTINE E3IEF2(VECA, VEC1, VEC2,ETRS,XINDX,ZYM1,ZYM2, 22 * DEN1,UDV,WORK,LWORK,KZYVR,KZYV1,KZYV2, 23 * IGRSYM,ISYMV1,ISYMV2,CMO,MJWOP) 24C 25C 26C Purpose: 27C Outer driver routine for IEF-PCM solvent contribution 28C to E[3] times two vectors. 29C Completely rewritten from scratch! 30C 31 32C now we build the V[3] term of the gradient 33#include "implicit.h" 34#include "dummy.h" 35C 36#include "maxorb.h" 37#include "mxcent.h" 38#include "priunit.h" 39#include "inforb.h" 40#include "infdim.h" 41#include "infinp.h" 42#include "infvar.h" 43#include "infrsp.h" 44#include "infpri.h" 45#include "rspprp.h" 46#include "infcr.h" 47#include "infspi.h" 48#include "infden.h" 49#include "pcmdef.h" 50#include "pcm.h" 51 PARAMETER ( D1 = 1.0D0, DM1 = -1.0D0, D2 = 2.0D0, DM2 = -2.0D0, 52 $ DP5= 0.5D0, DMP5= -0.5D0) 53 54 DIMENSION ETRS(KZYVR),XINDX(*) 55 DIMENSION UDV(NASHDI,NASHDI),DEN1(NASHDI,NASHDI) 56 DIMENSION ZYM1(*),ZYM2(*),WORK(*),CMO(*) 57 DIMENSION VEC1(KZYV1),VEC2(KZYV2),VECA(KZYVR) 58 DIMENSION MJWOP(2,MAXWOP,8) 59 INTEGER ADDFLG 60 61 LOGICAL DYNCHG,LORB,LCON,LREF 62 63 WRITE(LUPRI,*) 'INITIAL ETRS' 64 CALL OUTPUT(ETRS,1,KZYVR,1,1,KZYVR,1,1,LUPRI) 65 66 KFREE = 1 67 LFREE = LWORK 68 CALL MEMGET('REAL',KCHG, NTS, WORK,KFREE,LFREE) 69 CALL MEMGET('REAL',KQ12, NTS, WORK,KFREE,LFREE) 70 CALL MEMGET('REAL',KQ1, NTS, WORK,KFREE,LFREE) 71 CALL MEMGET('REAL',KQ2, NTS, WORK,KFREE,LFREE) 72 CALL MEMGET('REAL',KQD1, NTS, WORK,KFREE,LFREE) 73 CALL MEMGET('REAL',KQD2, NTS, WORK,KFREE,LFREE) 74 CALL MEMGET('REAL',KQ1D2, NTS, WORK,KFREE,LFREE) 75 CALL MEMGET('REAL',KQ2D1, NTS, WORK,KFREE,LFREE) 76 CALL MEMGET('REAL',KQD12, NTS, WORK,KFREE,LFREE) 77 CALL MEMGET('REAL',KCTS,3*NTS, WORK,KFREE,LFREE) 78C 79 CALL DZERO(WORK(KCHG), NTS) 80 CALL DZERO(WORK(KQ12), NTS) 81 CALL DZERO(WORK(KQ1), NTS) 82 CALL DZERO(WORK(KQ2), NTS) 83 CALL DZERO(WORK(KQD1), NTS) 84 CALL DZERO(WORK(KQD2), NTS) 85 CALL DZERO(WORK(KQ1D2), NTS) 86 CALL DZERO(WORK(KQ2D1), NTS) 87 CALL DZERO(WORK(KQD12), NTS) 88 CALL DZERO(WORK(KCTS),3*NTS) 89 90 DO I=1,NTS 91 K = 3 * (I-1) 92 WORK(KCTS + K) = XTSCOR(I) 93 WORK(KCTS + K + 1) = YTSCOR(I) 94 WORK(KCTS + K + 2) = ZTSCOR(I) 95 END DO 96 97C we fetch the transformation vectors 98 NSIM = 1 99 CALL GTZYMT(NSIM,VEC1,KZYV1,ISYMV1,ZYM1,MJWOP) 100 CALL GTZYMT(NSIM,VEC2,KZYV2,ISYMV2,ZYM2,MJWOP) 101 102c ISYMDN is the symmmetry of the density 103 ISYMDN = 1 104 OVLAP = D1 105 DYNCHG = .FALSE. 106 ISPING = ISPINA 107 ISPIN1 = ISPINB 108 ISPIN2 = ISPINC 109 ISPIN3 = 0 110 111 CALL DAXPY(NTS,D1,QSN,1,WORK(KCHG),1) 112 CALL DAXPY(NTS,D1,QSE,1,WORK(KCHG),1) 113C 114C Expectation values of transformed charges 115C 116 IF(ISPIN1.EQ.0) THEN 117 IKLVL = 1 118 CALL RSPASC(IKLVL,NTS,ISPIN1,IDUMMY,IDUMMY,ISYMV1,IDUMMY, 119 $ IDUMMY,ISYMDN,DYNCHG,OVLAP,ZYM1,DUMMY,DUMMY,WORK(KQ1), 120 $ WORK(KCTS),UDV,CMO,WORK(KFREE),LFREE) 121 END IF 122 123 IF(ISPIN2.EQ.0) THEN 124 IKLVL = 1 125 CALL RSPASC(IKLVL,NTS,ISPIN2,IDUMMY,IDUMMY,ISYMV2,IDUMMY, 126 $ IDUMMY,ISYMDN,DYNCHG,OVLAP,ZYM2,DUMMY,DUMMY,WORK(KQ2), 127 $ WORK(KCTS),UDV,CMO,WORK(KFREE),LFREE) 128 END IF 129 130 IF(ISPIN1.EQ.ISPIN2) THEN 131 IKLVL = 2 132 CALL RSPASC(IKLVL,NTS,ISPIN1,ISPIN2,IDUMMY,ISYMV1,ISYMV2, 133 $ IDUMMY,ISYMDN,DYNCHG,OVLAP,ZYM1,ZYM2,DUMMY,WORK(KQ12), 134 $ WORK(KCTS),UDV,CMO,WORK(KFREE),LFREE) 135 IKLVL = 2 136 CALL RSPASC(IKLVL,NTS,ISPIN2,ISPIN1,IDUMMY,ISYMV2,ISYMV1, 137 $ IDUMMY,ISYMDN,DYNCHG,OVLAP,ZYM2,ZYM1,DUMMY,WORK(KQ12), 138 $ WORK(KCTS),UDV,CMO,WORK(KFREE),LFREE) 139 END IF 140 141CLF IF (MZCONF(ISYMV1).GT.0) THEN 142 IF(.FALSE.) THEN 143 NCASE = 1 144 CALL DENGET(NCASE,DEN1) 145 IKLVL = 0 146 CALL RSPASC(IKLVL,NTS,IDUMMY,IDUMMY,IDUMMY,IDUMMY,IDUMMY, 147 $ IDUMMY,ISYMDN,DYNCHG,OVLAP,DUMMY,DUMMY,DUMMY,WORK(KQD1), 148 $ WORK(KCTS),DEN1,CMO,WORK(KFREE),LFREE) 149 END IF 150 151 IF(.FALSE.) THEN 152 WRITE(LUPRI,*) 'TRANSFORMED CHARGES' 153 WRITE(LUPRI,*) 'qsn+qse' 154 CALL OUTPUT(WORK(KCHG),1,NTS,1,1,NTS,1,1,LUPRI) 155 WRITE(LUPRI,*) 'q(1)' 156 CALL OUTPUT(WORK(KQ1),1,NTS,1,1,NTS,1,1,LUPRI) 157 WRITE(LUPRI,*) 'q(2)' 158 CALL OUTPUT(WORK(KQ2),1,NTS,1,1,NTS,1,1,LUPRI) 159 WRITE(LUPRI,*) 'q(12)' 160 CALL OUTPUT(WORK(KQ12),1,NTS,1,1,NTS,1,1,LUPRI) 161 END IF 162 163 IKLVL = 0 164 ADDFLG = 1 165C 166C Transformed potentials contracted with charges 167C 168 CALL RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS,ISPING, 169 $ IDUMMY,IDUMMY, 170 $ IDUMMY,IDUMMY,IDUMMY,IDUMMY,MJWOP, 171 $ ADDFLG,OVLAP,DMP5,DUMMY,DUMMY,DUMMY,UDV,ETRS,WORK(KQ12),CMO, 172 $ XINDX,WORK(KFREE),LFREE) 173 174 IKLVL = 1 175 ADDFLG = 1 176C 177C V(k2)*q(k1) 178C 179 CALL RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS,ISPING, 180 $ ISPIN1,IDUMMY, 181 $ IDUMMY,ISYMV1,ISYMV2,IDUMMY,MJWOP, 182 $ ADDFLG,OVLAP,DM1,ZYM1,DUMMY,DUMMY,UDV,ETRS,WORK(KQ2),CMO, 183 $ XINDX,WORK(KFREE),LFREE) 184C 185C V(k1)*q(k2) 186C 187 CALL RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS,ISPING, 188 $ ISPIN2,IDUMMY, 189 $ IDUMMY,ISYMV2,ISYMV1,IDUMMY,MJWOP, 190 $ ADDFLG,OVLAP,DM1,ZYM2,DUMMY,DUMMY,UDV,ETRS,WORK(KQ1),CMO, 191 $ XINDX,WORK(KFREE),LFREE) 192 193 IKLVL = 2 194 ADDFLG = 1 195 CALL RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS,ISPING, 196 $ ISPIN1,ISPIN2, 197 $ IDUMMY,ISYMV1,ISYMV2,IDUMMY,MJWOP, 198 $ ADDFLG,OVLAP,DMP5,ZYM1,ZYM2,DUMMY,UDV,ETRS,WORK(KCHG),CMO, 199 $ XINDX,WORK(KFREE),LFREE) 200 CALL RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS,ISPING, 201 $ ISPIN2,ISPIN1, 202 $ IDUMMY,ISYMV2,ISYMV1,IDUMMY,MJWOP, 203 $ ADDFLG,OVLAP,DMP5,ZYM2,ZYM1,DUMMY,UDV,ETRS,WORK(KCHG),CMO, 204 $ XINDX,WORK(KFREE),LFREE) 205 206Clf WRITE(LUPRI,*) 'FINAL ETRS' 207Clf CALL OUTPUT(ETRS,1,KZYVR,1,1,KZYVR,1,1,LUPRI) 208 RETURN 209 END 210 211C============================================================================= 212C /* Deck RSPPOT */ 213C============================================================================= 214 SUBROUTINE RSPPOT(NSIM,KZYVR,IGRSYM,ISYMDN,IKLVL,NTS, 215 $ ISPING,ISPIN1,ISPIN2,ISPIN3,ISYMV1,ISYMV2,ISYMV3, 216 $ MJWOP,ADDFLG,OVLAP,DFCTR, 217 $ ZYM1,ZYM2,ZYM3,UDV,ETRS,CHG,CMO,XINDX, 218 $ WORK,LWORK) 219C 220#include "implicit.h" 221#include "inforb.h" 222#include "infvar.h" 223#include "inftap.h" 224#include "priunit.h" 225#include "infrsp.h" 226#include "dummy.h" 227#include "infspi.h" 228C 229c Makes the potential part of the PCM solvent contribution to the 230c quadratic (and cubic) response gradient. 231C 232C 233 INTEGER ADDFLG 234 INTEGER MJWOP(2,MAXWOP,8) 235 LOGICAL FNDLAB,LORB,LCON,LREF 236 DIMENSION UDV(*),XINDX(*) 237 DIMENSION CMO(*),CHG(*),WORK(*),ETRS(*) 238 DIMENSION ZYM1(*),ZYM2(*),ZYM3(*) 239C 240 KFREE = 1 241 LFREE = LWORK 242 CALL MEMGET('REAL',KVAOT, NNBASX, WORK,KFREE,LFREE) 243 CALL MEMGET('REAL',KVMOT, NNORBX, WORK,KFREE,LFREE) 244 CALL MEMGET('REAL',KVMO, N2ORBX, WORK,KFREE,LFREE) 245 CALL MEMGET('REAL',KVMO2, N2ORBX, WORK,KFREE,LFREE) 246 CALL MEMGET('REAL',KVMO3, N2ORBX, WORK,KFREE,LFREE) 247 CALL MEMGET('REAL',KUCMO, NORBT*NBAST,WORK,KFREE,LFREE) 248C CALL MEMGET('REAL',KV3Q0,N2ORBX,WORK,KFREE,LFREE) 249 250 CALL DZERO(WORK(KVAOT),NNBASX) 251 CALL DZERO(WORK(KVMOT),NNORBX) 252 CALL DZERO(WORK(KVMO ),N2ORBX) 253 CALL DZERO(WORK(KVMO2),N2ORBX) 254 CALL DZERO(WORK(KVMO3),N2ORBX) 255 CALL DZERO(WORK(KUCMO),NORBT*NBAST) 256 257 CALL UPKCMO(CMO,WORK(KUCMO)) 258 259#ifdef PCM_DEBUG 260 print *,'beg of rsppot',iklvl,isymv1,isymv2,isymdn 261#endif 262 263Clf this initialization of symmetry needs to be rethought if one is to implement CR!!!! 264 ISYM = 1 265 IF (IKLVL.EQ.1) ISYM = ISYMV2 266 IF (IKLVL.EQ.2) ISYM = 1 267c IF (IKLVL.EQ.3) ISYM = MULD2H(MULD2H(ISYMV1,ISYMV2),ISYMV3) 268 ISYMT = MULD2H(ISYM,ISYMDN) 269C ISYMT = 1 270 271c 272C Calculate potentials in AO basis, mult. by the charges, transform to MO and unpack 273C 274 275Clf: NOTE: probably nosim needs to be changed according to symmetry. 276 NOSIM = 1 277 278 279#ifdef PCM_DEBUG 280 write(lupri,*) 'chg before pot',iklvl 281#endif 282 CALL OUTPUT(CHG,1,NTS,1,1,NTS,1,1,LUPRI) 283 CALL J1INT(CHG,.FALSE.,WORK(KVAOT),NOSIM,.FALSE.,'NPETES ', 284 & ISYMT,WORK(KFREE),LFREE) 285 CALL UTHU(WORK(KVAOT),WORK(KVMOT),WORK(KUCMO),WORK(KFREE), 286 $ NBAST,NORBT) 287 CALL DSPTSI(NORBT,WORK(KVMOT),WORK(KVMO2)) 288 289#ifdef PCM_DEBUG 290 write(lupri,*) 'potentialsXcharges', iklvl 291#endif 292 CALL OUTPUT(WORK(KVMO2),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 293 294 IGRSPI = 0 295C First transformation of potentials: V^e_{ab} --> V^e_{ab}({}^1\kappa) 296 IF (IKLVL.GE.1) THEN 297 IGRSPI = ISPIN1 298 CALL DZERO(WORK(KVMO),N2ORBX) 299 CALL OITH1(ISYMV1,ZYM1,WORK(KVMO2),WORK(KVMO),ISYM) 300 ISYM = ISYMV1 301 END IF 302C Second transformation of potentials: V^e_{ab}({}^1\kappa) --> V^e_{ab}({}^1\kappa {}^2\kappa) 303 IF (IKLVL.GE.2) THEN 304 IGRSPI = MULSP(IGRSPI,ISPIN2) 305 CALL DZERO(WORK(KVMO2),N2ORBX) 306 CALL OITH1(ISYMV2,ZYM2,WORK(KVMO),WORK(KVMO2),ISYM) 307 ISYM = MULD2H(ISYM,ISYMV2) 308 END IF 309C Third transformation of potentials..... 310 IF (IKLVL.GE.3) THEN 311 IGRSPI = MULSP(IGRSPI,ISPIN3) 312 CALL DZERO(WORK(KVMO),N2ORBX) 313 CALL OITH1(ISYMV3,ZYM3,WORK(KVMO2),WORK(KVMO),ISYM) 314 ISYM = MULD2H(ISYM,ISYMV3) 315 END IF 316 IF ((IKLVL.EQ.0).OR.(IKLVL.EQ.2)) THEN 317 CALL DCOPY(N2ORBX,WORK(KVMO2),1,WORK(KVMO),1) 318 END IF 319#ifdef PCM_DEBUG 320 write(lupri,*) 'potentialsXcharges after transformations', iklvl 321#endif 322 CALL OUTPUT(WORK(KVMO),1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 323 324C 325C ADD A SCALING FACTOR AND SUM THE OPERATOR (ADDFLG.LE.0 ONLY WHEN DEBUGGING) 326C 327 CALL DZERO(WORK(KVMO2),N2ORBX) 328 IF(ADDFLG.GT.0) THEN 329 CALL DAXPY(N2ORBX,DFCTR,WORK(KVMO),1,WORK(KVMO2),1) 330 END IF 331 332 333 ISYMV = IREFSY 334C ISPIN = 0 (ISPING) 335C THIS MUST BE SET HERE!!!!!!!!!!!!!!!! 336C IGRSPI = 0 337 NZYVEC = 0 338 NZCVEC = 0 339 LORB = .TRUE. 340 LCON = .FALSE. 341 LREF = .TRUE. 342 343 IF (.true.) THEN 344 WRITE(LUPRI,*) 'Norm of TOPGET ',DNRM2(N2ORBX,work(kvmo2),1) 345 END IF 346 347 IF(ADDFLG.GT.0) THEN 348 CALL PCM1GR(NSIM,KZYVR,IDUMMY,ISPING,IGRSYM,IGRSPI,ISYMV,ETRS, 349 $ DUMMY,NZYVEC,NZCVEC,OVLAP,ISYMDN,UDV,WORK(KVMO2),XINDX, 350 $ MJWOP,WORK(KFREE),LFREE,LORB,LCON,LREF) 351Clf WRITE(LUPRI,*) 'ETRS after the gradient' 352Clf CALL OUTPUT(ETRS,1,KZYVR,1,1,KZYVR,1,1,LUPRI) 353 END IF 354 355 RETURN 356 END 357 358C============================================================================= 359C /* Deck RSPASC */ 360C============================================================================= 361 SUBROUTINE RSPASC(IKLVL,NTS,ISPIN1,ISPIN2,ISPIN3,ISYMV1, 362 $ ISYMV2,ISYMV3,ISYMDN,DYNCHG,OVLAP,ZYM1,ZYM2,ZYM3,TCHG, 363 $ CTS,DEN,CMO,WORK,LWORK) 364C 365#include "implicit.h" 366#include "mxcent.h" 367#include "inforb.h" 368#include "inftap.h" 369#include "infrsp.h" 370#include "priunit.h" 371#include "orgcom.h" 372C 373C Makes the charge part of the PCM solvent contribution of the 374C quadraticn and (possibly) the cubic response. The subroutine 375C takes info about the n. of one-index transformations, spin and 376C symmetry. The required density is constructed inside. 377C 378C 379C 380 PARAMETER (D1=1.0D0) 381 LOGICAL FNDLAB 382 LOGICAL DYNCHG 383 LOGICAL TRIMAT,EXP1VL,TOFILE 384 DIMENSION CMO(*),DEN(*),TCHG(*),CTS(*),WORK(*) 385 DIMENSION ZYM1(*),ZYM2(*),ZYM3(*) 386 DIMENSION INTREP(9*MXCENT),INTADR(9*MXCENT) 387 CHARACTER*8 LABINT(9*MXCENT) 388 389 KFREE = 1 390 LFREE = LWORK 391Clf dirty fix!! 392 NTSIRR = NTS 393#ifdef PCM_DEBUG 394 print *,'rspasc',lwork,lfree,nnbasx,norbt,nbast,n2orbx,nts,ntsirr 395#endif 396 CALL MEMGET('REAL',KCHGAO,NNBASX, WORK,KFREE,LFREE) 397 CALL MEMGET('REAL',KUCMO ,NORBT*NBAST,WORK,KFREE,LFREE) 398 CALL MEMGET('REAL',KTELM ,N2ORBX, WORK,KFREE,LFREE) 399 CALL MEMGET('REAL',KTLMA ,N2ORBX, WORK,KFREE,LFREE) 400 CALL MEMGET('REAL',KTCHG ,NTS, WORK,KFREE,LFREE) 401 CALL MEMGET('REAL',KQCHG ,NTS, WORK,KFREE,LFREE) 402 CALL MEMGET('REAL',KQMAT ,NTS*NTSIRR, WORK,KFREE,LFREE) 403 404 CALL DZERO(WORK(KCHGAO),NNBASX) 405 CALL DZERO(WORK(KUCMO) ,NORBT*NBAST) 406 CALL DZERO(WORK(KTELM) ,N2ORBX) 407 CALL DZERO(WORK(KTLMA) ,N2ORBX) 408 CALL DZERO(WORK(KTCHG) ,NTS) 409 CALL DZERO(WORK(KQCHG) ,NTS) 410 CALL DZERO(WORK(KQMAT) ,NTS*NTSIRR) 411 412C CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ', 413C * 'UNFORMATTED',IDUMMY,.FALSE.) 414 415C 416C Unpack symmetry blocked CMO 417C 418 CALL UPKCMO(CMO,WORK(KUCMO)) 419 420 421 IF (IKLVL.LE.0) ISYM = 1 422 IF (IKLVL.EQ.1) ISYM = ISYMV1 423 IF (IKLVL.EQ.2) ISYM = MULD2H(ISYMV1,ISYMV2) 424 IF (IKLVL.EQ.3) ISYM = MULD2H(MULD2H(ISYMV1,ISYMV2),ISYMV3) 425 ISYMT = MULD2H(ISYM,ISYMDN) 426 427C REWIND(LUPROP) 428 DO ITS=1,NTS 429 ISYM = ISYMT 430 CALL DZERO(WORK(KCHGAO),NNBASX) 431C IF (DYNCHG) THEN 432C IF (FNDLAB('J3-PCMIN',LUPROP)) THEN 433C CALL READT(LUPROP,NNBASX,WORK(KCHGAO)) 434C ELSE 435C WRITE (LUPRI,'(/A)') ' Integral label J3-PCMIN not'// 436C & 'found on file AOPROPER' 437C CALL QUIT('Integral label not found in POPGET') 438C END IF 439C ELSE 440C IF (FNDLAB('J2-PCMIN',LUPROP)) THEN 441C CALL READT(LUPROP,NNBASX,WORK(KCHGAO)) 442C ELSE 443C WRITE (LUPRI,'(/A)') ' Integral label J2-PCMIN not'// 444C & 'found on file AOPROPER' 445C CALL QUIT('Integral label not found in POPGET') 446C END IF 447C END IF 448 L = 1 449 NCOMP = NSYM 450 KTS = 3*(ITS-1) 451 DIPORG(1) = CTS(KTS+1) 452 DIPORG(2) = CTS(KTS+2) 453 DIPORG(3) = CTS(KTS+3) 454 EXP1VL = .FALSE. 455 TOFILE = .FALSE. 456 KPATOM = 0 457 TRIMAT = .TRUE. 458 CALL GET1IN(WORK(KCHGAO),'NPETES ',NCOMP,WORK(KFREE),LWORK, 459 & LABINT,INTREP,INTADR,L,TOFILE,KPATOM,TRIMAT, 460 & DUMMY,EXP1VL,DUMMY,IPRRSP) 461 JCHGAO = KCHGAO 462 DO ILOP = 1, NSYM 463 ISYM = ILOP 464 JTS = (ILOP - 1)*NTSIRR + ITS 465 CALL UTHU(WORK(JCHGAO),WORK(KTLMA),WORK(KUCMO),WORK(KFREE), 466 $ NBAST,NORBT) 467 CALL DSPTSI(NORBT,WORK(KTLMA),WORK(KTELM)) 468C First transformation of charges: q^e_{ab} --> q^e_{ab}({}^1\kappa) 469 IF (IKLVL.GE.1) THEN 470 CALL DZERO(WORK(KTLMA),N2ORBX) 471 CALL OITH1(ISYMV1,ZYM1,WORK(KTELM),WORK(KTLMA),ISYM) 472 ISYM = MULD2H(ISYM,ISYMV1) 473 END IF 474C Second transformation of charges: q^e_{ab}({}^1\kappa) --> q^e_{ab}({}^1\kappa {}^2\kappa) 475 IF (IKLVL.GE.2) THEN 476 CALL DZERO(WORK(KTELM),N2ORBX) 477 CALL OITH1(ISYMV2,ZYM2,WORK(KTLMA),WORK(KTELM),ISYM) 478 ISYM = MULD2H(ISYM,ISYMV2) 479 END IF 480C Third transformation of charges: hope you can figure out the formula..... 481 IF (IKLVL.GE.3) THEN 482 CALL DZERO(WORK(KTLMA),N2ORBX) 483 CALL OITH1(ISYMV3,ZYM3,WORK(KTELM),WORK(KTLMA),ISYM) 484 ISYM = MULD2H(ISYM,ISYMV3) 485 END IF 486 IF ((IKLVL.EQ.1).OR.(IKLVL.EQ.3)) THEN 487 CALL DCOPY(N2ORBX,WORK(KTLMA),1,WORK(KTELM),1) 488 END IF 489C Contract transformed charges with the density 490 CALL MELONE(WORK(KTELM),ISYM,DEN,OVLAP,FACT,200,'RSPASC') 491#ifdef PCM_DEBUG 492 print *,its,fact,ilop,jts 493#endif 494 WORK(KTCHG + JTS - 1) = FACT 495 JCHGAO = JCHGAO + NNBASX 496 END DO 497 END DO 498 499 CALL V2Q(WORK(KQMAT),WORK(KTCHG),WORK(KQCHG),QTEXS,.false.) 500Clf CALL GPCLOSE(LUPCMD,'KEEP') 501 CALL DAXPY(NTS,D1,WORK(KQCHG),1,TCHG,1) 502 503 RETURN 504 END 505 506C============================================================================= 507C /* Deck PCMGR1 */ 508C============================================================================= 509 SUBROUTINE PCM1GR(NSIM,KZYVR,INTSYM,ISPIN,IGRSYM,IGRSPI, 510 * ISYMV,OTRS, 511 * VEC,NZYVEC,NZCVEC,OVLAP,ISYMDN,DEN1,OPMAT, 512 * XINDX,MJWOP,WORK,LWORK,LORB,LCON,LREFST) 513C 514C Compute the gradient resulting from multiplying the S matrix 515C with one ore more vectors. Copied from RSP1GR 516C 517#include "implicit.h" 518#include "infdim.h" 519#include "inforb.h" 520#include "priunit.h" 521#include "infrsp.h" 522#include "wrkrsp.h" 523#include "rspprp.h" 524#include "infhyp.h" 525#include "infvar.h" 526#include "qrinf.h" 527#include "infpri.h" 528C 529 LOGICAL LORB, LCON, LREFST 530C 531 DIMENSION OTRS (KZYVR), MJWOP(2,MAXWOP,8) 532 DIMENSION VEC(NZYVEC) 533 DIMENSION DEN1(NASHDI,NASHDI) 534 DIMENSION OPMAT(NORBT,NORBT) 535 DIMENSION XINDX(*) 536 DIMENSION WORK(*) 537C 538#ifdef PCM_DEBUG 539 print *,'lorb,lcon,lrefst',lorb,lcon,lrefst 540#endif 541clf IF ( IPRRSP .GT. 150 ) THEN 542 IF ( .true. ) THEN 543 WRITE(LUPRI,'(A)') ' Vector in PCM1GR' 544 IF (LREFST) WRITE(LUPRI,'(A)') ' (Reference state)' 545 CALL OUTPUT(VEC,1,NZYVEC,1,NSIM,NZYVEC,NSIM,1,LUPRI) 546 IF ( LORB ) THEN 547 WRITE(LUPRI,'(//A)') ' Density matrix in PCP1GR' 548 CALL OUTPUT(DEN1,1,NASHT,1,NASHT,NASHT,NASHT,1,LUPRI) 549 END IF 550 WRITE(LUPRI,'(//A)') ' One electron matrix in PCM1GR' 551 CALL OUTPUT(OPMAT,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 552 END IF 553C 554 IF ( LORB ) THEN 555 TRPLET = IGRSPI.NE.ISPIN 556 CALL ORBSX(NSIM,IGRSYM,KZYVR,OTRS,OPMAT,OVLAP, 557 * ISYMDN,DEN1,MJWOP,WORK,LWORK) 558 END IF 559C 560 IF ( LCON ) THEN 561 ISYMJ = MULD2H( IGRSYM, IREFSY ) 562 NZCSTJ = MZCONF( IGRSYM ) 563C 564 TRPLET = ISPIN.EQ.1 565 CALL CONSX(NSIM,KZYVR,IGRSYM,OPMAT,VEC,NZYVEC, 566 * NZCVEC,ISYMV,NZCSTJ,ISYMJ,LREFST,OTRS,XINDX, 567 * WORK,LWORK) 568 END IF 569C 570 RETURN 571 END 572 573C============================================================================= 574C /* Deck DENGET */ 575C============================================================================= 576 SUBROUTINE DENGET(NCASE,DEN1) 577#include "implicit.h" 578#include "dummy.h" 579 LOGICAL LREF 580C 581C L. Frediani, Nov 2005. Purpose: initialization of parameters for 582C different densities needed for singlet and triplet HF and MCSCF 583C quadratic response calculations 584C 585 586C 587C parameters initialization 588C 589 GO TO (101,102,103,104,105,106,107,108,109,110), NCASE 590 591 101 CONTINUE 592c ILSYM = MULD2H(IREFSY,ISYMV1) 593c IRSYM = MULD2H(IREFSY,ISYMV2) 594c NCL = MZCONF(ISYMV1) 595c NCR = MZCONF(ISYMV2) 596c KZVARL = MZYVAR(ISYMV1) 597c KZVARR = MZYVAR(ISYMV2) 598c LREF = .FALSE. 599cc ISYMDN = MULD2H(ILSYM,IRSYM) 600cC 601cC We get triplet density in case op. A is triplet 602cC 603c JSPIN1 = MULSP(ISPIN1,ISPIN2) 604c JSPIN2 = 0 605 606 GOTO 200 607 102 CONTINUE 608 GOTO 200 609 103 CONTINUE 610 GOTO 200 611 104 CONTINUE 612 GOTO 200 613 105 CONTINUE 614 GOTO 200 615 106 CONTINUE 616 GOTO 200 617 107 CONTINUE 618 GOTO 200 619 108 CONTINUE 620 GOTO 200 621 109 CONTINUE 622 GOTO 200 623 110 CONTINUE 624 GOTO 200 625C 626C Density calculation 627C 628 200 CONTINUE 629 630c CALL DZERO(DEN1,NASHT*NASHT) 631c CALL RSPGDM(NSIM,ILSYM,IRSYM,NCL,NCR,KZVARL,KZVARR, 632c * VEC1,VEC2,OVLAP,DEN1,DUMMY,JSPIN1,JSPIN2,TDM,NORHO2, 633c * XINDX,WRK,KFREE,LFREE,LREF) 634 RETURN 635 END 636