This code snippet shows how to build a bais set for a given molecule, simple and fast

from iodata import load_one
import basis_set_exchange as bse

mol = load_one("sample.sdf")

# get basis set information
selected_basis_set = bse.get_basis("def2-svpd",
                                   fmt="gaussian94",
                                   elements=mol.atnums.tolist(),
                                   header=False)

https://molssi-bse.github.io/basis_set_exchange/

Thanks for keeping track of this. I believe for Gaussian, when using the gen keyword, the basis set for each element should only be listed once. We may needs something like np.unique(mol.natoms).tolist().

Thanks for keeping track of this. I believe for Gaussian, when using the gen keyword, the basis set for each element should only be listed once. We may needs something like np.unique(mol.natoms).tolist().

Nice catch! We should use a list without duplicated elements (sort of set).

Good to know np.unique() instead of list(set(my_list)). Thanks! @tovrstra

I did some testing and extended the example to work out a specific case. This is just to get a better idea of how to proceed.

Python script:

from iodata import load_one, write_input
import basis_set_exchange as bse

mol = load_one("formamide.sdf")

basis_set = bse.get_basis(
    "def2-svpd",
    fmt="gaussian94",
    elements=mol.atnums.tolist(),
    header=False,
)

template = """\
#n {lot}/Gen {run_type}

{title}

{charge} {spinmult}
{geometry}

{basis_set}
"""

write_input(mol, "gaussian.com", "gaussian", template, basis_set=basis_set)

The SDF file formamide.sdf:

713
  -OEChem-03292115373D

  6  5  0     0  0  0  0  0  0999 v2000
    1.1280    0.2091    0.0000 O   0  0  0  0  0  0  0  0  0  0  0  0
   -1.1878    0.1791    0.0000 N   0  0  0  0  0  0  0  0  0  0  0  0
    0.0598   -0.3882    0.0000 C   0  0  0  0  0  0  0  0  0  0  0  0
   -1.3085    1.1864    0.0001 H   0  0  0  0  0  0  0  0  0  0  0  0
   -2.0305   -0.3861   -0.0001 H   0  0  0  0  0  0  0  0  0  0  0  0
   -0.0014   -1.4883   -0.0001 H   0  0  0  0  0  0  0  0  0  0  0  0
  1  3  2  0  0  0  0
  2  3  1  0  0  0  0
  2  4  1  0  0  0  0
  2  5  1  0  0  0  0
  3  6  1  0  0  0  0
M  END
$$$$

This results in a valid Gaussian input file. The bse Python package already knows that it should not repeat basis functions for the same element, which is convenient. @fwmeng88 Does this also work for you?

I'm not sure yet if we should build this into IOData, or just make sure it remains possible for users to use this approach in their own scripts. (We should at least add it as example to the documentation.) Implementing this inside IOData has the obvious advantage of convenience, but comes with some minor drawbacks as well: an extra dependency & (very minor) another option to control if built-in or bse basis sets should be used. To facilitate this decision: did I miss any pros or cons?

At least from my point of view, adding the increasingly standard dependency on the basis-set-exchange would be outweighed by increased convenience.

Good. It think it would become a much-loved feature. Should it be the default behavior?

By the way, something similar is possible with ORCA, see https://sites.google.com/site/orcainputlibrary/basis-sets (see newgto command).

I'd be in favour of using BSE to do this, especially since it already prepares the formatted lines, so we could just put it right into the template.

If we really wanted to be able to generate the non-standard basis from MolecularBasis, it wouldn't be terribly hard to add a load_one for their BSE JSON format (it's similar to QCSchema's basis schema), but then we'd need to build a converter for every format for which we want to support this feature, which seems like a lot of work without much benefit.

@wilhadams I hadn't thought yet of the second possibility you mention. It might be useful in some corner cases, e.g. when you have a MolecularBasis instance from a wavefunction file for which the name is not known or available, but I would indeed not worry about it for now.

Adding the basis_set_exchange would be a very nice feature to build Gaussian input files. I have tested the code you proposed and it worked perfectly on my end. @tovrstra

Speaking of putting an example in the documentation, maybe we need to put more than one? Such as one example for regular geometry optimization and another to show how to handle multi-step jobs (https://gaussian.com/input/). This can show how useful IOData can be for building Gaussian input files.

I've added it on the list for 1.0, to fix in this in relatively short term. A few related issues need to be handled before this one: #253 to avoid merge conflicts and #254 and #210 to provide a better infrastructure and place for more elaborate examples in the documentation.