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