1!--------------------------------------------------------------------------------------------------!
2!   CP2K: A general program to perform molecular dynamics simulations                              !
3!   Copyright (C) 2000 - 2019  CP2K developers group                                               !
4!--------------------------------------------------------------------------------------------------!
5
6! **************************************************************************************************
7!> \brief interface to use cp2k as library
8!> \note
9!>      useful additions for the future would be:
10!>      - string(path) based set/get of simple values (to change the new
11!>        input during the run and extract more data (energy types for example).
12!>      - set/get of a subset of atoms
13!> \par History
14!>      07.2004 created [fawzi]
15!>      11.2004 parallel version [fawzi]
16!> \author fawzi & Johanna
17! **************************************************************************************************
18MODULE f77_interface
19   USE base_hooks,                      ONLY: cp_abort_hook,&
20                                              cp_warn_hook,&
21                                              timeset_hook,&
22                                              timestop_hook
23   USE bibliography,                    ONLY: add_all_references
24   USE cell_types,                      ONLY: cell_type,&
25                                              init_cell
26   USE cp2k_info,                       ONLY: get_runtime_info
27   USE cp_error_handling,               ONLY: cp_error_handling_setup
28   USE cp_files,                        ONLY: init_preconnection_list,&
29                                              open_file
30   USE cp_log_handling,                 ONLY: &
31        cp_add_default_logger, cp_default_logger_stack_size, cp_failure_level, &
32        cp_get_default_logger, cp_logger_create, cp_logger_get_default_unit_nr, cp_logger_release, &
33        cp_logger_retain, cp_logger_type, cp_rm_default_logger, cp_to_string
34   USE cp_output_handling,              ONLY: cp_iterate
35   USE cp_para_env,                     ONLY: cp_para_env_create,&
36                                              cp_para_env_release,&
37                                              cp_para_env_retain
38   USE cp_para_types,                   ONLY: cp_para_env_type
39   USE cp_result_methods,               ONLY: get_results,&
40                                              test_for_result
41   USE cp_result_types,                 ONLY: cp_result_type
42   USE cp_subsys_types,                 ONLY: cp_subsys_get,&
43                                              cp_subsys_set,&
44                                              cp_subsys_type,&
45                                              unpack_subsys_particles
46   USE dbcsr_api,                       ONLY: dbcsr_acc_get_ndevices,&
47                                              dbcsr_acc_set_active_device,&
48                                              dbcsr_finalize_lib,&
49                                              dbcsr_init_lib
50   USE eip_environment,                 ONLY: eip_init
51   USE eip_environment_types,           ONLY: eip_env_create,&
52                                              eip_env_release,&
53                                              eip_environment_type
54   USE embed_main,                      ONLY: embed_create_force_env
55   USE embed_types,                     ONLY: embed_env_release,&
56                                              embed_env_type
57   USE environment,                     ONLY: cp2k_finalize,&
58                                              cp2k_init,&
59                                              cp2k_read,&
60                                              cp2k_setup
61   USE fist_main,                       ONLY: fist_create_force_env
62   USE force_env_methods,               ONLY: force_env_calc_energy_force,&
63                                              force_env_create
64   USE force_env_types,                 ONLY: &
65        force_env_get, force_env_get_frc, force_env_get_natom, force_env_get_nparticle, &
66        force_env_get_pos, force_env_get_vel, force_env_release, force_env_retain, force_env_set, &
67        force_env_type, multiple_fe_list
68   USE fp_types,                        ONLY: fp_env_create,&
69                                              fp_env_read,&
70                                              fp_env_release,&
71                                              fp_env_write,&
72                                              fp_type
73   USE global_types,                    ONLY: global_environment_type,&
74                                              globenv_create,&
75                                              globenv_release,&
76                                              globenv_retain
77   USE input_constants,                 ONLY: do_eip,&
78                                              do_embed,&
79                                              do_fist,&
80                                              do_mixed,&
81                                              do_qmmm,&
82                                              do_qmmmx,&
83                                              do_qs,&
84                                              do_sirius
85   USE input_cp2k_check,                ONLY: check_cp2k_input
86   USE input_cp2k_force_eval,           ONLY: create_force_eval_section
87   USE input_cp2k_read,                 ONLY: empty_initial_variables,&
88                                              read_input
89   USE input_enumeration_types,         ONLY: enum_i2c,&
90                                              enumeration_type
91   USE input_keyword_types,             ONLY: keyword_get,&
92                                              keyword_type
93   USE input_section_types,             ONLY: &
94        section_get_keyword, section_release, section_type, section_vals_duplicate, &
95        section_vals_get, section_vals_get_subs_vals, section_vals_release, &
96        section_vals_remove_values, section_vals_retain, section_vals_type, section_vals_val_get, &
97        section_vals_write
98   USE kinds,                           ONLY: default_path_length,&
99                                              default_string_length,&
100                                              dp
101   USE machine,                         ONLY: default_output_unit,&
102                                              m_chdir,&
103                                              m_getcwd,&
104                                              m_memory
105   USE message_passing,                 ONLY: &
106        add_mp_perf_env, get_mp_perf_env, mp_comm_world, mp_max, mp_perf_env_release, &
107        mp_perf_env_retain, mp_perf_env_type, mp_world_finalize, mp_world_init, rm_mp_perf_env
108   USE metadynamics_types,              ONLY: meta_env_release,&
109                                              meta_env_type
110   USE metadynamics_utils,              ONLY: metadyn_read
111   USE mixed_environment_types,         ONLY: mixed_env_release,&
112                                              mixed_environment_type
113   USE mixed_main,                      ONLY: mixed_create_force_env
114   USE periodic_table,                  ONLY: init_periodic_table
115   USE pw_cuda,                         ONLY: pw_cuda_finalize,&
116                                              pw_cuda_init
117   USE pwdft_environment,               ONLY: pwdft_init
118   USE pwdft_environment_types,         ONLY: pwdft_env_create,&
119                                              pwdft_env_release,&
120                                              pwdft_environment_type
121   USE qmmm_create,                     ONLY: qmmm_env_create
122   USE qmmm_types,                      ONLY: qmmm_env_release,&
123                                              qmmm_env_type
124   USE qmmmx_create,                    ONLY: qmmmx_env_create
125   USE qmmmx_types,                     ONLY: qmmmx_env_release,&
126                                              qmmmx_env_type
127   USE qs_environment,                  ONLY: qs_init
128   USE qs_environment_types,            ONLY: qs_env_create,&
129                                              qs_env_release,&
130                                              qs_environment_type
131   USE reference_manager,               ONLY: remove_all_references
132   USE string_table,                    ONLY: string_table_allocate,&
133                                              string_table_deallocate
134   USE timings,                         ONLY: add_timer_env,&
135                                              get_timer_env,&
136                                              rm_timer_env,&
137                                              timer_env_release,&
138                                              timer_env_retain,&
139                                              timings_register_hooks
140   USE timings_types,                   ONLY: timer_env_type
141   USE virial_types,                    ONLY: virial_type
142#include "./base/base_uses.f90"
143
144   IMPLICIT NONE
145   PRIVATE
146
147   LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
148   CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'f77_interface'
149
150! **************************************************************************************************
151   TYPE f_env_p_type
152      TYPE(f_env_type), POINTER :: f_env
153   END TYPE f_env_p_type
154
155! **************************************************************************************************
156   TYPE f_env_type
157      INTEGER :: id_nr
158      TYPE(force_env_type), POINTER      :: force_env
159      TYPE(cp_logger_type), POINTER      :: logger
160      TYPE(timer_env_type), POINTER      :: timer_env
161      TYPE(mp_perf_env_type), POINTER    :: mp_perf_env
162      CHARACTER(len=default_path_length) :: my_path, old_path
163   END TYPE f_env_type
164
165   TYPE(f_env_p_type), DIMENSION(:), POINTER, SAVE :: f_envs
166   TYPE(cp_para_env_type), POINTER, SAVE :: default_para_env
167   LOGICAL, SAVE :: module_initialized = .FALSE.
168   INTEGER, SAVE :: last_f_env_id = 0, n_f_envs = 0
169
170   PUBLIC :: default_para_env
171   PUBLIC :: init_cp2k, finalize_cp2k
172   PUBLIC :: create_force_env, destroy_force_env, set_pos, get_pos, &
173             get_force, calc_energy_force, get_energy, get_stress_tensor, &
174             calc_energy, calc_force, check_input, get_natom, get_nparticle, &
175             f_env_add_defaults, f_env_rm_defaults, f_env_type, &
176             f_env_get_from_id, &
177             set_vel, set_cell, get_cell, get_result_r1
178CONTAINS
179
180! **************************************************************************************************
181!> \brief returns the position of the force env corresponding to the given id
182!> \param env_id the id of the requested environment
183!> \return ...
184!> \author fawzi
185!> \note
186!>      private utility function
187! **************************************************************************************************
188   FUNCTION get_pos_of_env(env_id) RESULT(res)
189      INTEGER, INTENT(in)                                :: env_id
190      INTEGER                                            :: res
191
192      INTEGER                                            :: env_pos, isub
193
194      env_pos = -1
195      DO isub = 1, n_f_envs
196         IF (f_envs(isub)%f_env%id_nr == env_id) THEN
197            env_pos = isub
198         END IF
199      END DO
200      res = env_pos
201   END FUNCTION get_pos_of_env
202
203! **************************************************************************************************
204!> \brief initializes cp2k, needs to be called once before using any of the
205!>      other functions when using cp2k as library
206!> \param init_mpi if the mpi environment should be initialized
207!> \param ierr returns a number different from 0 if there was an error
208!> \author fawzi
209! **************************************************************************************************
210   SUBROUTINE init_cp2k(init_mpi, ierr)
211      LOGICAL, INTENT(in)                                :: init_mpi
212      INTEGER, INTENT(out)                               :: ierr
213
214      CHARACTER(len=*), PARAMETER :: routineN = 'init_cp2k', routineP = moduleN//':'//routineN
215
216      INTEGER                                            :: mpi_comm_default, unit_nr
217      TYPE(cp_logger_type), POINTER                      :: logger
218
219      IF (.NOT. module_initialized) THEN
220         ! install error handler hooks
221         CALL cp_error_handling_setup()
222
223         ! install timming handler hooks
224         CALL timings_register_hooks()
225
226         ! Initialise preconnection list
227         CALL init_preconnection_list()
228
229         ! get runtime information
230         CALL get_runtime_info
231
232         IF (init_mpi) THEN
233            ! get the default system wide communicator
234            CALL mp_world_init(mpi_comm_default)
235         ELSE
236            mpi_comm_default = mp_comm_world
237         END IF
238
239         CALL string_table_allocate()
240         CALL add_mp_perf_env()
241         CALL add_timer_env()
242
243         ! re-create the para_env and log with correct (reordered) ranks
244         NULLIFY (default_para_env)
245         CALL cp_para_env_create(default_para_env, group=mpi_comm_default, &
246                                 owns_group=.FALSE.)
247         IF (default_para_env%ionode) THEN
248            unit_nr = default_output_unit
249         ELSE
250            unit_nr = -1
251         END IF
252         NULLIFY (logger)
253
254         CALL cp_logger_create(logger, para_env=default_para_env, &
255                               default_global_unit_nr=unit_nr, &
256                               close_global_unit_on_dealloc=.FALSE.)
257         CALL cp_add_default_logger(logger)
258         CALL cp_logger_release(logger)
259
260         ALLOCATE (f_envs(0))
261         module_initialized = .TRUE.
262         ierr = 0
263
264         !   *** Initialize mathematical constants ***
265         CALL init_periodic_table()
266
267         !   *** init the bibliography ***
268         CALL add_all_references()
269
270         IF (dbcsr_acc_get_ndevices() > 0) &
271            CALL dbcsr_acc_set_active_device(MOD(default_para_env%mepos, dbcsr_acc_get_ndevices()))
272
273         CALL pw_cuda_init()
274
275         ! Initialize the DBCSR configuration
276         ! Attach the time handler hooks to DBCSR
277         CALL dbcsr_init_lib(default_para_env%group, timeset_hook, timestop_hook, &
278                             cp_abort_hook, cp_warn_hook, io_unit=unit_nr)
279      ELSE
280         ierr = cp_failure_level
281      END IF
282
283      !sample peak memory
284      CALL m_memory()
285
286   END SUBROUTINE init_cp2k
287
288! **************************************************************************************************
289!> \brief cleanup after you have finished using this interface
290!> \param finalize_mpi if the mpi environment should be finalized
291!> \param ierr returns a number different from 0 if there was an error
292!> \author fawzi
293! **************************************************************************************************
294   SUBROUTINE finalize_cp2k(finalize_mpi, ierr)
295      LOGICAL, INTENT(in)                                :: finalize_mpi
296      INTEGER, INTENT(out)                               :: ierr
297
298      CHARACTER(len=*), PARAMETER :: routineN = 'finalize_cp2k', routineP = moduleN//':'//routineN
299
300      INTEGER                                            :: ienv
301
302!sample peak memory
303
304      CALL m_memory()
305
306      IF (.NOT. module_initialized) THEN
307         ierr = cp_failure_level
308      ELSE
309         ! Finalize the DBCSR configuration
310         CALL dbcsr_finalize_lib()
311         CALL pw_cuda_finalize()
312         DO ienv = n_f_envs, 1, -1
313            CALL destroy_force_env(f_envs(ienv)%f_env%id_nr, ierr=ierr)
314            CPASSERT(ierr == 0)
315         END DO
316         DEALLOCATE (f_envs)
317         CALL cp_para_env_release(default_para_env)
318         ierr = 0
319         CALL cp_rm_default_logger()
320
321         !   *** deallocate the bibliography ***
322         CALL remove_all_references()
323         CALL rm_timer_env()
324         CALL rm_mp_perf_env()
325         CALL string_table_deallocate(0)
326         IF (finalize_mpi) THEN
327            CALL mp_world_finalize()
328         END IF
329      END IF
330   END SUBROUTINE finalize_cp2k
331
332! **************************************************************************************************
333!> \brief deallocates a f_env
334!> \param f_env the f_env to deallocate
335!> \author fawzi
336! **************************************************************************************************
337   RECURSIVE SUBROUTINE f_env_dealloc(f_env)
338      TYPE(f_env_type), POINTER                          :: f_env
339
340      CHARACTER(len=*), PARAMETER :: routineN = 'f_env_dealloc', routineP = moduleN//':'//routineN
341
342      INTEGER                                            :: ierr
343
344      CPASSERT(ASSOCIATED(f_env))
345      CALL force_env_release(f_env%force_env)
346      CALL cp_logger_release(f_env%logger)
347      CALL timer_env_release(f_env%timer_env)
348      CALL mp_perf_env_release(f_env%mp_perf_env)
349      IF (f_env%old_path /= f_env%my_path) THEN
350         CALL m_chdir(f_env%old_path, ierr)
351         CPASSERT(ierr == 0)
352      END IF
353   END SUBROUTINE f_env_dealloc
354
355! **************************************************************************************************
356!> \brief createates a f_env
357!> \param f_env the f_env to createate
358!> \param force_env the force_environment to be stored
359!> \param timer_env the timer env to be stored
360!> \param mp_perf_env the mp performance environment to be stored
361!> \param id_nr ...
362!> \param logger ...
363!> \param old_dir ...
364!> \author fawzi
365! **************************************************************************************************
366   SUBROUTINE f_env_create(f_env, force_env, timer_env, mp_perf_env, id_nr, logger, old_dir)
367      TYPE(f_env_type), POINTER                          :: f_env
368      TYPE(force_env_type), POINTER                      :: force_env
369      TYPE(timer_env_type), POINTER                      :: timer_env
370      TYPE(mp_perf_env_type), POINTER                    :: mp_perf_env
371      INTEGER, INTENT(in)                                :: id_nr
372      TYPE(cp_logger_type), POINTER                      :: logger
373      CHARACTER(len=*), INTENT(in)                       :: old_dir
374
375      CHARACTER(len=*), PARAMETER :: routineN = 'f_env_create', routineP = moduleN//':'//routineN
376
377      ALLOCATE (f_env)
378      f_env%force_env => force_env
379      CALL force_env_retain(f_env%force_env)
380      f_env%logger => logger
381      CALL cp_logger_retain(logger)
382      f_env%timer_env => timer_env
383      CALL timer_env_retain(f_env%timer_env)
384      f_env%mp_perf_env => mp_perf_env
385      CALL mp_perf_env_retain(f_env%mp_perf_env)
386      f_env%id_nr = id_nr
387      CALL m_getcwd(f_env%my_path)
388      f_env%old_path = old_dir
389   END SUBROUTINE f_env_create
390
391! **************************************************************************************************
392!> \brief ...
393!> \param f_env_id ...
394!> \param f_env ...
395! **************************************************************************************************
396   SUBROUTINE f_env_get_from_id(f_env_id, f_env)
397      INTEGER, INTENT(in)                                :: f_env_id
398      TYPE(f_env_type), POINTER                          :: f_env
399
400      CHARACTER(len=*), PARAMETER :: routineN = 'f_env_get_from_id', &
401         routineP = moduleN//':'//routineN
402
403      INTEGER                                            :: f_env_pos
404
405      NULLIFY (f_env)
406      f_env_pos = get_pos_of_env(f_env_id)
407      IF (f_env_pos < 1) THEN
408         CPABORT("invalid env_id "//cp_to_string(f_env_id))
409      ELSE
410         f_env => f_envs(f_env_pos)%f_env
411      END IF
412
413   END SUBROUTINE f_env_get_from_id
414
415! **************************************************************************************************
416!> \brief adds the default environments of the f_env to the stack of the
417!>      defaults, and returns a new error and sets failure to true if
418!>      something went wrong
419!> \param f_env_id the f_env from where to take the defaults
420!> \param f_env will contain the f_env corresponding to f_env_id
421!> \param handle ...
422!> \author fawzi
423!> \note
424!>      The following routines need to be synchronized wrt. adding/removing
425!>      of the default environments (logging, performance,error):
426!>      environment:cp2k_init, environment:cp2k_finalize,
427!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
428!>      f77_interface:create_force_env, f77_interface:destroy_force_env
429! **************************************************************************************************
430   SUBROUTINE f_env_add_defaults(f_env_id, f_env, handle)
431      INTEGER, INTENT(in)                                :: f_env_id
432      TYPE(f_env_type), POINTER                          :: f_env
433      INTEGER, INTENT(out), OPTIONAL                     :: handle
434
435      CHARACTER(len=*), PARAMETER :: routineN = 'f_env_add_defaults', &
436         routineP = moduleN//':'//routineN
437
438      INTEGER                                            :: f_env_pos, ierr
439      TYPE(cp_logger_type), POINTER                      :: logger
440
441      NULLIFY (f_env)
442      f_env_pos = get_pos_of_env(f_env_id)
443      IF (f_env_pos < 1) THEN
444         CPABORT("invalid env_id "//cp_to_string(f_env_id))
445      ELSE
446         f_env => f_envs(f_env_pos)%f_env
447         logger => f_env%logger
448         CPASSERT(ASSOCIATED(logger))
449         CALL m_getcwd(f_env%old_path)
450         IF (f_env%old_path /= f_env%my_path) THEN
451            CALL m_chdir(TRIM(f_env%my_path), ierr)
452            CPASSERT(ierr == 0)
453         END IF
454         CALL add_mp_perf_env(f_env%mp_perf_env)
455         CALL add_timer_env(f_env%timer_env)
456         CALL cp_add_default_logger(logger)
457         IF (PRESENT(handle)) handle = cp_default_logger_stack_size()
458      END IF
459   END SUBROUTINE f_env_add_defaults
460
461! **************************************************************************************************
462!> \brief removes the default environments of the f_env to the stack of the
463!>      defaults, and sets ierr accordingly to the failuers stored in error
464!>      It also releases the error
465!> \param f_env the f_env from where to take the defaults
466!> \param ierr variable that will be set to a number different from 0 if
467!>        error contains an error (otherwise it will be set to 0)
468!> \param handle ...
469!> \author fawzi
470!> \note
471!>      The following routines need to be synchronized wrt. adding/removing
472!>      of the default environments (logging, performance,error):
473!>      environment:cp2k_init, environment:cp2k_finalize,
474!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
475!>      f77_interface:create_force_env, f77_interface:destroy_force_env
476! **************************************************************************************************
477   SUBROUTINE f_env_rm_defaults(f_env, ierr, handle)
478      TYPE(f_env_type), POINTER                          :: f_env
479      INTEGER, INTENT(out), OPTIONAL                     :: ierr
480      INTEGER, INTENT(in), OPTIONAL                      :: handle
481
482      CHARACTER(len=*), PARAMETER :: routineN = 'f_env_rm_defaults', &
483         routineP = moduleN//':'//routineN
484
485      INTEGER                                            :: ierr2
486      TYPE(cp_logger_type), POINTER                      :: d_logger, logger
487      TYPE(mp_perf_env_type), POINTER                    :: d_mp_perf_env
488      TYPE(timer_env_type), POINTER                      :: d_timer_env
489
490      IF (ASSOCIATED(f_env)) THEN
491         IF (PRESENT(handle)) THEN
492            CPASSERT(handle == cp_default_logger_stack_size())
493         END IF
494
495         logger => f_env%logger
496         d_logger => cp_get_default_logger()
497         d_timer_env => get_timer_env()
498         d_mp_perf_env => get_mp_perf_env()
499         CPASSERT(ASSOCIATED(logger))
500         CPASSERT(ASSOCIATED(d_logger))
501         CPASSERT(ASSOCIATED(d_timer_env))
502         CPASSERT(ASSOCIATED(d_mp_perf_env))
503         CPASSERT(logger%id_nr == d_logger%id_nr)
504         ! CPASSERT(d_timer_env%id_nr==f_env%timer_env%id_nr)
505         CPASSERT(d_mp_perf_env%id_nr == f_env%mp_perf_env%id_nr)
506         IF (f_env%old_path /= f_env%my_path) THEN
507            CALL m_chdir(TRIM(f_env%old_path), ierr2)
508            CPASSERT(ierr2 == 0)
509         END IF
510         IF (PRESENT(ierr)) THEN
511            ierr = 0
512         ENDIF
513         CALL cp_rm_default_logger()
514         CALL rm_timer_env()
515         CALL rm_mp_perf_env()
516      ELSE
517         IF (PRESENT(ierr)) THEN
518            ierr = 0
519         ENDIF
520      END IF
521   END SUBROUTINE f_env_rm_defaults
522
523! **************************************************************************************************
524!> \brief creates a new force environment using the given input, and writing
525!>      the output to the given output unit
526!> \param new_env_id will contain the id of the newly created environment
527!> \param input_declaration ...
528!> \param input_path where to read the input (if the input is given it can
529!>        a virtual path)
530!> \param output_path filename (or name of the unit) for the output
531!> \param mpi_comm the mpi communicator to be used for this environment
532!>        it will not be freed when you get rid of the force_env
533!> \param output_unit if given it should be the unit for the output
534!>        and no file is open(should be valid on the processor with rank 0)
535!> \param owns_out_unit if the output unit should be closed upon destroing
536!>        of the force_env (defaults to true if not default_output_unit)
537!> \param input the parsed input, if given and valid it is used
538!>        instead of parsing from file
539!> \param ierr will return a number different from 0 if there was an error
540!> \param work_dir ...
541!> \param initial_variables ...
542!> \author fawzi
543!> \note
544!>      The following routines need to be synchronized wrt. adding/removing
545!>      of the default environments (logging, performance,error):
546!>      environment:cp2k_init, environment:cp2k_finalize,
547!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
548!>      f77_interface:create_force_env, f77_interface:destroy_force_env
549! **************************************************************************************************
550   RECURSIVE SUBROUTINE create_force_env(new_env_id, input_declaration, input_path, &
551                                         output_path, mpi_comm, output_unit, owns_out_unit, &
552                                         input, ierr, work_dir, initial_variables)
553      INTEGER, INTENT(out)                               :: new_env_id
554      TYPE(section_type), POINTER                        :: input_declaration
555      CHARACTER(len=*), INTENT(in)                       :: input_path
556      CHARACTER(len=*), INTENT(in), OPTIONAL             :: output_path
557      INTEGER, INTENT(in), OPTIONAL                      :: mpi_comm, output_unit
558      LOGICAL, INTENT(in), OPTIONAL                      :: owns_out_unit
559      TYPE(section_vals_type), OPTIONAL, POINTER         :: input
560      INTEGER, INTENT(out), OPTIONAL                     :: ierr
561      CHARACTER(len=*), INTENT(in), OPTIONAL             :: work_dir
562      CHARACTER(len=*), DIMENSION(:, :), OPTIONAL        :: initial_variables
563
564      CHARACTER(len=*), PARAMETER :: routineN = 'create_force_env', &
565         routineP = moduleN//':'//routineN
566
567      CHARACTER(len=default_path_length)                 :: old_dir, wdir
568      INTEGER :: handle, i, ierr2, iforce_eval, isubforce_eval, k, method_name_id, my_group, &
569         nforce_eval, ngroups, nsubforce_size, unit_nr
570      INTEGER, DIMENSION(:), POINTER                     :: group_distribution, i_force_eval, &
571                                                            lgroup_distribution
572      LOGICAL :: check, do_qmmm_force_mixing, multiple_subsys, my_echo, my_owns_out_unit, &
573         use_motion_section, use_multiple_para_env
574      TYPE(cp_logger_type), POINTER                      :: logger, my_logger
575      TYPE(cp_para_env_type), POINTER                    :: my_para_env, para_env
576      TYPE(eip_environment_type), POINTER                :: eip_env
577      TYPE(embed_env_type), POINTER                      :: embed_env
578      TYPE(enumeration_type), POINTER                    :: enum
579      TYPE(f_env_p_type), DIMENSION(:), POINTER          :: f_envs_old
580      TYPE(force_env_type), POINTER                      :: force_env, my_force_env
581      TYPE(fp_type), POINTER                             :: fp_env
582      TYPE(global_environment_type), POINTER             :: globenv
583      TYPE(keyword_type), POINTER                        :: keyword
584      TYPE(meta_env_type), POINTER                       :: meta_env
585      TYPE(mixed_environment_type), POINTER              :: mixed_env
586      TYPE(mp_perf_env_type), POINTER                    :: mp_perf_env
587      TYPE(pwdft_environment_type), POINTER              :: pwdft_env
588      TYPE(qmmm_env_type), POINTER                       :: qmmm_env
589      TYPE(qmmmx_env_type), POINTER                      :: qmmmx_env
590      TYPE(qs_environment_type), POINTER                 :: qs_env
591      TYPE(section_type), POINTER                        :: section
592      TYPE(section_vals_type), POINTER :: fe_section, force_env_section, force_env_sections, &
593         fp_section, input_file, qmmm_section, qmmmx_section, root_section, subsys_section, &
594         wrk_section
595      TYPE(timer_env_type), POINTER                      :: timer_env
596
597      CPASSERT(ASSOCIATED(input_declaration))
598      NULLIFY (para_env, force_env, timer_env, mp_perf_env, globenv, meta_env, &
599               fp_env, eip_env, pwdft_env, mixed_env, qs_env, qmmm_env, embed_env)
600      new_env_id = -1
601      IF (PRESENT(mpi_comm)) THEN
602         CALL cp_para_env_create(para_env, group=mpi_comm, owns_group=.FALSE.)
603      ELSE
604         para_env => default_para_env
605         CALL cp_para_env_retain(para_env)
606      END IF
607
608      CALL timeset(routineN, handle)
609
610      CALL m_getcwd(old_dir)
611      wdir = old_dir
612      IF (PRESENT(work_dir)) THEN
613         IF (work_dir /= " ") THEN
614            CALL m_chdir(work_dir, ierr2)
615            IF (ierr2 /= 0) THEN
616               IF (PRESENT(ierr)) ierr = ierr2
617               RETURN
618            END IF
619            wdir = work_dir
620         END IF
621      END IF
622
623      IF (PRESENT(output_unit)) THEN
624         unit_nr = output_unit
625      ELSE
626         IF (para_env%ionode) THEN
627            IF (output_path == "__STD_OUT__") THEN
628               unit_nr = default_output_unit
629            ELSE
630               CALL open_file(file_name=output_path, file_status="UNKNOWN", &
631                              file_action="WRITE", file_position="APPEND", &
632                              unit_number=unit_nr)
633            END IF
634         ELSE
635            unit_nr = -1
636         END IF
637      END IF
638      my_owns_out_unit = unit_nr /= default_output_unit
639      IF (PRESENT(owns_out_unit)) my_owns_out_unit = owns_out_unit
640      CALL globenv_create(globenv)
641      CALL cp2k_init(para_env, output_unit=unit_nr, globenv=globenv, input_file_name=input_path, &
642                     wdir=wdir)
643      logger => cp_get_default_logger()
644      ! warning this is dangerous, I did not check that all the subfunctions
645      ! support it, the program might crash upon error
646
647      NULLIFY (input_file)
648      IF (PRESENT(input)) input_file => input
649      IF (.NOT. ASSOCIATED(input_file)) THEN
650         IF (PRESENT(initial_variables)) THEN
651            input_file => read_input(input_declaration, input_path, initial_variables, para_env=para_env)
652         ELSE
653            input_file => read_input(input_declaration, input_path, empty_initial_variables, para_env=para_env)
654         ENDIF
655      ELSE
656         CALL section_vals_retain(input_file)
657      END IF
658      CALL section_vals_val_get(input_file, "GLOBAL%ECHO_INPUT", &
659                                l_val=my_echo)
660      ! echo after check?
661      IF (para_env%ionode .AND. my_echo) THEN
662         CALL section_vals_write(input_file, unit_nr=cp_logger_get_default_unit_nr(logger), &
663                                 hide_root=.TRUE., hide_defaults=.FALSE.)
664      END IF
665      ! XXXXXXXXXXXXXXXXXXXXXXXXXXX
666      ! root_section => input_file
667      ! XXXXXXXXXXXXXXXXXXXXXXXXXXX
668
669      CALL check_cp2k_input(input_declaration, input_file, para_env=para_env, output_unit=unit_nr)
670      ! XXXXXXXXXXXXXXXXXXXXXXXXXXX
671      ! NULLIFY(input_file)
672      ! XXXXXXXXXXXXXXXXXXXXXXXXXXX
673      root_section => input_file
674      CALL section_vals_retain(root_section)
675
676      IF (n_f_envs + 1 > SIZE(f_envs)) THEN
677         f_envs_old => f_envs
678         ALLOCATE (f_envs(n_f_envs + 10))
679         DO i = 1, n_f_envs
680            f_envs(i)%f_env => f_envs_old(i)%f_env
681         END DO
682         DO i = n_f_envs + 1, SIZE(f_envs)
683            NULLIFY (f_envs(i)%f_env)
684         END DO
685         DEALLOCATE (f_envs_old)
686      END IF
687
688      CALL cp2k_read(root_section, para_env, globenv)
689
690      CALL cp2k_setup(root_section, para_env, globenv)
691      ! Group Distribution
692      ALLOCATE (group_distribution(0:para_env%num_pe - 1))
693      group_distribution = 0
694      lgroup_distribution => group_distribution
695      ! Setup all possible force_env
696      force_env_sections => section_vals_get_subs_vals(root_section, "FORCE_EVAL")
697      CALL section_vals_val_get(root_section, "MULTIPLE_FORCE_EVALS%MULTIPLE_SUBSYS", &
698                                l_val=multiple_subsys)
699      CALL multiple_fe_list(force_env_sections, root_section, i_force_eval, nforce_eval)
700      ! Enforce the deletion of the subsys (unless not explicitly required)
701      IF (.NOT. multiple_subsys) THEN
702         DO iforce_eval = 2, nforce_eval
703            wrk_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
704                                                      i_rep_section=i_force_eval(iforce_eval))
705            CALL section_vals_remove_values(wrk_section)
706         END DO
707      END IF
708      nsubforce_size = nforce_eval - 1
709      use_multiple_para_env = .FALSE.
710      use_motion_section = .TRUE.
711      DO iforce_eval = 1, nforce_eval
712         NULLIFY (force_env_section, my_force_env, subsys_section)
713         ! Reference subsys from the first ordered force_eval
714         IF (.NOT. multiple_subsys) THEN
715            subsys_section => section_vals_get_subs_vals(force_env_sections, "SUBSYS", &
716                                                         i_rep_section=i_force_eval(1))
717         END IF
718         ! Handling para_env in case of multiple force_eval
719         IF (use_multiple_para_env) THEN
720            ! Check that the order of the force_eval is the correct one
721            CALL section_vals_val_get(force_env_sections, "METHOD", i_val=method_name_id, &
722                                      i_rep_section=i_force_eval(1))
723            IF ((method_name_id /= do_mixed) .AND. (method_name_id /= do_embed)) &
724               CALL cp_abort(__LOCATION__, &
725                             "In case of multiple force_eval the MAIN force_eval (the first in the list of FORCE_EVAL_ORDER or "// &
726                             "the one omitted from that order list) must be a MIXED_ENV type calculation. Please check your "// &
727                             "input file and possibly correct the MULTIPLE_FORCE_EVAL%FORCE_EVAL_ORDER. ")
728
729            IF (method_name_id .EQ. do_mixed) THEN
730               check = ASSOCIATED(force_env%mixed_env%sub_para_env)
731               CPASSERT(check)
732               ngroups = force_env%mixed_env%ngroups
733               my_group = lgroup_distribution(para_env%mepos)
734               isubforce_eval = iforce_eval - 1
735               ! If task not allocated on this procs skip setup..
736               IF (MODULO(isubforce_eval - 1, ngroups) /= my_group) CYCLE
737               my_para_env => force_env%mixed_env%sub_para_env(my_group + 1)%para_env
738               my_logger => force_env%mixed_env%sub_logger(my_group + 1)%p
739               CALL cp_rm_default_logger()
740               CALL cp_add_default_logger(my_logger)
741            ENDIF
742            IF (method_name_id .EQ. do_embed) THEN
743               check = ASSOCIATED(force_env%embed_env%sub_para_env)
744               CPASSERT(check)
745               ngroups = force_env%embed_env%ngroups
746               my_group = lgroup_distribution(para_env%mepos)
747               isubforce_eval = iforce_eval - 1
748               ! If task not allocated on this procs skip setup..
749               IF (MODULO(isubforce_eval - 1, ngroups) /= my_group) CYCLE
750               my_para_env => force_env%embed_env%sub_para_env(my_group + 1)%para_env
751               my_logger => force_env%embed_env%sub_logger(my_group + 1)%p
752               CALL cp_rm_default_logger()
753               CALL cp_add_default_logger(my_logger)
754            ENDIF
755         ELSE
756            my_para_env => para_env
757         END IF
758
759         ! Initialize force_env_section
760         ! No need to allocate one more force_env_section if only 1 force_eval
761         ! is provided.. this is in order to save memory..
762         IF (nforce_eval > 1) THEN
763            CALL section_vals_duplicate(force_env_sections, force_env_section, &
764                                        i_force_eval(iforce_eval), i_force_eval(iforce_eval))
765            IF (iforce_eval /= 1) use_motion_section = .FALSE.
766         ELSE
767            force_env_section => force_env_sections
768            use_motion_section = .TRUE.
769         END IF
770         CALL section_vals_val_get(force_env_section, "METHOD", i_val=method_name_id)
771
772         IF (method_name_id == do_qmmm) THEN
773            qmmmx_section => section_vals_get_subs_vals(force_env_section, "QMMM%FORCE_MIXING")
774            CALL section_vals_get(qmmmx_section, explicit=do_qmmm_force_mixing)
775            IF (do_qmmm_force_mixing) &
776               method_name_id = do_qmmmx ! QMMM Force-Mixing has its own (hidden) method_id
777         ENDIF
778
779         SELECT CASE (method_name_id)
780         CASE (do_fist)
781            CALL fist_create_force_env(my_force_env, root_section, my_para_env, globenv, &
782                                       force_env_section=force_env_section, subsys_section=subsys_section, &
783                                       use_motion_section=use_motion_section)
784
785         CASE (do_qs)
786            CALL qs_env_create(qs_env, globenv)
787            CALL qs_init(qs_env, my_para_env, root_section, globenv=globenv, force_env_section=force_env_section, &
788                         subsys_section=subsys_section, use_motion_section=use_motion_section)
789            CALL force_env_create(my_force_env, root_section, qs_env=qs_env, para_env=my_para_env, globenv=globenv, &
790                                  force_env_section=force_env_section)
791            CALL qs_env_release(qs_env)
792
793         CASE (do_qmmm)
794            qmmm_section => section_vals_get_subs_vals(force_env_section, "QMMM")
795            CALL qmmm_env_create(qmmm_env, root_section, my_para_env, globenv, &
796                                 force_env_section, qmmm_section, subsys_section, use_motion_section)
797            CALL force_env_create(my_force_env, root_section, qmmm_env=qmmm_env, para_env=my_para_env, &
798                                  globenv=globenv, force_env_section=force_env_section)
799            CALL qmmm_env_release(qmmm_env)
800
801         CASE (do_qmmmx)
802            CALL qmmmx_env_create(qmmmx_env, root_section, my_para_env, globenv, &
803                                  force_env_section, subsys_section, use_motion_section)
804            CALL force_env_create(my_force_env, root_section, qmmmx_env=qmmmx_env, para_env=my_para_env, &
805                                  globenv=globenv, force_env_section=force_env_section)
806            CALL qmmmx_env_release(qmmmx_env)
807
808         CASE (do_eip)
809            CALL eip_env_create(eip_env)
810            CALL eip_init(eip_env, root_section, my_para_env, force_env_section=force_env_section, &
811                          subsys_section=subsys_section)
812            CALL force_env_create(my_force_env, root_section, eip_env=eip_env, para_env=my_para_env, &
813                                  globenv=globenv, force_env_section=force_env_section)
814            CALL eip_env_release(eip_env)
815
816         CASE (do_sirius)
817            CALL pwdft_env_create(pwdft_env)
818            CALL pwdft_init(pwdft_env, root_section, my_para_env, force_env_section=force_env_section, &
819                            subsys_section=subsys_section, use_motion_section=use_motion_section)
820            CALL force_env_create(my_force_env, root_section, pwdft_env=pwdft_env, para_env=my_para_env, &
821                                  globenv=globenv, force_env_section=force_env_section)
822            CALL pwdft_env_release(pwdft_env)
823
824         CASE (do_mixed)
825            CALL mixed_create_force_env(mixed_env, root_section, my_para_env, &
826                                        force_env_section=force_env_section, n_subforce_eval=nsubforce_size, &
827                                        use_motion_section=use_motion_section)
828            CALL force_env_create(my_force_env, root_section, mixed_env=mixed_env, para_env=my_para_env, &
829                                  globenv=globenv, force_env_section=force_env_section)
830            CALL mixed_env_release(mixed_env)
831            !TODO: the sub_force_envs should really be created via recursion
832            use_multiple_para_env = .TRUE.
833            CALL cp_add_default_logger(logger) ! just to get the logger swapping started
834            lgroup_distribution => my_force_env%mixed_env%group_distribution
835
836         CASE (do_embed)
837            CALL embed_create_force_env(embed_env, root_section, my_para_env, &
838                                        force_env_section=force_env_section, n_subforce_eval=nsubforce_size, &
839                                        use_motion_section=use_motion_section)
840            CALL force_env_create(my_force_env, root_section, embed_env=embed_env, para_env=my_para_env, &
841                                  globenv=globenv, force_env_section=force_env_section)
842            CALL embed_env_release(embed_env)
843            !TODO: the sub_force_envs should really be created via recursion
844            use_multiple_para_env = .TRUE.
845            CALL cp_add_default_logger(logger) ! just to get the logger swapping started
846            lgroup_distribution => my_force_env%embed_env%group_distribution
847
848         CASE default
849            CALL create_force_eval_section(section)
850            keyword => section_get_keyword(section, "METHOD")
851            CALL keyword_get(keyword, enum=enum)
852            CALL cp_abort(__LOCATION__, &
853                          "Invalid METHOD <"//TRIM(enum_i2c(enum, method_name_id))// &
854                          "> was specified, ")
855            CALL section_release(section)
856         END SELECT
857
858         NULLIFY (meta_env, fp_env)
859         IF (use_motion_section) THEN
860            ! Metadynamics Setup
861            fe_section => section_vals_get_subs_vals(root_section, "MOTION%FREE_ENERGY")
862            CALL metadyn_read(meta_env, my_force_env, root_section, my_para_env, fe_section)
863            CALL force_env_set(my_force_env, meta_env=meta_env)
864            CALL meta_env_release(meta_env)
865            ! Flexible Partition Setup
866            fp_section => section_vals_get_subs_vals(root_section, "MOTION%FLEXIBLE_PARTITIONING")
867            CALL fp_env_create(fp_env)
868            CALL fp_env_read(fp_env, fp_section)
869            CALL fp_env_write(fp_env, fp_section)
870            CALL force_env_set(my_force_env, fp_env=fp_env)
871            CALL fp_env_release(fp_env)
872         END IF
873         ! Handle multiple force_eval
874         IF (nforce_eval > 1 .AND. iforce_eval == 1) THEN
875            ALLOCATE (my_force_env%sub_force_env(nsubforce_size))
876            ! Nullify subforce_env
877            DO k = 1, nsubforce_size
878               NULLIFY (my_force_env%sub_force_env(k)%force_env)
879            END DO
880         END IF
881         ! Reference the right force_env
882         IF (iforce_eval == 1) THEN
883            force_env => my_force_env
884         ELSE
885            force_env%sub_force_env(iforce_eval - 1)%force_env => my_force_env
886         END IF
887         ! Multiple para env for sub_force_eval
888         IF (.NOT. use_multiple_para_env) THEN
889            lgroup_distribution = iforce_eval
890         END IF
891         ! Release force_env_section
892         IF (nforce_eval > 1) CALL section_vals_release(force_env_section)
893      END DO
894      IF (use_multiple_para_env) &
895         CALL cp_rm_default_logger()
896      DEALLOCATE (group_distribution)
897      DEALLOCATE (i_force_eval)
898      timer_env => get_timer_env()
899      mp_perf_env => get_mp_perf_env()
900      CALL mp_max(last_f_env_id, para_env%group)
901      last_f_env_id = last_f_env_id + 1
902      new_env_id = last_f_env_id
903      n_f_envs = n_f_envs + 1
904      CALL f_env_create(f_envs(n_f_envs)%f_env, logger=logger, &
905                        timer_env=timer_env, mp_perf_env=mp_perf_env, force_env=force_env, &
906                        id_nr=last_f_env_id, old_dir=old_dir)
907      CALL force_env_release(force_env)
908      CALL globenv_release(globenv)
909      CALL section_vals_release(root_section)
910      CALL cp_para_env_release(para_env)
911      CALL f_env_rm_defaults(f_envs(n_f_envs)%f_env, ierr=ierr)
912      CALL timestop(handle)
913
914   END SUBROUTINE create_force_env
915
916! **************************************************************************************************
917!> \brief deallocates the force_env with the given id
918!> \param env_id the id of the force_env to remove
919!> \param ierr will contain a number different from 0 if
920!> \param q_finalize ...
921!> \author fawzi
922!> \note
923!>      The following routines need to be synchronized wrt. adding/removing
924!>      of the default environments (logging, performance,error):
925!>      environment:cp2k_init, environment:cp2k_finalize,
926!>      f77_interface:f_env_add_defaults, f77_interface:f_env_rm_defaults,
927!>      f77_interface:create_force_env, f77_interface:destroy_force_env
928! **************************************************************************************************
929   RECURSIVE SUBROUTINE destroy_force_env(env_id, ierr, q_finalize)
930      INTEGER, INTENT(in)                                :: env_id
931      INTEGER, INTENT(out)                               :: ierr
932      LOGICAL, INTENT(IN), OPTIONAL                      :: q_finalize
933
934      CHARACTER(len=*), PARAMETER :: routineN = 'destroy_force_env', &
935         routineP = moduleN//':'//routineN
936
937      INTEGER                                            :: env_pos, i
938      TYPE(cp_para_env_type), POINTER                    :: para_env
939      TYPE(f_env_type), POINTER                          :: f_env
940      TYPE(global_environment_type), POINTER             :: globenv
941      TYPE(section_vals_type), POINTER                   :: root_section
942
943      NULLIFY (f_env)
944      CALL f_env_add_defaults(env_id, f_env)
945      env_pos = get_pos_of_env(env_id)
946      n_f_envs = n_f_envs - 1
947      DO i = env_pos, n_f_envs
948         f_envs(i)%f_env => f_envs(i + 1)%f_env
949      END DO
950      NULLIFY (f_envs(n_f_envs + 1)%f_env)
951
952      CALL force_env_get(f_env%force_env, globenv=globenv, &
953                         root_section=root_section, para_env=para_env)
954
955      CPASSERT(ASSOCIATED(globenv))
956      CALL globenv_retain(globenv)
957      CALL f_env_dealloc(f_env)
958      IF (PRESENT(q_finalize)) THEN
959         CALL cp2k_finalize(root_section, para_env, globenv, f_env%old_path, q_finalize)
960      ELSE
961         CALL cp2k_finalize(root_section, para_env, globenv, f_env%old_path)
962      END IF
963      CALL section_vals_release(root_section)
964      CALL globenv_release(globenv)
965      DEALLOCATE (f_env)
966      ierr = 0
967   END SUBROUTINE destroy_force_env
968
969! **************************************************************************************************
970!> \brief returns the number of atoms in the given force env
971!> \param env_id id of the force_env
972!> \param n_atom ...
973!> \param ierr will return a number different from 0 if there was an error
974!> \date   22.11.2010 (MK)
975!> \author fawzi
976! **************************************************************************************************
977   SUBROUTINE get_natom(env_id, n_atom, ierr)
978
979      INTEGER, INTENT(IN)                                :: env_id
980      INTEGER, INTENT(OUT)                               :: n_atom, ierr
981
982      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_natom', routineP = moduleN//':'//routineN
983
984      TYPE(f_env_type), POINTER                          :: f_env
985
986      n_atom = 0
987      NULLIFY (f_env)
988      CALL f_env_add_defaults(env_id, f_env)
989      n_atom = force_env_get_natom(f_env%force_env)
990      CALL f_env_rm_defaults(f_env, ierr)
991
992   END SUBROUTINE get_natom
993
994! **************************************************************************************************
995!> \brief returns the number of particles in the given force env
996!> \param env_id id of the force_env
997!> \param n_particle ...
998!> \param ierr will return a number different from 0 if there was an error
999!> \author Matthias Krack
1000!>
1001! **************************************************************************************************
1002   SUBROUTINE get_nparticle(env_id, n_particle, ierr)
1003
1004      INTEGER, INTENT(IN)                                :: env_id
1005      INTEGER, INTENT(OUT)                               :: n_particle, ierr
1006
1007      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_nparticle', routineP = moduleN//':'//routineN
1008
1009      TYPE(f_env_type), POINTER                          :: f_env
1010
1011      n_particle = 0
1012      NULLIFY (f_env)
1013      CALL f_env_add_defaults(env_id, f_env)
1014      n_particle = force_env_get_nparticle(f_env%force_env)
1015      CALL f_env_rm_defaults(f_env, ierr)
1016
1017   END SUBROUTINE get_nparticle
1018
1019! **************************************************************************************************
1020!> \brief gets a cell
1021!> \param env_id id of the force_env
1022!> \param cell the array with the cell matrix
1023!> \param per periodicity
1024!> \param ierr will return a number different from 0 if there was an error
1025!> \author Joost VandeVondele
1026! **************************************************************************************************
1027   SUBROUTINE get_cell(env_id, cell, per, ierr)
1028
1029      INTEGER, INTENT(IN)                                :: env_id
1030      REAL(KIND=DP), DIMENSION(3, 3)                     :: cell
1031      INTEGER, DIMENSION(3), OPTIONAL                    :: per
1032      INTEGER, INTENT(OUT)                               :: ierr
1033
1034      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_cell', routineP = moduleN//':'//routineN
1035
1036      TYPE(cell_type), POINTER                           :: cell_full
1037      TYPE(f_env_type), POINTER                          :: f_env
1038
1039      NULLIFY (f_env)
1040      CALL f_env_add_defaults(env_id, f_env)
1041      NULLIFY (cell_full)
1042      CALL force_env_get(f_env%force_env, cell=cell_full)
1043      CPASSERT(ASSOCIATED(cell_full))
1044      cell = cell_full%hmat
1045      IF (PRESENT(per)) per(:) = cell_full%perd(:)
1046      CALL f_env_rm_defaults(f_env, ierr)
1047
1048   END SUBROUTINE get_cell
1049
1050! **************************************************************************************************
1051!> \brief gets a result from CP2K that is a real 1D array
1052!> \param env_id id of the force_env
1053!> \param description the tag of the result
1054!> \param N ...
1055!> \param RESULT ...
1056!> \param res_exist ...
1057!> \param ierr will return a number different from 0 if there was an error
1058!> \author Joost VandeVondele
1059! **************************************************************************************************
1060   SUBROUTINE get_result_r1(env_id, description, N, RESULT, res_exist, ierr)
1061      INTEGER                                            :: env_id
1062      CHARACTER(LEN=default_string_length)               :: description
1063      INTEGER                                            :: N
1064      REAL(KIND=dp), DIMENSION(1:N)                      :: RESULT
1065      LOGICAL, OPTIONAL                                  :: res_exist
1066      INTEGER                                            :: ierr
1067
1068      INTEGER                                            :: nres
1069      LOGICAL                                            :: exist_res
1070      TYPE(cp_result_type), POINTER                      :: results
1071      TYPE(cp_subsys_type), POINTER                      :: subsys
1072      TYPE(f_env_type), POINTER                          :: f_env
1073
1074      NULLIFY (f_env, subsys, results)
1075      CALL f_env_add_defaults(env_id, f_env)
1076
1077      CALL force_env_get(f_env%force_env, subsys=subsys)
1078      CALL cp_subsys_get(subsys, results=results)
1079      ! first test for the result
1080      IF (PRESENT(res_exist)) THEN
1081         res_exist = test_for_result(results, description=description)
1082         exist_res = res_exist
1083      ELSE
1084         exist_res = .TRUE.
1085      END IF
1086      ! if existing (or assuming the existance) read the results
1087      IF (exist_res) THEN
1088         CALL get_results(results, description=description, n_rep=nres)
1089         CALL get_results(results, description=description, values=RESULT, nval=nres)
1090      END IF
1091
1092      CALL f_env_rm_defaults(f_env, ierr)
1093
1094   END SUBROUTINE get_result_r1
1095
1096! **************************************************************************************************
1097!> \brief gets the forces of the particles
1098!> \param env_id id of the force_env
1099!> \param frc the array where to write the forces
1100!> \param n_el number of positions (3*nparticle) just to check
1101!> \param ierr will return a number different from 0 if there was an error
1102!> \date   22.11.2010 (MK)
1103!> \author fawzi
1104! **************************************************************************************************
1105   SUBROUTINE get_force(env_id, frc, n_el, ierr)
1106
1107      INTEGER, INTENT(IN)                                :: env_id, n_el
1108      REAL(KIND=dp), DIMENSION(1:n_el)                   :: frc
1109      INTEGER, INTENT(OUT)                               :: ierr
1110
1111      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_force', routineP = moduleN//':'//routineN
1112
1113      TYPE(f_env_type), POINTER                          :: f_env
1114
1115      NULLIFY (f_env)
1116      CALL f_env_add_defaults(env_id, f_env)
1117      CALL force_env_get_frc(f_env%force_env, frc, n_el)
1118      CALL f_env_rm_defaults(f_env, ierr)
1119
1120   END SUBROUTINE get_force
1121
1122! **************************************************************************************************
1123!> \brief gets the stress tensor
1124!> \param env_id id of the force_env
1125!> \param stress_tensor the array where to write the stress tensor
1126!> \param ierr will return a number different from 0 if there was an error
1127!> \author Ole Schuett
1128! **************************************************************************************************
1129   SUBROUTINE get_stress_tensor(env_id, stress_tensor, ierr)
1130
1131      INTEGER, INTENT(IN)                                :: env_id
1132      REAL(KIND=dp), DIMENSION(3, 3), INTENT(OUT)        :: stress_tensor
1133      INTEGER, INTENT(OUT)                               :: ierr
1134
1135      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_stress_tensor', &
1136         routineP = moduleN//':'//routineN
1137
1138      TYPE(cell_type), POINTER                           :: cell
1139      TYPE(cp_subsys_type), POINTER                      :: subsys
1140      TYPE(f_env_type), POINTER                          :: f_env
1141      TYPE(virial_type), POINTER                         :: virial
1142
1143      NULLIFY (f_env, subsys, virial, cell)
1144      stress_tensor(:, :) = 0.0_dp
1145
1146      CALL f_env_add_defaults(env_id, f_env)
1147      CALL force_env_get(f_env%force_env, subsys=subsys, cell=cell)
1148      CALL cp_subsys_get(subsys, virial=virial)
1149      IF (virial%pv_availability) THEN
1150         stress_tensor(:, :) = virial%pv_virial(:, :)/cell%deth
1151      ENDIF
1152      CALL f_env_rm_defaults(f_env, ierr)
1153
1154   END SUBROUTINE get_stress_tensor
1155
1156! **************************************************************************************************
1157!> \brief gets the positions of the particles
1158!> \param env_id id of the force_env
1159!> \param pos the array where to write the positions
1160!> \param n_el number of positions (3*nparticle) just to check
1161!> \param ierr will return a number different from 0 if there was an error
1162!> \date   22.11.2010 (MK)
1163!> \author fawzi
1164! **************************************************************************************************
1165   SUBROUTINE get_pos(env_id, pos, n_el, ierr)
1166
1167      INTEGER, INTENT(IN)                                :: env_id, n_el
1168      REAL(KIND=DP), DIMENSION(1:n_el)                   :: pos
1169      INTEGER, INTENT(OUT)                               :: ierr
1170
1171      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_pos', routineP = moduleN//':'//routineN
1172
1173      TYPE(f_env_type), POINTER                          :: f_env
1174
1175      NULLIFY (f_env)
1176      CALL f_env_add_defaults(env_id, f_env)
1177      CALL force_env_get_pos(f_env%force_env, pos, n_el)
1178      CALL f_env_rm_defaults(f_env, ierr)
1179
1180   END SUBROUTINE get_pos
1181
1182! **************************************************************************************************
1183!> \brief gets the velocities of the particles
1184!> \param env_id id of the force_env
1185!> \param vel the array where to write the velocities
1186!> \param n_el number of velocities (3*nparticle) just to check
1187!> \param ierr will return a number different from 0 if there was an error
1188!> \author fawzi
1189!> date    22.11.2010 (MK)
1190! **************************************************************************************************
1191   SUBROUTINE get_vel(env_id, vel, n_el, ierr)
1192
1193      INTEGER, INTENT(IN)                                :: env_id, n_el
1194      REAL(KIND=DP), DIMENSION(1:n_el)                   :: vel
1195      INTEGER, INTENT(OUT)                               :: ierr
1196
1197      CHARACTER(LEN=*), PARAMETER :: routineN = 'get_vel', routineP = moduleN//':'//routineN
1198
1199      TYPE(f_env_type), POINTER                          :: f_env
1200
1201      NULLIFY (f_env)
1202      CALL f_env_add_defaults(env_id, f_env)
1203      CALL force_env_get_vel(f_env%force_env, vel, n_el)
1204      CALL f_env_rm_defaults(f_env, ierr)
1205
1206   END SUBROUTINE get_vel
1207
1208! **************************************************************************************************
1209!> \brief sets a new cell
1210!> \param env_id id of the force_env
1211!> \param new_cell the array with the cell matrix
1212!> \param ierr will return a number different from 0 if there was an error
1213!> \author Joost VandeVondele
1214! **************************************************************************************************
1215   SUBROUTINE set_cell(env_id, new_cell, ierr)
1216
1217      INTEGER, INTENT(IN)                                :: env_id
1218      REAL(KIND=DP), DIMENSION(3, 3)                     :: new_cell
1219      INTEGER, INTENT(OUT)                               :: ierr
1220
1221      CHARACTER(LEN=*), PARAMETER :: routineN = 'set_cell', routineP = moduleN//':'//routineN
1222
1223      TYPE(cell_type), POINTER                           :: cell
1224      TYPE(cp_subsys_type), POINTER                      :: subsys
1225      TYPE(f_env_type), POINTER                          :: f_env
1226
1227      NULLIFY (f_env, cell, subsys)
1228      CALL f_env_add_defaults(env_id, f_env)
1229      NULLIFY (cell)
1230      CALL force_env_get(f_env%force_env, cell=cell)
1231      CPASSERT(ASSOCIATED(cell))
1232      cell%hmat = new_cell
1233      CALL init_cell(cell)
1234      CALL force_env_get(f_env%force_env, subsys=subsys)
1235      CALL cp_subsys_set(subsys, cell=cell)
1236      CALL f_env_rm_defaults(f_env, ierr)
1237
1238   END SUBROUTINE set_cell
1239
1240! **************************************************************************************************
1241!> \brief sets the positions of the particles
1242!> \param env_id id of the force_env
1243!> \param new_pos the array with the new positions
1244!> \param n_el number of positions (3*nparticle) just to check
1245!> \param ierr will return a number different from 0 if there was an error
1246!> \date   22.11.2010 updated (MK)
1247!> \author fawzi
1248! **************************************************************************************************
1249   SUBROUTINE set_pos(env_id, new_pos, n_el, ierr)
1250
1251      INTEGER, INTENT(IN)                                :: env_id, n_el
1252      REAL(KIND=dp), DIMENSION(1:n_el)                   :: new_pos
1253      INTEGER, INTENT(OUT)                               :: ierr
1254
1255      CHARACTER(LEN=*), PARAMETER :: routineN = 'set_pos', routineP = moduleN//':'//routineN
1256
1257      TYPE(cp_subsys_type), POINTER                      :: subsys
1258      TYPE(f_env_type), POINTER                          :: f_env
1259
1260      NULLIFY (f_env)
1261      CALL f_env_add_defaults(env_id, f_env)
1262      NULLIFY (subsys)
1263      CALL force_env_get(f_env%force_env, subsys=subsys)
1264      CALL unpack_subsys_particles(subsys=subsys, r=new_pos)
1265      CALL f_env_rm_defaults(f_env, ierr)
1266
1267   END SUBROUTINE set_pos
1268
1269! **************************************************************************************************
1270!> \brief sets the velocities of the particles
1271!> \param env_id id of the force_env
1272!> \param new_vel the array with the new velocities
1273!> \param n_el number of velocities (3*nparticle) just to check
1274!> \param ierr will return a number different from 0 if there was an error
1275!> \date   22.11.2010 updated (MK)
1276!> \author fawzi
1277! **************************************************************************************************
1278   SUBROUTINE set_vel(env_id, new_vel, n_el, ierr)
1279
1280      INTEGER, INTENT(IN)                                :: env_id, n_el
1281      REAL(kind=dp), DIMENSION(1:n_el)                   :: new_vel
1282      INTEGER, INTENT(OUT)                               :: ierr
1283
1284      CHARACTER(LEN=*), PARAMETER :: routineN = 'set_vel', routineP = moduleN//':'//routineN
1285
1286      TYPE(cp_subsys_type), POINTER                      :: subsys
1287      TYPE(f_env_type), POINTER                          :: f_env
1288
1289      NULLIFY (f_env)
1290      CALL f_env_add_defaults(env_id, f_env)
1291      NULLIFY (subsys)
1292      CALL force_env_get(f_env%force_env, subsys=subsys)
1293      CALL unpack_subsys_particles(subsys=subsys, v=new_vel)
1294      CALL f_env_rm_defaults(f_env, ierr)
1295
1296   END SUBROUTINE set_vel
1297
1298! **************************************************************************************************
1299!> \brief updates the energy and the forces of given force_env
1300!> \param env_id id of the force_env that you want to update
1301!> \param calc_force if the forces should be updated, if false the forces
1302!>        might be wrong.
1303!> \param ierr will return a number different from 0 if there was an error
1304!> \author fawzi
1305! **************************************************************************************************
1306   RECURSIVE SUBROUTINE calc_energy_force(env_id, calc_force, ierr)
1307
1308      INTEGER, INTENT(in)                                :: env_id
1309      LOGICAL, INTENT(in)                                :: calc_force
1310      INTEGER, INTENT(out)                               :: ierr
1311
1312      CHARACTER(len=*), PARAMETER :: routineN = 'calc_energy_force', &
1313         routineP = moduleN//':'//routineN
1314
1315      TYPE(cp_logger_type), POINTER                      :: logger
1316      TYPE(f_env_type), POINTER                          :: f_env
1317
1318      NULLIFY (f_env)
1319      CALL f_env_add_defaults(env_id, f_env)
1320      logger => cp_get_default_logger()
1321      CALL cp_iterate(logger%iter_info) ! add one to the iteration count
1322      CALL force_env_calc_energy_force(f_env%force_env, calc_force=calc_force)
1323      CALL f_env_rm_defaults(f_env, ierr)
1324
1325   END SUBROUTINE calc_energy_force
1326
1327! **************************************************************************************************
1328!> \brief returns the energy of the last configuration calculated
1329!> \param env_id id of the force_env that you want to update
1330!> \param e_pot the potential energy of the system
1331!> \param ierr will return a number different from 0 if there was an error
1332!> \author fawzi
1333! **************************************************************************************************
1334   SUBROUTINE get_energy(env_id, e_pot, ierr)
1335
1336      INTEGER, INTENT(in)                                :: env_id
1337      REAL(kind=dp), INTENT(out)                         :: e_pot
1338      INTEGER, INTENT(out)                               :: ierr
1339
1340      CHARACTER(len=*), PARAMETER :: routineN = 'get_energy', routineP = moduleN//':'//routineN
1341
1342      TYPE(f_env_type), POINTER                          :: f_env
1343
1344      NULLIFY (f_env)
1345      CALL f_env_add_defaults(env_id, f_env)
1346      CALL force_env_get(f_env%force_env, potential_energy=e_pot)
1347      CALL f_env_rm_defaults(f_env, ierr)
1348
1349   END SUBROUTINE get_energy
1350
1351! **************************************************************************************************
1352!> \brief returns the energy of the configuration given by the positions
1353!>      passed as argument
1354!> \param env_id id of the force_env that you want to update
1355!> \param pos array with the positions
1356!> \param n_el number of elements in pos (3*natom)
1357!> \param e_pot the potential energy of the system
1358!> \param ierr will return a number different from 0 if there was an error
1359!> \author fawzi
1360!> \note
1361!>      utility call
1362! **************************************************************************************************
1363   RECURSIVE SUBROUTINE calc_energy(env_id, pos, n_el, e_pot, ierr)
1364
1365      INTEGER, INTENT(IN)                                :: env_id, n_el
1366      REAL(KIND=dp), DIMENSION(1:n_el), INTENT(IN)       :: pos
1367      REAL(KIND=dp), INTENT(OUT)                         :: e_pot
1368      INTEGER, INTENT(OUT)                               :: ierr
1369
1370      REAL(KIND=dp), DIMENSION(1)                        :: dummy_f
1371
1372      CALL calc_force(env_id, pos, n_el, e_pot, dummy_f, 0, ierr)
1373
1374   END SUBROUTINE calc_energy
1375
1376! **************************************************************************************************
1377!> \brief returns the energy of the configuration given by the positions
1378!>      passed as argument
1379!> \param env_id id of the force_env that you want to update
1380!> \param pos array with the positions
1381!> \param n_el_pos number of elements in pos (3*natom)
1382!> \param e_pot the potential energy of the system
1383!> \param force array that will contain the forces
1384!> \param n_el_force number of elements in force (3*natom). If 0 the
1385!>        forces are not calculated
1386!> \param ierr will return a number different from 0 if there was an error
1387!> \author fawzi
1388!> \note
1389!>      utility call, but actually it could be a better and more efficient
1390!>      interface to connect to other codes if cp2k would be deeply
1391!>      refactored
1392! **************************************************************************************************
1393   RECURSIVE SUBROUTINE calc_force(env_id, pos, n_el_pos, e_pot, force, n_el_force, ierr)
1394
1395      INTEGER, INTENT(in)                                :: env_id, n_el_pos
1396      REAL(kind=dp), DIMENSION(1:n_el_pos), INTENT(in)   :: pos
1397      REAL(kind=dp), INTENT(out)                         :: e_pot
1398      INTEGER, INTENT(in)                                :: n_el_force
1399      REAL(kind=dp), DIMENSION(1:n_el_force), &
1400         INTENT(inout)                                   :: force
1401      INTEGER, INTENT(out)                               :: ierr
1402
1403      LOGICAL                                            :: calc_f
1404
1405      calc_f = (n_el_force /= 0)
1406      CALL set_pos(env_id, pos, n_el_pos, ierr)
1407      IF (ierr == 0) CALL calc_energy_force(env_id, calc_f, ierr)
1408      IF (ierr == 0) CALL get_energy(env_id, e_pot, ierr)
1409      IF (calc_f .AND. (ierr == 0)) CALL get_force(env_id, force, n_el_force, ierr)
1410
1411   END SUBROUTINE calc_force
1412
1413! **************************************************************************************************
1414!> \brief performs a check of the input
1415!> \param input_declaration ...
1416!> \param input_file_path the path of the input file to check
1417!> \param output_file_path path of the output file (to which it is appended)
1418!>        if it is "__STD_OUT__" the default_output_unit is used
1419!> \param echo_input if the parsed input should be written out with all the
1420!>        defaults made explicit
1421!> \param mpi_comm the mpi communicator (if not given it uses the default
1422!>        one)
1423!> \param ierr error control, if different from 0 there was an error
1424!> \author fawzi
1425! **************************************************************************************************
1426   SUBROUTINE check_input(input_declaration, input_file_path, output_file_path, &
1427                          echo_input, mpi_comm, ierr)
1428      TYPE(section_type), POINTER                        :: input_declaration
1429      CHARACTER(len=*), INTENT(in)                       :: input_file_path, output_file_path
1430      LOGICAL, INTENT(in), OPTIONAL                      :: echo_input
1431      INTEGER, INTENT(in), OPTIONAL                      :: mpi_comm
1432      INTEGER, INTENT(out)                               :: ierr
1433
1434      CHARACTER(len=*), PARAMETER :: routineN = 'check_input', routineP = moduleN//':'//routineN
1435
1436      INTEGER                                            :: unit_nr
1437      LOGICAL                                            :: my_echo_input
1438      TYPE(cp_logger_type), POINTER                      :: logger
1439      TYPE(cp_para_env_type), POINTER                    :: para_env
1440      TYPE(section_vals_type), POINTER                   :: input_file
1441
1442      my_echo_input = .FALSE.
1443      IF (PRESENT(echo_input)) my_echo_input = echo_input
1444
1445      IF (PRESENT(mpi_comm)) THEN
1446         NULLIFY (para_env)
1447         CALL cp_para_env_create(para_env, group=mpi_comm)
1448      ELSE
1449         para_env => default_para_env
1450         CALL cp_para_env_retain(para_env)
1451      END IF
1452      IF (para_env%ionode) THEN
1453         IF (output_file_path == "__STD_OUT__") THEN
1454            unit_nr = default_output_unit
1455         ELSE
1456            CALL open_file(file_name=output_file_path, file_status="UNKNOWN", &
1457                           file_action="WRITE", file_position="APPEND", &
1458                           unit_number=unit_nr)
1459         END IF
1460      ELSE
1461         unit_nr = -1
1462      END IF
1463
1464      NULLIFY (logger)
1465      CALL cp_logger_create(logger, para_env=para_env, &
1466                            default_global_unit_nr=unit_nr, &
1467                            close_global_unit_on_dealloc=.FALSE.)
1468      CALL cp_add_default_logger(logger)
1469      CALL cp_logger_release(logger)
1470
1471      input_file => read_input(input_declaration, input_file_path, initial_variables=empty_initial_variables, &
1472                               para_env=para_env)
1473      CALL check_cp2k_input(input_declaration, input_file, para_env=para_env, output_unit=unit_nr)
1474      IF (my_echo_input .AND. para_env%ionode) THEN
1475         CALL section_vals_write(input_file, &
1476                                 unit_nr=cp_logger_get_default_unit_nr(logger, local=.FALSE.), hide_root=.TRUE., &
1477                                 hide_defaults=.FALSE.)
1478      END IF
1479      CALL section_vals_release(input_file)
1480
1481      CALL cp_logger_release(logger)
1482      CALL cp_para_env_release(para_env)
1483      ierr = 0
1484      CALL cp_rm_default_logger()
1485   END SUBROUTINE check_input
1486
1487END MODULE f77_interface
1488