IRCCAR
The IRCCAR file defines an intrinsic reaction coordinate (IRC), similar to the Intrinsic-reaction-coordinate calculations using IBRION=40 for static calculations. However, this IRC is followed during a biased molecular dynamics (MD) simulation, similar to the Blue moon ensemble method.
Note, this how to follows the structure of [[1]].
Defining IRCCAR
Similarly to the IRC, the reaction coordinate can be defined following a ???
Part 1
This generates an XDATCAR file. Using the 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])
you can compute the values of the basis coordinates for the negative branch of the IRC:
python3 tellcoord.py m/XDATCAR > coords_m.dat
and the positive branch of the IRC:
python3 tellcoord.py p/XDATCAR > coords_p.dat
Part 2
The values of the two basis coordinates can then be separated:
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
Part 3
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
Part 4
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
Part 5
In this final step, the IRCCAR file itself is created using the cv_all.dat 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 numerical values to the problem at hand):
FBIAS_R0 = 0.2 FBIAS_D = 5.0 FBIAS_A = 5