h-CMD Tutorial

This tutorial offers detailed, step-by-step guidance for performing h-CMD simulations [11] with IBI fitting using DL_POLY Quantum 2.1. Follow these instructions to efficiently set up and execute your h-CMD simulation for accurate IR spectra calculations.

Step 1: Setup h-CMD

To perform h-CMD simulations, additional keywords are added to the CONTROL file. This CONTROL file is used for performing the initial PIMD reference calculations, where the keyword pimd indicates the integrator to be used.

CONTROL

# h-CMD of LiTFSI + q-TIP4P/F WATER PIMD Reference

# pimd simulation with 32 beads
pimd pile nbeads=32 taut=0.05

# invoking qTIP4P/f water potential model
qtip4pf

# SYSTEM TARGETED TEMPERATURE (K) 
temperature  298.0

# INTEGRATION TIMESTEP (ps)
timestep 0.0002

# SIMULATION & EQUILIBRATION LENGTH
equil	       250000
steps          1000000
trajectory     251000 5000 1

# OUTPUT
print          5000

# h-CMD Options
fqcmd pot=15 points=998 outpts=480 bnd=5 rmax=2.5 ang=8 wrt=100

# STATISTICS
stats          1000

# SYSTEM CUTOFFS FOR ELECTROSTATICS AND NON-BONDED
cutoff         10.0
rvdw           10.0
ewald precision   1.0e-6
delr width      0.5	! for Verlet list

# EXECUTION TIME
job time       9999999.0
close time     20.0

finish

1.1 h-CMD Specification

The simulation is specified as part of an h-CMD simulation using the fqcmd keyword. When this keyword is added for systems with multiple molecule types (as defined in the FIELD file), DL POLY Quantum automatically treats the simulation as h-CMD.

Note

In Version 2.1, the software assumes the final molecule type (typically water) is treated at the f-QCMD level, while other molecule types are treated at the f-CMD level.

1.2 FQCMD Keyword Options in `CONTROL` File

The options following the fqcmd keyword determine how the h-CMD simulation works:

pot Option: Specifies the number of non-bonded pair interactions using a correction potential.
points Option: Specifies the number of grid points in the TABLE files.
bnd Option: Specifies the total number of unique bonds.
rmax Option: Defines the maximum value for the bond distribution function.
ang Option: Specifies the total number of unique angles.
wrt Option: Defines the number of timesteps between distribution function calculations.

Note

Each non-bonded pair has an associated TABLE file named PAIR__POT.TABLE. The first line of the table contains the atom names as defined in the FIELD file. During PIMD reference simulations, potential and force values in the table are ignored; however, pair names are used to calculate RDFs.

Note

The maximum value for angle distribution is set to π.
The maximum for RDFs is the non-bonded cutoff value in the CONTROL file.

1.3 Output Files

h-CMD simulations generate three additional output files:

AVERAGE_INTRA.D: Includes intramolecular bond distributions.
AVERAGE_ANGLE.D: Includes intramolecular angle distributions.
AVERAGE_RDF.D: Includes the RDFs for the system.

Note

The number of points in these output files is defined by the outpts option in the CONTROL file.

1.4 FQCMD File

An additional input file, the FQCMD file, is required for h-CMD simulations. For the Li-TFSI system, an example is shown here.

FQCMD

FQCMD params
TFSI
bonds 4
  harm 4
  harm 6
  harm 2
  harm 2
angles 7
  harm 4
  harm 2
  harm 4
  harm 6
  harm 6
  harm 2
  harm 1
finish
Li-ion
bonds 0
angles 0
finish
H2O
bonds 1
  quar 2
angles 1
  harm 1
finish
CLOSE

Note

FQCMD file includes details about the bonds and angles for each molecule type in the same order as in the FIELD file:
bonds option specifies the number of unique bond types for a molecule and their occurrences.
angles option specifies the number of unique angle types and their occurrences.

Step 2: PIMD Reference Simulation

For h-CMD simulation, at first, a Reference PIMD Simulation must be carried out to compute the reference intra- and inter-molecular distribution functions. These distributions serve as the benchmark that the effective potentials will aim to replicate.

Note

It is often advantageous to perform multiple PIMD reference simulations starting from different initial configurations. This approach improves the statistical accuracy of the calculated distribution functions.

2.1 Prepare the Reference PIMD Simulation:

Set up the input files (CONFIG, CONTROL, FIELD, FQCMD, PAIR__POT.TABLE) of a PIMD simulation for DL_POLY Quantum.

2.2 Run the PIMD Simulation:

Execute the simulation to gather sufficient sampling of the reference system.
Simulation conditions typically include: - Canonical ensemble (NVT) with a quantum thermostat. - Number of beads (n_beads) to achieve convergence in quantum effects. - Adequate simulation time to ensure good sampling of the distributions.

2.3 Calculate Reference Distributions:

Average_Intra.d
Average_Angle.d
Average_RDF.d

Step 3: Classical Simulations (`iter0`)

To begin the IBI process, an iteration 0 is needed where the dynamics are performed classically using the base force field potential. This is equivalent to performing h-CMD with the correction potentials set to zero.

Most of the setup work done for the PIMD reference simulations can be reused here:

FQCMD file: No changes are required.
PAIR__POT.TABLE files: All values of the potential and force must be set to zero so that the dynamics are not altered from the base potential.
CONTROL file: This file remains very similar to the one used previously, with alter in integrator as:

.
integrator velocity verlet
ensemble nvt nhc 0.05 1 3
.
.
.
fqcmd  pot=15  points=998  outpts=480  bnd=5  rmax=2.5  ang=8  wrt=100
.

Note

Changes made in CONTROL file:

The integrator is changed to invoke the velocity Verlet algorithm.
The NVT ensemble is maintained through a NHC thermostat.
The simulation length can be shortened since RDFs are calculated using force sampling in classical dynamics.

This step performs classical dynamics without correction, providing a baseline, which serve as the starting point for subsequent iterations in the IBI process.

Step 4: IBI Iteration Workflow

After the PIMD reference and iteration 0 simulations are complete, the IBI process can begin to refine the effective ineraction potentials. Between each IBI iteration, the correction potentials are updated. This update is done using a separate Python script.

A small portion of the script, designed to work on high-performance computing clusters with SLURM schedulers, requires modification to detail the system of interest. Most of the detail parallels the FQCMD and CONTROL files, with some additional information required. For example, the pairtype array specifies the non-bonded pairs for which the correction potential is calculated.

Snippet of the IBI code detailing the non-bonded pairs and intramolecular details for the Li-TFSI system.

npair=15
cutoff=10.0
nrdf=480
pairtype=["C OW",
        "C HW",
        "N OW",
        "N HW",
        "O OW",
        "O HW",
        "F OW",
        "F HW",
        "S OW",
        "S HW",
        "Li OW",
        "Li HW",
        "OW OW",
        "OW HW",
        "HW HW"
          ]
nmols=2
nbonds=[4,1]
totbonds=sum(nbonds)
bonds=[
        ["harm",4],
        ["harm",6],
        ["harm",2],
        ["harm",2],
        ["quar",2]
        ]
nangs=[7,1]
totangs=sum(nangs)
angs=[
        ["harm",4],
        ["harm",2],
        ["harm",4],
        ["harm",6],
        ["harm",6],
        ["harm",2],
        ["harm",1],
        ["harm",1],
        ]

Note

pairtype: Holds the list of non-bonded pairs for correction potential calculations. This list should match the TABLE files.
nbonds: Specifies the number of unique bonds in each molecule type.
nangs: Specifies the number of unique angles in each molecule type.
bonds and angles: Store information about each unique bond and angle type, including occurrences.

Note

To execute the IBI script, two additional files are required:

OLD_FIELD: The FIELD file from iteration 0, renamed to OLD_FIELD. This file contains bond and angle force field parameters.
PREVIOUS_POTENTIAL.DAT: Combines the correction potentials for all pairs from the previous iteration. For iteration 0, these values will be zero, but the file is still necessary.

Note

Once the IBI script is executed, the following new files are created:

NEW_FIELD: Contains updated intramolecular force field parameters.
NEW_POTENTIAL.DAT: Includes updated intermolecular correction potentials.
PAIR__POT.TABLE files: Generated for each pair type, containing new potentials and forces.

The updated files are used to run simulations for the next iteration. Starting with iteration 1, the dynamics are now classical dynamics under the correction potentials, mimicking PI-based dynamics. The process is repeated until the distribution functions sufficiently match those of the PIMD reference simulations.

Note

A set of scripts is provided to automate the setup and execution of multiple IBI iterations, designed to work on high performance computing clusters with SLURM schedulers, with the jobs for later iterations requiring the previous ones to finish before running

Here is workflow of for single interation IBI fitting.

Step 4.1: IBI Iteration Setup:

Prepare the iter-setup.sh script to automate iteration setup: Each iteration runs 50 NVT simulations with 50 ps of equilibration and 50 ps of RDF gathering.
Edit the simulation parameters in the iterCONTROL file. Example settings include temperature, time steps, and RDF parameters.

Update the iteration number in the iter-setup.sh script:

i=<current_iteration_number>

iter-setup.sh

#!/bin/bash
let i=1
let j=i-1
let ntraj=50

if [ -d iter${i} ]; then
  echo iter${i} exists, change iteration number at first!!
  exit 1
elif [ ! -d iter${j} ]; then
  echo iter${j} does not exist, change iteration number at first!!
  exit 1
else
  echo "running simulation for iter${i}"
  echo "copying necessary files from iter${j}"
fi
mkdir iter${i}
mkdir iter${i}/run
cp base_scripts/run_dlpoly.sh iter${i}/run/
cp base_scripts/run_avg.sh iter${i}/run/
cp base_scripts/run_IBI.sh iter${i}/
cd iter${i}/run
mkdir base_input
cp ../../FQCMD base_input/
cp ../../iterCONTROL base_input/CONTROL
cp ../../avgCONTROL CONTROL
cp ../../iter${j}/*.table base_input/ 
cp ../../iter${j}/new_FIELD base_input/FIELD
cp ../../iter${j}/new_FIELD ../old_FIELD
cp ../../iter${j}/New_Potential.dat ../Previous_Potential.dat

for k in $(seq -f "%03g" 1 ${ntraj})
do
   mkdir traj${k}
   cp base_input/* traj${k}/
   cp run_dlpoly.sh traj${k}/
   cp ../../ref/run/traj${k}/CONFIG traj${k}/ 
   cd traj${k}
      sbatch -J itr${i}-${k} run_dlpoly.sh
   cd ../
done
cd ../..

Step 4.2: Running Simulations

Execute the iter-setup.sh script:

$ bash iter-setup.sh

This will start 50 classical-like dynamics under the correction potentials and generate RDFs for each trajectory. The CONTROL file looks like:

.
integrator velocity verlet
ensemble nvt nhc 0.05 1 3
.
.
.
fqcmd  pot=15  points=998  outpts=480  bnd=5  rmax=2.5  ang=8  wrt=100
.

Step 4.3: Averaging RDFs

Navigate to the simulation directory after the completion of simulation where the RDFs are stored:
```
$ cd iter<current_iteration_number>/run
```
Use the rdfAvg program (provided in utility folder) to compute the average RDFs from all trajectories:
```
$ rdfAvg
```
Note

If you change the number of trajecties or RDF points, update the avgCONTROL file to reflect these changes.

Step 4.4: Updating Potentials

Return to the iteration directory; Move up one directory to the main iteration forder:
```
$ cd ..
```
Ensure the following required files exist.
- OLD_FIELD (FIELD file from the previous iteration)
- PREVIOUS_POTENTIAL.DAT (potential corrections from the previous iteration)
Run the IBI Python script to generate updated potentials:
```
$ python IBI.py
```
The script generates the updated potential corrections:
- NEW_FIELD
- NEW_POTENTIAL.DAT.dat
- A series of PAIR__POT>TABLE files

Step 4.5: Evaluating Convergence

Compare the RDFs of the reference and the currect iteration.
If the RDFs do not match, repeat the process by starting a new iteration by updating the iteration number in iter-setup.sh and rerun the steps above.

Note

Production Simulations

After the correction potentials have converged, these corrected potentials are utilized for subsequent simulations of the same system. The CONTROL files for these simulations should retain the same options under the fqcmd keyword as in earlier simulations. However, the value of the wrt option can be increased to minimize computation time, as the distribution functions—though no longer needed—are still being calculated.

Key Files and Their Roles

This table describes the key files used during the IBI process and their roles.

File Name	Description
`iter-setup.sh`	Bash script to automate simulation setup and execution for the current iteration.
`iterCONTROL`	Configuration file for controlling simulation parameters.
`avgCONTROL`	Configuration file for averaging RDFs.
`OLD_FIELD`	Field file from the previous iteration.
`NEW_FIELD`	Updated field file generated by the IBI script.
`PREVIOUS_POTENTIAL.DAT`	Potential corrections from the previous iteration.
`NEW_POTENTIAL.DAT`	Updated potential corrections.
`PAIR_<i>_POT.TABLE`	Individual tables representing updated potentials for specific interactions.

h-CMD Tutorial

Step 1: Setup h-CMD

1.1 h-CMD Specification

1.2 FQCMD Keyword Options in CONTROL File

1.3 Output Files

1.4 FQCMD File

Step 2: PIMD Reference Simulation

Step 3: Classical Simulations (iter0)

Step 4: IBI Iteration Workflow

Step 4.1: IBI Iteration Setup:

Step 4.2: Running Simulations

Step 4.3: Averaging RDFs

Step 4.4: Updating Potentials

Step 4.5: Evaluating Convergence

Key Files and Their Roles

1.2 FQCMD Keyword Options in `CONTROL` File

Step 3: Classical Simulations (`iter0`)