Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Module performing a mdoe selective vibrational analysis
10 : !> \note
11 : !> Numerical accuracy for parallel runs:
12 : !> Each replica starts the SCF run from the one optimized
13 : !> in a previous run. It may happen then energies and derivatives
14 : !> of a serial run and a parallel run could be slightly different
15 : !> 'cause of a different starting density matrix.
16 : !> Exact results are obtained using:
17 : !> EXTRAPOLATION USE_GUESS in QS section (Teo 08.2006)
18 : !> \author Florian Schiffmann 08.2006
19 : ! **************************************************************************************************
20 : MODULE mode_selective
21 : USE cell_types, ONLY: cell_type
22 : USE cp_files, ONLY: close_file,&
23 : open_file
24 : USE cp_log_handling, ONLY: cp_get_default_logger,&
25 : cp_logger_get_default_io_unit,&
26 : cp_logger_type
27 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
28 : cp_print_key_unit_nr
29 : USE cp_result_methods, ONLY: get_results
30 : USE global_types, ONLY: global_environment_type
31 : USE input_constants, ONLY: ms_guess_atomic,&
32 : ms_guess_bfgs,&
33 : ms_guess_molden,&
34 : ms_guess_restart,&
35 : ms_guess_restart_vec
36 : USE input_section_types, ONLY: section_vals_get,&
37 : section_vals_get_subs_vals,&
38 : section_vals_type,&
39 : section_vals_val_get
40 : USE kinds, ONLY: default_path_length,&
41 : default_string_length,&
42 : dp,&
43 : max_line_length
44 : USE mathlib, ONLY: diamat_all
45 : USE message_passing, ONLY: mp_para_env_type
46 : USE molden_utils, ONLY: write_vibrations_molden
47 : USE particle_types, ONLY: particle_type
48 : USE physcon, ONLY: bohr,&
49 : debye,&
50 : massunit,&
51 : vibfac
52 : USE replica_methods, ONLY: rep_env_calc_e_f
53 : USE replica_types, ONLY: replica_env_type
54 : USE util, ONLY: sort
55 : #include "./base/base_uses.f90"
56 :
57 : IMPLICIT NONE
58 :
59 : PRIVATE
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mode_selective'
61 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
62 :
63 : TYPE ms_vib_type
64 : INTEGER :: mat_size = -1
65 : INTEGER :: select_id = -1
66 : INTEGER, DIMENSION(:), POINTER :: inv_atoms => NULL()
67 : REAL(KIND=dp) :: eps(2) = 0.0_dp
68 : REAL(KIND=dp) :: sel_freq = 0.0_dp
69 : REAL(KIND=dp) :: low_freq = 0.0_dp
70 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: b_vec => NULL()
71 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: delta_vec => NULL()
72 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ms_force => NULL()
73 : REAL(KIND=dp), DIMENSION(:), POINTER :: eig_bfgs => NULL()
74 : REAL(KIND=dp), DIMENSION(:), POINTER :: f_range => NULL()
75 : REAL(KIND=dp), DIMENSION(:), POINTER :: inv_range => NULL()
76 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_b => NULL()
77 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_r => NULL()
78 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: b_mat => NULL()
79 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_deriv => NULL()
80 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hes_bfgs => NULL()
81 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_mat => NULL()
82 : INTEGER :: initial_guess = -1
83 : END TYPE ms_vib_type
84 :
85 : PUBLIC :: ms_vb_anal
86 :
87 : CONTAINS
88 : ! **************************************************************************************************
89 : !> \brief Module performing a vibrational analysis
90 : !> \param input ...
91 : !> \param rep_env ...
92 : !> \param para_env ...
93 : !> \param globenv ...
94 : !> \param particles ...
95 : !> \param nrep ...
96 : !> \param calc_intens ...
97 : !> \param dx ...
98 : !> \param output_unit ...
99 : !> \param logger ...
100 : !> \param cell simulation cell
101 : !> \author Teodoro Laino 08.2006
102 : ! **************************************************************************************************
103 24 : SUBROUTINE ms_vb_anal(input, rep_env, para_env, globenv, particles, &
104 : nrep, calc_intens, dx, output_unit, logger, cell)
105 : TYPE(section_vals_type), POINTER :: input
106 : TYPE(replica_env_type), POINTER :: rep_env
107 : TYPE(mp_para_env_type), POINTER :: para_env
108 : TYPE(global_environment_type), POINTER :: globenv
109 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
110 : INTEGER :: nrep
111 : LOGICAL :: calc_intens
112 : REAL(KIND=dp) :: dx
113 : INTEGER :: output_unit
114 : TYPE(cp_logger_type), POINTER :: logger
115 : TYPE(cell_type), POINTER :: cell
116 :
117 : CHARACTER(len=*), PARAMETER :: routineN = 'ms_vb_anal'
118 :
119 : CHARACTER(LEN=default_string_length) :: description
120 : INTEGER :: handle, i, ip1, j, natoms, ncoord
121 : LOGICAL :: converged
122 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mass, pos0
123 24 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: tmp_deriv
124 24 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: tmp_dip
125 : TYPE(ms_vib_type) :: ms_vib
126 :
127 24 : CALL timeset(routineN, handle)
128 24 : converged = .FALSE.
129 24 : natoms = SIZE(particles)
130 24 : ncoord = 3*natoms
131 96 : ALLOCATE (mass(3*natoms))
132 886 : DO i = 1, natoms
133 3472 : DO j = 1, 3
134 2586 : mass((i - 1)*3 + j) = particles(i)%atomic_kind%mass
135 3448 : mass((i - 1)*3 + j) = SQRT(mass((i - 1)*3 + j))
136 : END DO
137 : END DO
138 : ! Allocate working arrays
139 96 : ALLOCATE (ms_vib%delta_vec(ncoord, nrep))
140 72 : ALLOCATE (ms_vib%b_vec(ncoord, nrep))
141 72 : ALLOCATE (ms_vib%step_r(nrep))
142 48 : ALLOCATE (ms_vib%step_b(nrep))
143 24 : IF (calc_intens) THEN
144 20 : description = '[DIPOLE]'
145 80 : ALLOCATE (tmp_dip(nrep, 3, 2))
146 60 : ALLOCATE (ms_vib%dip_deriv(3, nrep))
147 : END IF
148 : CALL MS_initial_moves(para_env, nrep, input, globenv, ms_vib, &
149 : particles, &
150 : mass, &
151 : dx, &
152 24 : calc_intens, logger)
153 24 : ncoord = 3*natoms
154 72 : ALLOCATE (pos0(ncoord))
155 96 : ALLOCATE (ms_vib%ms_force(ncoord, nrep))
156 886 : DO i = 1, natoms
157 3472 : DO j = 1, 3
158 3448 : pos0((i - 1)*3 + j) = particles((i))%r(j)
159 : END DO
160 : END DO
161 162 : ncoord = 3*natoms
162 : DO
163 23178 : ms_vib%ms_force = HUGE(0.0_dp)
164 336 : DO i = 1, nrep
165 23178 : DO j = 1, ncoord
166 23016 : rep_env%r(j, i) = pos0(j) + ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
167 : END DO
168 : END DO
169 162 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
170 :
171 336 : DO i = 1, nrep
172 174 : IF (calc_intens) THEN
173 : CALL get_results(results=rep_env%results(i)%results, &
174 : description=description, &
175 150 : n_rep=ip1)
176 : CALL get_results(results=rep_env%results(i)%results, &
177 : description=description, &
178 : values=tmp_dip(i, :, 1), &
179 150 : nval=ip1)
180 : END IF
181 23178 : DO j = 1, ncoord
182 23016 : ms_vib%ms_force(j, i) = rep_env%f(j, i)
183 : END DO
184 : END DO
185 336 : DO i = 1, nrep
186 23178 : DO j = 1, ncoord
187 23016 : rep_env%r(j, i) = pos0(j) - ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
188 : END DO
189 : END DO
190 162 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
191 162 : IF (calc_intens) THEN
192 300 : DO i = 1, nrep
193 : CALL get_results(results=rep_env%results(i)%results, &
194 : description=description, &
195 150 : n_rep=ip1)
196 : CALL get_results(results=rep_env%results(i)%results, &
197 : description=description, &
198 : values=tmp_dip(i, :, 2), &
199 150 : nval=ip1)
200 900 : ms_vib%dip_deriv(:, ms_vib%mat_size + i) = (tmp_dip(i, :, 1) - tmp_dip(i, :, 2))/(2*ms_vib%step_b(i))
201 : END DO
202 : END IF
203 :
204 : CALL evaluate_H_update_b(rep_env, ms_vib, input, nrep, &
205 : particles, &
206 : mass, &
207 : converged, &
208 : dx, calc_intens, &
209 162 : output_unit, logger, cell)
210 162 : IF (converged) EXIT
211 162 : IF (calc_intens) THEN
212 390 : ALLOCATE (tmp_deriv(3, ms_vib%mat_size))
213 3730 : tmp_deriv = ms_vib%dip_deriv
214 130 : DEALLOCATE (ms_vib%dip_deriv)
215 390 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
216 3730 : ms_vib%dip_deriv(:, 1:ms_vib%mat_size) = tmp_deriv(:, 1:ms_vib%mat_size)
217 130 : DEALLOCATE (tmp_deriv)
218 : END IF
219 : END DO
220 24 : DEALLOCATE (ms_vib%ms_force)
221 24 : DEALLOCATE (pos0)
222 24 : DEALLOCATE (ms_vib%step_r)
223 24 : DEALLOCATE (ms_vib%step_b)
224 24 : DEALLOCATE (ms_vib%b_vec)
225 24 : DEALLOCATE (ms_vib%delta_vec)
226 24 : DEALLOCATE (mass)
227 24 : DEALLOCATE (ms_vib%b_mat)
228 24 : DEALLOCATE (ms_vib%s_mat)
229 24 : IF (ms_vib%select_id == 3) THEN
230 8 : DEALLOCATE (ms_vib%inv_atoms)
231 : END IF
232 24 : IF (ASSOCIATED(ms_vib%eig_bfgs)) THEN
233 0 : DEALLOCATE (ms_vib%eig_bfgs)
234 : END IF
235 24 : IF (ASSOCIATED(ms_vib%hes_bfgs)) THEN
236 0 : DEALLOCATE (ms_vib%hes_bfgs)
237 : END IF
238 24 : IF (calc_intens) THEN
239 20 : DEALLOCATE (ms_vib%dip_deriv)
240 20 : DEALLOCATE (tmp_dip)
241 : END IF
242 24 : CALL timestop(handle)
243 72 : END SUBROUTINE ms_vb_anal
244 : ! **************************************************************************************************
245 : !> \brief Generates the first displacement vector for a mode selctive vibrational
246 : !> analysis. At the moment this is a random number for selected atoms
247 : !> \param para_env ...
248 : !> \param nrep ...
249 : !> \param input ...
250 : !> \param globenv ...
251 : !> \param ms_vib ...
252 : !> \param particles ...
253 : !> \param mass ...
254 : !> \param dx ...
255 : !> \param calc_intens ...
256 : !> \param logger ...
257 : !> \author Florian Schiffmann 11.2007
258 : ! **************************************************************************************************
259 24 : SUBROUTINE MS_initial_moves(para_env, nrep, input, globenv, ms_vib, particles, &
260 24 : mass, dx, &
261 : calc_intens, logger)
262 : TYPE(mp_para_env_type), POINTER :: para_env
263 : INTEGER :: nrep
264 : TYPE(section_vals_type), POINTER :: input
265 : TYPE(global_environment_type), POINTER :: globenv
266 : TYPE(ms_vib_type) :: ms_vib
267 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
268 : REAL(Kind=dp), DIMENSION(:) :: mass
269 : REAL(KIND=dp) :: dx
270 : LOGICAL :: calc_intens
271 : TYPE(cp_logger_type), POINTER :: logger
272 :
273 : CHARACTER(len=*), PARAMETER :: routineN = 'MS_initial_moves'
274 :
275 : INTEGER :: guess, handle, i, j, jj, k, m, &
276 : n_rep_val, natoms, ncoord
277 24 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map_atoms
278 24 : INTEGER, DIMENSION(:), POINTER :: tmplist
279 : LOGICAL :: do_involved_atoms, ionode
280 : REAL(KIND=dp) :: my_val, norm
281 : TYPE(section_vals_type), POINTER :: involved_at_section, ms_vib_section
282 :
283 24 : CALL timeset(routineN, handle)
284 24 : NULLIFY (ms_vib%eig_bfgs, ms_vib%f_range, ms_vib%hes_bfgs, ms_vib%inv_range)
285 24 : ms_vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
286 24 : CALL section_vals_val_get(ms_vib_section, "INITIAL_GUESS", i_val=guess)
287 24 : CALL section_vals_val_get(ms_vib_section, "EPS_MAX_VAL", r_val=ms_vib%eps(1))
288 24 : CALL section_vals_val_get(ms_vib_section, "EPS_NORM", r_val=ms_vib%eps(2))
289 24 : CALL section_vals_val_get(ms_vib_section, "RANGE", n_rep_val=n_rep_val)
290 24 : ms_vib%select_id = 0
291 24 : IF (n_rep_val /= 0) THEN
292 2 : CALL section_vals_val_get(ms_vib_section, "RANGE", r_vals=ms_vib%f_range)
293 2 : IF (ms_vib%f_range(1) > ms_vib%f_range(2)) THEN
294 0 : my_val = ms_vib%f_range(2)
295 0 : ms_vib%f_range(2) = ms_vib%f_range(1)
296 0 : ms_vib%f_range(1) = my_val
297 : END IF
298 2 : ms_vib%select_id = 2
299 : END IF
300 24 : CALL section_vals_val_get(ms_vib_section, "FREQUENCY", r_val=ms_vib%sel_freq)
301 24 : CALL section_vals_val_get(ms_vib_section, "LOWEST_FREQUENCY", r_val=ms_vib%low_freq)
302 24 : IF (ms_vib%sel_freq > 0._dp) ms_vib%select_id = 1
303 24 : involved_at_section => section_vals_get_subs_vals(ms_vib_section, "INVOLVED_ATOMS")
304 24 : CALL section_vals_get(involved_at_section, explicit=do_involved_atoms)
305 24 : IF (do_involved_atoms) THEN
306 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", n_rep_val=n_rep_val)
307 8 : jj = 0
308 16 : DO k = 1, n_rep_val
309 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=k, i_vals=tmplist)
310 32 : DO j = 1, SIZE(tmplist)
311 24 : jj = jj + 1
312 : END DO
313 : END DO
314 8 : IF (jj >= 1) THEN
315 8 : natoms = jj
316 24 : ALLOCATE (ms_vib%inv_atoms(natoms))
317 8 : jj = 0
318 16 : DO m = 1, n_rep_val
319 8 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=m, i_vals=tmplist)
320 32 : DO j = 1, SIZE(tmplist)
321 24 : ms_vib%inv_atoms(j) = tmplist(j)
322 : END DO
323 : END DO
324 8 : ms_vib%select_id = 3
325 : END IF
326 8 : CALL section_vals_val_get(involved_at_section, "RANGE", n_rep_val=n_rep_val)
327 8 : IF (n_rep_val /= 0) THEN
328 0 : CALL section_vals_val_get(involved_at_section, "RANGE", r_vals=ms_vib%inv_range)
329 0 : IF (ms_vib%inv_range(1) > ms_vib%inv_range(2)) THEN
330 0 : ms_vib%inv_range(2) = my_val
331 0 : ms_vib%inv_range(2) = ms_vib%inv_range(1)
332 0 : ms_vib%inv_range(1) = my_val
333 : END IF
334 : END IF
335 : END IF
336 24 : IF (ms_vib%select_id == 0) THEN
337 0 : CPABORT("no frequency, range or involved atoms specified ")
338 : END IF
339 24 : ionode = para_env%is_source()
340 12 : SELECT CASE (guess)
341 : CASE (ms_guess_atomic)
342 12 : ms_vib%initial_guess = 1
343 12 : CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
344 12 : jj = 0
345 22 : DO k = 1, n_rep_val
346 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
347 42 : DO j = 1, SIZE(tmplist)
348 30 : jj = jj + 1
349 : END DO
350 : END DO
351 12 : IF (jj < 1) THEN
352 2 : natoms = SIZE(particles)
353 6 : ALLOCATE (map_atoms(natoms))
354 14 : DO j = 1, natoms
355 14 : map_atoms(j) = j
356 : END DO
357 : ELSE
358 10 : natoms = jj
359 30 : ALLOCATE (map_atoms(natoms))
360 10 : jj = 0
361 20 : DO m = 1, n_rep_val
362 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
363 40 : DO j = 1, SIZE(tmplist)
364 30 : map_atoms(j) = tmplist(j)
365 : END DO
366 : END DO
367 : END IF
368 :
369 : ! apply random displacement along the mass weighted nuclear cartesian coordinates
370 526 : ms_vib%b_vec = 0._dp
371 526 : ms_vib%delta_vec = 0._dp
372 12 : jj = 0
373 :
374 28 : DO i = 1, nrep
375 56 : DO j = 1, natoms
376 176 : DO k = 1, 3
377 120 : jj = (map_atoms(j) - 1)*3 + k
378 160 : ms_vib%b_vec(jj, i) = ABS(globenv%gaussian_rng_stream%next())
379 : END DO
380 : END DO
381 514 : norm = NORM2(ms_vib%b_vec(:, i))
382 526 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
383 : END DO
384 :
385 12 : IF (nrep > 1) THEN
386 44 : DO k = 1, 10
387 124 : DO j = 1, nrep
388 280 : DO i = 1, nrep
389 240 : IF (i /= j) THEN
390 : ms_vib%b_vec(:, j) = &
391 1520 : ms_vib%b_vec(:, j) - DOT_PRODUCT(ms_vib%b_vec(:, j), ms_vib%b_vec(:, i))*ms_vib%b_vec(:, i)
392 : ms_vib%b_vec(:, j) = &
393 1520 : ms_vib%b_vec(:, j)/NORM2(ms_vib%b_vec(:, j))
394 : END IF
395 : END DO
396 : END DO
397 : END DO
398 : END IF
399 :
400 12 : ms_vib%mat_size = 0
401 474 : DO i = 1, SIZE(ms_vib%b_vec, 1)
402 972 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
403 : END DO
404 : CASE (ms_guess_bfgs)
405 :
406 4 : ms_vib%initial_guess = 2
407 4 : CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
408 4 : ms_vib%mat_size = 0
409 :
410 : CASE (ms_guess_restart_vec)
411 :
412 4 : ms_vib%initial_guess = 3
413 : ncoord = 3*SIZE(particles)
414 4 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
415 :
416 4 : ms_vib%mat_size = 0
417 : CASE (ms_guess_restart)
418 0 : ms_vib%initial_guess = 4
419 : ncoord = 3*SIZE(particles)
420 0 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
421 :
422 : CASE (ms_guess_molden)
423 4 : ms_vib%initial_guess = 5
424 4 : ncoord = 3*SIZE(particles)
425 4 : CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
426 28 : ms_vib%mat_size = 0
427 : END SELECT
428 5324 : CALL para_env%bcast(ms_vib%b_vec)
429 5324 : CALL para_env%bcast(ms_vib%delta_vec)
430 52 : DO i = 1, nrep
431 2650 : ms_vib%step_r(i) = dx/NORM2(ms_vib%delta_vec(:, i))
432 2674 : ms_vib%step_b(i) = NORM2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
433 : END DO
434 24 : CALL timestop(handle)
435 :
436 48 : END SUBROUTINE MS_initial_moves
437 :
438 : ! **************************************************************************************************
439 : !> \brief ...
440 : !> \param ms_vib_section ...
441 : !> \param ms_vib ...
442 : !> \param particles ...
443 : !> \param mass ...
444 : !> \param para_env ...
445 : !> \param nrep ...
446 : !> \author Florian Schiffmann 11.2007
447 : ! **************************************************************************************************
448 4 : SUBROUTINE bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
449 :
450 : TYPE(section_vals_type), POINTER :: ms_vib_section
451 : TYPE(ms_vib_type) :: ms_vib
452 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
453 : REAL(Kind=dp), DIMENSION(:) :: mass
454 : TYPE(mp_para_env_type), POINTER :: para_env
455 : INTEGER :: nrep
456 :
457 : CHARACTER(LEN=default_path_length) :: hes_filename
458 : INTEGER :: hesunit, i, istat, j, jj, k, natoms, &
459 : ncoord, output_unit, stat
460 4 : INTEGER, DIMENSION(:), POINTER :: tmplist
461 : REAL(KIND=dp) :: my_val, norm
462 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
463 : TYPE(cp_logger_type), POINTER :: logger
464 :
465 8 : logger => cp_get_default_logger()
466 4 : output_unit = cp_logger_get_default_io_unit(logger)
467 :
468 4 : natoms = SIZE(particles)
469 4 : ncoord = 3*natoms
470 :
471 16 : ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
472 12 : ALLOCATE (ms_vib%eig_bfgs(ncoord))
473 :
474 4 : IF (para_env%is_source()) THEN
475 2 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
476 2 : IF (hes_filename == "") hes_filename = "HESSIAN"
477 : CALL open_file(file_name=hes_filename, file_status="OLD", &
478 2 : file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
479 6 : ALLOCATE (tmp(ncoord))
480 6 : ALLOCATE (tmplist(ncoord))
481 :
482 : ! should use the cp_fm_read_unformatted...
483 2 : istat = 0
484 356 : DO i = 1, ncoord
485 354 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%hes_bfgs(:, i)
486 356 : istat = istat + stat
487 : END DO
488 2 : CALL close_file(hesunit)
489 2 : IF (output_unit > 0) THEN
490 2 : IF (istat /= 0) THEN
491 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading HESSIAN **"
492 : ELSE
493 : WRITE (output_unit, FMT="(/,T2,A)") &
494 2 : "*** Initial Hessian has been read successfully ***"
495 : END IF
496 : END IF
497 356 : DO i = 1, ncoord
498 63014 : DO j = 1, ncoord
499 63012 : ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
500 : END DO
501 : END DO
502 :
503 2 : CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
504 2 : tmp(:) = 0._dp
505 2 : IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
506 2 : IF (ms_vib%select_id == 2) my_val = (((ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp)/vibfac)**2/massunit
507 2 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
508 178 : DO i = 1, ncoord
509 178 : tmp(i) = ABS(my_val - ms_vib%eig_bfgs(i))
510 : END DO
511 1 : ELSE IF (ms_vib%select_id == 3) THEN
512 178 : DO i = 1, ncoord
513 531 : DO j = 1, SIZE(ms_vib%inv_atoms)
514 1593 : DO k = 1, 3
515 1062 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
516 1416 : tmp(i) = tmp(i) + SQRT(ms_vib%hes_bfgs(jj, i)**2)
517 : END DO
518 : END DO
519 178 : IF ((SIGN(1._dp, ms_vib%eig_bfgs(i))*SQRT(ABS(ms_vib%eig_bfgs(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
520 : END DO
521 178 : tmp(:) = -tmp(:)
522 : END IF
523 2 : CALL sort(tmp, ncoord, tmplist)
524 4 : DO i = 1, nrep
525 356 : ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
526 356 : norm = NORM2(ms_vib%b_vec(:, i))
527 358 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
528 : END DO
529 356 : DO i = 1, SIZE(ms_vib%b_vec, 1)
530 710 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
531 : END DO
532 2 : DEALLOCATE (tmp)
533 2 : DEALLOCATE (tmplist)
534 : END IF
535 :
536 1428 : CALL para_env%bcast(ms_vib%b_vec)
537 1428 : CALL para_env%bcast(ms_vib%delta_vec)
538 :
539 4 : DEALLOCATE (ms_vib%hes_bfgs)
540 4 : DEALLOCATE (ms_vib%eig_bfgs)
541 4 : ms_vib%mat_size = 0
542 :
543 4 : END SUBROUTINE bfgs_guess
544 :
545 : ! **************************************************************************************************
546 : !> \brief ...
547 : !> \param ms_vib_section ...
548 : !> \param para_env ...
549 : !> \param ms_vib ...
550 : !> \param mass ...
551 : !> \param ionode ...
552 : !> \param particles ...
553 : !> \param nrep ...
554 : !> \param calc_intens ...
555 : !> \author Florian Schiffmann 11.2007
556 : ! **************************************************************************************************
557 4 : SUBROUTINE rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
558 :
559 : TYPE(section_vals_type), POINTER :: ms_vib_section
560 : TYPE(mp_para_env_type), POINTER :: para_env
561 : TYPE(ms_vib_type) :: ms_vib
562 : REAL(Kind=dp), DIMENSION(:) :: mass
563 : LOGICAL :: ionode
564 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
565 : INTEGER :: nrep
566 : LOGICAL :: calc_intens
567 :
568 : CHARACTER(LEN=default_path_length) :: ms_filename
569 : INTEGER :: hesunit, i, j, mat, natoms, ncoord, &
570 : output_unit, stat, statint
571 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
572 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
573 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H
574 : TYPE(cp_logger_type), POINTER :: logger
575 :
576 4 : logger => cp_get_default_logger()
577 4 : output_unit = cp_logger_get_default_io_unit(logger)
578 :
579 4 : natoms = SIZE(particles)
580 4 : ncoord = 3*natoms
581 4 : IF (calc_intens) THEN
582 4 : DEALLOCATE (ms_vib%dip_deriv)
583 : END IF
584 :
585 4 : IF (ionode) THEN
586 :
587 2 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
588 2 : IF (ms_filename == "") ms_filename = "MS_RESTART"
589 : CALL open_file(file_name=ms_filename, &
590 : file_status="UNKNOWN", &
591 : file_form="UNFORMATTED", &
592 : file_action="READ", &
593 2 : unit_number=hesunit)
594 2 : READ (UNIT=hesunit, IOSTAT=stat) mat
595 2 : CPASSERT(stat == 0)
596 2 : ms_vib%mat_size = mat
597 : END IF
598 4 : CALL para_env%bcast(ms_vib%mat_size)
599 16 : ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
600 12 : ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
601 4 : IF (calc_intens) THEN
602 12 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
603 : END IF
604 4 : IF (ionode) THEN
605 2 : statint = 0
606 3384 : READ (UNIT=hesunit) ms_vib%b_mat
607 3384 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%s_mat
608 2 : IF (stat /= 0 .AND. output_unit > 0) THEN
609 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART **"
610 : END IF
611 2 : IF (calc_intens) THEN
612 2 : READ (UNIT=hesunit, IOSTAT=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
613 2 : IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART,", &
614 0 : "intensities are requested but not present in restart file **"
615 : END IF
616 2 : CALL close_file(hesunit)
617 2 : IF (stat == 0 .AND. statint == 0 .AND. output_unit > 0) THEN
618 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** MS_RESTART has been read successfully ***"
619 : END IF
620 : END IF
621 13532 : CALL para_env%bcast(ms_vib%b_mat)
622 13532 : CALL para_env%bcast(ms_vib%s_mat)
623 340 : IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
624 16 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
625 12 : ALLOCATE (eigenval(ms_vib%mat_size))
626 12 : ALLOCATE (ind(ms_vib%mat_size))
627 :
628 : CALL dgemm('T', 'N', ms_vib%mat_size, ms_vib%mat_size, SIZE(ms_vib%s_mat, 1), 1._dp, ms_vib%b_mat, SIZE(ms_vib%b_mat, 1), &
629 4 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
630 4 : CALL diamat_all(approx_H, eigenval)
631 :
632 4 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, ms_vib%b_vec)
633 4 : IF (ms_vib%initial_guess /= 4) THEN
634 :
635 716 : ms_vib%b_vec = 0._dp
636 8 : DO i = 1, nrep
637 42 : DO j = 1, ms_vib%mat_size
638 6768 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_H(j, ind(i))*ms_vib%b_mat(:, j)
639 : END DO
640 1424 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/NORM2(ms_vib%b_vec(:, i))
641 : END DO
642 :
643 4 : DEALLOCATE (ms_vib%s_mat)
644 4 : DEALLOCATE (ms_vib%b_mat)
645 4 : IF (calc_intens) THEN
646 4 : DEALLOCATE (ms_vib%dip_deriv)
647 12 : ALLOCATE (ms_vib%dip_deriv(3, nrep))
648 : END IF
649 : END IF
650 4 : DEALLOCATE (approx_H)
651 4 : DEALLOCATE (eigenval)
652 4 : DEALLOCATE (ind)
653 8 : DO i = 1, nrep
654 716 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
655 : END DO
656 :
657 4 : END SUBROUTINE rest_guess
658 :
659 : ! **************************************************************************************************
660 : !> \brief ...
661 : !> \param ms_vib_section ...
662 : !> \param input ...
663 : !> \param para_env ...
664 : !> \param ms_vib ...
665 : !> \param mass ...
666 : !> \param ncoord ...
667 : !> \param nrep ...
668 : !> \param logger ...
669 : !> \author Florian Schiffmann 11.2007
670 : ! **************************************************************************************************
671 4 : SUBROUTINE molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
672 : TYPE(section_vals_type), POINTER :: ms_vib_section, input
673 : TYPE(mp_para_env_type), POINTER :: para_env
674 : TYPE(ms_vib_type) :: ms_vib
675 : REAL(Kind=dp), DIMENSION(:) :: mass
676 : INTEGER :: ncoord, nrep
677 : TYPE(cp_logger_type), POINTER :: logger
678 :
679 : CHARACTER(LEN=2) :: at_name
680 : CHARACTER(LEN=default_path_length) :: ms_filename
681 : CHARACTER(LEN=max_line_length) :: info
682 : INTEGER :: i, istat, iw, j, jj, k, nvibs, &
683 : output_molden, output_unit, stat
684 4 : INTEGER, DIMENSION(:), POINTER :: tmplist
685 : LOGICAL :: reading_vib
686 : REAL(KIND=dp) :: my_val, norm
687 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
688 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
689 8 : REAL(KIND=dp), DIMENSION(3, ncoord/3) :: pos
690 :
691 8 : output_unit = cp_logger_get_default_io_unit(logger)
692 :
693 4 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
694 4 : IF (ms_filename == "") output_molden = &
695 : cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
696 : extension=".mol", file_status='UNKNOWN', &
697 0 : file_action="READ")
698 4 : IF (para_env%is_source()) THEN
699 :
700 2 : IF (ms_filename == "") THEN
701 0 : iw = output_molden
702 : ELSE
703 : CALL open_file(file_name=TRIM(ms_filename), &
704 : file_status="UNKNOWN", &
705 : file_form="FORMATTED", &
706 : file_action="READ", &
707 2 : unit_number=iw)
708 : END IF
709 2 : info = ""
710 2 : READ (iw, *) info
711 2 : READ (iw, *) info
712 2 : istat = 0
713 2 : nvibs = 0
714 2 : reading_vib = .FALSE.
715 140 : DO
716 142 : READ (iw, *, IOSTAT=stat) info
717 142 : istat = istat + stat
718 142 : IF (TRIM(ADJUSTL(info)) == "[FR-COORD]") EXIT
719 :
720 140 : CPASSERT(stat == 0)
721 :
722 140 : IF (reading_vib) nvibs = nvibs + 1
723 140 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") reading_vib = .TRUE.
724 : END DO
725 2 : REWIND (iw)
726 2 : istat = 0
727 2 : READ (iw, *, IOSTAT=stat) info
728 2 : istat = istat + stat
729 2 : READ (iw, *, IOSTAT=stat) info
730 2 : istat = istat + stat
731 : ! Skip [Atoms] section
732 118 : DO
733 120 : READ (iw, *, IOSTAT=stat) info
734 120 : istat = istat + stat
735 120 : CPASSERT(stat == 0)
736 120 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") EXIT
737 : END DO
738 : ! Read frequencies and modes
739 6 : ALLOCATE (freq(nvibs))
740 8 : ALLOCATE (modes(ncoord, nvibs))
741 :
742 22 : DO i = 1, nvibs
743 20 : READ (iw, *, IOSTAT=stat) freq(i)
744 22 : istat = istat + stat
745 : END DO
746 2 : READ (iw, *) info
747 120 : DO i = 1, ncoord/3
748 118 : READ (iw, *, IOSTAT=stat) at_name, pos(:, i)
749 120 : istat = istat + stat
750 : END DO
751 2 : READ (iw, *) info
752 22 : DO i = 1, nvibs
753 20 : READ (iw, *) info
754 20 : istat = istat + stat
755 1202 : DO j = 1, ncoord/3
756 1180 : k = (j - 1)*3 + 1
757 1180 : READ (iw, *, IOSTAT=stat) modes(k:k + 2, i)
758 1200 : istat = istat + stat
759 : END DO
760 : END DO
761 2 : IF (ms_filename /= "") CALL close_file(iw)
762 2 : IF (output_unit > 0) THEN
763 2 : IF (istat /= 0) THEN
764 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MOLDEN file **"
765 : ELSE
766 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
767 : END IF
768 : END IF
769 : !!!!!!! select modes !!!!!!
770 4 : ALLOCATE (tmp(nvibs))
771 2 : tmp(:) = 0.0_dp
772 6 : ALLOCATE (tmplist(nvibs))
773 2 : IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
774 2 : IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
775 2 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
776 11 : DO i = 1, nvibs
777 11 : tmp(i) = ABS(my_val - freq(i))
778 : END DO
779 1 : ELSE IF (ms_vib%select_id == 3) THEN
780 11 : DO i = 1, nvibs
781 30 : DO j = 1, SIZE(ms_vib%inv_atoms)
782 90 : DO k = 1, 3
783 60 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
784 80 : tmp(i) = tmp(i) + SQRT(modes(jj, i)**2)
785 : END DO
786 : END DO
787 11 : IF (freq(i) <= 400._dp) tmp(i) = 0._dp
788 : END DO
789 11 : tmp(:) = -tmp(:)
790 : END IF
791 2 : CALL sort(tmp, nvibs, tmplist)
792 4 : DO i = 1, nrep
793 356 : ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
794 356 : norm = NORM2(ms_vib%b_vec(:, i))
795 358 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
796 : END DO
797 4 : DO i = 1, nrep
798 358 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
799 : END DO
800 :
801 2 : DEALLOCATE (freq)
802 2 : DEALLOCATE (modes)
803 2 : DEALLOCATE (tmp)
804 2 : DEALLOCATE (tmplist)
805 :
806 : END IF
807 1428 : CALL para_env%bcast(ms_vib%b_vec)
808 1428 : CALL para_env%bcast(ms_vib%delta_vec)
809 :
810 4 : IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
811 0 : "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
812 4 : END SUBROUTINE molden_guess
813 :
814 : ! **************************************************************************************************
815 : !> \brief Davidson algorithm for to generate a approximate Hessian for mode
816 : !> selective vibrational analysis
817 : !> \param rep_env ...
818 : !> \param ms_vib ...
819 : !> \param input ...
820 : !> \param nrep ...
821 : !> \param particles ...
822 : !> \param mass ...
823 : !> \param converged ...
824 : !> \param dx ...
825 : !> \param calc_intens ...
826 : !> \param output_unit_ms ...
827 : !> \param logger ...
828 : !> \param cell simulation cell
829 : !> \author Florian Schiffmann 11.2007
830 : ! **************************************************************************************************
831 162 : SUBROUTINE evaluate_H_update_b(rep_env, ms_vib, input, nrep, &
832 : particles, &
833 162 : mass, &
834 : converged, dx, &
835 : calc_intens, output_unit_ms, logger, cell)
836 : TYPE(replica_env_type), POINTER :: rep_env
837 : TYPE(ms_vib_type) :: ms_vib
838 : TYPE(section_vals_type), POINTER :: input
839 : INTEGER :: nrep
840 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
841 : REAL(Kind=dp), DIMENSION(:) :: mass
842 : LOGICAL :: converged
843 : REAL(KIND=dp) :: dx
844 : LOGICAL :: calc_intens
845 : INTEGER :: output_unit_ms
846 : TYPE(cp_logger_type), POINTER :: logger
847 : TYPE(cell_type), POINTER :: cell
848 :
849 : INTEGER :: i, j, jj, k, natoms, ncoord
850 162 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
851 : LOGICAL :: dump_only_positive
852 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
853 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H, H_save, residuum, tmp_b, tmp_s
854 324 : REAL(KIND=dp), DIMENSION(2, nrep) :: criteria
855 162 : REAL(Kind=dp), DIMENSION(:), POINTER :: intensities
856 :
857 162 : natoms = SIZE(particles)
858 162 : ncoord = 3*natoms
859 162 : nrep = SIZE(rep_env%f, 2)
860 :
861 : !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
862 162 : IF (ms_vib%mat_size /= 0) THEN
863 :
864 690 : ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
865 414 : ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
866 :
867 153792 : tmp_b(:, :) = ms_vib%b_mat
868 153792 : tmp_s(:, :) = ms_vib%s_mat
869 :
870 138 : DEALLOCATE (ms_vib%b_mat)
871 138 : DEALLOCATE (ms_vib%s_mat)
872 : END IF
873 :
874 810 : ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
875 486 : ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
876 :
877 176832 : ms_vib%s_mat = 0.0_dp
878 :
879 22896 : DO i = 1, 3*natoms
880 22734 : IF (ms_vib%mat_size /= 0) THEN
881 172878 : DO j = 1, ms_vib%mat_size
882 152730 : ms_vib%b_mat(i, j) = tmp_b(i, j)
883 172878 : ms_vib%s_mat(i, j) = tmp_s(i, j)
884 : END DO
885 : END IF
886 45738 : DO j = 1, nrep
887 45576 : ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
888 : END DO
889 : END DO
890 :
891 162 : IF (ms_vib%mat_size /= 0) THEN
892 138 : DEALLOCATE (tmp_s)
893 138 : DEALLOCATE (tmp_b)
894 : END IF
895 :
896 162 : ms_vib%mat_size = ms_vib%mat_size + nrep
897 :
898 648 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
899 486 : ALLOCATE (H_save(ms_vib%mat_size, ms_vib%mat_size))
900 486 : ALLOCATE (eigenval(ms_vib%mat_size))
901 :
902 : !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
903 :
904 336 : DO i = 1, nrep
905 23178 : DO j = 1, 3*natoms
906 23016 : ms_vib%s_mat(j, ms_vib%mat_size - nrep + i) = -(ms_vib%ms_force(j, i) - rep_env%f(j, i))/(2*ms_vib%step_b(i)*mass(j))
907 : END DO
908 : END DO
909 :
910 : CALL dgemm('T', 'N', ms_vib%mat_size, ms_vib%mat_size, SIZE(ms_vib%s_mat, 1), 1._dp, ms_vib%b_mat, SIZE(ms_vib%b_mat, 1), &
911 162 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
912 14006 : H_save(:, :) = approx_H
913 :
914 162 : CALL diamat_all(approx_H, eigenval)
915 :
916 : !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
917 486 : ALLOCATE (ind(ms_vib%mat_size))
918 648 : ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
919 :
920 162 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
921 :
922 336 : DO i = 1, nrep
923 7950 : DO j = 1, natoms
924 30630 : DO k = 1, 3
925 22842 : jj = (j - 1)*3 + k
926 30456 : ms_vib%delta_vec(jj, i) = ms_vib%b_vec(jj, i)/mass(jj)
927 : END DO
928 : END DO
929 : END DO
930 :
931 336 : DO i = 1, nrep
932 23016 : ms_vib%step_r(i) = dx/NORM2(ms_vib%delta_vec(:, i))
933 23178 : ms_vib%step_b(i) = NORM2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
934 : END DO
935 162 : converged = .FALSE.
936 : IF (MAXVAL(criteria(1, :)) <= ms_vib%eps(1) .AND. MAXVAL(criteria(2, :)) &
937 510 : <= ms_vib%eps(2) .OR. ms_vib%mat_size >= ncoord) converged = .TRUE.
938 486 : ALLOCATE (freq(nrep))
939 336 : DO i = 1, nrep
940 336 : freq(i) = SQRT(ABS(eigenval(ind(i)))*massunit)*vibfac
941 : END DO
942 :
943 : !!! write information and output !!!
944 162 : IF (converged) THEN
945 198 : eigenval(:) = SIGN(1._dp, eigenval(:))*SQRT(ABS(eigenval(:))*massunit)*vibfac
946 96 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
947 24 : tmp_b = 0._dp
948 72 : ALLOCATE (tmp_s(3, ms_vib%mat_size))
949 24 : tmp_s = 0._dp
950 24 : IF (calc_intens) THEN
951 60 : ALLOCATE (intensities(ms_vib%mat_size))
952 170 : intensities = 0._dp
953 : END IF
954 198 : DO i = 1, ms_vib%mat_size
955 2268 : DO j = 1, ms_vib%mat_size
956 331218 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
957 : END DO
958 45882 : tmp_b(:, i) = tmp_b(:, i)/NORM2(tmp_b(:, i))
959 : END DO
960 24 : IF (calc_intens) THEN
961 170 : DO i = 1, ms_vib%mat_size
962 2100 : DO j = 1, ms_vib%mat_size
963 7950 : tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_H(j, i)
964 : END DO
965 620 : IF (calc_intens) intensities(i) = NORM2(tmp_s(:, i))
966 : END DO
967 : END IF
968 24 : IF (calc_intens) THEN
969 : CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
970 : input, nrep, approx_H, eigenval, calc_intens, &
971 20 : intensities=intensities, logger=logger)
972 : ELSE
973 : CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
974 4 : input, nrep, approx_H, eigenval, calc_intens, logger=logger)
975 : END IF
976 24 : dump_only_positive = ms_vib%low_freq > 0.0_dp
977 : CALL write_vibrations_molden(input, particles, eigenval, tmp_b, intensities, calc_intens, &
978 24 : dump_only_positive=dump_only_positive, logger=logger, cell=cell)
979 24 : IF (calc_intens) THEN
980 20 : DEALLOCATE (intensities)
981 : END IF
982 24 : DEALLOCATE (tmp_b)
983 24 : DEALLOCATE (tmp_s)
984 : END IF
985 :
986 162 : IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
987 138 : ms_vib, input, nrep, approx_H, eigenval, calc_intens, logger=logger)
988 :
989 162 : DEALLOCATE (freq)
990 162 : DEALLOCATE (approx_H)
991 162 : DEALLOCATE (eigenval)
992 162 : DEALLOCATE (residuum)
993 162 : DEALLOCATE (ind)
994 :
995 324 : END SUBROUTINE evaluate_H_update_b
996 :
997 : ! **************************************************************************************************
998 : !> \brief writes the output for a mode tracking calculation
999 : !> \param ms_vib ...
1000 : !> \param nrep ...
1001 : !> \param mass ...
1002 : !> \param ncoord ...
1003 : !> \param approx_H ...
1004 : !> \param eigenval ...
1005 : !> \param ind ...
1006 : !> \param residuum ...
1007 : !> \param criteria ...
1008 : !> \author Florian Schiffmann 11.2007
1009 : ! **************************************************************************************************
1010 166 : SUBROUTINE select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1011 :
1012 : TYPE(ms_vib_type) :: ms_vib
1013 : INTEGER :: nrep
1014 : REAL(Kind=dp), DIMENSION(:) :: mass
1015 : INTEGER :: ncoord
1016 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1017 : REAL(Kind=dp), DIMENSION(:) :: eigenval
1018 : INTEGER, DIMENSION(:) :: ind
1019 : REAL(KIND=dp), DIMENSION(:, :) :: residuum
1020 : REAL(KIND=dp), DIMENSION(2, nrep), OPTIONAL :: criteria
1021 :
1022 : INTEGER :: i, j, jj, k
1023 : REAL(KIND=dp) :: my_val, norm
1024 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1025 166 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1026 :
1027 498 : ALLOCATE (tmp(ms_vib%mat_size))
1028 :
1029 282 : SELECT CASE (ms_vib%select_id)
1030 : CASE (1)
1031 116 : my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1032 1006 : DO i = 1, ms_vib%mat_size
1033 1006 : tmp(i) = ABS(my_val - eigenval(i))
1034 : END DO
1035 116 : CALL sort(tmp, (ms_vib%mat_size), ind)
1036 15892 : residuum = 0._dp
1037 238 : DO j = 1, nrep
1038 1152 : DO i = 1, ms_vib%mat_size
1039 144040 : residuum(:, j) = residuum(:, j) + approx_H(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1040 : END DO
1041 : END DO
1042 : CASE (2)
1043 6 : CALL get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1044 : CASE (3)
1045 :
1046 176 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1047 44 : tmp_b = 0._dp
1048 :
1049 266 : DO i = 1, ms_vib%mat_size
1050 1728 : DO j = 1, ms_vib%mat_size
1051 268290 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
1052 : END DO
1053 78854 : tmp_b(:, i) = tmp_b(:, i)/NORM2(tmp_b(:, i))
1054 : END DO
1055 44 : tmp = 0._dp
1056 266 : DO i = 1, ms_vib%mat_size
1057 666 : DO j = 1, SIZE(ms_vib%inv_atoms)
1058 1998 : DO k = 1, 3
1059 1332 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1060 1776 : tmp(i) = tmp(i) + SQRT(tmp_b(jj, i)**2)
1061 : END DO
1062 : END DO
1063 266 : IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1064 222 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
1065 : ELSE
1066 0 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) <= ms_vib%inv_range(1)) tmp(i) = 0._dp
1067 0 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) >= ms_vib%inv_range(2)) tmp(i) = 0._dp
1068 : END IF
1069 : END DO
1070 266 : tmp(:) = -tmp(:)
1071 44 : CALL sort(tmp, (ms_vib%mat_size), ind)
1072 7876 : residuum(:, :) = 0._dp
1073 :
1074 88 : DO j = 1, nrep
1075 310 : DO i = 1, ms_vib%mat_size
1076 39560 : residuum(:, j) = residuum(:, j) + approx_H(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1077 : END DO
1078 : END DO
1079 210 : DEALLOCATE (tmp_b)
1080 : END SELECT
1081 :
1082 344 : DO j = 1, nrep
1083 1528 : DO i = 1, ms_vib%mat_size
1084 366822 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1085 : END DO
1086 : END DO
1087 166 : IF (PRESENT(criteria)) THEN
1088 336 : DO i = 1, nrep
1089 23016 : criteria(1, i) = MAXVAL((residuum(:, i)))
1090 23178 : criteria(2, i) = NORM2(residuum(:, i))
1091 : END DO
1092 : END IF
1093 :
1094 344 : DO i = 1, nrep
1095 23728 : norm = NORM2(residuum(:, i))
1096 23894 : residuum(:, i) = residuum(:, i)/norm
1097 : END DO
1098 :
1099 1826 : DO k = 1, 10
1100 3606 : DO j = 1, nrep
1101 13620 : DO i = 1, ms_vib%mat_size
1102 3666440 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1103 3668220 : residuum(:, j) = residuum(:, j)/NORM2(residuum(:, j))
1104 : END DO
1105 3440 : IF (nrep > 1) THEN
1106 720 : DO i = 1, nrep
1107 720 : IF (i /= j) THEN
1108 4560 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), residuum(:, i))*residuum(:, i)
1109 4560 : residuum(:, j) = residuum(:, j)/NORM2(residuum(:, j))
1110 : END IF
1111 : END DO
1112 : END IF
1113 : END DO
1114 : END DO
1115 23894 : ms_vib%b_vec = residuum
1116 166 : DEALLOCATE (tmp)
1117 166 : END SUBROUTINE select_vector
1118 :
1119 : ! **************************************************************************************************
1120 : !> \brief writes the output for a mode tracking calculation
1121 : !> \param iw ...
1122 : !> \param converged ...
1123 : !> \param freq ...
1124 : !> \param criter ...
1125 : !> \param ms_vib ...
1126 : !> \param input ...
1127 : !> \param nrep ...
1128 : !> \param approx_H ...
1129 : !> \param eigenval ...
1130 : !> \param calc_intens ...
1131 : !> \param intensities ...
1132 : !> \param logger ...
1133 : !> \author Florian Schiffmann 11.2007
1134 : ! **************************************************************************************************
1135 162 : SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1136 162 : approx_H, eigenval, calc_intens, intensities, logger)
1137 :
1138 : INTEGER :: iw
1139 : LOGICAL :: converged
1140 : REAL(KIND=dp), DIMENSION(:) :: freq
1141 : REAL(KIND=dp), DIMENSION(:, :) :: criter
1142 : TYPE(ms_vib_type) :: ms_vib
1143 : TYPE(section_vals_type), POINTER :: input
1144 : INTEGER :: nrep
1145 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1146 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1147 : LOGICAL :: calc_intens
1148 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: intensities
1149 : TYPE(cp_logger_type), POINTER :: logger
1150 :
1151 : INTEGER :: i, j, msunit
1152 : REAL(KIND=dp) :: crit_a, crit_b, fint, gintval
1153 162 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: residuum
1154 : TYPE(section_vals_type), POINTER :: ms_vib_section
1155 :
1156 : ms_vib_section => section_vals_get_subs_vals(input, &
1157 162 : "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1158 :
1159 162 : fint = 42.255_dp*massunit*debye**2*bohr**2
1160 :
1161 162 : IF (converged) THEN
1162 24 : IF (iw > 0) THEN
1163 12 : WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1164 26 : DO i = 1, nrep
1165 26 : WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1166 : END DO
1167 36 : ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1168 12 : WRITE (iw, '( /, 1X, 79("-") )')
1169 12 : WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1170 12 : IF (PRESENT(intensities)) THEN
1171 10 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1172 : ELSE
1173 2 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1174 : END IF
1175 99 : DO i = 1, SIZE(ms_vib%b_mat, 2)
1176 87 : residuum = 0._dp
1177 1134 : DO j = 1, SIZE(ms_vib%b_mat, 2)
1178 165609 : residuum(:) = residuum(:) + approx_H(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1179 : END DO
1180 1134 : DO j = 1, ms_vib%mat_size
1181 330084 : residuum(:) = residuum(:) - DOT_PRODUCT(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1182 : END DO
1183 11508 : crit_a = MAXVAL(residuum(:))
1184 11508 : crit_b = NORM2(residuum)
1185 99 : IF (PRESENT(intensities)) THEN
1186 75 : gintval = fint*intensities(i)**2
1187 75 : IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1188 28 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1189 26 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1190 : ELSE
1191 47 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,1X,F12.6,3X,E12.3,7X,E12.3,11X,A)') &
1192 47 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1193 : END IF
1194 : ELSE
1195 12 : IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1196 12 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1197 6 : 'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1198 : ELSE
1199 0 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1200 0 : 'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1201 : END IF
1202 : END IF
1203 : END DO
1204 12 : DEALLOCATE (residuum)
1205 :
1206 : msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1207 : "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1208 : file_status="REPLACE", file_form="UNFORMATTED", &
1209 12 : file_action="WRITE")
1210 :
1211 12 : IF (msunit > 0) THEN
1212 12 : WRITE (UNIT=msunit) ms_vib%mat_size
1213 11520 : WRITE (UNIT=msunit) ms_vib%b_mat
1214 11520 : WRITE (UNIT=msunit) ms_vib%s_mat
1215 312 : IF (calc_intens) WRITE (UNIT=msunit) ms_vib%dip_deriv
1216 : END IF
1217 :
1218 : CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1219 12 : "PRINT%MS_RESTART")
1220 : END IF
1221 : ELSE
1222 138 : IF (iw > 0) THEN
1223 : msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1224 : "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1225 : file_status="REPLACE", file_form="UNFORMATTED", &
1226 69 : file_action="WRITE")
1227 :
1228 69 : IF (msunit > 0) THEN
1229 69 : WRITE (UNIT=msunit) ms_vib%mat_size
1230 76896 : WRITE (UNIT=msunit) ms_vib%b_mat
1231 76896 : WRITE (UNIT=msunit) ms_vib%s_mat
1232 1869 : IF (calc_intens) WRITE (UNIT=msunit) ms_vib%dip_deriv
1233 : END IF
1234 :
1235 : CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1236 69 : "PRINT%MS_RESTART")
1237 :
1238 69 : WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1239 142 : DO i = 1, nrep
1240 142 : IF (criter(1, i) <= 1E-7 .AND. (criter(2, i)) <= 1E-6) THEN
1241 1 : WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1242 : ELSE
1243 72 : WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 NOT CONVERGED"
1244 : END IF
1245 : END DO
1246 : END IF
1247 : END IF
1248 :
1249 162 : END SUBROUTINE ms_out
1250 :
1251 : ! **************************************************************************************************
1252 : !> \brief ...
1253 : !> \param ms_vib ...
1254 : !> \param approx_H ...
1255 : !> \param eigenval ...
1256 : !> \param residuum ...
1257 : !> \param nrep ...
1258 : !> \param ind ...
1259 : !> \author Florian Schiffmann 11.2007
1260 : ! **************************************************************************************************
1261 6 : SUBROUTINE get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1262 :
1263 : TYPE(ms_vib_type) :: ms_vib
1264 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1265 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1266 : REAL(KIND=dp), DIMENSION(:, :) :: residuum
1267 : INTEGER :: nrep
1268 : INTEGER, DIMENSION(:) :: ind
1269 :
1270 : INTEGER :: count1, count2, i, j
1271 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1272 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1273 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1274 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1275 : REAL(KIND=dp), DIMENSION(2) :: myrange
1276 :
1277 18 : myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1278 6 : count1 = 0
1279 6 : count2 = 0
1280 126 : residuum = 0.0_dp
1281 6 : ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1282 18 : ALLOCATE (map1(SIZE(eigenval), 2))
1283 18 : ALLOCATE (tmp(SIZE(eigenval)))
1284 30 : DO i = 1, SIZE(eigenval)
1285 24 : IF (ABS(eigenval(i) - myrange(1)) + ABS(eigenval(i) - myrange(2)) <= &
1286 6 : ABS(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1287 0 : count1 = count1 + 1
1288 0 : map1(count1, 1) = i
1289 : ELSE
1290 24 : count2 = count2 + 1
1291 24 : map1(count2, 2) = i
1292 24 : tmp(count2) = MIN(ABS(eigenval(i) - myrange(1)), ABS(eigenval(i) - myrange(2)))
1293 : END IF
1294 : END DO
1295 :
1296 6 : IF (count1 == nrep) THEN
1297 0 : DO j = 1, count1
1298 0 : DO i = 1, ms_vib%mat_size
1299 0 : residuum(:, j) = residuum(:, j) + approx_H(i, map1(j, 1))*(ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1300 0 : ind(j) = map1(j, 1)
1301 : END DO
1302 : END DO
1303 6 : ELSE IF (count1 > nrep) THEN
1304 0 : ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1305 0 : ALLOCATE (tmp1(count1))
1306 0 : ALLOCATE (map2(count1))
1307 0 : tmp_resid = 0._dp
1308 0 : DO j = 1, count1
1309 0 : DO i = 1, ms_vib%mat_size
1310 : tmp_resid(:, j) = tmp_resid(:, j) + approx_H(i, map1(j, 1))* &
1311 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1312 : END DO
1313 : END DO
1314 :
1315 0 : DO j = 1, count1
1316 0 : DO i = 1, ms_vib%mat_size
1317 0 : tmp_resid(:, j) = tmp_resid(:, j) - DOT_PRODUCT(tmp_resid(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1318 : END DO
1319 0 : tmp(j) = MAXVAL(tmp_resid(:, j))
1320 : END DO
1321 0 : CALL sort(tmp, count1, map2)
1322 0 : DO j = 1, nrep
1323 0 : residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1324 0 : ind(j) = map1(map2(count1 + 1 - j), 1)
1325 : END DO
1326 0 : DEALLOCATE (tmp_resid)
1327 0 : DEALLOCATE (tmp1)
1328 0 : DEALLOCATE (map2)
1329 6 : ELSE IF (count1 < nrep) THEN
1330 :
1331 18 : ALLOCATE (map2(count2))
1332 6 : IF (count1 /= 0) THEN
1333 0 : DO j = 1, count1
1334 0 : DO i = 1, ms_vib%mat_size
1335 : residuum(:, j) = residuum(:, j) + approx_H(i, map1(j, 1))* &
1336 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1337 : END DO
1338 0 : ind(j) = map1(j, 1)
1339 : END DO
1340 : END IF
1341 6 : CALL sort(tmp, count2, map2)
1342 18 : DO j = 1, nrep - count1
1343 60 : DO i = 1, ms_vib%mat_size
1344 : residuum(:, count1 + j) = residuum(:, count1 + j) + approx_H(i, map1(map2(j), 2)) &
1345 492 : *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1346 : END DO
1347 18 : ind(count1 + j) = map1(map2(j), 2)
1348 : END DO
1349 :
1350 6 : DEALLOCATE (map2)
1351 : END IF
1352 :
1353 6 : DEALLOCATE (map1)
1354 6 : DEALLOCATE (tmp)
1355 :
1356 6 : END SUBROUTINE get_vibs_in_range
1357 0 : END MODULE mode_selective
|
