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