Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief - writing and printing the files, trajectory (pos, cell, dipoles) as
10 : !> well as restart files
11 : !> - usually just the Markov Chain elements are regarded, the elements
12 : !> beside this trajectory are neglected
13 : !> - futrthermore (by option) just the accepted configurations
14 : !> are print out to reduce the file sizes
15 : !> \par History
16 : !> 12.2012 created [Mandes Schoenherr]
17 : !> \author Mandes
18 : ! **************************************************************************************************
19 :
20 : MODULE tmc_file_io
21 : USE cp_files, ONLY: close_file,&
22 : open_file
23 : USE cp_log_handling, ONLY: cp_to_string
24 : USE kinds, ONLY: default_path_length,&
25 : default_string_length,&
26 : dp
27 : USE physcon, ONLY: au2a => angstrom
28 : USE tmc_analysis_types, ONLY: tmc_analysis_env
29 : USE tmc_calculations, ONLY: get_cell_scaling,&
30 : get_scaled_cell
31 : USE tmc_move_types, ONLY: nr_mv_types
32 : USE tmc_stati, ONLY: TMC_STATUS_FAILED,&
33 : TMC_STATUS_OK,&
34 : TMC_STATUS_WAIT_FOR_NEW_TASK,&
35 : tmc_default_restart_in_file_name,&
36 : tmc_default_restart_out_file_name,&
37 : tmc_default_trajectory_file_name
38 : USE tmc_tree_types, ONLY: elem_array_type,&
39 : tree_type
40 : USE tmc_types, ONLY: tmc_env_type,&
41 : tmc_param_type
42 : #include "../base/base_uses.f90"
43 :
44 : IMPLICIT NONE
45 :
46 : PRIVATE
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_file_io'
49 :
50 : ! filename manipulation
51 : PUBLIC :: expand_file_name_char, expand_file_name_temp, expand_file_name_int
52 : ! read/write restart file
53 : PUBLIC :: print_restart_file, read_restart_file
54 : ! write the configuration
55 : PUBLIC :: write_result_list_element
56 : PUBLIC :: write_element_in_file
57 : PUBLIC :: write_dipoles_in_file
58 : ! analysis read
59 : PUBLIC :: analyse_files_open, read_element_from_file, analyse_files_close
60 :
61 : CONTAINS
62 :
63 : !------------------------------------------------------------------------------
64 : ! routines for manipulating the file name
65 : !------------------------------------------------------------------------------
66 : ! **************************************************************************************************
67 : !> \brief placing a character string at the end of a file name
68 : !> (instead of the ending)
69 : !> \param file_name original file name
70 : !> \param extra string to be added before the file extension
71 : !> \return the new filename
72 : !> \author Mandes 11.2012
73 : ! **************************************************************************************************
74 2614 : FUNCTION expand_file_name_ending(file_name, extra) RESULT(result_file_name)
75 : CHARACTER(LEN=*) :: file_name, extra
76 : CHARACTER(LEN=default_path_length) :: result_file_name
77 :
78 : INTEGER :: ind
79 :
80 0 : CPASSERT(file_name .NE. "")
81 :
82 2614 : ind = INDEX(file_name, ".", BACK=.TRUE.)
83 2614 : IF (.NOT. ind .EQ. 0) THEN
84 2614 : WRITE (result_file_name, *) file_name(1:ind - 1), ".", &
85 5228 : TRIM(ADJUSTL(extra))
86 : ELSE
87 0 : WRITE (result_file_name, *) TRIM(file_name), ".", extra
88 : END IF
89 2614 : result_file_name = TRIM(ADJUSTL(result_file_name))
90 2614 : CPASSERT(result_file_name .NE. "")
91 2614 : END FUNCTION expand_file_name_ending
92 :
93 : ! **************************************************************************************************
94 : !> \brief placing a character string at the end of a file name
95 : !> (before the file extension)
96 : !> \param file_name original file name
97 : !> \param extra string to be added before the file extension
98 : !> \return the new filename
99 : !> \author Mandes 11.2012
100 : ! **************************************************************************************************
101 1427 : FUNCTION expand_file_name_char(file_name, extra) RESULT(result_file_name)
102 : CHARACTER(LEN=*) :: file_name, extra
103 : CHARACTER(LEN=default_path_length) :: result_file_name
104 :
105 : INTEGER :: ind
106 :
107 0 : CPASSERT(file_name .NE. "")
108 :
109 1427 : ind = INDEX(file_name, ".", BACK=.TRUE.)
110 1427 : IF (.NOT. ind .EQ. 0) THEN
111 1427 : WRITE (result_file_name, *) file_name(1:ind - 1), "_", &
112 2854 : TRIM(ADJUSTL(extra)), file_name(ind:LEN_TRIM(file_name))
113 : ELSE
114 0 : WRITE (result_file_name, *) TRIM(file_name), "_", extra
115 : END IF
116 1427 : result_file_name = TRIM(ADJUSTL(result_file_name))
117 1427 : CPASSERT(result_file_name .NE. "")
118 1427 : END FUNCTION expand_file_name_char
119 :
120 : ! **************************************************************************************************
121 : !> \brief placing the temperature at the end of a file name
122 : !> (before the file extension)
123 : !> \param file_name original file name
124 : !> \param rvalue temperature to be added
125 : !> \return the new filename
126 : !> \author Mandes 11.2012
127 : ! **************************************************************************************************
128 2870 : FUNCTION expand_file_name_temp(file_name, rvalue) RESULT(result_file_name)
129 : CHARACTER(LEN=*) :: file_name
130 : REAL(KIND=dp) :: rvalue
131 : CHARACTER(LEN=default_path_length) :: result_file_name
132 :
133 : CHARACTER(LEN=18) :: rval_to_string
134 : INTEGER :: ind
135 :
136 2870 : CPASSERT(file_name .NE. "")
137 :
138 2870 : rval_to_string = ""
139 :
140 2870 : WRITE (rval_to_string, "(F16.2)") rvalue
141 2870 : ind = INDEX(file_name, ".", BACK=.TRUE.)
142 2870 : IF (.NOT. ind .EQ. 0) THEN
143 2870 : WRITE (result_file_name, *) file_name(1:ind - 1), "_T", &
144 5740 : TRIM(ADJUSTL(rval_to_string)), file_name(ind:LEN_TRIM(file_name))
145 : ELSE
146 0 : IF (LEN(file_name) .EQ. 0) THEN
147 0 : WRITE (result_file_name, *) TRIM(file_name), "T", TRIM(ADJUSTL(rval_to_string)), &
148 0 : file_name(ind:LEN_TRIM(file_name))
149 : ELSE
150 0 : WRITE (result_file_name, *) TRIM(file_name), "_T", TRIM(ADJUSTL(rval_to_string))
151 : END IF
152 : END IF
153 2870 : result_file_name = TRIM(ADJUSTL(result_file_name))
154 2870 : CPASSERT(result_file_name .NE. "")
155 2870 : END FUNCTION expand_file_name_temp
156 :
157 : ! **************************************************************************************************
158 : !> \brief placing an integer at the end of a file name
159 : !> (before the file extension)
160 : !> \param file_name original file name
161 : !> \param ivalue number to be added
162 : !> \return the new filename
163 : !> \author Mandes 11.2012
164 : ! **************************************************************************************************
165 19 : FUNCTION expand_file_name_int(file_name, ivalue) RESULT(result_file_name)
166 : CHARACTER(LEN=*) :: file_name
167 : INTEGER :: ivalue
168 : CHARACTER(LEN=default_path_length) :: result_file_name
169 :
170 : CHARACTER(LEN=18) :: rval_to_string
171 : INTEGER :: ind
172 :
173 19 : CPASSERT(file_name .NE. "")
174 :
175 19 : rval_to_string = ""
176 :
177 19 : WRITE (rval_to_string, *) ivalue
178 19 : ind = INDEX(file_name, ".", BACK=.TRUE.)
179 19 : IF (.NOT. ind .EQ. 0) THEN
180 19 : WRITE (result_file_name, *) file_name(1:ind - 1), "_", &
181 38 : TRIM(ADJUSTL(rval_to_string)), file_name(ind:LEN_TRIM(file_name))
182 : ELSE
183 0 : IF (LEN(file_name) .EQ. 0) THEN
184 0 : WRITE (result_file_name, *) TRIM(file_name), "", TRIM(ADJUSTL(rval_to_string)), &
185 0 : file_name(ind:LEN_TRIM(file_name))
186 : ELSE
187 0 : WRITE (result_file_name, *) TRIM(file_name), "_", TRIM(ADJUSTL(rval_to_string)), &
188 0 : file_name(ind:LEN_TRIM(file_name))
189 : END IF
190 : END IF
191 19 : result_file_name = TRIM(ADJUSTL(result_file_name))
192 19 : CPASSERT(result_file_name .NE. "")
193 19 : END FUNCTION expand_file_name_int
194 :
195 : !------------------------------------------------------------------------------
196 : ! routines for reading and writing RESTART file
197 : !------------------------------------------------------------------------------
198 : ! **************************************************************************************************
199 : !> \brief prints out the TMC restart files with all last configurations and
200 : !> counters etc.
201 : !> \param tmc_env the tmc environment, storing result lists and counters an in
202 : !> temperatures
203 : !> \param job_counts the counters for counting the submitted different job types
204 : !> \param timings ...
205 : !> \author Mandes 11.2012
206 : ! **************************************************************************************************
207 3 : SUBROUTINE print_restart_file(tmc_env, job_counts, timings)
208 : TYPE(tmc_env_type), POINTER :: tmc_env
209 : INTEGER, DIMENSION(:) :: job_counts
210 : REAL(KIND=dp), DIMENSION(4) :: timings
211 :
212 : CHARACTER(LEN=default_path_length) :: c_tmp, file_name
213 : INTEGER :: f_unit, i
214 :
215 3 : c_tmp = ""
216 3 : CPASSERT(ASSOCIATED(tmc_env))
217 3 : CPASSERT(ASSOCIATED(tmc_env%m_env))
218 3 : CPASSERT(ASSOCIATED(tmc_env%params))
219 3 : CPASSERT(ASSOCIATED(tmc_env%m_env%gt_act))
220 :
221 3 : WRITE (c_tmp, FMT='(I9.9)') tmc_env%m_env%result_count(0)
222 : file_name = TRIM(expand_file_name_char( &
223 : file_name=tmc_default_restart_out_file_name, &
224 3 : extra=c_tmp))
225 : CALL open_file(file_name=file_name, file_status="REPLACE", &
226 : file_action="WRITE", file_form="UNFORMATTED", &
227 3 : unit_number=f_unit)
228 3 : WRITE (f_unit) SIZE(tmc_env%params%Temp)
229 3 : WRITE (f_unit) tmc_env%params%Temp(:), &
230 3 : tmc_env%m_env%gt_act%nr, &
231 3 : tmc_env%m_env%gt_act%rng_seed, &
232 3 : tmc_env%m_env%gt_act%rnd_nr, &
233 3 : tmc_env%m_env%gt_act%prob_acc, &
234 3 : tmc_env%m_env%gt_act%mv_conf, &
235 3 : tmc_env%m_env%gt_act%mv_next_conf, &
236 3 : tmc_env%m_env%result_count(0:), &
237 3 : tmc_env%params%move_types%mv_weight, &
238 3 : tmc_env%params%move_types%acc_count, &
239 3 : tmc_env%params%move_types%mv_count, &
240 3 : tmc_env%params%move_types%subbox_acc_count, &
241 3 : tmc_env%params%move_types%subbox_count, &
242 3 : tmc_env%params%cell%hmat, &
243 3 : job_counts, &
244 6 : timings
245 12 : DO i = 1, SIZE(tmc_env%params%Temp)
246 9 : WRITE (f_unit) tmc_env%m_env%result_list(i)%elem%nr, &
247 252 : tmc_env%m_env%result_list(i)%elem%rng_seed, &
248 576 : tmc_env%m_env%result_list(i)%elem%pos, &
249 576 : tmc_env%m_env%result_list(i)%elem%vel, &
250 36 : tmc_env%m_env%result_list(i)%elem%box_scale, &
251 9 : tmc_env%m_env%result_list(i)%elem%potential, &
252 9 : tmc_env%m_env%result_list(i)%elem%e_pot_approx, &
253 9 : tmc_env%m_env%result_list(i)%elem%ekin, &
254 9 : tmc_env%m_env%result_list(i)%elem%ekin_before_md, &
255 21 : tmc_env%m_env%result_list(i)%elem%temp_created
256 : END DO
257 3 : CALL close_file(unit_number=f_unit)
258 : ! write the file, where the restart file name is written in
259 : CALL open_file(file_name=tmc_default_restart_in_file_name, &
260 : file_action="WRITE", file_status="REPLACE", &
261 3 : unit_number=f_unit)
262 3 : WRITE (f_unit, *) TRIM(file_name)
263 3 : CALL close_file(unit_number=f_unit)
264 3 : END SUBROUTINE print_restart_file
265 :
266 : ! **************************************************************************************************
267 : !> \brief reads the TMC restart file with all last configurations and
268 : !> counters etc.
269 : !> \param tmc_env the tmc environment, storing result lists and counters an in
270 : !> temperatures
271 : !> \param job_counts the counters for counting the submitted different job types
272 : !> \param timings ...
273 : !> \param file_name the restart file name
274 : !> \author Mandes 11.2012
275 : ! **************************************************************************************************
276 2 : SUBROUTINE read_restart_file(tmc_env, job_counts, timings, file_name)
277 : TYPE(tmc_env_type), POINTER :: tmc_env
278 : INTEGER, DIMENSION(:) :: job_counts
279 : REAL(KIND=dp), DIMENSION(4) :: timings
280 : CHARACTER(LEN=*) :: file_name
281 :
282 : INTEGER :: file_ptr, i, temp_size
283 : LOGICAL :: flag
284 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp_temp
285 : REAL(KIND=dp), DIMENSION(nr_mv_types) :: mv_weight_tmp
286 :
287 2 : CPASSERT(ASSOCIATED(tmc_env))
288 2 : CPASSERT(ASSOCIATED(tmc_env%m_env))
289 2 : CPASSERT(ASSOCIATED(tmc_env%params))
290 2 : CPASSERT(ASSOCIATED(tmc_env%m_env%gt_act))
291 :
292 2 : IF (file_name .EQ. tmc_default_restart_in_file_name) THEN
293 2 : INQUIRE (FILE=tmc_default_restart_in_file_name, EXIST=flag)
294 2 : CPASSERT(flag)
295 : CALL open_file(file_name=tmc_default_restart_in_file_name, file_status="OLD", &
296 2 : file_action="READ", unit_number=file_ptr)
297 2 : READ (file_ptr, *) file_name
298 2 : CALL close_file(unit_number=file_ptr)
299 : END IF
300 :
301 : CALL open_file(file_name=file_name, file_status="OLD", file_form="UNFORMATTED", &
302 2 : file_action="READ", unit_number=file_ptr)
303 2 : READ (file_ptr) temp_size
304 2 : IF (temp_size .NE. SIZE(tmc_env%params%Temp)) &
305 : CALL cp_abort(__LOCATION__, &
306 : "the actual specified temperatures does not "// &
307 0 : "fit in amount with the one from restart file ")
308 6 : ALLOCATE (tmp_temp(temp_size))
309 2 : READ (file_ptr) tmp_temp(:), &
310 2 : tmc_env%m_env%gt_act%nr, &
311 2 : tmc_env%m_env%gt_act%rng_seed, &
312 2 : tmc_env%m_env%gt_act%rnd_nr, &
313 2 : tmc_env%m_env%gt_act%prob_acc, &
314 2 : tmc_env%m_env%gt_act%mv_conf, & !
315 2 : tmc_env%m_env%gt_act%mv_next_conf, & !
316 2 : tmc_env%m_env%result_count(0:), &
317 2 : mv_weight_tmp, & !
318 2 : tmc_env%params%move_types%acc_count, &
319 2 : tmc_env%params%move_types%mv_count, &
320 2 : tmc_env%params%move_types%subbox_acc_count, &
321 2 : tmc_env%params%move_types%subbox_count, & !
322 2 : tmc_env%params%cell%hmat, &
323 2 : job_counts, &
324 4 : timings
325 :
326 8 : IF (ANY(ABS(tmc_env%params%Temp(:) - tmp_temp(:)) .GE. 0.005)) &
327 : CALL cp_abort(__LOCATION__, "the temperatures differ from the previous calculation. "// &
328 0 : "There were the following temperatures used:")
329 22 : IF (ANY(mv_weight_tmp(:) .NE. tmc_env%params%move_types%mv_weight(:))) &
330 0 : CPWARN("The amount of mv types differs between the original and the restart run.")
331 :
332 8 : DO i = 1, SIZE(tmc_env%params%Temp)
333 6 : tmc_env%m_env%gt_act%conf(i)%elem => tmc_env%m_env%result_list(i)%elem
334 6 : READ (file_ptr) tmc_env%m_env%result_list(i)%elem%nr, &
335 168 : tmc_env%m_env%result_list(i)%elem%rng_seed, &
336 384 : tmc_env%m_env%result_list(i)%elem%pos, &
337 384 : tmc_env%m_env%result_list(i)%elem%vel, &
338 24 : tmc_env%m_env%result_list(i)%elem%box_scale, &
339 6 : tmc_env%m_env%result_list(i)%elem%potential, &
340 6 : tmc_env%m_env%result_list(i)%elem%e_pot_approx, &
341 6 : tmc_env%m_env%result_list(i)%elem%ekin, &
342 6 : tmc_env%m_env%result_list(i)%elem%ekin_before_md, &
343 14 : tmc_env%m_env%result_list(i)%elem%temp_created
344 : END DO
345 2 : CALL close_file(unit_number=file_ptr)
346 2 : END SUBROUTINE read_restart_file
347 :
348 : !----------------------------------------------------------------------------
349 : ! printing configuration in file
350 : !----------------------------------------------------------------------------
351 :
352 : ! **************************************************************************************************
353 : !> \brief select the correct configuration to print out the
354 : !> (coordinates, forces, cell ...)
355 : !> \param result_list list of configurations for each temperature
356 : !> \param result_count list with number of Markov Chain number
357 : !> for each teperature (index 0 for global tree)
358 : !> \param conf_updated index of the updated (modified element)
359 : !> \param accepted acceptance flag
360 : !> \param tmc_params TMC environment parameters
361 : !> \author Mandes 02.2013
362 : ! **************************************************************************************************
363 9296 : SUBROUTINE write_result_list_element(result_list, result_count, conf_updated, &
364 : accepted, tmc_params)
365 : TYPE(elem_array_type), DIMENSION(:), POINTER :: result_list
366 : INTEGER, DIMENSION(:), POINTER :: result_count
367 : INTEGER :: conf_updated
368 : LOGICAL, INTENT(IN) :: accepted
369 : TYPE(tmc_param_type), POINTER :: tmc_params
370 :
371 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_result_list_element'
372 :
373 : CHARACTER(LEN=default_path_length) :: file_name
374 : INTEGER :: handle, i
375 :
376 4648 : file_name = ""
377 :
378 4648 : CPASSERT(ASSOCIATED(result_list))
379 4648 : CPASSERT(ASSOCIATED(result_count))
380 4648 : CPASSERT(ASSOCIATED(tmc_params))
381 4648 : CPASSERT(ASSOCIATED(tmc_params%Temp))
382 4648 : CPASSERT(conf_updated .GE. 0)
383 4648 : CPASSERT(conf_updated .LE. SIZE(tmc_params%Temp))
384 :
385 : ! start the timing
386 4648 : CALL timeset(routineN, handle)
387 :
388 4648 : IF (conf_updated .EQ. 0) THEN
389 : ! for debugging print every configuration of every temperature
390 0 : DO i = 1, SIZE(tmc_params%Temp)
391 0 : WRITE (file_name, *) "every_step_", TRIM(tmc_default_trajectory_file_name)
392 : CALL write_element_in_file(elem=result_list(i)%elem, &
393 : tmc_params=tmc_params, conf_nr=result_count(0), &
394 0 : file_name=expand_file_name_temp(file_name=file_name, rvalue=tmc_params%Temp(i)))
395 : END DO
396 : ELSE
397 4648 : IF ((.NOT. tmc_params%print_only_diff_conf) .OR. &
398 : (tmc_params%print_only_diff_conf .AND. accepted)) THEN
399 : CALL write_element_in_file(elem=result_list(conf_updated)%elem, &
400 : tmc_params=tmc_params, conf_nr=result_count(conf_updated), &
401 : file_name=expand_file_name_temp(file_name=TRIM(tmc_default_trajectory_file_name), &
402 1024 : rvalue=tmc_params%Temp(conf_updated)))
403 : END IF
404 : END IF
405 : ! end the timing
406 4648 : CALL timestop(handle)
407 4648 : END SUBROUTINE write_result_list_element
408 :
409 : ! **************************************************************************************************
410 : !> \brief writes the trajectory element in a file from sub tree element
411 : !> \param elem actual tree element to be printed out
412 : !> \param tmc_params TMC environment parameters
413 : !> \param temp_index ...
414 : !> \param file_name file name will be extended by type of file (pos, cell,...)
415 : !> \param conf_nr Markov chain element number
416 : !> \param conf_info whole header line
417 : !> \author Mandes 11.2012
418 : ! **************************************************************************************************
419 1024 : SUBROUTINE write_element_in_file(elem, tmc_params, temp_index, file_name, conf_nr, &
420 : conf_info)
421 : TYPE(tree_type), POINTER :: elem
422 : TYPE(tmc_param_type), POINTER :: tmc_params
423 : INTEGER, OPTIONAL :: temp_index
424 : CHARACTER(LEN=*), OPTIONAL :: file_name
425 : INTEGER, OPTIONAL :: conf_nr
426 : CHARACTER(LEN=*), OPTIONAL :: conf_info
427 :
428 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_element_in_file'
429 :
430 : CHARACTER(LEN=default_path_length) :: file_name_act, tmp_name
431 : CHARACTER(LEN=default_string_length) :: header
432 : INTEGER :: file_ptr, handle, i, nr_atoms
433 : LOGICAL :: file_exists, print_it
434 : REAL(KIND=dp) :: vol
435 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat_scaled
436 :
437 1024 : file_name_act = ""
438 1024 : tmp_name = ""
439 1024 : header = ""
440 1024 : print_it = .TRUE.
441 :
442 0 : CPASSERT(ASSOCIATED(elem))
443 1024 : CPASSERT(ASSOCIATED(tmc_params))
444 1024 : CPASSERT(ASSOCIATED(tmc_params%atoms))
445 1024 : CPASSERT(PRESENT(conf_nr) .OR. PRESENT(conf_info))
446 :
447 1024 : IF (print_it) THEN
448 : ! start the timing
449 1024 : CALL timeset(routineN, handle)
450 :
451 : ! set default file name
452 1024 : IF (PRESENT(file_name)) THEN
453 1024 : CPASSERT(file_name .NE. "")
454 1024 : file_name_act = file_name
455 : ELSE
456 0 : CPASSERT(ASSOCIATED(tmc_params%Temp))
457 0 : CPASSERT(PRESENT(temp_index))
458 : file_name_act = expand_file_name_temp(file_name=tmc_default_trajectory_file_name, &
459 0 : rvalue=tmc_params%Temp(temp_index))
460 : END IF
461 :
462 1024 : nr_atoms = SIZE(elem%pos)/tmc_params%dim_per_elem
463 :
464 : ! set header (for coordinate or force file)
465 1024 : IF (tmc_params%print_trajectory .OR. tmc_params%print_forces) THEN
466 1024 : IF (PRESENT(conf_info)) THEN
467 0 : WRITE (header, *) TRIM(ADJUSTL(conf_info))
468 : ELSE
469 : !WRITE(header,FMT="(A,I8,A,F20.10)") " i = ", conf_nr,", E = ", elem%potential
470 1024 : WRITE (header, FMT="(A,I8,A,F20.10,F20.10,A,I8,I8)") "i =", conf_nr, " ,E =", &
471 2048 : elem%potential, elem%ekin, " st elem", elem%sub_tree_nr, elem%nr
472 : END IF
473 : END IF
474 :
475 : ! write the coordinates
476 1024 : IF (tmc_params%print_trajectory) THEN
477 1024 : tmp_name = expand_file_name_ending(file_name_act, "xyz")
478 : CALL open_file(file_name=tmp_name, file_status="UNKNOWN", &
479 : file_action="WRITE", file_position="APPEND", &
480 1024 : unit_number=file_ptr)
481 1024 : WRITE (file_ptr, FMT="(I8)") nr_atoms
482 1024 : WRITE (file_ptr, *) TRIM(header)
483 44775 : DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
484 : WRITE (file_ptr, FMT="(A4,1X,1000F20.10)") &
485 43751 : TRIM(tmc_params%atoms((i - 1)/tmc_params%dim_per_elem + 1)%name), &
486 219779 : elem%pos(i:i + tmc_params%dim_per_elem - 1)*au2a
487 : END DO
488 1024 : CALL close_file(unit_number=file_ptr)
489 : END IF
490 :
491 : ! write the forces
492 1024 : IF (tmc_params%print_forces) THEN
493 331 : tmp_name = expand_file_name_ending(file_name_act, "frc")
494 : CALL open_file(file_name=tmp_name, file_status="UNKNOWN", &
495 : file_action="WRITE", file_position="APPEND", &
496 331 : unit_number=file_ptr)
497 331 : WRITE (file_ptr, FMT="(I8)") nr_atoms
498 331 : WRITE (file_ptr, *) TRIM(header)
499 7282 : DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
500 : WRITE (file_ptr, FMT="(A4,1X,1000F20.10)") &
501 6951 : TRIM(tmc_params%atoms((i - 1)/tmc_params%dim_per_elem + 1)%name), &
502 14233 : elem%frc(i:i + tmc_params%dim_per_elem - 1)
503 : END DO
504 331 : CALL close_file(unit_number=file_ptr)
505 : END IF
506 :
507 : ! write the cell dipoles
508 1024 : IF (tmc_params%print_dipole) THEN
509 : CALL write_dipoles_in_file(file_name=file_name_act, &
510 0 : conf_nr=conf_nr, dip=elem%dipole)
511 : END IF
512 :
513 : ! write the cell file
514 1024 : IF (tmc_params%print_cell) THEN
515 392 : tmp_name = expand_file_name_ending(file_name_act, "cell")
516 : ! header
517 392 : INQUIRE (FILE=tmp_name, EXIST=file_exists) ! file_exists will be TRUE if the file exist
518 392 : IF (.NOT. file_exists) THEN
519 : CALL open_file(file_name=tmp_name, file_status="NEW", &
520 6 : file_action="WRITE", unit_number=file_ptr)
521 : WRITE (file_ptr, FMT='(A,9(7X,A2," [Angstrom]"),6X,A)') &
522 6 : "# MC step ", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
523 12 : "Volume [Angstrom^3]"
524 : ELSE
525 : CALL open_file(file_name=tmp_name, file_status="OLD", &
526 : file_action="WRITE", file_position="APPEND", &
527 386 : unit_number=file_ptr)
528 : END IF
529 : CALL get_scaled_cell(cell=tmc_params%cell, &
530 : box_scale=elem%box_scale, scaled_hmat=hmat_scaled, &
531 392 : vol=vol)
532 392 : WRITE (file_ptr, FMT="(I8,9(1X,F19.10),1X,F24.10)") conf_nr, &
533 5488 : hmat_scaled(:, :)*au2a, vol*au2a**3
534 : !TODO better cell output e.g. using cell_types routine
535 392 : CALL close_file(unit_number=file_ptr)
536 : END IF
537 :
538 : ! write the different energies
539 1024 : IF (tmc_params%print_energies) THEN
540 331 : tmp_name = expand_file_name_ending(file_name_act, "ener")
541 : ! header
542 331 : INQUIRE (FILE=tmp_name, EXIST=file_exists) ! file_exists will be TRUE if the file exist
543 331 : IF (.NOT. file_exists) THEN
544 : CALL open_file(file_name=tmp_name, file_status="NEW", &
545 3 : file_action="WRITE", unit_number=file_ptr)
546 : WRITE (file_ptr, FMT='(A,4A20)') &
547 3 : "# MC step ", " exact ", " approx ", " last SCF ", " kinetic "
548 : ELSE
549 : CALL open_file(file_name=tmp_name, file_status="OLD", &
550 : file_action="WRITE", file_position="APPEND", &
551 328 : unit_number=file_ptr)
552 : END IF
553 331 : WRITE (file_ptr, FMT="(I8,14F20.10)") conf_nr, elem%potential, elem%e_pot_approx, &
554 662 : elem%scf_energies(MOD(elem%scf_energies_count, 4) + 1), elem%ekin
555 331 : CALL close_file(unit_number=file_ptr)
556 : END IF
557 :
558 : ! end the timing
559 1024 : CALL timestop(handle)
560 : END IF
561 1024 : END SUBROUTINE write_element_in_file
562 :
563 : ! **************************************************************************************************
564 : !> \brief writes the cell dipoles in dipole trajectory file
565 : !> \param file_name ...
566 : !> \param conf_nr ...
567 : !> \param dip ...
568 : !> \param file_ext ...
569 : !> \param
570 : !> \author Mandes 11.2012
571 : ! **************************************************************************************************
572 500 : SUBROUTINE write_dipoles_in_file(file_name, conf_nr, dip, file_ext)
573 : CHARACTER(LEN=default_path_length) :: file_name
574 : INTEGER :: conf_nr
575 : REAL(KIND=dp), DIMENSION(:), POINTER :: dip
576 : CHARACTER(LEN=*), INTENT(in), OPTIONAL :: file_ext
577 :
578 : CHARACTER(LEN=default_path_length) :: file_name_tmp
579 : INTEGER :: file_ptr
580 : LOGICAL :: file_exists
581 :
582 500 : CPASSERT(ASSOCIATED(dip))
583 :
584 500 : IF (PRESENT(file_ext)) THEN
585 500 : CPASSERT(file_ext .NE. "")
586 500 : file_name_tmp = expand_file_name_ending(file_name, TRIM(file_ext))
587 : ELSE
588 0 : file_name_tmp = expand_file_name_ending(file_name, "dip")
589 : END IF
590 500 : INQUIRE (FILE=file_name_tmp, EXIST=file_exists)
591 500 : IF (.NOT. file_exists) THEN
592 : CALL open_file(file_name=file_name_tmp, file_status="NEW", &
593 3 : file_action="WRITE", unit_number=file_ptr)
594 3 : WRITE (file_ptr, FMT='(A8,10A20)') "# conf_nr", "dip_x [C Angstrom]", &
595 6 : "dip_y [C Angstrom]", "dip_z [C Angstrom]"
596 : ELSE
597 : CALL open_file(file_name=file_name_tmp, file_status="OLD", &
598 : file_action="WRITE", file_position="APPEND", &
599 497 : unit_number=file_ptr)
600 : END IF
601 500 : WRITE (file_ptr, FMT="(I8,10F20.10)") conf_nr, dip(:)
602 500 : CALL close_file(unit_number=file_ptr)
603 500 : END SUBROUTINE write_dipoles_in_file
604 :
605 : !----------------------------------------------------------------------------
606 : ! read configuration from file
607 : !----------------------------------------------------------------------------
608 :
609 : ! **************************************************************************************************
610 : !> \brief read the trajectory element from a file from sub tree element
611 : !> \param elem actual tree element to be printed out
612 : !> \param tmc_ana TMC analysis environment parameters
613 : !> \param conf_nr Markov chain element number
614 : !> (input the old number and read only if conf nr from file is greater
615 : !> \param stat ...
616 : !> \author Mandes 03.2013
617 : ! **************************************************************************************************
618 2098 : SUBROUTINE read_element_from_file(elem, tmc_ana, conf_nr, stat)
619 : TYPE(tree_type), POINTER :: elem
620 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
621 : INTEGER :: conf_nr, stat
622 :
623 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_element_from_file'
624 :
625 : INTEGER :: conf_nr_old, handle, i_tmp
626 : LOGICAL :: files_conf_missmatch
627 :
628 1049 : stat = TMC_STATUS_OK
629 1049 : conf_nr_old = conf_nr
630 1049 : files_conf_missmatch = .FALSE.
631 :
632 1049 : CPASSERT(ASSOCIATED(elem))
633 1049 : CPASSERT(ASSOCIATED(tmc_ana))
634 1049 : CPASSERT(ASSOCIATED(tmc_ana%atoms))
635 :
636 : ! start the timing
637 1049 : CALL timeset(routineN, handle)
638 :
639 : ! read the coordinates
640 1049 : IF (tmc_ana%id_traj .GT. 0) THEN
641 1049 : i_tmp = conf_nr_old
642 : CALL read_pos_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
643 1049 : conf_nr=i_tmp)
644 1049 : IF (stat .EQ. TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
645 : CALL cp_warn(__LOCATION__, &
646 : 'end of position file reached at line '// &
647 : cp_to_string(REAL(tmc_ana%lc_traj, KIND=dp))//", last element "// &
648 18 : cp_to_string(tmc_ana%last_elem%nr))
649 : ELSE
650 1031 : CPASSERT(i_tmp .GT. conf_nr_old)
651 1031 : conf_nr = i_tmp
652 1031 : elem%nr = i_tmp
653 : END IF
654 : END IF
655 :
656 : ! read the forces
657 : ! TODO if necessary
658 :
659 : ! read the dipoles file
660 1049 : IF (tmc_ana%id_dip .GT. 0 .AND. stat .EQ. TMC_STATUS_OK) THEN
661 0 : i_tmp = conf_nr_old
662 : search_conf_dip: DO
663 : CALL read_dipole_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
664 0 : conf_nr=i_tmp)
665 0 : IF (stat .EQ. TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
666 : CALL cp_warn(__LOCATION__, &
667 : 'end of dipole file reached at line'// &
668 0 : cp_to_string(REAL(tmc_ana%lc_dip, KIND=dp)))
669 0 : EXIT search_conf_dip
670 : END IF
671 : ! check consitence with pos file
672 0 : IF (tmc_ana%id_traj .GT. 0) THEN
673 0 : IF (i_tmp .EQ. conf_nr) THEN
674 : files_conf_missmatch = .FALSE.
675 : EXIT search_conf_dip
676 : ELSE
677 : ! the configuration numbering differ from the position file,
678 : ! but we keep on searching for the correct configuration
679 : files_conf_missmatch = .TRUE.
680 : END IF
681 : ! if no pos file, just take the next conf
682 0 : ELSE IF (i_tmp .GT. conf_nr_old) THEN
683 0 : conf_nr = i_tmp
684 0 : elem%nr = i_tmp
685 0 : EXIT search_conf_dip
686 : END IF
687 : END DO search_conf_dip
688 : END IF
689 :
690 : ! read the cell file
691 1049 : IF (tmc_ana%id_cell .GT. 0 .AND. stat .EQ. TMC_STATUS_OK) THEN
692 : search_conf_cell: DO
693 : CALL read_cell_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
694 1206 : conf_nr=i_tmp)
695 1206 : IF (stat .EQ. TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
696 : CALL cp_warn(__LOCATION__, &
697 : 'end of cell file reached at line at line'// &
698 0 : cp_to_string(REAL(tmc_ana%lc_cell, KIND=dp)))
699 0 : EXIT search_conf_cell
700 : END IF
701 : ! check consitence with pos file
702 1206 : IF (tmc_ana%id_traj .GT. 0) THEN
703 1206 : IF (i_tmp .EQ. conf_nr) THEN
704 : files_conf_missmatch = .FALSE.
705 : EXIT search_conf_cell
706 : ELSE
707 : ! the configuration numbering differ from the position file,
708 : ! but we keep on searching for the correct configuration
709 : files_conf_missmatch = .TRUE.
710 : END IF
711 : ! if no pos file, just take the next conf
712 0 : ELSE IF (i_tmp .GT. conf_nr_old) THEN
713 0 : conf_nr = i_tmp
714 0 : elem%nr = i_tmp
715 0 : EXIT search_conf_cell
716 : END IF
717 : END DO search_conf_cell
718 :
719 : END IF
720 :
721 : ! write the different energies
722 : ! TODO if necessary
723 :
724 1049 : IF (files_conf_missmatch) &
725 : CALL cp_warn(__LOCATION__, &
726 : 'there is a missmatch in the configuration numbering. '// &
727 : "Read number of lines (pos|cell|dip)"// &
728 : cp_to_string(tmc_ana%lc_traj)//"|"// &
729 : cp_to_string(tmc_ana%lc_cell)//"|"// &
730 0 : cp_to_string(tmc_ana%lc_dip))
731 :
732 : ! end the timing
733 1049 : CALL timestop(handle)
734 1049 : END SUBROUTINE read_element_from_file
735 :
736 : ! **************************************************************************************************
737 : !> \brief search for the next configurational position in file
738 : !> \param elem actual tree element to be read
739 : !> \param tmc_ana ...
740 : !> \param stat ...
741 : !> \param conf_nr Markov chain element number
742 : !> (input the old number and read only if conf nr from file is greater
743 : !> \param header_info ...
744 : !> \author Mandes 03.2013
745 : ! **************************************************************************************************
746 2098 : SUBROUTINE read_pos_from_file(elem, tmc_ana, stat, conf_nr, header_info)
747 : TYPE(tree_type), POINTER :: elem
748 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
749 : INTEGER :: stat, conf_nr
750 : CHARACTER(LEN=*), OPTIONAL :: header_info
751 :
752 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_pos_from_file'
753 :
754 : CHARACTER(LEN=default_string_length) :: c_tmp
755 : INTEGER :: handle, i, i_tmp, status
756 :
757 1049 : stat = TMC_STATUS_FAILED
758 :
759 0 : CPASSERT(ASSOCIATED(elem))
760 1049 : CPASSERT(ASSOCIATED(elem%pos))
761 1049 : CPASSERT(ASSOCIATED(tmc_ana))
762 1049 : CPASSERT(tmc_ana%id_traj .GT. 0)
763 :
764 : ! start the timing
765 1049 : CALL timeset(routineN, handle)
766 :
767 : search_next_conf: DO
768 6105 : c_tmp(:) = " "
769 6105 : tmc_ana%lc_traj = tmc_ana%lc_traj + 1
770 6105 : READ (tmc_ana%id_traj, '(A)', IOSTAT=status) c_tmp(:)
771 6105 : IF (status .GT. 0) &
772 : CALL cp_abort(__LOCATION__, &
773 : "configuration header read error at line: "// &
774 0 : cp_to_string(tmc_ana%lc_traj)//": "//c_tmp)
775 6105 : IF (status .LT. 0) THEN ! end of file reached
776 18 : stat = TMC_STATUS_WAIT_FOR_NEW_TASK
777 18 : EXIT search_next_conf
778 : END IF
779 6087 : IF (INDEX(c_tmp, "=") .GT. 0) THEN
780 1206 : READ (c_tmp(INDEX(c_tmp, "=") + 1:), *, IOSTAT=status) i_tmp ! read the configuration number
781 1206 : IF (status .NE. 0) &
782 : CALL cp_abort(__LOCATION__, &
783 : "configuration header read error (for conf nr) at line: "// &
784 0 : cp_to_string(tmc_ana%lc_traj))
785 1206 : IF (i_tmp .GT. conf_nr) THEN
786 : ! TODO we could also read the energy ...
787 1031 : conf_nr = i_tmp
788 1031 : IF (PRESENT(header_info)) header_info = c_tmp
789 1031 : stat = TMC_STATUS_OK
790 1031 : EXIT search_next_conf
791 : END IF
792 : END IF
793 : END DO search_next_conf
794 :
795 1049 : IF (stat .EQ. TMC_STATUS_OK) THEN
796 22682 : pos_loop: DO i = 1, SIZE(elem%pos), tmc_ana%dim_per_elem
797 21651 : tmc_ana%lc_traj = tmc_ana%lc_traj + 1
798 : READ (tmc_ana%id_traj, FMT="(A4,1X,1000F20.10)", IOSTAT=status) &
799 21651 : c_tmp, elem%pos(i:i + tmc_ana%dim_per_elem - 1)
800 22682 : IF (status .NE. 0) THEN
801 : CALL cp_abort(__LOCATION__, &
802 : "configuration pos read error at line: "// &
803 0 : cp_to_string(tmc_ana%lc_traj))
804 : END IF
805 : END DO pos_loop
806 65984 : elem%pos(:) = elem%pos(:)/au2a
807 : END IF
808 :
809 : ! end the timing
810 1049 : CALL timestop(handle)
811 1049 : END SUBROUTINE read_pos_from_file
812 :
813 : ! **************************************************************************************************
814 : !> \brief search for the dipole entry
815 : !> \param elem actual tree element to be read
816 : !> \param tmc_ana ...
817 : !> \param stat ...
818 : !> \param conf_nr Markov chain element number
819 : !> (input the old number and read only if conf nr from file is greater
820 : !> \author Mandes 03.2013
821 : ! **************************************************************************************************
822 0 : SUBROUTINE read_dipole_from_file(elem, tmc_ana, stat, conf_nr)
823 : TYPE(tree_type), POINTER :: elem
824 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
825 : INTEGER :: stat, conf_nr
826 :
827 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_dipole_from_file'
828 :
829 : CHARACTER(LEN=250) :: c_tmp
830 : INTEGER :: handle, status
831 :
832 0 : stat = TMC_STATUS_FAILED
833 :
834 0 : CPASSERT(ASSOCIATED(elem))
835 0 : CPASSERT(ASSOCIATED(elem%dipole))
836 0 : CPASSERT(ASSOCIATED(tmc_ana))
837 0 : CPASSERT(tmc_ana%id_dip .GT. 0)
838 :
839 : ! start the timing
840 0 : CALL timeset(routineN, handle)
841 0 : tmc_ana%lc_dip = tmc_ana%lc_dip + 1
842 0 : READ (tmc_ana%id_dip, FMT="(A)", IOSTAT=status) c_tmp
843 0 : IF (status .EQ. 0) THEN
844 : ! skip the initial line (header)
845 0 : IF (INDEX(c_tmp, "#") .GT. 0) THEN
846 0 : tmc_ana%lc_dip = tmc_ana%lc_dip + 1
847 0 : READ (tmc_ana%id_dip, FMT="(A)", IOSTAT=status) c_tmp
848 : END IF
849 : END IF
850 0 : IF (status .EQ. 0) THEN
851 : READ (c_tmp, FMT="(I8,10F20.10)", IOSTAT=status) &
852 0 : conf_nr, elem%dipole(:)
853 : END IF
854 0 : IF (status .EQ. 0) THEN ! success
855 0 : stat = TMC_STATUS_OK
856 0 : ELSE IF (status .LT. 0) THEN ! end of file reached
857 0 : stat = TMC_STATUS_WAIT_FOR_NEW_TASK
858 : ELSE
859 : IF (status .NE. 0) &
860 0 : CPWARN("configuration dipole read error at line: "//cp_to_string(tmc_ana%lc_dip))
861 0 : stat = TMC_STATUS_FAILED
862 : END IF
863 :
864 : ! end the timing
865 0 : CALL timestop(handle)
866 0 : END SUBROUTINE read_dipole_from_file
867 :
868 : ! **************************************************************************************************
869 : !> \brief search for the cell entry
870 : !> \param elem actual tree element to be read
871 : !> \param tmc_ana ...
872 : !> \param stat ...
873 : !> \param conf_nr Markov chain element number
874 : !> (input the old number and read only if conf nr from file is greater
875 : !> \author Mandes 03.2013
876 : ! **************************************************************************************************
877 2412 : SUBROUTINE read_cell_from_file(elem, tmc_ana, stat, conf_nr)
878 : TYPE(tree_type), POINTER :: elem
879 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
880 : INTEGER :: stat, conf_nr
881 :
882 : CHARACTER(LEN=*), PARAMETER :: routineN = 'read_cell_from_file'
883 :
884 : CHARACTER(LEN=250) :: c_tmp
885 : INTEGER :: handle, status
886 : REAL(KIND=dp) :: r_tmp
887 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
888 :
889 1206 : stat = TMC_STATUS_FAILED
890 :
891 1206 : CPASSERT(ASSOCIATED(elem))
892 1206 : CPASSERT(ASSOCIATED(tmc_ana))
893 1206 : CPASSERT(ASSOCIATED(tmc_ana%cell))
894 1206 : CPASSERT(tmc_ana%id_cell .GT. 0)
895 :
896 : ! start the timing
897 1206 : CALL timeset(routineN, handle)
898 :
899 1206 : tmc_ana%lc_cell = tmc_ana%lc_cell + 1
900 1206 : READ (tmc_ana%id_cell, FMT="(A)", IOSTAT=status) c_tmp
901 1206 : IF (status .EQ. 0) THEN
902 : ! skip the initial line (header)
903 1206 : IF (INDEX(c_tmp, "#") .GT. 0) THEN
904 18 : tmc_ana%lc_cell = tmc_ana%lc_cell + 1
905 18 : READ (tmc_ana%id_cell, FMT="(A)", IOSTAT=status) c_tmp
906 : END IF
907 : END IF
908 1206 : IF (status .EQ. 0) THEN
909 1206 : READ (c_tmp, FMT="(I8,9(1X,F19.10),1X,F24.10)", IOSTAT=status) conf_nr, &
910 2412 : hmat(:, :), r_tmp
911 : END IF
912 1206 : IF (status .LT. 0) THEN ! end of file reached
913 0 : stat = TMC_STATUS_WAIT_FOR_NEW_TASK
914 1206 : ELSE IF (status .GT. 0) THEN
915 : IF (status .NE. 0) &
916 0 : CPABORT("configuration cell read error at line: "//cp_to_string(tmc_ana%lc_cell))
917 0 : stat = TMC_STATUS_FAILED
918 : ELSE
919 1206 : IF (elem%nr .LT. 0) elem%nr = conf_nr
920 15678 : hmat(:, :) = hmat(:, :)/au2a
921 : ! get the box scaling
922 : CALL get_cell_scaling(cell=tmc_ana%cell, scaled_hmat=hmat, &
923 1206 : box_scale=elem%box_scale)
924 1206 : stat = TMC_STATUS_OK
925 : END IF
926 : ! end the timing
927 1206 : CALL timestop(handle)
928 1206 : END SUBROUTINE read_cell_from_file
929 :
930 : !----------------------------------------------------------------------------
931 : ! get the configurations from file and calc
932 : !----------------------------------------------------------------------------
933 :
934 : ! **************************************************************************************************
935 : !> \brief opens the files for reading configurations data to analyze
936 : !> \param tmc_ana ...
937 : !> \param stat ...
938 : !> \param dir_ind ...
939 : !> \param
940 : !> \author Mandes 02.2013
941 : ! **************************************************************************************************
942 36 : SUBROUTINE analyse_files_open(tmc_ana, stat, dir_ind)
943 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
944 : INTEGER :: stat
945 : INTEGER, OPTIONAL :: dir_ind
946 :
947 : CHARACTER(LEN=*), PARAMETER :: routineN = 'analyse_files_open'
948 :
949 : CHARACTER(LEN=default_path_length) :: dir_name, file_name_act, file_name_temp
950 : INTEGER :: handle
951 : LOGICAL :: file_exists
952 :
953 18 : CPASSERT(ASSOCIATED(tmc_ana))
954 :
955 18 : stat = TMC_STATUS_WAIT_FOR_NEW_TASK
956 :
957 : ! start the timing
958 18 : CALL timeset(routineN, handle)
959 :
960 18 : IF (PRESENT(dir_ind)) THEN
961 18 : CPASSERT(ASSOCIATED(tmc_ana%dirs))
962 18 : CPASSERT(dir_ind .GT. 0)
963 18 : CPASSERT(dir_ind .LE. SIZE(tmc_ana%dirs))
964 :
965 18 : IF (INDEX(tmc_ana%dirs(dir_ind), "/", BACK=.TRUE.) .EQ. &
966 : LEN_TRIM(tmc_ana%dirs(dir_ind))) THEN
967 18 : dir_name = TRIM(tmc_ana%dirs(dir_ind))
968 : ELSE
969 0 : dir_name = TRIM(tmc_ana%dirs(dir_ind))//"/"
970 : END IF
971 : ELSE
972 0 : dir_name = "./"
973 : END IF
974 :
975 : ! open the files
976 : file_name_temp = expand_file_name_temp( &
977 : file_name=tmc_default_trajectory_file_name, &
978 18 : rvalue=tmc_ana%temperature)
979 : ! position file
980 18 : IF (tmc_ana%costum_pos_file_name .NE. "") THEN
981 0 : file_name_act = TRIM(dir_name)//tmc_ana%costum_pos_file_name
982 : ELSE
983 : file_name_act = TRIM(dir_name)// &
984 18 : expand_file_name_ending(file_name_temp, "xyz")
985 : END IF
986 18 : INQUIRE (FILE=file_name_act, EXIST=file_exists)
987 18 : IF (file_exists) THEN
988 : CALL open_file(file_name=file_name_act, file_status="OLD", &
989 18 : file_action="READ", unit_number=tmc_ana%id_traj)
990 18 : WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
991 36 : "read xyz file", TRIM(file_name_act)
992 : END IF
993 :
994 : ! cell file
995 18 : IF (tmc_ana%costum_cell_file_name .NE. "") THEN
996 0 : file_name_act = TRIM(dir_name)//tmc_ana%costum_cell_file_name
997 : ELSE
998 : file_name_act = TRIM(dir_name)// &
999 18 : expand_file_name_ending(file_name_temp, "cell")
1000 : END IF
1001 18 : INQUIRE (FILE=file_name_act, EXIST=file_exists)
1002 18 : IF (file_exists) THEN
1003 : CALL open_file(file_name=file_name_act, file_status="OLD", &
1004 18 : file_action="READ", unit_number=tmc_ana%id_cell)
1005 18 : WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
1006 36 : "read cell file", TRIM(file_name_act)
1007 : END IF
1008 :
1009 : ! dipole file
1010 18 : IF (tmc_ana%costum_dip_file_name .NE. "") THEN
1011 18 : file_name_act = TRIM(dir_name)//tmc_ana%costum_dip_file_name
1012 : ELSE
1013 : file_name_act = TRIM(dir_name)// &
1014 0 : expand_file_name_ending(file_name_temp, "dip")
1015 : END IF
1016 18 : INQUIRE (FILE=file_name_act, EXIST=file_exists)
1017 18 : IF (file_exists) THEN
1018 : CALL open_file(file_name=file_name_act, file_status="OLD", &
1019 0 : file_action="READ", unit_number=tmc_ana%id_dip)
1020 0 : WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
1021 0 : "read dip file", TRIM(file_name_act)
1022 : END IF
1023 :
1024 18 : IF (tmc_ana%id_traj .GT. 0 .OR. tmc_ana%id_cell .GT. 0 .OR. &
1025 : tmc_ana%id_dip .GT. 0) THEN
1026 18 : stat = TMC_STATUS_OK
1027 : ELSE
1028 : CALL cp_warn(__LOCATION__, &
1029 : "There is no file to open for temperature "//cp_to_string(tmc_ana%temperature)// &
1030 0 : "K in directory "//TRIM(dir_name))
1031 : END IF
1032 : ! end the timing
1033 18 : CALL timestop(handle)
1034 18 : END SUBROUTINE analyse_files_open
1035 :
1036 : ! **************************************************************************************************
1037 : !> \brief close the files for reading configurations data to analyze
1038 : !> \param tmc_ana ...
1039 : !> \param
1040 : !> \author Mandes 02.2013
1041 : ! **************************************************************************************************
1042 36 : SUBROUTINE analyse_files_close(tmc_ana)
1043 : TYPE(tmc_analysis_env), POINTER :: tmc_ana
1044 :
1045 : CHARACTER(LEN=*), PARAMETER :: routineN = 'analyse_files_close'
1046 :
1047 : INTEGER :: handle
1048 :
1049 18 : CPASSERT(ASSOCIATED(tmc_ana))
1050 :
1051 : ! start the timing
1052 18 : CALL timeset(routineN, handle)
1053 :
1054 : ! position file
1055 18 : IF (tmc_ana%id_traj .GT. 0) CALL close_file(unit_number=tmc_ana%id_traj)
1056 :
1057 : ! cell file
1058 18 : IF (tmc_ana%id_cell .GT. 0) CALL close_file(unit_number=tmc_ana%id_cell)
1059 :
1060 : ! dipole file
1061 18 : IF (tmc_ana%id_dip .GT. 0) CALL close_file(unit_number=tmc_ana%id_dip)
1062 :
1063 : ! end the timing
1064 18 : CALL timestop(handle)
1065 18 : END SUBROUTINE analyse_files_close
1066 :
1067 : END MODULE tmc_file_io
|
