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