LCOV - code coverage report
Current view: top level - src/motion - reftraj_util.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5417965) Lines: 180 278 64.7 %
Date: 2024-02-24 07:43:42 Functions: 4 4 100.0 %

          Line data    Source code
       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          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          36 :       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          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          36 :       reftraj%isnap = nskip
     115          36 :       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          36 :       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          36 :       reftraj%natom = natom
     132          36 :       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          36 :       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          36 :    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         286 :    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         286 :       NULLIFY (logger)
     504         286 :       logger => cp_get_default_logger()
     505             : 
     506         286 :       NULLIFY (reftraj)
     507         286 :       NULLIFY (reftraj_section, root_section)
     508             : 
     509             :       CALL get_md_env(md_env=md_env, force_env=force_env, &
     510         286 :                       reftraj=reftraj)
     511             : 
     512         286 :       CALL force_env_get(force_env=force_env, root_section=root_section)
     513             : 
     514             :       reftraj_section => section_vals_get_subs_vals(root_section, &
     515         286 :                                                     "MOTION%MD%REFTRAJ")
     516             : 
     517         286 :       my_pos = "APPEND"
     518         286 :       my_act = "WRITE"
     519             : 
     520         286 :       IF (reftraj%init .AND. (reftraj%isnap == reftraj%info%first_snapshot)) THEN
     521          30 :          my_pos = "REWIND"
     522          30 :          first_entry = .TRUE.
     523             :       END IF
     524             : 
     525         286 :       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         286 :       reftraj%init = .FALSE.
     587             : 
     588         286 :    END SUBROUTINE write_output_reftraj
     589             : 
     590             : END MODULE reftraj_util
     591             : 

Generated by: LCOV version 1.15