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