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 19 SUBROUTINE CC_MOFCON2(XINT,OMEGA2, 20 * XLAMPA,XLAMHA,XLAMPB,XLAMHB,XLAMPC,XLAMHD, 21 * ISYMXLA,ISYMXLB,ISYMXLC,ISYMXLD, 22 * WORK,LWORK,IDEL,ISYMD,ISYOMEG,ISYHOP,IOPT) 23C 24C Written by Asger Halkier and Henrik Koch 3-5-95. 25C 26C Debugged By Ove Christiansen 25-7-1995 27C Generalized for cubic response calculations by C.H. in January 1997 28C Generalized for non-total symmetric integrals by C.H. in July 1998 29C 30C Purpose: To calculate the F-term's contribution to the vector 31C function and its derivatives using matrix vector routines. 32C 33C 34C IOPT=0: no symmetrization, returns: 35C F_{aibj} = (a^A i^B | b^C j^D) 36C 37C IOPT=1: symmetrization in (ai) <-> (bj), returns: 38C F_{aibj} = (a^A i^B | b^C j^D) + (b^A j^B | a^C i^D) 39C 40C IOPT=2: symmetrization in A <-> B, returns: 41C F_{aibj} = (a^A i^B | b^C j^D) + (a^B i^A | b^C j^D) 42C 43C IOPT=3: symmetrization in (ai) <-> (bj) and in A <-> B, returns: 44C F_{aibj} = (a^A i^B | b^C j^D) + (b^A j^B | a^C i^D) 45C + (a^B i^A | b^C j^D) + (b^B j^A | a^C i^D) 46C 47C N.B.: IOPT=0,2 assumes XLAMPA = XLAMPC and XLAMHB = XLAMHD 48C 49C 50C Symmetries: ISYOMEG -- OMEGA2 51C ISYMXLA -- XLAMPA, XLAMHA 52C ISYMXLB -- XLAMPB, XLAMHB 53C ISYMXLC -- XLAMPC 54C ISYMXLD -- XLAMHD 55C ISYMD -- AO IDEL 56C ISYHOP -- symmetry of integrals (all 4 indeces) 57C 58C 59C N.B. This routine assumes AO-symmetric integrals, and can therefor 60C not be used directly for calculations with London-orbitals!!! 61C 62#include "implicit.h" 63#include "priunit.h" 64#include "maxorb.h" 65 PARAMETER(ZERO = 0.0D0,ONE = 1.0D0,XMONE=-1.0D0,TWO = 2.0D0) 66 DIMENSION XINT(*),OMEGA2(*) 67 DIMENSION XLAMPA(*),XLAMHA(*),XLAMPB(*),XLAMHB(*), 68 * XLAMPC(*),XLAMHD(*) 69 DIMENSION WORK(LWORK) 70#include "ccorb.h" 71#include "symsq.h" 72#include "ccsdsym.h" 73C 74 INDEX(I,J) = MAX(I,J)*(MAX(I,J) - 3)/2 + I + J 75C 76 ISYDIS = MULD2H(ISYMD,ISYHOP) 77 ISYXAB = MULD2H(ISYMXLA,ISYMXLB) 78 ISYXCD = MULD2H(ISYMXLC,ISYMXLD) 79 80 IF (MULD2H(ISYXAB,ISYXCD) .NE. MULD2H(ISYOMEG,ISYHOP)) THEN 81 CALL QUIT('SYMMETRY MISMATCH IN CC_MOFCON2.') 82 END IF 83C 84 IF (IOPT.LT.0 .OR. IOPT.GT.3) THEN 85 CALL QUIT('CC_MOFCON2 called with an illegal value for IOPT.') 86 END IF 87 IF (IOPT.EQ.0 .OR. IOPT.EQ.2) THEN 88 IF (ISYMXLA.NE.ISYMXLC .OR. ISYMXLB.NE.ISYMXLD) THEN 89 CALL QUIT('CC_MOFCON2 called with inconsitent symmetries.') 90 END IF 91 END IF 92C 93 DO 100 ISYMG = 1,NSYM 94C 95 IF (NBAS(ISYMG) .EQ. 0) GOTO 100 96C 97 ISALBE = MULD2H(ISYMG,ISYDIS) 98 ISYMAI = MULD2H(ISALBE,ISYXAB) 99 ISYMJ = MULD2H(ISYMG,ISYMXLD) 100C 101C----------------------------------------- 102C Dynamic allocation of work space. 103C----------------------------------------- 104C 105 KSCR1 = 1 106 KSCR2 = KSCR1 + NNBST(ISALBE)*NRHF(ISYMJ) 107 KSCR3 = KSCR2 + N2BST(ISALBE) 108 KSCR4 = KSCR3 + NT1AM(ISYMAI) 109 KEND1 = KSCR4 + NT1AM(ISYMAI) 110 LWRK1 = LWORK - KEND1 111C 112 IF (LWRK1 .LT. 0) THEN 113 WRITE(LUPRI,*) 'Lwrk1 = ',LWRK1 114 CALL QUIT('Insufficient work space area in CC_MOFCON') 115 ENDIF 116C 117C-------------------------------- 118C Do first transformation. 119C-------------------------------- 120C 121 KOFF1 = IDSAOG(ISYMG,ISYDIS) + 1 122 KOFF2 = IGLMRH(ISYMG,ISYMJ) + 1 123C 124 NTALBE = MAX(NNBST(ISALBE),1) 125 NTOTG = MAX(NBAS(ISYMG),1) 126C 127 CALL DGEMM('N','N',NNBST(ISALBE),NRHF(ISYMJ),NBAS(ISYMG), 128 * ONE,XINT(KOFF1),NTALBE,XLAMHD(KOFF2),NTOTG, 129 * ZERO,WORK(KSCR1),NTALBE) 130C 131C----------------------------------- 132C Last index transformations. 133C----------------------------------- 134C 135 DO 110 J = 1,NRHF(ISYMJ) 136C 137 KOFF1 = KSCR1 + NNBST(ISALBE)*(J - 1) 138C 139 CALL CCSD_SYMSQ(WORK(KOFF1),ISALBE,WORK(KSCR2)) 140C 141 DO 120 ISYMI = 1,NSYM 142C 143 ISYMBE = MULD2H(ISYMI,ISYMXLB) 144 ISYMAL = MULD2H(ISYMBE,ISALBE) 145 ISYMA = MULD2H(ISYMAL,ISYMXLA) 146C 147 IF (LWRK1 .LT. NBAS(ISYMAL)*NRHF(ISYMI)) THEN 148 CALL QUIT('Insufficient space for 2. trf. '// 149 & 'in CC_MOFCON') 150 ENDIF 151C 152 KOFF2 = KSCR2 + IAODIS(ISYMAL,ISYMBE) 153 KOFF3 = IGLMRH(ISYMBE,ISYMI) + 1 154 KOFF4 = IGLMVI(ISYMAL,ISYMA) + 1 155 KOFF5 = KSCR3 + IT1AM(ISYMA,ISYMI) 156C 157 NTOTAL = MAX(NBAS(ISYMAL),1) 158 NTOTBE = MAX(NBAS(ISYMBE),1) 159 NTOTA = MAX(NVIR(ISYMA),1) 160C 161 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI),NBAS(ISYMBE), 162 * ONE,WORK(KOFF2),NTOTAL,XLAMHB(KOFF3),NTOTBE, 163 * ZERO,WORK(KEND1),NTOTAL) 164C 165 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI),NBAS(ISYMAL), 166 * ONE,XLAMPA(KOFF4),NTOTAL,WORK(KEND1),NTOTAL, 167 * ZERO,WORK(KOFF5),NTOTA) 168C 169C----------------------- 170C Symmetrize in A <-> B. 171C----------------------- 172 IF (IOPT .EQ. 2 .OR. IOPT .EQ. 3) THEN 173C 174 ISYMBE = MULD2H(ISYMI,ISYMXLA) 175 ISYMAL = MULD2H(ISYMBE,ISALBE) 176 ISYMA = MULD2H(ISYMAL,ISYMXLB) 177C 178 IF (LWRK1 .LT. NBAS(ISYMAL)*NRHF(ISYMI)) THEN 179 CALL QUIT('Insufficient space for 2. trf. '// 180 & 'in CC_MOFCON') 181 ENDIF 182C 183 KOFF2 = KSCR2 + IAODIS(ISYMAL,ISYMBE) 184 KOFF3 = IGLMRH(ISYMBE,ISYMI) + 1 185 KOFF4 = IGLMVI(ISYMAL,ISYMA) + 1 186 KOFF5 = KSCR3 + IT1AM(ISYMA,ISYMI) 187C 188 NTOTAL = MAX(NBAS(ISYMAL),1) 189 NTOTBE = MAX(NBAS(ISYMBE),1) 190 NTOTA = MAX(NVIR(ISYMA),1) 191C 192 CALL DGEMM('N','N',NBAS(ISYMAL),NRHF(ISYMI), 193 * NBAS(ISYMBE),ONE,WORK(KOFF2),NTOTAL, 194 * XLAMHA(KOFF3),NTOTBE,ZERO,WORK(KEND1), 195 * NTOTAL) 196C 197 CALL DGEMM('T','N',NVIR(ISYMA),NRHF(ISYMI), 198 * NBAS(ISYMAL),ONE,XLAMPB(KOFF4),NTOTAL, 199 * WORK(KEND1),NTOTAL,ONE,WORK(KOFF5),NTOTA) 200C 201 ENDIF 202C 203 204 120 CONTINUE 205C 206C-------------------------------------------------- 207C Storing the result in the omega2-array. 208C-------------------------------------------------- 209C 210 ISYMB = MULD2H(ISYMD,ISYMXLC) 211 ISYMBJ = MULD2H(ISYMB,ISYMJ) 212C 213 DO 130 B = 1,NVIR(ISYMB) 214C 215 NBJ = IT1AM(ISYMB,ISYMJ) + NVIR(ISYMB)*(J - 1) + B 216 NDB = IGLMVI(ISYMD,ISYMB) + NBAS(ISYMD)*(B - 1) 217 * + IDEL - IBAS(ISYMD) 218C 219 CALL DZERO(WORK(KSCR4),NT1AM(ISYMAI)) 220C 221 XLB = XLAMPC(NDB) 222C 223 CALL DAXPY(NT1AM(ISYMAI),XLB,WORK(KSCR3),1,WORK(KSCR4),1) 224C 225 IF (ISYMBJ .EQ. ISYMAI) THEN 226C 227 NTOTAI = NBJ 228C 229 IF (IOPT .EQ. 1 .OR. IOPT .EQ. 3) THEN 230 NTOTAI = NT1AM(ISYMAI) 231 WORK(KSCR4+NBJ-1) = TWO*WORK(KSCR4+NBJ-1) 232 ENDIF 233C 234 DO 140 NAI = 1,NTOTAI 235C 236 NAIBJ = IT2AM(ISYMAI,ISYMBJ) + INDEX(NAI,NBJ) 237C 238 OMEGA2(NAIBJ) = OMEGA2(NAIBJ) + WORK(KSCR4+NAI-1) 239C 240 140 CONTINUE 241C 242 ENDIF 243C 244 IF (ISYMAI .LT. ISYMBJ) THEN 245C 246 DO 150 NAI = 1,NT1AM(ISYMAI) 247C 248 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 249 * + NT1AM(ISYMAI)*(NBJ - 1) + NAI 250C 251 OMEGA2(NAIBJ) = OMEGA2(NAIBJ) + WORK(KSCR4+NAI-1) 252C 253 150 CONTINUE 254C 255 ENDIF 256C 257 IF ((ISYMBJ.LT.ISYMAI).AND.(IOPT.EQ.1.OR.IOPT.EQ.3))THEN 258C 259 DO 160 NAI = 1,NT1AM(ISYMAI) 260C 261 NAIBJ = IT2AM(ISYMAI,ISYMBJ) 262 * + NT1AM(ISYMBJ)*(NAI - 1) + NBJ 263C 264 OMEGA2(NAIBJ) = OMEGA2(NAIBJ) + WORK(KSCR4+NAI-1) 265C 266 160 CONTINUE 267C 268 ENDIF 269C 270 130 CONTINUE 271C 272 110 CONTINUE 273C 274 100 CONTINUE 275C 276 RETURN 277 END 278