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.
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 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
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.
PARENTS
anaddb
CHILDREN
asria_corr,matrginv,wrtout,zhpev
SOURCE
53 #if defined HAVE_CONFIG_H 54 #include "config.h" 55 #endif 56 57 #include "abi_common.h" 58 59 60 subroutine ddb_elast(anaddb_dtset,crystal,blkval,compl,compl_clamped,compl_stress,d2asr,& 61 & elast,elast_clamped,elast_stress,iblok,iblok_stress,& 62 & instrain,iout,mpert,& 63 ! &msize,& 64 & natom,nblok) 65 66 use defs_basis 67 use m_profiling_abi 68 use m_errors 69 use m_crystal 70 use m_ddb 71 72 use m_fstrings, only : itoa, sjoin 73 use m_abilasi, only : matrginv 74 use m_dynmat, only : asria_corr 75 use m_anaddb_dataset, only : anaddb_dataset_type 76 77 !This section has been created automatically by the script Abilint (TD). 78 !Do not modify the following lines by hand. 79 #undef ABI_FUNC 80 #define ABI_FUNC 'ddb_elast' 81 use interfaces_14_hidewrite 82 !End of the abilint section 83 84 implicit none 85 86 !Arguments ------------------------------------------- 87 !scalars 88 integer,intent(in) :: iblok,iblok_stress,iout,mpert,natom,nblok 89 !integer,intent(in) :: msize 90 type(crystal_t),intent(in) :: crystal 91 type(anaddb_dataset_type),intent(in) :: anaddb_dtset 92 !arrays 93 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),instrain(3*natom,6) 94 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 95 real(dp),intent(out) :: compl(6,6), compl_clamped(6,6),compl_stress(6,6) 96 real(dp),intent(out) :: elast(6,6), elast_clamped(6,6),elast_stress(6,6) 97 98 !Local variables------------------------------------ 99 !scalars 100 integer :: ier,ii1,ii2,ipert1,ipert2,ivarA,ivarB 101 real(dp) :: ucvol 102 logical :: iwrite 103 character(len=500) :: message 104 !arrays 105 real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom) 106 real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2) 107 real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3) 108 real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom) 109 real(dp) :: compl_relaxed(6,6),eigval(3*natom-3) 110 real(dp) :: eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3) 111 real(dp) :: eigvecp(2,3*natom,3*natom),elast_relaxed(6,6) 112 real(dp) :: kmatrix(3*natom,3*natom),new1(6,3*natom) 113 real(dp) :: new2(6,6),stress(6),zhpev1(2,2*3*natom-4),zhpev1p(2,2*3*natom-1) 114 real(dp) :: zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2) 115 real(dp) :: d2cart(2,3*natom,3*natom) 116 117 !*************************************************************************** 118 119 ucvol = crystal%ucvol 120 iwrite = iout > 0 121 122 !extraction of the elastic constants from the blkvals 123 124 do ivarA=1,6 125 do ivarB=1,6 126 ! because the elastic constant is 6*6, 127 ! so we should judge if the idir is larger than 3 or not 128 if(ivarA>3) then 129 ii1=ivarA-3 130 ipert1=natom+4 !for the shear modulus 131 else if(ivarA<=3) then 132 ii1=ivarA 133 ipert1=natom+3 !for the diagonal part 134 end if 135 if(ivarB>3) then 136 ii2=ivarB-3 137 ipert2=natom+4 !for the shear modulus 138 else if(ivarB<=3) then 139 ii2=ivarB 140 ipert2=natom+3 !for the diagonal part 141 end if 142 elast(ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok) 143 end do 144 end do 145 146 !then consider the volume, because the unit above is in 147 !Hartree, in fact the elastic constant should be in 148 !the units of presure, the energy/volume 149 !And then transform the unit to si unit using GPa 150 !from Hartree/Bohr^3 151 152 do ivarA=1,6 153 do ivarB=1,6 154 elast(ivarA,ivarB)=(elast(ivarA,ivarB)/ucvol)*HaBohr3_GPa 155 end do 156 end do 157 158 !then should consider the two situations:clamped and relaxed 159 !ions respectively,give the initial value of elast_clamped 160 elast_clamped(:,:)=elast(:,:) 161 elast_relaxed(:,:)=elast(:,:) 162 163 !then do the matrix mulplication of instrain*K*instrain to get the 164 !correction of the relaxed ion quantities, in case natom/=1 165 166 if( (anaddb_dtset%elaflag==2 .or. anaddb_dtset%elaflag==3& 167 & .or. anaddb_dtset%elaflag==4 .or. anaddb_dtset%elaflag==5) .and. natom/=1 )then 168 ! extracting force matrix at gamma 169 d2cart = zero 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 d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipert1,ii2,ipert2,iblok) 177 end do 178 end do 179 end do 180 end do 181 182 ! Eventually impose the acoustic sum rule 183 ! FIXME: this might depend on ifcflag: impose that it is 0 or generalize 184 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 185 call asria_corr(anaddb_dtset%asr,d2asr,d2cart,natom,natom) 186 kmatrix = d2cart(1,:,:) 187 188 ! DEBUG 189 ! write(std_out,'(/,a,/)')'the k matrix before inverse' 190 ! do ii1=1,3*natom 191 ! write(std_out,'(6es16.3)')kmatrix(ii1,1),kmatrix(ii1,2),kmatrix(ii1,3),& 192 ! & kmatrix(ii1,4),kmatrix(ii1,5),kmatrix(ii1,6) 193 ! end do 194 ! ENDDEBUG 195 196 ! according to formula, invert the kmatrix(3natom,3natom) 197 Apmatr(:,:)=kmatrix(:,:) 198 199 ! NOTE: MJV 13/3/2011 This is just the 3x3 unit matrix copied throughout the dynamical matrix 200 Nmatr(:,:)=zero 201 do ivarA=1,3*natom 202 do ivarB=1,3*natom 203 if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0) Nmatr(ivarA,ivarB)=one 204 if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1) Nmatr(ivarA,ivarB)=one 205 if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2) Nmatr(ivarA,ivarB)=one 206 end do 207 end do 208 209 ! DEBUG 210 ! The k matrix is not the inverse here - it has not been changed! 211 ! write(std_out,'(/,a,/)')'the direct inverse of the Kmatrix' 212 ! do ivarA=1,3*natom 213 ! write(std_out,'(/)') 214 ! do ivarB=1,3*natom 215 ! write(std_out,'(es16.6)')kmatrix(ivarA,ivarB) 216 ! end do 217 ! end do 218 ! ENDDEBUG 219 220 221 ! starting the pseudo-inverse processes 222 ! then get the eigenvectors of the big matrix,give values to matrixBp 223 ! Pack the Nmatr matrix in Hermitian form 224 ii1=1 225 do ivarA=1,3*natom 226 do ivarB=1,ivarA 227 Bpmatr(1,ii1)=Nmatr(ivarB,ivarA) 228 ii1=ii1+1 229 end do 230 end do 231 Bpmatr(2,:)=zero !the imaginary part of the force matrix 232 ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors 233 ! NOTE: MJV there is a huge indeterminacy in this matrix, which has all identical 3x3 block lines 234 ! this means the orientation of the 0-eigenvalue eigenvectors is kind of random... 235 ! Is the usage just to get out the translational modes? We know what the eigenvectors look like already! 236 ! The translational modes are the last 3 with eigenvalue 6 237 ! 238 call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier) 239 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 240 241 ! DEBUG 242 ! the eigenval and eigenvec 243 ! write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors' 244 ! do ivarA=1,3*natom 245 ! write(std_out,'(/)') 246 ! write(std_out,'(es16.6)')eigvalp(ivarA) 247 ! end do 248 ! do ivarA=1,3*natom 249 ! write(std_out,'(/)') 250 ! do ivarB=1,3*natom 251 ! write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA) 252 ! end do 253 ! end do 254 ! ENDDEBUG 255 256 ! do the multiplication to get the reduced matrix,in two steps 257 ! rotate to eigenbasis constructed above to isolate acoustic modes 258 Cpmatr(:,:)=zero 259 do ivarA=1,3*natom 260 do ivarB=1,3*natom 261 do ii1=1,3*natom 262 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*& 263 & Apmatr(ii1,ivarB) 264 end do 265 end do 266 end do 267 268 Apmatr(:,:)=zero 269 do ivarA=1,3*natom 270 do ivarB=1,3*natom 271 do ii1=1,3*natom 272 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*& 273 & eigvecp(1,ii1,ivarB) 274 end do 275 end do 276 end do 277 278 ! DEBUG 279 ! the blok diagonal parts 280 ! write(std_out,'(/,a,/)')'Apmatr' 281 ! do ivarA=1,3*natom 282 ! write(std_out,'(/)') 283 ! do ivarB=1,3*natom 284 ! write(std_out,'(es16.6)')Apmatr(ivarA,ivarB) 285 ! end do 286 ! end do 287 ! ENDDEBUG 288 289 290 ! check the last three eigenvalues whether too large 291 ivarB=0 292 do ivarA=3*natom-2,3*natom 293 if (ABS(Apmatr(ivarA,ivarA))>tol6) ivarB=1 294 end do 295 if(ivarB==1)then 296 write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,& 297 & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,& 298 & ' are too large at Gamma point',ch10,& 299 & ' increase cutoff energy or k-points sampling.',ch10,& 300 & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),& 301 & Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom) 302 MSG_WARNING(message) 303 call wrtout(iout,message,'COLL') 304 end if 305 ! then give the value of reduced matrix form Apmatr to Amatr 306 do ivarA=1,3*natom-3 307 do ivarB=1,3*natom-3 308 Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB) 309 end do 310 end do 311 312 ! now the reduced matrix is in the Amatr, the convert it 313 ! first give the give the value of Bmatr from Amatr 314 ii1=1 315 do ivarA=1,3*natom-3 316 do ivarB=1,ivarA 317 Bmatr(1,ii1)=Amatr(ivarB,ivarA) 318 ii1=ii1+1 319 end do 320 end do 321 Bmatr(2,:)=zero 322 ! then call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues 323 call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier) 324 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 325 326 ! check the unstable phonon modes, if the first is negative then print warning message 327 if(eigval(1)<-1.0*tol8)then 328 write(message,'(a,a,a,a)') ch10,& 329 & 'Unstable eigenvalue detected in force constant matrix at Gamma point.',ch10,& 330 & 'The system under calculation is physically unstable.' 331 MSG_WARNING(message) 332 call wrtout(iout,message,'COLL') 333 end if 334 335 ! the do the matrix muplication to get pseudoinverse inverse matrix 336 Cmatr(:,:)=zero 337 Amatr(:,:)=zero 338 do ivarA=1,3*natom-3 339 Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA) 340 end do 341 do ivarA=1,3*natom-3 342 do ivarB=1,3*natom-3 343 do ii1=1,3*natom-3 344 Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*& 345 & Cmatr(ii1,ivarB) 346 end do 347 end do 348 end do 349 350 ! then the second mulplication 351 Cmatr(:,:)=zero 352 do ivarA=1,3*natom-3 353 do ivarB=1,3*natom-3 354 do ii1=1,3*natom-3 355 Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+& 356 & Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1) 357 end do 358 end do 359 end do 360 361 ! DEBUG 362 ! write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix' 363 ! do ivarA=1,3*natom 364 ! write(std_out,'(/)') 365 ! do ivarB=1,3*natom 366 ! write(std_out,'(es16.6)')Cmatr(ivarA,ivarB) 367 ! end do 368 ! end do 369 ! ENDDEBUG 370 371 ! so now the inverse of the reduced matrix is in the matrixC 372 ! now do another mulplication to get the pseudoinverse of the original 373 Cpmatr(:,:)=zero 374 Apmatr(:,:)=zero 375 do ivarA=1,3*natom-3 376 do ivarB=1,3*natom-3 377 Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB) 378 end do 379 end do 380 381 ! times the eigvecp 382 do ivarA=1,3*natom 383 do ivarB=1,3*natom 384 do ii1=1,3*natom 385 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*& 386 & Cpmatr(ii1,ivarB) 387 end do 388 end do 389 end do 390 Cpmatr(:,:)=zero 391 do ivarA=1,3*natom 392 do ivarB=1,3*natom 393 do ii1=1,3*natom 394 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+& 395 & Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1) 396 end do 397 end do 398 end do 399 400 ! now the inverse in in Cpmatr 401 kmatrix(:,:)=Cpmatr(:,:) 402 ! transfer the inverse of k-matrix back to the k matrix 403 ! so now the inverse of k matrix is in the kmatrix 404 ! ending the part for pseudoinversing the K matrix 405 ! then do the first matrix mulplication 406 new1(:,:)=zero 407 do ii1=1,6 408 do ii2=1,3*natom 409 do ivarA=1,3*natom 410 new1(ii1,ii2)=new1(ii1,ii2)+instrain(ivarA,ii1)*& 411 & kmatrix(ivarA,ii2) 412 end do 413 end do 414 end do 415 ! then do the second matrix mulplication, and change the value of kmatrix 416 new2(:,:)=zero 417 do ii1=1,6 418 do ii2=1,6 419 do ivarB=1,3*natom 420 new2(ii1,ii2)=new2(ii1,ii2)+new1(ii1,ivarB)*& 421 & instrain(ivarB,ii2) 422 end do 423 end do 424 end do 425 ! then finish the matrix mupl., consider the unit cellvolume 426 ! and the unit change next step 427 do ivarA=1,6 428 do ivarB=1,6 429 new2(ivarA,ivarB)=(new2(ivarA,ivarB)/ucvol)*HaBohr3_GPa 430 end do 431 end do 432 ! then the relaxed one should be the previous one minus the new2 element 433 do ivarA=1,6 434 do ivarB=1,6 435 elast_relaxed(ivarA,ivarB)=elast_relaxed(ivarA,ivarB)-& 436 & new2(ivarA,ivarB) 437 end do 438 end do 439 end if 440 !the above end if end if for elaflag=2 or elafalg=3 or elafalg=4, 441 !or elafalg=5 in line 125 442 443 !DEBUG 444 !write(std_out,'(/,a,/)')'debug the unit cell volume' 445 !write(std_out,'(2es16.6)')ucvol,HaBohr3_GPa 446 !ENDDEBUG 447 448 !then give the initial value of the compl_relaxed(6,6) 449 compl_relaxed(:,:)=elast_relaxed(:,:) 450 451 !******************************************************************* 452 if(anaddb_dtset%elaflag==1.or. anaddb_dtset%elaflag==3)then 453 ! print out the clamped-ion elastic constants to output file 454 write(message,'(3a)')ch10,' Elastic Tensor (clamped ion) (unit:10^2GP):',ch10 455 call wrtout(std_out,message,'COLL') 456 do ivarA=1,6 457 write(std_out,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,& 458 & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,& 459 & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp 460 end do 461 462 call wrtout(iout,message,'COLL') 463 if (iwrite) then 464 do ivarA=1,6 465 write(iout,'(6f12.7)')elast(ivarA,1)/100.00_dp,elast(ivarA,2)/100.00_dp,& 466 & elast(ivarA,3)/100.00_dp,elast(ivarA,4)/100.00_dp,& 467 & elast(ivarA,5)/100.00_dp,elast(ivarA,6)/100.00_dp 468 end do 469 end if 470 end if 471 472 if(anaddb_dtset%elaflag==2.or.anaddb_dtset%elaflag==3& 473 & .or. anaddb_dtset%elaflag==4.or. anaddb_dtset%elaflag==5)then 474 if(anaddb_dtset%instrflag==0)then 475 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 476 & 'in order to get the elastic tensor(relaxed ion), ',ch10,& 477 & 'one needs information about internal strain ',ch10,& 478 & 'one should set instrflag==1;',ch10,& 479 & 'otherwise the program will continue but give wrong values.' 480 MSG_WARNING(message) 481 call wrtout(iout,message,'COLL') 482 end if 483 484 write(message,'(5a)')ch10,& 485 & ' Elastic Tensor (relaxed ion) (unit:10^2GP):',ch10,& 486 & ' (at fixed electric field boundary condition)',ch10 487 call wrtout(std_out,message,'COLL') 488 do ivarA=1,6 489 write(std_out,'(6f12.7)')elast_relaxed(ivarA,1)/100.00_dp,& 490 & elast_relaxed(ivarA,2)/100.00_dp,elast_relaxed(ivarA,3)/100.00_dp,& 491 & elast_relaxed(ivarA,4)/100.00_dp,elast_relaxed(ivarA,5)/100.00_dp,& 492 & elast_relaxed(ivarA,6)/100.00_dp 493 end do 494 495 call wrtout(iout,message,'COLL') 496 if (iwrite) then 497 do ivarA=1,6 498 write(iout,'(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 end if 504 end if 505 506 !then print the corresponding compliances 507 508 if(anaddb_dtset%elaflag==1.or.anaddb_dtset%elaflag==3)then 509 ! compl(:,:)=elast_clamped(:,:) !convert the elastic tensor 510 compl_clamped(:,:)=elast_clamped(:,:) 511 call matrginv(compl_clamped,6,6) 512 write(message,'(a,a,a)')ch10,' Compliance Tensor (clamped ion) (unit: 10^-2GP^-1):',ch10 513 call wrtout(std_out,message,'COLL') 514 515 do ivarB=1,6 516 write(std_out,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,& 517 & compl_clamped(ivarB,2)*100.00_dp,& 518 & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,& 519 & compl_clamped(ivarB,5)*100.00_dp,& 520 & compl_clamped(ivarB,6)*100.00_dp 521 end do 522 523 call wrtout(iout,message,'COLL') 524 525 if (iwrite) then 526 do ivarB=1,6 527 write(iout,'(6f12.7)')compl_clamped(ivarB,1)*100.00_dp,& 528 & compl_clamped(ivarB,2)*100.00_dp,& 529 & compl_clamped(ivarB,3)*100.00_dp,compl_clamped(ivarB,4)*100.00_dp,& 530 & compl_clamped(ivarB,5)*100.00_dp,& 531 & compl_clamped(ivarB,6)*100.00_dp 532 end do 533 end if 534 end if 535 536 if(anaddb_dtset%elaflag==2.or.anaddb_dtset%elaflag==3& 537 & .or. anaddb_dtset%elaflag==4 .or. anaddb_dtset%elaflag==5)then 538 ! compl(:,:)=elast_relaxed(:,:) 539 call matrginv(compl_relaxed,6,6) 540 if(anaddb_dtset%instrflag==0)then 541 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 542 & 'in order to get the compliance tensor(relaxed ion), ',ch10,& 543 & 'one needs information about internal strain ',ch10,& 544 & 'one should set instrflag==1;',ch10,& 545 & 'otherwise the program will continue but give wrong values.' 546 MSG_WARNING(message) 547 call wrtout(iout,message,'COLL') 548 end if 549 write(message,'(5a)')ch10,& 550 & ' Compliance Tensor (relaxed ion) (unit: 10^-2GP^-1):',ch10,& 551 & ' (at fixed electric field boundary condition)',ch10 552 call wrtout(std_out,message,'COLL') 553 554 do ivarB=1,6 555 write(std_out,'(6f12.7)')compl_relaxed(ivarB,1)*100.00_dp,& 556 & compl_relaxed(ivarB,2)*100.00_dp,& 557 & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,& 558 & compl_relaxed(ivarB,5)*100.00_dp,& 559 & compl_relaxed(ivarB,6)*100.00_dp 560 end do 561 call wrtout(iout,message,'COLL') 562 563 if (iwrite) then 564 do ivarB=1,6 565 write(iout,'(6f12.7)')compl_relaxed(ivarB,1)*100.00,& 566 & compl_relaxed(ivarB,2)*100.00_dp,& 567 & compl_relaxed(ivarB,3)*100.00_dp,compl_relaxed(ivarB,4)*100.00_dp,& 568 & compl_relaxed(ivarB,5)*100.00_dp,& 569 & compl_relaxed(ivarB,6)*100.00_dp 570 end do 571 end if 572 end if 573 574 !befor the end , make sure the tensor elast(6,6) 575 !will have the relaxed ion values and tensor elast_clamped has the clamped ion 576 !values, and similarily for the corresponding compliance tensors 577 elast_clamped(:,:)=elast(:,:) 578 elast(:,:)=elast_relaxed(:,:) 579 compl(:,:)=compl_relaxed(:,:) 580 581 !begin the part of computing stress corrected elastic tensors 582 if(anaddb_dtset%elaflag==5)then 583 584 ! DEBUG 585 ! check the iblok number of first derivative of energy 586 ! write(std_out,'(/,a,/)')'iblok number at 8:00Pm' 587 ! write(std_out,'(i)')iblok_stress 588 ! write(std_out,'(a,f12.7)')'the total energy', blkval(1,1,1) 589 ! write(std_out,*)'',blkval(1,:,:,:,:,iblok_stress) 590 ! write(std_out,*)'',blkval(1,:,7,1,1,iblok_stress) 591 ! ENDDEBUG 592 593 ! firts give the corect stress values diagonal parts 594 stress(1)=blkval(1,1,natom+3,1,1,iblok_stress) 595 stress(2)=blkval(1,2,natom+3,1,1,iblok_stress) 596 stress(3)=blkval(1,3,natom+3,1,1,iblok_stress) 597 ! the shear parts 598 stress(4)=blkval(1,1,natom+4,1,1,iblok_stress) 599 stress(5)=blkval(1,2,natom+4,1,1,iblok_stress) 600 stress(6)=blkval(1,3,natom+4,1,1,iblok_stress) 601 ! then convert the unit from atomic to the GPa unit 602 do ivarA=1,6 603 stress(ivarA)=stress(ivarA)*HaBohr3_GPa 604 end do 605 ! give the initial values of elast_stress tensor 606 elast_stress(:,:)=elast_relaxed(:,:) 607 ! notice that only the first three rows need to be corrected 608 do ivarA=1,3 609 do ivarB=1,6 610 elast_stress(ivarA,ivarB)=elast_stress(ivarA,ivarB)-stress(ivarB) 611 end do 612 end do 613 ! then compute the values of compliance tensor with stress correction 614 compl_stress(:,:)=elast_stress(:,:) 615 call matrginv(compl_stress,6,6) 616 ! then print out the results of stress corrected elastic and compliance tensors 617 if(anaddb_dtset%instrflag==0)then 618 write(message,'(a,a,a,a,a,a,a,a)' )ch10,& 619 & 'In order to get the elastic tensor (relaxed ion with stress correction), ',ch10,& 620 & 'one needs information about internal strain ',ch10,& 621 & 'one should set instrflag==1;',ch10,& 622 & 'otherwise the program will continue but give wrong values.' 623 MSG_WARNING(message) 624 call wrtout(iout,message,'COLL') 625 end if 626 write(message,'(5a)')ch10,& 627 & ' Elastic Tensor (relaxed ion with stress corrected) (unit:10^2GP)',ch10,& 628 & ' (at fixed electric field boundary condition)',ch10 629 call wrtout(std_out,message,'COLL') 630 call wrtout(iout,message,'COLL') 631 do ivarA=1,6 632 write(std_out,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,& 633 & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,& 634 & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp 635 end do 636 if (iwrite) then 637 do ivarA=1,6 638 write(iout,'(6f12.7)')elast_stress(ivarA,1)/100.00_dp,elast_stress(ivarA,2)/100.00_dp,& 639 & elast_stress(ivarA,3)/100.00_dp,elast_stress(ivarA,4)/100.00_dp,& 640 & elast_stress(ivarA,5)/100.00_dp,elast_stress(ivarA,6)/100.00_dp 641 end do 642 end if 643 644 ! then the complinace tensors with stress correction 645 write(message,'(5a)')ch10,& 646 & ' Compliance Tensor (relaxed ion with stress correction) (unit: 10^-2(GP)^-1):',ch10,& 647 & ' (at fixed electric field boundary condition)',ch10 648 call wrtout(std_out,message,'COLL') 649 call wrtout(iout,message,'COLL') 650 do ivarB=1,6 651 write(std_out,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,& 652 & compl_stress(ivarB,2)*100.00_dp,& 653 & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,& 654 & compl_stress(ivarB,5)*100.00_dp,& 655 & compl_stress(ivarB,6)*100.00_dp 656 end do 657 if (iwrite) then 658 do ivarB=1,6 659 write(iout,'(6f12.7)')compl_stress(ivarB,1)*100.00_dp,& 660 & compl_stress(ivarB,2)*100.00_dp,& 661 & compl_stress(ivarB,3)*100.00_dp,compl_stress(ivarB,4)*100.00_dp,& 662 & compl_stress(ivarB,5)*100.00_dp,& 663 & compl_stress(ivarB,6)*100.00_dp 664 end do 665 end if 666 end if 667 !end the if 510th line 668 !end the part of stress corrected elastic and compliance tensors 669 670 end subroutine ddb_elast