1!--------------------------------------------------------------------------------------------------! 2! CP2K: A general program to perform molecular dynamics simulations ! 3! Copyright (C) 2000 - 2019 CP2K developers group ! 4!--------------------------------------------------------------------------------------------------! 5 6! ************************************************************************************************** 7!> \brief Initialize the analysis of trajectories to be done 8!> by activating the REFTRAJ ensemble 9!> \par History 10!> Created 10-07 [MI] 11!> \author MI 12! ************************************************************************************************** 13MODULE reftraj_util 14 15 USE atomic_kind_list_types, ONLY: atomic_kind_list_type 16 USE atomic_kind_types, ONLY: atomic_kind_type,& 17 get_atomic_kind 18 USE cp_files, ONLY: close_file,& 19 open_file 20 USE cp_log_handling, ONLY: cp_get_default_logger,& 21 cp_logger_type,& 22 cp_to_string 23 USE cp_output_handling, ONLY: cp_print_key_finished_output,& 24 cp_print_key_unit_nr 25 USE cp_para_types, ONLY: cp_para_env_type 26 USE cp_parser_methods, ONLY: parser_get_next_line 27 USE cp_subsys_types, ONLY: cp_subsys_get,& 28 cp_subsys_type 29 USE cp_units, ONLY: cp_unit_to_cp2k 30 USE distribution_1d_types, ONLY: distribution_1d_type 31 USE force_env_types, ONLY: force_env_get,& 32 force_env_type 33 USE input_section_types, ONLY: section_vals_get_subs_vals,& 34 section_vals_type,& 35 section_vals_val_get 36 USE kinds, ONLY: default_path_length,& 37 default_string_length,& 38 dp 39 USE machine, ONLY: m_flush 40 USE md_environment_types, ONLY: get_md_env,& 41 md_environment_type 42 USE message_passing, ONLY: mp_bcast,& 43 mp_sum 44 USE molecule_kind_list_types, ONLY: molecule_kind_list_type 45 USE molecule_kind_types, ONLY: get_molecule_kind,& 46 molecule_kind_type 47 USE molecule_list_types, ONLY: molecule_list_type 48 USE molecule_types, ONLY: get_molecule,& 49 molecule_type 50 USE particle_list_types, ONLY: particle_list_type 51 USE particle_types, ONLY: particle_type 52 USE physcon, ONLY: angstrom,& 53 femtoseconds 54 USE reftraj_types, ONLY: reftraj_msd_type,& 55 reftraj_type 56 USE simpar_types, ONLY: simpar_type 57 USE util, ONLY: get_limit 58#include "../base/base_uses.f90" 59 60 IMPLICIT NONE 61 62 PRIVATE 63 64 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'reftraj_util' 65 66 PUBLIC :: initialize_reftraj, compute_msd_reftraj, write_output_reftraj 67 68CONTAINS 69 70! ************************************************************************************************** 71!> \brief ... 72!> \param reftraj ... 73!> \param reftraj_section ... 74!> \param md_env ... 75!> \par History 76!> 10.2007 created 77!> \author MI 78! ************************************************************************************************** 79 SUBROUTINE initialize_reftraj(reftraj, reftraj_section, md_env) 80 81 TYPE(reftraj_type), POINTER :: reftraj 82 TYPE(section_vals_type), POINTER :: reftraj_section 83 TYPE(md_environment_type), POINTER :: md_env 84 85 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_reftraj', & 86 routineP = moduleN//':'//routineN 87 88 INTEGER :: natom, nline_to_skip, nskip 89 LOGICAL :: my_end 90 TYPE(cp_para_env_type), POINTER :: para_env 91 TYPE(cp_subsys_type), POINTER :: subsys 92 TYPE(force_env_type), POINTER :: force_env 93 TYPE(particle_list_type), POINTER :: particles 94 TYPE(section_vals_type), POINTER :: msd_section 95 TYPE(simpar_type), POINTER :: simpar 96 97 NULLIFY (force_env, msd_section, particles, simpar, subsys) 98 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env, & 99 simpar=simpar) 100 CALL force_env_get(force_env=force_env, subsys=subsys) 101 CALL cp_subsys_get(subsys=subsys, particles=particles) 102 natom = particles%n_els 103 104 my_end = .FALSE. 105 nline_to_skip = 0 106 107 nskip = reftraj%info%first_snapshot - 1 108 CPASSERT(nskip >= 0) 109 110 IF (nskip > 0) THEN 111 nline_to_skip = (natom + 2)*nskip 112 CALL parser_get_next_line(reftraj%info%traj_parser, nline_to_skip, at_end=my_end) 113 END IF 114 115 reftraj%isnap = nskip 116 IF (my_end) & 117 CALL cp_abort(__LOCATION__, & 118 "Reached the end of the trajectory file for REFTRAJ. Number of steps skipped "// & 119 "equal to the number of steps present in the file.") 120 121 ! Cell File 122 IF (reftraj%info%variable_volume) THEN 123 IF (nskip > 0) THEN 124 CALL parser_get_next_line(reftraj%info%cell_parser, nskip, at_end=my_end) 125 END IF 126 IF (my_end) & 127 CALL cp_abort(__LOCATION__, & 128 "Reached the end of the cell file for REFTRAJ. Number of steps skipped "// & 129 "equal to the number of steps present in the file.") 130 END IF 131 132 reftraj%natom = natom 133 IF (reftraj%info%last_snapshot > 0) THEN 134 simpar%nsteps = (reftraj%info%last_snapshot - reftraj%info%first_snapshot + 1) 135 END IF 136 137 IF (reftraj%info%msd) THEN 138 msd_section => section_vals_get_subs_vals(reftraj_section, "MSD") 139 ! set up and printout 140 CALL initialize_msd_reftraj(reftraj%msd, msd_section, reftraj, md_env) 141 END IF 142 143 END SUBROUTINE initialize_reftraj 144 145! ************************************************************************************************** 146!> \brief ... 147!> \param msd ... 148!> \param msd_section ... 149!> \param reftraj ... 150!> \param md_env ... 151!> \par History 152!> 10.2007 created 153!> \author MI 154! ************************************************************************************************** 155 SUBROUTINE initialize_msd_reftraj(msd, msd_section, reftraj, md_env) 156 TYPE(reftraj_msd_type), POINTER :: msd 157 TYPE(section_vals_type), POINTER :: msd_section 158 TYPE(reftraj_type), POINTER :: reftraj 159 TYPE(md_environment_type), POINTER :: md_env 160 161 CHARACTER(len=*), PARAMETER :: routineN = 'initialize_msd_reftraj', & 162 routineP = moduleN//':'//routineN 163 164 CHARACTER(LEN=10) :: AA 165 CHARACTER(LEN=default_path_length) :: filename 166 CHARACTER(LEN=default_string_length) :: name, title 167 INTEGER :: first_atom, iatom, ikind, imol, & 168 last_atom, natom_read, nkind, nmol, & 169 nmolecule, nmolkind, npart 170 REAL(KIND=dp) :: com(3), mass, mass_mol, tol, x, y, z 171 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 172 TYPE(cp_para_env_type), POINTER :: para_env 173 TYPE(cp_subsys_type), POINTER :: subsys 174 TYPE(force_env_type), POINTER :: force_env 175 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 176 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 177 TYPE(molecule_kind_type), POINTER :: molecule_kind 178 TYPE(molecule_list_type), POINTER :: molecules 179 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 180 TYPE(molecule_type), POINTER :: molecule 181 TYPE(particle_list_type), POINTER :: particles 182 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 183 184 NULLIFY (molecule, molecules, molecule_kind, molecule_kind_set, & 185 molecule_kinds, molecule_set, subsys, force_env, particles, particle_set) 186 CPASSERT(.NOT. ASSOCIATED(msd)) 187 188 ALLOCATE (msd) 189 190 NULLIFY (msd%ref0_pos) 191 NULLIFY (msd%ref0_com_molecule) 192 NULLIFY (msd%val_msd_kind) 193 NULLIFY (msd%val_msd_molecule) 194 NULLIFY (msd%disp_atom_index) 195 NULLIFY (msd%disp_atom_dr) 196 197 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env) 198 CALL force_env_get(force_env=force_env, subsys=subsys) 199 CALL cp_subsys_get(subsys=subsys, particles=particles) 200 particle_set => particles%els 201 npart = SIZE(particle_set, 1) 202 203 msd%ref0_unit = -1 204 CALL section_vals_val_get(msd_section, "REF0_FILENAME", c_val=filename) 205 CALL open_file(TRIM(filename), unit_number=msd%ref0_unit) 206 207 ALLOCATE (msd%ref0_pos(3, reftraj%natom)) 208 msd%ref0_pos = 0.0_dp 209 210 IF (para_env%mepos == para_env%source) THEN 211 REWIND (msd%ref0_unit) 212 READ (msd%ref0_unit, *, ERR=999, END=998) natom_read 213 CPASSERT(natom_read == reftraj%natom) 214 READ (msd%ref0_unit, '(A)', ERR=999, END=998) title 215 msd%total_mass = 0.0_dp 216 msd%ref0_com = 0.0_dp 217 DO iatom = 1, natom_read 218 READ (msd%ref0_unit, *, ERR=999, END=998) AA, x, y, z 219 name = TRIM(particle_set(iatom)%atomic_kind%element_symbol) 220 CPASSERT((TRIM(AA) == name)) 221 222 x = cp_unit_to_cp2k(x, "angstrom") 223 y = cp_unit_to_cp2k(y, "angstrom") 224 z = cp_unit_to_cp2k(z, "angstrom") 225 msd%ref0_pos(1, iatom) = x 226 msd%ref0_pos(2, iatom) = y 227 msd%ref0_pos(3, iatom) = z 228 mass = particle_set(iatom)%atomic_kind%mass 229 msd%ref0_com(1) = msd%ref0_com(1) + x*mass 230 msd%ref0_com(2) = msd%ref0_com(2) + y*mass 231 msd%ref0_com(3) = msd%ref0_com(3) + z*mass 232 msd%total_mass = msd%total_mass + mass 233 END DO 234 msd%ref0_com = msd%ref0_com/msd%total_mass 235 END IF 236 CALL close_file(unit_number=msd%ref0_unit) 237 238 CALL mp_bcast(msd%total_mass, para_env%source, para_env%group) 239 CALL mp_bcast(msd%ref0_pos, para_env%source, para_env%group) 240 CALL mp_bcast(msd%ref0_com, para_env%source, para_env%group) 241 242 CALL section_vals_val_get(msd_section, "MSD_PER_KIND", l_val=msd%msd_kind) 243 CALL section_vals_val_get(msd_section, "MSD_PER_MOLKIND", l_val=msd%msd_molecule) 244 CALL section_vals_val_get(msd_section, "MSD_PER_REGION", l_val=msd%msd_region) 245 246 CALL section_vals_val_get(msd_section, "DISPLACED_ATOM", l_val=msd%disp_atom) 247 IF (msd%disp_atom) THEN 248 ALLOCATE (msd%disp_atom_index(npart)) 249 msd%disp_atom_index = 0 250 ALLOCATE (msd%disp_atom_dr(3, npart)) 251 msd%disp_atom_dr = 0.0_dp 252 msd%msd_kind = .TRUE. 253 END IF 254 CALL section_vals_val_get(msd_section, "DISPLACEMENT_TOL", r_val=tol) 255 msd%disp_atom_tol = tol*tol 256 257 IF (msd%msd_kind) THEN 258 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds) 259 nkind = atomic_kinds%n_els 260 261 ALLOCATE (msd%val_msd_kind(4, nkind)) 262 msd%val_msd_kind = 0.0_dp 263 END IF 264 265 IF (msd%msd_molecule) THEN 266 CALL cp_subsys_get(subsys=subsys, molecules=molecules, & 267 molecule_kinds=molecule_kinds) 268 nmolkind = molecule_kinds%n_els 269 ALLOCATE (msd%val_msd_molecule(4, nmolkind)) 270 271 molecule_kind_set => molecule_kinds%els 272 molecule_set => molecules%els 273 nmol = molecules%n_els 274 275 ALLOCATE (msd%ref0_com_molecule(3, nmol)) 276 277 DO ikind = 1, nmolkind 278 molecule_kind => molecule_kind_set(ikind) 279 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule) 280 DO imol = 1, nmolecule 281 molecule => molecule_set(molecule_kind%molecule_list(imol)) 282 CALL get_molecule(molecule=molecule, first_atom=first_atom, last_atom=last_atom) 283 com = 0.0_dp 284 mass_mol = 0.0_dp 285 DO iatom = first_atom, last_atom 286 mass = particle_set(iatom)%atomic_kind%mass 287 com(1) = com(1) + msd%ref0_pos(1, iatom)*mass 288 com(2) = com(2) + msd%ref0_pos(2, iatom)*mass 289 com(3) = com(3) + msd%ref0_pos(3, iatom)*mass 290 mass_mol = mass_mol + mass 291 ENDDO ! iatom 292 msd%ref0_com_molecule(1, molecule_kind%molecule_list(imol)) = com(1)/mass_mol 293 msd%ref0_com_molecule(2, molecule_kind%molecule_list(imol)) = com(2)/mass_mol 294 msd%ref0_com_molecule(3, molecule_kind%molecule_list(imol)) = com(3)/mass_mol 295 END DO ! imol 296 ENDDO ! ikind 297 END IF 298 299 IF (msd%msd_region) THEN 300 301 END IF 302 303 RETURN 304998 CONTINUE ! end of file 305 CPABORT("End of reference positions file reached") 306999 CONTINUE ! error 307 CPABORT("Error reading reference positions file") 308 309 END SUBROUTINE initialize_msd_reftraj 310 311! ************************************************************************************************** 312!> \brief ... 313!> \param reftraj ... 314!> \param md_env ... 315!> \param particle_set ... 316!> \par History 317!> 10.2007 created 318!> \author MI 319! ************************************************************************************************** 320 SUBROUTINE compute_msd_reftraj(reftraj, md_env, particle_set) 321 322 TYPE(reftraj_type), POINTER :: reftraj 323 TYPE(md_environment_type), POINTER :: md_env 324 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set 325 326 CHARACTER(len=*), PARAMETER :: routineN = 'compute_msd_reftraj', & 327 routineP = moduleN//':'//routineN 328 329 INTEGER :: atom, bo(2), first_atom, iatom, ikind, imol, imol_global, last_atom, mepos, & 330 natom_kind, nmol_per_kind, nmolecule, nmolkind, num_pe 331 INTEGER, DIMENSION(:), POINTER :: atom_list 332 REAL(KIND=dp) :: com(3), diff2_com(4), dr2, dx, dy, dz, & 333 mass, mass_mol, msd_mkind(4), rcom(3) 334 TYPE(atomic_kind_list_type), POINTER :: atomic_kinds 335 TYPE(atomic_kind_type), POINTER :: atomic_kind 336 TYPE(cp_para_env_type), POINTER :: para_env 337 TYPE(cp_subsys_type), POINTER :: subsys 338 TYPE(distribution_1d_type), POINTER :: local_molecules 339 TYPE(force_env_type), POINTER :: force_env 340 TYPE(molecule_kind_list_type), POINTER :: molecule_kinds 341 TYPE(molecule_kind_type), DIMENSION(:), POINTER :: molecule_kind_set 342 TYPE(molecule_kind_type), POINTER :: molecule_kind 343 TYPE(molecule_list_type), POINTER :: molecules 344 TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set 345 TYPE(molecule_type), POINTER :: molecule 346 347 NULLIFY (force_env, para_env, subsys) 348 NULLIFY (atomic_kind, atomic_kinds, atom_list) 349 NULLIFY (local_molecules, molecule, molecule_kind, molecule_kinds, & 350 molecule_kind_set, molecules, molecule_set) 351 352 CALL get_md_env(md_env=md_env, force_env=force_env, para_env=para_env) 353 CALL force_env_get(force_env=force_env, subsys=subsys) 354 CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds) 355 356 num_pe = para_env%num_pe 357 mepos = para_env%mepos 358 359 IF (reftraj%msd%msd_kind) THEN 360 reftraj%msd%val_msd_kind = 0.0_dp 361 reftraj%msd%num_disp_atom = 0 362 reftraj%msd%disp_atom_dr = 0.0_dp 363! compute com 364 rcom = 0.0_dp 365 DO ikind = 1, atomic_kinds%n_els 366 atomic_kind => atomic_kinds%els(ikind) 367 CALL get_atomic_kind(atomic_kind=atomic_kind, & 368 atom_list=atom_list, & 369 natom=natom_kind, mass=mass) 370 bo = get_limit(natom_kind, num_pe, mepos) 371 DO iatom = bo(1), bo(2) 372 atom = atom_list(iatom) 373 rcom(1) = rcom(1) + particle_set(atom)%r(1)*mass 374 rcom(2) = rcom(2) + particle_set(atom)%r(2)*mass 375 rcom(3) = rcom(3) + particle_set(atom)%r(3)*mass 376 END DO 377 END DO 378 CALL mp_sum(rcom, para_env%group) 379 rcom = rcom/reftraj%msd%total_mass 380 reftraj%msd%drcom(1) = rcom(1) - reftraj%msd%ref0_com(1) 381 reftraj%msd%drcom(2) = rcom(2) - reftraj%msd%ref0_com(2) 382 reftraj%msd%drcom(3) = rcom(3) - reftraj%msd%ref0_com(3) 383! IF(para_env%ionode) WRITE(*,'(A,T50,3f10.5)') ' COM displacement (dx,dy,dz) [angstrom]: ', & 384! drcom(1)*angstrom,drcom(2)*angstrom,drcom(3)*angstrom 385! compute_com 386 387 DO ikind = 1, atomic_kinds%n_els 388 atomic_kind => atomic_kinds%els(ikind) 389 CALL get_atomic_kind(atomic_kind=atomic_kind, & 390 atom_list=atom_list, & 391 natom=natom_kind) 392 bo = get_limit(natom_kind, num_pe, mepos) 393 DO iatom = bo(1), bo(2) 394 atom = atom_list(iatom) 395 dx = particle_set(atom)%r(1) - reftraj%msd%ref0_pos(1, atom) - & 396 reftraj%msd%drcom(1) 397 dy = particle_set(atom)%r(2) - reftraj%msd%ref0_pos(2, atom) - & 398 reftraj%msd%drcom(2) 399 dz = particle_set(atom)%r(3) - reftraj%msd%ref0_pos(3, atom) - & 400 reftraj%msd%drcom(3) 401 dr2 = dx*dx + dy*dy + dz*dz 402 403 reftraj%msd%val_msd_kind(1, ikind) = reftraj%msd%val_msd_kind(1, ikind) + dx*dx 404 reftraj%msd%val_msd_kind(2, ikind) = reftraj%msd%val_msd_kind(2, ikind) + dy*dy 405 reftraj%msd%val_msd_kind(3, ikind) = reftraj%msd%val_msd_kind(3, ikind) + dz*dz 406 reftraj%msd%val_msd_kind(4, ikind) = reftraj%msd%val_msd_kind(4, ikind) + dr2 407 408 IF (reftraj%msd%disp_atom) THEN 409 IF (dr2 > reftraj%msd%disp_atom_tol) THEN 410 reftraj%msd%num_disp_atom = reftraj%msd%num_disp_atom + 1 411 reftraj%msd%disp_atom_dr(1, atom) = dx 412 reftraj%msd%disp_atom_dr(2, atom) = dy 413 reftraj%msd%disp_atom_dr(3, atom) = dz 414 END IF 415 END IF 416 END DO !iatom 417 reftraj%msd%val_msd_kind(1:4, ikind) = & 418 reftraj%msd%val_msd_kind(1:4, ikind)/REAL(natom_kind, KIND=dp) 419 420 END DO ! ikind 421 ENDIF 422 CALL mp_sum(reftraj%msd%val_msd_kind, para_env%group) 423 CALL mp_sum(reftraj%msd%num_disp_atom, para_env%group) 424 CALL mp_sum(reftraj%msd%disp_atom_dr, para_env%group) 425 426 IF (reftraj%msd%msd_molecule) THEN 427 CALL cp_subsys_get(subsys=subsys, local_molecules=local_molecules, & 428 molecules=molecules, molecule_kinds=molecule_kinds) 429 430 nmolkind = molecule_kinds%n_els 431 molecule_kind_set => molecule_kinds%els 432 molecule_set => molecules%els 433 434 reftraj%msd%val_msd_molecule = 0.0_dp 435 DO ikind = 1, nmolkind 436 molecule_kind => molecule_kind_set(ikind) 437 CALL get_molecule_kind(molecule_kind=molecule_kind, nmolecule=nmolecule) 438 nmol_per_kind = local_molecules%n_el(ikind) 439 msd_mkind = 0.0_dp 440 DO imol = 1, nmol_per_kind 441 imol_global = local_molecules%list(ikind)%array(imol) 442 molecule => molecule_set(imol_global) 443 CALL get_molecule(molecule, first_atom=first_atom, last_atom=last_atom) 444 445 com = 0.0_dp 446 mass_mol = 0.0_dp 447 DO iatom = first_atom, last_atom 448 mass = particle_set(iatom)%atomic_kind%mass 449 com(1) = com(1) + particle_set(iatom)%r(1)*mass 450 com(2) = com(2) + particle_set(iatom)%r(2)*mass 451 com(3) = com(3) + particle_set(iatom)%r(3)*mass 452 mass_mol = mass_mol + mass 453 ENDDO ! iatom 454 com(1) = com(1)/mass_mol 455 com(2) = com(2)/mass_mol 456 com(3) = com(3)/mass_mol 457 diff2_com(1) = com(1) - reftraj%msd%ref0_com_molecule(1, imol_global) 458 diff2_com(2) = com(2) - reftraj%msd%ref0_com_molecule(2, imol_global) 459 diff2_com(3) = com(3) - reftraj%msd%ref0_com_molecule(3, imol_global) 460 diff2_com(1) = diff2_com(1)*diff2_com(1) 461 diff2_com(2) = diff2_com(2)*diff2_com(2) 462 diff2_com(3) = diff2_com(3)*diff2_com(3) 463 diff2_com(4) = diff2_com(1) + diff2_com(2) + diff2_com(3) 464 msd_mkind(1) = msd_mkind(1) + diff2_com(1) 465 msd_mkind(2) = msd_mkind(2) + diff2_com(2) 466 msd_mkind(3) = msd_mkind(3) + diff2_com(3) 467 msd_mkind(4) = msd_mkind(4) + diff2_com(4) 468 ENDDO ! imol 469 470 reftraj%msd%val_msd_molecule(1, ikind) = msd_mkind(1)/REAL(nmolecule, KIND=dp) 471 reftraj%msd%val_msd_molecule(2, ikind) = msd_mkind(2)/REAL(nmolecule, KIND=dp) 472 reftraj%msd%val_msd_molecule(3, ikind) = msd_mkind(3)/REAL(nmolecule, KIND=dp) 473 reftraj%msd%val_msd_molecule(4, ikind) = msd_mkind(4)/REAL(nmolecule, KIND=dp) 474 END DO ! ikind 475 CALL mp_sum(reftraj%msd%val_msd_molecule, para_env%group) 476 477 END IF 478 479 END SUBROUTINE compute_msd_reftraj 480 481! ************************************************************************************************** 482!> \brief ... 483!> \param md_env ... 484!> \par History 485!> 10.2007 created 486!> \author MI 487! ************************************************************************************************** 488 SUBROUTINE write_output_reftraj(md_env) 489 TYPE(md_environment_type), POINTER :: md_env 490 491 CHARACTER(len=*), PARAMETER :: routineN = 'write_output_reftraj', & 492 routineP = moduleN//':'//routineN 493 494 CHARACTER(LEN=default_string_length) :: my_act, my_mittle, my_pos 495 INTEGER :: iat, ikind, nkind, out_msd 496 LOGICAL, SAVE :: first_entry = .FALSE. 497 TYPE(cp_logger_type), POINTER :: logger 498 TYPE(force_env_type), POINTER :: force_env 499 TYPE(reftraj_type), POINTER :: reftraj 500 TYPE(section_vals_type), POINTER :: reftraj_section, root_section 501 502 NULLIFY (logger) 503 logger => cp_get_default_logger() 504 505 NULLIFY (reftraj) 506 NULLIFY (reftraj_section, root_section) 507 508 CALL get_md_env(md_env=md_env, force_env=force_env, & 509 reftraj=reftraj) 510 511 CALL force_env_get(force_env=force_env, root_section=root_section) 512 513 reftraj_section => section_vals_get_subs_vals(root_section, & 514 "MOTION%MD%REFTRAJ") 515 516 my_pos = "APPEND" 517 my_act = "WRITE" 518 519 IF (reftraj%init .AND. (reftraj%isnap == reftraj%info%first_snapshot)) THEN 520 my_pos = "REWIND" 521 first_entry = .TRUE. 522 END IF 523 524 IF (reftraj%info%msd) THEN 525 IF (reftraj%msd%msd_kind) THEN 526 nkind = SIZE(reftraj%msd%val_msd_kind, 2) 527 DO ikind = 1, nkind 528 my_mittle = "k"//TRIM(ADJUSTL(cp_to_string(ikind))) 529 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_KIND", & 530 extension=".msd", file_position=my_pos, file_action=my_act, & 531 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 532 IF (out_msd > 0) THEN 533 WRITE (UNIT=out_msd, FMT="(I8, F12.3,4F20.10)") reftraj%itimes, & 534 reftraj%time*femtoseconds, & 535 reftraj%msd%val_msd_kind(1:4, ikind)*angstrom*angstrom 536 CALL m_flush(out_msd) 537 END IF 538 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, & 539 "PRINT%MSD_KIND") 540 END DO 541 END IF 542 IF (reftraj%msd%msd_molecule) THEN 543 nkind = SIZE(reftraj%msd%val_msd_molecule, 2) 544 DO ikind = 1, nkind 545 my_mittle = "mk"//TRIM(ADJUSTL(cp_to_string(ikind))) 546 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%MSD_MOLECULE", & 547 extension=".msd", file_position=my_pos, file_action=my_act, & 548 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 549 IF (out_msd > 0) THEN 550 WRITE (UNIT=out_msd, FMT="(I8, F12.3,4F20.10)") reftraj%itimes, & 551 reftraj%time*femtoseconds, & 552 reftraj%msd%val_msd_molecule(1:4, ikind)*angstrom*angstrom 553 CALL m_flush(out_msd) 554 END IF 555 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, & 556 "PRINT%MSD_MOLECULE") 557 END DO 558 END IF 559 IF (reftraj%msd%disp_atom) THEN 560 561 IF (first_entry) my_pos = "REWIND" 562 my_mittle = "disp_at" 563 out_msd = cp_print_key_unit_nr(logger, reftraj_section, "PRINT%DISPLACED_ATOM", & 564 extension=".msd", file_position=my_pos, file_action=my_act, & 565 file_form="FORMATTED", middle_name=TRIM(my_mittle)) 566 IF (out_msd > 0 .AND. reftraj%msd%num_disp_atom > 0) THEN 567 IF (first_entry) THEN 568 first_entry = .FALSE. 569 END IF 570 WRITE (UNIT=out_msd, FMT="(A,T7,I8, A, T29, F12.3, A, T50, I10)") "# i = ", reftraj%itimes, " time (fs) = ", & 571 reftraj%time*femtoseconds, " nat = ", reftraj%msd%num_disp_atom 572 DO iat = 1, SIZE(reftraj%msd%disp_atom_dr, 2) 573 IF (ABS(reftraj%msd%disp_atom_dr(1, iat)) > 0.0_dp) THEN 574 WRITE (UNIT=out_msd, FMT="(I8, 3F20.10)") iat, & !reftraj%msd%disp_atom_index(iat),& 575 reftraj%msd%disp_atom_dr(1, iat)*angstrom, & 576 reftraj%msd%disp_atom_dr(2, iat)*angstrom, & 577 reftraj%msd%disp_atom_dr(3, iat)*angstrom 578 END IF 579 END DO 580 ENDIF 581 CALL cp_print_key_finished_output(out_msd, logger, reftraj_section, & 582 "PRINT%DISPLACED_ATOM") 583 END IF 584 ENDIF ! msd 585 reftraj%init = .FALSE. 586 587 END SUBROUTINE write_output_reftraj 588 589END MODULE reftraj_util 590 591