Line data Source code
1 : !!****m* ABINIT/m_strain
2 : !!
3 : !! NAME
4 : !! m_strain
5 : !!
6 : !! FUNCTION
7 : !! Module for get the strain
8 : !! Container type is defined
9 : !!
10 : !! COPYRIGHT
11 : !! Copyright (C) 2010-2026 ABINIT group (AM)
12 : !! This file is distributed under the terms of the
13 : !! GNU General Public Licence, see ~abinit/COPYING
14 : !! or http://www.gnu.org/copyleft/gpl.txt .
15 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
16 : !!
17 : !! SOURCE
18 :
19 : #if defined HAVE_CONFIG_H
20 : #include "config.h"
21 : #endif
22 :
23 : #include "abi_common.h"
24 :
25 : module m_strain
26 :
27 : use defs_basis
28 : use m_errors
29 : use m_abicore
30 : use m_xmpi
31 :
32 : use m_matrix, only : matr3inv
33 :
34 : implicit none
35 :
36 : private :: strain_def2strain
37 : private :: strain_strain2def
38 : public :: strain_print
39 : public :: strain_get
40 : public :: strain_init
41 : public :: strain_apply
42 : !!***
43 :
44 : !!****t* m_strain/strain_type
45 : !! NAME
46 : !! strain_type
47 : !!
48 : !! FUNCTION
49 : !! structure for a effective potential constructed.
50 : !!
51 : !! SOURCE
52 :
53 : type, public :: strain_type
54 : character(len=fnlen) :: name
55 : ! name of the strain (iso,uniaxial,shear...)
56 :
57 : real(dp) :: delta
58 : ! Value of the strain
59 :
60 : integer :: direction
61 : ! Direction of the strain (-1 if isostatic)
62 :
63 : real(dp) :: strain(3,3)
64 : ! Matrix representing the strain
65 :
66 : end type strain_type
67 : !!***
68 :
69 : CONTAINS !===========================================================================================
70 :
71 : !****f* m_strain/strain_init
72 : !!
73 : !! NAME
74 : !! strain_init
75 : !!
76 : !! FUNCTION
77 : !! routine to initialize strain structure
78 : !!
79 : !!
80 : !! INPUTS
81 : !! name = name of the perturbation
82 : !! direction = direction of the perturbation
83 : !! delta = delta to apply in the strain (in percent)
84 : !!
85 : !! OUTPUT
86 : !! strain = structure with all information of strain
87 : !!
88 : !! SOURCE
89 :
90 0 : subroutine strain_init(strain,delta,direction,name)
91 :
92 : !Arguments ------------------------------------
93 : !scalars
94 : character(len=fnlen),optional,intent(in) :: name
95 : real(dp),optional,intent(in) :: delta
96 : integer,optional,intent(in) :: direction
97 : !array
98 : type(strain_type),intent(out) :: strain
99 : !Local variables-------------------------------
100 : !scalar
101 : !arrays
102 : ! *************************************************************************
103 0 : if (present(name)) then
104 0 : strain%name = name
105 : else
106 0 : strain%name = ''
107 : end if
108 :
109 0 : if (present(delta)) then
110 0 : strain%delta = delta
111 : else
112 0 : strain%delta = zero
113 : end if
114 :
115 0 : if (present(direction)) then
116 0 : strain%direction = direction
117 : else
118 0 : strain%direction = 0
119 : end if
120 :
121 0 : call strain_strain2def(strain%strain,strain)
122 :
123 0 : end subroutine strain_init
124 : !!***
125 :
126 : !****f* m_strain/strain_free
127 : !!
128 : !! NAME
129 : !! strain_free
130 : !!
131 : !! FUNCTION
132 : !! routine to free strain structure
133 : !!
134 : !!
135 : !! INPUTS
136 : !! strain = structure with all information of strain
137 : !!
138 : !! OUTPUT
139 : !!
140 : !! SOURCE
141 :
142 0 : subroutine strain_free(strain)
143 :
144 : !Arguments ------------------------------------
145 : !scalars
146 : !array
147 : type(strain_type),intent(inout) :: strain
148 :
149 : !Local variables-------------------------------
150 : !scalar
151 : !arrays
152 : ! *************************************************************************
153 :
154 0 : strain%name = ''
155 0 : strain%delta = zero
156 0 : strain%direction = 0
157 0 : strain%strain = zero
158 :
159 0 : end subroutine strain_free
160 : !!***
161 :
162 :
163 : !****f* m_strain/strain_get
164 : !!
165 : !! NAME
166 : !! strain_get
167 : !!
168 : !! FUNCTION
169 : !! Get the strain for structure, compare to reference
170 : !! structure and fill strain type
171 : !!
172 : !!
173 : !! INPUTS
174 : !! symmetrized = (optional) symmetrize the output
175 : !!
176 : !! OUTPUT
177 : !! strain = structure with all information of strain
178 : !!
179 : !! SOURCE
180 :
181 27999 : subroutine strain_get(strain,rprim,rprim_def,mat_delta,symmetrized)
182 :
183 : !Arguments ------------------------------------
184 : !scalars
185 : !array
186 : type(strain_type),intent(inout) :: strain
187 : real(dp),optional,intent(in) :: rprim(3,3),rprim_def(3,3), mat_delta(3,3)
188 : logical,optional,intent(in) :: symmetrized
189 : !Local variables-------------------------------
190 : !scalar
191 : integer :: i,j
192 : logical :: symmetrized_in
193 : character(len=500) :: message
194 : !arrays
195 : real(dp) :: mat_delta_tmp(3,3),rprim_inv(3,3)
196 : real(dp) :: identity(3,3)
197 : ! *************************************************************************
198 :
199 : !check inputs
200 27999 : symmetrized_in = .FALSE.
201 27999 : if(present(symmetrized)) then
202 13519 : symmetrized_in = symmetrized
203 : end if
204 :
205 27999 : if((present(rprim_def).and..not.present(rprim)).or.&
206 27999 : & (present(rprim).and..not.present(rprim_def))) then
207 : write(message, '(a)' )&
208 0 : & ' strain_get: should give rprim_def and rprim as input of the routines'
209 0 : ABI_BUG(message)
210 : end if
211 :
212 27999 : if(present(rprim_def).and.present(rprim))then
213 : mat_delta_tmp = zero
214 : ! Fill the identity matrix
215 27038 : identity = zero
216 108152 : forall(i=1:3)identity(i,i)=1
217 :
218 27038 : call matr3inv(rprim,rprim_inv)
219 1433014 : mat_delta_tmp = matmul(rprim_def,transpose(rprim_inv))-identity
220 27038 : identity = zero
221 108152 : do i=1,3
222 351494 : do j=1,3
223 324456 : if (abs(mat_delta_tmp(i,j))>tol10) then
224 76532 : identity(i,j) = mat_delta_tmp(i,j)
225 : end if
226 : end do
227 : end do
228 :
229 27038 : mat_delta_tmp = identity
230 :
231 961 : else if (present(mat_delta)) then
232 961 : mat_delta_tmp = mat_delta
233 :
234 : else
235 : write(message, '(a)' )&
236 0 : & ' strain_get: should give rprim_def or mat_delta as input of the routines'
237 0 : ABI_BUG(message)
238 : end if
239 :
240 27999 : if(symmetrized_in)then
241 0 : mat_delta_tmp(2,3) = (mat_delta_tmp(2,3) + mat_delta_tmp(3,2)) / 2
242 0 : mat_delta_tmp(3,1) = (mat_delta_tmp(3,1) + mat_delta_tmp(1,3)) / 2
243 0 : mat_delta_tmp(2,1) = (mat_delta_tmp(2,1) + mat_delta_tmp(1,2)) / 2
244 :
245 0 : mat_delta_tmp(3,2) = mat_delta_tmp(2,3)
246 0 : mat_delta_tmp(1,3) = mat_delta_tmp(3,1)
247 0 : mat_delta_tmp(1,2) = mat_delta_tmp(2,1)
248 :
249 : end if
250 :
251 27999 : call strain_def2strain(mat_delta_tmp,strain)
252 :
253 27999 : end subroutine strain_get
254 : !!***
255 :
256 : !****f* m_strain/strain_apply
257 : !!
258 : !! NAME
259 : !! strain_get
260 : !!
261 : !! FUNCTION
262 : !! Get the strain for structure, compare to reference
263 : !! structure and fill strain type
264 : !!
265 : !!
266 : !! INPUTS
267 : !!
268 : !! OUTPUT
269 : !! strain = structure with all information of strain
270 : !!
271 : !! SOURCE
272 :
273 0 : subroutine strain_apply(rprim,rprim_def,strain)
274 :
275 : !Arguments ------------------------------------
276 : !scalars
277 : !array
278 : real(dp),intent(in) :: rprim(3,3)
279 : real(dp),intent(out) :: rprim_def(3,3)
280 : type(strain_type),intent(in) :: strain
281 : !Local variables-------------------------------
282 : !scalar
283 : !integer :: i
284 : !arrays
285 : ! *************************************************************************
286 :
287 0 : rprim_def(:,:) = zero
288 : ! Fill the identity matrix
289 0 : rprim_def(:,:) = matmul(strain%strain(:,:),transpose(rprim(:,:)))
290 :
291 0 : end subroutine strain_apply
292 : !!***
293 :
294 : !****f* m_strain/strain_def2strain
295 : !!
296 : !! NAME
297 : !! strain_matdef2strain
298 : !!
299 : !! FUNCTION
300 : !! transfer deformation matrix in structure strain
301 : !!
302 : !! INPUTS
303 : !! rprim = contains
304 : !!
305 : !! OUTPUT
306 : !!
307 : !!
308 : !! SOURCE
309 :
310 27999 : subroutine strain_def2strain(mat_strain,strain)
311 :
312 : !Arguments ------------------------------------
313 : !scalars
314 : !array
315 : real(dp),intent(in) :: mat_strain(3,3)
316 : type(strain_type),intent(inout) :: strain
317 : !Local variables-------------------------------
318 : !scalar
319 : !arrays
320 : ! *************************************************************************
321 27999 : strain%name = ""
322 27999 : strain%delta = zero
323 27999 : strain%direction = 0
324 363987 : strain%strain = mat_strain
325 :
326 241167 : if (all(abs(mat_strain)<tol10)) then
327 17764 : strain%name = "reference"
328 : strain%delta = zero
329 : strain%direction = 0
330 230932 : strain%strain = zero
331 : else
332 : if(abs(mat_strain(1,1))>tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))<tol10.and.&
333 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))>tol10.and.abs(mat_strain(2,3))<tol10.and.&
334 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))>tol10) then
335 1321 : if((mat_strain(1,1)-mat_strain(2,2))< tol10.and.&
336 : & (mat_strain(1,1)-mat_strain(3,3))< tol10) then
337 1321 : strain%name = "isostatic"
338 1321 : strain%delta = mat_strain(1,1)
339 1321 : strain%direction = -1
340 17173 : strain%strain = mat_strain
341 : end if
342 : end if
343 : if(abs(mat_strain(1,1))>tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))<tol10.and.&
344 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))<tol10.and.abs(mat_strain(2,3))<tol10.and.&
345 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))<tol10) then
346 0 : strain%name = "uniaxial"
347 0 : strain%delta = mat_strain(1,1)
348 0 : strain%direction = 1
349 0 : strain%strain = mat_strain
350 : end if
351 : if(abs(mat_strain(1,1))<tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))<tol10.and.&
352 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))>tol10.and.abs(mat_strain(2,3))<tol10.and.&
353 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))<tol10) then
354 0 : strain%name = "uniaxial"
355 0 : strain%delta = mat_strain(2,2)
356 0 : strain%direction = 2
357 0 : strain%strain = mat_strain
358 : end if
359 : if(abs(mat_strain(1,1))<tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))<tol10.and.&
360 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))<tol10.and.abs(mat_strain(2,3))<tol10.and.&
361 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))>tol10) then
362 0 : strain%name = "uniaxial"
363 0 : strain%delta = mat_strain(3,3)
364 0 : strain%direction = 3
365 0 : strain%strain = mat_strain
366 : end if
367 : if(abs(mat_strain(1,1))<tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))<tol10.and.&
368 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))<tol10.and.abs(mat_strain(2,3))>tol10.and.&
369 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))>tol10.and.abs(mat_strain(3,3))<tol10) then
370 0 : if (abs(mat_strain(3,2)-mat_strain(3,2))<tol10) then
371 0 : strain%name = "shear"
372 0 : strain%delta = mat_strain(3,2) * 2
373 0 : strain%direction = 4
374 0 : strain%strain = mat_strain
375 : end if
376 : end if
377 : if(abs(mat_strain(1,1))<tol10.and.abs(mat_strain(1,2))<tol10.and.abs(mat_strain(1,3))>tol10.and.&
378 : & abs(mat_strain(2,1))<tol10.and.abs(mat_strain(2,2))<tol10.and.abs(mat_strain(2,3))<tol10.and.&
379 10235 : & abs(mat_strain(3,1))>tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))<tol10) then
380 0 : if (abs(mat_strain(3,1)-mat_strain(1,3))<tol10) then
381 0 : strain%name = "shear"
382 0 : strain%delta = mat_strain(3,1) * 2
383 0 : strain%direction = 5
384 0 : strain%strain = mat_strain
385 : end if
386 : end if
387 : if(abs(mat_strain(1,1))<tol10.and.abs(mat_strain(1,2))>tol10.and.abs(mat_strain(1,3))<tol10.and.&
388 : & abs(mat_strain(2,1))>tol10.and.abs(mat_strain(2,2))<tol10.and.abs(mat_strain(2,3))<tol10.and.&
389 10235 : & abs(mat_strain(3,1))<tol10.and.abs(mat_strain(3,2))<tol10.and.abs(mat_strain(3,3))<tol10) then
390 0 : if (abs(mat_strain(1,2)-mat_strain(2,1))<tol10) then
391 0 : strain%name = "shear"
392 0 : strain%delta = mat_strain(2,1) * 2
393 0 : strain%direction = 6
394 0 : strain%strain = mat_strain
395 : end if
396 : end if
397 : end if
398 :
399 27999 : end subroutine strain_def2strain
400 : !!***
401 :
402 : !****f* m_strain/strain_strain2def
403 : !!
404 : !! NAME
405 : !! strain_matdef2strain
406 : !!
407 : !! FUNCTION
408 : !! transfer deformation matrix in structure strain
409 : !!
410 : !! INPUTS
411 : !! rprim = contains
412 : !!
413 : !! OUTPUT
414 : !!
415 : !!
416 : !! SOURCE
417 :
418 0 : subroutine strain_strain2def(mat_strain,strain)
419 :
420 : !Arguments ------------------------------------
421 : !scalars
422 : !array
423 : real(dp),intent(out) :: mat_strain(3,3)
424 : type(strain_type),intent(in) :: strain
425 : !Local variables-------------------------------
426 : !scalar
427 : integer :: i
428 : !arrays
429 : ! *************************************************************************
430 :
431 0 : mat_strain(:,:) = zero
432 0 : forall(i=1:3)mat_strain(i,i)=1
433 :
434 0 : select case(strain%direction)
435 : case(1)
436 0 : mat_strain(1,1) = mat_strain(1,1) + strain%delta
437 : case(2)
438 0 : mat_strain(2,2) = mat_strain(2,2) + strain%delta
439 : case(3)
440 0 : mat_strain(3,3) = mat_strain(3,3) + strain%delta
441 : case(4)
442 0 : mat_strain(3,2) = strain%delta / 2
443 0 : mat_strain(2,3) = strain%delta / 2
444 : case(5)
445 0 : mat_strain(3,1) = strain%delta / 2
446 0 : mat_strain(1,3) = strain%delta / 2
447 : case(6)
448 0 : mat_strain(2,1) = strain%delta / 2
449 0 : mat_strain(1,2) = strain%delta / 2
450 : end select
451 :
452 0 : end subroutine strain_strain2def
453 : !!***
454 :
455 : !****f* m_strain/strain_print
456 : !!
457 : !! NAME
458 : !! strain_print
459 : !!
460 : !! FUNCTION
461 : !! print the structure strain
462 : !!
463 : !! INPUTS
464 : !!
465 : !! OUTPUT
466 : !! eff_pot = supercell structure with data to be output
467 : !!
468 : !! SOURCE
469 :
470 2056 : subroutine strain_print(strain)
471 :
472 : !Arguments ------------------------------------
473 : !scalars
474 : !array
475 : type(strain_type),intent(in) :: strain
476 : !Local variables-------------------------------
477 : !scalar
478 : integer :: ii
479 : character(len=500) :: message
480 : !arrays
481 : ! *************************************************************************
482 :
483 2056 : if(strain%name == "reference") then
484 1081 : write(message,'(4a)') ch10,' no strain found:',&
485 2162 : & ' This structure is equivalent to the reference structure',ch10
486 1081 : call wrtout(std_out,message,'COLL')
487 : else
488 975 : if(strain%name /= "") then
489 : write(message,'(3a,I2,a,(ES10.2),a)') &
490 20 : & ' The strain is ',trim(strain%name),' type in the direction ',&
491 40 : & strain%direction,' with delta of ',strain%delta, ':'
492 20 : call wrtout(std_out,message,'COLL')
493 20 : call wrtout(ab_out,message,'COLL')
494 80 : do ii = 1,3
495 60 : write(message,'(3es17.8)') strain%strain(ii,1),strain%strain(ii,2),strain%strain(ii,3)
496 80 : call wrtout(std_out,message,'COLL')
497 : end do
498 : else
499 955 : write(message,'(a)') ' Strain does not correspond to standard strain:'
500 955 : call wrtout(std_out,message,'COLL')
501 3820 : do ii = 1,3
502 2865 : write(message,'(3es17.8)') strain%strain(ii,1),strain%strain(ii,2),strain%strain(ii,3)
503 3820 : call wrtout(std_out,message,'COLL')
504 : end do
505 : end if
506 : end if
507 2056 : end subroutine strain_print
508 : !!***
509 :
510 0 : end module m_strain
511 :
|
