Energy preservation of molecular dynamics

Total energy, geometry optimization, DFT+U, spin....

Moderator: bguster

Locked
tsukagoshi
Posts: 9
Joined: Fri Jun 15, 2012 2:19 am

Energy preservation of molecular dynamics

Post by tsukagoshi » Thu Jul 25, 2013 1:35 pm

When I calculated molecular dynamics of a H2 molecule using the verlet algorithm with ionmov=6 implemented in abinit-7.2.2, I found the sum of the electronic energy and kinetic energy of ions, which were written at the end of each step in a output file, was not preserved.
However, when I ignore the kinetic energy at a first step, the total energy was successfully preserved.
I think that this behavior comes from the difference between the timing when the electronic energy is calculated and when the kinetic energy is calculated.
From the source code of mover.F90 and pred_verlet.F90, it seems that the ionic velocity has not been calculated yet after the second step when the kinetic energy is calculated, although the electronic energy for the current ionic positions has been already calculated using scfcv routine.
As far as I know, this inconsistency does not occurred for the abinit-5.8.4, in which the kinetic energy is calculated after the ionic velocity is calculated in moldyn.F90. Would you reconsider about the source code ?

The following is my input file.
h2.in:
#molecular dynamics of a H2 molecule
#SCF
iscf 7
toldfe 1.d-8

#k-samplig
kptopt 0
nkpt 1

#electronic part
ixc 1
ecut 10.0
nsppol 1
nband 1

#atomic part
natom 2
ntypat 1
znucl 1
typat 2*1
ionmov 6
ntime 100
dtion 5

acell 3*10.0
rprim 1.0 0.0 0.0
0.0 1.0 0.0
0.0 0.0 1.0
xred
0.4 0.5 0.5
0.6 0.5 0.5


I can extract the electronic energy and the kinetic energy as follows.

$ grep "Total energy (etotal)" h2.out | cut -d"=" -f2 > etotal.dat
$ grep "Kinetic energy of ions (ekin)" h2.out | cut -d"=" -f2 > ekin.dat
$ paste etotal.dat ekin.dat > energy.dat

I can see the time evolution of the sum of the electronic energy and the kinetic energy using gnuplot as follows.

gnuplot> plot 'energy.dat' u 1 w lp, 'energy.dat' u ($1+$2) w lp

Locked