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_print_key_finished_output,&
36 : cp_print_key_generate_filename,&
37 : cp_print_key_unit_nr
38 : USE cp_units, ONLY: cp_unit_from_cp2k
39 : USE external_potential_types, ONLY: fist_potential_type,&
40 : get_potential
41 : USE input_constants, ONLY: dump_atomic,&
42 : dump_dcd,&
43 : dump_dcd_aligned_cell,&
44 : dump_extxyz,&
45 : dump_pdb,&
46 : dump_xmol
47 : USE input_cp2k_subsys, ONLY: create_cell_section
48 : USE input_enumeration_types, ONLY: enum_i2c,&
49 : enumeration_type
50 : USE input_keyword_types, ONLY: keyword_get,&
51 : keyword_type
52 : USE input_section_types, ONLY: section_get_keyword,&
53 : section_release,&
54 : section_type,&
55 : section_vals_get_subs_vals,&
56 : section_vals_type,&
57 : section_vals_val_get
58 : USE kinds, ONLY: default_path_length,&
59 : default_string_length,&
60 : dp,&
61 : sp
62 : USE machine, ONLY: m_timestamp,&
63 : timestamp_length
64 : USE mathconstants, ONLY: degree
65 : USE mathlib, ONLY: angle,&
66 : dihedral_angle,&
67 : gcd
68 : USE memory_utilities, ONLY: reallocate
69 : USE particle_types, ONLY: get_particle_pos_or_vel,&
70 : particle_type
71 : USE periodic_table, ONLY: nelem
72 : USE physcon, ONLY: massunit
73 : USE qmmm_ff_fist, ONLY: qmmm_ff_precond_only_qm
74 : USE qs_kind_types, ONLY: get_qs_kind,&
75 : qs_kind_type
76 : USE shell_potential_types, ONLY: get_shell,&
77 : shell_kind_type
78 : USE string_utilities, ONLY: uppercase
79 : USE util, ONLY: sort,&
80 : sort_unique
81 : #include "./base/base_uses.f90"
82 :
83 : IMPLICIT NONE
84 :
85 : PRIVATE
86 :
87 : ! Public subroutines
88 :
89 : PUBLIC :: write_fist_particle_coordinates, &
90 : write_qs_particle_coordinates, &
91 : write_particle_distances, &
92 : write_particle_coordinates, &
93 : write_structure_data, &
94 : get_particle_set, &
95 : write_particle_matrix, &
96 : write_final_structure
97 :
98 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'particle_methods'
99 :
100 : CONTAINS
101 :
102 : ! **************************************************************************************************
103 : !> \brief Get the components of a particle set.
104 : !> \param particle_set ...
105 : !> \param qs_kind_set ...
106 : !> \param first_sgf ...
107 : !> \param last_sgf ...
108 : !> \param nsgf ...
109 : !> \param nmao ...
110 : !> \param basis ...
111 : !> \param ncgf ...
112 : !> \date 14.01.2002
113 : !> \par History
114 : !> - particle type cleaned (13.10.2003,MK)
115 : !> - refactoring and add basis set option (17.08.2010,jhu)
116 : !> \author MK
117 : !> \version 1.0
118 : ! **************************************************************************************************
119 195219 : SUBROUTINE get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, &
120 195219 : nmao, basis, ncgf)
121 :
122 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
123 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
124 : INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: first_sgf, last_sgf, nsgf, nmao
125 : TYPE(gto_basis_set_p_type), DIMENSION(:), OPTIONAL :: basis
126 : INTEGER, DIMENSION(:), INTENT(INOUT), OPTIONAL :: ncgf
127 :
128 : INTEGER :: ikind, iparticle, isgf, nparticle, ns
129 :
130 195219 : CPASSERT(ASSOCIATED(particle_set))
131 :
132 195219 : nparticle = SIZE(particle_set)
133 195219 : IF (PRESENT(first_sgf)) THEN
134 42595 : CPASSERT(SIZE(first_sgf) >= nparticle)
135 : END IF
136 195219 : IF (PRESENT(last_sgf)) THEN
137 33941 : CPASSERT(SIZE(last_sgf) >= nparticle)
138 : END IF
139 195219 : IF (PRESENT(nsgf)) THEN
140 152176 : CPASSERT(SIZE(nsgf) >= nparticle)
141 : END IF
142 195219 : IF (PRESENT(nmao)) THEN
143 14 : CPASSERT(SIZE(nmao) >= nparticle)
144 : END IF
145 195219 : IF (PRESENT(ncgf)) THEN
146 4 : CPASSERT(SIZE(ncgf) >= nparticle)
147 : END IF
148 :
149 195219 : IF (PRESENT(first_sgf) .OR. PRESENT(last_sgf) .OR. PRESENT(nsgf)) THEN
150 : isgf = 0
151 1106372 : DO iparticle = 1, nparticle
152 911649 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
153 911649 : IF (PRESENT(basis)) THEN
154 681983 : IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
155 681979 : CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, nsgf=ns)
156 : ELSE
157 4 : ns = 0
158 : END IF
159 : ELSE
160 229666 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns)
161 : END IF
162 911649 : IF (PRESENT(nsgf)) nsgf(iparticle) = ns
163 911649 : IF (PRESENT(first_sgf)) first_sgf(iparticle) = isgf + 1
164 911649 : isgf = isgf + ns
165 2018517 : IF (PRESENT(last_sgf)) last_sgf(iparticle) = isgf
166 : END DO
167 : END IF
168 :
169 195219 : IF (PRESENT(ncgf)) THEN
170 12 : DO iparticle = 1, nparticle
171 8 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
172 8 : IF (PRESENT(basis)) THEN
173 8 : IF (ASSOCIATED(basis(ikind)%gto_basis_set)) THEN
174 8 : CALL get_gto_basis_set(gto_basis_set=basis(ikind)%gto_basis_set, ncgf=ns)
175 : ELSE
176 0 : ns = 0
177 : END IF
178 : ELSE
179 0 : CALL get_qs_kind(qs_kind_set(ikind), ncgf=ns)
180 : END IF
181 20 : ncgf(iparticle) = ns
182 : END DO
183 : END IF
184 :
185 195219 : IF (PRESENT(first_sgf)) THEN
186 42595 : IF (SIZE(first_sgf) > nparticle) first_sgf(nparticle + 1) = isgf + 1
187 : END IF
188 :
189 195219 : IF (PRESENT(nmao)) THEN
190 86 : DO iparticle = 1, nparticle
191 72 : CALL get_atomic_kind(particle_set(iparticle)%atomic_kind, kind_number=ikind)
192 72 : CALL get_qs_kind(qs_kind_set(ikind), mao=ns)
193 86 : nmao(iparticle) = ns
194 : END DO
195 : END IF
196 :
197 195219 : END SUBROUTINE get_particle_set
198 :
199 : ! **************************************************************************************************
200 : !> \brief Should be able to write a few formats e.g. xmol, and some binary
201 : !> format (dcd) some format can be used for x,v,f
202 : !>
203 : !> FORMAT CONTENT UNITS x, v, f
204 : !> XMOL POS, VEL, FORCE, POS_VEL, POS_VEL_FORCE Angstrom, a.u., a.u.
205 : !>
206 : !> \param particle_set ...
207 : !> \param iunit ...
208 : !> \param output_format ...
209 : !> \param content ...
210 : !> \param title ...
211 : !> \param cell ...
212 : !> \param array ...
213 : !> \param unit_conv ...
214 : !> \param charge_occup ...
215 : !> \param charge_beta ...
216 : !> \param charge_extended ...
217 : !> \param print_kind ...
218 : !> \date 14.01.2002
219 : !> \author MK
220 : !> \version 1.0
221 : ! **************************************************************************************************
222 26565 : SUBROUTINE write_particle_coordinates(particle_set, iunit, output_format, &
223 26565 : content, title, cell, array, unit_conv, &
224 : charge_occup, charge_beta, &
225 : charge_extended, print_kind)
226 :
227 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
228 : INTEGER :: iunit, output_format
229 : CHARACTER(LEN=*) :: content, title
230 : TYPE(cell_type), OPTIONAL, POINTER :: cell
231 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: array
232 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: unit_conv
233 : LOGICAL, INTENT(IN), OPTIONAL :: charge_occup, charge_beta, &
234 : charge_extended, print_kind
235 :
236 : CHARACTER(len=*), PARAMETER :: routineN = 'write_particle_coordinates'
237 :
238 : CHARACTER(LEN=120) :: line
239 : CHARACTER(LEN=2) :: element_symbol
240 : CHARACTER(LEN=4) :: name
241 : CHARACTER(LEN=default_string_length) :: atm_name, my_format
242 : INTEGER :: handle, iatom, natom
243 : LOGICAL :: dummy, my_charge_beta, &
244 : my_charge_extended, my_charge_occup, &
245 : my_print_kind
246 : REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma, &
247 : factor, qeff
248 26565 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: arr
249 : REAL(KIND=dp), DIMENSION(3) :: abc, angles, f, r, v
250 : REAL(KIND=dp), DIMENSION(3, 3) :: h
251 26565 : REAL(KIND=sp), ALLOCATABLE, DIMENSION(:) :: x4, y4, z4
252 : TYPE(cell_type), POINTER :: cell_dcd
253 : TYPE(fist_potential_type), POINTER :: fist_potential
254 : TYPE(shell_kind_type), POINTER :: shell
255 :
256 26565 : CALL timeset(routineN, handle)
257 :
258 26565 : natom = SIZE(particle_set)
259 26565 : IF (PRESENT(array)) THEN
260 1741 : SELECT CASE (TRIM(content))
261 : CASE ("POS_VEL", "POS_VEL_FORCE")
262 1741 : CPABORT("Illegal usage")
263 : END SELECT
264 : END IF
265 26565 : factor = 1.0_dp
266 26565 : IF (PRESENT(unit_conv)) THEN
267 26414 : factor = unit_conv
268 : END IF
269 53057 : SELECT CASE (output_format)
270 : CASE (dump_xmol, dump_extxyz)
271 26492 : my_print_kind = .FALSE.
272 26492 : IF (PRESENT(print_kind)) my_print_kind = print_kind
273 26492 : WRITE (iunit, "(I8)") natom
274 26492 : WRITE (iunit, "(A)") TRIM(title)
275 1322347 : DO iatom = 1, natom
276 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
277 1295855 : element_symbol=element_symbol)
278 1295855 : IF (LEN_TRIM(element_symbol) == 0 .OR. my_print_kind) THEN
279 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
280 24 : name=atm_name)
281 24 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
282 24 : my_format = "(A,"
283 24 : atm_name = TRIM(atm_name)
284 : ELSE
285 1295831 : my_format = "(T2,A2,"
286 1295831 : atm_name = TRIM(element_symbol)
287 : END IF
288 26492 : SELECT CASE (TRIM(content))
289 : CASE ("POS")
290 1180519 : IF (PRESENT(array)) THEN
291 47741 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
292 : ELSE
293 4531112 : r(:) = particle_set(iatom)%r(:)
294 : END IF
295 4722076 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), r(1:3)*factor
296 : CASE ("VEL")
297 85772 : IF (PRESENT(array)) THEN
298 0 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
299 : ELSE
300 343088 : v(:) = particle_set(iatom)%v(:)
301 : END IF
302 343088 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), v(1:3)*factor
303 : CASE ("FORCE")
304 20955 : IF (PRESENT(array)) THEN
305 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
306 : ELSE
307 83820 : f(:) = particle_set(iatom)%f(:)
308 : END IF
309 83820 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), f(1:3)*factor
310 : CASE ("FORCE_MIXING_LABELS")
311 8609 : IF (PRESENT(array)) THEN
312 34436 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
313 : ELSE
314 0 : f(:) = particle_set(iatom)%f(:)
315 : END IF
316 1330291 : WRITE (iunit, TRIM(my_format)//"1X,3F20.10)") TRIM(atm_name), f(1:3)*factor
317 : END SELECT
318 : END DO
319 : CASE (dump_atomic)
320 170 : DO iatom = 1, natom
321 10 : SELECT CASE (TRIM(content))
322 : CASE ("POS")
323 160 : IF (PRESENT(array)) THEN
324 0 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
325 : ELSE
326 640 : r(:) = particle_set(iatom)%r(:)
327 : END IF
328 640 : WRITE (iunit, "(3F20.10)") r(1:3)*factor
329 : CASE ("VEL")
330 0 : IF (PRESENT(array)) THEN
331 0 : v(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
332 : ELSE
333 0 : v(:) = particle_set(iatom)%v(:)
334 : END IF
335 0 : WRITE (iunit, "(3F20.10)") v(1:3)*factor
336 : CASE ("FORCE")
337 0 : IF (PRESENT(array)) THEN
338 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
339 : ELSE
340 0 : f(:) = particle_set(iatom)%f(:)
341 : END IF
342 0 : WRITE (iunit, "(3F20.10)") f(1:3)*factor
343 : CASE ("FORCE_MIXING_LABELS")
344 0 : IF (PRESENT(array)) THEN
345 0 : f(:) = array((iatom - 1)*3 + 1:(iatom - 1)*3 + 3)
346 : ELSE
347 0 : f(:) = particle_set(iatom)%f(:)
348 : END IF
349 160 : WRITE (iunit, "(3F20.10)") f(1:3)*factor
350 : END SELECT
351 : END DO
352 : CASE (dump_dcd, dump_dcd_aligned_cell)
353 4 : IF (.NOT. (PRESENT(cell))) THEN
354 0 : CPABORT("Cell is not present! Report this bug!")
355 : END IF
356 : CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, &
357 4 : abc=abc)
358 4 : IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN
359 : ! In the case of a non-orthorhombic cell adopt a common convention
360 : ! for the orientation of the cell with respect to the Cartesian axes:
361 : ! Cell vector a is aligned with the x axis and the cell vector b lies
362 : ! in the xy plane.
363 0 : NULLIFY (cell_dcd)
364 0 : CALL cell_create(cell_dcd)
365 0 : CALL cell_clone(cell, cell_dcd, tag="CELL_DCD")
366 0 : angles(1) = angle_alpha/degree
367 0 : angles(2) = angle_beta/degree
368 0 : angles(3) = angle_gamma/degree
369 : CALL set_cell_param(cell_dcd, abc, angles, &
370 0 : do_init_cell=.TRUE.)
371 0 : h(1:3, 1:3) = MATMUL(cell_dcd%hmat(1:3, 1:3), cell%h_inv(1:3, 1:3))
372 0 : CALL cell_release(cell_dcd)
373 : END IF
374 12 : ALLOCATE (arr(3, natom))
375 4 : IF (PRESENT(array)) THEN
376 0 : arr(1:3, 1:natom) = RESHAPE(array, [3, natom])
377 : ELSE
378 8 : SELECT CASE (TRIM(content))
379 : CASE ("POS")
380 1156 : DO iatom = 1, natom
381 4612 : arr(1:3, iatom) = particle_set(iatom)%r(1:3)
382 : END DO
383 : CASE ("VEL")
384 0 : DO iatom = 1, natom
385 0 : arr(1:3, iatom) = particle_set(iatom)%v(1:3)
386 : END DO
387 : CASE ("FORCE")
388 0 : DO iatom = 1, natom
389 0 : arr(1:3, iatom) = particle_set(iatom)%f(1:3)
390 : END DO
391 : CASE DEFAULT
392 4 : CPABORT("Illegal DCD dump type")
393 : END SELECT
394 : END IF
395 12 : ALLOCATE (x4(natom))
396 8 : ALLOCATE (y4(natom))
397 8 : ALLOCATE (z4(natom))
398 4 : IF (.NOT. cell%orthorhombic .AND. (output_format == dump_dcd_aligned_cell)) THEN
399 0 : x4(1:natom) = REAL(MATMUL(h(1, 1:3), arr(1:3, 1:natom)), KIND=sp)
400 0 : y4(1:natom) = REAL(MATMUL(h(2, 1:3), arr(1:3, 1:natom)), KIND=sp)
401 0 : z4(1:natom) = REAL(MATMUL(h(3, 1:3), arr(1:3, 1:natom)), KIND=sp)
402 : ELSE
403 1156 : x4(1:natom) = REAL(arr(1, 1:natom), KIND=sp)
404 1156 : y4(1:natom) = REAL(arr(2, 1:natom), KIND=sp)
405 1156 : z4(1:natom) = REAL(arr(3, 1:natom), KIND=sp)
406 : END IF
407 4 : WRITE (iunit) abc(1)*factor, angle_gamma, abc(2)*factor, &
408 8 : angle_beta, angle_alpha, abc(3)*factor
409 1156 : WRITE (iunit) x4*REAL(factor, KIND=sp)
410 1156 : WRITE (iunit) y4*REAL(factor, KIND=sp)
411 1156 : WRITE (iunit) z4*REAL(factor, KIND=sp)
412 : ! Release work storage
413 4 : DEALLOCATE (arr)
414 4 : DEALLOCATE (x4)
415 4 : DEALLOCATE (y4)
416 8 : DEALLOCATE (z4)
417 : CASE (dump_pdb)
418 59 : my_charge_occup = .FALSE.
419 59 : IF (PRESENT(charge_occup)) my_charge_occup = charge_occup
420 59 : my_charge_beta = .FALSE.
421 59 : IF (PRESENT(charge_beta)) my_charge_beta = charge_beta
422 59 : my_charge_extended = .FALSE.
423 59 : IF (PRESENT(charge_extended)) my_charge_extended = charge_extended
424 59 : IF (LEN_TRIM(title) > 0) THEN
425 : WRITE (UNIT=iunit, FMT="(A6,T11,A)") &
426 59 : "REMARK", TRIM(title)
427 : END IF
428 59 : CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, abc=abc)
429 : ! COLUMNS DATA TYPE CONTENTS
430 : ! --------------------------------------------------
431 : ! 1 - 6 Record name "CRYST1"
432 : ! 7 - 15 Real(9.3) a (Angstroms)
433 : ! 16 - 24 Real(9.3) b (Angstroms)
434 : ! 25 - 33 Real(9.3) c (Angstroms)
435 : ! 34 - 40 Real(7.2) alpha (degrees)
436 : ! 41 - 47 Real(7.2) beta (degrees)
437 : ! 48 - 54 Real(7.2) gamma (degrees)
438 : ! 56 - 66 LString Space group
439 : ! 67 - 70 Integer Z value
440 : WRITE (UNIT=iunit, FMT="(A6,3F9.3,3F7.2)") &
441 236 : "CRYST1", abc(1:3)*factor, angle_alpha, angle_beta, angle_gamma
442 59 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
443 2999 : DO iatom = 1, natom
444 2940 : line = ""
445 : ! COLUMNS DATA TYPE CONTENTS
446 : ! 1 - 6 Record name "ATOM "
447 : ! 7 - 11 Integer Atom serial number
448 : ! 13 - 16 Atom Atom name
449 : ! 17 Character Alternate location indicator
450 : ! 18 - 20 Residue name Residue name
451 : ! 22 Character Chain identifier
452 : ! 23 - 26 Integer Residue sequence number
453 : ! 27 AChar Code for insertion of residues
454 : ! 31 - 38 Real(8.3) Orthogonal coordinates for X in Angstrom
455 : ! 39 - 46 Real(8.3) Orthogonal coordinates for Y in Angstrom
456 : ! 47 - 54 Real(8.3) Orthogonal coordinates for Z in Angstrom
457 : ! 55 - 60 Real(6.2) Occupancy
458 : ! 61 - 66 Real(6.2) Temperature factor (Default = 0.0)
459 : ! 73 - 76 LString(4) Segment identifier, left-justified
460 : ! 77 - 78 LString(2) Element symbol, right-justified
461 : ! 79 - 80 LString(2) Charge on the atom
462 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
463 : element_symbol=element_symbol, name=atm_name, &
464 2940 : fist_potential=fist_potential, shell=shell)
465 2940 : IF (LEN_TRIM(element_symbol) == 0) THEN
466 0 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
467 : END IF
468 2940 : name = TRIM(atm_name)
469 2940 : IF (ASSOCIATED(fist_potential)) THEN
470 2940 : CALL get_potential(potential=fist_potential, qeff=qeff)
471 : ELSE
472 0 : qeff = 0.0_dp
473 : END IF
474 2940 : IF (ASSOCIATED(shell)) CALL get_shell(shell=shell, charge=qeff)
475 2940 : WRITE (UNIT=line(1:6), FMT="(A6)") "ATOM "
476 2940 : WRITE (UNIT=line(7:11), FMT="(I5)") MODULO(iatom, 100000)
477 2940 : WRITE (UNIT=line(13:16), FMT="(A4)") ADJUSTL(name)
478 : ! WRITE (UNIT=line(18:20),FMT="(A3)") TRIM(resname)
479 : ! WRITE (UNIT=line(23:26),FMT="(I4)") MODULO(idres,10000)
480 5880 : SELECT CASE (TRIM(content))
481 : CASE ("POS")
482 2940 : IF (PRESENT(array)) THEN
483 0 : r(1:3) = get_particle_pos_or_vel(iatom, particle_set, array)
484 : ELSE
485 11760 : r(:) = particle_set(iatom)%r(:)
486 : END IF
487 11760 : WRITE (UNIT=line(31:54), FMT="(3F8.3)") r(1:3)*factor
488 : CASE DEFAULT
489 2940 : CPABORT("PDB dump only for trajectory available")
490 : END SELECT
491 2940 : IF (my_charge_occup) THEN
492 2130 : WRITE (UNIT=line(55:60), FMT="(F6.2)") qeff
493 : ELSE
494 810 : WRITE (UNIT=line(55:60), FMT="(F6.2)") 0.0_dp
495 : END IF
496 2940 : IF (my_charge_beta) THEN
497 480 : WRITE (UNIT=line(61:66), FMT="(F6.2)") qeff
498 : ELSE
499 2460 : WRITE (UNIT=line(61:66), FMT="(F6.2)") 0.0_dp
500 : END IF
501 : ! WRITE (UNIT=line(73:76),FMT="(A4)") ADJUSTL(TRIM(molname))
502 2940 : WRITE (UNIT=line(77:78), FMT="(A2)") ADJUSTR(TRIM(element_symbol))
503 2940 : IF (my_charge_extended) THEN
504 330 : WRITE (UNIT=line(81:), FMT="(SP,F0.8)") qeff
505 : END IF
506 2999 : WRITE (UNIT=iunit, FMT="(A)") TRIM(line)
507 : END DO
508 59 : WRITE (UNIT=iunit, FMT="(A)") "END"
509 : CASE DEFAULT
510 26624 : CPABORT("Illegal dump type")
511 : END SELECT
512 :
513 26565 : CALL timestop(handle)
514 :
515 26565 : END SUBROUTINE write_particle_coordinates
516 :
517 : ! **************************************************************************************************
518 : !> \brief Write the atomic coordinates to the output unit.
519 : !> \param particle_set ...
520 : !> \param subsys_section ...
521 : !> \param charges ...
522 : !> \date 05.06.2000
523 : !> \author MK
524 : !> \version 1.0
525 : ! **************************************************************************************************
526 9689 : SUBROUTINE write_fist_particle_coordinates(particle_set, subsys_section, charges)
527 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
528 : TYPE(section_vals_type), POINTER :: subsys_section
529 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: charges
530 :
531 : CHARACTER(LEN=default_string_length) :: name, unit_str
532 : INTEGER :: iatom, ikind, iw, natom
533 : REAL(KIND=dp) :: conv, mass, qcore, qeff, qshell
534 : TYPE(cp_logger_type), POINTER :: logger
535 : TYPE(shell_kind_type), POINTER :: shell_kind
536 :
537 9689 : NULLIFY (logger)
538 9689 : NULLIFY (shell_kind)
539 :
540 9689 : logger => cp_get_default_logger()
541 : iw = cp_print_key_unit_nr(logger, subsys_section, &
542 9689 : "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
543 :
544 9689 : CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
545 9689 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
546 9689 : CALL uppercase(unit_str)
547 9689 : IF (iw > 0) THEN
548 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
549 2504 : "MODULE FIST: ATOMIC COORDINATES IN "//TRIM(unit_str)
550 : WRITE (UNIT=iw, FMT="(/,T4,A,T30,A,T44,A,T58,A,T66,A,T77,A)") &
551 2504 : "Atom Kind Name", "X", "Y", "Z", "q(eff)", "Mass"
552 2504 : natom = SIZE(particle_set)
553 362909 : DO iatom = 1, natom
554 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
555 : kind_number=ikind, &
556 : name=name, &
557 : mass=mass, &
558 : qeff=qeff, &
559 360405 : shell=shell_kind)
560 360405 : IF (PRESENT(charges)) qeff = charges(iatom)
561 360405 : IF (ASSOCIATED(shell_kind)) THEN
562 : CALL get_shell(shell=shell_kind, &
563 : charge_core=qcore, &
564 3426 : charge_shell=qshell)
565 3426 : qeff = qcore + qshell
566 : END IF
567 : WRITE (UNIT=iw, FMT="(T2,I6,1X,I4,1X,A7,3(1X,F13.6),2(1X,F8.4))") &
568 1804529 : iatom, ikind, name, particle_set(iatom)%r(1:3)*conv, qeff, mass/massunit
569 : END DO
570 2504 : WRITE (iw, "(A)") ""
571 : END IF
572 :
573 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
574 9689 : "PRINT%ATOMIC_COORDINATES")
575 :
576 9689 : END SUBROUTINE write_fist_particle_coordinates
577 :
578 : ! **************************************************************************************************
579 : !> \brief Write the atomic coordinates to the output unit.
580 : !> \param particle_set ...
581 : !> \param qs_kind_set ...
582 : !> \param subsys_section ...
583 : !> \param label ...
584 : !> \date 05.06.2000
585 : !> \author MK
586 : !> \version 1.0
587 : ! **************************************************************************************************
588 19138 : SUBROUTINE write_qs_particle_coordinates(particle_set, qs_kind_set, subsys_section, label)
589 :
590 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
591 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
592 : TYPE(section_vals_type), POINTER :: subsys_section
593 : CHARACTER(LEN=*), INTENT(IN) :: label
594 :
595 : CHARACTER(len=*), PARAMETER :: routineN = 'write_qs_particle_coordinates'
596 :
597 : CHARACTER(LEN=2) :: element_symbol
598 : CHARACTER(LEN=default_string_length) :: unit_str
599 : INTEGER :: handle, iatom, ikind, iw, natom, z
600 : REAL(KIND=dp) :: conv, mass, zeff
601 : TYPE(cp_logger_type), POINTER :: logger
602 :
603 19138 : CALL timeset(routineN, handle)
604 :
605 19138 : NULLIFY (logger)
606 19138 : logger => cp_get_default_logger()
607 : iw = cp_print_key_unit_nr(logger, subsys_section, &
608 19138 : "PRINT%ATOMIC_COORDINATES", extension=".coordLog")
609 :
610 19138 : CALL section_vals_val_get(subsys_section, "PRINT%ATOMIC_COORDINATES%UNIT", c_val=unit_str)
611 19138 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
612 19138 : CALL uppercase(unit_str)
613 19138 : IF (iw > 0) THEN
614 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
615 4439 : "MODULE "//TRIM(label)//": ATOMIC COORDINATES IN "//TRIM(unit_str)
616 : WRITE (UNIT=iw, FMT="(/,T4,A,T30,A,T44,A,T58,A,T66,A,T77,A)") &
617 4439 : "Atom Kind Element", "X", "Y", "Z", "Z(eff)", "Mass"
618 4439 : natom = SIZE(particle_set)
619 26233 : DO iatom = 1, natom
620 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
621 : kind_number=ikind, &
622 : element_symbol=element_symbol, &
623 : mass=mass, &
624 21794 : z=z)
625 21794 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff)
626 : WRITE (UNIT=iw, FMT="(T2,I6,1X,I4,1X,A2,1X,I4,3(1X,F13.6),2(1X,F8.4))") &
627 91615 : iatom, ikind, element_symbol, z, particle_set(iatom)%r(1:3)*conv, zeff, mass/massunit
628 : END DO
629 4439 : WRITE (iw, "(A)") ""
630 : END IF
631 :
632 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
633 19138 : "PRINT%ATOMIC_COORDINATES")
634 :
635 19138 : CALL timestop(handle)
636 :
637 19138 : END SUBROUTINE write_qs_particle_coordinates
638 :
639 : ! **************************************************************************************************
640 : !> \brief Write the matrix of the particle distances to the output unit.
641 : !> \param particle_set ...
642 : !> \param cell ...
643 : !> \param subsys_section ...
644 : !> \date 06.10.2000
645 : !> \author Matthias Krack
646 : !> \version 1.0
647 : ! **************************************************************************************************
648 11101 : SUBROUTINE write_particle_distances(particle_set, cell, subsys_section)
649 :
650 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
651 : TYPE(cell_type), POINTER :: cell
652 : TYPE(section_vals_type), POINTER :: subsys_section
653 :
654 : CHARACTER(len=*), PARAMETER :: routineN = 'write_particle_distances'
655 :
656 : CHARACTER(LEN=default_string_length) :: unit_str
657 : INTEGER :: handle, iatom, iw, jatom, natom
658 : INTEGER, DIMENSION(3) :: periodic
659 : LOGICAL :: explicit
660 : REAL(KIND=dp) :: conv, dab, dab_abort, dab_min, dab_warn
661 11101 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: distance_matrix
662 : REAL(KIND=dp), DIMENSION(3) :: rab
663 : TYPE(cp_logger_type), POINTER :: logger
664 :
665 11101 : CALL timeset(routineN, handle)
666 :
667 11101 : CPASSERT(ASSOCIATED(particle_set))
668 11101 : CPASSERT(ASSOCIATED(cell))
669 11101 : CPASSERT(ASSOCIATED(subsys_section))
670 :
671 11101 : NULLIFY (logger)
672 11101 : logger => cp_get_default_logger()
673 : iw = cp_print_key_unit_nr(logger, subsys_section, &
674 11101 : "PRINT%INTERATOMIC_DISTANCES", extension=".distLog")
675 :
676 11101 : CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%UNIT", c_val=unit_str)
677 11101 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
678 : CALL section_vals_val_get(subsys_section, "PRINT%INTERATOMIC_DISTANCES%CHECK_INTERATOMIC_DISTANCES", &
679 11101 : r_val=dab_min, explicit=explicit)
680 :
681 11101 : dab_abort = 0.0_dp
682 11101 : dab_warn = 0.0_dp
683 11101 : natom = SIZE(particle_set)
684 :
685 : ! Compute interatomic distances only if their printout or check is explicitly requested
686 : ! Disable the default check for systems with more than 3000 atoms
687 11101 : IF (explicit .OR. (iw > 0) .OR. (natom <= 2000)) THEN
688 11059 : IF (dab_min > 0.0_dp) THEN
689 11055 : dab_warn = dab_min*conv
690 4 : ELSE IF (dab_min < 0.0_dp) THEN
691 0 : dab_abort = ABS(dab_min)*conv
692 : END IF
693 : END IF
694 :
695 11101 : IF ((iw > 0) .OR. (dab_abort > 0.0_dp) .OR. (dab_warn > 0.0_dp)) THEN
696 11055 : CALL get_cell(cell=cell, periodic=periodic)
697 11055 : IF (iw > 0) THEN
698 132 : ALLOCATE (distance_matrix(natom, natom))
699 33 : distance_matrix(:, :) = 0.0_dp
700 : END IF
701 331960 : DO iatom = 1, natom
702 119868354 : DO jatom = iatom + 1, natom
703 : rab(:) = pbc(particle_set(iatom)%r(:), &
704 119536394 : particle_set(jatom)%r(:), cell)
705 119536394 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))*conv
706 119536394 : IF (dab_abort > 0.0_dp) THEN
707 : ! Stop the run for interatomic distances smaller than the requested threshold
708 0 : IF (dab < dab_abort) THEN
709 : CALL cp_abort(__LOCATION__, "The distance between the atoms "// &
710 : TRIM(ADJUSTL(cp_to_string(iatom, fmt="(I8)")))//" and "// &
711 : TRIM(ADJUSTL(cp_to_string(jatom, fmt="(I8)")))//" is only "// &
712 : TRIM(ADJUSTL(cp_to_string(dab, fmt="(F6.3)")))//" "// &
713 : TRIM(ADJUSTL(unit_str))//" and thus smaller than the requested threshold of "// &
714 : TRIM(ADJUSTL(cp_to_string(dab_abort, fmt="(F6.3)")))//" "// &
715 0 : TRIM(ADJUSTL(unit_str)))
716 : END IF
717 : END IF
718 119536394 : IF (dab < dab_warn) THEN
719 : ! Print warning for interatomic distances smaller than the requested threshold
720 : CALL cp_warn(__LOCATION__, "The distance between the atoms "// &
721 : TRIM(ADJUSTL(cp_to_string(iatom, fmt="(I8)")))//" and "// &
722 : TRIM(ADJUSTL(cp_to_string(jatom, fmt="(I8)")))//" is only "// &
723 : TRIM(ADJUSTL(cp_to_string(dab, fmt="(F6.3)")))//" "// &
724 : TRIM(ADJUSTL(unit_str))//" and thus smaller than the threshold of "// &
725 : TRIM(ADJUSTL(cp_to_string(dab_warn, fmt="(F6.3)")))//" "// &
726 910 : TRIM(ADJUSTL(unit_str)))
727 : END IF
728 119857299 : IF (iw > 0) THEN
729 35186 : distance_matrix(iatom, jatom) = dab
730 35186 : distance_matrix(jatom, iatom) = distance_matrix(iatom, jatom)
731 : END IF
732 : END DO
733 : END DO
734 11055 : IF (iw > 0) THEN
735 : ! Print the distance matrix
736 : WRITE (UNIT=iw, FMT="(/,/,T2,A)") &
737 33 : "INTERATOMIC DISTANCES IN "//TRIM(unit_str)
738 33 : CALL write_particle_matrix(distance_matrix, particle_set, iw)
739 33 : IF (ALLOCATED(distance_matrix)) DEALLOCATE (distance_matrix)
740 : END IF
741 : CALL cp_print_key_finished_output(iw, logger, subsys_section, &
742 11055 : "PRINT%INTERATOMIC_DISTANCES")
743 : END IF
744 :
745 11101 : CALL timestop(handle)
746 :
747 11101 : END SUBROUTINE write_particle_distances
748 :
749 : ! **************************************************************************************************
750 : !> \brief ...
751 : !> \param matrix ...
752 : !> \param particle_set ...
753 : !> \param iw ...
754 : !> \param el_per_part ...
755 : !> \param Ilist ...
756 : !> \param parts_per_line : number of particle columns to be printed in one line
757 : ! **************************************************************************************************
758 53 : SUBROUTINE write_particle_matrix(matrix, particle_set, iw, el_per_part, Ilist, parts_per_line)
759 : REAL(KIND=dp), DIMENSION(:, :) :: matrix
760 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
761 : INTEGER, INTENT(IN) :: iw
762 : INTEGER, INTENT(IN), OPTIONAL :: el_per_part
763 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: Ilist
764 : INTEGER, INTENT(IN), OPTIONAL :: parts_per_line
765 :
766 : CHARACTER(LEN=2) :: element_symbol
767 : CHARACTER(LEN=default_string_length) :: fmt_string1, fmt_string2
768 : INTEGER :: from, i, iatom, icol, jatom, katom, &
769 : my_el_per_part, my_parts_per_line, &
770 : natom, to
771 53 : INTEGER, DIMENSION(:), POINTER :: my_list
772 :
773 53 : my_el_per_part = 1
774 20 : IF (PRESENT(el_per_part)) my_el_per_part = el_per_part
775 53 : my_parts_per_line = 5
776 53 : IF (PRESENT(parts_per_line)) my_parts_per_line = MAX(parts_per_line, 1)
777 : WRITE (fmt_string1, FMT='(A,I0,A)') &
778 53 : "(/,T2,11X,", my_parts_per_line, "(4X,I5,4X))"
779 : WRITE (fmt_string2, FMT='(A,I0,A)') &
780 53 : "(T2,I5,2X,A2,2X,", my_parts_per_line, "(1X,F12.6))"
781 53 : IF (PRESENT(Ilist)) THEN
782 20 : natom = SIZE(Ilist)
783 : ELSE
784 33 : natom = SIZE(particle_set)
785 : END IF
786 159 : ALLOCATE (my_list(natom))
787 53 : IF (PRESENT(Ilist)) THEN
788 120 : my_list = Ilist
789 : ELSE
790 927 : DO i = 1, natom
791 927 : my_list(i) = i
792 : END DO
793 : END IF
794 53 : natom = natom*my_el_per_part
795 288 : DO jatom = 1, natom, my_parts_per_line
796 235 : from = jatom
797 235 : to = MIN(from + my_parts_per_line - 1, natom)
798 1279 : WRITE (UNIT=iw, FMT=TRIM(fmt_string1)) (icol, icol=from, to)
799 15446 : DO iatom = 1, natom
800 15158 : katom = iatom/my_el_per_part
801 15158 : IF (MOD(iatom, my_el_per_part) /= 0) katom = katom + 1
802 : CALL get_atomic_kind(atomic_kind=particle_set(my_list(katom))%atomic_kind, &
803 15158 : element_symbol=element_symbol)
804 : WRITE (UNIT=iw, FMT=TRIM(fmt_string2)) &
805 15158 : iatom, element_symbol, &
806 30551 : (matrix(iatom, icol), icol=from, to)
807 : END DO
808 : END DO
809 :
810 53 : DEALLOCATE (my_list)
811 :
812 53 : END SUBROUTINE write_particle_matrix
813 :
814 : ! **************************************************************************************************
815 : !> \brief Write structure data requested by a separate structure data input
816 : !> section to the output unit.
817 : !> input_section can be either motion_section or subsys_section.
818 : !>
819 : !> \param particle_set ...
820 : !> \param cell ...
821 : !> \param input_section ...
822 : !> \date 11.03.04
823 : !> \par History
824 : !> Recovered (23.03.06,MK)
825 : !> \author MK
826 : !> \version 1.0
827 : ! **************************************************************************************************
828 62193 : SUBROUTINE write_structure_data(particle_set, cell, input_section)
829 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
830 : TYPE(cell_type), POINTER :: cell
831 : TYPE(section_vals_type), POINTER :: input_section
832 :
833 : CHARACTER(LEN=*), PARAMETER :: routineN = 'write_structure_data'
834 :
835 : CHARACTER(LEN=default_string_length) :: string, unit_str
836 : INTEGER :: handle, i, i_rep, iw, n, n_rep, n_vals, &
837 : natom, new_size, old_size, wrk2(2), &
838 : wrk3(3), wrk4(4)
839 62193 : INTEGER, ALLOCATABLE, DIMENSION(:) :: work
840 62193 : INTEGER, DIMENSION(:), POINTER :: atomic_indices, index_list
841 : LOGICAL :: unique
842 : REAL(KIND=dp) :: conv, dab
843 : REAL(KIND=dp), DIMENSION(3) :: r, rab, rbc, rcd, s
844 : TYPE(cp_logger_type), POINTER :: logger
845 : TYPE(section_vals_type), POINTER :: section
846 :
847 62193 : CALL timeset(routineN, handle)
848 62193 : NULLIFY (atomic_indices)
849 62193 : NULLIFY (index_list)
850 62193 : NULLIFY (logger)
851 62193 : NULLIFY (section)
852 62193 : string = ""
853 :
854 62193 : logger => cp_get_default_logger()
855 : iw = cp_print_key_unit_nr(logger=logger, &
856 : basis_section=input_section, &
857 : print_key_path="PRINT%STRUCTURE_DATA", &
858 62193 : extension=".coordLog")
859 :
860 62193 : CALL section_vals_val_get(input_section, "PRINT%STRUCTURE_DATA%UNIT", c_val=unit_str)
861 62193 : conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
862 62193 : CALL uppercase(unit_str)
863 62193 : IF (iw > 0) THEN
864 568 : natom = SIZE(particle_set)
865 : section => section_vals_get_subs_vals(section_vals=input_section, &
866 568 : subsection_name="PRINT%STRUCTURE_DATA")
867 :
868 568 : WRITE (UNIT=iw, FMT="(/,T2,A)") "REQUESTED STRUCTURE DATA"
869 : ! Print the requested atomic position vectors
870 : CALL section_vals_val_get(section_vals=section, &
871 : keyword_name="POSITION", &
872 568 : n_rep_val=n_rep)
873 568 : IF (n_rep > 0) THEN
874 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
875 145 : "Position vectors r(i) of the atoms i in "//TRIM(unit_str)
876 145 : old_size = 0
877 848 : DO i_rep = 1, n_rep
878 : CALL section_vals_val_get(section_vals=section, &
879 : keyword_name="POSITION", &
880 : i_rep_val=i_rep, &
881 703 : i_vals=atomic_indices)
882 703 : n_vals = SIZE(atomic_indices)
883 703 : new_size = old_size + n_vals
884 703 : CALL reallocate(index_list, 1, new_size)
885 2903 : index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
886 848 : old_size = new_size
887 : END DO
888 435 : ALLOCATE (work(new_size))
889 145 : CALL sort(index_list, new_size, work)
890 145 : DEALLOCATE (work)
891 1245 : DO i = 1, new_size
892 1100 : WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")"
893 1100 : IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
894 : WRITE (UNIT=iw, FMT="(T3,A)") &
895 30 : "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored."
896 30 : CYCLE
897 : END IF
898 1070 : IF (i > 1) THEN
899 : ! Skip redundant indices
900 935 : IF (index_list(i) == index_list(i - 1)) CYCLE
901 : END IF
902 : WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") &
903 4425 : "r"//TRIM(string), "=", pbc(particle_set(index_list(i))%r(1:3), cell)*conv
904 : END DO
905 145 : DEALLOCATE (index_list)
906 : END IF
907 :
908 : ! Print the requested atomic position vectors in scaled coordinates
909 : CALL section_vals_val_get(section_vals=section, &
910 : keyword_name="POSITION_SCALED", &
911 568 : n_rep_val=n_rep)
912 568 : IF (n_rep > 0) THEN
913 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
914 27 : "Position vectors s(i) of the atoms i in scaled coordinates"
915 27 : old_size = 0
916 84 : DO i_rep = 1, n_rep
917 : CALL section_vals_val_get(section_vals=section, &
918 : keyword_name="POSITION_SCALED", &
919 : i_rep_val=i_rep, &
920 57 : i_vals=atomic_indices)
921 57 : n_vals = SIZE(atomic_indices)
922 57 : new_size = old_size + n_vals
923 57 : CALL reallocate(index_list, 1, new_size)
924 965 : index_list(old_size + 1:new_size) = atomic_indices(1:n_vals)
925 84 : old_size = new_size
926 : END DO
927 81 : ALLOCATE (work(new_size))
928 27 : CALL sort(index_list, new_size, work)
929 27 : DEALLOCATE (work)
930 481 : DO i = 1, new_size
931 454 : WRITE (UNIT=string, FMT="(A,I0,A)") "(", index_list(i), ")"
932 454 : IF ((index_list(i) < 1) .OR. (index_list(i) > natom)) THEN
933 : WRITE (UNIT=iw, FMT="(T3,A)") &
934 30 : "Invalid atomic index "//TRIM(string)//" specified. Print request is ignored."
935 30 : CYCLE
936 : END IF
937 424 : IF (i > 1) THEN
938 : ! Skip redundant indices
939 407 : IF (index_list(i) == index_list(i - 1)) CYCLE
940 : END IF
941 424 : r(1:3) = pbc(particle_set(index_list(i))%r(1:3), cell)
942 424 : CALL real_to_scaled(s, r, cell)
943 : WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6)") &
944 451 : "s"//TRIM(string), "=", s(1:3)
945 : END DO
946 27 : DEALLOCATE (index_list)
947 : END IF
948 :
949 : ! Print the requested distances
950 : CALL section_vals_val_get(section_vals=section, &
951 : keyword_name="DISTANCE", &
952 568 : n_rep_val=n)
953 568 : IF (n > 0) THEN
954 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
955 : "Distance vector r(i,j) between the atom i and j in "// &
956 129 : TRIM(unit_str)
957 355 : DO i = 1, n
958 : CALL section_vals_val_get(section_vals=section, &
959 : keyword_name="DISTANCE", &
960 : i_rep_val=i, &
961 226 : i_vals=atomic_indices)
962 226 : string = ""
963 : WRITE (UNIT=string, FMT="(A,2(I0,A))") &
964 226 : "(", atomic_indices(1), ",", atomic_indices(2), ")"
965 678 : wrk2 = atomic_indices
966 226 : CALL sort_unique(wrk2, unique)
967 355 : IF (((wrk2(1) >= 1) .AND. (wrk2(SIZE(wrk2)) <= natom)) .AND. unique) THEN
968 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
969 226 : particle_set(atomic_indices(2))%r(:), cell)
970 226 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
971 : WRITE (UNIT=iw, FMT="(T3,A,T20,A,3F13.6,3X,A,F13.6)") &
972 904 : "r"//TRIM(string), "=", rab(:)*conv, &
973 452 : "|r| =", dab*conv
974 : ELSE
975 : WRITE (UNIT=iw, FMT="(T3,A)") &
976 0 : "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored."
977 : END IF
978 : END DO
979 : END IF
980 :
981 : ! Print the requested angles
982 : CALL section_vals_val_get(section_vals=section, &
983 : keyword_name="ANGLE", &
984 568 : n_rep_val=n)
985 568 : IF (n > 0) THEN
986 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
987 : "Angle a(i,j,k) between the atomic distance vectors r(j,i) and "// &
988 67 : "r(j,k) in DEGREE"
989 139 : DO i = 1, n
990 : CALL section_vals_val_get(section_vals=section, &
991 : keyword_name="ANGLE", &
992 : i_rep_val=i, &
993 72 : i_vals=atomic_indices)
994 72 : string = ""
995 : WRITE (UNIT=string, FMT="(A,3(I0,A))") &
996 72 : "(", atomic_indices(1), ",", atomic_indices(2), ",", atomic_indices(3), ")"
997 288 : wrk3 = atomic_indices
998 72 : CALL sort_unique(wrk3, unique)
999 139 : IF (((wrk3(1) >= 1) .AND. (wrk3(SIZE(wrk3)) <= natom)) .AND. unique) THEN
1000 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1001 67 : particle_set(atomic_indices(2))%r(:), cell)
1002 : rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1003 67 : particle_set(atomic_indices(3))%r(:), cell)
1004 : WRITE (UNIT=iw, FMT="(T3,A,T26,A,F9.3)") &
1005 268 : "a"//TRIM(string), "=", angle(-rab, rbc)*degree
1006 : ELSE
1007 : WRITE (UNIT=iw, FMT="(T3,A)") &
1008 5 : "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored."
1009 : END IF
1010 : END DO
1011 : END IF
1012 :
1013 : ! Print the requested dihedral angles
1014 : CALL section_vals_val_get(section_vals=section, &
1015 : keyword_name="DIHEDRAL_ANGLE", &
1016 568 : n_rep_val=n)
1017 568 : IF (n > 0) THEN
1018 : WRITE (UNIT=iw, FMT="(/,T3,A,/)") &
1019 : "Dihedral angle d(i,j,k,l) between the planes (i,j,k) and (j,k,l) "// &
1020 5 : "in DEGREE"
1021 15 : DO i = 1, n
1022 : CALL section_vals_val_get(section_vals=section, &
1023 : keyword_name="DIHEDRAL_ANGLE", &
1024 : i_rep_val=i, &
1025 10 : i_vals=atomic_indices)
1026 10 : string = ""
1027 : WRITE (UNIT=string, FMT="(A,4(I0,A))") &
1028 10 : "(", atomic_indices(1), ",", atomic_indices(2), ",", &
1029 20 : atomic_indices(3), ",", atomic_indices(4), ")"
1030 50 : wrk4 = atomic_indices
1031 10 : CALL sort_unique(wrk4, unique)
1032 15 : IF (((wrk4(1) >= 1) .AND. (wrk4(SIZE(wrk4)) <= natom)) .AND. unique) THEN
1033 : rab(:) = pbc(particle_set(atomic_indices(1))%r(:), &
1034 0 : particle_set(atomic_indices(2))%r(:), cell)
1035 : rbc(:) = pbc(particle_set(atomic_indices(2))%r(:), &
1036 0 : particle_set(atomic_indices(3))%r(:), cell)
1037 : rcd(:) = pbc(particle_set(atomic_indices(3))%r(:), &
1038 0 : particle_set(atomic_indices(4))%r(:), cell)
1039 : WRITE (UNIT=iw, FMT="(T3,A,T26,A,F9.3)") &
1040 0 : "d"//TRIM(string), "=", dihedral_angle(rab, rbc, rcd)*degree
1041 : ELSE
1042 : WRITE (UNIT=iw, FMT="(T3,A)") &
1043 10 : "Invalid atomic indices "//TRIM(string)//" specified. Print request is ignored."
1044 : END IF
1045 : END DO
1046 : END IF
1047 : END IF
1048 : CALL cp_print_key_finished_output(iw, logger, input_section, &
1049 62193 : "PRINT%STRUCTURE_DATA")
1050 :
1051 62193 : CALL timestop(handle)
1052 :
1053 62193 : END SUBROUTINE write_structure_data
1054 :
1055 : ! **************************************************************************************************
1056 : !> \brief Write the final geometry and cell information to files
1057 : !> \param particle_set pointer to particles with atm_name, element_symbol and position
1058 : !> \param cell pointer to cell with abc, angle_alpha, angle_beta, angle_gamma and deth
1059 : !> \param input_section pointer to motion_section which has PRINT%FINAL_STRUCTURE
1060 : !> \param conv flag for whether convergence is achieved or not in optimization
1061 : !> \param keep_angles flag for whether cell optimization keeps initial angles
1062 : !> \param keep_symmetry flag for whether cell optimization keeps initial symmetry
1063 : !> \param keep_volume flag for whether cell optimization keeps initial volume
1064 : !> \param gopt_env_label the geometry optimization label "GEO_OPT", "CELL_OPT", ...
1065 : !> \param constraint_label label for directions with constraint in cell optimization
1066 : !> \par Intended to be invoked in gopt_f_methods:write_final_info.
1067 : !> This implementation does not consider higher space groups even if
1068 : !> one is detected, and the chemical formulae are neither written in
1069 : !> the sorted "Hill notation" nor expressed in groups of molecules.
1070 : !> Other potentially useful but yet to be written information includes:
1071 : !> the external pressure from CELL_OPT/EXTERNAL_POTENTIAL and the
1072 : !> stress tensor (virial) for CELL_OPT;
1073 : !> the fixed atoms from MOTION/CONSTRAINT/FIXED_ATOMS for all.
1074 : !>
1075 : !> History
1076 : !> 04.2026 - Created
1077 : !> \author HE Zilong
1078 : !> \version 1.0
1079 : ! **************************************************************************************************
1080 1073 : SUBROUTINE write_final_structure(particle_set, cell, input_section, conv, &
1081 : keep_angles, keep_symmetry, keep_volume, &
1082 : gopt_env_label, constraint_label)
1083 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1084 : TYPE(cell_type), INTENT(IN), POINTER :: cell
1085 : TYPE(section_vals_type), INTENT(IN), POINTER :: input_section
1086 : LOGICAL, INTENT(IN) :: conv, keep_angles, keep_symmetry, &
1087 : keep_volume
1088 : CHARACTER(LEN=default_string_length), INTENT(IN) :: gopt_env_label
1089 : CHARACTER(LEN=4), INTENT(IN) :: constraint_label
1090 :
1091 : CHARACTER(len=*), PARAMETER :: routineN = 'write_final_structure'
1092 :
1093 : CHARACTER(LEN=2) :: element_symbol
1094 1073 : CHARACTER(LEN=2), ALLOCATABLE :: element_list(:)
1095 : CHARACTER(LEN=5) :: perd_str
1096 1073 : CHARACTER(LEN=:), ALLOCATABLE :: formula_structural, formula_sum
1097 : CHARACTER(LEN=default_path_length) :: cell_str, record
1098 : CHARACTER(LEN=default_string_length) :: atm_name, f_cif, f_cif_label, &
1099 : f_cif_type_symbol, f_xyz
1100 1073 : CHARACTER(LEN=default_string_length), ALLOCATABLE :: cif_label(:), cif_type_symbol(:), &
1101 1073 : xyz_label(:)
1102 : CHARACTER(LEN=timestamp_length) :: timestamp
1103 : INTEGER :: elem_seen, file_unit, gcd_all, handle, i, iatom, ielem, natom, output_unit, &
1104 : symmetry_id, w_cif_label, w_cif_type_symbol, w_xyz_label
1105 1073 : INTEGER, ALLOCATABLE :: count_list(:)
1106 : INTEGER, DIMENSION(3) :: periodic
1107 : LOGICAL :: dummy, elem_in_list, orthorhombic, &
1108 : write_cif, write_xyz
1109 : REAL(KIND=dp) :: angle_alpha, angle_beta, angle_gamma, &
1110 : deth
1111 : REAL(KIND=dp), DIMENSION(3) :: abc, r, s
1112 : REAL(KIND=dp), DIMENSION(3, 3) :: hmat
1113 : TYPE(cp_logger_type), POINTER :: logger
1114 : TYPE(enumeration_type), POINTER :: enum
1115 : TYPE(keyword_type), POINTER :: symmetry_keyword
1116 : TYPE(section_type), POINTER :: tmp_cell_section
1117 : TYPE(section_vals_type), POINTER :: print_key
1118 :
1119 1073 : CALL timeset(routineN, handle)
1120 :
1121 1073 : NULLIFY (enum, logger, symmetry_keyword, print_key, tmp_cell_section)
1122 1073 : logger => cp_get_default_logger()
1123 1073 : output_unit = cp_logger_get_default_io_unit(logger)
1124 1073 : print_key => section_vals_get_subs_vals(input_section, "PRINT%FINAL_STRUCTURE")
1125 :
1126 : ! Collect cell information
1127 1073 : perd_str = "F F F"
1128 : CALL get_cell(cell, alpha=angle_alpha, beta=angle_beta, gamma=angle_gamma, &
1129 : deth=deth, orthorhombic=orthorhombic, abc=abc, periodic=periodic, &
1130 1073 : h=hmat, symmetry_id=symmetry_id)
1131 1073 : IF (periodic(1) == 1) perd_str(1:1) = "T"
1132 1073 : IF (periodic(2) == 1) perd_str(3:3) = "T"
1133 1073 : IF (periodic(3) == 1) perd_str(5:5) = "T"
1134 1073 : CALL create_cell_section(tmp_cell_section)
1135 1073 : symmetry_keyword => section_get_keyword(tmp_cell_section, "SYMMETRY")
1136 1073 : CALL keyword_get(symmetry_keyword, enum=enum)
1137 : ! cell_str is default_path_length which is longer
1138 : ! than default_string_length and should be enough
1139 : WRITE (UNIT=cell_str, FMT="(9(1X,F19.10))") &
1140 1073 : cp_unit_from_cp2k(hmat(1, 1), "angstrom"), &
1141 1073 : cp_unit_from_cp2k(hmat(2, 1), "angstrom"), &
1142 1073 : cp_unit_from_cp2k(hmat(3, 1), "angstrom"), &
1143 1073 : cp_unit_from_cp2k(hmat(1, 2), "angstrom"), &
1144 1073 : cp_unit_from_cp2k(hmat(2, 2), "angstrom"), &
1145 1073 : cp_unit_from_cp2k(hmat(3, 2), "angstrom"), &
1146 1073 : cp_unit_from_cp2k(hmat(1, 3), "angstrom"), &
1147 1073 : cp_unit_from_cp2k(hmat(2, 3), "angstrom"), &
1148 2146 : cp_unit_from_cp2k(hmat(3, 3), "angstrom")
1149 :
1150 : ! Collect atom information
1151 1073 : natom = SIZE(particle_set)
1152 1073 : ALLOCATE (element_list(nelem + 1), count_list(nelem + 1))
1153 1073 : count_list(:) = 0
1154 5365 : ALLOCATE (cif_label(natom), cif_type_symbol(natom), xyz_label(natom))
1155 1073 : elem_seen = 0
1156 1073 : w_cif_type_symbol = 0
1157 1073 : w_cif_label = 0
1158 1073 : w_xyz_label = 0
1159 37287 : atom_loop: DO iatom = 1, natom
1160 36214 : elem_in_list = .FALSE.
1161 : CALL get_atomic_kind(atomic_kind=particle_set(iatom)%atomic_kind, &
1162 36214 : name=atm_name, element_symbol=element_symbol)
1163 36214 : cif_type_symbol(iatom) = TRIM(atm_name)
1164 : ! From write_particle_coordinates above it seems possible
1165 : ! for some atoms to have empty element symbols; whatever
1166 : ! these are, do not count them in the chemical formula
1167 36214 : IF (LEN_TRIM(element_symbol) == 0) THEN
1168 0 : dummy = qmmm_ff_precond_only_qm(id1=atm_name)
1169 0 : cif_label(iatom) = TRIM(atm_name)//TRIM(ADJUSTL(cp_to_string(iatom)))
1170 0 : xyz_label(iatom) = TRIM(atm_name)
1171 : ELSE
1172 36214 : cif_label(iatom) = TRIM(element_symbol)//TRIM(ADJUSTL(cp_to_string(iatom)))
1173 36214 : xyz_label(iatom) = TRIM(element_symbol)
1174 60810 : elem_loop: DO ielem = 1, elem_seen
1175 60810 : IF (element_list(ielem) == element_symbol) THEN
1176 34139 : elem_in_list = .TRUE.
1177 34139 : count_list(ielem) = count_list(ielem) + 1
1178 : EXIT elem_loop
1179 : END IF
1180 : END DO elem_loop
1181 : IF (.NOT. elem_in_list) THEN
1182 2075 : elem_seen = elem_seen + 1
1183 2075 : element_list(elem_seen) = element_symbol
1184 2075 : count_list(elem_seen) = 1
1185 : END IF
1186 : END IF
1187 36214 : IF (LEN_TRIM(cif_type_symbol(iatom)) > w_cif_type_symbol) &
1188 : w_cif_type_symbol = LEN_TRIM(cif_type_symbol(iatom))
1189 36214 : IF (LEN_TRIM(cif_label(iatom)) > w_cif_label) &
1190 : w_cif_label = LEN_TRIM(cif_label(iatom))
1191 36214 : IF (LEN_TRIM(xyz_label(iatom)) > w_xyz_label) &
1192 1073 : w_xyz_label = LEN_TRIM(xyz_label(iatom))
1193 : END DO atom_loop
1194 :
1195 : ! Determine the format of each line in cif considering width of cif_type_symbol and cif_label
1196 : ! The fields are, in order:
1197 : ! _atom_site_type_symbol, _atom_site_label, _atom_site_symmetry_multiplicity,
1198 : ! _atom_site_fract_x, _atom_site_fract_y, _atom_site_fract_z, _atom_site_occupancy
1199 : ! in which:
1200 : ! _atom_site_type_symbol is taken as atm_name
1201 : ! _atom_site_label is taken as element_symbol//iatom
1202 : ! _atom_site_symmetry_multiplicity and _atom_site_occupancy are always 1
1203 1073 : f_cif_type_symbol = "A"//TRIM(ADJUSTL(cp_to_string(w_cif_type_symbol + 4)))
1204 1073 : f_cif_label = "A"//TRIM(ADJUSTL(cp_to_string(w_cif_label + 4)))
1205 1073 : f_cif = "(T3,"//TRIM(f_cif_type_symbol)//","//TRIM(f_cif_label)//",I4,3F14.8,F8.2)"
1206 :
1207 : ! Determine the format of each line in xyz
1208 1073 : f_xyz = "(T2,A"//TRIM(ADJUSTL(cp_to_string(w_xyz_label + 4)))//",1X,3F20.10)"
1209 :
1210 : ! Determine formula_sum
1211 1073 : CPASSERT(elem_seen > 0)
1212 1073 : CPASSERT(count_list(1) > 0)
1213 1073 : formula_sum = "'"
1214 3148 : DO ielem = 1, elem_seen
1215 2075 : formula_sum = formula_sum//TRIM(ADJUSTL(element_list(ielem)))
1216 2075 : formula_sum = formula_sum//TRIM(ADJUSTL(cp_to_string(count_list(ielem))))
1217 3148 : formula_sum = formula_sum//" "
1218 : END DO
1219 1073 : formula_sum = TRIM(ADJUSTL(formula_sum))//"'"
1220 :
1221 : ! Determine formula_structural and Z
1222 1073 : gcd_all = count_list(1)
1223 3148 : DO ielem = 1, elem_seen
1224 3148 : IF (count_list(ielem) /= 0) THEN
1225 2075 : gcd_all = gcd(gcd_all, count_list(ielem))
1226 : END IF
1227 : END DO
1228 63786 : IF (gcd_all > 1) count_list = count_list/gcd_all
1229 1073 : formula_structural = "'"
1230 3148 : DO ielem = 1, elem_seen
1231 2075 : formula_structural = formula_structural//TRIM(ADJUSTL(element_list(ielem)))
1232 2075 : formula_structural = formula_structural//TRIM(ADJUSTL(cp_to_string(count_list(ielem))))
1233 3148 : formula_structural = formula_structural//" "
1234 : END DO
1235 1073 : formula_structural = TRIM(ADJUSTL(formula_structural))//"'"
1236 :
1237 : ! Write XYZ
1238 1073 : CALL section_vals_val_get(print_key, "PRINT_XYZ", l_val=write_xyz)
1239 1073 : IF (write_xyz) THEN
1240 : ! Print a message to log
1241 : record = cp_print_key_generate_filename(logger, print_key, &
1242 : extension=".xyz", &
1243 1073 : my_local=.FALSE.)
1244 1073 : IF (output_unit > 0) THEN
1245 666 : IF (conv) THEN
1246 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1247 221 : routineN//": Optimization converged, writing XYZ file gladly:"
1248 : ELSE
1249 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1250 445 : routineN//": Optimization not yet converged, writing XYZ file anyway:"
1251 : END IF
1252 666 : WRITE (UNIT=output_unit, FMT="(T3,A)") TRIM(record)
1253 : END IF
1254 :
1255 : ! Make timestamp for the file
1256 1073 : CALL m_timestamp(timestamp)
1257 :
1258 : ! Prepare file unit and write to it
1259 : file_unit = cp_print_key_unit_nr(logger, input_section, "PRINT%FINAL_STRUCTURE", &
1260 1073 : file_status="REPLACE", extension=".xyz")
1261 1073 : IF (file_unit > 0) THEN
1262 : WRITE (UNIT=file_unit, FMT="(I8)") &
1263 579 : natom
1264 : WRITE (UNIT=file_unit, FMT="(A)") &
1265 : 'Lattice="'//TRIM(ADJUSTL(cell_str))// &
1266 : '" Properties=species:S:1:pos:R:3 pbc="'// &
1267 579 : TRIM(perd_str)//'"'
1268 19321 : DO iatom = 1, natom
1269 74968 : DO i = 1, 3
1270 74968 : r(i) = cp_unit_from_cp2k(particle_set(iatom)%r(i), "angstrom")
1271 : END DO
1272 : WRITE (UNIT=file_unit, FMT=TRIM(f_xyz)) &
1273 19321 : xyz_label(iatom), r(1:3)
1274 : END DO
1275 : END IF
1276 : END IF
1277 :
1278 : ! Write CIF
1279 1073 : CALL section_vals_val_get(print_key, "PRINT_CIF", l_val=write_cif)
1280 1073 : IF (write_cif) THEN
1281 : ! Print a message to log
1282 : record = cp_print_key_generate_filename(logger, print_key, &
1283 : extension=".cif", &
1284 1073 : my_local=.FALSE.)
1285 1073 : IF (output_unit > 0) THEN
1286 666 : IF (conv) THEN
1287 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1288 221 : routineN//": Optimization converged, writing CIF file gladly:"
1289 : ELSE
1290 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
1291 445 : routineN//": Optimization not yet converged, writing CIF file anyway:"
1292 : END IF
1293 666 : WRITE (UNIT=output_unit, FMT="(T3,A)") TRIM(record)
1294 : END IF
1295 :
1296 : ! Make timestamp for the file
1297 1073 : 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 1073 : file_status="REPLACE", extension=".cif")
1302 1073 : IF (file_unit > 0) THEN
1303 : ! Generic information
1304 : WRITE (UNIT=file_unit, FMT="(A)") &
1305 579 : "# CIF file created by CP2K "//TRIM(moduleN)//":"//TRIM(routineN)
1306 : WRITE (UNIT=file_unit, FMT="(A)") &
1307 579 : "data_"//TRIM(logger%iter_info%project_name)
1308 : WRITE (UNIT=file_unit, FMT="(A,T39,A)") &
1309 579 : "_audit_creation_date", timestamp(:10)
1310 : WRITE (UNIT=file_unit, FMT="(A,/,A,/,A)") &
1311 579 : "_audit_creation_method", ";", &
1312 1158 : TRIM(cp2k_version)//" (revision "//TRIM(compile_revision)//")"
1313 : WRITE (UNIT=file_unit, FMT="(A,/,A,/,A,/,A)") &
1314 579 : "Project name "//TRIM(logger%iter_info%project_name), &
1315 579 : "submitted by "//TRIM(r_user_name)//"@"//TRIM(r_host_name), &
1316 579 : "processed in "//TRIM(r_cwd), &
1317 1158 : "generated at "//TRIM(timestamp)
1318 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1319 579 : "- Optimization type: "//TRIM(gopt_env_label)
1320 579 : IF (conv) THEN
1321 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1322 221 : "- Optimization converged: TRUE"
1323 : ELSE
1324 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1325 358 : "- Optimization converged: FALSE"
1326 : END IF
1327 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1328 579 : "- Requested initial cell symmetry: "//TRIM(enum_i2c(enum, symmetry_id))
1329 579 : IF (orthorhombic) THEN
1330 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1331 468 : "- Cell is numerically orthorhombic: TRUE"
1332 : ELSE
1333 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1334 111 : "- Cell is numerically orthorhombic: FALSE"
1335 : END IF
1336 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1337 579 : "- Periodicity of cell: "//TRIM(perd_str)
1338 579 : IF (gopt_env_label == "CELL_OPT") THEN
1339 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1340 101 : "- Cell is subject to optimization: TRUE"
1341 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1342 101 : "- Cell has constraint on direction: "//TRIM(ADJUSTL(constraint_label))
1343 101 : IF (keep_angles) THEN
1344 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1345 15 : "- Keep angles between the cell vectors during optimization: TRUE"
1346 : ELSE
1347 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1348 86 : "- Keep angles between the cell vectors during optimization: FALSE"
1349 : END IF
1350 101 : IF (keep_symmetry) THEN
1351 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1352 21 : "- Keep initial cell symmetry during optimization: TRUE"
1353 : ELSE
1354 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1355 80 : "- Keep initial cell symmetry during optimization: FALSE"
1356 : END IF
1357 101 : IF (keep_volume) THEN
1358 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1359 3 : "- Keep initial cell volume during optimization: TRUE"
1360 : ELSE
1361 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1362 98 : "- Keep initial cell volume during optimization: FALSE"
1363 : END IF
1364 : ELSE
1365 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1366 478 : "- Cell is subject to optimization: FALSE"
1367 : END IF
1368 : WRITE (UNIT=file_unit, FMT="(T2,A)") &
1369 579 : "- Final cell vectors A, B, C by rows [angstrom]:"
1370 2316 : DO i = 1, 3
1371 : WRITE (UNIT=file_unit, FMT="(T3,3(1X,F19.10))") &
1372 1737 : cp_unit_from_cp2k(hmat(1, i), "angstrom"), &
1373 1737 : cp_unit_from_cp2k(hmat(2, i), "angstrom"), &
1374 4053 : cp_unit_from_cp2k(hmat(3, i), "angstrom")
1375 : END DO
1376 579 : WRITE (UNIT=file_unit, FMT="(A)") ";"
1377 : ! Data of cell and geometry
1378 : WRITE (UNIT=file_unit, FMT="(/,A,T44,A)") &
1379 579 : "_symmetry_space_group_name_H-M", "'P 1'"
1380 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1381 579 : "_cell_length_a", cp_unit_from_cp2k(abc(1), "angstrom")
1382 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1383 579 : "_cell_length_b", cp_unit_from_cp2k(abc(2), "angstrom")
1384 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1385 579 : "_cell_length_c", cp_unit_from_cp2k(abc(3), "angstrom")
1386 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1387 579 : "_cell_angle_alpha", angle_alpha
1388 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1389 579 : "_cell_angle_beta", angle_beta
1390 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1391 579 : "_cell_angle_gamma", angle_gamma
1392 : WRITE (UNIT=file_unit, FMT="(A,T48,A)") &
1393 579 : "_symmetry_Int_Tables_number", "1"
1394 : WRITE (UNIT=file_unit, FMT="(A,T36,A)") &
1395 579 : "_chemical_formula_structural", formula_structural
1396 : WRITE (UNIT=file_unit, FMT="(A,T36,A)") &
1397 579 : "_chemical_formula_sum", formula_sum
1398 : WRITE (UNIT=file_unit, FMT="(A,T31,F18.8)") &
1399 579 : "_cell_volume", cp_unit_from_cp2k(ABS(deth), "angstrom^3")
1400 : WRITE (UNIT=file_unit, FMT="(A,T41,I8)") &
1401 579 : "_cell_formula_units_Z", gcd_all
1402 : WRITE (UNIT=file_unit, FMT="(A,/,T2,A,/,T2,A,/,T3,A)") &
1403 579 : "loop_", "_symmetry_equiv_pos_site_id", &
1404 1158 : "_symmetry_equiv_pos_as_xyz", "1 'x, y, z'"
1405 : WRITE (UNIT=file_unit, FMT="(A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A,/,T2,A)") &
1406 579 : "loop_", "_atom_site_type_symbol", "_atom_site_label", &
1407 579 : "_atom_site_symmetry_multiplicity", "_atom_site_fract_x", &
1408 1158 : "_atom_site_fract_y", "_atom_site_fract_z", "_atom_site_occupancy"
1409 19321 : DO iatom = 1, natom
1410 18742 : r(1:3) = pbc(particle_set(iatom)%r(1:3), cell)
1411 18742 : CALL real_to_scaled(s, r, cell)
1412 74968 : DO i = 1, 3
1413 74968 : s(i) = MODULO(s(i), 1.0_dp)
1414 : END DO
1415 : WRITE (UNIT=file_unit, FMT=TRIM(f_cif)) &
1416 19321 : cif_type_symbol(iatom), cif_label(iatom), 1, s(1:3), 1.0_dp
1417 : END DO
1418 : END IF
1419 : END IF
1420 :
1421 : ! Finish
1422 0 : DEALLOCATE (element_list, count_list, formula_structural, &
1423 0 : formula_sum, cif_label, cif_type_symbol, &
1424 1073 : xyz_label)
1425 1073 : CALL section_release(tmp_cell_section)
1426 : CALL cp_print_key_finished_output(file_unit, logger, input_section, &
1427 1073 : "PRINT%FINAL_STRUCTURE")
1428 1073 : IF (output_unit > 0) &
1429 : WRITE (UNIT=output_unit, FMT='(/,T2,A)') &
1430 666 : routineN//": Done!"
1431 :
1432 1073 : CALL timestop(handle)
1433 :
1434 4292 : END SUBROUTINE write_final_structure
1435 :
1436 14596 : END MODULE particle_methods
|
