Simulations

A general procedure for using a force field generated by DFF is as follows.

Assigning a Force Field

Prepare or import a simulation model. DFF can be used to build three-dimensional cubic liquid boxes and liquid interfaces, though this functionality is intended primarily for parameterization and validation. For complex molecular systems, an external software package may be used to build the simulation models. DFF supports most popular file formats, including .pdb, .mol2, and .car files.

It is important to maintain a consistent assignment of topologies and formal charges in order to assign atom types. If atomic coordinates are provided with reasonable precision, DFF can automatically assign topologies and formal charges.

Once topologies and formal charges are assigned, the next step is to assign atom types and force field parameters using TEAMFF. Multiple force fields of the same type can be used together, but they are called in the order specified. All valence parameters for a given molecule must come from one force field; intra-molecular parameters cannot be mixed across force fields. Only inter-molecular parameters (nonbond terms) may be combined from different force fields using the combination rule. In practice, DFF scans all molecules against the selected force field. If any parameters are missing for a molecule, it is not assigned and is passed to the next force field in the list.

When all parameters are found, the atom types and parameters—possibly sourced from different force field tables—are unified into one consistent force field for the simulation model. A new set of atom types corresponding to this unified table is automatically generated and assigned to the simulation model. The model and the unified force field can then be used directly in simulations.

If parameters are still missing, TEAMFF will report the missing terms, which can be used for parametrization.

Native Functions

DFF includes a built-in simulation engine that can be used to optimize and validate parameters. Simulation jobs can be run interactively so that structures and energies are updated dynamically on the screen.

Structure Relaxation

For models with high internal tension due to distorted structures, it is often necessary to relax the system before performing molecular dynamics (MD) simulations or energy minimization. A typical example is tangled polymer chains, where poor initial construction leads to excessively high internal energy. To address such cases, a van der Waals (vdW) growing molecular dynamics method is implemented.

In this approach, a series of short NVT MD simulations are run sequentially. In each simulation, the vdW terms—Lennard-Jones (LJ) parameters and atomic charges—are scaled by a factor ranging from 0 to 1. At the start, vdW interactions are completely removed, so the dynamics are governed only by valence interactions, allowing intra-molecular structures to stabilize. As the strength of vdW interactions is gradually increased, both intra- and inter-molecular packing are progressively adjusted.

To apply this method, the key parameters to define are:

• the number of vdW growth steps, each step is a stage.
• the step-size and steps of molecular dynamics simulation at each stage.

Optimizers

Two different algorithms—the conjugate gradient method and the variable metric method—are implemented in DFF for energy minimization. The conjugate gradient method is based on the Polak–Ribiere approach, which is similar to the Fletcher–Reeves method. The variable metric method, also known as the quasi-Newton method, is an iterative approach that uses both the current gradient and the history of previous steps to determine the next move. For a quadratic functional form in N dimensions, both methods converge to the exact minimum in N steps.

These algorithms have different memory requirements. The conjugate gradient method requires storage on the order of N, while the variable metric method requires storage on the order of N × N.

Restrained energy minimization can also be performed to explore energy profiles. This is achieved by adding restraint terms to the total energy function on selected internal coordinates. The restraint function is harmonic, with user-defined reference values and force constants. In reported results, the restraint energy is subtracted from the total energy.

Molecular Dynamics

DFF implements the velocity Verlet algorithm for molecular dynamics. Choosing the appropriate integration time step is critical. A smaller step size increases accuracy but reduces the simulated time span per unit of execution time, while a larger step size does the opposite. In most molecular systems, the highest vibrational frequency is associated with X–H bond stretching, with a period of approximately 10⁻¹⁴ s. Therefore, a typical integration time step is about 1 femtosecond (fs).

A dynamic run consists of three stages: initialization, equilibration, and simulation. Initialization assigns initial positions and velocities to all atoms, though this step can be skipped in a continuous simulation. Equilibration allows the system to reach a configuration consistent with the target temperature and pressure. At equilibrium, system properties such as temperature, potential energy, and kinetic energy remain constant on average. Once equilibrium is established, a production simulation can be run, and the resulting data analyzed.

The available simulation types are:

• NVE (microcanonical) simulation: This ensemble is obtained by solving Newton’s equations without temperature or pressure controls. Except for rounding and truncation errors, the total energy is conserved.
• NVT (canonical) simulation: This ensemble controls temperature, either through direct velocity scaling or coupling to a thermal bath. Volume remains constant throughout the run.
• NPT (isothermal–isobaric) simulation: This ensemble controls both temperature and pressure. Unit cell vectors are allowed to change, and pressure is adjusted by varying volume. This ensemble is commonly used when accurate pressure, volume, and density are important. It can also be applied during equilibration before switching to a constant-volume or constant-energy ensemble for data collection.

Temperature can be controlled using two methods:

• Direct velocity scaling: Atom velocities are scaled uniformly whenever the instantaneous temperature falls outside a specified window, ensuring the system matches the target temperature.
• Andersen stochastic method: At intervals proportional to N²ᐟ³, where N is the number of atoms, each atom’s velocity is reassigned from a Maxwell–Boltzmann distribution, mimicking collisions with a heat bath. Between these stochastic events, the system evolves under constant energy, allowing time correlation functions to be computed.

Pressure is calculated using the virial theorem. Pressure control is achieved by adjusting particle coordinates and the unit cell size under periodic boundary conditions. The Berendsen method couples the system to a pressure “bath” that maintains the target pressure by scaling atomic coordinates and cell edges at each step. The strength of the coupling depends on the user-defined normalized compressibility. Note that this method scales the cell uniformly, changing its size but not its shape.

DFF applies periodic boundary conditions for condensed-phase simulations, using the minimum image convention. Nonbond interactions can be evaluated with atom-based or charge-group-based cutoffs, both of which can be supplemented with cutoff tail corrections. The charge-group cutoff method is particularly effective for organic molecules, where small neutral charge groups can be defined. As shown in the following table, this method yields results nearly identical to those obtained with Ewald summation, but with significantly reduced computational cost:

External Functions

Two external simulation packages, LAMMPS and GROMACS, are included with DFF distributions. These are multi-threaded parallel (MTP) versions released after 2021. Force field validation tasks such as energy minimization, and NVT, NPT, and NVE molecular dynamics can be launched directly from DFF using these programs.

Conversions

Because different software packages use different atom type conventions, DFF translates atom types into the formats recognized by the target program. This translation is a direct one-to-one mapping between DFF atom types and unique string identifiers used by the external software.

Since potential functions may be implemented differently from their published forms, some parameters must also be converted into the closest functional equivalents before export.

Dihedral in DREIDING

The dihedral potential in DREIDING is expressed as:

In LAMMPS, the closest function is the CHARMM dihedral_style:

The conversion between the two forms is given by:

Thus, the force constant remains unchanged (Eq. 3), the periodicity is preserved (Eq. 4), and the phase shift is adjusted by (Eq. 5).

Out-of-Plane in DREIDING

Two different functions are used for the out-of-plane (improper) potential in DREIDING.

For non-planar structures ():

For planar structures ():

LAMMPS uses a harmonic function (class2 improper style):

To achieve equivalence between the two functional forms, the following conversions are applied:

If :

If :

In summary, the reference angle remains unchanged (Eq. 9). For non-planar cases, the force constant is scaled by (Eq. 10), while for planar cases, it is halved (Eq. 11).