TABLE OF CONTENTS


ABINIT/ddb_piezo [ Functions ]

[ Top ] [ Functions ]

NAME

 ddb_piezo

FUNCTION

 Get the piezoelectric tensor (e-tensor), both clamped ion and relaxed ion;
 Compute physical(relaxed ion) piezoeletric (d, g, h) tensors;
 Compute relaxed ion and free stress dielectric tensor;
 Compute relaxed ion elastic and compliance tensors under fixed
 displacement field boundary conditions.

INPUTS

 inp= (derived datatype) contains all the input variables
 blkval(2,3,mpert,3,mpert,nblok)=
   second derivatives of total energy with respect to electric fields
   atom displacements,strain,...... all in cartesian coordinates
 dielt_rlx=relaxed ion dielectric tensor
 iblok= bolk number in DDB file contains 2 derivative of energy
 instrain=force response internal strain tensor
 iout=out file number
 mpert=maximum number of ipert
 natom=number of atoms in unit cell
 nblok=number of total bloks in DDB file
 ucvol=unit cell volume
 ncid=the id of the open NetCDF file. Set to nctk_noid if netcdf output is not wanted.

OUTPUT

 piezo = piezoelectric tensor

NOTES

 The elastic (compliance) tensors calculated here are under fixed D-field boundary
 condition, which include piezoelectric corrections to the elastic (compliance)
 tensors calculated in ddb_elast.F90 whose boundary condition is fixed E-field.

PARENTS

      anaddb

CHILDREN

      matrginv,wrtout,zhpev

SOURCE

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

ABINIT/m_ddb_piezo [ Modules ]

[ Top ] [ Modules ]

NAME

  m_ddb_piezo

FUNCTION

COPYRIGHT

  Copyright (C) 1999-2018 ABINIT group (XW)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_ddb_piezo
27 
28  use defs_basis
29  use m_abicore
30  use m_errors
31  use m_nctk
32 #ifdef HAVE_NETCDF
33  use netcdf
34 #endif
35 
36  use m_fstrings,       only : sjoin, itoa
37  use m_hide_lapack,    only : matrginv
38  use m_anaddb_dataset, only : anaddb_dataset_type
39 
40  implicit none
41 
42  private