Conformer Generation for Multiple Molecules

GFNFF Conformers

Apply AMS conformer generation to a system with more than one molecule.

Downloads: Notebook | Script ?

Requires: AMS2027 or later

Related examples

• Conformer Generation, Rescoring, and Filtering

• Generating COSMO-RS Compounds from XYZ Files and SMILES

• COSMO-RS Conformer Workflows for Acetic Acid

• Gasphase IR spectrum from Molecular Dynamics

Related tutorials

• Multiple molecules conformers tutorial

• Conformers tutorial

Related documentation

• AMS conformers utility

• GFNFF overview

Conformer generation for multiple molecules

Lets see how two alanine molecules orient themselves using CREST conformer generation. To do this we will constrain the system in a spherical region using the SphericalWall constraint. We start by setting up a system of two alanine molecules in a relatively small space.

Initial imports

import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scm.plams as plams
from scm.base import ChemicalSystem
from scm.conformers import ConformersJob
from scm.inputs import AMS, Conformers, GFNFF

Alanine dimer

smiles = "CC(N)C(=O)O.CC(N)C(=O)O"
dimer = ChemicalSystem.from_smiles(smiles)  # alternatively you may use packmol to build the dimer
dimer.translate(-dimer.center_of_mass())

Optimize the dimer structure prior to the conformer search.

s = AMS(Task="GeometryOptimization", Engine=GFNFF())
s.GeometryOptimization.Convergence.Quality = "VeryGood"
s.GeometryOptimization.MaxIterations = 1300
s.GeometryOptimization.Method = "Quasi-Newton"

job = plams.AMSJob(molecule=dimer, settings=s)
job.run()
dimer = job.results.get_main_system()
plams.view(dimer, height=300, width=300, direction="along_pca3")

[06.11|09:31:15] JOB plamsjob STARTED
[06.11|09:31:15] JOB plamsjob RUNNING
[06.11|09:31:17] JOB plamsjob FINISHED
[06.11|09:31:17] JOB plamsjob SUCCESSFUL

Calculation setup

To determine the radius of the SphericalWall we measure the size of the initial dimer.

dists = dimer.distance_matrix()
max_dist = np.max(dists)
diameter = 1.33 * max_dist
radius = diameter / 2
print(f"Largest distance between atoms: {max_dist:.3f} ang.")
print(f"Radius: {radius:.3f} ang.")

Largest distance between atoms: 8.460 ang.
Radius: 5.626 ang.

Now we can set up the Crest conformer generation job, with the appropriate spherical wall constraining the molecules close together.

# This seed will used for the starting MD velocities in conformers generation
# Due to numerical aspects it does not guarantee full reproducibility.
seed = random.randint(1, 10000000)
print(f"Seed used for stochastic aspects is {seed}.")

Seed used for stochastic aspects is 2371676.

settings = Conformers(Task="Generate", Engine=GFNFF())
settings.EngineAddons.WallPotential.Enabled = True
settings.EngineAddons.WallPotential.Radius = radius
settings.Generator.Method = "CREST"
settings.Output.KeepWorkDir = True
settings.GeometryOptimization.MaxConvergenceTime = "High"
settings.Generator.CREST.NCycles = 3  # at most 3 CREST cycles for this demo
settings.Generator.RNGSeed = seed
settings.Generator.CREST.NMolDynStepsFactor = 3

Run the conformers job

Now we can run the conformer generation job. This job will run somewhere between 30 minutes and 1 hour.

job = ConformersJob(molecule=dimer, settings=settings)
job.run();
# ConformersJob.load_external("plams_workdir/conformers/conformers.rkf")  # load from disk instead of running the job

[06.11|09:31:21] JOB conformers STARTED
[06.11|09:31:21] JOB conformers RUNNING
[06.11|09:56:26] JOB conformers FINISHED
[06.11|09:56:27] JOB conformers SUCCESSFUL

Now, remove the wall and reoptimize

opt_settings = settings.model_copy(deep=True)
del opt_settings.EngineAddons
opt_settings.Task = "Optimize"
opt_settings.InputConformersSet = [job.results.rkfpath()]
opt_job = ConformersJob(settings=opt_settings)
opt_job.run();

[06.11|09:56:27] JOB conformers STARTED
[06.11|09:56:27] Renaming job conformers to conformers.002
[06.11|09:56:27] JOB conformers.002 RUNNING
[06.11|09:59:15] JOB conformers.002 FINISHED
[06.11|09:59:16] JOB conformers.002 SUCCESSFUL

rkf = opt_job.results.rkfpath()
print(f"Conformers stored in {rkf}")

Conformers stored in /Users/ormrodmorley/Documents/code/ams/amshome_version_bump/userdoc/PythonExamples/conformers-multiple-molecules/plams_workdir/conformers.002/conformers.rkf

Results

relative_energies = opt_job.results.get_relative_energies("kcal/mol")
print(f"Total number of conformers: {len(relative_energies)}")
conformer_table = pd.DataFrame(
    {
        "Conformer Id": np.arange(1, len(relative_energies) + 1),
        "Rel. E (kcal/mol)": relative_energies,
        "Population at 298 K": opt_job.results.get_boltzmann_distribution(298),
    }
)

conformer_table.head(20)

Total number of conformers: 1681

	Conformer Id	Rel. E (kcal/mol)	Population at 298 K
0	1	0.000000	0.083733
1	2	0.007769	0.082642
2	3	0.011192	0.082165
3	4	0.071611	0.074196
4	5	0.103196	0.070342
5	6	0.205166	0.059215
6	7	0.227355	0.057038
7	8	0.305097	0.050020
8	9	0.314652	0.049220
9	10	0.329807	0.047976
10	11	0.470470	0.037833
11	12	0.494660	0.036318
12	13	0.543488	0.033444
13	14	0.934193	0.017290
14	15	0.948115	0.016888
15	16	0.967742	0.016337
16	17	1.091938	0.013246
17	18	1.370614	0.008274
18	19	1.374833	0.008215
19	20	1.400314	0.007869

Here we plot the three lowest-energy conformers.

from scm.conformers.plams.plot import plot_conformers

axes = plot_conformers(opt_job, width=300, height=300, direction="along_pca3")
axes;

You can also open the conformers in AMSmovie to browse the full conformer set:

!amsmovie "{opt_job.results.rkfpath()}"

Filtering out duplicates

The conformerset attribute holds information about the full set of conformers.

conformerset = opt_job.results.conformerset
print("Number of conformers: ", len(conformerset))

Number of conformers:  1681

This conformerset does not contain any duplicates. However, there is no single definition of the term duplicates when applied to conformers. By default, duplicates are defined based on energy and rotational constants, the latter quantifying the 3D shape of the system. The filtering method uses thresholds for energy and rotational constants, which have been thoroughly tested for (mostly small) single molecules systems. If the user feels that the stored conformers are too similar, these thresholds can of course be changed. Below we refilter the stored conformerset using more lenient thresholds (systems with greater differences will be considered duplicates).

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.CREST.EnergyThreshold = 0.2  # default is 0.05
s.Equivalence.CREST.ScaledRotationalConstantSettings.RotationalConstantThreshold = (
    0.01  # default is 0.003
)

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of conformers: ", len(conformerset))

[06.11|09:59:26] JOB conformers STARTED
[06.11|09:59:26] Renaming job conformers to conformers.003
[06.11|09:59:26] JOB conformers.003 RUNNING
[06.11|09:59:50] JOB conformers.003 FINISHED
[06.11|09:59:51] JOB conformers.003 SUCCESSFUL
Number of conformers:  1575

Alternatively, a filtering method can be applied that uses more local comparisons. The native AMS duplicate filter uses the interatomic distance matrix and torsion angles to determine if two structures are duplicates. This approach may be more appropriate for systems with multiple molecules, because the method only takes into account intra-molecular changes. On the other hand, this type of filtering is more time consuming, while in most cases the effect on the conformer set will not be extreme.

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.Method = "AMS"

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of AMS conformers: ", len(conformerset))

[06.11|09:59:51] JOB conformers STARTED
[06.11|09:59:51] Renaming job conformers to conformers.004
[06.11|09:59:51] JOB conformers.004 RUNNING
[06.11|10:00:44] JOB conformers.004 FINISHED
[06.11|10:00:44] JOB conformers.004 SUCCESSFUL
Number of AMS conformers:  1678

We also have a duplicate filter specifically for multiple molecules, where only interactions between molecules are compared, while differences in the molecules themselves are considered irrelevant.

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.Method = "Interactions"

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of differently interacting geometries: ", len(conformerset))

[06.11|10:00:44] JOB conformers STARTED
[06.11|10:00:44] Renaming job conformers to conformers.005
[06.11|10:00:44] JOB conformers.005 RUNNING
[06.11|10:00:56] JOB conformers.005 FINISHED
[06.11|10:00:56] JOB conformers.005 SUCCESSFUL
Number of differently interacting geometries:  209

See also

Related examples

• Conformer Generation, Rescoring, and Filtering

• Generating COSMO-RS Compounds from XYZ Files and SMILES

• COSMO-RS Conformer Workflows for Acetic Acid

• Gasphase IR spectrum from Molecular Dynamics

Related tutorials

• Multiple molecules conformers tutorial

• Conformers tutorial

Related documentation

• AMS conformers utility

• GFNFF overview

Python Script

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

# ## Conformer generation for multiple molecules
# 
# Lets see how two alanine molecules orient themselves using CREST conformer generation.
# To do this we will constrain the system in a spherical region using the `SphericalWall` constraint.
# We start by setting up a system of two alanine molecules in a relatively small space.

# ## Initial imports

import random
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import scm.plams as plams
from scm.base import ChemicalSystem
from scm.conformers import ConformersJob
from scm.inputs import AMS, Conformers, GFNFF


# ## Alanine dimer

smiles = "CC(N)C(=O)O.CC(N)C(=O)O"
dimer = ChemicalSystem.from_smiles(smiles)  # alternatively you may use packmol to build the dimer
dimer.translate(-dimer.center_of_mass())


# Optimize the dimer structure prior to the conformer search.

s = AMS(Task="GeometryOptimization", Engine=GFNFF())
s.GeometryOptimization.Convergence.Quality = "VeryGood"
s.GeometryOptimization.MaxIterations = 1300
s.GeometryOptimization.Method = "Quasi-Newton"

job = plams.AMSJob(molecule=dimer, settings=s)
job.run()
dimer = job.results.get_main_system()
plams.view(dimer, height=300, width=300, direction="along_pca3", picture_path="picture1.png")


# ## Calculation setup

# To determine the radius of the `SphericalWall` we measure the size of the initial dimer.

dists = dimer.distance_matrix()
max_dist = np.max(dists)
diameter = 1.33 * max_dist
radius = diameter / 2
print(f"Largest distance between atoms: {max_dist:.3f} ang.")
print(f"Radius: {radius:.3f} ang.")


# Now we can set up the Crest conformer generation job, with the appropriate spherical wall constraining the molecules close together. 

# This seed will used for the starting MD velocities in conformers generation
# Due to numerical aspects it does not guarantee full reproducibility.
seed = random.randint(1, 10000000)
print(f"Seed used for stochastic aspects is {seed}.")


settings = Conformers(Task="Generate", Engine=GFNFF())
settings.EngineAddons.WallPotential.Enabled = True
settings.EngineAddons.WallPotential.Radius = radius
settings.Generator.Method = "CREST"
settings.Output.KeepWorkDir = True
settings.GeometryOptimization.MaxConvergenceTime = "High"
settings.Generator.CREST.NCycles = 3  # at most 3 CREST cycles for this demo
settings.Generator.RNGSeed = seed
settings.Generator.CREST.NMolDynStepsFactor = 3


# ## Run the conformers job

# Now we can run the conformer generation job. This job will run somewhere between 30 minutes and 1 hour.

job = ConformersJob(molecule=dimer, settings=settings)
job.run();
# ConformersJob.load_external("plams_workdir/conformers/conformers.rkf")  # load from disk instead of running the job


# Now, remove the wall and reoptimize

opt_settings = settings.model_copy(deep=True)
del opt_settings.EngineAddons
opt_settings.Task = "Optimize"
opt_settings.InputConformersSet = [job.results.rkfpath()]
opt_job = ConformersJob(settings=opt_settings)
opt_job.run();


rkf = opt_job.results.rkfpath()
print(f"Conformers stored in {rkf}")


# ## Results

relative_energies = opt_job.results.get_relative_energies("kcal/mol")
print(f"Total number of conformers: {len(relative_energies)}")
conformer_table = pd.DataFrame(
    {
        "Conformer Id": np.arange(1, len(relative_energies) + 1),
        "Rel. E (kcal/mol)": relative_energies,
        "Population at 298 K": opt_job.results.get_boltzmann_distribution(298),
    }
)

print(conformer_table.head(20))


# Here we plot the three lowest-energy conformers.

from scm.conformers.plams.plot import plot_conformers

axes = plot_conformers(opt_job, width=300, height=300, direction="along_pca3")
axes;
import matplotlib.pyplot as plt
plt.gcf().savefig("picture2.png")


# You can also open the conformers in AMSmovie to browse the full conformer set:

# get_ipython().system('amsmovie "{opt_job.results.rkfpath()}"')


# ### Filtering out duplicates
# The `conformerset` attribute holds information about the full set of conformers.

conformerset = opt_job.results.conformerset
print("Number of conformers: ", len(conformerset))


# This conformerset does not contain any duplicates. However, there is no single definition of the term duplicates when applied to conformers. By default, duplicates are defined based on energy and rotational constants, the latter quantifying the 3D shape of the system. The filtering method uses thresholds for energy and rotational constants, which have been thoroughly tested for (mostly small) single molecules systems. If the user feels that the stored conformers are too similar, these thresholds can of course be changed. Below we refilter the stored conformerset using more lenient thresholds (systems with greater differences will be considered duplicates).

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.CREST.EnergyThreshold = 0.2  # default is 0.05
s.Equivalence.CREST.ScaledRotationalConstantSettings.RotationalConstantThreshold = (
    0.01  # default is 0.003
)

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of conformers: ", len(conformerset))


# Alternatively, a filtering method can be applied that uses more local comparisons. The native AMS duplicate filter uses the interatomic distance matrix and torsion angles to determine if two structures are duplicates. This approach may be more appropriate for systems with multiple molecules, because the method only takes into account intra-molecular changes. On the other hand, this type of filtering is more time consuming, while in most cases the effect on the conformer set will not be extreme.

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.Method = "AMS"

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of AMS conformers: ", len(conformerset))


# We also have a duplicate filter specifically for multiple molecules, where only interactions between molecules are compared, while differences in the molecules themselves are considered irrelevant.

s = Conformers(Task="Filter", InputConformersSet=[rkf])
s.Equivalence.Method = "Interactions"

filter_job = ConformersJob(settings=s)
filter_job.run()
conformerset = filter_job.results.conformerset
print("Number of differently interacting geometries: ", len(conformerset))