motion: Add space group symmetry

Add space group symmetry detection and preservation throughout geometry and cell optimization.

src/cryssym.F

@@ -165,9 +165,6 @@ SUBROUTINE crys_sym_gen(csym, scoor, types, hmat, delta, iounit)
          CALL cp_warn(__LOCATION__, "Symmetry library SPGLIB not available")
          spglib = .FALSE.
       ELSE
-         IF (major /= 1 .OR. minor /= 9 .OR. micro /= 9) THEN
-            CALL cp_warn(__LOCATION__, "Version of Symmetry Library SPGLIB not tested")
-         END IF
          spglib = .TRUE.
          CALL cite_reference(Togo2018)
          ierr = spg_get_international(csym%international_symbol, TRANSPOSE(hmat), scoor, types, nat, delta)

src/force_env_methods.F

@@ -1204,8 +1204,8 @@ SUBROUTINE mixed_energy_forces(force_env, calculate_forces)
          dip_exists = test_for_result(results=results_mix, description=description)
          IF (dip_exists) THEN
             CALL get_results(results=results_mix, description=description, values=dip_mix)
-            WRITE (unit_nr, '(/,1X,A,T48,3F11.6)') "MIXED ENV| DIPOLE  ( A.U.)|", dip_mix
-            WRITE (unit_nr, '(  1X,A,T48,3F11.6)') "MIXED ENV| DIPOLE  (Debye)|", dip_mix*debye
+            WRITE (unit_nr, '(/,1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE  ( A.U.)|", dip_mix
+            WRITE (unit_nr, '(  1X,A,T48,3F21.16)') "MIXED ENV| DIPOLE  (Debye)|", dip_mix*debye
          ELSE
             WRITE (unit_nr, *) "NO FORCE_EVAL section calculated the dipole"
          END IF

src/motion/bfgs_optimizer.F

@@ -7,6 +7,9 @@
 
 ! **************************************************************************************************
 !> \brief Routines for Geometry optimization using BFGS algorithm
+!> \par History
+!>      Module modified by Pierre-André Cazade [pcazade] 01.2020 - University of Limerick.
+!>      Modifications enable Space Group Symmetry.
 ! **************************************************************************************************
 MODULE bfgs_optimizer
 
@@ -72,6 +75,11 @@ MODULE bfgs_optimizer
    USE message_passing,                 ONLY: mp_min,&
                                               mp_sum
    USE particle_list_types,             ONLY: particle_list_type
+   USE space_groups,                    ONLY: identify_space_group,&
+                                              print_spgr,&
+                                              spgr_apply_rotations_coord,&
+                                              spgr_apply_rotations_force
+   USE space_groups_types,              ONLY: spgr_type
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -93,6 +101,8 @@ MODULE bfgs_optimizer
 !> \param geo_section ...
 !> \param gopt_env ...
 !> \param x0 ...
+!> \par History
+!>      01.2020 modified to perform Space Group Symmetry [pcazade]
 ! **************************************************************************************************
    RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, gopt_env, x0)
 
@@ -125,11 +135,13 @@ RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, go
       TYPE(cp_para_env_type), POINTER                    :: para_env
       TYPE(cp_subsys_type), POINTER                      :: subsys
       TYPE(section_vals_type), POINTER                   :: print_key, root_section
+      TYPE(spgr_type), POINTER                           :: spgr
 
-      NULLIFY (logger, g, blacs_env)
+      NULLIFY (logger, g, blacs_env, spgr)
       logger => cp_get_default_logger()
       para_env => force_env%para_env
       root_section => force_env%root_section
+      spgr => gopt_env%spgr
       t_old = m_walltime()
 
       CALL timeset(routineN, handle)
@@ -260,6 +272,14 @@ RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, go
       ALLOCATE (dr(ndf))
       dr(:) = 0.0_dp
 
+      ! find space_group
+      CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
+      IF (spgr%keep_space_group) THEN
+         CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
+         CALL spgr_apply_rotations_coord(spgr, x0)
+         CALL print_spgr(spgr)
+      END IF
+
       ! Geometry optimization starts now
       CALL cp_iterate(logger%iter_info, increment=0, iter_nr_out=iter_nr)
       CALL print_geo_opt_header(gopt_env, output_unit, wildcard)
@@ -268,6 +288,12 @@ RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, go
       CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
                       .FALSE., gopt_env%force_env%para_env)
 
+      ! Symmetrize coordinates and forces
+      IF (spgr%keep_space_group) THEN
+         CALL spgr_apply_rotations_coord(spgr, x0)
+         CALL spgr_apply_rotations_force(spgr, g)
+      END IF
+
       ! Print info at time 0
       emin = etot
       t_now = m_walltime()
@@ -290,20 +316,38 @@ RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, go
             IF (ionode) CALL close_file(unit_number=hesunit_read)
          ELSE
             IF ((its - iter_nr) > 1) THEN
+               ! Symmetrize old coordinates and old forces
+               IF (spgr%keep_space_group) THEN
+                  CALL spgr_apply_rotations_coord(spgr, xold)
+                  CALL spgr_apply_rotations_force(spgr, gold)
+               END IF
+
                DO indf = 1, ndf
                   dx(indf) = x0(indf) - xold(indf)
                   dg(indf) = g(indf) - gold(indf)
                END DO
 
                CALL bfgs(ndf, dx, dg, hess_mat, work, para_env)
 
+               ! Symmetrize coordinates and forces change
+               IF (spgr%keep_space_group) THEN
+                  CALL spgr_apply_rotations_force(spgr, dx)
+                  CALL spgr_apply_rotations_force(spgr, dg)
+               END IF
+
                !Possibly dump the Hessian file
                IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
                   CALL write_bfgs_hessian(geo_section, hess_mat, logger)
                ENDIF
             ENDIF
          END IF
 
+         ! Symmetrize coordinates and forces
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_coord(spgr, x0)
+            CALL spgr_apply_rotations_force(spgr, g)
+         END IF
+
          ! Setting the present positions & gradients as old
          xold(:) = x0
          gold(:) = g
@@ -328,19 +372,36 @@ RECURSIVE SUBROUTINE geoopt_bfgs(force_env, gopt_param, globenv, geo_section, go
             CALL rat_fun_opt(ndf, dg, eigval, work, eigvec_mat, g, para_env)
          END IF
          CALL geoopt_get_step(ndf, eigval, eigvec_mat, hess_tmp, dr, g, para_env, use_rfo)
+
+         ! Symmetrize dr
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_force(spgr, dr)
+         END IF
+
          CALL trust_radius(ndf, step, rad, rat, dr, output_unit)
 
          ! Update the atomic positions
          x0 = x0 + dr
 
+         ! Symmetrize coordinates
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_coord(spgr, x0)
+         END IF
+
          CALL energy_predict(ndf, work, hess_mat, dr, g, conv, pred, para_env)
          eold = etot
 
          ! Energy & Gradients at new step
          CALL cp_eval_at(gopt_env, x0, etot, g, gopt_env%force_env%para_env%mepos, &
                          .FALSE., gopt_env%force_env%para_env)
+
          ediff = etot - eold
 
+         ! Symmetrize forces
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_force(spgr, g)
+         END IF
+
          ! check for an external exit command
          CALL external_control(should_stop, "GEO", globenv=globenv)
          IF (should_stop) EXIT

src/motion/cg_optimizer.F

@@ -21,6 +21,7 @@ MODULE cg_optimizer
    USE cp_output_handling,              ONLY: cp_iterate,&
                                               cp_print_key_finished_output,&
                                               cp_print_key_unit_nr
+   USE cp_subsys_types,                 ONLY: cp_subsys_type
    USE force_env_types,                 ONLY: force_env_get,&
                                               force_env_type
    USE global_types,                    ONLY: global_environment_type
@@ -32,10 +33,22 @@ MODULE cg_optimizer
                                               print_geo_opt_nc
    USE gopt_f_types,                    ONLY: gopt_f_type
    USE gopt_param_types,                ONLY: gopt_param_type
+   USE input_constants,                 ONLY: default_cell_direct_id,&
+                                              default_cell_geo_opt_id,&
+                                              default_cell_md_id,&
+                                              default_cell_method_id,&
+                                              default_minimization_method_id,&
+                                              default_ts_method_id
    USE input_section_types,             ONLY: section_vals_type,&
+                                              section_vals_val_get,&
                                               section_vals_val_set
    USE kinds,                           ONLY: dp
    USE machine,                         ONLY: m_walltime
+   USE space_groups,                    ONLY: identify_space_group,&
+                                              print_spgr,&
+                                              spgr_apply_rotations_coord,&
+                                              spgr_apply_rotations_force
+   USE space_groups_types,              ONLY: spgr_type
 #include "../base/base_uses.f90"
 
    IMPLICIT NONE
@@ -77,16 +90,51 @@ RECURSIVE SUBROUTINE geoopt_cg(force_env, gopt_param, globenv, geo_section, &
       INTEGER                                            :: handle, output_unit
       LOGICAL                                            :: my_do_update
       TYPE(cp_logger_type), POINTER                      :: logger
+      TYPE(cp_subsys_type), POINTER                      :: subsys
+      TYPE(spgr_type), POINTER                           :: spgr
 
       CALL timeset(routineN, handle)
 
+      NULLIFY (spgr)
       logger => cp_get_default_logger()
+      spgr => gopt_env%spgr
 
       output_unit = cp_print_key_unit_nr(logger, geo_section, "PRINT%PROGRAM_RUN_INFO", &
                                          extension=".geoLog")
       CALL print_geo_opt_header(gopt_env, output_unit, "CONJUGATE GRADIENTS")
+
+      ! find space_group
+      CALL force_env_get(force_env, subsys=subsys)
+      CALL section_vals_val_get(geo_section, "KEEP_SPACE_GROUP", l_val=spgr%keep_space_group)
+      IF (spgr%keep_space_group) THEN
+         SELECT CASE (gopt_env%type_id)
+         CASE (default_minimization_method_id, default_ts_method_id)
+            CALL force_env_get(force_env, subsys=subsys)
+            CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
+            CALL spgr_apply_rotations_coord(spgr, x0)
+            CALL print_spgr(spgr)
+         CASE (default_cell_method_id)
+            SELECT CASE (gopt_env%cell_method_id)
+            CASE (default_cell_direct_id)
+               CALL force_env_get(force_env, subsys=subsys)
+               CALL identify_space_group(subsys, geo_section, gopt_env, output_unit)
+               CALL spgr_apply_rotations_coord(spgr, x0)
+               CALL print_spgr(spgr)
+            CASE (default_cell_geo_opt_id)
+               spgr%keep_space_group = .FALSE.
+            CASE (default_cell_md_id)
+               CPABORT("KEEP_SPACE_GROUP not implemented for motion method MD.")
+            CASE DEFAULT
+               spgr%keep_space_group = .FALSE.
+            END SELECT
+         CASE DEFAULT
+            spgr%keep_space_group = .FALSE.
+         END SELECT
+      END IF
+
       CALL cp_cg_main(force_env, x0, gopt_param, output_unit, globenv, &
                       gopt_env, do_update=my_do_update)
+
       CALL cp_print_key_finished_output(output_unit, logger, geo_section, &
                                         "PRINT%PROGRAM_RUN_INFO")
       IF (PRESENT(do_update)) do_update = my_do_update
@@ -132,11 +180,13 @@ RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv,
       REAL(KIND=dp), DIMENSION(:), POINTER               :: g, h, xi
       TYPE(cell_type), POINTER                           :: cell
       TYPE(cp_logger_type), POINTER                      :: logger
+      TYPE(cp_subsys_type), POINTER                      :: subsys
       TYPE(section_vals_type), POINTER                   :: root_section
+      TYPE(spgr_type), POINTER                           :: spgr
 
       CALL timeset(routineN, handle)
       t_old = m_walltime()
-      NULLIFY (logger, g, h, xi)
+      NULLIFY (logger, g, h, xi, spgr)
       root_section => force_env%root_section
       logger => cp_get_default_logger()
       conv = .FALSE.
@@ -148,16 +198,30 @@ RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv,
       ALLOCATE (h(SIZE(x0)))
       ALLOCATE (xi(SIZE(x0)))
       ALLOCATE (xold(SIZE(x0)))
-      CALL force_env_get(force_env, cell=cell)
+      CALL force_env_get(force_env, cell=cell, subsys=subsys)
+
+      spgr => gopt_env%spgr
+      ! applies rotation matrices to coordinates
+      IF (spgr%keep_space_group) THEN
+         CALL spgr_apply_rotations_coord(spgr, x0)
+      END IF
 
       ! Evaluate energy and forces at the first step
       ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
       save_consistent_energy_force = gopt_env%require_consistent_energy_force
       gopt_env%require_consistent_energy_force = .FALSE.
+
       CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
                       .FALSE., gopt_env%force_env%para_env)
+
       gopt_env%require_consistent_energy_force = save_consistent_energy_force
 
+      ! Symmetrize coordinates and forces
+      IF (spgr%keep_space_group) THEN
+         CALL spgr_apply_rotations_coord(spgr, x0)
+         CALL spgr_apply_rotations_force(spgr, xi)
+      END IF
+
       g = -xi
       h = g
       xi = h
@@ -175,11 +239,23 @@ RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv,
          CALL section_vals_val_set(gopt_env%geo_section, "STEP_START_VAL", i_val=its)
          CALL gopt_f_ii(its, output_unit)
 
+         ! Symmetrize coordinates and forces
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_coord(spgr, x0)
+            CALL spgr_apply_rotations_force(spgr, g)
+            CALL spgr_apply_rotations_force(spgr, xi)
+         END IF
+
          xold(:) = x0
 
          ! Line minimization
          CALL cg_linmin(gopt_env, x0, xi, g, opt_energy, output_unit, gopt_param, globenv)
 
+         ! Applies rotation matrices to coordinates
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_coord(spgr, x0)
+         END IF
+
          ! Check for an external exit command
          CALL external_control(should_stop, "GEO", globenv=globenv)
          IF (should_stop) EXIT
@@ -198,14 +274,27 @@ RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv,
          ![NB] consistent energies and forces not required for CG, but some line minimizers might set it
          save_consistent_energy_force = gopt_env%require_consistent_energy_force
          gopt_env%require_consistent_energy_force = .FALSE.
+
          CALL cp_eval_at(gopt_env, x0, opt_energy, xi, gopt_env%force_env%para_env%mepos, &
                          .FALSE., gopt_env%force_env%para_env)
+
          gopt_env%require_consistent_energy_force = save_consistent_energy_force
 
+         ! Symmetrize coordinates and forces
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_force(spgr, xi)
+         END IF
+
          ! Get Conjugate Directions:  updates the searching direction (h)
          wildcard = "   CG"
          CALL get_conjugate_direction(gopt_env, Fletcher_Reeves, g, xi, h)
 
+         ! Symmetrize coordinates and forces
+         IF (spgr%keep_space_group) THEN
+            CALL spgr_apply_rotations_force(spgr, g)
+            CALL spgr_apply_rotations_force(spgr, h)
+         END IF
+
          ! Reset Condition or Steepest Descent Requested
          ! ABS(DOT_PRODUCT(g, h))/SQRT((DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h))) > res_lim ...
          IF ((DOT_PRODUCT(g, h)*DOT_PRODUCT(g, h)) > (res_lim*res_lim*DOT_PRODUCT(g, g)*DOT_PRODUCT(h, h)) &
@@ -217,6 +306,7 @@ RECURSIVE SUBROUTINE cp_cg_main(force_env, x0, gopt_param, output_unit, globenv,
          g = -xi
          xi = h
       END DO
+
       IF (its == maxiter .AND. (.NOT. conv)) THEN
          CALL print_geo_opt_nc(gopt_env, output_unit)
       END IF