-
Notifications
You must be signed in to change notification settings - Fork 1
/
scine_utils.F
172 lines (157 loc) · 8.15 KB
/
scine_utils.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief SCINE interface
!> \author JGH - 01.2020
! **************************************************************************************************
MODULE scine_utils
USE cell_types, ONLY: cell_type
USE cp_control_types, ONLY: dft_control_type
USE force_env_types, ONLY: force_env_type
USE kinds, ONLY: dp
USE particle_types, ONLY: particle_type
USE physcon, ONLY: angstrom
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'scine_utils'
PUBLIC :: write_scine
! **************************************************************************************************
CONTAINS
! **************************************************************************************************
!> \brief Write SCINE interface file
!>
!> \param iounit ...
!> \param force_env ...
!> \param particles ...
!> \param energy ...
!> \param hessian ...
! **************************************************************************************************
SUBROUTINE write_scine(iounit, force_env, particles, energy, hessian)
INTEGER, INTENT(IN) :: iounit
TYPE(force_env_type), POINTER :: force_env
TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particles
REAL(KIND=dp), INTENT(IN) :: energy
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN), &
OPTIONAL :: hessian
REAL(KIND=dp), PARAMETER :: zero = 0.0_dp
CHARACTER(LEN=20) :: sunit
INTEGER :: i, j, natom, nc
LOGICAL :: nddo
REAL(KIND=dp) :: eout
REAL(KIND=dp), DIMENSION(5) :: fz
TYPE(cell_type), POINTER :: cell
TYPE(dft_control_type), POINTER :: dft_control
TYPE(qs_energy_type), POINTER :: qs_energy
IF (iounit > 0) THEN
! the units depend on the qs method!
CPASSERT(ASSOCIATED(force_env%qs_env))
CALL get_qs_env(force_env%qs_env, dft_control=dft_control)
nddo = dft_control%qs_control%semi_empirical
IF (nddo) THEN
CALL get_qs_env(force_env%qs_env, energy=qs_energy)
eout = energy + qs_energy%core_self
ELSE
eout = energy
END IF
!
CALL get_qs_env(force_env%qs_env, cell=cell)
natom = SIZE(particles)
WRITE (iounit, "(A,A)") "title: ", "'SCINE Interface File'"
sunit = "'hartree'"
WRITE (iounit, "(A,F18.12,A,A)") "energy: [", eout, ", "//TRIM(sunit), " ]"
sunit = "'angstrom'"
WRITE (iounit, "(A)") "system:"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
WRITE (iounit, "(T2,A)") "cell: "
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- A: [", cell%hmat(1:3, 1)*angstrom, " ]"
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- B: [", cell%hmat(1:3, 2)*angstrom, " ]"
WRITE (iounit, "(T4,A,F24.12,',',F24.12,',',F24.12,A)") &
"- C: [", cell%hmat(1:3, 3)*angstrom, " ]"
WRITE (iounit, "(T2,A,L1,', ',L1,', ',L1,' ]')") "periodicity: [ ", (cell%perd == 1)
WRITE (iounit, "(T2,A)") "coordinates: "
DO i = 1, natom
WRITE (iounit, "(T4,A,A2,A,F24.12,',',F24.12,',',F24.12,A)") &
"- ", TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":", &
" [", particles(i)%r(1:3)*angstrom, " ]"
END DO
WRITE (iounit, "(A)") "gradient:"
sunit = "'hartree/bohr'"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
WRITE (iounit, "(T2,A)") "values: "
DO i = 1, natom
WRITE (iounit, "(T4,A,A2,A,F24.12,',',F24.12,',',F24.12,A)") &
"- ", TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":", &
" [", particles(i)%f(1:3), " ]"
END DO
fz = zero
IF (PRESENT(hessian)) THEN
WRITE (iounit, "(A)") "hessian:"
sunit = "'hartree/bohr^2'"
WRITE (iounit, "(T2,A,A)") "unit: ", TRIM(sunit)
DO i = 1, natom
WRITE (iounit, "(T4,A)") TRIM(ADJUSTL(particles(i)%atomic_kind%element_symbol))//":"
nc = 3*(i - 1) + 1
WRITE (iounit, "(T6,'X: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
nc = 3*(i - 1) + 2
WRITE (iounit, "(T6,'Y: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
nc = 3*(i - 1) + 3
WRITE (iounit, "(T6,'Z: [')")
DO j = 1, nc, 5
IF (nc > j + 4) THEN
WRITE (iounit, "(T12,5(F20.12,','))") hessian(nc, j:j + 4)
ELSEIF (nc == j + 4) THEN
WRITE (iounit, "(T12,4(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 3) THEN
WRITE (iounit, "(T12,3(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 2) THEN
WRITE (iounit, "(T12,2(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j + 1) THEN
WRITE (iounit, "(T12,1(F20.12,','),F20.12,' ]')") hessian(nc, j:nc)
ELSEIF (nc == j) THEN
WRITE (iounit, "(T12,F20.12,' ]')") hessian(nc, j:nc)
END IF
END DO
END DO
END IF
END IF
END SUBROUTINE write_scine
END MODULE scine_utils