|
| 1 | +---------------------------- |
| 2 | +Running MRChem with QCEngine |
| 3 | +---------------------------- |
| 4 | + |
| 5 | +MRChem >=1.0 can be used as a computational engine with the `QCEngine |
| 6 | +<http://docs.qcarchive.molssi.org/projects/qcengine/>`_ program executor. |
| 7 | +QCEngine can be useful for running calculations on large sets of molecules and input parameters. |
| 8 | +The results are collected in standardised `QCScheme format |
| 9 | +<http://molssi-qc-schema.readthedocs.io/en/latest/index.html#>`_, which makes it |
| 10 | +easy to build post-processing pipelines and store data according to Findability, |
| 11 | +Accessibility, Interoperability, and Reuse (FAIR) of digital assets principles. |
| 12 | +Furthermore, QCEngine provides different geometry optimization drivers that can |
| 13 | +use the molecular gradient computed by MRChem for structural optimization. |
| 14 | + |
| 15 | +Installation |
| 16 | +------------ |
| 17 | + |
| 18 | +The easiest way is to install both QCEngine and MRChem in a Conda environment |
| 19 | +using the precompiled version: |
| 20 | + |
| 21 | +.. code-block:: bash |
| 22 | +
|
| 23 | + conda create -n mrchem-qcng mrchem qcengine qcelemental geometric optking pip -c conda-forge |
| 24 | + conda activate mrchem-qcng |
| 25 | + python -m pip install -U pyberny |
| 26 | +
|
| 27 | +It is also possible to use your own installation of MRChem: just make sure that |
| 28 | +the installation folder is in your ``PATH``. |
| 29 | + |
| 30 | + |
| 31 | +.. note:: |
| 32 | + |
| 33 | + If you want to use the precompiled, MPI-parallel version of MRChem with |
| 34 | + OpenMPI, install ``mrchem=*=*openmpi*`` insted of just ``mrchem``. |
| 35 | + A binary package compiled against MPICH is also available: |
| 36 | + ``mrchem=*=*mpich*``. |
| 37 | + |
| 38 | + |
| 39 | +Single compute |
| 40 | +~~~~~~~~~~~~~~ |
| 41 | + |
| 42 | +Calculations in QCEngine are defined in Python scripts. For example, the following runs MRChem to obtain the energy of water: |
| 43 | + |
| 44 | +.. code-block:: python |
| 45 | +
|
| 46 | + import qcelemental as qcel |
| 47 | + import qcengine as qcng |
| 48 | +
|
| 49 | +
|
| 50 | + mol = qcel.models.Molecule(geometry=[[0, 0, 0], [0, 1.5, 0], [0, 0, 1.5]], |
| 51 | + symbols=["O", "H", "H"], |
| 52 | + connectivity=[[0, 1, 1], [0, 2, 1]]) |
| 53 | + print(mol) |
| 54 | +
|
| 55 | + computation = { |
| 56 | + "molecule": mol, |
| 57 | + "driver": "energy", |
| 58 | + "model": {"method": "HF"}, |
| 59 | + "keywords": {"world_prec": 1.0e-3}, |
| 60 | + } |
| 61 | + ret = qcng.compute(computation, "mrchem") |
| 62 | +
|
| 63 | + print(f"E_HF = {ret.return_result} Hartree") |
| 64 | +
|
| 65 | +You can save this sample as `mrchem-run-hf.py` and execute it with: |
| 66 | + |
| 67 | +.. code-block:: bash |
| 68 | +
|
| 69 | + python mrchem-run-hf.py |
| 70 | +
|
| 71 | +Which will print to screen: |
| 72 | + |
| 73 | +.. code-block:: bash |
| 74 | +
|
| 75 | + Molecule(name='H2O', formula='H2O', hash='b41d0c5') |
| 76 | + E_HF = -75.9789291596064 Hartree |
| 77 | +
|
| 78 | +Note that: |
| 79 | + |
| 80 | +#. The molecule is specified, in Angstrom, using a `QCElemental |
| 81 | + <http://docs.qcarchive.molssi.org/projects/qcelemental/en/latest/>`_ object. |
| 82 | +#. The computation is described using a Python dictionary. |
| 83 | +#. The ``driver`` selects the kind of calculation you want to run with MRChem. |
| 84 | + Available drivers are: |
| 85 | + |
| 86 | + - ``energy``, for single-point energy calculations. |
| 87 | + - ``gradient``, for evaluation of the molecular gradient at a given |
| 88 | + geometry. |
| 89 | + - ``properties``, for the calculation of molecular properties. |
| 90 | +#. The ``model`` selects the wavefunction: ``HF`` for Hartree-Fock and any of |
| 91 | + the DFT functionals known to MRChem for a corresponding DFT calculation. |
| 92 | +#. The ``keywords`` key in the dictionary accepts a dictionary of MRChem |
| 93 | + options. Any of the options in the usual input file are recognized. |
| 94 | + |
| 95 | +Once you have a dictionary defining your computation, you can run it with: |
| 96 | + |
| 97 | +.. code-block:: python |
| 98 | +
|
| 99 | + ret = qcng.compute(computation, "mrchem") |
| 100 | +
|
| 101 | +You can reuse the same dictionary with *multiple* computational engine, *e.g.* |
| 102 | +other quantum chemistry programs that are recognized as executors by QCEngine. |
| 103 | +The return value from the ``compute`` function contains all data produced |
| 104 | +during the calculation in QCSchema format including, for example, the execution |
| 105 | +time elapsed. The full JSON output produced by MRChem is also available and can |
| 106 | +be inspected in Python as: |
| 107 | + |
| 108 | +.. code-block:: python |
| 109 | +
|
| 110 | + mrchem_json_out = ret.extras["raw_output"]["output"] |
| 111 | +
|
| 112 | +The full, human-readable input is saved as the ``stdout`` property of the |
| 113 | +object returned by ``compute``. |
| 114 | + |
| 115 | +Parallelism |
| 116 | +~~~~~~~~~~~ |
| 117 | + |
| 118 | +QCEngine allows you to exploit available parallel hardware. |
| 119 | +For example, to use 20 OpenMP threads in your MRChem calculation you would |
| 120 | +provide an additional task configuration dictionary as a ``task_config`` |
| 121 | +argument to ``compute``: |
| 122 | + |
| 123 | +.. code-block:: python |
| 124 | +
|
| 125 | + ret = qcng.compute( |
| 126 | + computation, |
| 127 | + "mrchem", |
| 128 | + task_config={"ncores": 20}) |
| 129 | +
|
| 130 | +You can inspect how the job was launched by printing out the ``provenance`` dictionary: |
| 131 | + |
| 132 | +.. code-block:: python |
| 133 | +
|
| 134 | + print(ret.extras["raw_output"]["output"]["provenance"]) |
| 135 | +
|
| 136 | +.. code-block:: bash |
| 137 | +
|
| 138 | + { |
| 139 | + "creator": "MRChem", |
| 140 | + "mpi_processes": 1, |
| 141 | + "routine": "/home/roberto/miniconda3/envs/mrchem-qcng/bin/mrchem.x", |
| 142 | + "total_cores": 1, |
| 143 | + "version": "1.1.0", |
| 144 | + "ncores": 12, |
| 145 | + "nnodes": 1, |
| 146 | + "ranks_per_node": 1, |
| 147 | + "cores_per_rank": 12, |
| 148 | + "total_ranks": 1 |
| 149 | + } |
| 150 | +
|
| 151 | + |
| 152 | +It is also possible to run MPI-parallel and hybrid MPI+OpenMP jobs. Assuming |
| 153 | +that you installed the MPICH version of the MRChem MPI-parallel Conda package, |
| 154 | +the basic ``task_config`` argument to ``compute`` would look like: |
| 155 | + |
| 156 | +.. code-block:: python |
| 157 | +
|
| 158 | + task = { |
| 159 | + "nnodes": 1, # number of nodes |
| 160 | + "ncores": 12, # number of cores per task on each node |
| 161 | + "cores_per_rank": 6, # number of cores per MPI rank |
| 162 | + "use_mpiexec": True, # launch with MPI |
| 163 | + "mpiexec_command": "mpiexec -n {total_ranks}", # the invocation of MPI |
| 164 | + } |
| 165 | +
|
| 166 | +This task configuration will launch a MPI job with 2 ranks on a single node. |
| 167 | +Each rank has access to 6 cores for OpenMP parallelization. The ``provenance`` |
| 168 | +dictionary now shows: |
| 169 | + |
| 170 | +.. code-block:: bash |
| 171 | +
|
| 172 | + { |
| 173 | + "creator": "MRChem", |
| 174 | + "mpi_processes": 2, |
| 175 | + "routine": "mpiexec -n 2 /home/roberto/miniconda3/envs/mrchem-qcng/bin/mrchem.x", |
| 176 | + "total_cores": 12, |
| 177 | + "version": "1.1.0", |
| 178 | + "ncores": 12, |
| 179 | + "nnodes": 1, |
| 180 | + "ranks_per_node": 2, |
| 181 | + "cores_per_rank": 6, |
| 182 | + "total_ranks": 2 |
| 183 | + } |
| 184 | +
|
| 185 | +
|
| 186 | +The ``mpiexec_command`` is a string that will be interpolated to provide the |
| 187 | +exact invocation. In the above example, MRChem will be run with: |
| 188 | + |
| 189 | +.. code-block:: bash |
| 190 | +
|
| 191 | + mpiexec -n 2 /home/roberto/miniconda3/envs/mrchem-qcng/bin/mrchem.x |
| 192 | +
|
| 193 | +The following interpolation parameters are understood by QCEngine when creating |
| 194 | +the MPI invocation: |
| 195 | + |
| 196 | +- ``{nnodes}``: number of nodes. |
| 197 | +- ``{cores_per_rank}``: number of cores to use for each MPI rank. |
| 198 | +- ``{ranks_per_node}``: number of MPI ranks per node. Computed as ``ncores // cores_per_rank``. |
| 199 | +- ``{total_ranks}``: total number of MPI ranks. Computed as ``nnodes * ranks_per_node``. |
| 200 | + |
| 201 | +More complex MPI invocations are possible by setting the appropriate |
| 202 | +``mpiexec_command`` in the task configuration. For usage with a scheduler, such |
| 203 | +as SLURM, you should refer to the documentation of your computing cluster and |
| 204 | +the documentation of QCEngine. |
| 205 | + |
| 206 | + |
| 207 | +Geometry optimizations |
| 208 | +~~~~~~~~~~~~~~~~~~~~~~ |
| 209 | + |
| 210 | +Running geometry optimizations is just as easy as single compute. The following |
| 211 | +example optimizes the structure of water using the SVWN5 functional with MW4. |
| 212 | +The `geomeTRIC <https://geometric.readthedocs.io/en/latest/>`_ package is used |
| 213 | +as optimization driver, but `pyberny |
| 214 | +<https://jhrmnn.github.io/pyberny/algorithm.html>`_ or `optking |
| 215 | +<https://optking.readthedocs.io/en/latest/?badge=latest>`_ would also work. |
| 216 | + |
| 217 | +.. warning:: |
| 218 | + |
| 219 | + The computation of the molecular gradient can be affected by significant |
| 220 | + numerical noise for MW3 and MW4, to the point that it can be impossible to |
| 221 | + converge a geometry optimization. Using a tighter precision might help, but |
| 222 | + the cost of the calculation might be prohibitively large. |
| 223 | + |
| 224 | +.. code-block:: python |
| 225 | +
|
| 226 | + import qcelemental as qcel |
| 227 | + import qcengine as qcng |
| 228 | +
|
| 229 | + mol = qcel.models.Molecule( |
| 230 | + geometry=[ |
| 231 | + [ 0.29127930, 3.00875625, 0.20308515], |
| 232 | + [-1.21253048, 1.95820900, 0.10303324], |
| 233 | + [ 0.10002049, 4.24958115,-1.10222079] |
| 234 | + ], |
| 235 | + symbols=["O", "H", "H"], |
| 236 | + fix_com=True, |
| 237 | + fix_orientation=True, |
| 238 | + fix_symmetry="c1") |
| 239 | +
|
| 240 | + opt_input = { |
| 241 | + "keywords": { |
| 242 | + "program": "mrchem", |
| 243 | + "maxiter": 70 |
| 244 | + }, |
| 245 | + "input_specification": { |
| 246 | + "driver": "gradient", |
| 247 | + "model": { |
| 248 | + "method": "SVWN5", |
| 249 | + }, |
| 250 | + "keywords": { |
| 251 | + "world_prec": 1.0e-4, |
| 252 | + "SCF": { |
| 253 | + "guess_type": "core_dz", |
| 254 | + } |
| 255 | + } |
| 256 | + }, |
| 257 | + "initial_molecule": mol, |
| 258 | + } |
| 259 | +
|
| 260 | + opt = qcng.compute_procedure( |
| 261 | + opt_input, |
| 262 | + "geometric", |
| 263 | + task_config={"ncores": 20}) |
| 264 | +
|
| 265 | + print(opt.stdout) |
| 266 | +
|
| 267 | + print("==> Optimized geometry <==") |
| 268 | + print(opt.final_molecule.pretty_print()) |
| 269 | +
|
| 270 | + print("==> Optimized geometric parameters <==") |
| 271 | + for m in [[0, 1], [0, 2], [1, 0, 2]]: |
| 272 | + opt_val = opt.final_molecule.measure(m) |
| 273 | + print(f"Internal degree of freedom {m} = {opt_val:.3f}") |
| 274 | +
|
| 275 | +Running this script will print all the steps taken during the structural optimization. |
| 276 | +The final printout contains the optimized geometry: |
| 277 | + |
| 278 | +.. code-block:: bash |
| 279 | +
|
| 280 | + Geometry (in Angstrom), charge = 0.0, multiplicity = 1: |
| 281 | +
|
| 282 | + Center X Y Z |
| 283 | + ------------ ----------------- ----------------- ----------------- |
| 284 | + O -4.146209038013 2.134923126314 -3.559202294678 |
| 285 | + H -4.906566693905 1.536801624016 -3.587431156799 |
| 286 | + H -4.270830051398 2.773072094238 -4.275607223691 |
| 287 | +
|
| 288 | +and the optimized values of bond distances and bond angle: |
| 289 | + |
| 290 | +.. code-block:: bash |
| 291 | +
|
| 292 | + Internal degree of freedom [0, 1] = 1.829 |
| 293 | + Internal degree of freedom [0, 2] = 1.828 |
| 294 | + Internal degree of freedom [1, 0, 2] = 106.549 |
0 commit comments