|
|
| Line 1: |
Line 1: |
| The {{FILE|IRCCAR}} file defines an intrinsic reaction coordinate (IRC), similar to the [[Intrinsic-reaction-coordinate calculations]] using {{TAG|IBRION}}=40 for static calculations. However, this IRC is followed during a biased molecular dynamics (MD) simulation, similar to the [[Blue moon ensemble]] method. | | The {{FILE|IRCCAR}} file defines an intrinsic reaction coordinate (IRC), similar to the [[Intrinsic-reaction-coordinate calculations]] using {{TAG|IBRION}}=40 for static calculations. However, this IRC is followed during a biased molecular dynamics (MD) simulation, similar to the [[Blue moon ensemble]] method. The Structure of the file is: |
| | |
| Note, this how to follows the structure of [[https://www.vasp.at/tutorials/latest/transition_states/part3/]].
| |
| | |
| == Defining IRCCAR ==
| |
| | |
| Similarly to the [[Intrinsic-reaction-coordinate calculations|IRC]], the reaction coordinate can be defined following a ???
| |
| | |
| === Part 1 ===
| |
| | |
| This generates an {{FILE|XDATCAR}} file. Using the <code>tellcoord.py</code> script, you can compute the values of the basis coordinates for the negative branch of the IRC:
| |
| | |
| <div class="toccolours mw-customtoggle-script">'''Click to show <code>tellcoord.py</code> script'''</div>
| |
| <div class="mw-collapsible mw-collapsed" id="mw-customcollapsible-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() | | N |
| xdat.readFile(fname, blocksize)
| | a1 b1 |
| numofatoms = xdat.numofatoms
| | a2 b2 |
| numconfig = xdat.numconfig
| | a3 b3 |
| katoms = xdat.atoms
| | ... |
| ntypes = len(katoms)
| | aN bN |
| 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])
| |
| </div>
| |
| | |
| 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 <code>m</code> and the <code>p</code> 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 {{FILE|IRCCAR}} file itself is created using the <code>ircprepare3.py</code> script.
| |
| | |
| <div class="toccolours mw-customtoggle-script">'''Click to show <code>ircprepare3.py</code> script'''</div>
| |
| <div class="mw-collapsible mw-collapsed" id="mw-customcollapsible-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)
| |
| | |
| </div>
| |
| | |
| It generates <code>7</code> equidistant points (within a tolerance <code>1e-3</code> 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 <code>1</code> 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 {{FILE|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 {{FILE|ICONST}} file is required. Instead of a simple sum <code>S</code> like in the [[Blue moon ensemble]], you can use <code>IS</code>, defining the weight using the <code>Suggested LAMBDA</code> 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 [[:Category:Biased_molecular_dynamics|step-shaped potential]], the following {{FILE|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 {{FILE|INCAR}} file should also have the following tags, adjusting numerical values to the problem at hand):
| |
|
| |
|
| FBIAS_R0 = 0.2
| | where N is the number of points in the file, a<sub>i</sub> is the first basis???, b<sub>i</sub> is the second basis. |
| FBIAS_D = 5.0
| |
| FBIAS_A = 5
| |
|
| |
|
| | There is a guide on how to use this tag in the How-to ADD LINK. |
|
| |
|
| <!--[[Category:Files]][[Category:Input files]][[Category:Advanced molecular-dynamics sampling]][[Category:Biased molecular dynamics]][[Category:Transition states]]--> | | <!--[[Category:Files]][[Category:Input files]][[Category:Advanced molecular-dynamics sampling]][[Category:Biased molecular dynamics]][[Category:Transition states]]--> |
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. The Structure of the file is:
N
a1 b1
a2 b2
a3 b3
...
aN bN
where N is the number of points in the file, ai is the first basis???, bi is the second basis.
There is a guide on how to use this tag in the How-to ADD LINK.
