Branch data 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 : : !!***
|
