← Back to branch summary

~albertog/siesta/trunk-nlso

« back to all changes in this revision

Viewing changes to Src/sankey_change_basis.F90

  • Committer: Alberto Garcia
  • Date: 2018-06-04 16:10:46 UTC
  • mfrom: (697.2.2 trunk)
  • Revision ID: albertog@icmab.es-20180604161046-fsufw3pl21dnb5g5
Sync to trunk-699

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/sankey_change_basis.F90
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_sankey_change_basis 
9
 
 
10
 
  IMPLICIT NONE
11
 
  PRIVATE
12
 
 
13
 
  PUBLIC :: sankey_change_basis
14
 
 
15
 
  CONTAINS
16
 
! *******************************************************************
17
 
 
18
 
     SUBROUTINE sankey_change_basis ( istpmove )
19
 
 
20
 
!******************************************************************************
21
 
! This subroutine transforms the TDKS wavefunctions into new-basis during
22
 
! Ehrenfest dynamics within TDDFT after atomic positions are changed. The 
23
 
! TDKS are transformed according to the transformation proposed by Tomfohr and 
24
 
! Sankey in Ref. phys. stat. sol. (b) 226 No.1, 115-123 (2001). 
25
 
! 
26
 
! After transforming TDKS states it calculates the density matrix in new basis set. 
27
 
! We use MatrixSwitch to manipulate matrices.
28
 
!
29
 
! This subroutine is roughly based on D. Sanchez-Portal's chgbasis subroutine from, 
30
 
! now absolete, serial version of tddft-siesta.
31
 
!
32
 
! Re-written for parallelization by Adiran Garaizar and Rafi Ullah, October 2015
33
 
! Modified by Rafi Ullah, February 2017
34
 
!********************************************************************************
35
 
 
36
 
  use precision
37
 
  use parallel,            only : Node, Nodes,BlockSize, IONode 
38
 
  use parallelsubs,        only : GlobalToLocalOrb, GetNodeOrbs,               &
39
 
                                  LocalToGlobalOrb
40
 
  use m_diag_option,       only : ParallelOverK
41
 
  use fdf
42
 
  use alloc
43
 
  use sys, only: die
44
 
#ifdef MPI
45
 
  use mpi_siesta,          only : mpi_bcast, mpi_comm_world, mpi_logical
46
 
#endif
47
 
  use m_spin,              only: nspin
48
 
  use m_gamma,             only: gamma
49
 
  use Kpoint_grid,         only: nkpnt, kpoint, kweight
50
 
  use atomlist,            only: no_u, indxuo
51
 
  use wavefunctions
52
 
  use sparse_matrices,     only : numh, listhptr, listh, S, xijo, Dscf
53
 
  use MatrixSwitch
54
 
  use matdiagon,           only: geteigen 
55
 
  !
56
 
  implicit none
57
 
  !
58
 
  integer, intent(in)     :: istpmove
59
 
  !
60
 
#ifdef MPI
61
 
  integer                 :: MPIerror,desch(9)
62
 
  external                :: diagkp
63
 
#endif
64
 
  !
65
 
  logical, save           :: frstme = .true.
66
 
  integer                 :: io, iuo, juo, nuo, jo, ind, ispin
67
 
  integer                 :: ik, j,ierror
68
 
  real(dp)                :: skxij,ckxij, kxij, qe 
69
 
  complex(dp)             :: cvar1
70
 
 
71
 
  type(matrix)                     :: Maux,invsqS,phi
72
 
  type(matrix)                     :: Sauxms
73
 
  type(matrix),allocatable,save    :: sqrtS(:)
74
 
  character(3)                     :: m_operation
75
 
  character(5)                     :: m_storage
76
 
  !
77
 
#ifdef DEBUG
78
 
  call write_debug( '    PRE sankey_change_basis' )
79
 
#endif
80
 
  !
81
 
#ifdef MPI
82
 
  call GetNodeOrbs(no_u,Node,Nodes,nuo)
83
 
  if (frstme) then
84
 
    if(ParallelOverK) then
85
 
      call die ( "chgbasis: tddft-siesta not parallelized over k-points." )
86
 
    endif
87
 
      call ms_scalapack_setup(mpi_comm_world,1,'c',BlockSize)
88
 
  endif
89
 
#else
90
 
  Node = 0
91
 
  Nodes = 1
92
 
  nuo = no_u
93
 
#endif
94
 
  call timer( 'sankey_change_basis', 1 )
95
 
  !
96
 
#ifdef MPI
97
 
  m_storage='pzdbc'
98
 
  m_operation='lap'
99
 
#else
100
 
  m_storage='szden'
101
 
  m_operation='lap'
102
 
#endif
103
 
  !
104
 
  IF (nspin .eq. 4) THEN
105
 
    call die ('chgbasis: ERROR: EID not yet prepared for non-collinear spin')
106
 
  END IF
107
 
    ! Allocate local arrays
108
 
  if(frstme) then
109
 
    allocate(sqrtS(nkpnt))
110
 
    do ik=1,nkpnt
111
 
      call m_allocate(sqrtS(ik),no_u,no_u,m_storage)
112
 
    end do
113
 
    frstme=.false.
114
 
  endif
115
 
  call m_allocate(Maux,no_u,no_u,m_storage)
116
 
  call m_allocate(invsqS,no_u,no_u,m_storage)
117
 
  !
118
 
  do ik = 1,nkpnt
119
 
  call m_allocate(Sauxms,no_u,no_u,m_storage)
120
 
    do iuo = 1,nuo
121
 
      call LocalToGlobalOrb(iuo, Node, Nodes, io)
122
 
      do j = 1,numh(iuo)
123
 
        ind = listhptr(iuo) + j
124
 
        juo = listh(ind)
125
 
        jo  = indxuo (juo)
126
 
        if(.not.gamma) then 
127
 
          kxij = kpoint(1,ik) * xijo(1,ind) +&
128
 
          kpoint(2,ik) * xijo(2,ind) +&
129
 
          kpoint(3,ik) * xijo(3,ind)
130
 
          ckxij = cos(kxij)
131
 
          skxij = -sin(kxij)
132
 
        else 
133
 
          ckxij=1.0_dp
134
 
          skxij=0.0_dp
135
 
        endif
136
 
        cvar1 =  cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
137
 
        call m_set_element(Sauxms, jo, io, cvar1, complx_1, m_operation)
138
 
      enddo
139
 
    enddo
140
 
    !
141
 
    if(istpmove.eq.1) then   ! istpmove 
142
 
      ! If first step calculate S0^1/2 and save for next step. 
143
 
      call calculatesqrtS(Sauxms,invsqS,sqrtS(ik),nuo,m_storage,m_operation)
144
 
    elseif(istpmove.gt.1) then 
145
 
      ! Calculate both Sn^1/2 and Sn^-1/2 where Sn^1/2 is used in n+1 step. 
146
 
      call calculatesqrtS(Sauxms,invsqS,Maux,nuo,m_storage,m_operation)
147
 
      !Saux= Sn-1^1/2*Sn^-1/2
148
 
      call mm_multiply(invsqS,'n',sqrtS(ik),'n',&
149
 
      Sauxms,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
150
 
      ! Passing Sn^1/2 from Maux to sqrtS for next step usage.
151
 
      call m_add ( Maux,'n',sqrtS(ik),cmplx(1.0_dp,0.0_dp,dp),              &
152
 
                  cmplx(0.0_dp,0.0_dp,dp),m_operation )
153
 
      ! C1=S0^1/2*S1^1/2*C0
154
 
      qe=2.0d0*kweight(ik)/dble(nspin)
155
 
      do ispin=1,nspin
156
 
        ! Cn = Saux*Cn-1 where Saux= Sn-1^1/2*Sn^-1/2
157
 
        call m_allocate ( phi,wavef_ms(ik,ispin)%dim1,                      &
158
 
                          wavef_ms(ik,ispin)%dim2,m_storage )
159
 
        call m_add( wavef_ms(ik,ispin),'n',phi,cmplx(1.0,0.0,dp),           &
160
 
                   cmplx(0.0_dp,0.0_dp,dp),m_operation )
161
 
        call mm_multiply(Sauxms,'n',phi,'n',               &
162
 
             wavef_ms(ik,ispin),cmplx(1.0_dp,0.0_dp,dp),     &
163
 
             cmplx(0.0_dp,0.0_dp,dp),m_operation)
164
 
        call m_deallocate(phi)
165
 
      enddo  
166
 
    endif   !istpmove 
167
 
  call m_deallocate(Sauxms)
168
 
  enddo          ! ik 
169
 
  !
170
 
  IF(istpmove.gt.1) THEN   ! istpmove 
171
 
    IF (IONode) THEN
172
 
      WRITE(*,'(a)') 'chgbasis: Computing DM in new basis'
173
 
    END IF
174
 
      call compute_tddm (Dscf)
175
 
  END IF
176
 
  call m_deallocate(Maux)
177
 
  call m_deallocate(invsqS)
178
 
 
179
 
  call timer('sankey_change_basis',2)
180
 
#ifdef DEBUG
181
 
  call write_debug( '    POS sankey_change_basis' )
182
 
#endif
183
 
  END SUBROUTINE sankey_change_basis 
184
 
 
185
 
  SUBROUTINE calculatesqrtS(S,invsqS,sqrtS,nu,m_storage,m_operation)
186
 
  
187
 
    use precision 
188
 
    use matdiagon
189
 
    use MatrixSwitch
190
 
    use parallelsubs,          only: LocalToGlobalOrb
191
 
    use parallel,              only: Node, Nodes
192
 
    use wavefunctions,         only: complx_0
193
 
    ! 
194
 
    implicit none
195
 
    ! 
196
 
    character(5), intent(in)                      :: m_storage
197
 
    character(3), intent(in)                      :: m_operation
198
 
    type(matrix), intent(inout)                   :: S,invsqS,sqrtS
199
 
    type(matrix)                                  :: SD01, SD02
200
 
    complex(dp)                                   :: varaux
201
 
    real(dp)                                      :: eig01, eig02
202
 
    real(dp), allocatable                         :: eigen(:)
203
 
    integer                                       :: no,nu, info, i, j,jo
204
 
    real(dp)  tiny
205
 
    data tiny  /1.0d-10/
206
 
    ! 
207
 
    no=S%dim1 
208
 
    allocate(eigen(no))
209
 
    call m_allocate(SD01,no,no,m_storage)
210
 
    call m_allocate(SD02,no,no,m_storage)
211
 
    ! Takes overlap matrix S in dense form and returns its eigenvalues
212
 
    ! in eigen(*) and eigenvectors in S.
213
 
    call geteigen(S,eigen,m_operation)
214
 
    !
215
 
    do j=1,nu
216
 
      call LocalToGlobalOrb(j,Node,Nodes,jo)
217
 
      eig01=dsqrt(dabs(eigen(jo)))
218
 
      eig02=1.0d0/(eig01+tiny)
219
 
      do i=1,no
220
 
        varaux = S%zval(i,j)
221
 
        call m_set_element(SD01,i,jo,eig01*varaux,complx_0,m_operation)
222
 
        call m_set_element(SD02,i,jo,eig02*varaux,complx_0,m_operation)
223
 
      enddo
224
 
    enddo 
225
 
    !
226
 
    deallocate(eigen)
227
 
    
228
 
    call mm_multiply(SD01,'n',S,'c',sqrtS,cmplx(1.0_dp,0.0_dp,dp),&
229
 
                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
230
 
    call mm_multiply(SD02,'n',S,'c',invsqS,cmplx(1.0_dp,0.0_dp,dp),&
231
 
                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
232
 
    call m_deallocate(SD01)
233
 
    call m_deallocate(SD02)
234
 
    !
235
 
  END SUBROUTINE calculatesqrtS
236
 
END MODULE m_sankey_change_basis
 
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_sankey_change_basis 
 
9
 
 
10
  IMPLICIT NONE
 
11
  PRIVATE
 
12
 
 
13
  PUBLIC :: sankey_change_basis
 
14
 
 
15
  CONTAINS
 
16
! *******************************************************************
 
17
 
 
18
     SUBROUTINE sankey_change_basis ( istpmove )
 
19
 
 
20
!******************************************************************************
 
21
! This subroutine transforms the TDKS wavefunctions into new-basis during
 
22
! Ehrenfest dynamics within TDDFT after atomic positions are changed. The 
 
23
! TDKS are transformed according to the transformation proposed by Tomfohr and 
 
24
! Sankey in Ref. phys. stat. sol. (b) 226 No.1, 115-123 (2001). 
 
25
! 
 
26
! After transforming TDKS states it calculates the density matrix in new basis set. 
 
27
! We use MatrixSwitch to manipulate matrices.
 
28
!
 
29
! This subroutine is roughly based on D. Sanchez-Portal's chgbasis subroutine from, 
 
30
! now absolete, serial version of tddft-siesta.
 
31
!
 
32
! Re-written for parallelization by Adiran Garaizar and Rafi Ullah, October 2015
 
33
! Modified by Rafi Ullah, February 2017
 
34
!********************************************************************************
 
35
 
 
36
  use precision
 
37
  use parallel,            only : Node, Nodes,BlockSize, IONode 
 
38
  use parallelsubs,        only : GlobalToLocalOrb, GetNodeOrbs,               &
 
39
                                  LocalToGlobalOrb
 
40
  use m_diag_option,       only : ParallelOverK
 
41
  use fdf
 
42
  use alloc
 
43
  use sys, only: die
 
44
#ifdef MPI
 
45
  use mpi_siesta,          only : mpi_bcast, mpi_comm_world, mpi_logical
 
46
#endif
 
47
  use m_spin,              only: nspin
 
48
  use m_gamma,             only: gamma
 
49
  use kpoint_scf_m,        only: kpoints_scf
 
50
  use atomlist,            only: no_u, indxuo
 
51
  use wavefunctions
 
52
  use sparse_matrices,     only : numh, listhptr, listh, S, xijo, Dscf
 
53
  use MatrixSwitch
 
54
  use matdiagon,           only: geteigen 
 
55
  !
 
56
  implicit none
 
57
  !
 
58
  integer, intent(in)     :: istpmove
 
59
  !
 
60
#ifdef MPI
 
61
  integer                 :: MPIerror,desch(9)
 
62
  external                :: diagkp
 
63
#endif
 
64
  !
 
65
  logical, save           :: frstme = .true.
 
66
  integer                 :: io, iuo, juo, nuo, jo, ind, ispin
 
67
  integer                 :: ik, j,ierror
 
68
  real(dp)                :: skxij,ckxij, kxij, qe 
 
69
  complex(dp)             :: cvar1
 
70
 
 
71
  type(matrix)                     :: Maux,invsqS,phi
 
72
  type(matrix)                     :: Sauxms
 
73
  type(matrix),allocatable,save    :: sqrtS(:)
 
74
  character(3)                     :: m_operation
 
75
  character(5)                     :: m_storage
 
76
  !
 
77
#ifdef DEBUG
 
78
  call write_debug( '    PRE sankey_change_basis' )
 
79
#endif
 
80
  !
 
81
#ifdef MPI
 
82
  call GetNodeOrbs(no_u,Node,Nodes,nuo)
 
83
  if (frstme) then
 
84
    if(ParallelOverK) then
 
85
      call die ( "chgbasis: tddft-siesta not parallelized over k-points." )
 
86
    endif
 
87
      call ms_scalapack_setup(mpi_comm_world,1,'c',BlockSize)
 
88
  endif
 
89
#else
 
90
  Node = 0
 
91
  Nodes = 1
 
92
  nuo = no_u
 
93
#endif
 
94
  call timer( 'sankey_change_basis', 1 )
 
95
  !
 
96
#ifdef MPI
 
97
  m_storage='pzdbc'
 
98
  m_operation='lap'
 
99
#else
 
100
  m_storage='szden'
 
101
  m_operation='lap'
 
102
#endif
 
103
  !
 
104
  IF (nspin .eq. 4) THEN
 
105
    call die ('chgbasis: ERROR: EID not yet prepared for non-collinear spin')
 
106
  END IF
 
107
    ! Allocate local arrays
 
108
  if(frstme) then
 
109
    allocate(sqrtS(kpoints_scf%N))
 
110
    do ik=1,kpoints_scf%N
 
111
      call m_allocate(sqrtS(ik),no_u,no_u,m_storage)
 
112
    end do
 
113
    frstme=.false.
 
114
  endif
 
115
  call m_allocate(Maux,no_u,no_u,m_storage)
 
116
  call m_allocate(invsqS,no_u,no_u,m_storage)
 
117
  !
 
118
  do ik = 1,kpoints_scf%N
 
119
  call m_allocate(Sauxms,no_u,no_u,m_storage)
 
120
    do iuo = 1,nuo
 
121
      call LocalToGlobalOrb(iuo, Node, Nodes, io)
 
122
      do j = 1,numh(iuo)
 
123
        ind = listhptr(iuo) + j
 
124
        juo = listh(ind)
 
125
        jo  = indxuo (juo)
 
126
        if(.not.gamma) then 
 
127
          kxij = kpoints_scf%k(1,ik) * xijo(1,ind) +&
 
128
          kpoints_scf%k(2,ik) * xijo(2,ind) +&
 
129
          kpoints_scf%k(3,ik) * xijo(3,ind)
 
130
          ckxij = cos(kxij)
 
131
          skxij = -sin(kxij)
 
132
        else 
 
133
          ckxij=1.0_dp
 
134
          skxij=0.0_dp
 
135
        endif
 
136
        cvar1 =  cmplx(S(ind)*ckxij,S(ind)*skxij,dp)
 
137
        call m_set_element(Sauxms, jo, io, cvar1, complx_1, m_operation)
 
138
      enddo
 
139
    enddo
 
140
    !
 
141
    if(istpmove.eq.1) then   ! istpmove 
 
142
      ! If first step calculate S0^1/2 and save for next step. 
 
143
      call calculatesqrtS(Sauxms,invsqS,sqrtS(ik),nuo,m_storage,m_operation)
 
144
    elseif(istpmove.gt.1) then 
 
145
      ! Calculate both Sn^1/2 and Sn^-1/2 where Sn^1/2 is used in n+1 step. 
 
146
      call calculatesqrtS(Sauxms,invsqS,Maux,nuo,m_storage,m_operation)
 
147
      !Saux= Sn-1^1/2*Sn^-1/2
 
148
      call mm_multiply(invsqS,'n',sqrtS(ik),'n',&
 
149
      Sauxms,cmplx(1.0_dp,0.0_dp,dp),cmplx(0.0_dp,0.0_dp,dp),m_operation)
 
150
      ! Passing Sn^1/2 from Maux to sqrtS for next step usage.
 
151
      call m_add ( Maux,'n',sqrtS(ik),cmplx(1.0_dp,0.0_dp,dp),              &
 
152
                  cmplx(0.0_dp,0.0_dp,dp),m_operation )
 
153
      ! C1=S0^1/2*S1^1/2*C0
 
154
      qe=2.0d0*kpoints_scf%w(ik)/dble(nspin)
 
155
      do ispin=1,nspin
 
156
        ! Cn = Saux*Cn-1 where Saux= Sn-1^1/2*Sn^-1/2
 
157
        call m_allocate ( phi,wavef_ms(ik,ispin)%dim1,                      &
 
158
                          wavef_ms(ik,ispin)%dim2,m_storage )
 
159
        call m_add( wavef_ms(ik,ispin),'n',phi,cmplx(1.0,0.0,dp),           &
 
160
                   cmplx(0.0_dp,0.0_dp,dp),m_operation )
 
161
        call mm_multiply(Sauxms,'n',phi,'n',               &
 
162
             wavef_ms(ik,ispin),cmplx(1.0_dp,0.0_dp,dp),     &
 
163
             cmplx(0.0_dp,0.0_dp,dp),m_operation)
 
164
        call m_deallocate(phi)
 
165
      enddo  
 
166
    endif   !istpmove 
 
167
  call m_deallocate(Sauxms)
 
168
  enddo          ! ik 
 
169
  !
 
170
  IF(istpmove.gt.1) THEN   ! istpmove 
 
171
    IF (IONode) THEN
 
172
      WRITE(*,'(a)') 'chgbasis: Computing DM in new basis'
 
173
    END IF
 
174
      call compute_tddm (Dscf)
 
175
  END IF
 
176
  call m_deallocate(Maux)
 
177
  call m_deallocate(invsqS)
 
178
 
 
179
  call timer('sankey_change_basis',2)
 
180
#ifdef DEBUG
 
181
  call write_debug( '    POS sankey_change_basis' )
 
182
#endif
 
183
  END SUBROUTINE sankey_change_basis 
 
184
 
 
185
  SUBROUTINE calculatesqrtS(S,invsqS,sqrtS,nu,m_storage,m_operation)
 
186
  
 
187
    use precision 
 
188
    use matdiagon
 
189
    use MatrixSwitch
 
190
    use parallelsubs,          only: LocalToGlobalOrb
 
191
    use parallel,              only: Node, Nodes
 
192
    use wavefunctions,         only: complx_0
 
193
    ! 
 
194
    implicit none
 
195
    ! 
 
196
    character(5), intent(in)                      :: m_storage
 
197
    character(3), intent(in)                      :: m_operation
 
198
    type(matrix), intent(inout)                   :: S,invsqS,sqrtS
 
199
    type(matrix)                                  :: SD01, SD02
 
200
    complex(dp)                                   :: varaux
 
201
    real(dp)                                      :: eig01, eig02
 
202
    real(dp), allocatable                         :: eigen(:)
 
203
    integer                                       :: no,nu, info, i, j,jo
 
204
    real(dp)  tiny
 
205
    data tiny  /1.0d-10/
 
206
    ! 
 
207
    no=S%dim1 
 
208
    allocate(eigen(no))
 
209
    call m_allocate(SD01,no,no,m_storage)
 
210
    call m_allocate(SD02,no,no,m_storage)
 
211
    ! Takes overlap matrix S in dense form and returns its eigenvalues
 
212
    ! in eigen(*) and eigenvectors in S.
 
213
    call geteigen(S,eigen,m_operation)
 
214
    !
 
215
    do j=1,nu
 
216
      call LocalToGlobalOrb(j,Node,Nodes,jo)
 
217
      eig01=dsqrt(dabs(eigen(jo)))
 
218
      eig02=1.0d0/(eig01+tiny)
 
219
      do i=1,no
 
220
        varaux = S%zval(i,j)
 
221
        call m_set_element(SD01,i,jo,eig01*varaux,complx_0,m_operation)
 
222
        call m_set_element(SD02,i,jo,eig02*varaux,complx_0,m_operation)
 
223
      enddo
 
224
    enddo 
 
225
    !
 
226
    deallocate(eigen)
 
227
    
 
228
    call mm_multiply(SD01,'n',S,'c',sqrtS,cmplx(1.0_dp,0.0_dp,dp),&
 
229
                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
 
230
    call mm_multiply(SD02,'n',S,'c',invsqS,cmplx(1.0_dp,0.0_dp,dp),&
 
231
                     cmplx(0.0_dp,0.0_dp,dp),m_operation)
 
232
    call m_deallocate(SD01)
 
233
    call m_deallocate(SD02)
 
234
    !
 
235
  END SUBROUTINE calculatesqrtS
 
236
END MODULE m_sankey_change_basis