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