How to optimize the structure with XYG3 functional?

[nesquik91]

Hello,

Thanks to google, at last, I found this nice program!
Unfortunately, all documents are written in Chinese and make me really hard to understand through the google translator.
So, it might be already in the documents but could you explain it once again how to optimize the geometry via XYG3 functional?
It will be great to use pyberny solver together.

If I have the H2O coordinate like this (xyz file)
O 0.0000 0.0000 0.0626
H -0.7920 0.0000 -0.4973
H 0.7920 0.0000 -0.4973
charge=0, spin=0, basis=ccpvtz,

then, which geometry is the correct one with the XYG3?
Thank you very much!

suhwan

[ajz34]

Hi~

It seems that I'm not notified by Github. So I saw this issue by your email. Sorry for that 😥
A possible optimized geometry could be below:

O  -0.         0.         0.0790668,
H  -0.7536111  0.        -0.5055334
H   0.7536111  0.        -0.5055334

where XYG3 energy of this molecule is -76.4225154492.
DFT grid of this calculation is (75, 302).

Actually geom optimize is not illustrated in document. These documents mostly emphasis on how gradients are derived and calculated by program.

The python code will be posted in next reply.

[ajz34]

# Import
from pyscf import gto, dft, lib
from pyxdh.DerivOnce import GradXDH
import numpy as np
from berny import Berny, geomlib
from pyscf.geomopt.berny_solver import to_berny_geom, _geom_to_atom
np.set_printoptions(7, suppress=True, linewidth=120)

# Mol Definition
mol = gto.Mole()
mol.atom = """
O 0.0000 0.0000 0.0626
H -0.7920 0.0000 -0.4973
H 0.7920 0.0000 -0.4973
"""
mol.basis = "cc-pVTZ"
mol.verbose = 0
mol.build()
geom = to_berny_geom(mol)

# Energy and Gradient
def solver(geom):
    # Update molecule
    mol.set_geom_(_geom_to_atom(mol, geom, False), unit="Bohr")
    # Generate (75, 302) grids
    grids = dft.Grids(mol)
    grids.atom_grid = (75, 302)
    grids.build()
    # Self-consistent part of XYG3
    scf_eng = dft.RKS(mol)
    scf_eng.xc = "B3LYPg"
    scf_eng.grids = grids
    # Non-self-consistent GGA part of XYG3
    nc_eng = dft.RKS(mol)
    nc_eng.xc = "0.8033*HF - 0.0140*LDA + 0.2107*B88, 0.6789*LYP"
    nc_eng.grids = grids
    # Dipole helper from pyxdh
    config = {
        "scf_eng": scf_eng,
        "nc_eng": nc_eng,
        "cc": 0.3211,
        "cphf_tol": 1e-10
    }
    grad_xDH = GradXDH(config)
    return grad_xDH.eng, grad_xDH.E_1

# Optimization
optimizer = Berny(geom)
for geom in optimizer:
    energy, gradients = solver(geom)
    print("Energy: {:16.10f}, Gradient Norm: {:16.6e}".format(energy, np.linalg.norm(gradients)))
    optimizer.send((energy, gradients))
>>> Energy:   -76.4212217055, Gradient Norm:     3.007911e-02
>>> Energy:   -76.4224930655, Gradient Norm:     6.383016e-03
>>> Energy:   -76.4225147540, Gradient Norm:     8.547240e-04
>>> Energy:   -76.4225154492, Gradient Norm:     6.597390e-06

# Final Coordinate
mol.atom_coords() * lib.param.BOHR
>>> array([[-0.       ,  0.       ,  0.0790668],
>>>        [-0.7536111,  0.       , -0.5055334],
>>>        [ 0.7536111,  0.       , -0.5055334]])

[nesquik91]

Thank you very much!
I got the very same result 👍
One more thing is, any specific reason for using 75,302 grids?

[ajz34]

No specific reason
Just because the efficiency of program is somehow poor, (99, 590) grids could be time consumable (even for water) 0.0

[nesquik91]

Hello ajz34,
Is it also possible to perform geometric optimization for the transition state with pyxdh?

[ajz34]

Actually I havn't tried that yet. Maybe in principle this could be possible.

[nesquik91]

got it. Thanks for the quick reply.