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