Correlation Plots for Computed Results
**************************************

Analysis Plotting Benchmark

This example illustrates how to plot, for example, forces predicted by
one engine against forces predicted by another.

**Requires:** AMS2026 or later

-[ Related tutorials ]-

* Generating a batch of jobs

* Python Getting Started tutorial

-[ Related documentation ]-

* PLAMS plotting tools

* PLAMS JobAnalysis

[image]


Initial imports
===============

   import scm.plams as plams
   from scm.plams.tools.plot import plot_correlation, get_correlation_xy
   import matplotlib.pyplot as plt


Define two engines to compare
=============================

Here we choose GFNFF and GFN1-xTB

   e1 = plams.Settings()
   e1.input.GFNFF

   e2 = plams.Settings()
   e2.input.DFTB.Model = "GFN1-xTB"

Let’s use a glycine molecule generated from SMILES:

   from scm.base import ChemicalSystem

   glycine = ChemicalSystem.from_smiles("C(C(=O)O)N")
   plams.view(glycine, width=200, height=200, padding=-0.5)

   [image: image generated from notebook][image]

Run a single-point calculation storing the Gradients (negative
forces):

   sp = plams.Settings()
   sp.input.ams.Task = "SinglePoint"
   sp.input.ams.Properties.Gradients = "Yes"
   sp.runscript.nproc = 1  # run in serial

   job1 = plams.AMSJob(settings=sp + e1, name="glycine-engine1", molecule=glycine)
   job2 = plams.AMSJob(settings=sp + e2, name="glycine-engine2", molecule=glycine)

   job1.run()
   job2.run()

   [18.03|15:53:50] JOB glycine-engine1 STARTED
   [18.03|15:53:50] JOB glycine-engine1 RUNNING
   [18.03|15:53:50] JOB glycine-engine1 FINISHED
   [18.03|15:53:50] JOB glycine-engine1 SUCCESSFUL
   [18.03|15:53:50] JOB glycine-engine2 STARTED
   [18.03|15:53:50] JOB glycine-engine2 RUNNING
   [18.03|15:53:50] JOB glycine-engine2 FINISHED
   [18.03|15:53:50] JOB glycine-engine2 SUCCESSFUL





   <scm.plams.interfaces.adfsuite.ams.AMSResults at 0x174a7ae20>

   plot_correlation(
       job1,
       job2,
       section="AMSResults",
       variable="Gradients",
       file="engine",
   );

   [image: image generated from notebook][image]

To get the actual numbers, use "get_correlation_xy":

   x, y = plams.tools.plot.get_correlation_xy(
       job1, job2, section="AMSResults", variable="Gradients", file="engine"
   )
   print("x")
   print(x)
   print("y")
   print(y)

   x
   [ 1.19424632e-03  1.92134706e-04  6.66726815e-04 -2.68722244e-02
    -9.34887191e-04 -4.01123910e-02  1.17988215e-02  4.31998219e-03
     5.43487308e-02  8.33720932e-03  5.26401890e-04  1.24250877e-02
     8.91893628e-03  8.12876763e-03  2.87360972e-02 -4.14954879e-05
    -3.23202276e-03 -9.44698138e-03 -3.49925920e-03  8.54464635e-03
     1.69353194e-03  1.06479919e-02 -5.23262379e-03 -2.43272087e-02
     7.04191867e-03 -1.95558968e-02 -1.18022032e-02 -1.75261448e-02
     7.24349781e-03 -1.21813902e-02]
   y
   [ 0.00326195 -0.00599574  0.02250686 -0.03506961 -0.00156447 -0.07804685
    -0.0008342   0.01094663  0.08968399  0.02986723 -0.00458719  0.01262981
     0.01335181  0.01483532  0.01659624 -0.00041231 -0.0029202  -0.00449173
    -0.00404233  0.01096039 -0.00074098  0.00750758 -0.006259   -0.03825885
     0.00589418 -0.01943988 -0.00943605 -0.01952429  0.00402414 -0.01044246]


Compare multiple jobs
=====================

   smiles_list = ["CC=C", "CCCO", "C(C(=O)O)N"]
   names = ["propene", "propanol", "glycine"]
   molecules = [ChemicalSystem.from_smiles(x) for x in smiles_list]
   imgs = {n: plams.view(mol, width=200, height=200) for n, mol in zip(names, molecules)}
   plams.plot_image_grid(imgs, rows=1);

   [image: image generated from notebook][image]

   jobs1 = [
       plams.AMSJob(settings=sp + e1, name="e1" + name, molecule=mol)
       for name, mol in zip(names, molecules)
   ]
   jobs2 = [
       plams.AMSJob(settings=sp + e2, name="e2" + name, molecule=mol)
       for name, mol in zip(names, molecules)
   ]

   for job in jobs1 + jobs2:
       job.run()

   [18.03|15:53:53] JOB e1propene STARTED
   [18.03|15:53:53] JOB e1propene RUNNING
   [18.03|15:53:53] JOB e1propene FINISHED
   [18.03|15:53:53] JOB e1propene SUCCESSFUL
   [18.03|15:53:53] JOB e1propanol STARTED
   [18.03|15:53:53] JOB e1propanol RUNNING
   [18.03|15:53:53] JOB e1propanol FINISHED
   [18.03|15:53:53] JOB e1propanol SUCCESSFUL
   [18.03|15:53:53] JOB e1glycine STARTED
   [18.03|15:53:53] Job e1glycine previously run as glycine-engine1, using old results
   ... output trimmed ....
   [18.03|15:53:53] JOB e2propene RUNNING
   [18.03|15:53:54] JOB e2propene FINISHED
   [18.03|15:53:54] JOB e2propene SUCCESSFUL
   [18.03|15:53:54] JOB e2propanol STARTED
   [18.03|15:53:54] JOB e2propanol RUNNING
   [18.03|15:53:54] JOB e2propanol FINISHED
   [18.03|15:53:54] JOB e2propanol SUCCESSFUL
   [18.03|15:53:54] JOB e2glycine STARTED
   [18.03|15:53:54] Job e2glycine previously run as glycine-engine2, using old results
   [18.03|15:53:54] JOB e2glycine COPIED

The correlation plot can be plotted as before. You can also add a unit
conversion to get your preferred units, and add custom xlabel and
ylabel:

   unit = "eV/angstrom"
   multiplier = plams.Units.convert(1.0, "hartree/bohr", unit)

   plot_correlation(
       jobs1,
       jobs2,
       section="AMSResults",
       variable="Gradients",
       file="engine",
       xlabel="Engine 1",
       ylabel="Engine 2",
       unit=unit,
       multiplier=multiplier,
   );

   [image: image generated from notebook][image]

   plot_correlation(
       jobs1,
       jobs2,
       section="AMSResults",
       variable="Charges",
       file="engine",
       xlabel="Engine 1",
       ylabel="Engine 2",
   );

   [image: image generated from notebook][image]


Use Task Replay to compare multiple frames from a trajectory
============================================================

The forces from an MD job can be stored with
"writeenginegradients=True"

   md = plams.AMSNVEJob(
       settings=e1,
       name="nve-md-e1",
       molecule=glycine,
       velocities=400,
       nsteps=100,
       samplingfreq=10,
       writeenginegradients=True,
   )
   md.run();

   [18.03|15:53:54] JOB nve-md-e1 STARTED
   [18.03|15:53:54] JOB nve-md-e1 RUNNING
   [18.03|15:53:55] JOB nve-md-e1 FINISHED
   [18.03|15:53:55] JOB nve-md-e1 SUCCESSFUL

When using the Replay task, set "Properties.Gradients" to get the
forces:

   replay_s = plams.Settings()
   replay_s.input.ams.Task = "Replay"
   replay_s.input.ams.Properties.Gradients = "Yes"
   replay_s.input.ams.Replay.File = md.results.rkfpath()
   replay = plams.AMSJob(settings=e2 + replay_s, name="replay-e2")
   replay.run();

   [18.03|15:53:55] JOB replay-e2 STARTED
   [18.03|15:53:55] JOB replay-e2 RUNNING
   [18.03|15:53:56] JOB replay-e2 FINISHED
   [18.03|15:53:56] JOB replay-e2 SUCCESSFUL

For the MD job the gradients (negative forces) are stored in
"History%EngineGradients", but for the Replay job they are stored in
"History%Gradients". Use the "alt_variable" to specify the variable
for the second job:

   plot_correlation(
       md, replay, section="History", variable="EngineGradients", alt_variable="Gradients", file="ams"
   );

   [image: image generated from notebook][image]


See also
========


Related tutorials
-----------------

* Generating a batch of jobs

* Python Getting Started tutorial


Related documentation
---------------------

* PLAMS plotting tools

* PLAMS JobAnalysis


Python Script
=============

   #!/usr/bin/env python
   # coding: utf-8

   # ## Initial imports

   import scm.plams as plams
   from scm.plams.tools.plot import plot_correlation, get_correlation_xy
   import matplotlib.pyplot as plt


   # ## Define two engines to compare
   #
   # Here we choose GFNFF and GFN1-xTB

   e1 = plams.Settings()
   e1.input.GFNFF

   e2 = plams.Settings()
   e2.input.DFTB.Model = "GFN1-xTB"


   # Let's use a glycine molecule generated from SMILES:

   from scm.base import ChemicalSystem

   glycine = ChemicalSystem.from_smiles("C(C(=O)O)N")
   plams.view(glycine, width=200, height=200, padding=-0.5)


   # Run a single-point calculation storing the Gradients (negative forces):

   sp = plams.Settings()
   sp.input.ams.Task = "SinglePoint"
   sp.input.ams.Properties.Gradients = "Yes"
   sp.runscript.nproc = 1  # run in serial


   job1 = plams.AMSJob(settings=sp + e1, name="glycine-engine1", molecule=glycine)
   job2 = plams.AMSJob(settings=sp + e2, name="glycine-engine2", molecule=glycine)


   job1.run()
   job2.run()


   plot_correlation(
       job1,
       job2,
       section="AMSResults",
       variable="Gradients",
       file="engine",
   )


   # To get the actual numbers, use ``get_correlation_xy``:

   x, y = plams.tools.plot.get_correlation_xy(job1, job2, section="AMSResults", variable="Gradients", file="engine")
   print("x")
   print(x)
   print("y")
   print(y)


   # ## Compare multiple jobs

   smiles_list = ["CC=C", "CCCO", "C(C(=O)O)N"]
   names = ["propene", "propanol", "glycine"]
   molecules = [ChemicalSystem.from_smiles(x) for x in smiles_list]
   imgs = {n: plams.view(mol, width=200, height=200) for n, mol in zip(names, molecules)}
   plams.plot_image_grid(imgs, rows=1)


   jobs1 = [plams.AMSJob(settings=sp + e1, name="e1" + name, molecule=mol) for name, mol in zip(names, molecules)]
   jobs2 = [plams.AMSJob(settings=sp + e2, name="e2" + name, molecule=mol) for name, mol in zip(names, molecules)]


   for job in jobs1 + jobs2:
       job.run()


   # The correlation plot can be plotted as before. You can also add a unit conversion to get your preferred units, and add custom xlabel and ylabel:

   unit = "eV/angstrom"
   multiplier = plams.Units.convert(1.0, "hartree/bohr", unit)

   plot_correlation(
       jobs1,
       jobs2,
       section="AMSResults",
       variable="Gradients",
       file="engine",
       xlabel="Engine 1",
       ylabel="Engine 2",
       unit=unit,
       multiplier=multiplier,
   )


   plot_correlation(
       jobs1,
       jobs2,
       section="AMSResults",
       variable="Charges",
       file="engine",
       xlabel="Engine 1",
       ylabel="Engine 2",
   )


   # ## Use Task Replay to compare multiple frames from a trajectory

   # The forces from an MD job can be stored with ``writeenginegradients=True``

   md = plams.AMSNVEJob(
       settings=e1,
       name="nve-md-e1",
       molecule=glycine,
       velocities=400,
       nsteps=100,
       samplingfreq=10,
       writeenginegradients=True,
   )
   md.run()


   # When using the Replay task, set ``Properties.Gradients`` to get the forces:

   replay_s = plams.Settings()
   replay_s.input.ams.Task = "Replay"
   replay_s.input.ams.Properties.Gradients = "Yes"
   replay_s.input.ams.Replay.File = md.results.rkfpath()
   replay = plams.AMSJob(settings=e2 + replay_s, name="replay-e2")
   replay.run()


   # For the MD job the gradients (negative forces) are stored in ``History%EngineGradients``, but for the Replay job they are stored in ``History%Gradients``. Use the ``alt_variable`` to specify the variable for the second job:

   plot_correlation(md, replay, section="History", variable="EngineGradients", alt_variable="Gradients", file="ams")
