LCOV - code coverage report
Current view: top level - src/tmc - tmc_analysis.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 72.7 % 719 523
Test Date: 2025-07-25 12:55:17 Functions: 85.0 % 20 17

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief module analyses element of the TMC tree element structure
      10              : !>        e.g. density, radial distribution function, dipole correlation,...
      11              : !> \par History
      12              : !>      02.2013 created [Mandes Schoenherr]
      13              : !> \author Mandes
      14              : ! **************************************************************************************************
      15              : 
      16              : MODULE tmc_analysis
      17              :    USE cell_types,                      ONLY: cell_type,&
      18              :                                               get_cell,&
      19              :                                               pbc
      20              :    USE cp_files,                        ONLY: close_file,&
      21              :                                               open_file
      22              :    USE cp_log_handling,                 ONLY: cp_to_string
      23              :    USE force_fields_input,              ONLY: read_chrg_section
      24              :    USE input_section_types,             ONLY: section_vals_get,&
      25              :                                               section_vals_get_subs_vals,&
      26              :                                               section_vals_type,&
      27              :                                               section_vals_val_get
      28              :    USE kinds,                           ONLY: default_path_length,&
      29              :                                               default_string_length,&
      30              :                                               dp
      31              :    USE mathconstants,                   ONLY: pi
      32              :    USE mathlib,                         ONLY: diag
      33              :    USE physcon,                         ONLY: a_mass,&
      34              :                                               au2a => angstrom,&
      35              :                                               boltzmann,&
      36              :                                               joule,&
      37              :                                               massunit
      38              :    USE tmc_analysis_types,              ONLY: &
      39              :         ana_type_default, ana_type_ice, ana_type_sym_xyz, atom_pairs_type, dipole_moment_type, &
      40              :         pair_correl_type, search_pair_in_list, tmc_ana_density_create, tmc_ana_density_file_name, &
      41              :         tmc_ana_dipole_analysis_create, tmc_ana_dipole_moment_create, tmc_ana_displacement_create, &
      42              :         tmc_ana_env_create, tmc_ana_pair_correl_create, tmc_ana_pair_correl_file_name, &
      43              :         tmc_analysis_env
      44              :    USE tmc_calculations,                ONLY: get_scaled_cell,&
      45              :                                               nearest_distance
      46              :    USE tmc_file_io,                     ONLY: analyse_files_close,&
      47              :                                               analyse_files_open,&
      48              :                                               expand_file_name_char,&
      49              :                                               expand_file_name_temp,&
      50              :                                               read_element_from_file,&
      51              :                                               write_dipoles_in_file
      52              :    USE tmc_stati,                       ONLY: TMC_STATUS_OK,&
      53              :                                               TMC_STATUS_WAIT_FOR_NEW_TASK,&
      54              :                                               tmc_default_restart_in_file_name,&
      55              :                                               tmc_default_restart_out_file_name,&
      56              :                                               tmc_default_trajectory_file_name,&
      57              :                                               tmc_default_unspecified_name
      58              :    USE tmc_tree_build,                  ONLY: allocate_new_sub_tree_node,&
      59              :                                               deallocate_sub_tree_node
      60              :    USE tmc_tree_types,                  ONLY: read_subtree_elem_unformated,&
      61              :                                               tree_type,&
      62              :                                               write_subtree_elem_unformated
      63              :    USE tmc_types,                       ONLY: tmc_atom_type,&
      64              :                                               tmc_param_type
      65              : #include "../base/base_uses.f90"
      66              : 
      67              :    IMPLICIT NONE
      68              : 
      69              :    PRIVATE
      70              : 
      71              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'tmc_analysis'
      72              : 
      73              :    PUBLIC :: tmc_read_ana_input
      74              :    PUBLIC :: analysis_init, do_tmc_analysis, analyze_file_configurations, finalize_tmc_analysis
      75              :    PUBLIC :: analysis_restart_print, analysis_restart_read
      76              : 
      77              : CONTAINS
      78              : 
      79              : ! **************************************************************************************************
      80              : !> \brief creates a new para environment for tmc analysis
      81              : !> \param tmc_ana_section ...
      82              : !> \param tmc_ana TMC analysis environment
      83              : !> \author Mandes 02.2013
      84              : ! **************************************************************************************************
      85           54 :    SUBROUTINE tmc_read_ana_input(tmc_ana_section, tmc_ana)
      86              :       TYPE(section_vals_type), POINTER                   :: tmc_ana_section
      87              :       TYPE(tmc_analysis_env), POINTER                    :: tmc_ana
      88              : 
      89              :       CHARACTER(LEN=default_path_length)                 :: c_tmp
      90           18 :       CHARACTER(LEN=default_string_length), POINTER      :: charge_atm(:)
      91              :       INTEGER                                            :: i_tmp, ntot
      92              :       INTEGER, DIMENSION(3)                              :: nr_bins
      93           18 :       INTEGER, DIMENSION(:), POINTER                     :: i_arr_tmp
      94              :       LOGICAL                                            :: explicit, explicit_key, flag
      95           18 :       REAL(KIND=dp), POINTER                             :: charge(:)
      96              :       TYPE(section_vals_type), POINTER                   :: tmp_section
      97              : 
      98           18 :       NULLIFY (tmp_section, charge_atm, i_arr_tmp, charge)
      99              : 
     100            0 :       CPASSERT(ASSOCIATED(tmc_ana_section))
     101           18 :       CPASSERT(.NOT. ASSOCIATED(tmc_ana))
     102              : 
     103           18 :       CALL section_vals_get(tmc_ana_section, explicit=explicit)
     104           18 :       IF (explicit) THEN
     105           18 :          CALL tmc_ana_env_create(tmc_ana=tmc_ana)
     106              :          ! restarting
     107           18 :          CALL section_vals_val_get(tmc_ana_section, "RESTART", l_val=tmc_ana%restart)
     108              :          ! file name prefix
     109              :          CALL section_vals_val_get(tmc_ana_section, "PREFIX_ANA_FILES", &
     110           18 :                                    c_val=tmc_ana%out_file_prefix)
     111           18 :          IF (tmc_ana%out_file_prefix .NE. "") THEN
     112            0 :             tmc_ana%out_file_prefix = TRIM(tmc_ana%out_file_prefix)//"_"
     113              :          END IF
     114              : 
     115              :          ! density calculation
     116           18 :          CALL section_vals_val_get(tmc_ana_section, "DENSITY", explicit=explicit_key)
     117           18 :          IF (explicit_key) THEN
     118            9 :             CALL section_vals_val_get(tmc_ana_section, "DENSITY", i_vals=i_arr_tmp)
     119              : 
     120            9 :             IF (SIZE(i_arr_tmp(:)) .EQ. 3) THEN
     121           36 :                IF (ANY(i_arr_tmp(:) .LE. 0)) &
     122              :                   CALL cp_abort(__LOCATION__, "The amount of intervals in each "// &
     123            0 :                                 "direction has to be greater than 0.")
     124           36 :                nr_bins(:) = i_arr_tmp(:)
     125            0 :             ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 1) THEN
     126            0 :                IF (ANY(i_arr_tmp(:) .LE. 0)) &
     127            0 :                   CPABORT("The amount of intervals has to be greater than 0.")
     128            0 :                nr_bins(:) = i_arr_tmp(1)
     129            0 :             ELSE IF (SIZE(i_arr_tmp(:)) .EQ. 0) THEN
     130            0 :                nr_bins(:) = 1
     131              :             ELSE
     132            0 :                CPABORT("unknown amount of dimensions for the binning.")
     133              :             END IF
     134            9 :             CALL tmc_ana_density_create(tmc_ana%density_3d, nr_bins)
     135              :          END IF
     136              : 
     137              :          ! radial distribution function calculation
     138           18 :          CALL section_vals_val_get(tmc_ana_section, "G_R", explicit=explicit_key)
     139           18 :          IF (explicit_key) THEN
     140            9 :             CALL section_vals_val_get(tmc_ana_section, "G_R", i_val=i_tmp)
     141              :             CALL tmc_ana_pair_correl_create(ana_pair_correl=tmc_ana%pair_correl, &
     142            9 :                                             nr_bins=i_tmp)
     143              :          END IF
     144              : 
     145              :          ! radial distribution function calculation
     146           18 :          CALL section_vals_val_get(tmc_ana_section, "CLASSICAL_DIPOLE_MOMENTS", explicit=explicit_key)
     147           18 :          IF (explicit_key) THEN
     148              :             ! charges for dipoles needed
     149            9 :             tmp_section => section_vals_get_subs_vals(tmc_ana_section, "CHARGE")
     150            9 :             CALL section_vals_get(tmp_section, explicit=explicit, n_repetition=i_tmp)
     151            9 :             IF (explicit) THEN
     152            9 :                ntot = 0
     153           27 :                ALLOCATE (charge_atm(i_tmp))
     154           27 :                ALLOCATE (charge(i_tmp))
     155            9 :                CALL read_chrg_section(charge_atm, charge, tmp_section, ntot)
     156              :             ELSE
     157              :                CALL cp_abort(__LOCATION__, &
     158              :                              "to calculate the classical cell dipole moment "// &
     159            0 :                              "the charges has to be specified")
     160              :             END IF
     161              : 
     162              :             CALL tmc_ana_dipole_moment_create(tmc_ana%dip_mom, charge_atm, charge, &
     163            9 :                                               tmc_ana%dim_per_elem)
     164              : 
     165            9 :             IF (ASSOCIATED(charge_atm)) DEALLOCATE (charge_atm)
     166            9 :             IF (ASSOCIATED(charge)) DEALLOCATE (charge)
     167              :          END IF
     168              : 
     169              :          ! dipole moment analysis
     170           18 :          CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", explicit=explicit_key)
     171           18 :          IF (explicit_key) THEN
     172            0 :             CALL tmc_ana_dipole_analysis_create(tmc_ana%dip_ana)
     173            0 :             CALL section_vals_val_get(tmc_ana_section, "DIPOLE_ANALYSIS", c_val=c_tmp)
     174            0 :             SELECT CASE (TRIM(c_tmp))
     175              :             CASE (TRIM(tmc_default_unspecified_name))
     176            0 :                tmc_ana%dip_ana%ana_type = ana_type_default
     177              :             CASE ("ICE")
     178            0 :                tmc_ana%dip_ana%ana_type = ana_type_ice
     179              :             CASE ("SYM_XYZ")
     180            0 :                tmc_ana%dip_ana%ana_type = ana_type_sym_xyz
     181              :             CASE DEFAULT
     182            0 :                CPWARN('unknown analysis type "'//TRIM(c_tmp)//'" specified. Set to default.')
     183            0 :                tmc_ana%dip_ana%ana_type = ana_type_default
     184              :             END SELECT
     185              :          END IF
     186              : 
     187              :       END IF
     188              : 
     189              :       ! cell displacement (deviation)
     190           18 :       CALL section_vals_val_get(tmc_ana_section, "DEVIATION", l_val=flag)
     191           18 :       IF (flag) THEN
     192              :          CALL tmc_ana_displacement_create(ana_disp=tmc_ana%displace, &
     193            9 :                                           dim_per_elem=tmc_ana%dim_per_elem)
     194              :       END IF
     195           18 :    END SUBROUTINE tmc_read_ana_input
     196              : 
     197              : ! **************************************************************************************************
     198              : !> \brief initialize all the necessarry analysis structures
     199              : !> \param ana_env ...
     200              : !> \param nr_dim dimension of the pos, frc etc. array
     201              : !> \author Mandes 02.2013
     202              : ! **************************************************************************************************
     203           18 :    SUBROUTINE analysis_init(ana_env, nr_dim)
     204              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     205              :       INTEGER                                            :: nr_dim
     206              : 
     207              :       CHARACTER(LEN=default_path_length)                 :: tmp_cell_file, tmp_dip_file, tmp_pos_file
     208              : 
     209           18 :       CPASSERT(ASSOCIATED(ana_env))
     210           18 :       CPASSERT(nr_dim > 0)
     211              : 
     212           18 :       ana_env%nr_dim = nr_dim
     213              : 
     214              :       ! save file names
     215           18 :       tmp_pos_file = ana_env%costum_pos_file_name
     216           18 :       tmp_cell_file = ana_env%costum_cell_file_name
     217           18 :       tmp_dip_file = ana_env%costum_dip_file_name
     218              : 
     219              :       ! unset all filenames
     220           18 :       ana_env%costum_pos_file_name = tmc_default_unspecified_name
     221           18 :       ana_env%costum_cell_file_name = tmc_default_unspecified_name
     222           18 :       ana_env%costum_dip_file_name = tmc_default_unspecified_name
     223              : 
     224              :       ! set the necessary files for ...
     225              :       ! density
     226           18 :       IF (ASSOCIATED(ana_env%density_3d)) THEN
     227            9 :          ana_env%costum_pos_file_name = tmp_pos_file
     228            9 :          ana_env%costum_cell_file_name = tmp_cell_file
     229              :       END IF
     230              :       ! pair correlation
     231           18 :       IF (ASSOCIATED(ana_env%pair_correl)) THEN
     232            9 :          ana_env%costum_pos_file_name = tmp_pos_file
     233            9 :          ana_env%costum_cell_file_name = tmp_cell_file
     234              :       END IF
     235              :       ! dipole moment
     236           18 :       IF (ASSOCIATED(ana_env%dip_mom)) THEN
     237            9 :          ana_env%costum_pos_file_name = tmp_pos_file
     238            9 :          ana_env%costum_cell_file_name = tmp_cell_file
     239              :       END IF
     240              :       ! dipole analysis
     241           18 :       IF (ASSOCIATED(ana_env%dip_ana)) THEN
     242            0 :          ana_env%costum_pos_file_name = tmp_pos_file
     243            0 :          ana_env%costum_cell_file_name = tmp_cell_file
     244            0 :          ana_env%costum_dip_file_name = tmp_dip_file
     245              :       END IF
     246              :       ! deviation / displacement
     247           18 :       IF (ASSOCIATED(ana_env%displace)) THEN
     248            9 :          ana_env%costum_pos_file_name = tmp_pos_file
     249            9 :          ana_env%costum_cell_file_name = tmp_cell_file
     250              :       END IF
     251              : 
     252              :       ! init radial distribution function
     253           18 :       IF (ASSOCIATED(ana_env%pair_correl)) &
     254              :          CALL ana_pair_correl_init(ana_pair_correl=ana_env%pair_correl, &
     255            9 :                                    atoms=ana_env%atoms, cell=ana_env%cell)
     256              :       ! init classical dipole moment calculations
     257           18 :       IF (ASSOCIATED(ana_env%dip_mom)) &
     258              :          CALL ana_dipole_moment_init(ana_dip_mom=ana_env%dip_mom, &
     259            9 :                                      atoms=ana_env%atoms)
     260           18 :    END SUBROUTINE analysis_init
     261              : 
     262              : ! **************************************************************************************************
     263              : !> \brief print analysis restart file
     264              : !> \param ana_env ...
     265              : !> \param
     266              : !> \author Mandes 02.2013
     267              : ! **************************************************************************************************
     268           18 :    SUBROUTINE analysis_restart_print(ana_env)
     269              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     270              : 
     271              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_tmp, &
     272              :                                                             restart_file_name
     273              :       INTEGER                                            :: file_ptr
     274              :       LOGICAL                                            :: l_tmp
     275              : 
     276           18 :       CPASSERT(ASSOCIATED(ana_env))
     277           18 :       CPASSERT(ASSOCIATED(ana_env%last_elem))
     278           18 :       IF (.NOT. ana_env%restart) RETURN
     279              : 
     280            6 :       WRITE (file_name, FMT='(I9.9)') ana_env%last_elem%nr
     281              :       file_name_tmp = TRIM(expand_file_name_temp(expand_file_name_char( &
     282              :                                                  TRIM(ana_env%out_file_prefix)// &
     283              :                                                  tmc_default_restart_out_file_name, &
     284            6 :                                                  "ana"), ana_env%temperature))
     285              :       restart_file_name = expand_file_name_char(file_name_tmp, &
     286            6 :                                                 file_name)
     287              :       CALL open_file(file_name=restart_file_name, file_status="REPLACE", &
     288              :                      file_action="WRITE", file_form="UNFORMATTED", &
     289            6 :                      unit_number=file_ptr)
     290            6 :       WRITE (file_ptr) ana_env%temperature
     291            6 :       CALL write_subtree_elem_unformated(ana_env%last_elem, file_ptr)
     292              : 
     293              :       ! first mention the different kind of anlysis types initialized
     294              :       ! then the variables for each calculation type
     295            6 :       l_tmp = ASSOCIATED(ana_env%density_3d)
     296            6 :       WRITE (file_ptr) l_tmp
     297            6 :       IF (l_tmp) THEN
     298            6 :          WRITE (file_ptr) ana_env%density_3d%conf_counter, &
     299            6 :             ana_env%density_3d%nr_bins, &
     300            6 :             ana_env%density_3d%sum_vol, &
     301            6 :             ana_env%density_3d%sum_vol2, &
     302            6 :             ana_env%density_3d%sum_box_length, &
     303            6 :             ana_env%density_3d%sum_box_length2, &
     304           30 :             ana_env%density_3d%sum_density, &
     305           36 :             ana_env%density_3d%sum_dens2
     306              :       END IF
     307              : 
     308            6 :       l_tmp = ASSOCIATED(ana_env%pair_correl)
     309            6 :       WRITE (file_ptr) l_tmp
     310            6 :       IF (l_tmp) THEN
     311            6 :          WRITE (file_ptr) ana_env%pair_correl%conf_counter, &
     312            6 :             ana_env%pair_correl%nr_bins, &
     313            6 :             ana_env%pair_correl%step_length, &
     314           24 :             ana_env%pair_correl%pairs, &
     315         5436 :             ana_env%pair_correl%g_r
     316              :       END IF
     317              : 
     318            6 :       l_tmp = ASSOCIATED(ana_env%dip_mom)
     319            6 :       WRITE (file_ptr) l_tmp
     320            6 :       IF (l_tmp) THEN
     321            6 :          WRITE (file_ptr) ana_env%dip_mom%conf_counter, &
     322          132 :             ana_env%dip_mom%charges, &
     323           30 :             ana_env%dip_mom%last_dip_cl
     324              :       END IF
     325              : 
     326            6 :       l_tmp = ASSOCIATED(ana_env%dip_ana)
     327            6 :       WRITE (file_ptr) l_tmp
     328            6 :       IF (l_tmp) THEN
     329            0 :          WRITE (file_ptr) ana_env%dip_ana%conf_counter, &
     330            0 :             ana_env%dip_ana%ana_type, &
     331            0 :             ana_env%dip_ana%mu2_pv_s, &
     332            0 :             ana_env%dip_ana%mu_psv, &
     333            0 :             ana_env%dip_ana%mu_pv, &
     334            0 :             ana_env%dip_ana%mu2_pv_mat, &
     335            0 :             ana_env%dip_ana%mu2_pv_mat
     336              :       END IF
     337              : 
     338            6 :       l_tmp = ASSOCIATED(ana_env%displace)
     339            6 :       WRITE (file_ptr) l_tmp
     340            6 :       IF (l_tmp) THEN
     341            6 :          WRITE (file_ptr) ana_env%displace%conf_counter, &
     342           12 :             ana_env%displace%disp
     343              :       END IF
     344              : 
     345            6 :       CALL close_file(unit_number=file_ptr)
     346              : 
     347              :       file_name_tmp = expand_file_name_char(TRIM(ana_env%out_file_prefix)// &
     348            6 :                                             tmc_default_restart_in_file_name, "ana")
     349              :       file_name = expand_file_name_temp(file_name_tmp, &
     350            6 :                                         ana_env%temperature)
     351              :       CALL open_file(file_name=file_name, &
     352              :                      file_action="WRITE", file_status="REPLACE", &
     353            6 :                      unit_number=file_ptr)
     354            6 :       WRITE (file_ptr, *) TRIM(restart_file_name)
     355            6 :       CALL close_file(unit_number=file_ptr)
     356              :    END SUBROUTINE analysis_restart_print
     357              : 
     358              : ! **************************************************************************************************
     359              : !> \brief read analysis restart file
     360              : !> \param ana_env ...
     361              : !> \param elem ...
     362              : !> \param
     363              : !> \author Mandes 02.2013
     364              : ! **************************************************************************************************
     365           18 :    SUBROUTINE analysis_restart_read(ana_env, elem)
     366              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     367              :       TYPE(tree_type), POINTER                           :: elem
     368              : 
     369              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_tmp
     370              :       INTEGER                                            :: file_ptr
     371              :       LOGICAL                                            :: l_tmp
     372              :       REAL(KIND=dp)                                      :: temp
     373              : 
     374           18 :       CPASSERT(ASSOCIATED(ana_env))
     375           18 :       CPASSERT(ASSOCIATED(elem))
     376           18 :       IF (.NOT. ana_env%restart) RETURN
     377              : 
     378              :       file_name_tmp = expand_file_name_char(TRIM(ana_env%out_file_prefix)// &
     379            6 :                                             tmc_default_restart_in_file_name, "ana")
     380              :       file_name = expand_file_name_temp(file_name_tmp, &
     381            6 :                                         ana_env%temperature)
     382            6 :       INQUIRE (FILE=file_name, EXIST=l_tmp)
     383            6 :       IF (l_tmp) THEN
     384              :          CALL open_file(file_name=file_name, file_status="OLD", &
     385            3 :                         file_action="READ", unit_number=file_ptr)
     386            3 :          READ (file_ptr, *) file_name_tmp
     387            3 :          CALL close_file(unit_number=file_ptr)
     388              : 
     389              :          CALL open_file(file_name=file_name_tmp, file_status="OLD", file_form="UNFORMATTED", &
     390            3 :                         file_action="READ", unit_number=file_ptr)
     391            3 :          READ (file_ptr) temp
     392            3 :          CPASSERT(ana_env%temperature .EQ. temp)
     393            3 :          ana_env%last_elem => elem
     394            3 :          CALL read_subtree_elem_unformated(elem, file_ptr)
     395              : 
     396              :          ! first mention the different kind of anlysis types initialized
     397              :          ! then the variables for each calculation type
     398            3 :          READ (file_ptr) l_tmp
     399            3 :          CPASSERT(ASSOCIATED(ana_env%density_3d) .EQV. l_tmp)
     400            3 :          IF (l_tmp) THEN
     401            3 :             READ (file_ptr) ana_env%density_3d%conf_counter, &
     402            3 :                ana_env%density_3d%nr_bins, &
     403            3 :                ana_env%density_3d%sum_vol, &
     404            3 :                ana_env%density_3d%sum_vol2, &
     405            3 :                ana_env%density_3d%sum_box_length, &
     406            3 :                ana_env%density_3d%sum_box_length2, &
     407           15 :                ana_env%density_3d%sum_density, &
     408           18 :                ana_env%density_3d%sum_dens2
     409              :          END IF
     410              : 
     411            3 :          READ (file_ptr) l_tmp
     412            3 :          CPASSERT(ASSOCIATED(ana_env%pair_correl) .EQV. l_tmp)
     413            3 :          IF (l_tmp) THEN
     414            3 :             READ (file_ptr) ana_env%pair_correl%conf_counter, &
     415            3 :                ana_env%pair_correl%nr_bins, &
     416            3 :                ana_env%pair_correl%step_length, &
     417           12 :                ana_env%pair_correl%pairs, &
     418         2718 :                ana_env%pair_correl%g_r
     419              :          END IF
     420              : 
     421            3 :          READ (file_ptr) l_tmp
     422            3 :          CPASSERT(ASSOCIATED(ana_env%dip_mom) .EQV. l_tmp)
     423            3 :          IF (l_tmp) THEN
     424            3 :             READ (file_ptr) ana_env%dip_mom%conf_counter, &
     425           66 :                ana_env%dip_mom%charges, &
     426           15 :                ana_env%dip_mom%last_dip_cl
     427              :          END IF
     428              : 
     429            3 :          READ (file_ptr) l_tmp
     430            3 :          CPASSERT(ASSOCIATED(ana_env%dip_ana) .EQV. l_tmp)
     431            3 :          IF (l_tmp) THEN
     432            0 :             READ (file_ptr) ana_env%dip_ana%conf_counter, &
     433            0 :                ana_env%dip_ana%ana_type, &
     434            0 :                ana_env%dip_ana%mu2_pv_s, &
     435            0 :                ana_env%dip_ana%mu_psv, &
     436            0 :                ana_env%dip_ana%mu_pv, &
     437            0 :                ana_env%dip_ana%mu2_pv_mat, &
     438            0 :                ana_env%dip_ana%mu2_pv_mat
     439              :          END IF
     440              : 
     441            3 :          READ (file_ptr) l_tmp
     442            3 :          CPASSERT(ASSOCIATED(ana_env%displace) .EQV. l_tmp)
     443            3 :          IF (l_tmp) THEN
     444            3 :             READ (file_ptr) ana_env%displace%conf_counter, &
     445            6 :                ana_env%displace%disp
     446              :          END IF
     447              : 
     448            3 :          CALL close_file(unit_number=file_ptr)
     449            3 :          elem => NULL()
     450              :       END IF
     451              :    END SUBROUTINE analysis_restart_read
     452              : 
     453              : ! **************************************************************************************************
     454              : !> \brief call all the necessarry analysis routines
     455              : !>         analysis the previous element with the weight of the different
     456              : !>        configuration numbers
     457              : !>        and stores the actual in the structur % last_elem
     458              : !>        afterwards the previous configuration can be deallocated (outside)
     459              : !> \param elem ...
     460              : !> \param ana_env ...
     461              : !> \param
     462              : !> \author Mandes 02.2013
     463              : ! **************************************************************************************************
     464         2062 :    SUBROUTINE do_tmc_analysis(elem, ana_env)
     465              :       TYPE(tree_type), POINTER                           :: elem
     466              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     467              : 
     468              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'do_tmc_analysis'
     469              : 
     470              :       INTEGER                                            :: handle, weight_act
     471              :       REAL(KIND=dp), DIMENSION(3)                        :: dip_tmp
     472              :       TYPE(tree_type), POINTER                           :: elem_tmp
     473              : 
     474         1031 :       CPASSERT(ASSOCIATED(elem))
     475         1031 :       CPASSERT(ASSOCIATED(ana_env))
     476              : 
     477              :       ! start the timing
     478         1031 :       CALL timeset(routineN, handle)
     479              : 
     480         1031 :       weight_act = 0
     481         1031 :       IF (ASSOCIATED(ana_env%last_elem)) THEN
     482         1016 :          weight_act = elem%nr - ana_env%last_elem%nr
     483              :       END IF
     484              : 
     485         1031 :       IF (weight_act .GT. 0) THEN
     486              :          ! calculates the 3 dimensional distributed density
     487         1016 :          IF (ASSOCIATED(ana_env%density_3d)) &
     488              :             CALL calc_density_3d(elem=ana_env%last_elem, &
     489              :                                  weight=weight_act, atoms=ana_env%atoms, &
     490          500 :                                  ana_env=ana_env)
     491              :          ! calculated the radial distribution function for each atom type
     492         1016 :          IF (ASSOCIATED(ana_env%pair_correl)) &
     493              :             CALL calc_paircorrelation(elem=ana_env%last_elem, weight=weight_act, &
     494          500 :                                       atoms=ana_env%atoms, ana_env=ana_env)
     495              :          ! calculates the classical dipole moments
     496         1016 :          IF (ASSOCIATED(ana_env%dip_mom)) &
     497              :             CALL calc_dipole_moment(elem=ana_env%last_elem, weight=weight_act, &
     498          500 :                                     ana_env=ana_env)
     499              :          ! calculates the dipole moments analysis and dielectric constant
     500         1016 :          IF (ASSOCIATED(ana_env%dip_ana)) THEN
     501              :             ! in symmetric case use also the dipoles
     502              :             !   (-x,y,z) .. .. (-x,-y,z).... (-x,-y-z) all have the same energy
     503            0 :             IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) THEN
     504              :                ! (-x,y,z)
     505            0 :                ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
     506            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     507            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     508            0 :                   ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
     509              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     510            0 :                                          ana_env=ana_env)
     511              :                ! (-x,-y,z)
     512            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     513            0 :                ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
     514            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     515            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     516            0 :                   ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
     517              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     518            0 :                                          ana_env=ana_env)
     519              :                ! (-x,-y,-z)
     520            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     521            0 :                ana_env%last_elem%dipole(3) = -ana_env%last_elem%dipole(3)
     522            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     523            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     524            0 :                   ana_env%dip_mom%last_dip_cl(3) = -ana_env%dip_mom%last_dip_cl(3)
     525              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     526            0 :                                          ana_env=ana_env)
     527              :                ! (x,-y,-z)
     528            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     529            0 :                ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
     530            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     531            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     532            0 :                   ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
     533              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     534            0 :                                          ana_env=ana_env)
     535              :                ! (x,y,-z)
     536            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     537            0 :                ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
     538            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     539            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     540            0 :                   ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
     541              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     542            0 :                                          ana_env=ana_env)
     543              :                ! (-x,y,-z)
     544            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     545            0 :                ana_env%last_elem%dipole(1) = -ana_env%last_elem%dipole(1)
     546            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     547            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     548            0 :                   ana_env%dip_mom%last_dip_cl(1) = -ana_env%dip_mom%last_dip_cl(1)
     549              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     550            0 :                                          ana_env=ana_env)
     551              :                ! (x,-y,z)
     552            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     553            0 :                ana_env%last_elem%dipole(:) = -ana_env%last_elem%dipole(:)
     554            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     555            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     556            0 :                   ana_env%dip_mom%last_dip_cl(:) = -ana_env%dip_mom%last_dip_cl(:)
     557              :                CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     558            0 :                                          ana_env=ana_env)
     559              :                ! back to (x,y,z)
     560            0 :                ana_env%last_elem%dipole(:) = dip_tmp(:)
     561            0 :                ana_env%last_elem%dipole(2) = -ana_env%last_elem%dipole(2)
     562            0 :                dip_tmp(:) = ana_env%last_elem%dipole(:)
     563            0 :                IF (ASSOCIATED(ana_env%dip_mom)) &
     564            0 :                   ana_env%dip_mom%last_dip_cl(2) = -ana_env%dip_mom%last_dip_cl(2)
     565              :             END IF
     566              :             CALL calc_dipole_analysis(elem=ana_env%last_elem, weight=weight_act, &
     567            0 :                                       ana_env=ana_env)
     568              :             CALL print_act_dipole_analysis(elem=ana_env%last_elem, &
     569            0 :                                            ana_env=ana_env)
     570              :          END IF
     571              : 
     572              :          ! calculates the cell displacement from last cell
     573         1016 :          IF (ASSOCIATED(ana_env%displace)) THEN
     574          500 :             CALL calc_displacement(elem=elem, ana_env=ana_env)
     575              :          END IF
     576              :       END IF
     577              :       ! swap elem with last elem, to delete original last element and store the actual one
     578         1031 :       elem_tmp => ana_env%last_elem
     579         1031 :       ana_env%last_elem => elem
     580         1031 :       elem => elem_tmp
     581              :       ! end the timing
     582         1031 :       CALL timestop(handle)
     583         1031 :    END SUBROUTINE do_tmc_analysis
     584              : 
     585              : ! **************************************************************************************************
     586              : !> \brief call all the necessarry analysis printing routines
     587              : !> \param ana_env ...
     588              : !> \param
     589              : !> \author Mandes 02.2013
     590              : ! **************************************************************************************************
     591           36 :    SUBROUTINE finalize_tmc_analysis(ana_env)
     592              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     593              : 
     594              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'finalize_tmc_analysis'
     595              : 
     596              :       INTEGER                                            :: handle
     597              : 
     598           18 :       CPASSERT(ASSOCIATED(ana_env))
     599              : 
     600              :       ! start the timing
     601           18 :       CALL timeset(routineN, handle)
     602           18 :       IF (ASSOCIATED(ana_env%density_3d)) THEN
     603            9 :          IF (ana_env%density_3d%conf_counter .GT. 0) &
     604            9 :             CALL print_density_3d(ana_env=ana_env)
     605              :       END IF
     606           18 :       IF (ASSOCIATED(ana_env%pair_correl)) THEN
     607            9 :          IF (ana_env%pair_correl%conf_counter .GT. 0) &
     608            9 :             CALL print_paircorrelation(ana_env=ana_env)
     609              :       END IF
     610           18 :       IF (ASSOCIATED(ana_env%dip_mom)) THEN
     611            9 :          IF (ana_env%dip_mom%conf_counter .GT. 0) &
     612            9 :             CALL print_dipole_moment(ana_env)
     613              :       END IF
     614           18 :       IF (ASSOCIATED(ana_env%dip_ana)) THEN
     615            0 :          IF (ana_env%dip_ana%conf_counter .GT. 0) &
     616            0 :             CALL print_dipole_analysis(ana_env)
     617              :       END IF
     618           18 :       IF (ASSOCIATED(ana_env%displace)) THEN
     619            9 :          IF (ana_env%displace%conf_counter .GT. 0) &
     620            9 :             CALL print_average_displacement(ana_env)
     621              :       END IF
     622              : 
     623              :       ! end the timing
     624           18 :       CALL timestop(handle)
     625           18 :    END SUBROUTINE finalize_tmc_analysis
     626              : 
     627              : ! **************************************************************************************************
     628              : !> \brief read the files and analyze the configurations
     629              : !> \param start_id ...
     630              : !> \param end_id ...
     631              : !> \param dir_ind ...
     632              : !> \param ana_env ...
     633              : !> \param tmc_params ...
     634              : !> \author Mandes 03.2013
     635              : ! **************************************************************************************************
     636           36 :    SUBROUTINE analyze_file_configurations(start_id, end_id, dir_ind, &
     637              :                                           ana_env, tmc_params)
     638              :       INTEGER                                            :: start_id, end_id
     639              :       INTEGER, OPTIONAL                                  :: dir_ind
     640              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     641              :       TYPE(tmc_param_type), POINTER                      :: tmc_params
     642              : 
     643              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'analyze_file_configurations'
     644              : 
     645              :       INTEGER                                            :: conf_nr, handle, nr_dim, stat
     646              :       TYPE(tree_type), POINTER                           :: elem
     647              : 
     648           18 :       NULLIFY (elem)
     649           18 :       conf_nr = -1
     650           18 :       stat = TMC_STATUS_WAIT_FOR_NEW_TASK
     651           18 :       CPASSERT(ASSOCIATED(ana_env))
     652           18 :       CPASSERT(ASSOCIATED(tmc_params))
     653              : 
     654              :       ! start the timing
     655           18 :       CALL timeset(routineN, handle)
     656              : 
     657              :       ! open the files
     658           18 :       CALL analyse_files_open(tmc_ana=ana_env, stat=stat, dir_ind=dir_ind)
     659              :       ! set the existence of exact dipoles (from file)
     660           18 :       IF (ana_env%id_dip .GT. 0) THEN
     661            0 :          tmc_params%print_dipole = .TRUE.
     662              :       ELSE
     663           18 :          tmc_params%print_dipole = .FALSE.
     664              :       END IF
     665              : 
     666              :       ! allocate the actual element structure
     667              :       CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
     668           18 :                                       nr_dim=ana_env%nr_dim)
     669              : 
     670           18 :       IF (ASSOCIATED(ana_env%last_elem)) conf_nr = ana_env%last_elem%nr
     671           18 :       nr_dim = SIZE(elem%pos)
     672              : 
     673           18 :       IF (stat .EQ. TMC_STATUS_OK) THEN
     674              :          conf_loop: DO
     675              :             CALL read_element_from_file(elem=elem, tmc_ana=ana_env, conf_nr=conf_nr, &
     676         1049 :                                         stat=stat)
     677         1049 :             IF (stat .EQ. TMC_STATUS_WAIT_FOR_NEW_TASK) THEN
     678           18 :                CALL deallocate_sub_tree_node(tree_elem=elem)
     679              :                EXIT conf_loop
     680              :             END IF
     681              :             ! if we want just a certain part of the trajectory
     682         1031 :             IF (start_id .LT. 0 .OR. conf_nr .GE. start_id) THEN
     683         1031 :                IF (end_id .LT. 0 .OR. conf_nr .LE. end_id) THEN
     684              :                   ! do the analysis calculations
     685         1031 :                   CALL do_tmc_analysis(elem=elem, ana_env=ana_env)
     686              :                END IF
     687              :             END IF
     688              : 
     689              :             ! clean temporary element (already analyzed)
     690         1031 :             IF (ASSOCIATED(elem)) THEN
     691         1016 :                CALL deallocate_sub_tree_node(tree_elem=elem)
     692              :             END IF
     693              :             ! if there was no previous element, create a new temp element to write in
     694         1031 :             IF (.NOT. ASSOCIATED(elem)) &
     695              :                CALL allocate_new_sub_tree_node(tmc_params=tmc_params, next_el=elem, &
     696         1031 :                                                nr_dim=nr_dim)
     697              :          END DO conf_loop
     698              :       END IF
     699              :       ! close the files
     700           18 :       CALL analyse_files_close(tmc_ana=ana_env)
     701              : 
     702           18 :       IF (ASSOCIATED(elem)) &
     703            0 :          CALL deallocate_sub_tree_node(tree_elem=elem)
     704              : 
     705              :       ! end the timing
     706           18 :       CALL timestop(handle)
     707           18 :    END SUBROUTINE analyze_file_configurations
     708              : 
     709              :    !============================================================================
     710              :    ! density calculations
     711              :    !============================================================================
     712              : 
     713              : ! **************************************************************************************************
     714              : !> \brief calculates the density in rectantangulares
     715              : !>        defined by the number of bins in each direction
     716              : !> \param elem ...
     717              : !> \param weight ...
     718              : !> \param atoms ...
     719              : !> \param ana_env ...
     720              : !> \param
     721              : !> \author Mandes 02.2013
     722              : ! **************************************************************************************************
     723          500 :    SUBROUTINE calc_density_3d(elem, weight, atoms, ana_env)
     724              :       TYPE(tree_type), POINTER                           :: elem
     725              :       INTEGER                                            :: weight
     726              :       TYPE(tmc_atom_type), DIMENSION(:), POINTER         :: atoms
     727              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     728              : 
     729              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_density_3d'
     730              : 
     731              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_tmp
     732              :       INTEGER                                            :: atom, bin_x, bin_y, bin_z, file_ptr, &
     733              :                                                             handle
     734              :       LOGICAL                                            :: flag
     735              :       REAL(KIND=dp)                                      :: mass_total, r_tmp, vol_cell, vol_sub_box
     736              :       REAL(KIND=dp), DIMENSION(3)                        :: atom_pos, cell_size, interval_size
     737          500 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: mass_bin
     738              : 
     739          500 :       NULLIFY (mass_bin)
     740              : 
     741            0 :       CPASSERT(ASSOCIATED(elem))
     742          500 :       CPASSERT(ASSOCIATED(elem%pos))
     743          500 :       CPASSERT(weight .GT. 0)
     744          500 :       CPASSERT(ASSOCIATED(atoms))
     745          500 :       CPASSERT(ASSOCIATED(ana_env))
     746          500 :       CPASSERT(ASSOCIATED(ana_env%cell))
     747          500 :       CPASSERT(ASSOCIATED(ana_env%density_3d))
     748          500 :       CPASSERT(ASSOCIATED(ana_env%density_3d%sum_density))
     749          500 :       CPASSERT(ASSOCIATED(ana_env%density_3d%sum_dens2))
     750              : 
     751              :       ! start the timing
     752          500 :       CALL timeset(routineN, handle)
     753              : 
     754          500 :       atom_pos(:) = 0.0_dp
     755          500 :       cell_size(:) = 0.0_dp
     756          500 :       interval_size(:) = 0.0_dp
     757          500 :       mass_total = 0.0_dp
     758              : 
     759          500 :       bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
     760          500 :       bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
     761          500 :       bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
     762         2500 :       ALLOCATE (mass_bin(bin_x, bin_y, bin_z))
     763         2500 :       mass_bin(:, :, :) = 0.0_dp
     764              : 
     765              :       ! if NPT -> box_scale.NE.1.0 use the scaled cell
     766              :       ! ATTENTION then the sub box middle points are not correct in the output
     767              :       !  espacially if we use multiple sub boxes
     768              :       CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
     769          500 :                            abc=cell_size, vol=vol_cell)
     770              :       ! volume summed over configurations for average volume [A]
     771              :       ana_env%density_3d%sum_vol = ana_env%density_3d%sum_vol + &
     772          500 :                                    vol_cell*(au2a)**3*weight
     773              :       ana_env%density_3d%sum_vol2 = ana_env%density_3d%sum_vol2 + &
     774          500 :                                     (vol_cell*(au2a)**3)**2*weight
     775              : 
     776              :       ana_env%density_3d%sum_box_length(:) = ana_env%density_3d%sum_box_length(:) &
     777         2000 :                                              + cell_size(:)*(au2a)*weight
     778              :       ana_env%density_3d%sum_box_length2(:) = ana_env%density_3d%sum_box_length2(:) &
     779         2000 :                                               + (cell_size(:)*(au2a))**2*weight
     780              : 
     781              :       ! sub interval length
     782          500 :       interval_size(1) = cell_size(1)/REAL(bin_x, dp)
     783          500 :       interval_size(2) = cell_size(2)/REAL(bin_y, dp)
     784          500 :       interval_size(3) = cell_size(3)/REAL(bin_z, dp)
     785              : 
     786              :       ! volume in [cm^3]
     787          500 :       vol_cell = vol_cell*(au2a*1E-8)**3
     788              :       vol_sub_box = interval_size(1)*interval_size(2)*interval_size(3)* &
     789          500 :                     (au2a*1E-8)**3
     790              : 
     791              :       ! count every atom
     792          500 :       DO atom = 1, SIZE(elem%pos), ana_env%dim_per_elem
     793              : 
     794        42000 :          atom_pos(:) = elem%pos(atom:atom + 2)
     795              :          ! fold into box
     796              :          CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
     797        10500 :                               vec=atom_pos)
     798              :          ! shifts the box to positive values (before 0,0,0 is the center)
     799        42000 :          atom_pos(:) = atom_pos(:) + 0.5_dp*cell_size(:)
     800              :          ! calculate the index of the sub box
     801        10500 :          bin_x = INT(atom_pos(1)/interval_size(1)) + 1
     802        10500 :          bin_y = INT(atom_pos(2)/interval_size(2)) + 1
     803        10500 :          bin_z = INT(atom_pos(3)/interval_size(3)) + 1
     804        10500 :          CPASSERT(bin_x .GT. 0 .AND. bin_y .GT. 0 .AND. bin_z .GT. 0)
     805        10500 :          CPASSERT(bin_x .LE. SIZE(ana_env%density_3d%sum_density(:, 1, 1)))
     806        10500 :          CPASSERT(bin_y .LE. SIZE(ana_env%density_3d%sum_density(1, :, 1)))
     807        10500 :          CPASSERT(bin_z .LE. SIZE(ana_env%density_3d%sum_density(1, 1, :)))
     808              : 
     809              :          ! sum mass in [g] (in bins and total)
     810              :          mass_bin(bin_x, bin_y, bin_z) = mass_bin(bin_x, bin_y, bin_z) + &
     811        10500 :                                          atoms(INT(atom/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%mass/massunit*1000*a_mass
     812              :          mass_total = mass_total + &
     813        10500 :                       atoms(INT(atom/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%mass/massunit*1000*a_mass
     814              :          !mass_bin(bin_x,bin_y,bin_z) = mass_bin(bin_x,bin_y,bin_z) + &
     815              :          !  atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
     816              :          !     massunit/n_avogadro
     817              :          !mass_total = mass_total + &
     818              :          !  atoms(INT(atom/REAL(ana_env%dim_per_elem,KIND=dp))+1)%mass/&
     819              :          !     massunit/n_avogadro
     820              :       END DO
     821              :       ! check total cell density
     822         4000 :       r_tmp = mass_total/vol_cell - SUM(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :))
     823          500 :       CPASSERT(ABS(r_tmp) .LT. 1E-5)
     824              : 
     825              :       ! calculate density (mass per volume) and sum up for average value
     826              :       ana_env%density_3d%sum_density(:, :, :) = ana_env%density_3d%sum_density(:, :, :) + &
     827         4500 :                                                 weight*mass_bin(:, :, :)/vol_sub_box
     828              : 
     829              :       ! calculate density squared ( (mass per volume)^2 ) for variance and sum up for average value
     830              :       ana_env%density_3d%sum_dens2(:, :, :) = ana_env%density_3d%sum_dens2(:, :, :) + &
     831         4500 :                                               weight*(mass_bin(:, :, :)/vol_sub_box)**2
     832              : 
     833          500 :       ana_env%density_3d%conf_counter = ana_env%density_3d%conf_counter + weight
     834              : 
     835              :       ! print out the actual and average density in file
     836          500 :       IF (ana_env%density_3d%print_dens) THEN
     837              :          file_name_tmp = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
     838              :                                                tmc_default_trajectory_file_name, &
     839          500 :                                                ana_env%temperature)
     840              :          file_name = TRIM(expand_file_name_char(file_name_tmp, &
     841          500 :                                                 "dens"))
     842          500 :          INQUIRE (FILE=file_name, EXIST=flag)
     843              :          CALL open_file(file_name=file_name, file_status="UNKNOWN", &
     844              :                         file_action="WRITE", file_position="APPEND", &
     845          500 :                         unit_number=file_ptr)
     846          500 :          IF (.NOT. flag) &
     847            3 :             WRITE (file_ptr, FMT='(A8,11A20)') "# conf_nr", "dens_act[g/cm^3]", &
     848            3 :             "dens_average[g/cm^3]", "density_variance", &
     849            3 :             "averages:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z", &
     850            6 :             "variances:volume", "box_lenth_x", "box_lenth_y", "box_lenth_z"
     851          500 :          WRITE (file_ptr, FMT="(I8,11F20.10)") ana_env%density_3d%conf_counter + 1 - weight, &
     852         4000 :             SUM(mass_bin(:, :, :))/vol_sub_box/SIZE(mass_bin(:, :, :)), &
     853              :             SUM(ana_env%density_3d%sum_density(:, :, :))/ &
     854              :             SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
     855         4000 :             REAL(ana_env%density_3d%conf_counter, KIND=dp), &
     856              :             SUM(ana_env%density_3d%sum_dens2(:, :, :))/ &
     857              :             SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
     858              :             REAL(ana_env%density_3d%conf_counter, KIND=dp) - &
     859              :             (SUM(ana_env%density_3d%sum_density(:, :, :))/ &
     860              :              SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
     861         7500 :              REAL(ana_env%density_3d%conf_counter, KIND=dp))**2, &
     862              :             ana_env%density_3d%sum_vol/ &
     863          500 :             REAL(ana_env%density_3d%conf_counter, KIND=dp), &
     864              :             ana_env%density_3d%sum_box_length(:)/ &
     865         2000 :             REAL(ana_env%density_3d%conf_counter, KIND=dp), &
     866              :             ana_env%density_3d%sum_vol2/ &
     867              :             REAL(ana_env%density_3d%conf_counter, KIND=dp) - &
     868              :             (ana_env%density_3d%sum_vol/ &
     869          500 :              REAL(ana_env%density_3d%conf_counter, KIND=dp))**2, &
     870              :             ana_env%density_3d%sum_box_length2(:)/ &
     871              :             REAL(ana_env%density_3d%conf_counter, KIND=dp) - &
     872              :             (ana_env%density_3d%sum_box_length(:)/ &
     873         2500 :              REAL(ana_env%density_3d%conf_counter, KIND=dp))**2
     874          500 :          CALL close_file(unit_number=file_ptr)
     875              :       END IF
     876              : 
     877          500 :       DEALLOCATE (mass_bin)
     878              :       ! end the timing
     879          500 :       CALL timestop(handle)
     880         1000 :    END SUBROUTINE calc_density_3d
     881              : 
     882              : ! **************************************************************************************************
     883              : !> \brief print the density in rectantangulares
     884              : !>        defined by the number of bins in each direction
     885              : !> \param ana_env ...
     886              : !> \param
     887              : !> \author Mandes 02.2013
     888              : ! **************************************************************************************************
     889           18 :    SUBROUTINE print_density_3d(ana_env)
     890              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
     891              : 
     892              :       CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA", &
     893              :          routineN = 'print_density_3d'
     894              : 
     895              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_vari
     896              :       INTEGER                                            :: bin_x, bin_y, bin_z, file_ptr_dens, &
     897              :                                                             file_ptr_vari, handle, i, j, k
     898              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_size, interval_size
     899              : 
     900            9 :       CPASSERT(ASSOCIATED(ana_env))
     901            9 :       CPASSERT(ASSOCIATED(ana_env%density_3d))
     902            9 :       CPASSERT(ASSOCIATED(ana_env%density_3d%sum_density))
     903            9 :       CPASSERT(ASSOCIATED(ana_env%density_3d%sum_dens2))
     904              : 
     905              :       ! start the timing
     906            9 :       CALL timeset(routineN, handle)
     907              : 
     908              :       file_name = ""
     909            9 :       file_name_vari = ""
     910              : 
     911            9 :       bin_x = SIZE(ana_env%density_3d%sum_density(:, 1, 1))
     912            9 :       bin_y = SIZE(ana_env%density_3d%sum_density(1, :, 1))
     913            9 :       bin_z = SIZE(ana_env%density_3d%sum_density(1, 1, :))
     914            9 :       CALL get_cell(cell=ana_env%cell, abc=cell_size)
     915            9 :       interval_size(1) = cell_size(1)/REAL(bin_x, KIND=dp)*au2a
     916            9 :       interval_size(2) = cell_size(2)/REAL(bin_y, KIND=dp)*au2a
     917            9 :       interval_size(3) = cell_size(3)/REAL(bin_z, KIND=dp)*au2a
     918              : 
     919              :       file_name = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
     920              :                                         tmc_ana_density_file_name, &
     921            9 :                                         ana_env%temperature)
     922              :       CALL open_file(file_name=file_name, file_status="REPLACE", &
     923              :                      file_action="WRITE", file_position="APPEND", &
     924            9 :                      unit_number=file_ptr_dens)
     925              :       WRITE (file_ptr_dens, FMT='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
     926            9 :          "# configurations", ana_env%density_3d%conf_counter, "bins", &
     927           18 :          ana_env%density_3d%nr_bins, "interval size", interval_size(:)
     928            9 :       WRITE (file_ptr_dens, FMT='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " density [g/cm^3] "
     929              : 
     930              :       file_name_vari = expand_file_name_temp(expand_file_name_char( &
     931              :                                              TRIM(ana_env%out_file_prefix)// &
     932              :                                              tmc_ana_density_file_name, "vari"), &
     933            9 :                                              ana_env%temperature)
     934              :       CALL open_file(file_name=file_name_vari, file_status="REPLACE", &
     935              :                      file_action="WRITE", file_position="APPEND", &
     936            9 :                      unit_number=file_ptr_vari)
     937              :       WRITE (file_ptr_vari, FMT='(A,1X,I0,1X,A,3(I0,1X),1X,A,1X,3F10.5)') &
     938            9 :          "# configurations", ana_env%density_3d%conf_counter, "bins", &
     939           18 :          ana_env%density_3d%nr_bins, "interval size", interval_size(:)
     940            9 :       WRITE (file_ptr_vari, FMT='(A,3A10,A20)') "#", " x [A] ", " y [A] ", " z [A] ", " variance"
     941              : 
     942           27 :       DO i = 1, SIZE(ana_env%density_3d%sum_density(:, 1, 1))
     943           45 :          DO j = 1, SIZE(ana_env%density_3d%sum_density(1, :, 1))
     944           54 :             DO k = 1, SIZE(ana_env%density_3d%sum_density(1, 1, :))
     945              :                WRITE (file_ptr_dens, FMT='(3F10.2,F20.10)') &
     946           18 :                   (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
     947           36 :                   ana_env%density_3d%sum_density(i, j, k)/REAL(ana_env%density_3d%conf_counter, KIND=dp)
     948              :                WRITE (file_ptr_vari, FMT='(3F10.2,F20.10)') &
     949           18 :                   (i - 0.5_dp)*interval_size(1), (j - 0.5_dp)*interval_size(2), (k - 0.5_dp)*interval_size(3), &
     950              :                   ana_env%density_3d%sum_dens2(i, j, k)/REAL(ana_env%density_3d%conf_counter, KIND=dp) - &
     951           54 :                   (ana_env%density_3d%sum_density(i, j, k)/REAL(ana_env%density_3d%conf_counter, KIND=dp))**2
     952              :             END DO
     953              :          END DO
     954              :       END DO
     955            9 :       CALL close_file(unit_number=file_ptr_dens)
     956            9 :       CALL close_file(unit_number=file_ptr_vari)
     957              : 
     958            9 :       WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
     959            9 :       WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "density calculation", "-"
     960            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
     961            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations", &
     962           18 :          cp_to_string(REAL(ana_env%density_3d%conf_counter, KIND=dp))
     963            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "average volume", &
     964              :          cp_to_string(ana_env%density_3d%sum_vol/ &
     965           18 :                       REAL(ana_env%density_3d%conf_counter, KIND=dp))
     966            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "average density in the cell: ", &
     967              :          cp_to_string(SUM(ana_env%density_3d%sum_density(:, :, :))/ &
     968              :                       SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
     969           81 :                       REAL(ana_env%density_3d%conf_counter, KIND=dp))
     970            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "density variance:", &
     971              :          cp_to_string(SUM(ana_env%density_3d%sum_dens2(:, :, :))/ &
     972              :                       SIZE(ana_env%density_3d%sum_dens2(:, :, :))/ &
     973              :                       REAL(ana_env%density_3d%conf_counter, KIND=dp) - &
     974              :                       (SUM(ana_env%density_3d%sum_density(:, :, :))/ &
     975              :                        SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
     976          144 :                        REAL(ana_env%density_3d%conf_counter, KIND=dp))**2)
     977            9 :       WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
     978            9 :       IF (ana_env%print_test_output) &
     979            9 :          WRITE (ana_env%io_unit, *) "TMC|ANALYSIS_CELL_DENSITY_X= ", &
     980              :          SUM(ana_env%density_3d%sum_density(:, :, :))/ &
     981              :          SIZE(ana_env%density_3d%sum_density(:, :, :))/ &
     982           81 :          REAL(ana_env%density_3d%conf_counter, KIND=dp)
     983              :       ! end the timing
     984            9 :       CALL timestop(handle)
     985            9 :    END SUBROUTINE print_density_3d
     986              : 
     987              :    !============================================================================
     988              :    ! radial distribution function
     989              :    !============================================================================
     990              : 
     991              : ! **************************************************************************************************
     992              : !> \brief init radial distribution function structures
     993              : !> \param ana_pair_correl ...
     994              : !> \param atoms ...
     995              : !> \param cell ...
     996              : !> \param
     997              : !> \author Mandes 02.2013
     998              : ! **************************************************************************************************
     999            9 :    SUBROUTINE ana_pair_correl_init(ana_pair_correl, atoms, cell)
    1000              :       TYPE(pair_correl_type), POINTER                    :: ana_pair_correl
    1001              :       TYPE(tmc_atom_type), DIMENSION(:), POINTER         :: atoms
    1002              :       TYPE(cell_type), POINTER                           :: cell
    1003              : 
    1004              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ana_pair_correl_init'
    1005              : 
    1006              :       INTEGER                                            :: counter, f_n, handle, list, list_ind, s_n
    1007              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_size
    1008            9 :       TYPE(atom_pairs_type), DIMENSION(:), POINTER       :: pairs_tmp
    1009              : 
    1010            0 :       CPASSERT(ASSOCIATED(ana_pair_correl))
    1011            9 :       CPASSERT(.NOT. ASSOCIATED(ana_pair_correl%g_r))
    1012            9 :       CPASSERT(.NOT. ASSOCIATED(ana_pair_correl%pairs))
    1013            9 :       CPASSERT(ASSOCIATED(atoms))
    1014            9 :       CPASSERT(SIZE(atoms) .GT. 1)
    1015            9 :       CPASSERT(ASSOCIATED(cell))
    1016              : 
    1017              :       ! start the timing
    1018            9 :       CALL timeset(routineN, handle)
    1019              : 
    1020            9 :       CALL get_cell(cell=cell, abc=cell_size)
    1021            9 :       IF (ana_pair_correl%nr_bins .LE. 0) THEN
    1022           45 :          ana_pair_correl%nr_bins = CEILING(MAXVAL(cell_size(:))/2.0_dp/(0.03/au2a))
    1023              :       END IF
    1024              :       ana_pair_correl%step_length = MAXVAL(cell_size(:))/2.0_dp/ &
    1025           45 :                                     ana_pair_correl%nr_bins
    1026            9 :       ana_pair_correl%conf_counter = 0
    1027              : 
    1028            9 :       counter = 1
    1029              :       ! initialise the atom pairs
    1030          216 :       ALLOCATE (pairs_tmp(SIZE(atoms)))
    1031          198 :       DO f_n = 1, SIZE(atoms)
    1032         2088 :          DO s_n = f_n + 1, SIZE(atoms)
    1033              :             ! search if atom pair is already selected
    1034              :             list_ind = search_pair_in_list(pair_list=pairs_tmp, n1=atoms(f_n)%name, &
    1035         1890 :                                            n2=atoms(s_n)%name, list_end=counter - 1)
    1036              :             ! add to list
    1037         2079 :             IF (list_ind .LT. 0) THEN
    1038           27 :                pairs_tmp(counter)%f_n = atoms(f_n)%name
    1039           27 :                pairs_tmp(counter)%s_n = atoms(s_n)%name
    1040           27 :                pairs_tmp(counter)%pair_count = 1
    1041           27 :                counter = counter + 1
    1042              :             ELSE
    1043         1863 :                pairs_tmp(list_ind)%pair_count = pairs_tmp(list_ind)%pair_count + 1
    1044              :             END IF
    1045              :          END DO
    1046              :       END DO
    1047              : 
    1048           54 :       ALLOCATE (ana_pair_correl%pairs(counter - 1))
    1049           36 :       DO list = 1, counter - 1
    1050           27 :          ana_pair_correl%pairs(list)%f_n = pairs_tmp(list)%f_n
    1051           27 :          ana_pair_correl%pairs(list)%s_n = pairs_tmp(list)%s_n
    1052           36 :          ana_pair_correl%pairs(list)%pair_count = pairs_tmp(list)%pair_count
    1053              :       END DO
    1054            9 :       DEALLOCATE (pairs_tmp)
    1055              : 
    1056           36 :       ALLOCATE (ana_pair_correl%g_r(SIZE(ana_pair_correl%pairs(:)), ana_pair_correl%nr_bins))
    1057         8145 :       ana_pair_correl%g_r = 0.0_dp
    1058              :       ! end the timing
    1059            9 :       CALL timestop(handle)
    1060           18 :    END SUBROUTINE ana_pair_correl_init
    1061              : 
    1062              : ! **************************************************************************************************
    1063              : !> \brief calculates the radial distribution function
    1064              : !> \param elem ...
    1065              : !> \param weight ...
    1066              : !> \param atoms ...
    1067              : !> \param ana_env ...
    1068              : !> \param
    1069              : !> \author Mandes 02.2013
    1070              : ! **************************************************************************************************
    1071         1000 :    SUBROUTINE calc_paircorrelation(elem, weight, atoms, ana_env)
    1072              :       TYPE(tree_type), POINTER                           :: elem
    1073              :       INTEGER                                            :: weight
    1074              :       TYPE(tmc_atom_type), DIMENSION(:), POINTER         :: atoms
    1075              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1076              : 
    1077              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_paircorrelation'
    1078              : 
    1079              :       INTEGER                                            :: handle, i, ind, j, pair_ind
    1080              :       REAL(KIND=dp)                                      :: dist
    1081              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_size
    1082              : 
    1083          500 :       CPASSERT(ASSOCIATED(elem))
    1084          500 :       CPASSERT(ASSOCIATED(elem%pos))
    1085         2000 :       CPASSERT(ALL(elem%box_scale(:) .GT. 0.0_dp))
    1086          500 :       CPASSERT(weight .GT. 0)
    1087          500 :       CPASSERT(ASSOCIATED(atoms))
    1088          500 :       CPASSERT(ASSOCIATED(ana_env))
    1089          500 :       CPASSERT(ASSOCIATED(ana_env%cell))
    1090          500 :       CPASSERT(ASSOCIATED(ana_env%pair_correl))
    1091          500 :       CPASSERT(ASSOCIATED(ana_env%pair_correl%g_r))
    1092          500 :       CPASSERT(ASSOCIATED(ana_env%pair_correl%pairs))
    1093              : 
    1094              :       ! start the timing
    1095          500 :       CALL timeset(routineN, handle)
    1096              : 
    1097          500 :       dist = -1.0_dp
    1098              : 
    1099        11000 :       first_elem_loop: DO i = 1, SIZE(elem%pos), ana_env%dim_per_elem
    1100       116000 :          second_elem_loop: DO j = i + 3, SIZE(elem%pos), ana_env%dim_per_elem
    1101              :             dist = nearest_distance(x1=elem%pos(i:i + ana_env%dim_per_elem - 1), &
    1102              :                                     x2=elem%pos(j:j + ana_env%dim_per_elem - 1), &
    1103       105000 :                                     cell=ana_env%cell, box_scale=elem%box_scale)
    1104       105000 :             ind = CEILING(dist/ana_env%pair_correl%step_length)
    1105       115500 :             IF (ind .LE. ana_env%pair_correl%nr_bins) THEN
    1106              :                pair_ind = search_pair_in_list(pair_list=ana_env%pair_correl%pairs, &
    1107              :                                               n1=atoms(INT(i/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%name, &
    1108        86745 :                                               n2=atoms(INT(j/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)%name)
    1109        86745 :                CPASSERT(pair_ind .GT. 0)
    1110              :                ana_env%pair_correl%g_r(pair_ind, ind) = &
    1111        86745 :                   ana_env%pair_correl%g_r(pair_ind, ind) + weight
    1112              :             END IF
    1113              :          END DO second_elem_loop
    1114              :       END DO first_elem_loop
    1115          500 :       ana_env%pair_correl%conf_counter = ana_env%pair_correl%conf_counter + weight
    1116          500 :       CALL get_cell(cell=ana_env%cell, abc=cell_size)
    1117              :       ana_env%pair_correl%sum_box_scale = ana_env%pair_correl%sum_box_scale + &
    1118         4000 :                                           (elem%box_scale(:)*weight)
    1119              :       ! end the timing
    1120          500 :       CALL timestop(handle)
    1121          500 :    END SUBROUTINE calc_paircorrelation
    1122              : 
    1123              : ! **************************************************************************************************
    1124              : !> \brief print the radial distribution function for each pair of atoms
    1125              : !> \param ana_env ...
    1126              : !> \param
    1127              : !> \author Mandes 02.2013
    1128              : ! **************************************************************************************************
    1129           18 :    SUBROUTINE print_paircorrelation(ana_env)
    1130              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1131              : 
    1132              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'print_paircorrelation'
    1133              : 
    1134              :       CHARACTER(LEN=default_path_length)                 :: file_name
    1135              :       INTEGER                                            :: bin, file_ptr, handle, pair
    1136              :       REAL(KIND=dp)                                      :: aver_box_scale(3), vol, voldr
    1137              :       REAL(KIND=dp), DIMENSION(3)                        :: cell_size
    1138              : 
    1139            9 :       CPASSERT(ASSOCIATED(ana_env))
    1140            9 :       CPASSERT(ASSOCIATED(ana_env%pair_correl))
    1141              : 
    1142              :       ! start the timing
    1143            9 :       CALL timeset(routineN, handle)
    1144              : 
    1145            9 :       CALL get_cell(cell=ana_env%cell, abc=cell_size)
    1146           36 :       aver_box_scale(:) = ana_env%pair_correl%sum_box_scale(:)/ana_env%pair_correl%conf_counter
    1147              :       vol = (cell_size(1)*aver_box_scale(1))* &
    1148              :             (cell_size(2)*aver_box_scale(2))* &
    1149            9 :             (cell_size(3)*aver_box_scale(3))
    1150              : 
    1151           36 :       DO pair = 1, SIZE(ana_env%pair_correl%pairs)
    1152              :          file_name = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
    1153              :                                            tmc_ana_pair_correl_file_name, &
    1154           27 :                                            ana_env%temperature)
    1155              :          CALL open_file(file_name=expand_file_name_char( &
    1156              :                         expand_file_name_char(file_name, &
    1157              :                                               ana_env%pair_correl%pairs(pair)%f_n), &
    1158              :                         ana_env%pair_correl%pairs(pair)%s_n), &
    1159              :                         file_status="REPLACE", &
    1160              :                         file_action="WRITE", file_position="APPEND", &
    1161           27 :                         unit_number=file_ptr)
    1162              :          WRITE (file_ptr, *) "# radial distribution function of "// &
    1163              :             TRIM(ana_env%pair_correl%pairs(pair)%f_n)//" and "// &
    1164           27 :             TRIM(ana_env%pair_correl%pairs(pair)%s_n)//" of ", &
    1165           54 :             ana_env%pair_correl%conf_counter, " configurations"
    1166           27 :          WRITE (file_ptr, *) "# using a bin size of ", &
    1167           27 :             ana_env%pair_correl%step_length*au2a, &
    1168           54 :             "[A] (for Vol changes: referring to the reference cell)"
    1169         6129 :          DO bin = 1, ana_env%pair_correl%nr_bins
    1170              :             voldr = 4.0/3.0*PI*ana_env%pair_correl%step_length**3* &
    1171         6102 :                     (REAL(bin, KIND=dp)**3 - REAL(bin - 1, KIND=dp)**3)
    1172         6102 :             WRITE (file_ptr, *) (bin - 0.5)*ana_env%pair_correl%step_length*au2a, &
    1173              :                (ana_env%pair_correl%g_r(pair, bin)/ana_env%pair_correl%conf_counter)/ &
    1174        12231 :                (voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
    1175              :          END DO
    1176           27 :          CALL close_file(unit_number=file_ptr)
    1177              : 
    1178           27 :          IF (ana_env%print_test_output) &
    1179              :             WRITE (*, *) "TMC|ANALYSIS_G_R_"// &
    1180              :             TRIM(ana_env%pair_correl%pairs(pair)%f_n)//"_"// &
    1181           27 :             TRIM(ana_env%pair_correl%pairs(pair)%s_n)//"_X= ", &
    1182              :             SUM(ana_env%pair_correl%g_r(pair, :)/ana_env%pair_correl%conf_counter/ &
    1183         6165 :                 voldr*ana_env%pair_correl%pairs(pair)%pair_count/vol)
    1184              :       END DO
    1185              : 
    1186              :       ! end the timing
    1187            9 :       CALL timestop(handle)
    1188            9 :    END SUBROUTINE print_paircorrelation
    1189              : 
    1190              :    !============================================================================
    1191              :    ! classical cell dipole moment
    1192              :    !============================================================================
    1193              : 
    1194              : ! **************************************************************************************************
    1195              : !> \brief init radial distribution function structures
    1196              : !> \param ana_dip_mom ...
    1197              : !> \param atoms ...
    1198              : !> \param
    1199              : !> \author Mandes 02.2013
    1200              : ! **************************************************************************************************
    1201            9 :    SUBROUTINE ana_dipole_moment_init(ana_dip_mom, atoms)
    1202              :       TYPE(dipole_moment_type), POINTER                  :: ana_dip_mom
    1203              :       TYPE(tmc_atom_type), DIMENSION(:), POINTER         :: atoms
    1204              : 
    1205              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'ana_dipole_moment_init'
    1206              : 
    1207              :       INTEGER                                            :: atom, charge, handle
    1208              : 
    1209            9 :       CPASSERT(ASSOCIATED(ana_dip_mom))
    1210            9 :       CPASSERT(ASSOCIATED(ana_dip_mom%charges_inp))
    1211            9 :       CPASSERT(ASSOCIATED(atoms))
    1212              : 
    1213              :       ! start the timing
    1214            9 :       CALL timeset(routineN, handle)
    1215              : 
    1216           27 :       ALLOCATE (ana_dip_mom%charges(SIZE(atoms)))
    1217          198 :       ana_dip_mom%charges = 0.0_dp
    1218              :       ! for every atom searcht the correct charge
    1219          198 :       DO atom = 1, SIZE(atoms)
    1220          324 :          charge_loop: DO charge = 1, SIZE(ana_dip_mom%charges_inp)
    1221          315 :             IF (atoms(atom)%name .EQ. ana_dip_mom%charges_inp(charge)%name) THEN
    1222          189 :                ana_dip_mom%charges(atom) = ana_dip_mom%charges_inp(charge)%mass
    1223          189 :                EXIT charge_loop
    1224              :             END IF
    1225              :          END DO charge_loop
    1226              :       END DO
    1227              : 
    1228            9 :       DEALLOCATE (ana_dip_mom%charges_inp)
    1229              :       ! end the timing
    1230            9 :       CALL timestop(handle)
    1231            9 :    END SUBROUTINE ana_dipole_moment_init
    1232              : 
    1233              : ! **************************************************************************************************
    1234              : !> \brief calculates the classical cell dipole moment
    1235              : !> \param elem ...
    1236              : !> \param weight ...
    1237              : !> \param ana_env ...
    1238              : !> \param
    1239              : !> \author Mandes 02.2013
    1240              : ! **************************************************************************************************
    1241          500 :    SUBROUTINE calc_dipole_moment(elem, weight, ana_env)
    1242              :       TYPE(tree_type), POINTER                           :: elem
    1243              :       INTEGER                                            :: weight
    1244              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1245              : 
    1246              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'calc_dipole_moment'
    1247              : 
    1248              :       CHARACTER(LEN=default_path_length)                 :: file_name
    1249              :       INTEGER                                            :: handle, i
    1250          500 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: dip_cl
    1251              : 
    1252            0 :       CPASSERT(ASSOCIATED(elem))
    1253          500 :       CPASSERT(ASSOCIATED(elem%pos))
    1254          500 :       CPASSERT(ASSOCIATED(ana_env))
    1255          500 :       CPASSERT(ASSOCIATED(ana_env%dip_mom))
    1256          500 :       CPASSERT(ASSOCIATED(ana_env%dip_mom%charges))
    1257              : 
    1258              :       ! start the timing
    1259          500 :       CALL timeset(routineN, handle)
    1260              : 
    1261         1500 :       ALLOCATE (dip_cl(ana_env%dim_per_elem))
    1262         2000 :       dip_cl(:) = 0.0_dp
    1263              : 
    1264          500 :       DO i = 1, SIZE(elem%pos, 1), ana_env%dim_per_elem
    1265              :          dip_cl(:) = dip_cl(:) + elem%pos(i:i + ana_env%dim_per_elem - 1)* &
    1266        84000 :                      ana_env%dip_mom%charges(INT(i/REAL(ana_env%dim_per_elem, KIND=dp)) + 1)
    1267              :       END DO
    1268              : 
    1269              :       ! if there are no exact dipoles save these ones in element structure
    1270          500 :       IF (.NOT. ASSOCIATED(elem%dipole)) THEN
    1271         1500 :          ALLOCATE (elem%dipole(ana_env%dim_per_elem))
    1272         4000 :          elem%dipole(:) = dip_cl(:)
    1273              :       END IF
    1274              : 
    1275          500 :       IF (ana_env%dip_mom%print_cl_dip) THEN
    1276              :          file_name = expand_file_name_temp(tmc_default_trajectory_file_name, &
    1277          500 :                                            ana_env%temperature)
    1278              :          CALL write_dipoles_in_file(file_name=file_name, &
    1279              :                                     conf_nr=ana_env%dip_mom%conf_counter + 1, dip=dip_cl, &
    1280          500 :                                     file_ext="dip_cl")
    1281              :       END IF
    1282          500 :       ana_env%dip_mom%conf_counter = ana_env%dip_mom%conf_counter + weight
    1283         4000 :       ana_env%dip_mom%last_dip_cl(:) = dip_cl
    1284              : 
    1285          500 :       DEALLOCATE (dip_cl)
    1286              : 
    1287              :       ! end the timing
    1288          500 :       CALL timestop(handle)
    1289         1000 :    END SUBROUTINE calc_dipole_moment
    1290              : 
    1291              : ! **************************************************************************************************
    1292              : !> \brief prints final values for classical cell dipole moment calculation
    1293              : !> \param ana_env ...
    1294              : !> \param
    1295              : !> \author Mandes 02.2013
    1296              : ! **************************************************************************************************
    1297            9 :    SUBROUTINE print_dipole_moment(ana_env)
    1298              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1299              : 
    1300            9 :       IF (ana_env%print_test_output) &
    1301            9 :          WRITE (*, *) "TMC|ANALYSIS_FINAL_CLASS_CELL_DIPOLE_MOMENT_X= ", &
    1302           18 :          ana_env%dip_mom%last_dip_cl(:)
    1303            9 :    END SUBROUTINE print_dipole_moment
    1304              : 
    1305              : ! **************************************************************************************************
    1306              : !> \brief calculates the dipole moment analysis
    1307              : !> \param elem ...
    1308              : !> \param weight ...
    1309              : !> \param ana_env ...
    1310              : !> \param
    1311              : !> \author Mandes 03.2013
    1312              : ! **************************************************************************************************
    1313            0 :    SUBROUTINE calc_dipole_analysis(elem, weight, ana_env)
    1314              :       TYPE(tree_type), POINTER                           :: elem
    1315              :       INTEGER                                            :: weight
    1316              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1317              : 
    1318              :       REAL(KIND=dp)                                      :: vol, weight_act
    1319              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: tmp_dip
    1320              :       TYPE(cell_type), POINTER                           :: scaled_cell
    1321              : 
    1322            0 :       NULLIFY (scaled_cell)
    1323              : 
    1324            0 :       CPASSERT(ASSOCIATED(elem))
    1325            0 :       CPASSERT(ASSOCIATED(elem%dipole))
    1326            0 :       CPASSERT(ASSOCIATED(ana_env))
    1327            0 :       CPASSERT(ASSOCIATED(ana_env%dip_ana))
    1328              : 
    1329            0 :       weight_act = weight
    1330            0 :       IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
    1331            0 :          weight_act = weight_act/REAL(8.0, KIND=dp)
    1332              : 
    1333              :       ! get the volume
    1334            0 :       ALLOCATE (scaled_cell)
    1335              :       CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, vol=vol, &
    1336            0 :                            scaled_cell=scaled_cell)
    1337              : 
    1338              :       ! fold exact dipole moments using the classical ones
    1339            0 :       IF (ASSOCIATED(ana_env%dip_mom)) THEN
    1340            0 :          IF (ALL(ana_env%dip_mom%last_dip_cl .NE. elem%dipole)) THEN
    1341              :             elem%dipole = pbc(r=elem%dipole(:) - ana_env%dip_mom%last_dip_cl, &
    1342            0 :                               cell=scaled_cell) + ana_env%dip_mom%last_dip_cl
    1343              :          END IF
    1344              :       END IF
    1345              : 
    1346            0 :       ana_env%dip_ana%conf_counter = ana_env%dip_ana%conf_counter + weight_act
    1347              : 
    1348              :       ! dipole sqared absolut value summed and weight_acted with volume and conf weight_act
    1349              :       ana_env%dip_ana%mu2_pv_s = ana_env%dip_ana%mu2_pv_s + &
    1350            0 :                                  DOT_PRODUCT(elem%dipole(:), elem%dipole(:))/vol*weight_act
    1351              : 
    1352            0 :       tmp_dip(:, :) = 0.0_dp
    1353            0 :       tmp_dip(:, 1) = elem%dipole(:)
    1354              : 
    1355              :       ! dipole sum, weight_acted with volume and conf weight_act
    1356              :       ana_env%dip_ana%mu_pv(:) = ana_env%dip_ana%mu_pv(:) + &
    1357            0 :                                  tmp_dip(:, 1)/vol*weight_act
    1358              : 
    1359              :       ! dipole sum, weight_acted with square root of volume and conf weight_act
    1360              :       ana_env%dip_ana%mu_psv(:) = ana_env%dip_ana%mu_psv(:) + &
    1361            0 :                                   tmp_dip(:, 1)/SQRT(vol)*weight_act
    1362              : 
    1363              :       ! dipole squared sum, weight_acted with volume and conf weight_act
    1364              :       ana_env%dip_ana%mu2_pv(:) = ana_env%dip_ana%mu2_pv(:) + &
    1365            0 :                                   tmp_dip(:, 1)**2/vol*weight_act
    1366              : 
    1367              :       ! calculate the directional average with componentwise correlation per volume
    1368            0 :       tmp_dip(:, :) = MATMUL(tmp_dip(:, :), TRANSPOSE(tmp_dip(:, :)))
    1369              :       ana_env%dip_ana%mu2_pv_mat(:, :) = ana_env%dip_ana%mu2_pv_mat(:, :) + &
    1370            0 :                                          tmp_dip(:, :)/vol*weight_act
    1371              : 
    1372            0 :    END SUBROUTINE calc_dipole_analysis
    1373              : 
    1374              : ! **************************************************************************************************
    1375              : !> \brief prints the actual dipole moment analysis (trajectories)
    1376              : !> \param elem ...
    1377              : !> \param ana_env ...
    1378              : !> \param
    1379              : !> \author Mandes 03.2013
    1380              : ! **************************************************************************************************
    1381            0 :    SUBROUTINE print_act_dipole_analysis(elem, ana_env)
    1382              :       TYPE(tree_type), POINTER                           :: elem
    1383              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1384              : 
    1385              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_tmp
    1386              :       INTEGER                                            :: counter_tmp, file_ptr
    1387              :       LOGICAL                                            :: flag
    1388              :       REAL(KIND=dp)                                      :: diel_const, diel_const_norm, &
    1389              :                                                             diel_const_sym, e0, kB
    1390              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: tmp_dip
    1391              : 
    1392            0 :       kB = boltzmann/joule
    1393            0 :       counter_tmp = INT(ana_env%dip_ana%conf_counter)
    1394              : 
    1395              :       ! TODO get correct constant using physcon
    1396            0 :       e0 = 0.07957747154594767_dp !e^2*a0*me*hbar^-2
    1397            0 :       diel_const_norm = 1/(3.0_dp*e0*kB*ana_env%temperature)
    1398              : 
    1399              :       file_name = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
    1400              :                                         tmc_default_trajectory_file_name, &
    1401            0 :                                         ana_env%temperature)
    1402              :       CALL write_dipoles_in_file(file_name=file_name, &
    1403              :                                  conf_nr=INT(ana_env%dip_ana%conf_counter) + 1, dip=elem%dipole, &
    1404            0 :                                  file_ext="dip_folded")
    1405              : 
    1406              :       ! set output file name
    1407              :       file_name_tmp = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
    1408              :                                             tmc_default_trajectory_file_name, &
    1409            0 :                                             ana_env%temperature)
    1410              : 
    1411            0 :       SELECT CASE (ana_env%dip_ana%ana_type)
    1412              :       CASE (ana_type_default)
    1413              :          file_name = TRIM(expand_file_name_char(file_name_tmp, &
    1414            0 :                                                 "diel_const"))
    1415              :          file_name_tmp = TRIM(expand_file_name_char(file_name_tmp, &
    1416            0 :                                                     "diel_const_tensor"))
    1417              :       CASE (ana_type_sym_xyz)
    1418              :          file_name = TRIM(expand_file_name_char(file_name_tmp, &
    1419            0 :                                                 "diel_const_sym"))
    1420              :          file_name_tmp = TRIM(expand_file_name_char(file_name_tmp, &
    1421            0 :                                                     "diel_const_tensor_sym"))
    1422              :       CASE DEFAULT
    1423            0 :          CPWARN('unknown analysis type "'//cp_to_string(ana_env%dip_ana%ana_type)//'" used.')
    1424              :       END SELECT
    1425              : 
    1426              :       ! calc the dielectric constant
    1427              :       ! 1+( <M^2> - <M>^2 ) / (3*e_0*V*k*T)
    1428              :       diel_const = 1.0_dp + (ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter) - &
    1429              :                              DOT_PRODUCT(ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter), &
    1430              :                                          ana_env%dip_ana%mu_psv(:)/(ana_env%dip_ana%conf_counter)))* &
    1431            0 :                    diel_const_norm
    1432              :       ! symmetrized dielctric constant
    1433              :       ! 1+( <M^2> ) / (3*e_0*V*k*T)
    1434              :       diel_const_sym = 1.0_dp + ana_env%dip_ana%mu2_pv_s/(ana_env%dip_ana%conf_counter)* &
    1435            0 :                        diel_const_norm
    1436              :       ! print dielectric constant trajectory
    1437              :       !  if szmetry used print only every 8th configuration, hence every different (not mirrowed)
    1438            0 :       INQUIRE (FILE=file_name, EXIST=flag)
    1439              :       CALL open_file(file_name=file_name, file_status="UNKNOWN", &
    1440              :                      file_action="WRITE", file_position="APPEND", &
    1441            0 :                      unit_number=file_ptr)
    1442            0 :       IF (.NOT. flag) THEN
    1443            0 :          WRITE (file_ptr, FMT='(A8,5A20)') "# conf", "diel_const", &
    1444            0 :             "diel_const_sym", "diel_const_sym_x", &
    1445            0 :             "diel_const_sym_y", "diel_const_sym_z"
    1446              :       END IF
    1447            0 :       WRITE (file_ptr, FMT="(I8,10F20.10)") counter_tmp, diel_const, &
    1448            0 :          diel_const_sym, &
    1449              :          4.0_dp*PI/(kB*ana_env%temperature)* &
    1450            0 :          ana_env%dip_ana%mu2_pv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
    1451            0 :       CALL close_file(unit_number=file_ptr)
    1452              : 
    1453              :       ! print dielectric constant tensor trajectory
    1454            0 :       INQUIRE (FILE=file_name_tmp, EXIST=flag)
    1455              :       CALL open_file(file_name=file_name_tmp, file_status="UNKNOWN", &
    1456              :                      file_action="WRITE", file_position="APPEND", &
    1457            0 :                      unit_number=file_ptr)
    1458            0 :       IF (.NOT. flag) THEN
    1459            0 :          WRITE (file_ptr, FMT='(A8,9A20)') "# conf", "xx", "xy", "xz", &
    1460            0 :             "yx", "yy", "yz", &
    1461            0 :             "zx", "zy", "zz"
    1462              :       END IF
    1463            0 :       tmp_dip(:, :) = 0.0_dp
    1464            0 :       tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
    1465              : 
    1466            0 :       WRITE (file_ptr, FMT="(I8,10F20.10)") counter_tmp, &
    1467              :          4.0_dp*PI/(kB*ana_env%temperature)* &
    1468              :          (ana_env%dip_ana%mu2_pv_mat(:, :)/REAL(ana_env%dip_ana%conf_counter, KIND=dp) - &
    1469            0 :           MATMUL(tmp_dip(:, :), TRANSPOSE(tmp_dip(:, :))))
    1470            0 :       CALL close_file(unit_number=file_ptr)
    1471            0 :    END SUBROUTINE print_act_dipole_analysis
    1472              : 
    1473              : ! **************************************************************************************************
    1474              : !> \brief prints the dipole moment analysis
    1475              : !> \param ana_env ...
    1476              : !> \param
    1477              : !> \author Mandes 03.2013
    1478              : ! **************************************************************************************************
    1479            0 :    SUBROUTINE print_dipole_analysis(ana_env)
    1480              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1481              : 
    1482              :       CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
    1483              : 
    1484              :       INTEGER                                            :: i
    1485              :       REAL(KIND=dp)                                      :: diel_const_scalar, kB
    1486              :       REAL(KIND=dp), DIMENSION(3)                        :: diel_const_sym, dielec_ev
    1487              :       REAL(KIND=dp), DIMENSION(3, 3)                     :: diel_const, tmp_dip, tmp_ev
    1488              : 
    1489            0 :       kB = boltzmann/joule
    1490              : 
    1491            0 :       CPASSERT(ASSOCIATED(ana_env))
    1492            0 :       CPASSERT(ASSOCIATED(ana_env%dip_ana))
    1493              : 
    1494            0 :       tmp_dip(:, :) = 0.0_dp
    1495              :       diel_const(:, :) = 0.0_dp
    1496            0 :       diel_const_scalar = 0.0_dp
    1497            0 :       diel_const_sym = 0.0_dp
    1498              : 
    1499              :       !dielectric constant
    1500            0 :       tmp_dip(:, 1) = ana_env%dip_ana%mu_psv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
    1501              :       diel_const(:, :) = 4.0_dp*PI/(kB*ana_env%temperature)* &
    1502              :                          (ana_env%dip_ana%mu2_pv_mat(:, :)/REAL(ana_env%dip_ana%conf_counter, KIND=dp) - &
    1503            0 :                           MATMUL(tmp_dip(:, :), TRANSPOSE(tmp_dip(:, :))))
    1504              : 
    1505              :       !dielectric constant for symmetric case
    1506              :       diel_const_sym(:) = 4.0_dp*PI/(kB*ana_env%temperature)* &
    1507            0 :                           ana_env%dip_ana%mu2_pv(:)/REAL(ana_env%dip_ana%conf_counter, KIND=dp)
    1508              : 
    1509            0 :       DO i = 1, 3
    1510            0 :          diel_const(i, i) = diel_const(i, i) + 1.0_dp ! +1 for unpolarizable models, 1.592 for polarizable
    1511            0 :          diel_const_scalar = diel_const_scalar + diel_const(i, i)
    1512              :       END DO
    1513            0 :       diel_const_scalar = diel_const_scalar/REAL(3, KIND=dp)
    1514              : 
    1515            0 :       tmp_dip(:, :) = diel_const
    1516            0 :       CALL diag(3, tmp_dip, dielec_ev, tmp_ev)
    1517              : 
    1518              :       ! print out results
    1519            0 :       WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
    1520            0 :       WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "average dipoles", "-"
    1521            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", cp_to_string(ana_env%temperature)
    1522            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations ", &
    1523            0 :          cp_to_string(REAL(ana_env%dip_ana%conf_counter, KIND=dp))
    1524            0 :       IF (ana_env%dip_ana%ana_type .EQ. ana_type_ice) &
    1525            0 :          WRITE (ana_env%io_unit, FMT='(T2,A,"| ",A)') plabel, &
    1526            0 :          "ice analysis with directions of hexagonal structure"
    1527            0 :       IF (ana_env%dip_ana%ana_type .EQ. ana_type_sym_xyz) &
    1528            0 :          WRITE (ana_env%io_unit, FMT='(T2,A,"| ",A)') plabel, &
    1529            0 :          "ice analysis with symmetrized dipoles in each direction."
    1530              : 
    1531            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "for product of 2 directions(per vol):"
    1532            0 :       DO i = 1, 3
    1533            0 :          WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", ana_env%dip_ana%mu2_pv_mat(i, :)/ &
    1534            0 :             REAL(ana_env%dip_ana%conf_counter, KIND=dp), " |"
    1535              :       END DO
    1536              : 
    1537            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric constant tensor:"
    1538            0 :       DO i = 1, 3
    1539            0 :          WRITE (ana_env%io_unit, '(A,3F16.8,A)') " |", diel_const(i, :), " |"
    1540              :       END DO
    1541              : 
    1542            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric tensor eigenvalues", &
    1543              :          cp_to_string(dielec_ev(1))//" "// &
    1544              :          cp_to_string(dielec_ev(2))//" "// &
    1545            0 :          cp_to_string(dielec_ev(3))
    1546            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric constant symm ", &
    1547              :          cp_to_string(diel_const_sym(1))//" | "// &
    1548              :          cp_to_string(diel_const_sym(2))//" | "// &
    1549            0 :          cp_to_string(diel_const_sym(3))
    1550            0 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "dielectric constant ", &
    1551            0 :          cp_to_string(diel_const_scalar)
    1552            0 :       WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
    1553              : 
    1554            0 :    END SUBROUTINE print_dipole_analysis
    1555              : 
    1556              :    !============================================================================
    1557              :    ! particle displacement in cell (from one configuration to the next)
    1558              :    !============================================================================
    1559              : 
    1560              : ! **************************************************************************************************
    1561              : !> \brief calculates the mean square displacement
    1562              : !> \param elem ...
    1563              : !> \param ana_env ...
    1564              : !> \param
    1565              : !> \author Mandes 02.2013
    1566              : ! **************************************************************************************************
    1567         1000 :    SUBROUTINE calc_displacement(elem, ana_env)
    1568              :       TYPE(tree_type), POINTER                           :: elem
    1569              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1570              : 
    1571              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'calc_displacement'
    1572              : 
    1573              :       CHARACTER(LEN=default_path_length)                 :: file_name, file_name_tmp
    1574              :       INTEGER                                            :: file_ptr, handle, ind
    1575              :       LOGICAL                                            :: flag
    1576              :       REAL(KIND=dp)                                      :: disp
    1577              :       REAL(KIND=dp), DIMENSION(3)                        :: atom_disp
    1578              : 
    1579          500 :       disp = 0.0_dp
    1580              : 
    1581          500 :       CPASSERT(ASSOCIATED(elem))
    1582          500 :       CPASSERT(ASSOCIATED(elem%pos))
    1583          500 :       CPASSERT(ASSOCIATED(ana_env))
    1584          500 :       CPASSERT(ASSOCIATED(ana_env%displace))
    1585          500 :       CPASSERT(ASSOCIATED(ana_env%last_elem))
    1586              : 
    1587              :       ! start the timing
    1588          500 :       CALL timeset(routineN, handle)
    1589              : 
    1590          500 :       DO ind = 1, SIZE(elem%pos), ana_env%dim_per_elem
    1591              :          ! fold into box
    1592        42000 :          atom_disp(:) = elem%pos(ind:ind + 2) - ana_env%last_elem%pos(ind:ind + 2)
    1593              :          CALL get_scaled_cell(cell=ana_env%cell, box_scale=elem%box_scale, &
    1594        10500 :                               vec=atom_disp)
    1595        42000 :          disp = disp + SUM((atom_disp(:)*au2a)**2)
    1596              :       END DO
    1597          500 :       ana_env%displace%disp = ana_env%displace%disp + disp
    1598          500 :       ana_env%displace%conf_counter = ana_env%displace%conf_counter + 1
    1599              : 
    1600          500 :       IF (ana_env%displace%print_disp) THEN
    1601              :          file_name_tmp = expand_file_name_temp(TRIM(ana_env%out_file_prefix)// &
    1602              :                                                tmc_default_trajectory_file_name, &
    1603          500 :                                                ana_env%temperature)
    1604              :          file_name = TRIM(expand_file_name_char(file_name_tmp, &
    1605          500 :                                                 "devi"))
    1606          500 :          INQUIRE (FILE=file_name, EXIST=flag)
    1607              :          CALL open_file(file_name=file_name, file_status="UNKNOWN", &
    1608              :                         file_action="WRITE", file_position="APPEND", &
    1609          500 :                         unit_number=file_ptr)
    1610          500 :          IF (.NOT. flag) &
    1611            3 :             WRITE (file_ptr, *) "# conf     squared deviation of the cell"
    1612          500 :          WRITE (file_ptr, *) elem%nr, disp
    1613          500 :          CALL close_file(unit_number=file_ptr)
    1614              :       END IF
    1615              : 
    1616              :       ! end the timing
    1617          500 :       CALL timestop(handle)
    1618              : 
    1619          500 :    END SUBROUTINE calc_displacement
    1620              : 
    1621              : ! **************************************************************************************************
    1622              : !> \brief prints final values for the displacement calculations
    1623              : !> \param ana_env ...
    1624              : !> \param
    1625              : !> \author Mandes 02.2013
    1626              : ! **************************************************************************************************
    1627            9 :    SUBROUTINE print_average_displacement(ana_env)
    1628              :       TYPE(tmc_analysis_env), POINTER                    :: ana_env
    1629              : 
    1630              :       CHARACTER(LEN=*), PARAMETER :: fmt_my = '(T2,A,"| ",A,T41,A40)', plabel = "TMC_ANA"
    1631              : 
    1632            9 :       WRITE (ana_env%io_unit, FMT="(/,T2,A)") REPEAT("-", 79)
    1633            9 :       WRITE (ana_env%io_unit, FMT="(T2,A,T35,A,T80,A)") "-", "average displacement", "-"
    1634            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "temperature ", &
    1635           18 :          cp_to_string(ana_env%temperature)
    1636            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "used configurations ", &
    1637           18 :          cp_to_string(REAL(ana_env%displace%conf_counter, KIND=dp))
    1638            9 :       WRITE (ana_env%io_unit, FMT=fmt_my) plabel, "cell root mean square deviation: ", &
    1639              :          cp_to_string(SQRT(ana_env%displace%disp/ &
    1640           18 :                            REAL(ana_env%displace%conf_counter, KIND=dp)))
    1641            9 :       IF (ana_env%print_test_output) &
    1642            9 :          WRITE (*, *) "TMC|ANALYSIS_AVERAGE_CELL_DISPLACEMENT_X= ", &
    1643              :          SQRT(ana_env%displace%disp/ &
    1644           18 :               REAL(ana_env%displace%conf_counter, KIND=dp))
    1645            9 :    END SUBROUTINE print_average_displacement
    1646              : END MODULE tmc_analysis
        

Generated by: LCOV version 2.0-1