TABLE OF CONTENTS
ABINIT/ddb_elast [ Functions ]
NAME
ddb_elast
FUNCTION
Get the elastic and compliance tensors, both clamped ion and relaxed ion, under the fixed electric field boundary condition; in which realxed ion tensors can generate two output tensors one is conventional, the other considers the sress correction.
INPUTS
inp= (derived datatype) contains all the input variables crystal<crystal_t>=Info on crystalline structure. blkval(2,3,mpert,3,mpert,nblok)= second derivatives of total energy with respect to electric fields atom displacements,strain,...... all in cartesian coordinates d2asr= ASR correction to the dynamical matrix at Gamma iblok= bolk number in DDB file iblok_stress= blok number which contain stress tensor 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 ncid=NC file handle (open in the caller)
OUTPUT
elast=relaxed-ion elastic tensor(without stress correction) (6*6) in Voigt notation elast_clamped=clamped-ion elastic tensor(without stress correction) (6*6) in Voigt notation elast_stress=relaxed-ion elastic tensor(with stress correction) (6*6) in Voigt notation
NOTES
The elastic (compliance) tensors calculated here are under boundary conditions of fixed Electric Field, different from those in ddb_piezo.F90 which are under fixed Displacement Field and incorporate piezoelectric corrections.
SOURCE
89 subroutine ddb_elast(inp,crystal,blkval,compl,compl_clamped,compl_stress,d2asr,& 90 & elast,elast_clamped,elast_stress,iblok,iblok_stress,& 91 & instrain,iout,mpert,natom,nblok,ncid) 92 93 !Arguments ------------------------------------------- 94 !scalars 95 integer,intent(in) :: iblok,iblok_stress,iout,mpert,natom,nblok,ncid 96 !integer,intent(in) :: msize 97 type(crystal_t),intent(in) :: crystal 98 type(anaddb_dataset_type),intent(in) :: inp 99 !arrays 100 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),instrain(3*natom,6) 101 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 102 real(dp),intent(out) :: compl(6,6), compl_clamped(6,6),compl_stress(6,6) 103 real(dp),intent(out) :: elast(6,6), elast_clamped(6,6),elast_stress(6,6) 104 105 !Local variables------------------------------------ 106 !scalars 107 integer :: ier,ii1,ii2,ipert1,ipert2,ivarA,ivarB,ncerr 108 real(dp) :: ucvol 109 logical :: iwrite 110 character(len=500) :: message 111 !arrays 112 real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom) 113 real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2) 114 real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3) 115 real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom) 116 real(dp) :: compl_relaxed(6,6),eigval(3*natom-3) 117 real(dp) :: eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3) 118 real(dp) :: eigvecp(2,3*natom,3*natom),elast_relaxed(6,6) 119 real(dp) :: kmatrix(3*natom,3*natom),new1(6,3*natom) 120 real(dp) :: new2(6,6),stress(6),zhpev1(2,2*3*natom-4),zhpev1p(2,2*3*natom-1) 121 real(dp) :: zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2) 122 real(dp) :: d2cart(2,3*natom,3*natom) 123 124 !*************************************************************************** 125 compl = zero; compl_clamped = zero; compl_stress = zero 126 elast = zero; elast_clamped = zero; elast_stress = zero 127 128 ucvol = crystal%ucvol 129 iwrite = iout > 0 130 131 !extraction of the elastic constants from the blkvals 132 133 do ivarA=1,6 134 do ivarB=1,6 135 ! because the elastic constant is 6*6, 136 ! so we should judge if the idir is larger than 3 or not 137 if(ivarA>3) then 138 ii1=ivarA-3 139 ipert1=natom+4 !for the shear modulus 140 else if(ivarA<=3) then 141 ii1=ivarA 142 ipert1=natom+3 !for the diagonal part 143 end if 144 if(ivarB>3) then 145 ii2=ivarB-3 146 ipert2=natom+4 !for the shear modulus 147 else if(ivarB<=3) then 148 ii2=ivarB 149 ipert2=natom+3 !for the diagonal part 150 end if 151 elast(ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok) 152 end do 153 end do 154 155 !then consider the volume, because the unit above is in 156 !Hartree, in fact the elastic constant should be in 157 !the units of presure, the energy/volume 158 !And then transform the unit to si unit using GPa 159 !from Hartree/Bohr^3 160 161 do ivarA=1,6 162 do ivarB=1,6 163 elast(ivarA,ivarB)=(elast(ivarA,ivarB)/ucvol)*HaBohr3_GPa 164 end do 165 end do 166 167 !then should consider the two situations:clamped and relaxed 168 !ions respectively,give the initial value of elast_clamped 169 elast_clamped(:,:)=elast(:,:) 170 elast_relaxed(:,:)=elast(:,:) 171 172 !then do the matrix mulplication of instrain*K*instrain to get the 173 !correction of the relaxed ion quantities, in case natom/=1 174 175 if( (inp%elaflag==2 .or. inp%elaflag==3& 176 & .or. inp%elaflag==4 .or. inp%elaflag==5) .and. natom/=1 )then 177 ! extracting force matrix at gamma 178 d2cart = zero 179 do ipert1=1,natom 180 do ii1=1,3 181 ivarA=ii1+3*(ipert1-1) 182 do ipert2=1,natom 183 do ii2=1,3 184 ivarB=ii2+3*(ipert2-1) 185 d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok) 186 end do 187 end do 188 end do 189 end do 190 191 ! Eventually impose the acoustic sum rule 192 ! FIXME: this might depend on ifcflag: impose that it is 0 or generalize 193 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 194 call asria_corr(inp%asr,d2asr,d2cart,natom,natom) 195 kmatrix = d2cart(1,:,:) 196 197 ! DEBUG 198 ! write(std_out,'(/,a,/)')'the k matrix before inverse' 199 ! do ii1=1,3*natom 200 ! write(std_out,'(6es16.3)')kmatrix(ii1,1),kmatrix(ii1,2),kmatrix(ii1,3),& 201 ! & kmatrix(ii1,4),kmatrix(ii1,5),kmatrix(ii1,6) 202 ! end do 203 ! ENDDEBUG 204 205 ! according to formula, invert the kmatrix(3natom,3natom) 206 Apmatr(:,:)=kmatrix(:,:) 207 208 ! NOTE: MJV 13/3/2011 This is just the 3x3 unit matrix copied throughout the dynamical matrix 209 Nmatr(:,:)=zero 210 do ivarA=1,3*natom 211 do ivarB=1,3*natom 212 if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0) Nmatr(ivarA,ivarB)=one 213 if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1) Nmatr(ivarA,ivarB)=one 214 if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2) Nmatr(ivarA,ivarB)=one 215 end do 216 end do 217 218 ! DEBUG 219 ! The k matrix is not the inverse here - it has not been changed! 220 ! write(std_out,'(/,a,/)')'the direct inverse of the Kmatrix' 221 ! do ivarA=1,3*natom 222 ! write(std_out,'(/)') 223 ! do ivarB=1,3*natom 224 ! write(std_out,'(es16.6)')kmatrix(ivarA,ivarB) 225 ! end do 226 ! end do 227 ! ENDDEBUG 228 229 230 ! starting the pseudo-inverse processes 231 ! then get the eigenvectors of the big matrix,give values to matrixBp 232 ! Pack the Nmatr matrix in Hermitian form 233 ii1=1 234 do ivarA=1,3*natom 235 do ivarB=1,ivarA 236 Bpmatr(1,ii1)=Nmatr(ivarB,ivarA) 237 ii1=ii1+1 238 end do 239 end do 240 Bpmatr(2,:)=zero !the imaginary part of the force matrix 241 ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors 242 ! NOTE: MJV there is a huge indeterminacy in this matrix, which has all identical 3x3 block lines 243 ! this means the orientation of the 0-eigenvalue eigenvectors is kind of random... 244 ! Is the usage just to get out the translational modes? We know what the eigenvectors look like already! 245 ! The translational modes are the last 3 with eigenvalue 6 246 ! 247 call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier) 248 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 249 250 ! DEBUG 251 ! the eigenval and eigenvec 252 ! write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors' 253 ! do ivarA=1,3*natom 254 ! write(std_out,'(/)') 255 ! write(std_out,'(es16.6)')eigvalp(ivarA) 256 ! end do 257 ! do ivarA=1,3*natom 258 ! write(std_out,'(/)') 259 ! do ivarB=1,3*natom 260 ! write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA) 261 ! end do 262 ! end do 263 ! ENDDEBUG 264 265 ! do the multiplication to get the reduced matrix,in two steps 266 ! rotate to eigenbasis constructed above to isolate acoustic modes 267 Cpmatr(:,:)=zero 268 do ivarA=1,3*natom 269 do ivarB=1,3*natom 270 do ii1=1,3*natom 271 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*& 272 & Apmatr(ii1,ivarB) 273 end do 274 end do 275 end do 276 277 Apmatr(:,:)=zero 278 do ivarA=1,3*natom 279 do ivarB=1,3*natom 280 do ii1=1,3*natom 281 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*& 282 & eigvecp(1,ii1,ivarB) 283 end do 284 end do 285 end do 286 287 ! DEBUG 288 ! the blok diagonal parts 289 ! write(std_out,'(/,a,/)')'Apmatr' 290 ! do ivarA=1,3*natom 291 ! write(std_out,'(/)') 292 ! do ivarB=1,3*natom 293 ! write(std_out,'(es16.6)')Apmatr(ivarA,ivarB) 294 ! end do 295 ! end do 296 ! ENDDEBUG 297 298 299 ! check the last three eigenvalues whether too large 300 ivarB=0 301 do ivarA=3*natom-2,3*natom 302 if (ABS(Apmatr(ivarA,ivarA))>tol6) ivarB=1 303 end do 304 if(ivarB==1)then 305 write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,& 306 & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,& 307 & ' are too large at Gamma point',ch10,& 308 & ' increase cutoff energy or k-points sampling.',ch10,& 309 & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),& 310 & Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom) 311 ABI_WARNING(message) 312 call wrtout(iout,message,'COLL') 313 end if 314 ! then give the value of reduced matrix form Apmatr to Amatr 315 do ivarA=1,3*natom-3 316 do ivarB=1,3*natom-3 317 Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB) 318 end do 319 end do 320 321 ! now the reduced matrix is in the Amatr, the convert it 322 ! first give the give the value of Bmatr from Amatr 323 ii1=1 324 do ivarA=1,3*natom-3 325 do ivarB=1,ivarA 326 Bmatr(1,ii1)=Amatr(ivarB,ivarA) 327 ii1=ii1+1 328 end do 329 end do 330 Bmatr(2,:)=zero 331 ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues 332 call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier) 333 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 334 335 ! check the unstable phonon modes, if the first is negative then print warning message 336 if(eigval(1)<-1.0*tol8)then 337 write(message,'(a,a,a,a)') ch10,& 338 & 'Unstable eigenvalue detected in force constant matrix at Gamma point.',ch10,& 339 & 'The system under calculation is physically unstable.' 340 ABI_WARNING(message) 341 call wrtout(iout,message,'COLL') 342 end if 343 344 ! the do the matrix muplication to get pseudoinverse inverse matrix 345 Cmatr(:,:)=zero 346 Amatr(:,:)=zero 347 do ivarA=1,3*natom-3 348 Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA) 349 end do 350 do ivarA=1,3*natom-3 351 do ivarB=1,3*natom-3 352 do ii1=1,3*natom-3 353 Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*& 354 & Cmatr(ii1,ivarB) 355 end do 356 end do 357 end do 358 359 ! then the second mulplication 360 Cmatr(:,:)=zero 361 do ivarA=1,3*natom-3 362 do ivarB=1,3*natom-3 363 do ii1=1,3*natom-3 364 Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+& 365 & Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1) 366 end do 367 end do 368 end do 369 370 ! DEBUG 371 ! write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix' 372 ! do ivarA=1,3*natom 373 ! write(std_out,'(/)') 374 ! do ivarB=1,3*natom 375 ! write(std_out,'(es16.6)')Cmatr(ivarA,ivarB) 376 ! end do 377 ! end do 378 ! ENDDEBUG 379 380 ! so now the inverse of the reduced matrix is in the matrixC 381 ! now do another mulplication to get the pseudoinverse of the original 382 Cpmatr(:,:)=zero 383 Apmatr(:,:)=zero 384 do ivarA=1,3*natom-3 385 do ivarB=1,3*natom-3 386 Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB) 387 end do 388 end do 389 390 ! times the eigvecp 391 do ivarA=1,3*natom 392 do ivarB=1,3*natom 393 do ii1=1,3*natom 394 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*& 395 & Cpmatr(ii1,ivarB) 396 end do 397 end do 398 end do 399 Cpmatr(:,:)=zero 400 do ivarA=1,3*natom 401 do ivarB=1,3*natom 402 do ii1=1,3*natom 403 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+& 404 & Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1) 405 end do 406 end do 407 end do 408 409 ! now the inverse in in Cpmatr 410 kmatrix(:,:)=Cpmatr(:,:) 411 ! transfer the inverse of k-matrix back to the k matrix 412 ! so now the inverse of k matrix is in the kmatrix 413 ! ending the part for pseudoinversing the K matrix 414 ! then do the first matrix mulplication 415 new1(:,:)=zero 416 do ii1=1,6 417 do ii2=1,3*natom 418 do ivarA=1,3*natom 419 new1(ii1,ii2)=new1(ii1,ii2)+instrain(ivarA,ii1)*& 420 & kmatrix(ivarA,ii2) 421 end do 422 end do 423 end do 424 ! then do the second matrix mulplication, and change the value of kmatrix 425 new2(:,:)=zero 426 do ii1=1,6 427 do ii2=1,6 428 do ivarB=1,3*natom 429 new2(ii1,ii2)=new2(ii1,ii2)+new1(ii1,ivarB)*& 430 & instrain(ivarB,ii2) 431 end do 432 end do 433 end do 434 ! then finish the matrix mupl., consider the unit cellvolume 435 ! and the unit change next step 436 do ivarA=1,6 437 do ivarB=1,6 438 new2(ivarA,ivarB)=(new2(ivarA,ivarB)/ucvol)*HaBohr3_GPa 439 end do 440 end do 441 ! then the relaxed one should be the previous one minus the new2 element 442 do ivarA=1,6 443 do ivarB=1,6 444 elast_relaxed(ivarA,ivarB)=elast_relaxed(ivarA,ivarB)-& 445 & new2(ivarA,ivarB) 446 end do 447 end do 448 end if 449 !the above end if end if for elaflag=2 or elafalg=3 or elafalg=4, 450 !or elafalg=5 in line 125 451 452 !DEBUG 453 !write(std_out,'(/,a,/)')'debug the unit cell volume' 454 !write(std_out,'(2es16.6)')ucvol,HaBohr3_GPa 455 !ENDDEBUG 456 457 !then give the initial value of the compl_relaxed(6,6) 458 compl_relaxed(:,:)=elast_relaxed(:,:) 459 460 !******************************************************************* 461 if(inp%elaflag==1.or. inp%elaflag==3)then 462 ! print out the clamped-ion elastic constants to output file 463 write(message,'(3a)')ch10,' Elastic Tensor (clamped ion) (unit:10^2GP):',ch10 464 call wrtout(std_out,message,'COLL') 465 do ivarA=1,6 466 write(std_out,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,& 467 & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,& 468 & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp 469 end do 470 471 call wrtout(iout,message,'COLL') 472 if (iwrite) then 473 do ivarA=1,6 474 write(iout,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,& 475 & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,& 476 & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp 477 end do 478 end if 479 end if 480 481 if(inp%elaflag==2.or.inp%elaflag==3& 482 & .or. inp%elaflag==4.or. inp%elaflag==5)then 483 if(inp%instrflag==0)then 484 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 485 & 'in order to get the elastic tensor(relaxed ion), ',ch10,& 486 & 'one needs information about internal strain ',ch10,& 487 & 'one should set instrflag==1;',ch10,& 488 & 'otherwise the program will continue but give wrong values.' 489 ABI_WARNING(message) 490 call wrtout(iout,message,'COLL') 491 end if 492 493 write(message,'(5a)')ch10,& 494 & ' Elastic Tensor (relaxed ion) (unit:10^2GP):',ch10,& 495 & ' (at fixed electric field boundary condition)',ch10 496 call wrtout(std_out,message,'COLL') 497 do ivarA=1,6 498 write(std_out,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,& 499 & elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,& 500 & elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,& 501 & elast_relaxed(ivarA,6)/100.00_dp 502 end do 503 504 call wrtout(iout,message,'COLL') 505 if (iwrite) then 506 do ivarA=1,6 507 write(iout,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,& 508 & elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,& 509 & elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,& 510 & elast_relaxed(ivarA,6)/100.00_dp 511 end do 512 end if 513 end if 514 515 !then print the corresponding compliances 516 517 if(inp%elaflag==1.or.inp%elaflag==3)then 518 ! compl(:,:)=elast_clamped(:,:) !convert the elastic tensor 519 compl_clamped(:,:)=elast_clamped(:,:) 520 call matrginv(compl_clamped,6,6) 521 write(message,'(a,a,a)')ch10,' Compliance Tensor (clamped ion) (unit: 10^-2GP^-1):',ch10 522 call wrtout(std_out,message,'COLL') 523 524 do ivarB=1,6 525 write(std_out,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,& 526 & compl_clamped(ivarB,2)*100.00_dp,& 527 & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,& 528 & compl_clamped(ivarB,5)*100.00_dp,& 529 & compl_clamped(ivarB,6)*100.00_dp 530 end do 531 532 call wrtout(iout,message,'COLL') 533 534 if (iwrite) then 535 do ivarB=1,6 536 write(iout,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,& 537 & compl_clamped(ivarB,2)*100.00_dp,& 538 & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,& 539 & compl_clamped(ivarB,5)*100.00_dp,& 540 & compl_clamped(ivarB,6)*100.00_dp 541 end do 542 end if 543 end if 544 545 if(inp%elaflag==2.or.inp%elaflag==3& 546 & .or. inp%elaflag==4 .or. inp%elaflag==5)then 547 ! compl(:,:)=elast_relaxed(:,:) 548 call matrginv(compl_relaxed,6,6) 549 if(inp%instrflag==0)then 550 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 551 & 'in order to get the compliance tensor(relaxed ion), ',ch10,& 552 & 'one needs information about internal strain ',ch10,& 553 & 'one should set instrflag==1;',ch10,& 554 & 'otherwise the program will continue but give wrong values.' 555 ABI_WARNING(message) 556 call wrtout(iout,message,'COLL') 557 end if 558 write(message,'(5a)')ch10,& 559 & ' Compliance Tensor (relaxed ion) (unit: 10^-2GP^-1):',ch10,& 560 & ' (at fixed electric field boundary condition)',ch10 561 call wrtout(std_out,message,'COLL') 562 563 do ivarB=1,6 564 write(std_out,'(6f12.7)')compl_relaxed(ivarB,1)*100.00_dp,& 565 & compl_relaxed(ivarB,2)*100.00_dp,& 566 & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,& 567 & compl_relaxed(ivarB,5)*100.00_dp,& 568 & compl_relaxed(ivarB,6)*100.00_dp 569 end do 570 call wrtout(iout,message,'COLL') 571 572 if (iwrite) then 573 do ivarB=1,6 574 write(iout,'(6f12.7)')compl_relaxed(ivarB,1)*100.00,& 575 & compl_relaxed(ivarB,2)*100.00_dp,& 576 & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,& 577 & compl_relaxed(ivarB,5)*100.00_dp,& 578 & compl_relaxed(ivarB,6)*100.00_dp 579 end do 580 end if 581 end if 582 583 !befor the end , make sure the tensor elast(6,6) 584 !will have the relaxed ion values and tensor elast_clamped has the clamped ion 585 !values, and similarily for the corresponding compliance tensors 586 elast_clamped(:,:)=elast(:,:) 587 elast(:,:)=elast_relaxed(:,:) 588 compl(:,:)=compl_relaxed(:,:) 589 590 !begin the part of computing stress corrected elastic tensors 591 if(inp%elaflag==5)then 592 593 ! DEBUG 594 ! check the iblok number of first derivative of energy 595 ! write(std_out,'(/,a,/)')'iblok number at 8:00Pm' 596 ! write(std_out,'(i)')iblok_stress 597 ! write(std_out,'(a,f12.7)')'the total energy', blkval(1,1,1) 598 ! write(std_out,*)'',blkval(1,:,:,:,:,iblok_stress) 599 ! write(std_out,*)'',blkval(1,:,7,1,1,iblok_stress) 600 ! ENDDEBUG 601 602 ! firts give the corect stress values diagonal parts 603 stress(1)=blkval(1,1,natom+3,1,1,iblok_stress) 604 stress(2)=blkval(1,2,natom+3,1,1,iblok_stress) 605 stress(3)=blkval(1,3,natom+3,1,1,iblok_stress) 606 ! the shear parts 607 stress(4)=blkval(1,1,natom+4,1,1,iblok_stress) 608 stress(5)=blkval(1,2,natom+4,1,1,iblok_stress) 609 stress(6)=blkval(1,3,natom+4,1,1,iblok_stress) 610 ! then convert the unit from atomic to the GPa unit 611 do ivarA=1,6 612 stress(ivarA)=stress(ivarA)*HaBohr3_GPa 613 end do 614 ! give the initial values of elast_stress tensor 615 elast_stress(:,:)=elast_relaxed(:,:) 616 ! notice that only the first three rows need to be corrected 617 do ivarA=1,3 618 do ivarB=1,6 619 elast_stress(ivarA,ivarB)=elast_stress(ivarA,ivarB)-stress(ivarB) 620 end do 621 end do 622 ! then compute the values of compliance tensor with stress correction 623 compl_stress(:,:)=elast_stress(:,:) 624 call matrginv(compl_stress,6,6) 625 ! then print out the results of stress corrected elastic and compliance tensors 626 if(inp%instrflag==0)then 627 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 628 & 'In order to get the elastic tensor (relaxed ion with stress correction), ',ch10,& 629 & 'one needs information about internal strain ',ch10,& 630 & 'one should set instrflag==1;',ch10,& 631 & 'otherwise the program will continue but give wrong values.' 632 ABI_WARNING(message) 633 call wrtout(iout,message,'COLL') 634 end if 635 write(message,'(5a)')ch10,& 636 & ' Elastic Tensor (relaxed ion with stress corrected) (unit:10^2GP)',ch10,& 637 & ' (at fixed electric field boundary condition)',ch10 638 call wrtout(std_out,message,'COLL') 639 call wrtout(iout,message,'COLL') 640 do ivarA=1,6 641 write(std_out,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,& 642 & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,& 643 & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp 644 end do 645 if (iwrite) then 646 do ivarA=1,6 647 write(iout,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,& 648 & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,& 649 & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp 650 end do 651 end if 652 653 ! then the complinace tensors with stress correction 654 write(message,'(5a)')ch10,& 655 & ' Compliance Tensor (relaxed ion with stress correction) (unit: 10^-2(GP)^-1):',ch10,& 656 & ' (at fixed electric field boundary condition)',ch10 657 call wrtout(std_out,message,'COLL') 658 call wrtout(iout,message,'COLL') 659 do ivarB=1,6 660 write(std_out,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,& 661 & compl_stress(ivarB,2)*100.00_dp,& 662 & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,& 663 & compl_stress(ivarB,5)*100.00_dp,& 664 & compl_stress(ivarB,6)*100.00_dp 665 end do 666 if (iwrite) then 667 do ivarB=1,6 668 write(iout,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,& 669 & compl_stress(ivarB,2)*100.00_dp,& 670 & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,& 671 & compl_stress(ivarB,5)*100.00_dp,& 672 & compl_stress(ivarB,6)*100.00_dp 673 end do 674 end if 675 end if 676 !end the if 510th line 677 !end the part of stress corrected elastic and compliance tensors 678 679 680 ! Writes the elastic constants tensors (clamped-ion, relaxed-ion with and 681 ! without stress corrections) to a netCDF file. 682 ! 683 ! compl=relaxed-ion compliance tensor(without stress correction) (6*6) in Voigt notation 684 ! compl_clamped=clamped-ion compliance tensor(without stress correction) (6*6) in Voigt notation 685 ! compl_stress=relaxed-ion compliance tensor(with stress correction) (6*6) in Voigt notation 686 ! elast=relaxed-ion elastic tensor(without stress correction) (6*6) in Voigt notation 687 ! elast_clamped=clamped-ion elastic tensor(without stress correction) (6*6) in Voigt notation 688 ! elast_stress=relaxed-ion elastic tensor(with stress correction) (6*6) in Voigt notation 689 ! Units are GPa for elastic constants and GPa^-1 for compliance constants 690 691 if (ncid /= nctk_noid) then 692 #ifdef HAVE_NETCDF 693 ! Define dimensions 694 NCF_CHECK(nctk_def_basedims(ncid, defmode=.True.)) 695 696 ncerr = nctk_def_iscalars(ncid, [character(len=nctk_slen) :: "asr", "elaflag", "instrflag"]) 697 NCF_CHECK(ncerr) 698 699 !arrays 700 ncerr = nctk_def_arrays(ncid, [& 701 nctkarr_t('internal_strain_tensor', "dp", 'natom3, six'), & 702 nctkarr_t('compliance_constants_relaxed_ion', "dp", 'six, six'), & 703 nctkarr_t('compliance_constants_clamped_ion', "dp", 'six, six'), & 704 nctkarr_t('compliance_constants_relaxed_ion_stress_corrected', "dp", "six, six"), & 705 nctkarr_t('elastic_constants_relaxed_ion', "dp", 'six, six'), & 706 nctkarr_t('elastic_constants_clamped_ion', "dp", 'six, six'), & 707 nctkarr_t('elastic_constants_relaxed_ion_stress_corrected', "dp", 'six, six')]) 708 NCF_CHECK(ncerr) 709 710 ! Write variables. 711 NCF_CHECK(nctk_set_datamode(ncid)) 712 NCF_CHECK(nf90_put_var(ncid, vid('asr'), inp%asr)) 713 NCF_CHECK(nf90_put_var(ncid, vid('elaflag'), inp%elaflag)) 714 NCF_CHECK(nf90_put_var(ncid, vid('instrflag'), inp%instrflag)) 715 NCF_CHECK(nf90_put_var(ncid, vid('internal_strain_tensor'), instrain)) 716 NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion'), compl)) 717 NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_clamped_ion'), compl_clamped)) 718 NCF_CHECK(nf90_put_var(ncid, vid('compliance_constants_relaxed_ion_stress_corrected'), compl_stress)) 719 NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion'), elast)) 720 NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_clamped_ion'), elast_clamped)) 721 NCF_CHECK(nf90_put_var(ncid, vid('elastic_constants_relaxed_ion_stress_corrected'), elast_stress)) 722 723 #else 724 ABI_ERROR("Netcdf support not enabled") 725 #endif 726 end if 727 728 contains 729 730 integer function vid(vname) 731 character(len=*),intent(in) :: vname 732 vid = nctk_idname(ncid, vname) 733 end function vid 734 735 end subroutine ddb_elast
ABINIT/m_ddb_elast [ Modules ]
NAME
m_ddb_elast
FUNCTION
COPYRIGHT
Copyright (C) 1999-2024 ABINIT group (XW, DW) 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
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_ddb_elast 23 24 use defs_basis 25 use m_abicore 26 use m_errors 27 use m_crystal 28 use m_ddb 29 use m_nctk 30 #ifdef HAVE_NETCDF 31 use netcdf 32 #endif 33 34 use m_fstrings, only : itoa, sjoin 35 use m_hide_lapack, only : matrginv 36 use m_dynmat, only : asria_corr 37 use m_anaddb_dataset, only : anaddb_dataset_type 38 39 implicit none 40 41 private