Skip to content

Commit

Permalink
Merge pull request merzlab#334 from Madu86/333-wrong-d-and-f-mo-coeff…
Browse files Browse the repository at this point in the history 
…icient-order-in-quick-generated-molden-files

Fix for issue 333
  • Loading branch information
@agoetz
agoetz authored Mar 15, 2024
2 parents 83ef917 + a9584de commit fd1ccc4
Show file tree
Hide file tree
Showing 3 changed files with 132 additions and 30 deletions.
6 changes: 6 additions & 0 deletions src/basis.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -444,6 +444,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
else
quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),0,0,0)
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=BB(k)

quick_basis%gcexpo(k,Ninitial)=AA(k)

Expand Down Expand Up @@ -486,6 +487,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
else
quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),1,0,0)
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=BB(k)
!quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),1,0,0)
quick_basis%gcexpo(k,Ninitial)=AA(k)
if(quick_basis%gcexpomin(jshell).gt.AA(k))quick_basis%gcexpomin(jshell)=AA(k)
Expand Down Expand Up @@ -529,6 +531,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
else
quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),0,0,0)
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=BB(k)
quick_basis%gcexpo(k,Ninitial)=AA(k)
if(quick_basis%gcexpomin(jshell).gt.AA(k))quick_basis%gcexpomin(jshell)=AA(k)
enddo
Expand All @@ -548,6 +551,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
else
quick_basis%gccoeff(k,Ninitial)=CC(k)*xnorm(AA(k),1,0,0)
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=CC(k)
quick_basis%gcexpo(k,Ninitial)=AA(k)
if(quick_basis%gcexpomin(jshell).gt.AA(k))quick_basis%gcexpomin(jshell)=AA(k)
enddo
Expand Down Expand Up @@ -613,6 +617,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),quick_basis%KLMN(1,Ninitial), &
quick_basis%KLMN(2,Ninitial),quick_basis%KLMN(3,Ninitial))
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=BB(k)

quick_basis%gcexpo(k,Ninitial)=AA(k)
if(quick_basis%gcexpomin(jshell).gt.AA(k))quick_basis%gcexpomin(jshell)=AA(k)
Expand Down Expand Up @@ -707,6 +712,7 @@ subroutine readbasis(natomxiao,natomstart,natomfinal,nbasisstart,nbasisfinal,ier
quick_basis%gccoeff(k,Ninitial)=BB(k)*xnorm(AA(k),quick_basis%KLMN(1,Ninitial), &
quick_basis%KLMN(2,Ninitial),quick_basis%KLMN(3,Ninitial))
endif
quick_basis%unnorm_gccoeff(k,Ninitial)=BB(k)
quick_basis%gcexpo(k,Ninitial)=AA(k)
if(quick_basis%gcexpomin(jshell).gt.AA(k))quick_basis%gcexpomin(jshell)=AA(k)
enddo
Expand Down
6 changes: 6 additions & 0 deletions src/modules/quick_basis_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -98,6 +98,9 @@ module quick_basis_module
! normalized coeffecient
double precision, allocatable, dimension(:,:) :: gccoeff

! unnormalized contraction coefficients
double precision, allocatable, dimension(:,:) :: unnorm_gccoeff

! basis set factor
double precision, allocatable, dimension(:) :: cons

Expand Down Expand Up @@ -255,11 +258,13 @@ subroutine allocate_quick_basis(self,natom_arg,nshell_arg,nbasis_arg)

if(.not. allocated(self%gcexpo)) allocate(self%gcexpo(MAXPRIM,nbasis_arg))
if(.not. allocated(self%gccoeff)) allocate(self%gccoeff(MAXPRIM,nbasis_arg))
if(.not. allocated(self%unnorm_gccoeff)) allocate(self%unnorm_gccoeff(MAXPRIM,nbasis_arg))
if(.not. allocated(self%cons)) allocate(self%cons(nbasis_arg))
do i = 1, MAXPRIM
do j = 1, nbasis_arg
self%gcexpo( i, j) = 0.0
self%gccoeff( i, j) = 0.0
self%unnorm_gccoeff( i, j) = 0.0
enddo
enddo

Expand Down Expand Up @@ -303,6 +308,7 @@ subroutine deallocate_quick_basis(self)
if (allocated(self%cons)) deallocate(self%cons)
if (allocated(self%gcexpo)) deallocate(self%gcexpo)
if (allocated(self%gccoeff)) deallocate(self%gccoeff)
if (allocated(self%unnorm_gccoeff)) deallocate(self%unnorm_gccoeff)
if (allocated(self%KLMN)) deallocate(self%KLMN)

#if defined CUDA_MPIV || defined HIP_MPIV
Expand Down
150 changes: 120 additions & 30 deletions src/modules/quick_molden_module.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,9 @@ module quick_molden_module

! counter to keep track of number of snapshots
integer :: iexport_snapshot

! temporary vector for reorganizing mo coefficients
double precision, dimension(:), allocatable :: reord_mo_vec

end type quick_molden_type

Expand Down Expand Up @@ -97,59 +100,94 @@ subroutine write_coordinates(self, ierr)
! we need to do this because optimizers may return the geometry
! for the next step which may be stored in xyz
k = self%iexport_snapshot - 1
write(self%iMoldenFile, '(3(4x,F10.4))') (self%xyz_snapshots(j,i,k),j=1,3)
write(self%iMoldenFile, '(3(2x,F20.10))') (self%xyz_snapshots(j,i,k),j=1,3)
else
! if it's a single point calculation we can use xyz
! we can't use xyz_snapshots because they have not been populated
write(self%iMoldenFile, '(3(4x,F10.4))') (xyz(j,i),j=1,3)
write(self%iMoldenFile, '(3(2x,F20.10))') (xyz(j,i),j=1,3)
endif
enddo

end subroutine write_coordinates

subroutine write_basis_info(self, ierr)

use quick_basis_module, only: quick_basis, nshell, nbasis, aexp, dcoeff, ncontract
use quick_basis_module, only: quick_basis, nshell, nbasis, ncontract
use quick_molspec_module, only: natom
implicit none
type (quick_molden_type), intent(in) :: self
integer, intent(out) :: ierr
integer :: iatom, ishell, ibas, iprim, nprim, j
integer :: iatom, ishell, ibas, iprim, nprim, j, ishell_idx
logical :: print_gto
double precision :: val_gccoeff, xnorm

! write basis function information
write(self%iMoldenFile, '("[GTO] (AU)")')
do iatom=1, natom
write(self%iMoldenFile, '(2x, I5)') iatom

! s basis functions
do ishell=1, nshell
if(quick_basis%katom(ishell) .eq. iatom) then
nprim = quick_basis%kprim(ishell)
if(quick_basis%ktype(ishell) .eq. 1) then
write(self%iMoldenFile, '(2x, "s", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 3) then
write(self%iMoldenFile, '(2x, "p", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 4) then
write(self%iMoldenFile, '(2x, "sp", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 6) then
write(self%iMoldenFile, '(2x, "d", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 10) then
write(self%iMoldenFile, '(2x, "f", 4x, I2)') nprim
endif

if(quick_basis%ktype(ishell) .eq. 4) then
do iprim=1, nprim
write(self%iMoldenFile, '(2x, E14.8, 2x, E14.8, 2x, E14.8)') &
aexp(iprim, quick_basis%ksumtype(ishell)), (dcoeff(iprim,quick_basis%ksumtype(ishell)+j), j=0,1)
ishell_idx=quick_basis%ksumtype(ishell)
write(self%iMoldenFile, '(2E20.10)') &
quick_basis%gcexpo(iprim,ishell_idx), quick_basis%unnorm_gccoeff(iprim,ishell_idx)
enddo
endif
endif
enddo

else
! s, p basis functions of sp shell
do ishell=1, nshell
if(quick_basis%katom(ishell) .eq. iatom) then
nprim = quick_basis%kprim(ishell)
if(quick_basis%ktype(ishell) .eq. 4) then
write(self%iMoldenFile, '(2x, "s", 4x, I2)') nprim
do iprim=1, nprim
write(self%iMoldenFile, '(2x, E14.8, 2x, E14.8)') &
aexp(iprim, quick_basis%ksumtype(ishell)), dcoeff(iprim,quick_basis%ksumtype(ishell))
ishell_idx=quick_basis%ksumtype(ishell)
write(self%iMoldenFile, '(2E20.10)') &
quick_basis%gcexpo(iprim,ishell_idx), quick_basis%unnorm_gccoeff(iprim,ishell_idx)
enddo
write(self%iMoldenFile, '(2x, "p", 4x, I2)') nprim
do iprim=1, nprim
ishell_idx=quick_basis%ksumtype(ishell)
write(self%iMoldenFile, '(2E20.10)') &
quick_basis%gcexpo(iprim,ishell_idx), (quick_basis%unnorm_gccoeff(iprim,ishell_idx+1))
enddo
endif
endif
enddo

! p, d, and f basis functions
do ishell=1, nshell
if(quick_basis%katom(ishell) .eq. iatom) then
nprim = quick_basis%kprim(ishell)
print_gto=.false.
if(quick_basis%ktype(ishell) .eq. 3) then
print_gto=.true.
write(self%iMoldenFile, '(2x, "p", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 6) then
print_gto=.true.
write(self%iMoldenFile, '(2x, "d", 4x, I2)') nprim
elseif(quick_basis%ktype(ishell) .eq. 10) then
print_gto=.true.
write(self%iMoldenFile, '(2x, "f", 4x, I2)') nprim
endif

do iprim=1, nprim
if(print_gto) then
ishell_idx=quick_basis%ksumtype(ishell)
write(self%iMoldenFile, '(2E20.10)') &
quick_basis%gcexpo(iprim,ishell_idx), quick_basis%unnorm_gccoeff(iprim,ishell_idx)
endif
enddo
endif
enddo

write(self%iMoldenFile, '("")')
enddo

Expand All @@ -163,12 +201,14 @@ subroutine write_mo(self, ierr)
use quick_molspec_module, only: quick_molspec
use quick_method_module, only: quick_method
implicit none
type (quick_molden_type), intent(in) :: self
type (quick_molden_type), intent(inout) :: self
integer, intent(out) :: ierr
integer :: i, j, k, neleca, nelecb
character(len=5) :: lbl1
double precision :: occnum, occval

if(.not. allocated(self%reord_mo_vec)) allocate(self%reord_mo_vec(nbasis))

write(self%iMoldenFile, '("[MO]")')

if(.not.quick_method%unrst) then
Expand All @@ -189,16 +229,19 @@ subroutine write_mo(self, ierr)
endif

write(lbl1,'(I5)') i
write(self%iMoldenFile, '(A11)') "Sym= a"//trim(adjustl(lbl1))
write(self%iMoldenFile, '(2x, "Ene= ", E16.10)') quick_qm_struct%E(i)
write(self%iMoldenFile, '(A11)') " Sym= a"//trim(adjustl(lbl1))
write(self%iMoldenFile, '(" Ene= ", E20.10)') quick_qm_struct%E(i)
write(self%iMoldenFile, '(2x, "Spin= Alpha" )')

! write orbital occupation numbers
write(self%iMoldenFile, '(2x, "Occup= ", F10.4)') occnum
write(self%iMoldenFile, '(" Occup= ", F10.8)') occnum

! reorder molecular orbital coefficients
call reorder_mo_coeffs(quick_qm_struct%co, quick_basis%KLMN, nbasis, i, self%reord_mo_vec, ierr)

! write molecular orbital coefficients
! write molecular orbital coefficients
do j=1, nbasis
write(self%iMoldenFile, '(2x, I5, 2x, E16.10)') j, quick_qm_struct%co(j,i)
write(self%iMoldenFile, '(I4,F15.10)') j, self%reord_mo_vec(j)
enddo
enddo

Expand All @@ -212,20 +255,25 @@ subroutine write_mo(self, ierr)
endif

write(lbl1,'(I5)') i
write(self%iMoldenFile, '(A11)') "Sym= b"//trim(adjustl(lbl1))
write(self%iMoldenFile, '(2x, "Ene= ", E16.10)') quick_qm_struct%Eb(i)
write(self%iMoldenFile, '(A11)') " Sym= b"//trim(adjustl(lbl1))
write(self%iMoldenFile, '(" Ene= ",E20.10)') quick_qm_struct%Eb(i)
write(self%iMoldenFile, '(2x, "Spin= Beta" )')

! write orbital occupation numbers
write(self%iMoldenFile, '(2x, "Occup= ", F10.4)') occnum
write(self%iMoldenFile, '(" Occup= ", F10.8)') occnum

! reorder molecular orbital coefficients
call reorder_mo_coeffs(quick_qm_struct%cob, quick_basis%KLMN, nbasis, i, self%reord_mo_vec, ierr)

! write molecular orbital coefficients
do j=1, nbasis
write(self%iMoldenFile, '(2x, I5, 2x, E16.10)') j, quick_qm_struct%cob(j,i)
write(self%iMoldenFile, '(I4,F15.10)') j, self%reord_mo_vec(j)
enddo
enddo
endif

if(allocated(self%reord_mo_vec)) deallocate(self%reord_mo_vec)

end subroutine write_mo

subroutine write_scf(self, ierr)
Expand Down Expand Up @@ -341,4 +389,46 @@ subroutine finalize_molden(self, ierr)

end subroutine finalize_molden

subroutine reorder_mo_coeffs(co, KLMN, nbasis, i, reord_mo_vec, ierr)

implicit none
double precision, intent(in) :: co(:,:)
integer, intent(in) :: KLMN(:,:)
integer, intent (in) :: nbasis, i
double precision, intent(inout) :: reord_mo_vec(:)
integer, intent(out) :: ierr
integer :: j

j=1
do while (j <= nbasis)
reord_mo_vec(j)=co(j,i)
! order d functions. xx is the first, followed by yy, zz, xy, xz, yz
if(KLMN(1,j) .eq. 2) then
reord_mo_vec(j+1)=co(j+2,i)
reord_mo_vec(j+2)=co(j+5,i)
reord_mo_vec(j+3)=co(j+1,i)
reord_mo_vec(j+4)=co(j+3,i)
reord_mo_vec(j+5)=co(j+4,i)
j=j+5
endif

! order f functions. xxx is the first, followed by yyy, zzz, xyy, xxy, xxz, xzz, yzz, yyz, xyz
if(KLMN(1,j) .eq. 3) then
reord_mo_vec(j+1)=co(j+3,i)
reord_mo_vec(j+2)=co(j+9,i)
reord_mo_vec(j+3)=co(j+2,i)
reord_mo_vec(j+4)=co(j+1,i)
reord_mo_vec(j+5)=co(j+4,i)
reord_mo_vec(j+6)=co(j+7,i)
reord_mo_vec(j+7)=co(j+8,i)
reord_mo_vec(j+8)=co(j+6,i)
reord_mo_vec(j+9)=co(j+5,i)
j=j+9
endif

j=j+1
enddo

end subroutine reorder_mo_coeffs

end module quick_molden_module

0 comments on commit fd1ccc4

Please sign in to comment.