1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief Work with symmetry
8!> \par History
9!> \author jgh
10! **************************************************************************************************
11MODULE cp_symmetry
12   USE atomic_kind_types,               ONLY: get_atomic_kind
13   USE cell_types,                      ONLY: cell_type,&
14                                              real_to_scaled
15   USE cp_log_handling,                 ONLY: cp_get_default_logger,&
16                                              cp_logger_type
17   USE cp_output_handling,              ONLY: cp_print_key_finished_output,&
18                                              cp_print_key_unit_nr
19   USE cryssym,                         ONLY: crys_sym_gen,&
20                                              csym_type,&
21                                              print_crys_symmetry,&
22                                              release_csym_type
23   USE input_section_types,             ONLY: section_vals_get_subs_vals,&
24                                              section_vals_type,&
25                                              section_vals_val_get
26   USE kinds,                           ONLY: dp
27   USE molsym,                          ONLY: molecular_symmetry,&
28                                              molsym_type,&
29                                              print_symmetry,&
30                                              release_molsym
31   USE particle_types,                  ONLY: particle_type
32   USE physcon,                         ONLY: massunit
33   USE string_utilities,                ONLY: uppercase
34#include "./base/base_uses.f90"
35
36   IMPLICIT NONE
37
38   PRIVATE
39
40   ! Global parameters (in this module)
41
42   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'cp_symmetry'
43
44   PUBLIC :: write_symmetry
45
46! **************************************************************************************************
47
48CONTAINS
49
50! **************************************************************************************************
51!> \brief Write symmetry information to output
52!> \param particle_set  Atom coordinates and types
53!> \param cell          Cell information
54!> \param input_section Input
55!> \par History
56!> \author jgh
57! **************************************************************************************************
58   SUBROUTINE write_symmetry(particle_set, cell, input_section)
59      TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
60      TYPE(cell_type), POINTER                           :: cell
61      TYPE(section_vals_type), POINTER                   :: input_section
62
63      CHARACTER(LEN=*), PARAMETER :: routineN = 'write_symmetry', routineP = moduleN//':'//routineN
64
65      CHARACTER(LEN=2), ALLOCATABLE, DIMENSION(:)        :: element
66      CHARACTER(LEN=8)                                   :: csymm, esymm
67      INTEGER                                            :: handle, i, iw, natom, plevel
68      INTEGER, ALLOCATABLE, DIMENSION(:)                 :: atype, z
69      LOGICAL                                            :: check, molecular, pall, pcoor, pinertia, &
70                                                            prmat, psymmele
71      REAL(KIND=dp)                                      :: eps_geo
72      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: weight
73      REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :)        :: coord, scoord
74      TYPE(cp_logger_type), POINTER                      :: logger
75      TYPE(csym_type)                                    :: crys_sym
76      TYPE(molsym_type), POINTER                         :: mol_sym
77      TYPE(section_vals_type), POINTER                   :: section
78
79      CALL timeset(routineN, handle)
80
81      NULLIFY (logger)
82      NULLIFY (section)
83
84      logger => cp_get_default_logger()
85      iw = cp_print_key_unit_nr(logger=logger, &
86                                basis_section=input_section, &
87                                print_key_path="PRINT%SYMMETRY", &
88                                extension=".symLog")
89
90      IF (iw > 0) THEN
91         section => section_vals_get_subs_vals(section_vals=input_section, &
92                                               subsection_name="PRINT%SYMMETRY")
93         CALL section_vals_val_get(section_vals=section, &
94                                   keyword_name="MOLECULE", l_val=molecular)
95         CALL section_vals_val_get(section_vals=section, &
96                                   keyword_name="EPS_GEO", r_val=eps_geo)
97         IF (molecular) THEN
98
99            NULLIFY (mol_sym)
100
101            natom = SIZE(particle_set)
102            ALLOCATE (coord(3, natom), z(natom), weight(natom), atype(natom), element(natom))
103
104            DO i = 1, natom
105               CALL get_atomic_kind(particle_set(i)%atomic_kind, z=z(i))
106               CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, &
107                                    kind_number=atype(i), element_symbol=element(i), mass=weight(i))
108               coord(1:3, i) = particle_set(i)%r(1:3)
109            END DO
110            weight(:) = weight(:)/massunit
111
112            CALL molecular_symmetry(mol_sym, eps_geo, coord, atype, weight)
113
114            CALL section_vals_val_get(section_vals=section, &
115                                      keyword_name="STANDARD_ORIENTATION", l_val=pcoor)
116            CALL section_vals_val_get(section_vals=section, &
117                                      keyword_name="INERTIA", l_val=pinertia)
118            CALL section_vals_val_get(section_vals=section, &
119                                      keyword_name="SYMMETRY_ELEMENTS", l_val=psymmele)
120            CALL section_vals_val_get(section_vals=section, &
121                                      keyword_name="ALL", l_val=pall)
122            plevel = 0
123            IF (pcoor) plevel = plevel + 1
124            IF (pinertia) plevel = plevel + 10
125            IF (psymmele) plevel = plevel + 100
126            IF (pall) plevel = 1111111111
127
128            CALL print_symmetry(mol_sym, coord, atype, element, z, weight, iw, plevel)
129
130            CALL section_vals_val_get(section_vals=section, &
131                                      keyword_name="CHECK_SYMMETRY", c_val=esymm)
132            CALL uppercase(esymm)
133            IF (TRIM(esymm) /= "NONE") THEN
134               csymm = mol_sym%point_group_symbol
135               CALL uppercase(csymm)
136               check = TRIM(ADJUSTL(csymm)) == TRIM(ADJUSTL(esymm))
137               IF (.NOT. check) THEN
138                  CALL cp_warn(__LOCATION__, "Symmetry check failed: "// &
139                               "Expected symmetry:"//TRIM(ADJUSTL(esymm))// &
140                               "Calculated symmetry:"//TRIM(ADJUSTL(csymm)))
141               END IF
142               CPASSERT(check)
143            END IF
144
145            DEALLOCATE (coord, z, weight, atype, element)
146
147            CALL release_molsym(mol_sym)
148
149         ELSE
150            ! Crystal symmetry
151
152            natom = SIZE(particle_set)
153            ALLOCATE (scoord(3, natom), atype(natom))
154
155            DO i = 1, natom
156               CALL get_atomic_kind(atomic_kind=particle_set(i)%atomic_kind, kind_number=atype(i))
157               CALL real_to_scaled(scoord(1:3, i), particle_set(i)%r(1:3), cell)
158            END DO
159
160            CALL crys_sym_gen(crys_sym, scoord, atype, cell%hmat, delta=eps_geo, iounit=iw)
161
162            CALL section_vals_val_get(section_vals=section, &
163                                      keyword_name="ROTATION_MATRICES", l_val=prmat)
164            CALL section_vals_val_get(section_vals=section, &
165                                      keyword_name="ALL", l_val=pall)
166            plevel = 0
167            IF (prmat) plevel = plevel + 1
168            IF (pall) plevel = 1111111111
169            crys_sym%plevel = plevel
170
171            CALL print_crys_symmetry(crys_sym)
172
173            DEALLOCATE (scoord, atype)
174
175            CALL release_csym_type(crys_sym)
176
177         END IF
178
179      END IF
180      CALL cp_print_key_finished_output(iw, logger, input_section, "PRINT%SYMMETRY")
181
182      CALL timestop(handle)
183
184   END SUBROUTINE write_symmetry
185
186! **************************************************************************************************
187
188END MODULE cp_symmetry
189