Context Navigation

restart.F90 @ 2908

Revision 2908, 19.7 KB (checked in by acastro, 3 years ago)
I moved the function restart_look to the states module, where it belongs. I also collected some code that was replicated in various places, and put it as a "restart_look_and_read" subroutine.
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`

Line
1	!! Copyright (C) 2002-2006 M. Marques, A. Castro, A. Rubio, G. Bertsch
2	!!
3	!! This program is free software; you can redistribute it and/or modify
4	!! it under the terms of the GNU General Public License as published by
5	!! the Free Software Foundation; either version 2, or (at your option)
6	!! any later version.
7	!!
8	!! This program is distributed in the hope that it will be useful,
9	!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10	!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11	!! GNU General Public License for more details.
12	!!
13	!! You should have received a copy of the GNU General Public License
14	!! along with this program; if not, write to the Free Software
15	!! Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
16	!! 02111-1307, USA.
17	!!
18	!! $Id$
19
20	#include "global.h"
21
22	module restart_m
23	use curvlinear_m
24	use datasets_m
25	use geometry_m
26	use global_m
27	use grid_m
28	use io_m
29	use lib_oct_m
30	use lib_oct_parser_m
31	use linear_response_m
32	use math_m
33	use mesh_function_m
34	use mesh_m
35	use messages_m
36	use mpi_m
37	use output_m
38	use simul_box_m
39	use states_m
40	use units_m
41	use varinfo_m
42
43	implicit none
44
45	private
46	public :: &
47	restart_init, &
48	clean_stop, &
49	restart_write, &
50	restart_read, &
51	restart_format, &
52	restart_look_and_read, &
53	drestart_write_function, &
54	zrestart_write_function, &
55	drestart_read_function, &
56	zrestart_read_function, &
57	drestart_write_lr_rho, &
58	zrestart_write_lr_rho, &
59	drestart_read_lr_rho, &
60	zrestart_read_lr_rho
61
62
63
64	integer, parameter :: &
65	RESTART_PLAIN = 1, &
66	RESTART_NETCDF = 2, &
67	RESTART_BINARY = 3
68
69	integer :: restart_format
70
71
72	contains
73
74	! returns true if a file named stop exists
75	function clean_stop()
76	logical clean_stop, file_exists
77
78	clean_stop = .false.
79	inquire(file='stop', exist=file_exists)
80	if(file_exists) then
81	message(1) = 'Clean STOP'
82	call write_warning(1)
83	clean_stop = .true.
84	#if defined(HAVE_MPI)
85	call MPI_Barrier(mpi_world%comm, mpi_err)
86	#endif
87	if(mpi_grp_is_root(mpi_world)) call loct_rm('stop')
88	end if
89
90	return
91	end function clean_stop
92
93
94	! ---------------------------------------------------------
95	! read restart format information
96	subroutine restart_init
97
98	integer :: default, parsed
99
100	call push_sub('restart.restart_init')
101
102	!%Variable RestartFileFormat
103	!%Type integer
104	!%Default restart_netcdf
105	!%Section Generalities::IO
106	!%Description
107	!% Determines in which format the restart files should be written. The default
108	!% is netcdf files (restart_netcdf) if the executable was compiled with netcdf
109	!% libraries; otherwise the default is plain binary format.
110	!%Option restart_plain 1
111	!% Binary (platform dependent) format
112	!%Option restart_netcdf 2
113	!% NetCDF (platform independent) format. This requires the NETCDF library.
114	!%Option restart_binary 3
115	!% Octopus Binary Format, the new and more flexible binary format
116	!% of Octopus. It is faster and produces smaller files than NetCDF
117	!% and it is compiler independent (in the future it will be
118	!% endianness independent).
119	!%End
120
121	#if defined(HAVE_NETCDF)
122	default = RESTART_NETCDF
123	#else
124	default = RESTART_BINARY
125	#endif
126
127	call loct_parse_int(check_inp('RestartFileFormat'), default, parsed)
128	#ifndef HAVE_NETCDF
129	if (parsed == RESTART_NETCDF) then
130	write(message(1),'(a)') 'Error: Octopus was compiled without NetCDF support but'
131	write(message(2),'(a)') 'NetCDF restart files were requested.'
132	call write_fatal(2)
133	end if
134	#endif
135	if(.not.varinfo_valid_option('RestartFileFormat', parsed)) call input_error('RestartFileFormat')
136
137	select case(parsed)
138	case(RESTART_NETCDF)
139	restart_format = output_fill_how("NETCDF")
140	case(RESTART_PLAIN)
141	restart_format = output_fill_how("Plain")
142	case(RESTART_BINARY)
143	restart_format = output_fill_how("Binary")
144	end select
145
146	call pop_sub()
147	end subroutine restart_init
148
149
150	! ---------------------------------------------------------
151	subroutine restart_look_and_read(dir, st, gr, geo, ierr)
152	character(len=*), intent(in) :: dir
153	type(states_t), intent(inout) :: st
154	type(grid_t), intent(in) :: gr
155	type(geometry_t), intent(in) :: geo
156	integer, intent(out) :: ierr
157
158	integer :: kpoints, dim, nst, j
159
160	call push_sub('restart.restart_look_and_read')
161
162	call states_look(trim(tmpdir)//'restart_gs', gr%m, kpoints, dim, nst, j)
163	if(j.ne.0) then
164	ierr = j
165	call pop_sub(); return
166	end if
167
168	st%nst = nst
169	st%st_end = nst
170	deallocate(st%eigenval, st%occ)
171
172	call states_allocate_wfns(st, gr%m)
173
174	ALLOCATE(st%eigenval(st%nst, st%d%nik), st%nst*st%d%nik)
175	ALLOCATE(st%occ(st%nst, st%d%nik), st%nst*st%d%nik)
176
177	if(st%d%ispin == SPINORS) then
178	ALLOCATE(st%spin(3, st%nst, st%d%nik), st%nstst%d%nik3)
179	st%spin = M_ZERO
180	end if
181	st%eigenval = huge(REAL_PRECISION)
182	st%occ = M_ZERO
183
184	call restart_read(trim(tmpdir)//'restart_gs', st, gr, geo, ierr)
185
186	call pop_sub()
187	end subroutine restart_look_and_read
188
189
190	! ---------------------------------------------------------
191	subroutine restart_write(dir, st, gr, ierr, iter, lr)
192	character(len=*), intent(in) :: dir
193	type(states_t), intent(in) :: st
194	type(grid_t), intent(in) :: gr
195	integer, intent(out) :: ierr
196	integer, optional, intent(in) :: iter
197	!if this next argument is present, the lr wfs are stored instead of the gs wfs
198	type(lr_t), optional, intent(in) :: lr
199
200	integer :: iunit, iunit2, iunit_mesh, iunit_states, err, ik, ist, idim, i
201	character(len=80) :: filename, mformat
202	logical :: wfns_are_associated, lr_wfns_are_associated
203
204	call push_sub('restart.restart_write')
205
206	wfns_are_associated = (associated(st%dpsi) .and. st%d%wfs_type == M_REAL) .or. &
207	(associated(st%zpsi) .and. st%d%wfs_type == M_CMPLX)
208	ASSERT(wfns_are_associated)
209
210	if(present(lr)) then
211	lr_wfns_are_associated = (associated(lr%ddl_psi) .and. st%d%wfs_type == M_REAL) .or. &
212	(associated(lr%zdl_psi) .and. st%d%wfs_type == M_CMPLX)
213	ASSERT(lr_wfns_are_associated)
214	endif
215
216	mformat = '(f12.8,a,f15.8,a,4(f12.8,a),i10.10,3(a,i8))'
217	ierr = 0
218
219	call block_signals()
220	if(mpi_grp_is_root(mpi_world)) then
221	call io_mkdir(dir, is_tmp=.true.)
222
223	iunit = io_open(trim(dir)//'/wfns', action='write', is_tmp=.true.)
224	write(iunit,'(a)') '# #kpoint #st #dim filename'
225	write(iunit,'(a)') '%Wavefunctions'
226
227	iunit2 = io_open(trim(dir)//'/occs', action='write', is_tmp=.true.)
228	write(iunit2,'(a)') '# occupations \| eigenvalue[a.u.] \| K-Points \| K-Weights \| filename \| ik \| ist \| idim'
229	write(iunit2,'(a)') '%Occupations_Eigenvalues_K-Points'
230
231	iunit_mesh = io_open(trim(dir)//'/mesh', action='write', is_tmp=.true.)
232	write(iunit_mesh,'(a)') '# This file contains the necessary information to generate the'
233	write(iunit_mesh,'(a)') '# mesh with which the functions in this directory were calculated,'
234	write(iunit_mesh,'(a)') '# except for the geometry of the system.'
235	call curvlinear_dump(gr%cv, iunit_mesh)
236	call simul_box_dump(gr%sb, iunit_mesh)
237	call mesh_dump(gr%m, iunit_mesh)
238	call io_close(iunit_mesh)
239
240	iunit_mesh = io_open(trim(dir)//'/Lxyz', action='write', is_tmp=.true.)
241	call mesh_Lxyz_dump(gr%m, iunit_mesh)
242	call io_close(iunit_mesh)
243
244	iunit_states = io_open(trim(dir)//'/states', action='write', is_tmp=.true.)
245	call states_dump(st, iunit_states)
246	call io_close(iunit_states)
247	end if
248
249	i = 1
250	do ik = 1, st%d%nik
251	do ist = 1, st%nst
252	do idim = 1, st%d%dim
253	write(filename,'(i10.10)') i
254
255	if(mpi_grp_is_root(mpi_world)) then
256	write(unit=iunit, fmt=*) ik, ' \| ', ist, ' \| ', idim, ' \| "', trim(filename), '"'
257	write(unit=iunit2, fmt=mformat) st%occ(ist,ik), ' \| ', st%eigenval(ist, ik), ' \| ', &
258	st%d%kpoints(1,ik), ' \| ', st%d%kpoints(2,ik), ' \| ', st%d%kpoints(3,ik) , ' \| ', &
259	st%d%kweights(ik), ' \| ', i, ' \| ', ik, ' \| ', ist, ' \| ', idim
260	end if
261
262	if(st%st_start <= ist .and. st%st_end >= ist) then
263	if( .not. present(lr) ) then
264	if (st%d%wfs_type == M_REAL) then
265	call drestart_write_function(dir, filename, gr, st%dpsi(:, idim, ist, ik), err, size(st%dpsi,1))
266	else
267	call zrestart_write_function(dir, filename, gr, st%zpsi(:, idim, ist, ik), err, size(st%zpsi,1))
268	end if
269	else
270	if (st%d%wfs_type == M_REAL) then
271	call drestart_write_function(dir, filename, gr, lr%ddl_psi(:, idim, ist, ik), err, size(st%dpsi,1))
272	else
273	call zrestart_write_function(dir, filename, gr, lr%zdl_psi(:, idim, ist, ik), err, size(st%zpsi,1))
274	end if
275	end if
276	if(err == 0) ierr = ierr + 1
277	end if
278	#if defined(HAVE_MPI)
279	call MPI_Barrier(MPI_COMM_WORLD, mpi_err) ! now we all wait
280	#endif
281	i = i + 1
282	end do
283	end do
284	end do
285
286	if(ierr == st%d%nik(st%st_end - st%st_start + 1)st%d%dim) ierr = 0 ! Alles OK
287
288	if(mpi_grp_is_root(mpi_world)) then
289	write(iunit,'(a)') '%'
290	if(present(iter)) write(iunit,'(a,i7)') 'Iter = ', iter
291	write(iunit2, '(a)') '%'
292
293	call io_close(iunit)
294	call io_close(iunit2)
295	end if
296
297	#if defined(HAVE_MPI)
298	call MPI_Barrier(MPI_COMM_WORLD, mpi_err) ! Since some processors did more than others...
299	#endif
300	call unblock_signals()
301
302	call pop_sub()
303	end subroutine restart_write
304
305
306	! ---------------------------------------------------------
307	! returns
308	! <0 => Fatal error
309	! =0 => read all wave-functions
310	! >0 => could only read x wavefunctions
311	subroutine restart_read(dir, st, gr, geo, ierr, iter, lr)
312	character(len=*), intent(in) :: dir
313	type(states_t), intent(inout) :: st
314	type(grid_t), intent(in) :: gr
315	type(geometry_t), intent(in) :: geo
316	integer, intent(out) :: ierr
317	integer, optional, intent(out) :: iter
318	!if this next argument is present, the lr wfs are read instead of the gs wfs
319	type(lr_t), optional, intent(inout) :: lr
320
321	integer :: iunit, iunit2, iunit_mesh, err, ik, ist, idim, i
322	character(len=12) :: filename
323	character(len=1) :: char
324	logical, allocatable :: filled(:, :, :)
325	character(len=256) :: line
326	character(len=50) :: str
327
328	FLOAT, allocatable :: dphi(:)
329	CMPLX, allocatable :: zphi(:)
330	type(mesh_t) :: old_mesh
331	type(curvlinear_t) :: old_cv
332	type(simul_box_t) :: old_sb
333	logical :: mesh_change, full_interpolation, boolean
334
335	call push_sub('restart.restart_read')
336
337	! sanity check
338	boolean = (associated(st%dpsi) .and. st%d%wfs_type == M_REAL) .or. &
339	(associated(st%zpsi) .and. st%d%wfs_type == M_CMPLX)
340	ASSERT(boolean)
341
342	if(present(lr)) then
343	boolean = (associated(lr%ddl_psi) .and. st%d%wfs_type == M_REAL) .or. &
344	(associated(lr%zdl_psi) .and. st%d%wfs_type == M_CMPLX)
345	ASSERT(boolean)
346	endif
347
348	ierr = 0
349
350	if(.not. present(lr)) then
351	write(str, '(a,i5)') 'Loading restart information'
352	else
353	write(str, '(a,i5)') 'Loading restart information for linear response'
354	end if
355
356	call messages_print_stress(stdout, trim(str))
357	! open files to read
358	iunit = io_open(trim(dir)//'/wfns', action='read', status='old', die=.false., is_tmp = .true., grp = gr%m%mpi_grp)
359	if(iunit < 0) then
360	ierr = -1
361	return
362	end if
363	iunit2 = io_open(trim(dir)//'/occs', action='read', status='old', die=.false., is_tmp = .true., grp = gr%m%mpi_grp)
364	if(iunit2 < 0) then
365	call io_close(iunit, grp = gr%m%mpi_grp)
366	ierr = -1
367	end if
368
369	if(ierr.ne.0) then
370	write(message(1),'(a)') 'Could not load any previous restart information.'
371	call write_info(1)
372	call messages_print_stress(stdout)
373	call pop_sub()
374	return
375	end if
376
377	! Reads out the previous mesh info.
378	call read_previous_mesh()
379
380	! now we really start
381	ALLOCATE(filled(st%d%dim, st%st_start:st%st_end, st%d%nik), st%d%dim(st%st_end-st%st_start+1)st%d%nik)
382	filled = .false.
383
384	! Skip two lines.
385	call iopar_read(gr%m%mpi_grp, iunit, line, err); call iopar_read(gr%m%mpi_grp, iunit, line, err)
386	call iopar_read(gr%m%mpi_grp, iunit2, line, err); call iopar_read(gr%m%mpi_grp, iunit2, line, err)
387
388	do
389	call iopar_read(gr%m%mpi_grp, iunit, line, i)
390	read(line, '(a)') char
391	if(i.ne.0.or.char=='%') exit
392
393	call iopar_backspace(gr%m%mpi_grp, iunit)
394
395	call iopar_read(gr%m%mpi_grp, iunit, line, err)
396	read(line, *) ik, char, ist, char, idim, char, filename
397	if(index_is_wrong()) then
398	call iopar_read(gr%m%mpi_grp, iunit2, line, err)
399	cycle
400	end if
401
402	call iopar_read(gr%m%mpi_grp, iunit2, line, err)
403	read(line, *) st%occ(ist, ik), char, st%eigenval(ist, ik)
404
405	if(ist >= st%st_start .and. ist <= st%st_end) then
406	if(.not.mesh_change) then
407	if( .not. present(lr) ) then
408	if (st%d%wfs_type == M_REAL) then
409	call drestart_read_function(dir, filename, gr%m, st%dpsi(:, idim, ist, ik), err)
410	else
411	call zrestart_read_function(dir, filename, gr%m, st%zpsi(:, idim, ist, ik), err)
412	end if
413	else
414	if (st%d%wfs_type == M_REAL) then
415	call drestart_read_function(dir, filename, gr%m, lr%ddl_psi(:, idim, ist, ik), err)
416	else
417	call zrestart_read_function(dir, filename, gr%m, lr%zdl_psi(:, idim, ist, ik), err)
418	end if
419	end if
420	if(err <= 0) then
421	filled(idim, ist, ik) = .true.
422	ierr = ierr + 1
423	end if
424	else
425	if (st%d%wfs_type == M_REAL) then
426	call drestart_read_function(dir, filename, old_mesh, dphi, err)
427	if( .not. present(lr) ) then
428	call dmf_interpolate(old_mesh, gr%m, full_interpolation, dphi, st%dpsi(:, idim, ist, ik))
429	else
430	call dmf_interpolate(old_mesh, gr%m, full_interpolation, dphi, lr%ddl_psi(:, idim, ist, ik))
431	end if
432	else
433	call zrestart_read_function(dir, filename, old_mesh, zphi, err)
434	if( .not. present(lr) ) then
435	call zmf_interpolate(old_mesh, gr%m, full_interpolation, zphi, st%zpsi(:, idim, ist, ik))
436	else
437	call zmf_interpolate(old_mesh, gr%m, full_interpolation, zphi, lr%zdl_psi(:, idim, ist, ik))
438	end if
439	end if
440	if(err <= 0) then
441	filled(idim, ist, ik) = .true.
442	ierr = ierr + 1
443	end if
444	endif
445	end if
446	end do
447
448	if(present(iter)) then
449	call iopar_read(gr%m%mpi_grp, iunit, line, err)
450	read(line, *) filename, filename, iter
451	end if
452
453	#if defined(HAVE_MPI)
454	if(st%parallel_in_states) then
455	call MPI_Allreduce(ierr, err, 1, MPI_INTEGER, MPI_SUM, st%mpi_grp%comm, mpi_err)
456	ierr = err
457	end if
458	#endif
459
460	call fill()
461	if(ierr == 0) then
462	ierr = -1 ! no files read
463	write(message(1),'(a)') 'No files could be read. No restart information can be used.'
464	call write_info(1)
465	else
466	! Everything o. k.
467	if(ierr == st%nstst%d%nikst%d%dim) then
468	ierr = 0
469	write(message(1),'(a)') 'All the needed files were succesfully read.'
470	call write_info(1)
471	else
472	write(message(1),'(a,i4,a,i4,a)') 'Only ', ierr,' files out of ', &
473	st%nstst%d%nikst%d%dim, ' could be read.'
474	call write_info(1)
475	endif
476	end if
477
478	deallocate(filled)
479	call io_close(iunit, grp = gr%m%mpi_grp)
480	call io_close(iunit2, grp = gr%m%mpi_grp)
481
482	if(mesh_change) call interpolation_end()
483
484	call messages_print_stress(stdout)
485	call pop_sub()
486
487	contains
488
489	! ---------------------------------------------------------
490	subroutine interpolation_init
491	call mesh_init_stage_1(old_sb, old_mesh, geo, old_cv, gr%f_der%n_ghost)
492	call mesh_init_stage_2(old_sb, old_mesh, geo, old_cv)
493	call mesh_init_stage_3(old_mesh, geo, old_cv)
494
495	if (st%d%wfs_type == M_REAL) then
496	ALLOCATE(dphi(old_mesh%np_global), old_mesh%np_global)
497	else
498	ALLOCATE(zphi(old_mesh%np_global), old_mesh%np_global)
499	end if
500	end subroutine interpolation_init
501
502
503	! ---------------------------------------------------------
504	subroutine interpolation_end
505	call mesh_end(old_mesh)
506	if (st%d%wfs_type == M_REAL) then
507	deallocate(dphi)
508	else
509	deallocate(zphi)
510	end if
511	end subroutine interpolation_end
512
513
514	! ---------------------------------------------------------
515	subroutine read_previous_mesh
516
517	mesh_change = .false.
518	full_interpolation = .true.
519
520	#if defined(HAVE_MPI)
521	return ! For the moment, avoid the complications of parallel stuff.
522	#endif
523	if(present(iter)) return ! No intepolation, in case we are in the td part.
524
525	iunit_mesh = io_open(trim(dir)//'/mesh', action='read', status='old', die=.false., grp = gr%m%mpi_grp)
526	if(iunit_mesh < 0) return
527
528	read(iunit_mesh, ); read(iunit_mesh, ); read(iunit_mesh, *)
529	call curvlinear_init_from_file(old_cv, iunit_mesh)
530	call simul_box_init_from_file(old_sb, iunit_mesh)
531	call io_close(iunit_mesh, grp = gr%m%mpi_grp)
532
533	if( .not. (old_cv.eq.gr%cv) ) then
534	mesh_change = .true.
535	return
536	end if
537	if( .not. (old_sb.eq.gr%sb) ) then
538	mesh_change = .true.
539	! First, check wether the spacings are the same.
540	if(old_sb%h .app. gr%sb%h) full_interpolation = .false.
541	end if
542
543	if(mesh_change) then
544	if(.not.full_interpolation) then
545	write(message(1),'(a)') 'The functions stored in "tmp/restart_gs" were obtained with'
546	write(message(2),'(a)') 'a different simulation box. The possible missing regions will be'
547	write(message(3),'(a)') 'padded with zeros.'
548	call write_info(3)
549	else
550	write(message(1),'(a)') 'The functions stored in "tmp/restart_gs" were obtained with'
551	write(message(2),'(a)') 'a different mesh. The values of the functions for the current'
552	write(message(3),'(a)') 'calculations will be interpolated/extrapolated.'
553	call write_info(3)
554	end if
555	call interpolation_init()
556	endif
557
558	end subroutine read_previous_mesh
559
560
561	! ---------------------------------------------------------
562	subroutine fill() ! Put random function in orbitals that could not be read.
563	do ik = 1, st%d%nik
564	do ist = st%st_start, st%st_end
565	do idim = 1, st%d%dim
566	if(filled(idim, ist, ik)) cycle
567	write(message(1),'(a,3i5)') 'Randomizing wavefunction: #dim, #ist, #ik = ', idim, ist, ik
568	call write_warning(1)
569
570	call states_generate_random(st, gr%m, ist, ist)
571	st%occ(ist, ik) = M_ZERO
572	end do
573	end do
574	end do
575	end subroutine fill
576
577
578	! ---------------------------------------------------------
579	logical function index_is_wrong() ! .true. if the index (idim, ist, ik) is not present in st structure...
580	if(idim > st%d%dim .or. idim < 1 .or. &
581	ist > st%nst .or. ist < 1 .or. &
582	ik > st%d%nik .or. ik < 1) then
583	index_is_wrong = .true.
584	else
585	index_is_wrong = .false.
586	end if
587	end function index_is_wrong
588
589	end subroutine restart_read
590
591	#include "undef.F90"
592	#include "real.F90"
593	#include "restart_inc.F90"
594
595	#include "undef.F90"
596	#include "complex.F90"
597	#include "restart_inc.F90"
598
599	end module restart_m
600
601	!! Local Variables:
602	!! mode: f90
603	!! coding: utf-8
604	!! End:

Note: See TracBrowser for help on using the browser.

Download in other formats:

Original Format