LCOV - code coverage report
Current view: top level - src/motion - gopt_f_methods.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:dc34ec9) Lines: 383 386 99.2 %
Date: 2023-03-24 20:09:49 Functions: 15 15 100.0 %

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

Generated by: LCOV version 1.15