diff --git a/src/tc_scf/EZFIO.cfg b/src/tc_scf/EZFIO.cfg new file mode 100644 index 00000000..313d6f2b --- /dev/null +++ b/src/tc_scf/EZFIO.cfg @@ -0,0 +1,4 @@ +[bitc_energy] +type: Threshold +doc: Energy bi-tc HF +interface: ezfio diff --git a/src/tc_scf/NEED b/src/tc_scf/NEED new file mode 100644 index 00000000..4e340cfe --- /dev/null +++ b/src/tc_scf/NEED @@ -0,0 +1,6 @@ +hartree_fock +bi_ortho_mos +three_body_ints +bi_ort_ints +tc_keywords +non_hermit_dav diff --git a/src/tc_scf/combine_lr_tcscf.irp.f b/src/tc_scf/combine_lr_tcscf.irp.f new file mode 100644 index 00000000..b257f4a5 --- /dev/null +++ b/src/tc_scf/combine_lr_tcscf.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f new file mode 100644 index 00000000..69083e33 --- /dev/null +++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f @@ -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 + diff --git a/src/tc_scf/fock_for_right.irp.f b/src/tc_scf/fock_for_right.irp.f new file mode 100644 index 00000000..5a51b324 --- /dev/null +++ b/src/tc_scf/fock_for_right.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/fock_tc.irp.f b/src/tc_scf/fock_tc.irp.f new file mode 100644 index 00000000..42b429d5 --- /dev/null +++ b/src/tc_scf/fock_tc.irp.f @@ -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) = +! +! 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 + +! --- + diff --git a/src/tc_scf/fock_three.irp.f b/src/tc_scf/fock_three.irp.f new file mode 100644 index 00000000..e4348892 --- /dev/null +++ b/src/tc_scf/fock_three.irp.f @@ -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) ! + 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 diff --git a/src/tc_scf/fock_three_bi_ortho.irp.f b/src/tc_scf/fock_three_bi_ortho.irp.f new file mode 100644 index 00000000..6960ebc2 --- /dev/null +++ b/src/tc_scf/fock_three_bi_ortho.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/fock_three_utils.irp.f b/src/tc_scf/fock_three_utils.irp.f new file mode 100644 index 00000000..5aec1d9e --- /dev/null +++ b/src/tc_scf/fock_three_utils.irp.f @@ -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 + diff --git a/src/tc_scf/integrals_in_r_stuff.irp.f b/src/tc_scf/integrals_in_r_stuff.irp.f new file mode 100644 index 00000000..3ce85a97 --- /dev/null +++ b/src/tc_scf/integrals_in_r_stuff.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/molden_lr_mos.irp.f b/src/tc_scf/molden_lr_mos.irp.f new file mode 100644 index 00000000..735349ba --- /dev/null +++ b/src/tc_scf/molden_lr_mos.irp.f @@ -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 + diff --git a/src/tc_scf/rotate_tcscf_orbitals.irp.f b/src/tc_scf/rotate_tcscf_orbitals.irp.f new file mode 100644 index 00000000..49f8bfa6 --- /dev/null +++ b/src/tc_scf/rotate_tcscf_orbitals.irp.f @@ -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 + +! --- diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f new file mode 100644 index 00000000..5af5206e --- /dev/null +++ b/src/tc_scf/tc_scf.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/tc_scf_dm.irp.f b/src/tc_scf/tc_scf_dm.irp.f new file mode 100644 index 00000000..f6ae3e1f --- /dev/null +++ b/src/tc_scf/tc_scf_dm.irp.f @@ -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 + + diff --git a/src/tc_scf/tc_scf_energy.irp.f b/src/tc_scf/tc_scf_energy.irp.f new file mode 100644 index 00000000..aa2a16ff --- /dev/null +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -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 + +! --- + diff --git a/src/tc_scf/tc_scf_utils.irp.f b/src/tc_scf/tc_scf_utils.irp.f new file mode 100644 index 00000000..09a4a1b9 --- /dev/null +++ b/src/tc_scf/tc_scf_utils.irp.f @@ -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 + +! --- diff --git a/src/utils/loc.f b/src/utils/loc.f new file mode 100644 index 00000000..32d2cd03 --- /dev/null +++ b/src/utils/loc.f @@ -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 +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 +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 +