(git:374b731)
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-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!>
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') .NE. 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') .NE. 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') .NE. 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') .NE. 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') .NE. 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') .NE. 0) THEN
164 l_fes_int = .true.
165 END IF
166
167 IF (index(wq_char, '-orac') .NE. 0) THEN
168 l_orac = .true.
169 END IF
170
171 IF (index(wq_char, '-cp2k') .NE. 0) THEN
172 l_cp2k = .true.
173 END IF
174
175 IF (index(wq_char, '-cpmd') .NE. 0) THEN
176 l_cpmd = .true.
177 END IF
178
179 IF (index(wq_char, '-find-minima') .NE. 0) THEN
180 l_fmin = .true.
181 END IF
182
183 IF (index(wq_char, '-find-path') .NE. 0) THEN
184 l_pmin = .true.
185 END IF
186
187 IF (index(wq_char, '-mathlab') .NE. 0) THEN
188 l_math = .true.
189 END IF
190
191 IF (index(wq_char, '-cube') .NE. 0) THEN
192 l_cube = .true.
193 END IF
194
195 IF (index(wq_char, '-opt-path') .NE. 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) &
204 cpabort("Error! You've to specify either ORAC, CP2K or CPMD!!")
205
206 ! For CPMD move filename to colvar_mtd
207 IF (l_cpmd) THEN
208 file = "colvar_mtd"
209 END IF
210
211 ! Initializing random numbers
212 CALL random_seed()
213 CALL random_number(dum)
214
215 ! Basic Allocation
216 ndw = ndim
217 ALLOCATE (ngrid(ndim), stat=stat)
218 IF (stat /= 0) cpabort("Allocation Error")
219 ALLOCATE (dp_grid(ndim), stat=stat)
220 IF (stat /= 0) cpabort("Allocation Error")
221 ALLOCATE (idw(ndw), stat=stat)
222 IF (stat /= 0) cpabort("Allocation Error")
223 ALLOCATE (iperd(ndim), stat=stat)
224 IF (stat /= 0) cpabort("Allocation Error")
225 ALLOCATE (iprd(nprd), stat=stat)
226 IF (stat /= 0) cpabort("Allocation Error")
227 DO i = 1, ndim
228 idw(i) = i
229 iperd(i) = 0
230 END DO
231
232 DO i = 1, command_argument_count()
233 CALL get_command_argument(i, wq_char, status=istat)
234 cpassert(istat == 0)
235
236 IF (index(wq_char, '-ndw') .NE. 0) THEN
237 DEALLOCATE (idw)
238
239 ndw = 0
240 ndw_loop: DO ix = i + 1, command_argument_count()
241 CALL get_command_argument(ix, wq_char, status=istat)
242 cpassert(istat == 0)
243 IF (index(wq_char, '-') .EQ. 0) THEN
244 ndw = ndw + 1
245 ELSE
246 EXIT ndw_loop
247 END IF
248 END DO ndw_loop
249
250 ALLOCATE (idw(ndw), stat=stat)
251 IF (stat /= 0) cpabort("Allocation Error")
252
253 DO id = 1, ndw
254 CALL get_command_argument(i + id, wq_char, status=istat)
255 cpassert(istat == 0)
256 READ (wq_char, *) idw(id)
257 END DO
258 END IF
259
260 IF (index(wq_char, '-periodic') .NE. 0) THEN
261 nprd = 0
262 nprd_loop: DO ix = i + 1, command_argument_count()
263 CALL get_command_argument(ix, wq_char, status=istat)
264 cpassert(istat == 0)
265 IF (index(wq_char, '-') .EQ. 0) THEN
266 nprd = nprd + 1
267 ELSE
268 EXIT nprd_loop
269 END IF
270 END DO nprd_loop
271
272 DEALLOCATE (iprd)
273 ALLOCATE (iprd(nprd), stat=stat)
274 IF (stat /= 0) cpabort("Allocation Error")
275
276 DO id = 1, nprd
277 CALL get_command_argument(i + id, wq_char, status=istat)
278 cpassert(istat == 0)
279 READ (wq_char, *) iprd(id)
280 END DO
281 END IF
282
283 IF (index(wq_char, '-ngrid') .NE. 0) THEN
284 DO ix = 1, ndim
285 CALL get_command_argument(i + ix, wq_char, status=istat)
286 cpassert(istat == 0)
287 READ (wq_char, *) ngrid(ix)
288 l_grid = .true.
289 END DO
290 END IF
291
292 IF (index(wq_char, '-dp') .NE. 0) THEN
293 l_dp = .true.
294 l_grid = .false.
295 DO ix = 1, ndim
296 CALL get_command_argument(i + ix, wq_char, status=istat)
297 cpassert(istat == 0)
298 READ (wq_char, *) dp_grid(ix)
299 END DO
300 END IF
301
302 IF (index(wq_char, '-fix') .NE. 0) THEN
303 fix = .true.
304 DO id = 1, ndw
305 CALL get_command_argument(i + 2*(id - 1) + 1, wq_char, status=istat)
306 cpassert(istat == 0)
307 READ (wq_char, *) x0w(id)
308 CALL get_command_argument(i + 2*(id - 1) + 2, wq_char, status=istat)
309 cpassert(istat == 0)
310 READ (wq_char, *) xfw(id)
311 END DO
312 END IF
313 END DO
314
315 IF (l_pmin) THEN
316 ALLOCATE (mep_input_data%minima(ndw, 2))
317 mep_input_data%minima = huge(0.0_dp)
318 mep_input_data%max_iter = 10000
319 mep_input_data%kb = 0.1_dp
320 mep_input_data%nreplica = 8
321 ! Read for starting point (a) and (b)
322 DO i = 1, command_argument_count()
323 CALL get_command_argument(i, wq_char, status=istat)
324 cpassert(istat == 0)
325
326 IF (index(wq_char, '-point-a') .NE. 0) THEN
327 DO id = 1, ndw
328 CALL get_command_argument(i + id, wq_char, status=istat)
329 cpassert(istat == 0)
330 READ (wq_char, *) mep_input_data%minima(id, 1)
331 END DO
332 END IF
333
334 IF (index(wq_char, '-point-b') .NE. 0) THEN
335 DO id = 1, ndw
336 CALL get_command_argument(i + id, wq_char, status=istat)
337 cpassert(istat == 0)
338 READ (wq_char, *) mep_input_data%minima(id, 2)
339 END DO
340 END IF
341
342 IF (index(wq_char, '-mep-iter') .NE. 0) THEN
343 CALL get_command_argument(i + 1, wq_char, status=istat)
344 cpassert(istat == 0)
345 READ (wq_char, *) mep_input_data%max_iter
346 END IF
347
348 IF (index(wq_char, '-mep-kb') .NE. 0) THEN
349 CALL get_command_argument(i + 1, wq_char, status=istat)
350 cpassert(istat == 0)
351 READ (wq_char, *) mep_input_data%kb
352 END IF
353
354 IF (index(wq_char, '-mep-nreplica') .NE. 0) THEN
355 CALL get_command_argument(i + 1, wq_char, status=istat)
356 cpassert(istat == 0)
357 READ (wq_char, *) mep_input_data%nreplica
358 END IF
359
360 END DO
361 IF (any(mep_input_data%minima == huge(0.0_dp))) &
362 cpabort("-find-path requires the specification of -point-a and -point-b !")
363 ELSE
364 ALLOCATE (mep_input_data%minima(0, 0))
365 END IF
366
367! Read parameters for Path_optimization
368 IF (l_popt) THEN
369 mep_input_data%nreplica = 0
370 mep_input_data%max_iter = 10000
371 mep_input_data%kb = 0.1_dp
372
373 DO i = 1, command_argument_count()
374 CALL get_command_argument(i, wq_char, status=istat)
375 cpassert(istat == 0)
376
377 IF (index(wq_char, '-mep-kb') .NE. 0) THEN
378 CALL get_command_argument(i + 1, wq_char, status=istat)
379 cpassert(istat == 0)
380 READ (wq_char, *) mep_input_data%kb
381 END IF
382
383 IF (index(wq_char, '-mep-iter') .NE. 0) THEN
384 CALL get_command_argument(i + 1, wq_char, status=istat)
385 cpassert(istat == 0)
386 READ (wq_char, *) mep_input_data%max_iter
387 END IF
388
389 IF (index(wq_char, '-mep-nreplica') .NE. 0) THEN
390 CALL get_command_argument(i + 1, wq_char, status=istat)
391 cpassert(istat == 0)
392 READ (wq_char, *) mep_input_data%nreplica
393 END IF
394 END DO
395
396 ALLOCATE (mep_input_data%minima(ndw, mep_input_data%nreplica))
397
398 CALL open_file(unit_number=unit_nr, file_name=path_file, file_status="OLD")
399 DO id = 1, mep_input_data%nreplica
400 READ (unit_nr, *) j, mep_input_data%minima(:, id)
401 END DO
402 CALL close_file(unit_nr)
403
404 DO id = 1, mep_input_data%nreplica
405 WRITE (*, *) mep_input_data%minima(:, id)
406 END DO
407 END IF
408
409 ! Defines the order of the collectiv var.: first the "wanted" ones, then the others
410 ALLOCATE (i_map(ndim), stat=stat)
411 IF (stat /= 0) cpabort("Allocation Error")
412 i_map = 0
413
414 DO id = 1, ndw
415 i_map(idw(id)) = id
416 END DO
417 ix = ndw
418 DO id = 1, ndim
419 IF (i_map(id) == 0) THEN
420 ix = ix + 1
421 i_map(id) = ix
422 END IF
423 END DO
424
425 ! Revert the order so we can perform averages (when projecting FES) more
426 ! efficiently
427 i_map = ndim - i_map + 1
428
429 ! Tag the periodic COLVAR according the new internal order
430 DO id = 1, nprd
431 iperd(i_map(iprd(id))) = 1
432 END DO
433
434 ! Grid size
435 IF (l_grid) THEN
436 ALLOCATE (tmp(ndim), stat=stat)
437 IF (stat /= 0) cpabort("Allocation Error")
438 tmp = ngrid
439 DO i = 1, ndim
440 ngrid(i_map(i)) = tmp(i)
441 END DO
442 DEALLOCATE (tmp)
443 ELSE
444 ngrid = 100
445 END IF
446
447 WRITE (iw, '(/,70("*"))')
448 WRITE (iw, '("FES|",T7,A,/)') "Parsing file: <"//trim(file)//">"
449
450 CALL open_file(unit_number=unit_nr, file_name=file, file_status="OLD")
451 IF (l_cp2k) THEN
452 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="NHILLS_START_VAL", i_val=nt)
453 ! These sections may not necessarily be present.. if not the values will be HUGE and negative..
454 ! If sc>0 but p and q are not defined, it fails miserably
455 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="HILL_TAIL_CUTOFF", r_val=sc)
456 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="P_EXPONENT", i_val=p)
457 CALL get_val_res(unit=unit_nr, section="&METADYN", keyword="Q_EXPONENT", i_val=q)
458 ELSE IF (l_orac .OR. l_cpmd) THEN
459 nt = 0
460 DO WHILE (.true.)
461 READ (unit_nr, *, END=100, ERR=100) dum
462 nt = nt + 1
463 END DO
464100 rewind(unit_nr)
465 END IF
466
467 ALLOCATE (x0(ndim), stat=stat)
468 IF (stat /= 0) cpabort("Allocation Error")
469 ALLOCATE (xf(ndim), stat=stat)
470 IF (stat /= 0) cpabort("Allocation Error")
471 ALLOCATE (ss0(ndim, nt), stat=stat)
472 IF (stat /= 0) cpabort("Allocation Error")
473 ALLOCATE (delta_s(ndim, nt), stat=stat)
474 IF (stat /= 0) cpabort("Allocation Error")
475 ALLOCATE (ww(nt), stat=stat)
476 IF (stat /= 0) cpabort("Allocation Error")
477 ALLOCATE (ind(ndim), stat=stat)
478 IF (stat /= 0) cpabort("Allocation Error")
479 ALLOCATE (inds(ndim), stat=stat)
480 IF (stat /= 0) cpabort("Allocation Error")
481 ALLOCATE (nn(ndim, nt), stat=stat)
482 IF (stat /= 0) cpabort("Allocation Error")
483 ALLOCATE (nn_max(ndim), stat=stat)
484 IF (stat /= 0) cpabort("Allocation Error")
485 ALLOCATE (dp_cut(ndim), stat=stat)
486 IF (stat /= 0) cpabort("Allocation Error")
487
488 IF (l_cp2k) THEN
489 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_POS")
490 DO i = 1, nt
491 READ (unit_nr, '(A120)') a
492 DO WHILE (a(len_trim(a):len_trim(a)) == '\\')
493 READ (unit_nr, '(A120)') b
494 a = a(1:len_trim(a) - 1)//b(1:len_trim(b))
495 END DO
496 READ (a, *) (ss0(i_map(id), i), id=1, ndim)
497 END DO
498 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_SCALE")
499 DO i = 1, nt
500 READ (unit_nr, '(A120)') a
501 DO WHILE (a(len_trim(a):len_trim(a)) == '\\')
502 READ (unit_nr, '(A120)') b
503 a = a(1:len_trim(a) - 1)//b(1:len_trim(b))
504 END DO
505 READ (a, *) (delta_s(i_map(id), i), id=1, ndim)
506 END DO
507 CALL get_val_res(unit=unit_nr, section="&METADYN", subsection="&SPAWNED_HILLS_HEIGHT")
508 DO i = 1, nt
509 READ (unit_nr, *) ww(i)
510 END DO
511 ELSE IF (l_orac) THEN
512 DO i = 1, nt
513 READ (unit_nr, *) dum, (ss0(i_map(id), i), id=1, ndim), (delta_s(i_map(id), i), id=1, ndim), ww(i)
514 END DO
515 ELSE IF (l_cpmd) THEN
516 CALL open_file(unit_number=unit_nr2, file_name="parvar_mtd", file_status="OLD")
517 DO i = 1, nt
518 READ (unit_nr, *) dum, (ss0(i_map(id), i), id=1, ndim), (delta_s(id, i), id=1, ndim)
519 READ (unit_nr2, *) dum, dum, delta_s_glob, ww(i)
520 delta_s(1:ndim, i) = delta_s_glob*delta_s(1:ndim, i)
521 END DO
522 CALL close_file(unit_nr2)
523 END IF
524 CALL close_file(unit_nr)
525
526 ! ORAC conversion factor
527 IF (l_orac) ww = ww*10000._dp/4.187_dp
528
529 ! Setting up the limit of definitions for the several colvars
530 DO id = 1, ndim
531 x0(id) = huge(1.0_dp)
532 xf(id) = -huge(1.0_dp)
533 END DO
534 IF (fix) THEN
535 DO it = 1, nt
536 DO id = 1, ndim - ndw
537 x0(id) = min(x0(id), ss0(id, it) - 3.*delta_s(id, it))
538 xf(id) = max(xf(id), ss0(id, it) + 3.*delta_s(id, it))
539 END DO
540 END DO
541 it = 0
542 DO id = ndim, ndim - ndw + 1, -1
543 it = it + 1
544 x0(id) = x0w(it)
545 xf(id) = xfw(it)
546 END DO
547 ELSE
548 DO it = 1, nt
549 DO id = ndim, 1, -1
550 IF (iperd(id) == 1) THEN
551 x0(id) = -pi
552 xf(id) = pi
553 ELSE
554 x0(id) = min(x0(id), ss0(id, it) - 3.*delta_s(id, it))
555 xf(id) = max(xf(id), ss0(id, it) + 3.*delta_s(id, it))
556 END IF
557 END DO
558 END DO
559 END IF
560
561 IF (l_dp) THEN
562 ALLOCATE (tmpr(ndim))
563 tmpr = dp_grid
564 DO i = 1, ndim
565 dp_grid(i_map(i)) = tmpr(i)
566 END DO
567 DEALLOCATE (tmpr)
568 ngrid = int((xf - x0)/dp_grid) + 1
569 ELSE
570 dp_grid = (xf - x0)/real(ngrid - 1, kind=dp)
571 END IF
572
573 WRITE (iw, '(70("*"))')
574 WRITE (iw, '("FES|",T7,A,/)') "Parameters for FES:"
575 WRITE (iw, '("FES|",T7,A15,5x,i7)') "NDIM ::", ndim
576 WRITE (iw, '("FES|",T7,A15,5x,i7)') "NWD ::", ndw
577 WRITE (iw, '("FES|",T7,A15,5x,i7)') "HILLS ::", nt
578 it = 0
579 DO i = ndim, 1, -1
580 it = it + 1
581 per_label = ""
582 active_label = "(NO MAPPED)"
583 IF (iperd(i) /= 0) per_label = "(PERIODIC)"
584 IF (it <= ndw) active_label = "( MAPPED)"
585 j = minloc((i_map - i)**2, 1)
586 WRITE (iw, '("FES|",T7,"COLVAR # ",i3," ::",5x,"(",f7.3," ,",f7.3,")",T48,A,T60,A)') &
587 j, x0(i), xf(i), trim(per_label), trim(active_label)
588 END DO
589 WRITE (iw, '("FES|",T7,a15,5x,7i7)') "NGRID ::", (ngrid(id), id=ndim, ndim - ndw + 1, -1)
590 WRITE (iw, '("FES|",T7,a15,5x,5f7.3)') "DX ::", (dp_grid(id), id=ndim, ndim - ndw + 1, -1)
591 WRITE (iw, '("FES|",T7,a15,5x,g10.5)') "CUTOFF ::", eps_cut
592 WRITE (iw, '(70("*"),/)')
593
594 nn_max = 0
595 DO i = 1, nt
596 dp_cut = sqrt(log(abs(ww(i))/eps_cut))*2.0_dp*delta_s(:, i)
597 nn(:, i) = int(dp_cut/dp_grid)
598 ww(i) = ww(i)**(1.0_dp/real(ndim, kind=dp))
599 END DO
600
601 nn_max = maxval(nn, dim=2)
602 ngauss = maxval(nn_max)*2 + 1
603 nfes = product(ngrid)
604
605 ALLOCATE (gauss(-maxval(nn_max):maxval(nn_max), ndim))
606 ALLOCATE (fes(nfes))
607 fes = 0.0_dp
608
609 nh = 1
610 nf = min(nh + nt_p - 1, nt)
611
612 IF (lstride) THEN
613 nwr = nt_p
614 ELSE
615 nwr = int(nt/10) + 1
616 END IF
617
618 ncount = 0
619 WRITE (iw, '(/,"FES|",T7,A)') "Computing Free Energy Surface"
620
621 stride: DO WHILE (nh <= nt)
622 hills: DO it = nh, nf
623 ind = int((ss0(:, it) - x0)/dp_grid) + 1
624 gauss = 0.0_dp
625
626 DO i = 1, ndim
627 coor = ind(i) - nn(i, it) - 1
628 ss = x0(i) + coor*dp_grid(i) - dp_grid(i)
629 DO ip = -nn(i, it), nn(i, it)
630 coor = coor + 1
631 ss = ss + dp_grid(i)
632 IF (iperd(i) == 0) THEN
633 IF (coor .GT. ngrid(i)) cycle
634 IF (coor .LT. 1) cycle
635 END IF
636 diff = ss - ss0(i, it)
637 dp2 = (diff/delta_s(i, it))**2
638 gauss(ip, i) = ww(it)*exp(-0.5_dp*dp2)
639 IF (sc > 0.0_dp .AND. p > 0.0_dp .AND. q > 0.0_dp .AND. q > p) THEN
640 gauss(ip, i) = gauss(ip, i)*(1 - (diff/sc*delta_s(i, it))**p)/(1 - (diff/sc*delta_s(i, it))**q)
641 END IF
642 END DO
643 END DO
644 inds = ind
645 CALL fes_compute_low(ndim, nn(:, it), fes, gauss, ind, inds, nfes, ndim, ngauss, ngrid, iperd)
646
647 IF (.NOT. lstride .AND. mod(it, nwr) == 0) THEN
648 WRITE (iw, '("FES|",T7,a,i4,a2)') "Mapping Gaussians ::", int(10*anint(10.*it/nt)), " %"
649 ELSEIF (.NOT. lstride .AND. it == nt) THEN
650 WRITE (iw, '("FES|",T7,a,i4,a2)') "Mapping Gaussians ::", int(10*anint(10.*it/nt)), " %"
651 END IF
652 END DO hills
653
654 IF (lstride) THEN
655 ncount = ncount + 1
656 WRITE (iw, '("FES|",T7,a13,i5," |-| Gaussians from ",i6," to",i6)') "Done frame ::", ncount, nh, nf
657 IF (l_fes_int) THEN
658 out3 = out2//"."
659 ELSE
660 out3 = out1//"."
661 END IF
662
663 IF (ncount < 10) THEN
664 WRITE (out3_stride, '(A,i1)') trim(out3), ncount
665 ELSEIF (ncount < 100) THEN
666 WRITE (out3_stride, '(A,i2)') trim(out3), ncount
667 ELSE
668 WRITE (out3_stride, '(A,i3)') trim(out3), ncount
669 END IF
670 CALL open_file(unit_number=unit_nr, file_name=out3_stride, file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
671 ind = 1
672 CALL fes_only_write(ndim, fes, ind, ndim, ngrid, dp_grid, ndw, l_fes_int, unit_nr)
673 CALL close_file(unit_nr)
674 END IF
675
676 nh = nh + nt_p
677 nf = min(nh + nt_p - 1, nt)
678 END DO stride
679 DEALLOCATE (gauss)
680
681 IF (l_fes_int) THEN
682 out3 = out2
683 ELSE
684 out3 = out1
685 END IF
686
687 WRITE (iw, '("FES|",T7,A)') "Dumping FES structure in file: < "//trim(out3)//" >"
688 CALL open_file(unit_number=unit_nr, file_name=out3, file_action="WRITE", file_status="UNKNOWN", file_form="FORMATTED")
689 IF (l_cube) THEN
690 ind = 1
691 CALL fes_cube_write(ndim, fes, ind, ndim, ngrid, dp_grid, x0, ndw, l_fes_int, file)
692 ELSE
693 ix = 0
694 IF (l_math) WRITE (unit_nr, '(10g12.5)') (ngrid(id), id=ndim, ndim - ndw + 1, -1), ix
695 ind = 1
696 CALL fes_write(unit_nr, ndim, fes, ind, ndim, ngrid, dp_grid, x0, ndw, l_fes_int)
697 END IF
698 CALL close_file(unit_nr)
699
700 ! If requested find minima
701 IF (l_fmin) CALL fes_min(fes, ndim, iperd, ngrid, dp_grid, x0, ndw)
702
703 ! If requested find or opt path
704 IF ((l_pmin) .AND. (l_popt)) cpabort("USE EITHER -find-path OR -opt-path")
705 IF (l_pmin) l_int = .true.
706 IF (l_popt) l_int = .false.
707
708 IF ((l_pmin) .OR. (l_popt)) CALL fes_path(fes, ndim, ngrid, dp_grid, iperd, x0, ndw, mep_input_data, l_int)
709
710 ! Free memory
711 DEALLOCATE (ngrid)
712 DEALLOCATE (dp_grid)
713 DEALLOCATE (idw)
714 DEALLOCATE (iperd)
715 DEALLOCATE (x0)
716 DEALLOCATE (xf)
717 DEALLOCATE (ss0)
718 DEALLOCATE (delta_s)
719 DEALLOCATE (ww)
720 DEALLOCATE (ind)
721 DEALLOCATE (inds)
722 DEALLOCATE (nn)
723 DEALLOCATE (nn_max)
724 DEALLOCATE (dp_cut)
725 DEALLOCATE (i_map)
726 DEALLOCATE (fes)
727 DEALLOCATE (iprd)
728 DEALLOCATE (mep_input_data%minima)
729
730 ! Terminate FES
731 WRITE (iw, '(/,A,/)') "FES| NORMAL FES TERMINATION."
732
733END 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: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::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