AuAg2 chain¶

We want to simulate external forces by the appication 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:
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.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¶ 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 propability distributions $$dp/dF$$ as shown in the figure above.

The figure was created with:

# creates: 1scogef.png
import pylab as plt
from ase.units import J, m

from cogef import COGEF
from cogef import Dissociation

cogef = COGEF(0, 2)

T = 300      # Temperature [K]
P = 101325.  # Pressure    [Pa]

# all forces in nN in the input and
# output of Dissociation methods
diss = Dissociation(cogef, force_unit='nN')

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

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(
p = plt.plot(F, DPDF, '-', label=None)
color = p.get_color()
else:
color = None
# Automatic limits
method='electronic')
DPDF, F = diss.probability_density(
p = plt.plot(F, DPDF, 'o-', color=color,
label='T={0}K'.format(T))
plt.axvline(x=cogef.get_maximum_force(method='use_energies') * f,
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

cogef = COGEF(0, 2)
diss = Dissociation(cogef, force_unit='nN')

P = 101325.  # Pressure    [Pa]

fstep = 0.003
for T in [5, 50, 100, 200, 300, 500]:  # Temperature [K]
method='electronic')
force, error = diss.rupture_force_and_uncertainty(
print('Rupture force ({0}K): ({1:4.2f} +- {2:4.2f}) nN'.format(
T, force, 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, force_unit='nN')
for step in range(max_steps + 1):
try:
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.'
cogef.pull(stepsize, 1, initialize, 'cogef.traj')
force, error = diss.rupture_force_and_uncertainty(
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

cogef = COGEF(0, 2)

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',
force_unit='nN')
diss.geometry = 'linear'  # Like in class IdealGasThermo