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_eandg 9 public :: eandg 10 CONTAINS 11 subroutine eandg(iopt,eta,qs,lam,nhmax,numh,listh,listhptr, 12 . ncmax,numc,listc,h,s,c,no_u,no_l,nbands, 13 . e3,e,grad,dm,edm,no_cl,nspin,Node,Nodes) 14C ************************************************************************** 15C Subroutine to link the CG algorithms to the calculation of 16C the functional energy and gradients, and to set up control 17C vectors for auxiliary sparse matrices. 18C It also computes the density matrix. 19C This routine works with the funcional of Kim et al. (PRB 52, 1640 (95)). 20C Written by P.Ordejon. October'96 21C ****************************** INPUT ************************************* 22C integer iopt : Input option parameter 23C iopt = 0 => Set up control vectors 24C iopt = 1 => Call energy routine for line. min. 25C iopt = 2 => Call gradient routine 26C iopt = 3 => Call density matrix routine 27C real*8 eta : Fermi level parameter of Kim et al 28C real*8 qs(nspin) : Total number of electrons 29C real*8 lam : Length of step for line minimization 30C integer nhmax : First dimension of listh and H, and maximum 31C number of nonzero elements of each row of H 32C integer numh(no_u) : Control vector of H matrix 33C (number of <>0 element of each row) 34C integer listh(nhmax) : Control vector of H matrix 35C (list of <>0 element of each row) 36C integer listhptr(no_u) : Control vector of H matrix 37C (pointer to start of row in listh/H) 38C integer ncmax : First dimension of listc and C, and maximum 39C number of nonzero elements of each row of C 40C integer numc(no_u) : Control vector of C matrix 41C (number of <>0 element of each row) 42C integer listc(ncmax,no_u): Control vector of C matrix 43C (list of <>0 element of each row) 44C real*8 h(nhmax) : Hamiltonian in sparse form 45C real*8 s(nhmax) : Overlap in sparse form 46C real*8 c(ncmax,no_u) : Current point (wave func. coeffs. in sparse) 47C integer no_u : Global number of atomic orbitals 48C integer no_l : Local number of atomic orbitals 49C integer nbands : Number of Localized Wave Functions 50C ******** INPUT OR OUTPUT (DEPENDING ON ARGUMENT IOPT) ******************** 51C real*8 grad(ncmax,no_u) : Gradient of the functional 52C (input if iopt = 1) 53C (output if iopt = 2) 54C ***************************** OUTPUT ************************************* 55C real*8 e3(3) : Value of the energy in three points C+LAM_i*GRAD 56C real*8 e : Value of the energy at point C 57C real*8 dm(nhmax) : Density matrix in sparse form 58C real*8 edm(nhmax) : Energy density matrix in sparse form 59C **************************** BEHAVIOUR *********************************** 60C The overlap matrix 'o' must be in the same sparse format as the 61C Hamiltonian matrix 'h', even if the overlap is more sparse than h 62C (as due to the KB projectors, for instance). It will, in general, 63C contain some zeros, therefore. 64C ************************************************************************** 65 66 use precision 67 use alloc 68 use fdf 69 use on_core, only : cttoc, fttof, listct, listf 70 use on_core, only : listft, listhij, numct, numf 71 use on_core, only : numft, numhij, indon, nindv 72 use on_core, only : f, fs 73 74 use on_main, only : nbL2G, nbG2L, nbandsloc, ncP2T 75 use sys, only : die 76 use m_gradient, only: gradient 77 use m_ener3, only: ener3 78 use m_denmat, only: denmat 79 use m_on_subs, only: axb_build1,axb_build2,ctrans1,ctrans2 80#ifdef MPI 81 use globalise, only : setglobaliseB, setglobalise, setglobaliseF 82 use globalise, only : globalinitB 83#endif 84 85 implicit none 86 87 integer, intent(in) :: 88 . iopt,nbands,no_u,no_l,no_cl,ncmax,nhmax,nspin, 89 . Node,Nodes 90 91 integer, intent(in) :: 92 . listc(ncmax,no_cl),listh(nhmax),listhptr(no_l), 93 . numc(no_u),numh(no_l) 94 95 real(dp), intent(in) :: 96 . c(ncmax,no_cl,nspin),qs(nspin),eta(nspin),h(nhmax,nspin), 97 . lam,s(nhmax) 98 99 real(dp), intent(inout) :: grad(ncmax,no_cl,nspin) 100 101 real(dp), intent(out) :: e,e3(3) 102 real(dp), intent(out) :: dm(nhmax,nspin), edm(nhmax,nspin) 103 104 105C Internal variables ..................................................... 106 integer, save :: maxnft = 1 107 integer, save :: maxnct = 1 108 integer, save :: maxnf = 1 109 integer, save :: maxnhij = 1 110 111 integer :: i 112 integer :: m, mt, n 113 114 logical, pointer, save :: lNeeded(:) 115 logical, save :: frstme = .true. 116 logical, save :: frstme2 = .true. 117 logical, save :: LoMem = .false. 118 119C ..................... 120 121* call timer('eandg',1) 122 123 if (frstme) then 124C Nullify pointers 125 nullify(cttoc) 126 nullify(fttof) 127 nullify(listct) 128 nullify(listf) 129 nullify(listft) 130 nullify(listhij) 131 nullify(numct) 132 nullify(numf) 133 nullify(numft) 134 nullify(numhij) 135 nullify(indon) 136 nullify(nindv) 137 nullify(nbL2G) 138 nullify(nbG2L) 139 nullify(f) 140 nullify(fs) 141 142C Set flag for memory algorithm choice 143 LoMem = fdf_boolean('ON.LowerMemory',.false.) 144 145 frstme = .false. 146 147 endif 148 149 if (iopt.eq.0) then 150C Find nbandsloc and pointers 151 call re_alloc( nbL2G, 1, nbands, 'nbL2G', 'eandg' ) 152 call re_alloc( nbG2L, 1, nbands, 'nbG2L', 'eandg' ) 153 call re_alloc( lNeeded, 1, nbands, 'lNeeded', 'eandg' ) 154 do i = 1,nbands 155 lNeeded(i) = .false. 156 enddo 157 do m = 1,no_l 158 mt = ncP2T(m) 159 do n = 1,numc(mt) 160 lNeeded(listc(n,mt)) = .true. 161 enddo 162 enddo 163 nbandsloc = 0 164 nbG2L(1:nbands) = 0 165 do i = 1,nbands 166 if (lNeeded(i)) then 167 nbandsloc = nbandsloc + 1 168 nbL2G(nbandsloc) = i 169 nbG2L(i) = nbandsloc 170 endif 171 enddo 172 deallocate(lNeeded) 173 174C Allocation of memory 175 call re_alloc(numct,1,nbandsloc,name='numct',copy=.true., 176 . shrink=.false.) 177 call re_alloc(numf,1,nbandsloc,name='numf',copy=.true., 178 . shrink=.false.) 179 call re_alloc(numft,1,no_cl,name='numft',copy=.true., 180 . shrink=.false.) 181 call re_alloc(numhij,1,nbandsloc,name='numhij',copy=.true., 182 . shrink=.false.) 183 call re_alloc(indon,1,no_u,name='indon',copy=.true., 184 . shrink=.false.) 185 call re_alloc(nindv,1,no_u,name='nindv',copy=.true., 186 . shrink=.false.) 187 call re_alloc(cttoc,1,maxnct,1,nbandsloc,name='cttoc', 188 . copy=.true.,shrink=.false.) 189 call re_alloc(fttof,1,maxnft,1,no_cl,name='fttof', 190 . copy=.true.,shrink=.false.) 191 call re_alloc(listct,1,maxnct,1,nbandsloc,name='listct', 192 . copy=.true.,shrink=.false.) 193 call re_alloc(listf,1,maxnf,1,nbandsloc,name='listf', 194 . copy=.true.,shrink=.false.) 195 call re_alloc(listft,1,maxnft,1,no_cl,name='listft', 196 . copy=.true.,shrink=.false.) 197 call re_alloc(listhij,1,maxnhij,1,nbandsloc,name='listhij', 198 . copy=.true.,shrink=.false.) 199 call re_alloc(f,1,maxnf,1,nbandsloc,name='f',copy=.true., 200 . shrink=.false.) 201 call re_alloc(fs,1,maxnf,1,nbandsloc,name='fs',copy=.true., 202 . shrink=.false.) 203 204 205 206C Set up index lists for sparse matrices .................................. 207 208C GET Ct LISTS 209 call ctrans1(no_cl,maxnct) 210 211C GET F LISTS 212 call axb_build1(no_u,maxnct,numct,listct,nhmax,numh,listh, 213 . listhptr,indon,maxnf) 214 215C Allocate arrays that depend on variables set above 216 call re_alloc( listf, 1, maxnf, 1, nbandsloc, 'listf', 'eandg' ) 217 218C GET Ft LISTS 219 call ctrans2(no_cl,maxnf,maxnft,numf,listf) 220 221C GET Hij LISTS 222 call axb_build2(nbands,maxnf,numf,listf, 223 . no_cl,ncmax,numc,listc,indon,maxnhij) 224 225C Allocate arrays that depend on variables set above 226 call re_alloc( f, 1, maxnf, 1, nbandsloc, 'f','eandg' ) 227 call re_alloc( fs, 1, maxnf, 1, nbandsloc, 'fs','eandg' ) 228 229#ifdef MPI 230C Set up communication information 231 call setglobalise(no_u,no_l,no_cl,nhmax,numh,listh, 232 . listhptr,Node,Nodes) 233 call setglobaliseB(nbands,Node) 234 call setglobaliseF(no_u,nbands,no_cl,maxnft,Node) 235#endif 236 goto 999 237 endif 238C......................... 239 240 if (frstme2) then 241#ifdef MPI 242C Initialise communication arrays - call here to avoid increasing size of arrays later 243 if (LoMem) then 244 call globalinitB(1) 245 else 246 call globalinitB(3) 247 endif 248#endif 249 frstme2 = .false. 250 endif 251 252C Calculate the energy at three points of the line, for the 253C CG line minimization ..................................................... 254 if (iopt .eq. 1) then 255 if (LoMem) then 256 call ener3lomem(c,grad,lam,eta,qs,h,s,no_u,no_l,nbands, 257 . ncmax,maxnct,maxnf,nhmax,maxnhij,numc,listc, 258 . numct,listct,cttoc,numf,listf,numh,listh, 259 . listhptr,numhij,listhij,e3,no_cl,nspin, 260 . Node) 261 else 262 call ener3(c,grad,lam,eta,qs,h,s,no_u,no_l,nbands, 263 . ncmax,maxnct,maxnf,nhmax,maxnhij,numc,listc,numct, 264 . listct,cttoc,numf,listf,numh,listh,listhptr, 265 . numhij,listhij,e3,no_cl,nspin,Node) 266 endif 267 goto 999 268 endif 269C......................... 270 271C Calculate the energy and Gradient at current point ....................... 272 if (iopt .eq. 2) then 273 if (LoMem) then 274 call gradientlomem(c,eta,qs,h,s,no_u,no_l,nbands,ncmax, 275 . maxnct,maxnf,maxnft,nhmax,maxnhij,numc, 276 . listc,numct,listct,cttoc,numf,listf,numft, 277 . listft,fttof,numh,listh,listhptr,numhij, 278 . listhij,f,fs,grad,e,no_cl,nspin,Node) 279 else 280 call gradient(c,eta,qs,h,s,no_u,no_l,nbands,ncmax, 281 . maxnct,maxnf,maxnft,nhmax,maxnhij,numc, 282 . listc,numct,listct,cttoc,numf,listf,numft, 283 . listft,fttof,numh,listh,listhptr,numhij, 284 . listhij,f,fs,grad,e,no_cl,nspin,Node) 285 endif 286 goto 999 287 endif 288C......................... 289 290C Calculate density matrix ................................................. 291C-JMS Modified denmat argument list 292 if (iopt .eq. 3) then 293 if (LoMem) then 294 call denmatlomem(c,eta,h,s,qs,no_u,no_l,nbands,ncmax, 295 . maxnct,maxnf,maxnft,nhmax,maxnhij,numc,listc, 296 . numct,listct,cttoc,numf,listf,numft,listft, 297 . numh,listh,listhptr,numhij,listhij,dm,edm, 298 . no_cl,nspin,Node) 299 else 300 call denmat(c,eta,h,s,qs,no_u,no_l,nbands,ncmax, 301 . maxnct,maxnf,maxnft,nhmax,maxnhij,numc,listc, 302 . numct,listct,cttoc,numf,listf,numft,listft, 303 . numh,listh,listhptr,numhij,listhij,dm,edm, 304 . no_cl,nspin,Node) 305 endif 306 goto 999 307 endif 308C......................... 309 call die('Error in eandg: incorrect iopt') 310 311 999 continue 312* call timer('eandg',2) 313 end subroutine eandg 314 end module m_eandg 315 316