diff --git a/examples/RKS_emuc_Goli2022/MuCN-euc.lowdin b/examples/RKS_emuc_Goli2022/MuCN-euc.lowdin new file mode 100644 index 00000000..6144f2a3 --- /dev/null +++ b/examples/RKS_emuc_Goli2022/MuCN-euc.lowdin @@ -0,0 +1,122 @@ +!This calculation reproduces the results for a molecule in Goli and Shahbazian muon correlation functional paper +!The paper reports total energy: -93.3882, mu kinetic: 0.0189, C-mu bond length: 1.072 +!From this calculation we get : -93.3908, : 0.0192, : 1.074 +!there is a small difference because we use cartesian and not spherical GTFs. + +SYSTEM_DESCRIPTION='Muon functional test - MuCN from 10.1063/5.0077179' + +!The geometry is taken from their supporting information +!We use the same basis sets +GEOMETRY +e-(C) PC-2 0.0 0.0 0.0 +e-(N) PC-2 0.0 0.0 -1.14598 +e-(H) PC-2 0.0 0.0 0.82596 +u+ 14S14P14D 0.0 0.0 0.82596 +C dirac 0.0 0.0 0.0 +N dirac 0.0 0.0 -1.14598 +END GEOMETRY + +TASKS +method = "RKS" +END TASKS + +!EMUC-1 is the functional defined in their paper +!For MCDFT calculations is better to start from reference vectors from HF, or DFT with point charges calculations +!But for this one, it manages to converge from the HCORE default guess +!To converge the SCF, we have to add a level shifting for the muons (this also applies to nuclei and positrons in MCDFT calculations) +CONTROL +electronExchangeCorrelationFunctional="B3LYP" +muonElectronCorrelationFunctional="EMUC-1" +readCoefficients=.F. +nonElectronicLevelShifting=0.10 +END CONTROL + +!Even-tempered parameters for this custom basis set are defined in the paper +BASIS 14S14P14D +O-U+ U+ (14S14P14D) BASIS TYPE: 2 +42 +0 0 1 +0.707106781187 1.0 +0 0 1 +1.000000000000 1.0 +2 0 1 +1.414213562373 1.0 +3 0 1 +2.000000000000 1.0 +4 0 1 +2.828427124746 1.0 +5 0 1 +4.000000000000 1.0 +6 0 1 +5.656854249492 1.0 +7 0 1 +8.000000000000 1.0 +8 0 1 +11.313708498985 1.0 +9 0 1 +16.000000000000 1.0 +10 0 1 +22.627416997970 1.0 +10 0 1 +32.000000000000 1.0 +12 0 1 +45.254833995939 1.0 +13 0 1 +64.000000000000 1.0 +0 1 1 +0.707106781187 1.0 +1 1 1 +1.000000000000 1.0 +2 1 1 +1.414213562373 1.0 +3 1 1 +2.000000000000 1.0 +4 1 1 +2.828427124746 1.0 +5 1 1 +4.000000000000 1.0 +6 1 1 +5.656854249492 1.0 +7 1 1 +8.000000000000 1.0 +8 1 1 +11.313708498985 1.0 +9 1 1 +16.000000000000 1.0 +10 1 1 +22.627416997970 1.0 +11 1 1 +32.000000000000 1.0 +12 1 1 +45.254833995939 1.0 +13 1 1 +64.000000000000 1.0 +0 2 1 +0.707106781187 1.0 +1 2 1 +1.000000000000 1.0 +2 2 1 +1.414213562373 1.0 +3 2 1 +2.000000000000 1.0 +4 2 1 +2.828427124746 1.0 +5 2 1 +4.000000000000 1.0 +6 2 1 +5.656854249492 1.0 +7 2 1 +8.000000000000 1.0 +8 2 1 +11.313708498985 1.0 +9 2 1 +16.000000000000 1.0 +10 2 1 +22.627416997970 1.0 +12 2 1 +32.000000000000 1.0 +12 2 1 +45.254833995939 1.0 +13 2 1 +64.000000000000 1.0 +END BASIS diff --git a/examples/RKS_emuc_Goli2022/MuCN-noeuc.lowdin b/examples/RKS_emuc_Goli2022/MuCN-noeuc.lowdin new file mode 100644 index 00000000..1d69c6da --- /dev/null +++ b/examples/RKS_emuc_Goli2022/MuCN-noeuc.lowdin @@ -0,0 +1,117 @@ +!This calculation reproduces the results for a molecule in Goli and Shahbazian muon correlation functional paper +!The paper reports total energy: -93.3077, mu kinetic: 0.0432, C-mu bond length: 1.135 +!From this calculation we get : -93.3087, : 0.0427, : 1.137 +!there is a small difference because we use cartesian and not spherical GTFs. + +SYSTEM_DESCRIPTION='No muon functional test - MuCN from 10.1063/5.0077179' + +!The geometry is taken from their supporting information +!We use the same basis sets +GEOMETRY +e-(C) PC-2 0.0 0.0 0.0 +e-(N) PC-2 0.0 0.0 -1.14737 +e-(H) PC-2 0.0 0.0 1.11091 +u+ 14S14P14D 0.0 0.0 1.11091 +C dirac 0.0 0.0 0.0 +N dirac 0.0 0.0 -1.14737 +END GEOMETRY + +TASKS +method = "RKS" +END TASKS + +!For DFT calculations without interspecies correlation functional, the SCF convergence is easier to achieve. +CONTROL +electronExchangeCorrelationFunctional="B3LYP" +readCoefficients=.F. +END CONTROL + +!Even-tempered parameters for this custom basis set are defined in the paper +BASIS 14S14P14D +O-U+ U+ (14S14P14D) BASIS TYPE: 2 +42 +0 0 1 +0.707106781187 1.0 +0 0 1 +1.000000000000 1.0 +2 0 1 +1.414213562373 1.0 +3 0 1 +2.000000000000 1.0 +4 0 1 +2.828427124746 1.0 +5 0 1 +4.000000000000 1.0 +6 0 1 +5.656854249492 1.0 +7 0 1 +8.000000000000 1.0 +8 0 1 +11.313708498985 1.0 +9 0 1 +16.000000000000 1.0 +10 0 1 +22.627416997970 1.0 +10 0 1 +32.000000000000 1.0 +12 0 1 +45.254833995939 1.0 +13 0 1 +64.000000000000 1.0 +0 1 1 +0.707106781187 1.0 +1 1 1 +1.000000000000 1.0 +2 1 1 +1.414213562373 1.0 +3 1 1 +2.000000000000 1.0 +4 1 1 +2.828427124746 1.0 +5 1 1 +4.000000000000 1.0 +6 1 1 +5.656854249492 1.0 +7 1 1 +8.000000000000 1.0 +8 1 1 +11.313708498985 1.0 +9 1 1 +16.000000000000 1.0 +10 1 1 +22.627416997970 1.0 +11 1 1 +32.000000000000 1.0 +12 1 1 +45.254833995939 1.0 +13 1 1 +64.000000000000 1.0 +0 2 1 +0.707106781187 1.0 +1 2 1 +1.000000000000 1.0 +2 2 1 +1.414213562373 1.0 +3 2 1 +2.000000000000 1.0 +4 2 1 +2.828427124746 1.0 +5 2 1 +4.000000000000 1.0 +6 2 1 +5.656854249492 1.0 +7 2 1 +8.000000000000 1.0 +8 2 1 +11.313708498985 1.0 +9 2 1 +16.000000000000 1.0 +10 2 1 +22.627416997970 1.0 +12 2 1 +32.000000000000 1.0 +12 2 1 +45.254833995939 1.0 +13 2 1 +64.000000000000 1.0 +END BASIS diff --git a/src/Core/CONTROL.f90 b/src/Core/CONTROL.f90 index 19a48939..c5206c31 100644 --- a/src/Core/CONTROL.f90 +++ b/src/Core/CONTROL.f90 @@ -267,6 +267,7 @@ module CONTROL_ character(50) :: ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL character(50) :: NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL character(50) :: POSITRON_ELECTRON_CORRELATION_FUNCTIONAL + character(50) :: MUON_ELECTRON_CORRELATION_FUNCTIONAL character(50) :: BETA_FUNCTION integer :: GRID_RADIAL_POINTS integer :: GRID_ANGULAR_POINTS @@ -618,6 +619,7 @@ module CONTROL_ character(50) :: LowdinParameters_electronExchangeCorrelationFunctional character(50) :: LowdinParameters_nuclearElectronCorrelationFunctional character(50) :: LowdinParameters_positronElectronCorrelationFunctional + character(50) :: LowdinParameters_muonElectronCorrelationFunctional character(50) :: LowdinParameters_betaFunction integer :: LowdinParameters_gridRadialPoints integer :: LowdinParameters_gridAngularPoints @@ -951,6 +953,7 @@ module CONTROL_ LowdinParameters_electronExchangeCorrelationFunctional, & LowdinParameters_nuclearElectronCorrelationFunctional, & LowdinParameters_positronElectronCorrelationFunctional, & + LowdinParameters_muonElectronCorrelationFunctional, & LowdinParameters_betaFunction, & LowdinParameters_gridRadialPoints, & LowdinParameters_gridAngularPoints, & @@ -1311,6 +1314,7 @@ subroutine CONTROL_start() LowdinParameters_electronExchangeCorrelationFunctional = "NONE" LowdinParameters_nuclearElectronCorrelationFunctional = "NONE" LowdinParameters_positronElectronCorrelationFunctional = "NONE" + LowdinParameters_muonElectronCorrelationFunctional = "NONE" LowdinParameters_betaFunction = "NONE" LowdinParameters_gridRadialPoints = 35 LowdinParameters_gridAngularPoints = 110 @@ -1658,6 +1662,7 @@ subroutine CONTROL_start() CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL = "NONE" CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL = "NONE" CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL = "NONE" + CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL = "NONE" CONTROL_instance%BETA_FUNCTION = "NONE" CONTROL_instance%GRID_RADIAL_POINTS = 35 CONTROL_instance%GRID_ANGULAR_POINTS = 110 @@ -2060,6 +2065,7 @@ subroutine CONTROL_load(unit) CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL = LowdinParameters_electronExchangeCorrelationFunctional CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL = LowdinParameters_nuclearElectronCorrelationFunctional CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL = LowdinParameters_positronElectronCorrelationFunctional + CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL = LowdinParameters_muonElectronCorrelationFunctional CONTROL_instance%BETA_FUNCTION = LowdinParameters_betaFunction CONTROL_instance%GRID_RADIAL_POINTS = LowdinParameters_gridRadialPoints CONTROL_instance%GRID_ANGULAR_POINTS = LowdinParameters_gridAngularPoints @@ -2428,6 +2434,7 @@ subroutine CONTROL_save(unit, lastStep, firstStep) LowdinParameters_electronExchangeCorrelationFunctional = CONTROL_instance%ELECTRON_EXCHANGE_CORRELATION_FUNCTIONAL LowdinParameters_nuclearElectronCorrelationFunctional = CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL LowdinParameters_positronElectronCorrelationFunctional = CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL + LowdinParameters_muonElectronCorrelationFunctional = CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL LowdinParameters_betaFunction = CONTROL_instance%BETA_FUNCTION LowdinParameters_gridRadialPoints = CONTROL_instance%GRID_RADIAL_POINTS LowdinParameters_gridAngularPoints = CONTROL_instance%GRID_ANGULAR_POINTS @@ -2571,6 +2578,7 @@ subroutine CONTROL_show() write (*, "(T10,A)") "ELECTRON-NUCLEAR CORRELATION FUNCTIONAL: "//trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) write (*, "(T10,A)") "ELECTRON-POSITRON CORRELATION FUNCTIONAL: "//trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + write (*, "(T10,A)") "ELECTRON-MUON CORRELATION FUNCTIONAL: "//trim(CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL) write (*, "(T10,A,I5,A,I5)") "SCF ATOMIC RADIALxANGULAR GRID SIZE:", CONTROL_instance%GRID_RADIAL_POINTS, "x", CONTROL_instance%GRID_ANGULAR_POINTS if (CONTROL_instance%FINAL_GRID_ANGULAR_POINTS*CONTROL_instance%FINAL_GRID_RADIAL_POINTS .gt. & CONTROL_instance%GRID_ANGULAR_POINTS*CONTROL_instance%GRID_RADIAL_POINTS) then diff --git a/src/Core/MolecularSystem.f90 b/src/Core/MolecularSystem.f90 index 0e99eb0b..eeaa0b41 100644 --- a/src/Core/MolecularSystem.f90 +++ b/src/Core/MolecularSystem.f90 @@ -200,6 +200,10 @@ subroutine MolecularSystem_build() MolecularSystem_instance%species(i)%ocupationNumber, "please check your input addParticles and multiplicity" call MolecularSystem_exception(ERROR, "Fractional ocupation number, imposible combination of charge and multiplicity", "MolecularSystem module at build function.") end if + + !!Check if a harmonic potential was loaded + if (MolecularSystem_getOmega(i, MolecularSystem_instance) .ne. 0.0_8) CONTROL_instance%ARE_THERE_QDO_POTENTIALS = .true. + end do do i = 1, MolecularSystem_instance%numberOfPointCharges @@ -1510,7 +1514,9 @@ function MolecularSystem_getPointChargesEnergy(this) result(output) output = 0.0_8 do i = 1, size(system%pointCharges) + if(system%pointCharges(i)%charge .eq. 0.0_8) cycle do j = i + 1, size(system%pointCharges) + if(system%pointCharges(j)%charge .eq. 0.0_8) cycle deltaOrigin = system%pointCharges(i)%origin & - system%pointCharges(j)%origin diff --git a/src/DFT/Functional.f90 b/src/DFT/Functional.f90 index 86483bbc..05d9d964 100644 --- a/src/DFT/Functional.f90 +++ b/src/DFT/Functional.f90 @@ -50,6 +50,7 @@ module Functional_ Functional_libxcEvaluate, & Functional_LDAEvaluate, & Functional_EPCEvaluate, & + Functional_EMUCEvaluate, & Functional_IKNEvaluate, & Functional_MLCSEvaluate, & Functional_MLCSAEvaluate, & @@ -286,6 +287,8 @@ subroutine Functional_constructor(this, speciesID, otherSpeciesID, molSys) auxstring = trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) else if (trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then auxstring = trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + else if (trim(CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring = trim(CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL) else auxstring = "NONE" end if @@ -353,15 +356,15 @@ subroutine Functional_show(these) if (this%correlationName .ne. "NONE") then - write (*, "(T5,A10,A10,A5,A12,A)") trim(this%symbol1), trim(this%symbol2), "", "exchange:", xc_f03_func_info_get_name(this%info1) + write (*, "(T5,A10,A10,A5,A12,A)") trim(this%symbol1), trim(this%symbol2), "", "exchange:", trim(xc_f03_func_info_get_name(this%info1)) ! print *, "family", xc_f03_func_info_get_family(this%info1), "shell", this%shell - write (*, "(T5,A10,A10,A5,A12,A)") trim(this%symbol1), trim(this%symbol2), "", "correlation:", xc_f03_func_info_get_name(this%info2) + write (*, "(T5,A10,A10,A5,A12,A)") trim(this%symbol1), trim(this%symbol2), "", "correlation:", trim(xc_f03_func_info_get_name(this%info2)) ! print *, "family", xc_f03_func_info_get_family(this%info2), "shell", this%shell else - write (*, "(T5,A10,A10,A5,A21,A)") trim(this%symbol1), trim(this%symbol2), "", "exchange-correlation:", xc_f03_func_info_get_name(this%info1) + write (*, "(T5,A10,A10,A5,A21,A)") trim(this%symbol1), trim(this%symbol2), "", "exchange-correlation:", trim(xc_f03_func_info_get_name(this%info1)) ! print *, "family", xc_f03_func_info_get_family(this%info1), "shell", this%shell @@ -421,8 +424,8 @@ subroutine Functional_show(these) else if (CONTROL_instance%BETA_FUNCTION .eq. "NEWNEWBETA") then if (this%mass2 .gt. 2.0) then !hydrogen - STOP "this beta function only works for electron-positron" - else !positron + call Exception_stopError("this beta function only works for electron-positron","At Functional_show") + else !positron a0 = 2.2919886876120283056 Eab = 0.25 Eab2 = 0.2620050702329801 @@ -572,7 +575,7 @@ subroutine Functional_EPCEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) c = 3.2 else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_EPCEvaluate") end if ec(1:n) = 0.0 @@ -598,6 +601,56 @@ subroutine Functional_EPCEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end subroutine Functional_EPCEvaluate + subroutine Functional_EMUCEvaluate(this, n, rhoABsum, rhoABdiff, rhoMu, ec, vcA, vcB, vcMu) + ! Evaluates Goli and Shabazian electron-positive muon correlation functional + ! Alpha and Beta electronic densities + ! Felix Moncada, 2026 + implicit none + type(Functional):: this !!type of functional + integer(8) :: n !!nuclear gridSize + real(8) :: rhoABsum(*), rhoABdiff(*), rhoMu(*) !! electron and muon Densities - input + real(8) :: ec(*) !! Energy density - output + real(8) :: vcA(*), vcB(*), vcMu(*) !! Potentials - output + + real(8) :: rhoA, rhoB, denominatorA, denominatorB, densityThreshold + integer :: i + + densityThreshold = CONTROL_instance%NUCLEAR_ELECTRON_DENSITY_THRESHOLD !TODO: add to other functionals + + if (this%name .ne. "correlation:EMUC-1") then + print *, this%name + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_EMUCEvaluate") + end if + + ec(1:n) = 0.0 + vcA(1:n) = 0.0 + vcB(1:n) = 0.0 + vcMu(1:n) = 0.0 + + ! print *, "i, rhoE, rhoN, denominator, energy density, potentialE, potentialN" + !$omp parallel private(rhoA,rhoB,denominatorA,denominatorB) + !$omp do schedule (dynamic) + do i = 1, n + if (rhoABsum(i) + rhoMu(i) .lt. densityThreshold) cycle + rhoA=(rhoABsum(i)+rhoABdiff(i))/2.0 + rhoB=(rhoABsum(i)-rhoABdiff(i))/2.0 + denominatorA = 1.0 + 4.0*rhoA*rhoMu(i)*(2.0+sqrt(rhoMu(i))) + denominatorB = 1.0 + 4.0*rhoB*rhoMu(i)*(2.0+sqrt(rhoMu(i))) +!!!Energy density per electron + ec(i) = (rhoA*rhoMu(i)*(-2.0+sqrt(rhoMu(i)))/denominatorA+& + rhoB*rhoMu(i)*(-2.0+sqrt(rhoMu(i)))/denominatorB)/rhoABsum(i) +!!!Potential + vcA(i) = rhoMu(i)*(-2.0+sqrt(rhoMu(i)))/denominatorA**2.0 + vcB(i) = rhoMu(i)*(-2.0+sqrt(rhoMu(i)))/denominatorB**2.0 + vcMu(i) = rhoA*(-4.0+sqrt(rhoMu(i))*(3.0+16.0*rhoA*rhoMu(i)))/2.0/denominatorA**2.0+& + rhoB*(-4.0+sqrt(rhoMu(i))*(3.0+16.0*rhoB*rhoMu(i)))/2.0/denominatorB**2.0 + ! write(*,"(I0.1,5ES16.6)") i, rhoE(i), rhoN(i), ec(i), vcE(i), vcN(i) + end do + !$omp end do + !$omp end parallel + + end subroutine Functional_EMUCEvaluate + subroutine Functional_IKNEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) ! Evaluates YUTAKA IMAMURA, HIROYOSHI KIRYU, HIROMI NAKAI Colle Salvetti nuclear electron correlation functional J Comput Chem 29: 735–740, 2008 ! Only works for Hydrogen - can be extended @@ -628,7 +681,7 @@ subroutine Functional_IKNEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_IKNEvaluate") end if do i = 1, n @@ -663,8 +716,6 @@ subroutine Functional_IKNEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if end do - ! STOP - end subroutine Functional_IKNEvaluate subroutine Functional_MLCSEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) @@ -699,7 +750,7 @@ subroutine Functional_MLCSEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_MLCSEvaluate") end if do i = 1, n @@ -734,8 +785,6 @@ subroutine Functional_MLCSEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if end do - ! STOP - end subroutine Functional_MLCSEvaluate subroutine Functional_MLCSAEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) @@ -770,7 +819,7 @@ subroutine Functional_MLCSAEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_MLCSAEvaluate") end if do i = 1, n @@ -805,8 +854,6 @@ subroutine Functional_MLCSAEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if end do - ! STOP - end subroutine Functional_MLCSAEvaluate subroutine Functional_MLCSANEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) @@ -841,7 +888,7 @@ subroutine Functional_MLCSANEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_MLCSANEvaluate") end if do i = 1, n @@ -878,8 +925,6 @@ subroutine Functional_MLCSANEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) end if end do - ! STOP - end subroutine Functional_MLCSANEvaluate subroutine Functional_myCSEvaluate(this, mass, npoints, rhoE, rhoN, ec, vcE, vcN) @@ -995,7 +1040,7 @@ subroutine Functional_myCSEvaluate(this, mass, npoints, rhoE, rhoN, ec, vcE, vcN bn = 5.580898227664 else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_MyCSEvaluate") end if if (CONTROL_instance%DUMMY_REAL(1) .ne. 0 .and. CONTROL_instance%DUMMY_REAL(2) .ne. 0) then @@ -1155,7 +1200,7 @@ subroutine Functional_expCSEvaluate(this, mass, npoints, rhoE, rhoP, ec, vcE, vc end if else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_expCSEvaluate") end if p = 1.0 @@ -1261,7 +1306,7 @@ subroutine Functional_expCSGGAEvaluate(this, mass, npoints, electronDensity, ele b0 = 1.0 end if - if (mass .gt. 2.0) STOP "the expCSGGA functional only works for positron-electron correlation at the moment" + if (mass .gt. 2.0) call Exception_stopError("the expCSGGA functional only works for positron-electron correlation at the moment","Functional_expCSGGA") p = 1.0 g1 = 1.0 @@ -1767,7 +1812,7 @@ subroutine Functional_PSNEvaluate(this, mass, n, rhoE, rhoP, ec, vcE, vcP) Cc = 5.21152*2.0_8 else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_PSN") end if ! densityThreshold=CONTROL_instance%NUCLEAR_ELECTRON_DENSITY_THRESHOLD @@ -1872,7 +1917,7 @@ subroutine Functional_PSNAPEvaluate(this, mass, n, rhoE, rhoP, ec, vcE, vcP) xcut = 6.0 else print *, this%name - STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_PSNAP") end if densityThreshold = CONTROL_instance%NUCLEAR_ELECTRON_DENSITY_THRESHOLD @@ -1965,7 +2010,7 @@ subroutine Functional_lowLimitEvaluate(this, mass, n, rhoE, rhoN, ec, vcE, vcN) energyDensity = -0.5_8*mass/(mass + 1.0_8) b = -0.5_8 else - ! STOP "The nuclear electron functional chosen is not implemented" + call Exception_stopError("The nuclear electron functional chosen is not implemented","Functional_LowLimit") end if do i = 1, n diff --git a/src/DFT/GridManager.f90 b/src/DFT/GridManager.f90 index c58ce7e9..a2abc45a 100644 --- a/src/DFT/GridManager.f90 +++ b/src/DFT/GridManager.f90 @@ -698,7 +698,7 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid(Grid_instance, Gr type(Vector) :: energyDensity type(Vector) :: sigma type(Vector) :: densityAB, potentialAB, sigmaAB, sigmaPotentialAB - type(Vector) :: electronicDensityAtOtherGrid, electronicGradientAtOtherGrid(3), electronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid(3) + type(Vector) :: electronicDensityAtOtherGrid, electronicGradientAtOtherGrid(3), spinDensityAtOtherGrid, spinGradientAtOtherGrid(3), electronicPotentialAtOtherGrid, otherElectronicPotentialAtOtherGrid, electronicGradientPotentialAtOtherGrid(3) integer :: i, j, dir integer(8) :: k @@ -715,21 +715,26 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid(Grid_instance, Gr call Vector_constructor(energyDensity, otherGridSize, 0.0_8) call Vector_constructor(electronicDensityAtOtherGrid, otherGridSize, 0.0_8) + call Vector_constructor(spinDensityAtOtherGrid, otherGridSize, 0.0_8) call Vector_constructor(electronicPotentialAtOtherGrid, otherGridSize, 0.0_8) + call Vector_constructor(otherElectronicPotentialAtOtherGrid, otherGridSize, 0.0_8) do dir = 1, 3 call Vector_constructor(electronicGradientAtOtherGrid(dir), otherGridSize, 0.0_8) + call Vector_constructor(spinGradientAtOtherGrid(dir), otherGridSize, 0.0_8) call Vector_constructor(electronicGradientPotentialAtOtherGrid(dir), otherGridSize, 0.0_8) end do !!This adds E-BETA density and gradient call GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommonPoints, speciesID, otherSpeciesID, & GridsCommonPoints(speciesID, otherSpeciesID)%totalSize, int(GridsCommonPoints(speciesID, otherSpeciesID)%points%values), & - electronicDensityAtOtherGrid, electronicGradientAtOtherGrid) + electronicDensityAtOtherGrid, electronicGradientAtOtherGrid,spinDensityAtOtherGrid, spinGradientAtOtherGrid) if (trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then auxstring = trim(CONTROL_instance%NUCLEAR_ELECTRON_CORRELATION_FUNCTIONAL) else if (trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then auxstring = trim(CONTROL_instance%POSITRON_ELECTRON_CORRELATION_FUNCTIONAL) + else if (trim(CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL) .ne. "NONE") then + auxstring = trim(CONTROL_instance%MUON_ELECTRON_CORRELATION_FUNCTIONAL) else auxstring = "NONE" end if @@ -760,6 +765,13 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid(Grid_instance, Gr electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & energyDensity%values, electronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values) + case ("EMUC-1") + call Functional_EMUCEvaluate(Functionals(speciesID, otherSpeciesID), otherGridSize, & + electronicDensityAtOtherGrid%values, spinDensityAtOtherGrid%values, & + Grid_instance(otherSpeciesID)%density%values, & + energyDensity%values, electronicPotentialAtOtherGrid%values, & + otherElectronicPotentialAtOtherGrid%values, Grid_instance(otherSpeciesID)%potential%values) + case ("IKN-NSF") call Functional_IKNEvaluate(Functionals(speciesID, otherSpeciesID), MolecularSystem_getMass(otherSpeciesID, Grid_instance(otherSpeciesID)%molSys), otherGridSize, & electronicDensityAtOtherGrid%values, Grid_instance(otherSpeciesID)%density%values, & @@ -824,7 +836,7 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid(Grid_instance, Gr case default print *, trim(auxstring) - call Exception_stopError("The "//otherNameOfSpecies//"electron functional chosen is not implemented", "at GridManager_getInterspeciesEnergyAndPotentialAtGrid") + call Exception_stopError("The "//trim(otherNameOfSpecies)//"electron functional chosen is not implemented", "at GridManager_getInterspeciesEnergyAndPotentialAtGrid") end select @@ -850,12 +862,17 @@ subroutine GridManager_getInterspeciesEnergyAndPotentialAtGrid(Grid_instance, Gr otherElectronExchangeCorrelationEnergy = otherElectronExchangeCorrelationEnergy + & energyDensity%values(j)*Grid_instance(otherElectronID)%density%values(i)*Grid_instance(otherSpeciesID)%points%values(j, 4) - Grid_instance(otherElectronID)%potential%values(i) = Grid_instance(otherElectronID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) - do dir = 1, 3 - Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) = Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) & + if (trim(auxstring) .eq. "EMUC-1") then + Grid_instance(otherElectronID)%potential%values(i) = Grid_instance(otherElectronID)%potential%values(i) + otherElectronicPotentialAtOtherGrid%values(j) + + else + Grid_instance(otherElectronID)%potential%values(i) = Grid_instance(otherElectronID)%potential%values(i) + electronicPotentialAtOtherGrid%values(j) + do dir = 1, 3 + Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) = Grid_instance(otherElectronID)%gradientPotential(dir)%values(i) & + electronicGradientPotentialAtOtherGrid(dir)%values(j) - end do + end do + end if end if end do @@ -974,7 +991,7 @@ subroutine GridManager_buildExchangeCorrelationMatrix(Grid_instance, GridsCommon end subroutine GridManager_buildExchangeCorrelationMatrix - subroutine GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommonPoints, electronicID, otherSpeciesID, commonGridSize, commonPoints, electronicDensityAtOtherGrid, electronicGradientAtOtherGrid) + subroutine GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommonPoints, electronicID, otherSpeciesID, commonGridSize, commonPoints, electronicDensityAtOtherGrid, electronicGradientAtOtherGrid, spinDensityAtOtherGrid, spinGradientAtOtherGrid) implicit none type(Grid) :: Grid_instance(:) type(Grid) :: GridsCommonPoints(:, :) @@ -983,6 +1000,8 @@ subroutine GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommo integer :: commonPoints(commonGridSize, 2) type(Vector) :: electronicDensityAtOtherGrid type(Vector) :: electronicGradientAtOtherGrid(3) + type(Vector) :: spinDensityAtOtherGrid + type(Vector) :: spinGradientAtOtherGrid(3) character(50) :: nameOfElectron integer :: otherElectronicID @@ -999,7 +1018,9 @@ subroutine GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommo if (nameOfElectron .eq. "E-ALPHA") otherElectronicID = MolecularSystem_getSpeciesID("E-BETA", Grid_instance(electronicID)%molSys) if (nameOfElectron .eq. "E-BETA") otherElectronicID = MolecularSystem_getSpeciesID("E-ALPHA", Grid_instance(electronicID)%molSys) + !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!Why is this not zero? FM 2026 call Vector_constructor(electronicDensityAtOtherGrid, otherGridSize, 1.0E-12_8) + call Vector_constructor(spinDensityAtOtherGrid, otherGridSize, 1.0E-12_8) time1 = omp_get_wtime() @@ -1013,6 +1034,10 @@ subroutine GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommo electronicGradientAtOtherGrid(1)%values(j) = Grid_instance(electronicID)%densityGradient(1)%values(i) + Grid_instance(otherElectronicID)%densityGradient(1)%values(i) electronicGradientAtOtherGrid(2)%values(j) = Grid_instance(electronicID)%densityGradient(2)%values(i) + Grid_instance(otherElectronicID)%densityGradient(2)%values(i) electronicGradientAtOtherGrid(3)%values(j) = Grid_instance(electronicID)%densityGradient(3)%values(i) + Grid_instance(otherElectronicID)%densityGradient(3)%values(i) + spinDensityAtOtherGrid%values(j) = Grid_instance(electronicID)%density%values(i) - Grid_instance(otherElectronicID)%density%values(i) + spinGradientAtOtherGrid(1)%values(j) = Grid_instance(electronicID)%densityGradient(1)%values(i) - Grid_instance(otherElectronicID)%densityGradient(1)%values(i) + spinGradientAtOtherGrid(2)%values(j) = Grid_instance(electronicID)%densityGradient(2)%values(i) - Grid_instance(otherElectronicID)%densityGradient(2)%values(i) + spinGradientAtOtherGrid(3)%values(j) = Grid_instance(electronicID)%densityGradient(3)%values(i) - Grid_instance(otherElectronicID)%densityGradient(3)%values(i) else electronicDensityAtOtherGrid%values(j) = Grid_instance(electronicID)%density%values(i) electronicGradientAtOtherGrid(1)%values(j) = Grid_instance(electronicID)%densityGradient(1)%values(i) @@ -1095,21 +1120,23 @@ subroutine GridManager_getContactDensity(Grid_instance, GridsCommonPoints, speci real(8) :: rhoE, rhoP, rhoTot, rhoDif, npos, densityThreshold real(8) :: beta, dBdE, dBdP, d2BdE2, d2BdP2, d2BdEP - integer :: n, nproc + integer :: n, nproc, dir real(8) :: contactDensity, overlapDensity character(50) :: auxstring - type(Vector) :: electronicDensityAtOtherGrid, electronicGradientAtOtherGrid(3), gfactor + type(Vector) :: electronicDensityAtOtherGrid, electronicGradientAtOtherGrid(3), spinDensityAtOtherGrid, spinGradientAtOtherGrid(3), gfactor gridSize = Grid_instance(otherSpeciesID)%totalSize call Vector_constructor(electronicDensityAtOtherGrid, gridSize, 0.0_8) - call Vector_constructor(electronicGradientAtOtherGrid(1), gridSize, 0.0_8) - call Vector_constructor(electronicGradientAtOtherGrid(2), gridSize, 0.0_8) - call Vector_constructor(electronicGradientAtOtherGrid(3), gridSize, 0.0_8) + call Vector_constructor(spinDensityAtOtherGrid, gridSize, 0.0_8) + do dir = 1, 3 + call Vector_constructor(electronicGradientAtOtherGrid(dir), gridSize, 0.0_8) + call Vector_constructor(spinGradientAtOtherGrid(dir), gridSize, 0.0_8) + end do !electrons go on the first position call GridManager_getElectronicDensityInOtherGrid(Grid_instance, GridsCommonPoints, speciesID, otherSpeciesID, & - GridsCommonPoints(speciesID, otherSpeciesID)%totalSize, int(GridsCommonPoints(speciesID, otherSpeciesID)%points%values), electronicDensityAtOtherGrid, electronicGradientAtOtherGrid) + GridsCommonPoints(speciesID, otherSpeciesID)%totalSize, int(GridsCommonPoints(speciesID, otherSpeciesID)%points%values), electronicDensityAtOtherGrid, electronicGradientAtOtherGrid,spinDensityAtOtherGrid, spinGradientAtOtherGrid) call Vector_constructor(gfactor, gridSize, 0.0_8) gfactor%values = 0.0_8 diff --git a/src/SCF/DensityMatrixSCFGuess.f90 b/src/SCF/DensityMatrixSCFGuess.f90 index 548d2721..e53624ce 100644 --- a/src/SCF/DensityMatrixSCFGuess.f90 +++ b/src/SCF/DensityMatrixSCFGuess.f90 @@ -38,10 +38,11 @@ module DensityMatrixSCFGuess_ !> !! @brief Obtiene la matriz de densidad inicial - subroutine DensityMatrixSCFGuess_getGuess(speciesID, hcoreMatrix, transformationMatrix, densityMatrix, orbitals, printInfo, system) + subroutine DensityMatrixSCFGuess_getGuess(speciesID, hcoreMatrix, overlapMatrix, transformationMatrix, densityMatrix, orbitals, printInfo, system) implicit none integer, intent(in) :: speciesID type(Matrix), intent(in) :: hcoreMatrix + type(Matrix), intent(in) :: overlapMatrix type(Matrix), intent(in) :: transformationMatrix type(Matrix), intent(inout) :: densityMatrix type(Matrix), intent(inout) :: orbitals @@ -50,7 +51,7 @@ subroutine DensityMatrixSCFGuess_getGuess(speciesID, hcoreMatrix, transformation type(MolecularSystem), pointer :: molSys - type(Matrix) :: auxMatrix + type(Matrix) :: guessMatrix, auxMatrix character(30) :: nameOfSpecies, symbolOfSpecies integer(8) :: orderOfMatrix, occupationNumber logical :: existPlain, existBinnary, readSuccess @@ -58,7 +59,8 @@ subroutine DensityMatrixSCFGuess_getGuess(speciesID, hcoreMatrix, transformation character(50) :: wfnFile character(50) :: arguments(20) integer :: wfnUnit - integer :: i, j, k + integer :: i, j, k, ktrial, kprev + real(8) :: normCheck if (present(system)) then molSys => system @@ -144,25 +146,44 @@ subroutine DensityMatrixSCFGuess_getGuess(speciesID, hcoreMatrix, transformation call Matrix_show(orbitals) end if - call Matrix_copyConstructor(auxMatrix, orbitals) + call Matrix_copyConstructor(guessMatrix, orbitals) !! Segment for fractional occupations: introduce fractional occupation if (trim(symbolOfSpecies) == trim(CONTROL_instance%IONIZE_SPECIES(1))) then do i = 1, size(CONTROL_instance%IONIZE_MO) if (CONTROL_instance%IONIZE_MO(i) .gt. 0 .and. CONTROL_instance%MO_FRACTION_OCCUPATION(i) .lt. 1.0_8) then if (printInfo) write (*, "(A,F6.2,A,I5,A,A)") "Removing ", (1.0 - CONTROL_instance%MO_FRACTION_OCCUPATION(i))*100, & " % of the density associated with orbital No. ", CONTROL_instance%IONIZE_MO(i), " of ", trim(symbolOfSpecies) - auxMatrix%values(:, CONTROL_instance%IONIZE_MO(i)) = auxMatrix%values(:, CONTROL_instance%IONIZE_MO(i))*sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION(i)) + guessMatrix%values(:, CONTROL_instance%IONIZE_MO(i)) = guessMatrix%values(:, CONTROL_instance%IONIZE_MO(i))*sqrt(CONTROL_instance%MO_FRACTION_OCCUPATION(i)) end if end do end if call Matrix_constructor(densityMatrix, int(orderOfMatrix, 8), int(orderOfMatrix, 8), 0.0_8) - do i = 1, orderOfMatrix - do j = 1, orderOfMatrix - do k = 1, occupationNumber - densityMatrix%values(i, j) = densityMatrix%values(i, j) + auxMatrix%values(i, k)*auxMatrix%values(j, k) - end do - end do + call Matrix_constructor(auxMatrix, int(orderOfMatrix, 8), int(orderOfMatrix, 8), 0.0_8) + + kprev=0 + do k = 1, occupationNumber + !Check if the orbital that's going to be added is not full of zeros! + !if it fails, try the next orbital + do ktrial = kprev+1, orderOfMatrix + do i = 1, orderOfMatrix + do j = 1, orderOfMatrix + auxMatrix%values(i, j) = guessMatrix%values(i, ktrial)*guessMatrix%values(j, ktrial) + end do + end do + normCheck = sum(transpose(auxMatrix%values)*overlapMatrix%values) + ! print *, "kt, normCheck", ktrial, normCheck + if (abs(normCheck) .gt. CONTROL_instance%DOUBLE_ZERO_THRESHOLD) then + kprev=ktrial + exit + end if + if (ktrial == orderOfMatrix) call Exception_sendWarning ("All the guess eigenvectors are zero! for species "//trim(nameOfSpecies), "at DensityMatrixSCFGuess_getGuess") + end do + do i = 1, orderOfMatrix + do j = 1, orderOfMatrix + densityMatrix%values(i, j) = densityMatrix%values(i, j) + guessMatrix%values(i, ktrial)*guessMatrix%values(j, ktrial) + end do + end do end do densityMatrix%values = densityMatrix%values*MolecularSystem_getEta(speciesID, molSys) diff --git a/src/SCF/MultiSCF.f90 b/src/SCF/MultiSCF.f90 index 449588c8..c7f89783 100644 --- a/src/SCF/MultiSCF.f90 +++ b/src/SCF/MultiSCF.f90 @@ -667,6 +667,7 @@ subroutine MultiSCF_getInitialGuess(this, wfObjects) !! do speciesID = 1, this%molSys%numberOfQuantumSpecies call DensityMatrixSCFGuess_getGuess(speciesID, wfObjects(speciesID)%HcoreMatrix, & + wfObjects(speciesID)%overlapMatrix, & wfObjects(speciesID)%transformationMatrix, & wfObjects(speciesID)%densityMatrix, & wfObjects(speciesID)%waveFunctionCoefficients, & @@ -1248,6 +1249,7 @@ subroutine MultiSCF_showResults(this, wfObjects) end if totalExternalPotentialEnergy = 0.0 + if (CONTROL_instance%IS_THERE_EXTERNAL_POTENTIAL .or. & sum(abs(CONTROL_instance%ELECTRIC_FIELD)) .ne. 0 .or. & CONTROL_instance%ARE_THERE_QDO_POTENTIALS) then diff --git a/src/SCF/WaveFunction.f90 b/src/SCF/WaveFunction.f90 index 15ff90b8..944135c2 100644 --- a/src/SCF/WaveFunction.f90 +++ b/src/SCF/WaveFunction.f90 @@ -558,7 +558,6 @@ subroutine WaveFunction_buildHCoreMatrix(this) auxOmega = MolecularSystem_getOmega(this%species, this%molSys) if (auxOmega .ne. 0.0_8) then - CONTROL_instance%ARE_THERE_QDO_POTENTIALS = .true. this%HCoreMatrix%values = this%HCoreMatrix%values + & (1.0/2.0)*MolecularSystem_getMass(this%species, this%molSys)*auxOmega**2*this%harmonic%values this%externalPotentialMatrix%values = this%externalPotentialMatrix%values + & diff --git a/test/Mu-DHT-euc.lowdin b/test/Mu-DHT-euc.lowdin new file mode 100644 index 00000000..360bd8cd --- /dev/null +++ b/test/Mu-DHT-euc.lowdin @@ -0,0 +1,122 @@ +SYSTEM_DESCRIPTION='Muon functional test - Double harmonic trap from 10.1103/PhysRevB.108.245155' +!This calculation reproduces very well the results reported in the paper +!Value paper openlowdin +!E_ground -0.4547 -0.454704 +!T_e 0.4072 0.407225 +!V_ext_e 0.0007 0.000679 +!T_u 0.0213 0.021334 +!V_exc_u 0.0108 0.010775 +!J_eu -0.8328 -0.832789 +!E_eu_c -0.0619 -0.061929 +!eps_e -0.47 -0.467260 +!eps_u -0.821 -0.820715 + +GEOMETRY +e-(H) 7S7P7D 0.0 0.0 0.0 multiplicity=2 omega = 0.02 +u+ 7S7P7D 0.0 0.0 0.0 omega = 0.02 +n0 dirac 0.0 0.0 0.0 qdoCenterOf=e-alpha +n0 dirac 0.0 0.0 0.0 qdoCenterOf=e-beta +n0 dirac 0.0 0.0 0.0 qdoCenterOf=u+ +END GEOMETRY + +TASKS +method = "UKS" +END TASKS + +CONTROL +readCoefficients=.F. +electronExchangeCorrelationFunctional="NONE" +muonElectronCorrelationFunctional="EMUC-1" +END CONTROL + +BASIS 7S7P7D +O-U+ U+ (SPD) BASIS TYPE: 2 +21 +1 0 1 +15.8560 1.0 +2 0 1 +10.4040 1.0 +3 0 1 +7.0649 1.0 +4 0 1 +5.2960 1.0 +5 0 1 +4.1400 1.0 +6 0 1 +2.0700 1.0 +7 0 1 +1.0350 1.0 +1 1 1 +15.8560 1.0 +2 1 1 +10.4040 1.0 +3 1 1 +7.0649 1.0 +4 1 1 +5.2960 1.0 +5 1 1 +4.1400 1.0 +6 1 1 +2.0700 1.0 +7 1 1 +1.0350 1.0 +1 2 1 +15.8560 1.0 +2 2 1 +10.4040 1.0 +3 2 1 +7.0649 1.0 +4 2 1 +5.2960 1.0 +5 2 1 +4.1400 1.0 +6 2 1 +2.0700 1.0 +7 2 1 +1.0350 1.0 +O-HYDROGEN H (SPD) BASIS TYPE: 1 +21 +1 0 1 +7.2212 1.0 +2 0 1 +3.2589 1.0 +3 0 1 +0.6472 1.0 +4 0 1 +0.0574 1.0 +5 0 1 +1.4523 1.0 +6 0 1 +0.1287 1.0 +7 0 1 +0.2879 1.0 +1 1 1 +7.2212 1.0 +2 1 1 +3.2589 1.0 +3 1 1 +0.6472 1.0 +4 1 1 +0.0574 1.0 +5 1 1 +1.4523 1.0 +6 1 1 +0.1287 1.0 +7 2 1 +0.2879 1.0 +1 2 1 +7.2212 1.0 +2 2 1 +3.2589 1.0 +3 2 1 +0.6472 1.0 +4 2 1 +0.0574 1.0 +5 2 1 +1.4523 1.0 +6 2 1 +0.1287 1.0 +7 2 1 +0.2879 1.0 +END BASIS + diff --git a/test/Mu-DHT-euc.py b/test/Mu-DHT-euc.py new file mode 100644 index 00000000..a3237348 --- /dev/null +++ b/test/Mu-DHT-euc.py @@ -0,0 +1,25 @@ +#!/usr/bin/env python +#The corresponding input file is testName.lowdin +#The functions setReferenceValues and getTestValues are specific for this test +#The common procedures are found in lowdinTestFunctions.py +import sys +import lowdinTestFunctions as test +def setReferenceValues(): + refValues = { +"KS energy" : [-0.454704120729,1E-6], +"U+/E- Corr energy" : [-0.061929601444,1E-4], +"U+/ext pot energy" : [0.010775182673,1E-4], +"E-/ext pot energy" : [0.000679766944,1E-4], +} + return refValues + +def getTestValues(testValues,testName): + testValues["KS energy"] = test.getSCFTotalEnergy(testName) + testValues["U+/E- Corr energy"] = test.getDFTCorrEnergy(testName,"E-ALPHA","U+") + testValues["U+/ext pot energy"] = test.getSCFExtPotEnergy(testName,"U+") + testValues["E-/ext pot energy"] = test.getSCFExtPotEnergy(testName,"E-ALPHA") + return + +if __name__ == '__main__': + testName = sys.argv[0][:-3] + test.performTest(testName,setReferenceValues,getTestValues) diff --git a/test/MuCN-euc-RKS.lowdin b/test/MuCN-euc-RKS.lowdin new file mode 100644 index 00000000..ce4719e5 --- /dev/null +++ b/test/MuCN-euc-RKS.lowdin @@ -0,0 +1,46 @@ +SYSTEM_DESCRIPTION='Molecula de H2' + +GEOMETRY +e-(C) PC-0 0.0 0.0 0.0 +e-(N) PC-0 0.0 0.0 -1.14598 +e-(H) PC-1 0.0 0.0 1.06602 +u+ 5S5P 0.0 0.0 1.06602 +C dirac 0.0 0.0 0.0 +N dirac 0.0 0.0 -1.14598 +END GEOMETRY + +TASKS +method = "RKS" +END TASKS + +CONTROL +readCoefficients=.F. +nonElectronicLevelShifting=0.10 +electronExchangeCorrelationFunctional="B3LYP" +nuclearElectronCorrelationFunctional="EMUC-1" +END CONTROL + +BASIS 5S5P +O-U+ U+ (SP) BASIS TYPE: 2 +10 +1 0 1 +10.4040 1.0 +2 0 1 +7.0649 1.0 +3 0 1 +5.2960 1.0 +4 0 1 +4.1400 1.0 +5 0 1 +2.0700 1.0 +1 1 1 +10.4040 1.0 +2 1 1 +7.0649 1.0 +3 1 1 +5.2960 1.0 +4 1 1 +4.1400 1.0 +5 1 1 +2.0700 1.0 +END BASIS diff --git a/test/MuCN-euc-RKS.py b/test/MuCN-euc-RKS.py new file mode 100644 index 00000000..5481422e --- /dev/null +++ b/test/MuCN-euc-RKS.py @@ -0,0 +1,21 @@ +#!/usr/bin/env python +#The corresponding input file is testName.lowdin +#The functions setReferenceValues and getTestValues are specific for this test +#The common procedures are found in lowdinTestFunctions.py +import sys +import lowdinTestFunctions as test +def setReferenceValues(): + refValues = { +"KS energy" : [-93.071088537483,1E-6], +"U+/E- Corr energy" : [-0.107102067712,1E-4] +} + return refValues + +def getTestValues(testValues,testName): + testValues["KS energy"] = test.getSCFTotalEnergy(testName) + testValues["U+/E- Corr energy"] = test.getDFTCorrEnergy(testName,"E-","U+") + return + +if __name__ == '__main__': + testName = sys.argv[0][:-3] + test.performTest(testName,setReferenceValues,getTestValues) diff --git a/test/MuCN-euc-UKS.lowdin b/test/MuCN-euc-UKS.lowdin new file mode 100644 index 00000000..41620dd5 --- /dev/null +++ b/test/MuCN-euc-UKS.lowdin @@ -0,0 +1,46 @@ +SYSTEM_DESCRIPTION='Molecula de H2' + +GEOMETRY +e-(C) PC-0 0.0 0.0 0.0 +e-(N) PC-0 0.0 0.0 -1.14598 +e-(H) PC-1 0.0 0.0 1.06602 +u+ 5S5P 0.0 0.0 1.06602 +C dirac 0.0 0.0 0.0 +N dirac 0.0 0.0 -1.14598 +END GEOMETRY + +TASKS +method = "UKS" +END TASKS + +CONTROL +readCoefficients=.F. +nonElectronicLevelShifting=0.10 +electronExchangeCorrelationFunctional="B3LYP" +nuclearElectronCorrelationFunctional="EMUC-1" +END CONTROL + +BASIS 5S5P +O-U+ U+ (SP) BASIS TYPE: 2 +10 +1 0 1 +10.4040 1.0 +2 0 1 +7.0649 1.0 +3 0 1 +5.2960 1.0 +4 0 1 +4.1400 1.0 +5 0 1 +2.0700 1.0 +1 1 1 +10.4040 1.0 +2 1 1 +7.0649 1.0 +3 1 1 +5.2960 1.0 +4 1 1 +4.1400 1.0 +5 1 1 +2.0700 1.0 +END BASIS diff --git a/test/MuCN-euc-UKS.py b/test/MuCN-euc-UKS.py new file mode 100644 index 00000000..e7b33037 --- /dev/null +++ b/test/MuCN-euc-UKS.py @@ -0,0 +1,23 @@ +#!/usr/bin/env python +#The corresponding input file is testName.lowdin +#The functions setReferenceValues and getTestValues are specific for this test +#The common procedures are found in lowdinTestFunctions.py +import sys +import lowdinTestFunctions as test +def setReferenceValues(): + refValues = { + "KS energy" : [-93.071088372363,1E-6], + "U+/E-A Corr energy" : [-0.053547483452,1E-4], + "U+/E-B Corr energy" : [-0.053547483452,1E-4] +} + return refValues + +def getTestValues(testValues,testName): + testValues["KS energy"] = test.getSCFTotalEnergy(testName) + testValues["U+/E-A Corr energy"] = test.getDFTCorrEnergy(testName,"E-ALPHA","U+") + testValues["U+/E-B Corr energy"] = test.getDFTCorrEnergy(testName,"E-BETA","U+") + return + +if __name__ == '__main__': + testName = sys.argv[0][:-3] + test.performTest(testName,setReferenceValues,getTestValues)