[CP2K:9634] Additional cell symmetry and bug fixed
Krack Matthias (PSI)
matthia... at psi.ch
Tue Nov 7 07:37:33 UTC 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.htm>
More information about the CP2K-user
mailing list