Multiple-molecule basis optimization

This example demonstrates how BasisOpt can be used to optimize a double-zeta molecule-optimized basis set. In this case five molecules (water, methane, methanol, formaldehyde and oxygen) are used for the optimization, resulting in a single set of exponents for hydrogen, carbon and oxygen. Please note that as the exponents are tightly coupled this optimization can take a long time to run. The full example script can be found in examples/multi-molecule/multi-molecule.py.

This example also shows the use of the logger functionality in BasisOpt.

Loading xyz geometries and choosing the method to use

In this example we use the Psi4 program to carry out the electronic structure calculations. The molecular geometries have been pre-optimized and xyz coordinate files can be found in examples/multi-molecule. The following code creates a set of molecule objects from the xyz files.

mb = MolecularBasis(name="double")
list_of_mols = ['water', 'methane', 'methanol', 'formaldehyde', 'oxygen']
mol_objs = [
bo.molecule.Molecule.from_xyz(mol+'.xyz', name=mol)
for mol in list_of_mols
        ]
mb = MolecularBasis(name="double", molecules=mol_objs)

We also set some calculation parameters that will be used in the optimization strategy. Specifically, we select the \(\omega\) B97X-D density functional and turn off the use of density fitting in Psi4.

params = {
'functional': "wb97x-d",
'scf_type': "pk"
}

Optimization strategy and first run

In this optimization we have elected to use the Karlsruhe def2-SVP basis as the starting guess for the exponents and use def2-QZVP to produce the reference energies. As part of the setup we uncontract the basis set and store the final value of the objective function for later re-use.

strategy = Strategy()
strategy.params = params
strategy.guess_params = {'name': 'def2-svp'}
mb.setup(method='dft', strategy=strategy, params=params, reference='def2-qzvp')
basis = mb.get_basis()
basis = bo.basis.uncontract(basis)

mb.optimize()
e_opt = []
e_opt.append(strategy.last_objective)
e_diff = e_opt[0]

Consistency iterations

As the exponents of the basis set are tightly coupled and because the optimal, for example, hydrogen exponents will depend upon the current carbon and oxygen exponents, we will need to iterate through the optimization until the change in energy is lower than some threshold. These are sometimes known as consistency iterations. In this case we set the threshold to be 1 micro hartree, but this will take a significant length of time to reach convergence.

Within these iterations we also use the logging capability in BasisOpt to produce some basic information that allows the user to monitor progress. This can be especially useful if BasisOpt is used with a batch scheduling system.

conv_crit = 1.0e-6
counter = 0

while e_diff > conv_crit:
    bo_logger.info("Starting consistency iteration %d", counter+1)
    mb.optimize()
    e_opt.append(strategy.last_objective)
    e_diff = strategy.delta_objective
    bo_logger.info("Objective function difference from previous iteration: %f\n", e_diff)
    counter += 1

Writing the basis exponents to file

As the final step in this example, we use the Basis Set Exchange wrapper to write out the optimized exponents in Molpro format.

filename = "opt_basis.txt"
bo_logger.info("Writing optimized basis to %s", filename)
f = open(filename, "x")
f.write(bo.bse_wrapper.internal_basis_converter(mb.get_basis(), fmt='molpro'))
f.close()

The exponents generated by this example can be found in examples/multi-molecule/opt_basis.txt.