Skip to content

muon-electron correlation functional from Goli and Shahbazian 2022 #177

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fsmoncadaa merged 1 commit into from
Aug 6, 2026
Merged

muon-electron correlation functional from Goli and Shahbazian 2022 #177

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
122 changes: 122 additions & 0 deletions examples/RKS_emuc_Goli2022/MuCN-euc.lowdin
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
117 changes: 117 additions & 0 deletions examples/RKS_emuc_Goli2022/MuCN-noeuc.lowdin
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
8 changes: 8 additions & 0 deletions src/Core/CONTROL.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -951,6 +953,7 @@ module CONTROL_
LowdinParameters_electronExchangeCorrelationFunctional, &
LowdinParameters_nuclearElectronCorrelationFunctional, &
LowdinParameters_positronElectronCorrelationFunctional, &
LowdinParameters_muonElectronCorrelationFunctional, &
LowdinParameters_betaFunction, &
LowdinParameters_gridRadialPoints, &
LowdinParameters_gridAngularPoints, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
6 changes: 6 additions & 0 deletions src/Core/MolecularSystem.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
Loading
Loading