(git:e7e05ae)
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 ! **************************************************************************************************
19 PROGRAM 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
464 100 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 
733 END 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