TABLE OF CONTENTS
- ABINIT/m_dynmat
- m_dynmat/asria_calc
- m_dynmat/asria_corr
- m_dynmat/asrif9
- m_dynmat/asrprs
- m_dynmat/axial9
- m_dynmat/bigbx9
- m_dynmat/canat9
- m_dynmat/canct9
- m_dynmat/cart29
- m_dynmat/cart39
- m_dynmat/chkph3
- m_dynmat/chkrp9
- m_dynmat/chneu9
- m_dynmat/d2cart_to_red
- m_dynmat/d2sym3
- m_dynmat/d3sym
- m_dynmat/dfpt_phfrq
- m_dynmat/dist9
- m_dynmat/dymfz9
- m_dynmat/dynmat_dq
- m_dynmat/ftgam
- m_dynmat/ftgam_init
- m_dynmat/ftifc_q2r
- m_dynmat/ftifc_r2q
- m_dynmat/gtdyn9
- m_dynmat/ifclo9
- m_dynmat/make_bigbox
- m_dynmat/massmult_and_breaksym
- m_dynmat/nanal9
- m_dynmat/q0dy3_apply
- m_dynmat/q0dy3_calc
- m_dynmat/symdm9
- m_dynmat/symdyma
- m_dynmat/sytens
- m_dynmat/wght9
- m_dynmat/wings3
ABINIT/m_dynmat [ Modules ]
NAME
m_dynmat
FUNCTION
This module provides low-level tools to operate on the dynamical matrix
COPYRIGHT
Copyright (C) 2014-2018 ABINIT group (XG, JCC, MJV, NH, RC, MVeithen, MM, MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
TODO
Use more explicative names for the procedures!
PARENTS
CHILDREN
SOURCE
24 #if defined HAVE_CONFIG_H 25 #include "config.h" 26 #endif 27 28 #include "abi_common.h" 29 30 module m_dynmat 31 32 use defs_basis 33 use m_profiling_abi 34 use m_errors 35 use m_linalg_interfaces 36 37 use m_fstrings, only : itoa, sjoin 38 use m_numeric_tools, only : wrap2_pmhalf, mkherm 39 use m_cgtools, only : fxphas_seq 40 use m_ewald, only : ewald9 41 42 implicit none 43 44 private 45 46 public :: asria_calc ! Calculate the correction for the Acoustic sum rule on 47 ! the InterAtomic Forces or on the dynamical matrix directly 48 public :: asria_corr ! Imposition of the Acoustic sum rule on the InterAtomic Forces 49 ! or on the dynamical matrix directly from the previously calculated d2asr 50 public :: asrprs ! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry 51 public :: cart29 ! Transform a second-derivative matrix from reduced coordinates to cartesian coordinates, and also 52 ! 1) add the ionic part of the effective charges, 53 ! 2) normalize the electronic dielectric tensor, and add the vacuum polarisation 54 public :: cart39 ! Transform a vector from reduced coordinates to cartesian coordinates, 55 ! taking into account the perturbation from which it was derived, 56 ! and also check the existence of the new values. 57 public :: d2cart_to_red ! Transform a second-derivative matrix 58 ! from cartesian to reduced coordinate. 59 public :: chkph3 ! Check the completeness of the dynamical matrix 60 public :: chneu9 ! Imposition of the Acoustic sum rule on the Effective charges 61 public :: d2sym3 ! Build (nearly) all the other matrix elements that can be build using symmetries. 62 public :: q0dy3_apply ! Takes care of the inclusion of the ewald q=0 term in the dynamical matrix 63 public :: q0dy3_calc ! Calculate the q=0 correction term to the dynamical matrix 64 public :: symdyma ! Symmetrize the dynamical matrices 65 public :: wings3 ! Suppress the wings of the cartesian 2DTE for which the diagonal element is not known 66 public :: asrif9 ! Imposes the Acoustic Sum Rule to Interatomic Forces 67 public :: make_bigbox ! Generates a Big Box of R points for the Fourier Transforms the dynamical matrix 68 public :: bigbx9 ! Helper functions that faciliates the generation of a Big Box containing 69 public :: canat9 ! From reduced to canonical coordinates 70 public :: canct9 ! Convert from canonical coordinates to cartesian coordinates 71 public :: chkrp9 ! Check if the rprim used for the definition of the unit cell (in the 72 ! inputs) are consistent with the rprim used in the routine generating the Big Box 73 public :: dist9 ! Compute the distance between atoms in the big box 74 public :: ftifc_q2r ! Fourier transform of the dynamical matrices to obtain interatomic forces (real space). 75 public :: ftifc_r2q ! Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space). 76 public :: dynmat_dq ! Compute the derivative D(q)/dq via Fourier transform of the interatomic forces 77 public :: ifclo9 ! Convert from cartesian coordinates to local coordinates 78 public :: wght9 ! Generates a weight to each R points of the Big Box and for each pair of atoms 79 public :: d3sym ! Given a set of calculated elements of the 3DTE matrix, 80 ! build (nearly) all the other matrix elements that can be build using symmetries. 81 public :: sytens ! Determines the set of irreductible elements of the non-linear optical susceptibility 82 ! and Raman tensors 83 public :: symdm9 ! Generate the dynamical matrix in the full BZ from the irreducible q-points. 84 public :: axial9 ! Generates the local coordinates system from the knowledge of the first vector (longitudinal) and 85 ! the ifc matrix in cartesian coordinates 86 public :: dymfz9 ! Multiply the dynamical matrix by a phase shift to account for normalized canonical coordinates. 87 public :: nanal9 ! Subtract/Add the non-analytical part from one dynamical matrix with number iqpt. 88 public :: gtdyn9 ! Generates a dynamical matrix from interatomic force constants and 89 ! long-range electrostatic interactions. 90 public :: dfpt_phfrq ! Diagonalize IFC(q), return phonon frequencies and eigenvectors. 91 ! If q is Gamma, the non-analytical behaviour can be included. 92 public :: massmult_and_breaksym ! Multiply IFC(q) by atomic masses. 93 94 ! TODO: Change name, 95 public :: ftgam 96 public :: ftgam_init 97 98 99 ! ************************************************************************* 100 101 contains
m_dynmat/asria_calc [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
asria_calc
FUNCTION
Calculate the correction for the Acoustic sum rule on the InterAtomic Forces or on the dynamical matrix directly
INPUTS
asr=(0 => no ASR, 1 or 2=> the diagonal element is modified to give the ASR,
5 => impose hermitian solution using lapack call)
d2cart=matrix of second derivatives of total energy, in cartesian coordinates
mpert =maximum number of ipert
natom=number of atom
OUTPUT
d2asr=matrix used to store the correction needed to fulfill the acoustic sum rule.
PARENTS
dfpt_gatherdy,m_ddb,m_effective_potential_file
CHILDREN
SOURCE
130 subroutine asria_calc(asr,d2asr,d2cart,mpert,natom) 131 132 133 !This section has been created automatically by the script Abilint (TD). 134 !Do not modify the following lines by hand. 135 #undef ABI_FUNC 136 #define ABI_FUNC 'asria_calc' 137 !End of the abilint section 138 139 implicit none 140 141 !Arguments ------------------------------- 142 !scalars 143 integer,intent(in) :: asr,mpert,natom 144 !arrays 145 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert) 146 real(dp),intent(out) :: d2asr(2,3,natom,3,natom) 147 148 !Local variables------------------------------- 149 !scalars 150 integer :: idir1,idir2,ii,ipert1,ipert2 151 integer :: constrank, imatelem, iconst, nconst, nd2_packed, info 152 !character(len=500) :: message 153 !arrays 154 integer, allocatable :: packingindex(:,:,:,:) 155 real(dp), allocatable :: constraints(:,:,:) 156 real(dp), allocatable :: d2cart_packed(:,:) 157 real(dp), allocatable :: singvals(:) 158 real(dp), allocatable :: constr_rhs(:,:) 159 real(dp), allocatable :: work(:,:),rwork(:) 160 161 ! ********************************************************************* 162 163 d2asr = zero 164 165 if (asr==0) return 166 167 !call wrtout(std_out,' asria_calc: calculation of the correction to the ASR for the interatomic forces.') 168 do ipert1=1,natom 169 do idir1=1,3 170 do idir2=1,3 171 172 ! Compute d2asr 173 do ipert2=1,natom 174 d2asr(:,idir1,ipert1,idir2,ipert1)=& 175 & d2asr(:,idir1,ipert1,idir2,ipert1)+& 176 & d2cart(:,idir1,ipert1,idir2,ipert2) 177 end do 178 end do 179 end do 180 end do 181 182 !holistic method: overwrite d2asr with hermitian solution 183 if (asr == 5) then 184 nconst = 9*natom 185 nd2_packed = 3*natom*(3*natom+1)/2 186 ABI_ALLOCATE(constraints,(2,nconst, nd2_packed)) 187 ABI_ALLOCATE(d2cart_packed,(2,nd2_packed)) 188 ABI_ALLOCATE(constr_rhs,(2,nd2_packed)) 189 ABI_ALLOCATE(singvals,(nconst)) 190 ABI_ALLOCATE(work,(2,3*nd2_packed)) 191 ABI_ALLOCATE(rwork,(5*nd2_packed)) 192 ABI_ALLOCATE(packingindex,(3,natom,3,natom)) 193 ii=1 194 packingindex=-1 195 do ipert2=1,natom 196 do idir2=1,3 197 do ipert1=1,ipert2-1 198 do idir1=1,3 199 packingindex(idir1,ipert1,idir2,ipert2) = ii 200 ii = ii+1 201 end do 202 end do 203 do idir1=1,idir2 204 packingindex(idir1,ipert2,idir2,ipert2) = ii 205 ii = ii+1 206 end do 207 end do 208 end do 209 ! setup constraint matrix 210 constraints = zero 211 do ipert1=1,natom 212 do idir1=1,3 213 do idir2=1,3 214 iconst = idir2+3*(idir1-1 + 3*(ipert1-1)) 215 ! set all atom forces, this component 216 do ipert2=1,natom 217 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 218 if (imatelem == -1) then 219 imatelem = packingindex(idir2,ipert2,idir1,ipert1) 220 end if 221 constraints(1,iconst,imatelem) = one 222 end do 223 end do 224 end do 225 end do 226 227 d2cart_packed = -999.0d0 228 do ipert2=1,natom 229 do idir2=1,3 230 do ipert1=1,natom 231 do idir1=1,3 232 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 233 if (imatelem == -1) cycle 234 d2cart_packed(:,imatelem) = d2cart(:,idir1,ipert1,idir2,ipert2) 235 end do 236 end do 237 end do 238 end do 239 constr_rhs = zero 240 constr_rhs(1,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(1,:)) 241 constr_rhs(2,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(2,:)) 242 243 ! lwork = 3*nd2_packed 244 call zgelss (nconst,nd2_packed,1,constraints,nconst,constr_rhs,nd2_packed,& 245 & singvals,-one,constrank,work,3*nd2_packed,rwork,info) 246 ABI_CHECK(info == 0, sjoin('zgelss returned:', itoa(info))) 247 248 ! unpack 249 do ipert2=1,natom 250 do idir2=1,3 251 do ipert1=1,natom 252 do idir1=1,3 253 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 254 if (imatelem == -1) then 255 imatelem = packingindex(idir2,ipert2,idir1,ipert1) 256 ! NOTE: should complex conjugate the correction below. 257 end if 258 d2asr(:,idir1,ipert1,idir2,ipert2) = constr_rhs(:,imatelem) 259 end do 260 end do 261 end do 262 end do 263 264 ABI_DEALLOCATE(constraints) 265 ABI_DEALLOCATE(d2cart_packed) 266 ABI_DEALLOCATE(singvals) 267 ABI_DEALLOCATE(constr_rhs) 268 ABI_DEALLOCATE(work) 269 ABI_DEALLOCATE(rwork) 270 ABI_DEALLOCATE(packingindex) 271 end if 272 273 end subroutine asria_calc
m_dynmat/asria_corr [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
asria_corr
FUNCTION
Imposition of the Acoustic sum rule on the InterAtomic Forces or on the dynamical matrix directly from the previously calculated d2asr
INPUTS
asr=(0 => no ASR, 1 or 2=> the diagonal element is modified to give the ASR,
5 => impose hermitian solution using lapack call)
d2asr=matrix used to store the correction needed to fulfill
the acoustic sum rule.
mpert =maximum number of ipert
natom=number of atom
OUTPUT
Input/Output: d2cart=matrix of second derivatives of total energy, in cartesian coordinates
PARENTS
ddb_elast,ddb_internalstr,dfpt_gatherdy,m_ddb,thmeig
CHILDREN
SOURCE
305 subroutine asria_corr(asr,d2asr,d2cart,mpert,natom) 306 307 308 !This section has been created automatically by the script Abilint (TD). 309 !Do not modify the following lines by hand. 310 #undef ABI_FUNC 311 #define ABI_FUNC 'asria_corr' 312 !End of the abilint section 313 314 implicit none 315 316 !Arguments ------------------------------- 317 !scalars 318 integer,intent(in) :: asr,mpert,natom 319 !arrays 320 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 321 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 322 323 !Local variables------------------------------- 324 !scalars 325 integer :: idir1,idir2,ipert1,ipert2 326 327 ! ********************************************************************* 328 329 if (asr==0) return 330 !call wrtout(std_out,' asria_corr: imposition of the ASR for the interatomic forces.','COLL') 331 332 ! Remove d2asr 333 do ipert2=1,natom 334 do idir2=1,3 335 do ipert1=1,natom 336 do idir1=1,3 337 d2cart(:,idir1,ipert1,idir2,ipert2)= d2cart(:,idir1,ipert1,idir2,ipert2) - d2asr(:,idir1,ipert1,idir2,ipert2) 338 end do 339 end do 340 end do 341 end do 342 343 end subroutine asria_corr
m_dynmat/asrif9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
asrif9
FUNCTION
Imposes the Acoustic Sum Rule to Interatomic Forces
INPUTS
asr= Option for the imposition of the ASR 0 => no ASR, 1 => modify "asymmetrically" the diagonal element 2 => modify "symmetrically" the diagonal element natom= Number of atoms in the unit cell nrpt= Number of R points in the Big Box rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) wghatm(natom,natom,nrpt)= Weight associated to the couple of atoms and the R vector atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces
OUTPUT
atmfrc(3,natom,3,natom,nrpt)= ASR-imposed Interatomic Forces
TODO
List of ouput should be included.
PARENTS
ddb_hybrid,m_ifc,m_tdep_abitypes
CHILDREN
SOURCE
2462 subroutine asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm) 2463 2464 2465 !This section has been created automatically by the script Abilint (TD). 2466 !Do not modify the following lines by hand. 2467 #undef ABI_FUNC 2468 #define ABI_FUNC 'asrif9' 2469 !End of the abilint section 2470 2471 implicit none 2472 2473 !Arguments ------------------------------- 2474 !scalars 2475 integer,intent(in) :: asr,natom,nrpt 2476 !arrays 2477 real(dp),intent(in) :: rpt(3,nrpt),wghatm(natom,natom,nrpt) 2478 real(dp),intent(inout) :: atmfrc(3,natom,3,natom,nrpt) 2479 2480 !Local variables ------------------------- 2481 !scalars 2482 integer :: found,ia,ib,irpt,izero,mu,nu 2483 real(dp) :: sumifc 2484 2485 ! ********************************************************************* 2486 2487 DBG_ENTER("COLL") 2488 2489 if(asr==1.or.asr==2)then 2490 2491 found=0 2492 ! Search for the R vector which is equal to ( 0 , 0 , 0 ) 2493 ! This vector leaves the atom a on itself ! 2494 do irpt=1,nrpt 2495 if (all(abs(rpt(:,irpt))<=1.0d-10)) then 2496 found=1 2497 izero=irpt 2498 end if 2499 if (found==1) exit 2500 end do 2501 2502 if(found==0)then 2503 MSG_BUG('Not able to find the vector R=(0,0,0).') 2504 end if 2505 2506 do mu=1,3 2507 do nu=1,3 2508 do ia=1,natom 2509 sumifc=zero 2510 do ib=1,natom 2511 2512 ! Get the sumifc of interatomic forces acting on the atom ia, 2513 ! either in a symmetrical manner, or an unsymmetrical one. 2514 if(asr==1)then 2515 do irpt=1,nrpt 2516 sumifc=sumifc+wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt) 2517 end do 2518 else if(asr==2)then 2519 do irpt=1,nrpt 2520 sumifc=sumifc+& 2521 & (wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)+& 2522 & wghatm(ia,ib,irpt)*atmfrc(nu,ia,mu,ib,irpt))/2 2523 end do 2524 end if 2525 end do 2526 2527 ! Correct the self-interaction in order to fulfill the ASR 2528 atmfrc(mu,ia,nu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)-sumifc 2529 if(asr==2)then 2530 atmfrc(nu,ia,mu,ia,izero)=atmfrc(mu,ia,nu,ia,izero) 2531 end if 2532 2533 end do 2534 end do 2535 end do 2536 end if 2537 2538 DBG_EXIT("COLL") 2539 2540 end subroutine asrif9
m_dynmat/asrprs [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
asrprs
FUNCTION
Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry
INPUTS
asr=(3 => 1D systems, all elements are modified to give ASR and
rotational symmetry)
(4 => 0D systems, all elements are modified to give ASR and
rotational symmetry)
asrflg=(1 => the correction to enforce asr is computed from
d2cart, but NOT applied;
2 => one uses the previously determined correction)
minvers=previously calculated inverted coefficient matrix
mpert =maximum number of ipert
natom=number of atom
rotinv=(1,2,3 => for linear systems along x,y,z
4 => non-linear molecule
xcart=cartesian coordinates of the ions
OUTPUT
(see side effects)
SIDE EFFECTS
Input/Output: d2cart=matrix of second derivatives of total energy, in cartesian coordinates minvers=inverse of the supermatrix for future application of the corrections
PARENTS
m_ddb
CHILDREN
SOURCE
385 subroutine asrprs(asr,asrflag,rotinv,uinvers,vtinvers,singular,d2cart,mpert,natom,xcart) 386 387 388 !This section has been created automatically by the script Abilint (TD). 389 !Do not modify the following lines by hand. 390 #undef ABI_FUNC 391 #define ABI_FUNC 'asrprs' 392 use interfaces_14_hidewrite 393 !End of the abilint section 394 395 implicit none 396 397 !Arguments ------------------------------------ 398 !scalars 399 integer,intent(in) :: asr,asrflag,mpert,natom,rotinv 400 !arrays 401 real(dp),intent(in) :: xcart(3,natom) 402 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 403 real(dp),intent(inout) :: singular(1:3*natom*(3*natom-1)/2) 404 real(dp),intent(inout) :: uinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2) 405 real(dp),intent(inout) :: vtinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2) 406 407 !Local variables------------------------------- 408 !scalars 409 integer :: column,idir1,idir2,ii,info,ipert1,ipert2,jj,n3,row,superdim 410 real(dp) :: rcond,test 411 ! real(dp) :: tau ! tau is present but commented out in this routine 412 character(len=500) :: message 413 !arrays 414 integer :: check(3,natom,3) 415 real(dp) :: tmp(natom,3,3),weightf(1:natom,1:natom) 416 real(dp),allocatable :: d2cartold(:,:,:,:,:),d2vecc(:),d2veccnew(:),d2vecr(:) 417 real(dp),allocatable :: d2vecrnew(:),superm(:,:),umatrix(:,:),vtmatrix(:) 418 real(dp),allocatable :: work(:) 419 420 ! ********************************************************************* 421 422 if(asr/=3 .and. asr/=4)then 423 write(message,'(3a,i0)')& 424 & 'The argument asr should be 3 or 4,',ch10,& 425 & 'however, asr = ',asr 426 MSG_BUG(message) 427 end if 428 429 if (asr==3.or.asr==4)then 430 write(message, '(a,a)' ) ch10, & 431 & 'asrprs: imposition of the ASR for the interatomic forces and rotational invariance' 432 call wrtout(std_out,message,'COLL') 433 end if 434 435 write(message,'(a,i6)')' asrflag is', asrflag 436 call wrtout(std_out,message,'COLL') 437 call wrtout(ab_out,message,'COLL') 438 439 !variables for the dimensions of the matrices 440 441 !n1=3*natom*(3*natom-1)/2 442 !n2=9*natom 443 n3=3*natom 444 445 superdim=9*natom*(natom-1)/2+n3 446 447 ABI_ALLOCATE(d2vecr,(1:superdim)) 448 ABI_ALLOCATE(d2vecc,(1:superdim)) 449 d2vecr=0d0 450 d2vecc=0d0 451 452 !should be changed set to delta function for debugging 453 weightf=1d0 454 !tau=1d-10 455 do ii=1, natom 456 ! do jj=1, ii-1 457 ! weightf(ii,jj)= & 458 ! & ((xcart(1,ii)-xcart(1,jj))**2+(xcart(2,ii)-xcart(2,jj))**2+(xcart(3,ii)-xcart(3,jj))**2)**tau 459 ! enddo 460 weightf(ii,ii)=0d0 461 end do 462 463 ABI_ALLOCATE(d2cartold,(2,3,mpert,3,mpert)) 464 465 d2cartold=d2cart 466 467 !setup vector with uncorrected derivatives 468 469 do ipert1=1, natom 470 do ipert2=1, ipert1-1 471 do idir1=1,3 472 do idir2=1,3 473 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 474 if(abs(d2cart(1,idir1,ipert1,idir2,ipert2))<1d-6)then 475 d2cart(1,idir1,ipert1,idir2,ipert2)=0d0 476 else 477 d2vecr(row)=4*weightf(ipert1,ipert2)*d2cart(1,idir1,ipert1,idir2,ipert2) 478 end if 479 if(abs(d2cart(2,idir1,ipert1,idir2,ipert2))<1d-6) then 480 d2cart(2,idir1,ipert1,idir2,ipert2)=0d0 481 else 482 d2vecc(row)=4*weightf(ipert1,ipert2)*d2cart(2,idir1,ipert1,idir2,ipert2) 483 end if 484 end do 485 end do 486 end do 487 end do 488 489 if(asrflag==1) then !calculate the pseudo-inverse of the supermatrix 490 ABI_ALLOCATE(superm,(1:superdim,1:superdim)) 491 492 superm=0d0 493 494 ! Setting up the supermatrix containing G, A, D 495 496 do ipert1=1, natom 497 do idir1=1, 3 498 ! Setting up G 499 idir2=mod(idir1,3)+1 500 row=3*(ipert1-1)+idir1 501 do ipert2=1, ipert1-1 502 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1 503 superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1) 504 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir2 505 superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2) 506 end do 507 do ipert2=ipert1+1, natom 508 column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir1-1)+rotinv 509 superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1) 510 column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir2-1)+rotinv 511 superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2) 512 end do 513 end do 514 do idir1=1, 3 515 ! Setting up D 516 idir2=mod(idir1,3)+1 517 ii=mod(idir1+1,3)+1 518 do ipert2=1, ipert1-1 519 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1 520 column=9*natom*(natom-1)/2+3*(ipert1-1)+idir1 521 superm(column,row)=superm(column,row)+xcart(idir2,ipert2)-xcart(idir2,ipert1) 522 column=9*natom*(natom-1)/2+3*(ipert1-1)+ii 523 superm(column,row)=superm(column,row)+xcart(ii,ipert1)-xcart(ii,ipert2) 524 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+rotinv 525 column=9*natom*(natom-1)/2+3*(ipert2-1)+idir1 526 superm(column,row)=superm(column,row)+xcart(idir2,ipert1)-xcart(idir2,ipert2) 527 column=9*natom*(natom-1)/2+3*(ipert2-1)+ii 528 superm(column,row)=superm(column,row)+xcart(ii,ipert2)-xcart(ii,ipert1) 529 end do 530 ! Setting up A 531 do idir2=1, 3 532 do ipert2=1, ipert1-1 533 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 534 row=n3+column 535 superm(column,row)=4*weightf(ipert1,ipert2) 536 end do 537 end do 538 end do 539 end do 540 541 ! calculate the pseudo-inverse of the supermatrix 542 543 ABI_ALLOCATE(work,(1:6*superdim)) 544 ABI_ALLOCATE(vtmatrix,(1:superdim)) 545 ABI_ALLOCATE(umatrix,(1:superdim,1:superdim)) 546 547 ! singular value decomposition of superm 548 549 call dgesvd('A','O',superdim,superdim,superm,superdim,singular,umatrix,superdim, & 550 & vtmatrix, 1, work,6*superdim,info) 551 ABI_CHECK(info == 0, sjoin('dgesvd returned:', itoa(info))) 552 553 ABI_DEALLOCATE(vtmatrix) 554 ABI_DEALLOCATE(work) 555 556 write(message, '(a,es16.8,es16.8)' )& 557 & ' Largest and smallest values from svd', singular(1), singular(superdim) 558 call wrtout(std_out,message,'COLL') 559 call wrtout(ab_out,message,'COLL') 560 561 ! Invert U and V**T, orthogonal matrices 562 563 do ii=1, superdim 564 do jj=1, superdim 565 uinvers(ii,jj)=umatrix(jj,ii) 566 vtinvers(ii,jj)=superm(jj,ii) 567 end do 568 end do 569 570 ABI_DEALLOCATE(umatrix) 571 ABI_DEALLOCATE(superm) 572 573 write(message,'(a,a)')' asrprs: done with asrflag 1', ch10 574 call wrtout(std_out,message,'COLL') 575 576 end if !asrflag=1 577 578 if(asrflag==2) then 579 580 ABI_ALLOCATE(d2vecrnew,(1:superdim)) 581 ABI_ALLOCATE(d2veccnew,(1:superdim)) 582 583 ! Calculate V**T**-1 Sigma**-1 U**-1 *rhs 584 585 do ii=1, superdim 586 d2vecrnew(ii)=0d0 587 d2veccnew(ii)=0d0 588 do jj=1, superdim 589 d2vecrnew(ii)=d2vecrnew(ii)+uinvers(ii,jj)*d2vecr(jj) 590 d2veccnew(ii)=d2veccnew(ii)+uinvers(ii,jj)*d2vecc(jj) 591 end do 592 end do 593 594 rcond=1d-10*singular(1) 595 do ii=1, superdim 596 if(singular(ii)>rcond) then 597 d2vecrnew(ii)=d2vecrnew(ii)/singular(ii) 598 d2veccnew(ii)=d2veccnew(ii)/singular(ii) 599 else 600 d2vecrnew(ii)=0d0 601 d2veccnew(ii)=0d0 602 end if 603 end do 604 605 do ii=1, superdim 606 d2vecr(ii)=0d0 607 d2vecc(ii)=0d0 608 do jj=1, superdim 609 d2vecr(ii)=d2vecr(ii)+vtinvers(ii,jj)*d2vecrnew(jj) 610 d2vecc(ii)=d2vecc(ii)+vtinvers(ii,jj)*d2veccnew(jj) 611 end do 612 end do 613 614 ! Store vector back into the matrix of 2nd order derivates 615 616 do ipert1=1, natom 617 do ipert2=1, ipert1-1 618 do idir1=1,3 619 do idir2=1,3 620 row=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 621 d2cart(1,idir1,ipert1,idir2,ipert2)=d2vecr(row) 622 d2cart(2,idir1,ipert1,idir2,ipert2)=d2vecc(row) 623 d2cart(1,idir2,ipert2,idir1,ipert1)=d2vecr(row) 624 d2cart(2,idir2,ipert2,idir1,ipert1)=d2vecc(row) 625 end do 626 end do 627 end do 628 end do 629 630 ABI_DEALLOCATE(d2vecrnew) 631 ABI_DEALLOCATE(d2veccnew) 632 633 check=0 634 635 do ipert1=1, natom 636 do idir1=1, 3 637 do idir2=1, 3 638 d2cart(1,idir1,ipert1,idir2,ipert1)=0d0 639 d2cart(2,idir1,ipert1,idir2,ipert1)=0d0 640 tmp(ipert1,idir1,idir2)=0d0 641 do ipert2=1, natom 642 if(ipert2/=ipert1) then 643 tmp(ipert1,idir1,idir2)=tmp(ipert1,idir1,idir2) & 644 & -d2cart(1,idir1,ipert1,idir2,ipert2) & 645 & -d2cart(1,idir2,ipert2,idir1,ipert1) 646 end if 647 end do 648 end do 649 end do 650 end do 651 652 do ipert1=1, natom 653 do idir1=1, 3 654 do idir2=1, 3 655 d2cart(1,idir1,ipert1,idir2,ipert1)=tmp(ipert1,idir1,idir2)/2 656 d2cart(1,idir2,ipert1,idir1,ipert1)=d2cart(1,idir1,ipert1,idir2,ipert1) 657 end do 658 end do 659 end do 660 661 write(std_out,*) 'this should all be zero' 662 663 do ipert1=1, natom 664 do idir1=1, 3 665 do idir2=1, 3 666 test=0d0 667 do ipert2=1, natom 668 test=test+d2cart(1,idir1,ipert1,idir2,ipert2)+d2cart(1,idir2,ipert2,idir1,ipert1) 669 end do 670 write(std_out,'(i3,i3,i3,es11.3)') idir1,ipert1,idir2,test 671 672 write(message, '(i3,i3,i3,es11.3)' ) idir1,ipert1,idir2,test 673 call wrtout(ab_out,message,'COLL') 674 675 end do 676 end do 677 end do 678 679 write(std_out,*) 'these as well' 680 do ipert2=1, natom 681 do idir1=1, 3 682 do idir2=1, 3 683 test=0d0 684 do ipert1=1, natom 685 test=test+d2cart(1,idir1,ipert1,idir2,ipert2) 686 end do 687 write(std_out,'(i3,i3,i3,i3,es11.3)') idir1,ipert1,idir2,ipert2,test 688 end do 689 end do 690 end do 691 692 write(message,'(a,a)')' asrprs: done with asrflag 2', ch10 693 call wrtout(std_out,message,'COLL') 694 695 end if !ends asrflag=2 696 697 ABI_DEALLOCATE(d2vecr) 698 ABI_DEALLOCATE(d2vecc) 699 ABI_DEALLOCATE(d2cartold) 700 701 end subroutine asrprs
m_dynmat/axial9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
axial9
FUNCTION
Generates the local coordinates system from the knowledge of the first vector (longitudinal) and the ifc matrix in cartesian coordinates
INPUTS
ifccar(3,3)= matrix of interatomic force constants in cartesian coordinates vect1(3)= cartesian coordinates of the first local vector
OUTPUT
vect2(3)= cartesian coordinates of the second local vector vect3(3)= cartesian coordinates of the third local vector
PARENTS
m_ifc
CHILDREN
SOURCE
5039 subroutine axial9(ifccar,vect1,vect2,vect3) 5040 5041 5042 !This section has been created automatically by the script Abilint (TD). 5043 !Do not modify the following lines by hand. 5044 #undef ABI_FUNC 5045 #define ABI_FUNC 'axial9' 5046 !End of the abilint section 5047 5048 implicit none 5049 5050 !Arguments ------------------------------- 5051 !arrays 5052 real(dp),intent(in) :: ifccar(3,3),vect1(3) 5053 real(dp),intent(out) :: vect2(3),vect3(3) 5054 5055 !Local variables ------------------------- 5056 !scalars 5057 integer :: flag,ii,itrial,jj 5058 real(dp) :: innorm,scprod 5059 !arrays 5060 real(dp) :: work(3) 5061 5062 ! ********************************************************************* 5063 5064 do jj=1,3 5065 work(jj)=zero 5066 do ii=1,3 5067 work(jj)=work(jj)+ifccar(jj,ii)*vect1(ii) 5068 end do 5069 end do 5070 5071 flag=0 5072 do itrial=1,4 5073 scprod=zero 5074 do ii=1,3 5075 scprod=scprod+work(ii)*vect1(ii) 5076 end do 5077 5078 do ii=1,3 5079 work(ii)=work(ii)-vect1(ii)*scprod 5080 end do 5081 5082 scprod=zero 5083 do ii=1,3 5084 scprod=scprod+work(ii)**2 5085 end do 5086 5087 if(scprod<1.0d-10)then 5088 work(1:3)=zero 5089 if(itrial>1)work(itrial-1)=1.0_dp 5090 else 5091 flag=1 5092 end if 5093 5094 if(flag==1)exit 5095 end do 5096 5097 innorm=scprod**(-0.5_dp) 5098 do ii=1,3 5099 vect2(ii)=work(ii)*innorm 5100 end do 5101 5102 vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2) 5103 vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3) 5104 vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1) 5105 5106 end subroutine axial9
m_dynmat/bigbx9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
bigbx9
FUNCTION
Generation of a Big Box containing all the R points in the cartesian real space needed to Fourier Transforms the dynamical matrix into its corresponding interatomic force.
INPUTS
brav= Bravais Lattice (1 or -1=S.C.;2=F.C.C.;3=BCC;4=Hex.) choice= if 0, simply count nrpt ; if 1, checks that the input mrpt is the same as nrpt, and generate rpt(3,mrpt) mrpt=dimension of rpt ngqpt(3)= Numbers used to generate the q points to sample the Brillouin zone using an homogeneous grid nqshft= number of q-points in the repeated cell for the Brillouin zone sampling When nqshft is not 1, but 2 or 4 (only other allowed values), the limits for the big box have to be extended by a factor of 2. rprim(3,3)= Normalized coordinates in real space !!! IS THIS CORRECT?
OUTPUT
cell= (nrpt,3) Give the index of the the cell and irpt nprt= Number of R points in the Big Box rpt(3,mrpt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) (output only if choice=1)
PARENTS
m_dynmat
CHILDREN
SOURCE
2661 subroutine bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt) 2662 2663 2664 !This section has been created automatically by the script Abilint (TD). 2665 !Do not modify the following lines by hand. 2666 #undef ABI_FUNC 2667 #define ABI_FUNC 'bigbx9' 2668 !End of the abilint section 2669 2670 implicit none 2671 2672 !Arguments ------------------------------- 2673 !scalars 2674 integer,intent(in) :: brav,choice,mrpt,nqshft 2675 integer,intent(out) :: nrpt 2676 !arrays 2677 integer,intent(in) :: ngqpt(3) 2678 real(dp),intent(in) :: rprim(3,3) 2679 real(dp),intent(out) :: rpt(3,mrpt) 2680 integer,intent(out) :: cell(3,mrpt) 2681 2682 !Local variables ------------------------- 2683 !In some cases, the atoms coordinates are not packed in the 2684 ! [0,1]^3 cube. Then, the parameter "buffer" might be increased, 2685 !to search relevant pairs of atoms in bigger boxes than usual. 2686 !scalars 2687 integer,parameter :: buffer=1 2688 integer :: irpt,lim1,lim2,lim3,lqshft,r1,r2,r3 2689 character(len=500) :: msg 2690 2691 ! ********************************************************************* 2692 2693 lqshft=1 2694 if(nqshft/=1)lqshft=2 2695 2696 2697 !Simple Cubic Lattice 2698 if (abs(brav)==1) then 2699 lim1=((ngqpt(1))+1)*lqshft+buffer 2700 lim2=((ngqpt(2))+1)*lqshft+buffer 2701 lim3=((ngqpt(3))+1)*lqshft+buffer 2702 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1) 2703 if(choice/=0)then 2704 if (nrpt/=mrpt) then 2705 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt 2706 MSG_BUG(msg) 2707 end if 2708 irpt=0 2709 do r1=-lim1,lim1 2710 do r2=-lim2,lim2 2711 do r3=-lim3,lim3 2712 irpt=irpt+1 2713 rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3) 2714 rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3) 2715 rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3) 2716 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2717 end do 2718 end do 2719 end do 2720 end if 2721 2722 ! Face Centered Cubic Lattice 2723 else if (brav==2) then 2724 lim1=((ngqpt(1)+3)/4)*lqshft+buffer 2725 lim2=((ngqpt(2)+3)/4)*lqshft+buffer 2726 lim3=((ngqpt(3)+3)/4)*lqshft+buffer 2727 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*4 2728 if(choice/=0)then 2729 if (nrpt/=mrpt) then 2730 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt 2731 MSG_BUG(msg) 2732 end if 2733 irpt=0 2734 do r1=-lim1,lim1 2735 do r2=-lim2,lim2 2736 do r3=-lim3,lim3 2737 irpt=irpt+4 2738 rpt(1,irpt-3)=r1 2739 rpt(2,irpt-3)=r2 2740 rpt(3,irpt-3)=r3 2741 rpt(1,irpt-2)=r1 2742 rpt(2,irpt-2)=r2+0.5 2743 rpt(3,irpt-2)=r3+0.5 2744 rpt(1,irpt-1)=r1+0.5 2745 rpt(2,irpt-1)=r2 2746 rpt(3,irpt-1)=r3+0.5 2747 rpt(1,irpt)=r1+0.5 2748 rpt(2,irpt)=r2+0.5 2749 rpt(3,irpt)=r3 2750 !TEST_AM 2751 ! cell(irpt-3,1)=r1;cell(irpt-3,2)=r2;cell(irpt-3,3)=r3 2752 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2753 end do 2754 end do 2755 end do 2756 end if 2757 2758 ! Body Centered Cubic Lattice 2759 else if (brav==3) then 2760 lim1=((ngqpt(1)+3)/4)*lqshft+buffer 2761 lim2=((ngqpt(2)+3)/4)*lqshft+buffer 2762 lim3=((ngqpt(3)+3)/4)*lqshft+buffer 2763 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*2 2764 if(choice/=0)then 2765 if(nrpt/=mrpt) then 2766 write(msg,'(2(a,i0))')' nrpt= ',nrpt,' is not equal to mrpt= ',mrpt 2767 MSG_BUG(msg) 2768 end if 2769 irpt=0 2770 do r1=-lim1,lim1 2771 do r2=-lim2,lim2 2772 do r3=-lim3,lim3 2773 irpt=irpt+2 2774 rpt(1,irpt-1)=r1 2775 rpt(2,irpt-1)=r2 2776 rpt(3,irpt-1)=r3 2777 rpt(1,irpt)=r1+0.5 2778 rpt(2,irpt)=r2+0.5 2779 rpt(3,irpt)=r3+0.5 2780 !TEST_AM 2781 ! cell(irpt-1,1)=r1;cell(irpt-1,2)=r2;cell(irpt-1,3)=r3 2782 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2783 end do 2784 end do 2785 end do 2786 end if 2787 2788 ! Hexagonal Lattice 2789 else if (brav==4) then 2790 lim1=(ngqpt(1)+1)*lqshft+buffer 2791 lim2=(ngqpt(2)+1)*lqshft+buffer 2792 lim3=((ngqpt(3)/2)+1)*lqshft+buffer 2793 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1) 2794 if(choice/=0)then 2795 if(nrpt/=mrpt)then 2796 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt=',mrpt 2797 MSG_BUG(msg) 2798 end if 2799 irpt=0 2800 do r1=-lim1,lim1 2801 do r2=-lim2,lim2 2802 do r3=-lim3,lim3 2803 irpt=irpt+1 2804 rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3) 2805 rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3) 2806 rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3) 2807 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2808 end do 2809 end do 2810 end do 2811 end if 2812 2813 else 2814 write(msg,'(a,i0,a)')' The value of brav= ',brav,' is not allowed (should be -1, 1, 2 or 4).' 2815 MSG_BUG(msg) 2816 end if 2817 2818 end subroutine bigbx9
m_dynmat/canat9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
canat9
FUNCTION
Transforms an atom whose coordinates (xred*rprim) would not be in the chosen unit cell used to generate the interatomic forces to its correspondent (rcan) in canonical coordinates.
INPUTS
brav= Bravais Lattice (1 or -1=S.C.;2=F.C.C.;3=BCC;4=Hex.) natom= Number of atoms in the unit cell rprim(3,3)= Normalized coordinates of primitive vectors
OUTPUT
rcan(3,natom) = Atomic position in canonical coordinates trans(3,natom) = Atomic translations : xred = rcan + trans
PARENTS
m_ifc
CHILDREN
SOURCE
2849 subroutine canat9(brav,natom,rcan,rprim,trans,xred) 2850 2851 2852 !This section has been created automatically by the script Abilint (TD). 2853 !Do not modify the following lines by hand. 2854 #undef ABI_FUNC 2855 #define ABI_FUNC 'canat9' 2856 use interfaces_14_hidewrite 2857 !End of the abilint section 2858 2859 implicit none 2860 2861 !Arguments ------------------------------- 2862 !scalars 2863 integer,intent(in) :: brav,natom 2864 !arrays 2865 real(dp),intent(in) :: rprim(3,3),xred(3,natom) 2866 real(dp),intent(out) :: rcan(3,natom),trans(3,natom) 2867 2868 !Local variables ------------------------- 2869 !scalars 2870 integer :: found,iatom,ii 2871 character(len=500) :: message 2872 !arrays 2873 real(dp) :: dontno(3,4),rec(3),rok(3),shift(3),tt(3) 2874 2875 ! ********************************************************************* 2876 2877 DBG_ENTER("COLL") 2878 2879 !Normalization of the cartesian atomic coordinates 2880 !If not normalized : rcan(i) <- rcan(i) * acell(i) 2881 do iatom=1,natom 2882 rcan(:,iatom)=xred(1,iatom)*rprim(:,1)+& 2883 & xred(2,iatom)*rprim(:,2)+& 2884 & xred(3,iatom)*rprim(:,3) 2885 end do 2886 2887 !Study of the different cases for the Bravais lattice: 2888 2889 !Simple Cubic Lattice 2890 if (abs(brav)==1) then 2891 2892 do iatom=1,natom 2893 2894 ! Canon will produces these coordinate transformations 2895 ! (Note : here we still use reduced coordinates ) 2896 call wrap2_pmhalf(xred(1,iatom),rok(1),shift(1)) 2897 call wrap2_pmhalf(xred(2,iatom),rok(2),shift(2)) 2898 call wrap2_pmhalf(xred(3,iatom),rok(3),shift(3)) 2899 2900 ! New coordinates : rcan 2901 rcan(:,iatom)=rok(1)*rprim(:,1)+& 2902 & rok(2)*rprim(:,2)+& 2903 & rok(3)*rprim(:,3) 2904 ! Translations between New and Old coordinates 2905 tt(:)=xred(1,iatom)*rprim(:,1)+& 2906 & xred(2,iatom)*rprim(:,2)+& 2907 & xred(3,iatom)*rprim(:,3) 2908 trans(:,iatom)=tt(:)-rcan(:,iatom) 2909 end do 2910 2911 else if (brav==2) then 2912 ! Face Centered Lattice 2913 ! Special possible translations in the F.C.C. case 2914 dontno(:,:)=zero 2915 dontno(2,2)=0.5_dp 2916 dontno(3,2)=0.5_dp 2917 dontno(1,3)=0.5_dp 2918 dontno(3,3)=0.5_dp 2919 dontno(1,4)=0.5_dp 2920 dontno(2,4)=0.5_dp 2921 do iatom=1,natom 2922 found=0 2923 do ii=1,4 2924 if (found==1) exit 2925 ! Canon will produces these coordinate transformations 2926 call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1)) 2927 call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2)) 2928 call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3)) 2929 ! In the F.C.C., ABS[ Ri ] + ABS[ Rj ] < or = 1/2 2930 ! The equal signs has been treated using a tolerance parameter 2931 ! not to have twice the same point in the unit cell ! 2932 rok(1)=rok(1)-1.0d-10 2933 rok(2)=rok(2)-2.0d-10 2934 rok(3)=rok(3)-5.0d-10 2935 if (abs(rok(1))+abs(rok(2))<=0.5_dp) then 2936 if (abs(rok(1))+abs(rok(3))<=0.5_dp) then 2937 if (abs(rok(2))+abs(rok(3))<=0.5_dp) then 2938 tt(:)=rcan(:,iatom) 2939 ! New coordinates : rcan 2940 rcan(1,iatom)=rok(1)+1.0d-10 2941 rcan(2,iatom)=rok(2)+2.0d-10 2942 rcan(3,iatom)=rok(3)+5.0d-10 2943 ! Translations between New and Old coordinates 2944 trans(:,iatom)=tt(:)-rcan(:,iatom) 2945 found=1 2946 end if 2947 end if 2948 end if 2949 end do 2950 end do 2951 2952 else if (brav==3) then 2953 ! Body Centered Cubic Lattice 2954 ! Special possible translations in the B.C.C. case 2955 2956 dontno(:,1)=zero 2957 dontno(:,2)=0.5_dp 2958 do iatom=1,natom 2959 found=0 2960 do ii=1,2 2961 if (found==1) exit 2962 ! Canon will produce these coordinate transformations 2963 call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1)) 2964 call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2)) 2965 call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3)) 2966 ! In the F.C.C., ABS[ Ri ] < or = 1/2 2967 ! and ABS[ R1 ] + ABS[ R2 ] + ABS[ R3 ] < or = 3/4 2968 ! The equal signs has been treated using a tolerance parameter 2969 ! not to have twice the same point in the unit cell ! 2970 rok(1)=rok(1)-1.0d-10 2971 rok(2)=rok(2)-2.0d-10 2972 rok(3)=rok(3)-5.0d-10 2973 if(abs(rok(1))+abs(rok(2))+abs(rok(3))<=0.75_dp)then 2974 if ( abs(rok(1))<=0.5_dp .and.& 2975 & abs(rok(2))<=0.5_dp .and.& 2976 & abs(rok(3))<=0.5_dp ) then 2977 tt(:)=rcan(:,iatom) 2978 ! New coordinates : rcan 2979 rcan(1,iatom)=rok(1)+1.0d-10 2980 rcan(2,iatom)=rok(2)+2.0d-10 2981 rcan(3,iatom)=rok(3)+5.0d-10 2982 ! Translations between New and Old coordinates 2983 trans(:,iatom)=tt(:)-rcan(:,iatom) 2984 found=1 2985 end if 2986 end if 2987 end do 2988 end do 2989 2990 else if (brav==4) then 2991 ! Hexagonal Lattice 2992 ! In this case, it is easier first to work in reduced coordinates space ! 2993 do iatom=1,natom 2994 ! Passage from the reduced space to the "lozenge" cell 2995 rec(1)=xred(1,iatom)-0.5_dp 2996 rec(2)=xred(2,iatom)-0.5_dp 2997 rec(3)=xred(3,iatom) 2998 ! Canon will produces these coordinate transformations 2999 call wrap2_pmhalf(rec(1),rok(1),shift(1)) 3000 call wrap2_pmhalf(rec(2),rok(2),shift(2)) 3001 call wrap2_pmhalf(rec(3),rok(3),shift(3)) 3002 rec(1)=rok(1)+0.5_dp 3003 rec(2)=rok(2)+0.5_dp 3004 rec(3)=rok(3) 3005 ! Passage in Cartesian Normalized Coordinates 3006 rcan(:,iatom)=rec(1)*rprim(:,1)+& 3007 & rec(2)*rprim(:,2)+& 3008 & rec(3)*rprim(:,3) 3009 ! Use of a tolerance parameter not to have twice the same point 3010 ! in the unit cell ! 3011 rcan(1,iatom)=rcan(1,iatom)-1.0d-10 3012 rcan(2,iatom)=rcan(2,iatom)-2.0d-10 3013 ! Passage to the honey-com hexagonal unit cell ! 3014 if (rcan(1,iatom)>0.5_dp) then 3015 rcan(1,iatom)=rcan(1,iatom)-1.0_dp 3016 end if 3017 if (rcan(1,iatom)>zero.and.rcan(1,iatom)+sqrt(3.0_dp)*rcan& 3018 & (2,iatom)>1.0_dp) then 3019 rcan(1,iatom)=rcan(1,iatom)-0.5_dp 3020 rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp 3021 end if 3022 if (rcan(1,iatom)<=zero.and.sqrt(3.0_dp)*rcan(2,iatom)-rcan& 3023 & (1,iatom)>1.0_dp) then 3024 rcan(1,iatom)=rcan(1,iatom)+0.5_dp 3025 rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp 3026 end if 3027 ! Translations between New and Old coordinates 3028 tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)& 3029 & +xred(3,iatom)*rprim(:,3) 3030 trans(:,iatom)=tt(:)-rcan(:,iatom) 3031 end do 3032 3033 ! End of the possible cases for brav : -1, 1, 2, 4. 3034 else 3035 3036 write(message, '(a,i0,a,a,a)' )& 3037 & 'The required value of brav=',brav,' is not available.',ch10,& 3038 & 'It should be -1, 1,2 or 4 .' 3039 MSG_BUG(message) 3040 end if 3041 3042 call wrtout(std_out,' Canonical Atomic Coordinates ',"COLL") 3043 3044 do iatom=1,natom 3045 write(message, '(a,i5,3es18.8)' )' atom',iatom,rcan(1,iatom),rcan(2,iatom),rcan(3,iatom) 3046 call wrtout(std_out,message,'COLL') 3047 end do 3048 3049 DBG_EXIT("COLL") 3050 3051 end subroutine canat9
m_dynmat/canct9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
canct9
FUNCTION
Convert from canonical coordinates to cartesian coordinates a vector defined by its index=ib+natom*(irpt-1)
INPUTS
acell(3)=length scales by which rprim is to be multiplied gprim(3,3)=dimensionless primitive translations in reciprocal space index= index of the atom natom=number of atoms in unit cell nrpt= Number of R points in the Big Box rcan(3,natom)=canonical coordinates of atoms rprim(3,3)=dimensionless primitive translations in real space rpt(3,nrpt)=canonical coordinates of the points in the BigBox.
OUTPUT
ib=number of the atom in the unit cell irpt= number of the unit cell to which belong the atom rcart(3)=cartesian coordinate of the atom indexed by index.
PARENTS
ddb_hybrid,m_ifc
CHILDREN
SOURCE
3087 subroutine canct9(acell,gprim,ib,index,irpt,natom,nrpt,rcan,rcart,rprim,rpt) 3088 3089 3090 !This section has been created automatically by the script Abilint (TD). 3091 !Do not modify the following lines by hand. 3092 #undef ABI_FUNC 3093 #define ABI_FUNC 'canct9' 3094 !End of the abilint section 3095 3096 implicit none 3097 3098 !Arguments ------------------------------- 3099 !scalars 3100 integer,intent(in) :: index,natom,nrpt 3101 integer,intent(out) :: ib,irpt 3102 !arrays 3103 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3) 3104 real(dp),intent(in) :: rpt(3,nrpt) 3105 real(dp),intent(out) :: rcart(3) 3106 3107 !Local variables ------------------------- 3108 !scalars 3109 integer :: jj 3110 !arrays 3111 real(dp) :: xred(3) 3112 3113 ! ********************************************************************* 3114 3115 irpt=(index-1)/natom+1 3116 ib=index-natom*(irpt-1) 3117 3118 !Transform the canonical coordinates to reduced coord. 3119 do jj=1,3 3120 xred(jj)=gprim(1,jj)*(rpt(1,irpt)+rcan(1,ib))& 3121 & +gprim(2,jj)*(rpt(2,irpt)+rcan(2,ib))& 3122 & +gprim(3,jj)*(rpt(3,irpt)+rcan(3,ib)) 3123 end do 3124 3125 !Then to cartesian coordinates (here the position of the atom b) 3126 do jj=1,3 3127 rcart(jj)=xred(1)*acell(1)*rprim(jj,1)+& 3128 & xred(2)*acell(2)*rprim(jj,2)+& 3129 & xred(3)*acell(3)*rprim(jj,3) 3130 end do 3131 3132 end subroutine canct9
m_dynmat/cart29 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
cart29
FUNCTION
Transform a second-derivative matrix from reduced coordinates to cartesian coordinates, and also 1) add the ionic part of the effective charges, 2) normalize the electronic dielectric tensor, and add the vacuum polarisation
INPUTS
blkflg(3,mpert,3,mpert,nblok)=
( 1 if the element of the dynamical matrix has been calculated ;
0 otherwise )
blkval(2,3,mpert,3,mpert,nblok)=DDB values
gprimd(3,3)=basis vector in the reciprocal space
iblok=number of the blok that will be transformed
mpert =maximum number of ipert
natom=number of atom
nblok=number of blocks (dimension of blkflg and blkval)
ntypat=number of atom types
rprimd(3,3)=basis vector in the real space
typat(natom)=integer label of each type of atom (1,2,...)
ucvol=unit cell volume
zion(ntypat)=charge corresponding to the atom type
OUTPUT
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian
2DTE matrix has been calculated correctly ; 0 otherwise )
d2cart(2,3,mpert,3,mpert)=
dynamical matrix, effective charges, dielectric tensor,....
all in cartesian coordinates
PARENTS
dfpt_gatherdy,m_ddb
CHILDREN
SOURCE
748 subroutine cart29(blkflg,blkval,carflg,d2cart,& 749 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 750 751 752 !This section has been created automatically by the script Abilint (TD). 753 !Do not modify the following lines by hand. 754 #undef ABI_FUNC 755 #define ABI_FUNC 'cart29' 756 !End of the abilint section 757 758 implicit none 759 760 !Arguments ------------------------------- 761 !scalars 762 integer,intent(in) :: iblok,mpert,natom,nblok,ntypat 763 real(dp),intent(in) :: ucvol 764 !arrays 765 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),typat(natom) 766 integer,intent(out) :: carflg(3,mpert,3,mpert) 767 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3) 768 real(dp),intent(in) :: zion(ntypat) 769 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 770 771 !Local variables ------------------------- 772 !scalars 773 integer :: idir1,idir2,ii,ipert1,ipert2 774 !arrays 775 integer :: flg1(3),flg2(3) 776 real(dp) :: vec1(3),vec2(3) 777 778 ! ********************************************************************* 779 780 !First, copy the data blok in place. 781 d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok) 782 783 !Cartesian coordinates transformation (in two steps) 784 !First step 785 do ipert1=1,mpert 786 do ipert2=1,mpert 787 do ii=1,2 788 do idir1=1,3 789 do idir2=1,3 790 vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2) 791 ! Note here blkflg 792 flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok) 793 end do 794 call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2) 795 do idir2=1,3 796 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) 797 ! And here carflg 798 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2) 799 end do 800 end do 801 end do 802 end do 803 end do 804 805 !Second step 806 do ipert1=1,mpert 807 do ipert2=1,mpert 808 do ii=1,2 809 do idir2=1,3 810 do idir1=1,3 811 vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2) 812 ! Note here carflg 813 flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2) 814 end do 815 call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2) 816 do idir1=1,3 817 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) 818 ! And here carflg again 819 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1) 820 end do 821 end do 822 end do 823 end do 824 end do 825 826 !For the dielectric tensor, takes into account the volume 827 !of the unit cell, and add the unit matrix (polarization of the vacuum) 828 do idir1=1,3 829 do idir2=1,3 830 do ii=1,2 831 d2cart(ii,idir1,natom+2,idir2,natom+2)=& 832 & -four_pi/ucvol*d2cart(ii,idir1,natom+2,idir2,natom+2) 833 end do 834 end do 835 end do 836 837 do idir1=1,3 838 d2cart(1,idir1,natom+2,idir1,natom+2)=& 839 & 1.0_dp+d2cart(1,idir1,natom+2,idir1,natom+2) 840 end do 841 842 !Add the ionic charges to delta z to get the effective charges 843 do ipert1=1,natom 844 do idir1=1,3 845 d2cart(1,idir1,ipert1,idir1,natom+2)=& 846 & zion(typat(ipert1))+d2cart(1,idir1,ipert1,idir1,natom+2) 847 end do 848 end do 849 do ipert2=1,natom 850 do idir2=1,3 851 d2cart(1,idir2,natom+2,idir2,ipert2)=& 852 & zion(typat(ipert2))+d2cart(1,idir2,natom+2,idir2,ipert2) 853 end do 854 end do 855 856 !For the piezoelectric tensor, takes into account the volume of the unit cell 857 do ipert2=natom+3,natom+4 858 do idir1=1,3 859 do idir2=1,3 860 do ii=1,2 861 d2cart(ii,idir1,natom+2,idir2,ipert2)=& 862 & (1.0_dp/ucvol)*d2cart(ii,idir1,natom+2,idir2,ipert2) 863 d2cart(ii,idir2,ipert2,idir1,natom+2)=& 864 & (1.0_dp/ucvol)*d2cart(ii,idir2,ipert2,idir1,natom+2) 865 end do 866 end do 867 end do 868 end do 869 870 end subroutine cart29
m_dynmat/cart39 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
cart39
FUNCTION
Transform a vector from reduced coordinates to cartesian coordinates, taking into account the perturbation from which it was derived, and also check the existence of the new values.
INPUTS
flg1(3)=tell if information of each component of vec1 is valid gprimd(3,3)=basis vector in the reciprocal space ipert=number of the perturbation natom=number of atom rprimd(3,3)=basis vector in the real space vec1(3)=input vector, in reduced coordinates
OUTPUT
flg2(3)=tell if information of each component of vec2 is valid vec2(3)=output vector, in cartesian coordinates
PARENTS
dfpt_gatherdy,m_ddb,m_dynmat
CHILDREN
SOURCE
904 subroutine cart39(flg1,flg2,gprimd,ipert,natom,rprimd,vec1,vec2) 905 906 907 !This section has been created automatically by the script Abilint (TD). 908 !Do not modify the following lines by hand. 909 #undef ABI_FUNC 910 #define ABI_FUNC 'cart39' 911 !End of the abilint section 912 913 implicit none 914 915 !Arguments ------------------------------- 916 !scalars 917 integer,intent(in) :: ipert,natom 918 !arrays 919 integer,intent(in) :: flg1(3) 920 integer,intent(out) :: flg2(3) 921 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3),vec1(3) 922 real(dp),intent(out) :: vec2(3) 923 924 !Local variables ------------------------- 925 !scalars 926 integer :: idir,ii 927 928 ! ********************************************************************* 929 930 !Treat phonon-type perturbation 931 if(ipert>=1.and.ipert<=natom)then 932 933 do idir=1,3 934 vec2(idir)=zero 935 flg2(idir)=1 936 do ii=1,3 937 if(abs(gprimd(idir,ii))>1.0d-10)then 938 if(flg1(ii)==1)then 939 vec2(idir)=vec2(idir)+gprimd(idir,ii)*vec1(ii) 940 else 941 flg2(idir)=0 942 end if 943 end if 944 end do 945 if(flg2(idir)==0)vec2(idir)=zero 946 end do 947 948 ! Treat electric field perturbation 949 else if(ipert==natom+2) then 950 ! OCL SCALAR 951 do idir=1,3 952 vec2(idir)=zero 953 flg2(idir)=1 954 ! OCL SCALAR 955 do ii=1,3 956 if(abs(rprimd(idir,ii))>1.0d-10)then 957 if(flg1(ii)==1)then 958 vec2(idir)=vec2(idir)+rprimd(idir,ii)*vec1(ii)/two_pi 959 else 960 flg2(idir)=0 961 end if 962 end if 963 end do 964 if(flg2(idir)==0)vec2(idir)=zero 965 end do 966 967 ! Treat other perturbations 968 else 969 do idir=1,3 970 vec2(idir)=vec1(idir) 971 flg2(idir)=flg1(idir) 972 end do 973 end if 974 975 end subroutine cart39
m_dynmat/chkph3 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
chkph3
FUNCTION
Check the completeness of the dynamical matrix and eventually send a warning
INPUTS
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian 2DTE matrix has been calculated correctly ; 0 otherwise ) idir = direction of the eventual electric field mpert =maximum number of ipert natom=number of atoms in unit cell
OUTPUT
eventually send a warning message
PARENTS
respfn
CHILDREN
SOURCE
1168 subroutine chkph3(carflg,idir,mpert,natom) 1169 1170 1171 !This section has been created automatically by the script Abilint (TD). 1172 !Do not modify the following lines by hand. 1173 #undef ABI_FUNC 1174 #define ABI_FUNC 'chkph3' 1175 use interfaces_14_hidewrite 1176 !End of the abilint section 1177 1178 implicit none 1179 1180 !Arguments ------------------------------- 1181 !scalars 1182 integer,intent(in) :: idir,mpert,natom 1183 !arrays 1184 integer,intent(in) :: carflg(3,mpert,3,mpert) 1185 1186 !Local variables ------------------------- 1187 !scalars 1188 integer :: idir1,idir2,ipert1,ipert2,send 1189 character(len=500) :: message 1190 1191 ! ********************************************************************* 1192 1193 send=0 1194 1195 !Check the elements of the analytical part of the dynamical matrix 1196 do ipert2=1,natom 1197 do idir2=1,3 1198 do ipert1=1,natom 1199 do idir1=1,3 1200 if(carflg(idir1,ipert1,idir2,ipert2)==0)then 1201 send=1 1202 end if 1203 end do 1204 end do 1205 end do 1206 end do 1207 1208 !If some electric field is present 1209 if(idir/=0)then 1210 1211 ! Check the dielectric constant 1212 if(carflg(idir,natom+2,idir,natom+2)==0)then 1213 send=1 1214 end if 1215 1216 ! Check the effective charges 1217 do ipert1=1,natom 1218 do idir1=1,3 1219 if(carflg(idir1,ipert1,idir,natom+2)==0)then 1220 send=1 1221 end if 1222 end do 1223 end do 1224 1225 end if 1226 1227 !If needed, send the message 1228 if(send==1)then 1229 write(message, '(a,a,a,a)' )& 1230 & ' chkph3 : WARNING -',ch10,& 1231 & ' The dynamical matrix was incomplete :',& 1232 & ' phonon frequencies may be wrong ...' 1233 call wrtout(std_out,message,'COLL') 1234 call wrtout(ab_out,message,'COLL') 1235 end if 1236 1237 end subroutine chkph3
m_dynmat/chkrp9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
chkrp9
FUNCTION
Check if the rprim used for the definition of the unit cell (in the inputs) are consistent with the rprim used in the routine generating the Big Box needed to generate the interatomic forces.
INPUTS
brav=bravais lattice (1 or -1=simple lattice,2=face centered lattice, 3=centered lattice,4=hexagonal lattice) rprimd(3,3)=dimensional primitive translations for real space (bohr)
OUTPUT
(only checking)
PARENTS
m_ifc
CHILDREN
SOURCE
3161 subroutine chkrp9(brav,rprim) 3162 3163 3164 !This section has been created automatically by the script Abilint (TD). 3165 !Do not modify the following lines by hand. 3166 #undef ABI_FUNC 3167 #define ABI_FUNC 'chkrp9' 3168 !End of the abilint section 3169 3170 implicit none 3171 3172 !Arguments ------------------------------- 3173 !scalars 3174 integer,intent(in) :: brav 3175 !arrays 3176 real(dp),intent(in) :: rprim(3,3) 3177 3178 !Local variables ------------------------- 3179 !scalars 3180 integer :: ii,jj 3181 character(len=500) :: message 3182 3183 ! ********************************************************************* 3184 3185 if (abs(brav)==1) then 3186 ! Simple Cubic Lattice No condition in this case ! 3187 continue 3188 3189 else if (brav==2) then 3190 ! Face Centered Lattice 3191 do ii=1,3 3192 do jj=1,3 3193 if ( ( ii==jj .and. abs(rprim(ii,jj))>tol10) .or.& 3194 & ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then 3195 write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3196 & 'The input variable rprim does not correspond to the',ch10,& 3197 & 'fixed rprim to be used with brav=2 and ifcflag=1 :',ch10,& 3198 & ' 0 1/2 1/2',ch10,& 3199 & ' 1/2 0 1/2',ch10,& 3200 & ' 1/2 1/2 0 ',ch10,& 3201 & 'Action: rebuild your DDB by using the latter rprim.' 3202 MSG_ERROR(message) 3203 end if 3204 end do 3205 end do 3206 3207 else if (brav==3) then 3208 ! Body Centered Cubic Lattice 3209 do ii=1,3 3210 do jj=1,3 3211 if ( ( ii==jj .and. abs(rprim(ii,jj)+.5_dp)>tol10) .or.& 3212 & ( ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then 3213 write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3214 & 'The input variable rprim does not correspond to the',ch10,& 3215 & 'fixed rprim to be used with brav=3 and ifcflag=1 :',ch10,& 3216 & ' -1/2 1/2 1/2',ch10,& 3217 & ' 1/2 -1/2 1/2',ch10,& 3218 & ' 1/2 1/2 -1/2',ch10,& 3219 & 'Action: rebuild your DDB by using the latter rprim.' 3220 MSG_ERROR(message) 3221 end if 3222 end do 3223 end do 3224 3225 else if (brav==4) then 3226 ! Hexagonal Lattice 3227 if (abs(rprim(1,1)-1.0_dp)>tol10 .or. & 3228 & abs(rprim(3,3)-1.0_dp)>tol10 .or. & 3229 & abs(rprim(2,1) )>tol10 .or. & 3230 & abs(rprim(3,1) )>tol10 .or. & 3231 & abs(rprim(1,3) )>tol10 .or. & 3232 & abs(rprim(2,3) )>tol10 .or. & 3233 & abs(rprim(3,2) )>tol10 .or. & 3234 & abs(rprim(1,2)+0.5_dp)>tol10 .or. & 3235 & abs(rprim(2,2)-0.5_dp*sqrt(3.0_dp))>tol10 ) then 3236 write(message, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3237 & 'The input variable rprim does not correspond to the',ch10,& 3238 & 'fixed rprim to be used with brav=4 and ifcflag=1 :',ch10,& 3239 & ' 1 0 0',ch10,& 3240 & ' -1/2 sqrt[3]/2 0',ch10,& 3241 & ' 0 0 1',ch10,& 3242 & 'Action: rebuild your DDB by using the latter rprim.' 3243 MSG_ERROR(message) 3244 end if 3245 3246 else 3247 3248 write(message, '(a,i4,a,a,a,a,a)' )& 3249 & 'The value of brav=',brav,' is not allowed.',ch10,& 3250 & 'Only -1, 1,2,3 or 4 are allowed.',ch10,& 3251 & 'Action: change the value of brav in your input file.' 3252 MSG_ERROR(message) 3253 end if 3254 3255 end subroutine chkrp9
m_dynmat/chneu9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
chneu9
FUNCTION
Imposition of the Acoustic sum rule on the Effective charges and suppress the imaginary part of the dynamical matrix
INPUTS
chneut=(0 => no ASR, 1 => equal repartition,2 => weighted repartition ) mpert =maximum number of ipert natom=number of atom ntypat=number of types of atoms in unit cell selectz=selection of some parts of the effective charge tensor attached to one atom. (0=> no selection, 1=> trace only, 2=> symmetric part only) typat(natom)=type of the atom zion(ntypat)=atomic charge for every type of atom
SIDE EFFECTS
Input/Output d2cart=matrix of second derivatives of total energy, in cartesian coordinates
PARENTS
dfpt_gatherdy,m_ddb
CHILDREN
SOURCE
1272 subroutine chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion) 1273 1274 1275 !This section has been created automatically by the script Abilint (TD). 1276 !Do not modify the following lines by hand. 1277 #undef ABI_FUNC 1278 #define ABI_FUNC 'chneu9' 1279 use interfaces_14_hidewrite 1280 !End of the abilint section 1281 1282 implicit none 1283 1284 !Arguments ------------------------------- 1285 !scalars 1286 integer,intent(in) :: chneut,mpert,natom,ntypat,selectz 1287 !arrays 1288 integer,intent(in) :: typat(natom) 1289 real(dp),intent(in) :: zion(ntypat) 1290 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 1291 1292 !Local variables ------------------------- 1293 !scalars 1294 integer :: idir1,idir2,ii,ipert1,ipert2 1295 character(len=500) :: message 1296 !arrays 1297 real(dp) :: sumwght(2) 1298 real(dp),allocatable :: wghtat(:) 1299 1300 ! ********************************************************************* 1301 1302 ABI_ALLOCATE(wghtat,(natom)) 1303 1304 !In case of acoustic sum rule imposition, compute the weights on each atom. 1305 if (chneut==1)then 1306 1307 ! The weight is the same for all atom 1308 do ipert1=1,natom 1309 wghtat(ipert1)=1./natom 1310 end do 1311 1312 else if (chneut==2) then 1313 1314 ! The weight is proportional to the diagonal electronic screening charge of the atom 1315 sumwght(1)=zero 1316 do ipert1=1,natom 1317 wghtat(ipert1)=zero 1318 do idir1=1,3 1319 wghtat(ipert1)=wghtat(ipert1)+& 1320 & d2cart(1,idir1,ipert1,idir1,natom+2)+& 1321 & d2cart(1,idir1,natom+2,idir1,ipert1)-2*zion(typat(ipert1)) 1322 end do 1323 sumwght(1)=sumwght(1)+wghtat(ipert1) 1324 end do 1325 1326 ! Normalize the weights to unity 1327 do ipert1=1,natom 1328 wghtat(ipert1)=wghtat(ipert1)/sumwght(1) 1329 end do 1330 end if 1331 1332 !Calculation of the violation of the charge neutrality 1333 !and imposition of the charge neutrality condition 1334 if (chneut/=0)then 1335 write(message, '(a,a,a,a,a,a,a)' )& 1336 & ' The violation of the charge neutrality conditions',ch10,& 1337 & ' by the effective charges is as follows :',ch10,& 1338 & ' atom electric field',ch10,& 1339 & ' displacement direction ' 1340 call wrtout(ab_out,message,'COLL') 1341 do idir1=1,3 1342 do idir2=1,3 1343 do ii=1,2 1344 sumwght(ii)=zero 1345 do ipert1=1,natom 1346 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir2,natom+2) 1347 end do 1348 do ipert1=1,natom 1349 d2cart(ii,idir1,ipert1,idir2,natom+2)=& 1350 & d2cart(ii,idir1,ipert1,idir2,natom+2)-sumwght(ii)*wghtat(ipert1) 1351 end do 1352 end do 1353 write(message, '(i8,i16,2f16.6)' ) idir1,idir2,sumwght(1),sumwght(2) 1354 call wrtout(ab_out,message,'COLL') 1355 end do 1356 end do 1357 write(message, '(a)' )' ' 1358 call wrtout(ab_out,message,'COLL') 1359 1360 ! The same for the symmetrical part 1361 do idir1=1,3 1362 do idir2=1,3 1363 do ii=1,2 1364 sumwght(ii)=zero 1365 do ipert2=1,natom 1366 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir2,ipert2) 1367 end do 1368 do ipert2=1,natom 1369 d2cart(ii,idir1,natom+2,idir2,ipert2)=& 1370 & d2cart(ii,idir1,natom+2,idir2,ipert2)-sumwght(ii)*wghtat(ipert2) 1371 end do 1372 end do 1373 end do 1374 end do 1375 end if 1376 1377 !Selection of the trace of the effective charge tensor attached to each atom 1378 if(selectz==1)then 1379 do ipert1=1,natom 1380 do ii=1,2 1381 sumwght(ii)=zero 1382 do idir1=1,3 1383 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir1,natom+2) 1384 end do 1385 do idir1=1,3 1386 do idir2=1,3 1387 d2cart(ii,idir1,ipert1,idir2,natom+2)=zero 1388 end do 1389 end do 1390 do idir1=1,3 1391 d2cart(ii,idir1,ipert1,idir1,natom+2)=sumwght(ii)/3.0_dp 1392 end do 1393 end do 1394 end do 1395 ! Do the same for the symmetrical part of d2cart 1396 do ipert2=1,natom 1397 do ii=1,2 1398 sumwght(ii)=zero 1399 do idir1=1,3 1400 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir1,ipert2) 1401 end do 1402 do idir1=1,3 1403 do idir2=1,3 1404 d2cart(ii,idir1,natom+2,idir2,ipert2)=zero 1405 end do 1406 end do 1407 do idir1=1,3 1408 d2cart(ii,idir1,natom+2,idir1,ipert2)=sumwght(ii)/3.0_dp 1409 end do 1410 end do 1411 end do 1412 end if 1413 1414 !Selection of the symmetric part of the effective charge tensor attached to each atom 1415 if(selectz==2)then 1416 do ipert1=1,natom 1417 do ii=1,2 1418 do idir1=1,3 1419 do idir2=1,3 1420 sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)& 1421 & +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp 1422 d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii) 1423 d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii) 1424 end do 1425 end do 1426 end do 1427 end do 1428 ! Do the same for the symmetrical part of d2cart 1429 do ipert1=1,natom 1430 do ii=1,2 1431 do idir1=1,3 1432 do idir2=1,3 1433 sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)& 1434 & +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp 1435 d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii) 1436 d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii) 1437 end do 1438 end do 1439 end do 1440 end do 1441 end if 1442 1443 !Write the effective charge tensor 1444 write(message, '(a,a,a,a,a,a,a)' )& 1445 & ' Effective charge tensors after ',ch10,& 1446 & ' imposition of the charge neutrality,',ch10,& 1447 & ' and eventual restriction to some part :',ch10,& 1448 & ' atom displacement ' 1449 call wrtout(ab_out,message,'COLL') 1450 1451 do ipert1=1,natom 1452 do idir1=1,3 1453 write(message, '(2i10,3es16.6)' )ipert1,idir1,& 1454 & (d2cart(1,idir1,ipert1,idir2,natom+2),idir2=1,3) 1455 call wrtout(ab_out,message,'COLL') 1456 end do 1457 end do 1458 1459 !Zero the imaginary part of the dynamical matrix 1460 write(message, '(a)' )' Now, the imaginary part of the dynamical matrix is zeroed ' 1461 call wrtout(ab_out,message,'COLL') 1462 call wrtout(std_out,message,'COLL') 1463 1464 do ipert1=1,natom 1465 do ipert2=1,natom 1466 do idir1=1,3 1467 do idir2=1,3 1468 d2cart(2,idir1,ipert1,idir2,ipert2)=zero 1469 end do 1470 end do 1471 end do 1472 end do 1473 1474 ABI_DEALLOCATE(wghtat) 1475 1476 end subroutine chneu9
m_dynmat/d2cart_to_red [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
d2cart_to_red
FUNCTION
Transform a second-derivative matrix from cartesian
coordinates to reduced coordinate. Also,
1) remove the ionic part of the effective charges,
2) remove the vacuum polarisation from the dielectric tensor
and scale it with the unit cell volume
In short, does the inverse operation of cart29.
INPUTS
d2cart(2,3,mpert,3,mpert)=
second-derivative matrix in cartesian coordinates
gprimd(3,3)=basis vector in the reciprocal space
rprimd(3,3)=basis vector in the real space
mpert =maximum number of ipert
natom=number of atom
OUTPUT
d2red(2,3,mpert,3,mpert)=
second-derivative matrix in reduced coordinates
PARENTS
ddb_interpolate
CHILDREN
SOURCE
1012 subroutine d2cart_to_red(d2cart, d2red, gprimd, rprimd, mpert, natom, & 1013 & ntypat,typat,ucvol,zion) 1014 1015 1016 !This section has been created automatically by the script Abilint (TD). 1017 !Do not modify the following lines by hand. 1018 #undef ABI_FUNC 1019 #define ABI_FUNC 'd2cart_to_red' 1020 !End of the abilint section 1021 1022 implicit none 1023 1024 !Arguments ------------------------------- 1025 !scalars 1026 integer,intent(in) :: mpert,natom,ntypat 1027 real(dp),intent(in) :: ucvol 1028 !arrays 1029 integer,intent(in) :: typat(natom) 1030 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert) 1031 real(dp),intent(out) :: d2red(2,3,mpert,3,mpert) 1032 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 1033 real(dp),intent(in) :: zion(ntypat) 1034 1035 !Local variables ------------------------- 1036 !scalars 1037 integer :: idir1,idir2,ii,ipert1,ipert2 1038 real(dp) :: fac 1039 !arrays 1040 integer :: flg1(3),flg2(3) 1041 real(dp) :: vec1(3),vec2(3) 1042 1043 ! ********************************************************************* 1044 1045 flg1 = one 1046 flg2 = one 1047 1048 d2red = d2cart 1049 1050 !Remove the ionic charges to z to get the change in effective charges 1051 do ipert1=1,natom 1052 do idir1=1,3 1053 d2red(1,idir1,ipert1,idir1,natom+2)=& 1054 & d2red(1,idir1,ipert1,idir1,natom+2) - zion(typat(ipert1)) 1055 end do 1056 end do 1057 do ipert2=1,natom 1058 do idir2=1,3 1059 d2red(1,idir2,natom+2,idir2,ipert2)=& 1060 & d2red(1,idir2,natom+2,idir2,ipert2) - zion(typat(ipert2)) 1061 end do 1062 end do 1063 1064 ! Remove the vacuum polarizability from the dielectric tensor 1065 do idir1=1,3 1066 d2red(1,idir1,natom+2,idir1,natom+2)=& 1067 & d2red(1,idir1,natom+2,idir1,natom+2) - 1.0_dp 1068 end do 1069 1070 ! Scale the dielectric tensor with the volue of the unit cell 1071 do idir1=1,3 1072 do idir2=1,3 1073 do ii=1,2 1074 d2red(ii,idir1,natom+2,idir2,natom+2)=& 1075 & - (ucvol / four_pi) * d2red(ii,idir1,natom+2,idir2,natom+2) 1076 end do 1077 end do 1078 end do 1079 1080 !For the piezoelectric tensor, takes into account the volume of the unit cell 1081 do ipert2=natom+3,natom+4 1082 do idir1=1,3 1083 do idir2=1,3 1084 do ii=1,2 1085 d2red(ii,idir1,natom+2,idir2,ipert2)=& 1086 & (ucvol)*d2red(ii,idir1,natom+2,idir2,ipert2) 1087 d2red(ii,idir2,ipert2,idir1,natom+2)=& 1088 & (ucvol)*d2red(ii,idir2,ipert2,idir1,natom+2) 1089 end do 1090 end do 1091 end do 1092 end do 1093 1094 ! Reduced coordinates transformation (in two steps) 1095 ! Note that rprimd and gprimd are swapped, compared to what cart39 expects 1096 ! A factor of (2pi) ** 2 is added to transform the electric field perturbations 1097 1098 !First step 1099 do ipert1=1,mpert 1100 fac = one; if (ipert1==natom+2) fac = two_pi ** 2 1101 1102 do ipert2=1,mpert 1103 do ii=1,2 1104 do idir1=1,3 1105 do idir2=1,3 1106 vec1(idir2)=d2red(ii,idir1,ipert1,idir2,ipert2) 1107 end do 1108 ! Transform vector from cartesian to reduced coordinates 1109 call cart39(flg1,flg2,transpose(rprimd),ipert1,natom,transpose(gprimd),vec1,vec2) 1110 do idir2=1,3 1111 d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) * fac 1112 end do 1113 end do 1114 end do 1115 end do 1116 end do 1117 1118 !Second step 1119 do ipert1=1,mpert 1120 do ipert2=1,mpert 1121 fac = one; if (ipert2==natom+2) fac = two_pi ** 2 1122 1123 do ii=1,2 1124 do idir2=1,3 1125 do idir1=1,3 1126 vec1(idir1)=d2red(ii,idir1,ipert1,idir2,ipert2) 1127 end do 1128 ! Transform vector from cartestian to reduced coordinates 1129 call cart39(flg1,flg2,transpose(rprimd),ipert2,natom,transpose(gprimd),vec1,vec2) 1130 do idir1=1,3 1131 d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) * fac 1132 end do 1133 end do 1134 end do 1135 end do 1136 end do 1137 1138 1139 end subroutine d2cart_to_red
m_dynmat/d2sym3 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
d2sym3
FUNCTION
Given a set of calculated elements of the 2DTE matrix d2, build (nearly) all the other matrix elements that can be build using symmetries. 1. Perform first some completion by symmetrisation (exchange) over the two defining perturbations 2. For each element, uses every symmetry, and build the element, in case EITHER all the needed elements are available, OR the only missing is itself OR the perturbation is the electric field, in a diamond symmetry (the last case was coded rather dirty)
INPUTS
indsym(4,nsym,natom)=indirect indexing array : for each isym,iatom, fourth element is label of atom into which iatom is sent by INVERSE of symmetry operation isym; first three elements are the primitive translations which must be subtracted after the transformation to get back to the original unit cell. mpert =maximum number of ipert natom= number of atoms nsym=number of space group symmetries qpt(3)=wavevector of the perturbation symq(4,2,nsym)= (integer) three first numbers define the G vector ; fourth number is zero if the q-vector is not preserved, is 1 otherwise second index is one without time-reversal symmetry, two with time-reversal symmetry symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space) timrev=1 if the time-reversal symmetry preserves the wavevector, modulo a reciprocal lattice vector, timrev=0 otherwise
SIDE EFFECTS
Input/Output
d2(2,3,mpert,3,mpert)= matrix of the 2DTE
blkflg(3,mpert,3,mpert)= ( 1 if the element of the dynamical
matrix has been calculated ; 0 otherwise)
NOTES
The complete search would be to have the possibility of a set of missing elements. See notes of July 2, 1994, in the blue notebook 'computer codes' The partial solution adopted here takes into account some mirror symmetries as well as the tetrahedral symmetry of the diamond lattice On 010331, replaced the loops up to mpert by loops up to natom+2, because of a crash bug under Windows. However, the problem lies likely in the use of the indsym array.
PARENTS
completeperts,m_ddb,respfn
CHILDREN
SOURCE
1541 subroutine d2sym3(blkflg,d2,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev) 1542 1543 1544 !This section has been created automatically by the script Abilint (TD). 1545 !Do not modify the following lines by hand. 1546 #undef ABI_FUNC 1547 #define ABI_FUNC 'd2sym3' 1548 !End of the abilint section 1549 1550 implicit none 1551 1552 !Arguments ------------------------------- 1553 !scalars 1554 integer,intent(in) :: mpert,natom,nsym,timrev 1555 !arrays 1556 integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym) 1557 integer,intent(in),target :: symrec(3,3,nsym),symrel(3,3,nsym) 1558 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 1559 real(dp),intent(in) :: qpt(3) 1560 real(dp),intent(inout) :: d2(2,3,mpert,3,mpert) 1561 1562 !Local variables ------------------------- 1563 !scalars 1564 logical, parameter :: do_final_sym=.true. 1565 logical :: qzero 1566 integer :: exch12,found,idir1,idir2,idisy1,idisy2,ipert1,ipert2 1567 integer :: ipesy1,ipesy2,isgn,isym,ithree,itirev,noccur,nsym_used,quit,quit1 1568 real(dp) :: arg1,arg2,im,norm,re,sumi,sumr,xi,xr 1569 !arrays 1570 integer,pointer :: sym1_(:,:,:),sym2_(:,:,:) 1571 real(dp),allocatable :: d2tmp1(:,:,:),d2tmp2(:,:,:),d2work(:,:,:,:,:) 1572 1573 ! ********************************************************************* 1574 1575 qzero=(qpt(1)**2+qpt(2)**2+qpt(3)**2<tol16) 1576 1577 !Here look after exchange of 1 and 2 axis, 1578 !for electric field in diamond symmetry 1579 exch12=0 1580 if (qzero) then 1581 do isym=1,nsym 1582 exch12=1 1583 if(symrel(1,1,isym)/=0)exch12=0 1584 if(symrel(1,2,isym)/=1)exch12=0 1585 if(symrel(1,3,isym)/=0)exch12=0 1586 if(symrel(2,1,isym)/=1)exch12=0 1587 if(symrel(2,2,isym)/=0)exch12=0 1588 if(symrel(2,3,isym)/=0)exch12=0 1589 if(symrel(3,1,isym)/=0)exch12=0 1590 if(symrel(3,2,isym)/=0)exch12=0 1591 if(symrel(3,3,isym)/=1)exch12=0 1592 ! if(exch12==1) write(std_out,*)' d2sym3 : found exchange 1 2 =',isym 1593 if(exch12==1)exit 1594 end do 1595 end if 1596 1597 !Exchange of perturbations 1598 1599 !Consider two cases : either time-reversal symmetry 1600 !conserves the wavevector, or not 1601 if(timrev==0)then 1602 1603 ! do ipert1=1,mpert See notes 1604 do ipert1=1,min(natom+2,mpert) 1605 do idir1=1,3 1606 1607 ! Since the matrix is hermitian, the diagonal elements are real 1608 d2(2,idir1,ipert1,idir1,ipert1)=zero 1609 1610 ! do ipert2=1,mpert See notes 1611 do ipert2=1,min(natom+2,mpert) 1612 do idir2=1,3 1613 1614 ! If an element exists 1615 if(blkflg(idir1,ipert1,idir2,ipert2)==1)then 1616 1617 ! Either complete the symmetric missing element 1618 if(blkflg(idir2,ipert2,idir1,ipert1)==0)then 1619 1620 d2(1,idir2,ipert2,idir1,ipert1)= d2(1,idir1,ipert1,idir2,ipert2) 1621 d2(2,idir2,ipert2,idir1,ipert1)=-d2(2,idir1,ipert1,idir2,ipert2) 1622 1623 blkflg(idir2,ipert2,idir1,ipert1)=1 1624 1625 ! Or symmetrize (the matrix is hermitian) in case both exists 1626 ! (Note : this opportunity has been disabled for more 1627 ! obvious search for bugs in the code ) 1628 ! else 1629 ! sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2) 1630 ! sumi=d2(1,idir2,ipert2,idir1,ipert1)-d2(1,idir1,ipert1,idir2,ipert2) 1631 ! d2(1,idir2,ipert2,idir1,ipert1)=half*sumr 1632 ! d2(1,idir1,ipert1,idir2,ipert2)=half*sumr 1633 ! d2(2,idir2,ipert2,idir1,ipert1)=half*sumi 1634 ! d2(2,idir1,ipert1,idir2,ipert2)=-half*sumi 1635 end if 1636 end if 1637 1638 end do 1639 end do 1640 1641 end do 1642 end do 1643 1644 ! Here, case with time-reversal symmetry 1645 else 1646 1647 ! do ipert1=1,mpert See notes 1648 do ipert1=1,min(natom+2,mpert) 1649 do idir1=1,3 1650 ! do ipert2=1,mpert See notes 1651 do ipert2=1,min(natom+2,mpert) 1652 do idir2=1,3 1653 d2(2,idir1,ipert1,idir2,ipert2)=zero 1654 1655 ! If an element exists 1656 if(blkflg(idir1,ipert1,idir2,ipert2)==1)then 1657 1658 ! Either complete the symmetric missing element 1659 if(blkflg(idir2,ipert2,idir1,ipert1)==0)then 1660 1661 d2(1,idir2,ipert2,idir1,ipert1)=d2(1,idir1,ipert1,idir2,ipert2) 1662 blkflg(idir2,ipert2,idir1,ipert1)=1 1663 1664 ! Or symmetrize (the matrix is hermitian) in case both exists 1665 ! (Note : this opportunity has been disabled for more 1666 ! obvious search for bugs in the code ) 1667 ! else 1668 ! sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2) 1669 ! d2(1,idir2,ipert2,idir1,ipert1)=half*sumr 1670 ! d2(1,idir1,ipert1,idir2,ipert2)=half*sumr 1671 end if 1672 1673 end if 1674 end do 1675 end do 1676 end do 1677 end do 1678 end if 1679 1680 !Big Big Loop : symmetrize three times, because 1681 !of some cases in which one element is not yet available 1682 !at the first pass, and even at the second one ! 1683 do ithree=1,3 1684 1685 ! Big loop on all elements 1686 ! do ipert1=1,mpert See notes 1687 do ipert1=1,min(natom+2,mpert) 1688 1689 ! Select the symmetries according to pertubation 1 1690 if (ipert1<=natom)then 1691 sym1_ => symrec 1692 else 1693 sym1_ => symrel 1694 end if 1695 1696 do idir1=1,3 1697 ! do ipert2=1,mpert See notes 1698 do ipert2=1,min(natom+2,mpert) 1699 1700 ! Select the symmetries according to pertubation 2 1701 if (ipert2<=natom)then 1702 sym2_ => symrec 1703 else 1704 sym2_ => symrel 1705 end if 1706 1707 do idir2=1,3 1708 1709 ! Will get element (idir1,ipert1,idir2,ipert2) 1710 ! so this element should not yet be present ... 1711 if(blkflg(idir1,ipert1,idir2,ipert2)/=1)then 1712 1713 d2(1,idir1,ipert1,idir2,ipert2)=zero 1714 d2(2,idir1,ipert1,idir2,ipert2)=zero 1715 1716 ! Loop on all symmetries, including time-reversal 1717 quit1=0 1718 do isym=1,nsym 1719 do itirev=1,2 1720 isgn=3-2*itirev 1721 1722 if(symq(4,itirev,isym)/=0)then 1723 found=1 1724 1725 ! Here select the symmetric of ipert1 1726 if(ipert1<=natom)then 1727 ipesy1=indsym(4,isym,ipert1) 1728 else if(ipert1==(natom+2).and.qzero)then 1729 ipesy1=ipert1 1730 else 1731 found=0 1732 end if 1733 1734 ! Here select the symmetric of ipert2 1735 if(ipert2<=natom)then 1736 ipesy2=indsym(4,isym,ipert2) 1737 else if(ipert2==(natom+2).and.qzero)then 1738 ipesy2=ipert2 1739 else 1740 found=0 1741 end if 1742 1743 ! Now that a symmetric perturbation has been obtained, 1744 ! including the expression of the symmetry matrix, see 1745 ! if the symmetric values are available 1746 if( found==1 ) then 1747 1748 sumr=zero 1749 sumi=zero 1750 noccur=0 1751 quit=0 1752 do idisy1=1,3 1753 do idisy2=1,3 1754 if(sym1_(idir1,idisy1,isym)/=0 .and. sym2_(idir2,idisy2,isym)/=0 )then 1755 if(blkflg(idisy1,ipesy1,idisy2,ipesy2)==1)then 1756 sumr=sumr+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*& 1757 & d2(1,idisy1,ipesy1,idisy2,ipesy2) 1758 sumi=sumi+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*& 1759 & d2(2,idisy1,ipesy1,idisy2,ipesy2) 1760 1761 ! Here, in case the symmetric of the element 1762 ! is the element, or the symmetric with 1763 ! respect to permutation of perturbations 1764 ! (some more conditions on the time-reversal 1765 ! symmetry must be fulfilled although) 1766 else if( idisy1==idir1 .and. ipesy1==ipert1& 1767 & .and. idisy2==idir2 .and. ipesy2==ipert2& 1768 & .and.(isgn==1 .or. timrev==1 & 1769 & .or. (idir1==idir2 .and. ipert1==ipert2)))& 1770 & then 1771 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1772 else if( idisy1==idir2 .and. ipesy1==ipert2& 1773 & .and. idisy2==idir1 .and. ipesy2==ipert1& 1774 & .and.(isgn==-1 .or. timrev==1& 1775 & .or. (idir1==idir2 .and. ipert1==ipert2)))& 1776 & then 1777 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1778 1779 ! Here, electric field case 1780 else if( exch12==1 .and. & 1781 & ipert1==natom+2 .and. ipert2==natom+2& 1782 & .and.(( idisy1+idir1 ==3 & 1783 & .and. idisy2==3 .and. idir2==3)& 1784 & .or. ( idisy1+idir2 ==3& 1785 & .and. idisy2==3 .and. idir1==3)& 1786 & .or. ( idisy2+idir2 ==3& 1787 & .and. idisy1==3 .and. idir1==3)& 1788 & .or. ( idisy2+idir1 ==3& 1789 & .and. idisy1==3 .and. idir2==3)))& 1790 & then 1791 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1792 1793 else 1794 ! Not found 1795 found=0 1796 quit=1 1797 exit 1798 end if 1799 1800 end if 1801 end do 1802 if(quit==1)exit 1803 end do 1804 end if 1805 1806 ! Now, if still found, put the correct value into array d2 1807 if(found==1)then 1808 1809 ! In case of phonons, need to take into account the 1810 ! time-reversal symmetry, and the shift back to the unit cell 1811 ! 1812 ! XG990712 : I am not sure this must be kept for electric field ... 1813 ! 1) Consider time-reversal symmetry 1814 sumi=isgn*sumi 1815 1816 if(ipert1<=natom .and. ipert2<=natom)then 1817 ! 2) Shift the atoms back to the unit cell. 1818 arg1=two_pi*( qpt(1)*indsym(1,isym,ipert1)& 1819 & +qpt(2)*indsym(2,isym,ipert1)& 1820 & +qpt(3)*indsym(3,isym,ipert1) ) 1821 arg2=two_pi*( qpt(1)*indsym(1,isym,ipert2)& 1822 & +qpt(2)*indsym(2,isym,ipert2)& 1823 & +qpt(3)*indsym(3,isym,ipert2) ) 1824 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 1825 ! XG010117 Must use isgn 1826 im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)) 1827 else 1828 re=one 1829 im=zero 1830 end if 1831 1832 ! Final check, could still fail if the 1833 ! element was its own symmetric 1834 if (abs(1.0_dp-re*noccur)< 1.0d-6.and.abs(im*noccur) <1.0d-6) then 1835 found=0 1836 end if 1837 1838 end if 1839 1840 if(found==1)then 1841 1842 if(noccur==0)then 1843 d2(1,idir1,ipert1,idir2,ipert2)=re*sumr-im*sumi 1844 d2(2,idir1,ipert1,idir2,ipert2)=re*sumi+im*sumr 1845 else 1846 ! See page July 2, 1994 in computer codes notebook 1847 xr=re*sumr-im*sumi 1848 xi=re*sumi+im*sumr 1849 norm=one+noccur**2-two*re*noccur 1850 xr=xr/norm 1851 xi=xi/norm 1852 d2(1,idir1,ipert1,idir2,ipert2)=& 1853 & (one-re*noccur)*xr-im*noccur*xi 1854 d2(2,idir1,ipert1,idir2,ipert2)=& 1855 & (one-re*noccur)*xi+im*noccur*xr 1856 end if 1857 1858 ! The element has been constructed ! 1859 blkflg(idir1,ipert1,idir2,ipert2)=1 1860 1861 quit1=1 1862 exit ! Exit loop on symmetry operations 1863 end if 1864 1865 ! End loop on all symmetries + time-reversal 1866 end if 1867 end do 1868 if(quit1==1)exit 1869 end do 1870 1871 end if 1872 end do ! End big loop on all elements 1873 end do 1874 end do 1875 end do 1876 1877 end do ! End Big Big Loop 1878 1879 !MT oct. 20, 2014: 1880 !Once the matrix has been built, it does not necessarily fulfill the correct symmetries. 1881 !It has just been filled up from rows or columns that only fulfill symmetries perserving 1882 !one particular perturbation. 1883 !An additional symmetrization might solve this (do not consider TR-symmetry) 1884 if (do_final_sym) then 1885 ABI_ALLOCATE(d2tmp1,(2,3,3)) 1886 ABI_ALLOCATE(d2tmp2,(2,3,3)) 1887 ABI_ALLOCATE(d2work,(2,3,mpert,3,mpert)) 1888 d2work(:,:,:,:,:)=d2(:,:,:,:,:) 1889 do ipert1=1,min(natom+2,mpert) 1890 if ((ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).or.(ipert1==natom+2.and.(.not.qzero))) cycle 1891 if (ipert1<=natom)then 1892 sym1_ => symrec 1893 else 1894 sym1_ => symrel 1895 end if 1896 do ipert2=1,min(natom+2,mpert) 1897 ! if (any(blkflg(:,ipert1,:,ipert2)==0)) cycle 1898 if ((ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11).or.(ipert2==natom+2.and.(.not.qzero))) cycle 1899 if (ipert2<=natom)then 1900 sym2_ => symrec 1901 else 1902 sym2_ => symrel 1903 end if 1904 nsym_used=0 1905 d2tmp2(:,:,:)=zero 1906 do isym=1,nsym 1907 if (symq(4,1,isym)==1) then 1908 ipesy1=ipert1;if (ipert1<=natom) ipesy1=indsym(4,isym,ipert1) 1909 ipesy2=ipert2;if (ipert2<=natom) ipesy2=indsym(4,isym,ipert2) 1910 if (all(blkflg(:,ipesy1,:,ipesy2)==1)) then 1911 nsym_used=nsym_used+1 1912 re=one;im=zero 1913 if (ipert1<=natom.and.ipert2<=natom.and.(.not.qzero)) then 1914 arg1=two_pi*(qpt(1)*(indsym(1,isym,ipert1)-indsym(1,isym,ipert2)) & 1915 & +qpt(2)*(indsym(2,isym,ipert1)-indsym(2,isym,ipert2)) & 1916 & +qpt(3)*(indsym(3,isym,ipert1)-indsym(3,isym,ipert2))) 1917 re=cos(arg1);im=sin(arg1) 1918 end if 1919 d2tmp1(:,:,:)=zero 1920 do idir2=1,3 !kappa 1921 do idir1=1,3 !mu 1922 do idisy1=1,3 !nu 1923 d2tmp1(:,idir1,idir2)=d2tmp1(:,idir1,idir2) & 1924 & +sym1_(idir1,idisy1,isym)*d2(:,idisy1,ipesy1,idir2,ipesy2) 1925 end do 1926 end do 1927 end do 1928 do idir2=1,3 !mu 1929 do idir1=1,3 !kappa 1930 do idisy2=1,3 !nu 1931 d2tmp2(1,idir1,idir2)=d2tmp2(1,idir1,idir2) & 1932 & +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*re-d2tmp1(2,idir1,idisy2)*im) 1933 d2tmp2(2,idir1,idir2)=d2tmp2(2,idir1,idir2) & 1934 & +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*im+d2tmp1(2,idir1,idisy2)*re) 1935 end do 1936 end do 1937 end do 1938 end if 1939 end if 1940 end do ! isym 1941 if (nsym_used>0) d2work(:,1:3,ipert1,1:3,ipert2)=d2tmp2(:,1:3,1:3)/dble(nsym_used) 1942 end do !ipert2 1943 end do !ipert1 1944 if (mpert>=natom) d2(:,1:3,1:natom,1:3,1:natom)=d2work(:,1:3,1:natom,1:3,1:natom) 1945 if (mpert>=natom+2) then 1946 d2(:,1:3,natom+2,1:3,1:natom)=d2work(:,1:3,natom+2,1:3,1:natom) 1947 d2(:,1:3,1:natom,1:3,natom+2)=d2work(:,1:3,1:natom,1:3,natom+2) 1948 d2(:,1:3,natom+2,1:3,natom+2)=d2work(:,1:3,natom+2,1:3,natom+2) 1949 end if 1950 ABI_DEALLOCATE(d2tmp1) 1951 ABI_DEALLOCATE(d2tmp2) 1952 ABI_DEALLOCATE(d2work) 1953 end if 1954 1955 end subroutine d2sym3
m_dynmat/d3sym [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
d3sym
FUNCTION
Given a set of calculated elements of the 3DTE matrix, build (nearly) all the other matrix elements that can be build using symmetries.
INPUTS
indsym(4,nsym,natom)=indirect indexing array : for each isym,iatom, fourth element is label of atom into which iatom is sent by INVERSE of symmetry operation isym; first three elements are the primitive translations which must be subtracted after the transformation to get back to the original unit cell. mpert =maximum number of ipert natom= number of atoms nsym=number of space group symmetries symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
SIDE EFFECTS
Input/Output blkflg(3,mpert,3,mpert,3,mpert)= matrix that indicates if an element of d3 is available (1 if available, 0 otherwise) d3(2,3,mpert,3,mpert,3,mpert)= matrix of the 3DTE
PARENTS
m_ddb,nonlinear
CHILDREN
SOURCE
4238 subroutine d3sym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel) 4239 4240 4241 !This section has been created automatically by the script Abilint (TD). 4242 !Do not modify the following lines by hand. 4243 #undef ABI_FUNC 4244 #define ABI_FUNC 'd3sym' 4245 !End of the abilint section 4246 4247 implicit none 4248 4249 !Arguments ------------------------------- 4250 !scalars 4251 integer,intent(in) :: mpert,natom,nsym 4252 !arrays 4253 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4254 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 4255 real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert) 4256 4257 !Local variables ------------------------- 4258 !scalars 4259 integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3 4260 integer :: ipesy1,ipesy2,ipesy3,isym,ithree 4261 real(dp) :: sumi,sumr 4262 !arrays 4263 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 4264 4265 ! ********************************************************************* 4266 4267 !DEBUG 4268 !write(std_out,*)'d3sym : enter' 4269 !do i1dir = 1, 3 4270 !do i2dir = 1, 3 4271 !do i3dir = 1, 3 4272 !write(std_out,*)i1dir,i2dir,i3dir,blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,natom+2) 4273 !end do 4274 !end do 4275 !end do 4276 !stop 4277 !ENDDEBUG 4278 4279 !First, take into account the permutations symmetry of 4280 !(i1pert,i1dir) and (i3pert,i3dir) 4281 4282 do i1pert = 1, mpert 4283 do i2pert = 1, mpert 4284 4285 do i3pert = 1, mpert 4286 4287 do i1dir = 1, 3 4288 do i2dir = 1, 3 4289 do i3dir = 1, 3 4290 4291 if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. & 4292 & (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=1)) then 4293 4294 d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = & 4295 & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4296 4297 blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = 1 4298 4299 end if 4300 4301 end do 4302 end do 4303 end do 4304 4305 end do 4306 end do 4307 end do 4308 4309 !Big Big Loop : symmetrize three times, because 4310 !of some cases in which one element is not yet available 4311 !at the first pass, and even at the second one ! 4312 4313 do ithree=1,3 4314 4315 ! Loop over perturbations 4316 do i1pert = 1, mpert 4317 do i2pert = 1, mpert 4318 do i3pert = 1, mpert 4319 4320 do i1dir = 1, 3 4321 do i2dir = 1, 3 4322 do i3dir = 1, 3 4323 4324 ! Will get element (idir1,ipert1,idir2,ipert2) 4325 ! so this element should not yet be present ... 4326 if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then 4327 4328 d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp 4329 4330 do isym = 1, nsym 4331 4332 found = 1 4333 4334 if (i1pert <= natom) then 4335 ipesy1 = indsym(4,isym,i1pert) 4336 sym1(:,:) = symrec(:,:,isym) 4337 else if (i1pert == natom + 2) then 4338 ipesy1 = i1pert 4339 sym1(:,:) = symrel(:,:,isym) 4340 else 4341 found = 0 4342 end if 4343 4344 if (i2pert <= natom) then 4345 ipesy2 = indsym(4,isym,i2pert) 4346 sym2(:,:) = symrec(:,:,isym) 4347 else if (i2pert == natom + 2) then 4348 ipesy2 = i2pert 4349 sym2(:,:) = symrel(:,:,isym) 4350 else 4351 found = 0 4352 end if 4353 4354 if (i3pert <= natom) then 4355 ipesy3 = indsym(4,isym,i3pert) 4356 sym3(:,:) = symrec(:,:,isym) 4357 else if (i3pert == natom + 2) then 4358 ipesy3 = i3pert 4359 sym3(:,:) = symrel(:,:,isym) 4360 else 4361 found = 0 4362 end if 4363 4364 sumr = 0_dp ; sumi = 0_dp; 4365 do idisy1 = 1, 3 4366 do idisy2 = 1, 3 4367 do idisy3 = 1, 3 4368 4369 if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) & 4370 & .and.(sym3(i3dir,idisy3) /=0)) then 4371 4372 if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then 4373 4374 sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4375 & sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4376 sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4377 & sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4378 4379 else 4380 4381 found = 0 4382 4383 end if 4384 4385 end if 4386 4387 end do 4388 end do 4389 end do 4390 4391 if (found == 1) then 4392 d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr 4393 d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 4394 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 4395 end if 4396 4397 end do ! isym 4398 4399 end if ! blkflg 4400 4401 ! Close loop over perturbations 4402 end do 4403 end do 4404 end do 4405 end do 4406 end do 4407 end do 4408 4409 end do ! close loop over ithree 4410 4411 end subroutine d3sym
m_dynmat/dfpt_phfrq [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dfpt_phfrq
FUNCTION
Get the phonon frequencies and eigenvectors (as well as the corresponding displacements) If q is at Gamma, the non-analytical behaviour can be included. Then, the effective dielectric tensor, the effective charges and oscillator strengths for the limiting direction are also returned
INPUTS
amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms) d2cart(2,3,mpert,3,mpert)=dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates indsym(4,msym*natom)=indirect indexing array : for each isym,iatom, fourth element is label of atom into which iatom is sent by INVERSE of symmetry operation isym; first three elements are the primitive translations which must be subtracted after the transformation to get back to the original unit cell. mpert =maximum number of ipert msym=maximum number of symmetries natom=number of atoms in unit cell nsym=number of space group symmetries ntypat=number of atom types qphnrm=(described above) qphon(3)= to be divided by qphnrm, give the phonon wavevector; if qphnrm==0.0_dp, then the wavevector is zero (Gamma point) and qphon gives the direction of the induced electric field in **CARTESIAN** coordinates. in the latter case, if qphon is zero, no non-analytical contribution is included. rprimd(3,3)=dimensional primitive translations (bohr) symdynmat=if 1, (re)symmetrize the dynamical matrix, except if Gamma wavevector with electric field added. symrel(3,3,nsym)=matrices of the group symmetries (real space) typat(natom)=integer label of each type of atom (1,2,...) ucvol=unit cell volume
OUTPUT
displ(2*3*natom*3*natom)= at the end, contains the displacements of atoms in cartesian coordinates.
The first index means either the real or the imaginary part,
The second index runs on the direction and the atoms displaced
The third index runs on the modes.
eigval(3*natom)=contains the eigenvalues of the dynamical matrix
eigvec(2*3*natom*3*natom)= at the end, contains the eigenvectors of the dynamical matrix in cartesian coordinates.
phfrq(3*natom)=phonon frequencies (square root of the dynamical matrix eigenvalues,
except if these are negative, and in this case, give minus the square root of the absolute value
of the matrix eigenvalues). Hartree units.
NOTES
1) One makes the dynamical matrix hermitian... 2) In case of q=Gamma, only the real part is used.
PARENTS
anaddb,m_ddb,m_effective_potential_file,m_ifc,m_phonons,respfn,thmeig
CHILDREN
SOURCE
5503 subroutine dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,& 5504 & mpert,msym,natom,nsym,ntypat,phfrq,qphnrm,qphon,rprimd,& 5505 & symdynmat,symrel,symafm,typat,ucvol) 5506 5507 5508 !This section has been created automatically by the script Abilint (TD). 5509 !Do not modify the following lines by hand. 5510 #undef ABI_FUNC 5511 #define ABI_FUNC 'dfpt_phfrq' 5512 !End of the abilint section 5513 5514 implicit none 5515 5516 !Arguments ------------------------------- 5517 !scalars 5518 integer,intent(in) :: mpert,msym,natom,nsym,ntypat,symdynmat 5519 real(dp),intent(in) :: qphnrm,ucvol 5520 !arrays 5521 integer,intent(in) :: indsym(4,msym*natom),symrel(3,3,nsym),typat(natom) 5522 integer,intent(in) :: symafm(nsym) 5523 real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),rprimd(3,3) 5524 real(dp),intent(inout) :: qphon(3) 5525 real(dp),intent(out) :: displ(2*3*natom*3*natom),eigval(3*natom) 5526 real(dp),intent(out) :: eigvec(2*3*natom*3*natom),phfrq(3*natom) 5527 5528 !Local variables ------------------------- 5529 !scalars 5530 integer :: analyt,i1,i2,idir1,idir2,ier,ii,imode,ipert1,ipert2 5531 integer :: jmode,indexi,indexj,index 5532 real(dp) :: epsq,norm,qphon2 5533 logical,parameter :: debug = .False. 5534 real(dp) :: sc_prod 5535 !arrays 5536 real(dp) :: qptn(3),dum(2,0) 5537 real(dp),allocatable :: matrx(:,:),zeff(:,:),zhpev1(:,:),zhpev2(:) 5538 5539 ! ********************************************************************* 5540 5541 !Prepare the diagonalisation: analytical part. 5542 !Note: displ is used as work space here 5543 i1=0 5544 do ipert1=1,natom 5545 do idir1=1,3 5546 i1=i1+1 5547 i2=0 5548 do ipert2=1,natom 5549 do idir2=1,3 5550 i2=i2+1 5551 index=i1+3*natom*(i2-1) 5552 displ(2*index-1)=d2cart(1,idir1,ipert1,idir2,ipert2) 5553 displ(2*index )=d2cart(2,idir1,ipert1,idir2,ipert2) 5554 end do 5555 end do 5556 end do 5557 end do 5558 5559 !Determine the analyticity of the matrix. 5560 analyt=1; if(abs(qphnrm)<tol8) analyt=0 5561 if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=2 5562 5563 !In case of q=Gamma, only the real part is used 5564 if(analyt==0 .or. analyt==2)then 5565 do i1=1,3*natom 5566 do i2=1,3*natom 5567 index=i1+3*natom*(i2-1) 5568 displ(2*index)=zero 5569 end do 5570 end do 5571 end if 5572 5573 !In the case the non-analyticity is required : 5574 ! MG: the tensor is in cartesian coordinates and this means that qphon must be in 5575 ! given in Cartesian coordinates. 5576 if(analyt==0)then 5577 5578 ! Normalize the limiting direction 5579 qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2 5580 qphon(:)=qphon(:)/sqrt(qphon2) 5581 5582 ! Get the dielectric constant for the limiting direction 5583 epsq=zero 5584 do idir1=1,3 5585 do idir2=1,3 5586 epsq= epsq + qphon(idir1)*qphon(idir2) * d2cart(1,idir1,natom+2,idir2,natom+2) 5587 end do 5588 end do 5589 5590 ABI_ALLOCATE(zeff,(3,natom)) 5591 5592 ! Get the effective charges for the limiting direction 5593 do idir1=1,3 5594 do ipert1=1,natom 5595 zeff(idir1,ipert1)=zero 5596 do idir2=1,3 5597 zeff(idir1,ipert1) = zeff(idir1,ipert1) + qphon(idir2)* d2cart(1,idir1,ipert1,idir2,natom+2) 5598 end do 5599 end do 5600 end do 5601 5602 ! Get the non-analytical part of the dynamical matrix, and suppress its imaginary part. 5603 i1=0 5604 do ipert1=1,natom 5605 do idir1=1,3 5606 i1=i1+1 5607 i2=0 5608 do ipert2=1,natom 5609 do idir2=1,3 5610 i2=i2+1 5611 index=i1+3*natom*(i2-1) 5612 displ(2*index-1)=displ(2*index-1)+four_pi/ucvol*zeff(idir1,ipert1)*zeff(idir2,ipert2)/epsq 5613 displ(2*index )=zero 5614 end do 5615 end do 5616 end do 5617 end do 5618 5619 ABI_DEALLOCATE(zeff) 5620 end if ! End of the non-analyticity treatment : 5621 5622 ! Multiply IFC(q) by masses 5623 call massmult_and_breaksym(natom, ntypat, typat, amu, displ) 5624 5625 !*********************************************************************** 5626 !Diagonalize the dynamical matrix 5627 5628 !Symmetrize the dynamical matrix 5629 !FIXME: swap the next 2 lines and update test files to include symmetrization for Gamma point too (except in non-analytic case) 5630 !if (symdynmat==1 .and. analyt > 0) then 5631 if (symdynmat==1 .and. analyt == 1) then 5632 qptn(:)=qphon(:) 5633 if (analyt==1) qptn(:)=qphon(:)/qphnrm 5634 call symdyma(displ,indsym,natom,nsym,qptn,rprimd,symrel,symafm) 5635 end if 5636 5637 ier=0; ii=1 5638 ABI_ALLOCATE(matrx,(2,(3*natom*(3*natom+1))/2)) 5639 do i2=1,3*natom 5640 do i1=1,i2 5641 matrx(1,ii)=displ(1+2*(i1-1)+2*(i2-1)*3*natom) 5642 matrx(2,ii)=displ(2+2*(i1-1)+2*(i2-1)*3*natom) 5643 ii=ii+1 5644 end do 5645 end do 5646 5647 ABI_ALLOCATE(zhpev1,(2,2*3*natom-1)) 5648 ABI_ALLOCATE(zhpev2,(3*3*natom-2)) 5649 5650 call ZHPEV ('V','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,zhpev2,ier) 5651 ABI_CHECK(ier == 0, sjoin('zhpev returned:', itoa(ier))) 5652 5653 ABI_DEALLOCATE(matrx) 5654 ABI_DEALLOCATE(zhpev1) 5655 ABI_DEALLOCATE(zhpev2) 5656 5657 if (debug) then 5658 ! Check the orthonormality of the eigenvectors 5659 do imode=1,3*natom 5660 do jmode=imode,3*natom 5661 indexi=2*3*natom*(imode-1) 5662 indexj=2*3*natom*(jmode-1) 5663 sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom)) 5664 write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod 5665 end do 5666 end do 5667 end if 5668 5669 !*********************************************************************** 5670 5671 !Get the phonon frequencies (negative by convention, if 5672 !the eigenvalue of the dynamical matrix is negative) 5673 do imode=1,3*natom 5674 if(eigval(imode)>=1.0d-16)then 5675 phfrq(imode)=sqrt(eigval(imode)) 5676 else if(eigval(imode)>=-1.0d-16)then 5677 phfrq(imode)=zero 5678 else 5679 phfrq(imode)=-sqrt(-eigval(imode)) 5680 end if 5681 end do 5682 5683 !Fix the phase of the eigenvectors 5684 call fxphas_seq(eigvec,dum,0,0,1,3*natom*3*natom,0,3*natom,3*natom,0) 5685 5686 !Normalise the eigenvectors 5687 do imode=1,3*natom 5688 norm=zero 5689 do idir1=1,3 5690 do ipert1=1,natom 5691 i1=idir1+(ipert1-1)*3 5692 index=i1+3*natom*(imode-1) 5693 norm=norm+eigvec(2*index-1)**2+eigvec(2*index)**2 5694 end do 5695 end do 5696 norm=sqrt(norm) 5697 do idir1=1,3 5698 do ipert1=1,natom 5699 i1=idir1+(ipert1-1)*3 5700 index=i1+3*natom*(imode-1) 5701 eigvec(2*index-1)=eigvec(2*index-1)/norm 5702 eigvec(2*index)=eigvec(2*index)/norm 5703 end do 5704 end do 5705 end do 5706 5707 !Get the phonon displacements 5708 do imode=1,3*natom 5709 do idir1=1,3 5710 do ipert1=1,natom 5711 i1=idir1+(ipert1-1)*3 5712 index=i1+3*natom*(imode-1) 5713 displ(2*index-1)=eigvec(2*index-1) / sqrt(amu(typat(ipert1))*amu_emass) 5714 displ(2*index )=eigvec(2*index ) / sqrt(amu(typat(ipert1))*amu_emass) 5715 end do 5716 end do 5717 end do 5718 5719 if (debug) then 5720 write(std_out,'(a)')' Phonon eigenvectors and displacements ' 5721 do imode=1,3*natom 5722 indexi=2*3*natom*(imode-1) 5723 write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' eigvec(1:6*natom)=',eigvec(indexi+1:indexi+6*natom) 5724 write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' displ(1:6*natom)=',displ(indexi+1:indexi+6*natom) 5725 end do 5726 5727 ! Check the orthonormality of the eigenvectors 5728 do imode=1,3*natom 5729 do jmode=imode,3*natom 5730 indexi=2*3*natom*(imode-1) 5731 indexj=2*3*natom*(jmode-1) 5732 sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom)) 5733 write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod 5734 end do 5735 end do 5736 end if 5737 5738 end subroutine dfpt_phfrq
m_dynmat/dist9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dist9
FUNCTION
Compute the distance between atoms
INPUTS
acell(3)=length scales by which rprim is to be multiplied dist(natom,natom,nrpt)=distances between atoms gprim(3,3)=dimensionless primitive translations in reciprocal space natom=number of atoms in unit cell nrpt= Number of R points in the Big Box rcan(3,natom)=canonical coordinates of atoms rprim(3,3)=dimensionless primitive translations in real space rpt(3,nrpt)=cartesian coordinates of the points in the BigBox.
OUTPUT
dist(natom,natom,nrpt)=distances between atoms
PARENTS
m_ifc
CHILDREN
SOURCE
3287 subroutine dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt) 3288 3289 3290 !This section has been created automatically by the script Abilint (TD). 3291 !Do not modify the following lines by hand. 3292 #undef ABI_FUNC 3293 #define ABI_FUNC 'dist9' 3294 !End of the abilint section 3295 3296 implicit none 3297 3298 !Arguments ------------------------------- 3299 !scalars 3300 integer,intent(in) :: natom,nrpt 3301 !arrays 3302 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3) 3303 real(dp),intent(in) :: rpt(3,nrpt) 3304 real(dp),intent(out) :: dist(natom,natom,nrpt) 3305 3306 !Local variables ------------------------- 3307 !scalars 3308 integer :: ia,ib,ii,irpt 3309 !arrays 3310 real(dp) :: ra(3),rb(3),rdiff(3),red(3),rptcar(3),xred(3) 3311 3312 ! ********************************************************************* 3313 3314 !BIG loop on all generic atoms 3315 do ia=1,natom 3316 3317 ! First transform canonical coordinates to reduced coordinates 3318 do ii=1,3 3319 xred(ii)=gprim(1,ii)*rcan(1,ia)+gprim(2,ii)*rcan(2,ia)& 3320 & +gprim(3,ii)*rcan(3,ia) 3321 end do 3322 ! Then to cartesian coordinates 3323 ra(:)=xred(1)*acell(1)*rprim(:,1)+xred(2)*acell(2)*rprim(:,2)+& 3324 & xred(3)*acell(3)*rprim(:,3) 3325 do ib=1,natom 3326 do ii=1,3 3327 xred(ii)=gprim(1,ii)*rcan(1,ib)+gprim(2,ii)*rcan(2,ib)& 3328 & +gprim(3,ii)*rcan(3,ib) 3329 end do 3330 do ii=1,3 3331 rb(ii)=xred(1)*acell(1)*rprim(ii,1)+& 3332 & xred(2)*acell(2)*rprim(ii,2)+& 3333 & xred(3)*acell(3)*rprim(ii,3) 3334 end do 3335 do irpt=1,nrpt 3336 ! First transform it to reduced coordinates 3337 do ii=1,3 3338 red(ii)=gprim(1,ii)*rpt(1,irpt)+gprim(2,ii)*rpt(2,irpt)& 3339 & +gprim(3,ii)*rpt(3,irpt) 3340 end do 3341 ! Then to cartesian coordinates 3342 do ii=1,3 3343 rptcar(ii)=red(1)*acell(1)*rprim(ii,1)+& 3344 & red(2)*acell(2)*rprim(ii,2)+& 3345 & red(3)*acell(3)*rprim(ii,3) 3346 end do 3347 do ii=1,3 3348 rdiff(ii)=-rptcar(ii)+ra(ii)-rb(ii) 3349 end do 3350 dist(ia,ib,irpt)=(rdiff(1)**2+rdiff(2)**2+rdiff(3)**2)**0.5 3351 3352 end do 3353 3354 end do 3355 end do 3356 3357 end subroutine dist9
m_dynmat/dymfz9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dymfz9
FUNCTION
As the subroutine canatm has transformed the coordinates of the atoms in normalized canonical coordinates, the corresponding dynamical matrix should be multiplied by a phase shift corresponding to the translation between New and Old coordinates of its two corresponding atoms.
INPUTS
dynmat = non-phase shifted dynamical matrices natom = number of atoms nqpt = number of qpoints gprim = reciprocal lattice vectors (cartesian but dimensionless) option=1 : the matrices are transformed from the old (tn) coordinate system to the new (normalized canonical) 2 : the matrices are restored from the normalized canonical coordinate system to the usual (tn) one... rcan = canonical coordinates of atoms spqpt = qpoint coordinates (reduced reciprocal) trans = Atomic translations : xred = rcan + trans
OUTPUT
dynmat = phase shifted dynamical matrices
PARENTS
m_dynmat,m_ifc
CHILDREN
SOURCE
5146 subroutine dymfz9(dynmat,natom,nqpt,gprim,option,spqpt,trans) 5147 5148 5149 !This section has been created automatically by the script Abilint (TD). 5150 !Do not modify the following lines by hand. 5151 #undef ABI_FUNC 5152 #define ABI_FUNC 'dymfz9' 5153 !End of the abilint section 5154 5155 implicit none 5156 5157 !Arguments ------------------------------- 5158 !scalars 5159 integer,intent(in) :: natom,nqpt,option 5160 !arrays 5161 real(dp),intent(in) :: gprim(3,3),spqpt(3,nqpt),trans(3,natom) 5162 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt) 5163 5164 !Local variables ------------------------- 5165 !scalars 5166 integer :: ia,ib,iqpt,mu,nu 5167 real(dp) :: im,ktrans,re 5168 !arrays 5169 real(dp) :: kk(3) 5170 5171 ! ********************************************************************* 5172 5173 DBG_ENTER("COLL") 5174 5175 do iqpt=1,nqpt 5176 ! Definition of q in normalized reciprocal space 5177 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 5178 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 5179 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 5180 5181 if(option==1)then 5182 kk(1)=-kk(1) 5183 kk(2)=-kk(2) 5184 kk(3)=-kk(3) 5185 end if 5186 5187 do ia=1,natom 5188 do ib=1,natom 5189 ! Product of q with the differences between the two atomic translations 5190 ktrans=kk(1)*(trans(1,ia)-trans(1,ib))+kk(2)*(trans(2,ia)-& 5191 & trans(2,ib))+kk(3)*(trans(3,ia)-trans(3,ib)) 5192 do mu=1,3 5193 do nu=1,3 5194 re=dynmat(1,mu,ia,nu,ib,iqpt) 5195 im=dynmat(2,mu,ia,nu,ib,iqpt) 5196 ! Transformation of the Old dynamical matrices by New ones by multiplication by a phase shift 5197 dynmat(1,mu,ia,nu,ib,iqpt)=re*cos(two_pi*ktrans)-im*sin(two_pi*ktrans) 5198 dynmat(2,mu,ia,nu,ib,iqpt)=re*sin(two_pi*ktrans)+im*cos(two_pi*ktrans) 5199 end do 5200 end do 5201 end do 5202 end do 5203 end do 5204 5205 DBG_EXIT("COLL") 5206 5207 end subroutine dymfz9
m_dynmat/dynmat_dq [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dynmat_dq
FUNCTION
Compute the derivative D(q)/dq of the dynamical matrix via Fourier transform of the interatomic forces
INPUTS
qpt(3)= Reduced coordinates of the q vector in reciprocal space natom= Number of atoms in the unit cell gprim(3,3)= Normalized coordinates in reciprocal space nrpt= Number of R points in the Big Box rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
OUTPUT
dddq(2,3,natom,3,natom,3)= Derivate of the dynamical matrix in cartesian coordinates. The tree directions are stored in the last dimension. These coordinates are normalized (=> * acell(3)!!)
PARENTS
m_ifc
CHILDREN
SOURCE
3610 subroutine dynmat_dq(qpt,natom,gprim,nrpt,rpt,atmfrc,wghatm,dddq) 3611 3612 3613 !This section has been created automatically by the script Abilint (TD). 3614 !Do not modify the following lines by hand. 3615 #undef ABI_FUNC 3616 #define ABI_FUNC 'dynmat_dq' 3617 !End of the abilint section 3618 3619 implicit none 3620 3621 !Arguments ------------------------------- 3622 !scalars 3623 integer,intent(in) :: natom,nrpt 3624 !arrays 3625 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt(3) 3626 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 3627 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 3628 real(dp),intent(out) :: dddq(2,3,natom,3,natom,3) 3629 3630 !Local variables ------------------------- 3631 !scalars 3632 integer :: ia,ib,irpt,mu,nu,ii 3633 real(dp) :: im,kr,re 3634 !arrays 3635 real(dp) :: kk(3),fact(2,3) 3636 3637 ! ********************************************************************* 3638 3639 dddq = zero 3640 3641 do irpt=1,nrpt 3642 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3643 kk(1) = qpt(1)*gprim(1,1)+ qpt(2)*gprim(1,2) + qpt(3)*gprim(1,3) 3644 kk(2) = qpt(1)*gprim(2,1)+ qpt(2)*gprim(2,2) + qpt(3)*gprim(2,3) 3645 kk(3) = qpt(1)*gprim(3,1)+ qpt(2)*gprim(3,2) + qpt(3)*gprim(3,3) 3646 3647 ! Product of k and r 3648 kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3649 3650 ! Get phase factor 3651 re=cos(two_pi*kr); im=sin(two_pi*kr) 3652 3653 ! Inner loop on atoms and directions 3654 do ib=1,natom 3655 do ia=1,natom 3656 if (abs(wghatm(ia,ib,irpt))>1.0d-10) then 3657 ! take into account rotation due to i. 3658 fact(1,:) = -im * wghatm(ia,ib,irpt) * rpt(:,irpt) 3659 fact(2,:) = re * wghatm(ia,ib,irpt) * rpt(:,irpt) 3660 do nu=1,3 3661 do mu=1,3 3662 ! Real and imaginary part of the dynamical matrices 3663 ! Atmfrc should be real 3664 do ii=1,3 3665 dddq(1,mu,ia,nu,ib,ii) = dddq(1,mu,ia,nu,ib,ii) + fact(1,ii) * atmfrc(mu,ia,nu,ib,irpt) 3666 dddq(2,mu,ia,nu,ib,ii) = dddq(2,mu,ia,nu,ib,ii) + fact(2,ii) * atmfrc(mu,ia,nu,ib,irpt) 3667 end do 3668 end do 3669 end do 3670 end if 3671 end do 3672 end do 3673 end do 3674 3675 end subroutine dynmat_dq
m_dynmat/ftgam [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ftgam
FUNCTION
If qtor=1 (q->r): Generates the Fourier transform of the recip space gkk matrices to obtain the real space ones. If qtor=0 (r->q): Generates the Fourier transform of the real space gkk matrices to obtain the reciprocal space ones.
INPUTS
natom= Number of atoms in the unit cell
nqpt= Number of q points in the Brillouin zone
if qtor=0 this number is read in the input file
nrpt= Number of R points in the Big Box
qtor= ( q to r : see above )
rpt(3,nprt)= Canonical coordinates of the R points in the unit cell
These coordinates are normalized (=> * acell(3)!!)
qpt_full(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
if qtor=0 these vectors are read in the input file
wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
OUTPUT
(see side effects)
SIDE EFFECTS
Input/output gam_qpt(2,3*natom*3*natom,nqpt) = gamma matrices in recip space coming from the Derivative Data Base gam_rpt(2,3*natom*3*natom,nrpt) = gamma matrices in real space stored in file unit_gkk_rpt
PARENTS
elphon,get_tau_k,integrate_gamma_alt,m_phgamma,mka2f,mka2f_tr mka2f_tr_lova,mkph_linwid
CHILDREN
NOTES
copied from ftiaf9.f recip to real space: real space is forced to disk file unit_gkk_rpt recip space depends on gkqwrite and unitgkq3 real to recip space: real space is forced to disk file unit_gkk_rpt recip space is necessarily in memory in gkk_qpt real space elements are complex, but could be reduced, as (-r) = (+r)*
SOURCE
5880 subroutine ftgam (wghatm,gam_qpt,gam_rpt,natom,nqpt,nrpt,qtor,coskr, sinkr) 5881 5882 5883 !This section has been created automatically by the script Abilint (TD). 5884 !Do not modify the following lines by hand. 5885 #undef ABI_FUNC 5886 #define ABI_FUNC 'ftgam' 5887 !End of the abilint section 5888 5889 implicit none 5890 5891 !Arguments ------------------------------- 5892 !scalars 5893 integer,intent(in) :: natom,nqpt,nrpt,qtor 5894 !arrays 5895 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 5896 real(dp),intent(inout) :: gam_qpt(2,3*natom*3*natom,nqpt) 5897 real(dp),intent(inout) :: gam_rpt(2,3*natom*3*natom,nrpt) 5898 real(dp),intent(in) :: coskr(nqpt,nrpt) 5899 real(dp),intent(in) :: sinkr(nqpt,nrpt) 5900 5901 !Local variables ------------------------- 5902 !scalars 5903 integer :: iatom,idir,ip,iqpt,irpt,jatom,jdir 5904 real(dp) :: im,re 5905 character(len=500) :: message 5906 5907 ! ********************************************************************* 5908 5909 select case (qtor) 5910 ! 5911 case (1) !Recip to real space 5912 gam_rpt(:,:,:) = zero 5913 do irpt=1,nrpt 5914 do iqpt=1,nqpt 5915 ! Get the phase factor with normalization! 5916 re=coskr(iqpt,irpt) 5917 im=sinkr(iqpt,irpt) 5918 do ip=1,3*natom*3*natom 5919 ! Real and imaginary part of the real-space gam matrices 5920 gam_rpt(1,ip,irpt) = gam_rpt(1,ip,irpt)& 5921 & +re*gam_qpt(1,ip,iqpt) & 5922 & +im*gam_qpt(2,ip,iqpt) 5923 gam_rpt(2,ip,irpt) = gam_rpt(2,ip,irpt)& 5924 & +re*gam_qpt(2,ip,iqpt) & 5925 & -im*gam_qpt(1,ip,iqpt) 5926 end do 5927 end do 5928 end do 5929 gam_rpt = gam_rpt/nqpt 5930 ! 5931 case (0) ! Recip space from real space 5932 5933 gam_qpt(:,:,:)=zero 5934 5935 do irpt=1,nrpt 5936 do iqpt=1,nqpt 5937 5938 do iatom=1,natom 5939 do jatom=1,natom 5940 re = coskr(iqpt,irpt)*wghatm(iatom,jatom,irpt) 5941 im = sinkr(iqpt,irpt)*wghatm(iatom,jatom,irpt) 5942 5943 do idir=1,3 5944 do jdir=1,3 5945 ! Get phase factor 5946 5947 ip= jdir + (jatom-1)*3 + (idir-1)*3*natom + (iatom-1)*9*natom 5948 ! Real and imaginary part of the interatomic forces 5949 gam_qpt(1,ip,iqpt)=& 5950 & gam_qpt(1,ip,iqpt)& 5951 & +re*gam_rpt(1,ip,irpt)& 5952 & -im*gam_rpt(2,ip,irpt) 5953 ! !DEBUG 5954 gam_qpt(2,ip,iqpt)=& 5955 & gam_qpt(2,ip,iqpt)& 5956 & +im*gam_rpt(1,ip,irpt)& 5957 & +re*gam_rpt(2,ip,irpt) 5958 ! !ENDDEBUG 5959 5960 end do ! end jdir 5961 end do ! end idir 5962 end do 5963 end do ! end iatom 5964 5965 end do ! end iqpt 5966 end do ! end irpt 5967 5968 case default ! There is no other space to Fourier transform from 5969 write(message,'(a,i0,a)' )' The only allowed values for qtor are 0 or 1, while qtor=',qtor,' has been required.' 5970 MSG_BUG(message) 5971 end select 5972 5973 end subroutine ftgam
m_dynmat/ftgam_init [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ftgam_init
FUNCTION
Generates the sin and cos phases for the Fourier transform of the gkk matrices
INPUTS
gprim = reciprocal space vectors to get cartesian coord for qpt nqpt= Number of q points in the Brillouin zone nrpt= Number of R points in the Big Box rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) qpt_full(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space if qtor=0 these vectors are read in the input file
OUTPUT
coskr, sinkr = cosine and sine of phase factors for given r and q points
PARENTS
elphon,get_tau_k,integrate_gamma_alt,m_phgamma,mka2f,mka2f_tr mka2f_tr_lova,mkph_linwid
CHILDREN
SOURCE
6004 subroutine ftgam_init (gprim,nqpt,nrpt,qpt_full,rpt,coskr, sinkr) 6005 6006 6007 !This section has been created automatically by the script Abilint (TD). 6008 !Do not modify the following lines by hand. 6009 #undef ABI_FUNC 6010 #define ABI_FUNC 'ftgam_init' 6011 !End of the abilint section 6012 6013 implicit none 6014 6015 !Arguments ------------------------------- 6016 !scalars 6017 integer,intent(in) :: nqpt,nrpt 6018 !arrays 6019 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt_full(3,nqpt) 6020 real(dp),intent(out) :: coskr(nqpt,nrpt) 6021 real(dp),intent(out) :: sinkr(nqpt,nrpt) 6022 6023 !Local variables ------------------------- 6024 !scalars 6025 integer :: iqpt,irpt 6026 real(dp) :: kr 6027 !arrays 6028 real(dp) :: kk(3) 6029 6030 ! ********************************************************************* 6031 6032 ! Prepare the phase factors 6033 do iqpt=1,nqpt 6034 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 6035 kk(1)= qpt_full(1,iqpt)*gprim(1,1)+& 6036 & qpt_full(2,iqpt)*gprim(1,2)+& 6037 & qpt_full(3,iqpt)*gprim(1,3) 6038 6039 kk(2)= qpt_full(1,iqpt)*gprim(2,1)+& 6040 & qpt_full(2,iqpt)*gprim(2,2)+& 6041 & qpt_full(3,iqpt)*gprim(2,3) 6042 6043 kk(3)= qpt_full(1,iqpt)*gprim(3,1)+& 6044 & qpt_full(2,iqpt)*gprim(3,2)+& 6045 & qpt_full(3,iqpt)*gprim(3,3) 6046 do irpt=1,nrpt 6047 ! Product of k and r 6048 kr = kk(1)*rpt(1,irpt)+ kk(2)*rpt(2,irpt)+ kk(3)*rpt(3,irpt) 6049 coskr(iqpt,irpt)=cos(two_pi*kr) 6050 sinkr(iqpt,irpt)=sin(two_pi*kr) 6051 end do 6052 end do 6053 6054 end subroutine ftgam_init
m_dynmat/ftifc_q2r [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ftifc_q2r
FUNCTION
Generates the Fourier transform of the dynamical matrices to obtain interatomic forces (real space).
INPUTS
dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base gprim(3,3)= Normalized coordinates in reciprocal space natom= Number of atoms in the unit cell nqpt= Number of q points in the Brillouin zone nrpt= Number of R points in the Big Box rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) spqpt(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space
OUTPUT
atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space !! We used the imaginary part just for debugging !
PARENTS
m_ifc
CHILDREN
SOURCE
3392 subroutine ftifc_q2r(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt) 3393 3394 3395 !This section has been created automatically by the script Abilint (TD). 3396 !Do not modify the following lines by hand. 3397 #undef ABI_FUNC 3398 #define ABI_FUNC 'ftifc_q2r' 3399 !End of the abilint section 3400 3401 implicit none 3402 3403 !Arguments ------------------------------- 3404 !scalars 3405 integer,intent(in) :: natom,nqpt,nrpt 3406 !arrays 3407 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt) 3408 real(dp),intent(out) :: atmfrc(3,natom,3,natom,nrpt) 3409 real(dp),intent(in) :: dynmat(2,3,natom,3,natom,nqpt) 3410 3411 !Local variables ------------------------- 3412 !scalars 3413 integer :: ia,ib,iqpt,irpt,mu,nu 3414 real(dp) :: im,kr,re 3415 !arrays 3416 real(dp) :: kk(3) 3417 3418 ! ********************************************************************* 3419 3420 DBG_ENTER("COLL") 3421 3422 !Interatomic Forces from Dynamical Matrices 3423 atmfrc = zero 3424 do irpt=1,nrpt 3425 do iqpt=1,nqpt 3426 3427 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3428 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 3429 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 3430 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 3431 3432 ! Product of k and r 3433 kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3434 3435 ! Get the phase factor 3436 re=cos(two_pi*kr) 3437 im=sin(two_pi*kr) 3438 3439 ! Now, big inner loops on atoms and directions 3440 3441 ! The indices are ordered to give better speed 3442 do ib=1,natom 3443 do nu=1,3 3444 do ia=1,natom 3445 do mu=1,3 3446 ! Real and imaginary part of the interatomic forces 3447 atmfrc(mu,ia,nu,ib,irpt)=atmfrc(mu,ia,nu,ib,irpt)& 3448 & +re*dynmat(1,mu,ia,nu,ib,iqpt)& 3449 & +im*dynmat(2,mu,ia,nu,ib,iqpt) 3450 ! The imaginary part should be equal to zero !!!!!! 3451 ! atmfrc(2,mu,ia,nu,ib,irpt)=atmfrc(2,mu,ia,nu,ib,irpt) 3452 ! & +re*dynmat(2,mu,ia,nu,ib,iqpt) 3453 ! & -im*dynmat(1,mu,ia,nu,ib,iqpt) 3454 end do 3455 end do 3456 end do 3457 end do 3458 3459 end do 3460 end do 3461 3462 !The sumifc has to be weighted by a normalization factor of 1/nqpt 3463 atmfrc = atmfrc/nqpt 3464 3465 DBG_EXIT("COLL") 3466 3467 end subroutine ftifc_q2r
m_dynmat/ftifc_r2q [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ftifc_r2q
FUNCTION
(r->q): Generates the Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space).
INPUTS
atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space gprim(3,3)= Normalized coordinates in reciprocal space natom= Number of atoms in the unit cell nqpt= Number of q points in the Brillouin zone nrpt= Number of R points in the Big Box rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) spqpt(3,nqpt)= Reduced coordinates of the q vectors in reciprocal space wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
OUTPUT
dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
PARENTS
m_dynmat
CHILDREN
SOURCE
3503 subroutine ftifc_r2q(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt,wghatm) 3504 3505 3506 !This section has been created automatically by the script Abilint (TD). 3507 !Do not modify the following lines by hand. 3508 #undef ABI_FUNC 3509 #define ABI_FUNC 'ftifc_r2q' 3510 !End of the abilint section 3511 3512 implicit none 3513 3514 !Arguments ------------------------------- 3515 !scalars 3516 integer,intent(in) :: natom,nqpt,nrpt 3517 !arrays 3518 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt) 3519 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 3520 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 3521 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt) 3522 3523 !Local variables ------------------------- 3524 !scalars 3525 integer :: ia,ib,iqpt,irpt,mu,nu 3526 real(dp) :: facti,factr,im,kr,re 3527 !arrays 3528 real(dp) :: kk(3) 3529 3530 ! ********************************************************************* 3531 3532 dynmat = zero 3533 3534 do iqpt=1,nqpt 3535 do irpt=1,nrpt 3536 3537 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3538 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)* gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 3539 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)* gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 3540 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)* gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 3541 3542 ! Product of k and r 3543 kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3544 3545 ! Get phase factor 3546 re=cos(two_pi*kr) 3547 im=sin(two_pi*kr) 3548 3549 ! Inner loop on atoms and directions 3550 do ib=1,natom 3551 do ia=1,natom 3552 if(abs(wghatm(ia,ib,irpt))>1.0d-10)then 3553 factr=re*wghatm(ia,ib,irpt) 3554 facti=im*wghatm(ia,ib,irpt) 3555 do nu=1,3 3556 do mu=1,3 3557 ! Real and imaginary part of the dynamical matrices 3558 dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt)& 3559 & +factr*atmfrc(mu,ia,nu,ib,irpt) 3560 ! Atmfrc should be real 3561 ! & -im*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt) 3562 dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt)& 3563 & +facti*atmfrc(mu,ia,nu,ib,irpt) 3564 ! Atmfrc should be real 3565 ! & +re*wghatm(ia,ib,irpt)*atmfrc(2,mu,ia,nu,ib,irpt) 3566 end do 3567 end do 3568 end if 3569 end do 3570 end do 3571 end do 3572 end do 3573 3574 end subroutine ftifc_r2q
m_dynmat/gtdyn9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
gtdyn9
FUNCTION
Generates a dynamical matrix from interatomic force constants and long-range electrostatic interactions.
INPUTS
acell(3)=length scales by which rprim is to be multiplied atmfrc(3,natom,3,natom,nrpt) = Interatomic Forces in real space (imaginary part only for debugging) dielt(3,3) = dielectric tensor dipdip= if 0, no dipole-dipole interaction was subtracted in atmfrc if 1, atmfrc has been build without dipole-dipole part dyewq0(3,3,natom)= Ewald part of the dynamical matrix, at q=0. gmet(3,3)= metric tensor in reciprocal space. gprim(3,3)= Normalized coordinates in reciprocal space mpert =maximum number of ipert natom= Number of atoms in the unit cell nrpt= Number of R points in the Big Box qphnrm= Normalisation coefficient for qpt qpt(3)= Reduced coordinates of the q vectors in reciprocal space rmet(3,3)= Metric tensor in real space. rprim(3,3)= dimensionless primitive translations in real space rpt(3,nprt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) trans(3,natom)= Atomic translations : xred = rcan + trans ucvol= unit cell volume wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector xred(3,natom)= relative coords of atoms in unit cell (dimensionless) zeff(3,3,natom)=effective charge on each atom, versus electric field and atomic displacement
OUTPUT
d2cart(2,3,mpert,3,mpert)=dynamical matrix obtained for the wavevector qpt (normalized using qphnrm)
PARENTS
anaddb,ddb_interpolate,m_effective_potential_file,m_ifc,m_phonons
CHILDREN
SOURCE
5349 subroutine gtdyn9(acell,atmfrc,dielt,dipdip,& 5350 & dyewq0,d2cart,gmet,gprim,mpert,natom,& 5351 & nrpt,qphnrm,qpt,rmet,rprim,rpt,& 5352 & trans,ucvol,wghatm,xred,zeff) 5353 5354 5355 !This section has been created automatically by the script Abilint (TD). 5356 !Do not modify the following lines by hand. 5357 #undef ABI_FUNC 5358 #define ABI_FUNC 'gtdyn9' 5359 !End of the abilint section 5360 5361 implicit none 5362 5363 !Arguments ------------------------------- 5364 !scalars 5365 integer,intent(in) :: dipdip,mpert,natom,nrpt 5366 real(dp),intent(in) :: qphnrm,ucvol 5367 !arrays 5368 real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qpt(3) 5369 real(dp),intent(in) :: rmet(3,3),rprim(3,3),rpt(3,nrpt) 5370 real(dp),intent(in) :: trans(3,natom),wghatm(natom,natom,nrpt),xred(3,natom) 5371 real(dp),intent(in) :: zeff(3,3,natom) 5372 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 5373 real(dp),intent(in) :: dyewq0(3,3,natom) 5374 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 5375 5376 !Local variables ------------------------- 5377 !scalars 5378 integer,parameter :: nqpt1=1,option2=2,sumg0=0,plus1=1,iqpt1=1 5379 integer :: i1,i2,ib,nsize 5380 !arrays 5381 real(dp) :: qphon(3) 5382 real(dp),allocatable :: dq(:,:,:,:,:),dyew(:,:,:,:,:) 5383 5384 ! ********************************************************************* 5385 5386 ABI_ALLOCATE(dq,(2,3,natom,3,natom)) 5387 5388 !Get the normalized wavevector 5389 if(abs(qphnrm)<1.0d-7)then 5390 qphon(1:3)=zero 5391 else 5392 qphon(1:3)=qpt(1:3)/qphnrm 5393 end if 5394 5395 !Generate the analytical part from the interatomic forces 5396 call ftifc_r2q (atmfrc,dq,gprim,natom,nqpt1,nrpt,rpt,qphon,wghatm) 5397 5398 !The analytical dynamical matrix dq has been generated 5399 !in the normalized canonical coordinate system. Now, the 5400 !phase is modified, in order to recover the usual (xred) coordinate of atoms. 5401 call dymfz9(dq,natom,nqpt1,gprim,option2,qphon,trans) 5402 5403 if (dipdip==1) then 5404 ! Add the non-analytical part 5405 ! Compute dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix, 5406 ! second energy derivative wrt xred(3,natom) in Hartrees 5407 ! (Denoted A-bar in the notes) 5408 ABI_ALLOCATE(dyew,(2,3,natom,3,natom)) 5409 5410 call ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff) 5411 call q0dy3_apply(natom,dyewq0,dyew) 5412 call nanal9(dyew,dq,iqpt1,natom,nqpt1,plus1) 5413 5414 ABI_DEALLOCATE(dyew) 5415 end if 5416 5417 !Copy the dynamical matrix in the proper location 5418 5419 !First zero all the elements 5420 nsize=2*(3*mpert)**2 5421 d2cart = zero 5422 5423 !Copy the elements from dq to d2cart 5424 d2cart(:,:,1:natom,:,1:natom)=dq(:,:,1:natom,:,1:natom) 5425 5426 !In case we have the gamma point, 5427 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)then 5428 ! Copy the effective charge and dielectric constant in the final array 5429 do i1=1,3 5430 do i2=1,3 5431 d2cart(1,i1,natom+2,i2,natom+2)=dielt(i1,i2) 5432 do ib=1,natom 5433 d2cart(1,i1,natom+2,i2,ib)=zeff(i1,i2,ib) 5434 d2cart(1,i2,ib,i1,natom+2)=zeff(i1,i2,ib) 5435 end do 5436 end do 5437 end do 5438 end if 5439 5440 ABI_DEALLOCATE(dq) 5441 5442 end subroutine gtdyn9
m_dynmat/ifclo9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ifclo9
FUNCTION
Convert from cartesian coordinates to local coordinates the 3*3 interatomic force constant matrix
INPUTS
ifccar(3,3)= matrix of interatomic force constants in cartesian coordinates vect1(3)= cartesian coordinates of the first local vector vect2(3)= cartesian coordinates of the second local vector vect3(3)= cartesian coordinates of the third local vector
OUTPUT
ifcloc(3,3)= matrix of interatomic force constants in local coordinates
PARENTS
m_ifc
CHILDREN
SOURCE
3705 subroutine ifclo9(ifccar,ifcloc,vect1,vect2,vect3) 3706 3707 3708 !This section has been created automatically by the script Abilint (TD). 3709 !Do not modify the following lines by hand. 3710 #undef ABI_FUNC 3711 #define ABI_FUNC 'ifclo9' 3712 !End of the abilint section 3713 3714 implicit none 3715 3716 !Arguments ------------------------------- 3717 !arrays 3718 real(dp),intent(in) :: ifccar(3,3),vect1(3),vect2(3),vect3(3) 3719 real(dp),intent(out) :: ifcloc(3,3) 3720 3721 !Local variables ------------------------- 3722 !scalars 3723 integer :: ii,jj 3724 !arrays 3725 real(dp) :: work(3,3) 3726 3727 ! ********************************************************************* 3728 3729 do jj=1,3 3730 do ii=1,3 3731 work(jj,ii)=zero 3732 end do 3733 do ii=1,3 3734 work(jj,1)=work(jj,1)+ifccar(jj,ii)*vect1(ii) 3735 work(jj,2)=work(jj,2)+ifccar(jj,ii)*vect2(ii) 3736 work(jj,3)=work(jj,3)+ifccar(jj,ii)*vect3(ii) 3737 end do 3738 end do 3739 3740 do jj=1,3 3741 do ii=1,3 3742 ifcloc(ii,jj)=zero 3743 end do 3744 do ii=1,3 3745 ifcloc(1,jj)=ifcloc(1,jj)+vect1(ii)*work(ii,jj) 3746 ifcloc(2,jj)=ifcloc(2,jj)+vect2(ii)*work(ii,jj) 3747 ifcloc(3,jj)=ifcloc(3,jj)+vect3(ii)*work(ii,jj) 3748 end do 3749 end do 3750 3751 end subroutine ifclo9
m_dynmat/make_bigbox [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
make_bigbox
FUNCTION
Helper functions that faciliates the generation of a Big Box containing all the R points in the cartesian real space needed to Fourier Transform the dynamical matrix into its corresponding interatomic force. See bigbx9 for the algorithm.
INPUTS
brav= Bravais Lattice (1 or -1=S.C.;2=F.C.C.;3=BCC;4=Hex.) ngqpt(3)= Numbers used to generate the q points to sample the Brillouin zone using an homogeneous grid nqshft= number of q-points in the repeated cell for the Brillouin zone sampling When nqshft is not 1, but 2 or 4 (only other allowed values), the limits for the big box have to be extended by a factor of 2. rprim(3,3)= Normalized coordinates in real space !!! IS THIS CORRECT?
OUTPUT
cell= (nrpt,3) Give the index of the the cell and irpt nprt= Number of R points in the Big Box rpt(3,nrpt)= Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)!!) The array is allocated here with the proper dimension. Client code is responsible for the deallocation.
PARENTS
m_ifc
CHILDREN
SOURCE
2579 subroutine make_bigbox(brav,cell,ngqpt,nqshft,rprim,nrpt,rpt) 2580 2581 2582 !This section has been created automatically by the script Abilint (TD). 2583 !Do not modify the following lines by hand. 2584 #undef ABI_FUNC 2585 #define ABI_FUNC 'make_bigbox' 2586 !End of the abilint section 2587 2588 implicit none 2589 2590 !Arguments ------------------------------- 2591 !scalars 2592 integer,intent(in) :: brav,nqshft 2593 integer,intent(out) :: nrpt 2594 !arrays 2595 integer,intent(in) :: ngqpt(3) 2596 real(dp),intent(in) :: rprim(3,3) 2597 real(dp),allocatable,intent(out) :: rpt(:,:) 2598 integer,allocatable,intent(out) :: cell(:,:) 2599 2600 !Local variables ------------------------- 2601 !scalars 2602 integer :: choice,mrpt 2603 !arrays 2604 real(dp) :: dummy_rpt(3,1) 2605 integer:: dummy_cell(1,3) 2606 2607 ! ********************************************************************* 2608 2609 ! Compute the number of points (cells) in real space 2610 choice=0 2611 call bigbx9(brav,dummy_cell,choice,1,ngqpt,nqshft,mrpt,rprim,dummy_rpt) 2612 2613 ! Now we can allocate and calculate the points and the weights. 2614 nrpt = mrpt 2615 ABI_ALLOCATE(rpt,(3,nrpt)) 2616 ABI_ALLOCATE(cell,(3,nrpt)) 2617 2618 choice=1 2619 call bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt) 2620 2621 end subroutine make_bigbox
m_dynmat/massmult_and_breaksym [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
mult_masses_and_break_symms
FUNCTION
Multiply the IFC(q) by the atomic masses, slightly break symmetry to make tests more portable and make the matrix hermitian before returning.
INPUTS
amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms) natom=number of atoms in unit cell ntypat=number of atom types typat(natom)=integer label of each type of atom (1,2,...)
SIDE EFFECTS
mat(2*3*natom*3*natom)=Multiplies by atomic masses in output.
PARENTS
CHILDREN
SOURCE
5765 subroutine massmult_and_breaksym(natom, ntypat, typat, amu, mat) 5766 5767 5768 !This section has been created automatically by the script Abilint (TD). 5769 !Do not modify the following lines by hand. 5770 #undef ABI_FUNC 5771 #define ABI_FUNC 'massmult_and_breaksym' 5772 !End of the abilint section 5773 5774 implicit none 5775 5776 !Arguments ------------------------------- 5777 !scalars 5778 integer,intent(in) :: natom,ntypat 5779 !arrays 5780 integer,intent(in) :: typat(natom) 5781 real(dp),intent(in) :: amu(ntypat) 5782 real(dp),intent(inout) :: mat(2*3*natom*3*natom) 5783 5784 !Local variables ------------------------- 5785 !scalars 5786 integer :: i1,i2,idir1,idir2,index,ipert1,ipert2 5787 real(dp),parameter :: break_symm=1.0d-12 5788 !real(dp),parameter :: break_symm=zero 5789 real(dp) :: fac 5790 !arrays 5791 real(dp) :: nearidentity(3,3) 5792 5793 ! ********************************************************************* 5794 5795 !This slight breaking of the symmetry allows the 5796 !results to be more portable between machines 5797 nearidentity(:,:)=one 5798 nearidentity(1,1)=one+break_symm 5799 nearidentity(3,3)=one-break_symm 5800 5801 !Include the masses in the dynamical matrix 5802 do ipert1=1,natom 5803 do ipert2=1,natom 5804 fac=1.0_dp/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass 5805 do idir1=1,3 5806 do idir2=1,3 5807 i1=idir1+(ipert1-1)*3 5808 i2=idir2+(ipert2-1)*3 5809 index=i1+3*natom*(i2-1) 5810 mat(2*index-1)=mat(2*index-1)*fac*nearidentity(idir1,idir2) 5811 mat(2*index )=mat(2*index )*fac*nearidentity(idir1,idir2) 5812 ! This is to break slightly the translation invariance, and make 5813 ! the automatic tests more portable 5814 if(ipert1==ipert2 .and. idir1==idir2)then 5815 mat(2*index-1)=mat(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp 5816 end if 5817 end do 5818 end do 5819 end do 5820 end do 5821 5822 ! Make the dynamical matrix hermitian 5823 call mkherm(mat,3*natom) 5824 5825 end subroutine massmult_and_breaksym
m_dynmat/nanal9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
nanal9
FUNCTION
If plus=0 then substracts the non-analytical part from one dynamical
matrices, with number iqpt.
If plus=1 then adds the non-analytical part to the dynamical
matrices, with number iqpt.
INPUTS
dyew(2,3,natom,3,natom)= Non-analytical part natom= Number of atoms in the unit cell iqpt= Referenced q point for the dynamical matrix nqpt= Number of q points plus= (see above)
OUTPUT
dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
PARENTS
m_dynmat,m_ifc
CHILDREN
SOURCE
5240 subroutine nanal9(dyew,dynmat,iqpt,natom,nqpt,plus) 5241 5242 5243 !This section has been created automatically by the script Abilint (TD). 5244 !Do not modify the following lines by hand. 5245 #undef ABI_FUNC 5246 #define ABI_FUNC 'nanal9' 5247 !End of the abilint section 5248 5249 implicit none 5250 5251 !Arguments ------------------------------- 5252 !scalars 5253 integer,intent(in) :: iqpt,natom,nqpt,plus 5254 !arrays 5255 real(dp),intent(in) :: dyew(2,3,natom,3,natom) 5256 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt) 5257 5258 !Local variables ------------------------- 5259 !scalars 5260 integer :: ia,ib,mu,nu 5261 character(len=500) :: message 5262 5263 ! ********************************************************************* 5264 5265 if (plus==0) then 5266 5267 do ia=1,natom 5268 do ib=1,natom 5269 do mu=1,3 5270 do nu=1,3 5271 ! The following four lines are OK 5272 dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) - dyew(1,mu,ia,nu,ib) 5273 dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) - dyew(2,mu,ia,nu,ib) 5274 end do 5275 end do 5276 end do 5277 end do 5278 5279 else if (plus==1) then 5280 ! write(std_out,*)' nanal9 : now, only the analytic part ' 5281 ! write(std_out,*)' nanal9 : now, only the ewald part ' 5282 do ia=1,natom 5283 do ib=1,natom 5284 do mu=1,3 5285 do nu=1,3 5286 dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) + dyew(1,mu,ia,nu,ib) 5287 dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) + dyew(2,mu,ia,nu,ib) 5288 end do 5289 end do 5290 end do 5291 end do 5292 5293 else 5294 write(message,'(a,a,a,i0,a)' )& 5295 & 'The argument "plus" must be equal to 0 or 1.',ch10,& 5296 & 'The value ',plus,' is not available.' 5297 MSG_BUG(message) 5298 end if 5299 5300 end subroutine nanal9
m_dynmat/q0dy3_apply [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
q0dy3_apply
FUNCTION
Takes care of the inclusion of the ewald q=0 term in the dynamical matrix - corrects the dyew matrix provided as input
INPUTS
dyewq0(3,3,natom) = part needed to correct
the dynamical matrix for atom self-interaction.
natom= number of atom in the unit cell
OUTPUT
SIDE EFFECTS
dyew(2,3,natom,3,natom)= dynamical matrix corrected on output
NOTES
Should be used just after each call to dfpt_ewald, for both q==0 and the real wavelength. The q0dy3_apply should be used in conjunction with the subroutine dfpt_ewald (or ewald9): First, the call of dfpt_ewald with q==0 should be done, then the call to q0dy3_calc will produce the dyewq0 matrix from the (q=0) dyew matrix Second, the call of dfpt_ewald with the real q (either =0 or diff 0) should be done, then the call to q0dy3_apply will produce the correct dynamical matrix dyew starting from the previously calculated dyewq0 and the bare(non-corrected) dyew matrix
PARENTS
m_dynmat,m_ifc,respfn
CHILDREN
SOURCE
2000 subroutine q0dy3_apply(natom,dyewq0,dyew) 2001 2002 2003 !This section has been created automatically by the script Abilint (TD). 2004 !Do not modify the following lines by hand. 2005 #undef ABI_FUNC 2006 #define ABI_FUNC 'q0dy3_apply' 2007 !End of the abilint section 2008 2009 implicit none 2010 2011 !Arguments ------------------------------- 2012 !scalars 2013 integer,intent(in) :: natom 2014 !arrays 2015 real(dp),intent(in) :: dyewq0(3,3,natom) 2016 real(dp),intent(inout) :: dyew(2,3,natom,3,natom) 2017 2018 !Local variables ------------------------- 2019 !scalars 2020 integer :: ia,mu,nu 2021 2022 ! ********************************************************************* 2023 2024 do mu=1,3 2025 do nu=1,3 2026 do ia=1,natom 2027 dyew(1,mu,ia,nu,ia)=dyew(1,mu,ia,nu,ia)-dyewq0(mu,nu,ia) 2028 end do 2029 end do 2030 end do 2031 2032 end subroutine q0dy3_apply
m_dynmat/q0dy3_calc [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
q0dy3_calc
FUNCTION
Calculate the q=0 correction term to the dynamical matrix
INPUTS
dyew(2,3,natom,3,natom)= dynamical matrix
input, non-corrected, for q=0 if option=1 or 2
natom= number of atom in the unit cell
option= either 1 or 2:
1: use dyew to calculate dyewq0 symmetrical form
2: use dyew to calculate dyewq0 symmetrical form
OUTPUT
dyewq0(3,3,natom) = part needed to correct
the dynamical matrix for atom self-interaction.
NOTES
Should be used just after each call to dfpt_ewald, for both q==0 and the real wavelength. If option=1 or 2, q0dy3_calc uses an Ewald dynamical matrix at q=0, called dyew, to produce a contracted form called dyewq0 : either: in an unsymmetrical form (if option=1), or in a symmetrical form (if option=2). The q0dy3_calc should be used in conjunction with the subroutine dfpt_ewald (or ewald9). First, the call of dfpt_ewald with q==0 should be done , then the call to q0dy3_calc will produce the dyewq0 matrix from the (q=0) dyew matrix Second, the call of dfpt_ewald with the real q (either =0 or diff 0) should be done, then the call to q0dy3_apply will produce the correct dynamical matrix dyew starting from the previously calculated dyewq0 and the bare(non-corrected) dyew matrix
PARENTS
ddb_hybrid,m_ifc,respfn
CHILDREN
SOURCE
2084 subroutine q0dy3_calc(natom,dyewq0,dyew,option) 2085 2086 2087 !This section has been created automatically by the script Abilint (TD). 2088 !Do not modify the following lines by hand. 2089 #undef ABI_FUNC 2090 #define ABI_FUNC 'q0dy3_calc' 2091 !End of the abilint section 2092 2093 implicit none 2094 2095 !Arguments ------------------------------- 2096 !scalars 2097 integer,intent(in) :: natom,option 2098 !arrays 2099 real(dp),intent(in) :: dyew(2,3,natom,3,natom) 2100 real(dp),intent(out) :: dyewq0(3,3,natom) 2101 2102 !Local variables ------------------------- 2103 !scalars 2104 integer :: ia,ib,mu,nu 2105 character(len=500) :: message 2106 2107 ! ********************************************************************* 2108 2109 if(option==1.or.option==2)then 2110 do mu=1,3 2111 do nu=1,3 2112 do ia=1,natom 2113 dyewq0(mu,nu,ia)=zero 2114 do ib=1,natom 2115 dyewq0(mu,nu,ia)=dyewq0(mu,nu,ia)+dyew(1,mu,ia,nu,ib) 2116 end do 2117 end do 2118 end do 2119 end do 2120 else 2121 write (message, '(3a)')& 2122 & 'option should be 1 or 2.',ch10,& 2123 & 'action: correct calling routine' 2124 MSG_BUG(message) 2125 end if 2126 2127 if(option==2)then 2128 do ia=1,natom 2129 do mu=1,3 2130 do nu=mu,3 2131 dyewq0(mu,nu,ia)=(dyewq0(mu,nu,ia)+dyewq0(nu,mu,ia))/2 2132 dyewq0(nu,mu,ia)=dyewq0(mu,nu,ia) 2133 end do 2134 end do 2135 end do 2136 end if 2137 2138 end subroutine q0dy3_calc
m_dynmat/symdm9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
symdm9
FUNCTION
Use the set of special k points calculated by the Monkhorst & Pack Technique. Check if all the informations for the k points are present in the DDB to determine their dynamical matrices. Generate the dynamical matrices of the set of k points which samples homogeneously the entire Brillouin zone.
INPUTS
blkflg(nsize,nblok)= flag of existence for each element of the DDB blknrm(1,nblok)=norm of qpt providing normalization blkqpt(1<ii<9,nblok)=q vector of a phonon mode (ii=1,2,3) blktyp(nblok)=1 or 2 depending on non-stationary or stationary block 3 for third order derivatives blkval(2,3*mpert*3*mpert,nblok)= all the dynamical matrices gprim(3,3)=dimensionless primitive translations in reciprocal space indsym = mapping of atoms under symops mpert =maximum number of ipert natom=number of atoms in unit cell nblok=number of blocks in the DDB nqpt=number of special q points nsym=number of space group symmetries rfmeth = 1 if non-stationary block 2 if stationary block 3 if third order derivatives rprim(3,3)=dimensionless primitive translations in real space spqpt(3,nqpt)=set of special q points generated by the Monkhorst & Pack Method symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
OUTPUT
dynmat(2,3,natom,3,natom,nqpt)=dynamical matrices relative to the q points of the B.Z. sampling
[qmissing]=Allocatable array with the indices of the q-points in the BZ that could not be obtained
by symmetry. If qmissing is present, the routine does not stop if the full BZ cannot be reconstructed.
The caller is responsible for filling the missing entries.
TODO
* A full description of the inputs should be included
NOTES
Time-reversal symmetry is always assumed
PARENTS
m_ifc
CHILDREN
SOURCE
4729 subroutine symdm9(blkflg,blknrm,blkqpt,blktyp,blkval,& 4730 & dynmat,gprim,indsym,mpert,natom,nblok,nqpt,nsym,rfmeth,& 4731 & rprim,spqpt,symrec,symrel,qmissing) 4732 4733 4734 !This section has been created automatically by the script Abilint (TD). 4735 !Do not modify the following lines by hand. 4736 #undef ABI_FUNC 4737 #define ABI_FUNC 'symdm9' 4738 use interfaces_14_hidewrite 4739 !End of the abilint section 4740 4741 implicit none 4742 4743 !Arguments ------------------------------- 4744 !scalars 4745 integer,intent(in) :: mpert,natom,nblok,nqpt,nsym,rfmeth 4746 !arrays 4747 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),blktyp(nblok) 4748 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4749 integer,allocatable,optional,intent(out) :: qmissing(:) 4750 real(dp),intent(in) :: blknrm(3,nblok),blkqpt(9,nblok) 4751 real(dp),intent(in) :: blkval(2,3*mpert*3*mpert,nblok),gprim(3,3),rprim(3,3) 4752 real(dp),intent(in) :: spqpt(3,nqpt) 4753 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt) 4754 4755 !Local variables ------------------------- 4756 !tol sets tolerance for equality of q points between those of 4757 !the DDB and those of the sampling grid 4758 !scalars 4759 integer :: ia,ib,iblok,idir1,idir2,ii,ipert1,ipert2,iqpt,isym,jj,kk,ll 4760 integer :: mu,nu,q1,q2,nqmiss 4761 real(dp),parameter :: tol=2.d-8 4762 real(dp) :: arg1,arg2,im,re,sumi,sumr 4763 logical :: allow_qmiss 4764 character(len=500) :: message 4765 !arrays 4766 integer,allocatable :: qtest(:,:) 4767 integer :: qmiss_(nqpt) 4768 real(dp) :: qq(3),qsym(6),ss(3,3) 4769 real(dp),allocatable :: ddd(:,:,:,:,:) 4770 4771 ! ********************************************************************* 4772 ! Initialize output (some q-points might not be reconstructed if qmissing is present) 4773 dynmat = huge(one) 4774 allow_qmiss = (present(qmissing)) 4775 4776 ABI_ALLOCATE(ddd,(2,3,natom,3,natom)) 4777 !Check if the blkqpt points and their symmetrics are sufficient 4778 !in the DDB to retrieve all the q points of the B.Z. sampling 4779 4780 !Initialization of a test variable 4781 ! qtest(iqpt,1)=iblok 4782 ! qtest(iqpt,2)=isym 4783 ! qtest(iqpt,3)=time_reversal 4784 ABI_ALLOCATE(qtest,(nqpt,3)) 4785 do iqpt=1,nqpt 4786 qtest(iqpt,1)=0 4787 end do 4788 4789 !Q points coming from the DDB 4790 !write(std_out,*)' Nbr. of Blocks -> ',nblok 4791 4792 do iblok=1,nblok 4793 4794 if (blktyp(iblok)==rfmeth) then 4795 qq(1)=blkqpt(1,iblok)/blknrm(1,iblok) 4796 qq(2)=blkqpt(2,iblok)/blknrm(1,iblok) 4797 qq(3)=blkqpt(3,iblok)/blknrm(1,iblok) 4798 4799 ! Calculation of the symmetric points (including Time Reversal) 4800 do isym=1,nsym 4801 qsym(1)=qq(1)*symrec(1,1,isym)+qq(2)*symrec(1,2,isym)+qq(3)*symrec(1,3,isym) 4802 qsym(2)=qq(1)*symrec(2,1,isym)+qq(2)*symrec(2,2,isym)+qq(3)*symrec(2,3,isym) 4803 qsym(3)=qq(1)*symrec(3,1,isym)+qq(2)*symrec(3,2,isym)+qq(3)*symrec(3,3,isym) 4804 4805 ! Dont forget the Time Reversal symmetry 4806 qsym(4)=-qq(1)*symrec(1,1,isym)-qq(2)*symrec(1,2,isym)-qq(3)*symrec(1,3,isym) 4807 qsym(5)=-qq(1)*symrec(2,1,isym)-qq(2)*symrec(2,2,isym)-qq(3)*symrec(2,3,isym) 4808 qsym(6)=-qq(1)*symrec(3,1,isym)-qq(2)*symrec(3,2,isym)-qq(3)*symrec(3,3,isym) 4809 4810 ! Comparison between the q points and their symmetric points 4811 ! and the set of q points which samples the entire Brillouin zone 4812 do iqpt=1,nqpt 4813 4814 if (mod(abs(spqpt(1,iqpt)-qsym(1))+tol,1._dp)<2*tol)then 4815 if (mod(abs(spqpt(2,iqpt)-qsym(2))+tol,1._dp)<2*tol)then 4816 if (mod(abs(spqpt(3,iqpt)-qsym(3))+tol,1._dp)<2*tol)then 4817 4818 ! write(std_out,*)' q point from the DDB ! ' 4819 ! write(std_out,*)' block -> ',iblok 4820 ! write(std_out,*)' sym. -> ',isym 4821 ! write(std_out,*)' No Time Reversal ' 4822 ! write(std_out,*)'(',qsym(1),',',qsym(2),',',qsym(3),')' 4823 ! write(std_out,*)' ' 4824 4825 qtest(iqpt,1)=iblok 4826 qtest(iqpt,2)=isym 4827 qtest(iqpt,3)=0 4828 end if 4829 end if 4830 end if 4831 4832 if (mod(abs(spqpt(1,iqpt)-qsym(4))+tol,1._dp)<2*tol)then 4833 if (mod(abs(spqpt(2,iqpt)-qsym(5))+tol,1._dp)<2*tol)then 4834 if (mod(abs(spqpt(3,iqpt)-qsym(6))+tol,1._dp)<2*tol)then 4835 4836 ! write(std_out,*)' q point from the DDB ! ' 4837 ! write(std_out,*)' block -> ',iblok 4838 ! write(std_out,*)' sym. -> ',isym 4839 ! write(std_out,*)' Time Reversal ' 4840 ! write(std_out,*)'(',qsym(4),',',qsym(5),',',qsym(6),')' 4841 ! write(std_out,*)' ' 4842 4843 qtest(iqpt,1)=iblok 4844 qtest(iqpt,2)=isym 4845 qtest(iqpt,3)=1 4846 end if 4847 end if 4848 end if 4849 4850 end do ! End of the loop on the q points of the sampling 4851 end do ! End of the loop on the symmetries 4852 4853 end if 4854 end do ! End of the loop on the q points of the DDB 4855 4856 ! Check if all the informations relatives to the q points sampling are found in the DDB; 4857 ! if not => stop message 4858 nqmiss = 0 4859 do iqpt=1,nqpt 4860 if (qtest(iqpt,1)==0) then 4861 nqmiss = nqmiss + 1 4862 qmiss_(nqmiss) = iqpt 4863 write(message, '(a,a,a)' )& 4864 & ' symdm9 : the bloks found in the DDB are characterized',ch10,& 4865 & ' by the following wavevectors :' 4866 call wrtout(std_out,message,'COLL') 4867 do iblok=1,nblok 4868 write(message, '(a,4d20.12)')& 4869 & ' ',blkqpt(1,iblok),blkqpt(2,iblok),blkqpt(3,iblok),blknrm(1,iblok) 4870 call wrtout(std_out,message,'COLL') 4871 end do 4872 write(message, '(a,a,a,i0,a,a,a,3es16.6,a,a,a,a)' )& 4873 & 'Information is missing in the DDB.',ch10,& 4874 & 'The dynamical matrix number ',iqpt,' cannot be built,',ch10,& 4875 & 'since no block with wavevector',spqpt(1:3,iqpt),ch10,& 4876 & 'has been found.',ch10,& 4877 & 'Action: add the required blok in the DDB, or modify your input file.' 4878 if (.not.allow_qmiss) then 4879 MSG_ERROR(message) 4880 else 4881 !continue 4882 MSG_COMMENT(message) 4883 end if 4884 end if 4885 end do 4886 4887 ! Will return a list with the index of the q-points that could not be symmetrized. 4888 if (allow_qmiss) then 4889 ABI_ALLOCATE(qmissing, (nqmiss)) 4890 if (nqmiss > 0) qmissing = qmiss_(1:nqmiss) 4891 end if 4892 4893 !Generation of the dynamical matrices relative to the q points 4894 !of the set which samples the entire Brillouin zone 4895 do iqpt=1,nqpt 4896 q1=qtest(iqpt,1) 4897 q2=qtest(iqpt,2) 4898 ! Skip this q-point if don't have enough info and allow_qmiss 4899 if (allow_qmiss .and. q1==0) cycle 4900 4901 ! Check if the symmetry accompagnied with time reversal : q <- -q 4902 if (qtest(iqpt,3)==0) then 4903 do ii=1,3 4904 qq(ii)=blkqpt(ii,q1)/blknrm(1,q1) 4905 do jj=1,3 4906 ss(ii,jj)=zero 4907 do kk=1,3 4908 do ll=1,3 4909 ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2) 4910 end do 4911 end do 4912 end do 4913 end do 4914 ! write(std_out,*)"ss",ss 4915 else 4916 do ii=1,3 4917 qq(ii)=-blkqpt(ii,q1)/blknrm(1,q1) 4918 do jj=1,3 4919 ss(ii,jj)=zero 4920 do kk=1,3 4921 do ll=1,3 4922 ss(ii,jj)=ss(ii,jj)+rprim(ii,kk)*gprim(jj,ll)*symrel(kk,ll,q2) 4923 end do 4924 end do 4925 end do 4926 end do 4927 end if 4928 ! 4929 ! Check whether all the information is contained in the DDB 4930 do ipert2=1,natom 4931 do idir2=1,3 4932 do ipert1=1,natom 4933 do idir1=1,3 4934 if(blkflg(idir1,ipert1,idir2,ipert2,q1)/=1)then 4935 write(message, '(a,a,a,i6,a,a,a,4i4,a,a,a,a)' )& 4936 & 'Informations are missing in the DDB.',ch10,& 4937 & 'In block',q1,' the following element is missing :',ch10,& 4938 & 'idir1,ipert1,idir2,ipert2=',idir1,ipert1,idir2,ipert2,ch10,& 4939 & 'Action: add the required information in the DDB,',ch10,& 4940 & 'or modify your input file.' 4941 MSG_ERROR(message) 4942 end if 4943 end do 4944 end do 4945 end do 4946 end do 4947 4948 ! Read the dynamical matrices in the DDB 4949 do ipert2=1,natom 4950 do idir2=1,3 4951 do ipert1=1,natom 4952 do idir1=1,3 4953 ddd(:,idir1,ipert1,idir2,ipert2)=blkval(:,idir1+3*(ipert1-1+mpert*(idir2-1+3*(ipert2-1))),q1) 4954 end do 4955 end do 4956 end do 4957 end do 4958 4959 ! Calculation of the dynamical matrix of a symmetrical q point 4960 do ia=1,natom 4961 do ib=1,natom 4962 ! write(std_out,*)'atom-> ',ia,indsym(4,q2,ia) 4963 ! write(std_out,*)'atom-> ',ib,indsym(4,q2,ib) 4964 arg1=two_pi*(qq(1)*indsym(1,q2,ia)+qq(2)*indsym(2,q2,ia)+qq(3)*indsym(3,q2,ia)) 4965 arg2=two_pi*(qq(1)*indsym(1,q2,ib)+qq(2)*indsym(2,q2,ib)+qq(3)*indsym(3,q2,ib)) 4966 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 4967 im=cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2) 4968 ! write(std_out,*)'re : ',re, 'im : ',im 4969 do mu=1,3 4970 do nu=1,3 4971 ! write(std_out,*)' ' 4972 sumr=zero 4973 sumi=zero 4974 do ii=1,3 4975 do jj=1,3 4976 ! If there is Time Reversal : D.M. <- Complex Conjugate D.M. 4977 if (qtest(iqpt,3)==0) then 4978 sumr=sumr+ss(mu,ii)*ss(nu,jj)*ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 4979 sumi=sumi+ss(mu,ii)*ss(nu,jj)*ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 4980 else 4981 sumr=sumr+ss(mu,ii)*ss(nu,jj)*ddd(1,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 4982 sumi=sumi-ss(mu,ii)*ss(nu,jj)*ddd(2,ii,indsym(4,q2,ia),jj,indsym(4,q2,ib)) 4983 end if 4984 end do 4985 end do 4986 4987 ! Dynmat -> Dynamical Matrix for the q point of the sampling 4988 ! write(std_out,*)' Sumr -> ',mu,nu,sumr 4989 ! write(std_out,*)' Sumi -> ',mu,nu,sumi 4990 dynmat(1,mu,ia,nu,ib,iqpt)=re*sumr-im*sumi 4991 dynmat(2,mu,ia,nu,ib,iqpt)=re*sumi+im*sumr 4992 4993 ! DEBUG 4994 ! if((ia==2 .or. ia==3) .and. ib==1)then 4995 ! write(std_out,'(5i3,2es16.8)' )mu,ia,nu,ib,iqpt,dynmat(1:2,mu,ia,nu,ib,iqpt) 4996 ! end if 4997 ! ENDDEBUG 4998 4999 end do ! End loop on the coordinates 5000 end do 5001 5002 end do ! End loop on the ia atoms 5003 end do ! End loop on the ib atoms 5004 end do ! End loop on the q points of the sampling 5005 5006 ABI_DEALLOCATE(ddd) 5007 ABI_DEALLOCATE(qtest) 5008 5009 end subroutine symdm9
m_dynmat/symdyma [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
symdyma
FUNCTION
Symmetrize the dynamical matrices
INPUTS
indsym(4,nsym*natom)=indirect indexing array : for each isym,iatom, fourth element is label of atom into which iatom is sent by INVERSE of symmetry operation isym; first three elements are the primitive translations which must be subtracted after the transformation to get back to the original unit cell. natom=number of atoms in unit cell nsym=number of space group symmetries qptn(3)=normalized phonon wavevector rprimd(3,3)=dimensional primitive translations (bohr) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
SIDE EFFECTS
Input/Output dmati(2*3*natom*3*natom)=dynamical matrices in cartesian coordinates relative to the q points of the B.Z. sampling
NOTES
the procedure of the symmetrization of the dynamical matrix follows the equations in: Hendrikse et al., Computer Phys. Comm. 86, 297 (1995)
TODO
A full description of the equations should be included
PARENTS
m_dynmat,m_phgamma,relaxpol
CHILDREN
SOURCE
2182 subroutine symdyma(dmati,indsym,natom,nsym,qptn,rprimd,symrel,symafm) 2183 2184 2185 !This section has been created automatically by the script Abilint (TD). 2186 !Do not modify the following lines by hand. 2187 #undef ABI_FUNC 2188 #define ABI_FUNC 'symdyma' 2189 use interfaces_32_util 2190 !End of the abilint section 2191 2192 implicit none 2193 2194 !Arguments ------------------------------- 2195 !scalars 2196 integer,intent(in) :: natom,nsym 2197 !arrays 2198 integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym) 2199 integer,intent(in) :: symafm(nsym) 2200 real(dp),intent(in) :: qptn(3),rprimd(3,3) 2201 real(dp),intent(inout) :: dmati(2*3*natom*3*natom) 2202 2203 !Local variables ------------------------- 2204 !scalars 2205 integer :: i1,i2,iat,idir,ii,index,isgn,isym,itirev,jat,jdir,jj,kk,ll 2206 integer :: niat,njat,timrev 2207 real(dp) :: arg1,arg2,dmint,im,re,sumi,sumr 2208 !arrays 2209 integer :: indij(natom,natom),symq(4,2,nsym),symrec(3,3,nsym) 2210 real(dp) :: TqR(3,3),TqS_(3,3),dynmat(2,3,natom,3,natom) 2211 real(dp) :: dynmatint(2*nsym,2,3,natom,3,natom),gprimd(3,3) 2212 real(dp) :: symcart(3,3,nsym) 2213 2214 ! ********************************************************************* 2215 2216 !0) initializations 2217 call matr3inv(rprimd,gprimd) 2218 do isym=1,nsym 2219 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 2220 end do 2221 2222 TqR=zero 2223 TqS_=zero 2224 dynmat=zero 2225 2226 !Note: dynmat is used as work space here 2227 i1=0 2228 do iat=1,natom 2229 do idir=1,3 2230 i1=i1+1 2231 i2=0 2232 do jat=1,natom 2233 do jdir=1,3 2234 i2=i2+1 2235 index=i1+3*natom*(i2-1) 2236 dynmat(1,idir,iat,jdir,jat)=dmati(2*index-1) 2237 dynmat(2,idir,iat,jdir,jat)=dmati(2*index ) 2238 end do 2239 end do 2240 end do 2241 end do 2242 2243 !Transform symrel to cartesian coordinates (RC coding) 2244 !do isym=1,nsym 2245 !symcart(:,:,isym)=matmul(rprimd,matmul(dble(symrel(:,:,isym)),gprimd)) 2246 !end do 2247 2248 !Coding from symdm9 2249 do isym=1,nsym 2250 do jj=1,3 2251 symcart(:,jj,isym)=zero 2252 do kk=1,3 2253 do ll=1,3 2254 symcart(:,jj,isym)=symcart(:,jj,isym)+rprimd(:,kk)*gprimd(jj,ll)*symrel(kk,ll,isym) 2255 end do 2256 end do 2257 end do 2258 end do 2259 2260 2261 !Get the symq of the CURRENT Q POINT 2262 !mjv: set prtvol=0 for production runs. 2263 call littlegroup_q(nsym,qptn,symq,symrec,symafm,timrev,prtvol=0) 2264 2265 indij(:,:)=0 2266 dynmatint=zero 2267 2268 do isym=1,nsym ! loop over all the symmetries 2269 ! write(std_out,*) 'current symmetry',isym 2270 do itirev=1,2 ! loop over the time-reversal symmetry 2271 isgn=3-2*itirev ! to take into accont the time-reversal 2272 ! write(std_out,*) 'timereversal',isgn 2273 if (symq(4,itirev,isym)==1) then ! isym belongs to the wave vector point group 2274 ! write(std_out,*) 'isym belongs to the wave vector point group' 2275 do iat=1,natom ! loop over the atoms 2276 do jat=1,natom ! loop over the atoms 2277 niat=indsym(4,isym,iat) ! niat={R|t}iat 2278 njat=indsym(4,isym,jat) ! njat={R|t}jat 2279 indij(niat,njat)=indij(niat,njat)+1 2280 ! write(std_out,'(a,5i5)') 'current status:',iat,jat,niat,njat,indij(niat,njat) 2281 ! phase calculation, arg1 and arg2 because of two-atom derivative 2282 arg1=two_pi*( qptn(1)*indsym(1,isym,iat)+& 2283 & qptn(2)*indsym(2,isym,iat)+& 2284 & qptn(3)*indsym(3,isym,iat) ) 2285 arg2=two_pi*( qptn(1)*indsym(1,isym,jat)+& 2286 & qptn(2)*indsym(2,isym,jat)+& 2287 & qptn(3)*indsym(3,isym,jat) ) 2288 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 2289 im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)) 2290 2291 do idir=1,3 ! loop over displacements 2292 do jdir=1,3 ! loop over displacements 2293 ! we pick the (iat,jat) (3x3) block of the dyn.mat. 2294 sumr=zero 2295 sumi=zero 2296 do ii=1,3 2297 do jj=1,3 2298 sumr=sumr+symcart(idir,ii,isym)*dynmat(1,ii,niat,jj,njat)*symcart(jdir,jj,isym) 2299 sumi=sumi+symcart(idir,ii,isym)*dynmat(2,ii,niat,jj,njat)*symcart(jdir,jj,isym) 2300 end do 2301 end do 2302 sumi=isgn*sumi 2303 2304 dynmatint(nsym*(itirev-1)+isym,1,idir,iat,jdir,jat)=re*sumr-im*sumi 2305 dynmatint(nsym*(itirev-1)+isym,2,idir,iat,jdir,jat)=re*sumi+im*sumr 2306 2307 end do 2308 end do 2309 end do 2310 end do ! end treatment of the (iat,jat) (3x3) block of dynmat 2311 end if ! symmetry check 2312 end do ! time-reversal 2313 end do ! symmetries 2314 2315 2316 !4) make the average, get the final symmetric dynamical matrix 2317 do iat=1,natom 2318 do jat=1,natom 2319 do idir=1,3 2320 do jdir=1,3 2321 dmint=zero 2322 do isym=1,2*nsym 2323 dmint=dmint+dynmatint(isym,1,idir,iat,jdir,jat) 2324 end do 2325 dynmat(1,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat)) 2326 dmint=zero 2327 do isym=1,2*nsym 2328 dmint=dmint+dynmatint(isym,2,idir,iat,jdir,jat) 2329 end do 2330 dynmat(2,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat)) 2331 end do 2332 end do 2333 end do 2334 end do 2335 2336 i1=0 2337 do iat=1,natom 2338 do idir=1,3 2339 i1=i1+1 2340 i2=0 2341 do jat=1,natom 2342 do jdir=1,3 2343 i2=i2+1 2344 index=i1+3*natom*(i2-1) 2345 dmati(2*index-1)=dynmat(1,idir,iat,jdir,jat) 2346 dmati(2*index )=dynmat(2,idir,iat,jdir,jat) 2347 end do 2348 end do 2349 end do 2350 end do 2351 2352 end subroutine symdyma
m_dynmat/sytens [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
sytens
FUNCTION
Determines the set of irreductible elements of the non-linear optical susceptibility and Raman tensors
INPUTS
indsym(4,nsym,natom)=indirect indexing array described above: for each isym,iatom, fourth element is label of atom into which iatom is sent by INVERSE of symmetry operation isym; first three elements are the primitive translations which must be subtracted after the transformation to get back to the original unit cell. mpert =maximum number of ipert natom= number of atoms nsym=number of space group symmetries symrec(3,3,nsym)=3x3 matrices of the group symmetries (reciprocal space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real space)
OUTPUT
(see side effects)
SIDE EFFECTS
rfpert(3,mpert,3,mpert,3,mpert) = array defining the type of perturbations
that have to be computed
At the input :
1 -> element has to be computed explicitely
At the output :
1 -> element has to be computed explicitely
-1 -> use symmetry operations to obtain the corresponding element
-2 -> element is zero by symmetry
PARENTS
m_ddb,nonlinear,respfn
CHILDREN
SOURCE
4457 subroutine sytens(indsym,mpert,natom,nsym,rfpert,symrec,symrel) 4458 4459 use defs_basis 4460 use m_profiling_abi 4461 4462 !This section has been created automatically by the script Abilint (TD). 4463 !Do not modify the following lines by hand. 4464 #undef ABI_FUNC 4465 #define ABI_FUNC 'sytens' 4466 !End of the abilint section 4467 4468 implicit none 4469 4470 !Arguments ------------------------------- 4471 !scalars 4472 integer,intent(in) :: mpert,natom,nsym 4473 !arrays 4474 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4475 integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert) 4476 4477 !Local variables ------------------------- 4478 !scalars 4479 integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_ 4480 integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2 4481 integer :: ipesy3,isym 4482 !arrays 4483 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 4484 integer,allocatable :: pertsy(:,:,:,:,:,:) 4485 4486 !*********************************************************************** 4487 4488 !DEBUG 4489 !write(std_out,*)'sytens : enter' 4490 !write(std_out,*)'indsym = ' 4491 !write(std_out,*)indsym 4492 !stop 4493 !ENDDEBUG 4494 4495 ABI_ALLOCATE(pertsy,(3,mpert,3,mpert,3,mpert)) 4496 pertsy(:,:,:,:,:,:) = 0 4497 4498 !Loop over perturbations 4499 4500 do i1pert_ = 1, mpert 4501 do i2pert_ = 1, mpert 4502 do i3pert_ = 1, mpert 4503 4504 do i1dir_ = 1, 3 4505 do i2dir_ = 1, 3 4506 do i3dir_ = 1, 3 4507 4508 i1pert = (mpert - i1pert_ + 1) 4509 if (i1pert <= natom) i1pert = natom + 1 - i1pert 4510 i2pert = (mpert - i2pert_ + 1) 4511 if (i2pert <= natom) i2pert = natom + 1 - i2pert 4512 i3pert = (mpert - i3pert_ + 1) 4513 if (i3pert <= natom) i3pert = natom + 1 - i3pert 4514 4515 if (i1pert <= natom) then 4516 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 4517 else if (i2pert <= natom) then 4518 i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_ 4519 else if (i3pert <= natom) then 4520 i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_ 4521 else 4522 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 4523 end if 4524 4525 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then 4526 4527 ! Loop over all symmetries 4528 4529 flag = 0 4530 do isym = 1, nsym 4531 4532 found = 1 4533 4534 ! Select the symmetric element of i1pert,i2pert,i3pert 4535 4536 if (i1pert <= natom) then 4537 ipesy1 = indsym(4,isym,i1pert) 4538 sym1(:,:) = symrec(:,:,isym) 4539 else if (i1pert == natom + 2) then 4540 ipesy1 = i1pert 4541 sym1(:,:) = symrel(:,:,isym) 4542 else 4543 found = 0 4544 end if 4545 4546 if (i2pert <= natom) then 4547 ipesy2 = indsym(4,isym,i2pert) 4548 sym2(:,:) = symrec(:,:,isym) 4549 else if (i2pert == natom + 2) then 4550 ipesy2 = i2pert 4551 sym2(:,:) = symrel(:,:,isym) 4552 else 4553 found = 0 4554 end if 4555 4556 if (i3pert <= natom) then 4557 ipesy3 = indsym(4,isym,i3pert) 4558 sym3(:,:) = symrec(:,:,isym) 4559 else if (i3pert == natom + 2) then 4560 ipesy3 = i3pert 4561 sym3(:,:) = symrel(:,:,isym) 4562 else 4563 found = 0 4564 end if 4565 4566 ! See if the symmetric element is available and check if some 4567 ! of the elements may be zeor. In the latter case, they do not need 4568 ! to be computed. 4569 4570 4571 if ((flag /= -1).and.& 4572 & (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then 4573 flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir) 4574 end if 4575 4576 4577 do idisy1 = 1, 3 4578 do idisy2 = 1, 3 4579 do idisy3 = 1, 3 4580 4581 if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.& 4582 & (sym3(i3dir,idisy3) /= 0)) then 4583 if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then 4584 found = 0 4585 ! exit ! exit loop over symmetries 4586 end if 4587 end if 4588 4589 4590 if ((flag == -1).and.& 4591 & ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then 4592 if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.& 4593 & (sym3(i3dir,idisy3)/=0)) then 4594 flag = 0 4595 end if 4596 end if 4597 4598 4599 4600 end do 4601 end do 4602 end do 4603 4604 4605 if (found == 1) then 4606 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1 4607 end if 4608 4609 4610 ! In case a symmetry operation only changes the sign of an 4611 ! element, this element has to be equal to zero 4612 4613 if (flag == -1) then 4614 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2 4615 exit 4616 end if 4617 4618 end do ! close loop on symmetries 4619 4620 4621 4622 ! If the elemetn i1pert,i2pert,i3pert is not symmetric 4623 ! to a basis element, it is a basis element 4624 4625 if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then 4626 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 4627 end if 4628 4629 end if ! rfpert /= 0 4630 4631 end do ! close loop over perturbations 4632 end do 4633 end do 4634 end do 4635 end do 4636 end do 4637 4638 4639 !Now, take into account the permutation of (i1pert,i1dir) 4640 !and (i3pert,i3dir) 4641 4642 do i1pert = 1, mpert 4643 do i2pert = 1, mpert 4644 do i3pert = 1, mpert 4645 4646 do i1dir = 1, 3 4647 do i2dir = 1, 3 4648 do i3dir = 1, 3 4649 4650 if ((i1pert /= i3pert).or.(i1dir /= i3dir)) then 4651 4652 if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.& 4653 & (pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) == 1)) then 4654 pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = -1 4655 end if 4656 4657 end if 4658 4659 end do 4660 end do 4661 end do 4662 4663 end do 4664 end do 4665 end do 4666 4667 rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:) 4668 4669 ABI_DEALLOCATE(pertsy) 4670 4671 4672 end subroutine sytens
m_dynmat/wght9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
wght9
FUNCTION
Generates a weight to each R points of the Big Box and for each pair of atoms For each R points included in the space generates by moving the unit cell around each atom; the weight will be one. Border conditions are provided. The R points outside the chosen space will have a 0 weight.
INPUTS
brav = Bravais lattice (1 or -1=S.C.;2=F.C.C.;4=Hex. -1 is for old algo to find weights, =1 is for Wigner-Seitz algo) gprim(3,3)= Normalized coordinates in reciprocal space natom= Number of atoms in the unit cell ngqpt(6)= Numbers used to sample the Brillouin zone nqpt= Number of q points used in the homogeneous grid sampling the Brillouin zone nqshft=number of shift vectors in the repeated cell nrpt=Number of R points in the Big Box qshft(3,nqshft)=vectors that will be used to determine the shifts from (0. 0. 0.) rcan(3,natom)=Atomic position in canonical coordinates rpt(3,nprt)=Canonical coordinates of the R points in the unit cell These coordinates are normalized (=> * acell(3)) rprimd(3,3)=dimensional primitive translations for real space (bohr)
OUTPUT
wghatm(natom,natom,nrpt)= Weight associated to the couple of atoms and the R vector The vector r(atom2)-r(atom1)+rpt should be inside the moving box ngqpt(6)= can be modified
PARENTS
m_ifc
CHILDREN
SOURCE
3795 subroutine wght9(brav,gprim,natom,ngqpt,nqpt,nqshft,nrpt,qshft,rcan,rpt,rprimd,r_inscribed_sphere,wghatm) 3796 3797 3798 !This section has been created automatically by the script Abilint (TD). 3799 !Do not modify the following lines by hand. 3800 #undef ABI_FUNC 3801 #define ABI_FUNC 'wght9' 3802 use interfaces_14_hidewrite 3803 !End of the abilint section 3804 3805 implicit none 3806 3807 !Arguments ------------------------------- 3808 !scalars 3809 integer,intent(in) :: brav,natom,nqpt,nqshft,nrpt 3810 real(dp),intent(out) :: r_inscribed_sphere 3811 !arrays 3812 integer,intent(inout) :: ngqpt(9) 3813 real(dp),intent(in) :: gprim(3,3),qshft(3,4),rcan(3,natom),rpt(3,nrpt),rprimd(3,3) 3814 real(dp),intent(out) :: wghatm(natom,natom,nrpt) 3815 3816 !Local variables ------------------------- 3817 !scalars 3818 integer :: ia,ib,ii,jj,kk,iqshft,irpt,jqshft,nbordh,tok,new_wght,nptws,nreq 3819 integer :: idir 3820 real(dp), parameter :: tolsym=tol8 3821 real(dp) :: factor,sumwght,normsq,proj 3822 character(len=500) :: message 3823 !arrays 3824 integer :: nbord(9) 3825 real(dp) :: rdiff(9),red(3,3),ptws(4, 729),pp(3),rdiff_tmp(3) 3826 3827 ! ********************************************************************* 3828 3829 DBG_ENTER("COLL") 3830 3831 !First analyze the vectors qshft 3832 if(nqshft/=1)then 3833 3834 if(brav==4)then 3835 write(message,'(a,a,a,i5,a,a,a)' )& 3836 & 'For the time being, only nqshft=1',ch10,& 3837 & 'is allowed with brav=4, while it is nqshft=',nqshft,'.',ch10,& 3838 & 'Action: in the input file, correct either brav or nqshft.' 3839 MSG_ERROR(message) 3840 end if 3841 3842 if(nqshft==2)then 3843 ! Make sure that the q vectors form a BCC lattice 3844 do ii=1,3 3845 if(abs(abs(qshft(ii,1)-qshft(ii,2))-.5_dp)>1.d-10)then 3846 write(message, '(a,a,a,a,a,a,a)' )& 3847 & 'The test of the q1shft vectors shows that they',ch10,& 3848 & 'do not generate a body-centered lattice, which',ch10,& 3849 & 'is mandatory for nqshft=2.',ch10,& 3850 & 'Action: change the q1shft vectors in your input file.' 3851 MSG_ERROR(message) 3852 end if 3853 end do 3854 else if(nqshft==4)then 3855 ! Make sure that the q vectors form a FCC lattice 3856 do iqshft=1,3 3857 do jqshft=iqshft+1,4 3858 tok=0 3859 do ii=1,3 3860 ! Test on the presence of a +-0.5 difference 3861 if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-.5_dp) <1.d-10) tok=tok+1 3862 ! Test on the presence of a 0 or +-1.0 difference 3863 if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-1._dp) <1.d-10 .or.& 3864 & abs(qshft(ii,iqshft)-qshft(ii,jqshft)) < 1.d-10) tok=tok+4 3865 end do 3866 ! Test 1 should be satisfied twice, and test 2 once 3867 if(tok/=6)then 3868 write(message, '(7a)' )& 3869 & 'The test of the q1shft vectors shows that they',ch10,& 3870 & 'do not generate a face-centered lattice, which',ch10,& 3871 & 'is mandatory for nqshft=4.',ch10,& 3872 & 'Action: change the q1shft vectors in your input file.' 3873 MSG_ERROR(message) 3874 end if 3875 end do 3876 end do 3877 else 3878 write(message, '(a,i4,3a)' )& 3879 & 'nqshft must be 1, 2 or 4. It is nqshft=',nqshft,'.',ch10,& 3880 & 'Action: change nqshft in your input file.' 3881 MSG_ERROR(message) 3882 end if 3883 3884 end if 3885 3886 factor=0.5_dp 3887 if(brav==2 .or. brav==3)then 3888 factor=0.25_dp 3889 end if 3890 if(nqshft/=1)factor=factor*2 3891 3892 if (brav==1) then 3893 3894 ! Does not support multiple shifts 3895 if (nqshft/=1) then 3896 write(message, '(a)' ) 'This version of the weights does not support nqshft/=1.' 3897 MSG_ERROR(message) 3898 end if 3899 3900 ! Find the points of the lattice given by ngqpt*acell. These are used to define 3901 ! a Wigner-Seitz cell around the origin. The origin is excluded from the list. 3902 ! TODO : in principle this should be only -1 to +1 for ii jj kk! 3903 nptws=0 3904 do ii=-2,2 3905 do jj=-2,2 3906 do kk=-2,2 3907 do idir=1,3 3908 pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3) 3909 end do 3910 normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) 3911 if (normsq > tol6) then 3912 nptws = nptws + 1 3913 ptws(:3,nptws) = pp(:) 3914 ptws(4,nptws) = half*normsq 3915 end if 3916 end do 3917 end do 3918 end do 3919 end if ! end new_wght 3920 3921 !DEBUG 3922 !write(std_out,*)'factor,ngqpt',factor,ngqpt(1:3) 3923 !ENDDEBUG 3924 3925 r_inscribed_sphere = sum((matmul(rprimd(:,:),ngqpt(1:3)))**2) 3926 do ii=-1,1 3927 do jj=-1,1 3928 do kk=-1,1 3929 if (ii==0 .and. jj==0 .and. kk==0) cycle 3930 do idir=1,3 3931 pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3) 3932 end do 3933 normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) 3934 r_inscribed_sphere = min(r_inscribed_sphere, normsq) 3935 end do 3936 end do 3937 end do 3938 r_inscribed_sphere = sqrt(r_inscribed_sphere) 3939 3940 3941 !Begin the big loop on ia and ib 3942 do ia=1,natom 3943 do ib=1,natom 3944 3945 ! Simple Lattice 3946 if (abs(brav)==1) then 3947 ! In this case, it is better to work in reduced coordinates 3948 ! As rcan is in canonical coordinates, => multiplication by gprim 3949 do ii=1,3 3950 red(1,ii)= rcan(1,ia)*gprim(1,ii) +rcan(2,ia)*gprim(2,ii) +rcan(3,ia)*gprim(3,ii) 3951 red(2,ii)= rcan(1,ib)*gprim(1,ii) +rcan(2,ib)*gprim(2,ii) +rcan(3,ib)*gprim(3,ii) 3952 end do 3953 end if 3954 3955 do irpt=1,nrpt 3956 3957 ! Initialization of the weights to 1.0 3958 wghatm(ia,ib,irpt)=1.0_dp 3959 3960 ! Compute the difference vector 3961 3962 ! Simple Cubic Lattice 3963 if (abs(brav)==1) then 3964 ! Change of rpt to reduced coordinates 3965 do ii=1,3 3966 red(3,ii)= rpt(1,irpt)*gprim(1,ii) +rpt(2,irpt)*gprim(2,ii) +rpt(3,irpt)*gprim(3,ii) 3967 rdiff(ii)=red(2,ii)-red(1,ii)+red(3,ii) 3968 end do 3969 if (brav==1) then 3970 ! rdiff in cartesian coordinates 3971 do ii=1,3 3972 rdiff_tmp(ii)=rdiff(1)*rprimd(ii,1)+rdiff(2)*rprimd(ii,2)+rdiff(3)*rprimd(ii,3) 3973 end do 3974 rdiff(1:3)=rdiff_tmp(1:3) 3975 end if 3976 ! Other lattices 3977 else 3978 do ii=1,3 3979 rdiff(ii)=rcan(ii,ib)-rcan(ii,ia)+rpt(ii,irpt) 3980 end do 3981 end if 3982 3983 ! *************************************************************** 3984 3985 ! Assignement of weights 3986 3987 if(nqshft==1 .and. brav/=4)then 3988 3989 if (brav/=1) then 3990 do ii=1,3 3991 ! If the rpt vector is greater than 3992 ! the allowed space => weight = 0.0 3993 if (abs(rdiff(ii))-tol10>factor*ngqpt(ii)) then 3994 wghatm(ia,ib,irpt)=zero 3995 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 3996 ! If the point is in a boundary position => weight/2 3997 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 3998 end if 3999 end do 4000 else 4001 ! new weights 4002 wghatm(ia,ib,irpt)=zero 4003 nreq = 1 4004 do ii=1,nptws 4005 proj = rdiff(1)*ptws(1,ii)+rdiff(2)*ptws(2,ii)+rdiff(3)*ptws(3,ii) 4006 ! if rdiff closer to ptws than the origin the weight is zero 4007 ! if rdiff close to the origin with respect to all the other ptws the weight is 1 4008 ! if rdiff is equidistant from the origin and N other ptws the weight is 1/(N+1) 4009 if (proj-ptws(4,ii)>tolsym) then 4010 nreq = 0 4011 EXIT 4012 else if(abs(proj-ptws(4,ii))<=tolsym) then 4013 nreq=nreq+1 4014 end if 4015 end do 4016 if (nreq>0) then 4017 wghatm(ia,ib,irpt)=one/DBLE(nreq) 4018 end if 4019 end if 4020 4021 else if(brav==4)then 4022 ! Hexagonal 4023 4024 ! Examination of the X and Y boundaries in order to form an hexagon 4025 ! First generate the relevant boundaries 4026 rdiff(4)=0.5_dp*( rdiff(1)+sqrt(3.0_dp)*rdiff(2) ) 4027 ngqpt(4)=ngqpt(1) 4028 rdiff(5)=0.5_dp*( rdiff(1)-sqrt(3.0_dp)*rdiff(2) ) 4029 ngqpt(5)=ngqpt(1) 4030 4031 ! Test the four inequalities 4032 do ii=1,5 4033 if(ii/=2)then 4034 4035 nbord(ii)=0 4036 ! If the rpt vector is greater than 4037 ! the allowed space => weight = 0.0 4038 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4039 wghatm(ia,ib,irpt)=zero 4040 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4041 ! If the point is in a boundary position increment nbord(ii) 4042 nbord(ii)=1 4043 end if 4044 4045 end if 4046 end do 4047 4048 ! Computation of weights 4049 nbordh=nbord(1)+nbord(4)+nbord(5) 4050 if (nbordh==1) then 4051 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4052 else if (nbordh==2) then 4053 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4054 else if (nbordh/=0) then 4055 message = 'There is a problem of borders and weights (hex).' 4056 MSG_BUG(message) 4057 end if 4058 if (nbord(3)==1)then 4059 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4060 end if 4061 4062 ! BCC packing of k-points 4063 else if(nqshft==2 .and. brav/=4)then 4064 4065 ! First, generate the relevant boundaries 4066 rdiff(4)= rdiff(1)+rdiff(2) 4067 rdiff(5)= rdiff(1)-rdiff(2) 4068 rdiff(6)= rdiff(1)+rdiff(3) 4069 rdiff(7)= rdiff(1)-rdiff(3) 4070 rdiff(8)= rdiff(3)+rdiff(2) 4071 rdiff(9)= rdiff(3)-rdiff(2) 4072 if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then 4073 write(message, '(a,a,a,3i6,a,a,a,a)' )& 4074 & 'In the BCC case, the three ngqpt numbers ',ch10,& 4075 & ' ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,& 4076 & 'should be equal.',ch10,& 4077 & 'Action: use identical ngqpt(1:3) in your input file.' 4078 MSG_ERROR(message) 4079 end if 4080 do ii=4,9 4081 ngqpt(ii)=ngqpt(1) 4082 end do 4083 4084 ! Test the relevant inequalities 4085 nbord(1)=0 4086 do ii=4,9 4087 ! If the rpt vector is greater than the allowed space => weight = 0.0 4088 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4089 wghatm(ia,ib,irpt)=zero 4090 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4091 ! If the point is in a boundary position increment nbord(1) 4092 nbord(1)=nbord(1)+1 4093 end if 4094 end do 4095 4096 ! Computation of weights 4097 if (nbord(1)==1) then 4098 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4099 else if (nbord(1)==2) then 4100 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4101 else if (nbord(1)==3) then 4102 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4 4103 else if (nbord(1)==4) then 4104 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/6 4105 else if (nbord(1)/=0) then 4106 message = ' There is a problem of borders and weights (BCC).' 4107 MSG_ERROR(message) 4108 end if 4109 4110 ! FCC packing of k-points 4111 else if(nqshft==4 .and. brav/=4)then 4112 4113 ! First, generate the relevant boundaries 4114 rdiff(4)= (rdiff(1)+rdiff(2)+rdiff(3))*2._dp/3._dp 4115 rdiff(5)= (rdiff(1)-rdiff(2)+rdiff(3))*2._dp/3._dp 4116 rdiff(6)= (rdiff(1)+rdiff(2)-rdiff(3))*2._dp/3._dp 4117 rdiff(7)= (rdiff(1)-rdiff(2)-rdiff(3))*2._dp/3._dp 4118 if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then 4119 write(message, '(a,a,a,3i6,a,a,a,a)' )& 4120 & 'In the FCC case, the three ngqpt numbers ',ch10,& 4121 & ' ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,& 4122 & 'should be equal.',ch10,& 4123 & 'Action: use identical ngqpt(1:3) in your input file.' 4124 MSG_ERROR(message) 4125 end if 4126 do ii=4,7 4127 ngqpt(ii)=ngqpt(1) 4128 end do 4129 4130 ! Test the relevant inequalities 4131 nbord(1)=0 4132 do ii=1,7 4133 4134 ! If the rpt vector is greater than the allowed space => weight = 0.0 4135 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4136 wghatm(ia,ib,irpt)=zero 4137 ! If the point is in a boundary position increment nbord(1) 4138 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4139 nbord(1)=nbord(1)+1 4140 end if 4141 end do 4142 4143 ! Computation of weights 4144 if (nbord(1)==1) then 4145 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4146 else if (nbord(1)==2) then 4147 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4148 else if (nbord(1)==3) then 4149 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4 4150 else if (nbord(1)/=0 .and. wghatm(ia,ib,irpt)>1.d-10) then 4151 ! Interestingly nbord(1)==4 happens for some points outside of the volume 4152 message = ' There is a problem of borders and weights (FCC).' 4153 MSG_BUG(message) 4154 end if 4155 4156 else 4157 write(message, '(3a,i0,a)' )& 4158 & 'One should not arrive here ... ',ch10,& 4159 & 'The value nqshft ',nqshft,' is not available' 4160 MSG_BUG(message) 4161 end if 4162 end do ! Assignement of weights is done 4163 4164 end do ! End of the double loop on ia and ib 4165 end do 4166 4167 ! Check the results 4168 do ia=1,natom 4169 do ib=1,natom 4170 sumwght=zero 4171 do irpt=1,nrpt 4172 ! Check if the sum of the weights is equal to the number of q points 4173 sumwght=sumwght+wghatm(ia,ib,irpt) 4174 ! write(std_out,'(a,3i5)' )' atom1, atom2, irpt ; rpt ; wghatm ',ia,ib,irpt 4175 ! write(std_out,'(3es16.6,es18.6)' )& 4176 ! & rpt(1,irpt),rpt(2,irpt),rpt(3,irpt),wghatm(ia,ib,irpt) 4177 end do 4178 if (abs(sumwght-nqpt)>tol10) then 4179 write(message, '(a,a,a,2i4,a,a,es14.4,a,a,i4)' )& 4180 & 'The sum of the weight is not equal to nqpt.',ch10,& 4181 & 'atoms :',ia,ib,ch10,& 4182 & 'The sum of the weights is : ',sumwght,ch10,& 4183 & 'The number of q points is : ',nqpt 4184 call wrtout(std_out,message,'COLL') 4185 write(message, '(13a)')& 4186 & 'This might have several sources.',ch10,& 4187 & 'If tolsym is larger than 1.0e-8, the atom positions might be loose',ch10,& 4188 & 'and the q point weights not computed properly.',ch10,& 4189 & 'Action: make input atomic positions more symmetric.',ch10,& 4190 & 'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine, and recompile.',ch10,& 4191 & 'Actually, this can also happen when ngqpt is 0 0 0,',ch10,& 4192 & 'if abs(brav)/=1, in which case you should change brav to 1.' 4193 MSG_BUG(message) 4194 end if 4195 end do 4196 end do 4197 4198 DBG_EXIT("COLL") 4199 4200 end subroutine wght9
m_dynmat/wings3 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
wings3
FUNCTION
Suppress the wings of the cartesian 2DTE for which the diagonal element is not known
INPUTS
carflg(3,mpert,3,mpert)= ( 1 if the element of the cartesian 2DTE matrix has been calculated correctly ; 0 otherwise ) d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates mpert =maximum number of ipert
OUTPUT
d2cart(2,3,mpert,3,mpert) without the wings
PARENTS
respfn
CHILDREN
SOURCE
2383 subroutine wings3(carflg,d2cart,mpert) 2384 2385 2386 !This section has been created automatically by the script Abilint (TD). 2387 !Do not modify the following lines by hand. 2388 #undef ABI_FUNC 2389 #define ABI_FUNC 'wings3' 2390 !End of the abilint section 2391 2392 implicit none 2393 2394 !Arguments ------------------------------- 2395 !scalars 2396 integer,intent(in) :: mpert 2397 !arrays 2398 integer,intent(inout) :: carflg(3,mpert,3,mpert) 2399 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 2400 2401 !Local variables ------------------------- 2402 !scalars 2403 integer :: idir,idir1,ipert,ipert1 2404 2405 ! ********************************************************************* 2406 2407 do ipert=1,mpert 2408 do idir=1,3 2409 if(carflg(idir,ipert,idir,ipert)==0)then 2410 do ipert1=1,mpert 2411 do idir1=1,3 2412 carflg(idir,ipert,idir1,ipert1)=0 2413 carflg(idir1,ipert1,idir,ipert)=0 2414 d2cart(1,idir,ipert,idir1,ipert1)=zero 2415 d2cart(2,idir,ipert,idir1,ipert1)=zero 2416 d2cart(1,idir1,ipert1,idir,ipert)=zero 2417 d2cart(2,idir1,ipert1,idir,ipert)=zero 2418 end do 2419 end do 2420 end if 2421 end do 2422 end do 2423 2424 end subroutine wings3