Line data Source code
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 : ! **************************************************************************************************
14 : MODULE 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 :
102 : CONTAINS
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 218214 : SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, &
122 218214 : 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 218214 : CPASSERT(ASSOCIATED(particle_set))
133 :
134 218214 : nparticle = SIZE(particle_set)
135 218214 : IF (PRESENT(first_sgf)) THEN
136 51933 : CPASSERT(SIZE(first_sgf) >= nparticle)
137 : END IF
138 218214 : IF (PRESENT(last_sgf)) THEN
139 41713 : CPASSERT(SIZE(last_sgf) >= nparticle)
140 : END IF
141 218214 : IF (PRESENT(nsgf)) THEN
142 165817 : CPASSERT(SIZE(nsgf) >= nparticle)
143 : END IF
144 218214 : IF (PRESENT(nmao)) THEN
145 14 : CPASSERT(SIZE(nmao) >= nparticle)
146 : END IF
147 218214 : IF (PRESENT(ncgf)) THEN
148 4 : CPASSERT(SIZE(ncgf) >= nparticle)
149 : END IF
150 :
151 218214 : IF (PRESENT(first_sgf) .OR. PRESENT(last_sgf) .OR. PRESENT(nsgf)) THEN
152 : isgf = 0
153 1191533 : DO iparticle = 1, nparticle
154 973831 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
155 973831 : IF (PRESENT(basis)) THEN
156 734273 : IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
157 734269 : CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
158 : ELSE
159 4 : ns = 0
160 : END IF
161 : ELSE
162 239558 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns)
163 : END IF
164 973831 : IF (PRESENT(nsgf)) nsgf(iparticle) = ns
165 973831 : IF (PRESENT(first_sgf)) first_sgf(iparticle) = isgf + 1
166 973831 : isgf = isgf + ns
167 2165876 : IF (PRESENT(last_sgf)) last_sgf(iparticle) = isgf
168 : END DO
169 : END IF
170 :
171 218214 : IF (PRESENT(ncgf)) THEN
172 12 : DO iparticle = 1, nparticle
173 8 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
174 8 : IF (PRESENT(basis)) THEN
175 8 : IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
176 8 : CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, ncgf=ns)
177 : ELSE
178 0 : ns = 0
179 : END IF
180 : ELSE
181 0 : CALL get_qs_kind(qs_kind_set(ikind), ncgf=ns)
182 : END IF
183 20 : ncgf(iparticle) = ns
184 : END DO
185 : END IF
186 :
187 218214 : IF (PRESENT(first_sgf)) THEN
188 51933 : IF (SIZE(first_sgf) > nparticle) first_sgf(nparticle + 1) = isgf + 1
189 : END IF
190 :
191 218214 : IF (PRESENT(nmao)) THEN
192 86 : DO iparticle = 1, nparticle
193 72 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
194 72 : CALL get_qs_kind(qs_kind_set(ikind), mao=ns)
195 86 : nmao(iparticle) = ns
196 : END DO
197 : END IF
198 :
199 218214 : 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 27597 : SUBROUTINE write_particle_coordinates(particle_set, iunit, output_format, &
225 27597 : 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 27597 : 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 27597 : 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 27597 : CALL timeset(routineN, handle)
259 :
260 27597 : natom = SIZE(particle_set)
261 27597 : IF (PRESENT(array)) THEN
262 1838 : SELECT CASE (TRIM(content))
263 : CASE ("POS_VEL", "POS_VEL_FORCE")
264 1838 : CPABORT("Illegal usage")
265 : END SELECT
266 : END IF
267 27597 : factor = 1.0_dp
268 27597 : IF (PRESENT(unit_conv)) THEN
269 27446 : factor = unit_conv
270 : END IF
271 55102 : SELECT CASE (output_format)
272 : CASE (dump_xmol, dump_extxyz)
273 27505 : my_print_kind = .FALSE.
274 27505 : IF (PRESENT(print_kind)) my_print_kind = print_kind
275 27505 : WRITE (iunit, "(I8)") natom
276 27505 : WRITE (iunit, "(A)") TRIM(title)
277 1384871 : DO iatom = 1, natom
278 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
279 1357366 : element_symbol=element_symbol)
280 1357366 : IF (LEN_TRIM(element_symbol) == 0 .OR. my_print_kind) THEN
281 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
282 24 : name=atm_name)
283 24 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
284 24 : my_format = "(A,"
285 24 : atm_name = TRIM(atm_name)
286 : ELSE
287 1357342 : my_format = "(T2,A2,"
288 1357342 : atm_name = TRIM(element_symbol)
289 : END IF
290 27505 : SELECT CASE (TRIM(content))
291 : CASE ("POS")
292 1242177 : IF (PRESENT(array)) THEN
293 53812 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
294 : ELSE
295 4753460 : r(:) = particle_set(iatom)%r(:)
296 : END IF
297 4968708 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), r(1:3)*factor
298 : CASE ("VEL")
299 85772 : IF (PRESENT(array)) THEN
300 0 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
301 : ELSE
302 343088 : v(:) = particle_set(iatom)%v(:)
303 : END IF
304 343088 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), v(1:3)*factor
305 : CASE ("FORCE")
306 20808 : IF (PRESENT(array)) THEN
307 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
308 : ELSE
309 83232 : f(:) = particle_set(iatom)%f(:)
310 : END IF
311 83232 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), f(1:3)*factor
312 : CASE ("FORCE_MIXING_LABELS")
313 8609 : IF (PRESENT(array)) THEN
314 34436 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
315 : ELSE
316 0 : f(:) = particle_set(iatom)%f(:)
317 : END IF
318 1391802 : 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 170 : DO iatom = 1, natom
323 10 : SELECT CASE (TRIM(content))
324 : CASE ("POS")
325 160 : IF (PRESENT(array)) THEN
326 0 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
327 : ELSE
328 640 : r(:) = particle_set(iatom)%r(:)
329 : END IF
330 640 : WRITE (iunit, "(3F20.10)") r(1:3)*factor
331 : CASE ("VEL")
332 0 : IF (PRESENT(array)) THEN
333 0 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
334 : ELSE
335 0 : v(:) = particle_set(iatom)%v(:)
336 : END IF
337 0 : WRITE (iunit, "(3F20.10)") v(1:3)*factor
338 : CASE ("FORCE")
339 0 : IF (PRESENT(array)) THEN
340 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
341 : ELSE
342 0 : f(:) = particle_set(iatom)%f(:)
343 : END IF
344 0 : WRITE (iunit, "(3F20.10)") f(1:3)*factor
345 : CASE ("FORCE_MIXING_LABELS")
346 0 : IF (PRESENT(array)) THEN
347 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
348 : ELSE
349 0 : f(:) = particle_set(iatom)%f(:)
350 : END IF
351 160 : WRITE (iunit, "(3F20.10)") f(1:3)*factor
352 : END SELECT
353 : END DO
354 : CASE (dump_dcd, dump_dcd_aligned_cell)
355 4 : IF (.NOT. (PRESENT(cell))) THEN
356 0 : 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 4 : abc=abc)
360 4 : 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 0 : NULLIFY (cell_dcd)
366 0 : CALL cell_create(cell_dcd)
367 0 : CALL cell_clone(cell, cell_dcd, tag="CELL_DCD")
368 0 : angles(1) = angle_alpha/degree
369 0 : angles(2) = angle_beta/degree
370 0 : angles(3) = angle_gamma/degree
371 : CALL set_cell_param(cell_dcd, abc, angles, &
372 0 : do_init_cell=.TRUE.)
373 0 : h(1:3, 1:3) = MATMUL(cell_dcd%hmat(1:3, 1:3), cell%h_inv(1:3, 1:3))
374 0 : CALL cell_release(cell_dcd)
375 : END IF
376 12 : ALLOCATE (arr(3, natom))
377 4 : IF (PRESENT(array)) THEN
378 0 : arr(1:3, 1:natom) = RESHAPE(array, [3, natom])
379 : ELSE
380 8 : SELECT CASE (TRIM(content))
381 : CASE ("POS")
382 1156 : DO iatom = 1, natom
383 4612 : arr(1:3, iatom) = particle_set(iatom)%r(1:3)
384 : END DO
385 : CASE ("VEL")
386 0 : DO iatom = 1, natom
387 0 : arr(1:3, iatom) = particle_set(iatom)%v(1:3)
388 : END DO
389 : CASE ("FORCE")
390 0 : DO iatom = 1, natom
391 0 : arr(1:3, iatom) = particle_set(iatom)%f(1:3)
392 : END DO
393 : CASE DEFAULT
394 4 : CPABORT("Illegal DCD dump type")
395 : END SELECT
396 : END IF
397 12 : ALLOCATE (x4(natom))
398 8 : ALLOCATE (y4(natom))
399 8 : ALLOCATE (z4(natom))
400 4 : IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN
401 0 : x4(1:natom) = REAL(MATMUL(h(1, 1:3), arr(1:3, 1:natom)), KIND=sp)
402 0 : y4(1:natom) = REAL(MATMUL(h(2, 1:3), arr(1:3, 1:natom)), KIND=sp)
403 0 : z4(1:natom) = REAL(MATMUL(h(3, 1:3), arr(1:3, 1:natom)), KIND=sp)
404 : ELSE
405 1156 : x4(1:natom) = REAL(arr(1, 1:natom), KIND=sp)
406 1156 : y4(1:natom) = REAL(arr(2, 1:natom), KIND=sp)
407 1156 : z4(1:natom) = REAL(arr(3, 1:natom), KIND=sp)
408 : END IF
409 4 : WRITE (iunit) abc(1)*factor, angle_gamma, abc(2)*factor, &
410 8 : angle_beta, angle_alpha, abc(3)*factor
411 1156 : WRITE (iunit) x4*REAL(factor, KIND=sp)
412 1156 : WRITE (iunit) y4*REAL(factor, KIND=sp)
413 1156 : WRITE (iunit) z4*REAL(factor, KIND=sp)
414 : ! Release work storage
415 4 : DEALLOCATE (arr)
416 4 : DEALLOCATE (x4)
417 4 : DEALLOCATE (y4)
418 8 : DEALLOCATE (z4)
419 : CASE (dump_pdb)
420 78 : my_charge_occup = .FALSE.
421 78 : IF (PRESENT(charge_occup)) my_charge_occup = charge_occup
422 78 : my_charge_beta = .FALSE.
423 78 : IF (PRESENT(charge_beta)) my_charge_beta = charge_beta
424 78 : my_charge_extended = .FALSE.
425 78 : IF (PRESENT(charge_extended)) my_charge_extended = charge_extended
426 78 : IF (LEN_TRIM(title) > 0) THEN
427 : WRITE (UNIT=iunit, FMT="(A6,T11,A)") &
428 78 : "REMARK", TRIM(title)
429 : END IF
430 78 : 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 312 : "CRYST1", abc(1:3)*factor, angle_alpha, angle_beta, angle_gamma
444 78 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
445 3588 : DO iatom = 1, natom
446 3510 : 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 3510 : fist_potential=fist_potential, shell=shell)
467 3510 : IF (LEN_TRIM(element_symbol) == 0) THEN
468 0 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
469 : END IF
470 3510 : name = TRIM(atm_name)
471 3510 : IF (ASSOCIATED(fist_potential)) THEN
472 3510 : CALL get_potential(potential=fist_potential, qeff=qeff)
473 : ELSE
474 0 : qeff = 0.0_dp
475 : END IF
476 3510 : IF (ASSOCIATED(shell)) CALL get_shell(shell=shell, charge=qeff)
477 3510 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
478 3510 : WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(iatom, 100000)
479 3510 : 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 7020 : SELECT CASE (TRIM(content))
483 : CASE ("POS")
484 3510 : IF (PRESENT(array)) THEN
485 0 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
486 : ELSE
487 14040 : r(:) = particle_set(iatom)%r(:)
488 : END IF
489 14040 : WRITE (UNIT=line(31:54), FMT="(3F8.3)") r(1:3)*factor
490 : CASE DEFAULT
491 3510 : CPABORT("PDB dump only for trajectory available")
492 : END SELECT
493 3510 : IF (my_charge_occup) THEN
494 2700 : WRITE (UNIT=line(55:60), FMT="(F6.2)") qeff
495 : ELSE
496 810 : WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp
497 : END IF
498 3510 : IF (my_charge_beta) THEN
499 480 : WRITE (UNIT=line(61:66), FMT="(F6.2)") qeff
500 : ELSE
501 3030 : 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 3510 : WRITE (UNIT=line(77:78), FMT="(A2)") ADJUSTR(TRIM(element_symbol))
505 3510 : IF (my_charge_extended) THEN
506 330 : WRITE (UNIT=line(81:), FMT="(SP,F0.8)") qeff
507 : END IF
508 3588 : WRITE (UNIT=iunit, FMT="(A)") TRIM(line)
509 : END DO
510 78 : WRITE (UNIT=iunit, FMT="(A)") "END"
511 : CASE DEFAULT
512 27675 : CPABORT("Illegal dump type")
513 : END SELECT
514 :
515 27597 : CALL timestop(handle)
516 :
517 27597 : 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 10233 : 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 10233 : NULLIFY (logger)
540 10233 : NULLIFY (shell_kind)
541 :
542 10233 : logger => cp_get_default_logger()
543 : iw = cp_print_key_unit_nr(logger, subsys_section, &
544 10233 : "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
545 :
546 10233 : CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
547 10233 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
548 10233 : CALL uppercase(unit_str)
549 10233 : IF (iw > 0) THEN
550 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
551 2520 : "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 2520 : "Atom Kind Name", "X", "Y", "Z", "q(eff)", "Mass"
554 2520 : natom = SIZE(particle_set)
555 366742 : 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 364222 : shell=shell_kind)
562 364222 : IF (PRESENT(charges)) qeff = charges(iatom)
563 364222 : IF (ASSOCIATED(shell_kind)) THEN
564 : CALL get_shell(shell=shell_kind, &
565 : charge_core=qcore, &
566 3426 : charge_shell=qshell)
567 3426 : 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 1823630 : iatom, ikind, name, particle_set(iatom)%r(1:3)*conv, qeff, mass/massunit
571 : END DO
572 2520 : WRITE (iw, "(A)") ""
573 : END IF
574 :
575 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
576 10233 : "PRINT%ATOMIC_COORDINATES")
577 :
578 10233 : 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 20733 : 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 20733 : CALL timeset(routineN, handle)
606 :
607 20733 : NULLIFY (logger)
608 20733 : logger => cp_get_default_logger()
609 : iw = cp_print_key_unit_nr(logger, subsys_section, &
610 20733 : "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
611 :
612 20733 : CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
613 20733 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
614 20733 : CALL uppercase(unit_str)
615 20733 : IF (iw > 0) THEN
616 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
617 4536 : "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 4536 : "Atom Kind Element", "X", "Y", "Z", "Z(eff)", "Mass"
620 4536 : natom = SIZE(particle_set)
621 26508 : 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 21972 : z=z)
627 21972 : 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 92424 : iatom, ikind, element_symbol, z, particle_set(iatom)%r(1:3)*conv, zeff, mass/massunit
630 : END DO
631 4536 : WRITE (iw, "(A)") ""
632 : END IF
633 :
634 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
635 20733 : "PRINT%ATOMIC_COORDINATES")
636 :
637 20733 : CALL timestop(handle)
638 :
639 20733 : 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 11962 : 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 11962 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: distance_matrix
664 : REAL(KIND=dp), DIMENSION(3) :: rab
665 : TYPE(cp_logger_type), POINTER :: logger
666 :
667 11962 : CALL timeset(routineN, handle)
668 :
669 11962 : CPASSERT(ASSOCIATED(particle_set))
670 11962 : CPASSERT(ASSOCIATED(cell))
671 11962 : CPASSERT(ASSOCIATED(subsys_section))
672 :
673 11962 : NULLIFY (logger)
674 11962 : logger => cp_get_default_logger()
675 : iw = cp_print_key_unit_nr(logger, subsys_section, &
676 11962 : "PRINT%INTERATOMIC_DISTANCES", extension=".distLog")
677 :
678 11962 : CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%UNIT", c_val=unit_str)
679 11962 : 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 11962 : r_val=dab_min, explicit=explicit)
682 :
683 11962 : dab_abort = 0.0_dp
684 11962 : dab_warn = 0.0_dp
685 11962 : 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 11962 : IF (explicit .OR. (iw > 0) .OR. (natom <= 2000)) THEN
690 11920 : IF (dab_min > 0.0_dp) THEN
691 11916 : dab_warn = dab_min*conv
692 4 : ELSE IF (dab_min < 0.0_dp) THEN
693 0 : dab_abort = ABS(dab_min)*conv
694 : END IF
695 : END IF
696 :
697 11962 : IF ((iw > 0) .OR. (dab_abort > 0.0_dp) .OR. (dab_warn > 0.0_dp)) THEN
698 11916 : CALL get_cell(cell=cell, periodic=periodic)
699 11916 : IF (iw > 0) THEN
700 132 : ALLOCATE (distance_matrix(natom, natom))
701 33 : distance_matrix(:, :) = 0.0_dp
702 : END IF
703 344679 : DO iatom = 1, natom
704 123382406 : DO jatom = iatom + 1, natom
705 : rab(:) = pbc(particle_set(iatom)%r(:), &
706 123037727 : particle_set(jatom)%r(:), cell)
707 123037727 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))*conv
708 123037727 : IF (dab_abort > 0.0_dp) THEN
709 : ! Stop the run for interatomic distances smaller than the requested threshold
710 0 : 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 0 : TRIM(ADJUSTL(unit_str)))
718 : END IF
719 : END IF
720 123037727 : 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 916 : TRIM(ADJUSTL(unit_str)))
729 : END IF
730 123370490 : IF (iw > 0) THEN
731 35186 : distance_matrix(iatom, jatom) = dab
732 35186 : distance_matrix(jatom, iatom) = distance_matrix(iatom, jatom)
733 : END IF
734 : END DO
735 : END DO
736 11916 : IF (iw > 0) THEN
737 : ! Print the distance matrix
738 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
739 33 : "INTERATOMIC DISTANCES IN "//TRIM(unit_str)
740 33 : CALL write_particle_matrix(distance_matrix, particle_set, iw)
741 33 : IF (ALLOCATED(distance_matrix)) DEALLOCATE (distance_matrix)
742 : END IF
743 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
744 11916 : "PRINT%INTERATOMIC_DISTANCES")
745 : END IF
746 :
747 11962 : CALL timestop(handle)
748 :
749 11962 : 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 63 : 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 63 : INTEGER, DIMENSION(:), POINTER :: my_list
774 :
775 63 : my_el_per_part = 1
776 30 : IF (PRESENT(el_per_part)) my_el_per_part = el_per_part
777 63 : my_parts_per_line = 5
778 63 : IF (PRESENT(parts_per_line)) my_parts_per_line = MAX(parts_per_line, 1)
779 : WRITE (fmt_string1, FMT='(A,I0,A)') &
780 63 : "(/,T2,9X,", my_parts_per_line, "(4X,I6,4X))"
781 : WRITE (fmt_string2, FMT='(A,I0,A)') &
782 63 : "(T2,I5,1X,A2,1X,", my_parts_per_line, "(1X,ES13.6E2))"
783 63 : IF (PRESENT(Ilist)) THEN
784 30 : natom = SIZE(Ilist)
785 : ELSE
786 33 : natom = SIZE(particle_set)
787 : END IF
788 189 : ALLOCATE (my_list(natom))
789 63 : IF (PRESENT(Ilist)) THEN
790 180 : my_list = Ilist
791 : ELSE
792 927 : DO i = 1, natom
793 927 : my_list(i) = i
794 : END DO
795 : END IF
796 63 : natom = natom*my_el_per_part
797 317 : DO jatom = 1, natom, my_parts_per_line
798 254 : from = jatom
799 254 : to = MIN(from + my_parts_per_line - 1, natom)
800 1373 : WRITE (UNIT=iw, FMT=TRIM(fmt_string1)) (icol, icol=from, to)
801 15622 : DO iatom = 1, natom
802 15305 : katom = iatom/my_el_per_part
803 15305 : 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 15305 : element_symbol=element_symbol)
806 : WRITE (UNIT=iw, FMT=TRIM(fmt_string2)) &
807 15305 : iatom, element_symbol, &
808 30864 : (matrix(iatom, icol), icol=from, to)
809 : END DO
810 : END DO
811 :
812 63 : DEALLOCATE (my_list)
813 :
814 63 : 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 62126 : 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 62126 : INTEGER, ALLOCATABLE, DIMENSION(:) :: work
842 62126 : 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 62126 : CALL timeset(routineN, handle)
850 62126 : NULLIFY (atomic_indices)
851 62126 : NULLIFY (index_list)
852 62126 : NULLIFY (logger)
853 62126 : NULLIFY (section)
854 62126 : string = ""
855 :
856 62126 : 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 62126 : extension=".coordLog")
861 :
862 62126 : CALL section_vals_val_get(input_section, "PRINT%STRUCTURE_DATA%UNIT", c_val=unit_str)
863 62126 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
864 62126 : CALL uppercase(unit_str)
865 62126 : IF (iw > 0) THEN
866 568 : natom = SIZE(particle_set)
867 : section => section_vals_get_subs_vals(section_vals=input_section, &
868 568 : subsection_name="PRINT%STRUCTURE_DATA")
869 :
870 568 : 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 568 : n_rep_val=n_rep)
875 568 : IF (n_rep > 0) THEN
876 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
877 164 : "Position vectors r(i) of the atoms i in "//TRIM(unit_str)
878 164 : old_size = 0
879 939 : 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 775 : i_vals=atomic_indices)
884 775 : n_vals = SIZE(atomic_indices)
885 775 : new_size = old_size + n_vals
886 775 : CALL reallocate(index_list, 1, new_size)
887 3337 : index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
888 939 : old_size = new_size
889 : END DO
890 492 : ALLOCATE (work(new_size))
891 164 : CALL sort(index_list, new_size, work)
892 164 : DEALLOCATE (work)
893 1445 : DO i = 1, new_size
894 1281 : WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")"
895 1281 : IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
896 : WRITE (UNIT=iw, FMT="(T3,A)") &
897 33 : "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored."
898 33 : CYCLE
899 : END IF
900 1248 : IF (i > 1) THEN
901 : ! Skip redundant indices
902 1095 : IF (index_list(i) == index_list(i - 1)) CYCLE
903 : END IF
904 : WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") &
905 5156 : "r"//TRIM(string), "=", pbc(particle_set(index_list(i))%r(1:3), cell)*conv
906 : END DO
907 164 : 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 568 : n_rep_val=n_rep)
914 568 : IF (n_rep > 0) THEN
915 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
916 36 : "Position vectors s(i) of the atoms i in scaled coordinates"
917 36 : old_size = 0
918 105 : 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 69 : i_vals=atomic_indices)
923 69 : n_vals = SIZE(atomic_indices)
924 69 : new_size = old_size + n_vals
925 69 : CALL reallocate(index_list, 1, new_size)
926 1219 : index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
927 105 : old_size = new_size
928 : END DO
929 108 : ALLOCATE (work(new_size))
930 36 : CALL sort(index_list, new_size, work)
931 36 : DEALLOCATE (work)
932 611 : DO i = 1, new_size
933 575 : WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")"
934 575 : IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
935 : WRITE (UNIT=iw, FMT="(T3,A)") &
936 33 : "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored."
937 33 : CYCLE
938 : END IF
939 542 : IF (i > 1) THEN
940 : ! Skip redundant indices
941 517 : IF (index_list(i) == index_list(i - 1)) CYCLE
942 : END IF
943 542 : r(1:3) = pbc(particle_set(index_list(i))%r(1:3), cell)
944 542 : CALL real_to_scaled(s, r, cell)
945 : WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") &
946 578 : "s"//TRIM(string), "=", s(1:3)
947 : END DO
948 36 : 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 568 : n_rep_val=n)
955 568 : 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 129 : TRIM(unit_str)
959 355 : DO i = 1, n
960 : CALL section_vals_val_get(section_vals=section, &
961 : keyword_name="DISTANCE", &
962 : i_rep_val=i, &
963 226 : i_vals=atomic_indices)
964 226 : string = ""
965 : WRITE (UNIT=string, FMT="(A,2(I0,A))") &
966 226 : "(", atomic_indices(1), ",", atomic_indices(2), ")"
967 678 : wrk2 = atomic_indices
968 226 : CALL sort_unique(wrk2, unique)
969 355 : IF (((wrk2(1) >= 1) .AND. (wrk2(SIZE(wrk2)) <= natom)) .AND. unique) THEN
970 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
971 226 : particle_set(atomic_indices(2))%r(:), cell)
972 226 : 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 904 : "r"//TRIM(string), "=", rab(:)*conv, &
975 452 : "|r| =", dab*conv
976 : ELSE
977 : WRITE (UNIT=iw, FMT="(T3,A)") &
978 0 : "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 568 : n_rep_val=n)
987 568 : 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 67 : "r(j,k) in DEGREE"
991 139 : DO i = 1, n
992 : CALL section_vals_val_get(section_vals=section, &
993 : keyword_name="ANGLE", &
994 : i_rep_val=i, &
995 72 : i_vals=atomic_indices)
996 72 : string = ""
997 : WRITE (UNIT=string, FMT="(A,3(I0,A))") &
998 72 : "(", atomic_indices(1), ",", atomic_indices(2), ",", atomic_indices(3), ")"
999 288 : wrk3 = atomic_indices
1000 72 : CALL sort_unique(wrk3, unique)
1001 139 : IF (((wrk3(1) >= 1) .AND. (wrk3(SIZE(wrk3)) <= natom)) .AND. unique) THEN
1002 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1003 67 : particle_set(atomic_indices(2))%r(:), cell)
1004 : rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1005 67 : particle_set(atomic_indices(3))%r(:), cell)
1006 : WRITE (UNIT=iw, FMT="(T3,A,T26,A,F9.3)") &
1007 268 : "a"//TRIM(string), "=", angle(-rab, rbc)*degree
1008 : ELSE
1009 : WRITE (UNIT=iw, FMT="(T3,A)") &
1010 5 : "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 568 : n_rep_val=n)
1019 568 : 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 6 : "in DEGREE"
1023 17 : DO i = 1, n
1024 : CALL section_vals_val_get(section_vals=section, &
1025 : keyword_name="DIHEDRAL_ANGLE", &
1026 : i_rep_val=i, &
1027 11 : i_vals=atomic_indices)
1028 11 : string = ""
1029 : WRITE (UNIT=string, FMT="(A,4(I0,A))") &
1030 11 : "(", atomic_indices(1), ",", atomic_indices(2), ",", &
1031 22 : atomic_indices(3), ",", atomic_indices(4), ")"
1032 55 : wrk4 = atomic_indices
1033 11 : CALL sort_unique(wrk4, unique)
1034 17 : IF (((wrk4(1) >= 1) .AND. (wrk4(SIZE(wrk4)) <= natom)) .AND. unique) THEN
1035 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1036 1 : particle_set(atomic_indices(2))%r(:), cell)
1037 : rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1038 1 : particle_set(atomic_indices(3))%r(:), cell)
1039 : rcd(:) = pbc(particle_set(atomic_indices(3))%r(:), &
1040 1 : 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 1 : "d"//TRIM(string), "=", -dihedral_angle(rab, rbc, rcd)*degree
1044 : ELSE
1045 : WRITE (UNIT=iw, FMT="(T3,A)") &
1046 10 : "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 62126 : "PRINT%STRUCTURE_DATA")
1053 :
1054 62126 : CALL timestop(handle)
1055 :
1056 62126 : 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 1119 : 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, constraint_label
1094 :
1095 : CHARACTER(len=*), PARAMETER :: routineN = 'write_final_structure'
1096 :
1097 : CHARACTER(LEN=1) :: conv_str
1098 : CHARACTER(LEN=2) :: element_symbol
1099 1119 : CHARACTER(LEN=2), ALLOCATABLE :: element_list(:)
1100 : CHARACTER(LEN=5) :: pbc_str
1101 1119 : CHARACTER(LEN=:), ALLOCATABLE :: formula_structural, formula_sum
1102 : CHARACTER(LEN=default_path_length) :: cell_str, record, title
1103 : CHARACTER(LEN=default_string_length) :: atm_name, f_cif, f_cif_label, &
1104 : f_cif_type_symbol
1105 1119 : CHARACTER(LEN=default_string_length), ALLOCATABLE :: cif_label(:), cif_type_symbol(:)
1106 : CHARACTER(LEN=timestamp_length) :: timestamp
1107 : INTEGER :: elem_seen, file_unit, gcd_all, handle, i, iatom, ielem, natom, output_unit, &
1108 : symmetry_id, w_cif_label, w_cif_type_symbol
1109 1119 : INTEGER, ALLOCATABLE :: count_list(:)
1110 : LOGICAL :: dummy, elem_in_list, orthorhombic, &
1111 : print_final_structure, print_kind, &
1112 : write_cif, write_xyz
1113 : REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma, &
1114 : deth, unit_conv
1115 : REAL(KIND=dp), DIMENSION(3) :: abc, r, s
1116 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1117 : TYPE(cp_logger_type), POINTER :: logger
1118 : TYPE(enumeration_type), POINTER :: enum
1119 : TYPE(keyword_type), POINTER :: symmetry_keyword
1120 : TYPE(section_type), POINTER :: tmp_cell_section
1121 : TYPE(section_vals_type), POINTER :: print_key
1122 :
1123 1119 : CALL timeset(routineN, handle)
1124 :
1125 1119 : NULLIFY (enum, logger, symmetry_keyword, print_key, tmp_cell_section)
1126 1119 : logger => cp_get_default_logger()
1127 1119 : output_unit = cp_logger_get_default_io_unit(logger)
1128 1119 : print_key => section_vals_get_subs_vals(input_section, "PRINT%FINAL_STRUCTURE")
1129 1119 : print_final_structure = BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)
1130 1119 : file_unit = -1
1131 1119 : conv_str = "F"
1132 1119 : IF (conv) conv_str(1:1) = "T"
1133 :
1134 : ! Collect cell information
1135 1119 : pbc_str = "F F F"
1136 : CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, &
1137 : deth=deth, orthorhombic=orthorhombic, abc=abc, h=hmat, &
1138 1119 : symmetry_id=symmetry_id)
1139 1119 : IF (cell%perd(1) == 1) pbc_str(1:1) = "T"
1140 1119 : IF (cell%perd(2) == 1) pbc_str(3:3) = "T"
1141 1119 : IF (cell%perd(3) == 1) pbc_str(5:5) = "T"
1142 1119 : CALL create_cell_section(tmp_cell_section)
1143 1119 : symmetry_keyword => section_get_keyword(tmp_cell_section, "SYMMETRY")
1144 1119 : CALL keyword_get(symmetry_keyword, enum=enum)
1145 : ! cell_str is default_path_length which is longer
1146 : ! than default_string_length and should be enough
1147 : WRITE (UNIT=cell_str, FMT="(9(1X,F19.10))") &
1148 1119 : cp_unit_from_cp2k(hmat(1, 1), "angstrom"), &
1149 1119 : cp_unit_from_cp2k(hmat(2, 1), "angstrom"), &
1150 1119 : cp_unit_from_cp2k(hmat(3, 1), "angstrom"), &
1151 1119 : cp_unit_from_cp2k(hmat(1, 2), "angstrom"), &
1152 1119 : cp_unit_from_cp2k(hmat(2, 2), "angstrom"), &
1153 1119 : cp_unit_from_cp2k(hmat(3, 2), "angstrom"), &
1154 1119 : cp_unit_from_cp2k(hmat(1, 3), "angstrom"), &
1155 1119 : cp_unit_from_cp2k(hmat(2, 3), "angstrom"), &
1156 2238 : cp_unit_from_cp2k(hmat(3, 3), "angstrom")
1157 :
1158 : ! Collect atom information
1159 1119 : natom = SIZE(particle_set)
1160 1119 : ALLOCATE (element_list(nelem + 1), count_list(nelem + 1))
1161 1119 : count_list(:) = 0
1162 4476 : ALLOCATE (cif_label(natom), cif_type_symbol(natom))
1163 1119 : elem_seen = 0
1164 1119 : w_cif_type_symbol = 0
1165 1119 : w_cif_label = 0
1166 44755 : atom_loop: DO iatom = 1, natom
1167 43636 : elem_in_list = .FALSE.
1168 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1169 43636 : name=atm_name, element_symbol=element_symbol)
1170 43636 : cif_type_symbol(iatom) = TRIM(atm_name)
1171 : ! From write_particle_coordinates above it seems possible
1172 : ! for some atoms to have empty element symbols; whatever
1173 : ! these are, do not count them in the chemical formula
1174 43636 : IF (LEN_TRIM(element_symbol) == 0) THEN
1175 0 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
1176 0 : cif_label(iatom) = TRIM(atm_name)//TRIM(ADJUSTL(cp_to_string(iatom)))
1177 : ELSE
1178 43636 : cif_label(iatom) = TRIM(element_symbol)//TRIM(ADJUSTL(cp_to_string(iatom)))
1179 73392 : elem_loop: DO ielem = 1, elem_seen
1180 73392 : IF (element_list(ielem) == element_symbol) THEN
1181 41467 : elem_in_list = .TRUE.
1182 41467 : count_list(ielem) = count_list(ielem) + 1
1183 : EXIT elem_loop
1184 : END IF
1185 : END DO elem_loop
1186 : IF (.NOT. elem_in_list) THEN
1187 2169 : elem_seen = elem_seen + 1
1188 2169 : element_list(elem_seen) = element_symbol
1189 2169 : count_list(elem_seen) = 1
1190 : END IF
1191 : END IF
1192 43636 : IF (LEN_TRIM(cif_type_symbol(iatom)) > w_cif_type_symbol) THEN
1193 : w_cif_type_symbol = LEN_TRIM(cif_type_symbol(iatom))
1194 : END IF
1195 44755 : IF (LEN_TRIM(cif_label(iatom)) > w_cif_label) THEN
1196 : w_cif_label = LEN_TRIM(cif_label(iatom))
1197 : END IF
1198 : END DO atom_loop
1199 :
1200 : ! Determine the format of each line in cif considering width of cif_type_symbol and cif_label
1201 : ! The fields are, in order:
1202 : ! _atom_site_type_symbol, _atom_site_label, _atom_site_symmetry_multiplicity,
1203 : ! _atom_site_fract_x, _atom_site_fract_y, _atom_site_fract_z, _atom_site_occupancy
1204 : ! in which:
1205 : ! _atom_site_type_symbol is taken as atm_name
1206 : ! _atom_site_label is taken as element_symbol//iatom
1207 : ! _atom_site_symmetry_multiplicity and _atom_site_occupancy are always 1
1208 1119 : f_cif_type_symbol = "A"//TRIM(ADJUSTL(cp_to_string(w_cif_type_symbol + 4)))
1209 1119 : f_cif_label = "A"//TRIM(ADJUSTL(cp_to_string(w_cif_label + 4)))
1210 1119 : f_cif = "(T3,"//TRIM(f_cif_type_symbol)//","//TRIM(f_cif_label)//",I4,3F14.8,F8.2)"
1211 :
1212 : ! Determine formula_sum
1213 1119 : CPASSERT(elem_seen > 0)
1214 1119 : CPASSERT(count_list(1) > 0)
1215 1119 : formula_sum = "'"
1216 3288 : DO ielem = 1, elem_seen
1217 2169 : formula_sum = formula_sum//TRIM(ADJUSTL(element_list(ielem)))
1218 2169 : formula_sum = formula_sum//TRIM(ADJUSTL(cp_to_string(count_list(ielem))))
1219 3288 : formula_sum = formula_sum//" "
1220 : END DO
1221 1119 : formula_sum = TRIM(ADJUSTL(formula_sum))//"'"
1222 :
1223 : ! Determine formula_structural and Z
1224 1119 : gcd_all = count_list(1)
1225 3288 : DO ielem = 1, elem_seen
1226 3288 : IF (count_list(ielem) /= 0) THEN
1227 2169 : gcd_all = gcd(gcd_all, count_list(ielem))
1228 : END IF
1229 : END DO
1230 66688 : IF (gcd_all > 1) count_list = count_list/gcd_all
1231 1119 : formula_structural = "'"
1232 3288 : DO ielem = 1, elem_seen
1233 2169 : formula_structural = formula_structural//TRIM(ADJUSTL(element_list(ielem)))
1234 2169 : formula_structural = formula_structural//TRIM(ADJUSTL(cp_to_string(count_list(ielem))))
1235 3288 : formula_structural = formula_structural//" "
1236 : END DO
1237 1119 : formula_structural = TRIM(ADJUSTL(formula_structural))//"'"
1238 :
1239 : ! Write XYZ
1240 1119 : CALL section_vals_val_get(print_key, "PRINT_XYZ", l_val=write_xyz)
1241 1119 : write_xyz = write_xyz .AND. print_final_structure
1242 1119 : CALL section_vals_val_get(print_key, "PRINT_ATOM_KIND", l_val=print_kind)
1243 1119 : IF (write_xyz) THEN
1244 : ! Print a message to log
1245 : record = cp_print_key_generate_filename(logger, print_key, &
1246 : extension=".xyz", &
1247 1106 : my_local=.FALSE.)
1248 1106 : IF (output_unit > 0) THEN
1249 583 : IF (conv) THEN
1250 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1251 214 : routineN//": Optimization converged, writing XYZ file gladly:"
1252 : ELSE
1253 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1254 369 : routineN//": Optimization not yet converged, writing XYZ file anyway:"
1255 : END IF
1256 583 : WRITE (UNIT=output_unit, FMT="(T3,A)") TRIM(record)
1257 : END IF
1258 :
1259 : ! Prepare title
1260 : WRITE (UNIT=title, FMT="(A)") &
1261 : 'Lattice="'//TRIM(ADJUSTL(cell_str))//'" '// &
1262 : 'Properties=species:S:1:pos:R:3 '// &
1263 : 'pbc="'//pbc_str//'" '// &
1264 1106 : 'Converged='//conv_str
1265 : ! Extended XYZ uses angstrom for positions
1266 1106 : unit_conv = cp_unit_from_cp2k(1.0_dp, "angstrom")
1267 : ! Prepare file unit and write to it
1268 : file_unit = cp_print_key_unit_nr(logger, input_section, "PRINT%FINAL_STRUCTURE", &
1269 : file_status="REPLACE", file_form="FORMATTED", &
1270 1106 : extension=".xyz")
1271 1106 : IF (file_unit > 0) THEN
1272 : CALL write_particle_coordinates(particle_set, file_unit, dump_extxyz, "POS", title, &
1273 583 : cell=cell, unit_conv=unit_conv, print_kind=print_kind)
1274 : END IF
1275 : END IF
1276 :
1277 : ! Write CIF
1278 1119 : CALL section_vals_val_get(print_key, "PRINT_CIF", l_val=write_cif)
1279 1119 : write_cif = write_cif .AND. print_final_structure
1280 1119 : IF (write_cif) THEN
1281 : ! Print a message to log
1282 : record = cp_print_key_generate_filename(logger, print_key, &
1283 : extension=".cif", &
1284 1106 : my_local=.FALSE.)
1285 1106 : IF (output_unit > 0) THEN
1286 583 : IF (conv) THEN
1287 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1288 214 : routineN//": Optimization converged, writing CIF file gladly:"
1289 : ELSE
1290 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1291 369 : routineN//": Optimization not yet converged, writing CIF file anyway:"
1292 : END IF
1293 583 : WRITE (UNIT=output_unit, FMT="(T3,A)") TRIM(record)
1294 : END IF
1295 :
1296 : ! Make timestamp for the file
1297 1106 : CALL m_timestamp(timestamp)
1298 :
1299 : ! Prepare file unit and write to it
1300 : file_unit = cp_print_key_unit_nr(logger, input_section, "PRINT%FINAL_STRUCTURE", &
1301 : file_status="REPLACE", file_form="FORMATTED", &
1302 1106 : extension=".cif")
1303 1106 : IF (file_unit > 0) THEN
1304 : ! Generic information
1305 : WRITE (UNIT=file_unit, FMT="(A)") &
1306 583 : "# CIF file created by CP2K "//TRIM(moduleN)//":"//TRIM(routineN)
1307 : WRITE (UNIT=file_unit, FMT="(A)") &
1308 583 : "data_"//TRIM(logger%iter_info%project_name)
1309 : WRITE (UNIT=file_unit, FMT="(A,T39,A)") &
1310 583 : "_audit_creation_date", timestamp(:10)
1311 : WRITE (UNIT=file_unit, FMT="(A,/,A,/,A)") &
1312 583 : "_audit_creation_method", ";", &
1313 1166 : TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")"
1314 : WRITE (UNIT=file_unit, FMT="(A,/,A,/,A,/,A)") &
1315 583 : "Project name "//TRIM(logger%iter_info%project_name), &
1316 583 : "submitted by "//TRIM(r_user_name)//"@"//TRIM(r_host_name), &
1317 583 : "processed in "//TRIM(r_cwd), &
1318 1166 : "generated at "//TRIM(timestamp)
1319 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1320 583 : REPEAT(" -*WARNING*- ", 6)
1321 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1322 583 : "WARNING| the lines enclosed in between -*WARNING*- contain metadata"
1323 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1324 583 : "WARNING| provided in a tentative format for the newly implemented"
1325 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1326 583 : "WARNING| CIF output from CP2K version 2026.2, which may be unstable"
1327 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1328 583 : "WARNING| and subject to overhaul in a future release. Contact the"
1329 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1330 583 : "WARNING| developers in case a stabilized format is needed for some"
1331 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1332 583 : "WARNING| downstream CIF parser with text pattern matchers (regexp)."
1333 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1334 583 : "- Optimization type: "//TRIM(gopt_env_label)
1335 583 : IF (conv) THEN
1336 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1337 214 : "- Optimization converged: TRUE"
1338 : ELSE
1339 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1340 369 : "- Optimization converged: FALSE"
1341 : END IF
1342 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1343 583 : "- Requested initial cell symmetry: "//TRIM(enum_i2c(enum, symmetry_id))
1344 583 : IF (orthorhombic) THEN
1345 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1346 465 : "- Cell is numerically orthorhombic: TRUE"
1347 : ELSE
1348 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1349 118 : "- Cell is numerically orthorhombic: FALSE"
1350 : END IF
1351 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1352 583 : "- Periodicity of cell: "//TRIM(pbc_str)
1353 583 : IF (gopt_env_label == "CELL_OPT") THEN
1354 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1355 117 : "- Cell is subject to optimization: TRUE"
1356 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1357 117 : "- Cell has constraint on direction: "//TRIM(ADJUSTL(constraint_label))
1358 117 : IF (keep_angles) THEN
1359 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1360 21 : "- Keep angles between the cell vectors during optimization: TRUE"
1361 : ELSE
1362 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1363 96 : "- Keep angles between the cell vectors during optimization: FALSE"
1364 : END IF
1365 117 : IF (keep_symmetry) THEN
1366 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1367 21 : "- Keep initial cell symmetry during optimization: TRUE"
1368 : ELSE
1369 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1370 96 : "- Keep initial cell symmetry during optimization: FALSE"
1371 : END IF
1372 117 : IF (keep_volume) THEN
1373 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1374 3 : "- Keep initial cell volume during optimization: TRUE"
1375 : ELSE
1376 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1377 114 : "- Keep initial cell volume during optimization: FALSE"
1378 : END IF
1379 : ELSE
1380 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1381 466 : "- Cell is subject to optimization: FALSE"
1382 : END IF
1383 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1384 583 : "- Final cell vectors A, B, C by rows [angstrom]:"
1385 2332 : DO i = 1, 3
1386 : WRITE (UNIT=file_unit, FMT="(T3,3(1X,F19.10))") &
1387 1749 : cp_unit_from_cp2k(hmat(1, i), "angstrom"), &
1388 1749 : cp_unit_from_cp2k(hmat(2, i), "angstrom"), &
1389 4081 : cp_unit_from_cp2k(hmat(3, i), "angstrom")
1390 : END DO
1391 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1392 583 : REPEAT("-*WARNING*- ", 6)
1393 583 : WRITE (UNIT=file_unit, FMT="(A)") ";"
1394 : ! Data of cell and geometry
1395 : WRITE (UNIT=file_unit, FMT="(/,A,T44,A)") &
1396 583 : "_symmetry_space_group_name_H-M", "'P 1'"
1397 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1398 583 : "_cell_length_a", cp_unit_from_cp2k(abc(1), "angstrom")
1399 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1400 583 : "_cell_length_b", cp_unit_from_cp2k(abc(2), "angstrom")
1401 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1402 583 : "_cell_length_c", cp_unit_from_cp2k(abc(3), "angstrom")
1403 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1404 583 : "_cell_angle_alpha", angle_alpha
1405 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1406 583 : "_cell_angle_beta", angle_beta
1407 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1408 583 : "_cell_angle_gamma", angle_gamma
1409 : WRITE (UNIT=file_unit, FMT="(A,T48,A)") &
1410 583 : "_symmetry_Int_Tables_number", "1"
1411 : WRITE (UNIT=file_unit, FMT="(A,T36,A)") &
1412 583 : "_chemical_formula_structural", formula_structural
1413 : WRITE (UNIT=file_unit, FMT="(A,T36,A)") &
1414 583 : "_chemical_formula_sum", formula_sum
1415 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1416 583 : "_cell_volume", cp_unit_from_cp2k(ABS(deth), "angstrom^3")
1417 : WRITE (UNIT=file_unit, FMT="(A,T41,I8)") &
1418 583 : "_cell_formula_units_Z", gcd_all
1419 : WRITE (UNIT=file_unit, FMT="(A,/,T2,A,/,T2,A,/,T3,A)") &
1420 583 : "loop_", "_symmetry_equiv_pos_site_id", &
1421 1166 : "_symmetry_equiv_pos_as_xyz", "1 'x, y, z'"
1422 : WRITE (UNIT=file_unit, FMT="(A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A)") &
1423 583 : "loop_", "_atom_site_type_symbol", "_atom_site_label", &
1424 583 : "_atom_site_symmetry_multiplicity", "_atom_site_fract_x", &
1425 1166 : "_atom_site_fract_y", "_atom_site_fract_z", "_atom_site_occupancy"
1426 22846 : DO iatom = 1, natom
1427 : ! positive_range=.TRUE. makes r(1:3) in [0, cell%hmat(i,i)] and
1428 : ! s(1:3) in [0, 1], so there is no need to MODULO s(1:3) by 1.0
1429 22263 : r(1:3) = pbc(particle_set(iatom)%r(1:3), cell, positive_range=.TRUE.)
1430 22263 : CALL real_to_scaled(s, r, cell)
1431 : WRITE (UNIT=file_unit, FMT=TRIM(f_cif)) &
1432 22846 : cif_type_symbol(iatom), cif_label(iatom), 1, s(1:3), 1.0_dp
1433 : END DO
1434 : END IF
1435 : END IF
1436 :
1437 : ! Finish
1438 0 : DEALLOCATE (element_list, count_list, formula_structural, &
1439 1119 : formula_sum, cif_label, cif_type_symbol)
1440 1119 : CALL section_release(tmp_cell_section)
1441 : CALL cp_print_key_finished_output(file_unit, logger, input_section, &
1442 1119 : "PRINT%FINAL_STRUCTURE")
1443 1119 : IF (output_unit > 0 .AND. (write_xyz .OR. write_cif)) THEN
1444 : WRITE (UNIT=output_unit, FMT='(/,T2,A)') &
1445 583 : routineN//": Done!"
1446 : END IF
1447 :
1448 1119 : CALL timestop(handle)
1449 :
1450 5595 : END SUBROUTINE write_final_structure
1451 :
1452 15252 : END MODULE particle_methods
|
