Molecular mechanics - VeloxChem Docs

import veloxchem as vlx

Let us assume we have a molecular for which we wish to construct a molecular mechanics (MM) force field.

molecule.show(atom_indices=True, width=600, height=450)

Force field generation¶

An initial force field is created based on calculated RESP charges and tabulated GAFF parameters.

ff_gen = vlx.MMForceFieldGenerator()

ff_gen.create_topology(molecule)

Force field reparameterization¶

Rotatable bonds have been identified and stored by the force field generator.

print(ff_gen.rotatable_bonds)

We will focus on the reparameterization of the barrier around the rotatable bond [21, 22], and make use of QM relaxed-scan reference data from file 16-21-22-27.xyz.

reparam_results = ff_gen.reparameterize_dihedrals(
    rotatable_bond=(21, 22),
    scan_file="../input_files/16-21-22-27.xyz",
    visualize=True,
    verbose=False,
)

The rotational barrier is now much improved but the relative energies between the two conformer minima are still not well reproduced, leading to significant errors in the statistical dihedral distributions. As a remedy to this situation, an additional dihedral potential can be added.

ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=180.0, periodicity=1)

We can now restart the reparametrization procedure with this added dihedral angle.

reparam_results = ff_gen.reparameterize_dihedrals(
    rotatable_bond=(21, 22),
    scan_file="../input_files/16-21-22-27.xyz",
    visualize=True,
    verbose=False,
    initial_validation=False,
)

Another dihedral angle can be added to reach an even better agreement with the QM reference data.

ff_gen.add_dihedral((16, 21, 22, 27), barrier=0.0, phase=0.0, periodicity=4)

reparam_results = ff_gen.reparameterize_dihedrals(
    rotatable_bond=(21, 22),
    scan_file="../input_files/16-21-22-27.xyz",
    visualize=True,
    verbose=False,
    initial_validation=False,
)

When a satisfactory quality has been reached in the force field, the files required to run GROMACS molecular dynamics simulations can be generated.

ff_gen.write_gromacs_files("HS-276_final", "MOL")