# Derivatives of total potential energy

Dear sir,
I need to calculate the first derivative of total/per atom potential energy in respect to size of the simulation box in Z direction. To do that I need to:

1. Run a NVT simulation and wait for the equilibrium
2. Calculate per atom potential energy and call it U1 (using compute all per/atom)
3. Change the simulation box the way I want (e.g. change_box all z scale 1.0001 y volume x volume remap units box)
4. Calculate, again, the per atom potential and name it U2
5. Switch back to the original box size. (change_box all z scale 0.9999 y volume x volume remap units box)
6. Use fix ave/spatial to sum up per atom quantities say each 100 time steps and print the results to a txt file
7. Run for 2 ns, for example.

Calculations of this type would give normal component of pressure and/or surface tension using perturbation theory.

The problem is when I do this, the simulation box changes twice before the first time step and never changes again. Consequently, U1 and U2 are the same.
I know I may calculate (U2-U1)/dV by using “run 0” command twice but I need to monitor this derivative during the course of simulation not just once. What would be the most convenient way to do this. Run every command maybe? Do I have to write my own compute or fix commands?

Sorry for the long post and thank you in advance for your support.

Dear sir,
I need to calculate the first derivative of total/per atom potential energy
in respect to size of the simulation box in Z direction. To do that I need
to:
1. Run a NVT simulation and wait for the equilibrium
2. Calculate per atom potential energy and call it U1 (using compute all
per/atom)
3. Change the simulation box the way I want (e.g. change_box all z scale
1.0001 y volume x volume remap units box)
4. Calculate, again, the per atom potential and name it U2
5. Switch back to the original box size. (change_box all z scale 0.9999 y
volume x volume remap units box)
6. Use fix ave/spatial to sum up per atom quantities say each 100 time steps
and print the results to a txt file
7. Run for 2 ns, for example.

Calculations of this type would give normal component of pressure and/or
surface tension using perturbation theory.

The problem is when I do this, the simulation box changes twice before the
first time step and never changes again. Consequently, U1 and U2 are the
same.
I know I may calculate (U2-U1)/dV by using "run 0" command twice but I need
to monitor this derivative during the course of simulation not just once.

well, then your input is not correct. LAMMPS commands
will always be executed immediately.

What would be the most convenient way to do this. Run every command maybe?

depends on what you call convenient.

Do I have to write my own compute or fix commands?

depends on how often you need to do this procedure.

axel.

It seems like the z-component of the virial tensor is (proportional
to) exactly what I think you are trying to calculate, except that it
looks like you are trying to approximate it numerically by rescaling
the coordinates and calculating the derivative numerically. There
should be no reason to do this.
Why not just use LAMMP's built-in command to calculate the pressure
tensor and take the z,z component?
http://lammps.sandia.gov/doc/compute_pressure.html
If you only want the contribution due to the potential energy (leaving
out kinetic energy) then add the "virial" keyword to the end of this
command.

On a related note, to compute only the z components of the kinetic
energy, take a look at:
http://lammps.sandia.gov/doc/compute_temp_partial.html
(Not sure if this is relevant or helpful.)

Cheers

Andrew