1************************************************************************ 2 subroutine gucc_fock(irefspc,itrefspc, 3 & ccvec1,ccvec2,ccvec3,ccvec4, 4 & civec1,civec2,c2vec, 5 & n_cc_typ,i_cc_typ, 6 & namp_cc_typ,ioff_cc_typ, 7 & iopsym,mxb_ci, 8 & luc,luamp,luomg,ludia, 9 & luec,luhc) 10************************************************************************ 11* 12* purpose : driver for tests of fock-space GUCC 13* 14* ak, early 2004 15* 16************************************************************************ 17* 18* units: 19* luc = definition of reference function 20* luamp = amplitude vectors (also output for most recent vector) 21* luampold = scratch containing old vectors from previous iterations 22* (on input it may also be a first trial vector) 23* luomg = error vectors 24* ludia = diagonal preconditioner 25* luec = scratch for exp(G)|ref> 26* luhc = scratch for H exp(G)|ref> 27 28* diverse inludes with commons and paramters 29c include 'implicit.inc' 30c include 'mxpdim.inc' 31 include 'wrkspc.inc' 32 include 'crun.inc' 33 include 'cstate.inc' 34 include 'cgas.inc' 35 include 'ctcc.inc' 36 include 'gasstr.inc' 37 include 'strinp.inc' 38 include 'orbinp.inc' 39 include 'lucinp.inc' 40 include 'cprnt.inc' 41 include 'corbex.inc' 42 include 'csm.inc' 43 include 'cecore.inc' 44 include 'gtbce.inc' 45 include 'cintfo.inc' 46 include 'glbbas.inc' 47 include 'opti.inc' 48* constants 49 integer, parameter :: 50 & ntest = 5 51 52* arrays 53 integer :: 54 & ioff_cc_typ(n_cc_typ), namp_cc_typ(n_cc_typ), 55 & i_cc_typ(4*ngas,n_cc_typ) 56 real*8 :: 57 & ccvec1(n_cc_amp), ccvec2(n_cc_amp), ccvec3(n_cc_amp) 58* local 59 logical :: 60 & calc_Omg, calc_gradE, tstgrad, tst_hss, comm_ops, 61 & do_eag, do_foo, do_hss 62 character*8 cctype 63 integer :: 64 & ictp(n_cc_typ) 65* external functions 66 real*8, external :: inprod, inprdd 67 68* ============================================================ 69* initialize : restart, set coefs to zero 70* ============================================================ 71 72 call atim(cpu0,wall0) 73 74 nprint = 1 75 lblk = -1 76 77 if (ntest.ge.5) then 78 write(6,*) '=======================' 79 write(6,*) 'entered gucc_fock with:' 80 write(6,*) '=======================' 81 write(6,*) ' iopsym = ',iopsym 82 end if 83 84* unit inits 85 lusc1 = iopen_nus('GTBSCR1') 86 lusc2 = iopen_nus('GTBSCR2') 87 lusc3 = iopen_nus('GTBSCR3') 88 lusc4 = iopen_nus('GTBSCR4') 89 lusc5 = iopen_nus('GTBSCR5') 90c lusc6 = iopen_nus('GTBSCR6') 91c lusc7 = iopen_nus('GTBSCR7') 92c lusc8 = iopen_nus('GTBSCR8') 93c lusc9 = iopen_nus('GTBSCR9') 94 95 call memman(idum,idum,'MARK ',idum,'FOCKSP') 96 97! be sure to have the correct h1 integrals: 98 call copvec(work(kint1o),work(kint1),nint1) 99 100 npairs = 4*ntoob**2 101 nheff = npairs*(npairs+1)/2 102 call memman(kheff,nheff,'ADDL ',2,'H_EFF ') 103 call memman(kheff2,npairs*npairs,'ADDL ',2,'H_EFF ') 104 nwmat = npairs*npairs 105 call memman(kwmat,nwmat,'ADDL ',2,'W_MAT ') 106 call memman(keig,npairs,'ADDL ',2,'EIGEN ') 107 108 npairs = (2*ntoob-1)*(2*ntoob)/2 109 nel = nactel 110 111 print *,' NEL = ', nel 112 113 call diag_effH(work(kheff),work(kheff2), 114 & work(kwmat),work(keig),nel) 115 116 call memchk2('a effH') 117 118c call logwmat(work(kwmat)) 119 120 call cleanvec(work(kwmat),npairs**2) 121 122 call rearr_g(work(kwmat),ccvec1, 123 & n_cc_typ,i_cc_typ,ioff_cc_typ) 124 125 igtbmod = 0 126 iopsym = 0 127 128 ! important: 129 call scalve(ccvec1,-1d0,n_cc_amp) 130 call gtbce_E(igtbmod,elen1,variance,ovl, 131 & ecore, 132 & ccvec1,iopsym,ccvec4, 133 & civec1,civec2,c2vec, 134 & n_cc_amp,mxb_ci, 135 & luc,luec,luhc,lusc1,lusc2) 136 h1 = elen1+ecore 137 138 ! effective Hamiltonian: 139 call rearr_g(work(kheff2),ccvec3, 140 & n_cc_typ,i_cc_typ,ioff_cc_typ) 141 142 call memchk2('a r g') 143 144 ! recalc energy with effective hamilt. 145 isigden = 1 146 call sigden_cc(civec1,civec2,luc,lusc1,ccvec3,isigden) 147 hf1 = ecore - inprdd(civec1,civec2,luc,lusc1,1,-1) 148 print *,'TEST: E(HF) = ',hf1 149 150 151 ! create a 2-fold excited determinant 152 ccvec2(ioff_cc_typ(25)+2) = 1d0 153 isigden = 1 154 call sigden_cc(civec1,civec2,luc,lusc1,ccvec2,isigden) 155 156 call gtbce_E(igtbmod,elen2,variance,ovl, 157 & ecore, 158 & ccvec1,iopsym,ccvec4, 159 & civec1,civec2,c2vec, 160 & n_cc_amp,mxb_ci, 161 & lusc1,lusc2,lusc3,lusc4,lusc5) 162 !ref eGref HeGref 163 164 h2 = elen2+ecore 165 ! <ex|e(-G)He(G)|0>: 166 h12 = inprdd(civec1,civec2,lusc2,luhc,1,-1) 167 h21 = inprdd(civec1,civec2,lusc3,luec,1,-1) 168 169 print *,' :: ',h1 , h21 170 print *,' :: ',h12, h2 171 172 ! Heff|expG 1> on lusc4 173 isigden = 1 174 call sigden_cc(civec1,civec2,luec,lusc4,ccvec3,isigden) 175 h1t = ecore - inprdd(civec1,civec2,luec,lusc4,1,-1) 176 ! Heff|expG 2> on lusc5 177 isigden = 1 178 call sigden_cc(civec1,civec2,lusc2,lusc5,ccvec3,isigden) 179 h2t = ecore - inprdd(civec1,civec2,lusc2,lusc5,1,-1) 180 181 h12t = -inprdd(civec1,civec2,lusc2,lusc4,1,-1) 182 h21t = -inprdd(civec1,civec2,luec,lusc5,1,-1) 183 184 print *,'TEST:' 185 print *,' :: ',h1t , h21t 186 print *,' :: ',h12t, h2t 187 188 189 call relunit(lusc1,'delete') 190 call relunit(lusc2,'delete') 191 call relunit(lusc3,'delete') 192 call relunit(lusc4,'delete') 193 call relunit(lusc5,'delete') 194c call relunit(lusc6,'delete') 195c call relunit(lusc7,'delete') 196c call relunit(lusc8,'delete') 197c call relunit(lusc9,'delete') 198 199 call memman(idum,idum,'FLUSM ',idum,'FOCKSP') 200 201 return 202 203 end 204 205************************************************************************ 206 207 subroutine diag_effH(h_eff,h_eff_out,wmat,eig,nel) 208************************************************************************ 209* 210* set up 211* H(pq,rs) = (pq|rs) + 1/(N-1)(h(pq)dlt(rs)+h(rs)dlt(pq)) 212* 213* and diagonalize in the sense: 214* 215* W(tu,pq)H(pq,rs)W(rs,vw) = e(tu)dlt(tu,vw) 216* 217************************************************************************ 218 219 include 'wrkspc.inc' 220 include 'crun.inc' 221 include 'cstate.inc' 222 include 'cgas.inc' 223 include 'ctcc.inc' 224 include 'gasstr.inc' 225 include 'strinp.inc' 226 include 'orbinp.inc' 227 include 'cprnt.inc' 228 include 'corbex.inc' 229 include 'csm.inc' 230 include 'cecore.inc' 231 include 'gtbce.inc' 232 include 'opti.inc' 233 include 'cintfo.inc' 234 include 'glbbas.inc' 235 include 'oper.inc' 236 237 dimension h_eff(*), wmat(*), eig(*), h_eff_out(*) 238 239 ! set up h_eff 240 241 npairs = 2*ntoob*(2*ntoob-1)/2 242 h_eff(1:npairs*(npairs+1)/2) = 0d0 243 xel = dble(nel) 244 245 ipr = 0 246 247 i_unrorb = 0 248 ispcas = 0 249 i_use_simtrh = 0 250 251 do mp = 1, 2 252c old: 253c do ip = 1, ntoob 254c do mr = 1, mp 255c 256c new: 257 do mr = 1, mp 258 do ip = 1, ntoob 259c 260 irmx = ntoob 261 if (mr.eq.mp) irmx = ip-1 262 do ir = 1, irmx 263 ipr = ipr + 1 264! ips = (imp-1)*2*ntoob + ims 265 266 iqs = 0 267 268 do mq = 1, 2 269c old: 270c do iq = 1, ntoob 271c do ms = 1, mq 272c 273c new: 274 do ms = 1, mq 275 do iq = 1, ntoob 276c 277 ismx = ntoob 278 if (ms.eq.mq) ismx = iq-1 279 do is = 1, ismx 280 iqs = iqs + 1 281! iqr = (imq-1)*2*ntoob + imr 282 283 if (ipr.lt.iqs) cycle ! hermitian operator 284 285 idxpqrs = ipr*(ipr-1)/2 + iqs 286 287c write(6,*) 'ip,iq,ir,is: ',ip,iq,ir,is 288c write(6,*) 'mp,mq,mr,ms: ',mp,mq,mr,ms 289cc 290c write(6,*) 'ipr,iqs,idxpqrs: ', ipr,iqs,idxpqrs 291cc 292 if (mp.eq.mq.and.mr.eq.ms) then 293c 294c write(6,*) ' before : ', h_eff(idxpqrs) 295c 296 h_eff(idxpqrs) = h_eff(idxpqrs) + gtijkl(ip,iq,ir,is) 297c 298c write(6,*) ' +2-el int : ', h_eff(idxpqrs) 299c 300 if (ip.eq.iq) 301 & h_eff(idxpqrs) = h_eff(idxpqrs) + 302 & 1d0/(xel-1d0)*geth1i_2(ir,is) 303 304 if (ir.eq.is) 305 & h_eff(idxpqrs) = h_eff(idxpqrs) + 306 & 1d0/(xel-1d0)*geth1i_2(ip,iq) 307 308c write(6,*) ' +1-el int : ', h_eff(idxpqrs) 309 end if 310 if (mp.eq.ms.and.mq.eq.mr) then 311c 312c write(6,*) ' before : ', h_eff(idxpqrs) 313c 314 h_eff(idxpqrs) = h_eff(idxpqrs) - gtijkl(ip,is,ir,iq) 315 316c write(6,*) ' +2-el int : ', h_eff(idxpqrs) 317 318 if (ip.eq.is) 319 & h_eff(idxpqrs) = h_eff(idxpqrs) - 320 & 1d0/(xel-1d0)*geth1i_2(ir,iq) 321 322 if (ir.eq.iq) 323 & h_eff(idxpqrs) = h_eff(idxpqrs) - 324 & 1d0/(xel-1d0)*geth1i_2(ip,is) 325 326c write(6,*) ' +1-el int : ', h_eff(idxpqrs) 327 328 end if 329 330 idx2pqrs = (ipr-1)*npairs+iqs 331 idx2qpsr = (iqs-1)*npairs+ipr 332 333 h_eff_out(idx2pqrs) = h_eff(idxpqrs) 334 h_eff_out(idx2qpsr) = h_eff(idxpqrs) 335 336 end do 337 end do 338 end do 339 end do 340 end do 341 end do 342 end do 343 end do 344 345 write(6,*) 'the hamiltonian:' 346c write(6,*) ' norm^2 = ',ddot(npairs**2,h_eff,1,h_eff,1) 347 348 call prtrlt(h_eff,npairs) 349 350 call copdia(h_eff,eig,npairs,1) 351 352 write(6,*) 'the diagonal of the hamiltonian:' 353 354 call wrtmat(eig,npairs,1,npairs,1) 355 356 ! call diagonalizer 357 358c call eigen(h_eff,wmat,npairs,0,1) 359 call eigen(h_eff,wmat,npairs,0,0) 360 361 call copdia(h_eff,eig,npairs,1) 362 363 write(6,*) 'the eigenvalues of the hamiltonian:' 364 365 call wrtmat(eig,npairs,1,npairs,1) 366 367 call memman(keigr,npairs,'ADDL ',2,'EIG RE') 368 call memman(keigi,npairs,'ADDL ',2,'EIG IM') 369 370 call dumsort(eig,npairs,npairs,work(keigr),work(keigi)) 371 call reo_vec(work(keigr),eig,npairs,work(keigi),1) 372 call copvec(work(keigi),eig,npairs) 373 374 write(6,*) 'the eigenvalues of the hamiltonian (sorted):' 375 376 call wrtmat(eig,npairs,1,npairs,1) 377 378 write(6,*) 'E_nuc = ',ecore 379 380 ener = ecore 381 do ii = 1, nel*(nel-1)/2 382 ener = ener + eig(ii) 383 print *,ii,eig(ii),ener 384 end do 385 386 write(6,*) 'the W-matrix: ',sqrt(ddot(npairs**2,wmat,1,wmat,1)) 387 388 call wrtmat2(wmat,npairs,npairs,npairs,npairs) 389 390c return 391 392c call memman(keigr,npairs,'ADDL ',2,'EIG RE') 393c call memman(keigi,npairs,'ADDL ',2,'EIG IM') 394 call memman(klneig,npairs,'ADDL ',2,'EIG IM') 395 call memman(keigv,npairs**2,'ADDL ',2,'EIGVEC') 396 call memman(keigvr,npairs**2,'ADDL ',2,'EIGVEC') 397 call memman(keigvi,npairs**2,'ADDL ',2,'EIGVEC') 398 call memman(keigvr2,npairs**2,'ADDL ',2,'EIGVEC') 399 call memman(keigvi2,npairs**2,'ADDL ',2,'EIGVEC') 400 call memman(kgmatr,npairs**2,'ADDL ',2,'GMAT R') 401 call memman(kgmati,npairs**2,'ADDL ',2,'GMAT I') 402 call memman(kiscr,npairs,'ADDL ',1,'INTSCR') 403 call memman(kscr,npairs,'ADDL ',2,'RELSCR') 404 405 call logumat(npairs,work(kgmatr),wmat, 406 & work(keigv),work(keigvr),work(keigvi)) 407c 408c call copvec(wmat,work(keigvr),npairs**2) 409c 410c call rg(npairs,npairs,work(keigvr), 411c & work(keigr),work(keigi),1,work(keigv), 412c & work(kiscr),work(kscr),ierr) 413c call nrmvec(npairs,work(keigv),work(keigi)) 414c 415c print *,'eigenvalues of W-matrix:' 416c 417c do ii = 1, npairs 418c eigr = work(keigr-1+ii) 419c eigi = work(keigi-1+ii) 420c rad = eigr*eigr+eigi*eigi 421c costau = eigr 422c tau = acos(eigr)*sign(1d0,eigi) 423c work(klneig-1+ii) = tau 424c print '(i4,2(2x,e20.10),2(2x,f15.10))',ii,eigr,eigi,rad, 425c & tau/3.1415926535897932d0*180d0 426c end do 427c 428cc print *,'the eigenvectors' 429cc call wrtmat2(work(keigv),npairs,npairs,npairs,npairs) 430c 431c ! sort components into real and imaginary matrices 432c 433c ii = 0 434c do while(ii.lt.npairs) 435c ii = ii+1 436c eigr = work(keigr-1+ii) 437c eigi = work(keigi-1+ii) 438c if (eigi.eq.0d0) then 439c call copvec(work(keigv+(ii-1)*npairs), 440c & work(keigvr+(ii-1)*npairs),npairs) 441c call setvec(work(keigvi+(ii-1)*npairs),0d0,npairs) 442c else 443c call copvec(work(keigv+(ii-1)*npairs), 444c & work(keigvr+(ii-1)*npairs),npairs) 445c call copvec(work(keigv+(ii-1)*npairs), 446c & work(keigvr+(ii)*npairs),npairs) 447c call copvec(work(keigv+(ii)*npairs), 448c & work(keigvi+(ii-1)*npairs),npairs) 449c call copvec(work(keigv+(ii)*npairs), 450c & work(keigvi+(ii)*npairs),npairs) 451c call scalve(work(keigvi+(ii)*npairs),-1d0,npairs) 452c ii = ii+1 ! additional increment 453c end if 454c 455c end do 456c 457cc print *,'the eigenvectors (re):' 458cc call wrtmat2(work(keigvr),npairs,npairs,npairs,npairs) 459cc print *,'the eigenvectors (im):' 460cc call wrtmat2(work(keigvi),npairs,npairs,npairs,npairs) 461c 462cc ! test the trafo matrices: 463cc 464cc ! A^T A - B^T B: 465cc call dgemm('t','n',npairs,npairs,npairs, 466cc & 1d0,work(keigvr),npairs, 467cc & work(keigvr),npairs, 468cc & 0d0,work(keigvr2),npairs) 469cc call dgemm('t','n',npairs,npairs,npairs, 470cc & 1d0,work(keigvi),npairs, 471cc & work(keigvi),npairs, 472cc & 1d0,work(keigvr2),npairs) 473cc 474cc ! A^T B + B^T A: 475cc call dgemm('t','n',npairs,npairs,npairs, 476cc & 1d0,work(keigvr),npairs, 477cc & work(keigvi),npairs, 478cc & 0d0,work(keigvi2),npairs) 479cc call dgemm('t','n',npairs,npairs,npairs, 480cc & -1d0,work(keigvi),npairs, 481cc & work(keigvr),npairs, 482cc & 1d0,work(keigvi2),npairs) 483cc 484cc print *,'U^T U, real part:' 485cc call wrtmat2(work(keigvr2),npairs,npairs,npairs,npairs) 486cc 487cc print *,'U^T U, imaginary part:' 488cc call wrtmat2(work(keigvi2),npairs,npairs,npairs,npairs) 489cc 490cc 491cc ! W A 492cc call dgemm('n','n',npairs,npairs,npairs, 493cc & 1d0,wmat,npairs, 494cc & work(keigvr),npairs, 495cc & 0d0,work(keigvr2),npairs) 496cc 497cc ! W B 498cc call dgemm('n','n',npairs,npairs,npairs, 499cc & 1d0,wmat,npairs, 500cc & work(keigvi),npairs, 501cc & 0d0,work(keigvi2),npairs) 502cc 503cc ! Re D = A^T W A - B^T W B 504cc call dgemm('t','n',npairs,npairs,npairs, 505cc & 1d0,work(keigvi),npairs, 506cc & work(keigvi2),npairs, 507cc & 0d0,work(kgmatr),npairs) 508cc call dgemm('t','n',npairs,npairs,npairs, 509cc & 1d0,work(keigvr),npairs, 510cc & work(keigvr2),npairs, 511cc & 1d0,work(kgmatr),npairs) 512cc 513cc ! Im D = A^T W B + B^T W A 514cc call dgemm('t','n',npairs,npairs,npairs, 515cc & 1d0,work(keigvr),npairs, 516cc & work(keigvi2),npairs, 517cc & 0d0,work(kgmati),npairs) 518cc call dgemm('t','n',npairs,npairs,npairs, 519cc & -1d0,work(keigvi),npairs, 520cc & work(keigvr2),npairs, 521cc & 1d0,work(kgmati),npairs) 522cc 523cc print *,'U^T W U, real part:' 524cc xnrm = sqrt(ddot(npairs**2,work(kgmatr),1,work(kgmatr),1 )) 525cc print *,'norm = ',xnrm 526cc 527cc call wrtmat2(work(kgmatr),npairs,npairs,npairs,npairs) 528cc 529cc print *,'U^T W U, imaginary part:' 530cc call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs) 531c 532c ! D A^+ 533c do ii = 1, npairs 534c tau = work(klneig-1+ii) 535c do jj = 1, npairs 536c ijadr = (jj-1)*npairs+ii-1 ! uahhh 537c jiadr = (ii-1)*npairs+jj-1 ! uahhh 538c work(keigvr2+ijadr) = tau * work(keigvr+jiadr) 539c end do 540c end do 541c 542c ! D B^+ 543c do ii = 1, npairs 544c tau = work(klneig-1+ii) 545c do jj = 1, npairs 546c ijadr = (jj-1)*npairs+ii-1 547c jiadr = (ii-1)*npairs+jj-1 548c work(keigvi2+ijadr) = tau * work(keigvi+jiadr) 549c end do 550c end do 551c 552c ! Re G = - B D A^+ + A D B^+ 553c call dgemm('n','n',npairs,npairs,npairs, 554c & -1d0,work(keigvi),npairs, 555c & work(keigvr2),npairs, 556c & 0d0,work(kgmatr),npairs) 557c call dgemm('n','n',npairs,npairs,npairs, 558c & +1d0,work(keigvr),npairs, 559c & work(keigvi2),npairs, 560c & 1d0,work(kgmatr),npairs) 561c 562c ! Im G = A D A^+ + B D B^+ 563c call dgemm('n','n',npairs,npairs,npairs, 564c & 1d0,work(keigvr),npairs, 565c & work(keigvr2),npairs, 566c & 0d0,work(kgmati),npairs) 567c call dgemm('n','n',npairs,npairs,npairs, 568c & 1d0,work(keigvi),npairs, 569c & work(keigvi2),npairs, 570c & 1d0,work(kgmati),npairs) 571c 572c print *,'G, real part:' 573c 574c call wrtmat2(work(kgmatr),npairs,npairs,npairs,npairs) 575c 576c print *,'G, imaginary part:' 577c xnrm = sqrt(ddot(npairs**2,work(kgmatr),1,work(kgmatr),1 )) 578c print *,'norm = ',xnrm 579c 580c if (xnrm.gt.1d-10) 581c & call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs) 582 583 call expmat(work(kgmati),work(kgmatr),work(keigv),work(keigvr), 584 & npairs,1d-20) 585 586 print *,'W from exp(G):' 587 call wrtmat2(work(kgmati),npairs,npairs,npairs,npairs) 588** 589 590 call copvec(work(kgmatr),wmat,npairs**2) 591 592c stop 'end of test' 593 594 end 595 596 subroutine expmat(expx,xmat,xscr1,xscr2,ndim,thrsh) 597 598* calculate exp(X), returned on expx, of (ndim,ndim)-matrix X, 599* input on xmat, by Taylor expansion (threshold thrsh) 600* xscr is a scratch matrix of the same dimensions as xmat, expx 601 602 implicit none 603 604 integer, parameter :: ntest = 100, maxn = 100 605 606 integer, intent(in) :: 607 & ndim 608 real(8), intent(in) :: 609 & thrsh 610 real(8), intent(inout) :: 611 & expx(ndim,ndim), xmat(ndim,ndim), 612 & xscr1(ndim,ndim), xscr2(ndim,ndim) 613 614 logical :: 615 & conv 616 integer :: 617 & n, ndim2, ii 618 real(8) :: 619 & xnrm, fac 620 621 real(8), external :: 622 & inprod 623 624 expx(1:ndim,1:ndim) = xmat(1:ndim,1:ndim) 625 xscr2(1:ndim,1:ndim) = xmat(1:ndim,1:ndim) 626 627 do ii = 1, ndim 628 expx(ii,ii) = expx(ii,ii) + 1d0 629 end do 630 631 ndim2 = ndim*ndim 632 n = 1 633 conv = .false. 634 635 do while (.not.conv) 636 n = n+1 637 if (n.gt.maxn) exit 638 639 fac = 1d0/dble(n) 640 641 ! Xscr = 1/N Xscr * X 642 call matml7(xscr1,xscr2,xmat, 643 & ndim,ndim, 644 & ndim,ndim, 645 & ndim,ndim,0d0,fac,0) 646 647 xnrm = sqrt(inprod(xscr1,xscr1,ndim2)) 648 if (xnrm.lt.thrsh) conv = .true. 649 650 if (ntest.ge.10) 651 & write(6,*) ' N = ',n,' |1/N! X^N| = ',xnrm 652 653 call vecsum(expx,expx,xscr1,1d0,1d0,ndim2) 654 655 call copvec(xscr1,xscr2,ndim2) 656 657 end do 658 659 if (.not.conv) then 660 write(6,*) ' Taylor expansion of exp(X) did not converge!' 661 stop 'expmat' 662 end if 663 664 return 665 end 666c 667 subroutine rearr_g(gmat_in,gmat_out,ntss_tp,itss_tp,ibtss_tp) 668 669 include 'implicit.inc' 670 include 'mxpdim.inc' 671 include 'cgas.inc' 672 include 'multd2h.inc' 673 include 'orbinp.inc' 674 include 'csm.inc' 675 include 'ctcc.inc' 676 include 'cc_exc.inc' 677 678 integer, parameter :: 679 & ntest = 1000 680 681* input 682 real(8), intent(inout) :: 683 & gmat_in(2*ntoob*(2*ntoob-1)/2,2*ntoob*(2*ntoob-1)/2), 684 & gmat_out(*) 685 686c input needed: itss_tp <-- work(klsobex), ntss_tp <-- nspobex_tp 687 integer, intent(in) :: 688 & ntss_tp, 689 & itss_tp(ngas,4,ntss_tp), 690 & ibtss_tp(ntss_tp) 691 692* local 693 integer :: 694 & igrp_ca(mxpngas), igrp_cb(mxpngas), 695 & igrp_aa(mxpngas), igrp_ab(mxpngas), 696 & iocc_ca(mx_st_tsoso_blk_mx), 697 & iocc_cb(mx_st_tsoso_blk_mx), 698 & iocc_aa(mx_st_tsoso_blk_mx), 699 & iocc_ab(mx_st_tsoso_blk_mx), 700 & idx_c(4), idx_s(4) 701 702 703 npairs = 2*ntoob*(2*ntoob-1)/2 704 705 if (ntest.ge.1000) then 706 write(6,*) ' input amplitudes: ' 707 call wrtmat2(gmat_in,npairs,npairs,npairs,npairs) 708 write(6,*) 'ibtss_tp:' 709 call iwrtma(ibtss_tp,1,ntss_tp,1,ntss_tp) 710 end if 711 712 713 ! loop over types 714 do itss = 1, ntss_tp 715 idx = ibtss_tp(itss) - 1 716 717 ! identify two-particle excitations: 718 nel_ca = ielsum(itss_tp(1,1,itss),ngas) 719 nel_cb = ielsum(itss_tp(1,2,itss),ngas) 720 nel_aa = ielsum(itss_tp(1,3,itss),ngas) 721 nel_ab = ielsum(itss_tp(1,4,itss),ngas) 722 nc = nel_ca + nel_cb 723 na = nel_aa + nel_ab 724 if (na.ne.2) cycle 725 ! transform occupations to groups 726 call occ_to_grp(itss_tp(1,1,itss),igrp_ca,1) 727 call occ_to_grp(itss_tp(1,2,itss),igrp_cb,1) 728 call occ_to_grp(itss_tp(1,3,itss),igrp_aa,1) 729 call occ_to_grp(itss_tp(1,4,itss),igrp_ab,1) 730 731 if (mscomb_cc.ne.0) then 732 call diag_exc_cc(itss_tp(1,1,itss),itss_tp(1,2,itss), 733 & itss_tp(1,3,itss),itss_tp(1,4,itss), 734 & ngas,idiag) 735 else 736 idiag = 0 737 end if 738 739 ! loop over symmetry blocks 740 ism = 1 ! totally symmetric operators 741 do ism_c = 1, nsmst 742 ism_a = multd2h(ism,ism_c) 743 do ism_ca = 1, nsmst 744 ism_cb = multd2h(ism_c,ism_ca) 745 do ism_aa = 1, nsmst 746 ism_ab = multd2h(ism_a,ism_aa) 747 ! get alpha and beta symmetry index 748 ism_alp = (ism_aa-1)*nsmst+ism_ca ! = (sym Ca,sym Aa) 749 ism_bet = (ism_ab-1)*nsmst+ism_cb ! = (sym Cb,sym Ab) 750 751 ! restrict to (sym Ca,sym Aa) >= (sym Cb,sym Ab) 752 if (idiag.eq.1.and.ism_bet.gt.ism_alp) cycle 753 if (idiag.eq.0.or.ism_alp.gt.ism_bet) then 754 irestr = 0 755 else 756 irestr = 1 757 end if 758 759 ! get the strings 760 call getstr2_totsm_spgp(igrp_ca,ngas,ism_ca,nel_ca, 761 & lca,iocc_ca,norb,0,idum,idum) 762 call getstr2_totsm_spgp(igrp_cb,ngas,ism_cb,nel_cb, 763 & lcb,iocc_cb,norb,0,idum,idum) 764 call getstr2_totsm_spgp(igrp_aa,ngas,ism_aa,nel_aa, 765 & laa,iocc_aa,norb,0,idum,idum) 766 call getstr2_totsm_spgp(igrp_ab,ngas,ism_ab,nel_ab, 767 & lab,iocc_ab,norb,0,idum,idum) 768 769 ! length of strings in this symmetry block 770 if (lca*lcb*laa*lab.eq.0) cycle 771 772 do iab = 1, lab 773 if (irestr.eq.1) then 774 iaa_min = iab 775 else 776 iaa_min = 1 777 end if 778 do iaa = iaa_min, laa 779 do icb = 1, lcb 780 if (irestr.eq.1.and.iaa.eq.iab) then 781 ica_min = icb 782 else 783 ica_min = 1 784 end if 785 do ica = ica_min, lca 786 idx = idx + 1 787 ! translate into canonical index quadrupel 788 ii = 0 789 do iel = 1, nel_ca 790 ii = ii + 1 791 idx_c(ii) = iocc_ca((ica-1)*nel_ca+iel) 792 idx_s(ii) = 1 793 end do 794 do iel = 1, nel_cb 795 ii = ii + 1 796 idx_c(ii) = iocc_cb((icb-1)*nel_cb+iel) 797 idx_s(ii) = 2 798 end do 799 do iel = 1, nel_aa 800 ii = ii + 1 801 idx_c(ii) = iocc_aa((iaa-1)*nel_aa+iel) 802 idx_s(ii) = 1 803 idx_s(ii) = 1 804 end do 805 do iel = 1, nel_ab 806 ii = ii + 1 807 idx_c(ii) = iocc_ab((iab-1)*nel_ab+iel) 808 idx_s(ii) = 2 809 end do 810 811 ! idx_c contains: 812 ! p1, p2, h1, h2 813 ! idx_s contains 814 ! mp1, mp2, mh1, mh2 815 816 print *,'p, r: ',idx_c(1),idx_c(2) 817 print *,' ',idx_s(1),idx_s(2) 818 print *,'q, s: ',idx_c(3),idx_c(4) 819 print *,' ',idx_s(3),idx_s(4) 820 821 idxpr = igtidx(idx_c(1),idx_c(2), 822 & idx_s(1),idx_s(2),ntoob) 823 idxqs = igtidx(idx_c(3),idx_c(4), 824 & idx_s(3),idx_s(4),ntoob) 825 826 print *,' >>> ',idx,' ---> ',idxpr,idxqs 827 828 gmat_out(idx) = gmat_in(idxpr,idxqs) 829 830c testing coverage: 831 if (gmat_in(idxpr,idxqs).ne.1d100) then 832 gmat_in(idxpr,idxqs) = 1d100 833 else 834 gmat_in(idxpr,idxqs) = 835 & gmat_in(idxpr,idxqs)+1d100 836 end if 837 838 end do ! ica 839 end do ! icb 840 end do ! iaa 841 end do ! iab 842 843 end do ! ism_aa 844 end do ! ism_ca 845 end do ! ism_c 846 847 end do ! itss 848 849 if (ntest.ge.1000) then 850 write(6,*) 'coverage:' 851 call wrtmat2(gmat_in,npairs,npairs,npairs,npairs) 852 write(6,*) ' output amplitudes: ' 853 call wrt_cc_vec2(gmat_out,6,'GEN_CC') 854 end if 855 856 return 857 end 858 859 integer function igtidx(ip,ir,mp,mr,ndim) 860 861 implicit none 862 863 integer, intent(in) :: 864 & ip, ir, 865 & mp, mr, ndim 866 867 integer :: 868 & ioff, idx1, idx2 869 870 ! get ipr: 871 if (mp.eq.1.and.mr.eq.1) then 872 ioff = 0 873 idx1 = max(ip,ir) 874 idx2 = min(ip,ir) 875 else if (mp.eq.1.and.mr.eq.2) then 876 ioff = ndim*(ndim-1)/2 ! NB: no diagonal in a/a block 877 idx1 = ir 878 idx2 = ip 879 else if (mp.eq.2.and.mr.eq.1) then 880 ioff = ndim*(ndim-1)/2 881 idx1 = ip 882 idx2 = ir 883 else if (mp.eq.2.and.mr.eq.2) then 884 ioff = ndim*(ndim-1)/2 + ndim*ndim 885 idx1 = max(ip,ir) 886 idx2 = min(ip,ir) 887 else 888 stop 'mp, mr ?' 889 end if 890 891 print *,'ioff,idx1,idx2: ',ioff,idx1,idx2 892 893 if (mp.eq.mr) then 894 if (ip.eq.ir) stop 'ip, ir ?' 895 if (idx1.eq.1) stop 'max(ip,ir) ?' 896 igtidx = ioff + (idx1-1)*(idx1-2)/2+idx2 897 else 898 igtidx = ioff + (idx1-1)*ndim+idx2 899 end if 900 901 print *,' >> gtidx = ',igtidx 902 903 return 904 905 end 906 907 subroutine cleanvec(vec,ndim) 908 909 implicit none 910 911 integer, intent(in) :: 912 & ndim 913 real(8), intent(inout) :: 914 & vec(ndim) 915 916 integer :: 917 & ii 918 real(8) :: 919 & thr 920 921 thr = 1000d0*epsilon(1d0) 922 923 print *,'thr = ',thr 924 print *,'ndim = ',ndim 925 926 do ii = 1, ndim 927 if (abs(vec(ii)).lt.thr) vec(ii)=0d0 928 end do 929 930 return 931 932 end 933c $Id$ 934