← Back to branch summary

~siesta-pseudos-bases/siesta/trunk-psml

« back to all changes in this revision

Viewing changes to Src/ts_energies.F90

  • Committer: Alberto Garcia
  • Date: 2019-09-02 14:09:43 UTC
  • mfrom: (427.6.323 trunk)
  • Revision ID: albertog@icmab.es-20190902140943-mzmbe1jacgefpgxw
Sync to trunk-776 (notably nc/soc wavefunction support)


Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/ts_energies.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
! This code segment has been fully created by:
 
9
! Nick Papior Andersen, 2019, nickpapior@gmail.com
 
10
! Please conctact the author, prior to re-using this code.
 
11
 
 
12
! Module for calculating energy contributions for the transiesta
 
13
! calculations.
 
14
! Effectively many of the energy contributions stem from Hamiltonian
 
15
! and/or DM elements.
 
16
! However, in TranSiesta we only update/use a subset of the full
 
17
! matrices. As such the energies should only be calculated
 
18
! on the used elements rather than the full.
 
19
!
 
20
! There are many more elements that need to be split, while currently
 
21
! we only calculate a subset of what is actually needed.
 
22
! The main problem in calculating NEGF energies is:
 
23
!  1. the open-boundary problem, and
 
24
!  2. how to handle real-space energies such as Exc from dhscf.
 
25
module ts_energies_m
 
26
 
 
27
  implicit none
 
28
 
 
29
contains
 
30
 
 
31
  subroutine ts_compute_energies()
 
32
 
 
33
    use precision, only: dp
 
34
    use m_ts_options, only: N_Elec, Elecs
 
35
    use m_spin, only: spin
 
36
 
 
37
    use sparse_matrices, only: dit => block_dist
 
38
    use sparse_matrices, only: n_nzs => maxnh
 
39
    use sparse_matrices, only: Dold, Dscf, H
 
40
    use sparse_matrices, only: H_kin_1D, H_vkb_1D
 
41
    use sparse_matrices, only: sp => sparse_pattern
 
42
    
 
43
    use m_ts_method
 
44
    use parallel, only : Node
 
45
#ifdef MPI
 
46
    use mpi_siesta
 
47
#endif
 
48
    use class_OrbitalDistribution
 
49
    use class_Sparsity
 
50
    use class_dSpData1D
 
51
    use class_dSpData2D
 
52
 
 
53
    use geom_helper, only : UCORB
 
54
    use m_ts_electype, only: Elec
 
55
    use m_energies, only: NEGF_Ebs, NEGF_Ekin, NEGF_Etot, NEGF_Enl
 
56
    use m_energies, only: NEGF_DEharr
 
57
 
 
58
! **********************
 
59
! * LOCAL variables    *
 
60
! **********************
 
61
    integer, pointer :: l_ncol(:), l_ptr(:), l_col(:)
 
62
    integer :: no_lo, no_u, lio, io, ind, jo, ir, jr, r, ispin
 
63
    real(dp), pointer :: H_vkb(:), H_kin(:)
 
64
    real(dp) :: Etmp(4,0:1+1+N_Elec*2)
 
65
#ifdef MPI
 
66
    real(dp) :: tmp(4,0:1+1+N_Elec*2)
 
67
    integer :: MPIerror
 
68
#endif
 
69
 
 
70
    if ( spin%NCol .or. spin%SO ) then
 
71
      call die('ts_energies: Not implemented for non-colinear or spin-orbit')
 
72
    end if
 
73
 
 
74
    H_vkb => val(H_vkb_1D)
 
75
    H_kin => val(H_kin_1D)
 
76
 
 
77
    ! Retrieve information about the sparsity pattern
 
78
    call attach(sp, &
 
79
        n_col=l_ncol,list_ptr=l_ptr,list_col=l_col, &
 
80
        nrows=no_lo,nrows_g=no_u)
 
81
    
 
82
    ! Initialize energies
 
83
    Etmp(:,:) = 0._dp
 
84
    
 
85
!$OMP parallel do default(shared), &
 
86
!$OMP&private(lio,io,ir,ind,jo,jr,r,ispin), &
 
87
!$OMP&reduction(+:Etmp)
 
88
    do lio = 1 , no_lo
 
89
 
 
90
      ! obtain the global index of the orbital.
 
91
      io = index_local_to_global(dit,lio,Node)
 
92
      ir = orb_type(io)
 
93
      
 
94
      ! Loop number of entries in the row... (index frame)
 
95
      do ind = l_ptr(lio) + 1 , l_ptr(lio) + l_ncol(lio)
 
96
        
 
97
        ! as the local sparsity pattern is a super-cell pattern,
 
98
        ! we need to check the unit-cell orbital
 
99
        ! The unit-cell column index
 
100
        jo = UCORB(l_col(ind),no_u)
 
101
        jr = orb_type(jo)
 
102
        
 
103
        if      ( all((/ir,jr/) == TYP_BUFFER) ) then
 
104
          r = 1 ! buffer
 
105
        else if ( any((/ir,jr/) == TYP_BUFFER) ) then
 
106
          r = 0 ! other
 
107
        else if ( all((/ir,jr/) == TYP_DEVICE) ) then
 
108
          r = 2 ! device
 
109
        else if ( any((/ir,jr/) == TYP_DEVICE) ) then
 
110
          r = 4+(ir+jr-1)*2 ! device/electrode
 
111
        else if ( ir == jr ) then
 
112
          r = 3+(ir-1)*2 ! electrode/electrode
 
113
        else
 
114
          r = 0 ! other
 
115
        end if
 
116
 
 
117
        ! Ebs
 
118
        Etmp(1,r) = Etmp(1,r) + sum(H(ind,:) * Dscf(ind,:) )
 
119
 
 
120
        do ispin = 1, spin%spinor
 
121
          ! Ekin
 
122
          Etmp(2,r) = Etmp(2,r) + H_kin(ind) * Dscf(ind,ispin)
 
123
          ! Enl
 
124
          Etmp(3,r) = Etmp(3,r) + H_vkb(ind) * Dscf(ind,ispin)
 
125
        end do
 
126
 
 
127
        ! DEharr
 
128
        do ispin = 1, spin%spinor
 
129
          Etmp(4,r) = Etmp(4,r) + H(ind,ispin) * (Dscf(ind,ispin) - Dold(ind,ispin))
 
130
        end do
 
131
 
 
132
      end do
 
133
    end do
 
134
!$OMP end parallel do
 
135
 
 
136
    ! Now add the *other* contributions
 
137
    do ir = 1, N_Elec
 
138
      if ( Elecs(ir)%DM_update >= 1 ) then
 
139
        ! add cross-terms
 
140
        r = 4 + (TYP_DEVICE+ir - 1) * 2
 
141
        Etmp(:,2) = Etmp(:,2) + Etmp(:,r)
 
142
      end if
 
143
      if ( Elecs(ir)%DM_update == 2 ) then
 
144
        r = 3 + (ir - 1) * 2
 
145
        ! add diagonal electrode terms
 
146
        Etmp(:,2) = Etmp(:,2) + Etmp(:,r)
 
147
      end if
 
148
    end do
 
149
 
 
150
#ifdef MPI
 
151
    call MPI_AllReduce(Etmp(1,2),tmp(1,2),size(Etmp, 1), &
 
152
        MPI_Double_Precision,MPI_SUM, MPI_Comm_World,MPIerror)
 
153
    Etmp(:,2) = tmp(:,2)
 
154
#endif
 
155
 
 
156
    ! Copy data over
 
157
    NEGF_Ebs = Etmp(1,2)
 
158
    NEGF_Ekin = Etmp(2,2)
 
159
    NEGF_Enl = Etmp(3,2)
 
160
    NEGF_DEharr = Etmp(4,2)
 
161
 
 
162
  end subroutine ts_compute_energies
 
163
 
 
164
end module ts_energies_m
 
165
 
 
166