plot_coulombmatrix.F90 9.05 KB
Newer Older
Daniel Wortmann's avatar
Daniel Wortmann committed
1 2 3
    !
    ! plot matrix
    !
Matthias Redies's avatar
Matthias Redies committed
4 5 6
    IF (l_plot) THEN
       IF (mpi%isize /= 1) STOP 'coulombmatrix: l_plot only works with one process'
       DO ikpt = 1, kpts%nkpt
Daniel Wortmann's avatar
Daniel Wortmann committed
7 8

          ic = 0
Matthias Redies's avatar
Matthias Redies committed
9 10
          DO i = 1, nbasm1(ikpt)
             DO j = 1, i
Daniel Wortmann's avatar
Daniel Wortmann committed
11
                ic = ic + 1
Matthias Redies's avatar
Matthias Redies committed
12 13 14
                IF (ABS(coulomb(ic, ikpt)) > 1E-8) THEN
                   WRITE (600 + ikpt, '(2i6)') i, j
                   WRITE (600 + ikpt, '(2i6)') j, i
Daniel Wortmann's avatar
Daniel Wortmann committed
15 16 17 18
                END IF
             END DO
          END DO

Matthias Redies's avatar
Matthias Redies committed
19 20
          ALLOCATE (coulhlp(nbasm1(ikpt), nbasm1(ikpt)))
          coulhlp = unpackmat(coulomb(:nbasm1(ikpt)*(nbasm1(ikpt) + 1)/2, ikpt))
Daniel Wortmann's avatar
Daniel Wortmann committed
21 22

          ic = 0
Matthias Redies's avatar
Matthias Redies committed
23 24 25 26 27 28 29
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)
                   DO M = -l, l
                      WRITE (700 + ikpt, *) l, M
                      DO n = 1, hybrid%nindxm1(l, itype)
                         WRITE (700 + ikpt, '(16f8.4)') coulhlp(ic + n, ic + 1:ic + hybrid%nindxm1(l, itype))
Daniel Wortmann's avatar
Daniel Wortmann committed
30
                      END DO
Matthias Redies's avatar
Matthias Redies committed
31
                      ic = ic + hybrid%nindxm1(l, itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
32 33 34 35
                   ENDDO
                END DO
             END DO
          END DO
Matthias Redies's avatar
Matthias Redies committed
36
          ALLOCATE (coulhlp1(nbasm1(ikpt), nbasm1(ikpt)))
Daniel Wortmann's avatar
Daniel Wortmann committed
37 38 39
          coulhlp1 = 0

          ic2 = 0
Matthias Redies's avatar
Matthias Redies committed
40 41 42 43 44
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)
                   DO M = -l, l
                      DO n = 1, hybrid%nindxm1(l, itype) - 1
Daniel Wortmann's avatar
Daniel Wortmann committed
45 46 47 48 49 50 51
                         ic2 = ic2 + 1
                      END DO
                   END DO
                END DO
             END DO
          END DO

Matthias Redies's avatar
Matthias Redies committed
52
          ic = 0
Daniel Wortmann's avatar
Daniel Wortmann committed
53
          ic1 = 0
Matthias Redies's avatar
Matthias Redies committed
54 55 56 57 58 59 60
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)
                   DO M = -l, l
                      DO n = 1, hybrid%nindxm1(l, itype) - 1
                         coulhlp1(ic + n, ic + 1:ic + hybrid%nindxm1(l, itype) - 1) &
                            = coulhlp(ic1 + n, ic1 + 1:ic1 + hybrid%nindxm1(l, itype) - 1)
Daniel Wortmann's avatar
Daniel Wortmann committed
61 62
                      END DO

Matthias Redies's avatar
Matthias Redies committed
63 64
                      coulhlp1(ic2 + 1, ic + 1:ic + hybrid%nindxm1(l, itype) - 1) &
                         = coulhlp(ic1 + hybrid%nindxm1(l, itype), ic1 + 1:ic1 + hybrid%nindxm1(l, itype) - 1)
Daniel Wortmann's avatar
Daniel Wortmann committed
65

Matthias Redies's avatar
Matthias Redies committed
66 67
                      ic = ic + hybrid%nindxm1(l, itype) - 1
                      ic1 = ic1 + hybrid%nindxm1(l, itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
68 69 70 71 72 73 74
                      ic2 = ic2 + 1

                   END DO
                END DO
             END DO
          END DO

Matthias Redies's avatar
Matthias Redies committed
75
          IF (ikpt == 1) THEN
Daniel Wortmann's avatar
Daniel Wortmann committed
76 77 78
             ! add contributions from s channel and G=0 component

             ic = 0; ic1 = 0
Matthias Redies's avatar
Matthias Redies committed
79 80 81 82 83 84 85 86 87 88 89 90
             DO itype = 1, atoms%ntype
                DO ineq = 1, atoms%neq(itype)
                   WRITE (*, *) ic + 1, ic + hybrid%nindxm1(0, itype) - 1, ic1 + 1, &
                      ic1 + hybrid%nindxm1(0, itype) - 1
                   coulhlp1(ic + 1:ic + hybrid%nindxm1(0, itype) - 1, nbasp + 1) &
                      = coulhlp(ic1 + 1:ic1 + hybrid%nindxm1(0, itype) - 1, nbasp + 1)
                   coulhlp1(nbasp + 1, ic + 1:ic + hybrid%nindxm1(0, itype) - 1) &
                      = coulhlp(nbasp + 1, ic1 + 1:ic1 + hybrid%nindxm1(0, itype) - 1)
                   ic = ic + SUM((/((2*l + 1)*(hybrid%nindxm1(l, itype) - 1), &
                                    l=0, hybrid%lcutm1(itype))/))
                   ic1 = ic1 + SUM((/((2*l + 1)*hybrid%nindxm1(l, itype), &
                                      l=0, hybrid%lcutm1(itype))/))
Daniel Wortmann's avatar
Daniel Wortmann committed
91 92 93 94
                END DO
             END DO

             ic2 = 0
Matthias Redies's avatar
Matthias Redies committed
95 96 97 98 99
             DO itype = 1, atoms%ntype
                DO ineq = 1, atoms%neq(itype)
                   DO l = 0, hybrid%lcutm1(itype)
                      DO M = -l, l
                         DO n = 1, hybrid%nindxm1(l, itype) - 1
Daniel Wortmann's avatar
Daniel Wortmann committed
100 101 102 103 104 105 106 107
                            ic2 = ic2 + 1
                         END DO
                      END DO
                   END DO
                END DO
             END DO

             ic1 = 0
Matthias Redies's avatar
Matthias Redies committed
108 109 110 111
             DO itype = 1, atoms%ntype
                DO ineq = 1, atoms%neq(itype)
                   DO l = 0, hybrid%lcutm1(itype)
                      DO M = -l, l
Daniel Wortmann's avatar
Daniel Wortmann committed
112 113
                         ic2 = ic2 + 1

Matthias Redies's avatar
Matthias Redies committed
114
                         ic1 = ic1 + hybrid%nindxm1(l, itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
115

Matthias Redies's avatar
Matthias Redies committed
116 117
                         IF (l /= 0) CYCLE
                         WRITE (900, *) ic2, ic1, itype, ineq
Daniel Wortmann's avatar
Daniel Wortmann committed
118 119 120

                         ic3 = 0
                         ic4 = 0
Matthias Redies's avatar
Matthias Redies committed
121 122 123 124 125 126
                         DO itype1 = 1, atoms%ntype
                            ishift = SUM((/((2*l2 + 1)*hybrid%nindxm1(l2, itype1), &
                                            l2=0, hybrid%lcutm1(itype1))/))
                            ishift1 = SUM((/((2*l2 + 1)*(hybrid%nindxm1(l2, itype1) - 1), &
                                             l2=0, hybrid%lcutm1(itype1))/))
                            DO ineq1 = 1, atoms%neq(itype1)
Daniel Wortmann's avatar
Daniel Wortmann committed
127
                               ic5 = ic3 + (ineq1 - 1)*ishift + 1
Matthias Redies's avatar
Matthias Redies committed
128
                               ic6 = ic5 + hybrid%nindxm1(0, itype1) - 2
Daniel Wortmann's avatar
Daniel Wortmann committed
129 130

                               ic7 = ic4 + (ineq1 - 1)*ishift1 + 1
Matthias Redies's avatar
Matthias Redies committed
131 132
                               ic8 = ic7 + hybrid%nindxm1(0, itype1) - 2
                               WRITE (901, *) ic2, ic7, ic8, ic1, ic5, ic6, itype, itype1
Daniel Wortmann's avatar
Daniel Wortmann committed
133

Matthias Redies's avatar
Matthias Redies committed
134 135 136
                               coulhlp1(ic2, ic7:ic8) = coulhlp(ic1, ic5:ic6)
                               IF (sym%invs) THEN
                                  coulhlp1(ic7:ic8, ic2) = coulhlp1(ic2, ic7:ic8)
Daniel Wortmann's avatar
Daniel Wortmann committed
137
                               ELSE
Matthias Redies's avatar
Matthias Redies committed
138
                                  coulhlp1(ic7:ic8, ic2) = CONJG(coulhlp1(ic2, ic7:ic8))
Daniel Wortmann's avatar
Daniel Wortmann committed
139 140
                               ENDIF
                            END DO
Matthias Redies's avatar
Matthias Redies committed
141
                            ic3 = ic3 + ishift*atoms%neq(itype1)
Daniel Wortmann's avatar
Daniel Wortmann committed
142 143 144 145 146 147 148 149 150 151 152
                            ic4 = ic4 + ishift1*atoms%neq(itype1)
                         END DO

                      END DO
                   END DO
                END DO
             END DO

          END IF

          ic2 = 0
Matthias Redies's avatar
Matthias Redies committed
153 154 155 156 157
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)
                   DO M = -l, l
                      DO n = 1, hybrid%nindxm1(l, itype) - 1
Daniel Wortmann's avatar
Daniel Wortmann committed
158 159 160 161 162 163 164 165
                         ic2 = ic2 + 1
                      END DO
                   END DO
                END DO
             END DO
          END DO

          ic1 = 0
Matthias Redies's avatar
Matthias Redies committed
166 167 168 169
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)
                   DO M = -l, l
Daniel Wortmann's avatar
Daniel Wortmann committed
170 171
                      ic2 = ic2 + 1

Matthias Redies's avatar
Matthias Redies committed
172
                      ic1 = ic1 + hybrid%nindxm1(l, itype)
Daniel Wortmann's avatar
Daniel Wortmann committed
173 174 175

                      ic3 = 0
                      ic4 = ic2
Matthias Redies's avatar
Matthias Redies committed
176 177 178 179 180 181 182 183 184 185
                      DO itype1 = 1, atoms%ntype
                         DO ineq1 = 1, atoms%neq(itype1)
                            DO l1 = 0, hybrid%lcutm1(itype1)
                               DO m1 = -l1, l1
                                  ic3 = ic3 + hybrid%nindxm1(l1, itype1)
                                  IF (ic3 < ic1) CYCLE
                                  WRITE (300, '(4i6,2f15.10)') ic2, ic4, ic1, ic3, coulhlp(ic1, ic3)
                                  coulhlp1(ic2, ic4) = coulhlp(ic1, ic3)
                                  IF (sym%invs) THEN
                                     coulhlp1(ic4, ic2) = coulhlp1(ic2, ic4)
Daniel Wortmann's avatar
Daniel Wortmann committed
186
                                  ELSE
Matthias Redies's avatar
Matthias Redies committed
187
                                     coulhlp1(ic4, ic2) = CONJG(coulhlp1(ic2, ic4))
Daniel Wortmann's avatar
Daniel Wortmann committed
188 189 190 191 192 193 194
                                  ENDIF
                                  ic4 = ic4 + 1
                               END DO
                            END DO
                         END DO
                      END DO

Matthias Redies's avatar
Matthias Redies committed
195 196 197 198
                      DO igpt = 1, hybrid%ngptm(ikpt)
                         coulhlp1(ic2, ic4) = coulhlp(ic1, nbasp + igpt)
                         IF (sym%invs) THEN
                            coulhlp1(ic4, ic2) = coulhlp1(ic2, ic4)
Daniel Wortmann's avatar
Daniel Wortmann committed
199
                         ELSE
Matthias Redies's avatar
Matthias Redies committed
200
                            coulhlp1(ic4, ic2) = CONJG(coulhlp1(ic2, ic4))
Daniel Wortmann's avatar
Daniel Wortmann committed
201 202 203 204 205 206 207 208 209
                         ENDIF
                         ic4 = ic4 + 1
                      END DO

                   END DO
                END DO
             END DO
          END DO

Matthias Redies's avatar
Matthias Redies committed
210
          coulhlp1(nbasp + 1:, nbasp + 1:) = coulhlp(nbasp + 1:, nbasp + 1:)
Daniel Wortmann's avatar
Daniel Wortmann committed
211 212

          ic = 0
Matthias Redies's avatar
Matthias Redies committed
213 214 215 216 217 218 219 220
          DO itype = 1, atoms%ntype
             DO ineq = 1, atoms%neq(itype)
                DO l = 0, hybrid%lcutm1(itype)

                   DO M = -l, l
                      WRITE (800 + ikpt, *) l, M
                      DO n = 1, hybrid%nindxm1(l, itype) - 1
                         WRITE (800 + ikpt, '(16f8.4)') coulhlp1(ic + n, ic + 1:ic + hybrid%nindxm1(l, itype) - 1)
Daniel Wortmann's avatar
Daniel Wortmann committed
221
                      END DO
Matthias Redies's avatar
Matthias Redies committed
222
                      ic = ic + hybrid%nindxm1(l, itype) - 1
Daniel Wortmann's avatar
Daniel Wortmann committed
223 224 225 226 227 228
                   ENDDO
                END DO
             END DO
          END DO

          ic = 0
Matthias Redies's avatar
Matthias Redies committed
229 230 231 232 233
          DO i = 1, nbasm1(ikpt)
             DO j = 1, i
                IF (ABS(coulhlp1(i, j)) > 1E-8) THEN
                   WRITE (850 + ikpt, '(2i6)') i, j
                   WRITE (850 + ikpt, '(2i6)') j, i
Daniel Wortmann's avatar
Daniel Wortmann committed
234 235 236
                END IF
             END DO
          END DO
Matthias Redies's avatar
Matthias Redies committed
237
          DEALLOCATE (coulhlp, coulhlp1)
Daniel Wortmann's avatar
Daniel Wortmann committed
238 239 240
       END DO
       STOP
    END IF ! lplot