LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-11-12 16:33:37 +01:00
Code Releases Activity

v0 of DIIS & level shift

Browse Source
This commit is contained in:
AbdAmmar 2022-12-07 11:26:34 +01:00
parent f6f1ca5fa4
commit f2c3c72978
11 changed files with 724 additions and 47 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
src/bi_ortho_mos/bi_density.irp.f
View File

@ -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

47
src/bi_ortho_mos/mos_rl.irp.f
View File

@ -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
! ---

24
src/tc_keywords/EZFIO.cfg
View File

@ -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

39
src/tc_scf/diago_bi_ort_tcfock.irp.f
View File

@ -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) ]

181
src/tc_scf/diis_tcscf.irp.f Normal file
View File

@ -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
! ---
~

40
src/tc_scf/fock_tc.irp.f
View File

@ -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
! ---

367
src/tc_scf/rh_tcscf.irp.f Normal file
View File

@ -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
! ---

24
src/tc_scf/tc_scf.irp.f
View File

@ -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
! ---

34
src/tc_scf/tc_scf_dm.irp.f
View File

@ -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

12
src/tc_scf/tc_scf_energy.irp.f
View File

@ -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

1
src/tc_scf/tc_scf_utils.irp.f
View File

@ -40,3 +40,4 @@ subroutine LTxSxR(n, m, L, S, R, C)
end subroutine LTxR
! ---