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.

Downloads: Notebook | Script ?

Requires: AMS2026 or later

Related examples
Related tutorials
Related documentation
../_images/conformers_generation_31_0_3770c31f.png

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 generated from notebook

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 generated from notebook

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 generated from notebook

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 generated from notebook

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 generated from notebook

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

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.