LCOV - code coverage report
Current view: top level - src/motion - gopt_f77_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:5911e77) Lines: 137 146 93.8 %
Date: 2024-12-12 07:20:04 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief evaluete the potential energy and its gradients using an array
      10             : !>      with same dimension as the particle_set
      11             : !> \param gopt_env the geometry optimization environment
      12             : !> \param x the position where the function should be evaluated
      13             : !> \param f the function value
      14             : !> \param gradient the value of its gradient
      15             : !> \param master ...
      16             : !> \param final_evaluation ...
      17             : !> \param para_env ...
      18             : !> \par History
      19             : !>       CELL OPTIMIZATION:  Teodoro Laino [tlaino] - University of Zurich - 03.2008
      20             : !>       07.2020 Pierre Cazade [pcazade] Space Group Symmetry
      21             : !> \author Teodoro Laino [tlaino] - University of Zurich - 01.2008
      22             : ! **************************************************************************************************
      23       20542 : SUBROUTINE cp_eval_at(gopt_env, x, f, gradient, master, &
      24             :                       final_evaluation, para_env)
      25             : 
      26             :    USE cp_log_handling, ONLY: cp_logger_type
      27             :    USE averages_types, ONLY: average_quantities_type, &
      28             :                              create_averages, &
      29             :                              release_averages
      30             :    USE bibliography, ONLY: Henkelman1999, &
      31             :                            cite_reference
      32             :    USE cell_opt_utils, ONLY: get_dg_dh, &
      33             :                              gopt_new_logger_create, &
      34             :                              gopt_new_logger_release
      35             :    USE cell_types, ONLY: cell_type
      36             :    USE cell_methods, ONLY: write_cell
      37             :    USE message_passing, ONLY: mp_para_env_type
      38             :    USE cp_subsys_types, ONLY: cp_subsys_get, &
      39             :                               cp_subsys_type, &
      40             :                               pack_subsys_particles, &
      41             :                               unpack_subsys_particles
      42             :    USE dimer_methods, ONLY: cp_eval_at_ts
      43             :    USE force_env_methods, ONLY: force_env_calc_energy_force
      44             :    USE force_env_types, ONLY: force_env_get, &
      45             :                               force_env_get_nparticle
      46             :    USE geo_opt, ONLY: cp_geo_opt
      47             :    USE gopt_f_types, ONLY: gopt_f_type
      48             :    USE gopt_f_methods, ONLY: apply_cell_change
      49             :    USE input_constants, ONLY: default_minimization_method_id, &
      50             :                               default_ts_method_id, &
      51             :                               default_cell_direct_id, &
      52             :                               default_cell_method_id, &
      53             :                               default_cell_geo_opt_id, &
      54             :                               default_cell_md_id, &
      55             :                               default_shellcore_method_id, &
      56             :                               nvt_ensemble, &
      57             :                               mol_dyn_run, &
      58             :                               geo_opt_run, &
      59             :                               cell_opt_run, &
      60             :                               fix_none
      61             :    USE input_section_types, ONLY: section_vals_get, &
      62             :                                   section_vals_get_subs_vals, &
      63             :                                   section_vals_type, &
      64             :                                   section_vals_val_get
      65             :    USE md_run, ONLY: qs_mol_dyn
      66             :    USE kinds, ONLY: dp, &
      67             :                     default_string_length
      68             :    USE particle_list_types, ONLY: particle_list_type
      69             :    USE particle_methods, ONLY: write_structure_data
      70             :    USE virial_methods, ONLY: virial_update
      71             :    USE virial_types, ONLY: virial_type
      72             :    USE cp_log_handling, ONLY: cp_add_default_logger, &
      73             :                               cp_rm_default_logger
      74             :    USE space_groups_types, ONLY: spgr_type
      75             :    USE space_groups, ONLY: spgr_apply_rotations_stress, &
      76             :                            spgr_apply_rotations_coord, &
      77             :                            spgr_apply_rotations_force, &
      78             :                            spgr_write_stress_tensor
      79             : 
      80             : #include "../base/base_uses.f90"
      81             :    IMPLICIT NONE
      82             :    TYPE(gopt_f_type), POINTER               :: gopt_env
      83             :    REAL(KIND=dp), DIMENSION(:), POINTER     :: x
      84             :    REAL(KIND=dp), INTENT(OUT), OPTIONAL     :: f
      85             :    REAL(KIND=dp), DIMENSION(:), OPTIONAL, &
      86             :       POINTER                                :: gradient
      87             :    INTEGER, INTENT(IN)                      :: master
      88             :    LOGICAL, INTENT(IN), OPTIONAL            :: final_evaluation
      89             :    TYPE(mp_para_env_type), POINTER          :: para_env
      90             : 
      91             :    CHARACTER(len=*), PARAMETER :: routineN = 'cp_eval_at'
      92             : 
      93             :    INTEGER                                  :: ensemble, handle, idg, idir, ip, &
      94             :                                                nparticle, nsize, shell_index
      95             :    LOGICAL                                  :: explicit, my_final_evaluation
      96             :    REAL(KIND=dp)                            :: f_ts, potential_energy
      97             :    REAL(KIND=dp), DIMENSION(3, 3)           :: av_ptens
      98       20542 :    REAL(KIND=dp), DIMENSION(:), POINTER     :: cell_gradient, gradient_ts
      99             :    TYPE(cell_type), POINTER                 :: cell
     100             :    TYPE(cp_subsys_type), POINTER            :: subsys
     101             :    TYPE(particle_list_type), POINTER        :: core_particles, particles, &
     102             :                                                shell_particles
     103             :    TYPE(virial_type), POINTER               :: virial
     104             :    TYPE(cp_logger_type), POINTER            :: new_logger
     105             :    CHARACTER(LEN=default_string_length)     :: project_name
     106             :    TYPE(average_quantities_type), POINTER   :: averages
     107             :    TYPE(section_vals_type), POINTER         :: work, avgs_section
     108             :    TYPE(spgr_type), POINTER                 :: spgr
     109             : 
     110       20542 :    NULLIFY (averages)
     111       20542 :    NULLIFY (cell)
     112       20542 :    NULLIFY (core_particles)
     113       20542 :    NULLIFY (gradient_ts)
     114       20542 :    NULLIFY (particles)
     115       20542 :    NULLIFY (shell_particles)
     116       20542 :    NULLIFY (subsys)
     117       20542 :    NULLIFY (virial)
     118       20542 :    NULLIFY (new_logger)
     119       20542 :    NULLIFY (spgr)
     120             : 
     121       20542 :    CALL timeset(routineN, handle)
     122             : 
     123       20542 :    CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
     124             :    CALL cp_subsys_get(subsys, &
     125             :                       core_particles=core_particles, &
     126             :                       particles=particles, &
     127             :                       shell_particles=shell_particles, &
     128       20542 :                       virial=virial)
     129             : 
     130       20542 :    spgr => gopt_env%spgr
     131             : 
     132       20542 :    my_final_evaluation = .FALSE.
     133       20542 :    IF (PRESENT(final_evaluation)) my_final_evaluation = final_evaluation
     134             : 
     135       34864 :    SELECT CASE (gopt_env%type_id)
     136             :    CASE (default_minimization_method_id, default_ts_method_id)
     137       14322 :       CALL unpack_subsys_particles(subsys=subsys, r=x)
     138       14322 :       CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
     139       14322 :       SELECT CASE (gopt_env%type_id)
     140             :       CASE (default_minimization_method_id)
     141             :          ! Geometry Minimization
     142             :          CALL force_env_calc_energy_force(gopt_env%force_env, &
     143             :                                           calc_force=PRESENT(gradient), &
     144       12506 :                                           require_consistent_energy_force=gopt_env%require_consistent_energy_force)
     145             :          ! Possibly take the potential energy
     146       12506 :          IF (PRESENT(f)) THEN
     147       12506 :             CALL force_env_get(gopt_env%force_env, potential_energy=f)
     148             :          END IF
     149             :          ! Possibly take the gradients
     150       12506 :          IF (PRESENT(gradient)) THEN
     151       11965 :             IF (master == para_env%mepos) THEN ! we are on the master
     152       10932 :                CALL pack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
     153       10932 :                IF (spgr%keep_space_group) THEN
     154           0 :                   CALL spgr_apply_rotations_force(spgr, gradient)
     155           0 :                   CALL unpack_subsys_particles(subsys=subsys, f=gradient, fscale=-1.0_dp)
     156             :                END IF
     157             :             END IF
     158             :          END IF
     159             :       CASE (default_ts_method_id)
     160             :          ! Transition State Optimization
     161        5448 :          ALLOCATE (gradient_ts(particles%n_els*3))
     162             :          ! Real calculation of energy and forces for transition state optimization:
     163             :          ! When doing dimer methods forces have to be always computed since the function
     164             :          ! to minimize is not the energy but the effective force
     165        1816 :          CALL cp_eval_at_ts(gopt_env, x, f_ts, gradient_ts, calc_force=.TRUE.)
     166        1816 :          CALL cite_reference(Henkelman1999)
     167             :          ! Possibly take the potential energy
     168        1816 :          IF (PRESENT(f)) f = f_ts
     169             :          ! Possibly take the gradients
     170        1816 :          IF (PRESENT(gradient)) THEN
     171         854 :             IF (master == para_env%mepos) THEN ! we are on the master
     172         854 :                CPASSERT(ASSOCIATED(gradient))
     173       29452 :                gradient = gradient_ts
     174             :             END IF
     175             :          END IF
     176        1816 :          DEALLOCATE (gradient_ts)
     177             :       END SELECT
     178             :       ! This call is necessary for QM/MM if a Translation is applied
     179             :       ! this makes the geometry optimizer consistent
     180       14322 :       CALL unpack_subsys_particles(subsys=subsys, r=x)
     181             :    CASE (default_cell_method_id)
     182             :       ! Check for VIRIAL
     183        5880 :       IF (.NOT. virial%pv_availability) &
     184             :          CALL cp_abort(__LOCATION__, &
     185             :                        "Cell optimization requested but FORCE_EVAL%STRESS_TENSOR was not defined! "// &
     186           0 :                        "Activate the evaluation of the stress tensor for cell optimization!")
     187       11160 :       SELECT CASE (gopt_env%cell_method_id)
     188             :       CASE (default_cell_direct_id)
     189        4940 :          CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
     190             :          ! Possibly output the new cell used for the next calculation
     191        4940 :          CALL write_cell(cell, gopt_env%geo_section)
     192             :          ! Compute the pressure tensor
     193        4940 :          BLOCK
     194             :             TYPE(virial_type) :: virial_avg
     195             :             CALL force_env_calc_energy_force(gopt_env%force_env, &
     196             :                                              calc_force=PRESENT(gradient), &
     197        4940 :                                              require_consistent_energy_force=gopt_env%require_consistent_energy_force)
     198             :             ! Possibly take the potential energy
     199        4940 :             CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
     200        4940 :             virial_avg = virial
     201        4940 :             CALL virial_update(virial_avg, subsys, para_env)
     202        4940 :             IF (PRESENT(f)) THEN
     203        4940 :                CALL force_env_get(gopt_env%force_env, potential_energy=f)
     204             :             END IF
     205             :             ! Possibly take the gradients
     206     1225120 :             IF (PRESENT(gradient)) THEN
     207        4526 :                CPASSERT(ANY(virial_avg%pv_total /= 0))
     208             :                ! Convert the average ptens
     209       58838 :                av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
     210        4526 :                IF (master == para_env%mepos) THEN ! we are on the master
     211        3014 :                   CPASSERT(ASSOCIATED(gradient))
     212        3014 :                   nparticle = force_env_get_nparticle(gopt_env%force_env)
     213        3014 :                   nsize = 3*nparticle
     214        3014 :                   CPASSERT((SIZE(gradient) == nsize + 6))
     215        3014 :                   CALL pack_subsys_particles(subsys=subsys, f=gradient(1:nsize), fscale=-1.0_dp)
     216        3014 :                   CALL apply_cell_change(gopt_env, cell, gradient, update_forces=.TRUE.)
     217        3014 :                   IF (spgr%keep_space_group) THEN
     218          20 :                      CALL spgr_apply_rotations_force(spgr, gradient)
     219          20 :                      CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
     220          20 :                      CALL spgr_write_stress_tensor(av_ptens, spgr)
     221             :                   END IF
     222        3014 :                   cell_gradient => gradient(nsize + 1:nsize + 6)
     223       21098 :                   cell_gradient = 0.0_dp
     224             :                   CALL get_dg_dh(cell_gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
     225             :                                  keep_angles=gopt_env%cell_env%keep_angles, &
     226             :                                  keep_symmetry=gopt_env%cell_env%keep_symmetry, &
     227             :                                  pres_int=gopt_env%cell_env%pres_int, &
     228             :                                  pres_constr=gopt_env%cell_env%pres_constr, &
     229        3014 :                                  constraint_id=gopt_env%cell_env%constraint_id)
     230             :                END IF
     231             :                ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
     232             :                ! Assume at least master==0
     233        4526 :                CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
     234        4526 :                IF (gopt_env%cell_env%constraint_id /= fix_none) &
     235          26 :                   CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
     236             :             END IF
     237             :          END BLOCK
     238             :       CASE (default_cell_geo_opt_id, default_cell_md_id)
     239         940 :          CALL apply_cell_change(gopt_env, cell, x, update_forces=.FALSE.)
     240             :          ! Possibly output the new cell used for the next calculation
     241         940 :          CALL write_cell(cell, gopt_env%geo_section)
     242             :          ! Compute the pressure tensor
     243         940 :          BLOCK
     244             :             TYPE(virial_type) :: virial_avg
     245         940 :             IF (my_final_evaluation) THEN
     246             :                CALL force_env_calc_energy_force(gopt_env%force_env, &
     247             :                                                 calc_force=PRESENT(gradient), &
     248          24 :                                                 require_consistent_energy_force=gopt_env%require_consistent_energy_force)
     249          24 :                IF (PRESENT(f)) THEN
     250          24 :                   CALL force_env_get(gopt_env%force_env, potential_energy=f)
     251             :                END IF
     252             :             ELSE
     253        1812 :                SELECT CASE (gopt_env%cell_method_id)
     254             :                CASE (default_cell_geo_opt_id)
     255         896 :                   work => section_vals_get_subs_vals(gopt_env%motion_section, "GEO_OPT")
     256         896 :                   CALL section_vals_get(work, explicit=explicit)
     257         896 :                   IF (.NOT. explicit) &
     258             :                      CALL cp_abort(__LOCATION__, &
     259           0 :                                    "Cell optimization at 0K was requested. GEO_OPT section MUST be provided in the input file!")
     260             :                   ! Perform a geometry optimization
     261             :                   CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
     262         896 :                                               project_name, id_run=geo_opt_run)
     263         896 :                   CALL cp_add_default_logger(new_logger)
     264         896 :                   CALL cp_geo_opt(gopt_env%force_env, gopt_env%globenv, eval_opt_geo=.FALSE.)
     265         896 :                   CALL force_env_get(gopt_env%force_env, potential_energy=potential_energy)
     266         896 :                   virial_avg = virial
     267             :                CASE (default_cell_md_id)
     268          20 :                   work => section_vals_get_subs_vals(gopt_env%motion_section, "MD")
     269          20 :                   avgs_section => section_vals_get_subs_vals(work, "AVERAGES")
     270          20 :                   CALL section_vals_get(work, explicit=explicit)
     271          20 :                   IF (.NOT. explicit) &
     272             :                      CALL cp_abort( &
     273             :                      __LOCATION__, &
     274           0 :                      "Cell optimization at finite temperature was requested. MD section MUST be provided in the input file!")
     275             :                   ! Only NVT ensemble is allowed..
     276          20 :                   CALL section_vals_val_get(gopt_env%motion_section, "MD%ENSEMBLE", i_val=ensemble)
     277          20 :                   IF (ensemble /= nvt_ensemble) &
     278             :                      CALL cp_abort(__LOCATION__, &
     279           0 :                                    "Cell optimization at finite temperature requires the NVT MD ensemble!")
     280             :                   ! Perform a molecular dynamics
     281             :                   CALL gopt_new_logger_create(new_logger, gopt_env%force_env%root_section, para_env, &
     282          20 :                                               project_name, id_run=mol_dyn_run)
     283          20 :                   CALL cp_add_default_logger(new_logger)
     284          20 :                   CALL create_averages(averages, avgs_section, virial_avg=.TRUE., force_env=gopt_env%force_env)
     285          20 :                   CALL qs_mol_dyn(gopt_env%force_env, gopt_env%globenv, averages, rm_restart_info=.FALSE.)
     286             :                   ! Retrieve the average of the stress tensor and the average of the potential energy
     287          20 :                   potential_energy = averages%avepot
     288          20 :                   virial_avg = averages%virial
     289          20 :                   CALL release_averages(averages)
     290             :                CASE DEFAULT
     291        2748 :                   CPABORT("")
     292             :                END SELECT
     293         916 :                CALL cp_rm_default_logger()
     294             :                CALL gopt_new_logger_release(new_logger, gopt_env%force_env%root_section, para_env, project_name, &
     295         916 :                                             cell_opt_run)
     296             :                ! Update the virial
     297         916 :                CALL virial_update(virial_avg, subsys, para_env)
     298             :                ! Possibly take give back the potential energy
     299         916 :                IF (PRESENT(f)) THEN
     300         916 :                   f = potential_energy
     301             :                END IF
     302             :             END IF
     303             :             ! Possibly give back the gradients
     304      233120 :             IF (PRESENT(gradient)) THEN
     305         916 :                CPASSERT(ANY(virial_avg%pv_total /= 0))
     306             :                ! Convert the average ptens
     307       11908 :                av_ptens(:, :) = virial_avg%pv_total(:, :)/cell%deth
     308         916 :                IF (master == para_env%mepos) THEN ! we are on the master
     309         885 :                   CPASSERT(ASSOCIATED(gradient))
     310         885 :                   IF (spgr%keep_space_group) THEN
     311           0 :                      CALL spgr_apply_rotations_stress(spgr, cell, av_ptens)
     312           0 :                      CALL spgr_write_stress_tensor(av_ptens, spgr)
     313             :                   END IF
     314             :                   ! Compute the gradients on the cell
     315             :                   CALL get_dg_dh(gradient, av_ptens, gopt_env%cell_env%pres_ext, cell, gopt_env%cell_env%mtrx, &
     316             :                                  keep_angles=gopt_env%cell_env%keep_angles, &
     317             :                                  keep_symmetry=gopt_env%cell_env%keep_symmetry, &
     318             :                                  pres_int=gopt_env%cell_env%pres_int, &
     319             :                                  pres_constr=gopt_env%cell_env%pres_constr, &
     320         885 :                                  constraint_id=gopt_env%cell_env%constraint_id)
     321             :                END IF
     322             :                ! some callers expect pres_int to be available on all ranks. Also, here master is not necessarily a single rank.
     323             :                ! Assume at least master==0
     324         916 :                CALL para_env%bcast(gopt_env%cell_env%pres_int, 0)
     325         916 :                IF (gopt_env%cell_env%constraint_id /= fix_none) &
     326           0 :                   CALL para_env%bcast(gopt_env%cell_env%pres_constr, 0)
     327             :             END IF
     328             :          END BLOCK
     329             :       CASE DEFAULT
     330        5880 :          CPABORT("")
     331             :       END SELECT
     332             :    CASE (default_shellcore_method_id)
     333         340 :       idg = 0
     334       32980 :       DO ip = 1, particles%n_els
     335       32640 :          shell_index = particles%els(ip)%shell_index
     336       32980 :          IF (shell_index /= 0) THEN
     337      130560 :             DO idir = 1, 3
     338       97920 :                idg = 3*(shell_index - 1) + idir
     339      130560 :                shell_particles%els(shell_index)%r(idir) = core_particles%els(ip)%r(idir) - x(idg)
     340             :             END DO
     341             :          END IF
     342             :       END DO
     343         340 :       CALL write_structure_data(particles%els, cell, gopt_env%motion_section)
     344             : 
     345             :       ! Shell-core optimization
     346             :       CALL force_env_calc_energy_force(gopt_env%force_env, &
     347             :                                        calc_force=PRESENT(gradient), &
     348         340 :                                        require_consistent_energy_force=gopt_env%require_consistent_energy_force)
     349             : 
     350             :       ! Possibly take the potential energy
     351         340 :       IF (PRESENT(f)) THEN
     352         340 :          CALL force_env_get(gopt_env%force_env, potential_energy=f)
     353             :       END IF
     354             : 
     355             :       ! Possibly take the gradients
     356         340 :       IF (PRESENT(gradient)) THEN
     357         340 :          IF (master == para_env%mepos) THEN ! we are on the master
     358         340 :             CPASSERT(ASSOCIATED(gradient))
     359         340 :             idg = 0
     360       32980 :             DO ip = 1, shell_particles%n_els
     361      130900 :                DO idir = 1, 3
     362       97920 :                   idg = idg + 1
     363      130560 :                   gradient(idg) = -(core_particles%els(ip)%f(idir) - shell_particles%els(ip)%f(idir))
     364             :                END DO
     365             :             END DO
     366             :          END IF
     367             :       END IF
     368             :    CASE DEFAULT
     369       20542 :       CPABORT("")
     370             :    END SELECT
     371             : 
     372       20542 :    CALL timestop(handle)
     373             : 
     374       41084 : END SUBROUTINE cp_eval_at

Generated by: LCOV version 1.15