diff --git a/src/mo_basis/mos.irp.f b/src/mo_basis/mos.irp.f index 6cf76c1b..04386e6b 100644 --- a/src/mo_basis/mos.irp.f +++ b/src/mo_basis/mos.irp.f @@ -93,82 +93,6 @@ BEGIN_PROVIDER [ double precision, mo_coef, (ao_num,mo_num) ] endif END_PROVIDER -BEGIN_PROVIDER [ double precision, mo_coef_imag, (ao_num,mo_num) ] - implicit none - BEGIN_DOC - ! Molecular orbital coefficients on |AO| basis set - ! - ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO| - ! - ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) - END_DOC - integer :: i, j - double precision, allocatable :: buffer(:,:) - logical :: exists - PROVIDE ezfio_filename - - - if (mpi_master) then - ! Coefs - call ezfio_has_mo_basis_mo_coef_imag(exists) - endif - IRP_IF MPI_DEBUG - print *, irp_here, mpi_rank - call MPI_BARRIER(MPI_COMM_WORLD, ierr) - IRP_ENDIF - IRP_IF MPI - include 'mpif.h' - integer :: ierr - call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read mo_coef_imag with MPI' - endif - IRP_ENDIF - - if (exists) then - if (mpi_master) then - call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag) - write(*,*) 'Read mo_coef_imag' - endif - IRP_IF MPI - call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) - if (ierr /= MPI_SUCCESS) then - stop 'Unable to read mo_coef_imag with MPI' - endif - IRP_ENDIF - else - ! Orthonormalized AO basis - do i=1,mo_num - do j=1,ao_num - mo_coef_imag(j,i) = 0.d0 - enddo - enddo - endif -END_PROVIDER - -BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] - implicit none - BEGIN_DOC - ! Molecular orbital coefficients on |AO| basis set - ! - ! mo_coef_complex(i,j) = coefficient of the i-th |AO| on the jth |MO| - ! - ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) - END_DOC - integer :: i, j - double precision, allocatable :: buffer(:,:) - logical :: exists - PROVIDE ezfio_filename - - provide mo_coef mo_coef_imag - - do i=1,mo_num - do j=1,ao_num - mo_coef_complex(j,i) = dcmplx(mo_coef(j,i),mo_coef_imag(j,i)) - enddo - enddo -END_PROVIDER - BEGIN_PROVIDER [ double precision, mo_coef_in_ao_ortho_basis, (ao_num, mo_num) ] implicit none BEGIN_DOC @@ -303,35 +227,6 @@ subroutine ao_to_mo(A_ao,LDA_ao,A_mo,LDA_mo) deallocate(T) end -subroutine ao_to_mo_complex(A_ao,LDA_ao,A_mo,LDA_mo) - implicit none - BEGIN_DOC - ! Transform A from the AO basis to the MO basis - ! where A is complex in the AO basis - ! - ! Ct.A_ao.C - END_DOC - integer, intent(in) :: LDA_ao,LDA_mo - complex*16, intent(in) :: A_ao(LDA_ao,ao_num) - complex*16, intent(out) :: A_mo(LDA_mo,mo_num) - complex*16, allocatable :: T(:,:) - - allocate ( T(ao_num,mo_num) ) - !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T - - call zgemm('N','N', ao_num, mo_num, ao_num, & - (1.d0,0.d0), A_ao,LDA_ao, & - mo_coef_complex, size(mo_coef_complex,1), & - (0.d0,0.d0), T, size(T,1)) - - call zgemm('C','N', mo_num, mo_num, ao_num, & - (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1), & - T, ao_num, & - (0.d0,0.d0), A_mo, size(A_mo,1)) - - deallocate(T) -end - subroutine mix_mo_jk(j,k) implicit none @@ -347,28 +242,43 @@ subroutine mix_mo_jk(j,k) ! by convention, the '+' |MO| is in the lowest index (min(j,k)) ! by convention, the '-' |MO| is in the highest index (max(j,k)) END_DOC - double precision :: array_tmp(ao_num,2),dsqrt_2 if(j==k)then print*,'You want to mix two orbitals that are the same !' print*,'It does not make sense ... ' print*,'Stopping ...' stop endif - array_tmp = 0.d0 + double precision :: dsqrt_2 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 + if (is_periodic) then + complex*16 :: array_tmp_c(ao_num,2) + array_tmp_c = (0.d0,0.d0) + do i = 1, ao_num + array_tmp_c(i,1) = dsqrt_2 * (mo_coef_complex(i,j) + mo_coef_complex(i,k)) + array_tmp_c(i,2) = dsqrt_2 * (mo_coef_complex(i,j) - mo_coef_complex(i,k)) + enddo + do i = 1, ao_num + mo_coef_complex(i,i_plus) = array_tmp_c(i,1) + mo_coef_complex(i,i_minus) = array_tmp_c(i,2) + enddo + else + double precision :: array_tmp(ao_num,2) + array_tmp = 0.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 + do i = 1, ao_num + mo_coef(i,i_plus) = array_tmp(i,1) + mo_coef(i,i_minus) = array_tmp(i,2) + enddo + endif end + subroutine ao_ortho_cano_to_ao(A_ao,LDA_ao,A,LDA) implicit none BEGIN_DOC diff --git a/src/mo_basis/mos_complex.irp.f b/src/mo_basis/mos_complex.irp.f new file mode 100644 index 00000000..54f98ef2 --- /dev/null +++ b/src/mo_basis/mos_complex.irp.f @@ -0,0 +1,163 @@ +BEGIN_PROVIDER [ double precision, mo_coef_imag, (ao_num,mo_num) ] + implicit none + BEGIN_DOC + ! Molecular orbital coefficients on |AO| basis set + ! + ! mo_coef_imag(i,j) = coefficient of the i-th |AO| on the jth |MO| + ! + ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) + END_DOC + integer :: i, j + double precision, allocatable :: buffer(:,:) + logical :: exists + PROVIDE ezfio_filename + + + if (mpi_master) then + ! Coefs + call ezfio_has_mo_basis_mo_coef_imag(exists) + endif + IRP_IF MPI_DEBUG + print *, irp_here, mpi_rank + call MPI_BARRIER(MPI_COMM_WORLD, ierr) + IRP_ENDIF + IRP_IF MPI + include 'mpif.h' + integer :: ierr + call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_imag with MPI' + endif + IRP_ENDIF + + if (exists) then + if (mpi_master) then + call ezfio_get_mo_basis_mo_coef_imag(mo_coef_imag) + write(*,*) 'Read mo_coef_imag' + endif + IRP_IF MPI + call MPI_BCAST( mo_coef_imag, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read mo_coef_imag with MPI' + endif + IRP_ENDIF + else + ! Orthonormalized AO basis + do i=1,mo_num + do j=1,ao_num + mo_coef_imag(j,i) = 0.d0 + enddo + enddo + endif +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, mo_coef_complex, (ao_num,mo_num) ] + implicit none + BEGIN_DOC + ! Molecular orbital coefficients on |AO| basis set + ! + ! mo_coef_complex(i,j) = coefficient of the i-th |AO| on the jth |MO| + ! + ! mo_label : Label characterizing the |MOs| (local, canonical, natural, etc) + END_DOC + integer :: i, j + PROVIDE ezfio_filename + + provide mo_coef mo_coef_imag + + do i=1,mo_num + do j=1,ao_num + mo_coef_complex(j,i) = dcmplx(mo_coef(j,i),mo_coef_imag(j,i)) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ complex*16, mo_coef_in_ao_ortho_basis_complex, (ao_num, mo_num) ] + implicit none + BEGIN_DOC + ! |MO| coefficients in orthogonalized |AO| basis + ! + ! $C^{-1}.C_{mo}$ + END_DOC + call zgemm('N','N',ao_num,mo_num,ao_num,(1.d0,0.d0), & + ao_ortho_canonical_coef_inv_complex, size(ao_ortho_canonical_coef_inv_complex,1),& + mo_coef_complex, size(mo_coef_complex,1), (0.d0,0.d0), & + mo_coef_in_ao_ortho_basis_complex, size(mo_coef_in_ao_ortho_basis_complex,1)) + +END_PROVIDER + + BEGIN_PROVIDER [ complex*16, mo_coef_transp_complex, (mo_num,ao_num) ] +&BEGIN_PROVIDER [ complex*16, mo_coef_transp_complex_conjg, (mo_num,ao_num) ] + implicit none + BEGIN_DOC + ! |MO| coefficients on |AO| basis set + END_DOC + integer :: i, j + + do j=1,ao_num + do i=1,mo_num + mo_coef_transp_complex(i,j) = mo_coef_complex(j,i) + mo_coef_transp_complex_conjg(i,j) = dconjg(mo_coef_complex(j,i)) + enddo + enddo + +END_PROVIDER + +subroutine ao_to_mo_complex(A_ao,LDA_ao,A_mo,LDA_mo) + implicit none + BEGIN_DOC + ! Transform A from the AO basis to the MO basis + ! where A is complex in the AO basis + ! + ! Ct.A_ao.C + END_DOC + integer, intent(in) :: LDA_ao,LDA_mo + complex*16, intent(in) :: A_ao(LDA_ao,ao_num) + complex*16, intent(out) :: A_mo(LDA_mo,mo_num) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num,mo_num) ) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: T + + call zgemm('N','N', ao_num, mo_num, ao_num, & + (1.d0,0.d0), A_ao,LDA_ao, & + mo_coef_complex, size(mo_coef_complex,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('C','N', mo_num, mo_num, ao_num, & + (1.d0,0.d0), mo_coef_complex,size(mo_coef_complex,1), & + T, ao_num, & + (0.d0,0.d0), A_mo, size(A_mo,1)) + + deallocate(T) +end + + +subroutine ao_ortho_cano_to_ao_complex(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}.C^{\dagger-1}$ + END_DOC + integer, intent(in) :: LDA_ao,LDA + complex*16, intent(in) :: A_ao(LDA_ao,*) + complex*16, intent(out) :: A(LDA,*) + complex*16, allocatable :: T(:,:) + + allocate ( T(ao_num,ao_num) ) + + call zgemm('C','N', ao_num, ao_num, ao_num, & + (1.d0,0.d0), & + ao_ortho_canonical_coef_inv_complex, size(ao_ortho_canonical_coef_inv_complex,1),& + A_ao,size(A_ao,1), & + (0.d0,0.d0), T, size(T,1)) + + call zgemm('N','N', ao_num, ao_num, ao_num, (1.d0,0.d0), & + T, size(T,1), & + ao_ortho_canonical_coef_inv_complex,size(ao_ortho_canonical_coef_inv_complex,1),& + (0.d0,0.d0), A, size(A,1)) + + deallocate(T) +end +