From 0242bf93031a0bb39b16736748d9583e33fec255 Mon Sep 17 00:00:00 2001 From: Manu Date: Wed, 7 Jan 2015 16:48:03 +0100 Subject: [PATCH] Molden interface --- src/AOs/aos.irp.f | 250 +++++++++++++++++++++++++++++++- src/CIS_dressed/CIS_full.irp.f | 4 +- src/MOs/NEEDED_MODULES | 4 +- src/MOs/mo_overlap.irp.f | 35 +++++ src/MOs/mos.irp.f | 119 +++++++++++++++ src/Molden/print_mo.irp.f | 131 +++++++++++++++++ src/MonoInts/ao_mono_ints.irp.f | 121 ---------------- src/MonoInts/mo_mono_ints.irp.f | 34 ----- src/MonoInts/pot_ao_ints.irp.f | 13 -- src/NEEDED_MODULES | 3 +- src/Nuclei/nuclei.irp.f | 42 +++++- src/Utils/integration.irp.f | 63 ++++++++ 12 files changed, 640 insertions(+), 179 deletions(-) create mode 100644 src/MOs/mo_overlap.irp.f create mode 100644 src/Molden/print_mo.irp.f diff --git a/src/AOs/aos.irp.f b/src/AOs/aos.irp.f index 613fdabd..1947867f 100644 --- a/src/AOs/aos.irp.f +++ b/src/AOs/aos.irp.f @@ -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 diff --git a/src/CIS_dressed/CIS_full.irp.f b/src/CIS_dressed/CIS_full.irp.f index 6be7da09..03194f7a 100644 --- a/src/CIS_dressed/CIS_full.irp.f +++ b/src/CIS_dressed/CIS_full.irp.f @@ -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) diff --git a/src/MOs/NEEDED_MODULES b/src/MOs/NEEDED_MODULES index 6e90bab2..5ca73603 100644 --- a/src/MOs/NEEDED_MODULES +++ b/src/MOs/NEEDED_MODULES @@ -1 +1,3 @@ -AOs Ezfio_files Nuclei Output Utils Electrons +AOs Electrons Ezfio_files Nuclei Output Utils + + diff --git a/src/MOs/mo_overlap.irp.f b/src/MOs/mo_overlap.irp.f new file mode 100644 index 00000000..9b8b48f0 --- /dev/null +++ b/src/MOs/mo_overlap.irp.f @@ -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 + diff --git a/src/MOs/mos.irp.f b/src/MOs/mos.irp.f index f1158530..17ca4cfd 100644 --- a/src/MOs/mos.irp.f +++ b/src/MOs/mos.irp.f @@ -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 diff --git a/src/Molden/print_mo.irp.f b/src/Molden/print_mo.irp.f new file mode 100644 index 00000000..30b14d92 --- /dev/null +++ b/src/Molden/print_mo.irp.f @@ -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 diff --git a/src/MonoInts/ao_mono_ints.irp.f b/src/MonoInts/ao_mono_ints.irp.f index 9a135dda..84e46f61 100644 --- a/src/MonoInts/ao_mono_ints.irp.f +++ b/src/MonoInts/ao_mono_ints.irp.f @@ -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 diff --git a/src/MonoInts/mo_mono_ints.irp.f b/src/MonoInts/mo_mono_ints.irp.f index e466c9cc..8d8fed6a 100644 --- a/src/MonoInts/mo_mono_ints.irp.f +++ b/src/MonoInts/mo_mono_ints.irp.f @@ -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 diff --git a/src/MonoInts/pot_ao_ints.irp.f b/src/MonoInts/pot_ao_ints.irp.f index ef68b0cb..d5d646ac 100644 --- a/src/MonoInts/pot_ao_ints.irp.f +++ b/src/MonoInts/pot_ao_ints.irp.f @@ -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) & diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index 40d37829..9eb1816c 100644 --- a/src/NEEDED_MODULES +++ b/src/NEEDED_MODULES @@ -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 diff --git a/src/Nuclei/nuclei.irp.f b/src/Nuclei/nuclei.irp.f index 98303e61..9fdc9c4d 100644 --- a/src/Nuclei/nuclei.irp.f +++ b/src/Nuclei/nuclei.irp.f @@ -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 diff --git a/src/Utils/integration.irp.f b/src/Utils/integration.irp.f index e1a6a59e..89dc559d 100644 --- a/src/Utils/integration.irp.f +++ b/src/Utils/integration.irp.f @@ -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