1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2020  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \par History
8!>      - taken out of input_cp2k_motion
9!> \author Ole Schuett
10! **************************************************************************************************
11
12MODULE input_cp2k_neb
13   USE bibliography,                    ONLY: Elber1987,&
14                                              Jonsson1998,&
15                                              Jonsson2000_1,&
16                                              Jonsson2000_2,&
17                                              Wales2004
18   USE cp_output_handling,              ONLY: add_last_numeric,&
19                                              cp_print_key_section_create,&
20                                              high_print_level,&
21                                              low_print_level,&
22                                              medium_print_level
23   USE cp_units,                        ONLY: cp_unit_to_cp2k
24   USE input_constants,                 ONLY: &
25        band_diis_opt, band_md_opt, do_b_neb, do_ci_neb, do_d_neb, do_eb, do_it_neb, &
26        do_rep_blocked, do_rep_interleaved, do_sm, pot_neb_fe, pot_neb_full, pot_neb_me
27   USE input_cp2k_thermostats,          ONLY: create_coord_section,&
28                                              create_velocity_section
29   USE input_keyword_types,             ONLY: keyword_create,&
30                                              keyword_release,&
31                                              keyword_type
32   USE input_section_types,             ONLY: section_add_keyword,&
33                                              section_add_subsection,&
34                                              section_create,&
35                                              section_release,&
36                                              section_type
37   USE input_val_types,                 ONLY: real_t
38   USE kinds,                           ONLY: dp
39   USE string_utilities,                ONLY: s2a
40#include "./base/base_uses.f90"
41
42   IMPLICIT NONE
43   PRIVATE
44
45   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
46   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'input_cp2k_neb'
47
48   PUBLIC :: create_band_section
49
50CONTAINS
51
52! **************************************************************************************************
53!> \brief creates the section for a BAND run
54!> \param section will contain the pint section
55!> \author Teodoro Laino 09.2006 [tlaino]
56! **************************************************************************************************
57   SUBROUTINE create_band_section(section)
58      TYPE(section_type), POINTER                        :: section
59
60      CHARACTER(len=*), PARAMETER :: routineN = 'create_band_section', &
61         routineP = moduleN//':'//routineN
62
63      TYPE(keyword_type), POINTER                        :: keyword
64      TYPE(section_type), POINTER                        :: print_key, subsection, subsubsection
65
66      CPASSERT(.NOT. ASSOCIATED(section))
67      CALL section_create(section, __LOCATION__, name="band", &
68                          description="The section that controls a BAND run", &
69                          n_keywords=1, n_subsections=0, repeats=.FALSE., &
70                          citations=(/Elber1987, Jonsson1998, Jonsson2000_1, Jonsson2000_2, Wales2004/))
71      NULLIFY (keyword, print_key, subsection, subsubsection)
72
73      CALL keyword_create(keyword, __LOCATION__, name="NPROC_REP", &
74                          description="Specify the number of processors to be used per replica "// &
75                          "environment (for parallel runs)", &
76                          default_i_val=1)
77      CALL section_add_keyword(section, keyword)
78      CALL keyword_release(keyword)
79
80      CALL keyword_create(keyword, __LOCATION__, name="PROC_DIST_TYPE", &
81                          description="Specify the topology of the mapping of processors into replicas.", &
82                          usage="PROC_DIST_TYPE (INTERLEAVED|BLOCKED)", &
83                          enum_c_vals=s2a("INTERLEAVED", &
84                                          "BLOCKED"), &
85                          enum_desc=s2a("Interleaved distribution", &
86                                        "Blocked distribution"), &
87                          enum_i_vals=(/do_rep_interleaved, do_rep_blocked/), &
88                          default_i_val=do_rep_blocked)
89      CALL section_add_keyword(section, keyword)
90      CALL keyword_release(keyword)
91
92      CALL keyword_create(keyword, __LOCATION__, name="BAND_TYPE", &
93                          description="Specifies the type of BAND calculation", &
94                          usage="BAND_TYPE (B-NEB|IT-NEB|CI-NEB|D-NEB|SM|EB)", &
95                          default_i_val=do_it_neb, &
96                          enum_c_vals=s2a("B-NEB", &
97                                          "IT-NEB", &
98                                          "CI-NEB", &
99                                          "D-NEB", &
100                                          "SM", &
101                                          "EB"), &
102                          enum_desc=s2a("Bisection nudged elastic band", &
103                                        "Improved tangent nudged elastic band", &
104                                        "Climbing image nudged elastic band", &
105                                        "Doubly nudged elastic band", &
106                                        "String Method", &
107                                        "Elastic band (Hamiltonian formulation)"), &
108                          enum_i_vals=(/do_b_neb, do_it_neb, do_ci_neb, do_d_neb, do_sm, do_eb/))
109      CALL section_add_keyword(section, keyword)
110      CALL keyword_release(keyword)
111
112      CALL keyword_create(keyword, __LOCATION__, name="NUMBER_OF_REPLICA", &
113                          description="Specify the number of Replica to use in the BAND", &
114                          default_i_val=10)
115      CALL section_add_keyword(section, keyword)
116      CALL keyword_release(keyword)
117
118      CALL keyword_create(keyword, __LOCATION__, name="USE_COLVARS", &
119                          description="Uses a version of the band scheme projected in a subspace of colvars.", &
120                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
121      CALL section_add_keyword(section, keyword)
122      CALL keyword_release(keyword)
123
124      CALL keyword_create(keyword, __LOCATION__, name="POT_TYPE", &
125                          description="Specifies the type of potential used in the BAND calculation", &
126                          usage="POT_TYPE (FULL|FE|ME)", &
127                          default_i_val=pot_neb_full, &
128                          enum_c_vals=s2a("FULL", &
129                                          "FE", &
130                                          "ME"), &
131                          enum_desc=s2a("Full potential (no projections in a subspace of colvars)", &
132                                        "Free energy (requires a projections in a subspace of colvars)", &
133                                        "Minimum energy (requires a projections in a subspace of colvars)"), &
134                          enum_i_vals=(/pot_neb_full, pot_neb_fe, pot_neb_me/))
135      CALL section_add_keyword(section, keyword)
136      CALL keyword_release(keyword)
137
138      CALL keyword_create(keyword, __LOCATION__, name="ROTATE_FRAMES", &
139                          description="Compute at each BAND step the RMSD and rotate the frames in order"// &
140                          " to minimize it.", &
141                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
142      CALL section_add_keyword(section, keyword)
143      CALL keyword_release(keyword)
144
145      CALL keyword_create(keyword, __LOCATION__, name="ALIGN_FRAMES", &
146                          description="Enables the alignment of the frames at the beginning of a BAND calculation. "// &
147                          "This keyword does not affect the rotation of the replicas during a BAND calculation.", &
148                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
149      CALL section_add_keyword(section, keyword)
150      CALL keyword_release(keyword)
151
152      CALL keyword_create(keyword, __LOCATION__, name="K_SPRING", &
153                          variants=(/"K"/), &
154                          description="Specify the value of the spring constant", &
155                          default_r_val=0.02_dp)
156      CALL section_add_keyword(section, keyword)
157      CALL keyword_release(keyword)
158
159      ! Convergence_control
160      CALL section_create(subsection, __LOCATION__, name="CONVERGENCE_CONTROL", &
161                          description="Setup parameters to control the convergence criteria for BAND", &
162                          repeats=.FALSE.)
163      CALL keyword_create(keyword, __LOCATION__, name="MAX_DR", &
164                          description="Tolerance on the maximum value of the displacement on the BAND.", &
165                          usage="MAX_DR {real}", &
166                          default_r_val=0.0002_dp)
167      CALL section_add_keyword(subsection, keyword)
168      CALL keyword_release(keyword)
169
170      CALL keyword_create(keyword, __LOCATION__, name="MAX_FORCE", &
171                          description="Tolerance on the maximum value of Forces on the BAND.", &
172                          usage="MAX_FORCE {real}", &
173                          default_r_val=0.00045_dp)
174      CALL section_add_keyword(subsection, keyword)
175      CALL keyword_release(keyword)
176
177      CALL keyword_create(keyword, __LOCATION__, name="RMS_DR", &
178                          description="Tolerance on RMS displacements on the BAND.", &
179                          usage="RMS_DR {real}", &
180                          default_r_val=0.0001_dp)
181      CALL section_add_keyword(subsection, keyword)
182      CALL keyword_release(keyword)
183
184      CALL keyword_create(keyword, __LOCATION__, name="RMS_FORCE", &
185                          description="Tolerance on RMS Forces on the BAND.", &
186                          usage="RMS_FORCE {real}", &
187                          default_r_val=0.00030_dp)
188      CALL section_add_keyword(subsection, keyword)
189      CALL keyword_release(keyword)
190      CALL section_add_subsection(section, subsection)
191      CALL section_release(subsection)
192
193      NULLIFY (subsection, subsubsection)
194      ! CI-NEB section
195      CALL section_create(subsection, __LOCATION__, name="CI_NEB", &
196                          description="Controls parameters for CI-NEB type calculation only.", &
197                          repeats=.FALSE.)
198      CALL keyword_create(keyword, __LOCATION__, name="NSTEPS_IT", &
199                          description="Specify the number of steps of IT-NEB to perform before "// &
200                          "switching on the CI algorithm", &
201                          default_i_val=5)
202      CALL section_add_keyword(subsection, keyword)
203      CALL keyword_release(keyword)
204      CALL section_add_subsection(section, subsection)
205      CALL section_release(subsection)
206
207      ! String Method section
208      CALL section_create(subsection, __LOCATION__, name="STRING_METHOD", &
209                          description="Controls parameters for String Method type calculation only.", &
210                          repeats=.FALSE.)
211
212      CALL keyword_create(keyword, __LOCATION__, name="SPLINE_ORDER", &
213                          description="Specify the oder of the spline used in the String Method.", &
214                          default_i_val=1)
215      CALL section_add_keyword(subsection, keyword)
216      CALL keyword_release(keyword)
217      CALL keyword_create(keyword, __LOCATION__, name="SMOOTHING", &
218                          description="Smoothing parameter for the reparametrization of the frames.", &
219                          default_r_val=0.2_dp)
220      CALL section_add_keyword(subsection, keyword)
221      CALL keyword_release(keyword)
222
223      CALL section_add_subsection(section, subsection)
224      CALL section_release(subsection)
225
226      ! Optimization section
227      CALL section_create(subsection, __LOCATION__, name="optimize_band", &
228                          description="Specify the optimization method for the band", &
229                          repeats=.TRUE.)
230      CALL create_opt_band_section(subsection)
231      CALL section_add_subsection(section, subsection)
232      CALL section_release(subsection)
233
234      ! replica section: to specify coordinates and velocities (possibly) of the
235      ! different replica used in the BAND
236      CALL section_create(subsection, __LOCATION__, name="replica", &
237                          description="Specify coordinates and velocities (possibly) of the replica", &
238                          repeats=.TRUE.)
239      ! Colvar
240      CALL keyword_create(keyword, __LOCATION__, name="COLLECTIVE", &
241                          description="Specifies the value of the collective variables used in the projected"// &
242                          " BAND method. The order of the values is the order of the COLLECTIVE section in the"// &
243                          " constraints/restraints section", &
244                          usage="COLLECTIVE {real} .. {real}", &
245                          type_of_var=real_t, n_var=-1)
246      CALL section_add_keyword(subsection, keyword)
247      CALL keyword_release(keyword)
248      ! Coordinates read through an external file
249      CALL keyword_create(keyword, __LOCATION__, name="COORD_FILE_NAME", &
250                          description="Name of the xyz file with coordinates (alternative to &COORD section)", &
251                          usage="COORD_FILE_NAME <CHAR>", &
252                          default_lc_val="")
253      CALL section_add_keyword(subsection, keyword)
254      CALL keyword_release(keyword)
255      ! Coordinates and velocities
256      CALL create_coord_section(subsubsection, "BAND")
257      CALL section_add_subsection(subsection, subsubsection)
258      CALL section_release(subsubsection)
259      CALL create_velocity_section(subsubsection, "BAND")
260      CALL section_add_subsection(subsection, subsubsection)
261      CALL section_release(subsubsection)
262
263      CALL section_add_subsection(section, subsection)
264      CALL section_release(subsection)
265
266      ! Print key section
267      CALL cp_print_key_section_create(print_key, __LOCATION__, "program_run_info", &
268                                       description="Controls the printing basic info about the BAND run", &
269                                       print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
270
271      CALL keyword_create(keyword, __LOCATION__, name="INITIAL_CONFIGURATION_INFO", &
272                          description="Print information for the setup of the initial configuration.", &
273                          usage="INITIAL_CONFIGURATION_INFO <LOGICAL>", &
274                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
275      CALL section_add_keyword(print_key, keyword)
276      CALL keyword_release(keyword)
277
278      CALL section_add_subsection(section, print_key)
279      CALL section_release(print_key)
280
281      CALL cp_print_key_section_create(print_key, __LOCATION__, "convergence_info", &
282                                       description="Controls the printing of the convergence criteria during a BAND run", &
283                                       print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
284      CALL section_add_subsection(section, print_key)
285      CALL section_release(print_key)
286
287      CALL cp_print_key_section_create(print_key, __LOCATION__, "replica_info", &
288                                       description="Controls the printing of each replica info during a BAND run", &
289                                       print_level=medium_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
290      CALL section_add_subsection(section, print_key)
291      CALL section_release(print_key)
292
293      CALL cp_print_key_section_create(print_key, __LOCATION__, "ENERGY", &
294                                       description="Controls the printing of the ENER file in a BAND run", &
295                                       print_level=low_print_level, common_iter_levels=1, &
296                                       filename="")
297      CALL section_add_subsection(section, print_key)
298      CALL section_release(print_key)
299
300      CALL cp_print_key_section_create(print_key, __LOCATION__, "BANNER", &
301                                       description="Controls the printing of the BAND banner", &
302                                       print_level=low_print_level, common_iter_levels=1, &
303                                       filename="__STD_OUT__")
304      CALL section_add_subsection(section, print_key)
305      CALL section_release(print_key)
306   END SUBROUTINE create_band_section
307
308! **************************************************************************************************
309!> \brief creates the optimization section for a BAND run
310!> \param section will contain the pint section
311!> \author Teodoro Laino 02.2007 [tlaino]
312! **************************************************************************************************
313   SUBROUTINE create_opt_band_section(section)
314      TYPE(section_type), POINTER                        :: section
315
316      CHARACTER(len=*), PARAMETER :: routineN = 'create_opt_band_section', &
317         routineP = moduleN//':'//routineN
318
319      TYPE(keyword_type), POINTER                        :: keyword
320      TYPE(section_type), POINTER                        :: print_key, subsection, subsubsection
321
322      CPASSERT(ASSOCIATED(section))
323      NULLIFY (keyword, print_key, subsection, subsubsection)
324
325      CALL keyword_create(keyword, __LOCATION__, name="OPT_TYPE", &
326                          description="Specifies the type optimizer used for the band", &
327                          usage="OPT_TYPE (MD|DIIS)", &
328                          default_i_val=band_diis_opt, &
329                          enum_c_vals=s2a("MD", &
330                                          "DIIS"), &
331                          enum_desc=s2a("Molecular dynamics-based optimizer", &
332                                        "Coupled steepest descent / direct inversion in the iterative subspace"), &
333                          enum_i_vals=(/band_md_opt, band_diis_opt/))
334      CALL section_add_keyword(section, keyword)
335      CALL keyword_release(keyword)
336
337      CALL keyword_create(keyword, __LOCATION__, name="OPTIMIZE_END_POINTS", &
338                          description="Performs also an optimization of the end points of the band.", &
339                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
340      CALL section_add_keyword(section, keyword)
341      CALL keyword_release(keyword)
342
343      ! MD optimization section
344      CALL section_create(subsection, __LOCATION__, name="MD", &
345                          description="Activate the MD based optimization procedure for BAND", &
346                          repeats=.FALSE.)
347
348      CALL keyword_create(keyword, __LOCATION__, name="MAX_STEPS", &
349                          description="Specify the maximum number of MD steps", &
350                          default_i_val=100)
351      CALL section_add_keyword(subsection, keyword)
352      CALL keyword_release(keyword)
353
354      CALL keyword_create( &
355         keyword, __LOCATION__, &
356         name="timestep", &
357         description="The length of an integration step", &
358         usage="timestep 1.0", &
359         default_r_val=cp_unit_to_cp2k(value=0.5_dp, &
360                                       unit_str="fs"), &
361         unit_str="fs")
362      CALL section_add_keyword(subsection, keyword)
363      CALL keyword_release(keyword)
364
365      CALL keyword_create(keyword, __LOCATION__, name="TEMPERATURE", &
366                          description="Specify the initial temperature", &
367                          default_r_val=cp_unit_to_cp2k(value=0.0_dp, &
368                                                        unit_str="K"), &
369                          unit_str="K")
370      CALL section_add_keyword(subsection, keyword)
371      CALL keyword_release(keyword)
372
373      ! Temp_control
374      CALL section_create(subsubsection, __LOCATION__, name="TEMP_CONTROL", &
375                          description="Setup parameters to control the temperature during a BAND MD run.", &
376                          repeats=.FALSE.)
377      CALL keyword_create(keyword, __LOCATION__, name="TEMPERATURE", &
378                          description="Specify the target temperature", &
379                          type_of_var=real_t, unit_str="K")
380      CALL section_add_keyword(subsubsection, keyword)
381      CALL keyword_release(keyword)
382
383      CALL keyword_create(keyword, __LOCATION__, name="TEMP_TOL", &
384                          description="Specify the tolerance on the temperature for rescaling", &
385                          default_r_val=cp_unit_to_cp2k(value=0.0_dp, &
386                                                        unit_str="K"), &
387                          unit_str="K")
388      CALL section_add_keyword(subsubsection, keyword)
389      CALL keyword_release(keyword)
390
391      CALL keyword_create(keyword, __LOCATION__, name="TEMP_TOL_STEPS", &
392                          description="Specify the number of steps to apply a temperature control", &
393                          default_i_val=0)
394      CALL section_add_keyword(subsubsection, keyword)
395      CALL keyword_release(keyword)
396      CALL section_add_subsection(subsection, subsubsection)
397      CALL section_release(subsubsection)
398
399      ! Vel_control
400      CALL section_create(subsubsection, __LOCATION__, name="VEL_CONTROL", &
401                          description="Setup parameters to control the velocity during a BAND MD run.", &
402                          repeats=.FALSE.)
403      CALL keyword_create(keyword, __LOCATION__, name="ANNEALING", &
404                          description="Specify the annealing coefficient", &
405                          default_r_val=1.0_dp)
406      CALL section_add_keyword(subsubsection, keyword)
407      CALL keyword_release(keyword)
408      CALL keyword_create(keyword, __LOCATION__, name="PROJ_VELOCITY_VERLET", &
409                          description="Uses a Projected Velocity Verlet instead of a normal Velocity Verlet."// &
410                          " Every time the cosine between velocities and forces is < 0 velocities are"// &
411                          " zeroed.", &
412                          usage="PROJ_VELOCITY_VERLET <LOGICAL>", &
413                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
414      CALL section_add_keyword(subsubsection, keyword)
415      CALL keyword_release(keyword)
416      CALL keyword_create(keyword, __LOCATION__, name="SD_LIKE", &
417                          description="Zeros velocity at each MD step emulating a steepest descent like"// &
418                          "(SD_LIKE) approach", &
419                          usage="SD_LIKE <LOGICAL>", &
420                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
421      CALL section_add_keyword(subsubsection, keyword)
422      CALL keyword_release(keyword)
423      CALL section_add_subsection(subsection, subsubsection)
424      CALL section_release(subsubsection)
425      ! End of MD
426      CALL section_add_subsection(section, subsection)
427      CALL section_release(subsection)
428
429      ! DIIS optimization section
430      CALL section_create(subsection, __LOCATION__, name="DIIS", &
431                          description="Activate the DIIS based optimization procedure for BAND", &
432                          repeats=.FALSE.)
433
434      CALL keyword_create(keyword, __LOCATION__, name="MAX_SD_STEPS", &
435                          description="Specify the maximum number of SD steps to perform"// &
436                          " before switching on DIIS (the minimum number will always be equal to N_DIIS).", &
437                          default_i_val=1)
438      CALL section_add_keyword(subsection, keyword)
439      CALL keyword_release(keyword)
440
441      CALL keyword_create(keyword, __LOCATION__, name="MAX_STEPS", &
442                          description="Specify the maximum number of optimization steps", &
443                          default_i_val=100)
444      CALL section_add_keyword(subsection, keyword)
445      CALL keyword_release(keyword)
446
447      CALL keyword_create(keyword, __LOCATION__, name="N_DIIS", &
448                          variants=(/"NDIIS"/), &
449                          description="Number of history vectors to be used with DIIS", &
450                          usage="N_DIIS 4", &
451                          default_i_val=7)
452      CALL section_add_keyword(subsection, keyword)
453      CALL keyword_release(keyword)
454
455      CALL keyword_create(keyword, __LOCATION__, name="STEPSIZE", &
456                          description="Initial stepsize used for the line search, sometimes this parameter"// &
457                          "can be reduced to stabilize DIIS", &
458                          usage="STEPSIZE <REAL>", &
459                          default_r_val=1.0_dp)
460      CALL section_add_keyword(subsection, keyword)
461      CALL keyword_release(keyword)
462
463      CALL keyword_create(keyword, __LOCATION__, name="MAX_STEPSIZE", &
464                          description="Maximum stepsize used for the line search, sometimes this parameter"// &
465                          "can be reduced to stabilize the LS for particularly difficult initial geometries", &
466                          usage="MAX_STEPSIZE <REAL>", &
467                          default_r_val=2.0_dp)
468      CALL section_add_keyword(subsection, keyword)
469      CALL keyword_release(keyword)
470
471      CALL keyword_create(keyword, __LOCATION__, name="NP_LS", &
472                          description="Number of points used in the line search SD.", &
473                          usage="NP_LS <INTEGER>", &
474                          default_i_val=2)
475      CALL section_add_keyword(subsection, keyword)
476      CALL keyword_release(keyword)
477
478      CALL keyword_create(keyword, __LOCATION__, name="NO_LS", &
479                          description="Does not perform LS during SD. Useful in combination with a proper STEPSIZE"// &
480                          " for particularly out of equilibrium starting geometries.", &
481                          default_l_val=.FALSE., lone_keyword_l_val=.TRUE.)
482      CALL section_add_keyword(subsection, keyword)
483      CALL keyword_release(keyword)
484
485      CALL keyword_create(keyword, __LOCATION__, name="CHECK_DIIS", &
486                          description="Performs a series of checks on the DIIS solution in order to accept the DIIS step."// &
487                          " If set to .FALSE. the only check performed is that the angle between the DIIS solution and the"// &
488                          " reference vector is less than Pi/2. Can be useful if many DIIS steps are rejected.", &
489                          default_l_val=.TRUE., lone_keyword_l_val=.TRUE.)
490      CALL section_add_keyword(subsection, keyword)
491      CALL keyword_release(keyword)
492
493      CALL cp_print_key_section_create(print_key, __LOCATION__, "diis_info", &
494                                       description="Controls the printing of DIIS info during a BAND run", &
495                                       print_level=high_print_level, add_last=add_last_numeric, filename="__STD_OUT__")
496      CALL section_add_subsection(subsection, print_key)
497      CALL section_release(print_key)
498
499      CALL section_add_subsection(section, subsection)
500      CALL section_release(subsection)
501   END SUBROUTINE create_opt_band_section
502
503END MODULE input_cp2k_neb
504