1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Methods and functions on the PWDFT environment
8!> \par History
9!>      07.2018 initial create
10!> \author JHU
11! **************************************************************************************************
12MODULE pwdft_environment
13   USE atomic_kind_types,               ONLY: atomic_kind_type
14   USE cell_types,                      ONLY: cell_type
15   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
16                                              cp_logger_get_default_io_unit,&
17                                              cp_logger_type
18   USE cp_para_types,                   ONLY: cp_para_env_type
19   USE cp_symmetry,                     ONLY: write_symmetry
20   USE distribution_1d_types,           ONLY: distribution_1d_release,&
21                                              distribution_1d_type
22   USE distribution_methods,            ONLY: distribute_molecules_1d
23   USE header,                          ONLY: sirius_header
24   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
25                                              section_vals_type
26   USE kinds,                           ONLY: dp
27   USE machine,                         ONLY: m_flush
28   USE molecule_kind_types,             ONLY: molecule_kind_type
29   USE molecule_types,                  ONLY: molecule_type
30   USE particle_methods,                ONLY: write_particle_distances,&
31                                              write_qs_particle_coordinates,&
32                                              write_structure_data
33   USE particle_types,                  ONLY: particle_type
34   USE pwdft_environment_types,         ONLY: pwdft_env_get,&
35                                              pwdft_env_set,&
36                                              pwdft_environment_type
37   USE qs_energy_types,                 ONLY: allocate_qs_energy,&
38                                              qs_energy_type
39   USE qs_kind_types,                   ONLY: qs_kind_type,&
40                                              write_qs_kind_set
41   USE qs_subsys_methods,               ONLY: qs_subsys_create
42   USE qs_subsys_types,                 ONLY: qs_subsys_get,&
43                                              qs_subsys_release,&
44                                              qs_subsys_set,&
45                                              qs_subsys_type
46   USE sirius_interface,                ONLY: cp_sirius_create_env,&
47                                              cp_sirius_energy_force,&
48                                              cp_sirius_update_context
49   USE virial_types,                    ONLY: virial_type
50#include "./base/base_uses.f90"
51
52   IMPLICIT NONE
53
54   PRIVATE
55
56! *** Global parameters ***
57
58   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pwdft_environment'
59
60! *** Public subroutines ***
61
62   PUBLIC :: pwdft_init, pwdft_calc_energy_force
63
64CONTAINS
65
66! **************************************************************************************************
67!> \brief Initialize the pwdft environment
68!> \param pwdft_env The pwdft environment to retain
69!> \param root_section ...
70!> \param para_env ...
71!> \param force_env_section ...
72!> \param subsys_section ...
73!> \param use_motion_section ...
74!> \par History
75!>      03.2018 initial create
76!> \author JHU
77! **************************************************************************************************
78   SUBROUTINE pwdft_init(pwdft_env, root_section, para_env, force_env_section, subsys_section, &
79                         use_motion_section)
80      TYPE(pwdft_environment_type), POINTER              :: pwdft_env
81      TYPE(section_vals_type), POINTER                   :: root_section
82      TYPE(cp_para_env_type), POINTER                    :: para_env
83      TYPE(section_vals_type), POINTER                   :: force_env_section, subsys_section
84      LOGICAL, INTENT(IN)                                :: use_motion_section
85
86      CHARACTER(len=*), PARAMETER :: routineN = 'pwdft_init', routineP = moduleN//':'//routineN
87
88      INTEGER                                            :: handle, iw, natom
89      LOGICAL                                            :: use_ref_cell
90      TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
91      TYPE(cell_type), POINTER                           :: my_cell, my_cell_ref
92      TYPE(cp_logger_type), POINTER                      :: logger
93      TYPE(distribution_1d_type), POINTER                :: local_molecules, local_particles
94      TYPE(molecule_kind_type), DIMENSION(:), POINTER    :: molecule_kind_set
95      TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
96      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
97      TYPE(qs_energy_type), POINTER                      :: energy
98      TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
99      TYPE(qs_subsys_type), POINTER                      :: qs_subsys
100      TYPE(section_vals_type), POINTER                   :: pwdft_section, xc_section
101
102      CALL timeset(routineN, handle)
103
104      CPASSERT(ASSOCIATED(pwdft_env))
105      CPASSERT(ASSOCIATED(force_env_section))
106
107      IF (.NOT. ASSOCIATED(subsys_section)) THEN
108         subsys_section => section_vals_get_subs_vals(force_env_section, "SUBSYS")
109      END IF
110      pwdft_section => section_vals_get_subs_vals(force_env_section, "PW_DFT")
111
112      ! retrieve the functionals parameters
113      xc_section => section_vals_get_subs_vals(force_env_section, "DFT%XC%XC_FUNCTIONAL")
114
115      CALL pwdft_env_set(pwdft_env=pwdft_env, pwdft_input=pwdft_section, &
116                         force_env_input=force_env_section, xc_input=xc_section)
117
118      ! parallel environment
119      CALL pwdft_env_set(pwdft_env=pwdft_env, para_env=para_env)
120
121      NULLIFY (qs_subsys)
122      CALL qs_subsys_create(qs_subsys, para_env, &
123                            force_env_section=force_env_section, &
124                            subsys_section=subsys_section, &
125                            use_motion_section=use_motion_section, &
126                            root_section=root_section)
127
128      ! Distribute molecules and atoms
129      NULLIFY (local_molecules)
130      NULLIFY (local_particles)
131      CALL qs_subsys_get(qs_subsys, particle_set=particle_set, &
132                         atomic_kind_set=atomic_kind_set, &
133                         molecule_set=molecule_set, &
134                         molecule_kind_set=molecule_kind_set)
135
136      CALL distribute_molecules_1d(atomic_kind_set=atomic_kind_set, &
137                                   particle_set=particle_set, &
138                                   local_particles=local_particles, &
139                                   molecule_kind_set=molecule_kind_set, &
140                                   molecule_set=molecule_set, &
141                                   local_molecules=local_molecules, &
142                                   force_env_section=force_env_section)
143
144      CALL qs_subsys_set(qs_subsys, local_molecules=local_molecules, &
145                         local_particles=local_particles)
146      CALL distribution_1d_release(local_particles)
147      CALL distribution_1d_release(local_molecules)
148
149      CALL pwdft_env_set(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
150      CALL qs_subsys_release(qs_subsys)
151
152      CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
153      CALL qs_subsys_get(qs_subsys, cell=my_cell, cell_ref=my_cell_ref, use_ref_cell=use_ref_cell)
154
155      NULLIFY (logger)
156      logger => cp_get_default_logger()
157      iw = cp_logger_get_default_io_unit(logger)
158      CALL sirius_header(iw)
159
160      NULLIFY (energy)
161      CALL allocate_qs_energy(energy)
162      CALL qs_subsys_set(qs_subsys, energy=energy)
163
164      CALL qs_subsys_get(qs_subsys, particle_set=particle_set, &
165                         qs_kind_set=qs_kind_set, &
166                         atomic_kind_set=atomic_kind_set, &
167                         molecule_set=molecule_set, &
168                         molecule_kind_set=molecule_kind_set)
169
170      ! init energy, force, stress arrays
171      CALL qs_subsys_get(qs_subsys, natom=natom)
172      ALLOCATE (pwdft_env%energy)
173      ALLOCATE (pwdft_env%forces(natom, 3))
174      pwdft_env%forces(1:natom, 1:3) = 0.0_dp
175      pwdft_env%stress(1:3, 1:3) = 0.0_dp
176
177      ! Print the atomic kind set
178      CALL write_qs_kind_set(qs_kind_set, subsys_section)
179      ! Print the atomic coordinates
180      CALL write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label="PWDFT/SIRIUS")
181      ! Print the interatomic distances
182      CALL write_particle_distances(particle_set, my_cell, subsys_section)
183      ! Print the requested structure data
184      CALL write_structure_data(particle_set, my_cell, subsys_section)
185      ! Print symmetry information
186      CALL write_symmetry(particle_set, my_cell, subsys_section)
187
188      IF (iw > 0) THEN
189         WRITE (iw, '(A,A,A)') " ==================================", " SIRIUS INIT ", &
190            "================================="
191      END IF
192      ! Sirius initialization
193      CALL cp_sirius_create_env(pwdft_env)
194
195      IF (iw > 0) THEN
196         WRITE (iw, '(A,A,A)') " =========================", " SIRIUS INIT FINISHED ", &
197            "================================="
198      END IF
199      IF (iw > 0) CALL m_flush(iw)
200
201      CALL timestop(handle)
202
203   END SUBROUTINE pwdft_init
204
205! **************************************************************************************************
206!> \brief Calculate energy and forces within the PWDFT/SIRIUS code
207!> \param pwdft_env The pwdft environment to retain
208!> \param calculate_forces ...
209!> \param calculate_stress ...
210!> \par History
211!>      03.2018 initial create
212!> \author JHU
213! **************************************************************************************************
214   SUBROUTINE pwdft_calc_energy_force(pwdft_env, calculate_forces, calculate_stress)
215      TYPE(pwdft_environment_type), POINTER              :: pwdft_env
216      LOGICAL, INTENT(IN)                                :: calculate_forces, calculate_stress
217
218      CHARACTER(len=*), PARAMETER :: routineN = 'pwdft_calc_energy_force', &
219         routineP = moduleN//':'//routineN
220
221      INTEGER                                            :: handle, iatom, iw, natom
222      REAL(KIND=dp), DIMENSION(1:3, 1:3)                 :: stress
223      REAL(KIND=dp), DIMENSION(:, :), POINTER            :: force
224      TYPE(cell_type), POINTER                           :: my_cell
225      TYPE(cp_logger_type), POINTER                      :: logger
226      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
227      TYPE(qs_subsys_type), POINTER                      :: qs_subsys
228      TYPE(virial_type), POINTER                         :: virial
229
230      CALL timeset(routineN, handle)
231
232      CPASSERT(ASSOCIATED(pwdft_env))
233
234      NULLIFY (logger)
235      logger => cp_get_default_logger()
236      iw = cp_logger_get_default_io_unit(logger)
237
238      ! update context for new positions/cell
239      IF (iw > 0) THEN
240         WRITE (iw, '(A,A,A)') " =============================", " SIRIUS UPDATE CONTEXT", &
241            "============================"
242      END IF
243      ! Sirius initialization
244      CALL cp_sirius_update_context(pwdft_env)
245
246      IF (iw > 0) THEN
247         WRITE (iw, '(A,A,A)') " ====================", " SIRIUS UPDATE CONTEXT FINISHED", &
248            "============================"
249      END IF
250      IF (iw > 0) CALL m_flush(iw)
251
252      ! calculate energy and forces/stress
253      CALL cp_sirius_energy_force(pwdft_env, calculate_forces, calculate_stress)
254
255      IF (calculate_forces) THEN
256         CALL pwdft_env_get(pwdft_env=pwdft_env, qs_subsys=qs_subsys)
257         CALL pwdft_env_get(pwdft_env=pwdft_env, forces=force)
258         CALL qs_subsys_get(qs_subsys, particle_set=particle_set, natom=natom)
259         DO iatom = 1, natom
260            particle_set(iatom)%f(1:3) = -force(iatom, 1:3)
261         END DO
262      END IF
263
264      IF (calculate_stress) THEN
265         ! i need to retrieve the volume of the unit cell for the stress tensor
266         CALL qs_subsys_get(qs_subsys, cell=my_cell, virial=virial)
267         CALL pwdft_env_get(pwdft_env=pwdft_env, stress=stress)
268         virial%pv_virial(1:3, 1:3) = -stress(1:3, 1:3)*my_cell%deth
269      END IF
270
271      CALL timestop(handle)
272
273   END SUBROUTINE pwdft_calc_energy_force
274END MODULE pwdft_environment
275