1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Routines that work on qs_subsys_type 8!> \author Ole Schuett 9! ************************************************************************************************** 10MODULE qs_subsys_methods 11 USE atom_types, ONLY: lmat 12 USE atomic_kind_types, ONLY: atomic_kind_type,& 13 get_atomic_kind 14 USE basis_set_types, ONLY: get_gto_basis_set,& 15 gto_basis_set_type 16 USE cell_methods, ONLY: read_cell,& 17 write_cell 18 USE cell_types, ONLY: cell_clone,& 19 cell_create,& 20 cell_release,& 21 cell_type 22 USE cp_para_types, ONLY: cp_para_env_type 23 USE cp_subsys_methods, ONLY: cp_subsys_create 24 USE cp_subsys_types, ONLY: cp_subsys_get,& 25 cp_subsys_release,& 26 cp_subsys_set,& 27 cp_subsys_type 28 USE external_potential_types, ONLY: all_potential_type,& 29 get_potential,& 30 gth_potential_type,& 31 sgp_potential_type 32 USE input_section_types, ONLY: section_vals_get_subs_vals,& 33 section_vals_type 34 USE kinds, ONLY: dp 35 USE molecule_kind_types, ONLY: get_molecule_kind,& 36 molecule_kind_type,& 37 set_molecule_kind 38 USE qs_kind_types, ONLY: create_qs_kind_set,& 39 get_qs_kind,& 40 init_atom_electronic_state,& 41 qs_kind_type 42 USE qs_subsys_types, ONLY: qs_subsys_set,& 43 qs_subsys_type 44#include "./base/base_uses.f90" 45 46 IMPLICIT NONE 47 PRIVATE 48 49 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_subsys_methods' 50 51 PUBLIC :: qs_subsys_create 52 53CONTAINS 54 55! ************************************************************************************************** 56!> \brief Creates a qs_subsys. Optionally an existsing cp_subsys is used. 57!> \param subsys ... 58!> \param para_env ... 59!> \param root_section ... 60!> \param force_env_section ... 61!> \param subsys_section ... 62!> \param use_motion_section ... 63!> \param cp_subsys ... 64!> \param cell ... 65!> \param cell_ref ... 66!> \param elkind ... 67! ************************************************************************************************** 68 SUBROUTINE qs_subsys_create(subsys, para_env, root_section, force_env_section, subsys_section, & 69 use_motion_section, cp_subsys, cell, cell_ref, elkind) 70 TYPE(qs_subsys_type), POINTER :: subsys 71 TYPE(cp_para_env_type), POINTER :: para_env 72 TYPE(section_vals_type), OPTIONAL, POINTER :: root_section 73 TYPE(section_vals_type), POINTER :: force_env_section, subsys_section 74 LOGICAL, INTENT(IN) :: use_motion_section 75 TYPE(cp_subsys_type), OPTIONAL, POINTER :: cp_subsys 76 TYPE(cell_type), OPTIONAL, POINTER :: cell, cell_ref 77 LOGICAL, INTENT(IN), OPTIONAL :: elkind 78 79 CHARACTER(len=*), PARAMETER :: routineN = 'qs_subsys_create', & 80 routineP = moduleN//':'//routineN 81 82 LOGICAL :: use_ref_cell 83 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 84 TYPE(cell_type), POINTER :: my_cell, my_cell_ref 85 TYPE(cp_subsys_type), POINTER :: my_cp_subsys 86 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 87 TYPE(section_vals_type), POINTER :: cell_section, kind_section 88 89 NULLIFY (atomic_kind_set, qs_kind_set, cell_section, kind_section, my_cell, my_cell_ref, my_cp_subsys) 90 91 IF (ASSOCIATED(subsys)) CPABORT("qs_subsys_create: subsys already associated") 92 93 ! create cp_subsys 94 IF (PRESENT(cp_subsys)) THEN 95 my_cp_subsys => cp_subsys 96 ELSE IF (PRESENT(root_section)) THEN 97 CALL cp_subsys_create(my_cp_subsys, para_env, root_section=root_section, & 98 force_env_section=force_env_section, & 99 subsys_section=subsys_section, & 100 use_motion_section=use_motion_section, & 101 elkind=elkind) 102 ELSE 103 CPABORT("qs_subsys_create: cp_subsys or root_section needed") 104 END IF 105 106 ! create cp_subsys%cell 107 !TODO: moved to cp_subsys_create(), needs further disentanglement of cell_ref. 108 IF (PRESENT(cell)) THEN 109 my_cell => cell 110 IF (PRESENT(cell_ref)) THEN 111 my_cell_ref => cell_ref 112 use_ref_cell = .TRUE. 113 ELSE 114 CALL cell_create(my_cell_ref) 115 CALL cell_clone(my_cell, my_cell_ref) 116 use_ref_cell = .FALSE. 117 END IF 118 ELSE 119 cell_section => section_vals_get_subs_vals(subsys_section, "CELL") 120 CALL read_cell(my_cell, my_cell_ref, use_ref_cell=use_ref_cell, & 121 cell_section=cell_section, para_env=para_env) 122 END IF 123 CALL cp_subsys_set(my_cp_subsys, cell=my_cell) 124 CALL write_cell(my_cell, subsys_section, cell_ref=my_cell_ref) 125 126 ! setup qs_kinds 127 CALL cp_subsys_get(my_cp_subsys, atomic_kind_set=atomic_kind_set) 128 kind_section => section_vals_get_subs_vals(subsys_section, "KIND") 129 CALL create_qs_kind_set(qs_kind_set, atomic_kind_set, kind_section, & 130 para_env, force_env_section) 131 132 CALL num_ao_el_per_molecule(my_cp_subsys%molecule_kinds%els, & 133 qs_kind_set) 134 135 ! finally create qs_subsys 136 ALLOCATE (subsys) 137 CALL qs_subsys_set(subsys, & 138 cp_subsys=my_cp_subsys, & 139 cell_ref=my_cell_ref, & 140 use_ref_cell=use_ref_cell, & 141 qs_kind_set=qs_kind_set) 142 143 IF (.NOT. PRESENT(cell)) CALL cell_release(my_cell) 144 IF (.NOT. PRESENT(cell_ref)) CALL cell_release(my_cell_ref) 145 IF (.NOT. PRESENT(cp_subsys)) CALL cp_subsys_release(my_cp_subsys) 146 END SUBROUTINE qs_subsys_create 147 148! ************************************************************************************************** 149!> \brief Read a molecule kind data set from the input file. 150!> \param molecule_kind_set ... 151!> \param qs_kind_set ... 152!> \date 22.11.2004 153!> \par History 154!> Rustam Z. Khaliullin 10.2014 - charges and electrons of molecules 155!> are now in agreement with atomic guess 156!> \author MI 157!> \version 1.0 158! ************************************************************************************************** 159 SUBROUTINE num_ao_el_per_molecule(molecule_kind_set, qs_kind_set) 160 161 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 162 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 163 164 CHARACTER(len=*), PARAMETER :: routineN = 'num_ao_el_per_molecule', & 165 routineP = moduleN//':'//routineN 166 167 INTEGER :: arbitrary_spin, iatom, ikind, imol, & 168 n_ao, natom, nmol_kind, nsgf, nspins, & 169 z_molecule 170 INTEGER, DIMENSION(0:lmat, 10) :: ne_core, ne_elem, ne_explicit 171 INTEGER, DIMENSION(2) :: n_occ_alpha_and_beta 172 REAL(KIND=dp) :: charge_molecule, zeff, zeff_correction 173 REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta 174 TYPE(all_potential_type), POINTER :: all_potential 175 TYPE(atomic_kind_type), POINTER :: atomic_kind 176 TYPE(gth_potential_type), POINTER :: gth_potential 177 TYPE(gto_basis_set_type), POINTER :: orb_basis_set 178 TYPE(molecule_kind_type), POINTER :: molecule_kind 179 TYPE(sgp_potential_type), POINTER :: sgp_potential 180 181 IF (ASSOCIATED(molecule_kind_set)) THEN 182 183 nspins = 2 184 nmol_kind = SIZE(molecule_kind_set, 1) 185 natom = 0 186 187 ! *** Initialize the molecule kind data structure *** 188 ARBITRARY_SPIN = 1 189 DO imol = 1, nmol_kind 190 191 molecule_kind => molecule_kind_set(imol) 192 CALL get_molecule_kind(molecule_kind=molecule_kind, & 193 natom=natom) 194 !nelectron = 0 195 n_ao = 0 196 n_occ_alpha_and_beta(1:nspins) = 0 197 z_molecule = 0 198 199 DO iatom = 1, natom 200 201 atomic_kind => molecule_kind%atom_list(iatom)%atomic_kind 202 CALL get_atomic_kind(atomic_kind, kind_number=ikind) 203 CALL get_qs_kind(qs_kind_set(ikind), & 204 basis_set=orb_basis_set, & 205 all_potential=all_potential, & 206 gth_potential=gth_potential, & 207 sgp_potential=sgp_potential) 208 209 ! Obtain the electronic state of the atom 210 ! The same state is used to calculate the ATOMIC GUESS 211 ! It is great that we are consistent with ATOMIC_GUESS 212 CALL init_atom_electronic_state(atomic_kind=atomic_kind, & 213 qs_kind=qs_kind_set(ikind), & 214 ncalc=ne_explicit, & 215 ncore=ne_core, & 216 nelem=ne_elem, & 217 edelta=edelta) 218 219 ! If &BS section is used ATOMIC_GUESS is calculated twice 220 ! for two separate wfns with their own alpha-beta combinations 221 ! This is done to break the spin symmetry of the initial wfn 222 ! For now, only alpha part of &BS is used to count electrons on 223 ! molecules 224 ! Get the number of explicit electrons (i.e. with orbitals) 225 ! For now, only the total number of electrons can be obtained 226 ! from init_atom_electronic_state 227 n_occ_alpha_and_beta(ARBITRARY_SPIN) = & 228 n_occ_alpha_and_beta(ARBITRARY_SPIN) + SUM(ne_explicit) + & 229 SUM(NINT(2*edelta(:, :, ARBITRARY_SPIN))) 230 ! We need a way to specify the number of alpha and beta electrons 231 ! on each molecule (i.e. multiplicity is not enough) 232 !n_occ(ispin) = n_occ(ispin) + SUM(ne_explicit) + SUM(NINT(2*edelta(:, :, ispin))) 233 234 IF (ASSOCIATED(all_potential)) THEN 235 CALL get_potential(potential=all_potential, zeff=zeff, & 236 zeff_correction=zeff_correction) 237 ELSE IF (ASSOCIATED(gth_potential)) THEN 238 CALL get_potential(potential=gth_potential, zeff=zeff, & 239 zeff_correction=zeff_correction) 240 ELSE IF (ASSOCIATED(sgp_potential)) THEN 241 CALL get_potential(potential=sgp_potential, zeff=zeff, & 242 zeff_correction=zeff_correction) 243 ELSE 244 zeff = 0.0_dp 245 zeff_correction = 0.0_dp 246 END IF 247 z_molecule = z_molecule + NINT(zeff - zeff_correction) 248 249 ! this one does not work because nelem is not adjusted in the symmetry breaking code 250 !CALL get_atomic_kind(atomic_kind,z=z) 251 !z_molecule=z_molecule+z 252 253 IF (ASSOCIATED(orb_basis_set)) THEN 254 CALL get_gto_basis_set(gto_basis_set=orb_basis_set, nsgf=nsgf) 255 ELSE 256 nsgf = 0 257 ENDIF 258 n_ao = n_ao + nsgf 259 260 END DO ! iatom 261 262 ! At this point we have the number of electrons (alpha+beta) on the molecule 263 ! as they are seen by the ATOMIC GUESS routines 264 charge_molecule = REAL(z_molecule - n_occ_alpha_and_beta(ARBITRARY_SPIN), dp) 265 CALL set_molecule_kind(molecule_kind=molecule_kind, & 266 nelectron=n_occ_alpha_and_beta(ARBITRARY_SPIN), & 267 charge=charge_molecule, & 268 nsgf=n_ao) 269 270 END DO ! imol 271 END IF 272 273 END SUBROUTINE num_ao_el_per_molecule 274 275END MODULE qs_subsys_methods 276