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