mirror of
https://github.com/LCPQ/quantum_package
synced 2025-01-03 10:05:57 +01:00
Almost working
This commit is contained in:
parent
f7da26ff09
commit
e15fbd2371
@ -1,4 +1,6 @@
|
||||
begin_provider [double precision, FPS_SPF_Matrix_AO, (AO_num_align, AO_num)]
|
||||
begin_template
|
||||
|
||||
begin_provider [double precision, FPS_SPF_Matrix_AO_$alpha, (AO_num, AO_num)]
|
||||
implicit none
|
||||
begin_doc
|
||||
! Commutator FPS - SPF
|
||||
@ -12,8 +14,8 @@ begin_provider [double precision, FPS_SPF_Matrix_AO, (AO_num_align, AO_num)]
|
||||
|
||||
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), &
|
||||
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))
|
||||
|
||||
@ -24,14 +26,14 @@ begin_provider [double precision, FPS_SPF_Matrix_AO, (AO_num_align, AO_num)]
|
||||
scratch,Size(scratch,1), &
|
||||
AO_Overlap,Size(AO_Overlap,1), &
|
||||
0.d0, &
|
||||
FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
|
||||
FPS_SPF_Matrix_AO_$alpha,Size(FPS_SPF_Matrix_AO_$alpha,1))
|
||||
|
||||
! Compute SP
|
||||
|
||||
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), &
|
||||
HF_Density_Matrix_AO_$alpha,Size(HF_Density_Matrix_AO_$alpha,1), &
|
||||
0.d0, &
|
||||
scratch,Size(scratch,1))
|
||||
|
||||
@ -40,12 +42,27 @@ begin_provider [double precision, FPS_SPF_Matrix_AO, (AO_num_align, AO_num)]
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
-1.d0, &
|
||||
scratch,Size(scratch,1), &
|
||||
Fock_Matrix_AO,Size(Fock_Matrix_AO,1), &
|
||||
Fock_Matrix_AO_$alpha,Size(Fock_Matrix_AO_$alpha,1), &
|
||||
1.d0, &
|
||||
FPS_SPF_Matrix_AO,Size(FPS_SPF_Matrix_AO,1))
|
||||
FPS_SPF_Matrix_AO_$alpha,Size(FPS_SPF_Matrix_AO_$alpha,1))
|
||||
|
||||
end_provider
|
||||
|
||||
begin_provider [double precision, FPS_SPF_Matrix_MO_$alpha, (AO_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
|
||||
|
||||
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) ]
|
||||
@ -122,12 +139,13 @@ BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num_align,AO_num) ]
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: num_linear_dependencies
|
||||
integer :: LDA, LDC
|
||||
double precision, allocatable :: U(:,:),Vt(:,:), D(:)
|
||||
integer :: info, i, j, k
|
||||
|
||||
LDA = size(AO_overlap,1)
|
||||
LDC = size(AO_overlap,1)
|
||||
LDC = size(X_matrix_AO,1)
|
||||
|
||||
allocate( &
|
||||
U(LDC,AO_num), &
|
||||
@ -141,9 +159,12 @@ BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num_align,AO_num) ]
|
||||
Vt,LDA, &
|
||||
AO_num,AO_num)
|
||||
|
||||
num_linear_dependencies = 0
|
||||
do i=1,AO_num
|
||||
print*,D(i)
|
||||
if(abs(D(i)) < threshold_overlap_AO_eigenvalues) then
|
||||
D(i) = 0.d0
|
||||
num_linear_dependencies += 1
|
||||
else
|
||||
D(i) = 1.d0/sqrt(D(i))
|
||||
endif
|
||||
@ -151,11 +172,13 @@ BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num_align,AO_num) ]
|
||||
X_matrix_AO(j,i) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
write(*,*) 'linear dependencies',num_linear_dependencies
|
||||
! stop
|
||||
|
||||
do k=1,AO_num
|
||||
if(D(k) /= 0.d0) then
|
||||
do j=1,AO_num
|
||||
do i=1,MO_tot_num
|
||||
do i=1,AO_num
|
||||
X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j)
|
||||
enddo
|
||||
enddo
|
||||
|
@ -14,13 +14,13 @@ default: 15
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the DIIS error vector during a Hartree-Fock calculation
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-5
|
||||
default: 1.e-6
|
||||
|
||||
[thresh_scf]
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the Hartree Fock energy
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-10
|
||||
default: 1.e-12
|
||||
|
||||
[n_it_scf_max]
|
||||
type: Strictly_positive_int
|
||||
@ -32,13 +32,13 @@ default: 128
|
||||
type: Positive_float
|
||||
doc: Energy shift on the virtual MOs to improve SCF convergence
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.5
|
||||
default: 0.4
|
||||
|
||||
[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 [ Damp | DIIS]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: damp
|
||||
default: Damp
|
||||
|
||||
[mo_guess_type]
|
||||
type: MO_guess
|
||||
|
@ -18,57 +18,57 @@
|
||||
END_DOC
|
||||
integer :: i,j,n
|
||||
if (elec_alpha_num == elec_beta_num) then
|
||||
Fock_matrix_mo = Fock_matrix_alpha_mo
|
||||
Fock_matrix_mo = Fock_matrix_mo_alpha
|
||||
else
|
||||
|
||||
do j=1,elec_beta_num
|
||||
! F-K
|
||||
do i=1,elec_beta_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
|
||||
- (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
|
||||
- (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
! F+K/2
|
||||
do i=elec_beta_num+1,elec_alpha_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
! F
|
||||
do i=elec_alpha_num+1, mo_tot_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do j=elec_beta_num+1,elec_alpha_num
|
||||
! F+K/2
|
||||
do i=1,elec_beta_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
! F
|
||||
do i=elec_beta_num+1,elec_alpha_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
|
||||
enddo
|
||||
! F-K/2
|
||||
do i=elec_alpha_num+1, mo_tot_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
|
||||
- 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
|
||||
- 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do j=elec_alpha_num+1, mo_tot_num
|
||||
! F
|
||||
do i=1,elec_beta_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
|
||||
enddo
|
||||
! F-K/2
|
||||
do i=elec_beta_num+1,elec_alpha_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j))&
|
||||
- 0.5d0*(Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
|
||||
- 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
! F+K
|
||||
do i=elec_alpha_num+1,mo_tot_num
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_alpha_mo(i,j)+Fock_matrix_beta_mo(i,j)) &
|
||||
+ (Fock_matrix_beta_mo(i,j) - Fock_matrix_alpha_mo(i,j))
|
||||
Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) &
|
||||
+ (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -81,8 +81,8 @@ END_PROVIDER
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_ao, (ao_num_align, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_beta_ao, (ao_num_align, ao_num) ]
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_ao_alpha, (ao_num_align, ao_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_ao_beta, (ao_num_align, ao_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Alpha Fock matrix in AO basis set
|
||||
@ -92,8 +92,8 @@ END_PROVIDER
|
||||
do j=1,ao_num
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,ao_num
|
||||
Fock_matrix_alpha_ao(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j)
|
||||
Fock_matrix_beta_ao (i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j)
|
||||
Fock_matrix_ao_alpha(i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_alpha(i,j)
|
||||
Fock_matrix_ao_beta (i,j) = ao_mono_elec_integral(i,j) + ao_bi_elec_integral_beta (i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -261,12 +261,7 @@ END_PROVIDER
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_tot_num) ]
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num_align,mo_tot_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Fock matrix on the MO basis
|
||||
@ -275,18 +270,18 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_alpha_mo, (mo_tot_num_align,mo_to
|
||||
allocate ( T(ao_num_align,mo_tot_num) )
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
|
||||
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
|
||||
1.d0, Fock_matrix_alpha_ao,size(Fock_matrix_alpha_ao,1), &
|
||||
1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
|
||||
mo_coef, size(mo_coef,1), &
|
||||
0.d0, T, ao_num_align)
|
||||
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
|
||||
1.d0, mo_coef,size(mo_coef,1), &
|
||||
T, size(T,1), &
|
||||
0.d0, Fock_matrix_alpha_mo, mo_tot_num_align)
|
||||
0.d0, Fock_matrix_mo_alpha, mo_tot_num_align)
|
||||
deallocate(T)
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot_num) ]
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num_align,mo_tot_num) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Fock matrix on the MO basis
|
||||
@ -295,13 +290,13 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_beta_mo, (mo_tot_num_align,mo_tot
|
||||
allocate ( T(ao_num_align,mo_tot_num) )
|
||||
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
|
||||
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
|
||||
1.d0, Fock_matrix_beta_ao,size(Fock_matrix_beta_ao,1), &
|
||||
1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
|
||||
mo_coef, size(mo_coef,1), &
|
||||
0.d0, T, ao_num_align)
|
||||
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
|
||||
1.d0, mo_coef,size(mo_coef,1), &
|
||||
T, size(T,1), &
|
||||
0.d0, Fock_matrix_beta_mo, mo_tot_num_align)
|
||||
0.d0, Fock_matrix_mo_beta, mo_tot_num_align)
|
||||
deallocate(T)
|
||||
END_PROVIDER
|
||||
|
||||
@ -316,8 +311,8 @@ BEGIN_PROVIDER [ double precision, HF_energy ]
|
||||
do j=1,ao_num
|
||||
do i=1,ao_num
|
||||
HF_energy += 0.5d0 * ( &
|
||||
(ao_mono_elec_integral(i,j) + Fock_matrix_alpha_ao(i,j) ) * HF_density_matrix_ao_alpha(i,j) +&
|
||||
(ao_mono_elec_integral(i,j) + Fock_matrix_beta_ao (i,j) ) * HF_density_matrix_ao_beta (i,j) )
|
||||
(ao_mono_elec_integral(i,j) + Fock_matrix_ao_alpha(i,j) ) * HF_density_matrix_ao_alpha(i,j) +&
|
||||
(ao_mono_elec_integral(i,j) + Fock_matrix_ao_beta (i,j) ) * HF_density_matrix_ao_beta (i,j) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -337,7 +332,7 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num_align, ao_num) ]
|
||||
do j=1,ao_num
|
||||
!DIR$ VECTOR ALIGNED
|
||||
do i=1,ao_num_align
|
||||
Fock_matrix_ao(i,j) = Fock_matrix_alpha_ao(i,j)
|
||||
Fock_matrix_ao(i,j) = Fock_matrix_ao_alpha(i,j)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
|
@ -6,16 +6,21 @@ END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
double precision :: energy_SCF,energy_SCF_previous,Delta_energy_SCF,max_error_DIIS
|
||||
double precision, allocatable :: Fock_matrix_DIIS(:,:,:),error_matrix_DIIS(:,:,:)
|
||||
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 (:,:,:)
|
||||
|
||||
integer :: iteration_SCF,dim_DIIS,index_dim_DIIS
|
||||
integer :: dim_DIIS_alpha, dim_DIIS_beta
|
||||
|
||||
integer :: i,j
|
||||
|
||||
allocate( &
|
||||
Fock_matrix_DIIS(AO_num,AO_num,max_dim_DIIS), &
|
||||
error_matrix_DIIS(AO_num,AO_num,max_dim_DIIS) &
|
||||
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) &
|
||||
)
|
||||
|
||||
call write_time(output_hartree_fock)
|
||||
@ -32,6 +37,8 @@ END_DOC
|
||||
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
|
||||
|
||||
@ -50,35 +57,44 @@ END_DOC
|
||||
|
||||
! Store Fock and error matrices at each iteration
|
||||
|
||||
do j=1,AO_num
|
||||
do i=1,AO_num
|
||||
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)
|
||||
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)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Compute the extrapolated Fock matrix
|
||||
|
||||
dim_DIIS_alpha = dim_DIIS
|
||||
call extrapolate_Fock_matrix( &
|
||||
error_matrix_DIIS,Fock_matrix_DIIS, &
|
||||
iteration_SCF,dim_DIIS &
|
||||
error_matrix_DIIS_alpha,Fock_matrix_DIIS_alpha, &
|
||||
Fock_matrix_AO_alpha,size(Fock_matrix_AO_alpha,1), &
|
||||
iteration_SCF,dim_DIIS_alpha &
|
||||
)
|
||||
|
||||
touch Fock_matrix_AO
|
||||
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 &
|
||||
)
|
||||
|
||||
MO_coef = eigenvectors_Fock_matrix_AO
|
||||
dim_DIIS = min(dim_DIIS_alpha,dim_DIIS_beta)
|
||||
touch Fock_matrix_AO_alpha Fock_matrix_AO_beta
|
||||
|
||||
! This algorithm still have an issue with linear dependencies
|
||||
! do i=1,AO_num
|
||||
! write(*,*) i,eigenvalues_Fock_matrix_AO(i)
|
||||
! enddo
|
||||
MO_coef = eigenvectors_Fock_matrix_MO
|
||||
|
||||
touch MO_coef
|
||||
|
||||
! Calculate error vectors
|
||||
|
||||
max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_AO))
|
||||
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)
|
||||
|
||||
! SCF energy
|
||||
|
||||
@ -88,8 +104,8 @@ END_DOC
|
||||
|
||||
! Print results at the end of each iteration
|
||||
|
||||
write(output_hartree_fock,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10)') &
|
||||
iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS
|
||||
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
|
||||
|
||||
enddo
|
||||
|
||||
@ -114,6 +130,7 @@ end
|
||||
|
||||
subroutine extrapolate_Fock_matrix( &
|
||||
error_matrix_DIIS,Fock_matrix_DIIS, &
|
||||
Fock_matrix_AO_,size_Fock_matrix_AO, &
|
||||
iteration_SCF,dim_DIIS &
|
||||
)
|
||||
|
||||
@ -123,11 +140,13 @@ END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
double precision,intent(in) :: Fock_matrix_DIIS(AO_num,AO_num,*),error_matrix_DIIS(AO_num,AO_num,*)
|
||||
integer,intent(in) :: iteration_SCF
|
||||
double precision,intent(in) :: Fock_matrix_DIIS(ao_num,ao_num,*),error_matrix_DIIS(ao_num,ao_num,*)
|
||||
integer,intent(in) :: iteration_SCF, size_Fock_matrix_AO
|
||||
double precision,intent(inout):: Fock_matrix_AO_(size_Fock_matrix_AO,ao_num)
|
||||
integer,intent(inout) :: dim_DIIS
|
||||
|
||||
double precision,allocatable :: B_matrix_DIIS(:,:),X_vector_DIIS(:)
|
||||
double precision,allocatable :: C_vector_DIIS(:)
|
||||
|
||||
double precision,allocatable :: scratch(:,:)
|
||||
integer :: i,j,k,i_DIIS,j_DIIS
|
||||
@ -135,7 +154,8 @@ END_DOC
|
||||
allocate( &
|
||||
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), &
|
||||
X_vector_DIIS(dim_DIIS+1), &
|
||||
scratch(AO_num,AO_num) &
|
||||
C_vector_DIIS(dim_DIIS+1), &
|
||||
scratch(ao_num,ao_num) &
|
||||
)
|
||||
|
||||
! Compute the matrices B and X
|
||||
@ -147,7 +167,7 @@ END_DOC
|
||||
|
||||
! Compute product of two errors vectors
|
||||
|
||||
call dgemm('N','N',AO_num,AO_num,AO_num, &
|
||||
call dgemm('N','N',ao_num,ao_num,ao_num, &
|
||||
1.d0, &
|
||||
error_matrix_DIIS(1,1,i_DIIS),size(error_matrix_DIIS,1), &
|
||||
error_matrix_DIIS(1,1,j_DIIS),size(error_matrix_DIIS,1), &
|
||||
@ -157,7 +177,7 @@ END_DOC
|
||||
! Compute Trace
|
||||
|
||||
B_matrix_DIIS(i,j) = 0.d0
|
||||
do k=1,AO_num
|
||||
do k=1,ao_num
|
||||
B_matrix_DIIS(i,j) += scratch(k,k)
|
||||
enddo
|
||||
enddo
|
||||
@ -168,10 +188,10 @@ END_DOC
|
||||
do i=1,dim_DIIS
|
||||
B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
|
||||
B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
|
||||
X_vector_DIIS(i) = 0.d0
|
||||
C_vector_DIIS(i) = 0.d0
|
||||
enddo
|
||||
B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1) = 0.d0
|
||||
X_vector_DIIS(dim_DIIS+1) = -1.d0
|
||||
C_vector_DIIS(dim_DIIS+1) = -1.d0
|
||||
|
||||
! Solve the linear system C = B.X
|
||||
|
||||
@ -182,24 +202,39 @@ END_DOC
|
||||
ipiv(dim_DIIS+1) &
|
||||
)
|
||||
|
||||
call dsysv('U',dim_DIIS+1,1, &
|
||||
double precision, allocatable :: AF(:,:)
|
||||
allocate (AF(dim_DIIS+1,dim_DIIS+1))
|
||||
double precision :: rcond, ferr, berr
|
||||
integer :: iwork(dim_DIIS+1)
|
||||
|
||||
call dsysvx('N','U',dim_DIIS+1,1, &
|
||||
B_matrix_DIIS,size(B_matrix_DIIS,1), &
|
||||
AF, size(AF,1), &
|
||||
ipiv, &
|
||||
C_vector_DIIS,size(C_vector_DIIS,1), &
|
||||
X_vector_DIIS,size(X_vector_DIIS,1), &
|
||||
rcond, &
|
||||
ferr, &
|
||||
berr, &
|
||||
scratch,size(scratch), &
|
||||
iwork, &
|
||||
info &
|
||||
)
|
||||
|
||||
if(info == 0) then
|
||||
if(info < 0) then
|
||||
stop 'bug in DIIS'
|
||||
endif
|
||||
|
||||
if (rcond > 1.d-8) then
|
||||
|
||||
! Compute extrapolated Fock matrix
|
||||
|
||||
Fock_matrix_AO(:,:) = 0.d0
|
||||
Fock_matrix_AO_(:,:) = 0.d0
|
||||
|
||||
do k=1,dim_DIIS
|
||||
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)
|
||||
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)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -209,9 +244,9 @@ END_DOC
|
||||
dim_DIIS = 0
|
||||
endif
|
||||
|
||||
! do i=1,AO_num
|
||||
! do j=1,AO_num
|
||||
! write(*,*) Fock_matrix_AO(i,j)
|
||||
! do i=1,ao_num
|
||||
! do j=1,ao_num
|
||||
! write(*,*) Fock_matrix_AO_(i,j)
|
||||
! enddo
|
||||
! enddo
|
||||
|
||||
|
@ -114,7 +114,6 @@ subroutine damping_SCF
|
||||
mo_coef = eigenvectors_fock_matrix_mo
|
||||
TOUCH mo_coef
|
||||
|
||||
|
||||
enddo
|
||||
write(output_hartree_fock,'(A4,1X,A16, 1X, A16, 1X, A16, 1X, A4 )') '====','================','================','================', '===='
|
||||
write(output_hartree_fock,*)
|
||||
|
@ -22,7 +22,7 @@ subroutine huckel_guess
|
||||
Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
|
||||
ao_mono_elec_integral_diag(j))
|
||||
enddo
|
||||
Fock_matrix_ao(j,j) = Fock_matrix_alpha_ao(j,j)
|
||||
Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j)
|
||||
enddo
|
||||
TOUCH Fock_matrix_ao
|
||||
mo_coef = eigenvectors_fock_matrix_mo
|
||||
|
@ -139,8 +139,6 @@ BEGIN_PROVIDER [ double precision, mo_occ, (mo_tot_num) ]
|
||||
endif
|
||||
END_PROVIDER
|
||||
|
||||
|
||||
|
||||
subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
|
@ -200,7 +200,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
|
||||
!
|
||||
! LDC : leftmost dimension of C
|
||||
!
|
||||
! m : Coefficients matrix is MxN, ( array is (LDC,N) )
|
||||
! M : Coefficients matrix is MxN, ( array is (LDC,N) )
|
||||
!
|
||||
END_DOC
|
||||
|
||||
@ -211,7 +211,6 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
|
||||
double precision, allocatable :: Vt(:,:)
|
||||
double precision, allocatable :: D(:)
|
||||
double precision, allocatable :: S_half(:,:)
|
||||
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
|
||||
integer :: info, i, j, k
|
||||
|
||||
if (n < 2) then
|
||||
|
Loading…
Reference in New Issue
Block a user