[ase-users] Dimer
Michael Joseph Waters
michael.j.waters at northwestern.edu
Fri Jan 24 22:44:21 CET 2020
Hi Zainab,
I am assuming that these calculations are ab initio calculations. Before we continue, I need to know are these spin polarized or magnetic? Which software are you using?
If you just want the exact barrier, consider doing a climbing image NEB calculation once the normal NEB has converged: https://wiki.fysik.dtu.dk/ase/ase/neb.html#climbing-image
I think you can use this snippet of my code to get what you want:
#### Set up the dimer method
d_control = DimerControl(initial_eigenmode_method = 'displacement',
displacement_method = 'vector', logfile = None,
mask = dimer_mask,
maximum_translation = 0.1,
dimer_separation = 0.1)
Most dimer searches start with a displacement from equilibrium. If you already know the rough saddle point, set the displacements to zero and feed in your rough saddle point here:
d_atoms = MinModeAtoms(atoms, d_control)
### Configure dimer atoms
# Displace the atoms
displacement_vector = [[0.0]*3]*(len(atoms))
displacement_vector[28] = [-0.2, 0.1, 0.1]
displacement_vector = fdstep *normalize(displacement_vector)
# set dimer direction, optional but helps reconverge
d_atoms.displace(displacement_vector = displacement_vector)
d_atoms.initialize_eigenmodes(method = 'displacement') # set by displacements
Since you already roughly know your saddle, you can orient your dimer direction vector to point across the saddle based on the neighboring NEB images. Set that direction vector as your eigenmode. Keep in mind that the eigenmodes must be initialized before overwriting them.
# Setting the eigenmode explicitly, must be normalized for good results,
# using a negative sign will find the saddle backwards
d_atoms.set_eigenmode(normalize(displacement_vector))
### Converge to a saddle point
dim_rlx = MinModeTranslate(d_atoms, trajectory = 'dimer_method.traj', logfile = logfile)
dim_rlx.run(fmax = fmax)
saddle = d_atoms.get_atoms()
# get dimer direction or eigenmode, returns normalized
dimer_direction = d_atoms.get_eigenmode()
## confirm saddle
saddle.set_calculator(EMT())
saddle_energy = saddle.get_potential_energy()
print('forces at saddle:', atoms.get_forces())
print('force at saddle:', np.linalg.norm(atoms.get_forces())) # a good check for convergence
Best,
-Mike
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk <ase-users-bounces at listserv.fysik.dtu.dk> on behalf of Alaithan, Zainab via ase-users <ase-users at listserv.fysik.dtu.dk>
Sent: Friday, January 24, 2020 7:09 AM
To: Michael Joseph Waters via ase-users <ase-users at listserv.fysik.dtu.dk>
Subject: Re: [ase-users] Dimer
Thank you some much for sharing your example. I think I am using the opposite approach. because my NEB calculation is converged and I was looking to refine the saddle point using the dimer method. My main problem is that I don't fully understand the dimer object d_atoms. The trajectory contains only 1 image per step, which I thought to contain two images. I am now assuming that the trajectory contains only the middle point forces and positions. I am running expensive calculation and I am limited by wall time. That's why I asked how to restart from previous dimer trajectory.
Should I keep an eigenmode log of the dimer and provide the initial eigenmodes ? I am thinking that it is best to start from exactly the same place the dimer calculation stopped but I am not sure how I can achieve this ? or there is not other option to restart?
No body is answering except for you. I am very grateful.
Regards,
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk <ase-users-bounces at listserv.fysik.dtu.dk> on behalf of Michael Joseph Waters via ase-users <ase-users at listserv.fysik.dtu.dk>
Sent: 23 January 2020 18:21
To: Alaithan, Zainab via ase-users <ase-users at listserv.fysik.dtu.dk>
Subject: Re: [ase-users] Dimer
Hi Zainab,
I am sorry that you are also experiencing a block on .py files. I like my example because it uses the saddle point as preconditioning to reduced the computational cost of the NEB by ~30% Here is the file that I sent earlier:
#!/usr/bin/env python3
from ase.visualize import view
from ase.build import fcc100, add_adsorbate
from ase.constraints import FixAtoms
from ase.calculators.emt import EMT
from ase.optimize import BFGS
from ase.dimer import DimerControl, MinModeAtoms, MinModeTranslate
import numpy as np
def normalize(a):
return a/np.linalg.norm(a)
logfile = '-'
fdstep = 0.2
fmax = 0.01
# Set up a small "slab" with an adatoms
atoms = fcc100('Pt', size = (3, 3, 2), vacuum = 10.0)
add_adsorbate(atoms, 'Cu', 1.611, 'hollow')
# Freeze the "slab"
mask = [atom.z < 11.0 for atom in atoms]
atoms.set_constraint(FixAtoms(mask = mask))
#dimer_mask = np.ones(len(atoms)) # doesn't work???
dimer_mask = len(atoms)*[1]
for flag,i in enumerate( mask):
if flag:
dimer_mask[i]=0
# Calculate using EMT
atoms.set_calculator(EMT())
#### relax initial structure
initial_energy = atoms.get_potential_energy()
dyn = BFGS(atoms, trajectory = 'bfgs.traj',logfile = logfile)
dyn.run(fmax=fmax)
relaxed_energy = atoms.get_potential_energy()
# we could compare this to minimum1 and minimum2 later
relaxed_structure = atoms.copy()
#### Set up the dimer method
d_control = DimerControl(initial_eigenmode_method = 'displacement',
displacement_method = 'vector', logfile = None,
mask = dimer_mask,
maximum_translation = 0.1,
dimer_separation = 0.1)
### Configure dimer atoms
d_atoms = MinModeAtoms(atoms, d_control)
# Displace the atoms
displacement_vector = [[0.0]*3]*(len(atoms))
displacement_vector[-1] = [0.1, 0.02, 0.1]
displacement_vector = fdstep *normalize(displacement_vector)
# set dimer direction, optional but helps reconverge
d_atoms.displace(displacement_vector = displacement_vector)
d_atoms.initialize_eigenmodes(method = 'displacement') # set by displacements
# Setting the eigenmode explicitly, must be normalized for good results,
# using a negative sign will find the saddle backwards
d_atoms.set_eigenmode(normalize(displacement_vector))
### Converge to a saddle point
dim_rlx = MinModeTranslate(d_atoms, trajectory = 'dimer_method.traj', logfile = logfile)
dim_rlx.run(fmax = fmax)
saddle = d_atoms.get_atoms()
# get dimer direction or eigenmode, returns normalized
dimer_direction = d_atoms.get_eigenmode()
## confirm saddle
saddle.set_calculator(EMT())
print('forces at saddle:', atoms.get_forces())
print('force at saddle:', np.linalg.norm(atoms.get_forces()))
##### relax forward through the saddle
displacement_vector = fdstep * normalize(dimer_direction)
atoms = saddle.copy()
atoms.set_positions(atoms.get_positions()+displacement_vector)
atoms.set_calculator(EMT())
dyn = BFGS(atoms, trajectory = 'bfgs_min1.traj',logfile = logfile)
dyn.run(fmax=fmax)
minimum1 = atoms.copy()
#### relax backwards through the saddle
displacement_vector = -fdstep * normalize(dimer_direction)
atoms = saddle.copy()
atoms.set_positions(atoms.get_positions()+displacement_vector)
atoms.set_calculator(EMT())
dyn = BFGS(atoms, trajectory = 'bfgs_min2.traj',logfile = logfile)
dyn.run(fmax=fmax)
minimum2 = atoms.copy()
###### neb
use_saddle = True
neb_images = 7
from ase.neb import NEB, SingleCalculatorNEB
from ase.optimize import MDMin
# Read initial and final states:
if use_saddle:
images = [minimum2.copy()]
images += [saddle.copy()]
images += [minimum1.copy()]
neb = NEB(images)
scneb=SingleCalculatorNEB(images)
# Linear Interpolation between minimum1 -> saddle
# and saddle -> minimum2 points
# why isn't this function also attached to NEB?
scneb.refine(steps = (neb_images-1)//2 )
neb = NEB(scneb.images)
else:
images = [minimum2.copy()]
images += [minimum2.copy() for i in range(neb_images)]
images += [minimum1.copy()]
neb = NEB(images)
# Linear Interpolation between minimum1 and minimum2 points:
neb.interpolate()
# Set calculators:
for image in neb.images:
image.set_calculator(EMT())
# Optimize:
optimizer = BFGS(neb, trajectory='neb.traj',logfile=logfile)
optimizer.run(fmax=fmax)
#### for viewing with enegies and forces
for image in neb.images:
image.set_calculator(EMT())
image.get_potential_energy()
view(neb.images)
###################
from ase.neb import NEBTools
from matplotlib import pyplot as plt
nebtools = NEBTools(neb.images)
# Get the calculated barrier and the energy change of the reaction.
Ef, dE = nebtools.get_barrier()
# Get the barrier without any interpolation between highest images.
Ef, dE = nebtools.get_barrier(fit=False)
# Get the actual maximum force at this point in the simulation.
max_force = nebtools.get_fmax()
# Create a figure with custom parameters.
fig2 = plt.figure(figsize=(5.5, 4.0))
ax = fig2.add_axes((0.15, 0.15, 0.8, 0.75))
nebtools.plot_band(ax)
fig2.savefig('diffusion-barrier.png')
plt.show()
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk <ase-users-bounces at listserv.fysik.dtu.dk> on behalf of Alaithan, Zainab via ase-users <ase-users at listserv.fysik.dtu.dk>
Sent: Thursday, January 23, 2020 10:33 AM
To: Michael Joseph Waters via ase-users <ase-users at listserv.fysik.dtu.dk>
Subject: Re: [ase-users] Dimer
Thanks a lot mike. it seems that I am blocked from downloading this file. If you don't mind can you try sending it with a text extension. Perhaps I can download it.
Regards,
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk <ase-users-bounces at listserv.fysik.dtu.dk> on behalf of Michael Joseph Waters via ase-users <ase-users at listserv.fysik.dtu.dk>
Sent: 23 January 2020 16:16
To: Offermans Willem via ase-users <ase-users at listserv.fysik.dtu.dk>
Subject: Re: [ase-users] Dimer
This email from ase-users at listserv.fysik.dtu.dk originates from outside Imperial. Do not click on links and attachments unless you recognise the sender. If you trust the sender, add them to your safe senders list<https://urldefense.proofpoint.com/v2/url?u=https-3A__spam.ic.ac.uk_SpamConsole_Senders.aspx&d=DwMFAw&c=yHlS04HhBraes5BQ9ueu5zKhE7rtNXt_d012z2PA6ws&r=IWIXV02QNefgPQmc-Ovk6evJ6KEOZjyl5g4atrVIlcJ4WAvRgH8dNpuAIeIrVHn4&m=nlpqyD29QIFfqC_jD8jnPFUfyBAYoydVGiR9JD5nVd4&s=_ShoFAwFAQMK2sNKxkK4cB6EqpB5altA-Kf8ZY4RI4s&e=> to disable email stamping for this address.
I agree,
The MinModeAtoms/DimerAtoms objects are not well explained. The things I think are missing is setting/getting the dimer location and direction/eigenmode. I've been building my personal example on finding reaction pathways which shows how to do these things. It takes about 20 - 30 seconds to run. Feel free to use it and acknowledge me when you cure cancer!
Cheers,
-Mike
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk <ase-users-bounces at listserv.fysik.dtu.dk> on behalf of Offermans Willem via ase-users <ase-users at listserv.fysik.dtu.dk>
Sent: Thursday, January 23, 2020 8:49 AM
To: ase users <ase-users at listserv.fysik.dtu.dk>
Subject: Re: [ase-users] Dimer
Dear Zainab and ASE friends,
To my opinion, Zainab is missing basic documentation on how to use the dimer method with ASE.
How to restart a dimer calculation?
How to monitor convergence?
How to identify the possible transitions state in the trajectory file?
How to retrieve the energy of the possible transition state?
How to retrieve the remaining forces of the possible transition state?
Met vriendelijke groeten,
Mit freundlichen Grüßen,
With kind regards,
Willem Offermans
Researcher Electrocatalysis SCT
VITO NV | Boeretang 200 | 2400 Mol
Phone:+32(0)14335263 Mobile:+32(0)492182073
Willem.Offermans at Vito.be<mailto:Willem.Offermans at Vito.be>
[cid:982BA063-B96A-4A1B-89AB-5A01CA9FC70D at vito.local]
On 23 Jan 2020, at 15:29, Michael Joseph Waters via ase-users <ase-users at listserv.fysik.dtu.dk<mailto:ase-users at listserv.fysik.dtu.dk>> wrote:
Hi Zainab,
You may want to try changing your optimizer. The documentation lists: BFGS, BFGSLineSearch, LBFGS, LBFGSLineSearch, GPMin, MDMin and FIRE as options. Many of the optimizers have options for maxstep/maxmove which limit large position jumps. Some optimizers have a damping option which can help.
Are you using molecular mechanics or an ab initio method like DFT for your dimer calculations?
If you are using molecular mechanics, try extending your cut-off distances.
If you are using an ab initio method, increasing your spatial sampling/plane wave cutoff and using tighter energy tolerances can help. If you are working with a spin-polarized ab initio method, things can get trickier due to the coexistence of multiple spin PES.
Best,
-Mike
________________________________
From: ase-users-bounces at listserv.fysik.dtu.dk<mailto:ase-users-bounces at listserv.fysik.dtu.dk> <ase-users-bounces at listserv.fysik.dtu.dk<mailto:ase-users-bounces at listserv.fysik.dtu.dk>> on behalf of Alaithan, Zainab via ase-users <ase-users at listserv.fysik.dtu.dk<mailto:ase-users at listserv.fysik.dtu.dk>>
Sent: Thursday, January 23, 2020 7:06 AM
To: Ase-users <ase-users at listserv.fysik.dtu.dk<mailto:ase-users at listserv.fysik.dtu.dk>>
Subject: [ase-users] Dimer
Dear Ase-user,
Can someone please share some tips about how to restart the dimer calculation from the dimer trajectory? and what is the appropriate displacement vector in this case? especially if the dimer calculation was converging and I just want to tighten the convergence criteria.
I will appreciate any help.
Regards,
_______________________________________________
ase-users mailing list
ase-users at listserv.fysik.dtu.dk<mailto:ase-users at listserv.fysik.dtu.dk>
https://listserv.fysik.dtu.dk/mailman/listinfo/ase-users<https://urldefense.proofpoint.com/v2/url?u=https-3A__listserv.fysik.dtu.dk_mailman_listinfo_ase-2Dusers&d=DwMGaQ&c=yHlS04HhBraes5BQ9ueu5zKhE7rtNXt_d012z2PA6ws&r=IWIXV02QNefgPQmc-Ovk6evJ6KEOZjyl5g4atrVIlcJ4WAvRgH8dNpuAIeIrVHn4&m=oHY7WPmo7YhF6NflQXTasP8vcx6DTay_8bDf9c9-jvo&s=8GefjROfBDg0fmQvHLziNv8tpfjIUCZH9yTHrP2bJHw&e=>
Indien u VITO Mol bezoekt, hou aub er dan rekening mee dat de hoofdingang voortaan enkel bereikbaar is vanuit de richting Dessel-Retie, niet vanuit richting Mol, zie vito.be/route.<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.vito.be_route&d=DwMGaQ&c=yHlS04HhBraes5BQ9ueu5zKhE7rtNXt_d012z2PA6ws&r=IWIXV02QNefgPQmc-Ovk6evJ6KEOZjyl5g4atrVIlcJ4WAvRgH8dNpuAIeIrVHn4&m=oHY7WPmo7YhF6NflQXTasP8vcx6DTay_8bDf9c9-jvo&s=Sb5zVibfDy5DBVtSIIXLd0KLlA-JK8lt6WlmL7FDHW0&e=>
If you plan to visit VITO at Mol, then please note that the main entrance can only be reached coming from Dessel-Retie and no longer coming from Mol, see vito.be/en/contact/locations.<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.vito.be_en_contact_locations&d=DwMGaQ&c=yHlS04HhBraes5BQ9ueu5zKhE7rtNXt_d012z2PA6ws&r=IWIXV02QNefgPQmc-Ovk6evJ6KEOZjyl5g4atrVIlcJ4WAvRgH8dNpuAIeIrVHn4&m=oHY7WPmo7YhF6NflQXTasP8vcx6DTay_8bDf9c9-jvo&s=LuSseMTxtADpSbPgDQkHn5_4-ce0xCSCAPRsRzhfqa4&e=>
VITO Disclaimer: http://www.vito.be/e-maildisclaimer<https://urldefense.proofpoint.com/v2/url?u=http-3A__www.vito.be_e-2Dmaildisclaimer&d=DwMGaQ&c=yHlS04HhBraes5BQ9ueu5zKhE7rtNXt_d012z2PA6ws&r=IWIXV02QNefgPQmc-Ovk6evJ6KEOZjyl5g4atrVIlcJ4WAvRgH8dNpuAIeIrVHn4&m=oHY7WPmo7YhF6NflQXTasP8vcx6DTay_8bDf9c9-jvo&s=ZmcLEylWWBmaVda1_1n-uQrz441t5gu7vwPod0TGfMY&e=>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20200124/742723dd/attachment-0001.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: vito.jpg
Type: image/jpeg
Size: 15232 bytes
Desc: vito.jpg
URL: <http://listserv.fysik.dtu.dk/pipermail/ase-users/attachments/20200124/742723dd/attachment-0001.jpg>
More information about the ase-users
mailing list