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