From d29f82c0800de5baf37047b010b8f868cf630cf5 Mon Sep 17 00:00:00 2001 From: Anthony Scemama Date: Mon, 24 Jun 2019 16:42:16 +0200 Subject: [PATCH] CAS-CI and wdens merged --- src/casscf/bavard.irp.f | 6 + src/casscf/bielec_create.irp.f | 118 +++++++ src/casscf/casscf.irp.f | 4 +- src/casscf/densities.irp.f | 177 +++++++++++ src/casscf/det_manip.irp.f | 131 ++++++++ src/casscf/driver_wdens.irp.f | 154 +++++++++ src/casscf/natorb.irp.f | 548 +++++++++++++++++++++++++++++++++ src/casscf/tot_en.irp.f | 122 ++++++++ 8 files changed, 1259 insertions(+), 1 deletion(-) create mode 100644 src/casscf/bavard.irp.f create mode 100644 src/casscf/bielec_create.irp.f create mode 100644 src/casscf/densities.irp.f create mode 100644 src/casscf/det_manip.irp.f create mode 100644 src/casscf/driver_wdens.irp.f create mode 100644 src/casscf/natorb.irp.f create mode 100644 src/casscf/tot_en.irp.f diff --git a/src/casscf/bavard.irp.f b/src/casscf/bavard.irp.f new file mode 100644 index 00000000..de71a346 --- /dev/null +++ b/src/casscf/bavard.irp.f @@ -0,0 +1,6 @@ +! -*- F90 -*- +BEGIN_PROVIDER [logical, bavard] + bavard=.true. + bavard=.false. +END_PROVIDER + diff --git a/src/casscf/bielec_create.irp.f b/src/casscf/bielec_create.irp.f new file mode 100644 index 00000000..7e6d16c8 --- /dev/null +++ b/src/casscf/bielec_create.irp.f @@ -0,0 +1,118 @@ +! -*- 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/casscf.irp.f b/src/casscf/casscf.irp.f index 28f56069..c08dd032 100644 --- a/src/casscf/casscf.irp.f +++ b/src/casscf/casscf.irp.f @@ -1,4 +1,4 @@ -program casscf_new +program casscf implicit none BEGIN_DOC ! TODO : Put the documentation of the program here @@ -11,4 +11,6 @@ end subroutine run implicit none call run_cipsi + call driver_wdens + end diff --git a/src/casscf/densities.irp.f b/src/casscf/densities.irp.f new file mode 100644 index 00000000..77f5593e --- /dev/null +++ b/src/casscf/densities.irp.f @@ -0,0 +1,177 @@ +! -*- F90 -*- +use bitmasks ! you need to include the bitmasks_module.f90 features + + 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 & + ,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 + D0tu(t,u)+=term +! 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 + D0tu(t,u)+=term + 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 + 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_PROVIDER diff --git a/src/casscf/det_manip.irp.f b/src/casscf/det_manip.irp.f new file mode 100644 index 00000000..c8e6c08a --- /dev/null +++ b/src/casscf/det_manip.irp.f @@ -0,0 +1,131 @@ +! -*- F90 -*- +use bitmasks ! you need to include the bitmasks_module.f90 features + + 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 + 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_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 + + end subroutine do_spinfree_mono_excitation + diff --git a/src/casscf/driver_wdens.irp.f b/src/casscf/driver_wdens.irp.f new file mode 100644 index 00000000..263e8441 --- /dev/null +++ b/src/casscf/driver_wdens.irp.f @@ -0,0 +1,154 @@ + subroutine driver_wdens + implicit none + integer :: istate,p,q,r,s,indx,i,j + + + write(6,*) ' total energy = ',eone+etwo+ecore + 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) + +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 + do q=1,mo_num + do p=q,mo_num + if (abs(onetrf(p,q)).gt.1.D-12) then + write(12,'(2i6,E20.12)') p,q,onetrf(p,q) + indx+=1 + end if + end do + end do + write(6,*) ' wrote ',indx,' mono-electronic integrals' + 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 + write(6,*) ' Orbital No ',i + write(6,'(5F14.6)') (NatOrbsFCI(j,i),j=1,mo_num) + write(6,*) + end do + + mo_label = "MCSCF" + mo_label = "Natural" + do i=1,mo_num + do j=1,ao_num + mo_coef(j,i)=NatOrbsFCI(j,i) + end do + end do + call save_mos + + write(6,*) ' ... done ' + + end + diff --git a/src/casscf/natorb.irp.f b/src/casscf/natorb.irp.f new file mode 100644 index 00000000..a903260c --- /dev/null +++ b/src/casscf/natorb.irp.f @@ -0,0 +1,548 @@ +! -*- 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) + + 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 + + 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 diff --git a/src/casscf/tot_en.irp.f b/src/casscf/tot_en.irp.f new file mode 100644 index 00000000..8734006e --- /dev/null +++ b/src/casscf/tot_en.irp.f @@ -0,0 +1,122 @@ +! -*- F90 -*- + BEGIN_PROVIDER [real*8, etwo] +&BEGIN_PROVIDER [real*8, eone] +&BEGIN_PROVIDER [real*8, eone_bis] +&BEGIN_PROVIDER [real*8, etwo_bis] +&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 + 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)) + 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 ((h1.ne.h2).or.(h1.ne.h3)) 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,*) + 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 + +