From ca3a985d50125a813418ddf209256e1cb60b1d2e Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 20 Mar 2019 13:38:42 +0100 Subject: [PATCH] create utils --- src/IntPak/CalcBoysF.f90 | 6 +- src/IntPak/CalcOm.f90 | 2 +- src/IntPak/CalcOm3e.f90 | 2 +- src/IntPak/CalcOmERI.f90 | 2 +- src/IntPak/CalcOmErf.f90 | 2 +- src/IntPak/CalcOmNuc.f90 | 2 +- src/IntPak/CalcOmYuk.f90 | 4 +- src/IntPak/Compute2eInt.f90 | 44 +-- src/IntPak/Compute3eInt.f90 | 56 +-- src/IntPak/Compute4eInt.f90 | 38 +- src/IntPak/ComputeKin.f90 | 22 +- src/IntPak/ComputeNuc.f90 | 24 +- src/IntPak/ComputeOv.f90 | 22 +- src/IntPak/FormVRR3e.f90 | 48 +-- src/IntPak/G2eInt.f90 | 30 +- src/IntPak/G3eInt.f90 | 18 +- src/IntPak/GF12Int.f90 | 2 +- src/IntPak/GenerateShell.f90 | 6 +- src/IntPak/HRR2e.f90 | 16 +- src/IntPak/HRR3e.f90 | 22 +- src/IntPak/HRRNuc.f90 | 8 +- src/IntPak/HRROv.f90 | 4 +- src/IntPak/IntPak.f90 | 26 +- src/IntPak/KinInt.f90 | 2 +- src/IntPak/NucInt.f90 | 8 +- src/IntPak/OvInt.f90 | 2 +- src/IntPak/S2eInt.f90 | 8 +- src/IntPak/S3eInt.f90 | 2 +- src/IntPak/VRR2e.f90 | 20 +- src/IntPak/VRR3e.f90 | 28 +- src/IntPak/VRRNuc.f90 | 8 +- src/IntPak/VRROv.f90 | 4 +- src/IntPak/obj/.gitignore | 1 + src/QuAcK/AO_values.f90 | 4 +- src/QuAcK/form_EOM_one_body.f90 | 44 +-- src/QuAcK/form_EOM_two_body.f90 | 112 +++--- src/QuAcK/obj/.gitingore | 1 + src/eDFT/AO_values_grid.f90 | 4 +- src/eDFT/obj/.gitingore | 1 + src/utils/elements.f90 | 170 +++++++++ src/utils/norm_coeff.f90 | 29 ++ src/utils/obj/.gitignore | 1 + src/utils/read_auxiliary_basis.f90 | 176 ++++++++++ src/utils/read_basis.f90 | 118 +++++++ src/utils/read_geometry.f90 | 68 ++++ src/utils/read_integrals.f90 | 119 +++++++ src/utils/read_molecule.f90 | 64 ++++ src/utils/utils.f90 | 533 +++++++++++++++++++++++++++++ src/utils/wrap_lapack.f90 | 207 +++++++++++ 49 files changed, 1816 insertions(+), 324 deletions(-) create mode 100644 src/IntPak/obj/.gitignore create mode 100644 src/QuAcK/obj/.gitingore create mode 100644 src/eDFT/obj/.gitingore create mode 100644 src/utils/elements.f90 create mode 100644 src/utils/norm_coeff.f90 create mode 100644 src/utils/obj/.gitignore create mode 100644 src/utils/read_auxiliary_basis.f90 create mode 100644 src/utils/read_basis.f90 create mode 100644 src/utils/read_geometry.f90 create mode 100644 src/utils/read_integrals.f90 create mode 100644 src/utils/read_molecule.f90 create mode 100644 src/utils/utils.f90 create mode 100644 src/utils/wrap_lapack.f90 diff --git a/src/IntPak/CalcBoysF.f90 b/src/IntPak/CalcBoysF.f90 index a9fc9ed..68fa43a 100644 --- a/src/IntPak/CalcBoysF.f90 +++ b/src/IntPak/CalcBoysF.f90 @@ -35,13 +35,13 @@ subroutine CalcBoysF(maxm,t,Fm) do m=0,maxm dm = dble(m) Fm(m) = 1d0/(2d0*dm+1d0) - enddo + end do else do m=0,maxm dm = dble(m) ! Fm(m) = gamma(dm+0.5d0)*gsl_sf_gamma_inc_P(dm+0.5d0,t)/(2d0*t**(dm+0.5d0)) Fm(m) = dgami(dm+0.5d0,t)/(2d0*t**(dm+0.5d0)) - enddo - endif + end do + end if end subroutine CalcBoysF diff --git a/src/IntPak/CalcOm.f90 b/src/IntPak/CalcOm.f90 index 0b88277..8b5c8e7 100644 --- a/src/IntPak/CalcOm.f90 +++ b/src/IntPak/CalcOm.f90 @@ -33,7 +33,7 @@ subroutine CalcOm(maxm,ExpPQi,NormPQSq,Om) do m=0,maxm dm =dble(m) Om(m) = (2d0/sqrt(pi))*(-1d0)**dm*(1d0/ExpPQi)**(dm+0.5d0)*Fm(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/CalcOm3e.f90 b/src/IntPak/CalcOm3e.f90 index aab4085..76a3525 100644 --- a/src/IntPak/CalcOm3e.f90 +++ b/src/IntPak/CalcOm3e.f90 @@ -37,7 +37,7 @@ subroutine CalcOm3e(maxm,delta0,delta1,Y1,Y0,Om) do m=0,maxm Om(m) = (2d0/sqrt(pi))*OG*sqrt(delta0/(delta1-delta0))*(delta1/(delta1-delta0))**m Om(m) = Om(m)*Fm(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/CalcOmERI.f90 b/src/IntPak/CalcOmERI.f90 index 6054b16..05c4abc 100644 --- a/src/IntPak/CalcOmERI.f90 +++ b/src/IntPak/CalcOmERI.f90 @@ -32,7 +32,7 @@ subroutine CalcOmERI(maxm,ExpY,NormYSq,Om) do m=0,maxm Om(m) = (2d0/sqrt(pi))*sqrt(ExpY)*Fm(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/CalcOmErf.f90 b/src/IntPak/CalcOmErf.f90 index 762cdb3..26adc72 100644 --- a/src/IntPak/CalcOmErf.f90 +++ b/src/IntPak/CalcOmErf.f90 @@ -32,7 +32,7 @@ subroutine CalcOmErf(maxm,ExpY,fG,NormYSq,Om) do m=0,maxm Om(m) = (2d0/sqrt(pi))*sqrt(fG)*(fG/ExpY)**m*Fm(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/CalcOmNuc.f90 b/src/IntPak/CalcOmNuc.f90 index dd953e9..5f5e730 100644 --- a/src/IntPak/CalcOmNuc.f90 +++ b/src/IntPak/CalcOmNuc.f90 @@ -33,7 +33,7 @@ subroutine CalcOmNuc(maxm,ExpPQi,NormPQSq,Om) do m=0,maxm dm =dble(m) Om(m) = (2d0/sqrt(pi))*(1d0/ExpPQi)**(dm+0.5d0)*Fm(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/CalcOmYuk.f90 b/src/IntPak/CalcOmYuk.f90 index 3836efa..739c10e 100644 --- a/src/IntPak/CalcOmYuk.f90 +++ b/src/IntPak/CalcOmYuk.f90 @@ -34,9 +34,9 @@ subroutine CalcOmYuk(maxm,ExpG,ExpY,fG,NormYSq,Om) Om(m) = 0d0 do k=0,m Om(m) = Om(m) + dbinom(m,k)*(ExpY/ExpG)**k*Fm(k) - enddo + end do Om(m) = (2d0/sqrt(pi))*sqrt(ExpY)*(fG/ExpG)*exp(-fG*NormYSq)*Om(m) - enddo + end do deallocate(Fm) diff --git a/src/IntPak/Compute2eInt.f90 b/src/IntPak/Compute2eInt.f90 index d208f11..5827d39 100644 --- a/src/IntPak/Compute2eInt.f90 +++ b/src/IntPak/Compute2eInt.f90 @@ -32,7 +32,7 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & integer,allocatable :: ShellFunctionB1(:,:),ShellFunctionB2(:,:) double precision :: ExpBra(2),ExpKet(2) double precision :: DBra(2),DKet(2) - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA1,iBasA2,iBasB1,iBasB2 integer :: iShA1,iShA2,iShB1,iShB2 @@ -110,7 +110,7 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & iFile = 24 open(unit=iFile,file='int/Erf.dat') - endif + end if !------------------------------------------------------------------------ ! Loops over shell A1 @@ -213,16 +213,16 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & do iKA1=1,KBra(1) ExpBra(1) = ExpShell(iShA1,iKA1) - DBra(1) = DShell(iShA1,iKA1)*NormCoeff(ExpBra(1),AngMomBra(1,1:3)) + DBra(1) = DShell(iShA1,iKA1)*norm_coeff(ExpBra(1),AngMomBra(1,1:3)) do iKA2=1,KBra(2) ExpBra(2) = ExpShell(iShA2,iKA2) - DBra(2) = DShell(iShA2,iKA2)*NormCoeff(ExpBra(2),AngMomBra(2,1:3)) + DBra(2) = DShell(iShA2,iKA2)*norm_coeff(ExpBra(2),AngMomBra(2,1:3)) do iKB1=1,KKet(1) ExpKet(1) = ExpShell(iShB1,iKB1) - DKet(1) = DShell(iShB1,iKB1)*NormCoeff(ExpKet(1),AngMomKet(1,1:3)) + DKet(1) = DShell(iShB1,iKB1)*norm_coeff(ExpKet(1),AngMomKet(1,1:3)) do iKB2=1,KKet(2) ExpKet(2) = ExpShell(iShB2,iKB2) - DKet(2) = DShell(iShB2,iKB2)*NormCoeff(ExpKet(2),AngMomKet(2,1:3)) + DKet(2) = DShell(iShB2,iKB2)*norm_coeff(ExpKet(2),AngMomKet(2,1:3)) call S2eInt(debug,iType,np2eInt,nSigp2eInt, & ExpS,KG,DG,ExpG, & @@ -232,10 +232,10 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & c2eInt = c2eInt + DBra(1)*DBra(2)*DKet(1)*DKet(2)*p2eInt - enddo - enddo - enddo - enddo + end do + end do + end do + end do call cpu_time(end_c2eInt) nc2eInt = nc2eInt + 1 @@ -252,7 +252,7 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & '(a1b1|a2b2) = ',c2eInt,iBasA1,iBasB1,iBasA2,iBasB2 - endif + end if else @@ -261,38 +261,38 @@ subroutine Compute2eInt(debug,chemist_notation,iType,nShell, & if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & ' = ',c2eInt,iBasA1,iBasA2,iBasB1,iBasB2 - endif + end if - endif - endif + end if + end if !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB2) - enddo + end do iBasB2 = 0 !------------------------------------------------------------------------ ! End loops over shell B2 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA2) - enddo + end do iBasA2 = 0 !------------------------------------------------------------------------ ! End loops over shell A2 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB1) - enddo + end do iBasB1 = 0 !------------------------------------------------------------------------ ! End loops over shell B1 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA1) - enddo + end do iBasA1 = 0 !------------------------------------------------------------------------ ! End loops over shell A1 diff --git a/src/IntPak/Compute3eInt.f90 b/src/IntPak/Compute3eInt.f90 index 008fdca..accf6de 100644 --- a/src/IntPak/Compute3eInt.f90 +++ b/src/IntPak/Compute3eInt.f90 @@ -32,7 +32,7 @@ subroutine Compute3eInt(debug,iType,nShell, & integer,allocatable :: ShellFunctionB1(:,:),ShellFunctionB2(:,:),ShellFunctionB3(:,:) double precision :: ExpBra(3),ExpKet(3) double precision :: DBra(3),DKet(3) - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3 integer :: iShA1,iShA2,iShA3,iShB1,iShB2,iShB3 @@ -80,7 +80,7 @@ subroutine Compute3eInt(debug,iType,nShell, & elseif(iType == 3) then iFile = 33 open(unit=iFile,file='int/3eInt_Type3.dat') - endif + end if !------------------------------------------------------------------------ ! Loops over shell A1 @@ -229,22 +229,22 @@ subroutine Compute3eInt(debug,iType,nShell, & do iKA1=1,KBra(1) ExpBra(1) = ExpShell(iShA1,iKA1) - DBra(1) = DShell(iShA1,iKA1)*NormCoeff(ExpBra(1),AngMomBra(1,1:3)) + DBra(1) = DShell(iShA1,iKA1)*norm_coeff(ExpBra(1),AngMomBra(1,1:3)) do iKA2=1,KBra(2) ExpBra(2) = ExpShell(iShA2,iKA2) - DBra(2) = DShell(iShA2,iKA2)*NormCoeff(ExpBra(2),AngMomBra(2,1:3)) + DBra(2) = DShell(iShA2,iKA2)*norm_coeff(ExpBra(2),AngMomBra(2,1:3)) do iKA3=1,KBra(3) ExpBra(3) = ExpShell(iShA3,iKA3) - DBra(3) = DShell(iShA3,iKA3)*NormCoeff(ExpBra(3),AngMomBra(3,1:3)) + DBra(3) = DShell(iShA3,iKA3)*norm_coeff(ExpBra(3),AngMomBra(3,1:3)) do iKB1=1,KKet(1) ExpKet(1) = ExpShell(iShB1,iKB1) - DKet(1) = DShell(iShB1,iKB1)*NormCoeff(ExpKet(1),AngMomKet(1,1:3)) + DKet(1) = DShell(iShB1,iKB1)*norm_coeff(ExpKet(1),AngMomKet(1,1:3)) do iKB2=1,KKet(2) ExpKet(2) = ExpShell(iShB2,iKB2) - DKet(2) = DShell(iShB2,iKB2)*NormCoeff(ExpKet(2),AngMomKet(2,1:3)) + DKet(2) = DShell(iShB2,iKB2)*norm_coeff(ExpKet(2),AngMomKet(2,1:3)) do iKB3=1,KKet(3) ExpKet(3) = ExpShell(iShB3,iKB3) - DKet(3) = DShell(iShB3,iKB3)*NormCoeff(ExpKet(3),AngMomKet(3,1:3)) + DKet(3) = DShell(iShB3,iKB3)*norm_coeff(ExpKet(3),AngMomKet(3,1:3)) call S3eInt(debug,iType,np3eInt,nSigp3eInt, & ExpS,KG,DG,ExpG, & @@ -254,12 +254,12 @@ subroutine Compute3eInt(debug,iType,nShell, & c3eInt = c3eInt + DBra(1)*DBra(2)*DBra(3)*DKet(1)*DKet(2)*DKet(3)*p3eInt - enddo - enddo - enddo - enddo - enddo - enddo + end do + end do + end do + end do + end do + end do call cpu_time(end_c3eInt) nc3eInt = nc3eInt + 1 @@ -271,50 +271,50 @@ subroutine Compute3eInt(debug,iType,nShell, & if(.true.) then write(*,'(A15,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6,1X,I6,1X,F16.10)') & '(a1a2a3|b1b2b3) = ',iBasA1,iBasA2,iBasA3,iBasB1,iBasB2,iBasB3,c3eInt - endif - endif + end if + end if !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB3) - enddo + end do iBasB3 = 0 !------------------------------------------------------------------------ ! End loops over shell B3 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB2) - enddo + end do iBasB2 = 0 !------------------------------------------------------------------------ ! End loops over shell B2 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB1) - enddo + end do iBasB1 = 0 !------------------------------------------------------------------------ ! End loops over shell B1 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA3) - enddo + end do iBasA3 = 0 !------------------------------------------------------------------------ ! End loops over shell A3 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA2) - enddo + end do iBasA2 = 0 !------------------------------------------------------------------------ ! End loops over shell A2 !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA1) - enddo + end do iBasA1 = 0 !------------------------------------------------------------------------ ! End loops over shell A1 diff --git a/src/IntPak/Compute4eInt.f90 b/src/IntPak/Compute4eInt.f90 index 0f21fef..3127942 100644 --- a/src/IntPak/Compute4eInt.f90 +++ b/src/IntPak/Compute4eInt.f90 @@ -29,7 +29,7 @@ subroutine Compute4eInt(debug,nEl,iType,nShell,ExpS, & ShellFunctionC(:,:),ShellFunctionD(:,:) double precision :: ExpA,ExpB,ExpC,ExpD double precision,allocatable :: DA,DB,DC,DD - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA,iBasB,iBasC,iBasD integer :: iShA,iShB,iShC,iShD @@ -166,16 +166,16 @@ subroutine Compute4eInt(debug,nEl,iType,nShell,ExpS, & do iKA=1,KA ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*NormCoeff(ExpA,AngMomA) + DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) do iKB=1,KB ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*NormCoeff(ExpB,AngMomB) + DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) do iKC=1,KC ExpC = ExpShell(iShC,iKC) - DC = DShell(iShC,iKC)*NormCoeff(ExpC,AngMomC) + DC = DShell(iShC,iKC)*norm_coeff(ExpC,AngMomC) do iKD=1,KD ExpD = ExpShell(iShD,iKD) - DD = DShell(iShD,iKD)*NormCoeff(ExpD,AngMomD) + DD = DShell(iShD,iKD)*norm_coeff(ExpD,AngMomD) ! Erf module ! call ErfInt(debug,npErf,nSigpErf, & @@ -188,10 +188,10 @@ subroutine Compute4eInt(debug,nEl,iType,nShell,ExpS, & ! cErf = cErf + DA*DB*DC*DD*pErf - enddo - enddo - enddo - enddo + end do + end do + end do + end do call cpu_time(end_cErf) ncErf = ncErf + 1 @@ -203,36 +203,36 @@ subroutine Compute4eInt(debug,nEl,iType,nShell,ExpS, & if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6,1X,I6,1X,I6)') & '(ab|erf(r)/r|cd) = ',cErf,iBasA,iBasB,iBasC,iBasD - endif - endif + end if + end if !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionD) - enddo + end do iBasD = 0 !------------------------------------------------------------------------ ! End loops over shell D !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionC) - enddo + end do iBasC = 0 !------------------------------------------------------------------------ ! End loops over shell C !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB) - enddo + end do iBasB = 0 !------------------------------------------------------------------------ ! End loops over shell B !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA) - enddo + end do iBasA = 0 !------------------------------------------------------------------------ ! End loops over shell A diff --git a/src/IntPak/ComputeKin.f90 b/src/IntPak/ComputeKin.f90 index 2a9fd0c..b50c4fb 100644 --- a/src/IntPak/ComputeKin.f90 +++ b/src/IntPak/ComputeKin.f90 @@ -26,7 +26,7 @@ subroutine ComputeKin(debug,nShell, & integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) double precision :: ExpA,ExpB double precision,allocatable :: DA,DB - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA,iBasB integer :: iShA,iShB @@ -115,10 +115,10 @@ subroutine ComputeKin(debug,nShell, & do iKA=1,KA ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*NormCoeff(ExpA,AngMomA) + DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) do iKB=1,KB ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*NormCoeff(ExpB,AngMomB) + DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) call KinInt(npKin,nSigpKin, & ExpA,CenterA,AngMomA, & @@ -127,8 +127,8 @@ subroutine ComputeKin(debug,nShell, & cKin = cKin + DA*DB*pKin - enddo - enddo + end do + end do call cpu_time(end_cKin) ncKin = ncKin + 1 @@ -138,21 +138,21 @@ subroutine ComputeKin(debug,nShell, & write(9,'(I6,I6,F20.15)') iBasA,iBasB,cKin if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|T|b) = ',cKin,iBasA,iBasB - endif - endif + end if + end if !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB) - enddo + end do iBasB = 0 !------------------------------------------------------------------------ ! End loops over shell B !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA) - enddo + end do iBasA = 0 !------------------------------------------------------------------------ ! End loops over shell A diff --git a/src/IntPak/ComputeNuc.f90 b/src/IntPak/ComputeNuc.f90 index 40e10f9..92b649f 100644 --- a/src/IntPak/ComputeNuc.f90 +++ b/src/IntPak/ComputeNuc.f90 @@ -29,7 +29,7 @@ subroutine ComputeNuc(debug,nShell, & integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) double precision :: ExpA,ExpB,ZC double precision,allocatable :: DA,DB - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA,iBasB integer :: iShA,iShB,iNucC @@ -131,10 +131,10 @@ subroutine ComputeNuc(debug,nShell, & do iKA=1,KA ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*NormCoeff(ExpA,AngMomA) + DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) do iKB=1,KB ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*NormCoeff(ExpB,AngMomB) + DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) call NucInt(debug,npNuc,nSigpNuc, & ExpA,CenterA,AngMomA, & @@ -144,12 +144,12 @@ subroutine ComputeNuc(debug,nShell, & cNuc = cNuc - DA*DB*ZC*pNuc - enddo - enddo + end do + end do !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do call cpu_time(end_cNuc) !------------------------------------------------------------------------ ! End loops over nuclear centers C @@ -163,19 +163,19 @@ subroutine ComputeNuc(debug,nShell, & if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|V|b) = ',cNuc,iBasA,iBasB write(*,*) - endif - endif + end if + end if - enddo + end do deallocate(ShellFunctionB) - enddo + end do iBasB = 0 !------------------------------------------------------------------------ ! End loops over shell B !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA) - enddo + end do iBasA = 0 !------------------------------------------------------------------------ ! End loops over shell A diff --git a/src/IntPak/ComputeOv.f90 b/src/IntPak/ComputeOv.f90 index 9c8c5d8..9995899 100644 --- a/src/IntPak/ComputeOv.f90 +++ b/src/IntPak/ComputeOv.f90 @@ -26,7 +26,7 @@ subroutine ComputeOv(debug,nBas,nShell, & integer,allocatable :: ShellFunctionA(:,:),ShellFunctionB(:,:) double precision :: ExpA,ExpB double precision,allocatable :: DA,DB - double precision :: NormCoeff + double precision :: norm_coeff integer :: iBasA,iBasB integer :: iShA,iShB @@ -117,10 +117,10 @@ subroutine ComputeOv(debug,nBas,nShell, & do iKA=1,KA ExpA = ExpShell(iShA,iKA) - DA = DShell(iShA,iKA)*NormCoeff(ExpA,AngMomA) + DA = DShell(iShA,iKA)*norm_coeff(ExpA,AngMomA) do iKB=1,KB ExpB = ExpShell(iShB,iKB) - DB = DShell(iShB,iKB)*NormCoeff(ExpB,AngMomB) + DB = DShell(iShB,iKB)*norm_coeff(ExpB,AngMomB) call OvInt(npOv,nSigpOv, & ExpA,CenterA,AngMomA, & @@ -129,8 +129,8 @@ subroutine ComputeOv(debug,nBas,nShell, & cOv = cOv + DA*DB*pOv - enddo - enddo + end do + end do call cpu_time(end_cOv) ncOv = ncOv + 1 @@ -141,22 +141,22 @@ subroutine ComputeOv(debug,nBas,nShell, & write(8,'(I6,I6,F20.15)') iBasA,iBasB,cOv if(debug) then write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '(a|b) = ',cOv,iBasA,iBasB - endif - endif + end if + end if !------------------------------------------------------------------------ ! End loops over contraction degrees !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionB) - enddo + end do iBasB = 0 !------------------------------------------------------------------------ ! End loops over shell B !------------------------------------------------------------------------ - enddo + end do deallocate(ShellFunctionA) - enddo + end do iBasA = 0 !------------------------------------------------------------------------ ! End loops over shell A diff --git a/src/IntPak/FormVRR3e.f90 b/src/IntPak/FormVRR3e.f90 index 124fbe9..97904b4 100644 --- a/src/IntPak/FormVRR3e.f90 +++ b/src/IntPak/FormVRR3e.f90 @@ -44,7 +44,7 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do i=1,3 ZetaMat(i,i) = ExpZ(i) - enddo + end do ! print*,'Zeta' ! call matout(3,3,ZetaMat) @@ -64,20 +64,20 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do i=1,3 do j=1,i-1 GMat(i,j) = - ExpG(j,i) - enddo + end do do j=i+1,3 GMat(i,j) = - ExpG(i,j) - enddo - enddo + end do + end do do i=1,3 do j=1,i-1 GMat(i,i) = GMat(i,i) + ExpG(j,i) - enddo + end do do j=i+1,3 GMat(i,i) = GMat(i,i) + ExpG(i,j) - enddo - enddo + end do + end do ! print*,'G' ! call matout(3,3,GMat) @@ -89,10 +89,10 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do k=1,3 CenterY(i,j,k) = CenterZ(i,k) - CenterZ(j,k) Y2Mat(i,j) = Y2Mat(i,j) + CenterY(i,j,k)**2 - enddo + end do YMat(i,j) = sqrt(Y2Mat(i,j)) - enddo - enddo + end do + end do ! print*,'Y' ! call matout(3,3,YMat) @@ -106,8 +106,8 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do j=1,3 Delta0Mat(i,j) = ZetaMat(i,j) + GMat(i,j) Delta1Mat(i,j) = Delta0Mat(i,j) + CMat(i,j) - enddo - enddo + end do + end do ! Form the DY and D2Y matrices @@ -117,10 +117,10 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) DYMat(i,j,k) = KappaCross(i,j,k)*YMat(j,k)/ExpZ(i) do l=1,3 D2YMat(i,j,k,l) = 0.5d0*KappaCross(i,k,l)*KappaCross(j,k,l)/(ExpZ(i)*ExpZ(j)) - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Compute the inverse of the Delta0 and Delta1 matrices @@ -130,8 +130,8 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do j=1,3 InvDelta0Mat(i,j) = Delta0Mat(i,j) InvDelta1Mat(i,j) = Delta1Mat(i,j) - enddo - enddo + end do + end do ! call amove(3,3,Delta0Mat,InvDelta0Mat) ! call amove(3,3,Delta1Mat,InvDelta1Mat) @@ -150,10 +150,10 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do l=1,3 D0Mat(i,j) = D0Mat(i,k) + ZetaMat(i,k)*InvDelta0Mat(k,l)*ZetaMat(l,j) D1Mat(i,j) = D1Mat(i,k) + ZetaMat(i,k)*InvDelta1Mat(k,l)*ZetaMat(l,j) - enddo - enddo - enddo - enddo + end do + end do + end do + end do ! Form the derivative matrices @@ -163,8 +163,8 @@ subroutine FormVRR3e(ExpZ,ExpG,CenterZ,DY0,DY1,D2Y0,D2Y1,delta0,delta1,Y0,Y1) do j=1,3 call CalcTrAB(3,D0Mat,D2YMat,D2Y0(i,j)) call CalcTrAB(3,D1Mat,D2YMat,D2Y1(i,j)) - enddo - enddo + end do + end do ! Compute Y0 and Y1 diff --git a/src/IntPak/G2eInt.f90 b/src/IntPak/G2eInt.f90 index bdccbd7..aa6e90d 100644 --- a/src/IntPak/G2eInt.f90 +++ b/src/IntPak/G2eInt.f90 @@ -40,7 +40,7 @@ function G2eInt(debug,iType, & do i=1,2 ExpZi(i) = 1d0/(ExpBra(i) + ExpKet(i)) - enddo + end do NormABSq = 0d0 do j=1,3 @@ -49,20 +49,20 @@ function G2eInt(debug,iType, & CenterAB(i,j) = CenterBra(i,j) - CenterKet(i,j) CenterZA(i,j) = CenterZ(i,j) - CenterBra(i,j) NormABSq(i) = NormABSq(i) + CenterAB(i,j)**2 - enddo - enddo + end do + end do do i=1,2 GAB(i) = (pi*ExpZi(i))**(1.5d0)*exp(-ExpBra(i)*ExpKet(i)*NormABSq(i)*ExpZi(i)) - enddo + end do ! Pre-computed shell-quartet quantities do i=1,2 do j=1,2 ExpY(i,j) = 1d0/(ExpZi(i) + ExpZi(j)) - enddo - enddo + end do + end do do i=1,2 do j=1,2 @@ -70,9 +70,9 @@ function G2eInt(debug,iType, & do k=1,3 CenterY(i,j,k) = CenterZ(i,k) - CenterZ(j,k) NormYSq(i,j) = NormYSq(i,j) + CenterY(i,j,k)**2 - enddo - enddo - enddo + end do + end do + end do ! fG = (ExpZ(1)*ExpZ(2)*ExpG)/(ExpZ(1)*ExpZ(2) + ExpZ(1)*ExpG + ExpZ(2)*ExpG) fG = 1d0/(ExpZi(1) + 1d0/ExpG + ExpZi(2)) @@ -84,7 +84,7 @@ function G2eInt(debug,iType, & TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) maxm = maxm + TotAngMomBra(i) + TotAngMomKet(i) - enddo + end do ! Pre-compute (00|00)^m @@ -97,7 +97,7 @@ function G2eInt(debug,iType, & call CalcOmYuk(maxm,ExpG,ExpY(1,2),fG,NormYSq(1,2),Om) elseif(iType == 4) then call CalcOmErf(maxm,ExpY(1,2),fG,NormYSq(1,2),Om) - endif + end if call cpu_time(finish_Om) @@ -107,9 +107,9 @@ function G2eInt(debug,iType, & write(*,*) '(00|00)^m' do i=0,maxm write(*,*) i,Om(i) - enddo + end do write(*,*) - endif + end if !------------------------------------------------------------------------ ! Launch reccurence relations! @@ -121,10 +121,10 @@ function G2eInt(debug,iType, & a1a2b1b2 = Om(0) else a1a2b1b2 = VRR2e(0,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - endif + end if else a1a2b1b2 = HRR2e(AngMomBra,AngMomKet,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) - endif + end if call cpu_time(finish_RR) diff --git a/src/IntPak/G3eInt.f90 b/src/IntPak/G3eInt.f90 index e91f721..d6926f5 100644 --- a/src/IntPak/G3eInt.f90 +++ b/src/IntPak/G3eInt.f90 @@ -49,7 +49,7 @@ function G3eInt(debug,iType, & do i=1,3 ExpZ(i) = ExpBra(i) + ExpKet(i) - enddo + end do NormABSq = 0d0 do i=1,3 @@ -58,12 +58,12 @@ function G3eInt(debug,iType, & CenterAB(i,j) = CenterBra(i,j) - CenterKet(i,j) CenterZA(i,j) = CenterZ(i,j) - CenterBra(i,j) NormABSq(i) = NormABSq(i) + CenterAB(i,j)**2 - enddo - enddo + end do + end do do i=1,3 GAB(i) = (pi/ExpZ(i))**(1.5d0)*exp(-ExpBra(i)*ExpKet(i)*NormABSq(i)/ExpZ(i)) - enddo + end do ! Pre-computed shell-sextet quantities @@ -76,7 +76,7 @@ function G3eInt(debug,iType, & TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) maxm = maxm + TotAngMomBra(i) + TotAngMomKet(i) - enddo + end do ! Pre-compute (000|000)^m @@ -91,9 +91,9 @@ function G3eInt(debug,iType, & write(*,*) '(000|000)^m' do i=0,maxm write(*,*) i,Om(i) - enddo + end do write(*,*) - endif + end if !------------------------------------------------------------------------ ! Launch reccurence relations! @@ -104,10 +104,10 @@ function G3eInt(debug,iType, & a1a2a3b1b2b3 = Om(0) else a1a2a3b1b2b3 = VRR3e(0,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if else a1a2a3b1b2b3 = HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if call cpu_time(finish_RR) diff --git a/src/IntPak/GF12Int.f90 b/src/IntPak/GF12Int.f90 index 38be824..ac1f555 100644 --- a/src/IntPak/GF12Int.f90 +++ b/src/IntPak/GF12Int.f90 @@ -87,7 +87,7 @@ function GF12Int(ExpG,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,ExpC,CenterC,Ang do i=1,3 CenterRA(i) = CenterP(i) - CenterA(i) + fP*(CenterQ(i) - CenterP(i)) CenterRC(i) = CenterQ(i) - CenterC(i) + fQ*(CenterP(i) - CenterQ(i)) - enddo + end do !------------------------------------------------------------------------ ! Launch reccurence relations! !------------------------------------------------------------------------ diff --git a/src/IntPak/GenerateShell.f90 b/src/IntPak/GenerateShell.f90 index 933d40b..7be3c00 100644 --- a/src/IntPak/GenerateShell.f90 +++ b/src/IntPak/GenerateShell.f90 @@ -23,8 +23,8 @@ subroutine GenerateShell(atot,nShellFunction,ShellFunction) ShellFunction(ia,1) = ax ShellFunction(ia,2) = ay ShellFunction(ia,3) = az - endif - enddo - enddo + end if + end do + end do end subroutine GenerateShell diff --git a/src/IntPak/HRR2e.f90 b/src/IntPak/HRR2e.f90 index bd22a34..1ab9126 100644 --- a/src/IntPak/HRR2e.f90 +++ b/src/IntPak/HRR2e.f90 @@ -31,7 +31,7 @@ recursive function HRR2e(AngMomBra,AngMomKet, & NegAngMomKet(i) = AngMomKet(i,1) < 0 .or. AngMomKet(i,2) < 0 .or. AngMomKet(i,3) < 0 TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - enddo + end do !------------------------------------------------------------------------ ! Termination condition @@ -52,8 +52,8 @@ recursive function HRR2e(AngMomBra,AngMomKet, & do j=1,3 a1p(i,j) = AngMomBra(i,j) b1m(i,j) = AngMomKet(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomKet(1,1) > 0) then @@ -64,7 +64,7 @@ recursive function HRR2e(AngMomBra,AngMomKet, & xyz = 3 else write(*,*) 'xyz = 0 in HRR2e!' - endif + end if ! End loop over cartesian directions a1p(1,xyz) = a1p(1,xyz) + 1 b1m(1,xyz) = b1m(1,xyz) - 1 @@ -78,8 +78,8 @@ recursive function HRR2e(AngMomBra,AngMomKet, & do j=1,3 a2p(i,j) = AngMomBra(i,j) b2m(i,j) = AngMomKet(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomKet(2,1) > 0) then @@ -90,12 +90,12 @@ recursive function HRR2e(AngMomBra,AngMomKet, & xyz = 3 else write(*,*) 'xyz = 0 in HRR2e!' - endif + end if ! End loop over cartesian directions a2p(2,xyz) = a2p(2,xyz) + 1 b2m(2,xyz) = b2m(2,xyz) - 1 a1a2b1b2 = HRR2e(a2p,b2m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) & + CenterAB(2,xyz)*HRR2e(AngMomBra,b2m,maxm,Om,ExpZi,ExpY,CenterAB,CenterZA,CenterY) - endif + end if end function HRR2e diff --git a/src/IntPak/HRR3e.f90 b/src/IntPak/HRR3e.f90 index 0e1fa21..6d597ca 100644 --- a/src/IntPak/HRR3e.f90 +++ b/src/IntPak/HRR3e.f90 @@ -30,7 +30,7 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, NegAngMomKet(i) = AngMomKet(i,1) < 0 .or. AngMomKet(i,2) < 0 .or. AngMomKet(i,3) < 0 TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) TotAngMomKet(i) = AngMomKet(i,1) + AngMomKet(i,2) + AngMomKet(i,3) - enddo + end do !------------------------------------------------------------------------ ! Termination condition @@ -50,8 +50,8 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, do j=1,3 a1p(i,j) = AngMomBra(i,j) b1m(i,j) = AngMomKet(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomKet(1,1) > 0) then @@ -62,7 +62,7 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, xyz = 3 else write(*,*) 'xyz = 0 in HRR3e!' - endif + end if ! End loop over cartesian directions a1p(1,xyz) = a1p(1,xyz) + 1 b1m(1,xyz) = b1m(1,xyz) - 1 @@ -77,8 +77,8 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, do j=1,3 a2p(i,j) = AngMomBra(i,j) b2m(i,j) = AngMomKet(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomKet(2,1) > 0) then @@ -89,7 +89,7 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, xyz = 3 else write(*,*) 'xyz = 0 in HRR3e!' - endif + end if ! End loop over cartesian directions a2p(2,xyz) = a2p(2,xyz) + 1 b2m(2,xyz) = b2m(2,xyz) - 1 @@ -104,8 +104,8 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, do j=1,3 a3p(i,j) = AngMomBra(i,j) b3m(i,j) = AngMomKet(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomKet(3,1) > 0) then @@ -116,13 +116,13 @@ recursive function HRR3e(AngMomBra,AngMomKet,maxm,Om,ExpZ,CenterAB,CenterZA,DY0, xyz = 3 else write(*,*) 'xyz = 0 in HRR3e!' - endif + end if ! End loop over cartesian directions a3p(3,xyz) = a3p(3,xyz) + 1 b3m(3,xyz) = b3m(3,xyz) - 1 a1a2a3b1b2b3 = HRR3e(a3p,b3m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) & + CenterAB(3,xyz)* & HRR3e(AngMomBra,b3m,maxm,Om,ExpZ,CenterAB,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if end function HRR3e diff --git a/src/IntPak/HRRNuc.f90 b/src/IntPak/HRRNuc.f90 index 822a9b5..50039ee 100644 --- a/src/IntPak/HRRNuc.f90 +++ b/src/IntPak/HRRNuc.f90 @@ -48,7 +48,7 @@ recursive function HRRNuc(AngMomA,AngMomB,maxm,Om,ExpPi,CenterAB,CenterPA,Center do i=1,3 ap(i) = AngMomA(i) bm(i) = AngMomB(i) - enddo + end do ! Loop over cartesian directions xyz = 0 if (AngMomB(1) > 0) then @@ -59,13 +59,13 @@ recursive function HRRNuc(AngMomA,AngMomB,maxm,Om,ExpPi,CenterAB,CenterPA,Center xyz = 3 else write(*,*) 'xyz = 0 in HRRNuc!' - endif + end if ! End loop over cartesian directions ap(xyz) = ap(xyz) + 1 bm(xyz) = bm(xyz) - 1 Gab = HRRNuc(ap,bm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & + CenterAB(xyz)*HRRNuc(AngMomA,bm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - endif - endif + end if + end if end function HRRNuc diff --git a/src/IntPak/HRROv.f90 b/src/IntPak/HRROv.f90 index 135140e..cc7311f 100644 --- a/src/IntPak/HRROv.f90 +++ b/src/IntPak/HRROv.f90 @@ -22,7 +22,7 @@ recursive function HRROv(AngMomA,AngMomB,ExpPi,CenterAB,CenterPA) & else Gab = HRROv(AngMomA+1,AngMomB-1,ExpPi,CenterAB,CenterPA) & + CenterAB*HRROv(AngMomA,AngMomB-1,ExpPi,CenterAB,CenterPA) - endif - endif + end if + end if end function HRROv diff --git a/src/IntPak/IntPak.f90 b/src/IntPak/IntPak.f90 index 8bab7eb..25ac6b2 100644 --- a/src/IntPak/IntPak.f90 +++ b/src/IntPak/IntPak.f90 @@ -109,7 +109,7 @@ program IntPak write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' write(*,*) - endif + end if !------------------------------------------------------------------------ ! Compute one-electron kinetic integrals @@ -137,7 +137,7 @@ program IntPak write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' write(*,*) - endif + end if !------------------------------------------------------------------------ ! Compute one-electron nuclear attraction integrals @@ -166,7 +166,7 @@ program IntPak write(*,'(A65,1X,F9.3,A8)') 'Total CPU time = ',t_1eInt(iType),' seconds' write(*,*) - endif + end if !------------------------------------------------------------------------ ! Compute ERIs @@ -201,7 +201,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute F12 two-electron integrals @@ -263,7 +263,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute Yukawa two-electron integrals @@ -298,7 +298,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute long-range Coulomb two-electron integrals @@ -334,7 +334,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute three-electron integrals: Type 1 => chain C12 S23 @@ -372,7 +372,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute three-electron integrals: Type 2 => cyclic C12 S13 S23 @@ -407,7 +407,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute three-electron integrals: Type 3 => chain S13 S23 @@ -442,7 +442,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute four-electron integrals: Type 1 => chain C12 S14 S23 @@ -477,7 +477,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute four-electron integrals: Type 2 => trident C12 S13 S14 @@ -511,7 +511,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! Compute four-electron integrals: Type 3 => chain C12 S13 S34 @@ -546,7 +546,7 @@ program IntPak deallocate(DG,ExpG) - endif + end if !------------------------------------------------------------------------ ! End of IntPak !------------------------------------------------------------------------ diff --git a/src/IntPak/KinInt.f90 b/src/IntPak/KinInt.f90 index e7394e8..c7f892d 100644 --- a/src/IntPak/KinInt.f90 +++ b/src/IntPak/KinInt.f90 @@ -71,6 +71,6 @@ subroutine KinInt(npKin,nSigpKin,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,pKin) npKin = npKin + 1 if(abs(pKin) > 1d-15) then nSigpKin = nSigpKin + 1 - endif + end if end subroutine KinInt diff --git a/src/IntPak/NucInt.f90 b/src/IntPak/NucInt.f90 index 36e678a..6e18276 100644 --- a/src/IntPak/NucInt.f90 +++ b/src/IntPak/NucInt.f90 @@ -57,7 +57,7 @@ subroutine NucInt(debug,npNuc,nSigpNuc, & CenterPC(i) = CenterP(i) - CenterC(i) NormABSq = NormABSq + CenterAB(i)**2 NormPCSq = NormPCSq + CenterPC(i)**2 - enddo + end do G = (pi*ExpPi)**(1.5d0)*exp(-NormABSq/(ExpAi+ExpBi)) @@ -81,9 +81,9 @@ subroutine NucInt(debug,npNuc,nSigpNuc, & write(*,*) '(0|V|0)^m' do i=0,maxm write(*,*) i,Om(i) - enddo + end do write(*,*) - endif + end if !------------------------------------------------------------------------ ! Launch reccurence relations! @@ -105,7 +105,7 @@ subroutine NucInt(debug,npNuc,nSigpNuc, & if(abs(pNuc) > 1d-15) then nSigpNuc = nSigpNuc + 1 ! write(*,'(A10,1X,F16.10,1X,I6,1X,I6)') '[a|V|b] = ',pNuc - endif + end if ! Deallocate arrays diff --git a/src/IntPak/OvInt.f90 b/src/IntPak/OvInt.f90 index 2426440..c8f930c 100644 --- a/src/IntPak/OvInt.f90 +++ b/src/IntPak/OvInt.f90 @@ -69,6 +69,6 @@ subroutine OvInt(npOv,nSigpOv,ExpA,CenterA,AngMomA,ExpB,CenterB,AngMomB,pOv) npOv = npOv + 1 if(abs(pOv) > 1d-15) then nSigpOv = nSigpOv + 1 - endif + end if end subroutine OvInt diff --git a/src/IntPak/S2eInt.f90 b/src/IntPak/S2eInt.f90 index 833677d..e3475d0 100644 --- a/src/IntPak/S2eInt.f90 +++ b/src/IntPak/S2eInt.f90 @@ -46,7 +46,7 @@ subroutine S2eInt(debug,iType,np2eInt,nSigp2eInt, & ExpKet(1),CenterKet(1,1:3),AngMomKet(1,1:3), & ExpBra(2),CenterBra(2,1:3),AngMomBra(2,1:3), & ExpKet(2),CenterKet(2,1:3),AngMomKet(2,1:3)) - enddo + end do else do k=1,KG ExpSG = ExpG(k)*ExpS**2 @@ -55,8 +55,8 @@ subroutine S2eInt(debug,iType,np2eInt,nSigp2eInt, & ExpSG, & ExpBra,CenterBra,AngMomBra, & ExpKet,CenterKet,AngMomKet) - enddo - endif + end do + end if ! Print result @@ -65,6 +65,6 @@ subroutine S2eInt(debug,iType,np2eInt,nSigp2eInt, & if(abs(p2eInt) > 1d-15) then nSigp2eInt = nSigp2eInt + 1 if(.false.) write(*,'(A15,1X,F16.10)') '[a1a2|b1b2] = ',p2eInt - endif + end if end subroutine S2eInt diff --git a/src/IntPak/S3eInt.f90 b/src/IntPak/S3eInt.f90 index 4c45eb3..997862c 100644 --- a/src/IntPak/S3eInt.f90 +++ b/src/IntPak/S3eInt.f90 @@ -71,6 +71,6 @@ subroutine S3eInt(debug,iType,np3eInt,nSigp3eInt, & if(abs(p3eInt) > 1d-15) then nSigp3eInt = nSigp3eInt + 1 if(.false.) write(*,'(A15,1X,F16.10)') '[a1a2a3|b1b2b3] = ',p3eInt - endif + end if end subroutine S3eInt diff --git a/src/IntPak/VRR2e.f90 b/src/IntPak/VRR2e.f90 index 0dde145..1b2189d 100644 --- a/src/IntPak/VRR2e.f90 +++ b/src/IntPak/VRR2e.f90 @@ -31,7 +31,7 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & do i=1,2 NegAngMomBra(i) = AngMomBra(i,1) < 0 .or. AngMomBra(i,2) < 0 .or. AngMomBra(i,3) < 0 TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - enddo + end do fZ(1) = ExpY(1,2)*ExpZi(1) fZ(2) = ExpY(1,2)*ExpZi(2) @@ -55,8 +55,8 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & do j=1,3 a1m(i,j) = AngMomBra(i,j) a1mm(i,j) = AngMomBra(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomBra(1,1) > 0) then @@ -67,7 +67,7 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & xyz = 3 else write(*,*) 'xyz = 0 in VRR2e!' - endif + end if ! End loop over cartesian directions a1m(1,xyz) = a1m(1,xyz) - 1 a1mm(1,xyz) = a1mm(1,xyz) - 2 @@ -82,7 +82,7 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & + 0.5d0*dble(AngMomBra(1,xyz)-1)*ExpZi(1)*( & VRR2e(m,a1mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - fZ(1)*VRR2e(m+1,a1mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY)) - endif + end if !------------------------------------------------------------------------ ! 2nd vertical recurrence relation (5 terms): (a0|c+0)^m !------------------------------------------------------------------------ @@ -92,8 +92,8 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & a2m(i,j) = AngMomBra(i,j) a2mm(i,j) = AngMomBra(i,j) a1m2m(i,j) = AngMomBra(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomBra(2,1) > 0) then @@ -104,7 +104,7 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & xyz = 3 else write(*,*) 'xyz = 0 in VRR2e!' - endif + end if ! End loop over cartesian directions a2m(2,xyz) = a2m(2,xyz) - 1 a2mm(2,xyz) = a2mm(2,xyz) - 2 @@ -121,10 +121,10 @@ recursive function VRR2e(m,AngMomBra,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & + 0.5d0*dble(AngMomBra(2,xyz)-1)*ExpZi(2)*( & VRR2e(m,a2mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) & - fZ(2)*VRR2e(m+1,a2mm,maxm,Om,ExpZi,ExpY,CenterZA,CenterY)) - endif + end if if(AngMomBra(1,xyz) > 0) & a1a2 = a1a2 & + 0.5d0*dble(AngMomBra(1,xyz))*fZ(2)*ExpZi(1)*VRR2e(m+1,a1m2m,maxm,Om,ExpZi,ExpY,CenterZA,CenterY) - endif + end if end function VRR2e diff --git a/src/IntPak/VRR3e.f90 b/src/IntPak/VRR3e.f90 index 3033bff..a866e97 100644 --- a/src/IntPak/VRR3e.f90 +++ b/src/IntPak/VRR3e.f90 @@ -31,7 +31,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & do i=1,3 NegAngMomBra(i) = AngMomBra(i,1) < 0 .or. AngMomBra(i,2) < 0 .or. AngMomBra(i,3) < 0 TotAngMomBra(i) = AngMomBra(i,1) + AngMomBra(i,2) + AngMomBra(i,3) - enddo + end do !------------------------------------------------------------------------ ! Termination condition @@ -51,8 +51,8 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & do j=1,3 a1m(i,j) = AngMomBra(i,j) a1mm(i,j) = AngMomBra(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomBra(1,1) > 0) then @@ -63,7 +63,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & xyz = 3 else write(*,*) 'xyz = 0 in VRR3e!' - endif + end if ! End loop over cartesian directions a1m(1,xyz) = a1m(1,xyz) - 1 a1mm(1,xyz) = a1mm(1,xyz) - 2 @@ -77,7 +77,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - (DY1(1) - DY0(1))*VRR3e(m+1,a1m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & + dble(AngMomBra(1,xyz)-1)*(0.5d0/ExpZ(1) - D2Y0(1,1))*VRR3e(m, a1mm,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - dble(AngMomBra(1,xyz)-1)*(D2Y1(1,1) - D2Y0(1,1))*VRR3e(m+1,a1mm,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if !------------------------------------------------------------------------ ! 2nd vertical recurrence relation (6 terms): (a1a2+0|000)^m !------------------------------------------------------------------------ @@ -87,8 +87,8 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & a2m(i,j) = AngMomBra(i,j) a2mm(i,j) = AngMomBra(i,j) a1m2m(i,j) = AngMomBra(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomBra(2,1) > 0) then @@ -99,7 +99,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & xyz = 3 else write(*,*) 'xyz = 0 in VRR3e!' - endif + end if ! End loop over cartesian directions a2m(2,xyz) = a2m(2,xyz) - 1 a2mm(2,xyz) = a2mm(2,xyz) - 2 @@ -115,7 +115,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - (DY1(2) - DY0(2))*VRR3e(m+1,a2m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & + dble(AngMomBra(2,xyz)-1)*(0.5d0/ExpZ(2) - D2Y0(2,2))*VRR3e(m, a2mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - dble(AngMomBra(2,xyz)-1)*(D2Y1(2,2) - D2Y0(2,2))*VRR3e(m+1,a2mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if if(AngMomBra(1,xyz) > 0) & a1a2a3 = a1a2a3 & + dble(AngMomBra(1,xyz))*(-D2Y0(2,1))*VRR3e(m, a1m2m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & @@ -130,8 +130,8 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & a3mm(i,j) = AngMomBra(i,j) a1m3m(i,j) = AngMomBra(i,j) a2m3m(i,j) = AngMomBra(i,j) - enddo - enddo + end do + end do ! Loop over cartesian directions xyz = 0 if (AngMomBra(3,1) > 0) then @@ -142,7 +142,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & xyz = 3 else write(*,*) 'xyz = 0 in VRR3e!' - endif + end if ! End loop over cartesian directions a3m(3,xyz) = a3m(3,xyz) - 1 a3mm(3,xyz) = a3mm(3,xyz) - 2 @@ -160,7 +160,7 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - (DY1(3) - DY0(3))*VRR3e(m+1,a3m, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & + dble(AngMomBra(3,xyz)-1)*(0.5d0/ExpZ(3) - D2Y0(3,3))*VRR3e(m, a3mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - dble(AngMomBra(3,xyz)-1)*(D2Y1(3,3) - D2Y0(3,3))*VRR3e(m+1,a3mm, maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if if(dble(AngMomBra(1,xyz)) > 0) & a1a2a3 = a1a2a3 & + dble(AngMomBra(1,xyz))*(-D2Y0(3,1))*VRR3e(m, a1m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & @@ -169,6 +169,6 @@ recursive function VRR3e(m,AngMomBra,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & a1a2a3 = a1a2a3 & + dble(AngMomBra(2,xyz))*(-D2Y0(3,2))*VRR3e(m, a2m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) & - dble(AngMomBra(2,xyz))*(D2Y1(3,2) - D2Y0(3,2))*VRR3e(m+1,a2m3m,maxm,Om,ExpZ,CenterZA,DY0,DY1,D2Y0,D2Y1) - endif + end if end function VRR3e diff --git a/src/IntPak/VRRNuc.f90 b/src/IntPak/VRRNuc.f90 index 82bd454..026e7a2 100644 --- a/src/IntPak/VRRNuc.f90 +++ b/src/IntPak/VRRNuc.f90 @@ -51,7 +51,7 @@ recursive function VRRNuc(m,AngMomA,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & do i=1,3 am(i) = AngMomA(i) amm(i) = AngMomA(i) - enddo + end do ! Loop over cartesian directions xyz = 0 if (AngMomA(1) > 0) then @@ -62,7 +62,7 @@ recursive function VRRNuc(m,AngMomA,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & xyz = 3 else write(*,*) 'xyz = 0 in VRRNuc!' - endif + end if ! End loop over cartesian directions am(xyz) = am(xyz) - 1 amm(xyz) = amm(xyz) - 2 @@ -70,7 +70,7 @@ recursive function VRRNuc(m,AngMomA,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & + 0.5d0*dble(am(xyz))*ExpPi*VRRNuc(m,amm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - CenterPC(xyz)*ExpPi*VRRNuc(m+1,am,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) & - 0.5d0*dble(am(xyz))*ExpPi**2*VRRNuc(m+1,amm,maxm,Om,ExpPi,CenterAB,CenterPA,CenterPC) - endif - endif + end if + end if end function VRRNuc diff --git a/src/IntPak/VRROv.f90 b/src/IntPak/VRROv.f90 index 0041318..10627ca 100644 --- a/src/IntPak/VRROv.f90 +++ b/src/IntPak/VRROv.f90 @@ -22,7 +22,7 @@ recursive function VRROv(AngMomA,ExpPi,CenterPA) & Ga = 1d0 else Ga = CenterPA*VRROv(AngMomA-1,ExpPi,CenterPA) + 0.5d0*dble(AngMomA-1)*ExpPi*VRROv(AngMomA-2,ExpPi,CenterPA) - endif - endif + end if + end if end function VRROv diff --git a/src/IntPak/obj/.gitignore b/src/IntPak/obj/.gitignore new file mode 100644 index 0000000..5761abc --- /dev/null +++ b/src/IntPak/obj/.gitignore @@ -0,0 +1 @@ +*.o diff --git a/src/QuAcK/AO_values.f90 b/src/QuAcK/AO_values.f90 index f57124d..0050f26 100644 --- a/src/QuAcK/AO_values.f90 +++ b/src/QuAcK/AO_values.f90 @@ -18,7 +18,7 @@ subroutine AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell integer :: atot,nShellFunction,a(3) integer,allocatable :: ShellFunction(:,:) - double precision :: rASq,xA,yA,zA,NormCoeff,prim + double precision :: rASq,xA,yA,zA,norm_coeff,prim integer :: iSh,iShF,iK,iW,iEl,iBas,ixyz @@ -66,7 +66,7 @@ subroutine AO_values(doDrift,nBas,nShell,nWalk,CenterShell,TotAngMomShell,KShell do iK=1,KShell(iSh) ! Calculate the exponential part - prim = DShell(iSh,iK)*NormCoeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq) + prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq) AO(iW,iEl,iBas) = AO(iW,iEl,iBas) + prim if(doDrift) then diff --git a/src/QuAcK/form_EOM_one_body.f90 b/src/QuAcK/form_EOM_one_body.f90 index 40ed0d0..58814e8 100644 --- a/src/QuAcK/form_EOM_one_body.f90 +++ b/src/QuAcK/form_EOM_one_body.f90 @@ -21,8 +21,8 @@ subroutine form_EOM_one_body(nO,nV,foo,fov,fvv,t1,t2,OOOV,OOVV,OVVV,cFoo,cFov,cF ! Local variables - integer :: i,j,k,l,m,n - integer :: a,b,c,d,e,f + integer :: i,j,k,l + integer :: a,b,c,d double precision,allocatable :: tau(:,:,:,:) ! Output variables @@ -54,20 +54,20 @@ subroutine form_EOM_one_body(nO,nV,foo,fov,fvv,t1,t2,OOOV,OOVV,OVVV,cFoo,cFov,cF do a=1,nV do b=1,nV - do m=1,nO - cFvv(a,b) = cFvv(a,b) - fov(m,b)*t1(m,a) + do i=1,nO + cFvv(a,b) = cFvv(a,b) - fov(i,b)*t1(i,a) end do - do m=1,nO - do f=1,nV - cFvv(a,b) = cFvv(a,b) + t1(m,f)*OVVV(m,a,f,b) + do i=1,nO + do c=1,nV + cFvv(a,b) = cFvv(a,b) + t1(i,c)*OVVV(i,a,c,b) end do end do - do m=1,nO - do n=1,nO - do e=1,nV - cFvv(a,b) = cFvv(a,b) - 0.5d0*tau(m,n,a,e)*OOVV(m,n,b,e) + do i=1,nO + do j=1,nO + do c=1,nV + cFvv(a,b) = cFvv(a,b) - 0.5d0*tau(i,j,a,c)*OOVV(i,j,b,c) end do end do end do @@ -79,23 +79,23 @@ subroutine form_EOM_one_body(nO,nV,foo,fov,fvv,t1,t2,OOOV,OOVV,OVVV,cFoo,cFov,cF cFoo(:,:) = foo(:,:) - do j=1,nO - do i=1,nO + do i=1,nO + do j=1,nO - do e=1,nV - cFoo(j,i) = cFoo(j,i) + t1(i,e)*fov(j,e) + do a=1,nV + cFoo(i,j) = cFoo(i,j) + t1(j,a)*fov(i,a) end do - do m=1,nO - do e=1,nV - cFoo(j,i) = cFoo(j,i) + t1(m,e)*OVVV(j,m,i,e) + do k=1,nO + do a=1,nV + cFoo(i,j) = cFoo(i,j) + t1(k,a)*OVVV(i,k,j,a) end do end do - do m=1,nO - do e=1,nV - do f=1,nV - cFoo(j,i) = cFoo(j,i) + 0.5d0*tau(i,m,e,f)*OOVV(j,m,e,f) + do k=1,nO + do a=1,nV + do b=1,nV + cFoo(i,j) = cFoo(i,j) + 0.5d0*tau(j,k,a,b)*OOVV(i,k,a,b) end do end do end do diff --git a/src/QuAcK/form_EOM_two_body.f90 b/src/QuAcK/form_EOM_two_body.f90 index 6b4c46e..cc29837 100644 --- a/src/QuAcK/form_EOM_two_body.f90 +++ b/src/QuAcK/form_EOM_two_body.f90 @@ -72,19 +72,19 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, cWoooo(:,:,:,:) = OOOO(:,:,:,:) - do k=1,nO - do l=1,nO - do i=1,nO - do j=1,nO + do i=1,nO + do j=1,nO + do k=1,nO + do l=1,nO do e=1,nV - cWoooo(k,l,i,j) = cWoooo(k,l,i,j) + t1(j,e)*OOOV(k,l,i,e) - cWoooo(k,l,i,j) = cWoooo(k,l,i,j) - t1(i,e)*OOOV(k,l,j,e) + cWoooo(i,j,k,l) = cWoooo(i,j,k,l) + t1(j,e)*OOOV(i,j,k,e) + cWoooo(i,j,k,l) = cWoooo(i,j,k,l) - t1(i,e)*OOOV(i,j,l,e) end do - do e=1,nV - do f=1,nV - cWoooo(k,l,i,j) = cWoooo(k,l,i,j) + 0.5d0*tau(i,j,e,f)*OOVV(k,l,e,f) + do a=1,nV + do b=1,nV + cWoooo(i,j,k,l) = cWoooo(i,j,k,l) + 0.5d0*tau(k,l,a,b)*OOVV(i,j,a,b) end do end do @@ -102,14 +102,14 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, do c=1,nV do d=1,nV - do m=1,nV - cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) - t1(m,b)*VOVV(a,m,c,d) - cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) + t1(m,a)*VOVV(b,m,c,d) + do i=1,nO + cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) - t1(i,b)*VOVV(a,i,c,d) + cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) + t1(i,a)*VOVV(b,i,c,d) end do - do m=1,nV - do n=1,nV - cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) + tau(m,n,a,b)*OOVV(m,n,c,d) + do i=1,nO + do j=1,nO + cWvvvv(a,b,c,d) = cWvvvv(a,b,c,d) + tau(i,j,a,b)*OOVV(i,j,c,d) end do end do @@ -127,8 +127,8 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, do b=1,nV do c=1,nV - do m=1,nV - cWvovv(a,i,b,c) = cWvovv(a,i,b,c) + t1(m,a)*OOVV(m,i,b,c) + do j=1,nO + cWvovv(a,i,b,c) = cWvovv(a,i,b,c) + t1(j,a)*OOVV(j,i,b,c) end do end do @@ -145,8 +145,8 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, do k=1,nO do a=1,nO - do e=1,nV - cWooov(i,j,k,a) = cWooov(i,j,k,a) + t1(i,e)*OOVV(j,k,e,a) + do b=1,nV + cWooov(i,j,k,a) = cWooov(i,j,k,a) + t1(i,b)*OOVV(j,k,b,a) end do end do @@ -195,36 +195,38 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, do c=1,nV do i=1,nO - do m=1,nO - do e=1,nV - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t2(i,m,b,e)*VOVV(a,m,c,e) - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t2(i,m,a,e)*VOVV(b,m,c,e) + do j=1,nO + do d=1,nV + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t2(i,j,b,d)*VOVV(a,j,c,d) + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t2(i,j,a,d)*VOVV(b,j,c,d) end do end do - do m=1,nO - do n=1,nO - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + 0.5d0*tau(m,n,a,b)*VOOO(c,i,m,n) + do j=1,nO + do k=1,nO + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + 0.5d0*tau(j,k,a,b)*VOOO(c,i,j,k) end do end do - do m=1,nO - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t2(i,m,a,b)*cFov(m,c) + do j=1,nO + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t2(i,j,a,b)*cFov(j,c) end do - do e=1,nV - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(i,e)*cWvvvv(a,b,c,e) + do d=1,nV + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(i,d)*cWvvvv(a,b,c,d) end do - do m=1,nO - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t1(m,a)*OVVO(m,b,c,i) - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(m,b)*OVVO(m,a,c,i) + do j=1,nO + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t1(j,a)*OVVO(j,b,c,i) + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(j,b)*OVVO(j,a,c,i) + do n=1,nO - do e=1,nV - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(m,a)*t2(n,i,b,e)*OOVV(m,n,c,e) - cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t1(m,b)*t2(n,i,a,e)*OOVV(m,n,c,e) + do d=1,nV + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) + t1(j,a)*t2(k,i,b,d)*OOVV(j,k,c,d) + cWvvvo(a,b,c,i) = cWvvvo(a,b,c,i) - t1(j,b)*t2(k,i,a,d)*OOVV(j,k,c,d) end do end do + end do end do @@ -241,36 +243,38 @@ subroutine form_EOM_two_body(nO,nV,foo,fov,fvv,t1,t2,cFoo,cFov,cFvv, do j=1,nO do k=1,nO - do m=1,nO - do e=1,nV - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t2(k,m,a,e)*OOOV(i,m,j,e) - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t2(k,m,a,e)*OOOV(j,m,i,e) + do l=1,nO + do b=1,nV + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t2(k,l,a,b)*OOOV(i,l,j,b) + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t2(k,l,a,b)*OOOV(j,l,i,b) end do end do - do e=1,nV - do f=1,nV - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + 0.5d0*tau(j,k,e,f)*OVVV(i,a,e,f) + do b=1,nV + do c=1,nV + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + 0.5d0*tau(j,k,b,c)*OVVV(i,a,b,c) end do end do - do e=1,nV - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t2(j,k,a,e)*cFov(i,e) + do b=1,nV + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t2(j,k,a,b)*cFov(i,b) end do - do m=1,nO - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(m,a)*cWoooo(i,m,j,k) + do l=1,nO + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(l,a)*cWoooo(i,l,j,k) end do - do e=1,nV - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t1(j,e)*OVVO(i,a,e,k) - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(i,e)*OVVO(j,a,e,k) - do m=1,nO - do f=1,nV - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(j,e)*t2(m,k,a,f)*OOVV(i,m,e,f) - cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t1(i,e)*t2(m,k,a,f)*OOVV(j,m,e,f) + do b=1,nV + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t1(j,b)*OVVO(i,a,b,k) + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(i,b)*OVVO(j,a,b,k) + + do l=1,nO + do c=1,nV + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) + t1(j,b)*t2(l,k,a,c)*OOVV(i,l,b,c) + cWovoo(i,a,j,k) = cWovoo(i,a,j,k) - t1(i,b)*t2(l,k,a,c)*OOVV(j,l,b,c) end do end do + end do end do diff --git a/src/QuAcK/obj/.gitingore b/src/QuAcK/obj/.gitingore new file mode 100644 index 0000000..5761abc --- /dev/null +++ b/src/QuAcK/obj/.gitingore @@ -0,0 +1 @@ +*.o diff --git a/src/eDFT/AO_values_grid.f90 b/src/eDFT/AO_values_grid.f90 index 75416e7..a5a72a6 100644 --- a/src/eDFT/AO_values_grid.f90 +++ b/src/eDFT/AO_values_grid.f90 @@ -21,7 +21,7 @@ subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,E integer :: atot,nShellFunction,a(3) integer,allocatable :: ShellFunction(:,:) - double precision :: rASq,xA,yA,zA,NormCoeff,prim + double precision :: rASq,xA,yA,zA,norm_coeff,prim integer :: iSh,iShF,iK,iG,iBas @@ -68,7 +68,7 @@ subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,E ! Calculate the exponential part - prim = DShell(iSh,iK)*NormCoeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq) + prim = DShell(iSh,iK)*norm_coeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq) AO(iBas,iG) = AO(iBas,iG) + prim prim = -2d0*ExpShell(iSh,iK)*prim diff --git a/src/eDFT/obj/.gitingore b/src/eDFT/obj/.gitingore new file mode 100644 index 0000000..5761abc --- /dev/null +++ b/src/eDFT/obj/.gitingore @@ -0,0 +1 @@ +*.o diff --git a/src/utils/elements.f90 b/src/utils/elements.f90 new file mode 100644 index 0000000..22953dc --- /dev/null +++ b/src/utils/elements.f90 @@ -0,0 +1,170 @@ +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/utils/norm_coeff.f90 b/src/utils/norm_coeff.f90 new file mode 100644 index 0000000..b28260f --- /dev/null +++ b/src/utils/norm_coeff.f90 @@ -0,0 +1,29 @@ +function norm_coeff(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 norm_coeff + + 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))) + + + norm_coeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot + norm_coeff = norm_coeff/(dfa(1)*dfa(2)*dfa(3)) + norm_coeff = sqrt(norm_coeff) + +end function norm_coeff diff --git a/src/utils/obj/.gitignore b/src/utils/obj/.gitignore new file mode 100644 index 0000000..5761abc --- /dev/null +++ b/src/utils/obj/.gitignore @@ -0,0 +1 @@ +*.o diff --git a/src/utils/read_auxiliary_basis.f90 b/src/utils/read_auxiliary_basis.f90 new file mode 100644 index 0000000..cf4f6ed --- /dev/null +++ b/src/utils/read_auxiliary_basis.f90 @@ -0,0 +1,176 @@ +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/utils/read_basis.f90 b/src/utils/read_basis.f90 new file mode 100644 index 0000000..075e1d6 --- /dev/null +++ b/src/utils/read_basis.f90 @@ -0,0 +1,118 @@ +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/utils/read_geometry.f90 b/src/utils/read_geometry.f90 new file mode 100644 index 0000000..a3a4cf9 --- /dev/null +++ b/src/utils/read_geometry.f90 @@ -0,0 +1,68 @@ +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/utils/read_integrals.f90 b/src/utils/read_integrals.f90 new file mode 100644 index 0000000..e6f9dd6 --- /dev/null +++ b/src/utils/read_integrals.f90 @@ -0,0 +1,119 @@ +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/utils/read_molecule.f90 b/src/utils/read_molecule.f90 new file mode 100644 index 0000000..5d8886e --- /dev/null +++ b/src/utils/read_molecule.f90 @@ -0,0 +1,64 @@ +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/utils/utils.f90 b/src/utils/utils.f90 new file mode 100644 index 0000000..3ff50d2 --- /dev/null +++ b/src/utils/utils.f90 @@ -0,0 +1,533 @@ +!------------------------------------------------------------------------ +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 + +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(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 + +!------------------------------------------------------------------------ + +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 diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 new file mode 100644 index 0000000..6c29ab7 --- /dev/null +++ b/src/utils/wrap_lapack.f90 @@ -0,0 +1,207 @@ +!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 +