AuAg2 chain

We want to simulate external forces by the application of geometry constraints. This “COnstraint GEometry to simulate Forces” is known as the COGEF method [Beyer2000].

1S-COGEF

The standard COGEF method constrains two atoms usually by their distance and relaxes all other degrees of freedom.

We first relax a linear Au-Ag-Ag chain with EMT:

# creates: 1scogef.png
from ase import io
from ase.atoms import Atoms
from ase.optimize import FIRE
from ase.calculators.emt import EMT

fmax = 0.01

fname = 'AuAg2.traj'
try:
    image = io.read(fname)
except FileNotFoundError:
    image = Atoms('AuAgAg', positions=((-1, 0, 0), (0, 0, 0), (1, 0, 0)))
    image.set_calculator(EMT())
    FIRE(image).run(fmax=fmax)
    image.write(fname)

Now we pull on the ends (atom indices 0, 2) of the chain:

# creates: 1scogef.png
from ase import io
from ase.optimize import FIRE
from ase.calculators.emt import EMT

from cogef import COGEF

image = io.read('AuAg2.traj')
image.calc = EMT()

fmax = 0.01
cogef = COGEF(0, 2, optimizer=FIRE, optimizer_logfile=None, fmax=fmax)

if not len(cogef.images):  # calculation was done already
    cogef.images = [image]
    
    stepsize = 0.1
    steps = 30
    cogef.pull(stepsize, steps)

which creates the directory cogef_0_2 and writes the file cogef_0_2/cogef.traj containing the corresponding trajectory. Observing the trajectory shows that eventually the Ag-Ag bond breaks. The maximum force needed for separation can be either obtained from a numerical derivative of the energy or directly from the forces:

from ase.units import J, m
from cogef import COGEF

f = m / J * 1e9  # eV / A -> nN

cogef = COGEF(0, 2)  # this reads cogef_0_2/cogef.traj if it is there

print('Maximal force: {0:4.2f} nN (from energies),'.format(
    cogef.get_maximum_force(method='use_energies') * f) +
    ' {0:4.2f} nN (from forces)'.format(
        cogef.get_maximum_force(method='use_forces') * f))

Dissociation

../../_images/1scogef.png

In reality, bonds can not withstand much smaller forces due to temperature effects. We may to include finite temperature into consideration. This is done using the Dissociation object that allows to obtain force-dependent probability distributions \(dp/dF\) as shown in the figure above.

The figure was created with:

# creates: 1scogef.png
import pylab as plt

from cogef import COGEF
from cogef import Dissociation
from cogef.units import nN

cogef = COGEF(0, 2)

alpha = 10  # [nN/s]
T = 300      # Temperature [K]
P = 101325.  # Pressure    [Pa]

# To have loading rate in nN/s and
# all forces in nN in the input and
# output of Dissociation methods
diss = Dissociation(cogef)


fstep = 0.003
for T in [5, 50, 100, 200, 300, 500]:
    if 0:
        # Manual limits
        fmin, fmax = 0.2, 2
        DPDF, F = diss.probability_density(
            T, P, loading_rate=alpha, force_min=fmin, force_max=fmax,
            force_step=fstep, method='electronic', probability_break_value=0.0)
        p = plt.plot(F / nN, DPDF, '-', label=None)
        color = p[0].get_color()
    else:
        color = None
    # Automatic limits
    fmin, fmax = diss.get_force_limits(T, P, alpha, force_step=fstep,
                                       method='electronic')
    DPDF, F = diss.probability_density(
        T, P, loading_rate=alpha, force_max=fmax, force_min=fmin,
        force_step=fstep, method='electronic')
    p = plt.plot(F / nN, DPDF, 'o-', color=color,
                 label='T={0}K'.format(T))
    plt.axvline(x=cogef.get_maximum_force(method='use_energies') / nN,
                ls='--', color='k')

plt.legend(loc=2)
plt.ylim(0, 110)
plt.xlabel('F [nN]')
plt.ylabel('dp/dF [nN$^{-1}$]')
plt.savefig('1scogef.png')

The peaks define the average rupture force and the standard deviation. These can be obtained using the Dissociation object:

from cogef import COGEF
from cogef import Dissociation
from cogef.units import nN

cogef = COGEF(0, 2)
diss = Dissociation(cogef)

LOADING_RATE = 10  # [nN/s]
P = 101325.  # Pressure    [Pa]

fstep = 0.003
for T in [5, 50, 100, 200, 300, 500]:  # Temperature [K]
    fmin, fmax = diss.get_force_limits(T, P, LOADING_RATE, force_step=fstep,
                                       method='electronic')
    force, error = diss.rupture_force_and_uncertainty(
        T, P, LOADING_RATE, fmax, fmin, fstep, method='electronic')
    print('Rupture force ({0}K): ({1:4.2f} +- {2:4.2f}) nN'.format(
        T, force / nN, error))

Rupture force calculation can be performed automatically, that is the number of needed steps for molecule stretching is obtained. The loop in the next example stops when there are enough images in order to obtain the rupture force:

T = 300
fmax = 0.01
stepsize = 0.1
max_steps = 100

image = Atoms('AuAgAg', positions=((-1, 0, 0), (0, 0, 0), (1, 0, 0)))
image.set_calculator(EMT())
FIRE(image).run(fmax=fmax)

cogef = COGEF([image], 0, 2, optimizer=FIRE, fmax=fmax)
diss = Dissociation(cogef)
for step in range(max_steps + 1):
    try:
        fmin, fmax = diss.get_force_limits(T, P, LOADING_RATE,
                                           force_step=fstep,
                                           method='electronic')
        break
    except ValueError as msg:
        if diss.error not in [1, 2, 3]:
            raise ValueError(msg)
        assert step < max_steps, 'Step limit reached.'
        # Add one image
        cogef.pull(stepsize, 1, initialize, 'cogef.traj')
force, error = diss.rupture_force_and_uncertainty(
    T, P, LOADING_RATE, fmax, fmin, fstep, method='electronic')
# in nN
force = force / nN
error = error / nN
print('Rupture force: ({0:4.3f} +- {1:4.3f}) nN'.format(force, error))

Gibbs energies can be used instead of electronic energies by a vibrational analysis and a calculation based on IdealGasThermo. The difference in rupture force is very small for this example:

from ase.calculators.emt import EMT
from cogef import COGEF
from cogef import Dissociation
from cogef.units import nN

cogef = COGEF(0, 2)

LOADING_RATE = 10  # [nN/s]
T = 300      # Temperature [K]
P = 101325.  # Pressure    [Pa]

def initialize(image, dummy_remove):  # XXX
    image.calc = EMT()
    return image

fstep = 0.003
diss = Dissociation(cogef, initialize, vib_method='frederiksen')
diss.geometry = 'linear'  # Like in class IdealGasThermo

fmin, fmax = diss.get_force_limits(T, P, LOADING_RATE, force_step=fstep,
                                   method='Gibbs')
force, error = diss.rupture_force_and_uncertainty(
    T, P, LOADING_RATE, fmax, fmin, fstep, method='Gibbs')
print('Rupture force (Gibbs): ({0:4.3f} +- {1:4.3f}) nN'.format(force / nN, error))

Most probable force

Rupture force or the most probable force can also be obtained analytically using Bell barrier or the General barrier. We can calculate the force dependent peak values of rupture force, the distribution and the range of most probable forces analytically.

from matplotlib import pyplot as plt

from cogef import COGEF, MPF
from cogef import Dissociation
from cogef.units import nN, _hplanck, _e
import numpy as np

h = _hplanck / _e

cogef = COGEF(0, 2)
energies, distances = cogef.get_energy_curve()
E = list(energies)
last_intact_bond = np.argmax(E)

diss = Dissociation(cogef)
cogef.last_intact_bond_image = last_intact_bond

D0 = diss.electronic_energy_barrier(0)
print('De=', D0, 'eV')

Fmax_cogef = cogef.get_maximum_force(method='use_energies')  # eV/A
print('Fmax = {0:4.2f} eV/A'.format(Fmax_cogef), '=', Fmax_cogef / nN, 'nN')

mpf = MPF(D0, Fmax_cogef)
fstep = 0.0001

# values for calculation
T = 300
anN = 10 # nN/s
print(' ')
print('For T =', T, 'K and alpha =', anN, 'nN/s')

################# Relative forces (f = F/Fmax) ###############################
fmp_bell = mpf.bell(T_K=T, alpha_nNs=anN)

delf = mpf.width_f(T_K=T)
b1 = (fmp_bell + delf / 2)
b2 = (fmp_bell - delf / 2)

fmp_gen = mpf.general(T_K=T, alpha_nNs=anN)
g1 = (fmp_gen + delf / 2)
g2 = (fmp_gen - delf / 2)

fmp_gen_ver = mpf.general_verify(fmp_gen, T_K=T, alpha_nNs=anN)

#################################################
dpdf_bell, f_dpdf_bell = mpf.bell_pd(T_K=T, alpha_nNs=anN, force_step=fstep,
                                     prob_break_value=0.)
idx_bell = np.argmax(dpdf_bell)
fmp_bell_dpdf = f_dpdf_bell[idx_bell]

################################################
dpdf_gen, f_dpdf_gen = mpf.general_pd(T_K=T, alpha_nNs=anN, force_step=fstep,
                                      prob_break_value=0.)
idx_gen = np.argmax(dpdf_gen)
fmp_gen_dpdf = f_dpdf_gen[idx_gen]

print('        Actual most probable force  (F) values in nN          ')
fmp_bell = round(fmp_bell * Fmax_cogef / nN, 4)
b1 = round(b1 * Fmax_cogef / nN, 4)
b2 = round(b2 * Fmax_cogef / nN, 4)
print('Most probable force_Bell is', fmp_bell, 'nN ranging '
                                               'from', b2, 'to', b1, 'nN')
fmp_bell_dpdf = round(fmp_bell_dpdf * Fmax_cogef / nN, 4)
print('From numerical, maximum force_Bell is', fmp_bell_dpdf, 'nN')

print(' ')
fmp_gen = round(fmp_gen * Fmax_cogef / nN, 4)
g1 = round(g1 * Fmax_cogef / nN, 4)
g2 = round(g2 * Fmax_cogef / nN, 4)
print('Most probable force_General:', fmp_gen, 'nN ranging '
                                               'from', g2, 'nN to', g1, 'nN')
fmp_gen_ver = round(fmp_gen_ver * Fmax_cogef / nN, 4)
print('From verification, MPF_General:', fmp_gen_ver, 'nN')

fmp_gen_dpdf = round(fmp_gen_dpdf * Fmax_cogef / nN, 4)
print('From numerical, maximum force_General is', fmp_gen_dpdf, 'nN')

#numerical values in nN
f_dpdf_bell = f_dpdf_bell * Fmax_cogef / nN
f_dpdf_gen = f_dpdf_gen * Fmax_cogef / nN

fig = plt.figure()
ax = fig.add_subplot()
plt.plot(f_dpdf_bell, dpdf_bell, 'o-', color='orange', label='Bell')
plt.plot(f_dpdf_gen, dpdf_gen, 'o-', color='cornflowerblue', label='General')
plt.axvline(x=b1, ls='--', color='orange')
plt.axvline(x=b2, ls='--', color='orange')
plt.axvline(x=g1, ls='--', color='cornflowerblue')
plt.axvline(x=g2, ls='--', color='cornflowerblue')
plt.xlim(0.60, 1.15)
plt.xlabel('Force (nN)')
plt.ylabel('dp/dF (nN$^{-1}$)')
plt.title(r'$T$={} $K$,  '.format(T) + r'$\alpha$={} $nN/s$'.format(anN))
plt.legend()
plt.savefig('mpf.png')
../../_images/mpf.png