IRCCAR

The IRCCAR file defines a discretized transformation path. Usually this is taken from an intrinsic reaction coordinate (IRC) - using IBRION=40, projected onto a small set of internal coordinates. It is required to define the path-based coordinates IS and IZ in ICONST .

However, this path is followed during a biased molecular dynamics (MD) simulation, similar to the Blue moon ensemble method, metadynamics, and other advanced MD aproaches . The structure of the file is:

M
chi_1(1) chi_2(1) ... chi(r)(1)
chi_1(2) chi_2(2) ... chi(r)(2)
...
chi_1(M) chi_2(M) ... chi(r)(M)

where M is the number of points in the file and chi_j(i) is [math]\displaystyle{ \tilde{\chi}_k(i) }[/math], the discretized transformation path in terms of internal coordinates, taking values between 0 and 1.

Using IRCCAR

In the context of following a reaction pathway, using IS and IZ represent one of the few ways to define collective variables for investigating chemical reactions. The following example covers defining IRCCAR and the corresponding ICONST file.

Defining IRCCAR

After performing an IRC calculation, a discretized transformation path can be defined using the XDATCAR files.

1. Extracting basis coordinates from XDATCAR

Using the tellcoord.py script, you can compute the values of the basis coordinates making up the path for the negative branch of the IRC:

Click to show tellcoord.py script

#!/usr/bin/python3

#from Numeric import *
from numpy import *
import sys
import string
import abstracts


def same_cell(x, y):
    r = 0.0
    for i in range(len(x)):
        for j in range(len(x[0])):
            while (x[i][j]-y[i][j]) > 0.5:
                x[i][j] -= 1
            while (x[i][j]-y[i][j]) <= -0.5:
                x[i][j] += 1
    return x

fname = sys.argv[1]

blocksize=1
if len(sys.argv)>2:
    blocksize = int(sys.argv[2])

xdat = abstracts.Xdatcarreader()
xdat.readFile(fname, blocksize)
numofatoms = xdat.numofatoms
numconfig = xdat.numconfig
katoms = xdat.atoms
ntypes = len(katoms)
katoms_kumul = 1*katoms
katoms_range = {}
flags = xdat.tags_atoms

for i in range(len(xdat.data)):
    if i > 0:
        xdat.data[i] = same_cell(array(xdat.data[i]), array(xdat.data[i-1]))
    coords_c = dot(xdat.data[i], xdat.lattmat[0])
    coord = abstracts.Coordinates()
    coord.getDefinition( coords_c, xdat.lattmat[0], "IRCCAR")
    coord.detect2(coords_c, xdat.lattmat[0])

python3 tellcoord.py m/XDATCAR > coords_m.dat

and the positive branch of the IRC:

python3 tellcoord.py p/XDATCAR > coords_p.dat

2. Separating the basis coordinates

The values of the two basis coordinates can then be separated from one another to generate four files:

awk '{if (NR%2==1) {print $0}}' coords_m.dat > cv1_m.dat
awk '{if (NR%2==1) {print $0}}' coords_p.dat > cv1_p.dat
awk '{if (NR%2==0) {print $0}}' coords_m.dat > cv2_m.dat
awk '{if (NR%2==0) {print $0}}' coords_p.dat > cv2_p.dat

3. Reordering and connecting each basis coordinate

Then merge the data for the first basis coordinate from the m and the p branches, reversing the order of the data for the negative branch.

tac cv1_m.dat > cv1.dat
cat cv1_p.dat >> cv1.dat

In a reaction pathway, this would result in the transition state (TS) to reactant (R) changing to be R to TS in the negative branch, which then seemlessly links TS to product (P) in the positive branch. Repeat this for the second basis coordinate:

tac cv2_m.dat > cv2.dat
cat cv2_p.dat >> cv2.dat

4. Combining the basis coordinates

Merge the resulting cv1.dat and cv2.dat into a single file such that each basis coordinate is presented in a single column:

paste cv1.dat cv2.dat > cv_all.dat

5. Creating IRCCAR

In this final step, the IRCCAR file itself is created using the ircprepare3.py script.

Click to show ircprepare3.py script

#!/usr/bin/python3.6

from numpy import *
import string
import sys 

def find_path(irc, np, incr, mytiny):
    ltest = 1
    dx = []
    irc2 = []

    indx = 0
    dx.append(indx+1)
    irc2.append(irc[indx])

    for i in range(1, np):
        indx, ltest_ = find_nearest(irc, indx, incr, mytiny)
        ltest *= ltest_
        irc2.append(irc[indx])
        dx.append(indx+1)
    return irc2, ltest, dx


def find_nearest(irc, istart, x, mytiny):
    ltest = 0
    target = 0
    mindist = abs(x)
    for i in range(istart, len(irc)):
        xx = (array(irc[i])-array(irc[istart]))
        xx = sum(xx*xx)**0.5
        dist = abs(xx-x)
        if dist < mindist:
            mindist = dist
            target = i
    if mindist < mytiny:
        ltest = 1
    return target, ltest


def give_alls(numpoints, istart, irc):
    s = 0.
    ilength = numpoints-istart
    alls = zeros(ilength, float)
    for i in range(istart+1, numpoints):
        count = i-istart
        ds = irc[i]-irc[i-1]
        ds = sum(ds*ds)**0.5
        s += ds
        alls[count] = s
    return alls, s


f  = sys.argv[1]
f  = open(f, 'r')

np = int(sys.argv[2])
npExtra = int(sys.argv[3])
mytiny = float(sys.argv[4])


irc = []


for line in f.readlines():
    line = line.split()
    num = len(line)
    r = line[:num]
    for j in range(len(r)):
        r[j] = float(r[j])
    irc.append(r)
f.close()

irc = array(irc)

numpoints = len(irc)
numcoords = len(irc[0])


alls, s = give_alls(numpoints, 0, irc)

print("Total length of path:", s)
incr = s/(np-1)

print("Estimated increment:", incr)

irc2, ltest, dx = find_path(irc, np, incr, mytiny)

while ltest != 1:
    incr *= 0.99
    irc2, ltest, dx = find_path(irc, np, incr, mytiny)


extra0 = []
extraN = []

for i in range(npExtra):
    dummy = irc2[0] + (npExtra-i)*(irc2[0] - irc2[1])
    extra0.append(dummy)

    dummy = irc2[-1] + (i+1)*(irc2[-1] - irc2[-2])
    extraN.append(dummy)


print("Increment in discretized path", incr)

f = 'IRCCAR'
f = open(f, 'w')
rrr = str(np+2*npExtra)
f.write(rrr+'\n')


for i in range(npExtra):
    rrr = 
    for j in range(numcoords):
        rrr += ' '+'%.6f' % (extra0[i][j])
    f.write(rrr+'\n')


for i in range(len(irc2)):
    rrr = 
    for j in range(numcoords):
        rrr += ' '+'%.6f' % (irc2[i][j])
    f.write(rrr+'\n')

rrr = 
for i in range(npExtra):
    rrr = 
    for j in range(numcoords):
        rrr += ' '+'%.6f' % (extraN[i][j])
    f.write(rrr+'\n')


irc2 = array(irc2)

for i in range(1, len(irc2)):
    dd = sum((irc2[i]-irc2[i-1])**2)**0.5
    print(dd)

av = 0.
for i in range(1, np):
    dist = irc2[i]-irc2[i-1]
    dist = sum(dist*dist)
    av += dist
av /= (np-1)
av = 1./av
print("")
print('Suggested LAMBDA:', av)

print(dx)

It generates 7 equidistant points (within a tolerance 1e-3 along the IRC path. There are also two additional points, one at the start and end of the sequence by linear extrapolation (if more than 1 extra point should be added to each side, replace 1 by another number). This creates well-behaved histograms for the stable states.

python3 ircprepare3.py cv_all.dat 7 1 1e-3

This will generate the following output:

Total length of path: 2.7295405195612816
Estimated increment: 0.4549234199268803
Increment in discretized path 0.44141173943163203
0.4412554263534783
0.4409102031186374
0.44160482930366773
0.4406936438932761
0.4416511471372773
0.4413125503811718

Suggested LAMBDA: 5.136342611992123
[1, 1154, 1567, 1799, 1956, 2164, 3275]

and an IRCCAR file:

9
 1.490808 3.877537
 1.487941 3.436290
 1.485073 2.995044
 1.502921 2.554496
 1.782918 2.213004
 2.087445 1.894454
 2.373883 1.558286
 2.781098 1.388190
 3.188313 1.218095

Defining ICONST

To run the biased MD, an ICONST file is required. Instead of a simple sum S like in the Blue moon ensemble, you can use IS, defining the weight using the Suggested LAMBDA from above:

R 3 25 0
R 3 2 0
IS 5.136342611992123 5.136342611992123 0

If you want to define also the IZ coordinate and use it to prevent side reaction by imposing a step-shaped potential, the following ICONST should be used:

R 3 25 0
R 3 2 0
IS 5.136342611992123 5.136342611992123 0
IZ 5.136342611992123 5.136342611992123 4

The INCAR file should also have the following tags, adjusting the values to your problem, depending on how you want to define the step-function:

FBIAS_R0 = 0.2
FBIAS_D = 5.0
FBIAS_A = 5

Related tags and articles

ICONST, FBIAS_R0, FBIAS_D, FBIAS_A