(git:6344619)
Loading...
Searching...
No Matches
reftraj_util.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \brief Initialize the analysis of trajectories to be done
10!> by activating the REFTRAJ ensemble
11!> \par History
12!> Created 10-07 [MI]
13!> \author MI
14! **************************************************************************************************
15MODULE reftraj_util
16
17 USE atomic_kind_list_types, ONLY: atomic_kind_list_type
18 USE atomic_kind_types, ONLY: atomic_kind_type,&
19 get_atomic_kind
20 USE cp_files, ONLY: close_file,&
21 open_file
22 USE cp_log_handling, ONLY: cp_get_default_logger,&
23 cp_logger_type,&
24 cp_to_string
25 USE cp_output_handling, ONLY: cp_print_key_finished_output,&
26 cp_print_key_unit_nr
27 USE cp_parser_methods, ONLY: parser_get_next_line
28 USE cp_subsys_types, ONLY: cp_subsys_get,&
29 cp_subsys_type
30 USE cp_units, ONLY: cp_unit_to_cp2k
31 USE distribution_1d_types, ONLY: distribution_1d_type
32 USE force_env_types, ONLY: force_env_get,&
33 force_env_type
34 USE input_section_types, ONLY: section_vals_get_subs_vals,&
35 section_vals_type,&
36 section_vals_val_get
37 USE kinds, ONLY: default_path_length,&
38 default_string_length,&
39 dp,&
40 max_line_length
41 USE machine, ONLY: m_flush
42 USE md_environment_types, ONLY: get_md_env,&
43 md_environment_type
44 USE message_passing, ONLY: mp_para_env_type
45 USE molecule_kind_list_types, ONLY: molecule_kind_list_type
46 USE molecule_kind_types, ONLY: get_molecule_kind,&
47 molecule_kind_type
48 USE molecule_list_types, ONLY: molecule_list_type
49 USE molecule_types, ONLY: get_molecule,&
50 molecule_type
51 USE particle_list_types, ONLY: particle_list_type
52 USE particle_types, ONLY: particle_type
53 USE physcon, ONLY: angstrom,&
54 femtoseconds
55 USE reftraj_types, ONLY: reftraj_msd_type,&
56 reftraj_type
57 USE simpar_types, ONLY: simpar_type
58 USE string_utilities, ONLY: uppercase
59 USE util, ONLY: get_limit
60#include "../base/base_uses.f90"
61
62 IMPLICIT NONE
63
64 PRIVATE
65
66 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'reftraj_util'
67
68 PUBLIC :: initialize_reftraj, compute_msd_reftraj, write_output_reftraj
69
70CONTAINS
71
72! **************************************************************************************************
73!> \brief ...
74!> \param reftraj ...
75!> \param reftraj_section ...
76!> \param md_env ...
77!> \par History
78!> 10.2007 created
79!> \author MI
80! **************************************************************************************************
81 SUBROUTINE initialize_reftraj(reftraj, reftraj_section, md_env)
82
83 TYPE(reftraj_type), POINTER :: reftraj
84 TYPE(section_vals_type), POINTER :: reftraj_section
85 TYPE(md_environment_type), POINTER :: md_env
86
87 INTEGER :: natom, nline_to_skip, nskip
88 LOGICAL :: my_end
89 TYPE(cp_subsys_type), POINTER :: subsys
90 TYPE(force_env_type), POINTER :: force_env
91 TYPE(mp_para_env_type), POINTER :: para_env
92 TYPE(particle_list_type), POINTER :: particles
93 TYPE(section_vals_type), POINTER :: msd_section
94 TYPE(simpar_type), POINTER :: simpar
95
96 NULLIFY (force_env, msd_section, particles, simpar, subsys)
97 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env, &
98 simpar=simpar)
99 CALL force_env_get(force_env=force_env, subsys=subsys)
100 CALL cp_subsys_get(subsys=subsys, particles=particles)
101 natom = particles%n_els
102
103 my_end = .false.
104 nline_to_skip = 0
105
106 nskip = reftraj%info%first_snapshot - 1
107 cpassert(nskip >= 0)
108
109 IF (nskip > 0) THEN
110 nline_to_skip = (natom + 2)*nskip
111 CALL parser_get_next_line(reftraj%info%traj_parser, nline_to_skip, at_end=my_end)
112 END IF
113
114 reftraj%isnap = nskip
115 IF (my_end) THEN
116 CALL cp_abort(__location__, &
117 "Reached the end of the trajectory file for REFTRAJ. Number of steps skipped "// &
118 "equal to the number of steps present in the file.")
119 END IF
120
121 ! Cell File
122 IF (reftraj%info%variable_volume) THEN
123 IF (nskip > 0) THEN
124 CALL parser_get_next_line(reftraj%info%cell_parser, nskip, at_end=my_end)
125 END IF
126 IF (my_end) THEN
127 CALL cp_abort(__location__, &
128 "Reached the end of the cell file for REFTRAJ. Number of steps skipped "// &
129 "equal to the number of steps present in the file.")
130 END IF
131 END IF
132
133 reftraj%natom = natom
134 IF (reftraj%info%last_snapshot > 0) THEN
135 simpar%nsteps = (reftraj%info%last_snapshot - reftraj%info%first_snapshot + 1)
136 END IF
137
138 IF (reftraj%info%msd) THEN
139 msd_section => section_vals_get_subs_vals(reftraj_section, "MSD")
140 ! set up and printout
141 CALL initialize_msd_reftraj(reftraj%msd, msd_section, reftraj, md_env)
142 END IF
143
144 END SUBROUTINE initialize_reftraj
145
146! **************************************************************************************************
147!> \brief ...
148!> \param msd ...
149!> \param msd_section ...
150!> \param reftraj ...
151!> \param md_env ...
152!> \par History
153!> 10.2007 created
154!> \author MI
155! **************************************************************************************************
156 SUBROUTINE initialize_msd_reftraj(msd, msd_section, reftraj, md_env)
157 TYPE(reftraj_msd_type), POINTER :: msd
158 TYPE(section_vals_type), POINTER :: msd_section
159 TYPE(reftraj_type), POINTER :: reftraj
160 TYPE(md_environment_type), POINTER :: md_env
161
162 CHARACTER(LEN=2) :: element_symbol, element_symbol_ref0
163 CHARACTER(LEN=default_path_length) :: filename
164 CHARACTER(LEN=default_string_length) :: title
165 CHARACTER(LEN=max_line_length) :: errmsg
166 INTEGER :: first_atom, iatom, ikind, imol, &
167 last_atom, natom_read, nkind, nmol, &
168 nmolecule, nmolkind, npart
169 REAL(kind=dp) :: com(3), mass, mass_mol, tol, x, y, z
170 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
171 TYPE(cp_subsys_type), POINTER :: subsys
172 TYPE(force_env_type), POINTER :: force_env
173 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
174 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
175 TYPE(molecule_kind_type), POINTER :: molecule_kind
176 TYPE(molecule_list_type), POINTER :: molecules
177 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
178 TYPE(molecule_type), POINTER :: molecule
179 TYPE(mp_para_env_type), POINTER :: para_env
180 TYPE(particle_list_type), POINTER :: particles
181 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
182
183 NULLIFY (molecule, molecules, molecule_kind, molecule_kind_set, &
184 molecule_kinds, molecule_set, subsys, force_env, particles, particle_set)
185 cpassert(.NOT. ASSOCIATED(msd))
186
187 ALLOCATE (msd)
188
189 NULLIFY (msd%ref0_pos)
190 NULLIFY (msd%ref0_com_molecule)
191 NULLIFY (msd%val_msd_kind)
192 NULLIFY (msd%val_msd_molecule)
193 NULLIFY (msd%disp_atom_index)
194 NULLIFY (msd%disp_atom_dr)
195
196 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
197 CALL force_env_get(force_env=force_env, subsys=subsys)
198 CALL cp_subsys_get(subsys=subsys, particles=particles)
199 particle_set => particles%els
200 npart = SIZE(particle_set, 1)
201
202 msd%ref0_unit = -1
203 CALL section_vals_val_get(msd_section, "REF0_FILENAME", c_val=filename)
204 CALL open_file(trim(filename), unit_number=msd%ref0_unit)
205
206 ALLOCATE (msd%ref0_pos(3, reftraj%natom))
207 msd%ref0_pos = 0.0_dp
208
209 IF (para_env%is_source()) THEN
210 rewind(msd%ref0_unit)
211 READ (msd%ref0_unit, *, err=999, END=998) natom_read
212 IF (natom_read /= reftraj%natom) THEN
213 errmsg = "The MSD reference configuration has a different number of atoms: "// &
214 trim(adjustl(cp_to_string(natom_read)))//" != "// &
215 trim(adjustl(cp_to_string(reftraj%natom)))
216 cpabort(errmsg)
217 END IF
218 READ (msd%ref0_unit, '(A)', err=999, END=998) title
219 msd%total_mass = 0.0_dp
220 msd%ref0_com = 0.0_dp
221 DO iatom = 1, natom_read
222 READ (msd%ref0_unit, *, err=999, END=998) element_symbol_ref0, x, y, z
223 CALL uppercase(element_symbol_ref0)
224 element_symbol = trim(particle_set(iatom)%atomic_kind%element_symbol)
225 CALL uppercase(element_symbol)
226 IF (element_symbol /= element_symbol_ref0) THEN
227 errmsg = "The MSD reference configuration shows a mismatch: Check atom "// &
228 trim(adjustl(cp_to_string(iatom)))
229 cpabort(errmsg)
230 END IF
231 x = cp_unit_to_cp2k(x, "angstrom")
232 y = cp_unit_to_cp2k(y, "angstrom")
233 z = cp_unit_to_cp2k(z, "angstrom")
234 msd%ref0_pos(1, iatom) = x
235 msd%ref0_pos(2, iatom) = y
236 msd%ref0_pos(3, iatom) = z
237 mass = particle_set(iatom)%atomic_kind%mass
238 msd%ref0_com(1) = msd%ref0_com(1) + x*mass
239 msd%ref0_com(2) = msd%ref0_com(2) + y*mass
240 msd%ref0_com(3) = msd%ref0_com(3) + z*mass
241 msd%total_mass = msd%total_mass + mass
242 END DO
243 msd%ref0_com = msd%ref0_com/msd%total_mass
244 END IF
245 CALL close_file(unit_number=msd%ref0_unit)
246
247 CALL para_env%bcast(msd%total_mass)
248 CALL para_env%bcast(msd%ref0_pos)
249 CALL para_env%bcast(msd%ref0_com)
250
251 CALL section_vals_val_get(msd_section, "MSD_PER_KIND", l_val=msd%msd_kind)
252 CALL section_vals_val_get(msd_section, "MSD_PER_MOLKIND", l_val=msd%msd_molecule)
253 CALL section_vals_val_get(msd_section, "MSD_PER_REGION", l_val=msd%msd_region)
254
255 CALL section_vals_val_get(msd_section, "DISPLACED_ATOM", l_val=msd%disp_atom)
256 IF (msd%disp_atom) THEN
257 ALLOCATE (msd%disp_atom_index(npart))
258 msd%disp_atom_index = 0
259 ALLOCATE (msd%disp_atom_dr(3, npart))
260 msd%disp_atom_dr = 0.0_dp
261 msd%msd_kind = .true.
262 END IF
263 CALL section_vals_val_get(msd_section, "DISPLACEMENT_TOL", r_val=tol)
264 msd%disp_atom_tol = tol*tol
265
266 IF (msd%msd_kind) THEN
267 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
268 nkind = atomic_kinds%n_els
269
270 ALLOCATE (msd%val_msd_kind(4, nkind))
271 msd%val_msd_kind = 0.0_dp
272 END IF
273
274 IF (msd%msd_molecule) THEN
275 CALL cp_subsys_get(subsys=subsys, molecules=molecules, &
276 molecule_kinds=molecule_kinds)
277 nmolkind = molecule_kinds%n_els
278 ALLOCATE (msd%val_msd_molecule(4, nmolkind))
279
280 molecule_kind_set => molecule_kinds%els
281 molecule_set => molecules%els
282 nmol = molecules%n_els
283
284 ALLOCATE (msd%ref0_com_molecule(3, nmol))
285
286 DO ikind = 1, nmolkind
287 molecule_kind => molecule_kind_set(ikind)
288 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
289 DO imol = 1, nmolecule
290 molecule => molecule_set(molecule_kind%molecule_list(imol))
291 CALL get_molecule(molecule=molecule, first_atom=first_atom, last_atom=last_atom)
292 com = 0.0_dp
293 mass_mol = 0.0_dp
294 DO iatom = first_atom, last_atom
295 mass = particle_set(iatom)%atomic_kind%mass
296 com(1) = com(1) + msd%ref0_pos(1, iatom)*mass
297 com(2) = com(2) + msd%ref0_pos(2, iatom)*mass
298 com(3) = com(3) + msd%ref0_pos(3, iatom)*mass
299 mass_mol = mass_mol + mass
300 END DO ! iatom
301 msd%ref0_com_molecule(1, molecule_kind%molecule_list(imol)) = com(1)/mass_mol
302 msd%ref0_com_molecule(2, molecule_kind%molecule_list(imol)) = com(2)/mass_mol
303 msd%ref0_com_molecule(3, molecule_kind%molecule_list(imol)) = com(3)/mass_mol
304 END DO ! imol
305 END DO ! ikind
306 END IF
307
308 IF (msd%msd_region) THEN
309
310 END IF
311
312 RETURN
313998 CONTINUE ! end of file
314 cpabort("End of reference positions file reached")
315999 CONTINUE ! error
316 cpabort("Error reading reference positions file")
317
318 END SUBROUTINE initialize_msd_reftraj
319
320! **************************************************************************************************
321!> \brief ...
322!> \param reftraj ...
323!> \param md_env ...
324!> \param particle_set ...
325!> \par History
326!> 10.2007 created
327!> \author MI
328! **************************************************************************************************
329 SUBROUTINE compute_msd_reftraj(reftraj, md_env, particle_set)
330
331 TYPE(reftraj_type), POINTER :: reftraj
332 TYPE(md_environment_type), POINTER :: md_env
333 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
334
335 INTEGER :: atom, bo(2), first_atom, iatom, ikind, imol, imol_global, last_atom, mepos, &
336 natom_kind, nmol_per_kind, nmolecule, nmolkind, num_pe
337 INTEGER, DIMENSION(:), POINTER :: atom_list
338 REAL(kind=dp) :: com(3), diff2_com(4), dr2, dx, dy, dz, &
339 mass, mass_mol, msd_mkind(4), rcom(3)
340 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
341 TYPE(atomic_kind_type), POINTER :: atomic_kind
342 TYPE(cp_subsys_type), POINTER :: subsys
343 TYPE(distribution_1d_type), POINTER :: local_molecules
344 TYPE(force_env_type), POINTER :: force_env
345 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
346 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
347 TYPE(molecule_kind_type), POINTER :: molecule_kind
348 TYPE(molecule_list_type), POINTER :: molecules
349 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
350 TYPE(molecule_type), POINTER :: molecule
351 TYPE(mp_para_env_type), POINTER :: para_env
352
353 NULLIFY (force_env, para_env, subsys)
354 NULLIFY (atomic_kind, atomic_kinds, atom_list)
355 NULLIFY (local_molecules, molecule, molecule_kind, molecule_kinds, &
356 molecule_kind_set, molecules, molecule_set)
357
358 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
359 CALL force_env_get(force_env=force_env, subsys=subsys)
360 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
361
362 num_pe = para_env%num_pe
363 mepos = para_env%mepos
364
365 IF (reftraj%msd%msd_kind) THEN
366 reftraj%msd%val_msd_kind = 0.0_dp
367 reftraj%msd%num_disp_atom = 0
368 reftraj%msd%disp_atom_dr = 0.0_dp
369! compute com
370 rcom = 0.0_dp
371 DO ikind = 1, atomic_kinds%n_els
372 atomic_kind => atomic_kinds%els(ikind)
373 CALL get_atomic_kind(atomic_kind=atomic_kind, &
374 atom_list=atom_list, &
375 natom=natom_kind, mass=mass)
376 bo = get_limit(natom_kind, num_pe, mepos)
377 DO iatom = bo(1), bo(2)
378 atom = atom_list(iatom)
379 rcom(1) = rcom(1) + particle_set(atom)%r(1)*mass
380 rcom(2) = rcom(2) + particle_set(atom)%r(2)*mass
381 rcom(3) = rcom(3) + particle_set(atom)%r(3)*mass
382 END DO
383 END DO
384 CALL para_env%sum(rcom)
385 rcom = rcom/reftraj%msd%total_mass
386 reftraj%msd%drcom(1) = rcom(1) - reftraj%msd%ref0_com(1)
387 reftraj%msd%drcom(2) = rcom(2) - reftraj%msd%ref0_com(2)
388 reftraj%msd%drcom(3) = rcom(3) - reftraj%msd%ref0_com(3)
389! IF(para_env%is_source()) WRITE(*,'(A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', &
390! drcom(1)*angstrom,drcom(2)*angstrom,drcom(3)*angstrom
391! compute_com
392
393 DO ikind = 1, atomic_kinds%n_els
394 atomic_kind => atomic_kinds%els(ikind)
395 CALL get_atomic_kind(atomic_kind=atomic_kind, &
396 atom_list=atom_list, &
397 natom=natom_kind)
398 bo = get_limit(natom_kind, num_pe, mepos)
399 DO iatom = bo(1), bo(2)
400 atom = atom_list(iatom)
401 dx = particle_set(atom)%r(1) - reftraj%msd%ref0_pos(1, atom) - &
402 reftraj%msd%drcom(1)
403 dy = particle_set(atom)%r(2) - reftraj%msd%ref0_pos(2, atom) - &
404 reftraj%msd%drcom(2)
405 dz = particle_set(atom)%r(3) - reftraj%msd%ref0_pos(3, atom) - &
406 reftraj%msd%drcom(3)
407 dr2 = dx*dx + dy*dy + dz*dz
408
409 reftraj%msd%val_msd_kind(1, ikind) = reftraj%msd%val_msd_kind(1, ikind) + dx*dx
410 reftraj%msd%val_msd_kind(2, ikind) = reftraj%msd%val_msd_kind(2, ikind) + dy*dy
411 reftraj%msd%val_msd_kind(3, ikind) = reftraj%msd%val_msd_kind(3, ikind) + dz*dz
412 reftraj%msd%val_msd_kind(4, ikind) = reftraj%msd%val_msd_kind(4, ikind) + dr2
413
414 IF (reftraj%msd%disp_atom) THEN
415 IF (dr2 > reftraj%msd%disp_atom_tol) THEN
416 reftraj%msd%num_disp_atom = reftraj%msd%num_disp_atom + 1
417 reftraj%msd%disp_atom_dr(1, atom) = dx
418 reftraj%msd%disp_atom_dr(2, atom) = dy
419 reftraj%msd%disp_atom_dr(3, atom) = dz
420 END IF
421 END IF
422 END DO !iatom
423 reftraj%msd%val_msd_kind(1:4, ikind) = &
424 reftraj%msd%val_msd_kind(1:4, ikind)/real(natom_kind, kind=dp)
425
426 END DO ! ikind
427 END IF
428 CALL para_env%sum(reftraj%msd%val_msd_kind)
429 CALL para_env%sum(reftraj%msd%num_disp_atom)
430 CALL para_env%sum(reftraj%msd%disp_atom_dr)
431
432 IF (reftraj%msd%msd_molecule) THEN
433 CALL cp_subsys_get(subsys=subsys, local_molecules=local_molecules, &
434 molecules=molecules, molecule_kinds=molecule_kinds)
435
436 nmolkind = molecule_kinds%n_els
437 molecule_kind_set => molecule_kinds%els
438 molecule_set => molecules%els
439
440 reftraj%msd%val_msd_molecule = 0.0_dp
441 DO ikind = 1, nmolkind
442 molecule_kind => molecule_kind_set(ikind)
443 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
444 nmol_per_kind = local_molecules%n_el(ikind)
445 msd_mkind = 0.0_dp
446 DO imol = 1, nmol_per_kind
447 imol_global = local_molecules%list(ikind)%array(imol)
448 molecule => molecule_set(imol_global)
449 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
450
451 com = 0.0_dp
452 mass_mol = 0.0_dp
453 DO iatom = first_atom, last_atom
454 mass = particle_set(iatom)%atomic_kind%mass
455 com(1) = com(1) + particle_set(iatom)%r(1)*mass
456 com(2) = com(2) + particle_set(iatom)%r(2)*mass
457 com(3) = com(3) + particle_set(iatom)%r(3)*mass
458 mass_mol = mass_mol + mass
459 END DO ! iatom
460 com(1) = com(1)/mass_mol
461 com(2) = com(2)/mass_mol
462 com(3) = com(3)/mass_mol
463 diff2_com(1) = com(1) - reftraj%msd%ref0_com_molecule(1, imol_global)
464 diff2_com(2) = com(2) - reftraj%msd%ref0_com_molecule(2, imol_global)
465 diff2_com(3) = com(3) - reftraj%msd%ref0_com_molecule(3, imol_global)
466 diff2_com(1) = diff2_com(1)*diff2_com(1)
467 diff2_com(2) = diff2_com(2)*diff2_com(2)
468 diff2_com(3) = diff2_com(3)*diff2_com(3)
469 diff2_com(4) = diff2_com(1) + diff2_com(2) + diff2_com(3)
470 msd_mkind(1) = msd_mkind(1) + diff2_com(1)
471 msd_mkind(2) = msd_mkind(2) + diff2_com(2)
472 msd_mkind(3) = msd_mkind(3) + diff2_com(3)
473 msd_mkind(4) = msd_mkind(4) + diff2_com(4)
474 END DO ! imol
475
476 reftraj%msd%val_msd_molecule(1, ikind) = msd_mkind(1)/real(nmolecule, kind=dp)
477 reftraj%msd%val_msd_molecule(2, ikind) = msd_mkind(2)/real(nmolecule, kind=dp)
478 reftraj%msd%val_msd_molecule(3, ikind) = msd_mkind(3)/real(nmolecule, kind=dp)
479 reftraj%msd%val_msd_molecule(4, ikind) = msd_mkind(4)/real(nmolecule, kind=dp)
480 END DO ! ikind
481 CALL para_env%sum(reftraj%msd%val_msd_molecule)
482
483 END IF
484
485 END SUBROUTINE compute_msd_reftraj
486
487! **************************************************************************************************
488!> \brief ...
489!> \param md_env ...
490!> \par History
491!> 10.2007 created
492!> \author MI
493! **************************************************************************************************
494 SUBROUTINE write_output_reftraj(md_env)
495 TYPE(md_environment_type), POINTER :: md_env
496
497 CHARACTER(LEN=default_string_length) :: my_act, my_mittle, my_pos
498 INTEGER :: iat, ikind, nkind, out_msd
499 LOGICAL, SAVE :: first_entry = .false.
500 TYPE(cp_logger_type), POINTER :: logger
501 TYPE(force_env_type), POINTER :: force_env
502 TYPE(reftraj_type), POINTER :: reftraj
503 TYPE(section_vals_type), POINTER :: reftraj_section, root_section
504
505 NULLIFY (logger)
506 logger => cp_get_default_logger()
507
508 NULLIFY (reftraj)
509 NULLIFY (reftraj_section, root_section)
510
511 CALL get_md_env(md_env=md_env, force_env=force_env, &
512 reftraj=reftraj)
513
514 CALL force_env_get(force_env=force_env, root_section=root_section)
515
516 reftraj_section => section_vals_get_subs_vals(root_section, &
517 "MOTION%MD%REFTRAJ")
518
519 my_pos = "APPEND"
520 my_act = "WRITE"
521
522 IF (reftraj%init .AND. (reftraj%isnap == reftraj%info%first_snapshot)) THEN
523 my_pos = "REWIND"
524 first_entry = .true.
525 END IF
526
527 IF (reftraj%info%msd) THEN
528 IF (reftraj%msd%msd_kind) THEN
529 nkind = SIZE(reftraj%msd%val_msd_kind, 2)
530 DO ikind = 1, nkind
531 my_mittle = "k"//trim(adjustl(cp_to_string(ikind)))
532 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_KIND", &
533 extension=".msd", file_position=my_pos, file_action=my_act, &
534 file_form="FORMATTED", middle_name=trim(my_mittle))
535 IF (out_msd > 0) THEN
536 WRITE (unit=out_msd, fmt="(I8, F12.3,4F20.10)") reftraj%itimes, &
537 reftraj%time*femtoseconds, &
538 reftraj%msd%val_msd_kind(1:4, ikind)*angstrom*angstrom
539 CALL m_flush(out_msd)
540 END IF
541 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
542 "PRINT%MSD_KIND")
543 END DO
544 END IF
545 IF (reftraj%msd%msd_molecule) THEN
546 nkind = SIZE(reftraj%msd%val_msd_molecule, 2)
547 DO ikind = 1, nkind
548 my_mittle = "mk"//trim(adjustl(cp_to_string(ikind)))
549 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_MOLECULE", &
550 extension=".msd", file_position=my_pos, file_action=my_act, &
551 file_form="FORMATTED", middle_name=trim(my_mittle))
552 IF (out_msd > 0) THEN
553 WRITE (unit=out_msd, fmt="(I8, F12.3,4F20.10)") reftraj%itimes, &
554 reftraj%time*femtoseconds, &
555 reftraj%msd%val_msd_molecule(1:4, ikind)*angstrom*angstrom
556 CALL m_flush(out_msd)
557 END IF
558 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
559 "PRINT%MSD_MOLECULE")
560 END DO
561 END IF
562 IF (reftraj%msd%disp_atom) THEN
563
564 IF (first_entry) my_pos = "REWIND"
565 my_mittle = "disp_at"
566 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%DISPLACED_ATOM", &
567 extension=".msd", file_position=my_pos, file_action=my_act, &
568 file_form="FORMATTED", middle_name=trim(my_mittle))
569 IF (out_msd > 0 .AND. reftraj%msd%num_disp_atom > 0) THEN
570 IF (first_entry) THEN
571 first_entry = .false.
572 END IF
573 WRITE (unit=out_msd, fmt="(A,T7,I8, A, T29, F12.3, A, T50, I10)") "# i = ", reftraj%itimes, " time (fs) = ", &
574 reftraj%time*femtoseconds, " nat = ", reftraj%msd%num_disp_atom
575 DO iat = 1, SIZE(reftraj%msd%disp_atom_dr, 2)
576 IF (abs(reftraj%msd%disp_atom_dr(1, iat)) > 0.0_dp) THEN
577 WRITE (unit=out_msd, fmt="(I8, 3F20.10)") iat, & !reftraj%msd%disp_atom_index(iat),&
578 reftraj%msd%disp_atom_dr(1, iat)*angstrom, &
579 reftraj%msd%disp_atom_dr(2, iat)*angstrom, &
580 reftraj%msd%disp_atom_dr(3, iat)*angstrom
581 END IF
582 END DO
583 END IF
584 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
585 "PRINT%DISPLACED_ATOM")
586 END IF
587 END IF ! msd
588 reftraj%init = .false.
589
590 END SUBROUTINE write_output_reftraj
591
592END MODULE reftraj_util
593
cp_log_handling::cp_to_string
Definition cp_log_handling.F:90
atom
Definition atom.F:9
atomic_kind_list_types
represent a simple array based list of the given type
Definition atomic_kind_list_types.F:15
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
cp_files
Utility routines to open and close files.
Definition cp_files.F:22
cp_files::open_file
subroutine, public open_file(file_name, unit_number, file_status, file_form, file_action, file_position, file_pad, file_access, debug)
Opens a file and returns its logical unit number.
Definition cp_files.F:197
cp_files::close_file
subroutine, public close_file(unit_number, file_status)
Closes an open file given by its logical unit number.
Definition cp_files.F:64
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:838
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1058
cp_parser_methods
Utility routines to read data from files. Kept as close as possible to the old parser because.
Definition cp_parser_methods.F:20
cp_parser_methods::parser_get_next_line
subroutine, public parser_get_next_line(parser, nline, at_end)
Read the next input line and broadcast the input information. Skip (nline-1) lines and skip also all ...
Definition cp_parser_methods.F:421
cp_subsys_types
types that represent a subsys, i.e. a part of the system
Definition cp_subsys_types.F:16
cp_subsys_types::cp_subsys_get
subroutine, public cp_subsys_get(subsys, ref_count, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:364
cp_units
unit conversion facility
Definition cp_units.F:30
cp_units::cp_unit_to_cp2k
real(kind=dp) function, public cp_unit_to_cp2k(value, unit_str, defaults, power)
converts to the internal cp2k units to the given unit
Definition cp_units.F:1282
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
force_env_types
Interface for the force calculations.
Definition force_env_types.F:18
force_env_types::force_env_get
recursive subroutine, public force_env_get(force_env, in_use, fist_env, qs_env, meta_env, fp_env, subsys, para_env, potential_energy, additional_potential, kinetic_energy, harmonic_shell, kinetic_shell, cell, sub_force_env, qmmm_env, qmmmx_env, eip_env, pwdft_env, globenv, input, force_env_section, method_name_id, root_section, mixed_env, nnp_env, embed_env, ipi_env)
returns various attributes about the force environment
Definition force_env_types.F:344
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition input_section_types.F:735
input_section_types::section_vals_val_get
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Definition input_section_types.F:1052
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::max_line_length
integer, parameter, public max_line_length
Definition kinds.F:59
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
kinds::default_path_length
integer, parameter, public default_path_length
Definition kinds.F:58
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_flush
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition machine.F:127
md_environment_types
Definition md_environment_types.F:15
md_environment_types::get_md_env
subroutine, public get_md_env(md_env, itimes, constant, used_time, cell, simpar, npt, force_env, para_env, reftraj, t, init, first_time, fe_env, thermostats, barostat, thermostat_coeff, thermostat_part, thermostat_shell, thermostat_baro, thermostat_fast, thermostat_slow, md_ener, averages, thermal_regions, ehrenfest_md)
get components of MD environment type
Definition md_environment_types.F:194
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:27
molecule_kind_list_types
represent a simple array based list of the given type
Definition molecule_kind_list_types.F:15
molecule_kind_types
Define the molecule kind structure types and the corresponding functionality.
Definition molecule_kind_types.F:16
molecule_kind_types::get_molecule_kind
subroutine, public get_molecule_kind(molecule_kind, atom_list, bond_list, bend_list, ub_list, impr_list, opbend_list, colv_list, fixd_list, g3x3_list, g4x6_list, vsite_list, torsion_list, shell_list, name, mass, charge, kind_number, natom, nbend, nbond, nub, nimpr, nopbend, nconstraint, nconstraint_fixd, nfixd, ncolv, ng3x3, ng4x6, nvsite, nfixd_restraint, ng3x3_restraint, ng4x6_restraint, nvsite_restraint, nrestraints, nmolecule, nsgf, nshell, ntorsion, molecule_list, nelectron, nelectron_alpha, nelectron_beta, bond_kind_set, bend_kind_set, ub_kind_set, impr_kind_set, opbend_kind_set, torsion_kind_set, molname_generated)
Get informations about a molecule kind.
Definition molecule_kind_types.F:518
molecule_list_types
represent a simple array based list of the given type
Definition molecule_list_types.F:15
molecule_types
Define the data structure for the molecule information.
Definition molecule_types.F:16
molecule_types::get_molecule
subroutine, public get_molecule(molecule, molecule_kind, lmi, lci, lg3x3, lg4x6, lcolv, first_atom, last_atom, first_shell, last_shell)
Get components from a molecule data set.
Definition molecule_types.F:333
particle_list_types
represent a simple array based list of the given type
Definition particle_list_types.F:15
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
physcon
Definition of physical constants:
Definition physcon.F:68
physcon::femtoseconds
real(kind=dp), parameter, public femtoseconds
Definition physcon.F:153
physcon::angstrom
real(kind=dp), parameter, public angstrom
Definition physcon.F:144
reftraj_types
initialization of the reftraj structure used to analyse previously generated trajectories
Definition reftraj_types.F:15
reftraj_util
Initialize the analysis of trajectories to be done by activating the REFTRAJ ensemble.
Definition reftraj_util.F:15
reftraj_util::write_output_reftraj
subroutine, public write_output_reftraj(md_env)
...
Definition reftraj_util.F:495
reftraj_util::compute_msd_reftraj
subroutine, public compute_msd_reftraj(reftraj, md_env, particle_set)
...
Definition reftraj_util.F:330
reftraj_util::initialize_reftraj
subroutine, public initialize_reftraj(reftraj, reftraj_section, md_env)
...
Definition reftraj_util.F:82
simpar_types
Type for storing MD parameters.
Definition simpar_types.F:14
string_utilities
Utilities for string manipulations.
Definition string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition string_utilities.F:3362
util
All kind of helpful little routines.
Definition util.F:14
util::get_limit
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition util.F:333
atomic_kind_list_types::atomic_kind_list_type
represent a list of objects
Definition atomic_kind_list_types.F:45
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
cp_subsys_types::cp_subsys_type
represents a system: atoms, molecules, their pos,vel,...
Definition cp_subsys_types.F:89
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
force_env_types::force_env_type
wrapper to abstract the force evaluation of the various methods
Definition force_env_types.F:145
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
md_environment_types::md_environment_type
Definition md_environment_types.F:52
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:768
molecule_kind_list_types::molecule_kind_list_type
represent a list of objects
Definition molecule_kind_list_types.F:45
molecule_kind_types::molecule_kind_type
Definition molecule_kind_types.F:153
molecule_list_types::molecule_list_type
represent a list of objects
Definition molecule_list_types.F:46
molecule_types::molecule_type
Definition molecule_types.F:160
particle_list_types::particle_list_type
represent a list of objects
Definition particle_list_types.F:46
particle_types::particle_type
Definition particle_types.F:35
reftraj_types::reftraj_msd_type
Definition reftraj_types.F:57
reftraj_types::reftraj_type
Definition reftraj_types.F:71
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30