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

Merge branch 'master' of github.com:LCPQ/quantum_package into LCPQ-master

Browse Source
This commit is contained in:
Anthony Scemama 2017-06-26 19:59:37 +02:00
commit 78741169d8
14 changed files with 462 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
ocaml/qp_run.ml
View File

@ -134,6 +134,7 @@ let run slave exe ezfio_file =
let duration = Time.diff (Time.now()) time_start
|> Core.Span.to_string in
Printf.printf "Wall time : %s\n\n" duration;
if (exit_code <> 0) then
exit exit_code
let spec =

12
plugins/Hartree_Fock_SlaterDressed/EZFIO.cfg
View File

@ -10,5 +10,17 @@ doc: Exponents of the additional Slater functions
size: (nuclei.nucl_num,mo_basis.mo_tot_num)
interface: ezfio, provider
[projector]
type: double precision
doc: Orthogonal AO basis
size: (ao_basis.ao_num,ao_basis.ao_num)
interface: ezfio
[ao_orthoSlaOverlap]
type: double precision
doc: Orthogonal AO basis
size: (ao_basis.ao_num,nuclei.nucl_num)
interface: ezfio

20
plugins/Hartree_Fock_SlaterDressed/LinearSystem.irp.f
View File

@ -13,7 +13,7 @@ BEGIN_PROVIDER [ double precision, cusp_A, (nucl_num, nucl_num) ]
cusp_A(A,B) -= slater_value_at_nucl(B,A)
! Projector
do mu=1,mo_tot_num
cusp_A(A,B) += MOSlaOverlap_matrix(mu,B) * mo_value_at_nucl(mu,A)
cusp_A(A,B) += AO_orthoSlaOverlap_matrix(mu,B) * ao_ortho_value_at_nucl(mu,A)
enddo
enddo
enddo
@ -41,13 +41,23 @@ BEGIN_PROVIDER [ double precision, cusp_C, (nucl_num, mo_tot_num) ]
! Equations to solve : A.C = B
END_DOC
double precision, allocatable :: AF(:,:)
integer :: info
integer :: ipiv(nucl_num)
double precision, allocatable :: AF(:,:)
allocate ( AF(nucl_num,nucl_num) )
call get_pseudo_inverse(cusp_A,nucl_num,nucl_num,AF,size(AF,1))
call dgemm('N','N',nucl_num,mo_tot_num,nucl_num,1.d0, &
AF,size(AF,1), cusp_B, size(cusp_B,1), 0.d0, cusp_C, size(cusp_C,1))
cusp_C(1:nucl_num,1:mo_tot_num) = cusp_B(1:nucl_num,1:mo_tot_num)
AF(1:nucl_num,1:nucl_num) = cusp_A(1:nucl_num,1:nucl_num)
call dgetrf(nucl_num,nucl_num,AF,size(AF,1),ipiv,info)
if (info /= 0) then
print *, info
stop 'dgetrf failed'
endif
call dgetrs('N',nucl_num,mo_tot_num,AF,size(AF,1),ipiv,cusp_C,size(cusp_C,1),info)
if (info /= 0) then
print *, info
stop 'dgetrs failed'
endif
END_PROVIDER

21
plugins/Hartree_Fock_SlaterDressed/SCF_dressed.irp.f
View File

@ -76,21 +76,28 @@ subroutine run
double precision :: EHF
integer :: i_it, i, j, k
mo_label = "CuspDressed"
mo_label = 'None'
print *, HF_energy
call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(cusp_C)
! print *, HF_energy
do i=1,ao_num
print *, mo_coef(i,1), cusp_corrected_mos(i,1)
enddo
mo_coef(1:ao_num,1:mo_tot_num) = cusp_corrected_mos(1:ao_num,1:mo_tot_num)
TOUCH mo_coef
SOFT_TOUCH mo_coef slater_coef
call ezfio_set_Hartree_Fock_SlaterDressed_slater_coef_ezfio(slater_coef)
call ezfio_set_Hartree_Fock_SlaterDressed_projector(ao_ortho_canonical_coef(1:ao_num,1:ao_num))
call ezfio_set_Hartree_Fock_SlaterDressed_ao_orthoSlaOverlap(AO_orthoSlaOverlap_matrix)
call save_mos
print *, 'ci'
print *, mo_coef(1:ao_num,1)
print *, 'cAi'
print *, slater_coef
EHF = HF_energy
print *, HF_energy
! EHF = HF_energy
! print *, HF_energy
! call Roothaan_Hall_SCF
end

20
plugins/Hartree_Fock_SlaterDressed/at_nucl.irp.f
View File

@ -28,14 +28,26 @@ BEGIN_PROVIDER [ double precision , ao_value_at_nucl, (ao_num,nucl_num) ]
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_value_at_nucl, (ao_num,nucl_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals at the nucleus
END_DOC
call dgemm('T','N',ao_num,nucl_num,ao_num,1.d0, &
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
ao_value_at_nucl, size(ao_value_at_nucl,1), &
0.d0, ao_ortho_value_at_nucl,size(ao_ortho_value_at_nucl,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_value_at_nucl, (mo_tot_num,nucl_num) ]
implicit none
BEGIN_DOC
! Values of the molecular orbitals at the nucleus
END_DOC
call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, &
mo_coef_transp, size(mo_coef_transp,1), &
call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0, &
mo_coef, size(mo_coef,1), &
ao_value_at_nucl, size(ao_value_at_nucl,1), &
0.d0, mo_value_at_nucl, size(mo_value_at_nucl,1))
END_PROVIDER
@ -54,10 +66,6 @@ BEGIN_PROVIDER [ double precision , slater_value_at_nucl, (nucl_num,nucl_num) ]
x = nucl_coord(i,1) - nucl_coord(k,1)
y = nucl_coord(i,2) - nucl_coord(k,2)
z = nucl_coord(i,3) - nucl_coord(k,3)
! expo = slater_expo(i)*slater_expo(i)*((x*x) + (y*y) + (z*z))
! if (expo > 160.d0) cycle
! expo = dsqrt(expo)
expo = slater_expo(i) * dsqrt((x*x) + (y*y) + (z*z))
slater_value_at_nucl(i,k) = dexp(-expo) * slater_normalization(i)
enddo

172
plugins/Hartree_Fock_SlaterDressed/dressing.irp.f Normal file
View File

@ -0,0 +1,172 @@
BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral_dressing, (ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the orthogonal AO basis set
END_DOC
integer :: i,j,k
integer :: mu, nu, lambda, A
double precision :: tmp
ao_ortho_mono_elec_integral_dressing = 0.d0
i = idx_dressing
do mu=1,ao_num
if (dabs(mo_coef_in_ao_ortho_basis(mu,i)) > 1.d-5) then
do A=1,nucl_num
tmp = 0.d0
do nu=1,ao_num
tmp += AO_orthoSlaOverlap_matrix(nu,A) * ao_ortho_mono_elec_integral(mu,nu)
enddo
ao_ortho_mono_elec_integral_dressing(mu,mu) += cusp_C(A,i) * (AO_orthoSlaH_matrix(mu,A) - tmp)
enddo
ao_ortho_mono_elec_integral_dressing(mu,mu) *= 1.d0/mo_coef_in_ao_ortho_basis(mu,i)
endif
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_mono_elec_integral, (ao_num_align, ao_num) ]
implicit none
BEGIN_DOC
! h core in orthogonal AO basis
END_DOC
double precision, allocatable :: T(:,:)
allocate(T(ao_num,ao_num))
call dgemm('T','N',ao_num,ao_num,ao_num,1.d0, &
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
ao_mono_elec_integral, size(ao_mono_elec_integral,1), &
0.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, size(ao_ortho_canonical_coef,1), &
0.d0, ao_ortho_mono_elec_integral, size(ao_ortho_mono_elec_integral,1))
deallocate(T)
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_dressing, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the AO basis set
END_DOC
call ao_ortho_cano_to_ao(ao_ortho_mono_elec_integral_dressing,size(ao_ortho_mono_elec_integral_dressing,1),&
ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral_dressing, (mo_tot_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Dressing of the core hamiltonian in the MO basis set
END_DOC
call ao_to_mo(ao_mono_elec_integral_dressing,size(ao_mono_elec_integral_dressing,1),&
mo_mono_elec_integral_dressing,size(mo_mono_elec_integral_dressing,1))
END_PROVIDER
BEGIN_PROVIDER [ integer, idx_dressing ]
implicit none
BEGIN_DOC
! Index of the MO which is being dressed
END_DOC
idx_dressing = 1
END_PROVIDER
BEGIN_PROVIDER [ double precision, cusp_corrected_mos, (ao_num_align,mo_tot_num) ]
implicit none
BEGIN_DOC
! Dressing core hamiltonian in the AO basis set
END_DOC
integer :: i,j
double precision, allocatable :: F(:,:), M(:,:)
allocate(F(mo_tot_num_align,mo_tot_num),M(ao_num,mo_tot_num))
logical :: oneshot
! oneshot = .True.
oneshot = .False.
if (oneshot) then
cusp_corrected_mos(1:ao_num,1:mo_tot_num) = mo_coef(1:ao_num,1:mo_tot_num)
slater_coef(1:nucl_num,1:mo_tot_num) = cusp_C(1:nucl_num,1:mo_tot_num)
return
else
do idx_dressing=1,mo_tot_num
if (idx_dressing>1) then
TOUCH idx_dressing
endif
do j=1,mo_tot_num
do i=1,mo_tot_num
F(i,j) = Fock_matrix_mo(i,j)
enddo
enddo
do j=1,mo_tot_num
do i=1,ao_num
M(i,j) = mo_coef(i,j)
enddo
enddo
integer :: it
do it=1,128
! print *, 'X', ao_ortho_canonical_coef(1:ao_num,1:ao_num)
! print *, 'C', mo_coef(1:ao_num,1:mo_tot_num)
! print *, 'Cp', mo_coef_in_ao_ortho_basis(1:ao_num,1:mo_tot_num)
! print *, 'cAi', cusp_C(1:nucl_num,1:mo_tot_num)
! print *, 'FmuA', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
! print *, 'Fock:', Fock_matrix_ao(1:ao_num,1:ao_num)
! print *, 'Diag Dressing:', ao_ortho_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'Dressing:', ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'Dressed Fock:', Fock_matrix_ao(1:ao_num,1:ao_num) + ao_mono_elec_integral_dressing(1:ao_num,1:ao_num)
! print *, 'AO_orthoSlaOverlap_matrix', AO_orthoSlaOverlap_matrix(1:ao_num,1:nucl_num)
! print *, 'AO_orthoSlaH_matrix', AO_orthoSlaH_matrix(1:ao_num,1:nucl_num)
! print *, 'ao_ortho_mono_elec_integral', ao_ortho_mono_elec_integral(1:ao_num,1:ao_num)
! print *, 'Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num)
do j=1,mo_tot_num
do i=1,mo_tot_num
Fock_matrix_mo(i,j) += mo_mono_elec_integral_dressing(i,j)
enddo
enddo
do i=1,mo_tot_num
Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)
enddo
! print *, 'Dressed Fock MO:', Fock_matrix_mo(1:mo_tot_num,1:mo_tot_num)
double precision :: conv
conv = 0.d0
do j=1,mo_tot_num
do i=1,mo_tot_num
if (i==j) cycle
conv = max(conv,Fock_matrix_mo(i,j))
enddo
enddo
TOUCH Fock_matrix_mo Fock_matrix_diag_mo
mo_coef(1:ao_num,1:mo_tot_num) = eigenvectors_fock_matrix_mo(1:ao_num,1:mo_tot_num)
TOUCH mo_coef
!print *, 'C', mo_coef(1:ao_num,1:mo_tot_num)
!print *, '-----'
print *, idx_dressing, it, real(mo_coef(1,idx_dressing)), real(conv)
if (conv < 1.d-5) exit
!stop
enddo
cusp_corrected_mos(1:ao_num,idx_dressing) = mo_coef(1:ao_num,idx_dressing)
slater_coef(1:nucl_num,idx_dressing) = cusp_C(1:nucl_num,idx_dressing)
enddo
idx_dressing = 1
mo_coef(1:ao_num,1:mo_tot_num) = M(1:ao_num,1:mo_tot_num)
soft_TOUCH mo_coef idx_dressing slater_coef
endif
END_PROVIDER

7
plugins/Hartree_Fock_SlaterDressed/integrals.irp.f
View File

@ -344,7 +344,7 @@ subroutine GauSlaNuclear(expGau,cGau,aGau,expSla,cSla,ZNuc,cNuc,result)
ss = k*ss
! Print result
write(*,*) ss
! write(*,*) ss
result = 0.d0
end
@ -502,8 +502,8 @@ BEGIN_PROVIDER [ double precision, MOSla$X_matrix, (mo_tot_num, nucl_num) ]
BEGIN_DOC
! <MO | Slater>
END_DOC
call dgemm('N','N',mo_tot_num,nucl_num,ao_num,1.d0, &
mo_coef_transp, size(mo_coef_transp,1), &
call dgemm('T','N',mo_tot_num,nucl_num,ao_num,1.d0, &
mo_coef, size(mo_coef,1), &
GauSla$X_matrix, size(GauSla$X_matrix,1), &
0.d0, MOSla$X_matrix, size(MOSla$X_matrix,1))
END_PROVIDER
@ -520,6 +520,7 @@ BEGIN_PROVIDER [ double precision, AO_orthoSla$X_matrix, (ao_num, nucl_num) ]
END_PROVIDER
SUBST [ X ]
Overlap ;;

5
plugins/Hartree_Fock_SlaterDressed/slater.irp.f
View File

@ -8,7 +8,10 @@ BEGIN_PROVIDER [ double precision, slater_expo, (nucl_num) ]
if (exists) then
slater_expo(1:nucl_num) = slater_expo_ezfio(1:nucl_num)
else
slater_expo(1:nucl_num) = nucl_charge(1:nucl_num)
integer :: i
do i=1,nucl_num
slater_expo(i) = nucl_charge(i)
enddo
call ezfio_set_Hartree_Fock_SlaterDressed_slater_expo_ezfio(slater_expo)
endif
END_PROVIDER

159
src/Determinants/mo_energy_expval.irp.f Normal file
View File

@ -0,0 +1,159 @@
BEGIN_PROVIDER [ double precision, mo_energy_expval, (N_states,mo_tot_num,2,2)]
use bitmasks
implicit none
BEGIN_DOC
! Third index is spin.
! Fourth index is 1:creation, 2:annihilation
END_DOC
integer :: i,j,k
integer :: ispin, istate
integer :: hp
double precision :: norm_out(N_states)
integer, parameter :: hole_particle(2) = (/ -1, 1 /)
double precision :: energies(n_states)
integer(bit_kind), allocatable :: psi_in_out(:,:,:)
double precision, allocatable :: psi_in_out_coef(:,:)
double precision :: E0(N_states), norm
double precision, parameter :: t=1.d-3
allocate (psi_in_out(N_int,2,N_det),psi_in_out_coef(N_det,N_states))
mo_energy_expval = 0.d0
psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states)
psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det)
! Truncate the wave function
do istate=1,N_states
norm = 0.d0
do k=1,N_det
if (dabs(psi_in_out_coef(k,istate)) < t) then
psi_in_out_coef(k,istate) = 0.d0
endif
norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate)
enddo
ASSERT (norm > 0.d0)
norm = 1.d0/dsqrt(norm)
psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm
call au0_h_au0(E0,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1))
enddo
do hp=1,2
do ispin=1,2
do i=1,mo_tot_num
psi_in_out_coef(1:N_det,1:N_states) = psi_coef(1:N_det,1:N_states)
psi_in_out(1:N_int,1:2,1:N_det) = psi_det(1:N_int,1:2,1:N_det)
call apply_exc_to_psi(i,hole_particle(hp),ispin, &
norm_out,psi_in_out,psi_in_out_coef, N_det,N_det,N_det,N_states)
! Truncate the wave function
do istate=1,N_states
norm = 0.d0
do k=1,N_det
if (dabs(psi_in_out_coef(k,istate)) < t) then
psi_in_out_coef(k,istate) = 0.d0
endif
norm = norm + psi_in_out_coef(k,istate)*psi_in_out_coef(k,istate)
enddo
if (norm == 0.d0) then
cycle
endif
norm = 1.d0/dsqrt(norm)
psi_in_out_coef(1:N_det,istate) = psi_in_out_coef(1:N_det,istate) * norm
enddo
call au0_h_au0(energies,psi_in_out,psi_in_out_coef,N_det,size(psi_in_out_coef,1))
mo_energy_expval(1:N_states,i,ispin,hp) = energies(1:N_states) - E0(1:N_states)
print *, i, ispin, real(energies(1)), real(E0(1))
enddo
enddo
enddo
mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1) = -mo_energy_expval(1:N_states,1:mo_tot_num,1:2,1)
END_PROVIDER
subroutine au0_h_au0(energies,psi_in,psi_in_coef,ndet,dim_psi_coef)
use bitmasks
implicit none
integer, intent(in) :: ndet,dim_psi_coef
integer(bit_kind), intent(in) :: psi_in(N_int,2,Ndet)
double precision, intent(in) :: psi_in_coef(dim_psi_coef,N_states)
double precision, intent(out) :: energies(N_states)
integer :: i,j, istate
double precision :: hij,accu
double precision, allocatable :: psi_coef_tmp(:)
energies(1:N_states) = 0.d0
do i = 1, Ndet
if(sum(dabs(psi_in_coef(i,1:N_states)))==0.d0) then
cycle
endif
call diag_H_mat_elem_au0_h_au0(psi_in(1,1,i),N_int,hij)
do istate=1,N_states
energies(istate) += psi_in_coef(i,istate) * psi_in_coef(i,istate) * hij
enddo
do j = i+1, Ndet
if(sum(dabs(psi_in_coef(j,1:N_states)))==0.d0) then
cycle
endif
call i_H_j(psi_in(1,1,i),psi_in(1,1,j),N_int,hij)
hij = hij+hij
do istate=1,N_states
energies(istate) = energies(istate) + psi_in_coef(i,istate) * psi_in_coef(j,istate) * hij
enddo
enddo
enddo
end
subroutine diag_H_mat_elem_au0_h_au0(det_in,Nint,hii)
use bitmasks
implicit none
BEGIN_DOC
! Computes <i|H|i> for any determinant i. Used for wave functions with an additional electron.
END_DOC
integer,intent(in) :: Nint
integer(bit_kind),intent(in) :: det_in(Nint,2)
double precision, intent(out) :: hii
integer :: i, j, iorb, jorb
integer :: occ(Nint*bit_kind_size,2)
integer :: elec_num_tab_local(2)
hii = 0.d0
call bitstring_to_list(det_in(1,1), occ(1,1), elec_num_tab_local(1), Nint)
call bitstring_to_list(det_in(1,2), occ(1,2), elec_num_tab_local(2), Nint)
! alpha - alpha
do i = 1, elec_num_tab_local(1)
iorb = occ(i,1)
hii += mo_mono_elec_integral(iorb,iorb)
do j = i+1, elec_num_tab_local(1)
jorb = occ(j,1)
hii += mo_bielec_integral_jj_anti(jorb,iorb)
enddo
enddo
! beta - beta
do i = 1, elec_num_tab_local(2)
iorb = occ(i,2)
hii += mo_mono_elec_integral(iorb,iorb)
do j = i+1, elec_num_tab_local(2)
jorb = occ(j,2)
hii += mo_bielec_integral_jj_anti(jorb,iorb)
enddo
enddo
! alpha - beta
do i = 1, elec_num_tab_local(2)
iorb = occ(i,2)
do j = 1, elec_num_tab_local(1)
jorb = occ(j,1)
hii += mo_bielec_integral_jj(jorb,iorb)
enddo
enddo
end

2
src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f
View File

@ -56,7 +56,6 @@ subroutine push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value,
stop 'error'
endif
! Activate is zmq_socket_push is a REQ
IRP_IF ZMQ_PUSH
IRP_ELSE
integer :: idummy
@ -189,7 +188,6 @@ subroutine ao_bielec_integrals_in_map_collector
rc = f77_zmq_recv( zmq_socket_pull, task_id, 4, 0)
! Activate if zmq_socket_pull is a REP
IRP_IF ZMQ_PUSH
IRP_ELSE
rc = f77_zmq_send( zmq_socket_pull, 0, 4, 0)

17
src/MO_Basis/ao_ortho_canonical.irp.f
View File

@ -76,19 +76,20 @@ BEGIN_PROVIDER [ double precision, ao_cart_to_sphe_inv, (ao_cart_to_sphe_num,ao_
! AO_cart_to_sphe_coef^(-1)
END_DOC
call get_pseudo_inverse(ao_cart_to_sphe_coef,ao_num,ao_cart_to_sphe_num, &
ao_cart_to_sphe_inv, size(ao_cart_to_sphe_coef,1))
call get_pseudo_inverse(ao_cart_to_sphe_coef,size(ao_cart_to_sphe_coef,1),&
ao_num,ao_cart_to_sphe_num, &
ao_cart_to_sphe_inv, size(ao_cart_to_sphe_inv,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num,ao_num)]
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef_inv, (ao_num_align,ao_num)]
implicit none
BEGIN_DOC
! ao_ortho_canonical_coef^(-1)
END_DOC
call get_pseudo_inverse(ao_ortho_canonical_coef,ao_num,ao_num, &
ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef,1))
call get_inverse(ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),&
ao_num, ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1))
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
@ -100,11 +101,15 @@ END_PROVIDER
! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
END_DOC
integer :: i
ao_ortho_canonical_coef(:,:) = 0.d0
ao_ortho_canonical_coef = 0.d0
do i=1,ao_num
ao_ortho_canonical_coef(i,i) = 1.d0
enddo
call ortho_lowdin(ao_overlap,size(ao_overlap,1),ao_num,ao_ortho_canonical_coef,size(ao_ortho_canonical_coef,1),ao_num)
ao_ortho_canonical_num=ao_num
return
if (ao_cartesian) then
ao_ortho_canonical_num = ao_num

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

@ -72,8 +72,10 @@ BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num_align, mo_
implicit none
BEGIN_DOC
! MO coefficients in orthogonalized AO basis
!
! C^(-1).C_mo
END_DOC
call dgemm('T','N',ao_num,mo_tot_num,ao_num,1.d0, &
call dgemm('N','N',ao_num,mo_tot_num,ao_num,1.d0, &
ao_ortho_canonical_coef_inv, size(ao_ortho_canonical_coef_inv,1),&
mo_coef, size(mo_coef,1), 0.d0, &
mo_coef_in_ao_ortho_basis, size(mo_coef_in_ao_ortho_basis,1))
@ -155,6 +157,8 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
implicit none
BEGIN_DOC
! Transform A from the AO basis to the MO basis
!
! C.A_ao.Ct
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_ao(LDA_ao)
@ -181,6 +185,8 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the AO basis
!
! (S.C).A_mo.(S.C)t
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo)
@ -273,6 +279,8 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
implicit none
BEGIN_DOC
! Transform A from the AO basis to the orthogonal AO basis
!
! C^(-1).A_ao.Ct^(-1)
END_DOC
integer, intent(in) :: LDA_ao,LDA
double precision, intent(in) :: A_ao(LDA_ao,*)
@ -296,35 +304,3 @@ subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA)
deallocate(T)
end
subroutine mo_to_ao_ortho_cano(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the AO orthogonal basis to the AO basis
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
double precision, allocatable :: T(:,:), SC(:,:)
allocate ( SC(ao_num_align,mo_tot_num) )
allocate ( T(mo_tot_num_align,ao_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
call dgemm('N','N', ao_num, mo_tot_num, ao_num, &
1.d0, ao_overlap,size(ao_overlap,1), &
ao_ortho_canonical_coef, size(ao_ortho_canonical_coef,1), &
0.d0, SC, ao_num_align)
call dgemm('N','T', mo_tot_num, ao_num, mo_tot_num, &
1.d0, A_mo,LDA_mo, &
SC, size(SC,1), &
0.d0, T, mo_tot_num_align)
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, SC,size(SC,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
deallocate(T,SC)
end

40
src/Utils/LinearAlgebra.irp.f
View File

@ -72,8 +72,6 @@ subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
double precision, allocatable :: S_half(:,:)
!DEC$ ATTRIBUTES ALIGN : 64 :: U, Vt, D
integer :: info, i, j
!call ortho_lowdin(overlap,LDA,N,C,LDC,m)
!return
if (n < 2) then
return
@ -270,14 +268,42 @@ end
subroutine get_pseudo_inverse(A,m,n,C,LDA)
subroutine get_inverse(A,LDA,m,C,LDC)
implicit none
BEGIN_DOC
! Returns the inverse of the square matrix A
END_DOC
integer, intent(in) :: m, LDA, LDC
double precision, intent(in) :: A(LDA,m)
double precision, intent(out) :: C(LDC,m)
integer :: info,lwork
integer, allocatable :: ipiv(:)
double precision,allocatable :: work(:)
allocate (ipiv(ao_num), work(ao_num*ao_num))
lwork = size(work)
C(1:m,1:m) = A(1:m,1:m)
call dgetrf(m,m,C,size(C,1),ipiv,info)
if (info /= 0) then
print *, info
stop 'error in inverse (dgetrf)'
endif
call dgetri(m,C,size(C,1),ipiv,work,lwork,info)
if (info /= 0) then
print *, info
stop 'error in inverse (dgetri)'
endif
deallocate(ipiv,work)
end
subroutine get_pseudo_inverse(A,LDA,m,n,C,LDC)
implicit none
BEGIN_DOC
! Find C = A^-1
END_DOC
integer, intent(in) :: m,n, LDA
integer, intent(in) :: m,n, LDA, LDC
double precision, intent(in) :: A(LDA,n)
double precision, intent(out) :: C(n,m)
double precision, intent(out) :: C(LDC,m)
double precision, allocatable :: U(:,:), D(:), Vt(:,:), work(:), A_tmp(:,:)
integer :: info, lwork
@ -304,7 +330,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
endif
do i=1,n
if (dabs(D(i)) > 1.d-6) then
if (D(i)/D(1) > 1.d-10) then
D(i) = 1.d0/D(i)
else
D(i) = 0.d0
@ -315,7 +341,7 @@ subroutine get_pseudo_inverse(A,m,n,C,LDA)
do i=1,m
do j=1,n
do k=1,n
C(j,i) += U(i,k) * D(k) * Vt(k,j)
C(j,i) = C(j,i) + U(i,k) * D(k) * Vt(k,j)
enddo
enddo
enddo