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
|
