Merge lp:~nickpapior/siesta/4.1-wfs into lp:~albertog/siesta/4.1-wfs
- 4.1-wfs
- Merge into 4.1-wfs
Status: | Merged |
---|---|
Merged at revision: | 1116 |
Proposed branch: | lp:~nickpapior/siesta/4.1-wfs |
Merge into: | lp:~albertog/siesta/4.1-wfs |
Diff against target: |
813 lines (+196/-188) 6 files modified
Src/Makefile (+1/-1) Src/bands.F (+25/-29) Src/print_spin.F90 (+51/-58) Src/siesta_analysis.F (+96/-73) Src/writewave.F (+22/-26) version.info (+1/-1) |
To merge this branch: | bzr merge lp:~nickpapior/siesta/4.1-wfs |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Alberto Garcia | Approve | ||
Review via email: mp+371082@code.launchpad.net |
Commit message
Implemented tSpin in bands/writewave and fixed a few output options
Using tSpin cleans up a bit some of the *magic* variables used
here and there. Actually the use of tSpin suggests that we need
a new variable in tSpin to clarify whether it is NCol .or. SO
which could help in certain places.
The output files are now a bit weird when comparing format
for NC.
- ioeig uses no_u * 2, 1
- bands uses no_u, 4 for WFS file
- bands uses no_u * 2, 1 for .bands
- wwave uses no_u, 4 for WFSX
It is not clear to me why they are different, and/or
if they *should* be different.
I would prefer if we can use the same convention all-over.
This may require changes all-around.
I don't have any preference. Probably I would prefer:
no_u, spin%H
which clarifies exactly which type of calculation it was.
This is probably opposite of my previous perception.
Also, I have cleaned up print_spin which did lots of
duplicate work for SO calculations.
Description of the change
Once this is in, I can approve the merge.
I think it would be best if we could clean up this once and for all.
Let me know what you think!
Alberto Garcia (albertog) wrote : | # |
More on the bands/EIG/WFS headings:
I agree that they could be improved, but doing so at this time would delay the pressing merge of the NC/SO wfs functionality. Besides, just changing a small set of numbers by another small set of numbers is probably not going to do a lot towards self-documentation of the files. To get out of the historical constraints we probably need a more thorough rethink of the files' contents.
Nick Papior (nickpapior) wrote : | # |
For EIG files etc. Agreed, it isn't totally necessary. I just wanted this on the table since it is rather confusing ;)
Once you know what they mean, it is easy to bypass the problems they yield.
Preview Diff
1 | === modified file 'Src/Makefile' |
2 | --- Src/Makefile 2019-06-20 13:50:42 +0000 |
3 | +++ Src/Makefile 2019-08-08 11:56:55 +0000 |
4 | @@ -1346,7 +1346,7 @@ |
5 | write_subs.o: m_steps.o m_stress.o m_ts_global_vars.o parallel.o precision.o |
6 | write_subs.o: siesta_cml.o siesta_geom.o siesta_options.o units.o zmatrix.o |
7 | writewave.o: alloc.o atmfuncs.o atomlist.o densematrix.o diag.o diag_option.o |
8 | -writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o parallel.o |
9 | +writewave.o: get_kpoints_scale.o kpoint_grid.o m_spin.o m_spin.o parallel.o |
10 | writewave.o: parallelsubs.o precision.o siesta_geom.o sys.o units.o |
11 | xml.o: precision.o |
12 | zm_broyden_optim.o: m_broyddj_nocomm.o parallel.o precision.o sys.o units.o |
13 | |
14 | === modified file 'Src/bands.F' |
15 | --- Src/bands.F 2019-02-24 16:24:33 +0000 |
16 | +++ Src/bands.F 2019-08-08 11:56:55 +0000 |
17 | @@ -296,7 +296,7 @@ |
18 | end subroutine initbands |
19 | |
20 | |
21 | - subroutine bands( no_s, h_spin_dim, spinor_dim, |
22 | + subroutine bands( no_s, spin, |
23 | . no_u, no_l,maxnh, maxk, |
24 | . numh, listhptr, listh, H, S, ef, xij, indxuo, |
25 | . writeb, nk, kpoint, ek, occtol, getPSI ) |
26 | @@ -306,11 +306,11 @@ |
27 | C Written by J.Soler, August 1997 and August 1998. |
28 | C Initialisation moved into a separate routine, JDG Jan 2000. |
29 | C WFS options by A. Garcia, April 2012 |
30 | +C Cleaned for tSpin by N. Papior, August 2019 |
31 | |
32 | C **************************** INPUT ********************************** |
33 | C integer no_s : Number of basis orbitals in supercell |
34 | -C integer h_spin_dim : |
35 | -C integer spinor_dim : |
36 | +C type(tSpin) spin : Containing all spin-information |
37 | C integer maxnh : Maximum number of orbitals interacting |
38 | C with any orbital |
39 | C integer maxk : Last dimension of kpoint and ek |
40 | @@ -320,7 +320,7 @@ |
41 | C hamiltonian matrix |
42 | C integer listh(maxlh) : Nonzero hamiltonian-matrix element |
43 | C column indexes for each matrix row |
44 | -C real*8 H(maxnh,h_spin_dim) : Hamiltonian in sparse form |
45 | +C real*8 H(maxnh,spin%H) : Hamiltonian in sparse form |
46 | C real*8 S(maxnh) : Overlap in sparse form |
47 | C real*8 ef : Fermi energy |
48 | C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) |
49 | @@ -338,7 +338,7 @@ |
50 | C real*8 kpoint(3,maxk) : k point vectors |
51 | C real*8 occtol : Occupancy threshold for DM build |
52 | C *************************** OUTPUT ********************************** |
53 | -C real*8 ek(no_u,spinor_dim,maxk) : Eigenvalues |
54 | +C real*8 ek(no_u,spin%spinor,maxk) : Eigenvalues |
55 | C *************************** UNITS *********************************** |
56 | C Lengths in atomic units (Bohr). |
57 | C k vectors in reciprocal atomic units. |
58 | @@ -366,7 +366,7 @@ |
59 | use atmfuncs, only : symfio, cnfigfio, labelfis, nofis |
60 | use writewave, only : wfs_filename |
61 | |
62 | - use m_spin, only: NoMagn, SPpol, NonCol, SpOrb |
63 | + use t_spin, only: tSpin |
64 | |
65 | use m_diag_option, only: ParallelOverK, Serial |
66 | use m_diag, only: diag_init |
67 | @@ -374,13 +374,13 @@ |
68 | |
69 | implicit none |
70 | |
71 | - integer :: h_spin_dim |
72 | - integer :: maxk, maxnh, spinor_dim, no_u, no_l, nk, no_s, |
73 | + type(tSpin), intent(in) :: spin |
74 | + integer :: maxk, maxnh, no_u, no_l, nk, no_s, |
75 | . indxuo(no_s), listh(maxnh), numh(no_l), |
76 | . listhptr(*) |
77 | logical :: writeb |
78 | - real(dp) :: ef, ek(no_u,spinor_dim,maxk), |
79 | - . H(maxnh,h_spin_dim), kpoint(3,maxk), |
80 | + real(dp) :: ef, ek(no_u,spin%spinor,maxk), |
81 | + . H(maxnh,spin%H), kpoint(3,maxk), |
82 | . S(maxnh), xij(3,maxnh), occtol |
83 | logical, intent(in) :: getPSI |
84 | |
85 | @@ -394,8 +394,7 @@ |
86 | logical, parameter :: fixspin = .false. |
87 | |
88 | integer :: ik, il, io, ispin, iu, iu_wfs, iuo, naux, nhs, j |
89 | - integer :: nspin_flag ! For wfs output |
90 | - logical :: SaveParallelOverK, non_coll |
91 | + logical :: SaveParallelOverK |
92 | |
93 | real(dp) |
94 | . Dnew, qs(2), e1, e2, efs(2), emax, emin, Enew, eV, qk, qtot, |
95 | @@ -424,15 +423,12 @@ |
96 | |
97 | C Allocate local arrays - only aux is relevant here |
98 | |
99 | - non_coll = (h_spin_dim >= 4) |
100 | - if ( non_coll ) then |
101 | - nspin_flag = 4 |
102 | - nhs = 2 * (2*no_u) * (2*no_l) |
103 | + if ( spin%Grid == 4 ) then |
104 | + nhs = 2 * (2*no_u) * (2*no_l) |
105 | else |
106 | - nspin_flag = h_spin_dim |
107 | nhs = 2 * no_u*no_l |
108 | endif |
109 | - naux = 2*no_u*5 |
110 | + naux = 2*no_u*5 |
111 | call allocDenseMatrix(nhs, nhs, nhs) |
112 | call re_alloc( aux, 1, naux, 'aux', 'bands' ) |
113 | |
114 | @@ -445,7 +441,7 @@ |
115 | rewind (iu_wfs) |
116 | |
117 | write(iu_wfs) nk, .false. ! nk, Gamma, same file-format in WFS as for Gamma-point |
118 | - write(iu_wfs) nspin_flag |
119 | + write(iu_wfs) spin%Grid |
120 | write(iu_wfs) no_u |
121 | write(iu_wfs) (iaorb(j),labelfis(isa(iaorb(j))), |
122 | . iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)), |
123 | @@ -455,7 +451,7 @@ |
124 | endif |
125 | |
126 | C Find the band energies |
127 | - if (NoMagn .or. SPpol) then |
128 | + if ( spin%Grid < 4 ) then |
129 | C fixspin and qs are not used in diagk, since getD=.false. ... |
130 | qs(1) = 0.0_dp |
131 | qs(2) = 0.0_dp |
132 | @@ -472,7 +468,7 @@ |
133 | call diag_init() |
134 | end if |
135 | |
136 | - call diagk( h_spin_dim, no_l, no_s, spinor_dim, |
137 | + call diagk( spin%H, no_l, no_s, spin%spinor, |
138 | . maxnh, maxnh, |
139 | . no_u, numh, listhptr, listh, numh, listhptr, |
140 | . listh, H, S, getD, getPSI, fixspin, qtot, qs, temp, |
141 | @@ -484,7 +480,7 @@ |
142 | ParallelOverK = SaveParallelOverK |
143 | if ( ParallelOverK ) Serial = .true. |
144 | |
145 | - elseif ( NonCol ) then |
146 | + elseif ( spin%NCol ) then |
147 | call diag2k(no_l, no_s, maxnh, maxnh, no_u, |
148 | . numh, listhptr, listh, numh, listhptr, |
149 | . listh, H, S, getD, getPSI, qtot, temp, e1, e2, |
150 | @@ -493,7 +489,7 @@ |
151 | . Haux, Saux, psi, Haux, Saux, aux, |
152 | . no_u, occtol, 1, no_u ) |
153 | |
154 | - elseif ( SpOrb ) then |
155 | + elseif ( spin%SO ) then |
156 | call diag3k(no_l, no_s, maxnh, maxnh, no_u, numh, |
157 | . listhptr, listh, numh, listhptr, listh, |
158 | . H, S, getD, getPSI, qtot, temp, e1, e2, xij, |
159 | @@ -502,7 +498,7 @@ |
160 | . Haux, Saux, psi, Haux, Saux, aux, |
161 | . no_u, occtol, 1, no_u ) |
162 | else |
163 | - call die( 'bands: ERROR: incorrect value of h_spin_dim') |
164 | + call die( 'bands: ERROR: incorrect value of spin%H') |
165 | endif |
166 | |
167 | C Write bands |
168 | @@ -525,7 +521,7 @@ |
169 | . path = path + sqrt( (kpoint(1,ik)-kpoint(1,ik-1))**2 + |
170 | . (kpoint(2,ik)-kpoint(2,ik-1))**2 + |
171 | . (kpoint(3,ik)-kpoint(3,ik-1))**2 ) |
172 | - do ispin = 1,spinor_dim |
173 | + do ispin = 1, spin%spinor |
174 | do io = 1, no_u |
175 | emax = max( emax, ek(io,ispin,ik) ) |
176 | emin = min( emin, ek(io,ispin,ik) ) |
177 | @@ -537,10 +533,10 @@ |
178 | write(iu,*) emin/eV, emax/eV |
179 | |
180 | C Write eigenvalues |
181 | - if ( non_coll) then |
182 | + if ( spin%Grid == 4 ) then |
183 | write(iu,*) 2*no_u, 1, nk ! A single spin channel, double number of bands |
184 | else |
185 | - write(iu,*) no_u, spinor_dim, nk |
186 | + write(iu,*) no_u, spin%spinor, nk |
187 | endif |
188 | |
189 | ! For non-collinear calculations, the ek array will contain the first no_u |
190 | @@ -554,11 +550,11 @@ |
191 | . (kpoint(3,ik)-kpoint(3,ik-1))**2 ) |
192 | write(iu,'(f10.6,10f12.4,/,(10x,10f12.4))') |
193 | . path, ((ek(io,ispin,ik)/eV,io=1,no_u), |
194 | - . ispin=1,spinor_dim) |
195 | + . ispin=1,spin%spinor) |
196 | else |
197 | write(iu,'(3f9.5,8f12.4,/,(27x,8f12.4))') |
198 | . kpoint(1,ik), kpoint(2,ik), kpoint(3,ik), |
199 | - . ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,spinor_dim) |
200 | + . ((ek(io,ispin,ik)/eV,io=1,no_u),ispin=1,spin%spinor) |
201 | endif |
202 | enddo |
203 | |
204 | |
205 | === modified file 'Src/print_spin.F90' |
206 | --- Src/print_spin.F90 2019-03-29 13:10:59 +0000 |
207 | +++ Src/print_spin.F90 2019-08-08 11:56:55 +0000 |
208 | @@ -2,7 +2,6 @@ |
209 | ! |
210 | ! Prints spin in output and CML files |
211 | ! |
212 | - use m_spin, only: nspin ! This is spin%grid (1, 2, or 4) |
213 | use m_spin, only: spin |
214 | use atomlist, only: no_l |
215 | use sparse_matrices, only: listhptr, numh |
216 | @@ -16,79 +15,73 @@ |
217 | |
218 | implicit none |
219 | |
220 | - real(dp), intent(out) :: qspin(4) |
221 | + real(dp), intent(out) :: qspin(spin%Grid) |
222 | |
223 | integer :: ispin, io, j, ind |
224 | real(dp) :: qaux |
225 | real(dp) :: Stot ! Total spin magnitude |
226 | real(dp) :: Svec(3) ! Total spin vector |
227 | - real(dp) :: dm_aux(4) ! Auxiliary array to hermitify Dscf for SOC |
228 | #ifdef MPI |
229 | - real(dp) :: qtmp(4) |
230 | + real(dp) :: qtmp(spin%Grid) |
231 | #endif |
232 | |
233 | - if (nspin < 2) RETURN |
234 | - |
235 | - if (spin%SO) then |
236 | + qspin(:) = 0.0_dp |
237 | + |
238 | + if ( spin%Grid < 2 ) return |
239 | + |
240 | + if ( spin%SO ) then |
241 | |
242 | - do ispin = 1,nspin |
243 | - qspin(ispin) = 0.0_dp |
244 | - do io = 1,no_l |
245 | - do j = 1,numh(io) |
246 | - ind = listhptr(io)+j |
247 | - ! In the SOC case, hermitify Dscf |
248 | - dm_aux(1:4) = dscf(ind,1:4) |
249 | - dm_aux(3) = 0.5_dp*(dm_aux(3)+dscf(ind,7)) |
250 | - dm_aux(4) = 0.5_dp*(dm_aux(4)+dscf(ind,8)) |
251 | - qspin(ispin) = qspin(ispin) + dm_aux(ispin) * S(ind) |
252 | - enddo |
253 | - enddo |
254 | - enddo |
255 | + do io = 1,no_l |
256 | + do j = 1,numh(io) |
257 | + ind = listhptr(io)+j |
258 | + ! In the SOC case, hermitify Dscf |
259 | + qspin(1:2) = qspin(1:2) + dscf(ind,1:2) * S(ind) |
260 | + qspin(3) = qspin(3) + 0.5_dp*(dscf(ind,3)+dscf(ind,7)) * S(ind) |
261 | + qspin(4) = qspin(4) + 0.5_dp*(dscf(ind,4)+dscf(ind,8)) * S(ind) |
262 | + end do |
263 | + end do |
264 | |
265 | else |
266 | |
267 | - do ispin = 1,nspin |
268 | - qspin(ispin) = 0.0_dp |
269 | - do io = 1,no_l |
270 | - do j = 1,numh(io) |
271 | - ind = listhptr(io)+j |
272 | - qspin(ispin) = qspin(ispin) + Dscf(ind,ispin) * S(ind) |
273 | - enddo |
274 | - enddo |
275 | - enddo |
276 | + do io = 1,no_l |
277 | + do j = 1,numh(io) |
278 | + ind = listhptr(io)+j |
279 | + qspin(:) = qspin(:) + Dscf(ind,:) * S(ind) |
280 | + end do |
281 | + end do |
282 | |
283 | endif |
284 | #ifdef MPI |
285 | - ! Global reduction of spin components |
286 | - call globalize_sum(qspin(1:nspin),qtmp(1:nspin)) |
287 | - qspin(1:nspin) = qtmp(1:nspin) |
288 | + ! Global reduction of spin components |
289 | + call globalize_sum(qspin(1:spin%Grid),qtmp(1:spin%Grid)) |
290 | + qspin(1:spin%Grid) = qtmp(1:spin%Grid) |
291 | #endif |
292 | |
293 | - if (nspin .eq. 2) then |
294 | - |
295 | - if (IOnode) then |
296 | - Svec(1) = 0.0_dp |
297 | - Svec(2) = 0.0_dp |
298 | - Svec(3) = qspin(1) - qspin(2) |
299 | - Stot = Svec(3) |
300 | - write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec |
301 | - endif |
302 | - if (cml_p) call cmlAddProperty(xf=mainXML, & |
303 | - value=qspin(1)-qspin(2), dictref='siesta:stot', & |
304 | - units='siestaUnits:spin') |
305 | - |
306 | - elseif (nspin .eq. 4) then |
307 | - |
308 | - call spnvec( nspin, qspin, qaux, Stot, Svec ) |
309 | - if (IONode) then |
310 | - write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec |
311 | - if (cml_p) then |
312 | - call cmlAddProperty(xf=mainXML, value=Stot, & |
313 | - dictref='siesta:stot', units='siestaUnits:spin') |
314 | - call cmlAddProperty(xf=mainXML, value=Svec, & |
315 | - dictref='siesta:svec', units='siestaUnits:spin') |
316 | - endif !cml_p |
317 | - endif |
318 | - endif |
319 | + ! We are left with printing out to stdout |
320 | + if ( .not. IONode ) return |
321 | + |
322 | + if ( spin%Grid == 2) then |
323 | + |
324 | + Svec(1) = 0.0_dp |
325 | + Svec(2) = 0.0_dp |
326 | + Svec(3) = qspin(1) - qspin(2) |
327 | + Stot = Svec(3) |
328 | + write(6,'(5x,a,f10.5,2f10.1,f10.5)') 'spin moment: S , {S} = ', Stot, Svec |
329 | + if (cml_p) call cmlAddProperty(xf=mainXML, & |
330 | + value=qspin(1)-qspin(2), dictref='siesta:stot', & |
331 | + units='siestaUnits:spin') |
332 | + |
333 | + else if ( spin%Grid == 4 ) then |
334 | + |
335 | + call spnvec( spin%Grid, qspin, qaux, Stot, Svec ) |
336 | + write(6,'(5x,a,4f10.5)') 'spin moment: S , {S} = ', Stot, Svec |
337 | + if (cml_p) then |
338 | + call cmlAddProperty(xf=mainXML, value=Stot, & |
339 | + dictref='siesta:stot', units='siestaUnits:spin') |
340 | + call cmlAddProperty(xf=mainXML, value=Svec, & |
341 | + dictref='siesta:svec', units='siestaUnits:spin') |
342 | + end if !cml_p |
343 | + |
344 | + end if |
345 | |
346 | end subroutine print_spin |
347 | |
348 | === modified file 'Src/siesta_analysis.F' |
349 | --- Src/siesta_analysis.F 2019-02-27 08:38:58 +0000 |
350 | +++ Src/siesta_analysis.F 2019-08-08 11:56:55 +0000 |
351 | @@ -46,8 +46,7 @@ |
352 | use m_energies |
353 | use m_steps, only: final |
354 | use m_ntm |
355 | - use m_spin, only: nspin, spinor_dim, h_spin_dim |
356 | - use m_spin, only: SpOrb, NonCol, SPpol, NoMagn |
357 | + use m_spin, only: spin |
358 | use m_dipol |
359 | use m_eo |
360 | use m_forces, only: fa |
361 | @@ -68,10 +67,10 @@ |
362 | |
363 | real(dp) :: dummy_str(3,3) |
364 | real(dp) :: dummy_strl(3,3) |
365 | - real(dp) :: qspin(4) ! Local |
366 | + real(dp) :: qspin(spin%Grid) ! Local |
367 | |
368 | - real(dp) :: polxyz(3,nspin) ! Autom., small |
369 | - real(dp) :: polR(3,nspin) ! Autom., small |
370 | + real(dp) :: polxyz(3,spin%Grid) ! Autom., small |
371 | + real(dp) :: polR(3,spin%Grid) ! Autom., small |
372 | real(dp) :: qaux |
373 | real(dp), pointer :: ebk(:,:,:) ! Band energies |
374 | |
375 | @@ -127,7 +126,7 @@ |
376 | |
377 | if (fdf_get("Read-H-from-file",.false.)) then |
378 | if (SIESTA_worker) then |
379 | - call read_spmatrix(maxnh, no_l, h_spin_dim, numh, |
380 | + call read_spmatrix(maxnh, no_l, spin%H, numh, |
381 | . listhptr, listh, H, found, userfile="H_IN") |
382 | if (.not. found) call die("Could not find H_IN") |
383 | current_ef = ef |
384 | @@ -137,7 +136,7 @@ |
385 | |
386 | #ifdef SIESTA__PEXSI |
387 | if (fdf_get("PEXSI.DOS",.false.)) then |
388 | - call pexsi_dos(no_u, no_l, spinor_dim, |
389 | + call pexsi_dos(no_u, no_l, spin%spinor, |
390 | $ maxnh, numh, listhptr, listh, H, S, qtot, ef) |
391 | |
392 | endif |
393 | @@ -263,6 +262,17 @@ |
394 | ! latest DM, although this is probably too much). |
395 | ! See the logic in siesta_forces. |
396 | |
397 | +! Setup the number of states depending on the used method |
398 | +! For polarized calculations we use no_u and separate spin-channels. |
399 | +! For NC/SO we have mixed spin-channels and separating makes no |
400 | +! sense. |
401 | + if ( spin%Grid == 4 ) then |
402 | +! We are NC/SO |
403 | + max_n_states = no_u * 2 |
404 | + else |
405 | + max_n_states = no_u |
406 | + end if |
407 | + |
408 | ! Find and print wavefunctions at selected k-points |
409 | !** This uses H,S, and xa |
410 | if (nwk.gt.0) then |
411 | @@ -270,7 +280,7 @@ |
412 | if (IONode) print "(a)", |
413 | $ "Writing WFSX for selected k-points in " |
414 | $ // trim(wfs_filename) |
415 | - call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, |
416 | + call wwave( no_s, spin, no_u, no_l, maxnh, |
417 | & nwk, |
418 | & numh, listhptr, listh, H, S, Ef, xijo, indxuo, |
419 | & gamma_wavefunctions, nwk, wfk, no_u, occtol ) |
420 | @@ -278,13 +288,10 @@ |
421 | |
422 | |
423 | !** This uses H,S, and xa |
424 | - if (write_coop) then |
425 | + if ( write_coop ) then |
426 | ! Output the wavefunctions for the kpoints in the SCF set |
427 | ! Note that we allow both a band number and an energy range |
428 | ! The user is responsible for using appropriate values. |
429 | - max_n_states = no_u |
430 | - if (h_spin_dim >= 4) max_n_states = 2*no_u |
431 | - |
432 | wfs_band_min = fdf_get("WFS.BandMin",1) |
433 | wfs_band_max = fdf_get("WFS.BandMax",max_n_states) |
434 | call setup_wfs_list(nkpnt,max_n_states, |
435 | @@ -294,7 +301,7 @@ |
436 | wfs_filename = trim(slabel)//".fullBZ.WFSX" |
437 | if (IONode) print "(a)", "Writing WFSX for COOP/COHP in " |
438 | $ // trim(wfs_filename) |
439 | - call wwave( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, |
440 | + call wwave( no_s, spin, no_u, no_l, maxnh, |
441 | . nkpnt, |
442 | . numh, listhptr, listh, H, S, Ef, xijo, indxuo, |
443 | . gamma_SCF, nkpnt, kpoint, no_u, occtol) |
444 | @@ -303,85 +310,101 @@ |
445 | ! Compute bands |
446 | !** This uses H,S, and xa |
447 | nullify( ebk ) |
448 | - call re_alloc( ebk, 1, no_u, 1, spinor_dim, 1, maxbk, |
449 | - & 'ebk', 'siesta_analysis' ) |
450 | - |
451 | - if (nbk.gt.0) then |
452 | + if ( spin%Grid == 4 ) then |
453 | + call re_alloc( ebk, 1, no_u*2, 1, 1, 1, maxbk, |
454 | + & 'ebk', 'siesta_analysis' ) |
455 | + else |
456 | + call re_alloc( ebk, 1, no_u, 1, spin%spinor, 1, maxbk, |
457 | + & 'ebk', 'siesta_analysis' ) |
458 | + end if |
459 | + |
460 | + if ( nbk > 0 ) then |
461 | if (IONode) print "(a)", "Computing bands..." |
462 | getPSI = fdf_get('WFS.Write.For.Bands', .false.) |
463 | if (getPSI) then |
464 | - wfs_filename = trim(slabel)//".bands.WFSX" |
465 | - if (IONode) print "(a)", "Writing WFSX for bands in " |
466 | - $ // trim(wfs_filename) |
467 | - max_n_states = no_u |
468 | - if (h_spin_dim >= 4) max_n_states = 2*no_u |
469 | - wfs_band_min = fdf_get("WFS.BandMin",1) |
470 | - wfs_band_max = fdf_get("WFS.BandMax",max_n_states) |
471 | - call setup_wfs_list(nbk,max_n_states, |
472 | - $ wfs_band_min,wfs_band_max, |
473 | - $ use_scf_weights=.false., |
474 | - $ use_energy_window=.false.) |
475 | - endif |
476 | - call bands( no_s, h_spin_dim, spinor_dim, no_u, no_l, maxnh, |
477 | + wfs_filename = trim(slabel)//".bands.WFSX" |
478 | + if (IONode) print "(a)", "Writing WFSX for bands in " |
479 | + $ // trim(wfs_filename) |
480 | + wfs_band_min = fdf_get("WFS.BandMin",1) |
481 | + wfs_band_max = fdf_get("WFS.BandMax",max_n_states) |
482 | + call setup_wfs_list(nbk,max_n_states, |
483 | + $ wfs_band_min,wfs_band_max, |
484 | + $ use_scf_weights=.false., |
485 | + $ use_energy_window=.false.) |
486 | + end if |
487 | + call bands( no_s, spin, no_u, no_l, maxnh, |
488 | . maxbk, |
489 | . numh, listhptr, listh, H, S, Ef, xijo, indxuo, |
490 | - . .true., nbk, bk, ebk, occtol, getPSI ) |
491 | - if (IOnode) then |
492 | + . .true., nbk, bk, ebk, occtol, getPSI ) |
493 | + |
494 | + if ( IOnode ) then |
495 | if ( writbk ) then |
496 | - write(6,'(/,a,/,a4,a12)') |
497 | - & 'siesta: Band k vectors (Bohr**-1):', 'ik', 'k' |
498 | + write(6,'(/,a,/,a4,tr1,a12)') |
499 | + & 'siesta: Band k vectors (Bohr**-1):', 'ik', 'k' |
500 | do ik = 1,nbk |
501 | - write(6,'(i4,3f12.6)') ik, (bk(ix,ik),ix=1,3) |
502 | - enddo |
503 | - endif |
504 | + write(6,'(i4,3(tr1,f12.6))') ik, (bk(ix,ik),ix=1,3) |
505 | + end do |
506 | + end if |
507 | |
508 | if ( writb ) then |
509 | - write(6,'(/,a,/,a4,a3,a7)') |
510 | - & 'siesta: Band energies (eV):', 'ik', 'is', 'eps' |
511 | - do ispin = 1, spinor_dim |
512 | + if ( spin%Grid == 4 ) then |
513 | + write(6,'(/,a,/,a4,tr1,a9)') |
514 | + & 'siesta: Band energies (eV):', 'ik', 'eps' |
515 | do ik = 1,nbk |
516 | - write(6,'(i4,i3,10f7.2)') |
517 | - & ik, ispin, (ebk(io,ispin,ik)/eV,io=1,min(10,no_u)) |
518 | - if (no_u.gt.10) write(6,'(7x,10f7.2)') |
519 | - & (ebk(io,ispin,ik)/eV,io=11,no_u) |
520 | + write(6,'(i4,10(tr1,f9.4))') |
521 | + & ik, (ebk(io,1,ik)/eV,io=1,min(10,2*no_u)) |
522 | + if ( 2*no_u > 10 ) write(6,'(4x,10(tr1,f9.4))') |
523 | + & (ebk(io,1,ik)/eV,io=11,no_u*2) |
524 | + end do |
525 | + |
526 | + else |
527 | + write(6,'(/,a,/,a4,a3,tr1,a9)') |
528 | + & 'siesta: Band energies (eV):', 'ik', 'is', 'eps' |
529 | + do ispin = 1, spin%spinor |
530 | + do ik = 1,nbk |
531 | + write(6,'(i4,i3,10(tr1,f9.4))') |
532 | + & ik, ispin, (ebk(io,ispin,ik)/eV,io=1,min(10,no_u)) |
533 | + if ( no_u > 10 ) write(6,'(7x,10(tr1,f9.4))') |
534 | + & (ebk(io,ispin,ik)/eV,io=11,no_u) |
535 | + enddo |
536 | enddo |
537 | - enddo |
538 | + end if |
539 | endif |
540 | + |
541 | endif |
542 | endif |
543 | |
544 | ! Print eigenvalues |
545 | - if (IOnode .and. writeig) then |
546 | - if ((isolve.eq.SOLVE_DIAGON .or. |
547 | - & ((isolve.eq.SOLVE_MINIM) .and. minim_calc_eigenvalues)) |
548 | - & .and. no_l.lt.1000) then |
549 | - if ( h_spin_dim <= 2 ) then |
550 | - write(6,'(/,a,/,a4,a3,a7)') |
551 | - & 'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps' |
552 | - do ik = 1,nkpnt |
553 | - do ispin = 1,spinor_dim |
554 | - write(6,'(i4,i3,10f7.2)') |
555 | + if ( IOnode .and. writeig .and. no_l < 1000 .and. |
556 | + & (isolve.eq.SOLVE_DIAGON .or. |
557 | + & (isolve.eq.SOLVE_MINIM .and. minim_calc_eigenvalues)) ) then |
558 | + |
559 | + if ( spin%Grid == 4 ) then |
560 | + write(6,'(/,a)') 'siesta: Eigenvalues (eV):' |
561 | + do ik = 1,nkpnt |
562 | + write(6,'(a,i6)') 'ik =', ik |
563 | + write(6,'(10(tr1,f9.4))') (eo(io,1,ik)/eV,io=1,2*neigwanted) |
564 | + end do |
565 | + |
566 | + else |
567 | + write(6,'(/,a,/,a4,a3,tr1,a9)') |
568 | + & 'siesta: Eigenvalues (eV):', 'ik', 'is', 'eps' |
569 | + do ik = 1,nkpnt |
570 | + do ispin = 1, spin%spinor |
571 | + write(6,'(i4,i3,10(tr1,f9.4))') |
572 | & ik,ispin,(eo(io,ispin,ik)/eV,io=1,min(10,neigwanted)) |
573 | - if (no_u.gt.10) write(6,'(7x,10f7.2)') |
574 | + if (no_u.gt.10) write(6,'(7x,10(tr1,f9.4))') |
575 | & (eo(io,ispin,ik)/eV,io=11,neigwanted) |
576 | - enddo |
577 | - enddo |
578 | - else |
579 | - write(6,'(/,a)') 'siesta: Eigenvalues (eV):' |
580 | - do ik = 1,nkpnt |
581 | - write(6,'(a,i6)') 'ik =', ik |
582 | - write(6,'(10f7.2)') |
583 | - & ((eo(io,ispin,ik)/eV,io=1,neigwanted),ispin=1,2) |
584 | - enddo |
585 | - endif |
586 | - write(6,'(a,f15.6,a)') 'siesta: Fermi energy =', ef/eV, ' eV' |
587 | - endif |
588 | + end do |
589 | + end do |
590 | + end if |
591 | + write(6,'(a,f15.6,a)') 'siesta: Fermi energy =', ef/eV, ' eV' |
592 | endif |
593 | |
594 | if (((isolve.eq.SOLVE_DIAGON).or. |
595 | & ((isolve.eq.SOLVE_MINIM).and.minim_calc_eigenvalues)) |
596 | & .and.IOnode) |
597 | - & call ioeig(eo,ef,neigwanted,nspin,nkpnt,no_u,spinor_dim, |
598 | + & call ioeig(eo,ef,neigwanted,spin%H,nkpnt,no_u,spin%spinor, |
599 | & nkpnt,kpoint, kweight) |
600 | |
601 | |
602 | @@ -437,7 +460,7 @@ |
603 | ! Calculation of the bulk polarization using the Berry phase |
604 | ! formulas by King-Smith and Vanderbilt |
605 | if (nkpol.gt.0 .and. .not.bornz) then |
606 | - if (NonCol .or. SpOrb) then |
607 | + if ( spin%NCol .or. spin%SO ) then |
608 | if (IOnode) then |
609 | write(6,'(/a)') |
610 | . 'siesta_analysis: bulk polarization implemented only for' |
611 | @@ -446,7 +469,7 @@ |
612 | endif |
613 | else |
614 | call KSV_pol(na_u, na_s, xa_last, rmaxo, scell_last, |
615 | - & ucell_last,no_u, no_l, no_s, nspin, qspin, |
616 | + & ucell_last,no_u, no_l, no_s, spin%Grid, qspin, |
617 | & maxnh, nkpol, numh, listhptr, listh, |
618 | & H, S, xijo, indxuo, isa, iphorb, |
619 | & iaorb, lasto, shape, |
620 | @@ -456,7 +479,7 @@ |
621 | |
622 | ! Calculation of the optical conductivity |
623 | call optical(na_u, na_s, xa_last, scell_last, ucell_last, |
624 | - & no_u, no_l, no_s, nspin, qspin, |
625 | + & no_u, no_l, no_s, spin%Grid, qspin, |
626 | & maxnh, numh, listhptr, listh, H, S, xijo, |
627 | $ indxuo, ebk, ef, temp, |
628 | & isa, iphorb, iphKB, lasto, lastkb, shape ) |
629 | @@ -486,7 +509,7 @@ |
630 | g2max = g2cut |
631 | dummy_str = 0.0 |
632 | dummy_strl = 0.0 |
633 | - call dhscf( nspin, no_s, iaorb, iphorb, no_l, |
634 | + call dhscf( spin%Grid, no_s, iaorb, iphorb, no_l, |
635 | . no_u, na_u, na_s, isa, xa_last, indxua, |
636 | & ntm, 0, 0, 0, filesOut, |
637 | & maxnh, numh, listhptr, listh, Dscf, Datm, |
638 | |
639 | === modified file 'Src/writewave.F' |
640 | --- Src/writewave.F 2019-06-17 08:37:13 +0000 |
641 | +++ Src/writewave.F 2019-08-08 11:56:55 +0000 |
642 | @@ -40,7 +40,7 @@ |
643 | USE alloc, only: re_alloc |
644 | USE sys, only: die |
645 | USE atomlist, only: no_u |
646 | - USE m_spin, only: h_spin_dim |
647 | + USE m_spin, only: spin |
648 | |
649 | implicit none |
650 | |
651 | @@ -51,11 +51,11 @@ |
652 | nullify(wfk) |
653 | |
654 | nwk = 0 |
655 | - if (h_spin_dim < 4) then |
656 | - max_n_states = no_u |
657 | + if ( spin%Grid == 4 ) then |
658 | + max_n_states = 2 * no_u |
659 | else |
660 | - max_n_states = 2*no_u |
661 | - endif |
662 | + max_n_states = no_u |
663 | + end if |
664 | call initwave( max_n_states, nwk, wfk ) |
665 | |
666 | if (nwk .eq. 0) then |
667 | @@ -391,7 +391,7 @@ |
668 | |
669 | end subroutine initwave |
670 | |
671 | - subroutine wwave( no, nspin, maxspn, maxo, maxuo, maxnh, |
672 | + subroutine wwave( no, spin, maxo, maxuo, maxnh, |
673 | . maxk, |
674 | . numh, listhptr, listh, H, S, ef, xij, indxuo, |
675 | . gamma, nk, kpoint, nuotot, occtol) |
676 | @@ -399,10 +399,10 @@ |
677 | C Finds wavefunctions at selected k-points. |
678 | C Written by P. Ordejon, June 2003 |
679 | C from routine 'bands' written by J.M.Soler |
680 | +C Changed to use tSpin, N. Papior, August 2019 |
681 | C **************************** INPUT ********************************** |
682 | C integer no : Number of basis orbitals |
683 | -C integer nspin = h_spin_dim : Number of spin components of H, D |
684 | -C integer maxspn = spinor_dim : Second dimension of ek |
685 | +C tSpin spin : Contains spin information |
686 | C integer maxo : First dimension of ek |
687 | C integer maxuo : Second dimension of H and S |
688 | C integer maxnh : Maximum number of orbitals interacting |
689 | @@ -414,7 +414,7 @@ |
690 | C hamiltonian matrix |
691 | C integer listh(maxlh) : Nonzero hamiltonian-matrix element |
692 | C column indexes for each matrix row |
693 | -C real*8 H(maxnh,nspin) : Hamiltonian in sparse form |
694 | +C real*8 H(maxnh,spin%H) : Hamiltonian in sparse form |
695 | C real*8 S(maxnh) : Overlap in sparse form |
696 | C real*8 ef : Fermi energy |
697 | C real*8 xij(3,maxnh) : Vectors between orbital centers (sparse) |
698 | @@ -446,7 +446,7 @@ |
699 | use atmfuncs, only : symfio, cnfigfio, labelfis, nofis |
700 | use atomlist, only : iaorb, iphorb |
701 | use siesta_geom, only : isa |
702 | - use m_spin, only : NoMagn, SPpol, NonCol, SpOrb |
703 | + use t_spin, only: tSpin |
704 | use units, only : eV |
705 | |
706 | #ifdef MPI |
707 | @@ -457,12 +457,12 @@ |
708 | |
709 | implicit none |
710 | |
711 | + type(tSpin), intent(in) :: spin |
712 | logical, intent(in) :: Gamma |
713 | integer maxk, maxnh, maxo, maxuo, nk, no, |
714 | - . h_spin_dim, spinor_dim, nspin, maxspn, |
715 | . nuotot, indxuo(no), listh(maxnh), numh(*), |
716 | . listhptr(*) |
717 | - real(dp) ef, H(:,:), kpoint(3,maxk), |
718 | + real(dp) ef, H(maxnh,spin%H), kpoint(3,maxk), |
719 | . S(maxnh), xij(3,maxnh), occtol |
720 | |
721 | external io_assign, io_close, memory |
722 | @@ -488,11 +488,7 @@ |
723 | data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/ |
724 | |
725 | logical :: SaveParallelOverK |
726 | - integer :: nspin_flag |
727 | |
728 | - h_spin_dim=size(H,dim=2) |
729 | - nspin_flag = h_spin_dim |
730 | - if (h_spin_dim == 8) nspin_flag = 4 |
731 | |
732 | C Get local number of orbitals |
733 | #ifdef MPI |
734 | @@ -517,15 +513,15 @@ |
735 | C Check spin |
736 | |
737 | C Allocate local arrays - only aux is relevant here |
738 | - if ( h_spin_dim >= 4 ) then |
739 | + if ( spin%Grid == 4 ) then |
740 | nhs = 2 * (2*nuotot) * (2*nuo) |
741 | else |
742 | nhs = 2 * nuotot * nuo |
743 | endif |
744 | call allocDenseMatrix(nhs, nhs, nhs) |
745 | |
746 | - allocate(ek(nuotot,maxspn,nk)) |
747 | - call memory('A','D',maxspn*nk*nuotot,'writewave') |
748 | + allocate(ek(nuotot,spin%spinor,nk)) |
749 | + call memory('A','D',nuotot*spin%spinor*nk,'writewave') |
750 | naux = 2*nuotot*5 |
751 | allocate(aux(naux)) |
752 | call memory('A','D',naux,'writewave') |
753 | @@ -544,12 +540,12 @@ |
754 | write(6,'(a)') 'writewave: Wave Functions Coefficients' |
755 | write(6,*) |
756 | write(6,'(a26,2x,i6)') 'Number of k-points = ', nk |
757 | - write(6,'(a26,2x,i6)') 'Number of Spins = ', nspin_flag |
758 | + write(6,'(a26,2x,i6)') 'Number of Spins = ', spin%Grid |
759 | write(6,'(a26,2x,i6)') 'Number of basis orbs = ',nuotot |
760 | write(6,*) |
761 | endif |
762 | write(iu) nk, gamma |
763 | - write(iu) nspin_flag ! This will be 1, 2, or 4 |
764 | + write(iu) spin%Grid ! This will be 1, 2, or 4 |
765 | write(iu) nuotot |
766 | write(iu) (iaorb(j),labelfis(isa(iaorb(j))), |
767 | . iphorb(j), cnfigfio(isa(iaorb(j)),iphorb(j)), |
768 | @@ -567,9 +563,9 @@ |
769 | |
770 | C Call appropriate diagonalization routine |
771 | |
772 | - if (NoMagn .or. SPpol) then |
773 | + if ( spin%Grid <= 2 ) then |
774 | if (gamma) then |
775 | - call diagg( h_spin_dim, nuo, maxuo, maxnh, maxnh, |
776 | + call diagg( spin%H, nuo, maxuo, maxnh, maxnh, |
777 | . maxo, numh, listhptr, listh, numh, listhptr, |
778 | . listh, H, S, getD, getPSI, |
779 | . fixspin, qtot, qs, temp, |
780 | @@ -589,7 +585,7 @@ |
781 | call diag_init() |
782 | end if |
783 | |
784 | - call diagk( h_spin_dim, nuo, no, maxspn, maxnh, maxnh, |
785 | + call diagk( spin%H, nuo, no, spin%spinor, maxnh, maxnh, |
786 | . maxo, numh, listhptr, listh, numh, listhptr, |
787 | . listh, H, S, getD, getPSI, |
788 | . fixspin, qtot, qs, temp, |
789 | @@ -603,7 +599,7 @@ |
790 | |
791 | endif |
792 | |
793 | - elseif (NonCol) then |
794 | + elseif ( spin%NCol) then |
795 | if (gamma) then |
796 | call diag2g( nuo, no, maxnh, maxnh, maxo, numh, |
797 | . listhptr, listh, numh, listhptr, listh, |
798 | @@ -620,7 +616,7 @@ |
799 | . Haux, Saux, psi, Haux, Saux, aux, |
800 | . nuotot, occtol, 1, nuotot ) |
801 | endif |
802 | - elseif (SpOrb) then |
803 | + elseif ( spin%SO ) then |
804 | if (gamma) then |
805 | call diag3g( nuo, no, maxnh, maxnh, maxo, numh, |
806 | . listhptr, listh, numh, listhptr, listh, |
807 | |
808 | === modified file 'version.info' |
809 | --- version.info 2019-06-20 13:50:42 +0000 |
810 | +++ version.info 2019-08-08 11:56:55 +0000 |
811 | @@ -1,1 +1,1 @@ |
812 | -siesta-4.1-1074--wfs-1115 |
813 | +siesta-4.1-1074--wfs-1115--nicpa-1 |
I agree with your tSpin and print_spin changes. They do make the code more clear. Thanks.
The new tSpin variable (equivalent to Ncol.or.SO), and other enhancements, can maybe be defined in the 4.2 branch.
Regarding the output files headings:
In bands and ioeig, no_u*2 refers to the number of states, and 1 to the
number of spin blocks (=1 since there is no spin distinction).
In the wavefunction files, no_u refers to the number of basis orbitals, and
4 to the number of expansions of length no_u needed at each band (4 for a complex
spinor). The number of states in the file is given by nk and the following records
for each k-point.
One possible source of fragility is that the "spinor_dim second-dimension" trick in ioeig and bands might not work when 'neigwanted' is not no_u. This could be also handled more robustly in 4.2, after the merge, and back-ported if deemed necessary for some use cases.