Merge lp:~albertog/siesta/4.0-images into lp:siesta/4.0
- 4.0-images
- Merge into rel-4.0
Status: | Superseded |
---|---|
Proposed branch: | lp:~albertog/siesta/4.0-images |
Merge into: | lp:siesta/4.0 |
Diff against target: |
839 lines (+478/-94) (has conflicts) 9 files modified
Src/Makefile (+7/-3) Src/diagk0.F (+295/-0) Src/diagon.F (+18/-2) Src/hsparse.F (+53/-76) Src/new_dm.F (+20/-1) Src/pdos.F (+9/-1) Src/state_init.F (+53/-8) Src/writewave.F (+14/-3) version.info (+9/-0) Text conflict in version.info |
To merge this branch: | bzr merge lp:~albertog/siesta/4.0-images |
Related bugs: |
Reviewer | Review Type | Date Requested | Status |
---|---|---|---|
Nick Papior | Approve | ||
Review via email: mp+344377@code.launchpad.net |
This proposal has been superseded by a proposal from 2018-05-08.
Commit message
(See individual commit messages)
Description of the change
Two independent changes:
1. The fix for the "severe folding" leading to wrong norm of the DM initialized by atomic data.
2. A saner hsparse.
I am not sure it is worth detecting the "severe folding" in hsparse (by checking for sameness of io and jo) for informative purposes.
Alberto Garcia (albertog) wrote : | # |
Nick Papior (nickpapior) wrote : | # |
Could we force the auxiliary cell for Harris functional?
1) I think very few people are using it,
2) The overhead is negligeble in this case.
Nick Papior (nickpapior) : | # |
- 576. By Alberto Garcia
-
Check for diagonal folding. Stop Harris runs if found
- 577. By Alberto Garcia
-
Add diagonalizer for gamma point with auxiliary supercell
The new routine 'diagk0' is a minor modification of 'diagg' that
can work with an auxiliary supercell while still using real
arithmetic. It is meant to be used when 'force-aux-cell T' is
specified in the input file, for gamma point calculations. (This
option can be set to carry out COOP/COHP analyses, or to avoid
severe folding in test cases.)The rest of the program logic is kept as it was. It might be worth
rethinking the use of the 'gamma' variable. - 578. By Alberto Garcia
-
Use 'diagk0' when appropriate in pdos and writewave
This might save time when computing the COOP-related wavefunctions,
and when carrying out a pDOS calculation with the 'gamma' sampling.
Unmerged revisions
- 578. By Alberto Garcia
-
Use 'diagk0' when appropriate in pdos and writewave
This might save time when computing the COOP-related wavefunctions,
and when carrying out a pDOS calculation with the 'gamma' sampling. - 577. By Alberto Garcia
-
Add diagonalizer for gamma point with auxiliary supercell
The new routine 'diagk0' is a minor modification of 'diagg' that
can work with an auxiliary supercell while still using real
arithmetic. It is meant to be used when 'force-aux-cell T' is
specified in the input file, for gamma point calculations. (This
option can be set to carry out COOP/COHP analyses, or to avoid
severe folding in test cases.)The rest of the program logic is kept as it was. It might be worth
rethinking the use of the 'gamma' variable.
Preview Diff
1 | === modified file 'Src/Makefile' |
2 | --- Src/Makefile 2018-04-11 11:00:02 +0000 |
3 | +++ Src/Makefile 2018-05-08 09:38:49 +0000 |
4 | @@ -84,7 +84,7 @@ |
5 | chempot.o coceri.o coxmol.o cross.o compute_norm.o\ |
6 | denmat.o denmatlomem.o detover.o dfscf.o m_diagon.o diagon.o digcel.o \ |
7 | fft.o dhscf.o constr.o diagk_file.o \ |
8 | - diagg.o diagk.o diagkp.o diag2g.o diag2k.o diagpol.o \ |
9 | + diagg.o diagk0.o diagk.o diagkp.o diag2g.o diag2k.o diagpol.o \ |
10 | diagsprl.o dipole.o dismin.o dnaefs.o doping_uniform.o dot.o \ |
11 | m_efield.o egandd.o ener3.o ener3lomem.o errorf.o extrapolon.o \ |
12 | fixed.o interpolation.o gradient.o gradientlomem.o grdsam.o \ |
13 | @@ -474,6 +474,8 @@ |
14 | diagg.o: writewave.o |
15 | diagk.o: compute_norm.o fermid.o parallel.o parallelsubs.o precision.o sys.o |
16 | diagk.o: writewave.o |
17 | +diagk0.o: alloc.o fermid.o parallel.o parallelsubs.o precision.o sys.o |
18 | +diagk0.o: writewave.o |
19 | diagk_file.o: fermid.o iowfs_netcdf.o parallel.o parallelsubs.o precision.o |
20 | diagk_file.o: sys.o writewave.o |
21 | diagkp.o: alloc.o fermid.o parallel.o parallelsubs.o precision.o sys.o |
22 | @@ -528,8 +530,8 @@ |
23 | grdsam.o: alloc.o dhscf.o files.o m_mpi_utils.o m_partial_charges.o parallel.o |
24 | grdsam.o: precision.o siesta_geom.o siesta_options.o sys.o units.o |
25 | hsparse.o: alloc.o atm_types.o atmfuncs.o atomlist.o listsc.o mneighb.o |
26 | -hsparse.o: parallel.o parallelsubs.o precision.o siesta_options.o sorting.o |
27 | -hsparse.o: sparse_matrices.o sys.o |
28 | +hsparse.o: parallel.o parallelsubs.o precision.o sorting.o sparse_matrices.o |
29 | +hsparse.o: sys.o |
30 | idiag.o: parallel.o sys.o |
31 | initatom.o: atmparams.o atom.o atom_options.o basis_io.o basis_specs.o |
32 | initatom.o: basis_types.o electrostatic.o old_atmfuncs.o precision.o |
33 | @@ -703,6 +705,8 @@ |
34 | ofc.o: alloc.o files.o precision.o units.o |
35 | old_atmfuncs.o: alloc.o atmparams.o precision.o sys.o |
36 | on_subs.o: alloc.o onmod.o onmod.o |
37 | +one.diagk0.o: compute_norm.o fermid.o parallel.o parallelsubs.o precision.o |
38 | +one.diagk0.o: sys.o writewave.o |
39 | onmod.o: precision.o |
40 | optical.o: alloc.o atomlist.o densematrix.o fermid.o files.o parallel.o |
41 | optical.o: parallelsubs.o precision.o sys.o units.o |
42 | |
43 | === added file 'Src/diagk0.F' |
44 | --- Src/diagk0.F 1970-01-01 00:00:00 +0000 |
45 | +++ Src/diagk0.F 2018-05-08 09:38:49 +0000 |
46 | @@ -0,0 +1,295 @@ |
47 | +! |
48 | +! Copyright (C) 1996-2016 The SIESTA group |
49 | +! This file is distributed under the terms of the |
50 | +! GNU General Public License: see COPYING in the top directory |
51 | +! or http://www.gnu.org/copyleft/gpl.txt. |
52 | +! See Docs/Contributors.txt for a list of contributors. |
53 | +! |
54 | + subroutine diagk0( nspin, nuo, maxuo, maxnh, maxnd, |
55 | + & maxo, numh, listhptr, listh, numd, |
56 | + & listdptr, listd, H, S, |
57 | + & getD, getPSI, fixspin, qtot, qs, temp, e1, e2, |
58 | + & no_s, indxuo, |
59 | + & eo, qo, Dnew, Enew, ef, efs, Entropy, |
60 | + & Haux, Saux, psi, nuotot, occtol, |
61 | + & iscf, neigwanted ) |
62 | +C ********************************************************************* |
63 | +C Subroutine to calculate the eigenvalues and eigenvectors, density |
64 | +C and energy-density matrices, and occupation weights of each |
65 | +C eigenvector, for given Hamiltonian and Overlap matrices (including |
66 | +C spin polarization). Gamma-point version. |
67 | +C Writen by J.Soler, August 1998. |
68 | +C **************************** INPUT ********************************** |
69 | +C integer nspin : Number of spin components (1 or 2) |
70 | +C integer nuo : Number of basis orbitals local to node |
71 | +C integer maxuo : Last dimension of xij |
72 | +C Must be at least max(indxuo) |
73 | +C integer maxnh : Maximum number of orbitals interacting |
74 | +C integer maxnd : Maximum number of nonzero elements of |
75 | +C each row of density matrix |
76 | +C integer maxo : First dimension of eo and qo |
77 | +C integer numh(nuo) : Number of nonzero elements of each row |
78 | +C of hamiltonian matrix |
79 | +C integer listhptr(nuo) : Pointer to each row (-1) of the |
80 | +C hamiltonian matrix |
81 | +C integer listh(maxnh) : Nonzero hamiltonian-matrix element |
82 | +C column indexes for each matrix row |
83 | +C integer numd(nuo) : Number of nonzero elements of each row |
84 | +C ofdensity matrix |
85 | +C integer listdptr(nuo) : Pointer to each row (-1) of the |
86 | +C density matrix |
87 | +C integer listd(maxnd) : Nonzero density-matrix element column |
88 | +C indexes for each matrix row |
89 | +C real*8 H(maxnh,nspin) : Hamiltonian in sparse form |
90 | +C real*8 S(maxnh) : Overlap in sparse form |
91 | +C logical getD : Find occupations and density matrices? |
92 | +C logical getPSI : Find and print wavefunctions? |
93 | +C logical fixspin : Fix the spin of the system? |
94 | +C real*8 qtot : Number of electrons in unit cell |
95 | +C real*8 qs(nspin) : Number of electrons in unit cell for each |
96 | +C spin component (if fixed spin option is used) |
97 | +C real*8 temp : Electronic temperature |
98 | +C real*8 e1, e2 : Energy range for density-matrix states |
99 | +C (to find local density of states) |
100 | +C Not used if e1 > e2 |
101 | +C integer no_s : Number of orbitals in auxiliary supercell |
102 | +C integer indxuo(no_s) : Map from aux supercell to unit cell |
103 | +C integer nuotot : total number of orbitals per unit cell |
104 | +C over all processors |
105 | +C integer iscf : SCF cycle number |
106 | +C real*8 occtol : Occupancy threshold for DM build |
107 | +C integer neigwanted : Number of eigenvalues wanted |
108 | +C *************************** OUTPUT ********************************** |
109 | +C real*8 eo(maxo,nspn) : Eigenvalues |
110 | +C ******************** OUTPUT (only if getD=.true.) ******************* |
111 | +C real*8 qo(maxo,nspn) : Occupations of eigenstates |
112 | +C real*8 Dnew(maxnd,nspin) : Output Density Matrix |
113 | +C real*8 Enew(maxnd,nspin) : Output Energy-Density Matrix |
114 | +C real*8 ef : Fermi energy |
115 | +C real*8 efs(nspin) : Fermi energy for each spin |
116 | +C (for fixed spin calculations) |
117 | +C real*8 Entropy : Electronic entropy |
118 | +C *************************** AUXILIARY ******************************* |
119 | +C real*8 Haux(nuotot,nuo) : Auxiliary space for the hamiltonian matrix |
120 | +C real*8 Saux(nuotot,nuo) : Auxiliary space for the overlap matrix |
121 | +C real*8 psi(nuotot,maxuo,nspin) : Auxiliary space for the eigenvectors |
122 | +C real*8 aux(nuotot) : Extra auxiliary space |
123 | +C *************************** UNITS *********************************** |
124 | +C eo, Enew and ef returned in the units of H. |
125 | +C *************************** PARALLEL ******************************** |
126 | +C The auxiliary arrays are now no longer symmetry and so the order |
127 | +C of referencing has been changed in several places to reflect this. |
128 | +C ********************************************************************* |
129 | +C |
130 | +C Modules |
131 | +C |
132 | + use precision |
133 | + use sys |
134 | + use parallel, only : Node, Nodes, BlockSize |
135 | + use parallelsubs, only : LocalToGlobalOrb |
136 | + use writewave, only : writew |
137 | + use m_fermid, only : fermid, fermispin, stepf |
138 | + use alloc |
139 | + |
140 | +#ifdef MPI |
141 | + use mpi_siesta |
142 | +#endif |
143 | + |
144 | + implicit none |
145 | + |
146 | +#ifdef MPI |
147 | + integer |
148 | + & MPIerror |
149 | +#endif |
150 | + |
151 | + integer |
152 | + & iscf, maxnd, maxnh, maxuo, maxo, nuo, nspin, nuotot, |
153 | + & neigwanted |
154 | + |
155 | + integer |
156 | + & listh(maxnh), numh(nuo), listhptr(nuo), |
157 | + & listd(maxnd), numd(nuo), listdptr(nuo) |
158 | + |
159 | + real(dp) |
160 | + & Dnew(maxnd,nspin), e1, e2, ef, Enew(maxnd,nspin), |
161 | + & Entropy, eo(maxo,nspin), H(maxnh,nspin), qo(maxo,nspin), |
162 | + & qtot, qs(nspin), S(maxnh), temp, efs(nspin), occtol |
163 | + |
164 | + real(dp) |
165 | + & Haux(nuotot,nuo), Saux(nuotot,nuo), psi(nuotot,maxuo,nspin) |
166 | + |
167 | + real(dp), dimension(1), parameter :: wk = (/ 1.0_dp /) |
168 | + integer, parameter :: nk = 1 |
169 | + |
170 | + integer no_s, indxuo(no_s) |
171 | + |
172 | + logical |
173 | + & getD, getPSI, fixspin |
174 | + |
175 | + external rdiag |
176 | + |
177 | +C Internal variables ............................................. |
178 | + integer ie, io, iio, ispin, ix, j, jo, BNode, iie, ind, |
179 | + & ierror, nd |
180 | + integer jo_s |
181 | + real(dp) ee, eei, qe, qei, rt, t, k(3) |
182 | + |
183 | + integer :: maxnuo, mm |
184 | + integer, pointer :: nuo_LOC(:) |
185 | + real(dp), pointer :: psi_tmp(:), paux(:) |
186 | + |
187 | +C Solve eigenvalue problem |
188 | + do ispin = 1,nspin |
189 | + call timer( 'r-eigvec', 1 ) |
190 | + call timer( 'r-buildHS', 1 ) |
191 | + do io = 1,nuo |
192 | + do jo = 1,nuotot |
193 | + Saux(jo,io) = 0.0d0 |
194 | + Haux(jo,io) = 0.0d0 |
195 | + enddo |
196 | + enddo |
197 | + do io = 1,nuo |
198 | + do j = 1,numh(io) |
199 | + ind = listhptr(io) + j |
200 | + jo_s = listh(ind) |
201 | + jo = indxuo(jo_s) ! In unit cell |
202 | + Saux(jo,io) = Saux(jo,io) + S(ind) |
203 | + Haux(jo,io) = Haux(jo,io) + H(ind,ispin) |
204 | + enddo |
205 | + enddo |
206 | + call timer( 'r-buildHS', 2 ) |
207 | + call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo(1,ispin), |
208 | + & psi(1,1,ispin), neigwanted, iscf, ierror ) |
209 | + call timer( 'r-eigvec', 2 ) |
210 | + |
211 | + |
212 | +C Check error flag and take appropriate action |
213 | + if (ierror.gt.0) then |
214 | + call die('Terminating due to failed diagonalisation') |
215 | + elseif (ierror.lt.0) then |
216 | +C Repeat diagonalisation with increased memory to handle clustering |
217 | + do io = 1,nuo |
218 | + do jo = 1,nuotot |
219 | + Saux(jo,io) = 0.0d0 |
220 | + Haux(jo,io) = 0.0d0 |
221 | + enddo |
222 | + enddo |
223 | + do io = 1,nuo |
224 | + do j = 1,numh(io) |
225 | + ind = listhptr(io) + j |
226 | + jo_s = listh(ind) |
227 | + jo = indxuo(jo_s) ! In unit cell |
228 | + Saux(jo,io) = Saux(jo,io) + S(ind) |
229 | + Haux(jo,io) = Haux(jo,io) + H(ind,ispin) |
230 | + enddo |
231 | + enddo |
232 | + call rdiag( Haux, Saux, nuotot, nuo, nuotot, eo(1,ispin), |
233 | + & psi(1,1,ispin), nuotot, iscf, ierror ) |
234 | + endif |
235 | + |
236 | + if (getPSI) then |
237 | + do ix = 1,3 |
238 | + k(ix) = 0.0d0 |
239 | + enddo |
240 | + call writew(nuotot,nuo,1,k,ispin, |
241 | + & eo(1,ispin),psi(1,1,ispin),.true.) |
242 | + endif |
243 | + |
244 | + enddo |
245 | +C Check if we are done ................................................ |
246 | + if (.not.getD) return |
247 | + |
248 | +C Find new Fermi energy and occupation weights ........................ |
249 | + if (fixspin) then |
250 | + call fermispin( nspin, nspin, nk, wk, maxo, neigwanted, |
251 | + & eo, temp, qs, qo, efs, Entropy ) |
252 | + else |
253 | + call fermid( nspin, nspin, nk, wk, maxo, neigwanted, eo, |
254 | + & temp, qtot, qo, ef, Entropy ) |
255 | + endif |
256 | + |
257 | +C Find weights for local density of states ............................ |
258 | + if (e1 .lt. e2) then |
259 | +* e1 = e1 - ef |
260 | +* e2 = e2 - ef |
261 | + t = max( temp, 1.d-6 ) |
262 | + rt = 1.0d0/t |
263 | + do ispin = 1,nspin |
264 | + do io = 1,nuotot |
265 | + qo(io,ispin) = ( stepf((eo(io,ispin)-e2)*rt) - |
266 | + & stepf((eo(io,ispin)-e1)*rt)) * |
267 | + & 2.0d0/dble(nspin) |
268 | + enddo |
269 | + enddo |
270 | + endif |
271 | + |
272 | + if (nuo.gt.0) then |
273 | +C New density and energy-density matrices of unit-cell orbitals ....... |
274 | + nd = listdptr(nuo) + numd(nuo) |
275 | + Dnew(1:nd,1:nspin) = 0.0d0 |
276 | + Enew(1:nd,1:nspin) = 0.0d0 |
277 | + endif |
278 | + |
279 | + call timer( 'r-buildD', 1 ) |
280 | +C Global operation to form new density matrix |
281 | + nullify( nuo_LOC ) |
282 | + call re_alloc( nuo_LOC, 0, Nodes-1, 'nuo_LOC', 'diagk0' ) |
283 | +#ifdef MPI |
284 | + call MPI_Allgather( nuo, 1, MPI_INTEGER, nuo_LOC, 1, |
285 | + & MPI_INTEGER, MPI_Comm_World, MPIerror ) |
286 | +#else |
287 | + nuo_LOC(0) = nuo |
288 | +#endif |
289 | + maxnuo = 0 |
290 | + do BNode=0, Nodes-1 |
291 | + maxnuo = max(nuo_LOC(BNode),maxnuo) |
292 | + enddo |
293 | + nullify( PSI_TMP ) |
294 | + call re_alloc( PSI_TMP, 1, nuotot*maxnuo*nspin, 'PSI_TMP', |
295 | + & 'diagk0' ) |
296 | + |
297 | + do Bnode=0, Nodes-1 |
298 | + if (BNode.eq.Node) then |
299 | + mm = 0 |
300 | + do ispin= 1, nspin |
301 | + do io= 1, nuo |
302 | + PSI_TMP(mm+1:mm+nuotot) = PSI(1:nuotot,io, ispin) |
303 | + mm = mm + nuotot |
304 | + enddo |
305 | + enddo |
306 | + endif |
307 | +#ifdef MPI |
308 | + call MPI_Bcast( PSI_TMP, nuotot*nuo_LOC(Bnode)*nspin, |
309 | + & MPI_double_precision, BNode, MPI_Comm_World, |
310 | + & MPIerror ) |
311 | +#endif |
312 | + do ispin = 1,nspin |
313 | + do io = 1,nuo |
314 | + call LocalToGlobalOrb(io,Node,Nodes,iio) |
315 | + mm = (ispin-1)*nuo_LOC(Bnode)*nuotot |
316 | + do iie = 1,nuo_LOC(Bnode) |
317 | + call LocalToGlobalOrb(iie,BNode,Nodes,ie) |
318 | + qe = qo(ie,ispin) |
319 | + if (abs(qe).gt.occtol) then |
320 | + paux => PSI_TMP(mm+1:mm+nuotot) |
321 | + ee = qe*eo(ie,ispin) |
322 | + qei = qe*paux(iio) |
323 | + eei = ee*paux(iio) |
324 | + do j = 1,numd(io) |
325 | + ind = listdptr(io) + j |
326 | + jo_s = listd(ind) |
327 | + jo = indxuo(jo_s) ! In unit cell |
328 | + Dnew(ind,ispin) = Dnew(ind,ispin) + qei*paux(jo) |
329 | + Enew(ind,ispin) = Enew(ind,ispin) + eei*paux(jo) |
330 | + enddo |
331 | + endif |
332 | + mm = mm + nuotot |
333 | + enddo |
334 | + enddo |
335 | + enddo |
336 | + enddo |
337 | + |
338 | + call timer( 'r-buildD', 2 ) |
339 | + call de_alloc( nuo_LOC, 'nuo_LOC', 'diagk0' ) |
340 | + call de_alloc( PSI_TMP, 'PSI_TMP', 'diagk0' ) |
341 | + end |
342 | |
343 | === modified file 'Src/diagon.F' |
344 | --- Src/diagon.F 2017-02-21 07:36:32 +0000 |
345 | +++ Src/diagon.F 2018-05-08 09:38:49 +0000 |
346 | @@ -126,7 +126,7 @@ |
347 | . fixspin, gamma, getD, getPSI |
348 | |
349 | external |
350 | - . diagg, diagk, diag2g, diag2k, readsp, diagsprl |
351 | + . diagg, diagk0, diagk, diag2g, diag2k, readsp, diagsprl |
352 | #ifdef CDF |
353 | external diagk_file |
354 | #endif |
355 | @@ -255,6 +255,7 @@ |
356 | . eo, qo, Dnew, Enew, ef, efs, Entropy, |
357 | . Haux, Saux, psi, Haux, Saux, aux, |
358 | . nuotot, occtol, iscf ) |
359 | + call output_kpoint_occupation() |
360 | else |
361 | #endif |
362 | #ifndef CDF |
363 | @@ -279,8 +280,22 @@ |
364 | . eo, qo, Dnew, Enew, ef, efs, Entropy, |
365 | . Haux, Saux, psi, Haux, Saux, aux, |
366 | . nuotot, occtol, iscf, neigwanted ) |
367 | + call output_kpoint_occupation() |
368 | #endif |
369 | else |
370 | + if ( (nk==1) .and. |
371 | + $ sum(kpoint(:,1)) == 0.0_dp) then |
372 | +! |
373 | +! Use special version (This is still gamma) |
374 | +! |
375 | + call diagk0( nspin, nuo, maxuo, maxnh, maxnd, maxo, |
376 | + . numh, listhptr, listh, numd, listdptr, listd, |
377 | + . H, S, getD, getPSI, fixspin, qtot, qs, temp, e1, e2, |
378 | + . no, indxuo, |
379 | + . eo, qo, Dnew, Enew, ef, efs, Entropy, |
380 | + . Haux, Saux, psi, nuotot, occtol, iscf, |
381 | + . neigwanted ) |
382 | + else |
383 | call diagk( nspin, nuo, no, maxspn, maxnh, maxnd, |
384 | . maxo, numh, listhptr, listh, numd, listdptr, |
385 | . listd, H, S, getD, getPSI, fixspin, qtot, qs, temp, |
386 | @@ -288,8 +303,9 @@ |
387 | . eo, qo, Dnew, Enew, ef, efs, Entropy, |
388 | . Haux, Saux, psi, Haux, Saux, aux, |
389 | . nuotot, occtol, iscf ) |
390 | + call output_kpoint_occupation() |
391 | + endif |
392 | endif |
393 | - call output_kpoint_occupation() |
394 | #ifdef MPI |
395 | endif |
396 | #endif |
397 | |
398 | === modified file 'Src/hsparse.F' |
399 | --- Src/hsparse.F 2016-01-25 16:00:16 +0000 |
400 | +++ Src/hsparse.F 2018-05-08 09:38:49 +0000 |
401 | @@ -48,7 +48,8 @@ |
402 | subroutine hsparse( negl, cell, nsc, na, isa, xa, |
403 | & lasto, lastkb, iphorb, iphkb, |
404 | & nlhmax, gamma, |
405 | - $ set_xijo, folding) |
406 | + $ set_xijo, folding, |
407 | + $ diagonal_folding, debug_folding) |
408 | |
409 | use precision |
410 | use parallel, only : Node, Nodes |
411 | @@ -61,7 +62,7 @@ |
412 | use neighbour, only : mneighb, reset_neighbour_arrays |
413 | use sys, only : die |
414 | use alloc, only : re_alloc, de_alloc |
415 | - use atomlist, only : no_l |
416 | + use atomlist, only : no_l, indxuo |
417 | use sparse_matrices, only : listhptr, numh, listh, xijo |
418 | implicit none |
419 | |
420 | @@ -76,6 +77,8 @@ |
421 | logical, intent(in) :: gamma |
422 | logical, intent(in) :: set_xijo |
423 | logical, intent(out) :: folding |
424 | + logical, intent(out) :: diagonal_folding |
425 | + logical, intent(in), optional :: debug_folding |
426 | |
427 | external timer |
428 | |
429 | @@ -99,8 +102,6 @@ |
430 | & rci, rcj, rck, rij, rik, rjk, |
431 | & rmax, rmaxkb, rmaxo |
432 | |
433 | - logical, save :: warn1 = .false. |
434 | - |
435 | real(dp), dimension(:), pointer :: rckb |
436 | logical, dimension(:), pointer :: conect |
437 | integer, dimension(:), pointer :: index |
438 | @@ -108,14 +109,18 @@ |
439 | integer, dimension(:), pointer :: listhtmp |
440 | |
441 | logical :: connected |
442 | - |
443 | + logical :: debug |
444 | + integer :: nprints ! for debugging |
445 | C ------------------------------------- |
446 | -#ifdef DEBUG |
447 | - call write_debug( ' PRE hsparse' ) |
448 | -#endif |
449 | + |
450 | C Start time counter |
451 | call timer( 'hsparse', 1 ) |
452 | |
453 | + debug = .false. |
454 | + if (present(debug_folding)) then |
455 | + debug = debug_folding |
456 | + endif |
457 | + |
458 | C Check size of internal arrays |
459 | ncells = nsc(1) * nsc(2) * nsc(3) |
460 | nua = na / ncells |
461 | @@ -192,6 +197,7 @@ |
462 | enddo |
463 | |
464 | folding = .false. |
465 | + diagonal_folding = .false. |
466 | |
467 | C------------------------------------C |
468 | C Find number of non-zeros in H C |
469 | @@ -324,6 +330,7 @@ |
470 | C------------------------------------C |
471 | C Find full H sparsity pattern C |
472 | C------------------------------------C |
473 | + nprints = 0 |
474 | C Loop on atoms in unit cell |
475 | do ia = 1,nua |
476 | C Find neighbour atoms within maximum range |
477 | @@ -396,62 +403,53 @@ |
478 | endif |
479 | |
480 | if (connected) then |
481 | - if (conect(jo)) then |
482 | - folding = .true. |
483 | - |
484 | - ! This test is now deferred to be able |
485 | - ! to catch multiple images while avoiding |
486 | - ! false positives (i.e., we test first |
487 | - ! whether there is indeed a connection). |
488 | - |
489 | - if (.true.) then |
490 | - ! If already connected and using supercell, |
491 | - ! the latter might not be big enough... |
492 | - ! We warn the user and keep the first instance |
493 | + if (conect(jo)) then |
494 | + folding = .true. |
495 | + if (io == indxuo(jo)) then |
496 | + diagonal_folding = .true. |
497 | + endif |
498 | + |
499 | + ! If already connected and using supercell, the |
500 | + ! latter might not be big enough... We defer |
501 | + ! the decision on whether to kill the program |
502 | + ! to the caller. Here keep the first instance |
503 | ! of xij (same behavior as the old xijorb, as |
504 | ! earlier jnats are closer) |
505 | - ! Warn also if Gamma-point calculation, just |
506 | - ! in case |
507 | - if (.not.warn1) then |
508 | - if (Node.eq.0) then |
509 | - if (gamma) then |
510 | - call check_cohp(io,jo) |
511 | - else |
512 | - write(6,'(/,a,2i6,a,/,a)') |
513 | - . 'WARNING: orbital pair ',io,jo, |
514 | + if (debug) then |
515 | + ! This might not be pretty in |
516 | + ! parallel. Better to build a list of a few |
517 | + ! cases, and pass it to the caller. |
518 | + if (nprints <= 20) then |
519 | + print "(a,2i6,a)", |
520 | + . 'WARNING: orbital pair ',io,indxuo(jo), |
521 | . ' is multiply connected' |
522 | - endif |
523 | + nprints = nprints + 1 |
524 | endif |
525 | - warn1 = .true. |
526 | - endif |
527 | - endif |
528 | - |
529 | - else |
530 | - |
531 | - conect(jo) = .true. |
532 | - numh(iio) = numh(iio) + 1 |
533 | - ind = listhptr(iio)+numh(iio) |
534 | - listh(ind) = jo |
535 | - if (set_xijo) then |
536 | - xijo(1:3,ind) = xij(1:3,jnat) |
537 | + endif |
538 | + |
539 | + else |
540 | + |
541 | + conect(jo) = .true. |
542 | + numh(iio) = numh(iio) + 1 |
543 | + ind = listhptr(iio)+numh(iio) |
544 | + listh(ind) = jo |
545 | + if (set_xijo) then |
546 | + xijo(1:3,ind) = xij(1:3,jnat) |
547 | + endif |
548 | endif |
549 | endif |
550 | - endif |
551 | + enddo |
552 | enddo |
553 | - enddo |
554 | |
555 | C Restore conect array for next orbital io |
556 | - do j = 1,numh(iio) |
557 | - jo = listh(listhptr(iio)+j) |
558 | - conect(jo) = .false. |
559 | - enddo |
560 | - endif ! iio > 0 |
561 | - enddo ! io |
562 | - enddo ! ia |
563 | + do j = 1,numh(iio) |
564 | + jo = listh(listhptr(iio)+j) |
565 | + conect(jo) = .false. |
566 | + enddo |
567 | + endif ! iio > 0 |
568 | + enddo ! io |
569 | + enddo ! ia |
570 | |
571 | -!! print "(a5,i3,a40,3i8)", |
572 | -!! $ "Node: ", Node, "in hsparse nuo, nuotot, nlhmax: ", |
573 | -!! $ nuo, nuotot, nlhmax |
574 | |
575 | C Initialize listsc |
576 | call LISTSC_INIT( nsc, nuotot ) |
577 | @@ -466,28 +464,7 @@ |
578 | call de_alloc( rckb, 'rckb', 'hsparse' ) |
579 | |
580 | call timer( 'hsparse', 2 ) |
581 | -#ifdef DEBUG |
582 | - call write_debug( ' POS hsparse' ) |
583 | -#endif |
584 | + |
585 | end subroutine hsparse |
586 | - |
587 | - subroutine check_cohp(io,jo) |
588 | - use siesta_options, only: write_coop |
589 | - |
590 | - integer, intent(in) :: io, jo |
591 | - |
592 | - if (write_coop) then |
593 | - write(6,'(/,a,2i6,a,/,a)') |
594 | - . 'NOTE: orbital pair ',io,jo, |
595 | - . ' (at least) is multiply connected.', |
596 | - . 'NOTE: Your COOP/COHP analysis might ' // |
597 | - $ 'be affected by folding.' |
598 | - write(0,'(/,a,2i6,a,/,a)') |
599 | - . 'NOTE: orbital pair ',io,jo, |
600 | - . ' (at least) is multiply connected.', |
601 | - . 'NOTE: Your COOP/COHP analysis might ' // |
602 | - $ 'be affected by folding.' |
603 | - endif |
604 | - end subroutine check_cohp |
605 | |
606 | end module m_hsparse |
607 | |
608 | === modified file 'Src/new_dm.F' |
609 | --- Src/new_dm.F 2018-04-06 12:32:04 +0000 |
610 | +++ Src/new_dm.F 2018-05-08 09:38:49 +0000 |
611 | @@ -485,6 +485,7 @@ |
612 | use mpi_siesta |
613 | #endif |
614 | use units, only : pi |
615 | + use sparse_matrices, only: S |
616 | |
617 | implicit none |
618 | |
619 | @@ -655,7 +656,7 @@ |
620 | endif |
621 | |
622 | |
623 | - else |
624 | + else ! Initialize with neutral atoms |
625 | |
626 | C See whether specific initial spins are given in a DM.InitSpin block |
627 | C and read them in a loop on atoms where lines are read and parsed |
628 | @@ -903,6 +904,24 @@ |
629 | |
630 | endif |
631 | |
632 | + ! We have initialized with atomic information. Correct in case we |
633 | + ! are using such a small cell that there are direct interactions |
634 | + ! of orbitals with their own images, and we insist on using the |
635 | + ! Gamma-point only. Otherwise S(diagonal) is always 1.0 and the |
636 | + ! simple atomic-orbital superposition works as intended. |
637 | + |
638 | + |
639 | + do io = 1, no_l |
640 | + call LocalToGlobalOrb(io,Node,Nodes,iio) |
641 | + do in = 1,numh(io) |
642 | + ind = listhptr(io)+in |
643 | + jo = listh(ind) |
644 | + if (iio .eq. jo) then ! diagonal element |
645 | + Dscf(ind,:) = Dscf(ind,:) / S(ind) |
646 | + endif |
647 | + enddo |
648 | + enddo |
649 | + |
650 | endif |
651 | |
652 | end subroutine initdm |
653 | |
654 | === modified file 'Src/pdos.F' |
655 | --- Src/pdos.F 2018-05-07 07:36:34 +0000 |
656 | +++ Src/pdos.F 2018-05-08 09:38:49 +0000 |
657 | @@ -180,11 +180,19 @@ |
658 | . XIJ, INDXUO, NK, KPOINT, WK, EO, |
659 | . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) |
660 | else |
661 | - call pdosk( NSPIN, NUO, NO, MAXSPN, MAXNH, |
662 | + if ((nk==1) .and. (sum(kpoint(:,1)) == 0.0d0)) then |
663 | + ! pdosg already setup for auxiliary supercell!! |
664 | + call pdosg( NSPIN, NUO, NO, MAXSPN, MAXNH, |
665 | + . MAXO, NUMH, LISTHPTR, LISTH, H, S, |
666 | + . E1, E2, NHIST, SIGMA, INDXUO, EO, |
667 | + . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) |
668 | + else |
669 | + call pdosk( NSPIN, NUO, NO, MAXSPN, MAXNH, |
670 | . MAXO, NUMH, LISTHPTR, LISTH, H, S, |
671 | . E1, E2, NHIST, SIGMA, |
672 | . XIJ, INDXUO, NK, KPOINT, WK, EO, |
673 | . HAUX, SAUX, PSI, DTOT, DPR, NO_U ) |
674 | + endif |
675 | endif |
676 | endif |
677 | |
678 | |
679 | === modified file 'Src/state_init.F' |
680 | --- Src/state_init.F 2018-04-11 09:23:27 +0000 |
681 | +++ Src/state_init.F 2018-05-08 09:38:49 +0000 |
682 | @@ -48,7 +48,7 @@ |
683 | use m_mpi_utils, only: globalize_or |
684 | use m_mpi_utils, only: globalize_sum |
685 | use domain_decom, only: domainDecom, use_dd, use_dd_perm |
686 | - |
687 | + use fdf, only: fdf_get |
688 | #ifdef TRANSIESTA |
689 | use m_ts_options, only : onlyS |
690 | use sys, only : bye |
691 | @@ -69,6 +69,7 @@ |
692 | integer :: i, ix, iadispl, ixdispl |
693 | logical :: auxchanged ! Auxiliary supercell changed? |
694 | logical :: folding, folding1 |
695 | + logical :: diag_folding, diag_folding1 |
696 | logical :: foundxv ! dummy for call to ioxv |
697 | real(dp) :: dummy_qspin(8) |
698 | external :: madelung, timer |
699 | @@ -255,7 +256,39 @@ |
700 | ! analyses. |
701 | call hsparse( negl, scell, nsc, na_s, isa, xa, lasto, |
702 | & lastkb, iphorb, iphKB, maxnh, gamma, |
703 | - $ set_xijo=.true., folding=folding1) |
704 | + $ set_xijo=.true., folding=folding1, |
705 | + $ diagonal_folding=diag_folding1, |
706 | + $ debug_folding=fdf_get('debug-folding',.false.)) |
707 | +! |
708 | + call globalize_or(diag_folding1,diag_folding) |
709 | + if (diag_folding .and. gamma) then |
710 | + if (IOnode) then |
711 | + write(6,"(a)") "Gamma-point calculation " // |
712 | + $ "with interaction between periodic images!" |
713 | + write(6,"(a)") "Some features might not work optimally" |
714 | + endif |
715 | + if (harrisfun) call die("Harris functional run needs " // |
716 | + $ "'force-aux-cell T'") |
717 | + endif |
718 | + ! |
719 | + call globalize_or(folding1,folding) |
720 | + if (folding) then |
721 | + if (gamma) then |
722 | + if (IOnode) then |
723 | + write(6,"(a)") "Gamma-point calculation " // |
724 | + $ "with multiply-connected orbital pairs" |
725 | + write(6,"(a)") "Folding of H and S implicitly performed" |
726 | + call check_cohp() |
727 | + endif |
728 | + else |
729 | + write(6,"(a,/,a)") "Non Gamma-point calculation " // |
730 | + $ "with multiply-connected orbital pairs " // |
731 | + $ "in auxiliary supercell.", |
732 | + $ "Possible internal error. " // |
733 | + $ "Use 'debug-folding T' to debug." |
734 | + call die("Inadequate auxiliary supercell") |
735 | + endif |
736 | + endif |
737 | ! |
738 | call globalize_sum(maxnh,nnz) |
739 | if (cml_p) then |
740 | @@ -269,12 +302,6 @@ |
741 | call cmlEndPropertyList(mainXML) |
742 | endif |
743 | ! |
744 | - call globalize_or(folding1,folding) |
745 | - if (folding) then |
746 | - if (IOnode) then |
747 | - print *, "Folding of H and S is implicitly performed" |
748 | - endif |
749 | - endif |
750 | ! |
751 | ! If using domain decomposition, redistribute orbitals |
752 | ! for this geometry, based on the hsparse info. |
753 | @@ -387,4 +414,22 @@ |
754 | |
755 | !--------------------------------------------------------------------------- END |
756 | END subroutine state_init |
757 | + |
758 | + subroutine check_cohp() |
759 | + use siesta_options, only: write_coop |
760 | + |
761 | + if (write_coop) then |
762 | + write(6,'(/,a,a,/,a)') |
763 | + . 'NOTE: There are multiply-connected orbitals.', |
764 | + . 'NOTE: Your COOP/COHP analysis might ' // |
765 | + $ 'be affected by folding.', |
766 | + . 'NOTE: Use "force-aux-cell T" or k-point sampling' |
767 | + write(0,'(/,a,a,/,a)') |
768 | + . 'NOTE: There are multiply-connected orbitals.', |
769 | + . 'NOTE: Your COOP/COHP analysis might ' // |
770 | + $ 'be affected by folding.', |
771 | + . 'NOTE: Use "force-aux-cell T" or k-point sampling' |
772 | + endif |
773 | + end subroutine check_cohp |
774 | + |
775 | END module m_state_init |
776 | |
777 | === modified file 'Src/writewave.F' |
778 | --- Src/writewave.F 2017-10-04 09:48:27 +0000 |
779 | +++ Src/writewave.F 2018-05-08 09:38:49 +0000 |
780 | @@ -471,6 +471,7 @@ |
781 | data Dnew, Enew, e1, e2, qk, qtot, temp, wk /8*0.d0/ |
782 | |
783 | logical SaveParallelOverK |
784 | + external :: diagg, diagk0, diagk |
785 | |
786 | C Get local number of orbitals |
787 | #ifdef MPI |
788 | @@ -600,8 +601,17 @@ |
789 | ParallelOverK = .false. |
790 | ResetFirstCall = .true. |
791 | |
792 | - |
793 | - call diagk( nspin, nuo, no, maxspn, maxnh, maxnh, |
794 | + if ( (nk==1) .and. (sum(kpoint(:,1)) == 0.0d0)) then |
795 | + call diagk0( nspin, nuo, maxuo, maxnh, maxnh, |
796 | + . maxo, numh, listhptr, listh, numh, listhptr, |
797 | + . listh, H, S, getD, getPSI, |
798 | + . fixspin, qtot, qs, temp, |
799 | + . e1, e2, |
800 | + $ no, indxuo, |
801 | + $ ek, qk, Dnew, Enew, ef, efs, Entropy, |
802 | + . Haux, Saux, psi, nuotot, occtol, 1, nuotot ) |
803 | + else |
804 | + call diagk( nspin, nuo, no, maxspn, maxnh, maxnh, |
805 | . maxo, numh, listhptr, listh, numh, listhptr, |
806 | . listh, H, S, getD, getPSI, |
807 | . fixspin, qtot, qs, temp, |
808 | @@ -609,6 +619,7 @@ |
809 | . ek, qk, Dnew, Enew, ef, efs, Entropy, |
810 | . Haux, Saux, psi, Haux, Saux, aux, nuotot, |
811 | . occtol, 1 ) |
812 | + endif |
813 | ParallelOverK = SaveParallelOverK |
814 | ResetFirstCall = .false. |
815 | |
816 | @@ -663,7 +674,7 @@ |
817 | real(SP), dimension(:,:), allocatable :: aux !! NOTE SP |
818 | |
819 | external io_assign, io_close |
820 | - |
821 | + |
822 | C ................... |
823 | |
824 | C Allocate auxiliary arrays |
825 | |
826 | === modified file 'version.info' |
827 | --- version.info 2018-05-07 12:07:29 +0000 |
828 | +++ version.info 2018-05-08 09:38:49 +0000 |
829 | @@ -1,1 +1,10 @@ |
830 | +<<<<<<< TREE |
831 | siesta-4.0--575 |
832 | +======= |
833 | +siesta-4.0--573--folding-578 |
834 | + |
835 | + |
836 | + |
837 | + |
838 | + |
839 | +>>>>>>> MERGE-SOURCE |
For Harris-functional calculations, severe folding is just plain crazy, as there is no time to smooth out the initial differences. We should remove the Harris option from the force_2 test.
In general, we should remove tests with non-physically small cells with gamma point (or at least flag them as showcases for the problems).