← Back to branch summary

~nickpapior/siesta/trunk-kpoint-dos

« back to all changes in this revision

Viewing changes to Src/rhoofd.f

  • Committer: Alberto Garcia
  • Date: 2004-11-25 18:49:43 UTC
  • Revision ID: Arch-1:siesta@uam.es--2004%siesta-devel--reference--0.11--patch-1
Siesta 0.11 -- imported from CVS
Import from cvs using date instead of siesta-0-11-release tag, since
the Pseudo structure was not properly integrated at that time and
did not get the tag.

Show diffs side-by-side

added added

removed removed

Lines of Context:
Src/rhoofd.f
1
 
      subroutine rhoofd( NO, endC,  listC,  C, NSPmax, NSP, 
 
1
      subroutine rhoofd( NO, indxuo, endC,  listC,  C, NSPmax, NSP, 
2
2
     .                   NP, endCt, listCt, CtoCt,
3
3
     .                   NDmax, numDs, listDs, Dscf,
4
 
     .                   RHOscf, MaxDi, Di )
 
4
     .                   RHOscf, MaxAux, ilocal )
5
5
 
6
6
C ********************************************************************
7
7
C Finds the SCF density at the mesh points from the density matrix.
8
8
C Written by P.Ordejon and J.M.Soler. May'95.
 
9
C Re-ordered so that mesh is the outer loop and the orbitals are
 
10
C handled as lower-half triangular. J.D.Gale and J.M.Soler, Feb'99
9
11
C *********************** INPUT **************************************
10
12
C integer NO              : Number of basis orbitals
 
13
C integer indxuo(NO)      : Index of equivalent atom in unit cell
11
14
C integer endC(NO)        : Accumulated umber of nonzero elements in
12
15
C                           each row of C
13
16
C integer listC(*)        : List of nonzero elements in each row of C
26
29
C integer numDs(NO)       : Number of nonzero elemts in each row of Dscf
27
30
C integer listDs(NVmax,NO): List of nonzero elements in each row of Dscf
28
31
C real*8  Dscf(NVmax,NO)  : Value of nonzero elemens in each row of Dscf
29
 
C integer MaxDi           : Dimension of auxiliary array Di
 
32
C integer MaxAux          : Dimension of auxiliary array ilocal
30
33
C *********************** OUTPUT **************************************
31
34
C real*4  RHOscf(NSP,NP)  : SCF density at mesh points
32
35
C ********************* AUXILIARY *************************************
33
 
C real*8  Di(MaxDi)       : Space that can be used freely between calls
34
 
C                           MaxDi must be at least NO
 
36
C integer ilocal(MaxAux)  : Space that can be used freely between calls
 
37
C                           MaxAux must be at least NO
35
38
C *********************************************************************
36
39
 
37
40
      implicit none
38
41
 
39
42
      integer
40
 
     .   NO, NP, NDmax, NSP, NSPmax, MaxDi,
41
 
     .   endC(0:NO),  listC(*),
 
43
     .   NO, NP, NDmax, NSP, NSPmax, MaxAux, 
 
44
     .   endC(0:NO),  listC(*), 
42
45
     .   endCt(0:NP), listCt(*), CtoCt(*),
43
 
     .   numDs(NO), listDs(NDmax,NO)
 
46
     .   numDs(NO), listDs(NDmax,NO), ilocal(MaxAux), indxuo(NO)
44
47
      real
45
48
     .   C(NSPmax,*), RHOscf(NSP,NP)
46
49
      double precision
47
 
     .   Dscf(NDmax,NO), Di(NO)
48
 
 
49
 
      integer i, in, ip, isp, j, jn, kn
50
 
      double precision Ci, Cj, Dij
 
50
     .   Dscf(NDmax,*)
 
51
 
 
52
 
 
53
      integer maxloc
 
54
      parameter ( maxloc = 300 )
 
55
 
 
56
      integer i, ii, il, imp, in, iorb(0:maxloc), ip, isp, iu,
 
57
     .        j, jl, jn, jmp, last
 
58
 
 
59
      double precision Ci, Cj, Dij, Dlocal(0:maxloc,0:maxloc)
51
60
      
52
61
      call timer('rhoofd',1)
53
62
 
54
 
C     Check size of auxiliary array
55
 
      if (NO .gt. MaxDi) stop 'RHOOFD: dimension MaxDi too small'
56
 
 
57
 
C     Initialize RHOscf
58
 
      do 80 ip = 1, NP
59
 
        do 70 isp = 1, NSP
 
63
C  Check size of auxiliary array
 
64
      if (NO .gt. MaxAux) then
 
65
        stop 'RHOOFD: dimension MaxAux too small'
 
66
      endif
 
67
 
 
68
C  Full initialization of arrays ilocal and iorb done only once
 
69
      do j = 1, NO
 
70
        ilocal(j) = 0
 
71
      enddo
 
72
      do il = 0,maxloc
 
73
        iorb(il) = 0
 
74
      enddo
 
75
      last = 0
 
76
 
 
77
C  Loop over grid points
 
78
      do ip = 1,np
 
79
 
 
80
C  Initialise RHOatm
 
81
        do isp = 1, nsp
60
82
          RHOscf(isp,ip) = 0.d0
61
 
   70   continue
62
 
   80 continue
63
 
 
64
 
C     Full initialization of array Di done only once
65
 
      do 90 j = 1, NO
66
 
        Di(j) = 0.d0
67
 
   90 continue
68
 
 
69
 
C     Loop on orbitals
70
 
      do 160 i = 1, NO
71
 
 
72
 
C       Copy row i of Dij from sparse format in Dscf to
73
 
C         full-row format in Di
74
 
        do 100 in = 1, numDs(i)
75
 
          j = listDs(in,i)
76
 
          Di(j) = Dscf(in,i)
77
 
  100   continue
78
 
 
79
 
C       Loop on mesh points of orbital i
80
 
        do 120 in = 1+endC(i-1), endC(i)
81
 
          ip = listC(in)
82
 
 
83
 
C         Loop over non-zero orbitals in mesh point
84
 
          do 110 kn = 1+endCt(ip-1), endCt(ip)
85
 
            j = listCt(kn)
86
 
            jn = CtoCt(kn)
87
 
            Dij = Di(j)
88
 
*           write(6,'(a,i3,5i6,i3,i6)')
89
 
*    .        'rhoofd: i,endC(i-1),in,ip,endCt(ip-1),kn,j,jn=',
90
 
*    .        i, endC(i-1), in, ip, endCt(ip-1), kn, j, jn
91
 
 
92
 
C           Loop over sub-points
93
 
            do 105 isp = 1, NSP
 
83
        enddo
 
84
 
 
85
C  Check that Dlocal can contain Dij for this point.
 
86
        if (endCt(ip)-endCt(ip-1) .gt. maxloc)
 
87
     .    stop 'RHOOFD: parameter maxloc too small'
 
88
 
 
89
C  iorb(il)>0 means that row il of Dlocal must not be overwritten
 
90
C  iorb(il)=0 means that row il of Dlocal is empty
 
91
C  iorb(il)<0 means that row il of Dlocal contains a valid row of 
 
92
C             Dscf, but which is not required at this point
 
93
        do imp = 1+endCt(ip-1), endCt(ip)
 
94
          i = listCt(imp)
 
95
          il = ilocal(i)
 
96
          if (il.gt.0) iorb(il) = i
 
97
        enddo
 
98
 
 
99
C  Look for required rows of Dscf not yet stored in Dlocal
 
100
        do imp = 1+endCt(ip-1), endCt(ip)
 
101
          i = listCt(imp)
 
102
          if (ilocal(i) .eq. 0) then
 
103
 
 
104
C           Look for an available row in Dlocal
 
105
            do il = 1,maxloc
 
106
C             last runs circularly over rows of Dlocal
 
107
              last = last + 1
 
108
              if (last .gt. maxloc) last = 1
 
109
              if (iorb(last) .le. 0) goto 10
 
110
            enddo
 
111
              stop 'rhoofd: no slot available in Dlocal'
 
112
   10       continue
 
113
 
 
114
C  Copy row i of Dscf into row last of Dlocal
 
115
            j = abs(iorb(last))
 
116
            if (j.ne.0) ilocal(j) = 0
 
117
            ilocal(i) = last
 
118
            iorb(last) = i
 
119
            il = last
 
120
            iu = indxuo(i)
 
121
            do ii = 1, numDs(i)
 
122
              j = listDs(ii,i)
 
123
              jl = ilocal(j)
 
124
              if (i .eq. j) then
 
125
                Dij = Dscf(ii,iu)
 
126
              else
 
127
                Dij = 2*Dscf(ii,iu)
 
128
              endif
 
129
              Dlocal(il,jl) = Dij
 
130
              Dlocal(jl,il) = Dij
 
131
            enddo
 
132
          endif
 
133
        enddo
 
134
 
 
135
C  Loop on first orbital of mesh point
 
136
        do imp = 1+endCt(ip-1), endCt(ip)
 
137
          i = listCt(imp)
 
138
          il = ilocal(i)
 
139
          iu = indxuo(i)
 
140
          in = CtoCt(imp)
 
141
 
 
142
C  Loop on second orbital of mesh point
 
143
          do jmp = 1+endCt(ip-1), imp
 
144
            j = listCt(jmp)
 
145
            jl = ilocal(j)
 
146
            jn = CtoCt(jmp)
 
147
            Dij = Dlocal(il,jl)
 
148
 
 
149
C  Loop over sub-points
 
150
            do isp = 1, nsp
94
151
              Ci = C(isp,in)
95
152
              Cj = C(isp,jn)
96
153
              RHOscf(isp,ip) = RHOscf(isp,ip) + Dij * Ci * Cj
97
 
  105       continue
98
 
 
99
 
  110     enddo
100
 
 
101
 
  120   enddo
102
 
 
103
 
C       Restore Di(j) for next orbital i
104
 
        do 150 in = 1, numDs(i)
105
 
          j = listDs(in,i)
106
 
          Di(j) = 0.d0
107
 
  150   continue
108
 
 
109
 
  160 enddo
 
154
            enddo
 
155
 
 
156
          enddo
 
157
        enddo
 
158
 
 
159
C  Restore iorb for next point
 
160
        do imp = 1+endCt(ip-1), endCt(ip)
 
161
          i = listCt(imp)
 
162
          il = ilocal(i)
 
163
          iorb(il) = -i
 
164
        enddo
 
165
 
 
166
      enddo
110
167
 
111
168
      call timer('rhoofd',2)
112
169
      end