Branch data Line data Source code
1 : : !!****m* ABINIT/m_tdep_latt
2 : : !! NAME
3 : : !! m_tdep_latt
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the TDEP Lattice data type, which describes the
7 : : !! lattice of the unitcell and of the supercell, and the routines building it.
8 : : !!
9 : : !! COPYRIGHT
10 : : !! Copyright (C) 2011-2026 ABINIT group (GA,FB,JB)
11 : : !! This file is distributed under the terms of the
12 : : !! GNU General Public License, see ~abinit/COPYING
13 : : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : : !!
15 : : !! SOURCE
16 : :
17 : : #if defined HAVE_CONFIG_H
18 : : #include "config.h"
19 : : #endif
20 : :
21 : : #include "abi_common.h"
22 : :
23 : : module m_tdep_latt
24 : :
25 : : use defs_basis
26 : : use m_abicore
27 : : use m_errors
28 : : use m_xmpi
29 : :
30 : : use m_matrix, only : matr3inv, mat33det
31 : : use m_geometry, only : metric
32 : : use m_crystal, only : crystal_t
33 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
34 : :
35 : : implicit none
36 : : !!***
37 : :
38 : : !!****t* m_tdep_latt/Lattice_type
39 : : !! NAME
40 : : !! Lattice_type
41 : : !!
42 : : !! FUNCTION
43 : : !! Holds the primitive vectors of the unitcell and of the supercell
44 : : !! (in various forms: transposed, inverted, dimensioned or not), the
45 : : !! lattice parameters, the Bravais lattice and the supercell multiplicity.
46 : : !!
47 : : !! SOURCE
48 : :
49 : : type Lattice_type
50 : :
51 : : double precision :: acell_unitcell(3)
52 : : ! Lattice parameters of the unitcell (bohr)
53 : :
54 : : double precision :: angle_alpha
55 : : ! Angle (in degrees) defining the lattice, for some Bravais lattices
56 : :
57 : : integer :: bravais(11)
58 : : ! Bravais lattice of the lattice alone (bravais(1): holohedry, bravais(2): centering), see symlatt
59 : :
60 : : integer :: line
61 : : ! How the lattice parameters multiply rprim (see tdep_make_latt):
62 : : ! 0=line or column, 1=line only, 2=column only, 3=neither
63 : :
64 : : double precision :: metmin (3,3)
65 : : ! Currently unused
66 : :
67 : : double precision :: minim (3,3)
68 : : ! Currently unused
69 : :
70 : : double precision :: multiplicity (3,3)
71 : : ! Supercell scaling matrix
72 : :
73 : : double precision :: multiplicitym1(3,3)
74 : : ! Inverse of multiplicity
75 : :
76 : : double precision :: gmet (3,3)
77 : : ! Reciprocal space metric of the unitcell (bohr^-2), computed from rprimdt
78 : :
79 : : double precision :: rmet (3,3)
80 : : ! Real space metric of the unitcell (bohr^2), computed from rprimdt
81 : :
82 : : double precision :: gprimd (3,3)
83 : : ! Dimensioned reciprocal primitive vectors of the unitcell (bohr^-1),
84 : : ! computed from rprimdt
85 : :
86 : : double precision :: gprim (3,3)
87 : : ! Inverse of rprim (dimensionless reciprocal primitive vectors of the unitcell)
88 : :
89 : : double precision :: gprimt (3,3)
90 : : ! Transpose of gprim
91 : :
92 : : double precision :: rprim (3,3)
93 : : ! Dimensionless primitive vectors of the unitcell
94 : :
95 : : double precision :: rprimt (3,3)
96 : : ! Transpose of rprim
97 : :
98 : : double precision :: rprimm1 (3,3)
99 : : ! Inverse of rprim
100 : :
101 : : double precision :: rprimd (3,3)
102 : : ! Dimensioned primitive vectors of the unitcell (bohr)
103 : :
104 : : double precision :: rprimdm1 (3,3)
105 : : ! Inverse of rprimd
106 : :
107 : : double precision :: rprimdt (3,3)
108 : : ! Transpose of rprimd
109 : :
110 : : double precision :: rprimdtm1 (3,3)
111 : : ! Inverse of rprimdt
112 : :
113 : : double precision :: rprimd_md (3,3)
114 : : ! Dimensioned primitive vectors of the supercell (bohr),
115 : : ! recomputed from rprim, acell_unitcell and multiplicity
116 : :
117 : : end type Lattice_type
118 : : !!***
119 : :
120 : : public :: tdep_make_inbox
121 : : public :: tdep_make_latt
122 : : public :: tdep_unitcell_rprimd
123 : :
124 : : contains
125 : :
126 : : !=====================================================================================================
127 : :
128 : : !!****f* ABINIT/m_tdep_latt/tdep_make_inbox
129 : : !! NAME
130 : : !! tdep_make_inbox
131 : : !!
132 : : !! FUNCTION
133 : : !! Wrap a set of vectors in reduced coordinates into [-0.5,0.5[, by removing
134 : : !! from each component the nearest integer (a small tolerance is used to
135 : : !! choose between the two boundaries). If temp is present, the integer shifts
136 : : !! are computed from tab but applied to temp, and tab is left unchanged.
137 : : !!
138 : : !! INPUTS
139 : : !! natom = Number of vectors.
140 : : !! tol = Tolerance used when rounding to the nearest integer.
141 : : !!
142 : : !! OUTPUT
143 : : !!
144 : : !! SIDE EFFECTS
145 : : !! tab(3,natom) = Vectors in reduced coordinates. Wrapped on exit,
146 : : !! if temp is not present.
147 : : !! temp(3,natom) = (optional) Vectors in reduced coordinates, shifted on exit
148 : : !! by the integer shifts computed from tab.
149 : : !!
150 : : !! SOURCE
151 : :
152 : 23145431 : subroutine tdep_make_inbox(tab,natom,tol,&
153 : : & temp) !Optional
154 : :
155 : : integer :: natom,ii,jj,iatom
156 : : double precision :: tol
157 : : double precision :: tab(3,natom) !Input (and ouput if temp is not present)
158 : : double precision,optional :: temp(3,natom) !Input/Output (if present)
159 : :
160 : : !GA: One could make use of function wrap2_pmhalf
161 : :
162 [ + + ]: 46290862 : do iatom=1,natom
163 [ + + ]: 115727155 : do ii=1,3
164 [ + + ]: 69436293 : if (tab(ii,iatom).lt.0.d0) then
165 : 21620643 : jj=int(tab(ii,iatom)-0.5+tol)
166 [ + - ]: 47815650 : else if (tab(ii,iatom).ge.0.d0) then
167 : 47815650 : jj=int(tab(ii,iatom)+0.5+tol)
168 : : end if
169 [ + + ]: 92581724 : if (present(temp)) then
170 : 69118575 : temp(ii,iatom)=temp(ii,iatom)-real(jj)
171 : : else
172 : 317718 : tab(ii,iatom)=tab(ii,iatom)-real(jj)
173 : : end if
174 : : end do
175 : : end do
176 : :
177 : 23145431 : end subroutine tdep_make_inbox
178 : : !!***
179 : :
180 : : !=====================================================================================================
181 : :
182 : : !!****f* ABINIT/m_tdep_latt/tdep_unitcell_rprimd
183 : : !! NAME
184 : : !! tdep_unitcell_rprimd
185 : : !!
186 : : !! FUNCTION
187 : : !! Return the primitive vectors of the unitcell given in input, deduced from
188 : : !! the primitive vectors of the supercell (Invar%rprimd_md) and the multiplicity:
189 : : !! rprimd_md = multiplicity * rprimd_unitcell
190 : : !! (with the vectors stored in lines, see tdep_make_latt).
191 : : !!
192 : : !! INPUTS
193 : : !! Invar = Input object containing the input variables.
194 : : !!
195 : : !! OUTPUT
196 : : !! rprimd(3,3) = Primitive vectors of the unitcell (bohr), in the ABINIT
197 : : !! convention (rprimd(:,ii) is the ii-th vector).
198 : : !!
199 : : !! SOURCE
200 : :
201 [ - + ]: 46 : function tdep_unitcell_rprimd(Invar) result(rprimd)
202 : :
203 : : type(atdep_dataset_type),intent(in) :: Invar
204 : : double precision :: rprimd(3,3)
205 : :
206 : : double precision :: multiplicitym1t(3,3)
207 : :
208 : 46 : call matr3inv(Invar%multiplicity, multiplicitym1t)
209 [ + + + + : 2438 : rprimd = MATMUL(TRANSPOSE(Invar%rprimd_md), multiplicitym1t)
+ + + + +
+ ]
210 : :
211 : 46 : end function tdep_unitcell_rprimd
212 : : !!***
213 : :
214 : : !=====================================================================================================
215 : :
216 : : !!****f* ABINIT/m_tdep_latt/tdep_make_latt
217 : : !! NAME
218 : : !! tdep_make_latt
219 : : !!
220 : : !! FUNCTION
221 : : !! Initialize the Lattice object. The dimensionless primitive vectors of
222 : : !! the unitcell (rprim) are built from the Bravais lattice of the lattice
223 : : !! alone (bravais_lattice, see Crystal%find_bravais) and, if needed, the angle
224 : : !! between the primitive vectors (Crystal%angdeg). The lattice parameters
225 : : !! (acell_unitcell) and the rotation matrix R are then deduced from the
226 : : !! primitive vectors of the supercell given in input (Invar%rprimd_md),
227 : : !! which are related to rprim by:
228 : : !! rprimd_md = R * A * multiplicity * rprim
229 : : !! where A is the diagonal matrix containing acell_unitcell.
230 : : !! The input supercell vectors are rotated accordingly, and checked against
231 : : !! the ones recomputed from rprim, acell_unitcell and multiplicity.
232 : : !! Finally, the dimensioned primitive vectors of the unitcell, their
233 : : !! transpose and inverse, and the metrics are computed.
234 : : !!
235 : : !! INPUTS
236 : : !! Crystal = Crystal object describing the unitcell given in input, before
237 : : !! rotation (see tdep_init_crystal).
238 : : !! bravais_lattice(11) = Bravais lattice of the lattice alone (holohedry, centering, ...),
239 : : !! see Crystal%find_bravais.
240 : : !!
241 : : !! OUTPUT
242 : : !! Lattice = Lattice object, initialized.
243 : : !! rotation_cart(3,3) = Rotation matrix, in cartesian coordinates, applied to
244 : : !! the input supercell primitive vectors (the cartesian
245 : : !! forces must be rotated accordingly, see
246 : : !! tdep_sampling_rotate).
247 : : !!
248 : : !! SIDE EFFECTS
249 : : !! Invar = Input object containing the input variables.
250 : : !! On exit, Invar%rprimd_md is rotated.
251 : : !!
252 : : !! SOURCE
253 : :
254 : 46 : subroutine tdep_make_latt(Invar,Crystal,bravais_lattice,Lattice,rotation_cart)
255 : :
256 : : type(Lattice_type),intent(out) :: Lattice
257 : : type(atdep_dataset_type),intent(inout) :: Invar
258 : : type(crystal_t),intent(in) :: Crystal
259 : : integer,intent(in) :: bravais_lattice(11)
260 : : double precision, intent(out) :: rotation_cart(3,3)
261 : :
262 : : integer :: ii,jj,line
263 : : integer :: bravais(11)
264 : : double precision :: angle_alpha
265 : : double precision :: acell_unitcell(3),multiplicity(3,3),multiplicitym1(3,3)
266 : : double precision :: rprimd(3,3),rprimdt(3,3),rprimd_md(3,3),rprimdm1(3,3)
267 : : double precision :: rprim(3,3),rprimm1(3,3),rprimt(3,3),rprimdtm1(3,3)
268 : : double precision :: rprimd_unitcell(3,3)
269 : : double precision :: rprimd_tmp(3,3), rprimdm1_tmp(3,3)
270 : : double precision :: RAmat(3,3),Amat2(3,3),Rmat(3,3)
271 : : double precision :: rotation(3,3)
272 : : double precision :: xi,hh
273 : : double precision :: ucvol
274 : : character(len=500) :: msg
275 : :
276 : : ! For bravais(1):
277 : : ! The holohedral groups are numbered as follows
278 : : ! (see international tables for crystallography (1983), p. 13)
279 : : ! iholohedry=1 triclinic 1bar
280 : : ! iholohedry=2 monoclinic 2/m
281 : : ! iholohedry=3 orthorhombic mmm
282 : : ! iholohedry=4 tetragonal 4/mmm
283 : : ! iholohedry=5 trigonal 3bar m
284 : : ! iholohedry=6 hexagonal 6/mmm
285 : : ! iholohedry=7 cubic m3bar m
286 : :
287 : : ! For bravais(2):
288 : : ! Centering
289 : : ! center=0 no centering
290 : : ! center=-1 body-centered
291 : : ! center=-3 face-centered
292 : : ! center=1 A-face centered
293 : : ! center=2 B-face centered
294 : : ! center=3 C-face centered
295 : :
296 : : ! The Bravais lattice of the lattice alone (the atoms may lower its symmetry, see
297 : : ! Sym%init), and the angle defining the primitive vectors of the monoclinic
298 : : ! (between the 1st and 3rd vectors) and trigonal (between any two vectors) lattices.
299 : 46 : bravais(:) = bravais_lattice(:)
300 : 46 : angle_alpha = 90.d0
301 [ + + ]: 46 : if (bravais(1).eq.2) angle_alpha = Crystal%angdeg(2)
302 [ + + ]: 46 : if (bravais(1).eq.5) angle_alpha = Crystal%angdeg(1)
303 : :
304 : : ! Initialize some local variables
305 : 46 : acell_unitcell(:)=zero; multiplicity(:,:)=zero; multiplicitym1(:,:)=zero
306 : 46 : rprimd(:,:)=zero; rprimdm1(:,:)=zero; rprimdt(:,:)=zero; rprimd_md(:,:)=zero
307 : 46 : rprim(:,:)=zero; rprimm1(:,:)=zero; rprimt(:,:)=zero; rprimdtm1(:,:)=zero
308 : 46 : rprimd_unitcell(:,:)=zero; RAmat(:,:)=zero; Amat2(:,:)=zero; Rmat(:,:)=zero
309 : :
310 : :
311 : : ! ---------------------------------------------------------------------------- !
312 : : ! Here, rprim defines a Primitive Lattice and NOT a Conventional lattice
313 : : ! The lattice parameters have to multiply rprim on:
314 : : ! 0/ line or column --> line 0
315 : : ! Cubic, Fcc, Bcc, Ortho, Tetra, Rhombo, Hexa
316 : : ! 1/ line only --> line=1 (see rprim using acell in ABINIT) :
317 : : ! Mono, Tri
318 : : ! 2/ column only --> line=2 (see rprim using scalecart in ABINIT) :
319 : : ! Bct, Face-centered-Ortho, Body-centered-Ortho, C-centered-Ortho
320 : : ! 3/ neither line or column --> line=3 (rprim has to be directly dimensioned in ABINIT) :
321 : : ! C-centered-Mono
322 : : ! ---------------------------------------------------------------------------- !
323 : 46 : line=0
324 : : ! For monoclinic: bravais(1)=2
325 [ + + + - ]: 46 : if (bravais(1).eq.2.and.bravais(2).eq.0) then !monoclinic
326 : 1 : line=1
327 : 1 : rprim(1,1)= 1.0d0 ; rprim(1,2)= 0.0d0 ; rprim(1,3)= 0.0d0
328 : 1 : rprim(2,1)= 0.0d0 ; rprim(2,2)= 1.0d0 ; rprim(2,3)= 0.0d0
329 : 1 : rprim(3,1)= dcos(angle_alpha*pi/180.d0) ; rprim(3,2)= 0.0d0 ; rprim(3,3)= dsin(angle_alpha*pi/180.d0)
330 : : ! For orthorhombic: bravais(1)=3
331 [ + + + + ]: 45 : else if (bravais(1).eq.3.and.bravais(2).eq.0) then !orthorhombic
332 : 2 : line=0
333 : 2 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
334 : 2 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
335 : 2 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
336 [ + + - + ]: 43 : else if (bravais(1).eq.3.and.bravais(2).eq.-3) then !face centered orthorhombic
337 : 0 : line=2
338 : 0 : rprim(1,1)=0.0d0 ; rprim(1,2)=0.5d0 ; rprim(1,3)=0.5d0
339 : 0 : rprim(2,1)=0.5d0 ; rprim(2,2)=0.0d0 ; rprim(2,3)=0.5d0
340 : 0 : rprim(3,1)=0.5d0 ; rprim(3,2)=0.5d0 ; rprim(3,3)=0.0d0
341 [ + + - + ]: 43 : else if (bravais(1).eq.3.and.bravais(2).eq.-1) then !body centered orthorhombic
342 : 0 : line=2
343 : 0 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
344 : 0 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
345 : 0 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
346 [ + + + - ]: 43 : else if (bravais(1).eq.3.and.bravais(2).eq.3) then !orthorombic C-face centered
347 : 7 : line=2
348 : 7 : rprim(1,1)= 0.5d0 ; rprim(1,2)=-0.5d0 ; rprim(1,3)=0.0d0
349 : 7 : rprim(2,1)= 0.5d0 ; rprim(2,2)= 0.5d0 ; rprim(2,3)=0.0d0
350 : 7 : rprim(3,1)= 0.0d0 ; rprim(3,2)= 0.0d0 ; rprim(3,3)=1.0d0
351 : : ! For tetragonal: bravais(1)=4
352 [ + + - + ]: 36 : else if (bravais(1).eq.4.and.bravais(2).eq.0) then !tetragonal
353 : 0 : line=2
354 : 0 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
355 : 0 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
356 : 0 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
357 [ + + + - ]: 36 : else if (bravais(1).eq.4.and.bravais(2).eq.-1) then !body centered tetragonal
358 : 2 : line=2
359 : 2 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
360 : 2 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
361 : 2 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
362 : : ! For trigonal: bravais(1)=5
363 [ + + + - ]: 34 : else if (bravais(1).eq.5.and.bravais(2).eq.0) then !rhombo
364 : 3 : line=0
365 : 3 : xi=dsin(angle_alpha*pi/180.d0/2d0)
366 : 3 : hh=dsqrt(1d0-4d0/3d0*xi**2)
367 : 3 : rprim(1,1)=xi ; rprim(1,2)=-xi/dsqrt(3d0) ; rprim(1,3)= hh
368 : 3 : rprim(2,1)= 0.d0 ; rprim(2,2)=2d0*xi/dsqrt(3d0) ; rprim(2,3)= hh
369 : 3 : rprim(3,1)=-xi ; rprim(3,2)=-xi/dsqrt(3d0) ; rprim(3,3)= hh
370 : : ! double precision :: AA,BB,CC,DD
371 : : ! AA=dcos(angle_alpha*pi/180.d0/2d0)
372 : : ! BB=dsin(angle_alpha*pi/180.d0/2d0)
373 : : ! CC=dcos(angle_alpha*pi/180.d0)
374 : : ! DD=dsqrt(1-CC**2/AA**2)
375 : : ! rprim(1,1)= AA ; rprim(1,2)= AA ; rprim(1,3)= CC/AA
376 : : ! rprim(2,1)=-BB ; rprim(2,2)= BB ; rprim(2,3)= 0.d0
377 : : ! rprim(3,1)= 0.d0 ; rprim(3,2)= 0.d0; rprim(3,3)= DD
378 : : ! rprim(1,1)= AA ; rprim(1,2)= -BB ; rprim(1,3)= 0.d0
379 : : ! rprim(2,1)= AA ; rprim(2,2)= BB ; rprim(2,3)= 0.d0
380 : : ! rprim(3,1)= CC/AA ; rprim(3,2)= 0.d0; rprim(3,3)= DD
381 : : ! For hexagonal: bravais(1)=6
382 [ + + + - ]: 31 : else if (bravais(1).eq.6.and.bravais(2).eq.0) then !hexagonal
383 : 5 : line=0
384 : : ! The following definition of the hcp is fixed in m_dynmat (chkrp9)
385 : 5 : rprim(1,1)= 1.0d0 ; rprim(1,2)= 0.0d0 ; rprim(1,3)= 0.0d0
386 : 5 : rprim(2,1)=-0.5d0 ; rprim(2,2)= dsqrt(3.d0)/2.d0 ; rprim(2,3)= 0.0d0
387 : 5 : rprim(3,1)= 0.0d0 ; rprim(3,2)= 0.0d0 ; rprim(3,3)= 1.0d0
388 : : ! For cubic: bravais(1)=7
389 [ + - + + ]: 26 : else if (bravais(1).eq.7.and.bravais(2).eq.0) then !simple cubic
390 : 6 : line=0
391 : 6 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
392 : 6 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
393 : 6 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
394 [ + - + + ]: 20 : else if (bravais(1).eq.7.and.bravais(2).eq.-3) then !face centered cubic
395 : 16 : line=0
396 : 16 : rprim(1,1)=0.0d0 ; rprim(1,2)=0.5d0 ; rprim(1,3)=0.5d0
397 : 16 : rprim(2,1)=0.5d0 ; rprim(2,2)=0.0d0 ; rprim(2,3)=0.5d0
398 : 16 : rprim(3,1)=0.5d0 ; rprim(3,2)=0.5d0 ; rprim(3,3)=0.0d0
399 [ + - + - ]: 4 : else if (bravais(1).eq.7.and.bravais(2).eq.-1) then !body centered cubic
400 : 4 : line=0
401 : 4 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
402 : 4 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
403 : 4 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
404 : : else
405 : 0 : write(msg,'(a,2(1x,i3),a)') 'The Bravais lattice of the unitcell, bravais(1:2)=',bravais(1:2),&
406 : 0 : ', is not available in aTDEP.'
407 : 0 : ABI_ERROR(msg)
408 : : end if
409 : :
410 : : ! ---------------------------------------------------------------------------- !
411 : :
412 : : ! Define inverse of the multiplicity
413 [ + + + + ]: 598 : multiplicity = Invar%multiplicity
414 : 46 : call matr3inv(multiplicity,multiplicitym1)
415 [ + + + + : 1150 : multiplicitym1 = TRANSPOSE(multiplicitym1)
+ + + + ]
416 : :
417 : : ! Compute gprim and (transpose of gprim) gprimt
418 : 46 : call matr3inv(rprim, Lattice%gprimt)
419 [ + + + + : 1150 : Lattice%gprim = TRANSPOSE(Lattice%gprimt)
+ + + + ]
420 : : ! GA: This is a weird convention, and is different from the rest of abinit.
421 : : ! See for example how the DDB has to be initialized with Lattice.
422 : :
423 : : ! Define transpose and inverse of rprim
424 [ + + + + ]: 598 : rprimt = TRANSPOSE(rprim)
425 : 46 : call matr3inv(rprimt, rprimm1)
426 : :
427 : : ! Compute acell_unitcell and a rotation matrix, given that
428 : : ! the supercell rprimd is related to the unitcell rprim according to
429 : : ! rprimd = R * A * multiplicity * rprim
430 : : ! where R is unitary and A is a diagonal matrix containing acell.
431 [ + + + + : 3634 : RAmat = MATMUL(MATMUL(Invar%rprimd_md, rprimm1), multiplicitym1)
+ + + + +
+ + + ]
432 : :
433 : : ! This matrix should be diagonal
434 [ + + + + : 1840 : Amat2 = MATMUL(TRANSPOSE(RAmat), RAmat)
+ + ]
435 [ + + ]: 184 : do ii=1,3
436 : 184 : acell_unitcell(ii) = sqrt(Amat2(ii,ii))
437 : : end do
438 [ + + ]: 184 : do ii=1,3
439 [ + + ]: 598 : do jj=1,3
440 : 552 : Rmat(ii,jj) = RAmat(ii,jj) / acell_unitcell(jj)
441 : : end do
442 : : end do
443 [ + + + + ]: 598 : rotation = TRANSPOSE(Rmat)
444 : :
445 : : ! Compute rotation matrix in cartesian coordinates
446 : 46 : call matr3inv(Invar%rprimd_md, rprimdm1_tmp)
447 [ + + + + : 1150 : rprimdm1_tmp = TRANSPOSE(rprimdm1_tmp)
+ + + + ]
448 [ + + + + : 1840 : rotation_cart = MATMUL(rotation, Invar%rprimd_md)
+ + ]
449 [ + + + + : 1840 : rotation_cart = MATMUL(rprimdm1_tmp, rotation_cart)
+ + ]
450 [ + + + + : 1150 : rotation_cart = TRANSPOSE(rotation_cart)
+ + + + ]
451 : :
452 : : ! Perform some checks
453 [ - + ]: 46 : if ((mat33det(rotation) - 1) .gt. tol8) then
454 [ # # # # : 0 : rprimd_tmp = MATMUL(multiplicitym1, Invar%rprimd_md)
# # ]
455 : : write(msg, '(6a,3(3f16.10,1x,a),2a,3(3f16.10,1x,a))')&
456 : 0 : 'The input primitive vectors cannot be aligned',ch10,&
457 : 0 : 'with the expected primitive vectors through rotation.',ch10,&
458 : 0 : 'Input unitcell primitive vectors:',ch10,&
459 : 0 : (rprimd_tmp(1,jj),jj=1,3),ch10,&
460 : 0 : (rprimd_tmp(2,jj),jj=1,3),ch10,&
461 : 0 : (rprimd_tmp(3,jj),jj=1,3),ch10,&
462 : 0 : 'Expected primitive vectors:',ch10,&
463 : 0 : (rprim(1,jj),jj=1,3),ch10,&
464 : 0 : (rprim(2,jj),jj=1,3),ch10,&
465 : 0 : (rprim(3,jj),jj=1,3),ch10
466 : 0 : ABI_ERROR(msg)
467 : : end if
468 : :
469 : : ! Apply rotation to rprim_md
470 [ + + + + : 2944 : Invar%rprimd_md = MATMUL(rotation, Invar%rprimd_md)
+ + + + +
+ + + +
+ ]
471 : :
472 : : ! ---------------------------------------------------------------------------- !
473 : :
474 : : ! Recompute dimensioned primitive vectors
475 [ + + + + : 1840 : rprimd_md = MATMUL(multiplicity, rprim)
+ + ]
476 [ + + ]: 46 : if (line==0.or.line==1) then
477 [ + + ]: 148 : do ii=1,3
478 [ + + ]: 481 : do jj=1,3
479 : 444 : rprimd_md(ii,jj) = acell_unitcell(ii) * rprimd_md(ii,jj)
480 : : end do
481 : : end do
482 : : else if (line==2) then
483 [ + + ]: 36 : do ii=1,3
484 [ + + ]: 117 : do jj=1,3
485 : 108 : rprimd_md(ii,jj) = acell_unitcell(jj) * rprimd_md(ii,jj)
486 : : end do
487 : : end do
488 : : end if
489 : :
490 : : ! Echo some (re)computed quantities
491 : 46 : write(ab_out,*) ' '
492 : 46 : write(ab_out,*) '#############################################################################'
493 : 46 : write(ab_out,*) '########################## Computed quantities ##############################'
494 : 46 : write(ab_out,*) '#############################################################################'
495 : :
496 : : ! Check the off-diagonal elements
497 [ + + ]: 184 : do ii=1,3
498 : 184 : write(std_out,'(a,1x,3(f16.10,1x))') 'The rprimd_md (computed)=',(rprimd_md(ii,jj),jj=1,3)
499 : : end do
500 : : if((abs(rprimd_md(1,2)-Invar%rprimd_md(1,2)).gt.tol5).or.(abs(rprimd_md(1,3)-Invar%rprimd_md(1,3)).gt.tol5).or.&
501 : : & (abs(rprimd_md(2,1)-Invar%rprimd_md(2,1)).gt.tol5).or.(abs(rprimd_md(3,1)-Invar%rprimd_md(3,1)).gt.tol5).or.&
502 : : & (abs(rprimd_md(2,3)-Invar%rprimd_md(2,3)).gt.tol5).or.(abs(rprimd_md(3,2)-Invar%rprimd_md(3,2)).gt.tol5).or.&
503 [ - + - + : 46 : & (abs(rprimd_md(1,1)-Invar%rprimd_md(1,1)).gt.tol5).or.(abs(rprimd_md(2,2)-Invar%rprimd_md(2,2)).gt.tol5).or.&
- + - + -
+ - + - +
- + - + ]
504 : : & (abs(rprimd_md(3,3)-Invar%rprimd_md(3,3)).gt.tol5)) then
505 [ # # ]: 0 : do ii=1,3
506 : 0 : write(std_out,'(a,1x,3(f16.10,1x))') 'The rprimd (from the input file or NetCDF file) is=',&
507 [ # # # # ]: 0 : & (Invar%rprimd_md(ii,jj),jj=1,3)
508 : : end do
509 [ # # ]: 0 : do ii=1,3
510 : 0 : write(std_out,'(a,1x,3(f16.10,1x))') 'However, using multiplicity (from the input file)=',&
511 : 0 : & (multiplicity(ii,jj),jj=1,3)
512 : : end do
513 [ # # ]: 0 : do ii=1,3
514 : 0 : write(std_out,'(a,1x,3(f16.10,1x))') 'rprim (from the aTDEP code)=',(rprim(ii,jj),jj=1,3)
515 : : end do
516 : 0 : write(std_out,'(a,1x,3(f16.10,1x))') 'and acell (from the calculation)=',(acell_unitcell(ii),ii=1,3)
517 : 0 : ABI_ERROR(' RPRIMD IS NOT RELATED TO RPRIM AND MULTIPLICITY. MODIFY YOUR RPRIMD.')
518 : : end if
519 : :
520 : : ! Check the diagonal elements
521 [ + + + - ]: 46 : if ((bravais(1).eq.2).and.(bravais(2).eq.0)) then !monoclinic
522 [ + - - + ]: 1 : if ((acell_unitcell(1).gt.acell_unitcell(3)).or.&
523 : : & (acell_unitcell(2).gt.acell_unitcell(3))) then
524 : 0 : ABI_ERROR('You must set a,b <= c in the conventional lattice')
525 : : end if
526 [ + + + + ]: 45 : else if ((bravais(1).eq.3).and.(bravais(2).eq.3)) then !C face centered orthorombique
527 [ - + ]: 7 : if (acell_unitcell(1).ge.acell_unitcell(2)) then
528 : 0 : ABI_ERROR('You must set a < b in the conventional lattice')
529 : : end if
530 [ + + ]: 38 : else if (bravais(1).eq.6) then !hexagonal
531 [ - + ]: 5 : if(abs(acell_unitcell(1)-acell_unitcell(2)).gt.tol8) then
532 : 0 : ABI_ERROR(' STOP: THE PRECISION ON THE LATTICE PARAMETERS IS NOT SUFFICIENT')
533 : : end if
534 : 5 : acell_unitcell(2)=acell_unitcell(1)
535 [ + + ]: 33 : else if (bravais(1).eq.7) then !cubic
536 [ + - - + ]: 26 : if((abs(acell_unitcell(1)-acell_unitcell(2)).gt.tol8).or.&
537 : : & (abs(acell_unitcell(2)-acell_unitcell(3)).gt.tol8)) then
538 : 0 : ABI_ERROR('THE PRECISION ON THE LATTICE PARAMETERS IS NOT SUFFICIENT')
539 : : end if
540 : 26 : acell_unitcell(2)=acell_unitcell(1)
541 : 26 : acell_unitcell(3)=acell_unitcell(1)
542 : : end if
543 : 46 : write(ab_out,'(a,1x,3(f16.10,1x))') ' acell_unitcell=',acell_unitcell(:)
544 : :
545 : : ! Recompute rprimd_md with the symmetrized acell,
546 : : ! in order to have a precision higher than 1.d-8.
547 [ + + + + : 1840 : rprimd_md = MATMUL(multiplicity, rprim)
+ + ]
548 [ + + ]: 46 : if (line==0.or.line==1) then
549 [ + + ]: 148 : do ii=1,3
550 [ + + ]: 481 : do jj=1,3
551 : 444 : rprimd_md(ii,jj) = acell_unitcell(ii) * rprimd_md(ii,jj)
552 : : end do
553 : : end do
554 : : else if (line==2) then
555 [ + + ]: 36 : do ii=1,3
556 [ + + ]: 117 : do jj=1,3
557 : 108 : rprimd_md(ii,jj) = acell_unitcell(jj) * rprimd_md(ii,jj)
558 : : end do
559 : : end do
560 : : end if
561 [ + + ]: 184 : do ii=1,3
562 : 184 : write(ab_out,'(a,1x,3(f16.10,1x))') ' rprimd_md=',(rprimd_md(ii,jj),jj=1,3)
563 : : end do
564 : :
565 : : ! Define the rprimd with respect to the acell_unitcell, according to the value of "line"
566 [ + + ]: 46 : if (line==2) then
567 : 9 : rprimd(1,1)=rprim(1,1)*acell_unitcell(1) ; rprimd(1,2)=rprim(1,2)*acell_unitcell(2) ; rprimd(1,3)=rprim(1,3)*acell_unitcell(3)
568 : 9 : rprimd(2,1)=rprim(2,1)*acell_unitcell(1) ; rprimd(2,2)=rprim(2,2)*acell_unitcell(2) ; rprimd(2,3)=rprim(2,3)*acell_unitcell(3)
569 : 9 : rprimd(3,1)=rprim(3,1)*acell_unitcell(1) ; rprimd(3,2)=rprim(3,2)*acell_unitcell(2) ; rprimd(3,3)=rprim(3,3)*acell_unitcell(3)
570 : : else if (line==0.or.line==1) then
571 : 37 : rprimd(1,1)=rprim(1,1)*acell_unitcell(1) ; rprimd(1,2)=rprim(1,2)*acell_unitcell(1) ; rprimd(1,3)=rprim(1,3)*acell_unitcell(1)
572 : 37 : rprimd(2,1)=rprim(2,1)*acell_unitcell(2) ; rprimd(2,2)=rprim(2,2)*acell_unitcell(2) ; rprimd(2,3)=rprim(2,3)*acell_unitcell(2)
573 : 37 : rprimd(3,1)=rprim(3,1)*acell_unitcell(3) ; rprimd(3,2)=rprim(3,2)*acell_unitcell(3) ; rprimd(3,3)=rprim(3,3)*acell_unitcell(3)
574 : : end if
575 : :
576 : : ! Define transpose and inverse of rprimd
577 [ + + + + ]: 598 : rprimdt = TRANSPOSE(rprimd)
578 : :
579 : : ! Compute gmet, rmet, gprimd
580 : : ! ucvol is discarded here: it is recomputed from Lattice%rprimd and stored
581 : : ! in Model%ucvol, once the Model object is available (see tdep_model_init).
582 : 46 : call metric(Lattice%gmet,Lattice%gprimd,std_out,Lattice%rmet,rprimdt,ucvol)
583 : :
584 : 46 : call matr3inv(rprimd, rprimdtm1)
585 : 46 : call matr3inv(rprimdt, rprimdm1)
586 : :
587 : : ! Store all these values in the 'Lattice' datatype
588 [ + + ]: 184 : Lattice%acell_unitcell(:) =acell_unitcell(:)
589 : 46 : Lattice%angle_alpha =angle_alpha
590 [ + + ]: 552 : Lattice%bravais(:) =bravais(:)
591 : 46 : Lattice%line =line
592 [ + + + + ]: 598 : Lattice%multiplicity (:,:)=multiplicity(:,:)
593 [ + + + + ]: 598 : Lattice%multiplicitym1(:,:)=multiplicitym1(:,:)
594 [ + + + + ]: 598 : Lattice%rprim (:,:)=rprim (:,:)
595 [ + + + + ]: 598 : Lattice%rprimt (:,:)=rprimt (:,:)
596 [ + + + + ]: 598 : Lattice%rprimm1 (:,:)=rprimm1 (:,:)
597 [ + + + + ]: 598 : Lattice%rprimd (:,:)=rprimd (:,:)
598 [ + + + + ]: 598 : Lattice%rprimdm1 (:,:)=rprimdm1 (:,:)
599 [ + + + + ]: 598 : Lattice%rprimdt (:,:)=rprimdt (:,:)
600 [ + + + + ]: 598 : Lattice%rprimdtm1 (:,:)=rprimdtm1 (:,:)
601 [ + + + + ]: 598 : Lattice%rprimd_md (:,:)=rprimd_md (:,:)
602 : :
603 : 46 : end subroutine tdep_make_latt
604 : : !!***
605 : :
606 : : !=====================================================================================================
607 : :
608 : 0 : end module m_tdep_latt
609 : : !!***
|
