Active Learning: Single-Molecule Production Simulation
******************************************************

MLPotential Active Learning ParAMS M3GNet Geometry Optimization

Reuse the retrained active-learning model in a production calculation,
then inspect optimized structures and vibrational properties for the
single-molecule test case.

**Requires:** AMS2026 or later

**Required AMS Packages:** m3gnet

-[ Related examples ]-

* Active Learning: Single-Molecule Setup and Run

* Active Learning: Single-Molecule Results

* Active Learning: Single-Molecule Comparison to M3GNet-UP-2022

* M3GNet Custom Model

-[ Related tutorials ]-

* Simple Active Learning tutorial

-[ Related documentation ]-

* Simple Active Learning Python API

* Geometry optimization

* ParAMSJob API

* ParAMSResults API

* ResultsImporter overview

[image]


Initial imports and load active learning job from disk
======================================================

   from scm.simple_active_learning import SimpleActiveLearningJob
   import scm.plams as plams
   import os
   import matplotlib.pyplot as plt

   plams.init()

   PLAMS working folder: /path/to/plams_workdir.003

   # replace the path with your own path !
   previous_sal_job_path = os.path.expandvars(
       "$AMSHOME/examples/SAL/Output/SingleMolecule/plams_workdir/sal"
   )
   sal_job = SimpleActiveLearningJob.load_external(previous_sal_job_path)
   params_job = sal_job.results.get_params_job()


Structure for production job
============================

We could initialize a PLAMS molecule in many different ways. Here, we
get the final frame from the final production simulation in the SAL
job, and preoptimize it with UFF.

   molecule = sal_job.results.get_main_molecule()
   molecule = plams.preoptimize(molecule, model="UFF")
   plams.plot_molecule(molecule)

   [image: image generated from notebook][image]


Settings for production job
===========================

Let’s run a geometry optimization + frequencies. But you could run any
AMS job!

The "MaxRestarts" option is useful when calculating normal modes. If
the geometry optimizer converges to a transition state, it will
continue until it finds a local minimum!

   s = plams.Settings()
   s.input.ams.Task = "GeometryOptimization"
   s.input.ams.Properties.NormalModes = "Yes"
   s.input.ams.GeometryOptimization.MaxRestarts = 5
   s.runscript.nproc = 1

   retrained_engine_settings = params_job.results.get_production_engine_settings()

   new_job = plams.AMSJob(
       settings=s + retrained_engine_settings, name="retrained_m3gnet", molecule=molecule
   )
   print(new_job.get_input())

   GeometryOptimization
     MaxRestarts 5
   End

   Properties
   ... output trimmed ....
     Model Custom
     ParameterDir /path/to/plams_workdir/sal/step4_attempt1_training/results/optimization/m3gnet/m3gnet
   EndEngine

   new_job.run();

   [02.02|12:39:32] JOB retrained_m3gnet STARTED
   [02.02|12:39:32] JOB retrained_m3gnet RUNNING
   [02.02|12:39:47] JOB retrained_m3gnet FINISHED
   [02.02|12:39:47] JOB retrained_m3gnet SUCCESSFUL


Optimized structure
-------------------

   plams.plot_molecule(new_job.results.get_main_molecule())

   [image: image generated from notebook][image]


Frequencies
-----------

   width = 50  # cm^-1
   x, y = new_job.results.get_ir_spectrum(
       broadening_width=width, post_process="all_intensities_to_1"
   )
   plt.plot(x, y, label="Retrained M3GNet")
   plt.xlabel("Frequency (cm^-1)")
   plt.ylabel("Normal mode count")
   plt.legend();

   [image: image generated from notebook][image]


Compare to reference UFF
========================

In this case, we can also calculate the normal modes with the
reference method (UFF) and compare:

   ref_engine_settings = plams.Settings()
   ref_engine_settings.input.ForceField.Type = "UFF"
   ref_job = plams.AMSJob(
       settings=s + ref_engine_settings, name="uff_ref", molecule=molecule
   )
   ref_job.run();

   [02.02|12:39:48] JOB uff_ref STARTED
   [02.02|12:39:48] JOB uff_ref RUNNING
   [02.02|12:39:48] JOB uff_ref FINISHED
   [02.02|12:39:48] JOB uff_ref SUCCESSFUL

   retrained_structure_rmsd = plams.Molecule.rmsd(
       ref_job.results.get_main_molecule(),
       new_job.results.get_main_molecule(),
       ignore_hydrogen=True,
   )
   print(f"Structural RMSD: {retrained_structure_rmsd:.2f} angstrom")

   Structural RMSD: 0.10 angstrom

   x_ref, y_ref = ref_job.results.get_ir_spectrum(
       broadening_width=width, post_process="all_intensities_to_1"
   )
   plt.plot(x, y, label="Retrained M3GNet")
   plt.plot(x_ref, y_ref, label="UFF (reference method)")
   plt.xlabel("Frequency (cm^-1)")
   plt.ylabel("Normal mode count")
   plt.legend();

   [image: image generated from notebook][image]

The agreement looks very good! The only significant difference is the
highest-frequency vibration (the O-H streching vibration). This
frequency is very sensitive to the calculated forces near the
equilibrium (minimum) structure. The agreement could have been
improved by

* having more training data, for example by setting a tighter success
  criterion for the forces and energy in the active learning

* running the active learning MD at a lower temperature (closer to the
  equilibrium structure, but this would mean less conformational
  sampling)

* training for more epochs

Tip: check if the vibrational frequencies with retrained M3GNet or
M3GNet-UP-2022 agree better or worse with the reference UFF
calculation.


See also
========


Related examples
----------------

* Active Learning: Single-Molecule Setup and Run

* Active Learning: Single-Molecule Results

* Active Learning: Single-Molecule Comparison to M3GNet-UP-2022

* M3GNet Custom Model


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

* Simple Active Learning tutorial


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

* Simple Active Learning Python API

* Geometry optimization

* ParAMSJob API

* ParAMSResults API

* ResultsImporter overview


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

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

   # ## Initial imports and load active learning job from disk

   from scm.simple_active_learning import SimpleActiveLearningJob
   import scm.plams as plams
   import os
   import matplotlib.pyplot as plt


   plams.init()


   # replace the path with your own path !
   previous_sal_job_path = os.path.expandvars("$AMSHOME/examples/SAL/Output/SingleMolecule/plams_workdir/sal")
   sal_job = SimpleActiveLearningJob.load_external(previous_sal_job_path)
   params_job = sal_job.results.get_params_job()


   # ## Structure for production job
   #
   # We could initialize a PLAMS molecule in many different ways. Here, we get the final frame from the final production simulation in the SAL job, and preoptimize it with UFF.

   molecule = sal_job.results.get_main_molecule()
   molecule = plams.preoptimize(molecule, model="UFF")
   plams.plot_molecule(molecule)


   # ## Settings for production job
   #
   # Let's run a geometry optimization + frequencies. But you could run any AMS job!
   #
   # The ``MaxRestarts`` option is useful when calculating normal modes. If the geometry optimizer converges to a transition state, it will continue until it finds a local minimum!

   s = plams.Settings()
   s.input.ams.Task = "GeometryOptimization"
   s.input.ams.Properties.NormalModes = "Yes"
   s.input.ams.GeometryOptimization.MaxRestarts = 5
   s.runscript.nproc = 1

   retrained_engine_settings = params_job.results.get_production_engine_settings()

   new_job = plams.AMSJob(settings=s + retrained_engine_settings, name="retrained_m3gnet", molecule=molecule)
   print(new_job.get_input())


   new_job.run()
   # ### Optimized structure

   plams.plot_molecule(new_job.results.get_main_molecule())


   # ### Frequencies

   width = 50  # cm^-1
   x, y = new_job.results.get_ir_spectrum(broadening_width=width, post_process="all_intensities_to_1")
   plt.plot(x, y, label="Retrained M3GNet")
   plt.xlabel("Frequency (cm^-1)")
   plt.ylabel("Normal mode count")
   plt.legend()
   plt.gcf().savefig("picture1.png")
   plt.clf()


   # ## Compare to reference UFF
   # In this case, we can also calculate the normal modes with the reference method (UFF) and compare:

   ref_engine_settings = plams.Settings()
   ref_engine_settings.input.ForceField.Type = "UFF"
   ref_job = plams.AMSJob(settings=s + ref_engine_settings, name="uff_ref", molecule=molecule)
   ref_job.run()
   retrained_structure_rmsd = plams.Molecule.rmsd(
       ref_job.results.get_main_molecule(),
       new_job.results.get_main_molecule(),
       ignore_hydrogen=True,
   )
   print(f"Structural RMSD: {retrained_structure_rmsd:.2f} angstrom")


   x_ref, y_ref = ref_job.results.get_ir_spectrum(broadening_width=width, post_process="all_intensities_to_1")
   plt.plot(x, y, label="Retrained M3GNet")
   plt.plot(x_ref, y_ref, label="UFF (reference method)")
   plt.xlabel("Frequency (cm^-1)")
   plt.ylabel("Normal mode count")
   plt.legend()
   plt.gcf().savefig("picture2.png")
   plt.clf()


   # The agreement looks very good! The only significant difference is the highest-frequency vibration (the O-H streching vibration). This frequency is very sensitive to the calculated forces near the equilibrium (minimum) structure. The agreement could have been improved by
   #
   # * having more training data, for example by setting a tighter success criterion for the forces and energy in the active learning
   # * running the active learning MD at a lower temperature (closer to the equilibrium structure, but this would mean less conformational sampling)
   # * training for more epochs
   #
   # Tip: check if the vibrational frequencies with retrained M3GNet or M3GNet-UP-2022 agree better or worse with the reference UFF calculation.
