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 Methods dealing with Neural Network interaction potential
10 : !> \author Laura Duran
11 : !> \date 2023-02-17
12 : ! **************************************************************************************************
13 : MODULE helium_nnp
14 :
15 : USE bibliography, ONLY: Behler2007,&
16 : Behler2011,&
17 : Schran2020a,&
18 : Schran2020b,&
19 : cite_reference
20 : USE cell_methods, ONLY: cell_create
21 : USE cell_types, ONLY: cell_release,&
22 : cell_type,&
23 : pbc
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_type
26 : USE cp_output_handling, ONLY: cp_p_file,&
27 : cp_print_key_finished_output,&
28 : cp_print_key_should_output,&
29 : cp_print_key_unit_nr
30 : USE cp_units, ONLY: cp_unit_from_cp2k
31 : USE helium_types, ONLY: helium_solvent_type
32 : USE input_section_types, ONLY: section_vals_get,&
33 : section_vals_get_subs_vals,&
34 : section_vals_type,&
35 : section_vals_val_get
36 : USE kinds, ONLY: default_path_length,&
37 : default_string_length,&
38 : dp
39 : USE nnp_environment, ONLY: nnp_init_model
40 : USE nnp_environment_types, ONLY: nnp_env_get,&
41 : nnp_env_set,&
42 : nnp_type
43 : USE periodic_table, ONLY: get_ptable_info
44 : USE physcon, ONLY: angstrom
45 : #include "../base/base_uses.f90"
46 :
47 : IMPLICIT NONE
48 :
49 : PRIVATE
50 :
51 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .TRUE.
52 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'helium_nnp'
53 :
54 : PUBLIC :: helium_init_nnp, &
55 : helium_nnp_print
56 :
57 : CONTAINS
58 :
59 : ! ***************************************************************************
60 : !> \brief Read and initialize all the information for neural network potentials
61 : !> \param helium ...
62 : !> \param nnp ...
63 : !> \param input ...
64 : !> \date 2023-02-21
65 : !> \author lduran
66 : ! **************************************************************************************************
67 1 : SUBROUTINE helium_init_nnp(helium, nnp, input)
68 : TYPE(helium_solvent_type), INTENT(INOUT) :: helium
69 : TYPE(nnp_type), POINTER :: nnp
70 : TYPE(section_vals_type), POINTER :: input
71 :
72 : CHARACTER(len=default_path_length) :: msg_str
73 : CHARACTER(len=default_string_length) :: elem
74 : INTEGER :: i, ig, is, j
75 : INTEGER, DIMENSION(3) :: periodicity
76 : LOGICAL :: found
77 : TYPE(cell_type), POINTER :: he_cell
78 : TYPE(cp_logger_type), POINTER :: logger
79 : TYPE(section_vals_type), POINTER :: sr_cutoff_section
80 :
81 1 : CALL cite_reference(Behler2007)
82 1 : CALL cite_reference(Behler2011)
83 1 : CALL cite_reference(Schran2020a)
84 1 : CALL cite_reference(Schran2020b)
85 :
86 1 : NULLIFY (logger)
87 1 : logger => cp_get_default_logger()
88 :
89 1 : CALL nnp_env_set(nnp_env=nnp, nnp_input=input)
90 :
91 1 : nnp%num_atoms = helium%solute_atoms + 1
92 :
93 1 : CALL nnp_init_model(nnp, "HELIUM NNP")
94 :
95 1 : periodicity = 0
96 1 : IF (helium%periodic) periodicity = 1
97 1 : NULLIFY (he_cell)
98 : CALL cell_create(he_cell, hmat=helium%cell_m, &
99 1 : periodic=periodicity, tag="HELIUM NNP")
100 1 : CALL nnp_env_set(nnp, cell=he_cell)
101 1 : CALL cell_release(he_cell)
102 :
103 : ! Set up arrays for calculation:
104 3 : ALLOCATE (nnp%ele_ind(nnp%num_atoms))
105 3 : ALLOCATE (nnp%nuc_atoms(nnp%num_atoms))
106 3 : ALLOCATE (nnp%coord(3, nnp%num_atoms))
107 3 : ALLOCATE (nnp%nnp_forces(3, nnp%num_atoms))
108 3 : ALLOCATE (nnp%atoms(nnp%num_atoms))
109 3 : ALLOCATE (nnp%sort(nnp%num_atoms - 1))
110 :
111 : !fill arrays, assume that order will not change during simulation
112 1 : ig = 1
113 1 : is = 1
114 4 : DO i = 1, nnp%n_ele
115 3 : IF (nnp%ele(i) == 'He') THEN
116 1 : nnp%atoms(ig) = 'He'
117 1 : CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
118 1 : nnp%ele_ind(ig) = i
119 1 : ig = ig + 1
120 : END IF
121 13 : DO j = 1, helium%solute_atoms
122 12 : IF (nnp%ele(i) == helium%solute_element(j)) THEN
123 3 : nnp%atoms(ig) = nnp%ele(i)
124 3 : CALL get_ptable_info(nnp%atoms(ig), number=nnp%nuc_atoms(ig))
125 3 : nnp%ele_ind(ig) = i
126 3 : nnp%sort(is) = j
127 3 : ig = ig + 1
128 3 : is = is + 1
129 : END IF
130 : END DO
131 : END DO
132 :
133 3 : ALLOCATE (helium%nnp_sr_cut(nnp%n_ele))
134 4 : helium%nnp_sr_cut = 0.0_dp
135 :
136 1 : sr_cutoff_section => section_vals_get_subs_vals(nnp%nnp_input, "SR_CUTOFF")
137 1 : CALL section_vals_get(sr_cutoff_section, n_repetition=is)
138 3 : DO i = 1, is
139 2 : CALL section_vals_val_get(sr_cutoff_section, "ELEMENT", c_val=elem, i_rep_section=i)
140 2 : found = .FALSE.
141 8 : DO ig = 1, nnp%n_ele
142 8 : IF (TRIM(nnp%ele(ig)) == TRIM(elem)) THEN
143 2 : found = .TRUE.
144 : CALL section_vals_val_get(sr_cutoff_section, "RADIUS", r_val=helium%nnp_sr_cut(ig), &
145 2 : i_rep_section=i)
146 : END IF
147 : END DO
148 3 : IF (.NOT. found) THEN
149 0 : msg_str = "SR_CUTOFF for element "//TRIM(elem)//" defined but not found in NNP"
150 0 : CPWARN(msg_str)
151 : END IF
152 : END DO
153 4 : helium%nnp_sr_cut(:) = helium%nnp_sr_cut(:)**2
154 :
155 1 : RETURN
156 :
157 1 : END SUBROUTINE helium_init_nnp
158 :
159 : ! **************************************************************************************************
160 : !> \brief Print properties according to the requests in input file
161 : !> \param nnp ...
162 : !> \param print_section ...
163 : !> \param ind_he ...
164 : !> \date 2023-07-31
165 : !> \author Laura Duran
166 : ! **************************************************************************************************
167 55172 : SUBROUTINE helium_nnp_print(nnp, print_section, ind_he)
168 : TYPE(nnp_type), INTENT(INOUT) :: nnp
169 : TYPE(section_vals_type), INTENT(IN), POINTER :: print_section
170 : INTEGER, INTENT(IN) :: ind_he
171 :
172 : INTEGER :: unit_nr
173 : LOGICAL :: file_is_new
174 : TYPE(cp_logger_type), POINTER :: logger
175 : TYPE(section_vals_type), POINTER :: print_key
176 :
177 55172 : NULLIFY (logger, print_key)
178 55172 : logger => cp_get_default_logger()
179 :
180 55172 : print_key => section_vals_get_subs_vals(print_section, "ENERGIES")
181 55172 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
182 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
183 0 : middle_name="helium-nnp-energies", is_new_file=file_is_new)
184 0 : IF (unit_nr > 0) CALL helium_nnp_print_energies(nnp, unit_nr, file_is_new)
185 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
186 : END IF
187 :
188 55172 : print_key => section_vals_get_subs_vals(print_section, "FORCES_SIGMA")
189 55172 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
190 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".data", &
191 0 : middle_name="helium-nnp-forces-std", is_new_file=file_is_new)
192 0 : IF (unit_nr > 0) CALL helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
193 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
194 : END IF
195 :
196 55172 : CALL logger%para_env%sum(nnp%output_expol)
197 55172 : IF (nnp%output_expol) THEN
198 0 : print_key => section_vals_get_subs_vals(print_section, "EXTRAPOLATION")
199 0 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
200 : unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".xyz", &
201 0 : middle_name="-NNP-He-extrapolation")
202 0 : IF (unit_nr > 0) CALL helium_nnp_print_expol(nnp, unit_nr, ind_he)
203 0 : CALL cp_print_key_finished_output(unit_nr, logger, print_key)
204 : END IF
205 : END IF
206 :
207 55172 : END SUBROUTINE helium_nnp_print
208 :
209 : ! **************************************************************************************************
210 : !> \brief Print NNP energies and standard deviation sigma
211 : !> \param nnp ...
212 : !> \param unit_nr ...
213 : !> \param file_is_new ...
214 : !> \date 2023-07-31
215 : !> \author Laura Duran
216 : ! **************************************************************************************************
217 0 : SUBROUTINE helium_nnp_print_energies(nnp, unit_nr, file_is_new)
218 : TYPE(nnp_type), INTENT(INOUT) :: nnp
219 : INTEGER, INTENT(IN) :: unit_nr
220 : LOGICAL, INTENT(IN) :: file_is_new
221 :
222 : CHARACTER(len=default_path_length) :: fmt_string
223 : INTEGER :: i
224 : REAL(KIND=dp) :: std
225 :
226 0 : IF (file_is_new) THEN
227 0 : WRITE (unit_nr, "(A1,1X,A20)", ADVANCE='no') "#", "NNP Average [a.u.],"
228 0 : WRITE (unit_nr, "(A20)", ADVANCE='no') "NNP sigma [a.u.]"
229 0 : DO i = 1, nnp%n_committee
230 0 : WRITE (unit_nr, "(A17,I3)", ADVANCE='no') "NNP", i
231 : END DO
232 0 : WRITE (unit_nr, "(A)") ""
233 : END IF
234 :
235 0 : fmt_string = "(2X,2(F20.9))"
236 0 : WRITE (fmt_string, "(A,I3,A)") "(2X", nnp%n_committee + 2, "(F20.9))"
237 0 : std = SUM((SUM(nnp%atomic_energy, 1) - nnp%nnp_potential_energy)**2)
238 0 : std = std/REAL(nnp%n_committee, dp)
239 0 : std = SQRT(std)
240 0 : WRITE (unit_nr, fmt_string) nnp%nnp_potential_energy, std, SUM(nnp%atomic_energy, 1)
241 :
242 0 : END SUBROUTINE helium_nnp_print_energies
243 :
244 : ! **************************************************************************************************
245 : !> \brief Print standard deviation sigma of NNP forces
246 : !> \param nnp ...
247 : !> \param unit_nr ...
248 : !> \param file_is_new ...
249 : !> \date 2023-07-31
250 : !> \author Laura Duran
251 : ! **************************************************************************************************
252 0 : SUBROUTINE helium_nnp_print_force_sigma(nnp, unit_nr, file_is_new)
253 : TYPE(nnp_type), INTENT(INOUT) :: nnp
254 : INTEGER, INTENT(IN) :: unit_nr
255 : LOGICAL, INTENT(IN) :: file_is_new
256 :
257 : INTEGER :: i, ig, j
258 : REAL(KIND=dp), DIMENSION(3) :: var
259 :
260 0 : IF (unit_nr > 0) THEN
261 0 : IF (file_is_new) THEN
262 0 : WRITE (unit_nr, "(A,1X,A)") "# NNP sigma of forces [a.u.] x, y, z coordinates"
263 : END IF
264 :
265 0 : ig = 1
266 0 : DO i = 1, nnp%num_atoms
267 0 : IF (nnp%ele(i) == 'He') THEN
268 0 : var = 0.0_dp
269 0 : DO j = 1, nnp%n_committee
270 0 : var = var + (nnp%committee_forces(:, i, j) - nnp%nnp_forces(:, i))**2
271 : END DO
272 0 : WRITE (unit_nr, "(A4,1X,3E20.10)") nnp%ele(i), var
273 : END IF
274 0 : ig = ig + 1
275 : END DO
276 : END IF
277 :
278 0 : END SUBROUTINE helium_nnp_print_force_sigma
279 :
280 : ! **************************************************************************************************
281 : !> \brief Print structures with extrapolation warning
282 : !> \param nnp ...
283 : !> \param unit_nr ...
284 : !> \param ind_he ...
285 : !> \date 2023-10-11
286 : !> \author Harald Forbert (harald.forbert@rub.de)
287 : ! **************************************************************************************************
288 0 : SUBROUTINE helium_nnp_print_expol(nnp, unit_nr, ind_he)
289 : TYPE(nnp_type), INTENT(INOUT) :: nnp
290 : INTEGER, INTENT(IN) :: unit_nr, ind_he
291 :
292 : CHARACTER(len=default_path_length) :: fmt_string
293 : INTEGER :: i
294 : REAL(KIND=dp) :: mass, unit_conv
295 : REAL(KIND=dp), DIMENSION(3) :: com
296 : TYPE(cell_type), POINTER :: cell
297 :
298 0 : NULLIFY (cell)
299 0 : CALL nnp_env_get(nnp_env=nnp, cell=cell)
300 :
301 0 : nnp%expol = nnp%expol + 1
302 0 : WRITE (unit_nr, *) nnp%num_atoms
303 0 : WRITE (unit_nr, "(A,1X,I6)") "HELIUM-NNP extrapolation point N =", nnp%expol
304 :
305 : ! move to COM of solute and wrap the box
306 : ! coord not needed afterwards, therefore manipulation ok
307 0 : com = 0.0_dp
308 0 : mass = 0.0_dp
309 0 : DO i = 1, nnp%num_atoms
310 0 : IF (i == ind_he) CYCLE
311 0 : CALL get_ptable_info(nnp%atoms(i), amass=unit_conv)
312 0 : com(:) = com(:) + nnp%coord(:, i)*unit_conv
313 0 : mass = mass + unit_conv
314 : END DO
315 0 : com(:) = com(:)/mass
316 :
317 0 : DO i = 1, nnp%num_atoms
318 0 : nnp%coord(:, i) = nnp%coord(:, i) - com(:)
319 0 : nnp%coord(:, i) = pbc(nnp%coord(:, i), cell)
320 : END DO
321 :
322 0 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
323 0 : fmt_string = "(A4,1X,3F20.10)"
324 0 : DO i = 1, nnp%num_atoms
325 : WRITE (unit_nr, fmt_string) &
326 0 : nnp%atoms(i), &
327 0 : nnp%coord(1, i)*unit_conv, &
328 0 : nnp%coord(2, i)*unit_conv, &
329 0 : nnp%coord(3, i)*unit_conv
330 : END DO
331 :
332 0 : END SUBROUTINE helium_nnp_print_expol
333 :
334 : END MODULE helium_nnp
|
