This tutorial will go over general QCFractal usage to give a feel for the ecosystem. In this tutorial, we employ Snowflake, a simple QCFractal stack which runs on a local machine for demonstration and exploration purposes.
To begin this quickstart tutorial, first install the QCArchive Snowflake environment from conda:
conda env create qcarchive/qcfractal-snowflake -n snowflake
conda activate snowflakeIf you have a pre-existing environment with qcfractal, ensure that rdkit and geometric are installed from the conda-forge channel and psi4 and dftd3 from the psi4 channel.
First let us import two items from the ecosystem:
FractalSnowflakeHandler - This is a FractalServer that is temporary and is used for trying out new things.
qcfractal.interface is the QCPortal module, but if using QCFractal it is best to import it locally.
Typically we alias qcportal as ptl. We will do the same for qcfractal.interface so that the code can be used anywhere.
In [1]:
from qcfractal import FractalSnowflakeHandler
import qcfractal.interface as ptl
We can now build a temporary server which acts just like a normal server, but we have a bit more direct control of it.
Warning! All data is lost when this notebook shuts down! This is for demonstration purposes only! For information about how to setup a permanent QCFractal server, see the Setup Quickstart Guide.
In [2]:
server = FractalSnowflakeHandler()
server
Out[2]:
We can then build a typical FractalClient to automatically connect to this server using the client() helper command. Note that the server names and addresses are identical in both the server and client.
In [3]:
client = server.client()
client
Out[3]:
A server starts with no data, so let's add some! We can do this by adding a water molecule at a poor geometry from XYZ coordinates.
Note that all internal QCFractal values are stored and used in atomic units;
whereas, the standard Molecule.from_data() assumes an input of Angstroms.
We can switch this back to Bohr by adding a units command in the text string.
In [4]:
mol = ptl.Molecule.from_data("""
O 0 0 0
H 0 0 2
H 0 2 0
units bohr
""")
mol
We can then measure various aspects of this molecule to determine its shape. Note that the measure command will provide a distance, angle, or dihedral depending if 2, 3, or 4 indices are passed in.
This molecule is quite far from optimal, so let's run an geometry optimization!
In [5]:
print(mol.measure([0, 1]))
print(mol.measure([1, 0, 2]))
We originally installed psi4 and geometric, so we can use these programs to perform a geometry optimization. In QCFractal, we call a geometry optimization a procedure, where procedure is a generic term for a higher level operation that will run multiple individual quantum chemistry energy, gradient, or Hessian evaluations. Other procedure examples are finite-difference computations, n-body computations, and torsiondrives.
We provide a JSON-like input to the client.add_procedure()
command to specify the method, basis, and program to be used.
The qc_spec field is used in all procedures to determine the underlying quantum chemistry method behind the individual procedure.
In this way, we can use any program or method that returns a energy or gradient quantity to run our geometry optimization!
(See also add_compute().)
In [6]:
spec = {
"keywords": None,
"qc_spec": {
"driver": "gradient",
"method": "b3lyp",
"basis": "6-31g",
"program": "psi4"
},
}
# Ask the server to compute a new computation
r = client.add_procedure("optimization", "geometric", spec, [mol])
print(r)
print(r.ids)
We can see that we submitted a single task to be evaluated and the server has not seen this particular procedure before.
The ids field returns the unique id of the procedure. Different procedures will always have a unique id, while identical procedures will always return the same id.
We can submit the same procedure again to see this effect:
In [7]:
r2 = client.add_procedure("optimization", "geometric", spec, [mol])
print(r)
print(r.ids)
Once a task is submitted, it will be placed in the compute queue and evaluated. In this particular case the FractalSnowflakeHandler uses your local hardware to evaluate these jobs. We recommend avoiding large tasks!
In general, the server can handle anywhere between laptop-scale resources to many hundreds of thousands of concurrent cores at many physical locations. The amount of resources to connect is up to you and the amount of compute that you require.
Since we did submit a very small job it is likely complete by now. Let us query this procedure from the server using its id like so:
In [15]:
proc = client.query_procedures(id=r.ids)[0]
proc
Out[15]:
This OptimizationRecord object has many different fields attached to it so that all quantities involved in the computation can be explored. For this example, let us pull the final molecule (optimized structure) and inspect the physical dimensions.
Note: if the status does not say COMPLETE, these fields will not be available. Try querying the procedure again in a few seconds to see if the task completed in the background.
In [16]:
final_mol = proc.get_final_molecule()
In [17]:
print(final_mol.measure([0, 1]))
print(final_mol.measure([1, 0, 2]))
final_mol
This water molecule has bond length and angle dimensions much closer to expected values. We can also plot the optimization history to see how each step in the geometry optimization affected the results. Though the chart is not too impressive for this simple molecule, it is hopefully illuminating and is available for any geometry optimization ever completed.
In [18]:
proc.show_history()
Submitting individual procedures or single quantum chemistry tasks is not typically done as it becomes hard to track individual tasks. To help resolve this, Collections are different ways of organizing standard computations so that many tasks can be referenced in a more human-friendly way. In this particular case, we will be exploring an intermolecular potential dataset.
To begin, we will create a new dataset and add a few intermolecular interactions to it.
In [19]:
ds = ptl.collections.ReactionDataset("My IE Dataset", ds_type="ie", client=client, default_program="psi4")
We can construct a water dimer that has fragments used in the intermolecular computation with the -- divider. A single water molecule with ghost atoms can be extracted like so:
In [20]:
water_dimer = ptl.Molecule.from_data("""
O 0.000000 0.000000 0.000000
H 0.758602 0.000000 0.504284
H 0.260455 0.000000 -0.872893
--
O 3.000000 0.500000 0.000000
H 3.758602 0.500000 0.504284
H 3.260455 0.500000 -0.872893
""")
water_dimer.get_fragment(0, 1)
Many molecular entries can be added to this dataset where each is entry is a given intermolecular complex that is given a unique name. In addition, the add_ie_rxn method to can automatically fragment molecules.
In [21]:
ds.add_ie_rxn("water dimer", water_dimer)
ds.add_ie_rxn("helium dimer", """
He 0 0 -3
--
He 0 0 3
""")
Out[21]:
Once the Collection is created it can be saved to the server so that it can always be retrived at a future date
In [22]:
ds.save()
Out[22]:
The client can list all Collections currently on the server and retrive collections to be manipulated:
In [23]:
client.list_collections()
Out[23]:
In [24]:
ds = client.get_collection("ReactionDataset", "My IE Dataset")
ds
Out[24]:
In [25]:
ds.compute("B3LYP-D3", "def2-SVP")
Out[25]:
By default this collection evaluates the non-counterpoise corrected interaction energy which typically requires three computations per entry (the complex and each monomer). In this case we compute the B3LYP and -D3 additive correction separately, nominally 12 total computations. However the collection is smart enough to understand that each Helium monomer is identical and does not need to be computed twice, reducing the total number of computations to 10 as shown here. We can continue to compute additional methods. Again, this is being evaluated on your computer! Be careful of the compute requirements.
In [26]:
ds.compute("PBE-D3", "def2-SVP")
Out[26]:
A list of all methods that have been computed for this dataset can also be shown:
In [28]:
ds.list_values()
Out[28]:
The above only shows what has been computed and does not pull this data from the server to your computer. To do so, the get_history command can be used:
In [32]:
print(f"DataFrame units: {ds.units}")
ds.get_values()
Out[32]:
You can also visualize results and more!
In [33]:
ds.visualize(["B3LYP-D3", "PBE-D3"], "def2-SVP", bench="B3LYP/def2-svp", kind="violin")
This is only the very beginning of what you can do with QCFractal! Explore the documentation to learn more capabilities. In particular, see the next section for a quickstart guide on how to set up QCFractal in production.
If you like the project, consider starring us on GitHub. If you have any questions, join our Slack channel.