From 2a386ffa419bb5f67b7b7bd075cbf3d839f5519b Mon Sep 17 00:00:00 2001 From: Kevin Gasperich Date: Mon, 27 Jan 2020 17:20:50 -0600 Subject: [PATCH] working on complex HF --- src/scf_utils/diagonalize_fock_complex.irp.f | 53 +++++++ src/scf_utils/fock_matrix.irp.f | 73 +++------ src/scf_utils/fock_matrix_complex.irp.f | 148 +++++++++++++++++++ 3 files changed, 221 insertions(+), 53 deletions(-) create mode 100644 src/scf_utils/diagonalize_fock_complex.irp.f create mode 100644 src/scf_utils/fock_matrix_complex.irp.f diff --git a/src/scf_utils/diagonalize_fock_complex.irp.f b/src/scf_utils/diagonalize_fock_complex.irp.f new file mode 100644 index 00000000..1150b773 --- /dev/null +++ b/src/scf_utils/diagonalize_fock_complex.irp.f @@ -0,0 +1,53 @@ +BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_complex, (ao_num,mo_num) ] + implicit none + BEGIN_DOC + ! Eigenvectors of the Fock matrix in the |MO| basis obtained with level shift. + END_DOC + + integer :: i,j + integer :: n + complex*16, allocatable :: F(:,:) + double precision, allocatable :: diag(:) + + + allocate( F(mo_num,mo_num) ) + allocate (diag(mo_num) ) + + do j=1,mo_num + do i=1,mo_num + F(i,j) = Fock_matrix_mo_complex(i,j) + enddo + enddo + + if(frozen_orb_scf)then + integer :: iorb,jorb + do i = 1, n_core_orb + iorb = list_core(i) + do j = 1, n_act_orb + jorb = list_act(j) + F(iorb,jorb) = (0.d0,0.d0) + F(jorb,iorb) = (0.d0,0.d0) + enddo + enddo + endif + + ! Insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F(i,i) += 0.5d0*level_shift + enddo + + do i = elec_alpha_num+1, mo_num + F(i,i) += level_shift + enddo + + n = mo_num + call lapack_diagd_diag_in_place_complex(diag,F,n,n) + + call zgemm('N','N',ao_num,mo_num,mo_num, (1.d0,0.d0), & + mo_coef_complex, size(mo_coef_complex,1), F, size(F,1), & + (0.d0,0.d0), eigenvectors_Fock_matrix_mo_complex, size(eigenvectors_Fock_matrix_mo_complex,1)) + deallocate(F, diag) + + +END_PROVIDER + diff --git a/src/scf_utils/fock_matrix.irp.f b/src/scf_utils/fock_matrix.irp.f index a49d51e5..0cecd7d4 100644 --- a/src/scf_utils/fock_matrix.irp.f +++ b/src/scf_utils/fock_matrix.irp.f @@ -110,15 +110,6 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_num,mo_num) ] endif END_PROVIDER -BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ] - implicit none - BEGIN_DOC - ! Fock matrix on the MO basis - END_DOC - call ao_to_mo_complex(Fock_matrix_ao_alpha_complex,size(Fock_matrix_ao_alpha_complex,1), & - Fock_matrix_mo_alpha_complex,size(Fock_matrix_mo_alpha_complex,1)) -END_PROVIDER - BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_num,mo_num) ] implicit none BEGIN_DOC @@ -133,16 +124,6 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_num,mo_num) ] endif END_PROVIDER -BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ] - implicit none - BEGIN_DOC - ! Fock matrix on the MO basis - END_DOC - call ao_to_mo_complex(Fock_matrix_ao_beta_complex,size(Fock_matrix_ao_beta_complex,1), & - Fock_matrix_mo_beta_complex,size(Fock_matrix_mo_beta_complex,1)) -END_PROVIDER - - BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] implicit none BEGIN_DOC @@ -169,34 +150,6 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ] endif END_PROVIDER - -BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] - implicit none - BEGIN_DOC - ! Fock matrix in AO basis set - END_DOC - - if(frozen_orb_scf)then - call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & - Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) - else - if ( (elec_alpha_num == elec_beta_num).and. & - (level_shift == 0.) ) & - then - integer :: i,j - do j=1,ao_num - do i=1,ao_num - Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(i,j) - enddo - enddo - else - call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & - Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) - endif - endif -END_PROVIDER - - BEGIN_PROVIDER [ double precision, SCF_energy ] implicit none BEGIN_DOC @@ -205,13 +158,27 @@ BEGIN_PROVIDER [ double precision, SCF_energy ] SCF_energy = nuclear_repulsion integer :: i,j - do j=1,ao_num - do i=1,ao_num - SCF_energy += 0.5d0 * ( & - (ao_one_e_integrals(i,j) + Fock_matrix_ao_alpha(i,j) ) * SCF_density_matrix_ao_alpha(i,j) +& - (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) * SCF_density_matrix_ao_beta (i,j) ) + if (is_periodic) then + complex*16 :: scf_e_tmp + scf_e_tmp = (0.d0,0.d0) + do j=1,ao_num + do i=1,ao_num + scf_e_tmp += 0.5d0 * ( & + (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_alpha_complex(i,j) ) * SCF_density_matrix_ao_alpha_complex(i,j) +& + (ao_one_e_integrals_complex(i,j) + Fock_matrix_ao_beta_complex (i,j) ) * SCF_density_matrix_ao_beta_complex (i,j) ) + enddo enddo - enddo + !TODO: add check for imaginary part? (should be zero) + SCF_energy = dble(scf_e_tmp) + else + do j=1,ao_num + do i=1,ao_num + SCF_energy += 0.5d0 * ( & + (ao_one_e_integrals(i,j) + Fock_matrix_ao_alpha(i,j) ) * SCF_density_matrix_ao_alpha(i,j) +& + (ao_one_e_integrals(i,j) + Fock_matrix_ao_beta (i,j) ) * SCF_density_matrix_ao_beta (i,j) ) + enddo + enddo + endif SCF_energy += extra_e_contrib_density END_PROVIDER diff --git a/src/scf_utils/fock_matrix_complex.irp.f b/src/scf_utils/fock_matrix_complex.irp.f new file mode 100644 index 00000000..290f9b9d --- /dev/null +++ b/src/scf_utils/fock_matrix_complex.irp.f @@ -0,0 +1,148 @@ + BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_complex, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_diag_mo_complex, (mo_num)] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is :: + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + integer :: i,j,n + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_mo_complex = Fock_matrix_mo_alpha_complex + else + + do j=1,elec_beta_num + ! F-K + do i=1,elec_beta_num !CC + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j))& + - (Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + ! F+K/2 + do i=elec_beta_num+1,elec_alpha_num !CA + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j))& + + 0.5d0*(Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + ! F + do i=elec_alpha_num+1, mo_num !CV + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j)) + enddo + enddo + + do j=elec_beta_num+1,elec_alpha_num + ! F+K/2 + do i=1,elec_beta_num !AC + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j))& + + 0.5d0*(Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + ! F + do i=elec_beta_num+1,elec_alpha_num !AA + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j)) + enddo + ! F-K/2 + do i=elec_alpha_num+1, mo_num !AV + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j))& + - 0.5d0*(Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + enddo + + do j=elec_alpha_num+1, mo_num + ! F + do i=1,elec_beta_num !VC + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j)) + enddo + ! F-K/2 + do i=elec_beta_num+1,elec_alpha_num !VA + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j))& + - 0.5d0*(Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + ! F+K + do i=elec_alpha_num+1,mo_num !VV + Fock_matrix_mo_complex(i,j) = 0.5d0*(Fock_matrix_mo_alpha_complex(i,j)+Fock_matrix_mo_beta_complex(i,j)) & + + (Fock_matrix_mo_beta_complex(i,j) - Fock_matrix_mo_alpha_complex(i,j)) + enddo + enddo + + endif + + do i = 1, mo_num + Fock_matrix_diag_mo_complex(i) = dble(Fock_matrix_mo_complex(i,i)) + if (dabs(dimag(Fock_matrix_mo_complex(i,i))) .gt. 1.0d-12) then + !stop 'diagonal elements of Fock matrix should be real' + print *, 'diagonal elements of Fock matrix should be real',i,Fock_matrix_mo_complex(i,i) + stop -1 + endif + enddo + + + if(frozen_orb_scf)then + integer :: iorb,jorb + do i = 1, n_core_orb + iorb = list_core(i) + do j = 1, n_act_orb + jorb = list_act(j) + Fock_matrix_mo_complex(iorb,jorb) = (0.d0,0.d0) + Fock_matrix_mo_complex(jorb,iorb) = (0.d0,0.d0) + enddo + enddo + endif + +END_PROVIDER + + + +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_alpha_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_alpha_complex,size(Fock_matrix_ao_alpha_complex,1), & + Fock_matrix_mo_alpha_complex,size(Fock_matrix_mo_alpha_complex,1)) +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, Fock_matrix_mo_beta_complex, (mo_num,mo_num) ] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis + END_DOC + call ao_to_mo_complex(Fock_matrix_ao_beta_complex,size(Fock_matrix_ao_beta_complex,1), & + Fock_matrix_mo_beta_complex,size(Fock_matrix_mo_beta_complex,1)) +END_PROVIDER + + +BEGIN_PROVIDER [ complex*16, Fock_matrix_ao_complex, (ao_num, ao_num) ] + implicit none + BEGIN_DOC + ! Fock matrix in AO basis set + END_DOC + + if(frozen_orb_scf)then + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + else + if ( (elec_alpha_num == elec_beta_num).and. & + (level_shift == 0.) ) & + then + integer :: i,j + do j=1,ao_num + do i=1,ao_num + Fock_matrix_ao_complex(i,j) = Fock_matrix_ao_alpha_complex(i,j) + enddo + enddo + else + call mo_to_ao_complex(Fock_matrix_mo_complex,size(Fock_matrix_mo_complex,1), & + Fock_matrix_ao_complex,size(Fock_matrix_ao_complex,1)) + endif + endif +END_PROVIDER +