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 GRRM interface
10 : !> \author JGH - 08.2019
11 : ! **************************************************************************************************
12 : MODULE grrm_utils
13 :
14 : USE cp_control_types, ONLY: dft_control_type
15 : USE force_env_types, ONLY: force_env_type
16 : USE kinds, ONLY: dp
17 : USE particle_types, ONLY: particle_type
18 : USE physcon, ONLY: angstrom
19 : USE qs_energy_types, ONLY: qs_energy_type
20 : USE qs_environment_types, ONLY: get_qs_env
21 : #include "./base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 :
25 : PRIVATE
26 :
27 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'grrm_utils'
28 :
29 : PUBLIC :: write_grrm
30 :
31 : ! **************************************************************************************************
32 :
33 : CONTAINS
34 :
35 : ! **************************************************************************************************
36 : !> \brief Write GRRM interface file
37 : !>
38 : !> \param iounit ...
39 : !> \param force_env ...
40 : !> \param particles ...
41 : !> \param energy ...
42 : !> \param dipole ...
43 : !> \param hessian ...
44 : !> \param dipder ...
45 : !> \param polar ...
46 : !> \param fixed_atoms ...
47 : ! **************************************************************************************************
48 40 : SUBROUTINE write_grrm(iounit, force_env, particles, energy, dipole, hessian, dipder, polar, &
49 : fixed_atoms)
50 :
51 : INTEGER, INTENT(IN) :: iounit
52 : TYPE(force_env_type), POINTER :: force_env
53 : TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particles
54 : REAL(KIND=dp), INTENT(IN) :: energy
55 : REAL(KIND=dp), DIMENSION(3), INTENT(IN), OPTIONAL :: dipole
56 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
57 : OPTIONAL :: hessian, dipder
58 : REAL(KIND=dp), DIMENSION(3, 3), INTENT(IN), &
59 : OPTIONAL :: polar
60 : INTEGER, INTENT(IN), OPTIONAL :: fixed_atoms
61 :
62 : REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
63 :
64 : INTEGER :: i, j, natom, nc
65 : LOGICAL :: nddo
66 : REAL(KIND=dp) :: eout
67 : REAL(KIND=dp), DIMENSION(5) :: fz
68 : TYPE(dft_control_type), POINTER :: dft_control
69 : TYPE(qs_energy_type), POINTER :: qs_energy
70 :
71 40 : IF (iounit > 0) THEN
72 : ! the units depend on the qs method!
73 40 : CPASSERT(ASSOCIATED(force_env%qs_env))
74 40 : CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
75 40 : nddo = dft_control%qs_control%semi_empirical
76 40 : IF (nddo) THEN
77 0 : CALL get_qs_env(force_env%qs_env, energy=qs_energy)
78 0 : eout = energy + qs_energy%core_self
79 : ELSE
80 40 : eout = energy
81 : END IF
82 : !
83 40 : natom = SIZE(particles)
84 40 : IF (PRESENT(fixed_atoms)) natom = natom - fixed_atoms
85 40 : WRITE (iounit, "(A7)") "RESULTS"
86 40 : WRITE (iounit, "(A18)") "CURRENT COORDINATE"
87 147 : DO i = 1, natom
88 107 : WRITE (iounit, "(A,3F24.12)") TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol)), &
89 575 : particles(i)%r(1:3)*angstrom
90 : END DO
91 40 : WRITE (iounit, "(A8,3F18.12)") "ENERGY =", eout, zero, zero
92 40 : WRITE (iounit, "(A8,3F18.12)") " =", zero, zero, zero
93 40 : WRITE (iounit, "(A8,F18.12)") "S**2 =", zero
94 40 : WRITE (iounit, "(A8)") "GRADIENT"
95 147 : DO i = 1, natom
96 468 : WRITE (iounit, "(F17.12)") - particles(i)%f(1:3)
97 : END DO
98 40 : IF (PRESENT(dipole)) THEN
99 0 : WRITE (iounit, "(A8,3F18.12)") "DIPOLE =", dipole(1:3)
100 : ELSE
101 40 : WRITE (iounit, "(A8,3F18.12)") "DIPOLE =", zero, zero, zero
102 : END IF
103 40 : fz = zero
104 40 : WRITE (iounit, "(A7)") "HESSIAN"
105 40 : IF (PRESENT(hessian)) THEN
106 2 : nc = 3*natom
107 6 : DO i = 1, nc, 5
108 26 : DO j = i, nc
109 24 : WRITE (iounit, "(5(F13.9,1X))") hessian(j, i:MIN(j, i + 4))
110 : END DO
111 : END DO
112 : ELSE
113 38 : nc = 3*natom
114 116 : DO i = 1, nc, 5
115 549 : DO j = i, nc
116 511 : WRITE (iounit, "(5(F13.9,1X))") fz(1:MIN(j - i + 1, 5))
117 : END DO
118 : END DO
119 : END IF
120 40 : WRITE (iounit, "(A18)") "DIPOLE DERIVATIVES"
121 40 : IF (PRESENT(dipder)) THEN
122 0 : DO i = 1, 3*natom
123 0 : WRITE (iounit, "(3(F17.12,7X))") dipder(1:3, i)
124 : END DO
125 : ELSE
126 361 : DO i = 1, 3*natom
127 361 : WRITE (iounit, "(3(F17.12,7X))") zero, zero, zero
128 : END DO
129 : END IF
130 40 : WRITE (iounit, "(A14)") "POLARIZABILITY"
131 40 : IF (PRESENT(polar)) THEN
132 0 : WRITE (iounit, "(1(F17.12,7X))") polar(1, 1)
133 0 : WRITE (iounit, "(2(F17.12,7X))") polar(2, 1), polar(2, 2)
134 0 : WRITE (iounit, "(3(F17.12,7X))") polar(3, 1), polar(3, 2), polar(3, 3)
135 : ELSE
136 40 : WRITE (iounit, "(1(F17.12,7X))") zero
137 40 : WRITE (iounit, "(2(F17.12,7X))") zero, zero
138 40 : WRITE (iounit, "(3(F17.12,7X))") zero, zero, zero
139 : END IF
140 : END IF
141 :
142 40 : END SUBROUTINE write_grrm
143 :
144 : END MODULE grrm_utils
|
