LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

added tc_scf

Browse Source
This commit is contained in:
eginer 2022-10-05 16:04:49 +02:00
parent bdce53d8b1
commit 8c51af6c5a
17 changed files with 2386 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
src/tc_scf/EZFIO.cfg Normal file
View File

@ -0,0 +1,4 @@
[bitc_energy]
type: Threshold
doc: Energy bi-tc HF
interface: ezfio

6
src/tc_scf/NEED Normal file
View File

@ -0,0 +1,6 @@
hartree_fock
bi_ortho_mos
three_body_ints
bi_ort_ints
tc_keywords
non_hermit_dav

74
src/tc_scf/combine_lr_tcscf.irp.f Normal file
View File

@ -0,0 +1,74 @@
! ---
program combine_lr_tcscf
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
bi_ortho = .True.
touch bi_ortho
call comb_orbitals()
end
! ---
subroutine comb_orbitals()
implicit none
integer :: i, m, n, nn, mm
double precision :: accu_d, accu_nd
double precision, allocatable :: R(:,:), L(:,:), Rnew(:,:), tmp(:,:), S(:,:)
n = ao_num
m = mo_num
nn = elec_alpha_num
mm = m - nn
allocate(L(n,m), R(n,m), Rnew(n,m), S(m,m))
L = mo_l_coef
R = mo_r_coef
call check_weighted_biorthog(n, m, ao_overlap, L, R, accu_d, accu_nd, S, .true.)
allocate(tmp(n,nn))
do i = 1, nn
tmp(1:n,i) = R(1:n,i)
enddo
call impose_weighted_orthog_svd(n, nn, ao_overlap, tmp)
do i = 1, nn
Rnew(1:n,i) = tmp(1:n,i)
enddo
deallocate(tmp)
allocate(tmp(n,mm))
do i = 1, mm
tmp(1:n,i) = L(1:n,i+nn)
enddo
call impose_weighted_orthog_svd(n, mm, ao_overlap, tmp)
do i = 1, mm
Rnew(1:n,i+nn) = tmp(1:n,i)
enddo
deallocate(tmp)
call check_weighted_biorthog(n, m, ao_overlap, Rnew, Rnew, accu_d, accu_nd, S, .true.)
mo_r_coef = Rnew
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
deallocate(L, R, Rnew, S)
end subroutine comb_orbitals
! ---

168
src/tc_scf/diago_bi_ort_tcfock.irp.f Normal file
View File

@ -0,0 +1,168 @@
BEGIN_PROVIDER [ double precision, fock_tc_reigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_tc_leigvec_mo, (mo_num, mo_num)]
&BEGIN_PROVIDER [ double precision, eigval_fock_tc_mo, (mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_fock_tc_eigvec_mo, (mo_num, mo_num)]
BEGIN_DOC
! EIGENVECTORS OF FOCK MATRIX ON THE MO BASIS and their OVERLAP
END_DOC
implicit none
integer :: n_real_tc
integer :: i, k, l
double precision :: accu_d, accu_nd, accu_tmp
double precision :: norm
double precision, allocatable :: eigval_right_tmp(:)
allocate( eigval_right_tmp(mo_num) )
PROVIDE Fock_matrix_tc_mo_tot
call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot &
, fock_tc_leigvec_mo, fock_tc_reigvec_mo &
, n_real_tc, eigval_right_tmp )
!if(max_ov_tc_scf)then
! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
!else
! call non_hrmt_diag_split_degen_bi_orthog( mo_num, Fock_matrix_tc_mo_tot &
! , fock_tc_leigvec_mo, fock_tc_reigvec_mo &
! , n_real_tc, eigval_right_tmp )
!endif
! if(n_real_tc .ne. mo_num)then
! print*,'n_real_tc ne mo_num ! ',n_real_tc
! stop
! endif
eigval_fock_tc_mo = eigval_right_tmp
! print*,'Eigenvalues of Fock_matrix_tc_mo_tot'
! do i = 1, mo_num
! print*, i, eigval_fock_tc_mo(i)
! enddo
! deallocate( eigval_right_tmp )
! L.T x R
call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0 &
, fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
, fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
, 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) )
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, mo_num
do k = 1, mo_num
if(i==k) then
accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
accu_d += dabs(accu_tmp )
else
accu_tmp = overlap_fock_tc_eigvec_mo(k,i)
accu_nd += accu_tmp * accu_tmp
if(dabs(overlap_fock_tc_eigvec_mo(k,i)).gt.1.d-10)then
print*,'k,i',k,i,overlap_fock_tc_eigvec_mo(k,i)
endif
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)/accu_d
if( accu_nd .gt. 1d-8 ) then
print *, ' bi-orthog failed'
print*,'accu_nd MO = ', accu_nd
print*,'overlap_fock_tc_eigvec_mo = '
do i = 1, mo_num
write(*,'(100(F16.10,X))') overlap_fock_tc_eigvec_mo(i,:)
enddo
stop
endif
if( dabs(accu_d - dble(mo_num)) .gt. 1e-7 ) then
print *, 'mo_num = ', mo_num
print *, 'accu_d MO = ', accu_d
print *, 'normalizing vectors ...'
do i = 1, mo_num
norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i)))
if( norm.gt.1e-7 ) then
do k = 1, mo_num
fock_tc_reigvec_mo(k,i) *= 1.d0/norm
fock_tc_leigvec_mo(k,i) *= 1.d0/norm
enddo
endif
enddo
call dgemm( "T", "N", mo_num, mo_num, mo_num, 1.d0 &
, fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
, fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
, 0.d0, overlap_fock_tc_eigvec_mo, size(overlap_fock_tc_eigvec_mo, 1) )
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_tc_reigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, fock_tc_leigvec_ao, (ao_num, mo_num)]
&BEGIN_PROVIDER [ double precision, overlap_fock_tc_eigvec_ao, (mo_num, mo_num) ]
BEGIN_DOC
! EIGENVECTORS OF FOCK MATRIX ON THE AO BASIS and their OVERLAP
!
! THE OVERLAP SHOULD BE THE SAME AS overlap_fock_tc_eigvec_mo
END_DOC
implicit none
integer :: i, j, k, q, p
double precision :: accu, accu_d
double precision, allocatable :: tmp(:,:)
! ! MO_R x R
call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
, mo_r_coef, size(mo_r_coef, 1) &
, fock_tc_reigvec_mo, size(fock_tc_reigvec_mo, 1) &
, 0.d0, fock_tc_reigvec_ao, size(fock_tc_reigvec_ao, 1) )
! MO_L x L
call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
, mo_l_coef, size(mo_l_coef, 1) &
, fock_tc_leigvec_mo, size(fock_tc_leigvec_mo, 1) &
, 0.d0, fock_tc_leigvec_ao, size(fock_tc_leigvec_ao, 1) )
allocate( tmp(mo_num,ao_num) )
! tmp <-- L.T x S_ao
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
, fock_tc_leigvec_ao, size(fock_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
, 0.d0, tmp, size(tmp, 1) )
! S <-- tmp x R
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
, tmp, size(tmp, 1), fock_tc_reigvec_ao, size(fock_tc_reigvec_ao, 1) &
, 0.d0, overlap_fock_tc_eigvec_ao, size(overlap_fock_tc_eigvec_ao, 1) )
deallocate( tmp )
! ---
double precision :: norm
do i = 1, mo_num
norm = 1.d0/dsqrt(dabs(overlap_fock_tc_eigvec_ao(i,i)))
do j = 1, mo_num
fock_tc_reigvec_ao(j,i) *= norm
fock_tc_leigvec_ao(j,i) *= norm
enddo
enddo
allocate( tmp(mo_num,ao_num) )
! tmp <-- L.T x S_ao
call dgemm( "T", "N", mo_num, ao_num, ao_num, 1.d0 &
, fock_tc_leigvec_ao, size(fock_tc_leigvec_ao, 1), ao_overlap, size(ao_overlap, 1) &
, 0.d0, tmp, size(tmp, 1) )
! S <-- tmp x R
call dgemm( "N", "N", mo_num, mo_num, ao_num, 1.d0 &
, tmp, size(tmp, 1), fock_tc_reigvec_ao, size(fock_tc_reigvec_ao, 1) &
, 0.d0, overlap_fock_tc_eigvec_ao, size(overlap_fock_tc_eigvec_ao, 1) )
deallocate( tmp )
END_PROVIDER

107
src/tc_scf/fock_for_right.irp.f Normal file
View File

@ -0,0 +1,107 @@
! ---
BEGIN_PROVIDER [ double precision, good_hermit_tc_fock_mat, (mo_num, mo_num)]
BEGIN_DOC
! good_hermit_tc_fock_mat = Hermitian Upper triangular Fock matrix
!
! The converged eigenvectors of such matrix yield to orthonormal vectors satisfying the left Brillouin theorem
END_DOC
implicit none
integer :: i, j
good_hermit_tc_fock_mat = Fock_matrix_tc_mo_tot
do j = 1, mo_num
do i = 1, j-1
good_hermit_tc_fock_mat(i,j) = Fock_matrix_tc_mo_tot(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, hermit_average_tc_fock_mat, (mo_num, mo_num)]
BEGIN_DOC
! hermit_average_tc_fock_mat = (F + F^\dagger)/2
END_DOC
implicit none
integer :: i, j
hermit_average_tc_fock_mat = Fock_matrix_tc_mo_tot
do j = 1, mo_num
do i = 1, mo_num
hermit_average_tc_fock_mat(i,j) = 0.5d0 * (Fock_matrix_tc_mo_tot(j,i) + Fock_matrix_tc_mo_tot(i,j))
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, grad_hermit]
implicit none
BEGIN_DOC
! square of gradient of the energy
END_DOC
if(symetric_fock_tc)then
grad_hermit = grad_hermit_average_tc_fock_mat
else
grad_hermit = grad_good_hermit_tc_fock_mat
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, grad_good_hermit_tc_fock_mat]
implicit none
BEGIN_DOC
! grad_good_hermit_tc_fock_mat = norm of gradients of the upper triangular TC fock
END_DOC
integer :: i, j
grad_good_hermit_tc_fock_mat = 0.d0
do i = 1, elec_alpha_num
do j = elec_alpha_num+1, mo_num
grad_good_hermit_tc_fock_mat += dabs(good_hermit_tc_fock_mat(i,j))
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, grad_hermit_average_tc_fock_mat]
implicit none
BEGIN_DOC
! grad_hermit_average_tc_fock_mat = norm of gradients of the upper triangular TC fock
END_DOC
integer :: i, j
grad_hermit_average_tc_fock_mat = 0.d0
do i = 1, elec_alpha_num
do j = elec_alpha_num+1, mo_num
grad_hermit_average_tc_fock_mat += dabs(hermit_average_tc_fock_mat(i,j))
enddo
enddo
END_PROVIDER
! ---
subroutine save_good_hermit_tc_eigvectors()
implicit none
integer :: sign
character*(64) :: label
logical :: output
sign = 1
label = "Canonical"
output = .False.
if(symetric_fock_tc)then
call mo_as_eigvectors_of_mo_matrix(hermit_average_tc_fock_mat, mo_num, mo_num, label, sign, output)
else
call mo_as_eigvectors_of_mo_matrix(good_hermit_tc_fock_mat, mo_num, mo_num, label, sign, output)
endif
end subroutine save_good_hermit_tc_eigvectors
! ---

133
src/tc_scf/fock_tc.irp.f Normal file
View File

@ -0,0 +1,133 @@
! ---
BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_alpha, (ao_num, ao_num)]
&BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_beta , (ao_num, ao_num)]
BEGIN_DOC
! two_e_tc_non_hermit_integral_alpha(k,i) = <k| F^tc_alpha |i>
!
! where F^tc is the two-body part of the TC Fock matrix and k,i are AO basis functions
END_DOC
implicit none
integer :: i, j, k, l
double precision :: density, density_a, density_b
two_e_tc_non_hermit_integral_alpha = 0.d0
two_e_tc_non_hermit_integral_beta = 0.d0
!! TODO :: parallelization properly done
do i = 1, ao_num
do k = 1, ao_num
!!$OMP PARALLEL &
!!$OMP DEFAULT (NONE) &
!!$OMP PRIVATE (j,l,density_a,density_b,density) &
!!$OMP SHARED (i,k,ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,ao_non_hermit_term_chemist) &
!!$OMP SHARED (two_e_tc_non_hermit_integral_alpha,two_e_tc_non_hermit_integral_beta)
!!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
do l = 1, ao_num
density_a = TCSCF_density_matrix_ao_alpha(l,j)
density_b = TCSCF_density_matrix_ao_beta (l,j)
density = density_a + density_b
! rho(l,j) * < k l| T | i j>
two_e_tc_non_hermit_integral_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i)
! rho(l,j) * < k l| T | i j>
two_e_tc_non_hermit_integral_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i)
! rho_a(l,j) * < l k| T | i j>
two_e_tc_non_hermit_integral_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i)
! rho_b(l,j) * < l k| T | i j>
two_e_tc_non_hermit_integral_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i)
enddo
enddo
!!$OMP END DO
!!$OMP END PARALLEL
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)]
implicit none
BEGIN_DOC
! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the AO basis
END_DOC
Fock_matrix_tc_ao_alpha = ao_one_e_integrals_tc_tot &
+ two_e_tc_non_hermit_integral_alpha
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_beta, (ao_num, ao_num)]
BEGIN_DOC
! Total beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
END_DOC
implicit none
Fock_matrix_tc_ao_beta = ao_one_e_integrals_tc_tot &
+ two_e_tc_non_hermit_integral_beta
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ]
implicit none
BEGIN_DOC
! Total alpha+beta TC Fock matrix : h_c + Two-e^TC terms on the AO basis
END_DOC
Fock_matrix_tc_ao_tot = 0.5d0 * (Fock_matrix_tc_ao_alpha + Fock_matrix_tc_ao_beta)
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
implicit none
BEGIN_DOC
! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis
END_DOC
if(bi_ortho)then
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
else
call ao_to_mo( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
endif
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
implicit none
BEGIN_DOC
! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis
END_DOC
if(bi_ortho)then
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
else
call ao_to_mo( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
endif
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
! ---

197
src/tc_scf/fock_three.irp.f Normal file
View File

@ -0,0 +1,197 @@
BEGIN_PROVIDER [ double precision, fock_3_mat, (mo_num, mo_num)]
implicit none
integer :: i,j
double precision :: contrib
fock_3_mat = 0.d0
if(.not.bi_ortho.and.three_body_h_tc)then
call give_fock_ia_three_e_total(1,1,contrib)
!! !$OMP PARALLEL &
!! !$OMP DEFAULT (NONE) &
!! !$OMP PRIVATE (i,j,m,integral) &
!! !$OMP SHARED (mo_num,three_body_3_index)
!! !$OMP DO SCHEDULE (guided) COLLAPSE(3)
do i = 1, mo_num
do j = 1, mo_num
call give_fock_ia_three_e_total(j,i,contrib)
fock_3_mat(j,i) = -contrib
enddo
enddo
!! !$OMP END DO
!! !$OMP END PARALLEL
!! do i = 1, mo_num
!! do j = 1, i-1
!! mat_three(j,i) = mat_three(i,j)
!! enddo
!! enddo
endif
END_PROVIDER
subroutine give_fock_ia_three_e_total(i,a,contrib)
implicit none
BEGIN_DOC
! contrib is the TOTAL (same spins / opposite spins) contribution from the three body term to the Fock operator
!
END_DOC
integer, intent(in) :: i,a
double precision, intent(out) :: contrib
double precision :: int_1, int_2, int_3
double precision :: mos_i, mos_a, w_ia
double precision :: mos_ia, weight
integer :: mm, ipoint,k,l
int_1 = 0.d0
int_2 = 0.d0
int_3 = 0.d0
do mm = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
mos_i = mos_in_r_array_transp(ipoint,i)
mos_a = mos_in_r_array_transp(ipoint,a)
mos_ia = mos_a * mos_i
w_ia = x_W_ij_erf_rk(ipoint,mm,i,a)
int_1 += weight * fock_3_w_kk_sum(ipoint,mm) * (4.d0 * fock_3_rho_beta(ipoint) * w_ia &
+ 2.0d0 * mos_ia * fock_3_w_kk_sum(ipoint,mm) &
- 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,i) * mos_a &
- 2.0d0 * fock_3_w_ki_mos_k(ipoint,mm,a) * mos_i )
int_2 += weight * (-1.d0) * ( 2.0d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * w_ia &
+ 2.0d0 * fock_3_rho_beta(ipoint) * fock_3_w_ki_wk_a(ipoint,mm,i,a) &
+ 1.0d0 * mos_ia * fock_3_trace_w_tilde(ipoint,mm) )
int_3 += weight * 1.d0 * (fock_3_w_kl_wla_phi_k(ipoint,mm,i) * mos_a + fock_3_w_kl_wla_phi_k(ipoint,mm,a) * mos_i &
+fock_3_w_ki_mos_k(ipoint,mm,i) * fock_3_w_ki_mos_k(ipoint,mm,a) )
enddo
enddo
contrib = int_1 + int_2 + int_3
end
BEGIN_PROVIDER [double precision, diag_three_elem_hf]
implicit none
integer :: i,j,k,ipoint,mm
double precision :: contrib,weight,four_third,one_third,two_third,exchange_int_231
if(.not.bi_ortho)then
if(three_body_h_tc)then
one_third = 1.d0/3.d0
two_third = 2.d0/3.d0
four_third = 4.d0/3.d0
diag_three_elem_hf = 0.d0
do i = 1, elec_beta_num
do j = 1, elec_beta_num
do k = 1, elec_beta_num
call give_integrals_3_body(k,j,i,j,i,k,exchange_int_231)
diag_three_elem_hf += two_third * exchange_int_231
enddo
enddo
enddo
do mm = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
contrib = 3.d0 * fock_3_w_kk_sum(ipoint,mm) * fock_3_rho_beta(ipoint) * fock_3_w_kk_sum(ipoint,mm) &
-2.d0 * fock_3_w_kl_mo_k_mo_l(ipoint,mm) * fock_3_w_kk_sum(ipoint,mm) &
-1.d0 * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm)
contrib *= four_third
contrib += -two_third * fock_3_rho_beta(ipoint) * fock_3_w_kl_w_kl(ipoint,mm) &
- four_third * fock_3_w_kk_sum(ipoint,mm) * fock_3_w_kl_mo_k_mo_l(ipoint,mm)
diag_three_elem_hf += weight * contrib
enddo
enddo
diag_three_elem_hf = - diag_three_elem_hf
else
diag_three_elem_hf = 0.D0
endif
else
diag_three_elem_hf = 0.D0
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh, (mo_num, mo_num)]
implicit none
integer :: h,p,i,j
double precision :: direct_int, exch_int, exchange_int_231, exchange_int_312
double precision :: exchange_int_23, exchange_int_12, exchange_int_13
fock_3_mat_a_op_sh = 0.d0
do h = 1, mo_num
do p = 1, mo_num
!F_a^{ab}(h,p)
do i = 1, elec_beta_num ! beta
do j = elec_beta_num+1, elec_alpha_num ! alpha
call give_integrals_3_body(h,j,i,p,j,i,direct_int) ! <hji|pji>
call give_integrals_3_body(h,j,i,j,p,i,exch_int)
fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
enddo
enddo
!F_a^{aa}(h,p)
do i = 1, elec_beta_num ! alpha
do j = elec_beta_num+1, elec_alpha_num ! alpha
direct_int = three_body_4_index(j,i,h,p)
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
call give_integrals_3_body(h,j,i,i,p,j,exchange_int_231)
call give_integrals_3_body(h,j,i,j,i,p,exchange_int_312)
call give_integrals_3_body(h,j,i,p,i,j,exchange_int_23)
call give_integrals_3_body(h,j,i,i,j,p,exchange_int_12)
call give_integrals_3_body(h,j,i,j,p,i,exchange_int_13)
fock_3_mat_a_op_sh(h,p) -= ( direct_int + exchange_int_231 + exchange_int_312 &
- exchange_int_23 & ! i <-> j
- exchange_int_12 & ! p <-> j
- exchange_int_13 )! p <-> i
enddo
enddo
enddo
enddo
! symmetrized
! do p = 1, elec_beta_num
! do h = elec_alpha_num +1, mo_num
! fock_3_mat_a_op_sh(h,p) = fock_3_mat_a_op_sh(p,h)
! enddo
! enddo
! do h = elec_beta_num+1, elec_alpha_num
! do p = elec_alpha_num +1, mo_num
! !F_a^{bb}(h,p)
! do i = 1, elec_beta_num
! do j = i+1, elec_beta_num
! call give_integrals_3_body(h,j,i,p,j,i,direct_int)
! call give_integrals_3_body(h,j,i,p,i,j,exch_int)
! fock_3_mat_a_op_sh(h,p) -= direct_int - exch_int
! enddo
! enddo
! enddo
! enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_mat_b_op_sh, (mo_num, mo_num)]
implicit none
integer :: h,p,i,j
double precision :: direct_int, exch_int
fock_3_mat_b_op_sh = 0.d0
do h = 1, elec_beta_num
do p = elec_alpha_num +1, mo_num
!F_b^{aa}(h,p)
do i = 1, elec_beta_num
do j = elec_beta_num+1, elec_alpha_num
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
call give_integrals_3_body(h,j,i,p,i,j,exch_int)
fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
enddo
enddo
!F_b^{ab}(h,p)
do i = elec_beta_num+1, elec_beta_num
do j = 1, elec_beta_num
call give_integrals_3_body(h,j,i,p,j,i,direct_int)
call give_integrals_3_body(h,j,i,j,p,i,exch_int)
fock_3_mat_b_op_sh(h,p) += direct_int - exch_int
enddo
enddo
enddo
enddo
END_PROVIDER

160
src/tc_scf/fock_three_bi_ortho.irp.f Normal file
View File

@ -0,0 +1,160 @@
! ---
BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh_bi_orth_old, (mo_num, mo_num)]
BEGIN_DOC
! Fock matrix for opposite spin contribution for bi ortho
END_DOC
implicit none
integer :: j, m, i, a
double precision :: direct_int, exch_int
fock_3_mat_a_op_sh_bi_orth_old = 0.d0
do i = 1, mo_num ! alpha single excitation
do a = 1, mo_num ! alpha single excitation
! ---
do j = 1, elec_beta_num
do m = 1, elec_beta_num
call give_integrals_3_body_bi_ort(a, m, j, i, m, j, direct_int)
fock_3_mat_a_op_sh_bi_orth_old(a,i) += 1.d0 * direct_int
call give_integrals_3_body_bi_ort(a, m, j, j, m, i, exch_int)
fock_3_mat_a_op_sh_bi_orth_old(a,i) += -1.d0 * exch_int
enddo
enddo
! ---
do j = 1, elec_beta_num ! beta
do m = j+1, elec_beta_num ! beta
call give_integrals_3_body_bi_ort(a, m, j, i, m, j, direct_int)
fock_3_mat_a_op_sh_bi_orth_old(a,i) += 1.d0 * direct_int
call give_integrals_3_body_bi_ort(a, m, j, i, j, m, exch_int)
fock_3_mat_a_op_sh_bi_orth_old(a,i) += -1.d0 * exch_int
enddo
enddo
! ---
enddo
enddo
fock_3_mat_a_op_sh_bi_orth_old = - fock_3_mat_a_op_sh_bi_orth_old
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, fock_3_mat_a_op_sh_bi_orth, (mo_num, mo_num)]
BEGIN_DOC
! Fock matrix for opposite spin contribution for bi ortho
END_DOC
implicit none
integer :: i, a
double precision :: integral1, integral2, integral3
fock_3_mat_a_op_sh_bi_orth = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, a, integral1, integral2, integral3) &
!$OMP SHARED (mo_num, fock_3_mat_a_op_sh_bi_orth)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num ! alpha single excitation
do a = 1, mo_num ! alpha single excitation
call direct_term_imj_bi_ortho(a, i, integral1)
call exch_term_jmi_bi_ortho (a, i, integral2)
call exch_term_ijm_bi_ortho (a, i, integral3)
fock_3_mat_a_op_sh_bi_orth(a,i) += 1.5d0 * integral1 - integral2 - 0.5d0 * integral3
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
fock_3_mat_a_op_sh_bi_orth = - fock_3_mat_a_op_sh_bi_orth
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, fock_3_mat_a_sa_sh_bi_orth_old, (mo_num, mo_num)]
BEGIN_DOC
! Fock matrix for same spin contribution for bi ortho
END_DOC
implicit none
integer :: j, m, i, a
double precision :: direct_int, cyclic_1, cyclic_2, non_cyclic_1, non_cyclic_2, non_cyclic_3
fock_3_mat_a_sa_sh_bi_orth_old = 0.d0
do i = 1, mo_num
do a = 1, mo_num
do j = 1, elec_beta_num
do m = j+1, elec_beta_num
call give_integrals_3_body_bi_ort(a, m, j, i, m, j, direct_int)
call give_integrals_3_body_bi_ort(a, m, j, j, i, m, cyclic_1)
call give_integrals_3_body_bi_ort(a, m, j, m, j, i, cyclic_2)
fock_3_mat_a_sa_sh_bi_orth_old(a,i) += direct_int + cyclic_1 + cyclic_2
call give_integrals_3_body_bi_ort(a, m, j, j, m, i, non_cyclic_1)
call give_integrals_3_body_bi_ort(a, m, j, i, j, m, non_cyclic_2)
call give_integrals_3_body_bi_ort(a, m, j, m, i, j, non_cyclic_3)
fock_3_mat_a_sa_sh_bi_orth_old(a,i) += -1.d0 * (non_cyclic_1 + non_cyclic_2 + non_cyclic_3)
enddo
enddo
enddo
enddo
fock_3_mat_a_sa_sh_bi_orth_old = -fock_3_mat_a_sa_sh_bi_orth_old
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, fock_3_mat_a_sa_sh_bi_orth, (mo_num, mo_num)]
BEGIN_DOC
! Fock matrix for same spin contribution for bi ortho
END_DOC
implicit none
integer :: j, m, i, a
double precision :: integral1, integral2, integral3, integral4
fock_3_mat_a_sa_sh_bi_orth = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, a, integral1, integral2, integral3, integral4) &
!$OMP SHARED (mo_num, fock_3_mat_a_sa_sh_bi_orth)
!$OMP DO SCHEDULE (dynamic)
do i = 1, mo_num
do a = 1, mo_num
call direct_term_imj_bi_ortho(a, i, integral1)
call cyclic_term_jim_bi_ortho(a, i, integral2)
call exch_term_jmi_bi_ortho (a, i, integral3)
call exch_term_ijm_bi_ortho (a, i, integral4)
fock_3_mat_a_sa_sh_bi_orth(a,i) += 0.5d0 * (integral1 - integral4) + integral2 - integral3
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
fock_3_mat_a_sa_sh_bi_orth = -fock_3_mat_a_sa_sh_bi_orth
END_PROVIDER
! ---

140
src/tc_scf/fock_three_utils.irp.f Normal file
View File

@ -0,0 +1,140 @@
BEGIN_PROVIDER [ double precision, fock_3_w_kk_sum, (n_points_final_grid,3)]
implicit none
integer :: mm, ipoint,k
double precision :: w_kk
fock_3_w_kk_sum = 0.d0
do k = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
w_kk = x_W_ij_erf_rk(ipoint,mm,k,k)
fock_3_w_kk_sum(ipoint,mm) += w_kk
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_w_ki_mos_k, (n_points_final_grid,3,mo_num)]
implicit none
integer :: mm, ipoint,k,i
double precision :: w_ki, mo_k
fock_3_w_ki_mos_k = 0.d0
do i = 1, mo_num
do k = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
w_ki = x_W_ij_erf_rk(ipoint,mm,k,i)
mo_k = mos_in_r_array(k,ipoint)
fock_3_w_ki_mos_k(ipoint,mm,i) += w_ki * mo_k
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_w_kl_w_kl, (n_points_final_grid,3)]
implicit none
integer :: k,j,ipoint,mm
double precision :: w_kj
fock_3_w_kl_w_kl = 0.d0
do j = 1, elec_beta_num
do k = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
w_kj = x_W_ij_erf_rk(ipoint,mm,k,j)
fock_3_w_kl_w_kl(ipoint,mm) += w_kj * w_kj
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_rho_beta, (n_points_final_grid)]
implicit none
integer :: ipoint,k
fock_3_rho_beta = 0.d0
do ipoint = 1, n_points_final_grid
do k = 1, elec_beta_num
fock_3_rho_beta(ipoint) += mos_in_r_array(k,ipoint) * mos_in_r_array(k,ipoint)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_w_kl_mo_k_mo_l, (n_points_final_grid,3)]
implicit none
integer :: ipoint,k,l,mm
double precision :: mos_k, mos_l, w_kl
fock_3_w_kl_mo_k_mo_l = 0.d0
do k = 1, elec_beta_num
do l = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
mos_k = mos_in_r_array_transp(ipoint,k)
mos_l = mos_in_r_array_transp(ipoint,l)
w_kl = x_W_ij_erf_rk(ipoint,mm,l,k)
fock_3_w_kl_mo_k_mo_l(ipoint,mm) += w_kl * mos_k * mos_l
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_w_ki_wk_a, (n_points_final_grid,3,mo_num, mo_num)]
implicit none
integer :: ipoint,i,a,k,mm
double precision :: w_ki,w_ka
fock_3_w_ki_wk_a = 0.d0
do i = 1, mo_num
do a = 1, mo_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
do k = 1, elec_beta_num
w_ki = x_W_ij_erf_rk(ipoint,mm,k,i)
w_ka = x_W_ij_erf_rk(ipoint,mm,k,a)
fock_3_w_ki_wk_a(ipoint,mm,a,i) += w_ki * w_ka
enddo
enddo
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_trace_w_tilde, (n_points_final_grid,3)]
implicit none
integer :: ipoint,k,mm
fock_3_trace_w_tilde = 0.d0
do k = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
fock_3_trace_w_tilde(ipoint,mm) += fock_3_w_ki_wk_a(ipoint,mm,k,k)
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, fock_3_w_kl_wla_phi_k, (n_points_final_grid,3,mo_num)]
implicit none
integer :: ipoint,a,k,mm,l
double precision :: w_kl,w_la, mo_k
fock_3_w_kl_wla_phi_k = 0.d0
do a = 1, mo_num
do k = 1, elec_beta_num
do l = 1, elec_beta_num
do mm = 1, 3
do ipoint = 1, n_points_final_grid
w_kl = x_W_ij_erf_rk(ipoint,mm,l,k)
w_la = x_W_ij_erf_rk(ipoint,mm,l,a)
mo_k = mos_in_r_array_transp(ipoint,k)
fock_3_w_kl_wla_phi_k(ipoint,mm,a) += w_kl * w_la * mo_k
enddo
enddo
enddo
enddo
enddo
END_PROVIDER

391
src/tc_scf/integrals_in_r_stuff.irp.f Normal file
View File

@ -0,0 +1,391 @@
! ---
BEGIN_PROVIDER [ double precision, tc_scf_dm_in_r, (n_points_final_grid) ]
implicit none
integer :: i, j
tc_scf_dm_in_r = 0.d0
do i = 1, n_points_final_grid
do j = 1, elec_beta_num
tc_scf_dm_in_r(i) += mos_r_in_r_array(j,i) * mos_l_in_r_array(j,i)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, w_sum_in_r, (n_points_final_grid, 3)]
implicit none
integer :: ipoint, j, xi
w_sum_in_r = 0.d0
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
!w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,j)
w_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ww_sum_in_r, (n_points_final_grid, 3)]
implicit none
integer :: ipoint, j, xi
double precision :: tmp
ww_sum_in_r = 0.d0
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
tmp = x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
ww_sum_in_r(ipoint,xi) += tmp * tmp
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_r_in_r, (n_points_final_grid, 3, mo_num)]
implicit none
integer :: i, j, xi, ipoint
! TODO: call lapack
W1_r_in_r = 0.d0
do i = 1, mo_num
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_r_in_r(ipoint,xi,i) += mos_r_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_l_in_r, (n_points_final_grid, 3, mo_num)]
implicit none
integer :: i, j, xi, ipoint
! TODO: call lapack
W1_l_in_r = 0.d0
do i = 1, mo_num
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_l_in_r(ipoint,xi,i) += mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_in_r, (n_points_final_grid, 3)]
implicit none
integer :: j, xi, ipoint
! TODO: call lapack
W1_in_r = 0.d0
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_in_r(ipoint,xi) += W1_l_in_r(ipoint,xi,j) * mos_r_in_r_array_transp(ipoint,j)
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_diag_in_r, (n_points_final_grid, 3)]
implicit none
integer :: j, xi, ipoint
! TODO: call lapack
W1_diag_in_r = 0.d0
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_diag_in_r(ipoint,xi) += mos_r_in_r_array_transp(ipoint,j) * mos_l_in_r_array_transp(ipoint,j) * x_W_ki_bi_ortho_erf_rk_diag(ipoint,xi,j)
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_sum_in_r, (n_points_final_grid, 3)]
implicit none
integer :: i, j, xi, ipoint
! TODO: call lapack
v_sum_in_r = 0.d0
do i = 1, elec_beta_num
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
v_sum_in_r(ipoint,xi) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_W1_r_in_r, (n_points_final_grid, 3, mo_num)]
implicit none
integer :: i, m, xi, ipoint
! TODO: call lapack
W1_W1_r_in_r = 0.d0
do i = 1, mo_num
do m = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_W1_r_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,m,i) * W1_r_in_r(ipoint,xi,m)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, W1_W1_l_in_r, (n_points_final_grid, 3, mo_num)]
implicit none
integer :: i, j, xi, ipoint
! TODO: call lapack
W1_W1_l_in_r = 0.d0
do i = 1, mo_num
do j = 1, elec_beta_num
do xi = 1, 3
do ipoint = 1, n_points_final_grid
W1_W1_l_in_r(ipoint,xi,i) += x_W_ki_bi_ortho_erf_rk(ipoint,xi,i,j) * W1_l_in_r(ipoint,xi,j)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
subroutine direct_term_imj_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j,m = 1, elec_beta_num) < a m j | i m j > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi
double precision :: weight, tmp
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
!integral += ( mos_l_in_r_array(a,ipoint) * mos_r_in_r_array(i,ipoint) * w_sum_in_r(ipoint,xi) * w_sum_in_r(ipoint,xi) &
! + 2.d0 * tc_scf_dm_in_r(ipoint) * w_sum_in_r(ipoint,xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) ) * weight
tmp = w_sum_in_r(ipoint,xi)
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * tmp * tmp &
+ 2.d0 * tc_scf_dm_in_r(ipoint) * tmp * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) &
) * weight
enddo
enddo
end
! ---
subroutine exch_term_jmi_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j,m = 1, elec_beta_num) < a m j | j m i > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi, j
double precision :: weight, tmp
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
tmp = 0.d0
do j = 1, elec_beta_num
tmp = tmp + x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,j) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,j,i)
enddo
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_r_in_r(ipoint,xi,i) * w_sum_in_r(ipoint,xi) &
+ tc_scf_dm_in_r(ipoint) * tmp &
+ mos_r_in_r_array_transp(ipoint,i) * W1_l_in_r(ipoint,xi,a) * w_sum_in_r(ipoint,xi) &
) * weight
enddo
enddo
end
! ---
subroutine exch_term_ijm_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j,m = 1, elec_beta_num) < a m j | i j m > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi
double precision :: weight
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * v_sum_in_r(ipoint,xi) &
+ 2.d0 * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) * W1_in_r(ipoint,xi) &
) * weight
enddo
enddo
end
! ---
subroutine direct_term_ijj_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j = 1, elec_beta_num) < a j j | i j j > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi
double precision :: weight
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += ( mos_l_in_r_array_transp(ipoint,a) * mos_r_in_r_array_transp(ipoint,i) * ww_sum_in_r(ipoint,xi) &
+ 2.d0 * W1_diag_in_r(ipoint, xi) * x_W_ki_bi_ortho_erf_rk(ipoint,xi,a,i) &
) * weight
enddo
enddo
end
! ---
subroutine cyclic_term_jim_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j,m = 1, elec_beta_num) < a m j | j i m > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi
double precision :: weight
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i) &
+ W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i) &
+ W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i) &
) * weight
enddo
enddo
end
! ---
subroutine cyclic_term_mji_bi_ortho(a, i, integral)
BEGIN_DOC
! computes sum_(j,m = 1, elec_beta_num) < a m j | m j i > with bi ortho mos
END_DOC
implicit none
integer, intent(in) :: i, a
double precision, intent(out) :: integral
integer :: ipoint, xi
double precision :: weight
integral = 0.d0
do xi = 1, 3
do ipoint = 1, n_points_final_grid
weight = final_weight_at_r_vector(ipoint)
integral += ( mos_l_in_r_array_transp(ipoint,a) * W1_W1_r_in_r(ipoint,xi,i) &
+ W1_l_in_r(ipoint,xi,a) * W1_r_in_r(ipoint,xi,i) &
+ W1_W1_l_in_r(ipoint,xi,a) * mos_r_in_r_array_transp(ipoint,i) &
) * weight
enddo
enddo
end
! ---

176
src/tc_scf/molden_lr_mos.irp.f Normal file
View File

@ -0,0 +1,176 @@
program molden
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
print *, 'starting ...'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
! my_n_pt_r_grid = 10 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call molden_lr
end
subroutine molden_lr
implicit none
BEGIN_DOC
! Produces a Molden file
END_DOC
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
integer :: i,j,k,l
double precision, parameter :: a0 = 0.529177249d0
PROVIDE ezfio_filename
output=trim(ezfio_filename)//'.mol'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
write(i_unit_output,'(A)') '[Molden Format]'
write(i_unit_output,'(A)') '[Atoms] Angs'
do i = 1, nucl_num
write(i_unit_output,'(A2,2X,I4,2X,I4,3(2X,F15.10))') &
trim(element_name(int(nucl_charge(i)))), &
i, &
int(nucl_charge(i)), &
nucl_coord(i,1)*a0, nucl_coord(i,2)*a0, nucl_coord(i,3)*a0
enddo
write(i_unit_output,'(A)') '[GTO]'
character*(1) :: character_shell
integer :: i_shell,i_prim,i_ao
integer :: iorder(ao_num)
integer :: nsort(ao_num)
i_shell = 0
i_prim = 0
do i=1,nucl_num
write(i_unit_output,*) i, 0
do j=1,nucl_num_shell_aos(i)
i_shell +=1
i_ao = nucl_list_shell_aos(i,j)
character_shell = trim(ao_l_char(i_ao))
write(i_unit_output,*) character_shell, ao_prim_num(i_ao), '1.00'
do k = 1, ao_prim_num(i_ao)
i_prim +=1
write(i_unit_output,'(E20.10,2X,E20.10)') ao_expo(i_ao,k), ao_coef(i_ao,k)
enddo
l = i_ao
do while ( ao_l(l) == ao_l(i_ao) )
nsort(l) = i*10000 + j*100
l += 1
if (l > ao_num) exit
enddo
enddo
write(i_unit_output,*)''
enddo
do i=1,ao_num
iorder(i) = i
! p
if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 3
! d
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
! f
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 10
! g
else if ((ao_power(i,1) == 4 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 1
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 4 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 2
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 4 )) then
nsort(i) += 3
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 4
else if ((ao_power(i,1) == 3 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 5
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 6
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 3 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 7
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 8
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 3 )) then
nsort(i) += 9
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 0 )) then
nsort(i) += 10
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 0 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 11
else if ((ao_power(i,1) == 0 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 12
else if ((ao_power(i,1) == 2 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 13
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 2 ).and.(ao_power(i,3) == 1 )) then
nsort(i) += 14
else if ((ao_power(i,1) == 1 ).and.(ao_power(i,2) == 1 ).and.(ao_power(i,3) == 2 )) then
nsort(i) += 15
endif
enddo
call isort(nsort,iorder,ao_num)
write(i_unit_output,'(A)') '[MO]'
do i=1,mo_num
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', mo_occ(i)
do j=1,ao_num
write(i_unit_output, '(I6,2X,E20.10)') j, mo_r_coef(iorder(j),i)
enddo
write (i_unit_output,*) 'Sym= 1'
write (i_unit_output,*) 'Ene=', Fock_matrix_tc_mo_tot(i,i)
write (i_unit_output,*) 'Spin= Alpha'
write (i_unit_output,*) 'Occup=', mo_occ(i)
do j=1,ao_num
write(i_unit_output, '(I6,2X,E20.10)') j, mo_l_coef(iorder(j),i)
enddo
enddo
close(i_unit_output)
end

248
src/tc_scf/rotate_tcscf_orbitals.irp.f Normal file
View File

@ -0,0 +1,248 @@
! ---
program rotate_tcscf_orbitals
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
bi_ortho = .True.
touch bi_ortho
call maximize_overlap()
end
! ---
subroutine maximize_overlap()
implicit none
integer :: i, m, n
double precision :: accu_d, accu_nd
double precision, allocatable :: C(:,:), R(:,:), L(:,:), W(:,:), e(:)
double precision, allocatable :: S(:,:)
n = ao_num
m = mo_num
allocate(L(n,m), R(n,m), C(n,m), W(n,n), e(m))
L = mo_l_coef
R = mo_r_coef
C = mo_coef
W = ao_overlap
print*, ' fock matrix diag elements'
do i = 1, m
e(i) = Fock_matrix_tc_mo_tot(i,i)
print*, e(i)
enddo
! ---
print *, ' overlap before :'
print *, ' '
allocate(S(m,m))
call LTxSxR(n, m, L, W, R, S)
!print*, " L.T x R"
!do i = 1, m
! write(*, '(100(F16.10,X))') S(i,i)
!enddo
call LTxSxR(n, m, L, W, C, S)
print*, " L.T x C"
do i = 1, m
write(*, '(100(F16.10,X))') S(i,:)
enddo
call LTxSxR(n, m, C, W, R, S)
print*, " C.T x R"
do i = 1, m
write(*, '(100(F16.10,X))') S(i,:)
enddo
deallocate(S)
! ---
call rotate_degen_eigvec(n, m, e, C, W, L, R)
! ---
print *, ' overlap after :'
print *, ' '
allocate(S(m,m))
call LTxSxR(n, m, L, W, R, S)
!print*, " L.T x R"
!do i = 1, m
! write(*, '(100(F16.10,X))') S(i,i)
!enddo
call LTxSxR(n, m, L, W, C, S)
print*, " L.T x C"
do i = 1, m
write(*, '(100(F16.10,X))') S(i,:)
enddo
call LTxSxR(n, m, C, W, R, S)
print*, " C.T x R"
do i = 1, m
write(*, '(100(F16.10,X))') S(i,:)
enddo
deallocate(S)
! ---
mo_l_coef = L
mo_r_coef = R
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
! ---
deallocate(L, R, C, W, e)
end subroutine maximize_overlap
! ---
subroutine rotate_degen_eigvec(n, m, e0, C0, W0, L0, R0)
implicit none
integer, intent(in) :: n, m
double precision, intent(in) :: e0(m), W0(n,n), C0(n,m)
double precision, intent(inout) :: L0(n,m), R0(n,m)
integer :: i, j, k, kk, mm, id1
double precision :: ei, ej, de, de_thr
integer, allocatable :: deg_num(:)
double precision, allocatable :: L(:,:), R(:,:), C(:,:), Lnew(:,:), Rnew(:,:), tmp(:,:)
!double precision, allocatable :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
double precision, allocatable :: S(:,:), Snew(:,:), T(:,:), Ttmp(:,:), Stmp(:,:)
!real*8 :: S(m,m), Snew(m,m), T(m,m)
id1 = 700
allocate(S(id1,id1), Snew(id1,id1), T(id1,id1))
! ---
allocate( deg_num(m) )
do i = 1, m
deg_num(i) = 1
enddo
de_thr = 1d-10
do i = 1, m-1
ei = e0(i)
! already considered in degen vectors
if(deg_num(i).eq.0) cycle
do j = i+1, m
ej = e0(j)
de = dabs(ei - ej)
if(de .lt. de_thr) then
deg_num(i) = deg_num(i) + 1
deg_num(j) = 0
endif
enddo
enddo
do i = 1, m
if(deg_num(i).gt.1) then
print *, ' degen on', i, deg_num(i)
endif
enddo
! ---
do i = 1, m
mm = deg_num(i)
if(mm .gt. 1) then
allocate(L(n,mm), R(n,mm), C(n,mm))
do j = 1, mm
L(1:n,j) = L0(1:n,i+j-1)
R(1:n,j) = R0(1:n,i+j-1)
C(1:n,j) = C0(1:n,i+j-1)
enddo
! ---
! C.T x W0 x R
allocate(tmp(mm,n), Stmp(mm,mm))
call dgemm( 'T', 'N', mm, n, n, 1.d0 &
, C, size(C, 1), W0, size(W0, 1) &
, 0.d0, tmp, size(tmp, 1) )
call dgemm( 'N', 'N', mm, mm, n, 1.d0 &
, tmp, size(tmp, 1), R, size(R, 1) &
, 0.d0, Stmp, size(Stmp, 1) )
deallocate(C, tmp)
S = 0.d0
do k = 1, mm
do kk = 1, mm
S(kk,k) = Stmp(kk,k)
enddo
enddo
deallocate(Stmp)
!print*, " overlap bef"
!do k = 1, mm
! write(*, '(100(F16.10,X))') (S(k,kk), kk=1, mm)
!enddo
T = 0.d0
Snew = 0.d0
call maxovl(mm, mm, S, T, Snew)
!print*, " overlap aft"
!do k = 1, mm
! write(*, '(100(F16.10,X))') (Snew(k,kk), kk=1, mm)
!enddo
allocate(Ttmp(mm,mm))
Ttmp(1:mm,1:mm) = T(1:mm,1:mm)
allocate(Lnew(n,mm), Rnew(n,mm))
call dgemm( 'N', 'N', n, mm, mm, 1.d0 &
, R, size(R, 1), Ttmp(1,1), size(Ttmp, 1) &
, 0.d0, Rnew, size(Rnew, 1) )
call dgemm( 'N', 'N', n, mm, mm, 1.d0 &
, L, size(L, 1), Ttmp(1,1), size(Ttmp, 1) &
, 0.d0, Lnew, size(Lnew, 1) )
deallocate(L, R)
deallocate(Ttmp)
! ---
do j = 1, mm
L0(1:n,i+j-1) = Lnew(1:n,j)
R0(1:n,i+j-1) = Rnew(1:n,j)
enddo
deallocate(Lnew, Rnew)
endif
enddo
deallocate(S, Snew, T)
end subroutine rotate_degen_eigvec
! ---

179
src/tc_scf/tc_scf.irp.f Normal file
View File

@ -0,0 +1,179 @@
program tc_scf
BEGIN_DOC
! TODO : Put the documentation of the program here
END_DOC
implicit none
print *, 'starting ...'
my_grid_becke = .True.
my_n_pt_r_grid = 30
my_n_pt_a_grid = 50
! my_n_pt_r_grid = 10 ! small grid for quick debug
! my_n_pt_a_grid = 26 ! small grid for quick debug
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
!call create_guess
!call orthonormalize_mos
call routine_scf()
end
! ---
subroutine create_guess
BEGIN_DOC
! Create a MO guess if no MOs are present in the EZFIO directory
END_DOC
implicit none
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_mo_basis_mo_coef(exists)
if (.not.exists) then
mo_label = 'Guess'
if (mo_guess_type == "HCore") then
mo_coef = ao_ortho_lowdin_coef
call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10)
TOUCH mo_coef
call mo_as_eigvectors_of_mo_matrix(mo_one_e_integrals, &
size(mo_one_e_integrals,1), &
size(mo_one_e_integrals,2), &
mo_label,1,.false.)
call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef,1), 1.d-10)
SOFT_TOUCH mo_coef
else if (mo_guess_type == "Huckel") then
call huckel_guess
else
print *, 'Unrecognized MO guess type : '//mo_guess_type
stop 1
endif
SOFT_TOUCH mo_label
endif
end subroutine create_guess
! ---
subroutine routine_scf()
implicit none
integer :: i, j, it
double precision :: e_save, e_delta, rho_delta
double precision, allocatable :: rho_old(:,:), rho_new(:,:)
allocate(rho_old(ao_num,ao_num), rho_new(ao_num,ao_num))
it = 0
print*,'iteration = ', it
!print*,'grad_hermit = ', grad_hermit
print*,'***'
print*,'TC HF total energy = ', TC_HF_energy
print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy
print*,'TC HF 2 e energy = ', TC_HF_two_e_energy
if(.not. bi_ortho)then
print*,'TC HF 3 body = ', diag_three_elem_hf
endif
print*,'***'
e_delta = 10.d0
e_save = 0.d0 !TC_HF_energy
rho_delta = 10.d0
if(bi_ortho)then
mo_l_coef = fock_tc_leigvec_ao
mo_r_coef = fock_tc_reigvec_ao
rho_old = TCSCF_bi_ort_dm_ao
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
TOUCH mo_l_coef mo_r_coef
else
print*,'grad_hermit = ',grad_hermit
call save_good_hermit_tc_eigvectors
TOUCH mo_coef
call save_mos
endif
! ---
if(bi_ortho) then
!do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. dsqrt(thresh_tcscf)) )
!do while( it .lt. n_it_tcscf_max .and. (e_delta .gt. thresh_tcscf) )
do while( it .lt. n_it_tcscf_max .and. (rho_delta .gt. thresh_tcscf) )
it += 1
print*,'iteration = ', it
print*,'***'
print*,'TC HF total energy = ', TC_HF_energy
print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy
print*,'TC HF 2 non hermit = ', TC_HF_two_e_energy
print*,'***'
e_delta = dabs( TC_HF_energy - e_save )
print*, 'it, delta E = ', it, e_delta
e_save = TC_HF_energy
mo_l_coef = fock_tc_leigvec_ao
mo_r_coef = fock_tc_reigvec_ao
rho_new = TCSCF_bi_ort_dm_ao
!print*, rho_new
rho_delta = 0.d0
do i = 1, ao_num
do j = 1, ao_num
rho_delta += dabs(rho_new(j,i) - rho_old(j,i))
enddo
enddo
print*, ' rho_delta =', rho_delta
rho_old = rho_new
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
TOUCH mo_l_coef mo_r_coef
call ezfio_set_tc_scf_bitc_energy(TC_HF_energy)
enddo
else
do while( (grad_hermit.gt.dsqrt(thresh_tcscf)) .and. it .lt. n_it_tcscf_max )
print*,'grad_hermit = ',grad_hermit
it += 1
print*,'iteration = ', it
print*,'***'
print*,'TC HF total energy = ', TC_HF_energy
print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy
print*,'TC HF 2 e energy = ', TC_HF_two_e_energy
print*,'TC HF 3 body = ', diag_three_elem_hf
print*,'***'
call save_good_hermit_tc_eigvectors
TOUCH mo_coef
call save_mos
enddo
endif
print*,'Energy converged !'
print*,'Diagonal Fock elements '
do i = 1, mo_num
print*,i,Fock_matrix_tc_mo_tot(i,i)
enddo
deallocate(rho_old, rho_new)
end subroutine routine_scf
! ---

25
src/tc_scf/tc_scf_dm.irp.f Normal file
View File

@ -0,0 +1,25 @@
BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) ]
implicit none
if(bi_ortho)then
TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta
else
TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num) ]
implicit none
if(bi_ortho)then
TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha
else
TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_tot, (ao_num, ao_num) ]
implicit none
TCSCF_density_matrix_ao_tot = TCSCF_density_matrix_ao_beta + TCSCF_density_matrix_ao_alpha
END_PROVIDER

32
src/tc_scf/tc_scf_energy.irp.f Normal file
View File

@ -0,0 +1,32 @@
BEGIN_PROVIDER [ double precision, TC_HF_energy]
&BEGIN_PROVIDER [ double precision, TC_HF_one_electron_energy]
&BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy]
BEGIN_DOC
! TC Hartree-Fock energy containing the nuclear repulsion, and its one- and two-body components.
END_DOC
implicit none
integer :: i, j
TC_HF_energy = nuclear_repulsion
TC_HF_one_electron_energy = 0.d0
TC_HF_two_e_energy = 0.d0
do j = 1, ao_num
do i = 1, ao_num
TC_HF_two_e_energy += 0.5d0 * ( two_e_tc_non_hermit_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
+ two_e_tc_non_hermit_integral_beta(i,j) * TCSCF_density_matrix_ao_beta(i,j) )
TC_HF_one_electron_energy += ao_one_e_integrals_tc_tot(i,j) &
* (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) )
enddo
enddo
TC_HF_energy += TC_HF_one_electron_energy + TC_HF_two_e_energy
TC_HF_energy += diag_three_elem_hf
END_PROVIDER
! ---

42
src/tc_scf/tc_scf_utils.irp.f Normal file
View File

@ -0,0 +1,42 @@
! ---
subroutine LTxSxR(n, m, L, S, R, C)
implicit none
integer, intent(in) :: n, m
double precision, intent(in) :: L(n,m), S(n,n), R(n,m)
double precision, intent(out) :: C(m,m)
integer :: i, j
double precision :: accu_d, accu_nd
double precision, allocatable :: tmp(:,:)
! L.T x S x R
allocate(tmp(m,n))
call dgemm( 'T', 'N', m, n, n, 1.d0 &
, L, size(L, 1), S, size(S, 1) &
, 0.d0, tmp, size(tmp, 1) )
call dgemm( 'N', 'N', m, m, n, 1.d0 &
, tmp, size(tmp, 1), R, size(R, 1) &
, 0.d0, C, size(C, 1) )
deallocate(tmp)
accu_d = 0.d0
accu_nd = 0.d0
do i = 1, m
do j = 1, m
if(j.eq.i) then
accu_d += dabs(C(j,i))
else
accu_nd += C(j,i) * C(j,i)
endif
enddo
enddo
accu_nd = dsqrt(accu_nd)
print*, ' accu_d = ', accu_d
print*, ' accu_nd = ', accu_nd
end subroutine LTxR
! ---

304
src/utils/loc.f Normal file
View File

@ -0,0 +1,304 @@
c************************************************************************
subroutine maxovl(n,m,s,t,w)
C
C This subprogram contains an iterative procedure to find the
C unitary transformation of a set of n vectors which maximizes
C the sum of their square overlaps with a set of m reference
C vectors (m.le.n)
C
C S: overlap matrix <ref|vec>
C T: rotation matrix
C W: new overlap matrix
C
C
implicit real*8(a-h,o-y),logical*1(z)
parameter (id1=700)
dimension s(id1,id1),t(id1,id1),w(id1,id1)
data small/1.d-6/
zprt=.true.
niter=1000000
conv=1.d-12
C niter=1000000
C conv=1.d-6
write (6,5) n,m,conv
5 format (//5x,'Unitary transformation of',i3,' vectors'/
* 5x,'following the principle of maximum overlap with a set of',
* i3,' reference vectors'/5x,'required convergence on rotation ',
* 'angle =',f13.10///5x,'Starting overlap matrix'/)
do 6 i=1,m
write (6,145) i
6 write (6,150) (s(i,j),j=1,n)
8 mm=m-1
if (m.lt.n) mm=m
iter=0
do 20 j=1,n
do 16 i=1,n
t(i,j)=0.d0
16 continue
do 18 i=1,m
18 w(i,j)=s(i,j)
20 t(j,j)=1.d0
sum=0.d0
do 10 i=1,m
sum=sum+s(i,i)*s(i,i)
10 continue
sum=sum/m
if (zprt) write (6,12) sum
12 format (//5x,'Average square overlap =',f10.6)
if (n.eq.1) goto 100
last=n
j=1
21 if (j.ge.last) goto 30
sum=0.d0
do 22 i=1,n
22 sum=sum+s(i,j)*s(i,j)
if (sum.gt.small) goto 28
do 24 i=1,n
sij=s(i,j)
s(i,j)=-s(i,last)
s(i,last)=sij
tij=t(i,j)
t(i,j)=-t(i,last)
t(i,last)=tij
24 continue
last=last-1
goto 21
28 j=j+1
goto 21
30 iter=iter+1
imax=0
jmax=0
dmax=0.d0
amax=0.d0
do 60 i=1,mm
ip=i+1
do 50 j=ip,n
a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
b=-s(i,i)*s(i,j)
if (j.gt.m) goto 31
a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
b=b+s(j,i)*s(j,j)
31 b=b+b
if (a.eq.0.d0) goto 32
ba=b/a
if (dabs(ba).gt.small) goto 32
if (a.gt.0.d0) goto 33
tang=-0.5d0*ba
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
32 tang=0.d0
if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
33 cosine=0.d0
sine=1.d0
34 delta=sine*(a*sine+b*cosine)
if (zprt.and.delta.lt.0.d0) write (6,71) i,j,a,b,sine,cosine,delta
do 35 k=1,m
p=s(k,i)*cosine-s(k,j)*sine
q=s(k,i)*sine+s(k,j)*cosine
s(k,i)=p
35 s(k,j)=q
do 40 k=1,n
p=t(k,i)*cosine-t(k,j)*sine
q=t(k,i)*sine+t(k,j)*cosine
t(k,i)=p
t(k,j)=q
40 continue
45 d=dabs(sine)
if (d.le.amax) goto 50
imax=i
jmax=j
amax=d
dmax=delta
50 continue
60 continue
if (zprt) write (6,70) iter,amax,imax,jmax,dmax
70 format (' iter=',i4,' largest rotation=',f12.8,
* ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
if (amax.lt.conv) goto 100
if (iter.lt.niter) goto 30
write (6,80)
write (6,*) 'niter=',niter
80 format (//5x,'*** maximum number of cycles exceeded ',
* 'in subroutine maxovl ***'//)
stop
100 continue
do 120 j=1,n
if (s(j,j).gt.0.d0) goto 120
do 105 i=1,m
105 s(i,j)=-s(i,j)
do 110 i=1,n
110 t(i,j)=-t(i,j)
120 continue
sum=0.d0
do 125 i=1,m
125 sum=sum+s(i,i)*s(i,i)
sum=sum/m
do 122 i=1,m
do 122 j=1,n
sw=s(i,j)
s(i,j)=w(i,j)
122 w(i,j)=sw
if (.not.zprt) return
write (6,12) sum
write (6,130)
130 format (//5x,'transformation matrix')
do 140 i=1,n
write (6,145) i
140 write (6,150) (t(i,j),j=1,n)
145 format (i8)
150 format (2x,10f12.8)
write (6,160)
160 format (//5x,'new overlap matrix'/)
do 170 i=1,m
write (6,145) i
170 write (6,150) (w(i,j),j=1,n)
return
end
c************************************************************************
subroutine maxovl_no_print(n,m,s,t,w)
C
C This subprogram contains an iterative procedure to find the
C unitary transformation of a set of n vectors which maximizes
C the sum of their square overlaps with a set of m reference
C vectors (m.le.n)
C
C S: overlap matrix <ref|vec>
C T: rotation matrix
C W: new overlap matrix
C
C
implicit real*8(a-h,o-y),logical*1(z)
parameter (id1=300)
dimension s(id1,id1),t(id1,id1),w(id1,id1)
data small/1.d-6/
zprt=.false.
niter=1000000
conv=1.d-8
C niter=1000000
C conv=1.d-6
8 mm=m-1
if (m.lt.n) mm=m
iter=0
do 20 j=1,n
do 16 i=1,n
t(i,j)=0.d0
16 continue
do 18 i=1,m
18 w(i,j)=s(i,j)
20 t(j,j)=1.d0
sum=0.d0
do 10 i=1,m
sum=sum+s(i,i)*s(i,i)
10 continue
sum=sum/m
12 format (//5x,'Average square overlap =',f10.6)
if (n.eq.1) goto 100
last=n
j=1
21 if (j.ge.last) goto 30
sum=0.d0
do 22 i=1,n
22 sum=sum+s(i,j)*s(i,j)
if (sum.gt.small) goto 28
do 24 i=1,n
sij=s(i,j)
s(i,j)=-s(i,last)
s(i,last)=sij
tij=t(i,j)
t(i,j)=-t(i,last)
t(i,last)=tij
24 continue
last=last-1
goto 21
28 j=j+1
goto 21
30 iter=iter+1
imax=0
jmax=0
dmax=0.d0
amax=0.d0
do 60 i=1,mm
ip=i+1
do 50 j=ip,n
a=s(i,j)*s(i,j)-s(i,i)*s(i,i)
b=-s(i,i)*s(i,j)
if (j.gt.m) goto 31
a=a+s(j,i)*s(j,i)-s(j,j)*s(j,j)
b=b+s(j,i)*s(j,j)
31 b=b+b
if (a.eq.0.d0) goto 32
ba=b/a
if (dabs(ba).gt.small) goto 32
if (a.gt.0.d0) goto 33
tang=-0.5d0*ba
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
32 tang=0.d0
if (b.ne.0.d0) tang=(a+dsqrt(a*a+b*b))/b
cosine=1.d0/dsqrt(1.d0+tang*tang)
sine=tang*cosine
goto 34
33 cosine=0.d0
sine=1.d0
34 delta=sine*(a*sine+b*cosine)
do 35 k=1,m
p=s(k,i)*cosine-s(k,j)*sine
q=s(k,i)*sine+s(k,j)*cosine
s(k,i)=p
35 s(k,j)=q
do 40 k=1,n
p=t(k,i)*cosine-t(k,j)*sine
q=t(k,i)*sine+t(k,j)*cosine
t(k,i)=p
t(k,j)=q
40 continue
45 d=dabs(sine)
if (d.le.amax) goto 50
imax=i
jmax=j
amax=d
dmax=delta
50 continue
60 continue
70 format (' iter=',i4,' largest rotation=',f12.8,
* ', vectors',i3,' and',i3,', incr. of diag. squares=',g12.5)
71 format (' i,j,a,b,sin,cos,delta =',2i3,5f10.5)
if (amax.lt.conv) goto 100
if (iter.lt.niter) goto 30
80 format (//5x,'*** maximum number of cycles exceeded ',
* 'in subroutine maxovl ***'//)
stop
100 continue
do 120 j=1,n
if (s(j,j).gt.0.d0) goto 120
do 105 i=1,m
105 s(i,j)=-s(i,j)
do 110 i=1,n
110 t(i,j)=-t(i,j)
120 continue
sum=0.d0
do 125 i=1,m
125 sum=sum+s(i,i)*s(i,i)
sum=sum/m
do 122 i=1,m
do 122 j=1,n
sw=s(i,j)
s(i,j)=w(i,j)
122 w(i,j)=sw
return
end