Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 0 : RECURSIVE SUBROUTINE fes_compute_low(idim, nn, fes, gauss, ind, ind0, nfes, ndim, &
69 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: ll, pos
81 : REAL(KIND=dp) :: prod
82 :
83 0 : ALLOCATE (pos(ndim), ll(ndim))
84 0 : pos = ind
85 0 : k = nn(idim)
86 :
87 0 : DO i = -k, k
88 0 : pos(idim) = ind(idim) + i
89 0 : IF (iperd(idim) == 0) THEN
90 0 : IF (pos(idim) > ngrid(idim)) CYCLE
91 0 : IF (pos(idim) < 1) CYCLE
92 : END IF
93 0 : IF (idim /= 1) THEN
94 0 : CALL fes_compute_low(idim - 1, nn, fes, gauss, pos, ind0, nfes, ndim, ngauss, ngrid, iperd)
95 : ELSE
96 0 : pnt = point_pbc(pos, iperd, ngrid, ndim)
97 0 : prod = 1.0_dp
98 0 : DO j = 1, ndim
99 0 : ll(j) = pos(j) - ind0(j)
100 0 : prod = prod*gauss(ll(j), j)
101 : END DO
102 0 : fes(pnt) = fes(pnt) + prod
103 : END IF
104 : END DO
105 0 : DEALLOCATE (pos, ll)
106 :
107 0 : 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 0 : RECURSIVE SUBROUTINE fes_write(unit_nr, idim, fes, pos, ndim, ngrid, &
128 0 : 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 0 : ALLOCATE (xx(ndim))
144 0 : xx = x0
145 0 : DO i = 1, ngrid(idim)
146 0 : pos(idim) = i
147 0 : IF (idim /= ndim - ndw + 1) THEN
148 0 : IF (PRESENT(array)) THEN
149 : CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, &
150 0 : x0, ndw, l_fes_int, array)
151 : ELSE
152 : CALL fes_write(unit_nr, idim - 1, fes, pos, ndim, ngrid, dp_grid, &
153 0 : x0, ndw, l_fes_int)
154 : END IF
155 : ELSE
156 0 : IF (PRESENT(array)) THEN
157 0 : ind = 1
158 0 : np = ngrid(ndim)*ngrid(ndim - 1)*ngrid(ndim - 2)
159 0 : DO is = 1, ndw
160 : itt = 1
161 0 : DO it = 1, is - 1
162 0 : itt = itt*ngrid(ndim - it)
163 : END DO
164 0 : ind = ind + (pos(ndim - is + 1) - 1)*itt
165 : END DO
166 0 : IF (ind > np) CPABORT("something wrong in indexing ..")
167 : END IF
168 0 : pnt = point_no_pbc(pos, ngrid, ndim)
169 0 : xx = x0 + dp_grid*(pos - 1)
170 0 : dimval = PRODUCT(ngrid(1:ndim - ndw))
171 :
172 0 : IF (.NOT. l_fes_int) THEN
173 0 : IF (PRESENT(array)) THEN
174 0 : array(ind) = MINVAL(-fes(pnt:pnt + dimval - 1))
175 : ELSE
176 0 : 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 0 : sum_fes = 0.0_dp
180 0 : dvol = 1.0_dp
181 0 : dvol = PRODUCT(dp_grid(1:ndim - ndw))
182 0 : DO is = pnt, pnt + dimval - 1
183 0 : sum_fes = sum_fes + fes(is)*dvol
184 : END DO
185 0 : IF (PRESENT(array)) THEN
186 0 : array(ind) = -sum_fes
187 : ELSE
188 0 : 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 0 : DEALLOCATE (xx)
194 :
195 0 : 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 0 : 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 0 : DO i = 1, ngrid(idim)
228 0 : pos(idim) = i
229 0 : IF (idim /= ndim - ndw + 1) THEN
230 0 : CALL fes_only_write(idim - 1, fes, pos, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
231 : ELSE
232 0 : pnt = point_no_pbc(pos, ngrid, ndim)
233 0 : dimval = PRODUCT(ngrid(1:ndim - ndw))
234 0 : IF (l_fes_int) THEN
235 0 : WRITE (unit_nr, '(1f12.5)') MINVAL(-fes(pnt:pnt + dimval - 1))
236 : ELSE
237 0 : sum_fes = 0.0_dp
238 0 : dvol = PRODUCT(dp_grid(1:ndim - ndw))
239 0 : DO is = pnt, pnt + dimval - 1
240 0 : sum_fes = sum_fes + fes(is)*dvol
241 : END DO
242 0 : WRITE (unit_nr, '(1f12.5)') - sum_fes
243 : END IF
244 : END IF
245 : END DO
246 :
247 0 : 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 0 : 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 0 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: history
273 0 : INTEGER, DIMENSION(:), POINTER :: pos, pos0
274 0 : INTEGER, DIMENSION(ndim) :: Dpos, ntrust
275 : LOGICAL :: do_save
276 : REAL(KIND=dp) :: fes_now, fes_old, norm_dx, resto
277 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: dx, rnd, xx
278 :
279 0 : IF (ndw /= ndim) CPABORT("Not implemented for projected FES!")
280 :
281 0 : ntrust = ngrid/10
282 0 : ntrials = PRODUCT(ngrid)
283 0 : WRITE (*, '(A,10I6)', ADVANCE="no") "FES| Trust hyper-radius ", ntrust
284 0 : WRITE (*, '(A,10F12.6)') " which is equivalent to: ", ntrust*dp_grid
285 :
286 0 : ALLOCATE (xx(ndim), dx(ndim), pos0(ndim), rnd(ndim), pos(ndim))
287 0 : ALLOCATE (history(ndim, ntrials))
288 0 : history = 0
289 0 : nacc = 0
290 0 : Trials: DO j = 1, ntrials
291 : ! Loop over all points
292 : pnt = j
293 0 : DO k = ndim, 2, -1
294 0 : pos0(k) = pnt/PRODUCT(ngrid(1:k - 1))
295 0 : resto = MOD(pnt, PRODUCT(ngrid(1:k - 1)))
296 0 : IF (resto /= 0) THEN
297 0 : pnt = pnt - pos0(k)*PRODUCT(ngrid(1:k - 1))
298 0 : pos0(k) = pos0(k) + 1
299 : ELSE
300 0 : pnt = PRODUCT(ngrid(1:k - 1))
301 : END IF
302 : END DO
303 0 : pos0(1) = pnt
304 :
305 : ! Loop over the frame points unless it is periodic
306 0 : DO k = 1, ndim
307 0 : IF ((iperd(k) == 0) .AND. (pos0(k) < ntrust(k))) CYCLE Trials
308 0 : IF ((iperd(k) == 0) .AND. (pos0(k) > ngrid(k) - ntrust(k))) CYCLE Trials
309 : END DO
310 :
311 : ! Evaluate position and derivative
312 0 : pos = pos0
313 0 : xx = x0 + dp_grid*(pos - 1)
314 0 : dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid)
315 :
316 : ! Integrate till derivative is small enough..
317 0 : pnt = point_no_pbc(pos, ngrid, ndim)
318 0 : fes_now = -fes(pnt)
319 0 : fes_old = HUGE(0.0_dp)
320 :
321 0 : i = 1
322 0 : DO WHILE ((i <= 100) .OR. (fes_now < fes_old))
323 0 : fes_old = fes_now
324 : !WRITE(10+j,'(10f20.10)')(xx(id),id=ndim,1,-1),-fes(pnt)
325 :
326 0 : norm_dx = SQRT(DOT_PRODUCT(dx, dx))
327 0 : IF (norm_dx == 0.0_dp) EXIT ! It is in a really flat region
328 0 : xx = xx - MIN(0.1_dp, norm_dx)*dx/norm_dx
329 : ! Re-evaluating pos
330 0 : pos = CEILING((xx - x0)/dp_grid) + 1
331 0 : CALL pbc(pos, iperd, ngrid, ndim)
332 :
333 : ! Incremental pos
334 0 : dx = derivative(fes, pos, iperd, ndim, ngrid, dp_grid)
335 0 : pnt = point_no_pbc(pos, ngrid, ndim)
336 0 : fes_now = -fes(pnt)
337 0 : i = i + 1
338 : END DO
339 0 : iter = i
340 :
341 : ! Compare with the available minima and if they are the same skip
342 : ! saving this position..
343 0 : do_save = fes(pnt) >= 1.0E-3_dp
344 0 : IF (do_save) THEN
345 0 : DO i = 1, nacc
346 0 : Dpos = pos - history(:, i)
347 0 : norm_dx = DOT_PRODUCT(Dpos, Dpos)
348 0 : max_ntrust = MAXVAL(ntrust)
349 : ! (SQRT(REAL(norm_dx, KIND=dp)) <= MAXVAL(ntrust)) ...
350 0 : 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 0 : IF (do_save) THEN
357 0 : pnt = point_no_pbc(pos, ngrid, ndim)
358 0 : xx = x0 + dp_grid*(pos - 1)
359 0 : WRITE (*, '(A,5F12.6)', ADVANCE="NO") "FES| Minimum found (", (xx(id), id=ndim, ndim - ndw + 1, -1)
360 0 : WRITE (*, '(A,F12.6,A,I6)') " ). FES value = ", -fes(pnt), " Hartree. Number of Iter: ", iter
361 0 : nacc = nacc + 1
362 0 : history(:, nacc) = pos
363 : END IF
364 : END DO Trials
365 0 : WRITE (*, '(A,I6,A)') "FES| Number of Minimum found: ", nacc, "."
366 :
367 0 : DEALLOCATE (xx, dx, pos0, rnd, pos)
368 0 : DEALLOCATE (history)
369 :
370 0 : 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 0 : 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 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: davg1, davg2, dxx, dyy, fes_rep, tang, &
405 0 : xx, yy
406 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: dx, pos, pos_old
407 :
408 0 : IF (ndw /= ndim) CPABORT("Not implemented for projected FES!")
409 0 : 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 0 : tang(ndim), davg2(ndim))
413 :
414 0 : IF (l_int) THEN
415 0 : id = 0
416 0 : DO i = ndim, ndim - ndw + 1, -1
417 0 : id = id + 1
418 0 : pos(i, 1) = mep_input_data%minima(id, 1)
419 0 : pos(i, nreplica) = mep_input_data%minima(id, 2)
420 : END DO
421 :
422 : ! Interpolate nreplica-2 points
423 0 : xx = (pos(:, nreplica) - pos(:, 1))/REAL(nreplica - 1, KIND=dp)
424 0 : DO irep = 2, nreplica - 1
425 0 : pos(:, irep) = pos(:, 1) + xx(:)*REAL(irep - 1, KIND=dp)
426 : END DO
427 :
428 : ELSE
429 0 : pos = mep_input_data%minima
430 : END IF
431 :
432 : ! Compute value and derivative in all replicas
433 0 : DO irep = 1, nreplica
434 0 : ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1
435 0 : pnt = point_no_pbc(ipos, ngrid, ndim)
436 0 : dx(:, irep) = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid)
437 0 : 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 0 : converged = .FALSE.
443 0 : pos_old = pos
444 0 : iter = 0
445 0 : DO WHILE ((.NOT. converged) .AND. (iter <= mep_input_data%max_iter))
446 0 : iter = iter + 1
447 0 : avg1 = 0.0_dp
448 : ! compute average length (distance 1)
449 0 : DO irep = 2, nreplica
450 0 : xx = pos(:, irep) - pos(:, irep - 1)
451 0 : avg1 = avg1 + SQRT(DOT_PRODUCT(xx, xx))
452 : END DO
453 0 : avg1 = avg1/REAL(nreplica - 1, KIND=dp)
454 :
455 0 : avg2 = 0.0_dp
456 : ! compute average length (distance 2)
457 0 : DO irep = 3, nreplica
458 0 : xx = pos(:, irep) - pos(:, irep - 2)
459 0 : avg2 = avg2 + SQRT(DOT_PRODUCT(xx, xx))
460 : END DO
461 0 : avg2 = avg2/REAL(nreplica - 2, KIND=dp)
462 :
463 : ! compute energy and derivatives
464 0 : dx = 0.0_dp
465 0 : ene = 0.0_dp
466 0 : ns = 1
467 0 : nf = nreplica
468 0 : DO irep = 1, nreplica
469 : ! compute energy and map point replica irep
470 0 : ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1
471 0 : pnt = point_no_pbc(ipos, ngrid, ndim)
472 0 : fes_rep(irep) = -fes(pnt)
473 0 : IF ((irep == 1) .OR. (irep == nreplica)) CYCLE
474 :
475 : ! -------------------------------------------------------------
476 : ! compute non-linear elastic terms : including only 2-d springs
477 : ! -------------------------------------------------------------
478 0 : davg2 = 0.0_dp
479 0 : IF (irep < nf - 1) THEN
480 0 : xx = pos(:, irep) - pos(:, irep + 2)
481 0 : xx0 = SQRT(DOT_PRODUCT(xx, xx))
482 0 : dxx = 1.0_dp/xx0*xx
483 0 : ene = ene + 0.25_dp*mep_input_data%kb*(xx0 - avg2)**2
484 0 : davg2 = davg2 + dxx
485 : END IF
486 :
487 0 : IF (irep > ns + 1) THEN
488 0 : xx = pos(:, irep) - pos(:, irep - 2)
489 0 : yy0 = SQRT(DOT_PRODUCT(xx, xx))
490 0 : dyy = 1.0_dp/yy0*xx
491 0 : davg2 = davg2 + dyy
492 : END IF
493 0 : davg2 = davg2/REAL(nreplica - 2, KIND=dp)
494 :
495 0 : IF (irep < nf - 1) THEN
496 0 : dx(:, irep) = dx(:, irep) + 0.5_dp*mep_input_data%kb*(xx0 - avg2)*(dxx - davg2)
497 : END IF
498 0 : IF (irep > ns + 1) THEN
499 0 : 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 0 : xx = pos(:, irep) - pos(:, irep + 1)
506 0 : yy0 = SQRT(DOT_PRODUCT(xx, xx))
507 0 : dyy = 1.0_dp/yy0*xx
508 :
509 0 : xx = pos(:, irep) - pos(:, irep - 1)
510 0 : xx0 = SQRT(DOT_PRODUCT(xx, xx))
511 0 : dxx = 1.0_dp/xx0*xx
512 0 : davg1 = (dxx + dyy)/REAL(nreplica - 1, KIND=dp)
513 :
514 0 : 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 0 : mep_input_data%kb*(yy0 - avg1)*(dyy - davg1)
517 :
518 : ! Evaluate the tangent
519 0 : xx = pos(:, irep + 1) - pos(:, irep)
520 0 : xx = xx/SQRT(DOT_PRODUCT(xx, xx))
521 0 : yy = pos(:, irep) - pos(:, irep - 1)
522 0 : yy = yy/SQRT(DOT_PRODUCT(yy, yy))
523 0 : tang = xx + yy
524 0 : tang = tang/SQRT(DOT_PRODUCT(tang, tang))
525 :
526 0 : xx = derivative(fes, ipos, iperd, ndim, ngrid, dp_grid)
527 : dx(:, irep) = DOT_PRODUCT(dx(:, irep), tang)*tang + &
528 0 : xx - DOT_PRODUCT(xx, tang)*tang
529 : END DO
530 0 : dx(:, 1) = 0.0_dp
531 0 : dx(:, nreplica) = 0.0_dp
532 :
533 : ! propagate the band with a SD step
534 0 : diff = 0.0_dp
535 0 : DO irep = 1, nreplica
536 0 : ene = ene + fes_rep(irep)
537 0 : IF ((irep == 1) .OR. (irep == nreplica)) CYCLE
538 :
539 0 : norm_dx = SQRT(DOT_PRODUCT(dx(:, irep), dx(:, irep)))
540 0 : IF (norm_dx /= 0.0_dp) THEN
541 0 : pos(:, irep) = pos(:, irep) - MIN(0.1_dp, norm_dx)*dx(:, irep)/norm_dx
542 : END IF
543 0 : xx = pos(:, irep) - pos_old(:, irep)
544 0 : diff = diff + DOT_PRODUCT(xx, xx)
545 : END DO
546 : ! SQRT(diff) <= 0.001_dp
547 0 : IF (diff <= 1.0e-6_dp) THEN
548 0 : converged = .TRUE.
549 : END IF
550 0 : pos_old = pos
551 0 : WRITE (*, *) "Iteration nr.", iter, SQRT(diff)
552 : END DO
553 :
554 0 : WRITE (*, *) "MEP saved on <mep.data> file."
555 0 : CALL open_file(unit_number=unit_nr, file_name="mep.data", file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
556 0 : DO irep = 1, nreplica
557 : ! compute energy and derivative for each single point of the replica
558 0 : ipos = FLOOR((pos(:, irep) - x0)/dp_grid) + 1
559 0 : pnt = point_no_pbc(ipos, ngrid, ndim)
560 0 : fes_rep(irep) = -fes(pnt)
561 0 : WRITE (unit_nr, *) irep, pos(:, nreplica - irep + 1), fes_rep(nreplica - irep + 1)
562 : END DO
563 0 : CALL close_file(unit_nr)
564 :
565 0 : DEALLOCATE (xx, dx, pos, fes_rep, ipos, pos_old, yy, dyy, dxx, davg1, tang, davg2)
566 0 : 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 0 : 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 0 : INTEGER, DIMENSION(:), POINTER :: izat
601 : REAL(KIND=dp) :: cell(3, 3), delta(3), dr(3), residual, &
602 : rt(3)
603 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: rho, rhot
604 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: xat
605 :
606 0 : CALL init_periodic_table()
607 0 : IF (ndw > 3) THEN
608 0 : WRITE (*, *)
609 0 : WRITE (*, *) 'ERROR: GAUSSIAN format can only handle FES on 3 CV !'
610 0 : CPABORT("")
611 : END IF
612 :
613 0 : OPEN (10, file=file, status='old')
614 0 : CALL get_val_res(unit=10, section="&SUBSYS", subsection="&CELL")
615 0 : READ (10, *) label, cell(1, 1), cell(2, 1), cell(3, 1)
616 0 : READ (10, *) label, cell(1, 2), cell(2, 2), cell(3, 2)
617 0 : READ (10, *) label, cell(1, 3), cell(2, 3), cell(3, 3)
618 0 : rt(1) = -(cell(1, 1)/2._dp)
619 0 : rt(2) = -(cell(2, 2)/2._dp)
620 0 : rt(3) = -(cell(3, 3)/2._dp)
621 :
622 0 : WRITE (*, *) 'Dumping GAUSSIAN CUBE format'
623 0 : WRITE (*, *) 'Cell vectors'
624 0 : WRITE (*, *)
625 :
626 0 : residual = 0.0d0
627 0 : DO ix = 1, 3
628 0 : DO iy = ix + 1, 3
629 0 : residual = residual + cell(ix, iy)**2
630 : END DO
631 : END DO
632 :
633 0 : IF (residual > 1.0d-6) THEN
634 0 : WRITE (*, *)
635 0 : WRITE (*, *) 'ERROR: this program can only handle orthogonal cells'
636 0 : WRITE (*, *) ' with vectors pointing in the X, Y and Z directions'
637 0 : CPABORT("")
638 : END IF
639 :
640 0 : WRITE (*, *)
641 0 : WRITE (*, *) 'Cube grid mesh: ', ngrid(1), 'x', ngrid(2), 'x', ngrid(3)
642 0 : WRITE (*, *) 'Origin in:', rt
643 0 : WRITE (*, *)
644 :
645 0 : DO ix = 1, 3
646 0 : dr(ix) = cell(ix, ix)/REAL(ngrid(ix), KIND=dp)
647 : END DO
648 :
649 0 : np = ngrid(1)*ngrid(2)*ngrid(3)
650 0 : ALLOCATE (rho(np), rhot(np))
651 0 : CALL fes_write(123, idim, fes, pos, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, rho)
652 0 : WRITE (*, *) 'Internal FES transfer completed!'
653 :
654 : ! translate cell
655 0 : DO ix = 1, 3
656 0 : delta(ix) = rt(ix)/dr(ix)
657 0 : id(ix) = INT(delta(ix))
658 0 : delta(ix) = rt(ix) - id(ix)*dr(ix)
659 : END DO
660 :
661 0 : DO iz = 1, ngrid(3)
662 0 : DO iy = 1, ngrid(2)
663 0 : DO ix = 1, ngrid(1)
664 0 : iix = ix + id(1)
665 0 : iiy = iy + id(2)
666 0 : iiz = iz + id(3)
667 0 : IF (iix < 1) iix = iix + ngrid(1)
668 0 : IF (iiy < 1) iiy = iiy + ngrid(2)
669 0 : IF (iiz < 1) iiz = iiz + ngrid(3)
670 0 : IF (iix > ngrid(1)) iix = iix - ngrid(1)
671 0 : IF (iiy > ngrid(2)) iiy = iiy - ngrid(2)
672 0 : IF (iiz > ngrid(3)) iiz = iiz - ngrid(3)
673 :
674 0 : IF (iix < 1) CPABORT("ix < 0")
675 0 : IF (iiy < 1) CPABORT("iy < 0")
676 0 : IF (iiz < 1) CPABORT("iz < 0")
677 0 : IF (iix > ngrid(1)) CPABORT("ix > cell")
678 0 : IF (iiy > ngrid(2)) CPABORT("iy > cell")
679 0 : IF (iiz > ngrid(3)) CPABORT("iz > cell")
680 0 : i = ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)*ngrid(2)
681 0 : ii = iix + (iiy - 1)*ngrid(1) + (iiz - 1)*ngrid(1)*ngrid(2)
682 0 : rhot(ii) = rho(i)
683 : END DO
684 : END DO
685 : END DO
686 :
687 0 : REWIND (10)
688 0 : CALL get_val_res(unit=10, section="&SUBSYS", subsection="&COORD")
689 0 : natoms = 0
690 0 : ALLOCATE (xat(1000, 3))
691 0 : ALLOCATE (izat(1000))
692 0 : DO WHILE (.TRUE.)
693 0 : READ (10, '(A)') line
694 0 : IF (INDEX(line, '&END') /= 0) EXIT
695 0 : natoms = natoms + 1
696 0 : READ (line, *) label, (xat(natoms, ix), ix=1, 3)
697 0 : IF (natoms == SIZE(xat, 1)) THEN
698 0 : CALL reallocate(xat, 1, SIZE(xat, 1)*2, 1, 3)
699 0 : CALL reallocate(izat, 1, SIZE(izat)*2)
700 : END IF
701 0 : CALL uppercase(label)
702 0 : DO i = 1, nelem
703 0 : labelp = ptable(i)%symbol
704 0 : CALL uppercase(labelp)
705 0 : IF (TRIM(label) == TRIM(labelp)) EXIT
706 : END DO
707 0 : IF (i == nelem + 1) THEN
708 0 : WRITE (*, *) TRIM(label), "In line: ", line
709 0 : CPABORT("Element not recognized!")
710 : END IF
711 0 : izat(natoms) = i
712 : END DO
713 0 : CALL reallocate(xat, 1, natoms, 1, 3)
714 0 : CALL reallocate(izat, 1, natoms)
715 :
716 0 : DO i = 1, natoms
717 0 : DO ix = 1, 3
718 0 : xat(i, ix) = xat(i, ix) + rt(ix) - delta(ix)
719 0 : IF (xat(i, ix) < rt(ix)) xat(i, ix) = xat(i, ix) + cell(ix, ix)
720 0 : IF (xat(i, ix) > -rt(ix)) xat(i, ix) = xat(i, ix) - cell(ix, ix)
721 : END DO
722 : END DO
723 :
724 0 : WRITE (123, *) "FES on CUBE"
725 0 : WRITE (123, *) "created by fes in CP2K"
726 0 : WRITE (123, '(i5,3f12.6)') natoms, rt(1:3)*bohr
727 :
728 0 : DO ix = 1, 3
729 0 : ii = ngrid(ix)
730 0 : WRITE (123, '(i5,4f12.6)') ii, (cell(ix, iy)/ii*bohr, iy=1, 3)
731 : END DO
732 :
733 0 : DO i = 1, natoms
734 0 : WRITE (123, '(i5,4f12.6)') izat(i), 0.0, (xat(i, ix)*bohr, ix=1, 3)
735 : END DO
736 :
737 0 : DO ix = 1, ngrid(1)
738 0 : DO iy = 1, ngrid(2)
739 :
740 : WRITE (123, '(6e13.5)') (rhot(ix + (iy - 1)*ngrid(1) + (iz - 1)*ngrid(1)&
741 0 : &*ngrid(2)), iz=1, ngrid(3))
742 : END DO
743 : END DO
744 0 : DEALLOCATE (xat, rho, rhot)
745 :
746 0 : END SUBROUTINE fes_cube_write
747 :
748 : END MODULE graph_methods
|
