LCOV - code coverage report
Current view: top level - src/motion - gopt_f_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:b195825) Lines: 383 386 99.2 %
Date: 2024-04-20 06:29:22 Functions: 15 15 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 contains a functional that calculates the energy and its derivatives
      10             : !>      for the geometry optimizer
      11             : !> \par History
      12             : !>      none
      13             : ! **************************************************************************************************
      14             : MODULE gopt_f_methods
      15             :    USE atomic_kind_list_types, ONLY: atomic_kind_list_type
      16             :    USE atomic_kind_types, ONLY: atomic_kind_type, &
      17             :                                 get_atomic_kind_set
      18             :    USE cell_methods, ONLY: cell_create, &
      19             :                            init_cell
      20             :    USE cell_types, ONLY: cell_copy, &
      21             :                          cell_release, &
      22             :                          cell_type, &
      23             :                          real_to_scaled, &
      24             :                          scaled_to_real
      25             :    USE cp_log_handling, ONLY: cp_to_string
      26             :    USE cp_subsys_types, ONLY: cp_subsys_get, &
      27             :                               cp_subsys_set, &
      28             :                               cp_subsys_type, &
      29             :                               pack_subsys_particles
      30             :    USE cp_units, ONLY: cp_unit_from_cp2k
      31             :    USE dimer_types, ONLY: dimer_env_type
      32             :    USE dimer_utils, ONLY: update_dimer_vec
      33             :    USE distribution_1d_types, ONLY: distribution_1d_type
      34             :    USE force_env_types, ONLY: force_env_get, &
      35             :                               force_env_get_natom, &
      36             :                               force_env_get_nparticle, &
      37             :                               force_env_type, &
      38             :                               use_qmmm, &
      39             :                               use_qmmmx
      40             :    USE gopt_f_types, ONLY: gopt_f_type
      41             :    USE gopt_param_types, ONLY: gopt_param_type
      42             :    USE input_constants, ONLY: default_cell_direct_id, &
      43             :                               default_cell_geo_opt_id, &
      44             :                               default_cell_md_id, &
      45             :                               default_cell_method_id, &
      46             :                               default_minimization_method_id, &
      47             :                               default_shellcore_method_id, &
      48             :                               default_ts_method_id, &
      49             :                               fix_none
      50             :    USE input_cp2k_restarts, ONLY: write_restart
      51             :    USE input_section_types, ONLY: section_vals_type, &
      52             :                                   section_vals_val_get
      53             :    USE kinds, ONLY: default_string_length, &
      54             :                     dp, &
      55             :                     int_8
      56             :    USE machine, ONLY: m_flush
      57             :    USE md_energies, ONLY: sample_memory
      58             :    USE message_passing, ONLY: mp_para_env_type
      59             :    USE motion_utils, ONLY: write_simulation_cell, &
      60             :                            write_stress_tensor, &
      61             :                            write_trajectory
      62             :    USE particle_list_types, ONLY: particle_list_type
      63             :    USE particle_methods, ONLY: write_structure_data
      64             :    USE particle_types, ONLY: particle_type
      65             :    USE qmmm_util, ONLY: apply_qmmm_translate
      66             :    USE qmmmx_util, ONLY: apply_qmmmx_translate
      67             :    USE virial_methods, ONLY: virial_evaluate
      68             :    USE virial_types, ONLY: virial_type
      69             : #include "../base/base_uses.f90"
      70             : 
      71             :    IMPLICIT NONE
      72             :    PRIVATE
      73             : 
      74             :    #:include "gopt_f77_methods.fypp"
      75             : 
      76             :    LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
      77             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gopt_f_methods'
      78             : 
      79             :    PUBLIC :: gopt_f_create_x0, &
      80             :              print_geo_opt_header, print_geo_opt_nc, &
      81             :              gopt_f_io_init, gopt_f_io, gopt_f_io_finalize, gopt_f_ii, &
      82             :              apply_cell_change
      83             : 
      84             : CONTAINS
      85             : 
      86             : ! **************************************************************************************************
      87             : !> \brief returns the value of the parameters for the actual configuration
      88             : !> \param gopt_env the geometry optimization environment you want the info about
      89             : !>      x0: the parameter vector (is allocated by this routine)
      90             : !> \param x0 ...
      91             : !> \par History
      92             : !>      - Cell optimization revised (06.11.2012,MK)
      93             : ! **************************************************************************************************
      94        1970 :    SUBROUTINE gopt_f_create_x0(gopt_env, x0)
      95             : 
      96             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
      97             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
      98             : 
      99             :       INTEGER                                            :: i, idg, j, nparticle
     100             :       TYPE(cell_type), POINTER                           :: cell
     101             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     102             : 
     103        1970 :       NULLIFY (cell)
     104        1970 :       NULLIFY (subsys)
     105             : 
     106        3730 :       SELECT CASE (gopt_env%type_id)
     107             :       CASE (default_minimization_method_id, default_ts_method_id)
     108        1760 :          CALL force_env_get(gopt_env%force_env, subsys=subsys)
     109             :          ! before starting we handle the case of translating coordinates (QM/MM)
     110        1760 :          IF (gopt_env%force_env%in_use == use_qmmm) &
     111          36 :             CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
     112        1760 :          IF (gopt_env%force_env%in_use == use_qmmmx) &
     113           0 :             CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
     114        1760 :          nparticle = force_env_get_nparticle(gopt_env%force_env)
     115        5280 :          ALLOCATE (x0(3*nparticle))
     116        1760 :          CALL pack_subsys_particles(subsys=subsys, r=x0)
     117             :       CASE (default_cell_method_id)
     118         378 :          SELECT CASE (gopt_env%cell_method_id)
     119             :          CASE (default_cell_direct_id)
     120         168 :             CALL force_env_get(gopt_env%force_env, subsys=subsys, cell=cell)
     121             :             ! Store reference cell
     122        4368 :             gopt_env%h_ref = cell%hmat
     123             :             ! before starting we handle the case of translating coordinates (QM/MM)
     124         168 :             IF (gopt_env%force_env%in_use == use_qmmm) &
     125           0 :                CALL apply_qmmm_translate(gopt_env%force_env%qmmm_env)
     126         168 :             IF (gopt_env%force_env%in_use == use_qmmmx) &
     127           0 :                CALL apply_qmmmx_translate(gopt_env%force_env%qmmmx_env)
     128         168 :             nparticle = force_env_get_nparticle(gopt_env%force_env)
     129         504 :             ALLOCATE (x0(3*nparticle + 6))
     130         168 :             CALL pack_subsys_particles(subsys=subsys, r=x0)
     131         168 :             idg = 3*nparticle
     132         672 :             DO i = 1, 3
     133        1680 :                DO j = 1, i
     134        1008 :                   idg = idg + 1
     135        1512 :                   x0(idg) = cell%hmat(j, i)
     136             :                END DO
     137             :             END DO
     138             :          CASE (default_cell_geo_opt_id, default_cell_md_id)
     139          42 :             CALL force_env_get(gopt_env%force_env, cell=cell)
     140          42 :             ALLOCATE (x0(6))
     141          42 :             idg = 0
     142         168 :             DO i = 1, 3
     143         420 :                DO j = 1, i
     144         252 :                   idg = idg + 1
     145         378 :                   x0(idg) = cell%hmat(j, i)
     146             :                END DO
     147             :             END DO
     148             :          CASE DEFAULT
     149         210 :             CPABORT("")
     150             :          END SELECT
     151             :       CASE DEFAULT
     152        1970 :          CPABORT("")
     153             :       END SELECT
     154             : 
     155        1970 :    END SUBROUTINE gopt_f_create_x0
     156             : 
     157             : ! **************************************************************************************************
     158             : !> \brief Prints iteration step of the optimization procedure on screen
     159             : !> \param its ...
     160             : !> \param output_unit ...
     161             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     162             : ! **************************************************************************************************
     163       13979 :    SUBROUTINE gopt_f_ii(its, output_unit)
     164             : 
     165             :       INTEGER, INTENT(IN)                                :: its, output_unit
     166             : 
     167       13979 :       IF (output_unit > 0) THEN
     168        6606 :          WRITE (UNIT=output_unit, FMT="(/,T2,26('-'))")
     169        6606 :          WRITE (UNIT=output_unit, FMT="(T2,A,I6)") "OPTIMIZATION STEP: ", its
     170        6606 :          WRITE (UNIT=output_unit, FMT="(T2,26('-'))")
     171        6606 :          CALL m_flush(output_unit)
     172             :       END IF
     173             : 
     174       13979 :    END SUBROUTINE gopt_f_ii
     175             : 
     176             : ! **************************************************************************************************
     177             : !> \brief Handles the Output during an optimization run
     178             : !> \param gopt_env ...
     179             : !> \param output_unit ...
     180             : !> \param opt_energy ...
     181             : !> \param wildcard ...
     182             : !> \param its ...
     183             : !> \param used_time ...
     184             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     185             : ! **************************************************************************************************
     186        1738 :    SUBROUTINE gopt_f_io_init(gopt_env, output_unit, opt_energy, wildcard, its, used_time)
     187             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     188             :       INTEGER, INTENT(IN)                                :: output_unit
     189             :       REAL(KIND=dp)                                      :: opt_energy
     190             :       CHARACTER(LEN=5)                                   :: wildcard
     191             :       INTEGER, INTENT(IN)                                :: its
     192             :       REAL(KIND=dp)                                      :: used_time
     193             : 
     194             :       REAL(KIND=dp)                                      :: pres_int
     195             : 
     196        3292 :       SELECT CASE (gopt_env%type_id)
     197             :       CASE (default_ts_method_id, default_minimization_method_id)
     198             :          ! Geometry Optimization (Minimization and Transition State Search)
     199        1554 :          IF (.NOT. gopt_env%dimer_rotation) THEN
     200             :             CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, wildcard=wildcard, &
     201        1422 :                                    used_time=used_time)
     202             :          ELSE
     203             :             CALL write_rot_cycle_infos(output_unit, it=its, etot=opt_energy, dimer_env=gopt_env%dimer_env, &
     204         132 :                                        wildcard=wildcard, used_time=used_time)
     205             :          END IF
     206             :       CASE (default_cell_method_id)
     207             :          ! Cell Optimization
     208         164 :          pres_int = gopt_env%cell_env%pres_int
     209             :          CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, pres_int=pres_int, wildcard=wildcard, &
     210         164 :                                 used_time=used_time)
     211             :       CASE (default_shellcore_method_id)
     212             :          CALL write_cycle_infos(output_unit, it=its, etot=opt_energy, wildcard=wildcard, &
     213        1738 :                                 used_time=used_time)
     214             :       END SELECT
     215             : 
     216        1738 :    END SUBROUTINE gopt_f_io_init
     217             : 
     218             : ! **************************************************************************************************
     219             : !> \brief Handles the Output during an optimization run
     220             : !> \param gopt_env ...
     221             : !> \param force_env ...
     222             : !> \param root_section ...
     223             : !> \param its ...
     224             : !> \param opt_energy ...
     225             : !> \param output_unit ...
     226             : !> \param eold ...
     227             : !> \param emin ...
     228             : !> \param wildcard ...
     229             : !> \param gopt_param ...
     230             : !> \param ndf ...
     231             : !> \param dx ...
     232             : !> \param xi ...
     233             : !> \param conv ...
     234             : !> \param pred ...
     235             : !> \param rat ...
     236             : !> \param step ...
     237             : !> \param rad ...
     238             : !> \param used_time ...
     239             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     240             : ! **************************************************************************************************
     241       27958 :    SUBROUTINE gopt_f_io(gopt_env, force_env, root_section, its, opt_energy, &
     242       13979 :                         output_unit, eold, emin, wildcard, gopt_param, ndf, dx, xi, conv, pred, rat, &
     243             :                         step, rad, used_time)
     244             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     245             :       TYPE(force_env_type), POINTER                      :: force_env
     246             :       TYPE(section_vals_type), POINTER                   :: root_section
     247             :       INTEGER, INTENT(IN)                                :: its
     248             :       REAL(KIND=dp), INTENT(IN)                          :: opt_energy
     249             :       INTEGER, INTENT(IN)                                :: output_unit
     250             :       REAL(KIND=dp)                                      :: eold, emin
     251             :       CHARACTER(LEN=5)                                   :: wildcard
     252             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     253             :       INTEGER, INTENT(IN), OPTIONAL                      :: ndf
     254             :       REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL  :: dx
     255             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: xi
     256             :       LOGICAL, OPTIONAL                                  :: conv
     257             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pred, rat, step, rad
     258             :       REAL(KIND=dp)                                      :: used_time
     259             : 
     260             :       INTEGER(KIND=int_8)                                :: max_memory
     261             :       LOGICAL                                            :: print_memory
     262             :       REAL(KIND=dp)                                      :: pres_diff, pres_diff_constr, pres_int, &
     263             :                                                             pres_tol
     264             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     265             : 
     266       13979 :       NULLIFY (para_env)
     267       13979 :       CALL section_vals_val_get(gopt_env%motion_section, "PRINT%MEMORY_INFO", l_val=print_memory)
     268       13979 :       max_memory = 0
     269       13979 :       IF (print_memory) THEN
     270       13979 :          CALL force_env_get(force_env, para_env=para_env)
     271       13979 :          max_memory = sample_memory(para_env)
     272             :       END IF
     273             : 
     274       23692 :       SELECT CASE (gopt_env%type_id)
     275             :       CASE (default_ts_method_id, default_minimization_method_id)
     276             :          ! Geometry Optimization (Minimization and Transition State Search)
     277        9713 :          IF (.NOT. gopt_env%dimer_rotation) THEN
     278             :             CALL geo_opt_io(force_env=force_env, root_section=root_section, &
     279        8987 :                             motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
     280             :             CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
     281        8987 :                                    pred=pred, rat=rat, step=step, rad=rad, emin=emin, wildcard=wildcard, used_time=used_time)
     282             :             ! Possibly check convergence
     283        8987 :             IF (PRESENT(conv)) THEN
     284        8987 :                CPASSERT(PRESENT(ndf))
     285        8987 :                CPASSERT(PRESENT(dx))
     286        8987 :                CPASSERT(PRESENT(xi))
     287        8987 :                CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory)
     288             :             END IF
     289             :          ELSE
     290         726 :             CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
     291         726 :             CALL write_restart(force_env=force_env, root_section=root_section)
     292             :             CALL write_rot_cycle_infos(output_unit, its, opt_energy, opt_energy - eold, emin, gopt_env%dimer_env, &
     293         726 :                                        used_time=used_time, wildcard=wildcard)
     294             :             ! Possibly check convergence
     295         726 :             IF (PRESENT(conv)) THEN
     296         726 :                CPASSERT(ASSOCIATED(gopt_env%dimer_env))
     297         726 :                CALL check_rot_conv(gopt_env%dimer_env, output_unit, conv)
     298             :             END IF
     299             :          END IF
     300             :       CASE (default_cell_method_id)
     301             :          ! Cell Optimization
     302        4096 :          pres_diff = gopt_env%cell_env%pres_int - gopt_env%cell_env%pres_ext
     303        4096 :          pres_int = gopt_env%cell_env%pres_int
     304        4096 :          pres_tol = gopt_env%cell_env%pres_tol
     305             :          CALL geo_opt_io(force_env=force_env, root_section=root_section, &
     306        4096 :                          motion_section=gopt_env%motion_section, its=its, opt_energy=opt_energy)
     307             :          CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
     308             :                                 pred=pred, rat=rat, step=step, rad=rad, emin=emin, pres_int=pres_int, wildcard=wildcard, &
     309        4096 :                                 used_time=used_time)
     310             :          ! Possibly check convergence
     311        4096 :          IF (PRESENT(conv)) THEN
     312        4096 :             CPASSERT(PRESENT(ndf))
     313        4096 :             CPASSERT(PRESENT(dx))
     314        4096 :             CPASSERT(PRESENT(xi))
     315        4096 :             IF (gopt_env%cell_env%constraint_id == fix_none) THEN
     316        4080 :                CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, pres_diff, pres_tol)
     317             :             ELSE
     318          16 :                pres_diff_constr = gopt_env%cell_env%pres_constr - gopt_env%cell_env%pres_ext
     319             :                CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory, pres_diff, &
     320          16 :                                   pres_tol, pres_diff_constr)
     321             :             END IF
     322             :          END IF
     323             :       CASE (default_shellcore_method_id)
     324             :          CALL write_cycle_infos(output_unit, its, etot=opt_energy, ediff=opt_energy - eold, &
     325         170 :                                 pred=pred, rat=rat, step=step, rad=rad, emin=emin, wildcard=wildcard, used_time=used_time)
     326             :          ! Possibly check convergence
     327       14149 :          IF (PRESENT(conv)) THEN
     328         170 :             CPASSERT(PRESENT(ndf))
     329         170 :             CPASSERT(PRESENT(dx))
     330         170 :             CPASSERT(PRESENT(xi))
     331         170 :             CALL check_converg(ndf, dx, xi, output_unit, conv, gopt_param, max_memory)
     332             :          END IF
     333             :       END SELECT
     334             : 
     335       13979 :    END SUBROUTINE gopt_f_io
     336             : 
     337             : ! **************************************************************************************************
     338             : !> \brief Handles the Output at the end of an optimization run
     339             : !> \param gopt_env ...
     340             : !> \param force_env ...
     341             : !> \param x0 ...
     342             : !> \param conv ...
     343             : !> \param its ...
     344             : !> \param root_section ...
     345             : !> \param para_env ...
     346             : !> \param master ...
     347             : !> \param output_unit ...
     348             : !> \author Teodoro Laino [tlaino] - University of Zurich - 03.2008
     349             : ! **************************************************************************************************
     350        2122 :    RECURSIVE SUBROUTINE gopt_f_io_finalize(gopt_env, force_env, x0, conv, its, root_section, &
     351             :                                            para_env, master, output_unit)
     352             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     353             :       TYPE(force_env_type), POINTER                      :: force_env
     354             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
     355             :       LOGICAL                                            :: conv
     356             :       INTEGER                                            :: its
     357             :       TYPE(section_vals_type), POINTER                   :: root_section
     358             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     359             :       INTEGER, INTENT(IN)                                :: master, output_unit
     360             : 
     361        2122 :       IF (gopt_env%eval_opt_geo) THEN
     362        1214 :          IF (.NOT. gopt_env%dimer_rotation) THEN
     363             :             CALL write_final_info(output_unit, conv, its, gopt_env, x0, master, &
     364        1082 :                                   para_env, force_env, gopt_env%motion_section, root_section)
     365             :          ELSE
     366         132 :             CALL update_dimer_vec(gopt_env%dimer_env, gopt_env%motion_section)
     367         132 :             CALL write_restart(force_env=force_env, root_section=root_section)
     368             :          END IF
     369             :       END IF
     370             : 
     371        2122 :    END SUBROUTINE gopt_f_io_finalize
     372             : 
     373             : ! **************************************************************************************************
     374             : !> \brief ...
     375             : !> \param output_unit ...
     376             : !> \param it ...
     377             : !> \param etot ...
     378             : !> \param ediff ...
     379             : !> \param pred ...
     380             : !> \param rat ...
     381             : !> \param step ...
     382             : !> \param rad ...
     383             : !> \param emin ...
     384             : !> \param pres_int ...
     385             : !> \param wildcard ...
     386             : !> \param used_time ...
     387             : ! **************************************************************************************************
     388       14859 :    SUBROUTINE write_cycle_infos(output_unit, it, etot, ediff, pred, rat, step, rad, emin, &
     389             :                                 pres_int, wildcard, used_time)
     390             : 
     391             :       INTEGER, INTENT(IN)                                :: output_unit, it
     392             :       REAL(KIND=dp), INTENT(IN)                          :: etot
     393             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: ediff, pred, rat, step, rad, emin, &
     394             :                                                             pres_int
     395             :       CHARACTER(LEN=5), OPTIONAL                         :: wildcard
     396             :       REAL(KIND=dp), INTENT(IN)                          :: used_time
     397             : 
     398             :       REAL(KIND=dp)                                      :: tmp_r1
     399             : 
     400       14859 :       IF (output_unit > 0) THEN
     401             :          WRITE (UNIT=output_unit, FMT="(/,T2,8('-'),A,I5,1X,12('-'))") &
     402        7046 :             "  Informations at step = ", it
     403             :          WRITE (UNIT=output_unit, FMT="(T2,A,T47,A)") &
     404        7046 :             " Optimization Method        = ", wildcard
     405             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     406        7046 :             " Total Energy               = ", etot
     407        7046 :          IF (PRESENT(pres_int)) THEN
     408        2130 :             tmp_r1 = cp_unit_from_cp2k(pres_int, "bar")
     409             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     410        2130 :                " Internal Pressure [bar]    = ", tmp_r1
     411             :          END IF
     412        7046 :          IF (PRESENT(ediff)) THEN
     413             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     414        6243 :                " Real energy change         = ", ediff
     415             :          END IF
     416        7046 :          IF (PRESENT(pred)) THEN
     417             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     418        3295 :                " Predicted change in energy = ", pred
     419             :          END IF
     420        7046 :          IF (PRESENT(rat)) THEN
     421             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     422        3295 :                " Scaling factor             = ", rat
     423             :          END IF
     424        7046 :          IF (PRESENT(step)) THEN
     425             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     426        3295 :                " Step size                  = ", step
     427             :          END IF
     428        7046 :          IF (PRESENT(rad)) THEN
     429             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     430        3295 :                " Trust radius               = ", rad
     431             :          END IF
     432        7046 :          IF (PRESENT(emin)) THEN
     433        6243 :             IF (etot < emin) THEN
     434             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     435        5677 :                   " Decrease in energy         =     ", &
     436       11354 :                   "             YES"
     437             :             ELSE
     438             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     439         566 :                   " Decrease in energy         =     ", &
     440        1132 :                   "              NO"
     441             :             END IF
     442             :          END IF
     443             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.3)") &
     444        7046 :             " Used time                  = ", used_time
     445        7046 :          IF (it == 0) WRITE (UNIT=output_unit, FMT="(T2,51('-'))")
     446             :       END IF
     447       14859 :    END SUBROUTINE write_cycle_infos
     448             : 
     449             : ! **************************************************************************************************
     450             : !> \brief ...
     451             : !> \param output_unit ...
     452             : !> \param it ...
     453             : !> \param etot ...
     454             : !> \param ediff ...
     455             : !> \param emin ...
     456             : !> \param dimer_env ...
     457             : !> \param used_time ...
     458             : !> \param wildcard ...
     459             : !> \date  01.2008
     460             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
     461             : ! **************************************************************************************************
     462         858 :    SUBROUTINE write_rot_cycle_infos(output_unit, it, etot, ediff, emin, dimer_env, used_time, wildcard)
     463             : 
     464             :       INTEGER, INTENT(IN)                                :: output_unit, it
     465             :       REAL(KIND=dp), INTENT(IN)                          :: etot
     466             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: ediff, emin
     467             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     468             :       REAL(KIND=dp)                                      :: used_time
     469             :       CHARACTER(LEN=5), OPTIONAL                         :: wildcard
     470             : 
     471         858 :       IF (output_unit > 0) THEN
     472             :          WRITE (UNIT=output_unit, FMT="(/,T2,4('-'),A,I5,1X,5('-'))") &
     473         429 :             "  Informations at rotational step = ", it
     474             :          WRITE (UNIT=output_unit, FMT="(T2,A,T47,A)") &
     475         429 :             " Optimization Method        = ", wildcard
     476             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     477         429 :             " Local Curvature            = ", dimer_env%rot%curvature
     478             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     479         429 :             " Total Rotational Force     = ", etot
     480         429 :          IF (PRESENT(ediff)) THEN
     481             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     482         363 :                " Real Force change          = ", ediff
     483             :          END IF
     484         429 :          IF (PRESENT(emin)) THEN
     485         363 :             IF (etot < emin) THEN
     486             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     487         157 :                   " Decrease in rotational force =   ", &
     488         314 :                   "             YES"
     489             :             ELSE
     490             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     491         206 :                   " Decrease in rotational force =   ", &
     492         412 :                   "              NO"
     493             :             END IF
     494             :          END IF
     495             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.3)") &
     496         429 :             " Used time                  = ", used_time
     497         429 :          IF (it == 0) WRITE (UNIT=output_unit, FMT="(T2,51('-'))")
     498             :       END IF
     499         858 :    END SUBROUTINE write_rot_cycle_infos
     500             : 
     501             : ! **************************************************************************************************
     502             : !> \brief ...
     503             : !> \param ndf ...
     504             : !> \param dr ...
     505             : !> \param g ...
     506             : !> \param output_unit ...
     507             : !> \param conv ...
     508             : !> \param gopt_param ...
     509             : !> \param max_memory ...
     510             : !> \param pres_diff ...
     511             : !> \param pres_tol ...
     512             : !> \param pres_diff_constr ...
     513             : ! **************************************************************************************************
     514       13253 :    SUBROUTINE check_converg(ndf, dr, g, output_unit, conv, gopt_param, max_memory, pres_diff, &
     515             :                             pres_tol, pres_diff_constr)
     516             : 
     517             :       INTEGER, INTENT(IN)                                :: ndf
     518             :       REAL(KIND=dp), INTENT(IN)                          :: dr(ndf), g(ndf)
     519             :       INTEGER, INTENT(IN)                                :: output_unit
     520             :       LOGICAL, INTENT(OUT)                               :: conv
     521             :       TYPE(gopt_param_type), POINTER                     :: gopt_param
     522             :       INTEGER(KIND=int_8), INTENT(IN)                    :: max_memory
     523             :       REAL(KIND=dp), INTENT(IN), OPTIONAL                :: pres_diff, pres_tol, pres_diff_constr
     524             : 
     525             :       INTEGER                                            :: indf
     526             :       LOGICAL                                            :: conv_dx, conv_g, conv_p, conv_rdx, &
     527             :                                                             conv_rg
     528             :       REAL(KIND=dp)                                      :: dumm, dxcon, gcon, maxdum(4), rmsgcon, &
     529             :                                                             rmsxcon, tmp_r1
     530             : 
     531       13253 :       dxcon = gopt_param%max_dr
     532       13253 :       gcon = gopt_param%max_force
     533       13253 :       rmsgcon = gopt_param%rms_force
     534       13253 :       rmsxcon = gopt_param%rms_dr
     535             : 
     536       13253 :       conv = .FALSE.
     537       13253 :       conv_dx = .TRUE.
     538       13253 :       conv_rdx = .TRUE.
     539       13253 :       conv_g = .TRUE.
     540       13253 :       conv_rg = .TRUE.
     541       13253 :       conv_p = .TRUE.
     542             : 
     543       13253 :       dumm = 0.0_dp
     544     2951318 :       DO indf = 1, ndf
     545     2938065 :          IF (indf == 1) maxdum(1) = ABS(dr(indf))
     546     2938065 :          dumm = dumm + dr(indf)**2
     547     2938065 :          IF (ABS(dr(indf)) > dxcon) conv_dx = .FALSE.
     548     2951318 :          IF (ABS(dr(indf)) > maxdum(1)) maxdum(1) = ABS(dr(indf))
     549             :       END DO
     550             :       ! SQRT(dumm/ndf) > rmsxcon
     551       13253 :       IF (dumm > (rmsxcon*rmsxcon*ndf)) conv_rdx = .FALSE.
     552       13253 :       maxdum(2) = SQRT(dumm/ndf)
     553             : 
     554       13253 :       dumm = 0.0_dp
     555     2951318 :       DO indf = 1, ndf
     556     2938065 :          IF (indf == 1) maxdum(3) = ABS(g(indf))
     557     2938065 :          dumm = dumm + g(indf)**2
     558     2938065 :          IF (ABS(g(indf)) > gcon) conv_g = .FALSE.
     559     2951318 :          IF (ABS(g(indf)) > maxdum(3)) maxdum(3) = ABS(g(indf))
     560             :       END DO
     561             :       ! SQRT(dumm/ndf) > rmsgcon
     562       13253 :       IF (dumm > (rmsgcon*rmsgcon*ndf)) conv_rg = .FALSE.
     563       13253 :       maxdum(4) = SQRT(dumm/ndf)
     564             : 
     565       13253 :       IF (PRESENT(pres_diff_constr) .AND. PRESENT(pres_tol)) THEN
     566          16 :          conv_p = ABS(pres_diff_constr) < ABS(pres_tol)
     567       13237 :       ELSEIF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
     568        4080 :          conv_p = ABS(pres_diff) < ABS(pres_tol)
     569             :       END IF
     570             : 
     571       13253 :       IF (output_unit > 0) THEN
     572             :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     573        6243 :             " Convergence check :"
     574             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     575        6243 :             " Max. step size             = ", maxdum(1)
     576             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     577        6243 :             " Conv. limit for step size  = ", dxcon
     578        6243 :          IF (conv_dx) THEN
     579             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     580        1435 :                " Convergence in step size   =     ", &
     581        2870 :                "             YES"
     582             :          ELSE
     583             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     584        4808 :                " Convergence in step size   =     ", &
     585        9616 :                "              NO"
     586             :          END IF
     587             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     588        6243 :             " RMS step size              = ", maxdum(2)
     589             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     590        6243 :             " Conv. limit for RMS step   = ", rmsxcon
     591        6243 :          IF (conv_rdx) THEN
     592             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     593        1569 :                " Convergence in RMS step    =     ", &
     594        3138 :                "             YES"
     595             :          ELSE
     596             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     597        4674 :                " Convergence in RMS step    =     ", &
     598        9348 :                "              NO"
     599             :          END IF
     600             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     601        6243 :             " Max. gradient              = ", maxdum(3)
     602             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     603        6243 :             " Conv. limit for gradients  = ", gcon
     604        6243 :          IF (conv_g) THEN
     605             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     606        1477 :                " Conv. in gradients         =     ", &
     607        2954 :                "             YES"
     608             :          ELSE
     609             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     610        4766 :                " Conv. for gradients        =     ", &
     611        9532 :                "              NO"
     612             :          END IF
     613             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     614        6243 :             " RMS gradient               = ", maxdum(4)
     615             :          WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     616        6243 :             " Conv. limit for RMS grad.  = ", rmsgcon
     617        6243 :          IF (conv_rg) THEN
     618             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     619        1771 :                " Conv. in RMS gradients     =     ", &
     620        3542 :                "             YES"
     621             : 
     622             :          ELSE
     623             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     624        4472 :                " Conv. for gradients        =     ", &
     625        8944 :                "              NO"
     626             :          END IF
     627        6243 :          IF (PRESENT(pres_diff) .AND. PRESENT(pres_tol)) THEN
     628        2048 :             tmp_r1 = cp_unit_from_cp2k(pres_diff, "bar")
     629        2048 :             IF (PRESENT(pres_diff_constr)) THEN
     630             :                WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     631           8 :                   " Pressure Deviation [bar]   = ", tmp_r1, &
     632          16 :                   " (without constraint)"
     633           8 :                tmp_r1 = cp_unit_from_cp2k(pres_diff_constr, "bar")
     634             :                WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     635           8 :                   " Pressure Deviation [bar]   = ", tmp_r1, &
     636          16 :                   " (with constraint)"
     637             :             ELSE
     638             :                WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     639        2040 :                   " Pressure Deviation [bar]   = ", tmp_r1
     640             :             END IF
     641        2048 :             tmp_r1 = cp_unit_from_cp2k(pres_tol, "bar")
     642             :             WRITE (UNIT=output_unit, FMT="(T2,A,F20.10)") &
     643        2048 :                " Pressure Tolerance [bar]   = ", tmp_r1
     644        2048 :             IF (conv_p) THEN
     645             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     646         311 :                   " Conv. for  PRESSURE        =     ", &
     647         622 :                   "             YES"
     648             :             ELSE
     649             :                WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     650        1737 :                   " Conv. for  PRESSURE        =     ", &
     651        3474 :                   "              NO"
     652             :             END IF
     653             :          END IF
     654        6243 :          WRITE (UNIT=output_unit, FMT="(T2,51('-'))")
     655        6243 :          IF (max_memory .NE. 0) THEN
     656             :             WRITE (output_unit, '(T2,A,T73,I8)') &
     657        6243 :                'Estimated peak process memory after this step [MiB]', &
     658       12486 :                (max_memory + (1024*1024) - 1)/(1024*1024)
     659             :          END IF
     660             :       END IF
     661             : 
     662       13253 :       IF (conv_dx .AND. conv_rdx .AND. conv_g .AND. conv_rg .AND. conv_p) conv = .TRUE.
     663             : 
     664       13253 :       IF ((conv) .AND. (output_unit > 0)) THEN
     665         649 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     666             :          WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
     667         649 :             "***", "GEOMETRY OPTIMIZATION COMPLETED", "***"
     668         649 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     669             :       END IF
     670             : 
     671       13253 :    END SUBROUTINE check_converg
     672             : 
     673             : ! **************************************************************************************************
     674             : !> \brief ...
     675             : !> \param dimer_env ...
     676             : !> \param output_unit ...
     677             : !> \param conv ...
     678             : !> \date  01.2008
     679             : !> \author Luca Bellucci and Teodoro Laino - created [tlaino]
     680             : ! **************************************************************************************************
     681         726 :    SUBROUTINE check_rot_conv(dimer_env, output_unit, conv)
     682             : 
     683             :       TYPE(dimer_env_type), POINTER                      :: dimer_env
     684             :       INTEGER, INTENT(IN)                                :: output_unit
     685             :       LOGICAL, INTENT(OUT)                               :: conv
     686             : 
     687         726 :       conv = (ABS(dimer_env%rot%angle2) < dimer_env%rot%angle_tol)
     688         726 :       IF (output_unit > 0) THEN
     689             :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     690         363 :             " Convergence check :"
     691             :          WRITE (UNIT=output_unit, FMT="(T2,A,F16.10)") &
     692         363 :             " Predicted angle step size  =     ", dimer_env%rot%angle1
     693             :          WRITE (UNIT=output_unit, FMT="(T2,A,F16.10)") &
     694         363 :             " Effective angle step size  =     ", dimer_env%rot%angle2
     695             :          WRITE (UNIT=output_unit, FMT="(T2,A,F16.10)") &
     696         363 :             " Conv. limit for angle step size =", dimer_env%rot%angle_tol
     697         363 :          IF (conv) THEN
     698             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     699          55 :                " Convergence in angle step size   =", &
     700         110 :                "            YES"
     701             :          ELSE
     702             :             WRITE (UNIT=output_unit, FMT="(T2,2A)") &
     703         308 :                " Convergence in angle step size   =", &
     704         616 :                "             NO"
     705             :          END IF
     706         363 :          WRITE (UNIT=output_unit, FMT="(T2,51('-'))")
     707             :       END IF
     708         726 :       IF ((conv) .AND. (output_unit > 0)) THEN
     709          55 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     710             :          WRITE (UNIT=output_unit, FMT="(T2,A,T25,A,T78,A)") &
     711          55 :             "***", "ROTATION OPTIMIZATION COMPLETED", "***"
     712          55 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     713             :       END IF
     714             : 
     715         726 :    END SUBROUTINE check_rot_conv
     716             : 
     717             : ! **************************************************************************************************
     718             : !> \brief ...
     719             : !> \param output_unit ...
     720             : !> \param conv ...
     721             : !> \param it ...
     722             : !> \param gopt_env ...
     723             : !> \param x0 ...
     724             : !> \param master ...
     725             : !> \param para_env ...
     726             : !> \param force_env ...
     727             : !> \param motion_section ...
     728             : !> \param root_section ...
     729             : !> \date  11.2007
     730             : !> \author Teodoro Laino [tlaino] - University of Zurich
     731             : ! **************************************************************************************************
     732        1082 :    RECURSIVE SUBROUTINE write_final_info(output_unit, conv, it, gopt_env, x0, master, para_env, force_env, &
     733             :                                          motion_section, root_section)
     734             :       INTEGER, INTENT(IN)                                :: output_unit
     735             :       LOGICAL, INTENT(IN)                                :: conv
     736             :       INTEGER, INTENT(INOUT)                             :: it
     737             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     738             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x0
     739             :       INTEGER, INTENT(IN)                                :: master
     740             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     741             :       TYPE(force_env_type), POINTER                      :: force_env
     742             :       TYPE(section_vals_type), POINTER                   :: motion_section, root_section
     743             : 
     744             :       REAL(KIND=dp)                                      :: etot
     745             :       TYPE(cell_type), POINTER                           :: cell
     746             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     747             :       TYPE(particle_list_type), POINTER                  :: particles
     748        1082 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     749             : 
     750        1082 :       CALL force_env_get(force_env, cell=cell, subsys=subsys)
     751        1082 :       CALL cp_subsys_get(subsys=subsys, particles=particles)
     752        1082 :       particle_set => particles%els
     753        1082 :       IF (conv) THEN
     754         403 :          it = it + 1
     755         403 :          CALL write_structure_data(particle_set, cell, motion_section)
     756         403 :          CALL write_restart(force_env=force_env, root_section=root_section)
     757             : 
     758         403 :          IF (output_unit > 0) &
     759         201 :             WRITE (UNIT=output_unit, FMT="(/,T20,' Reevaluating energy at the minimum')")
     760             : 
     761             :          CALL cp_eval_at(gopt_env, x0, f=etot, master=master, final_evaluation=.TRUE., &
     762         403 :                          para_env=para_env)
     763         403 :          CALL write_geo_traj(force_env, root_section, it, etot)
     764             :       END IF
     765             : 
     766        1082 :    END SUBROUTINE write_final_info
     767             : 
     768             : ! **************************************************************************************************
     769             : !> \brief  Specific driver for dumping trajectory during a GEO_OPT
     770             : !> \param force_env ...
     771             : !> \param root_section ...
     772             : !> \param it ...
     773             : !> \param etot ...
     774             : !> \date   11.2007
     775             : !> \par    History
     776             : !>         09.2010: Output of core and shell positions and forces (MK)
     777             : !> \author Teodoro Laino [tlaino] - University of Zurich
     778             : ! **************************************************************************************************
     779       26972 :    SUBROUTINE write_geo_traj(force_env, root_section, it, etot)
     780             : 
     781             :       TYPE(force_env_type), POINTER                      :: force_env
     782             :       TYPE(section_vals_type), POINTER                   :: root_section
     783             :       INTEGER, INTENT(IN)                                :: it
     784             :       REAL(KIND=dp), INTENT(IN)                          :: etot
     785             : 
     786             :       LOGICAL                                            :: shell_adiabatic, shell_present
     787             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     788       13486 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     789             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     790             :       TYPE(particle_list_type), POINTER                  :: core_particles, shell_particles
     791             : 
     792       13486 :       NULLIFY (atomic_kinds)
     793       13486 :       NULLIFY (atomic_kind_set)
     794       13486 :       NULLIFY (core_particles)
     795       13486 :       NULLIFY (shell_particles)
     796       13486 :       NULLIFY (subsys)
     797             : 
     798       13486 :       CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot)
     799             :       ! Print Force
     800       13486 :       CALL write_trajectory(force_env, root_section, it, 0.0_dp, 0.0_dp, etot, "FORCES", middle_name="frc")
     801       13486 :       CALL force_env_get(force_env, subsys=subsys)
     802       13486 :       CALL cp_subsys_get(subsys, atomic_kinds=atomic_kinds)
     803       13486 :       atomic_kind_set => atomic_kinds%els
     804             :       CALL get_atomic_kind_set(atomic_kind_set, &
     805             :                                shell_present=shell_present, &
     806       13486 :                                shell_adiabatic=shell_adiabatic)
     807       13486 :       IF (shell_present) THEN
     808             :          CALL cp_subsys_get(subsys, &
     809             :                             core_particles=core_particles, &
     810        4120 :                             shell_particles=shell_particles)
     811             :          CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
     812             :                                etot=etot, pk_name="SHELL_TRAJECTORY", middle_name="shpos", &
     813        4120 :                                particles=shell_particles)
     814        4120 :          IF (shell_adiabatic) THEN
     815             :             CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
     816             :                                   etot=etot, pk_name="SHELL_FORCES", middle_name="shfrc", &
     817        4120 :                                   particles=shell_particles)
     818             :             CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
     819             :                                   etot=etot, pk_name="CORE_TRAJECTORY", middle_name="copos", &
     820        4120 :                                   particles=core_particles)
     821             :             CALL write_trajectory(force_env, root_section, it=it, time=0.0_dp, dtime=0.0_dp, &
     822             :                                   etot=etot, pk_name="CORE_FORCES", middle_name="cofrc", &
     823        4120 :                                   particles=core_particles)
     824             :          END IF
     825             :       END IF
     826             : 
     827       13486 :    END SUBROUTINE write_geo_traj
     828             : 
     829             : ! **************************************************************************************************
     830             : !> \brief ...
     831             : !> \param gopt_env ...
     832             : !> \param output_unit ...
     833             : !> \param label ...
     834             : !> \date  01.2008
     835             : !> \author Teodoro Laino [tlaino] - University of Zurich
     836             : ! **************************************************************************************************
     837        2122 :    SUBROUTINE print_geo_opt_header(gopt_env, output_unit, label)
     838             : 
     839             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     840             :       INTEGER, INTENT(IN)                                :: output_unit
     841             :       CHARACTER(LEN=*), INTENT(IN)                       :: label
     842             : 
     843             :       CHARACTER(LEN=default_string_length)               :: my_format, my_label
     844             :       INTEGER                                            :: ix
     845             : 
     846        2122 :       IF (output_unit > 0) THEN
     847        1061 :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") REPEAT("*", 79)
     848        1061 :          IF (gopt_env%dimer_rotation) THEN
     849          66 :             my_label = "OPTIMIZING DIMER ROTATION"
     850             :          ELSE
     851         995 :             my_label = "STARTING "//gopt_env%tag(1:8)//" OPTIMIZATION"
     852             :          END IF
     853             : 
     854        1061 :          ix = (80 - 7 - LEN_TRIM(my_label))/2
     855        1061 :          ix = ix + 5
     856        1061 :          my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
     857        1061 :          WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(my_label), "***"
     858             : 
     859        1061 :          ix = (80 - 7 - LEN_TRIM(label))/2
     860        1061 :          ix = ix + 5
     861        1061 :          my_format = "(T2,A,T"//cp_to_string(ix)//",A,T78,A)"
     862        1061 :          WRITE (UNIT=output_unit, FMT=TRIM(my_format)) "***", TRIM(label), "***"
     863             : 
     864        1061 :          WRITE (UNIT=output_unit, FMT="(T2,A)") REPEAT("*", 79)
     865        1061 :          CALL m_flush(output_unit)
     866             :       END IF
     867        2122 :    END SUBROUTINE print_geo_opt_header
     868             : 
     869             : ! **************************************************************************************************
     870             : !> \brief ...
     871             : !> \param gopt_env ...
     872             : !> \param output_unit ...
     873             : !> \date  01.2008
     874             : !> \author Teodoro Laino [tlaino] - University of Zurich
     875             : ! **************************************************************************************************
     876         701 :    SUBROUTINE print_geo_opt_nc(gopt_env, output_unit)
     877             : 
     878             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     879             :       INTEGER, INTENT(IN)                                :: output_unit
     880             : 
     881         701 :       IF (output_unit > 0) THEN
     882             :          WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
     883         351 :             "*** MAXIMUM NUMBER OF OPTIMIZATION STEPS REACHED ***"
     884         351 :          IF (.NOT. gopt_env%dimer_rotation) THEN
     885             :             WRITE (UNIT=output_unit, FMT="(T2,A)") &
     886         340 :                "***        EXITING GEOMETRY OPTIMIZATION         ***"
     887             :          ELSE
     888             :             WRITE (UNIT=output_unit, FMT="(T2,A)") &
     889          11 :                "***        EXITING ROTATION OPTIMIZATION         ***"
     890             :          END IF
     891         351 :          CALL m_flush(output_unit)
     892             :       END IF
     893             : 
     894         701 :    END SUBROUTINE print_geo_opt_nc
     895             : 
     896             : ! **************************************************************************************************
     897             : !> \brief   Prints information during GEO_OPT common to all optimizers
     898             : !> \param force_env ...
     899             : !> \param root_section ...
     900             : !> \param motion_section ...
     901             : !> \param its ...
     902             : !> \param opt_energy ...
     903             : !> \date    02.2008
     904             : !> \author  Teodoro Laino [tlaino] - University of Zurich
     905             : !> \version 1.0
     906             : ! **************************************************************************************************
     907       13083 :    SUBROUTINE geo_opt_io(force_env, root_section, motion_section, its, opt_energy)
     908             : 
     909             :       TYPE(force_env_type), POINTER                      :: force_env
     910             :       TYPE(section_vals_type), POINTER                   :: root_section, motion_section
     911             :       INTEGER, INTENT(IN)                                :: its
     912             :       REAL(KIND=dp), INTENT(IN)                          :: opt_energy
     913             : 
     914             :       TYPE(atomic_kind_list_type), POINTER               :: atomic_kinds
     915       13083 :       TYPE(atomic_kind_type), DIMENSION(:), POINTER      :: atomic_kind_set
     916             :       TYPE(cell_type), POINTER                           :: cell
     917             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     918             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     919             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     920             :       TYPE(particle_list_type), POINTER                  :: particles
     921       13083 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     922             :       TYPE(virial_type), POINTER                         :: virial
     923             : 
     924       13083 :       NULLIFY (para_env, atomic_kind_set, subsys, particle_set, &
     925       13083 :                local_particles, atomic_kinds, particles)
     926             : 
     927             :       ! Write Restart File
     928       13083 :       CALL write_restart(force_env=force_env, root_section=root_section)
     929             : 
     930             :       ! Write Trajectory
     931       13083 :       CALL write_geo_traj(force_env, root_section, its, opt_energy)
     932             : 
     933             :       ! Write the stress Tensor
     934             :       CALL force_env_get(force_env, cell=cell, para_env=para_env, &
     935       13083 :                          subsys=subsys)
     936             :       CALL cp_subsys_get(subsys=subsys, atomic_kinds=atomic_kinds, local_particles=local_particles, &
     937       13083 :                          particles=particles, virial=virial)
     938       13083 :       atomic_kind_set => atomic_kinds%els
     939       13083 :       particle_set => particles%els
     940             :       CALL virial_evaluate(atomic_kind_set, particle_set, local_particles, &
     941       13083 :                            virial, para_env)
     942       13083 :       CALL write_stress_tensor(virial, cell, motion_section, its, 0.0_dp)
     943             : 
     944             :       ! Write the cell
     945       13083 :       CALL write_simulation_cell(cell, motion_section, its, 0.0_dp)
     946             : 
     947       13083 :    END SUBROUTINE geo_opt_io
     948             : 
     949             : ! **************************************************************************************************
     950             : !> \brief   Apply coordinate transformations after cell (shape) change
     951             : !> \param gopt_env ...
     952             : !> \param cell ...
     953             : !> \param x ...
     954             : !> \param update_forces ...
     955             : !> \date    05.11.2012 (revised version of unbiase_coordinates moved here, MK)
     956             : !> \author  Matthias Krack
     957             : !> \version 1.0
     958             : ! **************************************************************************************************
     959        8904 :    SUBROUTINE apply_cell_change(gopt_env, cell, x, update_forces)
     960             : 
     961             :       TYPE(gopt_f_type), POINTER                         :: gopt_env
     962             :       TYPE(cell_type), POINTER                           :: cell
     963             :       REAL(KIND=dp), DIMENSION(:), POINTER               :: x
     964             :       LOGICAL, INTENT(IN)                                :: update_forces
     965             : 
     966             :       INTEGER                                            :: i, iatom, idg, j, natom, nparticle, &
     967             :                                                             shell_index
     968             :       REAL(KIND=dp)                                      :: fc, fs, mass
     969             :       REAL(KIND=dp), DIMENSION(3)                        :: s
     970             :       TYPE(cell_type), POINTER                           :: cell_ref
     971             :       TYPE(cp_subsys_type), POINTER                      :: subsys
     972             :       TYPE(particle_list_type), POINTER                  :: core_particles, particles, &
     973             :                                                             shell_particles
     974             : 
     975        8904 :       NULLIFY (cell_ref)
     976        8904 :       NULLIFY (core_particles)
     977        8904 :       NULLIFY (particles)
     978        8904 :       NULLIFY (shell_particles)
     979        8904 :       NULLIFY (subsys)
     980             : 
     981        8904 :       natom = force_env_get_natom(gopt_env%force_env)
     982        8904 :       nparticle = force_env_get_nparticle(gopt_env%force_env)
     983             :       CALL force_env_get(gopt_env%force_env, &
     984        8904 :                          subsys=subsys)
     985             :       CALL cp_subsys_get(subsys=subsys, &
     986             :                          core_particles=core_particles, &
     987             :                          particles=particles, &
     988        8904 :                          shell_particles=shell_particles)
     989             : 
     990             :       ! Retrieve the reference cell
     991        8904 :       CALL cell_create(cell_ref)
     992        8904 :       CALL cell_copy(cell, cell_ref, tag="CELL_OPT_REF")
     993             : 
     994             :       ! Load the updated cell information
     995       16876 :       SELECT CASE (gopt_env%cell_method_id)
     996             :       CASE (default_cell_direct_id)
     997        7972 :          idg = 3*nparticle
     998        7972 :          CALL init_cell(cell_ref, hmat=gopt_env%h_ref)
     999             :       CASE (default_cell_geo_opt_id, default_cell_md_id)
    1000        8904 :          idg = 0
    1001             :       END SELECT
    1002        8904 :       CPASSERT((SIZE(x) == idg + 6))
    1003             : 
    1004        8904 :       IF (update_forces) THEN
    1005             : 
    1006             :          ! Transform particle forces back to reference cell
    1007             :          idg = 1
    1008      235244 :          DO iatom = 1, natom
    1009      232224 :             CALL real_to_scaled(s, x(idg:idg + 2), cell)
    1010      232224 :             CALL scaled_to_real(x(idg:idg + 2), s, cell_ref)
    1011      235244 :             idg = idg + 3
    1012             :          END DO
    1013             : 
    1014             :       ELSE
    1015             : 
    1016             :          ! Update cell
    1017       23536 :          DO i = 1, 3
    1018       58840 :             DO j = 1, i
    1019       35304 :                idg = idg + 1
    1020       52956 :                cell%hmat(j, i) = x(idg)
    1021             :             END DO
    1022             :          END DO
    1023        5884 :          CALL init_cell(cell)
    1024        5884 :          CALL cp_subsys_set(subsys, cell=cell)
    1025             : 
    1026             :          ! Retrieve particle coordinates for the current cell
    1027        5884 :          SELECT CASE (gopt_env%cell_method_id)
    1028             :          CASE (default_cell_direct_id)
    1029             :             idg = 1
    1030      441272 :             DO iatom = 1, natom
    1031      436320 :                CALL real_to_scaled(s, x(idg:idg + 2), cell_ref)
    1032      436320 :                shell_index = particles%els(iatom)%shell_index
    1033      436320 :                IF (shell_index == 0) THEN
    1034      144588 :                   CALL scaled_to_real(particles%els(iatom)%r, s, cell)
    1035             :                ELSE
    1036      291732 :                   CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
    1037      291732 :                   i = 3*(natom + shell_index - 1) + 1
    1038      291732 :                   CALL real_to_scaled(s, x(i:i + 2), cell_ref)
    1039      291732 :                   CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
    1040             :                   ! Update atomic position due to core and shell motion
    1041      291732 :                   mass = particles%els(iatom)%atomic_kind%mass
    1042      291732 :                   fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
    1043      291732 :                   fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
    1044             :                   particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
    1045     2333856 :                                                 fs*shell_particles%els(shell_index)%r(1:3)
    1046             :                END IF
    1047      441272 :                idg = idg + 3
    1048             :             END DO
    1049             :          CASE (default_cell_geo_opt_id, default_cell_md_id)
    1050       45852 :             DO iatom = 1, natom
    1051       39036 :                shell_index = particles%els(iatom)%shell_index
    1052       39968 :                IF (shell_index == 0) THEN
    1053       34932 :                   CALL real_to_scaled(s, particles%els(iatom)%r, cell_ref)
    1054       34932 :                   CALL scaled_to_real(particles%els(iatom)%r, s, cell)
    1055             :                ELSE
    1056        4104 :                   CALL real_to_scaled(s, core_particles%els(shell_index)%r, cell_ref)
    1057        4104 :                   CALL scaled_to_real(core_particles%els(shell_index)%r, s, cell)
    1058        4104 :                   i = 3*(natom + shell_index - 1) + 1
    1059        4104 :                   CALL real_to_scaled(s, shell_particles%els(shell_index)%r, cell_ref)
    1060        4104 :                   CALL scaled_to_real(shell_particles%els(shell_index)%r, s, cell)
    1061             :                   ! Update atomic position due to core and shell motion
    1062        4104 :                   mass = particles%els(iatom)%atomic_kind%mass
    1063        4104 :                   fc = core_particles%els(shell_index)%atomic_kind%shell%mass_core/mass
    1064        4104 :                   fs = shell_particles%els(shell_index)%atomic_kind%shell%mass_shell/mass
    1065             :                   particles%els(iatom)%r(1:3) = fc*core_particles%els(shell_index)%r(1:3) + &
    1066       32832 :                                                 fs*shell_particles%els(shell_index)%r(1:3)
    1067             :                END IF
    1068             :             END DO
    1069             :          END SELECT
    1070             : 
    1071             :       END IF
    1072             : 
    1073        8904 :       CALL cell_release(cell_ref)
    1074             : 
    1075        8904 :    END SUBROUTINE apply_cell_change
    1076             : 
    1077             : END MODULE gopt_f_methods

Generated by: LCOV version 1.15