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 bfgs_optimizer, ONLY: bfgs_read_restart_header
22 : USE cell_types, ONLY: cell_type
23 : USE cp_files, ONLY: close_file,&
24 : open_file
25 : USE cp_log_handling, ONLY: cp_get_default_logger,&
26 : cp_logger_get_default_io_unit,&
27 : cp_logger_type
28 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
29 : cp_print_key_unit_nr
30 : USE cp_result_methods, ONLY: get_results
31 : USE global_types, ONLY: global_environment_type
32 : USE ieee_arithmetic, ONLY: ieee_is_finite
33 : USE input_constants, ONLY: bfgs_direct,&
34 : bfgs_inverse,&
35 : ms_guess_atomic,&
36 : ms_guess_bfgs,&
37 : ms_guess_molden,&
38 : ms_guess_restart,&
39 : ms_guess_restart_vec
40 : USE input_section_types, ONLY: section_vals_get,&
41 : section_vals_get_subs_vals,&
42 : section_vals_type,&
43 : section_vals_val_get
44 : USE iso_fortran_env, ONLY: iostat_end
45 : USE kinds, ONLY: default_path_length,&
46 : default_string_length,&
47 : dp,&
48 : max_line_length
49 : USE mathlib, ONLY: diamat_all
50 : USE message_passing, ONLY: mp_para_env_type
51 : USE molden_utils, ONLY: write_vibrations_molden
52 : USE particle_types, ONLY: particle_type
53 : USE physcon, ONLY: bohr,&
54 : debye,&
55 : massunit,&
56 : vibfac
57 : USE replica_methods, ONLY: rep_env_calc_e_f
58 : USE replica_types, ONLY: replica_env_type
59 : USE util, ONLY: sort
60 : #include "./base/base_uses.f90"
61 :
62 : IMPLICIT NONE
63 :
64 : PRIVATE
65 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'mode_selective'
66 : LOGICAL, PARAMETER :: debug_this_module = .FALSE.
67 : REAL(KIND=dp), PARAMETER :: bfgs_guess_degeneracy_tol = 1.0E-8_dp
68 :
69 : TYPE ms_vib_type
70 : INTEGER :: mat_size = -1
71 : INTEGER :: select_id = -1
72 : INTEGER, DIMENSION(:), POINTER :: inv_atoms => NULL()
73 : REAL(KIND=dp) :: eps(2) = 0.0_dp
74 : REAL(KIND=dp) :: sel_freq = 0.0_dp
75 : REAL(KIND=dp) :: low_freq = 0.0_dp
76 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: b_vec => NULL()
77 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: delta_vec => NULL()
78 : REAL(KIND=dp), POINTER, DIMENSION(:, :) :: ms_force => NULL()
79 : REAL(KIND=dp), DIMENSION(:), POINTER :: eig_bfgs => NULL()
80 : REAL(KIND=dp), DIMENSION(:), POINTER :: f_range => NULL()
81 : REAL(KIND=dp), DIMENSION(:), POINTER :: inv_range => NULL()
82 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_b => NULL()
83 : REAL(KIND=dp), POINTER, DIMENSION(:) :: step_r => NULL()
84 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: b_mat => NULL()
85 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dip_deriv => NULL()
86 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: hes_bfgs => NULL()
87 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: s_mat => NULL()
88 : INTEGER :: initial_guess = -1
89 : END TYPE ms_vib_type
90 :
91 : PUBLIC :: ms_vb_anal
92 :
93 : CONTAINS
94 : ! **************************************************************************************************
95 : !> \brief Module performing a vibrational analysis
96 : !> \param input ...
97 : !> \param rep_env ...
98 : !> \param para_env ...
99 : !> \param globenv ...
100 : !> \param particles ...
101 : !> \param nrep ...
102 : !> \param calc_intens ...
103 : !> \param dx ...
104 : !> \param output_unit ...
105 : !> \param logger ...
106 : !> \param cell simulation cell
107 : !> \author Teodoro Laino 08.2006
108 : ! **************************************************************************************************
109 28 : SUBROUTINE ms_vb_anal(input, rep_env, para_env, globenv, particles, &
110 : nrep, calc_intens, dx, output_unit, logger, cell)
111 : TYPE(section_vals_type), POINTER :: input
112 : TYPE(replica_env_type), POINTER :: rep_env
113 : TYPE(mp_para_env_type), POINTER :: para_env
114 : TYPE(global_environment_type), POINTER :: globenv
115 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
116 : INTEGER :: nrep
117 : LOGICAL :: calc_intens
118 : REAL(KIND=dp) :: dx
119 : INTEGER :: output_unit
120 : TYPE(cp_logger_type), POINTER :: logger
121 : TYPE(cell_type), POINTER :: cell
122 :
123 : CHARACTER(len=*), PARAMETER :: routineN = 'ms_vb_anal'
124 :
125 : CHARACTER(LEN=default_string_length) :: description
126 : INTEGER :: handle, i, ip1, j, natoms, ncoord
127 : LOGICAL :: converged
128 28 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: mass, pos0
129 28 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: tmp_deriv
130 28 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: tmp_dip
131 : TYPE(ms_vib_type) :: ms_vib
132 :
133 28 : CALL timeset(routineN, handle)
134 28 : converged = .FALSE.
135 28 : natoms = SIZE(particles)
136 28 : ncoord = 3*natoms
137 112 : ALLOCATE (mass(3*natoms))
138 1068 : DO i = 1, natoms
139 4188 : DO j = 1, 3
140 3120 : mass((i - 1)*3 + j) = particles(i)%atomic_kind%mass
141 4160 : mass((i - 1)*3 + j) = SQRT(mass((i - 1)*3 + j))
142 : END DO
143 : END DO
144 : ! Allocate working arrays
145 112 : ALLOCATE (ms_vib%delta_vec(ncoord, nrep))
146 84 : ALLOCATE (ms_vib%b_vec(ncoord, nrep))
147 84 : ALLOCATE (ms_vib%step_r(nrep))
148 56 : ALLOCATE (ms_vib%step_b(nrep))
149 28 : IF (calc_intens) THEN
150 22 : description = '[DIPOLE]'
151 88 : ALLOCATE (tmp_dip(nrep, 3, 2))
152 66 : ALLOCATE (ms_vib%dip_deriv(3, nrep))
153 : END IF
154 : CALL MS_initial_moves(para_env, nrep, input, globenv, ms_vib, &
155 : particles, &
156 : mass, &
157 : dx, &
158 28 : calc_intens, logger)
159 28 : ncoord = 3*natoms
160 84 : ALLOCATE (pos0(ncoord))
161 112 : ALLOCATE (ms_vib%ms_force(ncoord, nrep))
162 1068 : DO i = 1, natoms
163 4188 : DO j = 1, 3
164 4160 : pos0((i - 1)*3 + j) = particles((i))%r(j)
165 : END DO
166 : END DO
167 186 : ncoord = 3*natoms
168 : DO
169 26952 : ms_vib%ms_force = HUGE(0.0_dp)
170 384 : DO i = 1, nrep
171 26952 : DO j = 1, ncoord
172 26766 : rep_env%r(j, i) = pos0(j) + ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
173 : END DO
174 : END DO
175 186 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
176 :
177 384 : DO i = 1, nrep
178 198 : IF (calc_intens) THEN
179 : CALL get_results(results=rep_env%results(i)%results, &
180 : description=description, &
181 168 : n_rep=ip1)
182 : CALL get_results(results=rep_env%results(i)%results, &
183 : description=description, &
184 : values=tmp_dip(i, :, 1), &
185 168 : nval=ip1)
186 : END IF
187 26952 : DO j = 1, ncoord
188 26766 : ms_vib%ms_force(j, i) = rep_env%f(j, i)
189 : END DO
190 : END DO
191 384 : DO i = 1, nrep
192 26952 : DO j = 1, ncoord
193 26766 : rep_env%r(j, i) = pos0(j) - ms_vib%step_r(i)*ms_vib%delta_vec(j, i)
194 : END DO
195 : END DO
196 186 : CALL rep_env_calc_e_f(rep_env, calc_f=.TRUE.)
197 186 : IF (calc_intens) THEN
198 336 : DO i = 1, nrep
199 : CALL get_results(results=rep_env%results(i)%results, &
200 : description=description, &
201 168 : n_rep=ip1)
202 : CALL get_results(results=rep_env%results(i)%results, &
203 : description=description, &
204 : values=tmp_dip(i, :, 2), &
205 168 : nval=ip1)
206 1008 : ms_vib%dip_deriv(:, ms_vib%mat_size + i) = (tmp_dip(i, :, 1) - tmp_dip(i, :, 2))/(2*ms_vib%step_b(i))
207 : END DO
208 : END IF
209 :
210 : CALL evaluate_H_update_b(rep_env, ms_vib, input, nrep, &
211 : particles, &
212 : mass, &
213 : converged, &
214 : dx, calc_intens, &
215 186 : output_unit, logger, cell)
216 186 : IF (converged) EXIT
217 186 : IF (calc_intens) THEN
218 438 : ALLOCATE (tmp_deriv(3, ms_vib%mat_size))
219 4034 : tmp_deriv = ms_vib%dip_deriv
220 146 : DEALLOCATE (ms_vib%dip_deriv)
221 438 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
222 4034 : ms_vib%dip_deriv(:, 1:ms_vib%mat_size) = tmp_deriv(:, 1:ms_vib%mat_size)
223 146 : DEALLOCATE (tmp_deriv)
224 : END IF
225 : END DO
226 28 : DEALLOCATE (ms_vib%ms_force)
227 28 : DEALLOCATE (pos0)
228 28 : DEALLOCATE (ms_vib%step_r)
229 28 : DEALLOCATE (ms_vib%step_b)
230 28 : DEALLOCATE (ms_vib%b_vec)
231 28 : DEALLOCATE (ms_vib%delta_vec)
232 28 : DEALLOCATE (mass)
233 28 : DEALLOCATE (ms_vib%b_mat)
234 28 : DEALLOCATE (ms_vib%s_mat)
235 28 : IF (ms_vib%select_id == 3) THEN
236 12 : DEALLOCATE (ms_vib%inv_atoms)
237 : END IF
238 28 : IF (ASSOCIATED(ms_vib%eig_bfgs)) THEN
239 0 : DEALLOCATE (ms_vib%eig_bfgs)
240 : END IF
241 28 : IF (ASSOCIATED(ms_vib%hes_bfgs)) THEN
242 0 : DEALLOCATE (ms_vib%hes_bfgs)
243 : END IF
244 28 : IF (calc_intens) THEN
245 22 : DEALLOCATE (ms_vib%dip_deriv)
246 22 : DEALLOCATE (tmp_dip)
247 : END IF
248 28 : CALL timestop(handle)
249 84 : END SUBROUTINE ms_vb_anal
250 : ! **************************************************************************************************
251 : !> \brief Generates the first displacement vector for a mode selctive vibrational
252 : !> analysis. At the moment this is a random number for selected atoms
253 : !> \param para_env ...
254 : !> \param nrep ...
255 : !> \param input ...
256 : !> \param globenv ...
257 : !> \param ms_vib ...
258 : !> \param particles ...
259 : !> \param mass ...
260 : !> \param dx ...
261 : !> \param calc_intens ...
262 : !> \param logger ...
263 : !> \author Florian Schiffmann 11.2007
264 : ! **************************************************************************************************
265 28 : SUBROUTINE MS_initial_moves(para_env, nrep, input, globenv, ms_vib, particles, &
266 28 : mass, dx, &
267 : calc_intens, logger)
268 : TYPE(mp_para_env_type), POINTER :: para_env
269 : INTEGER :: nrep
270 : TYPE(section_vals_type), POINTER :: input
271 : TYPE(global_environment_type), POINTER :: globenv
272 : TYPE(ms_vib_type) :: ms_vib
273 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
274 : REAL(Kind=dp), DIMENSION(:) :: mass
275 : REAL(KIND=dp) :: dx
276 : LOGICAL :: calc_intens
277 : TYPE(cp_logger_type), POINTER :: logger
278 :
279 : CHARACTER(len=*), PARAMETER :: routineN = 'MS_initial_moves'
280 :
281 : INTEGER :: guess, handle, i, j, jj, k, m, &
282 : n_rep_val, natoms, ncoord
283 28 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map_atoms
284 28 : INTEGER, DIMENSION(:), POINTER :: tmplist
285 : LOGICAL :: do_involved_atoms, ionode
286 : REAL(KIND=dp) :: my_val, norm
287 : TYPE(section_vals_type), POINTER :: involved_at_section, ms_vib_section
288 :
289 28 : CALL timeset(routineN, handle)
290 28 : NULLIFY (ms_vib%eig_bfgs, ms_vib%f_range, ms_vib%hes_bfgs, ms_vib%inv_range)
291 28 : ms_vib_section => section_vals_get_subs_vals(input, "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
292 28 : CALL section_vals_val_get(ms_vib_section, "INITIAL_GUESS", i_val=guess)
293 28 : CALL section_vals_val_get(ms_vib_section, "EPS_MAX_VAL", r_val=ms_vib%eps(1))
294 28 : CALL section_vals_val_get(ms_vib_section, "EPS_NORM", r_val=ms_vib%eps(2))
295 28 : CALL section_vals_val_get(ms_vib_section, "RANGE", n_rep_val=n_rep_val)
296 28 : ms_vib%select_id = 0
297 28 : IF (n_rep_val /= 0) THEN
298 2 : CALL section_vals_val_get(ms_vib_section, "RANGE", r_vals=ms_vib%f_range)
299 2 : IF (ms_vib%f_range(1) > ms_vib%f_range(2)) THEN
300 0 : my_val = ms_vib%f_range(2)
301 0 : ms_vib%f_range(2) = ms_vib%f_range(1)
302 0 : ms_vib%f_range(1) = my_val
303 : END IF
304 2 : ms_vib%select_id = 2
305 : END IF
306 28 : CALL section_vals_val_get(ms_vib_section, "FREQUENCY", r_val=ms_vib%sel_freq)
307 28 : CALL section_vals_val_get(ms_vib_section, "LOWEST_FREQUENCY", r_val=ms_vib%low_freq)
308 28 : IF (ms_vib%sel_freq > 0._dp) ms_vib%select_id = 1
309 28 : involved_at_section => section_vals_get_subs_vals(ms_vib_section, "INVOLVED_ATOMS")
310 28 : CALL section_vals_get(involved_at_section, explicit=do_involved_atoms)
311 28 : IF (do_involved_atoms) THEN
312 12 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", n_rep_val=n_rep_val)
313 12 : jj = 0
314 24 : DO k = 1, n_rep_val
315 12 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=k, i_vals=tmplist)
316 48 : DO j = 1, SIZE(tmplist)
317 36 : jj = jj + 1
318 : END DO
319 : END DO
320 12 : IF (jj >= 1) THEN
321 12 : natoms = jj
322 36 : ALLOCATE (ms_vib%inv_atoms(natoms))
323 12 : jj = 0
324 24 : DO m = 1, n_rep_val
325 12 : CALL section_vals_val_get(involved_at_section, "INVOLVED_ATOMS", i_rep_val=m, i_vals=tmplist)
326 48 : DO j = 1, SIZE(tmplist)
327 36 : ms_vib%inv_atoms(j) = tmplist(j)
328 : END DO
329 : END DO
330 12 : ms_vib%select_id = 3
331 : END IF
332 12 : CALL section_vals_val_get(involved_at_section, "RANGE", n_rep_val=n_rep_val)
333 12 : IF (n_rep_val /= 0) THEN
334 0 : CALL section_vals_val_get(involved_at_section, "RANGE", r_vals=ms_vib%inv_range)
335 0 : IF (ms_vib%inv_range(1) > ms_vib%inv_range(2)) THEN
336 0 : ms_vib%inv_range(2) = my_val
337 0 : ms_vib%inv_range(2) = ms_vib%inv_range(1)
338 0 : ms_vib%inv_range(1) = my_val
339 : END IF
340 : END IF
341 : END IF
342 28 : IF (ms_vib%select_id == 0) THEN
343 0 : CPABORT("no frequency, range or involved atoms specified ")
344 : END IF
345 28 : ionode = para_env%is_source()
346 12 : SELECT CASE (guess)
347 : CASE (ms_guess_atomic)
348 12 : ms_vib%initial_guess = 1
349 12 : CALL section_vals_val_get(ms_vib_section, "ATOMS", n_rep_val=n_rep_val)
350 12 : jj = 0
351 22 : DO k = 1, n_rep_val
352 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=k, i_vals=tmplist)
353 42 : DO j = 1, SIZE(tmplist)
354 30 : jj = jj + 1
355 : END DO
356 : END DO
357 12 : IF (jj < 1) THEN
358 2 : natoms = SIZE(particles)
359 6 : ALLOCATE (map_atoms(natoms))
360 14 : DO j = 1, natoms
361 14 : map_atoms(j) = j
362 : END DO
363 : ELSE
364 10 : natoms = jj
365 30 : ALLOCATE (map_atoms(natoms))
366 10 : jj = 0
367 20 : DO m = 1, n_rep_val
368 10 : CALL section_vals_val_get(ms_vib_section, "ATOMS", i_rep_val=m, i_vals=tmplist)
369 40 : DO j = 1, SIZE(tmplist)
370 30 : map_atoms(j) = tmplist(j)
371 : END DO
372 : END DO
373 : END IF
374 :
375 : ! apply random displacement along the mass weighted nuclear cartesian coordinates
376 526 : ms_vib%b_vec = 0._dp
377 526 : ms_vib%delta_vec = 0._dp
378 12 : jj = 0
379 :
380 28 : DO i = 1, nrep
381 56 : DO j = 1, natoms
382 176 : DO k = 1, 3
383 120 : jj = (map_atoms(j) - 1)*3 + k
384 160 : ms_vib%b_vec(jj, i) = ABS(globenv%gaussian_rng_stream%next())
385 : END DO
386 : END DO
387 514 : norm = NORM2(ms_vib%b_vec(:, i))
388 526 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
389 : END DO
390 :
391 12 : IF (nrep > 1) THEN
392 44 : DO k = 1, 10
393 124 : DO j = 1, nrep
394 280 : DO i = 1, nrep
395 240 : IF (i /= j) THEN
396 : ms_vib%b_vec(:, j) = &
397 1520 : ms_vib%b_vec(:, j) - DOT_PRODUCT(ms_vib%b_vec(:, j), ms_vib%b_vec(:, i))*ms_vib%b_vec(:, i)
398 : ms_vib%b_vec(:, j) = &
399 1520 : ms_vib%b_vec(:, j)/NORM2(ms_vib%b_vec(:, j))
400 : END IF
401 : END DO
402 : END DO
403 : END DO
404 : END IF
405 :
406 12 : ms_vib%mat_size = 0
407 474 : DO i = 1, SIZE(ms_vib%b_vec, 1)
408 972 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
409 : END DO
410 : CASE (ms_guess_bfgs)
411 :
412 8 : ms_vib%initial_guess = 2
413 8 : CALL bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
414 8 : ms_vib%mat_size = 0
415 :
416 : CASE (ms_guess_restart_vec)
417 :
418 4 : ms_vib%initial_guess = 3
419 : ncoord = 3*SIZE(particles)
420 4 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
421 :
422 4 : ms_vib%mat_size = 0
423 : CASE (ms_guess_restart)
424 0 : ms_vib%initial_guess = 4
425 : ncoord = 3*SIZE(particles)
426 0 : CALL rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
427 :
428 : CASE (ms_guess_molden)
429 4 : ms_vib%initial_guess = 5
430 4 : ncoord = 3*SIZE(particles)
431 4 : CALL molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
432 32 : ms_vib%mat_size = 0
433 : END SELECT
434 6404 : CALL para_env%bcast(ms_vib%b_vec)
435 6404 : CALL para_env%bcast(ms_vib%delta_vec)
436 60 : DO i = 1, nrep
437 3188 : ms_vib%step_r(i) = dx/NORM2(ms_vib%delta_vec(:, i))
438 3216 : ms_vib%step_b(i) = NORM2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
439 : END DO
440 28 : CALL timestop(handle)
441 :
442 56 : END SUBROUTINE MS_initial_moves
443 :
444 : ! **************************************************************************************************
445 : !> \brief Read a legacy or tagged BFGS restart as a direct ionic Hessian.
446 : !> \param unit File unit
447 : !> \param ncoord Number of ionic coordinates needed by mode-selective analysis
448 : !> \param hessian Direct ionic Hessian
449 : ! **************************************************************************************************
450 4 : SUBROUTINE read_bfgs_guess_hessian(unit, ncoord, hessian)
451 : INTEGER, INTENT(IN) :: unit, ncoord
452 : REAL(KIND=dp), DIMENSION(ncoord, ncoord), &
453 : INTENT(OUT) :: hessian
454 :
455 : INTEGER :: i, ios, stored_method, stored_ndf, &
456 : stored_type
457 : LOGICAL :: tagged
458 : REAL(KIND=dp) :: extra
459 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: column, eigenvalues
460 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: stored_matrix
461 : REAL(KIND=dp), DIMENSION(3, 3) :: cell_frame, reference_cell
462 :
463 : CALL bfgs_read_restart_header(unit, stored_method, stored_ndf, stored_type, tagged, &
464 4 : reference_cell, cell_frame)
465 :
466 4 : IF (.NOT. tagged) THEN
467 : ! Legacy restart: preserve the historical behaviour. Reading ncoord values from
468 : ! the first ncoord records also extracts the ionic block from a CELL_OPT Hessian.
469 0 : DO i = 1, ncoord
470 0 : READ (unit, IOSTAT=ios) hessian(:, i)
471 0 : IF (ios /= 0) CPABORT("Error while reading legacy BFGS Hessian")
472 : END DO
473 : RETURN
474 : END IF
475 :
476 4 : IF (stored_ndf < ncoord) THEN
477 0 : CPABORT("BFGS restart matrix is smaller than the mode-selective ionic Hessian")
478 : END IF
479 :
480 2 : SELECT CASE (stored_method)
481 : CASE (bfgs_direct)
482 6 : ALLOCATE (column(stored_ndf))
483 356 : DO i = 1, stored_ndf
484 354 : READ (unit, IOSTAT=ios) column(:)
485 354 : IF (ios /= 0) CPABORT("Truncated BFGS restart matrix")
486 63012 : IF (.NOT. ALL(ieee_is_finite(column(:)))) CPABORT("Non-finite BFGS restart matrix")
487 63014 : IF (i <= ncoord) hessian(:, i) = column(1:ncoord)
488 : END DO
489 2 : DEALLOCATE (column)
490 :
491 : CASE (bfgs_inverse)
492 12 : ALLOCATE (stored_matrix(stored_ndf, stored_ndf), eigenvalues(stored_ndf))
493 275 : DO i = 1, stored_ndf
494 273 : READ (unit, IOSTAT=ios) stored_matrix(:, i)
495 275 : IF (ios /= 0) CPABORT("Truncated inverse BFGS restart matrix")
496 : END DO
497 40820 : IF (.NOT. ALL(ieee_is_finite(stored_matrix(:, :)))) THEN
498 0 : CPABORT("Non-finite inverse BFGS restart matrix")
499 : END IF
500 :
501 : ! Invert the complete matrix before taking the ionic block. For CELL_OPT,
502 : ! inverse(H_RR) is not in general equal to the RR block of inverse(H).
503 2 : CALL diamat_all(stored_matrix, eigenvalues)
504 275 : IF (.NOT. ALL(ieee_is_finite(eigenvalues(:)))) THEN
505 0 : CPABORT("Non-finite inverse BFGS restart spectrum")
506 : END IF
507 275 : IF (MINVAL(eigenvalues(:)) <= 0.0_dp) THEN
508 0 : CPABORT("Inverse BFGS restart matrix is not positive definite")
509 : END IF
510 275 : eigenvalues(:) = 1.0_dp/eigenvalues(:)
511 275 : IF (.NOT. ALL(ieee_is_finite(eigenvalues(:)))) THEN
512 0 : CPABORT("Unrepresentable inverse BFGS restart matrix")
513 : END IF
514 275 : DO i = 1, stored_ndf
515 : stored_matrix(1:ncoord, i) = &
516 40244 : SQRT(eigenvalues(i))*stored_matrix(1:ncoord, i)
517 : END DO
518 4 : hessian(:, :) = MATMUL(stored_matrix(1:ncoord, :), &
519 2 : TRANSPOSE(stored_matrix(1:ncoord, :)))
520 2 : DEALLOCATE (stored_matrix, eigenvalues)
521 :
522 : CASE DEFAULT
523 4 : CPABORT("Unknown matrix representation in BFGS restart")
524 : END SELECT
525 :
526 : ! Tagged files have a known matrix dimension, so reject trailing/corrupt data.
527 4 : READ (unit, IOSTAT=ios) extra
528 4 : IF (ios /= iostat_end) CPABORT("Extra data in BFGS restart matrix")
529 : MARK_USED(stored_type)
530 : MARK_USED(reference_cell)
531 : MARK_USED(cell_frame)
532 4 : END SUBROUTINE read_bfgs_guess_hessian
533 :
534 : ! **************************************************************************************************
535 : !> \brief ...
536 : !> \param ms_vib_section ...
537 : !> \param ms_vib ...
538 : !> \param particles ...
539 : !> \param mass ...
540 : !> \param para_env ...
541 : !> \param nrep ...
542 : !> \author Florian Schiffmann 11.2007
543 : ! **************************************************************************************************
544 8 : SUBROUTINE bfgs_guess(ms_vib_section, ms_vib, particles, mass, para_env, nrep)
545 :
546 : TYPE(section_vals_type), POINTER :: ms_vib_section
547 : TYPE(ms_vib_type) :: ms_vib
548 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
549 : REAL(Kind=dp), DIMENSION(:) :: mass
550 : TYPE(mp_para_env_type), POINTER :: para_env
551 : INTEGER :: nrep
552 :
553 : CHARACTER(LEN=default_path_length) :: hes_filename
554 : INTEGER :: hesunit, i, j, jj, k, natoms, ncoord, &
555 : output_unit
556 8 : INTEGER, DIMENSION(:), POINTER :: tmplist
557 : REAL(KIND=dp) :: my_val, norm
558 8 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
559 : TYPE(cp_logger_type), POINTER :: logger
560 :
561 16 : logger => cp_get_default_logger()
562 8 : output_unit = cp_logger_get_default_io_unit(logger)
563 :
564 8 : natoms = SIZE(particles)
565 8 : ncoord = 3*natoms
566 :
567 32 : ALLOCATE (ms_vib%hes_bfgs(ncoord, ncoord))
568 24 : ALLOCATE (ms_vib%eig_bfgs(ncoord))
569 :
570 8 : IF (para_env%is_source()) THEN
571 4 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=hes_filename)
572 4 : IF (hes_filename == "") hes_filename = "HESSIAN"
573 : CALL open_file(file_name=hes_filename, file_status="OLD", &
574 4 : file_form="UNFORMATTED", file_action="READ", unit_number=hesunit)
575 12 : ALLOCATE (tmp(ncoord))
576 12 : ALLOCATE (tmplist(ncoord))
577 :
578 4 : CALL read_bfgs_guess_hessian(hesunit, ncoord, ms_vib%hes_bfgs)
579 4 : CALL close_file(hesunit)
580 4 : IF (output_unit > 0) THEN
581 : WRITE (output_unit, FMT="(/,T2,A)") &
582 4 : "*** Initial Hessian has been read successfully ***"
583 : END IF
584 625 : DO i = 1, ncoord
585 102712 : DO j = 1, ncoord
586 102708 : ms_vib%hes_bfgs(i, j) = ms_vib%hes_bfgs(i, j)/(mass(i)*mass(j))
587 : END DO
588 : END DO
589 :
590 4 : CALL diamat_all(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
591 4 : tmp(:) = 0._dp
592 4 : IF (ms_vib%select_id == 1) my_val = (ms_vib%sel_freq/vibfac)**2/massunit
593 4 : IF (ms_vib%select_id == 2) my_val = (((ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp)/vibfac)**2/massunit
594 4 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
595 178 : DO i = 1, ncoord
596 178 : tmp(i) = ABS(my_val - ms_vib%eig_bfgs(i))
597 : END DO
598 3 : ELSE IF (ms_vib%select_id == 3) THEN
599 3 : CALL bfgs_canonicalize_degenerate_modes(ms_vib%hes_bfgs, ms_vib%eig_bfgs)
600 447 : DO i = 1, ncoord
601 1332 : DO j = 1, SIZE(ms_vib%inv_atoms)
602 3996 : DO k = 1, 3
603 2664 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
604 3552 : tmp(i) = tmp(i) + SQRT(ms_vib%hes_bfgs(jj, i)**2)
605 : END DO
606 : END DO
607 447 : IF ((SIGN(1._dp, ms_vib%eig_bfgs(i))*SQRT(ABS(ms_vib%eig_bfgs(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
608 : END DO
609 447 : tmp(:) = -tmp(:)
610 : END IF
611 4 : CALL sort(tmp, ncoord, tmplist)
612 8 : DO i = 1, nrep
613 625 : ms_vib%b_vec(:, i) = ms_vib%hes_bfgs(:, tmplist(i))
614 625 : norm = NORM2(ms_vib%b_vec(:, i))
615 629 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
616 : END DO
617 625 : DO i = 1, SIZE(ms_vib%b_vec, 1)
618 1246 : ms_vib%delta_vec(i, :) = ms_vib%b_vec(i, :)/mass(i)
619 : END DO
620 4 : DEALLOCATE (tmp)
621 4 : DEALLOCATE (tmplist)
622 : END IF
623 :
624 2508 : CALL para_env%bcast(ms_vib%b_vec)
625 2508 : CALL para_env%bcast(ms_vib%delta_vec)
626 :
627 8 : DEALLOCATE (ms_vib%hes_bfgs)
628 8 : DEALLOCATE (ms_vib%eig_bfgs)
629 8 : ms_vib%mat_size = 0
630 :
631 8 : END SUBROUTINE bfgs_guess
632 :
633 : ! **************************************************************************************************
634 : !> \brief Replace arbitrary bases of near-degenerate BFGS eigenspaces by deterministic ones.
635 : !> This leaves the existing involved-atom scoring and mode ordering unchanged.
636 : !> \param modes Orthonormal Hessian eigenvectors
637 : !> \param eigenvalues Ascending Hessian eigenvalues
638 : ! **************************************************************************************************
639 3 : SUBROUTINE bfgs_canonicalize_degenerate_modes(modes, eigenvalues)
640 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: modes
641 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: eigenvalues
642 :
643 : INTEGER :: first, last, n
644 : REAL(KIND=dp) :: cutoff, scale
645 :
646 3 : n = SIZE(eigenvalues)
647 3 : CPASSERT(SIZE(modes, 1) == n .AND. SIZE(modes, 2) == n)
648 3 : cutoff = (400.0_dp/vibfac)**2/massunit
649 :
650 3 : first = 1
651 396 : DO WHILE (first <= n)
652 393 : IF (eigenvalues(first) <= cutoff) THEN
653 124 : first = first + 1
654 124 : CYCLE
655 : END IF
656 : last = first
657 320 : DO WHILE (last < n)
658 317 : scale = MAX(ABS(eigenvalues(first)), ABS(eigenvalues(last + 1)))
659 317 : IF (eigenvalues(last + 1) - eigenvalues(first) > bfgs_guess_degeneracy_tol*scale) EXIT
660 269 : last = last + 1
661 : END DO
662 269 : IF (last > first) CALL bfgs_canonical_modes(modes(:, first:last))
663 269 : first = last + 1
664 : END DO
665 3 : END SUBROUTINE bfgs_canonicalize_degenerate_modes
666 :
667 : ! **************************************************************************************************
668 : !> \brief Fix a subspace basis using projections of Cartesian coordinate axes.
669 : !> \param modes Orthonormal basis, replaced by a deterministic orthonormal basis
670 : ! **************************************************************************************************
671 3 : SUBROUTINE bfgs_canonical_modes(modes)
672 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: modes
673 :
674 : INTEGER :: i, j, k, m, n, pass, pivot
675 : REAL(KIND=dp) :: maximum, norm
676 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: norms, projection, vector
677 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: residual, rotation
678 :
679 3 : n = SIZE(modes, 1)
680 3 : m = SIZE(modes, 2)
681 36 : ALLOCATE (residual(n, m), rotation(m, m), norms(n), projection(n), vector(m))
682 4917 : residual(:, :) = modes(:, :)
683 3 : rotation(:, :) = 0.0_dp
684 57 : DO j = 1, m
685 133434 : norms(:) = SUM(residual(:, :)**2, DIM=2)
686 4914 : maximum = MAXVAL(norms)
687 2329 : pivot = 1
688 2329 : DO i = 1, n
689 2329 : IF (norms(i) >= maximum*(1.0_dp - bfgs_guess_degeneracy_tol)) THEN
690 : pivot = i
691 : EXIT
692 : END IF
693 : END DO
694 1482 : vector(:) = residual(pivot, :)
695 162 : DO pass = 1, 2
696 1536 : DO k = 1, j - 1
697 80178 : vector(:) = vector(:) - DOT_PRODUCT(rotation(:, k), vector)*rotation(:, k)
698 : END DO
699 : END DO
700 1482 : norm = NORM2(vector)
701 54 : CPASSERT(norm > SQRT(EPSILON(1.0_dp)))
702 1482 : rotation(:, j) = vector(:)/norm
703 130002 : projection(:) = MATMUL(residual, rotation(:, j))
704 1485 : DO k = 1, m
705 130002 : residual(:, k) = residual(:, k) - projection(:)*rotation(k, j)
706 : END DO
707 : END DO
708 369 : residual(:, :) = MATMUL(modes, rotation)
709 4917 : modes(:, :) = residual(:, :)
710 3 : DEALLOCATE (residual, rotation, norms, projection, vector)
711 3 : END SUBROUTINE bfgs_canonical_modes
712 :
713 : ! **************************************************************************************************
714 : !> \brief ...
715 : !> \param ms_vib_section ...
716 : !> \param para_env ...
717 : !> \param ms_vib ...
718 : !> \param mass ...
719 : !> \param ionode ...
720 : !> \param particles ...
721 : !> \param nrep ...
722 : !> \param calc_intens ...
723 : !> \author Florian Schiffmann 11.2007
724 : ! **************************************************************************************************
725 4 : SUBROUTINE rest_guess(ms_vib_section, para_env, ms_vib, mass, ionode, particles, nrep, calc_intens)
726 :
727 : TYPE(section_vals_type), POINTER :: ms_vib_section
728 : TYPE(mp_para_env_type), POINTER :: para_env
729 : TYPE(ms_vib_type) :: ms_vib
730 : REAL(Kind=dp), DIMENSION(:) :: mass
731 : LOGICAL :: ionode
732 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
733 : INTEGER :: nrep
734 : LOGICAL :: calc_intens
735 :
736 : CHARACTER(LEN=default_path_length) :: ms_filename
737 : INTEGER :: hesunit, i, j, mat, natoms, ncoord, &
738 : output_unit, stat, statint
739 4 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
740 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval
741 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H
742 : TYPE(cp_logger_type), POINTER :: logger
743 :
744 4 : logger => cp_get_default_logger()
745 4 : output_unit = cp_logger_get_default_io_unit(logger)
746 :
747 4 : natoms = SIZE(particles)
748 4 : ncoord = 3*natoms
749 4 : IF (calc_intens) THEN
750 4 : DEALLOCATE (ms_vib%dip_deriv)
751 : END IF
752 :
753 4 : IF (ionode) THEN
754 :
755 2 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
756 2 : IF (ms_filename == "") ms_filename = "MS_RESTART"
757 : CALL open_file(file_name=ms_filename, &
758 : file_status="UNKNOWN", &
759 : file_form="UNFORMATTED", &
760 : file_action="READ", &
761 2 : unit_number=hesunit)
762 2 : READ (UNIT=hesunit, IOSTAT=stat) mat
763 2 : CPASSERT(stat == 0)
764 2 : ms_vib%mat_size = mat
765 : END IF
766 4 : CALL para_env%bcast(ms_vib%mat_size)
767 16 : ALLOCATE (ms_vib%b_mat(ncoord, ms_vib%mat_size))
768 12 : ALLOCATE (ms_vib%s_mat(ncoord, ms_vib%mat_size))
769 4 : IF (calc_intens) THEN
770 12 : ALLOCATE (ms_vib%dip_deriv(3, ms_vib%mat_size + nrep))
771 : END IF
772 4 : IF (ionode) THEN
773 2 : statint = 0
774 3384 : READ (UNIT=hesunit) ms_vib%b_mat
775 3384 : READ (UNIT=hesunit, IOSTAT=stat) ms_vib%s_mat
776 2 : IF (stat /= 0 .AND. output_unit > 0) THEN
777 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART **"
778 : END IF
779 2 : IF (calc_intens) THEN
780 2 : READ (UNIT=hesunit, IOSTAT=statint) ms_vib%dip_deriv(:, 1:ms_vib%mat_size)
781 2 : IF (statint /= 0 .AND. output_unit > 0) WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MS_RESTART,", &
782 0 : "intensities are requested but not present in restart file **"
783 : END IF
784 2 : CALL close_file(hesunit)
785 2 : IF (stat == 0 .AND. statint == 0 .AND. output_unit > 0) THEN
786 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** MS_RESTART has been read successfully ***"
787 : END IF
788 : END IF
789 13532 : CALL para_env%bcast(ms_vib%b_mat)
790 13532 : CALL para_env%bcast(ms_vib%s_mat)
791 340 : IF (calc_intens) CALL para_env%bcast(ms_vib%dip_deriv)
792 16 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
793 12 : ALLOCATE (eigenval(ms_vib%mat_size))
794 12 : ALLOCATE (ind(ms_vib%mat_size))
795 :
796 : 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), &
797 4 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
798 4 : CALL diamat_all(approx_H, eigenval)
799 :
800 4 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, ms_vib%b_vec)
801 4 : IF (ms_vib%initial_guess /= 4) THEN
802 :
803 716 : ms_vib%b_vec = 0._dp
804 8 : DO i = 1, nrep
805 42 : DO j = 1, ms_vib%mat_size
806 6768 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i) + approx_H(j, ind(i))*ms_vib%b_mat(:, j)
807 : END DO
808 1424 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/NORM2(ms_vib%b_vec(:, i))
809 : END DO
810 :
811 4 : DEALLOCATE (ms_vib%s_mat)
812 4 : DEALLOCATE (ms_vib%b_mat)
813 4 : IF (calc_intens) THEN
814 4 : DEALLOCATE (ms_vib%dip_deriv)
815 12 : ALLOCATE (ms_vib%dip_deriv(3, nrep))
816 : END IF
817 : END IF
818 4 : DEALLOCATE (approx_H)
819 4 : DEALLOCATE (eigenval)
820 4 : DEALLOCATE (ind)
821 8 : DO i = 1, nrep
822 716 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
823 : END DO
824 :
825 4 : END SUBROUTINE rest_guess
826 :
827 : ! **************************************************************************************************
828 : !> \brief ...
829 : !> \param ms_vib_section ...
830 : !> \param input ...
831 : !> \param para_env ...
832 : !> \param ms_vib ...
833 : !> \param mass ...
834 : !> \param ncoord ...
835 : !> \param nrep ...
836 : !> \param logger ...
837 : !> \author Florian Schiffmann 11.2007
838 : ! **************************************************************************************************
839 4 : SUBROUTINE molden_guess(ms_vib_section, input, para_env, ms_vib, mass, ncoord, nrep, logger)
840 : TYPE(section_vals_type), POINTER :: ms_vib_section, input
841 : TYPE(mp_para_env_type), POINTER :: para_env
842 : TYPE(ms_vib_type) :: ms_vib
843 : REAL(Kind=dp), DIMENSION(:) :: mass
844 : INTEGER :: ncoord, nrep
845 : TYPE(cp_logger_type), POINTER :: logger
846 :
847 : CHARACTER(LEN=2) :: at_name
848 : CHARACTER(LEN=default_path_length) :: ms_filename
849 : CHARACTER(LEN=max_line_length) :: info
850 : INTEGER :: i, istat, iw, j, jj, k, nvibs, &
851 : output_molden, output_unit, stat
852 4 : INTEGER, DIMENSION(:), POINTER :: tmplist
853 : LOGICAL :: reading_vib
854 : REAL(KIND=dp) :: my_val, norm
855 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: freq, tmp
856 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: modes
857 8 : REAL(KIND=dp), DIMENSION(3, ncoord/3) :: pos
858 :
859 8 : output_unit = cp_logger_get_default_io_unit(logger)
860 :
861 4 : CALL section_vals_val_get(ms_vib_section, "RESTART_FILE_NAME", c_val=ms_filename)
862 4 : IF (ms_filename == "") output_molden = &
863 : cp_print_key_unit_nr(logger, input, "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB", &
864 : extension=".mol", file_status='UNKNOWN', &
865 0 : file_action="READ")
866 4 : IF (para_env%is_source()) THEN
867 :
868 2 : IF (ms_filename == "") THEN
869 0 : iw = output_molden
870 : ELSE
871 : CALL open_file(file_name=TRIM(ms_filename), &
872 : file_status="UNKNOWN", &
873 : file_form="FORMATTED", &
874 : file_action="READ", &
875 2 : unit_number=iw)
876 : END IF
877 2 : info = ""
878 2 : READ (iw, *) info
879 2 : READ (iw, *) info
880 2 : istat = 0
881 2 : nvibs = 0
882 2 : reading_vib = .FALSE.
883 140 : DO
884 142 : READ (iw, *, IOSTAT=stat) info
885 142 : istat = istat + stat
886 142 : IF (TRIM(ADJUSTL(info)) == "[FR-COORD]") EXIT
887 :
888 140 : CPASSERT(stat == 0)
889 :
890 140 : IF (reading_vib) nvibs = nvibs + 1
891 140 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") reading_vib = .TRUE.
892 : END DO
893 2 : REWIND (iw)
894 2 : istat = 0
895 2 : READ (iw, *, IOSTAT=stat) info
896 2 : istat = istat + stat
897 2 : READ (iw, *, IOSTAT=stat) info
898 2 : istat = istat + stat
899 : ! Skip [Atoms] section
900 118 : DO
901 120 : READ (iw, *, IOSTAT=stat) info
902 120 : istat = istat + stat
903 120 : CPASSERT(stat == 0)
904 120 : IF (TRIM(ADJUSTL(info)) == "[FREQ]") EXIT
905 : END DO
906 : ! Read frequencies and modes
907 6 : ALLOCATE (freq(nvibs))
908 8 : ALLOCATE (modes(ncoord, nvibs))
909 :
910 22 : DO i = 1, nvibs
911 20 : READ (iw, *, IOSTAT=stat) freq(i)
912 22 : istat = istat + stat
913 : END DO
914 2 : READ (iw, *) info
915 120 : DO i = 1, ncoord/3
916 118 : READ (iw, *, IOSTAT=stat) at_name, pos(:, i)
917 120 : istat = istat + stat
918 : END DO
919 2 : READ (iw, *) info
920 22 : DO i = 1, nvibs
921 20 : READ (iw, *) info
922 20 : istat = istat + stat
923 1202 : DO j = 1, ncoord/3
924 1180 : k = (j - 1)*3 + 1
925 1180 : READ (iw, *, IOSTAT=stat) modes(k:k + 2, i)
926 1200 : istat = istat + stat
927 : END DO
928 : END DO
929 2 : IF (ms_filename /= "") CALL close_file(iw)
930 2 : IF (output_unit > 0) THEN
931 2 : IF (istat /= 0) THEN
932 0 : WRITE (output_unit, FMT="(/,T2,A)") "** Error while reading MOLDEN file **"
933 : ELSE
934 2 : WRITE (output_unit, FMT="(/,T2,A)") "*** MOLDEN file has been read successfully ***"
935 : END IF
936 : END IF
937 : !!!!!!! select modes !!!!!!
938 4 : ALLOCATE (tmp(nvibs))
939 2 : tmp(:) = 0.0_dp
940 6 : ALLOCATE (tmplist(nvibs))
941 2 : IF (ms_vib%select_id == 1) my_val = ms_vib%sel_freq
942 2 : IF (ms_vib%select_id == 2) my_val = (ms_vib%f_range(2) + ms_vib%f_range(1))*0.5_dp
943 2 : IF (ms_vib%select_id == 1 .OR. ms_vib%select_id == 2) THEN
944 11 : DO i = 1, nvibs
945 11 : tmp(i) = ABS(my_val - freq(i))
946 : END DO
947 1 : ELSE IF (ms_vib%select_id == 3) THEN
948 11 : DO i = 1, nvibs
949 30 : DO j = 1, SIZE(ms_vib%inv_atoms)
950 90 : DO k = 1, 3
951 60 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
952 80 : tmp(i) = tmp(i) + SQRT(modes(jj, i)**2)
953 : END DO
954 : END DO
955 11 : IF (freq(i) <= 400._dp) tmp(i) = 0._dp
956 : END DO
957 11 : tmp(:) = -tmp(:)
958 : END IF
959 2 : CALL sort(tmp, nvibs, tmplist)
960 4 : DO i = 1, nrep
961 356 : ms_vib%b_vec(:, i) = modes(:, tmplist(i))*mass(:)
962 356 : norm = NORM2(ms_vib%b_vec(:, i))
963 358 : ms_vib%b_vec(:, i) = ms_vib%b_vec(:, i)/norm
964 : END DO
965 4 : DO i = 1, nrep
966 358 : ms_vib%delta_vec(:, i) = ms_vib%b_vec(:, i)/mass(:)
967 : END DO
968 :
969 2 : DEALLOCATE (freq)
970 2 : DEALLOCATE (modes)
971 2 : DEALLOCATE (tmp)
972 2 : DEALLOCATE (tmplist)
973 :
974 : END IF
975 1428 : CALL para_env%bcast(ms_vib%b_vec)
976 1428 : CALL para_env%bcast(ms_vib%delta_vec)
977 :
978 4 : IF (ms_filename == "") CALL cp_print_key_finished_output(output_molden, logger, input, &
979 0 : "VIBRATIONAL_ANALYSIS%PRINT%MOLDEN_VIB")
980 4 : END SUBROUTINE molden_guess
981 :
982 : ! **************************************************************************************************
983 : !> \brief Davidson algorithm for to generate a approximate Hessian for mode
984 : !> selective vibrational analysis
985 : !> \param rep_env ...
986 : !> \param ms_vib ...
987 : !> \param input ...
988 : !> \param nrep ...
989 : !> \param particles ...
990 : !> \param mass ...
991 : !> \param converged ...
992 : !> \param dx ...
993 : !> \param calc_intens ...
994 : !> \param output_unit_ms ...
995 : !> \param logger ...
996 : !> \param cell simulation cell
997 : !> \author Florian Schiffmann 11.2007
998 : ! **************************************************************************************************
999 186 : SUBROUTINE evaluate_H_update_b(rep_env, ms_vib, input, nrep, &
1000 : particles, &
1001 186 : mass, &
1002 : converged, dx, &
1003 : calc_intens, output_unit_ms, logger, cell)
1004 : TYPE(replica_env_type), POINTER :: rep_env
1005 : TYPE(ms_vib_type) :: ms_vib
1006 : TYPE(section_vals_type), POINTER :: input
1007 : INTEGER :: nrep
1008 : TYPE(particle_type), DIMENSION(:), POINTER :: particles
1009 : REAL(Kind=dp), DIMENSION(:) :: mass
1010 : LOGICAL :: converged
1011 : REAL(KIND=dp) :: dx
1012 : LOGICAL :: calc_intens
1013 : INTEGER :: output_unit_ms
1014 : TYPE(cp_logger_type), POINTER :: logger
1015 : TYPE(cell_type), POINTER :: cell
1016 :
1017 : INTEGER :: i, j, jj, k, natoms, ncoord
1018 186 : INTEGER, ALLOCATABLE, DIMENSION(:) :: ind
1019 : LOGICAL :: dump_only_positive
1020 186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: eigenval, freq
1021 186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: approx_H, H_save, residuum, tmp_b, tmp_s
1022 372 : REAL(KIND=dp), DIMENSION(2, nrep) :: criteria
1023 186 : REAL(Kind=dp), DIMENSION(:), POINTER :: intensities
1024 :
1025 186 : natoms = SIZE(particles)
1026 186 : ncoord = 3*natoms
1027 186 : nrep = SIZE(rep_env%f, 2)
1028 :
1029 : !!!!!!!! reallocate and update the davidson matrices !!!!!!!!!!
1030 186 : IF (ms_vib%mat_size /= 0) THEN
1031 :
1032 790 : ALLOCATE (tmp_b(3*natoms, ms_vib%mat_size))
1033 474 : ALLOCATE (tmp_s(3*natoms, ms_vib%mat_size))
1034 :
1035 167174 : tmp_b(:, :) = ms_vib%b_mat
1036 167174 : tmp_s(:, :) = ms_vib%s_mat
1037 :
1038 158 : DEALLOCATE (ms_vib%b_mat)
1039 158 : DEALLOCATE (ms_vib%s_mat)
1040 : END IF
1041 :
1042 930 : ALLOCATE (ms_vib%b_mat(3*natoms, ms_vib%mat_size + nrep))
1043 558 : ALLOCATE (ms_vib%s_mat(3*natoms, ms_vib%mat_size + nrep))
1044 :
1045 193968 : ms_vib%s_mat = 0.0_dp
1046 :
1047 26646 : DO i = 1, 3*natoms
1048 26460 : IF (ms_vib%mat_size /= 0) THEN
1049 189354 : DO j = 1, ms_vib%mat_size
1050 166014 : ms_vib%b_mat(i, j) = tmp_b(i, j)
1051 189354 : ms_vib%s_mat(i, j) = tmp_s(i, j)
1052 : END DO
1053 : END IF
1054 53214 : DO j = 1, nrep
1055 53028 : ms_vib%b_mat(i, ms_vib%mat_size + j) = ms_vib%b_vec(i, j)
1056 : END DO
1057 : END DO
1058 :
1059 186 : IF (ms_vib%mat_size /= 0) THEN
1060 158 : DEALLOCATE (tmp_s)
1061 158 : DEALLOCATE (tmp_b)
1062 : END IF
1063 :
1064 186 : ms_vib%mat_size = ms_vib%mat_size + nrep
1065 :
1066 744 : ALLOCATE (approx_H(ms_vib%mat_size, ms_vib%mat_size))
1067 558 : ALLOCATE (H_save(ms_vib%mat_size, ms_vib%mat_size))
1068 558 : ALLOCATE (eigenval(ms_vib%mat_size))
1069 :
1070 : !!!!!!!!!!!! calculate the new derivativ and the approximate hessian
1071 :
1072 384 : DO i = 1, nrep
1073 26952 : DO j = 1, 3*natoms
1074 26766 : 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))
1075 : END DO
1076 : END DO
1077 :
1078 : 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), &
1079 186 : ms_vib%s_mat, SIZE(ms_vib%s_mat, 1), 0._dp, approx_H, ms_vib%mat_size)
1080 14730 : H_save(:, :) = approx_H
1081 :
1082 186 : CALL diamat_all(approx_H, eigenval)
1083 :
1084 : !!!!!!!!!!!! select eigenvalue(s) and vector(s) and calculate the new displacement vector
1085 558 : ALLOCATE (ind(ms_vib%mat_size))
1086 744 : ALLOCATE (residuum(SIZE(ms_vib%s_mat, 1), nrep))
1087 :
1088 186 : CALL select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1089 :
1090 384 : DO i = 1, nrep
1091 9240 : DO j = 1, natoms
1092 35622 : DO k = 1, 3
1093 26568 : jj = (j - 1)*3 + k
1094 35424 : ms_vib%delta_vec(jj, i) = ms_vib%b_vec(jj, i)/mass(jj)
1095 : END DO
1096 : END DO
1097 : END DO
1098 :
1099 384 : DO i = 1, nrep
1100 26766 : ms_vib%step_r(i) = dx/NORM2(ms_vib%delta_vec(:, i))
1101 26952 : ms_vib%step_b(i) = NORM2(ms_vib%step_r(i)*ms_vib%b_vec(:, i))
1102 : END DO
1103 186 : converged = .FALSE.
1104 : IF (MAXVAL(criteria(1, :)) <= ms_vib%eps(1) .AND. MAXVAL(criteria(2, :)) &
1105 582 : <= ms_vib%eps(2) .OR. ms_vib%mat_size >= ncoord) converged = .TRUE.
1106 558 : ALLOCATE (freq(nrep))
1107 384 : DO i = 1, nrep
1108 384 : freq(i) = SQRT(ABS(eigenval(ind(i)))*massunit)*vibfac
1109 : END DO
1110 :
1111 : !!! write information and output !!!
1112 186 : IF (converged) THEN
1113 226 : eigenval(:) = SIGN(1._dp, eigenval(:))*SQRT(ABS(eigenval(:))*massunit)*vibfac
1114 112 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1115 28 : tmp_b = 0._dp
1116 84 : ALLOCATE (tmp_s(3, ms_vib%mat_size))
1117 28 : tmp_s = 0._dp
1118 28 : IF (calc_intens) THEN
1119 66 : ALLOCATE (intensities(ms_vib%mat_size))
1120 190 : intensities = 0._dp
1121 : END IF
1122 226 : DO i = 1, ms_vib%mat_size
1123 2472 : DO j = 1, ms_vib%mat_size
1124 361716 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
1125 : END DO
1126 53362 : tmp_b(:, i) = tmp_b(:, i)/NORM2(tmp_b(:, i))
1127 : END DO
1128 28 : IF (calc_intens) THEN
1129 190 : DO i = 1, ms_vib%mat_size
1130 2280 : DO j = 1, ms_vib%mat_size
1131 8616 : tmp_s(:, i) = tmp_s(:, i) + ms_vib%dip_deriv(:, j)*approx_H(j, i)
1132 : END DO
1133 694 : IF (calc_intens) intensities(i) = NORM2(tmp_s(:, i))
1134 : END DO
1135 : END IF
1136 28 : IF (calc_intens) THEN
1137 : CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
1138 : input, nrep, approx_H, eigenval, calc_intens, &
1139 22 : intensities=intensities, logger=logger)
1140 : ELSE
1141 : CALL ms_out(output_unit_ms, converged, freq, criteria, ms_vib, &
1142 6 : input, nrep, approx_H, eigenval, calc_intens, logger=logger)
1143 : END IF
1144 28 : dump_only_positive = ms_vib%low_freq > 0.0_dp
1145 : CALL write_vibrations_molden(input, particles, eigenval, tmp_b, intensities, calc_intens, &
1146 28 : dump_only_positive=dump_only_positive, logger=logger, cell=cell)
1147 28 : IF (calc_intens) THEN
1148 22 : DEALLOCATE (intensities)
1149 : END IF
1150 28 : DEALLOCATE (tmp_b)
1151 28 : DEALLOCATE (tmp_s)
1152 : END IF
1153 :
1154 186 : IF (.NOT. converged) CALL ms_out(output_unit_ms, converged, freq, criteria, &
1155 158 : ms_vib, input, nrep, approx_H, eigenval, calc_intens, logger=logger)
1156 :
1157 186 : DEALLOCATE (freq)
1158 186 : DEALLOCATE (approx_H)
1159 186 : DEALLOCATE (eigenval)
1160 186 : DEALLOCATE (residuum)
1161 186 : DEALLOCATE (ind)
1162 :
1163 372 : END SUBROUTINE evaluate_H_update_b
1164 :
1165 : ! **************************************************************************************************
1166 : !> \brief writes the output for a mode tracking calculation
1167 : !> \param ms_vib ...
1168 : !> \param nrep ...
1169 : !> \param mass ...
1170 : !> \param ncoord ...
1171 : !> \param approx_H ...
1172 : !> \param eigenval ...
1173 : !> \param ind ...
1174 : !> \param residuum ...
1175 : !> \param criteria ...
1176 : !> \author Florian Schiffmann 11.2007
1177 : ! **************************************************************************************************
1178 190 : SUBROUTINE select_vector(ms_vib, nrep, mass, ncoord, approx_H, eigenval, ind, residuum, criteria)
1179 :
1180 : TYPE(ms_vib_type) :: ms_vib
1181 : INTEGER :: nrep
1182 : REAL(Kind=dp), DIMENSION(:) :: mass
1183 : INTEGER :: ncoord
1184 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1185 : REAL(Kind=dp), DIMENSION(:) :: eigenval
1186 : INTEGER, DIMENSION(:) :: ind
1187 : REAL(KIND=dp), DIMENSION(:, :) :: residuum
1188 : REAL(KIND=dp), DIMENSION(2, nrep), OPTIONAL :: criteria
1189 :
1190 : INTEGER :: i, j, jj, k
1191 : REAL(KIND=dp) :: my_val, norm
1192 190 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp
1193 190 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_b
1194 :
1195 570 : ALLOCATE (tmp(ms_vib%mat_size))
1196 :
1197 306 : SELECT CASE (ms_vib%select_id)
1198 : CASE (1)
1199 116 : my_val = (ms_vib%sel_freq/(vibfac))**2/massunit
1200 1006 : DO i = 1, ms_vib%mat_size
1201 1006 : tmp(i) = ABS(my_val - eigenval(i))
1202 : END DO
1203 116 : CALL sort(tmp, (ms_vib%mat_size), ind)
1204 15892 : residuum = 0._dp
1205 238 : DO j = 1, nrep
1206 1152 : DO i = 1, ms_vib%mat_size
1207 144040 : residuum(:, j) = residuum(:, j) + approx_H(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1208 : END DO
1209 : END DO
1210 : CASE (2)
1211 6 : CALL get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1212 : CASE (3)
1213 :
1214 272 : ALLOCATE (tmp_b(ncoord, ms_vib%mat_size))
1215 68 : tmp_b = 0._dp
1216 :
1217 392 : DO i = 1, ms_vib%mat_size
1218 2428 : DO j = 1, ms_vib%mat_size
1219 372400 : tmp_b(:, i) = tmp_b(:, i) + approx_H(j, i)*ms_vib%b_mat(:, j)/mass(:)
1220 : END DO
1221 113000 : tmp_b(:, i) = tmp_b(:, i)/NORM2(tmp_b(:, i))
1222 : END DO
1223 68 : tmp = 0._dp
1224 392 : DO i = 1, ms_vib%mat_size
1225 972 : DO j = 1, SIZE(ms_vib%inv_atoms)
1226 2916 : DO k = 1, 3
1227 1944 : jj = (ms_vib%inv_atoms(j) - 1)*3 + k
1228 2592 : tmp(i) = tmp(i) + SQRT(tmp_b(jj, i)**2)
1229 : END DO
1230 : END DO
1231 392 : IF (.NOT. ASSOCIATED(ms_vib%inv_range)) THEN
1232 324 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) <= 400._dp) tmp(i) = 0._dp
1233 : ELSE
1234 0 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) <= ms_vib%inv_range(1)) tmp(i) = 0._dp
1235 0 : IF ((SIGN(1._dp, eigenval(i))*SQRT(ABS(eigenval(i))*massunit)*vibfac) >= ms_vib%inv_range(2)) tmp(i) = 0._dp
1236 : END IF
1237 : END DO
1238 392 : tmp(:) = -tmp(:)
1239 68 : CALL sort(tmp, (ms_vib%mat_size), ind)
1240 11650 : residuum(:, :) = 0._dp
1241 :
1242 136 : DO j = 1, nrep
1243 460 : DO i = 1, ms_vib%mat_size
1244 56696 : residuum(:, j) = residuum(:, j) + approx_H(i, ind(j))*(ms_vib%s_mat(:, i) - eigenval(ind(j))*ms_vib%b_mat(:, i))
1245 : END DO
1246 : END DO
1247 258 : DEALLOCATE (tmp_b)
1248 : END SELECT
1249 :
1250 392 : DO j = 1, nrep
1251 1678 : DO i = 1, ms_vib%mat_size
1252 400968 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1253 : END DO
1254 : END DO
1255 190 : IF (PRESENT(criteria)) THEN
1256 384 : DO i = 1, nrep
1257 26766 : criteria(1, i) = MAXVAL((residuum(:, i)))
1258 26952 : criteria(2, i) = NORM2(residuum(:, i))
1259 : END DO
1260 : END IF
1261 :
1262 392 : DO i = 1, nrep
1263 27478 : norm = NORM2(residuum(:, i))
1264 27668 : residuum(:, i) = residuum(:, i)/norm
1265 : END DO
1266 :
1267 2090 : DO k = 1, 10
1268 4110 : DO j = 1, nrep
1269 14880 : DO i = 1, ms_vib%mat_size
1270 4007660 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1271 4009680 : residuum(:, j) = residuum(:, j)/NORM2(residuum(:, j))
1272 : END DO
1273 3920 : IF (nrep > 1) THEN
1274 720 : DO i = 1, nrep
1275 720 : IF (i /= j) THEN
1276 4560 : residuum(:, j) = residuum(:, j) - DOT_PRODUCT(residuum(:, j), residuum(:, i))*residuum(:, i)
1277 4560 : residuum(:, j) = residuum(:, j)/NORM2(residuum(:, j))
1278 : END IF
1279 : END DO
1280 : END IF
1281 : END DO
1282 : END DO
1283 27668 : ms_vib%b_vec = residuum
1284 190 : DEALLOCATE (tmp)
1285 190 : END SUBROUTINE select_vector
1286 :
1287 : ! **************************************************************************************************
1288 : !> \brief writes the output for a mode tracking calculation
1289 : !> \param iw ...
1290 : !> \param converged ...
1291 : !> \param freq ...
1292 : !> \param criter ...
1293 : !> \param ms_vib ...
1294 : !> \param input ...
1295 : !> \param nrep ...
1296 : !> \param approx_H ...
1297 : !> \param eigenval ...
1298 : !> \param calc_intens ...
1299 : !> \param intensities ...
1300 : !> \param logger ...
1301 : !> \author Florian Schiffmann 11.2007
1302 : ! **************************************************************************************************
1303 186 : SUBROUTINE ms_out(iw, converged, freq, criter, ms_vib, input, nrep, &
1304 186 : approx_H, eigenval, calc_intens, intensities, logger)
1305 :
1306 : INTEGER :: iw
1307 : LOGICAL :: converged
1308 : REAL(KIND=dp), DIMENSION(:) :: freq
1309 : REAL(KIND=dp), DIMENSION(:, :) :: criter
1310 : TYPE(ms_vib_type) :: ms_vib
1311 : TYPE(section_vals_type), POINTER :: input
1312 : INTEGER :: nrep
1313 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1314 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1315 : LOGICAL :: calc_intens
1316 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: intensities
1317 : TYPE(cp_logger_type), POINTER :: logger
1318 :
1319 : REAL(KIND=dp), PARAMETER :: eps_criter1 = 1.0E-07_dp, &
1320 : eps_criter2 = 1.0E-06_dp
1321 :
1322 : INTEGER :: i, j, msunit
1323 : REAL(KIND=dp) :: crit_a, crit_b, fint, gintval
1324 186 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: residuum
1325 : TYPE(section_vals_type), POINTER :: ms_vib_section
1326 :
1327 : ms_vib_section => section_vals_get_subs_vals(input, &
1328 186 : "VIBRATIONAL_ANALYSIS%MODE_SELECTIVE")
1329 :
1330 186 : fint = 42.255_dp*massunit*debye**2*bohr**2
1331 :
1332 186 : IF (converged) THEN
1333 28 : IF (iw > 0) THEN
1334 14 : WRITE (iw, '(T2,A)') "MS| DAVIDSON ALGORITHM CONVERGED"
1335 30 : DO i = 1, nrep
1336 30 : WRITE (iw, '(T2,"MS| TRACKED FREQUENCY (",I0,") IS:",F12.6,3X,A)') i, freq(i), 'cm-1'
1337 : END DO
1338 42 : ALLOCATE (residuum(SIZE(ms_vib%b_mat, 1)))
1339 14 : WRITE (iw, '( /, 1X, 79("-") )')
1340 14 : WRITE (iw, '( 25X, A)') 'FREQUENCY AND CONVERGENCE LIST'
1341 14 : IF (PRESENT(intensities)) THEN
1342 11 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'INT[KM/Mole]', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1343 : ELSE
1344 3 : WRITE (iw, '(3X,5(4X, A))') 'FREQUENCY', 'MAXVAL CRITERIA', 'NORM CRITERIA', 'CONVERGENCE'
1345 : END IF
1346 113 : DO i = 1, SIZE(ms_vib%b_mat, 2)
1347 99 : residuum = 0._dp
1348 1236 : DO j = 1, SIZE(ms_vib%b_mat, 2)
1349 180858 : residuum(:) = residuum(:) + approx_H(j, i)*(ms_vib%s_mat(:, j) - eigenval(i)*ms_vib%b_mat(:, j))
1350 : END DO
1351 1236 : DO j = 1, ms_vib%mat_size
1352 360480 : residuum(:) = residuum(:) - DOT_PRODUCT(residuum(:), ms_vib%b_mat(:, j))*ms_vib%b_mat(:, j)
1353 : END DO
1354 13383 : crit_a = MAXVAL(residuum(:))
1355 13383 : crit_b = NORM2(residuum)
1356 113 : IF (PRESENT(intensities)) THEN
1357 84 : gintval = fint*intensities(i)**2
1358 84 : IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1359 30 : 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)') &
1360 28 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'YES'
1361 : ELSE
1362 54 : 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)') &
1363 54 : 'VIB|', eigenval(i), gintval, crit_a, crit_b, 'NO'
1364 : END IF
1365 : ELSE
1366 15 : IF (crit_a <= ms_vib%eps(1) .AND. crit_b <= ms_vib%eps(2)) THEN
1367 14 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1368 8 : 'VIB|', eigenval(i), crit_a, crit_b, 'YES'
1369 : ELSE
1370 1 : IF (eigenval(i) > ms_vib%low_freq) WRITE (iw, '(2X,A,2X,F9.3,5X,E12.6,5X,E12.3,11X,A)') &
1371 1 : 'VIB|', eigenval(i), crit_a, crit_b, 'NO'
1372 : END IF
1373 : END IF
1374 : END DO
1375 14 : DEALLOCATE (residuum)
1376 :
1377 : msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1378 : "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1379 : file_status="REPLACE", file_form="UNFORMATTED", &
1380 14 : file_action="WRITE")
1381 :
1382 14 : IF (msunit > 0) THEN
1383 14 : WRITE (UNIT=msunit) ms_vib%mat_size
1384 13397 : WRITE (UNIT=msunit) ms_vib%b_mat
1385 13397 : WRITE (UNIT=msunit) ms_vib%s_mat
1386 350 : IF (calc_intens) WRITE (UNIT=msunit) ms_vib%dip_deriv
1387 : END IF
1388 :
1389 : CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1390 14 : "PRINT%MS_RESTART")
1391 : END IF
1392 : ELSE
1393 158 : IF (iw > 0) THEN
1394 : msunit = cp_print_key_unit_nr(logger, ms_vib_section, &
1395 : "PRINT%MS_RESTART", extension=".bin", middle_name="MS_RESTART", &
1396 : file_status="REPLACE", file_form="UNFORMATTED", &
1397 79 : file_action="WRITE")
1398 :
1399 79 : IF (msunit > 0) THEN
1400 79 : WRITE (UNIT=msunit) ms_vib%mat_size
1401 83587 : WRITE (UNIT=msunit) ms_vib%b_mat
1402 83587 : WRITE (UNIT=msunit) ms_vib%s_mat
1403 2023 : IF (calc_intens) WRITE (UNIT=msunit) ms_vib%dip_deriv
1404 : END IF
1405 :
1406 : CALL cp_print_key_finished_output(msunit, logger, ms_vib_section, &
1407 79 : "PRINT%MS_RESTART")
1408 :
1409 79 : WRITE (iw, '(T2,A,3X,I6)') "MS| ITERATION STEP", ms_vib%mat_size/nrep
1410 162 : DO i = 1, nrep
1411 162 : IF (criter(1, i) <= eps_criter1 .AND. (criter(2, i)) <= eps_criter2) THEN
1412 1 : WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 IS CONVERGED"
1413 : ELSE
1414 82 : WRITE (iw, '(T2,A,3X,F12.6,A)') "MS| TRACKED MODE ", freq(i), "cm-1 NOT CONVERGED"
1415 : END IF
1416 : END DO
1417 : END IF
1418 : END IF
1419 :
1420 186 : END SUBROUTINE ms_out
1421 :
1422 : ! **************************************************************************************************
1423 : !> \brief ...
1424 : !> \param ms_vib ...
1425 : !> \param approx_H ...
1426 : !> \param eigenval ...
1427 : !> \param residuum ...
1428 : !> \param nrep ...
1429 : !> \param ind ...
1430 : !> \author Florian Schiffmann 11.2007
1431 : ! **************************************************************************************************
1432 6 : SUBROUTINE get_vibs_in_range(ms_vib, approx_H, eigenval, residuum, nrep, ind)
1433 :
1434 : TYPE(ms_vib_type) :: ms_vib
1435 : REAL(KIND=dp), DIMENSION(:, :) :: approx_H
1436 : REAL(KIND=dp), DIMENSION(:) :: eigenval
1437 : REAL(KIND=dp), DIMENSION(:, :) :: residuum
1438 : INTEGER :: nrep
1439 : INTEGER, DIMENSION(:) :: ind
1440 :
1441 : INTEGER :: count1, count2, i, j
1442 6 : INTEGER, ALLOCATABLE, DIMENSION(:) :: map2
1443 6 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: map1
1444 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tmp, tmp1
1445 6 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tmp_resid
1446 : REAL(KIND=dp), DIMENSION(2) :: myrange
1447 :
1448 18 : myrange(:) = (ms_vib%f_range(:)/(vibfac))**2/massunit
1449 6 : count1 = 0
1450 6 : count2 = 0
1451 126 : residuum = 0.0_dp
1452 6 : ms_vib%mat_size = SIZE(ms_vib%b_mat, 2)
1453 18 : ALLOCATE (map1(SIZE(eigenval), 2))
1454 18 : ALLOCATE (tmp(SIZE(eigenval)))
1455 30 : DO i = 1, SIZE(eigenval)
1456 24 : IF (ABS(eigenval(i) - myrange(1)) + ABS(eigenval(i) - myrange(2)) <= &
1457 6 : ABS(myrange(1) - myrange(2)) + myrange(1)*0.001_dp) THEN
1458 0 : count1 = count1 + 1
1459 0 : map1(count1, 1) = i
1460 : ELSE
1461 24 : count2 = count2 + 1
1462 24 : map1(count2, 2) = i
1463 24 : tmp(count2) = MIN(ABS(eigenval(i) - myrange(1)), ABS(eigenval(i) - myrange(2)))
1464 : END IF
1465 : END DO
1466 :
1467 6 : IF (count1 == nrep) THEN
1468 0 : DO j = 1, count1
1469 0 : DO i = 1, ms_vib%mat_size
1470 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))
1471 0 : ind(j) = map1(j, 1)
1472 : END DO
1473 : END DO
1474 6 : ELSE IF (count1 > nrep) THEN
1475 0 : ALLOCATE (tmp_resid(SIZE(ms_vib%b_mat, 1), count1))
1476 0 : ALLOCATE (tmp1(count1))
1477 0 : ALLOCATE (map2(count1))
1478 0 : tmp_resid = 0._dp
1479 0 : DO j = 1, count1
1480 0 : DO i = 1, ms_vib%mat_size
1481 : tmp_resid(:, j) = tmp_resid(:, j) + approx_H(i, map1(j, 1))* &
1482 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1483 : END DO
1484 : END DO
1485 :
1486 0 : DO j = 1, count1
1487 0 : DO i = 1, ms_vib%mat_size
1488 0 : tmp_resid(:, j) = tmp_resid(:, j) - DOT_PRODUCT(tmp_resid(:, j), ms_vib%b_mat(:, i))*ms_vib%b_mat(:, i)
1489 : END DO
1490 0 : tmp(j) = MAXVAL(tmp_resid(:, j))
1491 : END DO
1492 0 : CALL sort(tmp, count1, map2)
1493 0 : DO j = 1, nrep
1494 0 : residuum(:, j) = tmp_resid(:, map2(count1 + 1 - j))
1495 0 : ind(j) = map1(map2(count1 + 1 - j), 1)
1496 : END DO
1497 0 : DEALLOCATE (tmp_resid)
1498 0 : DEALLOCATE (tmp1)
1499 0 : DEALLOCATE (map2)
1500 6 : ELSE IF (count1 < nrep) THEN
1501 :
1502 18 : ALLOCATE (map2(count2))
1503 6 : IF (count1 /= 0) THEN
1504 0 : DO j = 1, count1
1505 0 : DO i = 1, ms_vib%mat_size
1506 : residuum(:, j) = residuum(:, j) + approx_H(i, map1(j, 1))* &
1507 0 : (ms_vib%s_mat(:, i) - eigenval(map1(j, 1))*ms_vib%b_mat(:, i))
1508 : END DO
1509 0 : ind(j) = map1(j, 1)
1510 : END DO
1511 : END IF
1512 6 : CALL sort(tmp, count2, map2)
1513 18 : DO j = 1, nrep - count1
1514 60 : DO i = 1, ms_vib%mat_size
1515 : residuum(:, count1 + j) = residuum(:, count1 + j) + approx_H(i, map1(map2(j), 2)) &
1516 492 : *(ms_vib%s_mat(:, i) - eigenval(map1(map2(j), 2))*ms_vib%b_mat(:, i))
1517 : END DO
1518 18 : ind(count1 + j) = map1(map2(j), 2)
1519 : END DO
1520 :
1521 6 : DEALLOCATE (map2)
1522 : END IF
1523 :
1524 6 : DEALLOCATE (map1)
1525 6 : DEALLOCATE (tmp)
1526 :
1527 6 : END SUBROUTINE get_vibs_in_range
1528 0 : END MODULE mode_selective
|
