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..726169d9 --- /dev/null +++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f @@ -0,0 +1,229 @@ +! --- + + 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, j, k, l + double precision :: accu_d, accu_nd, accu_tmp + double precision :: norm + double precision, allocatable :: eigval_right_tmp(:) + double precision, allocatable :: F_tmp(:,:) + + allocate( eigval_right_tmp(mo_num), F_tmp(mo_num,mo_num) ) + + PROVIDE Fock_matrix_tc_mo_tot + + do i = 1, mo_num + do j = 1, mo_num + F_tmp(j,i) = Fock_matrix_tc_mo_tot(j,i) + enddo + enddo + ! insert level shift here + do i = elec_beta_num+1, elec_alpha_num + F_tmp(i,i) += 0.5d0 * level_shift_tcscf + enddo + do i = elec_alpha_num+1, mo_num + F_tmp(i,i) += level_shift_tcscf + enddo + + call non_hrmt_bieig( mo_num, F_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag & + , 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, F_tmp, thresh_biorthog_diag, thresh_biorthog_nondiag & + ! , 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, F_tmp & + ! , fock_tc_leigvec_mo, fock_tc_reigvec_mo & + ! , n_real_tc, eigval_right_tmp ) + !endif + + deallocate(F_tmp) + + +! 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, elec_alpha_num +! print*, i, eigval_fock_tc_mo(i) +! enddo +! do i = elec_alpha_num+1, mo_num +! print*, i, eigval_fock_tc_mo(i) - level_shift_tcscf +! 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. thresh_biorthog_nondiag)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. thresh_biorthog_nondiag) then + print *, ' bi-orthog failed' + print *, ' accu_nd MO = ', accu_nd, thresh_biorthog_nondiag + 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))/dble(mo_num) .gt. thresh_biorthog_diag) then + + print *, ' mo_num = ', mo_num + print *, ' accu_d MO = ', accu_d, thresh_biorthog_diag + print *, ' normalizing vectors ...' + do i = 1, mo_num + norm = dsqrt(dabs(overlap_fock_tc_eigvec_mo(i,i))) + if(norm .gt. thresh_biorthog_diag) 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) ) + + 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. thresh_biorthog_nondiag)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. thresh_biorthog_diag) then + print *, ' bi-orthog failed' + print *, ' accu_nd MO = ', accu_nd, thresh_biorthog_nondiag + 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 + + 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(:,:) + + PROVIDE mo_l_coef mo_r_coef + +! ! 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/diis_tcscf.irp.f b/src/tc_scf/diis_tcscf.irp.f new file mode 100644 index 00000000..ff1077f5 --- /dev/null +++ b/src/tc_scf/diis_tcscf.irp.f @@ -0,0 +1,186 @@ +! --- + +BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero_TCSCF ] + + implicit none + + if(threshold_DIIS_TCSCF == 0.d0) then + threshold_DIIS_nonzero_TCSCF = dsqrt(thresh_tcscf) + else + threshold_DIIS_nonzero_TCSCF = threshold_DIIS_TCSCF + endif + ASSERT(threshold_DIIS_nonzero_TCSCF >= 0.d0) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, Q_alpha, (ao_num, ao_num) ] + + BEGIN_DOC + ! + ! Q_alpha = mo_r_coef x eta_occ_alpha x mo_l_coef.T + ! + ! [Q_alpha]_ij = \sum_{k=1}^{elec_alpha_num} [mo_r_coef]_ik [mo_l_coef]_jk + ! + END_DOC + + implicit none + + Q_alpha = 0.d0 + call dgemm( 'N', 'T', ao_num, ao_num, elec_alpha_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & + , 0.d0, Q_alpha, size(Q_alpha, 1) ) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Q_beta, (ao_num, ao_num) ] + + BEGIN_DOC + ! + ! Q_beta = mo_r_coef x eta_occ_beta x mo_l_coef.T + ! + ! [Q_beta]_ij = \sum_{k=1}^{elec_beta_num} [mo_r_coef]_ik [mo_l_coef]_jk + ! + END_DOC + + implicit none + + Q_beta = 0.d0 + call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1), mo_l_coef, size(mo_l_coef, 1) & + , 0.d0, Q_beta, size(Q_beta, 1) ) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Q_matrix, (ao_num, ao_num) ] + + BEGIN_DOC + ! + ! Q_matrix = 2 mo_r_coef x eta_occ x mo_l_coef.T + ! + ! with: + ! | 1 if i = j = 1, ..., nb of occ orbitals + ! [eta_occ]_ij = | + ! | 0 otherwise + ! + ! the diis error is defines as: + ! e = F_ao x Q x ao_overlap - ao_overlap x Q x F_ao + ! with: + ! mo_l_coef.T x ao_overlap x mo_r_coef = I + ! F_mo = mo_l_coef.T x F_ao x mo_r_coef + ! F_ao = (ao_overlap x mo_r_coef) x F_mo x (ao_overlap x mo_l_coef).T + ! + ! ==> e = 2 ao_overlap x mo_r_coef x [ F_mo x eta_occ - eta_occ x F_mo ] x (ao_overlap x mo_l_coef).T + ! + ! at convergence: + ! F_mo x eta_occ - eta_occ x F_mo = 0 + ! ==> [F_mo]_ij ([eta_occ]_ii - [eta_occ]_jj) = 0 + ! ==> [F_mo]_ia = [F_mo]_ai = 0 where: i = occ and a = vir + ! ==> Brillouin conditions + ! + END_DOC + + implicit none + + if(elec_alpha_num == elec_beta_num) then + Q_matrix = Q_alpha + Q_alpha + else + Q_matrix = Q_alpha + Q_beta + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)] + + implicit none + double precision, allocatable :: tmp(:,:) + + allocate(tmp(ao_num,ao_num)) + + ! F x Q + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & + , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), Q_matrix, size(Q_matrix, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! F x Q x S + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & + , tmp, size(tmp, 1), ao_overlap, size(ao_overlap, 1) & + , 0.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) ) + + ! S x Q + tmp = 0.d0 + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & + , ao_overlap, size(ao_overlap, 1), Q_matrix, size(Q_matrix, 1) & + , 0.d0, tmp, size(tmp, 1) ) + + ! F x Q x S - S x Q x F + call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 & + , tmp, size(tmp, 1), Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & + , 1.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) ) + + deallocate(tmp) + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, FQS_SQF_mo, (mo_num, mo_num)] + + implicit none + + call ao_to_mo_bi_ortho( FQS_SQF_ao, size(FQS_SQF_ao, 1) & + , FQS_SQF_mo, size(FQS_SQF_mo, 1) ) + +END_PROVIDER + +! --- + +! BEGIN_PROVIDER [ double precision, eigenval_Fock_tc_ao, (ao_num) ] +!&BEGIN_PROVIDER [ double precision, eigenvec_Fock_tc_ao, (ao_num,ao_num) ] +! +! BEGIN_DOC +! ! +! ! Eigenvalues and eigenvectors of the Fock matrix over the ao basis +! ! +! ! F' = X.T x F x X where X = ao_overlap^(-1/2) +! ! +! ! F' x Cr' = Cr' x E ==> F Cr = Cr x E with Cr = X x Cr' +! ! F'.T x Cl' = Cl' x E ==> F.T Cl = Cl x E with Cl = X x Cl' +! ! +! END_DOC +! +! implicit none +! double precision, allocatable :: tmp1(:,:), tmp2(:,:) +! +! ! --- +! ! Fock matrix in orthogonal basis: F' = X.T x F x X +! +! allocate(tmp1(ao_num,ao_num)) +! call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 & +! , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), S_half_inv, size(S_half_inv, 1) & +! , 0.d0, tmp1, size(tmp1, 1) ) +! +! allocate(tmp2(ao_num,ao_num)) +! call dgemm( 'T', 'N', ao_num, ao_num, ao_num, 1.d0 & +! , S_half_inv, size(S_half_inv, 1), tmp1, size(tmp1, 1) & +! , 0.d0, tmp2, size(tmp2, 1) ) +! +! ! --- +! +! ! Diagonalize F' to obtain eigenvectors in orthogonal basis C' and eigenvalues +! ! TODO +! +! ! Back-transform eigenvectors: C =X.C' +! +!END_PROVIDER + +! --- + +~ diff --git a/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f new file mode 100644 index 00000000..fccfd837 --- /dev/null +++ b/src/tc_scf/fock_3e_bi_ortho_uhf.irp.f @@ -0,0 +1,405 @@ + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_cs ...' + call wall_time(ti) + + fock_3e_uhf_mo_cs = 0.d0 + + do a = 1, mo_num + do b = 1, mo_num + + do j = 1, elec_beta_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_cs(b,a) -= 0.5d0 * ( 4.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - 2.d0 * I_bij_aji & + - 2.d0 * I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j, o + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_a ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_a = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - 2.d0 * I_bij_jia ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - 2.d0 * I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + ! --- + + do j = o, elec_alpha_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_a(b,a) -= 0.5d0 * ( I_bij_aij & + + I_bij_ija & + + I_bij_jai & + - I_bij_aji & + - I_bij_iaj & + - I_bij_jia ) + + enddo + enddo + + ! --- + + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_b, (mo_num, mo_num)] + + implicit none + integer :: a, b, i, j, o + double precision :: I_bij_aij, I_bij_ija, I_bij_jai, I_bij_aji, I_bij_iaj, I_bij_jia + double precision :: ti, tf + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' PROVIDING fock_3e_uhf_mo_b ...' + call wall_time(ti) + + o = elec_beta_num + 1 + + fock_3e_uhf_mo_b = fock_3e_uhf_mo_cs + + do a = 1, mo_num + do b = 1, mo_num + + ! --- + + do j = o, elec_alpha_num + do i = 1, elec_beta_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_iaj ) + + enddo + enddo + + ! --- + + do j = 1, elec_beta_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( 2.d0 * I_bij_aij & + - I_bij_aji & + - I_bij_jia ) + + enddo + enddo + + ! --- + + do j = o, elec_alpha_num + do i = o, elec_alpha_num + + call give_integrals_3_body_bi_ort(b, i, j, a, i, j, I_bij_aij) + call give_integrals_3_body_bi_ort(b, i, j, i, j, a, I_bij_ija) + call give_integrals_3_body_bi_ort(b, i, j, j, a, i, I_bij_jai) + call give_integrals_3_body_bi_ort(b, i, j, a, j, i, I_bij_aji) + call give_integrals_3_body_bi_ort(b, i, j, i, a, j, I_bij_iaj) + call give_integrals_3_body_bi_ort(b, i, j, j, i, a, I_bij_jia) + + fock_3e_uhf_mo_b(b,a) -= 0.5d0 * ( I_bij_aij & + - I_bij_aji ) + + enddo + enddo + + ! --- + + enddo + enddo + + call wall_time(tf) + !print *, ' total Wall time for fock_3e_uhf_mo_b =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_a, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + double precision, allocatable :: f_tmp(:,:) + + print *, ' PROVIDING fock_3e_uhf_ao_a ...' + call wall_time(ti) + + fock_3e_uhf_ao_a = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_a) + + allocate(f_tmp(ao_num,ao_num)) + f_tmp = 0.d0 + + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + do mu = 1, ao_num + do nu = 1, ao_num + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_a * dm_dk_a * i_mugd_eknu & + + dm_ge_a * dm_dk_a * i_mugd_knue & + - dm_ge_a * dm_dk * i_mugd_enuk & + - dm_ge * dm_dk_a * i_mugd_kenu & + - dm_ge_a * dm_dk_a * i_mugd_nuke & + - dm_ge_b * dm_dk_b * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do mu = 1, ao_num + do nu = 1, ao_num + fock_3e_uhf_ao_a(mu,nu) += f_tmp(mu,nu) + enddo + enddo + !$OMP END CRITICAL + + deallocate(f_tmp) + !$OMP END PARALLEL + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_a =', tf - ti + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_3e_uhf_ao_b, (ao_num, ao_num)] + + BEGIN_DOC + ! + ! Equations (B6) and (B7) + ! + ! g <--> gamma + ! d <--> delta + ! e <--> eta + ! k <--> kappa + ! + END_DOC + + implicit none + integer :: g, d, e, k, mu, nu + double precision :: dm_ge_a, dm_ge_b, dm_ge + double precision :: dm_dk_a, dm_dk_b, dm_dk + double precision :: i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu + double precision :: ti, tf + double precision, allocatable :: f_tmp(:,:) + + print *, ' PROVIDING fock_3e_uhf_ao_b ...' + call wall_time(ti) + + fock_3e_uhf_ao_b = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (g, e, d, k, mu, nu, dm_ge_a, dm_ge_b, dm_ge, dm_dk_a, dm_dk_b, dm_dk, f_tmp, & + !$OMP i_mugd_nuek, i_mugd_eknu, i_mugd_knue, i_mugd_nuke, i_mugd_enuk, i_mugd_kenu) & + !$OMP SHARED (ao_num, TCSCF_bi_ort_dm_ao_alpha, TCSCF_bi_ort_dm_ao_beta, fock_3e_uhf_ao_b) + + allocate(f_tmp(ao_num,ao_num)) + f_tmp = 0.d0 + + !$OMP DO + do g = 1, ao_num + do e = 1, ao_num + dm_ge_a = TCSCF_bi_ort_dm_ao_alpha(g,e) + dm_ge_b = TCSCF_bi_ort_dm_ao_beta (g,e) + dm_ge = dm_ge_a + dm_ge_b + do d = 1, ao_num + do k = 1, ao_num + dm_dk_a = TCSCF_bi_ort_dm_ao_alpha(d,k) + dm_dk_b = TCSCF_bi_ort_dm_ao_beta (d,k) + dm_dk = dm_dk_a + dm_dk_b + do mu = 1, ao_num + do nu = 1, ao_num + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, e, k, i_mugd_nuek) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, k, nu, i_mugd_eknu) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, nu, e, i_mugd_knue) + call give_integrals_3_body_bi_ort_ao(mu, g, d, nu, k, e, i_mugd_nuke) + call give_integrals_3_body_bi_ort_ao(mu, g, d, e, nu, k, i_mugd_enuk) + call give_integrals_3_body_bi_ort_ao(mu, g, d, k, e, nu, i_mugd_kenu) + f_tmp(mu,nu) -= 0.5d0 * ( dm_ge * dm_dk * i_mugd_nuek & + + dm_ge_b * dm_dk_b * i_mugd_eknu & + + dm_ge_b * dm_dk_b * i_mugd_knue & + - dm_ge_b * dm_dk * i_mugd_enuk & + - dm_ge * dm_dk_b * i_mugd_kenu & + - dm_ge_b * dm_dk_b * i_mugd_nuke & + - dm_ge_a * dm_dk_a * i_mugd_nuke ) + enddo + enddo + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do mu = 1, ao_num + do nu = 1, ao_num + fock_3e_uhf_ao_b(mu,nu) += f_tmp(mu,nu) + enddo + enddo + !$OMP END CRITICAL + + deallocate(f_tmp) + !$OMP END PARALLEL + + call wall_time(tf) + print *, ' total Wall time for fock_3e_uhf_ao_b =', tf - ti + +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..6796666d --- /dev/null +++ b/src/tc_scf/fock_tc.irp.f @@ -0,0 +1,307 @@ + +! --- + + BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_seq_alpha, (ao_num, ao_num)] +&BEGIN_PROVIDER [ double precision, two_e_tc_non_hermit_integral_seq_beta , (ao_num, ao_num)] + + BEGIN_DOC + ! + ! two_e_tc_non_hermit_integral_seq_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 + double precision :: t0, t1 + + !print*, ' providing two_e_tc_non_hermit_integral_seq ...' + !call wall_time(t0) + + two_e_tc_non_hermit_integral_seq_alpha = 0.d0 + two_e_tc_non_hermit_integral_seq_beta = 0.d0 + + do i = 1, ao_num + do k = 1, ao_num + 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_seq_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_seq_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_seq_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_seq_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) + + !! 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) + + ! rho(l,j) * < k l| T | i j> + two_e_tc_non_hermit_integral_seq_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j) + ! rho(l,j) * < k l| T | i j> + two_e_tc_non_hermit_integral_seq_beta (k,i) += density * ao_two_e_tc_tot(k,i,l,j) + ! rho_a(l,j) * < k l| T | j i> + two_e_tc_non_hermit_integral_seq_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i) + ! rho_b(l,j) * < k l| T | j i> + two_e_tc_non_hermit_integral_seq_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i) + + enddo + enddo + enddo + enddo + + !call wall_time(t1) + !print*, ' wall time for two_e_tc_non_hermit_integral_seq after = ', t1 - t0 + +END_PROVIDER + +! --- + + 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, I_coul, I_kjli + double precision :: t0, t1 + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + + !print*, ' providing two_e_tc_non_hermit_integral ...' + !call wall_time(t0) + + two_e_tc_non_hermit_integral_alpha = 0.d0 + two_e_tc_non_hermit_integral_beta = 0.d0 + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli) & + !$OMP SHARED (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_tc_tot, & + !$OMP two_e_tc_non_hermit_integral_alpha, two_e_tc_non_hermit_integral_beta) + + allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num)) + tmp_a = 0.d0 + tmp_b = 0.d0 + + !$OMP DO + 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 + do i = 1, ao_num + do k = 1, ao_num + + I_coul = density * ao_two_e_tc_tot(k,i,l,j) + I_kjli = ao_two_e_tc_tot(k,j,l,i) + + tmp_a(k,i) += I_coul - density_a * I_kjli + tmp_b(k,i) += I_coul - density_b * I_kjli + enddo + enddo + enddo + enddo + !$OMP END DO NOWAIT + + !$OMP CRITICAL + do i = 1, ao_num + do j = 1, ao_num + two_e_tc_non_hermit_integral_alpha(j,i) += tmp_a(j,i) + two_e_tc_non_hermit_integral_beta (j,i) += tmp_b(j,i) + enddo + enddo + !$OMP END CRITICAL + + deallocate(tmp_a, tmp_b) + !$OMP END PARALLEL + + !call wall_time(t1) + !print*, ' wall time for two_e_tc_non_hermit_integral after = ', t1 - t0 + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_alpha, (ao_num, ao_num)] + + BEGIN_DOC + ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the AO basis + END_DOC + + implicit none + + 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_mo_alpha, (mo_num, mo_num) ] + + BEGIN_DOC + ! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis + END_DOC + + implicit none + double precision, allocatable :: tmp(:,:) + + if(bi_ortho) then + + !allocate(tmp(ao_num,ao_num)) + !tmp = Fock_matrix_tc_ao_alpha + !if(three_body_h_tc) then + ! tmp += fock_3e_uhf_ao_a + !endif + !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1)) + !deallocate(tmp) + + 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) ) + if(three_body_h_tc) then + !Fock_matrix_tc_mo_alpha += fock_a_tot_3e_bi_orth + Fock_matrix_tc_mo_alpha += fock_3e_uhf_mo_a + endif + + 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) ] + + BEGIN_DOC + ! Total beta TC Fock matrix : h_c + Two-e^TC terms on the MO basis + END_DOC + + implicit none + double precision, allocatable :: tmp(:,:) + + if(bi_ortho) then + + !allocate(tmp(ao_num,ao_num)) + !tmp = Fock_matrix_tc_ao_beta + !if(three_body_h_tc) then + ! tmp += fock_3e_uhf_ao_b + !endif + !call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1)) + !deallocate(tmp) + + 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) ) + if(three_body_h_tc) then + !Fock_matrix_tc_mo_beta += fock_b_tot_3e_bi_orth + Fock_matrix_tc_mo_beta += fock_3e_uhf_mo_b + endif + + 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, grad_non_hermit_left] +&BEGIN_PROVIDER [ double precision, grad_non_hermit_right] +&BEGIN_PROVIDER [ double precision, grad_non_hermit] + + implicit none + integer :: i, k + + grad_non_hermit_left = 0.d0 + grad_non_hermit_right = 0.d0 + + do i = 1, elec_beta_num ! doc --> SOMO + do k = elec_beta_num+1, elec_alpha_num + grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) + grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) + !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) + !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) + !grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) + !grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) + enddo + enddo + + do i = 1, elec_beta_num ! doc --> virt + do k = elec_alpha_num+1, mo_num + grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) + grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) + !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) + !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) + grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) + grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) + enddo + enddo + + do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt + do k = elec_alpha_num+1, mo_num + grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i))) + grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k))) + !grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i)) + !grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k)) + grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i) + grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k) + enddo + enddo + + !grad_non_hermit = dsqrt(grad_non_hermit_left) + dsqrt(grad_non_hermit_right) + grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, Fock_matrix_tc_ao_tot, (ao_num, ao_num) ] + + implicit none + + call mo_to_ao_bi_ortho( Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) & + , Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) ) + +END_PROVIDER + +! --- + + diff --git a/src/tc_scf/fock_tc_mo_tot.irp.f b/src/tc_scf/fock_tc_mo_tot.irp.f new file mode 100644 index 00000000..2f33cd17 --- /dev/null +++ b/src/tc_scf/fock_tc_mo_tot.irp.f @@ -0,0 +1,144 @@ + + BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_tot, (mo_num,mo_num) ] +&BEGIN_PROVIDER [ double precision, Fock_matrix_tc_diag_mo_tot, (mo_num)] + implicit none + BEGIN_DOC + ! Fock matrix on the MO basis. + ! For open shells, the ROHF Fock Matrix is :: + ! + ! | F-K | F + K/2 | F | + ! |---------------------------------| + ! | F + K/2 | F | F - K/2 | + ! |---------------------------------| + ! | F | F - K/2 | F + K | + ! + ! + ! F = 1/2 (Fa + Fb) + ! + ! K = Fb - Fa + ! + END_DOC + integer :: i,j,n + if (elec_alpha_num == elec_beta_num) then + Fock_matrix_tc_mo_tot = Fock_matrix_tc_mo_alpha + else + + do j=1,elec_beta_num + ! F-K + do i=1,elec_beta_num !CC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F+K/2 + do i=elec_beta_num+1,elec_alpha_num !CA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F + do i=elec_alpha_num+1, mo_num !CV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + enddo + + do j=elec_beta_num+1,elec_alpha_num + ! F+K/2 + do i=1,elec_beta_num !AC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + + 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F + do i=elec_beta_num+1,elec_alpha_num !AA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_alpha_num+1, mo_num !AV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + enddo + + do j=elec_alpha_num+1, mo_num + ! F + do i=1,elec_beta_num !VC + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) + enddo + ! F-K/2 + do i=elec_beta_num+1,elec_alpha_num !VA + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j))& + - 0.5d0*(Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + ! F+K + do i=elec_alpha_num+1,mo_num !VV + Fock_matrix_tc_mo_tot(i,j) = 0.5d0*(Fock_matrix_tc_mo_alpha(i,j)+Fock_matrix_tc_mo_beta(i,j)) & + + (Fock_matrix_tc_mo_beta(i,j) - Fock_matrix_tc_mo_alpha(i,j)) + enddo + enddo + if(three_body_h_tc)then + ! C-O + do j = 1, elec_beta_num + do i = elec_beta_num+1, elec_alpha_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + ! C-V + do j = 1, elec_beta_num + do i = elec_alpha_num+1, mo_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + ! O-V + do j = elec_beta_num+1, elec_alpha_num + do i = elec_alpha_num+1, mo_num + Fock_matrix_tc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j)) + Fock_matrix_tc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i)) + enddo + enddo + endif + + endif + + do i = 1, mo_num + Fock_matrix_tc_diag_mo_tot(i) = Fock_matrix_tc_mo_tot(i,i) + enddo + + + if(frozen_orb_scf)then + integer :: iorb,jorb + do i = 1, n_core_orb + iorb = list_core(i) + do j = 1, n_act_orb + jorb = list_act(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + + if(no_oa_or_av_opt)then + do i = 1, n_act_orb + iorb = list_act(i) + do j = 1, n_inact_orb + jorb = list_inact(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_virt_orb + jorb = list_virt(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + do j = 1, n_core_orb + jorb = list_core(j) + Fock_matrix_tc_mo_tot(iorb,jorb) = 0.d0 + Fock_matrix_tc_mo_tot(jorb,iorb) = 0.d0 + enddo + enddo + endif + if(.not.bi_ortho .and. three_body_h_tc)then + Fock_matrix_tc_mo_tot += fock_3_mat + endif + +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..424eeffd --- /dev/null +++ b/src/tc_scf/fock_three.irp.f @@ -0,0 +1,229 @@ +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 + else if(bi_ortho.and.three_body_h_tc)then +!! !$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 + double precision :: integral_aaa, hthree, integral_aab, integral_abb, integral_bbb + + PROVIDE mo_l_coef mo_r_coef + + !print *, ' providing diag_three_elem_hf' + + if(.not. three_body_h_tc) then + + diag_three_elem_hf = 0.d0 + + else + + if(.not. bi_ortho) 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 + + provide mo_l_coef mo_r_coef + call give_aaa_contrib(integral_aaa) + call give_aab_contrib(integral_aab) + call give_abb_contrib(integral_abb) + call give_bbb_contrib(integral_bbb) + diag_three_elem_hf = integral_aaa + integral_aab + integral_abb + integral_bbb +! print*,'integral_aaa + integral_aab + integral_abb + integral_bbb' +! print*,integral_aaa , integral_aab , integral_abb , integral_bbb + + endif + + 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..279670b8 --- /dev/null +++ b/src/tc_scf/fock_three_bi_ortho.irp.f @@ -0,0 +1,178 @@ +BEGIN_PROVIDER [ double precision, fock_a_abb_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_abb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,beta,beta contribution + END_DOC + fock_a_abb_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_23_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + ! see contrib_3e_soo + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + fock_a_abb_3e_bi_orth_old(a,i) += direct_int - exch_23_int + enddo + enddo + + enddo + enddo + fock_a_abb_3e_bi_orth_old = - fock_a_abb_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_a_aba_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_aba_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,beta contribution + END_DOC + fock_a_aba_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num ! a + do k = 1, elec_beta_num ! b + ! a b a a b a + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + fock_a_aba_3e_bi_orth_old(a,i) += direct_int - exch_13_int + enddo + enddo + + enddo + enddo + fock_a_aba_3e_bi_orth_old = - fock_a_aba_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_a_aaa_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_a_aaa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution + END_DOC + fock_a_aaa_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + ! positive terms :: cycle contrib + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + fock_a_aaa_3e_bi_orth_old(a,i) += direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + fock_a_aaa_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int - exch_12_int + enddo + enddo + + enddo + enddo + fock_a_aaa_3e_bi_orth_old = - fock_a_aaa_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC + ! fock_a_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions + END_DOC + fock_a_tot_3e_bi_orth_old = fock_a_abb_3e_bi_orth_old + fock_a_aba_3e_bi_orth_old + fock_a_aaa_3e_bi_orth_old + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_baa_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_baa_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,alpha contribution + END_DOC + fock_b_baa_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_23_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + fock_b_baa_3e_bi_orth_old(a,i) += direct_int - exch_23_int + enddo + enddo + + enddo + enddo + fock_b_baa_3e_bi_orth_old = - fock_b_baa_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_bab_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_bab_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for beta electrons from beta,alpha,beta contribution + END_DOC + fock_b_bab_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = 1, elec_alpha_num + ! b a b b a b + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + fock_b_bab_3e_bi_orth_old(a,i) += direct_int - exch_13_int + enddo + enddo + + enddo + enddo + fock_b_bab_3e_bi_orth_old = - fock_b_bab_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_bbb_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC +! fock_b_bbb_3e_bi_orth_old(a,i) = bi-ortho 3-e Fock matrix for alpha electrons from alpha,alpha,alpha contribution + END_DOC + fock_b_bbb_3e_bi_orth_old = 0.d0 + integer :: i,a,j,k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + ! positive terms :: cycle contrib + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + fock_b_bbb_3e_bi_orth_old(a,i) += direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + fock_b_bbb_3e_bi_orth_old(a,i) += - exch_13_int - exch_23_int - exch_12_int + enddo + enddo + + enddo + enddo + fock_b_bbb_3e_bi_orth_old = - fock_b_bbb_3e_bi_orth_old +END_PROVIDER + +BEGIN_PROVIDER [ double precision, fock_b_tot_3e_bi_orth_old, (mo_num, mo_num)] + implicit none + BEGIN_DOC + ! fock_b_tot_3e_bi_orth_old = bi-ortho 3-e Fock matrix for alpha electrons from all possible spin contributions + END_DOC + fock_b_tot_3e_bi_orth_old = fock_b_bbb_3e_bi_orth_old + fock_b_bab_3e_bi_orth_old + fock_b_baa_3e_bi_orth_old + +END_PROVIDER diff --git a/src/tc_scf/fock_three_bi_ortho_new_new.irp.f b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f new file mode 100644 index 00000000..f73171a3 --- /dev/null +++ b/src/tc_scf/fock_three_bi_ortho_new_new.irp.f @@ -0,0 +1,286 @@ + +! --- + +BEGIN_PROVIDER [double precision, fock_a_tot_3e_bi_orth, (mo_num, mo_num)] + + implicit none + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_a_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp1_bi_ortho(a,i) + fock_a_tot_3e_bi_orth(a,i) += fock_a_tmp2_bi_ortho(a,i) + enddo + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_b_tot_3e_bi_orth, (mo_num, mo_num)] + + implicit none + integer :: i, a + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tot_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + fock_b_tot_3e_bi_orth(a,i) += fock_cs_3e_bi_orth (a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp2_bi_ortho(a,i) + fock_b_tot_3e_bi_orth(a,i) += fock_b_tmp1_bi_ortho(a,i) + enddo + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_cs_3e_bi_orth, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_cs_3e_bi_orth = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = 1, elec_beta_num + do k = 1, elec_beta_num + + !!call contrib_3e_sss(a,i,j,k,contrib_sss) + !!call contrib_3e_soo(a,i,j,k,contrib_soo) + !!call contrib_3e_sos(a,i,j,k,contrib_sos) + !!contrib = 0.5d0 * (contrib_sss + contrib_soo) + contrib_sos + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + new = 2.d0 * direct_int + 0.5d0 * (c_3_int + c_minus_3_int - exch_12_int) -1.5d0 * exch_13_int - exch_23_int + + fock_cs_3e_bi_orth(a,i) += new + enddo + enddo + enddo + enddo + + fock_cs_3e_bi_orth = - fock_cs_3e_bi_orth + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_a_tmp1_bi_ortho, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss, contrib_sos, contrib_soo, contrib + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_beta_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + + fock_a_tmp1_bi_ortho(a,i) += 1.5d0 * (direct_int - exch_13_int) + 0.5d0 * (c_3_int + c_minus_3_int - exch_23_int - exch_12_int) + enddo + enddo + enddo + enddo + + fock_a_tmp1_bi_ortho = - fock_a_tmp1_bi_ortho + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_a_tmp2_bi_ortho, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: contrib_sss + + PROVIDE mo_l_coef mo_r_coef + + fock_a_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_alpha_num + do k = elec_beta_num+1, elec_alpha_num + call contrib_3e_sss(a, i, j, k, contrib_sss) + + fock_a_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_sss + enddo + enddo + enddo + enddo + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_b_tmp1_bi_ortho, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int + double precision :: new + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tmp1_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = 1, elec_beta_num + do k = elec_beta_num+1, elec_alpha_num + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + + fock_b_tmp1_bi_ortho(a,i) += 1.5d0 * direct_int - 0.5d0 * exch_23_int - exch_13_int + enddo + enddo + enddo + enddo + + fock_b_tmp1_bi_ortho = - fock_b_tmp1_bi_ortho + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, fock_b_tmp2_bi_ortho, (mo_num, mo_num)] + + implicit none + integer :: i, a, j, k + double precision :: contrib_soo + + PROVIDE mo_l_coef mo_r_coef + + fock_b_tmp2_bi_ortho = 0.d0 + + do i = 1, mo_num + do a = 1, mo_num + do j = elec_beta_num + 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_soo(a, i, j, k, contrib_soo) + + fock_b_tmp2_bi_ortho(a,i) += 0.5d0 * contrib_soo + enddo + enddo + enddo + enddo + +END_PROVIDER + +! --- + +subroutine contrib_3e_sss(a, i, j, k, integral) + + BEGIN_DOC + ! returns the pure same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + + PROVIDE mo_l_coef mo_r_coef + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )!!! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, i, k, c_3_int) ! < a k j | j i k > + call give_integrals_3_body_bi_ort(a, k, j, k, j, i, c_minus_3_int)! < a k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)!!! < a k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)!!! < a k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(a, k, j, k, i, j, exch_12_int)!!! < a k j | k i j > : E_12 + integral += - exch_13_int - exch_23_int - exch_12_int + + integral = -integral + +end + +! --- + +subroutine contrib_3e_soo(a,i,j,k,integral) + + BEGIN_DOC + ! returns the same spin / opposite spin / opposite spin contribution to F(a,i) from two orbitals j,k + END_DOC + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_23_int + + PROVIDE mo_l_coef mo_r_coef + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int) ! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, i, j, k, exch_23_int)! < a k j | i j k > : E_23 + integral = direct_int - exch_23_int + + integral = -integral + +end + +! --- + +subroutine contrib_3e_sos(a, i, j, k, integral) + + BEGIN_DOC + ! returns the same spin / opposite spin / same spin contribution to F(a,i) from two orbitals j,k + END_DOC + + PROVIDE mo_l_coef mo_r_coef + + implicit none + integer, intent(in) :: a, i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int + + call give_integrals_3_body_bi_ort(a, k, j, i, k, j, direct_int )! < a k j | i k j > + call give_integrals_3_body_bi_ort(a, k, j, j, k, i, exch_13_int)! < a k j | j k i > : E_13 + integral = direct_int - exch_13_int + + integral = -integral + +end + +! --- + 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/minimize_tc_angles.irp.f b/src/tc_scf/minimize_tc_angles.irp.f new file mode 100644 index 00000000..cb729eb2 --- /dev/null +++ b/src/tc_scf/minimize_tc_angles.irp.f @@ -0,0 +1,12 @@ +program print_angles + implicit none + 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 = 14 ! small grid for quick debug + touch my_n_pt_r_grid my_n_pt_a_grid +! call sort_by_tc_fock + call minimize_tc_orb_angles +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/print_angle_tc_orb.irp.f b/src/tc_scf/print_angle_tc_orb.irp.f new file mode 100644 index 00000000..09260395 --- /dev/null +++ b/src/tc_scf/print_angle_tc_orb.irp.f @@ -0,0 +1,9 @@ +program print_angles + implicit none + 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 = 14 ! small grid for quick debug + call print_angles_tc +end diff --git a/src/tc_scf/print_fit_param.irp.f b/src/tc_scf/print_fit_param.irp.f new file mode 100644 index 00000000..f8bcfa7f --- /dev/null +++ b/src/tc_scf/print_fit_param.irp.f @@ -0,0 +1,60 @@ +program print_fit_param + + 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 +! 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 main() + +end + +! --- + +subroutine main() + + implicit none + integer :: i + + mu_erf = 1.d0 + touch mu_erf + + print *, ' fit for (1 - erf(x))^2' + do i = 1, n_max_fit_slat + print*, expo_gauss_1_erf_x_2(i), coef_gauss_1_erf_x_2(i) + enddo + + print *, '' + print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]' + do i = 1, n_max_fit_slat + print *, expo_gauss_j_mu_x(i), 2.d0 * coef_gauss_j_mu_x(i) + enddo + + print *, '' + print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)]^2' + do i = 1, n_max_fit_slat + print *, expo_gauss_j_mu_x_2(i), 4.d0 * coef_gauss_j_mu_x_2(i) + enddo + + print *, '' + print *, ' fit for [x * (1 - erf(x)) - 1/sqrt(pi) * exp(-x**2)] x [1 - erf(mu * r12)]' + do i = 1, n_max_fit_slat + print *, expo_gauss_j_mu_1_erf(i), 4.d0 * coef_gauss_j_mu_1_erf(i) + enddo + + return +end subroutine main + +! --- + diff --git a/src/tc_scf/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f new file mode 100644 index 00000000..0312df5f --- /dev/null +++ b/src/tc_scf/rh_tcscf.irp.f @@ -0,0 +1,336 @@ +! --- + +subroutine rh_tcscf() + + BEGIN_DOC + ! + ! Roothaan-Hall algorithm for TC-SCF calculation + ! + END_DOC + + implicit none + + integer :: i, j + integer :: iteration_TCSCF, dim_DIIS, index_dim_DIIS + double precision :: energy_TCSCF, energy_TCSCF_1e, energy_TCSCF_2e, energy_TCSCF_3e, gradie_TCSCF + double precision :: energy_TCSCF_previous, delta_energy_TCSCF + double precision :: gradie_TCSCF_previous, delta_gradie_TCSCF + double precision :: max_error_DIIS_TCSCF + double precision :: level_shift_save + double precision :: delta_energy_tmp, delta_gradie_tmp + double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:) + double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:) + + logical, external :: qp_stop + + + !PROVIDE ao_md5 mo_occ + PROVIDE level_shift_TCSCF + + allocate( mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num) & + , F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), e_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF) ) + + F_DIIS = 0.d0 + e_DIIS = 0.d0 + mo_l_coef_save = 0.d0 + mo_r_coef_save = 0.d0 + + call write_time(6) + + ! --- + ! Initialize energies and density matrices + + energy_TCSCF_previous = TC_HF_energy + energy_TCSCF_1e = TC_HF_one_e_energy + energy_TCSCF_2e = TC_HF_two_e_energy + energy_TCSCF_3e = 0.d0 + if(three_body_h_tc) then + energy_TCSCF_3e = diag_three_elem_hf + endif + gradie_TCSCF_previous = grad_non_hermit + delta_energy_TCSCF = 1.d0 + delta_gradie_TCSCF = 1.d0 + iteration_TCSCF = 0 + dim_DIIS = 0 + max_error_DIIS_TCSCF = 1.d0 + + ! --- + + ! Start of main SCF loop + + PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot + + do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. & + !(dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. & + (dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF)) ) + + iteration_TCSCF += 1 + if(iteration_TCSCF > n_it_TCSCF_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + stop + endif + + dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) + + ! --- + + if((tcscf_algorithm == 'DIIS') .and. (dabs(delta_energy_TCSCF) > 1.d-6)) then + + ! store Fock and error matrices at each iteration + index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1 + do j = 1, ao_num + do i = 1, ao_num + F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j) + e_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao(i,j) + enddo + enddo + + call extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), iteration_TCSCF, dim_DIIS) + + Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot + Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot + !TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta + + 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) ) + 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) ) + TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta + endif + + ! --- + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + TOUCH mo_l_coef mo_r_coef + + ! --- + + ! calculate error vectors + max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo)) + + ! --- + + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous + + ! --- + + do while((delta_gradie_tmp > 1.d-7) .and. (iteration_TCSCF > 1)) + !do while((dabs(delta_energy_tmp) > 0.5d0) .and. (iteration_TCSCF > 1)) + print *, ' very big or bad step : ', delta_energy_tmp, delta_gradie_tmp + print *, ' TC level shift = ', level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) + + if(level_shift_TCSCF <= .1d0) then + level_shift_TCSCF = 1.d0 + else + level_shift_TCSCF = level_shift_TCSCF * 3.0d0 + endif + TOUCH mo_l_coef mo_r_coef level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + TOUCH mo_l_coef mo_r_coef + + delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous + delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous + + if(level_shift_TCSCF - level_shift_save > 40.d0) then + level_shift_TCSCF = level_shift_save * 4.d0 + SOFT_TOUCH level_shift_TCSCF + exit + endif + + dim_DIIS = 0 + enddo +! print *, ' very big step : ', delta_energy_tmp +! print *, ' TC level shift = ', level_shift_TCSCF + + ! --- + + level_shift_TCSCF = 0.d0 + !level_shift_TCSCF = level_shift_TCSCF * 0.5d0 + SOFT_TOUCH level_shift_TCSCF + + gradie_TCSCF = grad_non_hermit + energy_TCSCF = TC_HF_energy + energy_TCSCF_1e = TC_HF_one_e_energy + energy_TCSCF_2e = TC_HF_two_e_energy + energy_TCSCF_3e = 0.d0 + if(three_body_h_tc) then + energy_TCSCF_3e = diag_three_elem_hf + endif + delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous + delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous + + energy_TCSCF_previous = energy_TCSCF + gradie_TCSCF_previous = gradie_TCSCF + + + level_shift_save = level_shift_TCSCF + mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) + mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) + + + print *, ' iteration = ', iteration_TCSCF + print *, ' total TC energy = ', energy_TCSCF + print *, ' 1-e TC energy = ', energy_TCSCF_1e + print *, ' 2-e TC energy = ', energy_TCSCF_2e + print *, ' 3-e TC energy = ', energy_TCSCF_3e + print *, ' |delta TC energy| = ', dabs(delta_energy_TCSCF) + print *, ' TC gradient = ', gradie_TCSCF + print *, ' delta TC gradient = ', delta_gradie_TCSCF + print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF + print *, ' TC DIIS dim = ', dim_DIIS + print *, ' TC level shift = ', level_shift_TCSCF + print *, ' ' + + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + + if(qp_stop()) exit + enddo + + ! --- + + print *, ' TCSCF DIIS converged !' + call print_energy_and_mos() + + call write_time(6) + + deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS) + +end + +! --- + +subroutine extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, F_ao, size_F_ao, iteration_TCSCF, dim_DIIS) + + BEGIN_DOC + ! + ! Compute the extrapolated Fock matrix using the DIIS procedure + ! + ! e = \sum_i c_i e_i and \sum_i c_i = 1 + ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1) + ! + END_DOC + + implicit none + + integer, intent(in) :: iteration_TCSCF, size_F_ao + integer, intent(inout) :: dim_DIIS + double precision, intent(in) :: F_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(in) :: e_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(inout) :: F_ao(size_F_ao,ao_num) + + double precision, allocatable :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:) + + integer :: i, j, k, l, i_DIIS, j_DIIS + integer :: lwork + double precision :: rcond, ferr, berr + integer, allocatable :: iwork(:) + double precision, allocatable :: scratch(:,:) + + if(dim_DIIS < 1) then + return + endif + + allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) & + , C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) ) + + ! Compute the matrices B and X + B_matrix_DIIS(:,:) = 0.d0 + do j = 1, dim_DIIS + j_DIIS = min(dim_DIIS, mod(iteration_TCSCF-j, max_dim_DIIS_TCSCF)+1) + + do i = 1, dim_DIIS + i_DIIS = min(dim_DIIS, mod(iteration_TCSCF-i, max_dim_DIIS_TCSCF)+1) + + ! Compute product of two errors vectors + do l = 1, ao_num + do k = 1, ao_num + B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + e_DIIS(k,l,i_DIIS) * e_DIIS(k,l,j_DIIS) + enddo + enddo + + enddo + enddo + + ! Pad B matrix and build the X matrix + + C_vector_DIIS(:) = 0.d0 + do i = 1, dim_DIIS + B_matrix_DIIS(i,dim_DIIS+1) = -1.d0 + B_matrix_DIIS(dim_DIIS+1,i) = -1.d0 + enddo + C_vector_DIIS(dim_DIIS+1) = -1.d0 + + deallocate(scratch) + + ! Estimate condition number of B + integer :: info + double precision :: anorm + integer, allocatable :: ipiv(:) + double precision, allocatable :: AF(:,:) + double precision, external :: dlange + + lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5) + allocate(AF(dim_DIIS+1,dim_DIIS+1)) + allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) ) + allocate(scratch(lwork,1)) + scratch(:,1) = 0.d0 + + anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1)) + + AF(:,:) = B_matrix_DIIS(:,:) + call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + if(rcond < 1.d-14) then + dim_DIIS = 0 + return + endif + + ! solve the linear system C = B x X + + X_vector_DIIS = C_vector_DIIS + call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info) + + deallocate(scratch, AF, iwork) + if(info < 0) then + stop ' bug in TC-DIIS' + endif + + ! Compute extrapolated Fock matrix + + !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) + do j = 1, ao_num + do i = 1, ao_num + F_ao(i,j) = 0.d0 + enddo + do k = 1, dim_DIIS + if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle + do i = 1,ao_num + ! FPE here + F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + +! --- + diff --git a/src/tc_scf/rh_tcscf_diis.irp.f b/src/tc_scf/rh_tcscf_diis.irp.f new file mode 100644 index 00000000..306c78b3 --- /dev/null +++ b/src/tc_scf/rh_tcscf_diis.irp.f @@ -0,0 +1,362 @@ +! --- + +subroutine rh_tcscf_diis() + + implicit none + + integer :: i, j, it + integer :: dim_DIIS, index_dim_DIIS + double precision :: etc_tot, etc_1e, etc_2e, etc_3e, e_save, e_delta + double precision :: tc_grad, g_save, g_delta, g_delta_th + double precision :: level_shift_save, rate_th + double precision :: t0, t1 + double precision :: er_DIIS, er_delta, er_save, er_delta_th + double precision, allocatable :: F_DIIS(:,:,:), E_DIIS(:,:,:) + double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:) + + logical, external :: qp_stop + + it = 0 + e_save = 0.d0 + dim_DIIS = 0 + g_delta_th = 1d0 + er_delta_th = 1d0 + rate_th = 100.d0 !0.01d0 !0.2d0 + + allocate(mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num)) + mo_l_coef_save = 0.d0 + mo_r_coef_save = 0.d0 + + allocate(F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), E_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF)) + F_DIIS = 0.d0 + E_DIIS = 0.d0 + + call write_time(6) + + ! --- + + PROVIDE level_shift_TCSCF + PROVIDE mo_l_coef mo_r_coef + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '================', '====', '========' + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & + , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '================', '====', '========' + + + ! first iteration (HF orbitals) + call wall_time(t0) + + etc_tot = TC_HF_energy + etc_1e = TC_HF_one_e_energy + etc_2e = TC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + tc_grad = grad_non_hermit + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + + e_save = etc_tot + g_save = tc_grad + er_save = er_DIIS + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + + ! --- + + PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot + + do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. threshold_DIIS_nonzero_TCSCF)) + + call wall_time(t0) + + it += 1 + if(it > n_it_TCSCF_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + stop + endif + + dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF) + + ! --- + + if(dabs(e_delta) > 1.d-12) then + + index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1 + do j = 1, ao_num + do i = 1, ao_num + F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j) + E_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao (i,j) + enddo + enddo + + call extrapolate_TC_Fock_matrix(E_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), it, dim_DIIS) + + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & + , Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) ) + TOUCH Fock_matrix_tc_mo_tot fock_matrix_tc_diag_mo_tot + endif + + ! --- + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + !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 + + ! --- + + g_delta = grad_non_hermit - g_save + er_delta = maxval(abs(FQS_SQF_mo)) - er_save + + !if((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) then + if((g_delta > rate_th * g_delta_th) .and. (it > 1)) then + !if((g_delta > 0.d0) .and. (it > 1)) then + + Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS) + call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) & + , Fock_matrix_tc_mo_tot, size(Fock_matrix_tc_mo_tot, 1) ) + TOUCH Fock_matrix_tc_mo_tot fock_matrix_tc_diag_mo_tot + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + !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 + + endif + + ! --- + + g_delta = grad_non_hermit - g_save + er_delta = maxval(abs(FQS_SQF_mo)) - er_save + + mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) + mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) + + !do while((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) + do while((g_delta > rate_th * g_delta_th) .and. (it > 1)) + print *, ' big or bad step : ', g_delta, rate_th * g_delta_th + + mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) + if(level_shift_TCSCF <= .1d0) then + level_shift_TCSCF = 1.d0 + else + level_shift_TCSCF = level_shift_TCSCF * 3.0d0 + endif + TOUCH mo_l_coef mo_r_coef level_shift_TCSCF + + mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num) + mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num) + !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 + + g_delta = grad_non_hermit - g_save + er_delta = maxval(abs(FQS_SQF_mo)) - er_save + + if(level_shift_TCSCF - level_shift_save > 40.d0) then + level_shift_TCSCF = level_shift_save * 4.d0 + SOFT_TOUCH level_shift_TCSCF + exit + endif + + dim_DIIS = 0 + enddo + + ! --- + + level_shift_TCSCF = level_shift_TCSCF * 0.5d0 + SOFT_TOUCH level_shift_TCSCF + + etc_tot = TC_HF_energy + etc_1e = TC_HF_one_e_energy + etc_2e = TC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + tc_grad = grad_non_hermit + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + g_delta = tc_grad - g_save + er_delta = er_DIIS - er_save + + e_save = etc_tot + g_save = tc_grad + level_shift_save = level_shift_TCSCF + er_save = er_DIIS + + g_delta_th = dabs(tc_grad) ! g_delta) + er_delta_th = dabs(er_DIIS) !er_delta) + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + + if(g_delta .lt. 0.d0) then + call ezfio_set_tc_scf_bitc_energy(etc_tot) + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + endif + + if(qp_stop()) exit + enddo + + ! --- + + print *, ' TCSCF DIIS converged !' + call print_energy_and_mos() + + call write_time(6) + + deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, E_DIIS) + + call ezfio_set_tc_scf_bitc_energy(TC_HF_energy) + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + +end + +! --- + +subroutine extrapolate_TC_Fock_matrix(E_DIIS, F_DIIS, F_ao, size_F_ao, it, dim_DIIS) + + BEGIN_DOC + ! + ! Compute the extrapolated Fock matrix using the DIIS procedure + ! + ! e = \sum_i c_i e_i and \sum_i c_i = 1 + ! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1) + ! + END_DOC + + implicit none + + integer, intent(in) :: it, size_F_ao + integer, intent(inout) :: dim_DIIS + double precision, intent(in) :: F_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(in) :: E_DIIS(ao_num,ao_num,dim_DIIS) + double precision, intent(inout) :: F_ao(size_F_ao,ao_num) + + double precision, allocatable :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:) + + integer :: i, j, k, l, i_DIIS, j_DIIS + integer :: lwork + double precision :: rcond, ferr, berr + integer, allocatable :: iwork(:) + double precision, allocatable :: scratch(:,:) + + if(dim_DIIS < 1) then + return + endif + + allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) & + , C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) ) + + ! Compute the matrices B and X + B_matrix_DIIS(:,:) = 0.d0 + do j = 1, dim_DIIS + j_DIIS = min(dim_DIIS, mod(it-j, max_dim_DIIS_TCSCF)+1) + + do i = 1, dim_DIIS + i_DIIS = min(dim_DIIS, mod(it-i, max_dim_DIIS_TCSCF)+1) + + ! Compute product of two errors vectors + do l = 1, ao_num + do k = 1, ao_num + B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + E_DIIS(k,l,i_DIIS) * E_DIIS(k,l,j_DIIS) + enddo + enddo + + enddo + enddo + + ! Pad B matrix and build the X matrix + + C_vector_DIIS(:) = 0.d0 + do i = 1, dim_DIIS + B_matrix_DIIS(i,dim_DIIS+1) = -1.d0 + B_matrix_DIIS(dim_DIIS+1,i) = -1.d0 + enddo + C_vector_DIIS(dim_DIIS+1) = -1.d0 + + deallocate(scratch) + + ! Estimate condition number of B + integer :: info + double precision :: anorm + integer, allocatable :: ipiv(:) + double precision, allocatable :: AF(:,:) + double precision, external :: dlange + + lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5) + allocate(AF(dim_DIIS+1,dim_DIIS+1)) + allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) ) + allocate(scratch(lwork,1)) + scratch(:,1) = 0.d0 + + anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1)) + + AF(:,:) = B_matrix_DIIS(:,:) + call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info) + if(info /= 0) then + dim_DIIS = 0 + return + endif + + if(rcond < 1.d-14) then + dim_DIIS = 0 + return + endif + + ! solve the linear system C = B x X + + X_vector_DIIS = C_vector_DIIS + call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info) + + deallocate(scratch, AF, iwork) + if(info < 0) then + stop ' bug in TC-DIIS' + endif + + ! Compute extrapolated Fock matrix + + !$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200) + do j = 1, ao_num + do i = 1, ao_num + F_ao(i,j) = 0.d0 + enddo + do k = 1, dim_DIIS + if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle + do i = 1,ao_num + ! FPE here + F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1) + enddo + enddo + enddo + !$OMP END PARALLEL DO + +end + +! --- + diff --git a/src/tc_scf/rh_tcscf_simple.irp.f b/src/tc_scf/rh_tcscf_simple.irp.f new file mode 100644 index 00000000..30798e3d --- /dev/null +++ b/src/tc_scf/rh_tcscf_simple.irp.f @@ -0,0 +1,129 @@ +! --- + +subroutine rh_tcscf_simple() + + implicit none + integer :: i, j, it, dim_DIIS + double precision :: t0, t1 + double precision :: e_save, e_delta, rho_delta + double precision :: etc_tot, etc_1e, etc_2e, etc_3e, tc_grad + double precision :: er_DIIS + double precision, allocatable :: rho_old(:,:), rho_new(:,:) + + allocate(rho_old(ao_num,ao_num), rho_new(ao_num,ao_num)) + + it = 0 + e_save = 0.d0 + dim_DIIS = 0 + + ! --- + + if(.not. bi_ortho) then + print *, ' grad_hermit = ', grad_hermit + call save_good_hermit_tc_eigvectors + TOUCH mo_coef + call save_mos + endif + + ! --- + + if(bi_ortho) then + + PROVIDE level_shift_tcscf + PROVIDE mo_l_coef mo_r_coef + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '================', '====', '========' + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' & + , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)' + + write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') & + '====', '================', '================', '================', '================', '================' & + , '================', '================', '================', '====', '========' + + + ! first iteration (HF orbitals) + call wall_time(t0) + + etc_tot = TC_HF_energy + etc_1e = TC_HF_one_e_energy + etc_2e = TC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + tc_grad = grad_non_hermit + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + e_save = etc_tot + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + + do while(tc_grad .gt. dsqrt(thresh_tcscf)) + call wall_time(t0) + + it += 1 + if(it > n_it_tcscf_max) then + print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max + stop + endif + + mo_l_coef = fock_tc_leigvec_ao + mo_r_coef = fock_tc_reigvec_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 + + etc_tot = TC_HF_energy + etc_1e = TC_HF_one_e_energy + etc_2e = TC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + etc_3e = diag_three_elem_hf + endif + tc_grad = grad_non_hermit + er_DIIS = maxval(abs(FQS_SQF_mo)) + e_delta = dabs(etc_tot - e_save) + e_save = etc_tot + + call ezfio_set_tc_scf_bitc_energy(etc_tot) + + call wall_time(t1) + write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') & + it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0 + 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_e_energy + print *, 'TC HF 2 e energy = ', TC_HF_two_e_energy + print *, 'TC HF 3 body = ', diag_three_elem_hf + print *, '***' + print *, '' + call save_good_hermit_tc_eigvectors + TOUCH mo_coef + call save_mos + enddo + + endif + + print *, ' TCSCF Simple converged !' + call print_energy_and_mos() + + deallocate(rho_old, rho_new) + +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..fc4a7935 --- /dev/null +++ b/src/tc_scf/rotate_tcscf_orbitals.irp.f @@ -0,0 +1,367 @@ + +! --- + +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_to_maximize_overlap(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_to_maximize_overlap(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, tot_deg + 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 = thr_degen_tc + + 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 + + tot_deg = 0 + do i = 1, m + if(deg_num(i).gt.1) then + print *, ' degen on', i, deg_num(i) + tot_deg = tot_deg + 1 + endif + enddo + + if(tot_deg .eq. 0) then + print *, ' no degen' + return + endif + + ! --- + + 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_to_maximize_overlap + +! --- + +subroutine fix_right_to_one() + + implicit none + integer :: i, j, m, n, mm, tot_deg + double precision :: accu_d, accu_nd + double precision :: de_thr, ei, ej, de + integer, allocatable :: deg_num(:) + double precision, allocatable :: R0(:,:), L0(:,:), W(:,:), e0(:) + double precision, allocatable :: R(:,:), L(:,:), S(:,:), Stmp(:,:), tmp(:,:) + + n = ao_num + m = mo_num + + allocate(L0(n,m), R0(n,m), W(n,n), e0(m)) + L0 = mo_l_coef + R0 = mo_r_coef + W = ao_overlap + + print*, ' fock matrix diag elements' + do i = 1, m + e0(i) = Fock_matrix_tc_mo_tot(i,i) + print*, e0(i) + enddo + + ! --- + + allocate( deg_num(m) ) + do i = 1, m + deg_num(i) = 1 + enddo + + de_thr = 1d-6 + + 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 + + deallocate(e0) + + tot_deg = 0 + do i = 1, m + if(deg_num(i).gt.1) then + print *, ' degen on', i, deg_num(i) + tot_deg = tot_deg + 1 + endif + enddo + + if(tot_deg .eq. 0) then + print *, ' no degen' + return + endif + + ! --- + + do i = 1, m + mm = deg_num(i) + + if(mm .gt. 1) then + + allocate(L(n,mm), R(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) + enddo + + ! --- + + call impose_weighted_orthog_svd(n, mm, W, R) + call impose_weighted_biorthog_qr(n, mm, thresh_biorthog_diag, thresh_biorthog_nondiag, R, W, L) + + ! --- + + do j = 1, mm + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + deallocate(L, R) + + endif + enddo + + call check_weighted_biorthog_binormalize(n, m, L0, W, R0, thresh_biorthog_diag, thresh_biorthog_nondiag, .true.) + + deallocate(W, deg_num) + + mo_l_coef = L0 + mo_r_coef = R0 + deallocate(L0, R0) + + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + print *, ' orbitals are rotated ' + + return +end subroutine fix_right_to_one + +! --- diff --git a/src/tc_scf/routines_rotates.irp.f b/src/tc_scf/routines_rotates.irp.f new file mode 100644 index 00000000..596ae500 --- /dev/null +++ b/src/tc_scf/routines_rotates.irp.f @@ -0,0 +1,359 @@ + +! --- + +subroutine minimize_tc_orb_angles() + + implicit none + logical :: good_angles + integer :: i + double precision :: thr_deg + + good_angles = .False. + thr_deg = thr_degen_tc + + call print_energy_and_mos() + + print *, ' Minimizing the angles between the TC orbitals' + i = 1 + do while (.not. good_angles) + print *, ' iteration = ', i + call routine_save_rotated_mos(thr_deg, good_angles) + thr_deg *= 10.d0 + i += 1 + if(i .gt. 100) then + print *, ' minimize_tc_orb_angles does not seem to converge ..' + print *, ' Something is weird in the tc orbitals ...' + print *, ' STOPPING' + stop + endif + enddo + print *, ' Converged ANGLES MINIMIZATION !!' + + call print_angles_tc() + call print_energy_and_mos() + +end + +! --- + +subroutine routine_save_rotated_mos(thr_deg, good_angles) + + implicit none + + double precision, intent(in) :: thr_deg + logical, intent(out) :: good_angles + + integer :: i, j, k, n_degen_list, m, n, n_degen, ilast, ifirst + double precision :: max_angle, norm + integer, allocatable :: list_degen(:,:) + double precision, allocatable :: new_angles(:) + double precision, allocatable :: mo_r_coef_good(:,:), mo_l_coef_good(:,:) + double precision, allocatable :: mo_r_coef_new(:,:) + double precision, allocatable :: fock_diag(:),s_mat(:,:) + double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:) + double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:) + + good_angles = .False. + + allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num)) + + print *, ' ***************************************' + print *, ' ***************************************' + print *, ' THRESHOLD FOR DEGENERACIES ::: ', thr_deg + print *, ' ***************************************' + print *, ' ***************************************' + print *, ' Starting with the following TC energy gradient :', grad_non_hermit + + mo_r_coef_good = mo_r_coef + mo_l_coef_good = mo_l_coef + + allocate(mo_r_coef_new(ao_num, mo_num)) + mo_r_coef_new = mo_r_coef + do i = 1, mo_num + norm = 1.d0/dsqrt(overlap_mo_r(i,i)) + do j = 1, ao_num + mo_r_coef_new(j,i) *= norm + enddo + enddo + + allocate(list_degen(mo_num,0:mo_num), s_mat(mo_num,mo_num), fock_diag(mo_num)) + do i = 1, mo_num + fock_diag(i) = Fock_matrix_tc_mo_tot(i,i) + enddo + + ! compute the overlap between the left and rescaled right + call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat) +! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list) + call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list) + print *, ' fock_matrix_mo' + do i = 1, mo_num + print *, i, fock_diag(i), angle_left_right(i) + enddo + + do i = 1, n_degen_list +! ifirst = list_degen(1,i) +! ilast = list_degen(2,i) +! n_degen = ilast - ifirst +1 + + n_degen = list_degen(i,0) + if(n_degen .eq. 1) cycle + + allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen)) + allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen)) + allocate(T(n_degen,n_degen), Snew(n_degen,n_degen)) + + do j = 1, n_degen + mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j)) + mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j)) + enddo + ! Orthogonalization of right functions + print *, ' Orthogonalization of RIGHT functions' + print *, ' ------------------------------------' + call orthog_functions(ao_num, n_degen, mo_r_coef_tmp, ao_overlap) + + ! Orthogonalization of left functions + print *, ' Orthogonalization of LEFT functions' + print *, ' ------------------------------------' + call orthog_functions(ao_num, n_degen, mo_l_coef_tmp, ao_overlap) + + print *, ' Overlap left-right ' + call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp) + do j = 1, n_degen + write(*,'(100(F8.4,X))') stmp(:,j) + enddo + call build_s_matrix(ao_num, n_degen, mo_l_coef_tmp, mo_l_coef_tmp, ao_overlap, stmp) + + !print*,'LEFT/LEFT OVERLAP ' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp) + !print*,'RIGHT/RIGHT OVERLAP ' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + + if(maxovl_tc) then + T = 0.d0 + Snew = 0.d0 + call maxovl(n_degen, n_degen, stmp, T, Snew) + !print*,'overlap after' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')Snew(:,j) + !enddo + call dgemm( 'N', 'N', ao_num, n_degen, n_degen, 1.d0 & + , mo_l_coef_tmp, size(mo_l_coef_tmp, 1), T(1,1), size(T, 1) & + , 0.d0, mo_l_coef_new, size(mo_l_coef_new, 1) ) + call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp) + !print*,'Overlap test' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + else + mo_l_coef_new = mo_l_coef_tmp + endif + + call impose_weighted_biorthog_svd(ao_num, n_degen, ao_overlap, mo_l_coef_new, mo_r_coef_tmp) + + !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_r_coef_tmp, ao_overlap, stmp) + !print*,'LAST OVERLAP ' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + !call build_s_matrix(ao_num, n_degen, mo_l_coef_new, mo_l_coef_new, ao_overlap, stmp) + !print*,'LEFT OVERLAP ' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + !call build_s_matrix(ao_num, n_degen, mo_r_coef_tmp, mo_r_coef_tmp, ao_overlap, stmp) + !print*,'RIGHT OVERLAP ' + !do j = 1, n_degen + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + do j = 1, n_degen +!!! mo_l_coef_good(1:ao_num,j+ifirst-1) = mo_l_coef_new(1:ao_num,j) +!!! mo_r_coef_good(1:ao_num,j+ifirst-1) = mo_r_coef_tmp(1:ao_num,j) + mo_l_coef_good(1:ao_num,list_degen(i,j)) = mo_l_coef_new(1:ao_num,j) + mo_r_coef_good(1:ao_num,list_degen(i,j)) = mo_r_coef_tmp(1:ao_num,j) + enddo + + deallocate(stmp, smat2) + deallocate(mo_r_coef_tmp, mo_l_coef_tmp, mo_l_coef_new) + deallocate(T, Snew) + enddo + + !allocate(stmp(mo_num, mo_num)) + !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_r_coef_good, ao_overlap, stmp) + !print*,'LEFT/RIGHT OVERLAP ' + !do j = 1, mo_num + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + !call build_s_matrix(ao_num, mo_num, mo_l_coef_good, mo_l_coef_good, ao_overlap, stmp) + !print*,'LEFT/LEFT OVERLAP ' + !do j = 1, mo_num + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + !call build_s_matrix(ao_num, mo_num, mo_r_coef_good, mo_r_coef_good, ao_overlap, stmp) + !print*,'RIGHT/RIGHT OVERLAP ' + !do j = 1, mo_num + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + + mo_r_coef = mo_r_coef_good + mo_l_coef = mo_l_coef_good + 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 + + allocate(new_angles(mo_num)) + new_angles(1:mo_num) = dabs(angle_left_right(1:mo_num)) + max_angle = maxval(new_angles) + good_angles = max_angle.lt.45.d0 + print *, ' max_angle = ', max_angle + +end + +! --- + +subroutine build_s_matrix(m, n, C1, C2, overlap, smat) + + implicit none + integer, intent(in) :: m, n + double precision, intent(in) :: C1(m,n), C2(m,n), overlap(m,m) + double precision, intent(out) :: smat(n,n) + integer :: i, j, k, l + double precision, allocatable :: S_tmp(:,:) + + smat = 0.d0 + + !do i = 1, n + ! do j = 1, n + ! do k = 1, m + ! do l = 1, m + ! smat(i,j) += C1(k,i) * overlap(l,k) * C2(l,j) + ! enddo + ! enddo + ! enddo + !enddo + + ! C1.T x overlap + allocate(S_tmp(n,m)) + call dgemm( 'T', 'N', n, m, m, 1.d0 & + , C1, size(C1, 1), overlap, size(overlap, 1) & + , 0.d0, S_tmp, size(S_tmp, 1) ) + ! C1.T x overlap x C2 + call dgemm( 'N', 'N', n, n, m, 1.d0 & + , S_tmp, size(S_tmp, 1), C2(1,1), size(C2, 1) & + , 0.d0, smat, size(smat, 1) ) + deallocate(S_tmp) + +end + +! --- + +subroutine orthog_functions(m, n, coef, overlap) + + implicit none + + integer, intent(in) :: m, n + double precision, intent(in) :: overlap(m,m) + double precision, intent(inout) :: coef(m,n) + double precision, allocatable :: stmp(:,:) + integer :: j, k + + allocate(stmp(n,n)) + call build_s_matrix(m, n, coef, coef, overlap, stmp) +! print*,'overlap before' +! do j = 1, n +! write(*,'(100(F16.10,X))')stmp(:,j) +! enddo + call impose_orthog_svd_overlap(m, n, coef, overlap) + call build_s_matrix(m, n, coef, coef, overlap, stmp) + do j = 1, n + ! --- + ! TODO: MANU check ici + !coef(1,:m) *= 1.d0/dsqrt(stmp(j,j)) + do k = 1, m + coef(k,j) *= 1.d0/dsqrt(stmp(j,j)) + enddo + ! --- + enddo + call build_s_matrix(m, n, coef, coef, overlap, stmp) + + !print*,'overlap after' + !do j = 1, n + ! write(*,'(100(F16.10,X))')stmp(:,j) + !enddo + + deallocate(stmp) + +end + +! --- + +subroutine print_angles_tc() + + implicit none + integer :: i, j + double precision :: left, right + + print *, ' product of norms, angle between vectors' + do i = 1, mo_num + left = overlap_mo_l(i,i) + right = overlap_mo_r(i,i) +! print*,Fock_matrix_tc_mo_tot(i,i),left*right,angle_left_right(i) + print *, left*right, angle_left_right(i) + enddo + +end + +! --- + +subroutine print_energy_and_mos() + + implicit none + integer :: i + + print *, ' ' + print *, ' TC energy = ', TC_HF_energy + print *, ' TC SCF energy gradient = ', grad_non_hermit + print *, ' Max angle Left/right = ', max_angle_left_right + + if(max_angle_left_right .lt. 45.d0) then + print *, ' Maximum angle BELOW 45 degrees, everthing is OK !' + else if(max_angle_left_right .gt. 45.d0 .and. max_angle_left_right .lt. 75.d0) then + print *, ' Maximum angle between 45 and 75 degrees, this is not the best for TC-CI calculations ...' + else if(max_angle_left_right .gt. 75.d0) then + print *, ' Maximum angle between ABOVE 75 degrees, YOU WILL CERTAINLY FIND TROUBLES IN TC-CI calculations ...' + endif + + print *, ' Diag Fock elem, product of left/right norm, angle left/right ' + do i = 1, mo_num + write(*, '(I3,X,100(F16.10,X))') i, Fock_matrix_tc_mo_tot(i,i), overlap_mo_l(i,i)*overlap_mo_r(i,i), angle_left_right(i) + enddo + +end + +! --- + +subroutine sort_by_tc_fock + implicit none + integer, allocatable :: iorder(:) + double precision, allocatable :: mo_l_tmp(:,:), mo_r_tmp(:,:),fock(:) + allocate(iorder(mo_num),fock(mo_num),mo_l_tmp(ao_num, mo_num),mo_r_tmp(ao_num,mo_num)) + integer :: i + mo_l_tmp = mo_l_coef + mo_r_tmp = mo_r_coef + do i = 1, mo_num + iorder(i) = i + fock(i) = Fock_matrix_tc_mo_tot(i,i) + enddo + call dsort(fock,iorder,mo_num) + do i = 1, mo_num + mo_l_coef(1:ao_num,i) = mo_l_tmp(1:ao_num,iorder(i)) + mo_r_coef(1:ao_num,i) = mo_r_tmp(1:ao_num,iorder(i)) + enddo + touch mo_l_coef mo_r_coef + +end + diff --git a/src/tc_scf/tc_petermann_factor.irp.f b/src/tc_scf/tc_petermann_factor.irp.f new file mode 100644 index 00000000..d3722098 --- /dev/null +++ b/src/tc_scf/tc_petermann_factor.irp.f @@ -0,0 +1,78 @@ + +! --- + +program tc_petermann_factor + + 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 +! 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 main() + +end + +! --- + +subroutine main() + + implicit none + integer :: i, j + double precision :: Pf_diag_av + double precision, allocatable :: Sl(:,:), Sr(:,:), Pf(:,:) + + allocate(Sl(mo_num,mo_num), Sr(mo_num,mo_num), Pf(mo_num,mo_num)) + + call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 & + , mo_l_coef, size(mo_l_coef, 1), mo_l_coef, size(mo_l_coef, 1) & + , 0.d0, Sl, size(Sl, 1) ) + + print *, '' + print *, ' left-orthog matrix:' + do i = 1, mo_num + write(*,'(100(F8.4,X))') Sl(:,i) + enddo + + call dgemm( "T", "N", mo_num, mo_num, ao_num, 1.d0 & + , mo_r_coef, size(mo_r_coef, 1), mo_r_coef, size(mo_r_coef, 1) & + , 0.d0, Sr, size(Sr, 1) ) + + print *, '' + print *, ' right-orthog matrix:' + do i = 1, mo_num + write(*,'(100(F8.4,X))') Sr(:,i) + enddo + + print *, '' + print *, ' Petermann matrix:' + do i = 1, mo_num + do j = 1, mo_num + Pf(j,i) = Sl(j,i) * Sr(j,i) + enddo + write(*,'(100(F8.4,X))') Pf(:,i) + enddo + + Pf_diag_av = 0.d0 + do i = 1, mo_num + Pf_diag_av = Pf_diag_av + Pf(i,i) + enddo + Pf_diag_av = Pf_diag_av / dble(mo_num) + + print *, '' + print *, ' mean of the diagonal Petermann factor = ', Pf_diag_av + + deallocate(Sl, Sr, Pf) + + return +end subroutine + +! --- + diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f new file mode 100644 index 00000000..deaf8d82 --- /dev/null +++ b/src/tc_scf/tc_scf.irp.f @@ -0,0 +1,75 @@ +! --- + +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 + + PROVIDE mu_erf + print *, ' mu = ', mu_erf + PROVIDE j1b_type + print *, ' j1b_type = ', j1b_type + print *, j1b_pen + + !call create_guess() + !call orthonormalize_mos() + + PROVIDE tcscf_algorithm + if(tcscf_algorithm == 'DIIS') then + call rh_tcscf_diis() + elseif(tcscf_algorithm == 'Simple') then + call rh_tcscf_simple() + else + print *, ' not implemented yet', tcscf_algorithm + stop + endif + + call minimize_tc_orb_angles() + call print_energy_and_mos() + +end + +! --- + +subroutine create_guess() + + implicit none + logical :: exists + + PROVIDE ezfio_filename + !call ezfio_has_mo_basis_mo_coef(exists) + exists = .false. + + 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 + elseif (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 + +! --- 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..90719f47 --- /dev/null +++ b/src/tc_scf/tc_scf_dm.irp.f @@ -0,0 +1,37 @@ +! --- + +BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) ] + + implicit none + + if(bi_ortho) then + PROVIDE mo_l_coef mo_r_coef + 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 + PROVIDE mo_l_coef mo_r_coef + 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..611b8b4c --- /dev/null +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -0,0 +1,34 @@ + + BEGIN_PROVIDER [ double precision, TC_HF_energy] +&BEGIN_PROVIDER [ double precision, TC_HF_one_e_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 + + PROVIDE mo_l_coef mo_r_coef + + TC_HF_energy = nuclear_repulsion + TC_HF_one_e_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_e_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_e_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..dde477c4 --- /dev/null +++ b/src/tc_scf/tc_scf_utils.irp.f @@ -0,0 +1,43 @@ + +! --- + +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/tc_scf/test_Ne.sh b/src/tc_scf/test_Ne.sh new file mode 100755 index 00000000..27ea73c2 --- /dev/null +++ b/src/tc_scf/test_Ne.sh @@ -0,0 +1,13 @@ +QP_ROOT=/home/eginer/new_qp2/qp2 +source ${QP_ROOT}/quantum_package.rc + echo Ne > Ne.xyz + echo $QP_ROOT + qp create_ezfio -b cc-pcvdz Ne.xyz + qp run scf + qp set tc_keywords bi_ortho True + qp set ao_two_e_erf_ints mu_erf 0.87 + qp set tc_keywords j1b_pen [1.5] + qp set tc_keywords j1b_type 3 + qp run tc_scf | tee Ne.ezfio.tc_scf.out + grep "TC energy =" Ne.ezfio.tc_scf.out | tail -1 + eref=-128.552134 diff --git a/src/tc_scf/test_int.irp.f b/src/tc_scf/test_int.irp.f new file mode 100644 index 00000000..a14c4126 --- /dev/null +++ b/src/tc_scf/test_int.irp.f @@ -0,0 +1,1003 @@ +program test_ints + + BEGIN_DOC +! TODO : Put the documentation of the program here + END_DOC + + implicit none + + print *, ' starting test_ints ...' + + my_grid_becke = .True. + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 +! my_n_pt_r_grid = 15 ! 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 + + my_extra_grid_becke = .True. + my_n_pt_r_extra_grid = 30 + my_n_pt_a_extra_grid = 50 ! small extra_grid for quick debug + touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid + +!! OK +!call routine_int2_u_grad1u_j1b2 +!! OK +!call routine_v_ij_erf_rk_cst_mu_j1b +!! OK +! call routine_x_v_ij_erf_rk_cst_mu_j1b +!! OK +! call routine_v_ij_u_cst_mu_j1b + +!! OK +!call routine_int2_u2_j1b2 + +!! OK +!call routine_int2_u_grad1u_x_j1b2 + +!! OK +! call routine_int2_grad1u2_grad2u2_j1b2 +! call routine_int2_u_grad1u_j1b2 +! call test_total_grad_lapl +! call test_total_grad_square +! call test_ao_tc_int_chemist +! call test_grid_points_ao +! call test_tc_scf + !call test_int_gauss + + !call test_fock_3e_uhf_ao() + !call test_fock_3e_uhf_mo() + + !call test_tc_grad_and_lapl_ao() + !call test_tc_grad_square_ao() + + call test_two_e_tc_non_hermit_integral() + +end + +! --- + +subroutine test_tc_scf + implicit none + integer :: i +! provide int2_u_grad1u_x_j1b2_test + provide x_v_ij_erf_rk_cst_mu_j1b_test +! provide x_v_ij_erf_rk_cst_mu_j1b_test +! print*,'TC_HF_energy = ',TC_HF_energy +! print*,'grad_non_hermit = ',grad_non_hermit +end + +subroutine test_ao_tc_int_chemist + implicit none + provide ao_tc_int_chemist +! provide ao_tc_int_chemist_test +! provide tc_grad_square_ao_test +! provide tc_grad_and_lapl_ao_test +end + +! --- + +subroutine routine_test_j1b + implicit none + integer :: i,icount,j + icount = 0 + do i = 1, List_all_comb_b3_size + if(dabs(List_all_comb_b3_coef(i)).gt.1.d-10)then + print*,'' + print*,List_all_comb_b3_expo(i),List_all_comb_b3_coef(i) + print*,List_all_comb_b3_cent(1:3,i) + print*,'' + icount += 1 + endif + + enddo + print*,'List_all_comb_b3_coef,icount = ',List_all_comb_b3_size,icount + do i = 1, ao_num + do j = 1, ao_num + do icount = 1, List_comb_thr_b3_size(j,i) + print*,'',j,i + print*,List_comb_thr_b3_expo(icount,j,i),List_comb_thr_b3_coef(icount,j,i) + print*,List_comb_thr_b3_cent(1:3,icount,j,i) + print*,'' + enddo +! enddo + enddo + enddo + print*,'max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size,List_all_comb_b3_size + +end + +subroutine routine_int2_u_grad1u_j1b2 + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_u_grad1u_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_u_grad1u_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_v_ij_erf_rk_cst_mu_j1b + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + + +subroutine routine_x_v_ij_erf_rk_cst_mu_j1b + implicit none + integer :: i,j,ipoint,k,l,m + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + do m = 1, 3 + array(j,i,l,k) += x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += x_v_ij_erf_rk_cst_mu_j1b (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + + + +subroutine routine_v_ij_u_cst_mu_j1b_test + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_int2_grad1u2_grad2u2_j1b2 + implicit none + integer :: i,j,ipoint,k,l + integer :: ii , jj + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + double precision, allocatable :: ints(:,:,:) + allocate(ints(ao_num, ao_num, n_points_final_grid)) +! do ipoint = 1, n_points_final_grid +! do i = 1, ao_num +! do j = 1, ao_num +! read(33,*)ints(j,i,ipoint) +! enddo +! enddo +! enddo + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! !array(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! array_ref(j,i,l,k) += int2_grad1u2_grad2u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += ints(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight +! if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then +! if(dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)).gt.1.d-6)then +! print*,j,i,ipoint +! print*,int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) , int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) - int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint)) +! print*,int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) , int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint), dabs(int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint) - int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)) +! stop +! endif +! endif + enddo + enddo + enddo + enddo + enddo + double precision :: e_ref, e_new + accu_relat = 0.d0 + accu_abs = 0.d0 + e_ref = 0.d0 + e_new = 0.d0 + do ii = 1, elec_alpha_num + do jj = ii, elec_alpha_num + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + e_ref += mo_coef(j,ii) * mo_coef(i,ii) * array_ref(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj) + e_new += mo_coef(j,ii) * mo_coef(i,ii) * array(j,i,l,k) * mo_coef(l,jj) * mo_coef(k,jj) + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib +! if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then +! accu_relat += contrib/dabs(array_ref(j,i,l,k)) +! endif + enddo + enddo + enddo + enddo + + enddo + enddo + print*,'e_ref = ',e_ref + print*,'e_new = ',e_new +! print*,'accu_abs = ',accu_abs/dble(ao_num)**4 +! print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_int2_u2_j1b2 + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += int2_u2_j1b2_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_u2_j1b2(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + + +subroutine routine_int2_u_grad1u_x_j1b2 + implicit none + integer :: i,j,ipoint,k,l,m + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) +! print*,'ao_overlap_abs = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_overlap_abs(i,:) +! enddo +! print*,'center = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_center(2,i,:) +! enddo +! print*,'sigma = ' +! do i = 1, ao_num +! write(*,'(100(F10.5,X))')ao_prod_sigma(i,:) +! enddo + + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + do m = 1, 3 + array(j,i,l,k) += int2_u_grad1u_x_j1b2_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += int2_u_grad1u_x_j1b2 (j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + + +end + +subroutine routine_v_ij_u_cst_mu_j1b + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:) + + allocate(array(ao_num, ao_num, ao_num, ao_num)) + array = 0.d0 + allocate(array_ref(ao_num, ao_num, ao_num, ao_num)) + array_ref = 0.d0 + do ipoint = 1, n_points_final_grid + weight = final_weight_at_r_vector(ipoint) + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + array(j,i,l,k) += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight + enddo + enddo + enddo + enddo + enddo + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k)) + accu_abs += contrib + if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(array_ref(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + +end + +! --- + +subroutine test_fock_3e_uhf_ao() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + double precision, allocatable :: fock_3e_uhf_ao_a_mo(:,:), fock_3e_uhf_ao_b_mo(:,:) + + thr_ih = 1d-7 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_ao_a fock_3e_uhf_ao_b + + ! --- + + allocate(fock_3e_uhf_ao_a_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_a , size(fock_3e_uhf_ao_a , 1) & + , fock_3e_uhf_ao_a_mo, size(fock_3e_uhf_ao_a_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_a_mo(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_a_mo (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' ' + + deallocate(fock_3e_uhf_ao_a_mo) + + ! --- + + allocate(fock_3e_uhf_ao_b_mo(mo_num,mo_num)) + call ao_to_mo_bi_ortho( fock_3e_uhf_ao_b , size(fock_3e_uhf_ao_b , 1) & + , fock_3e_uhf_ao_b_mo, size(fock_3e_uhf_ao_b_mo, 1) ) + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_ao_b_mo(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_ao_b_mo (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' ' + + deallocate(fock_3e_uhf_ao_b_mo) + + ! --- + +end subroutine test_fock_3e_uhf_ao() + +! --- + +subroutine test_fock_3e_uhf_mo() + + implicit none + integer :: i, j + double precision :: diff_tot, diff_ij, thr_ih, norm + + thr_ih = 1d-12 + + PROVIDE fock_a_tot_3e_bi_orth fock_b_tot_3e_bi_orth + PROVIDE fock_3e_uhf_mo_a fock_3e_uhf_mo_b + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_a(j,i) - fock_a_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_a_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_a (j,i) + !stop + endif + + norm += dabs(fock_a_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_a = ', diff_tot / norm + print *, ' norm_a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, mo_num + do j = 1, mo_num + + diff_ij = dabs(fock_3e_uhf_mo_b(j,i) - fock_b_tot_3e_bi_orth(j,i)) + if(diff_ij .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' MANU : ', fock_b_tot_3e_bi_orth(j,i) + print *, ' UHF : ', fock_3e_uhf_mo_b (j,i) + !stop + endif + + norm += dabs(fock_b_tot_3e_bi_orth(j,i)) + diff_tot += diff_ij + enddo + enddo + print *, ' diff on F_b = ', diff_tot/norm + print *, ' norm_b = ', norm + print *, ' ' + + ! --- + +end subroutine test_fock_3e_uhf_mo + +! --- + +subroutine test_total_grad_lapl + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(tc_grad_and_lapl_ao_test(j,i,l,k) - tc_grad_and_lapl_ao(j,i,l,k)) + accu_abs += contrib + if(dabs(tc_grad_and_lapl_ao(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(tc_grad_and_lapl_ao(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + +end + +subroutine test_total_grad_square + implicit none + integer :: i,j,ipoint,k,l + double precision :: weight,accu_relat, accu_abs, contrib + accu_relat = 0.d0 + accu_abs = 0.d0 + do k = 1, ao_num + do l = 1, ao_num + do i = 1, ao_num + do j = 1, ao_num + contrib = dabs(tc_grad_square_ao_test(j,i,l,k) - tc_grad_square_ao(j,i,l,k)) + accu_abs += contrib + if(dabs(tc_grad_square_ao(j,i,l,k)).gt.1.d-10)then + accu_relat += contrib/dabs(tc_grad_square_ao(j,i,l,k)) + endif + enddo + enddo + enddo + enddo + print*,'accu_abs = ',accu_abs/dble(ao_num)**4 + print*,'accu_relat = ',accu_relat/dble(ao_num)**4 + + +end + +subroutine test_grid_points_ao + implicit none + integer :: i,j,ipoint,icount,icount_good, icount_bad,icount_full + double precision :: thr + thr = 1.d-10 +! print*,'max_n_pts_grid_ao_prod = ',max_n_pts_grid_ao_prod +! print*,'n_pts_grid_ao_prod' + do i = 1, ao_num + do j = i, ao_num + icount = 0 + icount_good = 0 + icount_bad = 0 + icount_full = 0 + do ipoint = 1, n_points_final_grid +! if(dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,1)) & +! + dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,2)) & +! + dabs(int2_u_grad1u_x_j1b2_test(j,i,ipoint,3)) ) +! if(dabs(int2_u2_j1b2_test(j,i,ipoint)).gt.thr)then +! icount += 1 +! endif + if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then + icount_full += 1 + endif + if(dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)).gt.thr)then + icount += 1 + if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr*0.1d0)then + icount_good += 1 + else + print*,j,i,ipoint + print*,dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)),dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint))/dabs(v_ij_u_cst_mu_j1b_test(j,i,ipoint)) + icount_bad += 1 + endif + endif +! if(dabs(v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint)).gt.thr)then +! endif + enddo + print*,'' + print*,j,i + print*,icount,icount_full, icount_bad!,n_pts_grid_ao_prod(j,i) + print*,dble(icount)/dble(n_points_final_grid),dble(icount_full)/dble(n_points_final_grid) +! dble(n_pts_grid_ao_prod(j,i))/dble(n_points_final_grid) +! if(icount.gt.n_pts_grid_ao_prod(j,i))then +! print*,'pb !!' +! endif + enddo + enddo +end + +subroutine test_int_gauss + implicit none + integer :: i,j + print*,'center' + do i = 1, ao_num + do j = i, ao_num + print*,j,i + print*,ao_prod_sigma(j,i),ao_overlap_abs_grid(j,i) + print*,ao_prod_center(1:3,j,i) + enddo + enddo + print*,'' + double precision :: weight, r(3),integral_1,pi,center(3),f_r,alpha,distance,integral_2 + center = 0.d0 + pi = dacos(-1.d0) + integral_1 = 0.d0 + integral_2 = 0.d0 + alpha = 0.75d0 + do i = 1, n_points_final_grid + ! you get x, y and z of the ith grid point + r(1) = final_grid_points(1,i) + r(2) = final_grid_points(2,i) + r(3) = final_grid_points(3,i) + weight = final_weight_at_r_vector(i) + distance = dsqrt( (r(1) - center(1))**2 + (r(2) - center(2))**2 + (r(3) - center(3))**2 ) + f_r = dexp(-alpha * distance*distance) + ! you add the contribution of the grid point to the integral + integral_1 += f_r * weight + integral_2 += f_r * distance * weight + enddo + print*,'integral_1 =',integral_1 + print*,'(pi/alpha)**1.5 =',(pi / alpha)**1.5 + print*,'integral_2 =',integral_2 + print*,'(pi/alpha)**1.5 =',2.d0*pi / (alpha)**2 + + +end + +! --- + +subroutine test_tc_grad_and_lapl_ao() + + implicit none + integer :: i, j, k, l + double precision :: diff_tot, diff, thr_ih, norm + + thr_ih = 1d-10 + + PROVIDE tc_grad_and_lapl_ao tc_grad_and_lapl_ao_loop + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + + diff = dabs(tc_grad_and_lapl_ao_loop(l,k,j,i) - tc_grad_and_lapl_ao(l,k,j,i)) + if(diff .gt. thr_ih) then + print *, ' difference on ', l, k, j, i + print *, ' loops : ', tc_grad_and_lapl_ao_loop(l,k,j,i) + print *, ' lapack: ', tc_grad_and_lapl_ao (l,k,j,i) + !stop + endif + + norm += dabs(tc_grad_and_lapl_ao_loop(l,k,j,i)) + diff_tot += diff + enddo + enddo + enddo + enddo + + print *, ' diff tot = ', diff_tot / norm + print *, ' norm = ', norm + print *, ' ' + + return + +end + +! --- + +subroutine test_tc_grad_square_ao() + + implicit none + integer :: i, j, k, l + double precision :: diff_tot, diff, thr_ih, norm + + thr_ih = 1d-10 + + PROVIDE tc_grad_square_ao tc_grad_square_ao_loop + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + do k = 1, ao_num + do l = 1, ao_num + + diff = dabs(tc_grad_square_ao_loop(l,k,j,i) - tc_grad_square_ao(l,k,j,i)) + if(diff .gt. thr_ih) then + print *, ' difference on ', l, k, j, i + print *, ' loops : ', tc_grad_square_ao_loop(l,k,j,i) + print *, ' lapack: ', tc_grad_square_ao (l,k,j,i) + !stop + endif + + norm += dabs(tc_grad_square_ao_loop(l,k,j,i)) + diff_tot += diff + enddo + enddo + enddo + enddo + + print *, ' diff tot = ', diff_tot / norm + print *, ' norm = ', norm + print *, ' ' + + return + +end + +! --- + +subroutine test_two_e_tc_non_hermit_integral() + + implicit none + integer :: i, j + double precision :: diff_tot, diff, thr_ih, norm + + thr_ih = 1d-10 + + PROVIDE two_e_tc_non_hermit_integral_beta two_e_tc_non_hermit_integral_alpha + PROVIDE two_e_tc_non_hermit_integral_seq_beta two_e_tc_non_hermit_integral_seq_alpha + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + + diff = dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i) - two_e_tc_non_hermit_integral_alpha(j,i)) + if(diff .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' seq : ', two_e_tc_non_hermit_integral_seq_alpha(j,i) + print *, ' // : ', two_e_tc_non_hermit_integral_alpha (j,i) + !stop + endif + + norm += dabs(two_e_tc_non_hermit_integral_seq_alpha(j,i)) + diff_tot += diff + enddo + enddo + + print *, ' diff tot a = ', diff_tot / norm + print *, ' norm a = ', norm + print *, ' ' + + ! --- + + norm = 0.d0 + diff_tot = 0.d0 + do i = 1, ao_num + do j = 1, ao_num + + diff = dabs(two_e_tc_non_hermit_integral_seq_beta(j,i) - two_e_tc_non_hermit_integral_beta(j,i)) + if(diff .gt. thr_ih) then + print *, ' difference on ', j, i + print *, ' seq : ', two_e_tc_non_hermit_integral_seq_beta(j,i) + print *, ' // : ', two_e_tc_non_hermit_integral_beta (j,i) + !stop + endif + + norm += dabs(two_e_tc_non_hermit_integral_seq_beta(j,i)) + diff_tot += diff + enddo + enddo + + print *, ' diff tot b = ', diff_tot / norm + print *, ' norm b = ', norm + print *, ' ' + + ! --- + + return + +end + +! --- + +>>>>>>> 92a4e33f8a21717cab0c0e4f8412ed6903afb04a diff --git a/src/tc_scf/three_e_energy_bi_ortho.irp.f b/src/tc_scf/three_e_energy_bi_ortho.irp.f new file mode 100644 index 00000000..64212da8 --- /dev/null +++ b/src/tc_scf/three_e_energy_bi_ortho.irp.f @@ -0,0 +1,174 @@ + +subroutine contrib_3e_diag_sss(i,j,k,integral) + implicit none + integer, intent(in) :: i,j,k + BEGIN_DOC + ! returns the pure same spin contribution to diagonal matrix element of 3e term + END_DOC + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j > + call give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int) ! < i k j | j i k > + call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12 + integral += - exch_13_int - exch_23_int - exch_12_int + integral = -integral +end + +subroutine contrib_3e_diag_soo(i,j,k,integral) + implicit none + integer, intent(in) :: i,j,k + BEGIN_DOC + ! returns the pure same spin contribution to diagonal matrix element of 3e term + END_DOC + double precision, intent(out) :: integral + double precision :: direct_int, exch_23_int + call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int) ! < i k j | i k j > + call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)! < i k j | i j k > : E_23 + integral = direct_int - exch_23_int + integral = -integral +end + + +subroutine give_aaa_contrib_bis(integral_aaa) + implicit none + double precision, intent(out) :: integral_aaa + double precision :: integral + integer :: i,j,k + integral_aaa = 0.d0 + do i = 1, elec_alpha_num + do j = i+1, elec_alpha_num + do k = j+1, elec_alpha_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_aaa += integral + enddo + enddo + enddo + +end + +subroutine give_aaa_contrib(integral_aaa) + implicit none + double precision, intent(out) :: integral_aaa + double precision :: integral + integer :: i,j,k + integral_aaa = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_aaa += integral + enddo + enddo + enddo + integral_aaa *= 1.d0/6.d0 +end + + +subroutine give_aab_contrib(integral_aab) + implicit none + double precision, intent(out) :: integral_aab + double precision :: integral + integer :: i,j,k + integral_aab = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_aab += integral + enddo + enddo + enddo + integral_aab *= 0.5d0 +end + + +subroutine give_aab_contrib_bis(integral_aab) + implicit none + double precision, intent(out) :: integral_aab + double precision :: integral + integer :: i,j,k + integral_aab = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_alpha_num + do k = j+1, elec_alpha_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_aab += integral + enddo + enddo + enddo +end + + +subroutine give_abb_contrib(integral_abb) + implicit none + double precision, intent(out) :: integral_abb + double precision :: integral + integer :: i,j,k + integral_abb = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_abb += integral + enddo + enddo + enddo + integral_abb *= 0.5d0 +end + +subroutine give_abb_contrib_bis(integral_abb) + implicit none + double precision, intent(out) :: integral_abb + double precision :: integral + integer :: i,j,k + integral_abb = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_beta_num + do k = j+1, elec_beta_num + call contrib_3e_diag_soo(i,j,k,integral) + integral_abb += integral + enddo + enddo + enddo +end + +subroutine give_bbb_contrib_bis(integral_bbb) + implicit none + double precision, intent(out) :: integral_bbb + double precision :: integral + integer :: i,j,k + integral_bbb = 0.d0 + do i = 1, elec_beta_num + do j = i+1, elec_beta_num + do k = j+1, elec_beta_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_bbb += integral + enddo + enddo + enddo + +end + +subroutine give_bbb_contrib(integral_bbb) + implicit none + double precision, intent(out) :: integral_bbb + double precision :: integral + integer :: i,j,k + integral_bbb = 0.d0 + do i = 1, elec_beta_num + do j = 1, elec_beta_num + do k = 1, elec_beta_num + call contrib_3e_diag_sss(i,j,k,integral) + integral_bbb += integral + enddo + enddo + enddo + integral_bbb *= 1.d0/6.d0 +end + + diff --git a/src/utils/block_diag_degen.irp.f b/src/utils/block_diag_degen.irp.f new file mode 100644 index 00000000..188bfa58 --- /dev/null +++ b/src/utils/block_diag_degen.irp.f @@ -0,0 +1,218 @@ + +subroutine diag_mat_per_fock_degen(fock_diag, mat_ref, n, thr_d, thr_nd, thr_deg, leigvec, reigvec, eigval) + + + BEGIN_DOC + ! + ! subroutine that diagonalizes a matrix mat_ref BY BLOCK + ! + ! the blocks are defined by the elements having the SAME DEGENERACIES in the entries "fock_diag" + ! + ! examples : all elements having degeneracy 1 in fock_diag (i.e. not being degenerated) will be treated together + ! + ! : all elements having degeneracy 2 in fock_diag (i.e. two elements are equal) will be treated together + ! + ! : all elements having degeneracy 3 in fock_diag (i.e. two elements are equal) will be treated together + ! + ! etc... the advantage is to guarentee no spurious mixing because of numerical problems. + ! + END_DOC + + implicit none + integer, intent(in) :: n + double precision, intent(in) :: fock_diag(n), mat_ref(n,n), thr_d, thr_nd, thr_deg + double precision, intent(out) :: leigvec(n,n), reigvec(n,n), eigval(n) + + integer :: n_degen_list, n_degen,size_mat, i, j, k, icount, m, index_degen + integer :: ii, jj, i_good, j_good, n_real + integer :: icount_eigval + logical, allocatable :: is_ok(:) + integer, allocatable :: list_degen(:,:), list_same_degen(:) + integer, allocatable :: iorder(:), list_degen_sorted(:) + double precision, allocatable :: leigvec_unsrtd(:,:), reigvec_unsrtd(:,:), eigval_unsrtd(:) + double precision, allocatable :: mat_tmp(:,:), eigval_tmp(:), leigvec_tmp(:,:), reigvec_tmp(:,:) + + allocate(leigvec_unsrtd(n,n), reigvec_unsrtd(n,n), eigval_unsrtd(n)) + leigvec_unsrtd = 0.d0 + reigvec_unsrtd = 0.d0 + eigval_unsrtd = 0.d0 + + ! obtain degeneracies + allocate(list_degen(n,0:n)) + call give_degen_full_list(fock_diag, n, thr_deg, list_degen, n_degen_list) + + allocate(iorder(n_degen_list), list_degen_sorted(n_degen_list)) + do i = 1, n_degen_list + n_degen = list_degen(i,0) + list_degen_sorted(i) = n_degen + iorder(i) = i + enddo + + ! sort by number of degeneracies + call isort(list_degen_sorted, iorder, n_degen_list) + + allocate(is_ok(n_degen_list)) + is_ok = .True. + icount_eigval = 0 + + ! loop over degeneracies + do i = 1, n_degen_list + if(.not.is_ok(i)) cycle + + is_ok(i) = .False. + n_degen = list_degen_sorted(i) + + print *, ' diagonalizing for n_degen = ', n_degen + + k = 1 + + ! group all the entries having the same degeneracies +!! do while (list_degen_sorted(i+k)==n_degen) + do m = i+1, n_degen_list + if(list_degen_sorted(m)==n_degen) then + is_ok(i+k) = .False. + k += 1 + endif + enddo + + print *, ' number of identical degeneracies = ', k + size_mat = k*n_degen + print *, ' size_mat = ', size_mat + allocate(mat_tmp(size_mat,size_mat), list_same_degen(size_mat)) + allocate(eigval_tmp(size_mat), leigvec_tmp(size_mat,size_mat), reigvec_tmp(size_mat,size_mat)) + ! group all the elements sharing the same degeneracy + icount = 0 + do j = 1, k ! jth set of degeneracy + index_degen = iorder(i+j-1) + do m = 1, n_degen + icount += 1 + list_same_degen(icount) = list_degen(index_degen,m) + enddo + enddo + + print *, ' list of elements ' + do icount = 1, size_mat + print *, icount, list_same_degen(icount) + enddo + + ! you copy subset of matrix elements having all the same degeneracy in mat_tmp + do ii = 1, size_mat + i_good = list_same_degen(ii) + do jj = 1, size_mat + j_good = list_same_degen(jj) + mat_tmp(jj,ii) = mat_ref(j_good,i_good) + enddo + enddo + + call non_hrmt_bieig( size_mat, mat_tmp, thr_d, thr_nd & + , leigvec_tmp, reigvec_tmp & + , n_real, eigval_tmp ) + + do ii = 1, size_mat + icount_eigval += 1 + eigval_unsrtd(icount_eigval) = eigval_tmp(ii) ! copy eigenvalues + do jj = 1, size_mat ! copy the eigenvectors + j_good = list_same_degen(jj) + leigvec_unsrtd(j_good,icount_eigval) = leigvec_tmp(jj,ii) + reigvec_unsrtd(j_good,icount_eigval) = reigvec_tmp(jj,ii) + enddo + enddo + + deallocate(mat_tmp, list_same_degen) + deallocate(eigval_tmp, leigvec_tmp, reigvec_tmp) + enddo + + if(icount_eigval .ne. n) then + print *, ' pb !! (icount_eigval.ne.n)' + print *, ' icount_eigval,n', icount_eigval, n + stop + endif + + deallocate(iorder) + allocate(iorder(n)) + do i = 1, n + iorder(i) = i + enddo + call dsort(eigval_unsrtd, iorder, n) + + do i = 1, n + print*,'sorted eigenvalues ' + i_good = iorder(i) + eigval(i) = eigval_unsrtd(i) + print*,'i,eigval(i) = ',i,eigval(i) + do j = 1, n + leigvec(j,i) = leigvec_unsrtd(j,i_good) + reigvec(j,i) = reigvec_unsrtd(j,i_good) + enddo + enddo + + deallocate(leigvec_unsrtd, reigvec_unsrtd, eigval_unsrtd) + deallocate(list_degen) + deallocate(iorder, list_degen_sorted) + deallocate(is_ok) + +end + +! --- + +subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list) + + BEGIN_DOC + ! you enter with an array A(n) and spits out all the elements degenerated up to thr + ! + ! the elements of A(n) DON'T HAVE TO BE SORTED IN THE ENTRANCE: TOTALLY GENERAL + ! + ! list_degen(i,0) = number of degenerate entries + ! + ! list_degen(i,1) = index of the first degenerate entry + ! + ! list_degen(i,2:list_degen(i,0)) = list of all other dengenerate entries + ! + ! if list_degen(i,0) == 1 it means that there is no degeneracy for that element + END_DOC + + implicit none + + double precision, intent(in) :: A(n) + double precision, intent(in) :: thr + integer, intent(in) :: n + integer, intent(out) :: list_degen(n,0:n), n_degen_list + integer :: i, j, icount, icheck + logical, allocatable :: is_ok(:) + + + allocate(is_ok(n)) + n_degen_list = 0 + is_ok = .True. + do i = 1, n + if(.not.is_ok(i)) cycle + n_degen_list +=1 + is_ok(i) = .False. + list_degen(n_degen_list,1) = i + icount = 1 + do j = i+1, n + if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then + is_ok(j) = .False. + icount += 1 + list_degen(n_degen_list,icount) = j + endif + enddo + + list_degen(n_degen_list,0) = icount + enddo + + icheck = 0 + do i = 1, n_degen_list + icheck += list_degen(i,0) + enddo + + if(icheck.ne.n)then + print *, ' pb ! :: icheck.ne.n' + print *, icheck, n + stop + endif + +end + +! --- + diff --git a/src/utils/loc.f b/src/utils/loc.f new file mode 100644 index 00000000..02693281 --- /dev/null +++ b/src/utils/loc.f @@ -0,0 +1,327 @@ +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) + double precision, intent(inout) :: s(n,n) + double precision, intent(out) :: t(n,n),w(n,n) + 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 i=1,m + write (6,145) i + write (6,150) (s(i,j),j=1,n) + end do + 8 mm=m-1 + if (m.lt.n) mm=m + iter=0 + do j=1,n + do i=1,n + t(i,j)=0.d0 + end do + do i=1,m + w(i,j)=s(i,j) + enddo + t(j,j)=1.d0 + enddo + sum=0.d0 + do i=1,m + sum=sum+s(i,i)*s(i,i) + end do + 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 i=1,n + sum=sum+s(i,j)*s(i,j) + enddo + if (sum.gt.small) goto 28 + do 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 + end do + 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 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 + s(k,j)=q + enddo + do 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 + enddo + 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 j=1,n + if (s(j,j).gt.0.d0) cycle + do i=1,m + s(i,j)=-s(i,j) + enddo + do i=1,n + t(i,j)=-t(i,j) + enddo + enddo + sum=0.d0 + do i=1,m + sum=sum+s(i,i)*s(i,i) + enddo + sum=sum/m + do i=1,m + do j=1,n + sw=s(i,j) + s(i,j)=w(i,j) + w(i,j)=sw + enddo + enddo + if (.not.zprt) return + write (6,12) sum + write (6,130) + 130 format (//5x,'transformation matrix') + do i=1,n + write (6,145) i + write (6,150) (t(i,j),j=1,n) + enddo + 145 format (i8) + 150 format (2x,10f12.8) + write (6,160) + 160 format (//5x,'new overlap matrix'/) + do i=1,m + write (6,145) i + write (6,150) (w(i,j),j=1,n) + enddo + 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 j=1,n + do i=1,n + t(i,j)=0.d0 + enddo + do i=1,m + w(i,j)=s(i,j) + enddo + t(j,j)=1.d0 + enddo + sum=0.d0 + do i=1,m + sum=sum+s(i,i)*s(i,i) + enddo + 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 i=1,n + sum=sum+s(i,j)*s(i,j) + enddo + if (sum.gt.small) goto 28 + do 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 + end do + 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 i=1,mm + ip=i+1 + do 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 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 + s(k,j)=q + enddo + do 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 + enddo + 45 d=dabs(sine) + if (d.le.amax) goto 50 + imax=i + jmax=j + amax=d + dmax=delta + 50 continue + end do + end do + 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 j=1,n + if (s(j,j).gt.0.d0) cycle + do i=1,m + s(i,j)=-s(i,j) + enddo + do i=1,n + t(i,j)=-t(i,j) + enddo + enddo + sum=0.d0 + do i=1,m + sum=sum+s(i,i)*s(i,i) + enddo + sum=sum/m + do i=1,m + do j=1,n + sw=s(i,j) + s(i,j)=w(i,j) + w(i,j)=sw + enddo + enddo + return + end +