(git:d3d49ac)
Loading...
Searching...
No Matches
graph.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 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!>
16!> \par Note
17!> Please report any bug to the author
18! **************************************************************************************************
19PROGRAM graph
20
21 USE cp_files, ONLY: close_file,&
22 open_file
23 USE graph_methods, ONLY: fes_compute_low,&
24 fes_cube_write,&
25 fes_min,&
26 fes_only_write,&
27 fes_path,&
28 fes_write
29 USE graph_utils, ONLY: get_val_res,&
30 mep_input_data_type
31 USE kinds, ONLY: default_string_length,&
32 dp
33 USE mathconstants, ONLY: pi
34#include "../base/base_uses.f90"
35
36 IMPLICIT NONE
37
38 CHARACTER(LEN=80) :: file, out1, out2, out3, wq_char, &
39 path_file, out3_stride
40 CHARACTER(LEN=480) :: a, b
41 CHARACTER(LEN=default_string_length) :: active_label, per_label
42 INTEGER :: istat, coor, i, id, ip, &
43 it, iw, ix, j, ncount, ndim, &
44 ndw, nf, nfes, ngauss, nh, &
45 nprd, nt, nt_p, nwr, p, q, &
46 stat, unit_nr, unit_nr2
47 INTEGER, POINTER :: i_map(:), idw(:), ind(:), &
48 inds(:), iperd(:), iprd(:), &
49 ngrid(:), nn(:, :), nn_max(:), &
50 tmp(:)
51 LOGICAL :: fix, l_cp2k, l_cpmd, &
52 l_dp, l_fes_int, l_fmin, &
53 l_grid, l_math, l_orac, &
54 l_pmin, lstride, l_popt, l_int, &
55 l_cube
56 REAL(kind=dp) :: delta_s_glob, diff, dp2, dum, &
57 eps_cut, sc, ss, x0w(3), &
58 xfw(3)
59 REAL(kind=dp), POINTER :: delta_s(:, :), dp_cut(:), dp_grid(:), fes(:), &
60 gauss(:, :), ss0(:, :), tmpr(:), ww(:), x0(:), xf(:)
61 TYPE(mep_input_data_type) :: mep_input_data
62
63 ! Initialize variables
64 nprd = 0
65 ndim = 1
66 ndw = 1
67 nt_p = 9999999
68 eps_cut = 1e-6
69 file = 'HILLS'
70 out1 = 'fes.dat'
71 out2 = 'fes_int.dat'
72 fix = .false.
73 l_fes_int = .false.
74 lstride = .false.
75 l_grid = .false.
76 l_dp = .false.
77 l_orac = .false.
78 l_cp2k = .false.
79 l_cpmd = .false.
80 l_math = .false.
81 l_cube = .false.
82 l_fmin = .false.
83 l_pmin = .false.
84 l_popt = .false.
85 l_int = .false.
86 iw = 6
87
88 IF (command_argument_count() == 0) THEN
89 WRITE (iw, *) 'USAGE:'
90 WRITE (iw, *) 'graf '
91 WRITE (iw, *) '[-ngrid 50 .. ..] (Mesh dimension. Default :: 100)'
92 WRITE (iw, *) '[-dp 0.05 .. ..] (Alternative to -ngrid, allows the specification of the mesh dx)'
93 WRITE (iw, *) '[-ndim 3 ] (Number of collective variables NCV)'
94 WRITE (iw, *) '[-ndw 1 3 .. ] (CVs for the free energy surface)'
95 WRITE (iw, *) '[-periodic 2 3 ..] (CVs with periodic boundary conditions (-pi,pi] )'
96 WRITE (iw, *) '[-stride 10 ] (How often the FES is written)'
97 WRITE (iw, *) '[-fix 1.1 .. ..] (Define the region for the FES)'
98 WRITE (iw, *) ' (If omitted this is automatically calculated)'
99 WRITE (iw, *) '[-cutoff 2. ] (The hills are cutoffed at 2)'
100 WRITE (iw, *) '[-file filename]'
101 WRITE (iw, *) '[-out filename]'
102 WRITE (iw, *) '[-integrated_fes] (When projecting the FES print the integrated value, '
103 WRITE (iw, *) ' rather then the minimum value (minimum value is default))'
104 WRITE (iw, *) '[-orac] (If energies are written in orac intern units)'
105 WRITE (iw, *) '[-cp2k] (Specify if a CP2K restart file is provided)'
106 WRITE (iw, *) '[-cpmd] (Specify if CPMD colvar_mtd and parvar_mtd are provided)'
107 WRITE (iw, *) ' (With CPMD you do not need to specify -file, parvar_mtd and'
108 WRITE (iw, *) ' colvar_mtd are expected to be present in the working directory)'
109 WRITE (iw, *) '[-mathlab] (File storing FES in Mathlab format. Default format Gnuplot)'
110 WRITE (iw, *) '[-cube] (File storing FES in GAUSSIAN CUBE format. Default format Gnuplot)'
111 WRITE (iw, *) '[-find-minima] (Tries to finds all minima in the computed FES)'
112 WRITE (iw, *) '[-find-path] (Finds MEP between all minima (found) in the computed FES)'
113 WRITE (iw, *) '[-point-a] (Specifies point (a) when using -find-path option)'
114 WRITE (iw, *) '[-point-b] (Specifies point (b) when using -find-path option)'
115 WRITE (iw, *) '[-opt-path filename] (Optimize initial MEP of mep-nreplica points in the same format as mep.data)'
116 WRITE (iw, *) '[-mep-kb] (Specifies the value of the force constant for the MEP: default 0.1_dp)'
117 WRITE (iw, *) '[-mep-nreplica] (Specifies the number of replica points used in the MEP: default 8)'
118 WRITE (iw, *) '[-mep-iter] (Specifies the maximum number of iterations used in the MEP: default 10000)'
119 WRITE (iw, *) ''
120 WRITE (iw, *) 'DEFAULT OUTPUT: fes.dat'
121 WRITE (iw, *) ''
122 cpabort("Please provide arguments to run FES!")
123 END IF
124
125 DO i = 1, command_argument_count()
126 CALL get_command_argument(i, wq_char, status=istat)
127 cpassert(istat == 0)
128
129 IF (index(wq_char, '-file') /= 0) THEN
130 CALL get_command_argument(i + 1, wq_char, status=istat)
131 cpassert(istat == 0)
132 READ (wq_char, *) file
133 END IF
134
135 IF (index(wq_char, '-out') /= 0) THEN
136 CALL get_command_argument(i + 1, wq_char, status=istat)
137 cpassert(istat == 0)
138 READ (wq_char, *) out1
139 ! we read only 1 filename. If none is specified we differentiate between fes.dat and fes_int.dat
140 ! otherwise we use the one provided by the user
141 out2 = out1
142 END IF
143
144 IF (index(wq_char, '-ndim') /= 0) THEN
145 CALL get_command_argument(i + 1, wq_char, status=istat)
146 cpassert(istat == 0)
147 READ (wq_char, *) ndim
148 END IF
149
150 IF (index(wq_char, '-stride') /= 0) THEN
151 CALL get_command_argument(i + 1, wq_char, status=istat)
152 cpassert(istat == 0)
153 READ (wq_char, *) nt_p
154 lstride = .true.
155 END IF
156
157 IF (index(wq_char, '-cutoff') /= 0) THEN
158 CALL get_command_argument(i + 1, wq_char, status=istat)
159 cpassert(istat == 0)
160 READ (wq_char, *) eps_cut
161 END IF
162
163 IF (index(wq_char, '-integrated_fes') /= 0) THEN
164 l_fes_int = .true.
165 END IF
166
167 IF (index(wq_char, '-orac') /= 0) THEN
168 l_orac = .true.
169 END IF
170
171 IF (index(wq_char, '-cp2k') /= 0) THEN
172 l_cp2k = .true.
173 END IF
174
175 IF (index(wq_char, '-cpmd') /= 0) THEN
176 l_cpmd = .true.
177 END IF
178
179 IF (index(wq_char, '-find-minima') /= 0) THEN
180 l_fmin = .true.
181 END IF
182
183 IF (index(wq_char, '-find-path') /= 0) THEN
184 l_pmin = .true.
185 END IF
186
187 IF (index(wq_char, '-mathlab') /= 0) THEN
188 l_math = .true.
189 END IF
190
191 IF (index(wq_char, '-cube') /= 0) THEN
192 l_cube = .true.
193 END IF
194
195 IF (index(wq_char, '-opt-path') /= 0) THEN
196 l_popt = .true.
197 CALL get_command_argument(i + 1, wq_char, status=istat)
198 cpassert(istat == 0)
199 READ (wq_char, *) path_file
200 END IF
201
202 END DO
203 IF (count([l_orac, l_cp2k, l_cpmd]) /= 1) THEN
204 cpabort("Error! You've to specify either ORAC, CP2K or CPMD!")
205 END IF
206
207 ! For CPMD move filename to colvar_mtd
208 IF (l_cpmd) THEN
209 file = "colvar_mtd"
210 END IF
211
212 ! Initializing random numbers
213 CALL random_seed()
214 CALL random_number(dum)
215
216 ! Basic Allocation
217 ndw = ndim
218 ALLOCATE (ngrid(ndim), stat=stat)
219 IF (stat /= 0) cpabort("Allocation Error")
220 ALLOCATE (dp_grid(ndim), stat=stat)
221 IF (stat /= 0) cpabort("Allocation Error")
222 ALLOCATE (idw(ndw), stat=stat)
223 IF (stat /= 0) cpabort("Allocation Error")
224 ALLOCATE (iperd(ndim), stat=stat)
225 IF (stat /= 0) cpabort("Allocation Error")
226 ALLOCATE (iprd(nprd), stat=stat)
227 IF (stat /= 0) cpabort("Allocation Error")
228 DO i = 1, ndim
229 idw(i) = i
230 iperd(i) = 0
231 END DO
232
233 DO i = 1, command_argument_count()
234 CALL get_command_argument(i, wq_char, status=istat)
235 cpassert(istat == 0)
236
237 IF (index(wq_char, '-ndw') /= 0) THEN
238 DEALLOCATE (idw)
239
240 ndw = 0
241 ndw_loop: DO ix = i + 1, command_argument_count()
242 CALL get_command_argument(ix, wq_char, status=istat)
243 cpassert(istat == 0)
244 IF (index(wq_char, '-') == 0) THEN
245 ndw = ndw + 1
246 ELSE
247 EXIT ndw_loop
248 END IF
249 END DO ndw_loop
250
251 ALLOCATE (idw(ndw), stat=stat)
252 IF (stat /= 0) cpabort("Allocation Error")
253
254 DO id = 1, ndw
255 CALL get_command_argument(i + id, wq_char, status=istat)
256 cpassert(istat == 0)
257 READ (wq_char, *) idw(id)
258 END DO
259 END IF
260
261 IF (index(wq_char, '-periodic') /= 0) THEN
262 nprd = 0
263 nprd_loop: DO ix = i + 1, command_argument_count()
264 CALL get_command_argument(ix, wq_char, status=istat)
265 cpassert(istat == 0)
266 IF (index(wq_char, '-') == 0) THEN
267 nprd = nprd + 1
268 ELSE
269 EXIT nprd_loop
270 END IF
271 END DO nprd_loop
272
273 DEALLOCATE (iprd)
274 ALLOCATE (iprd(nprd), stat=stat)
275 IF (stat /= 0) cpabort("Allocation Error")
276
277 DO id = 1, nprd
278 CALL get_command_argument(i + id, wq_char, status=istat)
279 cpassert(istat == 0)
280 READ (wq_char, *) iprd(id)
281 END DO
282 END IF
283
284 IF (index(wq_char, '-ngrid') /= 0) THEN
285 DO ix = 1, ndim
286 CALL get_command_argument(i + ix, wq_char, status=istat)
287 cpassert(istat == 0)
288 READ (wq_char, *) ngrid(ix)
289 l_grid = .true.
290 END DO
291 END IF
292
293 IF (index(wq_char, '-dp') /= 0) THEN
294 l_dp = .true.
295 l_grid = .false.
296 DO ix = 1, ndim
297 CALL get_command_argument(i + ix, wq_char, status=istat)
298 cpassert(istat == 0)
299 READ (wq_char, *) dp_grid(ix)
300 END DO
301 END IF
302
303 IF (index(wq_char, '-fix') /= 0) THEN
304 fix = .true.
305 DO id = 1, ndw
306 CALL get_command_argument(i + 2*(id - 1) + 1, wq_char, status=istat)
307 cpassert(istat == 0)
308 READ (wq_char, *) x0w(id)
309 CALL get_command_argument(i + 2*(id - 1) + 2, wq_char, status=istat)
310 cpassert(istat == 0)
311 READ (wq_char, *) xfw(id)
312 END DO
313 END IF
314 END DO
315
316 IF (l_pmin) THEN
317 ALLOCATE (mep_input_data%minima(ndw, 2))
318 mep_input_data%minima = huge(0.0_dp)
319 mep_input_data%max_iter = 10000
320 mep_input_data%kb = 0.1_dp
321 mep_input_data%nreplica = 8
322 ! Read for starting point (a) and (b)
323 DO i = 1, command_argument_count()
324 CALL get_command_argument(i, wq_char, status=istat)
325 cpassert(istat == 0)
326
327 IF (index(wq_char, '-point-a') /= 0) THEN
328 DO id = 1, ndw
329 CALL get_command_argument(i + id, wq_char, status=istat)
330 cpassert(istat == 0)
331 READ (wq_char, *) mep_input_data%minima(id, 1)
332 END DO
333 END IF
334
335 IF (index(wq_char, '-point-b') /= 0) THEN
336 DO id = 1, ndw
337 CALL get_command_argument(i + id, wq_char, status=istat)
338 cpassert(istat == 0)
339 READ (wq_char, *) mep_input_data%minima(id, 2)
340 END DO
341 END IF
342
343 IF (index(wq_char, '-mep-iter') /= 0) THEN
344 CALL get_command_argument(i + 1, wq_char, status=istat)
345 cpassert(istat == 0)
346 READ (wq_char, *) mep_input_data%max_iter
347 END IF
348
349 IF (index(wq_char, '-mep-kb') /= 0) THEN
350 CALL get_command_argument(i + 1, wq_char, status=istat)
351 cpassert(istat == 0)
352 READ (wq_char, *) mep_input_data%kb
353 END IF
354
355 IF (index(wq_char, '-mep-nreplica') /= 0) THEN
356 CALL get_command_argument(i + 1, wq_char, status=istat)
357 cpassert(istat == 0)
358 READ (wq_char, *) mep_input_data%nreplica
359 END IF
360
361 END DO
362 IF (any(mep_input_data%minima == huge(0.0_dp))) THEN
363 CALL cp_abort(__location__, &
364 "-find-path requires the specification of -point-a and -point-b !")
365 END IF
366 ELSE
367 ALLOCATE (mep_input_data%minima(0, 0))
368 END IF
369
370! Read parameters for Path_optimization
371 IF (l_popt) THEN
372 mep_input_data%nreplica = 0
373 mep_input_data%max_iter = 10000
374 mep_input_data%kb = 0.1_dp
375
376 DO i = 1, command_argument_count()
377 CALL get_command_argument(i, wq_char, status=istat)
378 cpassert(istat == 0)
379
380 IF (index(wq_char, '-mep-kb') /= 0) THEN
381 CALL get_command_argument(i + 1, wq_char, status=istat)
382 cpassert(istat == 0)
383 READ (wq_char, *) mep_input_data%kb
384 END IF
385
386 IF (index(wq_char, '-mep-iter') /= 0) THEN
387 CALL get_command_argument(i + 1, wq_char, status=istat)
388 cpassert(istat == 0)
389 READ (wq_char, *) mep_input_data%max_iter
390 END IF
391
392 IF (index(wq_char, '-mep-nreplica') /= 0) THEN
393 CALL get_command_argument(i + 1, wq_char, status=istat)
394 cpassert(istat == 0)
395 READ (wq_char, *) mep_input_data%nreplica
396 END IF
397 END DO
398
399 ALLOCATE (mep_input_data%minima(ndw, mep_input_data%nreplica))
400
401 CALL open_file(unit_number=unit_nr, file_name=path_file, file_status="OLD")
402 DO id = 1, mep_input_data%nreplica
403 READ (unit_nr, *) j, mep_input_data%minima(:, id)
404 END DO
405 CALL close_file(unit_nr)
406
407 DO id = 1, mep_input_data%nreplica
408 WRITE (*, *) mep_input_data%minima(:, id)
409 END DO
410 END IF
411
412 ! Defines the order of the collectiv var.: first the "wanted" ones, then the others
413 ALLOCATE (i_map(ndim), stat=stat)
414 IF (stat /= 0) cpabort("Allocation Error")
415 i_map = 0
416
417 DO id = 1, ndw
418 i_map(idw(id)) = id
419 END DO
420 ix = ndw
421 DO id = 1, ndim
422 IF (i_map(id) == 0) THEN
423 ix = ix + 1
424 i_map(id) = ix
425 END IF
426 END DO
427
428 ! Revert the order so we can perform averages (when projecting FES) more
429 ! efficiently
430 i_map = ndim - i_map + 1
431
432 ! Tag the periodic COLVAR according the new internal order
433 DO id = 1, nprd
434 iperd(i_map(iprd(id))) = 1
435 END DO
436
437 ! Grid size
438 IF (l_grid) THEN
439 ALLOCATE (tmp(ndim), stat=stat)
440 IF (stat /= 0) cpabort("Allocation Error")
441 tmp = ngrid
442 DO i = 1, ndim
443 ngrid(i_map(i)) = tmp(i)
444 END DO
445 DEALLOCATE (tmp)
446 ELSE
447 ngrid = 100
448 END IF
449
450 WRITE (iw, '(/,70("*"))')
451 WRITE (iw, '("FES|",T7,A,/)') "Parsing file: <"//trim(file)//">"
452
453 CALL open_file(unit_number=unit_nr, file_name=file, file_status="OLD")
454 IF (l_cp2k) THEN
455 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="NHILLS_START_VAL", i_val=nt)
456 ! These sections may not necessarily be present.. if not the values will be HUGE and negative..
457 ! If sc>0 but p and q are not defined, it fails miserably
458 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="HILL_TAIL_CUTOFF", r_val=sc)
459 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="P_EXPONENT", i_val=p)
460 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="Q_EXPONENT", i_val=q)
461 ELSE IF (l_orac .OR. l_cpmd) THEN
462 nt = 0
463 DO WHILE (.true.)
464 READ (unit_nr, *, END=100, ERR=100) dum
465 nt = nt + 1
466 END DO
467100 rewind(unit_nr)
468 END IF
469
470 ALLOCATE (x0(ndim), stat=stat)
471 IF (stat /= 0) cpabort("Allocation Error")
472 ALLOCATE (xf(ndim), stat=stat)
473 IF (stat /= 0) cpabort("Allocation Error")
474 ALLOCATE (ss0(ndim, nt), stat=stat)
475 IF (stat /= 0) cpabort("Allocation Error")
476 ALLOCATE (delta_s(ndim, nt), stat=stat)
477 IF (stat /= 0) cpabort("Allocation Error")
478 ALLOCATE (ww(nt), stat=stat)
479 IF (stat /= 0) cpabort("Allocation Error")
480 ALLOCATE (ind(ndim), stat=stat)
481 IF (stat /= 0) cpabort("Allocation Error")
482 ALLOCATE (inds(ndim), stat=stat)
483 IF (stat /= 0) cpabort("Allocation Error")
484 ALLOCATE (nn(ndim, nt), stat=stat)
485 IF (stat /= 0) cpabort("Allocation Error")
486 ALLOCATE (nn_max(ndim), stat=stat)
487 IF (stat /= 0) cpabort("Allocation Error")
488 ALLOCATE (dp_cut(ndim), stat=stat)
489 IF (stat /= 0) cpabort("Allocation Error")
490
491 IF (l_cp2k) THEN
492 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_POS")
493 DO i = 1, nt
494 READ (unit_nr, '(A120)') a
495 DO WHILE (a(len_trim(a):len_trim(a)) == '\\')
496 READ (unit_nr, '(A120)') b
497 a = a(1:len_trim(a) - 1)//b(1:len_trim(b))
498 END DO
499 READ (a, *) (ss0(i_map(id), i), id=1, ndim)
500 END DO
501 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_SCALE")
502 DO i = 1, nt
503 READ (unit_nr, '(A120)') a
504 DO WHILE (a(len_trim(a):len_trim(a)) == '\\')
505 READ (unit_nr, '(A120)') b
506 a = a(1:len_trim(a) - 1)//b(1:len_trim(b))
507 END DO
508 READ (a, *) (delta_s(i_map(id), i), id=1, ndim)
509 END DO
510 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_HEIGHT")
511 DO i = 1, nt
512 READ (unit_nr, *) ww(i)
513 END DO
514 ELSE IF (l_orac) THEN
515 DO i = 1, nt
516 READ (unit_nr, *) dum, (ss0(i_map(id), i), id=1, ndim), (delta_s(i_map(id), i), id=1, ndim), ww(i)
517 END DO
518 ELSE IF (l_cpmd) THEN
519 CALL open_file(unit_number=unit_nr2, file_name="parvar_mtd", file_status="OLD")
520 DO i = 1, nt
521 READ (unit_nr, *) dum, (ss0(i_map(id), i), id=1, ndim), (delta_s(id, i), id=1, ndim)
522 READ (unit_nr2, *) dum, dum, delta_s_glob, ww(i)
523 delta_s(1:ndim, i) = delta_s_glob*delta_s(1:ndim, i)
524 END DO
525 CALL close_file(unit_nr2)
526 END IF
527 CALL close_file(unit_nr)
528
529 ! ORAC conversion factor
530 IF (l_orac) ww = ww*10000._dp/4.187_dp
531
532 ! Setting up the limit of definitions for the several colvars
533 DO id = 1, ndim
534 x0(id) = huge(1.0_dp)
535 xf(id) = -huge(1.0_dp)
536 END DO
537 IF (fix) THEN
538 DO it = 1, nt
539 DO id = 1, ndim - ndw
540 x0(id) = min(x0(id), ss0(id, it) - 3.*delta_s(id, it))
541 xf(id) = max(xf(id), ss0(id, it) + 3.*delta_s(id, it))
542 END DO
543 END DO
544 it = 0
545 DO id = ndim, ndim - ndw + 1, -1
546 it = it + 1
547 x0(id) = x0w(it)
548 xf(id) = xfw(it)
549 END DO
550 ELSE
551 DO it = 1, nt
552 DO id = ndim, 1, -1
553 IF (iperd(id) == 1) THEN
554 x0(id) = -pi
555 xf(id) = pi
556 ELSE
557 x0(id) = min(x0(id), ss0(id, it) - 3.*delta_s(id, it))
558 xf(id) = max(xf(id), ss0(id, it) + 3.*delta_s(id, it))
559 END IF
560 END DO
561 END DO
562 END IF
563
564 IF (l_dp) THEN
565 ALLOCATE (tmpr(ndim))
566 tmpr = dp_grid
567 DO i = 1, ndim
568 dp_grid(i_map(i)) = tmpr(i)
569 END DO
570 DEALLOCATE (tmpr)
571 ngrid = int((xf - x0)/dp_grid) + 1
572 ELSE
573 dp_grid = (xf - x0)/real(ngrid - 1, kind=dp)
574 END IF
575
576 WRITE (iw, '(70("*"))')
577 WRITE (iw, '("FES|",T7,A,/)') "Parameters for FES:"
578 WRITE (iw, '("FES|",T7,A15,5x,i7)') "NDIM ::", ndim
579 WRITE (iw, '("FES|",T7,A15,5x,i7)') "NWD ::", ndw
580 WRITE (iw, '("FES|",T7,A15,5x,i7)') "HILLS ::", nt
581 it = 0
582 DO i = ndim, 1, -1
583 it = it + 1
584 per_label = ""
585 active_label = "(NO MAPPED)"
586 IF (iperd(i) /= 0) per_label = "(PERIODIC)"
587 IF (it <= ndw) active_label = "( MAPPED)"
588 j = minloc((i_map - i)**2, 1)
589 WRITE (iw, '("FES|",T7,"COLVAR # ",i3," ::",5x,"(",f7.3," ,",f7.3,")",T48,A,T60,A)') &
590 j, x0(i), xf(i), trim(per_label), trim(active_label)
591 END DO
592 WRITE (iw, '("FES|",T7,a15,5x,7i7)') "NGRID ::", (ngrid(id), id=ndim, ndim - ndw + 1, -1)
593 WRITE (iw, '("FES|",T7,a15,5x,5f7.3)') "DX ::", (dp_grid(id), id=ndim, ndim - ndw + 1, -1)
594 WRITE (iw, '("FES|",T7,a15,5x,g10.5)') "CUTOFF ::", eps_cut
595 WRITE (iw, '(70("*"),/)')
596
597 nn_max = 0
598 DO i = 1, nt
599 dp_cut = sqrt(log(abs(ww(i))/eps_cut))*2.0_dp*delta_s(:, i)
600 nn(:, i) = int(dp_cut/dp_grid)
601 ww(i) = ww(i)**(1.0_dp/real(ndim, kind=dp))
602 END DO
603
604 nn_max = maxval(nn, dim=2)
605 ngauss = maxval(nn_max)*2 + 1
606 nfes = product(ngrid)
607
608 ALLOCATE (gauss(-maxval(nn_max):maxval(nn_max), ndim))
609 ALLOCATE (fes(nfes))
610 fes = 0.0_dp
611
612 nh = 1
613 nf = min(nh + nt_p - 1, nt)
614
615 IF (lstride) THEN
616 nwr = nt_p
617 ELSE
618 nwr = int(nt/10) + 1
619 END IF
620
621 ncount = 0
622 WRITE (iw, '(/,"FES|",T7,A)') "Computing Free Energy Surface"
623
624 stride: DO WHILE (nh <= nt)
625 hills: DO it = nh, nf
626 ind = int((ss0(:, it) - x0)/dp_grid) + 1
627 gauss = 0.0_dp
628
629 DO i = 1, ndim
630 coor = ind(i) - nn(i, it) - 1
631 ss = x0(i) + coor*dp_grid(i) - dp_grid(i)
632 DO ip = -nn(i, it), nn(i, it)
633 coor = coor + 1
634 ss = ss + dp_grid(i)
635 IF (iperd(i) == 0) THEN
636 IF (coor > ngrid(i)) cycle
637 IF (coor < 1) cycle
638 END IF
639 diff = ss - ss0(i, it)
640 dp2 = (diff/delta_s(i, it))**2
641 gauss(ip, i) = ww(it)*exp(-0.5_dp*dp2)
642 IF (sc > 0.0_dp .AND. p > 0.0_dp .AND. q > 0.0_dp .AND. q > p) THEN
643 gauss(ip, i) = gauss(ip, i)*(1 - (diff/sc*delta_s(i, it))**p)/(1 - (diff/sc*delta_s(i, it))**q)
644 END IF
645 END DO
646 END DO
647 inds = ind
648 CALL fes_compute_low(ndim, nn(:, it), fes, gauss, ind, inds, nfes, ndim, ngauss, ngrid, iperd)
649
650 IF (.NOT. lstride .AND. mod(it, nwr) == 0) THEN
651 WRITE (iw, '("FES|",T7,a,i4,a2)') "Mapping Gaussians ::", int(10*anint(10.*it/nt)), " %"
652 ELSE IF (.NOT. lstride .AND. it == nt) THEN
653 WRITE (iw, '("FES|",T7,a,i4,a2)') "Mapping Gaussians ::", int(10*anint(10.*it/nt)), " %"
654 END IF
655 END DO hills
656
657 IF (lstride) THEN
658 ncount = ncount + 1
659 WRITE (iw, '("FES|",T7,a13,i5," |-| Gaussians from ",i6," to",i6)') "Done frame ::", ncount, nh, nf
660 IF (l_fes_int) THEN
661 out3 = trim(out2)//"."
662 ELSE
663 out3 = trim(out1)//"."
664 END IF
665
666 IF (ncount < 10) THEN
667 WRITE (out3_stride, '(A,i1)') trim(out3), ncount
668 ELSE IF (ncount < 100) THEN
669 WRITE (out3_stride, '(A,i2)') trim(out3), ncount
670 ELSE
671 WRITE (out3_stride, '(A,i3)') trim(out3), ncount
672 END IF
673 CALL open_file(unit_number=unit_nr, file_name=out3_stride, file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
674 ind = 1
675 CALL fes_only_write(ndim, fes, ind, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
676 CALL close_file(unit_nr)
677 END IF
678
679 nh = nh + nt_p
680 nf = min(nh + nt_p - 1, nt)
681 END DO stride
682 DEALLOCATE (gauss)
683
684 IF (l_fes_int) THEN
685 out3 = out2
686 ELSE
687 out3 = out1
688 END IF
689
690 WRITE (iw, '("FES|",T7,A)') "Dumping FES structure in file: < "//trim(out3)//" >"
691 CALL open_file(unit_number=unit_nr, file_name=out3, file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
692 IF (l_cube) THEN
693 ind = 1
694 CALL fes_cube_write(ndim, fes, ind, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, file)
695 ELSE
696 ix = 0
697 IF (l_math) WRITE (unit_nr, '(10g12.5)') (ngrid(id), id=ndim, ndim - ndw + 1, -1), ix
698 ind = 1
699 CALL fes_write(unit_nr, ndim, fes, ind, ndim, ngrid, dp_grid, x0, ndw, l_fes_int)
700 END IF
701 CALL close_file(unit_nr)
702
703 ! If requested find minima
704 IF (l_fmin) CALL fes_min(fes, ndim, iperd, ngrid, dp_grid, x0, ndw)
705
706 ! If requested find or opt path
707 IF ((l_pmin) .AND. (l_popt)) cpabort("USE EITHER -find-path OR -opt-path")
708 IF (l_pmin) l_int = .true.
709 IF (l_popt) l_int = .false.
710
711 IF ((l_pmin) .OR. (l_popt)) CALL fes_path(fes, ndim, ngrid, dp_grid, iperd, x0, ndw, mep_input_data, l_int)
712
713 ! Free memory
714 DEALLOCATE (ngrid)
715 DEALLOCATE (dp_grid)
716 DEALLOCATE (idw)
717 DEALLOCATE (iperd)
718 DEALLOCATE (x0)
719 DEALLOCATE (xf)
720 DEALLOCATE (ss0)
721 DEALLOCATE (delta_s)
722 DEALLOCATE (ww)
723 DEALLOCATE (ind)
724 DEALLOCATE (inds)
725 DEALLOCATE (nn)
726 DEALLOCATE (nn_max)
727 DEALLOCATE (dp_cut)
728 DEALLOCATE (i_map)
729 DEALLOCATE (fes)
730 DEALLOCATE (iprd)
731 DEALLOCATE (mep_input_data%minima)
732
733 ! Terminate FES
734 WRITE (iw, '(/,A,/)') "FES| NORMAL FES TERMINATION."
735
736END PROGRAM graph
graph
program graph
Program to Map on grid the hills spawned during a metadynamics run.
Definition graph.F:19
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
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::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
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
graph_utils::mep_input_data_type
Definition graph_utils.F:25