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*=====================================================================* 20 SUBROUTINE CCFBINT1(ITRAN, LISTL, IDLSTL, LISTR, IDLSTR, 21 & XAINT, YAINT, MAINT, 22 & ZETA1, ZETA2, ISYCTR, 23 & TA1AMP, TA2AMP, ISYMTA, 24 & XLAMDP, XLAMDH, ISYLAM, 25 & XLAMDPQ, XLAMDHQ, ISYHOP, 26 & BZDENA, BZDENAB, 27 & FILPQAMO,LUPQAMO, IADRPQAMO,IADRPQA, 28 & FILPQA0, LUPQA0, IADRPQA0, IADRPQAI0, 29 & FILPQA1, LUPQA1, IADRPQA1, IADRPQAI1, 30 & LRELAX, LTWOEL, WORK, LWORK) 31*---------------------------------------------------------------------* 32* Purpose: 33* 34* Precalculate some intermediates for FbTa result vector depending 35* on ZETA, T^A and IOPER: 36* -- the XA and YA intermediates (in memory, allocated outside) 37* -- the effective density BZDENA for the rho^BZA intermediate (disk) 38* -- the effective density BZDENAB for the rho^BZQA intermediate (memory) 39* -- the PA and QA intermediates, both with 4 MO indices and 40* with 1 index in AO (0 and bar) 41* 42* Note: flags LRELAX & LTWOEL not yet carried through...? 43* 44* Sonia Coriani, February 1999 (based on CCETAINT1) 45* Version: 08/10-1999 46*---------------------------------------------------------------------* 47 IMPLICIT NONE 48#include "priunit.h" 49#include "ccsdinp.h" 50#include "ccsdsym.h" 51#include "maxorb.h" 52#include "ccorb.h" 53 54 LOGICAL LOCDBG 55 PARAMETER (LOCDBG = .FALSE.) 56 57 LOGICAL LRELAX, LTWOEL 58 CHARACTER*(*) LISTL, LISTR 59c 60 CHARACTER*(*) FILPQAMO, FILPQA0, FILPQA1 61c 62 INTEGER ITRAN, IDLSTL, ISYCTR, IDLSTR, ISYMTA 63 INTEGER ISYLAM, ISYHOP, LWORK 64c 65 INTEGER LUPQAMO, IADRPQA, IADRPQAMO(MXCORB_CC,*) 66 INTEGER LUPQA0, IADRPQAI0, IADRPQA0(MXCORB_CC,*) 67 INTEGER LUPQA1, IADRPQAI1, IADRPQA1(MXCORB_CC,*) 68 69#if defined (SYS_CRAY) 70 REAL XAINT(*), YAINT(*), MAINT(*), ZETA1(*), ZETA2(*) 71 REAL BZDENA(*), BZDENAB(*) 72 REAL XLAMDP(*), XLAMDH(*), XLAMDPQ(*), XLAMDHQ(*) 73 REAL TA1AMP(*), TA2AMP(*), WORK(*) 74 REAL DUMMY, ONE, TWO 75#else 76 DOUBLE PRECISION XAINT(*), YAINT(*), MAINT(*), ZETA1(*), ZETA2(*) 77 DOUBLE PRECISION BZDENA(*), BZDENAB(*) 78 DOUBLE PRECISION XLAMDP(*), XLAMDH(*), XLAMDPQ(*), XLAMDHQ(*) 79 DOUBLE PRECISION TA1AMP(*), TA2AMP(*), WORK(*) 80 DOUBLE PRECISION DUMMY, ONE, TWO 81#endif 82 PARAMETER (ONE = 1.0D0, TWO = 2.0D0) 83 84 CHARACTER MODEL*(10) 85 INTEGER IOPT, ISYINT, ISYINTQ, KCHI, KCHIQ, KEND1, LWRK1, IDEL 86 INTEGER IDUMMY, KZ1A, KCHIA 87 INTEGER ILLSTOLD 88 SAVE ILLSTOLD 89 INTEGER IRLSTOLD 90 SAVE IRLSTOLD 91 92*---------------------------------------------------------------------* 93* do some tests and set symmetries: 94*---------------------------------------------------------------------* 95* if ITRAN=1, make sure that everything will be initialzed: 96* 97 IF (ITRAN.EQ.1) THEN 98 ILLSTOLD = IDLSTL - 1 99 IRLSTOLD = IDLSTR - 1 100 END IF 101 102* set symmetries for intermediates: 103 104 ISYINT = MULD2H(ISYCTR,ISYMTA) !ZxTA 105 ISYINTQ = MULD2H(ISYINT,ISYHOP) !ZxTAxOperB 106 107*---------------------------------------------------------------------* 108* if only left vector (IDLSTL) has changed, 109* read new multipliers Z1+Z2 into memory and square Z2 up: 110*---------------------------------------------------------------------* 111 IF ((IDLSTL .NE. ILLSTOLD).AND.(IDLSTR.EQ.IRLSTOLD)) THEN 112 113 IF (LWORK .LT. NT2AM(ISYCTR)) THEN 114 CALL QUIT('Insufficient memory in CCFBINT1 (a1)') 115 END IF 116 117 IOPT = 3 118 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,ZETA1,WORK) 119 120 CALL CC_T2SQ(WORK,ZETA2,ISYCTR) 121 122 IF (LOCDBG) THEN 123 WRITE (LUPRI,*) 'CCFBINT1> Only NEW Zeta vector' 124 WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:' 125 CALL CC_PRP(WORK,WORK,ISYCTR,0,1) 126 WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:' 127 CALL CC_PRP(WORK,TA2AMP,ISYMTA,0,1) 128 END IF 129 END IF 130*---------------------------------------------------------------------* 131* if only right vector (IDLSTR) has changed, 132* read new response amplitudes TA1+TA2 into memory: 133*---------------------------------------------------------------------* 134 135 IF ((IDLSTL .EQ. ILLSTOLD).AND.(IDLSTR.NE.IRLSTOLD)) THEN 136 137 IF (LWORK .LT. NT2AM(ISYMTA)) THEN 138 CALL QUIT('Insufficient memory in CCFBINT1 (a2)') 139 END IF 140 141 IOPT = 3 142 CALL CC_RDRSP(LISTR,IDLSTR,ISYMTA,IOPT,MODEL,TA1AMP,TA2AMP) 143 144 IF (LOCDBG) THEN 145 WRITE (LUPRI,*) 'CCFBINT1> Only NEW T^A vector' 146 WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:' 147 CALL CC_PRP(WORK,WORK,ISYCTR,0,1) 148 WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:' 149 CALL CC_PRP(WORK,TA2AMP,ISYMTA,0,1) 150 END IF 151* 152* scale with 2 the tA_2 part 153* 154 IF ( (IOPT.EQ.3) .AND. (.NOT. LISTR(1:2).EQ.'R0') ) THEN 155 CALL CCLR_DIASCL(TA2AMP,TWO,ISYMTA) 156 END IF 157* 158 END IF 159*---------------------------------------------------------------------* 160* if both vectors (IDLSTR,IDLSTL) have changed, 161* read new Zeta and response amplitudes into memory: 162*---------------------------------------------------------------------* 163 164 IF ((IDLSTL .NE. ILLSTOLD).AND.(IDLSTR.NE.IRLSTOLD)) THEN 165 166c IF (LWORK .LT. (NT2AM(ISYCTR)+NT2AM(ISYMTA))) THEN 167 IF (LWORK .LT. NT2AM(ISYCTR)) THEN 168 CALL QUIT('Insufficient memory in CCFBINT1 (a3)') 169 END IF 170 171 IOPT = 3 !(both singles and doubles) 172 CALL CC_RDRSP(LISTR,IDLSTR,ISYMTA,IOPT,MODEL,TA1AMP,TA2AMP) 173 174c KSCR = 1 175c KEND1 = KSCR + NT2AM(ISYCTR) 176 177 IOPT = 3 178 CALL CC_RDRSP(LISTL,IDLSTL,ISYCTR,IOPT,MODEL,ZETA1,WORK) 179 180 CALL CC_T2SQ(WORK,ZETA2,ISYCTR) 181 182 IF (LOCDBG) THEN 183 WRITE (LUPRI,*) 'CCFBINT1> Both NEW T^A and Zeta vectors' 184 WRITE (LUPRI,*) 'CCFBINT1> the zeta2 vector:(packed)' 185 CALL CC_PRP(WORK,WORK,ISYCTR,0,1) 186 WRITE (LUPRI,*) 'CCFBINT1> the tA2amp vector:' 187 CALL CC_PRP(WORK(1),TA2AMP,ISYMTA,0,1) 188 END IF 189* 190* scale with 2 the tA_2 part 191* 192 IF ( (IOPT.EQ.3) .AND. (.NOT. LISTR(1:2).EQ.'R0') ) THEN 193 CALL CCLR_DIASCL(TA2AMP,TWO,ISYMTA) 194 END IF 195* 196 END IF 197*---------------------------------------------------------------------* 198* if either left or right vector changed, calculate 199* new XA, YA and MA intermediates: 200*---------------------------------------------------------------------* 201 IF (( IDLSTL .NE. ILLSTOLD ).OR.( IDLSTR .NE. IRLSTOLD )) THEN 202 203 CALL CC_XI(XAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK) 204 205 CALL CC_YI(YAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK) 206 207 IF (CCSD .AND. LRELAX) THEN 208 CALL CC_MI(MAINT,ZETA2,ISYCTR,TA2AMP,ISYMTA,WORK,LWORK) 209 END IF 210 211 IF (LOCDBG) THEN 212 WRITE (LUPRI,'(//A)') 'CCFBINT1> XA-intermediate:' 213 WRITE (LUPRI,'(5G15.6)') (XAINT(I),I=1,NMATIJ(ISYINT)) 214 WRITE (LUPRI,'(//A)') 'CCFBINT1> YA-intermediate:' 215 WRITE (LUPRI,'(5G15.6)') (YAINT(I),I=1,NMATAB(ISYINT)) 216 IF (CCSD .AND. LRELAX) THEN 217 WRITE (LUPRI,'(//A)') 'CCFBINT1> MA-intermediate:' 218 WRITE (LUPRI,'(5G15.6)') (MAINT(I),I=1,N3ORHF(ISYINT)) 219 END IF 220 END IF 221 222 END IF 223 224*---------------------------------------------------------------------* 225* calculate Chi^A matrices (calN)=(breve-Zeta1 - XA) intermediate 226* transform Zeta^ci to ZetaA^ki with T^A_ck 227* (cheap, do always!) 228*---------------------------------------------------------------------* 229 IF (CCSD) THEN 230 KCHIA = 1 231 KEND1 = KCHIA + NMATIJ(ISYINT) 232 LWRK1 = LWORK - KEND1 233 234 IF (LWRK1 .LT. NT2AM(ISYCTR)) THEN 235 CALL QUIT('Insufficient memory in CCFBINT1 (a)') 236 END IF 237 CALL CCLT_Z1A(ZETA1,ISYCTR,TA1AMP,ISYMTA, ISYINT,WORK(KCHIA) ) 238 IF (LOCDBG) THEN 239 WRITE (LUPRI,'(//A)') 'CCFBINT1> ChiA-intermediate 1:' 240 WRITE (LUPRI,'(5G15.6)') (WORK(KCHIA+I-1),I=1,NMATIJ(ISYINT)) 241 END IF 242 CALL DAXPY(NMATIJ(ISYINT),-ONE,XAINT,1,WORK(KCHIA),1) 243 IF (LOCDBG) THEN 244 WRITE (LUPRI,'(//A)') 'CCFBINT1> ChiA-intermediate:' 245 WRITE (LUPRI,'(5G15.6)') (WORK(KCHIA+I-1),I=1,NMATIJ(ISYINT)) 246 END IF 247 END IF 248 249*---------------------------------------------------------------------* 250* calculate the effective BZA density (only for new TA or ZETA) 251* D_{k\a}^{ij} (backtr. with Lambda^p) (does not depend on Lambda-bar) 252* D_{k\a}^{ij} -> BZDENA (in memory) 253*---------------------------------------------------------------------* 254 IF (CCSD .AND. LRELAX) THEN 255 IF ((IDLSTL .NE. ILLSTOLD).OR.(IDLSTR .NE. IRLSTOLD)) THEN 256 257 IF (LOCDBG) THEN 258 WRITE (LUPRI,*) 'CCFBINT1: Print LambdaP for symmetry 1' 259 CALL CC_PRLAM(XLAMDP,XLAMDH,ISYLAM) 260 WRITE (LUPRI,*) 'FBINT1: Print ZETA2 for symmetry 1' 261 CALL CC_PRSQ(ZETA1,ZETA2,ISYCTR,0,1) 262 CALL DZERO(BZDENA,NT2AOIJ(ISYINT)) 263 CALL DZERO(WORK(KEND1),LWRK1) 264 END IF 265 266 CALL CC_BFDENF( ZETA2, ISYCTR, MAINT, ISYINT, 267 * XLAMDP, ISYLAM, WORK(KCHIA), ISYINT, 268 * TA1AMP, ISYMTA, BZDENA, WORK(KEND1), LWRK1) 269 IF (LOCDBG) THEN 270 WRITE (LUPRI,*) 'FBINT1: Print BZDENA for symmetry 1' 271 CALL OUTPUT(BZDENA,1,NT1AO(1),1,NMATIJ(1), 272 * NT1AO(1),NMATIJ(1),1,LUPRI) 273 WRITE (LUPRI,*) 'returned 1 from CC_BFDENF... ' 274 CALL FLSHFO(LUPRI) 275 END IF 276 ENDIF 277 END IF 278 279*---------------------------------------------------------------------* 280* calculate the effective BZAB density which does depend on the 281* external perturbation for _each_ transformation (ITRAN, as it 282* depends on IOPER): (or pass IOPER?) 283* 2) Dbar_{k\a}^{ij} -> BZDENAB (send in LambdapQ zero) 284*---------------------------------------------------------------------* 285 IF (CCSD .AND. LRELAX) THEN 286 CALL CC_BFDENF( ZETA2, ISYCTR, MAINT, ISYINT, 287 * XLAMDPQ, ISYHOP, WORK(KCHIA), ISYINT, 288 * TA1AMP, ISYMTA, BZDENAB, WORK(KEND1), LWRK1) 289 IF (LOCDBG) THEN 290 WRITE (LUPRI,*) 'CCFBINT1: Print BZDENQA for symmetry 1' 291 CALL OUTPUT(BZDENAB,1,NT1AO(1),1,NMATIJ(1), 292 * NT1AO(1),NMATIJ(1),1,LUPRI) 293 WRITE (LUPRI,*) 'returned 2 from CC_BFDENF... ' 294 CALL FLSHFO(LUPRI) 295 END IF 296 ENDIF 297*---------------------------------------------------------------------* 298* calculate one-index backtransformed PA and QA intermediates used 299* in the F and G terms contribution to the result vector: 300*---------------------------------------------------------------------* 301 IF (CCSD .AND. LRELAX) THEN 302 303 IF (( IDLSTL .NE. ILLSTOLD ).OR.( IDLSTR .NE. IRLSTOLD )) THEN 304 305 IOPT = 1 306 CALL CC_PQIMO(ZETA2,ISYCTR,TA2AMP,ISYMTA,YAINT,IOPT, 307 * FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),IADRPQA, 308 * ITRAN,WORK(KEND1),LWRK1) 309 310 CALL CC_PQIAO(FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),ISYINT, 311 * FILPQA0, LUPQA0, IADRPQA0(1,ITRAN), IADRPQAI0, 312 * ITRAN, XLAMDH, ISYLAM,WORK(KEND1), LWRK1) 313 314 315 ELSE 316 DO IDEL = 1, NBAST 317 !dimensioned bigger that Nvir 318 IADRPQAMO(IDEL,ITRAN) = IADRPQAMO(IDEL,ITRAN-1) 319 IADRPQA0(IDEL,ITRAN) = IADRPQA0(IDEL,ITRAN-1) 320 END DO 321 END IF 322 323C Backtransformation with Lambda-bar always performed 324c this is not necessary. I could pass LNEWOP!!!!!!!!!!!!! 325 326 CALL CC_PQIAO(FILPQAMO,LUPQAMO,IADRPQAMO(1,ITRAN),ISYINT, 327 * FILPQA1, LUPQA1, IADRPQA1(1,ITRAN), IADRPQAI1, 328 * ITRAN,XLAMDHQ,ISYHOP,WORK(KEND1), LWRK1) 329 330 END IF 331*---------------------------------------------------------------------* 332* save present IDLSTL in IDLSTOLD and return: 333*---------------------------------------------------------------------* 334 ILLSTOLD = IDLSTL 335 IRLSTOLD = IDLSTR 336 337 RETURN 338 339 END 340*=====================================================================* 341