Line data Source code
1 : !!****m* ABINIT/m_ddb_piezo
2 : !! NAME
3 : !! m_ddb_piezo
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1999-2026 ABINIT group (XW)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_ddb_piezo
22 :
23 : use defs_basis
24 : use m_abicore
25 : use m_errors
26 : use m_nctk
27 : use netcdf
28 :
29 : use m_fstrings, only : sjoin, itoa
30 : use m_hide_lapack, only : matrginv
31 : use m_anaddb_dataset, only : anaddb_dataset_type
32 :
33 : implicit none
34 :
35 : private
36 : !!***
37 :
38 : public :: ddb_piezo
39 : !!***
40 :
41 : contains
42 : !!***
43 :
44 : !!****f* ABINIT/ddb_piezo
45 : !!
46 : !! NAME
47 : !! ddb_piezo
48 : !!
49 : !! FUNCTION
50 : !! Get the piezoelectric tensor (e-tensor), both clamped ion and relaxed ion;
51 : !! Compute physical(relaxed ion) piezoeletric (d, g, h) tensors;
52 : !! Compute relaxed ion and free stress dielectric tensor;
53 : !! Compute relaxed ion elastic and compliance tensors under fixed
54 : !! displacement field boundary conditions.
55 : !!
56 : !! INPUTS
57 : !! inp= (derived datatype) contains all the input variables
58 : !! blkval(2,3,mpert,3,mpert,nblok)=
59 : !! second derivatives of total energy with respect to electric fields
60 : !! atom displacements,strain,...... all in cartesian coordinates
61 : !! dielt_rlx=relaxed ion dielectric tensor
62 : !! iblok= bolk number in DDB file contains 2 derivative of energy
63 : !! instrain=force response internal strain tensor
64 : !! iout=out file number
65 : !! mpert=maximum number of ipert
66 : !! natom=number of atoms in unit cell
67 : !! nblok=number of total bloks in DDB file
68 : !! ucvol=unit cell volume
69 : !! ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.
70 : !!
71 : !! OUTPUT
72 : !! piezo = piezoelectric tensor
73 : !!
74 : !! NOTES
75 : !! The elastic (compliance) tensors calculated here are under fixed D-field boundary
76 : !! condition, which include piezoelectric corrections to the elastic (compliance)
77 : !! tensors calculated in ddb_elast.F90 whose boundary condition is fixed E-field.
78 : !!
79 : !! SOURCE
80 :
81 5 : subroutine ddb_piezo(inp,blkval,dielt_rlx,elast,iblok,instrain,iout,mpert,natom,nblok,piezo,ucvol,ncid)
82 :
83 : !Arguments-------------------------------------------
84 : !scalars
85 : integer,intent(in) :: iblok,iout,mpert,natom,nblok,ncid
86 : real(dp),intent(in) :: ucvol
87 : type(anaddb_dataset_type),intent(in) :: inp
88 : !arrays
89 : real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),dielt_rlx(3,3)
90 : real(dp),intent(in) :: elast(6,6),instrain(3*natom,6)
91 : real(dp),intent(out) :: piezo(6,3)
92 :
93 : !Local variables---------------------------------------
94 : !scalars
95 : integer :: idir1,idir2,ier,ii1,ii2,ipert1,ipert2,ivarA,ivarB,ncerr
96 : character(len=500) :: message
97 : logical :: iwrite
98 : !arrays
99 10 : real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom)
100 10 : real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2)
101 10 : real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3)
102 10 : real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom),beta_tensor(3,3)
103 : real(dp) :: compliance(6,6),compliance_dis(6,6)
104 10 : real(dp) :: d2cart_relaxed(2,3,mpert,3,mpert,nblok),d_tensor(6,3)
105 10 : real(dp) :: dielt_stress(3,3),eigval(3*natom-3),eigvalp(3*natom)
106 10 : real(dp) :: eigvec(2,3*natom-3,3*natom-3),eigvecp(2,3*natom,3*natom)
107 : real(dp) :: elast_dis(6,6),g_tensor(3,6),h_tensor(3,6)
108 10 : real(dp) :: kmatrix(3*natom,3*natom),new1(6,3*natom),piezo_clamped(6,3)
109 10 : real(dp) :: piezo_correction(6,3),piezo_relaxed(6,3),zhpev1(2,2*3*natom-4)
110 10 : real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2)
111 10 : real(dp) :: zstar1(3,3*natom),zstar2(3*natom,3)
112 :
113 : !****************************************************************
114 :
115 : !extraction of the clamped ion piezoelectric constants from blkvals
116 5 : iwrite = iout > 0
117 :
118 : !the six strain perturbations
119 35 : do ivarA=1,6
120 : ! the three E-field perturbations
121 125 : do ivarB=1,3
122 : ! judge if the ivarA>3 or not
123 90 : if(ivarA>3) then
124 45 : idir1=ivarA-3
125 45 : ipert1=natom+4
126 : ! for the shear part of the strain
127 : else if(ivarA<=3) then
128 45 : idir1=ivarA
129 45 : ipert1=natom+3
130 : ! for the diagonal part of strain
131 : end if
132 90 : idir2=ivarB
133 90 : ipert2=natom+2 !for the E-field perturbation only
134 120 : piezo(ivarA,ivarB)=blkval(1,idir2,ipert2,idir1,ipert1,iblok)
135 : end do
136 : end do
137 :
138 : !consider the volume and the -Qe before the piezo
139 : !according to the (30) in notes, the units are tranformed from atomic units to the SI units
140 35 : do ivarA=1,6
141 125 : do ivarB=1,3
142 120 : piezo(ivarA,ivarB)=piezo(ivarA,ivarB)*AmuBohr2_Cm2
143 : ! now it is in the SI unit
144 : end do
145 : end do
146 :
147 : !give the values of d2cart_relaxed as the same as blkval
148 : !and also give the initial values of piezo_clamped and piezo_relaxed
149 25911 : d2cart_relaxed(:,:,:,:,:,:)=blkval(:,:,:,:,:,:)
150 5 : piezo_clamped(:,:)=piezo(:,:)
151 :
152 : !********************************************************************
153 : !print the main results of the piezoelectric constants
154 5 : if(inp%piezoflag==1.or.inp%piezoflag==3 .or. inp%piezoflag==7)then
155 5 : write(message,'(3a)')ch10,' Proper piezoelectric constants (clamped ion) (unit:c/m^2)',ch10
156 5 : call wrtout(std_out,message,'COLL')
157 :
158 35 : do ivarA=1,6
159 35 : write(std_out,'(3f16.8)')piezo_clamped(ivarA,1),piezo_clamped(ivarA,2),piezo_clamped(ivarA,3)
160 : end do
161 :
162 5 : call wrtout(iout,message,'COLL')
163 5 : if (iwrite) then
164 35 : do ivarA=1,6
165 35 : write(iout,'(3f16.8)')piezo_clamped(ivarA,1),piezo_clamped(ivarA,2),piezo_clamped(ivarA,3)
166 : end do
167 : end if
168 : end if
169 :
170 : !the next is the calculation of the relaxed ion piezoelectric constants
171 : !first extract the K(force constant) matrix
172 :
173 : !if (piezoflag==2 .or. inp%piezoflag==3)then
174 : !extracting force matrix at gamma
175 19 : do ipert1=1,natom
176 61 : do ii1=1,3
177 42 : ivarA=ii1+3*(ipert1-1)
178 188 : do ipert2=1,natom
179 570 : do ii2=1,3
180 396 : ivarB=ii2+3*(ipert2-1)
181 528 : kmatrix(ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok)
182 : end do
183 : end do
184 : end do
185 : end do
186 :
187 443 : Apmatr(:,:)=kmatrix(:,:)
188 :
189 : !DEBUG
190 : !kmatrix values
191 : !write(std_out,'(/,a,/)')'the force constant matrix'
192 : !do ivarA=1,3*natom
193 : !write(std_out,'(/)')
194 : !do ivarB=1,3*natom
195 : !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA)
196 : !end do
197 : !end do
198 : !ENDDEBUG
199 :
200 443 : Nmatr(:,:)=zero
201 :
202 47 : do ivarA=1,3*natom
203 443 : do ivarB=1,3*natom
204 396 : if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0) Nmatr(ivarA,ivarB)=one
205 396 : if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1) Nmatr(ivarA,ivarB)=one
206 438 : if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2) Nmatr(ivarA,ivarB)=one
207 : end do
208 : end do
209 :
210 : !DEBUG
211 : !do ivarA=1,3*natom
212 : !write(std_out,'(/)')
213 : !do ivarB=1,3*natom
214 : !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA)
215 : !end do
216 : !end do
217 : !ENDDEBUG
218 :
219 : !starting the pseudoinversing processes
220 : !then get the eigenvectors of the big matrix,give values to matrixBp
221 : ii1=1
222 47 : do ivarA=1,3*natom
223 266 : do ivarB=1,ivarA
224 219 : Bpmatr(1,ii1)=Nmatr(ivarB,ivarA)
225 261 : ii1=ii1+1
226 : end do
227 : end do
228 :
229 : !the imaginary part of the force matrix
230 224 : Bpmatr(2,:)=zero
231 : !then call the subroutines CHPEV and ZHPEV to get the eigenvectors
232 5 : call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier)
233 5 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
234 :
235 : !DEBUG
236 : !the eigenval and eigenvec
237 : !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors'
238 : !do ivarA=1,3*natom
239 : !write(std_out,'(/)')
240 : !write(std_out,'(es16.6)')eigvalp(ivarA)
241 : !end do
242 : !do ivarA=1,3*natom
243 : !write(std_out,'(/)')
244 : !do ivarB=1,3*natom
245 : !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA)
246 : !end do
247 : !end do
248 : !ENDDEBUG
249 :
250 : !then do the muplication to get the reduced matrix,in two steps
251 : !After this the force constant matrix is decouple in two bloks,
252 : !accoustic and optical ones
253 443 : Cpmatr(:,:)=zero
254 47 : do ivarA=1,3*natom
255 443 : do ivarB=1,3*natom
256 4542 : do ii1=1,3*natom
257 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*&
258 4500 : & Apmatr(ii1,ivarB)
259 : end do
260 : end do
261 : end do
262 :
263 443 : Apmatr(:,:)=zero
264 47 : do ivarA=1,3*natom
265 443 : do ivarB=1,3*natom
266 4542 : do ii1=1,3*natom
267 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*&
268 4500 : & eigvecp(1,ii1,ivarB)
269 : end do
270 : end do
271 : end do
272 :
273 : !DEBUG
274 : !the blok diago
275 : !write(std_out,'(/,a,/)')'Apmatr'
276 : !do ivarA=1,3*natom
277 : !write(std_out,'(/)')
278 : !do ivarB=1,3*natom
279 : !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB)
280 : !end do
281 : !end do
282 : !ENDDEBUG
283 :
284 : !check the last three eigenvalues whether too large or not
285 : ivarB=0
286 20 : do ivarA=3*natom-2,3*natom
287 20 : if (ABS(Apmatr(ivarA,ivarA))>tol6) ivarB=1
288 : end do
289 5 : if(ivarB==1)then
290 4 : write(message,'(a,a,a,a,a,a,a,a,a,a,3es16.6)')ch10,&
291 4 : & ' ddb_piezo : WARNING -',ch10,&
292 4 : & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,&
293 4 : & ' are too large at Gamma point',ch10,&
294 4 : & ' Increase cutoff energy or k-points sampling.',ch10,&
295 8 : & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom)
296 4 : call wrtout(std_out, message, 'COLL')
297 4 : call wrtout(iout,message,'COLL')
298 : end if
299 :
300 : !give the value of reduced matrix form Apmatr to Amatr
301 32 : do ivarA=1,3*natom-3
302 221 : do ivarB=1,3*natom-3
303 216 : Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)
304 : end do
305 : end do
306 : !now the reduced matrix is in the matrixA, the convert it
307 : !first give the give the value of matixB from matrixA
308 : ii1=1
309 32 : do ivarA=1,3*natom-3
310 140 : do ivarB=1,ivarA
311 108 : Bmatr(1,ii1)=Amatr(ivarB,ivarA)
312 135 : ii1=ii1+1
313 : end do
314 : end do
315 113 : Bmatr(2,:)=zero
316 :
317 : !call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues
318 5 : call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier)
319 5 : ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier)))
320 :
321 : !check the unstable phonon modes, if the first is negative then print warning message
322 5 : if(eigval(1)<-1.0*tol8)then
323 0 : write(message,'(a,a,a,a,a,a)') ch10,&
324 0 : & ' ddb_piezo : WARNING -',ch10,&
325 0 : & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,&
326 0 : & ' The system under calculation is physically unstable.'
327 0 : call wrtout(std_out, message, 'COLL')
328 0 : call wrtout(iout,message,'COLL')
329 : end if
330 :
331 : !do the matrix multiplication to get pseudoinverse inverse matrix
332 221 : Cmatr(:,:)=zero
333 221 : Amatr(:,:)=zero
334 32 : do ivarA=1,3*natom-3
335 32 : Cmatr(ivarA,ivarA)=one/eigval(ivarA)
336 : end do
337 :
338 32 : do ivarA=1,3*natom-3
339 221 : do ivarB=1,3*natom-3
340 1755 : do ii1=1,3*natom-3
341 : Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*&
342 1728 : & Cmatr(ii1,ivarB)
343 : end do
344 : end do
345 : end do
346 :
347 : !the second mulplication
348 221 : Cmatr(:,:)=zero
349 32 : do ivarA=1,3*natom-3
350 221 : do ivarB=1,3*natom-3
351 1755 : do ii1=1,3*natom-3
352 : Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+&
353 1728 : & Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1)
354 : end do
355 : end do
356 : end do
357 :
358 : !DEBUG
359 : !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix'
360 : !do ivarA=1,3*natom
361 : !write(std_out,'(/)')
362 : !do ivarB=1,3*natom
363 : !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB)
364 : !end do
365 : !end do
366 : !ENDDEBUG
367 :
368 : !so now the inverse of the reduced matrix is in the matrixC
369 : !now do another mulplication to get the pseudoinverse of the original
370 443 : Cpmatr(:,:)=zero
371 443 : Apmatr(:,:)=zero
372 32 : do ivarA=1,3*natom-3
373 221 : do ivarB=1,3*natom-3
374 216 : Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)
375 : end do
376 : end do
377 :
378 : !now times the eigvecp
379 47 : do ivarA=1,3*natom
380 443 : do ivarB=1,3*natom
381 4542 : do ii1=1,3*natom
382 4500 : Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*Cpmatr(ii1,ivarB)
383 : end do
384 : end do
385 : end do
386 :
387 443 : Cpmatr(:,:)=zero
388 47 : do ivarA=1,3*natom
389 443 : do ivarB=1,3*natom
390 4542 : do ii1=1,3*natom
391 : Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+&
392 4500 : & Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1)
393 : end do
394 : end do
395 : end do
396 :
397 : !now the inverse in in Cpmatr
398 443 : kmatrix(:,:)=Cpmatr(:,:)
399 : !transfer the inverse of k-matrix back to the k matrix
400 : !so now the inverse of k matrix is in the kmatrix
401 : !ending the part for pseudoinversing the K matrix
402 :
403 : !we still need the z-star matrix
404 20 : do idir1=1,3
405 : d2cart_relaxed(1,idir1,natom+2,idir1,natom+2,iblok)=&
406 20 : & d2cart_relaxed(1,idir1,natom+2,idir1,natom+2,iblok)-1.0_dp
407 : end do
408 :
409 20 : do idir1=1,3
410 65 : do idir2=1,3
411 150 : do ii1=1,2
412 : d2cart_relaxed(ii1,idir1,natom+2,idir2,natom+2,iblok)=&
413 135 : & d2cart_relaxed(ii1,idir1,natom+2,idir2,natom+2,iblok)/four_pi
414 : end do
415 : end do
416 : end do
417 :
418 20 : do ivarA=1,3
419 15 : idir1=ivarA
420 15 : ipert1=natom+2
421 62 : do ipert2=1,natom
422 183 : do idir2=1,3
423 126 : ivarB=idir2+3*(ipert2-1)
424 168 : zstar1(ivarA,ivarB)=d2cart_relaxed(1,idir1,ipert1,idir2,ipert2,iblok)
425 : end do
426 : end do
427 : end do
428 :
429 : !then get the inverse of the zstar1 for zstar2(3*natom,3)
430 47 : do ivarA=1,3*natom
431 173 : do ivarB=1,3
432 168 : zstar2(ivarA,ivarB)=zstar1(ivarB,ivarA)
433 : end do
434 : end do
435 : !the the matrix I need for the multiplication is in kmatrix and zstar2
436 :
437 : !the first matrix mulplication
438 299 : new1(:,:)=zero
439 35 : do ii1=1,6
440 287 : do ii2=1,3*natom
441 2658 : do ivarA=1,3*natom
442 : new1(ii1,ii2)=new1(ii1,ii2)+instrain(ivarA,ii1)*&
443 2628 : & kmatrix(ivarA,ii2)
444 : end do
445 : end do
446 : end do
447 :
448 : !do the second matrix mulplication
449 5 : piezo_correction(:,:)=zero
450 35 : do ii1=1,6
451 125 : do ii2=1,3
452 876 : do ivarA=1,3*natom
453 : piezo_correction(ii1,ii2)=piezo_correction(ii1,ii2)+&
454 846 : & new1(ii1,ivarA)* zstar2(ivarA,ii2)
455 : end do
456 : end do
457 : end do
458 :
459 : !then consider the volume and the change the unit form atomic to SI
460 35 : do ii1=1,6
461 125 : do ii2=1,3
462 90 : piezo_correction(ii1,ii2)= (piezo_correction(ii1,ii2) / ucvol)*AmuBohr2_Cm2
463 120 : piezo_relaxed(ii1,ii2)=piezo_clamped(ii1,ii2)+ piezo_correction(ii1,ii2)
464 : end do
465 : end do
466 : !end the calculation of piezoelectric constants
467 :
468 : !then print out the relaxed ion piezoelectric constants
469 :
470 5 : if(inp%piezoflag==2.or.inp%piezoflag==3 .or. inp%piezoflag==7)then
471 5 : if(inp%instrflag==0)then
472 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
473 0 : & ' WARNING: in order to get the piezoelectric tensor (relaxed ion), ',ch10,&
474 0 : & ' one needs information about internal strain ',ch10,&
475 0 : & ' one should set instrflag==1;',ch10,&
476 0 : & ' otherwise the program will continue but will give wrong values.'
477 0 : call wrtout(std_out,message,'COLL')
478 0 : call wrtout(iout,message,'COLL')
479 : end if
480 5 : write(message,'(3a)')ch10,' Proper piezoelectric constants (relaxed ion) (unit:c/m^2)',ch10
481 5 : call wrtout(std_out,message,'COLL')
482 35 : do ivarA=1,6
483 35 : write(std_out,'(3f16.8)')piezo_relaxed(ivarA,1),piezo_relaxed(ivarA,2),piezo_relaxed(ivarA,3)
484 : end do
485 :
486 5 : call wrtout(iout,message,'COLL')
487 5 : if (iwrite) then
488 35 : do ivarA=1,6
489 35 : write(iout,'(3f16.8)')piezo_relaxed(ivarA,1),piezo_relaxed(ivarA,2),piezo_relaxed(ivarA,3)
490 : end do
491 : end if
492 : end if
493 :
494 : !DEBUG
495 : !check the values of the relaxed ion dielectric tensor
496 : !write(message,'(a,a,a,a)')ch10,'debug the dielt tensor values ',&
497 : !& '(unit:c/m^2)',ch10
498 : !call wrtout(std_out,message,'COLL')
499 : !do ivarA=1,3
500 : !write(std_out,'(3f16.8)')dielt_rlx(ivarA,1),dielt_rlx(ivarA,2),dielt_rlx(ivarA,3)
501 : !end do
502 : !END DEBUG
503 :
504 : !DEBUG
505 : !print the relaxed ion elast tensor
506 : !write(message,'(a,a,a,a)')ch10,' debugElastic Tensor(relaxed ion)',&
507 : !& '(unit:10^2GP,VOIGT notation):',ch10
508 : !call wrtout(std_out,message,'COLL')
509 : !do ivarA=1,6
510 : !write(std_out,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,&
511 : !& elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,&
512 : !& elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp
513 : !end do
514 : !ENDDEBUG
515 :
516 : !Start to compute the piezoelectric d tensors
517 : !first initialize the d_tensor values
518 : !first make sure the elastic tensor is not zero
519 5 : d_tensor(:,:)=zero
520 5 : if(inp%elaflag>1)then
521 : ! then get the relaxed ion compliance tensor
522 5 : compliance(:,:)=elast(:,:)
523 5 : call matrginv(compliance,6,6)
524 35 : do ivarA=1,6
525 125 : do ivarB=1,3
526 660 : do ii1=1,6
527 : d_tensor(ivarA,ivarB)=d_tensor(ivarA,ivarB)+compliance(ivarA,ii1)*&
528 630 : & piezo_relaxed(ii1,ivarB)
529 : end do
530 : end do
531 : end do
532 : ! then convert in to the right unit pc/N
533 35 : do ivarA=1,6
534 125 : do ivarB=1,3
535 120 : d_tensor(ivarA,ivarB)=1000*d_tensor(ivarA,ivarB)
536 : end do
537 : end do
538 : end if
539 :
540 : !then print out the results of d tensor in log and output files
541 5 : if(inp%piezoflag==4 .or. inp%piezoflag==7)then
542 2 : if(inp%instrflag==0 .or. inp%elaflag==0 .or. inp%elaflag==1)then
543 0 : write(message,'(12a)' )ch10,&
544 0 : & ' WARNING:in order to get the piezoelectric d tensor(relaxed ion),', ch10,&
545 0 : & ' one needs the elastic tensor(relaxed ion) and piezoelectric e tensor',ch10,&
546 0 : & ' the latter needs the information of internal strain;',ch10,&
547 0 : & ' please check that both instrflag and elaflag are set to correct numbers',ch10,&
548 0 : & ' (elaflag= 2,3,4, or 5; instrflag=1)',ch10,&
549 0 : & ' otherwise the program will continue but give wrong values.'
550 0 : call wrtout(std_out,message,'COLL')
551 0 : call wrtout(iout,message,'COLL')
552 : end if
553 2 : write(message,'(3a)')ch10,' Piezoelectric d tensor (relaxed ion) (unit:pc/N)',ch10
554 2 : call wrtout(std_out,message,'COLL')
555 14 : do ivarA=1,6
556 14 : write(std_out,'(3f16.8)')d_tensor(ivarA,1),d_tensor(ivarA,2),d_tensor(ivarA,3)
557 : end do
558 2 : call wrtout(iout,message,'COLL')
559 2 : if (iwrite) then
560 14 : do ivarA=1,6
561 14 : write(iout,'(3f16.8)')d_tensor(ivarA,1),d_tensor(ivarA,2),d_tensor(ivarA,3)
562 : end do
563 : end if
564 : end if
565 : !end the part of piezoelectric d tensor (relaxed ion only).
566 :
567 : !then start to compute the piezoelectric g tensor
568 : !according to the equation, we first need to know the information
569 : !of the free-stress dielectric tensor
570 : !first make sure dielt_rlx exits, so we do not invert zero matrix
571 5 : dielt_stress(:,:)=zero
572 5 : g_tensor(:,:)=zero
573 5 : if(inp%dieflag>0)then
574 12 : do ivarA=1,3
575 39 : do ivarB=1,3
576 27 : dielt_stress(ivarA,ivarB)=zero
577 198 : do ii1=1,6
578 : dielt_stress(ivarA,ivarB)=dielt_stress(ivarA,ivarB)+&
579 189 : & piezo_relaxed(ii1,ivarA)*d_tensor(ii1,ivarB)
580 : end do
581 : end do
582 : end do
583 :
584 : ! then combine the relaxed ion(fixed strain) dielectric
585 : ! tensor and also restore the unit
586 12 : do ivarA=1,3
587 39 : do ivarB=1,3
588 : dielt_stress(ivarA,ivarB)=dielt_rlx(ivarA,ivarB)+&
589 36 : & dielt_stress(ivarA,ivarB)/(eps0*1.0e12)
590 : end do
591 : end do
592 :
593 : ! DEBUG
594 : ! write(message,'(a,a,a,a)')ch10,'debug the free stress dielectric tensor ',&
595 : ! & '(unit:pc/N)',ch10
596 : ! call wrtout(std_out,message,'COLL')
597 : ! do ivarA=1,3
598 : ! write(std_out,'(3f16.8)')dielt_stress(ivarA,1),dielt_stress(ivarA,2),&
599 : ! & dielt_stress(ivarA,3)
600 : ! end do
601 : ! ENDDEBUG
602 :
603 : ! then get the g tensor
604 : beta_tensor(:,:)=0
605 3 : beta_tensor(:,:)=dielt_stress(:,:)
606 :
607 3 : call matrginv(beta_tensor,3,3)
608 12 : do ivarA=1,3
609 66 : do ivarB=1,6
610 54 : g_tensor(ivarA,ivarB)=zero
611 225 : do ii1=1,3
612 : g_tensor(ivarA,ivarB)=g_tensor(ivarA,ivarB)+beta_tensor(ivarA,ii1)*&
613 216 : & d_tensor(ivarB,ii1)
614 : end do
615 : end do
616 : end do
617 : ! then restore the unit to be m^2/C
618 12 : do ivarA=1,3
619 66 : do ivarB=1,6
620 63 : g_tensor(ivarA,ivarB)=g_tensor(ivarA,ivarB)/(eps0*1.0e12)
621 : end do
622 : end do
623 : end if
624 : !then print out the final results of the g tensors(relaxed ion)
625 5 : if(inp%piezoflag==5 .or. inp%piezoflag==7)then
626 : if(inp%instrflag==0 .or. inp%elaflag==0&
627 : & .or. inp%elaflag==1 .or. inp%elaflag==0&
628 2 : & .or. inp%dieflag==2 .or. inp%dieflag==1)then
629 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
630 0 : & ' WARNING:in order to get the piezoelectric g tensor(relaxed ion),',ch10,&
631 0 : & ' need internal strain, dielectric(relaxed-ion) and elastic(realxed ion)',ch10,&
632 0 : & ' please set instrflag==1, elaflag==2,3,4 or 5, dieflag==3 or 4',ch10,&
633 0 : & ' otherwise the program will still continue but give wrong values.'
634 0 : call wrtout(std_out,message,'COLL')
635 0 : call wrtout(iout,message,'COLL')
636 : end if
637 :
638 2 : write(message,'(3a)')ch10,' Piezoelectric g tensor (relaxed ion) (unit:m^2/c)',ch10
639 2 : call wrtout(std_out,message,'COLL')
640 14 : do ivarA=1,6
641 14 : write(std_out,'(3f16.8)')g_tensor(1,ivarA),g_tensor(2,ivarA),g_tensor(3,ivarA)
642 : end do
643 2 : call wrtout(iout,message,'COLL')
644 2 : if (iwrite) then
645 14 : do ivarA=1,6
646 14 : write(iout,'(3f16.8)')g_tensor(1,ivarA),g_tensor(2,ivarA),g_tensor(3,ivarA)
647 : end do
648 : end if
649 : end if
650 : !end the part of piezoelectric g tensor (relaxed ion only).
651 :
652 : !then start the part for computation of h tensor
653 5 : h_tensor(:,:)=zero
654 : !first make sure the dielt_rlx is not zero in the memory
655 5 : if(inp%dieflag>0)then
656 : beta_tensor(:,:)=0
657 3 : beta_tensor(:,:)=dielt_rlx(:,:)
658 : ! write(std_out,*)' call matrginv 3, dielt_rlx(:,:)= ',dielt_rlx(:,:)
659 :
660 3 : call matrginv(beta_tensor,3,3)
661 12 : do ivarA=1,3
662 66 : do ivarB=1,6
663 54 : h_tensor(ivarA,ivarB)=zero
664 225 : do ii1=1,3
665 : h_tensor(ivarA,ivarB)=h_tensor(ivarA,ivarB)+beta_tensor(ivarA,ii1)*&
666 216 : & piezo_relaxed(ivarB,ii1)
667 : end do
668 : end do
669 : end do
670 : ! then restore the unit to be N/c
671 12 : do ivarA=1,3
672 66 : do ivarB=1,6
673 63 : h_tensor(ivarA,ivarB)=1000.0*(h_tensor(ivarA,ivarB)/(eps0*1.0e12))
674 : end do
675 : end do
676 : end if
677 : !then print out the final results of h tensors
678 5 : if(inp%piezoflag==6 .or. inp%piezoflag==7)then
679 2 : if(inp%instrflag==0 .or. inp%dieflag==1 .or. &
680 : & inp%dieflag==2)then
681 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
682 0 : & ' WARNING: in order to get the h tensor, ',ch10,&
683 0 : & ' one needs information about internal strain and dielectric(relaxed ion)',ch10,&
684 0 : & ' one should set dieflag==3 or 4 and instrflag==1;',ch10,&
685 0 : & ' otherwise the program will continue but give wrong values.'
686 0 : call wrtout(std_out,message,'COLL')
687 0 : call wrtout(iout,message,'COLL')
688 : end if
689 2 : write(message,'(3a)')ch10,' Piezoelectric h tensor (relaxed ion) (unit:GN/c)',ch10
690 2 : call wrtout(std_out,message,'COLL')
691 14 : do ivarA=1,6
692 14 : write(std_out,'(3f16.8)')h_tensor(1,ivarA),h_tensor(2,ivarA),h_tensor(3,ivarA)
693 : end do
694 2 : call wrtout(iout,message,'COLL')
695 2 : if (iwrite) then
696 14 : do ivarA=1,6
697 14 : write(iout,'(3f16.8)')h_tensor(1,ivarA),h_tensor(2,ivarA),h_tensor(3,ivarA)
698 : end do
699 : end if
700 : end if
701 : !end the part of piezoelectric h tensor (relaxed ion only).
702 :
703 : !print the free stress dielectric tensor
704 5 : if(inp%dieflag==4)then
705 2 : write(message, '(a,a)')ch10,'************************************************'
706 2 : call wrtout(std_out,message,'COLL')
707 2 : call wrtout(iout,message,'COLL')
708 2 : if(inp%instrflag==0 .or. inp%elaflag==0 .or. inp%elaflag==1)then
709 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
710 0 : & ' WARNING: in order to get the free stress dielectric tensor,',ch10,&
711 0 : & ' one needs internal strain and elastic (relaxed ion)', ch10,&
712 0 : & ' we need set elaflag==2,3,4 or 5 and instrflag==1.',ch10,&
713 0 : & ' otherwise the program may continue but give wrong and nonsense values.'
714 0 : call wrtout(std_out,message,'COLL')
715 0 : call wrtout(iout,message,'COLL')
716 : end if
717 2 : write(message,'(a,a,a)')ch10,' Free stress dielectric tensor (dimensionless)',ch10
718 2 : call wrtout(std_out,message,'COLL')
719 8 : do ivarA=1,3
720 8 : write(std_out,'(3f16.8)')dielt_stress(ivarA,1),dielt_stress(ivarA,2),dielt_stress(ivarA,3)
721 : end do
722 2 : call wrtout(iout,message,'COLL')
723 2 : if (iwrite) then
724 8 : do ivarA=1,3
725 8 : write(iout,'(3f16.8)')dielt_stress(ivarA,1),dielt_stress(ivarA,2),dielt_stress(ivarA,3)
726 : end do
727 : end if
728 : end if
729 : !end the part of printing out the free stress dielectric tensor
730 :
731 : !then print out the fixed displacement elastic tensor
732 5 : elast_dis = zero; compliance_dis = zero
733 :
734 5 : if(inp%elaflag==4 .and. inp%dieflag>0)then
735 : if(inp%instrflag==0 .or. inp%dieflag==1 .or. &
736 2 : & inp%dieflag==2 .or. inp%dieflag==0)then
737 0 : write(message,'(a,a,a,a,a,a,a,a)' )ch10,&
738 0 : & ' WARNING: in order to get the elatic(fixed D field) tensor, ',ch10,&
739 0 : & ' one needs information about internal strain and dielectric(relaxed ion)',ch10,&
740 0 : & ' one should set dieflag==3 or 4 and instrflag==1;',ch10,&
741 0 : & ' otherwise the program will continue but give wrong values.'
742 0 : call wrtout(std_out,message,'COLL')
743 0 : call wrtout(iout,message,'COLL')
744 : end if
745 :
746 : !Then begin the computation of the fixed displacement elastic and compliance tensor(relaxed ion)
747 14 : do ivarA=1,6
748 86 : do ivarB=1,6
749 72 : elast_dis(ivarA,ivarB)=zero
750 300 : do ii1=1,3
751 : elast_dis(ivarA,ivarB)=elast_dis(ivarA,ivarB)+&
752 288 : & h_tensor(ii1,ivarA)*piezo_relaxed(ivarB,ii1)
753 : end do
754 : end do
755 : end do
756 : !Then should add the relaxed ion fixed E-field values
757 14 : do ivarA=1,6
758 86 : do ivarB=1,6
759 84 : elast_dis(ivarA,ivarB)=elast_dis(ivarA,ivarB)+elast(ivarA,ivarB)
760 : end do
761 : end do
762 :
763 2 : write(message, '(a,a)')ch10,'************************************************'
764 2 : call wrtout(std_out,message,'COLL')
765 2 : call wrtout(iout,message,'COLL')
766 2 : write(message,'(5a)')ch10,&
767 2 : & ' Elastic Tensor (relaxed ion) (unit:10^2GP)',ch10,&
768 4 : & ' (at fixed displacement field boundary condition)',ch10
769 2 : call wrtout(std_out,message,'COLL')
770 14 : do ivarA=1,6
771 12 : write(std_out,'(6f12.7)')elast_dis(ivarA,1)/100.00_dp,elast_dis(ivarA,2)/100.00_dp,&
772 12 : & elast_dis(ivarA,3)/100.00_dp,elast_dis(ivarA,4)/100.00_dp,&
773 26 : & elast_dis(ivarA,5)/100.00_dp,elast_dis(ivarA,6)/100.00_dp
774 : end do
775 2 : call wrtout(iout,message,'COLL')
776 2 : if (iwrite) then
777 14 : do ivarA=1,6
778 12 : write(iout,'(6f12.7)')elast_dis(ivarA,1)/100.00_dp,elast_dis(ivarA,2)/100.00_dp,&
779 12 : & elast_dis(ivarA,3)/100.00_dp,elast_dis(ivarA,4)/100.00_dp,&
780 26 : & elast_dis(ivarA,5)/100.00_dp,elast_dis(ivarA,6)/100.00_dp
781 : end do
782 : end if
783 : ! then invert the above to get the corresponding compliance tensor
784 : compliance_dis(:,:)=0
785 2 : compliance_dis(:,:)=elast_dis(:,:)
786 :
787 2 : call matrginv(compliance_dis,6,6)
788 : ! then print out the compliance tensor at fixed displacement field
789 2 : write(message,'(5a)')ch10,&
790 2 : & ' Compliance Tensor (relaxed ion) (unit: 10^-2(GP)^-1)',ch10,&
791 4 : & ' (at fixed displacement field boundary condition)',ch10
792 2 : call wrtout(std_out,message,'COLL')
793 14 : do ivarB=1,6
794 12 : write(std_out,'(6f12.7)')compliance_dis(ivarB,1)*100.00_dp,&
795 12 : & compliance_dis(ivarB,2)*100.00_dp,&
796 12 : & compliance_dis(ivarB,3)*100.00_dp,compliance_dis(ivarB,4)*100.00_dp,&
797 12 : & compliance_dis(ivarB,5)*100.00_dp,&
798 26 : & compliance_dis(ivarB,6)*100.00_dp
799 : end do
800 2 : call wrtout(iout,message,'COLL')
801 :
802 2 : if (iwrite) then
803 14 : do ivarB=1,6
804 12 : write(iout,'(6f12.7)')compliance_dis(ivarB,1)*100.00,&
805 12 : & compliance_dis(ivarB,2)*100.00_dp,&
806 12 : & compliance_dis(ivarB,3)*100.00_dp,compliance_dis(ivarB,4)*100.00_dp,&
807 12 : & compliance_dis(ivarB,5)*100.00_dp,&
808 26 : & compliance_dis(ivarB,6)*100.00_dp
809 : end do
810 : end if
811 : end if
812 : !end if the elaflag==4 for the fixed didplacement field elastic tensor
813 : !end the part for computation of elastic at fixed displacement field
814 :
815 : ! write tensors to netcdf file.
816 5 : if (ncid /= nctk_noid) then
817 : ncerr = nctk_def_arrays(ncid, [ &
818 : nctkarr_t("piezo_clamped_ion", "dp", "six, three"), &
819 : nctkarr_t("piezo_relaxed_ion", "dp", "six, three"), &
820 : nctkarr_t("d_tensor_relaxed_ion", "dp", "six, three"), &
821 : nctkarr_t("g_tensor_relaxed_ion", "dp", "three, six"), &
822 : nctkarr_t("h_tensor_relaxed_ion", "dp", "three, six"), &
823 : nctkarr_t("free_stress_dielectric_tensor", "dp", "three, three"), &
824 : nctkarr_t("elastic_tensor_relaxed_ion_fixed_D", "dp", "six, six"), &
825 : nctkarr_t("compliance_tensor_relaxed_ion_fixed_D", "dp", "six, six")], &
826 45 : defmode=.True.)
827 5 : NCF_CHECK(ncerr)
828 : ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: &
829 25 : "elaflag", "piezoflag", "instrflag", "dieflag"])
830 5 : NCF_CHECK(ncerr)
831 :
832 5 : NCF_CHECK(nctk_set_datamode(ncid))
833 : ncerr = nctk_write_iscalars(ncid, [character(len=nctk_slen) :: &
834 : "elaflag", "piezoflag", "instrflag", "dieflag"], &
835 45 : [inp%elaflag, inp%piezoflag, inp%instrflag, inp%dieflag])
836 5 : NCF_CHECK(ncerr)
837 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "piezo_clamped_ion"), piezo_clamped))
838 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "piezo_relaxed_ion"), piezo_relaxed))
839 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "d_tensor_relaxed_ion"), d_tensor))
840 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "g_tensor_relaxed_ion"), g_tensor))
841 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "h_tensor_relaxed_ion"), h_tensor))
842 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "free_stress_dielectric_tensor"), dielt_stress))
843 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "elastic_tensor_relaxed_ion_fixed_D"), elast_dis))
844 5 : NCF_CHECK(nf90_put_var(ncid, nctk_idname(ncid, "compliance_tensor_relaxed_ion_fixed_D"), compliance_dis))
845 : end if
846 :
847 5 : end subroutine ddb_piezo
848 : !!***
849 :
850 : end module m_ddb_piezo
851 : !!***
|
