Relaxing a Transition State

Greetings Lammpsters,
I am trying to optimize a force field with a large number of parameters against a set of known information obtained by dft. I need a good way to efficiently compare the transition state of a reaction obtained from dft (the known), with the transition state obtained by the force field. The reason it has to be efficient is that I need to be able to tweak parameters and quickly see the changes caused in the location and energy of the TS.

I have three questions in this message:

  1. Is there an efficient way to use neb in the replica package on a single processor? I have found when using user-reaxc that even a small group of atoms (20-30) that I need to specify at least 3gb memory in my PBS script in order to not get segmentation fault. Previously, I have tried to crowd 16 neb images onto one processor and found that the simulations are very slow with only 100 steps coming every minute or so. Is there a more efficient way of using neb on single processor? Why does reaxc need so much memory on my system?

  2. Another method I am considering is to fix the bonds involved in the reaction and relax the rest of the structure of the ts. Previously I have set the forces on these atoms to zero, but I have found that in doing so this neglects the angle and torsion strain in these atoms. I would like to be able to only fix the bonds in the reacting atoms (in order to prevent minimizing back down to reactants or products) without fixing their angles or torsions. It seems that fix shake is not compatible with reaxc as no actual “bonds” are specified in the data file. Is it possible to define a group of atoms with bonds of a specific bond type in order to use shake algorithm? Alternatively, I can use fix restrain with a stiff constant, but I would still like the minimum energy of the system to reflect the minimum energy of the force field and not the restraint+force field. What have others been doing?

  3. Are there any lammps packages for quickly finding a transition state from a guess state by doing something like eigenvector following etc.?

Best,

Josh

hi josh,

here is a suggestion for one part of your question and a comment on another.

[...]

2) Another method I am considering is to fix the bonds involved in the
reaction and relax the rest of the structure of the ts. Previously I have
set the forces on these atoms to zero, but I have found that in doing so
this neglects the angle and torsion strain in these atoms. I would like to
be able to only fix the bonds in the reacting atoms (in order to prevent
minimizing back down to reactants or products) without fixing their angles
or torsions. It seems that fix shake is not compatible with reaxc as no
actual "bonds" are specified in the data file. Is it possible to define a
group of atoms with bonds of a specific bond type in order to use shake
algorithm? Alternatively, I can use fix restrain with a stiff constant, but

yes, it should be possible to use fix shake with reaxc, but it
requires some trickery. you can define bonds, but to avoid having
those bonded pairs excluded from the neighborlist, you would have to
set special_bonds lj/coul to 1.0 1.0 1.0

mind you, fix shake only works for MD, not for minimization and thus
NEB (as a matter of principle). so setting restraints with stiff force
constants is the way to go. a minimizer can handle this much better
than and integrator. but please see the discussion on the restrain
force energy in the fix restrain docs. you still have access to the
individual energy components.

I would still like the minimum energy of the system to reflect the minimum
energy of the force field and not the restraint+force field. What have
others been doing?

since LAMMPS is primarily use as an MD engine, transition states are
(or at least have been) of lesser interest. what people primarily are
looking at are free energy barriers and for that there is a variety of
methods, many of which are available in LAMMPS. if your target is to
generate force field parameters for molecular dynamics runs, than you
may be missing significant contributions by only looking at minimum
and transition state energy structures.

axel.

Joshua,

Allow me to drop some comments in addition to the ones from Axel.

Below is a comment on your 1st Q. Ray may wish to
comment on memory usage within Reax/C.

Steve

Greetings Lammpsters,
I am trying to optimize a force field with a large number of parameters against a set of known information obtained by dft. I need a good way to efficiently compare the transition state of a reaction obtained from dft (the known), with the transition state obtained by the force field. The reason it has to be efficient is that I need to be able to tweak parameters and quickly see the changes caused in the location and energy of the TS.

I have three questions in this message:

  1. Is there an efficient way to use neb in the replica package on a single processor? I have found when using user-reaxc that even a small group of atoms (20-30) that I need to specify at least 3gb memory in my PBS script in order to not get segmentation fault. Previously, I have tried to crowd 16 neb images onto one processor and found that the simulations are very slow with only 100 steps coming every minute or so. Is there a more efficient way of using neb on single processor? Why does reaxc need so much memory on my system?

If you use neb on a single processor and run multiple replicas, then you
are presumably using MPI to create many virtual processes, one
for each replica, all of which run on the same physical core. This
will always run slower than running one replica per physical core.
Is there a reason you don’t want to run NEB in that more typicaly
fashion?

Josh,

It is extremely unlikely that a system of less than 30 atoms would require 3 GB of memory. The benchmark test of a 32,480 atom structure took only less than 1 GB (http://lammps.sandia.gov/bench.html). Something else might be happening in your input.

If you wish to optimize a ReaxFF force field, it might be faster and easier to acquire a copy of Prof. Adri van Duin’s fitting code. It does have the bond constraint feature you wish to use.

Ray