Linear H10 chain¶
1S-COGEF¶
We first pull on the outer atoms of a linear H10 chain kept together by Morse potentials. We let external forces act on indices 0 and 8 instead of 0 and 9 in order to break the symmetry and to break as single bond only
# creates: 1scogef_problem.png
from ase.atoms import Atoms
from ase.optimize import FIRE
from ase.calculators.morse import MorsePotential
from cogef import COGEF1D
fmax = 0.05
image = Atoms('H10', positions=[(i, 0, 0) for i in range(10)])
image.calc = MorsePotential()
FIRE(image, logfile=None).run(fmax=fmax)
cogef = COGEF1D(0, 8, optimizer=FIRE, optimizer_logfile=None, fmax=fmax)
stepsize = 0.03
if len(cogef.images) == 0:
cogef.images = [image]
steps = 35
cogef.move(stepsize, steps)
The transition state on the 1S-COGEF path can be wrong for long molecules as indicated by a discontinuity in energy [Brügner2018] in the figure above. The figure was created by
# creates: 1scogef_problem.png
import matplotlib.pyplot as plt
from cogef import COGEF1D
cogef = COGEF1D(0, 8)
energies = cogef.get_energies()
distances = cogef.get_distances()
plt.plot(distances, energies)
plt.xlabel(r'd [$\rm\AA$]')
plt.ylabel('U [eV]')
plt.savefig('1scogef_problem.png')
Let us obtain the corresponding electronic activation energies and rate constants, and save them:
# creates: 1scogef_error.png
from ase.io.jsonio import write_json
from cogef import COGEF, Dissociation
cogef = COGEF(0, 8)
T = 300 # Temperature [K]
P = 101325. # Pressure [Pa]
force_min = 0.1
force_max = 4.
force_step = 0.1
diss = Dissociation(cogef, force_unit='nN')
rates, forces = diss.get_rate_constants(
T, P, force_max, force_min, force_step, method='electronic')
barriers = []
for f_ext in forces:
barriers.append(diss.electronic_energy_barrier(f_ext))
with open('1scogef_barriers.json', 'w') as f:
write_json(f, {'forces': forces, 'barriers': barriers})
3S-COGEF¶
The 3S-COGEF method is based on a three-segmented path through a two-dimenisonal energy landscape. It can be used to get better transition states, activation energies and rate constants for the dissociation reaction in dependence of the external force. The following picture shows that 1S-COGEF overestimates the energy barrier. The script used
# creates: 1scogef_error.png
import matplotlib.pyplot as plt
from ase.units import J, m
from ase.io.jsonio import read_json
from ase.calculators.morse import MorsePotential
from cogef import COGEF
from cogef.cogef2d.fixed2d import FixedForce2D
cogef = COGEF(0, 8)
cogef2 = COGEF(0, 1, optimizer_logfile=None)
def initialize(image):
image.calc = MorsePotential()
return image
cogef3s = FixedForce2D(cogef, cogef2, initialize)
stepsize = 0.01
cogef3s.find_barrier(stepsize, range(1, 20))
forces, barriers = cogef3s.collect_barriers()
plt.plot(forces / J * m * 1e9, barriers, 'o-', label='3S-COGEF')
with open('1scogef_barriers.json') as f:
dct = read_json(f)
plt.plot(dct['forces'], dct['barriers'], '-', label='1S-COGEF')
plt.legend()
plt.xlabel('force [nN]')
plt.ylabel('barrier [eV]')
# plt.show()
plt.savefig('1scogef_error.png')
XXX fix following
and the error decreases with force compared to 3S-COGEF.
The class COGEF2D was used to apply the 3S-COGEF method. The first segment
of the 3S-COGEF path is the 1S-COGEF path up to bond breaking.
This path can also be obtained with COGEF2D, but we will
import it as we have it already. The breaking bond with length b must be set
([0, 1]) in addition to the end-to-end distance d defined like before.
Finally we choose between two different methods for variation of b.
Here, we want to fix the external force and not d which is often simpler to
use for chain molecules:
from cogef import COGEF2D, Dissociation2d
from os.path import join
def initialize2d(image, directory, imagenum, new_opt, get_filename):
"""Initialize the image or return the trajectory name.
"""
if initialize(image, -1, new_opt, get_filename) == 'Finished':
return 'Finished'
if get_filename:
return join(directory, 'cogef2d') + str(imagenum) + '.traj'
cogef = COGEF2D([0, 8], [0, 1], 'cogef1s.traj', optimizer=FIRE, fmax=fmax,
fix_force=True)
The discontinuity was already identified and we now the last image number, where the configuration of the 1S-COGEF path has an intact bond. This should be set to prevent error messages in the following:
cogef.set_last_intact_bond_image(26)
The second and the third segments of the 3S-COGEF method (maximum and minimum with respect to variation of b) should be obtained as far as we need them to obtain rate constants in the force range defined earlier:
energy_tolerance = 0.01
diss = Dissociation2d(cogef, force_unit='nN',
transition_into_two_fragments=True)
def get_rates():
# Use some global variables
while True:
try:
rates = diss.get_rate_constants(
T, P, force_max, force_min, force_step, method='electronic')
break
except ValueError:
assert diss.error in [1, 2]
# More images must be calculated
if len(diss.needed_max_images) > 0:
I = diss.needed_max_images[-1]
diss.clean_needed_images_list()
is_max = True
elif len(diss.needed_min_images) > 0:
I = diss.needed_min_images[-1]
diss.clean_needed_images_list()
is_max = False
else:
I = None
if (I is None or len(cogef.images) <= I):
# Calculate 1S-COGEF path further
cogef.pull(stepsize, 1, initialize, 'cogef1s.traj')
else:
if is_max:
cogef.calc_maximum_curve(
[I], stepsize, initialize2d=initialize2d,
max_trajectory='maximum.traj',
breakdirectory='pull',
min_trajectory='minimum.traj',
only_minimum_curve=not is_max,
energy_tolerance=energy_tolerance)
else:
cogef.calc_maximum_curve(
[I], stepsize, initialize2d=initialize2d,
max_trajectory='maximum.traj',
breakdirectory='pull',
min_trajectory='minimum.traj',
and_minimum_curve=not is_max,
energy_tolerance=energy_tolerance)
return rates
rates_3scogef, forces_3scogef = get_rates()
This will return an error message: “ValueError: Cannot find energy maximum. Maximum image number is too small or energy_tolerance is too large or too small. If you want to increase energy_tolerance, remove pull_2/pull.traj first.” Look into pull_2/pull.traj and you will see that the energy maximum cannot be identified within 0.01 eV energy tolerance. We do not want to decrease this tolerance, so let the program allow to calculate more images in b-direction:
cogef.max_image_number = 30
rates_3scogef, forces_3scogef = get_rates()
barriers_3scogef = []
for f_ext in forces_3scogef:
barriers_3scogef.append(diss.electronic_energy_barrier(f_ext))
Now we can obtain the picture from above for comparison of 1S- and 3S-COGEF:
plt.figure(1)
plt.plot(forces_1scogef, rates_1scogef, '--', color='blue', label='1S-COGEF',
lw=2)
plt.plot(forces_3scogef, rates_3scogef, color='blue', label='3S-COGEF', lw=2)
plt.xlabel('F [nN]')
plt.ylabel('Rate constant k [1/s]', color='blue')
plt.yscale('log')
ax1 = plt.gca()
for tl in ax1.get_yticklabels():
tl.set_color('blue')
plt.legend(loc=5)
ax2 = plt.twinx()
ax2.plot(forces_1scogef, barriers_1scogef, '--', color='green', lw=2)
ax2.plot(forces_3scogef, barriers_3scogef, color='green', lw=2)
try:
ddagger = unichr(int('2021', 16)) # Python 2
except NameError:
ddagger = chr(int('2021', 16)) # Python 3
plt.ylabel('Activation energy $\\Delta U^' + ddagger + '$ [eV]',
color='green')
for tl in ax2.get_yticklabels():
tl.set_color('green')
plt.xlim(1, 4)
plt.savefig('1scogef_error.png')
Rate constants can be saved with:
diss.save_rate_constants(
T, P, force_max, force_min, force_step, method='electronic',
fileout='bond1.dat')
Our python file for saving rate constants, and the trajectory file of the 1S-COGEF path can be copied in new directories to obtain rate constants for breaking of other bonds of the same molecule. Only the definition of b must be changed. The python file in one of the new directory may look like the following where we have to add the imports and the get_rates function:
cogef = COGEF2D([0, 8], [1, 2], 'cogef1s.traj', optimizer=FIRE, fmax=fmax,
fix_force=True)
cogef.set_last_intact_bond_image(26)
cogef.max_image_number = 30
energy_tolerance = 0.01
diss = Dissociation2d(cogef, force_unit='nN',
transition_into_two_fragments=True)
get_rates()
diss.save_rate_constants(
T, P, force_max, force_min, force_step, method='electronic',
fileout='bond2.dat')
Now, create a new directory where we add all the dat files and the following script. By running the script we can obtain rupture forces by taking multiple dissociation reactions into account. In the following we take only two (bond1 and bond2):
from numpy import array
from cogef import load_rate_constants, Minima, probability_density
from cogef import rupture_force_and_uncertainty_from_dpdf
loading_rate = 1000. # nN/s as defined in the Dissociation2d objects
# used for saving rate constants and associated forces
rates_bond1, forces_bond1 = load_rate_constants('bond1.dat')
rates_bond2, forces_bond2 = load_rate_constants('bond2.dat')
forces = forces_bond1
assert forces == forces_bond2
minima = Minima()
minima.add_destination(rates_bond1)
minima.add_destination(rates_bond2)
dpdf = probability_density(minima, forces, loading_rate)
print('Final probabilities:')
print(minima.prob)
f_rup, f_err = rupture_force_and_uncertainty_from_dpdf(
-array(dpdf[0]), forces)
print('Rupture force:')
print(str(f_rup) + 'nN')
print('Standard deviation:')
print(str(f_err) + 'nN')
In order to be sure that the force range is large enough, it is helpful to
plot the peak in -dpdf[0] over forces and check whether the final
probability of the reactant state minima.prob[0] is almost zero.