diff --git a/src/bi_ortho_mos/bi_density.irp.f b/src/bi_ortho_mos/bi_density.irp.f index 947be870..0de8ce69 100644 --- a/src/bi_ortho_mos/bi_density.irp.f +++ b/src/bi_ortho_mos/bi_density.irp.f @@ -22,7 +22,7 @@ BEGIN_PROVIDER [ double precision, TCSCF_bi_ort_dm_ao_beta, (ao_num, ao_num) ] ! ! This is the equivalent of the beta density of the HF Slater determinant, but with a couple of bi-orthonormal Slater determinant |Chi_0> and |Phi_0> END_DOC - call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & + call dgemm( 'N', 'T', ao_num, ao_num, elec_beta_num, 1.d0 & , mo_l_coef, size(mo_l_coef, 1), mo_r_coef, size(mo_r_coef, 1) & , 0.d0, TCSCF_bi_ort_dm_ao_beta, size(TCSCF_bi_ort_dm_ao_beta, 1) ) END_PROVIDER diff --git a/src/bi_ortho_mos/mos_rl.irp.f b/src/bi_ortho_mos/mos_rl.irp.f index 034a436e..9e3ed358 100644 --- a/src/bi_ortho_mos/mos_rl.irp.f +++ b/src/bi_ortho_mos/mos_rl.irp.f @@ -37,6 +37,52 @@ end subroutine ao_to_mo_bi_ortho ! --- +subroutine mo_to_ao_bi_ortho(A_mo, LDA_mo, A_ao, LDA_ao) + + BEGIN_DOC + ! + ! mo_l_coef.T x A_ao x mo_r_coef = A_mo + ! mo_l_coef.T x ao_overlap x mo_r_coef = I + ! + ! ==> A_ao = (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T + ! + END_DOC + + implicit none + integer, intent(in) :: LDA_ao, LDA_mo + double precision, intent(in) :: A_mo(LDA_mo,mo_num) + double precision, intent(out) :: A_ao(LDA_ao,ao_num) + double precision, allocatable :: tmp_1(:,:), tmp_2(:,:) + + ! ao_overlap x mo_r_coef + allocate( tmp_1(ao_num,mo_num) ) + call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 & + , ao_overlap, size(ao_overlap, 1), mo_r_coef, size(mo_r_coef, 1) & + , 0.d0, tmp_1, size(tmp_1, 1) ) + + ! (ao_overlap x mo_r_coef) x A_mo + allocate( tmp_1(ao_num,mo_num) ) + call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 & + , tmp_1, size(tmp_1, 1), A_mo, LDA_mo & + , 0.d0, tmp_2, size(tmp_2, 1) ) + + ! ao_overlap x mo_l_coef + tmp_1 = 0.d0 + call dgemm( 'N', 'N', ao_num, mo_num, ao_num, 1.d0 & + , ao_overlap, size(ao_overlap, 1), mo_l_coef, size(mo_l_coef, 1) & + , 0.d0, tmp_1, size(tmp_1, 1) ) + + ! (ao_overlap x mo_r_coef) x A_mo x (ao_overlap x mo_l_coef).T + call dgemm( 'N', 'T', ao_num, mo_num, mo_num, 1.d0 & + , tmp_2, size(tmp_2, 1), tmp_1, size(tmp_1, 1) & + , 0.d0, A_ao, LDA_ao ) + + deallocate(tmp_1, tmp_2) + +end subroutine mo_to_ao_bi_ortho + +! --- + BEGIN_PROVIDER [ double precision, mo_r_coef, (ao_num, mo_num) ] BEGIN_DOC @@ -175,3 +221,4 @@ END_PROVIDER ! --- + diff --git a/src/tc_keywords/EZFIO.cfg b/src/tc_keywords/EZFIO.cfg index 1afa26e9..8db73f9a 100644 --- a/src/tc_keywords/EZFIO.cfg +++ b/src/tc_keywords/EZFIO.cfg @@ -136,3 +136,27 @@ doc: nb of Gaussians used to fit Jastrow fcts interface: ezfio,provider,ocaml default: 6 +[max_dim_diis_tcscf] +type: integer +doc: Maximum size of the DIIS extrapolation procedure +interface: ezfio,provider,ocaml +default: 15 + +[threshold_diis_tcscf] +type: Threshold +doc: Threshold on the convergence of the DIIS error vector during a TCSCF calculation. If 0. is chosen, the square root of thresh_tcscf will be used. +interface: ezfio,provider,ocaml +default: 0. + +[level_shift_tcscf] +type: Positive_float +doc: Energy shift on the virtual MOs to improve TCSCF convergence +interface: ezfio,provider,ocaml +default: 0. + +[tcscf_algorithm] +type: character*(32) +doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS] +interface: ezfio,provider,ocaml +default: DIIS + diff --git a/src/tc_scf/diago_bi_ort_tcfock.irp.f b/src/tc_scf/diago_bi_ort_tcfock.irp.f index 856b7382..29ca0efe 100644 --- a/src/tc_scf/diago_bi_ort_tcfock.irp.f +++ b/src/tc_scf/diago_bi_ort_tcfock.irp.f @@ -1,3 +1,5 @@ +! --- + 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)] @@ -9,32 +11,50 @@ implicit none integer :: n_real_tc - integer :: i, k, l + integer :: i, j, k, l double precision :: accu_d, accu_nd, accu_tmp double precision :: thr_d, thr_nd double precision :: norm double precision, allocatable :: eigval_right_tmp(:) + double precision, allocatable :: F_tmp(:,:) thr_d = 1d-6 thr_nd = 1d-6 - allocate( eigval_right_tmp(mo_num) ) + allocate( eigval_right_tmp(mo_num), F_tmp(mo_num,mo_num) ) PROVIDE Fock_matrix_tc_mo_tot - call non_hrmt_bieig( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd & - , fock_tc_leigvec_mo, fock_tc_reigvec_mo & + 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, thr_d, thr_nd & + , fock_tc_leigvec_mo, fock_tc_reigvec_mo & , n_real_tc, eigval_right_tmp ) + !if(max_ov_tc_scf)then - ! call non_hrmt_fock_mat( mo_num, Fock_matrix_tc_mo_tot, thr_d, thr_nd & + ! call non_hrmt_fock_mat( mo_num, F_tmp, thr_d, thr_nd & ! , fock_tc_leigvec_mo, fock_tc_reigvec_mo & ! , n_real_tc, eigval_right_tmp ) !else - ! call non_hrmt_diag_split_degen_bi_orthog( mo_num, Fock_matrix_tc_mo_tot & + ! 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 @@ -42,9 +62,12 @@ eigval_fock_tc_mo = eigval_right_tmp ! print*,'Eigenvalues of Fock_matrix_tc_mo_tot' -! do i = 1, mo_num +! 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 @@ -102,6 +125,8 @@ 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) ] diff --git a/src/tc_scf/diis_tcscf.irp.f b/src/tc_scf/diis_tcscf.irp.f new file mode 100644 index 00000000..cf339175 --- /dev/null +++ b/src/tc_scf/diis_tcscf.irp.f @@ -0,0 +1,181 @@ +! --- + +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 + + 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 + + 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 + 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 P x S - S x P 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) ) + +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_tc.irp.f b/src/tc_scf/fock_tc.irp.f index 6b1c1d77..b9611836 100644 --- a/src/tc_scf/fock_tc.irp.f +++ b/src/tc_scf/fock_tc.irp.f @@ -141,27 +141,49 @@ 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 + + implicit none integer :: i, k - grad_non_hermit_left = 0.d0 + + 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+= dabs(Fock_matrix_tc_mo_tot(k,i)) - 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)) enddo enddo + do i = 1, elec_beta_num ! doc --> virt do k = elec_alpha_num+1, mo_num - 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 += dabs(Fock_matrix_tc_mo_tot(k,i)) + grad_non_hermit_right += dabs(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+= dabs(Fock_matrix_tc_mo_tot(k,i)) - 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)) enddo enddo - grad_non_hermit = grad_non_hermit_left + 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/rh_tcscf.irp.f b/src/tc_scf/rh_tcscf.irp.f new file mode 100644 index 00000000..dc7a34fc --- /dev/null +++ b/src/tc_scf/rh_tcscf.irp.f @@ -0,0 +1,367 @@ +! --- + +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_TCSCF_save + 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(delta_gradie_TCSCF) > 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 + exit + endif + + ! TODO + !if(frozen_orb_scf) then + ! call initialize_mo_coef_begin_iteration + !endif + + ! current size of the DIIS space + 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 + + ! Compute the extrapolated Fock matrix + 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 + + 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 + + ! --- + + ! TODO + !if(frozen_orb_scf) then + ! call reorder_core_orb + ! call initialize_mo_coef_begin_iteration + !endif + + ! calculate error vectors + max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo)) + + 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 + gradie_TCSCF = grad_non_hermit + delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous + delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous + + if((TCSCF_algorithm == 'DIIS') .and. (delta_gradie_TCSCF > 0.d0)) then + Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_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 + endif + + ! --- + + level_shift_TCSCF_save = level_shift_TCSCF + + mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num) + mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num) + + do while(delta_gradie_TCSCF > 0.d0) + + mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num) + mo_l_coef(1:ao_num,1:mo_num) = mo_l_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_r_coef mo_l_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) + + !if(frozen_orb_scf) then + ! call reorder_core_orb + ! call initialize_mo_coef_begin_iteration + !endif + TOUCH mo_l_coef mo_r_coef + + 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 + gradie_TCSCF = grad_non_hermit + delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous + delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous + + if(level_shift_TCSCF - level_shift_TCSCF_save > 40.d0) then + level_shift_TCSCF = level_shift_TCSCF_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 + + energy_TCSCF_previous = energy_TCSCF + 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 + + 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| = ', delta_energy_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 + + if(delta_gradie_TCSCF < 0.d0) then + call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef) + call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) + call ezfio_set_tc_scf_bitc_energy(energy_TCSCF) + endif + + if(qp_stop()) exit + enddo + + ! --- + + !if(iteration_TCSCF < n_it_TCSCF_max) then + ! mo_label = 'Canonical' + !endif + + !if(.not.frozen_orb_scf) then + ! call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo, size(Fock_matrix_mo,1), size(Fock_matrix_mo, 2), mo_label, 1, .true.) + ! call restore_symmetry(ao_num, mo_num, mo_coef, size(mo_coef, 1), 1.d-10) + ! call orthonormalize_mos + ! call save_mos + !endif + !call write_double(6, energy_TCSCF, 'TCSCF energy') + + call write_time(6) + +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/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 48cbbdc0..2b751e50 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -18,11 +18,16 @@ program tc_scf !call create_guess !call orthonormalize_mos - call routine_scf() + PROVIDE tcscf_algorithm + if(tcscf_algorithm == 'DIIS') then + call rh_tcscf() + else + call simple_tcscf() + endif + call minimize_tc_orb_angles() call print_energy_and_mos() - end ! --- @@ -64,7 +69,7 @@ end subroutine create_guess ! --- -subroutine routine_scf() +subroutine simple_tcscf() implicit none integer :: i, j, it @@ -79,9 +84,9 @@ subroutine routine_scf() !print*,'grad_hermit = ', grad_hermit print*,'***' print*,'TC HF total energy = ', TC_HF_energy - print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy + print*,'TC HF 1 e energy = ', TC_HF_one_e_energy print*,'TC HF 2 e energy = ', TC_HF_two_e_energy - if(three_body_h_tc)then + if(three_body_h_tc) then print*,'TC HF 3 body = ', diag_three_elem_hf endif print*,'***' @@ -99,7 +104,6 @@ subroutine routine_scf() call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef) TOUCH mo_l_coef mo_r_coef - else print*,'grad_hermit = ',grad_hermit @@ -122,7 +126,7 @@ subroutine routine_scf() print*,'iteration = ', it print*,'***' print*,'TC HF total energy = ', TC_HF_energy - print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy + print*,'TC HF 1 e energy = ', TC_HF_one_e_energy print*,'TC HF 2 non hermit = ', TC_HF_two_e_energy if(three_body_h_tc)then print*,'TC HF 3 body = ', diag_three_elem_hf @@ -161,7 +165,7 @@ subroutine routine_scf() print*,'iteration = ', it print*,'***' print*,'TC HF total energy = ', TC_HF_energy - print*,'TC HF 1 e energy = ', TC_HF_one_electron_energy + print*,'TC HF 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*,'***' @@ -174,11 +178,11 @@ subroutine routine_scf() endif print*,'Energy converged !' - call print_energy_and_mos + call print_energy_and_mos() deallocate(rho_old, rho_new) -end subroutine routine_scf +end subroutine simple_tcscf ! --- diff --git a/src/tc_scf/tc_scf_dm.irp.f b/src/tc_scf/tc_scf_dm.irp.f index f6ae3e1f..1f054a30 100644 --- a/src/tc_scf/tc_scf_dm.irp.f +++ b/src/tc_scf/tc_scf_dm.irp.f @@ -1,25 +1,31 @@ +! --- + BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_beta, (ao_num, ao_num) ] - implicit none - if(bi_ortho)then - TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta - else - TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta - endif + implicit none + if(bi_ortho) then + TCSCF_density_matrix_ao_beta = TCSCF_bi_ort_dm_ao_beta + else + TCSCF_density_matrix_ao_beta = SCF_density_matrix_ao_beta + endif END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_alpha, (ao_num, ao_num) ] - implicit none - if(bi_ortho)then - TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha - else - TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha - endif + implicit none + if(bi_ortho)then + TCSCF_density_matrix_ao_alpha = TCSCF_bi_ort_dm_ao_alpha + else + TCSCF_density_matrix_ao_alpha = SCF_density_matrix_ao_alpha + endif END_PROVIDER +! --- + BEGIN_PROVIDER [ double precision, TCSCF_density_matrix_ao_tot, (ao_num, ao_num) ] - implicit none - TCSCF_density_matrix_ao_tot = TCSCF_density_matrix_ao_beta + TCSCF_density_matrix_ao_alpha + 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 index aa2a16ff..c60ce761 100644 --- a/src/tc_scf/tc_scf_energy.irp.f +++ b/src/tc_scf/tc_scf_energy.irp.f @@ -1,6 +1,6 @@ BEGIN_PROVIDER [ double precision, TC_HF_energy] -&BEGIN_PROVIDER [ double precision, TC_HF_one_electron_energy] +&BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy] &BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy] BEGIN_DOC @@ -11,19 +11,19 @@ integer :: i, j TC_HF_energy = nuclear_repulsion - TC_HF_one_electron_energy = 0.d0 + 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_electron_energy += ao_one_e_integrals_tc_tot(i,j) & - * (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (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_electron_energy + TC_HF_two_e_energy + 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 index 09a4a1b9..dde477c4 100644 --- a/src/tc_scf/tc_scf_utils.irp.f +++ b/src/tc_scf/tc_scf_utils.irp.f @@ -40,3 +40,4 @@ subroutine LTxSxR(n, m, L, S, R, C) end subroutine LTxR ! --- +