1! ---
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt .
6! See Docs/Contributors.txt for a list of contributors.
7! ---
8      MODULE m_state_analysis
9      use write_subs
10
11      private
12      public :: state_analysis
13
14      CONTAINS
15
16      subroutine state_analysis( istep )
17      use siesta_cml
18      use m_born_charge, only : born_charge
19      use parallel,      only : IOnode
20      use m_wallclock,   only : wallclock
21      use zmatrix,       only : lUseZmatrix, iofaZmat,
22     &                          CartesianForce_to_ZmatForce
23      use atomlist,      only : iaorb, iphorb, amass, no_u, lasto
24      use atomlist,      only : indxuo, Qtot
25      use m_spin,        only : nspin, SpOrb
26      use m_fixed,       only : fixed
27      use sparse_matrices
28      use siesta_geom
29
30      USE siesta_options
31      use units, only: amu, eV
32      use m_stress
33      use m_energies, only: Etot, FreeE, Eharrs, FreeEHarris, Entropy
34      use m_energies, only: Ebs, Ef
35      use m_ntm
36      use m_forces
37      use m_energies, only: update_FreeE, update_FreeEHarris
38      use m_intramol_pressure, only: remove_intramol_pressure
39      use m_ts_global_vars, only: TSrun
40      use m_ts_options, only: N_elec, Elecs
41      use ts_charge_m, only: ts_charge_print, TS_Q_INFO_FULL
42
43#ifdef SIESTA__FLOOK
44      use flook_siesta, only : slua_call, LUA_FORCES
45#endif
46
47      implicit none
48
49      integer  :: istep
50      integer  :: ia, jx, ix
51      real(dp) :: volume
52      logical  :: eggbox_block=.true. ! Read eggbox info from data file?
53      real(dp) :: qspin
54
55      external :: eggbox, mulliken, moments
56      real(dp), external :: volcel
57
58!------------------------------------------------------------------------- BEGIN
59      call timer( 'state_analysis', 1 )
60#ifdef DEBUG
61      call write_debug( '  PRE state_analysis' )
62#endif
63
64      if (cml_p) then
65        call cmlStartModule(xf=mainXML, title='SCF Finalization')
66      endif
67
68!     Write final Kohn-Sham and Free Energy
69
70      FreeE       = Etot - Temp * Entropy
71      FreeEHarris = Eharrs - Temp * Entropy
72
73      if (cml_p) call cmlStartPropertyList(mainXML,
74     &                        title='Energies and spin')
75      if (IOnode) then
76        if ( .not. harrisfun)
77     &      write(6,"(/a,f14.4)")  'siesta: E_KS(eV) =        ', Etot/eV
78        if (cml_p) then
79           call cmlAddProperty(xf=mainXML, value=Etot/eV,
80     &       dictref='siesta:E_KS', units='siestaUnits:eV',
81     .       fmt='r6')
82           call cmlAddProperty(xf=mainXML, value=FreeE/eV,
83     &       dictref='siesta:FreeE', units='siestaUnits:eV',
84     .       fmt='r6')
85           call cmlAddProperty(xf=mainXML, value=Ebs/eV,
86     &       dictref='siesta:Ebs', units='siestaUnits:eV',
87     .       fmt='r6')
88           call cmlAddProperty(xf=mainXML, value=Ef/eV,
89     &       dictref='siesta:E_Fermi', units='siestaUnits:eV',
90     .       fmt='r6')
91        endif
92      endif
93
94!     Substract egg box effect from energy
95      if (eggbox_block) then
96        call eggbox( 'energy', ucell, na_u, isa, ntm, xa, fa, Etot,
97     &               eggbox_block )
98        FreeE  = Etot - Temp * Entropy
99        if (IOnode)
100     &    write(6,"(/a,f14.4)") 'siesta: E_KS - E_eggbox = ',Etot/eV
101        if (cml_p) call cmlAddProperty(xf=mainXML, value=Etot/eV,
102     &         dictref='siesta:E_KS_egg', units='siestaUnits:eV',
103     .         fmt='r6')
104      endif
105
106      call update_FreeE( Temp )
107      call update_FreeEHarris( Temp )
108      call print_spin(qspin)
109
110      if (cml_p) call cmlEndPropertyList( mainXML )
111
112!     Substract egg box effect from the forces
113      if (eggbox_block) then
114        call eggbox('forces',ucell,na_u,isa,ntm,xa,fa,Etot,eggbox_block)
115      endif
116
117      if (IOnode) call write_raw_efs(stress,na_u,fa,FreeE)
118
119!     Compute stress without internal molecular pressure
120      call remove_intramol_pressure(ucell,stress,na_u,xa,fa,mstress)
121
122!     Impose constraints to atomic movements by changing forces ...........
123      if (RemoveIntraMolecularPressure) then
124!        Consider intramolecular pressure-removal as another
125!        kind of constraint
126         call fixed( ucell, mstress, na_u, isa, amass, xa, fa,
127     &               cstress, cfa, ntcon ,
128     &               magnitude_usage = idyn==0 )
129      else
130         call fixed( ucell, stress, na_u, isa, amass, xa, fa,
131     &               cstress, cfa, ntcon ,
132     &               magnitude_usage = idyn==0 )
133      endif
134
135#ifdef SIESTA__FLOOK
136      ! We call it right after using the
137      ! geometry constraints.
138      ! In that way we can use both methods on top
139      ! of each other!
140      ! The easy, already implemented methods in fixed,
141      ! and custom ones in Lua :)
142      call slua_call(LUA, LUA_FORCES)
143#endif
144
145!     Calculate and output Zmatrix forces
146      if (lUseZmatrix .and. (idyn.eq.0)) then
147        call CartesianForce_to_ZmatForce(na_u,xa,fa)
148        if (IOnode) call iofaZmat()
149      endif
150
151!     Compute kinetic contribution to stress
152      kin_stress(1:3,1:3) = 0.0_dp
153      volume = volcel(ucell)
154      do ia = 1,na_u
155        do jx = 1,3
156          do ix = 1,3
157            kin_stress(ix,jx) = kin_stress(ix,jx) -
158     &             amu * amass(ia) * va(ix,ia) * va(jx,ia) / volume
159          enddo
160        enddo
161      enddo
162!     Add kinetic term to stress tensor
163      tstress = stress + kin_stress
164
165!     Force output
166      if (IOnode) then
167        call siesta_write_forces(istep)
168        if ( TSrun ) then
169          call transiesta_write_forces()
170        end if
171        call siesta_write_stress_pressure()
172        call wallclock('--- end of geometry step')
173      endif
174
175!     Population and moment analysis
176      if ( SpOrb .and. orbmoms) then
177         call moments( 1, na_u, no_u, maxnh, numh, listhptr,
178     .           listh, S, Dscf, isa, lasto, iaorb, iphorb,
179     .           indxuo )
180      endif
181      ! Call this unconditionally
182      call mulliken( mullipop, na_u, no_u, maxnh,
183     &               numh, listhptr, listh, S, Dscf, isa,
184     &               lasto, iaorb, iphorb )
185
186      ! Also write TS charges
187      if ( TSrun ) then
188        call ts_charge_print(N_Elec,Elecs, Qtot,
189     &      block_dist, sparse_pattern,
190     &      nspin, maxnh, Dscf, S, TS_Q_INFO_FULL)
191      end if
192
193!
194!     Call the born effective charge routine only in those steps (even)
195!     in which the dx  is positive.
196      if (bornz .and. (mod(istep,2) .eq. 0)) then
197         call born_charge()
198      endif
199
200!     End the xml module corresponding to the analysis
201      if (cml_p) then
202         call cmlEndModule(mainXML)
203      endif
204      call timer( 'state_analysis', 2 )
205
206!--------------------------------------------------------------------------- END
207      END subroutine state_analysis
208
209      END MODULE m_state_analysis
210