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 "stdio.fh" 17#include "global.fh" 18#include "util.fh" 19#include "detciP.fh" 20#include "detci.fh" 21#include "cdetci.fh" 22#include "cdetcistats.fh" 23c 24 integer norb 25 integer nela 26 integer nelb 27 integer nsym 28 integer symstate 29 integer osym(norb) 30 integer iprint 31 double precision eps(norb) 32 double precision h(*) 33 double precision g(*) 34c 35c 36c 37 double precision tx 38 integer nexa, nexb 39 integer nstra, nstrb 40 integer nekla, neklb 41 integer ntij, civlen 42 integer ijmap(detci_maxorb*detci_maxorb) 43 integer vtab((detci_maxorb+1)*(detci_maxelec+1)*(detci_maxsy)) 44 integer i, j, ii 45c$$$ integer ij 46 integer memuse03 47 integer mseg 48c 49c Reset statistics 50c 51 if (norb.gt.detci_maxorb) then 52 if (ga_nodeid().eq.0) then 53 write(LuOut,*)'*** detci_init: detci_maxorb = ',detci_maxorb 54 write(LuOut,*)'*** detci_init: norb = ',norb 55 write(LuOut,*)'*** maximum no. orbitals exceeded, edit ' 56 + //'detciP.fh and recompile' 57 endif 58 call errquit("detci_init: maximum no. orbitals exceeded",0, 59 + UERR) 60 endif 61 if (nela.gt.detci_maxelec.or.nelb.gt.detci_maxelec) then 62 if (ga_nodeid().eq.0) then 63 write(LuOut,*)'*** detci_init: detci_maxelec = ',detci_maxelec 64 write(LuOut,*)'*** detci_init: nela = ',nela 65 write(LuOut,*)'*** detci_init: nelb = ',nelb 66 write(LuOut,*)'*** maximum no. electrons exceeded, edit ' 67 + //'detciP.fh and recompile' 68 endif 69 call errquit("detci_init: maximum no. electrons exceeded",0, 70 + UERR) 71 endif 72 tx = util_cpusec() 73 detci_init_etime = 0.d0 74 detci_saa_etime = 0.d0 75 detci_sbb_etime = 0.d0 76 detci_sab_etime = 0.d0 77 detci_sigma_etime = 0.d0 78 detci_density_etime = 0.d0 79 detci_sigma_calls = 0 80 detci_density_calls = 0 81 detci_spinadp_etime = 0.d0 82 detci_symmadp_etime = 0.d0 83 detci_aaff_etime = 0.d0 84 detci_aadot_etime = 0.d0 85 detci_aagop_etime = 0.d0 86 detci_abstr_etime = 0.d0 87 detci_abgath_etime = 0.d0 88 detci_abdotab_etime = 0.d0 89 detci_abscat_etime = 0.d0 90 detci_absync_etime = 0.d0 91 detci_density_onepdm = 0.0d0 92 detci_density_twopdm = 0.0d0 93 detci_density_twopdmab = 0.0d0 94 cdetci_profprint = .true. 95c 96c Copy CI parameters 97c 98 cdetci_norb = norb 99 cdetci_nela = nela 100 cdetci_nelb = nelb 101 cdetci_nsym = nsym 102 cdetci_symstate = symstate 103c 104c Reorder orbital indices by symmetry 105c 106 ii = 0 107 do j=1,nsym 108 do i=1,norb 109 if (osym(i).eq.j) then 110 ii = ii + 1 111 cdetci_ixmap(ii) = i 112 cdetci_irmap(i) = ii 113 cdetci_osym(ii) = osym(i) 114 cdetci_eps(ii) = eps(i) 115 endif 116 enddo 117 enddo 118c 119c$$$ DO I=1,NORB 120c$$$ WRITE(6,'(I5,A5,I5,5X,F10.4)') I,' ->',CDETCI_IRMAP(I), 121c$$$ $ CDETCI_EPS(CDETCI_IRMAP(I)) 122c$$$ ENDDO 123c 124c Compute CI parameters 125c 126 call detci_ijmap( norb, nsym, cdetci_osym, ntij, ijmap ) 127 nexa = (norb-nela+1)*nela 128 nexb = (norb-nelb+1)*nelb 129 nstra = detci_binomial(norb,nela) 130 nstrb = detci_binomial(norb,nelb) 131 nekla = detci_binomial((norb-1),(nela-1)) 132 neklb = detci_binomial((norb-1),(nelb-1)) 133 civlen = nstra*nstrb 134 mseg = (10*ga_nnodes())/max(min(nexa,nexb),1) 135 mseg = min(mseg, 20) 136 MSEG = 1 ! Delete for better parallel effeciency 137c 138c Default reference energy is Aufbau ordering 139c ...although maybe reset by CI-guess routine! 140c 141 cdetci_eref = detci_refenergy( norb, nela, nelb, eps ) 142c 143c Toggle spin-adaption 144c Need some way of automatically turning this on and off 145c depending requested spin-state 146c 147 cdetci_spinadapt = .true. 148 cdetci_squant = (nela - nelb)/2.d0 149c 150c Info print 151c 152 if ((iprint.gt.0).and.(ga_nodeid().eq.0)) then 153 write(6,902) norb, nsym, civlen, 154 $ cdetci_symstate, 155 $ cdetci_spinadapt, 156 $ cdetci_squant, 157 $ nela, nelb, nstra, nstrb, 158 $ nexa, nexb, nekla, neklb 159 902 format(10x,'Active shells:',23x,i5,/, 160 $ 10x,'Irreps:',30x,i5,/, 161 $ 10x,'CI vector length:',5x,i20,/, 162 $ 10x,'State symmetry:',25x,i2,/, 163 $ 10x,'Spin adaption:',27x,l1,/, 164 $ 10x,'S quantum number:',15x,f10.3,//, 165 $ 37x,'Alpha',6x,'Beta',/, 166 $ 35x,2(7('-'),3x),/, 167 $ 10x,'Electrons:',12x,2i10,/, 168 $ 10x,'Strings:',14x,2i10,/, 169 $ 10x,'E_ij per string:',6x,2i10,/, 170 $ 10x,'Strings per E_ij:',5x,2i10,/) 171 call util_flush(6) 172 endif 173c 174c Construct arc weight tables 175c 176 call detci_vatable( norb, nela, nsym, cdetci_osym, 177 $ vtab, cdetci_ataba ) 178 call detci_vatable( norb, nelb, nsym, cdetci_osym, 179 $ vtab, cdetci_atabb ) 180c 181c Allocate and construct excitation operator table 182c 183 l_detci_exa = CDETCI_INVALID 184 l_detci_exb = CDETCI_INVALID 185 if (nexa.gt.0) then 186 if (.not.ma_push_get(MT_INT, (6*nexa*(nstra+1)), 'detci:exa', 187 $ l_detci_exa, k_detci_exa)) 188 $ call errquit('detci: cannot allocate',0, MA_ERR) 189 call detci_excit( norb, nela, nsym, nstra, nexa, cdetci_osym, 190 $ ijmap, cdetci_ataba, int_mb(k_detci_exa) ) 191 endif 192 if (nexb.gt.0) then 193 if (.not.ma_push_get(MT_INT, (6*nexb*(nstrb+1)), 'detci:exb', 194 $ l_detci_exb, k_detci_exb)) 195 $ call errquit('detci: cannot allocate',0, MA_ERR) 196 call detci_excit( norb, nelb, nsym, nstrb, nexb, cdetci_osym, 197 $ ijmap, cdetci_atabb, int_mb(k_detci_exb) ) 198 endif 199c$$$ IF (GA_NODEID().EQ.0) THEN 200c$$$ IJ = 6*(NEXA*(NSTRA+1) + NEXB*(NSTRB+1)) 201c$$$ I = MA_INQUIRE_AVAIL( MT_DBL ) 202c$$$ WRITE(6,913) IJ, I 203c$$$ 913 FORMAT('ALLOCATED AND GENERATED EXCITATION TABLE',/, 204c$$$ $ 'Requested ',I10,' Free',I10,' words') 205c$$$ CALL UTIL_FLUSH(6) 206c$$$ ENDIF 207c$$$ CALL GA_SYNC() 208c 209c Allocate ERI block 210c Copy integrals into internal block 211c 212 if (.not.ma_push_get(MT_DBL,ntij,'detci: h1', 213 $ l_detci_h, k_detci_h)) 214 $ call errquit('detci: cannot allocate',0, MA_ERR) 215 if (.not.ma_push_get( MT_DBL,(ntij*ntij),'detci: eri', 216 $ l_detci_g, k_detci_g)) 217 $ call errquit('detci: cannot allocate',0, MA_ERR) 218 call detci_moint_copy( norb, ntij, cdetci_ixmap, ijmap, 219 $ cdetci_osym, h, g, 220 $ dbl_mb(k_detci_h), dbl_mb(k_detci_g) ) 221c 222c Memory for scratch space should be allocated HERE! 223c 224 memuse03 = 2*max(nstra,nstrb)*mseg + 3*nekla + 4*nstrb*nekla 225 if (.not.ma_push_get(MT_DBL, (max(max(nstra,nstrb)*mseg,nexb)), 226 $ 'detci: fscr', l_detci_fscr, k_detci_fscr)) 227 $ call errquit('detci: cannot allocate fscr',0, MA_ERR) 228 229 if (.not.ma_push_get(MT_INT, nekla, 230 $ 'detci: rhsscr', l_detci_rhsscr, k_detci_rhsscr)) 231 $ call errquit('detci: cannot allocate rhs',0, MA_ERR) 232 233 if (.not.ma_push_get(MT_INT, max(nekla,neklb,nexb,nexa), 234 $ 'detci: lhsscr', l_detci_lhsscr, k_detci_lhsscr)) 235 $ call errquit('detci: cannot allocate lhs',0, MA_ERR) 236 237 if (.not.ma_push_get(MT_INT, nekla, 238 $ 'detci: iscr', l_detci_iscr, k_detci_iscr)) 239 $ call errquit('detci: cannot allocate iscr',0, MA_ERR) 240 241 if (.not.ma_push_get(MT_DBL, (nstrb*nekla), 242 $ 'detci: cprime', l_detci_cprime, k_detci_cprime)) 243 $ call errquit('detci: cannot allocate cprime',0, MA_ERR) 244 245 if (.not.ma_push_get(MT_DBL, (nstrb*nekla), 246 $ 'detci: sprime', l_detci_sprime, k_detci_sprime)) 247 $ call errquit('detci: cannot allocate sprime',0, MA_ERR) 248c 249c$$$ IF (GA_NODEID().EQ.0) THEN 250c$$$ IJ = MAX(NSTRA,NSTRB) + 3*NEKLA + 2*(NSTRB*NEKLA) 251c$$$ I = MA_INQUIRE_AVAIL( MT_DBL ) 252c$$$ WRITE(6,914) IJ, I 253c$$$ 914 FORMAT('ALLOCATED SCRATCH ARRAYS',/, 254c$$$ $ 'Requested ',I10,' Free',I10,' words') 255c$$$ CALL UTIL_FLUSH(6) 256c$$$ ENDIF 257c$$$ CALL GA_SYNC() 258c 259c Validate internals 260c 261 cdetci_valid = CDETCI_MAGIC 262 detci_init_etime = util_cpusec() - tx 263 return 264 end 265 266 267 268 269 270 271c 272c Compute internal memory required & CI-vector length 273c 274 275 integer function detci_memsiz( norb, nela, nelb, nsym, 276 $ osym, vlen ) 277 implicit none 278 integer norb 279 integer nela 280 integer nelb 281 integer nsym 282 integer osym(norb) 283 integer vlen 284c 285 integer nexa 286 integer nexb 287 integer nstra 288 integer nstrb 289 integer nekla 290 integer neklb 291 integer detci_binomial 292 external detci_binomial 293c 294c 295c 296 nexa = (norb-nela+1)*nela 297 nexb = (norb-nelb+1)*nelb 298 nstra = detci_binomial(norb,nela) 299 nstrb = detci_binomial(norb,nelb) 300 nekla = detci_binomial((norb-1),(nela-1)) 301 neklb = detci_binomial((norb-1),(nelb-1)) 302 303 vlen = nstra*nstrb 304 detci_memsiz = 0 305 return 306 end 307 308 309 310 311 312 subroutine detci_free() 313 implicit none 314#include "errquit.fh" 315c 316#include "mafdecls.fh" 317#include "global.fh" 318#include "msgids.fh" 319#include "detciP.fh" 320#include "detci.fh" 321#include "cdetci.fh" 322#include "cdetcistats.fh" 323c 324c 325c Check if initialized 326c 327 if (cdetci_valid .ne. CDETCI_MAGIC) 328 $ call errquit('detci_free with uninitialized internals',0, 329 & INPUT_ERR) 330c 331c Print stats 332c 333 if (cdetci_profprint) then 334 if (ga_nodeid().eq.0) write(6,901) detci_sigma_calls, 335 $ detci_saa_etime, 336 $ detci_sbb_etime, 337 $ detci_sab_etime, 338 $ detci_sigma_etime, 339 $ detci_aaff_etime, 340 $ detci_aagop_etime, 341 $ detci_aadot_etime, 342 $ detci_abstr_etime, 343 $ detci_abgath_etime, 344 $ detci_abdotab_etime, 345 $ detci_abscat_etime, 346 $ detci_absync_etime, 347 $ detci_density_etime, 348 $ detci_density_onepdm, 349 $ detci_density_twopdm, 350 $ detci_density_twopdmab, 351 $ detci_spinadp_etime, 352 $ detci_symmadp_etime 353 901 format(/,10x,'Number of sigma calls:',4x,i5, 354 $ /,23x,'o',5('<'),' (aa):',7x,f10.2, 355 $ /,23x,'o',5('<'),' (bb):',7x,f10.2, 356 $ /,23x,'o',5('<'),' (ab):',7x,f10.2, 357 $ /,23x,'o',5('<'),' (total)',5x,f10.2, 358 $ /,23x,'o',5('<'),' (aa) ff',5x,f10.2, 359 $ /,23x,'o',5('<'),' (aa) gop',4x,f10.2, 360 $ /,23x,'o',5('<'),' (aa) dot',4x,f10.2, 361 $ /,23x,'o',5('<'),' (ab) str',4x,f10.2, 362 $ /,23x,'o',5('<'),' (ab) gath',3x,f10.2, 363 $ /,23x,'o',5('<'),' (ab) dotab',2x,f10.2, 364 $ /,23x,'o',5('<'),' (ab) scat',3x,f10.2, 365 $ /,23x,'o',5('<'),' (ab) sync',3x,f10.2, 366 $ /,23x,'o',5('<'),' Density',5x,f10.2, 367 $ /,23x,'o',5('<'),' Density one',1x,f10.2, 368 $ /,23x,'o',5('<'),' Density two',1x,f10.2, 369 $ /,23x,'o',5('<'),' Density ab',2x,f10.2, 370 $ /,23x,'o',5('<'),' Spin adapt',2x,f10.2, 371 $ /,23x,'o',5('<'),' Symm adapt',2x,f10.2) 372 373 call ga_dgop(msg_detci_maxsync, detci_absync_etime, 374 $ 1, 'max' ) 375 if (ga_nodeid().eq.0) write(6,933) detci_absync_etime 376 933 format(/,23x,'o',5('<'),' (ab) max sync:',f10.2) 377 endif 378c 379c Free scratch memory 380c 381 if (.not.ma_pop_stack(l_detci_sprime)) 382 $ call errquit('detci: cannot pop stack',0, MA_ERR) 383 if (.not.ma_pop_stack(l_detci_cprime)) 384 $ call errquit('detci: cannot pop stack',0, MA_ERR) 385 if (.not.ma_pop_stack(l_detci_iscr)) 386 $ call errquit('detci: cannot pop stack',0, MA_ERR) 387 if (.not.ma_pop_stack(l_detci_lhsscr)) 388 $ call errquit('detci: cannot pop stack',0, MA_ERR) 389 if (.not.ma_pop_stack(l_detci_rhsscr)) 390 $ call errquit('detci: cannot pop stack',0, MA_ERR) 391 if (.not.ma_pop_stack(l_detci_fscr)) 392 $ call errquit('detci: cannot pop stack',0, MA_ERR) 393c 394c Free other internal memory blocks 395c 396 if (.not.ma_pop_stack(l_detci_g)) 397 $ call errquit('detci: cannot pop stack',0, MA_ERR) 398 if (.not.ma_pop_stack(l_detci_h)) 399 $ call errquit('detci: cannot pop stack',0, MA_ERR) 400 if (l_detci_exb.ne.CDETCI_INVALID) then 401 if (.not.ma_pop_stack(l_detci_exb)) 402 $ call errquit('detci: cannot pop stack',0, MA_ERR) 403 endif 404 if (l_detci_exa.ne.CDETCI_INVALID) then 405 if (.not.ma_pop_stack(l_detci_exa)) 406 $ call errquit('detci: cannot pop stack',0, MA_ERR) 407 endif 408c 409c 410c Invalidate common block 411c 412 cdetci_valid = 0 413c 414c 415 return 416 end 417 418 419 420 421 422c 423c Create a GA CI vector consistent 424c with previously initialized parameters 425c 426 logical function ga_civec_create( label, g_a ) 427 implicit none 428#include "errquit.fh" 429#include "mafdecls.fh" 430#include "global.fh" 431#include "detciP.fh" 432#include "detci.fh" 433#include "cdetci.fh" 434 character*(*) label 435 integer g_a 436c 437 integer nstra 438 integer nstrb 439c 440c Check if initialized 441c 442 if (cdetci_valid .ne. CDETCI_MAGIC) 443 $ call errquit('ga_civec_create called uninitialized',0, 444 & INPUT_ERR) 445c 446c 447c 448 nstra = detci_binomial( cdetci_norb, cdetci_nela ) 449 nstrb = detci_binomial( cdetci_norb, cdetci_nelb ) 450*ga:1:0 451 ga_civec_create = ga_create( MT_DBL, nstrb, nstra, label, 452 $ nstrb, 0, 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