I have been trying to replicate some Gaussian 16 DFT calculations in QCHEM with little success. The relative PE between two structures is off by 0.8 kcal/mol when using QCHEM.
My Gaussian 16 settings are: m062x int=grid=ultrafine 6-311G(d,p) scrf=(smd,solvent=tetrahydrofuran)
where ultrafine is a 99,590 grid
The 99,590 grid is specified by XC_GRID=000099000590 since XC_GRID=3 sets SG-3, which is an economized version of 99,590. However, I doubt that will account for all the difference.
Here’s a procedure that can be used to check results between any two codes. Begin by stripping the input from DFT and solvent model.
Compare nuclear repulsion energy to ensure it matches to machine precision between the outputs.
Use METHOD=HF and the same basis, set tight SCF convergence SCF_CONVERGENCE=10 or SCF_CONVERGENCE=12 and tight integral screening THRESH=14. Provide corresponding settings in G16 input and try to match Hartree-Fock energy between the codes. Small differences could indicate that there are additional tricks being used to speed up convergence (need to be disabled) or discrepancies in the definition of the basis set. The latter is a real possibility with Pople and Dunning basis sets, try to use an unambiguously defined basis, for example one from the Ahlrichs def2 family.
Once HF energy is reproduced exactly, turn on DFT by setting METHOD=M06-2X. Ensure that DFT grid settings and cutoffs are identical. Keep in mind that while the number of radial and angular points may match between two programs, different codes may use different radial quadratures. Maximize the density of the grid as much as possible and be prepared to be unable to exactly reproduce the energies with less than very tight grids. The other variable is the definition of the density functional, which may have subtle differences. Try another simpler functional, for example, PBE or PBE0, to confirm.
Finally, add the solvent model. Key parameters that may differ across codes here are the definitions of atomic radii and the number of tessellation points per atom, as well as solvent parametrization. For SMD those should be standard though, and both Q-Chem and Gaussian use the reference implementation of SMD from Truhlar’s group.
If you still have trouble reconciling the results, send full outputs to support@q-chem.com for inspection.
Note that SG-3 uses a different radial quadrature as compared to EML (99,590), so it’s not simply that the former is a subset of the latter. In addition, while SMD uses model-specific atomic radii (which should presumably be the same in both programs), it’s actually not true that the tessellation is the same because SMD uses PCM for its electrostatic component and the PCM discretization procedure does differ between programs.
I suggest trying to match the gas-phase calculation first, as 0.8 kcal/mol could easily be found in discretization error. That error is controllable (and can be systematically eliminated) in Q-Chem using the “HPoints” and “HeavyPoints” keywords in the $pcm section. I believe it’s controllable in the discretization scheme that is uses in Gaussian16 though I don’t know how. Older versions of Gaussian used GEPOL tessellation, which is not really systematically improvable.