1#define DETCI_INTFACE 2c 3c Need this macro to exclude declarations 4c in detci.fh 5c 6 7 subroutine detci_init( norb, nela, nelb, nsym, symstate, 8 $ osym, iprint, eps, h, g ) 9* 10* $Id$ 11* 12 implicit none 13#include "errquit.fh" 14c 15#include "mafdecls.fh" 16#include "global.fh" 17#include "util.fh" 18#include "detciP.fh" 19#include "detci.fh" 20#include "cdetci.fh" 21#include "cdetcistats.fh" 22c 23 integer norb 24 integer nela 25 integer nelb 26 integer nsym 27 integer symstate 28 integer osym(norb) 29 integer iprint 30 double precision eps(norb) 31 double precision h(*) 32 double precision g(*) 33c 34c 35c 36 double precision tx 37 integer nexa, nexb 38 integer nstra, nstrb 39 integer nekla, neklb 40 integer ntij, civlen 41 integer ijmap(detci_maxorb*detci_maxorb) 42 integer vtab((detci_maxorb+1)*(detci_maxelec+1)*(detci_maxsy)) 43 integer i, j, ii 44c$$$ integer ij 45 integer memuse03 46 integer mseg 47c 48c Reset statistics 49c 50 tx = util_cpusec() 51 detci_init_etime = 0.d0 52 detci_saa_etime = 0.d0 53 detci_sbb_etime = 0.d0 54 detci_sab_etime = 0.d0 55 detci_sigma_etime = 0.d0 56 detci_density_etime = 0.d0 57 detci_sigma_calls = 0 58 detci_density_calls = 0 59 detci_spinadp_etime = 0.d0 60 detci_symmadp_etime = 0.d0 61 detci_aaff_etime = 0.d0 62 detci_aadot_etime = 0.d0 63 detci_aagop_etime = 0.d0 64 detci_abstr_etime = 0.d0 65 detci_abgath_etime = 0.d0 66 detci_abdotab_etime = 0.d0 67 detci_abscat_etime = 0.d0 68 detci_absync_etime = 0.d0 69 detci_density_onepdm = 0.0d0 70 detci_density_twopdm = 0.0d0 71 detci_density_twopdmab = 0.0d0 72 cdetci_profprint = .true. 73c 74c Copy CI parameters 75c 76 cdetci_norb = norb 77 cdetci_nela = nela 78 cdetci_nelb = nelb 79 cdetci_nsym = nsym 80 cdetci_symstate = symstate 81c 82c Reorder orbital indices by symmetry 83c 84 ii = 0 85 do j=1,nsym 86 do i=1,norb 87 if (osym(i).eq.j) then 88 ii = ii + 1 89 cdetci_ixmap(ii) = i 90 cdetci_irmap(i) = ii 91 cdetci_osym(ii) = osym(i) 92 cdetci_eps(ii) = eps(i) 93 endif 94 enddo 95 enddo 96c 97c$$$ DO I=1,NORB 98c$$$ WRITE(6,'(I5,A5,I5,5X,F10.4)') I,' ->',CDETCI_IRMAP(I), 99c$$$ $ CDETCI_EPS(CDETCI_IRMAP(I)) 100c$$$ ENDDO 101c 102c Compute CI parameters 103c 104 call detci_ijmap( norb, nsym, cdetci_osym, ntij, ijmap ) 105 nexa = (norb-nela+1)*nela 106 nexb = (norb-nelb+1)*nelb 107 nstra = detci_binomial(norb,nela) 108 nstrb = detci_binomial(norb,nelb) 109 nekla = detci_binomial((norb-1),(nela-1)) 110 neklb = detci_binomial((norb-1),(nelb-1)) 111 civlen = nstra*nstrb 112 mseg = (10*ga_nnodes())/max(min(nexa,nexb),1) 113 mseg = min(mseg, 20) 114 MSEG = 1 ! Delete for better parallel effeciency 115c 116c Default reference energy is Aufbau ordering 117c ...although maybe reset by CI-guess routine! 118c 119 cdetci_eref = detci_refenergy( norb, nela, nelb, eps ) 120c 121c Toggle spin-adaption 122c Need some way of automatically turning this on and off 123c depending requested spin-state 124c 125 cdetci_spinadapt = .true. 126 cdetci_squant = (nela - nelb)/2.d0 127c 128c Info print 129c 130 if ((iprint.gt.0).and.(ga_nodeid().eq.0)) then 131 write(6,902) norb, nsym, civlen, 132 $ cdetci_symstate, 133 $ cdetci_spinadapt, 134 $ cdetci_squant, 135 $ nela, nelb, nstra, nstrb, 136 $ nexa, nexb, nekla, neklb 137 902 format(10x,'Active shells:',23x,i5,/, 138 $ 10x,'Irreps:',30x,i5,/, 139 $ 10x,'CI vector length:',5x,i20,/, 140 $ 10x,'State symmetry:',25x,i2,/, 141 $ 10x,'Spin adaption:',27x,l1,/, 142 $ 10x,'S quantum number:',15x,f10.3,//, 143 $ 37x,'Alpha',6x,'Beta',/, 144 $ 35x,2(7('-'),3x),/, 145 $ 10x,'Electrons:',12x,2i10,/, 146 $ 10x,'Strings:',14x,2i10,/, 147 $ 10x,'E_ij per string:',6x,2i10,/, 148 $ 10x,'Strings per E_ij:',5x,2i10,/) 149 call util_flush(6) 150 endif 151c 152c Construct arc weight tables 153c 154 call detci_vatable( norb, nela, nsym, cdetci_osym, 155 $ vtab, cdetci_ataba ) 156 call detci_vatable( norb, nelb, nsym, cdetci_osym, 157 $ vtab, cdetci_atabb ) 158c 159c Allocate and construct excitation operator table 160c 161 l_detci_exa = CDETCI_INVALID 162 l_detci_exb = CDETCI_INVALID 163 if (nexa.gt.0) then 164 if (.not.ma_push_get(MT_INT, (6*nexa*(nstra+1)), 'detci:exa', 165 $ l_detci_exa, k_detci_exa)) 166 $ call errquit('detci: cannot allocate',0, MA_ERR) 167 call detci_excit( norb, nela, nsym, nstra, nexa, cdetci_osym, 168 $ ijmap, cdetci_ataba, int_mb(k_detci_exa) ) 169 endif 170 if (nexb.gt.0) then 171 if (.not.ma_push_get(MT_INT, (6*nexb*(nstrb+1)), 'detci:exb', 172 $ l_detci_exb, k_detci_exb)) 173 $ call errquit('detci: cannot allocate',0, MA_ERR) 174 call detci_excit( norb, nelb, nsym, nstrb, nexb, cdetci_osym, 175 $ ijmap, cdetci_atabb, int_mb(k_detci_exb) ) 176 endif 177c$$$ IF (GA_NODEID().EQ.0) THEN 178c$$$ IJ = 6*(NEXA*(NSTRA+1) + NEXB*(NSTRB+1)) 179c$$$ I = MA_INQUIRE_AVAIL( MT_DBL ) 180c$$$ WRITE(6,913) IJ, I 181c$$$ 913 FORMAT('ALLOCATED AND GENERATED EXCITATION TABLE',/, 182c$$$ $ 'Requested ',I10,' Free',I10,' words') 183c$$$ CALL UTIL_FLUSH(6) 184c$$$ ENDIF 185c$$$ CALL GA_SYNC() 186c 187c Allocate ERI block 188c Copy integrals into internal block 189c 190 if (.not.ma_push_get(MT_DBL,ntij,'detci: h1', 191 $ l_detci_h, k_detci_h)) 192 $ call errquit('detci: cannot allocate',0, MA_ERR) 193 if (.not.ma_push_get( MT_DBL,(ntij*ntij),'detci: eri', 194 $ l_detci_g, k_detci_g)) 195 $ call errquit('detci: cannot allocate',0, MA_ERR) 196 call detci_moint_copy( norb, ntij, cdetci_ixmap, ijmap, 197 $ cdetci_osym, h, g, 198 $ dbl_mb(k_detci_h), dbl_mb(k_detci_g) ) 199c 200c Memory for scratch space should be allocated HERE! 201c 202 memuse03 = 2*max(nstra,nstrb)*mseg + 3*nekla + 4*nstrb*nekla 203 if (.not.ma_push_get(MT_DBL, (max(max(nstra,nstrb)*mseg,nexb)), 204 $ 'detci: fscr', l_detci_fscr, k_detci_fscr)) 205 $ call errquit('detci: cannot allocate fscr',0, MA_ERR) 206 207 if (.not.ma_push_get(MT_INT, nekla, 208 $ 'detci: rhsscr', l_detci_rhsscr, k_detci_rhsscr)) 209 $ call errquit('detci: cannot allocate rhs',0, MA_ERR) 210 211 if (.not.ma_push_get(MT_INT, max(nekla,neklb,nexb,nexa), 212 $ 'detci: lhsscr', l_detci_lhsscr, k_detci_lhsscr)) 213 $ call errquit('detci: cannot allocate lhs',0, MA_ERR) 214 215 if (.not.ma_push_get(MT_INT, nekla, 216 $ 'detci: iscr', l_detci_iscr, k_detci_iscr)) 217 $ call errquit('detci: cannot allocate iscr',0, MA_ERR) 218 219 if (.not.ma_push_get(MT_DBL, (nstrb*nekla), 220 $ 'detci: cprime', l_detci_cprime, k_detci_cprime)) 221 $ call errquit('detci: cannot allocate cprime',0, MA_ERR) 222 223 if (.not.ma_push_get(MT_DBL, (nstrb*nekla), 224 $ 'detci: sprime', l_detci_sprime, k_detci_sprime)) 225 $ call errquit('detci: cannot allocate sprime',0, MA_ERR) 226c 227c$$$ IF (GA_NODEID().EQ.0) THEN 228c$$$ IJ = MAX(NSTRA,NSTRB) + 3*NEKLA + 2*(NSTRB*NEKLA) 229c$$$ I = MA_INQUIRE_AVAIL( MT_DBL ) 230c$$$ WRITE(6,914) IJ, I 231c$$$ 914 FORMAT('ALLOCATED SCRATCH ARRAYS',/, 232c$$$ $ 'Requested ',I10,' Free',I10,' words') 233c$$$ CALL UTIL_FLUSH(6) 234c$$$ ENDIF 235c$$$ CALL GA_SYNC() 236c 237c Validate internals 238c 239 cdetci_valid = CDETCI_MAGIC 240 detci_init_etime = util_cpusec() - tx 241 return 242 end 243 244 245 246 247 248 249c 250c Compute internal memory required & CI-vector length 251c 252 253 integer function detci_memsiz( norb, nela, nelb, nsym, 254 $ osym, vlen ) 255 implicit none 256 integer norb 257 integer nela 258 integer nelb 259 integer nsym 260 integer osym(norb) 261 integer vlen 262c 263 integer nexa 264 integer nexb 265 integer nstra 266 integer nstrb 267 integer nekla 268 integer neklb 269 integer detci_binomial 270 external detci_binomial 271c 272c 273c 274 nexa = (norb-nela+1)*nela 275 nexb = (norb-nelb+1)*nelb 276 nstra = detci_binomial(norb,nela) 277 nstrb = detci_binomial(norb,nelb) 278 nekla = detci_binomial((norb-1),(nela-1)) 279 neklb = detci_binomial((norb-1),(nelb-1)) 280 281 vlen = nstra*nstrb 282 detci_memsiz = 0 283 return 284 end 285 286 287 288 289 290 subroutine detci_free() 291 implicit none 292#include "errquit.fh" 293c 294#include "mafdecls.fh" 295#include "global.fh" 296#include "msgids.fh" 297#include "detciP.fh" 298#include "detci.fh" 299#include "cdetci.fh" 300#include "cdetcistats.fh" 301c 302c 303c Check if initialized 304c 305 if (cdetci_valid .ne. CDETCI_MAGIC) 306 $ call errquit('detci_free with uninitialized internals',0, 307 & INPUT_ERR) 308c 309c Print stats 310c 311 if (cdetci_profprint) then 312 if (ga_nodeid().eq.0) write(6,901) detci_sigma_calls, 313 $ detci_saa_etime, 314 $ detci_sbb_etime, 315 $ detci_sab_etime, 316 $ detci_sigma_etime, 317 $ detci_aaff_etime, 318 $ detci_aagop_etime, 319 $ detci_aadot_etime, 320 $ detci_abstr_etime, 321 $ detci_abgath_etime, 322 $ detci_abdotab_etime, 323 $ detci_abscat_etime, 324 $ detci_absync_etime, 325 $ detci_density_etime, 326 $ detci_density_onepdm, 327 $ detci_density_twopdm, 328 $ detci_density_twopdmab, 329 $ detci_spinadp_etime, 330 $ detci_symmadp_etime 331 901 format(/,10x,'Number of sigma calls:',4x,i5, 332 $ /,23x,'o',5('<'),' (aa):',7x,f10.2, 333 $ /,23x,'o',5('<'),' (bb):',7x,f10.2, 334 $ /,23x,'o',5('<'),' (ab):',7x,f10.2, 335 $ /,23x,'o',5('<'),' (total)',5x,f10.2, 336 $ /,23x,'o',5('<'),' (aa) ff',5x,f10.2, 337 $ /,23x,'o',5('<'),' (aa) gop',4x,f10.2, 338 $ /,23x,'o',5('<'),' (aa) dot',4x,f10.2, 339 $ /,23x,'o',5('<'),' (ab) str',4x,f10.2, 340 $ /,23x,'o',5('<'),' (ab) gath',3x,f10.2, 341 $ /,23x,'o',5('<'),' (ab) dotab',2x,f10.2, 342 $ /,23x,'o',5('<'),' (ab) scat',3x,f10.2, 343 $ /,23x,'o',5('<'),' (ab) sync',3x,f10.2, 344 $ /,23x,'o',5('<'),' Density',5x,f10.2, 345 $ /,23x,'o',5('<'),' Density one',1x,f10.2, 346 $ /,23x,'o',5('<'),' Density two',1x,f10.2, 347 $ /,23x,'o',5('<'),' Density ab',2x,f10.2, 348 $ /,23x,'o',5('<'),' Spin adapt',2x,f10.2, 349 $ /,23x,'o',5('<'),' Symm adapt',2x,f10.2) 350 351 call ga_dgop(msg_detci_maxsync, detci_absync_etime, 352 $ 1, 'max' ) 353 if (ga_nodeid().eq.0) write(6,933) detci_absync_etime 354 933 format(/,23x,'o',5('<'),' (ab) max sync:',f10.2) 355 endif 356c 357c Free scratch memory 358c 359 if (.not.ma_pop_stack(l_detci_sprime)) 360 $ call errquit('detci: cannot pop stack',0, MA_ERR) 361 if (.not.ma_pop_stack(l_detci_cprime)) 362 $ call errquit('detci: cannot pop stack',0, MA_ERR) 363 if (.not.ma_pop_stack(l_detci_iscr)) 364 $ call errquit('detci: cannot pop stack',0, MA_ERR) 365 if (.not.ma_pop_stack(l_detci_lhsscr)) 366 $ call errquit('detci: cannot pop stack',0, MA_ERR) 367 if (.not.ma_pop_stack(l_detci_rhsscr)) 368 $ call errquit('detci: cannot pop stack',0, MA_ERR) 369 if (.not.ma_pop_stack(l_detci_fscr)) 370 $ call errquit('detci: cannot pop stack',0, MA_ERR) 371c 372c Free other internal memory blocks 373c 374 if (.not.ma_pop_stack(l_detci_g)) 375 $ call errquit('detci: cannot pop stack',0, MA_ERR) 376 if (.not.ma_pop_stack(l_detci_h)) 377 $ call errquit('detci: cannot pop stack',0, MA_ERR) 378 if (l_detci_exb.ne.CDETCI_INVALID) then 379 if (.not.ma_pop_stack(l_detci_exb)) 380 $ call errquit('detci: cannot pop stack',0, MA_ERR) 381 endif 382 if (l_detci_exa.ne.CDETCI_INVALID) then 383 if (.not.ma_pop_stack(l_detci_exa)) 384 $ call errquit('detci: cannot pop stack',0, MA_ERR) 385 endif 386c 387c 388c Invalidate common block 389c 390 cdetci_valid = 0 391c 392c 393 return 394 end 395 396 397 398 399 400c 401c Create a GA CI vector consistent 402c with previously initialized parameters 403c 404 logical function ga_civec_create( label, g_a ) 405 implicit none 406#include "errquit.fh" 407#include "mafdecls.fh" 408#include "global.fh" 409#include "detciP.fh" 410#include "detci.fh" 411#include "cdetci.fh" 412 character*(*) label 413 integer g_a 414c 415 integer nstra 416 integer nstrb 417c 418c Hardwired for now for 1024 processors 419c 420 integer map1(1024) 421c 422c Check if initialized 423c 424 if (cdetci_valid .ne. CDETCI_MAGIC) 425 $ call errquit('ga_civec_create called uninitialized',0, 426 & INPUT_ERR) 427c 428c Distribute both string_b and string_a such that different 429c symmetry blocks in string_b are not split between processors. 430c String_a will be split over nproc/nsym tasks 431c 432 nblock = 0 433 detstore = 0 434 ntasks = nproc/nsym 435 if (mod(nproc,nsym).ne.0) ntasks=ntasks+1 436 do isym=1,nsym 437 jsym=MULX(isym,ksym) 438 nblock=nblock+1 439 map1(nblock)=detstore+1 440 maxa=cdetci_ataba(norb,nela,jsym)-cdetci_ataba(0,0,jsym)+1 441 maxb=cdetci_atabb(norb,nelb,isym)-cdetci_atabb(0,0,isym)+1 442 do i=2,ntasks 443 nblock=nblock+1 444 map1(nblock)=detstore+(i-1)*(maxa/ntasks)+1 445 enddo 446 detstore=detstore+maxa*maxb 447 enddo 448c 449c Generate GA 450c 451 ga_civec_create = ga_create_irreg( MT_DBL, detstore, 1, label, 452 $ map1,nblock,1,1, g_a ) 453c 454 return 455 end 456 457 458 459 460 461 462c 463c Takes a list of configurations with coefficients 464c and generates an initial guess CI vector that 465c is both spin- and symmetry-adapted 466c Check to ensure initial guess has non-zero 467c component in requisite spin and symmetry state 468c 469 subroutine detci_guess( ngs, cfggs, cgs, g_civec, g_workvec ) 470 implicit none 471#include "errquit.fh" 472c 473#include "global.fh" 474#include "util.fh" 475#include "detciP.fh" 476#include "detci.fh" 477#include "cdetci.fh" 478#include "cdetcistats.fh" 479c 480 integer ngs 481 integer cfggs(*) 482 double precision cgs(ngs) 483 integer g_civec 484 integer g_workvec 485c 486 integer nstra 487 integer nstrb 488 integer rcfg(detci_maxelec_tot*detci_maxguess_cfg) 489 integer i, j, jj 490 double precision mxcnt, xx, tx 491c 492c 493c 494 if (cdetci_valid .ne. CDETCI_MAGIC) 495 $ call errquit('detci_guess with uninitialized internals',0, 496 & INPUT_ERR) 497 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 498 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 499c 500c Reorder guess indices to internal convention 501c 502 do j=1,ngs 503 jj = (j-1)*(cdetci_nela+cdetci_nelb) 504 do i=1,cdetci_nela+cdetci_nelb 505 rcfg(jj+i) = cdetci_irmap(cfggs(jj+i)) 506 enddo 507 call util_isort(cdetci_nela,rcfg(jj+1)) 508 call util_isort(cdetci_nelb,rcfg(jj+cdetci_nela+1)) 509 enddo 510c 511c Generate normalized guess CI vector 512c 513 call detci_ciguess( cdetci_norb, cdetci_nsym, cdetci_nela, 514 $ cdetci_nelb, nstra, nstrb, cdetci_osym, 515 $ cdetci_ataba, cdetci_atabb, 516 $ ngs, rcfg, cgs, g_civec ) 517c 518c Spin adapt 519c 520 if (cdetci_spinadapt) then 521 tx = util_cpusec() 522 call detci_spinadapt( g_civec, g_workvec ) 523 detci_spinadp_etime = detci_spinadp_etime + util_cpusec() - tx 524 endif 525c 526c Symmetry adapt 527c 528 if (cdetci_nsym.gt.1) then 529 tx = util_cpusec() 530 call detci_symmproject( cdetci_norb, cdetci_nsym, 531 $ cdetci_nela, cdetci_nelb, nstra, 532 $ nstrb, cdetci_osym, 533 $ cdetci_ataba, cdetci_atabb, 534 $ cdetci_symstate, .true., 535 $ mxcnt, g_civec ) 536 detci_symmadp_etime = detci_symmadp_etime + util_cpusec() - tx 537 endif 538c 539c Check for zero projected vector 540c 541 xx = ga_ddot(g_civec,g_civec) 542 if (xx.lt.0.0001d0) 543 $ call errquit('detci: initial guess has wrong spin/symmetry',0, 544 & INPUT_ERR) 545 return 546 end 547 548 549 550 551 552 553c 554c Easy interface to sigma vector product 555c 556c 557 subroutine detci_sigma( g_civec, g_sigma ) 558 implicit none 559#include "errquit.fh" 560c 561#include "mafdecls.fh" 562#include "global.fh" 563#include "util.fh" 564#include "detciP.fh" 565#include "detci.fh" 566#include "cdetci.fh" 567#include "cdetcistats.fh" 568c 569 integer g_civec 570 integer g_sigma 571c 572c$$$ double precision civec(*) 573c$$$ double precision sigma(*) 574c 575 double precision t0, tx 576 integer nexa, nexb 577 integer nstra, nstrb 578 integer nekla, neklb 579 integer ntij 580 integer ijmap(detci_maxorb*detci_maxorb) 581 integer g_ct, g_st 582 integer mseg 583c 584c Check internals 585c 586 if (cdetci_valid .ne. CDETCI_MAGIC) 587 $ call errquit('detci_sigma with uninitialized internals',0, 588 & INPUT_ERR) 589c 590c 591c 592 detci_sigma_calls = detci_sigma_calls + 1 593 t0 = util_cpusec() 594c 595c CI parameters 596c 597 call detci_ijmap( cdetci_norb, cdetci_nsym, cdetci_osym, 598 $ ntij, ijmap ) 599 nexa = (cdetci_norb-cdetci_nela+1)*cdetci_nela 600 nexb = (cdetci_norb-cdetci_nelb+1)*cdetci_nelb 601 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 602 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 603 nekla = detci_binomial( (cdetci_norb-1), (cdetci_nela-1) ) 604 neklb = detci_binomial( (cdetci_norb-1), (cdetci_nelb-1) ) 605 mseg = (10*ga_nnodes())/max(min(nexa,nexb),1) 606 mseg = min(mseg, 20) 607 MSEG = 1 ! Delete for better parallel effeciency 608c ! mseg = number of simultaneous F string contructed 609c 610c Sigma components 611c 612c beta-beta using transpose vectors 613c 614 tx = util_cpusec() 615*ga:1:0 616 if (.not.ga_create( MT_DBL, nstra, nstrb, 'transp CI', 617 $ nstra, 0, g_ct )) 618 $ call errquit('detci_sigma: cannot allocate transp ci',0, 619 & GA_ERR) 620*ga:1:0 621 if (.not.ga_create( MT_DBL, nstra, nstrb, 'transp sig', 622 $ nstra, 0, g_st )) 623 $ call errquit('detci_sigma: cannot allocate transp sig',0, 624 & GA_ERR) 625 call ga_transpose( g_civec, g_ct ) 626 call ga_zero(g_st) 627 628 call detci_sigmaaa( cdetci_norb, cdetci_nsym, cdetci_nelb, 629 $ cdetci_nela, nstrb, nstra, nexb, nexa, 630 $ neklb, nekla, cdetci_osym, 631 $ ijmap, int_mb(k_detci_exb), 632 $ int_mb(k_detci_exa), 633 $ cdetci_atabb, cdetci_ataba, ntij, mseg, 634 $ dbl_mb(k_detci_h), dbl_mb(k_detci_g), 635 $ dbl_mb(k_detci_fscr), 636 $ g_ct, g_st ) 637 638 call ga_transpose( g_st, g_sigma ) 639 if (.not.ga_destroy(g_ct)) 640 $ call errquit('detci_sigma: cannot destroy transp ci',0, 641 & GA_ERR) 642 if (.not.ga_destroy(g_st)) 643 $ call errquit('detci_sigma: cannot destroy transp sig',0, 644 & GA_ERR) 645 detci_sbb_etime = detci_sbb_etime + util_cpusec() - tx 646c 647c alpha-alpha 648c 649 tx = util_cpusec() 650 call detci_sigmaaa( cdetci_norb, cdetci_nsym, cdetci_nela, 651 $ cdetci_nelb, nstra, nstrb, nexa, nexb, 652 $ nekla, neklb, cdetci_osym, 653 $ ijmap, int_mb(k_detci_exa), 654 $ int_mb(k_detci_exb), 655 $ cdetci_ataba, cdetci_atabb, ntij, mseg, 656 $ dbl_mb(k_detci_h), dbl_mb(k_detci_g), 657 $ dbl_mb(k_detci_fscr), 658 $ g_civec, g_sigma ) 659 detci_saa_etime = detci_saa_etime + util_cpusec() - tx 660c 661c alpha-beta 662c 663 tx = util_cpusec() 664 call detci_sigmaab( cdetci_norb, cdetci_nsym, cdetci_nela, 665 $ cdetci_nelb, nstra, nstrb, nexa, nexb, 666 $ nekla, neklb, cdetci_osym, 667 $ ijmap, int_mb(k_detci_exa), 668 $ int_mb(k_detci_exb), 669 $ cdetci_ataba, cdetci_atabb, ntij, 670 $ dbl_mb(k_detci_g), 671 $ int_mb(k_detci_rhsscr), 672 $ int_mb(k_detci_lhsscr), 673 $ int_mb(k_detci_iscr), 674 $ dbl_mb(k_detci_fscr), 675 $ dbl_mb(k_detci_cprime), 676 $ dbl_mb(k_detci_sprime), 677 $ g_civec, g_sigma ) 678 detci_sab_etime = detci_sab_etime + util_cpusec() - tx 679 detci_sigma_etime = detci_sigma_etime + util_cpusec() - t0 680c 681c 682c$$$ CALL DETCI_CIVEC_PRINT( CDETCI_NORB, CDETCI_NSYM, CDETCI_NELA, CDETCI_NELB, 683c$$$ $ NSTRA, NSTRB, CDETCI_OSYM, 684c$$$ $ CDETCI_ATABA, CDETCI_ATABB, SIGMA, 1.D-5 ) 685c 686c 687c 688 return 689 end 690 691 692 693 694 695 696 697 698 subroutine detci_ciprecon( g_civec, g_workvec ) 699 implicit none 700#include "errquit.fh" 701#include "global.fh" 702#include "util.fh" 703#include "detciP.fh" 704#include "detci.fh" 705#include "cdetci.fh" 706#include "cdetcistats.fh" 707c 708 integer g_civec 709 integer g_workvec 710 integer nstra, nstrb 711c 712c$$$ double precision xx, yy 713 double precision tx 714 double precision mxcnt 715c 716c Check internals 717c 718 if (cdetci_valid .ne. CDETCI_MAGIC) 719 $ call errquit('detci_ciprecon with uninitialized common',0, 720 & INPUT_ERR) 721c 722c CI parameters 723c 724 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 725 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 726c 727c Apply preconditioning using orbital 728c energies (Moller-Plesset denominators) 729c 730 call detci_diagscale( cdetci_norb, cdetci_nsym, cdetci_nela, 731 $ cdetci_nelb, nstra, nstrb, cdetci_osym, 732 $ cdetci_ataba, cdetci_atabb, 733 $ cdetci_eref, cdetci_eps, g_civec ) 734c 735c Spin-adapt the vector 736c 737 if (cdetci_spinadapt) then 738c$$$ PRINT*,'=== ENTRY VECTOR ===' 739c$$$ CALL DETCI_PRINT(G_CIVEC,1.D-3) 740CCCCC CALL DETCI_RANDOM_ERROR( G_CIVEC ) ! INTRODUCE CONTAMINATION TO DEBUG 741c$$$ XX = SQRT( GA_DDOT( G_CIVEC, G_CIVEC ) ) 742 743 tx = util_cpusec() 744 call detci_spinadapt( g_civec, g_workvec ) 745 detci_spinadp_etime = detci_spinadp_etime + util_cpusec() - tx 746 747c$$$ YY = SQRT( GA_DDOT( G_CIVEC, G_CIVEC ) ) 748c$$$ WRITE(6,988) (1.d0 - YY/XX) 749c$$$ 988 FORMAT('Spin contamination:',E12.3) 750c$$$ PRINT*,'=== AFTER SPIN ADAPT ===' 751c$$$ CALL DETCI_PRINT(G_CIVEC,1.D-3) 752 endif 753c 754c Ensure vector is symmetry-adapted 755c 756c$$$ PRINT*,'=== AFTER PRECON ===' 757c$$$ CALL DETCI_PRINT(G_CIVEC,1.D-3) 758 if (cdetci_nsym.gt.1) then 759 tx = util_cpusec() 760 call detci_symmproject( cdetci_norb, cdetci_nsym, 761 $ cdetci_nela, cdetci_nelb, nstra, 762 $ nstrb, cdetci_osym, 763 $ cdetci_ataba, cdetci_atabb, 764 $ cdetci_symstate, .true., 765 $ mxcnt, g_civec ) 766 detci_symmadp_etime = detci_symmadp_etime + util_cpusec() - tx 767 endif 768 769c$$$ PRINT*,'=== AFTER SYMMPROJ ===' 770c$$$ CALL DETCI_PRINT(G_CIVEC,1.D-3) 771 772 return 773 end 774 775 776 777 778 779 780 781c 782c Print CI vector coefficients larger 783c than threshold 784c 785 subroutine detci_print( g_civec, thresh ) 786 implicit none 787#include "errquit.fh" 788c 789#include "global.fh" 790#include "detciP.fh" 791#include "detci.fh" 792#include "cdetci.fh" 793c 794 integer g_civec 795C double precision civec(*) 796 double precision thresh 797 integer nstra, nstrb 798c 799c Check internals 800c 801 if (cdetci_valid .ne. CDETCI_MAGIC) 802 $ call errquit('detci_print with uninitialized common',0, 803 & INPUT_ERR) 804c 805c CI parameters 806c 807 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 808 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 809c 810c 811c 812 call detci_civec_print( cdetci_norb, cdetci_nsym, cdetci_nela, 813 $ cdetci_nelb, nstra, nstrb, 814 $ cdetci_osym, cdetci_ixmap, 815 $ cdetci_ataba, cdetci_atabb, 816 $ g_civec, thresh ) 817c 818c 819c 820 return 821 end 822 823 824 825 826 827 828 829 830c 831c Generate 1- and 2-particle density matrices 832c for given CI vector (both RHS and LHS) 833c 834 subroutine detci_density( g_civec, onepdm, twopdm ) 835 implicit none 836#include "errquit.fh" 837c 838#include "global.fh" 839#include "mafdecls.fh" 840#include "util.fh" 841#include "detciP.fh" 842#include "detci.fh" 843#include "cdetci.fh" 844#include "cdetcistats.fh" 845c 846 integer g_civec 847 double precision onepdm(*) 848 double precision twopdm(*) 849c 850 integer nexa, nexb 851 integer nstra, nstrb 852 integer nekla, neklb 853 integer nn 854 integer g_civect 855 double precision tx, tx2 856c 857c Check internals 858c 859 if (cdetci_valid .ne. CDETCI_MAGIC) 860 $ call errquit('detci_density with uninitialized common',0, 861 & INPUT_ERR) 862c 863c CI parameters 864c 865 nexa = (cdetci_norb-cdetci_nela+1)*cdetci_nela 866 nexb = (cdetci_norb-cdetci_nelb+1)*cdetci_nelb 867 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 868 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 869 nekla = detci_binomial( (cdetci_norb-1),(cdetci_nela-1)) 870 neklb = detci_binomial( (cdetci_norb-1),(cdetci_nelb-1)) 871 nn = cdetci_norb*cdetci_norb 872c 873c 874c 875 tx = util_cpusec() 876*ga:1:0 877 if (.not.ga_create( MT_DBL, nstra, nstrb, 'civec transp', 878 $ nstra, 0, g_civect )) 879 $ call errquit('detci_density: cannot create CI transp.',0, 880 & GA_ERR) 881 call ga_transpose( g_civec, g_civect ) 882c 883c 884c 885 tx2 = util_cpusec() 886 call detci_onepdm( cdetci_norb, cdetci_nsym, cdetci_nela, 887 $ cdetci_nelb, nstra, nstrb, nexa, nexb, 888 $ cdetci_osym, int_mb(k_detci_exa), 889 $ int_mb(k_detci_exb), cdetci_ixmap, 890 $ g_civec, g_civect, onepdm ) 891 detci_density_onepdm = detci_density_onepdm + util_cpusec() - tx2 892c 893c 894c 895 tx2 = util_cpusec() 896 call dfill(nn*nn, 0.0d0, twopdm, 1) 897 call detci_twopdm( cdetci_norb, cdetci_nsym, cdetci_nela, 898 $ cdetci_nelb, nstra, nstrb, nexa, nexb, 899 $ cdetci_osym, int_mb(k_detci_exa), 900 $ int_mb(k_detci_exb), cdetci_ixmap, 901 $ g_civec, g_civect, twopdm ) 902 detci_density_twopdm = detci_density_twopdm + util_cpusec() - tx2 903c 904c 905c 906 tx2 = util_cpusec() 907 call detci_twopdm_ab( cdetci_norb, cdetci_nsym, cdetci_nela, 908 $ cdetci_nelb, nstra, nstrb, nexa, nexb, 909 $ nekla, neklb, cdetci_osym, cdetci_ixmap, 910 $ int_mb(k_detci_exa), int_mb(k_detci_exb), 911 $ cdetci_ataba, cdetci_atabb, 912 $ int_mb(k_detci_rhsscr), 913 $ int_mb(k_detci_lhsscr), 914 $ int_mb(k_detci_iscr), 915 $ dbl_mb(k_detci_cprime), 916 $ dbl_mb(k_detci_sprime), 917 $ g_civec, twopdm ) 918 detci_density_twopdmab=detci_density_twopdmab+util_cpusec()-tx2 919 call ga_dgop(1, twopdm, nn*nn, '+') 920c 921c 922c 923 if (.not.ga_destroy( g_civect )) 924 $ call errquit('detci_density: cannot destroy CI transp.',0, 925 & GA_ERR) 926 detci_density_etime = detci_density_etime + util_cpusec() - tx 927c 928c 929c 930 return 931 end 932 933 934 935 936 937c 938c Spin-adaption routine 939c Project off contaminants by Lowdin projection operator 940c S quantum number is determined by the highest component, M, 941c ((nela - nelb)/2) 942c 943 944 subroutine detci_spinadapt( g_civec, g_wvec ) 945 implicit none 946#include "errquit.fh" 947c 948#include "mafdecls.fh" 949#include "detciP.fh" 950#include "detci.fh" 951#include "cdetci.fh" 952#include "cdetcistats.fh" 953c 954c$$$ double precision civec(*) 955c$$$ double precision wvec(*) 956 integer g_civec 957 integer g_wvec 958c 959 integer nexa, nexb 960 integer nstra, nstrb 961 integer ntij 962 integer ijmap(detci_maxorb*detci_maxorb) 963c 964c Check internals 965c 966 if (cdetci_valid .ne. CDETCI_MAGIC) 967 $ call errquit('detci_spin with uninitialized internals',0, 968 & INPUT_ERR) 969c 970c CI parameters 971c 972 call detci_ijmap( cdetci_norb, cdetci_nsym, cdetci_osym, 973 $ ntij, ijmap ) 974 nexa = (cdetci_norb-cdetci_nela+1)*cdetci_nela 975 nexb = (cdetci_norb-cdetci_nelb+1)*cdetci_nelb 976 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 977 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 978c 979c 980c 981 call detci_spadpt( cdetci_norb, cdetci_nsym, cdetci_nela, 982 $ cdetci_nelb, nstra, nstrb, 983 $ cdetci_osym, cdetci_ataba, cdetci_atabb, 984 $ nexa, nexb, int_mb(k_detci_exa), 985 $ int_mb(k_detci_exb), g_civec, g_wvec ) 986c 987c 988 return 989 end 990 991 992 993 994 995 996 997 998 999 1000 subroutine detci_symmproj( g_civec, mxcnt, oscreen ) 1001 implicit none 1002#include "errquit.fh" 1003#include "global.fh" 1004#include "detciP.fh" 1005#include "detci.fh" 1006#include "cdetci.fh" 1007c 1008 integer g_civec 1009 logical oscreen 1010 double precision mxcnt 1011 integer nstra, nstrb 1012c 1013c Check internals 1014c 1015 if (cdetci_valid .ne. CDETCI_MAGIC) 1016 $ call errquit('detci_symmproj with uninitialized common',0, 1017 & INPUT_ERR) 1018c 1019c CI parameters 1020c 1021 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 1022 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 1023c 1024c Symmetry adapt 1025c 1026 mxcnt = 0.d0 1027 if (cdetci_nsym.gt.1) then 1028 call detci_symmproject( cdetci_norb, cdetci_nsym, 1029 $ cdetci_nela, cdetci_nelb, nstra, 1030 $ nstrb, cdetci_osym, 1031 $ cdetci_ataba, cdetci_atabb, 1032 $ cdetci_symstate, oscreen, 1033 $ mxcnt, g_civec ) 1034 endif 1035 1036 return 1037 end 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 subroutine detci_moint_copy( norb, ntij, irmap, ijmap, orbsym, 1048 $ h, g, h1, g1 ) 1049 implicit none 1050#include "errquit.fh" 1051#include "global.fh" 1052 integer norb 1053 integer ntij 1054 integer irmap(norb) 1055 integer ijmap(norb,norb) 1056 integer orbsym(norb) 1057 double precision h(ntij) 1058 double precision g(ntij,ntij) 1059 double precision h1(ntij) 1060 double precision g1(ntij,ntij) 1061c 1062 integer i, j, k, l 1063 integer ij, rij, kl, rkl 1064 integer ijsym, ijklsym, nmixed 1065 double precision TOL 1066#include "symmdef.fh" 1067#include "bitops.fh" 1068#include "symmmul.fh" 1069 data TOL/1.d-6/ ! rjh - was 1d-9 but this is not met - why? 1070c 1071 nmixed = 0 1072 do i=1,norb 1073 do j=1,i 1074 ij = ijmap(i,j) 1075 rij = ijmap(irmap(i),irmap(j)) 1076 h1(ij) = h(rij) 1077 ijsym = MULX(orbsym(i),orbsym(j)) 1078 if (ijsym.ne.1) then 1079 if (abs(h1(ij)).gt.TOL) then 1080 write(6,*) ' H1 contaminated ', i, j, h1(ij) 1081 nmixed = nmixed + 1 1082 endif 1083 h1(ij) = 0.0d0 1084 endif 1085 do k=1,norb 1086 do l=1,k 1087 kl = ijmap(k,l) 1088 rkl = ijmap(irmap(k),irmap(l)) 1089 g1(ij,kl) = g(rij,rkl) 1090 ijklsym = MULX(ijsym,MULX(orbsym(k),orbsym(l))) 1091 if (ijklsym.ne.1) then 1092 if (abs(g1(ij,kl)).gt.TOL) then 1093 nmixed = nmixed + 1 1094 write(6,*) ' G1 contaminated ', i,j,k,l,g1(ij,kl) 1095 endif 1096 g1(ij,kl) = 0.0d0 1097 endif 1098 enddo 1099 enddo 1100 enddo 1101 enddo 1102 if ((nmixed.gt.0).and.(ga_nodeid().eq.0)) then 1103 call errquit('DETCI: MO ints symmetry contaminated',nmixed, 1104 & GEOM_ERR) 1105 endif 1106c 1107 end 1108 1109