TABLE OF CONTENTS
ABINIT/ddb_internalstr [ Functions ]
NAME
ddb_internalstr
FUNCTION
Get the insternal strain tensors,both force response and displacement response ones.
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
blkval(2,3,mpert,3,mpert,nblok)= second derivatives of total energy with respect to electric fields atom displacements,strain,...... all in cartesian coordinates crystal<crystal_t>=Crystalline structure info. asrq0<asrq0_t>=Object for the treatment of the ASR based on the q=0 block found in the DDB file. iblok= bolk number in DDB file iout=out file number mpert=maximum number of ipert msize=Maximum size of dynamical matrices and other perturbations (ddk, dde...) natom=number of atoms in unit cell nblok=number of total bloks in DDB file
OUTPUT
instrain=force response internal strain tensor
NOTES
In output of internal strain tensor,column runs from strain1 to strain6(in Voigt notation),row runs from atom1x,atom1y,atom1z,atom2x,....... sum rule is applied on the internal strain tensor
PARENTS
anaddb,m_effective_potential_file
CHILDREN
asria_corr,wrtout,zhpev
SOURCE
47 #if defined HAVE_CONFIG_H 48 #include "config.h" 49 #endif 50 51 #include "abi_common.h" 52 53 54 subroutine ddb_internalstr(asr,& 55 !&crystal,& 56 & blkval,& 57 !&asrq0,& 58 & d2asr,iblok,instrain,iout,mpert,& 59 !&msize,& 60 natom,nblok) 61 62 use defs_basis 63 use m_profiling_abi 64 use m_errors 65 use m_crystal 66 use m_ddb 67 68 use m_fstrings, only : itoa, sjoin 69 use m_dynmat, only : asria_corr 70 71 !This section has been created automatically by the script Abilint (TD). 72 !Do not modify the following lines by hand. 73 #undef ABI_FUNC 74 #define ABI_FUNC 'ddb_internalstr' 75 use interfaces_14_hidewrite 76 !End of the abilint section 77 78 implicit none 79 80 !Arguments---------------------------------------------- 81 !scalars 82 integer,intent(in) :: asr,iblok,iout,mpert,natom,nblok 83 !integer,intent(in) :: msize 84 !type(crystal_t),intent(in) :: crystal 85 !type(asrq0_t),intent(inout) :: asrq0 86 !arrays 87 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 88 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok) 89 real(dp),intent(out) :: instrain(3*natom,6) 90 91 !Local variables------------------------------------ 92 !scalars 93 integer :: idirA,idirB,ier,ii1,ii2,ipertA,ipertB,ivarA,ivarB 94 character(len=500) :: direction,message 95 !arrays 96 real(dp) :: Amatr(3*natom-3,3*natom-3),Apmatr(3*natom,3*natom) 97 real(dp) :: Bmatr(2,((3*natom-3)*(3*natom-2))/2) 98 real(dp) :: Bpmatr(2,(3*natom*(3*natom+1))/2),Cmatr(3*natom-3,3*natom-3) 99 real(dp) :: Cpmatr(3*natom,3*natom),Nmatr(3*natom,3*natom),deviation(3,6) 100 real(dp) :: eigval(3*natom-3),eigvalp(3*natom),eigvec(2,3*natom-3,3*natom-3) 101 real(dp) :: eigvecp(2,3*natom,3*natom),instrain_dis(6,3*natom) 102 real(dp) :: kmatrix(3*natom,3*natom),zhpev1(2,2*3*natom-4) 103 real(dp) :: zhpev1p(2,2*3*natom-1),zhpev2(3*3*natom-5),zhpev2p(3*3*natom-2) 104 real(dp) :: d2cart(2,3*natom,3*natom) 105 106 !*************************************************************** 107 108 !extract internal strain from DDB matrix data 109 110 do ipertA=1,natom 111 do idirA=1,3 112 ivarA=idirA+3*(ipertA-1) 113 do ivarB=1,6 114 if(ivarB<=3) then 115 idirB=ivarB 116 ipertB=natom+3 117 ! for the diagonal modulus 118 else if(ivarB>3) then 119 idirB=ivarB-3 120 ipertB=natom+4 121 ! for the shear modulus 122 end if 123 instrain(ivarA,ivarB)=(-1.0_dp)*blkval(1,idirA,ipertA,idirB,ipertB,iblok) 124 ! write(std_out,'(es16.6)')blkval(1,idirA,ipertA,idirB,ipertB,iblok) 125 end do 126 end do 127 end do 128 !according to the definition there is a minus sign before the second derivative 129 130 !apply sum rule to the internal strain tensor 131 deviation(:,:)=zero 132 do ivarB=1,6 133 do ivarA=1,3*natom 134 if(mod(ivarA,3)==0)then 135 deviation(1,ivarB)=deviation(1,ivarB)+instrain(ivarA,ivarB) 136 end if 137 if(mod(ivarA,3)==1)then 138 deviation(2,ivarB)=deviation(2,ivarB)+instrain(ivarA,ivarB) 139 end if 140 if(mod(ivarA,3)==2)then 141 deviation(3,ivarB)=deviation(3,ivarB)+instrain(ivarA,ivarB) 142 end if 143 end do 144 end do 145 146 do ivarB=1,6 147 do ivarA=1,3*natom 148 if(mod(ivarA,3)==0)then 149 instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(1,ivarB)/natom 150 end if 151 if(mod(ivarA,3)==1)then 152 instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(2,ivarB)/natom 153 end if 154 if(mod(ivarA,3)==2)then 155 instrain(ivarA,ivarB)=instrain(ivarA,ivarB)-deviation(3,ivarB)/natom 156 end if 157 end do 158 end do 159 !ending the sum rule 160 161 !print the force response internal strain constants into the output file 162 163 write(message,'(a,a,a,a)')ch10,& 164 & ' Force-response internal strain tensor','(Unit:Hartree/bohr)',ch10 165 call wrtout(std_out,message,'COLL') 166 call wrtout(iout,message,'COLL') 167 168 write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',& 169 & 'strainyy','strainzz','strainyz','strainxz','strainxy' 170 call wrtout(std_out,message,'COLL') 171 do ii1=1,3*natom 172 if(mod(ii1,3)==1)then 173 direction='x' 174 elseif(mod(ii1,3)==2)then 175 direction='y' 176 elseif(mod(ii1,3)==0)then 177 direction='z' 178 end if 179 write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,& 180 & instrain(ii1,1),instrain(ii1,2),instrain(ii1,3),& 181 & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6) 182 call wrtout(std_out,message,'COLL') 183 end do 184 185 !now write into the ddb output file 186 write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',& 187 & 'strainyy','strainzz','strainyz','strainxz','strainxy' 188 call wrtout(iout,message,'COLL') 189 do ii1=1,3*natom 190 if(mod(ii1,3)==1)then 191 direction='x' 192 elseif(mod(ii1,3)==2)then 193 direction='y' 194 elseif(mod(ii1,3)==0)then 195 direction='z' 196 end if 197 write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ii1-1)/3)+1,' ',direction,& 198 & instrain(ii1,1),& 199 & instrain(ii1,2),instrain(ii1,3),& 200 & instrain(ii1,4),instrain(ii1,5),instrain(ii1,6) 201 call wrtout(iout,message,'COLL') 202 end do 203 204 !try to get the displacement response internal strain tensor 205 !first need the inverse of force constant matrix 206 d2cart = zero 207 do ipertA=1,natom 208 do ii1=1,3 209 ivarA=ii1+3*(ipertA-1) 210 do ipertB=1,natom 211 do ii2=1,3 212 ivarB=ii2+3*(ipertB-1) 213 d2cart(1,ivarA,ivarB)=blkval(1,ii1,ipertA,ii2,ipertB,iblok) 214 end do 215 end do 216 end do 217 end do 218 219 !Eventually impose the acoustic sum rule 220 !FIXME: this might depend on ifcflag: impose that it is 0 or generalize 221 call asria_corr(asr,d2asr,d2cart,natom,natom) 222 !call asrq0_apply(asrq0, natom, mpert, msize, crystal%xcart, d2cart) 223 kmatrix = d2cart(1,:,:) 224 Apmatr(:,:)=kmatrix(:,:) 225 226 !DEBUG 227 !write(std_out,'(/,a,/)')'the force constant matrix' 228 !do ivarA=1,3*natom 229 !write(std_out,'(/)') 230 !do ivarB=1,3*natom 231 !write(std_out,'(es16.6)')kmatrix(ivarB,ivarA) 232 !end do 233 !end do 234 !ENDDEBUG 235 236 Nmatr(:,:)=0.0_dp 237 do ivarA=1,3*natom 238 do ivarB=1,3*natom 239 if (mod(ivarA,3)==0 .and. mod(ivarB,3)==0)then 240 Nmatr(ivarA,ivarB)=one 241 end if 242 if (mod(ivarA,3)==1 .and. mod(ivarB,3)==1)then 243 Nmatr(ivarA,ivarB)=one 244 end if 245 if (mod(ivarA,3)==2 .and. mod(ivarB,3)==2)then 246 Nmatr(ivarA,ivarB)=one 247 end if 248 end do 249 end do 250 251 !DEBUG 252 !do ivarA=1,3*natom 253 !write(std_out,'(/)') 254 !do ivarB=1,3*natom 255 !write(std_out,'(es16.6)')Nmatr(ivarB,ivarA) 256 !end do 257 !end do 258 !ENDDEBUG 259 260 !starting the pseudoinervering processes 261 !then get the eigenvectors of the big matrix,give values to matrixBp 262 Bpmatr=0.0_dp 263 ii1=1 264 do ivarA=1,3*natom 265 do ivarB=1,ivarA 266 Bpmatr(1,ii1)=Nmatr(ivarB,ivarA) 267 ii1=ii1+1 268 end do 269 end do 270 271 !Bpmatr(2,:) is the imaginary part of the force matrix 272 !then call the subroutines CHPEV and ZHPEV to get the eigenvectors 273 call ZHPEV ('V','U',3*natom,Bpmatr,eigvalp,eigvecp,3*natom,zhpev1p,zhpev2p,ier) 274 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 275 276 !DEBUG 277 !the eigenval and eigenvec 278 !write(std_out,'(/,a,/)')'the eigenvalues and eigenvectors' 279 !do ivarA=1,3*natom 280 !write(std_out,'(/)') 281 !write(std_out,'(es16.6)')eigvalp(ivarA) 282 !end do 283 !do ivarA=1,3*natom 284 !write(std_out,'(/)') 285 !do ivarB=1,3*natom 286 !write(std_out,'(es16.6)')eigvecp(1,ivarB,ivarA) 287 !end do 288 !end do 289 !ENDDEBUG 290 291 !Then do the multiplication to get the reduced matrix,in two steps 292 !After this the force constant matrix is decouple in two bloks, 293 !acoustic and optical ones 294 Cpmatr(:,:)=0.0_dp 295 do ivarA=1,3*natom 296 do ivarB=1,3*natom 297 do ii1=1,3*natom 298 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+eigvecp(1,ii1,ivarA)*Apmatr(ii1,ivarB) 299 end do 300 end do 301 end do 302 303 Apmatr(:,:)=0.0_dp 304 do ivarA=1,3*natom 305 do ivarB=1,3*natom 306 do ii1=1,3*natom 307 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+Cpmatr(ivarA,ii1)*eigvecp(1,ii1,ivarB) 308 end do 309 end do 310 end do 311 312 !DEBUG 313 !the blok diago 314 !write(std_out,'(/,a,/)')'matrixAp' 315 !do ivarA=1,3*natom 316 !write(std_out,'(/)') 317 !do ivarB=1,3*natom 318 !write(std_out,'(es16.6)')Apmatr(ivarA,ivarB) 319 !end do 320 !end do 321 !ENDDEBUG 322 323 !Check the last three eigenvalues whether too large or not 324 ivarB=0 325 do ivarA=3*natom-2,3*natom 326 if (ABS(Apmatr(ivarA,ivarA))>tol6)then 327 ivarB=1 328 end if 329 end do 330 331 if(ivarB==1)then 332 write(message,'(a,a,a,a,a,a,a,a,3es16.6)')ch10,& 333 & ' Acoustic sum rule violation met : the eigenvalues of accoustic mode',ch10,& 334 & ' are too large at Gamma point.',ch10,& 335 & ' Increase cutoff energy or k-points sampling.',ch10,& 336 & ' The three eigenvalues are:',Apmatr(3*natom-2,3*natom-2),Apmatr(3*natom-1,natom-1),Apmatr(3*natom,3*natom) 337 MSG_WARNING(message) 338 call wrtout(iout,message,'COLL') 339 end if 340 341 !Give the value of reduced matrix form Apmatr to Amatr 342 do ivarA=1,3*natom-3 343 do ivarB=1,3*natom-3 344 Amatr(ivarA,ivarB)=Apmatr(ivarA,ivarB) 345 end do 346 end do 347 348 !Now the reduced matrix is in the matrixA, the convert it 349 !first give the give the value of matixB from matrixA 350 ii1=1 351 do ivarA=1,3*natom-3 352 do ivarB=1,ivarA 353 Bmatr(1,ii1)=Amatr(ivarB,ivarA) 354 ii1=ii1+1 355 end do 356 end do 357 Bmatr(2,:)=0.0_dp 358 359 !Call the subroutines CHPEV and ZHPEV to get the eigenvectors and the eigenvalues 360 call ZHPEV ('V','U',3*natom-3,Bmatr,eigval,eigvec,3*natom-3,zhpev1,zhpev2,ier) 361 ABI_CHECK(ier == 0, sjoin("ZHPEV returned:", itoa(ier))) 362 363 !Check the unstable phonon modes, if the first is negative then print 364 !warning message 365 if(eigval(1)<-1.0*tol8)then 366 write(message,'(9a)') ch10,& 367 & ' --- !WARNING',ch10,& 368 & ' Unstable eigenvalue detected in force constant matrix at Gamma point',ch10,& 369 & ' The system under calculation is physically unstable.',ch10,& 370 & ' ---',ch10 371 call wrtout(std_out,message,'COLL') 372 end if 373 374 !Do the matrix muplication to get pseudoinverse inverse matrix 375 Cmatr(:,:)=0.0_dp 376 Amatr(:,:)=0.0_dp 377 do ivarA=1,3*natom-3 378 Cmatr(ivarA,ivarA)=1.0_dp/eigval(ivarA) 379 end do 380 381 do ivarA=1,3*natom-3 382 do ivarB=1,3*natom-3 383 do ii1=1,3*natom-3 384 Amatr(ivarA,ivarB)=Amatr(ivarA,ivarB)+eigvec(1,ivarA,ii1)*Cmatr(ii1,ivarB) 385 end do 386 end do 387 end do 388 389 390 !The second mulplication 391 Cmatr(:,:)=0.0_dp 392 do ivarA=1,3*natom-3 393 do ivarB=1,3*natom-3 394 do ii1=1,3*natom-3 395 Cmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB)+ Amatr(ivarA,ii1)*eigvec(1,ivarB,ii1) 396 end do 397 end do 398 end do 399 400 !DEBUG 401 !write(std_out,'(/,a,/)')'the pseudo inverse of the force matrix' 402 !do ivarA=1,3*natom 403 !write(std_out,'(/)') 404 !do ivarB=1,3*natom 405 !write(std_out,'(es16.6)')Cmatr(ivarA,ivarB) 406 !end do 407 !end do 408 !ENDDEBUG 409 410 !So now the inverse of the reduced matrix is in the matrixC 411 !now do another mulplication to get the pseudoinverse of the original 412 Cpmatr(:,:)=0.0_dp 413 Apmatr(:,:)=0.0_dp 414 do ivarA=1,3*natom-3 415 do ivarB=1,3*natom-3 416 Cpmatr(ivarA,ivarB)=Cmatr(ivarA,ivarB) 417 end do 418 end do 419 420 !Now times the eigvecp 421 do ivarA=1,3*natom 422 do ivarB=1,3*natom 423 do ii1=1,3*natom 424 Apmatr(ivarA,ivarB)=Apmatr(ivarA,ivarB)+eigvecp(1,ivarA,ii1)*& 425 & Cpmatr(ii1,ivarB) 426 end do 427 end do 428 end do 429 Cpmatr(:,:)=0.0_dp 430 do ivarA=1,3*natom 431 do ivarB=1,3*natom 432 do ii1=1,3*natom 433 Cpmatr(ivarA,ivarB)=Cpmatr(ivarA,ivarB)+ Apmatr(ivarA,ii1)*eigvecp(1,ivarB,ii1) 434 end do 435 end do 436 end do 437 438 !Now the inverse is in Cpmatr 439 kmatrix(:,:)=Cpmatr(:,:) 440 !transfer the inverse of k-matrix back to the k matrix 441 !so now the inverse of k matrix is in the kmatrix 442 !ending the part for pseudoinversing the K matrix 443 444 !Now do simple mulplication to obtain the displacement response 445 !internal strain tensor 446 447 instrain_dis(:,:)=0.0_dp 448 do ivarA=1,6 449 do ivarB=1,3*natom 450 do ii1=1,3*natom 451 instrain_dis(ivarA,ivarB)=instrain_dis(ivarA,ivarB)+& 452 & instrain(ii1,ivarA)*kmatrix(ii1,ivarB) 453 end do 454 end do 455 end do 456 457 !Print out the results 458 write(message,'(a,a,a,a)')ch10,& 459 & ' Displacement-response internal strain ', 'tensor (Unit:Bohr)',ch10 460 call wrtout(std_out,message,'COLL') 461 call wrtout(iout,message,'COLL') 462 write(message,'(a5,a4,a11,a12,a12,a12,a12,a12)')' Atom',' dir','strainxx',& 463 & 'strainyy','strainzz','strainyz','strainxz','strainxy' 464 call wrtout(std_out,message,'COLL') 465 call wrtout(iout,message,'COLL') 466 do ivarA=1,3*natom 467 if(mod(ivarA,3)==1)then 468 direction='x' 469 elseif(mod(ivarA,3)==2)then 470 direction='y' 471 elseif(mod(ivarA,3)==0)then 472 direction='z' 473 end if 474 write(message,'(a1,i2,a2,a3,6f12.7)')' ',int((ivarA-1)/3)+1,' ',direction,& 475 & instrain_dis(1,ivarA),instrain_dis(2,ivarA),& 476 & instrain_dis(3,ivarA),instrain_dis(4,ivarA),instrain_dis(5,ivarA),& 477 & instrain_dis(6,ivarA) 478 call wrtout(std_out,message,'COLL') 479 call wrtout(iout,message,'COLL') 480 end do 481 482 end subroutine ddb_internalstr