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 cc_aodens3 */ 20 SUBROUTINE CC_AODENS3(XLAMDP1,ISYMP1,XLAMDH2,ISYMH2, 21 & XLAMDP3,ISYMP3,XLAMDH4,ISYMH4, 22 & DENS,ISYDEN,ICORE,IOPT,WORK,LWORK) 23C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 24C 25C Calculate (special) AO-density matrix used in contructing 26C the (special) AO Fock matrix. 27C 28C IOPT = 1 : XLAMDP1 x XLAMDH2 29C IOPT = 2 : XLAMDP1 x XLAMDH2 + XLAMDP3 x XLAMDH4 30C 31C Sonia Coriani 10-Feb-1999, based on AODENS2 32C 33C Careful: No special treatment of ICORE. 34C ORDER IS IMPORTANT: always P1*H2 (+ P3*H4) 35C Debug 12.8.99 OK 36C Reinsert in AODENS2 at some point 37C&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& 38#include "implicit.h" 39 PARAMETER (ZERO = 0.0D0, ONE = 1.0D0, TWO = 2.0D0) 40 DIMENSION XLAMDP1(*), XLAMDH2(*), XLAMDP3(*), XLAMDH4(*) 41 DIMENSION DENS(*), WORK(LWORK) 42#include "ccorb.h" 43#include "ccsdinp.h" 44#include "ccsdsym.h" 45#include "priunit.h" 46#include "dummy.h" 47C 48* 49* Symmetry test 50* 51 ISYTOT1 = MULD2H(ISYMP1,ISYMH2) 52 IF (ISYTOT1.NE.ISYDEN) 53 * CALL QUIT('Symmetry mismatch 1 in AODENS3' ) 54 IF (IOPT.EQ.2) THEN 55 ISYTOT2 = MULD2H(ISYMP3,ISYMH4) 56 IF (ISYTOT1.NE.ISYTOT2) 57 * CALL QUIT('Symmetry mismatch 1 in AODENS3' ) 58 END IF 59* 60 IF (IOPT.GE.1) THEN 61 CALL DZERO(DENS,N2BST(ISYDEN)) 62 END IF 63* 64 DO ISYMK = 1,NSYM 65C 66 IF (IOPT.GE.1) THEN 67C 68 ISYMA = MULD2H(ISYMP1,ISYMK) ! XLAMDP1 69 ISYMB = MULD2H(ISYMH2,ISYMK) ! XLAMDH2 70C 71 KOFF1 = 1 + IGLMRH(ISYMA,ISYMK) ! XLAMDP1 72 KOFF2 = 1 + IGLMRH(ISYMB,ISYMK) ! XLAMDH2 73 KOFF3 = 1 + IAODIS(ISYMA,ISYMB) 74 NBASA = MAX(NBAS(ISYMA),1) 75 NBASB = MAX(NBAS(ISYMB),1) 76C 77 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYMK),ONE, 78 * XLAMDP1(KOFF1),NBASA,XLAMDH2(KOFF2),NBASB,ONE, 79 * DENS(KOFF3),NBASA) 80C 81 END IF 82C 83 IF (IOPT.EQ.2) THEN !2nd contribution 84C 85 ISYMA = MULD2H(ISYMP3,ISYMK) ! XLAMDP3 86 ISYMB = MULD2H(ISYMH4,ISYMK) ! XLAMDH4 87C 88 KOFF1 = 1 + IGLMRH(ISYMA,ISYMK) ! XLAMDP3 89 KOFF2 = 1 + IGLMRH(ISYMB,ISYMK) ! XLAMDH4 90 KOFF3 = 1 + IAODIS(ISYMA,ISYMB) 91 NBASA = MAX(NBAS(ISYMA),1) 92 NBASB = MAX(NBAS(ISYMB),1) 93C 94 CALL DGEMM('N','T',NBAS(ISYMA),NBAS(ISYMB),NRHF(ISYMK),ONE, 95 * XLAMDP3(KOFF1),NBASA,XLAMDH4(KOFF2),NBASB,ONE, 96 * DENS(KOFF3),NBASA) 97C 98 END IF !result has been added to previous (beta=1) 99C 100 END DO 101C 102C 103C----------------------------- 104C Include frozen orbitals. 105C----------------------------- 106C 107 IF ( (FROIMP.OR.FROEXP) .AND. (ICORE .EQ. 1) ) THEN 108C 109 IF (IOPT.NE.0) THEN 110 WRITE (LUPRI,*) 111 * 'CC_AODENS3: ICORE=1 not yet available for IOPT>0.' 112 CALL QUIT( 113 * 'CC_AODENS3: ICORE=1 not yet available for IOPT>0.') 114 END IF 115C 116 IF (LWORK .LT. NLAMDS) THEN 117 CALL QUIT('Insufficient space in CCSD_AODENS') 118 ENDIF 119C 120C------------------------------------------------- 121C Read MO-coefficients from interface file. 122C------------------------------------------------- 123C 124 LUSIFC = -1 125 CALL GPOPEN(LUSIFC,'SIRIFC','OLD',' ','UNFORMATTED', 126 * IDUMMY,.FALSE.) 127 REWIND LUSIFC 128C 129 CALL MOLLAB('TRCCINT ',LUSIFC,LUPRI) 130 READ (LUSIFC) 131C 132 READ (LUSIFC) 133 READ (LUSIFC) (WORK(I), I=1,NLAMDS) 134C 135 CALL GPCLOSE(LUSIFC,'KEEP') 136C 137C------------------------------------------------------- 138C Add contribution from frozen occupied orbitals. 139C------------------------------------------------------- 140C 141 KOFF1 = 0 142 KOFF2 = 0 143 DO ISYMK = 1,NSYM 144C 145CCH the following has been changed because ISYMP0 is not defined...? 146CCH ISYMA = MULD2H(ISYMP0,ISYMK) 147CCH ISYMB = MULD2H(ISYMP0,ISYMK) 148CCH 149CCH Sonia is this correct? 150 ISYMA = ISYMK 151 ISYMB = ISYMK 152CCH 153C 154 DO II = 1,NRHFFR(ISYMK) 155C 156 K = KFRRHF(II,ISYMK) 157C 158 DO B = 1,NBAS(ISYMB) 159 DO A = 1,NBAS(ISYMA) 160C 161 NAK = KOFF1 + NBAS(ISYMA)*(K - 1) + A 162 NBK = KOFF1 + NBAS(ISYMB)*(K - 1) + B 163 NAB = KOFF2 + NBAS(ISYMA)*(B - 1) + A 164C 165 DENS(NAB) = DENS(NAB) + WORK(NAK)*WORK(NBK) 166C 167 END DO 168 END DO 169C 170 END DO 171C 172 KOFF1 = KOFF1 + NBAS(ISYMK)*NORBS(ISYMK) 173 KOFF2 = KOFF2 + NBAS(ISYMA)*NBAS(ISYMB) 174C 175 END DO 176C 177 ENDIF 178C 179 END 180