LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-11-06 22:24:00 +01:00
Code Issues Releases Wiki Activity

Fixed

Browse Source
This commit is contained in:
Anthony Scemama 2017-09-13 09:06:32 +02:00
parent 485ffb4bef
commit eed7cc8c14
25 changed files with 175 additions and 281 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
plugins/FOBOCI/diag_fock_inactiv_virt.irp.f
View File

@ -29,7 +29,7 @@ subroutine diag_inactive_virt_and_update_mos
label = "Canonical"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.)
soft_touch mo_coef
@ -76,7 +76,7 @@ subroutine diag_inactive_virt_new_and_update_mos
label = "Canonical"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.false.)
soft_touch mo_coef

2
plugins/FOBOCI/new_approach.irp.f
View File

@ -655,7 +655,7 @@ subroutine new_approach
integer :: sign
sign = -1
call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign)
call mo_as_eigvectors_of_mo_matrix(one_body_dm_mo,size(one_body_dm_mo,1),size(one_body_dm_mo,2),label,sign,.true.)
soft_touch mo_coef
call save_mos

4
plugins/FOBOCI/routines_foboci.irp.f
View File

@ -449,7 +449,7 @@ subroutine save_osoci_natural_mos
label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.)
!if(disk_access_ao_integrals == "None" .or. disk_access_ao_integrals == "Write" )then
! disk_access_ao_integrals = "Read"
! touch disk_access_ao_integrals
@ -584,7 +584,7 @@ subroutine set_osoci_natural_mos
enddo
label = "Natural"
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1)
call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label,1,.true.)
soft_touch mo_coef
deallocate(tmp,occ)

63
plugins/Hartree_Fock/DIIS.irp.f
View File

@ -94,7 +94,7 @@ END_PROVIDER
do i=1,AO_num
do j=1,AO_num
Xt(i,j) = X_Matrix_AO(j,i)
Xt(i,j) = S_half_inv(j,i)
enddo
enddo
@ -103,7 +103,7 @@ END_PROVIDER
call dgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, &
Fock_matrix_AO,size(Fock_matrix_AO,1), &
X_Matrix_AO,size(X_Matrix_AO,1), &
S_half_inv,size(S_half_inv,1), &
0.d0, &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))
@ -130,67 +130,10 @@ END_PROVIDER
call dgemm('N','N',AO_num,AO_num,AO_num, &
1.d0, &
X_matrix_AO,size(X_matrix_AO,1), &
S_half_inv,size(S_half_inv,1), &
scratch,size(scratch,1), &
0.d0, &
eigenvectors_Fock_matrix_AO,size(eigenvectors_Fock_matrix_AO,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, X_matrix_AO, (AO_num,AO_num) ]
BEGIN_DOC
! Matrix X = S^{-1/2} obtained by SVD
END_DOC
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(X_matrix_AO,1)
allocate( &
U(LDC,AO_num), &
Vt(LDA,AO_num), &
D(AO_num))
call svd( &
AO_overlap,LDA, &
U,LDC, &
D, &
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
ASSERT (D(i) > 0.d0)
D(i) = 1.d0/sqrt(D(i))
endif
do j=1,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,AO_num
X_matrix_AO(i,j) = X_matrix_AO(i,j) + U(i,k)*D(k)*Vt(k,j)
enddo
enddo
endif
enddo
END_PROVIDER

118
plugins/Hartree_Fock/Fock_matrix.irp.f
View File

@ -263,17 +263,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_alpha, (mo_tot_num,mo_tot_num)
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate ( T(ao_num,mo_tot_num) )
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, size(T,1))
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_mo_alpha, size(Fock_matrix_mo_alpha,1))
deallocate(T)
call ao_to_mo(Fock_matrix_ao_alpha,size(Fock_matrix_ao_alpha,1), &
Fock_matrix_mo_alpha,size(Fock_matrix_mo_alpha,1))
END_PROVIDER
@ -282,17 +273,8 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_mo_beta, (mo_tot_num,mo_tot_num)
BEGIN_DOC
! Fock matrix on the MO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate ( T(ao_num,mo_tot_num) )
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
mo_coef, size(mo_coef,1), &
0.d0, T, size(T,1))
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_mo_beta, size(Fock_matrix_mo_beta,1))
deallocate(T)
call ao_to_mo(Fock_matrix_ao_beta,size(Fock_matrix_ao_beta,1), &
Fock_matrix_mo_beta,size(Fock_matrix_mo_beta,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_energy ]
@ -330,97 +312,9 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_ao, (ao_num, ao_num) ]
enddo
enddo
else
double precision, allocatable :: T(:,:), M(:,:)
integer :: ierr
! F_ao = S C F_mo C^t S
allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
if (ierr /=0 ) then
print *, irp_here, ' : allocation failed'
endif
! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
! -> M(ao_num,mo_tot_num)
call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, &
ao_overlap, size(ao_overlap,1), &
mo_coef, size(mo_coef,1), &
0.d0, &
M, size(M,1))
! M(ao_num,mo_tot_num) . Fock_matrix_mo (mo_tot_num,mo_tot_num)
! -> T(ao_num,mo_tot_num)
call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, &
M, size(M,1), &
Fock_matrix_mo, size(Fock_matrix_mo,1), &
0.d0, &
T, size(T,1))
! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
! -> M(ao_num,ao_num)
call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, &
T, size(T,1), &
mo_coef, size(mo_coef,1), &
0.d0, &
M, size(M,1))
! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
! -> Fock_matrix_ao(ao_num,ao_num)
call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, &
M, size(M,1), &
ao_overlap, size(ao_overlap,1), &
0.d0, &
Fock_matrix_ao, size(Fock_matrix_ao,1))
deallocate(T)
call mo_to_ao(Fock_matrix_mo,size(Fock_matrix_mo,1), &
Fock_matrix_ao,size(Fock_matrix_ao,1))
endif
END_PROVIDER
subroutine Fock_mo_to_ao(FMO,LDFMO,FAO,LDFAO)
implicit none
integer, intent(in) :: LDFMO ! size(FMO,1)
integer, intent(in) :: LDFAO ! size(FAO,1)
double precision, intent(in) :: FMO(LDFMO,*)
double precision, intent(out) :: FAO(LDFAO,*)
double precision, allocatable :: T(:,:), M(:,:)
integer :: ierr
! F_ao = S C F_mo C^t S
allocate (T(ao_num,ao_num),M(ao_num,ao_num),stat=ierr)
if (ierr /=0 ) then
print *, irp_here, ' : allocation failed'
endif
! ao_overlap (ao_num,ao_num) . mo_coef (ao_num,mo_tot_num)
! -> M(ao_num,mo_tot_num)
call dgemm('N','N', ao_num,mo_tot_num,ao_num, 1.d0, &
ao_overlap, size(ao_overlap,1), &
mo_coef, size(mo_coef,1), &
0.d0, &
M, size(M,1))
! M(ao_num,mo_tot_num) . FMO (mo_tot_num,mo_tot_num)
! -> T(ao_num,mo_tot_num)
call dgemm('N','N', ao_num,mo_tot_num,mo_tot_num, 1.d0, &
M, size(M,1), &
FMO, size(FMO,1), &
0.d0, &
T, size(T,1))
! T(ao_num,mo_tot_num) . mo_coef^T (mo_tot_num,ao_num)
! -> M(ao_num,ao_num)
call dgemm('N','T', ao_num,ao_num,mo_tot_num, 1.d0, &
T, size(T,1), &
mo_coef, size(mo_coef,1), &
0.d0, &
M, size(M,1))
! M(ao_num,ao_num) . ao_overlap (ao_num,ao_num)
! -> Fock_matrix_ao(ao_num,ao_num)
call dgemm('N','N', ao_num,ao_num,ao_num, 1.d0, &
M, size(M,1), &
ao_overlap, size(ao_overlap,1), &
0.d0, &
FAO, size(FAO,1))
deallocate(T,M)
end

3
plugins/Hartree_Fock/Roothaan_Hall_SCF.irp.f
View File

@ -83,7 +83,6 @@ END_DOC
! Calculate error vectors
max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_MO))
! max_error_DIIS = maxval(Abs(FPS_SPF_Matrix_AO))
! SCF energy
@ -134,7 +133,7 @@ END_DOC
write(output_hartree_fock,*)
if(.not.no_oa_or_av_opt)then
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
endif
call write_double(output_hartree_fock, Energy_SCF, 'Hartree-Fock energy')

2
plugins/Hartree_Fock/SCF.irp.f
View File

@ -23,7 +23,7 @@ subroutine create_guess
mo_coef = ao_ortho_lowdin_coef
TOUCH mo_coef
mo_label = 'Guess'
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label)
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,size(mo_mono_elec_integral,1),size(mo_mono_elec_integral,2),mo_label,.false.)
SOFT_TOUCH mo_coef mo_label
else if (mo_guess_type == "Huckel") then
call huckel_guess

2
plugins/Hartree_Fock/damping_SCF.irp.f
View File

@ -119,7 +119,7 @@ subroutine damping_SCF
write(output_hartree_fock,*)
if(.not.no_oa_or_av_opt)then
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1)
call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1),size(Fock_matrix_mo,2),mo_label,1,.true.)
endif
call write_double(output_hartree_fock, E_min, 'Hartree-Fock energy')

13
plugins/Hartree_Fock/huckel.irp.f
View File

@ -11,21 +11,24 @@ subroutine huckel_guess
label = "Guess"
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, &
size(mo_mono_elec_integral,1), &
size(mo_mono_elec_integral,2),label,1)
size(mo_mono_elec_integral,2),label,1,.false.)
TOUCH mo_coef
c = 0.5d0 * 1.75d0
do j=1,ao_num
do i=1,ao_num
Fock_matrix_ao(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
ao_mono_elec_integral_diag(j))
Fock_matrix_ao_alpha(i,j) = c*ao_overlap(i,j)*(ao_mono_elec_integral_diag(i) + &
ao_mono_elec_integral_diag(j))
Fock_matrix_ao_beta (i,j) = Fock_matrix_ao_alpha(i,j)
enddo
Fock_matrix_ao(j,j) = Fock_matrix_ao_alpha(j,j)
Fock_matrix_ao_alpha(j,j) = ao_mono_elec_integral(j,j) + ao_bi_elec_integral_alpha(j,j)
Fock_matrix_ao_beta (j,j) = Fock_matrix_ao_alpha(j,j)
enddo
TOUCH Fock_matrix_ao
TOUCH Fock_matrix_ao_alpha Fock_matrix_ao_beta
mo_coef = eigenvectors_fock_matrix_mo
SOFT_TOUCH mo_coef
call save_mos
print *, 'E=', HF_energy
end

71
src/AO_Basis/ao_overlap.irp.f
View File

@ -1,7 +1,7 @@
BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num,ao_num) ]
BEGIN_PROVIDER [ double precision, ao_overlap,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_x,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_y,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z,(ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between atomic basis functions:
@ -69,7 +69,7 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between absolute value of atomic basis functions:
@ -123,11 +123,68 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_inv,(ao_num,ao_num) ]
BEGIN_PROVIDER [ double precision, S_inv,(ao_num,ao_num) ]
implicit none
BEGIN_DOC
! S^-1
END_DOC
call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,ao_overlap_inv,size(ao_overlap_inv,1))
call get_pseudo_inverse(ao_overlap,size(ao_overlap,1),ao_num,ao_num,S_inv,size(S_inv,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, S_half_inv, (AO_num,AO_num) ]
BEGIN_DOC
! Matrix X = S^{-1/2} obtained by SVD
END_DOC
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(S_half_inv,1)
allocate( &
U(LDC,AO_num), &
Vt(LDA,AO_num), &
D(AO_num))
call svd( &
AO_overlap,LDA, &
U,LDC, &
D, &
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
ASSERT (D(i) > 0.d0)
D(i) = 1.d0/sqrt(D(i))
endif
do j=1,AO_num
S_half_inv(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,AO_num
S_half_inv(i,j) = S_half_inv(i,j) + U(i,k)*D(k)*Vt(k,j)
enddo
enddo
endif
enddo
END_PROVIDER

13
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -337,8 +337,8 @@ end
END_DOC
mo_integrals_cache_min_8 = max(1_8,elec_alpha_num - 63_8)
mo_integrals_cache_max_8 = min(int(mo_tot_num,8),mo_integrals_cache_min+127_8)
mo_integrals_cache_min = mo_integrals_cache_min_8
mo_integrals_cache_max = mo_integrals_cache_max_8
mo_integrals_cache_min = max(1,elec_alpha_num - 63)
mo_integrals_cache_max = min(mo_tot_num,mo_integrals_cache_min+127)
END_PROVIDER
@ -349,17 +349,22 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:128_8*128_8*128_8*12
END_DOC
PROVIDE mo_bielec_integrals_in_map
integer*8 :: i,j,k,l
integer*8 :: i4,j4,k4,l4
integer*8 :: ii
integer(key_kind) :: idx
real(integral_kind) :: integral
FREE ao_integrals_cache
!$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
!$OMP PARALLEL DO PRIVATE (i,j,k,l,i4,j4,k4,l4,idx,ii,integral)
do l=mo_integrals_cache_min_8,mo_integrals_cache_max_8
do k=mo_integrals_cache_min_8,mo_integrals_cache_max_8
do j=mo_integrals_cache_min_8,mo_integrals_cache_max_8
do i=mo_integrals_cache_min_8,mo_integrals_cache_max_8
i4 = int(i,4)
j4 = int(j,4)
k4 = int(k,4)
l4 = int(l,4)
!DIR$ FORCEINLINE
call bielec_integrals_index(i,j,k,l,idx)
call bielec_integrals_index(i4,j4,k4,l4,idx)
!DIR$ FORCEINLINE
call map_get(mo_integrals_map,idx,integral)
ii = l-mo_integrals_cache_min_8

2
src/Integrals_Monoelec/ao_mono_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_diag,(ao_num)]
implicit none
integer :: i,j,n,l

13
src/Integrals_Monoelec/kin_ao_ints.irp.f
View File

@ -1,6 +1,6 @@
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num,ao_num) ]
BEGIN_PROVIDER [ double precision, ao_deriv2_x,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_y,(ao_num_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_z,(ao_num_align,ao_num) ]
implicit none
integer :: i,j,n,l
double precision :: f
@ -119,7 +119,7 @@
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)]
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! array of the priminitve basis kinetic integrals
@ -134,11 +134,14 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)]
else
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(ao_num, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
!$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
do j = 1, ao_num
do i = 1, ao_num
ao_kinetic_integral(i,j) = -0.5d0 * (ao_deriv2_x(i,j) + ao_deriv2_y(i,j) + ao_deriv2_z(i,j) )
enddo
do i = ao_num +1,ao_num_align
ao_kinetic_integral(i,j) = 0.d0
enddo
enddo
!$OMP END PARALLEL DO
endif

2
src/Integrals_Monoelec/kin_mo_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! Kinetic energy integrals in the MO basis

2
src/Integrals_Monoelec/mo_mono_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l
BEGIN_DOC

4
src/Integrals_Monoelec/pot_ao_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
BEGIN_DOC
! interaction nuclear electron
END_DOC
@ -80,7 +80,7 @@ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num,ao_num,nucl_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)]
BEGIN_DOC
! ao_nucl_elec_integral_per_atom(i,j,k) = -<AO(i)|1/|r-Rk|AO(j)>
! where Rk is the geometry of the kth atom

6
src/Integrals_Monoelec/pot_ao_pseudo_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! Pseudo-potential integrals
@ -29,7 +29,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! Local pseudo-potential
@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! Local pseudo-potential

4
src/Integrals_Monoelec/pot_mo_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! interaction nuclear electron on the MO basis
@ -25,7 +25,7 @@ BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num,mo_tot_num)
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num,mo_tot_num,nucl_num)]
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num_align,mo_tot_num,nucl_num)]
implicit none
BEGIN_DOC
! mo_nucl_elec_integral_per_atom(i,j,k) = -<MO(i)|1/|r-Rk|MO(j)>

2
src/Integrals_Monoelec/pot_mo_pseudo_ints.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num_align,mo_tot_num)]
implicit none
BEGIN_DOC
! interaction nuclear electron on the MO basis

18
src/Integrals_Monoelec/spread_dipole_ao.irp.f
View File

@ -1,6 +1,6 @@
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_spread_x, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_y, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_spread_z, (ao_num_align,ao_num)]
BEGIN_DOC
! array of the integrals of AO_i * x^2 AO_j
! array of the integrals of AO_i * y^2 AO_j
@ -69,9 +69,9 @@
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_dipole_x, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_y, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_dipole_z, (ao_num_align,ao_num)]
BEGIN_DOC
! array of the integrals of AO_i * x AO_j
! array of the integrals of AO_i * y AO_j
@ -139,9 +139,9 @@
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_deriv_1_x, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_y, (ao_num_align,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_deriv_1_z, (ao_num_align,ao_num)]
BEGIN_DOC
! array of the integrals of AO_i * d/dx AO_j
! array of the integrals of AO_i * d/dy AO_j

12
src/Integrals_Monoelec/spread_dipole_mo.irp.f
View File

@ -1,6 +1,6 @@
BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_dipole_x , (mo_tot_num_align,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_dipole_y , (mo_tot_num_align,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_dipole_z , (mo_tot_num_align,mo_tot_num)]
BEGIN_DOC
! array of the integrals of MO_i * x MO_j
! array of the integrals of MO_i * y MO_j
@ -29,9 +29,9 @@
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_spread_x , (mo_tot_num_align,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_spread_y , (mo_tot_num_align,mo_tot_num)]
&BEGIN_PROVIDER [double precision, mo_spread_z , (mo_tot_num_align,mo_tot_num)]
BEGIN_DOC
! array of the integrals of MO_i * x^2 MO_j
! array of the integrals of MO_i * y^2 MO_j

8
src/MOGuess/H_CORE_guess.irp.f
View File

@ -5,11 +5,5 @@ program H_CORE_guess
END_DOC
implicit none
character*(64) :: label
label = "Guess"
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, &
size(mo_mono_elec_integral,1), &
size(mo_mono_elec_integral,2),label,1)
print *, 'save mos'
call save_mos
call h_core_guess
end

3
src/MOGuess/h_core_guess_routine.irp.f
View File

@ -7,8 +7,7 @@ subroutine hcore_guess
label = "Guess"
call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral, &
size(mo_mono_elec_integral,1), &
size(mo_mono_elec_integral,2),label,1)
print *, 'save mos'
size(mo_mono_elec_integral,2),label,1,.false.)
call save_mos
SOFT_TOUCH mo_coef mo_label
end

56
src/MO_Basis/mos.irp.f
View File

@ -26,7 +26,7 @@ BEGIN_PROVIDER [ integer, mo_tot_num_align ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ]
BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Molecular orbital coefficients on AO basis set
@ -35,7 +35,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ]
END_DOC
integer :: i, j
double precision, allocatable :: buffer(:,:)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: buffer
logical :: exists
PROVIDE ezfio_filename
@ -50,9 +49,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ]
do j=1,ao_num
mo_coef(j,i) = buffer(j,i)
enddo
do j=ao_num+1,ao_num_align
mo_coef(j,i) = 0.d0
enddo
enddo
deallocate(buffer)
else
@ -61,9 +57,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num_align,mo_tot_num) ]
do j=1,ao_num
mo_coef(j,i) = ao_ortho_canonical_coef(j,i)
enddo
do j=ao_num+1,ao_num_align
mo_coef(j,i) = 0.d0
enddo
enddo
endif
END_PROVIDER
@ -125,7 +118,7 @@ BEGIN_PROVIDER [ double precision, S_mo_coef, (ao_num_align, mo_tot_num) ]
END_DOC
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap,size(ao_overlap,1), &
1.d0, ao_overlap,size(ao_overlap,1), &
mo_coef, size(mo_coef,1), &
0.d0, S_mo_coef, size(S_mo_coef,1))
@ -159,24 +152,25 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
BEGIN_DOC
! Transform A from the AO basis to the MO basis
!
! (S.C)t.A_ao.S.C
! Ct.A_ao.C
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_ao(LDA_ao,ao_num)
double precision, intent(out) :: A_mo(LDA_mo,mo_tot_num)
double precision, allocatable :: T(:,:)
allocate ( T(mo_tot_num,ao_num) )
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('T','N', mo_tot_num, ao_num, ao_num, &
1.d0, S_mo_coef, size(S_mo_coef,1), &
A_ao, size(A_ao,1), &
0.d0, T, size(T,1))
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, A_ao,LDA_ao, &
mo_coef, size(mo_coef,1), &
0.d0, T, ao_num_align)
call dgemm('N','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, T,size(T,1), &
S_mo_coef, size(S_mo_coef,1), &
0.d0, A_mo, size(A_mo,1))
call dgemm('T','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, ao_num_align, &
0.d0, A_mo, LDA_mo)
deallocate(T)
end
@ -186,7 +180,7 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
BEGIN_DOC
! Transform A from the MO basis to the AO basis
!
! C.A_mo.Ct
! (S.C).A_mo.(S.C)t
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo,mo_tot_num)
@ -194,15 +188,14 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
double precision, allocatable :: T(:,:)
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
mo_coef, size(mo_coef,1), &
1.d0, A_mo,size(A_mo,1), &
S_mo_coef, size(S_mo_coef,1), &
0.d0, T, size(T,1))
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, mo_coef,size(mo_coef,1), &
1.d0, S_mo_coef, size(S_mo_coef,1), &
T, size(T,1), &
0.d0, A_ao, size(A_ao,1))
@ -257,18 +250,17 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
double precision, allocatable :: T(:,:)
allocate ( T(ao_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('T','N', ao_num, ao_num, ao_num, &
call dgemm('T','N', ao_num, ao_num, ao_num, &
1.d0, &
ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1), &
A_ao,LDA_ao, &
0.d0, T, ao_num_align)
ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),&
A_ao,size(A_ao,1), &
0.d0, T, size(T,1))
call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, &
T, size(T,1), &
call dgemm('N','N', ao_num, ao_num, ao_num, 1.d0, &
T, size(T,1), &
ao_ortho_canonical_coef_inv,size(ao_ortho_canonical_coef_inv,1),&
0.d0, A, LDA)
0.d0, A, size(A,1))
deallocate(T)
end

29
src/MO_Basis/utils.irp.f
View File

@ -44,11 +44,12 @@ subroutine save_mos_truncated(n)
end
subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign,output)
implicit none
integer,intent(in) :: n,m, sign
character*(64), intent(in) :: label
double precision, intent(in) :: matrix(n,m)
logical, intent(in) :: output
integer :: i,j
double precision, allocatable :: mo_coef_new(:,:), R(:,:),eigvalues(:), A(:,:)
@ -76,22 +77,26 @@ subroutine mo_as_eigvectors_of_mo_matrix(matrix,n,m,label,sign)
mo_coef_new = mo_coef
call lapack_diag(eigvalues,R,A,n,m)
write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**'
write (output_mo_basis,'(A)') ''
write (output_mo_basis,'(A)') 'Eigenvalues'
write (output_mo_basis,'(A)') '-----------'
write (output_mo_basis,'(A)') ''
write (output_mo_basis,'(A)') '======== ================'
if (output) then
write (output_mo_basis,'(A)') 'MOs are now **'//trim(label)//'**'
write (output_mo_basis,'(A)') ''
write (output_mo_basis,'(A)') 'Eigenvalues'
write (output_mo_basis,'(A)') '-----------'
write (output_mo_basis,'(A)') ''
write (output_mo_basis,'(A)') '======== ================'
endif
if (sign == -1) then
do i=1,m
eigvalues(i) = -eigvalues(i)
enddo
endif
do i=1,m
write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i)
enddo
write (output_mo_basis,'(A)') '======== ================'
write (output_mo_basis,'(A)') ''
if (output) then
do i=1,m
write (output_mo_basis,'(I8,1X,F16.10)') i,eigvalues(i)
enddo
write (output_mo_basis,'(A)') '======== ================'
write (output_mo_basis,'(A)') ''
endif
call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),R,size(R,1),0.d0,mo_coef,size(mo_coef,1))
deallocate(A,mo_coef_new,R,eigvalues)