LIBEFP/EFPMD user manual

General functionality

EfpMD is a molecular simulation program for LibEFP. It supports EFP-only (no QM region) calculations. The following functionality is currently available:

  • single point energy

  • gradient

  • geometry optimization

  • semi-numerical Hessian and normal mode analysis

  • molecular dynamics simulations in microcanonical (NVE), canonical (NVT), and isobaric-isothermal (NPT) ensembles

Additioanally, EfpMD provides calculations of

  • electrostatic potential on fragment atoms

  • electric field on fragment atoms

  • pairwise interaction energy decomposition analysis (PIEDA)

EfpMD supports calculations with periodic boundary conditions in arbitrary parallelopiped unit cells. Switching functions can be used for treating long-range interactions.

Since version 2.0, EfpMD can perform energy and geoemtry optimization calculations utilizing neural network (NN) potentials through interface with LibTorch.

Additional examples of input files can be found in the tests directory in source code archive.

Input file format

The input file generally consists of two parts: a part containing parameters of the requested calculation (EFPMD parameters), and the second part containing coordinates of a molecular system (Molecular geometry).

An example EfpMD input file for performing geometry optimization of a 4-water cluster opt_1.in.

 1# water tetramer
 2
 3run_type opt
 4terms elec pol disp xr
 5max_steps 100
 6coord points
 7fraglib_path ../fraglib
 8
 9fragment h2o_l
10   -6.6939 -0.7053  0.2031
11   -7.5290 -0.1798  0.2855
12   -6.9161 -1.6539  0.0273
13fragment h2o_l
14   -3.0446  1.4217  0.1994
15   -3.8439  0.8389  0.2372
16   -3.3220  2.3687  0.2792
17fragment h2o_l
18   -1.9443 -2.0764  0.1287
19   -0.9588 -1.9865  0.1596
20   -2.3592 -1.3491  0.6570
21fragment h2o_l
22   -5.5235 -4.1554  0.2711
23   -5.4795 -4.8791 -0.4031
24   -4.6654 -3.6617  0.2824

Note

Lines beginning with the # symbol are ignored during input parsing.

EFPMD parameters

Basic parameters

Type of the simulation

run_type [sp|grad|hess|opt|md|efield|elpot|gtest|etest]

  • sp - single point energy calculation.

  • grad - energy gradient calculation.

  • hess - semi-numerical Hessian calculation and normal mode analysis.

  • opt - geometry optimization.

  • md - molecular dynamics simulation.

  • efield - compute and print electric field on all atoms.

  • elpot - compute and print electrostatic potential on all atoms.

  • frag_elpot - compute and print electrostatic potential on all atoms of a special fragment.

  • gtest - compute and compare numerical and analytical gradients.

  • etest - compute and compare total energy.

Default value: sp

Format of fragment input

coord [xyzabc|points|rotmat|atoms]

  • xyzabc - Coordinates of the center of mass and Euler angles.

  • points - Coordinates of three atoms for each fragment.

  • rotmat - Coordinates of the center of mass and rotation matrix.

  • atoms - Atom names and cartesian coordinates of all atoms.

Default value: points

Note

See Fragment input specification for more details.

Warning

Default changed to points from xyzabc in version 1.7.3

EFP energy terms

terms [elec pol disp xr qq lj]

  • elec - Include electrostatics energy.

  • pol - Include polarization energy.

  • disp - Include dispersion energy.

  • xr - Include exchange repulsion energy.

  • qq - Include point-charge Coulomb energy due to MM charges (MM_CHARGE section in .efp file is required).

  • lj - Include Lennard-Jones interaction term (MM_LJ section in .efp file is required).

Default value: elec pol disp xr

Short-range damping parameters

Electrostatic damping type

elec_damp [screen|overlap|off]

  • screen - Damping formula based on SCREEN group in the EFP potential.

  • overlap - Overlap-based damping formula. This damping correction is printed

as charge penetration energy. - off - No electrostatic damping.

Default value: screen

Dispersion damping type

disp_damp [tt|overlap|off]

  • tt - Damping based on the formula by Tang and Toennies, with hard-coded coefficient \(b = 1.5\).

  • overlap - Overlap-based dispersion damping.

  • off - No dispersion damping.

Default value: overlap

Polarization damping type

pol_damp [tt|off]

  • tt - Tang and Toennies like damping formula.

  • off - No polarization damping.

Default value: tt

Long-range cutoffs

Enable cutoff for fragment/fragment interactions

enable_cutoff [true|false]

Default value: false

Cutoff distance for fragment/fragment interactions

swf_cutoff <value>

Default value: 10.0

Unit: Angstrom

Cutoff distance for exchange-repulsion interactions between fragments

xr_cutoff <value>

Default value: swf_cutoff

Unit: Angstrom

Polarization solver

Polarization solver

pol_driver [iterative|direct]

  • iterative - Iterative solution of system of linear equations for polarization induced dipoles.

  • direct - Direct solution of system of linear equations for polarization induced dipoles. This solver does not have convergence issues but is unsuitable for large systems (more than 2000 polarizable points). The direct solver is not parallelized.

Default value: iterative

Geometry optimization

Maximum number of optimization and MD steps

max_steps <number>

Default value: 100

Note

max_steps specifies maximum number of steps for both geometry optimization and MD.

Gradient tolerance

opt_tol <value>

Default value: 0.0003

Unit: Hartree/Bohr

Energy tolerance

opt_energy_tol <value>

Default value: 0.000001

Unit: Hartree

Note

Optimization will stop when maximum gradient component is less than opt_tol and energy change is less than opt_energy_tol.

Hessian calculation

Hessian accuracy

hess_central [true|false]

Default value: false

If hess_central is true, then the more accurate central differences will be used for numerical Hessian calculation. Otherwise, only forward differences will be used. Note that central differences require twice as many gradient calculations.

Numerical differentiation step length for distances

num_step_dist <value>

Default value: 0.001

Unit: Angstrom

Numerical differentiation step length for angles

num_step_angle <value>

Default value: 0.01

Unit: Radian

Molecular dynamics simulations

Ensemble

ensemble [nve|nvt|npt]

  • nve - Microcanonical ensemble.

  • nvt - Canonical ensemble with Nose-Hoover thermostat. For the description of the algorithm, see _Phys. Rev. A 31, 1695 (1985)_.

  • npt - Isobaric-isothermal ensemble. This also sets enable_pbc to `true. For the description of the algorithm, see _Mol. Phys. 78, 533 (1993)_.

Default value: nve

Time step

time_step <value>

Unit: Femtosecond

Default value: 1.0

Print step: Number of steps between outputs of the system state.

print_step <value>

Default value: 1

Assign initial velocities

velocitize [true|false]

Default value: false

If true, random initial velocities will be assigned to fragments using Gaussian distribution. Velocity magnitudes are chosen so that initial temperature of the system is approximately equal to the target simulation temperature.

Simulation temperature: target simulation temperature for NVT and NPT thermostats.

temperature <value>

Unit: Kelvin

Default value: 300.0

Simulation pressure

Target simulation pressure for NPT barostat. Note that whether or not pressure coupling is used, the pressure value will oscillate significantly. Fluctuations on the order of hundreds of bar are typical. This variation is entirely normal due to the fact that pressure is a macroscopic property and can only be measured properly as time average. The actual variations of instantaneous pressure depend on the size of the system and the values of barostat parameters.

pressure <value>

Unit: Bar

Default value: 1.0

Thermostat parameter: temperature relaxation time parameter of the Nose-Hoover thermostat

thermostat_tau <value>

Unit: Femtosecond

Default value: 1000.0

Barostat parameter: pressure relaxation time parameter of the barostat

barostat_tau <value>

Units: Femtosecond

Default value: 10000.0

Enable multistep MD

enable_multistep [true|false]

Default value: false

Number of steps between recomputation of slow terms in multistep MD

multistep_steps <number>

Default value: 1

Currently only exchange-repulsion EFP term is affected.

Periodic Boundary Conditions (PBC)

Enable/Disable PBC

enable_pbc [true|false]

Default value: false

Note

Setting enable_pbc to true also sets enable_cutoff to true.

Dimensions of periodic cell

Units: Angstroms, degrees.

periodic_box <x> <y> <z> <alpha> <beta> <gamma>

Default value: 30.0 30.0 30.0 90.0 90.0 90.0

Note

If only three values are given, the angles are set to 90 degrees (orthogonal box). Non-orthogonal PBC are implemented only for single-point energy calculations. The smallest box dimension must be greater than 2 * swf_cutoff.

Print PBC coordinates

Prints coordinates of the system contained in a single periodic cell around a fragment specified by ligand keyword.

print_pbc [true/false]

Default value: false

Pairwise interaction energy decomposition analysis

Enable/Disable pairwise analysis

The PIEDA analysis will be performed with respect to the fragment specified by the ligand keyword.

enable_pairwise [true|false]

Default value: false

Specify ligand

ligand [integer]

Default value: 0 (the first fragment in the system)

Symmetric molecular crystals

Use symmetry

If true, effectively performs calculations only on symmetry-unique fragments, which speeds up calculations of symmetric crystal systems with PBC. .. warning:: Implemented for single-point energy calculations. Not parallelized. See symm_frag keyword for specifying symmetry-identical fragments.

symmetry [true|false]

Default value: false

Specifying symmetry-identical fragments

symm_frag [frag | list]

  • frag - assumes that all fragments of the same type are identical.

  • list - fragments provided in a list are symmetric (not fully implemented).

Default value: frag

Calculations with NN potentials

Starting from version 2.0, calculations with NN potentials are enabled through interface with LibTorch library. Standard ANI potentials (TorchANI) and custom models with embedded electrostatic potentials (see enable_elpot) can be used. Single point energy and geometry optimization are currently implemented, see opt_special_frag for detail.

Enable NN calculations

enable_torch [true|false]

Default value: false

Specify NN fragment: no default value, need to specify the fragment to be treated with NN.

special_fragment <value>

Default value: none

Enable electrostatic potential embedding in NN calculations, i.e., NN/EFP

enable_elpot [true|false]

  • true - provide custom NN using custom_nn and aev_nn

  • false - provide one of ANI potentials with torch_nn

Default value: false

Optimization protocol in NN calculations

opt_special_frag [0|1]

  • 0 - optimize special (NN) fragment only

  • 1 - optimize the whole system

Default value: 0

Atomic gradients in NN calculations

atom_gradient [mm|frag]

  • mm - use direct atomic gradients from MM model

  • frag - gradient on special fragment atoms computed by distributing the COM force/torque gradient

Default value: frag

Note

mm algorithm can be used only with MM (qq and lj) potentials.

Note

frag algorithm distributes EFP center of mass force and torque into gradients on atoms. This decomposition is not unique resulting in non-exact gradients. If you experience convergence issues in geometry optimization, loosen optimization tolerance criteria opt_tol.

Name of the NN potential for standard ANI models

torch_nn <value>

Default value: ani.pt

Name of custom NN potential for NN/EFP embedded model

custom_nn <value>

Default value: custom_model_script.pt

Name of file containing atomic environment vector (AEV) for NN/EFP embedded model

aev_nn <value>

Default value: aev_scripted.pt

Path to the directory with provided NN potentials

ml_path <value>

Default value: “$(prefix)/share/libefp” (data install directory)

Path to the directory with user NN potentials

userml_path <value>

Default value: “.” (current directory)

Various other keywords

Print additional information for debugging

print <value>

  • 0 - standard output

  • 1 - detailed information on gradients

  • 2 - additional information on fragments; polarization convergence

  • 3 - various information for debugging

Default value: 0

Use single EFP parameters file

single_params_file [true|false]

Default value: false

Single EFP parameters file path

efp_params_file <path>

Default value: params.efp

The path to the directory with fragment library

fraglib_path <path>

Default value: “$(prefix)/share/libefp” (data install directory)

The <path> parameter should not contain spaces or should be inside double quotes otherwise.

The path to the directory with user-created fragments

userlib_path <path>

Default value: "." (current directory)

The <path> parameter should not contain spaces or should be in double quotes otherwise.

Parameters for test jobs

See also num_step_dist and num_step_angle in Hessian calculation section for controlling step size for numerical gradient testing.

Test tolerance

gtest_tol <value>

Default value: 1.0e-6

Unit: Hartree/Bohr

Reference energy value

ref_energy <value>

Default value: 0.0

Unit: Hartree

Depricated Parameters

Geometry of the molecular-mechanics part

ff_geometry <path>

Default value: ff.xyz

Molecular-mechanics force-field parameters file

ff_parameters <path>

Default value: fraglib/params/amber99.prm

Enable molecular-mechanics force-field for flexible EFP links

enable_ff [true|false]

Default value: false

Molecular geometry

Fragment name

LibEFP/EfpMD distribution comes with a library of fragment potentials used for testing and debugging purposes. Those fragments are stored in the fraglib directory.

Warning

.efp potentials from fraglib directory should not be used for any practical calculations as they might produce erroneous results; they are used for code testing only.

Note

In order to use the .efp potentials from fraglib directory (defined by fraglib_path parameter), the fragment <name> should contain an _l suffix, e.g., h2o_l, c6h6_l, etc.

User-defined parameters (or parameters from EFP parameter databases) do not need suffix _l, e.g., h2o, c6h6. Location of the parameter files is defined by userlib_path keyword, with the default being the working directory ./.

The name of the .efp file must be the same as the name of the fragment (without an _l suffix). For example for the fragment named H2O_L the parameters must be in the fraglib_path directory in the file named h2o.efp. For the fragment named NH3 the parameters must be in the userlib_path directory in the file named nh3.efp.

See $FRAGNAME for additional information.

Fragment input

Fragment input should contain one or more groups starting with fragment <name> line (see Fragment name) and followed by fragment coordinates in one of four possible formats controlled by coord keyword (see Basic parameters).

``xyzabc`` format

fragment h2o
0.0 0.0 0.0 0.0 0.0 0.0

The numbers are coordinates of the center of mass of a fragment in Angstroms and three Euler rotation angels in Radians.

``points`` format

fragment h2o
        0.0 0.0 0.0
        1.0 0.0 0.0
        0.0 1.0 0.0

Exactly three lines of three numbers is expected. These are x,y,z coordinates of the first three points (i.e., the first three atoms specified in the .efp potential) of a fragment in Angstroms.

``rotmat`` format

fragment h2o
        0.0 0.0 0.0
        1.0 0.0 0.0
        0.0 1.0 0.0
        0.0 0.0 1.0

The first line contain x,y,z coordinates of fragment’s center of mass in Angstroms. The next three lines contain a 3 x 3 rotation matrix determined with respect to the reference fragment stored in .efp file.

``atoms`` format

fragment nh3_l
A01N1                0.235191     4.772273    -0.103166
A02H2               -0.523189     4.322511     0.366654
A03H3                0.517310     5.557741     0.446020
A04H4                0.999952     4.130518    -0.141193

Atom name (as defined in .efp) and x,y,z coordinates in Angstroms for all fragment atoms should be given.

Additional fragment input

Fragment velocities

Each fragment input can be followed by fragment velocities (for restarting MD simulations) initiated with the velocity keyword, with the center of mass velocity and angular velocity in atomic units specified on the next line. An example showing the usage:

 1run_type md
 2coord xyzabc
 3ensemble nve
 4time_step 0.5
 5max_steps 50
 6fraglib_path ../fraglib
 7
 8fragment h2o_l
 9   0.0   0.0   0.0     0.0   0.0   0.0
10velocity
11   0.0   0.0   5.0e-4  0.0   0.0   0.0
12
13fragment nh3_l
14   0.0   0.0   5.0     0.0   0.0   0.0
15velocity
16   0.0   0.0  -7.0e-4  0.0   0.0   0.0

Fragment constraints

Quadratic constraint on the fragment center of mass can be specified using constraint keyword with the force constant k (in a.u.) and constraint position xyz (in Angstroms) specified on the next line. Example:

 1run_type opt
 2coord points
 3elec_damp screen
 4disp_damp tt
 5pol_damp tt
 6fraglib_path ../fraglib
 7
 8fragment h2o_l
 9  -3.394  -1.900  -3.700
10  -3.524  -1.089  -3.147
11  -2.544  -2.340  -3.445
12constraint
13  0.1  -3.0 -2.0 -3.0
14
15fragment nh3_l
16  -5.515   1.083   0.968
17  -5.161   0.130   0.813
18  -4.833   1.766   0.609
19
20fragment nh3_l
21   1.848   0.114   0.130
22   1.966   0.674  -0.726
23   0.909   0.273   0.517

Point charges

Additionally to fragments, a system can contain a set of point charges. They can be specified using the following format for each charge:

charge <q> <x> <y> <z>