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