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