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