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:
- Run a NVT simulation and wait for the equilibrium
- Calculate per atom potential energy and call it U1 (using compute all per/atom)
- Change the simulation box the way I want (e.g. change_box all z scale 1.0001 y volume x volume remap units box)
- Calculate, again, the per atom potential and name it U2
- Switch back to the original box size. (change_box all z scale 0.9999 y volume x volume remap units box)
- Use fix ave/spatial to sum up per atom quantities say each 100 time steps and print the results to a txt file
- 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
Thanks for your quick reply.
I have to do it numerically for two reasons:
-
I have an elongated simulation box with Lx=a, Ly=a and Lz=3a where there is a stable liquid slab in the middle and 2 vapor phases on each side. Using built-in functions in Lammps will cause some problems since I have to be worry about tail corrections. I can not simply use “pair_modify tail yes” because that is just good for homogenous systems. When I calculate that numerically, the tail correction will be canceled out in a formula like (U2-U1)/dV
-
My final goal is to calculate surface tension and I don’t think there is a built-in command in Lammps to do that. Reproducing the Pz was just and example and also good practice before I tackle surface tension. To get surface tension, I should calculate dU/dA where A is the cross-section area of the simulation box while I keep the volume constant.
Have you got a suggestion for my problem yet?