(git:e7e05ae)
reftraj_util.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 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 ! **************************************************************************************************
15 MODULE 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 
70 CONTAINS
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
311 998 CONTINUE ! end of file
312  cpabort("End of reference positions file reached")
313 999 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 
590 END MODULE reftraj_util
591 
uppercase
subroutine uppercase(string)
...
Definition: dumpdcd.F:1376
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:803
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:1013
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)
returns various attributes about the force environment
Definition: force_env_types.F:329
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:724
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:1040
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:106
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:3361
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