1
MODULE m_state_analysis
5
public :: state_analysis
9
subroutine state_analysis( istep )
11
use m_born_charge, only: born_charge
12
use parallel, only: IOnode
13
use m_wallclock, only : wallclock
14
use zmatrix, only: lUseZmatrix, iofaZmat, write_Zmatrix,
15
. CartesianForce_to_ZmatForce,
16
. write_canonical_ucell_and_Zmatrix
17
use m_iostruct, only: write_struct
19
use atomlist, only: iaorb, iphorb, iza, amass, no_u, lasto, superx
23
use m_energies, only: Etot
25
use m_spin, only: nspin
33
logical :: eggbox_block=.true. ! Read eggbox info from data file?
35
external :: eggbox, mulliken, fixed
36
real(dp), external :: volcel
38
!------------------------------------------------------------------------- BEGIN
39
! Write final Kohn-Sham Energy
40
if (cml_p) call cmlStartPropertyList(mainXML,
41
. title='Final KS Energy')
44
. write(6,"(/a,f14.4)") 'siesta: E_KS(eV) = ', Etot/eV
45
if (cml_p) call cmlAddProperty(xf=mainXML, property=Etot/eV,
46
. dictref='siesta:E_KS', units='siestaUnits:eV',
50
! Substract egg box effect form energy
51
if (eggbox_block) then
52
call eggbox('energy',ucell,na_u,isa,ntm,xa,fa,Etot,
55
. write(6,"(/a,f14.4)") 'siesta: E_KS - E_eggbox = ',Etot/eV
56
if (cml_p) call cmlAddProperty(xf=mainXML, property=Etot/eV,
57
. dictref='siesta:E_KS_egg', units='siestaUnits:eV',
60
if (cml_p) call cmlEndPropertyList(mainXML)
62
! Substract egg box effect from the forces
63
if (eggbox_block) then
64
call eggbox('forces',ucell,na_u,isa,ntm,xa,fa,Etot,eggbox_block)
67
! Compute stress without internal molecular pressure
68
call remove_intramol_pressure(ucell,stress,na_u,xa,fa,mstress)
70
C Impose constraints to atomic movements by changing forces ...........
71
if (RemoveIntraMolecularPressure) then
72
! Consider intramolecular pressure-removal as another
74
call fixed(ucell,mstress,na_u,isa, amass, xa, fa,
75
$ cstress, cfa, ntcon )
77
call fixed(ucell,stress,na_u,isa, amass, xa, fa,
78
$ cstress, cfa, ntcon )
82
! Calculate and output Zmatrix forces
84
call CartesianForce_to_ZmatForce(na_u,xa,fa)
85
if (IOnode) call iofaZmat()
88
! Compute kinetic contribution to stress
89
kin_stress(1:3,1:3) = 0.0_dp
90
volume = volcel(ucell)
94
kin_stress(ix,jx) = kin_stress(ix,jx) -
95
. amu * amass(ia) * va(ix,ia) * va(jx,ia) / volume
99
! Add kinetic term to stress tensor
100
tstress = stress + kin_stress
104
call siesta_write_forces()
105
call siesta_write_stress_pressure()
106
call wallclock('--- end of geometry step')
109
! Mulliken population analysis
110
call mulliken( mullipop, nspin, na_u, no_u, maxnh,
111
. numh, listhptr, listh, S, Dscf, isa,
112
. lasto, iaorb, iphorb )
114
! Save structural information in crystallographic format,
115
! and canonical Zmatrix (if applicable) before moving the atoms
118
call write_struct( ucell, na_u, isa, iza, xa )
119
if (lUseZmatrix) then
120
call write_canonical_ucell_and_Zmatrix()
124
! Call the born effective charge routine only in those steps (even)
125
! in which the dx is positive.
126
if (bornz .and. (mod(istep,2) .eq. 0)) then
129
!--------------------------------------------------------------------------- END
130
END subroutine state_analysis
132
END MODULE m_state_analysis