TABLE OF CONTENTS
ABINIT/ddb_hybrid [ Functions ]
NAME
ddb_hybrid
FUNCTION
Modify selected elements in atmfrc, dielt, zeff as specified in an input file named: "modifs.ddb" If this file does not exist, return immediately
COPYRIGHT
Copyright (C) 2000-2018 ABINIT group (PhG,XG) 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
acell = lengths of lattice vectors
asr = flag to impose the acoustic sum rule
atmfrc(2,3,natom,3,natom,nrpt)= Analytical part of the
Interatonmic Force Constants in real space.
dielt(3,3)=dielectric tensor
dipdip = flag to include dipole-dipole contribution
dyew = full Ewald matrix
dyewq0(3,3,natom)=contraction of the Ewald matrix at q=0
modifications to the IFCs.
gmet = reciprocal space metric
gprim = reciprocal lattice vectors
iout = unit number for main output file. If -1 we are in a child mpi process which should not write
natom=number of atom in unit cell
nrpt = number of points in real space for the FT on the dynamical matrices
rcan = canonical positions of atoms
rmet = real-space metric
rprim = unit cell vectors
rpt = positions of points in real space for the FT on the dynamical matrices
ucvol = unit cell volume
wghatm = wieghts for pairs of atoms, in FT interpolations of dyn mat
xred = reduced positions of atoms
zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement.
OUTPUT
atmfrc(2,3,natom,3,natom,nrpt)= modified dielt(3,3)=modified zeff(3,3,natom)=modified
PARENTS
CHILDREN
asrif9,canct9,ewald9,matr3inv,q0dy3_calc
SOURCE
55 #if defined HAVE_CONFIG_H 56 #include "config.h" 57 #endif 58 59 #include "abi_common.h" 60 61 62 subroutine ddb_hybrid(acell,asr,atmfrc,dielt,dipdip,dyew,dyewq0,& 63 & gmet,gprim,iout,natom,nrpt,rcan,rmet,& 64 & rprim,rpt,ucvol,wghatm,xred,zeff) 65 66 use defs_basis 67 use m_profiling_abi 68 use m_errors 69 70 use m_io_tools, only : open_file 71 use m_dynmat, only : q0dy3_calc, asrif9, canct9 72 use m_ewald, only : ewald9 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'ddb_hybrid' 78 use interfaces_32_util 79 !End of the abilint section 80 81 implicit none 82 83 !Arguments ------------------------------------ 84 !scalars 85 integer,intent(in) :: asr,dipdip,iout,natom,nrpt 86 real(dp),intent(in) :: ucvol 87 !arrays 88 real(dp),intent(in) :: acell(3),gmet(3,3),gprim(3,3),rcan(3,natom),rmet(3,3) 89 real(dp),intent(in) :: rprim(3,3),rpt(3,nrpt),wghatm(natom,natom,nrpt) 90 real(dp),intent(in) :: xred(3,natom) 91 real(dp),intent(inout) :: atmfrc(2,3,natom,3,natom,nrpt),dielt(3,3) 92 real(dp),intent(inout) :: dyew(2,3,natom,3,natom),dyewq0(3,3,natom) 93 real(dp),intent(inout) :: zeff(3,3,natom) 94 95 !Local variables------------------------------- 96 !Allocate should be used, instead of these fixed values for msrd and mddd 97 !scalars 98 integer,parameter :: hyun=12,mddd=2,msrd=5000 99 integer :: chk,d1,d2,flgr,genat1,genat2,ia,ib,id1,index,irpt,isrd,jj,kk,mu 100 integer :: nbr,nddd,nsrd,nu,option,sumg0 101 real(dp) :: dd,dd2,detdlt,detdlt2,normr2,rsq,rsq2 102 logical :: ex 103 character(len=500) :: msg 104 !arrays 105 integer :: at1(msrd),at2(msrd),sta(natom,natom,nrpt) 106 real(dp) :: del(3),delta(3),delta2(3),dielt2(3,3),difr(3),dyew2q0(3,3,natom) 107 real(dp) :: ewab(3,3),ewab2(3,3),ewiaf(3,3),ewiaf2(3,3),ifcsr(3,3,msrd) 108 real(dp) :: invdlt(3,3),invdlt2(3,3),qpt(3),ra(3),rb(3) 109 real(dp) :: rpt2(3,msrd),xreda(3),zeff2(3,3,natom) 110 111 ! ****************************************************************** 112 113 #if defined HAVE_OS_MACOSX 114 return ! macOSX seem to have problem with the inquire statement below 115 #endif 116 117 !Do the modifications only if the 'modifs.ddb' file exist. 118 inquire(file='modifs.ddb',exist=ex) 119 if(.not.ex)return 120 121 !PART 1: read input file for modifications 122 !+++++++++++++++++++++++++++++++++++++++++ 123 124 if (iout > 0) then 125 write(iout,*) 126 write(iout,*) 127 write(iout,*)' **WARNING** some interatomic force constants' 128 write(iout,*)' will be artificially modified' 129 end if 130 131 if (open_file('modifs.ddb',msg,unit=hyun,status='old',form='formatted') /= 0) then 132 MSG_ERROR(msg) 133 end if 134 135 read(hyun,*) nsrd 136 if (nsrd<0)then 137 write(std_out,'(a,i8,a)' )& 138 & 'ddb_hybrid : nsrd is',nsrd,', which is lower than 0 => stop' 139 MSG_ERROR("Aborting now") 140 end if 141 if (nsrd>msrd)then 142 write(std_out,'(a,i8,a,a,i8,a)' )& 143 & ' ddb_hybrid : nsrd is',nsrd,', which is larger than',ch10,& 144 & msrd,' , the maximum number allowed => stop' 145 MSG_ERROR("Aborting now") 146 end if 147 if (nsrd/=0)then 148 isrd=0 149 do ia=1,natom 150 read(hyun,*) genat1 151 do ib=1,(nsrd/natom) 152 isrd=isrd+1 153 at1(isrd)=genat1 154 read(hyun,*) genat2 155 at2(isrd)=genat2 156 read(hyun,*) rpt2(1:3,isrd) 157 read(hyun,*) ifcsr(1:3,1,isrd) 158 read(hyun,*) ifcsr(1:3,2,isrd) 159 read(hyun,*) ifcsr(1:3,3,isrd) 160 end do 161 end do 162 end if 163 164 read(hyun,*) nddd 165 if (nddd<0)then 166 write(std_out,'(a,i8,a)' )& 167 & ' ddb_hybrid : nddd is',nddd,', which is lower than 0 => stop' 168 MSG_ERROR("Aborting now") 169 end if 170 if (nddd>mddd)then 171 write(std_out,'(a,i8,a,a,a)' )& 172 & ' ddb_hybrid : nddd is',nddd,', which is not',ch10,& 173 & ' one of the allowed values (0-1-2) => stop' 174 MSG_ERROR("Aborting now") 175 end if 176 if (nddd/=0)then 177 do ia=1,natom 178 do id1=1,3 179 read(hyun,*) zeff2(id1,1:3,ia) 180 end do 181 end do 182 do id1=1,3 183 read(hyun,*) dielt2(id1,1),dielt2(id1,2),dielt2(id1,3) 184 end do 185 end if 186 187 close(hyun) 188 189 !PART 2: include modifications in the SR part 190 !++++++++++++++++++++++++++++++++++++++++++++ 191 192 if (nsrd/=0) then 193 194 ! write(iout,*)' List of SR modifications:' 195 ! write(iout,*)' dir, at, dir, at, R-vec, old-value, new-value' 196 ! write(iout,*)' old-value, new-value, weight' 197 198 chk=0 199 nbr=0 200 201 do ia=1,natom 202 do ib=1,natom 203 do irpt=1,nrpt 204 sta(ia,ib,irpt)=0 205 if (wghatm(ia,ib,irpt)/=0) sta(ia,ib,irpt)=1 206 end do 207 end do 208 end do 209 210 ! BIG loop on the N SR-modifications 211 do isrd=1,nsrd 212 213 ! identify the R-vector 214 flgr=0 215 do irpt=1,nrpt 216 difr(1:3)= abs(rpt2(1:3,isrd)-rpt(1:3,irpt)) 217 if (difr(1)<1.d-4) then 218 if (difr(2)<1.d-4) then 219 if (difr(3)<1.d-4) then 220 flgr=irpt 221 end if 222 end if 223 end if 224 end do 225 226 if (flgr==0)then 227 write(std_out,'(a,a,a,i8,a)' )& 228 & ' ddb_hybrid : not able to identify the',ch10,& 229 & ' R-vector associated to SR-data',isrd,'=> stop' 230 cycle 231 end if 232 233 ! include the modifications in atmfrc 234 if (wghatm(at1(isrd),at2(isrd),flgr)/=0.0_dp) then 235 do d1=1,3 236 do d2=1,3 237 atmfrc(1,d1,at1(isrd),d2,at2(isrd),flgr) & 238 & = ifcsr(d1,d2,isrd)/wghatm(at1(isrd),at2(isrd),flgr) 239 end do 240 end do 241 sta(at1(isrd),at2(isrd),flgr)=0 242 nbr=nbr+1 243 end if 244 245 chk=chk+1 246 247 end do 248 249 ! end if nsrd><0 250 end if 251 252 if (iout > 0) then 253 write(iout,*)' Number of SR modifications:',chk 254 write(iout,*)' Modifications really included:',nbr 255 256 write(iout,*)' Eventual ifc missing:' 257 write(iout,*)' (ia,ib,irpt -- rpt1,rpt2,rpt3)' 258 chk=0 259 do ia=1,natom 260 do ib=1,natom 261 do irpt=1,nrpt 262 if ((wghatm(ia,ib,irpt)/=0).and.(sta(ia,ib,irpt)==1)) then 263 write(iout,*) ia,ib,irpt 264 write(iout,*) rpt(1,irpt),rpt(2,irpt),rpt(2,irpt) 265 chk=1 266 end if 267 end do 268 end do 269 end do 270 if (chk==0) write(iout,*)' -no problem detected-' 271 end if 272 273 274 !PART 3: include modifications in the DD part 275 !++++++++++++++++++++++++++++++++++++++++++++ 276 277 if (nddd/=0) then 278 279 if (iout > 0) then 280 write(iout,*)' The DD interaction has also been modified:' 281 write(iout,*)' The Born effective chages are now:' 282 do ia=1,natom 283 write(iout,'(a,i4)' )' atom',ia 284 do id1=1,3 285 write(iout,*)zeff2(id1,1:3,ia) 286 end do 287 end do 288 write(iout,*)' The dielectric tensor is now:' 289 do id1=1,3 290 write(iout,*) dielt2(id1,1:3) 291 end do 292 end if 293 294 ! The former values are replaced by the new ones in part 4 295 296 ! Modify dyew2q0 accordingly to zeff2 and dielt2 (if dip-dip non-zero) 297 if (dipdip==1) then 298 sumg0=0 299 qpt(1)=0.0_dp 300 qpt(2)=0.0_dp 301 qpt(3)=0.0_dp 302 call ewald9(acell,dielt2,dyew,gmet,gprim,natom,& 303 & qpt,rmet,rprim,sumg0,ucvol,xred,zeff2) 304 if (asr==1.or.asr==2) then 305 option=asr 306 call q0dy3_calc(natom,dyew2q0,dyew,option) 307 else if (asr==0) then 308 dyew2q0(:,:,:)=0.0_dp 309 end if 310 end if 311 312 ! end if nddd/=0 313 end if 314 315 !Eventually keep the total ifc within the box unchanged 316 !This basically corresponds to modify the SR part accordingly 317 !to the change of the DD part within the box... 318 319 if (nddd==2) then 320 321 ! Store the interatomic distances 322 ! call dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt) 323 324 ! calculating the inverse (transpose) of the dielectric tensor 325 call matr3inv(dielt,invdlt) 326 call matr3inv(dielt2,invdlt2) 327 ! calculating the determinant of the dielectric tensor 328 detdlt=dielt(1,1)*dielt(2,2)*dielt(3,3)+dielt(1,3)*dielt(2,1)*& 329 & dielt(3,2)+dielt(1,2)*dielt(2,3)*dielt(3,1)-dielt(1,3)*& 330 & dielt(2,2)*dielt(3,1)-dielt(1,1)*dielt(2,3)*dielt(3,2)-& 331 & dielt(1,2)*dielt(2,1)*dielt(3,3) 332 detdlt2=dielt2(1,1)*dielt2(2,2)*dielt2(3,3)+dielt2(1,3)*& 333 & dielt2(2,1)*& 334 & dielt2(3,2)+dielt2(1,2)*dielt2(2,3)*dielt2(3,1)-& 335 & dielt2(1,3)*& 336 & dielt2(2,2)*dielt2(3,1)-dielt2(1,1)*dielt2(2,3)*& 337 & dielt2(3,2)-& 338 & dielt2(1,2)*dielt2(2,1)*dielt2(3,3) 339 340 ! Big loop on first atom ia 341 do ia=1,natom 342 343 ! First transform canonical coordinates to reduced coordinates 344 xreda(:)=gprim(1,:)*rcan(1,ia)+gprim(2,:)*rcan(2,ia)& 345 & +gprim(3,:)*rcan(3,ia) 346 347 ! Then to cartesian coordinates 348 ra(:)=xreda(1)*acell(1)*rprim(:,1)+& 349 & xreda(2)*acell(2)*rprim(:,2)+& 350 & xreda(3)*acell(3)*rprim(:,3) 351 352 ! Big intra-loop on the atoms in the box ib 353 do index=1,(natom*nrpt) 354 355 call canct9(acell,gprim,ib,index,irpt,natom,nrpt,& 356 & rcan,rb,rprim,rpt) 357 358 if (wghatm(ia,ib,irpt)/=0)then 359 360 del(:)=ra(:)-rb(:) 361 rsq=0.0_dp 362 rsq2=0.0_dp 363 delta(:)=0.0_dp 364 delta2(:)=0.0_dp 365 do jj=1,3 366 do kk=1,3 367 ewab(jj,kk)=0.0_dp 368 ewab2(jj,kk)=0 369 rsq=rsq+del(jj)*invdlt(kk,jj)*del(kk) 370 rsq2=rsq2+del(jj)*invdlt2(kk,jj)*del(kk) 371 delta(kk)=delta(kk)+invdlt(kk,jj)*del(jj) 372 delta2(kk)=delta2(kk)+invdlt2(kk,jj)*del(jj) 373 end do 374 end do 375 dd=sqrt(rsq) 376 dd2=sqrt(rsq2) 377 378 ! Avoid zero denominators in 'term': 379 if (sqrt(rsq)>=1.0d-12) then 380 do mu=1,3 381 do nu=1,3 382 ewab(mu,nu)=(-3*delta(nu)*delta(mu)+invdlt(nu,mu)*dd**2)& 383 & /dd**5/sqrt(detdlt) 384 end do 385 end do 386 else 387 if (ia/=ib)then 388 write(std_out,*)' ddb_hybrid : interatomic distance vanishes ',& 389 & ' Check the input for the following atoms :' 390 write(std_out,*)ia,ib 391 MSG_ERROR("Aborting now") 392 end if 393 end if 394 395 ! Avoid zero denominators in 'term': 396 if (sqrt(rsq)>=1.0d-12) then 397 do mu=1,3 398 do nu=1,3 399 ewab2(mu,nu)=& 400 & (-3*delta2(nu)*delta2(mu)+invdlt2(nu,mu)*dd2**2)& 401 & /dd2**5/sqrt(detdlt2) 402 end do 403 end do 404 else 405 if (ia/=ib)then 406 write(std_out,*)' ddb_hybrid : inter-atomic distance vanishes ',& 407 & ' Check the input for the following atoms :' 408 write(std_out,*)ia,ib 409 MSG_ERROR("Aborting now") 410 end if 411 end if 412 413 ! Take into account the effective charge tensor 414 normr2=rpt(1,irpt)**2+rpt(2,irpt)**2+rpt(3,irpt)**2 415 do mu=1,3 416 do nu=1,3 417 ewiaf(mu,nu)=0.0_dp 418 ewiaf2(mu,nu)=0.0_dp 419 if((ia==ib).and.(normr2<1.d-7))then 420 ewiaf(mu,nu)=-dyewq0(mu,nu,ia) 421 ewiaf2(mu,nu)=-dyew2q0(mu,nu,ia) 422 end if 423 do jj=1,3 424 do kk=1,3 425 ewiaf(mu,nu)=ewiaf(mu,nu)& 426 & +zeff(jj,mu,ia)*zeff(kk,nu,ib)*& 427 & ewab(jj,kk) 428 ewiaf2(mu,nu)=ewiaf2(mu,nu)& 429 & +zeff2(jj,mu,ia)*zeff2(kk,nu,ib)*& 430 & ewab2(jj,kk) 431 end do 432 end do 433 end do 434 end do 435 436 ! add DD(old)-DD(new) to the SR part... 437 do mu=1,3 438 do nu=1,3 439 atmfrc(1,mu,ia,nu,ib,irpt)=atmfrc(1,mu,ia,nu,ib,irpt)& 440 & +(ewiaf(mu,nu)-ewiaf2(mu,nu))/wghatm(ia,ib,irpt) 441 end do 442 end do 443 444 ! end if wghatm><0 445 end if 446 447 ! End loop ia 448 end do 449 450 ! End loop index 451 end do 452 453 ! end if nddd=2 454 end if 455 456 !PART 4: actualize the matrices before exiting 457 !+++++++++++++++++++++++++++++++++++++++++++++ 458 459 if (nddd>0) then 460 461 ! Update of dielt and zeff 462 zeff(:,:,:)=zeff2(:,:,:) 463 dielt(:,:)=dielt2(:,:) 464 465 ! Update of dyewq0 466 dyewq0(:,:,:)=dyew2q0(:,:,:) 467 468 end if 469 470 if ((nsrd>0).or.(nddd==2)) then 471 ! Reimpose the ASR on the ifc that have been modified... 472 if(asr>0)then 473 write(std_out,*)' ddb_hybrid : enter asrif9 ' 474 call asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm) 475 write(std_out,*)' ddb_hybrid : exit asrif9 ' 476 end if 477 end if 478 479 end subroutine ddb_hybrid