search

dftlibs / xcfun

XCFun: A library of exchange-correlation functionals with arbitrary-order derivatives
https://dftlibs.org/xcfun/
Mozilla Public License 2.0
57 stars 32 forks source link

VWN3 Open-Shell could not reproduce Gaussian form of VWN3 #132

Closed ajz34 closed 4 years ago

ajz34 commented 4 years ago

When I was using PySCF with xcfun engine (default PySCF engine is libxc), the closed-shell B3LYP energy (Gaussian convention) can be reproduced, but open-shell could not. I think this problem originates from xcfun's definition of VWN3, so I raised issue here. Code below is almost PySCF. I'm not accustomed to xcfun debugging (T.T)

Possible Solution

In file vwn.hpp, possibly change definition of dd as

num dd = g * ((vwn_f(s, ferro) - vwn_f(s, para));

This could possibly reproduce VWN3 of Gaussian, or VWN_RPA of libxc. (So I don't know why it seems that libxc's VWN3 is different to that of Gaussian's, but libxc's VWN_RPA is the same to Gaussian's VWN3 ...)

Reproduction details

Below is some lengthy reproduction details. For radical CH3, Gaussian input card

#p B3LYP/6-31G Int(Grid=99590)

CH3

0 2
C  0. 0. 0.
H  0. 0. 1.
H  0. 1. 0.
H  1. 0. 0.

should give SCF Done value like -39.7370512431. Almost same value can be reproduced by PySCF with default dft engine (libxc):

def calc_eng(xc):
    from pyscf import gto, scf, dft
    mol = gto.Mole()
    mol.atom = """
    C  0. 0. 0.
    H  0. 0. 1.
    H  0. 1. 0.
    H  1. 0. 0.
    """
    mol.basis = "6-31G"
    mol.spin = 1
    mol.verbose = 0
    mol.build()
    grids = dft.Grids(mol)
    grids.atom_grid = (99, 590)
    grids.prune = None
    grids.radi_method = dft.gauss_chebyshev
    grids.becke_scheme = dft.gen_grid.stratmann
    grids.build()
    scf_eng = dft.UKS(mol)
    scf_eng.xc = xc
    scf_eng.grids = grids
    return scf_eng.run().e_tot
calc_eng("B3LYPg")  # -39.73705123814413

But change numint.py in PySCF as

# --- originally ---
# try:
#     from pyscf.dft import libxc
# --- change to ---
try:
    from pyscf.dft import xcfun
    libxc = xcfun

Then execuate the python code above. The returned energy is about -39.737219767297034.

Above is unrestricted calculation. If molecule's alpha density equals to beta density, Gaussian, xcfun, libxc should be the same.

Your Environment

Configuration is the same as https://github.com/dftlibs/xcfun/issues/131.

uekstrom commented 4 years ago

Dear Zhenyu Zhu, there are many implementations of the VWN functional. What is the reference for your solution?

Regards, Ulf

Den lör 6 juni 2020 kl 07:04 skrev Zhenyu Zhu ajz34 < notifications@github.com>:

When I was using PySCF with xcfun engine (default PySCF engine is libxc), the closed-shell B3LYP energy (Gaussian convention) can be reproduced, but open-shell could not. I think this problem originates from xcfun's definition of VWN3, so I raised issue here. Code below is almost PySCF. I'm not accustomed to xcfun debugging (T.T) Possible Solution

In file vwn.hpp https://github.com/dftlibs/xcfun/blob/master/src/functionals/vwn.hpp#L74, possibly change definition of dd as

num dd = g * ((vwn_f(s, ferro) - vwn_f(s, para));

This could possibly reproduce VWN3 of Gaussian, or VWN_RPA of libxc https://gitlab.com/libxc/libxc/-/blob/master/maple/lda_exc/lda_c_vwn_rpa.mpl . (So I don't know why it seems that libxc's VWN3 is different to that of Gaussian's, but libxc's VWN_RPA is the same to Gaussian's VWN3 ...) Reproduction details

Below is some lengthy reproduction details. For radical CH3, Gaussian input card

p B3LYP/6-31G Int(Grid=99590)

CH3

0 2 C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0.

should give SCF Done value like -39.7370512431. Almost same value can be reproduced by PySCF with default dft engine (libxc):

def calc_eng(xc): from pyscf import gto, scf, dft mol = gto.Mole() mol.atom = """ C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0. """ mol.basis = "6-31G" mol.spin = 1 mol.verbose = 0 mol.build() grids = dft.Grids(mol) grids.atom_grid = (99, 590) grids.prune = None grids.radi_method = dft.gauss_chebyshev grids.becke_scheme = dft.gen_grid.stratmann grids.build() scf_eng = dft.UKS(mol) scf_eng.xc = xc scf_eng.grids = grids return scf_eng.run().e_totcalc_eng("B3LYPg") # -39.73705123814413

But change numint.py https://github.com/pyscf/pyscf/blob/master/pyscf/dft/numint.py#L25 in PySCF as

--- originally ---# try:# from pyscf.dft import libxc# --- change to ---try:

from pyscf.dft import xcfun
libxc = xcfun

Then execuate the python code above. The returned energy is about -39.737219767297034.

Above is unrestricted calculation. If molecule's alpha density equals to beta density, Gaussian, xcfun, libxc should be the same. Your Environment

Configuration is the same as #131 https://github.com/dftlibs/xcfun/issues/131.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/dftlibs/xcfun/issues/132, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADD25GONPONUL5FKMBXU5ZDRVHE6RANCNFSM4NVS5L2Q .

uekstrom commented 4 years ago

Dear Zhenyu Zhu, I looked a bit at the Molpro documentation and it looks like you may have a point, but I don't know if your fix is correct. Could one of the current maintainers check with the VWN3 paper if the implementation is correct? Here's the Molpro formula: https://www.molpro.net/info/release/doc/manual/node1065.html If XCFun is wrong then Dalton (at least the version I have) is also wrong. Possibly we only tested using VWN5. Here's a link about using that functional (the recommended one) in Gaussian: https://sites.google.com/site/rangsiman1993/comp-chem/tips-and-tricks/gaussian-b3lyp

Regards, Ulf

Den lör 6 juni 2020 kl 21:19 skrev uekstrom@gmail.com uekstrom@gmail.com:

Dear Zhenyu Zhu, there are many implementations of the VWN functional. What is the reference for your solution?

Regards, Ulf

Den lör 6 juni 2020 kl 07:04 skrev Zhenyu Zhu ajz34 < notifications@github.com>:

When I was using PySCF with xcfun engine (default PySCF engine is libxc), the closed-shell B3LYP energy (Gaussian convention) can be reproduced, but open-shell could not. I think this problem originates from xcfun's definition of VWN3, so I raised issue here. Code below is almost PySCF. I'm not accustomed to xcfun debugging (T.T) Possible Solution

In file vwn.hpp https://github.com/dftlibs/xcfun/blob/master/src/functionals/vwn.hpp#L74, possibly change definition of dd as

num dd = g * ((vwn_f(s, ferro) - vwn_f(s, para));

This could possibly reproduce VWN3 of Gaussian, or VWN_RPA of libxc https://gitlab.com/libxc/libxc/-/blob/master/maple/lda_exc/lda_c_vwn_rpa.mpl . (So I don't know why it seems that libxc's VWN3 is different to that of Gaussian's, but libxc's VWN_RPA is the same to Gaussian's VWN3 ...) Reproduction details

Below is some lengthy reproduction details. For radical CH3, Gaussian input card

p B3LYP/6-31G Int(Grid=99590)

CH3

0 2 C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0.

should give SCF Done value like -39.7370512431. Almost same value can be reproduced by PySCF with default dft engine (libxc):

def calc_eng(xc): from pyscf import gto, scf, dft mol = gto.Mole() mol.atom = """ C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0. """ mol.basis = "6-31G" mol.spin = 1 mol.verbose = 0 mol.build() grids = dft.Grids(mol) grids.atom_grid = (99, 590) grids.prune = None grids.radi_method = dft.gauss_chebyshev grids.becke_scheme = dft.gen_grid.stratmann grids.build() scf_eng = dft.UKS(mol) scf_eng.xc = xc scf_eng.grids = grids return scf_eng.run().e_totcalc_eng("B3LYPg") # -39.73705123814413

But change numint.py https://github.com/pyscf/pyscf/blob/master/pyscf/dft/numint.py#L25 in PySCF as

--- originally ---# try:# from pyscf.dft import libxc# --- change to ---try:

from pyscf.dft import xcfun
libxc = xcfun

Then execuate the python code above. The returned energy is about -39.737219767297034.

Above is unrestricted calculation. If molecule's alpha density equals to beta density, Gaussian, xcfun, libxc should be the same. Your Environment

Configuration is the same as #131 https://github.com/dftlibs/xcfun/issues/131.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/dftlibs/xcfun/issues/132, or unsubscribe https://github.com/notifications/unsubscribe-auth/ADD25GONPONUL5FKMBXU5ZDRVHE6RANCNFSM4NVS5L2Q .

robertodr commented 4 years ago

I can have a look in roughly a week's time

On Sat, Jun 6, 2020, 21:32 Ulf Ekström notifications@github.com wrote:

Dear Zhenyu Zhu, I looked a bit at the Molpro documentation and it looks like you may have a point, but I don't know if your fix is correct. Could one of the current maintainers check with the VWN3 paper if the implementation is correct? Here's the Molpro formula: https://www.molpro.net/info/release/doc/manual/node1065.html If XCFun is wrong then Dalton (at least the version I have) is also wrong. Possibly we only tested using VWN5. Here's a link about using that functional (the recommended one) in Gaussian:

https://sites.google.com/site/rangsiman1993/comp-chem/tips-and-tricks/gaussian-b3lyp

Regards, Ulf

Den lör 6 juni 2020 kl 21:19 skrev uekstrom@gmail.com <uekstrom@gmail.com

:

Dear Zhenyu Zhu, there are many implementations of the VWN functional. What is the reference for your solution?

Regards, Ulf

Den lör 6 juni 2020 kl 07:04 skrev Zhenyu Zhu ajz34 < notifications@github.com>:

When I was using PySCF with xcfun engine (default PySCF engine is libxc), the closed-shell B3LYP energy (Gaussian convention) can be reproduced, but open-shell could not. I think this problem originates from xcfun's definition of VWN3, so I raised issue here. Code below is almost PySCF. I'm not accustomed to xcfun debugging (T.T) Possible Solution

In file vwn.hpp < https://github.com/dftlibs/xcfun/blob/master/src/functionals/vwn.hpp#L74>, possibly change definition of dd as

num dd = g * ((vwn_f(s, ferro) - vwn_f(s, para));

This could possibly reproduce VWN3 of Gaussian, or VWN_RPA of libxc < https://gitlab.com/libxc/libxc/-/blob/master/maple/lda_exc/lda_c_vwn_rpa.mpl

. (So I don't know why it seems that libxc's VWN3 is different to that of Gaussian's, but libxc's VWN_RPA is the same to Gaussian's VWN3 ...) Reproduction details

Below is some lengthy reproduction details. For radical CH3, Gaussian input card

p B3LYP/6-31G Int(Grid=99590)

CH3

0 2 C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0.

should give SCF Done value like -39.7370512431. Almost same value can be reproduced by PySCF with default dft engine (libxc):

def calc_eng(xc): from pyscf import gto, scf, dft mol = gto.Mole() mol.atom = """ C 0. 0. 0. H 0. 0. 1. H 0. 1. 0. H 1. 0. 0. """ mol.basis = "6-31G" mol.spin = 1 mol.verbose = 0 mol.build() grids = dft.Grids(mol) grids.atom_grid = (99, 590) grids.prune = None grids.radi_method = dft.gauss_chebyshev grids.becke_scheme = dft.gen_grid.stratmann grids.build() scf_eng = dft.UKS(mol) scf_eng.xc = xc scf_eng.grids = grids return scf_eng.run().e_totcalc_eng("B3LYPg") # -39.73705123814413

But change numint.py https://github.com/pyscf/pyscf/blob/master/pyscf/dft/numint.py#L25 in PySCF as

--- originally ---# try:# from pyscf.dft import libxc# --- change to

---try: from pyscf.dft import xcfun libxc = xcfun

Then execuate the python code above. The returned energy is about -39.737219767297034.

Above is unrestricted calculation. If molecule's alpha density equals to beta density, Gaussian, xcfun, libxc should be the same. Your Environment

Configuration is the same as #131 https://github.com/dftlibs/xcfun/issues/131.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/dftlibs/xcfun/issues/132, or unsubscribe < https://github.com/notifications/unsubscribe-auth/ADD25GONPONUL5FKMBXU5ZDRVHE6RANCNFSM4NVS5L2Q

.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/dftlibs/xcfun/issues/132#issuecomment-640107660, or unsubscribe https://github.com/notifications/unsubscribe-auth/AA4JOEND3NJT52DPHTMF7E3RVKKVTANCNFSM4NVS5L2Q .

bast commented 4 years ago

I will also have a look (I am now in the context because of a little side-project where the definition of VWN).

bast commented 4 years ago

I am finally looking at this issue ...

bast commented 4 years ago

After some tests I also confirm that with that suggested change I reproduce lda_c_vwn_rpa of libxc in the spin-polarized case. As Zhenyu Zhu indicated, the problem does not show up if na == nb.

I will some more testing but will commit a patch soon.