1! 2! Dalton, a molecular electronic structure program 3! Copyright (C) by the authors of Dalton. 4! 5! This program is free software; you can redistribute it and/or 6! modify it under the terms of the GNU Lesser General Public 7! License version 2.1 as published by the Free Software Foundation. 8! 9! This program is distributed in the hope that it will be useful, 10! but WITHOUT ANY WARRANTY; without even the implied warranty of 11! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 12! Lesser General Public License for more details. 13! 14! If a copy of the GNU LGPL v2.1 was not distributed with this 15! code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html. 16! 17! 18C 19C /* Deck dftdns */ 20 SUBROUTINE DFTDNS(DMAT,WORK,LWORK,IPRINT) 21C 22C T. Helgaker Feb 01 23C 24#include "implicit.h" 25#include "priunit.h" 26#include "dummy.h" 27#include "mxcent.h" 28#include "iratdef.h" 29#include "inforb.h" 30C 31 DIMENSION DMAT(NBAST,NBAST), WORK(LWORK) 32C 33#include "inftap.h" 34#include "infvar.h" 35#include "nuclei.h" 36#include "dftcom.h" 37C 38 IF (NASHT .GT. 0) 39 & CALL QUIT('DFTDNS ERROR: open-shell DFT not implemented') 40 IF (NCMOT .GT. LWORK) CALL STOPIT('DFTDNS',' ',NCMOT,LWORK) 41 REWIND LUSIFC 42 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 43 READ (LUSIFC) 44 READ (LUSIFC) 45 CALL READI(LUSIFC,IRAT*NCMOT,WORK) 46 JKEEP = JWOPSY 47 JWOPSY = 1 48 CALL FCKDEN(.TRUE.,.FALSE.,DMAT,DUMMY,WORK,DUMMY,WORK,LWORK) 49 JWOPSY = JKEEP 50 IF (IPRINT.GT.100) THEN 51 CALL HEADER('AO density matrix in DFTDNS ',-1) 52 CALL OUTPUT(DMAT,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 53 END IF 54 RETURN 55 END 56C /* Deck getden */ 57 SUBROUTINE GETDEN(DMAT,WORK,LWORK,NODC,NODV,IPRINT) 58C 59C T. Helgaker Sep 1999 60C 61#include "implicit.h" 62#include "priunit.h" 63C 64 PARAMETER (DP5 = 0.50D0, D1 = 1.0D0) 65C 66#include "inforb.h" 67C 68 LOGICAL NODC, NODV 69 DIMENSION DMAT(NBAST,NBAST), WORK(LWORK) 70C 71 IF (IPRINT .GE. 5) CALL TITLER('Output from GETDEN','*',103) 72C 73 KDSO = 1 74 KDNS = KDSO + NNBAST 75 KLAST = KDNS + NNBASX 76 LWRK = LWORK - KLAST + 1 77 IF (KLAST.GT.LWORK) CALL STOPIT('GETDEN','QDRDSO',KLAST,LWORK) 78 CALL QDRDSO(WORK(KDSO),WORK(KLAST),LWRK,IPRINT,NODC,NODV) 79 CALL QDRDAO(WORK(KDNS),WORK(KDSO),NBAST,IPRINT) 80C 81 IJ = 0 82 DO 100 I = 1, NBAST 83 DO 100 J = 1, I 84 FAC = D1 85 IF (I.NE.J) FAC = DP5 86 DMAT(I,J) = FAC*WORK(KDNS+IJ) 87 DMAT(J,I) = FAC*WORK(KDNS+IJ) 88 IJ = IJ + 1 89 100 CONTINUE 90 IF (IPRINT .GT. 5) THEN 91 CALL HEADER('Total density matrix from GETDEN',-1) 92 CALL OUTPUT(DMAT,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 93 END IF 94 RETURN 95 END 96C 97C /* Deck qdrdso */ 98 SUBROUTINE QDRDSO(DSO,WORK,LWORK,IPRINT,NODC,NODV) 99#include "implicit.h" 100#include "priunit.h" 101#include "inforb.h" 102 LOGICAL NODC, NODV 103 DIMENSION DSO(NNBAST), WORK(LWORK) 104C 105 KCMO = 1 106 KDV = KCMO + NCMOT 107 KLAST = KDV + NNASHX 108 IF (KLAST .GT. LWORK) CALL STOPIT('QDRDSO',' ',KLAST,LWORK) 109 CALL QDRDS1(WORK(KCMO),WORK(KDV),DSO,IPRINT,NODC,NODV) 110 RETURN 111 END 112C /* Deck qdrds1 */ 113 SUBROUTINE QDRDS1(CMO,DV,DSO,IPRINT,NODC,NODV) 114C 115#include "implicit.h" 116 PARAMETER (D0 = 0.0D0, TWO = 2.0D0, HALF = 0.5D0) 117#include "maxorb.h" 118#include "inforb.h" 119 LOGICAL NODC, NODV 120 DIMENSION CMO(*), DV(*), DSO(NNBAST) 121#include "iratdef.h" 122#include "mxcent.h" 123#include "priunit.h" 124C 125#include "abainf.h" 126#include "inftap.h" 127 INTEGER R, S, RS, U, V, UV 128C 129C 130C ***** Read input from LUSIFC ***** 131C 132 REWIND LUSIFC 133 CALL MOLLAB('SIR IPH ',LUSIFC,LUPRI) 134 READ (LUSIFC) 135 READ (LUSIFC) NISHT,NASHT,NOCCT,NORBT,NBAST,NCONF,NWOPT,NWOPH, 136 * NCMOT,NNASHX,NNASHY,NNORBT,N2ORBT 137 NSSHT = NORBT - NOCCT 138 CALL READI(LUSIFC,IRAT*NCMOT,CMO) 139 READ (LUSIFC) 140 IF (NASHT .GT. 0) THEN 141 CALL READI(LUSIFC,IRAT*NNASHX,DV) 142 ELSE 143 READ (LUSIFC) 144 END IF 145C 146C ***** Print Section ***** 147C 148 IF (IPRINT .GT. 05) THEN 149 WRITE (LUPRI, '(//A/)') ' ----- SUBROUTINE QDRDS1 ------' 150 WRITE (LUPRI, '(A,8I5)') ' NISH ', (NISH(I),I = 1,NSYM) 151 WRITE (LUPRI, '(A,8I5)') ' NASH ', (NASH(I),I = 1,NSYM) 152 WRITE (LUPRI, '(A,8I5)') ' NOCC ', (NOCC(I),I = 1,NSYM) 153 WRITE (LUPRI, '(A,8I5)') ' NORB ', (NORB(I),I = 1,NSYM) 154 WRITE (LUPRI, '(A,8I5)') ' NBAS ', (NBAS(I),I = 1,NSYM) 155 IF (IPRINT .GE. 10) THEN 156 CALL HEADER('Occupied molecular orbitals',0) 157 IEND = 0 158 DO 1000 ISYM = 1,NSYM 159 IF (NBAS(ISYM) .EQ. 0) GOTO 1000 160 IF (NOCC(ISYM) .EQ. 0) GOTO 1100 161 WRITE (LUPRI, '(//,A,I5,/)') ' Symmetry ', ISYM 162 IENDI = IEND 163 DO 1200 I = 1, NOCC(ISYM) 164 WRITE (LUPRI,'(/,A,I5,/)') ' Molecular orbital ', I 165 WRITE (LUPRI,'(6F12.6)') (CMO(IENDI+J),J=1,NBAS(ISYM)) 166 IENDI = IENDI + NBAS(ISYM) 1671200 CONTINUE 1681100 CONTINUE 169 IEND = IEND + NORB(ISYM)*NBAS(ISYM) 1701000 CONTINUE 171 CALL HEADER('Active density matrix (MO basis)',-1) 172 CALL OUTPAK(DV,NASHT,1,LUPRI) 173 END IF 174 END IF 175C 176C ***** Construct contravariant SO matrices ***** 177C 178 ISEND = 0 179 ICEND = 0 180 DO 110 ISYM = 1,NSYM 181 NORBI = NORB(ISYM) 182 NISHI = NISH(ISYM) 183 NASHI = NASH(ISYM) 184 IASHI = IASH(ISYM) 185 NBASI = NBAS(ISYM) 186 IF (NBASI .EQ. 0) GOTO 120 187 IF (NOCC(ISYM) .EQ. 0) THEN 188 CALL DZERO(DSO(ISEND+1),NNBAS(ISYM)) 189 GO TO 120 190 END IF 191 RS = 0 192 DO 100 R = 1, NBASI 193 DO 200 S = 1,R 194 RS = RS + 1 195C 196 DTRS = D0 197C 198C (I) Inactive contribution 199C 200 IF (NISHI .GT. 0) THEN 201 ICENDI = ICEND 202 DO 300 I = 1, NISHI 203 DTRS = DTRS + CMO(ICENDI+R)*CMO(ICENDI+S) 204 ICENDI = ICENDI + NBASI 205 300 CONTINUE 206 DTRS = DTRS + DTRS 207 END IF 208 IF (NODC) DTRS = D0 209C 210C (II) Active contribution 211C 212 IF (.NOT. NODV) THEN 213 IF (NASHI .GT. 0) THEN 214 UV = ((IASHI + 1)*(IASHI + 2))/2 215 IDVEND = ICEND + NISHI*NBASI 216 ICENDU = IDVEND 217 DO 400 U = 1,NASHI 218 ICENDV = IDVEND 219 DO 410 V = 1, U 220 DUV = DV(UV) 221 IF (ABS(DUV) .GT. D0) THEN 222 TEMP = CMO(ICENDU+R)*CMO(ICENDV+S) 223 IF (U .NE. V) TEMP = TEMP 224 * + CMO(ICENDU+S)*CMO(ICENDV+R) 225 DTRS = DTRS + DUV*TEMP 226 END IF 227 UV = UV + 1 228 ICENDV = ICENDV + NBASI 229 410 CONTINUE 230 UV = UV + IASHI 231 ICENDU = ICENDU + NBASI 232 400 CONTINUE 233 END IF 234 END IF 235 IF (R .NE. S) DTRS = DTRS + DTRS 236 DSO(ISEND+RS) = DTRS 237C 238 200 CONTINUE 239 100 CONTINUE 240C 241C ***** Print Section ***** 242C 243 IF (IPRINT .GE. 10) THEN 244 WRITE (LUPRI,'(1X,A,I5)') ' Symmetry', ISYM 245 CALL HEADER('Total density matrix (SO basis)',-1) 246 CALL OUTPAK(DSO(ISEND+1),NBASI,1,LUPRI) 247 END IF 248120 CONTINUE 249 ISEND = ISEND + (NBASI*(NBASI + 1))/2 250 ICEND = ICEND + NORBI*NBASI 251110 CONTINUE 252 RETURN 253 END 254C /* Deck qdrdao */ 255 SUBROUTINE QDRDAO(DENMAT,DSO,NBAST,IPRINT) 256C 257#include "implicit.h" 258#include "priunit.h" 259#include "maxaqn.h" 260#include "maxorb.h" 261#include "mxcent.h" 262 DIMENSION DSO(*), DENMAT(*) 263#include "shells.h" 264#include "pincom.h" 265#include "symmet.h" 266 267C 268 IF (IPRINT .GT. 10) CALL HEADER('Subroutine QDRDAO',-1) 269C 270C Loop over all irreps in molecule 271C 272 ISOFF = 0 273 ISTR = 1 274 NNBASX = NBAST*(NBAST + 1)/2 275 CALL DZERO(DENMAT,NNBASX) 276 DO 100 IREP = 0, MAXREP 277 NORBI = NAOS(IREP+1) 278 IF (NORBI .EQ. 0) GOTO 110 279 DO 200 I = ISTR,ISTR + NORBI - 1 280 IA = IAND(ISHFT(IPIND(I),-16),65535) 281 NA = IAND(ISHFT(IPIND(I), -8), 255) 282 IOFF = KSTRT(IA) 283 MULA = ISTBAO(IA) 284 INDA = IOFF + NA 285 DO 300 J = ISTR,I 286 IB = IAND(ISHFT(IPIND(J),-16),65535) 287 NB = IAND(ISHFT(IPIND(J), -8), 255) 288 JOFF = KSTRT(IB) 289 NHKTB = NHKT(IB) 290 KHKTB = KHKT(IB) 291 MULB = ISTBAO(IB) 292 MAB = IOR(MULA,MULB) 293 KAB = IAND(MULA,MULB) 294 HKAB = FMULT(KAB) 295 ISOFF = ISOFF + 1 296 DSYMIJ = DSO(ISOFF) 297 INDB = JOFF + NB - KHKTB 298 DO 400 ISYMOP = 0, MAXOPR 299 IF (IAND(ISYMOP,MAB) .NE. 0) GOTO 400 300 INDB = INDB + KHKTB 301C 302C Weight and parity factor 303C 304 FAC = HKAB* 305 * PT(IAND(ISYMOP,IEOR(IREP,ISYMAO(NHKTB,NB)))) 306 INDM = MAX(INDA,INDB) 307 IND = (INDM*(INDM - 3))/2 + INDA + INDB 308 DENMAT(IND) = DENMAT(IND) + FAC*DSYMIJ 309400 CONTINUE 310300 CONTINUE 311200 CONTINUE 312110 CONTINUE 313 ISTR = ISTR + NORBI 314100 CONTINUE 315 IF (IPRINT .GT. 10) THEN 316 CALL HEADER('Total density matrix (sym. distinct AO basis)',-1) 317 CALL OUTPAK(DENMAT,NBAST,1,LUPRI) 318 END IF 319 RETURN 320 END 321C /* Deck dftdnsab */ 322 SUBROUTINE DFTDNSAB(DMATA,DMATB,WORK,LWORK,IPRINT) 323C 324C DFTDNS adaptation for open shell case ZR 325C 326#include "implicit.h" 327 DIMENSION DMATA(*), DMATB(*) 328 DIMENSION WORK(LWORK) 329 INTEGER IPRINT 330#include "inforb.h" 331#include "inftap.h" 332#include "dummy.h" 333#include "priunit.h" 334 PARAMETER ( D1 = 1.0D0, DP5 =0.5D0) 335 LOGICAL CLOSED 336C 337C used from infvar.h: JWOPSY 338#include "infvar.h" 339C 340 CALL QENTER('DFTDNSAB') 341C 342 KFREE = 1 343 LFREE = LWORK 344 CALL MEMGET('REAL',KCMO,NCMOT,WORK,KFREE,LFREE) 345 CALL MEMGET('REAL',KUDV,NNASHX,WORK,KFREE,LFREE) 346C 347 CLOSED=LUSIFC.LT.0 348 IF (CLOSED) CALL GPOPEN( 349 & LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED',IDUMMY,.FALSE.) 350 REWIND(LUSIFC) 351 CALL MOLLAB(LBSIFC,LUSIFC,LUPRI) 352 READ(LUSIFC) 353 READ(LUSIFC) 354 CALL READT(LUSIFC,NCMOT,WORK(KCMO)) 355 READ(LUSIFC) 356 CALL READT(LUSIFC,NNASHX,WORK(KUDV)) 357C 358 JKEEP = JWOPSY 359 JWOPSY = 1 360 IF (NASHT.EQ.0) THEN 361 CALL FCKDEN(.TRUE.,.FALSE.,DMATB,DUMMY,WORK(KCMO),DUMMY, 362 & WORK(KFREE),LFREE) 363 CALL DZERO(DMATA,N2BASX) 364 ELSE 365 CALL FCKDEN(.TRUE.,.TRUE.,DMATB,DMATA,WORK(KCMO),WORK(KUDV), 366 & WORK(KFREE),LFREE) 367 END IF 368 JWOPSY = JKEEP 369 CALL DSCAL(N2BASX,DP5,DMATB,1) 370 CALL DAXPY(N2BASX,D1,DMATB,1,DMATA,1) 371C 372 IF (CLOSED) CALL GPCLOSE(LUSIFC,'KEEP') 373C 374 IF (IPRINT.GT.100) THEN 375 CALL HEADER('AO alpha density matrix in DFTDNSAB ',-1) 376 CALL OUTPUT(DMATA,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 377 CALL HEADER('AO beta density matrix in DFTDNSAB ',-1) 378 CALL OUTPUT(DMATB,1,NBAST,1,NBAST,NBAST,NBAST,1,LUPRI) 379 END IF 380C 381 CALL MEMREL('DFTDNSAB',WORK,1,1,KFREE,LFREE) 382C 383 CALL QEXIT('DFTDNSAB') 384C 385 END 386 387C -- end of dft_den.F -- 388