LCOV - code coverage report
Current view: top level - src/motion - helium_io.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 68.6 % 809 555
Test Date: 2025-07-25 12:55:17 Functions: 92.9 % 14 13

            Line data    Source code
       1              : !--------------------------------------------------------------------------------------------------!
       2              : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3              : !   Copyright 2000-2025 CP2K developers group <https://cp2k.org>                                   !
       4              : !                                                                                                  !
       5              : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6              : !--------------------------------------------------------------------------------------------------!
       7              : 
       8              : ! **************************************************************************************************
       9              : !> \brief  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 2.0-1