PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-23 04:43:45 +01:00
Code Releases Activity

testing for real kpts; not clean

Browse Source
This commit is contained in:
Kevin Gasperich 2020-07-13 18:24:37 -05:00
parent 75dbda613a
commit 4d9299ad7c
13 changed files with 669 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
src/ao_one_e_ints/ao_overlap.irp.f
View File

@ -122,6 +122,22 @@ BEGIN_PROVIDER [ complex*16, ao_overlap_kpts, (ao_num_per_kpt, ao_num_per_kpt, k
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_kpts_real, (ao_num_per_kpt, ao_num_per_kpt, kpt_num) ]
implicit none
BEGIN_DOC
! Overlap for complex AOs
END_DOC
integer :: i,j,k
do k=1,kpt_num
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
ao_overlap_kpts_real(i,j,k) = dble(ao_overlap_kpts(i,j,k))
enddo
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]

92
src/hartree_fock/scf_k_real.irp.f Normal file
View File

@ -0,0 +1,92 @@
program scf_k_real
BEGIN_DOC
!
! The :ref:`scf` program performs *Restricted* Hartree-Fock
! calculations (the spatial part of the |MOs| is common for alpha and beta
! spinorbitals).
!
! It performs the following actions:
!
! #. Compute/Read all the one- and two-electron integrals, and store them
! in memory
! #. Check in the |EZFIO| database if there is a set of |MOs|.
! If there is, it will read them as initial guess. Otherwise, it will
! create a guess.
! #. Perform the |SCF| iterations
!
! For the keywords related to the |SCF| procedure, see the ``scf_utils``
! directory where you will find all options.
!
! At each iteration, the |MOs| are saved in the |EZFIO| database. Hence,
! if the calculation crashes for any unexpected reason, the calculation
! can be restarted by running again the |SCF| with the same |EZFIO|
! database.
!
! To start again a fresh |SCF| calculation, the |MOs| can be reset by
! running the :ref:`qp_reset` command.
!
! The `DIIS`_ algorithm is implemented, as well as the `level-shifting`_
! method. If the |SCF| does not converge, try again with a higher value of
! :option:`level_shift`.
!
! .. _DIIS: https://en.wikipedia.org/w/index.php?title=DIIS
! .. _level-shifting: https://doi.org/10.1002/qua.560070407
!
END_DOC
call create_guess_k_real
call orthonormalize_mos_k_real
call run_k_real
end
subroutine create_guess_k_real
implicit none
BEGIN_DOC
! Create a MO guess if no MOs are present in the EZFIO directory
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_mo_basis_mo_coef_kpts(exists)
if (.not.exists) then
if (mo_guess_type == "HCore") then
!mo_coef_complex = ao_ortho_lowdin_coef_complex
mo_coef_kpts = ao_ortho_lowdin_coef_kpts_real
TOUCH mo_coef_kpts
mo_label = 'Guess'
!call mo_as_eigvectors_of_mo_matrix_complex(mo_one_e_integrals_kpts, &
call mo_as_eigvectors_of_mo_matrix_kpts_real(mo_one_e_integrals_kpts_real, &
size(mo_one_e_integrals_kpts_real,1), &
size(mo_one_e_integrals_kpts_real,2), &
size(mo_one_e_integrals_kpts_real,3), &
mo_label,1,.false.)
SOFT_TOUCH mo_coef_kpts mo_label
else if (mo_guess_type == "Huckel") then
call huckel_guess_kpts_real
else
print *, 'Unrecognized MO guess type : '//mo_guess_type
stop 1
endif
endif
end
subroutine run_k_real
BEGIN_DOC
! Run SCF calculation
END_DOC
use bitmasks
implicit none
integer :: i_it, i, j, k
mo_label = "Orthonormalized"
call roothaan_hall_scf_kpts_real
call ezfio_set_hartree_fock_energy(SCF_energy)
print*,'hf 1e,2e,total energy'
print*,hf_one_electron_energy
print*,hf_two_electron_energy
print*,hf_energy
end

79
src/mo_basis/utils_cplx.irp.f
View File

@ -327,6 +327,85 @@ subroutine mo_as_eigvectors_of_mo_matrix_kpts(matrix,n,m,nk,label,sign,output)
endif
end
subroutine mo_as_eigvectors_of_mo_matrix_kpts_real(matrix,n,m,nk,label,sign,output)
!TODO: test this
implicit none
integer,intent(in) :: n,m,nk, sign
character*(64), intent(in) :: label
double precision, intent(in) :: matrix(n,m,nk)
logical, intent(in) :: output
integer :: i,j,k
double precision, allocatable :: eigvalues(:)
!complex*16, allocatable :: mo_coef_new(:,:)
double precision, allocatable :: mo_coef_new(:,:),mo_coef_tmp(:,:),R(:,:), A(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: mo_coef_new, R
call write_time(6)
if (m /= mo_num_per_kpt) then
print *, irp_here, ': Error : m/= mo_num_per_kpt'
stop 1
endif
if (nk /= kpt_num) then
print *, irp_here, ': Error : nk/= kpt_num'
stop 1
endif
allocate(A(n,m),R(n,m),mo_coef_tmp(ao_num_per_kpt,m),mo_coef_new(ao_num_per_kpt,m),eigvalues(m))
do k=1,nk
if (sign == -1) then
do j=1,m
do i=1,n
A(i,j) = -matrix(i,j,k)
enddo
enddo
else
do j=1,m
do i=1,n
A(i,j) = matrix(i,j,k)
enddo
enddo
endif
mo_coef_new = dble(mo_coef_kpts(:,:,k))
call lapack_diag(eigvalues,R,A,n,m)
if (sign == -1) then
do i=1,m
eigvalues(i) = -eigvalues(i)
enddo
endif
if (output) then
do i=1,m
write (6,'(2(I8),1X,F16.10)') k,i,eigvalues(i)
enddo
write (6,'(A)') '======== ================'
write (6,'(A)') ''
!write (6,'(A)') 'Fock Matrix'
!write (6,'(A)') '-----------'
!do i=1,n
! write(*,'(200(E24.15))') A(i,:)
!enddo
endif
call dgemm('N','N',ao_num_per_kpt,m,m,1.d0, &
mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0, &
mo_coef_tmp,size(mo_coef_tmp,1))
call zlacp2('N',ao_num_per_kpt,m,mo_coef_tmp,size(mo_coef_tmp,1), &
mo_coef_kpts(:,:,k),size(mo_coef_kpts,1))
enddo
deallocate(A,mo_coef_new,mo_coef_tmp,R,eigvalues)
call write_time(6)
mo_label = label
if (output) then
write (6,'(A)') 'MOs are now **'//trim(label)//'**'
write (6,'(A)') ''
write (6,'(A)') 'Eigenvalues'
write (6,'(A)') '-----------'
write (6,'(A)') ''
write (6,'(A)') '======== ================'
endif
end
subroutine mo_as_svd_vectors_of_mo_matrix_kpts(matrix,lda,m,n,label)
!TODO: implement
print *, irp_here, ' not implemented for kpts'

32
src/mo_guess/mo_ortho_lowdin_cplx.irp.f
View File

@ -107,3 +107,35 @@ BEGIN_PROVIDER [complex*16, ao_ortho_lowdin_overlap_kpts, (ao_num_per_kpt,ao_num
enddo
enddo
END_PROVIDER
!============================================!
! !
! kpts_real !
! !
!============================================!
BEGIN_PROVIDER [ double precision, ao_ortho_lowdin_coef_kpts_real, (ao_num_per_kpt,ao_num_per_kpt,kpt_num)]
implicit none
BEGIN_DOC
! matrix of the coefficients of the mos generated by the
! orthonormalization by the S^{-1/2} canonical transformation of the aos
! ao_ortho_lowdin_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_lowdin orbital
END_DOC
integer :: i,j,k,l
double precision, allocatable :: tmp_matrix(:,:)
allocate (tmp_matrix(ao_num,ao_num))
do k=1,kpt_num
tmp_matrix(:,:) = 0.d0
do j=1, ao_num
tmp_matrix(j,j) = 1.d0
enddo
call ortho_lowdin(ao_overlap_kpts_real(:,:,k),ao_num_per_kpt,ao_num_per_kpt,tmp_matrix,ao_num_per_kpt,ao_num_per_kpt,lin_dep_cutoff)
do i=1, ao_num_per_kpt
do j=1, ao_num_per_kpt
ao_ortho_lowdin_coef_kpts_real(j,i,k) = tmp_matrix(i,j)
enddo
enddo
enddo
deallocate(tmp_matrix)
END_PROVIDER

17
src/mo_one_e_ints/mo_one_e_ints_cplx.irp.f
View File

@ -59,3 +59,20 @@ BEGIN_PROVIDER [ complex*16, mo_one_e_integrals_kpts,(mo_num_per_kpt,mo_num_per_
print*,'Provided the one-electron integrals'
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_one_e_integrals_kpts_real,(mo_num_per_kpt,mo_num_per_kpt,kpt_num)]
implicit none
BEGIN_DOC
! array of the one-electron Hamiltonian on the |MO| basis :
! sum of the kinetic and nuclear electronic potentials (and pseudo potential if needed)
END_DOC
integer :: i,j,k
do k=1,kpt_num
do j=1,mo_num_per_kpt
do i=1,mo_num_per_kpt
mo_one_e_integrals_kpts_real(i,j,k) = dble(mo_one_e_integrals_kpts(i,j,k))
enddo
enddo
enddo
END_PROVIDER

15
src/mo_one_e_ints/mo_overlap.irp.f
View File

@ -128,3 +128,18 @@ BEGIN_PROVIDER [ complex*16, mo_overlap_kpts,(mo_num_per_kpt,mo_num_per_kpt,kpt_
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_overlap_kpts_real, (mo_num_per_kpt, mo_num_per_kpt, kpt_num) ]
implicit none
BEGIN_DOC
! Overlap for complex MOs
END_DOC
integer :: i,j,k
do k=1,kpt_num
do j=1,mo_num_per_kpt
do i=1,mo_num_per_kpt
mo_overlap_kpts_real(i,j,k) = dble(mo_overlap_kpts(i,j,k))
enddo
enddo
enddo
END_PROVIDER

20
src/mo_one_e_ints/orthonormalize.irp.f
View File

@ -19,3 +19,23 @@ subroutine orthonormalize_mos
end
subroutine orthonormalize_mos_k_real
implicit none
integer :: m,p,s,k
double precision, allocatable :: mo_coef_tmp(:,:)
allocate(mo_coef_tmp(ao_num_per_kpt,mo_num_per_kpt))
do k=1,kpt_num
m = size(mo_coef_kpts,1)
p = size(mo_overlap_kpts,1)
mo_coef_tmp = dble(mo_coef_kpts(:,:,k))
call ortho_lowdin(mo_overlap_kpts_real(1,1,k),p,mo_num_per_kpt,mo_coef_tmp,m,ao_num_per_kpt,lin_dep_cutoff)
call zlacp2('X',ao_num_per_kpt,mo_num_per_kpt,mo_coef_tmp,size(mo_coef_tmp,1), &
mo_coef_kpts(1,1,k),size(mo_coef_kpts,1))
enddo
deallocate(mo_coef_tmp)
mo_label = 'Orthonormalized'
SOFT_TOUCH mo_coef_kpts mo_label
end

67
src/scf_utils/diagonalize_fock_cplx.irp.f
View File

@ -112,4 +112,71 @@ BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts, (ao_num_per_kpt,m
deallocate(F, diag)
END_PROVIDER
BEGIN_PROVIDER [ complex*16, eigenvectors_Fock_matrix_mo_kpts_real, (ao_num_per_kpt,mo_num_per_kpt,kpt_num) ]
implicit none
BEGIN_DOC
! Eigenvectors of the Fock matrix in the |MO| basis obtained with level shift.
END_DOC
integer :: i,j,k
integer :: n
!complex*16, allocatable :: F(:,:)
double precision, allocatable :: F(:,:)
double precision, allocatable :: diag(:), mo_coef_tmp(:,:), eigvecs_tmp(:,:)
allocate( F(mo_num_per_kpt,mo_num_per_kpt) )
allocate (diag(mo_num_per_kpt) )
allocate (mo_coef_tmp(ao_num_per_kpt,mo_num_per_kpt) )
allocate (eigvecs_tmp(ao_num_per_kpt,mo_num_per_kpt) )
do k=1,kpt_num
do j=1,mo_num_per_kpt
do i=1,mo_num_per_kpt
!F(i,j) = fock_matrix_mo_complex(i,j)
F(i,j) = dble(fock_matrix_mo_kpts(i,j,k))
enddo
enddo
if(frozen_orb_scf)then
integer :: iorb,jorb
!todo: core/act per kpt
do i = 1, n_core_orb
iorb = list_core(i)
do j = 1, n_act_orb
jorb = list_act(j)
F(iorb,jorb) = 0.d0
F(jorb,iorb) = 0.d0
enddo
enddo
endif
! Insert level shift here
!todo: elec per kpt
do i = elec_beta_num_kpts(k)+1, elec_alpha_num_kpts(k)
F(i,i) += 0.5d0*level_shift
enddo
do i = elec_alpha_num_kpts(k)+1, mo_num_per_kpt
F(i,i) += level_shift
enddo
n = mo_num_per_kpt
call lapack_diagd_diag_in_place(diag,F,n,n)
mo_coef_tmp = dble(mo_coef_kpts(:,:,k))
call dgemm('N','N',ao_num_per_kpt,mo_num_per_kpt,mo_num_per_kpt, 1.d0, &
mo_coef_tmp, size(mo_coef_tmp,1), F, size(F,1), &
0.d0, eigvecs_tmp, size(eigvecs_tmp,1))
call zlacp2('X',ao_num_per_kpt,mo_num_per_kpt,eigvecs_tmp,size(eigvecs_tmp,1), &
eigenvectors_fock_matrix_mo_kpts_real(:,:,k), size(eigenvectors_Fock_matrix_mo_kpts_real,1))
! call zgemm('N','N',ao_num_per_kpt,mo_num_per_kpt,mo_num_per_kpt, (1.d0,0.d0), &
! mo_coef_kpts(:,:,k), size(mo_coef_kpts,1), F, size(F,1), &
! (0.d0,0.d0), eigenvectors_Fock_matrix_mo_kpts(:,:,k), size(eigenvectors_Fock_matrix_mo_kpts,1))
enddo
deallocate(F, diag,mo_coef_tmp,eigvecs_tmp)
END_PROVIDER

2
src/scf_utils/diis_cplx.irp.f
View File

@ -164,7 +164,7 @@ BEGIN_PROVIDER [complex*16, FPS_SPF_Matrix_AO_kpts, (AO_num_per_kpt, AO_num_per_
call zgemm('N','N',AO_num_per_kpt,AO_num_per_kpt,AO_num_per_kpt, &
(1.d0,0.d0), &
Fock_Matrix_AO_kpts(1,1,k),Size(Fock_Matrix_AO_kpts,1), &
SCF_Density_Matrix_AO_kpts(1,1,k),Size(SCF_Density_Matrix_AO_kpts,1), &
scf_density_matrix_ao_kpts(1,1,k),Size(SCF_Density_Matrix_AO_kpts,1), &
(0.d0,0.d0), &
scratch,Size(scratch,1))

18
src/scf_utils/fock_matrix_cplx.irp.f
View File

@ -360,6 +360,24 @@ END_PROVIDER
END_PROVIDER
!============================================!
! !
! kpts_real !
! !
!============================================!
BEGIN_PROVIDER [ double precision, Fock_matrix_mo_kpts_real, (mo_num_per_kpt,mo_num_per_kpt,kpt_num) ]
implicit none
integer :: i,j,k
do k=1,kpt_num
do j=1,mo_num_per_kpt
do i=1,mo_num_per_kpt
fock_matrix_mo_kpts_real(i,j,k) = dble(fock_matrix_mo_kpts(i,j,k))
enddo
enddo
enddo
END_PROVIDER
!============================================!
! !
! kpts !

49
src/scf_utils/huckel_cplx.irp.f
View File

@ -89,3 +89,52 @@ subroutine huckel_guess_kpts
deallocate(A)
end
subroutine huckel_guess_kpts_real
implicit none
BEGIN_DOC
! Build the MOs using the extended Huckel model
END_DOC
integer :: i,j,k
double precision :: accu
double precision :: c
character*(64) :: label
!complex*16, allocatable :: A(:,:)
double precision, allocatable :: A(:,:)
label = "Guess"
c = 0.5d0 * 1.75d0
allocate (A(ao_num_per_kpt, ao_num_per_kpt))
do k=1,kpt_num
A = 0.d0
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
A(i,j) = c * ao_overlap_kpts_real(i,j,k) * (ao_one_e_integrals_diag_kpts(i,k) + ao_one_e_integrals_diag_kpts(j,k))
enddo
A(j,j) = ao_one_e_integrals_diag_kpts(j,k) + dble(ao_two_e_integral_alpha_kpts(j,j,k))
if (dabs(dimag(ao_two_e_integral_alpha_kpts(j,j,k))) .gt. 1.0d-10) then
stop 'diagonal elements of ao_two_e_integral_alpha should be real'
endif
enddo
! Fock_matrix_ao_alpha(1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
! Fock_matrix_ao_beta (1:ao_num,1:ao_num) = A(1:ao_num,1:ao_num)
call zlacp2('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), &
Fock_matrix_ao_alpha_kpts(:,:,k), size(Fock_matrix_ao_alpha_kpts,1))
call zlacp2('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), &
Fock_matrix_ao_beta_kpts(:,:,k), size(Fock_matrix_ao_beta_kpts, 1))
!call zlacpy('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), &
! Fock_matrix_ao_alpha_kpts(:,:,k), size(Fock_matrix_ao_alpha_kpts,1))
!call zlacpy('X', ao_num_per_kpt, ao_num_per_kpt, A, size(A,1), &
! Fock_matrix_ao_beta_kpts(:,:,k), size(Fock_matrix_ao_beta_kpts, 1))
enddo
! TOUCH mo_coef
!TOUCH fock_matrix_ao_alpha_complex fock_matrix_ao_beta_kpts
TOUCH fock_matrix_ao_alpha_kpts fock_matrix_ao_beta_kpts
mo_coef_kpts = eigenvectors_fock_matrix_mo_kpts_real
SOFT_TOUCH mo_coef_kpts
call save_mos
deallocate(A)
end

189
src/scf_utils/roothaan_hall_scf_cplx.irp.f
View File

@ -653,3 +653,192 @@ END_DOC
endif
end
!============================================!
! !
! kpts_real !
! !
!============================================!
subroutine Roothaan_Hall_SCF_kpts_real
BEGIN_DOC
! Roothaan-Hall algorithm for SCF Hartree-Fock calculation
END_DOC
implicit none
double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF
double precision :: max_error_DIIS,max_error_DIIS_alpha,max_error_DIIS_beta
complex*16, allocatable :: Fock_matrix_DIIS(:,:,:,:),error_matrix_DIIS(:,:,:,:)
integer :: iteration_SCF,dim_DIIS,index_dim_DIIS
integer :: i,j,k,kk
logical, external :: qp_stop
complex*16, allocatable :: mo_coef_save(:,:,:)
PROVIDE ao_md5 mo_occ_kpts level_shift
allocate(mo_coef_save(ao_num_per_kpt,mo_num_per_kpt,kpt_num), &
Fock_matrix_DIIS (ao_num_per_kpt,ao_num_per_kpt,max_dim_DIIS,kpt_num), &
error_matrix_DIIS(ao_num_per_kpt,ao_num_per_kpt,max_dim_DIIS,kpt_num) &
)
!todo: add kpt_num dim to diis mats? (3 or 4)
call write_time(6)
print*,'Energy of the guess = ',scf_energy
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
' N ', 'Energy ', 'Energy diff ', 'DIIS error ', 'Level shift '
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
! Initialize energies and density matrices
energy_SCF_previous = SCF_energy
Delta_energy_SCF = 1.d0
iteration_SCF = 0
dim_DIIS = 0
max_error_DIIS = 1.d0
!
! Start of main SCF loop
!
!PROVIDE fps_spf_matrix_ao_complex fock_matrix_ao_complex
PROVIDE fps_spf_matrix_ao_kpts fock_matrix_ao_kpts
do while ( &
( (max_error_DIIS > threshold_DIIS_nonzero) .or. &
(dabs(Delta_energy_SCF) > thresh_SCF) &
) .and. (iteration_SCF < n_it_SCF_max) )
! Increment cycle number
iteration_SCF += 1
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)
if (scf_algorithm == 'DIIS') then
do kk=1,kpt_num
! Store Fock and error matrices at each iteration
do j=1,ao_num_per_kpt
do i=1,ao_num_per_kpt
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
Fock_matrix_DIIS (i,j,index_dim_DIIS,kk) = fock_matrix_ao_kpts(i,j,kk)
error_matrix_DIIS(i,j,index_dim_DIIS,kk) = fps_spf_matrix_ao_kpts(i,j,kk)
enddo
enddo
! Compute the extrapolated Fock matrix
call extrapolate_fock_matrix_kpts( &
error_matrix_DIIS(1,1,1,kk),Fock_matrix_DIIS(1,1,1,kk), &
Fock_matrix_AO_kpts(1,1,kk),size(Fock_matrix_AO_kpts,1), &
iteration_SCF,dim_DIIS &
)
enddo
Fock_matrix_AO_alpha_kpts = Fock_matrix_AO_kpts*0.5d0
Fock_matrix_AO_beta_kpts = Fock_matrix_AO_kpts*0.5d0
TOUCH Fock_matrix_AO_alpha_kpts Fock_matrix_AO_beta_kpts
endif
mo_coef_kpts = eigenvectors_fock_matrix_mo_kpts
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_kpts
! Calculate error vectors
max_error_DIIS = maxval(cdabs(FPS_SPF_Matrix_MO_kpts))
! SCF energy
! call print_debug_scf_complex
energy_SCF = scf_energy
Delta_Energy_SCF = energy_SCF - energy_SCF_previous
if ( (SCF_algorithm == 'DIIS').and.(Delta_Energy_SCF > 0.d0) ) then
do kk=1,kpt_num
Fock_matrix_AO_kpts(1:ao_num_per_kpt,1:ao_num_per_kpt,kk) = &
Fock_matrix_DIIS (1:ao_num_per_kpt,1:ao_num_per_kpt,index_dim_DIIS,kk)
enddo
Fock_matrix_AO_alpha_kpts = Fock_matrix_AO_kpts*0.5d0
Fock_matrix_AO_beta_kpts = Fock_matrix_AO_kpts*0.5d0
TOUCH fock_matrix_ao_alpha_kpts Fock_matrix_AO_beta_kpts
endif
double precision :: level_shift_save
level_shift_save = level_shift
mo_coef_save(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num)
do while (Delta_energy_SCF > 0.d0)
mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = mo_coef_save
if (level_shift <= .1d0) then
level_shift = 1.d0
else
level_shift = level_shift * 3.0d0
endif
TOUCH mo_coef_kpts level_shift
mo_coef_kpts(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num) = &
eigenvectors_fock_matrix_mo_kpts_real(1:ao_num_per_kpt,1:mo_num_per_kpt,1:kpt_num)
if(frozen_orb_scf)then
call reorder_core_orb
call initialize_mo_coef_begin_iteration
endif
TOUCH mo_coef_kpts
Delta_Energy_SCF = SCF_energy - energy_SCF_previous
energy_SCF = SCF_energy
if (level_shift-level_shift_save > 40.d0) then
level_shift = level_shift_save * 4.d0
SOFT_TOUCH level_shift
exit
endif
dim_DIIS=0
enddo
level_shift = level_shift * 0.5d0
SOFT_TOUCH level_shift
energy_SCF_previous = energy_SCF
! Print results at the end of each iteration
write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') &
iteration_SCF, energy_scf, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS
if (Delta_energy_SCF < 0.d0) then
call save_mos
endif
if (qp_stop()) exit
enddo
if (iteration_SCF < n_it_SCF_max) then
mo_label = "Canonical"
endif
!
! End of Main SCF loop
!
write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)') &
'====','================','================','================','================'
write(6,*)
if(.not.frozen_orb_scf)then
call mo_as_eigvectors_of_mo_matrix_kpts_real(fock_matrix_mo_kpts_real,size(Fock_matrix_mo_kpts_real,1),size(Fock_matrix_mo_kpts_real,2),size(Fock_matrix_mo_kpts_real,3),mo_label,1,.true.)
call save_mos
endif
call write_double(6, Energy_SCF, 'SCF energy')
call write_time(6)
end

74
src/utils/linear_algebra.irp.f
View File

@ -1277,3 +1277,77 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
deallocate(A,eigenvalues)
end
subroutine lapack_diagd_diag_in_place(eigvalues,eigvectors,nmax,n)
implicit none
BEGIN_DOC
! Diagonalize matrix H(complex)
!
! H is untouched between input and ouptut
!
! eigevalues(i) = ith lowest eigenvalue of the H matrix
!
! eigvectors(i,j) = <i|psi_j> where i is the basis function and psi_j is the j th eigenvector
!
END_DOC
integer, intent(in) :: n,nmax
double precision, intent(out) :: eigvectors(nmax,n)
! complex*16, intent(inout) :: eigvectors(nmax,n)
double precision, intent(out) :: eigvalues(n)
! double precision, intent(in) :: H(nmax,n)
double precision,allocatable :: work(:)
integer ,allocatable :: iwork(:)
! complex*16,allocatable :: A(:,:)
integer :: lwork, info, i,j,l,k, liwork
! print*,'Diagonalization by jacobi'
! print*,'n = ',n
lwork = 2*n*n + 6*n + 1
liwork = 5*n + 3
allocate (work(lwork),iwork(liwork))
lwork = -1
liwork = -1
! get optimal work size
call DSYEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, &
iwork, liwork, info )
if (info < 0) then
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
endif
lwork = int( real(work(1)))
liwork = iwork(1)
deallocate (work,iwork)
allocate (work(lwork),iwork(liwork))
call DSYEVD( 'V', 'U', n, eigvectors, nmax, eigvalues, work, lwork, &
iwork, liwork, info )
deallocate(work,iwork)
if (info < 0) then
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
else if( info > 0 ) then
write(*,*)'DSYEVD Failed; calling DSYEV'
lwork = 3*n - 1
allocate(work(lwork))
lwork = -1
call DSYEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,info)
if (info < 0) then
print *, irp_here, ': DSYEV: the ',-info,'-th argument had an illegal value'
stop 2
endif
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call DSYEV('V','L',n,eigvectors,nmax,eigvalues,work,lwork,info)
if (info /= 0 ) then
write(*,*)'DSYEV Failed'
stop 1
endif
deallocate(work)
end if
end