Implementation of a numerical Hessian via a new fix

Dear users,

I want to implement the computation of a finite difference based Hessian in LAMMPS. My Idea is to write a new fix. The fix should use the methods initial_integrate and final_integrate. In the initial_integrate method the an atoms coordinate is displaced, and in the final_integrate method the actual finite differencing should be done. By storing the current atom and coordinate (x,y or z) as private variables of the fix class, it is always known which is the next atom to process.

Since this is the first time I write an extension to LAMMPS I want to ask two questions:

  • Does this approach make sense, or is there a more clever/easier way of doing this?
  • How do I have to setup the actual Hessian matrix (2d array) so that I can access it after the computation of the Hessian? The ideal case would be to extract it directly with the extract_fix method of the python interface.
    Best

Johannes Dürholt

Dear users,

I want to implement the computation of a finite difference based Hessian
in LAMMPS. My Idea is to write a new fix. The fix should use the methods
initial_integrate and final_integrate. In the initial_integrate method the
an atoms coordinate is displaced, and in the final_integrate method the
actual finite differencing should be done. By storing the current atom and
coordinate (x,y or z) as private variables of the fix class, it is always
known which is the next atom to process.

Since this is the first time I write an extension to LAMMPS I want to ask
two questions:

   - Does this approach make sense, or is there a more clever/easier way
   of doing this?
   - How do I have to setup the actual Hessian matrix (2d array) so that
   I can access it after the computation of the Hessian? The ideal case would
   be to extract it directly with the extract_fix method of the python
   interface.

​if you are already planning to do the processing from the python

interface, why not do the bulk of the steps in python and only ​run LAMMPS
in between with "run 0" to compute the forces?

axel.