diff --git a/src/IntPak/NormCoeff.f90 b/src/IntPak/NormCoeff.f90 deleted file mode 100644 index 9e6cabf..0000000 --- a/src/IntPak/NormCoeff.f90 +++ /dev/null @@ -1,29 +0,0 @@ -function NormCoeff(alpha,a) - - implicit none - -! Input variables - - double precision,intent(in) :: alpha - integer,intent(in) :: a(3) - -! local variable - double precision :: pi,dfa(3),dfac - integer :: atot - -! Output variable - double precision NormCoeff - - pi = 4d0*atan(1d0) - atot = a(1) + a(2) + a(3) - - dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1))) - dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2))) - dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3))) - - - NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot - NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3)) - NormCoeff = sqrt(NormCoeff) - -end function NormCoeff diff --git a/src/IntPak/elements.f90 b/src/IntPak/elements.f90 deleted file mode 100644 index 22953dc..0000000 --- a/src/IntPak/elements.f90 +++ /dev/null @@ -1,170 +0,0 @@ -function element_number(element_name) - - implicit none - - integer,parameter :: nelement_max = 103 - character(len=2),intent(in) :: element_name - integer :: element_number - character(len=2),parameter :: element_list(nelement_max) = & - (/' H', 'He', & ! 2 - 'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10 - 'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18 - ' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36 - 'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54 - 'Cs','Ba', & ! 56 - 'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70 - 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86 - 'Fr','Ra', & ! 88 - 'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102 - 'Lr' & ! 103 - /) - -!===== - integer :: ielement -!===== - - ielement=1 - do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) ) - if( ielement == nelement_max ) then - write(*,'(a,a)') ' Input symbol ',element_name - write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements' - write(*,*) '!!! element symbol not understood !!!' - stop - endif - ielement = ielement + 1 - enddo - - element_number = ielement - -end function element_number - -function element_core(zval,zatom) - implicit none - double precision,intent(in) :: zval - double precision,intent(in) :: zatom - integer :: element_core -!===== - - ! - ! If zval /= zatom, this is certainly an effective core potential - ! and no core states should be frozen. - if( ABS(zval - zatom) > 1d0-3 ) then - element_core = 0 - else - - if( zval <= 4.00001d0 ) then ! up to Be - element_core = 0 - else if( zval <= 12.00001d0 ) then ! up to Mg - element_core = 1 - else if( zval <= 30.00001d0 ) then ! up to Ca - element_core = 5 - else if( zval <= 48.00001d0 ) then ! up to Sr - element_core = 9 - else - write(*,*) '!!! not imlemented in element_core !!!' - stop - endif - - endif - - -end function element_core - -function element_covalent_radius(zatom) - -! Return covalent radius of an atom - - implicit none - include 'parameters.h' - - integer,intent(in) :: zatom - double precision :: element_covalent_radius - - ! - ! Data from Cambridge Structural Database - ! http://en.wikipedia.org/wiki/Covalent_radius - ! - ! Values are first given in picometer - ! They will be converted in bohr later on - select case(zatom) - case( 1) - element_covalent_radius = 31. - case( 2) - element_covalent_radius = 28. - case( 3) - element_covalent_radius = 128. - case( 4) - element_covalent_radius = 96. - case( 5) - element_covalent_radius = 84. - case( 6) - element_covalent_radius = 73. - case( 7) - element_covalent_radius = 71. - case( 8) - element_covalent_radius = 66. - case( 9) - element_covalent_radius = 57. - case(10) ! Ne. - element_covalent_radius = 58. - case(11) - element_covalent_radius = 166. - case(12) - element_covalent_radius = 141. - case(13) - element_covalent_radius = 121. - case(14) - element_covalent_radius = 111. - case(15) - element_covalent_radius = 107. - case(16) - element_covalent_radius = 105. - case(17) - element_covalent_radius = 102. - case(18) ! Ar. - element_covalent_radius = 106. - case(19) - element_covalent_radius = 203. - case(20) - element_covalent_radius = 176. - case(21) - element_covalent_radius = 170. - case(22) - element_covalent_radius = 160. - case(23) - element_covalent_radius = 153. - case(24) - element_covalent_radius = 139. - case(25) - element_covalent_radius = 145. - case(26) - element_covalent_radius = 145. - case(27) - element_covalent_radius = 140. - case(28) - element_covalent_radius = 124. - case(29) - element_covalent_radius = 132. - case(30) - element_covalent_radius = 122. - case(31) - element_covalent_radius = 120. - case(32) - element_covalent_radius = 119. - case(34) - element_covalent_radius = 120. - case(35) - element_covalent_radius = 120. - case(36) ! Kr. - element_covalent_radius = 116. - case default - write(*,*) '!!! covalent radius not available !!!' - stop - end select - - ! pm to bohr conversion - element_covalent_radius = element_covalent_radius*pmtoau - - -end function element_covalent_radius - diff --git a/src/IntPak/read_auxiliary_basis.f90 b/src/IntPak/read_auxiliary_basis.f90 deleted file mode 100644 index cf4f6ed..0000000 --- a/src/IntPak/read_auxiliary_basis.f90 +++ /dev/null @@ -1,176 +0,0 @@ -subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, & - TotAngMomShell,KShell,DShell,ExpShell) - -! Read auxiliary basis set information - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: NAtoms - double precision,intent(in) :: XYZAtoms(NAtoms,3) - -! Local variables - - integer :: nShAt,iAt - integer :: i,j,k - character :: shelltype - -! Output variables - - integer,intent(out) :: nShell - double precision,intent(out) :: CenterShell(maxShell,3) - integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - -!------------------------------------------------------------------------ -! Primary basis set information -!------------------------------------------------------------------------ - -! Open file with basis set specification - - open(unit=2,file='input/basis') - -! Read basis information - - write(*,'(A28)') 'Gaussian basis set' - write(*,'(A28)') '------------------' - - nShell = 0 - do i=1,NAtoms - read(2,*) iAt,nShAt - write(*,'(A28,1X,I16)') 'Atom n. ',iAt - write(*,'(A28,1X,I16)') 'number of shells ',nShAt - write(*,'(A28)') '------------------' - -! Basis function centers - - do j=1,nShAt - nShell = nShell + 1 - do k=1,3 - CenterShell(nShell,k) = XYZAtoms(iAt,k) - enddo - -! Shell type and contraction degree - - read(2,*) shelltype,KShell(nShell) - if(shelltype == "S") then - TotAngMomShell(nShell) = 0 - write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - elseif(shelltype == "P") then - TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - elseif(shelltype == "D") then - TotAngMomShell(nShell) = 2 - write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - elseif(shelltype == "F") then - TotAngMomShell(nShell) = 3 - write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - elseif(shelltype == "G") then - TotAngMomShell(nShell) = 4 - write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - elseif(shelltype == "H") then - TotAngMomShell(nShell) = 5 - write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - elseif(shelltype == "I") then - TotAngMomShell(nShell) = 6 - write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif - -! Read exponents and contraction coefficients - - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(2,*) ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo - write(*,'(A28)') '------------------' - enddo - -! Total number of shells - - write(*,'(A28,1X,I16)') 'Number of shells in OBS',nShell - write(*,'(A28)') '------------------' - write(*,*) - -! Close file with basis set specification - - close(unit=2) - -!------------------------------------------------------------------------ -! Auxiliary basis set information -!------------------------------------------------------------------------ - -! Open file with auxilairy basis specification - - open(unit=3,file='input/auxbasis') - -! Read basis information - - write(*,'(A28)') 'Auxiliary basis set' - write(*,'(A28)') '-------------------' - - do i=1,NAtoms - read(3,*) iAt,nShAt - write(*,'(A28,1X,I16)') 'Atom n. ',iAt - write(*,'(A28,1X,I16)') 'number of shells ',nShAt - write(*,'(A28)') '------------------' - -! Basis function centers - - do j=1,nShAt - nShell = nShell + 1 - do k=1,3 - CenterShell(nShell,k) = XYZAtoms(iAt,k) - enddo - -! Shell type and contraction degree - - read(3,*) shelltype,KShell(nShell) - if(shelltype == "S") then - TotAngMomShell(nShell) = 0 - write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - elseif(shelltype == "P") then - TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - elseif(shelltype == "D") then - TotAngMomShell(nShell) = 2 - write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - elseif(shelltype == "F") then - TotAngMomShell(nShell) = 3 - write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - elseif(shelltype == "G") then - TotAngMomShell(nShell) = 4 - write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - elseif(shelltype == "H") then - TotAngMomShell(nShell) = 5 - write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - elseif(shelltype == "I") then - TotAngMomShell(nShell) = 6 - write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif - -! Read exponents and contraction coefficients - - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(3,*) ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo - write(*,'(A28)') '------------------' - enddo - -! Total number of shells - - write(*,'(A28,1X,I16)') 'Number of shells in ABS',nShell - write(*,'(A28)') '------------------' - write(*,*) - -! Close file with basis set specification - - close(unit=3) - -end subroutine read_auxiliary_basis diff --git a/src/IntPak/read_basis.f90 b/src/IntPak/read_basis.f90 deleted file mode 100644 index 075e1d6..0000000 --- a/src/IntPak/read_basis.f90 +++ /dev/null @@ -1,118 +0,0 @@ -subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) - -! Read basis set information - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nNuc,nO - double precision,intent(in) :: rNuc(nNuc,ncart) - -! Local variables - - integer :: nShAt,iNuc,iShell - integer :: i,j,k - character :: shelltype - -! Output variables - - integer,intent(out) :: nShell,nBas - double precision,intent(out) :: CenterShell(maxShell,ncart) - integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - integer,intent(out) :: nV - -!------------------------------------------------------------------------ -! Primary basis set information -!------------------------------------------------------------------------ - -! Open file with basis set specification - - open(unit=2,file='input/basis') - -! Read basis information - - write(*,'(A28)') 'Gaussian basis set' - write(*,'(A28)') '------------------' - - nShell = 0 - do i=1,nNuc - read(2,*) iNuc,nShAt - write(*,'(A28,1X,I16)') 'Atom n. ',iNuc - write(*,'(A28,1X,I16)') 'number of shells ',nShAt - write(*,'(A28)') '------------------' - -! Basis function centers - - do j=1,nShAt - nShell = nShell + 1 - do k=1,ncart - CenterShell(nShell,k) = rNuc(iNuc,k) - enddo - -! Shell type and contraction degree - - read(2,*) shelltype,KShell(nShell) - if(shelltype == "S") then - TotAngMomShell(nShell) = 0 - write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - elseif(shelltype == "P") then - TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - elseif(shelltype == "D") then - TotAngMomShell(nShell) = 2 - write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - elseif(shelltype == "F") then - TotAngMomShell(nShell) = 3 - write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - elseif(shelltype == "G") then - TotAngMomShell(nShell) = 4 - write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - elseif(shelltype == "H") then - TotAngMomShell(nShell) = 5 - write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - elseif(shelltype == "I") then - TotAngMomShell(nShell) = 6 - write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif - -! Read exponents and contraction coefficients - - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(2,*) ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo - write(*,'(A28)') '------------------' - enddo - -! Total number of shells - - write(*,'(A28,1X,I16)') 'Number of shells',nShell - write(*,'(A28)') '------------------' - write(*,*) - -! Close file with basis set specification - - close(unit=2) - -! Calculate number of basis functions - - nBas = 0 - do iShell=1,nShell - nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 - enddo - - write(*,'(A28)') '------------------' - write(*,'(A28,1X,I16)') 'Number of basis functions',NBas - write(*,'(A28)') '------------------' - write(*,*) - -! Number of virtual orbitals - - nV = nBas - nO - -end subroutine read_basis diff --git a/src/IntPak/read_geometry.f90 b/src/IntPak/read_geometry.f90 deleted file mode 100644 index 60c60b8..0000000 --- a/src/IntPak/read_geometry.f90 +++ /dev/null @@ -1,68 +0,0 @@ -subroutine read_geometry(nAt,ZNuc,rA,ENuc) - -! Read molecular geometry - - implicit none - - include 'parameters.h' - -! Ouput variables - - integer,intent(in) :: nAt - -! Local variables - - integer :: i,j - double precision :: RAB - character(len=2) :: El - integer,external :: element_number - -! Ouput variables - - double precision,intent(out) :: ZNuc(NAt),rA(nAt,ncart),ENuc - -! Open file with geometry specification - - open(unit=1,file='input/molecule') - -! Read geometry - - read(1,*) - read(1,*) - read(1,*) - - do i=1,nAt - read(1,*) El,rA(i,1),rA(i,2),rA(i,3) - ZNuc(i) = element_number(El) - enddo - -! Compute nuclear repulsion energy - - ENuc = 0 - - do i=1,nAt-1 - do j=i+1,nAt - RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2 - ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB) - enddo - enddo - -! Close file with geometry specification - close(unit=1) - -! Print geometry - write(*,'(A28)') '------------------' - write(*,'(A28)') 'Molecular geometry' - write(*,'(A28)') '------------------' - do i=1,NAt - write(*,'(A28,1X,I16)') 'Atom n. ',i - write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i) - write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,ncart) - enddo - write(*,*) - write(*,'(A28)') '------------------' - write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc - write(*,'(A28)') '------------------' - write(*,*) - -end subroutine read_geometry diff --git a/src/IntPak/read_molecule.f90 b/src/IntPak/read_molecule.f90 deleted file mode 100644 index 4d0ef7c..0000000 --- a/src/IntPak/read_molecule.f90 +++ /dev/null @@ -1,39 +0,0 @@ -subroutine read_molecule(nAt,nEl,nO,nC,nR) - -! Read number of atoms nAt and number of electrons nEl - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(out) :: nAt,nEl,nO,nC,nR - -! Open file with geometry specification - - open(unit=1,file='input/molecule') - -! Read number of atoms and number of electrons - - read(1,*) - read(1,*) nAt,nEl,nC,nR - - nO = nEl/2 - -! Print results - - write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of atoms',nAt - write(*,'(A28)') '----------------------' - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of electrons',nEl - write(*,'(A28)') '----------------------' - write(*,*) - -! Close file with geometry specification - - close(unit=1) - -end subroutine read_molecule diff --git a/src/IntPak/utils.f90 b/src/IntPak/utils.f90 deleted file mode 100644 index 17d2103..0000000 --- a/src/IntPak/utils.f90 +++ /dev/null @@ -1,385 +0,0 @@ -!------------------------------------------------------------------------ - -function KroneckerDelta(i,j) result(delta) - -! Kronecker Delta - - implicit none - -! Input variables - - integer,intent(in) :: i,j - - -! Output variables - - integer :: delta - - if(i == j) then - delta = 1 - else - delta = 0 - endif - -end function KroneckerDelta - -!------------------------------------------------------------------------ - -subroutine matout(m,n,A) - -! Print the MxN array A - - implicit none - - integer,parameter :: ncol = 5 - double precision,parameter :: small = 1d-10 - integer,intent(in) :: m,n - double precision,intent(in) :: A(m,n) - double precision :: B(ncol) - integer :: ilower,iupper,num,i,j - - do ilower=1,n,ncol - iupper = min(ilower + ncol - 1,n) - num = iupper - ilower + 1 - write(*,'(3X,10(7X,I6))') (j,j=ilower,iupper) - do i=1,m - do j=ilower,iupper - B(j-ilower+1) = A(i,j) - enddo - do j=1,num - if(abs(B(j)) < small) B(j) = 0d0 - enddo - write(*,'(I7,10F15.8)') i,(B(j),j=1,num) - enddo - enddo - -end subroutine matout - -!------------------------------------------------------------------------ - -subroutine CalcTrAB(n,A,B,Tr) - -! Calculate the trace of the square matrix A.B - - implicit none - -! Input variables - - integer,intent(in) :: n - double precision,intent(in) :: A(n,n),B(n,n) - -! Local variables - - integer :: i,j - -! Output variables - - double precision,intent(out) :: Tr - - Tr = 0d0 - do i=1,n - do j=1,n - Tr = Tr + A(i,j)*B(j,i) - enddo - enddo - -end subroutine CalcTrAB - -!------------------------------------------------------------------------ - -function EpsilonSwitch(i,j) result(delta) - -! Epsilon function - - implicit none - -! Input variables - - integer,intent(in) :: i,j - integer :: delta - - if(i <= j) then - delta = 1 - else - delta = -1 - endif - -end function EpsilonSwitch - -!------------------------------------------------------------------------ - -function KappaCross(i,j,k) result(kappa) - -! kappa(i,j,k) = eps(i,j) delta(i,k) - eps(k,i) delta(i,j) - - implicit none - -! Input variables - - integer,intent(in) :: i,j,k - -! Local variables - - integer :: EpsilonSwitch,KroneckerDelta - double precision :: kappa - - kappa = dble(EpsilonSwitch(i,j)*KroneckerDelta(i,k) - EpsilonSwitch(k,i)*KroneckerDelta(i,j)) - -end function KappaCross - -!------------------------------------------------------------------------ - -subroutine CalcInv3(a,det) - -! Calculate the inverse and the determinant of a 3x3 matrix - - implicit none - - double precision,intent(inout) :: a(3,3) - double precision, intent(inout) :: det - double precision :: b(3,3) - integer :: i,j - - det = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) & - - a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) & - + a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1)) - - do i=1,3 - b(i,1) = a(i,1) - b(i,2) = a(i,2) - b(i,3) = a(i,3) - enddo - - a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2) - a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3) - a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1) - - a(1,2) = b(1,3)*b(3,2) - b(1,2)*b(3,3) - a(2,2) = b(1,1)*b(3,3) - b(1,3)*b(3,1) - a(3,2) = b(1,2)*b(3,1) - b(1,1)*b(3,2) - - a(1,3) = b(1,2)*b(2,3) - b(1,3)*b(2,2) - a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3) - a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1) - - do i=1,3 - do j=1,3 - a(i,j) = a(i,j)/det - enddo - enddo - -end subroutine CalcInv3 - -!------------------------------------------------------------------------ - -subroutine CalcInv4(a,det) - - implicit none - - double precision,intent(inout) :: a(4,4) - double precision,intent(inout) :: det - double precision :: b(4,4) - integer :: i,j - - det = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) & - -a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) & - +a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) & - - a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) & - -a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) & - +a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) & - + a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) & - -a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) & - +a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) & - - a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) & - -a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) & - +a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) - do i=1,4 - b(1,i) = a(1,i) - b(2,i) = a(2,i) - b(3,i) = a(3,i) - b(4,i) = a(4,i) - enddo - - a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)) - a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)) - a(3,1) = b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) - a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) - - a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2)) - a(2,2) = b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1)) - a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) - a(4,2) = b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1)) - - a(1,3) = b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2)) - a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1)) - a(3,3) = b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1)) - a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1)) - - a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2)) - a(2,4) = b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*b(3,1)) - a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) - a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1)) - - do i=1,4 - do j=1,4 - a(i,j) = a(i,j)/det - enddo - enddo - -end subroutine CalcInv4 - - -!double precision function binom(i,j) -! implicit none -! integer,intent(in) :: i,j -! double precision :: logfact -! integer, save :: ifirst -! double precision, save :: memo(0:15,0:15) -! integer :: k,l -! if (ifirst == 0) then -! ifirst = 1 -! do k=0,15 -! do l=0,15 -! memo(k,l) = dexp( logfact(k)-logfact(l)-logfact(k-l) ) -! enddo -! enddo -! endif -! if ( (i<=15).and.(j<=15) ) then -! binom = memo(i,j) -! else -! binom = dexp( logfact(i)-logfact(j)-logfact(i-j) ) -! endif -!end -! -!double precision function fact(n) -! implicit none -! integer :: n -! double precision, save :: memo(1:100) -! integer, save :: memomax = 1 -! -! if (n<=memomax) then -! if (n<2) then -! fact = 1.d0 -! else -! fact = memo(n) -! endif -! return -! endif -! -! integer :: i -! memo(1) = 1.d0 -! do i=memomax+1,min(n,100) -! memo(i) = memo(i-1)*dble(i) -! enddo -! memomax = min(n,100) -! double precision :: logfact -! fact = dexp(logfact(n)) -!end function -! -!double precision function logfact(n) -! implicit none -! integer :: n -! double precision, save :: memo(1:100) -! integer, save :: memomax = 1 -! -! if (n<=memomax) then -! if (n<2) then -! logfact = 0.d0 -! else -! logfact = memo(n) -! endif -! return -! endif -! -! integer :: i -! memo(1) = 0.d0 -! do i=memomax+1,min(n,100) -! memo(i) = memo(i-1)+dlog(dble(i)) -! enddo -! memomax = min(n,100) -! logfact = memo(memomax) -! do i=101,n -! logfact += dlog(dble(i)) -! enddo -!end function -! -!double precision function dble_fact(n) -! implicit none -! integer :: n -! double precision :: dble_fact_even, dble_fact_odd -! -! dble_fact = 1.d0 -! -! if(n.lt.0) return -! -! if(iand(n,1).eq.0)then -! dble_fact = dble_fact_even(n) -! else -! dble_fact= dble_fact_odd(n) -! endif -! -!end function -! -!double precision function dble_fact_even(n) result(fact2) -! implicit none -! integer :: n,k -! double precision, save :: memo(0:100) -! integer, save :: memomax = 0 -! double precision :: prod -! -! -! if (n <= memomax) then -! if (n < 2) then -! fact2 = 1.d0 -! else -! fact2 = memo(n) -! endif -! return -! endif -! -! integer :: i -! memo(0)=1.d0 -! memo(1)=1.d0 -! do i=memomax+2,min(n,100),2 -! memo(i) = memo(i-2)* dble(i) -! enddo -! memomax = min(n,100) -! fact2 = memo(memomax) -! -! if (n > 100) then -! double precision :: dble_logfact -! fact2 = dexp(dble_logfact(n)) -! endif -! -!end function -! -!double precision function dble_fact_odd(n) result(fact2) -! implicit none -! integer :: n -! double precision, save :: memo(1:100) -! integer, save :: memomax = 1 -! -! if (n<=memomax) then -! if (n<3) then -! fact2 = 1.d0 -! else -! fact2 = memo(n) -! endif -! return -! endif -! -! integer :: i -! memo(1) = 1.d0 -! do i=memomax+2,min(n,99),2 -! memo(i) = memo(i-2)* dble(i) -! enddo -! memomax = min(n,99) -! fact2 = memo(memomax) -! -! if (n > 99) then -! double precision :: dble_logfact -! fact2 = dexp(dble_logfact(n)) -! endif -! -!end function - diff --git a/src/QuAcK/NormCoeff.f90 b/src/QuAcK/NormCoeff.f90 deleted file mode 100644 index 9e6cabf..0000000 --- a/src/QuAcK/NormCoeff.f90 +++ /dev/null @@ -1,29 +0,0 @@ -function NormCoeff(alpha,a) - - implicit none - -! Input variables - - double precision,intent(in) :: alpha - integer,intent(in) :: a(3) - -! local variable - double precision :: pi,dfa(3),dfac - integer :: atot - -! Output variable - double precision NormCoeff - - pi = 4d0*atan(1d0) - atot = a(1) + a(2) + a(3) - - dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1))) - dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2))) - dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3))) - - - NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot - NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3)) - NormCoeff = sqrt(NormCoeff) - -end function NormCoeff diff --git a/src/QuAcK/elements.f90 b/src/QuAcK/elements.f90 deleted file mode 100644 index 22953dc..0000000 --- a/src/QuAcK/elements.f90 +++ /dev/null @@ -1,170 +0,0 @@ -function element_number(element_name) - - implicit none - - integer,parameter :: nelement_max = 103 - character(len=2),intent(in) :: element_name - integer :: element_number - character(len=2),parameter :: element_list(nelement_max) = & - (/' H', 'He', & ! 2 - 'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10 - 'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18 - ' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36 - 'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54 - 'Cs','Ba', & ! 56 - 'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70 - 'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86 - 'Fr','Ra', & ! 88 - 'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102 - 'Lr' & ! 103 - /) - -!===== - integer :: ielement -!===== - - ielement=1 - do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) ) - if( ielement == nelement_max ) then - write(*,'(a,a)') ' Input symbol ',element_name - write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements' - write(*,*) '!!! element symbol not understood !!!' - stop - endif - ielement = ielement + 1 - enddo - - element_number = ielement - -end function element_number - -function element_core(zval,zatom) - implicit none - double precision,intent(in) :: zval - double precision,intent(in) :: zatom - integer :: element_core -!===== - - ! - ! If zval /= zatom, this is certainly an effective core potential - ! and no core states should be frozen. - if( ABS(zval - zatom) > 1d0-3 ) then - element_core = 0 - else - - if( zval <= 4.00001d0 ) then ! up to Be - element_core = 0 - else if( zval <= 12.00001d0 ) then ! up to Mg - element_core = 1 - else if( zval <= 30.00001d0 ) then ! up to Ca - element_core = 5 - else if( zval <= 48.00001d0 ) then ! up to Sr - element_core = 9 - else - write(*,*) '!!! not imlemented in element_core !!!' - stop - endif - - endif - - -end function element_core - -function element_covalent_radius(zatom) - -! Return covalent radius of an atom - - implicit none - include 'parameters.h' - - integer,intent(in) :: zatom - double precision :: element_covalent_radius - - ! - ! Data from Cambridge Structural Database - ! http://en.wikipedia.org/wiki/Covalent_radius - ! - ! Values are first given in picometer - ! They will be converted in bohr later on - select case(zatom) - case( 1) - element_covalent_radius = 31. - case( 2) - element_covalent_radius = 28. - case( 3) - element_covalent_radius = 128. - case( 4) - element_covalent_radius = 96. - case( 5) - element_covalent_radius = 84. - case( 6) - element_covalent_radius = 73. - case( 7) - element_covalent_radius = 71. - case( 8) - element_covalent_radius = 66. - case( 9) - element_covalent_radius = 57. - case(10) ! Ne. - element_covalent_radius = 58. - case(11) - element_covalent_radius = 166. - case(12) - element_covalent_radius = 141. - case(13) - element_covalent_radius = 121. - case(14) - element_covalent_radius = 111. - case(15) - element_covalent_radius = 107. - case(16) - element_covalent_radius = 105. - case(17) - element_covalent_radius = 102. - case(18) ! Ar. - element_covalent_radius = 106. - case(19) - element_covalent_radius = 203. - case(20) - element_covalent_radius = 176. - case(21) - element_covalent_radius = 170. - case(22) - element_covalent_radius = 160. - case(23) - element_covalent_radius = 153. - case(24) - element_covalent_radius = 139. - case(25) - element_covalent_radius = 145. - case(26) - element_covalent_radius = 145. - case(27) - element_covalent_radius = 140. - case(28) - element_covalent_radius = 124. - case(29) - element_covalent_radius = 132. - case(30) - element_covalent_radius = 122. - case(31) - element_covalent_radius = 120. - case(32) - element_covalent_radius = 119. - case(34) - element_covalent_radius = 120. - case(35) - element_covalent_radius = 120. - case(36) ! Kr. - element_covalent_radius = 116. - case default - write(*,*) '!!! covalent radius not available !!!' - stop - end select - - ! pm to bohr conversion - element_covalent_radius = element_covalent_radius*pmtoau - - -end function element_covalent_radius - diff --git a/src/QuAcK/read_basis.f90 b/src/QuAcK/read_basis.f90 deleted file mode 100644 index 2f3768c..0000000 --- a/src/QuAcK/read_basis.f90 +++ /dev/null @@ -1,118 +0,0 @@ -subroutine read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell) - -! Read basis set information - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nAt,nO - double precision,intent(in) :: rAt(nAt,ncart) - -! Local variables - - integer :: nShAt,iAt,iShell - integer :: i,j,k - character :: shelltype - -! Output variables - - integer,intent(out) :: nShell,nBas - double precision,intent(out) :: CenterShell(maxShell,ncart) - integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell) - double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK) - integer,intent(out) :: nV - -!------------------------------------------------------------------------ -! Primary basis set information -!------------------------------------------------------------------------ - -! Open file with basis set specification - - open(unit=2,file='input/basis') - -! Read basis information - - write(*,'(A28)') 'Gaussian basis set' - write(*,'(A28)') '------------------' - - nShell = 0 - do i=1,nAt - read(2,*) iAt,nShAt - write(*,'(A28,1X,I16)') 'Atom n. ',iAt - write(*,'(A28,1X,I16)') 'number of shells ',nShAt - write(*,'(A28)') '------------------' - -! Basis function centers - - do j=1,nShAt - nShell = nShell + 1 - do k=1,ncart - CenterShell(nShell,k) = rAt(iAt,k) - enddo - -! Shell type and contraction degree - - read(2,*) shelltype,KShell(nShell) - if(shelltype == "S") then - TotAngMomShell(nShell) = 0 - write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell) - elseif(shelltype == "P") then - TotAngMomShell(nShell) = 1 - write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell) - elseif(shelltype == "D") then - TotAngMomShell(nShell) = 2 - write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell) - elseif(shelltype == "F") then - TotAngMomShell(nShell) = 3 - write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell) - elseif(shelltype == "G") then - TotAngMomShell(nShell) = 4 - write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell) - elseif(shelltype == "H") then - TotAngMomShell(nShell) = 5 - write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell) - elseif(shelltype == "I") then - TotAngMomShell(nShell) = 6 - write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell) - endif - -! Read exponents and contraction coefficients - - write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction' - do k=1,Kshell(nShell) - read(2,*) ExpShell(nShell,k),DShell(nShell,k) - write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k) - enddo - enddo - write(*,'(A28)') '------------------' - enddo - -! Total number of shells - - write(*,'(A28,1X,I16)') 'Number of shells',nShell - write(*,'(A28)') '------------------' - write(*,*) - -! Close file with basis set specification - - close(unit=2) - -! Calculate number of basis functions - - nBas = 0 - do iShell=1,nShell - nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2 - enddo - - write(*,'(A28)') '------------------' - write(*,'(A28,1X,I16)') 'Number of basis functions',NBas - write(*,'(A28)') '------------------' - write(*,*) - -! Number of virtual orbitals - - nV = nBas - nO - -end subroutine read_basis diff --git a/src/QuAcK/read_geometry.f90 b/src/QuAcK/read_geometry.f90 deleted file mode 100644 index a3a4cf9..0000000 --- a/src/QuAcK/read_geometry.f90 +++ /dev/null @@ -1,68 +0,0 @@ -subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc) - -! Read molecular geometry - - implicit none - - include 'parameters.h' - -! Ouput variables - - integer,intent(in) :: nNuc - -! Local variables - - integer :: i,j - double precision :: RAB - character(len=2) :: El - integer,external :: element_number - -! Ouput variables - - double precision,intent(out) :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc - -! Open file with geometry specification - - open(unit=1,file='input/molecule') - -! Read geometry - - read(1,*) - read(1,*) - read(1,*) - - do i=1,nNuc - read(1,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3) - ZNuc(i) = element_number(El) - enddo - -! Compute nuclear repulsion energy - - ENuc = 0 - - do i=1,nNuc-1 - do j=i+1,nNuc - RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(j,3))**2 - ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB) - enddo - enddo - -! Close file with geometry specification - close(unit=1) - -! Print geometry - write(*,'(A28)') '------------------' - write(*,'(A28)') 'Molecular geometry' - write(*,'(A28)') '------------------' - do i=1,nNuc - write(*,'(A28,1X,I16)') 'Atom n. ',i - write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i) - write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rNuc(i,j),j=1,ncart) - enddo - write(*,*) - write(*,'(A28)') '------------------' - write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc - write(*,'(A28)') '------------------' - write(*,*) - -end subroutine read_geometry diff --git a/src/QuAcK/read_integrals.f90 b/src/QuAcK/read_integrals.f90 deleted file mode 100644 index e6f9dd6..0000000 --- a/src/QuAcK/read_integrals.f90 +++ /dev/null @@ -1,119 +0,0 @@ -subroutine read_integrals(nBas,S,T,V,Hc,G) - -! Read one- and two-electron integrals from files - - implicit none - -! Input variables - - integer,intent(in) :: nBas - -! Local variables - - logical :: debug - integer :: mu,nu,la,si - double precision :: Ov,Kin,Nuc,ERI - double precision :: scale - -! Output variables - - double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas) - -! Open file with integrals - - debug = .false. - - scale = 1d0 - - open(unit=8 ,file='int/Ov.dat') - open(unit=9 ,file='int/Kin.dat') - open(unit=10,file='int/Nuc.dat') - open(unit=11,file='int/ERI.dat') - -! Read overlap integrals - - S = 0d0 - do - read(8,*,end=8) mu,nu,Ov - S(mu,nu) = Ov - enddo - 8 close(unit=8) - -! Read kinetic integrals - - T = 0d0 - do - read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Kin/scale**2 - enddo - 9 close(unit=9) - -! Read nuclear integrals - - V = 0d0 - do - read(10,*,end=10) mu,nu,Nuc - V(mu,nu) = Nuc - enddo - 10 close(unit=10) - -! Define core Hamiltonian - - Hc = T + V - -! Read nuclear integrals - - G = 0d0 - do - read(11,*,end=11) mu,nu,la,si,ERI - - ERI = ERI/scale -! <12|34> - G(mu,nu,la,si) = ERI -! <32|14> - G(la,nu,mu,si) = ERI -! <14|32> - G(mu,si,la,nu) = ERI -! <34|12> - G(la,si,mu,nu) = ERI -! <41|23> - G(si,mu,nu,la) = ERI -! <23|41> - G(nu,la,si,mu) = ERI -! <21|43> - G(nu,mu,si,la) = ERI -! <43|21> - G(si,la,nu,mu) = ERI - enddo - 11 close(unit=11) - - -! Print results - if(debug) then - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Overlap integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,S) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Kinetic integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,T) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Nuclear integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,V) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Electron repulsion integrals' - write(*,'(A28)') '----------------------' - do la=1,nBas - do si=1,nBas - call matout(nBas,nBas,G(1,1,la,si)) - enddo - enddo - write(*,*) - endif - -end subroutine read_integrals diff --git a/src/QuAcK/read_molecule.f90 b/src/QuAcK/read_molecule.f90 deleted file mode 100644 index 5d8886e..0000000 --- a/src/QuAcK/read_molecule.f90 +++ /dev/null @@ -1,64 +0,0 @@ -subroutine read_molecule(nNuc,nEl,nO,nC,nR) - -! Read number of atoms and number of electrons - - implicit none - - include 'parameters.h' - -! Input variables - - integer,intent(out) :: nNuc - integer,intent(out) :: nEl(nspin) - integer,intent(out) :: nO(nspin) - integer,intent(out) :: nC(nspin) - integer,intent(out) :: nR(nspin) - -! Local variables - - integer :: nCore - integer :: nRyd - -! Open file with geometry specification - - open(unit=1,file='input/molecule') - -! Read number of atoms and number of electrons - - read(1,*) - read(1,*) nNuc,nEl(1),nEl(2),nCore,nRyd - - if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then - - print*, 'The number of core and Rydberg electrons must be even!' - stop - - end if - - nO(:) = nEl(:) - nC(:) = nCore/2 - nR(:) = nRyd/2 - -! Print results - - write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of atoms',nNuc - write(*,'(A28)') '----------------------' - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nEl(1) - write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nEl(2) - write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nEl(:)) - write(*,'(A28)') '----------------------' - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28,1X,I16)') 'Number of core electrons',sum(nC(:)) - write(*,'(A28,1X,I16)') 'Number of Rydberg electrons',sum(nR(:)) - write(*,'(A28)') '----------------------' - write(*,*) - -! Close file with geometry specification - - close(unit=1) - -end subroutine read_molecule diff --git a/src/QuAcK/utils.f90 b/src/QuAcK/utils.f90 deleted file mode 100644 index 8669f34..0000000 --- a/src/QuAcK/utils.f90 +++ /dev/null @@ -1,339 +0,0 @@ -!------------------------------------------------------------------------ -function Kronecker_delta(i,j) result(delta) - -! Kronecker Delta - - implicit none - -! Input variables - - integer,intent(in) :: i,j - -! Output variables - - double precision :: delta - - if(i == j) then - delta = 1d0 - else - delta = 0d0 - endif - -end function Kronecker_delta - -!------------------------------------------------------------------------ -subroutine matout(m,n,A) - -! Print the MxN array A - - implicit none - - integer,parameter :: ncol = 5 - double precision,parameter :: small = 1d-10 - integer,intent(in) :: m,n - double precision,intent(in) :: A(m,n) - double precision :: B(ncol) - integer :: ilower,iupper,num,i,j - - do ilower=1,n,ncol - iupper = min(ilower + ncol - 1,n) - num = iupper - ilower + 1 - write(*,'(3X,10(9X,I6))') (j,j=ilower,iupper) - do i=1,m - do j=ilower,iupper - B(j-ilower+1) = A(i,j) - enddo - do j=1,num - if(abs(B(j)) < small) B(j) = 0d0 - enddo - write(*,'(I7,10F15.8)') i,(B(j),j=1,num) - enddo - enddo - -end subroutine matout - -!------------------------------------------------------------------------ -subroutine trace_vector(n,v,Tr) - -! Calculate the trace of the vector v of length n -!!! Please use the intrinsic fortran sum() !!! - - implicit none - -! Input variables - - integer,intent(in) :: n - double precision,intent(in) :: v(n) - -! Local variables - - integer :: i - -! Output variables - - double precision,intent(out) :: Tr - - Tr = 0d0 - do i=1,n - Tr = Tr + v(i) - enddo - -end subroutine trace_vector - -!------------------------------------------------------------------------ -function trace_matrix(n,A) result(Tr) - -! Calculate the trace of the square matrix A - - implicit none - -! Input variables - - integer,intent(in) :: n - double precision,intent(in) :: A(n,n) - -! Local variables - - integer :: i - -! Output variables - - double precision :: Tr - - Tr = 0d0 - do i=1,n - Tr = Tr + A(i,i) - enddo - -end function trace_matrix - -!------------------------------------------------------------------------ -subroutine compute_error(nData,Mean,Var,Error) - -! Calculate the statistical error - - implicit none - -! Input variables - - double precision,intent(in) :: nData,Mean(3) - -! Output variables - - double precision,intent(out) :: Error(3) - double precision,intent(inout):: Var(3) - - Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0)) - -end subroutine compute_error - -!------------------------------------------------------------------------ -subroutine identity_matrix(N,A) - -! Set the matrix A to the identity matrix - - implicit none - -! Input variables - - integer,intent(in) :: N - -! Local viaruabkes - - integer :: i - -! Output variables - - double precision,intent(out) :: A(N,N) - - A = 0d0 - - do i=1,N - A(i,i) = 1d0 - enddo - -end subroutine identity_matrix - -!------------------------------------------------------------------------ -subroutine prepend(N,M,A,b) - -! Prepend the vector b of size N into the matrix A of size NxM - - implicit none - -! Input variables - - integer,intent(in) :: N,M - double precision,intent(in) :: b(N) - -! Local viaruabkes - - integer :: i,j - -! Output variables - - double precision,intent(out) :: A(N,M) - - -! print*,'b in append' -! call matout(N,1,b) - - do i=1,N - do j=M-1,1,-1 - A(i,j+1) = A(i,j) - enddo - A(i,1) = b(i) - enddo - -end subroutine prepend - -!------------------------------------------------------------------------ -subroutine append(N,M,A,b) - -! Append the vector b of size N into the matrix A of size NxM - - implicit none - -! Input variables - - integer,intent(in) :: N,M - double precision,intent(in) :: b(N) - -! Local viaruabkes - - integer :: i,j - -! Output variables - - double precision,intent(out) :: A(N,M) - - do i=1,N - do j=2,M - A(i,j-1) = A(i,j) - enddo - A(i,M) = b(i) - enddo - -end subroutine append - -!------------------------------------------------------------------------ -subroutine AtDA(N,A,D,B) - -! Perform B = At.D.A where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),D(N) - -! Local viaruabkes - - integer :: i,j,k - -! Output variables - - double precision,intent(out) :: B(N,N) - - B = 0d0 - - do i=1,N - do j=1,N - do k=1,N - B(i,k) = B(i,k) + A(j,i)*D(j)*A(j,k) - enddo - enddo - enddo - -end subroutine AtDA - -!------------------------------------------------------------------------ -subroutine ADAt(N,A,D,B) - -! Perform B = A.D.At where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),D(N) - -! Local viaruabkes - - integer :: i,j,k - -! Output variables - - double precision,intent(out) :: B(N,N) - - B = 0d0 - - do i=1,N - do j=1,N - do k=1,N - B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j) - enddo - enddo - enddo - -end subroutine ADAt -!------------------------------------------------------------------------ -subroutine DA(N,D,A) - -! Perform A <- D.A where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - - integer,intent(in) :: N - integer :: i,j,k - double precision,intent(in) :: D(N) - double precision,intent(inout):: A(N,N) - - do i=1,N - do j=1,N - A(i,j) = D(i)*A(i,j) - enddo - enddo - -end subroutine DA - -!------------------------------------------------------------------------ -subroutine AD(N,A,D) - -! Perform A <- A.D where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - - integer,intent(in) :: N - integer :: i,j,k - double precision,intent(in) :: D(N) - double precision,intent(inout):: A(N,N) - - do i=1,N - do j=1,N - A(i,j) = A(i,j)*D(j) - enddo - enddo - -end subroutine AD - -!------------------------------------------------------------------------ -subroutine print_warning(message) - -! Print warning - - implicit none - - character(len=*),intent(in) :: message - - write(*,*) message - -end subroutine print_warning - - diff --git a/src/QuAcK/wrap_lapack.f90 b/src/QuAcK/wrap_lapack.f90 deleted file mode 100644 index 6c29ab7..0000000 --- a/src/QuAcK/wrap_lapack.f90 +++ /dev/null @@ -1,207 +0,0 @@ -!subroutine eigenvalues_non_symmetric_matrix(N,A,e) -! -!! Diagonalize a square matrix -! -! implicit none -! -!! Input variables -! -! integer,intent(in) :: N -! double precision,intent(inout):: A(N,N) -! double precision,intent(out) :: e(N) -! -!! Local variables -! -! integer :: lwork,info -! double precision,allocatable :: work(:) -! -!! Memory allocation -! -! allocate(eRe(N),eIm(N),work(3*N)) -! lwork = size(work) -! -! call DGEEV('N','N',N,A,N, eRe, eIm, 0d0,1, VR,LDVR, WORK, LWORK, INFO ) -! -! if(info /= 0) then -! print*,'Problem in diagonalize_matrix (dseev)!!' -! stop -! endif -! -!end subroutine eigenvalues_non_symmetric_matrix - -subroutine diagonalize_matrix(N,A,e) - -! Diagonalize a square matrix - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(inout):: A(N,N) - double precision,intent(out) :: e(N) - -! Local variables - - integer :: lwork,info - double precision,allocatable :: work(:) - -! Memory allocation - - allocate(work(3*N)) - lwork = size(work) - - call dsyev('V','U',N,A,N,e,work,lwork,info) - - if(info /= 0) then - print*,'Problem in diagonalize_matrix (dsyev)!!' - endif - -end subroutine diagonalize_matrix - -subroutine svd(N,A,U,D,Vt) - - ! Compute A = U.D.Vt - ! Dimension of A is NxN - - implicit none - - integer, intent(in) :: N - double precision,intent(in) :: A(N,N) - double precision,intent(out) :: U(N,N) - double precision,intent(out) :: Vt(N,N) - double precision,intent(out) :: D(N) - double precision,allocatable :: work(:) - integer :: info,lwork - - double precision,allocatable :: scr(:,:) - - allocate (scr(N,N)) - - scr(:,:) = A(:,:) - - ! Find optimal size for temporary arrays - - allocate(work(1)) - - lwork = -1 - call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info) - lwork = int(work(1)) - - deallocate(work) - - allocate(work(lwork)) - - call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info) - - deallocate(work,scr) - - if (info /= 0) then - print *, info, ': SVD failed' - stop - endif - -end - -subroutine inverse_matrix(N,A,B) - -! Returns the inverse of the square matrix A in B - - implicit none - - integer,intent(in) :: N - double precision, intent(in) :: A(N,N) - double precision, intent(out) :: B(N,N) - - integer :: info,lwork - integer, allocatable :: ipiv(:) - double precision,allocatable :: work(:) - - allocate (ipiv(N),work(N*N)) - lwork = size(work) - - B(1:N,1:N) = A(1:N,1:N) - - call dgetrf(N,N,B,N,ipiv,info) - - if (info /= 0) then - - print*,info - stop 'error in inverse (dgetrf)!!' - - endif - - call dgetri(N,B,N,ipiv,work,lwork,info) - - if (info /= 0) then - - print *, info - stop 'error in inverse (dgetri)!!' - - endif - - deallocate(ipiv,work) - -end subroutine inverse_matrix - -subroutine linear_solve(N,A,b,x,rcond) - -! Solve the linear system A.x = b where A is a NxN matrix -! and x and x are vectors of size N - - implicit none - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),b(N),rcond - double precision,intent(out) :: x(N) - - integer :: info,lwork - double precision :: ferr,berr - integer,allocatable :: ipiv(:),iwork(:) - double precision,allocatable :: AF(:,:),work(:) - - lwork = 3*N - allocate(AF(N,N),ipiv(N),work(lwork),iwork(N)) - - call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info) - -! if (info /= 0) then - -! print *, info -! stop 'error in linear_solve (dsysvx)!!' - -! endif - -end subroutine linear_solve - -subroutine easy_linear_solve(N,A,b,x) - -! Solve the linear system A.x = b where A is a NxN matrix -! and x and x are vectors of size N - - implicit none - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),b(N) - double precision,intent(out) :: x(N) - - integer :: info,lwork - integer,allocatable :: ipiv(:) - double precision,allocatable :: work(:) - - allocate(ipiv(N),work(N*N)) - lwork = size(work) - - x = b - - call dsysv('U',N,1,A,N,ipiv,x,N,work,lwork,info) - - if (info /= 0) then - - print *, info - stop 'error in linear_solve (dsysv)!!' - - endif - -end subroutine easy_linear_solve - diff --git a/src/eDFT/utils.f90 b/src/eDFT/utils.f90 deleted file mode 100644 index 4142b10..0000000 --- a/src/eDFT/utils.f90 +++ /dev/null @@ -1,400 +0,0 @@ -!------------------------------------------------------------------------ -function Kronecker_delta(i,j) result(delta) - -! Kronecker Delta - - implicit none - -! Input variables - - integer,intent(in) :: i,j - -! Output variables - - double precision :: delta - - if(i == j) then - delta = 1d0 - else - delta = 0d0 - endif - -end function Kronecker_delta - -!------------------------------------------------------------------------ -subroutine matout(m,n,A) - -! Print the MxN array A - - implicit none - - integer,parameter :: ncol = 5 - double precision,parameter :: small = 1d-10 - integer,intent(in) :: m,n - double precision,intent(in) :: A(m,n) - double precision :: B(ncol) - integer :: ilower,iupper,num,i,j - - do ilower=1,n,ncol - iupper = min(ilower + ncol - 1,n) - num = iupper - ilower + 1 - write(*,'(3X,10(9X,I6))') (j,j=ilower,iupper) - do i=1,m - do j=ilower,iupper - B(j-ilower+1) = A(i,j) - enddo - do j=1,num - if(abs(B(j)) < small) B(j) = 0d0 - enddo - write(*,'(I7,10F15.8)') i,(B(j),j=1,num) - enddo - enddo - -end subroutine matout - -!------------------------------------------------------------------------ -subroutine trace_vector(n,v,Tr) - -! Calculate the trace of the vector v of length n -!!! Please use the intrinsic fortran sum() !!! - - implicit none - -! Input variables - - integer,intent(in) :: n - double precision,intent(in) :: v(n) - -! Local variables - - integer :: i - -! Output variables - - double precision,intent(out) :: Tr - - Tr = 0d0 - do i=1,n - Tr = Tr + v(i) - enddo - -end subroutine trace_vector - -!------------------------------------------------------------------------ -function trace_matrix(n,A) result(Tr) - -! Calculate the trace of the square matrix A - - implicit none - -! Input variables - - integer,intent(in) :: n - double precision,intent(in) :: A(n,n) - -! Local variables - - integer :: i - -! Output variables - - double precision :: Tr - - Tr = 0d0 - do i=1,n - Tr = Tr + A(i,i) - enddo - -end function trace_matrix - -!------------------------------------------------------------------------ -subroutine compute_error(nData,Mean,Var,Error) - -! Calculate the statistical error - - implicit none - -! Input variables - - double precision,intent(in) :: nData,Mean(3) - -! Output variables - - double precision,intent(out) :: Error(3) - double precision,intent(inout):: Var(3) - - Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0)) - -end subroutine compute_error - -!------------------------------------------------------------------------ -subroutine identity_matrix(N,A) - -! Set the matrix A to the identity matrix - - implicit none - -! Input variables - - integer,intent(in) :: N - -! Local viaruabkes - - integer :: i - -! Output variables - - double precision,intent(out) :: A(N,N) - - A = 0d0 - - do i=1,N - A(i,i) = 1d0 - enddo - -end subroutine identity_matrix - -!------------------------------------------------------------------------ -subroutine prepend(N,M,A,b) - -! Prepend the vector b of size N into the matrix A of size NxM - - implicit none - -! Input variables - - integer,intent(in) :: N,M - double precision,intent(in) :: b(N) - -! Local viaruabkes - - integer :: i,j - -! Output variables - - double precision,intent(out) :: A(N,M) - - -! print*,'b in append' -! call matout(N,1,b) - - do i=1,N - do j=M-1,1,-1 - A(i,j+1) = A(i,j) - enddo - A(i,1) = b(i) - enddo - -end subroutine prepend - -!------------------------------------------------------------------------ -subroutine append(N,M,A,b) - -! Append the vector b of size N into the matrix A of size NxM - - implicit none - -! Input variables - - integer,intent(in) :: N,M - double precision,intent(in) :: b(N) - -! Local viaruabkes - - integer :: i,j - -! Output variables - - double precision,intent(out) :: A(N,M) - - do i=1,N - do j=2,M - A(i,j-1) = A(i,j) - enddo - A(i,M) = b(i) - enddo - -end subroutine append - -!------------------------------------------------------------------------ -subroutine AtDA(N,A,D,B) - -! Perform B = At.D.A where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),D(N) - -! Local viaruabkes - - integer :: i,j,k - -! Output variables - - double precision,intent(out) :: B(N,N) - - B = 0d0 - - do i=1,N - do j=1,N - do k=1,N - B(i,k) = B(i,k) + A(j,i)*D(j)*A(j,k) - enddo - enddo - enddo - -end subroutine AtDA - -!------------------------------------------------------------------------ -subroutine ADAt(N,A,D,B) - -! Perform B = A.D.At where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),D(N) - -! Local viaruabkes - - integer :: i,j,k - -! Output variables - - double precision,intent(out) :: B(N,N) - - B = 0d0 - - do i=1,N - do j=1,N - do k=1,N - B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j) - enddo - enddo - enddo - -end subroutine ADAt -!------------------------------------------------------------------------ -subroutine DA(N,D,A) - -! Perform A <- D.A where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - - integer,intent(in) :: N - integer :: i,j,k - double precision,intent(in) :: D(N) - double precision,intent(inout):: A(N,N) - - do i=1,N - do j=1,N - A(i,j) = D(i)*A(i,j) - enddo - enddo - -end subroutine DA - -!------------------------------------------------------------------------ -subroutine AD(N,A,D) - -! Perform A <- A.D where A is a NxN matrix and D is a diagonal matrix given -! as a vector of length N - - implicit none - - integer,intent(in) :: N - integer :: i,j,k - double precision,intent(in) :: D(N) - double precision,intent(inout):: A(N,N) - - do i=1,N - do j=1,N - A(i,j) = A(i,j)*D(j) - enddo - enddo - -end subroutine AD - -!------------------------------------------------------------------------ - -subroutine print_warning(message) - -! Print warning - - implicit none - - character(len=*),intent(in) :: message - - write(*,*) message - -end subroutine print_warning - -!------------------------------------------------------------------------ -recursive function fac(n) result(fact) - - implicit none - integer :: fact - integer, intent(in) :: n - - if (n == 0) then - fact = 1 - else - fact = n * fac(n-1) - end if - -end function fac - -!------------------------------------------------------------------------ - -function dfac(n) result(fact) - - implicit none - double precision :: fact - integer, intent(in) :: n - integer :: fac - - fact = dble(fac(n)) - -end function dfac - -!------------------------------------------------------------------------ - -function NormCoeff(alpha,a) - -! Compute normalization coefficients for cartesian gaussians - - implicit none - -! Input variables - - double precision,intent(in) :: alpha - integer,intent(in) :: a(3) - -! local variable - double precision :: pi,dfa(3),dfac - integer :: atot - -! Output variable - double precision NormCoeff - - pi = 4d0*atan(1d0) - atot = a(1) + a(2) + a(3) - - dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1))) - dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2))) - dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3))) - - - NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot - NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3)) - NormCoeff = sqrt(NormCoeff) - -end function NormCoeff diff --git a/src/eDFT/wrap_lapack.f90 b/src/eDFT/wrap_lapack.f90 deleted file mode 100644 index 6c29ab7..0000000 --- a/src/eDFT/wrap_lapack.f90 +++ /dev/null @@ -1,207 +0,0 @@ -!subroutine eigenvalues_non_symmetric_matrix(N,A,e) -! -!! Diagonalize a square matrix -! -! implicit none -! -!! Input variables -! -! integer,intent(in) :: N -! double precision,intent(inout):: A(N,N) -! double precision,intent(out) :: e(N) -! -!! Local variables -! -! integer :: lwork,info -! double precision,allocatable :: work(:) -! -!! Memory allocation -! -! allocate(eRe(N),eIm(N),work(3*N)) -! lwork = size(work) -! -! call DGEEV('N','N',N,A,N, eRe, eIm, 0d0,1, VR,LDVR, WORK, LWORK, INFO ) -! -! if(info /= 0) then -! print*,'Problem in diagonalize_matrix (dseev)!!' -! stop -! endif -! -!end subroutine eigenvalues_non_symmetric_matrix - -subroutine diagonalize_matrix(N,A,e) - -! Diagonalize a square matrix - - implicit none - -! Input variables - - integer,intent(in) :: N - double precision,intent(inout):: A(N,N) - double precision,intent(out) :: e(N) - -! Local variables - - integer :: lwork,info - double precision,allocatable :: work(:) - -! Memory allocation - - allocate(work(3*N)) - lwork = size(work) - - call dsyev('V','U',N,A,N,e,work,lwork,info) - - if(info /= 0) then - print*,'Problem in diagonalize_matrix (dsyev)!!' - endif - -end subroutine diagonalize_matrix - -subroutine svd(N,A,U,D,Vt) - - ! Compute A = U.D.Vt - ! Dimension of A is NxN - - implicit none - - integer, intent(in) :: N - double precision,intent(in) :: A(N,N) - double precision,intent(out) :: U(N,N) - double precision,intent(out) :: Vt(N,N) - double precision,intent(out) :: D(N) - double precision,allocatable :: work(:) - integer :: info,lwork - - double precision,allocatable :: scr(:,:) - - allocate (scr(N,N)) - - scr(:,:) = A(:,:) - - ! Find optimal size for temporary arrays - - allocate(work(1)) - - lwork = -1 - call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info) - lwork = int(work(1)) - - deallocate(work) - - allocate(work(lwork)) - - call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info) - - deallocate(work,scr) - - if (info /= 0) then - print *, info, ': SVD failed' - stop - endif - -end - -subroutine inverse_matrix(N,A,B) - -! Returns the inverse of the square matrix A in B - - implicit none - - integer,intent(in) :: N - double precision, intent(in) :: A(N,N) - double precision, intent(out) :: B(N,N) - - integer :: info,lwork - integer, allocatable :: ipiv(:) - double precision,allocatable :: work(:) - - allocate (ipiv(N),work(N*N)) - lwork = size(work) - - B(1:N,1:N) = A(1:N,1:N) - - call dgetrf(N,N,B,N,ipiv,info) - - if (info /= 0) then - - print*,info - stop 'error in inverse (dgetrf)!!' - - endif - - call dgetri(N,B,N,ipiv,work,lwork,info) - - if (info /= 0) then - - print *, info - stop 'error in inverse (dgetri)!!' - - endif - - deallocate(ipiv,work) - -end subroutine inverse_matrix - -subroutine linear_solve(N,A,b,x,rcond) - -! Solve the linear system A.x = b where A is a NxN matrix -! and x and x are vectors of size N - - implicit none - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),b(N),rcond - double precision,intent(out) :: x(N) - - integer :: info,lwork - double precision :: ferr,berr - integer,allocatable :: ipiv(:),iwork(:) - double precision,allocatable :: AF(:,:),work(:) - - lwork = 3*N - allocate(AF(N,N),ipiv(N),work(lwork),iwork(N)) - - call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info) - -! if (info /= 0) then - -! print *, info -! stop 'error in linear_solve (dsysvx)!!' - -! endif - -end subroutine linear_solve - -subroutine easy_linear_solve(N,A,b,x) - -! Solve the linear system A.x = b where A is a NxN matrix -! and x and x are vectors of size N - - implicit none - - integer,intent(in) :: N - double precision,intent(in) :: A(N,N),b(N) - double precision,intent(out) :: x(N) - - integer :: info,lwork - integer,allocatable :: ipiv(:) - double precision,allocatable :: work(:) - - allocate(ipiv(N),work(N*N)) - lwork = size(work) - - x = b - - call dsysv('U',N,1,A,N,ipiv,x,N,work,lwork,info) - - if (info /= 0) then - - print *, info - stop 'error in linear_solve (dsysv)!!' - - endif - -end subroutine easy_linear_solve -