def run_pyscf(xyzfile, charge, mult, do_grad,hess=False):
    import numpy as np
    from pyscf.dft import rks
    from pyscf import  grad

    atoms = read_xyz(xyzfile)

    mol = gto.M(
        atom=atoms,
        basis={
            'C': 'aug-cc-pvdz',
            'N': 'aug-cc-pvdz',
            'H': 'aug-cc-pvdz',
            'Br': 'aug-cc-pvdz',
            'O': 'aug-cc-pvdz',
            'Mn': gto.basis.load('stuttgart_rsc', 'Mn')
        },
        ecp={'Mn': gto.load_ecp('stuttgart_rsc', 'Mn')},
        unit="angstrom",
        charge=charge,
        spin=mult-1,
        verbose=3
    )

    # Create GPU RKS object
    mf = rks.RKS(mol, xc='M06L').density_fit()

    # Optional solvent (may need CPU fallback if not supported on GPU)
    mf = mf.SMD()
    mf.with_solvent.method = 'SMD'
    mf.with_solvent.solvent = '1,4-dioxane'

    mf.max_cycle = 1000
    mf.grids.level = 8
    mf.conv_tol = 1e-8

    # Energy calculation
    energy = mf.kernel()

    # Gradient calculation using GPU4PySCF
    if do_grad:
        grad_calc = mf.nuc_grad_method().kernel()  # GPU-compatible
    else:
        grad_calc = np.zeros((mol.natm, 3))

    if hess:
        hess_matx=mf.Hessian().kernel()
    else:
        hess_matx = np.zeros((3*mol.natm, 3*mol.natm))
    return energy, grad_calc, hess_matx