Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Interpolation mechanics

  • Explain what is interpolation mechanics

  • Deriving the Force-Field

  • Running a dynamics

Force field generation

We define first a VeloxChem Molecule for which the interpolation mechaics force field (IM-MM) will be constructed.

import veloxchem as vlx
# Initialize the Molecule
molecule_xyz = '''16

C              1.354430000000         0.723890000000         0.062850000000
C              0.075030000000         1.193180000000        -0.131980000000
S             -1.063880000000        -0.076840000000        -0.238910000000
C              0.145050000000        -1.275960000000        -0.023130000000
C              1.392450000000        -0.699780000000         0.133940000000
C             -0.161190000000        -2.684380000000         0.016020000000
C              0.347150000000        -3.612800000000         0.907410000000
C             -0.139800000000        -4.931010000000         0.664450000000
C             -0.992970000000        -4.972470000000        -0.413030000000
S             -1.212210000000        -3.430270000000        -1.117120000000
H              2.219640000000         1.368740000000         0.151510000000
H             -0.248470000000         2.221020000000        -0.218570000000
H              2.299910000000        -1.273800000000         0.280250000000
H              1.028740000000        -3.357480000000         1.709320000000
H              0.125220000000        -5.802850000000         1.250990000000
H             -1.502580000000        -5.833680000000        -0.824080000000
'''


molecule = vlx.Molecule.from_xyz_string(molecule_xyz)
molecule.show(atom_indices=False)
Loading...

We then initialize the QM-Driver determining the underlying quality of the interpolation data.

# Initialize the basis set and the QM driver
basis_set_label = "def2-svp"
basis = vlx.MolecularBasis.read(molecule, basis_set_label)
qm_driver = vlx.ScfRestrictedDriver()
qm_driver.xcfun = "b3lyp"

Define dyanmics and interpolation settings which are crucial for the quality in timing of the database. Important variables:

.desired_point_density -> how many datapoints should the database contain for each conformer.
.converged_cycle -> after how many cycles the construction for 1 conformer terminates (1 cycle: 3000 steps without database expansion).
.energy_thrsh -> allowed devaition of the IM energy and QM reference calculations.

As this notebook is intended to be a practival example on how to use the IMFroceFieldGenerator but keeping the computational effort at a reasoable level, the shown settings have been used. For production run we recommend using the following settings:

.nsteps = 50000
.snapshots = 10000
.desired_point_density = 80
.converged_cycle = 5
.energy_threshold = 2.0 kcal/mol
# Initialize the IM force field generator with:
# dihedrals: list of dihedral angles to sample
# sampling_structures: number of structures to sample
# ensemble: NVT or NVE
# temperature: temperature of the ensemble
# nsteps: number of steps to run the dynamics
# snapshots: number of snapshots to save

imffgen_drv = vlx.IMForceFieldGenerator(qm_driver)
imffgen_drv.dihedrals_dict = [((2, 3, 5, 9), 10)]
imffgen_drv.ensemble = 'NVT'
imffgen_drv.temperature = 250
imffgen_drv.nsteps = 10000
imffgen_drv.snapshots = 5000
imffgen_drv.converged_cycle = 3
imffgen_drv.energy_threshold = 2.0
imffgen_drv.desired_point_density = 10
imffgen_drv.nstruc_to_confirm_database_quality = 10
imffgen_drv.minimize = False
IMPORTANT: IM ForceFieldFile is initalized as im_database.h5

The compute function starts the construction of a new or given database. The construction run ends by confirming the qulity of the database and the result are stored.

Note that here we provide the constructed database (‘im_database.h5’) as the computational cost of the second derivatives are very demanding.

# Run the database consruction
imffgen_drv.imforcefieldfile = "im_database.h5"
im_results = imffgen_drv.compute(molecule, basis)

ff_gen = vlx.MMForceFieldGenerator()
ff_gen.create_topology(molecule)

To use the constructed IM force field one has to initialize the InterpolationDriver() with user-defined settings (recommended). The current implementation is restricted to well seperated IM and MM regions and is not yet included into a linking atom framework.

z_matrix = imffgen_drv.define_z_matrix(molecule)
im_driver = vlx.InterpolationDriver(z_matrix)
interpolation_settings = { 'interpolation_type':'shepard', 
                            'exponent_p':'2',
                            'exponent_q':'2', 
                            'confidence_radius':'0.5',
                            'imforcefieldfile':'im_database.h5'}
im_driver.update_settings(interpolation_settings)

openmmdyn = vlx.OpenMMDynamics()
openmmdyn.create_system_from_molecule(molecule, ff_gen, qm_atoms='all')
openmmdyn.run_qmmm(im_driver, im_driver)