(git:ed6f26b)
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-2025 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) &
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
120 ! Cell File
121 IF (reftraj%info%variable_volume) THEN
122 IF (nskip > 0) THEN
123 CALL parser_get_next_line(reftraj%info%cell_parser, nskip, at_end=my_end)
124 END IF
125 IF (my_end) &
126 CALL cp_abort(__location__, &
127 "Reached the end of the cell file for REFTRAJ. Number of steps skipped "// &
128 "equal to the number of steps present in the file.")
129 END IF
130
131 reftraj%natom = natom
132 IF (reftraj%info%last_snapshot > 0) THEN
133 simpar%nsteps = (reftraj%info%last_snapshot - reftraj%info%first_snapshot + 1)
134 END IF
135
136 IF (reftraj%info%msd) THEN
137 msd_section => section_vals_get_subs_vals(reftraj_section, "MSD")
138 ! set up and printout
139 CALL initialize_msd_reftraj(reftraj%msd, msd_section, reftraj, md_env)
140 END IF
141
142 END SUBROUTINE initialize_reftraj
143
144! **************************************************************************************************
145!> \brief ...
146!> \param msd ...
147!> \param msd_section ...
148!> \param reftraj ...
149!> \param md_env ...
150!> \par History
151!> 10.2007 created
152!> \author MI
153! **************************************************************************************************
154 SUBROUTINE initialize_msd_reftraj(msd, msd_section, reftraj, md_env)
155 TYPE(reftraj_msd_type), POINTER :: msd
156 TYPE(section_vals_type), POINTER :: msd_section
157 TYPE(reftraj_type), POINTER :: reftraj
158 TYPE(md_environment_type), POINTER :: md_env
159
160 CHARACTER(LEN=2) :: element_symbol, element_symbol_ref0
161 CHARACTER(LEN=default_path_length) :: filename
162 CHARACTER(LEN=default_string_length) :: title
163 CHARACTER(LEN=max_line_length) :: errmsg
164 INTEGER :: first_atom, iatom, ikind, imol, &
165 last_atom, natom_read, nkind, nmol, &
166 nmolecule, nmolkind, npart
167 REAL(kind=dp) :: com(3), mass, mass_mol, tol, x, y, z
168 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
169 TYPE(cp_subsys_type), POINTER :: subsys
170 TYPE(force_env_type), POINTER :: force_env
171 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
172 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
173 TYPE(molecule_kind_type), POINTER :: molecule_kind
174 TYPE(molecule_list_type), POINTER :: molecules
175 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
176 TYPE(molecule_type), POINTER :: molecule
177 TYPE(mp_para_env_type), POINTER :: para_env
178 TYPE(particle_list_type), POINTER :: particles
179 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
180
181 NULLIFY (molecule, molecules, molecule_kind, molecule_kind_set, &
182 molecule_kinds, molecule_set, subsys, force_env, particles, particle_set)
183 cpassert(.NOT. ASSOCIATED(msd))
184
185 ALLOCATE (msd)
186
187 NULLIFY (msd%ref0_pos)
188 NULLIFY (msd%ref0_com_molecule)
189 NULLIFY (msd%val_msd_kind)
190 NULLIFY (msd%val_msd_molecule)
191 NULLIFY (msd%disp_atom_index)
192 NULLIFY (msd%disp_atom_dr)
193
194 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
195 CALL force_env_get(force_env=force_env, subsys=subsys)
196 CALL cp_subsys_get(subsys=subsys, particles=particles)
197 particle_set => particles%els
198 npart = SIZE(particle_set, 1)
199
200 msd%ref0_unit = -1
201 CALL section_vals_val_get(msd_section, "REF0_FILENAME", c_val=filename)
202 CALL open_file(trim(filename), unit_number=msd%ref0_unit)
203
204 ALLOCATE (msd%ref0_pos(3, reftraj%natom))
205 msd%ref0_pos = 0.0_dp
206
207 IF (para_env%is_source()) THEN
208 rewind(msd%ref0_unit)
209 READ (msd%ref0_unit, *, err=999, END=998) natom_read
210 IF (natom_read /= reftraj%natom) THEN
211 errmsg = "The MSD reference configuration has a different number of atoms: "// &
212 trim(adjustl(cp_to_string(natom_read)))//" != "// &
213 trim(adjustl(cp_to_string(reftraj%natom)))
214 cpabort(errmsg)
215 END IF
216 READ (msd%ref0_unit, '(A)', err=999, END=998) title
217 msd%total_mass = 0.0_dp
218 msd%ref0_com = 0.0_dp
219 DO iatom = 1, natom_read
220 READ (msd%ref0_unit, *, err=999, END=998) element_symbol_ref0, x, y, z
221 CALL uppercase(element_symbol_ref0)
222 element_symbol = trim(particle_set(iatom)%atomic_kind%element_symbol)
223 CALL uppercase(element_symbol)
224 IF (element_symbol /= element_symbol_ref0) THEN
225 errmsg = "The MSD reference configuration shows a mismatch: Check atom "// &
226 trim(adjustl(cp_to_string(iatom)))
227 cpabort(errmsg)
228 END IF
229 x = cp_unit_to_cp2k(x, "angstrom")
230 y = cp_unit_to_cp2k(y, "angstrom")
231 z = cp_unit_to_cp2k(z, "angstrom")
232 msd%ref0_pos(1, iatom) = x
233 msd%ref0_pos(2, iatom) = y
234 msd%ref0_pos(3, iatom) = z
235 mass = particle_set(iatom)%atomic_kind%mass
236 msd%ref0_com(1) = msd%ref0_com(1) + x*mass
237 msd%ref0_com(2) = msd%ref0_com(2) + y*mass
238 msd%ref0_com(3) = msd%ref0_com(3) + z*mass
239 msd%total_mass = msd%total_mass + mass
240 END DO
241 msd%ref0_com = msd%ref0_com/msd%total_mass
242 END IF
243 CALL close_file(unit_number=msd%ref0_unit)
244
245 CALL para_env%bcast(msd%total_mass)
246 CALL para_env%bcast(msd%ref0_pos)
247 CALL para_env%bcast(msd%ref0_com)
248
249 CALL section_vals_val_get(msd_section, "MSD_PER_KIND", l_val=msd%msd_kind)
250 CALL section_vals_val_get(msd_section, "MSD_PER_MOLKIND", l_val=msd%msd_molecule)
251 CALL section_vals_val_get(msd_section, "MSD_PER_REGION", l_val=msd%msd_region)
252
253 CALL section_vals_val_get(msd_section, "DISPLACED_ATOM", l_val=msd%disp_atom)
254 IF (msd%disp_atom) THEN
255 ALLOCATE (msd%disp_atom_index(npart))
256 msd%disp_atom_index = 0
257 ALLOCATE (msd%disp_atom_dr(3, npart))
258 msd%disp_atom_dr = 0.0_dp
259 msd%msd_kind = .true.
260 END IF
261 CALL section_vals_val_get(msd_section, "DISPLACEMENT_TOL", r_val=tol)
262 msd%disp_atom_tol = tol*tol
263
264 IF (msd%msd_kind) THEN
265 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
266 nkind = atomic_kinds%n_els
267
268 ALLOCATE (msd%val_msd_kind(4, nkind))
269 msd%val_msd_kind = 0.0_dp
270 END IF
271
272 IF (msd%msd_molecule) THEN
273 CALL cp_subsys_get(subsys=subsys, molecules=molecules, &
274 molecule_kinds=molecule_kinds)
275 nmolkind = molecule_kinds%n_els
276 ALLOCATE (msd%val_msd_molecule(4, nmolkind))
277
278 molecule_kind_set => molecule_kinds%els
279 molecule_set => molecules%els
280 nmol = molecules%n_els
281
282 ALLOCATE (msd%ref0_com_molecule(3, nmol))
283
284 DO ikind = 1, nmolkind
285 molecule_kind => molecule_kind_set(ikind)
286 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
287 DO imol = 1, nmolecule
288 molecule => molecule_set(molecule_kind%molecule_list(imol))
289 CALL get_molecule(molecule=molecule, first_atom=first_atom, last_atom=last_atom)
290 com = 0.0_dp
291 mass_mol = 0.0_dp
292 DO iatom = first_atom, last_atom
293 mass = particle_set(iatom)%atomic_kind%mass
294 com(1) = com(1) + msd%ref0_pos(1, iatom)*mass
295 com(2) = com(2) + msd%ref0_pos(2, iatom)*mass
296 com(3) = com(3) + msd%ref0_pos(3, iatom)*mass
297 mass_mol = mass_mol + mass
298 END DO ! iatom
299 msd%ref0_com_molecule(1, molecule_kind%molecule_list(imol)) = com(1)/mass_mol
300 msd%ref0_com_molecule(2, molecule_kind%molecule_list(imol)) = com(2)/mass_mol
301 msd%ref0_com_molecule(3, molecule_kind%molecule_list(imol)) = com(3)/mass_mol
302 END DO ! imol
303 END DO ! ikind
304 END IF
305
306 IF (msd%msd_region) THEN
307
308 END IF
309
310 RETURN
311998 CONTINUE ! end of file
312 cpabort("End of reference positions file reached")
313999 CONTINUE ! error
314 cpabort("Error reading reference positions file")
315
316 END SUBROUTINE initialize_msd_reftraj
317
318! **************************************************************************************************
319!> \brief ...
320!> \param reftraj ...
321!> \param md_env ...
322!> \param particle_set ...
323!> \par History
324!> 10.2007 created
325!> \author MI
326! **************************************************************************************************
327 SUBROUTINE compute_msd_reftraj(reftraj, md_env, particle_set)
328
329 TYPE(reftraj_type), POINTER :: reftraj
330 TYPE(md_environment_type), POINTER :: md_env
331 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
332
333 INTEGER :: atom, bo(2), first_atom, iatom, ikind, imol, imol_global, last_atom, mepos, &
334 natom_kind, nmol_per_kind, nmolecule, nmolkind, num_pe
335 INTEGER, DIMENSION(:), POINTER :: atom_list
336 REAL(kind=dp) :: com(3), diff2_com(4), dr2, dx, dy, dz, &
337 mass, mass_mol, msd_mkind(4), rcom(3)
338 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds
339 TYPE(atomic_kind_type), POINTER :: atomic_kind
340 TYPE(cp_subsys_type), POINTER :: subsys
341 TYPE(distribution_1d_type), POINTER :: local_molecules
342 TYPE(force_env_type), POINTER :: force_env
343 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds
344 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set
345 TYPE(molecule_kind_type), POINTER :: molecule_kind
346 TYPE(molecule_list_type), POINTER :: molecules
347 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
348 TYPE(molecule_type), POINTER :: molecule
349 TYPE(mp_para_env_type), POINTER :: para_env
350
351 NULLIFY (force_env, para_env, subsys)
352 NULLIFY (atomic_kind, atomic_kinds, atom_list)
353 NULLIFY (local_molecules, molecule, molecule_kind, molecule_kinds, &
354 molecule_kind_set, molecules, molecule_set)
355
356 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
357 CALL force_env_get(force_env=force_env, subsys=subsys)
358 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
359
360 num_pe = para_env%num_pe
361 mepos = para_env%mepos
362
363 IF (reftraj%msd%msd_kind) THEN
364 reftraj%msd%val_msd_kind = 0.0_dp
365 reftraj%msd%num_disp_atom = 0
366 reftraj%msd%disp_atom_dr = 0.0_dp
367! compute com
368 rcom = 0.0_dp
369 DO ikind = 1, atomic_kinds%n_els
370 atomic_kind => atomic_kinds%els(ikind)
371 CALL get_atomic_kind(atomic_kind=atomic_kind, &
372 atom_list=atom_list, &
373 natom=natom_kind, mass=mass)
374 bo = get_limit(natom_kind, num_pe, mepos)
375 DO iatom = bo(1), bo(2)
376 atom = atom_list(iatom)
377 rcom(1) = rcom(1) + particle_set(atom)%r(1)*mass
378 rcom(2) = rcom(2) + particle_set(atom)%r(2)*mass
379 rcom(3) = rcom(3) + particle_set(atom)%r(3)*mass
380 END DO
381 END DO
382 CALL para_env%sum(rcom)
383 rcom = rcom/reftraj%msd%total_mass
384 reftraj%msd%drcom(1) = rcom(1) - reftraj%msd%ref0_com(1)
385 reftraj%msd%drcom(2) = rcom(2) - reftraj%msd%ref0_com(2)
386 reftraj%msd%drcom(3) = rcom(3) - reftraj%msd%ref0_com(3)
387! IF(para_env%is_source()) WRITE(*,'(A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', &
388! drcom(1)*angstrom,drcom(2)*angstrom,drcom(3)*angstrom
389! compute_com
390
391 DO ikind = 1, atomic_kinds%n_els
392 atomic_kind => atomic_kinds%els(ikind)
393 CALL get_atomic_kind(atomic_kind=atomic_kind, &
394 atom_list=atom_list, &
395 natom=natom_kind)
396 bo = get_limit(natom_kind, num_pe, mepos)
397 DO iatom = bo(1), bo(2)
398 atom = atom_list(iatom)
399 dx = particle_set(atom)%r(1) - reftraj%msd%ref0_pos(1, atom) - &
400 reftraj%msd%drcom(1)
401 dy = particle_set(atom)%r(2) - reftraj%msd%ref0_pos(2, atom) - &
402 reftraj%msd%drcom(2)
403 dz = particle_set(atom)%r(3) - reftraj%msd%ref0_pos(3, atom) - &
404 reftraj%msd%drcom(3)
405 dr2 = dx*dx + dy*dy + dz*dz
406
407 reftraj%msd%val_msd_kind(1, ikind) = reftraj%msd%val_msd_kind(1, ikind) + dx*dx
408 reftraj%msd%val_msd_kind(2, ikind) = reftraj%msd%val_msd_kind(2, ikind) + dy*dy
409 reftraj%msd%val_msd_kind(3, ikind) = reftraj%msd%val_msd_kind(3, ikind) + dz*dz
410 reftraj%msd%val_msd_kind(4, ikind) = reftraj%msd%val_msd_kind(4, ikind) + dr2
411
412 IF (reftraj%msd%disp_atom) THEN
413 IF (dr2 > reftraj%msd%disp_atom_tol) THEN
414 reftraj%msd%num_disp_atom = reftraj%msd%num_disp_atom + 1
415 reftraj%msd%disp_atom_dr(1, atom) = dx
416 reftraj%msd%disp_atom_dr(2, atom) = dy
417 reftraj%msd%disp_atom_dr(3, atom) = dz
418 END IF
419 END IF
420 END DO !iatom
421 reftraj%msd%val_msd_kind(1:4, ikind) = &
422 reftraj%msd%val_msd_kind(1:4, ikind)/real(natom_kind, kind=dp)
423
424 END DO ! ikind
425 END IF
426 CALL para_env%sum(reftraj%msd%val_msd_kind)
427 CALL para_env%sum(reftraj%msd%num_disp_atom)
428 CALL para_env%sum(reftraj%msd%disp_atom_dr)
429
430 IF (reftraj%msd%msd_molecule) THEN
431 CALL cp_subsys_get(subsys=subsys, local_molecules=local_molecules, &
432 molecules=molecules, molecule_kinds=molecule_kinds)
433
434 nmolkind = molecule_kinds%n_els
435 molecule_kind_set => molecule_kinds%els
436 molecule_set => molecules%els
437
438 reftraj%msd%val_msd_molecule = 0.0_dp
439 DO ikind = 1, nmolkind
440 molecule_kind => molecule_kind_set(ikind)
441 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
442 nmol_per_kind = local_molecules%n_el(ikind)
443 msd_mkind = 0.0_dp
444 DO imol = 1, nmol_per_kind
445 imol_global = local_molecules%list(ikind)%array(imol)
446 molecule => molecule_set(imol_global)
447 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
448
449 com = 0.0_dp
450 mass_mol = 0.0_dp
451 DO iatom = first_atom, last_atom
452 mass = particle_set(iatom)%atomic_kind%mass
453 com(1) = com(1) + particle_set(iatom)%r(1)*mass
454 com(2) = com(2) + particle_set(iatom)%r(2)*mass
455 com(3) = com(3) + particle_set(iatom)%r(3)*mass
456 mass_mol = mass_mol + mass
457 END DO ! iatom
458 com(1) = com(1)/mass_mol
459 com(2) = com(2)/mass_mol
460 com(3) = com(3)/mass_mol
461 diff2_com(1) = com(1) - reftraj%msd%ref0_com_molecule(1, imol_global)
462 diff2_com(2) = com(2) - reftraj%msd%ref0_com_molecule(2, imol_global)
463 diff2_com(3) = com(3) - reftraj%msd%ref0_com_molecule(3, imol_global)
464 diff2_com(1) = diff2_com(1)*diff2_com(1)
465 diff2_com(2) = diff2_com(2)*diff2_com(2)
466 diff2_com(3) = diff2_com(3)*diff2_com(3)
467 diff2_com(4) = diff2_com(1) + diff2_com(2) + diff2_com(3)
468 msd_mkind(1) = msd_mkind(1) + diff2_com(1)
469 msd_mkind(2) = msd_mkind(2) + diff2_com(2)
470 msd_mkind(3) = msd_mkind(3) + diff2_com(3)
471 msd_mkind(4) = msd_mkind(4) + diff2_com(4)
472 END DO ! imol
473
474 reftraj%msd%val_msd_molecule(1, ikind) = msd_mkind(1)/real(nmolecule, kind=dp)
475 reftraj%msd%val_msd_molecule(2, ikind) = msd_mkind(2)/real(nmolecule, kind=dp)
476 reftraj%msd%val_msd_molecule(3, ikind) = msd_mkind(3)/real(nmolecule, kind=dp)
477 reftraj%msd%val_msd_molecule(4, ikind) = msd_mkind(4)/real(nmolecule, kind=dp)
478 END DO ! ikind
479 CALL para_env%sum(reftraj%msd%val_msd_molecule)
480
481 END IF
482
483 END SUBROUTINE compute_msd_reftraj
484
485! **************************************************************************************************
486!> \brief ...
487!> \param md_env ...
488!> \par History
489!> 10.2007 created
490!> \author MI
491! **************************************************************************************************
492 SUBROUTINE write_output_reftraj(md_env)
493 TYPE(md_environment_type), POINTER :: md_env
494
495 CHARACTER(LEN=default_string_length) :: my_act, my_mittle, my_pos
496 INTEGER :: iat, ikind, nkind, out_msd
497 LOGICAL, SAVE :: first_entry = .false.
498 TYPE(cp_logger_type), POINTER :: logger
499 TYPE(force_env_type), POINTER :: force_env
500 TYPE(reftraj_type), POINTER :: reftraj
501 TYPE(section_vals_type), POINTER :: reftraj_section, root_section
502
503 NULLIFY (logger)
504 logger => cp_get_default_logger()
505
506 NULLIFY (reftraj)
507 NULLIFY (reftraj_section, root_section)
508
509 CALL get_md_env(md_env=md_env, force_env=force_env, &
510 reftraj=reftraj)
511
512 CALL force_env_get(force_env=force_env, root_section=root_section)
513
514 reftraj_section => section_vals_get_subs_vals(root_section, &
515 "MOTION%MD%REFTRAJ")
516
517 my_pos = "APPEND"
518 my_act = "WRITE"
519
520 IF (reftraj%init .AND. (reftraj%isnap == reftraj%info%first_snapshot)) THEN
521 my_pos = "REWIND"
522 first_entry = .true.
523 END IF
524
525 IF (reftraj%info%msd) THEN
526 IF (reftraj%msd%msd_kind) THEN
527 nkind = SIZE(reftraj%msd%val_msd_kind, 2)
528 DO ikind = 1, nkind
529 my_mittle = "k"//trim(adjustl(cp_to_string(ikind)))
530 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_KIND", &
531 extension=".msd", file_position=my_pos, file_action=my_act, &
532 file_form="FORMATTED", middle_name=trim(my_mittle))
533 IF (out_msd > 0) THEN
534 WRITE (unit=out_msd, fmt="(I8, F12.3,4F20.10)") reftraj%itimes, &
535 reftraj%time*femtoseconds, &
536 reftraj%msd%val_msd_kind(1:4, ikind)*angstrom*angstrom
537 CALL m_flush(out_msd)
538 END IF
539 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
540 "PRINT%MSD_KIND")
541 END DO
542 END IF
543 IF (reftraj%msd%msd_molecule) THEN
544 nkind = SIZE(reftraj%msd%val_msd_molecule, 2)
545 DO ikind = 1, nkind
546 my_mittle = "mk"//trim(adjustl(cp_to_string(ikind)))
547 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_MOLECULE", &
548 extension=".msd", file_position=my_pos, file_action=my_act, &
549 file_form="FORMATTED", middle_name=trim(my_mittle))
550 IF (out_msd > 0) THEN
551 WRITE (unit=out_msd, fmt="(I8, F12.3,4F20.10)") reftraj%itimes, &
552 reftraj%time*femtoseconds, &
553 reftraj%msd%val_msd_molecule(1:4, ikind)*angstrom*angstrom
554 CALL m_flush(out_msd)
555 END IF
556 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
557 "PRINT%MSD_MOLECULE")
558 END DO
559 END IF
560 IF (reftraj%msd%disp_atom) THEN
561
562 IF (first_entry) my_pos = "REWIND"
563 my_mittle = "disp_at"
564 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%DISPLACED_ATOM", &
565 extension=".msd", file_position=my_pos, file_action=my_act, &
566 file_form="FORMATTED", middle_name=trim(my_mittle))
567 IF (out_msd > 0 .AND. reftraj%msd%num_disp_atom > 0) THEN
568 IF (first_entry) THEN
569 first_entry = .false.
570 END IF
571 WRITE (unit=out_msd, fmt="(A,T7,I8, A, T29, F12.3, A, T50, I10)") "# i = ", reftraj%itimes, " time (fs) = ", &
572 reftraj%time*femtoseconds, " nat = ", reftraj%msd%num_disp_atom
573 DO iat = 1, SIZE(reftraj%msd%disp_atom_dr, 2)
574 IF (abs(reftraj%msd%disp_atom_dr(1, iat)) > 0.0_dp) THEN
575 WRITE (unit=out_msd, fmt="(I8, 3F20.10)") iat, & !reftraj%msd%disp_atom_index(iat),&
576 reftraj%msd%disp_atom_dr(1, iat)*angstrom, &
577 reftraj%msd%disp_atom_dr(2, iat)*angstrom, &
578 reftraj%msd%disp_atom_dr(3, iat)*angstrom
579 END IF
580 END DO
581 END IF
582 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
583 "PRINT%DISPLACED_ATOM")
584 END IF
585 END IF ! msd
586 reftraj%init = .false.
587
588 END SUBROUTINE write_output_reftraj
589
590END MODULE reftraj_util
591
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. Tracking of preconnections.
Definition cp_files.F:16
cp_files::open_file
subroutine, public open_file(file_name, file_status, file_form, file_action, file_position, file_pad, unit_number, debug, skip_get_unit_number, file_access)
Opens the requested file using a free unit number.
Definition cp_files.F:308
cp_files::close_file
subroutine, public close_file(unit_number, file_status, keep_preconnection)
Close an open file given by its logical unit number. Optionally, keep the file and unit preconnected.
Definition cp_files.F:119
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:853
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:1063
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)
returns information about various attributes of the given subsys
Definition cp_subsys_types.F:345
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:1150
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:342
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:731
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:1047
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:130
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:195
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
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:510
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:277
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:493
reftraj_util::compute_msd_reftraj
subroutine, public compute_msd_reftraj(reftraj, md_env, particle_set)
...
Definition reftraj_util.F:328
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:53
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
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:161
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:110
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:53
reftraj_types::reftraj_type
Definition reftraj_types.F:67
simpar_types::simpar_type
Simulation parameter type for molecular dynamics.
Definition simpar_types.F:30