From 26be853c18189096dea671da1766724540f0a859 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Tue, 25 Jun 2019 16:46:14 +0200 Subject: [PATCH] Cleaning --- src/casscf/bielec.irp.f | 237 ++++--- src/casscf/bielec_create.irp.f | 118 ---- src/casscf/bielec_natorb.irp.f | 273 ++++++++ src/casscf/densities.irp.f | 377 +++++----- src/casscf/det_manip.irp.f | 249 ++++--- src/casscf/driver_wdens.irp.f | 104 +-- src/casscf/gradient.irp.f | 460 ++++++------ src/casscf/hessian.irp.f | 1204 ++++++++++++++++---------------- src/casscf/mcscf_fock.irp.f | 141 ++-- src/casscf/natorb.irp.f | 915 ++++++++++-------------- src/casscf/natorb_casscf.irp.f | 65 -- src/casscf/tot_en.irp.f | 203 +++--- 12 files changed, 2126 insertions(+), 2220 deletions(-) delete mode 100644 src/casscf/bielec_create.irp.f create mode 100644 src/casscf/bielec_natorb.irp.f delete mode 100644 src/casscf/natorb_casscf.irp.f diff --git a/src/casscf/bielec.irp.f b/src/casscf/bielec.irp.f index a1ec155d..9bb953f8 100644 --- a/src/casscf/bielec.irp.f +++ b/src/casscf/bielec.irp.f @@ -1,104 +1,151 @@ -! -*- F90 -*- - BEGIN_PROVIDER[real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] -&BEGIN_PROVIDER[real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] -BEGIN_DOC -! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active -! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active -! indices are unshifted orbital numbers -! all integrals are read from files -END_DOC - implicit none - integer :: i,j,p,q,indx,kk - real*8 :: hhh - logical :: lread + BEGIN_PROVIDER [real*8, bielec_PQxx, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] + BEGIN_DOC + ! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active + ! indices are unshifted orbital numbers + END_DOC + implicit none + integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 + double precision, allocatable :: integrals_array(:,:) + real*8 :: mo_two_e_integral + + allocate(integrals_array(mo_num,mo_num)) + + bielec_PQxx = 0.d0 + + do i=1,n_core_orb + ii=list_core(i) + do j=i,n_core_orb + jj=list_core(j) + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PQxx(p,q,i,j)=integrals_array(p,q) + bielec_PQxx(p,q,j,i)=integrals_array(p,q) + end do + end do + end do + do j=1,n_act_orb + jj=list_act(j) + j3=j+n_core_orb + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PQxx(p,q,i,j3)=integrals_array(p,q) + bielec_PQxx(p,q,j3,i)=integrals_array(p,q) + end do + end do + end do + end do - do i=1,n_core_orb+n_act_orb - do j=1,n_core_orb+n_act_orb - do p=1,mo_num - do q=1,mo_num - bielec_PQxx(p,q,i,j)=0.D0 - bielec_PxxQ(p,i,j,q)=0.D0 - end do - end do - end do - end do - open(unit=12,form='formatted',status='old',file='bielec_PQxx.tmp') - lread=.true. - indx=0 - do while (lread) - read(12,*,iostat=kk) p,q,i,j,hhh - if (kk.ne.0) then - lread=.false. - else -! stored with p.le.q, i.le.j - bielec_PQxx(p,q,i,j)=hhh - bielec_PQxx(q,p,i,j)=hhh - bielec_PQxx(q,p,j,i)=hhh - bielec_PQxx(p,q,j,i)=hhh - indx+=1 - end if - end do - close(12) - write(6,*) ' read ',indx,' integrals PQxx into core ' - - open(unit=12,form='formatted',status='old',file='bielec_PxxQ.tmp') - lread=.true. - indx=0 - do while (lread) - read(12,*,iostat=kk) p,i,j,q,hhh - if (kk.ne.0) then - lread=.false. - else -! stored with (ip).le.(jq) - bielec_PxxQ(p,i,j,q)=hhh - bielec_PxxQ(q,j,i,p)=hhh - indx+=1 - end if - end do - write(6,*) ' read ',indx,' integrals PxxQ into core ' - close(12) - write(6,*) ' provided integrals (PQ|xx) and (Px|xQ) ' + ! (ij|pq) + do i=1,n_act_orb + ii=list_act(i) + i3=i+n_core_orb + do j=i,n_act_orb + jj=list_act(j) + j3=j+n_core_orb + call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PQxx(p,q,i3,j3)=integrals_array(p,q) + bielec_PQxx(p,q,j3,i3)=integrals_array(p,q) + end do + end do + end do + end do + + write(6,*) ' provided integrals (PQ|xx) ' END_PROVIDER -BEGIN_PROVIDER[real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] -BEGIN_DOC -! bielecCI : integrals (tu|vp) with p arbitrary, tuv active -! index p runs over the whole basis, t,u,v only over the active orbitals -! inegrals read from file -END_DOC - implicit none - integer :: i,j,k,p,t,u,v,kk,indx - real*8 :: hhh - logical :: lread - write(6,*) ' reading integrals bielecCI ' - do i=1,n_act_orb - do j=1,n_act_orb - do k=1,n_act_orb - do p=1,mo_num - bielecCI(i,k,j,p)=0.D0 - end do - end do - end do - end do +BEGIN_PROVIDER [real*8, bielec_PxxQ, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] + BEGIN_DOC + ! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active + ! indices are unshifted orbital numbers + END_DOC + implicit none + integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 + double precision, allocatable :: integrals_array(:,:) + real*8 :: mo_two_e_integral + + allocate(integrals_array(mo_num,mo_num)) + + bielec_PxxQ = 0.d0 + + do i=1,n_core_orb + ii=list_core(i) + do j=i,n_core_orb + jj=list_core(j) + call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PxxQ(p,i,j,q)=integrals_array(p,q) + bielec_PxxQ(p,j,i,q)=integrals_array(q,p) + end do + end do + end do + do j=1,n_act_orb + jj=list_act(j) + j3=j+n_core_orb + call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PxxQ(p,i,j3,q)=integrals_array(p,q) + bielec_PxxQ(p,j3,i,q)=integrals_array(q,p) + end do + end do + end do + end do - open(unit=12,form='formatted',status='old',file='bielecCI.tmp') - lread=.true. - indx=0 - do while (lread) - read(12,*,iostat=kk) i,j,k,p,hhh - if (kk.ne.0) then - lread=.false. - else - bielecCI(i,j,k,p)=hhh - bielecCI(j,i,k,p)=hhh - indx+=1 - end if - end do - write(6,*) ' read ',indx,' integrals (aa|aP) into core ' - close(12) - write(6,*) ' provided integrals (tu|xP) ' + + ! (ip|qj) + do i=1,n_act_orb + ii=list_act(i) + i3=i+n_core_orb + do j=i,n_act_orb + jj=list_act(j) + j3=j+n_core_orb + call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + do q=1,mo_num + bielec_PxxQ(p,i3,j3,q)=integrals_array(p,q) + bielec_PxxQ(p,j3,i3,q)=integrals_array(q,p) + end do + end do + end do + end do + write(6,*) ' provided integrals (Px|xQ) ' +END_PROVIDER + + +BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)] + BEGIN_DOC + ! bielecCI : integrals (tu|vp) with p arbitrary, tuv active + ! index p runs over the whole basis, t,u,v only over the active orbitals + END_DOC + implicit none + integer :: i,j,k,p,t,u,v + double precision, allocatable :: integrals_array(:) + real*8 :: mo_two_e_integral + + allocate(integrals_array(mo_num)) + + do i=1,n_act_orb + t=list_act(i) + do j=1,n_act_orb + u=list_act(j) + do k=1,n_act_orb + v=list_act(k) + ! (tu|vp) + call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array,mo_integrals_map) + do p=1,mo_num + bielecCI(i,k,j,p)=integrals_array(p) + end do + end do + end do + end do + write(6,*) ' provided integrals (tu|xP) ' END_PROVIDER diff --git a/src/casscf/bielec_create.irp.f b/src/casscf/bielec_create.irp.f deleted file mode 100644 index 7e6d16c8..00000000 --- a/src/casscf/bielec_create.irp.f +++ /dev/null @@ -1,118 +0,0 @@ -! -*- F90 -*- - BEGIN_PROVIDER[real*8, bielec_PQxxtmp, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] -&BEGIN_PROVIDER[real*8, bielec_PxxQtmp, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] -&BEGIN_PROVIDER[integer, num_PQxx_written] -&BEGIN_PROVIDER[integer, num_PxxQ_written] -BEGIN_DOC -! bielec_PQxx : integral (pq|xx) with p,q arbitrary, x core or active -! bielec_PxxQ : integral (px|xq) with p,q arbitrary, x core or active -! indices are unshifted orbital numbers -END_DOC - implicit none - integer :: i,j,ii,jj,p,q,i3,j3,t3,v3 - double precision, allocatable :: integrals_array1(:,:) - double precision, allocatable :: integrals_array2(:,:) - real*8 :: mo_two_e_integral - - allocate(integrals_array1(mo_num,mo_num)) - allocate(integrals_array2(mo_num,mo_num)) - - do i=1,n_core_orb+n_act_orb - do j=1,n_core_orb+n_act_orb - do p=1,mo_num - do q=1,mo_num - bielec_PQxxtmp(p,q,i,j)=0.D0 - bielec_PxxQtmp(p,i,j,q)=0.D0 - end do - end do - end do - end do - - do i=1,n_core_orb - ii=list_core(i) - do j=i,n_core_orb - jj=list_core(j) -! (ij|pq) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map) -! (ip|qj) - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxxtmp(p,q,i,j)=integrals_array1(p,q) - bielec_PQxxtmp(p,q,j,i)=integrals_array1(p,q) - bielec_PxxQtmp(p,i,j,q)=integrals_array2(p,q) - bielec_PxxQtmp(p,j,i,q)=integrals_array2(q,p) - end do - end do - end do - do j=1,n_act_orb - jj=list_act(j) - j3=j+n_core_orb -! (ij|pq) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map) -! (ip|qj) - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxxtmp(p,q,i,j3)=integrals_array1(p,q) - bielec_PQxxtmp(p,q,j3,i)=integrals_array1(p,q) - bielec_PxxQtmp(p,i,j3,q)=integrals_array2(p,q) - bielec_PxxQtmp(p,j3,i,q)=integrals_array2(q,p) - end do - end do - end do - end do - do i=1,n_act_orb - ii=list_act(i) - i3=i+n_core_orb - do j=i,n_act_orb - jj=list_act(j) - j3=j+n_core_orb -! (ij|pq) - call get_mo_two_e_integrals_i1j1(ii,jj,mo_num,integrals_array1,mo_integrals_map) -! (ip|qj) - call get_mo_two_e_integrals_ij (ii,jj,mo_num,integrals_array2,mo_integrals_map) - do p=1,mo_num - do q=1,mo_num - bielec_PQxxtmp(p,q,i3,j3)=integrals_array1(p,q) - bielec_PQxxtmp(p,q,j3,i3)=integrals_array1(p,q) - bielec_PxxQtmp(p,i3,j3,q)=integrals_array2(p,q) - bielec_PxxQtmp(p,j3,i3,q)=integrals_array2(q,p) - end do - end do - end do - end do - write(6,*) ' provided integrals (PQ|xx) ' - write(6,*) ' provided integrals (Px|xQ) ' -!!$ write(6,*) ' 1 1 1 2 = ',bielec_PQxxtmp(2,2,2,3),bielec_PxxQtmp(2,2,2,3) -END_PROVIDER - -BEGIN_PROVIDER[real*8, bielecCItmp, (n_act_orb,n_act_orb,n_act_orb, mo_num)] -BEGIN_DOC -! bielecCI : integrals (tu|vp) with p arbitrary, tuv active -! index p runs over the whole basis, t,u,v only over the active orbitals -END_DOC - implicit none - integer :: i,j,k,p,t,u,v - double precision, allocatable :: integrals_array1(:) - real*8 :: mo_two_e_integral - - allocate(integrals_array1(mo_num)) - - do i=1,n_act_orb - t=list_act(i) - do j=1,n_act_orb - u=list_act(j) - do k=1,n_act_orb - v=list_act(k) -! (tu|vp) - call get_mo_two_e_integrals(t,u,v,mo_num,integrals_array1,mo_integrals_map) - do p=1,mo_num - bielecCItmp(i,k,j,p)=integrals_array1(p) - end do - end do - end do - end do - write(6,*) ' provided integrals (tu|xP) ' -END_PROVIDER - diff --git a/src/casscf/bielec_natorb.irp.f b/src/casscf/bielec_natorb.irp.f new file mode 100644 index 00000000..2f1e43eb --- /dev/null +++ b/src/casscf/bielec_natorb.irp.f @@ -0,0 +1,273 @@ + BEGIN_PROVIDER [real*8, bielec_PQxx_no, (mo_num, mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb)] + BEGIN_DOC + ! integral (pq|xx) in the basis of natural MOs + ! indices are unshifted orbital numbers + END_DOC + implicit none + integer :: i,j,k,l,t,u,p,q,pp + real*8 :: d(n_act_orb) + + bielec_PQxx_no(:,:,:,:) = bielec_PQxx(:,:,:,:) + + do j=1,mo_num + do k=1,n_core_orb+n_act_orb + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PQxx_no(list_act(q),j,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PQxx_no(list_act(p),j,k,l)=d(p) + end do + end do + end do + end do + ! 2nd quarter + do j=1,mo_num + do k=1,n_core_orb+n_act_orb + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PQxx_no(j,list_act(q),k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PQxx_no(j,list_act(p),k,l)=d(p) + end do + end do + end do + end do + ! 3rd quarter + do j=1,mo_num + do k=1,mo_num + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PQxx_no(j,k,n_core_orb+q,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PQxx_no(j,k,n_core_orb+p,l)=d(p) + end do + end do + end do + end do + ! 4th quarter + do j=1,mo_num + do k=1,mo_num + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PQxx_no(j,k,l,n_core_orb+q)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PQxx_no(j,k,l,n_core_orb+p)=d(p) + end do + end do + end do + end do + write(6,*) ' transformed PQxx' + +END_PROVIDER + + + +BEGIN_PROVIDER [real*8, bielec_PxxQ_no, (mo_num,n_core_orb+n_act_orb,n_core_orb+n_act_orb, mo_num)] + BEGIN_DOC + ! integral (px|xq) in the basis of natural MOs + ! indices are unshifted orbital numbers + END_DOC + implicit none + integer :: i,j,k,l,t,u,p,q,pp + real*8 :: d(n_act_orb) + + bielec_PxxQ_no(:,:,:,:) = bielec_PxxQ(:,:,:,:) + + do j=1,mo_num + do k=1,n_core_orb+n_act_orb + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PxxQ_no(list_act(q),k,l,j)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PxxQ_no(list_act(p),k,l,j)=d(p) + end do + end do + end do + end do + ! 2nd quarter + do j=1,mo_num + do k=1,n_core_orb+n_act_orb + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PxxQ_no(j,k,l,list_act(q))*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PxxQ_no(j,k,l,list_act(p))=d(p) + end do + end do + end do + end do + ! 3rd quarter + do j=1,mo_num + do k=1,mo_num + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PxxQ_no(j,n_core_orb+q,l,k)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PxxQ_no(j,n_core_orb+p,l,k)=d(p) + end do + end do + end do + end do + ! 4th quarter + do j=1,mo_num + do k=1,mo_num + do l=1,n_core_orb+n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielec_PxxQ_no(j,l,n_core_orb+q,k)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielec_PxxQ_no(j,l,n_core_orb+p,k)=d(p) + end do + end do + end do + end do + write(6,*) ' transformed PxxQ ' + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, bielecCI_no, (n_act_orb,n_act_orb,n_act_orb, mo_num)] + BEGIN_DOC + ! integrals (tu|vp) in the basis of natural MOs + ! index p runs over the whole basis, t,u,v only over the active orbitals + END_DOC + implicit none + integer :: i,j,k,l,t,u,p,q,pp + real*8 :: d(n_act_orb) + + bielecCI_no(:,:,:,:) = bielecCI(:,:,:,:) + + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielecCI_no(q,j,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielecCI_no(p,j,k,l)=d(p) + end do + end do + end do + end do + ! 2nd quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielecCI_no(j,q,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielecCI_no(j,p,k,l)=d(p) + end do + end do + end do + end do + ! 3rd quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielecCI_no(j,k,q,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielecCI_no(j,k,p,l)=d(p) + end do + end do + end do + end do + ! 4th quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=bielecCI_no(j,k,l,list_act(q))*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + bielecCI_no(j,k,l,list_act(p))=d(p) + end do + end do + end do + end do + write(6,*) ' transformed tuvP ' + +END_PROVIDER + diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f index 77f5593e..6e8065e2 100644 --- a/src/casscf/densities.irp.f +++ b/src/casscf/densities.irp.f @@ -1,177 +1,216 @@ -! -*- F90 -*- -use bitmasks ! you need to include the bitmasks_module.f90 features +use bitmasks - BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ] -&BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] -BEGIN_DOC -! the first-order density matrix in the basis of the starting MOs -! the second-order density matrix in the basis of the starting MOs -! matrices are state averaged -! -! we use the spin-free generators of mono-excitations -! E_pq destroys q and creates p -! D_pq = <0|E_pq|0> = D_qp -! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0> -! -END_DOC - implicit none - integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart - integer :: ierr - integer(bit_kind), allocatable :: det_mu(:,:) - integer(bit_kind), allocatable :: det_mu_ex(:,:) - integer(bit_kind), allocatable :: det_mu_ex1(:,:) - integer(bit_kind), allocatable :: det_mu_ex11(:,:) - integer(bit_kind), allocatable :: det_mu_ex12(:,:) - integer(bit_kind), allocatable :: det_mu_ex2(:,:) - integer(bit_kind), allocatable :: det_mu_ex21(:,:) - integer(bit_kind), allocatable :: det_mu_ex22(:,:) - real*8 :: phase1,phase11,phase12,phase2,phase21,phase22 - integer :: nu1,nu2,nu11,nu12,nu21,nu22 - integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22 - real*8 :: cI_mu(N_states),term - allocate(det_mu(N_int,2)) - allocate(det_mu_ex(N_int,2)) - allocate(det_mu_ex1(N_int,2)) - allocate(det_mu_ex11(N_int,2)) - allocate(det_mu_ex12(N_int,2)) - allocate(det_mu_ex2(N_int,2)) - allocate(det_mu_ex21(N_int,2)) - allocate(det_mu_ex22(N_int,2)) - - write(6,*) ' providing density matrices D0 and P0 ' - -! set all to zero - do t=1,n_act_orb - do u=1,n_act_orb - D0tu(u,t)=0.D0 - do v=1,n_act_orb - do x=1,n_act_orb - P0tuvx(x,v,u,t)=0.D0 - end do - end do - end do - end do - -! first loop: we apply E_tu, once for D_tu, once for -P_tvvu - do mu=1,n_det - call det_extract(det_mu,mu,N_int) - do istate=1,n_states - cI_mu(istate)=psi_coef(mu,istate) - end do - do t=1,n_act_orb - ipart=list_act(t) - do u=1,n_act_orb - ihole=list_act(u) -! apply E_tu - call det_copy(det_mu,det_mu_ex1,N_int) - call det_copy(det_mu,det_mu_ex2,N_int) - call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & +BEGIN_PROVIDER [real*8, D0tu, (n_act_orb,n_act_orb) ] + BEGIN_DOC + ! the first-order density matrix in the basis of the starting MOs + ! matrices are state averaged + ! + ! we use the spin-free generators of mono-excitations + ! E_pq destroys q and creates p + ! D_pq = <0|E_pq|0> = D_qp + ! + END_DOC + implicit none + integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart + integer :: ierr + integer(bit_kind) :: det_mu(N_int,2) + integer(bit_kind) :: det_mu_ex(N_int,2) + integer(bit_kind) :: det_mu_ex1(N_int,2) + integer(bit_kind) :: det_mu_ex2(N_int,2) + real*8 :: phase1,phase2,term + integer :: nu1,nu2 + integer :: ierr1,ierr2 + real*8 :: cI_mu(N_states) + + write(6,*) ' providing density matrices D0 and P0 ' + + D0tu = 0.d0 + + ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu + do mu=1,n_det + call det_extract(det_mu,mu,N_int) + do istate=1,n_states + cI_mu(istate)=psi_coef(mu,istate) + end do + do t=1,n_act_orb + ipart=list_act(t) + do u=1,n_act_orb + ihole=list_act(u) + ! apply E_tu + call det_copy(det_mu,det_mu_ex1,N_int) + call det_copy(det_mu,det_mu_ex2,N_int) + call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) -! det_mu_ex1 is in the list - if (nu1.ne.-1) then + ! det_mu_ex1 is in the list + if (nu1.ne.-1) then do istate=1,n_states - term=cI_mu(istate)*psi_coef(nu1,istate)*phase1 - D0tu(t,u)+=term -! and we fill P0_tvvu - do v=1,n_act_orb - P0tuvx(t,v,v,u)-=term - end do + term=cI_mu(istate)*psi_coef(nu1,istate)*phase1 + D0tu(t,u)+=term end do - end if -! det_mu_ex2 is in the list - if (nu2.ne.-1) then + end if + ! det_mu_ex2 is in the list + if (nu2.ne.-1) then do istate=1,n_states - term=cI_mu(istate)*psi_coef(nu2,istate)*phase2 - D0tu(t,u)+=term - do v=1,n_act_orb - P0tuvx(t,v,v,u)-=term - end do + term=cI_mu(istate)*psi_coef(nu2,istate)*phase2 + D0tu(t,u)+=term end do - end if - end do - end do + end if end do -! now we do the double excitation E_tu E_vx |0> - do mu=1,n_det - call det_extract(det_mu,mu,N_int) - do istate=1,n_states - cI_mu(istate)=psi_coef(mu,istate) - end do - do v=1,n_act_orb - ipart=list_act(v) - do x=1,n_act_orb - ihole=list_act(x) -! apply E_vx - call det_copy(det_mu,det_mu_ex1,N_int) - call det_copy(det_mu,det_mu_ex2,N_int) - call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & - ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) -! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0> - if (ierr1.eq.1) then - do t=1,n_act_orb - jpart=list_act(t) - do u=1,n_act_orb - jhole=list_act(u) - call det_copy(det_mu_ex1,det_mu_ex11,N_int) - call det_copy(det_mu_ex1,det_mu_ex12,N_int) - call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11 & - ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12) - if (nu11.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate) & - *phase11*phase1 - end do - end if - if (nu12.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate) & - *phase12*phase1 - end do - end if - end do - end do - end if - -! we apply E_tu to the second resultant determinant - if (ierr2.eq.1) then - do t=1,n_act_orb - jpart=list_act(t) - do u=1,n_act_orb - jhole=list_act(u) - call det_copy(det_mu_ex2,det_mu_ex21,N_int) - call det_copy(det_mu_ex2,det_mu_ex22,N_int) - call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21 & - ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22) - if (nu21.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate) & - *phase21*phase2 - end do - end if - if (nu22.ne.-1) then - do istate=1,n_states - P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate) & - *phase22*phase2 - end do - end if - end do - end do - end if - - end do - end do - end do - -! we average by just dividing by the number of states - do x=1,n_act_orb - do v=1,n_act_orb - D0tu(v,x)*=1.0D0/dble(N_states) - do u=1,n_act_orb - do t=1,n_act_orb - P0tuvx(t,u,v,x)*=0.5D0/dble(N_states) - end do - end do - end do - end do - + end do + end do + + ! we average by just dividing by the number of states + do x=1,n_act_orb + do v=1,n_act_orb + D0tu(v,x)*=1.0D0/dble(N_states) + end do + end do + +END_PROVIDER + +BEGIN_PROVIDER [real*8, P0tuvx, (n_act_orb,n_act_orb,n_act_orb,n_act_orb) ] + BEGIN_DOC + ! the second-order density matrix in the basis of the starting MOs + ! matrices are state averaged + ! + ! we use the spin-free generators of mono-excitations + ! E_pq destroys q and creates p + ! D_pq = <0|E_pq|0> = D_qp + ! P_pqrs = 1/2 <0|E_pq E_rs - delta_qr E_ps|0> + ! + END_DOC + implicit none + integer :: t,u,v,x,mu,nu,istate,ispin,jspin,ihole,ipart,jhole,jpart + integer :: ierr + real*8 :: phase1,phase11,phase12,phase2,phase21,phase22 + integer :: nu1,nu2,nu11,nu12,nu21,nu22 + integer :: ierr1,ierr2,ierr11,ierr12,ierr21,ierr22 + real*8 :: cI_mu(N_states),term + integer(bit_kind), dimension(N_int,2) :: det_mu, det_mu_ex + integer(bit_kind), dimension(N_int,2) :: det_mu_ex1, det_mu_ex11, det_mu_ex12 + integer(bit_kind), dimension(N_int,2) :: det_mu_ex2, det_mu_ex21, det_mu_ex22 + + write(6,*) ' providing density matrices D0 and P0 ' + + P0tuvx = 0.d0 + + ! first loop: we apply E_tu, once for D_tu, once for -P_tvvu + do mu=1,n_det + call det_extract(det_mu,mu,N_int) + do istate=1,n_states + cI_mu(istate)=psi_coef(mu,istate) + end do + do t=1,n_act_orb + ipart=list_act(t) + do u=1,n_act_orb + ihole=list_act(u) + ! apply E_tu + call det_copy(det_mu,det_mu_ex1,N_int) + call det_copy(det_mu,det_mu_ex2,N_int) + call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & + ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) + ! det_mu_ex1 is in the list + if (nu1.ne.-1) then + do istate=1,n_states + term=cI_mu(istate)*psi_coef(nu1,istate)*phase1 + ! and we fill P0_tvvu + do v=1,n_act_orb + P0tuvx(t,v,v,u)-=term + end do + end do + end if + ! det_mu_ex2 is in the list + if (nu2.ne.-1) then + do istate=1,n_states + term=cI_mu(istate)*psi_coef(nu2,istate)*phase2 + do v=1,n_act_orb + P0tuvx(t,v,v,u)-=term + end do + end do + end if + end do + end do + end do + ! now we do the double excitation E_tu E_vx |0> + do mu=1,n_det + call det_extract(det_mu,mu,N_int) + do istate=1,n_states + cI_mu(istate)=psi_coef(mu,istate) + end do + do v=1,n_act_orb + ipart=list_act(v) + do x=1,n_act_orb + ihole=list_act(x) + ! apply E_vx + call det_copy(det_mu,det_mu_ex1,N_int) + call det_copy(det_mu,det_mu_ex2,N_int) + call do_spinfree_mono_excitation(det_mu,det_mu_ex1 & + ,det_mu_ex2,nu1,nu2,ihole,ipart,phase1,phase2,ierr1,ierr2) + ! we apply E_tu to the first resultant determinant, thus E_tu E_vx |0> + if (ierr1.eq.1) then + do t=1,n_act_orb + jpart=list_act(t) + do u=1,n_act_orb + jhole=list_act(u) + call det_copy(det_mu_ex1,det_mu_ex11,N_int) + call det_copy(det_mu_ex1,det_mu_ex12,N_int) + call do_spinfree_mono_excitation(det_mu_ex1,det_mu_ex11& + ,det_mu_ex12,nu11,nu12,jhole,jpart,phase11,phase12,ierr11,ierr12) + if (nu11.ne.-1) then + do istate=1,n_states + P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu11,istate)& + *phase11*phase1 + end do + end if + if (nu12.ne.-1) then + do istate=1,n_states + P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu12,istate)& + *phase12*phase1 + end do + end if + end do + end do + end if + + ! we apply E_tu to the second resultant determinant + if (ierr2.eq.1) then + do t=1,n_act_orb + jpart=list_act(t) + do u=1,n_act_orb + jhole=list_act(u) + call det_copy(det_mu_ex2,det_mu_ex21,N_int) + call det_copy(det_mu_ex2,det_mu_ex22,N_int) + call do_spinfree_mono_excitation(det_mu_ex2,det_mu_ex21& + ,det_mu_ex22,nu21,nu22,jhole,jpart,phase21,phase22,ierr21,ierr22) + if (nu21.ne.-1) then + do istate=1,n_states + P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu21,istate)& + *phase21*phase2 + end do + end if + if (nu22.ne.-1) then + do istate=1,n_states + P0tuvx(t,u,v,x)+=cI_mu(istate)*psi_coef(nu22,istate)& + *phase22*phase2 + end do + end if + end do + end do + end if + + end do + end do + end do + + ! we average by just dividing by the number of states + do x=1,n_act_orb + do v=1,n_act_orb + do u=1,n_act_orb + do t=1,n_act_orb + P0tuvx(t,u,v,x)*=0.5D0/dble(N_states) + end do + end do + end do + end do + END_PROVIDER diff --git a/src/casscf/det_manip.irp.f b/src/casscf/det_manip.irp.f index c8e6c08a..adf90196 100644 --- a/src/casscf/det_manip.irp.f +++ b/src/casscf/det_manip.irp.f @@ -1,131 +1,130 @@ -! -*- F90 -*- -use bitmasks ! you need to include the bitmasks_module.f90 features +use bitmasks - subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, & - ispin,phase,ierr) -BEGIN_DOC -! we create the mono-excitation, and determine, if possible, -! the phase and the number in the list of determinants -END_DOC - implicit none - integer(bit_kind) :: key1(N_int,2),key2(N_int,2) - integer(bit_kind), allocatable :: keytmp(:,:) - integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin - real*8 :: phase - logical :: found - allocate(keytmp(N_int,2)) - - nu=-1 - phase=1.D0 - ierr=0 - call det_copy(key1,key2,N_int) -! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin -! call print_det(key2,N_int) - call do_single_excitation(key2,ihole,ipart,ispin,ierr) -! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin -! call print_det(key2,N_int) -! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr - if (ierr.eq.1) then -! excitation is possible -! get the phase - call get_single_excitation(key1,key2,exc,phase,N_int) -! get the number in the list - found=.false. - nu=0 - do while (.not.found) - nu+=1 - if (nu.gt.N_det) then -! the determinant is possible, but not in the list - found=.true. - nu=-1 - else - call det_extract(keytmp,nu,N_int) -integer :: i,ii - found=.true. - do ii=1,2 - do i=1,N_int - if (keytmp(i,ii).ne.key2(i,ii)) then - found=.false. - end if - end do - end do - end if +subroutine do_signed_mono_excitation(key1,key2,nu,ihole,ipart, & + ispin,phase,ierr) + BEGIN_DOC + ! we create the mono-excitation, and determine, if possible, + ! the phase and the number in the list of determinants + END_DOC + implicit none + integer(bit_kind) :: key1(N_int,2),key2(N_int,2) + integer(bit_kind), allocatable :: keytmp(:,:) + integer :: exc(0:2,2,2),ihole,ipart,ierr,nu,ispin + real*8 :: phase + logical :: found + allocate(keytmp(N_int,2)) + + nu=-1 + phase=1.D0 + ierr=0 + call det_copy(key1,key2,N_int) + ! write(6,*) ' key2 before excitation ',ihole,' -> ',ipart,' spin = ',ispin + ! call print_det(key2,N_int) + call do_single_excitation(key2,ihole,ipart,ispin,ierr) + ! write(6,*) ' key2 after ',ihole,' -> ',ipart,' spin = ',ispin + ! call print_det(key2,N_int) + ! write(6,*) ' excitation ',ihole,' -> ',ipart,' gives ierr = ',ierr + if (ierr.eq.1) then + ! excitation is possible + ! get the phase + call get_single_excitation(key1,key2,exc,phase,N_int) + ! get the number in the list + found=.false. + nu=0 + do while (.not.found) + nu+=1 + if (nu.gt.N_det) then + ! the determinant is possible, but not in the list + found=.true. + nu=-1 + else + call det_extract(keytmp,nu,N_int) + integer :: i,ii + found=.true. + do ii=1,2 + do i=1,N_int + if (keytmp(i,ii).ne.key2(i,ii)) then + found=.false. + end if end do -! if (found) then -! if (nu.eq.-1) then -! write(6,*) ' image not found in the list, thus nu = ',nu -! else -! write(6,*) ' found in the list as No ',nu,' phase = ',phase -! end if -! end if - end if -! -! we found the new string, the phase, and possibly the number in the list -! - end subroutine do_signed_mono_excitation + end do + end if + end do + ! if (found) then + ! if (nu.eq.-1) then + ! write(6,*) ' image not found in the list, thus nu = ',nu + ! else + ! write(6,*) ' found in the list as No ',nu,' phase = ',phase + ! end if + ! end if + end if + ! + ! we found the new string, the phase, and possibly the number in the list + ! +end subroutine do_signed_mono_excitation - subroutine det_extract(key,nu,Nint) -BEGIN_DOC -! extract a determinant from the list of determinants -END_DOC - implicit none - integer :: ispin,i,nu,Nint - integer(bit_kind) :: key(Nint,2) - do ispin=1,2 - do i=1,Nint - key(i,ispin)=psi_det(i,ispin,nu) - end do - end do - end subroutine det_extract +subroutine det_extract(key,nu,Nint) + BEGIN_DOC + ! extract a determinant from the list of determinants + END_DOC + implicit none + integer :: ispin,i,nu,Nint + integer(bit_kind) :: key(Nint,2) + do ispin=1,2 + do i=1,Nint + key(i,ispin)=psi_det(i,ispin,nu) + end do + end do +end subroutine det_extract - subroutine det_copy(key1,key2,Nint) - use bitmasks ! you need to include the bitmasks_module.f90 features -BEGIN_DOC -! copy a determinant from key1 to key2 -END_DOC - implicit none - integer :: ispin,i,Nint - integer(bit_kind) :: key1(Nint,2),key2(Nint,2) - do ispin=1,2 - do i=1,Nint - key2(i,ispin)=key1(i,ispin) - end do - end do - end subroutine det_copy - - subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 & - ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr) -BEGIN_DOC -! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q) -! we may create two determinants as result -! -END_DOC - implicit none - integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2) - integer(bit_kind) :: key_out2(N_int,2) - integer :: ihole,ipart,ierr,jerr,nu1,nu2 - integer :: ispin - real*8 :: phase1,phase2 - -! write(6,*) ' applying E_',ipart,ihole,' on determinant ' -! call print_det(key_in,N_int) - -! spin alpha - ispin=1 - call do_signed_mono_excitation(key_in,key_out1,nu1,ihole & - ,ipart,ispin,phase1,ierr) -! if (ierr.eq.1) then -! write(6,*) ' 1 result is ',nu1,phase1 -! call print_det(key_out1,N_int) -! end if -! spin beta - ispin=2 - call do_signed_mono_excitation(key_in,key_out2,nu2,ihole & - ,ipart,ispin,phase2,jerr) -! if (jerr.eq.1) then -! write(6,*) ' 2 result is ',nu2,phase2 -! call print_det(key_out2,N_int) -! end if +subroutine det_copy(key1,key2,Nint) + use bitmasks ! you need to include the bitmasks_module.f90 features + BEGIN_DOC + ! copy a determinant from key1 to key2 + END_DOC + implicit none + integer :: ispin,i,Nint + integer(bit_kind) :: key1(Nint,2),key2(Nint,2) + do ispin=1,2 + do i=1,Nint + key2(i,ispin)=key1(i,ispin) + end do + end do +end subroutine det_copy - end subroutine do_spinfree_mono_excitation +subroutine do_spinfree_mono_excitation(key_in,key_out1,key_out2 & + ,nu1,nu2,ihole,ipart,phase1,phase2,ierr,jerr) + BEGIN_DOC + ! we create the spin-free mono-excitation E_pq=(a^+_p a_q + a^+_P a_Q) + ! we may create two determinants as result + ! + END_DOC + implicit none + integer(bit_kind) :: key_in(N_int,2),key_out1(N_int,2) + integer(bit_kind) :: key_out2(N_int,2) + integer :: ihole,ipart,ierr,jerr,nu1,nu2 + integer :: ispin + real*8 :: phase1,phase2 + + ! write(6,*) ' applying E_',ipart,ihole,' on determinant ' + ! call print_det(key_in,N_int) + + ! spin alpha + ispin=1 + call do_signed_mono_excitation(key_in,key_out1,nu1,ihole & + ,ipart,ispin,phase1,ierr) + ! if (ierr.eq.1) then + ! write(6,*) ' 1 result is ',nu1,phase1 + ! call print_det(key_out1,N_int) + ! end if + ! spin beta + ispin=2 + call do_signed_mono_excitation(key_in,key_out2,nu2,ihole & + ,ipart,ispin,phase2,jerr) + ! if (jerr.eq.1) then + ! write(6,*) ' 2 result is ',nu2,phase2 + ! call print_det(key_out2,N_int) + ! end if + +end subroutine do_spinfree_mono_excitation diff --git a/src/casscf/driver_wdens.irp.f b/src/casscf/driver_wdens.irp.f index 263e8441..5a3863a3 100644 --- a/src/casscf/driver_wdens.irp.f +++ b/src/casscf/driver_wdens.irp.f @@ -7,58 +7,13 @@ write(6,*) ' generating natural orbitals ' write(6,*) write(6,*) - call trf_to_natorb write(6,*) ' all data available ! ' write(6,*) ' writing out files ' - open(unit=12,file='D0tu.dat',form='formatted',status='unknown') - do p=1,n_act_orb - do q=1,n_act_orb - if (abs(D0tu(p,q)).gt.1.D-12) then - write(12,'(2i8,E20.12)') p,q,D0tu(p,q) - end if - end do - end do - close(12) - + call trf_to_natorb real*8 :: approx,np,nq,nr,ns logical :: lpq,lrs,lps,lqr - open(unit=12,file='P0tuvx.dat',form='formatted',status='unknown') - do p=1,n_act_orb - np=D0tu(p,p) - do q=1,n_act_orb - lpq=p.eq.q - nq=D0tu(q,q) - do r=1,n_act_orb - lqr=q.eq.r - nr=D0tu(r,r) - do s=1,n_act_orb - lrs=r.eq.s - lps=p.eq.s - approx=0.D0 - if (lpq.and.lrs) then - if (lqr) then -! pppp - approx=0.5D0*np*(np-1.D0) - else -! pprr - approx=0.5D0*np*nr - end if - else - if (lps.and.lqr.and..not.lpq) then -! pqqp - approx=-0.25D0*np*nq - end if - end if - if (abs(P0tuvx(p,q,r,s)).gt.1.D-12) then - write(12,'(4i4,2E20.12)') p,q,r,s,P0tuvx(p,q,r,s),approx - end if - end do - end do - end do - end do - close(12) open(unit=12,form='formatted',status='unknown',file='onetrf.tmp') indx=0 @@ -74,63 +29,6 @@ logical :: lpq,lrs,lps,lqr close(12) - open(unit=12,form='formatted',status='unknown',file='bielec_PQxx.tmp') - indx=0 - do p=1,mo_num - do q=p,mo_num - do r=1,n_core_orb+n_act_orb - do s=r,n_core_orb+n_act_orb - if (abs(bielec_PQxxtmp(p,q,r,s)).gt.1.D-12) then - write(12,'(4i8,E20.12)') p,q,r,s,bielec_PQxxtmp(p,q,r,s) - indx+=1 - end if - end do - end do - end do - end do - write(6,*) ' wrote ',indx,' integrals (PQ|xx)' - close(12) - - open(unit=12,form='formatted',status='unknown',file='bielec_PxxQ.tmp') - indx=0 - do p=1,mo_num - do q=1,n_core_orb+n_act_orb - do r=q,n_core_orb+n_act_orb -integer ::s_start - if (q.eq.r) then - s_start=p - else - s_start=1 - end if - do s=s_start,mo_num - if (abs(bielec_PxxQtmp(p,q,r,s)).gt.1.D-12) then - write(12,'(4i8,E20.12)') p,q,r,s,bielec_PxxQtmp(p,q,r,s) - indx+=1 - end if - end do - end do - end do - end do - write(6,*) ' wrote ',indx,' integrals (Px|xQ)' - close(12) - - open(unit=12,form='formatted',status='unknown',file='bielecCI.tmp') - indx=0 - do p=1,n_act_orb - do q=p,n_act_orb - do r=1,n_act_orb - do s=1,mo_num - if (abs(bielecCItmp(p,q,r,s)).gt.1.D-12) then - write(12,'(4i8,E20.12)') p,q,r,s,bielecCItmp(p,q,r,s) - indx+=1 - end if - end do - end do - end do - end do - write(6,*) ' wrote ',indx,' integrals (tu|xP)' - close(12) - write(6,*) write(6,*) ' creating new orbitals ' do i=1,mo_num diff --git a/src/casscf/gradient.irp.f b/src/casscf/gradient.irp.f index d35d96ed..606bf12b 100644 --- a/src/casscf/gradient.irp.f +++ b/src/casscf/gradient.irp.f @@ -1,251 +1,249 @@ -! -*- F90 -*- - -use bitmasks ! you need to include the bitmasks_module.f90 features +use bitmasks BEGIN_PROVIDER [ integer, nMonoEx ] -BEGIN_DOC -! -END_DOC - implicit none - nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb - write(6,*) ' nMonoEx = ',nMonoEx + BEGIN_DOC + ! Number of single excitations + END_DOC + implicit none + nMonoEx=n_core_orb*n_act_orb+n_core_orb*n_virt_orb+n_act_orb*n_virt_orb + write(6,*) ' nMonoEx = ',nMonoEx END_PROVIDER BEGIN_PROVIDER [integer, excit, (2,nMonoEx)] &BEGIN_PROVIDER [character*3, excit_class, (nMonoEx)] -BEGIN_DOC -! a list of the orbitals involved in the excitation -END_DOC - - implicit none - integer :: i,t,a,ii,tt,aa,indx - indx=0 - do ii=1,n_core_orb - i=list_core(ii) - do tt=1,n_act_orb - t=list_act(tt) - indx+=1 - excit(1,indx)=i - excit(2,indx)=t - excit_class(indx)='c-a' - end do - end do - - do ii=1,n_core_orb - i=list_core(ii) - do aa=1,n_virt_orb - a=list_virt(aa) - indx+=1 - excit(1,indx)=i - excit(2,indx)=a - excit_class(indx)='c-v' - end do - end do - - do tt=1,n_act_orb + BEGIN_DOC + ! a list of the orbitals involved in the excitation + END_DOC + + implicit none + integer :: i,t,a,ii,tt,aa,indx + indx=0 + do ii=1,n_core_orb + i=list_core(ii) + do tt=1,n_act_orb t=list_act(tt) - do aa=1,n_virt_orb - a=list_virt(aa) - indx+=1 - excit(1,indx)=t - excit(2,indx)=a - excit_class(indx)='a-v' - end do - end do - - if (bavard) then - write(6,*) ' Filled the table of the Monoexcitations ' - do indx=1,nMonoEx - write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & - ,excit(2,indx),' ',excit_class(indx) - end do - end if - + indx+=1 + excit(1,indx)=i + excit(2,indx)=t + excit_class(indx)='c-a' + end do + end do + + do ii=1,n_core_orb + i=list_core(ii) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=i + excit(2,indx)=a + excit_class(indx)='c-v' + end do + end do + + do tt=1,n_act_orb + t=list_act(tt) + do aa=1,n_virt_orb + a=list_virt(aa) + indx+=1 + excit(1,indx)=t + excit(2,indx)=a + excit_class(indx)='a-v' + end do + end do + + if (bavard) then + write(6,*) ' Filled the table of the Monoexcitations ' + do indx=1,nMonoEx + write(6,*) ' ex ',indx,' : ',excit(1,indx),' -> ' & + ,excit(2,indx),' ',excit_class(indx) + end do + end if + END_PROVIDER BEGIN_PROVIDER [real*8, gradvec, (nMonoEx)] -BEGIN_DOC -! calculate the orbital gradient by hand, i.e. for -! each determinant I we determine the string E_pq |I> (alpha and beta -! separately) and generate -! sum_I c_I is then the pq component of the orbital -! gradient -! E_pq = a^+_pa_q + a^+_Pa_Q -END_DOC - implicit none - integer :: ii,tt,aa,indx,ihole,ipart,istate - real*8 :: res - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - call calc_grad_elem(ihole,ipart,res) - gradvec(indx)=res - end do - -real*8 :: norm_grad - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec(indx)*gradvec(indx) - end do - norm_grad=sqrt(norm_grad) - write(6,*) - write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad - write(6,*) - - + BEGIN_DOC + ! calculate the orbital gradient by hand, i.e. for + ! each determinant I we determine the string E_pq |I> (alpha and beta + ! separately) and generate + ! sum_I c_I is then the pq component of the orbital + ! gradient + ! E_pq = a^+_pa_q + a^+_Pa_Q + END_DOC + implicit none + integer :: ii,tt,aa,indx,ihole,ipart,istate + real*8 :: res + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + call calc_grad_elem(ihole,ipart,res) + gradvec(indx)=res + end do + + real*8 :: norm_grad + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec(indx)*gradvec(indx) + end do + norm_grad=sqrt(norm_grad) + write(6,*) + write(6,*) ' Norm of the orbital gradient (via <0|EH|0>) : ', norm_grad + write(6,*) + + END_PROVIDER - subroutine calc_grad_elem(ihole,ipart,res) -BEGIN_DOC -! eq 18 of Siegbahn et al, Physica Scripta 1980 -! we calculate 2 , q=hole, p=particle -END_DOC - implicit none - integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) - real*8 :: i_H_psi_array(N_states),phase - allocate(det_mu(N_int,2)) - allocate(det_mu_ex(N_int,2)) - - res=0.D0 - - do mu=1,n_det -! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 -! do the monoexcitation on it - call det_copy(det_mu,det_mu_ex,N_int) - call do_signed_mono_excitation(det_mu,det_mu_ex,nu & - ,ihole,ipart,ispin,phase,ierr) - if (ierr.eq.1) then -! write(6,*) -! write(6,*) ' mu = ',mu -! call print_det(det_mu,N_int) -! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin -! call print_det(det_mu_ex,N_int) - call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase - end do -! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res - end if +subroutine calc_grad_elem(ihole,ipart,res) + BEGIN_DOC + ! eq 18 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 , q=hole, p=particle + END_DOC + implicit none + integer :: ihole,ipart,mu,iii,ispin,ierr,nu,istate + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:),det_mu_ex(:,:) + real*8 :: i_H_psi_array(N_states),phase + allocate(det_mu(N_int,2)) + allocate(det_mu_ex(N_int,2)) + + res=0.D0 + + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation on it + call det_copy(det_mu,det_mu_ex,N_int) + call do_signed_mono_excitation(det_mu,det_mu_ex,nu & + ,ihole,ipart,ispin,phase,ierr) + if (ierr.eq.1) then + ! write(6,*) + ! write(6,*) ' mu = ',mu + ! call print_det(det_mu,N_int) + ! write(6,*) ' generated nu = ',nu,' for excitation ',ihole,' -> ',ipart,' ierr = ',ierr,' phase = ',phase,' ispin = ',ispin + ! call print_det(det_mu_ex,N_int) + call i_H_psi(det_mu_ex,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase end do - end do - -! state-averaged gradient - res*=2.D0/dble(N_states) - - end subroutine calc_grad_elem + ! write(6,*) ' contribution = ',i_H_psi_array(1)*psi_coef(mu,1)*phase,res + end if + end do + end do + + ! state-averaged gradient + res*=2.D0/dble(N_states) + +end subroutine calc_grad_elem BEGIN_PROVIDER [real*8, gradvec2, (nMonoEx)] -BEGIN_DOC -! calculate the orbital gradient from density -! matrices and integrals; Siegbahn et al, Phys Scr 1980 -! eqs 14 a,b,c -END_DOC - implicit none - integer :: i,t,a,indx - real*8 :: gradvec_it,gradvec_ia,gradvec_ta - real*8 :: norm_grad - - indx=0 - do i=1,n_core_orb - do t=1,n_act_orb - indx+=1 - gradvec2(indx)=gradvec_it(i,t) - end do - end do - - do i=1,n_core_orb - do a=1,n_virt_orb - indx+=1 - gradvec2(indx)=gradvec_ia(i,a) - end do - end do - - do t=1,n_act_orb - do a=1,n_virt_orb - indx+=1 - gradvec2(indx)=gradvec_ta(t,a) - end do - end do - - norm_grad=0.d0 - do indx=1,nMonoEx - norm_grad+=gradvec2(indx)*gradvec2(indx) - end do - norm_grad=sqrt(norm_grad) - write(6,*) - write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad - write(6,*) - + BEGIN_DOC + ! calculate the orbital gradient from density + ! matrices and integrals; Siegbahn et al, Phys Scr 1980 + ! eqs 14 a,b,c + END_DOC + implicit none + integer :: i,t,a,indx + real*8 :: gradvec_it,gradvec_ia,gradvec_ta + real*8 :: norm_grad + + indx=0 + do i=1,n_core_orb + do t=1,n_act_orb + indx+=1 + gradvec2(indx)=gradvec_it(i,t) + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + indx+=1 + gradvec2(indx)=gradvec_ia(i,a) + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + indx+=1 + gradvec2(indx)=gradvec_ta(t,a) + end do + end do + + norm_grad=0.d0 + do indx=1,nMonoEx + norm_grad+=gradvec2(indx)*gradvec2(indx) + end do + norm_grad=sqrt(norm_grad) + write(6,*) + write(6,*) ' Norm of the orbital gradient (via D, P and integrals): ', norm_grad + write(6,*) + END_PROVIDER - real*8 function gradvec_it(i,t) -BEGIN_DOC -! the orbital gradient core -> active -! we assume natural orbitals -END_DOC - implicit none - integer :: i,t +real*8 function gradvec_it(i,t) + BEGIN_DOC + ! the orbital gradient core -> active + ! we assume natural orbitals + END_DOC + implicit none + integer :: i,t + + integer :: ii,tt,v,vv,x,y + integer :: x3,y3 + + ii=list_core(i) + tt=list_act(t) + gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) + gradvec_it-=occnum(tt)*Fipq(ii,tt) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + x3=x+n_core_orb + do y=1,n_act_orb + y3=y+n_core_orb + gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx_no(ii,vv,x3,y3) + end do + end do + end do + gradvec_it*=2.D0 +end function gradvec_it - integer :: ii,tt,v,vv,x,y - integer :: x3,y3 +real*8 function gradvec_ia(i,a) + BEGIN_DOC + ! the orbital gradient core -> virtual + END_DOC + implicit none + integer :: i,a,ii,aa + + ii=list_core(i) + aa=list_virt(a) + gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) + gradvec_ia*=2.D0 + +end function gradvec_ia - ii=list_core(i) - tt=list_act(t) - gradvec_it=2.D0*(Fipq(tt,ii)+Fapq(tt,ii)) - gradvec_it-=occnum(tt)*Fipq(ii,tt) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb - x3=x+n_core_orb - do y=1,n_act_orb - y3=y+n_core_orb - gradvec_it-=2.D0*P0tuvx_no(t,v,x,y)*bielec_PQxx(ii,vv,x3,y3) - end do - end do - end do - gradvec_it*=2.D0 - end function gradvec_it - - real*8 function gradvec_ia(i,a) -BEGIN_DOC -! the orbital gradient core -> virtual -END_DOC - implicit none - integer :: i,a,ii,aa - - ii=list_core(i) - aa=list_virt(a) - gradvec_ia=2.D0*(Fipq(aa,ii)+Fapq(aa,ii)) - gradvec_ia*=2.D0 - - end function gradvec_ia - - real*8 function gradvec_ta(t,a) -BEGIN_DOC -! the orbital gradient active -> virtual -! we assume natural orbitals -END_DOC - implicit none - integer :: t,a,tt,aa,v,vv,x,y - - tt=list_act(t) - aa=list_virt(a) - gradvec_ta=0.D0 - gradvec_ta+=occnum(tt)*Fipq(aa,tt) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa) - end do - end do - end do - gradvec_ta*=2.D0 - - end function gradvec_ta +real*8 function gradvec_ta(t,a) + BEGIN_DOC + ! the orbital gradient active -> virtual + ! we assume natural orbitals + END_DOC + implicit none + integer :: t,a,tt,aa,v,vv,x,y + + tt=list_act(t) + aa=list_virt(a) + gradvec_ta=0.D0 + gradvec_ta+=occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + gradvec_ta+=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) + end do + end do + end do + gradvec_ta*=2.D0 + +end function gradvec_ta diff --git a/src/casscf/hessian.irp.f b/src/casscf/hessian.irp.f index 4603d11e..65734a25 100644 --- a/src/casscf/hessian.irp.f +++ b/src/casscf/hessian.irp.f @@ -1,639 +1,637 @@ -! -*- F90 -*- - -use bitmasks ! you need to include the bitmasks_module.f90 features +use bitmasks BEGIN_PROVIDER [real*8, hessmat, (nMonoEx,nMonoEx)] -BEGIN_DOC -! calculate the orbital hessian 2 -! + + by hand, -! determinant per determinant, as for the gradient -! -! we assume that we have natural active orbitals -END_DOC - implicit none - integer :: indx,ihole,ipart - integer :: jndx,jhole,jpart - character*3 :: iexc,jexc - real*8 :: res - - write(6,*) ' providing Hessian matrix hessmat ' - write(6,*) ' nMonoEx = ',nMonoEx - - do indx=1,nMonoEx - do jndx=1,nMonoEx - hessmat(indx,jndx)=0.D0 - end do - end do - - do indx=1,nMonoEx - ihole=excit(1,indx) - ipart=excit(2,indx) - iexc=excit_class(indx) - do jndx=indx,nMonoEx - jhole=excit(1,jndx) - jpart=excit(2,jndx) - jexc=excit_class(jndx) - call calc_hess_elem(ihole,ipart,jhole,jpart,res) -! write(6,*) ' Hessian ',ihole,'->',ipart & -! ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res - hessmat(indx,jndx)=res - hessmat(jndx,indx)=res - end do - end do - + BEGIN_DOC + ! calculate the orbital hessian 2 + ! + + by hand, + ! determinant per determinant, as for the gradient + ! + ! we assume that we have natural active orbitals + END_DOC + implicit none + integer :: indx,ihole,ipart + integer :: jndx,jhole,jpart + character*3 :: iexc,jexc + real*8 :: res + + write(6,*) ' providing Hessian matrix hessmat ' + write(6,*) ' nMonoEx = ',nMonoEx + + do indx=1,nMonoEx + do jndx=1,nMonoEx + hessmat(indx,jndx)=0.D0 + end do + end do + + do indx=1,nMonoEx + ihole=excit(1,indx) + ipart=excit(2,indx) + iexc=excit_class(indx) + do jndx=indx,nMonoEx + jhole=excit(1,jndx) + jpart=excit(2,jndx) + jexc=excit_class(jndx) + call calc_hess_elem(ihole,ipart,jhole,jpart,res) + ! write(6,*) ' Hessian ',ihole,'->',ipart & + ! ,' (',iexc,')',jhole,'->',jpart,' (',jexc,')',res + hessmat(indx,jndx)=res + hessmat(jndx,indx)=res + end do + end do + END_PROVIDER - subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) -BEGIN_DOC -! eq 19 of Siegbahn et al, Physica Scripta 1980 -! we calculate 2 -! + + -! average over all states is performed. -! no transition between states. -END_DOC - implicit none - integer :: ihole,ipart,ispin,mu,istate - integer :: jhole,jpart,jspin - integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu - real*8 :: res - integer(bit_kind), allocatable :: det_mu(:,:) - integer(bit_kind), allocatable :: det_nu(:,:) - integer(bit_kind), allocatable :: det_mu_pq(:,:) - integer(bit_kind), allocatable :: det_mu_rs(:,:) - integer(bit_kind), allocatable :: det_nu_rs(:,:) - integer(bit_kind), allocatable :: det_mu_pqrs(:,:) - integer(bit_kind), allocatable :: det_mu_rspq(:,:) - real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 - real*8 :: i_H_j_element - allocate(det_mu(N_int,2)) - allocate(det_nu(N_int,2)) - allocate(det_mu_pq(N_int,2)) - allocate(det_mu_rs(N_int,2)) - allocate(det_nu_rs(N_int,2)) - allocate(det_mu_pqrs(N_int,2)) - allocate(det_mu_rspq(N_int,2)) - integer :: mu_pq_possible - integer :: mu_rs_possible - integer :: nu_rs_possible - integer :: mu_pqrs_possible - integer :: mu_rspq_possible - - res=0.D0 - -! the terms <0|E E H |0> - do mu=1,n_det -! get the string of the determinant - call det_extract(det_mu,mu,N_int) - do ispin=1,2 -! do the monoexcitation pq on it - call det_copy(det_mu,det_mu_pq,N_int) - call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & - ,ihole,ipart,ispin,phase,mu_pq_possible) - if (mu_pq_possible.eq.1) then -! possible, but not necessarily in the list -! do the second excitation - do jspin=1,2 - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs & +subroutine calc_hess_elem(ihole,ipart,jhole,jpart,res) + BEGIN_DOC + ! eq 19 of Siegbahn et al, Physica Scripta 1980 + ! we calculate 2 + ! + + + ! average over all states is performed. + ! no transition between states. + END_DOC + implicit none + integer :: ihole,ipart,ispin,mu,istate + integer :: jhole,jpart,jspin + integer :: mu_pq, mu_pqrs, mu_rs, mu_rspq, nu_rs,nu + real*8 :: res + integer(bit_kind), allocatable :: det_mu(:,:) + integer(bit_kind), allocatable :: det_nu(:,:) + integer(bit_kind), allocatable :: det_mu_pq(:,:) + integer(bit_kind), allocatable :: det_mu_rs(:,:) + integer(bit_kind), allocatable :: det_nu_rs(:,:) + integer(bit_kind), allocatable :: det_mu_pqrs(:,:) + integer(bit_kind), allocatable :: det_mu_rspq(:,:) + real*8 :: i_H_psi_array(N_states),phase,phase2,phase3 + real*8 :: i_H_j_element + allocate(det_mu(N_int,2)) + allocate(det_nu(N_int,2)) + allocate(det_mu_pq(N_int,2)) + allocate(det_mu_rs(N_int,2)) + allocate(det_nu_rs(N_int,2)) + allocate(det_mu_pqrs(N_int,2)) + allocate(det_mu_rspq(N_int,2)) + integer :: mu_pq_possible + integer :: mu_rs_possible + integer :: nu_rs_possible + integer :: mu_pqrs_possible + integer :: mu_rspq_possible + + res=0.D0 + + ! the terms <0|E E H |0> + do mu=1,n_det + ! get the string of the determinant + call det_extract(det_mu,mu,N_int) + do ispin=1,2 + ! do the monoexcitation pq on it + call det_copy(det_mu,det_mu_pq,N_int) + call do_signed_mono_excitation(det_mu,det_mu_pq,mu_pq & + ,ihole,ipart,ispin,phase,mu_pq_possible) + if (mu_pq_possible.eq.1) then + ! possible, but not necessarily in the list + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& ,jhole,jpart,jspin,phase2,mu_pqrs_possible) -! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 end do - end if -! try the de-excitation with opposite sign - call det_copy(det_mu_pq,det_mu_pqrs,N_int) - call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs & + end if + ! try the de-excitation with opposite sign + call det_copy(det_mu_pq,det_mu_pqrs,N_int) + call do_signed_mono_excitation(det_mu_pq,det_mu_pqrs,mu_pqrs& ,jpart,jhole,jspin,phase2,mu_pqrs_possible) - phase2=-phase2 -! excitation possible - if (mu_pqrs_possible.eq.1) then - call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) + phase2=-phase2 + ! excitation possible + if (mu_pqrs_possible.eq.1) then + call i_H_psi(det_mu_pqrs,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase*phase2 end do - end if - end do - end if -! exchange the notion of pq and rs -! do the monoexcitation rs on the initial determinant - call det_copy(det_mu,det_mu_rs,N_int) - call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & - ,jhole,jpart,ispin,phase2,mu_rs_possible) - if (mu_rs_possible.eq.1) then -! do the second excitation - do jspin=1,2 - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq & - ,ihole,ipart,jspin,phase3,mu_rspq_possible) -! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if -! we may try the de-excitation, with opposite sign - call det_copy(det_mu_rs,det_mu_rspq,N_int) - call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq & - ,ipart,ihole,jspin,phase3,mu_rspq_possible) - phase3=-phase3 -! excitation possible (of course, the result is outside the CAS) - if (mu_rspq_possible.eq.1) then - call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & - ,N_det,N_det,N_states,i_H_psi_array) - do istate=1,N_states - res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 - end do - end if - end do - end if -! -! the operator E H E, we have to do a double loop over the determinants -! we still have the determinant mu_pq and the phase in memory - if (mu_pq_possible.eq.1) then - do nu=1,N_det - call det_extract(det_nu,nu,N_int) - do jspin=1,2 - call det_copy(det_nu,det_nu_rs,N_int) - call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & - ,jhole,jpart,jspin,phase2,nu_rs_possible) -! excitation possible ? - if (nu_rs_possible.eq.1) then - call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) - do istate=1,N_states - res+=2.D0*i_H_j_element*psi_coef(mu,istate) & - *psi_coef(nu,istate)*phase*phase2 - end do - end if - end do - end do - end if + end if end do - end do - -! state-averaged Hessian - res*=1.D0/dble(N_states) - - end subroutine calc_hess_elem + end if + ! exchange the notion of pq and rs + ! do the monoexcitation rs on the initial determinant + call det_copy(det_mu,det_mu_rs,N_int) + call do_signed_mono_excitation(det_mu,det_mu_rs,mu_rs & + ,jhole,jpart,ispin,phase2,mu_rs_possible) + if (mu_rs_possible.eq.1) then + ! do the second excitation + do jspin=1,2 + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ihole,ipart,jspin,phase3,mu_rspq_possible) + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + ! we may try the de-excitation, with opposite sign + call det_copy(det_mu_rs,det_mu_rspq,N_int) + call do_signed_mono_excitation(det_mu_rs,det_mu_rspq,mu_rspq& + ,ipart,ihole,jspin,phase3,mu_rspq_possible) + phase3=-phase3 + ! excitation possible (of course, the result is outside the CAS) + if (mu_rspq_possible.eq.1) then + call i_H_psi(det_mu_rspq,psi_det,psi_coef,N_int & + ,N_det,N_det,N_states,i_H_psi_array) + do istate=1,N_states + res+=i_H_psi_array(istate)*psi_coef(mu,istate)*phase2*phase3 + end do + end if + end do + end if + ! + ! the operator E H E, we have to do a double loop over the determinants + ! we still have the determinant mu_pq and the phase in memory + if (mu_pq_possible.eq.1) then + do nu=1,N_det + call det_extract(det_nu,nu,N_int) + do jspin=1,2 + call det_copy(det_nu,det_nu_rs,N_int) + call do_signed_mono_excitation(det_nu,det_nu_rs,nu_rs & + ,jhole,jpart,jspin,phase2,nu_rs_possible) + ! excitation possible ? + if (nu_rs_possible.eq.1) then + call i_H_j(det_mu_pq,det_nu_rs,N_int,i_H_j_element) + do istate=1,N_states + res+=2.D0*i_H_j_element*psi_coef(mu,istate) & + *psi_coef(nu,istate)*phase*phase2 + end do + end if + end do + end do + end if + end do + end do + + ! state-averaged Hessian + res*=1.D0/dble(N_states) + +end subroutine calc_hess_elem BEGIN_PROVIDER [real*8, hessmat2, (nMonoEx,nMonoEx)] -BEGIN_DOC -! explicit hessian matrix from density matrices and integrals -! of course, this will be used for a direct Davidson procedure later -! we will not store the matrix in real life -! formulas are broken down as functions for the 6 classes of matrix elements -! -END_DOC - implicit none - integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart - - real*8 :: hessmat_itju - real*8 :: hessmat_itja - real*8 :: hessmat_itua - real*8 :: hessmat_iajb - real*8 :: hessmat_iatb - real*8 :: hessmat_taub - - write(6,*) ' providing Hessian matrix hessmat2 ' - write(6,*) ' nMonoEx = ',nMonoEx - - indx=1 - do i=1,n_core_orb - do t=1,n_act_orb - jndx=indx - do j=i,n_core_orb - if (i.eq.j) then - ustart=t - else - ustart=1 - end if - do u=ustart,n_act_orb - hessmat2(indx,jndx)=hessmat_itju(i,t,j,u) - hessmat2(jndx,indx)=hessmat2(indx,jndx) -! write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx) - jndx+=1 - end do - end do - do j=1,n_core_orb - do a=1,n_virt_orb - hessmat2(indx,jndx)=hessmat_itja(i,t,j,a) - hessmat2(jndx,indx)=hessmat2(indx,jndx) - jndx+=1 - end do - end do - do u=1,n_act_orb - do a=1,n_virt_orb - hessmat2(indx,jndx)=hessmat_itua(i,t,u,a) - hessmat2(jndx,indx)=hessmat2(indx,jndx) - jndx+=1 - end do - end do - indx+=1 + BEGIN_DOC + ! explicit hessian matrix from density matrices and integrals + ! of course, this will be used for a direct Davidson procedure later + ! we will not store the matrix in real life + ! formulas are broken down as functions for the 6 classes of matrix elements + ! + END_DOC + implicit none + integer :: i,j,t,u,a,b,indx,jndx,bstart,ustart + + real*8 :: hessmat_itju + real*8 :: hessmat_itja + real*8 :: hessmat_itua + real*8 :: hessmat_iajb + real*8 :: hessmat_iatb + real*8 :: hessmat_taub + + write(6,*) ' providing Hessian matrix hessmat2 ' + write(6,*) ' nMonoEx = ',nMonoEx + + indx=1 + do i=1,n_core_orb + do t=1,n_act_orb + jndx=indx + do j=i,n_core_orb + if (i.eq.j) then + ustart=t + else + ustart=1 + end if + do u=ustart,n_act_orb + hessmat2(indx,jndx)=hessmat_itju(i,t,j,u) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + ! write(6,*) ' result I :',i,t,j,u,indx,jndx,hessmat(indx,jndx),hessmat2(indx,jndx) + jndx+=1 end do - end do - - do i=1,n_core_orb + end do + do j=1,n_core_orb do a=1,n_virt_orb - jndx=indx - do j=i,n_core_orb - if (i.eq.j) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(indx,jndx)=hessmat_iajb(i,a,j,b) - hessmat2(jndx,indx)=hessmat2(indx,jndx) - jndx+=1 - end do - end do - do t=1,n_act_orb - do b=1,n_virt_orb - hessmat2(indx,jndx)=hessmat_iatb(i,a,t,b) - hessmat2(jndx,indx)=hessmat2(indx,jndx) - jndx+=1 - end do - end do - indx+=1 + hessmat2(indx,jndx)=hessmat_itja(i,t,j,a) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 end do - end do - - do t=1,n_act_orb + end do + do u=1,n_act_orb do a=1,n_virt_orb - jndx=indx - do u=t,n_act_orb - if (t.eq.u) then - bstart=a - else - bstart=1 - end if - do b=bstart,n_virt_orb - hessmat2(indx,jndx)=hessmat_taub(t,a,u,b) - hessmat2(jndx,indx)=hessmat2(indx,jndx) - jndx+=1 - end do - end do - indx+=1 + hessmat2(indx,jndx)=hessmat_itua(i,t,u,a) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 end do - end do - + end do + indx+=1 + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + jndx=indx + do j=i,n_core_orb + if (i.eq.j) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat2(indx,jndx)=hessmat_iajb(i,a,j,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + do t=1,n_act_orb + do b=1,n_virt_orb + hessmat2(indx,jndx)=hessmat_iatb(i,a,t,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + indx+=1 + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + jndx=indx + do u=t,n_act_orb + if (t.eq.u) then + bstart=a + else + bstart=1 + end if + do b=bstart,n_virt_orb + hessmat2(indx,jndx)=hessmat_taub(t,a,u,b) + hessmat2(jndx,indx)=hessmat2(indx,jndx) + jndx+=1 + end do + end do + indx+=1 + end do + end do + END_PROVIDER - real*8 function hessmat_itju(i,t,j,u) -BEGIN_DOC -! the orbital hessian for core->act,core->act -! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu -! -! we assume natural orbitals -END_DOC - implicit none - integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj - real*8 :: term,t2 - -! write(6,*) ' hessmat_itju ',i,t,j,u - ii=list_core(i) - tt=list_act(t) - if (i.eq.j) then - if (t.eq.u) then -! diagonal element - term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & - -2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) - term+=2.D0*(3.D0*bielec_pxxq(tt,i,i,tt)-bielec_pqxx(tt,tt,i,i)) - term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq(tt,i,i,tt) & - -bielec_pqxx(tt,tt,i,i)) - term-=occnum(tt)*Fipq(tt,tt) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb - xx=list_act(x) - term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(vv,xx,i,i) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq(vv,i,i,xx)) - do y=1,n_act_orb - term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx) - end do - end do - end do - else -! it/iu, t != u - uu=list_act(u) - term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu)) - term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & - -bielec_PQxx(tt,uu,i,j)) - term-=occnum(tt)*Fipq(uu,tt) - term-=(occnum(tt)+occnum(uu)) & - *(3.D0*bielec_PxxQ(tt,i,i,uu)-bielec_PQxx(uu,tt,i,i)) - do v=1,n_act_orb - vv=list_act(v) -! term-=D0tu(u,v)*Fipq(tt,vv) ! published, but inverting t and u seems more correct - do x=1,n_act_orb - xx=list_act(x) - term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,i) & - +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & - *bielec_pxxq(vv,i,i,xx)) - do y=1,n_act_orb - term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(u,v,y,xx) - end do - end do - end do -!!! write(6,*) ' direct diff ',i,t,j,u,term,term2 -!!! term=term2 - end if - else -! it/ju - jj=list_core(j) - uu=list_act(u) - if (t.eq.u) then - term=occnum(tt)*Fipq(ii,jj) - term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) - else - term=0.D0 - end if - term+=2.D0*(4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & - -bielec_PQxx(tt,uu,i,j)) - term-=(occnum(tt)+occnum(uu))* & - (4.D0*bielec_PxxQ(tt,i,j,uu)-bielec_PxxQ(uu,i,j,tt) & - -bielec_PQxx(uu,tt,i,j)) - do v=1,n_act_orb - vv=list_act(v) - do x=1,n_act_orb +real*8 function hessmat_itju(i,t,j,u) + BEGIN_DOC + ! the orbital hessian for core->act,core->act + ! i, t, j, u are list indices, the corresponding orbitals are ii,tt,jj,uu + ! + ! we assume natural orbitals + END_DOC + implicit none + integer :: i,t,j,u,ii,tt,uu,v,vv,x,xx,y,jj + real*8 :: term,t2 + + ! write(6,*) ' hessmat_itju ',i,t,j,u + ii=list_core(i) + tt=list_act(t) + if (i.eq.j) then + if (t.eq.u) then + ! diagonal element + term=occnum(tt)*Fipq(ii,ii)+2.D0*(Fipq(tt,tt)+Fapq(tt,tt)) & + -2.D0*(Fipq(ii,ii)+Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq_no(tt,i,i,tt)-bielec_pqxx_no(tt,tt,i,i)) + term-=2.D0*occnum(tt)*(3.D0*bielec_pxxq_no(tt,i,i,tt) & + -bielec_pqxx_no(tt,tt,i,i)) + term-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb xx=list_act(x) - term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx(vv,xx,i,j) & - +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & - *bielec_pxxq(vv,i,j,xx)) - end do - end do - end if - - term*=2.D0 - hessmat_itju=term - - end function hessmat_itju - - real*8 function hessmat_itja(i,t,j,a) -BEGIN_DOC -! the orbital hessian for core->act,core->virt -END_DOC - implicit none - integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y - real*8 :: term - -! write(6,*) ' hessmat_itja ',i,t,j,a -! it/ja - ii=list_core(i) - tt=list_act(t) - jj=list_core(j) - aa=list_virt(a) - term=2.D0*(4.D0*bielec_pxxq(aa,j,i,tt) & - -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt)) - term-=occnum(tt)*(4.D0*bielec_pxxq(aa,j,i,tt) & - -bielec_pqxx(aa,tt,i,j) -bielec_pxxq(aa,i,j,tt)) - if (i.eq.j) then - term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt)) - term-=0.5D0*occnum(tt)*Fipq(aa,tt) - do v=1,n_act_orb - do x=1,n_act_orb + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq_no(vv,i,i,xx)) do y=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,aa) + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) end do - end do end do - end if - term*=2.D0 - hessmat_itja=term + end do + else + ! it/iu, t != u + uu=list_act(u) + term=2.D0*(Fipq(tt,uu)+Fapq(tt,uu)) + term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(tt,uu,i,j)) + term-=occnum(tt)*Fipq(uu,tt) + term-=(occnum(tt)+occnum(uu)) & + *(3.D0*bielec_PxxQ_no(tt,i,i,uu)-bielec_PQxx_no(uu,tt,i,i)) + do v=1,n_act_orb + vv=list_act(v) + ! term-=D0tu(u,v)*Fipq(tt,vv) ! published, but inverting t and u seems more correct + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,i) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq_no(vv,i,i,xx)) + do y=1,n_act_orb + term-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(u,v,y,xx) + end do + end do + end do + !!! write(6,*) ' direct diff ',i,t,j,u,term,term2 + !!! term=term2 + end if + else + ! it/ju + jj=list_core(j) + uu=list_act(u) + if (t.eq.u) then + term=occnum(tt)*Fipq(ii,jj) + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + else + term=0.D0 + end if + term+=2.D0*(4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(tt,uu,i,j)) + term-=(occnum(tt)+occnum(uu))* & + (4.D0*bielec_PxxQ_no(tt,i,j,uu)-bielec_PxxQ_no(uu,i,j,tt) & + -bielec_PQxx_no(uu,tt,i,j)) + do v=1,n_act_orb + vv=list_act(v) + do x=1,n_act_orb + xx=list_act(x) + term+=2.D0*(P0tuvx_no(u,t,v,x)*bielec_pqxx_no(vv,xx,i,j) & + +(P0tuvx_no(u,x,v,t)+P0tuvx_no(u,x,t,v)) & + *bielec_pxxq_no(vv,i,j,xx)) + end do + end do + end if + + term*=2.D0 + hessmat_itju=term + +end function hessmat_itju - end function hessmat_itja +real*8 function hessmat_itja(i,t,j,a) + BEGIN_DOC + ! the orbital hessian for core->act,core->virt + END_DOC + implicit none + integer :: i,t,j,a,ii,tt,jj,aa,v,vv,x,y + real*8 :: term + + ! write(6,*) ' hessmat_itja ',i,t,j,a + ! it/ja + ii=list_core(i) + tt=list_act(t) + jj=list_core(j) + aa=list_virt(a) + term=2.D0*(4.D0*bielec_pxxq_no(aa,j,i,tt) & + -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) + term-=occnum(tt)*(4.D0*bielec_pxxq_no(aa,j,i,tt) & + -bielec_pqxx_no(aa,tt,i,j) -bielec_pxxq_no(aa,i,j,tt)) + if (i.eq.j) then + term+=2.D0*(Fipq(aa,tt)+Fapq(aa,tt)) + term-=0.5D0*occnum(tt)*Fipq(aa,tt) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,aa) + end do + end do + end do + end if + term*=2.D0 + hessmat_itja=term + +end function hessmat_itja - real*8 function hessmat_itua(i,t,u,a) -BEGIN_DOC -! the orbital hessian for core->act,act->virt -END_DOC - implicit none - integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 - real*8 :: term +real*8 function hessmat_itua(i,t,u,a) + BEGIN_DOC + ! the orbital hessian for core->act,act->virt + END_DOC + implicit none + integer :: i,t,u,a,ii,tt,uu,aa,v,vv,x,xx,u3,t3,v3 + real*8 :: term + + ! write(6,*) ' hessmat_itua ',i,t,u,a + ii=list_core(i) + tt=list_act(t) + t3=t+n_core_orb + uu=list_act(u) + u3=u+n_core_orb + aa=list_virt(a) + if (t.eq.u) then + term=-occnum(tt)*Fipq(aa,ii) + else + term=0.D0 + end if + term-=occnum(uu)*(bielec_pqxx_no(aa,ii,t3,u3)-4.D0*bielec_pqxx_no(aa,uu,t3,i)& + +bielec_pxxq_no(aa,t3,u3,ii)) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + integer :: x3 + xx=list_act(x) + x3=x+n_core_orb + term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,ii,v3,x3) & + +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u)) & + *bielec_pqxx_no(aa,xx,v3,i)) + end do + end do + if (t.eq.u) then + term+=Fipq(aa,ii)+Fapq(aa,ii) + end if + term*=2.D0 + hessmat_itua=term + +end function hessmat_itua -! write(6,*) ' hessmat_itua ',i,t,u,a - ii=list_core(i) - tt=list_act(t) - t3=t+n_core_orb - uu=list_act(u) - u3=u+n_core_orb - aa=list_virt(a) - if (t.eq.u) then - term=-occnum(tt)*Fipq(aa,ii) - else - term=0.D0 - end if - term-=occnum(uu)*(bielec_pqxx(aa,ii,t3,u3)-4.D0*bielec_pqxx(aa,uu,t3,i) & - +bielec_pxxq(aa,t3,u3,ii)) - do v=1,n_act_orb +real*8 function hessmat_iajb(i,a,j,b) + BEGIN_DOC + ! the orbital hessian for core->virt,core->virt + END_DOC + implicit none + integer :: i,a,j,b,ii,aa,jj,bb + real*8 :: term + ! write(6,*) ' hessmat_iajb ',i,a,j,b + + ii=list_core(i) + aa=list_virt(a) + if (i.eq.j) then + if (a.eq.b) then + ! ia/ia + term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii)) + term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,aa)-bielec_pqxx_no(aa,aa,i,i)) + else + bb=list_virt(b) + ! ia/ib + term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb)) + term+=2.D0*(3.D0*bielec_pxxq_no(aa,i,i,bb)-bielec_pqxx_no(aa,bb,i,i)) + end if + else + ! ia/jb + jj=list_core(j) + bb=list_virt(b) + term=2.D0*(4.D0*bielec_pxxq_no(aa,i,j,bb)-bielec_pqxx_no(aa,bb,i,j) & + -bielec_pxxq_no(aa,j,i,bb)) + if (a.eq.b) then + term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) + end if + end if + term*=2.D0 + hessmat_iajb=term + +end function hessmat_iajb + +real*8 function hessmat_iatb(i,a,t,b) + BEGIN_DOC + ! the orbital hessian for core->virt,act->virt + END_DOC + implicit none + integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 + real*8 :: term + + ! write(6,*) ' hessmat_iatb ',i,a,t,b + ii=list_core(i) + aa=list_virt(a) + tt=list_act(t) + bb=list_virt(b) + t3=t+n_core_orb + term=occnum(tt)*(4.D0*bielec_pxxq_no(aa,i,t3,bb)-bielec_pxxq_no(aa,t3,i,bb)& + -bielec_pqxx_no(aa,bb,i,t3)) + if (a.eq.b) then + term-=Fipq(tt,ii)+Fapq(tt,ii) + term-=0.5D0*occnum(tt)*Fipq(tt,ii) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,ii) + end do + end do + end do + end if + term*=2.D0 + hessmat_iatb=term + +end function hessmat_iatb + +real*8 function hessmat_taub(t,a,u,b) + BEGIN_DOC + ! the orbital hessian for act->virt,act->virt + END_DOC + implicit none + integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y + integer :: v3,x3 + real*8 :: term,t1,t2,t3 + + tt=list_act(t) + aa=list_virt(a) + if (t.eq.u) then + if (a.eq.b) then + ! ta/ta + t1=occnum(tt)*Fipq(aa,aa) + t2=0.D0 + t3=0.D0 + t1-=occnum(tt)*Fipq(tt,tt) + do v=1,n_act_orb vv=list_act(v) v3=v+n_core_orb do x=1,n_act_orb -integer :: x3 - xx=list_act(x) - x3=x+n_core_orb - term-=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,ii,v3,x3) & - +(P0tuvx_no(t,v,u,x)+P0tuvx_no(t,v,x,u)) & - *bielec_pqxx(aa,xx,v3,i)) - end do - end do - if (t.eq.u) then - term+=Fipq(aa,ii)+Fapq(aa,ii) - end if - term*=2.D0 - hessmat_itua=term - - end function hessmat_itua - - real*8 function hessmat_iajb(i,a,j,b) -BEGIN_DOC -! the orbital hessian for core->virt,core->virt -END_DOC - implicit none - integer :: i,a,j,b,ii,aa,jj,bb - real*8 :: term -! write(6,*) ' hessmat_iajb ',i,a,j,b - - ii=list_core(i) - aa=list_virt(a) - if (i.eq.j) then - if (a.eq.b) then -! ia/ia - term=2.D0*(Fipq(aa,aa)+Fapq(aa,aa)-Fipq(ii,ii)-Fapq(ii,ii)) - term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,aa)-bielec_pqxx(aa,aa,i,i)) - else - bb=list_virt(b) -! ia/ib - term=2.D0*(Fipq(aa,bb)+Fapq(aa,bb)) - term+=2.D0*(3.D0*bielec_pxxq(aa,i,i,bb)-bielec_pqxx(aa,bb,i,i)) - end if - else -! ia/jb - jj=list_core(j) - bb=list_virt(b) - term=2.D0*(4.D0*bielec_pxxq(aa,i,j,bb)-bielec_pqxx(aa,bb,i,j) & - -bielec_pxxq(aa,j,i,bb)) - if (a.eq.b) then - term-=2.D0*(Fipq(ii,jj)+Fapq(ii,jj)) - end if - end if - term*=2.D0 - hessmat_iajb=term - - end function hessmat_iajb - - real*8 function hessmat_iatb(i,a,t,b) -BEGIN_DOC -! the orbital hessian for core->virt,act->virt -END_DOC - implicit none - integer :: i,a,t,b,ii,aa,tt,bb,v,vv,x,y,v3,t3 - real*8 :: term - -! write(6,*) ' hessmat_iatb ',i,a,t,b - ii=list_core(i) - aa=list_virt(a) - tt=list_act(t) - bb=list_virt(b) - t3=t+n_core_orb - term=occnum(tt)*(4.D0*bielec_pxxq(aa,i,t3,bb)-bielec_pxxq(aa,t3,i,bb) & - -bielec_pqxx(aa,bb,i,t3)) - if (a.eq.b) then - term-=Fipq(tt,ii)+Fapq(tt,ii) - term-=0.5D0*occnum(tt)*Fipq(tt,ii) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,ii) - end do - end do - end do - end if - term*=2.D0 - hessmat_iatb=term - - end function hessmat_iatb - - real*8 function hessmat_taub(t,a,u,b) -BEGIN_DOC -! the orbital hessian for act->virt,act->virt -END_DOC - implicit none - integer :: t,a,u,b,tt,aa,uu,bb,v,vv,x,xx,y - integer :: v3,x3 - real*8 :: term,t1,t2,t3 - - tt=list_act(t) - aa=list_virt(a) - if (t.eq.u) then - if (a.eq.b) then -! ta/ta - t1=occnum(tt)*Fipq(aa,aa) - t2=0.D0 - t3=0.D0 - t1-=occnum(tt)*Fipq(tt,tt) - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_orb - t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,aa,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & - bielec_pxxq(aa,x3,v3,aa)) - do y=1,n_act_orb - t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI(t,v,y,xx) - end do - end do - end do - term=t1+t2+t3 -! write(6,*) ' Hess taub ',t,a,t1,t2,t3 - else - bb=list_virt(b) -! ta/tb b/=a - term=occnum(tt)*Fipq(aa,bb) - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_orb - term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & - *bielec_pxxq(aa,x3,v3,bb)) - end do - end do - end if - else -! ta/ub t/=u - uu=list_act(u) - bb=list_virt(b) - term=0.D0 - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_orb - do x=1,n_act_orb xx=list_act(x) x3=x+n_core_orb - term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx(aa,bb,v3,x3) & - +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & - *bielec_pxxq(aa,x3,v3,bb)) - end do - end do - if (a.eq.b) then - term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) - do v=1,n_act_orb - do x=1,n_act_orb - do y=1,n_act_orb - term-=P0tuvx_no(t,v,x,y)*bielecCI(x,y,v,uu) - term-=P0tuvx_no(u,v,x,y)*bielecCI(x,y,v,tt) - end do + t2+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,aa,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v))* & + bielec_pxxq_no(aa,x3,v3,aa)) + do y=1,n_act_orb + t3-=2.D0*P0tuvx_no(t,v,x,y)*bielecCI_no(t,v,y,xx) end do - end do - end if - - end if - - term*=2.D0 - hessmat_taub=term - - end function hessmat_taub + end do + end do + term=t1+t2+t3 + ! write(6,*) ' Hess taub ',t,a,t1,t2,t3 + else + bb=list_virt(b) + ! ta/tb b/=a + term=occnum(tt)*Fipq(aa,bb) + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + term+=2.D0*(P0tuvx_no(t,t,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,t)+P0tuvx_no(t,x,t,v)) & + *bielec_pxxq_no(aa,x3,v3,bb)) + end do + end do + end if + else + ! ta/ub t/=u + uu=list_act(u) + bb=list_virt(b) + term=0.D0 + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + term+=2.D0*(P0tuvx_no(t,u,v,x)*bielec_pqxx_no(aa,bb,v3,x3) & + +(P0tuvx_no(t,x,v,u)+P0tuvx_no(t,x,u,v)) & + *bielec_pxxq_no(aa,x3,v3,bb)) + end do + end do + if (a.eq.b) then + term-=0.5D0*(occnum(tt)*Fipq(uu,tt)+occnum(uu)*Fipq(tt,uu)) + do v=1,n_act_orb + do x=1,n_act_orb + do y=1,n_act_orb + term-=P0tuvx_no(t,v,x,y)*bielecCI_no(x,y,v,uu) + term-=P0tuvx_no(u,v,x,y)*bielecCI_no(x,y,v,tt) + end do + end do + end do + end if + + end if + + term*=2.D0 + hessmat_taub=term + +end function hessmat_taub BEGIN_PROVIDER [real*8, hessdiag, (nMonoEx)] -BEGIN_DOC -! the diagonal of the Hessian, needed for the Davidson procedure -END_DOC - implicit none - integer :: i,t,a,indx - real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub - - indx=0 - do i=1,n_core_orb - do t=1,n_act_orb - indx+=1 - hessdiag(indx)=hessmat_itju(i,t,i,t) - end do - end do - - do i=1,n_core_orb - do a=1,n_virt_orb - indx+=1 - hessdiag(indx)=hessmat_iajb(i,a,i,a) - end do - end do - - do t=1,n_act_orb - do a=1,n_virt_orb - indx+=1 - hessdiag(indx)=hessmat_taub(t,a,t,a) - end do - end do - + BEGIN_DOC + ! the diagonal of the Hessian, needed for the Davidson procedure + END_DOC + implicit none + integer :: i,t,a,indx + real*8 :: hessmat_itju,hessmat_iajb,hessmat_taub + + indx=0 + do i=1,n_core_orb + do t=1,n_act_orb + indx+=1 + hessdiag(indx)=hessmat_itju(i,t,i,t) + end do + end do + + do i=1,n_core_orb + do a=1,n_virt_orb + indx+=1 + hessdiag(indx)=hessmat_iajb(i,a,i,a) + end do + end do + + do t=1,n_act_orb + do a=1,n_virt_orb + indx+=1 + hessdiag(indx)=hessmat_taub(t,a,t,a) + end do + end do + END_PROVIDER diff --git a/src/casscf/mcscf_fock.irp.f b/src/casscf/mcscf_fock.irp.f index 301b1418..68845eb4 100644 --- a/src/casscf/mcscf_fock.irp.f +++ b/src/casscf/mcscf_fock.irp.f @@ -1,67 +1,80 @@ -! -*- F90 -*- - BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] -&BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] -BEGIN_DOC -! the inactive and the active Fock matrices, in molecular -! orbitals -! we create them in MOs, quite expensive -! -! for an implementation in AOs we need first the natural orbitals -! for forming an active density matrix in AOs -! -END_DOC - implicit none - double precision, allocatable :: integrals_array1(:,:) - double precision, allocatable :: integrals_array2(:,:) - integer :: p,q,k,kk,t,tt,u,uu - allocate(integrals_array1(mo_num,mo_num)) - allocate(integrals_array2(mo_num,mo_num)) - +BEGIN_PROVIDER [real*8, Fipq, (mo_num,mo_num) ] + BEGIN_DOC + ! the inactive Fock matrix, in molecular orbitals + END_DOC + implicit none + integer :: p,q,k,kk,t,tt,u,uu + + do q=1,mo_num + do p=1,mo_num + Fipq(p,q)=one_ints(p,q) + end do + end do + + ! the inactive Fock matrix + do k=1,n_core_orb + kk=list_core(k) + do q=1,mo_num do p=1,mo_num - do q=1,mo_num - Fipq(p,q)=one_ints(p,q) - Fapq(p,q)=0.D0 - end do + Fipq(p,q)+=2.D0*bielec_pqxx_no(p,q,k,k) -bielec_pxxq_no(p,k,k,q) end do - -! the inactive Fock matrix - do k=1,n_core_orb - kk=list_core(k) - do p=1,mo_num - do q=1,mo_num - Fipq(p,q)+=2.D0*bielec_pqxx(p,q,k,k) -bielec_pxxq(p,k,k,q) - end do - end do - end do - -! the active Fock matrix, D0tu is diagonal - do t=1,n_act_orb - tt=list_act(t) - do p=1,mo_num - do q=1,mo_num - Fapq(p,q)+=occnum(tt) & - *(bielec_pqxx(p,q,tt,tt)-0.5D0*bielec_pxxq(p,tt,tt,q)) - end do - end do - end do - -if (bavard) then -integer :: i - write(6,*) - write(6,*) ' the effective Fock matrix over MOs' - write(6,*) - - write(6,*) - write(6,*) ' the diagonal of the inactive effective Fock matrix ' - write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) - write(6,*) - write(6,*) - write(6,*) ' the diagonal of the active Fock matrix ' - write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) - write(6,*) -end if - - + end do + end do + + if (bavard) then + integer :: i + write(6,*) + write(6,*) ' the diagonal of the inactive effective Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) + write(6,*) + end if + + END_PROVIDER - - + + +BEGIN_PROVIDER [real*8, Fapq, (mo_num,mo_num) ] + BEGIN_DOC + ! the active active Fock matrix, in molecular orbitals + ! we create them in MOs, quite expensive + ! + ! for an implementation in AOs we need first the natural orbitals + ! for forming an active density matrix in AOs + ! + END_DOC + implicit none + integer :: p,q,k,kk,t,tt,u,uu + + Fapq = 0.d0 + + ! the active Fock matrix, D0tu is diagonal + do t=1,n_act_orb + tt=list_act(t) + do q=1,mo_num + do p=1,mo_num + Fapq(p,q)+=occnum(tt) & + *(bielec_pqxx_no(p,q,tt,tt)-0.5D0*bielec_pxxq_no(p,tt,tt,q)) + end do + end do + end do + + if (bavard) then + integer :: i + write(6,*) + write(6,*) ' the effective Fock matrix over MOs' + write(6,*) + + write(6,*) + write(6,*) ' the diagonal of the inactive effective Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fipq(i,i),i=1,mo_num) + write(6,*) + write(6,*) + write(6,*) ' the diagonal of the active Fock matrix ' + write(6,'(5(i3,F12.5))') (i,Fapq(i,i),i=1,mo_num) + write(6,*) + end if + + +END_PROVIDER + + diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f index a903260c..d2cc6736 100644 --- a/src/casscf/natorb.irp.f +++ b/src/casscf/natorb.irp.f @@ -1,548 +1,373 @@ -! -*- F90 -*- -! diagonalize D0tu -! save the diagonal somewhere, in inverse order -! 4-index-transform the 2-particle density matrix over active orbitals -! correct the bielectronic integrals -! correct the monoelectronic integrals -! put integrals on file, as well orbitals, and the density matrices -! - subroutine trf_to_natorb - implicit none - integer :: i,j,k,l,t,u,p,q,pp - real*8 :: eigval(n_act_orb),natorbsCI(n_act_orb,n_act_orb) - real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb) + BEGIN_PROVIDER [real*8, occnum, (mo_num)] + implicit none + BEGIN_DOC + ! MO occupation numbers + END_DOC + + integer :: i + occnum=0.D0 + do i=1,n_core_orb + occnum(list_core(i))=2.D0 + end do + + do i=1,n_act_orb + occnum(list_act(i))=occ_act(n_act_orb-i+1) + end do - call lapack_diag(eigval,natorbsCI,D0tu,n_act_orb,n_act_orb) - write(6,*) ' found occupation numbers as ' - do i=1,n_act_orb - write(6,*) i,eigval(i) - end do + write(6,*) ' occupation numbers ' + do i=1,mo_num + write(6,*) i,occnum(i) + end do - if (bavard) then -! - -integer :: nmx -real*8 :: xmx - do i=1,n_act_orb -! largest element of the eigenvector should be positive - xmx=0.D0 - nmx=0 - do j=1,n_act_orb - if (abs(natOrbsCI(j,i)).gt.xmx) then - nmx=j - xmx=abs(natOrbsCI(j,i)) - end if - end do - xmx=sign(1.D0,natOrbsCI(nmx,i)) - do j=1,n_act_orb - natOrbsCI(j,i)*=xmx - end do - - - write(6,*) ' Eigenvector No ',i - write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb) - end do - end if - - do i=1,n_act_orb - do j=1,n_act_orb - D0tu(i,j)=0.D0 - end do -! fill occupation numbers in descending order - D0tu(i,i)=eigval(n_act_orb-i+1) - end do -! -! 4-index transformation of 2part matrices -! -! index per index -! first quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=P0tuvx(q,j,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx(p,j,k,l)=d(p) - end do - end do - end do - end do -! 2nd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=P0tuvx(j,q,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx(j,p,k,l)=d(p) - end do - end do - end do - end do -! 3rd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=P0tuvx(j,k,q,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx(j,k,p,l)=d(p) - end do - end do - end do - end do -! 4th quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=P0tuvx(j,k,l,q)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - P0tuvx(j,k,l,p)=d(p) - end do - end do - end do - end do - write(6,*) ' transformed P0tuvx ' -! -! one-electron integrals -! - do i=1,mo_num - do j=1,mo_num - onetrf(i,j)=mo_one_e_integrals(i,j) - end do - end do -! 1st half-trf - do j=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=onetrf(list_act(q),j)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - onetrf(list_act(p),j)=d(p) - end do - end do -! 2nd half-trf - do j=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=onetrf(j,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - onetrf(j,list_act(p))=d(p) - end do - end do - write(6,*) ' transformed onetrf ' -! -! Orbitals -! - do j=1,ao_num - do i=1,mo_num - NatOrbsFCI(j,i)=mo_coef(j,i) - end do - end do - - do j=1,ao_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - NatOrbsFCI(j,list_act(p))=d(p) - end do - end do - write(6,*) ' transformed orbitals ' -! -! now the bielectronic integrals -! -!!$ write(6,*) ' before the transformation ' -!!$integer :: kk,ll,ii,jj -!!$real*8 :: h1,h2,h3 -!!$ do i=1,n_act_orb -!!$ ii=list_act(i) -!!$ do j=1,n_act_orb -!!$ jj=list_act(j) -!!$ do k=1,n_act_orb -!!$ kk=list_act(k) -!!$ do l=1,n_act_orb -!!$ ll=list_act(l) -!!$ h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb) -!!$ h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll) -!!$ h3=bielecCItmp(i,j,k,ll) -!!$ if ((h1.ne.h2).or.(h1.ne.h3)) then -!!$ write(6,9901) i,j,k,l,h1,h2,h3 -!!$9901 format(' aie ',4i4,3E20.12) -!!$9902 format('correct',4i4,3E20.12) -!!$ else -!!$ write(6,9902) i,j,k,l,h1,h2,h3 -!!$ end if -!!$ end do -!!$ end do -!!$ end do -!!$ end do - - do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d1(p)=0.D0 - d2(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d1(pp)+=bielec_PQxxtmp(list_act(q),j,k,l)*natorbsCI(q,p) - d2(pp)+=bielec_PxxQtmp(list_act(q),k,l,j)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxxtmp(list_act(p),j,k,l)=d1(p) - bielec_PxxQtmp(list_act(p),k,l,j)=d2(p) - end do - end do - end do - end do -! 2nd quarter - do j=1,mo_num - do k=1,n_core_orb+n_act_orb - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d1(p)=0.D0 - d2(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d1(pp)+=bielec_PQxxtmp(j,list_act(q),k,l)*natorbsCI(q,p) - d2(pp)+=bielec_PxxQtmp(j,k,l,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxxtmp(j,list_act(p),k,l)=d1(p) - bielec_PxxQtmp(j,k,l,list_act(p))=d2(p) - end do - end do - end do - end do -! 3rd quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d1(p)=0.D0 - d2(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d1(pp)+=bielec_PQxxtmp(j,k,n_core_orb+q,l)*natorbsCI(q,p) - d2(pp)+=bielec_PxxQtmp(j,n_core_orb+q,l,k)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxxtmp(j,k,n_core_orb+p,l)=d1(p) - bielec_PxxQtmp(j,n_core_orb+p,l,k)=d2(p) - end do - end do - end do - end do -! 4th quarter - do j=1,mo_num - do k=1,mo_num - do l=1,n_core_orb+n_act_orb - do p=1,n_act_orb - d1(p)=0.D0 - d2(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d1(pp)+=bielec_PQxxtmp(j,k,l,n_core_orb+q)*natorbsCI(q,p) - d2(pp)+=bielec_PxxQtmp(j,l,n_core_orb+q,k)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielec_PQxxtmp(j,k,l,n_core_orb+p)=d1(p) - bielec_PxxQtmp(j,l,n_core_orb+p,k)=d2(p) - end do - end do - end do - end do - write(6,*) ' transformed PQxx and PxxQ ' -! -! and finally the bielecCI integrals -! - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCItmp(q,j,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCItmp(p,j,k,l)=d(p) - end do - end do - end do - end do -! 2nd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCItmp(j,q,k,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCItmp(j,p,k,l)=d(p) - end do - end do - end do - end do -! 3rd quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,mo_num - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCItmp(j,k,q,l)*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCItmp(j,k,p,l)=d(p) - end do - end do - end do - end do -! 4th quarter - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - do p=1,n_act_orb - d(p)=0.D0 - end do - do p=1,n_act_orb - pp=n_act_orb-p+1 - do q=1,n_act_orb - d(pp)+=bielecCItmp(j,k,l,list_act(q))*natorbsCI(q,p) - end do - end do - do p=1,n_act_orb - bielecCItmp(j,k,l,list_act(p))=d(p) - end do - end do - end do - end do - write(6,*) ' transformed tuvP ' -! -! that's all -! -!!$ -!!$! test coherence of the bielectronic integals -!!$! PQxx = PxxQ = tuvP for some of the indices -!!$ write(6,*) ' after the transformation ' -!!$ do i=1,n_act_orb -!!$ ii=list_act(i) -!!$ do j=1,n_act_orb -!!$ jj=list_act(j) -!!$ do k=1,n_act_orb -!!$ kk=list_act(k) -!!$ do l=1,n_act_orb -!!$ ll=list_act(l) -!!$ h1=bielec_PQxxtmp(ii,jj,k+n_core_orb,l+n_core_orb) -!!$ h2=bielec_PxxQtmp(ii,j+n_core_orb,k+n_core_orb,ll) -!!$ h3=bielecCItmp(i,j,k,ll) -!!$ if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then -!!$ write(6,9901) i,j,k,l,h1,h1-h2,h1-h3 -!!$ else -!!$ write(6,9902) i,j,k,l,h1,h2,h3 -!!$ end if -!!$ end do -!!$ end do -!!$ end do -!!$ end do - -! we recalculate total energies - write(6,*) - write(6,*) ' recalculating energies after the transformation ' - write(6,*) - write(6,*) - real*8 :: e_one_all - real*8 :: e_two_all - integer :: ii - integer :: jj - integer :: t3 - integer :: tt - integer :: u3 - integer :: uu - integer :: v - integer :: v3 - integer :: vv - integer :: x - integer :: x3 - integer :: xx - - e_one_all=0.D0 - e_two_all=0.D0 - do i=1,n_core_orb - ii=list_core(i) - e_one_all+=2.D0*onetrf(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i) - end do - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_orb - e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) & - -bielec_PQxxtmp(tt,ii,i,u3)) - end do - end do - end do - do t=1,n_act_orb - tt=list_act(t) - do u=1,n_act_orb - uu=list_act(u) - e_one_all+=D0tu(t,u)*onetrf(tt,uu) - do v=1,n_act_orb - v3=v+n_core_orb - do x=1,n_act_orb - x3=x+n_core_orb - e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3) - end do - end do - end do - end do - write(6,*) ' e_one_all = ',e_one_all - write(6,*) ' e_two_all = ',e_two_all - ecore =nuclear_repulsion - ecore_bis=nuclear_repulsion - do i=1,n_core_orb - ii=list_core(i) - ecore +=2.D0*onetrf(ii,ii) - ecore_bis+=2.D0*onetrf(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - ecore +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i) - ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii) - end do - end do - eone =0.D0 - eone_bis=0.D0 - etwo =0.D0 - etwo_bis=0.D0 - etwo_ter=0.D0 - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_orb - eone +=D0tu(t,u)*onetrf(tt,uu) - eone_bis+=D0tu(t,u)*onetrf(tt,uu) - do i=1,n_core_orb - ii=list_core(i) - eone +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) & - -bielec_PQxxtmp(tt,ii,i,u3)) - eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) & - -bielec_PxxQtmp(tt,i,i,uu)) - end do - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_orb - do x=1,n_act_orb - xx=list_act(x) - x3=x+n_core_orb -real*8 :: h1,h2,h3 - h1=bielec_PQxxtmp(tt,uu,v3,x3) - h2=bielec_PxxQtmp(tt,u3,v3,xx) - h3=bielecCItmp(t,u,v,xx) - etwo +=P0tuvx(t,u,v,x)*h1 - etwo_bis+=P0tuvx(t,u,v,x)*h2 - etwo_ter+=P0tuvx(t,u,v,x)*h3 - if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then - write(6,9901) t,u,v,x,h1,h2,h3 -9901 format('aie: ',4I4,3E20.12) - end if - end do - end do - end do - end do - - write(6,*) ' energy contributions ' - write(6,*) ' core energy = ',ecore,' using PQxx integrals ' - write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' - write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' - write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' - write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' - write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' - write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' - write(6,*) ' ----------------------------------------- ' - write(6,*) ' sum of all = ',eone+etwo+ecore - write(6,*) - - end subroutine trf_to_natorb - - BEGIN_PROVIDER [real*8, onetrf, (mo_num,mo_num)] -&BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)] END_PROVIDER + + + BEGIN_PROVIDER [ real*8, natorbsCI, (n_act_orb,n_act_orb) ] +&BEGIN_PROVIDER [ real*8, occ_act, (n_act_orb) ] + implicit none + BEGIN_DOC + ! Natural orbitals of CI + END_DOC + integer :: i, j + + call lapack_diag(occ_act,natorbsCI,D0tu,n_act_orb,n_act_orb) + + write(6,*) ' found occupation numbers as ' + do i=1,n_act_orb + write(6,*) i,occ_act(i) + end do + + if (bavard) then + ! + + integer :: nmx + real*8 :: xmx + do i=1,n_act_orb + ! largest element of the eigenvector should be positive + xmx=0.D0 + nmx=0 + do j=1,n_act_orb + if (abs(natOrbsCI(j,i)).gt.xmx) then + nmx=j + xmx=abs(natOrbsCI(j,i)) + end if + end do + xmx=sign(1.D0,natOrbsCI(nmx,i)) + do j=1,n_act_orb + natOrbsCI(j,i)*=xmx + end do + + write(6,*) ' Eigenvector No ',i + write(6,'(5(I3,F12.5))') (j,natOrbsCI(j,i),j=1,n_act_orb) + end do + end if + +END_PROVIDER + + +BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] + implicit none + BEGIN_DOC + ! 4-index transformation of 2part matrices + END_DOC + integer :: i,j,k,l,p,q,pp + real*8 :: d(n_act_orb) + + ! index per index + ! first quarter + P0tuvx_no(:,:,:,:) = P0tuvx(:,:,:,:) + + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=P0tuvx_no(q,j,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(p,j,k,l)=d(p) + end do + end do + end do + end do + ! 2nd quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=P0tuvx_no(j,q,k,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,p,k,l)=d(p) + end do + end do + end do + end do + ! 3rd quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=P0tuvx_no(j,k,q,l)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,k,p,l)=d(p) + end do + end do + end do + end do + ! 4th quarter + do j=1,n_act_orb + do k=1,n_act_orb + do l=1,n_act_orb + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=P0tuvx_no(j,k,l,q)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + P0tuvx_no(j,k,l,p)=d(p) + end do + end do + end do + end do + write(6,*) ' transformed P0tuvx ' + +END_PROVIDER + + + +BEGIN_PROVIDER [real*8, onetrf, (mo_num,mo_num)] + implicit none + BEGIN_DOC + ! Transformed one-e integrals + END_DOC + integer :: i,j, p, pp, q + real*8 :: d(n_act_orb) + onetrf(:,:)=mo_one_e_integrals(:,:) + + ! 1st half-trf + do j=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=onetrf(list_act(q),j)*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + onetrf(list_act(p),j)=d(p) + end do + end do + + ! 2nd half-trf + do j=1,mo_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=onetrf(j,list_act(q))*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + onetrf(j,list_act(p))=d(p) + end do + end do + write(6,*) ' transformed onetrf ' +END_PROVIDER + + +BEGIN_PROVIDER [real*8, NatOrbsFCI, (ao_num,mo_num)] + implicit none + BEGIN_DOC +! FCI natural orbitals + END_DOC + integer :: i,j, p, pp, q + real*8 :: d(n_act_orb) + + NatOrbsFCI(:,:)=mo_coef(:,:) + + do j=1,ao_num + do p=1,n_act_orb + d(p)=0.D0 + end do + do p=1,n_act_orb + pp=n_act_orb-p+1 + do q=1,n_act_orb + d(pp)+=NatOrbsFCI(j,list_act(q))*natorbsCI(q,p) + end do + end do + do p=1,n_act_orb + NatOrbsFCI(j,list_act(p))=d(p) + end do + end do + write(6,*) ' transformed orbitals ' +END_PROVIDER + + + + + + +subroutine trf_to_natorb() + implicit none + BEGIN_DOC + ! save the diagonal somewhere, in inverse order + ! 4-index-transform the 2-particle density matrix over active orbitals + ! correct the bielectronic integrals + ! correct the monoelectronic integrals + ! put integrals on file, as well orbitals, and the density matrices + ! + END_DOC + integer :: i,j,k,l,t,u,p,q,pp + real*8 :: d(n_act_orb),d1(n_act_orb),d2(n_act_orb) + + ! we recalculate total energies + write(6,*) + write(6,*) ' recalculating energies after the transformation ' + write(6,*) + write(6,*) + real*8 :: e_one_all + real*8 :: e_two_all + integer :: ii + integer :: jj + integer :: t3 + integer :: tt + integer :: u3 + integer :: uu + integer :: v + integer :: v3 + integer :: vv + integer :: x + integer :: x3 + integer :: xx + + e_one_all=0.D0 + e_two_all=0.D0 + do i=1,n_core_orb + ii=list_core(i) + e_one_all+=2.D0*onetrf(ii,ii) + do j=1,n_core_orb + jj=list_core(j) + e_two_all+=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i) + end do + do t=1,n_act_orb + tt=list_act(t) + t3=t+n_core_orb + e_two_all += occnum(list_act(t)) * & + (2.d0*bielec_PQxx_no(tt,tt,i,i) - bielec_PQxx_no(tt,ii,i,t3)) + end do + end do + + + + do t=1,n_act_orb + tt=list_act(t) + e_one_all += occnum(list_act(t))*onetrf(tt,tt) + do u=1,n_act_orb + uu=list_act(u) + do v=1,n_act_orb + v3=v+n_core_orb + do x=1,n_act_orb + x3=x+n_core_orb + e_two_all +=P0tuvx_no(t,u,v,x)*bielec_PQxx_no(tt,uu,v3,x3) + end do + end do + end do + end do + write(6,*) ' e_one_all = ',e_one_all + write(6,*) ' e_two_all = ',e_two_all + ecore =nuclear_repulsion + ecore_bis=nuclear_repulsion + do i=1,n_core_orb + ii=list_core(i) + ecore +=2.D0*onetrf(ii,ii) + ecore_bis+=2.D0*onetrf(ii,ii) + do j=1,n_core_orb + jj=list_core(j) + ecore +=2.D0*bielec_PQxx_no(ii,ii,j,j)-bielec_PQxx_no(ii,jj,j,i) + ecore_bis+=2.D0*bielec_PxxQ_no(ii,i,j,jj)-bielec_PxxQ_no(ii,j,j,ii) + end do + end do + eone =0.D0 + eone_bis=0.D0 + etwo =0.D0 + etwo_bis=0.D0 + etwo_ter=0.D0 + do t=1,n_act_orb + tt=list_act(t) + t3=t+n_core_orb + eone += occnum(list_act(t))*onetrf(tt,tt) + eone_bis += occnum(list_act(t))*onetrf(tt,tt) + do i=1,n_core_orb + ii=list_core(i) + eone += occnum(list_act(t)) * & + (2.D0*bielec_PQxx_no(tt,tt,i,i ) - bielec_PQxx_no(tt,ii,i,t3)) + eone_bis += occnum(list_act(t)) * & + (2.D0*bielec_PxxQ_no(tt,t3,i,ii) - bielec_PxxQ_no(tt,i ,i,tt)) + end do + do u=1,n_act_orb + uu=list_act(u) + u3=u+n_core_orb + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb + xx=list_act(x) + x3=x+n_core_orb + real*8 :: h1,h2,h3 + h1=bielec_PQxx_no(tt,uu,v3,x3) + h2=bielec_PxxQ_no(tt,u3,v3,xx) + h3=bielecCI_no(t,u,v,xx) + etwo +=P0tuvx_no(t,u,v,x)*h1 + etwo_bis+=P0tuvx_no(t,u,v,x)*h2 + etwo_ter+=P0tuvx_no(t,u,v,x)*h3 + if ((abs(h1-h2).gt.1.D-14).or.(abs(h1-h3).gt.1.D-14)) then + write(6,9901) t,u,v,x,h1,h2,h3 + 9901 format('aie: ',4I4,3E20.12) + end if + end do + end do + end do + end do + + write(6,*) ' energy contributions ' + write(6,*) ' core energy = ',ecore,' using PQxx integrals ' + write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' + write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' + write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' + write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' + write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' + write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' + write(6,*) ' ----------------------------------------- ' + write(6,*) ' sum of all = ',eone+etwo+ecore + write(6,*) + SOFT_TOUCH ecore ecore_bis eone eone_bis etwo etwo_bis etwo_ter + +end subroutine trf_to_natorb + diff --git a/src/casscf/natorb_casscf.irp.f b/src/casscf/natorb_casscf.irp.f deleted file mode 100644 index 0a818a34..00000000 --- a/src/casscf/natorb_casscf.irp.f +++ /dev/null @@ -1,65 +0,0 @@ -! -*- F90 -*- -BEGIN_PROVIDER [real*8, occnum, (mo_num)] - implicit none - integer :: i,kk,j - logical :: lread - real*8 :: rdum - do i=1,mo_num - occnum(i)=0.D0 - end do - do i=1,n_core_orb - occnum(list_core(i))=2.D0 - end do - - open(unit=12,file='D0tu.dat',form='formatted',status='old') - lread=.true. - do while (lread) - read(12,*,iostat=kk) i,j,rdum - if (kk.ne.0) then - lread=.false. - else - if (i.eq.j) then - occnum(list_act(i))=rdum - else - write(6,*) ' WARNING - no natural orbitals !' - write(6,*) i,j,rdum - end if - end if - end do - close(12) - write(6,*) ' read occupation numbers ' - do i=1,mo_num - write(6,*) i,occnum(i) - end do - -END_PROVIDER - -BEGIN_PROVIDER [real*8, P0tuvx_no, (n_act_orb,n_act_orb,n_act_orb,n_act_orb)] - implicit none - integer :: i,j,k,l,kk - real*8 :: rdum - logical :: lread - - do i=1,n_act_orb - do j=1,n_act_orb - do k=1,n_act_orb - do l=1,n_act_orb - P0tuvx_no(l,k,j,i)=0.D0 - end do - end do - end do - end do - - open(unit=12,file='P0tuvx.dat',form='formatted',status='old') - lread=.true. - do while (lread) - read(12,*,iostat=kk) i,j,k,l,rdum - if (kk.ne.0) then - lread=.false. - else - P0tuvx_no(i,j,k,l)=rdum - end if - end do - close(12) - write(6,*) ' read the 2-particle density matrix ' -END_PROVIDER diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f index 8734006e..780cd543 100644 --- a/src/casscf/tot_en.irp.f +++ b/src/casscf/tot_en.irp.f @@ -1,4 +1,3 @@ -! -*- F90 -*- BEGIN_PROVIDER [real*8, etwo] &BEGIN_PROVIDER [real*8, eone] &BEGIN_PROVIDER [real*8, eone_bis] @@ -6,117 +5,117 @@ &BEGIN_PROVIDER [real*8, etwo_ter] &BEGIN_PROVIDER [real*8, ecore] &BEGIN_PROVIDER [real*8, ecore_bis] - implicit none - integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3 -real*8 :: e_one_all,e_two_all - e_one_all=0.D0 - e_two_all=0.D0 - do i=1,n_core_orb - ii=list_core(i) - e_one_all+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - e_two_all+=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i) - end do - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - do u=1,n_act_orb - uu=list_act(u) - u3=u+n_core_orb - e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) & - -bielec_PQxxtmp(tt,ii,i,u3)) - end do - end do - end do - do t=1,n_act_orb - tt=list_act(t) - do u=1,n_act_orb - uu=list_act(u) - e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu) - do v=1,n_act_orb - v3=v+n_core_orb - do x=1,n_act_orb - x3=x+n_core_orb - e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxxtmp(tt,uu,v3,x3) - end do - end do - end do - end do - write(6,*) ' e_one_all = ',e_one_all - write(6,*) ' e_two_all = ',e_two_all - ecore =nuclear_repulsion - ecore_bis=nuclear_repulsion - do i=1,n_core_orb - ii=list_core(i) - ecore +=2.D0*mo_one_e_integrals(ii,ii) - ecore_bis+=2.D0*mo_one_e_integrals(ii,ii) - do j=1,n_core_orb - jj=list_core(j) - ecore +=2.D0*bielec_PQxxtmp(ii,ii,j,j)-bielec_PQxxtmp(ii,jj,j,i) - ecore_bis+=2.D0*bielec_PxxQtmp(ii,i,j,jj)-bielec_PxxQtmp(ii,j,j,ii) - end do - end do - eone =0.D0 - eone_bis=0.D0 - etwo =0.D0 - etwo_bis=0.D0 - etwo_ter=0.D0 - do t=1,n_act_orb - tt=list_act(t) - t3=t+n_core_orb - do u=1,n_act_orb + implicit none + integer :: t,u,v,x,i,ii,tt,uu,vv,xx,j,jj,t3,u3,v3,x3 + real*8 :: e_one_all,e_two_all + e_one_all=0.D0 + e_two_all=0.D0 + do i=1,n_core_orb + ii=list_core(i) + e_one_all+=2.D0*mo_one_e_integrals(ii,ii) + do j=1,n_core_orb + jj=list_core(j) + e_two_all+=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) + end do + do t=1,n_act_orb + tt=list_act(t) + t3=t+n_core_orb + do u=1,n_act_orb uu=list_act(u) u3=u+n_core_orb - eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu) - eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu) - do i=1,n_core_orb - ii=list_core(i) - eone +=D0tu(t,u)*(2.D0*bielec_PQxxtmp(tt,uu,i,i) & - -bielec_PQxxtmp(tt,ii,i,u3)) - eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQtmp(tt,u3,i,ii) & - -bielec_PxxQtmp(tt,i,i,uu)) + e_two_all+=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & + -bielec_PQxx(tt,ii,i,u3)) + end do + end do + end do + do t=1,n_act_orb + tt=list_act(t) + do u=1,n_act_orb + uu=list_act(u) + e_one_all+=D0tu(t,u)*mo_one_e_integrals(tt,uu) + do v=1,n_act_orb + v3=v+n_core_orb + do x=1,n_act_orb + x3=x+n_core_orb + e_two_all +=P0tuvx(t,u,v,x)*bielec_PQxx(tt,uu,v3,x3) end do - do v=1,n_act_orb - vv=list_act(v) - v3=v+n_core_orb - do x=1,n_act_orb + end do + end do + end do + write(6,*) ' e_one_all = ',e_one_all + write(6,*) ' e_two_all = ',e_two_all + ecore =nuclear_repulsion + ecore_bis=nuclear_repulsion + do i=1,n_core_orb + ii=list_core(i) + ecore +=2.D0*mo_one_e_integrals(ii,ii) + ecore_bis+=2.D0*mo_one_e_integrals(ii,ii) + do j=1,n_core_orb + jj=list_core(j) + ecore +=2.D0*bielec_PQxx(ii,ii,j,j)-bielec_PQxx(ii,jj,j,i) + ecore_bis+=2.D0*bielec_PxxQ(ii,i,j,jj)-bielec_PxxQ(ii,j,j,ii) + end do + end do + eone =0.D0 + eone_bis=0.D0 + etwo =0.D0 + etwo_bis=0.D0 + etwo_ter=0.D0 + do t=1,n_act_orb + tt=list_act(t) + t3=t+n_core_orb + do u=1,n_act_orb + uu=list_act(u) + u3=u+n_core_orb + eone +=D0tu(t,u)*mo_one_e_integrals(tt,uu) + eone_bis+=D0tu(t,u)*mo_one_e_integrals(tt,uu) + do i=1,n_core_orb + ii=list_core(i) + eone +=D0tu(t,u)*(2.D0*bielec_PQxx(tt,uu,i,i) & + -bielec_PQxx(tt,ii,i,u3)) + eone_bis+=D0tu(t,u)*(2.D0*bielec_PxxQ(tt,u3,i,ii) & + -bielec_PxxQ(tt,i,i,uu)) + end do + do v=1,n_act_orb + vv=list_act(v) + v3=v+n_core_orb + do x=1,n_act_orb xx=list_act(x) x3=x+n_core_orb -real*8 :: h1,h2,h3 - h1=bielec_PQxxtmp(tt,uu,v3,x3) - h2=bielec_PxxQtmp(tt,u3,v3,xx) - h3=bielecCItmp(t,u,v,xx) + real*8 :: h1,h2,h3 + h1=bielec_PQxx(tt,uu,v3,x3) + h2=bielec_PxxQ(tt,u3,v3,xx) + h3=bielecCI(t,u,v,xx) etwo +=P0tuvx(t,u,v,x)*h1 etwo_bis+=P0tuvx(t,u,v,x)*h2 etwo_ter+=P0tuvx(t,u,v,x)*h3 if ((h1.ne.h2).or.(h1.ne.h3)) then - write(6,9901) t,u,v,x,h1,h2,h3 -9901 format('aie: ',4I4,3E20.12) + write(6,9901) t,u,v,x,h1,h2,h3 + 9901 format('aie: ',4I4,3E20.12) end if - end do end do - end do end do - - write(6,*) ' energy contributions ' - write(6,*) ' core energy = ',ecore,' using PQxx integrals ' - write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' - write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' - write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' - write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' - write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' - write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' - write(6,*) ' ----------------------------------------- ' - write(6,*) ' sum of all = ',eone+etwo+ecore - write(6,*) - write(6,*) ' nuclear (qp) = ',nuclear_repulsion - write(6,*) ' core energy (qp) = ',core_energy - write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1) - write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1) - write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1) - write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1) - + end do + end do + + write(6,*) ' energy contributions ' + write(6,*) ' core energy = ',ecore,' using PQxx integrals ' + write(6,*) ' core energy (bis) = ',ecore,' using PxxQ integrals ' + write(6,*) ' 1el energy = ',eone ,' using PQxx integrals ' + write(6,*) ' 1el energy (bis) = ',eone ,' using PxxQ integrals ' + write(6,*) ' 2el energy = ',etwo ,' using PQxx integrals ' + write(6,*) ' 2el energy (bis) = ',etwo_bis,' using PxxQ integrals ' + write(6,*) ' 2el energy (ter) = ',etwo_ter,' using tuvP integrals ' + write(6,*) ' ----------------------------------------- ' + write(6,*) ' sum of all = ',eone+etwo+ecore + write(6,*) + write(6,*) ' nuclear (qp) = ',nuclear_repulsion + write(6,*) ' core energy (qp) = ',core_energy + write(6,*) ' 1el energy (qp) = ',psi_energy_h_core(1) + write(6,*) ' 2el energy (qp) = ',psi_energy_two_e(1) + write(6,*) ' nuc + 1 + 2 (qp) = ',nuclear_repulsion+psi_energy_h_core(1)+psi_energy_two_e(1) + write(6,*) ' <0|H|0> (qp) = ',psi_energy_with_nucl_rep(1) + END_PROVIDER - - + +