Conformer Generation, Rescoring, and Filtering
**********************************************

GFNFF DFTB Conformers

Generate conformers, reoptimize them at higher levels of theory, and
filter the final ensemble to keep only the most relevant structures.

**Requires:** AMS2026 or later

-[ Related examples ]-

* Conformer Generation for Multiple Molecules

* Automated Transition-State Workflow for a Diels-Alder Addition

* Screening Molecules for Promising Electronic Excitations

-[ Related tutorials ]-

* Conformers tutorial

* Multiple molecules conformers tutorial

-[ Related documentation ]-

* AMS conformers utility

* GFNFF overview

* DFTB tasks and properties

[image]


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

   import os
   import sys

   import matplotlib.pyplot as plt
   import numpy as np

   from scm.conformers import ConformersJob
   from scm.plams import Settings, view


Initial structure
=================

   from scm.base import ChemicalSystem

   molecule = ChemicalSystem.from_smiles("OC(CC1c2ccccc2Sc2ccccc21)CN1CCCC1")
   view(molecule, width=300, height=300, direction="along_pca3")

   [image: image generated from notebook][image]


Generate conformers with RDKit and UFF
======================================

The fastest way to generate conformers is to use RDKit with the UFF
force field.

Below we specify to generate 16 initial conformers. The final number
of conformers may be smaller, as the geometry optimization may cause
several structures to enter the same minimum.


Conformer generation settings
-----------------------------

   s = Settings()
   s.input.ams.Task = "Generate"  # default
   s.input.ams.Generator.Method = "RDKit"  # default
   s.input.ams.Generator.RDKit.InitialNConformers = 16  # optional, non-default
   s.input.ForceField.Type = "UFF"  # default


Conformer generation input file
-------------------------------

   print(ConformersJob(settings=s).get_input())

   Generator
     Method RDKit
     RDKit
       InitialNConformers 16
     End
   End

   Task Generate


   Engine ForceField
     Type UFF
   EndEngine


Run conformer generation
------------------------

   generate_job = ConformersJob(name="generate", molecule=molecule, settings=s)
   generate_job.run();

   [09.04|12:59:17] JOB generate STARTED
   [09.04|12:59:17] JOB generate RUNNING
   [09.04|12:59:27] JOB generate FINISHED
   [09.04|12:59:27] JOB generate SUCCESSFUL


Conformer generation results
============================


Some helper functions
---------------------

   import pandas as pd


   def get_df(job: ConformersJob, temperature=298, unit="kcal/mol") -> pd.DataFrame:
       relative_energies = job.results.get_relative_energies(unit)
       pop = job.results.get_boltzmann_distribution(temperature)
       df = pd.DataFrame(
           {
               "Conformer Id": np.arange(1, len(relative_energies) + 1),
               f"Rel. E ({unit})": relative_energies,
               f"Population at {temperature} K": pop,
           }
       )
       return df


Actual results
--------------

Below we see that the **conformer generation gave 15 distinct
conformers**, where the highest-energy conformer is 18 kcal/mol higher
in energy than the lowest energy conformer.

You can also see the **relative populations** of these conformers at
the specified temperature. The populations are calculated from the
**Boltzmann distribution** and the relative energies.

   unit = "kcal/mol"
   temperature = 298

   df = get_df(generate_job, temperature, unit)
   df

Conformer Id

Rel. E (kcal/mol)

Population at 298 K

0

1

0.000000

3.631715e-01

1

2

0.006054

3.594777e-01

2

3

0.472937

1.634079e-01

3

4

0.867098

8.398580e-02

4

5

1.764069

1.846683e-02

5

6

2.111678

1.026751e-02

6

7

3.675484

7.321730e-04

7

8

3.957073

4.550966e-04

8

9

5.730983

2.275975e-05

9

10

6.072706

1.278078e-05

10

11

15.374290

1.927672e-12

11

12

17.740336

3.546829e-14

12

13

18.151045

1.772709e-14

13

14

19.330941

2.417323e-15

14

15

24.757089

2.534824e-19

   from scm.conformers.plams.plot import plot_conformers

   axes = plot_conformers(
       generate_job,
       4,
       temperature=temperature,
       unit=unit,
       lowest=True,
       width=300,
       height=300,
       direction="along_pca3",
   )
   axes;

   [image: image generated from notebook][image]


Re-optimize conformers with GFNFF
=================================

The UFF force field is not very accurate for geometries and energies.
From an initial conformer set you can reoptimize it with a better
level of theory.

The **Optimize** task performs **GeometryOptimization** jobs on each
conformer in a set.

Below, the most stable conformers (within 8 kcal/mol of the most
stable conformer) at the UFF level of theory are re-optimized with
GFNFF, which gives more accurate geometries.

   s = Settings()
   s.input.ams.Task = "Optimize"
   # must be absolute path
   s.input.ams.InputConformersSet = os.path.abspath(generate_job.results.rkfpath())
   # only conformers within 8 kcal/mol at the PREVIOUS level of theory
   s.input.ams.InputMaxEnergy = 8.0
   s.input.GFNFF  # or choose a different engine if you don't have a GFNFF license

   reoptimize_job = ConformersJob(settings=s, name="reoptimize")
   print(reoptimize_job.get_input())

   InputConformersSet /Users/ormrodmorley/Documents/code/ams/amshome/userdoc/PythonExamples/conformers-generation/plams_workdir.005/generate/conformers.rkf

   InputMaxEnergy 8.0

   Task Optimize


   Engine GFNFF
   EndEngine

   reoptimize_job.run();

   [09.04|12:59:32] JOB reoptimize STARTED
   [09.04|12:59:32] JOB reoptimize RUNNING
   [09.04|12:59:34] JOB reoptimize FINISHED
   [09.04|12:59:34] JOB reoptimize SUCCESSFUL

   df = get_df(reoptimize_job, temperature=temperature, unit=unit)
   df

Conformer Id

Rel. E (kcal/mol)

Population at 298 K

0

1

0.000000

0.455124

1

2

0.432698

0.219180

2

3

0.861018

0.106337

3

4

0.901472

0.099315

4

5

1.311437

0.049700

5

6

1.469166

0.038079

6

7

1.867458

0.019435

7

8

2.476361

0.006951

8

9

2.693401

0.004818

9

10

3.589730

0.001061

   axes = plot_conformers(reoptimize_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes;

   [image: image generated from notebook][image]


Score conformers with DFTB
==========================

If you have many conformers or a very large molecule, it can be
computationally expensive to do the conformer generation or
reoptimization and a high level of theory.

The **Score** task runs **SinglePoint** jobs on the conformers in a
set. This lets you use a more computationally expensive method. Here,
we choose DFTB, although normally you may choose some DFT method.

   s = Settings()
   s.input.ams.Task = "Score"
   # must be absolute path
   s.input.ams.InputConformersSet = os.path.abspath(reoptimize_job.results.rkfpath())
   # only conformers within 4 kcal/mol at the PREVIOUS level of theory
   s.input.ams.InputMaxEnergy = 4.0
   s.input.DFTB.Model = "GFN1-xTB"  # or choose a different engine if you don't have a DFTB license
   # s.input.adf.XC.GGA = 'PBE'                       # to use ADF PBE
   # s.input.adf.XC.DISPERSION = 'GRIMME3 BJDAMP'     # to use ADF PBE with Grimme D3(BJ) dispersion

   score_job = ConformersJob(settings=s, name="score")
   score_job.run();

   [09.04|12:59:39] JOB score STARTED
   [09.04|12:59:39] JOB score RUNNING
   [09.04|12:59:41] JOB score FINISHED
   [09.04|12:59:41] JOB score SUCCESSFUL

   df = get_df(score_job, temperature=temperature, unit=unit)
   df

Conformer Id

Rel. E (kcal/mol)

Population at 298 K

0

1

0.000000

0.674970

1

2

0.791452

0.177360

2

3

1.415553

0.061824

3

4

1.852765

0.029547

4

5

2.001944

0.022968

5

6

2.241456

0.015327

6

7

2.465533

0.010499

7

8

2.818300

0.005787

8

9

3.800985

0.001101

9

10

4.142879

0.000618

   axes = plot_conformers(score_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes;

   [image: image generated from notebook][image]

Here, you see that from the conformers in the set, **DFTB predicts a
different lowest-energy conformer than GFNFF** (compare to previous
figure).


Filter a conformer set
======================

In practice, you may have generated thousands of conformers for a
particular structure. Many of those conformers may be so high in
energy that their Boltzmann weights are very small.

The **Filter** task only filters the conformers, it does not perform
any additional calculations. It can be used to reduce a conformer set
so that it is more convenient to work with.

Below, we filter the conformers set to only the conformers within 2
kcal/mol of the minimum.

   s = Settings()
   s.input.ams.Task = "Filter"
   s.input.ams.InputConformersSet = os.path.abspath(score_job.results.rkfpath())
   s.input.ams.InputMaxEnergy = 2.0

   filter_job = ConformersJob(settings=s, name="filter")
   filter_job.run();

   [09.04|12:59:47] JOB filter STARTED
   [09.04|12:59:47] JOB filter RUNNING
   [09.04|12:59:48] JOB filter FINISHED
   [09.04|12:59:48] JOB filter SUCCESSFUL

   df = get_df(filter_job, temperature=temperature, unit=unit)
   df

Conformer Id

Rel. E (kcal/mol)

Population at 298 K

0

1

0.000000

0.715237

1

2

0.791452

0.187941

2

3

1.415553

0.065512

3

4

1.852765

0.031310

   axes = plot_conformers(filter_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes;

   [image: image generated from notebook][image]

The structures and energies are identical to before. However, the
relative populations changed slightly as there are now fewer
conformers in the set.


More about conformers
=====================

* Try **CREST** instead of RDKit to generate the initial conformer set

* The **Expand** task can be used to expand a set of conformers.


See also
========


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

* Conformer Generation for Multiple Molecules

* Automated Transition-State Workflow for a Diels-Alder Addition

* Screening Molecules for Promising Electronic Excitations


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

* Conformers tutorial

* Multiple molecules conformers tutorial


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

* AMS conformers utility

* GFNFF overview

* DFTB tasks and properties


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

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

   # ## Initial imports

   import os
   import sys

   import matplotlib.pyplot as plt
   import numpy as np

   from scm.conformers import ConformersJob
   from scm.plams import Settings, view


   # ## Initial structure

   from scm.base import ChemicalSystem

   molecule = ChemicalSystem.from_smiles("OC(CC1c2ccccc2Sc2ccccc21)CN1CCCC1")
   view(molecule, width=300, height=300, direction="along_pca3", picture_path="picture1.png")


   # ## Generate conformers with RDKit and UFF
   # The fastest way to generate conformers is to use RDKit with the UFF force field.
   #
   # Below we specify to generate 16 initial conformers. The final number of conformers may be smaller, as the geometry optimization may cause several structures to enter the same minimum.

   # ### Conformer generation settings

   s = Settings()
   s.input.ams.Task = "Generate"  # default
   s.input.ams.Generator.Method = "RDKit"  # default
   s.input.ams.Generator.RDKit.InitialNConformers = 16  # optional, non-default
   s.input.ForceField.Type = "UFF"  # default


   # ### Conformer generation input file

   print(ConformersJob(settings=s).get_input())


   # ### Run conformer generation

   generate_job = ConformersJob(name="generate", molecule=molecule, settings=s)
   generate_job.run()
   # ## Conformer generation results

   # ### Some helper functions

   import pandas as pd


   def get_df(job: ConformersJob, temperature=298, unit="kcal/mol") -> pd.DataFrame:
       relative_energies = job.results.get_relative_energies(unit)
       pop = job.results.get_boltzmann_distribution(temperature)
       df = pd.DataFrame(
           {
               "Conformer Id": np.arange(1, len(relative_energies) + 1),
               f"Rel. E ({unit})": relative_energies,
               f"Population at {temperature} K": pop,
           }
       )
       return df


   # ### Actual results
   #
   # Below we see that the **conformer generation gave 15 distinct conformers**, where the highest-energy conformer is 18 kcal/mol higher in energy than the lowest energy conformer.
   #
   # You can also see the **relative populations** of these conformers at the specified temperature. The populations are calculated from the **Boltzmann distribution** and the relative energies.

   unit = "kcal/mol"
   temperature = 298


   df = get_df(generate_job, temperature, unit)
   print(df)


   from scm.conformers.plams.plot import plot_conformers

   axes = plot_conformers(
       generate_job,
       4,
       temperature=temperature,
       unit=unit,
       lowest=True,
       width=300,
       height=300,
       direction="along_pca3",
   )
   axes
   import matplotlib.pyplot as plt

   plt.gcf().savefig("picture2.png")


   # ## Re-optimize conformers with GFNFF
   #
   # The UFF force field is not very accurate for geometries and energies. From an initial conformer set you can reoptimize it with a better level of theory.
   #
   # The **Optimize** task performs **GeometryOptimization** jobs on each conformer in a set.
   #
   # Below, the most stable conformers (within 8 kcal/mol of the most stable conformer) at the UFF level of theory are re-optimized with GFNFF, which gives more accurate geometries.

   s = Settings()
   s.input.ams.Task = "Optimize"
   # must be absolute path
   s.input.ams.InputConformersSet = os.path.abspath(generate_job.results.rkfpath())
   # only conformers within 8 kcal/mol at the PREVIOUS level of theory
   s.input.ams.InputMaxEnergy = 8.0
   s.input.GFNFF  # or choose a different engine if you don't have a GFNFF license

   reoptimize_job = ConformersJob(settings=s, name="reoptimize")
   print(reoptimize_job.get_input())


   reoptimize_job.run()
   df = get_df(reoptimize_job, temperature=temperature, unit=unit)
   print(df)


   axes = plot_conformers(reoptimize_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes
   import matplotlib.pyplot as plt

   plt.gcf().savefig("picture3.png")


   # ## Score conformers with DFTB
   #
   # If you have many conformers or a very large molecule, it can be computationally expensive to do the conformer generation or reoptimization and a high level of theory.
   #
   # The **Score** task runs **SinglePoint** jobs on the conformers in a set. This lets you use a more computationally expensive method. Here, we choose DFTB, although normally you may choose some DFT method.

   s = Settings()
   s.input.ams.Task = "Score"
   # must be absolute path
   s.input.ams.InputConformersSet = os.path.abspath(reoptimize_job.results.rkfpath())
   # only conformers within 4 kcal/mol at the PREVIOUS level of theory
   s.input.ams.InputMaxEnergy = 4.0
   s.input.DFTB.Model = "GFN1-xTB"  # or choose a different engine if you don't have a DFTB license
   # s.input.adf.XC.GGA = 'PBE'                       # to use ADF PBE
   # s.input.adf.XC.DISPERSION = 'GRIMME3 BJDAMP'     # to use ADF PBE with Grimme D3(BJ) dispersion

   score_job = ConformersJob(settings=s, name="score")
   score_job.run()
   df = get_df(score_job, temperature=temperature, unit=unit)
   print(df)


   axes = plot_conformers(score_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes
   import matplotlib.pyplot as plt

   plt.gcf().savefig("picture4.png")


   # Here, you see that from the conformers in the set, **DFTB predicts a different lowest-energy conformer than GFNFF** (compare to previous figure).

   # ## Filter a conformer set
   #
   # In practice, you may have generated thousands of conformers for a particular structure. Many of those conformers may be so high in energy that their Boltzmann weights are very small.
   #
   # The **Filter** task only filters the conformers, it does not perform any additional calculations. It can be used to reduce a conformer set so that it is more convenient to work with.
   #
   # Below, we filter the conformers set to only the conformers within 2 kcal/mol of the minimum.

   s = Settings()
   s.input.ams.Task = "Filter"
   s.input.ams.InputConformersSet = os.path.abspath(score_job.results.rkfpath())
   s.input.ams.InputMaxEnergy = 2.0

   filter_job = ConformersJob(settings=s, name="filter")
   filter_job.run()
   df = get_df(filter_job, temperature=temperature, unit=unit)
   print(df)


   axes = plot_conformers(filter_job, 4, temperature=temperature, unit=unit, lowest=True)
   axes
   import matplotlib.pyplot as plt

   plt.gcf().savefig("picture5.png")


   # The structures and energies are identical to before. However, the relative populations changed slightly as there are now fewer conformers in the set.

   # ## More about conformers
   #
   # * Try **CREST** instead of RDKit to generate the initial conformer set
   #
   # * The **Expand** task can be used to expand a set of conformers.
