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