1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief GRRM interface
8!> \author JGH - 08.2019
9! **************************************************************************************************
10MODULE grrm_utils
11
12   USE kinds,                           ONLY: dp
13   USE particle_types,                  ONLY: particle_type
14   USE physcon,                         ONLY: angstrom
15#include "./base/base_uses.f90"
16
17   IMPLICIT NONE
18
19   PRIVATE
20
21   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'grrm_utils'
22
23   PUBLIC :: write_grrm
24
25! **************************************************************************************************
26
27CONTAINS
28
29! **************************************************************************************************
30!> \brief Write GRRM interface file
31!>
32!> \param iounit ...
33!> \param particles ...
34!> \param energy ...
35!> \param dipole ...
36!> \param hessian ...
37!> \param dipder ...
38!> \param polar ...
39! **************************************************************************************************
40   SUBROUTINE write_grrm(iounit, particles, energy, dipole, hessian, dipder, polar)
41
42      INTEGER, INTENT(IN)                                :: iounit
43      TYPE(particle_type), DIMENSION(:), INTENT(IN)      :: particles
44      REAL(KIND=dp), INTENT(IN)                          :: energy
45      REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL  :: dipole
46      REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
47         OPTIONAL                                        :: hessian, dipder
48      REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
49         OPTIONAL                                        :: polar
50
51      CHARACTER(LEN=*), PARAMETER :: routineN = 'write_grrm', routineP = moduleN//':'//routineN
52      REAL(KIND=dp), PARAMETER                           :: zero = 0.0_dp
53
54      INTEGER                                            :: i, j, natom, nc
55      REAL(KIND=dp), DIMENSION(5)                        :: fz
56
57      IF (iounit > 0) THEN
58         natom = SIZE(particles)
59         WRITE (iounit, "(A7)") "RESULTS"
60         WRITE (iounit, "(A18)") "CURRENT COORDINATE"
61         DO i = 1, natom
62            WRITE (iounit, "(A,3F24.12)") TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol)), &
63               particles(i)%r(1:3)*angstrom
64         END DO
65         WRITE (iounit, "(A8,3F18.12)") "ENERGY =", energy, zero, zero
66         WRITE (iounit, "(A8,3F18.12)") "       =", zero, zero, zero
67         WRITE (iounit, "(A8,F18.12)") "S**2   =", zero
68         WRITE (iounit, "(A8)") "GRADIENT"
69         DO i = 1, natom
70            WRITE (iounit, "(F17.12)") particles(i)%f(1:3)
71         END DO
72         IF (PRESENT(dipole)) THEN
73            WRITE (iounit, "(A8,3F18.12)") "DIPOLE =", dipole(1:3)
74         ELSE
75            WRITE (iounit, "(A8,3F18.12)") "DIPOLE =", zero, zero, zero
76         END IF
77         fz = zero
78         WRITE (iounit, "(A7)") "HESSIAN"
79         IF (PRESENT(hessian)) THEN
80            nc = SIZE(hessian, 1)
81            CPASSERT(nc == 3*natom)
82            DO i = 1, nc, 5
83               DO j = i, nc
84                  WRITE (iounit, "(5(F13.9,1X))") hessian(j, i:MIN(j, i + 4))
85               END DO
86            END DO
87         ELSE
88            nc = 3*natom
89            DO i = 1, nc, 5
90               DO j = i, nc
91                  WRITE (iounit, "(5(F13.9,1X))") fz(1:MIN(j - i + 1, 5))
92               END DO
93            END DO
94         END IF
95         WRITE (iounit, "(A18)") "DIPOLE DERIVATIVES"
96         IF (PRESENT(dipder)) THEN
97            DO i = 1, 3*natom
98               WRITE (iounit, "(3(F17.12,7X))") dipder(1:3, i)
99            END DO
100         ELSE
101            DO i = 1, 3*natom
102               WRITE (iounit, "(3(F17.12,7X))") zero, zero, zero
103            END DO
104         END IF
105         WRITE (iounit, "(A14)") "POLARIZABILITY"
106         IF (PRESENT(polar)) THEN
107            WRITE (iounit, "(1(F17.12,7X))") polar(1, 1)
108            WRITE (iounit, "(2(F17.12,7X))") polar(2, 1), polar(2, 2)
109            WRITE (iounit, "(3(F17.12,7X))") polar(3, 1), polar(3, 2), polar(3, 3)
110         ELSE
111            WRITE (iounit, "(1(F17.12,7X))") zero
112            WRITE (iounit, "(2(F17.12,7X))") zero, zero
113            WRITE (iounit, "(3(F17.12,7X))") zero, zero, zero
114         END IF
115      END IF
116
117   END SUBROUTINE write_grrm
118
119END MODULE grrm_utils
120