1 subroutine hv_write(lun,ix,len,x) 2* 3* $Id$ 4* 5 implicit none 6 integer lun, ix, len 7 double precision x(len) 8 integer i, k, khi 9 10 write(lun,901) ix,(x(i),i=1,8) 11 901 format(i5,8f12.6) 12 do k=9,len,8 13 khi = min((k+7),len) 14 write(lun,902) (x(i),i=k,khi) 15 902 format(5x,8f12.6) 16 enddo 17 return 18 end 19 20 21 22 23 subroutine hv_reorder( nbf, nclosed, nact, g_x, x, y ) 24 implicit none 25 integer nbf, nclosed, nact 26 integer g_x 27 double precision x(*), y(*) 28c 29 integer nvir, vlen, voff, aoff, aend, xoff, xend, ii 30 integer nca, i, j 31c 32c 33 nvir = nbf - nclosed - nact 34 vlen = (nclosed+nact)*nvir + nclosed*nact 35 voff = nclosed + nact + 1 36 aoff = nclosed + 1 37 aend = nclosed + nact 38 nca = nclosed + nact 39c 40c 41 xoff = nvir*(nclosed+nact) + 1 42 xend = vlen 43 call ga_get(g_x,xoff,xend,1,1,x,(nclosed*nact)) 44 call dcopy((nclosed*nact),x,1,y,1) 45 46 xoff = 1 47 xend = nvir*nclosed 48 call ga_get(g_x,xoff,xend,1,1,x,(nclosed*nvir)) 49 ii = nact*nclosed 50 do i=1,nvir 51 do j=1,nclosed 52 y(ii+(i-1)*nca+j) = x((j-1)*nvir+i) 53 enddo 54 enddo 55 56 xoff = xend + 1 57 xend = nvir*(nclosed+nact) 58 call ga_get(g_x,xoff,xend,1,1,x,(nact*nvir)) 59 do i=1,nvir 60 do j=1,nact 61 y(ii+(i-1)*nca+nclosed+j) = x((j-1)*nvir+i) 62 enddo 63 enddo 64 65 call dscal(vlen,0.5d0,y,1) 66c 67c 68 return 69 end 70 71 72 73c$$$ 74c$$$ subroutine mattr(n,m,x,y) 75c$$$ implicit none 76c$$$ integer n,m 77c$$$ double precision x(n,m),y(m,n) 78c$$$ integer i,j 79c$$$ 80c$$$ do i=1,n 81c$$$ do j=1,m 82c$$$ y(j,i) = x(i,j) 83c$$$ enddo 84c$$$ enddo 85c$$$ return 86c$$$ end 87c$$$ 88 89 90 91 92 subroutine hv_writev(lun,len,x) 93 implicit none 94 integer lun, len 95 double precision x(len) 96 97 write(lun) x 98 return 99 end 100 101 102 103 subroutine hv_readv(lun,len,x) 104 implicit none 105 integer lun, len 106 double precision x(len) 107 108 read(lun) x 109 return 110 end 111 112 113 114 115 subroutine ga_rowprint( text, g_a ) 116 implicit none 117#include "errquit.fh" 118#include "mafdecls.fh" 119#include "global.fh" 120 character*(*) text 121 integer g_a 122 integer gtype,vlen,rlen 123 integer l_tmp, k_tmp 124 integer i 125 126 if (ga_nodeid().eq.0) then 127 call ga_inquire(g_a,gtype,vlen,rlen) 128 if (.not.ma_push_get(MT_DBL,vlen,'tmp',l_tmp,k_tmp)) 129 $ call errquit('ga_rowprint: cannot allocate',0, MA_ERR) 130 call ga_get(g_a,1,vlen,1,1,dbl_mb(k_tmp),1) 131 write(6,909) text 132 909 format(/,a) 133 write(6,910) (dbl_mb(k_tmp+i-1),i=1,vlen) 134 910 format(6F12.8) 135 if (.not.ma_pop_stack(l_tmp)) 136 $ call errquit('ga_rowprint: cannot pop stack',0, MA_ERR) 137 endif 138 call ga_sync() 139 return 140 end 141 142 143 144 145 146 subroutine mcscf_twopdm_print(n,dm2) 147 implicit none 148 integer n 149 double precision dm2(n,n,n,n) 150 double precision TOL 151 integer i,j,k,l,ltop 152 data TOL/2.d-1/ 153 154 do i=1,n 155 do j=1,n 156 do k=1,n 157 ltop = k 158 if (k.gt.j) ltop = j 159 do l=1,n 160 if (abs(dm2(i,j,k,l)).gt.TOL) then 161 write(6,771) i,j,k,l,dm2(i,j,k,l) 162 771 format(4i5,f18.10) 163 endif 164 enddo 165 enddo 166 enddo 167 enddo 168 return 169 end 170 171 172 173 174 175 subroutine ga_print_x( g_a ) 176 implicit none 177#include "global.fh" 178 integer g_a 179 integer type, n, m 180 integer dim 181 parameter(dim=40) 182 double precision xtmp(dim*dim) 183 integer chi, clo, i, j 184 185 call ga_inquire(g_a, type, n, m ) 186 if ((n*m).gt.(dim*dim)) return 187 call ga_get(g_a,1,n,1,n,xtmp,n) 188 chi = 0 189 11 clo = chi + 1 190 chi = min((clo + 7),m) 191 do j=1,n 192 write(6,771) (xtmp((i-1)*n + j),i=clo,chi) 193 771 format(8f12.6) 194 enddo 195 write(6,*) 196 if (chi.lt.n) goto 11 197 return 198 end 199 200 201 202 subroutine mcscf_trprint(n,x) 203 implicit none 204 integer n 205 double precision x(*) 206 integer i, ii, j, itop 207 208 do i=1,n 209 ii = (i*(i-1))/2 210 itop = min(10,i) 211 write(6,922) (x(ii+j),j=1,itop) 212 922 format(10f12.6) 213 enddo 214 return 215 end 216 217 218 219 220 221 222 223 224 225 226 227 228#ifdef MCSCF_DEBUGGER 229c 230c 231c This is the main routine for debugging 232c the MCSCF (this is not distributed) 233c 234c 235 subroutine mcscf_debugger( rtdb, basis, geom, nbf, nclosed, nact, 236 $ g_movecs, evals ) 237 implicit none 238#include "errquit.fh" 239#include "mafdecls.fh" 240#include "global.fh" 241#include "bas.fh" 242#include "geom.fh" 243#include "rtdb.fh" 244#include "util.fh" 245#include "sym.fh" 246#include "pstat.fh" 247#include "mcscfprof.fh" 248c 249 integer rtdb 250 integer geom, basis 251 integer nbf, nclosed, nact 252 integer g_movecs 253 double precision evals(*) 254c 255 integer nvir, noper, nsym 256 integer nactel, nela, nelb, multip, orlen 257 integer l_occ, k_occ, l_sym, k_sym 258 integer l_dm1, k_dm1, l_dm2, k_dm2 259 integer l_tmp, k_tmp 260 integer g_coul, g_exch 261 integer g_grad, g_prod, g_x 262 integer mo_lo, mo_hi 263 integer i, j, nmixed, clo, chi 264 integer blen 265 double precision pfac 266 double precision eone, etwo, energy, enrep, etrace 267 double precision ecore, eci, e0 268 double precision citol 269 double precision tol2e, gnorm, xx 270 logical oskel 271 logical ohalf, oblk 272c 273c--------------------- 274c Debugging variables 275c 276 integer mclosed, mopen ! ROHF occupation 277 integer g_fcv, g_fpv, g_fcp ! ROHF Fock matrices 278 integer g_afock, g_ifock, g_gfock 279 integer g_u, g_b, g_newgfock, g_tmp, g_tmp2 280 integer g_coul2, g_exch2, g_grad2 281 integer info 282 283 integer hdim 284 parameter(hdim=100) 285 double precision crap(1000), hh(hdim*hdim) 286 double precision scr(4*hdim), ev(hdim) 287 double precision edel4, edel, theta, gg1, gg2, edelF 288 double precision phi, xxstep, etgt, ezzz 289 290 double precision cjgtol 291 integer iii 292 double precision phimin,phimax,phiinc 293C data phimin,phimax,phiinc/-6.3d0,6.3d0,0.1d0/ 294 data phimin,phimax,phiinc/1.0d0,1.0d0,0.2d0/ 295 296c--------------------- 297c 298 integer ga_create_atom_blocked, ga_create_JKblocked 299 external ga_create_atom_blocked, ga_create_JKblocked 300 integer mcscf_rohf_den2occ 301 external mcscf_rohf_den2occ 302 double precision ga_trace_diag 303 external ga_trace_diag 304c 305 data ohalf/.true./ 306 data blen/16/ 307c 308c 309c 310c Get w.f. parameters 311c 312 nvir = nbf - nclosed - nact 313 orlen = (nclosed*nvir) + (nact*nvir) + (nclosed*nact) 314 nsym = sym_number_ops(geom)+1 315 if (.not.rtdb_get(rtdb,'scf:skeleton',MT_LOG,1,oskel)) 316 $ oskel = sym_number_ops(geom).gt.0 317 if (.not.geom_nuc_rep_energy( geom, enrep )) 318 $ call errquit('mcscf: cannot retrieve nuclear repulsion',0, 319 & GEOM_ERR) 320 if (.not.rtdb_get(rtdb,'mcscf:tol2e',MT_DBL,1,tol2e)) 321 $ tol2e = 1.d-12 ! Redundant recovered later 322c 323c Get electron and spin multiplicity (NB: for info only) 324c Active elec and multiplicity must be set --- no defaults 325c 326 if (.not.rtdb_get(rtdb,'mcscf:nactelec',MT_INT,1,nactel)) 327 $ call errquit('number of active electrons not set',0, 328 & RTDB_ERR) 329 if (.not.rtdb_get(rtdb,'mcscf:multiplicity',MT_INT,1,multip)) 330 $ call errquit('spin multiplicity not set',0, RTDB_ERR) 331 nela = (nactel + multip - 1)/2 332 nelb = nactel - nela 333 if ((mod((nactel + multip - 1),2).ne.0).or. 334 $ (nela.lt.0).or.(nelb.lt.0)) 335 $ call errquit('mcscf: incompatible elec and spin',0, 336 & INPUT_ERR) 337 MOPEN = 2 338 MCLOSED = (NACTEL - MOPEN)/2 339c 340c Print info 341c 342 write(6,900) 343 900 format(///,10x,40('='), 344 $ /,17x,'MCSCF Debug Section', 345 $ /,10x,40('='),//) 346 write(6,901) nbf, nclosed, nact, nactel, 347 $ multip, orlen 348 write(6,902) (nclosed*nvir),(nact*nvir), 349 $ (nclosed*nact) 350 901 format(18x,35('-'),/, 351 $ 20x,'Basis functions:',10x,i5,/, 352 $ 20x,'Inactive shells:',10x,i5,/, 353 $ 20x,'Active shells:',12x,i5,/, 354 $ 20x,'Active electrons:',9x,i5,/, 355 $ 20x,'Multiplicity:',13x,i5,/, 356 $ 20x,'Orbital rotations:',8x,i5) 357 902 format(25x,'Inact - Virt',9x,i5,/, 358 $ 25x,'Act - Virt',11x,i5,/, 359 $ 25x,'Inact - Act',10x,i5) 360 write(6,903) 361 903 format(18x,35('-')) 362c 363c Allocate Fock and gradient matrices 364c 365*ga:1:0 366 if (.not.ga_create(MT_DBL,nbf,nbf,'Act Fock',nbf,0,g_afock)) 367 $ call errquit('mcscf: cannot allocate active Fock',0, GA_ERR) 368*ga:1:0 369 if (.not.ga_create(MT_DBL,nbf,nbf,'In Fock',nbf,0,g_ifock)) 370 $ call errquit('mcscf: cannot allocate inactive Fock',0, 371 & GA_ERR) 372*ga:1:0 373 if (.not.ga_create(MT_DBL,nbf,nbf,'Gen Fock',nbf,0,g_gfock)) 374 $ call errquit('mcscf: cannot allocate general Fock',0, 375 & GA_ERR) 376*ga:1:0 377 if (.not.ga_create(MT_DBL,orlen,1,'Gradient',0,0,g_grad)) 378 $ call errquit('rohf_head: cannot allocate',0, GA_ERR) 379*ga:1:0 380 if (.not.ga_create(MT_DBL,orlen,1,'Product',0,0,g_prod)) 381 $ call errquit('rohf_head: cannot allocate',0, GA_ERR) 382*ga:1:0 383 if (.not.ga_create(MT_DBL,orlen,1,'Arg vec',0,0,g_x)) 384 $ call errquit('rohf_head: cannot allocate',0, GA_ERR) 385c 386c Create occupation vectors 387c 388 if (.not.ma_push_get(MT_DBL, nbf, 'MO occ', l_occ, k_occ)) 389 $ call errquit('mcscf: cannot allocate',0, MA_ERR) 390c 391c Allocate 1- & 2-PDM 392c 393 if (.not.ma_push_get(MT_DBL, (nact*nact*nact*nact), 394 $ '2P density', l_dm2, k_dm2)) 395 $ call errquit('mcscf: cannot allocate MO density',0, 396 & MA_ERR) 397 if (.not.ma_push_get(MT_DBL, (nact*nact), 398 $ '1P density', l_dm1, k_dm1)) 399 $ call errquit('mcscf: cannot allocate MO density',0, MA_ERR) 400c 401c Get orbital symmetries 402c 403 if (.not.ma_push_get(MT_INT, nbf, 'MO sym', l_sym, k_sym)) 404 $ call errquit('mcscf: cannot allocate symmetry',0, MA_ERR) 405 call sym_movecs_adapt( basis, 1.d-8, g_movecs, 406 $ int_mb(k_sym), nmixed ) 407 if (nmixed .ne. 0) call errquit( 408 $ 'mcscf: symmetry contamination in starting MOs', nmixed, 409 & GEOM_ERR) 410c 411c Print orbital info 412c 413 write(6,550) 414 550 format(/,2x,'Starting Orbital Energies') 415 write(6,551) (evals(i),i=1,nbf) 416 551 format(7f12.6) 417 write(6,887) 418 887 format(/,2x,'Orbital Symmetry Irreps') 419 write(6,888) (int_mb(k_sym+i),i=0,nbf-1) 420 888 format(16i3) 421c 422c Allocate operator matrices 423c Memory test required here! 424c 425 mo_lo = nclosed + 1 426 mo_hi = nclosed + nact 427 noper = (nact*(nact+1))/2 428 g_coul = ga_create_JKblocked(noper,nbf,nbf,'Coulomb Oper') 429 g_exch = ga_create_JKblocked(noper,nbf,nbf,'X Oper') 430 if (.not.rtdb_get(rtdb,'mcscf:aoblock',MT_LOG,1,oblk)) then 431 if (.not.rtdb_get(rtdb,'fourindex:aoblock',MT_LOG,1,oblk)) 432 $ oblk = .false. 433 endif 434c 435c Initial 4-Index Tranformation 436c 437 call moints_build_2x( basis, ohalf, oskel, 438 $ mo_lo, mo_lo, mo_hi, 1, nbf, 439 $ g_movecs, g_coul, .true., 440 $ g_exch, .true., blen, oblk ) 441c 442c Initial core energy 443c 444 call dfill((nact*nact),0.d0,dbl_mb(k_dm1),1) 445 call dfill((nact*nact*nact*nact),0.d0,dbl_mb(k_dm2),1) 446 call mcscf_etrace( geom, basis, nbf, nclosed, nact, 447 $ .false., oskel, tol2e, dbl_mb(k_dm1), 448 $ dbl_mb(k_dm2), g_movecs, g_coul, 449 $ eone, etwo, ecore ) 450c 451c 452c 453 cjgtol = 1.d-3 454 gnorm = 4.d0 455 e0 = ecore + enrep 456c 457c Solve CI to set density 458c 459 citol = 1.d-8 460 call mcscf_ifock( geom, basis, nbf, nclosed, nact, 461 $ oskel, tol2e, g_movecs, eone, etwo, 462 $ ecore, g_ifock ) 463 e0 = ecore + enrep 464 call mcscf_cisolve( rtdb, geom, basis, nbf, nclosed, nact, 465 $ nsym, int_mb(k_sym), e0, 466 $ evals, g_ifock, g_coul, 467 $ citol, .false., .true., .true., 468 $ dbl_mb(k_dm1), dbl_mb(k_dm2), eci ) 469c 470c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 471c 472c Natural Orbital Section 473c 474c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 475c 476c 477c$$$ call mcscf_naturalorb( nbf, nclosed, nact, dbl_mb(k_dm1), 478c$$$ $ dbl_mb(k_occ), evals, g_movecs ) 479c 480c Print natural orbitals 481c 482c$$$ write(6,880) 483c$$$ 880 format(///,10x,'Natural Orbitals and Occupation') 484c$$$ if (.not.ma_push_get(MT_DBL, (nbf*(nclosed+nact)), 'tmp', 485c$$$ $ l_tmp, k_tmp)) 486c$$$ $ call errquit('mcscf: cannot allocate local MO',0) 487c$$$ call ga_get(g_movecs, 1, nbf, 1, (nclosed+nact), 488c$$$ $ dbl_mb(k_tmp), nbf) 489c$$$ chi = 0 490c$$$ 33 clo = chi + 1 491c$$$ chi = min((clo + 7),(nact+nclosed)) 492c$$$ write(6,*) 493c$$$ write(6,881) (dbl_mb(k_occ+i-1),i=clo,chi) 494c$$$ 881 format(8f12.6) 495c$$$ write(6,*) 496c$$$ do i=1,nbf 497c$$$ write(6,881) (dbl_mb(k_tmp+(j-1)*nbf+i-1),j=clo,chi) 498c$$$ enddo 499c$$$ if (chi.ne.(nact+nclosed)) goto 33 500c$$$ if (.not.ma_pop_stack(l_tmp)) 501c$$$ $ call errquit('mcscf: cannot pop local MO',0) 502c 503c Resolve CI for Natural Orbitals 504c 505c$$$ write(6,228) 506c$$$ 228 format(/,'Resolving CI for natural orbitals',/) 507c$$$ call moints_build_2x( basis, ohalf, oskel, 508c$$$ $ mo_lo, mo_lo, mo_hi, 1, nbf, 509c$$$ $ g_movecs, g_coul, .true., 510c$$$ $ g_exch, .true., blen, oblk ) 511c$$$ call mcscf_fcore( basis, nbf, nclosed, nact, g_movecs, 512c$$$ $ g_coul, g_exch, g_ifock ) 513c$$$ call mcscf_cisolve( rtdb, geom, basis, nbf, nclosed, nact, 514c$$$ $ nsym, int_mb(k_sym), e0, evals, 515c$$$ $ g_ifock, g_coul, 516c$$$ $ citol, .false., .false., .false., 517c$$$ $ dbl_mb(k_dm1), dbl_mb(k_dm2), eci ) 518c 519c Print 1-pdm and 2-pdm 520c 521 if (util_print('density matrix',print_debug)) then 522 if (ga_nodeid().eq.0) then 523 write(6,671) 524 671 format(/,'<<<<<<< 1pdm density matrix >>>>>>>>>') 525 call moints_matprint( nact, nact, dbl_mb(k_dm1) ) 526 write(6,672) 527 672 format(/,'<<<<<<< 2pdm density matrix >>>>>>>>>') 528 call mcscf_twopdm_print(nact,dbl_mb(k_dm2)) 529 write(6,673) 530 673 format(/,'<<<<<<< symm. 2pdm density matrix >>>>>>>>>') 531 call mcscf_symmetrize_2pdm( nact, dbl_mb(k_dm2), crap ) 532 call mcscf_twopdm_print(nact,crap) 533 endif 534 endif 535c 536c 537c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 538c 539c MCSCF Energy Trace and Fock Section 540c 541c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 542c 543c Check energy trace 544c 545 call mcscf_etrace( geom, basis, nbf, nclosed, nact, 546 $ .true., oskel, tol2e, dbl_mb(k_dm1), 547 $ dbl_mb(k_dm2), g_movecs, g_coul, 548 $ eone, etwo, ecore ) 549 etrace = eone + etwo + enrep 550 write(6,674) etrace 551 674 format(20x,'Trace Energy: ',f20.14) 552c 553c Check Fock energy and gradient 554c 555 call ga_zero(g_grad) 556 call mcscf_fock( geom, basis, nbf, nclosed, nact, 557 $ oskel, tol2e, dbl_mb(k_dm1), dbl_mb(k_dm2), 558 $ g_movecs, g_coul, eone, etwo, e0, 559 $ g_ifock, g_afock, g_gfock ) 560 call mcscf_gfock2grad( nbf, nclosed, nact, g_gfock, g_grad ) 561 gnorm = sqrt(ga_ddot(g_grad,g_grad)) 562 energy = eone + etwo + enrep 563 write(6,675) energy 564 675 format(20x,'Fock Energy: ',f20.14) 565 write(6,676) gnorm 566 676 format(20x,'Gradient norm:',e20.8) 567c 568c 569c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 570c 571c MCSCF Hessian vector product Test Section 572c 573c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 574c 575c Regenerate MO integrals 576c 577 call moints_build_2x( basis, ohalf, oskel, 578 $ mo_lo, mo_lo, mo_hi, 1, nbf, 579 $ g_movecs, g_coul, .true., 580 $ g_exch, .true., blen, oblk ) 581c 582c 583 call mcscf_fcore( basis, nbf, nclosed, nact, g_movecs, 584 $ g_coul, g_exch, g_afock ) 585c 586c Finite difference gradient 587c 588c$$$ call ga_rowprint( ' ==== Analytical Gradient ====', g_grad ) 589c$$$ call mcscf_fdiff_grad( geom, basis, nbf, nclosed, nact, 590c$$$ $ oskel, tol2e, dbl_mb(k_dm1), 591c$$$ $ dbl_mb(k_dm2), g_movecs, g_coul, g_grad) 592c$$$ call ga_rowprint( ' ==== Finite Diff Gradient ====', g_grad ) 593c$$$ gnorm = sqrt(ga_ddot(g_grad,g_grad)) 594c$$$ write(6,773) gnorm 595c$$$ 773 format(/,10x,'Finite diff gradient norm:',e12.4,/) 596c 597c MCSCF Fock matrices and gradient 598c 599 call mcscf_fock( geom, basis, nbf, nclosed, nact, 600 $ oskel, tol2e, dbl_mb(k_dm1), dbl_mb(k_dm2), 601 $ g_movecs, g_coul, eone, etwo, e0, 602 $ g_ifock, g_afock, g_gfock ) 603 call mcscf_gfock2grad( nbf, nclosed, nact, g_gfock, g_grad ) 604c 605c 606c 607c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 608c 609c B Matrix Test Section 610c 611c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 612c 613c 614c Make G matrix on disk 615c (doesn't work..yet) 616c 617c$$$ call makeJK( basis, nbf, nclosed, nact, g_movecs ) 618c$$$ call mcscf_g( nbf, nclosed, nact, dbl_mb(k_dm1), dbl_mb(k_dm2), 619c$$$ $ g_ifock, g_afock ) 620c 621c 622c 623c$$$ if (.not.ga_duplicate(g_coul, g_coul2, 'Coul Copy')) 624c$$$ $ call errquit('mcscf: cannot allocate Coulomb copy',0) 625c$$$ if (.not.ga_duplicate(g_exch, g_exch2, 'Exch Copy')) 626c$$$ $ call errquit('mcscf: cannot allocate Exch copy',0) 627c$$$ if (.not.ga_duplicate(g_grad, g_grad2, 'Grad Copy')) 628c$$$ $ call errquit('mcscf: cannot allocate Exch copy',0) 629*ga:1:0 630c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'Unitary',nbf,0,g_u)) 631c$$$ $ call errquit('mcscf: cannot allocate unitrary',0) 632*ga:1:0 633c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'Temp',nbf,0,g_tmp)) 634c$$$ $ call errquit('mcscf: cannot allocate tmp',0) 635*ga:1:0 636c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'Temp 2',nbf,0,g_tmp2)) 637c$$$ $ call errquit('mcscf: cannot allocate tmp',0) 638*ga:1:0 639c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'B',nbf,0,g_b)) 640c$$$ $ call errquit('mcscf: cannot allocate B',0) 641*ga:1:0 642c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'new G Fock',nbf,0,g_newgfock)) 643c$$$ $ call errquit('mcscf: cannot allocate new GFock',0) 644c 645c Step and generate U = exp(K) (2x2 rotation) 646c 647c$$$ call ga_zero(g_b) 648c$$$ iii = (nclosed+nact-1)*nvir + 1 649c$$$ do phi=phimin,phimax,phiinc 650c$$$ xxstep = sqrt((1.d0 - cos(phi))) 651c$$$ call ga_copy( g_coul, g_coul2 ) 652c$$$ call ga_copy( g_exch, g_exch2 ) 653c$$$ call ga_copy( g_grad, g_grad2 ) 654c$$$ call ga_copy( g_grad2, g_x) 655c$$$ call ga_zero(g_u) 656c$$$ call theta2u( nbf, nclosed+nact, nclosed+nact+1, phi, g_u ) 657 658 659c$$$ PRINT*,'Unitrary transformation' 660c$$$ CALL GA_GET(G_U,1,NBF,1,NBF,CRAP,NBF) 661c$$$ CALL MOINTS_MATPRINT(NBF,NBF,CRAP) 662c 663c Generate B (Gradient response) 664c 665c$$$ call ga_zero(g_b) 666c$$$ call mcscf_b( geom, basis, nbf, nclosed, nact, int_mb(k_sym), 667c$$$ $ dbl_mb(k_dm1), dbl_mb(k_dm2), oskel, tol2e, 668c$$$ $ g_movecs, g_ifock, g_afock, g_coul2, g_exch2, 669c$$$ $ g_u, g_b ) 670c 671c ~ t 672c A = U B 673c 674c$$$ call ga_dgemm( 't', 'n', nbf, nbf, nbf, 1.d0, g_b, g_u, 675c$$$ $ 0.d0, g_newgfock ) 676c$$$ call mcscf_gfock2grad(nbf, nclosed, nact, g_newgfock,g_grad2) 677C CALL GA_ROWPRINT('Approximate New Grad',g_grad2) 678c$$$ gnorm = sqrt(ga_ddot(g_grad2,g_grad2)) 679c$$$ call ga_get(g_grad2,iii,iii,1,1,gg2,1) 680c 681c t 682c dE = tr( T (A + B) ) 683c 684c$$$ call ga_zero(g_tmp) 685c$$$ do i=1+ga_nodeid(),nbf,ga_nnodes() 686c$$$ call ga_put(g_tmp,i,i,i,i,1.d0,1) 687c$$$ enddo 688c$$$ call ga_dadd( 1.d0, g_u, -1.d0, g_tmp, g_tmp ) 689c$$$ edelF = ga_ddot( g_tmp, g_gfock ) 690c$$$ 691c$$$ call ga_transpose( g_b, g_tmp2 ) 692c$$$ call ga_dadd( 1.d0, g_gfock, 1.d0, g_tmp2, g_newgfock ) 693c$$$ call ga_transpose( g_newgfock, g_tmp2 ) 694c$$$ edel4 = ga_ddot( g_tmp, g_tmp2 ) 695c 696c G operator code is broken! 697c 698c (2) t ij 699c e = sum (T G T ) 700c ij ij 701c 702c (1) t 703c e = tr( T A ) 704c 705c$$$ call ga_get(g_tmp,1,nbf,1,nbf,crap,nbf) 706c$$$ call mcscf_gt( nbf, nclosed, nact, crap, etgt ) 707c$$$ ezzz = 2.d0*ga_ddot( g_tmp, g_gfock ) + etgt 708c 709c Compute exact new gradient for comparison 710c 711c$$$ call ga_dgemm( 'n', 'n', nbf, nbf, nbf, 1.d0, g_movecs, g_u, 712c$$$ $ 0.d0, g_newgfock ) 713c$$$ call moints_build_2x( basis, ohalf, oskel, 714c$$$ $ mo_lo, mo_lo, mo_hi, 1, nbf, 715c$$$ $ g_newgfock, g_coul2, .true., 716c$$$ $ g_exch2, .true., blen, oblk ) 717c$$$ call mcscf_fock( geom, basis, nbf, nclosed, nact, 718c$$$ $ oskel, tol2e, dbl_mb(k_dm1), dbl_mb(k_dm2), 719c$$$ $ g_newgfock, g_coul2, eone, etwo, e0, 720c$$$ $ g_u, g_tmp, g_tmp2 ) 721c$$$ call mcscf_gfock2grad( nbf, nclosed, nact, g_tmp2, g_grad2 ) 722c$$$ CALL GA_ROWPRINT('Exact New Grad',g_grad2) 723c$$$ gnorm = sqrt(ga_ddot(g_grad2,g_grad2)) 724c$$$ edel = eone + etwo + enrep - energy 725c$$$ call ga_get(g_grad2,iii,iii,1,1,gg1,1) 726c$$$ write(6,939) ((phi*180.d0)/3.1415629d0), 727c$$$ $ edel,edel4,ezzz,edelF,(edel4-edelF),etgt 728c$$$ 939 format('qqqq',f10.2,6e15.6) 729c$$$ enddo 730c 731c 732c 733c$$$ if (.not.ga_destroy(g_newgfock)) 734c$$$ $ call errquit('mcscf: cannot destroy unitrary',0) 735c$$$ if (.not.ga_destroy(g_u)) 736c$$$ $ call errquit('mcscf: cannot destroy unitrary',0) 737c$$$ if (.not.ga_destroy(g_tmp)) 738c$$$ $ call errquit('mcscf: cannot destroy unitrary',0) 739c$$$ if (.not.ga_destroy(g_tmp2)) 740c$$$ $ call errquit('mcscf: cannot destroy unitrary',0) 741c$$$ if (.not.ga_destroy(g_b)) 742c$$$ $ call errquit('mcscf: cannot destroy B',0) 743c$$$ if (.not.ga_destroy(g_coul2)) 744c$$$ $ call errquit('mcscf: cannot destroy coul copy',0) 745c$$$ if (.not.ga_destroy(g_exch2)) 746c$$$ $ call errquit('mcscf: cannot destroy exch copy',0) 747c$$$ if (.not.ga_destroy(g_grad2)) 748c$$$ $ call errquit('mcscf: cannot destroy grad copy',0) 749c 750c 751c 752c 753c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 754c 755c Explicit Hessian contruction and Eigenvalue Test Section 756c 757c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 758c 759c$$$ if (orlen.le.hdim) then 760c$$$ call makeJK( basis, nbf, nclosed, nact, g_movecs ) 761c$$$ call dfill((orlen*orlen),0.d0,hh,1) 762c$$$ call hmat( nbf, nclosed, nact, orlen, dbl_mb(k_dm1), 763c$$$ $ dbl_mb(k_dm2), g_ifock, g_afock, g_gfock, 764c$$$ $ g_coul, g_exch, hh) 765c$$$ 766c$$$ write(6,951) 767c$$$ 951 format('Hessian diagonal') 768c$$$ write(6,222) (hh((i-1)*orlen+i),i=1,orlen) 769c$$$ 222 format(8f16.11) 770c$$$ call dsyev('V','L',orlen,hh,orlen,ev,scr,(4*hdim),info) 771c$$$ write(6,953) 772c$$$ 953 format('Hessian eigenvalues') 773c$$$ write(6,222) (ev(i),i=1,orlen) 774c$$$ endif 775c 776c 777c Make explicit Hessian matrix from vector products 778c 779c$$$ call dfill((orlen*orlen),0.d0,hh,1) 780c$$$ call mcscf_hessmake( geom, basis, nbf, nclosed, nact, 781c$$$ $ oskel, orlen, g_movecs, dbl_mb(k_dm1), 782c$$$ $ dbl_mb(k_dm2), g_ifock, g_afock, 783c$$$ $ g_gfock, g_coul, g_exch, g_x, g_prod, 784c$$$ $ hh ) 785c 786c 787c Finite difference Hessian 788c 789c$$$ call mcscf_fdiff_hess( geom, basis, nbf, nclosed, nact, 790c$$$ $ oskel, tol2e, dbl_mb(k_dm1), 791c$$$ $ dbl_mb(k_dm2), g_movecs, g_coul, 792c$$$ $ g_grad ) 793c 794c Free temporaries for debugging section 795c 796c 797c 798c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 799c 800c ROHF Test Section 801c 802c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 803c 804c ROHF occupation and trivial density 805c 806c$$$ call mcscf_occ2int( nbf, dbl_mb(k_occ), mclosed, mopen ) 807c$$$ write(6,903) mclosed, mopen 808c$$$ 903 format(' ROHF Occupation',6x,'(closed):',i3,5x,'(open):',i3) 809c$$$ call mcscf_rohf_modens(mopen,nact,dbl_mb(k_dm1),dbl_mb(k_dm2)) 810c 811c 812c Fock build and 1e-Hessian vector product 813c 814*ga:1:0 815c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'CV Fock',nbf,0,g_fcv)) 816c$$$ $ call errquit('mcscf: cannot allocate active Fock',0) 817*ga:1:0 818c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'PV Fock',nbf,0,g_fpv)) 819c$$$ $ call errquit('mcscf: cannot allocate active Fock',0) 820*ga:1:0 821c$$$ if (.not.ga_create(MT_DBL,nbf,nbf,'CP Fock',nbf,0,g_fcp)) 822c$$$ $ call errquit('mcscf: cannot allocate active Fock',0) 823c 824c ROHF Gradient 825c 826c$$$ call rohf_fock( geom, basis, nclosed, nopen, tol2e, g_movecs, 827c$$$ $ eone, etwo, g_fcv, g_fpv, g_fcp, oskel ) 828c$$$ call rohf_fock2grad( nbf, nclosed, nact, 829c$$$ $ g_fcv, g_fpv, g_fcp, g_grad ) 830c$$$ call ga_rowprint( '==== ROHF Gradient ====', g_grad ) 831c 832c ROHF Hessian vector product 833c 834c$$$ call ga_copy(g_grad, g_x) 835c$$$ call ga_zero(g_prod) 836c$$$ pflg = 2 837c$$$ lshift = 0.d0 838c$$$ call rohf_hessv_xx( basis, geom, nbf, nclosed, nact, pflg, 839c$$$ $ g_movecs, oskel, g_fcv, g_fpv, g_fcp, 840c$$$ $ tol2e, lshift, g_x, g_prod ) 841c$$$ 842c$$$ call ga_rowprint( '==== ROHF Product', g_prod ) 843c 844c Make ROHF Hessian 845c 846c$$$ call rohf_hessmake( basis, geom, nbf, nclosed, nact, 847c$$$ $ g_movecs, oskel, g_fcv, g_fpv, g_fcp, 848c$$$ $ tol2e, g_x, g_prod ) 849c$$$ 850c 851c Create ROHF Hessian diagonal to compare with 852c 853c$$$ call rohf_hxxx( nbf, nclosed, nact, 0.d0, g_fcv, g_fpv, 854c$$$ $ g_fcp, g_prod ) 855c 856c Deallocate stuff 857c 858c$$$ if (.not.ga_destroy(g_fcv)) 859c$$$ $ call errquit('mcscf: cannot destroy MO vectors',0) 860c$$$ if (.not.ga_destroy(g_fpv)) 861c$$$ $ call errquit('mcscf: cannot destroy MO vectors',0) 862c$$$ if (.not.ga_destroy(g_fcp)) 863c$$$ $ call errquit('mcscf: cannot destroy MO vectors',0) 864c$$$ if (.not.ga_destroy(g_grad)) 865c$$$ $ call errquit('mcscf: cannot destroy gradient',0) 866c$$$ if (.not.ga_destroy(g_prod)) 867c$$$ $ call errquit('mcscf: cannot destroy product',0) 868c$$$ if (.not.ga_destroy(g_x)) 869c$$$ $ call errquit('mcscf: cannot destroy product',0) 870c 871c 872c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 873c 874c Final Cleanup 875c 876c %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 877c 878 401 continue 879c 880 if (.not.ma_pop_stack(l_sym)) 881 $ call errquit('mcscf: cannot pop stack?',0, MA_ERR) 882 if (.not.ma_pop_stack(l_dm1)) 883 $ call errquit('mcscf: cannot pop stack?',0, MA_ERR) 884 if (.not.ma_pop_stack(l_dm2)) 885 $ call errquit('mcscf: cannot pop stack?',0, MA_ERR) 886 if (.not.ma_pop_stack(l_occ)) 887 $ call errquit('mcscf: cannot pop stack?',0, MA_ERR) 888 if (.not.ga_destroy(g_exch)) 889 $ call errquit('mcscf: cannot destroy exchange',0, GA_ERR) 890 if (.not.ga_destroy(g_coul)) 891 $ call errquit('mcscf: cannot destroy Coulomb',0, GA_ERR) 892 893 if (.not.ga_destroy(g_grad)) 894 $ call errquit('rohf_head: cannot destroy gradient',0, GA_ERR) 895 if (.not.ga_destroy(g_prod)) 896 $ call errquit('rohf_head: cannot destroy product',0, GA_ERR) 897 if (.not.ga_destroy(g_x)) 898 $ call errquit('rohf_head: cannot destroy product',0, GA_ERR) 899 if (.not.ga_destroy(g_afock)) 900 $ call errquit('rohf_head: cannot destroy afock',0, GA_ERR) 901 if (.not.ga_destroy(g_ifock)) 902 $ call errquit('rohf_head: cannot destroy ifock',0, GA_ERR) 903 if (.not.ga_destroy(g_gfock)) 904 $ call errquit('rohf_head: cannot destroy gfock',0, GA_ERR) 905 906 write(6,771) 907 771 format(///,10x,40('='), 908 $ /,17x,'END DEBUG SECTION', 909 $ /,10x,40('='),//) 910c 911c 912c 913 return 914 end 915 916 917 918 919 920 921 922 923 924 subroutine mcscf_hessmake( geom, basis, nbf, nclosed, nact, 925 $ oskel, orlen, g_movecs, dm1, dm2, 926 $ g_ifock, g_afock, g_gfock, 927 $ g_coul, g_exch, 928 $ g_x, g_prod, hh ) 929 implicit none 930#include "errquit.fh" 931#include "mafdecls.fh" 932 integer geom, basis 933 integer nbf, nclosed, nact 934 logical oskel 935 integer orlen 936 double precision dm1(*), dm2(*) 937 integer g_movecs, g_ifock, g_afock, g_gfock 938 integer g_coul, g_exch 939 integer g_x, g_prod 940 double precision hh(orlen,orlen) 941c 942c 943C INTEGER L_HD, K_HD 944 integer l_hv, k_hv, l_hy, k_hy 945 integer nvir, incr, xoff, i, j, ii, pflg 946 double precision xx, lshift, tol2e 947 data pflg/2/ 948 data lshift/0.d0/ 949 data tol2e/1.d-12/ 950c 951c 952c 953 if (.not.ma_push_get(MT_DBL, orlen, 'H', l_hv, k_hv)) 954 $ call errquit('mcscf: cannot allocate',0, MA_ERR) 955 if (.not.ma_push_get(MT_DBL, orlen, 'Hy', l_hy, k_hy)) 956 $ call errquit('mcscf: cannot allocate',0, MA_ERR) 957 nvir = nbf - nclosed - nact 958C goto 10000 959 open(unit=88,file='hess',form='unformatted', 960 $ status='unknown') 961 write(88) nclosed,nact,nvir 962 write(88) (nclosed*nvir),((nclosed+nact)*nvir),orlen 963 xx = 1.d0 964 incr = 0 965c 966c 967c 968 do j=1,nclosed+nact 969 do i=1,nvir 970 ii = (j-1)*nvir + i 971 call ga_zero(g_x) 972 call ga_put(g_x,ii,ii,1,1,xx,1) 973 call mcscf_hessv( geom, basis, nbf, nclosed, nact, 974 $ tol2e, .false., pflg, lshift, dm1, dm2, 975 $ g_movecs, g_ifock, g_afock, g_gfock, 976 $ g_coul, g_exch, g_x, g_prod ) 977 call ga_get(g_prod,1,orlen,1,1,dbl_mb(k_hv),orlen) 978 call hv_writev(88,orlen,dbl_mb(k_hv)) 979 call dcopy(orlen,dbl_mb(k_hv),1,hh(1,ii),1) 980 enddo 981 enddo 982c 983 xoff = nvir*(nclosed+nact) 984 do j=1,nclosed 985 do i=1,nact 986 ii = xoff + (j-1)*nact + i 987 call ga_zero(g_x) 988 call ga_put(g_x,ii,ii,1,1,xx,1) 989 call mcscf_hessv( geom, basis, nbf, nclosed, nact, 990 $ tol2e, .false., pflg, lshift, dm1, dm2, 991 $ g_movecs, g_ifock, g_afock, g_gfock, 992 $ g_coul, g_exch, g_x, g_prod ) 993 call ga_get(g_prod,1,orlen,1,1,dbl_mb(k_hv),orlen) 994 call hv_writev(88,orlen,dbl_mb(k_hv)) 995 call dcopy(orlen,dbl_mb(k_hv),1,hh(1,ii),1) 996 enddo 997 enddo 998c 999c 1000c 1001c$$$ if (.not.ma_push_get(MT_DBL, orlen, 'Hy', l_hd, k_hd)) 1002c$$$ $ call errquit('mcscf: cannot allocate',0) 1003c$$$c 1004c$$$c 1005c$$$c 1006c$$$ open(unit=11,file='hessian.ascii',form='formatted', 1007c$$$ $ status='unknown') 1008c$$$ call hv_write(11,incr,orlen,dbl_mb(k_hv)) 1009c$$$ call hv_writev(12,orlen,dbl_mb(k_hv)) 1010c$$$ open(unit=12,file='hessian',form='unformatted', 1011c$$$ $ status='unknown') 1012c$$$ close(11) 1013c$$$ close(12) 1014 1015c$$$10000 continue 1016c$$$ xx = 1.d0 1017c$$$ do i=1,orlen 1018c$$$ call ga_zero(g_x) 1019c$$$ call ga_put(g_x,i,i,1,1,xx,1) 1020c$$$ call mcscf_hessv( geom, basis, nbf, nclosed, nact, 1021c$$$ $ tol2e, oskel, dm1, dm2, g_movecs, 1022c$$$ $ g_ifock, g_afock, g_gfock, 1023c$$$ $ g_coul, g_exch, g_x, g_prod ) 1024c$$$ call ga_get(g_prod,i,i,1,1,dbl_mb(k_hd+i-1),1) 1025c$$$ enddo 1026c$$$ write(6,900) 1027c$$$ 900 format('Exact Hessian diagonal') 1028c$$$ write(6,901) (dbl_mb(k_hd+i-1),i=1,orlen) 1029c$$$ 901 format(10f12.6) 1030 1031c$$$ if (.not.ma_pop_stack(l_hd)) 1032c$$$ $ call errquit('mcscf: damn',0) 1033c 1034c Clean up 1035c 1036 if (.not.ma_pop_stack(l_hy)) 1037 $ call errquit('mcscf: damn',0, MA_ERR) 1038 if (.not.ma_pop_stack(l_hv)) 1039 $ call errquit('mcscf: damn',0, MA_ERR) 1040 1041 return 1042 end 1043 1044 1045 1046 1047c 1048c Make J and K integrals for closed + active operators 1049c and dump to disk 1050c 1051 subroutine makeJK( basis, nbf, nclosed, nact, g_movecs ) 1052 implicit none 1053#include "errquit.fh" 1054#include "global.fh" 1055 integer basis, nbf, nclosed, nact, g_movecs 1056 integer mo_lo, mo_hi, noper, nn, i 1057 integer g_coul, g_exch 1058 double precision tmp(1000) 1059 logical oskel 1060 integer ga_create_JKblocked 1061 external ga_create_JKblocked 1062 data oskel/.false./ 1063 1064 mo_lo = 1 1065 mo_hi = (nclosed+nact) 1066 noper = (mo_hi*(mo_hi+1))/2 1067 nn = nbf*nbf 1068 g_coul = ga_create_JKblocked(noper,nbf,nbf,'Coulomb Oper') 1069 g_exch = ga_create_JKblocked(noper,nbf,nbf,'X Oper') 1070 call moints_build_6x( basis, oskel, 1071 $ mo_lo, mo_lo, mo_hi, 1, nbf, 1072 $ g_movecs, g_coul, .true., 1073 $ g_exch, .true., 16, .false., .false. ) 1074 open(unit=88,file='JKints',status='unknown',form='unformatted') 1075 do i=1,noper 1076 call ga_get(g_coul,1,nn,i,i,tmp,1) 1077 call hv_writev(88,nn,tmp) 1078 enddo 1079 do i=1,noper 1080 call ga_get(g_exch,1,nn,i,i,tmp,1) 1081 call hv_writev(88,nn,tmp) 1082 enddo 1083 close(88) 1084 if (.not. ga_destroy(g_coul)) call errquit('mcscf: ga?',0, 1085 & GA_ERR) 1086 if (.not. ga_destroy(g_exch)) call errquit('mcscf: ga?',0, 1087 & GA_ERR) 1088 return 1089 end 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100#endif 1101 1102