LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-01 19:05:25 +02:00
Code Issues Releases Wiki Activity

ao_to_mo is broken

Browse Source
This commit is contained in:
Anthony Scemama 2017-09-12 01:57:45 +02:00
parent fa58b656f8
commit 485ffb4bef
15 changed files with 62 additions and 84 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
plugins/Hartree_Fock/EZFIO.cfg
View File

@ -32,7 +32,7 @@ default: 500
type: Positive_float
doc: Energy shift on the virtual MOs to improve SCF convergence
interface: ezfio,provider,ocaml
default: 0.0
default: 0.1
[scf_algorithm]
type: character*(32)

6
plugins/Hartree_Fock/HF_density_matrix_ao.irp.f
View File

@ -1,7 +1,7 @@
BEGIN_PROVIDER [double precision, HF_density_matrix_ao_alpha, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! S^-1 x Alpha density matrix in the AO basis x S^-1
! Alpha density matrix in the AO basis
END_DOC
call dgemm('N','T',ao_num,ao_num,elec_alpha_num,1.d0, &
@ -14,7 +14,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao_beta, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! S^-1 Beta density matrix in the AO basis x S^-1
! Beta density matrix in the AO basis
END_DOC
call dgemm('N','T',ao_num,ao_num,elec_beta_num,1.d0, &
@ -27,7 +27,7 @@ END_PROVIDER
BEGIN_PROVIDER [ double precision, HF_density_matrix_ao, (ao_num,ao_num) ]
implicit none
BEGIN_DOC
! S^-1 Density matrix in the AO basis S^-1
! Density matrix in the AO basis
END_DOC
ASSERT (size(HF_density_matrix_ao,1) == size(HF_density_matrix_ao_alpha,1))
if (elec_alpha_num== elec_beta_num) then

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

@ -83,6 +83,7 @@ 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

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

@ -1,7 +1,7 @@
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) ]
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) ]
implicit none
BEGIN_DOC
! Overlap between atomic basis functions:
@ -34,8 +34,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_overlap(i,j)= 0.d0
ao_overlap_x(i,j)= 0.d0
@ -49,7 +47,6 @@
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)
@ -72,7 +69,7 @@
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between absolute value of atomic basis functions:
@ -103,8 +100,6 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_overlap_abs(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
@ -115,7 +110,6 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1)
@ -129,3 +123,11 @@ BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_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))
END_PROVIDER

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

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_mono_elec_integral_diag,(ao_num)]
implicit none
integer :: i,j,n,l
@ -7,7 +7,6 @@
! : sum of the kinetic and nuclear electronic potential
END_DOC
do j = 1, ao_num
!DIR$ VECTOR ALIGNED
do i = 1, ao_num
ao_mono_elec_integral(i,j) = ao_nucl_elec_integral(i,j) + ao_kinetic_integral(i,j) + ao_pseudo_integral(i,j)
enddo

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

@ -1,6 +1,6 @@
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) ]
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) ]
implicit none
integer :: i,j,n,l
double precision :: f
@ -45,8 +45,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
ao_deriv2_x(i,j)= 0.d0
ao_deriv2_y(i,j)= 0.d0
@ -59,7 +57,6 @@
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1)
@ -122,7 +119,7 @@
END_PROVIDER
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! array of the priminitve basis kinetic integrals
@ -137,16 +134,11 @@ BEGIN_PROVIDER [double precision, ao_kinetic_integral, (ao_num_align,ao_num)]
else
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i,j) &
!$OMP SHARED(ao_num, ao_num_align, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
!$OMP SHARED(ao_num, ao_kinetic_integral,ao_deriv2_x,ao_deriv2_y,ao_deriv2_z)
do j = 1, ao_num
!DEC$ VECTOR ALWAYS
!DEC$ VECTOR ALIGNED
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_align,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_kinetic_integral, (mo_tot_num,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_align,mo_tot_num)]
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num,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_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num,ao_num)]
BEGIN_DOC
! interaction nuclear electron
END_DOC
@ -80,7 +80,7 @@ BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral, (ao_num_align,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num_align,ao_num,nucl_num)]
BEGIN_PROVIDER [ double precision, ao_nucl_elec_integral_per_atom, (ao_num,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_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! Pseudo-potential integrals
@ -29,7 +29,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral, (ao_num_align,ao_num)]
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num,ao_num)]
implicit none
BEGIN_DOC
! Local pseudo-potential
@ -128,7 +128,7 @@ BEGIN_PROVIDER [ double precision, ao_pseudo_integral_local, (ao_num_align,ao_nu
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num_align,ao_num)]
BEGIN_PROVIDER [ double precision, ao_pseudo_integral_non_local, (ao_num,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_align,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral, (mo_tot_num,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_align,mo_to
END_PROVIDER
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num_align,mo_tot_num,nucl_num)]
BEGIN_PROVIDER [double precision, mo_nucl_elec_integral_per_atom, (mo_tot_num,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_align,mo_tot_num)]
BEGIN_PROVIDER [double precision, mo_pseudo_integral, (mo_tot_num,mo_tot_num)]
implicit none
BEGIN_DOC
! interaction nuclear electron on the MO basis

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

@ -1,6 +1,6 @@
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_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_DOC
! array of the integrals of AO_i * x^2 AO_j
! array of the integrals of AO_i * y^2 AO_j
@ -35,8 +35,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -49,7 +47,6 @@
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
c = ao_coef_normalized_ordered_transp(n,j)*ao_coef_normalized_ordered_transp(l,i)
beta = ao_expo_ordered_transp(l,i)
@ -72,9 +69,9 @@
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_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_DOC
! array of the integrals of AO_i * x AO_j
! array of the integrals of AO_i * y AO_j
@ -109,8 +106,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -123,7 +118,6 @@
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
c = ao_coef_normalized_ordered_transp(l,i)*ao_coef_normalized_ordered_transp(n,j)
@ -145,9 +139,9 @@
!$OMP END PARALLEL DO
END_PROVIDER
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_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_DOC
! array of the integrals of AO_i * d/dx AO_j
! array of the integrals of AO_i * d/dy AO_j
@ -183,8 +177,6 @@
power_A(1) = ao_power( j, 1 )
power_A(2) = ao_power( j, 2 )
power_A(3) = ao_power( j, 3 )
!DEC$ VECTOR ALIGNED
!DEC$ VECTOR ALWAYS
do i= 1,ao_num
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
@ -197,7 +189,6 @@
accu_z = 0.d0
do n = 1,ao_prim_num(j)
alpha = ao_expo_ordered_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(l,i)
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1)

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_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_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_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_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_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_DOC
! array of the integrals of MO_i * x^2 MO_j
! array of the integrals of MO_i * y^2 MO_j

33
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
@ -125,7 +125,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))
@ -162,30 +162,23 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo)
! (S.C)t.A_ao.S.C
END_DOC
integer, intent(in) :: LDA_ao,LDA_mo
double precision, intent(in) :: A_ao(LDA_ao)
double precision, intent(out) :: A_mo(LDA_mo)
double precision, allocatable :: T(:,:), SC(:,:)
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 ( SC(ao_num_align,mo_tot_num) )
allocate ( T(ao_num_align,mo_tot_num) )
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T
allocate ( T(mo_tot_num,ao_num) )
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, ao_overlap,size(ao_overlap,1), &
mo_coef, size(mo_coef,1), &
0.d0, SC, ao_num_align)
call dgemm('T','N', ao_num, mo_tot_num, ao_num, &
1.d0, SC, size(SC,1), &
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, &
call dgemm('N','N', mo_tot_num, mo_tot_num, ao_num, &
1.d0, T,size(T,1), &
SC, size(SC,1), &
S_mo_coef, size(S_mo_coef,1), &
0.d0, A_mo, size(A_mo,1))
deallocate(T,SC)
deallocate(T)
end
subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
@ -196,8 +189,8 @@ subroutine mo_to_ao(A_mo,LDA_mo,A_ao,LDA_ao)
! C.A_mo.Ct
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, intent(in) :: A_mo(LDA_mo,mo_tot_num)
double precision, intent(out) :: A_ao(LDA_ao,ao_num)
double precision, allocatable :: T(:,:)
allocate ( T(mo_tot_num_align,ao_num) )