foct_tc_mo_tot is OK

2024-06-26 15:12:19 +02:00 · 2022-11-03 16:15:39 +01:00 · 2022-11-03 16:15:39 +01:00 · 392b5436c2
commit 392b5436c2
parent 0bea180fc6
5 changed files with 173 additions and 15 deletions
--- 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 
--- 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) = <phi_j| |r-R_nucl| |phi_j>
+  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 
--- 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 

 ! ---

--- a/src/tc_scf/fock_tc_mo_tot.irp.f
+++ 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
+
--- 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
+