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 Interface for Voronoi Integration and output of BQB files
10 : !> \par History
11 : !> 11/2020 created [mbrehm]
12 : !> \author Martin Brehm
13 : ! **************************************************************************************************
14 : MODULE voronoi_interface
15 : USE input_section_types, ONLY: section_vals_type, &
16 : section_vals_val_get
17 : USE cp_log_handling, ONLY: cp_get_default_logger, &
18 : cp_logger_get_default_io_unit, &
19 : cp_logger_type
20 : USE bibliography, ONLY: Rycroft2009, Thomas2015, Brehm2018, Brehm2020, &
21 : Brehm2021, cite_reference
22 : USE kinds, ONLY: dp, default_path_length
23 : USE cell_types, ONLY: cell_type, pbc
24 : USE pw_types, ONLY: pw_r3d_rs_type
25 : USE physcon, ONLY: bohr, debye
26 : USE mathconstants, ONLY: fourpi
27 : USE orbital_pointers, ONLY: indco
28 : USE qs_environment_types, ONLY: get_qs_env, &
29 : qs_environment_type
30 : USE molecule_kind_types, ONLY: molecule_kind_type, &
31 : write_molecule_kind_set
32 : USE molecule_types, ONLY: molecule_type
33 : USE qs_rho_types, ONLY: qs_rho_get, &
34 : qs_rho_type
35 : USE atomic_kind_types, ONLY: atomic_kind_type, &
36 : get_atomic_kind
37 : USE particle_list_types, ONLY: particle_list_type
38 : USE particle_types, ONLY: particle_type
39 : USE cp_files, ONLY: file_exists, close_file, open_file
40 : USE qs_kind_types, ONLY: get_qs_kind, &
41 : qs_kind_type
42 : USE message_passing, ONLY: mp_para_env_type
43 : USE qs_subsys_types, ONLY: qs_subsys_get, &
44 : qs_subsys_type
45 : USE cp_control_types, ONLY: dft_control_type
46 : USE pw_grid_types, ONLY: PW_MODE_LOCAL
47 : USE physcon, ONLY: angstrom, femtoseconds
48 : USE message_passing, ONLY: mp_comm_type
49 : USE input_constants, ONLY: &
50 : voro_radii_unity, voro_radii_vdw, voro_radii_cov, voro_radii_user, &
51 : bqb_opt_off, bqb_opt_quick, bqb_opt_normal, bqb_opt_patient, bqb_opt_exhaustive, do_method_gapw
52 : USE qs_rho0_types, ONLY: rho0_atom_type, &
53 : rho0_mpole_type
54 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE, C_CHAR
55 :
56 : #if defined(__HAS_IEEE_EXCEPTIONS)
57 : USE ieee_exceptions, ONLY: ieee_get_halting_mode, &
58 : ieee_set_halting_mode, &
59 : ieee_all
60 : #endif
61 :
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 : PRIVATE
66 :
67 : ! Global parameters
68 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'voronoi_interface'
69 : PUBLIC :: entry_voronoi_or_bqb, finalize_libvori
70 : INTEGER :: step_count = 0
71 :
72 : #if defined(__LIBVORI)
73 :
74 : ! The C interface to libvori
75 :
76 : INTERFACE
77 :
78 : INTEGER(C_INT) FUNCTION libvori_setBQBSkipFirst(i) BIND(C, NAME='libvori_setBQBSkipFirst')
79 : USE ISO_C_BINDING, ONLY: C_INT
80 : INTEGER(C_INT), VALUE :: i
81 :
82 : END FUNCTION libvori_setBQBSkipFirst
83 :
84 : INTEGER(C_INT) FUNCTION libvori_setBQBStoreStep(i) BIND(C, NAME='libvori_setBQBStoreStep')
85 : USE ISO_C_BINDING, ONLY: C_INT
86 : INTEGER(C_INT), VALUE :: i
87 :
88 : END FUNCTION libvori_setBQBStoreStep
89 :
90 : INTEGER(C_INT) FUNCTION libvori_setVoronoiSkipFirst(i) BIND(C, NAME='libvori_setVoronoiSkipFirst')
91 : USE ISO_C_BINDING, ONLY: C_INT
92 : INTEGER(C_INT), VALUE :: i
93 :
94 : END FUNCTION libvori_setVoronoiSkipFirst
95 :
96 : INTEGER(C_INT) FUNCTION libvori_setBQBCheck(i) BIND(C, NAME='libvori_setBQBCheck')
97 : USE ISO_C_BINDING, ONLY: C_INT
98 : INTEGER(C_INT), VALUE :: i
99 :
100 : END FUNCTION libvori_setBQBCheck
101 :
102 : INTEGER(C_INT) FUNCTION libvori_setBQBFilename(len, s) BIND(C, NAME='libvori_setBQBFilename')
103 : USE ISO_C_BINDING, ONLY: C_INT, C_CHAR
104 : INTEGER(C_INT), VALUE :: len
105 : CHARACTER(C_CHAR) :: s(*)
106 :
107 : END FUNCTION libvori_setBQBFilename
108 :
109 : INTEGER(C_INT) FUNCTION libvori_setBQBParmString(len, s) BIND(C, NAME='libvori_setBQBParmString')
110 : USE ISO_C_BINDING, ONLY: C_INT, C_CHAR
111 : INTEGER(C_INT), VALUE :: len
112 : CHARACTER(C_CHAR) :: s(*)
113 :
114 : END FUNCTION libvori_setBQBParmString
115 :
116 : INTEGER(C_INT) FUNCTION libvori_setBQBHistory(i) BIND(C, NAME='libvori_setBQBHistory')
117 : USE ISO_C_BINDING, ONLY: C_INT
118 : INTEGER(C_INT), VALUE :: i
119 :
120 : END FUNCTION libvori_setBQBHistory
121 :
122 : INTEGER(C_INT) FUNCTION libvori_setBQBOptimization(i) BIND(C, NAME='libvori_setBQBOptimization')
123 : USE ISO_C_BINDING, ONLY: C_INT
124 : INTEGER(C_INT), VALUE :: i
125 :
126 : END FUNCTION libvori_setBQBOptimization
127 :
128 : INTEGER(C_INT) FUNCTION libvori_processBQBFrame(step, t) BIND(C, NAME='libvori_processBQBFrame')
129 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
130 : INTEGER(C_INT), VALUE :: step
131 : REAL(C_DOUBLE), VALUE :: t
132 :
133 : END FUNCTION libvori_processBQBFrame
134 :
135 : INTEGER(C_INT) FUNCTION libvori_setPrefix_Voronoi() BIND(C, NAME='libvori_setPrefix_Voronoi')
136 : USE ISO_C_BINDING, ONLY: C_INT
137 : END FUNCTION libvori_setPrefix_Voronoi
138 :
139 : INTEGER(C_INT) FUNCTION libvori_setPrefix_BQB() BIND(C, NAME='libvori_setPrefix_BQB')
140 : USE ISO_C_BINDING, ONLY: C_INT
141 : END FUNCTION libvori_setPrefix_BQB
142 :
143 : INTEGER(C_INT) FUNCTION libvori_setRefinementFactor(i) BIND(C, NAME='libvori_setRefinementFactor')
144 : USE ISO_C_BINDING, ONLY: C_INT
145 : INTEGER(C_INT), VALUE :: i
146 :
147 : END FUNCTION libvori_setRefinementFactor
148 :
149 : INTEGER(C_INT) FUNCTION libvori_setJitter(i) BIND(C, NAME='libvori_setJitter')
150 : USE ISO_C_BINDING, ONLY: C_INT
151 : INTEGER(C_INT), VALUE :: i
152 :
153 : END FUNCTION libvori_setJitter
154 :
155 : INTEGER(C_INT) FUNCTION libvori_setJitterAmplitude(amp) BIND(C, NAME='libvori_setJitterAmplitude')
156 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
157 : REAL(C_DOUBLE), VALUE :: amp
158 :
159 : END FUNCTION libvori_setJitterAmplitude
160 :
161 : INTEGER(C_INT) FUNCTION libvori_setJitterSeed(seed) BIND(C, NAME='libvori_setJitterSeed')
162 : USE ISO_C_BINDING, ONLY: C_INT
163 : INTEGER(C_INT), VALUE :: seed
164 :
165 : END FUNCTION libvori_setJitterSeed
166 :
167 : INTEGER(C_INT) FUNCTION libvori_setEMPOutput(i) BIND(C, NAME='libvori_setEMPOutput')
168 : USE ISO_C_BINDING, ONLY: C_INT
169 : INTEGER(C_INT), VALUE :: i
170 :
171 : END FUNCTION libvori_setEMPOutput
172 :
173 : INTEGER(C_INT) FUNCTION libvori_setPrintLevel_Verbose() BIND(C, NAME='libvori_setPrintLevel_Verbose')
174 : USE ISO_C_BINDING, ONLY: C_INT
175 : END FUNCTION libvori_setPrintLevel_Verbose
176 :
177 : INTEGER(C_INT) FUNCTION libvori_setRadii_Unity() BIND(C, NAME='libvori_setRadii_Unity')
178 : USE ISO_C_BINDING, ONLY: C_INT
179 : END FUNCTION libvori_setRadii_Unity
180 :
181 : INTEGER(C_INT) FUNCTION libvori_setRadii_Covalent() BIND(C, NAME='libvori_setRadii_Covalent')
182 : USE ISO_C_BINDING, ONLY: C_INT
183 : END FUNCTION libvori_setRadii_Covalent
184 :
185 : INTEGER(C_INT) FUNCTION libvori_setRadii_User(factor, rad) BIND(C, NAME='libvori_setRadii_User')
186 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
187 : REAL(C_DOUBLE), VALUE :: factor
188 : REAL(C_DOUBLE) :: rad(*)
189 :
190 : END FUNCTION libvori_setRadii_User
191 :
192 : INTEGER(C_INT) FUNCTION libvori_step(step, t) BIND(C, NAME='libvori_step')
193 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
194 : INTEGER(C_INT), VALUE :: step
195 : REAL(C_DOUBLE), VALUE :: t
196 :
197 : END FUNCTION libvori_step
198 :
199 : INTEGER(C_INT) FUNCTION libvori_sanitycheck(step, t) BIND(C, NAME='libvori_sanitycheck')
200 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
201 : INTEGER(C_INT), VALUE :: step
202 : REAL(C_DOUBLE), VALUE :: t
203 :
204 : END FUNCTION libvori_sanitycheck
205 :
206 : INTEGER(C_INT) FUNCTION libvori_setGrid(rx, ry, rz, ax, ay, az, bx, by, bz, cx, cy, cz, tax, tay, taz, tbx, tby, tbz, &
207 : tcx, tcy, tcz) BIND(C, NAME='libvori_setGrid')
208 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
209 : INTEGER(C_INT), VALUE :: rx, ry, rz
210 : REAL(C_DOUBLE), VALUE :: ax, ay, az, bx, by, bz, cx, cy, cz, tax, &
211 : tay, taz, tbx, tby, tbz, tcx, tcy, tcz
212 :
213 : END FUNCTION libvori_setGrid
214 :
215 : INTEGER(C_INT) FUNCTION libvori_pushAtoms(n, pord, pchg, posx, posy, posz) BIND(C, NAME='libvori_pushAtoms')
216 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
217 : INTEGER(C_INT), VALUE :: n
218 : INTEGER(C_INT) :: pord(*)
219 : REAL(C_DOUBLE) :: pchg(*), posx(*), posy(*), posz(*)
220 :
221 : END FUNCTION libvori_pushAtoms
222 :
223 : INTEGER(C_INT) FUNCTION libvori_push_rho_zrow(ix, iy, length, buf) BIND(C, NAME='libvori_push_rho_zrow')
224 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
225 : INTEGER(C_INT), VALUE :: ix, iy, length
226 : REAL(C_DOUBLE) :: buf(*)
227 :
228 : END FUNCTION libvori_push_rho_zrow
229 :
230 : INTEGER(C_INT) FUNCTION libvori_setBQBOverwrite(i) BIND(C, NAME='libvori_setBQBOverwrite')
231 : USE ISO_C_BINDING, ONLY: C_INT
232 : INTEGER(C_INT), VALUE :: i
233 :
234 : END FUNCTION libvori_setBQBOverwrite
235 :
236 : INTEGER(C_INT) FUNCTION libvori_setVoriOverwrite(i) BIND(C, NAME='libvori_setVoriOverwrite')
237 : USE ISO_C_BINDING, ONLY: C_INT
238 : INTEGER(C_INT), VALUE :: i
239 :
240 : END FUNCTION libvori_setVoriOverwrite
241 :
242 : INTEGER(C_INT) FUNCTION libvori_get_radius(length, buf) BIND(C, NAME='libvori_get_radius')
243 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
244 : INTEGER(C_INT), VALUE :: length
245 : REAL(C_DOUBLE) :: buf(*)
246 :
247 : END FUNCTION libvori_get_radius
248 :
249 : INTEGER(C_INT) FUNCTION libvori_get_volume(length, buf) BIND(C, NAME='libvori_get_volume')
250 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
251 : INTEGER(C_INT), VALUE :: length
252 : REAL(C_DOUBLE) :: buf(*)
253 :
254 : END FUNCTION libvori_get_volume
255 :
256 : INTEGER(C_INT) FUNCTION libvori_get_charge(length, buf) BIND(C, NAME='libvori_get_charge')
257 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
258 : INTEGER(C_INT), VALUE :: length
259 : REAL(C_DOUBLE) :: buf(*)
260 :
261 : END FUNCTION libvori_get_charge
262 :
263 : INTEGER(C_INT) FUNCTION libvori_get_dipole(component, length, buf) BIND(C, NAME='libvori_get_dipole')
264 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
265 : INTEGER(C_INT), VALUE :: component, length
266 : REAL(C_DOUBLE) :: buf(*)
267 :
268 : END FUNCTION libvori_get_dipole
269 :
270 : INTEGER(C_INT) FUNCTION libvori_get_quadrupole(component, length, buf) BIND(C, NAME='libvori_get_quadrupole')
271 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
272 : INTEGER(C_INT), VALUE :: component, length
273 : REAL(C_DOUBLE) :: buf(*)
274 :
275 : END FUNCTION libvori_get_quadrupole
276 :
277 : INTEGER(C_INT) FUNCTION libvori_get_wrapped_pos(component, length, buf) BIND(C, NAME='libvori_get_wrapped_pos')
278 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
279 : INTEGER(C_INT), VALUE :: component, length
280 : REAL(C_DOUBLE) :: buf(*)
281 :
282 : END FUNCTION libvori_get_wrapped_pos
283 :
284 : INTEGER(C_INT) FUNCTION libvori_get_charge_center(component, length, buf) BIND(C, NAME='libvori_get_charge_center')
285 : USE ISO_C_BINDING, ONLY: C_INT, C_DOUBLE
286 : INTEGER(C_INT), VALUE :: component, length
287 : REAL(C_DOUBLE) :: buf(*)
288 :
289 : END FUNCTION libvori_get_charge_center
290 :
291 : INTEGER(C_INT) FUNCTION libvori_finalize() BIND(C, NAME='libvori_finalize')
292 : USE ISO_C_BINDING, ONLY: C_INT
293 : END FUNCTION libvori_finalize
294 :
295 : END INTERFACE
296 :
297 : #endif
298 :
299 : ! **************************************************************************************************
300 :
301 : CONTAINS
302 :
303 : ! **************************************************************************************************
304 : !> \brief Does a Voronoi integration of density or stores the density to compressed BQB format
305 : !> \param do_voro whether a Voronoi integration shall be performed
306 : !> \param do_bqb whether the density shall be written to compressed BQB format
307 : !> \param input_voro the input section for Voronoi integration
308 : !> \param input_bqb the input section for the BQB compression
309 : !> \param unit_voro the output unit number for the Voronoi integration results
310 : !> \param qs_env the qs_env where to calculate the charges
311 : !> \param rspace_pw the grid with the real-space electron density to integrate/compress
312 : !> \author Martin Brehm
313 : ! **************************************************************************************************
314 30 : SUBROUTINE entry_voronoi_or_bqb(do_voro, do_bqb, input_voro, input_bqb, unit_voro, qs_env, rspace_pw)
315 : INTEGER :: do_voro, do_bqb
316 : TYPE(section_vals_type), POINTER :: input_voro, input_bqb
317 : INTEGER, INTENT(IN) :: unit_voro
318 : TYPE(qs_environment_type), POINTER :: qs_env
319 : TYPE(pw_r3d_rs_type) :: rspace_pw
320 :
321 : #if defined(__LIBVORI)
322 :
323 : CHARACTER(len=*), PARAMETER :: routineN = 'entry_voronoi_or_bqb'
324 : INTEGER :: handle, iounit, &
325 : ret, i, tag, &
326 : nkind, natom, ikind, iat, ord, source, dest, &
327 : ip, i1, i2, reffac, radius_type, bqb_optimize, &
328 : bqb_history, nspins, jitter_seed
329 : LOGICAL :: outemp, bqb_skip_first, voro_skip_first, &
330 : bqb_store_step, bqb_check, voro_sanity, &
331 : bqb_overwrite, vori_overwrite, molprop, &
332 : gapw, jitter
333 : REAL(KIND=dp) :: zeff, qa, fn0, fn1, jitter_amplitude
334 : TYPE(qs_rho_type), POINTER :: rho
335 : TYPE(cp_logger_type), POINTER :: logger
336 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: buf
337 30 : INTEGER, ALLOCATABLE, DIMENSION(:) :: particles_z
338 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: particles_r
339 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: particles_c
340 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: particles_radius
341 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
342 : TYPE(atomic_kind_type), POINTER :: atomic_kind
343 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
344 : TYPE(mp_para_env_type), POINTER :: para_env
345 : TYPE(particle_list_type), POINTER :: particles
346 : REAL(KIND=dp) :: r
347 : TYPE(qs_subsys_type), POINTER :: subsys
348 30 : INTEGER, DIMENSION(:), POINTER :: atom_list
349 : TYPE(dft_control_type), POINTER :: dft_control
350 : TYPE(cell_type), POINTER :: cell
351 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: voro_radii
352 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: voro_charge
353 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: voro_volume
354 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: voro_dipole
355 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: voro_quadrupole
356 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: voro_buffer
357 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: voro_wrapped_pos
358 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: voro_charge_center
359 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: user_radii
360 : CHARACTER(len=default_path_length) :: bqb_file_name, mp_file_name
361 : CHARACTER(len=128) :: bqb_parm_string
362 : TYPE(rho0_mpole_type), POINTER :: rho0_mpole
363 : #if defined(__HAS_IEEE_EXCEPTIONS)
364 : LOGICAL, DIMENSION(5) :: halt
365 : #endif
366 :
367 30 : CALL timeset(routineN, handle)
368 30 : NULLIFY (logger)
369 30 : logger => cp_get_default_logger()
370 30 : iounit = cp_logger_get_default_io_unit(logger)
371 :
372 : CALL get_qs_env(qs_env=qs_env, rho=rho, qs_kind_set=qs_kind_set, &
373 : atomic_kind_set=atomic_kind_set, para_env=para_env, &
374 : nkind=nkind, natom=natom, subsys=subsys, dft_control=dft_control, &
375 30 : cell=cell)
376 :
377 30 : tag = 1
378 :
379 30 : IF (do_voro /= 0) THEN
380 28 : CALL section_vals_val_get(input_voro, "REFINEMENT_FACTOR", i_val=reffac)
381 28 : CALL section_vals_val_get(input_voro, "OUTPUT_EMP", l_val=outemp)
382 28 : CALL section_vals_val_get(input_voro, "JITTER", l_val=jitter)
383 28 : CALL section_vals_val_get(input_voro, "JITTER_AMPLITUDE", r_val=jitter_amplitude)
384 28 : CALL section_vals_val_get(input_voro, "JITTER_SEED", i_val=jitter_seed)
385 28 : CALL section_vals_val_get(input_voro, "VORONOI_RADII", i_val=radius_type)
386 28 : CALL section_vals_val_get(input_voro, "SKIP_FIRST", l_val=voro_skip_first)
387 28 : CALL section_vals_val_get(input_voro, "SANITY_CHECK", l_val=voro_sanity)
388 28 : CALL section_vals_val_get(input_voro, "OVERWRITE", l_val=vori_overwrite)
389 28 : IF (radius_type == voro_radii_user) THEN
390 0 : CALL section_vals_val_get(input_voro, "USER_RADII", r_vals=user_radii)
391 : END IF
392 28 : IF (qs_env%single_point_run) THEN
393 20 : voro_skip_first = .FALSE.
394 : END IF
395 28 : CALL cite_reference(Rycroft2009)
396 28 : CALL cite_reference(Thomas2015)
397 28 : CALL cite_reference(Brehm2018)
398 28 : CALL cite_reference(Brehm2020)
399 28 : CALL cite_reference(Brehm2021)
400 : !
401 28 : CALL section_vals_val_get(input_voro, "MOLECULAR_PROPERTIES", l_val=molprop)
402 28 : CALL section_vals_val_get(input_voro, "MOLPROP_FILE_NAME", c_val=mp_file_name)
403 : END IF
404 :
405 30 : IF (do_bqb /= 0) THEN
406 2 : CALL section_vals_val_get(input_bqb, "HISTORY", i_val=bqb_history)
407 2 : CALL section_vals_val_get(input_bqb, "OPTIMIZE", i_val=bqb_optimize)
408 2 : CALL section_vals_val_get(input_bqb, "FILENAME", c_val=bqb_file_name)
409 2 : CALL section_vals_val_get(input_bqb, "SKIP_FIRST", l_val=bqb_skip_first)
410 2 : CALL section_vals_val_get(input_bqb, "STORE_STEP_NUMBER", l_val=bqb_store_step)
411 2 : CALL section_vals_val_get(input_bqb, "CHECK", l_val=bqb_check)
412 2 : CALL section_vals_val_get(input_bqb, "OVERWRITE", l_val=bqb_overwrite)
413 2 : CALL section_vals_val_get(input_bqb, "PARAMETER_KEY", c_val=bqb_parm_string)
414 2 : IF (qs_env%single_point_run) THEN
415 2 : bqb_skip_first = .FALSE.
416 2 : bqb_history = 1
417 : END IF
418 2 : IF (bqb_history < 1) THEN
419 0 : bqb_history = 1
420 : END IF
421 2 : CALL cite_reference(Brehm2018)
422 : END IF
423 :
424 30 : CALL qs_subsys_get(subsys, particles=particles)
425 :
426 : ! Temporarily disable floating point traps because libvori raise IEEE754 exceptions.
427 : #if defined(__HAS_IEEE_EXCEPTIONS)
428 : CALL ieee_get_halting_mode(IEEE_ALL, halt)
429 : CALL ieee_set_halting_mode(IEEE_ALL, .FALSE.)
430 : #endif
431 :
432 : ASSOCIATE (ionode => para_env%is_source(), my_rank => para_env%mepos, &
433 : num_pe => para_env%num_pe, &
434 : sim_step => qs_env%sim_step, sim_time => qs_env%sim_time, &
435 : L1 => rspace_pw%pw_grid%bounds(1, 1), L2 => rspace_pw%pw_grid%bounds(1, 2), &
436 : L3 => rspace_pw%pw_grid%bounds(1, 3), U1 => rspace_pw%pw_grid%bounds(2, 1), &
437 : U2 => rspace_pw%pw_grid%bounds(2, 2), U3 => rspace_pw%pw_grid%bounds(2, 3))
438 30 : IF (ionode) THEN
439 :
440 15 : IF (iounit > 0) THEN
441 15 : WRITE (iounit, *) ""
442 : END IF
443 :
444 15 : IF (do_voro /= 0) THEN
445 14 : ret = libvori_setPrefix_Voronoi()
446 14 : ret = libvori_setRefinementFactor(reffac)
447 14 : ret = libvori_setJitter(MERGE(1, 0, jitter))
448 14 : ret = libvori_setJitterAmplitude(jitter_amplitude*angstrom)
449 14 : ret = libvori_setJitterSeed(jitter_seed)
450 28 : ret = libvori_setVoronoiSkipFirst(MERGE(1, 0, voro_skip_first))
451 28 : ret = libvori_setVoriOverwrite(MERGE(1, 0, vori_overwrite))
452 23 : ret = libvori_setEMPOutput(MERGE(1, 0, outemp))
453 : ELSE
454 1 : ret = libvori_setPrefix_BQB()
455 : END IF
456 :
457 15 : IF (do_bqb /= 0) THEN
458 1 : ret = libvori_setBQBFilename(default_path_length, bqb_file_name)
459 1 : ret = libvori_setBQBParmString(128, bqb_parm_string)
460 0 : SELECT CASE (bqb_optimize)
461 : CASE (bqb_opt_off)
462 0 : bqb_optimize = 0
463 : CASE (bqb_opt_quick)
464 1 : bqb_optimize = 1
465 : CASE (bqb_opt_normal)
466 0 : bqb_optimize = 2
467 : CASE (bqb_opt_patient)
468 0 : bqb_optimize = 3
469 : CASE (bqb_opt_exhaustive)
470 1 : bqb_optimize = 4
471 : END SELECT
472 1 : ret = libvori_setBQBOptimization(bqb_optimize)
473 1 : ret = libvori_setBQBHistory(bqb_history)
474 2 : ret = libvori_setBQBSkipFirst(MERGE(1, 0, bqb_skip_first))
475 1 : ret = libvori_setBQBCheck(MERGE(1, 0, bqb_check))
476 2 : ret = libvori_setBQBOverwrite(MERGE(1, 0, bqb_overwrite))
477 1 : ret = libvori_setBQBStoreStep(MERGE(1, 0, bqb_store_step))
478 : END IF
479 :
480 : ret = libvori_setgrid( &
481 : U1 - L1 + 1, &
482 : U2 - L2 + 1, &
483 : U3 - L3 + 1, &
484 : rspace_pw%pw_grid%dh(1, 1)*(U1 - L1 + 1), &
485 : rspace_pw%pw_grid%dh(2, 1)*(U1 - L1 + 1), &
486 : rspace_pw%pw_grid%dh(3, 1)*(U1 - L1 + 1), &
487 : rspace_pw%pw_grid%dh(1, 2)*(U2 - L2 + 1), &
488 : rspace_pw%pw_grid%dh(2, 2)*(U2 - L2 + 1), &
489 : rspace_pw%pw_grid%dh(3, 2)*(U2 - L2 + 1), &
490 : rspace_pw%pw_grid%dh(1, 3)*(U3 - L3 + 1), &
491 : rspace_pw%pw_grid%dh(2, 3)*(U3 - L3 + 1), &
492 : rspace_pw%pw_grid%dh(3, 3)*(U3 - L3 + 1), &
493 : cell%hmat(1, 1), &
494 : cell%hmat(2, 1), &
495 : cell%hmat(3, 1), &
496 : cell%hmat(1, 2), &
497 : cell%hmat(2, 2), &
498 : cell%hmat(3, 2), &
499 : cell%hmat(1, 3), &
500 : cell%hmat(2, 3), &
501 : cell%hmat(3, 3) &
502 15 : )
503 :
504 15 : IF (ret /= 0) THEN
505 0 : CPABORT("The library returned an error. Aborting.")
506 : END IF
507 :
508 45 : ALLOCATE (particles_z(natom))
509 45 : ALLOCATE (particles_c(natom))
510 45 : ALLOCATE (particles_r(3, natom))
511 30 : ALLOCATE (particles_radius(natom))
512 :
513 46 : DO ikind = 1, nkind
514 31 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, vdw_radius=r)
515 31 : r = r*angstrom
516 31 : atomic_kind => atomic_kind_set(ikind)
517 31 : CALL get_atomic_kind(atomic_kind, atom_list=atom_list, z=ord)
518 105 : DO iat = 1, SIZE(atom_list)
519 59 : i = atom_list(iat)
520 59 : particles_c(i) = zeff
521 59 : particles_z(i) = ord
522 236 : particles_r(:, i) = particles%els(i)%r(:)
523 90 : particles_radius(i) = r
524 : END DO
525 : END DO
526 :
527 369 : ret = libvori_pushatoms(natom, particles_z, particles_c, particles_r(1, :), particles_r(2, :), particles_r(3, :))
528 :
529 15 : IF (ret /= 0) THEN
530 0 : CPABORT("The library returned an error. Aborting.")
531 : END IF
532 :
533 : END IF
534 :
535 30 : IF (iounit > 0) THEN
536 15 : IF (do_voro /= 0) THEN
537 14 : WRITE (iounit, *) "VORONOI| Collecting electron density from MPI ranks and sending to library..."
538 : ELSE
539 1 : WRITE (iounit, *) "BQB| Collecting electron density from MPI ranks and sending to library..."
540 : END IF
541 : END IF
542 :
543 90 : ALLOCATE (buf(L3:U3))
544 :
545 30 : dest = 0
546 :
547 2222 : DO I1 = L1, U1
548 169922 : DO I2 = L2, U2
549 :
550 : ! cycling through the CPUs, check if the current ray (I1,I2) is local to that CPU
551 167700 : IF (rspace_pw%pw_grid%para%mode /= PW_MODE_LOCAL) THEN
552 503100 : DO ip = 0, num_pe - 1
553 : IF (rspace_pw%pw_grid%para%bo(1, 1, ip, 1) <= I1 - L1 + 1 &
554 : .AND. rspace_pw%pw_grid%para%bo(2, 1, ip, 1) >= I1 - L1 + 1 .AND. &
555 : rspace_pw%pw_grid%para%bo(1, 2, ip, 1) <= I2 - L2 + 1 &
556 503100 : .AND. rspace_pw%pw_grid%para%bo(2, 2, ip, 1) >= I2 - L2 + 1) THEN
557 167700 : source = ip
558 : END IF
559 : END DO
560 : ELSE
561 0 : source = dest
562 : END IF
563 :
564 167700 : IF (source == dest) THEN
565 84300 : IF (my_rank == source) THEN
566 3368006 : buf(:) = rspace_pw%array(I1, I2, :)
567 : END IF
568 : ELSE
569 83400 : IF (my_rank == source) THEN
570 3333806 : buf(:) = rspace_pw%array(I1, I2, :)
571 41700 : CALL para_env%send(buf, dest, tag)
572 : END IF
573 83400 : IF (my_rank == dest) THEN
574 41700 : CALL para_env%recv(buf, source, tag)
575 : END IF
576 : END IF
577 :
578 167700 : IF (my_rank == dest) THEN
579 83850 : ret = libvori_push_rho_zrow(I1 - L1, I2 - L2, U3 - L3 + 1, buf)
580 83850 : IF (ret /= 0) THEN
581 0 : CPABORT("The library returned an error. Aborting.")
582 : END IF
583 : END IF
584 :
585 : ! this double loop generates so many messages that it can overload
586 : ! the message passing system, e.g. on XT3
587 : ! we therefore put a barrier here that limits the amount of message
588 : ! that flies around at any given time.
589 : ! if ever this routine becomes a bottleneck, we should go for a
590 : ! more complicated rewrite
591 169892 : CALL para_env%sync()
592 :
593 : END DO
594 : END DO
595 :
596 30 : DEALLOCATE (buf)
597 :
598 30 : IF (ionode) THEN
599 :
600 15 : gapw = .FALSE.
601 15 : IF (dft_control%qs_control%method_id == do_method_gapw) gapw = .TRUE.
602 :
603 15 : IF (do_voro /= 0) THEN
604 :
605 14 : IF (radius_type == voro_radii_unity) THEN
606 0 : ret = libvori_setRadii_Unity()
607 14 : ELSE IF (radius_type == voro_radii_cov) THEN
608 : ! Use the covalent radii from LIBVORI
609 0 : ret = libvori_setRadii_Covalent()
610 14 : ELSE IF (radius_type == voro_radii_vdw) THEN
611 : ! Use the van der Waals radii from CP2K
612 14 : ret = libvori_setRadii_User(100.0_dp, particles_radius)
613 0 : ELSE IF (radius_type == voro_radii_user) THEN
614 : ! Use the user defined atomic radii
615 0 : IF (natom /= SIZE(user_radii)) THEN
616 : CALL cp_abort(__LOCATION__, &
617 : "Length of keyword VORONOI\USER_RADII does not "// &
618 0 : "match number of atoms in the input coordinate file.")
619 : END IF
620 0 : ret = libvori_setRadii_User(100.0_dp, user_radii)
621 : ELSE
622 0 : CPABORT("No valid radius type was specified for VORONOI")
623 : END IF
624 :
625 14 : IF (voro_sanity) THEN
626 :
627 9 : ret = libvori_sanitycheck(sim_step, sim_time)
628 :
629 9 : IF (ret /= 0) THEN
630 0 : CPABORT("The library returned an error. Aborting.")
631 : END IF
632 :
633 : END IF
634 :
635 14 : ret = libvori_step(sim_step, sim_time)
636 :
637 14 : step_count = step_count + 1
638 :
639 14 : IF (ret /= 0) THEN
640 0 : CPABORT("The library returned an error. Aborting.")
641 : END IF
642 :
643 14 : IF ((step_count > 1) .OR. (.NOT. voro_skip_first)) THEN
644 :
645 42 : ALLOCATE (voro_radii(natom))
646 28 : ALLOCATE (voro_charge(natom))
647 28 : ALLOCATE (voro_volume(natom))
648 42 : ALLOCATE (voro_dipole(natom, 3))
649 42 : ALLOCATE (voro_quadrupole(natom, 9))
650 28 : ALLOCATE (voro_buffer(natom))
651 28 : ALLOCATE (voro_wrapped_pos(natom, 3))
652 28 : ALLOCATE (voro_charge_center(natom, 3))
653 :
654 14 : ret = libvori_get_radius(natom, voro_radii)
655 :
656 14 : ret = libvori_get_charge(natom, voro_charge)
657 :
658 14 : ret = libvori_get_volume(natom, voro_volume)
659 :
660 56 : DO i1 = 1, 3
661 42 : ret = libvori_get_dipole(i1, natom, voro_buffer)
662 215 : voro_dipole(:, i1) = voro_buffer(:)
663 : END DO
664 :
665 140 : DO i1 = 1, 9
666 126 : ret = libvori_get_quadrupole(i1, natom, voro_buffer)
667 617 : voro_quadrupole(:, i1) = voro_buffer(:)
668 : END DO
669 :
670 56 : DO i1 = 1, 3
671 42 : ret = libvori_get_wrapped_pos(i1, natom, voro_buffer)
672 215 : voro_wrapped_pos(:, i1) = voro_buffer(:)
673 : END DO
674 :
675 56 : DO i1 = 1, 3
676 42 : ret = libvori_get_charge_center(i1, natom, voro_buffer)
677 215 : voro_charge_center(:, i1) = voro_buffer(:)
678 : END DO
679 :
680 14 : IF (gapw) THEN
681 2 : CALL get_qs_env(qs_env=qs_env, rho0_mpole=rho0_mpole)
682 2 : nspins = dft_control%nspins
683 6 : DO i1 = 1, natom
684 8 : voro_charge(i1) = voro_charge(i1) - SUM(rho0_mpole%mp_rho(i1)%Q0(1:nspins))
685 4 : fn0 = rho0_mpole%norm_g0l_h(1)/bohr*100._dp
686 16 : voro_dipole(i1, 1:3) = voro_dipole(i1, 1:3) + rho0_mpole%mp_rho(i1)%Qlm_car(2:4)/fn0
687 4 : qa = voro_charge(i1) - particles_c(i1)
688 16 : voro_charge_center(i1, 1:3) = voro_dipole(i1, 1:3)/qa
689 4 : fn1 = rho0_mpole%norm_g0l_h(2)/bohr/bohr*10000._dp
690 4 : voro_quadrupole(i1, 1) = voro_quadrupole(i1, 1) + rho0_mpole%mp_rho(i1)%Qlm_car(5)/fn1
691 4 : voro_quadrupole(i1, 2) = voro_quadrupole(i1, 2) + rho0_mpole%mp_rho(i1)%Qlm_car(6)/fn1
692 4 : voro_quadrupole(i1, 3) = voro_quadrupole(i1, 3) + rho0_mpole%mp_rho(i1)%Qlm_car(7)/fn1
693 4 : voro_quadrupole(i1, 4) = voro_quadrupole(i1, 4) + rho0_mpole%mp_rho(i1)%Qlm_car(6)/fn1
694 4 : voro_quadrupole(i1, 5) = voro_quadrupole(i1, 5) + rho0_mpole%mp_rho(i1)%Qlm_car(8)/fn1
695 4 : voro_quadrupole(i1, 6) = voro_quadrupole(i1, 6) + rho0_mpole%mp_rho(i1)%Qlm_car(9)/fn1
696 4 : voro_quadrupole(i1, 7) = voro_quadrupole(i1, 7) + rho0_mpole%mp_rho(i1)%Qlm_car(7)/fn1
697 4 : voro_quadrupole(i1, 8) = voro_quadrupole(i1, 8) + rho0_mpole%mp_rho(i1)%Qlm_car(9)/fn1
698 6 : voro_quadrupole(i1, 9) = voro_quadrupole(i1, 9) + rho0_mpole%mp_rho(i1)%Qlm_car(10)/fn1
699 : END DO
700 : END IF
701 :
702 14 : IF (unit_voro > 0) THEN
703 14 : WRITE (unit_voro, FMT="(T2,I0)") natom
704 14 : WRITE (unit_voro, FMT="(A,I8,A,F12.4,A)") "# Step ", sim_step, ", Time ", &
705 28 : sim_time*femtoseconds, " fs"
706 14 : WRITE (unit_voro, FMT="(A,9F20.10)") "# Cell ", &
707 14 : cell%hmat(1, 1)*angstrom, cell%hmat(2, 1)*angstrom, cell%hmat(3, 1)*angstrom, &
708 14 : cell%hmat(1, 2)*angstrom, cell%hmat(2, 2)*angstrom, cell%hmat(3, 2)*angstrom, &
709 28 : cell%hmat(1, 3)*angstrom, cell%hmat(2, 3)*angstrom, cell%hmat(3, 3)*angstrom
710 14 : WRITE (unit_voro, FMT="(A,22A20)") "# Atom Z", &
711 14 : "Radius", "Position(X)", "Position(Y)", "Position(Z)", &
712 14 : "Voronoi_Volume", "Z(eff)", "Charge", "Dipole(X)", "Dipole(Y)", "Dipole(Z)", &
713 14 : "ChargeCenter(X)", "ChargeCenter(Y)", "ChargeCenter(Z)", &
714 14 : "Quadrupole(XX)", "Quadrupole(XY)", "Quadrupole(XZ)", &
715 14 : "Quadrupole(YX)", "Quadrupole(YY)", "Quadrupole(YZ)", &
716 28 : "Quadrupole(ZX)", "Quadrupole(ZY)", "Quadrupole(ZZ)"
717 67 : DO i1 = 1, natom
718 : WRITE (unit_voro, FMT="(2I6,22F20.10)") &
719 53 : i1, &
720 53 : particles_z(i1), &
721 53 : voro_radii(i1)/100.0_dp, &
722 212 : particles_r(1:3, i1)*angstrom, &
723 53 : voro_volume(i1)/1000000.0_dp, &
724 53 : particles_c(i1), &
725 53 : voro_charge(i1), &
726 53 : voro_dipole(i1, 1:3), &
727 212 : voro_charge_center(i1, 1:3)/100.0_dp, &
728 120 : voro_quadrupole(i1, 1:9)
729 : END DO
730 : END IF
731 :
732 14 : IF (molprop) THEN
733 : CALL molecular_properties(subsys, cell, sim_step, sim_time, iounit, &
734 : particles_r, particles_c, &
735 8 : voro_charge, voro_charge_center, mp_file_name)
736 : END IF
737 :
738 14 : DEALLOCATE (voro_radii)
739 14 : DEALLOCATE (voro_charge)
740 14 : DEALLOCATE (voro_volume)
741 14 : DEALLOCATE (voro_dipole)
742 14 : DEALLOCATE (voro_quadrupole)
743 14 : DEALLOCATE (voro_buffer)
744 14 : DEALLOCATE (voro_wrapped_pos)
745 14 : DEALLOCATE (voro_charge_center)
746 :
747 : END IF ! not skip_first
748 :
749 14 : IF (iounit > 0) THEN
750 14 : WRITE (iounit, *) "VORONOI| Voronoi integration finished."
751 : END IF
752 :
753 : END IF ! do_voro
754 :
755 15 : IF (do_bqb /= 0) THEN
756 :
757 1 : ret = libvori_processBQBFrame(sim_step, sim_time*femtoseconds)
758 :
759 1 : IF (ret /= 0) THEN
760 0 : CPABORT("The library returned an error. Aborting.")
761 : END IF
762 :
763 1 : IF (do_voro /= 0) THEN
764 0 : IF (iounit > 0) THEN
765 0 : WRITE (iounit, *) "VORONOI| BQB compression finished."
766 : END IF
767 : ELSE
768 1 : IF (iounit > 0) THEN
769 1 : WRITE (iounit, *) "BQB| BQB compression finished."
770 : END IF
771 : END IF
772 :
773 : END IF ! do_bqb
774 :
775 : END IF
776 :
777 30 : IF (ionode) THEN
778 15 : DEALLOCATE (particles_z)
779 15 : DEALLOCATE (particles_c)
780 15 : DEALLOCATE (particles_r)
781 15 : DEALLOCATE (particles_radius)
782 : END IF
783 : END ASSOCIATE
784 :
785 : #if defined(__HAS_IEEE_EXCEPTIONS)
786 : CALL ieee_set_halting_mode(IEEE_ALL, halt)
787 : #endif
788 :
789 30 : CALL timestop(handle)
790 :
791 : #else
792 :
793 : MARK_USED(do_voro)
794 : MARK_USED(do_bqb)
795 : MARK_USED(input_voro)
796 : MARK_USED(input_bqb)
797 : MARK_USED(unit_voro)
798 : MARK_USED(qs_env)
799 : MARK_USED(rspace_pw)
800 :
801 : CALL cp_warn(__LOCATION__, &
802 : "Voronoi integration and BQB output require CP2k to be compiled"// &
803 : " with the -D__LIBVORI preprocessor option.")
804 :
805 : #endif
806 :
807 60 : END SUBROUTINE entry_voronoi_or_bqb
808 :
809 : ! **************************************************************************************************
810 : !> \brief Call libvori's finalize if support is compiled in
811 : ! **************************************************************************************************
812 9671 : SUBROUTINE finalize_libvori()
813 : #if defined(__LIBVORI)
814 : INTEGER(KIND=C_INT) :: ret
815 9671 : ret = libvori_finalize()
816 : #endif
817 9671 : END SUBROUTINE finalize_libvori
818 :
819 : ! **************************************************************************************************
820 : !> \brief ...
821 : !> \param subsys ...
822 : !> \param cell ...
823 : !> \param sim_step ...
824 : !> \param sim_time ...
825 : !> \param iounit ...
826 : !> \param particles_r ...
827 : !> \param particles_c ...
828 : !> \param voro_charge ...
829 : !> \param voro_charge_center ...
830 : !> \param mp_file_name ...
831 : ! **************************************************************************************************
832 8 : SUBROUTINE molecular_properties(subsys, cell, sim_step, sim_time, iounit, &
833 16 : particles_r, particles_c, voro_charge, &
834 8 : voro_charge_center, mp_file_name)
835 : TYPE(qs_subsys_type), POINTER :: subsys
836 : TYPE(cell_type), POINTER :: cell
837 : INTEGER, INTENT(IN) :: sim_step
838 : REAL(KIND=dp), INTENT(IN) :: sim_time
839 : INTEGER, INTENT(IN) :: iounit
840 : REAL(KIND=dp), DIMENSION(:, :) :: particles_r
841 : REAL(KIND=dp), DIMENSION(:) :: particles_c, voro_charge
842 : REAL(KIND=dp), DIMENSION(:, :) :: voro_charge_center
843 : CHARACTER(len=default_path_length) :: mp_file_name
844 :
845 : CHARACTER(len=3) :: fstatus
846 : CHARACTER(len=default_path_length) :: fname
847 : INTEGER :: ia, imol, mk, mpunit, na, na1, na2, &
848 : nmolecule
849 : REAL(KIND=dp) :: cm, ddip
850 : REAL(KIND=dp), DIMENSION(3) :: dipm, posa, posc, ref
851 : TYPE(molecule_kind_type), POINTER :: molecule_kind
852 8 : TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
853 :
854 8 : IF (iounit > 0) THEN
855 8 : WRITE (iounit, *) "VORONOI| Start Calculation of Molecular Properties from Voronoi Integration"
856 : END IF
857 8 : CALL qs_subsys_get(subsys, molecule_set=molecule_set)
858 :
859 8 : IF (INDEX(mp_file_name, "__STD_OUT__") /= 0) THEN
860 8 : mpunit = iounit
861 : ELSE
862 0 : fname = ADJUSTL(mp_file_name)
863 0 : IF (fname(1:2) /= "./") THEN
864 0 : fname = TRIM(fname)//".molprop"
865 : END IF
866 0 : IF (file_exists(fname)) THEN
867 0 : fstatus = "old"
868 : ELSE
869 0 : fstatus = "new"
870 : END IF
871 : CALL open_file(file_name=fname, file_status=fstatus, file_action="write", &
872 0 : file_position="append", unit_number=mpunit)
873 : END IF
874 8 : nmolecule = SIZE(molecule_set)
875 8 : WRITE (mpunit, FMT="(T2,I0)") nmolecule
876 8 : WRITE (mpunit, FMT="(A,I8,A,F12.4,A)") " # Step ", sim_step, ", Time ", &
877 16 : sim_time*femtoseconds, " fs"
878 8 : WRITE (mpunit, FMT="(A,T25,A)") " # Mol Type Charge", &
879 16 : " Dipole[Debye] Total Dipole[Debye]"
880 18 : DO imol = 1, nmolecule
881 10 : molecule_kind => molecule_set(imol)%molecule_kind
882 10 : mk = molecule_kind%kind_number
883 10 : na1 = molecule_set(imol)%first_atom
884 10 : na2 = molecule_set(imol)%last_atom
885 10 : na = na2 - na1 + 1
886 10 : ref(1:3) = 0.0_dp
887 31 : DO ia = na1, na2
888 94 : ref(1:3) = ref(1:3) + pbc(particles_r(1:3, ia), cell)
889 : END DO
890 40 : ref(1:3) = ref(1:3)/REAL(na, KIND=dp)
891 10 : dipm = 0.0_dp
892 31 : DO ia = na1, na2
893 84 : posa(1:3) = particles_r(1:3, ia) - ref(1:3)
894 84 : posa(1:3) = pbc(posa, cell)
895 84 : posc(1:3) = posa(1:3) + bohr*voro_charge_center(ia, 1:3)/100.0_dp
896 84 : posc(1:3) = pbc(posc, cell)
897 21 : cm = -particles_c(ia) + voro_charge(ia)
898 94 : dipm(1:3) = dipm(1:3) + posa(1:3)*particles_c(ia) + posc(1:3)*cm
899 : END DO
900 40 : dipm(1:3) = dipm(1:3)*debye
901 40 : ddip = SQRT(SUM(dipm**2))
902 31 : cm = SUM(voro_charge(na1:na2))
903 18 : WRITE (mpunit, FMT="(I8,I6,F12.4,T25,3F12.4,8X,F12.4)") imol, mk, cm, dipm(1:3), ddip
904 : END DO
905 8 : IF (mpunit /= iounit) THEN
906 0 : CALL close_file(mpunit)
907 : END IF
908 :
909 8 : END SUBROUTINE molecular_properties
910 :
911 : END MODULE voronoi_interface
912 :
|
