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