Native Engine

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 number of steps of the 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: