Example 7: Example of multiple system parametrization¶
Mapping of different topologies¶
ParaMol is able to parametrize multiple systems at the same time. Even though this is possible to do without applying any symmetry constraint, in the context of force field development it is often desirable to derive the parameters for a given term type and for a set of molecules. Therefore, in this example we are going to concomitantly optimize the barrier height of the dihedral type hc-c3-c3-hc for ethane and propane.
NOTE: When optimizing more than once system at once, the mapping of the symmetry groups of different systems has to be performed manually so that there is a one-to-one correspondence between symmetry group and term type. The automatic mapping of symmetries for different systems is a feature that will likely be implemented in future ParaMol versions.
The first step is to write out the ParaMol Force Field files so that the correct symmetries are set, which can be done using the following code:
# ParaMol imports
from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
# --------------------------------------------------------- #
# Preparation #
# --------------------------------------------------------- #
# Create the OpenMM engine for ethane
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='ethane.prmtop', crd_format='AMBER', crd_file='ethane.inpcrd')
# Create ethane ParaMol System
ethane = ParaMolSystem(name="ethane", engine=openmm_system, n_atoms=8)
# Create ParaMol's force field representation
ethane.force_field.create_force_field()
# Create the OpenMM engine for propane
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='propane.prmtop', crd_format='AMBER', crd_file='propane.inpcrd')
# Create propane ParaMol System
propane = ParaMolSystem(name="propane", engine=openmm_system, n_atoms=11)
# Create ParaMol's force field representation
propane.force_field.create_force_field()
# --------------------------------------------------------- #
# Save ParaMol ForceField #
# --------------------------------------------------------- #
ethane.force_field.write_ff_file("ethane.ff")
propane.force_field.write_ff_file("propane.ff")
Original ethane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 1 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
0 0 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
2 1 0 4 7 3.00000000 0.00000000 0.62760000 0 0 0 X
3 2 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
4 2 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
5 2 0 4 7 3.00000000 0.00000000 0.62760000 0 0 0 X
6 3 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
7 3 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
8 3 0 4 7 3.00000000 0.00000000 0.62760000 0 0 0 X
...
Original propane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 0 4 7 8 3.00000000 0.00000000 0.66944000 0 0 0 X
1 0 4 7 9 3.00000000 0.00000000 0.66944000 0 0 0 X
2 0 4 7 10 3.00000000 0.00000000 0.66944000 0 0 0 X
3 1 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
4 1 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
5 1 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
6 2 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
7 2 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
8 2 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
9 3 0 4 5 3.00000000 0.00000000 0.62760000 0 0 0 X
10 3 0 4 6 3.00000000 0.00000000 0.62760000 0 0 0 X
11 3 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
12 5 4 7 8 3.00000000 0.00000000 0.62760000 0 0 0 X
13 5 4 7 9 3.00000000 0.00000000 0.62760000 0 0 0 X
14 5 4 7 10 3.00000000 0.00000000 0.62760000 0 0 0 X
15 6 4 7 8 3.00000000 0.00000000 0.62760000 0 0 0 X
16 6 4 7 9 3.00000000 0.00000000 0.62760000 0 0 0 X
17 6 4 7 10 3.00000000 0.00000000 0.62760000 0 0 0 X
...
As can be seen in the ParaMol Force Field files, all term types belong to the default symmetry group, i.e., “X”, which means that they do not possess any symmetry. We will set the symmetry label of the dihedral types hc-c3-c3-hc to “T0” (see modified ParaMol Force Fields below). In this way, we will be able to find the best barrier height that minimizes the objective function for both ethane and propane.
Modified ethane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 1 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
1 1 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
2 1 0 4 7 3.00000000 0.00000000 0.62760000 0 0 1 T0
3 2 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
4 2 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
5 2 0 4 7 3.00000000 0.00000000 0.62760000 0 0 1 T0
6 3 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
7 3 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
8 3 0 4 7 3.00000000 0.00000000 0.62760000 0 0 1 T0
...
Modified propane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 0 4 7 8 3.00000000 0.00000000 0.66944000 0 0 0 X
1 0 4 7 9 3.00000000 0.00000000 0.66944000 0 0 0 X
2 0 4 7 10 3.00000000 0.00000000 0.66944000 0 0 0 X
3 1 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
4 1 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
5 1 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
6 2 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
7 2 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
8 2 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
9 3 0 4 5 3.00000000 0.00000000 0.62760000 0 0 1 T0
10 3 0 4 6 3.00000000 0.00000000 0.62760000 0 0 1 T0
11 3 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
12 5 4 7 8 3.00000000 0.00000000 0.62760000 0 0 1 T0
13 5 4 7 9 3.00000000 0.00000000 0.62760000 0 0 1 T0
14 5 4 7 10 3.00000000 0.00000000 0.62760000 0 0 1 T0
15 6 4 7 8 3.00000000 0.00000000 0.62760000 0 0 1 T0
16 6 4 7 9 3.00000000 0.00000000 0.62760000 0 0 1 T0
17 6 4 7 10 3.00000000 0.00000000 0.62760000 0 0 1 T0
...
Simultaneous parametrization of multiple systems¶
Now that we have done the necessary changes in the ParaMol Force Field files, we are ready to perform the conformational sampling, ab initio properties calculation and subsequent parameters’ optimization. This can be done using the following script:
# ParaMol imports
from ParaMol.System.system import *
from ParaMol.MM_engines.openmm import *
from ParaMol.QM_engines.qm_engine import *
# ParaMol Tasks imports
from ParaMol.Tasks.parametrization import *
from ParaMol.Tasks.torsions_scan import *
from ParaMol.Utils.settings import *
from ParaMol.Utils.conformational_sampling import *
# --------------------------------------------------------- #
# Settings #
# --------------------------------------------------------- #
# Create ParaMol settings instance
paramol_settings = Settings()
# The objective function will contain an energy and a regularization term
paramol_settings.properties["include_energies"] = True
paramol_settings.properties["include_forces"] = False
paramol_settings.properties["include_regularization"] = True
# --------------------------------------------------------- #
# Preparation #
# --------------------------------------------------------- #
# Create the OpenMM engine for ethane
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='ethane.prmtop', crd_format='AMBER', crd_file='ethane.inpcrd')
# Create ethane ParaMol System
ethane = ParaMolSystem(name="ethane", engine=openmm_system, n_atoms=8)
# Create ParaMol's force field representation
ethane.force_field.create_force_field(ff_file="ethane_sym.ff")
# Create the OpenMM engine for propane
openmm_system = OpenMMEngine(init_openmm=True, topology_format='AMBER', top_file='propane.prmtop', crd_format='AMBER', crd_file='propane.inpcrd')
# Create propane ParaMol System
propane = ParaMolSystem(name="propane", engine=openmm_system, n_atoms=11)
# Create ParaMol's force field representation
propane.force_field.create_force_field(ff_file="propane_sym.ff")
# --------------------------------------------------------- #
# Set the QM Engine #
# --------------------------------------------------------- #
# Create the ASE calculator
from ase.calculators.dftb import *
calc = Dftb(Hamiltonian_='DFTB', # line is included by default
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_H='s',
Hamiltonian_MaxAngularMomentum_O='p',
Hamiltonian_MaxAngularMomentum_C='p',
Hamiltonian_MaxAngularMomentum_N="p",
Hamiltonian_Dispersion="DftD3 { \n s6=1.000 \n s8=0.5883 \n Damping = BeckeJohnson { \n a1=0.5719 \n a2=3.6017 \n } \n }",
Hamiltonian_SCC='Yes',
Hamiltonian_SCCTolerance=1e-8, )
# Alternative, we could set the calculator in the settings
paramol_settings.qm_engine["ase"]["calculator"] = calc
# --------------------------------------------------------- #
# Conformational Sampling #
# --------------------------------------------------------- #
systems = [ethane, propane]
# Perform conformational sampling
conformational_sampling(paramol_settings, systems, 1000, 1000)
# --------------------------------------------------------- #
# Parameter Optimization #
# --------------------------------------------------------- #
parametrization = Parametrization()
systems, parameter_space, objective_function, optimizer = parametrization.run_task(paramol_settings, systems)
# Write final ParaMol FFs
ethane.force_field.write_ff_file("ethane_optimized.ff")
propane.force_field.write_ff_file("propane_optimized.ff")
Finally, as can be seen in the final ParaMol Force Field files, a new barrier height value was found that best suits both systems. Specifically, GAFF slightly overestimates the barrier height of this torsions with respect to the SCC-DFTB-D3 level of theory, since, after re-parametrization, its values has decreased from 0.6276 kj/mol to 0.4690 kj/mol.
Final ethane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 1 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
1 1 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
2 1 0 4 7 3.00000000 0.00000000 0.46899172 0 0 1 T0
3 2 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
4 2 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
5 2 0 4 7 3.00000000 0.00000000 0.46899172 0 0 1 T0
6 3 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
7 3 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
8 3 0 4 7 3.00000000 0.00000000 0.46899172 0 0 1 T0
...
Final propane ParaMol Force Field file.¶
...
PeriodicTorsionForce 2
0 0 4 7 8 3.00000000 0.00000000 0.66944000 0 0 0 X
1 0 4 7 9 3.00000000 0.00000000 0.66944000 0 0 0 X
2 0 4 7 10 3.00000000 0.00000000 0.66944000 0 0 0 X
3 1 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
4 1 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
5 1 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
6 2 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
7 2 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
8 2 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
9 3 0 4 5 3.00000000 0.00000000 0.46899172 0 0 1 T0
10 3 0 4 6 3.00000000 0.00000000 0.46899172 0 0 1 T0
11 3 0 4 7 3.00000000 0.00000000 0.66944000 0 0 0 X
12 5 4 7 8 3.00000000 0.00000000 0.46899172 0 0 1 T0
13 5 4 7 9 3.00000000 0.00000000 0.46899172 0 0 1 T0
14 5 4 7 10 3.00000000 0.00000000 0.46899172 0 0 1 T0
15 6 4 7 8 3.00000000 0.00000000 0.46899172 0 0 1 T0
16 6 4 7 9 3.00000000 0.00000000 0.46899172 0 0 1 T0
17 6 4 7 10 3.00000000 0.00000000 0.46899172 0 0 1 T0
...