LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-03 18:16:12 +01:00
Code Issues Releases Wiki Activity

Molden interface

Browse Source
This commit is contained in:
Manu 2015-01-07 16:48:03 +01:00
parent 1deea047f9
commit 0242bf9303
12 changed files with 640 additions and 179 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

250
src/AOs/aos.irp.f
View File

@ -19,13 +19,10 @@ END_PROVIDER
BEGIN_PROVIDER [ integer, ao_power, (ao_num_align,3) ]
&BEGIN_PROVIDER [ double precision, ao_expo, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_coef, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
implicit none
BEGIN_DOC
! Coefficients, exponents and powers of x,y and z
! ao_coef(i,j) = coefficient of the jth primitive on the ith ao
! ao_l = l value of the AO: a+b+c in x^a y^b z^c
END_DOC
PROVIDE ezfio_filename
@ -33,6 +30,7 @@ END_PROVIDER
allocate ( buffer(ao_num,ao_prim_num_max) )
integer :: ibuffer(ao_num,3)
integer :: i,j,k
character*(128) :: give_ao_character_space
ibuffer = 0
call ezfio_get_ao_basis_ao_power(ibuffer)
ao_power = 0
@ -49,6 +47,7 @@ END_PROVIDER
ao_expo(i,j) = buffer(i,j)
enddo
enddo
ao_coef = 0.d0
buffer = 0.d0
call ezfio_get_ao_basis_ao_coef(buffer)
@ -57,6 +56,7 @@ END_PROVIDER
ao_coef(i,j) = buffer(i,j)
enddo
enddo
deallocate(buffer)
! Normalization of the AO coefficients
@ -93,12 +93,66 @@ END_PROVIDER
ao_coef(i,j) = d(j,2)
enddo
enddo
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
implicit none
BEGIN_DOC
! Transposed ao_coef and ao_expo
END_DOC
integer :: i,j
do j=1, ao_num
do i=1, ao_prim_num_max
ao_coef_transp(i,j) = ao_coef(j,i)
ao_expo_transp(i,j) = ao_expo(j,i)
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_unnormalized, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ double precision, ao_expo_unsorted, (ao_num_align,ao_prim_num_max) ]
&BEGIN_PROVIDER [ integer, ao_l, (ao_num) ]
&BEGIN_PROVIDER [ character*(128), ao_l_char, (ao_num) ]
implicit none
BEGIN_DOC
! Coefficients, exponents and powers of x,y and z as in the EZFIO file
! ao_coef(i,j) = coefficient of the jth primitive on the ith ao
! ao_l = l value of the AO: a+b+c in x^a y^b z^c
END_DOC
PROVIDE ezfio_filename
double precision, allocatable :: buffer(:,:)
allocate ( buffer(ao_num,ao_prim_num_max) )
integer :: i,j,k
character*(128) :: give_ao_character_space
buffer = 0.d0
call ezfio_get_ao_basis_ao_expo(buffer)
do j = 1, ao_prim_num_max
do i = 1, ao_num
ao_expo_unsorted(i,j) = buffer(i,j)
enddo
enddo
buffer = 0.d0
call ezfio_get_ao_basis_ao_coef(buffer)
do j = 1, ao_prim_num_max
do i = 1, ao_num
ao_coef_unnormalized(i,j) = buffer(i,j)
enddo
enddo
deallocate(buffer)
do i=1,ao_num
ao_l(i) = ao_power(i,1) + ao_power(i,2) + ao_power(i,3)
ao_l_char(i) = l_to_charater(ao_l(i))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_coef_transp, (ao_prim_num_max_align,ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_expo_transp, (ao_prim_num_max_align,ao_num) ]
@ -159,3 +213,187 @@ BEGIN_PROVIDER [ integer, ao_nucl, (ao_num)]
call ezfio_get_ao_basis_ao_nucl(ao_nucl)
END_PROVIDER
BEGIN_PROVIDER [ character*(128), l_to_charater, (0:4)]
BEGIN_DOC
! character corresponding to the "L" value of an AO orbital
END_DOC
implicit none
l_to_charater(0)='S'
l_to_charater(1)='P'
l_to_charater(2)='D'
l_to_charater(3)='F'
l_to_charater(4)='G'
END_PROVIDER
BEGIN_PROVIDER [ integer, Nucl_N_Aos, (nucl_num)]
&BEGIN_PROVIDER [ integer, N_AOs_max ]
implicit none
integer :: i
BEGIN_DOC
! Number of AOs per atom
END_DOC
Nucl_N_Aos = 0
do i = 1, ao_num
Nucl_N_Aos(ao_nucl(i)) +=1
enddo
N_AOs_max = maxval(Nucl_N_Aos)
END_PROVIDER
BEGIN_PROVIDER [ integer, Nucl_Aos, (nucl_num,N_AOs_max)]
implicit none
BEGIN_DOC
! List of AOs attached on each atom
END_DOC
integer :: i
integer, allocatable :: nucl_tmp(:)
allocate(nucl_tmp(nucl_num))
nucl_tmp = 0
Nucl_Aos = 0
do i = 1, ao_num
nucl_tmp(ao_nucl(i))+=1
Nucl_Aos(ao_nucl(i),nucl_tmp(ao_nucl(i))) = i
enddo
deallocate(nucl_tmp)
END_PROVIDER
BEGIN_PROVIDER [ integer, Nucl_list_shell_Aos, (nucl_num,N_AOs_max)]
&BEGIN_PROVIDER [ integer, Nucl_num_shell_Aos, (nucl_num)]
implicit none
integer :: i,j,k
BEGIN_DOC
! Index of the shell type Aos and of the corresponding Aos
! Per convention, for P,D,F and G AOs, we take the index
! of the AO with the the corresponding power in the "X" axis
END_DOC
do i = 1, nucl_num
Nucl_num_shell_Aos(i) = 0
do j = 1, Nucl_N_Aos(i)
if(ao_l(Nucl_Aos(i,j))==0)then
! S type function
Nucl_num_shell_Aos(i)+=1
Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
elseif(ao_l(Nucl_Aos(i,j))==1)then
! P type function
if(ao_power(Nucl_Aos(i,j),1)==1)then
Nucl_num_shell_Aos(i)+=1
Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
endif
elseif(ao_l(Nucl_Aos(i,j))==2)then
! D type function
if(ao_power(Nucl_Aos(i,j),1)==2)then
Nucl_num_shell_Aos(i)+=1
Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
endif
elseif(ao_l(Nucl_Aos(i,j))==3)then
! F type function
if(ao_power(Nucl_Aos(i,j),1)==3)then
Nucl_num_shell_Aos(i)+=1
Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
endif
elseif(ao_l(Nucl_Aos(i,j))==4)then
! G type function
if(ao_power(Nucl_Aos(i,j),1)==4)then
Nucl_num_shell_Aos(i)+=1
Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
endif
endif
enddo
enddo
END_PROVIDER
BEGIN_PROVIDER [ character*(4), ao_l_char_space, (ao_num) ]
implicit none
integer :: i
character*(4) :: give_ao_character_space
do i=1,ao_num
if(ao_l(i)==0)then
! S type AO
give_ao_character_space = 'S '
elseif(ao_l(i) == 1)then
! P type AO
if(ao_power(i,1)==1)then
give_ao_character_space = 'X '
elseif(ao_power(i,2) == 1)then
give_ao_character_space = 'Y '
else
give_ao_character_space = 'Z '
endif
elseif(ao_l(i) == 2)then
! D type AO
if(ao_power(i,1)==2)then
give_ao_character_space = 'XX '
elseif(ao_power(i,2) == 2)then
give_ao_character_space = 'YY '
elseif(ao_power(i,3) == 2)then
give_ao_character_space = 'ZZ '
elseif(ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'XY '
elseif(ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'XZ '
else
give_ao_character_space = 'YZ '
endif
elseif(ao_l(i) == 3)then
! F type AO
if(ao_power(i,1)==3)then
give_ao_character_space = 'XXX '
elseif(ao_power(i,2) == 3)then
give_ao_character_space = 'YYY '
elseif(ao_power(i,3) == 3)then
give_ao_character_space = 'ZZZ '
elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'XXY '
elseif(ao_power(i,1) == 2 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'XXZ '
elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1)then
give_ao_character_space = 'YYX '
elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'YYZ '
elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1)then
give_ao_character_space = 'ZZX '
elseif(ao_power(i,3) == 2 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'ZZY '
elseif(ao_power(i,3) == 1 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'XYZ '
endif
elseif(ao_l(i) == 4)then
! G type AO
if(ao_power(i,1)==4)then
give_ao_character_space = 'XXXX'
elseif(ao_power(i,2) == 4)then
give_ao_character_space = 'YYYY'
elseif(ao_power(i,3) == 4)then
give_ao_character_space = 'ZZZZ'
elseif(ao_power(i,1) == 3 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'XXXY'
elseif(ao_power(i,1) == 3 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'XXXZ'
elseif(ao_power(i,2) == 3 .and. ao_power(i,1) == 1)then
give_ao_character_space = 'YYYX'
elseif(ao_power(i,2) == 3 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'YYYZ'
elseif(ao_power(i,3) == 3 .and. ao_power(i,1) == 1)then
give_ao_character_space = 'ZZZX'
elseif(ao_power(i,3) == 3 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'ZZZY'
elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 2)then
give_ao_character_space = 'XXYY'
elseif(ao_power(i,2) == 2 .and. ao_power(i,3) == 2)then
give_ao_character_space = 'YYZZ'
elseif(ao_power(i,1) == 2 .and. ao_power(i,2) == 1 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'XXYZ'
elseif(ao_power(i,2) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,3) == 1)then
give_ao_character_space = 'YYXZ'
elseif(ao_power(i,3) == 2 .and. ao_power(i,1) == 1 .and. ao_power(i,2) == 1)then
give_ao_character_space = 'ZZXY'
endif
endif
ao_l_char_space(i) = give_ao_character_space
enddo
END_PROVIDER

4
src/CIS_dressed/CIS_full.irp.f
View File

@ -6,7 +6,9 @@ program CIS_full
end
subroutine save_cis
do i = 1, n_state_cis
print*,'eigenvalues_CIS(i) = ',eigenvalues_CIS(i)
print*,''
print*,'eigenvalues_CIS(i) = ',eigenvalues_CIS(i)+nuclear_repulsion
print*,'s2 = ',s_2_CIS(i)
enddo
call save_wavefunction_general(size_psi_CIS,n_state_cis,psi_CIS,coefs_CIS)

4
src/MOs/NEEDED_MODULES
View File

@ -1 +1,3 @@
AOs Ezfio_files Nuclei Output Utils Electrons
AOs Electrons Ezfio_files Nuclei Output Utils

35
src/MOs/mo_overlap.irp.f Normal file
View File

@ -0,0 +1,35 @@
BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l
double precision :: f
integer :: lmax
lmax = iand(ao_num,-4)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
!$OMP PRIVATE(i,j,n,l) &
!$OMP SHARED(mo_overlap,mo_coef,ao_overlap, &
!$OMP mo_tot_num,ao_num,lmax)
do j=1,mo_tot_num
do i= 1,mo_tot_num
mo_overlap(i,j) = 0.d0
do n = 1, lmax,4
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(l,i) * &
( mo_coef(n ,j) * ao_overlap(l,n ) &
+ mo_coef(n+1,j) * ao_overlap(l,n+1) &
+ mo_coef(n+2,j) * ao_overlap(l,n+2) &
+ mo_coef(n+3,j) * ao_overlap(l,n+3) )
enddo
enddo
do n = lmax+1, ao_num
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(n,j) * mo_coef(l,i) * ao_overlap(l,n)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER

119
src/MOs/mos.irp.f
View File

@ -121,3 +121,122 @@ 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
! Transform A from the AO basis to the MO basis
END_DOC
double precision, intent(in) :: A_ao(LDA_ao)
double precision, intent(out) :: A_mo(LDA_mo)
integer, intent(in) :: LDA_ao,LDA_mo
double precision, allocatable :: T(:,:)
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, A_ao,LDA_ao, &
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, ao_num_align, &
0.d0, A_mo, LDA_mo)
deallocate(T)
end
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
END_DOC
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
integer, intent(in) :: LDA_ao,LDA_mo
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), &
mo_coef, size(mo_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
subroutine mo_to_ao_no_overlap(A_mo,LDA_mo,A_ao,LDA_ao)
implicit none
BEGIN_DOC
! Transform A from the MO basis to the S^-1 AO basis
END_DOC
double precision, intent(in) :: A_mo(LDA_mo)
double precision, intent(out) :: A_ao(LDA_ao)
integer, intent(in) :: LDA_ao,LDA_mo
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), &
0.d0, T, mo_tot_num_align)
call dgemm('N','N', ao_num, ao_num, mo_tot_num, &
1.d0, mo_coef,size(mo_coef,1), &
T, mo_tot_num_align, &
0.d0, A_ao, LDA_ao)
deallocate(T)
end
subroutine mix_mo_jk(j,k)
implicit none
integer, intent(in) :: j,k
integer :: i,i_plus,i_minus
BEGIN_DOC
! subroutine that rotates the jth MO with the kth MO
! to give two new MO's that are
! '+' = 1/sqrt(2) (|j> + |k>)
! '-' = 1/sqrt(2) (|j> - |k>)
! by convention, the '+' MO is in the lower index (min(j,k))
! by convention, the '-' MO is in the greater index (max(j,k))
END_DOC
double precision :: array_tmp(ao_num,2),dsqrt_2
if(j==k)then
print*,'Youwant to mix two orbitals that are the same !'
print*,'It does not make sense ... '
print*,'Stopping ...'
stop
endif
array_tmp = 0.d0
dsqrt_2 = 1.d0/dsqrt(2.d0)
do i = 1, ao_num
array_tmp(i,1) = dsqrt_2 * (mo_coef(i,j) + mo_coef(i,k))
array_tmp(i,2) = dsqrt_2 * (mo_coef(i,j) - mo_coef(i,k))
enddo
i_plus = min(j,k)
i_minus = max(j,k)
do i = 1, ao_num
mo_coef(i,i_plus) = array_tmp(i,1)
mo_coef(i,i_minus) = array_tmp(i,2)
enddo
end

131
src/Molden/print_mo.irp.f Normal file
View File

@ -0,0 +1,131 @@
program print_mos
implicit none
character*(128) :: output
integer :: i_unit_output,getUnitAndOpen
provide ezfio_filename
integer :: i
print*,trim(ezfio_filename)
output=trim(ezfio_filename)//'.mol'
print*,'output = ',trim(output)
i_unit_output = getUnitAndOpen(output,'w')
print*,'i_unit_output = ',i_unit_output
call write_intro_gamess(i_unit_output)
call write_geometry(i_unit_output)
call write_Ao_basis(i_unit_output)
call write_Mo_basis(i_unit_output)
write(i_unit_output,*),''
write(i_unit_output,*),''
write(i_unit_output,*)' ------------------------'
close(i_unit_output)
end
subroutine write_intro_gamess(i_unit_output)
implicit none
integer, intent(in) :: i_unit_output
integer :: i,j,k,l
write(i_unit_output,*)' * GAMESS VERSION = 22 FEB 2006 (R5) *'
write(i_unit_output,*)' * FROM IOWA STATE UNIVERSITY *'
write(i_unit_output,*)' * M.W.SCHMIDT, K.K.BALDRIDGE, J.A.BOATZ, S.T.ELBERT, *'
write(i_unit_output,*)' * M.S.GORDON, J.H.JENSEN, S.KOSEKI, N.MATSUNAGA, *'
write(i_unit_output,*)' * K.A.NGUYEN, S.J.SU, T.L.WINDUS, *'
write(i_unit_output,*)' * TOGETHER WITH M.DUPUIS, J.A.MONTGOMERY *'
write(i_unit_output,*)' * J.COMPUT.CHEM. 14, 1347-1363(1993) *'
write(i_unit_output,*)''
end
subroutine write_geometry(i_unit_output)
implicit none
integer, intent(in) :: i_unit_output
integer :: i,j,k,l, getUnitAndOpen
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(i_unit_output,*)'ATOM ATOMIC COORDINATES (BOHR) '
write(i_unit_output,*)' CHARGE X Y Z'
do i = 1, nucl_num
! write(i_unit_output,'(A2 I3 X F3.1 X 3(F16.10))') trim(element_name(int(nucl_charge(i)))),i,(nucl_charge(i)), nucl_coord(i,1), nucl_coord(i,2), nucl_coord(i,3)
write(i_unit_output,'(A2,I1, 9X F5.1 X 3(F16.10 ,4X))') trim(element_name(int(nucl_charge(i)))),i,(nucl_charge(i)), nucl_coord(i,1), nucl_coord(i,2), nucl_coord(i,3)
enddo
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
end
subroutine write_Ao_basis(i_unit_output)
implicit none
integer, intent(in) :: i_unit_output
integer :: i,j,k,l, getUnitAndOpen
character*(128) :: character_shell
integer :: i_shell,i_prim,i_ao
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
write(i_unit_output,*)''
write(i_unit_output,*)''
write(i_unit_output,*)' ATOMIC BASIS SET'
write(i_unit_output,*)' ----------------'
write(i_unit_output,*)'THE CONTRACTED PRIMITIVE FUNCTIONS HAVE BEEN UNNORMALIZED'
write(i_unit_output,*)'THE CONTRACTED BASIS FUNCTIONS ARE NOW NORMALIZED TO UNITY'
write(i_unit_output,*)''
write(i_unit_output,*)'SHELL TYPE PRIMITIVE EXPONENT CONTRACTION COEFFICIENT(S)'
write(i_unit_output,*)''
write(i_unit_output,*)''
i_shell = 0
i_prim = 0
do i = 1, Nucl_num
write(i_unit_output,'(A2,I1)') trim(element_name(int(nucl_charge(i)))),i
write(i_unit_output,*)' '
! write(i_unit_output,*)'Nucl_num_shell_Aos(i) = ',Nucl_num_shell_Aos(i)
do j = 1, Nucl_num_shell_Aos(i)
i_shell +=1
i_ao = Nucl_list_shell_Aos(i,j)
character_shell = trim(ao_l_char(i_ao))
! write(i_unit_output,*),j,i_shell,i_ao!trim(character_shell)
do k = 1, ao_prim_num(i_ao)
i_prim +=1
write(i_unit_output,'(4X,I3,3X,A1,6X,I2,6X,F16.7,2X,F16.12)')i_shell,character_shell,i_prim,ao_expo_unsorted(i_ao,k),ao_coef_unnormalized(i_ao,k)
enddo
write(i_unit_output,*)''
enddo
enddo
write(i_unit_output,*)''
write(i_unit_output,'(A47,2X,I3)')'TOTAL NUMBER OF BASIS SET SHELLS =', i_shell
write(i_unit_output,'(A47,2X,I3)')'NUMBER OF CARTESIAN GAUSSIAN BASIS FUNCTIONS =', ao_num
write(i_unit_output,*)''
end
subroutine write_Mo_basis(i_unit_output)
implicit none
integer, intent(in) :: i_unit_output
integer :: i,j,k,l, getUnitAndOpen
integer :: i_5,i_mod
write(i_unit_output,*),' ----------------------'
write(i_unit_output,*),' MCSCF NATURAL ORBITALS'
write(i_unit_output,*),' ----------------------'
write(i_unit_output,*),' '
do j = 1, mo_tot_num
write(i_unit_output,'(18X,I3)')j
write(i_unit_output,*)''
write(i_unit_output,'(18X,F8.5)')-1.d0
write(i_unit_output,*)''
do i = 1, ao_num
write(i_unit_output,'(2X,I3, 2X A1, I3, 2X A4 , F9.6)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i)),mo_coef(i,j)
! write(i_unit_output,'(I3, X A1, X I3, X A4 X F16.8)')i,trim(element_name(int(nucl_charge(ao_nucl(i))))),ao_nucl(i),(ao_l_char_space(i))
enddo
write(i_unit_output,*)''
enddo
end

121
src/MonoInts/ao_mono_ints.irp.f
View File

@ -1,124 +1,3 @@
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:
! :math:`\int \chi_i(r) \chi_j(r) dr)`
END_DOC
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta, c
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
dim1=100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_transp,ao_nucl, &
!$OMP ao_expo_transp,dim1)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
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
ao_overlap_y(i,j)= 0.d0
ao_overlap_z(i,j)= 0.d0
B_center(1) = nucl_coord( ao_nucl(i), 1 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_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)
c = ao_coef_transp(n,j) * ao_coef_transp(l,i)
ao_overlap(i,j) += c * overlap
ao_overlap_x(i,j) += c * overlap_x
ao_overlap_y(i,j) += c * overlap_y
ao_overlap_z(i,j) += c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_overlap_abs,(ao_num_align,ao_num) ]
implicit none
BEGIN_DOC
! Overlap between absolute value of atomic basis functions:
! :math:`\int |\chi_i(r)| |\chi_j(r)| dr)`
END_DOC
integer :: i,j,n,l
double precision :: f
integer :: dim1
double precision :: overlap, overlap_x, overlap_y, overlap_z
double precision :: alpha, beta
double precision :: A_center(3), B_center(3)
integer :: power_A(3), power_B(3)
double precision :: lower_exp_val, dx
dim1=100
lower_exp_val = 40.d0
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,dx) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_abs,ao_num,ao_coef_transp,ao_nucl, &
!$OMP ao_expo_transp,dim1,lower_exp_val)
do j=1,ao_num
A_center(1) = nucl_coord( ao_nucl(j), 1 )
A_center(2) = nucl_coord( ao_nucl(j), 2 )
A_center(3) = nucl_coord( ao_nucl(j), 3 )
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 )
B_center(2) = nucl_coord( ao_nucl(i), 2 )
B_center(3) = nucl_coord( ao_nucl(i), 3 )
power_B(1) = ao_power( i, 1 )
power_B(2) = ao_power( i, 2 )
power_B(3) = ao_power( i, 3 )
do n = 1,ao_prim_num(j)
alpha = ao_expo_transp(n,j)
!DEC$ VECTOR ALIGNED
do l = 1, ao_prim_num(i)
beta = ao_expo_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)
call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1)
call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1)
ao_overlap_abs(i,j) += abs(ao_coef_transp(n,j) * ao_coef_transp(l,i)) * overlap_x * overlap_y * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, ao_mono_elec_integral,(ao_num_align,ao_num)]
implicit none
integer :: i,j,n,l

34
src/MonoInts/mo_mono_ints.irp.f
View File

@ -1,37 +1,3 @@
BEGIN_PROVIDER [ double precision, mo_overlap,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l
double precision :: f
integer :: lmax
lmax = iand(ao_num,-4)
!$OMP PARALLEL DO SCHEDULE(STATIC) DEFAULT(NONE) &
!$OMP PRIVATE(i,j,n,l) &
!$OMP SHARED(mo_overlap,mo_coef,ao_overlap, &
!$OMP mo_tot_num,ao_num,lmax)
do j=1,mo_tot_num
do i= 1,mo_tot_num
mo_overlap(i,j) = 0.d0
do n = 1, lmax,4
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(l,i) * &
( mo_coef(n ,j) * ao_overlap(l,n ) &
+ mo_coef(n+1,j) * ao_overlap(l,n+1) &
+ mo_coef(n+2,j) * ao_overlap(l,n+2) &
+ mo_coef(n+3,j) * ao_overlap(l,n+3) )
enddo
enddo
do n = lmax+1, ao_num
!DIR$ VECTOR ALIGNED
do l = 1, ao_num
mo_overlap(i,j) = mo_overlap(i,j) + mo_coef(n,j) * mo_coef(l,i) * ao_overlap(l,n)
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
BEGIN_PROVIDER [ double precision, mo_mono_elec_integral,(mo_tot_num_align,mo_tot_num)]
implicit none
integer :: i,j,n,l

13
src/MonoInts/pot_ao_ints.irp.f
View File

@ -14,19 +14,6 @@
ao_nucl_elec_integral = 0.d0
! -- Dummy run to provide everything called by NAI_pol_mult
!power_A(:)= ao_power(1,:)
!A_center(:) = nucl_coord(1,:)
!B_center(:) = nucl_coord(1,:)+0.5d0
!C_center(:) = 0.5d0*(B_center(:)-A_center(:))
!alpha = 1.d0
!beta = 0.1d0
!Z = 1.d0
!n_pt_in = n_pt_max_integrals
!c = NAI_pol_mult(A_center,B_center,power_A,power_A,alpha,beta,C_center,n_pt_in)
!c = NAI_pol_mult(A_center,A_center,power_A,power_A,alpha,beta,C_center,n_pt_in)
!c = NAI_pol_mult(A_center,A_center,power_A,power_A,alpha,beta,A_center,n_pt_in)
! --
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &

3
src/NEEDED_MODULES
View File

@ -1,2 +1 @@
AOs BiInts Bitmask CASSD DDCI Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full Utils
AOs BiInts Bitmask CASSD DDCI Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MonoInts MOs Nuclei Output Perturbation Selectors_full Utils Good_states Molden Primitive_basis Loc_MOs

42
src/Nuclei/nuclei.irp.f
View File

@ -196,4 +196,44 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ]
'Nuclear repulsion energy')
END_PROVIDER
BEGIN_PROVIDER [ character*(128), element_name, (36)]
BEGIN_DOC
! Array of the name of element, sorted by nuclear charge (integer)
END_DOC
element_name(1) = 'H'
element_name(2) = 'He'
element_name(3) = 'Li'
element_name(4) = 'Be'
element_name(5) = 'B'
element_name(6) = 'C'
element_name(7) = 'N'
element_name(8) = 'O'
element_name(9) = 'F'
element_name(10) = 'Ne'
element_name(11) = 'Na'
element_name(12) = 'Mg'
element_name(13) = 'Al'
element_name(14) = 'Si'
element_name(15) = 'P'
element_name(16) = 'S'
element_name(17) = 'Cl'
element_name(18) = 'Ar'
element_name(19) = 'K'
element_name(20) = 'Ca'
element_name(21) = 'Sc'
element_name(22) = 'Ti'
element_name(23) = 'V'
element_name(24) = 'Cr'
element_name(25) = 'Mn'
element_name(26) = 'Fe'
element_name(27) = 'Co'
element_name(28) = 'Ni'
element_name(29) = 'Cu'
element_name(30) = 'Zn'
element_name(31) = 'Ga'
element_name(32) = 'Ge'
element_name(33) = 'As'
element_name(34) = 'Se'
element_name(35) = 'Br'
element_name(36) = 'Kr'
END_PROVIDER

63
src/Utils/integration.irp.f
View File

@ -118,6 +118,69 @@ subroutine give_explicit_poly_and_gaussian(P_new,P_center,p,fact_k,iorder,alpha,
end
subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
BEGIN_DOC
! Transforms the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3)
! exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) exp(-(r-Nucl_center)^2 gama
!
! into
! fact_k * [ sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x ] exp (- p (x-P_center(1))^2 )
! * [ sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y ] exp (- p (y-P_center(2))^2 )
! * [ sum (l_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
END_DOC
implicit none
include 'constants.F'
integer, intent(in) :: dim
integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
double precision, intent(in) :: alpha, beta, gama ! exponents
double precision, intent(in) :: A_center(3) ! A center
double precision, intent(in) :: B_center (3) ! B center
double precision, intent(in) :: Nucl_center(3) ! B center
double precision, intent(out) :: P_center(3) ! new center
double precision, intent(out) :: p ! new exponent
double precision, intent(out) :: fact_k ! constant factor
double precision, intent(out) :: P_new(0:max_dim,3)! polynomial
integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials
double precision :: P_center_tmp(3) ! new center
double precision :: p_tmp ! new exponent
double precision :: fact_k_tmp,fact_k_bis ! constant factor
double precision :: P_new_tmp(0:max_dim,3)! polynomial
integer :: i,j
double precision :: binom_func
! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp
call give_explicit_poly_and_gaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim)
! Then you create the new gaussian from the product of the new one per the Nuclei one
call gaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center)
fact_k = fact_k_bis * fact_k_tmp
! Then you build the coefficient of the new polynom
do i = 0, iorder(1)
P_new(i,1) = 0.d0
do j = i,iorder(1)
P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i)
enddo
enddo
do i = 0, iorder(2)
P_new(i,2) = 0.d0
do j = i,iorder(2)
P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i)
enddo
enddo
do i = 0, iorder(3)
P_new(i,3) = 0.d0
do j = i,iorder(3)
P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i)
enddo
enddo
end
subroutine gaussian_product(a,xa,b,xb,k,p,xp)
implicit none
BEGIN_DOC