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 CC_PQI(CTR2,ISYCTR,T2AM,ISYTAMP, 21 & FILNAM,LUPQIM,IADRPQ,IADR,IV,WORK,LWORK) 22*---------------------------------------------------------------------- 23* 24* Purpose: Calculate the P and Q intermediates from the 25* Lagrangian multipiers CTR2 26* and the amplitude vector T2AM 27* and write them to a file direct access file FILNAM 28* 29* P{aik,del} = sum_{dl} (2 t(dl,k;del) - t(dk,l;del)) Zeta(dl;ai) 30* P{aik,del} = sum_{dl} t(dl,k;del) (2Zeta(di;al) + Zeta(dl;ai)) 31* 32* Christof Haettig & Asger Halkier, August 1998 33* 34*====================================================================== 35#if defined (IMPLICIT_NONE) 36 IMPLICIT NONE 37#else 38# include "implicit.h" 39#endif 40#include "priunit.h" 41#include "ccorb.h" 42#include "maxorb.h" 43#include "ccsdsym.h" 44#include "ccsdinp.h" 45#include "iratdef.h" 46#include "second.h" 47 48 LOGICAL LOCDBG 49 PARAMETER (LOCDBG = .FALSE.) 50 INTEGER LWORK, LUPQIM 51 CHARACTER*(*) FILNAM 52 53#if defined (SYS_CRAY) 54 REAL ONE, TWO, THREE, HALF, XNORM, DDOT, ZERO 55 REAL WORK(*), CTR2(*), T2AM(*) 56 REAL TIMET, TIMIO, TIMZWV, TIMTAM, DTIME, TIMTI 57 REAL TIMZET,TIMSCL 58#else 59 DOUBLE PRECISION ONE, TWO, THREE, HALF, DDOT, ZERO 60 DOUBLE PRECISION WORK(*), CTR2(*), T2AM(*) 61 DOUBLE PRECISION TIMET,TIMIO,TIMTI,TIMZWV,TIMTAM,DTIME 62 DOUBLE PRECISION TIMZET,TIMSCL 63#endif 64 PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,HALF=0.5D0,ZERO=0.0D0) 65 66 INTEGER IV, IADR, IADRQ 67 INTEGER IADRPQ(MXCORB_CC,IV) 68 INTEGER ISYTAMP, ISYCTR, ISYMD, ISYTIN, ISYAIK 69 INTEGER KEND1, LWRK1, KEND2, LWRK2 70 INTEGER IOPT, IDEL, ILLL 71 INTEGER KTINT, KTJNT, KPINT, KQINT, KLAMDH, KLAMDP, KT1AM 72 73 74 CALL QENTER('CC_PQI') 75*---------------------------------------------------------------------- 76* set symmetries and allocate work space: 77*---------------------------------------------------------------------- 78 KLAMDH = 1 79 KLAMDP = KLAMDH + NLAMDT 80 KT1AM = KLAMDP + NLAMDT 81 KEND1 = KT1AM + NT1AMX 82 LWRK1 = LWORK - KEND1 83 84 IF ( LWRK1 .LT. 0 ) THEN 85 WRITE (LUPRI,*) 'Insufficient memory in CC_PQI.' 86 CALL QUIT('Insufficient memory in CC_PQI.') 87 END IF 88 89 TIMET = SECOND() 90 TIMIO = ZERO 91 TIMTI = ZERO 92 TIMZWV = ZERO 93 TIMTAM = ZERO 94 TIMZET = ZERO 95 TIMSCL = ZERO 96 97 IF (LOCDBG) THEN 98 WRITE (LUPRI,*) 'Norm of T2AMP:', 99 & DDOT(NT2AM(ISYTAMP),T2AM,1,T2AM,1) 100 WRITE (LUPRI,*) 'T2AMP:' 101 CALL CC_PRP(WORK,T2AM,ISYTAMP,0,1) 102 WRITE (LUPRI,*) 'Norm of CTR2 :', 103 & DDOT(NT2SQ(ISYCTR),CTR2,1,CTR2,1) 104 WRITE (LUPRI,*) 'CTR2:' 105 CALL CC_PRSQ(WORK,CTR2,ISYCTR,0,1) 106 END IF 107*---------------------------------------------------------------------- 108* get XLAMD matrices from zero T1 amplitudes: 109*---------------------------------------------------------------------- 110 CALL DZERO(WORK(KT1AM),NT1AMX) 111 112 CALL LAMMAT(WORK(KLAMDP),WORK(KLAMDH),WORK(KT1AM), 113 & WORK(KEND1),LWRK1) 114 115*---------------------------------------------------------------------- 116* calculate the P intermediate in a loop over AO index delta: 117*---------------------------------------------------------------------- 118 DO ISYMD = 1, NSYM 119 DO ILLL = 1, NBAS(ISYMD) 120 IDEL = IBAS(ISYMD) + ILLL 121 122 ISYTIN = MULD2H(ISYMD,ISYTAMP) 123 ISYAIK = MULD2H(ISYTIN,ISYCTR) 124 125 KTINT = KEND1 126 KTJNT = KTINT + NT2BCD(ISYTIN) 127 KPINT = KTJNT + NT2BCD(ISYTIN) 128 KQINT = KPINT + NT2BCD(ISYAIK) 129 KEND2 = KQINT + NT2BCD(ISYAIK) 130 LWRK2 = LWORK - KEND2 131 132 IF (LWRK2 .LT. 0) THEN 133 WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.' 134 CALL QUIT('Insufficient memory in CC_PQIM.') 135 END IF 136C 137C calculate delta batch of backtransformed amplitudes: 138C 139 DTIME = SECOND() 140 CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,WORK(KLAMDH),1, 141 & WORK(KEND2),LWRK2,IDEL,ISYMD) 142 TIMTI = TIMTI + SECOND() - DTIME 143C 144C calculate 2 x Coul - Exc combination of backtransformed T: 145C 146 DTIME = SECOND() 147 CALL DCOPY(NT2BCD(ISYTIN),WORK(KTINT),1,WORK(KTJNT),1) 148 149 CALL DSCAL(NT2BCD(ISYTIN),-ONE,WORK(KTJNT),1) 150 151 CALL CCLT_P21I(WORK(KTJNT),ISYTIN,WORK(KEND2),LWRK2, 152 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 153 154 CALL DAXPY(NT2BCD(ISYTIN),TWO,WORK(KTINT),1,WORK(KTJNT),1) 155 TIMTAM = TIMTAM + SECOND() - DTIME 156C 157C calculate the P intermediate: 158C 159 IOPT = 1 160 DTIME = SECOND() 161 CALL CC_ZWVI(WORK(KPINT),CTR2,ISYCTR,WORK(KTJNT), 162 & ISYTIN,WORK(KEND2),LWRK2,IOPT) 163 TIMZWV = TIMZWV + SECOND() - DTIME 164 165 DTIME = SECOND() 166 CALL DSCAL(NT2BCD(ISYAIK),HALF,WORK(KPINT),1) 167 TIMSCL = TIMSCL + SECOND() - DTIME 168C 169C write the intermediate to file: 170C 171 IADRPQ(IDEL,IV) = IADR 172 173 DTIME = SECOND() 174 CALL PUTWA2(LUPQIM,FILNAM,WORK(KPINT),IADR,NT2BCD(ISYAIK)) 175 TIMIO = TIMIO + SECOND() - DTIME 176 177 IADR = IADR + 2*NT2BCD(ISYAIK) 178 179 IF (LOCDBG) THEN 180 WRITE (LUPRI,*) 'CC_PQI> P interm. in AO for IDEL = ',IDEL 181 WRITE (LUPRI,'(5(F12.8))') 182 & (WORK(KPINT+I),I=0,NT2BCD(ISYAIK)-1) 183 END IF 184 185 END DO 186 END DO 187 188*---------------------------------------------------------------------- 189* calculate (2 x Exc + Coul)/3 combination of ZETA: 190* (we interrupt here the loop over AO to calculate the modified 191* ZETA and accept that we recalculate the backtransformed 192* amplitude since it the transformation of the amplitudes is 193* less expansive than the transformation and restruction of ZETA 194* for each delta. both the transformation of TAMP and the 195* transformation/restruction of ZETA scale with N V^2 O^2. ) 196*---------------------------------------------------------------------- 197 DTIME = SECOND() 198 CALL CCRHS_T2BT(CTR2,WORK(KEND1),LWRK1,ISYCTR) 199 CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYCTR) 200 TIMZET = TIMZET + SECOND() - DTIME 201 202*---------------------------------------------------------------------- 203* calculate the Q intermediate in a loop over AO index delta: 204*---------------------------------------------------------------------- 205 DO ISYMD = 1, NSYM 206 DO ILLL = 1, NBAS(ISYMD) 207 IDEL = IBAS(ISYMD) + ILLL 208 209 ISYTIN = MULD2H(ISYMD,ISYTAMP) 210 ISYAIK = MULD2H(ISYTIN,ISYCTR) 211 212 KTINT = KEND1 213 KTJNT = KTINT + NT2BCD(ISYTIN) 214 KPINT = KTJNT + NT2BCD(ISYTIN) 215 KQINT = KPINT + NT2BCD(ISYAIK) 216 KEND2 = KQINT + NT2BCD(ISYAIK) 217 LWRK2 = LWORK - KEND2 218 219 IF (LWRK2 .LT. 0) THEN 220 WRITE (LUPRI,*) 'Insufficient memory in CC_PQIM.' 221 CALL QUIT('Insufficient memory in CC_PQIM.') 222 END IF 223C 224C calculate delta batch of backtransformed amplitudes: 225C 226 DTIME = SECOND() 227 CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,WORK(KLAMDH),1, 228 & WORK(KEND2),LWRK2,IDEL,ISYMD) 229 TIMTI = TIMTI + SECOND() - DTIME 230C 231C calculate Q intermediate: 232C 233 IOPT = 2 234 DTIME = SECOND() 235 CALL CC_ZWVI(WORK(KQINT),CTR2,ISYCTR,WORK(KTINT), 236 & ISYTIN,WORK(KEND2),LWRK2,IOPT) 237 TIMZWV = TIMZWV + SECOND() - DTIME 238 239 DTIME = SECOND() 240 CALL DSCAL(NT2BCD(ISYAIK),THREE*HALF,WORK(KQINT),1) 241 TIMSCL = TIMSCL + SECOND() - DTIME 242C 243C write the intermediate to file: 244C 245 IADRQ = IADRPQ(IDEL,IV) + NT2BCD(ISYAIK) 246 247 DTIME = SECOND() 248 CALL PUTWA2(LUPQIM,FILNAM,WORK(KQINT),IADRQ,NT2BCD(ISYAIK)) 249 TIMIO = TIMIO + SECOND() - DTIME 250 251 IF (LOCDBG) THEN 252 WRITE (LUPRI,*) 'CC_PQI> Q interm. in AO for IDEL = ',IDEL 253 WRITE (LUPRI,'(5(F12.8))') 254 & (WORK(KQINT+I),I=0,NT2BCD(ISYAIK)-1) 255 END IF 256 257 END DO 258 END DO 259 260*---------------------------------------------------------------------- 261* recover Zeta vector: 262*---------------------------------------------------------------------- 263 DTIME = SECOND() 264 CALL CCSD_T2TP(CTR2,WORK(KEND1),LWRK1,ISYCTR) 265 CALL CCRHS_T2TR(CTR2,WORK(KEND1),LWRK1,ISYCTR) 266 TIMZET = TIMZET + SECOND() - DTIME 267 268*---------------------------------------------------------------------- 269* return: 270*---------------------------------------------------------------------- 271 IF (IPRINT.GE.10) THEN 272 TIMET = SECOND() - TIMET 273 WRITE(LUPRI,*) ' Timings of CC_PQIM:' 274 WRITE(LUPRI,1) 'I/O ', TIMIO 275 WRITE(LUPRI,1) 'CC_TI ', TIMTI 276 WRITE(LUPRI,1) '2C-E of T2 ', TIMTAM 277 WRITE(LUPRI,1) '2E+C of L2 ', TIMZET 278 WRITE(LUPRI,1) 'CC_ZWV ', TIMZWV 279 WRITE(LUPRI,1) 'scaling etc. ', TIMSCL 280 WRITE(LUPRI,1) 'total CC_PQIM: ', TIMET 281 END IF 282 283 1 FORMAT(1x,'Time used for',2x,A18,2x,': ',f10.2,' seconds') 284 285 CALL QEXIT('CC_PQI') 286 RETURN 287 END 288C 289*====================================================================== 290 SUBROUTINE CC_PQI3(CTR2,ISYCTR,T2AM,ISYTAMP, 291 & FILNAM,LUPQIM,IADRPQ,IADR,IV, 292 & XLAMDH,ICAL,WORK,LWORK) 293*---------------------------------------------------------------------- 294* 295* Purpose: Calculate the P intermediates from the 296* Lagrangian multipiers CTR2 297* and the amplitude vector T2AM 298* and write them to a file direct access file FILNAM 299* 300* P{aik,del} = sum_{dl} (2 t(dl,k;del) - t(dk,l;del)) Ctr2(dl;ai) 301C where in the triplet case, ctr2 should contain 302C (+)L(dl;ai) - (-)L(dl;ai) 303* 304C Based on CCPQI31 by 305* Christof Haettig & Asger Halkier, August 1998 306* 307*====================================================================== 308 IMPLICIT NONE 309#include "priunit.h" 310#include "ccorb.h" 311#include "maxorb.h" 312#include "ccsdsym.h" 313#include "ccsdinp.h" 314#include "iratdef.h" 315#include "second.h" 316 317 CHARACTER, PARAMETER :: MYNAME = 'CC_PQI3' 318 INTEGER, PARAMETER :: NCAL = 3 ! Number of times this is called 319 LOGICAL LOCDBG 320 PARAMETER (LOCDBG = .FALSE.) 321 322 323 INTEGER, INTENT(IN) :: LWORK, LUPQIM 324 INTEGER, INTENT(INOUT) :: IADRPQ(MXCORB_CC,IV), IADR 325 INTEGER, INTENT(IN) :: ISYCTR, ISYTAMP, IV, ICAL 326 CHARACTER*(*),INTENT(IN) :: FILNAM 327 328 329 DOUBLE PRECISION, INTENT(IN) :: WORK(*), CTR2(*), T2AM(*), 330 & XLAMDH(*) 331 DOUBLE PRECISION TIMET,TIMIO,TIMTI,TIMZWV,TIMTAM,DTIME 332 DOUBLE PRECISION TIMZET,TIMSCL 333 334 DOUBLE PRECISION ONE, TWO, THREE, HALF, DDOT, ZERO 335 PARAMETER(ONE=1.0D0,TWO=2.0D0,THREE=3.0D0,HALF=0.5D0,ZERO=0.0D0) 336 337 INTEGER :: IADRQ 338 INTEGER ISYMD, ISYTIN, ISYAIK 339 INTEGER KEND1, LWRK1, KEND2, LWRK2 340 INTEGER IOPT, IDEL, ILLL 341 INTEGER KTINT, KTJNT, KPINT, KQINT, KLAMDH, KLAMDP, KT1AM 342 343 344 CALL QENTER(myname) 345*---------------------------------------------------------------------- 346* set symmetries and allocate work space: 347*---------------------------------------------------------------------- 348 KEND1 = 1 349 LWRK1 = LWORK 350 351 IF ( LWRK1 .LT. 0 ) THEN 352 WRITE (LUPRI,*) 'Insufficient memory in '//myname//'.' 353 CALL QUIT('Insufficient memory in '//myname//'.') 354 END IF 355 356 TIMET = SECOND() 357 TIMIO = ZERO 358 TIMTI = ZERO 359 TIMZWV = ZERO 360 TIMTAM = ZERO 361 TIMZET = ZERO 362 TIMSCL = ZERO 363 364 IF (LOCDBG) THEN 365 WRITE (LUPRI,*) 'Norm of T2AMP:', 366 & DDOT(NT2AM(ISYTAMP),T2AM,1,T2AM,1) 367 WRITE (LUPRI,*) 'T2AMP:' 368 CALL CC_PRP(WORK,T2AM,ISYTAMP,0,1) 369 WRITE (LUPRI,*) 'Norm of CTR2 :', 370 & DDOT(NT2SQ(ISYCTR),CTR2,1,CTR2,1) 371 WRITE (LUPRI,*) 'CTR2:' 372 CALL CC_PRSQ(WORK,CTR2,ISYCTR,0,1) 373 END IF 374 375*---------------------------------------------------------------------- 376* calculate the P intermediate in a loop over AO index delta: 377*---------------------------------------------------------------------- 378 DO ISYMD = 1, NSYM 379 DO ILLL = 1, NBAS(ISYMD) 380 IDEL = IBAS(ISYMD) + ILLL 381 382 ISYTIN = MULD2H(ISYMD,ISYTAMP) 383 ISYAIK = MULD2H(ISYTIN,ISYCTR) 384 385 KTINT = KEND1 386 KTJNT = KTINT + NT2BCD(ISYTIN) 387 KPINT = KTJNT + NT2BCD(ISYTIN) 388 KEND2 = KPINT + NT2BCD(ISYAIK) 389 LWRK2 = LWORK - KEND2 390 391 IF (LWRK2 .LT. 0) THEN 392 WRITE (LUPRI,*) 'Insufficient memory in '//myname//'.' 393 CALL QUIT('Insufficient memory in '//myname//'.') 394 END IF 395C 396C calculate delta batch of backtransformed amplitudes: 397C 398 DTIME = SECOND() 399 CALL CC_TI(WORK(KTINT),ISYTIN,T2AM,ISYTAMP,XLAMDH,1, 400 & WORK(KEND2),LWRK2,IDEL,ISYMD) 401 TIMTI = TIMTI + SECOND() - DTIME 402C 403C calculate 2 x Coul - Exc combination of backtransformed T: 404C 405 IF (ICAL .EQ. 0) THEN 406 DTIME = SECOND() 407 CALL DCOPY(NT2BCD(ISYTIN),WORK(KTINT),1,WORK(KTJNT),1) 408 409 CALL DSCAL(NT2BCD(ISYTIN),TWO,WORK(KTINT),1) 410 411 CALL CCLT_P21I(WORK(KTJNT),ISYTIN,WORK(KEND2),LWRK2, 412 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 413 414 CALL DAXPY(NT2BCD(ISYTIN),-ONE,WORK(KTJNT),1,WORK(KTINT),1) 415 416 IADRPQ(IDEL,IV) = IADR 417 IADR = IADR + NCAL*NT2BCD(ISYAIK) 418 ELSE 419 CALL CCLT_P21I(WORK(KTINT),ISYTIN,WORK(KEND2),LWRK2, 420 & IT2BCD,NT2BCD,IT1AM,NT1AM,NVIR) 421 END IF 422 IADRQ = IADRPQ(IDEL,IV) + NT2BCD(ISYAIK)*ICAL 423 424 TIMTAM = TIMTAM + SECOND() - DTIME 425C 426C calculate the P intermediate: 427C 428 IOPT = 1 429 DTIME = SECOND() 430 CALL CC_ZWVI3(WORK(KPINT),CTR2,ISYCTR,WORK(KTINT), 431 & ISYTIN,WORK(KEND2),LWRK2) 432 TIMZWV = TIMZWV + SECOND() - DTIME 433 434 DTIME = SECOND() 435 CALL DSCAL(NT2BCD(ISYAIK),HALF,WORK(KPINT),1) 436 TIMSCL = TIMSCL + SECOND() - DTIME 437C 438C write the intermediate to file: 439C 440 441 DTIME = SECOND() 442 CALL PUTWA2(LUPQIM,FILNAM,WORK(KPINT),IADRQ,NT2BCD(ISYAIK)) 443 TIMIO = TIMIO + SECOND() - DTIME 444 445 446 IF (LOCDBG) THEN 447 WRITE (LUPRI,*) 'CC_PQI> P interm. in AO for IDEL = ',IDEL 448 WRITE (LUPRI,'(5(F12.8))') 449 & (WORK(KPINT+I),I=0,NT2BCD(ISYAIK)-1) 450 END IF 451 452 END DO 453 END DO 454 455 CALL QEXIT(myname) 456 RETURN 457 END 458 459