(git:34ef472)
graph_methods.F
Go to the documentation of this file.
1 !--------------------------------------------------------------------------------------------------!
2 ! CP2K: A general program to perform molecular dynamics simulations !
3 ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 ! !
5 ! SPDX-License-Identifier: GPL-2.0-or-later !
6 !--------------------------------------------------------------------------------------------------!
7 
8 ! **************************************************************************************************
9 !> \brief Program to Map on grid the hills spawned during a metadynamics run
10 !> \author Teodoro Laino [tlaino] - 06.2009
11 !> \par History
12 !> 03.2006 created [tlaino]
13 !> teodoro.laino .at. gmail.com
14 !> 11.2007 - tlaino (University of Zurich): Periodic COLVAR - cleaning.
15 !> 12.2010 - teodoro.laino@gmail.com: addition of the MEP for FES
16 !>
17 !> \par Note
18 !> Please report any bug to the author
19 ! **************************************************************************************************
20 MODULE graph_methods
21 
22  USE cp_files, ONLY: close_file,&
23  open_file
24  USE graph_utils, ONLY: derivative,&
25  get_val_res,&
26  mep_input_data_type,&
27  pbc,&
28  point_no_pbc,&
29  point_pbc
30  USE kinds, ONLY: dp
31  USE memory_utilities, ONLY: reallocate
32  USE periodic_table, ONLY: init_periodic_table,&
33  nelem,&
34  ptable
35  USE physcon, ONLY: bohr
36  USE string_utilities, ONLY: uppercase
37 #include "../base/base_uses.f90"
38 
39  IMPLICIT NONE
40  PRIVATE
41 
42  PUBLIC :: fes_compute_low, &
43  fes_write, &
44  fes_only_write, &
45  fes_min, &
46  fes_path, &
47  fes_cube_write
48 
49 CONTAINS
50 ! **************************************************************************************************
51 !> \brief Efficiently map the gaussians on the grid
52 !> \param idim ...
53 !> \param nn ...
54 !> \param fes ...
55 !> \param gauss ...
56 !> \param ind ...
57 !> \param ind0 ...
58 !> \param nfes ...
59 !> \param ndim ...
60 !> \param ngauss ...
61 !> \param ngrid ...
62 !> \param iperd ...
63 !> \par History
64 !> 03.2006 created [tlaino]
65 !> teodoro.laino .at. gmail.com
66 !> \author Teodoro Laino
67 ! **************************************************************************************************
68  RECURSIVE SUBROUTINE fes_compute_low(idim, nn, fes, gauss, ind, ind0, nfes, ndim, &
69  ngauss, ngrid, iperd)
70  INTEGER, INTENT(in) :: idim
71  INTEGER, DIMENSION(:) :: nn
72  REAL(kind=dp), DIMENSION(:), POINTER :: fes
73  REAL(kind=dp), DIMENSION(:, :), POINTER :: gauss
74  INTEGER, DIMENSION(:) :: ind, ind0
75  INTEGER, INTENT(in) :: nfes, ndim, ngauss
76  INTEGER, DIMENSION(:), POINTER :: ngrid
77  INTEGER, DIMENSION(:) :: iperd
78 
79  INTEGER :: i, j, k, pnt
80  INTEGER, DIMENSION(:), POINTER :: ll, pos
81  REAL(kind=dp) :: prod
82 
83  ALLOCATE (pos(ndim), ll(ndim))
84  pos = ind
85  k = nn(idim)
86 
87  DO i = -k, k
88  pos(idim) = ind(idim) + i
89  IF (iperd(idim) == 0) THEN
90  IF (pos(idim) .GT. ngrid(idim)) cycle
91  IF (pos(idim) .LT. 1) cycle
92  END IF
93  IF (idim /= 1) THEN
94  CALL fes_compute_low(idim - 1, nn, fes, gauss, pos, ind0, nfes, ndim, ngauss, ngrid, iperd)
95  ELSE
96  pnt = point_pbc(pos, iperd, ngrid, ndim)
97  prod = 1.0_dp
98  DO j = 1, ndim
99  ll(j) = pos(j) - ind0(j)
100  prod = prod*gauss(ll(j), j)
101  END DO
102  fes(pnt) = fes(pnt) + prod
103  END IF
104  END DO
105  DEALLOCATE (pos, ll)
106 
107  END SUBROUTINE fes_compute_low
108 
109 ! **************************************************************************************************
110 !> \brief Writes the FES on the file
111 !> \param unit_nr ...
112 !> \param idim ...
113 !> \param fes ...
114 !> \param pos ...
115 !> \param ndim ...
116 !> \param ngrid ...
117 !> \param dp_grid ...
118 !> \param x0 ...
119 !> \param ndw ...
120 !> \param l_fes_int ...
121 !> \param array ...
122 !> \par History
123 !> 03.2006 created [tlaino]
124 !> teodoro.laino .at. gmail.com
125 !> \author Teodoro Laino
126 ! **************************************************************************************************
127  RECURSIVE SUBROUTINE fes_write(unit_nr, idim, fes, pos, ndim, ngrid, &
128  dp_grid, x0, ndw, l_fes_int, array)
129  INTEGER, INTENT(IN) :: unit_nr, idim
130  REAL(kind=dp), DIMENSION(:), POINTER :: fes
131  INTEGER, DIMENSION(:), POINTER :: pos
132  INTEGER, INTENT(IN) :: ndim
133  INTEGER, DIMENSION(:), POINTER :: ngrid
134  REAL(kind=dp), DIMENSION(:), POINTER :: dp_grid, x0
135  INTEGER, INTENT(IN) :: ndw
136  LOGICAL, INTENT(IN) :: l_fes_int
137  REAL(kind=dp), DIMENSION(:), OPTIONAL :: array
138 
139  INTEGER :: dimval, i, id, ind, is, it, itt, np, pnt
140  REAL(kind=dp) :: dvol, sum_fes
141  REAL(kind=dp), DIMENSION(:), POINTER :: xx
142 
143  ALLOCATE (xx(ndim))
144  xx = x0
145  DO i = 1, ngrid(idim)
146  pos(idim) = i
147  IF (idim /= ndim - ndw + 1) THEN
148  IF (PRESENT(array)) THEN
149  CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, &
150  x0, ndw, l_fes_int, array)
151  ELSE
152  CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, &
153  x0, ndw, l_fes_int)
154  END IF
155  ELSE
156  IF (PRESENT(array)) THEN
157  ind = 1
158  np = ngrid(ndim)*ngrid(ndim - 1)*ngrid(ndim - 2)
159  DO is = 1, ndw
160  itt = 1
161  DO it = 1, is - 1
162  itt = itt*ngrid(ndim - it)
163  END DO
164  ind = ind + (pos(ndim - is + 1) - 1)*itt
165  END DO
166  IF (ind > np) cpabort("something wrong in indexing ..")
167  END IF
168  pnt = point_no_pbc(pos, ngrid, ndim)
169  xx = x0 + dp_grid*(pos - 1)
170  dimval = product(ngrid(1:ndim - ndw))
171 
172  IF (.NOT. l_fes_int) THEN
173  IF (PRESENT(array)) THEN
174  array(ind) = minval(-fes(pnt:pnt + dimval - 1))
175  ELSE
176  WRITE (unit_nr, '(10f20.10)') (xx(id), id=ndim, ndim - ndw + 1, -1), minval(-fes(pnt:pnt + dimval - 1))
177  END IF
178  ELSE
179  sum_fes = 0.0_dp
180  dvol = 1.0_dp
181  dvol = product(dp_grid(1:ndim - ndw))
182  DO is = pnt, pnt + dimval - 1
183  sum_fes = sum_fes + fes(is)*dvol
184  END DO
185  IF (PRESENT(array)) THEN
186  array(ind) = -sum_fes
187  ELSE
188  WRITE (unit_nr, '(10f20.10)') (xx(id), id=ndim, ndim - ndw + 1, -1), -sum_fes
189  END IF
190  END IF
191  END IF
192  END DO
193  DEALLOCATE (xx)
194 
195  END SUBROUTINE fes_write
196 
197 ! **************************************************************************************************
198 !> \brief Writes the FES on the file when stride is requested
199 !> \param idim ...
200 !> \param fes ...
201 !> \param pos ...
202 !> \param ndim ...
203 !> \param ngrid ...
204 !> \param dp_grid ...
205 !> \param ndw ...
206 !> \param l_fes_int ...
207 !> \param unit_nr ...
208 !> \par History
209 !> 03.2006 created [tlaino]
210 !> teodoro.laino .at. gmail.com
211 !> \author Teodoro Laino
212 ! **************************************************************************************************
213  RECURSIVE SUBROUTINE fes_only_write(idim, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
214  INTEGER, INTENT(IN) :: idim
215  REAL(kind=dp), DIMENSION(:), POINTER :: fes
216  INTEGER, DIMENSION(:), POINTER :: pos
217  INTEGER, INTENT(IN) :: ndim
218  INTEGER, DIMENSION(:), POINTER :: ngrid
219  REAL(kind=dp), DIMENSION(:), POINTER :: dp_grid
220  INTEGER, INTENT(IN) :: ndw
221  LOGICAL, INTENT(IN) :: l_fes_int
222  INTEGER :: unit_nr
223 
224  INTEGER :: dimval, i, is, pnt
225  REAL(kind=dp) :: dvol, sum_fes
226 
227  DO i = 1, ngrid(idim)
228  pos(idim) = i
229  IF (idim /= ndim - ndw + 1) THEN
230  CALL fes_only_write(idim - 1, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
231  ELSE
232  pnt = point_no_pbc(pos, ngrid, ndim)
233  dimval = product(ngrid(1:ndim - ndw))
234  IF (l_fes_int) THEN
235  WRITE (unit_nr, '(1f12.5)') minval(-fes(pnt:pnt + dimval - 1))
236  ELSE
237  sum_fes = 0.0_dp
238  dvol = product(dp_grid(1:ndim - ndw))
239  DO is = pnt, pnt + dimval - 1
240  sum_fes = sum_fes + fes(is)*dvol
241  END DO
242  WRITE (unit_nr, '(1f12.5)') - sum_fes
243  END IF
244  END IF
245  END DO
246 
247  END SUBROUTINE fes_only_write
248 
249 ! **************************************************************************************************
250 !> \brief Finds minima of the FES
251 !> \param fes ...
252 !> \param ndim ...
253 !> \param iperd ...
254 !> \param ngrid ...
255 !> \param dp_grid ...
256 !> \param x0 ...
257 !> \param ndw ...
258 !> \par History
259 !> 06.2009 created [tlaino]
260 !> teodoro.laino .at. gmail.com
261 !> \author Teodoro Laino
262 ! **************************************************************************************************
263  SUBROUTINE fes_min(fes, ndim, iperd, ngrid, dp_grid, x0, ndw)
264  REAL(kind=dp), DIMENSION(:), POINTER :: fes
265  INTEGER, INTENT(IN) :: ndim
266  INTEGER, DIMENSION(:), POINTER :: iperd, ngrid
267  REAL(kind=dp), DIMENSION(:), POINTER :: dp_grid, x0
268  INTEGER, INTENT(IN) :: ndw
269 
270  INTEGER :: i, id, iter, j, k, max_ntrust, nacc, &
271  ntrials, pnt
272  INTEGER, ALLOCATABLE, DIMENSION(:, :) :: history
273  INTEGER, DIMENSION(:), POINTER :: pos, pos0
274  INTEGER, DIMENSION(ndim) :: dpos, ntrust
275  LOGICAL :: do_save
276  REAL(kind=dp) :: fes_now, fes_old, norm_dx, resto
277  REAL(kind=dp), DIMENSION(:), POINTER :: dx, rnd, xx
278 
279  IF (ndw /= ndim) cpabort("Not implemented for projected FES!")
280 
281  ntrust = ngrid/10
282  ntrials = product(ngrid)
283  WRITE (*, '(A,10I6)', advance="no") "FES| Trust hyper-radius ", ntrust
284  WRITE (*, '(A,10F12.6)') " which is equivalent to: ", ntrust*dp_grid
285 
286  ALLOCATE (xx(ndim), dx(ndim), pos0(ndim), rnd(ndim), pos(ndim))
287  ALLOCATE (history(ndim, ntrials))
288  history = 0
289  nacc = 0
290  trials: DO j = 1, ntrials
291  ! Loop over all points
292  pnt = j
293  DO k = ndim, 2, -1
294  pos0(k) = pnt/product(ngrid(1:k - 1))
295  resto = mod(pnt, product(ngrid(1:k - 1)))
296  IF (resto /= 0) THEN
297  pnt = pnt - pos0(k)*product(ngrid(1:k - 1))
298  pos0(k) = pos0(k) + 1
299  ELSE
300  pnt = product(ngrid(1:k - 1))
301  END IF
302  END DO
303  pos0(1) = pnt
304 
305  ! Loop over the frame points unless it is periodic
306  DO k = 1, ndim
307  IF ((iperd(k) == 0) .AND. (pos0(k) < ntrust(k))) cycle trials
308  IF ((iperd(k) == 0) .AND. (pos0(k) > ngrid(k) - ntrust(k))) cycle trials
309  END DO
310 
311  ! Evaluate position and derivative
312  pos = pos0
313  xx = x0 + dp_grid*(pos - 1)
314  dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid)
315 
316  ! Integrate till derivative is small enough..
317  pnt = point_no_pbc(pos, ngrid, ndim)
318  fes_now = -fes(pnt)
319  fes_old = huge(0.0_dp)
320 
321  i = 1
322  DO WHILE ((i <= 100) .OR. (fes_now < fes_old))
323  fes_old = fes_now
324  !WRITE(10+j,'(10f20.10)')(xx(id),id=ndim,1,-1),-fes(pnt)
325 
326  norm_dx = sqrt(dot_product(dx, dx))
327  IF (norm_dx == 0.0_dp) EXIT ! It is in a really flat region
328  xx = xx - min(0.1_dp, norm_dx)*dx/norm_dx
329  ! Re-evaluating pos
330  pos = ceiling((xx - x0)/dp_grid) + 1
331  CALL pbc(pos, iperd, ngrid, ndim)
332 
333  ! Incremental pos
334  dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid)
335  pnt = point_no_pbc(pos, ngrid, ndim)
336  fes_now = -fes(pnt)
337  i = i + 1
338  END DO
339  iter = i
340 
341  ! Compare with the available minima and if they are the same skip
342  ! saving this position..
343  do_save = fes(pnt) >= 1.0e-3_dp
344  IF (do_save) THEN
345  DO i = 1, nacc
346  dpos = pos - history(:, i)
347  norm_dx = dot_product(dpos, dpos)
348  max_ntrust = maxval(ntrust)
349  ! (SQRT(REAL(norm_dx, KIND=dp)) <= MAXVAL(ntrust)) ...
350  IF ((norm_dx <= real(max_ntrust*max_ntrust, kind=dp)) .OR. (fes(pnt) < 1.0e-3_dp)) THEN
351  do_save = .false.
352  EXIT
353  END IF
354  END DO
355  END IF
356  IF (do_save) THEN
357  pnt = point_no_pbc(pos, ngrid, ndim)
358  xx = x0 + dp_grid*(pos - 1)
359  WRITE (*, '(A,5F12.6)', advance="NO") "FES| Minimum found (", (xx(id), id=ndim, ndim - ndw + 1, -1)
360  WRITE (*, '(A,F12.6,A,I6)') " ). FES value = ", -fes(pnt), " Hartree. Number of Iter: ", iter
361  nacc = nacc + 1
362  history(:, nacc) = pos
363  END IF
364  END DO trials
365  WRITE (*, '(A,I6,A)') "FES| Number of Minimum found: ", nacc, "."
366 
367  DEALLOCATE (xx, dx, pos0, rnd, pos)
368  DEALLOCATE (history)
369 
370  END SUBROUTINE fes_min
371 
372 ! **************************************************************************************************
373 !> \brief Finds path between two points (a) and (b)
374 !> \param fes ...
375 !> \param ndim ...
376 !> \param ngrid ...
377 !> \param dp_grid ...
378 !> \param iperd ...
379 !> \param x0 ...
380 !> \param ndw ...
381 !> \param mep_input_data ...
382 !> \param l_int ...
383 !> \par History
384 !> 12.2010 created [tlaino]
385 !> teodoro.laino .at. gmail.com
386 !> \author Teodoro Laino
387 ! **************************************************************************************************
388  SUBROUTINE fes_path(fes, ndim, ngrid, dp_grid, iperd, x0, ndw, mep_input_data, l_int)
389  REAL(kind=dp), DIMENSION(:), POINTER :: fes
390  INTEGER, INTENT(IN) :: ndim
391  INTEGER, DIMENSION(:), POINTER :: ngrid
392  REAL(kind=dp), DIMENSION(:), POINTER :: dp_grid
393  INTEGER, DIMENSION(:), POINTER :: iperd
394  REAL(kind=dp), DIMENSION(:), POINTER :: x0
395  INTEGER, INTENT(IN) :: ndw
396  TYPE(mep_input_data_type), INTENT(IN) :: mep_input_data
397  LOGICAL :: l_int
398 
399  INTEGER :: i, id, irep, iter, nf, nreplica, ns, &
400  pnt, unit_nr
401  INTEGER, DIMENSION(:), POINTER :: ipos
402  LOGICAL :: converged
403  REAL(kind=dp) :: avg1, avg2, diff, ene, norm_dx, xx0, yy0
404  REAL(kind=dp), DIMENSION(:), POINTER :: davg1, davg2, dxx, dyy, fes_rep, tang, &
405  xx, yy
406  REAL(kind=dp), DIMENSION(:, :), POINTER :: dx, pos, pos_old
407 
408  IF (ndw /= ndim) cpabort("Not implemented for projected FES!")
409  nreplica = mep_input_data%nreplica
410  ALLOCATE (xx(ndim), dx(ndim, nreplica), pos_old(ndim, nreplica), pos(ndim, nreplica), &
411  ipos(ndim), fes_rep(nreplica), dxx(ndim), dyy(ndim), yy(ndim), davg1(ndim), &
412  tang(ndim), davg2(ndim))
413 
414  IF (l_int) THEN
415  id = 0
416  DO i = ndim, ndim - ndw + 1, -1
417  id = id + 1
418  pos(i, 1) = mep_input_data%minima(id, 1)
419  pos(i, nreplica) = mep_input_data%minima(id, 2)
420  END DO
421 
422  ! Interpolate nreplica-2 points
423  xx = (pos(:, nreplica) - pos(:, 1))/real(nreplica - 1, kind=dp)
424  DO irep = 2, nreplica - 1
425  pos(:, irep) = pos(:, 1) + xx(:)*real(irep - 1, kind=dp)
426  END DO
427 
428  ELSE
429  pos = mep_input_data%minima
430  END IF
431 
432  ! Compute value and derivative in all replicas
433  DO irep = 1, nreplica
434  ipos = floor((pos(:, irep) - x0)/dp_grid) + 1
435  pnt = point_no_pbc(ipos, ngrid, ndim)
436  dx(:, irep) = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid)
437  fes_rep(irep) = -fes(pnt)
438  END DO
439 
440  ! Implement a simple elastic band method (Hamiltonian): definitely not the best
441  ! method, but for such a simple task it should be more than enough
442  converged = .false.
443  pos_old = pos
444  iter = 0
445  DO WHILE ((.NOT. converged) .AND. (iter <= mep_input_data%max_iter))
446  iter = iter + 1
447  avg1 = 0.0_dp
448  ! compute average length (distance 1)
449  DO irep = 2, nreplica
450  xx = pos(:, irep) - pos(:, irep - 1)
451  avg1 = avg1 + sqrt(dot_product(xx, xx))
452  END DO
453  avg1 = avg1/real(nreplica - 1, kind=dp)
454 
455  avg2 = 0.0_dp
456  ! compute average length (distance 2)
457  DO irep = 3, nreplica
458  xx = pos(:, irep) - pos(:, irep - 2)
459  avg2 = avg2 + sqrt(dot_product(xx, xx))
460  END DO
461  avg2 = avg2/real(nreplica - 2, kind=dp)
462 
463  ! compute energy and derivatives
464  dx = 0.0_dp
465  ene = 0.0_dp
466  ns = 1
467  nf = nreplica
468  DO irep = 1, nreplica
469  ! compute energy and map point replica irep
470  ipos = floor((pos(:, irep) - x0)/dp_grid) + 1
471  pnt = point_no_pbc(ipos, ngrid, ndim)
472  fes_rep(irep) = -fes(pnt)
473  IF ((irep == 1) .OR. (irep == nreplica)) cycle
474 
475  ! -------------------------------------------------------------
476  ! compute non-linear elastic terms : including only 2-d springs
477  ! -------------------------------------------------------------
478  davg2 = 0.0_dp
479  IF (irep < nf - 1) THEN
480  xx = pos(:, irep) - pos(:, irep + 2)
481  xx0 = sqrt(dot_product(xx, xx))
482  dxx = 1.0_dp/xx0*xx
483  ene = ene + 0.25_dp*mep_input_data%kb*(xx0 - avg2)**2
484  davg2 = davg2 + dxx
485  END IF
486 
487  IF (irep > ns + 1) THEN
488  xx = pos(:, irep) - pos(:, irep - 2)
489  yy0 = sqrt(dot_product(xx, xx))
490  dyy = 1.0_dp/yy0*xx
491  davg2 = davg2 + dyy
492  END IF
493  davg2 = davg2/real(nreplica - 2, kind=dp)
494 
495  IF (irep < nf - 1) THEN
496  dx(:, irep) = dx(:, irep) + 0.5_dp*mep_input_data%kb*(xx0 - avg2)*(dxx - davg2)
497  END IF
498  IF (irep > ns + 1) THEN
499  dx(:, irep) = dx(:, irep) + 0.5_dp*mep_input_data%kb*(yy0 - avg2)*(dyy - davg2)
500  END IF
501 
502  ! -------------------------------------------------------------
503  ! Evaluation of the elastic term
504  ! -------------------------------------------------------------
505  xx = pos(:, irep) - pos(:, irep + 1)
506  yy0 = sqrt(dot_product(xx, xx))
507  dyy = 1.0_dp/yy0*xx
508 
509  xx = pos(:, irep) - pos(:, irep - 1)
510  xx0 = sqrt(dot_product(xx, xx))
511  dxx = 1.0_dp/xx0*xx
512  davg1 = (dxx + dyy)/real(nreplica - 1, kind=dp)
513 
514  ene = ene + 0.5_dp*mep_input_data%kb*(xx0 - avg1)**2
515  dx(:, irep) = dx(:, irep) + mep_input_data%kb*(xx0 - avg1)*(dxx - davg1) + &
516  mep_input_data%kb*(yy0 - avg1)*(dyy - davg1)
517 
518  ! Evaluate the tangent
519  xx = pos(:, irep + 1) - pos(:, irep)
520  xx = xx/sqrt(dot_product(xx, xx))
521  yy = pos(:, irep) - pos(:, irep - 1)
522  yy = yy/sqrt(dot_product(yy, yy))
523  tang = xx + yy
524  tang = tang/sqrt(dot_product(tang, tang))
525 
526  xx = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid)
527  dx(:, irep) = dot_product(dx(:, irep), tang)*tang + &
528  xx - dot_product(xx, tang)*tang
529  END DO
530  dx(:, 1) = 0.0_dp
531  dx(:, nreplica) = 0.0_dp
532 
533  ! propagate the band with a SD step
534  diff = 0.0_dp
535  DO irep = 1, nreplica
536  ene = ene + fes_rep(irep)
537  IF ((irep == 1) .OR. (irep == nreplica)) cycle
538 
539  norm_dx = sqrt(dot_product(dx(:, irep), dx(:, irep)))
540  IF (norm_dx /= 0.0_dp) THEN
541  pos(:, irep) = pos(:, irep) - min(0.1_dp, norm_dx)*dx(:, irep)/norm_dx
542  END IF
543  xx = pos(:, irep) - pos_old(:, irep)
544  diff = diff + dot_product(xx, xx)
545  END DO
546  ! SQRT(diff) <= 0.001_dp
547  IF (diff <= 1.0e-6_dp) THEN
548  converged = .true.
549  END IF
550  pos_old = pos
551  WRITE (*, *) "Iteration nr.", iter, sqrt(diff)
552  END DO
553 
554  WRITE (*, *) "MEP saved on <mep.data> file."
555  CALL open_file(unit_number=unit_nr, file_name="mep.data", file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
556  DO irep = 1, nreplica
557  ! compute energy and derivative for each single point of the replica
558  ipos = floor((pos(:, irep) - x0)/dp_grid) + 1
559  pnt = point_no_pbc(ipos, ngrid, ndim)
560  fes_rep(irep) = -fes(pnt)
561  WRITE (unit_nr, *) irep, pos(:, nreplica - irep + 1), fes_rep(nreplica - irep + 1)
562  END DO
563  CALL close_file(unit_nr)
564 
565  DEALLOCATE (xx, dx, pos, fes_rep, ipos, pos_old, yy, dyy, dxx, davg1, tang, davg2)
566  END SUBROUTINE fes_path
567 
568 ! **************************************************************************************************
569 !> \brief Dump FES with a GAUSSIAN cube format - Useful for multidimensional FES
570 !> \param idim ...
571 !> \param fes ...
572 !> \param pos ...
573 !> \param ndim ...
574 !> \param ngrid ...
575 !> \param dp_grid ...
576 !> \param x0 ...
577 !> \param ndw ...
578 !> \param l_fes_int ...
579 !> \param file ...
580 !> \par History
581 !> 12.2013 created [tlaino]
582 !> teodoro.laino .at. gmail.com
583 !> \author Teodoro Laino
584 ! **************************************************************************************************
585  SUBROUTINE fes_cube_write(idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, file)
586  INTEGER, INTENT(IN) :: idim
587  REAL(kind=dp), DIMENSION(:), POINTER :: fes
588  INTEGER, DIMENSION(:), POINTER :: pos
589  INTEGER, INTENT(IN) :: ndim
590  INTEGER, DIMENSION(:), POINTER :: ngrid
591  REAL(kind=dp), DIMENSION(:), POINTER :: dp_grid, x0
592  INTEGER, INTENT(IN) :: ndw
593  LOGICAL, INTENT(IN) :: l_fes_int
594  CHARACTER(LEN=80) :: file
595 
596  CHARACTER(LEN=120) :: line
597  CHARACTER(LEN=5) :: label, labelp
598  INTEGER :: i, id(3), ii, iix, iiy, iiz, ix, iy, iz, &
599  natoms, np
600  INTEGER, DIMENSION(:), POINTER :: izat
601  REAL(kind=dp) :: cell(3, 3), delta(3), dr(3), residual, &
602  rt(3)
603  REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: rho, rhot
604  REAL(kind=dp), DIMENSION(:, :), POINTER :: xat
605 
606  CALL init_periodic_table()
607  IF (ndw .GT. 3) THEN
608  WRITE (*, *)
609  WRITE (*, *) 'ERROR: GAUSSIAN format can only handle FES on 3 CV !'
610  cpabort("")
611  END IF
612 
613  OPEN (10, file=file, status='old')
614  CALL get_val_res(unit=10, section="&SUBSYS", subsection="&CELL")
615  READ (10, *) label, cell(1, 1), cell(2, 1), cell(3, 1)
616  READ (10, *) label, cell(1, 2), cell(2, 2), cell(3, 2)
617  READ (10, *) label, cell(1, 3), cell(2, 3), cell(3, 3)
618  rt(1) = -(cell(1, 1)/2._dp)
619  rt(2) = -(cell(2, 2)/2._dp)
620  rt(3) = -(cell(3, 3)/2._dp)
621 
622  WRITE (*, *) 'Dumping GAUSSIAN CUBE format'
623  WRITE (*, *) 'Cell vectors'
624  WRITE (*, *)
625 
626  residual = 0.0d0
627  DO ix = 1, 3
628  DO iy = ix + 1, 3
629  residual = residual + cell(ix, iy)**2
630  END DO
631  END DO
632 
633  IF (residual .GT. 1.0d-6) THEN
634  WRITE (*, *)
635  WRITE (*, *) 'ERROR: this program can only handle orthogonal cells'
636  WRITE (*, *) ' with vectors pointing in the X, Y and Z directions'
637  cpabort("")
638  END IF
639 
640  WRITE (*, *)
641  WRITE (*, *) 'Cube grid mesh: ', ngrid(1), 'x', ngrid(2), 'x', ngrid(3)
642  WRITE (*, *) 'Origin in:', rt
643  WRITE (*, *)
644 
645  DO ix = 1, 3
646  dr(ix) = cell(ix, ix)/real(ngrid(ix), kind=dp)
647  END DO
648 
649  np = ngrid(1)*ngrid(2)*ngrid(3)
650  ALLOCATE (rho(np), rhot(np))
651  CALL fes_write(123, idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, rho)
652  WRITE (*, *) 'Internal FES transfer completed!'
653 
654  ! translate cell
655  DO ix = 1, 3
656  delta(ix) = rt(ix)/dr(ix)
657  id(ix) = int(delta(ix))
658  delta(ix) = rt(ix) - id(ix)*dr(ix)
659  END DO
660 
661  DO iz = 1, ngrid(3)
662  DO iy = 1, ngrid(2)
663  DO ix = 1, ngrid(1)
664  iix = ix + id(1)
665  iiy = iy + id(2)
666  iiz = iz + id(3)
667  IF (iix .LT. 1) iix = iix + ngrid(1)
668  IF (iiy .LT. 1) iiy = iiy + ngrid(2)
669  IF (iiz .LT. 1) iiz = iiz + ngrid(3)
670  IF (iix .GT. ngrid(1)) iix = iix - ngrid(1)
671  IF (iiy .GT. ngrid(2)) iiy = iiy - ngrid(2)
672  IF (iiz .GT. ngrid(3)) iiz = iiz - ngrid(3)
673 
674  IF (iix .LT. 1) cpabort("ix < 0")
675  IF (iiy .LT. 1) cpabort("iy < 0")
676  IF (iiz .LT. 1) cpabort("iz < 0")
677  IF (iix .GT. ngrid(1)) cpabort("ix > cell")
678  IF (iiy .GT. ngrid(2)) cpabort("iy > cell")
679  IF (iiz .GT. ngrid(3)) cpabort("iz > cell")
680  i = ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)*ngrid(2)
681  ii = iix + (iiy - 1)*ngrid(1) + (iiz - 1)*ngrid(1)*ngrid(2)
682  rhot(ii) = rho(i)
683  END DO
684  END DO
685  END DO
686 
687  rewind(10)
688  CALL get_val_res(unit=10, section="&SUBSYS", subsection="&COORD")
689  natoms = 0
690  ALLOCATE (xat(1000, 3))
691  ALLOCATE (izat(1000))
692  DO WHILE (.true.)
693  READ (10, '(A)') line
694  IF (index(line, '&END') /= 0) EXIT
695  natoms = natoms + 1
696  READ (line, *) label, (xat(natoms, ix), ix=1, 3)
697  IF (natoms == SIZE(xat, 1)) THEN
698  CALL reallocate(xat, 1, SIZE(xat, 1)*2, 1, 3)
699  CALL reallocate(izat, 1, SIZE(izat)*2)
700  END IF
701  CALL uppercase(label)
702  DO i = 1, nelem
703  labelp = ptable(i)%symbol
704  CALL uppercase(labelp)
705  IF (trim(label) == trim(labelp)) EXIT
706  END DO
707  IF (i == nelem + 1) THEN
708  WRITE (*, *) trim(label), "In line: ", line
709  cpabort("Element not recognized!")
710  END IF
711  izat(natoms) = i
712  END DO
713  CALL reallocate(xat, 1, natoms, 1, 3)
714  CALL reallocate(izat, 1, natoms)
715 
716  DO i = 1, natoms
717  DO ix = 1, 3
718  xat(i, ix) = xat(i, ix) + rt(ix) - delta(ix)
719  IF (xat(i, ix) .LT. rt(ix)) xat(i, ix) = xat(i, ix) + cell(ix, ix)
720  IF (xat(i, ix) .GT. -rt(ix)) xat(i, ix) = xat(i, ix) - cell(ix, ix)
721  END DO
722  END DO
723 
724  WRITE (123, *) "FES on CUBE"
725  WRITE (123, *) "created by fes in CP2K"
726  WRITE (123, '(i5,3f12.6)') natoms, rt(1:3)*bohr
727 
728  DO ix = 1, 3
729  ii = ngrid(ix)
730  WRITE (123, '(i5,4f12.6)') ii, (cell(ix, iy)/ii*bohr, iy=1, 3)
731  END DO
732 
733  DO i = 1, natoms
734  WRITE (123, '(i5,4f12.6)') izat(i), 0.0, (xat(i, ix)*bohr, ix=1, 3)
735  END DO
736 
737  DO ix = 1, ngrid(1)
738  DO iy = 1, ngrid(2)
739 
740  WRITE (123, '(6e13.5)') (rhot(ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)&
741  &*ngrid(2)), iz=1, ngrid(3))
742  END DO
743  END DO
744  DEALLOCATE (xat, rho, rhot)
745 
746  END SUBROUTINE fes_cube_write
747 
748 END MODULE graph_methods
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:308
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:119
graph_methods
Program to Map on grid the hills spawned during a metadynamics run.
Definition: graph_methods.F:20
graph_methods::fes_min
subroutine, public fes_min(fes, ndim, iperd, ngrid, dp_grid, x0, ndw)
Finds minima of the FES.
Definition: graph_methods.F:264
graph_methods::fes_compute_low
recursive subroutine, public fes_compute_low(idim, nn, fes, gauss, ind, ind0, nfes, ndim, ngauss, ngrid, iperd)
Efficiently map the gaussians on the grid.
Definition: graph_methods.F:70
graph_methods::fes_cube_write
subroutine, public fes_cube_write(idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, file)
Dump FES with a GAUSSIAN cube format - Useful for multidimensional FES.
Definition: graph_methods.F:586
graph_methods::fes_write
recursive subroutine, public fes_write(unit_nr, idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, array)
Writes the FES on the file.
Definition: graph_methods.F:129
graph_methods::fes_only_write
recursive subroutine, public fes_only_write(idim, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
Writes the FES on the file when stride is requested.
Definition: graph_methods.F:214
graph_methods::fes_path
subroutine, public fes_path(fes, ndim, ngrid, dp_grid, iperd, x0, ndw, mep_input_data, l_int)
Finds path between two points (a) and (b)
Definition: graph_methods.F:389
graph_utils
Module containing utils for mapping FESs.
Definition: graph_utils.F:18
graph_utils::derivative
real(kind=dp) function, dimension(ndim), public derivative(fes, pos0, iperd, ndim, ngrid, dp_grid)
computes the derivative of the FES w.r.t CVs
Definition: graph_utils.F:56
graph_utils::get_val_res
subroutine, public get_val_res(unit, section, keyword, subsection, i_val, r_val)
Parser informations from the cp2k input/restart.
Definition: graph_utils.F:193
graph_utils::point_no_pbc
integer function, public point_no_pbc(pos, ngrid, ndim)
Computes the pointer to the 1D array given the n-dimensional position non-PBC version.
Definition: graph_utils.F:166
graph_utils::pbc
subroutine, public pbc(pos, iperd, ngrid, ndim)
Computes the pointer to the 1D array given the n-dimensional position PBC version.
Definition: graph_utils.F:136
graph_utils::point_pbc
integer function, public point_pbc(pos, iperd, ngrid, ndim)
Computes the pointer to the 1D array given the n-dimensional position PBC version.
Definition: graph_utils.F:96
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
memory_utilities
Utility routines for the memory handling.
Definition: memory_utilities.F:14
periodic_table
Periodic Table related data definitions.
Definition: periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition: periodic_table.F:65
periodic_table::init_periodic_table
subroutine, public init_periodic_table()
Initialization of Periodic Table related data.
Definition: periodic_table.F:137
periodic_table::nelem
integer, parameter, public nelem
Definition: periodic_table.F:63
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147
string_utilities
Utilities for string manipulations.
Definition: string_utilities.F:16
string_utilities::uppercase
elemental subroutine, public uppercase(string)
Convert all lower case characters in a string to upper case.
Definition: string_utilities.F:3361