Branch data Line data Source code
1 : :
2 : : #if defined HAVE_CONFIG_H
3 : : #include "config.h"
4 : : #endif
5 : :
6 : : #include "abi_common.h"
7 : :
8 : : module m_tdep_latt
9 : :
10 : : use defs_basis
11 : : use m_abicore
12 : : use m_errors
13 : : use m_xmpi
14 : :
15 : : use m_matrix, only : matr3inv, mat33det
16 : : use m_geometry, only : metric
17 : : use m_tdep_dataset, only : atdep_dataset_type, MPI_enreg_type
18 : :
19 : : implicit none
20 : :
21 : : type Lattice_type
22 : :
23 : : double precision :: acell_unitcell(3)
24 : : double precision :: angle_alpha
25 : : integer :: brav
26 : : integer :: bravais(11)
27 : : integer :: line
28 : : double precision :: metmin (3,3)
29 : : double precision :: minim (3,3)
30 : : double precision :: multiplicity (3,3)
31 : : double precision :: multiplicitym1(3,3)
32 : : double precision :: gmet (3,3)
33 : : double precision :: rmet (3,3)
34 : : double precision :: gprimd (3,3)
35 : : double precision :: gprim (3,3)
36 : : double precision :: gprimt (3,3)
37 : : double precision :: rprim (3,3)
38 : : double precision :: rprimt (3,3)
39 : : double precision :: rprimm1 (3,3)
40 : : double precision :: rprimd (3,3)
41 : : double precision :: rprimdm1 (3,3)
42 : : double precision :: rprimdt (3,3)
43 : : double precision :: rprimdtm1 (3,3)
44 : : double precision :: rprimd_md (3,3)
45 : : double precision :: Sij (6,6)
46 : : double precision :: ucvol
47 : : double precision :: BulkModulus_T
48 : : double precision :: BulkModulus_S
49 : : double precision :: HeatCapa_V
50 : : double precision :: HeatCapa_P
51 : : double precision :: Shear
52 : : double precision :: Density
53 : :
54 : : end type Lattice_type
55 : :
56 : : public :: tdep_make_inbox
57 : : public :: tdep_make_latt
58 : :
59 : : contains
60 : :
61 : : !=====================================================================================================
62 : 23774195 : subroutine tdep_make_inbox(tab,natom,tol,&
63 : : & temp) !Optional
64 : :
65 : : integer :: natom,ii,jj,iatom
66 : : double precision :: tol
67 : : double precision :: tab(3,natom) !Input (and ouput if temp is not present)
68 : : double precision,optional :: temp(3,natom) !Input/Output (if present)
69 : :
70 : : !GA: One could make use of function wrap2_pmhalf
71 : :
72 [ + + ]: 47548390 : do iatom=1,natom
73 [ + + ]: 118870975 : do ii=1,3
74 [ + + ]: 71322585 : if (tab(ii,iatom).lt.0.d0) then
75 : 22489363 : jj=int(tab(ii,iatom)-0.5+tol)
76 [ + - ]: 48833222 : else if (tab(ii,iatom).ge.0.d0) then
77 : 48833222 : jj=int(tab(ii,iatom)+0.5+tol)
78 : : end if
79 [ + + ]: 95096780 : if (present(temp)) then
80 : 66880881 : temp(ii,iatom)=temp(ii,iatom)-real(jj)
81 : : else
82 : 4441704 : tab(ii,iatom)=tab(ii,iatom)-real(jj)
83 : : end if
84 : : end do
85 : : end do
86 : :
87 : 23774195 : end subroutine tdep_make_inbox
88 : :
89 : : !=====================================================================================================
90 : 44 : subroutine tdep_make_latt(Invar,Lattice,rotation_cart)
91 : :
92 : : type(Lattice_type),intent(out) :: Lattice
93 : : type(atdep_dataset_type),intent(inout) :: Invar
94 : : double precision, intent(out) :: rotation_cart(3,3)
95 : :
96 : : integer :: brav,ii,jj,line
97 : : double precision :: acell_unitcell(3),multiplicity(3,3),multiplicitym1(3,3)
98 : : double precision :: rprimd(3,3),rprimdt(3,3),rprimd_md(3,3),rprimdm1(3,3)
99 : : double precision :: rprim(3,3),rprimm1(3,3),rprimt(3,3),rprimdtm1(3,3)
100 : : double precision :: rprimd_unitcell(3,3)
101 : : double precision :: rprimd_tmp(3,3), rprimdm1_tmp(3,3)
102 : : double precision :: RAmat(3,3),Amat2(3,3),Rmat(3,3)
103 : : double precision :: rotation(3,3)
104 : : double precision :: xi,hh
105 : : character(len=500) :: msg
106 : :
107 : : ! For bravais(1):
108 : : ! The holohedral groups are numbered as follows
109 : : ! (see international tables for crystallography (1983), p. 13)
110 : : ! iholohedry=1 triclinic 1bar
111 : : ! iholohedry=2 monoclinic 2/m
112 : : ! iholohedry=3 orthorhombic mmm
113 : : ! iholohedry=4 tetragonal 4/mmm
114 : : ! iholohedry=5 trigonal 3bar m
115 : : ! iholohedry=6 hexagonal 6/mmm
116 : : ! iholohedry=7 cubic m3bar m
117 : :
118 : : ! For bravais(2):
119 : : ! Centering
120 : : ! center=0 no centering
121 : : ! center=-1 body-centered
122 : : ! center=-3 face-centered
123 : : ! center=1 A-face centered
124 : : ! center=2 B-face centered
125 : : ! center=3 C-face centered
126 : :
127 : : ! Correspondency between brav and bravais: brav=1-S.C., 2-F.C., 3-B.C., 4-Hex.)
128 : :
129 : : ! Initialize some local variables
130 : 44 : acell_unitcell(:)=zero; multiplicity(:,:)=zero; multiplicitym1(:,:)=zero
131 : 44 : rprimd(:,:)=zero; rprimdm1(:,:)=zero; rprimdt(:,:)=zero; rprimd_md(:,:)=zero
132 : 44 : rprim(:,:)=zero; rprimm1(:,:)=zero; rprimt(:,:)=zero; rprimdtm1(:,:)=zero
133 : 44 : rprimd_unitcell(:,:)=zero; RAmat(:,:)=zero; Amat2(:,:)=zero; Rmat(:,:)=zero
134 : :
135 : :
136 : : ! ---------------------------------------------------------------------------- !
137 : : ! Here, rprim defines a Primitive Lattice and NOT a Conventional lattice
138 : : ! The lattice parameters have to multiply rprim on:
139 : : ! 0/ line or column --> line 0
140 : : ! Cubic, Fcc, Bcc, Ortho, Tetra, Rhombo, Hexa
141 : : ! 1/ line only --> line=1 (see rprim using acell in ABINIT) :
142 : : ! Mono, Tri
143 : : ! 2/ column only --> line=2 (see rprim using scalecart in ABINIT) :
144 : : ! Bct, Face-centered-Ortho, Body-centered-Ortho, C-centered-Ortho
145 : : ! 3/ neither line or column --> line=3 (rprim has to be directly dimensioned in ABINIT) :
146 : : ! C-centered-Mono
147 : : ! ---------------------------------------------------------------------------- !
148 : 44 : line=0
149 : : ! For monoclinic: bravais(1)=2
150 [ + + + - ]: 44 : if (Invar%bravais(1).eq.2.and.Invar%bravais(2).eq.0) then !monoclinic
151 : 1 : brav=1
152 : 1 : line=1
153 : 1 : rprim(1,1)= 1.0d0 ; rprim(1,2)= 0.0d0 ; rprim(1,3)= 0.0d0
154 : 1 : rprim(2,1)= 0.0d0 ; rprim(2,2)= 1.0d0 ; rprim(2,3)= 0.0d0
155 : 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)
156 : : ! For orthorhombic: bravais(1)=3
157 [ + + + + ]: 43 : else if (Invar%bravais(1).eq.3.and.Invar%bravais(2).eq.0) then !orthorhombic
158 : 2 : brav=1
159 : 2 : line=0
160 : 2 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
161 : 2 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
162 : 2 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
163 [ + + - + ]: 41 : else if (Invar%bravais(1).eq.3.and.Invar%bravais(2).eq.-3) then !face centered orthorhombic
164 : 0 : brav=1
165 : 0 : line=2
166 : 0 : rprim(1,1)=0.0d0 ; rprim(1,2)=0.5d0 ; rprim(1,3)=0.5d0
167 : 0 : rprim(2,1)=0.5d0 ; rprim(2,2)=0.0d0 ; rprim(2,3)=0.5d0
168 : 0 : rprim(3,1)=0.5d0 ; rprim(3,2)=0.5d0 ; rprim(3,3)=0.0d0
169 [ + + - + ]: 41 : else if (Invar%bravais(1).eq.3.and.Invar%bravais(2).eq.-1) then !body centered orthorhombic
170 : 0 : brav=1
171 : 0 : line=2
172 : 0 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
173 : 0 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
174 : 0 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
175 [ + + + - ]: 41 : else if (Invar%bravais(1).eq.3.and.Invar%bravais(2).eq.3) then !orthorombic C-face centered
176 : 7 : brav=1
177 : 7 : line=2
178 : 7 : rprim(1,1)= 0.5d0 ; rprim(1,2)=-0.5d0 ; rprim(1,3)=0.0d0
179 : 7 : rprim(2,1)= 0.5d0 ; rprim(2,2)= 0.5d0 ; rprim(2,3)=0.0d0
180 : 7 : rprim(3,1)= 0.0d0 ; rprim(3,2)= 0.0d0 ; rprim(3,3)=1.0d0
181 : : ! For tetragonal: bravais(1)=4
182 [ + + - + ]: 34 : else if (Invar%bravais(1).eq.4.and.Invar%bravais(2).eq.0) then !tetragonal
183 : 0 : brav=1
184 : 0 : line=2
185 : 0 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
186 : 0 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
187 : 0 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
188 [ + + + - ]: 34 : else if (Invar%bravais(1).eq.4.and.Invar%bravais(2).eq.-1) then !body centered tetragonal
189 : 2 : brav=3
190 : 2 : line=2
191 : 2 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
192 : 2 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
193 : 2 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
194 : : ! For trigonal: bravais(1)=5
195 [ + + + - ]: 32 : else if (Invar%bravais(1).eq.5.and.Invar%bravais(2).eq.0) then !rhombo
196 : 3 : brav=1
197 : 3 : line=0
198 : 3 : xi=dsin(Invar%angle_alpha*pi/180.d0/2d0)
199 : 3 : hh=dsqrt(1d0-4d0/3d0*xi**2)
200 : 3 : rprim(1,1)=xi ; rprim(1,2)=-xi/dsqrt(3d0) ; rprim(1,3)= hh
201 : 3 : rprim(2,1)= 0.d0 ; rprim(2,2)=2d0*xi/dsqrt(3d0) ; rprim(2,3)= hh
202 : 3 : rprim(3,1)=-xi ; rprim(3,2)=-xi/dsqrt(3d0) ; rprim(3,3)= hh
203 : : ! double precision :: AA,BB,CC,DD
204 : : ! AA=dcos(Invar%angle_alpha*pi/180.d0/2d0)
205 : : ! BB=dsin(Invar%angle_alpha*pi/180.d0/2d0)
206 : : ! CC=dcos(Invar%angle_alpha*pi/180.d0)
207 : : ! DD=dsqrt(1-CC**2/AA**2)
208 : : ! rprim(1,1)= AA ; rprim(1,2)= AA ; rprim(1,3)= CC/AA
209 : : ! rprim(2,1)=-BB ; rprim(2,2)= BB ; rprim(2,3)= 0.d0
210 : : ! rprim(3,1)= 0.d0 ; rprim(3,2)= 0.d0; rprim(3,3)= DD
211 : : ! rprim(1,1)= AA ; rprim(1,2)= -BB ; rprim(1,3)= 0.d0
212 : : ! rprim(2,1)= AA ; rprim(2,2)= BB ; rprim(2,3)= 0.d0
213 : : ! rprim(3,1)= CC/AA ; rprim(3,2)= 0.d0; rprim(3,3)= DD
214 : : ! For hexagonal: bravais(1)=6
215 [ + + + - ]: 29 : else if (Invar%bravais(1).eq.6.and.Invar%bravais(2).eq.0) then !hexagonal
216 : 4 : brav=4
217 : 4 : line=0
218 : : ! The following definition of the hcp is fixed in m_dynmat (chkrp9)
219 : 4 : rprim(1,1)= 1.0d0 ; rprim(1,2)= 0.0d0 ; rprim(1,3)= 0.0d0
220 : 4 : rprim(2,1)=-0.5d0 ; rprim(2,2)= dsqrt(3.d0)/2.d0 ; rprim(2,3)= 0.0d0
221 : 4 : rprim(3,1)= 0.0d0 ; rprim(3,2)= 0.0d0 ; rprim(3,3)= 1.0d0
222 : : ! For cubic: bravais(1)=7
223 [ + - + + ]: 25 : else if (Invar%bravais(1).eq.7.and.Invar%bravais(2).eq.0) then !simple cubic
224 : 6 : brav=1
225 : 6 : line=0
226 : 6 : rprim(1,1)=1.0d0 ; rprim(1,2)=0.0d0 ; rprim(1,3)=0.0d0
227 : 6 : rprim(2,1)=0.0d0 ; rprim(2,2)=1.0d0 ; rprim(2,3)=0.0d0
228 : 6 : rprim(3,1)=0.0d0 ; rprim(3,2)=0.0d0 ; rprim(3,3)=1.0d0
229 [ + - + + ]: 19 : else if (Invar%bravais(1).eq.7.and.Invar%bravais(2).eq.-3) then !face centered cubic
230 : 15 : brav=2
231 : 15 : line=0
232 : 15 : rprim(1,1)=0.0d0 ; rprim(1,2)=0.5d0 ; rprim(1,3)=0.5d0
233 : 15 : rprim(2,1)=0.5d0 ; rprim(2,2)=0.0d0 ; rprim(2,3)=0.5d0
234 : 15 : rprim(3,1)=0.5d0 ; rprim(3,2)=0.5d0 ; rprim(3,3)=0.0d0
235 [ + - + - ]: 4 : else if (Invar%bravais(1).eq.7.and.Invar%bravais(2).eq.-1) then !body centered cubic
236 : 4 : brav=3
237 : 4 : line=0
238 : 4 : rprim(1,1)=-0.5d0 ; rprim(1,2)= 0.5d0 ; rprim(1,3)= 0.5d0
239 : 4 : rprim(2,1)= 0.5d0 ; rprim(2,2)=-0.5d0 ; rprim(2,3)= 0.5d0
240 : 4 : rprim(3,1)= 0.5d0 ; rprim(3,2)= 0.5d0 ; rprim(3,3)=-0.5d0
241 : : else
242 : 0 : ABI_ERROR('THIS BRAVAIS IS NOT DEFINED')
243 : : end if
244 : :
245 : : ! ---------------------------------------------------------------------------- !
246 : :
247 : : ! Define inverse of the multiplicity
248 [ + + + + ]: 572 : multiplicity = Invar%multiplicity
249 : 44 : call matr3inv(multiplicity,multiplicitym1)
250 [ + + + + : 1100 : multiplicitym1 = TRANSPOSE(multiplicitym1)
+ + + + ]
251 : :
252 : : ! Compute gprim and (transpose of gprim) gprimt
253 : 44 : call matr3inv(rprim, Lattice%gprimt)
254 [ + + + + : 1100 : Lattice%gprim = TRANSPOSE(Lattice%gprimt)
+ + + + ]
255 : : ! GA: This is a weird convention, and is different from the rest of abinit.
256 : : ! See for example how the DDB has to be initialized with Lattice.
257 : :
258 : : ! Define transpose and inverse of rprim
259 [ + + + + ]: 572 : rprimt = TRANSPOSE(rprim)
260 : 44 : call matr3inv(rprimt, rprimm1)
261 : :
262 : : ! Compute acell_unitcell and a rotation matrix, given that
263 : : ! the supercell rprimd is related to the unitcell rprim according to
264 : : ! rprimd = R * A * multiplicity * rprim
265 : : ! where R is unitary and A is a diagonal matrix containing acell.
266 [ + + + + : 3476 : RAmat = MATMUL(MATMUL(Invar%rprimd_md, rprimm1), multiplicitym1)
+ + + + +
+ + + ]
267 : :
268 : : ! This matrix should be diagonal
269 [ + + + + : 1760 : Amat2 = MATMUL(TRANSPOSE(RAmat), RAmat)
+ + ]
270 [ + + ]: 176 : do ii=1,3
271 : 176 : acell_unitcell(ii) = sqrt(Amat2(ii,ii))
272 : : end do
273 [ + + ]: 176 : do ii=1,3
274 [ + + ]: 572 : do jj=1,3
275 : 528 : Rmat(ii,jj) = RAmat(ii,jj) / acell_unitcell(jj)
276 : : end do
277 : : end do
278 [ + + + + ]: 572 : rotation = TRANSPOSE(Rmat)
279 : :
280 : : ! Compute rotation matrix in cartesian coordinates
281 : 44 : call matr3inv(Invar%rprimd_md, rprimdm1_tmp)
282 [ + + + + : 1100 : rprimdm1_tmp = TRANSPOSE(rprimdm1_tmp)
+ + + + ]
283 [ + + + + : 1760 : rotation_cart = MATMUL(rotation, Invar%rprimd_md)
+ + ]
284 [ + + + + : 1760 : rotation_cart = MATMUL(rprimdm1_tmp, rotation_cart)
+ + ]
285 [ + + + + : 1100 : rotation_cart = TRANSPOSE(rotation_cart)
+ + + + ]
286 : :
287 : : ! Perform some checks
288 [ - + ]: 44 : if ((mat33det(rotation) - 1) .gt. tol8) then
289 [ # # # # : 0 : rprimd_tmp = MATMUL(multiplicitym1, Invar%rprimd_md)
# # ]
290 : : write(msg, '(6a,3(3f16.10,1x,a),2a,3(3f16.10,1x,a))')&
291 : 0 : 'The input primitive vectors cannot be aligned',ch10,&
292 : 0 : 'with the expected primitive vectors through rotation.',ch10,&
293 : 0 : 'Input unitcell primitive vectors:',ch10,&
294 : 0 : (rprimd_tmp(1,jj),jj=1,3),ch10,&
295 : 0 : (rprimd_tmp(2,jj),jj=1,3),ch10,&
296 : 0 : (rprimd_tmp(3,jj),jj=1,3),ch10,&
297 : 0 : 'Expected primitive vectors:',ch10,&
298 : 0 : (rprim(1,jj),jj=1,3),ch10,&
299 : 0 : (rprim(2,jj),jj=1,3),ch10,&
300 : 0 : (rprim(3,jj),jj=1,3),ch10
301 : 0 : ABI_ERROR(msg)
302 : : end if
303 : :
304 : : ! Apply rotation to rprim_md
305 [ + + + + : 2816 : Invar%rprimd_md = MATMUL(rotation, Invar%rprimd_md)
+ + + + +
+ + + +
+ ]
306 : :
307 : : ! ---------------------------------------------------------------------------- !
308 : :
309 : : ! Recompute dimensioned primitive vectors
310 [ + + + + : 1760 : rprimd_md = MATMUL(multiplicity, rprim)
+ + ]
311 [ + + ]: 44 : if (line==0.or.line==1) then
312 [ + + ]: 140 : do ii=1,3
313 [ + + ]: 455 : do jj=1,3
314 : 420 : rprimd_md(ii,jj) = acell_unitcell(ii) * rprimd_md(ii,jj)
315 : : end do
316 : : end do
317 : : else if (line==2) then
318 [ + + ]: 36 : do ii=1,3
319 [ + + ]: 117 : do jj=1,3
320 : 108 : rprimd_md(ii,jj) = acell_unitcell(jj) * rprimd_md(ii,jj)
321 : : end do
322 : : end do
323 : : end if
324 : :
325 : : ! Echo some (re)computed quantities
326 : 44 : write(Invar%stdout,*) ' '
327 : 44 : write(Invar%stdout,*) '#############################################################################'
328 : 44 : write(Invar%stdout,*) '########################## Computed quantities ##############################'
329 : 44 : write(Invar%stdout,*) '#############################################################################'
330 : :
331 : : ! Check the off-diagonal elements
332 [ + + ]: 176 : do ii=1,3
333 : 176 : write(Invar%stdlog,'(a,1x,3(f16.10,1x))') 'The rprimd_md (computed)=',(rprimd_md(ii,jj),jj=1,3)
334 : : end do
335 : : 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.&
336 : : & (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.&
337 : : & (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.&
338 [ - + - + : 44 : & (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.&
- + - + -
+ - + - +
- + - + ]
339 : : & (abs(rprimd_md(3,3)-Invar%rprimd_md(3,3)).gt.tol5)) then
340 [ # # ]: 0 : do ii=1,3
341 : 0 : write(Invar%stdlog,'(a,1x,3(f16.10,1x))') 'The rprimd (from the input file or NetCDF file) is=',&
342 [ # # # # ]: 0 : & (Invar%rprimd_md(ii,jj),jj=1,3)
343 : : end do
344 [ # # ]: 0 : do ii=1,3
345 : 0 : write(Invar%stdlog,'(a,1x,3(f16.10,1x))') 'However, using multiplicity (from the input file)=',&
346 : 0 : & (multiplicity(ii,jj),jj=1,3)
347 : : end do
348 [ # # ]: 0 : do ii=1,3
349 : 0 : write(Invar%stdlog,'(a,1x,3(f16.10,1x))') 'rprim (from the aTDEP code)=',(rprim(ii,jj),jj=1,3)
350 : : end do
351 : 0 : write(Invar%stdlog,'(a,1x,3(f16.10,1x))') 'and acell (from the calculation)=',(acell_unitcell(ii),ii=1,3)
352 : 0 : ABI_ERROR(' RPRIMD IS NOT RELATED TO RPRIM AND MULTIPLICITY. MODIFY YOUR RPRIMD.')
353 : : end if
354 : :
355 : : ! Check the diagonal elements
356 [ + + + - ]: 44 : if ((Invar%bravais(1).eq.2).and.(Invar%bravais(2).eq.0)) then !monoclinic
357 [ + - - + ]: 1 : if ((acell_unitcell(1).gt.acell_unitcell(3)).or.&
358 : : & (acell_unitcell(2).gt.acell_unitcell(3))) then
359 : 0 : ABI_ERROR('You must set a,b <= c in the conventional lattice')
360 : : end if
361 [ + + + + ]: 43 : else if ((Invar%bravais(1).eq.3).and.(Invar%bravais(2).eq.3)) then !C face centered orthorombique
362 [ - + ]: 7 : if (acell_unitcell(1).ge.acell_unitcell(2)) then
363 : 0 : ABI_ERROR('You must set a < b in the conventional lattice')
364 : : end if
365 [ + + ]: 36 : else if (Invar%bravais(1).eq.6) then !hexagonal
366 [ - + ]: 4 : if(abs(acell_unitcell(1)-acell_unitcell(2)).gt.tol8) then
367 : 0 : ABI_ERROR(' STOP: THE PRECISION ON THE LATTICE PARAMETERS IS NOT SUFFICIENT')
368 : : end if
369 : 4 : acell_unitcell(2)=acell_unitcell(1)
370 [ + + ]: 32 : else if (Invar%bravais(1).eq.7) then !cubic
371 [ + - - + ]: 25 : if((abs(acell_unitcell(1)-acell_unitcell(2)).gt.tol8).or.&
372 : : & (abs(acell_unitcell(2)-acell_unitcell(3)).gt.tol8)) then
373 : 0 : ABI_ERROR('THE PRECISION ON THE LATTICE PARAMETERS IS NOT SUFFICIENT')
374 : : end if
375 : 25 : acell_unitcell(2)=acell_unitcell(1)
376 : 25 : acell_unitcell(3)=acell_unitcell(1)
377 : : end if
378 : 44 : write(Invar%stdout,'(a,1x,3(f16.10,1x))') ' acell_unitcell=',acell_unitcell(:)
379 : :
380 : : ! Recompute rprimd_md with the symmetrized acell,
381 : : ! in order to have a precision higher than 1.d-8.
382 [ + + + + : 1760 : rprimd_md = MATMUL(multiplicity, rprim)
+ + ]
383 [ + + ]: 44 : if (line==0.or.line==1) then
384 [ + + ]: 140 : do ii=1,3
385 [ + + ]: 455 : do jj=1,3
386 : 420 : rprimd_md(ii,jj) = acell_unitcell(ii) * rprimd_md(ii,jj)
387 : : end do
388 : : end do
389 : : else if (line==2) then
390 [ + + ]: 36 : do ii=1,3
391 [ + + ]: 117 : do jj=1,3
392 : 108 : rprimd_md(ii,jj) = acell_unitcell(jj) * rprimd_md(ii,jj)
393 : : end do
394 : : end do
395 : : end if
396 [ + + ]: 176 : do ii=1,3
397 : 176 : write(Invar%stdout,'(a,1x,3(f16.10,1x))') ' rprimd_md=',(rprimd_md(ii,jj),jj=1,3)
398 : : end do
399 : :
400 : : ! Define the rprimd with respect to the acell_unitcell, according to the value of "line"
401 [ + + ]: 44 : if (line==2) then
402 : 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)
403 : 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)
404 : 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)
405 : : else if (line==0.or.line==1) then
406 : 35 : 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)
407 : 35 : 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)
408 : 35 : 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)
409 : : end if
410 : :
411 : : ! Define transpose and inverse of rprimd
412 [ + + + + ]: 572 : rprimdt = TRANSPOSE(rprimd)
413 : :
414 : : ! Compute gmet, rmet, gprimd
415 : 44 : call metric(Lattice%gmet,Lattice%gprimd,Invar%stdlog,Lattice%rmet,rprimdt,Lattice%ucvol)
416 : :
417 : 44 : call matr3inv(rprimd, rprimdtm1)
418 : 44 : call matr3inv(rprimdt, rprimdm1)
419 : :
420 : : ! Store all these values in the 'Lattice' datatype
421 [ + + ]: 176 : Lattice%acell_unitcell(:) =acell_unitcell(:)
422 : 44 : Lattice%angle_alpha =Invar%angle_alpha
423 : 44 : Lattice%brav =brav
424 [ + + ]: 528 : Lattice%bravais(:) =Invar%bravais(:)
425 : 44 : Lattice%line =line
426 [ + + + + ]: 572 : Lattice%multiplicity (:,:)=multiplicity(:,:)
427 [ + + + + ]: 572 : Lattice%multiplicitym1(:,:)=multiplicitym1(:,:)
428 [ + + + + ]: 572 : Lattice%rprim (:,:)=rprim (:,:)
429 [ + + + + ]: 572 : Lattice%rprimt (:,:)=rprimt (:,:)
430 [ + + + + ]: 572 : Lattice%rprimm1 (:,:)=rprimm1 (:,:)
431 [ + + + + ]: 572 : Lattice%rprimd (:,:)=rprimd (:,:)
432 [ + + + + ]: 572 : Lattice%rprimdm1 (:,:)=rprimdm1 (:,:)
433 [ + + + + ]: 572 : Lattice%rprimdt (:,:)=rprimdt (:,:)
434 [ + + + + ]: 572 : Lattice%rprimdtm1 (:,:)=rprimdtm1 (:,:)
435 [ + + + + ]: 572 : Lattice%rprimd_md (:,:)=rprimd_md (:,:)
436 : :
437 : 44 : end subroutine tdep_make_latt
438 : :
439 : : !=====================================================================================================
440 : :
441 : 0 : end module m_tdep_latt
|
