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! 18 subroutine ccsdt_t32_noddy(tab3am,listab,idlstab,freqab, 19 & xint1s0,xint2s0, 20 & xlamdp0,xlamdh0,fock0,fockd, 21 & fieldao,field, 22 & scr1,work,lwork) 23*---------------------------------------------------------------------* 24* Purpose: compute triples part of second-order response amplitudes 25* 26* Input: listab, idlstab 27* xlamdp0, xlamdh0, fock0, fockd 28* xint1s, xint2s 29* work, lwork 30* 31* Output: tab3am, freqab 32* 33* Written by Christof Haettig, Mai 2003, based on ccsdt_t31_noddy. 34*---------------------------------------------------------------------* 35 implicit none 36#include "priunit.h" 37#include "ccsdinp.h" 38#include "maxorb.h" 39#include "ccsdsym.h" 40#include "ccorb.h" 41#include "cco2rsp.h" 42#include "ccr2rsp.h" 43#include "ccer1rsp.h" 44#include "ccexci.h" 45#include "ccfield.h" 46#include "dummy.h" 47 48 49 logical locdbg, print_t3 50 parameter (locdbg=.false., print_t3=.false.) 51 52 integer isym0 53 parameter ( isym0 = 1 ) 54 55 character*3 listab 56 integer idlstab, lwork 57 58#if defined (SYS_CRAY) 59 real tab3am(*), freqab, freqa, freqb 60 real scr1(*), work(*), xlamdp0(*), xlamdh0(*) 61 real fock0(*), fockd(*), xint1s0(*), xint2s0(*) 62 real field(*), fieldao(*) 63 real two, one, ddot 64#else 65 double precision tab3am(*), freqab, freqa, freqb 66 double precision scr1(*), work(*), xlamdp0(*), xlamdh0(*) 67 double precision fock0(*), fockd(*), xint1s0(*), xint2s0(*) 68 double precision field(*), fieldao(*) 69 double precision two, one, ddot 70#endif 71 parameter ( one = 1.0d0, two = 2.0d0 ) 72 73 character lista*3, listb*3, labela*8, labelb*8, model*10 74 logical lorxa, lorxb, rhs_only 75 integer isyma, isymb, idlsta, idlstb, ierr, irrep, kend1, lwrk1, 76 & kta1am, kta2am, klampa, klamha, kfocka_ao, kfocka, 77 & ktb1am, ktb2am, klampb, klamhb, kfockb_ao, kfockb, 78 & kint1sa, kint2sa, kint1sb, kint2sb, kint1sab, kint2sab, 79 & kta3am, ktb3am, kt03am, iopt, ktab1am, ktab2am, isymab, 80 & kt02am, kfockab, kfckbuf 81 82* external functions: 83 integer ir1tamp, ilstsym 84 85 call qenter('CCT32NOD') 86*---------------------------------------------------------------------* 87* Do some initializations and checks: 88*---------------------------------------------------------------------* 89 if (listab(1:3).eq.'O2 ') then 90 ! we compute the rhs vectors for second-order amplitude resp. 91 rhs_only = .true. 92 93 lista = 'R1 ' 94 labela = lblo2(idlstab,1) 95 isyma = isyo2(idlstab,1) 96 freqa = frqo2(idlstab,1) 97 lorxa = lorxo2(idlstab,1) 98 idlsta = ir1tamp(labela,lorxa,freqa,isyma) 99 100 listb = 'R1 ' 101 labelb = lblo2(idlstab,2) 102 isymb = isyo2(idlstab,2) 103 freqb = frqo2(idlstab,2) 104 lorxb = lorxo2(idlstab,2) 105 idlstb = ir1tamp(labelb,lorxb,freqb,isymb) 106 107 else if (listab(1:3).eq.'R2 ') then 108 ! we compute the second-order amplitude response vectors 109 rhs_only = .false. 110 111 lista = 'R1 ' 112 labela = lblr2t(idlstab,1) 113 isyma = isyr2t(idlstab,1) 114 freqa = frqr2t(idlstab,1) 115 lorxa = lorxr2t(idlstab,1) 116 idlsta = ir1tamp(labela,lorxa,freqa,isyma) 117 118 listb = 'R1 ' 119 labelb = lblr2t(idlstab,2) 120 isymb = isyr2t(idlstab,2) 121 freqb = frqr2t(idlstab,2) 122 lorxb = lorxr2t(idlstab,2) 123 idlstb = ir1tamp(labelb,lorxb,freqb,isymb) 124 125 else if (listab(1:3).eq.'EO1' .or. listab(1:3).eq.'ER1') then 126 ! ER1 - first-order response of excited states 127 ! EO1 - rhs vectors for ER1 equations 128 rhs_only = listab(1:3) .eq. 'EO1' 129 130 lista = 'R1 ' 131 labela = lbler1(idlstab) 132 isyma = isyoer1(idlstab) 133 freqa = frqer1(idlstab) 134 lorxa = lorxer1(idlstab) 135 idlsta = ir1tamp(labela,lorxa,freqa,isyma) 136 137 listb = 'RE ' 138 labelb = '-- XX --' 139 isymb = isyser1(idlstab) 140 freqb = eiger1(idlstab) 141 lorxb = .false. 142 idlstb = ister1(idlstab) 143 144 else 145 call quit('Unknown LISTAB in CCSDT_T32_NODDY.') 146 end if 147 148 isymab = ilstsym(listab,idlstab) 149 freqab = freqa + freqb 150 151 152 if (lorxa.or.lorxb) then 153 call quit('Orbital relaxation not allowed in CCSDT_T32_NODDY.') 154 end if 155 156 if (listb(1:3).ne.'R1 ' .and. listb(1:3).NE.'RE ') then 157 call quit('Unknown LISTB in CCSDT_T32_NODDY.') 158 end if 159 160 if (lista(1:3).ne.'R1 ' .and. lista(1:3).NE.'RE ') then 161 call quit('Unknown LISTA in CCSDT_T32_NODDY.') 162 end if 163 164 call dzero(tab3am,nt1amx*nt1amx*nt1amx) 165 166*---------------------------------------------------------------------* 167* Memory allocation: 168*---------------------------------------------------------------------* 169 kend1 = 1 170 171 kta1am = kend1 172 klampa = kta1am + nt1amx 173 klamha = klampa + nlamdt 174 kend1 = klamha + nlamdt 175 176 ktb1am = kend1 177 klampb = ktb1am + nt1amx 178 klamhb = klampb + nlamdt 179 kend1 = klamhb + nlamdt 180 181 kt02am = kend1 182 kta2am = kt02am + nt1amx*nt1amx 183 ktb2am = kta2am + nt1amx*nt1amx 184 kend1 = ktb2am + nt1amx*nt1amx 185 186 if (lista.eq.'R1 ') then 187 kfocka_ao = kend1 188 kfocka = kfocka_ao + nbast*nbast 189 kend1 = kfocka + nbast*nbast 190 else 191 kfocka_ao = -999 999 192 kfocka = -999 999 193 end if 194 195 if (listb.eq.'R1 ') then 196 kfockb_ao = kend1 197 kfockb = kfockb_ao + nbast*nbast 198 kend1 = kfockb + nbast*nbast 199 else 200 kfockb_ao = -999 999 201 kfockb = -999 999 202 end if 203 204 if (lista.eq.'R1 '.or.listb.eq.'R1 ') then 205 kfockab = kend1 206 kfckbuf = kfockab + nbast*nbast 207 kend1 = kfckbuf + nbast*nbast 208 else 209 kfockab = -999 999 210 kfckbuf = -999 999 211 end if 212 213 kint1sa = kend1 214 kint2sa = kint1sa + nt1amx*nvirt*nvirt 215 kend1 = kint2sa + nrhft*nrhft*nt1amx 216 217 kint1sb = kend1 218 kint2sb = kint1sb + nt1amx*nvirt*nvirt 219 kend1 = kint2sb + nrhft*nrhft*nt1amx 220 221 kint1sab = kend1 222 kint2sab = kint1sab + nt1amx*nvirt*nvirt 223 kend1 = kint2sab + nrhft*nrhft*nt1amx 224 225 kta3am = kend1 226 kend1 = kta3am + nt1amx*nt1amx*nt1amx 227 ktb3am = kta3am 228 kt03am = kta3am 229 230 lwrk1 = lwork - kend1 231 if (lwrk1 .lt. 0) then 232 call quit('Insufficient space in CCSDT_T32_NODDY') 233 endif 234 235*---------------------------------------------------------------------* 236* Read zeroth-order amplitudes and response vectors into memory 237*---------------------------------------------------------------------* 238 if (lwrk1.lt.nt2amx) then 239 call quit('Insufficient space in CCSDT_T32_NODDY.') 240 end if 241 242 iopt = 2 243 call cc_rdrsp('R0 ',0,isym0,iopt,model,dummy,work(kend1)) 244 call cc_t2sq(work(kend1),work(kt02am),isym0) 245 246 iopt = 3 247 call cc_rdrsp(lista,idlsta,isyma,iopt,model, 248 & work(kta1am),work(kend1)) 249 call cclr_diascl(work(kend1),two,isyma) 250 call cc_t2sq(work(kend1),work(kta2am),isyma) 251 252 iopt = 3 253 call cc_rdrsp(listb,idlstb,isymb,iopt,model, 254 & work(ktb1am),work(kend1)) 255 call cclr_diascl(work(kend1),two,isymb) 256 call cc_t2sq(work(kend1),work(ktb2am),isymb) 257 258*---------------------------------------------------------------------* 259* Get property matrices A and B and the matrix A^B+B^A: 260*---------------------------------------------------------------------* 261 if (lista(1:3).eq.'R1 ') then 262 ! read property integrals from file: 263 call ccprpao(labela,.TRUE.,work(kfocka_ao),irrep,isyma,ierr, 264 & work(kend1),lwrk1) 265 if ((ierr.gt.0) .or. (ierr.eq.0 .and. irrep.ne.isyma)) then 266 call quit('CCSDT_T32_NODDY: error reading operator '//LABELA) 267 else if (ierr.lt.0) then 268 call dzero(work(kfocka_ao),n2bst(isyma)) 269 end if 270 call dcopy(nbast*nbast,work(kfocka_ao),1,work(kfocka),1) 271 272 ! transform property integrals to Lambda-MO basis 273 call cc_fckmo(work(kfocka),xlamdp0,xlamdh0, 274 & work(kend1),lwrk1,isyma,1,1) 275 end if 276 277 if (listb(1:3).eq.'R1 ') then 278 ! read property integrals from file: 279 call ccprpao(labelb,.TRUE.,work(kfockb_ao),irrep,isymb,ierr, 280 & work(kend1),lwrk1) 281 if ((ierr.gt.0) .or. (ierr.eq.0 .and. irrep.ne.isymb)) then 282 call quit('CCSDT_T32_NODDY: error reading operator '//LABELB) 283 else if (ierr.lt.0) then 284 call dzero(work(kfockb_ao),n2bst(isymb)) 285 end if 286 call dcopy(nbast*nbast,work(kfockb_ao),1,work(kfockb),1) 287 288 ! transform property integrals to Lambda-MO basis 289 call cc_fckmo(work(kfockb),xlamdp0,xlamdh0, 290 & work(kend1),lwrk1,isymb,1,1) 291 end if 292 293 294*---------------------------------------------------------------------* 295* Compute contributions 296* 297* <mu_3| [B,T^A_3] + [[B,T^A_2],T^0_2] |HF> 298* 299* Compute corrections to triples vector T^A_3, and corresponding 300* lambda matrices and the XINT1SA,XINT2SA integrals and set FREQA, 301*---------------------------------------------------------------------* 302 if (listb(1:3).eq.'R1 ') then 303 304 if (lista(1:3).eq.'R1 ' .or. lista(1:3).eq.'RE ' .or. 305 & lista(1:3).eq.'RC ' ) then 306 call ccsdt_t31_noddy(work(kta3am),lista,idlsta,freqa,.false., 307 & .false.,xint1s0,xint2s0, 308 & .false.,dummy,dummy, 309 & .false.,dummy,dummy, 310 & work(kint1sa),work(kint2sa), 311 & work(klampa),work(klamha),work(kfocka), 312 & xlamdp0,xlamdh0,fock0,dummy,fockd, 313 & work(kend1),lwrk1) 314 call dscal(nt1amx*nt1amx*nt1amx,-1.0d0,work(kta3am),1) 315 else 316 call quit('Unknown LISTA in CCSDT_T32_NODDY.') 317 end if 318 319c write(lupri,*) 'norm^2(tab3arm) before cont. of t^B_3:', 320c & ddot(nt1amx**3,tab3am,1,tab3am,1) 321 call ccsdt_xksi3_2(tab3am,work(kfockb),work(kta3am)) 322c write(lupri,*) 'norm^2(tab3arm) after cont. of t^B_3-1:', 323c & ddot(nt1amx**3,tab3am,1,tab3am,1) 324 325 call ccsdt_xksi3_1(tab3am,work(kfockb), 326 & work(kt02am),work(kta2am),one) 327 call ccsdt_xksi3_1(tab3am,work(kfockb), 328 & work(kta2am),work(kt02am),one) 329c write(lupri,*) 'norm^2(tab3arm) after cont. of t^B_3-2:', 330c & ddot(nt1amx**3,tab3am,1,tab3am,1) 331 332 else 333 ! we don't need the triples vector, but still we need in 334 ! the following the response Lambda matrices and the 335 ! one-index transformed integrals 336 337 call cclr_lamtra(xlamdp0,work(klampa),xlamdh0,work(klamha), 338 & work(kta1am),isyma) 339 340 call ccsdt_ints1_noddy(.true.,work(kint1sa),work(kint2sa), 341 & .false.,dummy,dummy, 342 & xlamdp0,xlamdh0, 343 & work(klampa),work(klamha), 344 & work(kend1),lwrk1) 345 346 end if 347 348 if (locdbg) then 349 write(lupri,*) 'norm^2(tab3am) after cont. of t^B_3:', 350 & ddot(nt1amx**3,tab3am,1,tab3am,1) 351 end if 352 353*---------------------------------------------------------------------* 354* Compute the contributions 355* 356* <mu_3| [A,T^B_3] + [[A,T^B_2],T^0_2] |HF> 357* 358* Compute corrections to triples vector T^B_3, and corresponding 359* lambda matrices and the XINT1SB,XINT2SB integrals and set FREQB: 360*---------------------------------------------------------------------* 361 if (lista(1:3).eq.'R1 ') then 362 363 if (listb(1:3).eq.'R1 ' .or. listb(1:3).eq.'RE ' .or. 364 & listb(1:3).eq.'RC ' ) then 365 call ccsdt_t31_noddy(work(ktb3am),listb,idlstb,freqb,.false., 366 & .false.,xint1s0,xint2s0, 367 & .false.,dummy,dummy, 368 & .false.,dummy,dummy, 369 & work(kint1sb),work(kint2sb), 370 & work(klampb),work(klamhb),work(kfockb), 371 & xlamdp0,xlamdh0,fock0,dummy,fockd, 372 & work(kend1),lwrk1) 373 call dscal(nt1amx*nt1amx*nt1amx,-1.0d0,work(ktb3am),1) 374 else 375 call quit('Unknown LISTB in CCSDT_T32_NODDY.') 376 end if 377 378 call ccsdt_xksi3_2(tab3am,work(kfocka),work(ktb3am)) 379 380 call ccsdt_xksi3_1(tab3am,work(kfocka), 381 & work(kt02am),work(ktb2am),one) 382 call ccsdt_xksi3_1(tab3am,work(kfocka), 383 & work(ktb2am),work(kt02am),one) 384 385 else 386 ! we don't need the triples vector, but still we need in 387 ! the following the response Lambda matrices and the 388 ! one-index transformed integrals 389 390 call cclr_lamtra(xlamdp0,work(klampb),xlamdh0,work(klamhb), 391 & work(ktb1am),isymb) 392 393 call ccsdt_ints1_noddy(.true.,work(kint1sb),work(kint2sb), 394 & .false.,dummy,dummy, 395 & xlamdp0,xlamdh0, 396 & work(klampb),work(klamhb), 397 & work(kend1),lwrk1) 398 399 end if 400 401 if (locdbg) then 402 write(lupri,*) 'norm^2(tab3arm) after cont. of t^A_3:', 403 & ddot(nt1amx**3,tab3am,1,tab3am,1) 404 end if 405 406*---------------------------------------------------------------------* 407* Compute contributions: 408* 409* <mu_3| [ ([A,T^B_1] + [B,T^A_1]) , T^0_3 ] |HF> 410* 411* Compute the matrix with one-index transformed property 412* integrals FOCKAB = [A,T^B_1] + [B,T^A_1] 413*---------------------------------------------------------------------* 414 if (lista(1:3).eq.'R1 ' .or. listb(1:3).eq.'R1 ') then 415 call dzero(work(kfockab),nbast*nbast) 416 417 if (lista(1:3).eq.'R1 ') then 418 ! add [A,T^B_1] 419 call dcopy(nbast*nbast,work(kfocka_ao),1,work(kfckbuf),1) 420 call cc_fckmo(work(kfckbuf),work(klampb),xlamdh0, 421 & work(kend1),lwrk1,isyma,isymb,isym0) 422 call daxpy(nbast*nbast,one,work(kfckbuf),1,work(kfockab),1) 423 424 call dcopy(nbast*nbast,work(kfocka_ao),1,work(kfckbuf),1) 425 call cc_fckmo(work(kfckbuf),xlamdp0,work(klamhb), 426 & work(kend1),lwrk1,isyma,isym0,isymb) 427 call daxpy(nbast*nbast,one,work(kfckbuf),1,work(kfockab),1) 428 end if 429 430 if (listb(1:3).eq.'R1 ') then 431 ! add [B,T^A_1] 432 call dcopy(nbast*nbast,work(kfockb_ao),1,work(kfckbuf),1) 433 call cc_fckmo(work(kfckbuf),work(klampa),xlamdh0, 434 & work(kend1),lwrk1,isymb,isyma,isym0) 435 call daxpy(nbast*nbast,one,work(kfckbuf),1,work(kfockab),1) 436 437 call dcopy(nbast*nbast,work(kfockb_ao),1,work(kfckbuf),1) 438 call cc_fckmo(work(kfckbuf),xlamdp0,work(klamha), 439 & work(kend1),lwrk1,isymb,isym0,isyma) 440 call daxpy(nbast*nbast,one,work(kfckbuf),1,work(kfockab),1) 441 end if 442 443 if (nonhf .and. (nfield.gt.0) .and. 444 & (lwrk1.lt.nt1amx*nt1amx*nt1amx)) then 445 call quit('Insufficient space in CCSDT_T32_NODDY') 446 endif 447 448 call dzero(work(kt03am),nt1amx*nt1amx*nt1amx) 449 call ccsdt_t03am(work(kt03am),xint1s0,xint2s0, 450 & work(kt02am),scr1,fockd, 451 & field,work(kend1)) 452 453 call ccsdt_xksi3_2(tab3am,work(kfockab),work(kt03am)) 454 end if 455 456 if (locdbg) then 457 call ccsdt_clean_t3(tab3am,nt1amx,nvirt,nrhft) 458 write(lupri,*) 'norm^2(tab3am) after A{O} contrib.:', 459 & ddot(nt1amx**3,tab3am,1,tab3am,1) 460 call print_pt3_noddy(tab3am) 461 end if 462 463*---------------------------------------------------------------------* 464* Compute double one-index transformed integrals and evaluate 465* the B matrix contributions 466* 467* <mu_3| [H^AB,T^0_2] + [H^A,T^B_2] + [H^B,T^A_2] |HF> 468* 469*---------------------------------------------------------------------* 470 call ccsdt_ints2_noddy(work(kint1sab),work(kint2sab), 471 & xlamdp0,xlamdh0, 472 & work(klampb),work(klamhb), 473 & work(klampa),work(klamha), 474 & work(kend1),lwrk1) 475 476 if (locdbg) then 477 write(lupri,*) 'norm^2(tab3am) before b3am:', 478 & ddot(nt1amx**3,tab3am,1,tab3am,1) 479 end if 480 481 ! note: this routine overwrites work(kt02am) 482 call ccsdt_b3am(tab3am, 483 & work(kint1sab),work(kint2sab),fockd, 484 & lista,idlsta,work(kint1sa),work(kint2sa), 485 & listb,idlstb,work(kint1sb),work(kint2sb), 486 & scr1,work(kt02am),work(kend1),lwrk1) 487 if (nonhf .and. (nfield.gt.0)) 488 & call quit('Finite Field unfinished in CCSDT_T32_NODDY.') 489cch 490c call dscal(nt1amx*nt1amx*nt1amx,-1.0d0,tab3am,1) 491cch 492c call dscal(nt1amx**3,0.5d0,tab3am,1) 493 494 if (locdbg) then 495 call ccsdt_clean_t3(tab3am,nt1amx,nvirt,nrhft) 496 write(lupri,*) 'norm^2(tab3am) after b3am:', 497 & ddot(nt1amx**3,tab3am,1,tab3am,1) 498 end if 499 500*---------------------------------------------------------------------* 501* if solution vector requested add contribution from jacobian: 502* 503* <mu_3| [[H,T^AB_1],T^0_2] + [H,T^AB_2] |HF> 504* 505* and solve the triples equations: 506* 507* 508*---------------------------------------------------------------------* 509 if (.not. rhs_only) then 510 ktab1am = ktb1am 511 ktab2am = ktb2am 512 513 if (lwrk1.lt.max(nt2amx,nt1amx*nt1amx*nt1amx)) then 514 call quit('Insufficient space in CCSDT_T32_NODDY.') 515 end if 516 517 iopt = 2 518 call cc_rdrsp('R0 ',0,isym0,iopt,model,dummy,work(kend1)) 519 call cc_t2sq(work(kend1),work(kt02am),isym0) 520 521 iopt = 3 522 call cc_rdrsp(listab,idlstab,isymab,iopt,model, 523 & work(ktab1am),work(kend1)) 524 call cclr_diascl(work(kend1),two,isymab) 525 call cc_t2sq(work(kend1),work(ktab2am),isymab) 526 527 ! note: here the doubly one-index transformed integrals 528 ! on kint1sab and kint2sab are overwritten by the 529 ! one t^ab one-index transformed integrals 530 call ccsdt_a3am(tab3am,work(ktab1am),work(ktab2am), 531 & freqab,isymab,work(kt02am),work(kt03am), 532 & xint1s0,xint2s0, 533 & work(kint1sab),work(kint2sab), 534 & fockd,xlamdp0,xlamdh0, 535 & fieldao,field,scr1,work(kend1),lwrk1) 536 call ccsdt_clean_t3(tab3am,nt1amx,nvirt,nrhft) 537 538 if (locdbg) then 539 write(lupri,*) 'norm^2(tab3am) after a3am:', 540 & ddot(nt1amx**3,tab3am,1,tab3am,1) 541 end if 542 543 call ccsdt_3am(tab3am,freqab,scr1,fockd, 544 & nonhf,field,.false.,work(kend1)) 545 546 call dscal(nt1amx*nt1amx*nt1amx,-1.0d0,tab3am,1) 547 548 end if 549 550 call ccsdt_clean_t3(tab3am,nt1amx,nvirt,nrhft) 551 552 if (print_t3) then 553 write(lupri,*)'CCSDT_T32_AM> vector type:',listab 554 write(lupri,*)'CCSDT_T32_AM> list,idlst:',listab,idlstab 555 write(lupri,*)'CCSDT_T32_AM> freq:',freqab 556 call print_pt3_noddy(tab3am) 557 end if 558 559 call qexit('CCT32NOD') 560 return 561 end 562*---------------------------------------------------------------------* 563* end of subroutine ccsdt_t32_noddy 564*---------------------------------------------------------------------* 565