scnc
/
siesta
mirror of https://gitlab.com/siesta-project/siesta.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

maint: added CODATA units for vmb, dynamics and fire code 

I still need to figure out exactly how amu is calculated.

Signed-off-by: Nick Papior <nickpapior@gmail.com>
This commit is contained in:
Nick Papior 2020-10-08 12:42:48 +02:00
parent 21d7ea9dc5
commit 9ab8543b9b
4 changed files with 152 additions and 20 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

145
Src/dynamics.f → Src/dynamics.F
View File

@ -16,6 +16,7 @@
use m_mpi_utils, only : broadcast
use files, only : slabel
use alloc, only : re_alloc, de_alloc
use units, only : Kelvin, eV, Ang, amu
implicit none
@ -161,17 +162,34 @@ C Define constants and conversion factors ......................................
twodt = dt*2.0d0
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
C convert target temperature into target kinetic energy
C Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in eV if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * 8.617d-5 * tt
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038
fovermp = 0.009579038d0
#else
! convert target temperature into target kinetic energy
! Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in eV if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * Kelvin * tt / eV
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert target temperature into target kinetic energy
C Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in Ry if Temp is in Kelvin)
tekin = eV * 0.5d0 * (3.d0 * natoms - ct) * 8.617d-5 * tt
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038 * Ang**2 / eV
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! convert target temperature into target kinetic energy
! Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in Ry if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * Kelvin * tt
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C calculate cell volume at current time
@ -541,9 +559,17 @@ C Potential energy of Parrinello-Rahman variables
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5/eV
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin
#endif
endif
C .....................
@ -676,12 +702,24 @@ C Define constants and conversion factors ......................................
twodt = dt*2.0d0
if (iunit .eq. 1) then
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038d0
#else
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038 * Ang**2 / eV
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C calculate cell volume at current time
vol = volcel( h )
C ........................
@ -1063,9 +1101,17 @@ C Potential energy of Parrinello-Rahman variables
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5/eV
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin
#endif
endif
C .....................
@ -1344,9 +1390,17 @@ C Potential energy of Nose variable (in eV)
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = kin / (0.5d0 * (3.d0 * natoms - ct) * 8.617d-5)
#else
temp = kin / (0.5d0 * (3.d0 * natoms - ct) * Kelvin / eV)
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = kin / (0.5d0 * (3.d0 * natoms - ct) * 8.617d-5 * eV)
#else
temp = kin / (0.5d0 * (3.d0 * natoms - ct) * Kelvin)
#endif
endif
if (Node .eq. 0) then
@ -1482,17 +1536,34 @@ C Define constants and conversion factors ......................................
have_p_target = (ianneal .eq. 2 .or. ianneal .eq. 3)
if (iunit .eq. 1) then
C convert target ionic temperature into target kinetic energy
#ifdef SIESTA__UNITS_ORIGINAL
C convert target temperature into target kinetic energy
C Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in eV if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * 8.617d-5 * tt
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038
fovermp = 0.009579038d0
#else
! convert target temperature into target kinetic energy
! Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in eV if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * Kelvin * tt / eV
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert target temperature into target kinetic energy
C Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in Ry if Temp is in Kelvin)
tekin = eV * 0.5d0 * (3.d0 * natoms - ct) * 8.617d-5 * tt
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038 * Ang**2 / eV
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! convert target temperature into target kinetic energy
! Ekin=1/2*(3N-3)*kB*Temp (yields Ekin in Ry if Temp is in Kelvin)
tekin = 0.5d0 * (3.d0 * natoms - ct) * Kelvin * tt
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C calculate cell volume at current time
@ -1716,9 +1787,17 @@ C Kinetic energy of atoms
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/8.617d-5/eV
#else
temp = 3.0d0*vol*pgas/(3.0d0*natoms-ct)/Kelvin
#endif
endif
C .....................
@ -1823,11 +1902,22 @@ C Define constants and conversion factors .....................................
twodt = dt*2.0d0
if (iunit .eq. 1) then
C convert F/m in (eV/Amstrong)/amu to Amstrong/fs**2
fovermp = 0.009579038
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038d0
#else
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038 * Ang**2 / eV
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C ........................
@ -1898,9 +1988,17 @@ C Kinetic energy of atoms
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = 2.0d0*kin/(3.0d0*natoms-ct)/8.617d-5
#else
temp = 2.0d0*kin/(3.0d0*natoms-ct)/Kelvin * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = 2.0d0*kin/(3.0d0*natoms-ct)/8.617d-5/eV
#else
temp = 2.0d0*kin/(3.0d0*natoms-ct)/Kelvin
#endif
endif
C .....................
@ -2009,11 +2107,22 @@ C Define constants and conversion factors .....................................
dtby2 = dt/2.0d0
if (iunit .eq. 1) then
C convert F/m in (eV/Amstrong)/amu to Amstrong/fs**2
fovermp = 0.009579038
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038d0
#else
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038 * Ang**2 / eV
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C ........................
@ -2223,9 +2332,17 @@ C Kinetic energy of atoms
C Instantaneous temperature (Kelvin)
if (iunit .eq. 1) then
#ifdef SIESTA__UNITS_ORIGINAL
temp = 2.0d0*kin/(3.0d0*natoms-ct)/8.617d-5
#else
temp = 2.0d0*kin/(3.0d0*natoms-ct)/Kelvin * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
temp = 2.0d0*kin/(3.0d0*natoms-ct)/8.617d-5/eV
#else
temp = 2.0d0*kin/(3.0d0*natoms-ct)/Kelvin
#endif
endif
!

6
Src/fire_optim.F
View File

@ -12,7 +12,7 @@
use precision, only : dp
use fdf, only: fdf_block
use m_mpi_utils, only: broadcast
use units, only: kBar, eV, Ang
use units, only: kBar, eV, Ang, amu
use alloc, only: re_alloc, de_alloc
use fdf, only: fdf_get
@ -95,7 +95,11 @@ c Internal variables and arrays
$ fire_alpha_factor, fire_dtmax
integer :: fire_nmin
#ifdef SIESTA__UNITS_ORIGINAL
real(dp), parameter :: fovermp = 0.009579038 * Ang**2 / eV
#else
real(dp), parameter :: fovermp = 1._dp / amu
#endif
real(dp), parameter :: dstrain_max = 0.1_dp
real(dp) :: volcel

6
Src/units.F90
View File

@ -68,7 +68,6 @@ module units
real(dp), parameter, public :: Ang = 1._dp / 0.529177_dp
real(dp), parameter, public :: eV = 1._dp / 13.60580_dp
real(dp), parameter, public :: kBar = 1._dp / 1.47108e5_dp
real(dp), parameter, public :: GPa = kBar * 10
real(dp), parameter, public :: Kelvin = eV / 11604.45_dp
real(dp), parameter, public :: Debye = 0.393430_dp
real(dp), parameter, public :: amu = 2.133107_dp
@ -77,15 +76,16 @@ module units
real(dp), parameter, public :: Ang = 1.88972612462577017_dp
real(dp), parameter, public :: eV = 7.34986443513115789e-2_dp
real(dp), parameter, public :: kBar = 6.79786184348648780e-6_dp
real(dp), parameter, public :: GPa = kBar * 10._dp
real(dp), parameter, public :: Kelvin = 6.33362312691136091e-6_dp
real(dp), parameter, public :: Debye = 3.93430269519899511e-1_dp
real(dp), parameter, public :: amu = 2.133107_dp
real(dp), parameter, public :: Ryd_time = 1._dp/0.04837769_dp
real(dp), parameter, public :: Ryd_time = 1._dp/0.048377686531714_dp
#endif
real(dp), parameter, public :: GPa = kBar * 10
! pi to 50 digits
real(dp), parameter, public :: pi = 3.14159265358979323846264338327950288419716939937510_dp
real(dp), parameter, public :: pi2 = pi * 2._dp
real(dp), parameter, public :: deg = pi / 180.0_dp

15
Src/vmb.F
View File

@ -239,7 +239,7 @@ C real*8 tempe : Instantaneous system temperature
C *****************************************************************************
use precision, only : dp
use units, only: eV, Ang
use units, only: eV, Ang, amu
use sys, only : die
implicit none
@ -266,11 +266,22 @@ C ........................
C Define constants and conversion factors .....................................
if (iunit .eq. 1) then
C convert F/m in (eV/Amstrong)/amu to Amstrong/fs**2
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (eV/Angstrom)/amu to Angstrom/fs**2
fovermp = 0.009579038d0
#else
! The default units of Siesta are Bohr and fs
! Hence any conversion is simpler from that respect
fovermp = 1._dp / amu / Ang ** 2 * eV
#endif
else
#ifdef SIESTA__UNITS_ORIGINAL
C convert F/m in (Ry/Bohr)/amu to Bohr/fs**2
fovermp = 0.009579038d0 *Ang**2 / eV
#else
! The default units of Siesta are Bohr and fs, so no conversion needed
fovermp = 1._dp / amu
#endif
endif
C ........................