From 392b5436c2c11991b80c72d1bd393e2c59143de6 Mon Sep 17 00:00:00 2001 From: eginer Date: Thu, 3 Nov 2022 16:15:39 +0100 Subject: [PATCH] foct_tc_mo_tot is OK --- .../grid_becke_per_atom.irp.f | 4 - src/dft_utils_in_r/mo_in_r.irp.f | 27 ++++ src/tc_scf/fock_tc.irp.f | 22 ++-- src/tc_scf/fock_tc_mo_tot.irp.f | 118 ++++++++++++++++++ src/tools/print_wf.irp.f | 17 +++ 5 files changed, 173 insertions(+), 15 deletions(-) create mode 100644 src/tc_scf/fock_tc_mo_tot.irp.f diff --git a/src/becke_numerical_grid/grid_becke_per_atom.irp.f b/src/becke_numerical_grid/grid_becke_per_atom.irp.f index 6026350b..75610090 100644 --- a/src/becke_numerical_grid/grid_becke_per_atom.irp.f +++ b/src/becke_numerical_grid/grid_becke_per_atom.irp.f @@ -46,8 +46,4 @@ END_PROVIDER enddo enddo enddo - - - - END_PROVIDER diff --git a/src/dft_utils_in_r/mo_in_r.irp.f b/src/dft_utils_in_r/mo_in_r.irp.f index 0a8b4d52..2bc6fc30 100644 --- a/src/dft_utils_in_r/mo_in_r.irp.f +++ b/src/dft_utils_in_r/mo_in_r.irp.f @@ -181,3 +181,30 @@ END_PROVIDER + BEGIN_PROVIDER [double precision, mo_abs_dist_per_nuclei, (mo_num,nucl_num)] + implicit none + BEGIN_DOC + ! mo_abs_dist_per_nuclei(j,i_nucl) = + END_DOC + integer :: ipoint,i_nucl,m,j + double precision :: weight, r(3),r_nucl(3),dist + mo_abs_dist_per_nuclei = 0.d0 + do ipoint = 1, n_points_final_grid + r(1) = final_grid_points(1,ipoint) + r(2) = final_grid_points(2,ipoint) + r(3) = final_grid_points(3,ipoint) + weight = final_weight_at_r_vector(ipoint) + do i_nucl = 1, nucl_num + dist = 0.d0 + do m = 1, 3 + r_nucl(m) = r(m) - nucl_coord_transp(m,i_nucl) + r_nucl(m) *= r_nucl(m) + dist += r_nucl(m) + enddo + dist = dsqrt(dist) + do j = 1, mo_num + mo_abs_dist_per_nuclei(j,i_nucl) += weight * mos_in_r_array(j,ipoint)*mos_in_r_array(j,ipoint) * dist + enddo + enddo + enddo + END_PROVIDER diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 4907f2c8..7eba9231 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -117,17 +117,17 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ] END_PROVIDER -BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)] - implicit none - BEGIN_DOC - ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis - END_DOC - Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta) - if(three_body_h_tc) then - Fock_matrix_tc_mo_tot += fock_3_mat - endif - !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10) -END_PROVIDER +!BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num, mo_num)] +! implicit none +! BEGIN_DOC +! ! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis +! END_DOC +! Fock_matrix_tc_mo_tot = 0.5d0 * (Fock_matrix_tc_mo_alpha + Fock_matrix_tc_mo_beta) +! if(three_body_h_tc) then +! Fock_matrix_tc_mo_tot += fock_3_mat +! endif +! !call restore_symmetry(mo_num, mo_num, Fock_matrix_tc_mo_tot, mo_num, 1.d-10) +!END_PROVIDER ! --- diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f new file mode 100644 index 00000000..5eaa4726 --- /dev/null +++ b/src/tc_scf/fock_tc_mo_tot.irp.f @@ -0,0 +1,118 @@ + + BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (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_tc_mo_tot = Fock_matrix_tc_mo_alpha + else + + do j=1,elec_beta_num + ! F-K + do i=1,elec_beta_num !CC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F+K/2 + do i=elec_beta_num+1,elec_alpha_num !CA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F + do i=elec_alpha_num+1, mo_num !CV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(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_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F + do i=elec_beta_num+1,elec_alpha_num !AA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_alpha_num+1, mo_num !AV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + enddo + + do j=elec_alpha_num+1, mo_num + ! F + do i=1,elec_beta_num !VC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_beta_num+1,elec_alpha_num !VA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F+K + do i=elec_alpha_num+1,mo_num !VV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) & + + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + enddo + + endif + + do i = 1, mo_num + Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i) + 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_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + + if(no_oa_or_av_opt)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + +END_PROVIDER + diff --git a/src/tools/print_wf.irp.f b/src/tools/print_wf.irp.f index 7e51caaf..3b9f6978 100644 --- a/src/tools/print_wf.irp.f +++ b/src/tools/print_wf.irp.f @@ -17,6 +17,7 @@ program print_wf ! psi_coef_sorted are the wave function stored in the |EZFIO| directory. read_wf = .True. touch read_wf + call write_wf call routine end @@ -120,3 +121,19 @@ subroutine routine print*,'L2 norm of pert beta = ',norm_mono_b_pert_2 end + +subroutine write_wf + implicit none + character*(128) :: output + integer :: i_unit_output,getUnitAndOpen + output=trim(ezfio_filename)//'.wf' + i_unit_output = getUnitAndOpen(output,'w') + integer :: i + print*,'Writing the sorted wf' + do i = 1, N_det + write(i_unit_output,*)i,psi_coef_sorted(i,1)/psi_coef_sorted(1,1) + enddo + + +end +