1C /* Deck pcmltr */ 2 SUBROUTINE PCMLTR(NCSIM,NOSIM,BCVECS,BOVECS,CREF,CMO,INDXCI, 3 & UDV,DV,UDVTR,DVTR,DTV,DTVTR,SCVECS,SOVECS, 4 & WRK,LWRK) 5C 6C LF BM RC 24 oct 01 based on PCMLIN 7C 8C Common driver for PCMLNC and PCMLNO 9C 10#include "implicit.h" 11#include "dummy.h" 12 13 DIMENSION BCVECS(*),BOVECS(*),CREF(*),CMO(*),INDXCI(*) 14 DIMENSION UDV(*),DV(*),DTV(*),SCVECS(*),SOVECS(*),WRK(LWRK) 15 DIMENSION UDVTR(*),DVTR(*),DTVTR(*) 16C 17C 18#include "maxorb.h" 19#include "priunit.h" 20#include "infpri.h" 21#include "infrsp.h" 22#include "inftap.h" 23#include "wrkrsp.h" 24C 25 CALL QENTER('PCMLTR') 26C 27 CALL GPOPEN(LUPROP,'AOPROPER','OLD',' ', 28 * 'UNFORMATTED',IDUMMY,.FALSE.) 29C if SOPRPA, then only put PCM contribs into RPA part 30C IF (NCSIM .GT. 0) THEN 31 IF ((NCSIM .GT. 0) .and. (.not. SOPRPA)) THEN 32 IF (IPRRSP.GT.101) THEN 33C IF (.true.) THEN 34 WRITE(LUPRI,*)' LINEAR TRANSFORMED CONFIGURATION VECTOR' 35 WRITE(LUPRI,*)' **** BEFORE IEFLNC **** ' 36 CALL OUTPUT(SCVECS,1,KZYVAR,1,NCSIM,KZYVAR,NCSIM,1,LUPRI) 37 END IF 38 CALL IEFLNC(NCSIM,BCVECS,CREF,CMO,INDXCI, 39 * UDV,DV,UDVTR,DVTR,DTV,DTVTR,SCVECS,WRK,LWRK) 40 IF (IPRRSP.GT.101) THEN 41C IF (.true.) THEN 42 WRITE(LUPRI,*)' LINEAR TRANSFORMED CONFIGURATION VECTOR' 43 WRITE(LUPRI,*)' **** AFTER IEFLNC **** ' 44 CALL OUTPUT(SCVECS,1,KZYVAR,1,NCSIM,KZYVAR,NCSIM,1,LUPRI) 45 END IF 46 END IF 47 IF ( NOSIM .GT.0 ) THEN 48 IF (IPRRSP.GT.101) THEN 49 WRITE(LUPRI,*)' LINEAR TRANSFORMED ORBITAL VECTOR' 50 WRITE(LUPRI,*)' **** BEFORE IEFLNO **** ',kzyvar,nosim 51 CALL OUTPUT(SOVECS,1,KZYVAR,1,NOSIM,KZYVAR,NOSIM,1,LUPRI) 52 END IF 53 CALL IEFLNO(NCSIM,NOSIM,BOVECS,CREF,CMO,INDXCI, 54 * UDV,DV,UDVTR,DVTR,SOVECS,WRK,LWRK) 55 IF (IPRRSP.GT.101) THEN 56 WRITE(LUPRI,*)' LINEAR TRANSFORMED ORBITAL VECTOR' 57 WRITE(LUPRI,*)' **** AFTER IEFLNO **** ' 58 CALL OUTPUT(SOVECS,1,KZYVAR,1,NOSIM,KZYVAR,NOSIM,1,LUPRI) 59 END IF 60 END IF 61 IF (LUPROP .GT. 0) CALL GPCLOSE(LUPROP,'KEEP') 62 CALL QEXIT('PCMLTR') 63 RETURN 64 END 65C /* Deck ieflnc */ 66 SUBROUTINE IEFLNC(NCSIM,BCVEC,CREF,CMO,INDXCI, 67 * UDV,DV,UDVTR,DVTR,DTV,DTVTR,SVEC, 68 * WRK,LFREE) 69C 70C 24 OCT 01 71C 72C Purpose: Calculate MCSCF E2 contribution from a 73C pcm-ief surrounding medium to a csf trial vector. 74C 75#include "implicit.h" 76 DIMENSION BCVEC(*), CREF(*), CMO(*) 77 DIMENSION INDXCI(*), UDV(*), DV(*), DTV(N2ASHX,*) 78 DIMENSION SVEC(KZYVAR,*), WRK(*) 79 DIMENSION UDVTR(*),DVTR(*),DTVTR(N2ASHX,*) 80C 81#include "iratdef.h" 82C 83 PARAMETER ( D0 = 0.0D0 , D2 = 2.0D0 , THRZER = 1.0D-14 ) 84 LOGICAL FNDLAB 85C 86C 87C Used from common blocks: 88C INFINP : ? 89C INFORB : NNASHX, NNORBX, NNBASX, etc. 90C INFTAP : LUSOL, LBSYMB 91C 92#include "maxash.h" 93#include "maxorb.h" 94#include "mxcent.h" 95#include "pcmdef.h" 96#include "priunit.h" 97#include "orgcom.h" 98#include "infinp.h" 99#include "inforb.h" 100#include "infrsp.h" 101#include "wrkrsp.h" 102#include "inftap.h" 103#include "infpri.h" 104#include "pcm.h" 105#include "pcmlog.h" 106 107C 108Clf this is wrong!!! Must be fixed!!!!!11 109 DIMENSION VTEX(NTS,NCSIM) 110C 111 CALL QENTER('IEFLNC') 112C 113C Core allocation 114C 115 KUCMO = 1 116Clara define KJPXAO 117 KJPXAO = KUCMO + NORBT*NBAST 118 KJPX = KJPXAO + NNBASX 119Clara 120C KJPX = KUCMO + NORBT*NBAST 121 KJPXAC = KJPX + NNORBX 122 KPCMXB = KJPXAC + NNASHX 123 KXBAC = KPCMXB + NCSIM*NNORBX 124 KEXBAC = KXBAC + NCSIM*NNASHX 125 KW10 = KEXBAC + NCSIM 126C 127 KCHRG = KW10 128 KW20 = KCHRG + NTS 129C 2.1 read rlmao in ao basis and transform to rlm in mo basis 130 KDV = KW20 131Clara add definition of KDW 132 KDW = KDV + NCSIM*NNBASX 133 KW30 = KDW + NCSIM*N2BASX 134C 135 KJPCMAO = KW10 136 KW50 = KJPCMAO + NNBASX 137 LW30 = LFREE + 1 - KW30 138C 139C 4.0 rspsor 140 KURXC = KW10 + NNORBX 141 KURYC = KURXC + N2ORBX 142 KW40 = KURYC + N2ORBX 143 KNEED = MAX(KW20,KW30) 144 KNEED = MAX(KNEED,KW40) 145 IF (KNEED .GT. LFREE) CALL ERRWRK('IEFLNC',KNEED,LFREE) 146C 147C 148C 149 CALL DZERO(WRK(KJPX),NNORBX) 150Clara clear space for charges, jpxao 151 CALL DZERO(WRK(KCHRG),NTS) 152 CALL DZERO(WRK(KJPXAO),NNBASX) 153C 154C Unpack symmetry blocked CMO 155C 156 CALL UPKCMO(CMO,WRK(KUCMO)) 157C 158CLF Read Potential integrals on tesserae 159C 160 CALL DCOPY(NTS,QSN,1,WRK(KCHRG),1) 161 CALL DAXPY(NTS,1.0D0,QSE,1,WRK(KCHRG),1) 162Clara multiply by -1 the tot charges 163 CALL DSCAL(NTS,-1.0D0,WRK(KCHRG),1) 164C write (lupri,*),'qn' 165C call output(qsn,1,nts,1,1,nts,1,1,lupri) 166C write (lupri,*),'qe' 167C call output(qse,1,nts,1,1,nts,1,1,lupri) 168 CALL J1INT(WRK(KCHRG),.FALSE.,WRK(KJPXAO),1,.FALSE.,'NPETES ', 169 & 1,WRK(KW30),LW30) 170 CALL UTHU(WRK(KJPXAO),WRK(KJPX),WRK(KUCMO),WRK(KW30),NBAST,NORBT) 171C 172 IF (NASHT .GT. 0) THEN 173 CALL GETAC2(WRK(KJPX),WRK(KJPXAC)) 174 END IF 175 IF (IPRRSP .GE. 15) THEN 176 WRITE (LUPRI,'(/A)') ' J+X_mo matrix:' 177 CALL OUTPAK(WRK(KJPX), NORBT,1,LUPRI) 178 IF (NASHT .GT. 0) THEN 179 WRITE (LUPRI,'(/A)') ' J+X_ac matrix:' 180 CALL OUTPAK(WRK(KJPXAC),NASHT,1,LUPRI) 181 END IF 182 END IF 183C 184C Expectation value of J+X 185C 186 EXPJPX = SOLELM(DV,WRK(KJPXAC),WRK(KJPX),ACJPX) 187 IF (IPRRSP .GE. 6) THEN 188 WRITE (LUPRI,'(A,F17.8)') 189 * ' --- J+X expectation value :',EXPJPX 190 WRITE (LUPRI,'(A,F17.8)') 191 * ' --- active part of J1X :',ACJPX 192 END IF 193C 194 DO ICSIM = 1, NCSIM 195Clara first clear for KDV, KDW and charges 196 CALL DZERO(WRK(KCHRG),NTS) 197 CALL DZERO(WRK(KDW),N2BASX) 198 CALL DZERO(WRK(KDV),NNBASX) 199Clara 200#ifdef PCM_DEBUG 201 write (lupri,*) 'lara dtv=' 202 call output(dtv(1,icsim),1,nasht,1,nasht,nasht,nasht,1,lupri) 203#endif 204clara 205Clara change call to fckden with fckden2 and KDV with KDW 206 CALL FCKDEN2(.FALSE.,.TRUE.,DUMMY,WRK(KDW),CMO,DTV(1,ICSIM), 207 & WRK(KW30),LW30) 208Clara add call to dgfsp 209Clara 210#ifdef PCM_DEBUG 211 write (lupri,*) 'lara kdw=' 212 call output(wrk(kdw),1,nbast,1,nbast,nbast,nbast,1,lupri) 213#endif 214clara 215 CALL DGEFSP(NBAST,WRK(KDW),WRK(KDV)) 216clara 217#ifdef PCM_DEBUG 218 write (lupri,*) 'lara kdv=' 219 call outpak (wrk(kdv),nbast,1,lupri) 220 write (lupri,*) 'lara kdv after=' 221 call outpak (wrk(kdv),nbast,1,lupri) 222#endif 223clara 224 CALL J1INT(WRK(KCHRG),.TRUE.,WRK(KDV),1,.FALSE.,'NPETES ', 225 & KSYMOP,WRK(KW30),LW30) 226 CALL DSCAL(NTS,-1.0D0,WRK(KCHRG),1) 227Clara 228#ifdef PCM_DEBUG 229 write (lupri,*) 'lara kchrg=' 230 call output (wrk(kchrg),1,nts,1,1,nts,1,1,lupri) 231#endif 232clara 233 CALL V2Q(WRK(KW30),WRK(KCHRG),VTEX(1,ICSIM),QET,NEQRSP) 234 IF (KSYMOP .NE. 1) THEN 235 CALL DCOPY(NTSIRR,VTEX((KSYMOP - 1)*NTSIRR + 1,ICSIM),1, 236 & VTEX(1,ICSIM),1) 237 CALL DZERO(VTEX(NTSIRR+1,ICSIM),NTS-NTSIRR) 238 END IF 239 END DO 240 CALL GPCLOSE(LUPCMD,'KEEP') 241C 242 CALL DZERO(WRK(KPCMXB),NCSIM*NNORBX) 243C 244 DO ICSIM = 1, NCSIM 245 CALL DZERO(WRK(KJPCMAO),NNBASX) 246Clara 247#ifdef PCM_DEBUG 248 write (lupri,*) 'lara vtex=' 249 call output(vtex,1,nts,1,ncsim,nts,ncsim,1,lupri) 250C write (lupri,'(<1+(nts-1)/4>(4f18.12/))') vtex(:,icsim) 251#endif 252clara 253 CALL J1INT(VTEX(1,ICSIM),.FALSE.,WRK(KJPCMAO),1,.FALSE., 254 & 'NPETES ',KSYMOP,WRK(KW30),LW30) 255 JPCMXB = KPCMXB + (ICSIM - 1)*NNORBX 256Clara 257#ifdef PCM_DEBUG 258 write (lupri,*) 'lara icsim=', icsim,' nnbasx=', nnbasx, 'jpcmao=' 259 call outpak(wrk(kjpcmao),nbast,1,lupri) 260#endif 261clara 262 CALL UTHU(WRK(KJPCMAO),WRK(JPCMXB),WRK(KUCMO),WRK(KW30), 263 & NBAST,NORBT) 264 END DO 265C 266Clara clear for kxbac 267 CALL DZERO(WRK(KXBAC),NNASHX*NCSIM) 268 DO ICSIM = 1, NCSIM 269 JPCMXB = KPCMXB + (ICSIM - 1)*NNORBX 270 JXBAC = KXBAC + (ICSIM - 1)*NNASHX 271 CALL GETAC2(WRK(JPCMXB),WRK(JXBAC)) 272Clara 273#ifdef PCM_DEBUG 274 write (lupri,*) 'lara icsim=',icsim,' nnashx=',nnashx, 'jxbac=' 275 call outpak(wrk(jxbac),nbast,1,lupri) 276#endif 277clara 278 IF (TRPLET) THEN 279 TEMP = SOLELM(DVTR,WRK(JXBAC),WRK(JPCMXB),XBAC) 280 TEMP = XBAC 281 ELSE 282 TEMP = SOLELM(DV,WRK(JXBAC),WRK(JPCMXB),XBAC) 283 ENDIF 284 WRK(KEXBAC - 1 + ICSIM) = XBAC 285 END DO 286C 287 IF (IPRRSP.GT.101) THEN 288 WRITE(LUPRI,*)' LINEAR TRANSFORMED CONFIGURATION VECTOR' 289 WRITE(LUPRI,*)' **** BEFORE SLVSC in IEFLNC **** ' 290 CALL OUTPUT(SVEC,1,KZYVAR,1,NCSIM,KZYVAR,NCSIM,1,LUPRI) 291 END IF 292C 293 CALL SLVSC(NCSIM,0,NNASHX,BCVEC,CREF,SVEC,WRK(KXBAC),WRK(KJPXAC), 294 * WRK(KEXBAC),ACJPX,INDXCI,WRK(KW30),LW30) 295C CALL SLVSC(NCSIM,NOSIM,NNASHX,BCVECS,CREF,SVECS, 296C * RXAC,RYAC,TRXAC,TRYAC,INDXCI,WRK,LWRK) 297C 298 IF (IPRRSP.GT.101) THEN 299 WRITE(LUPRI,*)' LINEAR TRANSFORMED CONFIGURATION VECTOR' 300 WRITE(LUPRI,*)' **** AFTER SLVSC in IEFLNC **** ' 301 CALL OUTPUT(SVEC,1,KZYVAR,1,NCSIM,KZYVAR,NCSIM,1,LUPRI) 302 END IF 303C 304C ... orbital part of sigma vector(s) 305C 306 IF (KZWOPT .GT. 0) THEN 307 JPCMXB = KPCMXB 308C 309 CALL DSPTSI(NORBT,WRK(KJPX),WRK(KURYC)) 310 DO 800 ICSIM = 1,NCSIM 311C 312 CALL DSPTSI(NORBT,WRK(JPCMXB),WRK(KURXC)) 313 IF (TRPLET) THEN 314 CALL SLVSOR(.TRUE.,.FALSE.,1,UDVTR, 315 * SVEC(1,ICSIM),WRK(KURXC)) 316 ELSE 317 CALL SLVSOR(.TRUE.,.TRUE.,1,UDV, 318 * SVEC(1,ICSIM),WRK(KURXC)) 319 END IF 320 IF (IPRRSP.GT.101) THEN 321 WRITE(LUPRI,*)' **** AFTER SLVSOR in IEFLNC **** ' 322 WRITE(LUPRI,*) 323 * ' orbital part of LINEAR TRANSFORMED CONF VEC No',ICSIM 324 CALL OUTPUT(SVEC(1,ICSIM),1,KZYVAR,1,1,KZYVAR,1,1,LUPRI) 325 END IF 326 IF (TRPLET) THEN 327 CALL SLVSOR(.FALSE.,.FALSE.,1,DTVTR(1,ICSIM), 328 * SVEC(1,ICSIM),WRK(KURYC)) 329 ELSE 330 CALL SLVSOR(.FALSE.,.FALSE.,1,DTV(1,ICSIM), 331 * SVEC(1,ICSIM),WRK(KURYC)) 332 END IF 333 IF (IPRRSP.GT.101) THEN 334 WRITE(LUPRI,*) 335 * ' orbital part of LINEAR TRANSFORMED CONF VEC No',ICSIM 336 CALL OUTPUT(SVEC(1,ICSIM),1,KZYVAR,1,1,KZYVAR,1,1,LUPRI) 337 END IF 338 339 JPCMXB = JPCMXB + NNORBX 340 800 CONTINUE 341 IF (IPRRSP.GT.101) THEN 342 WRITE(LUPRI,*)' LINEAR TRANSFIRMED CONFIGURATION VECTOR' 343 WRITE(LUPRI,*)' **** AFTER SLVSOR in IEFLNC **** ' 344 CALL OUTPUT(SVEC,1,KZYVAR,1,NCSIM,KZYVAR,NCSIM,1,LUPRI) 345 END IF 346 END IF 347C 348 CALL QEXIT('IEFLNC') 349 RETURN 350C end of IEFLNC. 351 END 352C /* Deck ieflno */ 353 SUBROUTINE IEFLNO(NCSIM,NOSIM,BOVECS,CREF,CMO,INDXCI, 354 * UDV,DV,UDVTR,DVTR,SVEC, 355 * WRK,LFREE) 356C 357C 358C Purpose: Calculate MCSCF E2 contribution from a 359C surrounding ief-pcm medium to an orbital trial vector. 360C 361C 362#include "implicit.h" 363 DIMENSION BOVECS(*), CREF(*), CMO(*) 364 DIMENSION INDXCI(*), UDV(*), DV(*) 365 DIMENSION SVEC(KZYVAR,*), WRK(*) 366 DIMENSION UDVTR(*),DVTR(*) 367C 368#include "iratdef.h" 369C 370 PARAMETER ( D0 = 0.0D0 , D1 = 1.0D0, D2 = 2.0D0, DP5 = 0.50D0) 371 LOGICAL FNDLAB, TOFILE, EXP1VL, TRIMAT 372#include "dummy.h" 373C 374C 375C Used from common blocks: 376C INFINP : ? 377C INFORB : NNASHX, NNORBX, NNBASX, etc. 378C INFVAR : JWOP 379C INFRSP : 380C WRKRSP : 381C INFTAP : LUSOL, LBSYMB 382C 383#include "pcmdef.h" 384#include "maxash.h" 385#include "maxorb.h" 386#include "mxcent.h" 387Ckr VTEX and QTEX should be brought into this subroutine through 388Ckr the call 389 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 390 CHARACTER*8 LABINT(9*MXCENT) 391 LOGICAL FILE_EXIST 392#include "orgcom.h" 393#include "priunit.h" 394#include "infinp.h" 395#include "pcm.h" 396#include "pcmlog.h" 397#include "inforb.h" 398#include "infvar.h" 399#include "infrsp.h" 400#include "wrkrsp.h" 401#include "inftap.h" 402#include "infpri.h" 403#include "infpar.h" 404#include "qm3.h" 405C 406C 407C 408 CALL QENTER('IEFLNO') 409C 410C Determine if full Hessian or only orbital Hessian 411C 412C 413 IF (IPRRSP .GE. 40) THEN 414 WRITE (LUPRI,'(//A)') ' --- TEST OUTPUT FROM IEFLNO ---' 415 END IF 416 IF (IPRRSP .GE. 140) THEN 417 WRITE (LUPRI,'(/A)') ' --- IEFLNO - svec(1,nosim) on entry' 418 CALL OUTPUT(SVEC,1,KZYVAR,1,NOSIM,KZYVAR,NOSIM,1,LUPRI) 419 END IF 420C 421 IF (.NOT. MMPCM) LADDMM = .TRUE. 422C 423C Core allocation 424C 425 KINTRP = 1 426 KINTAD = KINTRP + (3*MXCOOR + 1)/IRAT 427 KVTEX = KINTAD + (3*MXCOOR + 1)/IRAT 428 KQTEX = KVTEX + NTS*NOSIM 429 KUBO = KQTEX + NTS*NOSIM 430 KW10 = KUBO + NOSIM*N2ORBX 431C 432 KUCMO = KW10 433 KPCMX = KUCMO + NORBT*NBAST 434 435 IF (TRPLET) THEN 436 KPCMXT = KPCMX + NOSIM*N2ORBX 437 KPCMXA = KPCMXT + NOSIM*N2ORBX 438 KPCMXAT = KPCMXA + NOSIM*N2ASHX 439 KJPXAC = KPCMXAT + NOSIM*N2ASHX 440 ELSE 441 KPCMXA = KPCMX + NOSIM*N2ORBX 442 KJPXAC = KPCMXA + NOSIM*N2ASHX 443 ENDIF 444 IF (TRPLET) THEN 445 KJPXACT = KJPXAC + NOSIM 446 KW20 = KJPXACT + NOSIM 447 ELSE 448 KW20 = KJPXAC + NOSIM 449 ENDIF 450C 451 KJ1AO = KW20 452 KJ1SQ = KJ1AO + NNBASX*NSYM 453 KJ1 = KJ1SQ + N2ORBX 454 KJ1XSQ = KJ1 + NNORBX 455 KJ1XAC = KJ1XSQ + N2ORBX*NOSIM 456 IF (TRPLET) THEN 457 KJ1XACT = KJ1XAC + NOSIM * N2ASHX 458 KW30 = KJ1XACT + NOSIM * N2ASHX 459 ELSE 460 KW30 = KJ1XAC + NOSIM * N2ASHX 461 ENDIF 462C 463C 3.0 SOLSC 464 KOVLP = KW30 465 KW50 = KOVLP + NOSIM 466 LW50 = LFREE + 1 - KW50 467C 468 KNEED = MAX(KW30,KW50) 469 IF (KNEED .GT. LFREE) CALL ERRWRK('IEFLNO',KNEED,LFREE) 470C 471C Unpack symmetry blocked CMO 472C 473 CALL UPKCMO(CMO,WRK(KUCMO)) 474C 475C Calculate unpacked orbital trial vectors in UBO 476C 477 IF (NOSIM.GT.0) THEN 478 CALL RSPZYM(NOSIM,BOVECS,WRK(KUBO)) 479 CALL DSCAL(NOSIM*N2ORBX,-1.0D0,WRK(KUBO),1) 480 IF (IPRRSP .GE. 55) THEN 481 DO 210 IOSIM = 1,NOSIM 482 JUBO = KUBO + (IOSIM-1)*N2ORBX 483 WRITE (LUPRI,2110) IOSIM,NOSIM 484 CALL OUTPUT(WRK(JUBO),1,NORBT,1,NORBT,NORBT,NORBT,1, 485 & LUPRI) 486 210 CONTINUE 487 END IF 488 END IF 489 2110 FORMAT (/,'ABC Orbital trial vector unpacked to matrix form (no.', 490 * I3,' of',I3,')') 491C 492C Contributions from all tesserae are included. 493C 494C 495 CALL DZERO(WRK(KPCMX),NOSIM*N2ORBX) 496 IF (TRPLET) CALL DZERO(WRK(KPCMXT),NOSIM*N2ORBX) 497 REWIND LUPROP 498 XI = DIPORG(1) 499 YI = DIPORG(2) 500 ZI = DIPORG(3) 501#if defined (VAR_MPI) 502 IF (NODTOT .GE. 1 .AND. .NOT. MMPCM) THEN 503 CALL J1XP(NOSIM,WRK(KVTEX),WRK(KUCMO),WRK(KUBO),UDV,UDVTR, 504 & WRK(KPCMX),WRK(KPCMXT),.FALSE.,TRPLET,JWOPSY, 505 & WRK(KW50),LW50) 506 ELSE 507#endif 508 DO I = 1, NTSIRR 509C Read AO potential integrals on tesserae 510 L = 1 511 NCOMP = NSYM 512 DIPORG(1) = XTSCOR(I) 513 DIPORG(2) = YTSCOR(I) 514 DIPORG(3) = ZTSCOR(I) 515 EXP1VL = .FALSE. 516 TOFILE = .FALSE. 517 KPATOM = 0 518 TRIMAT = .TRUE. 519 CALL GET1IN(WRK(KJ1AO),'NPETES ',NCOMP,WRK(KW50),LW50, 520 & LABINT,WRK(KINTRP),WRK(KINTAD),L,TOFILE,KPATOM, 521 & TRIMAT,DUMMY,EXP1VL,DUMMY,IPRRSP) 522 JJ1AO = KJ1AO 523 CALL UTHU(WRK(JJ1AO),WRK(KJ1),WRK(KUCMO),WRK(KW30), 524 * NBAST,NORBT) 525C Transform V(MO) from triangular to square format 526 CALL DSPTSI(NORBT,WRK(KJ1),WRK(KJ1SQ)) 527 CALL DZERO(WRK(KJ1XSQ),N2ORBX*NOSIM) 528 DO IOSIM = 1, NOSIM 529 JUBO = KUBO + (IOSIM - 1) * N2ORBX 530 JJ1X = KJ1XSQ + (IOSIM - 1) * N2ORBX 531 JJ1XAC = KJ1XAC + (IOSIM - 1) * N2ASHX 532 CALL ONEXH1(WRK(JUBO),WRK(KJ1SQ),WRK(JJ1X)) 533C 534 IF (NASHT .GT. 0) THEN 535 CALL GETACQ(WRK(JJ1X),WRK(JJ1XAC)) 536 END IF 537 IF (IPRRSP .GE. 15) THEN 538 WRITE (LUPRI,'(/A,I5)') 'J1X_mo matrix: tess',I 539 CALL OUTPUT(WRK(JJ1X),1,NORBT,1,NORBT, 540 & NORBT,NORBT,1,LUPRI) 541 IF (NASHT .GT. 0) THEN 542 WRITE (LUPRI,'(/A)') ' J1X_ac matrix:' 543 CALL OUTPUT(WRK(JJ1XAC),1,NASHT,1,NASHT, 544 & NASHT,NASHT,1,LUPRI) 545 END IF 546 END IF 547C 548C Expectation value of transformed potential on tesserae: 549C <0|\tilde{V}|0> 550C 551 IF (KSYMOP .EQ. 1) THEN 552 IF (TRPLET) THEN 553 FACTOR = SLVTLM(UDVTR,WRK(JJ1XAC),WRK(JJ1X),TJ1XAC) 554 ELSE 555 FACTOR = SLVQLM(UDV,WRK(JJ1XAC),WRK(JJ1X),TJ1XAC) 556 END IF 557 WRK(KVTEX + (I-1) + (IOSIM-1)*NTS) = FACTOR 558 IF (IPRRSP .GE. 6) THEN 559 WRITE (LUPRI,'(A,F17.8)') 560 * ' --- J1X expectation value :', 561 & WRK(KVTEX + (I-1) + (IOSIM-1)*NTS) 562 WRITE (LUPRI,'(A,F17.8)') 563 * ' --- active part of J1X :',TJ1XAC 564 END IF 565 END IF 566 ENDDO 567 QFACTOR = QSN(I)+QSE(I) 568 IF (MMPCM) QFACTOR = QSN(I)+QSE(I)+QSMM(I) 569Cjje 570 IF(.not. SOPRPA) then 571 IF (TRPLET) THEN 572 CALL DAXPY(NOSIM*N2ORBX,QFACTOR,WRK(KJ1XSQ), 573 & 1,WRK(KPCMXT),1) 574 ELSE 575 CALL DAXPY(NOSIM*N2ORBX,QFACTOR,WRK(KJ1XSQ), 576 & 1,WRK(KPCMX),1) 577 ENDIF 578 END IF 579Cjje 580C KPCMX: \tilde{J} + \tilde{X} 581C 582C For non-totally symmetric perturbation operators 583C 584 IF (KSYMOP .NE. 1) THEN 585 ITS = (KSYMOP - 1)*NTSIRR + I 586C Transform AO pot. int. into MO basis V(AO) --> V(MO) 587 JJ1AO = KJ1AO + (KSYMOP - 1)*NNBASX 588 CALL UTHU(WRK(JJ1AO),WRK(KJ1),WRK(KUCMO),WRK(KW50), 589 * NBAST,NORBT) 590C Transform V(MO) from triangular to square format 591 CALL DSPTSI(NORBT,WRK(KJ1),WRK(KJ1SQ)) 592 CALL DZERO(WRK(KJ1XSQ),N2ORBX*NOSIM) 593 DO IOSIM = 1, NOSIM 594 JUBO = KUBO + (IOSIM - 1) * N2ORBX 595 JJ1X = KJ1XSQ + (IOSIM - 1) * N2ORBX 596 JJ1XAC = KJ1XAC + (IOSIM - 1) * N2ASHX 597 CALL ONEXH1(WRK(JUBO),WRK(KJ1SQ),WRK(JJ1X)) 598C 599 IF (NASHT .GT. 0) THEN 600 CALL GETACQ(WRK(JJ1X),WRK(JJ1XAC)) 601 END IF 602 IF (IPRRSP .GE. 15) THEN 603 WRITE (LUPRI,'(/A)') ' J1X_mo matrix :' 604 CALL OUTPUT(WRK(JJ1X),1,NORBT,1,NORBT, 605 & NORBT,NORBT,1,LUPRI) 606 IF (NASHT .GT. 0) THEN 607 WRITE (LUPRI,'(/A)') ' J1X_ac matrix :' 608 CALL OUTPUT(WRK(JJ1XAC),1,NASHT,1,NASHT, 609 & NASHT,NASHT,1,LUPRI) 610 END IF 611 END IF 612C 613C Expectation value of transformed potential on tesserae: 614C <0|\tilde{V}|0> 615C 616 IF (TRPLET) THEN 617 FACTOR = SLVTLM(UDVTR,WRK(JJ1XAC),WRK(JJ1X),TJ1XAC) 618 ELSE 619 FACTOR = SLVQLM(UDV,WRK(JJ1XAC),WRK(JJ1X),TJ1XAC) 620 END IF 621 WRK(KVTEX + (ITS-1) + (IOSIM-1)*NTS) = FACTOR 622 IF (IPRRSP .GE. 6) THEN 623 WRITE (LUPRI,'(A,F17.8)') 624 * ' --- J1X expectation value :', 625 & WRK(KVTEX + (ITS-1) + (IOSIM-1)*NTS) 626 WRITE (LUPRI,'(A,F17.8)') 627 * ' --- active part of J1X :',TJ1XAC 628 END IF 629 ENDDO 630C KPCMX: \tilde{J} + \tilde{X} 631 END IF 632 ENDDO 633 634#if defined (VAR_MPI) 635 END IF 636#endif 637 IF (IPRRSP .GE. 50) THEN 638 DO IOSIM = 1,NOSIM 639 JPCMX = KPCMX + (IOSIM-1)*N2ORBX 640 WRITE (LUPRI,'(/A,I3,A,I3)') 641 * ' --- IEFLNO - (JPCMX) half.',IOSIM,' of',NOSIM 642 CALL OUTPUT(WRK(JPCMX),1,NORBT,1,NORBT,NORBT,NORBT, 643 * 1,LUPRI) 644 END DO 645 END IF 646 DO IOSIM = 1, NOSIM 647 CALL V2Q(WRK(KW50),WRK(KVTEX + NTS*(IOSIM-1)), 648 & WRK(KQTEX + NTS*(IOSIM-1)), 649 & QTEXS,NEQRSP) 650 IF (IPRRSP .GE. 6) THEN 651 DO I = 1, NTS 652 WRITE (LUPRI,'(A,F17.8)') 653 * ' --- J2X expectation value :', 654 & WRK(KQTEX + NTS*(IOSIM-1) + I - 1) 655 END DO 656 END IF 657 ENDDO 658 CALL GPCLOSE(LUPCMD,'KEEP') 659C 660 IF (MMPCM) THEN 661 DO I = 1, NOSIM 662 DO ITS = 1, NTS 663 INDEX = NTS * (I-1) + ITS - 1 664 WRK(KQTEX + INDEX) = 665 & QSMM1(1 + INDEX) + WRK(KQTEX + INDEX) 666 ENDDO 667 ENDDO 668 ENDIF 669C 670C TRANSFORMED CHARGES MULTIPIED TO POTENTIALS. 671C There are only one-index-transformed charges for the totally 672C symmetric irrep. 673C 674 CALL DZERO(WRK(KJ1XSQ),NOSIM*NNBASX) 675 IF (KSYMOP .GT. 1) THEN 676 DO IOSIM = 1, NOSIM 677 KTSSIM = (IOSIM-1) * NTS 678 KTSSYMOP = (KSYMOP-1) * NTSIRR 679 CALL DCOPY(NTSIRR,WRK(KQTEX + KTSSIM + KTSSYMOP),1, 680 & WRK(KQTEX + KTSSIM),1) 681 CALL DZERO(WRK(KQTEX + KTSSIM + KTSSYMOP), NTSIRR) 682 END DO 683 END IF 684 CALL J1INT(WRK(KQTEX),.FALSE.,WRK(KJ1XSQ),NOSIM,.FALSE., 685 & 'NPETES ',KSYMOP,WRK(KW50),LW50) 686 DO IOSIM = 1, NOSIM 687 JAOP = KJ1XSQ + (IOSIM - 1) * NNBASX 688 JPCMX = KPCMX + (IOSIM - 1) * N2ORBX 689 CALL UTHU(WRK(JAOP),WRK(KJ1),WRK(KUCMO),WRK(KW50), 690 * NBAST,NORBT) 691 CALL DSPTSI(NORBT,WRK(KJ1),WRK(KJ1SQ)) 692 CALL DAXPY(N2ORBX,-D1,WRK(KJ1SQ),1, 693 & WRK(JPCMX),1) 694 END DO 695C 696C Restore dipole origin 697C 698 DIPORG(1) = XI 699 DIPORG(2) = YI 700 DIPORG(3) = ZI 701C 702 IF (IPRRSP .GE. 50) THEN 703 DO IOSIM = 1,NOSIM 704 JPCMX = KPCMX + (IOSIM-1)*N2ORBX 705 WRITE (LUPRI,'(/A,I3,A,I3)') 706 * ' --- IEFLNO - (JPCMX) matrix no.',IOSIM,' of',NOSIM 707 CALL OUTPUT(WRK(JPCMX),1,NORBT,1,NORBT,NORBT,NORBT, 708 * 1,LUPRI) 709 END DO 710 END IF 711C 712Clara add the case of KZCONF.GT.0 713 IF ((.NOT.TDHF).AND.(KZCONF.GT.0)) THEN 714C 715C ... CSF part of sigma vectors 716C 717 DO 700 IOSIM = 1,NOSIM 718 JPCMX = KPCMX + (IOSIM-1)*N2ORBX 719 JPCMXA = KPCMXA + (IOSIM-1)*N2ASHX 720 IF (TRPLET) THEN 721 JPCMXT = KPCMXT + (IOSIM-1)*N2ORBX 722 JPCMXAT = KPCMXAT + (IOSIM-1)*N2ASHX 723 ENDIF 724 IF (IREFSY .EQ. KSYMST) THEN 725 WRK(KOVLP-1+IOSIM) = DDOT(KZCONF,CREF,1,SVEC(1,IOSIM),1) 726 ELSE 727 WRK(KOVLP-1+IOSIM) = D0 728 END IF 729 CALL GETACQ(WRK(JPCMX),WRK(JPCMXA)) 730C <0|\tilde{Z}|0> (Z=J+X+V*Q just to remind) 731 TJPX = SLVQLM(UDV,WRK(JPCMXA),WRK(JPCMX),TJPXAC) 732 WRK(KJPXAC-1+IOSIM) = TJPXAC 733 IF (TRPLET) THEN 734 CALL GETACQ(WRK(JPCMXT),WRK(JPCMXAT)) 735 TJPXT = SLVQLM(UDVTR,WRK(JPCMXAT),WRK(JPCMXT),TJPXACT) 736 WRK(KJPXACT-1+IOSIM) = TJPXACT 737 ELSE 738 IF (IPRRSP .GE. 40) WRITE (LUPRI,'(/A,I5,3F15.10)') 739 * ' IOSIM, C_OVLP, TJPXAC, TJPX :', 740 * IOSIM,WRK(KOVLP-1+IOSIM),TJPXAC,TJPX 741 ENDIF 742 700 CONTINUE 743C \sigma_{co} = <i|\tilde{Z}|0> 744C KW30??? 745 CALL SLVSC(0,NOSIM,N2ASHX,DUMMY,CREF,SVEC,WRK(KPCMXA), 746 * WRK(KPCMXAT), 747 * WRK(KJPXAC),DUMMY,INDXCI,WRK(KW50),LW50) 748C CALL SLVSC(NCSIM,NOSIM,BCVECS,CREF,SVECS, 749C * RXAC,RYAC,TRXAC,TRYAC,INDXCI,WRK,LWRK) 750 IF (IPRRSP.GT.101) THEN 751 WRITE(LUPRI,*)' LINEAR TRANSFORMED ORBITAL VECTOR' 752 WRITE(LUPRI,*)' **** AFTER SLVSC in IEFLNO **** ' 753 CALL OUTPUT(SVEC,1,KZYVAR,1,NOSIM,KZYVAR,NOSIM,1,LUPRI) 754 END IF 755 END IF 756C \sigma_{oo} = <0|[q_j,\tilde{Z}]|0> 757 IF (.NOT. LADDMM) THEN 758 CALL DZERO(WRK(KPCMX),NOSIM*N2ORBX) 759 IF (TRPLET) CALL DZERO(WRK(KPCMXT),NOSIM*N2ORBX) 760 ENDIF 761 762 IF(TRPLET) THEN 763 CALL SLVSOR(.TRUE.,.FALSE.,NOSIM,UDVTR,SVEC(1,1),WRK(KPCMX)) 764 CALL SLVSOR(.TRUE.,.TRUE., NOSIM,UDV, SVEC(1,1),WRK(KPCMXT)) 765 ELSE 766 CALL SLVSOR(.TRUE.,.TRUE., NOSIM,UDV, SVEC(1,1),WRK(KPCMX)) 767 ENDIF 768C allow for SOPRPA 769C IF (KZCONF.GT.0 .AND. IREFSY.EQ.KSYMST) THEN 770 IF ((KZCONF.GT.0) .AND. (IREFSY.EQ.KSYMST) 771 * .AND. (.NOT. SOPRPA)) THEN 772 IF (NCREF .NE. KZCONF) CALL QUIT('IEFLNO: NCREF .ne. KZCONF') 773C 774C ... test orthogonality 775C 776 DO 900 IOSIM = 1,NOSIM 777 TEST = DDOT(KZCONF,CREF,1,SVEC(1,IOSIM),1) 778 * - WRK(KOVLP-1+IOSIM) 779 IF (ABS(TEST) .GT. 1.D-8) THEN 780 NWARN = NWARN + 1 781 WRITE (LUPRI,'(/A,I5,/A,1P,D12.4)') 782 * ' --- IEFLNO WARNING, for IOSIM =',IOSIM, 783 * ' <CREF | SVEC_solvent(iosim) > =',TEST 784 END IF 785 900 CONTINUE 786 END IF 787C 788CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 789C 790C ... test print 791C 792 IF (IPRRSP .GE. 140) THEN 793 WRITE (LUPRI,'(/A)') ' --- IEFLNO - svec(ci,1) on exit' 794 DO 930 I = 1,KZCONF 795 IF (SVEC(I,1) .NE. D0) WRITE (LUPRI,'(A,I10,F15.10)') 796 * ' conf #',I,SVEC(I,1) 797 930 CONTINUE 798 END IF 799 IF (IPRRSP .GE. 140) THEN 800 WRITE (LUPRI,'(/A)') ' --- IEFLNO - svec(orb,1) on exit' 801 WRITE (LUPRI,'(/A)') ' Z - PART OR VECTOR ' 802 CALL OUTPUT(SVEC(KZCONF+1,1),1,KZWOPT,1,1,KZWOPT,1,1,LUPRI) 803 WRITE (LUPRI,'(/A)') ' Y - PART OR VECTOR ' 804 CALL OUTPUT(SVEC(KZVAR+KZCONF+1,1),1,KZWOPT,1,1,KZWOPT, 805 * 1,1,LUPRI) 806 END IF 807C 808 CALL QEXIT('IEFLNO') 809 RETURN 810C ... end of IEFLNO. 811 END 812 813C /* Deck pcm_comp_pot */ 814 SUBROUTINE PCM_COMP_POT(MOCOEF, DENMAT, DENTRP, BOVECS, 815 & XJSQ, XJSQT, NOSIM, KSYMP, WORK, LWORK) 816#include "implicit.h" 817#include "priunit.h" 818#include "dummy.h" 819#include "thrzer.h" 820#include "iratdef.h" 821#include "mxcent.h" 822#include "pcmdef.h" 823#include "pcm.h" 824#include "pcmlog.h" 825#include "qm3.h" 826#include "orgcom.h" 827#include "dftcom.h" 828#include "inforb.h" 829#include "infrsp.h" 830#include "infvar.h" 831 832 DIMENSION MOCOEF(*), DENMAT(*), DENTRP(*), BOVECS(*), 833 & XJSQ(*), XJSQT(*), WORK(*) 834 IF (VSNFLG .EQ. -1) THEN 835 CALL DZERO(VSN, NTS) 836 CALL PCM_COMP_NUC_POT(VSN) 837 VSNFLG = 1 838 END IF 839 IF (VSEFLG .EQ. -1) THEN 840 CALL DZERO(VSE, NTS) 841 CALL PCM_COMP_EL_POT(VSE, DENMAT, NOSIM, NBAS, NSYM, NNBASX, 842 & KSYMP, WORK, LWORK) 843 VSEFLG = 1 844 END IF 845 IF (VSE1FLG .EQ. -1) THEN 846 CALL DZERO(VSE1, NTS*NOSIM) 847 CALL PCM_COMP_EL_POT1(VSE1, NOSIM, DENMAT, DENTRP, JWOPSY, 848 & MOCOEF, BOVECS, XJSQ, XJSQT, 849 & WORK, LWORK) 850 VSE1FLG = 1 851 END IF 852 IF (VSMMFLG .EQ. -1) THEN 853 CALL DZERO(VSMM, NTS) 854 CALL PCM_COMP_MM_POT(VSMM) 855 VSMMFLG = 1 856 END IF 857 IF (VSMM1FLG .EQ. -1) THEN 858 CALL DZERO(VSMM1, NTS*NOSIM) 859 CALL PCM_COMP_MM_POT1(VSMM1, NTS, XTSCOR, YTSCOR, ZTSCOR, 860 & WORK, LWORK) 861 VSMM1FLG = 1 862 END IF 863 RETURN 864 END 865 866C /* Deck pcm_comp_nuc_pot */ 867 SUBROUTINE PCM_COMP_NUC_POT(VSN) 868#include "implicit.h" 869 CALL COMP_NUC_POT_CAV(VSN, .TRUE., .FALSE.,.FALSE.) 870 RETURN 871 END 872 873C /* Deck pcm_comp_el_pot */ 874 SUBROUTINE PCM_COMP_EL_POT(VSE, DENMAT, NOSIM, NBAS, NSYM, NNBASX, 875 & KSYMP, WORK, LWORK) 876#include "implicit.h" 877 DIMENSION WORK(*), DENMAT(*), VSE(*) 878 KDEN = 1 879 KFREE = KDEN + NNBASX 880 LFREE = LWORK - KFREE + 1 881 IF (LFREE .LT. 0) THEN 882 CALL QUIT('Not enough memory in pcm_comp_el_pot') 883 ENDIF 884 885 CALL PKSYM1(WORK(KDEN),DENMAT,NBAS,NSYM,-1) 886 887 CALL J1INT(VSE, .TRUE., WORK(KDEN), NOSIM, .FALSE., 'NPETES ', 888 & KSYMP, WORK(KFREE), LFREE) 889 RETURN 890 END 891 892C /* Deck pcm_comp_el_pot1 */ 893 SUBROUTINE PCM_COMP_EL_POT1(VSE1, NOSIM, UDV, UDVTR, JWOPSY, 894 & CMO, BOVECS, XJSQ, XJSQT, 895 & WORK, LWORK) 896#include "implicit.h" 897#include "inforb.h" 898#include "priunit.h" 899#include "infrsp.h" 900#include "maxorb.h" 901#include "infpar.h" 902#include "qm3.h" 903 904 LOGICAL TOFILE 905 DIMENSION BOVECS(*), UDV(*), UDVTR(*), CMO(*), WORK(*) 906 907C iprrsp = 1 908 KUBO = 1 909 KUCMO = KUBO + NOSIM * N2ORBX 910 KFREE = KUCMO + NORBT * NBAST 911 LFREE = LWORK - KFREE + 1 912 IF (LFREE .LT. 0) THEN 913 CALL QUIT('PCM_COMP_EL_POT1: insufficient memory') 914 END IF 915 CALL UPKCMO(CMO,WORK(KUCMO)) 916 917 IF (NOSIM.GT.0) THEN 918 CALL RSPZYM(NOSIM,BOVECS,WORK(KUBO)) 919 CALL DSCAL(NOSIM*N2ORBX,-1.0D0,WORK(KUBO),1) 920 IF (IPRRSP .GE. 55) THEN 921 DO 210 IOSIM = 1,NOSIM 922 JUBO = KUBO + (IOSIM-1)*N2ORBX 923 WRITE (LUPRI,2110) IOSIM,NOSIM 924 CALL OUTPUT(WORK(JUBO),1,NORBT,1,NORBT,NORBT,NORBT,1, 925 & LUPRI) 926 210 CONTINUE 927 END IF 928 END IF 929 2110 FORMAT (/,'Orbital trial vector unpacked to matrix form (no.', 930 * I3,' of',I3,')') 931 932#if defined (VAR_MPI) 933 IF (NODTOT .GE. 1 .AND. .NOT. MMPCM) THEN 934 CALL J1XP(NOSIM,VSE1,WORK(KUCMO),WORK(KUBO),UDV,UDVTR, 935 & XJ1SQ,XJ1SQT, 936 & .FALSE.,TRPLET,JWOPSY,WORK(KFREE),LFREE) 937 ELSE 938#endif 939 CALL J1X(NOSIM,VSE1,WORK(KUCMO),WORK(KUBO),UDV,UDVTR, 940 & TOFILE,JWOPSY,WORK(KFREE),LFREE) 941#if defined (VAR_MPI) 942 ENDIF 943#endif 944 RETURN 945 END 946 947C /* Deck pcm_comp_mm_pot */ 948 SUBROUTINE PCM_COMP_MM_POT(VCAV, NTS, XTS, YTS, ZTS, XMM, YMM, 949 & ZMM, SOURCE, NSOURCES, TYPE) 950#include "implicit.h" 951 CHARACTER*3 TYPE 952 IF (TYPE .EQ. 'CHG') THEN 953 CALL PCM_COMP_POT_CHG(VCAV, NTS, XTS, YTS, ZTS, XMM, YMM, 954 & ZMM, SOURCE, NSOURCES) 955 ELSE IF (TYPE .EQ. 'DIP') THEN 956 CALL PCM_COMP_POT_DIP(VCAV, NTS, XTS, YTS, ZTS, XMM, YMM, 957 & ZMM, SOURCE, NSOURCES) 958Clf quadrupole not yet imlemented 959Clf ELSE IF (TYPE .EQ. 'QDR') THEN 960Clf PCM_COMP_POT_QDR(VCAV, NTS, COORD, SOURCE, NSOURCES) 961Clf quadrupole components order xx, xy, xz, yy, yz, zz 962Clf quadrupole potential: sum_{a,b} 1/2*(3*Ra*Rb-R^2*delta_ab)/R^5 Q_ab 963Clf 964 ELSE 965 CALL QUIT('Wrong type in PCM_COMP_MM_POT') 966 END IF 967 RETURN 968 END 969 970C /* Deck pcm_comp_mm_pot1 */ 971 SUBROUTINE PCM_COMP_MM_POT1(VCAV, NTS, XTS, YTS, ZTS, 972 & WORK, LWORK) 973#include "implicit.h" 974 CHARACTER*3 TYPE 975 DIMENSION WORK(*), VCAV(*), XTS(*), YTS(*), ZTS(*) 976 977 CALL PCM_READ_MM_INDIP(NOSIMMM, KXMM, KYMM, KZMM, KSOURCES, 978 & NSOURCES, WORK, KFREE, LWORK) 979 DO I = 1, NOSIMMM 980 JXMM = KXMM + NSOURCES * (I-1) 981 JYMM = KYMM + NSOURCES * (I-1) 982 JZMM = KZMM + NSOURCES * (I-1) 983 JSOURCES = KSOURCES + 3 * (I-1) * NSOURCES 984 JPOT = 1 + NTS * (I-1) 985 CALL PCM_COMP_POT_DIP(VCAV(JPOT), NTS, XTS, YTS, ZTS, 986 & WORK(JXMM), WORK(JYMM), WORK(JZMM), 987 & WORK(JSOURCES), NSOURCES) 988 END DO 989 RETURN 990 END 991 992C /* Deck pcm_comp_mm_pot */ 993 SUBROUTINE PCM_COMP_POT_CHG(VCAV, NTS, XTS, YTS, ZTS, XMM, YMM, 994 & ZMM, SOURCE, NSITES) 995#include "implicit.h" 996 DIMENSION VCAV(*), XTS(*), YTS(*), ZTS(*),XMM(*),YMM(*),ZMM(*), 997 & SOURCE(*) 998 DO IATOM = 1, NSITES 999 XL = XMM(IATOM) 1000 YL = YMM(IATOM) 1001 ZL = ZMM(IATOM) 1002 CHGMM = SOURCE(IATOM) 1003 DO ITS = 1, NTS 1004 RIL=DSQRT((XTS(ITS)-XL)**2+(YTS(ITS)-YL)**2 1005 & +(ZTS(ITS)-ZL)**2) 1006 VCAV(ITS) = VCAV(ITS) + CHGMM/RIL 1007 ENDDO 1008 ENDDO 1009 RETURN 1010 END 1011 1012C /* Deck pcm_comp_pot_dip */ 1013 SUBROUTINE PCM_COMP_POT_DIP(VCAV, NTS, XTS, YTS, ZTS, XMM, YMM, 1014 & ZMM, SOURCE, NSITES) 1015#include "implicit.h" 1016 DIMENSION VCAV(*), XTS(*), YTS(*), ZTS(*), 1017 & XMM(*), YMM(*),ZMM(*), 1018 & SOURCE(*) 1019 DO IATOM = 1, NSITES 1020 XDIP = SOURCE(3*IATOM-2) 1021 YDIP = SOURCE(3*IATOM-1) 1022 ZDIP = SOURCE(3*IATOM) 1023 DO ITS = 1, NTS 1024 XL = XTS(ITS) - XMM(IATOM) 1025 YL = YTS(ITS) - YMM(IATOM) 1026 ZL = ZTS(ITS) - ZMM(IATOM) 1027 SCALAR = XL * XDIP + YL * YDIP + ZL * ZDIP 1028 R2 = XL ** 2 + YL ** 2 + ZL ** 2 1029 R3 = (DSQRT(R2)) ** 3 1030 VCAV(ITS) = VCAV(ITS) + SCALAR/R3 1031 ENDDO 1032 ENDDO 1033 RETURN 1034 END 1035 1036C 1037C /* Deck pcm_read_mm_indip */ 1038C 1039 SUBROUTINE PCM_READ_MM_INDIP(NOSIMMM, KXMM, KYMM, KZMM, KSOURCES, 1040 & NSOURCES, WORK, KFREE, LWORK) 1041#include "implicit.h" 1042#include "priunit.h" 1043 LOGICAL FILE_EXIST 1044 DIMENSION WORK(*) 1045 INQUIRE(FILE='MYFPCM',EXIST=FILE_EXIST) 1046 IF (FILE_EXIST) THEN 1047 LQM3PCM = -1 1048 CALL GPOPEN(LQM3PCM,'MYFPCM','OLD',' ','FORMATTED ',IDUMMY, 1049 & .FALSE.) 1050 REWIND(LQM3PCM) 1051 READ(LQM3PCM,*) NOSIMMM,NSOURCES 1052 ELSE 1053 CALL QUIT('There is no input from QM/MM response 1054 & to PCM response') 1055 ENDIF 1056 KXMM = 1 1057 KYMM = KXMM + NSOURCES * NOSIMMM 1058 KZMM = KYMM + NSOURCES * NOSIMMM 1059 KSOURCES = KZMM + NSOURCES * NOSIMMM 1060 KFREE = KSOURCES + 3 * NSOURCES * NOSIMMM 1061 LFREE = LWORK - KFREE + 1 1062 L = 0 1063 DO I = 1, NOSIMMM 1064 DO J = 1, NSOURCES 1065 L = L + 1 1066 READ(LQM3PCM,'(6(E25.15,2x))') 1067 & WORK(KXMM + L - 1), 1068 & WORK(KYMM + L - 1), 1069 & WORK(KZMM + L - 1), 1070 & WORK(KSOURCES + 3 * L - 3 ), 1071 & WORK(KSOURCES + 3 * L - 2 ), 1072 & WORK(KSOURCES + 3 * L - 1 ) 1073 1074 ENDDO 1075 ENDDO 1076 CALL GPCLOSE(LQM3PCM,'KEEP') 1077 LWORK = LFREE 1078 RETURN 1079 END 1080C 1081C /* Deck pcm_write_pcm2mm */ 1082C 1083 SUBROUTINE PCM_WRITE_PCM2MM(NOSIMMM,NTS,XTSCOR,YTSCOR,ZTSCOR, 1084 & QSE1,QSMM1) 1085C Arnfinn Apr. -09 1086C Make the file QFQM3 with induced charges for use 1087C by QM3 response calculations. 1088#include "implicit.h" 1089#include "priunit.h" 1090C#include "mxcent.h" 1091C#include "pcmdef.h" 1092C#include "pcm.h" 1093 1094 LOGICAL FILE_EXIST 1095 DIMENSION XTSCOR(*),YTSCOR(*),ZTSCOR(*),QSE1(*),QSMM1(*) 1096 1097 INQUIRE(FILE='QFQM3',EXIST=FILE_EXIST) 1098 IF (FILE_EXIST) THEN 1099 LPCMQM3 = -1 1100 CALL GPOPEN(LPCMQM3,'QFQM3','OLD',' ','FORMATTED', 1101 & IDUMMY,.FALSE.) 1102 CALL GPCLOSE(LPCMQM3,'DELETE') 1103 ENDIF 1104 LPCMQM3 = -1 1105 CALL GPOPEN(LPCMQM3,'QFQM3','NEW',' ','FORMATTED', 1106 & IDUMMY,.FALSE.) 1107 WRITE (LPCMQM3,*) NOSIMMM, NTS 1108 DO I = 1, NOSIMMM 1109 DO ITS = 1, NTS 1110 INDEX = NTS * (I-1) + ITS 1111 CHG1 = QSE1(INDEX)+QSMM1(INDEX) 1112 WRITE (LPCMQM3,'(4(E25.15,2x))') XTSCOR(ITS),YTSCOR(ITS), 1113 & ZTSCOR(ITS), CHG1 1114 ENDDO 1115 ENDDO 1116 CALL GPCLOSE(LPCMQM3,'KEEP') 1117 RETURN 1118 END 1119