(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
20MODULE 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
49CONTAINS
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
748END MODULE graph_methods
memory_utilities::reallocate
Definition memory_utilities.F:27
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
graph_utils::mep_input_data_type
Definition graph_utils.F:25