<html xmlns:v="urn:schemas-microsoft-com:vml" xmlns:o="urn:schemas-microsoft-com:office:office" xmlns:w="urn:schemas-microsoft-com:office:word" xmlns:m="http://schemas.microsoft.com/office/2004/12/omml" xmlns="http://www.w3.org/TR/REC-html40"> <head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8"> <meta name="Generator" content="Microsoft Word 14 (filtered medium)"> <style></style> </head> <body lang="EN-GB" link="blue" vlink="purple"> <div class="WordSection1"> Hi Tobias<o:p></o:p> <o:p> </o:p> Thanks for the bug fix and the additional cell symmetry. I will commit it after checking.<o:p></o:p> <o:p> </o:p> Best regards<o:p></o:p> <o:p> </o:p> Matthias<o:p></o:p> <o:p> </o:p> From: cp...@googlegroups.com [mailto:cp...@googlegroups.com] On Behalf Of Tobias Binninger Sent: 06 November 2017 17:41 To: cp2k Subject: [CP2K:9634] Additional cell symmetry and bug fixed<o:p></o:p> <o:p> </o:p> <div> 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 <o:p></o:p> </div> -- 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 <a href="mailto:cp2k+unsu...@googlegroups.com">cp2k+unsu...@googlegroups.com</a>. To post to this group, send email to <a href="mailto:cp...@googlegroups.com">cp...@googlegroups.com</a>. Visit this group at <a href="https://groups.google.com/group/cp2k">https://groups.google.com/group/cp2k</a>. For more options, visit <a href="https://groups.google.com/d/optout">https://groups.google.com/d/optout</a>.<o:p></o:p> </div> </body> </html>