← Back to branch summary

~dikan.va/siesta/thtr-func-self-test-example-Jks

« back to all changes in this revision

Viewing changes to Src/rmatrix.f

  • Committer: Alberto Garcia
  • Date: 2019-05-09 14:35:26 UTC
  • Revision ID: albertog@icmab.es-20190509143526-patkkeluhx730ltp
First installment of computation of R_mu_nu elements for thermal transport

(Work in progress)

routine rmatrix in module m_rmatrix in rmatrix.f is modelled after 'overlap'.

It computes Rmat, wich is a pointer to the value section of the dSp2D object Rmat_2D.
Both live in module sparse_matrices.

TODO:

-- Use the correct origin for the X,Y,Z matrix elements (matel puts the origin in
the second orbital).
-- Debug obscure error coming from the spher_harm module.
-- Debug parallel operation
-- "Deletion" of Rmat_2D in reset_sparse_matrices might need to be wrapped.


Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/rmatrix.f
 
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_rmatrix
 
9
 
 
10
      use precision,     only : dp
 
11
      use parallel,      only : Node, Nodes
 
12
      use parallelsubs,  only : GlobalToLocalOrb
 
13
      use atmfuncs,      only : rcut, orb_gindex
 
14
      use neighbour,     only : jna=>jan, r2ij, xij, mneighb,
 
15
     &                          reset_neighbour_arrays
 
16
      use alloc,         only : re_alloc, de_alloc
 
17
      use m_new_matel,   only : new_matel
 
18
      use m_iodm_old,    only : write_dm
 
19
      use m_matio,       only : write_mat
 
20
      use atomlist, only: no_l
 
21
      use fdf
 
22
 
 
23
      implicit none
 
24
 
 
25
      public :: rmatrix
 
26
      private
 
27
 
 
28
      CONTAINS
 
29
 
 
30
      subroutine rmatrix(nua, na, no, scell, xa, indxua, rmaxo,
 
31
     &                   maxnh, lasto, iphorb, isa, 
 
32
     &                   numh, listhptr, listh, Rmat)
 
33
C *********************************************************************
 
34
C Computes the "R" matrix < phi_1 | R | phi_2>
 
35
C Energies in Ry. Lengths in Bohr.
 
36
C Based on 'overlap'
 
37
C **************************** INPUT **********************************
 
38
C integer nua              : Number of atoms in unit cell
 
39
C integer na               : Number of atoms in supercell
 
40
C integer no               : Number of orbitals in supercell
 
41
C real*8  scell(3,3)       : Supercell vectors SCELL(IXYZ,IVECT)
 
42
C real*8  xa(3,na)         : Atomic positions in cartesian coordinates
 
43
c integer indxua(na)       : Index of equivalent atom in unit cell
 
44
C real*8  rmaxo            : Maximum cutoff for atomic orbitals
 
45
C integer maxnh            : First dimension of S and listh
 
46
C integer lasto(0:na)      : Last orbital index of each atom
 
47
C integer iphorb(no)       : Orbital index of each orbital in its atom
 
48
C integer isa(na)          : Species index of each atom
 
49
C integer numh(nuotot)     : Number of nonzero elements of each row
 
50
C                            of the overlap matrix
 
51
C integer listhptr(nuotot) : Pointer to start of rows (-1) of overlap
 
52
C                            matrix
 
53
C integer listh(maxnh)     : Column indexes of the nonzero elements  
 
54
C                            of each row of the overlap matrix
 
55
C **************************** OUTPUT *********************************
 
56
C real*8  Rmat(maxnh,3)    : Matrix elements of "R"
 
57
 
 
58
      integer, intent(in)   ::  maxnh, na, no, nua
 
59
      integer, intent(in)   :: indxua(na), iphorb(no), isa(na),
 
60
     &                         lasto(0:na), listh(maxnh), numh(*),
 
61
     &                         listhptr(*)
 
62
      real(dp) , intent(in) :: scell(3,3), rmaxo, xa(3,na)
 
63
      real(dp), intent(out) :: Rmat(maxnh,3)
 
64
C Internal variables ......................................................
 
65
      integer               :: ia, ind, io, ioa, is,  iio, j, ja, jn,
 
66
     &                         jo, joa, js, jua, nnia, ig, jg
 
67
      real(dp)              :: grSij(3) , rij, Sij
 
68
      real(dp),     pointer :: Ri(:,:)
 
69
      external  timer
 
70
 
 
71
C     Start timer
 
72
      call timer( 'rmatrix', 1 )
 
73
 
 
74
C     Initialize neighb subroutine 
 
75
      call mneighb( scell, 2.d0*rmaxo, na, xa, 0, 0, nnia )
 
76
 
 
77
C     Allocate local memory
 
78
      nullify(Ri)
 
79
      call re_alloc( Ri, 1, no, 1, 3, 'Ri', 'rmatrix' )
 
80
 
 
81
      do ia = 1,nua
 
82
        is = isa(ia)
 
83
        call mneighb( scell, 2.d0*rmaxo, na, xa, ia, 0, nnia )
 
84
        do io = lasto(ia-1)+1,lasto(ia)
 
85
 
 
86
C         Is this orbital on this Node?
 
87
          call GlobalToLocalOrb(io,Node,Nodes,iio)
 
88
          if (iio.gt.0) then
 
89
 
 
90
C           Valid orbital
 
91
            ioa = iphorb(io)
 
92
            ig = orb_gindex(is,ioa)
 
93
            do jn = 1,nnia
 
94
              ja = jna(jn)
 
95
              jua = indxua(ja)
 
96
              rij = sqrt( r2ij(jn) )
 
97
              do jo = lasto(ja-1)+1,lasto(ja)
 
98
                joa = iphorb(jo)
 
99
                js = isa(ja)
 
100
                !
 
101
                ! Use global indexes for new version of matel
 
102
                !
 
103
                jg = orb_gindex(js,joa)
 
104
                if (rcut(is,ioa)+rcut(js,joa) .gt. rij) then
 
105
                  call new_MATEL( 'X', ig, jg, xij(1:3,jn),
 
106
     &                        Sij, grSij )
 
107
                  Ri(1,jo) = Ri(1,jo) + Sij
 
108
                  call new_MATEL( 'Y', ig, jg, xij(1:3,jn),
 
109
     &                        Sij, grSij )
 
110
                  Ri(2,jo) = Ri(2,jo) + Sij
 
111
                  call new_MATEL( 'Z', ig, jg, xij(1:3,jn),
 
112
     &                        Sij, grSij )
 
113
                  Ri(3,jo) = Ri(3,jo) + Sij
 
114
                endif
 
115
              enddo
 
116
            enddo
 
117
            do j = 1,numh(iio)
 
118
              ind = listhptr(iio)+j
 
119
              jo = listh(ind)
 
120
              Rmat(ind,:) = Ri(jo,:)
 
121
              Ri(jo,:) = 0.0d0
 
122
            enddo
 
123
          endif
 
124
        enddo
 
125
      enddo
 
126
 
 
127
C     Deallocate local memory
 
128
      call reset_neighbour_arrays( )
 
129
      call de_alloc( Ri, 'Ri', 'rmatrix' )
 
130
 
 
131
C     Finish timer
 
132
      call timer( 'rmatrix', 2 )
 
133
      end subroutine rmatrix
 
134
      end module m_rmatrix