mirror of
https://github.com/LCPQ/quantum_package
synced 2024-12-22 20:35:19 +01:00
Finalized DIIS
This commit is contained in:
parent
54b5e86608
commit
de4ba1961f
@ -1,68 +1,74 @@
|
||||
begin_template
|
||||
BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! If threshold_DIIS is zero, choose sqrt(thresh_scf)
|
||||
END_DOC
|
||||
if (threshold_DIIS == 0.d0) then
|
||||
threshold_DIIS_nonzero = dsqrt(thresh_scf)
|
||||
else
|
||||
threshold_DIIS_nonzero = threshold_DIIS
|
||||
endif
|
||||
ASSERT (threshold_DIIS_nonzero >= 0.d0)
|
||||
|
||||
begin_provider [double precision, FPS_SPF_Matrix_AO_$alpha, (AO_num, AO_num)]
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [double precision, FPS_SPF_Matrix_AO, (AO_num, AO_num)]
|
||||
implicit none
|
||||
begin_doc
|
||||
! Commutator FPS - SPF
|
||||
end_doc
|
||||
double precision, allocatable:: scratch(:,:)
|
||||
allocate( &
|
||||
scratch(AO_num_align, AO_num) &
|
||||
)
|
||||
BEGIN_DOC
|
||||
! Commutator FPS - SPF
|
||||
END_DOC
|
||||
double precision, allocatable :: scratch(:,:)
|
||||
allocate( &
|
||||
scratch(AO_num_align, AO_num) &
|
||||
)
|
||||
|
||||
! Compute FP
|
||||
! Compute FP
|
||||
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
Fock_Matrix_AO_$alpha,Size(Fock_Matrix_AO_$alpha,1), &
|
||||
HF_Density_Matrix_AO_$alpha,Size(HF_Density_Matrix_AO_$alpha,1), &
|
||||
0.d0, &
|
||||
scratch,Size(scratch,1))
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
Fock_Matrix_AO,Size(Fock_Matrix_AO,1), &
|
||||
HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), &
|
||||
0.d0, &
|
||||
scratch,Size(scratch,1))
|
||||
|
||||
! Compute FPS
|
||||
! Compute FPS
|
||||
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
scratch,Size(scratch,1), &
|
||||
AO_Overlap,Size(AO_Overlap,1), &
|
||||
0.d0, &
|
||||
FPS_SPF_Matrix_AO_$alpha,Size(FPS_SPF_Matrix_AO_$alpha,1))
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
scratch,Size(scratch,1), &
|
||||
AO_Overlap,Size(AO_Overlap,1), &
|
||||
0.d0, &
|
||||
FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
|
||||
|
||||
! Compute SP
|
||||
! Compute SP
|
||||
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
AO_Overlap,Size(AO_Overlap,1), &
|
||||
HF_Density_Matrix_AO_$alpha,Size(HF_Density_Matrix_AO_$alpha,1), &
|
||||
0.d0, &
|
||||
scratch,Size(scratch,1))
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
1.d0, &
|
||||
AO_Overlap,Size(AO_Overlap,1), &
|
||||
HF_Density_Matrix_AO,Size(HF_Density_Matrix_AO,1), &
|
||||
0.d0, &
|
||||
scratch,Size(scratch,1))
|
||||
|
||||
! Compute FPS - SPF
|
||||
! Compute FPS - SPF
|
||||
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
-1.d0, &
|
||||
scratch,Size(scratch,1), &
|
||||
Fock_Matrix_AO_$alpha,Size(Fock_Matrix_AO_$alpha,1), &
|
||||
1.d0, &
|
||||
FPS_SPF_Matrix_AO_$alpha,Size(FPS_SPF_Matrix_AO_$alpha,1))
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
-1.d0, &
|
||||
scratch,Size(scratch,1), &
|
||||
Fock_Matrix_AO,Size(Fock_Matrix_AO,1), &
|
||||
1.d0, &
|
||||
FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
|
||||
|
||||
end_provider
|
||||
END_PROVIDER
|
||||
|
||||
begin_provider [double precision, FPS_SPF_Matrix_MO_$alpha, (mo_tot_num, mo_tot_num)]
|
||||
bEGIN_PROVIDER [double precision, FPS_SPF_Matrix_MO, (mo_tot_num, mo_tot_num)]
|
||||
implicit none
|
||||
begin_doc
|
||||
! Commutator FPS - SPF in MO basis
|
||||
end_doc
|
||||
call ao_to_mo(FPS_SPF_Matrix_AO_$alpha, size(FPS_SPF_Matrix_AO_$alpha,1), &
|
||||
FPS_SPF_Matrix_MO_$alpha, size(FPS_SPF_Matrix_MO_$alpha,1))
|
||||
end_provider
|
||||
call ao_to_mo(FPS_SPF_Matrix_AO, size(FPS_SPF_Matrix_AO,1), &
|
||||
FPS_SPF_Matrix_MO, size(FPS_SPF_Matrix_MO,1))
|
||||
END_PROVIDER
|
||||
|
||||
subst [alpha]
|
||||
|
||||
alpha ;;
|
||||
beta ;;
|
||||
|
||||
end_template
|
||||
|
||||
BEGIN_PROVIDER [ double precision, eigenvalues_Fock_matrix_AO, (AO_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, eigenvectors_Fock_matrix_AO, (AO_num_align,AO_num) ]
|
||||
@ -78,11 +84,11 @@ end_template
|
||||
integer :: i,j
|
||||
|
||||
lwork = 3*AO_num - 1
|
||||
allocate( &
|
||||
scratch(AO_num_align,AO_num), &
|
||||
work(lwork), &
|
||||
Xt(AO_num,AO_num) &
|
||||
)
|
||||
allocate( &
|
||||
scratch(AO_num_align,AO_num), &
|
||||
work(lwork), &
|
||||
Xt(AO_num,AO_num) &
|
||||
)
|
||||
|
||||
! Calculate Xt
|
||||
|
||||
@ -162,10 +168,11 @@ BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num_align,AO_num) ]
|
||||
num_linear_dependencies = 0
|
||||
do i=1,AO_num
|
||||
print*,D(i)
|
||||
if(abs(D(i)) < threshold_overlap_AO_eigenvalues) then
|
||||
if(abs(D(i)) <= threshold_overlap_AO_eigenvalues) then
|
||||
D(i) = 0.d0
|
||||
num_linear_dependencies += 1
|
||||
else
|
||||
ASSERT (D(i) > 0.d0)
|
||||
D(i) = 1.d0/sqrt(D(i))
|
||||
endif
|
||||
do j=1,AO_num
|
||||
|
||||
@ -12,13 +12,13 @@ default: 15
|
||||
|
||||
[threshold_diis]
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the DIIS error vector during a Hartree-Fock calculation
|
||||
doc: Threshold on the convergence of the DIIS error vector during a Hartree-Fock calculation. If 0. is chosen, the square root of thresh_scf will be used.
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-6
|
||||
default: 0.
|
||||
|
||||
[thresh_scf]
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the Hartree Fock energy
|
||||
doc: Threshold on the convergence of the Hartree Fock energy.
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-12
|
||||
|
||||
@ -32,11 +32,11 @@ default: 128
|
||||
type: Positive_float
|
||||
doc: Energy shift on the virtual MOs to improve SCF convergence
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.4
|
||||
default: 0.0
|
||||
|
||||
[scf_algorithm]
|
||||
type: character*(32)
|
||||
doc: Type of SCF algorithm used. Possible choices are [ Damp | DIIS]
|
||||
doc: Type of SCF algorithm used. Possible choices are [ Simple | DIIS]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: DIIS
|
||||
|
||||
|
||||
@ -8,20 +8,16 @@ END_DOC
|
||||
|
||||
double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF
|
||||
double precision :: max_error_DIIS,max_error_DIIS_alpha,max_error_DIIS_beta
|
||||
double precision, allocatable :: Fock_matrix_DIIS_alpha(:,:,:),error_matrix_DIIS_alpha(:,:,:)
|
||||
double precision, allocatable :: Fock_matrix_DIIS_beta (:,:,:),error_matrix_DIIS_beta (:,:,:)
|
||||
double precision, allocatable :: Fock_matrix_DIIS(:,:,:),error_matrix_DIIS(:,:,:)
|
||||
|
||||
integer :: iteration_SCF,dim_DIIS,index_dim_DIIS
|
||||
integer :: dim_DIIS_alpha, dim_DIIS_beta
|
||||
|
||||
integer :: i,j
|
||||
|
||||
allocate( &
|
||||
Fock_matrix_DIIS_alpha(ao_num,ao_num,max_dim_DIIS), &
|
||||
Fock_matrix_DIIS_beta (ao_num,ao_num,max_dim_DIIS), &
|
||||
error_matrix_DIIS_alpha(ao_num,ao_num,max_dim_DIIS), &
|
||||
error_matrix_DIIS_beta (ao_num,ao_num,max_dim_DIIS) &
|
||||
)
|
||||
allocate( &
|
||||
Fock_matrix_DIIS (ao_num,ao_num,max_dim_DIIS), &
|
||||
error_matrix_DIIS(ao_num,ao_num,max_dim_DIIS) &
|
||||
)
|
||||
|
||||
call write_time(output_hartree_fock)
|
||||
|
||||
@ -35,17 +31,16 @@ END_DOC
|
||||
! Initialize energies and density matrices
|
||||
|
||||
energy_SCF_previous = HF_energy
|
||||
Delta_energy_SCF = 0.d0
|
||||
iteration_SCF = 0
|
||||
dim_DIIS_alpha = 0
|
||||
dim_DIIS_beta = 0
|
||||
dim_DIIS = 0
|
||||
max_error_DIIS = 1.d0
|
||||
Delta_energy_SCF = 1.d0
|
||||
iteration_SCF = 0
|
||||
dim_DIIS = 0
|
||||
max_error_DIIS = 1.d0
|
||||
|
||||
|
||||
!
|
||||
! Start of main SCF loop
|
||||
!
|
||||
do while((max_error_DIIS > threshold_DIIS) .and. (iteration_SCF < n_it_SCF_max))
|
||||
do while((max_error_DIIS > threshold_DIIS_nonzero) .and. (iteration_SCF < n_it_SCF_max) .and. dabs(Delta_energy_SCF) > thresh_SCF)
|
||||
|
||||
! Increment cycle number
|
||||
|
||||
@ -55,46 +50,38 @@ END_DOC
|
||||
|
||||
dim_DIIS = min(dim_DIIS+1,max_dim_DIIS)
|
||||
|
||||
! Store Fock and error matrices at each iteration
|
||||
if (scf_algorithm == 'DIIS') then
|
||||
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
|
||||
Fock_matrix_DIIS_alpha (i,j,index_dim_DIIS) = Fock_matrix_AO_alpha(i,j)
|
||||
Fock_matrix_DIIS_beta (i,j,index_dim_DIIS) = Fock_matrix_AO_beta (i,j)
|
||||
error_matrix_DIIS_alpha(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO_alpha(i,j)
|
||||
error_matrix_DIIS_beta (i,j,index_dim_DIIS) = FPS_SPF_matrix_AO_beta (i,j)
|
||||
! Store Fock and error matrices at each iteration
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
index_dim_DIIS = mod(dim_DIIS-1,max_dim_DIIS)+1
|
||||
Fock_matrix_DIIS (i,j,index_dim_DIIS) = Fock_matrix_AO(i,j)
|
||||
error_matrix_DIIS(i,j,index_dim_DIIS) = FPS_SPF_matrix_AO(i,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Compute the extrapolated Fock matrix
|
||||
! Compute the extrapolated Fock matrix
|
||||
|
||||
dim_DIIS_alpha = dim_DIIS
|
||||
call extrapolate_Fock_matrix( &
|
||||
error_matrix_DIIS_alpha,Fock_matrix_DIIS_alpha, &
|
||||
Fock_matrix_AO_alpha,size(Fock_matrix_AO_alpha,1), &
|
||||
iteration_SCF,dim_DIIS_alpha &
|
||||
)
|
||||
call extrapolate_Fock_matrix( &
|
||||
error_matrix_DIIS,Fock_matrix_DIIS, &
|
||||
Fock_matrix_AO,size(Fock_matrix_AO,1), &
|
||||
iteration_SCF,dim_DIIS &
|
||||
)
|
||||
|
||||
dim_DIIS_beta = dim_DIIS
|
||||
call extrapolate_Fock_matrix( &
|
||||
error_matrix_DIIS_beta,Fock_matrix_DIIS_beta, &
|
||||
Fock_matrix_AO_beta,size(Fock_matrix_AO_beta,1), &
|
||||
iteration_SCF,dim_DIIS_beta &
|
||||
)
|
||||
Fock_matrix_AO_alpha = Fock_matrix_AO*0.5d0
|
||||
Fock_matrix_AO_Beta = Fock_matrix_AO*0.5d0
|
||||
touch Fock_matrix_AO_alpha Fock_matrix_AO_beta
|
||||
|
||||
dim_DIIS = min(dim_DIIS_alpha,dim_DIIS_beta)
|
||||
touch Fock_matrix_AO_alpha Fock_matrix_AO_beta
|
||||
endif
|
||||
|
||||
MO_coef = eigenvectors_Fock_matrix_MO
|
||||
|
||||
touch MO_coef
|
||||
|
||||
! Calculate error vectors
|
||||
! Calculate error vectors
|
||||
|
||||
max_error_DIIS_alpha = maxval(Abs(FPS_SPF_Matrix_MO_alpha))
|
||||
max_error_DIIS_beta = maxval(Abs(FPS_SPF_Matrix_MO_beta ))
|
||||
max_error_DIIS = max(max_error_DIIS_alpha,max_error_DIIS_beta)
|
||||
max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO))
|
||||
|
||||
! SCF energy
|
||||
|
||||
@ -102,11 +89,15 @@ END_DOC
|
||||
Delta_Energy_SCF = energy_SCF - energy_SCF_previous
|
||||
energy_SCF_previous = energy_SCF
|
||||
|
||||
! Print results at the end of each iteration
|
||||
! Print results at the end of each iteration
|
||||
|
||||
write(output_hartree_fock,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)') &
|
||||
iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, dim_DIIS
|
||||
|
||||
if (Delta_energy_SCF < 0.d0) then
|
||||
call save_mos
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
!
|
||||
@ -178,7 +169,7 @@ END_DOC
|
||||
|
||||
B_matrix_DIIS(i,j) = 0.d0
|
||||
do k=1,ao_num
|
||||
B_matrix_DIIS(i,j) += scratch(k,k)
|
||||
B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + scratch(k,k)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -225,29 +216,27 @@ END_DOC
|
||||
stop 'bug in DIIS'
|
||||
endif
|
||||
|
||||
if (rcond > 1.d-8) then
|
||||
if (rcond > 1.d-14) then
|
||||
|
||||
! Compute extrapolated Fock matrix
|
||||
! Compute extrapolated Fock matrix
|
||||
|
||||
Fock_matrix_AO_(:,:) = 0.d0
|
||||
|
||||
do k=1,dim_DIIS
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED)
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
Fock_matrix_AO_(i,j) += X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)
|
||||
Fock_matrix_AO_(i,j) = 0.d0
|
||||
enddo
|
||||
do k=1,dim_DIIS
|
||||
do i=1,ao_num
|
||||
Fock_matrix_AO_(i,j) = Fock_matrix_AO_(i,j) + &
|
||||
X_vector_DIIS(k)*Fock_matrix_DIIS(i,j,dim_DIIS-k+1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
else
|
||||
write(*,*) 'Re-initialize DIIS!!'
|
||||
dim_DIIS = 0
|
||||
endif
|
||||
|
||||
! do i=1,ao_num
|
||||
! do j=1,ao_num
|
||||
! write(*,*) Fock_matrix_AO_(i,j)
|
||||
! enddo
|
||||
! enddo
|
||||
|
||||
end
|
||||
|
||||
@ -56,24 +56,9 @@ subroutine run
|
||||
|
||||
! Choose SCF algorithm
|
||||
|
||||
if(scf_algorithm == 'Damp') then
|
||||
call damping_SCF
|
||||
else if(scf_algorithm == 'DIIS') then
|
||||
call Roothaan_Hall_SCF
|
||||
else
|
||||
write(*,*) 'Unrecognized SCF algorithm: '//scf_algorithm
|
||||
stop 1
|
||||
endif
|
||||
! call damping_SCF
|
||||
call Roothaan_Hall_SCF
|
||||
|
||||
end
|
||||
|
||||
|
||||
subroutine debug
|
||||
implicit none
|
||||
integer :: i,j
|
||||
do i=1,mo_tot_num
|
||||
do j=1,mo_tot_num
|
||||
print *, i, j, fps_spf_matrix_mo_alpha(i,j)
|
||||
enddo
|
||||
enddo
|
||||
end
|
||||
|
||||
Loading…
Reference in New Issue
Block a user