[CP2K:9634] Additional cell symmetry and bug fixed

Krack Matthias (PSI) matthia... at psi.ch
Tue Nov 7 08:37:33 CET 2017


Hi Tobias

Thanks for the bug fix and the additional cell symmetry. I will commit it after checking.

Best regards

Matthias

From: cp... at googlegroups.com [mailto:cp... at googlegroups.com] On Behalf Of Tobias Binninger
Sent: 06 November 2017 17:41
To: cp2k
Subject: [CP2K:9634] Additional cell symmetry and bug fixed

Hello,

I implemented an additional cell symmetry in the cp2k_trunk version. It is characterized by a=b unequal c and alpha=beta=90deg unequal gamma.
It can be invoked with "SYMMETRY MONOCLINIC_GAMMA_AB" in the CELL section of the input file.

In the course of the implementation, I also fixed a bug in the "init_cell" subroutine in "cell_types.F". Here, the determinant cell%deth and the inverse cell%h_inv of cell%hmat had been computed before enforcing the cell symmetry on cell%hmat with the consequence that, in case the symmetry enforcement actually did change the cell h-matrix, cell%deth and cell%h_inv still corresponded to the unchanged h-matrix. Now, cell%deth and cell%h_inv are recomputed at the end of the subroutine so that they always correspond to the updated cell%hmat.

The entire patch is given here in terms of the printout of svn diff:

Index: cp2k/src/input_cp2k_subsys.F
===================================================================
--- cp2k/src/input_cp2k_subsys.F    (revision 18110)
+++ cp2k/src/input_cp2k_subsys.F    (working copy)
@@ -16,7 +16,7 @@
                                               Krack2005,&
                                               VandeVondele2005a,&
                                               VandeVondele2007
-   USE cell_types,                      ONLY: &
+   USE cell_types,                      ONLY: cell_sym_monoclinic_gamma_ab, &
         cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_none, &
         cell_sym_orthorhombic, cell_sym_rhombohedral, cell_sym_tetragonal_ab, &
         cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, cell_sym_triclinic, use_perd_none, &
@@ -202,12 +202,14 @@
                        "Tetragonal (alias for TETRAGONAL_AB)", &
                        "Rhombohedral (a = b = c, α = β = γ ≠ 90°)", &
                        "Hexagonal (a = b ≠ c, α = β = 90°, γ = 60°)", &
-                       "Cubic (a = b = c, α = β = γ = 90°)"), &
+                       "Cubic (a = b = c, α = β = γ = 90°)", &
+                       "Monoclinic_gamma_ab (a = b ≠ c, α = β = 90°, γ ≠ 90°)"), &
          enum_c_vals=s2a("NONE", "TRICLINIC", "MONOCLINIC", "ORTHORHOMBIC", "TETRAGONAL_AB", "TETRAGONAL_AC", &
-                         "TETRAGONAL_BC", "TETRAGONAL", "RHOMBOHEDRAL", "HEXAGONAL", "CUBIC"), &
+                         "TETRAGONAL_BC", "TETRAGONAL", "RHOMBOHEDRAL", "HEXAGONAL", "CUBIC", "MONOCLINIC_GAMMA_AB"), &
          enum_i_vals=(/cell_sym_none, cell_sym_triclinic, cell_sym_monoclinic, cell_sym_orthorhombic, &
                        cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, &
-                       cell_sym_tetragonal_ab, cell_sym_rhombohedral, cell_sym_hexagonal, cell_sym_cubic/), &
+                       cell_sym_tetragonal_ab, cell_sym_rhombohedral, cell_sym_hexagonal, cell_sym_cubic, &
+                       cell_sym_monoclinic_gamma_ab/), &
          default_i_val=cell_sym_none)
       CALL section_add_keyword(section, keyword)
       CALL keyword_release(keyword)
Index: cp2k/src/motion/cell_opt_utils.F
===================================================================
--- cp2k/src/motion/cell_opt_utils.F    (revision 18110)
+++ cp2k/src/motion/cell_opt_utils.F    (working copy)
@@ -10,7 +10,7 @@
 !>      03.2008 - Teodoro Laino [tlaino] - University of Zurich - Cell Optimization
 ! **************************************************************************************************
 MODULE cell_opt_utils
-   USE cell_types,                      ONLY: &
+   USE cell_types,                      ONLY: cell_sym_monoclinic_gamma_ab, &
         cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_orthorhombic, &
         cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &
         cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type, get_cell_param, set_cell_param
@@ -337,7 +337,8 @@
          routineP = moduleN//':'//routineN

       REAL(KIND=dp)                                      :: a, cosa, cosah, g, norm, norm_b, norm_c, &
-                                                            sina, sinah
+                                                            sina, sinah, &
+                                                            a_length, b_length, ab_length, gamma, cosgamma, singamma, deriv_gamma

       IF (keep_angles) THEN
          ! If we want to keep the angles constant we have to project out the
@@ -421,6 +422,26 @@
          CASE (cell_sym_monoclinic)
             gradient(2) = 0.0_dp
             gradient(5) = 0.0_dp
+! Case cell_sym_monoclinic_gamma_ab for cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg
+         CASE (cell_sym_monoclinic_gamma_ab)
+            a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1)+ &
+            cell%hmat(2, 1)*cell%hmat(2, 1)+ &
+            cell%hmat(3, 1)*cell%hmat(3, 1))
+            b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2)+ &
+            cell%hmat(2, 2)*cell%hmat(2, 2)+ &
+            cell%hmat(3, 2)*cell%hmat(3, 2))
+            ab_length = 0.5_dp*(a_length+b_length)
+            gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
+            cosgamma = COS(gamma)
+            singamma = SIN(gamma)
+! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma
+            g = 0.5_dp*(gradient(1)+cosgamma*gradient(2)+singamma*gradient(3))
+            deriv_gamma = (gradient(3)*cosgamma-gradient(2)*singamma)/b_length
+            gradient(1) = g
+            gradient(2) = g*cosgamma-ab_length*singamma*deriv_gamma
+            gradient(3) = g*singamma+ab_length*cosgamma*deriv_gamma
+            gradient(4) = 0.0_dp
+            gradient(5) = 0.0_dp
          CASE (cell_sym_triclinic)
             ! Nothing to do
          END SELECT
Index: cp2k/src/subsys/cell_types.F
===================================================================
--- cp2k/src/subsys/cell_types.F    (revision 18110)
+++ cp2k/src/subsys/cell_types.F    (working copy)
@@ -39,7 +39,8 @@
                                                cell_sym_tetragonal_bc = 6, &
                                                cell_sym_rhombohedral = 7, &
                                                cell_sym_hexagonal = 8, &
-                                               cell_sym_cubic = 9
+                                               cell_sym_cubic = 9, &
+                                               cell_sym_monoclinic_gamma_ab = 10

    INTEGER, PARAMETER, PUBLIC               :: use_perd_x = 0, &
                                                use_perd_y = 1, &
@@ -343,7 +344,8 @@

       INTEGER                                            :: dim
       REAL(KIND=dp)                                      :: a, acosa, acosah, alpha, asina, asinah, &
-                                                            beta, norm, norm_c
+                                                            beta, norm, norm_c, &
+                                                            gamma, acosgamma, asingamma
       REAL(KIND=dp), DIMENSION(3)                        :: abc

       IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)
@@ -356,7 +358,6 @@
                        "An invalid set of cell vectors was specified. "// &
                        "The determinant det(h) is too small")
       END IF
-      cell%h_inv = inv_3x3(cell%hmat)

       SELECT CASE (cell%symmetry_id)
       CASE (cell_sym_cubic, &
@@ -420,6 +421,16 @@
          cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)
          cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp
          cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)
+! Case cell_sym_monoclinic_gamma_ab for cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg
+      CASE (cell_sym_monoclinic_gamma_ab)
+         CALL get_cell(cell=cell, abc=abc)
+         a = 0.5_dp*(abc(1)+abc(2))
+         gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))
+         acosgamma = a*COS(gamma)
+         asingamma = a*SIN(gamma)
+         cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosgamma; cell%hmat(1, 3) = 0.0_dp
+         cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asingamma; cell%hmat(2, 3) = 0.0_dp
+         cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)
       CASE (cell_sym_triclinic)
          ! Nothing to do
       END SELECT
@@ -449,6 +460,10 @@
          CPABORT("Non-orthorhombic and not periodic")
       END IF

+! Update deth and hmat_inv with enforced symmetry
+      cell%deth = ABS(det_3x3(cell%hmat))
+      cell%h_inv = inv_3x3(cell%hmat)
+
    END SUBROUTINE init_cell

 ! **************************************************************************************************


A regtest input file for the new symmetry and the corresponding calculation output file are attached.

Best regards,
Tobias


--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns... at googlegroups.com<mailto:cp2k+uns... at googlegroups.com>.
To post to this group, send email to cp... at googlegroups.com<mailto:cp... at googlegroups.com>.
Visit this group at https://groups.google.com/group/cp2k.
For more options, visit https://groups.google.com/d/optout.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20171107/3b22322b/attachment.html>


More information about the CP2K-user mailing list