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