1c -*-mode:fortran; fortran-continuation-string: "&" -*- 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 20C /* Deck getfd */ 21 SUBROUTINE GETFD(IATOM,NOORBI,TWOH1,NOH1,NOH2,NODC,NODV,IPRINT, 22 & DERIV,EMYD,CMO,DV,FCD,FVD,FDFTD,WORK,LWORK,AOMAT, 23 & DOREAL,INTPRI) 24C 25C 3-Nov-1988 hjaaj -- interface routine for GETFD 26C 27#include "implicit.h" 28 LOGICAL TWOH1, NOH1, NOH2, NODC, NODV, DERIV(3), 29 * NOORBI, AOMAT, DOREAL 30 DIMENSION CMO(*), DV(*), FCD(N2BASX,3), FVD(N2ORBX,3), 31 * EMYD(3), FDFTD(N2BASX,3), WORK(LWORK) 32#include "priunit.h" 33#include "iratdef.h" 34C 35C Used from common blocks: 36C INFORB : NNBASX,... 37C 38#include "inforb.h" 39C 40 CALL QENTER('GETFD') 41 KFCAOX = 1 42 KFCAOY = KFCAOX + NNBASX 43 KFCAOZ = KFCAOY + NNBASX 44 IF (NODV) THEN 45 KFVAOX = KFCAOZ + NNBASX 46 KFVAOY = KFVAOX 47 KFVAOZ = KFVAOX 48 KDVAO = KFVAOX 49 KDCAO = KFVAOX 50 ELSE 51 KFVAOX = KFCAOZ + NNBASX 52 KFVAOY = KFVAOX + NNBASX 53 KFVAOZ = KFVAOY + NNBASX 54 KDVAO = KFVAOZ + NNBASX 55 KDCAO = KDVAO + NNBASX 56 END IF 57 KDCFO = KDCAO + NNBASX 58 KIINDX4= KDCFO + NNBASX 59 KGETFD = KIINDX4+ 4*600/IRAT 60 LGETFD = LWORK - KGETFD + 1 61 IF (KGETFD.GT.LWORK) CALL STOPIT('GETFD',' ',KGETFD,LWORK) 62C 63 CALL GETFD1(IATOM,NOORBI,TWOH1,NOH1,NOH2,NODC,NODV,IPRINT, 64 & DERIV,EMYD,CMO,DV,FCD,FVD,FDFTD, 65 & WORK(KFCAOX),WORK(KFCAOY),WORK(KFCAOZ), 66 & WORK(KFVAOX),WORK(KFVAOY),WORK(KFVAOZ), 67 & WORK(KDCAO), WORK(KDVAO), WORK(KDCFO), 68 & WORK(KIINDX4),WORK(KGETFD),LGETFD,AOMAT,DOREAL, 69 & INTPRI) 70 CALL QEXIT('GETFD') 71 RETURN 72 END 73C /* Deck getfd1 */ 74 SUBROUTINE GETFD1(IATOM,NOORBI,TWOH1,NOH1,NOH2,NODC,NODV,IPRINT, 75 & DERIV,EMYD,CMO,DV,FCD,FVD,FDFTD,FCSOX,FCSOY, 76 & FCSOZ,FVSOX,FVSOY,FVSOZ,DCSO,DVSO,DCFOLD, 77 & IINDX4,WORK,LWORK,AOMAT,DOREAL,INTPRI) 78C 79C This subroutine calculates the AO-differentiated Fock (inactive 80C and active) in MO basis. 81C 82C tuh Sep 1988 83C 84#include "implicit.h" 85#include "priunit.h" 86#include "mxcent.h" 87#include "maxorb.h" 88#include "maxaqn.h" 89#include "iratdef.h" 90 PARAMETER (D0 = 0.00 D00, D1 = 1.0D0, DP5 = 0.50D00) 91 LOGICAL TWOH1, NOH1, NOH2, NODC, NODV, DERIV(3), DERX, DERY, DERZ, 92 * NOORBI, AOMAT, DOFV, DOREAL 93#include "inforb.h" 94 DIMENSION IINDX4(4,600), CMO(*), DV(*), 95 * FCD(N2BASX,3), FVD(N2ORBX,3), FDFTD(N2BASX,3), 96 * FCSOX(NNBASX), FCSOY(NNBASX), FCSOZ(NNBASX), 97 * FVSOX(NNBASX), FVSOY(NNBASX), FVSOZ(NNBASX), 98 * DCSO(NNBASX), DVSO(NNBASX), 99 * EMYD(3), DCFOLD(*), WORK(LWORK) 100#include "ccom.h" 101#include "abainf.h" 102#include "dftcom.h" 103#include "nuclei.h" 104#include "inftap.h" 105#include "eribuf.h" 106 107 DOFV = (NISHT .GT. 0) .AND. (NASHT .GT. 0) .AND. (.NOT.NOORBI) 108 DERX = DERIV(1) 109 DERY = DERIV(2) 110 DERZ = DERIV(3) 111 IF (IPRINT .GT. 05 .OR. AOMAT) THEN 112 CALL TITLER('Output from GETFD','*',103) 113 END IF 114 IF (IPRINT .GT. 05) THEN 115 WRITE (LUPRI, '(A,3L5)')' DERIV ', DERX, DERY, DERZ 116 WRITE (LUPRI, '(A,L5)') ' NOH1 ', NOH1 117 WRITE (LUPRI, '(A,L5)') ' NOH2 ', NOH2 118 WRITE (LUPRI, '(A,L5)') ' NODC ', NODC 119 WRITE (LUPRI, '(A,L5)') ' NODV ', NODV 120 WRITE (LUPRI, '(A,L5)') ' AOMAT ', AOMAT 121 END IF 122 CALL ERIBUF_INI 123 LBUF = 600 124C 125 LENINT4 = 2*LBUF + NIBUF*LBUF + 1 ! length in integer*4 units 126C 127C 128C *************************************************** 129C ***** Set up active and inactive one-electron ***** 130C ***** density matrices in SO basis ***** 131C *************************************************** 132C 133 CALL DCDVSO(CMO,DV,DCSO,DVSO,DCFOLD,IPRINT,NODC,NODV) 134C 135C ******************************************************* 136C ***** One-electron contributions to Fock matrices ***** 137C ******************************************************* 138C 139C Inactive FOCK matrix 140C 141 CALL DZERO(FCD,3*N2BASX) 142 IF (NOH1) THEN 143 EMY1X = D0 144 EMY1Y = D0 145 EMY1Z = D0 146 ELSE 147C 148C Due to changes in GETSD, UFCSOX/Y/Z are now square, K.Ruud, Nov.-93 149C Note that in order to save memory, we use FCD for UFCDSOX/Y/Z 150C 151 CALL ONEDRL('HMAT',FCD(1,1),FCD(1,2),FCD(1,3),IATOM,DERX, 152 & DERY,DERZ,WORK,LWORK,IPRINT,INTPRI) 153C 154C DFT contribution 155C 156 IF (DFTADD) THEN 157 CALL DZERO(FDFTD,3*N2BASX) 158 IF (IATOM.GT.0) THEN 159 CALL DFTHED(IATOM,FCD,FDFTD,WORK,LWORK,IPRINT) 160 ELSE IF (IATOM .EQ. 0 .AND. .NOT.NOLOND) THEN 161 CALL DFTBRHS(FCD,WORK,LWORK,IPRINT) 162 END IF 163 END IF 164C 165c IF (DFTADD .AND. IATOM .EQ. 0 .AND. .NOT.NOLOND) THEN 166c CALL DFTBRHS(FCD,WORK,LWORK,IPRINT) 167c END IF 168C 169 IF (IATOM .GT. 0) THEN ! nuclear derivative - symmetric matrix 170 CALL DSITSP(NBAST,FCD(1,1),FCSOX) 171 CALL DSITSP(NBAST,FCD(1,2),FCSOY) 172 CALL DSITSP(NBAST,FCD(1,3),FCSOZ) 173 IF (DERX) EMY1X = DDOT(NNBASX,DCFOLD,1,FCSOX,1) 174 IF (DERY) EMY1Y = DDOT(NNBASX,DCFOLD,1,FCSOY,1) 175 IF (DERZ) EMY1Z = DDOT(NNBASX,DCFOLD,1,FCSOZ,1) 176 ELSE ! magnetic field - antisymmetric matrix 177 CALL DGETAP(NBAST,FCD(1,1),FCSOX) 178 CALL DGETAP(NBAST,FCD(1,2),FCSOY) 179 CALL DGETAP(NBAST,FCD(1,3),FCSOZ) 180 END IF 181 END IF 182C 183C Active Fock matrix 184C 185 IF (DOFV) THEN 186 CALL DZERO(FVSOX,NNBASX) 187 CALL DZERO(FVSOY,NNBASX) 188 CALL DZERO(FVSOZ,NNBASX) 189 END IF 190C 191C ******************************************************* 192C ***** Two-electron contributions to Fock matrices ***** 193C ******************************************************* 194C 195 IF ((NISHT.GT.0) .AND. (.NOT.(NOH2.OR.NOORBI)) 196 & .AND. (.NOT.(IATOM.EQ.0 .AND. NOLOND))) THEN 197C 198C Field Derivatives 199C ================= 200C 201 IF (IATOM.EQ.0) THEN 202C 203 IF (FCKDDR) THEN 204 KREAD = 1 205 LLAST = KREAD + NNBASX 206 IF (LLAST.GT.LWORK) CALL STOPIT('GETFD1','1',LLAST,LWORK) 207 IREC = 0 208 IF (DOREAL) THEN 209 IREC = 3 210 IF (EXPFCK) IREC = 3*NUCIND 211 END IF 212 IF (DOFV) IREC = 2*IREC 213 LNGTH = IRAT*NNBASX 214 CALL READDX(LUDFCK,IREC+1,LNGTH,WORK(KREAD)) 215 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOX,1) 216 CALL READDX(LUDFCK,IREC+2,LNGTH,WORK(KREAD)) 217 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOY,1) 218 CALL READDX(LUDFCK,IREC+3,LNGTH,WORK(KREAD)) 219 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOZ,1) 220 IF (DOFV) THEN 221 CALL READDX(LUDFCK,IREC+4,LNGTH,FVSOX) 222 CALL READDX(LUDFCK,IREC+5,LNGTH,FVSOY) 223 CALL READDX(LUDFCK,IREC+6,LNGTH,FVSOZ) 224 END IF 225 ELSE 226 KFCSQ = 1 227 KDCSQ = KFCSQ + 3*NBAST*NBAST 228 LLAST = KDCSQ + NBAST*NBAST 229 IF (NASHT .GT. 0) THEN 230 KFVSQ = LLAST 231 KDVSQ = KFVSQ + 3*NBAST*NBAST 232 LLAST = KDVSQ + NBAST*NBAST 233 END IF 234 LVLAST = LLAST + 600*(IRAT + 1) + 1 235C 236 IF (LVLAST.GT.LWORK) 237 & CALL STOPIT('GETFD1',' ',LVLAST,LWORK) 238C 239 CALL FCKMG2(FCSOX,FVSOX,DCSO,DVSO,WORK(KFCSQ), 240 & WORK(KFVSQ),WORK(KDCSQ),WORK(KDVSQ),IINDX4, 241 & IPRINT,.TRUE.,WORK(LLAST)) 242 END IF 243C 244C Geometrical Derivatives 245C ======================= 246C 247 ELSE 248 IF (FCKDDR) THEN 249 IREC = 0 250 IF (EXPFCK) IREC = 3*(IATOM - 1) 251 IF (DOFV) IREC = 2*IREC 252 KREAD = 1 253 LLAST = KREAD + NNBASX 254 IF (LLAST.GT.LWORK) CALL STOPIT('GETFD1','2',LLAST,LWORK) 255 LNGTH = IRAT*NNBASX 256 CALL READDX(LUDFCK,IREC+1,LNGTH,WORK(KREAD)) 257 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOX,1) 258 CALL READDX(LUDFCK,IREC+2,LNGTH,WORK(KREAD)) 259 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOY,1) 260 CALL READDX(LUDFCK,IREC+3,LNGTH,WORK(KREAD)) 261 CALL DAXPY(NNBASX,D1,WORK(KREAD),1,FCSOZ,1) 262 IF (DOFV) THEN 263 CALL READDX(LUDFCK,IREC+4,LNGTH,FVSOX) 264 CALL READDX(LUDFCK,IREC+5,LNGTH,FVSOY) 265 CALL READDX(LUDFCK,IREC+6,LNGTH,FVSOZ) 266 END IF 267 ELSE 268 CALL REWSPL(LU2DER) 269 KINT = 1 270 KIINT = KINT + LBUF 271 150 CONTINUE 272 CALL READI4(LU2DER,LENINT4,WORK(KINT)) 273 CALL AOLAB4(WORK(KIINT),LBUF,NIBUF,NBITS,IINDX4,LENGTH) 274 IF (IPRINT .GE. 05) WRITE (LUPRI, '(/,A,I5)') 275 & ' Two-electron integral buffer read, LENGTH =', 276 & LENGTH 277 IF (LENGTH .GT. 0) THEN 278 CALL FCKTWO(FCSOX,FVSOX,DCSO,DVSO,WORK(KINT), 279 & IINDX4,LENGTH,IPRINT) 280 ELSE IF (LENGTH .LT. 0 ) THEN 281 GO TO 300 282 END IF 283 GO TO 150 284C 285 300 CONTINUE 286 END IF 287 END IF 288 END IF 289C 290C ******************************************************** 291C ***** Transform matrices to MO basis and calculate ***** 292C ***** inactive contributions to molecular gradient ***** 293C ******************************************************** 294C 295C x direction 296C =========== 297C 298 KFSOSQ = 1 299 KLAST = KFSOSQ + N2BASX 300 LFREE = LWORK - KLAST + 1 301 IF (KLAST .GT. LWORK) CALL STOPIT('GETFD1',' ',KLAST,LWORK) 302 IF (DERX) THEN 303 IF (IPRINT .GE. 5 .OR. AOMAT) CALL TITLER 304 * ('Differentiation in x direction',' ',200) 305 ISIM = 1 306 IF (IATOM .GT. 0) THEN 307 IF (AOMAT) THEN 308 CALL AROUND('Inactive Fock matrix diff. in x direction') 309 CALL TRSOAO(FCSOX,WORK,LWORK,NBAST,IATOM,IPRINT) 310 IF (DOFV) THEN 311 CALL AROUND('Active Fock matrix diff. in x direction') 312 CALL TRSOAO(FVSOX,WORK,LWORK,NBAST,IATOM,IPRINT) 313 END IF 314 END IF 315 EMYD(ISIM) = DP5*(EMY1X + DDOT(NNBASX,DCFOLD,1,FCSOX,1)) 316 END IF 317 IF (IATOM .EQ. 0) THEN 318 CALL DAPTGE(NBAST,FCSOX,WORK(KFSOSQ)) 319 ELSE 320 CALL DSPTSI(NBAST,FCSOX,WORK(KFSOSQ)) 321 END IF 322 CALL ONETRA(CMO,FCD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 323 & IATOM,1,IPRINT) 324 IF (DOFV) THEN 325 IF (IATOM .GT. 0) THEN 326 CALL DSPTSI(NBAST,FVSOX,WORK(KFSOSQ)) 327 ELSE 328 CALL DAPTGE(NBAST,FVSOX,WORK(KFSOSQ)) 329 END IF 330 CALL ONETRA(CMO,FVD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 331 & IATOM,1,IPRINT) 332 END IF 333 IF (IPRINT .GE. 5) THEN 334 WRITE (LUPRI,'(A,F24.10,/)') 335 * ' Inactive contribution to molecular gradient:',EMYD(ISIM) 336 CALL HEADER('Inactive SO Fock matrix',-1) 337 CALL OUTPAK(FCSOX,NBAST,1,LUPRI) 338 CALL HEADER('Inactive MO Fock matrix',-1) 339 CALL OUTPUT(FCD(1,ISIM),1,NORBT,1,NORBT,NBAST,NBAST,1,LUPRI) 340 IF (DOFV) THEN 341 CALL HEADER('Active SO Fock matrix',-1) 342 CALL OUTPAK(FVSOX,NBAST,1,LUPRI) 343 CALL HEADER('Active MO Fock matrix',-1) 344 CALL OUTPUT(FVD(1,ISIM),1,NORBT,1,NORBT,NORBT,NORBT,1, 345 & LUPRI) 346 END IF 347 END IF 348 END IF 349C 350C y direction 351C =========== 352C 353 IF (DERY) THEN 354 IF (IPRINT .GE. 5 .OR. AOMAT) CALL TITLER 355 * ('Differentiation in y direction',' ',200) 356 ISIM = 2 357 IF (IATOM .GT. 0) THEN 358 IF (AOMAT) THEN 359 CALL AROUND('Inactive Fock matrix diff. in y direction') 360 CALL TRSOAO(FCSOY,WORK,LWORK,NBAST,IATOM,IPRINT) 361 IF (DOFV) THEN 362 CALL AROUND('Active Fock matrix diff. in y direction') 363 CALL TRSOAO(FVSOY,WORK,LWORK,NBAST,IATOM,IPRINT) 364 END IF 365 END IF 366 EMYD(ISIM) = DP5*(EMY1Y + DDOT(NNBASX,DCFOLD,1,FCSOY,1)) 367 END IF 368 IF (IATOM .GT. 0) THEN 369 CALL DSPTSI(NBAST,FCSOY,WORK(KFSOSQ)) 370 ELSE 371 CALL DAPTGE(NBAST,FCSOY,WORK(KFSOSQ)) 372 END IF 373 CALL ONETRA(CMO,FCD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 374 & IATOM,2,IPRINT) 375 IF (DOFV) THEN 376 IF (IATOM .GT. 0) THEN 377 CALL DSPTSI(NBAST,FVSOY,WORK(KFSOSQ)) 378 ELSE 379 CALL DAPTGE(NBAST,FVSOY,WORK(KFSOSQ)) 380 END IF 381 CALL ONETRA(CMO,FVD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 382 & IATOM,2,IPRINT) 383 END IF 384 IF (IPRINT .GE. 5) THEN 385 WRITE (LUPRI,'(A,F24.10,/)') 386 * ' Inactive contribution to molecular gradient:',EMYD(ISIM) 387 CALL HEADER('Inactive SO Fock matrix',-1) 388 CALL OUTPAK(FCSOY,NBAST,1,LUPRI) 389 CALL HEADER('Inactive MO Fock matrix',-1) 390 CALL OUTPUT(FCD(1,ISIM),1,NORBT,1,NORBT,NBAST,NBAST,1,LUPRI) 391 IF (DOFV) THEN 392 CALL HEADER('Active SO Fock matrix',-1) 393 CALL OUTPAK(FVSOY,NBAST,1,LUPRI) 394 CALL HEADER('Active MO Fock matrix',-1) 395 CALL OUTPUT(FVD(1,ISIM),1,NORBT,1,NORBT,NORBT,NORBT,1, 396 & LUPRI) 397 END IF 398 END IF 399 END IF 400C 401C z direction 402C =========== 403C 404 IF (DERZ) THEN 405 IF (IPRINT .GE. 5 .OR. AOMAT) CALL TITLER 406 * ('Differentiation in z direction',' ',200) 407 ISIM = 3 408 IF (IATOM .GT. 0) THEN 409 IF (AOMAT) THEN 410 CALL AROUND('Inactive Fock matrix diff. in z direction') 411 CALL TRSOAO(FCSOZ,WORK,LWORK,NBAST,IATOM,IPRINT) 412 IF (DOFV) THEN 413 CALL AROUND('Active Fock matrix diff. in z direction') 414 CALL TRSOAO(FVSOZ,WORK,LWORK,NBAST,IATOM,IPRINT) 415 END IF 416 END IF 417 EMYD(ISIM) = DP5*(EMY1Z + DDOT(NNBASX,DCFOLD,1,FCSOZ,1)) 418 END IF 419 IF (IATOM .GT. 0) THEN 420 CALL DSPTSI(NBAST,FCSOZ,WORK(KFSOSQ)) 421 ELSE 422 CALL DAPTGE(NBAST,FCSOZ,WORK(KFSOSQ)) 423 END IF 424 CALL ONETRA(CMO,FCD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 425 & IATOM,3,IPRINT) 426 IF (DOFV) THEN 427 IF (IATOM .EQ. 0) THEN 428 CALL DAPTGE(NBAST,FVSOZ,WORK(KFSOSQ)) 429 ELSE 430 CALL DSPTSI(NBAST,FVSOZ,WORK(KFSOSQ)) 431 END IF 432 CALL ONETRA(CMO,FVD(1,ISIM),WORK(KFSOSQ),WORK(KLAST),LFREE, 433 & IATOM,3,IPRINT) 434 END IF 435 IF (IPRINT .GE. 5) THEN 436 WRITE (LUPRI,'(A,F24.10,/)') 437 * ' Inactive contribution to molecular gradient:',EMYD(ISIM) 438 CALL HEADER('Inactive SO Fock matrix',-1) 439 CALL OUTPAK(FCSOZ,NBAST,1,LUPRI) 440 CALL HEADER('Inactive MO Fock matrix',-1) 441 CALL OUTPUT(FCD(1,ISIM),1,NORBT,1,NORBT,NBAST,NBAST,1,LUPRI) 442 IF (DOFV) THEN 443 CALL HEADER('Active SO Fock matrix',-1) 444 CALL OUTPAK(FVSOZ,NBAST,1,LUPRI) 445 CALL HEADER('Active MO Fock matrix',-1) 446 CALL OUTPUT(FVD(1,ISIM),1,NORBT,1,NORBT,NORBT,NORBT,1, 447 & LUPRI) 448 END IF 449 END IF 450 END IF 451 RETURN 452 END 453C /* Deck dcdvso */ 454 SUBROUTINE DCDVSO(CMO,DV,DCSO,DVSO,DCFOLD,IPRINT,NODC,NODV) 455C 456C This subroutine calculates the inactive and active one-electron 457C density matrices in SO basis (contravariant). 458C 459C tuh 310888 460C 461#include "implicit.h" 462#include "priunit.h" 463#include "maxorb.h" 464 PARAMETER (D0 = 0.0D0) 465 LOGICAL NODC, NODV 466 INTEGER R, S, U, V, UV, UR, US, VR, VS, RS 467 DIMENSION CMO(*), DV(*), DCSO(*), DVSO(*), DCFOLD(*) 468#include "inforb.h" 469C 470C ***** Print Section ***** 471C 472 IF (IPRINT .GT. 03) THEN 473 WRITE (LUPRI, '(//,A)') ' ----- SUBROUTINE DCDVSO ----- ' 474 WRITE (LUPRI, '(A,L5)') ' NODC ', NODC 475 WRITE (LUPRI, '(A,L5)') ' NODV ', NODV 476 WRITE (LUPRI, '(A,8I5)') ' NISH ', (NISH(I),I = 1,NSYM) 477 WRITE (LUPRI, '(A,8I5)') ' NASH ', (NASH(I),I = 1,NSYM) 478 WRITE (LUPRI, '(A,8I5)') ' NOCC ', (NOCC(I),I = 1,NSYM) 479 WRITE (LUPRI, '(A,8I5)') ' NORB ', (NORB(I),I = 1,NSYM) 480 WRITE (LUPRI, '(A,8I5)') ' NBAS ', (NBAS(I),I = 1,NSYM) 481 IF (IPRINT .GE. 05) THEN 482 CALL HEADER('Occupied molecular orbitals',0) 483 IEND = 0 484 DO 1000 ISYM = 1,NSYM 485 IF (NBAS(ISYM) .EQ. 0) GOTO 1000 486 IF (NOCC(ISYM) .EQ. 0) GOTO 1100 487 WRITE (LUPRI, '(//,A,I5,/)') ' Symmetry ', ISYM 488 IENDI = 0 489 DO 1200 I = 1, NOCC(ISYM) 490 WRITE (LUPRI, '(/,A,I5,/)') 491 * ' Molecular orbital ', I 492 WRITE (LUPRI, '(6F12.6)') 493 * (CMO(IEND+IENDI+J), J = 1, NBAS(ISYM)) 494 IENDI = IENDI + NBAS(ISYM) 4951200 CONTINUE 4961100 CONTINUE 497 IEND = IEND + NORB(ISYM)*NBAS(ISYM) 4981000 CONTINUE 499 CALL HEADER('Active density matrix (MO basis)',-1) 500 CALL OUTPAK(DV,NASHT,1,LUPRI) 501 END IF 502 END IF 503C 504C ***** Construct contravariant SO matrices ***** 505C 506 CALL DZERO(DCSO,NNBASX) 507 CALL DZERO(DCFOLD,NNBASX) 508 CALL DZERO(DVSO,NNBASX) 509 ISEND = 0 510 ICEND = 0 511 DO 110 ISYM = 1,NSYM 512 IOFF = IBAS(ISYM) 513 DO 100 R = 1, NBAS(ISYM) 514 DO 200 S = 1,R 515C 516C (I) Inactive contribution DCRS 517C 518 DCRS = D0 519 ICENDI = 0 520 DO 300 I = 1, NISH(ISYM) 521 DCRS = DCRS + CMO(ICEND+ICENDI+R)*CMO(ICEND+ICENDI+S) 522 ICENDI = ICENDI + NBAS(ISYM) 523 300 CONTINUE 524 DCRS = DCRS + DCRS 525 IF (NODC) DCRS = D0 526C 527C (II) Active contribution DVRS 528C 529 DVRS = D0 530 IF (.NOT. NODV) THEN 531 IASHI = IASH(ISYM) 532 UV = ((IASHI + 1)*(IASHI + 2))/2 533 IDVEND = ICEND + NBAS(ISYM)*NISH(ISYM) 534 ICENDU = IDVEND 535 DO 400 U = 1,NASH(ISYM) 536 ICENDV = IDVEND 537 DO 410 V = 1, U 538 DUV = DV(UV) 539 IF (ABS(DUV) .GT. D0) THEN 540 TMP = CMO(ICENDU+R)*CMO(ICENDV+S) 541 IF (U .NE. V) TMP = TMP 542 * + CMO(ICENDU+S)*CMO(ICENDV+R) 543 DVRS = DVRS + DUV*TMP 544 END IF 545 UV = UV + 1 546 ICENDV = ICENDV + NBAS(ISYM) 547 410 CONTINUE 548 UV = UV + IASHI 549 ICENDU = ICENDU + NBAS(ISYM) 550 400 CONTINUE 551 END IF 552 RS = (IOFF + R)*(IOFF + R - 1)/2 + IOFF + S 553 DCSO(RS) = DCRS 554 DVSO(RS) = DVRS 555 IF (R .EQ. S) THEN 556 DCFOLD(RS) = DCRS 557 ELSE 558 DCFOLD(RS) = DCRS + DCRS 559 END IF 560 200 CONTINUE 561 100 CONTINUE 562 ISEND = ISEND + (NBAS(ISYM)*(NBAS(ISYM) + 1))/2 563 ICEND = ICEND + NORB(ISYM)*NBAS(ISYM) 564110 CONTINUE 565C 566C ***** Print Section ***** 567C 568 IF (IPRINT .GE. 10) THEN 569 CALL HEADER('Inactive SO density matrix',-1) 570 CALL OUTPAK(DCSO,NBAST,1,LUPRI) 571 CALL HEADER('Inactive SO density matrix (folded)',-1) 572 CALL OUTPAK(DCFOLD,NBAST,1,LUPRI) 573 IF (.NOT. NODV) THEN 574 CALL HEADER('Active SO density matrix',-1) 575 CALL OUTPAK(DVSO,NBAST,1,LUPRI) 576 END IF 577 END IF 578 RETURN 579 END 580C /* Deck fcktwo */ 581 SUBROUTINE FCKTWO(FC,FV,DC,DV,BUF,IINDX4,LENGTH,IPRINT) 582C 583C This subroutine adds derivative two-electron integrals to the 584C active and inactive Fock matrices. 585C 586C The 14 classes of integral labels have been taken from Duke, 587C Chem. Phys. Lett. 13 (1972) 76 588C 589C Note: Canonical ordering of indices assumed 590C 591C tuh 120988 592C 593#include "implicit.h" 594#include "priunit.h" 595#include "maxorb.h" 596 PARAMETER (TWO = 2.00 D00, THRHLF = 1.50 D00, 597 * HALF = 0.50D00, IBIT08 = 255) 598 INTEGER P, Q, R, S 599 LOGICAL FIRST, ACTIVE 600#include "inforb.h" 601 DIMENSION BUF(600), IINDX4(4,600), ITRI(MXCORB), 602 * FC(3*NNBASX), FV(3*NNBASX), 603 * DC(NNBASX), DV(NNBASX) 604 SAVE FIRST, ITRI, ACTIVE, IOFF 605 DATA FIRST /.TRUE./ 606 607C 608 IF (FIRST) THEN 609 DO 100 I = 1, NBAST 610 ITRI(I) = I*(I - 1)/2 611 100 CONTINUE 612 ACTIVE = NASHT .GT. 0 613 FIRST = .FALSE. 614 END IF 615C 616 DO 200 I = 1, LENGTH 617 S = IINDX4(4,I) 618 IF (S .EQ. 0) THEN 619 IOFF = (IINDX4(3,I) - 1)*NNBASX 620 IF (IPRINT .GE. 25) THEN 621 ICOOR = IINDX4(3,I) 622 IREPE = IINDX4(2,I) 623 WRITE (LUPRI,'(A,2I5)') ' ICOOR, IREPE', ICOOR,IREPE 624 END IF 625 ELSE 626 DINT = BUF(I) 627 P = IINDX4(1,I) 628 Q = IINDX4(2,I) 629 R = IINDX4(3,I) 630 IF (IPRINT .GE. 25) WRITE (LUPRI,'(1X,F12.6,5X,4I5)') 631 * DINT,P,Q,R,S 632 IF (Q .LT. R) THEN 633 IF (Q .LT. S) THEN 634 IF (R .EQ. S) THEN 635C 636C p > r = s > q 637C 638 II = ITRI(R) + R 639 IK = ITRI(P) + R 640 IL = ITRI(R) + Q 641 KL = ITRI(P) + Q 642 FC(II+IOFF) = FC(II+IOFF) + TWO*DC(KL)*DINT 643 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(IK)*DINT 644 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(IL)*DINT 645 FC(KL+IOFF) = FC(KL+IOFF) + DC(II)*DINT 646 IF (ACTIVE) THEN 647 FV(II+IOFF) = FV(II+IOFF) + TWO*DV(KL)*DINT 648 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(IK)*DINT 649 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(IL)*DINT 650 FV(KL+IOFF) = FV(KL+IOFF) + DV(II)*DINT 651 END IF 652 ELSE 653C 654C p > r > s > q 655C 656 IJ = ITRI(P) + Q 657 IK = ITRI(P) + R 658 IL = ITRI(P) + S 659 JK = ITRI(R) + Q 660 JL = ITRI(S) + Q 661 KL = ITRI(R) + S 662 FC(IJ+IOFF) = FC(IJ+IOFF) + TWO*DC(KL)*DINT 663 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(JL)*DINT 664 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JK)*DINT 665 FC(JK+IOFF) = FC(JK+IOFF) - HALF*DC(IL)*DINT 666 FC(JL+IOFF) = FC(JL+IOFF) - HALF*DC(IK)*DINT 667 FC(KL+IOFF) = FC(KL+IOFF) + TWO*DC(IJ)*DINT 668 IF (ACTIVE) THEN 669 FV(IJ+IOFF) = FV(IJ+IOFF) + TWO*DV(KL)*DINT 670 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(JL)*DINT 671 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JK)*DINT 672 FV(JK+IOFF) = FV(JK+IOFF) - HALF*DV(IL)*DINT 673 FV(JL+IOFF) = FV(JL+IOFF) - HALF*DV(IK)*DINT 674 FV(KL+IOFF) = FV(KL+IOFF) + TWO*DV(IJ)*DINT 675 END IF 676 END IF 677 ELSE IF (Q .EQ. S) THEN 678 IF (P .EQ. R) THEN 679C 680C p = r > q = s 681C 682 II = ITRI(P) + P 683 IJ = ITRI(P) + Q 684 JJ = ITRI(Q) + Q 685 FC(II+IOFF) = FC(II+IOFF) - HALF*DC(JJ)*DINT 686 FC(IJ+IOFF) = FC(IJ+IOFF) + THRHLF*DC(IJ)*DINT 687 FC(JJ+IOFF) = FC(JJ+IOFF) - HALF*DC(II)*DINT 688 IF (ACTIVE) THEN 689 FV(II+IOFF) = FV(II+IOFF) - HALF*DV(JJ)*DINT 690 FV(IJ+IOFF) = FV(IJ+IOFF) + THRHLF*DV(IJ)*DINT 691 FV(JJ+IOFF) = FV(JJ+IOFF) - HALF*DV(II)*DINT 692 END IF 693 ELSE 694C 695C p > r > q = s 696C 697 IJ = ITRI(P) + Q 698 IL = ITRI(P) + R 699 JJ = ITRI(Q) + Q 700 JL = ITRI(R) + Q 701 FC(IJ+IOFF) = FC(IJ+IOFF) + THRHLF*DC(JL)*DINT 702 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JJ)*DINT 703 FC(JJ+IOFF) = FC(JJ+IOFF) - DC(IL)*DINT 704 FC(JL+IOFF) = FC(JL+IOFF) + THRHLF*DC(IJ)*DINT 705 IF (ACTIVE) THEN 706 FV(IJ+IOFF) = FV(IJ+IOFF) + THRHLF*DV(JL)*DINT 707 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JJ)*DINT 708 FV(JJ+IOFF) = FV(JJ+IOFF) - DV(IL)*DINT 709 FV(JL+IOFF) = FV(JL+IOFF) + THRHLF*DV(IJ)*DINT 710 END IF 711 END IF 712 ELSE 713 IF (P .EQ. R) THEN 714C 715C p = r > q > s 716C 717 IJ = ITRI(P) + Q 718 IL = ITRI(Q) + S 719 JJ = ITRI(P) + P 720 JL = ITRI(P) + S 721 FC(IJ+IOFF) = FC(IJ+IOFF) + THRHLF*DC(JL)*DINT 722 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JJ)*DINT 723 FC(JJ+IOFF) = FC(JJ+IOFF) - DC(IL)*DINT 724 FC(JL+IOFF) = FC(JL+IOFF) + THRHLF*DC(IJ)*DINT 725 IF (ACTIVE) THEN 726 FV(IJ+IOFF) = FV(IJ+IOFF) + THRHLF*DV(JL)*DINT 727 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JJ)*DINT 728 FV(JJ+IOFF) = FV(JJ+IOFF) - DV(IL)*DINT 729 FV(JL+IOFF) = FV(JL+IOFF) + THRHLF*DV(IJ)*DINT 730 END IF 731 ELSE 732C 733C p > r > q > s 734C 735 IJ = ITRI(P) + Q 736 IK = ITRI(P) + R 737 IL = ITRI(P) + S 738 JK = ITRI(R) + Q 739 JL = ITRI(Q) + S 740 KL = ITRI(R) + S 741 FC(IJ+IOFF) = FC(IJ+IOFF) + TWO*DC(KL)*DINT 742 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(JL)*DINT 743 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JK)*DINT 744 FC(JK+IOFF) = FC(JK+IOFF) - HALF*DC(IL)*DINT 745 FC(JL+IOFF) = FC(JL+IOFF) - HALF*DC(IK)*DINT 746 FC(KL+IOFF) = FC(KL+IOFF) + TWO*DC(IJ)*DINT 747 IF (ACTIVE) THEN 748 FV(IJ+IOFF) = FV(IJ+IOFF) + TWO*DV(KL)*DINT 749 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(JL)*DINT 750 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JK)*DINT 751 FV(JK+IOFF) = FV(JK+IOFF) - HALF*DV(IL)*DINT 752 FV(JL+IOFF) = FV(JL+IOFF) - HALF*DV(IK)*DINT 753 FV(KL+IOFF) = FV(KL+IOFF) + TWO*DV(IJ)*DINT 754 END IF 755 END IF 756 END IF 757 ELSE IF (Q .EQ. R) THEN 758 IF (P .EQ. Q) THEN 759 IF (R .EQ. S) THEN 760C 761C p = q = r = s 762C 763 II = ITRI(P) + P 764 FC(II+IOFF) = FC(II+IOFF) + HALF*DC(II)*DINT 765 FV(II+IOFF) = FV(II+IOFF) + HALF*DV(II)*DINT 766 ELSE 767C 768C p = q = r > s 769C 770 II = ITRI(P) + P 771 IL = ITRI(P) + S 772 FC(II+IOFF) = FC(II+IOFF) + DC(IL)*DINT 773 FC(IL+IOFF) = FC(IL+IOFF) + HALF*DC(II)*DINT 774 IF (ACTIVE) THEN 775 FV(II+IOFF) = FV(II+IOFF) + DV(IL)*DINT 776 FV(IL+IOFF) = FV(IL+IOFF) + HALF*DV(II)*DINT 777 END IF 778 END IF 779 ELSE 780 IF (R .EQ. S) THEN 781C 782C p > q = r = s 783C 784 II = ITRI(Q) + Q 785 IL = ITRI(P) + Q 786 FC(II+IOFF) = FC(II+IOFF) + DC(IL)*DINT 787 FC(IL+IOFF) = FC(IL+IOFF) + HALF*DC(II)*DINT 788 IF (ACTIVE) THEN 789 FV(II+IOFF) = FV(II+IOFF) + DV(IL)*DINT 790 FV(IL+IOFF) = FV(IL+IOFF) + HALF*DV(II)*DINT 791 END IF 792 ELSE 793C 794C p > q = r > s 795C 796 IJ = ITRI(P) + Q 797 IL = ITRI(P) + S 798 JJ = ITRI(Q) + Q 799 JL = ITRI(Q) + S 800 FC(IJ+IOFF) = FC(IJ+IOFF) + THRHLF*DC(JL)*DINT 801 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JJ)*DINT 802 FC(JJ+IOFF) = FC(JJ+IOFF) - DC(IL)*DINT 803 FC(JL+IOFF) = FC(JL+IOFF) + THRHLF*DC(IJ)*DINT 804 IF (ACTIVE) THEN 805 FV(IJ+IOFF) = FV(IJ+IOFF) + THRHLF*DV(JL)*DINT 806 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JJ)*DINT 807 FV(JJ+IOFF) = FV(JJ+IOFF) - DV(IL)*DINT 808 FV(JL+IOFF) = FV(JL+IOFF) + THRHLF*DV(IJ)*DINT 809 END IF 810 END IF 811 END IF 812 ELSE 813 IF (P .EQ. Q) THEN 814 IF (R .EQ. S) THEN 815C 816C p = q > r = s 817C 818 II = ITRI(P) + P 819 IK = ITRI(P) + R 820 KK = ITRI(R) + R 821 FC(II+IOFF) = FC(II+IOFF) + DC(KK)*DINT 822 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(IK)*DINT 823 FC(KK+IOFF) = FC(KK+IOFF) + DC(II)*DINT 824 IF (ACTIVE) THEN 825 FV(II+IOFF) = FV(II+IOFF) + DV(KK)*DINT 826 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(IK)*DINT 827 FV(KK+IOFF) = FV(KK+IOFF) + DV(II)*DINT 828 END IF 829 ELSE 830C 831C p = q > r > s 832C 833 II = ITRI(P) + P 834 IK = ITRI(P) + R 835 IL = ITRI(P) + S 836 KL = ITRI(R) + S 837 FC(II+IOFF) = FC(II+IOFF) + TWO*DC(KL)*DINT 838 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(IK)*DINT 839 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(IL)*DINT 840 FC(KL+IOFF) = FC(KL+IOFF) + DC(II)*DINT 841 IF (ACTIVE) THEN 842 FV(II+IOFF) = FV(II+IOFF) + TWO*DV(KL)*DINT 843 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(IK)*DINT 844 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(IL)*DINT 845 FV(KL+IOFF) = FV(KL+IOFF) + DV(II)*DINT 846 END IF 847 END IF 848 ELSE 849 IF (R .EQ. S) THEN 850C 851C p > q > r = s 852C 853 II = ITRI(R) + R 854 IK = ITRI(P) + R 855 IL = ITRI(Q) + R 856 KL = ITRI(P) + Q 857 FC(II+IOFF) = FC(II+IOFF) + TWO*DC(KL)*DINT 858 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(IK)*DINT 859 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(IL)*DINT 860 FC(KL+IOFF) = FC(KL+IOFF) + DC(II)*DINT 861 IF (ACTIVE) THEN 862 FV(II+IOFF) = FV(II+IOFF) + TWO*DV(KL)*DINT 863 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(IK)*DINT 864 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(IL)*DINT 865 FV(KL+IOFF) = FV(KL+IOFF) + DV(II)*DINT 866 END IF 867 ELSE 868C 869C p > q > r > s 870C 871 IJ = ITRI(P) + Q 872 IK = ITRI(P) + R 873 IL = ITRI(P) + S 874 JK = ITRI(Q) + R 875 JL = ITRI(Q) + S 876 KL = ITRI(R) + S 877 FC(IJ+IOFF) = FC(IJ+IOFF) + TWO*DC(KL)*DINT 878 FC(IK+IOFF) = FC(IK+IOFF) - HALF*DC(JL)*DINT 879 FC(IL+IOFF) = FC(IL+IOFF) - HALF*DC(JK)*DINT 880 FC(JK+IOFF) = FC(JK+IOFF) - HALF*DC(IL)*DINT 881 FC(JL+IOFF) = FC(JL+IOFF) - HALF*DC(IK)*DINT 882 FC(KL+IOFF) = FC(KL+IOFF) + TWO*DC(IJ)*DINT 883 IF (ACTIVE) THEN 884 FV(IJ+IOFF) = FV(IJ+IOFF) + TWO*DV(KL)*DINT 885 FV(IK+IOFF) = FV(IK+IOFF) - HALF*DV(JL)*DINT 886 FV(IL+IOFF) = FV(IL+IOFF) - HALF*DV(JK)*DINT 887 FV(JK+IOFF) = FV(JK+IOFF) - HALF*DV(IL)*DINT 888 FV(JL+IOFF) = FV(JL+IOFF) - HALF*DV(IK)*DINT 889 FV(KL+IOFF) = FV(KL+IOFF) + TWO*DV(IJ)*DINT 890 END IF 891 END IF 892 END IF 893 END IF 894 END IF 895 200 CONTINUE 896 RETURN 897 END 898C /* Deck fckmg2 */ 899 SUBROUTINE FCKMG2(FC,FV,DC,DV,FCSQ,FVSQ,DCSQ,DVSQ,IINDX4, 900 & IPRINT,SYMMET,WORK) 901C 902C ================================================================ 903C | This subroutine adds magnetic-field derivative two-electron | 904C | integrals to the active and inactive Fock matrices. | 905C ================================================================ 906C 907#include "implicit.h" 908#include "priunit.h" 909#include "iratdef.h" 910#include "maxorb.h" 911#include "mxcent.h" 912 PARAMETER (DP5 = 0.50D00, IBIT08 = 255) 913 LOGICAL SYMMET 914 INTEGER P, Q, R, S, X 915 DIMENSION FC(3*NNBASX), FV(3*NNBASX), DC(NNBASX), DV(NNBASX), 916 & FCSQ(NBAST,NBAST,3), FVSQ(NBAST,NBAST,3), 917 & DCSQ(NBAST,NBAST), DVSQ(NBAST,NBAST), 918 & IINDX4(4,600), WORK(*) 919#include "dftcom.h" 920#include "nuclei.h" 921#include "inftap.h" 922#include "eribuf.h" 923#include "inforb.h" 924 925C 926 LBUF = 600 927 CALL ERIBUF_INI 928C 929 LENINT4 = 2*LBUF + NIBUF*LBUF + 1 ! length in integer*4 units 930 KINT = 1 931 KIINT = KINT + LBUF 932C 933 CALL DZERO(FCSQ,3*N2BASX) 934 IF (NASHT .GE. 1) CALL DZERO(FVSQ,3*N2BASX) 935C 936C Form full (squared) active and inactive density matrices 937C ======================================================== 938C 939 IF (SYMMET) THEN 940 CALL DSPTSI(NBAST,DC,DCSQ) 941 IF (NASHT .GT. 0) THEN 942 CALL DSPTSI(NBAST,DV,DVSQ) 943 END IF 944 END IF 945C 946C Read integrals into buffers and process the integrals 947C ===================================================== 948C 949 CALL REWSPL(LU2DER) 950 CALL MOLLAB('AO2MGINT',LU2DER,LUPRI) 951C 952 300 CONTINUE 953 CALL READI4(LU2DER,LENINT4,WORK(KINT)) 954 CALL AOLAB4(WORK(KIINT),LBUF,NIBUF,NBITS,IINDX4,LENGTH) 955 IF (IPRINT .GE. 05) WRITE (LUPRI, '(/,A,I5)') 956 * ' Two-electron integral buffer read, LENGTH =', LENGTH 957 IF (LENGTH .GT. 0) THEN 958C 959 DO 400 I = 1, LENGTH 960 S = IINDX4(4,I) 961 IF (S .EQ. 0) THEN 962 X = IINDX4(3,I) 963 IF (IPRINT .GE. 25) THEN 964 IREPE = IINDX4(2,I) 965 WRITE (LUPRI,'(A,2I5)') ' ICOOR, IREPE', X,IREPE 966 END IF 967 ELSE 968 CINT = WORK(KINT + I - 1) 969 P = IINDX4(1,I) 970 Q = IINDX4(2,I) 971 R = IINDX4(3,I) 972C IF (IPRINT .GE. 25) WRITE (LUPRI,'(1X,F12.6,5X,4I5)') 973C & CINT,P,Q,R,S 974C 975 IF (P.EQ.R .AND. S.EQ.Q) CINT = DP5*CINT 976 EINT = DP5*HFXFAC*CINT 977C 978C Add Coulomb contributions to inactive fock matrix 979C ================================================= 980C 981 FCSQ(P,Q,X) = FCSQ(P,Q,X) - CINT*DCSQ(R,S) 982 FCSQ(Q,P,X) = FCSQ(Q,P,X) + CINT*DCSQ(S,R) 983 FCSQ(R,S,X) = FCSQ(R,S,X) - CINT*DCSQ(P,Q) 984 FCSQ(S,R,X) = FCSQ(S,R,X) + CINT*DCSQ(Q,P) 985C 986C Add Exchange contributions to inactive fock matrix 987C ================================================== 988C 989 FCSQ(P,S,X) = FCSQ(P,S,X) + EINT*DCSQ(Q,R) 990 FCSQ(S,P,X) = FCSQ(S,P,X) - EINT*DCSQ(R,Q) 991 FCSQ(R,Q,X) = FCSQ(R,Q,X) + EINT*DCSQ(S,P) 992 FCSQ(Q,R,X) = FCSQ(Q,R,X) - EINT*DCSQ(P,S) 993C 994 IF (NASHT.GE.1) THEN 995C 996C Add Coulomb contributions to active fock matrix 997C =============================================== 998C 999 FVSQ(P,Q,X) = FVSQ(P,Q,X) - CINT*DVSQ(R,S) 1000 FVSQ(Q,P,X) = FVSQ(Q,P,X) + CINT*DVSQ(S,R) 1001 FVSQ(R,S,X) = FVSQ(R,S,X) - CINT*DVSQ(P,Q) 1002 FVSQ(S,R,X) = FVSQ(S,R,X) + CINT*DVSQ(Q,P) 1003C 1004C Add Exchange contributions to active fock matrix 1005C ================================================ 1006C 1007 FVSQ(P,S,X) = FVSQ(P,S,X) + EINT*DVSQ(Q,R) 1008 FVSQ(S,P,X) = FVSQ(S,P,X) - EINT*DVSQ(R,Q) 1009 FVSQ(R,Q,X) = FVSQ(R,Q,X) + EINT*DVSQ(S,P) 1010 FVSQ(Q,R,X) = FVSQ(Q,R,X) - EINT*DVSQ(P,S) 1011C 1012 END IF 1013C 1014 END IF 1015 400 CONTINUE 1016 ELSE IF (LENGTH .LT. 0 ) THEN 1017 GO TO 500 1018 END IF 1019 GO TO 300 1020C 1021C Pack active and inactive fock matrices in lower triangular form 1022C =============================================================== 1023C 1024 500 CONTINUE 1025C 1026 IF (SYMMET) THEN 1027 IPCK = 0 1028 IF (NASHT .GT. 0) THEN 1029 DO 800 X = 1,3 1030 DO 700 J = 1, NBAST 1031 DO 600 I = 1, J 1032 IPCK = IPCK + 1 1033 FC(IPCK) = FC(IPCK) + FCSQ(I,J,X) 1034 FV(IPCK) = FV(IPCK) + FVSQ(I,J,X) 1035 600 CONTINUE 1036 700 CONTINUE 1037 800 CONTINUE 1038 ELSE 1039 DO 810 X = 1, 3 1040 DO 710 J = 1, NBAST 1041 DO 610 I = 1, J 1042 IPCK = IPCK + 1 1043 FC(IPCK) = FC(IPCK) + FCSQ(I,J,X) 1044 610 CONTINUE 1045 710 CONTINUE 1046 810 CONTINUE 1047 END IF 1048 END IF 1049C 1050 RETURN 1051 END 1052C /* Deck trsoao */ 1053 SUBROUTINE TRSOAO(FCSO,WORK,LWORK,IDIM,IATOM,IPRINT) 1054#include "implicit.h" 1055#include "priunit.h" 1056 DIMENSION FCSO(*), WORK(LWORK) 1057C 1058 IF (IPRINT .GT. 2) CALL TITLER('Subroutine TRSOAO','*',103) 1059C 1060C Allocation 1061C 1062 KFCAO = 1 1063 KTRMT = KFCAO + IDIM*(IDIM + 1)/2 1064 KWRK = KTRMT + IDIM*IDIM 1065 LMAX = KWRK + IDIM 1066 IF (LMAX.GT.LWORK) CALL STOPIT('TRSOAO',' ',LMAX,LWORK) 1067C 1068C Set up AO to SO transformation matrix 1069C 1070 CALL TRAORB(WORK(KTRMT),IDIM,IATOM,IPRINT) 1071C 1072C Transform one-electron matrix 1073C 1074 CALL DZERO(WORK(KFCAO),IDIM*(IDIM + 1)/2) 1075 CALL UTHU(FCSO,WORK(KFCAO),WORK(KTRMT), 1076 * WORK(KWRK),IDIM,IDIM) 1077 IF (IPRINT .GT. 6) THEN 1078 CALL HEADER('Matrix transformed from SO to AO basis',-1) 1079 CALL OUTPAK(WORK(KFCAO),IDIM,1,LUPRI) 1080 END IF 1081 RETURN 1082 END 1083C /* Deck traoso */ 1084 SUBROUTINE TRAOSO(FCAO,WORK,LWORK,IDIM,IATOM,IPRINT) 1085#include "implicit.h" 1086#include "priunit.h" 1087#include "maxorb.h" 1088 DIMENSION FCAO(*), WORK(LWORK) 1089#include "aosotr.h" 1090C 1091 IF (IPRINT .GT. 2) CALL TITLER('Subroutine TRAOSO','*',103) 1092C 1093C Allocation 1094C 1095 KFCSO = 1 1096 KTRMT = KFCSO + IDIM*(IDIM + 1)/2 1097 KWRK = KTRMT + IDIM*IDIM 1098 LMAX = KWRK + IDIM 1099 IF (LMAX.GT.LWORK) CALL STOPIT('TRAOSO',' ',LMAX,LWORK) 1100C 1101C Construct transformation matrix 1102C 1103 CALL TRAORB(WORK(KTRMT),IDIM,IATOM,IPRINT) 1104 CALL DGETRN(WORK(KTRMT),IDIM,IDIM) 1105C 1106C Transform one-electron matrix 1107C 1108 CALL DZERO(WORK(KFCSO),IDIM*(IDIM + 1)/2) 1109 CALL UTHU(FCAO,WORK(KFCSO),WORK(KTRMT), 1110 * WORK(KWRK),IDIM,IDIM) 1111 IF (IPRINT .GT. 6) THEN 1112 CALL HEADER('Matrix transformed from AO to SO basis',-1) 1113 CALL OUTPAK(WORK(KFCSO),IDIM,1,LUPRI) 1114 END IF 1115 RETURN 1116 END 1117C /* Deck traorb */ 1118 SUBROUTINE TRAORB(ORBTRA,IDIM,IATOM,IPRINT) 1119C 1120C Sets up transformation matrices between SO and AO orbitals 1121C tuh Sep 07 1988 1122C 1123C IATOM > 0 : SO -> AO transformation, modified for atom IATOM 1124C IATOM = 0 : SO -> AO transformation 1125C IATOM = -1 : AO -> SO transformation 1126C 1127#include "implicit.h" 1128#include "priunit.h" 1129#include "maxaqn.h" 1130#include "mxcent.h" 1131#include "maxorb.h" 1132 PARAMETER (D0 = 0.0D0) 1133#include "nuclei.h" 1134#include "shells.h" 1135#include "symmet.h" 1136#include "aosotr.h" 1137 DIMENSION ORBTRA(IDIM,IDIM) 1138 1139C 1140 IF (IATOM .GT. 0) THEN 1141 FACTOR = SQRT(FMULT(ISTBNU(IATOM))) 1142 ELSE IF (IATOM .EQ. 0 .OR. IATOM .EQ. -1) THEN 1143 FACTOR = 1.0D0 1144 ELSE 1145 CALL QENTER('TRAORB') 1146 CALL QUIT('IATOM .lt. -1') 1147 END IF 1148 CALL DZERO(ORBTRA,IDIM*IDIM) 1149 DO 100 I = 1, IDIM 1150 DO 200 J = 1, MAXREP + 1 1151 IF (IATOM .LE. 0) THEN 1152 K = ITRAN(I,J) 1153 ELSE 1154 K = IAOAO(ITRAN(I,J)) 1155 END IF 1156 IF (K .GT. 0) THEN 1157 ORBTRA(I,K) = CTRAN(I,J) 1158 END IF 1159 200 CONTINUE 1160 100 CONTINUE 1161 IF (IATOM .GE. 0) THEN 1162 DO 300 I = 1, IDIM 1163 SUM = D0 1164 DO 400 J = 1, IDIM 1165 SUM = SUM + ABS(ORBTRA(I,J)) 1166 400 CONTINUE 1167 SUM = 1.0D0/(FACTOR*SUM) 1168 DO 500 J = 1, IDIM 1169 ORBTRA(I,J) = ORBTRA(I,J)*SUM 1170 500 CONTINUE 1171 300 CONTINUE 1172 END IF 1173C 1174 IF (IPRINT .GE. 15) THEN 1175 CALL HEADER('SO to AO transformation matrix',-1) 1176 WRITE(LUPRI,*) ' for IATOM = ',IATOM 1177 CALL OUTPUT(ORBTRA,1,IDIM,1,IDIM,IDIM,IDIM,-1,LUPRI) 1178 END IF 1179 RETURN 1180 END 1181C /* Deck testfd */ 1182 SUBROUTINE TESTFD(IATOM,ICOOR,IPRINT,NOH1,NOH2,NODC,NODV,NOPV, 1183 * FTD,FCD,FVD,QD,DV) 1184C 1185C This subroutine calculates the total AO differentiated Fock 1186C matrix in MO basis and tests by taking the trace of the Fock 1187C matrix and comparing it with the appropriate integral 1188C contributions to the molecular gradient. 1189C 1190#include "implicit.h" 1191#include "priunit.h" 1192#include "mxcent.h" 1193 PARAMETER (D0 = 0.D00, TWO = 2.D00) 1194 INTEGER P, Q, T, U, TU, QU 1195 LOGICAL NOH1, NOH2, NODC, NODV, NOPV 1196 DIMENSION FTD(NOCCT,NORBT), FCD(NNORBT), FVD(NNORBT), 1197 * QD(NASHDI,NORBT), DV(*) 1198C 1199C Used from common blocks: 1200C INFDIM: NASHDI 1201C 1202#include "inforb.h" 1203#include "infdim.h" 1204#include "energy.h" 1205 ITRI(I) = I*(I - 1)/2 1206C 1207 CALL DZERO(FTD(1,1),NOCCT*NORBT) 1208C 1209C Inactive block 1210C 1211 IF (.NOT. NODC) THEN 1212 DO 200 I = 1, NISHT 1213 DO 210 Q = 1, NORBT 1214 IF (I .GE. Q) THEN 1215 IQ = ITRI(I) + Q 1216 ELSE 1217 IQ = ITRI(Q) + I 1218 END IF 1219 FTD(I,Q) = TWO*(FCD(IQ) + FVD(IQ)) 1220 210 CONTINUE 1221 200 CONTINUE 1222 END IF 1223C 1224C Active block 1225C 1226 DO 300 T = 1, NASHT 1227 DO 310 Q = 1, NORBT 1228 FTDTQ = D0 1229 IF (.NOT. NODV) THEN 1230 DO 400 U = 1, NASHT 1231 IU = NISHT + U 1232 IF (Q .GE. IU) THEN 1233 QU = ITRI(Q) + IU 1234 ELSE 1235 QU = ITRI(IU) + Q 1236 END IF 1237 IF (T .GE. U) THEN 1238 TU = ITRI(T) + U 1239 ELSE 1240 TU = ITRI(U) + T 1241 END IF 1242 FTDTQ = FTDTQ + DV(TU)*FCD(QU) 1243 400 CONTINUE 1244 END IF 1245 IF (.NOT. (NOH2 .OR. NOPV)) THEN 1246 FTDTQ = FTDTQ + QD(T,Q) 1247 END IF 1248 FTD(NISHT + T,Q) = FTDTQ 1249 310 CONTINUE 1250 300 CONTINUE 1251 IF (IPRINT .GE. 10) THEN 1252 CALL HEADER('AO derivative of total Fock matrix',-1) 1253 CALL OUTPUT(FTD(1,1),1,NOCCT,1,NORBT,NOCCT,NORBT,1,LUPRI) 1254 END IF 1255 TRACE = D0 1256 DO 500 P = 1, NOCCT 1257 TRACE = TRACE + FTD(P,P) 1258 500 CONTINUE 1259 IAC = 3*(IATOM - 1) + ICOOR 1260 GRDINT = D0 1261 IF (.NOT. NOH1) THEN 1262 GRDINT = GRDINT + GRADKE(IAC) + GRADNA(IAC) 1263 END IF 1264 IF (.NOT. NOH2) THEN 1265 GRDINT = GRDINT + TWO*GRADEE(IAC) 1266 END IF 1267 CALL HEADER('FOCK MATRIX TEST:',0) 1268 WRITE (LUPRI,'(8X,A,/)') 1269 * ' GRD(KE+NA+2*EE) TRACE(FTD) DIFFERENCE' 1270 WRITE (LUPRI,'(5X,3F16.10)') GRDINT, TRACE, GRDINT - TRACE 1271 RETURN 1272 END 1273C /* Deck testtr */ 1274 SUBROUTINE TESTTR(IATOM,ICOOR,IOPSYM,PV,H2D,QD,IPRINT) 1275#include "implicit.h" 1276#include "priunit.h" 1277C 1278C ***************************************************** 1279C *** Take trace of Q matrices and PV.H2D product *** 1280C *** as a check of transformation. This version *** 1281C *** handles symmetry. *** 1282C ***************************************************** 1283C 1284#include "maxaqn.h" 1285#include "maxorb.h" 1286#include "mxcent.h" 1287#include "iratdef.h" 1288 PARAMETER (D0 = 0.D00, D1 = 1.D00, TWO = 2.D00) 1289 PARAMETER (HALF = 0.5D00) 1290 INTEGER T, U, V, X, TU, VX 1291 DIMENSION PV(NNASHX,NNASHX), H2D(NNASHX,NNASHX) 1292 DIMENSION QD(NORBT,NASHT) 1293#include "symmet.h" 1294#include "inforb.h" 1295#include "energy.h" 1296 ITRI(I) = I*(I - 1)/2 1297C 1298 WRITE (LUPRI,'(/,A,/)') 1299 * ' ---------- Output from TESTTR ---------- ' 1300 IAC = IPTCNT(3*(IATOM-1)+ICOOR,IOPSYM-1,1) 1301 GRDINT = GRADEE(IAC) 1302C 1303C ***** H2D elements ***** 1304C 1305 IF (IPRINT .GT. 60) THEN 1306 CALL HEADER('Two-electron density elements',0) 1307 CALL OUTPUT(PV,1,NNASHX,1,NNASHX,NNASHX,NNASHX,1,LUPRI) 1308 CALL HEADER('Derivative MO two-electron integrals',0) 1309 CALL OUTPUT(H2D(1,1),1,NNASHX,1,NNASHX, 1310 * NNASHX,NNASHX,1,LUPRI) 1311 END IF 1312 TRACE = D0 1313C 1314C Extra factor of half for new definition of 2-matrix 1315C 1316 DO 100 T = 1, NASHT 1317 DO 110 U = 1, T 1318 TU = ITRI(T) + U 1319 FACTU = HALF 1320 IF (T .NE. U) FACTU = TWO*FACTU 1321 DO 120 V = 1, NASHT 1322 DO 130 X = 1, V 1323 VX = ITRI(V) + X 1324 FACVX = D1 1325 IF (V .NE. X) FACVX = TWO 1326 TRACE = TRACE + FACTU*FACVX*PV(TU,VX)*H2D(TU,VX) 1327130 CONTINUE 1328120 CONTINUE 1329110 CONTINUE 1330100 CONTINUE 1331 CALL HEADER('Test on H2D matrix:',0) 1332 WRITE (LUPRI,'(9X,A,/)') 1333 * ' GRADEE(PV) PV*H2D Difference ' 1334 WRITE (LUPRI,'(5X,3F15.10)') GRDINT, TRACE, GRDINT - TRACE 1335 WRITE (LUPRI,'(/,A,/)') 1336 * ' Note: Only PV must be included in TWOEXP!' 1337C 1338C ***** QD elements ***** 1339C 1340 IF (IPRINT .GT. 60) THEN 1341 CALL HEADER('QD matrix',-1) 1342 CALL OUTPUT(QD(1,1),1,NORBT,1,NASHT,NORBT,NASHT,1,LUPRI) 1343 END IF 1344C 1345C Loop over symmetries here and accumulate only active-active terms 1346C 1347 TRACE = D0 1348 IOFF = 0 1349 JOFF = 0 1350 DO 300 ISYM = 1,NSYM 1351 NASHI = NASH(ISYM) 1352 NISHI = NISH(ISYM) 1353 NORBI = NORB(ISYM) 1354 IF (NASHI .NE. 0) THEN 1355 DO 310 I = 1, NASHI 1356 TRACE = TRACE + QD(IOFF+NISHI+I,JOFF+I) 1357310 CONTINUE 1358 ENDIF 1359 IOFF = IOFF + NORBI 1360 JOFF = JOFF + NASHI 1361300 CONTINUE 1362 GRDINT = TWO*GRDINT 1363 CALL HEADER('Test on QD matrix:',0) 1364 WRITE (LUPRI,'(9X,A,/)') 1365 * ' 2*GRADEE(PV) Trace(QD) Difference ' 1366 WRITE (LUPRI,'(5X,3F15.10)') GRDINT, TRACE, GRDINT - TRACE 1367 WRITE (LUPRI,'(/,A,/)') 1368 * ' Note: Only PV must be included in TWOEXP!' 1369 RETURN 1370 END 1371C /* Deck testsd */ 1372 SUBROUTINE TESTSD(TD,FT,IATOM,DERIV,IPRINT) 1373C 1374C This routine performs a test on differentiated overlap matrices by 1375C calculating the gradient reorthonormalization in MO basis and 1376C comparing it with the results obtained in AO basis. 1377C 1378#include "implicit.h" 1379#include "priunit.h" 1380#include "maxaqn.h" 1381#include "mxcent.h" 1382#include "maxorb.h" 1383#include "iratdef.h" 1384 PARAMETER (D0 = 0.D00, TWO = 2.D00) 1385 LOGICAL DERIV(3) 1386 DIMENSION TD(NNORBX,3), FT(*), TRACE(3) 1387#include "symmet.h" 1388#include "inforb.h" 1389#include "energy.h" 1390C 1391 CALL HEADER('Reorthonormalization of gradient from'// 1392 * ' overlap matrices in MO basis',0) 1393 IF (IPRINT .GT. 5) THEN 1394 DO 50 ISYM = 1, NSYM 1395 WRITE (LUPRI,'(A,I5)') ' Irrep', ISYM 1396 NORBI = NORB(ISYM) 1397 ISTRI = I2ORB(ISYM) + 1 1398 CALL HEADER('Total Fock matrix',-1) 1399 CALL OUTPUT(FT(ISTRI),1,NORBI,1,NORBI,NORBI,NORBI,1,LUPRI) 1400 50 CONTINUE 1401 END IF 1402 DO 100 ICOOR = 1, 3 1403 IF (DERIV(ICOOR)) THEN 1404 IF (IPRINT .GT. 5) THEN 1405 CALL HEADER('TD matrix',-1) 1406 CALL OUTPAK(TD(1,ICOOR),NORBT,1,LUPRI) 1407 END IF 1408 TRAC = D0 1409 DO 150 ISYM = 1, NSYM 1410 ISTR = IORB(ISYM) + 1 1411 IEND = IORB(ISYM) + NORB(ISYM) 1412 IJ = I2ORB(ISYM) 1413 DO 200 I = ISTR, IEND 1414 DO 300 J = ISTR, IEND 1415 IJ = IJ + 1 1416 IJTD = (MAX(I,J)*(MAX(I,J) - 3))/2 + I + J 1417 TRAC = TRAC + FT(IJ)*TD(IJTD,ICOOR) 1418 300 CONTINUE 1419 200 CONTINUE 1420 150 CONTINUE 1421 TRACE(ICOOR) = TWO*TRAC 1422 END IF 1423 100 CONTINUE 1424 IOFF = 3*(IATOM - 1) 1425 CALL HEADER 1426 * ('Coordinate AO basis MO basis Difference',10) 1427 DO 500 ICOOR = 1, 3 1428 IF (DERIV(ICOOR)) THEN 1429 JCOOR = IPTCNT(IOFF + ICOOR,0,1) 1430 IF (JCOOR .GT. 0) THEN 1431 FS = GRADFS(JCOOR) 1432 TRAC = TRACE(ICOOR) 1433 WRITE (LUPRI,'(5X,I10,5X,3F15.10)') 1434 * JCOOR, FS, TRAC, FS - TRAC 1435 END IF 1436 END IF 1437 500 CONTINUE 1438 RETURN 1439 END 1440C /* Deck onetra */ 1441 SUBROUTINE ONETRA(CMO,SQMO,SQSO,WORK,LWORK,IATOM,ICOOR,IPRINT) 1442C 1443C This routine transforms derivative SO matrices to MO basis. 1444C 1445C tuh 010988 1446C 1447#include "implicit.h" 1448#include "priunit.h" 1449#include "maxaqn.h" 1450#include "mxcent.h" 1451#include "maxorb.h" 1452#include "symmet.h" 1453#include "inforb.h" 1454#include "infdim.h" 1455 DIMENSION CMO(NCMOT), WORK(LWORK), 1456 * SQMO(NORBT,NORBT), SQSO(N2BASX) 1457 LOGICAL NONSYM 1458 PARAMETER (D1 =1.0D0, D2 = 2.0D0, DP5 = 0.5D0) 1459C 1460C ***** Print Section ***** 1461C 1462 IF (IPRINT .GT. 5) THEN 1463 CALL HEADER('Output from ONETRA',-1) 1464 WRITE (LUPRI,'(A,I5)') ' IATOM ', IATOM 1465 WRITE (LUPRI,'(A,I5)') ' ICOOR ', ICOOR 1466 WRITE (LUPRI,'(A,I5)') ' NORBT ', NORBT 1467 WRITE (LUPRI,'(A,I5)') ' NBAST ', NBAST 1468 WRITE (LUPRI,'(A,2I5)')' NNBASX, N2BASX ', NNBASX, N2BASX 1469 WRITE (LUPRI,'(A,2I5)')' NNORBX, N2ORBX ', NNORBX, N2ORBX 1470 WRITE (LUPRI,'(A,I5)') ' NCMOT ', NCMOT 1471 WRITE (LUPRI,'(A,I10)')' LWORK ', LWORK 1472 END IF 1473C 1474C ***** Make square matrix ***** 1475C 1476 IF (IPRINT .GT. 15) THEN 1477 CALL HEADER('Square SO matrix',-1) 1478 CALL OUTPUT(SQSO,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1479 END IF 1480C 1481C ***** Transform matrix from SO to MO basis ***** 1482C 1483 CALL DZERO(SQMO,N2ORBX) 1484C 1485C Loop over irreps for orbitals 1486C 1487 NONSYM = .FALSE. 1488 JCOOR = MAX(1,3*(IATOM - 1) + ICOOR) 1489C ... hjaaj: to avoid referring to IPTCNT(-2,...) below when IATOM.eq.0 1490 DO 100 IREPA = 1,NSYM 1491 IF (NORB(IREPA) .GT. 0) THEN 1492 DO 200 IREPB = 1, IREPA 1493 IF (NORB(IREPB) .GT. 0) THEN 1494C 1495C Check if operator belongs to this symmetry 1496C 1497 IREPO = MULD2H(IREPA,IREPB) 1498 IF ((IATOM.EQ.0 .AND. (IREPO-1).EQ.ISYMAX(ICOOR,2)) .OR. 1499 & (IATOM.GT.0 .AND. IPTCNT(JCOOR,IREPO-1,1).GT.0)) THEN 1500 NONSYM = NONSYM .OR. IREPO .GT. 1 1501C 1502C Print symmetries and orbitals 1503C 1504 IF (IPRINT.GT.15) THEN 1505 WRITE (LUPRI,'(A,3I3)') ' IREPA/B/O',IREPA,IREPB,IREPO 1506 WRITE (LUPRI,'(/A,I5,A)') 1507 * ' MO coefficients for symmetry', IREPA 1508 CALL OUTPUT(CMO(ICMO(IREPA)+1),1,NBAS(IREPA),1, 1509 * NORB(IREPA),NBAS(IREPA),NORB(IREPA),1,LUPRI) 1510 WRITE (LUPRI,'(/A,I5,A)') 1511 * ' MO coefficients for symmetry', IREPB 1512 CALL OUTPUT(CMO(ICMO(IREPB)+1),1,NBAS(IREPB),1, 1513 * NORB(IREPB),NBAS(IREPB),NORB(IREPB),1,LUPRI) 1514 END IF 1515C 1516C Transform matrix block(s) 1517C 1518 CALL UTHV(CMO(ICMO(IREPA)+1),SQSO,CMO(ICMO(IREPB)+1), 1519 * IREPA,IREPB,NBAS(IREPA),NBAS(IREPB),SQMO,WORK) 1520 IF (NONSYM) THEN 1521 CALL UTHV(CMO(ICMO(IREPB)+1),SQSO,CMO(ICMO(IREPA)+1), 1522 & IREPB,IREPA,NBAS(IREPB),NBAS(IREPA),SQMO,WORK) 1523 END IF 1524C 1525C Print transformed matrix thus far 1526C 1527 IF (IPRINT.GT.25) THEN 1528 WRITE(LUPRI,'(/4A)')' Unfinished matrix in MO basis' 1529 CALL OUTPUT(SQMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1530 END IF 1531 END IF 1532 END IF 1533 200 CONTINUE 1534 END IF 1535 100 CONTINUE 1536 IF (IPRINT.GT.15) THEN 1537 CALL HEADER('Matrix in MO basis',-1) 1538 CALL OUTPUT(SQMO,1,NORBT,1,NORBT,NORBT,NORBT,1,LUPRI) 1539 END IF 1540 RETURN 1541 END 1542C /* Deck transx */ 1543 SUBROUTINE TRANSX(A,B,NA,NB,ANTSYM,IPRINT) 1544C 1545C Modified from TRANSA, tuh 2-Sep-1988 1546C TRANSA written 1-Feb-1984 PJ 1547C 1548C Put the lower triangle of A into the upper triangle of B 1549C (if A = B, the matrix will become (anti)symmetric 1550C based on the lower triangle) 1551C 1552#include "implicit.h" 1553#include "priunit.h" 1554C 1555C 1556 DIMENSION A(NA,*), B(NB,*) 1557 IF (IPRINT.GT.15) THEN 1558 WRITE(LUPRI,150) NA,NB,ANTSYM 1559 150 FORMAT(/' Dimension of matrix A is ',I5,I5, 1560 & /' ANTSYM =',F10.2) 1561 WRITE(LUPRI,'(/A)')' Matrix to be (anti)symmetrized ' 1562 CALL OUTPUT(A(1,1),1,NA,1,NB,NA,NB,1,LUPRI) 1563 END IF 1564 DO 200 I=1,NB 1565 DO 100 J=1,I 1566 B(J,I)=A(I,J)*ANTSYM 1567 100 CONTINUE 1568 200 CONTINUE 1569 IF (IPRINT.GT.15) THEN 1570 WRITE(LUPRI,'(/A)')' Matrix after (anti)symmetrization ' 1571 CALL OUTPUT(B(1,1),1,NB,1,NA,NB,NA,1,LUPRI) 1572 END IF 1573 RETURN 1574 END 1575C /* Deck getsd */ 1576 SUBROUTINE GETSD(DERIV,CMO,TD,WORK,LWORK,IATOM,DIAGTD, 1577 & KEY,DOTD,IPRINT,INTPRI) 1578C 1579C This subroutine calculates the AO-differentiated overlap matrices 1580C in MO basis. 1581C 1582C 1583C Due to the need of square matrices in constructing right-hand 1584C sides when using differentiated creation operators, all matrices 1585C in this routine are square, all packing is to be done outside this 1586C routine, K.Ruud, Nov.-93 1587C 1588C tuh 051088 1589C 1590#include "implicit.h" 1591#include "priunit.h" 1592#include "mxcent.h" 1593 LOGICAL DERIV(3), DIAGTD, DIVIDE, DOTD 1594 PARAMETER (D1 = 1.0D0) 1595 CHARACTER*4 KEY 1596#include "abainf.h" 1597#include "inforb.h" 1598 DIMENSION CMO(*), TD(N2BASX,3), WORK(LWORK) 1599 CALL QENTER('GETSD') 1600 IF (IPRINT .GT. 05) THEN 1601 CALL TITLER('Output from GETSD','*',103) 1602 WRITE (LUPRI, '(A,3L5)') ' DERIV ', (DERIV(I),I=1,3) 1603 WRITE (LUPRI, '(A,I10)') ' LWORK ', LWORK 1604 WRITE (LUPRI, '(A,I5)') ' IATOM ', IATOM 1605 WRITE (LUPRI, '(A,L5)') ' DIAGTD', DIAGTD 1606 END IF 1607 CALL DZERO(TD,3*N2BASX) 1608 CALL ONEDRL(KEY,TD(1,1),TD(1,2),TD(1,3),IATOM, 1609 & DERIV(1),DERIV(2),DERIV(3),WORK,LWORK,IPRINT,INTPRI) 1610 IF (IPRINT .GT. 10) THEN 1611 CALL HEADER('SO overlap matrix - x direction',-1) 1612 CALL OUTPUT(TD(1,1),1,NBAST,1,NBAST,NBAST,NBAST, 1613 & 1,LUPRI) 1614 CALL HEADER('SO overlap matrix - y direction',-1) 1615 CALL OUTPUT(TD(1,2),1,NBAST,1,NBAST,NBAST,NBAST, 1616 & 1,LUPRI) 1617 CALL HEADER('SO overlap matrix - z direction',-1) 1618 CALL OUTPUT(TD(1,3),1,NBAST,1,NBAST,NBAST,NBAST, 1619 & 1,LUPRI) 1620 END IF 1621 KSMOSQ = 1 1622 KGETSD = KSMOSQ + N2BASX 1623 LGETSD = LWORK - KGETSD + 1 1624 IF (KGETSD.GT.LWORK) CALL STOPIT('GETSD',' ',KGETSD,LWORK) 1625 DO 100 ICOOR = 1, 3 1626 IF (DERIV(ICOOR)) THEN 1627 IF (IPRINT .GE. 10) THEN 1628 IF (ICOOR. EQ. 1) THEN 1629 CALL TITLER('Differentiation in x direction',' ',200) 1630 ELSE IF (ICOOR. EQ. 2) THEN 1631 CALL TITLER('Differentiation in y direction',' ',200) 1632 ELSE 1633 CALL TITLER('Differentiation in z direction',' ',200) 1634 END IF 1635 END IF 1636 CALL ONETRA(CMO,WORK(KSMOSQ),TD(1,ICOOR),WORK(KGETSD), 1637 & LGETSD,IATOM,ICOOR,IPRINT) 1638 CALL DCOPY(N2ORBX,WORK(KSMOSQ),1,TD(1,ICOOR),1) 1639 IF (DOTD) THEN 1640 IF (KEY .EQ. 'SMAT' .OR. KEY .EQ. 'OMAT') THEN 1641 DIVIDE = .TRUE. 1642 ELSE 1643 DIVIDE = .FALSE. 1644 END IF 1645 CALL GETTD(TD(1,ICOOR),NORBT,DIAGTD,DIVIDE,IPRINT) 1646 IF (IPRINT .GE. 10) THEN 1647 CALL HEADER('TD matrix',-1) 1648 CALL OUTPUT(TD(1,ICOOR),1,NORBT,1,NORBT,NORBT,NORBT, 1649 & 1,LUPRI) 1650 END IF 1651 ELSE 1652 IF (IPRINT .GE. 10) THEN 1653 CALL HEADER('SO overlap matrix',-1) 1654 CALL OUTPUT(TD(1,ICOOR),1,NORBT,1,NORBT,NORBT,NORBT, 1655 & 1,LUPRI) 1656 END IF 1657 END IF 1658 END IF 1659 100 CONTINUE 1660 CALL QEXIT('GETSD') 1661 RETURN 1662 END 1663C /* Deck onedrl */ 1664 SUBROUTINE ONEDRL(KEY,OMATX,OMATY,OMATZ,IATOM,DERX,DERY,DERZ,WORK, 1665 & LWORK,IPRINT,INTPRI) 1666C 1667 use pelib_interface, only: use_pelib, pelib_ifc_london 1668#include "implicit.h" 1669#include "priunit.h" 1670#include "dummy.h" 1671#include "mxcent.h" 1672#include "maxorb.h" 1673 PARAMETER (D0 = 0.0D0, DP5 = 0.5D0, D1 = 1.0D0) 1674 INTEGER TMAT 1675 LOGICAL DERX, DERY, DERZ 1676 CHARACTER*4 KEY 1677 CHARACTER*7 KEY2 1678#include "inforb.h" 1679#include "gnrinf.h" 1680#include "abainf.h" 1681#include "cbisol.h" 1682#include "infinp.h" 1683#include "pcmlog.h" 1684#include "qmmm.h" 1685#include "efield.h" 1686#include "infpar.h" 1687 DIMENSION OMATX(N2BASX), OMATY(N2BASX), OMATZ(N2BASX), 1688 & WORK(LWORK) 1689 IF (IPRINT .GT. 5) CALL TITLER('Output from ONEDRL','*',103) 1690C 1691C Field derivatives 1692C ================= 1693C 1694 IF (IATOM .EQ. 0) THEN 1695 KMATS = 1 1696 KWRK = KMATS + 3*N2BASX 1697 LWRK = LWORK - KWRK + 1 1698 IF (KWRK .GT. LWORK) CALL STOPIT('ONEDRL','GET1PR',KWRK,LWORK) 1699 KEY2 = 'S1MAG ' 1700 IF (KEY .EQ. 'BMAT') KEY2 = 'S1MAGR ' 1701 IF ((KEY .EQ. 'HMAT') .AND. (.NOT. NOLOND)) KEY2 = 'MAGMOM ' 1702 IF ((KEY .EQ. 'HMAT') .AND. NOLOND) KEY2 = 'ANGMOM ' 1703 NCOMP = 3 1704 IF (KEY2 .EQ. 'S1MAGR ') THEN 1705 CALL GET1PR(WORK(KMATS),KEY2,NCOMP,'SQUARE',.FALSE., 1706 & WORK(KWRK),LWRK,IDUMMY,INTPRI) 1707 ELSE 1708 CALL GET1PR(WORK(KMATS),KEY2,NCOMP,'SQUARE',.TRUE., 1709 & WORK(KWRK),LWRK,IDUMMY,INTPRI) 1710 END IF 1711 IF (KEY2.EQ. 'ANGMOM') CALL DSCAL(3*N2BASX,DP5,WORK(KMATS),1) 1712C 1713C Add contribution from electric field 1714C 1715 IF (NFIELD .GT. 0 .AND. (KEY .EQ. 'HMAT') 1716 & .AND. (.NOT. NOLOND)) THEN 1717 KTMAT = KWRK 1718 MWRK = KTMAT + 3*N2BASX 1719 IF (MWRK .GT. LWRK) CALL STOPIT('ONEDRL','CM1 ',MWRK, 1720 & LWRK) 1721 LLEFT = LWRK - MWRK + 1 1722 DO 556 IFIELD = 1, NFIELD 1723 IF (EFIELD(IFIELD) .NE. D0) THEN 1724 IF (LFIELD(IFIELD)(2:7) .EQ. 'DIPLEN') THEN 1725 IF (LFIELD(IFIELD) .EQ. 'XDIPLEN ') THEN 1726 FIELD1 = 'X-FIELD' 1727 ELSE IF (LFIELD(IFIELD) .EQ. 'YDIPLEN ') THEN 1728 FIELD1 = 'Y-FIELD' 1729 ELSE IF (LFIELD(IFIELD) .EQ. 'ZDIPLEN ') THEN 1730 FIELD1 = 'Z-FIELD' 1731 ELSE 1732 WRITE (LUPRI,'(/,3A,/)') 'Field type ', 1733 & LFIELD(IFIELD), 1734 & ' not implemented for magnetic properties' 1735 CALL QUIT('Illegal field type for '// 1736 & 'magnetic properties') 1737 END IF 1738 NCOMP = 3 1739 CALL GET1PR(WORK(KTMAT),'CM1 ',NCOMP,'SQUARE', 1740 & .TRUE.,WORK(MWRK),LLEFT,IDUMMY,INTPRI) 1741 ELSE IF (LFIELD(IFIELD)(3:7) .EQ. 'THETA' .OR. 1742 & LFIELD(IFIELD)(3:8) .EQ. 'SECMOM') THEN 1743 IF (LFIELD(IFIELD) .EQ. 'XXTHETA ') THEN 1744 FIELD3 = 'XX-FGRD' 1745 ELSE IF (LFIELD(IFIELD) .EQ. 'YYTHETA ') THEN 1746 FIELD3 = 'YY-FGRD' 1747 ELSE IF (LFIELD(IFIELD) .EQ. 'ZZTHETA ') THEN 1748 FIELD3 = 'ZZ-FGRD' 1749 ELSE IF (LFIELD(IFIELD) .EQ. 'XYTHETA ') THEN 1750 FIELD3 = 'XY-FGRD' 1751 ELSE IF (LFIELD(IFIELD) .EQ. 'XZTHETA ') THEN 1752 FIELD3 = 'XZ-FGRD' 1753 ELSE IF (LFIELD(IFIELD) .EQ. 'YZTHETA ') THEN 1754 FIELD3 = 'YZ-FGRD' 1755 ELSE IF (LFIELD(IFIELD) .EQ. 'XXSECMOM') THEN 1756 FIELD3 = 'XX-2MOM' 1757 ELSE IF (LFIELD(IFIELD) .EQ. 'YYSECMOM') THEN 1758 FIELD3 = 'YY-2MOM' 1759 ELSE IF (LFIELD(IFIELD) .EQ. 'ZZSECMOM') THEN 1760 FIELD3 = 'ZZ-2MOM' 1761 ELSE IF (LFIELD(IFIELD) .EQ. 'XYSECMOM') THEN 1762 FIELD3 = 'XY-2MOM' 1763 ELSE IF (LFIELD(IFIELD) .EQ. 'XZSECMOM') THEN 1764 FIELD3 = 'XZ-2MOM' 1765 ELSE IF (LFIELD(IFIELD) .EQ. 'YZSECMOM') THEN 1766 FIELD3 = 'YZ-2MOM' 1767 ELSE 1768 WRITE (LUPRI,'(/,3A,/)') 'Field type ', 1769 & LFIELD(IFIELD), 1770 & ' not implemented for magnetic properties' 1771 CALL QUIT('Illegal field type for '// 1772 & 'magnetic properties') 1773 END IF 1774 NCOMP = 3 1775 CALL GET1PR(WORK(KTMAT),'QDBINT ',NCOMP,'SQUARE', 1776 & .TRUE.,WORK(MWRK),LLEFT,IDUMMY,INTPRI) 1777 END IF 1778 CALL DAXPY(N2BASX*3,EFIELD(IFIELD),WORK(KTMAT),1, 1779 & WORK(KMATS),1) 1780 END IF 1781 556 CONTINUE 1782 END IF 1783C 1784 IF ((SOLVNT .OR. PCM) .AND. (KEY .EQ. 'HMAT') 1785 & .AND. (.NOT. NOLOND)) THEN 1786 KTMAT = KWRK 1787 MWRK = KTMAT + 3*N2BASX 1788 IF (MWRK .GT. LWRK) CALL STOPIT('ONEDRL','MAGSOL',MWRK, 1789 & LWRK) 1790 LLEFT = LWRK - MWRK + 1 1791 IF (SOLVNT) THEN 1792 CALL MAGSOL(WORK(KTMAT),3,IPRINT,INTPRI, 1793 & WORK(MWRK),LLEFT) 1794 ELSE 1795 CALL MAGPCMFIRST(WORK(KTMAT),3,IPRINT,INTPRI,WORK(MWRK), 1796 & LLEFT) 1797 END IF 1798 CALL DAXPY(N2BASX*3,D1,WORK(KTMAT),1,WORK(KMATS),1) 1799 END IF 1800C 1801 IF (QM3 .AND. (KEY .EQ. 'HMAT') .AND. .NOT. NOLOND) THEN 1802 KTMAT = KWRK 1803 MWRK = KTMAT + 3*N2BASX 1804 IF (MWRK .GT. LWRK) CALL STOPIT('ONEDRL','QM3-1',MWRK,LWRK) 1805 LLEFT = LWRK - MWRK + 1 1806C Calls QM3FIRST_P if parallel and with #nodes>1, Arnfinn nov. -08 1807 IF (NODTOT.GE.1) THEN 1808 CALL QM3FIRST_P(WORK(KTMAT),3,IPRINT,INTPRI,WORK(MWRK), 1809 & LLEFT) 1810 ELSE 1811 CALL QM3FIRST(WORK(KTMAT),3,IPRINT,INTPRI,WORK(MWRK), 1812 & LLEFT) 1813 ENDIF 1814 CALL DAXPY(N2BASX*3,D1,WORK(KTMAT),1,WORK(KMATS),1) 1815 END IF 1816 1817 IF (QMMM .AND. (KEY .EQ. 'HMAT') .AND. .NOT. NOLOND) THEN 1818 IF (IPQMMM .GT. 1) CALL TIMER('START ',TIMSTR,TIMEND) 1819 KTMAT = KWRK 1820 MWRK = KTMAT + 3*N2BASX 1821 IF (MWRK .GT. LWRK) CALL STOPIT('ONEDRL','QMMM-1',MWRK,LWRK) 1822 LLEFT = LWRK - MWRK + 1 1823 CALL QMMMFIRST(WORK(KTMAT),IPRINT,INTPRI,WORK(MWRK),LLEFT) 1824 CALL DAXPY(N2BASX*3,D1,WORK(KTMAT),1,WORK(KMATS),1) 1825 IF (IPQMMM .GT. 1) CALL TIMER('QMMMFIRST',TIMSTR,TIMEND) 1826 END IF 1827 1828 IF (USE_PELIB() .AND. (KEY.EQ.'HMAT') .AND. .NOT. NOLOND) THEN 1829 KTMAT = KWRK 1830 MWRK = KTMAT + 3 * N2BASX 1831 IF (MWRK > LWRK) CALL STOPIT('ONEDRL', 'PELIB', MWRK, LWRK) 1832 CALL PELIB_IFC_LONDON(WORK(KTMAT)) 1833 CALL DAXPY(3*N2BASX, D1, WORK(KTMAT), 1, WORK(KMATS), 1) 1834 END IF 1835 1836 IF (DERX) CALL DCOPY(N2BASX,WORK(KMATS ),1,OMATX,1) 1837 IF (DERY) CALL DCOPY(N2BASX,WORK(KMATS+ N2BASX),1,OMATY,1) 1838 IF (DERZ) CALL DCOPY(N2BASX,WORK(KMATS+2*N2BASX),1,OMATZ,1) 1839C 1840C Geometrical derivatives 1841C ======================= 1842C 1843 ELSE 1844 IF (LWORK.LT.(NNBASX + N2BASX)) 1845 & CALL STOPIT('ONEDRL','ONEDER',LWORK,N2BASX) 1846 IF (DERX) THEN 1847 CALL ONEDER(KEY,OMATX,WORK(1),WORK(1+NNBASX), 1848 & IATOM,1,IPRINT) 1849 END IF 1850 IF (DERY) THEN 1851 CALL ONEDER(KEY,OMATY,WORK(1),WORK(1+NNBASX), 1852 & IATOM,2,IPRINT) 1853 END IF 1854 IF (DERZ) THEN 1855 CALL ONEDER(KEY,OMATZ,WORK(1),WORK(1+NNBASX), 1856 & IATOM,3,IPRINT) 1857 END IF 1858 END IF 1859C 1860C Print 1861C ===== 1862C 1863 IF (IPRINT .GE. 15) THEN 1864 WRITE (LUPRI,'(2X,2A,I3)') 'Integral type and atom:',KEY,IATOM 1865 IF (DERX) THEN 1866 CALL HEADER('X component in ONEDRL',-1) 1867 CALL OUTPUT(OMATX,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1868 END IF 1869 IF (DERY) THEN 1870 CALL HEADER('Y component in ONEDRL',-1) 1871 CALL OUTPUT(OMATY,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1872 END IF 1873 IF (DERZ) THEN 1874 CALL HEADER('Z component in ONEDRL',-1) 1875 CALL OUTPUT(OMATZ,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 1876 END IF 1877 END IF 1878 RETURN 1879 END 1880C /* Deck oneder */ 1881 SUBROUTINE ONEDER(KEY,ONEMAT,TRONEM,WORK,ICENT,ICOOR,IPRINT) 1882C 1883C 130988 tuh 1884C 1885C This subroutine retrieves one-electron differentiated matrices 1886C (one-electron Hamiltonian or overlap) 1887C 1888#include "implicit.h" 1889#include "priunit.h" 1890#include "maxaqn.h" 1891#include "maxorb.h" 1892#include "mxcent.h" 1893C 1894 LOGICAL DNULL, TOTSYM, FNDLAB 1895 CHARACTER*4 KEY 1896 CHARACTER*8 LABEL 1897#include "nuclei.h" 1898#include "symmet.h" 1899#include "inforb.h" 1900#include "inftap.h" 1901#include "infinp.h" 1902#include "chrnos.h" 1903 DIMENSION WORK(N2BASX), ONEMAT(N2BASX), TRONEM(NNBASX) 1904 1905C 1906 IF (IPRINT .GE. 10) CALL TITLER('Output from ONEDER','*',103) 1907 CALL DZERO(TRONEM,NNBASX) 1908 CALL DZERO(WORK,N2BASX) 1909 DO 100 IREPO = 0, MAXREP 1910 TOTSYM = IREPO .EQ. 0 1911 ISCOOR = IPTCNT(3*(ICENT-1)+ICOOR,IREPO,1) 1912 IF (ISCOOR .GT. 0) THEN 1913 IF (KEY .EQ. 'OMAT' .OR. KEY .EQ. 'HMAT') THEN 1914 DO 200 IREPA = 0, MAXREP 1915 IREPB = IEOR(IREPO,IREPA) 1916 IF (IREPA .GE. IREPB) THEN 1917 NDIMA = NBAS(IREPA + 1) 1918 NDIMB = NBAS(IREPB + 1) 1919 NDIM = NDIMA*NDIMB 1920 IF (NDIM .GT. 0) THEN 1921 IF (TOTSYM) NDIM = NDIMA*(NDIMA + 1)/2 1922C 1923C ***** Get elements ***** 1924C 1925 CALL GETONE(KEY,ISCOOR,IREPA,NDIM,WORK,DNULL) 1926C 1927C ***** Print section ***** 1928C 1929 IF (IPRINT .GE. 15) THEN 1930 WRITE (LUPRI,'(//A,2X,A,I5)') 1931 * ' Coordinate and symmetry: ', 1932 * NAMEX(3*ICENT - 3 + ICOOR), IREPO 1933 WRITE (LUPRI,'(/A,2I5)') 1934 * ' Symmetry of orbitals:', IREPA, IREPB 1935 WRITE (LUPRI,'(A,2I5)') 1936 * ' Number of orbitals: ', NDIMA, NDIMB 1937 IF (.NOT.DNULL) THEN 1938 IF (TOTSYM) THEN 1939 CALL OUTPAK(WORK,NDIMA,1,LUPRI) 1940 ELSE 1941 CALL OUTPUT(WORK,1,NDIMB,1,NDIMA, 1942 * NDIMB,NDIMA,1,LUPRI) 1943 END IF 1944 END IF 1945 END IF 1946C 1947C ***** Add elements to matrix ***** 1948C 1949 IF (.NOT.DNULL) THEN 1950 IJ = 0 1951 IOFFA = IBAS(IREPA + 1) 1952 IOFFB = IBAS(IREPB + 1) 1953 DO 300 I = IOFFA + 1, IOFFA + NDIMA 1954 DO 400 J = IOFFB + 1, MIN(IOFFB + NDIMB,I) 1955 IJ = IJ + 1 1956 TRONEM((I*(I - 1))/2 + J) = WORK(IJ) 1957 400 CONTINUE 1958 300 CONTINUE 1959 END IF 1960 END IF 1961 END IF 1962 200 CONTINUE 1963 IF (NFIELD .GT. 0 .AND. KEY .EQ.'HMAT') THEN 1964 CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ', 1965 & 'UNFORMATTED',IDUMMY,.FALSE.) 1966 DO 556 IFIELD = 1, NFIELD 1967 IF (EFIELD(IFIELD) .NE. 0.D0) THEN 1968 IF (LFIELD(IFIELD)(2:8) .NE. 'DIPLEN ') THEN 1969 WRITE (LUPRI,'(/,3A,/)') 'Field type ', 1970 & LFIELD(IFIELD), 1971 & ' not implemented for geometric derivatives' 1972 CALL QUIT('Illegal field type for '// 1973 & 'geometric derivative') 1974 END IF 1975C 1976C Determine label to search for 1977C 1978 LABEL = CHRNOS(ISCOOR/100) 1979 & //CHRNOS(MOD(ISCOOR,100)/10) 1980 & //CHRNOS(MOD((MOD(ISCOOR,100)),10)) 1981 & //'DPG '// LFIELD(IFIELD)(1:1) 1982 REWIND LUPROP 1983 IF (.NOT. FNDLAB(LABEL,LUPROP)) THEN 1984 WRITE(LUPRI,'(1X,3A)')'Label ', 1985 & LABEL, ' not found on file AOPROPER' 1986 CALL QUIT('Error in ONEDER') 1987 END IF 1988 CALL READT(LUPROP,NNBASX,WORK) 1989 CALL DAXPY(NNBASX,-EFIELD(IFIELD), 1990 & WORK,1,TRONEM,1) 1991 END IF 1992 556 CONTINUE 1993 CALL GPCLOSE(LUPROP,'KEEP') 1994 END IF 1995 ELSE 1996 CALL GETDMA(ISCOOR,ONEMAT,WORK,DNULL) 1997 END IF 1998 END IF 1999 100 CONTINUE 2000 IF (KEY .EQ. 'OMAT' .OR. KEY .EQ. 'HMAT') THEN 2001 CALL DSPTSI(NBAST,TRONEM,ONEMAT) 2002 END IF 2003 IF (IPRINT .GE. 10) THEN 2004 IF (KEY .EQ. 'HMAT') THEN 2005 CALL HEADER('One-electron differentiated Hamiltonian',-1) 2006 ELSE 2007 CALL HEADER('Differentiated overlap matrix',-1) 2008 END IF 2009 CALL OUTPUT(ONEMAT,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 2010 END IF 2011 RETURN 2012 END 2013C /* Deck gettd */ 2014 SUBROUTINE GETTD(TD,NORBT,DIAGTD,DIVIDE,IPRINT) 2015#include "implicit.h" 2016#include "priunit.h" 2017 PARAMETER (D1 = 1.0D0, DP5 = 0.5D00, DP25 = 0.25D00) 2018 LOGICAL DIAGTD, DIVIDE 2019 DIMENSION TD(NORBT,NORBT) 2020C 2021 IF (DIAGTD) THEN 2022 CALL DZERO(TD,NORBT*NORBT) 2023 DO 100 I = 1, NORBT 2024 TD(I,I) = DP25 2025 100 CONTINUE 2026 ELSE 2027 FACTOR = - D1 2028 IF (DIVIDE) FACTOR = - DP5 2029 CALL DSCAL(NORBT*NORBT,FACTOR,TD,1) 2030 END IF 2031 RETURN 2032 END 2033C /* Deck getdma */ 2034 SUBROUTINE GETDMA(ISCOOR,ONEMAT,WORK,DNULL) 2035C 2036C Collect the half differentiated overlap matrice needed when using 2037C new connections. K.Ruud June-94 2038C 2039#include "implicit.h" 2040#include "dummy.h" 2041#include "priunit.h" 2042#include "iratdef.h" 2043#include "maxorb.h" 2044#include "mxcent.h" 2045 PARAMETER (D1 = 1.0D0) 2046 CHARACTER*8 LABEL 2047 DIMENSION WORK(N2BASX), ONEMAT(N2BASX) 2048 LOGICAL DNULL, FNDLAB, OPTST 2049#include "inforb.h" 2050#include "inftap.h" 2051#include "chrnos.h" 2052C 2053 CALL GPOPEN(LUPROP,'AOPROPER','UNKNOWN',' ', 2054 & 'UNFORMATTED',IDUMMY,.FALSE.) 2055C 2056C Determine label to search for 2057C 2058 LABEL = 'SQHDR'//CHRNOS(ISCOOR/100) 2059 & //CHRNOS(MOD(ISCOOR,100)/10) 2060 & //CHRNOS(MOD((MOD(ISCOOR,100)),10)) 2061 REWIND LUPROP 2062 IF (.NOT. FNDLAB(LABEL,LUPROP)) THEN 2063 WRITE(LUPRI,'(1X,3A)') 2064 & 'Label ', LABEL, ' not found on file AOPROPER' 2065 CALL QUIT('Error in GETSD') 2066 END IF 2067 CALL READT(LUPROP,N2BASX,WORK) 2068 CALL DAXPY(N2BASX,-D1,WORK,1,ONEMAT,1) 2069 IF (DNRM2(N2BASX,WORK,1) .GT. 0.0D0) THEN 2070 DNULL = .FALSE. 2071 ELSE 2072 DNULL = .TRUE. 2073 END IF 2074 CALL GPCLOSE(LUPROP,'KEEP') 2075 RETURN 2076 END 2077