(git:24d69ee)
Loading...
Searching...
No Matches
particle_methods.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Define methods related to particle_type
10!> \par History
11!> 10.2014 Move routines out of particle_types.F [Ole Schuett]
12!> \author Ole Schuett
13! **************************************************************************************************
14MODULE particle_methods
15 USE atomic_kind_types, ONLY: get_atomic_kind
16 USE basis_set_types, ONLY: get_gto_basis_set,&
17 gto_basis_set_p_type
18 USE cell_methods, ONLY: cell_create,&
19 set_cell_param
20 USE cell_types, ONLY: cell_clone,&
21 cell_release,&
22 cell_type,&
23 get_cell,&
24 pbc,&
25 real_to_scaled
26 USE cp2k_info, ONLY: compile_revision,&
27 cp2k_version,&
28 r_cwd,&
29 r_host_name,&
30 r_user_name
31 USE cp_log_handling, ONLY: cp_get_default_logger,&
32 cp_logger_get_default_io_unit,&
33 cp_logger_type,&
34 cp_to_string
35 USE cp_output_handling, ONLY: cp_p_file,&
36 cp_print_key_finished_output,&
37 cp_print_key_generate_filename,&
38 cp_print_key_should_output,&
39 cp_print_key_unit_nr
40 USE cp_units, ONLY: cp_unit_from_cp2k
41 USE external_potential_types, ONLY: fist_potential_type,&
42 get_potential
43 USE input_constants, ONLY: dump_atomic,&
44 dump_dcd,&
45 dump_dcd_aligned_cell,&
46 dump_extxyz,&
47 dump_pdb,&
48 dump_xmol
49 USE input_cp2k_subsys, ONLY: create_cell_section
50 USE input_enumeration_types, ONLY: enum_i2c,&
51 enumeration_type
52 USE input_keyword_types, ONLY: keyword_get,&
53 keyword_type
54 USE input_section_types, ONLY: section_get_keyword,&
55 section_release,&
56 section_type,&
57 section_vals_get_subs_vals,&
58 section_vals_type,&
59 section_vals_val_get
60 USE kinds, ONLY: default_path_length,&
61 default_string_length,&
62 dp,&
63 sp
64 USE machine, ONLY: m_timestamp,&
65 timestamp_length
66 USE mathconstants, ONLY: degree
67 USE mathlib, ONLY: angle,&
68 dihedral_angle,&
69 gcd
70 USE memory_utilities, ONLY: reallocate
71 USE particle_types, ONLY: get_particle_pos_or_vel,&
72 particle_type
73 USE periodic_table, ONLY: nelem
74 USE physcon, ONLY: massunit
75 USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
76 USE qs_kind_types, ONLY: get_qs_kind,&
77 qs_kind_type
78 USE shell_potential_types, ONLY: get_shell,&
79 shell_kind_type
80 USE string_utilities, ONLY: uppercase
81 USE util, ONLY: sort,&
82 sort_unique
83#include "./base/base_uses.f90"
84
85 IMPLICIT NONE
86
87 PRIVATE
88
89 ! Public subroutines
90
91 PUBLIC :: write_fist_particle_coordinates, &
92 write_qs_particle_coordinates, &
93 write_particle_distances, &
94 write_particle_coordinates, &
95 write_structure_data, &
96 get_particle_set, &
97 write_particle_matrix, &
98 write_final_structure
99
100 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'particle_methods'
101
102CONTAINS
103
104! **************************************************************************************************
105!> \brief Get the components of a particle set.
106!> \param particle_set ...
107!> \param qs_kind_set ...
108!> \param first_sgf ...
109!> \param last_sgf ...
110!> \param nsgf ...
111!> \param nmao ...
112!> \param basis ...
113!> \param ncgf ...
114!> \date 14.01.2002
115!> \par History
116!> - particle type cleaned (13.10.2003,MK)
117!> - refactoring and add basis set option (17.08.2010,jhu)
118!> \author MK
119!> \version 1.0
120! **************************************************************************************************
121 SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, &
122 nmao, basis, ncgf)
123
124 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
125 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
126 INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: first_sgf, last_sgf, nsgf, nmao
127 TYPE(gto_basis_set_p_type), DIMENSION(:), OPTIONAL :: basis
128 INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: ncgf
129
130 INTEGER :: ikind, iparticle, isgf, nparticle, ns
131
132 cpassert(ASSOCIATED(particle_set))
133
134 nparticle = SIZE(particle_set)
135 IF (PRESENT(first_sgf)) THEN
136 cpassert(SIZE(first_sgf) >= nparticle)
137 END IF
138 IF (PRESENT(last_sgf)) THEN
139 cpassert(SIZE(last_sgf) >= nparticle)
140 END IF
141 IF (PRESENT(nsgf)) THEN
142 cpassert(SIZE(nsgf) >= nparticle)
143 END IF
144 IF (PRESENT(nmao)) THEN
145 cpassert(SIZE(nmao) >= nparticle)
146 END IF
147 IF (PRESENT(ncgf)) THEN
148 cpassert(SIZE(ncgf) >= nparticle)
149 END IF
150
151 IF (PRESENT(first_sgf) .OR. PRESENT(last_sgf) .OR. PRESENT(nsgf)) THEN
152 isgf = 0
153 DO iparticle = 1, nparticle
154 CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
155 IF (PRESENT(basis)) THEN
156 IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
157 CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
158 ELSE
159 ns = 0
160 END IF
161 ELSE
162 CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns)
163 END IF
164 IF (PRESENT(nsgf)) nsgf(iparticle) = ns
165 IF (PRESENT(first_sgf)) first_sgf(iparticle) = isgf + 1
166 isgf = isgf + ns
167 IF (PRESENT(last_sgf)) last_sgf(iparticle) = isgf
168 END DO
169 END IF
170
171 IF (PRESENT(ncgf)) THEN
172 DO iparticle = 1, nparticle
173 CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
174 IF (PRESENT(basis)) THEN
175 IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
176 CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, ncgf=ns)
177 ELSE
178 ns = 0
179 END IF
180 ELSE
181 CALL get_qs_kind(qs_kind_set(ikind), ncgf=ns)
182 END IF
183 ncgf(iparticle) = ns
184 END DO
185 END IF
186
187 IF (PRESENT(first_sgf)) THEN
188 IF (SIZE(first_sgf) > nparticle) first_sgf(nparticle + 1) = isgf + 1
189 END IF
190
191 IF (PRESENT(nmao)) THEN
192 DO iparticle = 1, nparticle
193 CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
194 CALL get_qs_kind(qs_kind_set(ikind), mao=ns)
195 nmao(iparticle) = ns
196 END DO
197 END IF
198
199 END SUBROUTINE get_particle_set
200
201! **************************************************************************************************
202!> \brief Should be able to write a few formats e.g. xmol, and some binary
203!> format (dcd) some format can be used for x,v,f
204!>
205!> FORMAT CONTENT UNITS x, v, f
206!> XMOL POS, VEL, FORCE, POS_VEL, POS_VEL_FORCE Angstrom, a.u., a.u.
207!>
208!> \param particle_set ...
209!> \param iunit ...
210!> \param output_format ...
211!> \param content ...
212!> \param title ...
213!> \param cell ...
214!> \param array ...
215!> \param unit_conv ...
216!> \param charge_occup ...
217!> \param charge_beta ...
218!> \param charge_extended ...
219!> \param print_kind ...
220!> \date 14.01.2002
221!> \author MK
222!> \version 1.0
223! **************************************************************************************************
224 SUBROUTINE write_particle_coordinates(particle_set, iunit, output_format, &
225 content, title, cell, array, unit_conv, &
226 charge_occup, charge_beta, &
227 charge_extended, print_kind)
228
229 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
230 INTEGER :: iunit, output_format
231 CHARACTER(LEN=*) :: content, title
232 TYPE(cell_type), OPTIONAL, POINTER :: cell
233 REAL(kind=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: array
234 REAL(kind=dp), INTENT(IN), OPTIONAL :: unit_conv
235 LOGICAL, INTENT(IN), OPTIONAL :: charge_occup, charge_beta, &
236 charge_extended, print_kind
237
238 CHARACTER(len=*), PARAMETER :: routinen = 'write_particle_coordinates'
239
240 CHARACTER(LEN=120) :: line
241 CHARACTER(LEN=2) :: element_symbol
242 CHARACTER(LEN=4) :: name
243 CHARACTER(LEN=default_string_length) :: atm_name, my_format
244 INTEGER :: handle, iatom, natom
245 LOGICAL :: dummy, my_charge_beta, &
246 my_charge_extended, my_charge_occup, &
247 my_print_kind
248 REAL(kind=dp) :: angle_alpha, angle_beta, angle_gamma, &
249 factor, qeff
250 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: arr
251 REAL(kind=dp), DIMENSION(3) :: abc, angles, f, r, v
252 REAL(kind=dp), DIMENSION(3, 3) :: h
253 REAL(kind=sp), ALLOCATABLE, DIMENSION(:) :: x4, y4, z4
254 TYPE(cell_type), POINTER :: cell_dcd
255 TYPE(fist_potential_type), POINTER :: fist_potential
256 TYPE(shell_kind_type), POINTER :: shell
257
258 CALL timeset(routinen, handle)
259
260 natom = SIZE(particle_set)
261 IF (PRESENT(array)) THEN
262 SELECT CASE (trim(content))
263 CASE ("POS_VEL", "POS_VEL_FORCE")
264 cpabort("Illegal usage")
265 END SELECT
266 END IF
267 factor = 1.0_dp
268 IF (PRESENT(unit_conv)) THEN
269 factor = unit_conv
270 END IF
271 SELECT CASE (output_format)
272 CASE (dump_xmol, dump_extxyz)
273 my_print_kind = .false.
274 IF (PRESENT(print_kind)) my_print_kind = print_kind
275 WRITE (iunit, "(I8)") natom
276 WRITE (iunit, "(A)") trim(title)
277 DO iatom = 1, natom
278 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
279 element_symbol=element_symbol)
280 IF (len_trim(element_symbol) == 0 .OR. my_print_kind) THEN
281 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
282 name=atm_name)
283 dummy = qmmm_ff_precond_only_qm(id1=atm_name)
284 my_format = "(A,"
285 atm_name = trim(atm_name)
286 ELSE
287 my_format = "(T2,A2,"
288 atm_name = trim(element_symbol)
289 END IF
290 SELECT CASE (trim(content))
291 CASE ("POS")
292 IF (PRESENT(array)) THEN
293 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
294 ELSE
295 r(:) = particle_set(iatom)%r(:)
296 END IF
297 WRITE (iunit, trim(my_format)//"1X,3F20.10)") trim(atm_name), r(1:3)*factor
298 CASE ("VEL")
299 IF (PRESENT(array)) THEN
300 v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
301 ELSE
302 v(:) = particle_set(iatom)%v(:)
303 END IF
304 WRITE (iunit, trim(my_format)//"1X,3F20.10)") trim(atm_name), v(1:3)*factor
305 CASE ("FORCE")
306 IF (PRESENT(array)) THEN
307 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
308 ELSE
309 f(:) = particle_set(iatom)%f(:)
310 END IF
311 WRITE (iunit, trim(my_format)//"1X,3F20.10)") trim(atm_name), f(1:3)*factor
312 CASE ("FORCE_MIXING_LABELS")
313 IF (PRESENT(array)) THEN
314 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
315 ELSE
316 f(:) = particle_set(iatom)%f(:)
317 END IF
318 WRITE (iunit, trim(my_format)//"1X,3F20.10)") trim(atm_name), f(1:3)*factor
319 END SELECT
320 END DO
321 CASE (dump_atomic)
322 DO iatom = 1, natom
323 SELECT CASE (trim(content))
324 CASE ("POS")
325 IF (PRESENT(array)) THEN
326 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
327 ELSE
328 r(:) = particle_set(iatom)%r(:)
329 END IF
330 WRITE (iunit, "(3F20.10)") r(1:3)*factor
331 CASE ("VEL")
332 IF (PRESENT(array)) THEN
333 v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
334 ELSE
335 v(:) = particle_set(iatom)%v(:)
336 END IF
337 WRITE (iunit, "(3F20.10)") v(1:3)*factor
338 CASE ("FORCE")
339 IF (PRESENT(array)) THEN
340 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
341 ELSE
342 f(:) = particle_set(iatom)%f(:)
343 END IF
344 WRITE (iunit, "(3F20.10)") f(1:3)*factor
345 CASE ("FORCE_MIXING_LABELS")
346 IF (PRESENT(array)) THEN
347 f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
348 ELSE
349 f(:) = particle_set(iatom)%f(:)
350 END IF
351 WRITE (iunit, "(3F20.10)") f(1:3)*factor
352 END SELECT
353 END DO
354 CASE (dump_dcd, dump_dcd_aligned_cell)
355 IF (.NOT. (PRESENT(cell))) THEN
356 cpabort("Cell is not present! Report this bug!")
357 END IF
358 CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, &
359 abc=abc)
360 IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN
361 ! In the case of a non-orthorhombic cell adopt a common convention
362 ! for the orientation of the cell with respect to the Cartesian axes:
363 ! Cell vector a is aligned with the x axis and the cell vector b lies
364 ! in the xy plane.
365 NULLIFY (cell_dcd)
366 CALL cell_create(cell_dcd)
367 CALL cell_clone(cell, cell_dcd, tag="CELL_DCD")
368 angles(1) = angle_alpha/degree
369 angles(2) = angle_beta/degree
370 angles(3) = angle_gamma/degree
371 CALL set_cell_param(cell_dcd, abc, angles, &
372 do_init_cell=.true.)
373 h(1:3, 1:3) = matmul(cell_dcd%hmat(1:3, 1:3), cell%h_inv(1:3, 1:3))
374 CALL cell_release(cell_dcd)
375 END IF
376 ALLOCATE (arr(3, natom))
377 IF (PRESENT(array)) THEN
378 arr(1:3, 1:natom) = reshape(array, [3, natom])
379 ELSE
380 SELECT CASE (trim(content))
381 CASE ("POS")
382 DO iatom = 1, natom
383 arr(1:3, iatom) = particle_set(iatom)%r(1:3)
384 END DO
385 CASE ("VEL")
386 DO iatom = 1, natom
387 arr(1:3, iatom) = particle_set(iatom)%v(1:3)
388 END DO
389 CASE ("FORCE")
390 DO iatom = 1, natom
391 arr(1:3, iatom) = particle_set(iatom)%f(1:3)
392 END DO
393 CASE DEFAULT
394 cpabort("Illegal DCD dump type")
395 END SELECT
396 END IF
397 ALLOCATE (x4(natom))
398 ALLOCATE (y4(natom))
399 ALLOCATE (z4(natom))
400 IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN
401 x4(1:natom) = real(matmul(h(1, 1:3), arr(1:3, 1:natom)), kind=sp)
402 y4(1:natom) = real(matmul(h(2, 1:3), arr(1:3, 1:natom)), kind=sp)
403 z4(1:natom) = real(matmul(h(3, 1:3), arr(1:3, 1:natom)), kind=sp)
404 ELSE
405 x4(1:natom) = real(arr(1, 1:natom), kind=sp)
406 y4(1:natom) = real(arr(2, 1:natom), kind=sp)
407 z4(1:natom) = real(arr(3, 1:natom), kind=sp)
408 END IF
409 WRITE (iunit) abc(1)*factor, angle_gamma, abc(2)*factor, &
410 angle_beta, angle_alpha, abc(3)*factor
411 WRITE (iunit) x4*real(factor, kind=sp)
412 WRITE (iunit) y4*real(factor, kind=sp)
413 WRITE (iunit) z4*real(factor, kind=sp)
414 ! Release work storage
415 DEALLOCATE (arr)
416 DEALLOCATE (x4)
417 DEALLOCATE (y4)
418 DEALLOCATE (z4)
419 CASE (dump_pdb)
420 my_charge_occup = .false.
421 IF (PRESENT(charge_occup)) my_charge_occup = charge_occup
422 my_charge_beta = .false.
423 IF (PRESENT(charge_beta)) my_charge_beta = charge_beta
424 my_charge_extended = .false.
425 IF (PRESENT(charge_extended)) my_charge_extended = charge_extended
426 IF (len_trim(title) > 0) THEN
427 WRITE (unit=iunit, fmt="(A6,T11,A)") &
428 "REMARK", trim(title)
429 END IF
430 CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc)
431 ! COLUMNS DATA TYPE CONTENTS
432 ! --------------------------------------------------
433 ! 1 - 6 Record name "CRYST1"
434 ! 7 - 15 Real(9.3) a (Angstroms)
435 ! 16 - 24 Real(9.3) b (Angstroms)
436 ! 25 - 33 Real(9.3) c (Angstroms)
437 ! 34 - 40 Real(7.2) alpha (degrees)
438 ! 41 - 47 Real(7.2) beta (degrees)
439 ! 48 - 54 Real(7.2) gamma (degrees)
440 ! 56 - 66 LString Space group
441 ! 67 - 70 Integer Z value
442 WRITE (unit=iunit, fmt="(A6,3F9.3,3F7.2)") &
443 "CRYST1", abc(1:3)*factor, angle_alpha, angle_beta, angle_gamma
444 WRITE (unit=line(1:6), fmt="(A6)") "ATOM "
445 DO iatom = 1, natom
446 line = ""
447 ! COLUMNS DATA TYPE CONTENTS
448 ! 1 - 6 Record name "ATOM "
449 ! 7 - 11 Integer Atom serial number
450 ! 13 - 16 Atom Atom name
451 ! 17 Character Alternate location indicator
452 ! 18 - 20 Residue name Residue name
453 ! 22 Character Chain identifier
454 ! 23 - 26 Integer Residue sequence number
455 ! 27 AChar Code for insertion of residues
456 ! 31 - 38 Real(8.3) Orthogonal coordinates for X in Angstrom
457 ! 39 - 46 Real(8.3) Orthogonal coordinates for Y in Angstrom
458 ! 47 - 54 Real(8.3) Orthogonal coordinates for Z in Angstrom
459 ! 55 - 60 Real(6.2) Occupancy
460 ! 61 - 66 Real(6.2) Temperature factor (Default = 0.0)
461 ! 73 - 76 LString(4) Segment identifier, left-justified
462 ! 77 - 78 LString(2) Element symbol, right-justified
463 ! 79 - 80 LString(2) Charge on the atom
464 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
465 element_symbol=element_symbol, name=atm_name, &
466 fist_potential=fist_potential, shell=shell)
467 IF (len_trim(element_symbol) == 0) THEN
468 dummy = qmmm_ff_precond_only_qm(id1=atm_name)
469 END IF
470 name = trim(atm_name)
471 IF (ASSOCIATED(fist_potential)) THEN
472 CALL get_potential(potential=fist_potential, qeff=qeff)
473 ELSE
474 qeff = 0.0_dp
475 END IF
476 IF (ASSOCIATED(shell)) CALL get_shell(shell=shell, charge=qeff)
477 WRITE (unit=line(1:6), fmt="(A6)") "ATOM "
478 WRITE (unit=line(7:11), fmt="(I5)") modulo(iatom, 100000)
479 WRITE (unit=line(13:16), fmt="(A4)") adjustl(name)
480 ! WRITE (UNIT=line(18:20),FMT="(A3)") TRIM(resname)
481 ! WRITE (UNIT=line(23:26),FMT="(I4)") MODULO(idres,10000)
482 SELECT CASE (trim(content))
483 CASE ("POS")
484 IF (PRESENT(array)) THEN
485 r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
486 ELSE
487 r(:) = particle_set(iatom)%r(:)
488 END IF
489 WRITE (unit=line(31:54), fmt="(3F8.3)") r(1:3)*factor
490 CASE DEFAULT
491 cpabort("PDB dump only for trajectory available")
492 END SELECT
493 IF (my_charge_occup) THEN
494 WRITE (unit=line(55:60), fmt="(F6.2)") qeff
495 ELSE
496 WRITE (unit=line(55:60), fmt="(F6.2)") 0.0_dp
497 END IF
498 IF (my_charge_beta) THEN
499 WRITE (unit=line(61:66), fmt="(F6.2)") qeff
500 ELSE
501 WRITE (unit=line(61:66), fmt="(F6.2)") 0.0_dp
502 END IF
503 ! WRITE (UNIT=line(73:76),FMT="(A4)") ADJUSTL(TRIM(molname))
504 WRITE (unit=line(77:78), fmt="(A2)") adjustr(trim(element_symbol))
505 IF (my_charge_extended) THEN
506 WRITE (unit=line(81:), fmt="(SP,F0.8)") qeff
507 END IF
508 WRITE (unit=iunit, fmt="(A)") trim(line)
509 END DO
510 WRITE (unit=iunit, fmt="(A)") "END"
511 CASE DEFAULT
512 cpabort("Illegal dump type")
513 END SELECT
514
515 CALL timestop(handle)
516
517 END SUBROUTINE write_particle_coordinates
518
519! **************************************************************************************************
520!> \brief Write the atomic coordinates to the output unit.
521!> \param particle_set ...
522!> \param subsys_section ...
523!> \param charges ...
524!> \date 05.06.2000
525!> \author MK
526!> \version 1.0
527! **************************************************************************************************
528 SUBROUTINE write_fist_particle_coordinates(particle_set, subsys_section, charges)
529 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
530 TYPE(section_vals_type), POINTER :: subsys_section
531 REAL(kind=dp), DIMENSION(:), OPTIONAL :: charges
532
533 CHARACTER(LEN=default_string_length) :: name, unit_str
534 INTEGER :: iatom, ikind, iw, natom
535 REAL(kind=dp) :: conv, mass, qcore, qeff, qshell
536 TYPE(cp_logger_type), POINTER :: logger
537 TYPE(shell_kind_type), POINTER :: shell_kind
538
539 NULLIFY (logger)
540 NULLIFY (shell_kind)
541
542 logger => cp_get_default_logger()
543 iw = cp_print_key_unit_nr(logger, subsys_section, &
544 "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
545
546 CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
547 conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
548 CALL uppercase(unit_str)
549 IF (iw > 0) THEN
550 WRITE (unit=iw, fmt="(/,/,T2,A)") &
551 "MODULE FIST: ATOMIC COORDINATES IN "//trim(unit_str)
552 WRITE (unit=iw, fmt="(/,T4,A,T30,A,T44,A,T58,A,T66,A,T77,A)") &
553 "Atom Kind Name", "X", "Y", "Z", "q(eff)", "Mass"
554 natom = SIZE(particle_set)
555 DO iatom = 1, natom
556 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
557 kind_number=ikind, &
558 name=name, &
559 mass=mass, &
560 qeff=qeff, &
561 shell=shell_kind)
562 IF (PRESENT(charges)) qeff = charges(iatom)
563 IF (ASSOCIATED(shell_kind)) THEN
564 CALL get_shell(shell=shell_kind, &
565 charge_core=qcore, &
566 charge_shell=qshell)
567 qeff = qcore + qshell
568 END IF
569 WRITE (unit=iw, fmt="(T2,I6,1X,I4,1X,A7,3(1X,F13.6),2(1X,F8.4))") &
570 iatom, ikind, name, particle_set(iatom)%r(1:3)*conv, qeff, mass/massunit
571 END DO
572 WRITE (iw, "(A)") ""
573 END IF
574
575 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
576 "PRINT%ATOMIC_COORDINATES")
577
578 END SUBROUTINE write_fist_particle_coordinates
579
580! **************************************************************************************************
581!> \brief Write the atomic coordinates to the output unit.
582!> \param particle_set ...
583!> \param qs_kind_set ...
584!> \param subsys_section ...
585!> \param label ...
586!> \date 05.06.2000
587!> \author MK
588!> \version 1.0
589! **************************************************************************************************
590 SUBROUTINE write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
591
592 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
593 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
594 TYPE(section_vals_type), POINTER :: subsys_section
595 CHARACTER(LEN=*), INTENT(IN) :: label
596
597 CHARACTER(len=*), PARAMETER :: routinen = 'write_qs_particle_coordinates'
598
599 CHARACTER(LEN=2) :: element_symbol
600 CHARACTER(LEN=default_string_length) :: unit_str
601 INTEGER :: handle, iatom, ikind, iw, natom, z
602 REAL(kind=dp) :: conv, mass, zeff
603 TYPE(cp_logger_type), POINTER :: logger
604
605 CALL timeset(routinen, handle)
606
607 NULLIFY (logger)
608 logger => cp_get_default_logger()
609 iw = cp_print_key_unit_nr(logger, subsys_section, &
610 "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
611
612 CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
613 conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
614 CALL uppercase(unit_str)
615 IF (iw > 0) THEN
616 WRITE (unit=iw, fmt="(/,/,T2,A)") &
617 "MODULE "//trim(label)//": ATOMIC COORDINATES IN "//trim(unit_str)
618 WRITE (unit=iw, fmt="(/,T4,A,T30,A,T44,A,T58,A,T66,A,T77,A)") &
619 "Atom Kind Element", "X", "Y", "Z", "Z(eff)", "Mass"
620 natom = SIZE(particle_set)
621 DO iatom = 1, natom
622 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
623 kind_number=ikind, &
624 element_symbol=element_symbol, &
625 mass=mass, &
626 z=z)
627 CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
628 WRITE (unit=iw, fmt="(T2,I6,1X,I4,1X,A2,1X,I4,3(1X,F13.6),2(1X,F8.4))") &
629 iatom, ikind, element_symbol, z, particle_set(iatom)%r(1:3)*conv, zeff, mass/massunit
630 END DO
631 WRITE (iw, "(A)") ""
632 END IF
633
634 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
635 "PRINT%ATOMIC_COORDINATES")
636
637 CALL timestop(handle)
638
639 END SUBROUTINE write_qs_particle_coordinates
640
641! **************************************************************************************************
642!> \brief Write the matrix of the particle distances to the output unit.
643!> \param particle_set ...
644!> \param cell ...
645!> \param subsys_section ...
646!> \date 06.10.2000
647!> \author Matthias Krack
648!> \version 1.0
649! **************************************************************************************************
650 SUBROUTINE write_particle_distances(particle_set, cell, subsys_section)
651
652 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
653 TYPE(cell_type), POINTER :: cell
654 TYPE(section_vals_type), POINTER :: subsys_section
655
656 CHARACTER(len=*), PARAMETER :: routinen = 'write_particle_distances'
657
658 CHARACTER(LEN=default_string_length) :: unit_str
659 INTEGER :: handle, iatom, iw, jatom, natom
660 INTEGER, DIMENSION(3) :: periodic
661 LOGICAL :: explicit
662 REAL(kind=dp) :: conv, dab, dab_abort, dab_min, dab_warn
663 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: distance_matrix
664 REAL(kind=dp), DIMENSION(3) :: rab
665 TYPE(cp_logger_type), POINTER :: logger
666
667 CALL timeset(routinen, handle)
668
669 cpassert(ASSOCIATED(particle_set))
670 cpassert(ASSOCIATED(cell))
671 cpassert(ASSOCIATED(subsys_section))
672
673 NULLIFY (logger)
674 logger => cp_get_default_logger()
675 iw = cp_print_key_unit_nr(logger, subsys_section, &
676 "PRINT%INTERATOMIC_DISTANCES", extension=".distLog")
677
678 CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%UNIT", c_val=unit_str)
679 conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
680 CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%CHECK_INTERATOMIC_DISTANCES", &
681 r_val=dab_min, explicit=explicit)
682
683 dab_abort = 0.0_dp
684 dab_warn = 0.0_dp
685 natom = SIZE(particle_set)
686
687 ! Compute interatomic distances only if their printout or check is explicitly requested
688 ! Disable the default check for systems with more than 3000 atoms
689 IF (explicit .OR. (iw > 0) .OR. (natom <= 2000)) THEN
690 IF (dab_min > 0.0_dp) THEN
691 dab_warn = dab_min*conv
692 ELSE IF (dab_min < 0.0_dp) THEN
693 dab_abort = abs(dab_min)*conv
694 END IF
695 END IF
696
697 IF ((iw > 0) .OR. (dab_abort > 0.0_dp) .OR. (dab_warn > 0.0_dp)) THEN
698 CALL get_cell(cell=cell, periodic=periodic)
699 IF (iw > 0) THEN
700 ALLOCATE (distance_matrix(natom, natom))
701 distance_matrix(:, :) = 0.0_dp
702 END IF
703 DO iatom = 1, natom
704 DO jatom = iatom + 1, natom
705 rab(:) = pbc(particle_set(iatom)%r(:), &
706 particle_set(jatom)%r(:), cell)
707 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))*conv
708 IF (dab_abort > 0.0_dp) THEN
709 ! Stop the run for interatomic distances smaller than the requested threshold
710 IF (dab < dab_abort) THEN
711 CALL cp_abort(__location__, "The distance between the atoms "// &
712 trim(adjustl(cp_to_string(iatom, fmt="(I8)")))//" and "// &
713 trim(adjustl(cp_to_string(jatom, fmt="(I8)")))//" is only "// &
714 trim(adjustl(cp_to_string(dab, fmt="(F6.3)")))//" "// &
715 trim(adjustl(unit_str))//" and thus smaller than the requested threshold of "// &
716 trim(adjustl(cp_to_string(dab_abort, fmt="(F6.3)")))//" "// &
717 trim(adjustl(unit_str)))
718 END IF
719 END IF
720 IF (dab < dab_warn) THEN
721 ! Print warning for interatomic distances smaller than the requested threshold
722 CALL cp_warn(__location__, "The distance between the atoms "// &
723 trim(adjustl(cp_to_string(iatom, fmt="(I8)")))//" and "// &
724 trim(adjustl(cp_to_string(jatom, fmt="(I8)")))//" is only "// &
725 trim(adjustl(cp_to_string(dab, fmt="(F6.3)")))//" "// &
726 trim(adjustl(unit_str))//" and thus smaller than the threshold of "// &
727 trim(adjustl(cp_to_string(dab_warn, fmt="(F6.3)")))//" "// &
728 trim(adjustl(unit_str)))
729 END IF
730 IF (iw > 0) THEN
731 distance_matrix(iatom, jatom) = dab
732 distance_matrix(jatom, iatom) = distance_matrix(iatom, jatom)
733 END IF
734 END DO
735 END DO
736 IF (iw > 0) THEN
737 ! Print the distance matrix
738 WRITE (unit=iw, fmt="(/,/,T2,A)") &
739 "INTERATOMIC DISTANCES IN "//trim(unit_str)
740 CALL write_particle_matrix(distance_matrix, particle_set, iw)
741 IF (ALLOCATED(distance_matrix)) DEALLOCATE (distance_matrix)
742 END IF
743 CALL cp_print_key_finished_output(iw, logger, subsys_section, &
744 "PRINT%INTERATOMIC_DISTANCES")
745 END IF
746
747 CALL timestop(handle)
748
749 END SUBROUTINE write_particle_distances
750
751! **************************************************************************************************
752!> \brief ...
753!> \param matrix ...
754!> \param particle_set ...
755!> \param iw ...
756!> \param el_per_part ...
757!> \param Ilist ...
758!> \param parts_per_line : number of particle columns to be printed in one line
759! **************************************************************************************************
760 SUBROUTINE write_particle_matrix(matrix, particle_set, iw, el_per_part, Ilist, parts_per_line)
761 REAL(kind=dp), DIMENSION(:, :) :: matrix
762 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
763 INTEGER, INTENT(IN) :: iw
764 INTEGER, INTENT(IN), OPTIONAL :: el_per_part
765 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: ilist
766 INTEGER, INTENT(IN), OPTIONAL :: parts_per_line
767
768 CHARACTER(LEN=2) :: element_symbol
769 CHARACTER(LEN=default_string_length) :: fmt_string1, fmt_string2
770 INTEGER :: from, i, iatom, icol, jatom, katom, &
771 my_el_per_part, my_parts_per_line, &
772 natom, to
773 INTEGER, DIMENSION(:), POINTER :: my_list
774
775 my_el_per_part = 1
776 IF (PRESENT(el_per_part)) my_el_per_part = el_per_part
777 my_parts_per_line = 5
778 IF (PRESENT(parts_per_line)) my_parts_per_line = max(parts_per_line, 1)
779 WRITE (fmt_string1, fmt='(A,I0,A)') &
780 "(/,T2,9X,", my_parts_per_line, "(4X,I6,4X))"
781 WRITE (fmt_string2, fmt='(A,I0,A)') &
782 "(T2,I5,1X,A2,1X,", my_parts_per_line, "(1X,ES13.6E2))"
783 IF (PRESENT(ilist)) THEN
784 natom = SIZE(ilist)
785 ELSE
786 natom = SIZE(particle_set)
787 END IF
788 ALLOCATE (my_list(natom))
789 IF (PRESENT(ilist)) THEN
790 my_list = ilist
791 ELSE
792 DO i = 1, natom
793 my_list(i) = i
794 END DO
795 END IF
796 natom = natom*my_el_per_part
797 DO jatom = 1, natom, my_parts_per_line
798 from = jatom
799 to = min(from + my_parts_per_line - 1, natom)
800 WRITE (unit=iw, fmt=trim(fmt_string1)) (icol, icol=from, to)
801 DO iatom = 1, natom
802 katom = iatom/my_el_per_part
803 IF (mod(iatom, my_el_per_part) /= 0) katom = katom + 1
804 CALL get_atomic_kind(atomic_kind=particle_set(my_list(katom))%atomic_kind, &
805 element_symbol=element_symbol)
806 WRITE (unit=iw, fmt=trim(fmt_string2)) &
807 iatom, element_symbol, &
808 (matrix(iatom, icol), icol=from, to)
809 END DO
810 END DO
811
812 DEALLOCATE (my_list)
813
814 END SUBROUTINE write_particle_matrix
815
816! **************************************************************************************************
817!> \brief Write structure data requested by a separate structure data input
818!> section to the output unit.
819!> input_section can be either motion_section or subsys_section.
820!>
821!> \param particle_set ...
822!> \param cell ...
823!> \param input_section ...
824!> \date 11.03.04
825!> \par History
826!> Recovered (23.03.06,MK)
827!> \author MK
828!> \version 1.0
829! **************************************************************************************************
830 SUBROUTINE write_structure_data(particle_set, cell, input_section)
831 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
832 TYPE(cell_type), POINTER :: cell
833 TYPE(section_vals_type), POINTER :: input_section
834
835 CHARACTER(LEN=*), PARAMETER :: routinen = 'write_structure_data'
836
837 CHARACTER(LEN=default_string_length) :: string, unit_str
838 INTEGER :: handle, i, i_rep, iw, n, n_rep, n_vals, &
839 natom, new_size, old_size, wrk2(2), &
840 wrk3(3), wrk4(4)
841 INTEGER, ALLOCATABLE, DIMENSION(:) :: work
842 INTEGER, DIMENSION(:), POINTER :: atomic_indices, index_list
843 LOGICAL :: unique
844 REAL(kind=dp) :: conv, dab
845 REAL(kind=dp), DIMENSION(3) :: r, rab, rbc, rcd, s
846 TYPE(cp_logger_type), POINTER :: logger
847 TYPE(section_vals_type), POINTER :: section
848
849 CALL timeset(routinen, handle)
850 NULLIFY (atomic_indices)
851 NULLIFY (index_list)
852 NULLIFY (logger)
853 NULLIFY (section)
854 string = ""
855
856 logger => cp_get_default_logger()
857 iw = cp_print_key_unit_nr(logger=logger, &
858 basis_section=input_section, &
859 print_key_path="PRINT%STRUCTURE_DATA", &
860 extension=".coordLog")
861
862 CALL section_vals_val_get(input_section, "PRINT%STRUCTURE_DATA%UNIT", c_val=unit_str)
863 conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
864 CALL uppercase(unit_str)
865 IF (iw > 0) THEN
866 natom = SIZE(particle_set)
867 section => section_vals_get_subs_vals(section_vals=input_section, &
868 subsection_name="PRINT%STRUCTURE_DATA")
869
870 WRITE (unit=iw, fmt="(/,T2,A)") "REQUESTED STRUCTURE DATA"
871 ! Print the requested atomic position vectors
872 CALL section_vals_val_get(section_vals=section, &
873 keyword_name="POSITION", &
874 n_rep_val=n_rep)
875 IF (n_rep > 0) THEN
876 WRITE (unit=iw, fmt="(/,T3,A,/)") &
877 "Position vectors r(i) of the atoms i in "//trim(unit_str)
878 old_size = 0
879 DO i_rep = 1, n_rep
880 CALL section_vals_val_get(section_vals=section, &
881 keyword_name="POSITION", &
882 i_rep_val=i_rep, &
883 i_vals=atomic_indices)
884 n_vals = SIZE(atomic_indices)
885 new_size = old_size + n_vals
886 CALL reallocate(index_list, 1, new_size)
887 index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
888 old_size = new_size
889 END DO
890 ALLOCATE (work(new_size))
891 CALL sort(index_list, new_size, work)
892 DEALLOCATE (work)
893 DO i = 1, new_size
894 WRITE (unit=string, fmt="(A,I0,A)") "(", index_list(i), ")"
895 IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
896 WRITE (unit=iw, fmt="(T3,A)") &
897 "Invalid atomic index "//trim(string)//" specified. Print request is ignored."
898 cycle
899 END IF
900 IF (i > 1) THEN
901 ! Skip redundant indices
902 IF (index_list(i) == index_list(i - 1)) cycle
903 END IF
904 WRITE (unit=iw, fmt="(T3,A,T20,A,3F13.6)") &
905 "r"//trim(string), "=", pbc(particle_set(index_list(i))%r(1:3), cell)*conv
906 END DO
907 DEALLOCATE (index_list)
908 END IF
909
910 ! Print the requested atomic position vectors in scaled coordinates
911 CALL section_vals_val_get(section_vals=section, &
912 keyword_name="POSITION_SCALED", &
913 n_rep_val=n_rep)
914 IF (n_rep > 0) THEN
915 WRITE (unit=iw, fmt="(/,T3,A,/)") &
916 "Position vectors s(i) of the atoms i in scaled coordinates"
917 old_size = 0
918 DO i_rep = 1, n_rep
919 CALL section_vals_val_get(section_vals=section, &
920 keyword_name="POSITION_SCALED", &
921 i_rep_val=i_rep, &
922 i_vals=atomic_indices)
923 n_vals = SIZE(atomic_indices)
924 new_size = old_size + n_vals
925 CALL reallocate(index_list, 1, new_size)
926 index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
927 old_size = new_size
928 END DO
929 ALLOCATE (work(new_size))
930 CALL sort(index_list, new_size, work)
931 DEALLOCATE (work)
932 DO i = 1, new_size
933 WRITE (unit=string, fmt="(A,I0,A)") "(", index_list(i), ")"
934 IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
935 WRITE (unit=iw, fmt="(T3,A)") &
936 "Invalid atomic index "//trim(string)//" specified. Print request is ignored."
937 cycle
938 END IF
939 IF (i > 1) THEN
940 ! Skip redundant indices
941 IF (index_list(i) == index_list(i - 1)) cycle
942 END IF
943 r(1:3) = pbc(particle_set(index_list(i))%r(1:3), cell)
944 CALL real_to_scaled(s, r, cell)
945 WRITE (unit=iw, fmt="(T3,A,T20,A,3F13.6)") &
946 "s"//trim(string), "=", s(1:3)
947 END DO
948 DEALLOCATE (index_list)
949 END IF
950
951 ! Print the requested distances
952 CALL section_vals_val_get(section_vals=section, &
953 keyword_name="DISTANCE", &
954 n_rep_val=n)
955 IF (n > 0) THEN
956 WRITE (unit=iw, fmt="(/,T3,A,/)") &
957 "Distance vector r(i,j) between the atom i and j in "// &
958 trim(unit_str)
959 DO i = 1, n
960 CALL section_vals_val_get(section_vals=section, &
961 keyword_name="DISTANCE", &
962 i_rep_val=i, &
963 i_vals=atomic_indices)
964 string = ""
965 WRITE (unit=string, fmt="(A,2(I0,A))") &
966 "(", atomic_indices(1), ",", atomic_indices(2), ")"
967 wrk2 = atomic_indices
968 CALL sort_unique(wrk2, unique)
969 IF (((wrk2(1) >= 1) .AND. (wrk2(SIZE(wrk2)) <= natom)) .AND. unique) THEN
970 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
971 particle_set(atomic_indices(2))%r(:), cell)
972 dab = sqrt(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
973 WRITE (unit=iw, fmt="(T3,A,T20,A,3F13.6,3X,A,F13.6)") &
974 "r"//trim(string), "=", rab(:)*conv, &
975 "|r| =", dab*conv
976 ELSE
977 WRITE (unit=iw, fmt="(T3,A)") &
978 "Invalid atomic indices "//trim(string)//" specified. Print request is ignored."
979 END IF
980 END DO
981 END IF
982
983 ! Print the requested angles
984 CALL section_vals_val_get(section_vals=section, &
985 keyword_name="ANGLE", &
986 n_rep_val=n)
987 IF (n > 0) THEN
988 WRITE (unit=iw, fmt="(/,T3,A,/)") &
989 "Angle a(i,j,k) between the atomic distance vectors r(j,i) and "// &
990 "r(j,k) in DEGREE"
991 DO i = 1, n
992 CALL section_vals_val_get(section_vals=section, &
993 keyword_name="ANGLE", &
994 i_rep_val=i, &
995 i_vals=atomic_indices)
996 string = ""
997 WRITE (unit=string, fmt="(A,3(I0,A))") &
998 "(", atomic_indices(1), ",", atomic_indices(2), ",", atomic_indices(3), ")"
999 wrk3 = atomic_indices
1000 CALL sort_unique(wrk3, unique)
1001 IF (((wrk3(1) >= 1) .AND. (wrk3(SIZE(wrk3)) <= natom)) .AND. unique) THEN
1002 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1003 particle_set(atomic_indices(2))%r(:), cell)
1004 rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1005 particle_set(atomic_indices(3))%r(:), cell)
1006 WRITE (unit=iw, fmt="(T3,A,T26,A,F9.3)") &
1007 "a"//trim(string), "=", angle(-rab, rbc)*degree
1008 ELSE
1009 WRITE (unit=iw, fmt="(T3,A)") &
1010 "Invalid atomic indices "//trim(string)//" specified. Print request is ignored."
1011 END IF
1012 END DO
1013 END IF
1014
1015 ! Print the requested dihedral angles
1016 CALL section_vals_val_get(section_vals=section, &
1017 keyword_name="DIHEDRAL_ANGLE", &
1018 n_rep_val=n)
1019 IF (n > 0) THEN
1020 WRITE (unit=iw, fmt="(/,T3,A,/)") &
1021 "Dihedral angle d(i,j,k,l) between the planes (i,j,k) and (j,k,l) "// &
1022 "in DEGREE"
1023 DO i = 1, n
1024 CALL section_vals_val_get(section_vals=section, &
1025 keyword_name="DIHEDRAL_ANGLE", &
1026 i_rep_val=i, &
1027 i_vals=atomic_indices)
1028 string = ""
1029 WRITE (unit=string, fmt="(A,4(I0,A))") &
1030 "(", atomic_indices(1), ",", atomic_indices(2), ",", &
1031 atomic_indices(3), ",", atomic_indices(4), ")"
1032 wrk4 = atomic_indices
1033 CALL sort_unique(wrk4, unique)
1034 IF (((wrk4(1) >= 1) .AND. (wrk4(SIZE(wrk4)) <= natom)) .AND. unique) THEN
1035 rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1036 particle_set(atomic_indices(2))%r(:), cell)
1037 rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1038 particle_set(atomic_indices(3))%r(:), cell)
1039 rcd(:) = pbc(particle_set(atomic_indices(3))%r(:), &
1040 particle_set(atomic_indices(4))%r(:), cell)
1041 ! Match the sign convention of the TORSION collective variable.
1042 WRITE (unit=iw, fmt="(T3,A,T26,A,F9.3)") &
1043 "d"//trim(string), "=", -dihedral_angle(rab, rbc, rcd)*degree
1044 ELSE
1045 WRITE (unit=iw, fmt="(T3,A)") &
1046 "Invalid atomic indices "//trim(string)//" specified. Print request is ignored."
1047 END IF
1048 END DO
1049 END IF
1050 END IF
1051 CALL cp_print_key_finished_output(iw, logger, input_section, &
1052 "PRINT%STRUCTURE_DATA")
1053
1054 CALL timestop(handle)
1055
1056 END SUBROUTINE write_structure_data
1057
1058! **************************************************************************************************
1059!> \brief Write the final geometry and cell information to files
1060!> \param particle_set pointer to particles with atm_name, element_symbol and position
1061!> \param cell pointer to cell with abc, angle_alpha, angle_beta, angle_gamma and deth
1062!> \param input_section pointer to motion_section which has PRINT%FINAL_STRUCTURE
1063!> \param conv flag for whether convergence is achieved or not in optimization
1064!> \param keep_angles flag for whether cell optimization keeps initial angles
1065!> \param keep_symmetry flag for whether cell optimization keeps initial symmetry
1066!> \param keep_volume flag for whether cell optimization keeps initial volume
1067!> \param gopt_env_label the geometry optimization label "GEO_OPT", "CELL_OPT", ...
1068!> \param constraint_label label for directions with constraint in cell optimization
1069!> \par Intended to be invoked in gopt_f_methods:write_final_info.
1070!> This implementation does not consider higher space groups even if
1071!> one is detected, and the chemical formulae are neither written in
1072!> the sorted "Hill notation" nor expressed in groups of molecules.
1073!> Other potentially useful but yet to be written information includes:
1074!> the external pressure from CELL_OPT/EXTERNAL_POTENTIAL and the
1075!> stress tensor (virial) for CELL_OPT;
1076!> the fixed atoms from MOTION/CONSTRAINT/FIXED_ATOMS for all.
1077!>
1078!> History
1079!> 04.2026 - Created as write_final_cif
1080!> 05.2026 - Generalized to write_final_structure and enable extxyz
1081!> 06.2026 - Adopted write_particle_coordinates for handling kind in extxyz
1082!> \author HE Zilong
1083!> \version 1.0
1084! **************************************************************************************************
1085 SUBROUTINE write_final_structure(particle_set, cell, input_section, conv, &
1086 keep_angles, keep_symmetry, keep_volume, &
1087 gopt_env_label, constraint_label)
1088 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1089 TYPE(cell_type), INTENT(IN), POINTER :: cell
1090 TYPE(section_vals_type), INTENT(IN), POINTER :: input_section
1091 LOGICAL, INTENT(IN) :: conv, keep_angles, keep_symmetry, &
1092 keep_volume
1093 CHARACTER(LEN=default_string_length), INTENT(IN) :: gopt_env_label
1094 CHARACTER(LEN=4), INTENT(IN) :: constraint_label
1095
1096 CHARACTER(len=*), PARAMETER :: routinen = 'write_final_structure'
1097
1098 CHARACTER(LEN=1) :: conv_str
1099 CHARACTER(LEN=2) :: element_symbol
1100 CHARACTER(LEN=2), ALLOCATABLE :: element_list(:)
1101 CHARACTER(LEN=5) :: pbc_str
1102 CHARACTER(LEN=:), ALLOCATABLE :: formula_structural, formula_sum
1103 CHARACTER(LEN=default_path_length) :: cell_str, record, title
1104 CHARACTER(LEN=default_string_length) :: atm_name, f_cif, f_cif_label, &
1105 f_cif_type_symbol
1106 CHARACTER(LEN=default_string_length), ALLOCATABLE :: cif_label(:), cif_type_symbol(:)
1107 CHARACTER(LEN=timestamp_length) :: timestamp
1108 INTEGER :: elem_seen, file_unit, gcd_all, handle, i, iatom, ielem, natom, output_unit, &
1109 symmetry_id, w_cif_label, w_cif_type_symbol
1110 INTEGER, ALLOCATABLE :: count_list(:)
1111 LOGICAL :: dummy, elem_in_list, orthorhombic, &
1112 print_final_structure, print_kind, &
1113 write_cif, write_xyz
1114 REAL(kind=dp) :: angle_alpha, angle_beta, angle_gamma, &
1115 deth, unit_conv
1116 REAL(kind=dp), DIMENSION(3) :: abc, r, s
1117 REAL(kind=dp), DIMENSION(3, 3) :: hmat
1118 TYPE(cp_logger_type), POINTER :: logger
1119 TYPE(enumeration_type), POINTER :: enum
1120 TYPE(keyword_type), POINTER :: symmetry_keyword
1121 TYPE(section_type), POINTER :: tmp_cell_section
1122 TYPE(section_vals_type), POINTER :: print_key
1123
1124 CALL timeset(routinen, handle)
1125
1126 NULLIFY (enum, logger, symmetry_keyword, print_key, tmp_cell_section)
1127 logger => cp_get_default_logger()
1128 output_unit = cp_logger_get_default_io_unit(logger)
1129 print_key => section_vals_get_subs_vals(input_section, "PRINT%FINAL_STRUCTURE")
1130 print_final_structure = btest(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
1131 file_unit = -1
1132 conv_str = "F"
1133 IF (conv) conv_str(1:1) = "T"
1134
1135 ! Collect cell information
1136 pbc_str = "F F F"
1137 CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, &
1138 deth=deth, orthorhombic=orthorhombic, abc=abc, h=hmat, &
1139 symmetry_id=symmetry_id)
1140 IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
1141 IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
1142 IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
1143 CALL create_cell_section(tmp_cell_section)
1144 symmetry_keyword => section_get_keyword(tmp_cell_section, "SYMMETRY")
1145 CALL keyword_get(symmetry_keyword, enum=enum)
1146 ! cell_str is default_path_length which is longer
1147 ! than default_string_length and should be enough
1148 WRITE (unit=cell_str, fmt="(9(1X,F19.10))") &
1149 cp_unit_from_cp2k(hmat(1, 1), "angstrom"), &
1150 cp_unit_from_cp2k(hmat(2, 1), "angstrom"), &
1151 cp_unit_from_cp2k(hmat(3, 1), "angstrom"), &
1152 cp_unit_from_cp2k(hmat(1, 2), "angstrom"), &
1153 cp_unit_from_cp2k(hmat(2, 2), "angstrom"), &
1154 cp_unit_from_cp2k(hmat(3, 2), "angstrom"), &
1155 cp_unit_from_cp2k(hmat(1, 3), "angstrom"), &
1156 cp_unit_from_cp2k(hmat(2, 3), "angstrom"), &
1157 cp_unit_from_cp2k(hmat(3, 3), "angstrom")
1158
1159 ! Collect atom information
1160 natom = SIZE(particle_set)
1161 ALLOCATE (element_list(nelem + 1), count_list(nelem + 1))
1162 count_list(:) = 0
1163 ALLOCATE (cif_label(natom), cif_type_symbol(natom))
1164 elem_seen = 0
1165 w_cif_type_symbol = 0
1166 w_cif_label = 0
1167 atom_loop: DO iatom = 1, natom
1168 elem_in_list = .false.
1169 CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1170 name=atm_name, element_symbol=element_symbol)
1171 cif_type_symbol(iatom) = trim(atm_name)
1172 ! From write_particle_coordinates above it seems possible
1173 ! for some atoms to have empty element symbols; whatever
1174 ! these are, do not count them in the chemical formula
1175 IF (len_trim(element_symbol) == 0) THEN
1176 dummy = qmmm_ff_precond_only_qm(id1=atm_name)
1177 cif_label(iatom) = trim(atm_name)//trim(adjustl(cp_to_string(iatom)))
1178 ELSE
1179 cif_label(iatom) = trim(element_symbol)//trim(adjustl(cp_to_string(iatom)))
1180 elem_loop: DO ielem = 1, elem_seen
1181 IF (element_list(ielem) == element_symbol) THEN
1182 elem_in_list = .true.
1183 count_list(ielem) = count_list(ielem) + 1
1184 EXIT elem_loop
1185 END IF
1186 END DO elem_loop
1187 IF (.NOT. elem_in_list) THEN
1188 elem_seen = elem_seen + 1
1189 element_list(elem_seen) = element_symbol
1190 count_list(elem_seen) = 1
1191 END IF
1192 END IF
1193 IF (len_trim(cif_type_symbol(iatom)) > w_cif_type_symbol) THEN
1194 w_cif_type_symbol = len_trim(cif_type_symbol(iatom))
1195 END IF
1196 IF (len_trim(cif_label(iatom)) > w_cif_label) THEN
1197 w_cif_label = len_trim(cif_label(iatom))
1198 END IF
1199 END DO atom_loop
1200
1201 ! Determine the format of each line in cif considering width of cif_type_symbol and cif_label
1202 ! The fields are, in order:
1203 ! _atom_site_type_symbol, _atom_site_label, _atom_site_symmetry_multiplicity,
1204 ! _atom_site_fract_x, _atom_site_fract_y, _atom_site_fract_z, _atom_site_occupancy
1205 ! in which:
1206 ! _atom_site_type_symbol is taken as atm_name
1207 ! _atom_site_label is taken as element_symbol//iatom
1208 ! _atom_site_symmetry_multiplicity and _atom_site_occupancy are always 1
1209 f_cif_type_symbol = "A"//trim(adjustl(cp_to_string(w_cif_type_symbol + 4)))
1210 f_cif_label = "A"//trim(adjustl(cp_to_string(w_cif_label + 4)))
1211 f_cif = "(T3,"//trim(f_cif_type_symbol)//","//trim(f_cif_label)//",I4,3F14.8,F8.2)"
1212
1213 ! Determine formula_sum
1214 cpassert(elem_seen > 0)
1215 cpassert(count_list(1) > 0)
1216 formula_sum = "'"
1217 DO ielem = 1, elem_seen
1218 formula_sum = formula_sum//trim(adjustl(element_list(ielem)))
1219 formula_sum = formula_sum//trim(adjustl(cp_to_string(count_list(ielem))))
1220 formula_sum = formula_sum//" "
1221 END DO
1222 formula_sum = trim(adjustl(formula_sum))//"'"
1223
1224 ! Determine formula_structural and Z
1225 gcd_all = count_list(1)
1226 DO ielem = 1, elem_seen
1227 IF (count_list(ielem) /= 0) THEN
1228 gcd_all = gcd(gcd_all, count_list(ielem))
1229 END IF
1230 END DO
1231 IF (gcd_all > 1) count_list = count_list/gcd_all
1232 formula_structural = "'"
1233 DO ielem = 1, elem_seen
1234 formula_structural = formula_structural//trim(adjustl(element_list(ielem)))
1235 formula_structural = formula_structural//trim(adjustl(cp_to_string(count_list(ielem))))
1236 formula_structural = formula_structural//" "
1237 END DO
1238 formula_structural = trim(adjustl(formula_structural))//"'"
1239
1240 ! Write XYZ
1241 CALL section_vals_val_get(print_key, "PRINT_XYZ", l_val=write_xyz)
1242 write_xyz = write_xyz .AND. print_final_structure
1243 CALL section_vals_val_get(print_key, "PRINT_ATOM_KIND", l_val=print_kind)
1244 IF (write_xyz) THEN
1245 ! Print a message to log
1246 record = cp_print_key_generate_filename(logger, print_key, &
1247 extension=".xyz", &
1248 my_local=.false.)
1249 IF (output_unit > 0) THEN
1250 IF (conv) THEN
1251 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1252 routinen//": Optimization converged, writing XYZ file gladly:"
1253 ELSE
1254 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1255 routinen//": Optimization not yet converged, writing XYZ file anyway:"
1256 END IF
1257 WRITE (unit=output_unit, fmt="(T3,A)") trim(record)
1258 END IF
1259
1260 ! Prepare title
1261 WRITE (unit=title, fmt="(A)") &
1262 'Lattice="'//trim(adjustl(cell_str))//'" '// &
1263 'Properties=species:S:1:pos:R:3 '// &
1264 'pbc="'//pbc_str//'" '// &
1265 'Converged='//conv_str
1266 ! Extended XYZ uses angstrom for positions
1267 unit_conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1268 ! Prepare file unit and write to it
1269 file_unit = cp_print_key_unit_nr(logger, input_section, "PRINT%FINAL_STRUCTURE", &
1270 file_status="REPLACE", file_form="FORMATTED", &
1271 extension=".xyz")
1272 IF (file_unit > 0) THEN
1273 CALL write_particle_coordinates(particle_set, file_unit, dump_extxyz, "POS", title, &
1274 cell=cell, unit_conv=unit_conv, print_kind=print_kind)
1275 END IF
1276 END IF
1277
1278 ! Write CIF
1279 CALL section_vals_val_get(print_key, "PRINT_CIF", l_val=write_cif)
1280 write_cif = write_cif .AND. print_final_structure
1281 IF (write_cif) THEN
1282 ! Print a message to log
1283 record = cp_print_key_generate_filename(logger, print_key, &
1284 extension=".cif", &
1285 my_local=.false.)
1286 IF (output_unit > 0) THEN
1287 IF (conv) THEN
1288 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1289 routinen//": Optimization converged, writing CIF file gladly:"
1290 ELSE
1291 WRITE (unit=output_unit, fmt="(/,T2,A)") &
1292 routinen//": Optimization not yet converged, writing CIF file anyway:"
1293 END IF
1294 WRITE (unit=output_unit, fmt="(T3,A)") trim(record)
1295 END IF
1296
1297 ! Make timestamp for the file
1298 CALL m_timestamp(timestamp)
1299
1300 ! Prepare file unit and write to it
1301 file_unit = cp_print_key_unit_nr(logger, input_section, "PRINT%FINAL_STRUCTURE", &
1302 file_status="REPLACE", file_form="FORMATTED", &
1303 extension=".cif")
1304 IF (file_unit > 0) THEN
1305 ! Generic information
1306 WRITE (unit=file_unit, fmt="(A)") &
1307 "# CIF file created by CP2K "//trim(modulen)//":"//trim(routinen)
1308 WRITE (unit=file_unit, fmt="(A)") &
1309 "data_"//trim(logger%iter_info%project_name)
1310 WRITE (unit=file_unit, fmt="(A,T39,A)") &
1311 "_audit_creation_date", timestamp(:10)
1312 WRITE (unit=file_unit, fmt="(A,/,A,/,A)") &
1313 "_audit_creation_method", ";", &
1314 trim(cp2k_version)//" (revision "//trim(compile_revision)//")"
1315 WRITE (unit=file_unit, fmt="(A,/,A,/,A,/,A)") &
1316 "Project name "//trim(logger%iter_info%project_name), &
1317 "submitted by "//trim(r_user_name)//"@"//trim(r_host_name), &
1318 "processed in "//trim(r_cwd), &
1319 "generated at "//trim(timestamp)
1320 WRITE (unit=file_unit, fmt="(T2,A)") &
1321 repeat(" -*WARNING*- ", 6)
1322 WRITE (unit=file_unit, fmt="(T2,A)") &
1323 "WARNING| the lines enclosed in between -*WARNING*- contain metadata"
1324 WRITE (unit=file_unit, fmt="(T2,A)") &
1325 "WARNING| provided in a tentative format for the newly implemented"
1326 WRITE (unit=file_unit, fmt="(T2,A)") &
1327 "WARNING| CIF output from CP2K version 2026.2, which may be unstable"
1328 WRITE (unit=file_unit, fmt="(T2,A)") &
1329 "WARNING| and subject to overhaul in a future release. Contact the"
1330 WRITE (unit=file_unit, fmt="(T2,A)") &
1331 "WARNING| developers in case a stabilized format is needed for some"
1332 WRITE (unit=file_unit, fmt="(T2,A)") &
1333 "WARNING| downstream CIF parser with text pattern matchers (regexp)."
1334 WRITE (unit=file_unit, fmt="(T2,A)") &
1335 "- Optimization type: "//trim(gopt_env_label)
1336 IF (conv) THEN
1337 WRITE (unit=file_unit, fmt="(T2,A)") &
1338 "- Optimization converged: TRUE"
1339 ELSE
1340 WRITE (unit=file_unit, fmt="(T2,A)") &
1341 "- Optimization converged: FALSE"
1342 END IF
1343 WRITE (unit=file_unit, fmt="(T2,A)") &
1344 "- Requested initial cell symmetry: "//trim(enum_i2c(enum, symmetry_id))
1345 IF (orthorhombic) THEN
1346 WRITE (unit=file_unit, fmt="(T2,A)") &
1347 "- Cell is numerically orthorhombic: TRUE"
1348 ELSE
1349 WRITE (unit=file_unit, fmt="(T2,A)") &
1350 "- Cell is numerically orthorhombic: FALSE"
1351 END IF
1352 WRITE (unit=file_unit, fmt="(T2,A)") &
1353 "- Periodicity of cell: "//trim(pbc_str)
1354 IF (gopt_env_label == "CELL_OPT") THEN
1355 WRITE (unit=file_unit, fmt="(T2,A)") &
1356 "- Cell is subject to optimization: TRUE"
1357 WRITE (unit=file_unit, fmt="(T2,A)") &
1358 "- Cell has constraint on direction: "//trim(adjustl(constraint_label))
1359 IF (keep_angles) THEN
1360 WRITE (unit=file_unit, fmt="(T2,A)") &
1361 "- Keep angles between the cell vectors during optimization: TRUE"
1362 ELSE
1363 WRITE (unit=file_unit, fmt="(T2,A)") &
1364 "- Keep angles between the cell vectors during optimization: FALSE"
1365 END IF
1366 IF (keep_symmetry) THEN
1367 WRITE (unit=file_unit, fmt="(T2,A)") &
1368 "- Keep initial cell symmetry during optimization: TRUE"
1369 ELSE
1370 WRITE (unit=file_unit, fmt="(T2,A)") &
1371 "- Keep initial cell symmetry during optimization: FALSE"
1372 END IF
1373 IF (keep_volume) THEN
1374 WRITE (unit=file_unit, fmt="(T2,A)") &
1375 "- Keep initial cell volume during optimization: TRUE"
1376 ELSE
1377 WRITE (unit=file_unit, fmt="(T2,A)") &
1378 "- Keep initial cell volume during optimization: FALSE"
1379 END IF
1380 ELSE
1381 WRITE (unit=file_unit, fmt="(T2,A)") &
1382 "- Cell is subject to optimization: FALSE"
1383 END IF
1384 WRITE (unit=file_unit, fmt="(T2,A)") &
1385 "- Final cell vectors A, B, C by rows [angstrom]:"
1386 DO i = 1, 3
1387 WRITE (unit=file_unit, fmt="(T3,3(1X,F19.10))") &
1388 cp_unit_from_cp2k(hmat(1, i), "angstrom"), &
1389 cp_unit_from_cp2k(hmat(2, i), "angstrom"), &
1390 cp_unit_from_cp2k(hmat(3, i), "angstrom")
1391 END DO
1392 WRITE (unit=file_unit, fmt="(T2,A)") &
1393 repeat("-*WARNING*- ", 6)
1394 WRITE (unit=file_unit, fmt="(A)") ";"
1395 ! Data of cell and geometry
1396 WRITE (unit=file_unit, fmt="(/,A,T44,A)") &
1397 "_symmetry_space_group_name_H-M", "'P 1'"
1398 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1399 "_cell_length_a", cp_unit_from_cp2k(abc(1), "angstrom")
1400 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1401 "_cell_length_b", cp_unit_from_cp2k(abc(2), "angstrom")
1402 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1403 "_cell_length_c", cp_unit_from_cp2k(abc(3), "angstrom")
1404 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1405 "_cell_angle_alpha", angle_alpha
1406 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1407 "_cell_angle_beta", angle_beta
1408 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1409 "_cell_angle_gamma", angle_gamma
1410 WRITE (unit=file_unit, fmt="(A,T48,A)") &
1411 "_symmetry_Int_Tables_number", "1"
1412 WRITE (unit=file_unit, fmt="(A,T36,A)") &
1413 "_chemical_formula_structural", formula_structural
1414 WRITE (unit=file_unit, fmt="(A,T36,A)") &
1415 "_chemical_formula_sum", formula_sum
1416 WRITE (unit=file_unit, fmt="(A,T31,F18.8)") &
1417 "_cell_volume", cp_unit_from_cp2k(abs(deth), "angstrom^3")
1418 WRITE (unit=file_unit, fmt="(A,T41,I8)") &
1419 "_cell_formula_units_Z", gcd_all
1420 WRITE (unit=file_unit, fmt="(A,/,T2,A,/,T2,A,/,T3,A)") &
1421 "loop_", "_symmetry_equiv_pos_site_id", &
1422 "_symmetry_equiv_pos_as_xyz", "1 'x, y, z'"
1423 WRITE (unit=file_unit, fmt="(A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A)") &
1424 "loop_", "_atom_site_type_symbol", "_atom_site_label", &
1425 "_atom_site_symmetry_multiplicity", "_atom_site_fract_x", &
1426 "_atom_site_fract_y", "_atom_site_fract_z", "_atom_site_occupancy"
1427 DO iatom = 1, natom
1428 ! positive_range=.TRUE. makes r(1:3) in [0, cell%hmat(i,i)] and
1429 ! s(1:3) in [0, 1], so there is no need to MODULO s(1:3) by 1.0
1430 r(1:3) = pbc(particle_set(iatom)%r(1:3), cell, positive_range=.true.)
1431 CALL real_to_scaled(s, r, cell)
1432 WRITE (unit=file_unit, fmt=trim(f_cif)) &
1433 cif_type_symbol(iatom), cif_label(iatom), 1, s(1:3), 1.0_dp
1434 END DO
1435 END IF
1436 END IF
1437
1438 ! Finish
1439 DEALLOCATE (element_list, count_list, formula_structural, &
1440 formula_sum, cif_label, cif_type_symbol)
1441 CALL section_release(tmp_cell_section)
1442 CALL cp_print_key_finished_output(file_unit, logger, input_section, &
1443 "PRINT%FINAL_STRUCTURE")
1444 IF (output_unit > 0 .AND. (write_xyz .OR. write_cif)) THEN
1445 WRITE (unit=output_unit, fmt='(/,T2,A)') &
1446 routinen//": Done!"
1447 END IF
1448
1449 CALL timestop(handle)
1450
1451 END SUBROUTINE write_final_structure
1452
1453END MODULE particle_methods
modulo
static GRID_HOST_DEVICE int modulo(int a, int m)
Equivalent of Fortran's MODULO, which always return a positive number. https://gcc....
Definition grid_common.h:125
cell_types::pbc
Definition cell_types.F:104
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
external_potential_types::get_potential
Definition external_potential_types.F:278
memory_utilities::reallocate
Definition memory_utilities.F:27
util::sort_unique
Definition util.F:36
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_methods
Handles all functions related to the CELL.
Definition cell_methods.F:15
cell_methods::set_cell_param
subroutine, public set_cell_param(cell, cell_length, cell_angle, periodic, do_init_cell)
Sets the cell using the internal parameters (a,b,c) (alpha,beta,gamma) using the convention: a parall...
Definition cell_methods.F:672
cell_methods::cell_create
subroutine, public cell_create(cell, hmat, periodic, tag)
allocates and initializes a cell
Definition cell_methods.F:95
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cell_types::real_to_scaled
subroutine, public real_to_scaled(s, r, cell)
Transform real to scaled cell coordinates. s=h_inv*r.
Definition cell_types.F:595
cell_types::cell_release
subroutine, public cell_release(cell)
releases the given cell (see doc/ReferenceCounting.html)
Definition cell_types.F:668
cell_types::cell_clone
subroutine, public cell_clone(cell_in, cell_out, tag)
Clone cell variable.
Definition cell_types.F:141
cell_types::get_cell
subroutine, public get_cell(cell, alpha, beta, gamma, deth, orthorhombic, abc, periodic, h, h_inv, symmetry_id, tag)
Get informations about a simulation cell.
Definition cell_types.F:233
cp2k_info
some minimal info about CP2K, including its version and license
Definition cp2k_info.F:22
cp2k_info::r_host_name
character(len=default_string_length), public r_host_name
Definition cp2k_info.F:74
cp2k_info::compile_revision
character(len= *), parameter, public compile_revision
Definition cp2k_info.F:45
cp2k_info::cp2k_version
character(len= *), parameter, public cp2k_version
Definition cp2k_info.F:49
cp2k_info::r_cwd
character(len=default_path_length), public r_cwd
Definition cp2k_info.F:73
cp2k_info::r_user_name
character(len=default_string_length), public r_user_name
Definition cp2k_info.F:74
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_generate_filename
character(len=default_path_length) function, public cp_print_key_generate_filename(logger, print_key, middle_name, extension, my_local)
Utility function that returns a unit number to write the print key. Might open a file with a unique f...
Definition cp_output_handling.F:754
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1073
cp_output_handling::cp_p_file
integer, parameter, public cp_p_file
Definition cp_output_handling.F:103
cp_output_handling::cp_print_key_should_output
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
Definition cp_output_handling.F:345
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_from_cp2k
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition cp_units.F:1251
external_potential_types
Definition of the atomic potential types.
Definition external_potential_types.F:14
gamma
Calculation of the incomplete Gamma function F_n(t) for multi-center integrals over Cartesian Gaussia...
Definition gamma.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:18
input_constants::dump_xmol
integer, parameter, public dump_xmol
Definition input_constants.F:389
input_constants::dump_pdb
integer, parameter, public dump_pdb
Definition input_constants.F:389
input_constants::dump_extxyz
integer, parameter, public dump_extxyz
Definition input_constants.F:389
input_constants::dump_atomic
integer, parameter, public dump_atomic
Definition input_constants.F:389
input_constants::dump_dcd_aligned_cell
integer, parameter, public dump_dcd_aligned_cell
Definition input_constants.F:389
input_constants::dump_dcd
integer, parameter, public dump_dcd
Definition input_constants.F:389
input_cp2k_subsys
builds the subsystem section of the input
Definition input_cp2k_subsys.F:14
input_cp2k_subsys::create_cell_section
subroutine, public create_cell_section(section, periodic)
creates the cell section
Definition input_cp2k_subsys.F:75
input_enumeration_types
represents an enumeration, i.e. a mapping between integers and strings
Definition input_enumeration_types.F:14
input_enumeration_types::enum_i2c
character(len=default_string_length) function, public enum_i2c(enum, i)
maps an integer to a string
Definition input_enumeration_types.F:152
input_keyword_types
represents keywords in an input
Definition input_keyword_types.F:14
input_keyword_types::keyword_get
subroutine, public keyword_get(keyword, names, usage, description, type_of_var, n_var, default_value, lone_keyword_value, repeats, enum, citations)
...
Definition input_keyword_types.F:489
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_release
recursive subroutine, public section_release(section)
releases the given keyword list (see doc/ReferenceCounting.html)
Definition input_section_types.F:238
input_section_types::section_get_keyword
recursive type(keyword_type) function, pointer, public section_get_keyword(section, keyword_name)
returns the requested keyword
Definition input_section_types.F:481
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
kinds::sp
integer, parameter, public sp
Definition kinds.F:33
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::timestamp_length
integer, parameter, public timestamp_length
Definition machine.F:46
machine::m_timestamp
subroutine, public m_timestamp(timestamp)
Returns a human readable timestamp.
Definition machine.F:381
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::degree
real(kind=dp), parameter, public degree
Definition mathconstants.F:139
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::angle
pure real(kind=dp) function, public angle(a, b)
Calculation of the angle between the vectors a and b. The angle is returned in radians.
Definition mathlib.F:184
mathlib::gcd
elemental integer function, public gcd(a, b)
computes the greatest common divisor of two number
Definition mathlib.F:1289
mathlib::dihedral_angle
pure real(kind=dp) function, public dihedral_angle(ab, bc, cd)
Returns the dihedral angle, i.e. the angle between the planes defined by the vectors (-ab,...
Definition mathlib.F:476
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
particle_methods
Define methods related to particle_type.
Definition particle_methods.F:14
particle_methods::write_qs_particle_coordinates
subroutine, public write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
Write the atomic coordinates to the output unit.
Definition particle_methods.F:591
particle_methods::write_fist_particle_coordinates
subroutine, public write_fist_particle_coordinates(particle_set, subsys_section, charges)
Write the atomic coordinates to the output unit.
Definition particle_methods.F:529
particle_methods::write_final_structure
subroutine, public write_final_structure(particle_set, cell, input_section, conv, keep_angles, keep_symmetry, keep_volume, gopt_env_label, constraint_label)
Write the final geometry and cell information to files.
Definition particle_methods.F:1088
particle_methods::write_particle_matrix
subroutine, public write_particle_matrix(matrix, particle_set, iw, el_per_part, ilist, parts_per_line)
...
Definition particle_methods.F:761
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis, ncgf)
Get the components of a particle set.
Definition particle_methods.F:123
particle_methods::write_structure_data
subroutine, public write_structure_data(particle_set, cell, input_section)
Write structure data requested by a separate structure data input section to the output unit....
Definition particle_methods.F:831
particle_methods::write_particle_distances
subroutine, public write_particle_distances(particle_set, cell, subsys_section)
Write the matrix of the particle distances to the output unit.
Definition particle_methods.F:651
particle_methods::write_particle_coordinates
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Definition particle_methods.F:228
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
particle_types::get_particle_pos_or_vel
pure real(kind=dp) function, dimension(3), public get_particle_pos_or_vel(iatom, particle_set, vector)
Return the atomic position or velocity of atom iatom in x from a packed vector even if core-shell par...
Definition particle_types.F:182
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::nelem
integer, parameter, public nelem
Definition periodic_table.F:63
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::massunit
real(kind=dp), parameter, public massunit
Definition physcon.F:141
qmmm_ff_fist
Definition qmmm_ff_fist.F:11
qmmm_ff_fist::qmmm_ff_precond_only_qm
logical function, public qmmm_ff_precond_only_qm(id1, id2, id3, id4, is_link)
This function handles the atom names and modifies the "_QM_" prefix, in order to find the parameters ...
Definition qmmm_ff_fist.F:39
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
shell_potential_types
Definition shell_potential_types.F:11
shell_potential_types::get_shell
elemental subroutine, public get_shell(shell, charge, charge_core, charge_shell, mass_core, mass_shell, k2_spring, k4_spring, max_dist, shell_cutoff)
...
Definition shell_potential_types.F:71
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
util
All kind of helpful little routines.
Definition util.F:14
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:94
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
external_potential_types::fist_potential_type
Definition external_potential_types.F:85
input_enumeration_types::enumeration_type
Definition input_enumeration_types.F:43
input_keyword_types::keyword_type
represent a keyword in the input
Definition input_keyword_types.F:89
input_section_types::section_type
represent a section of the input file
Definition input_section_types.F:102
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
particle_types::particle_type
Definition particle_types.F:35
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28