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      subroutine extrapolon(istep,iord,nspin,nbasis,nbasisCloc,maxnc,
9     .                      numc,listc,aux,numcold,listcold,cold,c,
10     .                      nbasisloc)
11C ******************************************************************************
12C Subroutine to extrapolate a given matrix M (like the coefficients of the
13C wave functions, or the density matrix) for the next MD step.
14C The matrix M is given in sparse form.
15C
16C Writen by P.Ordejon, November'96.
17C ******************************* INPUT ***************************************
18C integer istep                : Time step of the simulation
19C integer iord                 : Extrapolation order (0 or 1)
20C                                0 = 0th order;  1 = 1st order
21C integer nspin                : Number of spin polarizations (1 or 2)
22C integer nbasis               : Number of rows of matrix M
23C integer nbasisloc            : Number of rows of matrix M held locally
24C integer nbasisCloc           : Maximum number of rows of matrix M (dimension)
25C integer maxnc                : First dimension of M matrix, and maximum
26C                                number of nonzero elements of each column of M
27C integer numc(maxnc)          : Control vector 1 of M matrix at t
28C integer listc(maxnc,nbasisCloc) : Control vector 2 of M matrix at t
29C real*8 aux(2,nbasis)         : Auxiliary storage array
30C ************************** INPUT AND OUTPUT *********************************
31C integer numcold(maxnc)         : Input: Control vector 1 of M matrix at t-dt
32C                                       (if istep .ne. 1)
33C                                Output: Control vector 1 of M matrix at t
34C integer listcold(maxnc,nbasisCloc) : Input: Control vector 2 of M matrix at t-dt
35C                                       (if istep .ne. 1)
36C                                Output: Control vector 2 of M matrix at t
37C real*8 cold(maxnc,nbasisCloc,nspin): Input: matrix M at t-2dt
38C                                Output: matrix M at t-dt
39C real*8 c(maxnc,nbasisCloc,nspin): New matrix M (extrapolated)
40C                                Input: matrix at t-dt
41C                                Output: matrix at t
42C                                If istep = 1, c returned uncahanged
43C **************************** SCRATCH SPACE **********************************
44C real*8 aux(2,nbasis)         : Auxiliary storage array
45C **************************** BEHAVIOUR **************************************
46C The routine allows for the sparse structure of the matrix M to change
47C between MD time steps. On input, the matrices of former steps (c and cold)
48C have the structure of last step (t-dt): numold and listold; whereas the new
49C (extrapolated) matrix has the structure of the current time step (which
50C must be determined before calling this routine!!): num and list.
51C On output, the routine updates the structure of c and cold, to that
52C at the current (t) time steps respectively. Same with numold and listold
53C
54C For the first MD time step (istep = 1), there is no extrapolation.
55C In that case, c is returned unchanged.
56C Also, in that case numold and listold are only an output, and are set equal
57C to num and list
58C *****************************************************************************
59
60      use precision, only : dp
61      use parallel,  only : IONode
62      use sys,       only : die
63
64      implicit none
65
66      integer, intent(in) :: iord
67      integer, intent(in) :: istep
68      integer, intent(in) :: nspin
69      integer, intent(in) :: nbasis
70      integer, intent(in) :: nbasisloc
71      integer, intent(in) :: nbasisCloc
72      integer, intent(in) :: maxnc
73      integer, intent(in) :: numc(nbasisCloc)
74      integer, intent(in) :: listc(maxnc, nbasisCloc)
75
76      integer, intent(inout) :: numcold(nbasisloc)
77      integer, intent(inout) :: listcold(maxnc, nbasisloc)
78      real(dp), intent(inout) :: cold(maxnc, nbasisloc, nspin)
79      real(dp), intent(inout) :: c(maxnc, nbasisCloc, nspin)
80
81      real(dp), intent(out) :: aux(2,nbasis)
82C  Internal variables .......................................................
83
84      integer
85     .  i,in,ispin,j
86
87      real(dp) ::
88     .  msave
89
90      logical ::
91     .  changed
92C ...........................................................................
93
94      if (iord /= 0 .and. iord /= 1) then
95        if (IONode) then
96          call die ('extrapolon: Wrong iord: '//
97     .              'only 0 and 1 order available')
98        endif
99      endif
100
101C Just initialize numcold and listcold if istep = 1 ...........................
102      if (istep .eq. 1) then
103        do i = 1,nbasisloc
104          numcold(i) = numc(i)
105          do in = 1,numc(i)
106            listcold(in,i) = listc(in,i)
107            do ispin = 1,nspin
108              cold(in,i,ispin) = 0.0_dp
109            enddo
110          enddo
111        enddo
112
113      else
114
115C Check if sparse structure has changed .....................................
116        changed = .false.
117
118        do i = 1,nbasisloc
119          if (numcold(i).ne.numc(i)) then
120            changed = .true.
121            exit
122          endif
123          do in = 1,numc(i)
124            if (listcold(in,i).ne.listc(in,i)) then
125              changed = .true.
126              exit
127            endif
128          enddo
129        enddo
130
131C If sparse structure has changed, re-order c and cold
132C and change numcold and listcold to current ones .............................
133
134        if (changed) then
135          aux = 0.0_dp
136
137          do i = 1,nbasisloc
138            do ispin = 1,nspin
139              do in = 1,numcold(i)
140                j = listcold(in,i)
141                aux(1,j) = c(in,i,ispin)
142                aux(2,j) = cold(in,i,ispin)
143              enddo
144              do in = 1,numc(i)
145                j = listc(in,i)
146                c(in,i,ispin) = aux(1,j)
147                cold(in,i,ispin) = aux(2,j)
148              enddo
149              do in = 1,numcold(i)
150                j = listcold(in,i)
151                aux(1,j) = 0.0_dp
152                aux(2,j) = 0.0_dp
153              enddo
154            enddo
155            numcold(i) = numc(i)
156            do in = 1,numc(i)
157              listcold(in,i) = listc(in,i)
158            enddo
159          enddo
160        endif !changed
161
162C Extrapolate matrix M ......................................................
163
164        do ispin = 1,nspin
165          do i = 1,nbasisloc
166            do in = 1,numc(i)
167              msave = c(in,i,ispin)
168              if (iord == 1) then
169                c(in,i,ispin) = 2.0_dp*c(in,i,ispin) - cold(in,i,ispin)
170              endif
171              cold(in,i,ispin) = msave
172            enddo
173          enddo
174        enddo
175
176      endif !istep==1
177
178      end
179