1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Driver for the localization that should be general 8!> for all the methods available and all the definition of the 9!> spread functional 10!> Write centers, spread and cubes only if required and for the 11!> selected states 12!> The localized functions are copied in the standard mos array 13!> for the next use 14!> \par History 15!> 01.2008 Teodoro Laino [tlaino] - University of Zurich 16!> - Merging the two localization codes and updating to new structures 17!> \author MI (04.2005) 18! ************************************************************************************************** 19MODULE qs_loc_methods 20 USE atomic_kind_types, ONLY: atomic_kind_type,& 21 deallocate_atomic_kind_set,& 22 set_atomic_kind 23 USE cell_types, ONLY: cell_type,& 24 pbc 25 USE cp_control_types, ONLY: dft_control_type 26 USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm 27 USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,& 28 fm_pool_create_fm 29 USE cp_fm_struct, ONLY: cp_fm_struct_create,& 30 cp_fm_struct_release,& 31 cp_fm_struct_type 32 USE cp_fm_types, ONLY: & 33 cp_fm_create, cp_fm_get_element, cp_fm_get_info, cp_fm_get_submatrix, cp_fm_p_type, & 34 cp_fm_release, cp_fm_set_all, cp_fm_set_submatrix, cp_fm_to_fm, cp_fm_type 35 USE cp_gemm_interface, ONLY: cp_gemm 36 USE cp_log_handling, ONLY: cp_get_default_logger,& 37 cp_logger_type,& 38 cp_to_string 39 USE cp_output_handling, ONLY: cp_iter_string,& 40 cp_p_file,& 41 cp_print_key_finished_output,& 42 cp_print_key_should_output,& 43 cp_print_key_unit_nr 44 USE cp_para_types, ONLY: cp_para_env_type 45 USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube 46 USE cp_units, ONLY: cp_unit_from_cp2k 47 USE dbcsr_api, ONLY: dbcsr_p_type 48 USE input_constants, ONLY: & 49 do_loc_crazy, do_loc_direct, do_loc_jacobi, do_loc_l1_norm_sd, do_loc_none, do_loc_scdm, & 50 dump_dcd, dump_dcd_aligned_cell, dump_xmol, op_loc_berry, op_loc_boys, op_loc_pipek, & 51 state_loc_list 52 USE input_section_types, ONLY: section_get_ival,& 53 section_get_ivals,& 54 section_get_lval,& 55 section_vals_get_subs_vals,& 56 section_vals_type,& 57 section_vals_val_get 58 USE kinds, ONLY: default_path_length,& 59 default_string_length,& 60 dp,& 61 sp 62 USE machine, ONLY: m_flush 63 USE mathconstants, ONLY: twopi 64 USE memory_utilities, ONLY: reallocate 65 USE motion_utils, ONLY: get_output_format 66 USE particle_list_types, ONLY: particle_list_type 67 USE particle_methods, ONLY: get_particle_set,& 68 write_particle_coordinates 69 USE particle_types, ONLY: allocate_particle_set,& 70 deallocate_particle_set,& 71 particle_type 72 USE physcon, ONLY: angstrom 73 USE pw_env_types, ONLY: pw_env_get,& 74 pw_env_type 75 USE pw_pool_types, ONLY: pw_pool_create_pw,& 76 pw_pool_give_back_pw,& 77 pw_pool_type 78 USE pw_types, ONLY: COMPLEXDATA1D,& 79 REALDATA3D,& 80 REALSPACE,& 81 RECIPROCALSPACE,& 82 pw_p_type 83 USE qs_collocate_density, ONLY: calculate_wavefunction 84 USE qs_environment_types, ONLY: get_qs_env,& 85 qs_environment_type 86 USE qs_kind_types, ONLY: qs_kind_type 87 USE qs_loc_types, ONLY: get_qs_loc_env,& 88 localized_wfn_control_type,& 89 qs_loc_env_new_type 90 USE qs_loc_utils, ONLY: jacobi_rotation_pipek 91 USE qs_localization_methods, ONLY: approx_l1_norm_sd,& 92 crazy_rotations,& 93 direct_mini,& 94 jacobi_rotations,& 95 scdm_qrfact,& 96 zij_matrix 97 USE qs_matrix_pools, ONLY: mpools_get 98 USE qs_mo_types, ONLY: get_mo_set,& 99 mo_set_p_type 100 USE qs_subsys_types, ONLY: qs_subsys_get,& 101 qs_subsys_type 102 USE string_utilities, ONLY: xstring 103#include "./base/base_uses.f90" 104 105 IMPLICIT NONE 106 107 PRIVATE 108 109! *** Global parameters *** 110 111 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_methods' 112 113! *** Public *** 114 PUBLIC :: qs_loc_driver, qs_print_cubes, centers_spreads_berry 115 116CONTAINS 117 118! ************************************************************************************************** 119!> \brief Calculate and optimize the spread functional as calculated from 120!> the selected mos and by the definition using the berry phase 121!> as given by silvestrelli et al 122!> If required the centers and the spreads for each mos selected 123!> are calculated from z_ij and printed to file. 124!> The centers files is appended if already exists 125!> \param method indicates localization algorithm 126!> \param qs_loc_env new environment for the localization calculations 127!> \param vectors selected mos to be localized 128!> \param op_sm_set sparse matrices containing the integrals of the kind <mi e{iGr} nu> 129!> \param zij_fm_set set of full matrix of size nmoloc x nmoloc, will contain the z_ij numbers 130!> as defined by Silvestrelli et al 131!> \param para_env ... 132!> \param cell ... 133!> \param weights ... 134!> \param ispin ... 135!> \param print_loc_section ... 136!> \par History 137!> 04.2005 created [MI] 138!> \author MI 139!> \note 140!> This definition need the use of complex numbers, therefore the 141!> optimization routines are specific for this case 142!> The file for the centers and the spreads have a xyz format 143! ************************************************************************************************** 144 SUBROUTINE optimize_loc_berry(method, qs_loc_env, vectors, op_sm_set, & 145 zij_fm_set, para_env, cell, weights, ispin, print_loc_section) 146 147 INTEGER, INTENT(IN) :: method 148 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 149 TYPE(cp_fm_type), POINTER :: vectors 150 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set 151 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: zij_fm_set 152 TYPE(cp_para_env_type), POINTER :: para_env 153 TYPE(cell_type), POINTER :: cell 154 REAL(dp), DIMENSION(:) :: weights 155 INTEGER, INTENT(IN) :: ispin 156 TYPE(section_vals_type), POINTER :: print_loc_section 157 158 CHARACTER(len=*), PARAMETER :: routineN = 'optimize_loc_berry', & 159 routineP = moduleN//':'//routineN 160 161 INTEGER :: handle, max_iter, nao, nmoloc, out_each, & 162 output_unit, sweeps 163 LOGICAL :: converged, crazy_use_diag, & 164 do_jacobi_refinement 165 REAL(dp) :: crazy_scale, eps_localization, & 166 max_crazy_angle, start_time, & 167 target_time 168 TYPE(cp_logger_type), POINTER :: logger 169 170! INTEGER :: i,j 171 172 CALL timeset(routineN, handle) 173 logger => cp_get_default_logger() 174 output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", & 175 extension=".locInfo") 176 177 ! get rows and cols of the input 178 CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc) 179 180 CALL zij_matrix(vectors, op_sm_set, zij_fm_set) 181 182 max_iter = qs_loc_env%localized_wfn_control%max_iter 183 max_crazy_angle = qs_loc_env%localized_wfn_control%max_crazy_angle 184 crazy_use_diag = qs_loc_env%localized_wfn_control%crazy_use_diag 185 crazy_scale = qs_loc_env%localized_wfn_control%crazy_scale 186 eps_localization = qs_loc_env%localized_wfn_control%eps_localization 187 out_each = qs_loc_env%localized_wfn_control%out_each 188 target_time = qs_loc_env%target_time 189 start_time = qs_loc_env%start_time 190 do_jacobi_refinement = qs_loc_env%localized_wfn_control%jacobi_refinement 191 CALL centers_spreads_berry(qs_loc_env, zij_fm_set, nmoloc, cell, weights, & 192 ispin, print_loc_section, only_initial_out=.TRUE.) 193 SELECT CASE (method) 194 CASE (do_loc_jacobi) 195 CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, & 196 eps_localization=eps_localization, sweeps=sweeps, & 197 out_each=out_each, target_time=target_time, start_time=start_time) 198 CASE (do_loc_scdm) 199 ! Decomposition 200 CALL scdm_qrfact(vectors) 201 ! Calculation of Zij 202 CALL zij_matrix(vectors, op_sm_set, zij_fm_set) 203 IF (do_jacobi_refinement) THEN 204 ! Intermediate spread and centers 205 CALL centers_spreads_berry(qs_loc_env, zij_fm_set, nmoloc, cell, weights, & 206 ispin, print_loc_section, only_initial_out=.TRUE.) 207 CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, & 208 eps_localization=eps_localization, sweeps=sweeps, & 209 out_each=out_each, target_time=target_time, start_time=start_time) 210 END IF 211 CASE (do_loc_crazy) 212 CALL crazy_rotations(weights, zij_fm_set, vectors, max_iter=max_iter, max_crazy_angle=max_crazy_angle, & 213 crazy_scale=crazy_scale, crazy_use_diag=crazy_use_diag, & 214 eps_localization=eps_localization, iterations=sweeps, converged=converged) 215 ! Possibly fallback to jacobi if the crazy rotation fails 216 IF (.NOT. converged) THEN 217 IF (qs_loc_env%localized_wfn_control%jacobi_fallback) THEN 218 IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") & 219 " Crazy Wannier localization not converged after ", sweeps, & 220 " iterations, switching to jacobi rotations" 221 CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, & 222 eps_localization=eps_localization, sweeps=sweeps, & 223 out_each=out_each, target_time=target_time, start_time=start_time) 224 ELSE 225 IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") & 226 " Crazy Wannier localization not converged after ", sweeps, & 227 " iterations, and jacobi_fallback switched off " 228 ENDIF 229 END IF 230 CASE (do_loc_direct) 231 CALL direct_mini(weights, zij_fm_set, vectors, max_iter=max_iter, & 232 eps_localization=eps_localization, iterations=sweeps) 233 CASE (do_loc_l1_norm_sd) 234 IF (.NOT. cell%orthorhombic) THEN 235 CPABORT("Non-orthorhombic cell with the selected method NYI") 236 ELSE 237 CALL approx_l1_norm_sd(vectors, max_iter, eps_localization, converged, sweeps) 238 ! here we need to set zij for the computation of the centers and spreads 239 CALL zij_matrix(vectors, op_sm_set, zij_fm_set) 240 END IF 241 CASE (do_loc_none) 242 IF (output_unit > 0) THEN 243 WRITE (output_unit, '(T4,A,I6,A)') " No MOS localization applied " 244 ENDIF 245 CASE DEFAULT 246 CPABORT("Unknown localization method") 247 END SELECT 248 IF (output_unit > 0) THEN 249 IF (sweeps <= max_iter) WRITE (output_unit, '(T4,A,I3,A,I6,A)') " Localization for spin ", ispin, & 250 " converged in ", sweeps, " iterations" 251 ENDIF 252 253 CALL centers_spreads_berry(qs_loc_env, zij_fm_set, nmoloc, cell, weights, & 254 ispin, print_loc_section) 255 CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO") 256 257 CALL timestop(handle) 258 259 END SUBROUTINE optimize_loc_berry 260 261! ************************************************************************************************** 262!> \brief ... 263!> \param qs_env ... 264!> \param method ... 265!> \param qs_loc_env ... 266!> \param vectors ... 267!> \param zij_fm_set ... 268!> \param ispin ... 269!> \param print_loc_section ... 270!> \par History 271!> 04.2005 created [MI] 272!> \author MI 273! ************************************************************************************************** 274 SUBROUTINE optimize_loc_pipek(qs_env, method, qs_loc_env, vectors, zij_fm_set, & 275 ispin, print_loc_section) 276 TYPE(qs_environment_type), POINTER :: qs_env 277 INTEGER, INTENT(IN) :: method 278 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 279 TYPE(cp_fm_type), POINTER :: vectors 280 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: zij_fm_set 281 INTEGER, INTENT(IN) :: ispin 282 TYPE(section_vals_type), POINTER :: print_loc_section 283 284 CHARACTER(len=*), PARAMETER :: routineN = 'optimize_loc_pipek', & 285 routineP = moduleN//':'//routineN 286 287 INTEGER :: handle, iatom, isgf, ldz, nao, natom, & 288 ncol, nmoloc, output_unit, sweeps 289 INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgf 290 TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools 291 TYPE(cp_fm_type), POINTER :: opvec, ov_fm 292 TYPE(cp_logger_type), POINTER :: logger 293 TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s 294 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 295 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 296 297 CALL timeset(routineN, handle) 298 logger => cp_get_default_logger() 299 output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", & 300 extension=".locInfo") 301 302 NULLIFY (particle_set) 303 ! get rows and cols of the input 304 CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc) 305 306 ! replicate the input kind of matrix 307 CALL cp_fm_create(opvec, vectors%matrix_struct) 308 CALL cp_fm_set_all(opvec, 0.0_dp) 309 310 CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, & 311 particle_set=particle_set, qs_kind_set=qs_kind_set) 312 313 natom = SIZE(particle_set, 1) 314 ALLOCATE (first_sgf(natom)) 315 ALLOCATE (last_sgf(natom)) 316 ALLOCATE (nsgf(natom)) 317 318 ! construction of 319 CALL get_particle_set(particle_set, qs_kind_set, & 320 first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf) 321 322 ! Copy the overlap sparse matrix in a full matrix 323 CALL mpools_get(qs_env%mpools, ao_ao_fm_pools=ao_ao_fm_pools) 324 CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, ov_fm, name=" ") 325 CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, ov_fm) 326 327 ! Compute zij here 328 DO iatom = 1, natom 329 CALL cp_fm_set_all(zij_fm_set(iatom, 1)%matrix, 0.0_dp) 330 CALL cp_fm_get_info(zij_fm_set(iatom, 1)%matrix, ncol_global=ldz) 331 isgf = first_sgf(iatom) 332 ncol = nsgf(iatom) 333 334 ! multiply fmxfm, using only part of the ao : Ct x S 335 CALL cp_gemm('N', 'N', nao, nmoloc, nao, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, & 336 a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1) 337 338 CALL cp_gemm('T', 'N', nmoloc, nmoloc, ncol, 0.5_dp, vectors, opvec, & 339 0.0_dp, zij_fm_set(iatom, 1)%matrix, & 340 a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf) 341 342 CALL cp_gemm('N', 'N', nao, nmoloc, ncol, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, & 343 a_first_col=isgf, a_first_row=1, b_first_col=1, b_first_row=isgf) 344 345 CALL cp_gemm('T', 'N', nmoloc, nmoloc, nao, 0.5_dp, vectors, opvec, & 346 1.0_dp, zij_fm_set(iatom, 1)%matrix, & 347 a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1) 348 349 END DO ! iatom 350 351 ! And now perform the optimization and rotate the orbitals 352 SELECT CASE (method) 353 CASE (do_loc_jacobi) 354 CALL jacobi_rotation_pipek(zij_fm_set, vectors, sweeps) 355 CASE (do_loc_crazy) 356 CPABORT("Crazy and Pipek not implemented.") 357 CASE (do_loc_l1_norm_sd) 358 CPABORT("L1 norm and Pipek not implemented.") 359 CASE (do_loc_direct) 360 CPABORT("Direct and Pipek not implemented.") 361 CASE (do_loc_none) 362 IF (output_unit > 0) WRITE (output_unit, '(A,I6,A)') " No MOS localization applied " 363 CASE DEFAULT 364 CPABORT("Unknown localization method") 365 END SELECT 366 367 IF (output_unit > 0) WRITE (output_unit, '(A,I3,A,I6,A)') " Localization for spin ", ispin, & 368 " converged in ", sweeps, " iterations" 369 370 CALL centers_spreads_pipek(qs_loc_env, zij_fm_set, particle_set, ispin, & 371 print_loc_section) 372 373 DEALLOCATE (first_sgf, last_sgf, nsgf) 374 375 CALL cp_fm_release(opvec) 376 CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO") 377 378 CALL timestop(handle) 379 380 END SUBROUTINE optimize_loc_pipek 381 382! ************************************************************************************************** 383!> \brief ... 384!> \param qs_loc_env ... 385!> \param zij ... 386!> \param nmoloc ... 387!> \param cell ... 388!> \param weights ... 389!> \param ispin ... 390!> \param print_loc_section ... 391!> \param only_initial_out ... 392!> \par History 393!> 04.2005 created [MI] 394!> \author MI 395! ************************************************************************************************** 396 SUBROUTINE centers_spreads_berry(qs_loc_env, zij, nmoloc, cell, weights, ispin, & 397 print_loc_section, only_initial_out) 398 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 399 TYPE(cp_fm_p_type), INTENT(INOUT) :: zij(:, :) 400 INTEGER, INTENT(IN) :: nmoloc 401 TYPE(cell_type), POINTER :: cell 402 REAL(dp), DIMENSION(:) :: weights 403 INTEGER, INTENT(IN) :: ispin 404 TYPE(section_vals_type), POINTER :: print_loc_section 405 LOGICAL, INTENT(IN), OPTIONAL :: only_initial_out 406 407 CHARACTER(len=*), PARAMETER :: routineN = 'centers_spreads_berry', & 408 routineP = moduleN//':'//routineN 409 410 CHARACTER(len=default_path_length) :: file_tmp, iter 411 COMPLEX(KIND=dp) :: z 412 INTEGER :: idir, istate, jdir, nstates, & 413 output_unit, unit_out_s 414 LOGICAL :: my_only_init 415 REAL(dp) :: spread_i, spread_ii, sum_spread_i, & 416 sum_spread_ii 417 REAL(dp), DIMENSION(3) :: c, c2, cpbc 418 REAL(dp), DIMENSION(:, :), POINTER :: centers 419 REAL(KIND=dp) :: imagpart, realpart 420 TYPE(cp_logger_type), POINTER :: logger 421 TYPE(section_vals_type), POINTER :: print_key 422 423 NULLIFY (centers, logger, print_key) 424 logger => cp_get_default_logger() 425 my_only_init = .FALSE. 426 IF (PRESENT(only_initial_out)) my_only_init = only_initial_out 427 428 file_tmp = TRIM(qs_loc_env%tag_mo)//"_spreads_s"//TRIM(ADJUSTL(cp_to_string(ispin))) 429 output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", & 430 extension=".locInfo") 431 unit_out_s = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", & 432 middle_name=file_tmp, extension=".data") 433 iter = cp_iter_string(logger%iter_info) 434 IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(i6,/,A)') nmoloc, TRIM(iter) 435 436 CALL cp_fm_get_info(zij(1, 1)%matrix, nrow_global=nstates) 437 CPASSERT(nstates >= nmoloc) 438 439 centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array 440 CPASSERT(ASSOCIATED(centers)) 441 CPASSERT(SIZE(centers, 2) == nmoloc) 442 sum_spread_i = 0.0_dp 443 sum_spread_ii = 0.0_dp 444 DO istate = 1, nmoloc 445 c = 0.0_dp 446 c2 = 0.0_dp 447 spread_i = 0.0_dp 448 spread_ii = 0.0_dp 449 DO jdir = 1, SIZE(zij, 2) 450 CALL cp_fm_get_element(zij(1, jdir)%matrix, istate, istate, realpart) 451 CALL cp_fm_get_element(zij(2, jdir)%matrix, istate, istate, imagpart) 452 z = CMPLX(realpart, imagpart, dp) 453 spread_i = spread_i - weights(jdir)* & 454 LOG(realpart*realpart + imagpart*imagpart)/twopi/twopi 455 spread_ii = spread_ii + weights(jdir)* & 456 (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi 457 IF (jdir < 4) THEN 458 DO idir = 1, 3 459 c(idir) = c(idir) + & 460 (cell%hmat(idir, jdir)/twopi)*AIMAG(LOG(z)) 461 ENDDO 462 END IF 463 END DO 464 cpbc = pbc(c, cell) 465 centers(1:3, istate) = cpbc(1:3) 466 centers(4, istate) = spread_i 467 centers(5, istate) = spread_ii 468 sum_spread_i = sum_spread_i + centers(4, istate) 469 sum_spread_ii = sum_spread_ii + centers(5, istate) 470 IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(I6,2F16.8)') istate, centers(4:5, istate) 471 ENDDO 472 473 ! Print of wannier centers 474 print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS") 475 IF (.NOT. my_only_init) CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin) 476 477 IF (output_unit > 0) THEN 478 WRITE (output_unit, '(T4, A, 2x, A26, A26)') " Spread Functional ", "sum_in -w_i ln(|z_in|^2)", & 479 "sum_in w_i(1-|z_in|^2)" 480 IF (my_only_init) THEN 481 WRITE (output_unit, '(T4,A,T38,2F20.10)') " Initial Spread (Berry) : ", sum_spread_i, sum_spread_ii 482 ELSE 483 WRITE (output_unit, '(T4,A,T38,2F20.10)') " Total Spread (Berry) : ", sum_spread_i, sum_spread_ii 484 END IF 485 END IF 486 487 IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(A,2F16.10)') " Total ", sum_spread_i, sum_spread_ii 488 489 CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS") 490 CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO") 491 492 END SUBROUTINE centers_spreads_berry 493 494! ************************************************************************************************** 495!> \brief define and print the centers and spread 496!> when the pipek operator is used 497!> \param qs_loc_env ... 498!> \param zij_fm_set matrix elements that define the populations on atoms 499!> \param particle_set ... 500!> \param ispin spin 1 or 2 501!> \param print_loc_section ... 502!> \par History 503!> 05.2005 created [MI] 504! ************************************************************************************************** 505 SUBROUTINE centers_spreads_pipek(qs_loc_env, zij_fm_set, particle_set, ispin, & 506 print_loc_section) 507 508 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 509 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: zij_fm_set 510 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 511 INTEGER, INTENT(IN) :: ispin 512 TYPE(section_vals_type), POINTER :: print_loc_section 513 514 CHARACTER(len=*), PARAMETER :: routineN = 'centers_spreads_pipek', & 515 routineP = moduleN//':'//routineN 516 517 CHARACTER(len=default_path_length) :: file_tmp, iter 518 INTEGER :: iatom, istate, natom, nstate, unit_out_s 519 INTEGER, DIMENSION(:), POINTER :: atom_of_state 520 REAL(dp) :: r(3) 521 REAL(dp), ALLOCATABLE, DIMENSION(:) :: Qii, ziimax 522 REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: diag 523 REAL(dp), DIMENSION(:, :), POINTER :: centers 524 TYPE(cp_logger_type), POINTER :: logger 525 TYPE(section_vals_type), POINTER :: print_key 526 527 NULLIFY (centers, logger, print_key) 528 logger => cp_get_default_logger() 529 530 CALL cp_fm_get_info(zij_fm_set(1, 1)%matrix, nrow_global=nstate) 531 natom = SIZE(zij_fm_set, 1) 532 533 centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array 534 CPASSERT(ASSOCIATED(centers)) 535 CPASSERT(SIZE(centers, 2) == nstate) 536 537 file_tmp = TRIM(qs_loc_env%tag_mo)//"_spreads_s"//TRIM(ADJUSTL(cp_to_string(ispin))) 538 unit_out_s = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", & 539 middle_name=file_tmp, extension=".data", log_filename=.FALSE.) 540 iter = cp_iter_string(logger%iter_info) 541 IF (unit_out_s > 0) WRITE (unit_out_s, '(i6,/,A)') TRIM(iter) 542 543 ALLOCATE (atom_of_state(nstate)) 544 atom_of_state = 0 545 546 ALLOCATE (diag(nstate, natom)) 547 diag = 0.0_dp 548 549 DO iatom = 1, natom 550 DO istate = 1, nstate 551 CALL cp_fm_get_element(zij_fm_set(iatom, 1)%matrix, istate, istate, diag(istate, iatom)) 552 END DO 553 END DO 554 555 ALLOCATE (Qii(nstate), ziimax(nstate)) 556 ziimax = 0.0_dp 557 Qii = 0.0_dp 558 559 DO iatom = 1, natom 560 DO istate = 1, nstate 561 Qii(istate) = Qii(istate) + diag(istate, iatom)*diag(istate, iatom) 562 IF (ABS(diag(istate, iatom)) > ziimax(istate)) THEN 563 ziimax(istate) = ABS(diag(istate, iatom)) 564 atom_of_state(istate) = iatom 565 END IF 566 END DO 567 END DO 568 569 DO istate = 1, nstate 570 iatom = atom_of_state(istate) 571 r(1:3) = particle_set(iatom)%r(1:3) 572 centers(1:3, istate) = r(1:3) 573 centers(4, istate) = 1.0_dp/Qii(istate) 574 IF (unit_out_s > 0) WRITE (unit_out_s, '(I6,F16.8)') istate, angstrom*centers(4, istate) 575 END DO 576 577 ! Print the wannier centers 578 print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS") 579 CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin) 580 581 CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS") 582 583 DEALLOCATE (Qii, ziimax, atom_of_state, diag) 584 585 END SUBROUTINE centers_spreads_pipek 586 587! ************************************************************************************************** 588!> \brief set up the calculation of localized orbitals 589!> \param qs_env ... 590!> \param qs_loc_env ... 591!> \param print_loc_section ... 592!> \param myspin ... 593!> \param ext_mo_coeff ... 594!> \par History 595!> 04.2005 created [MI] 596!> \author MI 597! ************************************************************************************************** 598 SUBROUTINE qs_loc_driver(qs_env, qs_loc_env, print_loc_section, myspin, & 599 ext_mo_coeff) 600 601 TYPE(qs_environment_type), POINTER :: qs_env 602 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 603 TYPE(section_vals_type), POINTER :: print_loc_section 604 INTEGER, INTENT(IN), OPTIONAL :: myspin 605 TYPE(cp_fm_type), OPTIONAL, POINTER :: ext_mo_coeff 606 607 CHARACTER(len=*), PARAMETER :: routineN = 'qs_loc_driver', routineP = moduleN//':'//routineN 608 609 CHARACTER(LEN=default_string_length) :: my_pos 610 INTEGER :: dim_op, handle, i, imo, imoloc, ir, & 611 ispin, istate, j, jstate, l_spin, lb, & 612 loc_method, n_rep, nao, ncubes, nmo, & 613 nmosub, s_spin, ub 614 INTEGER, DIMENSION(:), POINTER :: bounds, list, list_cubes 615 LOGICAL :: append_cube, list_cubes_setup 616 LOGICAL, SAVE :: first_time = .TRUE. 617 REAL(dp), DIMENSION(6) :: weights 618 REAL(KIND=dp), DIMENSION(:, :), POINTER :: centers, vecbuffer 619 TYPE(cell_type), POINTER :: cell 620 TYPE(cp_fm_p_type), DIMENSION(:), POINTER :: moloc_coeff 621 TYPE(cp_fm_p_type), DIMENSION(:, :), POINTER :: op_fm_set 622 TYPE(cp_fm_struct_type), POINTER :: tmp_fm_struct 623 TYPE(cp_fm_type), POINTER :: mo_coeff 624 TYPE(cp_para_env_type), POINTER :: para_env 625 TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set 626 TYPE(dft_control_type), POINTER :: dft_control 627 TYPE(localized_wfn_control_type), POINTER :: localized_wfn_control 628 TYPE(mo_set_p_type), DIMENSION(:), POINTER :: mos 629 TYPE(section_vals_type), POINTER :: print_key 630 631 CALL timeset(routineN, handle) 632 NULLIFY (para_env, mos, dft_control, list_cubes) 633 NULLIFY (cell, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set) 634 qs_loc_env%first_time = first_time 635 qs_loc_env%target_time = qs_env%target_time 636 qs_loc_env%start_time = qs_env%start_time 637 638 CALL get_qs_loc_env(qs_loc_env=qs_loc_env, & 639 localized_wfn_control=localized_wfn_control, & 640 moloc_coeff=moloc_coeff, op_sm_set=op_sm_set, op_fm_set=op_fm_set, cell=cell, & 641 weights=weights, dim_op=dim_op) 642 643 CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, & 644 para_env=para_env, mos=mos) 645 646 s_spin = 1 647 l_spin = dft_control%nspins 648 IF (PRESENT(myspin)) THEN 649 s_spin = myspin 650 l_spin = myspin 651 END IF 652 653 DO ispin = s_spin, l_spin 654 655 CALL get_mo_set(mo_set=mos(ispin)%mo_set, nao=nao, nmo=nmo) 656 loc_method = localized_wfn_control%localization_method 657 658 SELECT CASE (localized_wfn_control%operator_type) 659 660 CASE (op_loc_berry) 661 ! Here we allocate op_fm_set with the RIGHT size for uks 662 NULLIFY (tmp_fm_struct, mo_coeff) 663 nmosub = localized_wfn_control%nloc_states(ispin) 664 IF (PRESENT(ext_mo_coeff)) THEN 665 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmosub, & 666 ncol_global=nmosub, para_env=para_env, context=ext_mo_coeff%matrix_struct%context) 667 ELSE 668 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 669 CALL cp_fm_struct_create(tmp_fm_struct, nrow_global=nmosub, & 670 ncol_global=nmosub, para_env=para_env, context=mo_coeff%matrix_struct%context) 671 END IF 672 ! 673 ALLOCATE (op_fm_set(2, dim_op)) 674 DO i = 1, dim_op 675 DO j = 1, SIZE(op_fm_set, 1) 676 NULLIFY (op_fm_set(j, i)%matrix) 677 CALL cp_fm_create(op_fm_set(j, i)%matrix, tmp_fm_struct) 678 CALL cp_fm_get_info(op_fm_set(j, i)%matrix, nrow_global=nmosub) 679 CPASSERT(nmo >= nmosub) 680 CALL cp_fm_set_all(op_fm_set(j, i)%matrix, 0.0_dp) 681 END DO 682 END DO 683 CALL cp_fm_struct_release(tmp_fm_struct) 684 685 CALL optimize_loc_berry(loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, op_sm_set, & 686 op_fm_set, para_env, cell, weights, ispin, print_loc_section) 687 688 ! Here we dealloctate op_fm_set 689 IF (ASSOCIATED(op_fm_set)) THEN 690 DO i = 1, dim_op 691 DO j = 1, SIZE(op_fm_set, 1) 692 CALL cp_fm_release(op_fm_set(j, i)%matrix) 693 END DO 694 END DO 695 DEALLOCATE (op_fm_set) 696 END IF 697 698 CASE (op_loc_boys) 699 CPABORT("Boys localization not implemented") 700 701 CASE (op_loc_pipek) 702 703 CALL optimize_loc_pipek(qs_env, loc_method, qs_loc_env, moloc_coeff(ispin)%matrix, & 704 op_fm_set, ispin, print_loc_section) 705 706 END SELECT 707 708 ! give back the localized orbitals 709 IF (.NOT. PRESENT(ext_mo_coeff)) THEN 710 CALL get_mo_set(mo_set=mos(ispin)%mo_set, mo_coeff=mo_coeff) 711 END IF 712 713 lb = localized_wfn_control%lu_bound_states(1, ispin) 714 ub = localized_wfn_control%lu_bound_states(2, ispin) 715 716 IF (localized_wfn_control%set_of_states == state_loc_list) THEN 717 ALLOCATE (vecbuffer(1, nao)) 718 nmosub = SIZE(localized_wfn_control%loc_states, 1) 719 imoloc = 0 720 DO i = lb, ub 721 ! Get the index in the subset 722 imoloc = imoloc + 1 723 ! Get the index in the full set 724 imo = localized_wfn_control%loc_states(i, ispin) 725 726 CALL cp_fm_get_submatrix(moloc_coeff(ispin)%matrix, vecbuffer, 1, imoloc, & 727 nao, 1, transpose=.TRUE.) 728 IF (PRESENT(ext_mo_coeff)) THEN 729 CALL cp_fm_set_submatrix(ext_mo_coeff, vecbuffer, 1, imo, & 730 nao, 1, transpose=.TRUE.) 731 ELSE 732 CALL cp_fm_set_submatrix(mo_coeff, vecbuffer, 1, imo, & 733 nao, 1, transpose=.TRUE.) 734 END IF 735 END DO 736 DEALLOCATE (vecbuffer) 737 738 ELSE 739 nmosub = localized_wfn_control%nloc_states(ispin) 740 IF (PRESENT(ext_mo_coeff)) THEN 741 CALL cp_fm_to_fm(moloc_coeff(ispin)%matrix, ext_mo_coeff, nmosub, 1, lb) 742 ELSE 743 CALL cp_fm_to_fm(moloc_coeff(ispin)%matrix, mo_coeff, nmosub, 1, lb) 744 END IF 745 END IF 746 747 ! Write cube files if required 748 IF (localized_wfn_control%print_cubes) THEN 749 list_cubes_setup = .FALSE. 750 NULLIFY (bounds, list, list_cubes) 751 752 ! Provides boundaries of MOs 753 CALL section_vals_val_get(print_loc_section, "WANNIER_CUBES%CUBES_LU_BOUNDS", & 754 i_vals=bounds) 755 ncubes = bounds(2) - bounds(1) + 1 756 IF (ncubes > 0) THEN 757 list_cubes_setup = .TRUE. 758 ALLOCATE (list_cubes(ncubes)) 759 DO ir = 1, ncubes 760 list_cubes(ir) = bounds(1) + (ir - 1) 761 END DO 762 END IF 763 764 ! Provides the list of MOs 765 CALL section_vals_val_get(print_loc_section, "WANNIER_CUBES%CUBES_LIST", & 766 n_rep_val=n_rep) 767 IF (.NOT. list_cubes_setup) THEN 768 ncubes = 0 769 DO ir = 1, n_rep 770 CALL section_vals_val_get(print_loc_section, "WANNIER_CUBES%CUBES_LIST", & 771 i_rep_val=ir, i_vals=list) 772 IF (ASSOCIATED(list)) THEN 773 CALL reallocate(list_cubes, 1, ncubes + SIZE(list)) 774 DO i = 1, SIZE(list) 775 list_cubes(i + ncubes) = list(i) 776 END DO 777 ncubes = ncubes + SIZE(list) 778 END IF 779 END DO 780 IF (ncubes > 0) list_cubes_setup = .TRUE. 781 END IF 782 783 ! Full list of Mos 784 IF (.NOT. list_cubes_setup) THEN 785 list_cubes_setup = .TRUE. 786 ncubes = localized_wfn_control%nloc_states(1) 787 IF (ncubes > 0) THEN 788 ALLOCATE (list_cubes(ncubes)) 789 END IF 790 DO i = 1, ncubes 791 list_cubes(i) = i 792 END DO 793 END IF 794 795 ncubes = SIZE(list_cubes) 796 ncubes = MIN(ncubes, nmo) 797 ALLOCATE (centers(6, ncubes)) 798 DO i = 1, ncubes 799 istate = list_cubes(i) 800 DO j = 1, localized_wfn_control%nloc_states(ispin) 801 jstate = localized_wfn_control%loc_states(j, ispin) 802 IF (istate == jstate) THEN 803 centers(1:6, i) = localized_wfn_control%centers_set(ispin)%array(1:6, j) 804 EXIT 805 ENDIF 806 END DO 807 END DO ! ncubes 808 809 ! Real call for dumping the cube files 810 print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CUBES") 811 append_cube = section_get_lval(print_loc_section, "WANNIER_CUBES%APPEND") 812 my_pos = "REWIND" 813 IF (append_cube) THEN 814 my_pos = "APPEND" 815 END IF 816 817 CALL qs_print_cubes(qs_env, moloc_coeff(ispin)%matrix, ncubes, list_cubes, centers, & 818 print_key, "loc"//TRIM(ADJUSTL(qs_loc_env%tag_mo)), & 819 ispin=ispin, file_position=my_pos) 820 821 DEALLOCATE (centers) 822 DEALLOCATE (list_cubes) 823 END IF 824 END DO ! ispin 825 first_time = .FALSE. 826 CALL timestop(handle) 827 END SUBROUTINE qs_loc_driver 828 829! ************************************************************************************************** 830!> \brief write the cube files for a set of selected states 831!> \param qs_env ... 832!> \param mo_coeff set mos from which the states to be printed are extracted 833!> \param nstates number of states to be printed 834!> \param state_list list of the indexes of the states to be printed 835!> \param centers centers and spread, all=0 if they hva not been calculated 836!> \param print_key ... 837!> \param root initial part of the cube file names 838!> \param ispin ... 839!> \param idir ... 840!> \param state0 ... 841!> \param file_position ... 842!> \par History 843!> 08.2005 created [MI] 844!> \author MI 845!> \note 846!> This routine should not be in this module 847! ************************************************************************************************** 848 SUBROUTINE qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, & 849 print_key, root, ispin, idir, state0, file_position) 850 851 TYPE(qs_environment_type), POINTER :: qs_env 852 TYPE(cp_fm_type), POINTER :: mo_coeff 853 INTEGER, INTENT(IN) :: nstates 854 INTEGER, DIMENSION(:), POINTER :: state_list 855 REAL(dp), DIMENSION(:, :), POINTER :: centers 856 TYPE(section_vals_type), POINTER :: print_key 857 CHARACTER(LEN=*) :: root 858 INTEGER, INTENT(IN), OPTIONAL :: ispin, idir 859 INTEGER, OPTIONAL :: state0 860 CHARACTER(LEN=default_string_length), OPTIONAL :: file_position 861 862 CHARACTER(len=*), PARAMETER :: routineN = 'qs_print_cubes', routineP = moduleN//':'//routineN 863 CHARACTER, DIMENSION(3), PARAMETER :: labels = (/'x', 'y', 'z'/) 864 865 CHARACTER :: label 866 CHARACTER(LEN=default_path_length) :: file_tmp, filename, my_pos 867 CHARACTER(LEN=default_string_length) :: title 868 INTEGER :: handle, ia, ie, istate, ivector, & 869 my_ispin, my_state0, unit_out_c 870 LOGICAL :: add_idir, add_spin, mpi_io 871 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 872 TYPE(cell_type), POINTER :: cell 873 TYPE(cp_logger_type), POINTER :: logger 874 TYPE(dft_control_type), POINTER :: dft_control 875 TYPE(particle_list_type), POINTER :: particles 876 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 877 TYPE(pw_env_type), POINTER :: pw_env 878 TYPE(pw_p_type) :: wf_g, wf_r 879 TYPE(pw_pool_type), POINTER :: auxbas_pw_pool 880 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 881 TYPE(qs_subsys_type), POINTER :: subsys 882 883 CALL timeset(routineN, handle) 884 NULLIFY (auxbas_pw_pool, pw_env, logger) 885 logger => cp_get_default_logger() 886 887 CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys) 888 CALL qs_subsys_get(subsys, particles=particles) 889 890 CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool) 891 892 CALL pw_pool_create_pw(auxbas_pw_pool, wf_r%pw, & 893 use_data=REALDATA3D, & 894 in_space=REALSPACE) 895 CALL pw_pool_create_pw(auxbas_pw_pool, wf_g%pw, & 896 use_data=COMPLEXDATA1D, & 897 in_space=RECIPROCALSPACE) 898 899 my_state0 = 0 900 IF (PRESENT(state0)) my_state0 = state0 901 902 add_spin = .FALSE. 903 my_ispin = 1 904 IF (PRESENT(ispin)) THEN 905 add_spin = .TRUE. 906 my_ispin = ispin 907 END IF 908 add_idir = .FALSE. 909 IF (PRESENT(idir)) THEN 910 add_idir = .TRUE. 911 label = labels(idir) 912 END IF 913 914 my_pos = "REWIND" 915 IF (PRESENT(file_position)) THEN 916 my_pos = file_position 917 END IF 918 919 DO istate = 1, nstates 920 ivector = state_list(istate) - my_state0 921 CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, & 922 dft_control=dft_control, particle_set=particle_set, pw_env=pw_env) 923 924 CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, & 925 qs_kind_set, cell, dft_control, particle_set, pw_env) 926 927 ! Formatting the middle part of the name 928 ivector = state_list(istate) 929 CALL xstring(root, ia, ie) 930 IF (add_idir) THEN 931 filename = root(ia:ie)//"_"//label//"_w"//TRIM(ADJUSTL(cp_to_string(ivector))) 932 ELSE 933 filename = root(ia:ie)//"_w"//TRIM(ADJUSTL(cp_to_string(ivector))) 934 END IF 935 IF (add_spin) THEN 936 file_tmp = filename 937 CALL xstring(file_tmp, ia, ie) 938 filename = file_tmp(ia:ie)//"_s"//TRIM(ADJUSTL(cp_to_string(ispin))) 939 END IF 940 941 ! Using the print_key tools to open the file with the proper name 942 mpi_io = .TRUE. 943 unit_out_c = cp_print_key_unit_nr(logger, print_key, "", middle_name=filename, & 944 extension=".cube", file_position=my_pos, log_filename=.FALSE., & 945 mpi_io=mpi_io) 946 IF (SIZE(centers, 1) == 6) THEN 947 WRITE (title, '(A7,I5.5,A2,I1.1,A1,6F10.4)') "WFN ", ivector, "_s", my_ispin, " ", & 948 centers(1:3, istate)*angstrom, centers(4:6, istate)*angstrom 949 ELSE 950 WRITE (title, '(A7,I5.5,A2,I1.1,A1,3F10.4)') "WFN ", ivector, "_s", my_ispin, " ", & 951 centers(1:3, istate)*angstrom 952 END IF 953 CALL cp_pw_to_cube(wf_r%pw, unit_out_c, title, & 954 particles=particles, & 955 stride=section_get_ivals(print_key, "STRIDE"), mpi_io=mpi_io) 956 CALL cp_print_key_finished_output(unit_out_c, logger, print_key, "", mpi_io=mpi_io) 957 END DO 958 959 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_r%pw) 960 CALL pw_pool_give_back_pw(auxbas_pw_pool, wf_g%pw) 961 CALL timestop(handle) 962 END SUBROUTINE qs_print_cubes 963 964! ************************************************************************************************** 965!> \brief Prints wannier centers 966!> \param qs_loc_env ... 967!> \param print_key ... 968!> \param center ... 969!> \param logger ... 970!> \param ispin ... 971! ************************************************************************************************** 972 SUBROUTINE print_wannier_centers(qs_loc_env, print_key, center, logger, ispin) 973 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 974 TYPE(section_vals_type), POINTER :: print_key 975 REAL(KIND=dp), INTENT(IN) :: center(:, :) 976 TYPE(cp_logger_type), POINTER :: logger 977 INTEGER, INTENT(IN) :: ispin 978 979 CHARACTER(len=*), PARAMETER :: routineN = 'print_wannier_centers', & 980 routineP = moduleN//':'//routineN 981 982 CHARACTER(default_string_length) :: iter, unit_str 983 CHARACTER(LEN=default_string_length) :: my_ext, my_form 984 INTEGER :: iunit, l, nstates 985 LOGICAL :: first_time, init_traj 986 REAL(KIND=dp) :: unit_conv 987 988 nstates = SIZE(center, 2) 989 my_form = "formatted" 990 my_ext = ".data" 991 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, first_time=first_time) & 992 , cp_p_file)) THEN 993 ! Find out if we want to print IONS+CENTERS or ONLY CENTERS.. 994 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS") & 995 , cp_p_file)) THEN 996 CALL get_output_format(print_key, my_form=my_form, my_ext=my_ext) 997 END IF 998 IF (first_time .OR. (.NOT. qs_loc_env%first_time)) THEN 999 iunit = cp_print_key_unit_nr(logger, print_key, "", extension=my_ext, file_form=my_form, & 1000 middle_name=TRIM(qs_loc_env%tag_mo)//"_centers_s"//TRIM(ADJUSTL(cp_to_string(ispin))), & 1001 log_filename=.FALSE., on_file=.TRUE., is_new_file=init_traj) 1002 IF (iunit > 0) THEN 1003 ! Gather units of measure for output (if available) 1004 CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str) 1005 unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 1006 1007 IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS"), cp_p_file)) THEN 1008 ! Different possible formats 1009 CALL print_wannier_traj(qs_loc_env, print_key, center, iunit, init_traj, unit_conv) 1010 ELSE 1011 ! Default print format 1012 iter = cp_iter_string(logger%iter_info) 1013 WRITE (iunit, '(i6,/,a)') nstates, TRIM(iter) 1014 DO l = 1, nstates 1015 WRITE (iunit, '(A,4F16.8)') "X", unit_conv*center(1:4, l) 1016 END DO 1017 END IF 1018 END IF 1019 CALL cp_print_key_finished_output(iunit, logger, print_key, on_file=.TRUE.) 1020 END IF 1021 END IF 1022 END SUBROUTINE print_wannier_centers 1023 1024! ************************************************************************************************** 1025!> \brief computes spread and centers 1026!> \param qs_loc_env ... 1027!> \param print_key ... 1028!> \param center ... 1029!> \param iunit ... 1030!> \param init_traj ... 1031!> \param unit_conv ... 1032! ************************************************************************************************** 1033 SUBROUTINE print_wannier_traj(qs_loc_env, print_key, center, iunit, init_traj, unit_conv) 1034 TYPE(qs_loc_env_new_type), POINTER :: qs_loc_env 1035 TYPE(section_vals_type), POINTER :: print_key 1036 REAL(KIND=dp), INTENT(IN) :: center(:, :) 1037 INTEGER, INTENT(IN) :: iunit 1038 LOGICAL, INTENT(IN) :: init_traj 1039 REAL(KIND=dp), INTENT(IN) :: unit_conv 1040 1041 CHARACTER(len=*), PARAMETER :: routineN = 'print_wannier_traj', & 1042 routineP = moduleN//':'//routineN 1043 1044 CHARACTER(len=default_string_length) :: iter, remark1, remark2, title 1045 INTEGER :: i, iskip, natom, ntot, outformat 1046 TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set 1047 TYPE(atomic_kind_type), POINTER :: atomic_kind 1048 TYPE(cp_logger_type), POINTER :: logger 1049 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 1050 1051 NULLIFY (particle_set, atomic_kind_set, atomic_kind, logger) 1052 logger => cp_get_default_logger() 1053 natom = SIZE(qs_loc_env%particle_set) 1054 ntot = natom + SIZE(center, 2) 1055 CALL allocate_particle_set(particle_set, ntot) 1056 ALLOCATE (atomic_kind_set(1)) 1057 atomic_kind => atomic_kind_set(1) 1058 CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=0, & 1059 name="X", element_symbol="X", mass=0.0_dp) 1060 ! Particles 1061 DO i = 1, natom 1062 particle_set(i)%atomic_kind => qs_loc_env%particle_set(i)%atomic_kind 1063 particle_set(i)%r = pbc(qs_loc_env%particle_set(i)%r, qs_loc_env%cell) 1064 END DO 1065 ! Wannier Centers 1066 DO i = natom + 1, ntot 1067 particle_set(i)%atomic_kind => atomic_kind 1068 particle_set(i)%r = pbc(center(1:3, i - natom), qs_loc_env%cell) 1069 END DO 1070 ! Dump the structure 1071 CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat) 1072 1073 ! Header file 1074 SELECT CASE (outformat) 1075 CASE (dump_dcd, dump_dcd_aligned_cell) 1076 IF (init_traj) THEN 1077 !Lets write the header for the coordinate dcd 1078 ! Note (TL) : even the new DCD format is unfortunately too poor 1079 ! for our capabilities.. for example here the printing 1080 ! of the geometry could be nested inside several iteration 1081 ! levels.. this cannot be exactly reproduce with DCD. 1082 ! Just as a compromise let's pick-up the value of the MD iteration 1083 ! level. In any case this is not any sensible information for the standard.. 1084 iskip = section_get_ival(print_key, "EACH%MD") 1085 WRITE (iunit) "CORD", 0, -1, iskip, & 1086 0, 0, 0, 0, 0, 0, REAL(0, KIND=sp), 1, 0, 0, 0, 0, 0, 0, 0, 0, 24 1087 remark1 = "REMARK FILETYPE CORD DCD GENERATED BY CP2K" 1088 remark2 = "REMARK Support new DCD format with cell information" 1089 WRITE (iunit) 2, remark1, remark2 1090 WRITE (iunit) ntot 1091 CALL m_flush(iunit) 1092 END IF 1093 CASE (dump_xmol) 1094 iter = cp_iter_string(logger%iter_info) 1095 WRITE (UNIT=title, FMT="(A)") " Particles+Wannier centers. Iteration:"//TRIM(iter) 1096 CASE DEFAULT 1097 title = "" 1098 END SELECT 1099 CALL write_particle_coordinates(particle_set, iunit, outformat, "POS", title, qs_loc_env%cell, & 1100 unit_conv=unit_conv) 1101 CALL m_flush(iunit) 1102 CALL deallocate_particle_set(particle_set) 1103 CALL deallocate_atomic_kind_set(atomic_kind_set) 1104 END SUBROUTINE print_wannier_traj 1105 1106END MODULE qs_loc_methods 1107