The INPUT Files¶

Fig. 15 DL_POLY_5 input (left) and output (right) files. **Note**: files marked with an asterisk are non-mandatory.¶

DL_POLY_5 may require many input files. However, only CONTROL, CONFIG and FIELD are mandatory. The MPOLES and TAB* are complimentary to FIELD and are only required when opted within. HISTORY is required when an old trajectory is opted for re-play in CONTROL. REFERENCE is optionally required when defect detection is switched on in CONTROL. REVOLD is required only when the job represents a continuation of a previous job. In the following sections we describe the form and content of these files.

It is worth noting that historically DL_POLY used hard-coded names for different I/O files as shown in Figure (15). Upon instructions in the CONTROL many I/O file name can be overridden with specific, user-defined filenames (see ). Even the CONTROL file can be named differently but in this case the alternative name must be passed as a command line argument to the DL_POLY_5 executable (usually named DLPOLY.Z). Thus the DL_POLY_5 engine can be efficiently embedded and utilised within external frameworks for user-definable work-flows.

The CONTROL File¶

The CONTROL file is read by the subroutine read_new_control and defines the control variables for running a DL_POLY_5 job. Keywords are character strings that appear as the first entry on a data record (or line) and which invoke a particular operation or provide numerical parameters. Keywords can appear in any order in the CONTROL file and some of the directives are mandatory (for example the timestep directive that defines the timestep), others are optional.

This way of constructing the file is very convenient, but it has inherent dangers. It is, for example, quite easy to specify contradictory directives, or invoke algorithms that do not work together. By large DL_POLY_5 tries to sort out these difficulties and print helpful error messages, but it does not claim to be fully foolproof. It is important to think carefully about a simulation beforehand and ensure that DL_POLY_5 is being asked to do something that is physically reasonable. It should also be remembered that the present capabilities the package may not allow the simulation required and it may be necessary for you yourself to add new features.

An example CONTROL file appears below. The directives and keywords appearing are described in the following section.

title DL_POLY_5 CONTROL DIRECTIVES

# I/O REDIRECT
io_file_output my-dl_poly-run
io_file_field FIELD-dl_field.out
io_file_config CONFIG-CDS-generated
io_file_history new-trajectory
io_file_revive REVIVE-2ps-run
io_file_revcon REVCON-100k-output
io_file_revold REVIVE-1ps-run

# SYSTEM REPLICATION & IMPACT OPTION
nfold [ 10 10 10 ]
impact_part_index 1
impact_time  2000.0 steps
impact_energy  7.5 ke.V
impact_direction [  1.0  2.0  3.0 ]

# DENSITY VARIATION ARRAY BOOST
density_variance  10.0 %

# INDEX AND VERIFICATION BYPASS AND NO TOPOLOGY REPORTING
ignore_config_indices ON
strict_checks OFF
print_topology_info OFF

# INTERACTIONS BYPASS
vdw_method OFF
coul_method OFF

# APPLY MIXING TO ALLOWED & AVAILABLE VDW CROSS INTERACTIONS
vdw_mix_method Lorentz-Berthelot

# DIRECT CALCULATION OF VDW/METAL INTERACTIONS INSTEAD OF
# EVALUATION BY SPLINING OVER TABULATED VALUES IN MEMORY
vdw_method direct
metal_direct ON

# FORCE-SHIFT VDW INTERACTIONS SO THAT ENERGY AND FORCE
# CONTRIBUTIONS FALL SMOOTHLY TO ZERO WHEN APPROACHING R_CUT
vdw_force_shift ON

# RANDOM NUMBER GENERATOR SEEDING
random_seed [  100  200  300 ]

# RESTART OPTIONS
restart noscale
data_dump_frequency  1000 steps

# SYSTEM TARGET TEMPERATURE AND PRESSURE
pressure_hydrostatic  0.001 katm
temperature  300.0 K

# SYSTEM CUTOFFS AND ELECTROSTATICS
vdw_cutoff  8.0 ang
padding  0.35 ang
cutoff  10.0 ang
subcelling_threshold_density   50.0 %
coul_extended_exclusion ON
coul_dielectric_constant 1.0
coul_method spme
spme_precision 1e-05

# RELAXED SHELL MODEL TOLERANCE
rlx_tol  1.0

# CONSTRANTS ITERATION LENGTH and TOLERANCE
shake_max_iter 250
shake_tolerance  1e-05 ang

# INTEGRATION FLAVOUR, ENSEMBLE AND PSEUDO THERMOSTAT
ensemble nst
ensemble_method berendsen
ensemble_thermostat_coupling  0.5 ps
ensemble_barostat_coupling  1.5 ps
pseudo_thermostat_method langevin
pseudo_thermostat_width  2.0 ang
pseudo_thermostat_temperature  150.0 K

# INTEGRATION TIMESTEP
timestep  0.001 ps
timestep_variable ON
timestep_variable_min_dist  0.03 ang
timestep_variable_max_dist  0.1 ang
timestep_variable_max_delta  0.005 ps

# SIMULATION & EQUILIBRATION LENGTH
time_run  10000.0 steps
time_equilibration  1000.0 steps

# EQUILIBRATION DIRECTIVES
reset_temperature_interval  1.0 steps
equilibration_force_cap  500.0 k_B.temp/ang
rescale_frequency  5.0 steps
regauss_frequency 3.0 steps
minimisation_criterion energy
minimisation_tolerance  0.001 internal_e
minimisation_frequency  20.0 steps
minimisation_step_length  1e-05 ang

# STATISTICS
record_equilibration ON
stack_size  50.0 steps
stats_frequency  10.0 steps

# OUTPUT
print_frequency  20.0 steps

# HISTORY
traj_calculate ON
traj_start  20.0 steps
traj_interval  30.0 steps
traj_key pos

# DEFECTS TRAJECTORY - DEFECTS
defects_calculate ON
defects_start  40.0 steps
defects_interval  15.0 steps
defects_distance  0.75 ang

# DISPLACEMENTS TRAJECTORY - RSDDAT
displacements_calculate ON
displacements_start  70.0 steps
displacements_interval  10.0 steps
displacements_distance  0.25 ang

# MSDTMP
msd_calculate ON
msd_start  1000.0 steps
msd_frequency  100.0 steps

# INTRAMOLECULAR PDF ANALYSIS BY TYPE IF PRESENT
analyse_bonds ON
analyse_frequency_bonds 100 steps
analyse_num_bins_bonds 250

analyse_angles ON
analyse_frequency_angles 100 steps
analyse_num_bins_angles 360

analyse_dihedrals ON
analyse_frequency_dihedrals 100 steps
analyse_num_bins_dihedrals 720

analyse_inversions ON
analyse_frequency_inversions 100 steps
analyse_num_bins_inversions 360

analyse_all ON
analyse_frequency 100 steps
analyse_num_bins 1000
analyse_max_dist 5.0 ang


# RDF & Z-DENSITY
rdf_print ON
rdf_calculate ON
rdf_frequency  7.0 steps
rdf_binsize 0.05
zden_print ON
zden_calculate ON
zden_frequency  7.0 steps
zden_binsize 0.05

# EMPIRICAL VALENCE BOND (EVB) DIRECTIVE
evb 2

# EXECUTION TIME
time_job  1000.0 s
time_close  10.0 s

The CONTROL File Format¶

The file is free-formatted and not case-sensitive. Every line is treated as a command sentence (record). Commented records (beginning with a # or !) and blank lines are not processed and may be added to aid legibility (see example above). Records must be limited in length to 200 characters. Records are read in as up to three words “keyword value unit”. A word must not exceed 256 characters in length.

It is possible to split a record across multiple lines using “&”. See below for an example. Note any spacing must come on the new line

pressure_tensor [ 6.666 6.666 6.666&
 6.666 6.666 6.666 ] katm
correlation_observable [v_x&
-v_x&
 v_y-v_y]

Additional annotation must be rendered as a comment.

The first keyword (i.e. not comment) in the CONTROL file is must be title and a suiltable title.

The CONTROL File Directives¶

Note that in some cases additional keywords, shown in brackets “(…)”, may also be supplied in the directives, or directives may be used in a long form. However, it is strongly recommended that the user uses only the bold part of these directives.

The of directives available is as follows:

title	STRING	Run title
simulation_method	OPTION	Set simulation method, options: MD, EVB, FFS	md
random_seed	DATA_INT_VECTOR	Set random seed	1 2 3
density_variance	FLOAT	Set expected density variance for determining maximum array sizes	0.0	%
data_dump_frequency	FLOAT	Set data dumping frequency for restarts	1000	steps
subcell_threshold	FLOAT	Set subcelling threshold density for setting minimum particles per link-cell	50.0
evb_num_ff	INT	Set number of forcefields to be used in EVB	1
time_run	FLOAT	Set duration of simulation (inc. equilibration)	0	steps
time_equilibration	FLOAT	Set equilibration duration	0	steps
time_job	FLOAT	Set total job time before attempted safe closure	-1.0	hr
time_close	FLOAT	Estimated closure time for finite-time jobs	-1.0	min
stats_frequency	FLOAT	Set frequency of stats sampling to statis file	0	steps
stack_size	FLOAT	Set rolling average stack to n timesteps	0	steps
record_equilibration	BOOL	Include equilibration in output	off
print_probability_distribution	BOOL	Calculate and print probability distribution (enforces RDF print)	off
analyse_all	BOOL	Enable analysis for all bonds, angles, dihedrals and inversions	off
analyse_angles	BOOL	Enable analysis for all angles	off
analyse_bonds	BOOL	Enable analysis for all bonds	off
analyse_dihedrals	BOOL	Enable analysis for all dihedrals	off
analyse_inversions	BOOL	Enable analysis for all inversions	off
analyse_frequency	FLOAT	Set global frequency of data analysis	1	steps
analyse_frequency_bonds	FLOAT	Set frequency of bonds data analysis	1	steps
analyse_frequency_angles	FLOAT	Set frequency of angles data analysis	1	steps
analyse_frequency_dihedrals	FLOAT	Set frequency of dihedrals data analysis	1	steps
analyse_frequency_inversions	FLOAT	Set frequency of inversions data analysis	1	steps
analyse_max_dist	FLOAT	Set cutoff for bonds analysis	2.0	ang
analyse_num_bins	INT	Set global number of bins to be used in bonding analysis	-1
analyse_num_bins_bonds	INT	Set number of bins to be used in bond analysis	0
analyse_num_bins_angles	INT	Set number of bins to be used in angle analysis	0
analyse_num_bins_dihedrals	INT	Set number of bins to be used in dihedral analysis	0
analyse_num_bins_inversions	INT	Set number of bins to be used in inversion analysis	0
msd_calculate	BOOL	Enable calculation of MSD	off
msd_print	BOOL	Enable printing of MSD	off
msd_start	FLOAT	Start timestep for dumping MSD configurations	0	steps
msd_frequency	FLOAT	Interval between dumping MSD configurations	1	steps
correlation_observable	DATA_STRING_VECTOR	Enable calculation of correlations
correlation_blocks	DATA_INT_VECTOR	Correlation blocks
correlation_block_points	DATA_INT_VECTOR	Correlation blocks
correlation_window	DATA_INT_VECTOR	correlation window averaging
correlation_update_frequency	DATA_INT_VECTOR	correlation update frequency		steps
correlation_dump_frequency	FLOAT	frequency to dump correlation data	0	steps
traj_calculate	BOOL	Enable calculation of trajectory	off
traj_key	OPTION	Set trajectory output, options: pos, pos-vel, pos-vel-force, compressed	pos
traj_start	FLOAT	Start timestep for dumping trajectory configurations	0	steps
traj_interval	FLOAT	Interval between dumping trajectory configurations	1	steps
defects_calculate	BOOL	Enable calculation of defects	off
defects_start	FLOAT	Start timestep for dumping defects configurations	0	steps
defects_interval	FLOAT	Interval between dumping defects configurations	1	steps
defects_distance	FLOAT	Set cutoff for deviation to be considered by defects as interstitial	0.75	ang
defects_backup	BOOL	Enable defects backup	off
displacements_calculate	BOOL	Enable calculation of displacements	off
displacements_start	FLOAT	Start timestep for dumping displacements configurations	0	steps
displacements_interval	FLOAT	Interval between dumping displacements configurations	1	steps
displacements_distance	FLOAT	Set cutoff for qualifying as displacement	0.75	ang
coord_calculate	BOOL	Enable calculation of coordination	off
coord_ops	OPTION	Set Coordops, options: icoord: only dumps the coordination of each atom at i; CCOORD: only dumps coordination of each atom at i; FULL: dumps the coordination of each atom every j steps	icoord
coord_start	FLOAT	Start timestep for dumping coordination configurations	0	steps
coord_interval	FLOAT	Interval between dumping coordination configurations	100	steps
adf_calculate	BOOL	Enable calculation of ADF	off
adf_frequency	FLOAT	Set frequency of ADF sampling	100	steps
adf_precision	FLOAT	Set precision of angular distribution bins in ADF analysis	0.0
rdf_calculate	BOOL	Enable calculation of RDF	off
rdf_print	BOOL	Enable printing of RDF	on
rdf_frequency	FLOAT	Set frequency of RDF sampling	1	steps
rdf_start	FLOAT	Set the start of RDF sampling	0	steps
rdf_binsize	FLOAT	Set number of bins to be used in RDF analysis	0.05	ang
rdf_error_analysis	OPTION	Enable RDF error analysis, options: Off, Jackknife, Block	off
rdf_error_analysis_blocks	INT	Set number of RDF error analysis blocks	1
zden_calculate	BOOL	Enable calculation of ZDen	off
zden_print	BOOL	Enable printing of ZDen	on
zden_frequency	FLOAT	Set frequency of ZDen sampling	1	steps
zden_binsize	FLOAT	Set number of bins to be used in ZDen analysis	0.05	ang
vaf_calculate	BOOL	Enable calculation of VAF	off
vaf_print	BOOL	Enable printing of VAF	on
vaf_frequency	FLOAT	Set frequency of VAF sampling	1	steps
vaf_binsize	INT	Set number of bins to be used in VAF analysis	0
vaf_averaging	BOOL	Ignore time-averaging of VAF, report all calculated VAF to VAFDAT files and final profile to OUTPUT	on
currents_calculate	BOOL	Enable calculation of currents	off
energy_stress_currents	BOOL	Enable calculation of energy currents and k-dependent stress	off
heat_flux	BOOL	Enable calculation of heat flux	off
momentum_density	DATA_STRING_VECTOR	Enable calculation of momentum_density
write_per_particle	BOOL	Enable dumping of per-particle information	off
elastic_constants	BOOL	Enable calculation of elastic constants	off
print_frequency	FLOAT	Set frequency of printing results to output	0	steps
io_units_scheme	OPTION	Set I/O units scheme, options: internal, si, atomic, hartree, kj, kcal, boltzmann, dpd	internal
io_units_length	OPTION	Set I/O units for length (unused)	internal_l
io_units_time	OPTION	Set I/O units for time (unused)	internal_t
io_units_mass	OPTION	Set I/O units for mass (unused)	internal_m
io_units_charge	OPTION	Set I/O units for charge (unused)	internal_q
io_units_energy	OPTION	Set I/O units for energy (unused)	internal_e
io_units_pressure	OPTION	Set I/O units for pressure (unused)	internal_p
io_units_force	OPTION	Set I/O units for force (unused)	internal_f
io_units_velocity	OPTION	Set I/O units for velocity (unused)	internal_v
io_units_power	OPTION	Set I/O units for power (unused)	internal_e/internal_t
io_units_surface_tension	OPTION	Set I/O units for surface tension (unused)	internal_f/internal_l
io_units_emf	OPTION	Set I/O units for electromotive force (unused)	internal_e/internal_q
io_read_method	OPTION	Set I/O read method, possible read methods: mpiio, direct, master	mpiio
io_read_readers	INT	Set number of parallel I/O readers	0	(Automatic)
io_read_batch_size	INT	Set I/O reader batch size	0	(Automatic)
io_read_buffer_size	INT	Set I/O reader buffer size	0	(Automatic)
io_read_error_check	BOOL	Enable extended error checking on read	off
io_read_ascii_revold	BOOL	Read human-readable (ASCII) REVOLD file	off
io_write_method	OPTION	Set I/O write method, possible write methods: mpiio, direct, master	mpiio
io_write_writers	INT	Set number of parallel I/O writers	0	(Automatic)
io_write_batch_size	INT	Set I/O writer batch size	0	(Automatic)
io_write_buffer_size	INT	Set I/O writer buffer size	0	(Automatic)
io_write_sorted	BOOL	Enable sorted output for atomic data	on
io_write_error_check	BOOL	Enable extended error checking on write	off
io_write_ascii_revive	BOOL	Write REVIVE as a human-readable (ASCII) file	off
io_file_output	STRING	Set output filepath, special options: SCREEN, NONE	OUTPUT
io_file_config	STRING	Set input configuration filepath	CONFIG
io_file_field	STRING	Set input field filepath	FIELD
io_file_field_2	STRING	Set input field filepath for evb second state	FIELD2
io_file_field_3	STRING	Set input field filepath for evb third state	FIELD3
io_file_statis	STRING	Set output statistics filepath	STATIS
io_file_heatflux	STRING	Set output heatflux filepath	HEATFLUX
io_file_history	STRING	Set output history filepath	HISTORY
io_file_historf	STRING	Set output historf filepath	HISTORF
io_file_revive	STRING	Set output revive filepath, special options: NONE	REVIVE
io_file_revold	STRING	Set output revold filepath	REVOLD
io_file_revcon	STRING	Set output revcon filepath, special options: NONE	REVCON
io_file_rdf	STRING	Set output RDF filepath, special options: NONE	RDFDAT
io_file_msd	STRING	Set output MSD filepath, special options: NONE	MSDTMP
io_file_currents	STRING	Set output CURRENTS filepath	CURRENTS
io_file_tabbnd	STRING	Set input TABBND filepath	TABBND
io_file_tabang	STRING	Set input TABANG filepath	TABANG
io_file_tabdih	STRING	Set input TABDIH filepath	TABDIH
io_file_tabinv	STRING	Set input TABINV filepath	TABINV
io_file_tabvdw	STRING	Set input TABVDW filepath	TABVDW
io_file_tabeam	STRING	Set input TABEAM filepath	TABEAM
io_file_cor	STRING	Set output COR filepath, special options: NONE	COR
io_file_setevb	STRING	Set filepath for SETEVB	SETEVB
io_file_popevb	STRING	Set filepath for POPEVB	POPEVB
io_statis_yaml	BOOL	write statis file in yaml format	off
io_rdf_yaml	BOOL	write rdf file in yaml format	off
io_file_revcon_2	STRING	Set output revcon filepath, second evb	REVCON2
io_file_revcon_3	STRING	Set output revcon filepath, third evb	REVCON3
io_file_config_2	STRING	Set input configuration filepath for second evb	CONFIG2
io_file_config_3	STRING	Set input configuration filepath for third EVB	CONFIG3
output_energy	BOOL	Output final energy e_tot in output file	off
output_std_dev	BOOL	Output rolling standard deviations in OUTPUT	off
ignore_config_indices	BOOL	Ignore indices as defined in CONFIG and use read order instead	off
print_topology_info	BOOL	Print topology information in output file	off
print_level	INT	Disable unnecessary printing, levels: 0 - silent, 1 - quiet, 2 - standard, 3 - full	1
timer_depth	INT	Do not display timers beyond this many levels in final timer output	4
timer_yaml_file	BOOL	Output timers as a yaml file	off
timer_per_mpi	BOOL	Write timings for each MPI process individually	off
timestep	FLOAT	Set calculation timestep or initial timestep for variable timestep calculations	0.0	internal_t
timestep_variable	BOOL	Enable variable timestep	off
timestep_variable_min_dist	FLOAT	Set minimum permissible distance for variable timestep	0.03	ang
timestep_variable_max_dist	FLOAT	Set maximum permissible distance for variable timestep	0.1	ang
timestep_variable_max_delta	FLOAT	Set maximum timestep delta for variable timestep	0.0	internal_t
reference_scaling_matrix	DATA_FLOAT_VECTOR	Set the reference simulation cell for strain calculations (column major) defaults to initial cell.		ang
ensemble	OPTION	Set ensemble constraints, options: NVE, PMF, NVT, NPT, NST	NVE
ensemble_method	OPTION	Set ensemble method, options NVT: Evans, Langevin, Andersen, Berendsen, Hoover, gentle, ttm, dpds1, dpds2. NP\|ST: Langevin, Berendsen, Hoover, MTK.
ensemble_thermostat_coupling	FLOAT	Set thermostat relaxation/decorrelation times (use ensemble_thermostat_friction for langevin)	0.0	ps
ensemble_dpd_order	OPTION	Set dpd method, options: off, first, second	off
ensemble_dpd_drag	FLOAT	Set DPD drag coefficient	0.0	Da/ps
ensemble_thermostat_friction	FLOAT	Set thermostat friction for langevin and gentle stochastic thermostats	0.0	ps^-1
ensemble_thermostat_softness	FLOAT	Set thermostat softness for Andersen thermostat	0.0
ensemble_barostat_coupling	FLOAT	Set barostat relaxation/decorrelation times (use ensemble_barostat_friction for langevin)	0.0	ps
ensemble_barostat_friction	FLOAT	Set barostat friction	0.0	ps^-1
ensemble_semi_isotropic	OPTION	Enable semi-isotropic barostat constraints, options: area, tension, orthorhombic	off
ensemble_semi_orthorhombic	BOOL	Enable semi-orthorhombic barostat constraints	off
ensemble_tension	FLOAT	Set tension in NPngT calctulation	0.0	N/m
pressure_tensor	DATA_FLOAT_VECTOR	Set the target pressure tensor for NsT calculations	0.0 0.0 0.0 0.0 0.0 0.0	katm
pressure_hydrostatic	FLOAT	Set the target hydrostatic pressure (1/3Tr[P]) for NPT calculations	0.0	katm
pressure_perpendicular	DATA_FLOAT_VECTOR	Set the target pressure as x, y, z perpendicular to cell faces for NPT calculations	0.0 0.0 0.0	katm
temperature	FLOAT	Set the initial temperature or target temperature (for thermostats)	0.0	K
pseudo_thermostat_method	OPTION	Set pseudo thermostat method, possible options: Off, Langevin-Direct, Langevin, Gauss, Direct	off
pseudo_thermostat_width	FLOAT	Set the width of thermostatted boundaries for pseudo thermostats	2.0	ang
pseudo_thermostat_temperature	FLOAT	Set the temperature of the pseudo thermostat	0.0	K
impact_part_index	INT	Set particle index for impact simulations	0
impact_time	FLOAT	Set time for impact in impact simulations	0.0	internal_t
impact_energy	FLOAT	Set impact energy for impact simulations	0.0	ke.V
impact_direction	DATA_FLOAT_VECTOR	Direction vector for impact simulations	1.0 1.0 1.0
ttm_calculate	BOOL	Enable calculation of two-temperature model	off
ttm_num_ion_cells	INT	Set number of coarse-grained ion temperature cells (CIT)	10
ttm_num_elec_cells	DATA_FLOAT_VECTOR	Set number of coarse-grained electronic temperature cells (CET)	50 50 50
ttm_metal	BOOL	Specifies parameters for metallic system are required for two-temperature model, i.e. thermal conductivity	off
ttm_heat_cap_model	OPTION	Sets model for specific heat capacity in TTM, options: const, linear, tabulated, tanh
ttm_heat_cap	FLOAT	Sets constant, scale or maximum (per atom) specific heat capacity in TTM	0.0	k_B
ttm_temp_term	FLOAT	Set Fermi temperature in TTM, for tanh	0.0	K^-1
ttm_fermi_temp	FLOAT	Set Fermi temperature in TTM, for linear	0.0	K
ttm_elec_cond_model	OPTION	Set electronic conductivity model in TTM, options: Infinite, constant, drude, tabulated
ttm_elec_cond	FLOAT	Set electronic conductivity in TTM	0.0	W/m/K
ttm_diff_model	OPTION	Set diffusion model in TTM, options: constant, recip, tabulated
ttm_diff	FLOAT	Set TTM thermal diffusivity	0.0	m^2/s
ttm_dens_model	OPTION	Set density model in TTM, options are: constant, dynamic
ttm_dens	FLOAT	Set constant density in TTM	0.0	ang^-3
ttm_min_atoms	INT	Minimum number of atoms needed per ionic temperature cell	0
ttm_stopping_power	FLOAT	Electronic stopping power of projectile entering electronic system	0.0	e.V/nm
ttm_spatial_dist	OPTION	Set the spatial distribution of TTM, options: flat, gaussian, flat-laser, exp-laser
ttm_spatial_sigma	FLOAT	Set the sigma for spatial distributions of TTM	1.0	nm
ttm_spatial_cutoff	FLOAT	Set the cutoff for spatial distributions of TTM as a multiple of sigma	5.0
ttm_fluence	FLOAT	Initial energy deposition into electronic system by laser for TTM	0.0	mJ/cm^2
ttm_penetration_depth	FLOAT	Set laser penetration depth for TTM	0.0	nm
ttm_laser_type	OPTION	Set laser deposition type. options: flat, exponential	flat
ttm_temporal_dist	OPTION	Set temporal distribution for TTM, options: gaussian, exponential, delta, square
ttm_temporal_duration	FLOAT	Set duration of energy deposition for TTM (gaussian, exponential, square)	0.001	ps
ttm_temporal_cutoff	FLOAT	Set temporal cutoff for TTM (gaussian, exponential) as a multiple of its duration	5.0
ttm_variable_ep	OPTION	Set electron-phonon coupling for TTM, options: homo, hetero
ttm_boundary_condition	OPTION	Set boundary conditions for TTM, options: periodic, dirichlet, neumann, robin
ttm_boundary_xy	BOOL	Fix Neumann (zero-flux) boundary in Z	off
ttm_boundary_heat_flux	FLOAT	Set boundary heat flux in Robin boundaries for TTM	96	%
ttm_time_offset	FLOAT	Set electron-ion coupling offset for TTM	0.0	ps
ttm_oneway	BOOL	Enable one-way electron-phonon coupling when electronic temperature is greater than ionic temperature	off
ttm_stats_frequency	FLOAT	Frequency of write to TTM PEAK_E and PEAK_I	0	steps
ttm_traj_frequency	FLOAT	Frequency of write to TTM LATS_E and LATS_I	0	steps
ttm_com_correction	OPTION	Apply inhomogeneous Langevin thermostat to selected directions in TTM, options: full, zdir, off	full
ttm_redistribute	BOOL	Redistribute electronic energy in newly-deactivated temperature cells to nearest active neighbours
ttm_e-phonon_friction	FLOAT	Set TTM electron-phonon friction	0.0	ps^-1
ttm_e-stopping_friction	FLOAT	Set TTM electron-stopping friction	0.0	ps^-1
ttm_e-stopping_velocity	FLOAT	Set TTM electron-stopping velocity	0.0	ang/ps
ttm_e-phonon_cutoff_velocity	FLOAT	Set TTM electron-phonon coupling cutoff velocity	0.0	ang/ps
rlx_cgm_step	FLOAT	Set CGM stepping for relaxed shell model	-1.0	ang
rlx_tol	FLOAT	Set force tolerance for relaxed shell model	1.0	internal_f
shake_max_iter	INT	Set maximum number of SHAKE/RATTLE iterations	250
shake_tolerance	FLOAT	Set accepted SHAKE/RATTLE tolerance	1e-6	ang
dftb	BOOL	Enable DFTB	off
fixed_com	BOOL	Remove net centre of mass momentum	on
reset_temperature_interval	FLOAT	Interval between temperature zeroing during equilibration for minimisation	-1	steps
regauss_frequency	FLOAT	Set the frequency of temperature regaussing	-1	steps
rescale_frequency	FLOAT	Set the frequency of temperature rescaling	-1	steps
temperature_increment_frequency	FLOAT	Set the frequency of temperature incrementing	-1	steps
temperature_increment_start	FLOAT	Set the start of temperature incrementing	0	steps
temperature_increment_stop	FLOAT	Set the target heating/cooling temperature	0	K
temperature_increment	FLOAT	Set the (positive) temperature increment	0	K
equilibration_force_cap	FLOAT	Set force cap clamping maximum force during equilibration	1000.0	k_b.temp/ang
minimisation_criterion	OPTION	Set minimisation criterion, options: off, force, energy, distance	off
minimisation_tolerance	FLOAT	Set minimisation tolerance, units: determined by criterion	0.0
minimisation_step_length	FLOAT	Set minimisation tolerance	-1.0	ang
minimisation_frequency	FLOAT	Set minimisation frequency	0	steps
initial_minimum_separation	FLOAT	Turn on the check on minimum separation distance between VNL pairs at re/start	-1.0	internal_l
restart	OPTION	Set restart settings, possible options: Clean, Continue, Rescale, Noscale	clean
nfold	DATA_INT_VECTOR	Expand cell before running	1 1 1
cutoff	FLOAT	Set the global cutoff for real-speace potentials	1.0	internal_l
padding	FLOAT	Set padding for sizing of Verlet neighbour lists	0.0	internal_l
coul_damping	FLOAT	Calculate electrostatics using Fennell damping (Ewald-like) with given alpha	0.0	1/ang
coul_dielectric_constant	FLOAT	Set dielectric constant relative to vacuum	1.0
coul_bjerrum_length	FLOAT	Set Bjerrum length (use instead of dielectric constant relative to vacuum)	0.0	ang
coul_extended_exclusion	BOOL	Enable extended coulombic exclusion affecting intra-molecular interactions	off
coul_method	OPTION	Set method for electrostatics method, options: off, spme, dddp, pairwise, reaction_field, force_shifted	off
coul_precision	FLOAT	Calculate electrostatics using Fennell damping (Ewald-like) with given precision	0.0
charge_smearing	OPTION	Set method for charge smearing, options: off, linear, slater_approx, slater, gaussian, gaussian_equal	off
charge_smearing_length	FLOAT	Set length parameter for charge smearing	0.0	ang
charge_smearing_beta	OPTION	Set beta-lambda correspondence for Slater-type charge smearing,options: original, overlap, distribution	original
spme_precision	FLOAT	Set spme parameters to calculate within given precision for Spme calculations	1.0e-6
spme_alpha	FLOAT	Set real-recip changeover location for Spme calculations	0.0	ang^-1
spme_kvec	DATA_INT_VECTOR	Set number of k-space samples for SPME calculations	0 0 0
spme_kvec_spacing	FLOAT	Calculate k-vector samples for an even sampling of given spacing in spme calculations	0.0	ang^-1
spme_nsplines	INT	Set number of B-Splines for Ewald SPME calculations, min=3	8
polarisation_model	OPTION	Enable polarisation, options: default, CHARMM	default
polarisation_thole	FLOAT	Set global atomic damping factor	1.3
metal_direct	BOOL	Enable direct (non-tabulated) calculation of metallic forces	off
metal_sqrtrho	BOOL	Enable metal sqrtrho interpolation option for EAM embeding function in TABEAM	off
vdw_method	OPTION	Set method for Van der Waal’s calculations, options: off, direct, tabulated, spme	tabulated
vdw_cutoff	FLOAT	Set cut-off for Van der Waal’s potentials	-1.0	internal_l
vdw_mix_method	OPTION	Enable VdW mixing, possible mixing schemes: Off, Lorentz-Berthelot, Fender-Hasley, Hogervorst, Waldman-Hagler, Tang-Toennies, Functional	off
vdw_force_shift	BOOL	Enable force shift corrections to Van der Waals’ forces	off
plumed	BOOL	Enabled plumed dynamics	off
plumed_input	STRING	Set plumed input file
plumed_log	STRING	Set plumed log file
plumed_precision	INT	Set plumed numerical precision (4=single, 8=double)	8
plumed_restart	BOOL	Restart plumed dynamics	on
strict_checks	BOOL	Enforce strict checks such as: good system cutoff, particle index contiguity, disable non-error warnings, minimisation information	on
unsafe_comms	BOOL	Do not ensure checks of logicals are enforced in parallel	off
dftb_test	BOOL	Do not perform a DLPOLY run, instead run dftb tests	off
replay	BOOL	Don’t perform simulation, instead replay history	off
replay_calculate_forces	BOOL	On history replay recalculate forces	on

Further Comments on the CONTROL File¶

All keywords in new-style control files must be unique, double specification will result in an error. This is to protect users from unexpected behaviour.
The following parameters are mandatory:
1. ensemble
2. ensemble_method
3. cutoff
4. timestep
Some directives are optional. If not specified DL_POLY_5 may give default values if necessary. (defaults are specified above in the list of directives) However fail-safe DL_POLY_5 is, ensure parameters are appropriate for the system of interest and defaults should only be used after it is known that they are valid.
The time_run and time_equilibration directives have a default of zero. If not used or used with their default values a “dry run” is performed. This includes force generation and system dump (REVCON and REVIVE) and, depending on the rest of the options, may include; velocity generation, force capping, application of the CGM minimiser, application of the pseudo thermostat, and dumps of HISTORY, DEFECTS, RDFDAT, ZDNDAT and MSDTMP. Note that, since no actual dynamics is to be performed, the temperature and pressure directives do not play any role and are therefore not necessary.
If the CGM minimiser, minimise_frequency, is specified with zero frequency, it is only applied at timestep zero if $\textbf{time\_equilibration}~\ge~\textbf{time\_run}$ (i.e. optimise structure at start only!). In this way it can be used as a configuration optimiser at the beginning of the equilibration period or when a “dry run” (time_run$= 0$) is performed (i.e. equilibrate without any actual dynamics!). Note that the default CGM search algorithm stepping uses a step that is proportional to the square of the instantaneous value of the timestep and thus its usage in may lead to algorithm’s instability and failure in the cases when optimisation is applied before any dynamics occurs in a model system in an unphysical state (i.e. much away from equilibrium) and/or with an ill/badly defined forcefield for the state. Hence, special care (specifying the optional CGM stepping for example) should be taken when the option is used, especially in a “dry run” mode with a timestep value too large for the model system state.
The timestep_variable option requires the user to specify an initial guess for a reasonable timestep for the system. The simulation is unlikely to retain this as the operational timestep however, as the latter may change in response to the dynamics of the system. The option is used in conjunction with the default values of timestep_variable_max_dist ($0.10$ Å) and timestep_variable_min_dist ($0.03$ Å), which can also be optionally altered if used as directives (note the rule that timestep_variable_max_dist $>~2.5$timestep_variable_min_dist applies). Also, an additional timestep_variable_max_delta (in ps) control can be applied. These serve as control values in the variable timestep algorithm, which calculates the greatest distance a particle has travelled in any timestep during the simulation. If the maximum distance is exceeded, the timestep variable is halved and the step repeated. If the greatest move is less than the minimum allowed, the timestep variable is doubled and the step repeated provided it does not exceed the user specified timestep_variable_max_delta. If it does then it scales to timestep_variable_max_delta and the step is repeated. In this way the integration timestep self-adjusts in response to the dynamics of the system. Note that this option is abandoned when used in conjunction with DPD thermostats (ensemble_method dpd), since it cannot be applied efficiently and furthermore it is not well defined in a DPD sense either.
The starting options for a simulation are governed by the keyword restart. If this is not specified in the control file, the simulation will start as new. When specified, it will continue a previous simulation (restart) provided all needed restart files are in place and not corrupted. If they are not in place or are found corrupted, it will start a new simulation without initial temperature scaling of the previous configuration (restart noscale).
The ensemble nst keyword is also used in the N$\sigma$T ensembles extension to NP$_{n}$AT and NP$_{n}\gamma$T ones. Note that these semi-isotropic ensembles are only correct for infinite interfaces placed perpendicularly to the z axis! This means that the interface is homogeneous (unbroken) and continuous in the (x,y) plane of the MD cell, which assumes that that two of the cell vectors have a cross product only in the z direction. (For example, if the MD box is defined by its lattice vectors $(\underline{a},\underline{b},\underline{c})$ then $\underline{a} \times \underline{b} = \pm (0,0,1)$.) It is the users’ responsibility to ensure this holds for their model system.
The reset_temperature_interval directive, enables a “zero temperature” optimisation. In this case a crude energy minimiser is used to help the system relax before each integration of the equations of motion. The function of the minimiser can be summarised as

\[\begin{split}\underline{v}_{i} \leftarrow \left\{ \begin{array} {l@{\quad:\quad}l} 0 & \underline{v}_{i} \cdot \underline{f}_{i} < 0 \\ \underline{f}_{i} \; \frac{\underline{v}_{i}~\cdot~\underline{f}_{i}} {\underline{f}_{i}~\cdot~\underline{f}_{i}} & \underline{v}_{i} \cdot \underline{f}_{i} \ge 0 \end{array} \right.\end{split}\]

for systems with free particles only, where $\underline{v}_{i}$ and $\underline{f}_{i}$ are the force and velocity of particle $i$. The algorithm is extended in the case of RBs by including

\[\begin{split}\begin{aligned} \underline{V}_{j} \leftarrow \left\{ \begin{array} {l@{\quad:\quad}l} 0 & \underline{V}_{j} \cdot \underline{F}_{j} < 0 \\ \underline{F}_{j} \; \frac{\underline{V}_{j}~\cdot~\underline{F}_{j}} {\underline{F}_{j}~\cdot~\underline{F}_{j}} & \underline{V}_{j} \cdot \underline{F}_{j} \ge 0 \end{array} \right. \\ \nonumber \underline{\omega}_{j} \leftarrow \left\{ \begin{array} {l@{\quad:\quad}l} 0 & \underline{\omega}_{j} \cdot \underline{\tau}_{j} < 0 \\ \underline{\tau}_{j} \; \frac{\underline{\omega}_{j}~\cdot~\underline{\tau}_{j}} {\underline{\tau}_{j}~\cdot~\underline{\tau}_{j}} & \underline{\omega}_{j} \cdot \underline{\tau}_{j} \ge 0 \end{array} \right.~~,\end{aligned}\end{split}\]

where $\underline{V}_{j}$ and $\underline{F}_{j}$ are the velocity and force of the RB’s centre of mass, and $\underline{\omega}_{j}$ and $\underline{\tau}_{j}$ are the angular velocity and torque of RB $j$. Measures are taken to conserve the MD cell momentum and the thermostat’s instantaneous kinetic energy.

This must not be thought of as a true energy minimization method. Note that this optimisation is only applied when the simulation runs in equilibration mode.

The algorithm is developed in the DL_POLY_5 routine zero_k_optimise.
The impact directives will not be activated if the particle index is beyond the one of the last particle. The option will fail in a controlled manner at application time if the particle is found to be in a frozen state or the shell of an ion or part of a rigid body. During application the center of mass momentum is re-zeroed to prevent any drifts. The user must take care to have the impactinitiated after any possible equilibration. Otherwise, the system will be thermostatted and the impact energy dissipated during the equilibration.
The pseudo_thermostat directives are intended to be used in highly non-equilibrium simulations when users are primarily interested in the structural changes in the core of the simulated system as the the MD cell boundaries of the system are coupled to a thermal bath.

The thermal bath can be used with several types of temperature scaling algorithms
1. Langevin (stochastic thermostat),
2. Gauss
3. Direct (direct thermostat).
The user is also required to specify the width of the pseudo thermostat, pseudo_thermostat_width, which must be larger than 2 Å and less than or equal to a quarter of minimum width of the MD cell. The thermostat is an pseudo_thermostat_width thick buffer layer attached on the inside at the MD cell boundaries.

The temperature of the bath is specified by the user, pseudo_thermostat_temperature, which must be larger than 1 Kelvin.
- langevin
  
  The stochasticity of the Langevin thermostat emulates an infinite environment around the MD cell, providing a means for “natural” heat exchange between the MD system and the heath bath thus aiding possible heat build up in the system. In this way the instantaneous temperature of the system is driven naturally towards the bath temperature. Every particle within the thermostat buffer layer is coupled to a viscous background and a stochastic heat bath, such that
  
  \[\begin{split}\begin{aligned} {d \underline{r}_{i}(t) \over d t} =& \underline{v}_{i}(t) \nonumber \\ {d \underline{v}_{i}(t) \over d t} =& {{\underline{f}_{i}(t)+\underline{R}_{i}(t)} \over m_{i}} - \chi (t) \; \underline{v}_{i}(t)~~,\end{aligned}\end{split}\]
  
  where $\chi (t)$ is the friction parameter from the dynamics in the the MD cell and $R(t)$ is stochastic force with zero mean that satisfies the fluctuation-dissipation theorem:
  
  \[\left< R^{\alpha}_{i}(t)~R^{\beta}_{j}(t^\prime) \right> = 2~\chi(t)~m_{i}~k_{B}T~\delta_{ij}~\delta_{\alpha \beta}~\delta(t-t^\prime)~~,\]
  
  where superscripts denote Cartesian indices, subscripts particle indices, $k_{B}$ is the Boltzmann constant, $T$ the bath temperature and $m_{i}$ the particle’s mass. The algorithm is implemented in routine pseudo and has two stages:
  - Generate random forces on all particles within the thermostat. Here, care must be exercised to prevent introduction of non-zero net force when the random forces are added to the system force field.
  - Rescale the kinetic energy of the thermostat bath so that particles within have Gaussian distributed kinetic energy with respect to the target temperature and determine the (Gaussian constraint) friction within the thermostat:
    
    \[\chi(t) = Max \left( 0, \frac {\sum_{i} [\vec{f}_{i}(t) + \vec{R}_{i}(t)] \cdot \vec{v}_{i}(t)} {\sum_{i} m_{i}~\vec{v}_{i}^{2}(t)} \right)~~.\]
    
    Care must be exercised to prevent introduction of non-zero net momentum. (Users are reminded to use for target temperature the temperature at which the original system was equilibrated in order to avoid simulation instabilities.)
  The effect of this algorithm is to relax the buffer region of the system on a local scale and to effectively dissipate the incomingexcess kinetic energy from the rest of the system, thus emulating an infinite-like environment surrounding the MD cell. The thermostat width matters as the more violent the events on the inside of the MD cell, the bigger width may be needed in order to ensure safe dissipation of the excess kinetic energy.
- Gaussian
  - Rescale the kinetic energy of the thermostat bath so that particles within have Gaussian distributed kinetic energy with respect to the target temperature.
- direct
  
  The Direct thermostat is the simplest possible model allowing for heat exchange between the MD system and the heath bath. All (mass, non-frozen) particles within the bath have their kinetic energy scaled to $1.5~k_{B}T~$ at the end of each time step during the simulation. Care is exercised to prevent introduction of non-zero net momentum when scaling velocities. (Users are reminded to use for target temperature the temperature at which the original system was equilibrated in order to avoid simulation instabilities.) Due to the “unphysical” nature of this temperature control the thermostat width does not matter to the same extent as in the case of the Langevin thermostat.
Note that embedding a thermostat in the MD cell walls is bound to produce wrong ensemble averages, and instantaneous pressure and stress build-ups at the thermostat boundary. Therefore, ensembles lose their meaning as such and so does the conserved quantity for true ensembles.

The algorithms are developed in the DL_POLY_5 routine pseudo_vv.
The defects_calculate option will trigger reading of REFERENCE (see Section The REFERENCE File), which defines a reference MD cell with particles’ positions defining the crystalline lattice sites. If REFERENCE is not found the simulation will either:
- halt if the simulation has been restarted, i.e. is a continuation of an old one - the restart option is used in CONTROL and the REVOLD (see Section The REVOLD File) file has been provided.
- recover using CONFIG (see Section The CONFIG File) if it is a new simulation run, i.e restart option is not used in CONTROL or REVOLD has not been provided.
The actual defect detection is based on comparison of the simulated MD cell to the reference MD cell based on a user defined site-interstitial cutoff, $R_{def}$,

\[{\rm Min}\left[0.3,r_{\rm cut}/3\right]~\AA~\le~R_{def}~\le~{\rm Min}\left[1.2,r_{\rm cut}/2\right]~\AA \nonumber\]

with a default value of ${\rm Min}\left[0.75,r_{\rm cut}/3\right]$ Å. (If the supplied value exceeds the limits the simulation execution will halt). If a particle, $p$, is located in the vicinity of a site, $s$, defined by a sphere with its centre at this site and a radius, $R_{def}$, then the particle is a first hand claimee of $s$, and the site is not vacant. Otherwise, the site is presumed vacant and the particle is presumed a general interstitial. If a site, $s$, is claimed and another particle, $p\prime$, is located within the sphere around it, then $p\prime$ becomes an interstitial associated with $s$. After all particles and all sites are considered, it is clear which sites are vacancies. Finally, for every claimed site, distances between the site and its first hand claimee and interstitials are compared and the particle with the shortest one becomes the real claimee. If a first hand claimee of $s$ is not the real claimee it becomes an interstitial associated with $s$. At this stage it is clear which particles are interstitials. The sum of interstitials and vacancies gives the total number of defects in the simulated MD cell.

Frozen particles and particles detected to be shells of polarisable ions are not considered in the defect detection.

Note that the algorithm cannot be applied safely if $R_{def}$ is larger than half the shortest interatomic distance within the reference MD cell since a particle may; (i) claim more than one site, (ii) be an interstitial associated with more than one site, or both (i) and (ii). On the other hand, low values of $R_{def}$ are likely to lead to slight overestimation of defects.

If the simulation and reference MD cell have the same number of atoms then the total number of interstitials is always equal to the total number of defects.
The displacements_calculate option will trigger dump of atom displacements based on a qualifying cutoff in a trajectory like manner. Displacements of atoms from their original position at the end of equilibration (the start of statistics), $t~=0~$, is carried out at each timestep.
The tolerance and stepping for relaxed shell model rlx_tol, is a last resort option to aid shell relaxation of systems with very energetic and/or rough potential surface. Users are advised to use it with caution, should there really need be, as the use of high values for the tolerance (default of 1) may result in physically incorrect dynamics and small stepping (default $Max(k_{core-shell})/2$) to expensive force evaluations.
The dielectric constant relative to vacuum, $\epsilon$, is set using the coul_dielectric_constant directive and defaults to a value of 1, i.e. the value for vacuum. An alternative approach (and preferred for simulations specifying quantities in DPD units) is to specify a non-zero Bjerrum length $\lambda_B$ using the coul_bjerrum_length directive, which is used along with the system temperature to determine the dielectric constant [1].
The choice of reaction field electrostatics (directive coul_method reaction_field) relies on the specification of the relative dielectric constant external to the cavity. This is specified by either one of the coul_dielectric_constant or coul_bjerrum_length directives.
DL_POLY_5 uses two different potential cutoffs. These are as follows:
1. $r_{\rm cut}$ - the universal cutoff set by cutoff. It applies to the real space part of the electrostatics calculations and to the van der Waals potentials if no other cutoff is applied.
2. $r_{\rm vdw}$ - the user-specified cutoff for the van der Waals potentials set by vdw_cutoff. If not specified its value defaults to $r_{\rm cut}$.
Constraint algorithms in , SHAKE/RATTLE (see Section Bond Constraints), use default iteration precision of 10^-6 and limit of iteration cycles of 250. Users may experience that during optimisation of a new built system containing constraints simulation may fail prematurely since a constraint algorithm failed to converge. In such cases directives shake_max_iter (to increase) and shake_tolerance (to decrease) may be used to decrease the strain in the system and stablise the simulation numerics until equilibration is achieved.
DL_POLY_5’s DD strategy assumes that the local (per domain/node or link cell) density of various system entities (i.e. atoms, bonds, angles, etc.) does not vary much during a simulation and some limits for these are assumed empirically. This may not the case in extremely non-equilibrium simulations, where the assumed limits are prone to be exceeded or in some specific systems where these do not hold from the start. A way to tackle such circumstances and avoid simulations crash (by controlled termination) is to use the density_variance $f$ option. In the set_bounds subroutine DL_POLY_5 makes assumptions at the beginning of the simulation and corrects the lengths of bonded-like interaction lists arrays (mxshl, mxcons, mxrgd, mxteth, mxbond, mxangl, mxdihd, mxinv) as well as the lengths of link-cell (mxlist) and domain (mxatms, mxatdm) lists arrays when the option is activated with $f > 0$. Greater values of $f$ will correspond to allocation bigger global arrays and larger memory consumption by DL_POLY_5 during the simulation. Note that this option may demand more memory than available on the computer architecture. In such cases DL_POLY_5 will terminate with an array allocation failure message.
As a default, DL_POLY_5 does not store statistical data during the equilibration period. If the directive record_equilibration is used, equilibration data will be incorporated into the overall statistics.
The vaf_calculate directive switches on velocity autocorrelation function (VAF) calculations for individual atomic species in DL_POLY_5 after equilibration or immediately at start if the directive record_equilibration is used. It controls how often VAF profiles are started what the size of each profile (in timesteps). Overlapping profiles are possible and require more memory to store them (and initial velocities) while they are being calculated. By default DL_POLY_5 will report time-averaged VAF profiles. This can be overridden using the vaf_averaging directive, which will instead report individual ‘instantaneous’ VAF profiles.
The fixed_com option will trigger a default bypass of the getvom routine which will return zero and thus no COM removal will happen. Note that this will lead to COM momentum accumulation for many though not all ensembles!. Such accumulation will propagate to the generation of flow in the MD cell and ultimately suppress the thermal motion of the particles in the system, leading to the so called “frozen ice cube effect”! It is worth nothing that this option must be turned on for the correct application of stochastic dynamics via the langevin temperature control (NVT Langevin)! If the option is not applied then the dynamics will lead to peculiar thermalisation of different atomic species to mass- and system size-dependent temperatures.
The padding option will add extra distance, $r_{\rm pad}$, if larger than $f \ge {\rm Min}[0.05,0.5\%.r_{\rm cut}]$ Å, to the major cutoff, $r_{\rm cut}$, to construct a larger link-cell width, $r_{\rm lnk}=r_{\rm cut}+r_{\rm pad}$, which will trigger a construction of a larger Verlet neighbour list (VNL) while at the same time facilitate its conditional update, rather at every timestpe. The VNL conditaional update is check at the end of each tiemstep and triggered only when the most travelled particle has moved a distance larger than $r_{\rm pad}/2$. It is worth noting that padding is at expense of extra memory but if used wisely it could improve time to solution from 10% to 100% depending on force-field complexity. If it is too large or too small (that is the reason for the $f \ge {\rm Min}[0.05,0.5\%.r_{\rm cut}]$ Å limit) it will lead to performace degradation. It is recomended that $r_{\rm pad}$ is set up at a value of $\approx 1 \div 5\%$ of the cutoff, $r_{\rm cut}$, as long as the major link-cell algorithm uses a link-cell decomposition that not worse than $4 \otimes 4 \otimes 4$ per domain. For such setups, in practice, one may expect average compute speedups [2] of the order of $10 \div 30\%$ for force-fields involving the Ewald summation methodology and $60 \div 100\%$ for force-fields without electrostatics evaluations involving the Ewald summation methodology.
The subcelling_threshold_density option will set the threshold density of particles per link cell below which subcelling (decreasing link-cell size) stops. The default is 50 although experimenting with values between 10 and 100 may provide sufficient information about a performance sweetspot. The comparison condition is carried out on the average density per sub-link cell of each domain. Thus subcelling takes into account only the per domain density variation. Obviously, the subcelling threshhold is limited to 1 particle per cell for the sake of safety and sanity reasons.
The coul_extended_exclusion option will make sure that for all conventional (no distance restraints) intra-molecular interactions (bonds, angles, dihedrals, inversions) as well as for CB and RB units any intra-core-shell interactions fall within the list of excluded interactions. This is not a default behaviour. The option is also triggered by the polarisation directvies.
The ttm_calculate options has the effect of switching on the two-temperature model (TTM), which assumes the use of the inhomogeneous Langevin NVT ensemble. It is not possible to specify any other ensemble for TTM-based systems, so the inhomogeneous Langevin NVT ensemble will always be used in this case, with default values for the friction terms and cut-off velocity if these are not specified. No implementation of the inhomogeneous Langevin ensemble currently exists for rigid bodies. To model metallic systems, the thermal conductivity must be specified either as infinitely large, a constant value, a linear function of electronic temperature (based on the Drude model) or as a tabulated function of temperature, while non-metallic systems require a constant thermal diffusivity. Any energy deposition is applied after equilibration. An additional restart file, DUMP_E, consisting of electronic temperatures for each electronic temperature cell (CET), is produced along with REVCON and REVIVE when TTM is switched on.

Users are advised to study the example CONTROL files appearing in the data sub-directory to see how different files are constructed.

The CONFIG File¶

The CONFIG file contains the dimensions of the unit cell, the key for periodic boundary conditions and the atomic labels, coordinates, velocities and forces. This file is read by the subroutine read_config (optionally by scan_config) in the set_bounds routine. The first few records of a typical CONFIG file are shown below:

IceI structure 6x6x6 unit cells with proton disorder
         2         3               276
  26.988000000000000   0.000000000000000   0.000000000000000
 -13.494000000000000  23.372293600000000   0.000000000000000
   0.000000000000000   0.000000000000000  44.028000000000000
      OW         1
    -2.505228382        -1.484234330        -7.274585343
    0.5446573999        -1.872177437       -0.7702718106
     3515.939287         13070.74357         4432.030587
      HW         2
    -1.622622646        -1.972916834        -7.340573742
     1.507099154        -1.577400769         4.328786484
     7455.527553        -4806.880540        -1255.814536
      HW         3
    -3.258494716        -2.125627191        -7.491549620
     2.413871957        -4.336956694         2.951142896
    -7896.278327        -8318.045939        -2379.766752
      OW         4
    0.9720599243E-01    -2.503798635        -3.732081894
     1.787340483        -1.021777575        0.5473436377
     9226.455153         9445.662860         5365.202509

etc.

The CONFIG File Format¶

The file is free-formatted and case sensitive for the atom species names. Every line is treated as a command sentence (record). However, line records are limited to 72 characters in length. Records are read in words, as a word must not exceed 40 characters in length. Words are recognised as such by separation by one or more space characters. The first record in the CONFIG file is a header (up to 72 characters long) to aid identification of the file. Blank and commented lines are not allowed.

Definitions of Variables in the CONFIG File¶

record 1
header      a72       title line

record 2
levcfg      integer   CONFIG file key. See Table (6) for permitted values
imcon       integer   Periodic boundary key. See Table (7) for permitted values
megatm      integer   Total number of particles (crystalographic entities)

record 3 omitted if imcon = 0
cell(1)     real      x component of the a cell vector in Å (or DPD length units)
cell(2)     real      y component of the a cell vector in Å (or DPD length units)
cell(3)     real      z component of the a cell vector in Å (or DPD length units)

record 4 omitted if imcon = 0
cell(4)     real      x component of the b cell vector in Å (or DPD length units)
cell(5)     real      y component of the b cell vector in Å (or DPD length units)
cell(6)     real      z component of the b cell vector in Å (or DPD length units)

record 5 omitted if imcon = 0
cell(7)     real      x component of the c cell vector in Å (or DPD length units)
cell(8)     real      y component of the c cell vector in Å (or DPD length units)
cell(9)     real      z component of the c cell vector in Å (or DPD length units)

Note that record 2 may contain more information apart from the mandatory as listed above. If the file has been produced by DL_POLY_5 then it also contains other items intended to help possible parallel I/O reading. Also, it is worth mentioning that the periodic boundary conditions (PBC), as specified in Table (7) and described in detail in Appendix B, refer generally to a triclinic type of super-cell, for which there are no symmetry assumptions! Records 3, 4 and 5 contain the Cartesian components of the super-cell’s lattice vectors in Å. DL_POLY_5 can only tract triclinic type of super-cells as the only types of super-cell shapes that are commensurate with the domain decomposition (DD) parallelisation strategy of it. However, this is not a restriction for the replicated data (RD) parallelisation that adopts and thus it can also accept truncated octahedral and rhombic dodecahedral periodic boundaries.

Subsequent records consists of blocks of between 2 and 4 records depending on the value of the levcfg variable. Each block refers to one atom. The atoms do not need to be listed sequentially in order of increasing index. Within each block the data are as follows:

record i
atmnam   a8       atom name
index    integer  atom index

record ii
xxx      real     x coordinate in Å (or DPD length units)
yyy      real     y coordinate in Å (or DPD length units)
zzz      real     z coordinate in Å (or DPD length units)

record iii included only if levcfg \(>0\)
vxx      real     x component of velocity in Å/picosecond (or DPD velocity units)
vyy      real     y component of velocity in Å/picosecond (or DPD velocity units)
vzz      real     x component of velocity in Å/picosecond (or DPD velocity units)

record iv included only if levcfg \(>1\)
fxx      real     x component of force in Å\(\cdot\)Dalton/picosecond² (or DPD force units)
fyy      real     y component of force in Å\(\cdot\)Dalton/picosecond² (or DPD force units)
fzz      real     z component of force in Å\(\cdot\)Dalton/picosecond² (or DPD force units)

Note that on record i only the atom name is strictly mandatory, any other items are not read by DL_POLY_CLASSIC but may be added to aid alternative uses of the file, for example alike DL_POLY GUI ¹¹⁸, DL_POLY_CLASSIC assume that the atoms’ indices are in a sequentially ascending order starting form 1. However, needs the index or the no index option needs to be specified in the CONTROL file! It is worth mentioning that DL_POLY_5 (as well as DL_POLY_CLASSIC) assumes that the origin of Cartesian system with respect to which the particle positions are specified is the middle of MD cell. Also, as both the cell vectors and the particles’ positions are specified in Å, there is a fine connection between them! This would not be the case if the particles’ positions were kept in reduced space with fractional coordinates. Last but not least, it is worth pointing out that composite entities, such as velocities and forces, have their units expressed as composites of the default DL_POLY units as shown in Section [units].

Table 6 CONFIG File Key (record 2)¶
`levcfg`
0	coordinates included in file
1	coordinates and velocities included in file
2	coordinates, velocities and forces included in file

Table 7 Periodic Boundary Key (record 2)¶
`imcon`
0	no periodic boundaries
1	cubic boundary conditions
2	orthorhombic boundary conditions
3	parallelepiped boundary conditions
6	x-y parallelogram boundary conditions with
	no periodicity in the z direction

Further Comments on the CONFIG File¶

The CONFIG file has the same format as the output file REVCON (Section The REVCON File). When restarting from a previous run of DL_POLY_5 (i.e. using the restart, restart noscale or restart scale directives in the CONTROL file - above), the CONFIG file must be replaced by the REVCON file, which is renamed as the CONFIG file. The copy macro in the execute sub-directory of DL_POLY_5 does this for you.

The CONFIG file has the same format as the optional output file CFGMIN, which is only produced when the minimise (optimise) option has been used during an equilibration simulation or a “dry run”.

The FIELD File¶

The FIELD file contains the force field information defining the nature of the molecular forces. This information explicitly includes the (site) topology of the system which sequence must be matched (implicitly) in the crystallographic description of the system in the CONFIG file. The FIELD file is read by the subroutine read_field. (It is also read by the subroutine scan_field in the set_bounds routine.) Excerpts from a force field file are shown below. The example is the antibiotic Valinomycin in a cluster of 146 water molecules:

Valinomycin Molecule with 146 SPC Waters
UNITS kcal

MOLECULES    2
Valinomycin
NUMMOLS 1
ATOMS 168
    O        16.0000     -0.4160         1
    OS       16.0000     -0.4550         1
    "          "           "        "
    "          "           "        "
    HC        1.0080      0.0580         1
    C        12.0100      0.4770         1
BONDS 78
harm   31   19 674.000     1.44900
harm   33   31 620.000     1.52600
    "    "    "    "         "
    "    "    "    "         "
harm  168   19 980.000     1.33500
harm  168  162 634.000     1.52200
CONSTRAINTS 90
   20   19    1.000017
   22   21    1.000032
    "    "     "
    "    "     "
  166  164    1.000087
  167  164    0.999968
ANGLES 312
harm   43    2   44  200.00      116.40
harm   69    5   70  200.00      116.40
  "    "    "    "     "           "
  "    "    "    "     "           "
harm   18  168  162  160.00      120.40
harm   19  168  162  140.00      116.60
DIHEDRALS 371
harm    1   43    2   44  2.3000      180.00
harm   31   43    2   44  2.3000      180.00
  "    "    "    "    "     "           "
  "    "    "    "    "     "           "
cos   149   17  161   16  10.500      180.00
cos   162   19  168   18  10.500      180.00
FINISH
SPC Water
NUMMOLS 146
ATOMS  3
    OW       16.0000     -0.8200
    HW        1.0080      0.4100
    HW        1.0080      0.4100
CONSTRAINTS  3
    1    2    1.0000
    1    3    1.0000
    2    3    1.63299
FINISH
VDW   45
C       C         lj   0.12000      3.2963
C       CT        lj   0.08485      3.2518
"       "         "     "            "
"       "         "     "            "
"       "         "     "            "
OW      OS        lj   0.15100      3.0451
OS      OS        lj   0.15000      2.9400
CLOSE

The FIELD File Format¶

The file is free-formatted and not case-sensitive (except for the site names). Every line is treated as a command sentence (record). Commented records (beginning with a #) and blank lines are not processed and may be added to aid legibility (see example above). Records must be limited in length to 200 characters. Records are read in words, as a word must not exceed 40 characters in length. Words are recognised as such by separation by one or more space characters. The contents of the file are variable and are defined by the use of directives. Additional information is associated with the directives.

Definitions of Variables in the FIELD File¶

The file divides into three sections: general information, molecular descriptions, and non-bonded interaction descriptions, appearing in that order in the file.

General information¶

The first viable record in the FIELD file is the title. The second is the units directive. Both of these are mandatory.

record 1
header      a200    field file header

record 2
units       a40     Unit of energy used for input and output

The energy units on the units directive are described by additional keywords:

eV, for electron-Volts
kcal/mol, for k-calories per mol
kJ/mol, for k-Joules per mol
Kelvin/Boltzmann, for Kelvin per Boltzmann
dpd, for DPD-based units (see DPD units)
internal, for DL_POLY internal units (10 Joules per mol).

If no units keyword is entered, DL_POLY internal units are assumed for both input and output. The units directive only affects the input and output interfaces, all internal calculations are handled using DL_POLY units. System input and output energies are read in units per MD cell.

Note that all energy bearing potential parameters are read in terms of the specified energy units. If such a parameter depends on an angle then the dependence is read in terms of radians although the following angle in the parameter sequence is read in terms of degrees. As to any rule, there is an exception - the electrostatic bond potential has no energy bearing parameter!

A third optional record can then be supplied before any molecular description:

record 3
multipolar order \(n\)  a50, integer    Electrostatics
evaluation to order \(n\) poles

This will later trigger the parsing of the MPOLES file (see Section The MPOLES File) which supplies the multipolar momenta values of the molecular sites specified in FIELD. Sites with specified pole orders, $m$, smaller than the one required, $n$, will have their $m+1$ to $n$ order poles’ momenta zeroed. Similarly, if sites have momenta of poles of higher order than the one required, $n$, these will not be processed by .

Note

Although algorithms in DL_POLY_5 could in principle handle any high pole order summation, in practice, however, DL_POLY_5 will abort if the order is higher than hexadecapole (order 4)! For more information on this functionality refer to Section The MPOLES File.

Molecular details¶

It is important for the user to understand that there is an organisational correspondence between the FIELD file and the CONFIG file described above. It is required that the order of specification of molecular types and their atomic constituents in the FIELD file follows the order of indices in which they appear in the CONFIG file. Failure to adhere to this common sequence will be detected by DL_POLY_5 and result in premature termination of the job. It is therefore essential to work from the CONFIG file when constructing the FIELD file. It is not as difficult as it sounds!

The entry of the molecular details begins with the mandatory directive:

molecules n

where n is an integer specifying the number of different types of molecule appearing in the FIELD file. Once this directive has been encountered, DL_POLY_5 enters the molecular description environment in which only molecular description keywords and data are valid.

Immediately following the molecules directive, are the records defining individual molecules:

name-of-molecule which can be any character string up to 200 characters in length. (Note: this is not a directive, just a simple character string.)
nummols n where n is the number of times a molecule of this type appears in the simulated system. The molecular data then follow in subsequent records:

atoms n where n indicates the number of atoms in this type of molecule. A number of records follow, each giving details of the atoms in the molecule i.e. site names, masses and charges. Each record carries the entries:

sitnam  a8      atomic site name
weight  real    atomic site mass (in Daltons or DPD mass units)
chge    real    atomic site charge (in protons)
nrept   integer repeat counter
ifrz    integer ‘frozen’ atom (if ifrz > 0`)

The integer nrept need not be specified if the atom/site is not frozen (in which case a value of 1 is assumed.) A number greater than 1 specified here indicates that the next (nrept-1) entries in the CONFIG file are ascribed the atomic characteristics given in the current record. The sum of the repeat numbers for all atoms in a molecule should equal the number specified by the atoms directive.

shell n where n is the number of core-shell units. Each of the subsequent n records contains:
```
index 1 (i)     integer   site index of core
index 2 (j)     integer   site index of shell
k_2             real      force constant of core-shell spring
k_4             real      quartic (anharmonic) force constant of spring
```
The spring potential is

\[U(r)=\frac{1}{2}k_{2} r_{ij}^{2}+\frac{1}{4!}k_{4} r_{ij}^{4}~~,\]

with the force constant $k_{2}$ entered in units of engunit$\times$Å^-2 and $k_{4}$ in engunit Å^-4, where usually $k_{2} >> k_{4}$. The engunit is the energy unit specified in the units directive. (In the case of DPD units, substitute Å with the DPD length unit.)

Note that the atomic site indices referred to above are indices arising from numbering each atom in the molecule from 1 to the number specified in the atoms directive for this molecule. This same numbering scheme should be used for all descriptions of this molecule, including the constraints, pmf, rigid, teth, bonds, angles, dihedrals and inversions entries described below. DL_POLY_5 will itself construct the global indices for all atoms in the systems.

Note that DL_POLY_5 determines which shell model to use by scanning shells’ weights provided the FIELD file. If all shells have zero weight the DL_POLY_5 will choose the relaxed shell model. If no shell has zero weight then DL_POLY_5 will choose the dynamical one. In case when some shells are massless and some are not DL_POLY_5 will terminate execution controllably and provide information about the error and possible possible choices of action in the OUTPUT file (see Section The OUTPUT Files). Shell models’ extensions are dealt according to the user specifications in the FIELD files.

This directive (and associated data records) need not be specified if the molecule contains no core-shell units.

constraints n where n is the number of constraint bonds in the molecule. Each of the following n records contains:

index 1     integer   first atomic site index
index 2     integer   second atomic site index
bondlength  real      constraint bond length

This directive (and associated data records) need not be specified if the molecule contains no constraint bonds. See the note on the atomic indices appearing under the shell directive above.

pmf b where b is the potential of mean force bondlength (Å). There follows the definitions of two PMF units:
1. pmf unit n1 where n1 is the number of sites in the first unit. The subsequent n1 records provide the site indices and weighting. Each record contains:
```
index     integer   atomic site index
weight    real      site weighting
```
2. pmf unit n2 where n2 is the number of sites in the second unit. The subsequent n2 records provide the site indices and weighting. Each record contains:
```
index     integer     atomic site index
weight    real        site weighting
```
This directive (and associated data records) need not be specified if no PMF constraints are present. See the note on the atomic indices appearing under the shell directive.

Note that if a site weighting is not supplied DL_POLY_5 will assume it is zero. However, DL_POLY_5 detects that all sites in a PMF unit have zero weighting then the PMF unit sites will be assigned the masses of the original atomic sites.

The PMF bondlength applies to the distance between the centres of the two PMF units. The centre, $\vec{R}_{i}$, of each unit is given by

\[\underline{R}_{i}~=~\frac{\sum_{j=1}^{n_{i}} w_{j}~\vec{r}_{j}}{\sum_{j=1}^{n_{j}} w_{j}}~~,\]

where $r_{j}$ is a site position and $w_{j}$ the site weighting.

Note that the PMF constraint is intramolecular. To define a constraint between two molecules, the molecules must be described as part of the same DL_POLY_5 “molecule”. DL_POLY_5 allows only one type of PMF constraint per system. The value of nummols for this molecule determines the number of PMF constraint in the system.

Note that in DL_POLY_5 PMF constraints are handled in every available ensemble.
rigid n where n is the number of basic rigid units in the molecule. It is followed by at least n records, each specifying the sites in a rigid unit:
```
m         integer     number of sites in rigid unit
site 1    integer     first site atomic index
site 2    integer     second site atomic index
site 3    integer     third site atomic index

etc.
site m integer m’th site atomic index
```
Up to 15 sites can be specified on the first record. Additional records can be used if necessary. Up to 16 sites are specified per record thereafter.

This directive (and associated data records) need not be specified if the molecule contains no rigid units. See the note on the atomic indices appearing under the shell directive above.

teth n where n is the number of tethered atoms in the molecule. It is followed n records specifying the tehered sites in the molecule:

tether key      a4        potential key, see Table (8)
index 1 (i)     integer   atomic site index
variable 1      real      potential parameter, see Table (8)
variable 2      real      potential parameter, see Table (8)

The meaning of these variables is given in Table (8).

This directive (and associated data records) need not be specified if the molecule contains no flexible chemical bonds. See the note on the atomic indices appearing under the shell directive above.

Table 8 Tether Potentials¶
Key	Potential Type	Variables	Functional Form
harm	Harmonic	$k$	$U(r) = \frac{1}{2}~k~(r_{i}-r_{i}^{t=0})^{2}$
rhrm	Restraint	$k$, $r_{c}$	$U(r) = \frac{1}{2}~k~(r_{i}-r_{i}^{t=0})^{2}$ : $\|r_{i}-r_{i}^{t=0}\| \le r_{c}$ $U(r) = \frac{1}{2}~k~r_{c}^{2}+k~r_{c}(\|r_{i}-r_{i}^{t=0}\|-r_{c})~:~\|r_{i}-r_{i}^{t=0}\|>r_{c}$
quar	Quartic	$k$, $k'$, $k''$	$U(r) = \frac{k}{2}~(r_{i}-r_{i}^{t=0})^{2}+\frac{k'}{3}~(r_{i}-r_{i}^{t=0})^{3}$ $+\frac{k''}{4}~(r_{i}-r_{i}^{t=0})^{4}$

bonds n where n is the number of flexible chemical bonds in the molecule. Each of the subsequent n records contains:

Bond key      a4        potential key, see Table (9)
index 1 (i)   integer   first atomic site index in bond
index 2 (j)   integer   second atomic site index in bond
variable 1    real      potential parameter, see Table (9)
variable 2    real      potential parameter, see Table (9)
variable 3    real      potential parameter, see Table (9)
variable 4    real      potential parameter, see Table (9)

The meaning of these variables is given in Table (9).

Table 9 Bond Potentials¶
Key	Potential Type	Variables	Functional Form
tab-tab	Tabulation		See tabulated potential Section The TABBND, TABANG, TABDIH & TABINV Files in TABBND file.
harm-hrm	Harmonic	$k$, $r_{0}$	$U(r) = \frac{1}{2}~k~(r_{ij}-r_{0})^{2}$
mors-mrs	Morse	$E_{0}$, $r_{0}$, $k$	$U(r) = E_{0}~[\{1-\exp(-k~(r_{ij}-r_{0}))\}^{2}-1]$
12-6 - 126	12-6	$A$, $B$	$U(r) = \left (\frac{A}{r_{ij}^{12}}\right)-\left(\frac{B}{r_{ij}^{6}}\right)$
lj	Lennard-Jones	$\epsilon$, $\sigma$	$U(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r_{ij}} \right)^{12} - \left( \frac{\sigma}{r_{ij}}\right)^{6} \right]$
rhrm-rhm	Restraint	$k$, $r_{0}$, $r_{c}$	$U(r) = \frac{1}{2}~k~(r_{ij}-r_{0})^{2}~:~\|r_{ij}-r_{0}\| \le r_{c}$ $U(r) = \frac{1}{2}~k~r_{c}^{2}+k~r_{c}(\|r_{ij}-r_{0}\|-r_{c})~:~\|r_{ij}-r_{0}\|>r_{c}$
quar-qur	Quartic	$k$, $r_{0}$, $k'$	$U(r) = \frac{k}{2}~(r_{ij}-r_{0})^{2}+\frac{k'}{3}~(r_{ij}-r_{0})^{3}$ $+\frac{k''}{4}~(r_{ij}-r_{0})^{4}$
buck-bck	Buckingham	$A$, $\rho$, $C$	$U(r) = A~\exp\left(-\frac{r_{ij}}{\rho}\right)-\frac{C}{r_{ij}^{6}}$
coul-cul	Coulomb	$k$	$U(r) = k \cdot U^{Electrostatics}(r_{ij}) \; \left(= \frac{k}{4\pi\epsilon_{0}\epsilon}\frac{q_{i}q_{j}}{r_{ij}}\right)$
fene-fne	Shifted $^{*}$ FENE ^12,44,143	$k$, $R_{o}$, $\Delta$	$U(r) = -0.5~k~R_{o}~ln\left[1-\left(\frac{r_{ij}-\Delta}{R_{o}^{2}}\right)^{2}\right]~:~r_{ij} < R_{o} + \Delta$ $U(r) = \infty~:~r_{ij} \ge R_{o} + \Delta$
mmst-mst	MM3 bond stretch ⁵	$k$, $r_{o}$	$U(r) = k~\delta^{2}\left[1-2.55~\delta+(7/12)~2.55^{2}~\delta^{2}\right]~;~\delta~ = r-r_{o}$

$^{*}$ Note: $\Delta$ defaults to zero if $|\Delta| > 0.5~R_{o}$ or if it is not specified in the FIELD file.

Note: Bond potentials with a dash (-) as the first character of the keyword, do not contribute to the excluded atoms list (see Section Force Fields). In this case DL_POLY_5 will also calculate the non-bonded pair potentials between the described atoms, unless these are deactivated by another potential specification.

angles n where n is the number of valence angle bonds in the molecule. Each of the n records following contains:

Angle key     a4        potential key, see Table (10)
index 1 (i)   integer   first atomic site index
index 2 (j)   integer   second atomic site index (central site)
index 3 (k)   integer   third atomic site index
variable 1    real      potential parameter, see Table (10)
variable 2    real      potential parameter, see Table (10)
variable 3    real      potential parameter, see Table (10)
variable 4    real      potential parameter, see Table (10)

The meaning of these variables is given in Table (10).

This directive (and associated data records) need not be specified if the molecule contains no angular terms. See the note on the atomic indices appearing under the shell directive above.

Table 10 Bond Angle Potentials¶
Key	Potential Type	Variables	Functional Form
tab-tab	Tabulation		see tabulated potential Section The TABBND, TABANG, TABDIH & TABINV Files in TABANG file
harm-hrm	Harmonic	$k$	$\theta_{0}$ $U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}$
quar-qur	Quartic	$k$, $\theta_{0}$, $k'$, $k''$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2} + {k'\over 3}(\theta - \theta_{0})^{3} + {k''\over 4}(\theta - \theta_{0})^{4}$
thrm-thm	Truncated harmonic	$k$, $\theta_{0}$, $\rho$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2} \exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]$
shrm-shm	Screened harmonic	$k$, $\theta_{0}$, $\rho_{1}$, $\rho_{2}$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]$
bvs1-bv1	Screened Vessal ¹⁴¹	$k$, $\theta_{0}$, $\rho_{1}$, $\rho_{2}$	$U(\theta) = {k \over 8(\theta_{0}-\pi)^{2}} \left\{ \left[ (\theta_{0} - \pi)^{2} -(\theta-\pi)^{2}\right]^{2}\right\} \times$ $\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]$
bvs2-bv2	Truncated Vessal ¹²³	$k$, $\theta_{0}$, $a$, $\rho$	$U(\theta) = k~(\theta-\theta_{0})^{2}~\big[ \theta^a (\theta+\theta_{0}-2\pi)^{2}$ $+ {a \over 2} \pi^{a-1} (\theta_{0}-\pi)^{3}\big] \exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]$
hcos-hcs	Harmonic Cosine	$k$, $\theta_{0}$	$U(\theta) = {k \over 2}~(\cos(\theta) - \cos(\theta_{0}))^{2}$
cos-cos	Cosine	$A$, $\delta$, $m$	$U(\theta) = A~[1+\cos(m~\theta-\delta)]$
mmsb-msb	MM3 stretch-bend ⁵	$A$, $\theta_{0}$, $r_{ij}^{o}$, $r_{jk}^{o}$	$U(\theta) = A~(\theta-\theta_{0})~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o})$
stst-sts	Compass stretch-stretch ¹²⁶	$A$, $r_{ij}^{o}$, $r_{jk}^{o}$	$U(\theta) = A~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o})$
stbe-stb	Compass stretch-bend ¹²⁶	$A$, $\theta_{0}$, $r_{ij}^{o}$	$U(\theta) = A~(\theta-\theta_{0})~(r_{ij} - r_{ij}^{o})$
cmps-cmp	Compass all terms ¹²⁶	$A$, $B$, $C$, $\theta_{0}$, $r_{ij}^{o}$, $r_{jk}^{o}$	$U(\theta) = A~(r_{ij} - r_{ij}^{o})~(r_{ik} - r_{ik}^{o}) + (\theta-\theta_{0}) \times$ $[B~(r_{ij} - r_{ij}^{o})+C~(r_{ik} - r_{ik}^{o})]$
mmbd-mbd	MM3 angle bend ⁵	$k$, $\theta_{0}$	$U(\theta) = k~\Delta^{2}[1 - 1.4 \cdot 10^{-2}\Delta + 5.6 \cdot 10^{-5}\Delta^{2}$ $- 7.0 \cdot 10^{-7}\Delta^{3}+2.2 \cdot 10^{-8}\Delta^{4})]~;~\Delta = \theta-\theta_{0}$
kky-kky	KKY ⁸⁹	$f_{k}$, $\theta_{0}$, $g_{r}$, $r_{o}$	$U(\theta) = 2~f_{k}~\sqrt{K_{ij} \cdot K_{ik}}~\sin^{2}\left[(\theta-\theta_{0})\right];$ $K_{ij} = 1/\left[\exp\left[g_{r}(r_{ij}-r_{o})\right]+1\right]$

$\theta$ is the $i$-$j$-$k$ angle.

Note: valence angle potentials with a dash (-) as the first character of the keyword, do not contribute to the excluded atoms list (see Section Force Fields). In this case DL_POLY_5 will calculate the non-bonded pair potentials between the described atoms.

dihedrals n where n is the number of dihedral interactions present in the molecule. Each of the following n records contains:

Dihedral key    a4        potential key, see Table (11)
index 1 (i)     integer   first atomic site index
index 2 (j)     integer   second atomic site index (central site)
index 3 (k)     integer   third atomic site index
index 4 (l)     integer   fourth atomic site index
variable 1      real      first potential parameter, see  Table (11)
variable 2      real      second potential parameter, see  Table (11)
variable 3      real      third potential parameter, see  Table (11)
variable 4      real      1-4 electrostatic interaction scale factor
variable 5      real      1-4 van der Waals interaction scale factor
variable 6      real      fourth potential parameter, see Table (11)
variable 7      real      fifth potential parameter, see Table (11)

The meaning of the variables 1-3,6-7 is given in Table (11). The variables 4 and 5 specify the scaling factor for the 1-4 electrostatic and van der Waals non-bonded interactions respectively.

This directive (and associated data records) need not be specified if the molecule contains no dihedral angle terms. See the note on the atomic indices appearing under the shell directive above.

Table 11 Dihedral Potentials¶
Key	Potential Type	Variables	Functional Form $^\ddagger$
tab	Tabulation		see tabulated potential Section The TABBND, TABANG, TABDIH & TABINV Files in TABDIH file
cos	Cosine	$A$, $\delta$, $m$	$U(\phi) = A~\left[ 1 + \cos(m\phi - \delta) \right]$
harm	Harmonic	$k$, $\phi_{0}$	$U(\phi) = {k \over 2}~(\phi - \phi_{0})^{2}$
hcos	Harmonic cosine	$k$, $\phi_{0}$	$U(\phi) = {k \over 2}~(\cos(\phi) - \cos(\phi_{0}))^{2}$
cos3	Triple cosine	$A_{1}$, $A_{2}$, $A_{3}$	$U(\phi) = {1 \over 2}~\{ A_{1}~(1+\cos(\phi))~+$ $A_{2}~(1-\cos(2\phi))~+$ $A_{3}~(1+\cos(3\phi)) \}$
ryck	Ryckaert-Bellemans ¹⁰¹	$A$	$U(\phi) = A~\{a+b~\cos(\phi)+c~\cos^{2}(\phi)~+$ $d~\cos^{3}(\phi)+e~\cos^{4}(\phi)+f~\cos^{5}(\phi) \}$
rbf	Fluorinated Ryckaert-Bellemans ¹⁰⁷	$A$	$U(\phi) = A~\{ a+b~\cos(\phi)+c~\cos^{2}(\phi)~+$ $d~\cos^{3}(\phi)+e~\cos^{4}(\phi)+f~\cos^{5}(\phi))~+$ $g~\exp(-h(\phi-\pi)^{2})) \}$
opls	OPLS torsion	$A_{0}$, $A_{1}$, $A_{2}$, $A_{3}$, $\phi_{0}$	$U(\phi) = A_{0} + \frac{1}{2}~\{ A_{1}~(1+\cos(\phi-\phi_{0}))~+$ $A_{2}~(1-\cos(2(\phi-\phi_{0})))~+$ $A_{3}~(1+\cos(3(\phi-\phi_{0}))) \}$

$^\ddagger$ $\phi$ is the $i$-$j$-$k$-$l$ dihedral angle.

inversions n where n is the number of inversion interactions present in the molecule. Each of the following n records contains:

Inversion key     a4        potential key, see Table (12)
index 1 (i)       integer   first atomic site index (central site)
index 2 (j)       integer   second atomic site index
index 3 (k)       integer   third atomic site index
index 4 (l)       integer   fourth atomic site index
variable 1        real      potential parameter, see Table (12)
variable 2        real      potential parameter, see Table (12)
variable 3        real      potential parameter, see Table (12)

The meaning of the variables 1-2 is given in Table (12).

This directive (and associated data records) need not be specified if the molecule contains no inversion angle terms. See the note on the atomic indices appearing under the shell directive above.

Table 12 Inversion Angle Potentials¶
Key	Potential Type	Variables	Functional Form $^\ddagger$
tab	Tabulation		see tabulated potential Section The TABBND, TABANG, TABDIH & TABINV Files in TABINV file
harm	Harmonic	$k$, $\phi_{0}$	$U(\phi) = {k \over 2}~(\phi - \phi_{0})^{2}$
hcos	Harmonic cosine	$k$, $\phi_{0}$	$U(\phi)={k \over 2}~(\cos(\phi) - \cos(\phi_{0}))^{2}$
plan	Planar	$A$	$U(\phi) = A~[1 - \cos(\phi)]$
xpln	Extended planar	$k$, $m$, $\phi_{0}$	$U(\phi) = {k \over 2}~[ 1 - \cos(m~\phi - \phi_{0}) ]$
calc	Calcite [calcite]	$A$, $B$	$U(u)= Au^{2}+Bu^{4}$

$^\ddagger$ $\phi$ is the $i$-$j$-$k$-$l$ inversion angle.

Note that the calcite potential is not dependent on an angle $\phi$, but on a displacement $u$. See Section The Calcite Four-Body Potential for details.

finish This directive is entered to signal to DL_POLY_5 that the entry of the details of a molecule has been completed.

The entries for a second molecule may now be entered, beginning with the name-of-molecule record and ending with the finish directive.

The cycle is repeated until all the types of molecules indicated by the molecules directive have been entered.

The user is recommended to look at the example FIELD files in the data directory to see how typical FIELD files are constructed.

Non-bonded Interactions¶

Non-bonded interactions are identified by atom types as opposed to specific atomic indices. The following different types of non-bonded potentials are available in ; vdw - van der Waals pair, metal - metal, tersoff - Tersoff, tbp - three-body and fbp - four-body. Each of these types is specified by a specific keyword as described bellow.

When DL_POLY_5 is cross-compiled with an OpenKIM functionality (see Section Open Knowledgebase of Interatomic Models - OpenKIM), it is possible to specify a complete model of inter-molecular interactions, by calling the OpenKIM model name provided it available in your local KIM library. Two keywords are employed when using OpenKIM IMs, one to select the IM and perform necessary initialisation (kim_init), and the other (kim_interactions) to define a mapping between atom types in DL_POLY_5 to the available species in the OpenKIM IM.

kim_init model_name where the model_name is the OpenKIM model identifier, it uses the information retrieved from the OpenKIM repository to initialise and activate OpenKIM IMs for use in the .
kim_interactions site_names where site_names defines a list of unique atomic names. It defines a mapping between atom types in DL_POLY_5 to the available species in the OpenKIM IM. For example, consider an OpenKIM IM that supports Si and C species. If the DL_POLY_5 simulation has four atoms, where the first three are Si, and the fourth is C, the kim_interactions would be used as:
```
kim_interactions Si C
```
The Si and C arguments map the DL_POLY_5 atom types to the Si and C species as defined within KIM PM.

In addition to the usual DL_POLY_5 error messages, the KIM library itself may generate errors, which should be printed to the screen. In this case, it is also useful to check the kim.log file for additional error information. The file kim.log should be generated in the same directory where DL_POLY_5 is running.

Note that although a KIM model fully describes a model system, it is still possible to specify further, complementing intra- and inter-molecular interactions in the FIELD! This is a viable option only when the model system is extended beyond what the specific KIM model is intended to describe.

By default, all the species in the DL_POLY_5 system should match to the available species in the OpenKIM IM.

Note that there is an experimental feature, implemented in the , where one can use a KIM IM in a hybrid style. In this case, one can create a system where part of species are interacting using a KIM IM (e.g., a machine learning model in KIM), and the rest of the species are interacting with the internal DL_POLY_5 intra- and inter-molecular interactions provided in the FIELD. While this new feature provides great flexibility, it is not fully compliant with the KIM API interface standard. Thus using this feature some of the KIM IMs might fail with the error message unexpected species code detected or unknown species detected.

vdw n where n is the number of pair potentials to be entered. It is followed by n records, each specifying a particular pair potential in the following manner:

atmnam 1      a8        first atom type
atmnam 2      a8        second atom type
key           a4        potential key, see Table (13)
variable 1    real      potential parameter, see Table (13)
variable 2    real      potential parameter, see Table (13)
variable 3    real      potential parameter, see Table (13)
variable 4    real      potential parameter, see Table (13)
variable 5    real      potential parameter, see Table (13)
variable 6    real      potential parameter, see Table (13)
variable 7    real      potential parameter, see Table (13)
variable 8    real      potential parameter, see Table (13)
variable 9    real      potential parameter, see Table (13)

The variables pertaining to each potential are described in Table (13). If any of the DPD thermostats are in use, two additional parameters for the atom pair’s dissipative force parameter $\gamma_{ij}$ and the thermostat cutoff distance $r_{t,ij}$ can be included after the parameters for the pair potential and used instead of any default values. If $\gamma_{ij}$ is not specified, DL_POLY_5 will either use mixing rules or a global value specified by ensemble_dpd_drag in the CONTROL file. If $r_{t,ij}$ is not supplied, DL_POLY_5 will use either the pair potential’s interaction cutoff distance (e.g. $r_c$ for standard DPD) or, failing that, the maximum van der Waals cutoff distance rvdw.

Note that any pair potential not specified in the FIELD file, will be assumed to be zero.

Table 13 Pair Potentials¶
Key	Potential Type	Variables	Functional Form
tab	Tabulation		see tabulated potential Sections Short Ranged (van der Waals) Potentials and The TABBND, TABANG, TABDIH & TABINV Files in TABVDW file
12-6	12-6	$A$, $B$	$U(r) = \left( \frac{A}{r^{12}} \right) - \left( \frac{B}{r^{6}} \right)$
lj	Lennard-Jones	$\epsilon$, $\sigma$	$U(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right]$
ljc	LJ cohesive ¹⁰	$\epsilon$, $\sigma$, $c$	$U(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - c~\left( \frac{\sigma}{r}\right)^{6} \right]$
ljf	LJ Frenkel ¹⁴²	$\epsilon$, $\sigma$, $r_c$	see Eq. (53)
nm	n-m ^18,78	$E_{o}$, $n$, $m$, $r_{0}$	$U(r) = \frac{E_{o}}{(n-m)} \left[ m \left( \frac{r_{o}}{r} \right)^{n}-n \left( \frac{r_{o}}{r} \right)^{m} \right]$
buck	Buckingham	$A$, $\rho$, $C$	$U(r) = A~\exp \left( -\frac{r}{\rho} \right) - \frac{C}{r^{6}}$
bhm	Born-Huggins-Meyer	$A$, $B$, $\sigma$, $C$, $D$	$U(r)=A~\exp[B(\sigma-r)]-\frac{C}{r^{6}}-\frac{D}{r^{8}}$
hbnd	12-10 H-bond	$A$, $B$	$U(r) = \left( \frac{A}{r^{12}} \right) - \left( \frac{B}{r^{10}} \right)$
snm	Shifted force$^{\dagger}$ n-m ^18,78	$E_{o}$, $n$, $m$, $r_{0}$, $r_{c}{}^{\ddagger}$	$U(r) = \frac{\alpha E_{o}}{(n-m)}\times$ $r_{c}$$^{\ddagger}$ $\left[ m\beta^{n} \left\{ \left( \frac{r_{o}}{r} \right)^{n} - \left( \frac{1}{\gamma} \right)^{n} \right\} - n\beta^{m} \left\{\left( \frac{r_{o}}{r} \right)^{m} - \left( \frac{1}{\gamma}\right)^{m}\right\} \right]$ $+ \frac{nm\alpha E_{o}}{(n-m)} \left( \frac{r-\gamma r_{o}}{\gamma r_{o}} \right) \left\{ \left( \frac{\beta}{\gamma}\right)^{n} - \left( \frac{\beta}{\gamma} \right)^{m} \right\}$
mors	Morse	$E_{0}$, $r_{0}$, $k$	$U(r) = E_{0}[\{1-\exp(-k(r-r_{0}))\}^{2}-1]$
wca	Shifted$^{*}$ ¹⁴⁸ Weeks-Chandler-Andersen	$\epsilon$, $\sigma^{\ddagger}$, $\Delta$	$U(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r-\Delta} \right)^{12} - \left( \frac{\sigma}{r-\Delta}\right)^{6} \right] + \epsilon : r_{ij} < 2^{1 \over 6}~\sigma + \Delta$ $U(r) = 0 : r_{ij} \ge 2^{1 \over 6}~\sigma + \Delta$
dpd	Standard DPD ⁴⁷ (Groot-Warren) $^{\dagger}$	$A$, $r_{c}{}^{\ddagger}$	$U(r) = \frac{A}{2}~r_{c}~\left(1-\frac{r}{r_{c}}\right)^{2}~:~r < r_{c}$ $U(r) = 0~:~r \ge r_{c}$
ndpd	$n$DPD ¹²⁴	$A$, $b$, $n$, $r_{c}{}^{\ddagger}$	$U(r) = \frac{Ab}{n+1}~r_{c}~\left(1-\frac{r}{r_{c}}\right)^{n+1}-\frac{A}{2}~r_{c}~~\left(1-\frac{r}{r_{c}}\right)^{2}~:~r < r_{c}$ $U(r) = 0~~~~~~~~~~~~~~~~~~~:~r \ge r_{c}$
mdpd	Generalised many-body DPD ¹⁴⁰	$A$, $B$, $n$, $m$, $r_{c}{}^{\ddagger}$, $r_{d}$	See Dissipative Particle Dynamics (DPD)
14-7	14-7 buffered AMOEBA FF ⁸⁹	$\epsilon$, $r_{o}$	$U(r) = \epsilon\left(\frac{1.07}{{(r_{ij}/r_{o})+0.07}^{\vphantom{7}}}\right)^{7}\left(\frac{1.12}{{(r_{ij}/r_{o})^{7}+0.12}^{\vphantom{7}}}-2\right)$
mstw	Morse modified	$E_{0}$, $r_{0}$, $k$, $c$	$U(r) = E_{0}\{\left[1-\exp\left(-k\left(r-r_{0}\right)\right)\right]^{2}-1\}+ \frac{c}{r^{12}}$
ryd	Rydberg	$a$, $b$, $\rho$	$U(r) = (a+br)\exp\left(-r/\rho\right)$
zbl	ZBL	$Z_{1}$, $Z_{2}$	$U(r) = \frac{Z_{1}~Z_{2}~e^{2}}{4\pi\varepsilon_0\varepsilon_r}\sum_{i=1}^{4}b_{i} \exp\left(-c_{i}r/a\right)$
zbls	ZBL mixed with Morse	$Z_{1}$, $Z_{2}$, $r_{m}$, $\xi$, $E_{0}$, $r_{0}$, $k$	$U(r) = f\left(r\right)U_{ZBL}\left(r\right)~~+~~$ $~~~~~~~~+~~\left(1-f\left(r\right)\right)U_{morse}\left(r\right)$
zblb	ZBL mixed with Buckingham	$Z_{1}$, $Z_{2}$, $r_{m}$, $\xi$, $A$, $\rho$, $C$	$U(r) = f\left(r\right)U_{ZBL}\left(r\right)~~+~~$ $~~~~~+~~\left(1-f\left(r\right)\right)U_{buckingham}\left(r\right)$
mlj	Lennard-Jones tapered with MDF	$\epsilon$, $\sigma$, $r_i$	$U(r) = 4 \epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r}\right)^{6} \right]f(r)$
mbuc	Buckingham tapered with MDF	$A$, $\rho$, $C$, $r_i$	$U(r) = \left[A~\exp \left( -\frac{r}{\rho} \right) - \frac{C}{r^{6}}\right]f(r)$
m126	12-6 Lennard-Jones tapered with MDF	$A$, $B$, $r_i$	$U(r) = \left( \frac{A}{r^{12}} - \frac{B}{r^{6}} \right)f(r)$

$^{\dagger}$ Note: in this formula the terms $\alpha$, $\beta$ and $\gamma$ are compound expressions involving the variables $E_{o},n,m,r_{0}$ and $r_{c}$. See Section Short Ranged (van der Waals) Potentials for further details.

$^{\ddagger}$ Note: All local potential cutoffs, $r_{c}$, default to the general van der Waals cutoff, rvdw, or the general domain decomposition cutoff, rcut, if all are unspecified or set to zero in the FIELD file! Any non-zero values of $r_{c}$ in the FIELD file will take precedence over any value specified in the CONTROL file when determining the general van der Waals cutoff rvdw, with two notable consequences: (1) a value for rvdw given in the CONTROL file cannot override any non-zero values given in the FIELD file, and (2) the largest non-zero cutoff value given in the FIELD file (including the WCA equivalent $2^{1 \over 6}~\sigma + \Delta$) will always be used as the general van der Waals cutoff. The general domain decomposition cutoff rcut will only be set to the general van der Waals cutoff if it is not otherwise specified in the CONTROL file or it is smaller. Therefore, we strongly recommend users check the OUTPUT file for the van der Waals cutoff DL_POLY_5 has determined for a given simulation, as this value will be used for all potentials involving force-shifting or long-range corrections!

$^{*}$ Note: $\Delta$ defaults to zero if $|\Delta| > 0.5~\sigma$ or it is not specified in the FIELD file.

metal n where n is the number of metal potentials to be entered. It is followed by n records, each specifying a particular metal potential in the following manner:

atmnam 1      a8      first atom type
atmnam 2      a8      second atom type
key           a4      potential key, see Table (14)
variable 1    real    potential parameter, see Table (14)
variable 2    real    potential parameter, see Table (14)
variable 3    real    potential parameter, see Table (14)
variable 4    real    potential parameter, see Table (14)
variable 5    real    potential parameter, see Table (14)
variable 6    real    potential parameter, see Table (14)
variable 7    real    potential parameter, see Table (14)
variable 8    real    potential parameter, see Table (14)
variable 9    real    potential parameter, see Table (14)

The variables pertaining to each potential are described in Table (14).

Table 14 Metal Potential¶
Key	Potential Type	Variables	Functional Form
eam	EAM		tabulated potential in TABEAM
eeam	EEAM		tabulated potential in TABEAM
2bea	2BEAM		tabulated potential in TABEAM
2bee	2BEEAM		tabulated potential in TABEAM
fnsc	Finnis-Sinclair	$c_{0}$, $c_{1}$, $c_{2}$, $c$, $A$, $d$, $\beta$	$U_{i}(r) = {1 \over 2} \sum\limits_{j \ne i} (r_{ij}-c)^{2} (c_{0}+c_{1}r_{ij}+c_{2}r_{ij}^{2}) - A \sqrt{\rho_{i}}~;$ $\rho_{i} = \sum\limits_{j\ne i} \left[(r_{ij}-d)^{2} + \beta \frac{(r_{ij}-d)^{3}}{d}\right]$
exfs	Extended Finnis-Sinclair	$c_{0}$, $c_{1}$, $c_{2}$, $c_{3}$, $c_{4}$, $c$, $A$, $d$, $B$	$U_{i}(r) = {1 \over 2} \sum\limits_{j \ne i} (r_{ij}-c)^{2} (c_{0}+c_{1}r_{ij}+c_{2}r_{ij}^{2}+c_{3}r_{ij}^{3}+c_{4}r_{ij}^{4})$ $- A \sqrt{\rho_{i}}~~~;~~~\rho_{i} = \sum\limits_{j\ne i} \left[(r_{ij}-d)^{2} + B^{2} (r_{ij}-d)^{4}\right]$
stch	Sutton-Chen	$\epsilon$, $a$, $n$, $m$, $c$	$U_{i}(r) = \epsilon \left[ {1 \over 2} \sum\limits_{j \ne i} \left(\frac{a}{r_{ij}} \right)^{n} - c \sqrt{\rho_{i}} \right]~;~\rho_{i} = \sum\limits_{j\ne i} \left(\frac{a}{r_{ij}}\right)^{m}$
gupt	Gupta	$A$, $r_{0}$, $p$, $B$, $q_{ij}$	$U_{i}(r) = \sum\limits_{j \ne i} A \exp \left(-p \frac{r_{ij}-r_{0}}{r_{0}}\right) - B \sqrt{\rho_{i}}~;$ $\rho_{i} = \sum\limits_{j\ne i} \exp \left(-2 q_{ij} \frac{r_{ij}-r_{0}}{r_{0}}\right)$
mbpc	MBPC$^{\dagger}$	$\epsilon$, $a$, $m$, $\alpha$, $r_{\rm o}$	$U_{i}(r) = -\epsilon \sqrt{\rho_{i}}~;~\rho_{i} = \sum\limits_{j\ne i} \left(\frac{a}{r_{ij}^{m}}\right) \frac{1}{2}\left[1+{\rm erf}\left(\alpha(r_{ij}-r_{\rm o})\right)\right]$

$^{\dagger}$Note that the parameters $\alpha$ and $r_{\rm o}$ must be the same for all defined potentials of this type. DL_POLY_5 will set $\alpha={\rm Max}(0,\alpha_{pq})$ and $r_{\rm o}={\rm Max}(0,r_{{\rm o}\_pq})$ for all defined interactions of this type between species $p$ and $q$. If after this any is left undefined, i.e. zero, the undefined entities will be set to their defaults: $\alpha=20$ and $r_{\rm o}={\rm Min}(1.5,0.2~r_{\rm cut})$.

rdf n where n is the number of RDF pairs to be entered. It is followed by n records, each specifying a particular RDF pair in the following manner:
```
atmnam 1    a8    first atom type
atmnam 2    a8    second atom type
```
By default in and DL_POLY_5 every vdw and met potential specifies an RDF pair. If the control option rdf f is specified in the CONTROL file then all pairs defined in vdw and/or met potentials sections will also have their RDF calculated. The user has two choices to enable the calculation of RDFs in systems with force fields that do not have vdw and/or met potentials: (i) to define fictitious potentials with zero contributions or (ii) to use rdf n option - which not only provides a neater way for specification of RDF pairs but also better memory efficiency since DL_POLY_5 will not allocate (additional) potential arrays for fictitious interactions that will not be used. (This option is not available in .)

Note that rdf and vdw/met are not complementary - i.e. if the former is used in FIELD none of the pairs defined by the latter will be considered for RDF calculations.

The selected RDFs are calculated in the rdf_collect, rdf_excl_collect, rdf_frzn_collect and rdf_compute by collecting distance information from all two-body pairs as encountered in the Verlet neighbour list created in the link_cell_pairs routine within the two_body_forces routine. In the construction of the Verlet neighbour list, pairs of particles (part of the exclusion list) are excluded. The exclusion list contains particles that are part of:
- core-shell units
- bond constraints
- chemical bonds, that are NOT distance restraints
- valence angles, that are NOT distance restraints
- dihedrals
- inversions
- frozen particles
RDF pairs containing type(s) of particles that fall in this list will be polluted. However, there are many ways to overcome such effects.

tersoff n where n is the number of specified Tersoff potentials. There are two types of Tersoff potential forms that cannot be mixed (used simultaneously). They are shorthanded as ters and kihs atomkeys.

ters atomkey expects $2n$ records specifying $n$ particular Tersoff single atom type parameter sets and $n(n+1)/2$ records specifying cross atom type parameter sets in the following manner:

potential 1 : record 1
atmnam        a8      atom type
key           a4      potential key, see Table (15)
variable 1    real    potential parameter, see Table (15)
variable 2    real    potential parameter, see Table (15)
variable 3    real    potential parameter, see Table (15)
variable 4    real    potential parameter, see Table (15)
variable 5    real    cutoff range for this potential (Å or DPD length units)
potential 1 : record 2
variable 6    real    potential parameter, see Table (15)
variable 7    real    potential parameter, see Table (15)
variable 8    real    potential parameter, see Table (15)
variable 9    real    potential parameter, see Table (15)
variable 10   real    potential parameter, see Table (15)
variable 11   real    potential parameter, see Table (15)
    ...       ...       ...
    ...       ...       ...
potential n : record 2n-1
    ...       ...       ...
potential n : record 2n
    ...       ...       ...
cross term 1 : record 2n+1
atmnam 1      a8      first atom type
atmnam 2      a8      second atom type
variable a    real    potential parameter, see Table (15)
variable b    real    potential parameter, see Table (15)
variable c    real    potential parameter, see Table (15)
    ...       ...       ...
    ...       ...       ...
cross term n(n+1)/2 : record 2n+n(n+1)/2`
    ...       ...       ...

kihs atomkey expects $3n$ records specifying $n$ particular Tersoff single atom type parameter sets in the following manner:

potential 1 : record 1
atmnam        a8        atom type
key           a4        potential key, see Table (15)
variable 1    real      potential parameter, see Table (15)
variable 2    real      potential parameter, see Table (15)
variable 3    real      potential parameter, see Table (15)
variable 4    real      potential parameter, see Table (15)
variable 5    real      cutoff range for this potential (Å or DPD length units)
potential 1 : record 2
variable 6    real      potential parameter, see Table (15)
variable 7    real      potential parameter, see Table (15)
variable 8    real      potential parameter, see Table (15)
variable 9    real      potential parameter, see Table (15)
variable 10   real      potential parameter, see Table (15)
variable 11   real      potential parameter, see Table (15)
potential 1 : record 3
variable 12   real      potential parameter, see Table (15)
variable 13   real      potential parameter, see Table (15)
variable 14   real      potential parameter, see Table (15)
variable 15   real      potential parameter, see Table (15)
variable 16   real      potential parameter, see Table (15)
    ...       ...         ...
    ...       ...         ...
potential n : record 2n-1
    ...       ...         ...
potential n : record 2n
    ...       ...         ...

The variables pertaining to each potential are described in Table (15).

Note that the fifth variable is the range at which the particular tersoff potential is truncated. The distance is in Å (or DPD length units).

Table 15 Tersoff Potential¶
Key	Potential Type	Variables	Functional Form
ters	Tersoff	$A$, $a$, $B$, $b$, $R$	as shown in Section Tersoff Potentials
	(single)	$S$, $\beta$, $\eta$, $c$, $d$, $h$
	(cross)	$\chi$, $\omega$, $\delta$
kihs	KIHS	$A$, $a$, $B$, $b$, $R$, $S$, $\eta$, $\delta$, $c_{1}$, $c_{2}$ $c_{3}$, $c_{4}$, $c_{5}$, $h$, $\alpha$, $\beta$	Section Tersoff Potentials

tbp n where n is the number of three-body potentials to be entered. It is followed by n records, each specifying a particular three-body potential in the following manner:

atmnam 1 (i)      a8      first atom type
atmnam 2 (j)      a8      second (central) atom type
atmnam 3 (k)      a8      third atom type
key               a4      potential key, see Table (16)
variable 1        real    potential parameter, see Table (16)
variable 2        real    potential parameter, see Table (16)
variable 3        real    potential parameter, see Table (16)
variable 4        real    potential parameter, see Table (16)
variable 5        real    cutoff range for this potential (Å or DPD length units)

The variables pertaining to each potential are described in Table (16).

Note that the fifth variable is the range at which the three body potential is truncated. The distance is in Å, measured from the central atom. Note that interactions defined with less than four potential parameters must provide zeros for the ones up to the cutoff one.

Table 16 Three-body Potentials¶
Key	Potential Type	Variables	Functional Form
harm	Harmonic	$k$, $\theta_{0}$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}$
thrm	Truncated harmonic	$k$, $\theta_{0}$, $\rho$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}~\exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]$
shrm	Screened harmonic	$k$, $\theta_{0}$, $\rho_{1}$, $\rho_{2}$	$U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]$
bvs1	Screened Vessal ¹⁴¹	$k$ $\theta_{0}$, $\rho_{1}$, $\rho_{2}$	$U(\theta) = {k \over 8(\theta_{0}-\pi)^{2}} \left\{ \left[ (\theta_{0} -\pi)^{2} -(\theta-\pi)^{2}\right]^{2}\right\} \times$ $\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]$
bvs2	Truncated Vessal ¹²³	$k$, $\theta_{0}$, $a$, $\rho$	$U(\theta) = k~(\theta-\theta_{0})^{2}~\big[ \theta^a (\theta-\theta_{0})^{2} (\theta+\theta_{0}-2\pi){^{2}}$ $+ {a \over 2} \pi^{a-1}(\theta_{0}-\pi)^{3}\big]~\exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]$
hbnd	H-bond ⁷⁵	$D_{hb}$, $R_{hb}$	$U(\theta) = D_{hb}~\cos^{4}(\theta) \times$ $[5(R_{hb}/r_{jk})^{12} - 6(R_{hb}/r_{jk})^{10}]$

$\theta$ is the $i$-$j$-$k$ angle.

fbp n where n is the number of four-body potentials to be entered. It is followed by n records, each specifying a particular four-body potential in the following manner:

atmnam 1 (i)    a8      first (central) atom type
atmnam 2 (j)    a8      second atom type
atmnam 3 (k)    a8      third atom type
atmnam 4 (l)    a8      fourth atom type
key             a4      potential key, see Table (17)
variable 1      real    potential parameter, see Table (17)
variable 2      real    potential parameter, see Table (17)
variable 3      real    cutoff range for this potential (Å or DPD length units)

The variables pertaining to each potential are described in Table (17).

Note that the third variable is the range at which the four-body potential is truncated. The distance is in Å, measured from the central atom. Note that interactions idefined with less than two potential parameters must provide zeros for the ones up to the cutoff one.

Table 17 Four-body Potentials¶
Key	Potential Type	Variables	Functional Form $^\ddagger$
harm	Harmonic	$k$, $\phi_{0}$	$U(\phi) = {k \over 2}~(\phi - \phi_{0})^{2}$
hcos	Harmonic cosine	$k$, $\phi_{0}$	$U(\phi) = {k \over 2}~(\cos(\phi) - \cos(\phi_{0}))^{2}$
plan	Planar	$A$	$U(\phi) = A~[1 - \cos(\phi)]$

$^\ddagger$ $\phi$ is the $i$-$j$-$k$-$l$ four-body angle.

External Field¶

The presence of an external field is flagged by the directive: extern The following line in the FIELD file must contain another directiveindicating what type of field is to be applied, followed by the field parameters in the following manner:

field key     a4      external field key, see Table (18)
variable 1    real    potential parameter, see Table (18)
variable 2    real    potential parameter, see Table (18)
variable 3    real    potential parameter, see Table (18)
variable 4    real    potential parameter, see Table (18)
variable 5    real    potential parameter, see Table (18)
variable 6    real    potential parameter, see Table (18)

The variables pertaining to each field potential are described in Table (18).

Note: only one type of field can be applied at a time.

Table 18 External Fields¶
Key	Potential Type	Variables	Functional Form
elec	Electric Field	$E_{x}$, $E_{y}$, $E_{z}$	$\underline{F} = q \; \underline{E}$
oshr	Oscillating Shear	$A$, $n$	$\underline{F}_x = A \; \cos(2 n\pi \cdot z/L_{z})$
shrx	Continuous Shear	$A$, $z_{0}$	$\underline{v}_x = \frac{A}{2} \frac{\|z\|}{z}~:~\|z\|~>~z_{0}$
grav	Gravitational Field	$G_{x}$, $G_{y}$, $G_{z}$	$\underline{F} = m \; \underline{G}$
magn	Magnetic Field	$H_{x}$, $H_{y}$, $H_{z}$	$\underline{F} = q \; (\underline{v} \times \underline{H})$
sphr	Containing Sphere	$A$, $R_{0}$, $n$, $R_{\rm cut}$	$\underline{F} = A \; (R_{0} - r)^{-n}~:~r > R_{\rm cut}$
zbnd	Repulsive Wall	$A$, $z_{0}$, $f=\pm 1$	$\underline{F} = A \; (z_{0} - z) ~:~f \cdot z > f \cdot z_{0}$
xpis	X-Piston	$i_{gid}$, $j_{gid}$, $P_{k {\rm atm}}$	$\underline{F}_{x} = \frac{P \cdot \underline{\textrm{Area}}(\perp\textrm{X-}dir)}{\left[\sum_{k=i}^{j}m_{k}\right] \textrm{/} m_{k}} : \forall~k=i,..,j$
zres	Molecule in HR Zone	$i_{gid}$, $j_{gid}$, $A$, $z_{mn}$, $z_{mx}$	$\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z_{cm}-z_{mx}) & z_{cm} > z_{mx} \\A(z_{mn}-z_{cm}) & z_{cm} < z_{mn}\end{array} \right.$
zrs-	HR Zone (push out)	$i_{gid}$, $j_{gid}$, $A$, $z_{mn}$, $z_{mx}$	$\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z-z_{mx}) & z \ge \frac{z_{mx}+z_{mn}}{2} \\A(z_{mn}-z) & z < \frac{z_{mx}+z_{mn}}{2}\end{array} \right.$
zrs+	HR Zone (pull in)	$i_{gid}$, $j_{gid}$, $A$, $z_{mn}$, $z_{mx}$	$\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z-z_{mx}) & z > z_{mx} \\A(z_{mn}-z) & z < z_{mn}\end{array} \right.$
osel	Osc. Electric Field	$E_{x}$, $E_{y}$, $E_{z}$, $\omega_{{\rm ps}^{-1}}$	$\underline{F} = q \; \underline{E} \; \sin(2 \pi \omega t)$
ushr	umbrella sampling ¹³⁶	$i^{A}_{gid}$, $j^{A}_{gid}$, $i^{B}_{gid}$, $j^{B}_{gid}$, $k$, $R_{0}$	$U_{AB} = \frac{k}{2} (R_{AB}-R_{0})^{2}$
	harm. constraint ⁵⁸

Note

External force parameters are read in terms of the specified energy units and the general DL_POLY units so that the two sides of the equation defining the field are balanced. For example, the magnetic field units, $\underline{H} = (H_{1},H_{2},H_{3})$, in the DL_POLY FIELD scope will follow from the interaction definition as seen in Table (18):

(214)¶\[\begin{split} \begin{aligned} \underline{F} =& q \; (\underline{v} \times \underline{H})~~~~~~\texttt{therefore}, \nonumber \\ \left[H\right] =& \frac{\left[F\right]}{\left[q\right]~\left[v\right]} = \frac{\left[m\right]~\left[a\right]}{\left[q\right]~\left[v\right]} \nonumber \\ \left[H\right] =& \frac{\texttt{Dalton~Å/ps$^{2}$}}{\texttt{proton~Å/ps}} = \frac{\texttt{Dalton}}{\texttt{proton~ps}} \\ \left[H\right] =& 1.037837512 \times 10^{4}~\texttt{Tesla} \\ H(DL\_POLY) =& H(MKS)~1.037837512 \times 10^{4}~~. \nonumber \end{aligned}\end{split}\]

Thus to apply a magnetic field of 1 Tesla along the $y$ axis, one could specify in FIELD the following:

UNITS internal
...
external
magnetic 0  1/1.037837512e04   0
...

when working in DL_POLY internal units. If we worked in unit units then

\[\begin{split}\begin{aligned} H &\propto& Energy \nonumber \\ Energy(DL\_POLY) &=& Energy(unit) ~ k_{unit \to DL\_POLY} \nonumber \\ H(DL\_POLY) &=& H(MKS)~1.037837512 \times 10^{4} \\ H(unit) &=& H(MKS)~\frac{1.037837512 \times 10^{4}}{k_{unit \to DL\_POLY}}~~, \nonumber\end{aligned}\end{split}\]

with the following conversion factors values:

\[\begin{split}\begin{aligned} k_{eV \to DL\_POLY} &=& 9648.530821 \nonumber \\ k_{kcal/mol \to DL\_POLY} &=& 418.4 \nonumber \\ k_{kJ/mol \to DL\_POLY} &=& 100.0 \\ k_{K/Boltz \to DL\_POLY} &=& 0.831451115 \nonumber \\ k_{DPD \to DL\_POLY} &=& 1.0 \nonumber \\ k_{DL\_POLY \to DL\_POLY} &=& 1.0~~. \nonumber\end{aligned}\end{split}\]

Obviously, for $eV$ units

\[\begin{split}\begin{aligned} H(unit) &=& H(MKS)~\frac{1.037837512 \times 10^{4}}{k_{unit \to DL\_POLY}} \nonumber \\ k_{eV \to DL\_POLY} &=& 9648.530821 \\ H(eV) &=& H(MKS)~1.07564305~~, \nonumber\end{aligned}\end{split}\]

the FIELD file should be amended to read:

UNITS eV
...
external
magnetic 0  1/1.07564305   0
...

Note that if $DPD$ units are in use, the ‘real-life’ conversion factors will differ depending on the choices for the DPD mass, length and energy units. (The DPD charge unit is equal to the charge magnitude of an electron.) As such, any external field parameters will always be given in DPD units, i.e. electric fields in $\frac{[E]}{[q] [L]}$ rather than in volts per angstrom, magnetic fields in $\frac{[E] [t]}{[q] [L]}$ instead of teslas, and piston pressures in $\frac{[E]}{[L]^3}$ in place of kilo-atmospheres.

crd¶

crd i j k

This section defines the atomic pairs to use for both coordination and angular distribution calculations. The start of this section is given by the keyword crd followed by the i j k records:

i = The number of unique atomic pairs

i number of the following pattern: atom1 atom2 bond_length

j = The number of unique atoms you are building the pairs from

j number of the following pattern: atom minimun_displacement

k = The number of lists to build coordination and angular distribution between

k number of the following pattern: atom atom … ‘-‘ atom atom …

Example in the FIELD File

crd 2 3 1
B O 1.9
Si O 2.0
B 0.9
Si 1.0
O 1.0
Si B - O

Closing the FIELD File¶

The FIELD file must be closed with the directive: close which signals the end of the force field data. Without this directive DL_POLY_5 will abort.

The MPOLES File¶

The MPOLES file serves to define the multipolar momenta of the sites defined in FIELD (see Section The FIELD File). It is only read by DL_POLY_5 when the directive is specified at the top of FIELD. The file is read by the subroutine read_mpoles. The MPOLES file and has the same molecular topology syntax and rules as FIELD except that (i) it only understands the stoichiometry information; and (ii) will abort should more FIELD like detail is supplied. Thus the molecular topology/stoichiometry information must explicitly match that in FIELD.

Example from a water system with multipolar momenta beyond the single point charges is shown below:

Puddle Test Case (32000 AMOEBA Waters)

MOLECULES 1

SPC WATER
NUMMOLS   32000

ATOMS     3

O  2  1  0.837  0.39
-0.51966
 0.00000   0.00000   0.14279
 0.37928   0.00000   0.00000  -0.41809   0.00000   0.03881

H  2  2  0.496  0.39
0.25983
-0.03859   0.00000  -0.05818
-0.03673   0.00000  -0.00203  -0.10739   0.00000   0.14412

FINISH

CLOSE

The MPOLES File Format¶

Definitions of Variables in the MPOLES File¶

The file divides into two sections: general information, molecular descriptions.

General information¶

The first viable record in the MPOLES file is the title, which is mandatory.

**record 1**
header    a200      field file header

Molecular details¶

The entry of the molecular details begins with the mandatory directive:

molecules n

Immediately following the molecules directive, are the records defining individual molecules:

name-of-molecule which can be any character string up to 200 characters in length. (Note: this is not a directive, just a simple character string.)
nummols n where n is the number of times a molecule of this type appears in the simulated system. The molecular data then follow in subsequent records:
atoms n where n indicates the number of atoms in this type of molecule. A number of records follow, each giving details of the atoms in the molecule i.e. site names, masses and charges. Each record carries the entries:
```
sitnam    a8        atomic site name
order     integer   multipolar order supplied
nrept     integer   repeat counter
\alpha    real(1)   optional atomic polarisability (in Å³ or cubed DPD length units)
a         real(1)   optional Thole dumping factor
```
The integer nrept may be omitted (in which case a value or 1 is assumed) only if no further optional directives are provided! A number greater than 1 specified here indicates that the next (nrept-1) entries in the CONFIG file are ascribed the atomic characteristics given in the current record. The sum of the repeat numbers for all atoms in a molecule should equal the number specified by the atoms directive.

The atomic polarisability, $\alpha$, and the Thole ¹²⁹ dumping factor, $a$, are optional and are only parsed for the core (non-Druder) particles of core-shell units. If all polarisabilities as well as core and shell charges and associated force constants are well defined then $a$ may default for (Druder nuclei) core particles if no Thole dumping is specified in both MPOLES and CONTROL.

For each pole order specified DL_POLY_5 will read a new line that will specify the pole order momenta. If a site has a specified pole order, $m$, smaller than the one specified in FIELD, $n$, then the $m+1$ to $n$ order poles’ momenta will be zero. Similarly, if a site has momenta of poles of higher order than the one specified in FIELD, $n$, these will not be processed by .

DL_POLY_5 will read same order momenta in their natural vector format:
```
charge        real(1)     scalar q (in protons)
dipole        real(3)     vector (x, y, z) (in protons/Å or protons/DPD length unit)
quadrupole    real(6)     vector (xx, xy, xz, yy, yz, zz) (in protons/Å² or protons/squared DPD length unit)
octupole      real(10)    vector (xxx, xxy, xxz, xyy, xyz, xzz, yyy, yyz, yzz, zzz) (in protons/Å³ or protons/cubed DPD length unit)
hexadecapole  real(15)    vector (xxxx, xxxy, xxxz,...., yyyz, yyzz, yzzz, zzzz)` (in protons/Å⁴ or protons/quartic DPD length unit)
```
Note that the charge values supplied in FIELD will be overwritten with those supplied here!

Note, although algorithms in DL_POLY_5 could in principle handle any high pole order summation, in practice, however, DL_POLY_5 will abort if the order is higher than hexadecapole (order 4)!
finish This directive is entered to signal DL_POLY_5 that the entry of the details of a molecule has been completed.

The entries for a second molecule may now be entered, beginning with the name-of-molecule record and ending with the finish directive.

The cycle is repeated until all the types of molecules indicated by the molecules directive have been entered.

Closing the MPOLES File¶

The MPOLES file must be closed with the directive: close after which DL_POLY_5 will resume will processing the intermolecular part of FIELD.

The REFERENCE File¶

The REFERENCE has the same format and structure as CONFIG (see Section The CONFIG File) file with the exception that $\texttt{imcon}~\textbf{MUST~BE} \ne 0$. REFERENCE may contain more or less particles than CONFIG does and may have particles with identities that are not defined in FIELD (see Section The FIELD File). The positions of these particles are used to define the crystalline lattice sites to which the particles in CONFIG compare during simulation when the defect detection option, defects, is used. REFERENCE is read by the subroutine defects_reference_read.

The REVOLD File¶

This file contains statistics arrays from a previous job. It is not required if the current job is not a continuation of a previous run (i.e. if the restart directive is not present in the CONTROL file - see above). The file is unformatted and therefore not human readable. DL_POLY_5 normally produces the file REVIVE (see Section The REVIVE File) at the end of a job which contains the statistics data. REVIVE should be copied to REVOLD before a continuation run commences. This may be done by the copy macro supplied in the execute sub-directory of .

Format¶

The REVOLD file is unformatted. All variables appearing are written in native working precision (see Section Choosing Ewald Sum Variables) real representation. Nominally, integer quantities (e.g. the timestep number nstep) are represented by the nearest real number. The contents are as follows (the dimensions of array variables are given in brackets, in terms of parameters from the setup_module file - see Section Comments on setup_module).

record 1:
nstep       timestep of final configuration
numacc      number of configurations used in averages
numrdf      number of configurations used in RDF averages
numzdn      number of configurations used in Z-density averages
time        elapsed simulation time
tmst        elapsed simulation before averages were switched on
chit        thermostat related quantity (first)
chip        barostat related quantity
cint        thermostat related quantity (second)
record 2:
eta         scaling factors for simulation cell matrix elements (9)
record 3:
stpval      instantaneous values of thermodynamic variables (mxnstk)
record 4:
sumval      average values of thermodynamic variables (mxnstk)
record 5:
ssqval      fluctuation (squared) of thermodynamic variables (mxnstk)
record 6:
zumval      running totals of thermodynamic variables (mxnstk)
record 7:
ravval      rolling averages of thermodynamic variables (mxnstk)
record 8:
stkval      stacked values of thermodynamic variables (mxstak\(\times\)mxnstk)
record 9:
strcon      constraint bond stress (9)
record 10:
strpmf      PMF constraint stress (9)
record 11:
stress      atomic stress (9)
record 12: (Optional)
rdf         RDF arrays (mxgrdf\(\times\)mxrdf)
record 13: (Optional)
usr         umbrella sampling RDF array (mxgusr)
record 14: (Optional)
zdens       Z-density array (mxgrdf\(\times\)mxatyp)
record 15: (Optional)
vaf         VAF arrays (sizes dependent on sampling frequency and VAF bin size)

Further Comments¶

Note that different versions of DL_POLY_5 may have a different order of the above parameters or include more or less such. Therefore different versions of DL_POLY_5 may render any existing REVOLD file unreadable by the code.

The TABVDW File¶

The TABVDW file provides an alternative way of reading in the short range potentials - in tabular form. This is particularly useful if an analytical form of the potential does not exist or is too complicated to specify in the vdw_generate subroutine. The table file is read by the subroutine vdw_table_read (see Chapter Source Code).

The option of using tabulated potentials is specified in the FIELD file (see above). The specific potentials that are to be tabulated are indicated by the use of the tab keyword on the record defining the short range potential (see Table (13)).

The TABVDW File Format¶

The file is free-formatted but blank and commented lines are not allowed.

Definitions of Variables¶

record 1
header    a200      file header
record 2
delpot    real      mesh resolution in Å or DPD length units (delpot = cutpot / (ngrid-4))
cutpot    real      cutoff used to define tables in Å or DPD length units
ngrid     integer   number of grid points in tables

The subsequent records define each tabulated potential in turn, in the order indicated by the specification in the FIELD file. Each potential is defined by a header record and a set of data records with the potential and force tables.

header record:
atom 1      a8      first atom type
atom 2      a8      second atom type
potential data records: (number of data records = Int((ngrid+3)/4))
data 1      real    data item 1
data 2      real    data item 2
data 3      real    data item 3
data 4      real    data item 4
force data records: (number of data records = Int((ngrid+3)/4))
data 1      real    data item 1
data 2      real    data item 2
data 3      real    data item 3
data 4      real    data item 4

Further Comments¶

It should be noted that the number of grid points in the TABVDW file should not be less than the number of grid points DL_POLY_5 is expecting. (This number is given by the parameter mxgvdw calculated in the setup_module file - see Section Note on the Interpolation Scheme and Comments on setup_module.) DL_POLY_5 will re-interpolate the tables if $\texttt{delpot} = \frac{\texttt{cutpot}}{{\texttt{ngrid}-4}}~<~\texttt{dlrvdw} = \frac{\texttt{rvdw}}{{\texttt{mxgvdw}-4}}$ (usually when $\texttt{ngrid}>\texttt{mxgvdw}$), but will abort if $\texttt{delpot}>\texttt{dlrvdw}$.

The potential and force tables are used to fill the internal arrays vvdw and gvdw respectively (see Section Short Ranged (van der Waals) Potentials). The contents of force arrays are derived from the potential via the formula:

\[G(r)=-r\frac{\partial}{\partial r}U(r)~~.\]

Note

this is not the same as the true force.

During simulation, interactions beyond distance $\texttt{cutpot}$ are discarded.

The TABEAM File¶

The TABEAM file contains the tabulated potential functions (no explicit analytic form) describing the EAM or EEAM metal interactions in the MD system. This file is read by the subroutine metal_table_read (see Chapter Source Code).

The TABEAM File Format¶

The file is free-formatted but blank and commented lines are not allowed.

Definitions of Variables¶

record 1
header    a200      file header
record 2
numpot    integer   number of potential functions in file

For an $n$ component alloy, numpot is

$n(n+5)/2$ for the EAM potential or
$3n(n+1)/2$ for the EEAM potential or
$n(n+4)$ for the 2BEAM potential or
$5n(n+1)/2$ for the 2BEEAM potential.

The subsequent records for an $n$ component alloy define $n(n+1)/2$ cross pair potential functions - pairs keyword and

EAM: $n$ embedding functions (one for each atom type) - embedding keyword, $n$ electron density functions (one for each atom type) - density keyword;
EEAM: $n$ embedding functions (one for each atom type) and $n^{2}$ for the EEAM potential (one for each non-commuting pair of atoms types);
2BEAM: $n$ $s$-band embedding functions (one for each atom type) - sembedding keyword, $n(n+1)/2$ $s$-band density functions (one for each cross-pair of atoms types) - sdensity keyword, $n$ $d$-band embedding functions (one for each atom type) - dembedding or embeedding keyword, and $n$ $d$-band density functions (one for each atom types) - ddensity or density keyword;
2BEEAM: $n$ $s$-band embedding functions (one for each atom type) - sembedding keyword, $n^{2}$ $s$-band density functions (one for each non-commuting pair of atoms types) - sdensity keyword, $n$ $d$-band embedding functions (one for each atom type) - dembedding or embeedding keyword, and $n^{2}$ $d$-band density functions (one for each non-commuting pair of atoms types) - ddensity or density keyword.

The functions may appear in any random order in TABEAM as their identification is based on their unique keyword, defined first in the function’s header record. The header record is followed by predefined number of data records as a maximum of four data per record are read in - allowing for incompletion of the very last record.

header record:
keyword     a4        type of EAM function: pair, embed or density, with 2B extension alternatives for the s-band - [sembed and sdensity] and d-band - dembed = embed and ddensity = density
atom 1      a8        first atom type
atom 2      a8        second atom type - only specified for pair potential functions and for the (i) density functions in the EEAM potential case or (ii) sdensity functions in the 2BEAM potential case or (iii) sden and dden functions in the 2BEEAM potential case
ngrid       integer   number of function data points to read in
limit 1     real      lower interpolation limit in Å or DPD length units for dens/sden/dden and pair or in density units for embe/semb/demb
limit 2     real      upper interpolation limit in Å or DPD length units for dens/sden/dden and pair or in density units for embe/semb/demb
function data records: (number of data records = Int((ngrid+3)/4))
data 1      real      data item 1
data 2      real      data item 2
data 3      real      data item 3
data 4      real      data item 4

Further Comments¶

The tabled data are used to fill the internal arrays vmet, dmet and fmet, and optionally dmes and fmes for the 2B extensions of EAM and EEAM (see Section Metal Potentials). The force arrays are generated from these (by the metal_table_derivatives routine) using a five point interpolation procedure. During simulation, interactions beyond distance $Min(r_{\rm cut},\texttt{limit~2})$ are discarded, whereas interactions at distances shorter than limit 1 will cause the simulation to abort. For the purpose of extrapolating the embedding functions $F(\rho)$ beyond its limit 2 specified in the tabulated array, it is assumed that

\[F(\rho > \texttt{limit~2}) = F(\rho = \texttt{limit~2})~~.\]

The simulation will however abort if any local density is less than the limit 1 for its corresponding embedding function.

It is worth noting that in the 2BEAM and 2BEEAM the $s$-band contribution is usually only for the alloy component, so that local concentrations of a single element revert to the standard EAM or EEAM! In such case, the densities functions must be zeroed in the DL_POLY_5 TABEAM file. A convenient way to do this, for example, will be data record of the type:

SDEN   Atom1   Atom1   1   0   1   0

The TABBND, TABANG, TABDIH & TABINV Files¶

DL_POLY_5 allows the specification of tabulated data for intramolecular interactions:

TABBND - for chemical bonds potentials - distance dependent
TABANG - for bond angles potentials - angle dependent
TABDIH - for dihedrals (torsional) potentials - angle dependent
TABINV - for inversions potentials - angle dependent.

The files have the same formatting rules with examples shown in Section Setting up Tabulated Intramolecular Force-Field Files. Refer to Section User-Defined Coarse-Grain Models with Tabulated Force-Fields for their derivation and usage in coarse grained model systems.

Definitions of Variables¶

record 1
header    a200      file header
record 2
#         a1        a hash (#) symbol
cutpot    real      cutoff in Å or DPD length units - only expected in TABBND as the cutoff ranges are known for TABANG, TABDIH & TABINV
ngrid     integer   number of grid points in table for all potentials
record 3
#         a1        a hash (#) symbol

empty record:
id record:
#         a1      a hash (#) symbol
atom 1    a8      first atom type
atom 2    a8      second atom type
atom 3    a8      third atom type - only required for TABANG
atom 4    a8      forth atom type - only required for TABDIH & TABINV
interaction data records 0/1-ngrid:
abscissa  real    consecutive value over the full cutoff/range in Å or DPD length units for TABBND and degrees for TABANG, TABDIH & TABINV
potential real    potential at the abscissa grid point in units as specified in FIELD
force     real    complementary force (virial for TABBND) value

Further Comments¶

It should be noted that the number of grid points in the table files should not be less than the number of grid points DL_POLY_5 is expecting. For more information the reader is advised to examine setup_module and inspect the mxgint variables, where int refers to bnd for bonds, ang for angles, dih for dihedrals and inv for inversions.

The potential and force tables are used to fill the internal arrays vint and gint for the respective intrmolecular potential (see Chapter Force Fields).

The contents of force arrays for TABBND are derived from the potential via the formula:

\[G(r)=-r\frac{\partial}{\partial r}U(r)~~.\]

Note

this is not the same as the true force. During simulation, interactions beyond distance $\texttt{cutpot}$ will bring the run to a controlled termination.

The DUMP_E File¶

The DUMP_E file contains the values of coarse-grained electronic temperature (CET) cells from a previous job with the two-temperature model (TTM). It is not required if the current job is not a continuation of a previous run (i.e. if the restart directive is not present in the CONTROL file - see above). The two-line header consists of the following records:

record 1:
eltsys(1)     integer     number of CET cells in x-direction
eltsys(2)     integer     number of CET cells in y-direction
eltsys(3)     integer     number of CET cells in z-direction
record 2:
nstep         integer     timestep of current electronic temperature profile
time          real        elapsed simulation time
depostart     real        time when energy deposition started
depoend       real        time when energy deposition ended or is due to end

and each subsequent record is formatted as follows:

x       integer     x-coordinate of CET cell
y       integer     y-coordinate of CET cell
z       integer     z-coordinate of CET cell
eltemp  real        electronic temperature at current CET cell (k)

with the origin of CET cell coordinates at the centre of the grid.

The file is read by the subroutine ttm_system_init. It will not be accepted by DL_POLY_5 if the number of CET cells in each direction does not match the values given in the CONTROL file or insufficient electronic temperatures are supplied. No matching up of restart timestep or elapsed simulation time between REVOLD and DUMP_E is absolutely necessary, but the user will be warned if they are not the same.

The Ce.dat, Ke.dat, De.dat and g.dat Files¶

The two-temperature model (TTM) implementation in DL_POLY_5 allows specification of the following tabulated data:

Ce.dat - for electronic volumetric heat capacity - temperature dependent
Ke.dat - for metallic thermal conductivity - temperature dependent
De.dat - for non-metallic thermal diffusivity - temperature dependent
g.dat - for electron-phonon coupling constants - temperature dependent

Each file is free-formatted and only consist of two columns: the first column gives temperature in K, while the second gives electronic volumetric heat capacity in J m^-3 K^-1 (Ce.dat), thermal conductivity in W m^-1 K^-1 (Ke.dat), thermal diffusivity in m² s^-1 (De.dat) or the electron-phonon coupling constant in W m^-3 K^-1 (g.dat). These files are read by the subroutine ttm_table_read if the ttm cetab, ttm ketab, ttm detab and/or ttm gvar directives are included in the CONTROL file.

The HISTORY/HISTORF File¶

The HISTORY file is usually an output file (see Section The HISTORY File). However, upon specifying the replay option a HISTORY trajectory can be replayed and various observables recreated should there is enough information within to recover the specific observable (e.g. velocities supplied within HISTORY for regeneration of velocity autocorrelation data). Thus it can be used as input as well.

If the replay force option is used in CONTROL it enforce a new trajectory generation, i.e. generation of HISTORY, which necessitates the trajectory reading routine to default to reading a differently named, HISTORF (copy of HISTORY), basic trajectory file instead. However, using this option expects a full force evaluation driven by expectedly a different FIELD field file from the one used for the HISTORF generation. For more information do refer to CONTROL file directives and options (see Section The CONTROL File).

The SETEVB File¶

The file SETEVB is needed for EVB simulations. If this file is not found, the execution of DL_POLY_5 is aborted. See section Setting EVB calculations for a detailed explanation of the input parameters for EVB calculations.

The KPOINTS File¶

The KPOINTS file is needed for currents to be calculated. See section Currents. The file format is

record 1:
n     integer count of kpoints
record 2:
x     real     1st kpoint x-component
y     real     1st kpoint y-component
z     real     1st kpoint z-component
record 3:
x     real     2nd kpoint x-component
y     real     2nd kpoint y-component
z     real     2nd kpoint z-component
...
record n:
x     real     nth kpoint x-component
y     real     nth kpoint y-component
z     real     nth kpoint z-component

Key	Potential Type	Variables	Functional Form
harm	Harmonic	\(k\)	\(U(r) = \frac{1}{2}~k~(r_{i}-r_{i}^{t=0})^{2}\)
rhrm	Restraint	\(k\), \(r_{c}\)	\(U(r) = \frac{1}{2}~k~(r_{i}-r_{i}^{t=0})^{2}\) : \(\|r_{i}-r_{i}^{t=0}\| \le r_{c}\) \(U(r) = \frac{1}{2}~k~r_{c}^{2}+k~r_{c}(\|r_{i}-r_{i}^{t=0}\|-r_{c})~:~\|r_{i}-r_{i}^{t=0}\|>r_{c}\)
quar	Quartic	\(k\), \(k'\), \(k''\)	\(U(r) = \frac{k}{2}~(r_{i}-r_{i}^{t=0})^{2}+\frac{k'}{3}~(r_{i}-r_{i}^{t=0})^{3}\) \(+\frac{k''}{4}~(r_{i}-r_{i}^{t=0})^{4}\)

Key	Potential Type	Variables	Functional Form
ters	Tersoff	\(A\), \(a\), \(B\), \(b\), \(R\)	as shown in Section Tersoff Potentials
	(single)	\(S\), \(\beta\), \(\eta\), \(c\), \(d\), \(h\)
	(cross)	\(\chi\), \(\omega\), \(\delta\)
kihs	KIHS	\(A\), \(a\), \(B\), \(b\), \(R\), \(S\), \(\eta\), \(\delta\), \(c_{1}\), \(c_{2}\) \(c_{3}\), \(c_{4}\), \(c_{5}\), \(h\), \(\alpha\), \(\beta\)	Section Tersoff Potentials

Key	Potential Type	Variables	Functional Form
harm	Harmonic	\(k\), \(\theta_{0}\)	\(U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}\)
thrm	Truncated harmonic	\(k\), \(\theta_{0}\), \(\rho\)	\(U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}~\exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]\)
shrm	Screened harmonic	\(k\), \(\theta_{0}\), \(\rho_{1}\), \(\rho_{2}\)	\(U(\theta) = {k \over 2}~(\theta - \theta_{0})^{2}\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]\)
bvs1	Screened Vessal ¹⁴¹	\(k\) \(\theta_{0}\), \(\rho_{1}\), \(\rho_{2}\)	\(U(\theta) = {k \over 8(\theta_{0}-\pi)^{2}} \left\{ \left[ (\theta_{0} -\pi)^{2} -(\theta-\pi)^{2}\right]^{2}\right\} \times\) \(\exp[-(r_{ij}/\rho_{1} + r_{ik}/\rho_{2})]\)
bvs2	Truncated Vessal ¹²³	\(k\), \(\theta_{0}\), \(a\), \(\rho\)	\(U(\theta) = k~(\theta-\theta_{0})^{2}~\big[ \theta^a (\theta-\theta_{0})^{2} (\theta+\theta_{0}-2\pi){^{2}}\) \(+ {a \over 2} \pi^{a-1}(\theta_{0}-\pi)^{3}\big]~\exp[-(r_{ij}^{8} + r_{ik}^{8})/\rho^{8}]\)
hbnd	H-bond ⁷⁵	\(D_{hb}\), \(R_{hb}\)	\(U(\theta) = D_{hb}~\cos^{4}(\theta) \times\) \([5(R_{hb}/r_{jk})^{12} - 6(R_{hb}/r_{jk})^{10}]\)

Key	Potential Type	Variables	Functional Form \(^\ddagger\)
harm	Harmonic	\(k\), \(\phi_{0}\)	\(U(\phi) = {k \over 2}~(\phi - \phi_{0})^{2}\)
hcos	Harmonic cosine	\(k\), \(\phi_{0}\)	\(U(\phi) = {k \over 2}~(\cos(\phi) - \cos(\phi_{0}))^{2}\)
plan	Planar	\(A\)	\(U(\phi) = A~[1 - \cos(\phi)]\)

Key	Potential Type	Variables	Functional Form
elec	Electric Field	\(E_{x}\), \(E_{y}\), \(E_{z}\)	\(\underline{F} = q \; \underline{E}\)
oshr	Oscillating Shear	\(A\), \(n\)	\(\underline{F}_x = A \; \cos(2 n\pi \cdot z/L_{z})\)
shrx	Continuous Shear	\(A\), \(z_{0}\)	\(\underline{v}_x = \frac{A}{2} \frac{\|z\|}{z}~:~\|z\|~>~z_{0}\)
grav	Gravitational Field	\(G_{x}\), \(G_{y}\), \(G_{z}\)	\(\underline{F} = m \; \underline{G}\)
magn	Magnetic Field	\(H_{x}\), \(H_{y}\), \(H_{z}\)	\(\underline{F} = q \; (\underline{v} \times \underline{H})\)
sphr	Containing Sphere	\(A\), \(R_{0}\), \(n\), \(R_{\rm cut}\)	\(\underline{F} = A \; (R_{0} - r)^{-n}~:~r > R_{\rm cut}\)
zbnd	Repulsive Wall	\(A\), \(z_{0}\), \(f=\pm 1\)	\(\underline{F} = A \; (z_{0} - z) ~:~f \cdot z > f \cdot z_{0}\)
xpis	X-Piston	\(i_{gid}\), \(j_{gid}\), \(P_{k {\rm atm}}\)	\(\underline{F}_{x} = \frac{P \cdot \underline{\textrm{Area}}(\perp\textrm{X-}dir)}{\left[\sum_{k=i}^{j}m_{k}\right] \textrm{/} m_{k}} : \forall~k=i,..,j\)
zres	Molecule in HR Zone	\(i_{gid}\), \(j_{gid}\), \(A\), \(z_{mn}\), \(z_{mx}\)	\(\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z_{cm}-z_{mx}) & z_{cm} > z_{mx} \\A(z_{mn}-z_{cm}) & z_{cm} < z_{mn}\end{array} \right.\)
zrs-	HR Zone (push out)	\(i_{gid}\), \(j_{gid}\), \(A\), \(z_{mn}\), \(z_{mx}\)	\(\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z-z_{mx}) & z \ge \frac{z_{mx}+z_{mn}}{2} \\A(z_{mn}-z) & z < \frac{z_{mx}+z_{mn}}{2}\end{array} \right.\)
zrs+	HR Zone (pull in)	\(i_{gid}\), \(j_{gid}\), \(A\), \(z_{mn}\), \(z_{mx}\)	\(\underline{F}_{z} = \left\{ \begin{array} {l@{~:~}l}A(z-z_{mx}) & z > z_{mx} \\A(z_{mn}-z) & z < z_{mn}\end{array} \right.\)
osel	Osc. Electric Field	\(E_{x}\), \(E_{y}\), \(E_{z}\), \(\omega_{{\rm ps}^{-1}}\)	\(\underline{F} = q \; \underline{E} \; \sin(2 \pi \omega t)\)
ushr	umbrella sampling ¹³⁶	\(i^{A}_{gid}\), \(j^{A}_{gid}\), \(i^{B}_{gid}\), \(j^{B}_{gid}\), \(k\), \(R_{0}\)	\(U_{AB} = \frac{k}{2} (R_{AB}-R_{0})^{2}\)
	harm. constraint ⁵⁸