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