1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2020 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Define methods related to particle_type 8!> \par History 9!> 10.2014 Move routines out of particle_types.F [Ole Schuett] 10!> \author Ole Schuett 11! ************************************************************************************************** 12MODULE particle_methods 13 USE atomic_kind_types, ONLY: get_atomic_kind 14 USE basis_set_types, ONLY: get_gto_basis_set,& 15 gto_basis_set_p_type 16 USE cell_types, ONLY: cell_clone,& 17 cell_create,& 18 cell_release,& 19 cell_type,& 20 get_cell,& 21 pbc,& 22 real_to_scaled,& 23 set_cell_param 24 USE cp_log_handling, ONLY: cp_get_default_logger,& 25 cp_logger_type 26 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 27 cp_print_key_unit_nr 28 USE cp_units, ONLY: cp_unit_from_cp2k 29 USE external_potential_types, ONLY: fist_potential_type,& 30 get_potential 31 USE input_constants, ONLY: dump_atomic,& 32 dump_dcd,& 33 dump_dcd_aligned_cell,& 34 dump_pdb,& 35 dump_xmol 36 USE input_section_types, ONLY: section_vals_get_subs_vals,& 37 section_vals_type,& 38 section_vals_val_get 39 USE kinds, ONLY: default_string_length,& 40 dp,& 41 sp 42 USE mathconstants, ONLY: degree 43 USE mathlib, ONLY: angle,& 44 dihedral_angle 45 USE memory_utilities, ONLY: reallocate 46 USE particle_types, ONLY: get_particle_pos_or_vel,& 47 particle_type 48 USE physcon, ONLY: massunit 49 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm 50 USE qs_kind_types, ONLY: get_qs_kind,& 51 qs_kind_type 52 USE shell_potential_types, ONLY: get_shell,& 53 shell_kind_type 54 USE string_utilities, ONLY: uppercase 55 USE util, ONLY: sort,& 56 sort_unique 57#include "./base/base_uses.f90" 58 59 IMPLICIT NONE 60 61 PRIVATE 62 63 ! Public subroutines 64 65 PUBLIC :: write_fist_particle_coordinates, & 66 write_qs_particle_coordinates, & 67 write_particle_distances, & 68 write_particle_coordinates, & 69 write_structure_data, & 70 get_particle_set, & 71 write_particle_matrix 72 73 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'particle_methods' 74 75CONTAINS 76 77! ************************************************************************************************** 78!> \brief Get the components of a particle set. 79!> \param particle_set ... 80!> \param qs_kind_set ... 81!> \param first_sgf ... 82!> \param last_sgf ... 83!> \param nsgf ... 84!> \param nmao ... 85!> \param basis ... 86!> \date 14.01.2002 87!> \par History 88!> - particle type cleaned (13.10.2003,MK) 89!> - refactoring and add basis set option (17.08.2010,jhu) 90!> \author MK 91!> \version 1.0 92! ************************************************************************************************** 93 SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, & 94 nmao, basis) 95 96 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 97 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 98 INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: first_sgf, last_sgf, nsgf, nmao 99 TYPE(gto_basis_set_p_type), DIMENSION(:), OPTIONAL :: basis 100 101 CHARACTER(len=*), PARAMETER :: routineN = 'get_particle_set', & 102 routineP = moduleN//':'//routineN 103 104 INTEGER :: ikind, iparticle, isgf, nparticle, ns 105 106 CPASSERT(ASSOCIATED(particle_set)) 107 108 nparticle = SIZE(particle_set) 109 IF (PRESENT(first_sgf)) THEN 110 CPASSERT(SIZE(first_sgf) >= nparticle) 111 END IF 112 IF (PRESENT(last_sgf)) THEN 113 CPASSERT(SIZE(last_sgf) >= nparticle) 114 END IF 115 IF (PRESENT(nsgf)) THEN 116 CPASSERT(SIZE(nsgf) >= nparticle) 117 END IF 118 IF (PRESENT(nmao)) THEN 119 CPASSERT(SIZE(nmao) >= nparticle) 120 END IF 121 122 IF (PRESENT(first_sgf) .OR. PRESENT(last_sgf) .OR. PRESENT(nsgf)) THEN 123 isgf = 0 124 DO iparticle = 1, nparticle 125 CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind) 126 IF (PRESENT(basis)) THEN 127 CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns) 128 ELSE 129 CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns) 130 END IF 131 IF (PRESENT(nsgf)) nsgf(iparticle) = ns 132 IF (PRESENT(first_sgf)) first_sgf(iparticle) = isgf + 1 133 isgf = isgf + ns 134 IF (PRESENT(last_sgf)) last_sgf(iparticle) = isgf 135 END DO 136 END IF 137 138 IF (PRESENT(first_sgf)) THEN 139 IF (SIZE(first_sgf) > nparticle) first_sgf(nparticle + 1) = isgf + 1 140 END IF 141 142 IF (PRESENT(nmao)) THEN 143 DO iparticle = 1, nparticle 144 CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind) 145 CALL get_qs_kind(qs_kind_set(ikind), mao=ns) 146 nmao(iparticle) = ns 147 END DO 148 END IF 149 150 END SUBROUTINE get_particle_set 151 152! ************************************************************************************************** 153!> \brief Should be able to write a few formats e.g. xmol, and some binary 154!> format (dcd) some format can be used for x,v,f 155!> 156!> FORMAT CONTENT UNITS x, v, f 157!> XMOL POS, VEL, FORCE, POS_VEL, POS_VEL_FORCE Angstrom, a.u., a.u. 158!> 159!> \param particle_set ... 160!> \param iunit ... 161!> \param output_format ... 162!> \param content ... 163!> \param title ... 164!> \param cell ... 165!> \param array ... 166!> \param unit_conv ... 167!> \param charge_occup ... 168!> \param charge_beta ... 169!> \param charge_extended ... 170!> \date 14.01.2002 171!> \author MK 172!> \version 1.0 173! ************************************************************************************************** 174 SUBROUTINE write_particle_coordinates(particle_set, iunit, output_format, & 175 content, title, cell, array, unit_conv, & 176 charge_occup, charge_beta, & 177 charge_extended) 178 179 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 180 INTEGER :: iunit, output_format 181 CHARACTER(LEN=*) :: content, title 182 TYPE(cell_type), OPTIONAL, POINTER :: cell 183 REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: array 184 REAL(KIND=dp), INTENT(IN), OPTIONAL :: unit_conv 185 LOGICAL, INTENT(IN), OPTIONAL :: charge_occup, charge_beta, & 186 charge_extended 187 188 CHARACTER(len=*), PARAMETER :: routineN = 'write_particle_coordinates', & 189 routineP = moduleN//':'//routineN 190 191 CHARACTER(LEN=120) :: line 192 CHARACTER(LEN=2) :: element_symbol 193 CHARACTER(LEN=4) :: name 194 CHARACTER(LEN=default_string_length) :: atm_name, my_format 195 INTEGER :: handle, iatom, natom 196 LOGICAL :: dummy, my_charge_beta, & 197 my_charge_extended, my_charge_occup 198 REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma, & 199 factor, qeff 200 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: arr 201 REAL(KIND=dp), DIMENSION(3) :: abc, angles, f, r, v 202 REAL(KIND=dp), DIMENSION(3, 3) :: h 203 REAL(KIND=sp), ALLOCATABLE, DIMENSION(:) :: x4, y4, z4 204 TYPE(cell_type), POINTER :: cell_dcd 205 TYPE(fist_potential_type), POINTER :: fist_potential 206 TYPE(shell_kind_type), POINTER :: shell 207 208 CALL timeset(routineN, handle) 209 210 natom = SIZE(particle_set) 211 IF (PRESENT(array)) THEN 212 SELECT CASE (TRIM(content)) 213 CASE ("POS_VEL", "POS_VEL_FORCE") 214 CPABORT("Illegal usage") 215 END SELECT 216 END IF 217 factor = 1.0_dp 218 IF (PRESENT(unit_conv)) THEN 219 factor = unit_conv 220 END IF 221 SELECT CASE (output_format) 222 CASE (dump_xmol) 223 WRITE (iunit, "(I8)") natom 224 WRITE (iunit, "(A)") TRIM(title) 225 DO iatom = 1, natom 226 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 227 element_symbol=element_symbol) 228 IF (LEN_TRIM(element_symbol) == 0) THEN 229 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 230 name=atm_name) 231 dummy = qmmm_ff_precond_only_qm(id1=atm_name) 232 my_format = "(A4," 233 name = TRIM(atm_name) 234 ELSE 235 my_format = "(T2,A2," 236 name = TRIM(element_symbol) 237 END IF 238 SELECT CASE (TRIM(content)) 239 CASE ("POS") 240 IF (PRESENT(array)) THEN 241 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array) 242 ELSE 243 r(:) = particle_set(iatom)%r(:) 244 END IF 245 WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(name), r(1:3)*factor 246 CASE ("VEL") 247 IF (PRESENT(array)) THEN 248 v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array) 249 ELSE 250 v(:) = particle_set(iatom)%v(:) 251 END IF 252 WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(name), v(1:3)*factor 253 CASE ("FORCE") 254 IF (PRESENT(array)) THEN 255 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3) 256 ELSE 257 f(:) = particle_set(iatom)%f(:) 258 END IF 259 WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(name), f(1:3)*factor 260 CASE ("FORCE_MIXING_LABELS") 261 IF (PRESENT(array)) THEN 262 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3) 263 ELSE 264 f(:) = particle_set(iatom)%f(:) 265 END IF 266 WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(name), f(1:3)*factor 267 END SELECT 268 END DO 269 CASE (dump_atomic) 270 DO iatom = 1, natom 271 SELECT CASE (TRIM(content)) 272 CASE ("POS") 273 IF (PRESENT(array)) THEN 274 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array) 275 ELSE 276 r(:) = particle_set(iatom)%r(:) 277 END IF 278 WRITE (iunit, "(3F20.10)") r(1:3)*factor 279 CASE ("VEL") 280 IF (PRESENT(array)) THEN 281 v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array) 282 ELSE 283 v(:) = particle_set(iatom)%v(:) 284 END IF 285 WRITE (iunit, "(3F20.10)") v(1:3)*factor 286 CASE ("FORCE") 287 IF (PRESENT(array)) THEN 288 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3) 289 ELSE 290 f(:) = particle_set(iatom)%f(:) 291 END IF 292 WRITE (iunit, "(3F20.10)") f(1:3)*factor 293 CASE ("FORCE_MIXING_LABELS") 294 IF (PRESENT(array)) THEN 295 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3) 296 ELSE 297 f(:) = particle_set(iatom)%f(:) 298 END IF 299 WRITE (iunit, "(3F20.10)") f(1:3)*factor 300 END SELECT 301 END DO 302 CASE (dump_dcd, dump_dcd_aligned_cell) 303 IF (.NOT. (PRESENT(cell))) THEN 304 CPABORT("Cell is not present! Report this bug!") 305 END IF 306 CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, & 307 abc=abc) 308 IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN 309 ! In the case of a non-orthorhombic cell adopt a common convention 310 ! for the orientation of the cell with respect to the Cartesian axes: 311 ! Cell vector a is aligned with the x axis and the cell vector b lies 312 ! in the xy plane. 313 NULLIFY (cell_dcd) 314 CALL cell_create(cell_dcd) 315 CALL cell_clone(cell, cell_dcd) 316 angles(1) = angle_alpha/degree 317 angles(2) = angle_beta/degree 318 angles(3) = angle_gamma/degree 319 CALL set_cell_param(cell_dcd, abc, angles, & 320 do_init_cell=.TRUE.) 321 h(1:3, 1:3) = MATMUL(cell_dcd%hmat(1:3, 1:3), cell%h_inv(1:3, 1:3)) 322 CALL cell_release(cell_dcd) 323 END IF 324 ALLOCATE (arr(3, natom)) 325 IF (PRESENT(array)) THEN 326 arr(1:3, 1:natom) = RESHAPE(array, (/3, natom/)) 327 ELSE 328 SELECT CASE (TRIM(content)) 329 CASE ("POS") 330 DO iatom = 1, natom 331 arr(1:3, iatom) = particle_set(iatom)%r(1:3) 332 END DO 333 CASE ("VEL") 334 DO iatom = 1, natom 335 arr(1:3, iatom) = particle_set(iatom)%v(1:3) 336 END DO 337 CASE ("FORCE") 338 DO iatom = 1, natom 339 arr(1:3, iatom) = particle_set(iatom)%f(1:3) 340 END DO 341 CASE DEFAULT 342 CPABORT("Illegal DCD dump type") 343 END SELECT 344 END IF 345 ALLOCATE (x4(natom)) 346 ALLOCATE (y4(natom)) 347 ALLOCATE (z4(natom)) 348 IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN 349 x4(1:natom) = REAL(MATMUL(h(1, 1:3), arr(1:3, 1:natom)), KIND=sp) 350 y4(1:natom) = REAL(MATMUL(h(2, 1:3), arr(1:3, 1:natom)), KIND=sp) 351 z4(1:natom) = REAL(MATMUL(h(3, 1:3), arr(1:3, 1:natom)), KIND=sp) 352 ELSE 353 x4(1:natom) = REAL(arr(1, 1:natom), KIND=sp) 354 y4(1:natom) = REAL(arr(2, 1:natom), KIND=sp) 355 z4(1:natom) = REAL(arr(3, 1:natom), KIND=sp) 356 END IF 357 WRITE (iunit) abc(1)*factor, angle_gamma, abc(2)*factor, & 358 angle_beta, angle_alpha, abc(3)*factor 359 WRITE (iunit) x4*REAL(factor, KIND=sp) 360 WRITE (iunit) y4*REAL(factor, KIND=sp) 361 WRITE (iunit) z4*REAL(factor, KIND=sp) 362 ! Release work storage 363 DEALLOCATE (arr) 364 DEALLOCATE (x4) 365 DEALLOCATE (y4) 366 DEALLOCATE (z4) 367 CASE (dump_pdb) 368 my_charge_occup = .FALSE. 369 IF (PRESENT(charge_occup)) my_charge_occup = charge_occup 370 my_charge_beta = .FALSE. 371 IF (PRESENT(charge_beta)) my_charge_beta = charge_beta 372 my_charge_extended = .FALSE. 373 IF (PRESENT(charge_extended)) my_charge_extended = charge_extended 374 IF (LEN_TRIM(title) > 0) THEN 375 WRITE (UNIT=iunit, FMT="(A6,T11,A)") & 376 "REMARK", TRIM(title) 377 END IF 378 CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc) 379 ! COLUMNS DATA TYPE CONTENTS 380 ! -------------------------------------------------- 381 ! 1 - 6 Record name "CRYST1" 382 ! 7 - 15 Real(9.3) a (Angstroms) 383 ! 16 - 24 Real(9.3) b (Angstroms) 384 ! 25 - 33 Real(9.3) c (Angstroms) 385 ! 34 - 40 Real(7.2) alpha (degrees) 386 ! 41 - 47 Real(7.2) beta (degrees) 387 ! 48 - 54 Real(7.2) gamma (degrees) 388 ! 56 - 66 LString Space group 389 ! 67 - 70 Integer Z value 390 WRITE (UNIT=iunit, FMT="(A6,3F9.3,3F7.2)") & 391 "CRYST1", abc(1:3)*factor, angle_alpha, angle_beta, angle_gamma 392 WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM " 393 DO iatom = 1, natom 394 line = "" 395 ! COLUMNS DATA TYPE CONTENTS 396 ! 1 - 6 Record name "ATOM " 397 ! 7 - 11 Integer Atom serial number 398 ! 13 - 16 Atom Atom name 399 ! 17 Character Alternate location indicator 400 ! 18 - 20 Residue name Residue name 401 ! 22 Character Chain identifier 402 ! 23 - 26 Integer Residue sequence number 403 ! 27 AChar Code for insertion of residues 404 ! 31 - 38 Real(8.3) Orthogonal coordinates for X in Angstrom 405 ! 39 - 46 Real(8.3) Orthogonal coordinates for Y in Angstrom 406 ! 47 - 54 Real(8.3) Orthogonal coordinates for Z in Angstrom 407 ! 55 - 60 Real(6.2) Occupancy 408 ! 61 - 66 Real(6.2) Temperature factor (Default = 0.0) 409 ! 73 - 76 LString(4) Segment identifier, left-justified 410 ! 77 - 78 LString(2) Element symbol, right-justified 411 ! 79 - 80 LString(2) Charge on the atom 412 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 413 element_symbol=element_symbol, name=atm_name, & 414 fist_potential=fist_potential, shell=shell) 415 IF (LEN_TRIM(element_symbol) == 0) THEN 416 dummy = qmmm_ff_precond_only_qm(id1=atm_name) 417 END IF 418 name = TRIM(atm_name) 419 IF (ASSOCIATED(fist_potential)) THEN 420 CALL get_potential(potential=fist_potential, qeff=qeff) 421 ELSE 422 qeff = 0.0_dp 423 END IF 424 IF (ASSOCIATED(shell)) CALL get_shell(shell=shell, charge=qeff) 425 WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM " 426 WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(iatom, 100000) 427 WRITE (UNIT=line(13:16), FMT="(A4)") ADJUSTL(name) 428 ! WRITE (UNIT=line(18:20),FMT="(A3)") TRIM(resname) 429 ! WRITE (UNIT=line(23:26),FMT="(I4)") MODULO(idres,10000) 430 SELECT CASE (TRIM(content)) 431 CASE ("POS") 432 IF (PRESENT(array)) THEN 433 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array) 434 ELSE 435 r(:) = particle_set(iatom)%r(:) 436 END IF 437 WRITE (UNIT=line(31:54), FMT="(3F8.3)") r(1:3)*factor 438 CASE DEFAULT 439 CPABORT("PDB dump only for trajectory available") 440 END SELECT 441 IF (my_charge_occup) THEN 442 WRITE (UNIT=line(55:60), FMT="(F6.2)") qeff 443 ELSE 444 WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp 445 END IF 446 IF (my_charge_beta) THEN 447 WRITE (UNIT=line(61:66), FMT="(F6.2)") qeff 448 ELSE 449 WRITE (UNIT=line(61:66), FMT="(F6.2)") 0.0_dp 450 END IF 451 ! WRITE (UNIT=line(73:76),FMT="(A4)") ADJUSTL(TRIM(molname)) 452 WRITE (UNIT=line(77:78), FMT="(A2)") ADJUSTR(TRIM(element_symbol)) 453 IF (my_charge_extended) THEN 454 WRITE (UNIT=line(81:), FMT="(SP,F0.8)") qeff 455 END IF 456 WRITE (UNIT=iunit, FMT="(A)") TRIM(line) 457 END DO 458 WRITE (UNIT=iunit, FMT="(A)") "END" 459 CASE DEFAULT 460 CPABORT("Illegal dump type") 461 END SELECT 462 463 CALL timestop(handle) 464 465 END SUBROUTINE write_particle_coordinates 466 467! ************************************************************************************************** 468!> \brief Write the atomic coordinates to the output unit. 469!> \param particle_set ... 470!> \param subsys_section ... 471!> \param charges ... 472!> \date 05.06.2000 473!> \author MK 474!> \version 1.0 475! ************************************************************************************************** 476 SUBROUTINE write_fist_particle_coordinates(particle_set, subsys_section, & 477 charges) 478 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 479 TYPE(section_vals_type), POINTER :: subsys_section 480 REAL(KIND=dp), DIMENSION(:), POINTER :: charges 481 482 CHARACTER(len=*), PARAMETER :: routineN = 'write_fist_particle_coordinates', & 483 routineP = moduleN//':'//routineN 484 485 CHARACTER(LEN=default_string_length) :: name, unit_str 486 INTEGER :: iatom, ikind, iw, natom 487 REAL(KIND=dp) :: conv, mass, qcore, qeff, qshell 488 TYPE(cp_logger_type), POINTER :: logger 489 TYPE(shell_kind_type), POINTER :: shell_kind 490 491 NULLIFY (logger) 492 NULLIFY (shell_kind) 493 494 logger => cp_get_default_logger() 495 iw = cp_print_key_unit_nr(logger, subsys_section, & 496 "PRINT%ATOMIC_COORDINATES", extension=".coordLog") 497 498 CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str) 499 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 500 CALL uppercase(unit_str) 501 IF (iw > 0) THEN 502 WRITE (UNIT=iw, FMT="(/,/,T2,A)") & 503 "MODULE FIST: ATOMIC COORDINATES IN "//TRIM(unit_str) 504 WRITE (UNIT=iw, & 505 FMT="(/,T3,A,7X,2(A1,11X),A1,8X,A8,5X,A6,/)") & 506 "Atom Kind ATM_TYP", "X", "Y", "Z", " q(eff)", " Mass" 507 natom = SIZE(particle_set) 508 DO iatom = 1, natom 509 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 510 kind_number=ikind, & 511 name=name, & 512 mass=mass, & 513 qeff=qeff, & 514 shell=shell_kind) 515 IF (ASSOCIATED(charges)) qeff = charges(iatom) 516 IF (ASSOCIATED(shell_kind)) THEN 517 CALL get_shell(shell=shell_kind, & 518 charge_core=qcore, & 519 charge_shell=qshell) 520 qeff = qcore + qshell 521 END IF 522 WRITE (UNIT=iw, & 523 FMT="(T2,I5,1X,I4,3X,A4,3X,3F12.6,4X,F6.2,2X,F11.4)") & 524 iatom, ikind, name, & 525 particle_set(iatom)%r(1:3)*conv, qeff, mass/massunit 526 END DO 527 WRITE (iw, '(/)') 528 END IF 529 530 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 531 "PRINT%ATOMIC_COORDINATES") 532 533 END SUBROUTINE write_fist_particle_coordinates 534 535! ************************************************************************************************** 536!> \brief Write the atomic coordinates to the output unit. 537!> \param particle_set ... 538!> \param qs_kind_set ... 539!> \param subsys_section ... 540!> \param label ... 541!> \date 05.06.2000 542!> \author MK 543!> \version 1.0 544! ************************************************************************************************** 545 SUBROUTINE write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label) 546 547 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 548 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set 549 TYPE(section_vals_type), POINTER :: subsys_section 550 CHARACTER(LEN=*), INTENT(IN) :: label 551 552 CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates', & 553 routineP = moduleN//':'//routineN 554 555 CHARACTER(LEN=2) :: element_symbol 556 CHARACTER(LEN=default_string_length) :: unit_str 557 INTEGER :: handle, iatom, ikind, iw, natom, z 558 REAL(KIND=dp) :: conv, mass, zeff 559 TYPE(cp_logger_type), POINTER :: logger 560 561 CALL timeset(routineN, handle) 562 563 NULLIFY (logger) 564 logger => cp_get_default_logger() 565 iw = cp_print_key_unit_nr(logger, subsys_section, & 566 "PRINT%ATOMIC_COORDINATES", extension=".coordLog") 567 568 CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str) 569 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 570 IF (iw > 0) THEN 571 WRITE (UNIT=iw, FMT="(/,/,T2,A)") & 572 "MODULE "//TRIM(label)//": ATOMIC COORDINATES IN "//TRIM(unit_str) 573 WRITE (UNIT=iw, & 574 FMT="(/,T3,A,7X,2(A1,11X),A1,8X,A8,5X,A6,/)") & 575 "Atom Kind Element", "X", "Y", "Z", " Z(eff)", " Mass" 576 577 natom = SIZE(particle_set) 578 DO iatom = 1, natom 579 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, & 580 kind_number=ikind, & 581 element_symbol=element_symbol, & 582 mass=mass, & 583 z=z) 584 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff) 585 WRITE (UNIT=iw, & 586 FMT="(T2,I7,1X,I5,1X,A2,1X,I3,3F12.6,4X,F8.4,2X,F11.4)") & 587 iatom, ikind, element_symbol, z, & 588 particle_set(iatom)%r(1:3)*conv, zeff, mass/massunit 589 END DO 590 WRITE (iw, '(/)') 591 END IF 592 593 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 594 "PRINT%ATOMIC_COORDINATES") 595 596 CALL timestop(handle) 597 598 END SUBROUTINE write_qs_particle_coordinates 599 600! ************************************************************************************************** 601!> \brief Write the matrix of the particle distances to the output unit. 602!> \param particle_set ... 603!> \param cell ... 604!> \param subsys_section ... 605!> \date 06.10.2000 606!> \author Matthias Krack 607!> \version 1.0 608! ************************************************************************************************** 609 SUBROUTINE write_particle_distances(particle_set, cell, subsys_section) 610 611 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 612 TYPE(cell_type), POINTER :: cell 613 TYPE(section_vals_type), POINTER :: subsys_section 614 615 CHARACTER(len=*), PARAMETER :: routineN = 'write_particle_distances', & 616 routineP = moduleN//':'//routineN 617 618 CHARACTER(LEN=default_string_length) :: unit_str 619 INTEGER :: handle, iatom, iw, jatom, natom 620 INTEGER, DIMENSION(3) :: periodic 621 REAL(KIND=dp) :: conv, dab 622 REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: distance_matrix 623 REAL(KIND=dp), DIMENSION(3) :: rab 624 TYPE(cp_logger_type), POINTER :: logger 625 626 CALL timeset(routineN, handle) 627 628 NULLIFY (logger) 629 logger => cp_get_default_logger() 630 iw = cp_print_key_unit_nr(logger, subsys_section, & 631 "PRINT%INTERATOMIC_DISTANCES", extension=".distLog") 632 633 CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%UNIT", c_val=unit_str) 634 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 635 IF (iw > 0) THEN 636 CALL get_cell(cell=cell, periodic=periodic) 637 natom = SIZE(particle_set) 638 ALLOCATE (distance_matrix(natom, natom)) 639 distance_matrix(:, :) = 0.0_dp 640 DO iatom = 1, natom 641 DO jatom = iatom + 1, natom 642 rab(:) = pbc(particle_set(iatom)%r(:), & 643 particle_set(jatom)%r(:), cell) 644 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 645 distance_matrix(iatom, jatom) = dab*conv 646 distance_matrix(jatom, iatom) = distance_matrix(iatom, jatom) 647 END DO 648 END DO 649 650 ! Print the distance matrix 651 WRITE (UNIT=iw, FMT="(/,/,T2,A)") & 652 "INTERATOMIC DISTANCES IN "//TRIM(unit_str) 653 654 CALL write_particle_matrix(distance_matrix, particle_set, iw) 655 END IF 656 657 CALL cp_print_key_finished_output(iw, logger, subsys_section, & 658 "PRINT%INTERATOMIC_DISTANCES") 659 660 CALL timestop(handle) 661 662 END SUBROUTINE write_particle_distances 663 664! ************************************************************************************************** 665!> \brief ... 666!> \param matrix ... 667!> \param particle_set ... 668!> \param iw ... 669!> \param el_per_part ... 670!> \param Ilist ... 671!> \param parts_per_line : number of particle columns to be printed in one line 672! ************************************************************************************************** 673 SUBROUTINE write_particle_matrix(matrix, particle_set, iw, el_per_part, Ilist, parts_per_line) 674 REAL(KIND=dp), DIMENSION(:, :) :: matrix 675 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 676 INTEGER, INTENT(IN) :: iw 677 INTEGER, INTENT(IN), OPTIONAL :: el_per_part 678 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist 679 INTEGER, INTENT(IN), OPTIONAL :: parts_per_line 680 681 CHARACTER(LEN=2) :: element_symbol 682 CHARACTER(LEN=default_string_length) :: fmt_string1, fmt_string2 683 INTEGER :: from, i, iatom, icol, jatom, katom, & 684 my_el_per_part, my_parts_per_line, & 685 natom, to 686 INTEGER, DIMENSION(:), POINTER :: my_list 687 688 my_el_per_part = 1 689 IF (PRESENT(el_per_part)) my_el_per_part = el_per_part 690 my_parts_per_line = 5 691 IF (PRESENT(parts_per_line)) my_parts_per_line = MAX(parts_per_line, 1) 692 WRITE (fmt_string1, FMT='(A,I0,A)') & 693 "(/,T2,11X,", my_parts_per_line, "(4X,I5,4X))" 694 WRITE (fmt_string2, FMT='(A,I0,A)') & 695 "(T2,I5,2X,A2,2X,", my_parts_per_line, "(1X,F12.6))" 696 IF (PRESENT(Ilist)) THEN 697 natom = SIZE(Ilist) 698 ELSE 699 natom = SIZE(particle_set) 700 END IF 701 ALLOCATE (my_list(natom)) 702 IF (PRESENT(Ilist)) THEN 703 my_list = Ilist 704 ELSE 705 DO i = 1, natom 706 my_list(i) = i 707 END DO 708 END IF 709 natom = natom*my_el_per_part 710 DO jatom = 1, natom, my_parts_per_line 711 from = jatom 712 to = MIN(from + my_parts_per_line - 1, natom) 713 WRITE (UNIT=iw, FMT=TRIM(fmt_string1)) (icol, icol=from, to) 714 DO iatom = 1, natom 715 katom = iatom/my_el_per_part 716 IF (MOD(iatom, my_el_per_part) /= 0) katom = katom + 1 717 CALL get_atomic_kind(atomic_kind=particle_set(my_list(katom))%atomic_kind, & 718 element_symbol=element_symbol) 719 WRITE (UNIT=iw, FMT=TRIM(fmt_string2)) & 720 iatom, element_symbol, & 721 (matrix(iatom, icol), icol=from, to) 722 END DO 723 END DO 724 DEALLOCATE (my_list) 725 END SUBROUTINE write_particle_matrix 726 727! ************************************************************************************************** 728!> \brief Write structure data requested by a separate structure data input 729!> section to the output unit. 730!> input_section can be either motion_section or subsys_section. 731!> 732!> \param particle_set ... 733!> \param cell ... 734!> \param input_section ... 735!> \date 11.03.04 736!> \par History 737!> Recovered (23.03.06,MK) 738!> \author MK 739!> \version 1.0 740! ************************************************************************************************** 741 SUBROUTINE write_structure_data(particle_set, cell, input_section) 742 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 743 TYPE(cell_type), POINTER :: cell 744 TYPE(section_vals_type), POINTER :: input_section 745 746 CHARACTER(LEN=*), PARAMETER :: routineN = 'write_structure_data', & 747 routineP = moduleN//':'//routineN 748 749 CHARACTER(LEN=default_string_length) :: string, unit_str 750 INTEGER :: handle, i, i_rep, iw, n, n_rep, n_vals, & 751 natom, new_size, old_size, wrk2(2), & 752 wrk3(3), wrk4(4) 753 INTEGER, ALLOCATABLE, DIMENSION(:) :: work 754 INTEGER, DIMENSION(:), POINTER :: atomic_indices, index_list 755 LOGICAL :: unique 756 REAL(KIND=dp) :: conv, dab 757 REAL(KIND=dp), DIMENSION(3) :: r, rab, rbc, rcd, s 758 TYPE(cp_logger_type), POINTER :: logger 759 TYPE(section_vals_type), POINTER :: section 760 761 CALL timeset(routineN, handle) 762 NULLIFY (atomic_indices) 763 NULLIFY (index_list) 764 NULLIFY (logger) 765 NULLIFY (section) 766 string = "" 767 768 logger => cp_get_default_logger() 769 iw = cp_print_key_unit_nr(logger=logger, & 770 basis_section=input_section, & 771 print_key_path="PRINT%STRUCTURE_DATA", & 772 extension=".coordLog") 773 774 CALL section_vals_val_get(input_section, "PRINT%STRUCTURE_DATA%UNIT", c_val=unit_str) 775 conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str)) 776 CALL uppercase(unit_str) 777 IF (iw > 0) THEN 778 natom = SIZE(particle_set) 779 section => section_vals_get_subs_vals(section_vals=input_section, & 780 subsection_name="PRINT%STRUCTURE_DATA") 781 782 WRITE (UNIT=iw, FMT="(/,T2,A)") "REQUESTED STRUCTURE DATA" 783 ! Print the requested atomic position vectors 784 CALL section_vals_val_get(section_vals=section, & 785 keyword_name="POSITION", & 786 n_rep_val=n_rep) 787 IF (n_rep > 0) THEN 788 WRITE (UNIT=iw, FMT="(/,T3,A,/)") & 789 "Position vectors r(i) of the atoms i in "//TRIM(unit_str) 790 old_size = 0 791 DO i_rep = 1, n_rep 792 CALL section_vals_val_get(section_vals=section, & 793 keyword_name="POSITION", & 794 i_rep_val=i_rep, & 795 i_vals=atomic_indices) 796 n_vals = SIZE(atomic_indices) 797 new_size = old_size + n_vals 798 CALL reallocate(index_list, 1, new_size) 799 index_list(old_size + 1:new_size) = atomic_indices(1:n_vals) 800 old_size = new_size 801 END DO 802 ALLOCATE (work(new_size)) 803 CALL sort(index_list, new_size, work) 804 DEALLOCATE (work) 805 DO i = 1, new_size 806 WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")" 807 IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN 808 WRITE (UNIT=iw, FMT="(T3,A)") & 809 "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored." 810 CYCLE 811 END IF 812 IF (i > 1) THEN 813 ! Skip redundant indices 814 IF (index_list(i) == index_list(i - 1)) CYCLE 815 END IF 816 WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") & 817 "r"//TRIM(string), "=", pbc(particle_set(index_list(i))%r(1:3), cell)*conv 818 END DO 819 DEALLOCATE (index_list) 820 END IF 821 822 ! Print the requested atomic position vectors in scaled coordinates 823 CALL section_vals_val_get(section_vals=section, & 824 keyword_name="POSITION_SCALED", & 825 n_rep_val=n_rep) 826 IF (n_rep > 0) THEN 827 WRITE (UNIT=iw, FMT="(/,T3,A,/)") & 828 "Position vectors s(i) of the atoms i in scaled coordinates" 829 old_size = 0 830 DO i_rep = 1, n_rep 831 CALL section_vals_val_get(section_vals=section, & 832 keyword_name="POSITION_SCALED", & 833 i_rep_val=i_rep, & 834 i_vals=atomic_indices) 835 n_vals = SIZE(atomic_indices) 836 new_size = old_size + n_vals 837 CALL reallocate(index_list, 1, new_size) 838 index_list(old_size + 1:new_size) = atomic_indices(1:n_vals) 839 old_size = new_size 840 END DO 841 ALLOCATE (work(new_size)) 842 CALL sort(index_list, new_size, work) 843 DEALLOCATE (work) 844 DO i = 1, new_size 845 WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")" 846 IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN 847 WRITE (UNIT=iw, FMT="(T3,A)") & 848 "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored." 849 CYCLE 850 END IF 851 IF (i > 1) THEN 852 ! Skip redundant indices 853 IF (index_list(i) == index_list(i - 1)) CYCLE 854 END IF 855 r(1:3) = pbc(particle_set(index_list(i))%r(1:3), cell) 856 CALL real_to_scaled(s, r, cell) 857 WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") & 858 "s"//TRIM(string), "=", s(1:3) 859 END DO 860 DEALLOCATE (index_list) 861 END IF 862 863 ! Print the requested distances 864 CALL section_vals_val_get(section_vals=section, & 865 keyword_name="DISTANCE", & 866 n_rep_val=n) 867 IF (n > 0) THEN 868 WRITE (UNIT=iw, FMT="(/,T3,A,/)") & 869 "Distance vector r(i,j) between the atom i and j in "// & 870 TRIM(unit_str) 871 DO i = 1, n 872 CALL section_vals_val_get(section_vals=section, & 873 keyword_name="DISTANCE", & 874 i_rep_val=i, & 875 i_vals=atomic_indices) 876 string = "" 877 WRITE (UNIT=string, FMT="(A,2(I0,A))") & 878 "(", atomic_indices(1), ",", atomic_indices(2), ")" 879 wrk2 = atomic_indices 880 CALL sort_unique(wrk2, unique) 881 IF (((wrk2(1) >= 1) .AND. (wrk2(SIZE(wrk2)) <= natom)) .AND. unique) THEN 882 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), & 883 particle_set(atomic_indices(2))%r(:), cell) 884 dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)) 885 WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6,3X,A,F13.6)") & 886 "r"//TRIM(string), "=", rab(:)*conv, & 887 "|r| =", dab*conv 888 ELSE 889 WRITE (UNIT=iw, FMT="(T3,A)") & 890 "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored." 891 END IF 892 END DO 893 END IF 894 895 ! Print the requested angles 896 CALL section_vals_val_get(section_vals=section, & 897 keyword_name="ANGLE", & 898 n_rep_val=n) 899 IF (n > 0) THEN 900 WRITE (UNIT=iw, FMT="(/,T3,A,/)") & 901 "Angle a(i,j,k) between the atomic distance vectors r(j,i) and "// & 902 "r(j,k) in DEGREE" 903 DO i = 1, n 904 CALL section_vals_val_get(section_vals=section, & 905 keyword_name="ANGLE", & 906 i_rep_val=i, & 907 i_vals=atomic_indices) 908 string = "" 909 WRITE (UNIT=string, FMT="(A,3(I0,A))") & 910 "(", atomic_indices(1), ",", atomic_indices(2), ",", atomic_indices(3), ")" 911 wrk3 = atomic_indices 912 CALL sort_unique(wrk3, unique) 913 IF (((wrk3(1) >= 1) .AND. (wrk3(SIZE(wrk3)) <= natom)) .AND. unique) THEN 914 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), & 915 particle_set(atomic_indices(2))%r(:), cell) 916 rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), & 917 particle_set(atomic_indices(3))%r(:), cell) 918 WRITE (UNIT=iw, FMT="(T3,A,T26,A,F9.3)") & 919 "a"//TRIM(string), "=", angle(-rab, rbc)*degree 920 ELSE 921 WRITE (UNIT=iw, FMT="(T3,A)") & 922 "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored." 923 END IF 924 END DO 925 END IF 926 927 ! Print the requested dihedral angles 928 CALL section_vals_val_get(section_vals=section, & 929 keyword_name="DIHEDRAL_ANGLE", & 930 n_rep_val=n) 931 IF (n > 0) THEN 932 WRITE (UNIT=iw, FMT="(/,T3,A,/)") & 933 "Dihedral angle d(i,j,k,l) between the planes (i,j,k) and (j,k,l) "// & 934 "in DEGREE" 935 DO i = 1, n 936 CALL section_vals_val_get(section_vals=section, & 937 keyword_name="DIHEDRAL_ANGLE", & 938 i_rep_val=i, & 939 i_vals=atomic_indices) 940 string = "" 941 WRITE (UNIT=string, FMT="(A,4(I0,A))") & 942 "(", atomic_indices(1), ",", atomic_indices(2), ",", & 943 atomic_indices(3), ",", atomic_indices(4), ")" 944 wrk4 = atomic_indices 945 CALL sort_unique(wrk4, unique) 946 IF (((wrk4(1) >= 1) .AND. (wrk4(SIZE(wrk4)) <= natom)) .AND. unique) THEN 947 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), & 948 particle_set(atomic_indices(2))%r(:), cell) 949 rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), & 950 particle_set(atomic_indices(3))%r(:), cell) 951 rcd(:) = pbc(particle_set(atomic_indices(3))%r(:), & 952 particle_set(atomic_indices(4))%r(:), cell) 953 WRITE (UNIT=iw, FMT="(T3,A,T26,A,F9.3)") & 954 "d"//TRIM(string), "=", dihedral_angle(rab, rbc, rcd)*degree 955 ELSE 956 WRITE (UNIT=iw, FMT="(T3,A)") & 957 "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored." 958 END IF 959 END DO 960 END IF 961 END IF 962 CALL cp_print_key_finished_output(iw, logger, input_section, & 963 "PRINT%STRUCTURE_DATA") 964 965 CALL timestop(handle) 966 967 END SUBROUTINE write_structure_data 968 969END MODULE particle_methods 970