<div dir="ltr">Hello,<br><br>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.<br>It can be invoked with "SYMMETRY MONOCLINIC_GAMMA_AB" in the CELL section of the input file.<br><br>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.<br><br>The entire patch is given here in terms of the printout of svn diff:<br><br>Index: cp2k/src/input_cp2k_subsys.F<br>===================================================================<br>--- cp2k/src/input_cp2k_subsys.F    (revision 18110)<br>+++ cp2k/src/input_cp2k_subsys.F    (working copy)<br>@@ -16,7 +16,7 @@<br>                                               Krack2005,&<br>                                               VandeVondele2005a,&<br>                                               VandeVondele2007<br>-   USE cell_types,                      ONLY: &<br>+   USE cell_types,                      ONLY: cell_sym_monoclinic_gamma_ab, &<br>         cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_none, &<br>         cell_sym_orthorhombic, cell_sym_rhombohedral, cell_sym_tetragonal_ab, &<br>         cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, cell_sym_triclinic, use_perd_none, &<br>@@ -202,12 +202,14 @@<br>                        "Tetragonal (alias for TETRAGONAL_AB)", &<br>                        "Rhombohedral (a = b = c, &alpha; = &beta; = &gamma; &ne; 90&deg;)", &<br>                        "Hexagonal (a = b &ne; c, &alpha; = &beta; = 90&deg;, &gamma; = 60&deg;)", &<br>-                       "Cubic (a = b = c, &alpha; = &beta; = &gamma; = 90&deg;)"), &<br>+                       "Cubic (a = b = c, &alpha; = &beta; = &gamma; = 90&deg;)", &<br>+                       "Monoclinic_gamma_ab (a = b &ne; c, &alpha; = &beta; = 90&deg;, &gamma; &ne; 90&deg;)"), &<br>          enum_c_vals=s2a("NONE", "TRICLINIC", "MONOCLINIC", "ORTHORHOMBIC", "TETRAGONAL_AB", "TETRAGONAL_AC", &<br>-                         "TETRAGONAL_BC", "TETRAGONAL", "RHOMBOHEDRAL", "HEXAGONAL", "CUBIC"), &<br>+                         "TETRAGONAL_BC", "TETRAGONAL", "RHOMBOHEDRAL", "HEXAGONAL", "CUBIC", "MONOCLINIC_GAMMA_AB"), &<br>          enum_i_vals=(/cell_sym_none, cell_sym_triclinic, cell_sym_monoclinic, cell_sym_orthorhombic, &<br>                        cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, cell_sym_tetragonal_bc, &<br>-                       cell_sym_tetragonal_ab, cell_sym_rhombohedral, cell_sym_hexagonal, cell_sym_cubic/), &<br>+                       cell_sym_tetragonal_ab, cell_sym_rhombohedral, cell_sym_hexagonal, cell_sym_cubic, &<br>+                       cell_sym_monoclinic_gamma_ab/), &<br>          default_i_val=cell_sym_none)<br>       CALL section_add_keyword(section, keyword)<br>       CALL keyword_release(keyword)<br>Index: cp2k/src/motion/cell_opt_utils.F<br>===================================================================<br>--- cp2k/src/motion/cell_opt_utils.F    (revision 18110)<br>+++ cp2k/src/motion/cell_opt_utils.F    (working copy)<br>@@ -10,7 +10,7 @@<br> !>      03.2008 - Teodoro Laino [tlaino] - University of Zurich - Cell Optimization<br> ! **************************************************************************************************<br> MODULE cell_opt_utils<br>-   USE cell_types,                      ONLY: &<br>+   USE cell_types,                      ONLY: cell_sym_monoclinic_gamma_ab, &<br>         cell_sym_cubic, cell_sym_hexagonal, cell_sym_monoclinic, cell_sym_orthorhombic, &<br>         cell_sym_rhombohedral, cell_sym_tetragonal_ab, cell_sym_tetragonal_ac, &<br>         cell_sym_tetragonal_bc, cell_sym_triclinic, cell_type, get_cell_param, set_cell_param<br>@@ -337,7 +337,8 @@<br>          routineP = moduleN//':'//routineN<br> <br>       REAL(KIND=dp)                                      :: a, cosa, cosah, g, norm, norm_b, norm_c, &<br>-                                                            sina, sinah<br>+                                                            sina, sinah, &<br>+                                                            a_length, b_length, ab_length, gamma, cosgamma, singamma, deriv_gamma<br> <br>       IF (keep_angles) THEN<br>          ! If we want to keep the angles constant we have to project out the<br>@@ -421,6 +422,26 @@<br>          CASE (cell_sym_monoclinic)<br>             gradient(2) = 0.0_dp<br>             gradient(5) = 0.0_dp<br>+! Case cell_sym_monoclinic_gamma_ab for cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg<br>+         CASE (cell_sym_monoclinic_gamma_ab)<br>+            a_length = SQRT(cell%hmat(1, 1)*cell%hmat(1, 1)+ &<br>+            cell%hmat(2, 1)*cell%hmat(2, 1)+ &<br>+            cell%hmat(3, 1)*cell%hmat(3, 1))<br>+            b_length = SQRT(cell%hmat(1, 2)*cell%hmat(1, 2)+ &<br>+            cell%hmat(2, 2)*cell%hmat(2, 2)+ &<br>+            cell%hmat(3, 2)*cell%hmat(3, 2))<br>+            ab_length = 0.5_dp*(a_length+b_length)<br>+            gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))<br>+            cosgamma = COS(gamma)<br>+            singamma = SIN(gamma)<br>+! Here, g is the average derivative of the cell vector length ab_length, and deriv_gamma is the derivative of the angle gamma<br>+            g = 0.5_dp*(gradient(1)+cosgamma*gradient(2)+singamma*gradient(3))<br>+            deriv_gamma = (gradient(3)*cosgamma-gradient(2)*singamma)/b_length<br>+            gradient(1) = g<br>+            gradient(2) = g*cosgamma-ab_length*singamma*deriv_gamma<br>+            gradient(3) = g*singamma+ab_length*cosgamma*deriv_gamma<br>+            gradient(4) = 0.0_dp<br>+            gradient(5) = 0.0_dp<br>          CASE (cell_sym_triclinic)<br>             ! Nothing to do<br>          END SELECT<br>Index: cp2k/src/subsys/cell_types.F<br>===================================================================<br>--- cp2k/src/subsys/cell_types.F    (revision 18110)<br>+++ cp2k/src/subsys/cell_types.F    (working copy)<br>@@ -39,7 +39,8 @@<br>                                                cell_sym_tetragonal_bc = 6, &<br>                                                cell_sym_rhombohedral = 7, &<br>                                                cell_sym_hexagonal = 8, &<br>-                                               cell_sym_cubic = 9<br>+                                               cell_sym_cubic = 9, &<br>+                                               cell_sym_monoclinic_gamma_ab = 10<br> <br>    INTEGER, PARAMETER, PUBLIC               :: use_perd_x = 0, &<br>                                                use_perd_y = 1, &<br>@@ -343,7 +344,8 @@<br> <br>       INTEGER                                            :: dim<br>       REAL(KIND=dp)                                      :: a, acosa, acosah, alpha, asina, asinah, &<br>-                                                            beta, norm, norm_c<br>+                                                            beta, norm, norm_c, &<br>+                                                            gamma, acosgamma, asingamma<br>       REAL(KIND=dp), DIMENSION(3)                        :: abc<br> <br>       IF (PRESENT(hmat)) cell%hmat(:, :) = hmat(:, :)<br>@@ -356,7 +358,6 @@<br>                        "An invalid set of cell vectors was specified. "// &<br>                        "The determinant det(h) is too small")<br>       END IF<br>-      cell%h_inv = inv_3x3(cell%hmat)<br> <br>       SELECT CASE (cell%symmetry_id)<br>       CASE (cell_sym_cubic, &<br>@@ -420,6 +421,16 @@<br>          cell%hmat(1, 1) = abc(1); cell%hmat(1, 2) = 0.0_dp; cell%hmat(1, 3) = abc(3)*COS(beta)<br>          cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = abc(2); cell%hmat(2, 3) = 0.0_dp<br>          cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)*SIN(beta)<br>+! Case cell_sym_monoclinic_gamma_ab for cell symmetry with a=b, alpha=beta=90deg and gammma unequal 90deg<br>+      CASE (cell_sym_monoclinic_gamma_ab)<br>+         CALL get_cell(cell=cell, abc=abc)<br>+         a = 0.5_dp*(abc(1)+abc(2))<br>+         gamma = angle(cell%hmat(:, 1), cell%hmat(:, 2))<br>+         acosgamma = a*COS(gamma)<br>+         asingamma = a*SIN(gamma)<br>+         cell%hmat(1, 1) = a; cell%hmat(1, 2) = acosgamma; cell%hmat(1, 3) = 0.0_dp<br>+         cell%hmat(2, 1) = 0.0_dp; cell%hmat(2, 2) = asingamma; cell%hmat(2, 3) = 0.0_dp<br>+         cell%hmat(3, 1) = 0.0_dp; cell%hmat(3, 2) = 0.0_dp; cell%hmat(3, 3) = abc(3)<br>       CASE (cell_sym_triclinic)<br>          ! Nothing to do<br>       END SELECT<br>@@ -449,6 +460,10 @@<br>          CPABORT("Non-orthorhombic and not periodic")<br>       END IF<br> <br>+! Update deth and hmat_inv with enforced symmetry<br>+      cell%deth = ABS(det_3x3(cell%hmat))<br>+      cell%h_inv = inv_3x3(cell%hmat)<br>+<br>    END SUBROUTINE init_cell<br> <br> ! **************************************************************************************************<br><br><br>A regtest input file for the new symmetry and the corresponding calculation output file are attached.<br><br>Best regards,<br>Tobias<br><br><br><br></div>