1 2! 3! Dalton, a molecular electronic structure program 4! Copyright (C) by the authors of Dalton. 5! 6! This program is free software; you can redistribute it and/or 7! modify it under the terms of the GNU Lesser General Public 8! License version 2.1 as published by the Free Software Foundation. 9! 10! This program is distributed in the hope that it will be useful, 11! but WITHOUT ANY WARRANTY; without even the implied warranty of 12! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13! Lesser General Public License for more details. 14! 15! If a copy of the GNU LGPL v2.1 was not distributed with this 16! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 17! 18! 19C CBN+JK 20C************************************************************************* 21C /* Deck qm3fck */ 22 SUBROUTINE QM3FCK(DCAO,DVAO,FSOL,EQM3T,ENSEL,EPOL,EELEL, 23 * WRK,LWRK,IPRINT) 24 25#include "implicit.h" 26#include "dummy.h" 27#include "inftap.h" 28#include "priunit.h" 29#include "qm3.h" 30#include "mxcent.h" 31#include "thrzer.h" 32#include "iratdef.h" 33#include "codata.h" 34#include "maxash.h" 35#include "maxorb.h" 36#include "infinp.h" 37#include "inforb.h" 38#include "infpri.h" 39#include "scbrhf.h" 40#include "maxaqn.h" 41#include "symmet.h" 42#include "orgcom.h" 43#include "infpar.h" 44 45 DIMENSION DCAO(*), DVAO(*) 46 DIMENSION FSOL(*), WRK(LWRK) 47 DIMENSION DFTNS(MXQM3) 48 49 50 CALL QENTER('QM3FCK') 51 KDTAO = 1 52 KTAO = KDTAO + NNBASX 53 KEND = KTAO + NNBASX 54 LWRK1 = LWRK - KEND 55 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3FCK',-KEND,LWRK) 56 57 58C Get total density 59 IF (NASHT .EQ. 0) THEN 60 CALL PKSYM1(WRK(KDTAO),DCAO,NBAS,NSYM,-1) 61 ELSE 62 DO I = 1,NNBAST 63 IF (HSROHF) THEN 64 WRK(KTAO-1+I) = DCAO(I) 65 ELSE 66 WRK(KTAO-1+I) = DCAO(I) + DVAO(I) 67 END IF 68 END DO 69 CALL PKSYM1(WRK(KDTAO),WRK(KTAO),NBAS,NSYM,-1) 70 END IF 71 72C From COMMON 73 ENSQM3 = 0.0D0 74 EPOQM3 = 0.0D0 75 76C Calculate RRa, Obar and ENSA and write these to file 77 CALL QM3_RRA(WRK(KDTAO),WRK(KEND),LWRK1,IPRINT) 78 79C Calculate Ns and keep in memory (DFTNS). 80C (Run QM3_NSP instead of QM3_NS if #nodes > 1, Arnfinn nov. -08) 81#if defined(VAR_MPI) 82 IF (NODTOT.GE.1 .AND. .NOT. OLDTG .AND. .NOT. MMPCM) THEN 83 CALL QM3_NSP(WRK(KDTAO),DFTNS,WRK(KEND),LWRK1,IPRINT) 84 ELSE 85#endif 86 CALL QM3_NS(WRK(KDTAO),DFTNS,WRK(KEND),LWRK1,IPRINT) 87#if defined(VAR_MPI) 88 END IF 89#endif 90 91C Modify the fock operator. Modification returned in FSOL. 92 CALL QM3_FMO(FSOL,WRK(KEND),LWRK1,IPRINT) 93 94C Calculate the QM3 contribution to the energy (returned in EQM3T) 95 CALL QM3_ENERGY(DFTNS,ENSEL,EPOL,EELEL,EQM3T,WRK(KEND),LWRK1, 96 & IPRINT) 97 98 CALL QEXIT('QM3FCK') 99 RETURN 100 END 101C *********************************************************************** 102C /* Deck qm3_rra */ 103 SUBROUTINE QM3_RRA(DCAO,WRK,LWRK,IPRINT) 104 105#include "implicit.h" 106#include "priunit.h" 107#include "dummy.h" 108#include "qm3.h" 109#include "mxcent.h" 110#include "iratdef.h" 111#include "maxash.h" 112#include "maxorb.h" 113#include "inforb.h" 114#include "inftap.h" 115#include "infpri.h" 116#include "scbrhf.h" 117#include "maxaqn.h" 118#include "symmet.h" 119#include "orgcom.h" 120#include "infinp.h" 121#include "nuclei.h" 122 123 124 DIMENSION WRK(LWRK) 125 DIMENSION EELEC(3,MXQM3) 126 DIMENSION FFJ(3) 127 LOGICAL LOINDM 128 CHARACTER*8 LABINT(9*MXCENT) 129 LOGICAL TOFILE, TRIMAT, EXP1VL 130 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 131 132C From here. 133C This is NOT good. In fact very bad. Needs to be fixed. 134C Neede for f.ex. spin-spin couplings since nucdep is used here for allocation. 135C Can only be done AFTER coordinates for charges and polarizabilities have been defined 136C since nudep defines no. of coordinates !! Only active when requsted by the user. 137 138 IF (REDCNT) NUCDEP=NUCIND 139 140C To here 141 142 IF (LOSPC) RETURN 143 144 CALL QENTER('QM3_RRA') 145 146 KRRAOx = 1 147 KRRAOy = KRRAOx + NNBASX 148 KRRAOz = KRRAOy + NNBASX 149 KRAx = KRRAOz + NNBASX 150 KRAy = KRAx + NCOMS 151 KRAz = KRAy + NCOMS 152 KWRK1 = KRAz + NCOMS 153 LWRK1 = LWRK - KWRK1 154 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3_RRA',-KWRK1,LWRK) 155 156 CALL DZERO(WRK(KRRAOx),3*NNBASX) 157 CALL DZERO(WRK(KRAx),3*NCOMS) 158 159 DO 879 I = 1, MXQM3 160 DO 880 J = 1, 3 161 EELEC(J,I) = 0.0D0 162 880 CONTINUE 163 879 CONTINUE 164 165 IF (.NOT. INTDIR) THEN 166 IF (IQM3PR .GT. 15) THEN 167 WRITE(LUPRI,*) 'QM3RAINT: Read in integrals' 168 ENDIF 169 170 CALL GPOPEN(LUQM3E,'ELFDMM','OLD',' ', 171 & 'UNFORMATTED',IDUMMY,.FALSE.) 172 REWIND(LUQM3E) 173 END IF 174 175 LM = 0 176 177 IF (INTDIR) THEN 178 L = NUSITE + NUALIS(0) 179 OBKPX = DIPORG(1) 180 OBKPY = DIPORG(2) 181 OBKPZ = DIPORG(3) 182 END IF 183 184 DO 520 I = 1, ISYTP 185 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 186 DO 521 J = NSYSBG(I), NSYSED(I) 187 DO 522 K = 1, NUALIS(I) 188 LM = LM + 1 189 190 IF (INTDIR) THEN 191 KMAT = KWRK1 192 KLAST = KMAT + 3*NNBASX 193 LWRK2 = LWRK - KLAST + 1 194 IATNOW = NUCIND + L + LM 195 196 KPATOM = 0 197 NOSIMI = 3 198 TOFILE = .FALSE. 199 TRIMAT = .TRUE. 200 EXP1VL = .FALSE. 201 DIPORG(1) = CORD(1,IATNOW) 202 DIPORG(2) = CORD(2,IATNOW) 203 DIPORG(3) = CORD(3,IATNOW) 204 205 RUNQM3 = .TRUE. 206 CALL GET1IN(WRK(KMAT),'NEFIELD',NOSIMI,WRK(KLAST), 207 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 208 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 209 RUNQM3 = .FALSE. 210 211 IF (IPRINT.GT.15) THEN 212 WRITE (LUPRI,'(/A)') ' Rra_ao_x matrix:' 213 CALL OUTPAK(WRK(KMAT),NBAST,1,LUPRI) 214 215 WRITE (LUPRI,'(/A)') ' Rra_ao_y matrix:' 216 CALL OUTPAK(WRK(KMAT+NNBASX),NBAST,1,LUPRI) 217 218 WRITE (LUPRI,'(/A)') ' Rra_ao_z matrix:' 219 CALL OUTPAK(WRK(KMAT+2*NNBASX),NBAST,1,LUPRI) 220 END IF 221 222 IF (QMDAMP) THEN 223 DIST = (DIPORG(1)-QMCOM(1))**2 + 224 & (DIPORG(2)-QMCOM(2))**2 + 225 & (DIPORG(3)-QMCOM(3))**2 226 DIST = SQRT(DIST) 227 DFACT = (1-exp(-ADAMP*DIST))**3 228 CALL DSCAL(3*NNBASX,DFACT,WRK(KMAT),1) 229 ENDIF 230 231 WRK(KRAx+LM-1) = DDOT(NNBASX,WRK(KMAT),1,DCAO,1) 232 WRK(KRAy+LM-1) = DDOT(NNBASX,WRK(KMAT+NNBASX),1,DCAO,1) 233 WRK(KRAz+LM-1) = DDOT(NNBASX,WRK(KMAT+2*NNBASX), 234 * 1,DCAO,1) 235 ELSE 236 237 CALL READT(LUQM3E,NNBASX,WRK(KRRAOx)) 238 WRK(KRAx+LM-1) = DDOT(NNBASX,WRK(KRRAOx),1,DCAO,1) 239 240 IF (IPRINT.GT.15) THEN 241 CALL AROUND('Rra_ao_x matrix:') 242 CALL OUTPAK(WRK(KRRAOx),NBAST,1,LUPRI) 243 END IF 244 245 CALL READT(LUQM3E,NNBASX,WRK(KRRAOy)) 246 WRK(KRAy+LM-1) = DDOT(NNBASX,WRK(KRRAOy),1,DCAO,1) 247 248 IF (IPRINT.GT.15) THEN 249 CALL AROUND('Rra_ao_y matrix:') 250 CALL OUTPAK(WRK(KRRAOy),NBAST,1,LUPRI) 251 END IF 252 253 CALL READT(LUQM3E,NNBASX,WRK(KRRAOz)) 254 WRK(KRAz+LM-1) = DDOT(NNBASX,WRK(KRRAOz),1,DCAO,1) 255 256 IF (IPRINT.GT.15) THEN 257 CALL AROUND('Rra_ao_z matrix:') 258 CALL OUTPAK(WRK(KRRAOz),NBAST,1,LUPRI) 259 END IF 260 261 END IF 262 522 CONTINUE 263 521 CONTINUE 264 END IF 265 520 CONTINUE 266 267 IF (INTDIR) THEN 268 DIPORG(1) = OBKPX 269 DIPORG(2) = OBKPY 270 DIPORG(3) = OBKPZ 271 END IF 272 273 274 IF (IQM3PR .GE. 5) THEN 275 WRITE(LUPRI,'(/A/A/A)') 276 * ' +==========================================================+', 277 * ' | COM| <Rra>_x | <Rra>_y | <Rra>_z |', 278 * ' +==========================================================+' 279 280 NUM = 0 281 282 DO 215 I = 1, ISYTP 283 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 284 DO 216 J = NSYSBG(I), NSYSED(I) 285 DO 217 K=1,NUALIS(I) 286 287 NUM = NUM + 1 288 289 WRITE(LUPRI,'(A,I3,A,3(F16.10,A)/A)') 290 * ' | ', J,'|', WRK(KRAx + NUM - 1),' |', 291 * WRK(KRAy + NUM - 1),' |', WRK(KRAz + NUM - 1),' |', 292 * ' +----------------------------------------------------------+' 293 217 CONTINUE 294 216 CONTINUE 295 END IF 296 215 CONTINUE 297 WRITE(LUPRI,'(//,A)') 298 END IF 299 300 IF (LM .NE. NCOMS) THEN 301 CALL QUIT('Error in no. of center of masses in QM3RAINT') 302 END IF 303 304 DO 534 LM = 1,NCOMS 305 EELEC(1,LM) = WRK(KRAx+LM-1) 306 EELEC(2,LM) = WRK(KRAy+LM-1) 307 EELEC(3,LM) = WRK(KRAz+LM-1) 308 534 CONTINUE 309 310C If RELMOM is true we want to include the external field(s) in 311C the determination of the induced dipole moments 312 313 FFJ(1) = 0.0D0 314 FFJ(2) = 0.0D0 315 FFJ(3) = 0.0D0 316 317 IF (RELMOM) THEN 318 DO 330 IF =1, NFIELD 319 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(1) = FFJ(1) + EFIELD(IF) 320 IF (LFIELD(IF) .EQ. 'YDIPLEN') FFJ(2) = FFJ(2) + EFIELD(IF) 321 IF (LFIELD(IF) .EQ. 'ZDIPLEN') FFJ(3) = FFJ(3) + EFIELD(IF) 322 330 CONTINUE 323 END IF 324 325 IF (FIXMOM) THEN 326 WRITE(LUPRI,'(/A)')'FIXMOM: NO ITER. DETERM. OF MYIND' 327 ELSE IF (LOSPC) THEN 328 WRITE(LUPRI,'(/A)')'All MM models are SPC, INDMOM not called' 329 ELSE 330 LOINDM = .FALSE. 331 CALL INDMOM(EELEC,LOINDM,FFJ) 332 END IF 333 334 IF (FIXMOM .AND. MMPCM) CALL MMPCM_FIXMOM() 335 336 CALL QM3_OBAR(FFJ) 337 338 CALL CC_PUT31('CC_RA',NCOMS,WRK(KRAx),WRK(KRAy),WRK(KRAz)) 339 340 IF (.NOT. INTDIR) CALL GPCLOSE(LUQM3E,'KEEP') 341 342 CALL QEXIT('QM3_RRA') 343 RETURN 344 END 345C********************************************************************* 346C /* Deck qm3_ns */ 347 SUBROUTINE QM3_NS(DCAO,DFTNS,WRK,LWRK,IPRINT) 348 349#include "implicit.h" 350#include "priunit.h" 351#include "dummy.h" 352#include "qm3.h" 353#include "mxcent.h" 354#include "iratdef.h" 355#include "maxash.h" 356#include "maxorb.h" 357#include "inforb.h" 358#include "inftap.h" 359#include "infpri.h" 360#include "scbrhf.h" 361#include "maxaqn.h" 362#include "symmet.h" 363#include "orgcom.h" 364#include "nuclei.h" 365 366 367 DIMENSION WRK(LWRK) 368 DIMENSION DFTNS(*), DCAO(NNBASX) 369 CHARACTER*8 LABINT(9*MXCENT) 370 LOGICAL TOFILE, TRIMAT, EXP1VL 371 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 372 373 374 CALL QENTER('QM3_NS') 375 376C From here. 377C This is NOT good. In fact very bad. Needs to be fixed. 378C Neede for f.ex. spin-spin couplings since nucdep is used here for allocation. 379C Can only be done AFTER coordinates for charges and polarizabilities have been defined 380C since nudep defines no. of coordinates !! Only active when requsted by the user. 381 382 IF (REDCNT) NUCDEP=NUCIND 383 384C To here 385 386 387 KNSAO = 1 388 KWRK1 = KNSAO + NNBASX 389 LWRK1 = LWRK - KWRK1 390 391 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3_NS',-KWRK1,LWRK) 392 393 CALL DZERO(WRK(KNSAO),NNBASX) 394 395 ENSEL = 0.0D0 396 397 IF (.NOT. INTDIR) THEN 398 CALL GPOPEN(LUQM3P,'POTMM','UNKNOWN',' ', 399 & 'UNFORMATTED',IDUMMY,.FALSE.) 400 REWIND (LUQM3P) 401 ENDIF 402 403 FAC1 = -1.0D0 404 405 L = 0 406 407 IF (INTDIR) THEN 408 OBKPX = DIPORG(1) 409 OBKPY = DIPORG(2) 410 OBKPZ = DIPORG(3) 411 ENDIF 412 413 DO 525 I = 1, ISYTP 414 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 415 DO 526 J = NSYSBG(I), NSYSED(I) 416 DO 527 K = 1,NSISY(I) 417 418 L = L +1 419 420 IF (INTDIR) THEN 421 422 KMAT = KWRK1 423 KLAST = KMAT + NNBASX 424 LWRK2 = LWRK - KLAST + 1 425 IATNOW = NUCIND + L 426 427 KPATOM = 0 428 NOSIMI = 1 429 TOFILE = .FALSE. 430 TRIMAT = .TRUE. 431 EXP1VL = .FALSE. 432 DIPORG(1) = CORD(1,IATNOW) 433 DIPORG(2) = CORD(2,IATNOW) 434 DIPORG(3) = CORD(3,IATNOW) 435 436 RUNQM3=.TRUE. 437 CALL GET1IN(WRK(KMAT),'NPETES ',NOSIMI,WRK(KLAST), 438 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 439 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 440 RUNQM3=.FALSE. 441 IF (OLDTG) THEN 442 FACQ = -1.0D0*CHAOLD(IATNOW) 443 ELSE 444 FACQ = -1.0D0*CHARGE(IATNOW) 445 ENDIF 446 447 CALL DSCAL(NNBASX,FACQ,WRK(KMAT),1) 448 449 IF (IPRINT.GT.15) THEN 450 WRITE (LUPRI,'(/A)') 'NUCPOT matrix' 451 CALL OUTPAK(WRK(KMAT),NBAST,1,LUPRI) 452 ENDIF 453 454 DFTNS(L) = DDOT(NNBASX,DCAO,1,WRK(KMAT),1) 455C MM/PCM interface. Coordinates and charges are fixed 456C here so in fact on needed on the first call. Leave it 457C for now. 458 IF (MMPCM) THEN 459 XMMQ(L) = CORD(1,IATNOW) 460 YMMQ(L) = CORD(2,IATNOW) 461 ZMMQ(L) = CORD(3,IATNOW) 462 MMQ(L) = CHARGE(IATNOW) 463 IF (L.GT.MXQ) THEN 464 CALL QUIT('Too many charges in MM.'// 465 & ' Increase MXQ') 466 ENDIF 467 ENDIF 468 ELSE 469 CALL READT(LUQM3P,NNBASX,WRK(KNSAO)) 470 IF (IPRINT.GT.15) THEN 471 WRITE (LUPRI,'(/A,I3,A)') 472 * ' N(',L,')_ao matrix: ' 473 CALL OUTPAK(WRK(KNSAO),NBAST,1,LUPRI) 474 ENDIF 475 476 DFTNS(L) = DDOT(NNBASX,DCAO,1,WRK(KNSAO),1) 477 ENSEL = ENSEL - DFTNS(L) 478 ENDIF 479 527 CONTINUE 480 526 CONTINUE 481 END IF 482 525 CONTINUE 483 484 IF (INTDIR) THEN 485 DIPORG(1) = OBKPX 486 DIPORG(2) = OBKPY 487 DIPORG(3) = OBKPZ 488 ENDIF 489 490 IF (.NOT. INTDIR) CALL GPCLOSE(LUQM3P,'KEEP') 491 492C------------------- 493C Print section: 494C------------------- 495C 496 IF (IQM3PR .GT. 3) THEN 497 WRITE(LUPRI,'(/A/A/A)') 498 * ' +======================+', 499 * ' |Site| <N_s> |', 500 * ' +======================+' 501 502 LS = 0 503 504 DO 215 I = 1, ISYTP 505 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 506 DO 216 J = NSYSBG(I), NSYSED(I) 507 DO 217 K = 1, NSISY(I) 508 509 LS = LS + 1 510 511 WRITE(LUPRI,'(A,I3,A,F16.10,A/A)') 512 * ' | ', LS,'|', DFTNS(LS),' |', 513 * ' +----------------------+' 514 217 CONTINUE 515 216 CONTINUE 516 END IF 517 215 CONTINUE 518 WRITE(LUPRI,'()') 519 END IF 520 CALL QEXIT('QM3_NS') 521 RETURN 522 END 523 524 525#if defined(VAR_MPI) 526C********************************************************************* 527C /* Deck qm3_nsp */ 528 SUBROUTINE QM3_NSP(DCAO,DFTNS,WRK,LWRK,IPRINT) 529C 530C The parallel version of the routine QM3_NS, Arnfinn nov. -08 531C 532#include "implicit.h" 533#include "priunit.h" 534#include "dummy.h" 535#include "qm3.h" 536#include "mxcent.h" 537#include "iratdef.h" 538#include "maxash.h" 539#include "maxorb.h" 540#include "inforb.h" 541#include "inftap.h" 542#include "infpri.h" 543#include "scbrhf.h" 544#include "maxaqn.h" 545#include "symmet.h" 546#include "orgcom.h" 547#include "nuclei.h" 548#include "mtags.h" 549#include "infpar.h" 550 551 DIMENSION WRK(LWRK) 552 DIMENSION DFTNS(MXQM3), DCAO(NNBASX) 553 CHARACTER*8 LABINT(9*MXCENT) 554 LOGICAL TOFILE, TRIMAT, EXP1VL 555 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 556 557 558C defined parallel calculation types 559#include "iprtyp.h" 560 IPRTYP = QM3_NSP_WORK 561 562 CALL QENTER('QM3_NSP') 563 564C From here. TODO TODO TODO 565C This is NOT good. In fact very bad. Needs to be fixed. 566C Needed for f.ex. spin-spin couplings since nucdep is used here for allocation. 567C Can only be done AFTER coordinates for charges and polarizabilities have been defined 568C since nudep defines no. of coordinates !! Only active when requsted by the user. 569 570 IF (REDCNT) NUCDEP=NUCIND 571 572C To here 573 574CHJAaJ-b KDFTNS not used /June 09 575C KDFTNS = 1 576C KWRK1 = KDFTNS + MXQM3 577CHJAaJ-e 578 KWRK1 = 1 579 LWRK1 = LWRK - KWRK1 + 1 580 581 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3_NSP',-KWRK1,LWRK) 582 583 ENSEL = 0.0D0 584 585 FAC1 = -1.0D0 586 587 L = 0 588 589 OBKPX = DIPORG(1) 590 OBKPY = DIPORG(2) 591 OBKPZ = DIPORG(3) 592 593 CALL DZERO(DFTNS,MXQM3) 594 595 DO 525 I = 1, ISYTP 596 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 597 CALL QM3_NSPM(IPRTYP,DCAO,DFTNS,WRK(KWRK1), 598 & LWRK1,IPRINT,EXP1VL,.FALSE.,ENSEL,FAC1,L,I) 599 600 END IF 601 525 CONTINUE 602 603 DIPORG(1) = OBKPX 604 DIPORG(2) = OBKPY 605 DIPORG(3) = OBKPZ 606 607C------------------- 608C Print section: 609C------------------- 610C 611 IF (IQM3PR .GT. 3) THEN 612 WRITE(LUPRI,'(/A/A/A)') 613 * ' +======================+', 614 * ' |Site| <N_s> |', 615 * ' +======================+' 616 617 LS = 0 618 619 DO 215 I = 1, ISYTP 620 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 621 DO 216 J = NSYSBG(I), NSYSED(I) 622 DO 217 K = 1, NSISY(I) 623 624 LS = LS + 1 625 626 WRITE(LUPRI,'(A,I3,A,F16.10,A/A)') 627 * ' | ', LS,'|', DFTNS(LS),' |', 628 * ' +----------------------+' 629 217 CONTINUE 630 216 CONTINUE 631 END IF 632 215 CONTINUE 633 WRITE(LUPRI,'()') 634 END IF 635 CALL QEXIT('QM3_NSP') 636 RETURN 637 END 638C********************************************************************* 639 SUBROUTINE QM3_NSPM(IPRTYP,DCAO,DFTNS,WRK,LWRK,IPRINT, 640 & EXP1VL,TOFILE,ENSEL,FAC1,LNUM,INUM) 641#include "implicit.h" 642#include "priunit.h" 643#include "dummy.h" 644#include "qm3.h" 645#include "mxcent.h" 646#include "iratdef.h" 647#include "maxash.h" 648#include "maxorb.h" 649#include "inforb.h" 650#include "inftap.h" 651#include "infpri.h" 652#include "scbrhf.h" 653#include "maxaqn.h" 654#include "symmet.h" 655#include "orgcom.h" 656#include "nuclei.h" 657#include "infpar.h" 658#include "mtags.h" 659#if defined(VAR_MPI) 660#include "mpif.h" 661#endif 662 663#include "cbiher.h" 664 665C The master routine 666 667 LOGICAL TOFILE, EXP1VL 668 DIMENSION WRK(LWRK) 669 DIMENSION DCAO(NNBASX) 670 DIMENSION DFTNS(MXQM3) 671 672 IF (TOFILE) CALL QUIT('Parallel calculations do not allow '// 673 & 'for storing NS-integrals on disk') 674 675 KDCAO=1 676 KDFTNS1=KDCAO+NNBASX 677 KDFTNS2=KDFTNS1+MXQM3 678 KLAST=KDFTNS2+MXQM3 679 LWRK1=LWRK-KLAST+1 680 CALL MPIXBCAST(IPRTYP,1,'INTEGER',MASTER) 681 CALL MPIXBCAST(IQM3PR,1,'INTEGER',MASTER) 682C 683 CALL MPIXBCAST(NBAST,1,'INTEGER',MASTER) 684 CALL MPIXBCAST(INTDIR,1,'LOGICAL',MASTER) 685 CALL MPIXBCAST(ISYTP,1,'INTEGER',MASTER) 686 CALL MPIXBCAST(NCTOT,1,'INTEGER',MASTER) 687 CALL MPIXBCAST(DCAO,NNBASX,'DOUBLE',MASTER) 688 CALL MPIXBCAST(CORD,3*MXCENT,'DOUBLE',MASTER) 689 CALL MPIXBCAST(CHARGE,MXCENT,'DOUBLE',MASTER) 690 CALL MPIXBCAST(TOFILE,1,'LOGICAL',MASTER) 691 CALL MPIXBCAST(ENSEL,1,'DOUBLE',MASTER) 692 CALL MPIXBCAST(FAC1,1,'DOUBLE',MASTER) 693 CALL MPIXBCAST(NUCIND,1,'INTEGER',MASTER) 694 CALL MPIXBCAST(NUCDEP,1,'INTEGER',MASTER) 695 CALL MPIXBCAST(LUPROP,1,'INTEGER',MASTER) 696 CALL MPIXBCAST(NPATOM,1,'INTEGER',MASTER) 697 CALL MPIXBCAST(IPATOM,MXCENT,'INTEGER',MASTER) 698 699C Loop over all MM nuclei 700 701C ENSEL = 0.0D0 702C FAC1 = -1.0D0 703 704 DO J = NSYSBG(INUM), NSYSED(INUM) 705 DO K = 1, NSISY(INUM) 706 LNUM = LNUM + 1 707 IWHO = -1 708 CALL MPIXRECV(NWHO, 1, 'INTEGER', IWHO, MPTAG1) 709 CALL MPIXSEND(LNUM,1,'INTEGER',NWHO,MPTAG2) 710 END DO 711 END DO 712 713C Send end message to all slaves 714 715 LEND = -1 716 DO ISLAVE = 1, NODTOT 717 IWHO = -1 718 CALL MPIXRECV(NWHO,1,'INTEGER',IWHO,MPTAG1) 719 CALL MPIXSEND(LEND ,1,'INTEGER',NWHO,MPTAG2) 720 END DO 721 722C Collect data from all slaves 723 CALL DZERO(WRK(KDFTNS1),MXQM3) 724 CALL DZERO(WRK(KDFTNS2),MXQM3) 725 CALL MPI_REDUCE(WRK(KDFTNS1),WRK(KDFTNS2),MXQM3, 726 & MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD, 727 & IERR) 728 DO I = 1, MXQM3 729 DFTNS(I) = DFTNS(I) + WRK(KDFTNS2 + I - 1) 730 END DO 731 732 RETURN 733 END 734 735C********************************************************************* 736 737 SUBROUTINE QM3_NSPS(WRK,LWRK,IPRTMP) 738#include "implicit.h" 739#include "priunit.h" 740#include "dummy.h" 741#include "gnrinf.h" 742#include "qm3.h" 743#include "mxcent.h" 744#include "iratdef.h" 745#include "maxash.h" 746#include "maxorb.h" 747#include "inforb.h" 748#include "inftap.h" 749#include "infpri.h" 750#include "scbrhf.h" 751#include "maxaqn.h" 752#include "symmet.h" 753#include "orgcom.h" 754#include "nuclei.h" 755#include "infpar.h" 756#include "mtags.h" 757#if defined(VAR_MPI) 758#include "mpif.h" 759#endif 760 761#include "cbiher.h" 762 763C The slave routine 764 765 CHARACTER*8 LABINT(9*MXCENT) 766 LOGICAL TOFILE, EXP1VL, TRIMAT 767 DIMENSION WRK(LWRK) 768 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 769 770 IQM3PR = IPRTMP 771 QM3 = .TRUE. 772 773 CALL MPIXBCAST(NBAST,1,'INTEGER',MASTER) 774 CALL MPIXBCAST(INTDIR,1,'LOGICAL',MASTER) 775 CALL MPIXBCAST(ISYTP,1,'INTEGER',MASTER) 776 CALL MPIXBCAST(NCTOT,1,'INTEGER',MASTER) 777 KDCAO = 1 778 KDFTNS = KDCAO + NNBASX 779 KLAST1 = KDFTNS + MXQM3 780 LWRK1 = LWRK - KLAST1 + 1 781 CALL MPIXBCAST(WRK(KDCAO),NNBASX,'DOUBLE',MASTER) 782 CALL MPIXBCAST(CORD,3*MXCENT,'DOUBLE',MASTER) 783 CALL MPIXBCAST(CHARGE,MXCENT,'DOUBLE',MASTER) 784 CALL MPIXBCAST(TOFILE,1,'LOGICAL',MASTER) 785 CALL MPIXBCAST(ENSEL,1,'DOUBLE',MASTER) 786 CALL MPIXBCAST(FAC1,1,'DOUBLE',MASTER) 787 CALL MPIXBCAST(NUCIND,1,'INTEGER',MASTER) 788 CALL MPIXBCAST(NUCDEP,1,'INTEGER',MASTER) 789 CALL MPIXBCAST(LUPROP,1,'INTEGER',MASTER) 790 CALL MPIXBCAST(NPATOM,1,'INTEGER',MASTER) 791 CALL MPIXBCAST(IPATOM,MXCENT,'INTEGER',MASTER) 792 793 KMAT = KLAST1 794 KLAST2 = KMAT + NNBASX 795 LWRK2 = LWRK1 - KLAST2 + KLAST1 796 797 CALL DZERO(WRK(KDFTNS),MXQM3) 798 799C Run loop over MM nuclear charges 800 801 20 CONTINUE 802 803 CALL MPIXSEND(MYNUM,1,'INTEGER',MASTER,MPTAG1) 804 CALL MPIXRECV(L,1,'INTEGER',MASTER,MPTAG2) 805 CALL DZERO(WRK(KMAT),NNBASX) 806 IF (L .GT. 0) THEN 807 IATNOW = NUCIND + L 808 KPATOM = 0 809 NOSIMI = 1 810 TOFILE = .FALSE. 811 TRIMAT = .TRUE. 812 EXP1VL = .FALSE. 813 DIPORG(1) = CORD(1,IATNOW) 814 DIPORG(2) = CORD(2,IATNOW) 815 DIPORG(3) = CORD(3,IATNOW) 816 817 RUNQM3=.TRUE. 818 CALL GET1IN(WRK(KMAT),'NPETES ',NOSIMI,WRK(KLAST2), 819 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 820 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 821 822 RUNQM3=.FALSE. 823 IF (OLDTG) THEN 824 FACQ = -1.0D0*CHAOLD(IATNOW) 825 ELSE 826 FACQ = -1.0D0*CHARGE(IATNOW) 827 ENDIF 828 CALL DSCAL(NNBASX,FACQ,WRK(KMAT),1) 829 WRK(KDFTNS + L - 1) = DDOT(NNBASX,WRK(KDCAO),1,WRK(KMAT),1) 830 ENSEL = ENSEL - WRK(KDFTNS + L - 1) 831 GO TO 20 832 END IF 833 834C No more Ns to calculate 835 CALL MPI_REDUCE(WRK(KDFTNS),MPI_IN_PLACE,MXQM3, 836 & MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD, 837 & IERR) 838 839 RETURN 840 END 841#endif 842C********************************************************************* 843C /* Deck qm3_fmo */ 844 SUBROUTINE QM3_FMO(FSOL,WRK,LWRK,IPRINT) 845 846#include "implicit.h" 847#include "priunit.h" 848#include "dummy.h" 849#include "qm3.h" 850#include "mxcent.h" 851#include "iratdef.h" 852#include "maxash.h" 853#include "maxorb.h" 854#include "inforb.h" 855#include "inftap.h" 856#include "infpri.h" 857#include "scbrhf.h" 858#include "maxaqn.h" 859#include "symmet.h" 860#include "orgcom.h" 861#include "infinp.h" 862#include "nuclei.h" 863 864 865 DIMENSION WRK(LWRK) 866 DIMENSION FSOL(*) 867 CHARACTER*8 LABINT(9*MXCENT) 868 LOGICAL TOFILE, TRIMAT, EXP1VL 869 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 870 871 872 IF (LOSPC .AND. (.NOT. OLDTG)) RETURN 873 874 CALL QENTER('QM3_FMO') 875 876 KRRAOx = 1 877 KRRAOy = KRRAOx + NNBASX 878 KRRAOz = KRRAOy + NNBASX 879 KRAx = KRRAOz + NNBASX 880 KRAy = KRAx + NCOMS 881 KRAz = KRAy + NCOMS 882 KENSAx = KRAz + NCOMS 883 KENSAy = KENSAx + NCOMS 884 KENSAz = KENSAy + NCOMS 885 KTAO = KENSAz + NCOMS 886 KWRK1 = KTAO + NNBASX 887 LWRK1 = LWRK - KWRK1 888 889 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3_FMO',-KWRK1,LWRK) 890 891 CALL DZERO(WRK(KRRAOx),3*NNBASX) 892 CALL DZERO(WRK(KRAx),6*NCOMS) 893 CALL DZERO(WRK(KTAO),NNBASX) 894 895 IF (.NOT. LOSPC) THEN 896 897 CALL CC_GET31('CC_RA',NCOMS, 898 * WRK(KRAx),WRK(KRAy),WRK(KRAz)) 899 900 CALL CC_GET31('ENSAFILE',NCOMS, 901 * WRK(KENSAx),WRK(KENSAy),WRK(KENSAz)) 902 903 IF (IQM3PR .GT. 5) THEN 904 WRITE(LUPRI,'(/A/A/A)') 905 * ' +==========================================================+', 906 * ' | COM| <Rra>_x | <Rra>_y | <Rra>_z |', 907 * ' +==========================================================+' 908 909 NUM = 0 910 911 DO 215 I = 1, ISYTP 912 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 913 DO 216 J = NSYSBG(I), NSYSED(I) 914 DO 217 K=1,NUALIS(I) 915 916 NUM = NUM + 1 917 918 WRITE(LUPRI,'(A,I3,A,F16.10,A,F16.10,A,F16.10,A/A)') 919 * ' | ', J,'|', WRK(KRAx + NUM - 1),' |', 920 * WRK(KRAy + NUM - 1),' |', WRK(KRAz + NUM - 1),' |', 921 * '+----------------------------------------------------------+' 922 217 CONTINUE 923 216 CONTINUE 924 END IF 925 215 CONTINUE 926 927 WRITE(LUPRI,'(/A/A/A)') 928 * ' +==========================================================+', 929 * ' | COM| ENSA_x | ENSA_y | ENSA_z |', 930 * ' +==========================================================+' 931 932 NUM = 0 933 934 DO 415 I = 1, ISYTP 935 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 936 937 NUM = NUM + 1 938 939 DO 416 J = NSYSBG(I), NSYSED(I) 940 DO 417 K=1,NUALIS(I) 941 WRITE(LUPRI,'(A,I3,A,F16.10,A, 942 & F16.10,A,F16.10,A/A)') 943 * ' | ', J,'|', WRK(KENSAx + NUM - 1),' |', 944 * WRK(KENSAy + NUM - 1),' |', WRK(KENSAz + NUM - 1),' |', 945 * ' +----------------------------------------------------------+' 946 417 CONTINUE 947 416 CONTINUE 948 END IF 949 415 CONTINUE 950 END IF 951 952 IF (.NOT. INTDIR) THEN 953 CALL GPOPEN(LUQM3E,'ELFDMM','OLD',' ', 954 & 'UNFORMATTED',IDUMMY,.FALSE.) 955 REWIND(LUQM3E) 956 ENDIF 957 958 LM = 0 959 960 IF (INTDIR) THEN 961 L = NUSITE + NUALIS(0) 962 OBKPX = DIPORG(1) 963 OBKPY = DIPORG(2) 964 OBKPZ = DIPORG(3) 965 ENDIF 966 967 DO 520 I = 1, ISYTP 968 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 969 DO 521 J = NSYSBG(I), NSYSED(I) 970 DO 522 K = 1, NUALIS(I) 971 LM = LM + 1 972 973 IF (INTDIR) THEN 974 KMAT = KWRK1 975 KLAST = KMAT + 3*NNBASX 976 LWRK2 = LWRK - KLAST + 1 977 IATNOW = NUCIND + L + LM 978 979 KPATOM = 0 980 NOSIMI = 3 981 TOFILE = .FALSE. 982 TRIMAT = .TRUE. 983 EXP1VL = .FALSE. 984 DIPORG(1) = CORD(1,IATNOW) 985 DIPORG(2) = CORD(2,IATNOW) 986 DIPORG(3) = CORD(3,IATNOW) 987 988 RUNQM3=.TRUE. 989 CALL GET1IN(WRK(KMAT),'NEFIELD',NOSIMI,WRK(KLAST), 990 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 991 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 992 RUNQM3=.FALSE. 993 994 IF (QMDAMP) THEN 995 DIST = (DIPORG(1)-QMCOM(1))**2 + 996 & (DIPORG(2)-QMCOM(2))**2 + 997 & (DIPORG(3)-QMCOM(3))**2 998 DIST = SQRT(DIST) 999 DFACT = (1-exp(-ADAMP*DIST))**3 1000 CALL DSCAL(3*NNBASX,DFACT,WRK(KMAT),1) 1001 ENDIF 1002 1003 ELSE 1004 1005 CALL READT(LUQM3E,NNBASX,WRK(KRRAOx)) 1006 CALL READT(LUQM3E,NNBASX,WRK(KRRAOy)) 1007 CALL READT(LUQM3E,NNBASX,WRK(KRRAOz)) 1008 1009 ENDIF 1010 1011 IF (MDLWRD(I) .EQ. 'SPC_E01') THEN 1012 1013 FACx = - ALPIMM(I,K) * (WRK(KRAx + LM - 1) 1014 & + 0.5D0 * WRK(KENSAx + LM - 1)) 1015 FACy = - ALPIMM(I,K) * (WRK(KRAy + LM - 1) 1016 & + 0.5D0 * WRK(KENSAy + LM - 1)) 1017 FACz = - ALPIMM(I,K) * (WRK(KRAz + LM - 1) 1018 & + 0.5D0 * WRK(KENSAz + LM - 1)) 1019 1020 IF (INTDIR) THEN 1021 CALL DAXPY(NNBASX,FACx,WRK(KMAT),1,WRK(KTAO),1) 1022 CALL DAXPY(NNBASX,FACy,WRK(KMAT+NNBASX),1, 1023 * WRK(KTAO),1) 1024 CALL DAXPY(NNBASX,FACz,WRK(KMAT+2*NNBASX),1, 1025 * WRK(KTAO),1) 1026 ELSE 1027 CALL DAXPY(NNBASX,FACx,WRK(KRRAOx),1,WRK(KTAO),1) 1028 CALL DAXPY(NNBASX,FACy,WRK(KRRAOy),1,WRK(KTAO),1) 1029 CALL DAXPY(NNBASX,FACz,WRK(KRRAOz),1,WRK(KTAO),1) 1030 ENDIF 1031 1032 ENDIF 1033 522 CONTINUE 1034 521 CONTINUE 1035 END IF 1036 520 CONTINUE 1037 1038 IF (.NOT. INTDIR) CALL GPCLOSE(LUQM3E,'KEEP') 1039 1040 IF (IQM3PR .GE. 15) THEN 1041 WRITE (LUPRI,'(/A)') ' QM/MM contribution:' 1042 CALL OUTPAK(WRK(KTAO),NBAST,1,LUPRI) 1043 END IF 1044 1045 END IF ! LOSPC 1046 1047 IF (OLDTG) THEN 1048 1049C We use RRAOx in this case for the pot.energy integrals 1050 CALL DZERO(WRK(KRRAOx),NNBASX) 1051 1052 IF (.NOT. INTDIR) THEN 1053 CALL GPOPEN(LUQM3P,'POTMM','UNKNOWN',' ', 1054 & 'UNFORMATTED',IDUMMY,.FALSE.) 1055 REWIND (LUQM3P) 1056 ENDIF 1057 1058 FAC1 = -1.0D0 1059 1060 L = 0 1061 1062 IF (INTDIR) THEN 1063 OBKPX = DIPORG(1) 1064 OBKPY = DIPORG(2) 1065 OBKPZ = DIPORG(3) 1066 ENDIF 1067 1068 DO 525 I = 1, ISYTP 1069 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 1070 DO 526 J = NSYSBG(I), NSYSED(I) 1071 DO 527 K = 1,NSISY(I) 1072 1073 L = L +1 1074 1075 IF (INTDIR) THEN 1076 1077 KMAT = KWRK1 1078 KLAST = KMAT + NNBASX 1079 LWRK2 = LWRK - KLAST + 1 1080 1081 IATNOW = NUCIND + L 1082 1083 CALL DZERO(WRK(KMAT),NNBASX) 1084 1085 KPATOM = 0 1086 NOSIMI = 1 1087 TOFILE = .FALSE. 1088 TRIMAT = .TRUE. 1089 EXP1VL = .FALSE. 1090 DIPORG(1) = CORD(1,IATNOW) 1091 DIPORG(2) = CORD(2,IATNOW) 1092 DIPORG(3) = CORD(3,IATNOW) 1093 1094 RUNQM3=.TRUE. 1095 CALL GET1IN(WRK(KMAT),'NPETES ',NOSIMI,WRK(KLAST), 1096 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 1097 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 1098 RUNQM3=.FALSE. 1099 1100 FACQ = -1.0D0*CHAOLD(IATNOW) 1101 1102 CALL DSCAL(NNBASX,FACQ,WRK(KMAT),1) 1103 1104 IF (IPRINT.GT.15) THEN 1105 WRITE (LUPRI,'(/A)') 'NUCPOT matrix in QM3_FMO' 1106 CALL OUTPAK(WRK(KMAT),NBAST,1,LUPRI) 1107 ENDIF 1108 1109 CALL DAXPY(NNBASX,FAC1,WRK(KMAT), 1110 * 1,WRK(KTAO),1) 1111 1112 ELSE 1113 1114 CALL READT(LUQM3P,NNBASX,WRK(KRRAOx)) 1115 1116 IF (IPRINT.GT.15) THEN 1117 WRITE (LUPRI,'(/A,I3,A)') 1118 * ' N(',L,')_ao matrix: ' 1119 CALL OUTPAK(WRK(KRRAOx),NBAST,1,LUPRI) 1120 ENDIF 1121 1122 CALL DAXPY(NNBASX,FAC1,WRK(KRRAOx), 1123 * 1,WRK(KTAO),1) 1124 1125 ENDIF 1126 1127 527 CONTINUE 1128 526 CONTINUE 1129 END IF 1130 525 CONTINUE 1131 1132 IF (INTDIR) THEN 1133 DIPORG(1) = OBKPX 1134 DIPORG(2) = OBKPY 1135 DIPORG(3) = OBKPZ 1136 ENDIF 1137 1138 IF (.NOT. INTDIR) CALL GPCLOSE(LUQM3P,'KEEP') 1139 1140 END IF 1141 1142 CALL PKSYM1(WRK(KTAO),FSOL,NBAS,NSYM,+1) 1143 1144 CALL QEXIT('QM3_FMO') 1145 RETURN 1146 END 1147C********************************************************************* 1148C /* Deck qm3_energy */ 1149 SUBROUTINE QM3_ENERGY(DFTNS,EEL,EPOL,EELEL,EQM3T,WRK,LWRK,IPRINT) 1150 1151#include "implicit.h" 1152#include "priunit.h" 1153#include "dummy.h" 1154#include "qm3.h" 1155#include "mxcent.h" 1156#include "iratdef.h" 1157#include "maxash.h" 1158#include "maxorb.h" 1159#include "inforb.h" 1160#include "inftap.h" 1161#include "infpri.h" 1162#include "scbrhf.h" 1163#include "maxaqn.h" 1164#include "symmet.h" 1165#include "orgcom.h" 1166#include "infinp.h" 1167 1168 DIMENSION WRK(LWRK) 1169 DIMENSION FFJ(3), DFTNS(*) 1170 1171 CALL QENTER('QM3_ENERGY') 1172 IF ( .NOT. (LOSPC) ) THEN 1173 1174 KOMMSx = 1 1175 KOMMSy = KOMMSx + NCOMS 1176 KOMMSz = KOMMSy + NCOMS 1177 KRAx = KOMMSz + NCOMS 1178 KRAy = KRAx + NCOMS 1179 KRAz = KRAy + NCOMS 1180 KWRK1 = KRAz + NCOMS 1181 LWRK1 = LWRK - KWRK1 1182 1183 IF (LWRK1.LT.0) CALL QUIT( 'Too little work in QM3_ENERGY') 1184 1185 CALL DZERO(WRK(KOMMSx),6*NCOMS) 1186 1187 END IF 1188 1189C First we see if any static fields are to be added 1190 1191 FFJ(1) = 0.0D0 1192 FFJ(2) = 0.0D0 1193 FFJ(3) = 0.0D0 1194 1195 IF (RELMOM) THEN 1196 DO 330 IF =1, NFIELD 1197 IF (LFIELD(IF) .EQ. 'XDIPLEN') FFJ(1) = FFJ(1) + EFIELD(IF) 1198 IF (LFIELD(IF) .EQ. 'YDIPLEN') FFJ(2) = FFJ(2) + EFIELD(IF) 1199 IF (LFIELD(IF) .EQ. 'ZDIPLEN') FFJ(3) = FFJ(3) + EFIELD(IF) 1200 330 CONTINUE 1201 END IF 1202 1203C------------------------------------- 1204C Add up the energy contributions: 1205C 1) Epol: 1206C------------------------------------- 1207 1208 EQM3T = 0.0D0 1209 1210 IF (.NOT. LOSPC) THEN 1211 CALL CC_GET31('OBARFILE',NCOMS, 1212 * WRK(KOMMSx),WRK(KOMMSy),WRK(KOMMSz)) 1213 1214 CALL CC_GET31('CC_RA',NCOMS, 1215 * WRK(KRAx),WRK(KRAy),WRK(KRAz)) 1216 1217 KS = 0 1218 DO 110 I = 1, ISYTP 1219 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1220 DO 120 J = NSYSBG(I), NSYSED(I) 1221 DO 121 K = 1, NUALIS(I) 1222 KS = KS+1 1223 RAT = 0.0D0 1224 RAT = 0.5D0 * ALPIMM(I,K) * WRK(KRAx + KS - 1) * 1225 & (WRK(KRAx + KS - 1) + WRK(KOMMSx + KS - 1)) 1226 & + 0.5D0 * ALPIMM(I,K) * WRK(KRAy + KS - 1) * 1227 & (WRK(KRAy + KS - 1) + WRK(KOMMSy + KS - 1)) 1228 & + 0.5D0 * ALPIMM(I,K) * WRK(KRAz + KS - 1) * 1229 & (WRK(KRAz + KS - 1) + WRK(KOMMSz + KS - 1)) 1230 EQM3T = EQM3T - RAT 1231 121 CONTINUE 1232 120 CONTINUE 1233 END IF 1234 110 CONTINUE 1235 1236 TEMP1 = EQM3T 1237 1238 CALL QM3_OTILDE(OTILDE,FFJ) 1239 1240 ELSE 1241 1242 OTILDE = 0.0D0 1243 TEMP1 = 0.0D0 1244 1245 END IF 1246 1247 EQM3T = EQM3T + OTILDE 1248 EPOL = EQM3T 1249C-------- 1250C 2) Eel: 1251C-------- 1252C 1253 ETEMP = 0.0D0 1254 L = 0 1255 1256 DO 130 I = 1, ISYTP 1257 IF (MDLWRD(I)(1:3) .EQ. 'SPC') THEN 1258 DO 140 J = NSYSBG(I), NSYSED(I) 1259 DO 150 K = 1, NSISY(I) 1260 L = L + 1 1261 ETEMP = ETEMP - DFTNS(L) 1262 150 CONTINUE 1263 140 CONTINUE 1264 END IF 1265 130 CONTINUE 1266 1267 EELEL = ETEMP 1268 1269 EQM3T = EQM3T + ENUQM3 + EELEL 1270 1271 EEL = ENUQM3 + EELEL 1272 1273C------------- 1274C 3) E(QM/MM): 1275C------------- 1276 1277 EQM3T = EQM3T + ECLVDW 1278 1279C----------------------------------------------- 1280C Testing the energy contributions one by one: 1281C----------------------------------------------- 1282 1283 IF (IQM3PR .GT. 1) THEN 1284 WRITE(LUPRI,'(A)') 1285 * ' |xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx ' 1286 * //'The different energy contributions outlined' 1287 * //' xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx|' 1288 WRITE(LUPRI,'(A)') 1289 * ' | Evdw | E-nuclear |' 1290 * //' -Sum_s <N_s> |' 1291 * //' -alpha*Sum_a[ > | O-tilde |', 1292 1293 * ' | contribution | contribution |' 1294 * //' |' 1295 * //'<Rra>*{<Rra>+OmmS}] | contribution |', 1296 1297 * ' +-------------------------------------' 1298 * //'--------------------------------------' 1299 * //'----------------------------------+' 1300 WRITE(LUPRI,'(A,F20.16,A,F20.16,A,F20.16,A, 1301 & F20.16,A,F20.16,A)') 1302 * ' |',ECLVDW,' |',ENUQM3,' |',EELEL,' |',TEMP1,' |', 1303 * OTILDE,' |' 1304 END IF 1305 1306 CALL QEXIT('QM3_ENERGY') 1307 RETURN 1308 END 1309C************************************************************************* 1310C /* Deck qm3grad */ 1311 SUBROUTINE QM3GRAD(CREF,CMO,INDXCI,G,EQM3T,WRK,LWRK,IPRINT) 1312 1313#include "implicit.h" 1314#include "priunit.h" 1315#include "dummy.h" 1316#include "qm3.h" 1317#include "mxcent.h" 1318#include "maxash.h" 1319#include "maxorb.h" 1320#include "infinp.h" 1321#include "infvar.h" 1322#include "inforb.h" 1323#include "infind.h" 1324#include "inftap.h" 1325#include "infpri.h" 1326#include "inflin.h" 1327 1328 DIMENSION CREF(*),CMO(*),INDXCI(*) 1329 DIMENSION G(*),WRK(LWRK) 1330 LOGICAL FIRST, FNDLAB 1331 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0, DCVAL = D2 ) 1332 CHARACTER*8 STAR8, SOLVDI, EODATA 1333 SAVE FIRST 1334 DATA FIRST/.TRUE./, STAR8/'********'/ 1335 DATA SOLVDI/'SOLVDIAG'/, EODATA/'EODATA '/ 1336 1337 CALL QENTER('QM3GRAD') 1338 1339 KDIASH = 1 1340 KGRDLM = KDIASH + NVAR 1341 KUCMO = KGRDLM + NVARH 1342 KFSOLAO = KUCMO + NORBT*NBAST 1343 KFSOLMO = KFSOLAO + NNBASX 1344 KDIALM = KFSOLMO + NNORBX 1345 KDV = KDIALM + NVAR 1346 KDENS = KDV + N2BASX 1347 KDVS = KDENS + NNBASX 1348 KWRK = KDVS + NNBASX 1349 LWRK1 = LWRK - KWRK 1350 1351 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3GRAD',-KWRK,LWRK) 1352 1353C Construct the ao density matrix from the molecular orbial 1354C coefficients. 1355 1356 CALL FCKDEN((NISHT.GT.0),.FALSE.,WRK(KDV), 1357 * DUMMY,CMO,DUMMY,WRK(KWRK),LWRK1) 1358 1359 CALL DGEFSP(NBAST,WRK(KDV),WRK(KDVS)) 1360 CALL PKSYM1(WRK(KDVS),WRK(KDENS),NBAS,NSYM,1) 1361 1362 1363 CALL DZERO(WRK(KDIASH),NVAR) 1364 CALL DZERO(WRK(KGRDLM),NVARH) 1365 1366 EQM3T = 0.0D0 1367 ENSEL = 0.0D0 1368 EPOL = 0.0D0 1369 EELEL = 0.0D0 1370 1371C Make effective solvent operator and transform to mo basis 1372 CALL QM3FCK(WRK(KDENS),WRK(KDENS),WRK(KFSOLAO),EQM3T,ENSEL, 1373 & EPOL,EELEL,WRK(KWRK),LWRK1,IPRINT) 1374 CALL UPKCMO(CMO,WRK(KUCMO)) 1375 CALL UTHU(WRK(KFSOLAO),WRK(KFSOLMO),WRK(KUCMO),WRK(KWRK), 1376 & NBAST,NORBT) 1377 1378 IF (.NOT. OLDTG) EQM3T = EQM3T - EELEL 1379 1380C Calculate gradients 1381 IF (NWOPT .GT. 0) THEN 1382 CALL SOLGO(D2,WRK(KDENS),WRK(KFSOLMO),WRK(KGRDLM+NCONF)) 1383 ENDIF 1384 1385 IF (NWOPPT .GT. 0) THEN 1386 CALL OR1DIA(WRK(KDENS),WRK(KFSOLMO),WRK(KDIALM+NCONST), 1387 & IPRINT) 1388 END IF 1389 1390 DO 420 I = 0,(NVAR-1) 1391 WRK(KDIASH+I) = WRK(KDIASH+I) 1392 & + WRK(KDIALM+I) 1393 & + WRK(KGRDLM+I) * WRK(KGRDLM+I) 1394 420 CONTINUE 1395 1396 FAC=1.0D0 1397 CALL DAXPY(NVARH,FAC,WRK(KGRDLM),1,G,1) 1398 1399 IF (LUIT2 .GT. 0) THEN 1400 NW4 = MAX(NWOPT,4) 1401 REWIND LUIT2 1402 IF (FNDLAB(EODATA,LUIT2)) BACKSPACE LUIT2 1403 WRITE (LUIT2) STAR8,STAR8,STAR8,SOLVDI 1404 IF (NWOPT .GT. 0) CALL WRITT(LUIT2,NW4,WRK(KDIASH+NCONF)) 1405 WRITE (LUIT2) STAR8,STAR8,STAR8,EODATA 1406 END IF 1407 1408 FIRST = .FALSE. 1409 CALL QEXIT('QM3GRAD') 1410 RETURN 1411 END 1412C 1413C************************************************************************* 1414C /* Deck qm3lin */ 1415 SUBROUTINE QM3LIN(NOSIM,BOVECS,CREF,CMO,INDXCI, 1416 & SOVECS,ORBLIN,WRK,LWRK) 1417C 1418#include "implicit.h" 1419#include "inflin.h" 1420#include "inforb.h" 1421#include "dummy.h" 1422#include "priunit.h" 1423#include "infpri.h" 1424 1425 DIMENSION BOVECS(*),CREF(*),CMO(*),INDXCI(*) 1426 DIMENSION SOVECS(*),WRK(LWRK) 1427 LOGICAL ORBLIN 1428 1429 CALL QENTER('QM3LIN') 1430 1431 KDV = 1 1432 KDVS = KDV + N2BASX 1433 KDENS = KDVS + NNBASX 1434 KWRK = KDENS + NNBASX 1435 LWRK1 = LWRK - KWRK 1436 1437 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3LIN',-KWRK,LWRK) 1438 1439C Construct the ao density matrix from the molecular orbial 1440C coefficients. 1441 CALL FCKDEN((NISHT.GT.0),.FALSE.,WRK(KDV), 1442 * DUMMY,CMO,DUMMY,WRK(KWRK),LWRK1) 1443 1444 CALL DGEFSP(NBAST,WRK(KDV),WRK(KDVS)) 1445 CALL PKSYM1(WRK(KDVS),WRK(KDENS),NBAS,NSYM,1) 1446 1447C 1448 IF ( NOSIM .GT.0 ) THEN 1449 IF ( .NOT.ORBLIN ) THEN 1450 NSVAR = NVARPT 1451 ELSE 1452 NSVAR = NWOPPT 1453 END IF 1454 1455 1456 CALL QM3HESS(NOSIM,BOVECS,CREF,CMO,INDXCI, 1457 & WRK(KDENS),SOVECS,NSVAR,WRK(KWRK),LWRK1) 1458 1459 END IF 1460 CALL QEXIT('QM3LIN') 1461 RETURN 1462 END 1463C 1464C************************************************************************* 1465C /* Deck qm3hess */ 1466 SUBROUTINE QM3HESS(NOSIM,BOVECS,CREF,CMO,INDXCI, 1467 & DV,SVEC,NSVEC,WRK,LWRK) 1468C 1469#include "implicit.h" 1470#include "maxorb.h" 1471#include "inflin.h" 1472#include "inforb.h" 1473#include "infvar.h" 1474#include "maxash.h" 1475#include "infind.h" 1476#include "qm3.h" 1477#include "mxcent.h" 1478#include "priunit.h" 1479#include "dummy.h" 1480#include "inftap.h" 1481#include "orgcom.h" 1482#include "infinp.h" 1483#include "nuclei.h" 1484C 1485 DIMENSION BOVECS(NWOPPT,*),CREF(*),CMO(*),INDXCI(*) 1486 DIMENSION DV(*),SVEC(NSVEC,*),WRK(LWRK) 1487 LOGICAL FULHES 1488 PARAMETER ( D0 = 0.0D0, D2 = 2.0D0 ) 1489 1490 CHARACTER*8 LABINT(9*MXCENT) 1491 LOGICAL TOFILE, TRIMAT, EXP1VL 1492 DIMENSION INTREP(9*MXCENT), INTADR(9*MXCENT) 1493C 1494 CALL QENTER('QM3HESS') 1495 1496 CALL QUIT('QM3HESS needs further debuging. Sorry !') 1497C 1498 KRXYO = 1 1499 KUCMO = KRXYO + NOSIM*NNORBX 1500 KUBO = KUCMO + NORBT*NBAST 1501 KFSOLAO = KUBO + NOSIM*N2ORBX 1502 KFSOLMO = KFSOLAO + NNBASX 1503 KUFSOL = KFSOLMO + NNORBX 1504 KTGX = KUFSOL + N2ORBX 1505 KRRAOx = KTGX + NNORBX 1506 KRRAOy = KRRAOx + NNBASX 1507 KRRAOz = KRRAOy + NNBASX 1508 KTRAO = KRRAOz + NNBASX 1509 KTRMO = KTRAO + NNBASX 1510 KUTRX = KTRMO + NNORBX 1511 KUTR = KUTRX + N2ORBX 1512 KTRX = KUTR + N2ORBX 1513 KWRK = KTRX + NNORBX 1514 LWRK1 = LWRK - KWRK 1515 1516 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3HESS',-KWRK,LWRK) 1517 1518 IF (JWOPSY .NE. LSYMPT .OR. NWOPPT .NE. NWOPT) THEN 1519 WRITE (LUPRI,*) 'QM3HESS ERROR: JWOPSY .ne. LSYMPT .or.'// 1520 * ' NWOPPT .ne. NWOPT' 1521 WRITE (LUPRI,*) 'LSYMPT,JWOPSY:',LSYMPT,JWOPSY 1522 WRITE (LUPRI,*) 'NWOPPT,NWOPT :',NWOPPT,NWOPT 1523 CALL QUIT('QM3HESS ERROR,lsympt.ne.jwopsy or nwoppt.ne.nwopt') 1524 END IF 1525 1526C Determine if full Hessian or only orbital Hessian 1527 1528 FULHES = (NSVEC .EQ. NVARPT) 1529 IF (FULHES) THEN 1530 JSOVEC = 1 + NCONST 1531 ELSE 1532 JSOVEC = 1 1533 END IF 1534 1535 IF (IQM3PR .GT. 15) THEN 1536 WRITE (LUPRI,'(/A)') ' SOLLNO - svec(orb) on entry' 1537 CALL OUTPUT(SVEC(JSOVEC,1),1,NWOPPT,1,NOSIM, 1538 * NSVEC,NOSIM,1,LUPRI) 1539 END IF 1540 1541C 1. Construct effective operators Tgxo and RRayo 1542 1543 CALL DZERO(WRK(KRXYO),NOSIM*NNORBX ) 1544 CALL DZERO(WRK(KUBO),NOSIM*N2ORBX) 1545 1546C 1.a. Unpack symmetry blocked CMO 1547 CALL UPKCMO(CMO,WRK(KUCMO)) 1548C 1.b. Calculate unpacked orbital trial vectors in UBO 1549 1550 IF (NOSIM.GT.0) THEN 1551 DO 210 IOSIM = 1,NOSIM 1552 JUBO = KUBO + (IOSIM-1)*N2ORBX 1553 CALL UPKWOP(NWOPPT,JWOP,BOVECS(1,IOSIM),WRK(JUBO)) 1554 210 CONTINUE 1555 END IF 1556 1557 1558 CALL DZERO(WRK(KFSOLMO),NNORBX) 1559 CALL DZERO(WRK(KUFSOL),N2ORBX) 1560 1561C 1.c. Construct the Tg operator in ao and transform to mo 1562 CALL QM3FCKMO(CMO,WRK(KFSOLMO),WRK(KWRK),LWRK1,IPRINT) 1563 1564 CALL DSPTSI(NORBT,WRK(KFSOLMO),WRK(KUFSOL)) 1565 1566 DO 200 IOSIM = 1,NOSIM 1567 JRXYO = KRXYO + (IOSIM-1)*NNORBX 1568 JUBO = KUBO + (IOSIM-1)*N2ORBX 1569 1570C 1.d. Calculate one-index transformed Tg integrals 1571 1572 CALL DZERO(WRK(KUTRX),N2ORBX) 1573 CALL TR1UH1(WRK(JUBO),WRK(KUFSOL),WRK(KUTRX),1) 1574 CALL DGETSP(NORBT,WRK(KUTRX),WRK(KTGX)) 1575 1576 FAC = 1.0D0 1577 CALL DAXPY(NNORBX,FAC,WRK(KTGX),1,WRK(JRXYO),1) 1578 1579 1580 WRITE (LUPRI,'(/A,I3,A,I3)') 1581 * ' --- SOLLNO - Rxo matrix no.',IOSIM,' of',NOSIM 1582 CALL OUTPAK(WRK(JRXYO),NORBT,1,LUPRI) 1583 1584 1585 IF (.NOT. INTDIR) THEN 1586 CALL GPOPEN(LUQM3E,'ELFDMM','OLD',' ', 1587 & 'UNFORMATTED',IDUMMY,.FALSE.) 1588 REWIND(LUQM3E) 1589 ENDIF 1590 1591C 1.e. Readin Rra integrals and transform to mo 1592 1593 LM = 0 1594 1595 IF (INTDIR) THEN 1596 L = NUSITE + NUALIS(0) 1597 OBKPX = DIPORG(1) 1598 OBKPY = DIPORG(2) 1599 OBKPZ = DIPORG(3) 1600 ENDIF 1601 1602 DO 520 I = 1, ISYTP 1603 IF (MDLWRD(I)(1:5) .EQ. 'SPC_E') THEN 1604 DO 521 J = NSYSBG(I), NSYSED(I) 1605 DO 522 K = 1, NUALIS(I) 1606 LM = LM + 1 1607 1608 FAC = -ALPIMM(I,K) 1609 1610C Calculate one-index transformed Rra.x integrals 1611 1612 IF (INTDIR) THEN 1613 KMAT = KWRK 1614 KLAST = KMAT + 3*NNBASX 1615 LWRK2 = LWRK - KLAST + 1 1616 IATNOW = NUCIND + L + LM 1617 1618 CALL DZERO(WRK(KMAT),3*NNBASX) 1619 1620 KPATOM = 0 1621 NOSIMI = 3 1622 TOFILE = .FALSE. 1623 TRIMAT = .TRUE. 1624 EXP1VL = .FALSE. 1625 DIPORG(1) = CORD(1,IATNOW) 1626 DIPORG(2) = CORD(2,IATNOW) 1627 DIPORG(3) = CORD(3,IATNOW) 1628 1629 RUNQM3=.TRUE. 1630 CALL GET1IN(WRK(KMAT),'NEFIELD',NOSIMI,WRK(KLAST), 1631 & LWRK2,LABINT,INTREP,INTADR,IATNOW,TOFILE, 1632 & KPATOM,TRIMAT,DUMMY,EXP1VL,DUMMY,IQM3PR) 1633 RUNQM3=.FALSE. 1634 1635 IF (QMDAMP) THEN 1636 DIST = (DIPORG(1)-QMCOM(1))**2 + 1637 & (DIPORG(2)-QMCOM(2))**2 + 1638 & (DIPORG(3)-QMCOM(3))**2 1639 DIST = SQRT(DIST) 1640 DFACT = (1-exp(-ADAMP*DIST))**3 1641 CALL DSCAL(3*NNBASX,DFACT,WRK(KMAT),1) 1642 ENDIF 1643 1644 CALL DZERO(WRK(KTRMO),NNORBX) 1645 CALL UTHU(WRK(KMAT),WRK(KTRMO),WRK(KUCMO),WRK(KLAST), 1646 & NBAST,NORBT) 1647 ELSE 1648 CALL DZERO(WRK(KRRAOx),NNBASX) 1649 CALL READT(LUQM3E,NNBASX,WRK(KRRAOx)) 1650 CALL DZERO(WRK(KTRMO),NNORBX) 1651 CALL UTHU(WRK(KRRAOx),WRK(KTRMO),WRK(KUCMO), 1652 & WRK(KWRK),NBAST,NORBT) 1653 ENDIF 1654 1655 CALL DZERO(WRK(KUTR),N2ORBX) 1656 CALL DZERO(WRK(KUTRX),N2ORBX) 1657 CALL DZERO(WRK(KTRX),NNORBX) 1658 CALL DSPTSI(NORBT,WRK(KTRMO),WRK(KUTR)) 1659 CALL TR1UH1(WRK(JUBO),WRK(KUTR),WRK(KUTRX),1) 1660 CALL DGETSP(NORBT,WRK(KUTRX),WRK(KTRX)) 1661 1662 FACx = 0.0D0 1663 DO 523 IW = 1,NISHT 1664 IG = ISX(IW) 1665 FACx = FACx + 2.0D0*WRK(KTRX -1 + IROW(IG+1)) 1666 523 CONTINUE 1667 FACx = FACx * FAC 1668 1669 CALL DAXPY(NNORBX,FACx,WRK(KTRMO),1,WRK(JRXYO),1) 1670 1671C Calculate one-index transformed Rra.y integrals 1672 1673 IF (INTDIR) THEN 1674 CALL DZERO(WRK(KTRMO),NNORBX) 1675 CALL UTHU(WRK(KMAT+NNBASX),WRK(KTRMO),WRK(KUCMO), 1676 & WRK(KLAST),NBAST,NORBT) 1677 ELSE 1678 CALL DZERO(WRK(KRRAOy),NNBASX) 1679 CALL READT(LUQM3E,NNBASX,WRK(KRRAOy)) 1680 CALL DZERO(WRK(KTRMO),NNORBX) 1681 CALL UTHU(WRK(KRRAOy),WRK(KTRMO),WRK(KUCMO), 1682 & WRK(KWRK),NBAST,NORBT) 1683 ENDIF 1684 1685 CALL DZERO(WRK(KUTR),N2ORBX) 1686 CALL DZERO(WRK(KUTRX),N2ORBX) 1687 CALL DZERO(WRK(KTRX),NNORBX) 1688 CALL DSPTSI(NORBT,WRK(KTRMO),WRK(KUTR)) 1689 CALL TR1UH1(WRK(JUBO),WRK(KUTR),WRK(KUTRX),1) 1690 CALL DGETSP(NORBT,WRK(KUTRX),WRK(KTRX)) 1691 1692 FACy = 0.0D0 1693 DO 524 IW = 1,NISHT 1694 IG = ISX(IW) 1695 FACy = FACy + 2.0D0*WRK(KTRX -1 + IROW(IG+1)) 1696 524 CONTINUE 1697 FACy = FACy * FAC 1698 1699 CALL DAXPY(NNORBX,FACy,WRK(KTRMO),1,WRK(JRXYO),1) 1700 1701C Calculate one-index transformed Rra.z integrals 1702 1703 IF (INTDIR) THEN 1704 CALL DZERO(WRK(KTRMO),NNORBX) 1705 CALL UTHU(WRK(KMAT+2*NNBASX),WRK(KTRMO),WRK(KUCMO), 1706 & WRK(KLAST),NBAST,NORBT) 1707 ELSE 1708 CALL DZERO(WRK(KRRAOz),NNBASX) 1709 CALL READT(LUQM3E,NNBASX,WRK(KRRAOz)) 1710 CALL DZERO(WRK(KTRMO),NNORBX) 1711 CALL UTHU(WRK(KRRAOz),WRK(KTRMO),WRK(KUCMO), 1712 & WRK(KWRK),NBAST,NORBT) 1713 ENDIF 1714 1715 CALL DZERO(WRK(KUTR),N2ORBX) 1716 CALL DZERO(WRK(KUTRX),N2ORBX) 1717 CALL DZERO(WRK(KTRX),NNORBX) 1718 CALL DSPTSI(NORBT,WRK(KTRMO),WRK(KUTR)) 1719 CALL TR1UH1(WRK(JUBO),WRK(KUTR),WRK(KUTRX),1) 1720 CALL DGETSP(NORBT,WRK(KUTRX),WRK(KTRX)) 1721 1722 FACz = 0.0D0 1723 DO 525 IW = 1,NISHT 1724 IG = ISX(IW) 1725 FACz = FACz + 2.0D0*WRK(KTRX -1 + IROW(IG+1)) 1726 525 CONTINUE 1727 FACz = FACz * FAC 1728 1729 CALL DAXPY(NNORBX,FACz,WRK(KTRMO),1,WRK(JRXYO),1) 1730 1731 522 CONTINUE 1732 521 CONTINUE 1733 END IF 1734 520 CONTINUE 1735C 1736 IF (.NOT. INTDIR) CALL GPCLOSE(LUQM3E,'KEEP') 1737 1738 IF (INTDIR) THEN 1739 DIPORG(1) = OBKPX 1740 DIPORG(2) = OBKPY 1741 DIPORG(3) = OBKPZ 1742 END IF 1743 1744 200 CONTINUE 1745 1746 IF (IQM3PR .GT. 15) THEN 1747 DO 600 IOSIM = 1,NOSIM 1748 JRXYO = KRXYO + (IOSIM-1)*NNORBX 1749 WRITE (LUPRI,'(/A,I3,A,I3)') 1750 * ' SOLLNO - (Rxo + Ryo) matrix no.',IOSIM,' of',NOSIM 1751 CALL OUTPAK(WRK(JRXYO),NORBT,1,LUPRI) 1752 600 CONTINUE 1753 END IF 1754 1755C 2. Add effective operators Tgxo and RRayo contribution to SVEC(NSVEC,NOSIM) 1756C Tell SOLGO only to use the NWOPPT first JWOP entries 1757 DO 800 IOSIM = 1,NOSIM 1758 JRXYO = KRXYO + (IOSIM-1)*NNORBX 1759 CALL SOLGO(D2,DV,WRK(JRXYO),SVEC(JSOVEC,IOSIM)) 1760 800 CONTINUE 1761C 1762 IF (IQM3PR .GT. 15) THEN 1763 WRITE (LUPRI,'(/A)') ' SOLLNO - svec(orb) on exit' 1764 CALL OUTPUT(SVEC(JSOVEC,1),1,NWOPPT,1,NOSIM, 1765 * NSVEC,NOSIM,1,LUPRI) 1766 END IF 1767 1768 CALL QEXIT('QM3HESS') 1769 RETURN 1770 END 1771C************************************************************************* 1772C /* Deck qm3fck */ 1773 SUBROUTINE QM3FCKMO(CMO,FSOL,WRK,LWRK,IPRINT) 1774C 1775C Construct the QM/MM contribution to the Fock-matrix in MO basis 1776C 1777#include "implicit.h" 1778#include "priunit.h" 1779#include "dummy.h" 1780#include "qm3.h" 1781#include "inforb.h" 1782#include "infopt.h" 1783C 1784 DIMENSION CMO(*), FSOL(*), WRK(LWRK) 1785C 1786 CALL QENTER('QM3FCKMO') 1787C 1788 KDV = 1 1789 KDENS = KDV + N2BASX 1790 KDVS = KDENS + NNBASX 1791 KFSOLAO = KDVS + NNBASX 1792 KUCMO = KFSOLAO + NNBASX 1793 KWRK = KUCMO + NORBT*NBAST 1794 LWRK1 = LWRK - KWRK 1795 1796 IF (LWRK1 .LT. 0) CALL ERRWRK('QM3FCKMO',-KWRK,LWRK) 1797 1798C Construct the density matrix 1799 CALL FCKDEN((NISHT.GT.0),.FALSE.,WRK(KDV), 1800 * DUMMY,CMO,DUMMY,WRK(KWRK),LWRK1) 1801 1802 CALL DGEFSP(NBAST,WRK(KDV),WRK(KDVS)) 1803 CALL PKSYM1(WRK(KDVS),WRK(KDENS),NBAS,NSYM,1) 1804 1805C Construct the QM/MM contribution to the Fock-matrix in AO 1806 CALL QM3FCK(WRK(KDENS),WRK(KDENS),WRK(KFSOLAO),ESOLT, 1807 * ENSQM3,EPOQM3,EELEL,WRK(KWRK),LWRK1,IPRINT) 1808 1809C Transform to mo 1810 CALL UPKCMO(CMO,WRK(KUCMO)) 1811 CALL UTHU(WRK(KFSOLAO),FSOL,WRK(KUCMO),WRK(KWRK), 1812 & NBAST,NORBT) 1813C 1814 CALL QEXIT('QM3FCKMO') 1815 RETURN 1816 END 1817