AuAg₂ 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 want to perform this now for a linear Au-Ag-Ag chain using the ASE built in EMT potential. Optimization/relaxation is done with the FIRE (Fast Inertial Relaxation Engine) algorithm. The convergence criterion is that the force on all individual atoms should be less than f_max

We first relax the chain:

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.calc = EMT()
    FIRE(image).run(fmax=fmax)
    image.write(fname)

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

from ase import io
from ase.calculators.emt import EMT

from cogef import COGEF1D

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

fmax = 0.01
cogef = COGEF1D(0, 2, fmax=fmax)

if not len(cogef.images):  # calculation was done already
    cogef.images = [image]

    stepsize = 0.1
    steps = 30
    cogef.move(stepsize, steps)

which creates the directory cogef_0_2 and writes the file cogef1d_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 cogef import COGEF1D
from cogef.units import nN

cogef = COGEF1D(0, 2)  # this reads cogef1d_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') / nN) +
    ' {0:4.2f} nN (from forces)'.format(
        cogef.get_maximum_force(method='use_forces') / nN))

Dissociation

The figure shows the force-dependent probability distributions dp/dF for different temperatures. The mean rupture force and standard deviation are defined by the peaks.

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 ** link for the dissociation file probability distributions \(dp/dF\) as shown in the figure above.

The figure was created with:

import pylab as plt
from cogef.generalized import COGEF1D
from cogef import Dissociation
from cogef.dissociation.barrier import ElectronicBarrier

from cogef.units import nN

cogef = COGEF1D(0, 2)

# The force is assumed to increase uniformly by this loading rate.
alpha_nNs = 10  # [nN/s]
alpha = alpha_nNs * nN

diss = Dissociation(ElectronicBarrier(cogef))

fstep = 0.001
for T in [5, 50, 100, 200, 300, 500]:  # Temperature [K]
    # Automatic limits
    fmin, fmax = diss.get_force_limits(T, alpha, force_step=fstep,
                                       method='electronic')
    dpdf, F = diss.probability_density(
        T, loading_rate=alpha, force_max=fmax, force_min=fmin,
        force_step=fstep)
    p = plt.plot(F / nN, dpdf * nN, 'o-',
                 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.plot()

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:

 ### Rupture force calculation
 
  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.

# creates: mpf.png

from matplotlib import pyplot as plt

from cogef import COGEF1D, MPF
from cogef import Dissociation
from cogef.dissociation.diss_old import Dissociation as Dissociation_old
from cogef.units import nN#, _hplanck, _e
import numpy as np

#h = _hplanck / _e

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

diss_old = Dissociation_old(cogef)
cogef.last_intact_bond_image = last_intact_bond

D0 = diss_old.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')