Additional cell symmetry and bug fixed
Tobias Binninger
t... at zurich.ibm.com
Mon Nov 6 16:41:18 UTC 2017
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
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20171106/5c70c411/attachment.htm>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: cell_sym_monoclinic_gamma_ab.inp
Type: chemical/x-gamess-input
Size: 8134 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20171106/5c70c411/attachment.inp>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: output_cell_sym_monoclinic_gamma_ab
Type: application/octet-stream
Size: 292537 bytes
Desc: not available
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20171106/5c70c411/attachment.obj>
More information about the CP2K-user
mailing list