1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Utilities for Molecular Dynamics
8!> \author Teodoro Laino [tlaino] - University of Zurich - 09.2007
9! **************************************************************************************************
10MODULE md_util
11
12   USE cp_files,                        ONLY: close_file,&
13                                              open_file
14   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
15                                              cp_logger_type
16   USE cp_output_handling,              ONLY: cp_print_key_generate_filename
17   USE cp_para_types,                   ONLY: cp_para_env_type
18   USE input_cp2k_restarts,             ONLY: write_restart
19   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
20                                              section_vals_type,&
21                                              section_vals_val_get
22   USE kinds,                           ONLY: default_path_length,&
23                                              dp
24   USE md_energies,                     ONLY: md_write_output
25   USE md_environment_types,            ONLY: md_environment_type
26   USE message_passing,                 ONLY: mp_bcast
27#include "../base/base_uses.f90"
28
29   IMPLICIT NONE
30
31   PRIVATE
32
33! *** Global parameters ***
34
35   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'md_util'
36
37   PUBLIC :: md_output, &
38             read_vib_eigs_unformatted
39
40CONTAINS
41
42! **************************************************************************************************
43!> \brief collects the part of the MD that, basically, does the output
44!> \param md_env ...
45!> \param md_section ...
46!> \param root_section ...
47!> \param forced_io ...
48!> \par History
49!>      03.2006 created [Joost VandeVondele]
50! **************************************************************************************************
51   SUBROUTINE md_output(md_env, md_section, root_section, forced_io)
52      TYPE(md_environment_type), POINTER                 :: md_env
53      TYPE(section_vals_type), POINTER                   :: md_section, root_section
54      LOGICAL, INTENT(IN)                                :: forced_io
55
56      CHARACTER(LEN=*), PARAMETER :: routineN = 'md_output', routineP = moduleN//':'//routineN
57
58      INTEGER                                            :: handle
59      LOGICAL                                            :: do_print
60      TYPE(section_vals_type), POINTER                   :: print_section
61
62      CALL timeset(routineN, handle)
63      do_print = .TRUE.
64      IF (forced_io) THEN
65         print_section => section_vals_get_subs_vals(md_section, "PRINT")
66         CALL section_vals_val_get(print_section, "FORCE_LAST", l_val=do_print)
67      END IF
68      IF (do_print) THEN
69         ! Dumps all files related to the MD run
70         CALL md_write_output(md_env)
71         CALL write_restart(md_env=md_env, root_section=root_section)
72      END IF
73      CALL timestop(handle)
74
75   END SUBROUTINE md_output
76
77! **************************************************************************************************
78!> \brief read eigenvalues and eigenvectors of Hessian from vibrational analysis results, for use
79!>        of initialising MD simulations. Expects to read an unformatted binary file
80!> \param md_section : input section object containing MD subsections and keywords. This should
81!>                     provide the filename to read vib analysis eigenvalues and eigenvectors.
82!>                     If the filename is not explicitly specified by the user in the input, then
83!>                     it will use the default CARTESIAN_EIGS print key filename defined in the
84!>                     vibrational analysis input section as the filename.
85!> \param vib_section : input section object containing vibrational analysis subsections
86!>                      and keywords
87!> \param para_env : cp2k mpi environment object, needed for IO in parallel computations
88!> \param dof : outputs the total number of eigenvalues (no. degrees of freedom) read from the file
89!> \param eigenvalues : outputs the eigenvalues (Cartesian frequencies) read from the file
90!> \param eigenvectors : outputs the corresponding eigenvectors read from the file
91!> \author Lianheng Tong, lianheng.tong@kcl.ac.uk
92! **************************************************************************************************
93   SUBROUTINE read_vib_eigs_unformatted(md_section, &
94                                        vib_section, &
95                                        para_env, &
96                                        dof, &
97                                        eigenvalues, &
98                                        eigenvectors)
99      TYPE(section_vals_type), POINTER                   :: md_section, vib_section
100      TYPE(cp_para_env_type), POINTER                    :: para_env
101      INTEGER, INTENT(OUT)                               :: dof
102      REAL(KIND=dp), DIMENSION(:), INTENT(OUT)           :: eigenvalues
103      REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT)        :: eigenvectors
104
105      CHARACTER(len=*), PARAMETER :: routineN = 'read_vib_eigs_unformatted', &
106         routineP = moduleN//':'//routineN
107
108      CHARACTER(LEN=default_path_length)                 :: filename
109      INTEGER                                            :: jj, n_rep_val, unit_nr
110      LOGICAL                                            :: exist
111      TYPE(cp_logger_type), POINTER                      :: logger
112      TYPE(section_vals_type), POINTER                   :: print_key
113
114      logger => cp_get_default_logger()
115      dof = 0
116      eigenvalues = 0.0_dp
117      eigenvectors = 0.0_dp
118      ! obtain file name
119      CALL section_vals_val_get(md_section, "INITIAL_VIBRATION%VIB_EIGS_FILE_NAME", &
120                                n_rep_val=n_rep_val)
121      IF (n_rep_val > 0) THEN
122         CALL section_vals_val_get(md_section, "VIB_EIGS_FILE_NAME", c_val=filename)
123      ELSE
124         print_key => section_vals_get_subs_vals(vib_section, "PRINT%CARTESIAN_EIGS")
125         filename = cp_print_key_generate_filename(logger, print_key, extension="eig", &
126                                                   my_local=.FALSE.)
127      END IF
128      ! read file
129      IF (para_env%ionode) THEN
130         INQUIRE (FILE=filename, exist=exist)
131         IF (.NOT. exist) THEN
132            CPABORT("File "//filename//" is not found.")
133         END IF
134         CALL open_file(file_name=filename, &
135                        file_action="READ", &
136                        file_form="UNFORMATTED", &
137                        file_status="OLD", &
138                        unit_number=unit_nr)
139         ! the first record contains one integer giving degrees of freedom
140         READ (unit_nr) dof
141         IF (dof .GT. SIZE(eigenvalues)) THEN
142            CPABORT("Too many DoFs found in "//filename)
143         END IF
144         ! the second record contains the eigenvalues
145         READ (unit_nr) eigenvalues(1:dof)
146         ! the rest of the records contain the eigenvectors
147         DO jj = 1, dof
148            READ (unit_nr) eigenvectors(1:dof, jj)
149         END DO
150      END IF
151      ! broadcast to all compulational nodes. note that it is assumed
152      ! that source is the ionode
153      CALL mp_bcast(dof, para_env%source, para_env%group)
154      CALL mp_bcast(eigenvalues, para_env%source, para_env%group)
155      CALL mp_bcast(eigenvectors, para_env%source, para_env%group)
156      ! close file
157      IF (para_env%ionode) THEN
158         CALL close_file(unit_number=unit_nr)
159      END IF
160   END SUBROUTINE read_vib_eigs_unformatted
161
162END MODULE md_util
163