Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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 36 : 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 36 : 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 36 : simpar=simpar)
99 36 : CALL force_env_get(force_env=force_env, subsys=subsys)
100 36 : CALL cp_subsys_get(subsys=subsys, particles=particles)
101 36 : natom = particles%n_els
102 :
103 36 : my_end = .FALSE.
104 : nline_to_skip = 0
105 :
106 36 : nskip = reftraj%info%first_snapshot - 1
107 36 : CPASSERT(nskip >= 0)
108 :
109 36 : IF (nskip > 0) THEN
110 10 : nline_to_skip = (natom + 2)*nskip
111 10 : CALL parser_get_next_line(reftraj%info%traj_parser, nline_to_skip, at_end=my_end)
112 : END IF
113 :
114 36 : reftraj%isnap = nskip
115 36 : IF (my_end) THEN
116 : CALL cp_abort(__LOCATION__, &
117 : "Reached the end of the trajectory file for REFTRAJ. Number of steps skipped "// &
118 0 : "equal to the number of steps present in the file.")
119 : END IF
120 :
121 : ! Cell File
122 36 : IF (reftraj%info%variable_volume) THEN
123 6 : IF (nskip > 0) THEN
124 6 : CALL parser_get_next_line(reftraj%info%cell_parser, nskip, at_end=my_end)
125 : END IF
126 6 : IF (my_end) THEN
127 : CALL cp_abort(__LOCATION__, &
128 : "Reached the end of the cell file for REFTRAJ. Number of steps skipped "// &
129 0 : "equal to the number of steps present in the file.")
130 : END IF
131 : END IF
132 :
133 36 : reftraj%natom = natom
134 36 : IF (reftraj%info%last_snapshot > 0) THEN
135 10 : simpar%nsteps = (reftraj%info%last_snapshot - reftraj%info%first_snapshot + 1)
136 : END IF
137 :
138 36 : IF (reftraj%info%msd) THEN
139 2 : msd_section => section_vals_get_subs_vals(reftraj_section, "MSD")
140 : ! set up and printout
141 2 : CALL initialize_msd_reftraj(reftraj%msd, msd_section, reftraj, md_env)
142 : END IF
143 :
144 36 : 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 2 : 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 2 : 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 2 : 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 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
182 :
183 2 : NULLIFY (molecule, molecules, molecule_kind, molecule_kind_set, &
184 2 : molecule_kinds, molecule_set, subsys, force_env, particles, particle_set)
185 0 : CPASSERT(.NOT. ASSOCIATED(msd))
186 :
187 14 : 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 2 : CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
197 2 : CALL force_env_get(force_env=force_env, subsys=subsys)
198 2 : CALL cp_subsys_get(subsys=subsys, particles=particles)
199 2 : particle_set => particles%els
200 2 : npart = SIZE(particle_set, 1)
201 :
202 2 : msd%ref0_unit = -1
203 2 : CALL section_vals_val_get(msd_section, "REF0_FILENAME", c_val=filename)
204 2 : CALL open_file(TRIM(filename), unit_number=msd%ref0_unit)
205 :
206 6 : ALLOCATE (msd%ref0_pos(3, reftraj%natom))
207 770 : msd%ref0_pos = 0.0_dp
208 :
209 2 : IF (para_env%is_source()) THEN
210 1 : REWIND (msd%ref0_unit)
211 1 : READ (msd%ref0_unit, *, ERR=999, END=998) natom_read
212 1 : 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 0 : TRIM(ADJUSTL(cp_to_string(reftraj%natom)))
216 0 : CPABORT(errmsg)
217 : END IF
218 1 : READ (msd%ref0_unit, '(A)', ERR=999, END=998) title
219 1 : msd%total_mass = 0.0_dp
220 4 : msd%ref0_com = 0.0_dp
221 97 : DO iatom = 1, natom_read
222 96 : READ (msd%ref0_unit, *, ERR=999, END=998) element_symbol_ref0, x, y, z
223 96 : CALL uppercase(element_symbol_ref0)
224 96 : element_symbol = TRIM(particle_set(iatom)%atomic_kind%element_symbol)
225 96 : CALL uppercase(element_symbol)
226 96 : IF (element_symbol /= element_symbol_ref0) THEN
227 : errmsg = "The MSD reference configuration shows a mismatch: Check atom "// &
228 0 : TRIM(ADJUSTL(cp_to_string(iatom)))
229 0 : CPABORT(errmsg)
230 : END IF
231 96 : x = cp_unit_to_cp2k(x, "angstrom")
232 96 : y = cp_unit_to_cp2k(y, "angstrom")
233 96 : z = cp_unit_to_cp2k(z, "angstrom")
234 96 : msd%ref0_pos(1, iatom) = x
235 96 : msd%ref0_pos(2, iatom) = y
236 96 : msd%ref0_pos(3, iatom) = z
237 96 : mass = particle_set(iatom)%atomic_kind%mass
238 96 : msd%ref0_com(1) = msd%ref0_com(1) + x*mass
239 96 : msd%ref0_com(2) = msd%ref0_com(2) + y*mass
240 96 : msd%ref0_com(3) = msd%ref0_com(3) + z*mass
241 97 : msd%total_mass = msd%total_mass + mass
242 : END DO
243 4 : msd%ref0_com = msd%ref0_com/msd%total_mass
244 : END IF
245 2 : CALL close_file(unit_number=msd%ref0_unit)
246 :
247 2 : CALL para_env%bcast(msd%total_mass)
248 1538 : CALL para_env%bcast(msd%ref0_pos)
249 2 : CALL para_env%bcast(msd%ref0_com)
250 :
251 2 : CALL section_vals_val_get(msd_section, "MSD_PER_KIND", l_val=msd%msd_kind)
252 2 : CALL section_vals_val_get(msd_section, "MSD_PER_MOLKIND", l_val=msd%msd_molecule)
253 2 : CALL section_vals_val_get(msd_section, "MSD_PER_REGION", l_val=msd%msd_region)
254 :
255 2 : CALL section_vals_val_get(msd_section, "DISPLACED_ATOM", l_val=msd%disp_atom)
256 2 : IF (msd%disp_atom) THEN
257 6 : ALLOCATE (msd%disp_atom_index(npart))
258 194 : msd%disp_atom_index = 0
259 6 : ALLOCATE (msd%disp_atom_dr(3, npart))
260 770 : msd%disp_atom_dr = 0.0_dp
261 2 : msd%msd_kind = .TRUE.
262 : END IF
263 2 : CALL section_vals_val_get(msd_section, "DISPLACEMENT_TOL", r_val=tol)
264 2 : msd%disp_atom_tol = tol*tol
265 :
266 2 : IF (msd%msd_kind) THEN
267 2 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
268 2 : nkind = atomic_kinds%n_els
269 :
270 6 : ALLOCATE (msd%val_msd_kind(4, nkind))
271 22 : msd%val_msd_kind = 0.0_dp
272 : END IF
273 :
274 2 : IF (msd%msd_molecule) THEN
275 : CALL cp_subsys_get(subsys=subsys, molecules=molecules, &
276 0 : molecule_kinds=molecule_kinds)
277 0 : nmolkind = molecule_kinds%n_els
278 0 : ALLOCATE (msd%val_msd_molecule(4, nmolkind))
279 :
280 0 : molecule_kind_set => molecule_kinds%els
281 0 : molecule_set => molecules%els
282 0 : nmol = molecules%n_els
283 :
284 0 : ALLOCATE (msd%ref0_com_molecule(3, nmol))
285 :
286 0 : DO ikind = 1, nmolkind
287 0 : molecule_kind => molecule_kind_set(ikind)
288 0 : CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
289 0 : DO imol = 1, nmolecule
290 0 : molecule => molecule_set(molecule_kind%molecule_list(imol))
291 0 : CALL get_molecule(molecule=molecule, first_atom=first_atom, last_atom=last_atom)
292 0 : com = 0.0_dp
293 0 : mass_mol = 0.0_dp
294 0 : DO iatom = first_atom, last_atom
295 0 : mass = particle_set(iatom)%atomic_kind%mass
296 0 : com(1) = com(1) + msd%ref0_pos(1, iatom)*mass
297 0 : com(2) = com(2) + msd%ref0_pos(2, iatom)*mass
298 0 : com(3) = com(3) + msd%ref0_pos(3, iatom)*mass
299 0 : mass_mol = mass_mol + mass
300 : END DO ! iatom
301 0 : msd%ref0_com_molecule(1, molecule_kind%molecule_list(imol)) = com(1)/mass_mol
302 0 : msd%ref0_com_molecule(2, molecule_kind%molecule_list(imol)) = com(2)/mass_mol
303 0 : 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 2 : RETURN
313 : 998 CONTINUE ! end of file
314 0 : CPABORT("End of reference positions file reached")
315 : 999 CONTINUE ! error
316 0 : CPABORT("Error reading reference positions file")
317 :
318 4 : 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 14 : 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 14 : 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 14 : 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 14 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
350 : TYPE(molecule_type), POINTER :: molecule
351 : TYPE(mp_para_env_type), POINTER :: para_env
352 :
353 14 : NULLIFY (force_env, para_env, subsys)
354 14 : NULLIFY (atomic_kind, atomic_kinds, atom_list)
355 14 : NULLIFY (local_molecules, molecule, molecule_kind, molecule_kinds, &
356 14 : molecule_kind_set, molecules, molecule_set)
357 :
358 14 : CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env)
359 14 : CALL force_env_get(force_env=force_env, subsys=subsys)
360 14 : CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds)
361 :
362 14 : num_pe = para_env%num_pe
363 14 : mepos = para_env%mepos
364 :
365 14 : IF (reftraj%msd%msd_kind) THEN
366 154 : reftraj%msd%val_msd_kind = 0.0_dp
367 14 : reftraj%msd%num_disp_atom = 0
368 5390 : reftraj%msd%disp_atom_dr = 0.0_dp
369 : ! compute com
370 14 : rcom = 0.0_dp
371 42 : DO ikind = 1, atomic_kinds%n_els
372 28 : atomic_kind => atomic_kinds%els(ikind)
373 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
374 : atom_list=atom_list, &
375 28 : natom=natom_kind, mass=mass)
376 28 : bo = get_limit(natom_kind, num_pe, mepos)
377 742 : DO iatom = bo(1), bo(2)
378 672 : atom = atom_list(iatom)
379 672 : rcom(1) = rcom(1) + particle_set(atom)%r(1)*mass
380 672 : rcom(2) = rcom(2) + particle_set(atom)%r(2)*mass
381 700 : rcom(3) = rcom(3) + particle_set(atom)%r(3)*mass
382 : END DO
383 : END DO
384 14 : CALL para_env%sum(rcom)
385 56 : rcom = rcom/reftraj%msd%total_mass
386 14 : reftraj%msd%drcom(1) = rcom(1) - reftraj%msd%ref0_com(1)
387 14 : reftraj%msd%drcom(2) = rcom(2) - reftraj%msd%ref0_com(2)
388 14 : 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 42 : DO ikind = 1, atomic_kinds%n_els
394 28 : atomic_kind => atomic_kinds%els(ikind)
395 : CALL get_atomic_kind(atomic_kind=atomic_kind, &
396 : atom_list=atom_list, &
397 28 : natom=natom_kind)
398 28 : bo = get_limit(natom_kind, num_pe, mepos)
399 700 : DO iatom = bo(1), bo(2)
400 672 : atom = atom_list(iatom)
401 : dx = particle_set(atom)%r(1) - reftraj%msd%ref0_pos(1, atom) - &
402 672 : reftraj%msd%drcom(1)
403 : dy = particle_set(atom)%r(2) - reftraj%msd%ref0_pos(2, atom) - &
404 672 : reftraj%msd%drcom(2)
405 : dz = particle_set(atom)%r(3) - reftraj%msd%ref0_pos(3, atom) - &
406 672 : reftraj%msd%drcom(3)
407 672 : dr2 = dx*dx + dy*dy + dz*dz
408 :
409 672 : reftraj%msd%val_msd_kind(1, ikind) = reftraj%msd%val_msd_kind(1, ikind) + dx*dx
410 672 : reftraj%msd%val_msd_kind(2, ikind) = reftraj%msd%val_msd_kind(2, ikind) + dy*dy
411 672 : reftraj%msd%val_msd_kind(3, ikind) = reftraj%msd%val_msd_kind(3, ikind) + dz*dz
412 672 : reftraj%msd%val_msd_kind(4, ikind) = reftraj%msd%val_msd_kind(4, ikind) + dr2
413 :
414 700 : IF (reftraj%msd%disp_atom) THEN
415 672 : IF (dr2 > reftraj%msd%disp_atom_tol) THEN
416 0 : reftraj%msd%num_disp_atom = reftraj%msd%num_disp_atom + 1
417 0 : reftraj%msd%disp_atom_dr(1, atom) = dx
418 0 : reftraj%msd%disp_atom_dr(2, atom) = dy
419 0 : 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 182 : reftraj%msd%val_msd_kind(1:4, ikind)/REAL(natom_kind, KIND=dp)
425 :
426 : END DO ! ikind
427 : END IF
428 294 : CALL para_env%sum(reftraj%msd%val_msd_kind)
429 14 : CALL para_env%sum(reftraj%msd%num_disp_atom)
430 10766 : CALL para_env%sum(reftraj%msd%disp_atom_dr)
431 :
432 14 : IF (reftraj%msd%msd_molecule) THEN
433 : CALL cp_subsys_get(subsys=subsys, local_molecules=local_molecules, &
434 0 : molecules=molecules, molecule_kinds=molecule_kinds)
435 :
436 0 : nmolkind = molecule_kinds%n_els
437 0 : molecule_kind_set => molecule_kinds%els
438 0 : molecule_set => molecules%els
439 :
440 0 : reftraj%msd%val_msd_molecule = 0.0_dp
441 0 : DO ikind = 1, nmolkind
442 0 : molecule_kind => molecule_kind_set(ikind)
443 0 : CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule)
444 0 : nmol_per_kind = local_molecules%n_el(ikind)
445 0 : msd_mkind = 0.0_dp
446 0 : DO imol = 1, nmol_per_kind
447 0 : imol_global = local_molecules%list(ikind)%array(imol)
448 0 : molecule => molecule_set(imol_global)
449 0 : CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom)
450 :
451 0 : com = 0.0_dp
452 0 : mass_mol = 0.0_dp
453 0 : DO iatom = first_atom, last_atom
454 0 : mass = particle_set(iatom)%atomic_kind%mass
455 0 : com(1) = com(1) + particle_set(iatom)%r(1)*mass
456 0 : com(2) = com(2) + particle_set(iatom)%r(2)*mass
457 0 : com(3) = com(3) + particle_set(iatom)%r(3)*mass
458 0 : mass_mol = mass_mol + mass
459 : END DO ! iatom
460 0 : com(1) = com(1)/mass_mol
461 0 : com(2) = com(2)/mass_mol
462 0 : com(3) = com(3)/mass_mol
463 0 : diff2_com(1) = com(1) - reftraj%msd%ref0_com_molecule(1, imol_global)
464 0 : diff2_com(2) = com(2) - reftraj%msd%ref0_com_molecule(2, imol_global)
465 0 : diff2_com(3) = com(3) - reftraj%msd%ref0_com_molecule(3, imol_global)
466 0 : diff2_com(1) = diff2_com(1)*diff2_com(1)
467 0 : diff2_com(2) = diff2_com(2)*diff2_com(2)
468 0 : diff2_com(3) = diff2_com(3)*diff2_com(3)
469 0 : diff2_com(4) = diff2_com(1) + diff2_com(2) + diff2_com(3)
470 0 : msd_mkind(1) = msd_mkind(1) + diff2_com(1)
471 0 : msd_mkind(2) = msd_mkind(2) + diff2_com(2)
472 0 : msd_mkind(3) = msd_mkind(3) + diff2_com(3)
473 0 : msd_mkind(4) = msd_mkind(4) + diff2_com(4)
474 : END DO ! imol
475 :
476 0 : reftraj%msd%val_msd_molecule(1, ikind) = msd_mkind(1)/REAL(nmolecule, KIND=dp)
477 0 : reftraj%msd%val_msd_molecule(2, ikind) = msd_mkind(2)/REAL(nmolecule, KIND=dp)
478 0 : reftraj%msd%val_msd_molecule(3, ikind) = msd_mkind(3)/REAL(nmolecule, KIND=dp)
479 0 : reftraj%msd%val_msd_molecule(4, ikind) = msd_mkind(4)/REAL(nmolecule, KIND=dp)
480 : END DO ! ikind
481 0 : CALL para_env%sum(reftraj%msd%val_msd_molecule)
482 :
483 : END IF
484 :
485 14 : 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 282 : 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 282 : NULLIFY (logger)
506 282 : logger => cp_get_default_logger()
507 :
508 282 : NULLIFY (reftraj)
509 282 : NULLIFY (reftraj_section, root_section)
510 :
511 : CALL get_md_env(md_env=md_env, force_env=force_env, &
512 282 : reftraj=reftraj)
513 :
514 282 : CALL force_env_get(force_env=force_env, root_section=root_section)
515 :
516 : reftraj_section => section_vals_get_subs_vals(root_section, &
517 282 : "MOTION%MD%REFTRAJ")
518 :
519 282 : my_pos = "APPEND"
520 282 : my_act = "WRITE"
521 :
522 282 : IF (reftraj%init .AND. (reftraj%isnap == reftraj%info%first_snapshot)) THEN
523 32 : my_pos = "REWIND"
524 32 : first_entry = .TRUE.
525 : END IF
526 :
527 282 : IF (reftraj%info%msd) THEN
528 14 : IF (reftraj%msd%msd_kind) THEN
529 14 : nkind = SIZE(reftraj%msd%val_msd_kind, 2)
530 42 : DO ikind = 1, nkind
531 28 : 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 28 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
535 28 : IF (out_msd > 0) THEN
536 14 : WRITE (UNIT=out_msd, FMT="(I8, F12.3,4F20.10)") reftraj%itimes, &
537 14 : reftraj%time*femtoseconds, &
538 84 : reftraj%msd%val_msd_kind(1:4, ikind)*angstrom*angstrom
539 14 : CALL m_flush(out_msd)
540 : END IF
541 : CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
542 42 : "PRINT%MSD_KIND")
543 : END DO
544 : END IF
545 14 : IF (reftraj%msd%msd_molecule) THEN
546 0 : nkind = SIZE(reftraj%msd%val_msd_molecule, 2)
547 0 : DO ikind = 1, nkind
548 0 : 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 0 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
552 0 : IF (out_msd > 0) THEN
553 0 : WRITE (UNIT=out_msd, FMT="(I8, F12.3,4F20.10)") reftraj%itimes, &
554 0 : reftraj%time*femtoseconds, &
555 0 : reftraj%msd%val_msd_molecule(1:4, ikind)*angstrom*angstrom
556 0 : CALL m_flush(out_msd)
557 : END IF
558 : CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, &
559 0 : "PRINT%MSD_MOLECULE")
560 : END DO
561 : END IF
562 14 : IF (reftraj%msd%disp_atom) THEN
563 :
564 14 : IF (first_entry) my_pos = "REWIND"
565 14 : 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 14 : file_form="FORMATTED", middle_name=TRIM(my_mittle))
569 14 : IF (out_msd > 0 .AND. reftraj%msd%num_disp_atom > 0) THEN
570 0 : IF (first_entry) THEN
571 0 : first_entry = .FALSE.
572 : END IF
573 0 : WRITE (UNIT=out_msd, FMT="(A,T7,I8, A, T29, F12.3, A, T50, I10)") "# i = ", reftraj%itimes, " time (fs) = ", &
574 0 : reftraj%time*femtoseconds, " nat = ", reftraj%msd%num_disp_atom
575 0 : DO iat = 1, SIZE(reftraj%msd%disp_atom_dr, 2)
576 0 : IF (ABS(reftraj%msd%disp_atom_dr(1, iat)) > 0.0_dp) THEN
577 0 : WRITE (UNIT=out_msd, FMT="(I8, 3F20.10)") iat, & !reftraj%msd%disp_atom_index(iat),&
578 0 : reftraj%msd%disp_atom_dr(1, iat)*angstrom, &
579 0 : reftraj%msd%disp_atom_dr(2, iat)*angstrom, &
580 0 : 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 14 : "PRINT%DISPLACED_ATOM")
586 : END IF
587 : END IF ! msd
588 282 : reftraj%init = .FALSE.
589 :
590 282 : END SUBROUTINE write_output_reftraj
591 :
592 : END MODULE reftraj_util
593 :
|
