LCOV - code coverage report
Current view: top level - src/tmc - tmc_file_io.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:2c0d679) Lines: 80.2 % 419 336
Test Date: 2026-09-25 00:58:37 Functions: 93.3 % 15 14

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2026 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         2608 :    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 /= "")
      81              : 
      82         2608 :       ind = INDEX(file_name, ".", BACK=.TRUE.)
      83         2608 :       IF (.NOT. ind == 0) THEN
      84         2608 :          WRITE (result_file_name, *) file_name(1:ind - 1), ".", &
      85         5216 :             TRIM(ADJUSTL(extra))
      86              :       ELSE
      87            0 :          WRITE (result_file_name, *) TRIM(file_name), ".", extra
      88              :       END IF
      89         2608 :       result_file_name = TRIM(ADJUSTL(result_file_name))
      90         2608 :       CPASSERT(result_file_name /= "")
      91         2608 :    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 /= "")
     108              : 
     109         1427 :       ind = INDEX(file_name, ".", BACK=.TRUE.)
     110         1427 :       IF (.NOT. ind == 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 /= "")
     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         2864 :    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         2864 :       CPASSERT(file_name /= "")
     137              : 
     138         2864 :       rval_to_string = ""
     139              : 
     140         2864 :       WRITE (rval_to_string, "(F16.2)") rvalue
     141         2864 :       ind = INDEX(file_name, ".", BACK=.TRUE.)
     142         2864 :       IF (.NOT. ind == 0) THEN
     143         2864 :          WRITE (result_file_name, *) file_name(1:ind - 1), "_T", &
     144         5728 :             TRIM(ADJUSTL(rval_to_string)), file_name(ind:LEN_TRIM(file_name))
     145              :       ELSE
     146            0 :          IF (LEN(file_name) == 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         2864 :       result_file_name = TRIM(ADJUSTL(result_file_name))
     154         2864 :       CPASSERT(result_file_name /= "")
     155         2864 :    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 /= "")
     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 == 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) == 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 /= "")
     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           12 :       WRITE (f_unit) tmc_env%params%Temp(:), &
     230            3 :          tmc_env%m_env%gt_act%nr, &
     231           84 :          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           15 :          tmc_env%m_env%result_count(0:), &
     237           33 :          tmc_env%params%move_types%mv_weight, &
     238          111 :          tmc_env%params%move_types%acc_count, &
     239          111 :          tmc_env%params%move_types%mv_count, &
     240          102 :          tmc_env%params%move_types%subbox_acc_count, &
     241          102 :          tmc_env%params%move_types%subbox_count, &
     242           39 :          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              :       REAL(KIND=dp), PARAMETER                           :: eps_temp = 5.0E-03_dp
     283              : 
     284              :       INTEGER                                            :: file_ptr, i, temp_size
     285              :       LOGICAL                                            :: flag
     286            2 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: tmp_temp
     287              :       REAL(KIND=dp), DIMENSION(nr_mv_types)              :: mv_weight_tmp
     288              : 
     289            2 :       CPASSERT(ASSOCIATED(tmc_env))
     290            2 :       CPASSERT(ASSOCIATED(tmc_env%m_env))
     291            2 :       CPASSERT(ASSOCIATED(tmc_env%params))
     292            2 :       CPASSERT(ASSOCIATED(tmc_env%m_env%gt_act))
     293              : 
     294            2 :       IF (file_name == tmc_default_restart_in_file_name) THEN
     295            2 :          INQUIRE (FILE=tmc_default_restart_in_file_name, EXIST=flag)
     296            2 :          CPASSERT(flag)
     297              :          CALL open_file(file_name=tmc_default_restart_in_file_name, file_status="OLD", &
     298            2 :                         file_action="READ", unit_number=file_ptr)
     299            2 :          READ (file_ptr, *) file_name
     300            2 :          CALL close_file(unit_number=file_ptr)
     301              :       END IF
     302              : 
     303              :       CALL open_file(file_name=file_name, file_status="OLD", file_form="UNFORMATTED", &
     304            2 :                      file_action="READ", unit_number=file_ptr)
     305            2 :       READ (file_ptr) temp_size
     306            2 :       IF (temp_size /= SIZE(tmc_env%params%Temp)) THEN
     307              :          CALL cp_abort(__LOCATION__, &
     308              :                        "the actual specified temperatures does not "// &
     309            0 :                        "fit in amount with the one from restart file ")
     310              :       END IF
     311            6 :       ALLOCATE (tmp_temp(temp_size))
     312            2 :       READ (file_ptr) tmp_temp(:), &
     313            2 :          tmc_env%m_env%gt_act%nr, &
     314           56 :          tmc_env%m_env%gt_act%rng_seed, &
     315            2 :          tmc_env%m_env%gt_act%rnd_nr, &
     316            2 :          tmc_env%m_env%gt_act%prob_acc, &
     317            2 :          tmc_env%m_env%gt_act%mv_conf, & !
     318            2 :          tmc_env%m_env%gt_act%mv_next_conf, & !
     319           10 :          tmc_env%m_env%result_count(0:), &
     320            2 :          mv_weight_tmp, & !
     321           74 :          tmc_env%params%move_types%acc_count, &
     322           74 :          tmc_env%params%move_types%mv_count, &
     323           68 :          tmc_env%params%move_types%subbox_acc_count, &
     324           68 :          tmc_env%params%move_types%subbox_count, & !
     325           26 :          tmc_env%params%cell%hmat, &
     326            2 :          job_counts, &
     327            4 :          timings
     328              : 
     329            8 :       IF (ANY(ABS(tmc_env%params%Temp(:) - tmp_temp(:)) >= eps_temp)) THEN
     330              :          CALL cp_abort(__LOCATION__, "the temperatures differ from the previous calculation. "// &
     331            0 :                        "There were the following temperatures used:")
     332              :       END IF
     333           22 :       IF (ANY(mv_weight_tmp(:) /= tmc_env%params%move_types%mv_weight(:))) THEN
     334            0 :          CPWARN("The amount of mv types differs between the original and the restart run.")
     335              :       END IF
     336              : 
     337            8 :       DO i = 1, SIZE(tmc_env%params%Temp)
     338            6 :          tmc_env%m_env%gt_act%conf(i)%elem => tmc_env%m_env%result_list(i)%elem
     339            6 :          READ (file_ptr) tmc_env%m_env%result_list(i)%elem%nr, &
     340          168 :             tmc_env%m_env%result_list(i)%elem%rng_seed, &
     341          384 :             tmc_env%m_env%result_list(i)%elem%pos, &
     342          384 :             tmc_env%m_env%result_list(i)%elem%vel, &
     343           24 :             tmc_env%m_env%result_list(i)%elem%box_scale, &
     344            6 :             tmc_env%m_env%result_list(i)%elem%potential, &
     345            6 :             tmc_env%m_env%result_list(i)%elem%e_pot_approx, &
     346            6 :             tmc_env%m_env%result_list(i)%elem%ekin, &
     347            6 :             tmc_env%m_env%result_list(i)%elem%ekin_before_md, &
     348           14 :             tmc_env%m_env%result_list(i)%elem%temp_created
     349              :       END DO
     350            2 :       CALL close_file(unit_number=file_ptr)
     351            2 :    END SUBROUTINE read_restart_file
     352              : 
     353              :    !----------------------------------------------------------------------------
     354              :    ! printing configuration in file
     355              :    !----------------------------------------------------------------------------
     356              : 
     357              : ! **************************************************************************************************
     358              : !> \brief select the correct configuration to print out the
     359              : !>        (coordinates, forces, cell ...)
     360              : !> \param result_list list of configurations for each temperature
     361              : !> \param result_count list with number of Markov Chain number
     362              : !>          for each teperature (index 0 for global tree)
     363              : !> \param conf_updated index of the updated (modified element)
     364              : !> \param accepted acceptance flag
     365              : !> \param tmc_params TMC environment parameters
     366              : !> \author Mandes 02.2013
     367              : ! **************************************************************************************************
     368         9086 :    SUBROUTINE write_result_list_element(result_list, result_count, conf_updated, &
     369              :                                         accepted, tmc_params)
     370              :       TYPE(elem_array_type), DIMENSION(:), POINTER       :: result_list
     371              :       INTEGER, DIMENSION(:), POINTER                     :: result_count
     372              :       INTEGER                                            :: conf_updated
     373              :       LOGICAL, INTENT(IN)                                :: accepted
     374              :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     375              : 
     376              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_result_list_element'
     377              : 
     378              :       CHARACTER(LEN=default_path_length)                 :: file_name
     379              :       INTEGER                                            :: handle, i
     380              : 
     381         4543 :       file_name = ""
     382              : 
     383         4543 :       CPASSERT(ASSOCIATED(result_list))
     384         4543 :       CPASSERT(ASSOCIATED(result_count))
     385         4543 :       CPASSERT(ASSOCIATED(tmc_params))
     386         4543 :       CPASSERT(ASSOCIATED(tmc_params%Temp))
     387         4543 :       CPASSERT(conf_updated >= 0)
     388         4543 :       CPASSERT(conf_updated <= SIZE(tmc_params%Temp))
     389              : 
     390              :       ! start the timing
     391         4543 :       CALL timeset(routineN, handle)
     392              : 
     393         4543 :       IF (conf_updated == 0) THEN
     394              :          ! for debugging print every configuration of every temperature
     395            0 :          DO i = 1, SIZE(tmc_params%Temp)
     396            0 :             WRITE (file_name, *) "every_step_", TRIM(tmc_default_trajectory_file_name)
     397              :             CALL write_element_in_file(elem=result_list(i)%elem, &
     398              :                                        tmc_params=tmc_params, conf_nr=result_count(0), &
     399            0 :                                        file_name=expand_file_name_temp(file_name=file_name, rvalue=tmc_params%Temp(i)))
     400              :          END DO
     401              :       ELSE
     402         4543 :          IF ((.NOT. tmc_params%print_only_diff_conf) .OR. &
     403              :              (tmc_params%print_only_diff_conf .AND. accepted)) THEN
     404              :             CALL write_element_in_file(elem=result_list(conf_updated)%elem, &
     405              :                                        tmc_params=tmc_params, conf_nr=result_count(conf_updated), &
     406              :                                        file_name=expand_file_name_temp(file_name=TRIM(tmc_default_trajectory_file_name), &
     407         1018 :                                                                        rvalue=tmc_params%Temp(conf_updated)))
     408              :          END IF
     409              :       END IF
     410              :       ! end the timing
     411         4543 :       CALL timestop(handle)
     412         4543 :    END SUBROUTINE write_result_list_element
     413              : 
     414              : ! **************************************************************************************************
     415              : !> \brief writes the trajectory element in a file from sub tree element
     416              : !> \param elem actual tree element to be printed out
     417              : !> \param tmc_params TMC environment parameters
     418              : !> \param temp_index ...
     419              : !> \param file_name file name will be extended by type of file (pos, cell,...)
     420              : !> \param conf_nr Markov chain element number
     421              : !> \param conf_info whole header line
     422              : !> \author Mandes 11.2012
     423              : ! **************************************************************************************************
     424         1018 :    SUBROUTINE write_element_in_file(elem, tmc_params, temp_index, file_name, conf_nr, &
     425              :                                     conf_info)
     426              :       TYPE(tree_type), POINTER                           :: elem
     427              :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     428              :       INTEGER, OPTIONAL                                  :: temp_index
     429              :       CHARACTER(LEN=*), OPTIONAL                         :: file_name
     430              :       INTEGER, OPTIONAL                                  :: conf_nr
     431              :       CHARACTER(LEN=*), OPTIONAL                         :: conf_info
     432              : 
     433              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'write_element_in_file'
     434              : 
     435              :       CHARACTER(LEN=default_path_length)                 :: file_name_act, tmp_name
     436              :       CHARACTER(LEN=default_string_length)               :: header
     437              :       INTEGER                                            :: file_ptr, handle, i, nr_atoms
     438              :       LOGICAL                                            :: file_exists, print_it
     439              :       REAL(KIND=dp)                                      :: vol
     440              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat_scaled
     441              : 
     442         1018 :       file_name_act = ""
     443         1018 :       tmp_name = ""
     444         1018 :       header = ""
     445         1018 :       print_it = .TRUE.
     446              : 
     447            0 :       CPASSERT(ASSOCIATED(elem))
     448         1018 :       CPASSERT(ASSOCIATED(tmc_params))
     449         1018 :       CPASSERT(ASSOCIATED(tmc_params%atoms))
     450         1018 :       CPASSERT(PRESENT(conf_nr) .OR. PRESENT(conf_info))
     451              : 
     452         1018 :       IF (print_it) THEN
     453              :          ! start the timing
     454         1018 :          CALL timeset(routineN, handle)
     455              : 
     456              :          ! set default file name
     457         1018 :          IF (PRESENT(file_name)) THEN
     458         1018 :             CPASSERT(file_name /= "")
     459         1018 :             file_name_act = file_name
     460              :          ELSE
     461            0 :             CPASSERT(ASSOCIATED(tmc_params%Temp))
     462            0 :             CPASSERT(PRESENT(temp_index))
     463              :             file_name_act = expand_file_name_temp(file_name=tmc_default_trajectory_file_name, &
     464            0 :                                                   rvalue=tmc_params%Temp(temp_index))
     465              :          END IF
     466              : 
     467         1018 :          nr_atoms = SIZE(elem%pos)/tmc_params%dim_per_elem
     468              : 
     469              :          ! set header (for coordinate or force file)
     470         1018 :          IF (tmc_params%print_trajectory .OR. tmc_params%print_forces) THEN
     471         1018 :             IF (PRESENT(conf_info)) THEN
     472            0 :                WRITE (header, *) TRIM(ADJUSTL(conf_info))
     473              :             ELSE
     474              :                !WRITE(header,FMT="(A,I8,A,F20.10)") " i = ", conf_nr,", E = ", elem%potential
     475         1018 :                WRITE (header, FMT="(A,I8,A,F20.10,F20.10,A,I8,I8)") "i =", conf_nr, " ,E =", &
     476         2036 :                   elem%potential, elem%ekin, " st elem", elem%sub_tree_nr, elem%nr
     477              :             END IF
     478              :          END IF
     479              : 
     480              :          ! write the coordinates
     481         1018 :          IF (tmc_params%print_trajectory) THEN
     482         1018 :             tmp_name = expand_file_name_ending(file_name_act, "xyz")
     483              :             CALL open_file(file_name=tmp_name, file_status="UNKNOWN", &
     484              :                            file_action="WRITE", file_position="APPEND", &
     485         1018 :                            unit_number=file_ptr)
     486         1018 :             WRITE (file_ptr, FMT="(I8)") nr_atoms
     487         1018 :             WRITE (file_ptr, *) TRIM(header)
     488        44643 :             DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
     489              :                WRITE (file_ptr, FMT="(A4,1X,1000F20.10)") &
     490        43625 :                   TRIM(tmc_params%atoms((i - 1)/tmc_params%dim_per_elem + 1)%name), &
     491       219143 :                   elem%pos(i:i + tmc_params%dim_per_elem - 1)*au2a
     492              :             END DO
     493         1018 :             CALL close_file(unit_number=file_ptr)
     494              :          END IF
     495              : 
     496              :          ! write the forces
     497         1018 :          IF (tmc_params%print_forces) THEN
     498          331 :             tmp_name = expand_file_name_ending(file_name_act, "frc")
     499              :             CALL open_file(file_name=tmp_name, file_status="UNKNOWN", &
     500              :                            file_action="WRITE", file_position="APPEND", &
     501          331 :                            unit_number=file_ptr)
     502          331 :             WRITE (file_ptr, FMT="(I8)") nr_atoms
     503          331 :             WRITE (file_ptr, *) TRIM(header)
     504         7282 :             DO i = 1, SIZE(elem%pos), tmc_params%dim_per_elem
     505              :                WRITE (file_ptr, FMT="(A4,1X,1000F20.10)") &
     506         6951 :                   TRIM(tmc_params%atoms((i - 1)/tmc_params%dim_per_elem + 1)%name), &
     507        35086 :                   elem%frc(i:i + tmc_params%dim_per_elem - 1)
     508              :             END DO
     509          331 :             CALL close_file(unit_number=file_ptr)
     510              :          END IF
     511              : 
     512              :          ! write the cell dipoles
     513         1018 :          IF (tmc_params%print_dipole) THEN
     514              :             CALL write_dipoles_in_file(file_name=file_name_act, &
     515            0 :                                        conf_nr=conf_nr, dip=elem%dipole)
     516              :          END IF
     517              : 
     518              :          ! write the cell file
     519         1018 :          IF (tmc_params%print_cell) THEN
     520          392 :             tmp_name = expand_file_name_ending(file_name_act, "cell")
     521              :             ! header
     522          392 :             INQUIRE (FILE=tmp_name, EXIST=file_exists) ! file_exists will be TRUE if the file exist
     523          392 :             IF (.NOT. file_exists) THEN
     524              :                CALL open_file(file_name=tmp_name, file_status="NEW", &
     525            6 :                               file_action="WRITE", unit_number=file_ptr)
     526              :                WRITE (file_ptr, FMT='(A,9(7X,A2," [Angstrom]"),6X,A)') &
     527            6 :                   "# MC step ", "Ax", "Ay", "Az", "Bx", "By", "Bz", "Cx", "Cy", "Cz", &
     528           12 :                   "Volume [Angstrom^3]"
     529              :             ELSE
     530              :                CALL open_file(file_name=tmp_name, file_status="OLD", &
     531              :                               file_action="WRITE", file_position="APPEND", &
     532          386 :                               unit_number=file_ptr)
     533              :             END IF
     534              :             CALL get_scaled_cell(cell=tmc_params%cell, &
     535              :                                  box_scale=elem%box_scale, scaled_hmat=hmat_scaled, &
     536          392 :                                  vol=vol)
     537          392 :             WRITE (file_ptr, FMT="(I8,9(1X,F19.10),1X,F24.10)") conf_nr, &
     538         5488 :                hmat_scaled(:, :)*au2a, vol*au2a**3
     539              :             !TODO better cell output e.g. using cell_types routine
     540          392 :             CALL close_file(unit_number=file_ptr)
     541              :          END IF
     542              : 
     543              :          ! write the different energies
     544         1018 :          IF (tmc_params%print_energies) THEN
     545          331 :             tmp_name = expand_file_name_ending(file_name_act, "ener")
     546              :             ! header
     547          331 :             INQUIRE (FILE=tmp_name, EXIST=file_exists) ! file_exists will be TRUE if the file exist
     548          331 :             IF (.NOT. file_exists) THEN
     549              :                CALL open_file(file_name=tmp_name, file_status="NEW", &
     550            3 :                               file_action="WRITE", unit_number=file_ptr)
     551              :                WRITE (file_ptr, FMT='(A,4A20)') &
     552            3 :                   "# MC step ", " exact ", " approx ", " last SCF ", " kinetic "
     553              :             ELSE
     554              :                CALL open_file(file_name=tmp_name, file_status="OLD", &
     555              :                               file_action="WRITE", file_position="APPEND", &
     556          328 :                               unit_number=file_ptr)
     557              :             END IF
     558          331 :             WRITE (file_ptr, FMT="(I8,14F20.10)") conf_nr, elem%potential, elem%e_pot_approx, &
     559          662 :                elem%scf_energies(MOD(elem%scf_energies_count, 4) + 1), elem%ekin
     560          331 :             CALL close_file(unit_number=file_ptr)
     561              :          END IF
     562              : 
     563              :          ! end the timing
     564         1018 :          CALL timestop(handle)
     565              :       END IF
     566         1018 :    END SUBROUTINE write_element_in_file
     567              : 
     568              : ! **************************************************************************************************
     569              : !> \brief writes the cell dipoles in dipole trajectory file
     570              : !> \param file_name ...
     571              : !> \param conf_nr ...
     572              : !> \param dip ...
     573              : !> \param file_ext ...
     574              : !> \param
     575              : !> \author Mandes 11.2012
     576              : ! **************************************************************************************************
     577          500 :    SUBROUTINE write_dipoles_in_file(file_name, conf_nr, dip, file_ext)
     578              :       CHARACTER(LEN=default_path_length)                 :: file_name
     579              :       INTEGER                                            :: conf_nr
     580              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: dip
     581              :       CHARACTER(LEN=*), INTENT(in), OPTIONAL             :: file_ext
     582              : 
     583              :       CHARACTER(LEN=default_path_length)                 :: file_name_tmp
     584              :       INTEGER                                            :: file_ptr
     585              :       LOGICAL                                            :: file_exists
     586              : 
     587          500 :       CPASSERT(ASSOCIATED(dip))
     588              : 
     589          500 :       IF (PRESENT(file_ext)) THEN
     590          500 :          CPASSERT(file_ext /= "")
     591          500 :          file_name_tmp = expand_file_name_ending(file_name, TRIM(file_ext))
     592              :       ELSE
     593            0 :          file_name_tmp = expand_file_name_ending(file_name, "dip")
     594              :       END IF
     595          500 :       INQUIRE (FILE=file_name_tmp, EXIST=file_exists)
     596          500 :       IF (.NOT. file_exists) THEN
     597              :          CALL open_file(file_name=file_name_tmp, file_status="NEW", &
     598            3 :                         file_action="WRITE", unit_number=file_ptr)
     599            3 :          WRITE (file_ptr, FMT='(A8,10A20)') "# conf_nr", "dip_x [C Angstrom]", &
     600            6 :             "dip_y [C Angstrom]", "dip_z [C Angstrom]"
     601              :       ELSE
     602              :          CALL open_file(file_name=file_name_tmp, file_status="OLD", &
     603              :                         file_action="WRITE", file_position="APPEND", &
     604          497 :                         unit_number=file_ptr)
     605              :       END IF
     606         2000 :       WRITE (file_ptr, FMT="(I8,10F20.10)") conf_nr, dip(:)
     607          500 :       CALL close_file(unit_number=file_ptr)
     608          500 :    END SUBROUTINE write_dipoles_in_file
     609              : 
     610              :    !----------------------------------------------------------------------------
     611              :    ! read configuration from file
     612              :    !----------------------------------------------------------------------------
     613              : 
     614              : ! **************************************************************************************************
     615              : !> \brief read the trajectory element from a file from sub tree element
     616              : !> \param elem actual tree element to be printed out
     617              : !> \param tmc_ana TMC analysis environment parameters
     618              : !> \param conf_nr Markov chain element number
     619              : !>        (input the old number and read only if conf nr from file is greater
     620              : !> \param stat ...
     621              : !> \author Mandes 03.2013
     622              : ! **************************************************************************************************
     623         2098 :    SUBROUTINE read_element_from_file(elem, tmc_ana, conf_nr, stat)
     624              :       TYPE(tree_type), POINTER                           :: elem
     625              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
     626              :       INTEGER                                            :: conf_nr, stat
     627              : 
     628              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_element_from_file'
     629              : 
     630              :       INTEGER                                            :: conf_nr_old, handle, i_tmp
     631              :       LOGICAL                                            :: files_conf_missmatch
     632              : 
     633         1049 :       stat = TMC_STATUS_OK
     634         1049 :       conf_nr_old = conf_nr
     635         1049 :       files_conf_missmatch = .FALSE.
     636              : 
     637         1049 :       CPASSERT(ASSOCIATED(elem))
     638         1049 :       CPASSERT(ASSOCIATED(tmc_ana))
     639         1049 :       CPASSERT(ASSOCIATED(tmc_ana%atoms))
     640              : 
     641              :       ! start the timing
     642         1049 :       CALL timeset(routineN, handle)
     643              : 
     644              :       ! read the coordinates
     645         1049 :       IF (tmc_ana%id_traj > 0) THEN
     646         1049 :          i_tmp = conf_nr_old
     647              :          CALL read_pos_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
     648         1049 :                                  conf_nr=i_tmp)
     649         1049 :          IF (stat == TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
     650              :             CALL cp_warn(__LOCATION__, &
     651              :                          'end of position file reached at line '// &
     652              :                          cp_to_string(REAL(tmc_ana%lc_traj, KIND=dp))//", last element "// &
     653           18 :                          cp_to_string(tmc_ana%last_elem%nr))
     654              :          ELSE
     655         1031 :             CPASSERT(i_tmp > conf_nr_old)
     656         1031 :             conf_nr = i_tmp
     657         1031 :             elem%nr = i_tmp
     658              :          END IF
     659              :       END IF
     660              : 
     661              :       ! read the forces
     662              :       ! TODO if necessary
     663              : 
     664              :       ! read the dipoles file
     665         1049 :       IF (tmc_ana%id_dip > 0 .AND. stat == TMC_STATUS_OK) THEN
     666            0 :          i_tmp = conf_nr_old
     667              :          search_conf_dip: DO
     668              :             CALL read_dipole_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
     669            0 :                                        conf_nr=i_tmp)
     670            0 :             IF (stat == TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
     671              :                CALL cp_warn(__LOCATION__, &
     672              :                             'end of dipole file reached at line'// &
     673            0 :                             cp_to_string(REAL(tmc_ana%lc_dip, KIND=dp)))
     674            0 :                EXIT search_conf_dip
     675              :             END IF
     676              :             ! check consitence with pos file
     677            0 :             IF (tmc_ana%id_traj > 0) THEN
     678            0 :                IF (i_tmp == conf_nr) THEN
     679              :                   files_conf_missmatch = .FALSE.
     680              :                   EXIT search_conf_dip
     681              :                ELSE
     682              :                   ! the configuration numbering differ from the position file,
     683              :                   !  but we keep on searching for the correct configuration
     684              :                   files_conf_missmatch = .TRUE.
     685              :                END IF
     686              :                ! if no pos file, just take the next conf
     687            0 :             ELSE IF (i_tmp > conf_nr_old) THEN
     688            0 :                conf_nr = i_tmp
     689            0 :                elem%nr = i_tmp
     690            0 :                EXIT search_conf_dip
     691              :             END IF
     692              :          END DO search_conf_dip
     693              :       END IF
     694              : 
     695              :       ! read the cell file
     696         1049 :       IF (tmc_ana%id_cell > 0 .AND. stat == TMC_STATUS_OK) THEN
     697              :          search_conf_cell: DO
     698              :             CALL read_cell_from_file(elem=elem, tmc_ana=tmc_ana, stat=stat, &
     699         1206 :                                      conf_nr=i_tmp)
     700         1206 :             IF (stat == TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
     701              :                CALL cp_warn(__LOCATION__, &
     702              :                             'end of cell file reached at line at line'// &
     703            0 :                             cp_to_string(REAL(tmc_ana%lc_cell, KIND=dp)))
     704            0 :                EXIT search_conf_cell
     705              :             END IF
     706              :             ! check consitence with pos file
     707         1206 :             IF (tmc_ana%id_traj > 0) THEN
     708         1206 :                IF (i_tmp == conf_nr) THEN
     709              :                   files_conf_missmatch = .FALSE.
     710              :                   EXIT search_conf_cell
     711              :                ELSE
     712              :                   ! the configuration numbering differ from the position file,
     713              :                   !  but we keep on searching for the correct configuration
     714              :                   files_conf_missmatch = .TRUE.
     715              :                END IF
     716              :                ! if no pos file, just take the next conf
     717            0 :             ELSE IF (i_tmp > conf_nr_old) THEN
     718            0 :                conf_nr = i_tmp
     719            0 :                elem%nr = i_tmp
     720            0 :                EXIT search_conf_cell
     721              :             END IF
     722              :          END DO search_conf_cell
     723              : 
     724              :       END IF
     725              : 
     726              :       ! write the different energies
     727              :       ! TODO if necessary
     728              : 
     729         1049 :       IF (files_conf_missmatch) THEN
     730              :          CALL cp_warn(__LOCATION__, &
     731              :                       'there is a missmatch in the configuration numbering. '// &
     732              :                       "Read number of lines (pos|cell|dip)"// &
     733              :                       cp_to_string(tmc_ana%lc_traj)//"|"// &
     734              :                       cp_to_string(tmc_ana%lc_cell)//"|"// &
     735            0 :                       cp_to_string(tmc_ana%lc_dip))
     736              :       END IF
     737              : 
     738              :       ! end the timing
     739         1049 :       CALL timestop(handle)
     740         1049 :    END SUBROUTINE read_element_from_file
     741              : 
     742              : ! **************************************************************************************************
     743              : !> \brief search for the next configurational position in file
     744              : !> \param elem actual tree element to be read
     745              : !> \param tmc_ana ...
     746              : !> \param stat ...
     747              : !> \param conf_nr Markov chain element number
     748              : !>        (input the old number and read only if conf nr from file is greater
     749              : !> \param header_info ...
     750              : !> \author Mandes 03.2013
     751              : ! **************************************************************************************************
     752         2098 :    SUBROUTINE read_pos_from_file(elem, tmc_ana, stat, conf_nr, header_info)
     753              :       TYPE(tree_type), POINTER                           :: elem
     754              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
     755              :       INTEGER                                            :: stat, conf_nr
     756              :       CHARACTER(LEN=*), OPTIONAL                         :: header_info
     757              : 
     758              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_pos_from_file'
     759              : 
     760              :       CHARACTER(LEN=default_string_length)               :: c_tmp
     761              :       INTEGER                                            :: handle, i, i_tmp, status
     762              : 
     763         1049 :       stat = TMC_STATUS_FAILED
     764              : 
     765            0 :       CPASSERT(ASSOCIATED(elem))
     766         1049 :       CPASSERT(ASSOCIATED(elem%pos))
     767         1049 :       CPASSERT(ASSOCIATED(tmc_ana))
     768         1049 :       CPASSERT(tmc_ana%id_traj > 0)
     769              : 
     770              :       ! start the timing
     771         1049 :       CALL timeset(routineN, handle)
     772              : 
     773              :       search_next_conf: DO
     774         6105 :          c_tmp(:) = " "
     775         6105 :          tmc_ana%lc_traj = tmc_ana%lc_traj + 1
     776         6105 :          READ (tmc_ana%id_traj, '(A)', IOSTAT=status) c_tmp(:)
     777         6105 :          IF (status > 0) THEN
     778              :             CALL cp_abort(__LOCATION__, &
     779              :                           "configuration header read error at line: "// &
     780            0 :                           cp_to_string(tmc_ana%lc_traj)//": "//c_tmp)
     781              :          END IF
     782         6105 :          IF (status < 0) THEN ! end of file reached
     783           18 :             stat = TMC_STATUS_WAIT_FOR_NEW_TASK
     784           18 :             EXIT search_next_conf
     785              :          END IF
     786         6087 :          IF (INDEX(c_tmp, "=") > 0) THEN
     787         1206 :             READ (c_tmp(INDEX(c_tmp, "=") + 1:), *, IOSTAT=status) i_tmp ! read the configuration number
     788         1206 :             IF (status /= 0) THEN
     789              :                CALL cp_abort(__LOCATION__, &
     790              :                              "configuration header read error (for conf nr) at line: "// &
     791            0 :                              cp_to_string(tmc_ana%lc_traj))
     792              :             END IF
     793         1206 :             IF (i_tmp > conf_nr) THEN
     794              :                ! TODO we could also read the energy ...
     795         1031 :                conf_nr = i_tmp
     796         1031 :                IF (PRESENT(header_info)) header_info = c_tmp
     797         1031 :                stat = TMC_STATUS_OK
     798         1031 :                EXIT search_next_conf
     799              :             END IF
     800              :          END IF
     801              :       END DO search_next_conf
     802              : 
     803         1049 :       IF (stat == TMC_STATUS_OK) THEN
     804        22682 :          pos_loop: DO i = 1, SIZE(elem%pos), tmc_ana%dim_per_elem
     805        21651 :             tmc_ana%lc_traj = tmc_ana%lc_traj + 1
     806              :             READ (tmc_ana%id_traj, FMT="(A4,1X,1000F20.10)", IOSTAT=status) &
     807        86604 :                c_tmp, elem%pos(i:i + tmc_ana%dim_per_elem - 1)
     808        22682 :             IF (status /= 0) THEN
     809              :                CALL cp_abort(__LOCATION__, &
     810              :                              "configuration pos read error at line: "// &
     811            0 :                              cp_to_string(tmc_ana%lc_traj))
     812              :             END IF
     813              :          END DO pos_loop
     814        65984 :          elem%pos(:) = elem%pos(:)/au2a
     815              :       END IF
     816              : 
     817              :       ! end the timing
     818         1049 :       CALL timestop(handle)
     819         1049 :    END SUBROUTINE read_pos_from_file
     820              : 
     821              : ! **************************************************************************************************
     822              : !> \brief search for the dipole entry
     823              : !> \param elem actual tree element to be read
     824              : !> \param tmc_ana ...
     825              : !> \param stat ...
     826              : !> \param conf_nr Markov chain element number
     827              : !>        (input the old number and read only if conf nr from file is greater
     828              : !> \author Mandes 03.2013
     829              : ! **************************************************************************************************
     830            0 :    SUBROUTINE read_dipole_from_file(elem, tmc_ana, stat, conf_nr)
     831              :       TYPE(tree_type), POINTER                           :: elem
     832              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
     833              :       INTEGER                                            :: stat, conf_nr
     834              : 
     835              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_dipole_from_file'
     836              : 
     837              :       CHARACTER(LEN=250)                                 :: c_tmp
     838              :       INTEGER                                            :: handle, status
     839              : 
     840            0 :       stat = TMC_STATUS_FAILED
     841              : 
     842            0 :       CPASSERT(ASSOCIATED(elem))
     843            0 :       CPASSERT(ASSOCIATED(elem%dipole))
     844            0 :       CPASSERT(ASSOCIATED(tmc_ana))
     845            0 :       CPASSERT(tmc_ana%id_dip > 0)
     846              : 
     847              :       ! start the timing
     848            0 :       CALL timeset(routineN, handle)
     849            0 :       tmc_ana%lc_dip = tmc_ana%lc_dip + 1
     850            0 :       READ (tmc_ana%id_dip, FMT="(A)", IOSTAT=status) c_tmp
     851            0 :       IF (status == 0) THEN
     852              :          ! skip the initial line (header)
     853            0 :          IF (INDEX(c_tmp, "#") > 0) THEN
     854            0 :             tmc_ana%lc_dip = tmc_ana%lc_dip + 1
     855            0 :             READ (tmc_ana%id_dip, FMT="(A)", IOSTAT=status) c_tmp
     856              :          END IF
     857              :       END IF
     858            0 :       IF (status == 0) THEN
     859              :          READ (c_tmp, FMT="(I8,10F20.10)", IOSTAT=status) &
     860            0 :             conf_nr, elem%dipole(:)
     861              :       END IF
     862            0 :       IF (status == 0) THEN ! success
     863            0 :          stat = TMC_STATUS_OK
     864            0 :       ELSE IF (status < 0) THEN ! end of file reached
     865            0 :          stat = TMC_STATUS_WAIT_FOR_NEW_TASK
     866              :       ELSE
     867              :          IF (status /= 0) THEN
     868            0 :             CPWARN("configuration dipole read error at line: "//cp_to_string(tmc_ana%lc_dip))
     869              :          END IF
     870            0 :          stat = TMC_STATUS_FAILED
     871              :       END IF
     872              : 
     873              :       ! end the timing
     874            0 :       CALL timestop(handle)
     875            0 :    END SUBROUTINE read_dipole_from_file
     876              : 
     877              : ! **************************************************************************************************
     878              : !> \brief search for the cell entry
     879              : !> \param elem actual tree element to be read
     880              : !> \param tmc_ana ...
     881              : !> \param stat ...
     882              : !> \param conf_nr Markov chain element number
     883              : !>        (input the old number and read only if conf nr from file is greater
     884              : !> \author Mandes 03.2013
     885              : ! **************************************************************************************************
     886         2412 :    SUBROUTINE read_cell_from_file(elem, tmc_ana, stat, conf_nr)
     887              :       TYPE(tree_type), POINTER                           :: elem
     888              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
     889              :       INTEGER                                            :: stat, conf_nr
     890              : 
     891              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'read_cell_from_file'
     892              : 
     893              :       CHARACTER(LEN=250)                                 :: c_tmp
     894              :       INTEGER                                            :: handle, status
     895              :       REAL(KIND=dp)                                      :: r_tmp
     896              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: hmat
     897              : 
     898         1206 :       stat = TMC_STATUS_FAILED
     899              : 
     900         1206 :       CPASSERT(ASSOCIATED(elem))
     901         1206 :       CPASSERT(ASSOCIATED(tmc_ana))
     902         1206 :       CPASSERT(ASSOCIATED(tmc_ana%cell))
     903         1206 :       CPASSERT(tmc_ana%id_cell > 0)
     904              : 
     905              :       ! start the timing
     906         1206 :       CALL timeset(routineN, handle)
     907              : 
     908         1206 :       tmc_ana%lc_cell = tmc_ana%lc_cell + 1
     909         1206 :       READ (tmc_ana%id_cell, FMT="(A)", IOSTAT=status) c_tmp
     910         1206 :       IF (status == 0) THEN
     911              :          ! skip the initial line (header)
     912         1206 :          IF (INDEX(c_tmp, "#") > 0) THEN
     913           18 :             tmc_ana%lc_cell = tmc_ana%lc_cell + 1
     914           18 :             READ (tmc_ana%id_cell, FMT="(A)", IOSTAT=status) c_tmp
     915              :          END IF
     916              :       END IF
     917         1206 :       IF (status == 0) THEN
     918         1206 :          READ (c_tmp, FMT="(I8,9(1X,F19.10),1X,F24.10)", IOSTAT=status) conf_nr, &
     919         2412 :             hmat(:, :), r_tmp
     920              :       END IF
     921         1206 :       IF (status < 0) THEN ! end of file reached
     922            0 :          stat = TMC_STATUS_WAIT_FOR_NEW_TASK
     923         1206 :       ELSE IF (status > 0) THEN
     924              :          IF (status /= 0) THEN
     925            0 :             CPABORT("configuration cell read error at line: "//cp_to_string(tmc_ana%lc_cell))
     926              :          END IF
     927            0 :          stat = TMC_STATUS_FAILED
     928              :       ELSE
     929         1206 :          IF (elem%nr < 0) elem%nr = conf_nr
     930        15678 :          hmat(:, :) = hmat(:, :)/au2a
     931              :          ! get the box scaling
     932              :          CALL get_cell_scaling(cell=tmc_ana%cell, scaled_hmat=hmat, &
     933         1206 :                                box_scale=elem%box_scale)
     934         1206 :          stat = TMC_STATUS_OK
     935              :       END IF
     936              :       ! end the timing
     937         1206 :       CALL timestop(handle)
     938         1206 :    END SUBROUTINE read_cell_from_file
     939              : 
     940              :    !----------------------------------------------------------------------------
     941              :    ! get the configurations from file and calc
     942              :    !----------------------------------------------------------------------------
     943              : 
     944              : ! **************************************************************************************************
     945              : !> \brief opens the files for reading configurations data to analyze
     946              : !> \param tmc_ana ...
     947              : !> \param stat ...
     948              : !> \param dir_ind ...
     949              : !> \param
     950              : !> \author Mandes 02.2013
     951              : ! **************************************************************************************************
     952           36 :    SUBROUTINE analyse_files_open(tmc_ana, stat, dir_ind)
     953              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
     954              :       INTEGER                                            :: stat
     955              :       INTEGER, OPTIONAL                                  :: dir_ind
     956              : 
     957              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'analyse_files_open'
     958              : 
     959              :       CHARACTER(LEN=default_path_length)                 :: dir_name, file_name_act, file_name_temp
     960              :       INTEGER                                            :: handle
     961              :       LOGICAL                                            :: file_exists
     962              : 
     963           18 :       CPASSERT(ASSOCIATED(tmc_ana))
     964              : 
     965           18 :       stat = TMC_STATUS_WAIT_FOR_NEW_TASK
     966              : 
     967              :       ! start the timing
     968           18 :       CALL timeset(routineN, handle)
     969              : 
     970           18 :       IF (PRESENT(dir_ind)) THEN
     971           18 :          CPASSERT(ASSOCIATED(tmc_ana%dirs))
     972           18 :          CPASSERT(dir_ind > 0)
     973           18 :          CPASSERT(dir_ind <= SIZE(tmc_ana%dirs))
     974              : 
     975           18 :          IF (INDEX(tmc_ana%dirs(dir_ind), "/", BACK=.TRUE.) == &
     976              :              LEN_TRIM(tmc_ana%dirs(dir_ind))) THEN
     977           18 :             dir_name = TRIM(tmc_ana%dirs(dir_ind))
     978              :          ELSE
     979            0 :             dir_name = TRIM(tmc_ana%dirs(dir_ind))//"/"
     980              :          END IF
     981              :       ELSE
     982            0 :          dir_name = "./"
     983              :       END IF
     984              : 
     985              :       ! open the files
     986              :       file_name_temp = expand_file_name_temp( &
     987              :                        file_name=tmc_default_trajectory_file_name, &
     988           18 :                        rvalue=tmc_ana%temperature)
     989              :       ! position file
     990           18 :       IF (tmc_ana%costum_pos_file_name /= "") THEN
     991            0 :          file_name_act = TRIM(dir_name)//tmc_ana%costum_pos_file_name
     992              :       ELSE
     993              :          file_name_act = TRIM(dir_name)// &
     994           18 :                          expand_file_name_ending(file_name_temp, "xyz")
     995              :       END IF
     996           18 :       INQUIRE (FILE=file_name_act, EXIST=file_exists)
     997           18 :       IF (file_exists) THEN
     998              :          CALL open_file(file_name=file_name_act, file_status="OLD", &
     999           18 :                         file_action="READ", unit_number=tmc_ana%id_traj)
    1000           18 :          WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
    1001           36 :             "read xyz file", TRIM(file_name_act)
    1002              :       END IF
    1003              : 
    1004              :       ! cell file
    1005           18 :       IF (tmc_ana%costum_cell_file_name /= "") THEN
    1006            0 :          file_name_act = TRIM(dir_name)//tmc_ana%costum_cell_file_name
    1007              :       ELSE
    1008              :          file_name_act = TRIM(dir_name)// &
    1009           18 :                          expand_file_name_ending(file_name_temp, "cell")
    1010              :       END IF
    1011           18 :       INQUIRE (FILE=file_name_act, EXIST=file_exists)
    1012           18 :       IF (file_exists) THEN
    1013              :          CALL open_file(file_name=file_name_act, file_status="OLD", &
    1014           18 :                         file_action="READ", unit_number=tmc_ana%id_cell)
    1015           18 :          WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
    1016           36 :             "read cell file", TRIM(file_name_act)
    1017              :       END IF
    1018              : 
    1019              :       ! dipole file
    1020           18 :       IF (tmc_ana%costum_dip_file_name /= "") THEN
    1021           18 :          file_name_act = TRIM(dir_name)//tmc_ana%costum_dip_file_name
    1022              :       ELSE
    1023              :          file_name_act = TRIM(dir_name)// &
    1024            0 :                          expand_file_name_ending(file_name_temp, "dip")
    1025              :       END IF
    1026           18 :       INQUIRE (FILE=file_name_act, EXIST=file_exists)
    1027           18 :       IF (file_exists) THEN
    1028              :          CALL open_file(file_name=file_name_act, file_status="OLD", &
    1029            0 :                         file_action="READ", unit_number=tmc_ana%id_dip)
    1030            0 :          WRITE (tmc_ana%io_unit, FMT='(T2,A,"| ",A,T41,A40)') "TMC_ANA", &
    1031            0 :             "read dip file", TRIM(file_name_act)
    1032              :       END IF
    1033              : 
    1034           18 :       IF (tmc_ana%id_traj > 0 .OR. tmc_ana%id_cell > 0 .OR. &
    1035              :           tmc_ana%id_dip > 0) THEN
    1036           18 :          stat = TMC_STATUS_OK
    1037              :       ELSE
    1038              :          CALL cp_warn(__LOCATION__, &
    1039              :                       "There is no file to open for temperature "//cp_to_string(tmc_ana%temperature)// &
    1040            0 :                       "K in directory "//TRIM(dir_name))
    1041              :       END IF
    1042              :       ! end the timing
    1043           18 :       CALL timestop(handle)
    1044           18 :    END SUBROUTINE analyse_files_open
    1045              : 
    1046              : ! **************************************************************************************************
    1047              : !> \brief close the files for reading configurations data to analyze
    1048              : !> \param tmc_ana ...
    1049              : !> \param
    1050              : !> \author Mandes 02.2013
    1051              : ! **************************************************************************************************
    1052           36 :    SUBROUTINE analyse_files_close(tmc_ana)
    1053              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
    1054              : 
    1055              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'analyse_files_close'
    1056              : 
    1057              :       INTEGER                                            :: handle
    1058              : 
    1059           18 :       CPASSERT(ASSOCIATED(tmc_ana))
    1060              : 
    1061              :       ! start the timing
    1062           18 :       CALL timeset(routineN, handle)
    1063              : 
    1064              :       ! position file
    1065           18 :       IF (tmc_ana%id_traj > 0) CALL close_file(unit_number=tmc_ana%id_traj)
    1066              : 
    1067              :       ! cell file
    1068           18 :       IF (tmc_ana%id_cell > 0) CALL close_file(unit_number=tmc_ana%id_cell)
    1069              : 
    1070              :       ! dipole file
    1071           18 :       IF (tmc_ana%id_dip > 0) CALL close_file(unit_number=tmc_ana%id_dip)
    1072              : 
    1073              :       ! end the timing
    1074           18 :       CALL timestop(handle)
    1075           18 :    END SUBROUTINE analyse_files_close
    1076              : 
    1077              : END MODULE tmc_file_io
        

Generated by: LCOV version 2.0-1