LCOV - code coverage report
Current view: top level - src/motion - helium_io.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 555 809 68.6 %
Date: 2024-12-12 07:20:04 Functions: 13 14 92.9 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief  I/O subroutines for helium
      10             : !> \author Lukasz Walewski
      11             : !> \date   2009-06-08
      12             : ! **************************************************************************************************
      13             : MODULE helium_io
      14             : 
      15             :    USE cell_types,                      ONLY: get_cell
      16             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      17             :                                               cp_logger_get_default_unit_nr,&
      18             :                                               cp_logger_type
      19             :    USE cp_output_handling,              ONLY: cp_p_file,&
      20             :                                               cp_print_key_finished_output,&
      21             :                                               cp_print_key_should_output,&
      22             :                                               cp_print_key_unit_nr
      23             :    USE cp_parser_methods,               ONLY: parser_get_next_line,&
      24             :                                               parser_get_object
      25             :    USE cp_parser_types,                 ONLY: cp_parser_type,&
      26             :                                               parser_create,&
      27             :                                               parser_release
      28             :    USE cp_units,                        ONLY: cp_unit_from_cp2k,&
      29             :                                               cp_unit_to_cp2k
      30             :    USE helium_common,                   ONLY: helium_cycle_number,&
      31             :                                               helium_cycle_of,&
      32             :                                               helium_is_winding,&
      33             :                                               helium_path_length,&
      34             :                                               helium_pbc
      35             :    USE helium_interactions,             ONLY: helium_total_inter_action,&
      36             :                                               helium_total_link_action,&
      37             :                                               helium_total_pair_action
      38             :    USE helium_types,                    ONLY: &
      39             :         e_id_interact, e_id_kinetic, e_id_potential, e_id_thermo, e_id_total, e_id_virial, &
      40             :         helium_solvent_p_type, helium_solvent_type, rho_num
      41             :    USE input_constants,                 ONLY: &
      42             :         fmt_id_pdb, fmt_id_xyz, helium_cell_shape_cube, helium_cell_shape_octahedron, &
      43             :         helium_sampling_ceperley, helium_sampling_worm, helium_solute_intpot_mwater, &
      44             :         helium_solute_intpot_nnp, helium_solute_intpot_none, perm_cycle, perm_plain
      45             :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      46             :                                               section_vals_type,&
      47             :                                               section_vals_val_get
      48             :    USE kinds,                           ONLY: default_path_length,&
      49             :                                               default_string_length,&
      50             :                                               dp,&
      51             :                                               int_8
      52             :    USE machine,                         ONLY: m_flush
      53             :    USE memory_utilities,                ONLY: reallocate
      54             :    USE message_passing,                 ONLY: mp_para_env_type
      55             :    USE physcon,                         ONLY: angstrom,&
      56             :                                               massunit
      57             :    USE pint_types,                      ONLY: pint_env_type
      58             : #include "../base/base_uses.f90"
      59             : 
      60             :    IMPLICIT NONE
      61             : 
      62             :    PRIVATE
      63             : 
      64             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      65             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_io'
      66             : 
      67             :    CHARACTER(len=*), DIMENSION(6), PARAMETER, PRIVATE :: m_dist_name = (/ &
      68             :                                                          "SINGLEV    ", &
      69             :                                                          "UNIFORM    ", &
      70             :                                                          "LINEAR     ", &
      71             :                                                          "QUADRATIC  ", &
      72             :                                                          "EXPONENTIAL", &
      73             :                                                          "GAUSSIAN   "/)
      74             : 
      75             :    PUBLIC :: helium_read_xyz
      76             :    PUBLIC :: helium_print_rdf
      77             :    PUBLIC :: helium_print_rho
      78             :    PUBLIC :: helium_write_line
      79             :    PUBLIC :: helium_write_setup
      80             :    PUBLIC :: helium_print_energy
      81             :    PUBLIC :: helium_print_vector
      82             :    PUBLIC :: helium_print_plength
      83             :    PUBLIC :: helium_print_coordinates
      84             :    PUBLIC :: helium_print_force
      85             :    PUBLIC :: helium_print_accepts
      86             :    PUBLIC :: helium_print_perm
      87             :    PUBLIC :: helium_print_action
      88             :    PUBLIC :: helium_write_cubefile
      89             : 
      90             : CONTAINS
      91             : 
      92             : ! ***************************************************************************
      93             : !> \brief  Read XYZ coordinates from file
      94             : !> \param coords ...
      95             : !> \param file_name ...
      96             : !> \param para_env ...
      97             : !> \date   2009-06-03
      98             : !> \author Lukasz Walewski
      99             : !> \note   This is a parallel routine, all ranks get coords defined
     100             : ! **************************************************************************************************
     101           0 :    SUBROUTINE helium_read_xyz(coords, file_name, para_env)
     102             : 
     103             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: coords
     104             :       CHARACTER(LEN=default_path_length)                 :: file_name
     105             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     106             : 
     107             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_read_xyz'
     108             : 
     109             :       CHARACTER(LEN=default_string_length)               :: strtmp
     110             :       INTEGER                                            :: frame, handle, istat, j, natom
     111             :       LOGICAL                                            :: found, my_end
     112             :       TYPE(cp_parser_type)                               :: parser
     113             : 
     114           0 :       CALL timeset(routineN, handle)
     115             : 
     116             :       ! check if the file can be accessed
     117           0 :       INQUIRE (FILE=file_name, EXIST=found, IOSTAT=istat)
     118           0 :       IF (istat /= 0) THEN
     119             :          WRITE (UNIT=strtmp, FMT="(A,I0,A)") &
     120             :             "An error occurred inquiring the file <"// &
     121           0 :             TRIM(file_name)//">"
     122           0 :          CPABORT(TRIM(strtmp))
     123             :       END IF
     124           0 :       IF (.NOT. found) THEN
     125           0 :          CPABORT("Coordinate file <"//TRIM(file_name)//"> not found.")
     126             :       END IF
     127             : 
     128             :       CALL parser_create( &
     129             :          parser, &
     130             :          file_name, &
     131             :          para_env=para_env, &
     132           0 :          parse_white_lines=.TRUE.)
     133             : 
     134             :       natom = 0
     135           0 :       frame = 0
     136           0 :       CALL parser_get_next_line(parser, 1)
     137           0 :       Frames: DO
     138             :          ! Atom number
     139           0 :          CALL parser_get_object(parser, natom)
     140           0 :          frame = frame + 1
     141           0 :          IF (frame == 1) THEN
     142           0 :             ALLOCATE (coords(3*natom))
     143             :          ELSE
     144           0 :             strtmp = "Warning: more than one frame on file <"//TRIM(file_name)//">"
     145           0 :             CALL helium_write_line(strtmp)
     146           0 :             strtmp = "Warning: using the first frame only!"
     147           0 :             CALL helium_write_line(strtmp)
     148           0 :             EXIT Frames
     149             :          END IF
     150             :          ! Dummy line
     151           0 :          CALL parser_get_next_line(parser, 2)
     152           0 :          DO j = 1, natom
     153             :             ! Atom coordinates
     154           0 :             READ (parser%input_line, *) strtmp, &
     155           0 :                coords(3*(j - 1) + 1), &
     156           0 :                coords(3*(j - 1) + 2), &
     157           0 :                coords(3*(j - 1) + 3)
     158           0 :             coords(3*(j - 1) + 1) = cp_unit_to_cp2k(coords(3*(j - 1) + 1), "angstrom")
     159           0 :             coords(3*(j - 1) + 2) = cp_unit_to_cp2k(coords(3*(j - 1) + 2), "angstrom")
     160           0 :             coords(3*(j - 1) + 3) = cp_unit_to_cp2k(coords(3*(j - 1) + 3), "angstrom")
     161             :             ! If there's a white line or end of file exit.. otherwise go on
     162           0 :             CALL parser_get_next_line(parser, 1, at_end=my_end)
     163           0 :             my_end = my_end .OR. (LEN_TRIM(parser%input_line) == 0)
     164           0 :             IF (my_end) THEN
     165           0 :                IF (j /= natom) THEN
     166           0 :                   CPABORT("Error in XYZ format.")
     167             :                END IF
     168             :                EXIT Frames
     169             :             END IF
     170             :          END DO
     171             :       END DO Frames
     172           0 :       CALL parser_release(parser)
     173             : 
     174           0 :       CALL timestop(handle)
     175             : 
     176           0 :    END SUBROUTINE helium_read_xyz
     177             : 
     178             : ! ***************************************************************************
     179             : !> \brief  Write helium parameters to the output unit
     180             : !> \param helium ...
     181             : !> \date   2009-06-03
     182             : !> \par    History
     183             : !>         2023-07-23 Modified to work with NNP solute-solvent interactions [lduran]
     184             : !> \author Lukasz Walewski
     185             : ! **************************************************************************************************
     186          25 :    SUBROUTINE helium_write_setup(helium)
     187             : 
     188             :       TYPE(helium_solvent_type), POINTER                 :: helium
     189             : 
     190             :       CHARACTER(len=default_string_length)               :: msg_str, my_label, stmp, stmp1, stmp2, &
     191             :                                                             unit_str
     192             :       INTEGER                                            :: i, itmp, unit_nr
     193             :       INTEGER(KIND=int_8)                                :: i8tmp
     194             :       REAL(KIND=dp)                                      :: rtmp, v1, v2, v3
     195             :       REAL(KIND=dp), DIMENSION(3)                        :: my_abc
     196             :       TYPE(cp_logger_type), POINTER                      :: logger
     197             : 
     198          25 :       NULLIFY (logger)
     199          25 :       logger => cp_get_default_logger()
     200          25 :       my_label = "HELIUM| "
     201             : 
     202          25 :       IF (logger%para_env%is_source()) THEN
     203          13 :          unit_nr = cp_logger_get_default_unit_nr(logger)
     204             : 
     205          13 :          WRITE (unit_nr, *)
     206             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     207          13 :             " Number of helium environments:     ", helium%num_env
     208             : 
     209             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     210          13 :             " Number of solvent atoms:           ", helium%atoms
     211             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     212          13 :             " Number of solvent beads:           ", helium%beads
     213             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     214          13 :             " Total number of solvent particles: ", helium%atoms*helium%beads
     215             : 
     216          13 :          unit_str = "angstrom^-3"
     217             :          rtmp = cp_unit_from_cp2k(helium%density, &
     218          13 :                                   unit_str)
     219             :          WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" Density   ["// &
     220          13 :             TRIM(unit_str)//"]:", rtmp
     221             : 
     222          13 :          unit_str = "a.m.u."
     223          13 :          rtmp = helium%he_mass_au/massunit
     224             :          WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" He Mass   ["// &
     225          13 :             TRIM(unit_str)//"]: ", rtmp
     226             : 
     227          13 :          unit_str = "angstrom"
     228             :          rtmp = cp_unit_from_cp2k(helium%cell_size, &
     229          13 :                                   unit_str)
     230             :          WRITE (unit_nr, '(T2,A,F12.6)') TRIM(my_label)//" Cell size ["// &
     231          13 :             TRIM(unit_str)//"]:   ", rtmp
     232             : 
     233          13 :          IF (helium%periodic) THEN
     234           8 :             SELECT CASE (helium%cell_shape)
     235             :             CASE (helium_cell_shape_cube)
     236           0 :                CALL helium_write_line("PBC cell shape: CUBE.")
     237             :             CASE (helium_cell_shape_octahedron)
     238           8 :                CALL helium_write_line("PBC cell shape: TRUNCATED OCTAHEDRON.")
     239             :             CASE DEFAULT
     240           8 :                CALL helium_write_line("*** Warning: unknown cell shape.")
     241             :             END SELECT
     242             :          ELSE
     243           5 :             CALL helium_write_line("PBC turned off.")
     244             :          END IF
     245             : 
     246          13 :          IF (helium%droplet_radius .LT. HUGE(1.0_dp)) THEN
     247           5 :             WRITE (stmp, *) helium%droplet_radius*angstrom
     248           5 :             CALL helium_write_line("Droplet radius: "//TRIM(ADJUSTL(stmp))//" [angstrom]")
     249             :          END IF
     250             : 
     251             :          ! first step gets incremented during first iteration
     252             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     253          13 :             " First MC step                      :", helium%first_step + 1
     254             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     255          13 :             " Last MC step                       :", helium%last_step
     256             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     257          13 :             " Total number of MC steps           :", helium%num_steps
     258             :          WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     259          13 :             " Number of outer MC trials per step :", helium%iter_rot
     260          13 :          i8tmp = INT(helium%iter_rot, int_8)
     261          13 :          IF (helium%sampling_method == helium_sampling_ceperley) THEN
     262             :             WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     263           8 :                " Number of inner MC trials per step :", helium%iter_norot
     264           8 :             i8tmp = i8tmp*INT(helium%iter_norot, int_8)
     265             :          END IF
     266          13 :          stmp = ""
     267          13 :          WRITE (stmp, *) i8tmp
     268             :          WRITE (unit_nr, '(T2,A)') TRIM(my_label)// &
     269          13 :             " Total number of MC trials per step : "//TRIM(ADJUSTL(stmp))
     270          13 :          i8tmp = INT(helium%num_steps, int_8)
     271          13 :          i8tmp = i8tmp*INT(helium%iter_rot, int_8)
     272          13 :          IF (helium%sampling_method == helium_sampling_ceperley) THEN
     273           8 :             i8tmp = i8tmp*INT(helium%iter_norot, int_8)
     274             :          END IF
     275          13 :          stmp = ""
     276          13 :          WRITE (stmp, *) i8tmp
     277             :          WRITE (unit_nr, '(T2,A)') TRIM(my_label)// &
     278          13 :             " Total number of MC trials          : "//TRIM(ADJUSTL(stmp))
     279             : 
     280          21 :          SELECT CASE (helium%sampling_method)
     281             : 
     282             :          CASE (helium_sampling_ceperley)
     283             : 
     284             :             ! permutation cycle length sampling
     285           8 :             stmp = ""
     286           8 :             CALL helium_write_line(stmp)
     287           8 :             WRITE (stmp, *) helium%maxcycle
     288           8 :             stmp2 = ""
     289             :             WRITE (stmp2, *) "Using maximum permutation cycle length: "// &
     290           8 :                TRIM(ADJUSTL(stmp))
     291           8 :             CALL helium_write_line(stmp2)
     292           8 :             stmp = ""
     293             :             WRITE (stmp, *) "Permutation cycle length distribution: "// &
     294           8 :                TRIM(ADJUSTL(m_dist_name(helium%m_dist_type)))
     295           8 :             CALL helium_write_line(stmp)
     296           8 :             stmp = ""
     297           8 :             stmp1 = ""
     298           8 :             WRITE (stmp1, *) helium%m_ratio
     299           8 :             stmp2 = ""
     300           8 :             WRITE (stmp2, *) helium%m_value
     301             :             WRITE (stmp, *) "Using ratio "//TRIM(ADJUSTL(stmp1))// &
     302           8 :                " for M = "//TRIM(ADJUSTL(stmp2))
     303           8 :             CALL helium_write_line(stmp)
     304           8 :             stmp = ""
     305           8 :             CALL helium_write_line(stmp)
     306             : 
     307             :          CASE (helium_sampling_worm)
     308             : 
     309           5 :             stmp1 = ""
     310           5 :             stmp2 = ""
     311           5 :             CALL helium_write_line(stmp1)
     312           5 :             WRITE (stmp1, *) helium%worm_centroid_drmax
     313             :             WRITE (stmp2, *) "WORM| Centroid move max. displacement: "// &
     314           5 :                TRIM(ADJUSTL(stmp1))
     315           5 :             CALL helium_write_line(stmp2)
     316           5 :             stmp1 = ""
     317           5 :             stmp2 = ""
     318           5 :             WRITE (stmp1, *) helium%worm_staging_l
     319           5 :             WRITE (stmp2, *) "WORM| Maximal staging length: "//TRIM(ADJUSTL(stmp1))
     320           5 :             CALL helium_write_line(stmp2)
     321           5 :             stmp1 = ""
     322           5 :             stmp2 = ""
     323           5 :             WRITE (stmp1, *) helium%worm_open_close_scale
     324           5 :             WRITE (stmp2, *) "WORM| Open/Close scaling: "//TRIM(ADJUSTL(stmp1))
     325           5 :             CALL helium_write_line(stmp2)
     326           5 :             stmp1 = ""
     327           5 :             stmp2 = ""
     328           5 :             WRITE (stmp1, *) helium%worm_allow_open
     329           5 :             WRITE (stmp2, *) "WORM| Open worm sector: "//TRIM(ADJUSTL(stmp1))
     330           5 :             CALL helium_write_line(stmp2)
     331           5 :             stmp1 = ""
     332           5 :             stmp2 = ""
     333           5 :             WRITE (stmp1, *) helium%worm_show_statistics
     334           5 :             WRITE (stmp2, *) "WORM| Print statistics: "//TRIM(ADJUSTL(stmp1))
     335           5 :             CALL helium_write_line(stmp2)
     336           5 :             IF (helium%worm_max_open_cycles > 0 .AND. helium%worm_allow_open) THEN
     337           0 :                stmp1 = ""
     338           0 :                stmp2 = ""
     339           0 :                WRITE (stmp1, *) helium%worm_max_open_cycles
     340           0 :                WRITE (stmp2, *) "WORM| Max. nb of MC cycles in open sector: "//TRIM(ADJUSTL(stmp1))
     341           0 :                CALL helium_write_line(stmp2)
     342             :             END IF
     343           5 :             stmp1 = ""
     344           5 :             CALL helium_write_line(stmp1)
     345             : 
     346             :          CASE DEFAULT
     347           0 :             WRITE (msg_str, *) helium%sampling_method
     348           0 :             msg_str = "Sampling method ("//TRIM(ADJUSTL(msg_str))//") not supported."
     349          13 :             CPABORT(msg_str)
     350             : 
     351             :          END SELECT
     352             : 
     353             :          ! solute related data
     354          13 :          IF (helium%solute_present) THEN
     355             :             WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     356           8 :                " Number of solute atoms:            ", helium%solute_atoms
     357             :             WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     358           8 :                " Number of solute beads:            ", helium%solute_beads
     359             :             WRITE (unit_nr, '(T2,A,1X,I0)') TRIM(my_label)// &
     360           8 :                " Total number of solute particles:  ", helium%solute_atoms* &
     361          16 :                helium%solute_beads
     362             : 
     363           8 :             stmp1 = ""
     364           8 :             SELECT CASE (helium%solute_interaction)
     365             :             CASE (helium_solute_intpot_none)
     366           0 :                WRITE (stmp1, *) "NONE"
     367             :             CASE (helium_solute_intpot_mwater)
     368           7 :                WRITE (stmp1, *) "MWATER"
     369             :             CASE (helium_solute_intpot_nnp)
     370           1 :                WRITE (stmp1, *) "NNP"
     371             :             CASE DEFAULT
     372           8 :                WRITE (stmp1, *) "***UNKNOWN***"
     373             :             END SELECT
     374             :             WRITE (unit_nr, '(T2,A,A,A)') &
     375           8 :                TRIM(my_label), &
     376           8 :                " Solute interaction type: ", &
     377          16 :                TRIM(ADJUSTL(stmp1))
     378             : 
     379           8 :             CALL get_cell(helium%solute_cell, abc=my_abc)
     380           8 :             unit_str = "angstrom"
     381           8 :             v1 = cp_unit_from_cp2k(my_abc(1), unit_str)
     382           8 :             v2 = cp_unit_from_cp2k(my_abc(2), unit_str)
     383           8 :             v3 = cp_unit_from_cp2k(my_abc(3), unit_str)
     384             :             WRITE (unit_nr, '(T2,A,F12.6,1X,F12.6,1X,F12.6)') &
     385             :                TRIM(my_label)//" Solute cell size ["// &
     386           8 :                TRIM(unit_str)//"]:   ", v1, v2, v3
     387             :          ELSE
     388           5 :             WRITE (unit_nr, '(T2,A)') TRIM(my_label)//" Solute is NOT present"
     389             :          END IF
     390             :       END IF
     391          25 :       CALL helium_write_line("")
     392             : 
     393             :       ! RDF related settings
     394          25 :       IF (helium%rdf_present) THEN
     395           2 :          WRITE (stmp, '(I6)') helium%rdf_num
     396           2 :          CALL helium_write_line("RDF| number of centers: "//TRIM(stmp))
     397           2 :          rtmp = cp_unit_from_cp2k(helium%rdf_delr, "angstrom")
     398           2 :          WRITE (stmp, '(1X,F12.6)') rtmp
     399           2 :          CALL helium_write_line("RDF| delr [angstrom]  : "//TRIM(stmp))
     400           2 :          rtmp = cp_unit_from_cp2k(helium%rdf_maxr, "angstrom")
     401           2 :          WRITE (stmp, '(1X,F12.6)') rtmp
     402           2 :          CALL helium_write_line("RDF| maxr [angstrom]  : "//TRIM(stmp))
     403           2 :          itmp = helium%rdf_nbin
     404           2 :          WRITE (stmp, '(I6)') itmp
     405           2 :          CALL helium_write_line("RDF| nbin             : "//TRIM(stmp))
     406           2 :          CALL helium_write_line("")
     407             :       ELSE
     408          23 :          CALL helium_write_line("RDF| radial distributions will NOT be calculated.")
     409          23 :          CALL helium_write_line("")
     410             :       END IF
     411             : 
     412             :       ! RHO related settings
     413          25 :       IF (helium%rho_present) THEN
     414           2 :          CALL helium_write_line('RHO| The following densities will be calculated:')
     415          12 :          DO i = 1, rho_num
     416          12 :             IF (helium%rho_property(i)%is_calculated) THEN
     417           2 :                WRITE (stmp, '(A)') 'RHO|    '//TRIM(helium%rho_property(i)%name)
     418           2 :                CALL helium_write_line(TRIM(stmp))
     419             :             END IF
     420             :          END DO
     421           2 :          CALL helium_write_line('RHO|')
     422           2 :          rtmp = cp_unit_from_cp2k(helium%rho_delr, "angstrom")
     423           2 :          WRITE (stmp, '(1X,F12.6)') rtmp
     424           2 :          CALL helium_write_line("RHO| delr [angstrom]  : "//TRIM(stmp))
     425           2 :          rtmp = cp_unit_from_cp2k(helium%rho_maxr, "angstrom")
     426           2 :          WRITE (stmp, '(1X,F12.6)') rtmp
     427           2 :          CALL helium_write_line("RHO| maxr [angstrom]  : "//TRIM(stmp))
     428           2 :          itmp = helium%rho_nbin
     429           2 :          WRITE (stmp, '(I6)') itmp
     430           2 :          CALL helium_write_line("RHO| nbin             : "//TRIM(stmp))
     431           2 :          CALL helium_write_line("")
     432             :       ELSE
     433          23 :          CALL helium_write_line("RHO| Density distributions will NOT be calculated.")
     434          23 :          CALL helium_write_line("")
     435             :       END IF
     436             : 
     437          25 :       RETURN
     438             :    END SUBROUTINE helium_write_setup
     439             : 
     440             : ! ***************************************************************************
     441             : !> \brief  Writes out a line of text to the default output unit
     442             : !> \param line ...
     443             : !> \date   2009-07-10
     444             : !> \author Lukasz Walewski
     445             : ! **************************************************************************************************
     446         897 :    SUBROUTINE helium_write_line(line)
     447             : 
     448             :       CHARACTER(len=*), INTENT(IN)                       :: line
     449             : 
     450             :       CHARACTER(len=default_string_length)               :: my_label
     451             :       INTEGER                                            :: unit_nr
     452             :       TYPE(cp_logger_type), POINTER                      :: logger
     453             : 
     454         897 :       NULLIFY (logger)
     455         897 :       logger => cp_get_default_logger()
     456         897 :       my_label = "HELIUM|"
     457             : 
     458         897 :       IF (logger%para_env%is_source()) THEN
     459         537 :          unit_nr = cp_logger_get_default_unit_nr(logger)
     460         537 :          WRITE (unit_nr, '(T2,A)') TRIM(my_label)//" "//TRIM(line)
     461             :       END IF
     462             : 
     463         897 :       RETURN
     464             :    END SUBROUTINE helium_write_line
     465             : 
     466             : ! ***************************************************************************
     467             : !> \brief  Print energies according to HELIUM%PRINT%ENERGY
     468             : !> \param helium_env ...
     469             : !> \date   2009-06-08
     470             : !> \par    History
     471             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     472             : !> \author Lukasz Walewski
     473             : ! **************************************************************************************************
     474         105 :    SUBROUTINE helium_print_energy(helium_env)
     475             : 
     476             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     477             : 
     478             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_energy'
     479             : 
     480             :       INTEGER                                            :: handle, iteration, k, m, unit_nr
     481             :       LOGICAL                                            :: cepsample, file_is_new, should_output
     482             :       REAL(KIND=dp)                                      :: naccptd, ntot
     483         105 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: my_energy
     484             :       TYPE(cp_logger_type), POINTER                      :: logger
     485             :       TYPE(section_vals_type), POINTER                   :: print_key
     486             : 
     487         105 :       CALL timeset(routineN, handle)
     488             : 
     489         105 :       NULLIFY (logger, print_key)
     490         105 :       logger => cp_get_default_logger()
     491             : 
     492         105 :       IF (logger%para_env%is_source()) THEN
     493             :          print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
     494          55 :                                                  "MOTION%PINT%HELIUM%PRINT%ENERGY")
     495             :          should_output = BTEST(cp_print_key_should_output( &
     496             :                                iteration_info=logger%iter_info, &
     497          55 :                                basis_section=print_key), cp_p_file)
     498          55 :          cepsample = helium_env(1)%helium%sampling_method == helium_sampling_ceperley
     499             :       END IF
     500         105 :       CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
     501         105 :       CALL helium_env(1)%comm%bcast(cepsample, logger%para_env%source)
     502             : 
     503         105 :       IF (should_output) THEN
     504             : 
     505         105 :          naccptd = 0.0_dp
     506         105 :          IF (cepsample) THEN
     507          80 :             ntot = 0.0_dp
     508         190 :             DO k = 1, SIZE(helium_env)
     509         110 :                ntot = ntot + helium_env(1)%helium%iter_norot*helium_env(1)%helium%iter_rot
     510         480 :                DO m = 1, helium_env(k)%helium%maxcycle
     511         400 :                   naccptd = naccptd + helium_env(k)%helium%num_accepted(helium_env(k)%helium%bisctlog2 + 2, m)
     512             :                END DO
     513             :             END DO
     514             :          ELSE !(wormsample)
     515          25 :             ntot = 0.0_dp
     516          50 :             DO k = 1, SIZE(helium_env)
     517          25 :                naccptd = naccptd + helium_env(k)%helium%num_accepted(1, 1)
     518          50 :                ntot = ntot + helium_env(k)%helium%num_accepted(2, 1)
     519             :             END DO
     520             :          END IF
     521         105 :          CALL helium_env(1)%comm%sum(naccptd)
     522         105 :          CALL helium_env(1)%comm%sum(ntot)
     523             : 
     524         105 :          IF (logger%para_env%is_source()) THEN
     525          55 :             my_energy => helium_env(1)%helium%energy_avrg
     526          55 :             iteration = logger%iter_info%iteration(2)
     527             : 
     528             :             unit_nr = cp_print_key_unit_nr( &
     529             :                       logger, &
     530             :                       print_key, &
     531             :                       middle_name="helium-energy", &
     532             :                       extension=".dat", &
     533          55 :                       is_new_file=file_is_new)
     534             : 
     535          55 :             IF (file_is_new) THEN
     536             :                WRITE (unit_nr, '(A9,7(1X,A20))') &
     537          11 :                   "#    Step", &
     538          11 :                   "            AccRatio", &
     539          11 :                   "               E_pot", &
     540          11 :                   "               E_kin", &
     541          11 :                   "            E_thermo", &
     542          11 :                   "            E_virial", &
     543          11 :                   "             E_inter", &
     544          22 :                   "               E_tot"
     545             :             END IF
     546             : 
     547             :             WRITE (unit_nr, "(I9,7(1X,F20.9))") &
     548          55 :                iteration, &
     549          55 :                naccptd/ntot, &
     550          55 :                my_energy(e_id_potential), &
     551          55 :                my_energy(e_id_kinetic), &
     552          55 :                my_energy(e_id_thermo), &
     553          55 :                my_energy(e_id_virial), &
     554          55 :                my_energy(e_id_interact), &
     555         110 :                my_energy(e_id_total)
     556             : 
     557          55 :             CALL m_flush(unit_nr)
     558          55 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     559             : 
     560             :          END IF
     561             :       END IF
     562             : 
     563         105 :       CALL timestop(handle)
     564             : 
     565         105 :       RETURN
     566         105 :    END SUBROUTINE helium_print_energy
     567             : 
     568             : ! ***************************************************************************
     569             : !> \brief  Print a 3D real vector according to printkey <pkey>
     570             : !> \param helium_env ...
     571             : !> \param pkey ...
     572             : !> \param DATA ...
     573             : !> \param uconv ...
     574             : !> \param col_label ...
     575             : !> \param cmmnt ...
     576             : !> \param fname ...
     577             : !> \param fpos ...
     578             : !> \param avg ...
     579             : !> \date   2014-09-09
     580             : !> \par    History
     581             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     582             : !> \author Lukasz Walewski
     583             : ! **************************************************************************************************
     584         630 :    SUBROUTINE helium_print_vector(helium_env, pkey, DATA, uconv, col_label, cmmnt, fname, fpos, avg)
     585             : 
     586             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     587             :       CHARACTER(len=*)                                   :: pkey
     588             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: DATA
     589             :       REAL(KIND=dp)                                      :: uconv
     590             :       CHARACTER(len=*), DIMENSION(3)                     :: col_label
     591             :       CHARACTER(len=*)                                   :: cmmnt, fname
     592             :       CHARACTER(len=*), OPTIONAL                         :: fpos
     593             :       LOGICAL, OPTIONAL                                  :: avg
     594             : 
     595             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_vector'
     596             : 
     597             :       CHARACTER(len=default_string_length)               :: my_fpos
     598             :       INTEGER                                            :: handle, i, irank, msglen, nenv, offset, &
     599             :                                                             unit_nr
     600             :       LOGICAL                                            :: is_new, my_avg, should_output
     601             :       REAL(KIND=dp), DIMENSION(3)                        :: avg_data
     602         630 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: data_p
     603             :       TYPE(cp_logger_type), POINTER                      :: logger
     604             :       TYPE(section_vals_type), POINTER                   :: psection
     605             : 
     606         630 :       CALL timeset(routineN, handle)
     607             : 
     608         630 :       IF (PRESENT(avg)) THEN
     609         315 :          my_avg = avg
     610             :       ELSE
     611             :          my_avg = .FALSE.
     612             :       END IF
     613             : 
     614         630 :       IF (PRESENT(fpos)) THEN
     615         315 :          my_fpos = fpos
     616             :       ELSE
     617         315 :          my_fpos = "APPEND"
     618             :       END IF
     619             : 
     620         630 :       NULLIFY (logger, psection)
     621         630 :       logger => cp_get_default_logger()
     622             : 
     623         630 :       psection => section_vals_get_subs_vals(helium_env(1)%helium%input, pkey)
     624             :       should_output = BTEST(cp_print_key_should_output( &
     625             :                             iteration_info=logger%iter_info, &
     626         630 :                             basis_section=psection), cp_p_file)
     627             : 
     628         630 :       IF (.NOT. should_output) THEN
     629         315 :          CALL timestop(handle)
     630         315 :          RETURN
     631             :       END IF
     632             : 
     633         315 :       IF (my_avg) THEN
     634             :          ! average data over all processors and environments
     635         315 :          avg_data(:) = 0.0_dp
     636         315 :          msglen = SIZE(avg_data)
     637         720 :          DO i = 0, SIZE(helium_env) - 1
     638        1935 :             avg_data(:) = avg_data(:) + DATA(msglen*i + 1:msglen*(i + 1))
     639             :          END DO
     640         315 :          CALL helium_env(1)%comm%sum(avg_data)
     641         315 :          nenv = helium_env(1)%helium%num_env
     642        1260 :          avg_data(:) = avg_data(:)/REAL(nenv, dp)
     643             :       ELSE
     644             :          ! gather data from all processors
     645           0 :          offset = 0
     646           0 :          DO i = 1, logger%para_env%mepos
     647           0 :             offset = offset + helium_env(1)%env_all(i)
     648             :          END DO
     649             : 
     650           0 :          helium_env(1)%helium%rtmp_3_np_1d = 0.0_dp
     651           0 :          msglen = SIZE(avg_data)
     652           0 :          DO i = 0, SIZE(helium_env) - 1
     653           0 :             helium_env(1)%helium%rtmp_3_np_1d(msglen*(offset + i) + 1:msglen*(offset + i + 1)) = DATA(msglen*i + 1:msglen*(i + 1))
     654             :          END DO
     655           0 :          CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_np_1d)
     656             :       END IF
     657             : 
     658             :       unit_nr = cp_print_key_unit_nr( &
     659             :                 logger, &
     660             :                 psection, &
     661             :                 middle_name=fname, &
     662             :                 extension=".dat", &
     663             :                 file_position=my_fpos, &
     664         315 :                 is_new_file=is_new)
     665             : 
     666         315 :       IF (logger%para_env%is_source()) THEN
     667             : 
     668         165 :          IF (is_new) THEN
     669             :             ! comment
     670         165 :             IF (cmmnt .NE. "") THEN
     671         165 :                WRITE (unit_nr, '(A)') "# "//cmmnt
     672             :             END IF
     673             :             ! column labels
     674             :             WRITE (unit_nr, '(A2,A18,1X,A20,1X,A20)') &
     675         165 :                "# ", &
     676         165 :                col_label(1), &
     677         165 :                col_label(2), &
     678         330 :                col_label(3)
     679             :          END IF
     680             : 
     681         165 :          IF (my_avg) THEN
     682         660 :             DO i = 1, 3
     683         495 :                WRITE (unit_nr, '(E27.20)', ADVANCE='NO') uconv*avg_data(i)
     684         660 :                IF (i .LT. 3) THEN
     685         330 :                   WRITE (unit_nr, '(1X)', ADVANCE='NO')
     686             :                END IF
     687             :             END DO
     688         165 :             WRITE (unit_nr, '(A)') ""
     689             :          ELSE
     690             :             ! iterate over processors/helium environments
     691           0 :             DO irank = 1, helium_env(1)%helium%num_env
     692             :                ! unpack data (actually point to the right fragment only)
     693           0 :                msglen = SIZE(avg_data)
     694           0 :                offset = (irank - 1)*msglen
     695           0 :                data_p => helium_env(1)%helium%rtmp_3_np_1d(offset + 1:offset + msglen)
     696             :                ! write out the data
     697           0 :                DO i = 1, 3
     698           0 :                   WRITE (unit_nr, '(E27.20)', ADVANCE='NO') uconv*data_p(i)
     699           0 :                   IF (i .LT. 3) THEN
     700           0 :                      WRITE (unit_nr, '(1X)', ADVANCE='NO')
     701             :                   END IF
     702             :                END DO
     703           0 :                WRITE (unit_nr, '(A)') ""
     704             :             END DO
     705             :          END IF
     706             : 
     707         165 :          CALL m_flush(unit_nr)
     708         165 :          CALL cp_print_key_finished_output(unit_nr, logger, psection)
     709             : 
     710             :       END IF
     711             : 
     712         315 :       CALL timestop(handle)
     713             : 
     714         315 :       RETURN
     715         630 :    END SUBROUTINE helium_print_vector
     716             : 
     717             : ! ***************************************************************************
     718             : !> \brief  Print acceptance counts according to HELIUM%PRINT%ACCEPTS
     719             : !> \param helium_env ...
     720             : !> \date   2010-05-27
     721             : !> \par    History
     722             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     723             : !> \author Lukasz Walewski
     724             : ! **************************************************************************************************
     725         105 :    SUBROUTINE helium_print_accepts(helium_env)
     726             : 
     727             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     728             : 
     729             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_accepts'
     730             : 
     731             :       INTEGER                                            :: handle, i, j, unit_nr
     732             :       LOGICAL                                            :: file_is_new, should_output
     733             :       TYPE(cp_logger_type), POINTER                      :: logger
     734             :       TYPE(section_vals_type), POINTER                   :: print_key
     735             : 
     736         105 :       CALL timeset(routineN, handle)
     737             : 
     738         105 :       NULLIFY (logger, print_key)
     739         105 :       logger => cp_get_default_logger()
     740             : 
     741         105 :       IF (logger%para_env%is_source()) THEN
     742             :          print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
     743          55 :                                                  "MOTION%PINT%HELIUM%PRINT%ACCEPTS")
     744             :          should_output = BTEST(cp_print_key_should_output( &
     745             :                                iteration_info=logger%iter_info, &
     746          55 :                                basis_section=print_key), cp_p_file)
     747             : 
     748          55 :          IF (should_output) THEN
     749             :             unit_nr = cp_print_key_unit_nr( &
     750             :                       logger, &
     751             :                       print_key, &
     752             :                       middle_name="helium-accepts", &
     753             :                       extension=".dat", &
     754           0 :                       is_new_file=file_is_new)
     755             : 
     756           0 :             IF (file_is_new) THEN
     757             :                WRITE (unit_nr, '(A8,1X,A15,1X,A20)', ADVANCE='NO') &
     758           0 :                   "# Length", &
     759           0 :                   "         Trials", &
     760           0 :                   "            Selected"
     761           0 :                DO j = 1, helium_env(1)%helium%bisctlog2
     762           0 :                   WRITE (unit_nr, '(A17,1X,I3)', ADVANCE='NO') "            Level", j
     763             :                END DO
     764           0 :                WRITE (unit_nr, '(A)') ""
     765             :             END IF
     766             : 
     767           0 :             DO i = 1, helium_env(1)%helium%maxcycle
     768           0 :                WRITE (unit_nr, '(I3)', ADVANCE='NO') i
     769           0 :                DO j = 1, helium_env(1)%helium%bisctlog2 + 2
     770           0 :                   WRITE (unit_nr, '(1X,F20.2)', ADVANCE='NO') helium_env(1)%helium%num_accepted(j, i)
     771             :                END DO
     772           0 :                WRITE (unit_nr, '(A)') ""
     773             :             END DO
     774           0 :             WRITE (unit_nr, '(A)') "&"
     775             : 
     776           0 :             CALL m_flush(unit_nr)
     777           0 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     778             : 
     779             :          END IF
     780             :       END IF
     781             : 
     782         105 :       CALL timestop(handle)
     783         105 :       RETURN
     784             :    END SUBROUTINE helium_print_accepts
     785             : 
     786             : ! ***************************************************************************
     787             : !> \brief  Print permutation state according to HELIUM%PRINT%PERM
     788             : !> \param helium_env ...
     789             : !> \date   2010-06-07
     790             : !> \par    History
     791             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     792             : !> \author Lukasz Walewski
     793             : ! **************************************************************************************************
     794         105 :    SUBROUTINE helium_print_perm(helium_env)
     795             : 
     796             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     797             : 
     798             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_print_perm'
     799             : 
     800             :       CHARACTER                                          :: left_delim, right_delim
     801             :       CHARACTER(len=default_string_length)               :: msg_str, my_middle_name, stmp
     802             :       INTEGER                                            :: curat, handle, i, irank, j, lb, msglen, &
     803             :                                                             nused, offset, outformat, ub, unit_nr
     804         105 :       INTEGER, DIMENSION(:), POINTER                     :: my_cycle, my_perm, used_indices
     805             :       LOGICAL                                            :: first, first_cycle, should_output
     806             :       TYPE(cp_logger_type), POINTER                      :: logger
     807             :       TYPE(section_vals_type), POINTER                   :: print_key
     808             : 
     809         105 :       CALL timeset(routineN, handle)
     810             : 
     811         105 :       NULLIFY (logger, print_key)
     812         105 :       NULLIFY (used_indices)
     813         105 :       logger => cp_get_default_logger()
     814             : 
     815             :       ! determine offset for arrays
     816         105 :       offset = 0
     817         155 :       DO i = 1, logger%para_env%mepos
     818         155 :          offset = offset + helium_env(1)%env_all(i)
     819             :       END DO
     820             : 
     821             :       print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
     822         105 :                                               "MOTION%PINT%HELIUM%PRINT%PERM")
     823             :       should_output = BTEST(cp_print_key_should_output( &
     824             :                             iteration_info=logger%iter_info, &
     825         105 :                             basis_section=print_key), cp_p_file)
     826             : 
     827         105 :       IF (.NOT. should_output) THEN
     828         105 :          CALL timestop(handle)
     829         105 :          RETURN
     830             :       END IF
     831             : 
     832             :       ! get the output type
     833           0 :       CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
     834           0 :       IF (outformat .EQ. perm_cycle) THEN
     835             :          ! gather positions from all processors
     836           0 :          helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
     837           0 :          j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
     838           0 :          DO i = 1, SIZE(helium_env)
     839             :             helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = &
     840           0 :                PACK(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .TRUE.)
     841             :          END DO
     842           0 :          CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d)
     843             :          ! set logical mask for unpacking coordinates gathered from other ranks
     844           0 :          helium_env(1)%helium%ltmp_3_atoms_beads_3d(:, :, :) = .TRUE.
     845             :       END IF
     846             : 
     847             :       ! gather permutation state from all processors to logger%para_env%source
     848           0 :       helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
     849           0 :       msglen = SIZE(helium_env(1)%helium%permutation)
     850           0 :       DO i = 1, SIZE(helium_env)
     851           0 :          helium_env(1)%helium%itmp_atoms_np_1d(msglen*(offset + i - 1) + 1:msglen*(offset + i)) = helium_env(i)%helium%permutation
     852             :       END DO
     853             : 
     854           0 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%itmp_atoms_np_1d)
     855             : 
     856           0 :       IF (logger%para_env%is_source()) THEN
     857             : 
     858             :          ! iterate over helium environments
     859           0 :          DO irank = 1, helium_env(1)%helium%num_env
     860             : 
     861             :             ! generate one file per environment
     862           0 :             stmp = ""
     863           0 :             WRITE (stmp, *) irank
     864           0 :             my_middle_name = "helium-perm-"//TRIM(ADJUSTL(stmp))
     865             :             unit_nr = cp_print_key_unit_nr( &
     866             :                       logger, &
     867             :                       print_key, &
     868             :                       middle_name=TRIM(my_middle_name), &
     869           0 :                       extension=".dat")
     870             : 
     871             :             ! unpack permutation state (actually point to the right section only)
     872           0 :             lb = (irank - 1)*helium_env(1)%helium%atoms + 1
     873           0 :             ub = irank*helium_env(1)%helium%atoms
     874           0 :             my_perm => helium_env(1)%helium%itmp_atoms_np_1d(lb:ub)
     875             : 
     876           0 :             SELECT CASE (outformat)
     877             : 
     878             :             CASE (perm_cycle)
     879             :                ! write the permutation grouped into cycles
     880             : 
     881             :                ! unpack coordinates (necessary only for winding path delimiters)
     882           0 :                msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
     883           0 :                offset = (irank - 1)*msglen
     884             :                helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
     885             :                   UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
     886           0 :                          MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
     887             : 
     888           0 :                curat = 1
     889           0 :                nused = 0
     890           0 :                first_cycle = .TRUE.
     891           0 :                DO WHILE (curat .LE. helium_env(1)%helium%atoms)
     892             : 
     893             :                   ! get the permutation cycle the current atom belongs to
     894           0 :                   my_cycle => helium_cycle_of(curat, my_perm)
     895             : 
     896             :                   ! include the current cycle in the pool of "used" indices
     897           0 :                   nused = nused + SIZE(my_cycle)
     898           0 :                   CALL reallocate(used_indices, 1, nused)
     899           0 :                   used_indices(nused - SIZE(my_cycle) + 1:nused) = my_cycle(1:SIZE(my_cycle))
     900             : 
     901             :                   ! select delimiters according to the cycle's winding state
     902           0 :                   IF (helium_is_winding(helium_env(1)%helium, curat, helium_env(1)%helium%work, my_perm)) THEN
     903           0 :                      left_delim = "["
     904           0 :                      right_delim = "]"
     905             :                   ELSE
     906           0 :                      left_delim = "("
     907           0 :                      right_delim = ")"
     908             :                   END IF
     909             : 
     910             :                   ! cycle delimiter
     911           0 :                   IF (first_cycle) THEN
     912             :                      first_cycle = .FALSE.
     913             :                   ELSE
     914           0 :                      WRITE (unit_nr, '(1X)', ADVANCE='NO')
     915             :                   END IF
     916             : 
     917             :                   ! write out the current cycle
     918           0 :                   WRITE (unit_nr, '(A1)', ADVANCE='NO') left_delim
     919           0 :                   first = .TRUE.
     920           0 :                   DO i = 1, SIZE(my_cycle)
     921           0 :                      IF (first) THEN
     922             :                         first = .FALSE.
     923             :                      ELSE
     924           0 :                         WRITE (unit_nr, '(1X)', ADVANCE='NO')
     925             :                      END IF
     926           0 :                      WRITE (unit_nr, '(I0)', ADVANCE='NO') my_cycle(i)
     927             :                   END DO
     928           0 :                   WRITE (unit_nr, '(A1)', ADVANCE='NO') right_delim
     929             : 
     930             :                   ! clean up
     931           0 :                   DEALLOCATE (my_cycle)
     932           0 :                   NULLIFY (my_cycle)
     933             : 
     934             :                   ! try to increment the current atom index
     935           0 :                   DO WHILE (ANY(used_indices .EQ. curat))
     936           0 :                      curat = curat + 1
     937             :                   END DO
     938             : 
     939             :                END DO
     940           0 :                WRITE (unit_nr, *)
     941             : 
     942           0 :                DEALLOCATE (used_indices)
     943           0 :                NULLIFY (used_indices)
     944             : 
     945             :             CASE (perm_plain)
     946             :                ! write the plain permutation
     947             : 
     948           0 :                first = .TRUE.
     949           0 :                DO i = 1, helium_env(1)%helium%atoms
     950           0 :                   IF (first) THEN
     951             :                      first = .FALSE.
     952             :                   ELSE
     953           0 :                      WRITE (unit_nr, '(1X)', ADVANCE='NO')
     954             :                   END IF
     955           0 :                   WRITE (unit_nr, '(I0)', ADVANCE='NO') my_perm(i)
     956             :                END DO
     957           0 :                WRITE (unit_nr, '(A)') ""
     958             : 
     959             :             CASE default
     960             : 
     961           0 :                msg_str = "Format for permutation output unknown."
     962           0 :                CPABORT(msg_str)
     963             : 
     964             :             END SELECT
     965             : 
     966           0 :             CALL m_flush(unit_nr)
     967           0 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     968             : 
     969             :          END DO
     970             : 
     971             :       END IF
     972             : 
     973           0 :       CALL timestop(handle)
     974             : 
     975           0 :       RETURN
     976         105 :    END SUBROUTINE helium_print_perm
     977             : 
     978             : ! **************************************************************************************************
     979             : !> \brief Print helium action file to HELIUM%PRINT%ACTION
     980             : !> \param pint_env ...
     981             : !> \param helium_env ...
     982             : !> \date   2016-06-07
     983             : !> \par    History
     984             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
     985             : !> \author Felix Uhl
     986             : ! **************************************************************************************************
     987         105 :    SUBROUTINE helium_print_action(pint_env, helium_env)
     988             : 
     989             :       TYPE(pint_env_type), INTENT(IN)                    :: pint_env
     990             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
     991             : 
     992             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_action'
     993             : 
     994             :       CHARACTER(len=default_string_length)               :: my_middle_name, stmp
     995             :       INTEGER                                            :: handle, i, irank, iteration, k, offset, &
     996             :                                                             unit_nr
     997             :       LOGICAL                                            :: file_is_new, should_output
     998         105 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: tmp_inter_action, tmp_link_action, &
     999         105 :                                                             tmp_pair_action
    1000             :       TYPE(cp_logger_type), POINTER                      :: logger
    1001             :       TYPE(section_vals_type), POINTER                   :: print_key
    1002             : 
    1003         105 :       CALL timeset(routineN, handle)
    1004             : 
    1005         105 :       NULLIFY (logger, print_key)
    1006         105 :       logger => cp_get_default_logger()
    1007             : 
    1008             :       ! determine offset for arrays
    1009         105 :       offset = 0
    1010         155 :       DO i = 1, logger%para_env%mepos
    1011         155 :          offset = offset + helium_env(1)%env_all(i)
    1012             :       END DO
    1013             : 
    1014             :       print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1015         105 :                                               "MOTION%PINT%HELIUM%PRINT%ACTION")
    1016             :       should_output = BTEST(cp_print_key_should_output( &
    1017             :                             iteration_info=logger%iter_info, &
    1018         105 :                             basis_section=print_key), cp_p_file)
    1019             : 
    1020         105 :       IF (.NOT. should_output) THEN
    1021          85 :          CALL timestop(handle)
    1022          85 :          RETURN
    1023             :       END IF
    1024             : 
    1025          70 :       DO k = 1, SIZE(helium_env)
    1026          50 :          helium_env(k)%helium%link_action = helium_total_link_action(helium_env(k)%helium)
    1027          50 :          helium_env(k)%helium%pair_action = helium_total_pair_action(helium_env(k)%helium)
    1028          70 :          helium_env(k)%helium%inter_action = helium_total_inter_action(pint_env, helium_env(k)%helium)
    1029             :       END DO
    1030             : 
    1031          20 :       NULLIFY (tmp_link_action)
    1032          20 :       NULLIFY (tmp_pair_action)
    1033          20 :       NULLIFY (tmp_inter_action)
    1034          60 :       ALLOCATE (tmp_link_action(helium_env(1)%helium%num_env))
    1035          40 :       ALLOCATE (tmp_pair_action(helium_env(1)%helium%num_env))
    1036          40 :       ALLOCATE (tmp_inter_action(helium_env(1)%helium%num_env))
    1037         120 :       tmp_link_action(:) = 0.0_dp
    1038         120 :       tmp_pair_action(:) = 0.0_dp
    1039         120 :       tmp_inter_action(:) = 0.0_dp
    1040             :       ! gather Action from all processors to logger%para_env%source
    1041          70 :       DO k = 1, SIZE(helium_env)
    1042          50 :          tmp_link_action(offset + k) = helium_env(k)%helium%link_action
    1043          50 :          tmp_pair_action(offset + k) = helium_env(k)%helium%pair_action
    1044          70 :          tmp_inter_action(offset + k) = helium_env(k)%helium%inter_action
    1045             :       END DO
    1046         220 :       CALL helium_env(1)%comm%sum(tmp_link_action)
    1047         220 :       CALL helium_env(1)%comm%sum(tmp_pair_action)
    1048         220 :       CALL helium_env(1)%comm%sum(tmp_inter_action)
    1049             : 
    1050          20 :       IF (logger%para_env%is_source()) THEN
    1051          10 :          iteration = logger%iter_info%iteration(2)
    1052             :          ! iterate over processors/helium environments
    1053          60 :          DO irank = 1, helium_env(1)%helium%num_env
    1054             : 
    1055             :             ! generate one file per helium_env
    1056          50 :             stmp = ""
    1057          50 :             WRITE (stmp, *) irank
    1058          50 :             my_middle_name = "helium-action-"//TRIM(ADJUSTL(stmp))
    1059             :             unit_nr = cp_print_key_unit_nr( &
    1060             :                       logger, &
    1061             :                       print_key, &
    1062             :                       middle_name=TRIM(my_middle_name), &
    1063             :                       extension=".dat", &
    1064          50 :                       is_new_file=file_is_new)
    1065             : 
    1066          50 :             IF (file_is_new) THEN
    1067             :                WRITE (unit_nr, '(A9,3(1X,A25))') &
    1068           5 :                   "#    Step", &
    1069           5 :                   "     He_Total_Link_Action", &
    1070           5 :                   "     He_Total_Pair_Action", &
    1071          10 :                   "     He_Total_Interaction"
    1072             :             END IF
    1073             : 
    1074             :             WRITE (unit_nr, "(I9,3(1X,F25.14))") &
    1075          50 :                iteration, &
    1076          50 :                tmp_link_action(irank), &
    1077          50 :                tmp_pair_action(irank), &
    1078         100 :                tmp_inter_action(irank)
    1079             : 
    1080          50 :             CALL m_flush(unit_nr)
    1081          60 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1082             : 
    1083             :          END DO
    1084             :       END IF
    1085             : 
    1086          20 :       DEALLOCATE (tmp_link_action)
    1087          20 :       DEALLOCATE (tmp_pair_action)
    1088          20 :       DEALLOCATE (tmp_inter_action)
    1089          20 :       NULLIFY (tmp_link_action)
    1090          20 :       NULLIFY (tmp_pair_action)
    1091          20 :       NULLIFY (tmp_inter_action)
    1092             : 
    1093          20 :       CALL timestop(handle)
    1094             : 
    1095          20 :       RETURN
    1096         105 :    END SUBROUTINE helium_print_action
    1097             : 
    1098             : ! ***************************************************************************
    1099             : !> \brief  Print coordinates according to HELIUM%PRINT%COORDINATES
    1100             : !> \param helium_env ...
    1101             : !> \date   2009-07-16
    1102             : !> \par    History
    1103             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1104             : !> \author Lukasz Walewski
    1105             : ! **************************************************************************************************
    1106         105 :    SUBROUTINE helium_print_coordinates(helium_env)
    1107             : 
    1108             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1109             : 
    1110             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_coordinates'
    1111             : 
    1112             :       CHARACTER(3)                                       :: resName
    1113             :       CHARACTER(len=default_string_length)               :: fmt_string, my_middle_name, stmp
    1114             :       INTEGER                                            :: handle, i, ia, ib, ib1, ib2, ic, icycle, &
    1115             :                                                             irank, j, k, msglen, offset, &
    1116             :                                                             outformat, tmp1, tmp2, unit_nr
    1117         105 :       INTEGER, DIMENSION(:), POINTER                     :: my_perm
    1118             :       LOGICAL                                            :: are_connected, is_winding, ltmp, &
    1119             :                                                             should_output
    1120             :       REAL(KIND=dp)                                      :: xtmp, ytmp, ztmp
    1121             :       REAL(KIND=dp), DIMENSION(3)                        :: r0, r1, r2
    1122             :       TYPE(cp_logger_type), POINTER                      :: logger
    1123             :       TYPE(section_vals_type), POINTER                   :: print_key
    1124             : 
    1125         105 :       CALL timeset(routineN, handle)
    1126             : 
    1127         105 :       NULLIFY (logger, print_key)
    1128         105 :       logger => cp_get_default_logger()
    1129             : 
    1130             :       ! determine offset for arrays
    1131         105 :       offset = 0
    1132         155 :       DO i = 1, logger%para_env%mepos
    1133         155 :          offset = offset + helium_env(1)%env_all(i)
    1134             :       END DO
    1135             : 
    1136             :       print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1137         105 :                                               "MOTION%PINT%HELIUM%PRINT%COORDINATES")
    1138             :       should_output = BTEST(cp_print_key_should_output( &
    1139             :                             iteration_info=logger%iter_info, &
    1140         105 :                             basis_section=print_key), cp_p_file)
    1141             : 
    1142         105 :       IF (.NOT. should_output) THEN
    1143          85 :          CALL timestop(handle)
    1144          85 :          RETURN
    1145             :       END IF
    1146             : 
    1147             :       ! prepare the coordinates for output (use unit cell centered around r0)
    1148          70 :       DO k = 1, SIZE(helium_env)
    1149         200 :          r0(:) = helium_env(k)%helium%center(:)
    1150         320 :          DO ia = 1, helium_env(k)%helium%atoms
    1151        4300 :             DO ib = 1, helium_env(k)%helium%beads
    1152       16000 :                r1(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
    1153       16000 :                r2(:) = helium_env(k)%helium%pos(:, ia, ib) - r0(:)
    1154        4000 :                CALL helium_pbc(helium_env(k)%helium, r2)
    1155        4000 :                ltmp = .FALSE.
    1156       16000 :                DO ic = 1, 3
    1157       16000 :                   IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
    1158           0 :                      ltmp = .TRUE.
    1159           0 :                      CYCLE
    1160             :                   END IF
    1161             :                END DO
    1162        4250 :                IF (ltmp) THEN
    1163           0 :                   helium_env(k)%helium%work(:, ia, ib) = r0(:) + r2(:)
    1164             :                ELSE
    1165       16000 :                   helium_env(k)%helium%work(:, ia, ib) = helium_env(k)%helium%pos(:, ia, ib)
    1166             :                END IF
    1167             :             END DO
    1168             :          END DO
    1169             :       END DO
    1170             : 
    1171             :       ! gather positions from all processors to logger%para_env%source
    1172       24020 :       helium_env(1)%helium%rtmp_3_atoms_beads_np_1d = 0.0_dp
    1173          20 :       j = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
    1174          70 :       DO i = 1, SIZE(helium_env)
    1175             :          helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = &
    1176       12070 :             PACK(helium_env(i)%helium%pos(:, :, 1:helium_env(1)%helium%beads), .TRUE.)
    1177             :       END DO
    1178       48020 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d)
    1179             : 
    1180             :       ! gather permutation state from all processors to logger%para_env%source
    1181         520 :       helium_env(1)%helium%itmp_atoms_np_1d(:) = 0
    1182          20 :       j = SIZE(helium_env(1)%helium%permutation)
    1183          70 :       DO i = 1, SIZE(helium_env)
    1184         320 :          helium_env(1)%helium%itmp_atoms_np_1d(j*(offset + i - 1) + 1:j*(offset + i)) = helium_env(i)%helium%permutation
    1185             :       END DO
    1186             : 
    1187        1020 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%itmp_atoms_np_1d)
    1188             : 
    1189             :       ! set logical mask for unpacking coordinates gathered from other ranks
    1190        6740 :       helium_env(1)%helium%ltmp_3_atoms_beads_3d(:, :, :) = .TRUE.
    1191             : 
    1192          20 :       IF (logger%para_env%is_source()) THEN
    1193             : 
    1194          10 :          CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
    1195             : 
    1196             :          ! iterate over helium environments
    1197          60 :          DO irank = 1, helium_env(1)%helium%num_env
    1198          60 :             IF (outformat == fmt_id_pdb) THEN
    1199             :                ! generate one file per environment
    1200          50 :                stmp = ""
    1201          50 :                WRITE (stmp, *) irank
    1202          50 :                my_middle_name = "helium-pos-"//TRIM(ADJUSTL(stmp))
    1203             :                unit_nr = cp_print_key_unit_nr( &
    1204             :                          logger, &
    1205             :                          print_key, &
    1206             :                          middle_name=TRIM(my_middle_name), &
    1207          50 :                          extension=".pdb")
    1208             : 
    1209             :                ! write out the unit cell parameters
    1210          50 :                fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
    1211          50 :                xtmp = helium_env(1)%helium%cell_size
    1212          50 :                xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
    1213         100 :                SELECT CASE (helium_env(1)%helium%cell_shape)
    1214             :                CASE (helium_cell_shape_cube)
    1215          50 :                   stmp = "C          "
    1216             :                CASE (helium_cell_shape_octahedron)
    1217           0 :                   stmp = "O          "
    1218             :                CASE DEFAULT
    1219          50 :                   stmp = "UNKNOWN    "
    1220             :                END SELECT
    1221          50 :                WRITE (unit_nr, fmt_string) "CRYST1", &
    1222          50 :                   xtmp, xtmp, xtmp, &
    1223          50 :                   90.0_dp, 90.0_dp, 90.0_dp, &
    1224         100 :                   stmp, helium_env(1)%helium%beads
    1225             : 
    1226             :                ! unpack coordinates
    1227          50 :                msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
    1228          50 :                offset = (irank - 1)*msglen
    1229             :                helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
    1230             :                   UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
    1231       16850 :                          MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
    1232             : 
    1233             :                ! unpack permutation state (actually point to the right section only)
    1234          50 :                msglen = SIZE(helium_env(1)%helium%permutation)
    1235          50 :                offset = (irank - 1)*msglen
    1236          50 :                my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
    1237             : 
    1238             :                ! write out coordinates
    1239             :                fmt_string = &
    1240          50 :                   "(A6,I5,1X,A4,A1,A3,1X,A1,I4,A1,3X,3F8.3,2F6.2,10X,A2,A2)"
    1241         300 :                DO ia = 1, helium_env(1)%helium%atoms
    1242         250 :                   icycle = helium_cycle_number(helium_env(1)%helium, ia, my_perm)
    1243         250 :                   is_winding = helium_is_winding(helium_env(1)%helium, ia, helium_env(1)%helium%work, my_perm)
    1244         250 :                   IF (is_winding) THEN
    1245           0 :                      resName = "SPR"
    1246             :                   ELSE
    1247         250 :                      resName = "NRM"
    1248             :                   END IF
    1249        4300 :                   DO ib = 1, helium_env(1)%helium%beads
    1250        4000 :                      xtmp = helium_env(1)%helium%work(1, ia, ib)
    1251        4000 :                      xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
    1252        4000 :                      ytmp = helium_env(1)%helium%work(2, ia, ib)
    1253        4000 :                      ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
    1254        4000 :                      ztmp = helium_env(1)%helium%work(3, ia, ib)
    1255        4000 :                      ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
    1256        4000 :                      WRITE (unit_nr, fmt_string) "ATOM  ", &
    1257        4000 :                         (ia - 1)*helium_env(1)%helium%beads + ib, &
    1258        4000 :                         " He ", " ", resName, "X", &
    1259        4000 :                         icycle, &
    1260        4000 :                         " ", &
    1261        4000 :                         xtmp, ytmp, ztmp, &
    1262        8250 :                         1.0_dp, 0.0_dp, "HE", "  "
    1263             :                   END DO
    1264             :                END DO
    1265             : 
    1266             :                ! write out the bead connectivity information
    1267         300 :                DO ia = 1, helium_env(1)%helium%atoms
    1268             : 
    1269             :                   ! write connectivity records for this atom only if the path
    1270             :                   ! it belongs to is longer than 1.
    1271         250 :                   IF (helium_path_length(helium_env(1)%helium, ia, my_perm) .LE. 1) THEN
    1272             :                      CYCLE
    1273             :                   END IF
    1274             : 
    1275           0 :                   DO ib = 1, helium_env(1)%helium%beads - 1
    1276             :                      ! check whether the consecutive beads belong to the same box
    1277           0 :                      r1(:) = helium_env(1)%helium%work(:, ia, ib) - helium_env(1)%helium%work(:, ia, ib + 1)
    1278           0 :                      r2(:) = r1(:)
    1279           0 :                      CALL helium_pbc(helium_env(1)%helium, r2)
    1280           0 :                      are_connected = .TRUE.
    1281           0 :                      DO ic = 1, 3
    1282           0 :                         IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
    1283             :                            ! if the distance betw ib and ib+1 changes upon applying
    1284             :                            ! PBC do not connect them
    1285           0 :                            are_connected = .FALSE.
    1286           0 :                            CYCLE
    1287             :                         END IF
    1288             :                      END DO
    1289           0 :                      IF (are_connected) THEN
    1290           0 :                         tmp1 = (ia - 1)*helium_env(1)%helium%beads + ib
    1291           0 :                         tmp2 = (ia - 1)*helium_env(1)%helium%beads + ib + 1
    1292             :                         ! smaller value has to go first
    1293             :                         IF (tmp1 .LT. tmp2) THEN
    1294           0 :                            ib1 = tmp1
    1295           0 :                            ib2 = tmp2
    1296             :                         ELSE
    1297             :                            ib1 = tmp2
    1298             :                            ib2 = tmp1
    1299             :                         END IF
    1300           0 :                         WRITE (unit_nr, '(A6,2I5)') "CONECT", ib1, ib2
    1301             :                      END IF
    1302             :                   END DO
    1303             : 
    1304             :                   ! last bead of atom <ia> connects to the first bead
    1305             :                   ! of the next atom in the permutation cycle
    1306           0 :                  r1(:) = helium_env(1)%helium%work(:, ia, helium_env(1)%helium%beads) - helium_env(1)%helium%work(:, my_perm(ia), 1)
    1307           0 :                   r2(:) = r1(:)
    1308           0 :                   CALL helium_pbc(helium_env(1)%helium, r2)
    1309           0 :                   are_connected = .TRUE.
    1310           0 :                   DO ic = 1, 3
    1311           0 :                      IF (ABS(r1(ic) - r2(ic)) .GT. 100.0_dp*EPSILON(0.0_dp)) THEN
    1312             :                         ! if the distance betw ib and ib+1 changes upon applying
    1313             :                         ! PBC do not connect them
    1314           0 :                         are_connected = .FALSE.
    1315           0 :                         CYCLE
    1316             :                      END IF
    1317             :                   END DO
    1318          50 :                   IF (are_connected) THEN
    1319           0 :                      tmp1 = ia*helium_env(1)%helium%beads
    1320           0 :                      tmp2 = (my_perm(ia) - 1)*helium_env(1)%helium%beads + 1
    1321           0 :                      IF (tmp1 .LT. tmp2) THEN
    1322           0 :                         ib1 = tmp1
    1323           0 :                         ib2 = tmp2
    1324             :                      ELSE
    1325           0 :                         ib1 = tmp2
    1326           0 :                         ib2 = tmp1
    1327             :                      END IF
    1328           0 :                      WRITE (unit_nr, '(A6,2I5)') "CONECT", ib1, ib2
    1329             :                   END IF
    1330             :                END DO
    1331          50 :                WRITE (unit_nr, '(A)') "END"
    1332             : 
    1333          50 :                CALL m_flush(unit_nr)
    1334          50 :                CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1335           0 :             ELSE IF (outformat == fmt_id_xyz) THEN
    1336             :                ! generate one file per environment and bead
    1337           0 :                DO ib = 1, helium_env(1)%helium%beads
    1338           0 :                   stmp = ""
    1339           0 :                   WRITE (stmp, *) irank
    1340           0 :                   my_middle_name = "helium-pos-"//TRIM(ADJUSTL(stmp))
    1341           0 :                   WRITE (stmp, *) ib
    1342           0 :                   my_middle_name = TRIM(my_middle_name)//"-"//TRIM(ADJUSTL(stmp))
    1343             :                   unit_nr = cp_print_key_unit_nr( &
    1344             :                             logger, &
    1345             :                             print_key, &
    1346             :                             middle_name=TRIM(my_middle_name), &
    1347           0 :                             extension=".xyz")
    1348             :                   ! write out xyz header
    1349           0 :                   WRITE (unit_nr, *) helium_env(1)%helium%atoms
    1350           0 :                   stmp = ""
    1351           0 :                   WRITE (stmp, *) logger%iter_info%n_rlevel
    1352           0 :                   fmt_string = "(A6,"//TRIM(ADJUSTL(stmp))//"I12)"
    1353           0 :                   WRITE (unit_nr, fmt_string) "iter= ", logger%iter_info%iteration(:)
    1354           0 :                   fmt_string = "(A6,3F9.3,3F7.2,1X,A11,1X,I3)"
    1355             : 
    1356             :                   ! unpack coordinates
    1357           0 :                   msglen = SIZE(helium_env(1)%helium%rtmp_3_atoms_beads_1d)
    1358           0 :                   offset = (irank - 1)*msglen
    1359             :                   helium_env(1)%helium%work(:, :, 1:helium_env(1)%helium%beads) = &
    1360             :                      UNPACK(helium_env(1)%helium%rtmp_3_atoms_beads_np_1d(offset + 1:offset + msglen), &
    1361           0 :                             MASK=helium_env(1)%helium%ltmp_3_atoms_beads_3d, FIELD=0.0_dp)
    1362             : 
    1363             :                   ! unpack permutation state (actually point to the right section only)
    1364           0 :                   msglen = SIZE(helium_env(1)%helium%permutation)
    1365           0 :                   offset = (irank - 1)*msglen
    1366           0 :                   my_perm => helium_env(1)%helium%itmp_atoms_np_1d(offset + 1:offset + msglen)
    1367             : 
    1368             :                   ! write out coordinates
    1369           0 :                   fmt_string = "(A2,3(1X,F20.10))"
    1370           0 :                   DO ia = 1, helium_env(1)%helium%atoms
    1371           0 :                      xtmp = helium_env(1)%helium%work(1, ia, ib)
    1372           0 :                      xtmp = cp_unit_from_cp2k(xtmp, "angstrom")
    1373           0 :                      ytmp = helium_env(1)%helium%work(2, ia, ib)
    1374           0 :                      ytmp = cp_unit_from_cp2k(ytmp, "angstrom")
    1375           0 :                      ztmp = helium_env(1)%helium%work(3, ia, ib)
    1376           0 :                      ztmp = cp_unit_from_cp2k(ztmp, "angstrom")
    1377           0 :                      WRITE (unit_nr, fmt_string) "He", xtmp, ytmp, ztmp
    1378             :                   END DO
    1379           0 :                   CALL m_flush(unit_nr)
    1380           0 :                   CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1381             :                END DO
    1382             :             ELSE
    1383           0 :                CPABORT("")
    1384             :             END IF
    1385             :          END DO
    1386             : 
    1387             :       END IF
    1388             : 
    1389          20 :       CALL timestop(handle)
    1390             : 
    1391          20 :       RETURN
    1392         105 :    END SUBROUTINE helium_print_coordinates
    1393             : 
    1394             : ! ***************************************************************************
    1395             : !> \brief  Print radial distribution functions according to HELIUM%PRINT%RDF
    1396             : !> \param helium_env ...
    1397             : !> \date   2009-07-23
    1398             : !> \par    History
    1399             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1400             : !> \author Lukasz Walewski
    1401             : ! **************************************************************************************************
    1402          10 :    SUBROUTINE helium_print_rdf(helium_env)
    1403             : 
    1404             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1405             : 
    1406             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_print_rdf'
    1407             : 
    1408             :       CHARACTER(len=default_string_length)               :: stmp
    1409             :       INTEGER                                            :: handle, ia, ic, id, itmp, iweight, k, &
    1410             :                                                             nsteps, unit_nr
    1411             :       LOGICAL                                            :: should_output
    1412             :       REAL(KIND=dp)                                      :: inv_norm, rtmp, rtmp2
    1413          10 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: message
    1414             :       TYPE(cp_logger_type), POINTER                      :: logger
    1415             :       TYPE(section_vals_type), POINTER                   :: print_key
    1416             : 
    1417          10 :       CALL timeset(routineN, handle)
    1418             : 
    1419          10 :       NULLIFY (logger, print_key)
    1420          10 :       logger => cp_get_default_logger()
    1421             : 
    1422          10 :       IF (logger%para_env%is_source()) THEN
    1423             :          print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1424           5 :                                                  "MOTION%PINT%HELIUM%PRINT%RDF")
    1425             :          should_output = BTEST(cp_print_key_should_output( &
    1426             :                                iteration_info=logger%iter_info, &
    1427           5 :                                basis_section=print_key), cp_p_file)
    1428             :       END IF
    1429          10 :       CALL helium_env(1)%comm%bcast(should_output, logger%para_env%source)
    1430             : 
    1431          10 :       IF (should_output) THEN
    1432             :          ! work on the temporary array so that accumulated data remains intact
    1433             :          ! save accumulated data of different env on same core in first temp
    1434        5010 :          helium_env(1)%helium%rdf_inst(:, :) = 0.0_dp
    1435          20 :          DO k = 1, SIZE(helium_env)
    1436             :             helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :) + &
    1437        5020 :                                                   helium_env(k)%helium%rdf_accu(:, :)
    1438             :          END DO
    1439             : 
    1440             :          ! average over processors
    1441          10 :          NULLIFY (message)
    1442          10 :          message => helium_env(1)%helium%rdf_inst(:, :)
    1443       10010 :          CALL helium_env(1)%comm%sum(message)
    1444          10 :          itmp = helium_env(1)%helium%num_env
    1445          10 :          inv_norm = 1.0_dp/REAL(itmp, dp)
    1446        5010 :          helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm
    1447             : 
    1448             :          ! average over steps performed so far in this run
    1449          10 :          nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
    1450          10 :          inv_norm = 1.0_dp/REAL(nsteps, dp)
    1451        5010 :          helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)*inv_norm
    1452             : 
    1453          10 :          iweight = helium_env(1)%helium%rdf_iweight
    1454             :          ! average over the old and the current density (observe the weights!)
    1455             :          helium_env(1)%helium%rdf_inst(:, :) = nsteps*helium_env(1)%helium%rdf_inst(:, :) + &
    1456        5010 :                                                iweight*helium_env(1)%helium%rdf_rstr(:, :)
    1457        5010 :          helium_env(1)%helium%rdf_inst(:, :) = helium_env(1)%helium%rdf_inst(:, :)/REAL(nsteps + iweight, dp)
    1458             : 
    1459          10 :          IF (logger%para_env%is_source()) THEN
    1460             : 
    1461           5 :             ia = 0
    1462           5 :             rtmp = cp_unit_from_cp2k(1.0_dp, "angstrom")
    1463           5 :             rtmp2 = 1.0_dp
    1464           5 :             IF (.NOT. helium_env(1)%helium%periodic) THEN
    1465             :                ! RDF in non-periodic case has unit 1/bohr**3, convert to Angstrom:
    1466           0 :                rtmp2 = rtmp**(-3)
    1467             :             END IF
    1468             : 
    1469           5 :             IF (helium_env(1)%helium%rdf_he_he) THEN
    1470             :                ! overwrite RDF file each time it is written
    1471           5 :                ia = 1
    1472           5 :                stmp = ""
    1473           5 :                WRITE (stmp, *) "He-He"
    1474             :                unit_nr = cp_print_key_unit_nr( &
    1475             :                          logger, &
    1476             :                          print_key, &
    1477             :                          middle_name="helium-rdf-"//TRIM(ADJUSTL(stmp)), &
    1478             :                          extension=".dat", &
    1479             :                          file_position="REWIND", &
    1480           5 :                          do_backup=.FALSE.)
    1481             : 
    1482        1255 :                DO ic = 1, helium_env(1)%helium%rdf_nbin
    1483        1250 :                   WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
    1484        1250 :                   WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(ia, ic)*rtmp2
    1485        1255 :                   WRITE (unit_nr, *)
    1486             :                END DO
    1487             : 
    1488           5 :                CALL m_flush(unit_nr)
    1489           5 :                CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1490             :             END IF
    1491             : 
    1492           5 :             IF (helium_env(1)%helium%rdf_sol_he) THEN
    1493             :                ! overwrite RDF file each time it is written
    1494           0 :                stmp = ""
    1495           0 :                WRITE (stmp, *) "Solute-He"
    1496             :                unit_nr = cp_print_key_unit_nr( &
    1497             :                          logger, &
    1498             :                          print_key, &
    1499             :                          middle_name="helium-rdf-"//TRIM(ADJUSTL(stmp)), &
    1500             :                          extension=".dat", &
    1501             :                          file_position="REWIND", &
    1502           0 :                          do_backup=.FALSE.)
    1503             : 
    1504           0 :                DO ic = 1, helium_env(1)%helium%rdf_nbin
    1505           0 :                   WRITE (unit_nr, '(F20.10)', ADVANCE='NO') (REAL(ic, dp) - 0.5_dp)*helium_env(1)%helium%rdf_delr*rtmp
    1506           0 :                   DO id = 1 + ia, helium_env(1)%helium%rdf_num
    1507           0 :                      WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%rdf_inst(id, ic)*rtmp2
    1508             :                   END DO
    1509           0 :                   WRITE (unit_nr, *)
    1510             :                END DO
    1511             : 
    1512           0 :                CALL m_flush(unit_nr)
    1513           0 :                CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1514             :             END IF
    1515             : 
    1516             :          END IF
    1517             :       END IF
    1518             : 
    1519          10 :       CALL timestop(handle)
    1520             : 
    1521          10 :       RETURN
    1522          10 :    END SUBROUTINE helium_print_rdf
    1523             : 
    1524             : ! ***************************************************************************
    1525             : !> \brief  Print densities according to HELIUM%PRINT%RHO
    1526             : !> \param helium_env ...
    1527             : !> \date   2011-06-21
    1528             : !> \par History
    1529             : !>      08.2015 cleaned code from unneeded arrays
    1530             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1531             : !> \author Lukasz Walewski
    1532             : ! **************************************************************************************************
    1533          10 :    SUBROUTINE helium_print_rho(helium_env)
    1534             : 
    1535             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1536             : 
    1537             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'helium_print_rho'
    1538             : 
    1539             :       CHARACTER(len=default_string_length)               :: comment, fname
    1540             :       INTEGER                                            :: handle, ic, id, itmp, iweight, k, &
    1541             :                                                             nsteps, unit_nr
    1542             :       LOGICAL                                            :: should_output
    1543             :       REAL(KIND=dp)                                      :: inv_norm, invproc
    1544             :       REAL(KIND=dp), DIMENSION(3)                        :: center
    1545          10 :       REAL(KIND=dp), DIMENSION(:, :, :), POINTER         :: cubdata
    1546             :       TYPE(cp_logger_type), POINTER                      :: logger
    1547             :       TYPE(section_vals_type), POINTER                   :: print_key
    1548             : 
    1549          10 :       CALL timeset(routineN, handle)
    1550             : 
    1551          10 :       NULLIFY (logger, print_key)
    1552          10 :       logger => cp_get_default_logger()
    1553             : 
    1554             :       print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1555          10 :                                               "MOTION%PINT%HELIUM%PRINT%RHO")
    1556             :       should_output = BTEST(cp_print_key_should_output( &
    1557             :                             iteration_info=logger%iter_info, &
    1558          10 :                             basis_section=print_key), cp_p_file)
    1559             : 
    1560          10 :       IF (.NOT. should_output) THEN
    1561           8 :          CALL timestop(handle)
    1562           8 :          RETURN
    1563             :       END IF
    1564             : 
    1565             :       ! work on temporary array so that the average remains intact
    1566             :       ! save accumulated data of different env on same core in first temp
    1567        4222 :       helium_env(1)%helium%rho_inst(:, :, :, :) = 0.0_dp
    1568           4 :       DO k = 1, SIZE(helium_env)
    1569             :          helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :) + &
    1570        4224 :                                                      helium_env(k)%helium%rho_accu(:, :, :, :)
    1571             :       END DO
    1572             : 
    1573             :       ! average over processors
    1574        8442 :       CALL helium_env(1)%comm%sum(helium_env(1)%helium%rho_inst)
    1575           2 :       itmp = helium_env(1)%helium%num_env
    1576           2 :       invproc = 1.0_dp/REAL(itmp, dp)
    1577        4222 :       helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*invproc
    1578             : 
    1579             :       ! average over steps performed so far in this run
    1580           2 :       nsteps = helium_env(1)%helium%current_step - helium_env(1)%helium%first_step
    1581           2 :       inv_norm = 1.0_dp/REAL(nsteps, dp)
    1582        4222 :       helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)*inv_norm
    1583             : 
    1584           2 :       iweight = helium_env(1)%helium%rho_iweight
    1585             :       ! average over the old and the current density (observe the weights!)
    1586             :       helium_env(1)%helium%rho_inst(:, :, :, :) = nsteps*helium_env(1)%helium%rho_inst(:, :, :, :) + &
    1587        4222 :                                                   iweight*helium_env(1)%helium%rho_rstr(:, :, :, :)
    1588        4222 :       helium_env(1)%helium%rho_inst(:, :, :, :) = helium_env(1)%helium%rho_inst(:, :, :, :)/REAL(nsteps + iweight, dp)
    1589             : 
    1590             :       ! set center of the cubefile
    1591           2 :       IF (helium_env(1)%helium%solute_present) THEN
    1592             :          ! should be set to solute's COM
    1593           0 :          center(:) = helium_env(1)%helium%center(:)
    1594             :       ELSE
    1595             :          ! regardless of whether we are periodic or not use the origin, since
    1596             :          ! pure cluster's COM might drift, but we want the cube fixed (note that
    1597             :          ! the densities are correctly calculated wrt to COM in such case)
    1598           2 :          center(:) = (/0.0_dp, 0.0_dp, 0.0_dp/)
    1599             :       END IF
    1600             : 
    1601          12 :       DO id = 1, rho_num ! loop over densities ---
    1602             : 
    1603          10 :          IF (.NOT. helium_env(1)%helium%rho_property(id)%is_calculated) THEN
    1604             :             ! skip densities that are not requested by the user
    1605             :             CYCLE
    1606             :          END IF
    1607             : 
    1608           6 :          DO ic = 1, helium_env(1)%helium%rho_property(id)%num_components ! loop over components
    1609             : 
    1610             :             WRITE (fname, '(A)') "helium-rho-"// &
    1611           2 :                TRIM(ADJUSTL(helium_env(1)%helium%rho_property(id)%filename_suffix(ic)))
    1612           2 :             IF (helium_env(1)%helium%rho_property(id)%component_name(ic) .EQ. "") THEN
    1613           2 :                WRITE (comment, '(A)') TRIM(helium_env(1)%helium%rho_property(id)%name)
    1614             :             ELSE
    1615             :                WRITE (comment, '(A)') TRIM(helium_env(1)%helium%rho_property(id)%name)// &
    1616             :                   ", "// &
    1617           0 :                   TRIM(helium_env(1)%helium%rho_property(id)%component_name(ic))
    1618             :             END IF
    1619           2 :             cubdata => helium_env(1)%helium%rho_inst(helium_env(1)%helium%rho_property(id)%component_index(ic), :, :, :)
    1620             : 
    1621             :             unit_nr = cp_print_key_unit_nr( &
    1622             :                       logger, &
    1623             :                       print_key, &
    1624             :                       middle_name=TRIM(ADJUSTL(fname)), &
    1625             :                       extension=".cube", &
    1626             :                       file_position="REWIND", &
    1627           2 :                       do_backup=.FALSE.)
    1628             : 
    1629          12 :             IF (logger%para_env%is_source()) THEN
    1630             :                CALL helium_write_cubefile( &
    1631             :                   unit_nr, &
    1632             :                   comment, &
    1633             :                   center - 0.5_dp*(helium_env(1)%helium%rho_maxr - helium_env(1)%helium%rho_delr), &
    1634             :                   helium_env(1)%helium%rho_delr, &
    1635             :                   helium_env(1)%helium%rho_nbin, &
    1636           4 :                   cubdata)
    1637           1 :                CALL m_flush(unit_nr)
    1638           1 :                CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1639             :             END IF
    1640             : 
    1641             :          END DO ! loop over components
    1642             : 
    1643             :       END DO ! loop over densities ---
    1644             : 
    1645           2 :       CALL timestop(handle)
    1646             : 
    1647          10 :    END SUBROUTINE helium_print_rho
    1648             : 
    1649             : ! ***************************************************************************
    1650             : !> \brief Write volumetric data to an orthorhombic cubefile
    1651             : !> \param   unit   unit number to which output will be sent
    1652             : !> \param   comment   description of the data stored in the cubefile
    1653             : !> \param   origin   position of the cubefile origin
    1654             : !> \param   deltar   voxel side length
    1655             : !> \param   ndim   number of voxels in each dimension
    1656             : !> \param   DATA   array (ndim x ndim x ndim) with the data for output
    1657             : !> \date 2013-11-25
    1658             : !> \author Lukasz Walewski
    1659             : ! **************************************************************************************************
    1660           1 :    SUBROUTINE helium_write_cubefile(unit, comment, origin, deltar, ndim, DATA)
    1661             : 
    1662             :       INTEGER, INTENT(IN)                                :: unit
    1663             :       CHARACTER(len=default_string_length), INTENT(IN)   :: comment
    1664             :       REAL(KIND=dp), DIMENSION(3), INTENT(IN)            :: origin
    1665             :       REAL(KIND=dp), INTENT(IN)                          :: deltar
    1666             :       INTEGER, INTENT(IN)                                :: ndim
    1667             :       REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN), &
    1668             :          POINTER                                         :: DATA
    1669             : 
    1670             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_write_cubefile'
    1671             : 
    1672             :       INTEGER                                            :: handle, ix, iy, iz, nw
    1673             :       REAL(KIND=dp)                                      :: delr, inva3
    1674             :       REAL(KIND=dp), DIMENSION(3)                        :: orig
    1675             : 
    1676           1 :       CALL timeset(routineN, handle)
    1677             : 
    1678             :       ! convert cubefile data to the proper units of measure
    1679           1 :       delr = angstrom*deltar
    1680           4 :       orig(:) = angstrom*origin(:)
    1681             : 
    1682             :       ! output cube file header
    1683           1 :       WRITE (unit, '(A)') comment
    1684           1 :       WRITE (unit, '(A)') "Volumetric data in cubefile format generated by CP2K"
    1685           1 :       WRITE (unit, '(I5,3(1X,F12.8))') 0, orig(1), orig(2), orig(3)
    1686           1 :       WRITE (unit, '(I5,3(1X,F12.8))') ndim, delr, 0.0_dp, 0.0_dp
    1687           1 :       WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, delr, 0.0_dp
    1688           1 :       WRITE (unit, '(I5,3(1X,F12.8))') ndim, 0.0_dp, 0.0_dp, delr
    1689             : 
    1690             :       ! output cube file data
    1691           1 :       nw = 0
    1692           1 :       inva3 = 1.0_dp/(angstrom*angstrom*angstrom)
    1693          11 :       DO ix = 1, ndim
    1694         111 :          DO iy = 1, ndim
    1695        1110 :             DO iz = 1, ndim
    1696        1000 :                WRITE (unit, '(1X,E13.5)', ADVANCE='NO') inva3*DATA(ix, iy, iz)
    1697        1000 :                nw = nw + 1
    1698        1100 :                IF (MOD(nw, 6) .EQ. 0) THEN
    1699         166 :                   nw = 0
    1700         166 :                   WRITE (unit, *)
    1701             :                END IF
    1702             :             END DO
    1703             :          END DO
    1704             :       END DO
    1705             :       ! some compilers write over the first entry on a line losing all previous
    1706             :       ! values written on that line unless line terminator is written at the end
    1707             :       ! so make sure that the last WRITE statement runs without ADVANCE='NO'
    1708             :       ! (observed for ifort 14.0.1 and 14.0.2 but not for gfortran 4.8.2)
    1709           1 :       IF (nw .NE. 0) THEN
    1710           1 :          WRITE (unit, *)
    1711             :       END IF
    1712             : 
    1713           1 :       CALL timestop(handle)
    1714             : 
    1715           1 :    END SUBROUTINE helium_write_cubefile
    1716             : 
    1717             : ! ***************************************************************************
    1718             : !> \brief  Print permutation length according to HELIUM%PRINT%PLENGTH
    1719             : !> \param helium_env ...
    1720             : !> \date   2010-06-07
    1721             : !> \par    History
    1722             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1723             : !> \author Lukasz Walewski
    1724             : ! **************************************************************************************************
    1725         105 :    SUBROUTINE helium_print_plength(helium_env)
    1726             : 
    1727             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1728             : 
    1729             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_plength'
    1730             : 
    1731             :       INTEGER                                            :: handle, i, unit_nr
    1732             :       LOGICAL                                            :: is_new, should_output
    1733             :       TYPE(cp_logger_type), POINTER                      :: logger
    1734             :       TYPE(section_vals_type), POINTER                   :: print_key
    1735             : 
    1736         105 :       CALL timeset(routineN, handle)
    1737             : 
    1738         105 :       NULLIFY (logger, print_key)
    1739         105 :       logger => cp_get_default_logger()
    1740             : 
    1741         105 :       IF (logger%para_env%is_source()) THEN
    1742             :          print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1743          55 :                                                  "MOTION%PINT%HELIUM%PRINT%PLENGTH")
    1744             :          should_output = BTEST(cp_print_key_should_output( &
    1745             :                                iteration_info=logger%iter_info, &
    1746          55 :                                basis_section=print_key), cp_p_file)
    1747             : 
    1748          55 :          IF (should_output) THEN
    1749             : 
    1750             :             unit_nr = cp_print_key_unit_nr( &
    1751             :                       logger, &
    1752             :                       print_key, &
    1753             :                       middle_name="helium-plength", &
    1754             :                       extension=".dat", &
    1755           0 :                       is_new_file=is_new)
    1756             : 
    1757           0 :             DO i = 1, helium_env(1)%helium%atoms
    1758           0 :                WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%plength_avrg(i)
    1759           0 :                IF (i .LT. helium_env(1)%helium%atoms) THEN
    1760           0 :                   WRITE (unit_nr, '(1X)', ADVANCE='NO')
    1761             :                END IF
    1762             :             END DO
    1763           0 :             WRITE (unit_nr, '(A)') ""
    1764             : 
    1765           0 :             CALL m_flush(unit_nr)
    1766           0 :             CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1767             : 
    1768             :          END IF
    1769             :       END IF
    1770             : 
    1771         105 :       CALL timestop(handle)
    1772             : 
    1773         105 :       RETURN
    1774             :    END SUBROUTINE helium_print_plength
    1775             : 
    1776             : ! ***************************************************************************
    1777             : !> \brief  Print helium force according to HELIUM%PRINT%FORCE
    1778             : !> \param helium_env ...
    1779             : !> \date   2010-01-27
    1780             : !> \par    History
    1781             : !>         2016-07-14 Modified to work with independent helium_env [cschran]
    1782             : !> \author Lukasz Walewski
    1783             : ! **************************************************************************************************
    1784         105 :    SUBROUTINE helium_print_force(helium_env)
    1785             : 
    1786             :       TYPE(helium_solvent_p_type), DIMENSION(:), POINTER :: helium_env
    1787             : 
    1788             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_force'
    1789             : 
    1790             :       CHARACTER(len=default_string_length)               :: msgstr
    1791             :       INTEGER                                            :: handle, ia, ib, ic, idim, unit_nr
    1792             :       LOGICAL                                            :: should_output
    1793             :       TYPE(cp_logger_type), POINTER                      :: logger
    1794             :       TYPE(section_vals_type), POINTER                   :: print_key
    1795             : 
    1796         105 :       CALL timeset(routineN, handle)
    1797             : 
    1798         105 :       NULLIFY (logger, print_key)
    1799         105 :       logger => cp_get_default_logger()
    1800             : 
    1801             :       print_key => section_vals_get_subs_vals(helium_env(1)%helium%input, &
    1802         105 :                                               "MOTION%PINT%HELIUM%PRINT%FORCES")
    1803             :       should_output = BTEST(cp_print_key_should_output( &
    1804             :                             logger%iter_info, &
    1805         105 :                             basis_section=print_key), cp_p_file)
    1806             : 
    1807         105 :       IF (.NOT. should_output) THEN
    1808         105 :          CALL timestop(handle)
    1809         105 :          RETURN
    1810             :       END IF
    1811             : 
    1812             :       ! check if there is anything to be printed out
    1813           0 :       IF (.NOT. helium_env(1)%helium%solute_present) THEN
    1814           0 :          msgstr = "Warning: force printout requested but there is no solute!"
    1815           0 :          CALL helium_write_line(msgstr)
    1816           0 :          CALL timestop(handle)
    1817           0 :          RETURN
    1818             :       END IF
    1819             : 
    1820           0 :       IF (logger%para_env%is_source()) THEN
    1821             : 
    1822             :          unit_nr = cp_print_key_unit_nr( &
    1823             :                    logger, &
    1824             :                    print_key, &
    1825             :                    middle_name="helium-force", &
    1826           0 :                    extension=".dat")
    1827             : 
    1828             :          ! print all force components in one line
    1829           0 :          DO ib = 1, helium_env(1)%helium%solute_beads
    1830           0 :             idim = 0
    1831           0 :             DO ia = 1, helium_env(1)%helium%solute_atoms
    1832           0 :                DO ic = 1, 3
    1833           0 :                   idim = idim + 1
    1834           0 :                   WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium_env(1)%helium%force_avrg(ib, idim)
    1835             :                END DO
    1836             :             END DO
    1837             :          END DO
    1838           0 :          WRITE (unit_nr, *)
    1839             : 
    1840           0 :          CALL m_flush(unit_nr)
    1841           0 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1842             : 
    1843             :       END IF
    1844             : 
    1845           0 :       CALL timestop(handle)
    1846             : 
    1847           0 :       RETURN
    1848             :    END SUBROUTINE helium_print_force
    1849             : 
    1850             : #if 0
    1851             : 
    1852             : ! ***************************************************************************
    1853             : !> \brief  Print instantaneous force according to HELIUM%PRINT%FORCES_INST
    1854             : !> \param helium ...
    1855             : !> \date   2010-01-29
    1856             : !> \author Lukasz Walewski
    1857             : !> \note   Collects instantaneous helium forces from all processors on
    1858             : !> logger%para_env%source and writes them to files - one file per processor.
    1859             : !>         !TODO: Generalize for helium_env
    1860             : ! **************************************************************************************************
    1861             :    SUBROUTINE helium_print_force_inst(helium)
    1862             : 
    1863             :       TYPE(helium_solvent_type), POINTER                 :: helium
    1864             : 
    1865             :       CHARACTER(len=*), PARAMETER :: routineN = 'helium_print_force_inst', &
    1866             :          routineP = moduleN//':'//routineN
    1867             : 
    1868             :       CHARACTER(len=default_string_length)               :: my_middle_name, stmp
    1869             :       INTEGER                                            :: handle, ia, ib, ic, idim, irank, offset, &
    1870             :                                                             unit_nr
    1871             :       LOGICAL                                            :: should_output
    1872             :       TYPE(cp_logger_type), POINTER                      :: logger
    1873             :       TYPE(section_vals_type), POINTER                   :: print_key
    1874             : 
    1875             :       CALL timeset(routineN, handle)
    1876             : 
    1877             :       NULLIFY (logger, print_key)
    1878             :       logger => cp_get_default_logger()
    1879             :       print_key => section_vals_get_subs_vals(helium%input, &
    1880             :                                               "MOTION%PINT%HELIUM%PRINT%FORCES_INST")
    1881             :       should_output = BTEST(cp_print_key_should_output( &
    1882             :                             logger%iter_info, &
    1883             :                             basis_section=print_key), cp_p_file)
    1884             : 
    1885             :       IF (should_output) THEN
    1886             : 
    1887             :          ! check if there is anything to be printed out
    1888             :          IF (.NOT. helium%solute_present) THEN
    1889             :             stmp = "Warning: force printout requested but there is no solute!"
    1890             :             CALL helium_write_line(stmp)
    1891             :             CALL timestop(handle)
    1892             :             RETURN
    1893             :          END IF
    1894             : 
    1895             :          ! fill the tmp buffer with instantaneous helium forces at each proc
    1896             :          helium%rtmp_p_ndim_1d(:) = PACK(helium%force_inst, .TRUE.)
    1897             : 
    1898             :          ! pass the message from all processors to logger%para_env%source
    1899             :          helium%rtmp_p_ndim_np_1d(:) = 0.0_dp
    1900             :          CALL logger%para_env%gather(helium%rtmp_p_ndim_1d, helium%rtmp_p_ndim_np_1d)
    1901             : 
    1902             :          IF (logger%para_env%is_source()) THEN
    1903             : 
    1904             :             ! iterate over processors/helium environments
    1905             :             DO irank = 1, helium%num_env
    1906             : 
    1907             :                ! generate one file per processor
    1908             :                stmp = ""
    1909             :                WRITE (stmp, *) irank
    1910             :                my_middle_name = "helium-force-inst-"//TRIM(ADJUSTL(stmp))
    1911             :                unit_nr = cp_print_key_unit_nr( &
    1912             :                          logger, &
    1913             :                          print_key, &
    1914             :                          middle_name=TRIM(my_middle_name), &
    1915             :                          extension=".dat")
    1916             : 
    1917             :                ! unpack and actually print the forces - all components in one line
    1918             :                offset = (irank - 1)*SIZE(helium%rtmp_p_ndim_1d)
    1919             :                idim = 0
    1920             :                DO ib = 1, helium%solute_beads
    1921             :                   DO ia = 1, helium%solute_atoms
    1922             :                      DO ic = 1, 3
    1923             :                         idim = idim + 1
    1924             :                         WRITE (unit_nr, '(F20.10)', ADVANCE='NO') helium%rtmp_p_ndim_np_1d(offset + idim)
    1925             :                      END DO
    1926             :                   END DO
    1927             :                END DO
    1928             :                WRITE (unit_nr, *)
    1929             : 
    1930             :                CALL m_flush(unit_nr)
    1931             :                CALL cp_print_key_finished_output(unit_nr, logger, print_key)
    1932             : 
    1933             :             END DO
    1934             : 
    1935             :          END IF
    1936             :       END IF
    1937             : 
    1938             :       CALL timestop(handle)
    1939             : 
    1940             :    END SUBROUTINE helium_print_force_inst
    1941             : 
    1942             : #endif
    1943             : 
    1944             : END MODULE helium_io

Generated by: LCOV version 1.15