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/d3lwsym
- m_dynmat/d3sym
- m_dynmat/dfpt_phfrq
- m_dynmat/dfpt_prtph
- m_dynmat/dfpt_sydy
- m_dynmat/dfpt_sygra
- 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/get_bigbox_and_weights
- m_dynmat/gtdyn9
- m_dynmat/ifclo9
- m_dynmat/make_bigbox
- m_dynmat/massmult_and_breaksym
- m_dynmat/massmult_and_breaksym_cplx
- m_dynmat/nanal9
- m_dynmat/phdispl_from_eigvec
- m_dynmat/pheigvec_normalize
- m_dynmat/q0dy3_apply
- m_dynmat/q0dy3_calc
- m_dynmat/sylwtens
- 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-2022 ABINIT group (XG, JCC, MJV, NH, RC, MVeithen, MM, MG, MT, DCA) 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!
SOURCE
19 #if defined HAVE_CONFIG_H 20 #include "config.h" 21 #endif 22 23 #include "abi_common.h" 24 25 module m_dynmat 26 27 use, intrinsic :: iso_c_binding 28 use defs_basis 29 use m_abicore 30 use m_errors 31 use m_linalg_interfaces 32 use m_xmpi 33 34 use m_fstrings, only : itoa, sjoin 35 use m_numeric_tools, only : wrap2_pmhalf, mkherm 36 use m_symtk, only : mati3inv, matr3inv, littlegroup_q 37 use m_cgtools, only : fxphas_seq 38 use m_ewald, only : ewald9 39 use m_time, only : timab 40 41 implicit none 42 43 private 44 45 public :: asria_calc ! Calculate the correction for the Acoustic sum rule on 46 ! the InterAtomic Forces or on the dynamical matrix directly 47 public :: asria_corr ! Imposition of the Acoustic sum rule on the InterAtomic Forces 48 ! or on the dynamical matrix directly from the previously calculated d2asr 49 public :: asrprs ! Imposition of the Acoustic sum rule on the InterAtomic Forces Plus Rotational Symmetry 50 public :: cart29 ! Transform a second-derivative matrix from reduced coordinates to cartesian coordinates, and also 51 ! 1) add the ionic part of the effective charges, 52 ! 2) normalize the electronic dielectric tensor, and add the vacuum polarisation 53 public :: cart39 ! Transform a vector from reduced coordinates to cartesian coordinates, 54 ! taking into account the perturbation from which it was derived, 55 ! and also check the existence of the new values. 56 public :: d2cart_to_red ! Transform a second-derivative matrix 57 ! from cartesian to reduced coordinate. 58 public :: chkph3 ! Check the completeness of the dynamical matrix 59 public :: chneu9 ! Imposition of the charge neutrality sum rule on the Effective charges 60 public :: d2sym3 ! Build (nearly) all the other matrix elements that can be build using symmetries. 61 public :: q0dy3_apply ! Takes care of the inclusion of the ewald q=0 term in the dynamical matrix 62 public :: q0dy3_calc ! Calculate the q=0 correction term to the dynamical matrix 63 ! TODO: 3 routines to symmetrize. Clarify different use cases 64 public :: symdyma ! Symmetrize the dynamical matrices 65 public :: dfpt_sygra ! Symmetrize derivatives of energy with respect to coordinates, 66 public :: dfpt_sydy ! Symmetrize dynamical matrix (eventually diagonal wrt to the atoms) 67 public :: wings3 ! Suppress the wings of the cartesian 2DTE for which the diagonal element is not known 68 public :: asrif9 ! Imposes the Acoustic Sum Rule to Interatomic Forces 69 public :: get_bigbox_and_weights ! Compute 70 public :: bigbx9 ! Generates a Big Box of R points for the Fourier Transforms the dynamical matrix 71 public :: make_bigbox ! Helper functions that faciliates the generation of a Big Box containing 72 public :: canat9 ! From reduced to canonical coordinates 73 public :: canct9 ! Convert from canonical coordinates to cartesian coordinates 74 public :: chkrp9 ! Check if the rprim used for the definition of the unit cell (in the 75 ! inputs) are consistent with the rprim used in the routine generating the Big Box 76 public :: dist9 ! Compute the distance between atoms in the big box 77 public :: ftifc_q2r ! Fourier transform of the dynamical matrices to obtain interatomic forces (real space). 78 private :: ftifc_r2q ! Fourier transform of the interatomic forces to obtain dynamical matrices (reciprocal space). 79 public :: dynmat_dq ! Compute the derivative D(q)/dq via Fourier transform of the interatomic forces 80 public :: ifclo9 ! Convert from cartesian coordinates to local coordinates 81 public :: wght9 ! Generates a weight to each R points of the Big Box and for each pair of atoms 82 public :: d3sym ! Given a set of calculated elements of the 3DTE matrix, 83 ! build (nearly) all the other matrix elements that can be build using symmetries. 84 public :: d3lwsym 85 public :: sylwtens ! Determines the set of irreductible elements of the spatial-dispersion tensors 86 public :: sytens ! Determines the set of irreductible elements of the nonlinear optical susceptibility 87 ! and Raman tensors 88 public :: axial9 ! Generates the local coordinates system from the knowledge of the first vector (longitudinal) and 89 ! the ifc matrix in cartesian coordinates 90 public :: dymfz9 ! Multiply the dynamical matrix by a phase shift to account for normalized canonical coordinates. 91 public :: nanal9 ! Subtract/Add the non-analytical part from one dynamical matrix with number iqpt. 92 public :: gtdyn9 ! Generates a dynamical matrix from interatomic force constants and 93 ! long-range electrostatic interactions. 94 public :: dfpt_phfrq ! Diagonalize IFC(q), return phonon frequencies and eigenvectors. 95 ! If q is Gamma, the non-analytical behaviour can be included. 96 public :: pheigvec_normalize ! Normalize input eigenvectors in cartesian coordinates. 97 public :: phdispl_from_eigvec ! Phonon displacements from eigenvectors 98 public :: dfpt_prtph ! Print phonon frequencies 99 public :: massmult_and_breaksym ! Multiply IFC(q) by atomic masses. 100 public :: massmult_and_breaksym_cplx ! Version for complex array 101 102 ! TODO: Change name, 103 public :: ftgam 104 public :: ftgam_init 105 106 107 ! ************************************************************************* 108 109 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.
SOURCE
133 subroutine asria_calc(asr,d2asr,d2cart,mpert,natom) 134 135 !Arguments ------------------------------- 136 !scalars 137 integer,intent(in) :: asr,mpert,natom 138 !arrays 139 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert) 140 real(dp),intent(out) :: d2asr(2,3,natom,3,natom) 141 142 !Local variables------------------------------- 143 !scalars 144 integer :: idir1,idir2,ii,ipert1,ipert2 145 integer :: constrank, imatelem, iconst, nconst, nd2_packed, info 146 !character(len=500) :: msg 147 !arrays 148 integer, allocatable :: packingindex(:,:,:,:) 149 real(dp), allocatable :: constraints(:,:,:) 150 real(dp), allocatable :: d2cart_packed(:,:) 151 real(dp), allocatable :: singvals(:) 152 real(dp), allocatable :: constr_rhs(:,:) 153 real(dp), allocatable :: work(:,:),rwork(:) 154 155 ! ********************************************************************* 156 157 d2asr = zero 158 159 if (asr==0) return 160 161 !call wrtout(std_out,' asria_calc: calculation of the correction to the ASR for the interatomic forces.') 162 do ipert1=1,natom 163 do idir1=1,3 164 do idir2=1,3 165 166 ! Compute d2asr 167 do ipert2=1,natom 168 d2asr(:,idir1,ipert1,idir2,ipert1)=& 169 & d2asr(:,idir1,ipert1,idir2,ipert1)+& 170 & d2cart(:,idir1,ipert1,idir2,ipert2) 171 end do 172 end do 173 end do 174 end do 175 176 !holistic method: overwrite d2asr with hermitian solution 177 if (asr == 5) then 178 nconst = 9*natom 179 nd2_packed = 3*natom*(3*natom+1)/2 180 ABI_MALLOC(constraints,(2,nconst, nd2_packed)) 181 ABI_MALLOC(d2cart_packed,(2,nd2_packed)) 182 ABI_MALLOC(constr_rhs,(2,nd2_packed)) 183 ABI_MALLOC(singvals,(nconst)) 184 ABI_MALLOC(work,(2,3*nd2_packed)) 185 ABI_MALLOC(rwork,(5*nd2_packed)) 186 ABI_MALLOC(packingindex,(3,natom,3,natom)) 187 ii=1 188 packingindex=-1 189 do ipert2=1,natom 190 do idir2=1,3 191 do ipert1=1,ipert2-1 192 do idir1=1,3 193 packingindex(idir1,ipert1,idir2,ipert2) = ii 194 ii = ii+1 195 end do 196 end do 197 do idir1=1,idir2 198 packingindex(idir1,ipert2,idir2,ipert2) = ii 199 ii = ii+1 200 end do 201 end do 202 end do 203 ! setup constraint matrix 204 constraints = zero 205 do ipert1=1,natom 206 do idir1=1,3 207 do idir2=1,3 208 iconst = idir2+3*(idir1-1 + 3*(ipert1-1)) 209 ! set all atom forces, this component 210 do ipert2=1,natom 211 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 212 if (imatelem == -1) then 213 imatelem = packingindex(idir2,ipert2,idir1,ipert1) 214 end if 215 constraints(1,iconst,imatelem) = one 216 end do 217 end do 218 end do 219 end do 220 221 d2cart_packed = -999.0d0 222 do ipert2=1,natom 223 do idir2=1,3 224 do ipert1=1,natom 225 do idir1=1,3 226 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 227 if (imatelem == -1) cycle 228 d2cart_packed(:,imatelem) = d2cart(:,idir1,ipert1,idir2,ipert2) 229 end do 230 end do 231 end do 232 end do 233 constr_rhs = zero 234 constr_rhs(1,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(1,:)) 235 constr_rhs(2,1:nconst) = matmul(constraints(1,:,:),d2cart_packed(2,:)) 236 237 ! lwork = 3*nd2_packed 238 call zgelss (nconst,nd2_packed,1,constraints,nconst,constr_rhs,nd2_packed,& 239 & singvals,-one,constrank,work,3*nd2_packed,rwork,info) 240 ABI_CHECK(info == 0, sjoin('zgelss returned:', itoa(info))) 241 242 ! unpack 243 do ipert2=1,natom 244 do idir2=1,3 245 do ipert1=1,natom 246 do idir1=1,3 247 imatelem = packingindex(idir1,ipert1,idir2,ipert2) 248 if (imatelem == -1) then 249 imatelem = packingindex(idir2,ipert2,idir1,ipert1) 250 ! NOTE: should complex conjugate the correction below. 251 end if 252 d2asr(:,idir1,ipert1,idir2,ipert2) = constr_rhs(:,imatelem) 253 end do 254 end do 255 end do 256 end do 257 258 ABI_FREE(constraints) 259 ABI_FREE(d2cart_packed) 260 ABI_FREE(singvals) 261 ABI_FREE(constr_rhs) 262 ABI_FREE(work) 263 ABI_FREE(rwork) 264 ABI_FREE(packingindex) 265 end if 266 267 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
SOURCE
294 subroutine asria_corr(asr,d2asr,d2cart,mpert,natom) 295 296 !Arguments ------------------------------- 297 !scalars 298 integer,intent(in) :: asr,mpert,natom 299 !arrays 300 real(dp),intent(in) :: d2asr(2,3,natom,3,natom) 301 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 302 303 !Local variables------------------------------- 304 !scalars 305 integer :: idir1,idir2,ipert1,ipert2 306 307 ! ********************************************************************* 308 309 if (asr==0) return 310 !call wrtout(std_out,' asria_corr: imposition of the ASR for the interatomic forces.') 311 312 ! Remove d2asr 313 do ipert2=1,natom 314 do idir2=1,3 315 do ipert1=1,natom 316 do idir1=1,3 317 d2cart(:,idir1,ipert1,idir2,ipert2)= d2cart(:,idir1,ipert1,idir2,ipert2) - d2asr(:,idir1,ipert1,idir2,ipert2) 318 end do 319 end do 320 end do 321 end do 322 323 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.
SOURCE
2571 subroutine asrif9(asr,atmfrc,natom,nrpt,rpt,wghatm) 2572 2573 !Arguments ------------------------------- 2574 !scalars 2575 integer,intent(in) :: asr,natom,nrpt 2576 !arrays 2577 real(dp),intent(in) :: rpt(3,nrpt),wghatm(natom,natom,nrpt) 2578 real(dp),intent(inout) :: atmfrc(3,natom,3,natom,nrpt) 2579 2580 !Local variables ------------------------- 2581 !scalars 2582 integer :: found,ia,ib,irpt,izero,mu,nu 2583 real(dp) :: sumifc 2584 2585 ! ********************************************************************* 2586 2587 if(asr==1.or.asr==2)then 2588 found=0 2589 ! Search for the R vector which is equal to ( 0 , 0 , 0 ) 2590 ! This vector leaves the atom a on itself ! 2591 do irpt=1,nrpt 2592 if (all(abs(rpt(:,irpt))<=1.0d-10)) then 2593 found=1 2594 izero=irpt 2595 end if 2596 if (found==1) exit 2597 end do 2598 2599 if(found==0)then 2600 ABI_BUG('Not able to find the vector R=(0,0,0).') 2601 end if 2602 2603 do mu=1,3 2604 do nu=1,3 2605 do ia=1,natom 2606 sumifc=zero 2607 do ib=1,natom 2608 2609 ! Get the sumifc of interatomic forces acting on the atom ia, 2610 ! either in a symmetrical manner, or an unsymmetrical one. 2611 if(asr==1)then 2612 do irpt=1,nrpt 2613 sumifc=sumifc+wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt) 2614 end do 2615 else if(asr==2)then 2616 do irpt=1,nrpt 2617 sumifc=sumifc+& 2618 (wghatm(ia,ib,irpt)*atmfrc(mu,ia,nu,ib,irpt)+& 2619 wghatm(ia,ib,irpt)*atmfrc(nu,ia,mu,ib,irpt))/2 2620 end do 2621 end if 2622 end do 2623 2624 ! Correct the self-interaction in order to fulfill the ASR 2625 atmfrc(mu,ia,nu,ia,izero)=atmfrc(mu,ia,nu,ia,izero)-sumifc 2626 if (asr==2) atmfrc(nu,ia,mu,ia,izero)=atmfrc(mu,ia,nu,ia,izero) 2627 end do 2628 end do 2629 end do 2630 end if 2631 2632 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
SOURCE
360 subroutine asrprs(asr,asrflag,rotinv,uinvers,vtinvers,singular,d2cart,mpert,natom,xcart) 361 362 !Arguments ------------------------------------ 363 !scalars 364 integer,intent(in) :: asr,asrflag,mpert,natom,rotinv 365 !arrays 366 real(dp),intent(in) :: xcart(3,natom) 367 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 368 real(dp),intent(inout) :: singular(1:3*natom*(3*natom-1)/2) 369 real(dp),intent(inout) :: uinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2) 370 real(dp),intent(inout) :: vtinvers(1:3*natom*(3*natom-1)/2,1:3*natom*(3*natom-1)/2) 371 372 !Local variables------------------------------- 373 !scalars 374 integer :: column,idir1,idir2,ii,info,ipert1,ipert2,jj,n3,row,superdim 375 real(dp) :: rcond,test 376 ! real(dp) :: tau ! tau is present but commented out in this routine 377 character(len=500) :: msg 378 !arrays 379 integer :: check(3,natom,3) 380 real(dp) :: tmp(natom,3,3),weightf(1:natom,1:natom) 381 real(dp),allocatable :: d2cartold(:,:,:,:,:),d2vecc(:),d2veccnew(:),d2vecr(:) 382 real(dp),allocatable :: d2vecrnew(:),superm(:,:),umatrix(:,:),vtmatrix(:) 383 real(dp),allocatable :: work(:) 384 385 ! ********************************************************************* 386 387 if(asr/=3 .and. asr/=4)then 388 write(msg,'(3a,i0)')& 389 'The argument asr should be 3 or 4,',ch10, 'however, asr = ',asr 390 ABI_BUG(msg) 391 end if 392 393 if (asr==3.or.asr==4)then 394 write(msg, '(a,a)' ) ch10, & 395 'asrprs: imposition of the ASR for the interatomic forces and rotational invariance' 396 call wrtout(std_out,msg) 397 end if 398 399 write(msg,'(a,i0)')' asrflag is ', asrflag 400 call wrtout([std_out, ab_out], msg) 401 402 !variables for the dimensions of the matrices 403 404 !n1=3*natom*(3*natom-1)/2 405 !n2=9*natom 406 n3=3*natom 407 408 superdim=9*natom*(natom-1)/2+n3 409 410 ABI_MALLOC(d2vecr,(1:superdim)) 411 ABI_MALLOC(d2vecc,(1:superdim)) 412 d2vecr=0d0 413 d2vecc=0d0 414 415 !should be changed set to delta function for debugging 416 weightf=1d0 417 !tau=1d-10 418 do ii=1, natom 419 ! do jj=1, ii-1 420 ! weightf(ii,jj)= & 421 ! & ((xcart(1,ii)-xcart(1,jj))**2+(xcart(2,ii)-xcart(2,jj))**2+(xcart(3,ii)-xcart(3,jj))**2)**tau 422 ! enddo 423 weightf(ii,ii)=0d0 424 end do 425 426 ABI_MALLOC(d2cartold,(2,3,mpert,3,mpert)) 427 428 d2cartold=d2cart 429 430 !setup vector with uncorrected derivatives 431 432 do ipert1=1, natom 433 do ipert2=1, ipert1-1 434 do idir1=1,3 435 do idir2=1,3 436 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 437 if(abs(d2cart(1,idir1,ipert1,idir2,ipert2))<1d-6)then 438 d2cart(1,idir1,ipert1,idir2,ipert2)=0d0 439 else 440 d2vecr(row)=4*weightf(ipert1,ipert2)*d2cart(1,idir1,ipert1,idir2,ipert2) 441 end if 442 if(abs(d2cart(2,idir1,ipert1,idir2,ipert2))<1d-6) then 443 d2cart(2,idir1,ipert1,idir2,ipert2)=0d0 444 else 445 d2vecc(row)=4*weightf(ipert1,ipert2)*d2cart(2,idir1,ipert1,idir2,ipert2) 446 end if 447 end do 448 end do 449 end do 450 end do 451 452 if(asrflag==1) then !calculate the pseudo-inverse of the supermatrix 453 ABI_MALLOC(superm,(1:superdim,1:superdim)) 454 455 superm=0d0 456 457 ! Setting up the supermatrix containing G, A, D 458 459 do ipert1=1, natom 460 do idir1=1, 3 461 ! Setting up G 462 idir2=mod(idir1,3)+1 463 row=3*(ipert1-1)+idir1 464 do ipert2=1, ipert1-1 465 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1 466 superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1) 467 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir2 468 superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2) 469 end do 470 do ipert2=ipert1+1, natom 471 column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir1-1)+rotinv 472 superm(column,row)=xcart(idir2,ipert2)-xcart(idir2,ipert1) 473 column=9*(ipert2-1)*(ipert2-2)/2+9*(ipert1-1)+3*(idir2-1)+rotinv 474 superm(column,row)=xcart(idir1,ipert1)-xcart(idir1,ipert2) 475 end do 476 end do 477 do idir1=1, 3 478 ! Setting up D 479 idir2=mod(idir1,3)+1 480 ii=mod(idir1+1,3)+1 481 do ipert2=1, ipert1-1 482 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(rotinv-1)+idir1 483 column=9*natom*(natom-1)/2+3*(ipert1-1)+idir1 484 superm(column,row)=superm(column,row)+xcart(idir2,ipert2)-xcart(idir2,ipert1) 485 column=9*natom*(natom-1)/2+3*(ipert1-1)+ii 486 superm(column,row)=superm(column,row)+xcart(ii,ipert1)-xcart(ii,ipert2) 487 row=n3+9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+rotinv 488 column=9*natom*(natom-1)/2+3*(ipert2-1)+idir1 489 superm(column,row)=superm(column,row)+xcart(idir2,ipert1)-xcart(idir2,ipert2) 490 column=9*natom*(natom-1)/2+3*(ipert2-1)+ii 491 superm(column,row)=superm(column,row)+xcart(ii,ipert2)-xcart(ii,ipert1) 492 end do 493 ! Setting up A 494 do idir2=1, 3 495 do ipert2=1, ipert1-1 496 column=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 497 row=n3+column 498 superm(column,row)=4*weightf(ipert1,ipert2) 499 end do 500 end do 501 end do 502 end do 503 504 ! calculate the pseudo-inverse of the supermatrix 505 506 ABI_MALLOC(work,(1:6*superdim)) 507 ABI_MALLOC(vtmatrix,(1:superdim)) 508 ABI_MALLOC(umatrix,(1:superdim,1:superdim)) 509 510 ! singular value decomposition of superm 511 512 call dgesvd('A','O',superdim,superdim,superm,superdim,singular,umatrix,superdim, & 513 vtmatrix, 1, work,6*superdim,info) 514 ABI_CHECK(info == 0, sjoin('dgesvd returned:', itoa(info))) 515 516 ABI_FREE(vtmatrix) 517 ABI_FREE(work) 518 519 write(msg, '(a,es16.8,es16.8)' )' Largest and smallest values from svd', singular(1), singular(superdim) 520 call wrtout([std_out, ab_out], msg) 521 522 ! Invert U and V**T, orthogonal matrices 523 524 do ii=1, superdim 525 do jj=1, superdim 526 uinvers(ii,jj)=umatrix(jj,ii) 527 vtinvers(ii,jj)=superm(jj,ii) 528 end do 529 end do 530 531 ABI_FREE(umatrix) 532 ABI_FREE(superm) 533 534 write(msg,'(a,a)')' asrprs: done with asrflag 1', ch10 535 call wrtout(std_out,msg) 536 537 end if !asrflag=1 538 539 if(asrflag==2) then 540 541 ABI_MALLOC(d2vecrnew,(1:superdim)) 542 ABI_MALLOC(d2veccnew,(1:superdim)) 543 544 ! Calculate V**T**-1 Sigma**-1 U**-1 *rhs 545 546 do ii=1, superdim 547 d2vecrnew(ii)=0d0 548 d2veccnew(ii)=0d0 549 do jj=1, superdim 550 d2vecrnew(ii)=d2vecrnew(ii)+uinvers(ii,jj)*d2vecr(jj) 551 d2veccnew(ii)=d2veccnew(ii)+uinvers(ii,jj)*d2vecc(jj) 552 end do 553 end do 554 555 rcond=1d-10*singular(1) 556 do ii=1, superdim 557 if(singular(ii)>rcond) then 558 d2vecrnew(ii)=d2vecrnew(ii)/singular(ii) 559 d2veccnew(ii)=d2veccnew(ii)/singular(ii) 560 else 561 d2vecrnew(ii)=0d0 562 d2veccnew(ii)=0d0 563 end if 564 end do 565 566 do ii=1, superdim 567 d2vecr(ii)=0d0 568 d2vecc(ii)=0d0 569 do jj=1, superdim 570 d2vecr(ii)=d2vecr(ii)+vtinvers(ii,jj)*d2vecrnew(jj) 571 d2vecc(ii)=d2vecc(ii)+vtinvers(ii,jj)*d2veccnew(jj) 572 end do 573 end do 574 575 ! Store vector back into the matrix of 2nd order derivates 576 577 do ipert1=1, natom 578 do ipert2=1, ipert1-1 579 do idir1=1,3 580 do idir2=1,3 581 row=9*(ipert1-1)*(ipert1-2)/2+9*(ipert2-1)+3*(idir1-1)+idir2 582 d2cart(1,idir1,ipert1,idir2,ipert2)=d2vecr(row) 583 d2cart(2,idir1,ipert1,idir2,ipert2)=d2vecc(row) 584 d2cart(1,idir2,ipert2,idir1,ipert1)=d2vecr(row) 585 d2cart(2,idir2,ipert2,idir1,ipert1)=d2vecc(row) 586 end do 587 end do 588 end do 589 end do 590 591 ABI_FREE(d2vecrnew) 592 ABI_FREE(d2veccnew) 593 594 check=0 595 596 do ipert1=1, natom 597 do idir1=1, 3 598 do idir2=1, 3 599 d2cart(1,idir1,ipert1,idir2,ipert1)=0d0 600 d2cart(2,idir1,ipert1,idir2,ipert1)=0d0 601 tmp(ipert1,idir1,idir2)=0d0 602 do ipert2=1, natom 603 if(ipert2/=ipert1) then 604 tmp(ipert1,idir1,idir2)=tmp(ipert1,idir1,idir2) & 605 & -d2cart(1,idir1,ipert1,idir2,ipert2) & 606 & -d2cart(1,idir2,ipert2,idir1,ipert1) 607 end if 608 end do 609 end do 610 end do 611 end do 612 613 do ipert1=1, natom 614 do idir1=1, 3 615 do idir2=1, 3 616 d2cart(1,idir1,ipert1,idir2,ipert1)=tmp(ipert1,idir1,idir2)/2 617 d2cart(1,idir2,ipert1,idir1,ipert1)=d2cart(1,idir1,ipert1,idir2,ipert1) 618 end do 619 end do 620 end do 621 622 write(std_out,*) 'this should all be zero' 623 624 do ipert1=1, natom 625 do idir1=1, 3 626 do idir2=1, 3 627 test=0d0 628 do ipert2=1, natom 629 test=test+d2cart(1,idir1,ipert1,idir2,ipert2)+d2cart(1,idir2,ipert2,idir1,ipert1) 630 end do 631 write(std_out,'(i3,i3,i3,es11.3)') idir1,ipert1,idir2,test 632 633 write(msg, '(i3,i3,i3,es11.3)' ) idir1,ipert1,idir2,test 634 call wrtout(ab_out,msg) 635 end do 636 end do 637 end do 638 639 write(std_out,*) 'these as well' 640 do ipert2=1, natom 641 do idir1=1, 3 642 do idir2=1, 3 643 test=0d0 644 do ipert1=1, natom 645 test=test+d2cart(1,idir1,ipert1,idir2,ipert2) 646 end do 647 write(std_out,'(i3,i3,i3,i3,es11.3)') idir1,ipert1,idir2,ipert2,test 648 end do 649 end do 650 end do 651 652 write(msg,'(a,a)')' asrprs: done with asrflag 2', ch10 653 call wrtout(std_out,msg) 654 655 end if !ends asrflag=2 656 657 ABI_FREE(d2vecr) 658 ABI_FREE(d2vecc) 659 ABI_FREE(d2cartold) 660 661 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
SOURCE
5244 subroutine axial9(ifccar,vect1,vect2,vect3) 5245 5246 !Arguments ------------------------------- 5247 !arrays 5248 real(dp),intent(in) :: ifccar(3,3),vect1(3) 5249 real(dp),intent(out) :: vect2(3),vect3(3) 5250 5251 !Local variables ------------------------- 5252 !scalars 5253 integer :: flag,ii,itrial,jj 5254 real(dp) :: innorm,scprod 5255 !arrays 5256 real(dp) :: work(3) 5257 5258 ! ********************************************************************* 5259 5260 do jj=1,3 5261 work(jj)=zero 5262 do ii=1,3 5263 work(jj)=work(jj)+ifccar(jj,ii)*vect1(ii) 5264 end do 5265 end do 5266 5267 flag=0 5268 do itrial=1,4 5269 scprod=zero 5270 do ii=1,3 5271 scprod=scprod+work(ii)*vect1(ii) 5272 end do 5273 5274 do ii=1,3 5275 work(ii)=work(ii)-vect1(ii)*scprod 5276 end do 5277 5278 scprod=zero 5279 do ii=1,3 5280 scprod=scprod+work(ii)**2 5281 end do 5282 5283 if(scprod<1.0d-10)then 5284 work(1:3)=zero 5285 if(itrial>1)work(itrial-1)=1.0_dp 5286 else 5287 flag=1 5288 end if 5289 5290 if(flag==1)exit 5291 end do 5292 5293 innorm=scprod**(-0.5_dp) 5294 do ii=1,3 5295 vect2(ii)=work(ii)*innorm 5296 end do 5297 5298 vect3(1)=vect1(2)*vect2(3)-vect1(3)*vect2(2) 5299 vect3(2)=vect1(3)*vect2(1)-vect1(1)*vect2(3) 5300 vect3(3)=vect1(1)*vect2(2)-vect1(2)*vect2(1) 5301 5302 end subroutine axial9
m_dynmat/bigbx9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
bigbx9
FUNCTION
Generation of a Big Box containing all the R points (cells) 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(3,nrpt)= integer coordinates of the cells (R points) in the rprim basis nprt= Number of cells (R points) in the Big Box rpt(3,mrpt)= canonical coordinates of the cells (R points) These coordinates are normalized (=> * acell(3)!!) (output only if choice=1)
SOURCE
2885 subroutine bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt) 2886 2887 !Arguments ------------------------------- 2888 !scalars 2889 integer,intent(in) :: brav,choice,mrpt,nqshft 2890 integer,intent(out) :: nrpt 2891 !arrays 2892 integer,intent(in) :: ngqpt(3) 2893 real(dp),intent(in) :: rprim(3,3) 2894 real(dp),intent(out) :: rpt(3,mrpt) 2895 integer,intent(out) :: cell(3,mrpt) 2896 2897 !Local variables ------------------------- 2898 !In some cases, the atoms coordinates are not packed in the 2899 ! [0,1]^3 cube. Then, the parameter "buffer" might be increased, 2900 !to search relevant pairs of atoms in bigger boxes than usual. 2901 !scalars 2902 integer,parameter :: buffer=1 2903 integer :: irpt,lim1,lim2,lim3,lqshft,r1,r2,r3 2904 character(len=500) :: msg 2905 2906 ! ********************************************************************* 2907 2908 lqshft=1 2909 if(nqshft/=1)lqshft=2 2910 2911 2912 !Simple Cubic Lattice 2913 if (abs(brav)==1) then 2914 lim1=((ngqpt(1))+1)*lqshft+buffer 2915 lim2=((ngqpt(2))+1)*lqshft+buffer 2916 lim3=((ngqpt(3))+1)*lqshft+buffer 2917 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1) 2918 if(choice/=0)then 2919 if (nrpt/=mrpt) then 2920 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt 2921 ABI_BUG(msg) 2922 end if 2923 irpt=0 2924 do r1=-lim1,lim1 2925 do r2=-lim2,lim2 2926 do r3=-lim3,lim3 2927 irpt=irpt+1 2928 rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3) 2929 rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3) 2930 rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3) 2931 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2932 end do 2933 end do 2934 end do 2935 end if 2936 2937 ! Face Centered Cubic Lattice 2938 else if (brav==2) then 2939 lim1=((ngqpt(1)+3)/4)*lqshft+buffer 2940 lim2=((ngqpt(2)+3)/4)*lqshft+buffer 2941 lim3=((ngqpt(3)+3)/4)*lqshft+buffer 2942 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*4 2943 if(choice/=0)then 2944 if (nrpt/=mrpt) then 2945 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt= ',mrpt 2946 ABI_BUG(msg) 2947 end if 2948 irpt=0 2949 do r1=-lim1,lim1 2950 do r2=-lim2,lim2 2951 do r3=-lim3,lim3 2952 irpt=irpt+4 2953 rpt(1,irpt-3)=r1 2954 rpt(2,irpt-3)=r2 2955 rpt(3,irpt-3)=r3 2956 rpt(1,irpt-2)=r1 2957 rpt(2,irpt-2)=r2+0.5 2958 rpt(3,irpt-2)=r3+0.5 2959 rpt(1,irpt-1)=r1+0.5 2960 rpt(2,irpt-1)=r2 2961 rpt(3,irpt-1)=r3+0.5 2962 rpt(1,irpt)=r1+0.5 2963 rpt(2,irpt)=r2+0.5 2964 rpt(3,irpt)=r3 2965 !TEST_AM 2966 ! cell(irpt-3,1)=r1;cell(irpt-3,2)=r2;cell(irpt-3,3)=r3 2967 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2968 end do 2969 end do 2970 end do 2971 end if 2972 2973 ! Body Centered Cubic Lattice 2974 else if (brav==3) then 2975 lim1=((ngqpt(1)+3)/4)*lqshft+buffer 2976 lim2=((ngqpt(2)+3)/4)*lqshft+buffer 2977 lim3=((ngqpt(3)+3)/4)*lqshft+buffer 2978 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1)*2 2979 if(choice/=0)then 2980 if(nrpt/=mrpt) then 2981 write(msg,'(2(a,i0))')' nrpt= ',nrpt,' is not equal to mrpt= ',mrpt 2982 ABI_BUG(msg) 2983 end if 2984 irpt=0 2985 do r1=-lim1,lim1 2986 do r2=-lim2,lim2 2987 do r3=-lim3,lim3 2988 irpt=irpt+2 2989 rpt(1,irpt-1)=r1 2990 rpt(2,irpt-1)=r2 2991 rpt(3,irpt-1)=r3 2992 rpt(1,irpt)=r1+0.5 2993 rpt(2,irpt)=r2+0.5 2994 rpt(3,irpt)=r3+0.5 2995 !TEST_AM 2996 ! cell(irpt-1,1)=r1;cell(irpt-1,2)=r2;cell(irpt-1,3)=r3 2997 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 2998 end do 2999 end do 3000 end do 3001 end if 3002 3003 ! Hexagonal Lattice 3004 else if (brav==4) then 3005 lim1=(ngqpt(1)+1)*lqshft+buffer 3006 lim2=(ngqpt(2)+1)*lqshft+buffer 3007 lim3=((ngqpt(3)/2)+1)*lqshft+buffer 3008 nrpt=(2*lim1+1)*(2*lim2+1)*(2*lim3+1) 3009 if(choice/=0)then 3010 if(nrpt/=mrpt)then 3011 write(msg,'(2(a,i0))')' nrpt=',nrpt,' is not equal to mrpt=',mrpt 3012 ABI_BUG(msg) 3013 end if 3014 irpt=0 3015 do r1=-lim1,lim1 3016 do r2=-lim2,lim2 3017 do r3=-lim3,lim3 3018 irpt=irpt+1 3019 rpt(1,irpt)=r1*rprim(1,1)+r2*rprim(1,2)+r3*rprim(1,3) 3020 rpt(2,irpt)=r1*rprim(2,1)+r2*rprim(2,2)+r3*rprim(2,3) 3021 rpt(3,irpt)=r1*rprim(3,1)+r2*rprim(3,2)+r3*rprim(3,3) 3022 cell(1,irpt)=r1;cell(2,irpt)=r2;cell(3,irpt)=r3 3023 end do 3024 end do 3025 end do 3026 end if 3027 3028 else 3029 write(msg,'(a,i0,a)')' The value of brav= ',brav,' is not allowed (should be -1, 1, 2 or 4).' 3030 ABI_BUG(msg) 3031 end if 3032 3033 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
SOURCE
3059 subroutine canat9(brav,natom,rcan,rprim,trans,xred) 3060 3061 !Arguments ------------------------------- 3062 !scalars 3063 integer,intent(in) :: brav,natom 3064 !arrays 3065 real(dp),intent(in) :: rprim(3,3),xred(3,natom) 3066 real(dp),intent(out) :: rcan(3,natom),trans(3,natom) 3067 3068 !Local variables ------------------------- 3069 !scalars 3070 integer :: found,iatom,ii 3071 character(len=500) :: msg 3072 !arrays 3073 real(dp) :: dontno(3,4),rec(3),rok(3),shift(3),tt(3) 3074 3075 ! ********************************************************************* 3076 3077 !Normalization of the cartesian atomic coordinates 3078 !If not normalized : rcan(i) <- rcan(i) * acell(i) 3079 do iatom=1,natom 3080 rcan(:,iatom)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3) 3081 end do 3082 3083 !Study of the different cases for the Bravais lattice: 3084 if (abs(brav)==1) then 3085 !Simple Cubic Lattice 3086 3087 do iatom=1,natom 3088 ! Canon will produces these coordinate transformations 3089 ! (Note: here we still use reduced coordinates ) 3090 call wrap2_pmhalf(xred(1,iatom),rok(1),shift(1)) 3091 call wrap2_pmhalf(xred(2,iatom),rok(2),shift(2)) 3092 call wrap2_pmhalf(xred(3,iatom),rok(3),shift(3)) 3093 3094 ! New coordinates : rcan 3095 rcan(:,iatom)=rok(1)*rprim(:,1)+rok(2)*rprim(:,2)+rok(3)*rprim(:,3) 3096 ! Translations between New and Old coordinates 3097 tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3) 3098 trans(:,iatom)=tt(:)-rcan(:,iatom) 3099 end do 3100 3101 else if (brav==2) then 3102 ! Face Centered Lattice 3103 ! Special possible translations in the F.C.C. case 3104 dontno(:,:)=zero 3105 dontno(2,2)=0.5_dp 3106 dontno(3,2)=0.5_dp 3107 dontno(1,3)=0.5_dp 3108 dontno(3,3)=0.5_dp 3109 dontno(1,4)=0.5_dp 3110 dontno(2,4)=0.5_dp 3111 do iatom=1,natom 3112 found=0 3113 do ii=1,4 3114 if (found==1) exit 3115 ! Canon will produce these coordinate transformations 3116 call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1)) 3117 call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2)) 3118 call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3)) 3119 ! In the F.C.C., ABS[ Ri ] + ABS[ Rj ] < or = 1/2 3120 ! The equal sign hase been treated using a tolerance parameter 3121 ! not to have twice the same point in the unit cell ! 3122 rok(1)=rok(1)-1.0d-10 3123 rok(2)=rok(2)-2.0d-10 3124 rok(3)=rok(3)-5.0d-10 3125 if (abs(rok(1))+abs(rok(2))<=0.5_dp) then 3126 if (abs(rok(1))+abs(rok(3))<=0.5_dp) then 3127 if (abs(rok(2))+abs(rok(3))<=0.5_dp) then 3128 tt(:)=rcan(:,iatom) 3129 ! New coordinates : rcan 3130 rcan(1,iatom)=rok(1)+1.0d-10 3131 rcan(2,iatom)=rok(2)+2.0d-10 3132 rcan(3,iatom)=rok(3)+5.0d-10 3133 ! Translations between New and Old coordinates 3134 trans(:,iatom)=tt(:)-rcan(:,iatom) 3135 found=1 3136 end if 3137 end if 3138 end if 3139 end do 3140 end do 3141 3142 else if (brav==3) then 3143 ! Body Centered Cubic Lattice 3144 ! Special possible translations in the B.C.C. case 3145 dontno(:,1)=zero 3146 dontno(:,2)=0.5_dp 3147 do iatom=1,natom 3148 found=0 3149 do ii=1,2 3150 if (found==1) exit 3151 ! Canon will produce these coordinate transformations 3152 call wrap2_pmhalf(rcan(1,iatom)+dontno(1,ii),rok(1),shift(1)) 3153 call wrap2_pmhalf(rcan(2,iatom)+dontno(2,ii),rok(2),shift(2)) 3154 call wrap2_pmhalf(rcan(3,iatom)+dontno(3,ii),rok(3),shift(3)) 3155 ! In the F.C.C., ABS[ Ri ] < or = 1/2 3156 ! and ABS[ R1 ] + ABS[ R2 ] + ABS[ R3 ] < or = 3/4 3157 ! The equal signs have been treated using a tolerance parameter 3158 ! not to have twice the same point in the unit cell ! 3159 rok(1)=rok(1)-1.0d-10 3160 rok(2)=rok(2)-2.0d-10 3161 rok(3)=rok(3)-5.0d-10 3162 if(abs(rok(1))+abs(rok(2))+abs(rok(3))<=0.75_dp)then 3163 if ( abs(rok(1))<=0.5_dp .and. abs(rok(2))<=0.5_dp .and. abs(rok(3))<=0.5_dp) then 3164 tt(:)=rcan(:,iatom) 3165 ! New coordinates : rcan 3166 rcan(1,iatom)=rok(1)+1.0d-10 3167 rcan(2,iatom)=rok(2)+2.0d-10 3168 rcan(3,iatom)=rok(3)+5.0d-10 3169 ! Translations between New and Old coordinates 3170 trans(:,iatom)=tt(:)-rcan(:,iatom) 3171 found=1 3172 end if 3173 end if 3174 end do 3175 end do 3176 3177 else if (brav==4) then 3178 ! Hexagonal Lattice 3179 ! In this case, it is easier first to work in reduced coordinates space ! 3180 do iatom=1,natom 3181 ! Passage from the reduced space to the "lozenge" cell 3182 rec(1)=xred(1,iatom)-0.5_dp 3183 rec(2)=xred(2,iatom)-0.5_dp 3184 rec(3)=xred(3,iatom) 3185 ! Canon will produces these coordinate transformations 3186 call wrap2_pmhalf(rec(1),rok(1),shift(1)) 3187 call wrap2_pmhalf(rec(2),rok(2),shift(2)) 3188 call wrap2_pmhalf(rec(3),rok(3),shift(3)) 3189 rec(1)=rok(1)+0.5_dp 3190 rec(2)=rok(2)+0.5_dp 3191 rec(3)=rok(3) 3192 ! Passage in Cartesian Normalized Coordinates 3193 rcan(:,iatom)=rec(1)*rprim(:,1)+rec(2)*rprim(:,2)+rec(3)*rprim(:,3) 3194 ! Use of a tolerance parameter not to have twice the same point in the unit cell ! 3195 rcan(1,iatom)=rcan(1,iatom)-1.0d-10 3196 rcan(2,iatom)=rcan(2,iatom)-2.0d-10 3197 ! Passage to the honeycomb hexagonal unit cell ! 3198 if (rcan(1,iatom)>0.5_dp) then 3199 rcan(1,iatom)=rcan(1,iatom)-1.0_dp 3200 end if 3201 if (rcan(1,iatom)>zero.and.rcan(1,iatom)+sqrt(3.0_dp)*rcan(2,iatom)>1.0_dp) then 3202 rcan(1,iatom)=rcan(1,iatom)-0.5_dp 3203 rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp 3204 end if 3205 if (rcan(1,iatom)<=zero.and.sqrt(3.0_dp)*rcan(2,iatom)-rcan(1,iatom)>1.0_dp) then 3206 rcan(1,iatom)=rcan(1,iatom)+0.5_dp 3207 rcan(2,iatom)=rcan(2,iatom)-sqrt(3.0_dp)*0.5_dp 3208 end if 3209 ! Translations between New and Old coordinates 3210 tt(:)=xred(1,iatom)*rprim(:,1)+xred(2,iatom)*rprim(:,2)+xred(3,iatom)*rprim(:,3) 3211 trans(:,iatom)=tt(:)-rcan(:,iatom) 3212 end do 3213 3214 ! End of the possible cases for brav : -1, 1, 2, 4. 3215 else 3216 write(msg, '(a,i0,a,a,a)' )& 3217 'The required value of brav=',brav,' is not available.',ch10,& 3218 'It should be -1, 1,2 or 4 .' 3219 ABI_BUG(msg) 3220 end if 3221 3222 call wrtout(std_out,' Canonical Atomic Coordinates ') 3223 do iatom=1,natom 3224 write(msg, '(a,i5,3es18.8)' )' atom',iatom,rcan(1,iatom),rcan(2,iatom),rcan(3,iatom) 3225 call wrtout(std_out,msg) 3226 end do 3227 3228 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.
SOURCE
3259 subroutine canct9(acell,gprim,ib,index,irpt,natom,nrpt,rcan,rcart,rprim,rpt) 3260 3261 !Arguments ------------------------------- 3262 !scalars 3263 integer,intent(in) :: index,natom,nrpt 3264 integer,intent(out) :: ib,irpt 3265 !arrays 3266 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3) 3267 real(dp),intent(in) :: rpt(3,nrpt) 3268 real(dp),intent(out) :: rcart(3) 3269 3270 !Local variables ------------------------- 3271 !scalars 3272 integer :: jj 3273 !arrays 3274 real(dp) :: xred(3) 3275 3276 ! ********************************************************************* 3277 3278 irpt=(index-1)/natom+1 3279 ib=index-natom*(irpt-1) 3280 3281 !Transform the canonical coordinates to reduced coord. 3282 do jj=1,3 3283 xred(jj)=gprim(1,jj)*(rpt(1,irpt)+rcan(1,ib))& 3284 & +gprim(2,jj)*(rpt(2,irpt)+rcan(2,ib))& 3285 & +gprim(3,jj)*(rpt(3,irpt)+rcan(3,ib)) 3286 end do 3287 3288 !Then to cartesian coordinates (here the position of the atom b) 3289 do jj=1,3 3290 rcart(jj)=xred(1)*acell(1)*rprim(jj,1)+& 3291 & xred(2)*acell(2)*rprim(jj,2)+& 3292 & xred(3)*acell(3)*rprim(jj,3) 3293 end do 3294 3295 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
SOURCE
703 subroutine cart29(blkflg,blkval,carflg,d2cart,& 704 & gprimd,iblok,mpert,natom,nblok,ntypat,rprimd,typat,ucvol,zion) 705 706 !Arguments ------------------------------- 707 !scalars 708 integer,intent(in) :: iblok,mpert,natom,nblok,ntypat 709 real(dp),intent(in) :: ucvol 710 !arrays 711 integer,intent(in) :: blkflg(3,mpert,3,mpert,nblok),typat(natom) 712 integer,intent(out) :: carflg(3,mpert,3,mpert) 713 real(dp),intent(in) :: blkval(2,3,mpert,3,mpert,nblok),gprimd(3,3),rprimd(3,3) 714 real(dp),intent(in) :: zion(ntypat) 715 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 716 717 !Local variables ------------------------- 718 !scalars 719 integer :: idir1,idir2,ii,ipert1,ipert2 720 !arrays 721 integer :: flg1(3),flg2(3) 722 real(dp) :: vec1(3),vec2(3) 723 724 ! ********************************************************************* 725 726 !First, copy the data blok in place. 727 d2cart(:,:,:,:,:)=blkval(:,:,:,:,:,iblok) 728 729 !Cartesian coordinates transformation (in two steps) 730 !First step 731 do ipert1=1,mpert 732 do ipert2=1,mpert 733 do ii=1,2 734 do idir1=1,3 735 do idir2=1,3 736 vec1(idir2)=d2cart(ii,idir1,ipert1,idir2,ipert2) 737 ! Note here blkflg 738 flg1(idir2)=blkflg(idir1,ipert1,idir2,ipert2,iblok) 739 end do 740 call cart39(flg1,flg2,gprimd,ipert2,natom,rprimd,vec1,vec2) 741 do idir2=1,3 742 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) 743 ! And here carflg 744 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir2) 745 end do 746 end do 747 end do 748 end do 749 end do 750 751 !Second step 752 do ipert1=1,mpert 753 do ipert2=1,mpert 754 do ii=1,2 755 do idir2=1,3 756 do idir1=1,3 757 vec1(idir1)=d2cart(ii,idir1,ipert1,idir2,ipert2) 758 ! Note here carflg 759 flg1(idir1)=carflg(idir1,ipert1,idir2,ipert2) 760 end do 761 call cart39(flg1,flg2,gprimd,ipert1,natom,rprimd,vec1,vec2) 762 do idir1=1,3 763 d2cart(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) 764 ! And here carflg again 765 carflg(idir1,ipert1,idir2,ipert2)=flg2(idir1) 766 end do 767 end do 768 end do 769 end do 770 end do 771 772 !For the dielectric tensor, takes into account the volume 773 !of the unit cell, and add the unit matrix (polarization of the vacuum) 774 do idir1=1,3 775 do idir2=1,3 776 do ii=1,2 777 d2cart(ii,idir1,natom+2,idir2,natom+2)=& 778 & -four_pi/ucvol*d2cart(ii,idir1,natom+2,idir2,natom+2) 779 end do 780 end do 781 end do 782 783 do idir1=1,3 784 d2cart(1,idir1,natom+2,idir1,natom+2)=& 785 & 1.0_dp+d2cart(1,idir1,natom+2,idir1,natom+2) 786 end do 787 788 !Add the ionic charges to delta z to get the effective charges 789 do ipert1=1,natom 790 do idir1=1,3 791 d2cart(1,idir1,ipert1,idir1,natom+2)=& 792 & zion(typat(ipert1))+d2cart(1,idir1,ipert1,idir1,natom+2) 793 end do 794 end do 795 do ipert2=1,natom 796 do idir2=1,3 797 d2cart(1,idir2,natom+2,idir2,ipert2)=& 798 & zion(typat(ipert2))+d2cart(1,idir2,natom+2,idir2,ipert2) 799 end do 800 end do 801 802 !For the piezoelectric tensor, takes into account the volume of the unit cell 803 do ipert2=natom+3,natom+4 804 do idir1=1,3 805 do idir2=1,3 806 do ii=1,2 807 d2cart(ii,idir1,natom+2,idir2,ipert2)=& 808 & (1.0_dp/ucvol)*d2cart(ii,idir1,natom+2,idir2,ipert2) 809 d2cart(ii,idir2,ipert2,idir1,natom+2)=& 810 & (1.0_dp/ucvol)*d2cart(ii,idir2,ipert2,idir1,natom+2) 811 end do 812 end do 813 end do 814 end do 815 816 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
SOURCE
845 subroutine cart39(flg1,flg2,gprimd,ipert,natom,rprimd,vec1,vec2) 846 847 !Arguments ------------------------------- 848 !scalars 849 integer,intent(in) :: ipert,natom 850 !arrays 851 integer,intent(in) :: flg1(3) 852 integer,intent(out) :: flg2(3) 853 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3),vec1(3) 854 real(dp),intent(out) :: vec2(3) 855 856 !Local variables ------------------------- 857 !scalars 858 integer :: idir,ii 859 860 ! ********************************************************************* 861 862 !Treat phonon-type perturbation 863 if(ipert>=1.and.ipert<=natom)then 864 865 do idir=1,3 866 vec2(idir)=zero 867 flg2(idir)=1 868 do ii=1,3 869 if(abs(gprimd(idir,ii))>1.0d-10)then 870 if(flg1(ii)==1)then 871 vec2(idir)=vec2(idir)+gprimd(idir,ii)*vec1(ii) 872 else 873 flg2(idir)=0 874 end if 875 end if 876 end do 877 if(flg2(idir)==0)vec2(idir)=zero 878 end do 879 880 ! Treat electric field and qvec perturbations 881 else if(ipert==natom+2.or.ipert==natom+8) then 882 ! OCL SCALAR 883 do idir=1,3 884 vec2(idir)=zero 885 flg2(idir)=1 886 ! OCL SCALAR 887 do ii=1,3 888 if(abs(rprimd(idir,ii))>1.0d-10)then 889 if(flg1(ii)==1)then 890 vec2(idir)=vec2(idir)+rprimd(idir,ii)*vec1(ii)/two_pi 891 else 892 flg2(idir)=0 893 end if 894 end if 895 end do 896 if(flg2(idir)==0)vec2(idir)=zero 897 end do 898 899 ! Treat other perturbations 900 else 901 do idir=1,3 902 vec2(idir)=vec1(idir) 903 flg2(idir)=flg1(idir) 904 end do 905 end if 906 907 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
SOURCE
1081 subroutine chkph3(carflg,idir,mpert,natom) 1082 1083 !Arguments ------------------------------- 1084 !scalars 1085 integer,intent(in) :: idir,mpert,natom 1086 !arrays 1087 integer,intent(in) :: carflg(3,mpert,3,mpert) 1088 1089 !Local variables ------------------------- 1090 !scalars 1091 integer :: idir1,idir2,ipert1,ipert2,send 1092 character(len=500) :: msg 1093 1094 ! ********************************************************************* 1095 1096 send=0 1097 1098 !Check the elements of the analytical part of the dynamical matrix 1099 do ipert2=1,natom 1100 do idir2=1,3 1101 do ipert1=1,natom 1102 do idir1=1,3 1103 if(carflg(idir1,ipert1,idir2,ipert2)==0)then 1104 send=1 1105 end if 1106 end do 1107 end do 1108 end do 1109 end do 1110 1111 !If some electric field is present 1112 if(idir/=0)then 1113 1114 ! Check the dielectric constant 1115 if(carflg(idir,natom+2,idir,natom+2)==0)then 1116 send=1 1117 end if 1118 1119 ! Check the effective charges 1120 do ipert1=1,natom 1121 do idir1=1,3 1122 if(carflg(idir1,ipert1,idir,natom+2)==0)then 1123 send=1 1124 end if 1125 end do 1126 end do 1127 1128 end if 1129 1130 ! If needed, send the message 1131 if(send==1)then 1132 write(msg, '(a,a,a,a)' )& 1133 ' chkph3 : WARNING -',ch10,& 1134 ' Dynamical matrix incomplete, phonon frequencies may be wrong, check input variables rfatpol and rfdir.' 1135 call wrtout([std_out, ab_out],msg) 1136 end if 1137 1138 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)
SOURCE
3319 subroutine chkrp9(brav,rprim) 3320 3321 !Arguments ------------------------------- 3322 !scalars 3323 integer,intent(in) :: brav 3324 !arrays 3325 real(dp),intent(in) :: rprim(3,3) 3326 3327 !Local variables ------------------------- 3328 !scalars 3329 integer :: ii,jj 3330 character(len=500) :: msg 3331 3332 ! ********************************************************************* 3333 3334 if (abs(brav)==1) then 3335 ! Simple Cubic Lattice No condition in this case ! 3336 continue 3337 3338 else if (brav==2) then 3339 ! Face Centered Lattice 3340 do ii=1,3 3341 do jj=1,3 3342 if ( ( ii==jj .and. abs(rprim(ii,jj))>tol10) .or. (ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then 3343 write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3344 'The input variable rprim does not correspond to the',ch10,& 3345 'fixed rprim to be used with brav=2 and ifcflag=1 :',ch10,& 3346 ' 0 1/2 1/2',ch10,& 3347 ' 1/2 0 1/2',ch10,& 3348 ' 1/2 1/2 0 ',ch10,& 3349 'Action: rebuild your DDB by using the latter rprim.' 3350 ABI_ERROR(msg) 3351 end if 3352 end do 3353 end do 3354 3355 else if (brav==3) then 3356 ! Body Centered Cubic Lattice 3357 do ii=1,3 3358 do jj=1,3 3359 if ( ( ii==jj .and. abs(rprim(ii,jj)+.5_dp)>tol10) .or. (ii/=jj .and. abs(rprim(ii,jj)-.5_dp)>tol10) ) then 3360 write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3361 'The input variable rprim does not correspond to the',ch10,& 3362 'fixed rprim to be used with brav=3 and ifcflag=1 :',ch10,& 3363 ' -1/2 1/2 1/2',ch10,& 3364 ' 1/2 -1/2 1/2',ch10,& 3365 ' 1/2 1/2 -1/2',ch10,& 3366 'Action: rebuild your DDB by using the latter rprim.' 3367 ABI_ERROR(msg) 3368 end if 3369 end do 3370 end do 3371 3372 else if (brav==4) then 3373 ! Hexagonal Lattice 3374 if (abs(rprim(1,1)-1.0_dp)>tol10 .or. & 3375 abs(rprim(3,3)-1.0_dp)>tol10 .or. & 3376 abs(rprim(2,1) )>tol10 .or. & 3377 abs(rprim(3,1) )>tol10 .or. & 3378 abs(rprim(1,3) )>tol10 .or. & 3379 abs(rprim(2,3) )>tol10 .or. & 3380 abs(rprim(3,2) )>tol10 .or. & 3381 abs(rprim(1,2)+0.5_dp)>tol10 .or. & 3382 abs(rprim(2,2)-0.5_dp*sqrt(3.0_dp))>tol10 ) then 3383 write(msg, '(a,a,a,a,a,a,a,a,a,a,a)' )& 3384 'The input variable rprim does not correspond to the',ch10,& 3385 'fixed rprim to be used with brav=4 and ifcflag=1 :',ch10,& 3386 ' 1 0 0',ch10,& 3387 ' -1/2 sqrt[3]/2 0',ch10,& 3388 ' 0 0 1',ch10,& 3389 'Action: rebuild your DDB by using the latter rprim.' 3390 ABI_ERROR(msg) 3391 end if 3392 3393 else 3394 write(msg, '(a,i4,a,a,a,a,a)' )& 3395 'The value of brav=',brav,' is not allowed.',ch10,& 3396 'Only -1, 1,2,3 or 4 are allowed.',ch10,& 3397 'Action: change the value of brav in your input file.' 3398 ABI_ERROR(msg) 3399 end if 3400 3401 end subroutine chkrp9
m_dynmat/chneu9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
chneu9
FUNCTION
Imposition of the charge neutrality 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
SOURCE
1168 subroutine chneu9(chneut,d2cart,mpert,natom,ntypat,selectz,typat,zion) 1169 1170 !Arguments ------------------------------- 1171 !scalars 1172 integer,intent(in) :: chneut,mpert,natom,ntypat,selectz 1173 !arrays 1174 integer,intent(in) :: typat(natom) 1175 real(dp),intent(in) :: zion(ntypat) 1176 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 1177 1178 !Local variables ------------------------- 1179 !scalars 1180 integer :: idir1,idir2,ii,ipert1,ipert2 1181 character(len=500) :: msg 1182 !arrays 1183 real(dp) :: sumwght(2) 1184 real(dp),allocatable :: wghtat(:) 1185 1186 ! ********************************************************************* 1187 1188 ABI_MALLOC(wghtat,(natom)) 1189 1190 !In case of acoustic sum rule imposition, compute the weights on each atom. 1191 if (chneut==1)then 1192 1193 ! The weight is the same for all atom 1194 do ipert1=1,natom 1195 wghtat(ipert1)=1./natom 1196 end do 1197 1198 else if (chneut==2) then 1199 1200 ! The weight is proportional to the diagonal electronic screening charge of the atom 1201 sumwght(1)=zero 1202 do ipert1=1,natom 1203 wghtat(ipert1)=zero 1204 do idir1=1,3 1205 wghtat(ipert1)=wghtat(ipert1)+& 1206 & d2cart(1,idir1,ipert1,idir1,natom+2)+& 1207 & d2cart(1,idir1,natom+2,idir1,ipert1)-2*zion(typat(ipert1)) 1208 end do 1209 sumwght(1)=sumwght(1)+wghtat(ipert1) 1210 end do 1211 1212 ! Normalize the weights to unity 1213 do ipert1=1,natom 1214 wghtat(ipert1)=wghtat(ipert1)/sumwght(1) 1215 end do 1216 end if 1217 1218 !Calculation of the violation of the charge neutrality 1219 !and imposition of the charge neutrality condition 1220 if (chneut/=0)then 1221 write(msg, '(a,a,a,a,a,a,a)' )& 1222 ' The violation of the charge neutrality conditions',ch10,& 1223 ' by the effective charges is as follows :',ch10,& 1224 ' atom electric field',ch10,& 1225 ' displacement direction ' 1226 call wrtout(ab_out,msg) 1227 do idir1=1,3 1228 do idir2=1,3 1229 do ii=1,2 1230 sumwght(ii)=zero 1231 do ipert1=1,natom 1232 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir2,natom+2) 1233 end do 1234 do ipert1=1,natom 1235 d2cart(ii,idir1,ipert1,idir2,natom+2)=& 1236 d2cart(ii,idir1,ipert1,idir2,natom+2)-sumwght(ii)*wghtat(ipert1) 1237 end do 1238 end do 1239 write(msg, '(i8,i16,2f16.6)' ) idir1,idir2,sumwght(1),sumwght(2) 1240 call wrtout(ab_out,msg) 1241 end do 1242 end do 1243 write(msg, '(a)' )' ' 1244 call wrtout(ab_out,msg) 1245 1246 ! The same for the symmetrical part 1247 do idir1=1,3 1248 do idir2=1,3 1249 do ii=1,2 1250 sumwght(ii)=zero 1251 do ipert2=1,natom 1252 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir2,ipert2) 1253 end do 1254 do ipert2=1,natom 1255 d2cart(ii,idir1,natom+2,idir2,ipert2)=& 1256 d2cart(ii,idir1,natom+2,idir2,ipert2)-sumwght(ii)*wghtat(ipert2) 1257 end do 1258 end do 1259 end do 1260 end do 1261 end if 1262 1263 !Selection of the trace of the effective charge tensor attached to each atom 1264 if(selectz==1)then 1265 do ipert1=1,natom 1266 do ii=1,2 1267 sumwght(ii)=zero 1268 do idir1=1,3 1269 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,ipert1,idir1,natom+2) 1270 end do 1271 do idir1=1,3 1272 do idir2=1,3 1273 d2cart(ii,idir1,ipert1,idir2,natom+2)=zero 1274 end do 1275 end do 1276 do idir1=1,3 1277 d2cart(ii,idir1,ipert1,idir1,natom+2)=sumwght(ii)/3.0_dp 1278 end do 1279 end do 1280 end do 1281 ! Do the same for the symmetrical part of d2cart 1282 do ipert2=1,natom 1283 do ii=1,2 1284 sumwght(ii)=zero 1285 do idir1=1,3 1286 sumwght(ii)=sumwght(ii)+d2cart(ii,idir1,natom+2,idir1,ipert2) 1287 end do 1288 do idir1=1,3 1289 do idir2=1,3 1290 d2cart(ii,idir1,natom+2,idir2,ipert2)=zero 1291 end do 1292 end do 1293 do idir1=1,3 1294 d2cart(ii,idir1,natom+2,idir1,ipert2)=sumwght(ii)/3.0_dp 1295 end do 1296 end do 1297 end do 1298 end if 1299 1300 !Selection of the symmetric part of the effective charge tensor attached to each atom 1301 if(selectz==2)then 1302 do ipert1=1,natom 1303 do ii=1,2 1304 do idir1=1,3 1305 do idir2=1,3 1306 sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)& 1307 & +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp 1308 d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii) 1309 d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii) 1310 end do 1311 end do 1312 end do 1313 end do 1314 ! Do the same for the symmetrical part of d2cart 1315 do ipert1=1,natom 1316 do ii=1,2 1317 do idir1=1,3 1318 do idir2=1,3 1319 sumwght(ii)=(d2cart(ii,idir1,ipert1,idir2,natom+2)& 1320 & +d2cart(ii,idir2,ipert1,idir1,natom+2))/2.0_dp 1321 d2cart(ii,idir1,ipert1,idir2,natom+2)=sumwght(ii) 1322 d2cart(ii,idir2,ipert1,idir1,natom+2)=sumwght(ii) 1323 end do 1324 end do 1325 end do 1326 end do 1327 end if 1328 1329 !Write the effective charge tensor 1330 write(msg, '(a,a,a,a,a,a,a)' )& 1331 ' Effective charge tensors after ',ch10,& 1332 ' imposition of the charge neutrality,',ch10,& 1333 ' and eventual restriction to some part :',ch10,& 1334 ' atom displacement ' 1335 call wrtout(ab_out,msg) 1336 1337 do ipert1=1,natom 1338 do idir1=1,3 1339 write(msg, '(2i10,3es16.6)' )ipert1,idir1,(d2cart(1,idir1,ipert1,idir2,natom+2),idir2=1,3) 1340 call wrtout(ab_out,msg) 1341 end do 1342 end do 1343 1344 !Zero the imaginary part of the dynamical matrix 1345 write(msg, '(a)' )' Now, the imaginary part of the dynamical matrix is zeroed ' 1346 call wrtout(ab_out,msg) 1347 call wrtout(std_out,msg) 1348 1349 do ipert1=1,natom 1350 do ipert2=1,natom 1351 do idir1=1,3 1352 do idir2=1,3 1353 d2cart(2,idir1,ipert1,idir2,ipert2)=zero 1354 end do 1355 end do 1356 end do 1357 end do 1358 1359 ABI_FREE(wghtat) 1360 1361 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
SOURCE
939 subroutine d2cart_to_red(d2cart, d2red, gprimd, rprimd, mpert, natom, & 940 & ntypat,typat,ucvol,zion) 941 942 !Arguments ------------------------------- 943 !scalars 944 integer,intent(in) :: mpert,natom,ntypat 945 real(dp),intent(in) :: ucvol 946 !arrays 947 integer,intent(in) :: typat(natom) 948 real(dp),intent(in) :: d2cart(2,3,mpert,3,mpert) 949 real(dp),intent(out) :: d2red(2,3,mpert,3,mpert) 950 real(dp),intent(in) :: gprimd(3,3),rprimd(3,3) 951 real(dp),intent(in) :: zion(ntypat) 952 953 !Local variables ------------------------- 954 !scalars 955 integer :: idir1,idir2,ii,ipert1,ipert2 956 real(dp) :: fac 957 !arrays 958 integer :: flg1(3),flg2(3) 959 real(dp) :: vec1(3),vec2(3) 960 961 ! ********************************************************************* 962 963 flg1 = one 964 flg2 = one 965 966 d2red = d2cart 967 968 !Remove the ionic charges to z to get the change in effective charges 969 do ipert1=1,natom 970 do idir1=1,3 971 d2red(1,idir1,ipert1,idir1,natom+2)=& 972 & d2red(1,idir1,ipert1,idir1,natom+2) - zion(typat(ipert1)) 973 end do 974 end do 975 do ipert2=1,natom 976 do idir2=1,3 977 d2red(1,idir2,natom+2,idir2,ipert2)=& 978 & d2red(1,idir2,natom+2,idir2,ipert2) - zion(typat(ipert2)) 979 end do 980 end do 981 982 ! Remove the vacuum polarizability from the dielectric tensor 983 do idir1=1,3 984 d2red(1,idir1,natom+2,idir1,natom+2)=& 985 & d2red(1,idir1,natom+2,idir1,natom+2) - 1.0_dp 986 end do 987 988 ! Scale the dielectric tensor with the volue of the unit cell 989 do idir1=1,3 990 do idir2=1,3 991 do ii=1,2 992 d2red(ii,idir1,natom+2,idir2,natom+2)=& 993 & - (ucvol / four_pi) * d2red(ii,idir1,natom+2,idir2,natom+2) 994 end do 995 end do 996 end do 997 998 !For the piezoelectric tensor, takes into account the volume of the unit cell 999 do ipert2=natom+3,natom+4 1000 do idir1=1,3 1001 do idir2=1,3 1002 do ii=1,2 1003 d2red(ii,idir1,natom+2,idir2,ipert2)=& 1004 & (ucvol)*d2red(ii,idir1,natom+2,idir2,ipert2) 1005 d2red(ii,idir2,ipert2,idir1,natom+2)=& 1006 & (ucvol)*d2red(ii,idir2,ipert2,idir1,natom+2) 1007 end do 1008 end do 1009 end do 1010 end do 1011 1012 ! Reduced coordinates transformation (in two steps) 1013 ! Note that rprimd and gprimd are swapped, compared to what cart39 expects 1014 ! A factor of (2pi) ** 2 is added to transform the electric field perturbations 1015 1016 !First step 1017 do ipert1=1,mpert 1018 fac = one; if (ipert1==natom+2) fac = two_pi ** 2 1019 1020 do ipert2=1,mpert 1021 do ii=1,2 1022 do idir1=1,3 1023 do idir2=1,3 1024 vec1(idir2)=d2red(ii,idir1,ipert1,idir2,ipert2) 1025 end do 1026 ! Transform vector from cartesian to reduced coordinates 1027 call cart39(flg1,flg2,transpose(rprimd),ipert1,natom,transpose(gprimd),vec1,vec2) 1028 do idir2=1,3 1029 d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir2) * fac 1030 end do 1031 end do 1032 end do 1033 end do 1034 end do 1035 1036 !Second step 1037 do ipert1=1,mpert 1038 do ipert2=1,mpert 1039 fac = one; if (ipert2==natom+2) fac = two_pi ** 2 1040 1041 do ii=1,2 1042 do idir2=1,3 1043 do idir1=1,3 1044 vec1(idir1)=d2red(ii,idir1,ipert1,idir2,ipert2) 1045 end do 1046 ! Transform vector from cartesian to reduced coordinates 1047 call cart39(flg1,flg2,transpose(rprimd),ipert2,natom,transpose(gprimd),vec1,vec2) 1048 do idir1=1,3 1049 d2red(ii,idir1,ipert1,idir2,ipert2)=vec2(idir1) * fac 1050 end do 1051 end do 1052 end do 1053 end do 1054 end do 1055 1056 1057 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 zero_by_symm= if 1, set blkflg to 1 for the elements that must be zero by symmetry, and zero them. This has the indirect effect of being able to resymmetrize the whole matrix, thus enforcing better the symmetry for the 2DTE.
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.
SOURCE
1424 subroutine d2sym3(blkflg,d2,indsym,mpert,natom,nsym,qpt,symq,symrec,symrel,timrev,zero_by_symm) 1425 1426 !Arguments ------------------------------- 1427 !scalars 1428 integer,intent(in) :: mpert,natom,nsym,timrev,zero_by_symm 1429 !arrays 1430 integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym) 1431 integer,intent(in),target :: symrec(3,3,nsym),symrel(3,3,nsym) 1432 integer,intent(inout) :: blkflg(3,mpert,3,mpert) 1433 real(dp),intent(in) :: qpt(3) 1434 real(dp),intent(inout) :: d2(2,3,mpert,3,mpert) 1435 1436 !Local variables ------------------------- 1437 !scalars 1438 logical, parameter :: do_final_sym=.true. 1439 logical :: qzero 1440 integer :: exch12,found,idir1,idir2,idisy1,idisy2,ipert1,ipert2 1441 integer :: ipesy1,ipesy2,isgn,isym,ithree,itirev,nblkflg_is_one,noccur,nsym_used,quit,quit1 1442 real(dp) :: arg1,arg2,im,norm,re,sumi,sumr,xi,xr 1443 !arrays 1444 integer,pointer :: sym1_(:,:,:),sym2_(:,:,:) 1445 real(dp),allocatable :: d2tmp1(:,:,:),d2tmp2(:,:,:),d2work(:,:,:,:,:) 1446 1447 ! ********************************************************************* 1448 1449 qzero=(qpt(1)**2+qpt(2)**2+qpt(3)**2<tol16) 1450 1451 !Here look after exchange of 1 and 2 axis, 1452 !for electric field in diamond symmetry 1453 exch12=0 1454 if (qzero) then 1455 do isym=1,nsym 1456 exch12=1 1457 if(symrel(1,1,isym)/=0)exch12=0 1458 if(symrel(1,2,isym)/=1)exch12=0 1459 if(symrel(1,3,isym)/=0)exch12=0 1460 if(symrel(2,1,isym)/=1)exch12=0 1461 if(symrel(2,2,isym)/=0)exch12=0 1462 if(symrel(2,3,isym)/=0)exch12=0 1463 if(symrel(3,1,isym)/=0)exch12=0 1464 if(symrel(3,2,isym)/=0)exch12=0 1465 if(symrel(3,3,isym)/=1)exch12=0 1466 ! if(exch12==1) write(std_out,*)' d2sym3 : found exchange 1 2 =',isym 1467 if(exch12==1)exit 1468 end do 1469 end if 1470 1471 !Exchange of perturbations 1472 1473 !Consider two cases : either time-reversal symmetry 1474 !conserves the wavevector, or not 1475 if(timrev==0)then 1476 1477 ! do ipert1=1,mpert See notes 1478 do ipert1=1,min(natom+2,mpert) 1479 do idir1=1,3 1480 1481 ! Since the matrix is hermitian, the diagonal elements are real 1482 d2(2,idir1,ipert1,idir1,ipert1)=zero 1483 1484 ! do ipert2=1,mpert See notes 1485 do ipert2=1,min(natom+2,mpert) 1486 do idir2=1,3 1487 1488 ! If an element exists 1489 if(blkflg(idir1,ipert1,idir2,ipert2)==1)then 1490 1491 ! Either complete the symmetric missing element 1492 if(blkflg(idir2,ipert2,idir1,ipert1)==0)then 1493 1494 d2(1,idir2,ipert2,idir1,ipert1)= d2(1,idir1,ipert1,idir2,ipert2) 1495 d2(2,idir2,ipert2,idir1,ipert1)=-d2(2,idir1,ipert1,idir2,ipert2) 1496 1497 blkflg(idir2,ipert2,idir1,ipert1)=1 1498 1499 ! Or symmetrize (the matrix is hermitian) in case both exists 1500 ! (Note : this opportunity has been disabled for more 1501 ! obvious search for bugs in the code ) 1502 ! else 1503 ! sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2) 1504 ! sumi=d2(1,idir2,ipert2,idir1,ipert1)-d2(1,idir1,ipert1,idir2,ipert2) 1505 ! d2(1,idir2,ipert2,idir1,ipert1)=half*sumr 1506 ! d2(1,idir1,ipert1,idir2,ipert2)=half*sumr 1507 ! d2(2,idir2,ipert2,idir1,ipert1)=half*sumi 1508 ! d2(2,idir1,ipert1,idir2,ipert2)=-half*sumi 1509 end if 1510 end if 1511 1512 end do 1513 end do 1514 1515 end do 1516 end do 1517 1518 ! Here, case with time-reversal symmetry 1519 else 1520 1521 ! do ipert1=1,mpert See notes 1522 do ipert1=1,min(natom+2,mpert) 1523 do idir1=1,3 1524 ! do ipert2=1,mpert See notes 1525 do ipert2=1,min(natom+2,mpert) 1526 do idir2=1,3 1527 d2(2,idir1,ipert1,idir2,ipert2)=zero 1528 1529 ! If an element exists 1530 if(blkflg(idir1,ipert1,idir2,ipert2)==1)then 1531 1532 ! Either complete the symmetric missing element 1533 if(blkflg(idir2,ipert2,idir1,ipert1)==0)then 1534 1535 d2(1,idir2,ipert2,idir1,ipert1)=d2(1,idir1,ipert1,idir2,ipert2) 1536 blkflg(idir2,ipert2,idir1,ipert1)=1 1537 1538 ! Or symmetrize (the matrix is hermitian) in case both exists 1539 ! (Note : this opportunity has been disabled for more 1540 ! obvious search for bugs in the code ) 1541 ! else 1542 ! sumr=d2(1,idir2,ipert2,idir1,ipert1)+d2(1,idir1,ipert1,idir2,ipert2) 1543 ! d2(1,idir2,ipert2,idir1,ipert1)=half*sumr 1544 ! d2(1,idir1,ipert1,idir2,ipert2)=half*sumr 1545 end if 1546 1547 end if 1548 end do 1549 end do 1550 end do 1551 end do 1552 end if 1553 1554 !Big Big Loop : symmetrize three times, because 1555 !of some cases in which one element is not yet available 1556 !at the first pass, and even at the second one ! 1557 do ithree=1,3 1558 1559 ! Big loop on all elements 1560 ! do ipert1=1,mpert See notes 1561 do ipert1=1,min(natom+2,mpert) 1562 1563 ! Select the symmetries according to pertubation 1 1564 if (ipert1<=natom)then 1565 sym1_ => symrec 1566 else 1567 sym1_ => symrel 1568 end if 1569 1570 do idir1=1,3 1571 ! do ipert2=1,mpert See notes 1572 do ipert2=1,min(natom+2,mpert) 1573 1574 ! Select the symmetries according to pertubation 2 1575 if (ipert2<=natom)then 1576 sym2_ => symrec 1577 else 1578 sym2_ => symrel 1579 end if 1580 1581 do idir2=1,3 1582 1583 ! Will get element (idir1,ipert1,idir2,ipert2) 1584 ! so this element should not yet be present ... 1585 if(blkflg(idir1,ipert1,idir2,ipert2)/=1)then 1586 1587 d2(1,idir1,ipert1,idir2,ipert2)=zero 1588 d2(2,idir1,ipert1,idir2,ipert2)=zero 1589 1590 ! Loop on all symmetries, including time-reversal 1591 quit1=0 1592 do isym=1,nsym 1593 do itirev=1,2 1594 isgn=3-2*itirev 1595 1596 if(symq(4,itirev,isym)/=0)then 1597 found=1 1598 1599 ! Here select the symmetric of ipert1 1600 if(ipert1<=natom)then 1601 ipesy1=indsym(4,isym,ipert1) 1602 else if(ipert1==(natom+2).and.qzero)then 1603 ipesy1=ipert1 1604 else 1605 found=0 1606 end if 1607 1608 ! Here select the symmetric of ipert2 1609 if(ipert2<=natom)then 1610 ipesy2=indsym(4,isym,ipert2) 1611 else if(ipert2==(natom+2).and.qzero)then 1612 ipesy2=ipert2 1613 else 1614 found=0 1615 end if 1616 1617 ! Now that a symmetric perturbation has been obtained, 1618 ! including the expression of the symmetry matrix, see 1619 ! if the symmetric values are available 1620 if( found==1 ) then 1621 1622 sumr=zero 1623 sumi=zero 1624 noccur=0 1625 nblkflg_is_one=0 1626 quit=0 1627 do idisy1=1,3 1628 do idisy2=1,3 1629 if(sym1_(idir1,idisy1,isym)/=0 .and. sym2_(idir2,idisy2,isym)/=0 )then 1630 if(blkflg(idisy1,ipesy1,idisy2,ipesy2)==1)then 1631 nblkflg_is_one=nblkflg_is_one+1 1632 sumr=sumr+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*& 1633 & d2(1,idisy1,ipesy1,idisy2,ipesy2) 1634 sumi=sumi+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym)*& 1635 & d2(2,idisy1,ipesy1,idisy2,ipesy2) 1636 1637 ! Here, in case the symmetric of the element 1638 ! is the element, or the symmetric with 1639 ! respect to permutation of perturbations 1640 ! (some more conditions on the time-reversal 1641 ! symmetry must be fulfilled although) 1642 else if( idisy1==idir1 .and. ipesy1==ipert1& 1643 & .and. idisy2==idir2 .and. ipesy2==ipert2& 1644 & .and.(isgn==1 .or. timrev==1 & 1645 & .or. (idir1==idir2 .and. ipert1==ipert2)))& 1646 & then 1647 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1648 else if( idisy1==idir2 .and. ipesy1==ipert2& 1649 & .and. idisy2==idir1 .and. ipesy2==ipert1& 1650 & .and.(isgn==-1 .or. timrev==1& 1651 & .or. (idir1==idir2 .and. ipert1==ipert2)))& 1652 & then 1653 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1654 1655 ! Here, electric field case 1656 else if( exch12==1 .and. & 1657 & ipert1==natom+2 .and. ipert2==natom+2& 1658 & .and.(( idisy1+idir1 ==3 & 1659 & .and. idisy2==3 .and. idir2==3)& 1660 & .or. ( idisy1+idir2 ==3& 1661 & .and. idisy2==3 .and. idir1==3)& 1662 & .or. ( idisy2+idir2 ==3& 1663 & .and. idisy1==3 .and. idir1==3)& 1664 & .or. ( idisy2+idir1 ==3& 1665 & .and. idisy1==3 .and. idir2==3)))& 1666 & then 1667 noccur=noccur+sym1_(idir1,idisy1,isym)*sym2_(idir2,idisy2,isym) 1668 1669 else 1670 ! Not found 1671 found=0 1672 quit=1 1673 exit 1674 end if 1675 1676 end if 1677 end do 1678 if(quit==1)exit 1679 end do 1680 end if 1681 1682 ! In case zero_by_symm==0, the computed materix element must be associated to at least one really computed matrix element 1683 if(zero_by_symm==0 .and. nblkflg_is_one==0)then 1684 found=0 1685 endif 1686 1687 ! Now, if still found and associated to at least one really computed matrix element, put the correct value into array d2 1688 if(found==1)then 1689 1690 ! In case of phonons, need to take into account the 1691 ! time-reversal symmetry, and the shift back to the unit cell 1692 ! 1693 ! XG990712 : I am not sure this must be kept for electric field ... 1694 ! 1) Consider time-reversal symmetry 1695 sumi=isgn*sumi 1696 1697 if(ipert1<=natom .and. ipert2<=natom)then 1698 ! 2) Shift the atoms back to the unit cell. 1699 arg1=two_pi*( qpt(1)*indsym(1,isym,ipert1)& 1700 & +qpt(2)*indsym(2,isym,ipert1)& 1701 & +qpt(3)*indsym(3,isym,ipert1) ) 1702 arg2=two_pi*( qpt(1)*indsym(1,isym,ipert2)& 1703 & +qpt(2)*indsym(2,isym,ipert2)& 1704 & +qpt(3)*indsym(3,isym,ipert2) ) 1705 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 1706 ! XG010117 Must use isgn 1707 im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)) 1708 else 1709 re=one 1710 im=zero 1711 end if 1712 1713 ! Final check, could still fail if the 1714 ! element was its own symmetric 1715 if (abs(1.0_dp-re*noccur)< 1.0d-6.and.abs(im*noccur) <1.0d-6) then 1716 found=0 1717 end if 1718 1719 end if 1720 1721 if(found==1)then 1722 1723 if(noccur==0)then 1724 d2(1,idir1,ipert1,idir2,ipert2)=re*sumr-im*sumi 1725 d2(2,idir1,ipert1,idir2,ipert2)=re*sumi+im*sumr 1726 else 1727 ! See page July 2, 1994 in computer codes notebook 1728 xr=re*sumr-im*sumi 1729 xi=re*sumi+im*sumr 1730 norm=one+noccur**2-two*re*noccur 1731 xr=xr/norm 1732 xi=xi/norm 1733 d2(1,idir1,ipert1,idir2,ipert2)=& 1734 & (one-re*noccur)*xr-im*noccur*xi 1735 d2(2,idir1,ipert1,idir2,ipert2)=& 1736 & (one-re*noccur)*xi+im*noccur*xr 1737 end if 1738 1739 ! The element has been constructed ! 1740 blkflg(idir1,ipert1,idir2,ipert2)=1 1741 1742 quit1=1 1743 exit ! Exit loop on symmetry operations 1744 end if 1745 1746 ! End loop on all symmetries + time-reversal 1747 end if 1748 end do 1749 if(quit1==1)exit 1750 end do 1751 1752 end if 1753 end do ! End big loop on all elements 1754 end do 1755 end do 1756 end do 1757 1758 end do ! End Big Big Loop 1759 1760 !MT oct. 20, 2014: 1761 !Once the matrix has been built, it does not necessarily fulfill the correct symmetries. 1762 !It has just been filled up from rows or columns that only fulfill symmetries preserving 1763 !one particular perturbation. 1764 !An additional symmetrization might solve this (do not consider TR-symmetry) 1765 if (do_final_sym) then 1766 ABI_MALLOC(d2tmp1,(2,3,3)) 1767 ABI_MALLOC(d2tmp2,(2,3,3)) 1768 ABI_MALLOC(d2work,(2,3,mpert,3,mpert)) 1769 d2work(:,:,:,:,:)=d2(:,:,:,:,:) 1770 do ipert1=1,min(natom+2,mpert) 1771 if ((ipert1==natom+1.or.ipert1==natom+10.or.ipert1==natom+11).or.(ipert1==natom+2.and.(.not.qzero))) cycle 1772 if (ipert1<=natom)then 1773 sym1_ => symrec 1774 else 1775 sym1_ => symrel 1776 end if 1777 do ipert2=1,min(natom+2,mpert) 1778 ! if (any(blkflg(:,ipert1,:,ipert2)==0)) cycle 1779 if ((ipert2==natom+1.or.ipert2==natom+10.or.ipert2==natom+11).or.(ipert2==natom+2.and.(.not.qzero))) cycle 1780 if (ipert2<=natom)then 1781 sym2_ => symrec 1782 else 1783 sym2_ => symrel 1784 end if 1785 nsym_used=0 1786 d2tmp2(:,:,:)=zero 1787 do isym=1,nsym 1788 if (symq(4,1,isym)==1) then 1789 ipesy1=ipert1;if (ipert1<=natom) ipesy1=indsym(4,isym,ipert1) 1790 ipesy2=ipert2;if (ipert2<=natom) ipesy2=indsym(4,isym,ipert2) 1791 ! The condition on next line is too severe, since some elements of sym1_ or sym2_ might be zero, 1792 ! which means not all blkflg(:,ipesy1,:,ipesy2) would need to be 1 to symmetrize the matrix. 1793 ! However, coding something more refined is really more difficult. 1794 ! This condition then has the side effect that more symmetries can be applied when zero_by_symm==1, 1795 ! since blkflg can be set to 1 when the symmetries guarantee the matrix element to be zero. 1796 if (all(blkflg(:,ipesy1,:,ipesy2)==1)) then 1797 nsym_used=nsym_used+1 1798 re=one;im=zero 1799 if (ipert1<=natom.and.ipert2<=natom.and.(.not.qzero)) then 1800 arg1=two_pi*(qpt(1)*(indsym(1,isym,ipert1)-indsym(1,isym,ipert2)) & 1801 & +qpt(2)*(indsym(2,isym,ipert1)-indsym(2,isym,ipert2)) & 1802 & +qpt(3)*(indsym(3,isym,ipert1)-indsym(3,isym,ipert2))) 1803 re=cos(arg1);im=sin(arg1) 1804 end if 1805 d2tmp1(:,:,:)=zero 1806 do idir2=1,3 !kappa 1807 do idir1=1,3 !mu 1808 do idisy1=1,3 !nu 1809 d2tmp1(:,idir1,idir2)=d2tmp1(:,idir1,idir2) & 1810 & +sym1_(idir1,idisy1,isym)*d2(:,idisy1,ipesy1,idir2,ipesy2) 1811 end do 1812 end do 1813 end do 1814 do idir2=1,3 !mu 1815 do idir1=1,3 !kappa 1816 do idisy2=1,3 !nu 1817 d2tmp2(1,idir1,idir2)=d2tmp2(1,idir1,idir2) & 1818 & +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*re-d2tmp1(2,idir1,idisy2)*im) 1819 d2tmp2(2,idir1,idir2)=d2tmp2(2,idir1,idir2) & 1820 & +sym2_(idir2,idisy2,isym)*(d2tmp1(1,idir1,idisy2)*im+d2tmp1(2,idir1,idisy2)*re) 1821 end do 1822 end do 1823 end do 1824 end if 1825 end if 1826 end do ! isym 1827 if (nsym_used>0) d2work(:,1:3,ipert1,1:3,ipert2)=d2tmp2(:,1:3,1:3)/dble(nsym_used) 1828 end do !ipert2 1829 end do !ipert1 1830 if (mpert>=natom) d2(:,1:3,1:natom,1:3,1:natom)=d2work(:,1:3,1:natom,1:3,1:natom) 1831 if (mpert>=natom+2) then 1832 d2(:,1:3,natom+2,1:3,1:natom)=d2work(:,1:3,natom+2,1:3,1:natom) 1833 d2(:,1:3,1:natom,1:3,natom+2)=d2work(:,1:3,1:natom,1:3,natom+2) 1834 d2(:,1:3,natom+2,1:3,natom+2)=d2work(:,1:3,natom+2,1:3,natom+2) 1835 end if 1836 ABI_FREE(d2tmp1) 1837 ABI_FREE(d2tmp2) 1838 ABI_FREE(d2work) 1839 end if 1840 1841 end subroutine d2sym3
m_dynmat/d3lwsym [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
d3lwsym
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
has_strain = if .true. i2pert includes strain perturbation 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 reduced space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space) symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian 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,m_nonlinear
CHILDREN
SOURCE
4466 !subroutine d3lwsym(blkflg,d3,has_strain,indsym,mpert,natom,nsym,symrec,symrel,symrel_cart) 4467 subroutine d3lwsym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel) 4468 4469 !Arguments ------------------------------- 4470 !scalars 4471 integer,intent(in) :: mpert,natom,nsym 4472 ! logical,intent(in) :: has_strain 4473 !arrays 4474 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4475 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 4476 real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert) 4477 ! real(dp),intent(in) :: symrel_cart(3,3,nsym) 4478 4479 !Local variables ------------------------- 4480 !scalars 4481 integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3 4482 integer :: ipesy1,ipesy2,ipesy3,isym,ithree 4483 !integer :: istr,i2dir_a,i2dir_b,disy2_a,idisy2_b 4484 real(dp) :: sumi,sumr 4485 logical :: is_strain 4486 !arrays 4487 ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 4488 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 4489 ! integer :: strflg(3,mpert,3,3,3,mpert),strflg_car(3,mpert,3,3,3,mpert) 4490 ! real(dp) :: d3str(2,3,mpert,3,3,3,mpert) 4491 4492 ! ********************************************************************* 4493 4494 !First, take into account the permutations symmetry of 4495 !(i1pert,i1dir) and (i2pert,i2dir) 4496 do i1pert = 1, mpert 4497 do i2pert = 1, mpert 4498 do i3pert = 1, mpert 4499 4500 do i1dir = 1, 3 4501 do i2dir = 1, 3 4502 do i3dir = 1, 3 4503 4504 if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. & 4505 (blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert)/=1)) then 4506 4507 d3(1,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = & 4508 d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4509 d3(2,i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = & 4510 -d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4511 4512 blkflg(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = 1 4513 4514 end if 4515 4516 end do 4517 end do 4518 end do 4519 4520 end do 4521 end do 4522 end do 4523 4524 !For strain perturbation we need an array with the two strain indexes 4525 ! if (has_strain) then 4526 ! strflg=0 4527 ! d3str=zero 4528 ! do i3pert=1, mpert 4529 ! do i3dir=1,3 4530 ! do i2pert=natom+3,natom+4 4531 ! do i2dir=1,3 4532 ! if (i2pert==natom+3) istr=i2dir 4533 ! if (i2pert==natom+4) istr=3+i2dir 4534 ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr) 4535 ! do i1pert=1,mpert 4536 ! do i1dir=1,3 4537 ! if (blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1) then 4538 ! strflg(i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)=1 4539 ! d3str(:,i1dir,i1pert,i2dir_a,i2dir_b,i3dir,i3pert)= & 4540 ! & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4541 ! if (i2pert==natom+4) then 4542 ! strflg(i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)=1 4543 ! d3str(:,i1dir,i1pert,i2dir_b,i2dir_a,i3dir,i3pert)= & 4544 ! & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4545 ! end if 4546 ! end if 4547 ! end do 4548 ! end do 4549 ! end do 4550 ! end do 4551 ! end do 4552 ! end do 4553 ! end if 4554 4555 !Big Big Loop : symmetrize three times, because 4556 !of some cases in which one element is not yet available 4557 !at the first pass, and even at the second one ! 4558 4559 do ithree=1,3 4560 4561 ! Loop over perturbations 4562 do i1pert = 1, mpert 4563 do i2pert = 1, mpert 4564 is_strain=.false. 4565 do i3pert = 1, mpert 4566 4567 do i1dir = 1, 3 4568 do i2dir = 1, 3 4569 do i3dir = 1, 3 4570 4571 ! Will get element (idir1,ipert1,idir2,ipert2) 4572 ! so this element should not yet be present ... 4573 if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then 4574 4575 d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp 4576 4577 do isym = 1, nsym 4578 4579 found = 1 4580 4581 if (i1pert <= natom) then 4582 ipesy1 = indsym(4,isym,i1pert) 4583 sym1(:,:) = symrec(:,:,isym) 4584 else if (i1pert == natom + 2) then 4585 ipesy1 = i1pert 4586 sym1(:,:) = symrel(:,:,isym) 4587 else 4588 found = 0 4589 end if 4590 4591 if (i2pert <= natom) then 4592 ipesy2 = indsym(4,isym,i2pert) 4593 sym2(:,:) = symrec(:,:,isym) 4594 else if (i2pert == natom + 2) then 4595 ipesy2 = i2pert 4596 sym2(:,:) = symrel(:,:,isym) 4597 else if (i2pert == natom + 3.or. i2pert == natom + 4) then 4598 !TODO: Symmetries on strain perturbation do not work yet. 4599 found = 0 4600 is_strain=.true. 4601 4602 ! ipesy2 = i2pert 4603 ! sym2(:,:) = NINT(symrel_cart(:,:,isym)) 4604 ! if (i2pert==natom+3) istr=i2dir 4605 ! if (i2pert==natom+4) istr=3+i2dir 4606 ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr) 4607 else 4608 found = 0 4609 end if 4610 4611 if (i3pert <= natom) then 4612 ipesy3 = indsym(4,isym,i3pert) 4613 sym3(:,:) = symrec(:,:,isym) 4614 else if (i3pert == natom + 2.or.i3pert == natom + 8) then 4615 ipesy3 = i3pert 4616 sym3(:,:) = symrel(:,:,isym) 4617 else 4618 found = 0 4619 end if 4620 4621 sumr = 0_dp ; sumi = 0_dp; 4622 if (.not.is_strain) then 4623 do idisy1 = 1, 3 4624 do idisy2 = 1, 3 4625 do idisy3 = 1, 3 4626 4627 if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) & 4628 & .and.(sym3(i3dir,idisy3) /=0)) then 4629 4630 if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then 4631 4632 sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4633 & sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4634 sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4635 & sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4636 4637 else 4638 4639 found = 0 4640 4641 end if 4642 4643 end if 4644 4645 end do 4646 end do 4647 end do 4648 else 4649 ! do idisy1 = 1, 3 4650 ! !do idisy2_a = 1, 3 4651 ! ! do idisy2_b = 1, 3 4652 ! do idisy2 = 1, 3 4653 ! if (ipesy2==natom+3) istr=idisy2 4654 ! if (ipesy2==natom+4) istr=3+idisy2 4655 ! idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr) 4656 ! do idisy3 = 1, 3 4657 ! 4658 ! if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir_a,idisy2_a) /=0) & 4659 !& .and.(sym2(i2dir_b,idisy2_b) /=0).and.(sym3(i3dir,idisy3) /=0)) then 4660 ! 4661 ! if (strflg(idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3) == 1) then 4662 ! 4663 ! sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_a)* & 4664 !& sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*& 4665 !& d3str(1,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3) 4666 ! sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir_a,idisy2_b)*& 4667 !& sym2(i2dir_b,idisy2_b)*sym3(i3dir,idisy3)*& 4668 !& d3str(2,idisy1,ipesy1,idisy2_a,idisy2_b,idisy3,ipesy3) 4669 ! 4670 ! else 4671 ! 4672 ! found = 0 4673 ! 4674 ! end if 4675 ! 4676 ! end if 4677 ! 4678 ! end do 4679 ! !end do 4680 ! end do 4681 ! end do 4682 end if 4683 4684 if (found == 1) then 4685 d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr 4686 d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 4687 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 4688 end if 4689 4690 end do ! isym 4691 4692 end if ! blkflg 4693 4694 ! Close loop over perturbations 4695 end do 4696 end do 4697 end do 4698 end do 4699 end do 4700 end do 4701 4702 end do ! close loop over ithree 4703 4704 end subroutine d3lwsym
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
SOURCE
4266 subroutine d3sym(blkflg,d3,indsym,mpert,natom,nsym,symrec,symrel) 4267 4268 !Arguments ------------------------------- 4269 !scalars 4270 integer,intent(in) :: mpert,natom,nsym 4271 !arrays 4272 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4273 integer,intent(inout) :: blkflg(3,mpert,3,mpert,3,mpert) 4274 real(dp),intent(inout) :: d3(2,3,mpert,3,mpert,3,mpert) 4275 4276 !Local variables ------------------------- 4277 !scalars 4278 integer :: found,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert,idisy1,idisy2,idisy3 4279 integer :: ipesy1,ipesy2,ipesy3,isym,ithree 4280 real(dp) :: sumi,sumr 4281 !arrays 4282 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 4283 4284 ! ********************************************************************* 4285 4286 !DEBUG 4287 !write(std_out,*)'d3sym : enter' 4288 !do i1dir = 1, 3 4289 !do i2dir = 1, 3 4290 !do i3dir = 1, 3 4291 !write(std_out,*)i1dir,i2dir,i3dir,blkflg(i1dir,natom+2,i2dir,natom+2,i3dir,natom+2) 4292 !end do 4293 !end do 4294 !end do 4295 !stop 4296 !ENDDEBUG 4297 4298 !First, take into account the permutations symmetry of 4299 !(i1pert,i1dir) and (i3pert,i3dir) 4300 do i1pert = 1, mpert 4301 do i2pert = 1, mpert 4302 do i3pert = 1, mpert 4303 4304 do i1dir = 1, 3 4305 do i2dir = 1, 3 4306 do i3dir = 1, 3 4307 4308 if ((blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)==1).and. & 4309 & (blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert)/=1)) then 4310 4311 d3(:,i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = & 4312 & d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) 4313 4314 blkflg(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = 1 4315 4316 end if 4317 4318 end do 4319 end do 4320 end do 4321 4322 end do 4323 end do 4324 end do 4325 4326 !Big Big Loop : symmetrize three times, because 4327 !of some cases in which one element is not yet available 4328 !at the first pass, and even at the second one ! 4329 4330 do ithree=1,3 4331 4332 ! Loop over perturbations 4333 do i1pert = 1, mpert 4334 do i2pert = 1, mpert 4335 do i3pert = 1, mpert 4336 4337 do i1dir = 1, 3 4338 do i2dir = 1, 3 4339 do i3dir = 1, 3 4340 4341 ! Will get element (idir1,ipert1,idir2,ipert2) 4342 ! so this element should not yet be present ... 4343 if(blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert)/=1)then 4344 4345 d3(:,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 0_dp 4346 4347 do isym = 1, nsym 4348 4349 found = 1 4350 4351 if (i1pert <= natom) then 4352 ipesy1 = indsym(4,isym,i1pert) 4353 sym1(:,:) = symrec(:,:,isym) 4354 else if (i1pert == natom + 2) then 4355 ipesy1 = i1pert 4356 sym1(:,:) = symrel(:,:,isym) 4357 else 4358 found = 0 4359 end if 4360 4361 if (i2pert <= natom) then 4362 ipesy2 = indsym(4,isym,i2pert) 4363 sym2(:,:) = symrec(:,:,isym) 4364 else if (i2pert == natom + 2) then 4365 ipesy2 = i2pert 4366 sym2(:,:) = symrel(:,:,isym) 4367 else 4368 found = 0 4369 end if 4370 4371 if (i3pert <= natom) then 4372 ipesy3 = indsym(4,isym,i3pert) 4373 sym3(:,:) = symrec(:,:,isym) 4374 else if (i3pert == natom + 2) then 4375 ipesy3 = i3pert 4376 sym3(:,:) = symrel(:,:,isym) 4377 else 4378 found = 0 4379 end if 4380 4381 sumr = 0_dp ; sumi = 0_dp; 4382 do idisy1 = 1, 3 4383 do idisy2 = 1, 3 4384 do idisy3 = 1, 3 4385 4386 if ((sym1(i1dir,idisy1) /=0).and.(sym2(i2dir,idisy2) /=0) & 4387 & .and.(sym3(i3dir,idisy3) /=0)) then 4388 4389 if (blkflg(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 1) then 4390 4391 sumr = sumr + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4392 & sym3(i3dir,idisy3)*d3(1,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4393 sumi = sumi + sym1(i1dir,idisy1)*sym2(i2dir,idisy2)*& 4394 & sym3(i3dir,idisy3)*d3(2,idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) 4395 4396 else 4397 4398 found = 0 4399 4400 end if 4401 4402 end if 4403 4404 end do 4405 end do 4406 end do 4407 4408 if (found == 1) then 4409 d3(1,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumr 4410 d3(2,i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = sumi 4411 blkflg(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 4412 end if 4413 4414 end do ! isym 4415 4416 end if ! blkflg 4417 4418 ! Close loop over perturbations 4419 end do 4420 end do 4421 end do 4422 end do 4423 end do 4424 end do 4425 4426 end do ! close loop over ithree 4427 4428 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.
SOURCE
5660 subroutine dfpt_phfrq(amu,displ,d2cart,eigval,eigvec,indsym,& 5661 & mpert,msym,natom,nsym,ntypat,phfrq,qphnrm,qphon,rprimd,& 5662 & symdynmat,symrel,symafm,typat,ucvol) 5663 5664 !Arguments ------------------------------- 5665 !scalars 5666 integer,intent(in) :: mpert,msym,natom,nsym,ntypat,symdynmat 5667 real(dp),intent(in) :: qphnrm,ucvol 5668 !arrays 5669 integer,intent(in) :: indsym(4,msym*natom),symrel(3,3,nsym),typat(natom) 5670 integer,intent(in) :: symafm(nsym) 5671 real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),rprimd(3,3) 5672 real(dp),intent(inout) :: qphon(3) 5673 real(dp),intent(out) :: displ(2*3*natom*3*natom),eigval(3*natom) 5674 real(dp),intent(out) :: eigvec(2*3*natom*3*natom),phfrq(3*natom) 5675 5676 !Local variables ------------------------- 5677 !scalars 5678 integer :: analyt,i1,i2,idir1,idir2,ier,ii,imode,ipert1,ipert2 5679 integer :: jmode,indexi,indexj,index 5680 real(dp) :: epsq,qphon2 5681 logical,parameter :: debug = .False. 5682 real(dp) :: sc_prod 5683 !arrays 5684 real(dp) :: qptn(3),dum(2,0) !, tsec(2) 5685 real(dp),allocatable :: matrx(:,:),zeff(:,:),zhpev1(:,:),zhpev2(:) 5686 5687 ! ********************************************************************* 5688 5689 ! Keep track of time spent in dfpt_phfrq 5690 !call timab(1751, 1, tsec) 5691 5692 ! Prepare the diagonalisation: analytical part. 5693 ! Note: displ is used as work space here 5694 i1=0 5695 do ipert1=1,natom 5696 do idir1=1,3 5697 i1=i1+1 5698 i2=0 5699 do ipert2=1,natom 5700 do idir2=1,3 5701 i2=i2+1 5702 index=i1+3*natom*(i2-1) 5703 displ(2*index-1)=d2cart(1,idir1,ipert1,idir2,ipert2) 5704 displ(2*index )=d2cart(2,idir1,ipert1,idir2,ipert2) 5705 end do 5706 end do 5707 end do 5708 end do 5709 5710 ! Determine the analyticity of the matrix. 5711 analyt=1; if(abs(qphnrm)<tol8) analyt=0 5712 if(abs(qphon(1))<tol8.and.abs(qphon(2))<tol8.and.abs(qphon(3))<tol8) analyt=2 5713 5714 ! In case of q=Gamma, only the real part is used 5715 if(analyt==0 .or. analyt==2)then 5716 do i1=1,3*natom 5717 do i2=1,3*natom 5718 index=i1+3*natom*(i2-1) 5719 displ(2*index)=zero 5720 end do 5721 end do 5722 end if 5723 5724 ! In the case the non-analyticity is required: 5725 ! the tensor is in cartesian coordinates and this means that qphon must be in given in Cartesian coordinates. 5726 if(analyt==0)then 5727 5728 ! Normalize the limiting direction 5729 qphon2=qphon(1)**2+qphon(2)**2+qphon(3)**2 5730 qphon(:)=qphon(:)/sqrt(qphon2) 5731 5732 ! Get the dielectric constant for the limiting direction 5733 epsq=zero 5734 do idir1=1,3 5735 do idir2=1,3 5736 epsq= epsq + qphon(idir1)*qphon(idir2) * d2cart(1,idir1,natom+2,idir2,natom+2) 5737 end do 5738 end do 5739 5740 ABI_MALLOC(zeff,(3,natom)) 5741 5742 ! Get the effective charges for the limiting direction 5743 do idir1=1,3 5744 do ipert1=1,natom 5745 zeff(idir1,ipert1)=zero 5746 do idir2=1,3 5747 zeff(idir1,ipert1) = zeff(idir1,ipert1) + qphon(idir2)* d2cart(1,idir1,ipert1,idir2,natom+2) 5748 end do 5749 end do 5750 end do 5751 5752 ! Get the non-analytical part of the dynamical matrix, and suppress its imaginary part. 5753 i1=0 5754 do ipert1=1,natom 5755 do idir1=1,3 5756 i1=i1+1 5757 i2=0 5758 do ipert2=1,natom 5759 do idir2=1,3 5760 i2=i2+1 5761 index=i1+3*natom*(i2-1) 5762 displ(2*index-1)=displ(2*index-1)+four_pi/ucvol*zeff(idir1,ipert1)*zeff(idir2,ipert2)/epsq 5763 displ(2*index )=zero 5764 end do 5765 end do 5766 end do 5767 end do 5768 5769 ABI_FREE(zeff) 5770 end if ! End of the non-analyticity treatment 5771 5772 ! Multiply IFC(q) by masses 5773 call massmult_and_breaksym(natom, ntypat, typat, amu, displ) 5774 5775 ! *********************************************************************** 5776 ! Diagonalize the dynamical matrix 5777 5778 !Symmetrize the dynamical matrix 5779 !FIXME: swap the next 2 lines and update test files to include symmetrization 5780 ! for Gamma point too (except in non-analytic case) 5781 !if (symdynmat==1 .and. analyt > 0) then 5782 if (symdynmat==1 .and. analyt == 1) then 5783 qptn(:)=qphon(:) 5784 if (analyt==1) qptn(:)=qphon(:)/qphnrm 5785 call symdyma(displ,indsym,natom,nsym,qptn,rprimd,symrel,symafm) 5786 end if 5787 5788 ii=1 5789 ABI_MALLOC(matrx,(2,(3*natom*(3*natom+1))/2)) 5790 do i2=1,3*natom 5791 do i1=1,i2 5792 matrx(1,ii)=displ(1+2*(i1-1)+2*(i2-1)*3*natom) 5793 matrx(2,ii)=displ(2+2*(i1-1)+2*(i2-1)*3*natom) 5794 ii=ii+1 5795 end do 5796 end do 5797 5798 ABI_MALLOC(zhpev1,(2,2*3*natom-1)) 5799 ABI_MALLOC(zhpev2,(3*3*natom-2)) 5800 5801 call ZHPEV ('V','U',3*natom,matrx,eigval,eigvec,3*natom,zhpev1,zhpev2,ier) 5802 ABI_CHECK(ier == 0, sjoin('zhpev returned:', itoa(ier))) 5803 5804 ABI_FREE(matrx) 5805 ABI_FREE(zhpev1) 5806 ABI_FREE(zhpev2) 5807 5808 if (debug) then 5809 ! Check the orthonormality of the eigenvectors 5810 do imode=1,3*natom 5811 do jmode=imode,3*natom 5812 indexi=2*3*natom*(imode-1) 5813 indexj=2*3*natom*(jmode-1) 5814 sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom)) 5815 write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod 5816 end do 5817 end do 5818 end if 5819 5820 !*********************************************************************** 5821 5822 ! Get the phonon frequencies (negative by convention, if the eigenvalue of the dynamical matrix is negative) 5823 do imode=1,3*natom 5824 if(eigval(imode)>=1.0d-16)then 5825 phfrq(imode)=sqrt(eigval(imode)) 5826 else if(eigval(imode)>=-1.0d-16)then 5827 phfrq(imode)=zero 5828 else 5829 phfrq(imode)=-sqrt(-eigval(imode)) 5830 end if 5831 end do 5832 5833 ! Fix the phase of the eigenvectors 5834 call fxphas_seq(eigvec,dum, 0, 0, 1, 3*natom*3*natom, 0, 3*natom, 3*natom, 0) 5835 5836 ! Normalise the eigenvectors 5837 call pheigvec_normalize(natom, eigvec) 5838 5839 ! Get the phonon displacements 5840 call phdispl_from_eigvec(natom, ntypat, typat, amu, eigvec, displ) 5841 5842 if (debug) then 5843 write(std_out,'(a)')' Phonon eigenvectors and displacements ' 5844 do imode=1,3*natom 5845 indexi=2*3*natom*(imode-1) 5846 write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' eigvec(1:6*natom)=',eigvec(indexi+1:indexi+6*natom) 5847 write(std_out,'(a,i4,a,12es16.6)')' imode=',imode,' displ(1:6*natom)=',displ(indexi+1:indexi+6*natom) 5848 end do 5849 5850 ! Check the orthonormality of the eigenvectors 5851 do imode=1,3*natom 5852 do jmode=imode,3*natom 5853 indexi=2*3*natom*(imode-1) 5854 indexj=2*3*natom*(jmode-1) 5855 sc_prod=sum(eigvec(indexi+1:indexi+6*natom)*eigvec(indexj+1:indexj+6*natom)) 5856 write(std_out,'(a,2i4,a,es16.6)')' imode,jmode=',imode,jmode,' real scalar product =',sc_prod 5857 end do 5858 end do 5859 end if 5860 5861 !call timab(1751, 2, tsec) 5862 5863 end subroutine dfpt_phfrq
m_dynmat/dfpt_prtph [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dfpt_prtph
FUNCTION
Print the phonon frequencies, on unit 6 as well as the printing unit (except if the associated number -iout- is negative), and for the latter, in Hartree, meV, Thz, Kelvin or cm-1. If eivec==1,2, also print the eigenmodes : displacements in cartesian coordinates. If eivec==4, generate output files for band2eps (drawing tool for the phonon band structure
INPUTS
displ(2,3*natom,3*natom)= 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. eivec=(if eivec==0, the eigendisplacements are not printed, if eivec==1,2, the eigendisplacements are printed, if eivec==4, files for band2eps enunit=units for output of the phonon frequencies : 0=> Hartree and cm-1, 1=> eV and Thz, other=> Ha,Thz,eV,cm-1 and K iout= unit for long print (if negative, the routine only print on unit 6, and in Hartree only). natom= number of atom phfreq(3*natom)= phonon frequencies in Hartree qphnrm=phonon wavevector normalisation factor qphon(3)=phonon wavevector
OUTPUT
Only printing
NOTES
called by one processor only
SOURCE
6014 subroutine dfpt_prtph(displ,eivec,enunit,iout,natom,phfrq,qphnrm,qphon) 6015 6016 !Arguments ------------------------------- 6017 !scalars 6018 integer,intent(in) :: eivec,enunit,iout,natom 6019 real(dp),intent(in) :: qphnrm 6020 !arrays 6021 real(dp),intent(in) :: displ(2,3*natom,3*natom),phfrq(3*natom),qphon(3) 6022 6023 !Local variables ------------------------- 6024 !scalars 6025 integer :: i,idir,ii,imode,jj 6026 real(dp) :: tolerance 6027 logical :: t_degenerate 6028 character(len=500) :: msg 6029 !arrays 6030 real(dp) :: vecti(3),vectr(3) 6031 character(len=1) :: metacharacter(3*natom) 6032 6033 ! ********************************************************************* 6034 6035 !Check the value of eivec 6036 if (all(eivec /= [0,1,2,4])) then 6037 write(msg, '(a,i0,a,a)' )& 6038 'In the calling subroutine, eivec is',eivec,ch10,& 6039 'but allowed values are between 0 and 4.' 6040 ABI_BUG(msg) 6041 end if 6042 6043 !write the phonon frequencies on unit std_out 6044 write(msg,'(4a)' )' ',ch10,' phonon wavelength (reduced coordinates) , ','norm, and energies in hartree' 6045 call wrtout(std_out,msg) 6046 6047 !The next format should be rewritten 6048 write(msg,'(a,4f5.2)' )' ',(qphon(i),i=1,3),qphnrm 6049 call wrtout(std_out,msg) 6050 do jj=1,3*natom,5 6051 if (3*natom-jj<5) then 6052 write(msg,'(5es17.9)') (phfrq(ii),ii=jj,3*natom) 6053 else 6054 write(msg,'(5es17.9)') (phfrq(ii),ii=jj,jj+4) 6055 end if 6056 call wrtout(std_out,msg) 6057 end do 6058 write(msg,'(a,a,es17.9)') ch10,' Zero Point Motion energy (sum of freqs/2)=',sum(phfrq(1:3*natom))/2 6059 call wrtout(std_out,msg) 6060 6061 !Put the wavevector in nice format 6062 if(iout>=0)then 6063 call wrtout(iout,' ') 6064 if(qphnrm/=0.0_dp)then 6065 write(msg, '(a,3f9.5)' )& 6066 ' Phonon wavevector (reduced coordinates) :',(qphon(i)/qphnrm+tol10,i=1,3) 6067 else 6068 write(msg, '(3a,3f9.5)' )& 6069 ' Phonon at Gamma, with non-analyticity in the',ch10,& 6070 ' direction (cartesian coordinates)',qphon(1:3)+tol10 6071 end if 6072 call wrtout(iout,msg) 6073 6074 ! Write it, in different units. 6075 if(enunit/=1)then 6076 write(iout, '(a)' )' Phonon energies in Hartree :' 6077 do jj=1,3*natom,5 6078 if (3*natom-jj<5) then 6079 write(msg, '(1x,5es14.6)') (phfrq(ii),ii=jj,3*natom) 6080 else 6081 write(msg, '(1x,5es14.6)') (phfrq(ii),ii=jj,jj+4) 6082 end if 6083 call wrtout(iout,msg) 6084 end do 6085 end if 6086 if(enunit/=0)then 6087 write(iout, '(a)' )' Phonon energies in meV :' 6088 do jj=1,3*natom,5 6089 if (3*natom-jj<5) then 6090 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,3*natom) 6091 else 6092 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_eV*1.0d3,ii=jj,jj+4) 6093 end if 6094 call wrtout(iout,msg) 6095 end do 6096 end if 6097 if(enunit/=1)then 6098 write(iout, '(a)' )' Phonon frequencies in cm-1 :' 6099 do jj=1,3*natom,5 6100 if (3*natom-jj<5) then 6101 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,3*natom) 6102 else 6103 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_cmm1,ii=jj,jj+4) 6104 end if 6105 call wrtout(iout,msg) 6106 end do 6107 end if 6108 if(enunit/=0)then 6109 write(iout, '(a)' )' Phonon frequencies in Thz :' 6110 do jj=1,3*natom,5 6111 if (3*natom-jj<5) then 6112 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,3*natom) 6113 else 6114 write(msg, '("-",5es14.6)') (phfrq(ii)*Ha_THz,ii=jj,jj+4) 6115 end if 6116 call wrtout(iout,msg) 6117 end do 6118 end if 6119 if(enunit/=0.and.enunit/=1)then 6120 write(iout, '(a)' )' Phonon energies in Kelvin :' 6121 do jj=1,3*natom,5 6122 if (3*natom-jj<5) then 6123 write(msg, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,3*natom) 6124 else 6125 write(msg, '("-",5es14.6)') (phfrq(ii)/kb_HaK,ii=jj,jj+4) 6126 end if 6127 call wrtout(iout,msg) 6128 end do 6129 end if 6130 end if 6131 6132 !Take care of the eigendisplacements 6133 if(eivec==1 .or. eivec==2)then 6134 write(msg, '(a,a,a,a,a,a,a,a)' ) ch10,& 6135 ' Eigendisplacements ',ch10,& 6136 ' (will be given, for each mode : in cartesian coordinates',ch10,& 6137 ' for each atom the real part of the displacement vector,',ch10,& 6138 ' then the imaginary part of the displacement vector - absolute values smaller than 1.0d-7 are set to zero)' 6139 call wrtout(std_out,msg) 6140 if(iout>=0) then 6141 call wrtout(iout,msg) 6142 end if 6143 6144 ! Examine the degeneracy of each mode. The portability of the echo of the eigendisplacements 6145 ! is very hard to obtain, and has not been attempted. 6146 do imode=1,3*natom 6147 ! The degenerate modes are not portable 6148 t_degenerate=.false. 6149 if(imode>1)then 6150 if(phfrq(imode)-phfrq(imode-1)<tol6)t_degenerate=.true. 6151 end if 6152 if(imode<3*natom)then 6153 if(phfrq(imode+1)-phfrq(imode)<tol6)t_degenerate=.true. 6154 end if 6155 metacharacter(imode)=';'; if(t_degenerate)metacharacter(imode)='-' 6156 end do 6157 6158 do imode=1,3*natom 6159 write(msg,'(a,i4,a,es16.6)' )' Mode number ',imode,' Energy',phfrq(imode) 6160 call wrtout(std_out,msg) 6161 if(iout>=0)then 6162 write(msg, '(a,i4,a,es16.6)' )' Mode number ',imode,' Energy',phfrq(imode) 6163 call wrtout(iout,msg) 6164 end if 6165 tolerance=1.0d-7 6166 if(abs(phfrq(imode))<1.0d-5)tolerance=2.0d-7 6167 if(phfrq(imode)<1.0d-5)then 6168 write(msg,'(3a)' )' Attention : low frequency mode.',ch10,& 6169 ' (Could be unstable or acoustic mode)' 6170 call wrtout(std_out,msg) 6171 if(iout>=0)then 6172 write(iout, '(3a)' )' Attention : low frequency mode.',ch10,& 6173 ' (Could be unstable or acoustic mode)' 6174 end if 6175 end if 6176 do ii=1,natom 6177 do idir=1,3 6178 vectr(idir)=displ(1,idir+(ii-1)*3,imode) 6179 if(abs(vectr(idir))<tolerance)vectr(idir)=0.0_dp 6180 vecti(idir)=displ(2,idir+(ii-1)*3,imode) 6181 if(abs(vecti(idir))<tolerance)vecti(idir)=0.0_dp 6182 end do 6183 write(msg,'(i4,3es16.8,a,4x,3es16.8)' ) ii,vectr(:),ch10,vecti(:) 6184 call wrtout(std_out,msg) 6185 if(iout>=0)then 6186 write(msg,'(a,i3,3es16.8,2a,3x,3es16.8)') metacharacter(imode),ii,vectr(:),ch10,& 6187 metacharacter(imode), vecti(:) 6188 call wrtout(iout,msg) 6189 end if 6190 end do 6191 end do 6192 end if 6193 6194 end subroutine dfpt_prtph
m_dynmat/dfpt_sydy [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dfpt_sydy
FUNCTION
Symmetrize dynamical matrix (eventually diagonal wrt to the atoms) Unsymmetrized dynamical matrix is input as dyfrow; symmetrized dynamical matrix is then placed in sdyfro. If nsym=1 simply copy dyfrow into sdyfro.
INPUTS
cplex=1 if dynamical matrix is real, 2 if it is complex dyfrow(3,3,natom,1+(natom-1)*nondiag)=unsymmetrized dynamical matrix 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. natom=number of atoms in cell. nondiag=0 if dynamical matrix is diagonal with respect to atoms nsym=number of symmetry operators in group. qphon(3)=wavevector of the phonon symq(4,2,nsym)=1 if symmetry preserves present qpoint. From littlegroup_q symrec(3,3,nsym)=symmetries of group in terms of operations on real space primitive translations (integers).
OUTPUT
sdyfro(3,3,natom,1+(natom-1)*nondiag)=symmetrized dynamical matrix
NOTES
Symmetrization of gradients with respect to reduced coordinates tn is conducted according to the expression $[d(e)/d(t(n,a))]_{symmetrized} = (1/Nsym)*Sum(S)*symrec(n,m,S)* [d(e)/d(t(m,b))]_{unsymmetrized}$ where $t(m,b)= (symrel^{-1})(m,n)*(t(n,a)-tnons(n))$ and tnons is a possible nonsymmorphic translation. The label "b" here refers to the atom which gets rotated into "a" under symmetry "S". symrel is the symmetry matrix in real space, which is the inverse transpose of symrec. symrec is the symmetry matrix in reciprocal space. $sym_{cartesian} = R * symrel * R^{-1} = G * symrec * G^{-1}$ where the columns of R and G are the dimensional primitive translations in real and reciprocal space respectively. Note the use of "symrec" in the symmetrization expression above.
SOURCE
2359 subroutine dfpt_sydy(cplex,dyfrow,indsym,natom,nondiag,nsym,qphon,sdyfro,symq,symrec) 2360 2361 !Arguments ------------------------------- 2362 !scalars 2363 integer,intent(in) :: cplex,natom,nondiag,nsym 2364 !arrays 2365 integer,intent(in) :: indsym(4,nsym,natom),symq(4,2,nsym),symrec(3,3,nsym) 2366 real(dp),intent(in) :: dyfrow(cplex,3,3,natom,1+(natom-1)*nondiag),qphon(3) 2367 real(dp),intent(out) :: sdyfro(cplex,3,3,natom,1+(natom-1)*nondiag) 2368 2369 !Local variables ------------------------- 2370 !scalars 2371 integer :: ia,indi,indj,isym,ja,kappa,mu,natom_nondiag,nsym_used,nu 2372 logical :: qeq0 2373 real(dp) :: arg,div,phasei,phaser 2374 !arrays 2375 real(dp) :: work(cplex,3,3) 2376 2377 ! ********************************************************************* 2378 2379 if (nsym==1) then 2380 2381 ! Only symmetry is identity so simply copy 2382 sdyfro(:,:,:,:,:)=dyfrow(:,:,:,:,:) 2383 2384 else 2385 2386 ! Actually carry out symmetrization 2387 sdyfro(:,:,:,:,:)=zero 2388 qeq0=(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14) 2389 ! === Diagonal dyn. matrix OR q=0 2390 if (nondiag==0.or.qeq0) then 2391 natom_nondiag=1;if (nondiag==1) natom_nondiag=natom 2392 do ja=1,natom_nondiag 2393 do ia=1,natom 2394 do isym=1,nsym 2395 indi=indsym(4,isym,ia) 2396 indj=1;if (nondiag==1) indj=indsym(4,isym,ja) 2397 work(:,:,:)=zero 2398 do mu=1,3 2399 do nu=1,3 2400 do kappa=1,3 2401 work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj) 2402 end do 2403 end do 2404 end do 2405 do mu=1,3 2406 do nu=1,3 2407 do kappa=1,3 2408 sdyfro(:,kappa,mu,ia,ja)=sdyfro(:,kappa,mu,ia,ja)+symrec(mu,nu,isym)*work(:,kappa,nu) 2409 end do 2410 end do 2411 end do 2412 end do 2413 end do 2414 end do 2415 div=one/dble(nsym) 2416 sdyfro(:,:,:,:,:)=div*sdyfro(:,:,:,:,:) 2417 ! === Non diagonal dyn. matrix AND q<>0 2418 else 2419 do ja=1,natom 2420 do ia=1,natom 2421 nsym_used=0 2422 do isym=1,nsym 2423 if (symq(4,1,isym)==1) then 2424 arg=two_pi*(qphon(1)*(indsym(1,isym,ia)-indsym(1,isym,ja)) & 2425 & +qphon(2)*(indsym(2,isym,ia)-indsym(2,isym,ja)) & 2426 & +qphon(3)*(indsym(3,isym,ia)-indsym(3,isym,ja))) 2427 phaser=cos(arg);phasei=sin(arg) 2428 nsym_used=nsym_used+1 2429 indi=indsym(4,isym,ia) 2430 indj=indsym(4,isym,ja) 2431 work(:,:,:)=zero 2432 do mu=1,3 2433 do nu=1,3 2434 do kappa=1,3 2435 work(:,mu,kappa)=work(:,mu,kappa)+symrec(mu,nu,isym)*dyfrow(:,nu,kappa,indi,indj) 2436 end do 2437 end do 2438 end do 2439 do mu=1,3 2440 do nu=1,3 2441 do kappa=1,3 2442 sdyfro(1,kappa,mu,ia,ja)=sdyfro(1,kappa,mu,ia,ja) & 2443 & +symrec(mu,nu,isym)*(work(1,kappa,nu)*phaser-work(2,kappa,nu)*phasei) 2444 end do 2445 end do 2446 end do 2447 if (cplex==2) then 2448 do mu=1,3 2449 do nu=1,3 2450 do kappa=1,3 2451 sdyfro(2,kappa,mu,ia,ja)=sdyfro(2,kappa,mu,ia,ja) & 2452 & +symrec(mu,nu,isym)*(work(1,kappa,nu)*phasei+work(2,kappa,nu)*phaser) 2453 end do 2454 end do 2455 end do 2456 end if 2457 end if 2458 end do 2459 div=one/dble(nsym_used) 2460 sdyfro(:,:,:,ia,ja)=div*sdyfro(:,:,:,ia,ja) 2461 end do 2462 end do 2463 end if 2464 2465 end if 2466 2467 end subroutine dfpt_sydy
m_dynmat/dfpt_sygra [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
dfpt_sygra
FUNCTION
Symmetrize derivatives of energy with respect to coordinates, as appearing in phonon calculations. Unsymmetrized gradients are input as deunsy; symmetrized grads are then placed in desym. If nsym=1 simply copy deunsy into desym (only symmetry is identity). The index of the initial perturbation is needed, in case there is a change of atom position (moved in another cell) due to the symmetry operation.
INPUTS
natom=number of atoms in cell deunsy(2,3,natom)=unsymmetrized gradients wrt dimensionless tn (hartree) note: there is a real and a imaginary part ... indsym(4,nsym,natom)=label given by subroutine symatm, indicating atom label which gets rotated into given atom by given symmetry (first three elements are related primitive translation-- see symatm where this is computed) nsym=number of symmetry operators in group ipert=index of the initial perturbation qpt(3)= wavevector of the phonon, in reduced coordinates symrec(3,3,nsym)=symmetries of group in terms of operations on reciprocal space primitive translations--see comments below
OUTPUT
desym(2,3,natom)=symmetrized gradients wrt dimensionless tn (hartree)
NOTES
Written by X. Gonze starting from sygrad, written by D.C. Allan: introduction of the q vector for phonon symmetrization This routine should once be merged with sygrad...
SOURCE
2232 subroutine dfpt_sygra(natom,desym,deunsy,indsym,ipert,nsym,qpt,symrec) 2233 2234 !Arguments ------------------------------- 2235 !scalars 2236 integer,intent(in) :: ipert,natom,nsym 2237 !arrays 2238 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym) 2239 real(dp),intent(in) :: deunsy(2,3,natom),qpt(3) 2240 real(dp),intent(out) :: desym(2,3,natom) 2241 2242 !Local variables ------------------------- 2243 !scalars 2244 integer :: ia,ind,isym,mu 2245 real(dp) :: arg,im,re,sumi,sumr 2246 2247 ! ********************************************************************* 2248 2249 !DEBUG 2250 !write(std_out,*)' dfpt_sygra : enter ' 2251 !write(std_out,*)' dfpt_sygra : qpt(:)',qpt(:) 2252 !do ia=1,natom 2253 !do mu=1,3 2254 !write(std_out,*)' dfpt_sygra : deunsy(:2,mu,ia)',deunsy(:2,mu,ia) 2255 !enddo 2256 !enddo 2257 !ENDDEBUG 2258 2259 if (nsym==1) then 2260 2261 ! Only symmetry is identity so simply copy 2262 desym(:,:,:)=deunsy(:,:,:) 2263 2264 else 2265 2266 ! Actually conduct symmetrization 2267 ! write(std_out,*)' dfpt_sygra : desym(:2,:3,:natom),qpt(:)',desym(:2,:3,:natom),qpt(:) 2268 do ia=1,natom 2269 ! write(std_out,*)' dfpt_sygra : ia=',ia 2270 do mu=1,3 2271 sumr=zero 2272 sumi=zero 2273 ! write(std_out,*)' dfpt_sygra : mu=',mu 2274 do isym=1,nsym 2275 ind=indsym(4,isym,ia) 2276 ! Must shift the atoms back to the unit cell. 2277 ! arg=two_pi*( qpt(1)*indsym(1,isym,ia)& 2278 ! & +qpt(2)*indsym(2,isym,ia)& 2279 ! & +qpt(3)*indsym(3,isym,ia) ) 2280 ! Selection of non-zero q point, to avoid ipert being outside the 1 ... natom range 2281 if(qpt(1)**2+qpt(2)**2+qpt(3)**2 > tol16)then 2282 arg=two_pi*( qpt(1)*(indsym(1,isym,ia)-indsym(1,isym,ipert))& 2283 & +qpt(2)* (indsym(2,isym,ia)-indsym(2,isym,ipert))& 2284 & +qpt(3)* (indsym(3,isym,ia)-indsym(3,isym,ipert))) 2285 else 2286 arg=zero 2287 end if 2288 2289 re=dble(symrec(mu,1,isym))*deunsy(1,1,ind)+& 2290 & dble(symrec(mu,2,isym))*deunsy(1,2,ind)+& 2291 & dble(symrec(mu,3,isym))*deunsy(1,3,ind) 2292 im=dble(symrec(mu,1,isym))*deunsy(2,1,ind)+& 2293 & dble(symrec(mu,2,isym))*deunsy(2,2,ind)+& 2294 & dble(symrec(mu,3,isym))*deunsy(2,3,ind) 2295 sumr=sumr+re*cos(arg)-im*sin(arg) 2296 sumi=sumi+re*sin(arg)+im*cos(arg) 2297 ! sumr=sumr+re 2298 ! sumi=sumi+im 2299 ! write(std_out,*)' dfpt_sygra : isym,indsym(4,isym,ia),arg,re,im,sumr,sumi',& 2300 ! & isym,indsym(4,isym,ia),arg,re,im,sumr,sumi 2301 end do 2302 desym(1,mu,ia)=sumr/dble(nsym) 2303 desym(2,mu,ia)=sumi/dble(nsym) 2304 ! write(std_out,*)' dfpt_sygra : desym(:,mu,ia)',desym(:,mu,ia) 2305 end do 2306 end do 2307 end if 2308 2309 end subroutine dfpt_sygra
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
SOURCE
3428 subroutine dist9(acell,dist,gprim,natom,nrpt,rcan,rprim,rpt) 3429 3430 !Arguments ------------------------------- 3431 !scalars 3432 integer,intent(in) :: natom,nrpt 3433 !arrays 3434 real(dp),intent(in) :: acell(3),gprim(3,3),rcan(3,natom),rprim(3,3) 3435 real(dp),intent(in) :: rpt(3,nrpt) 3436 real(dp),intent(out) :: dist(natom,natom,nrpt) 3437 3438 !Local variables ------------------------- 3439 !scalars 3440 integer :: ia,ib,ii,irpt 3441 !arrays 3442 real(dp) :: ra(3),rb(3),rdiff(3),red(3),rptcar(3),xred(3) 3443 3444 ! ********************************************************************* 3445 3446 !BIG loop on all generic atoms 3447 do ia=1,natom 3448 ! First transform canonical coordinates to reduced coordinates 3449 do ii=1,3 3450 xred(ii)=gprim(1,ii)*rcan(1,ia)+gprim(2,ii)*rcan(2,ia)+gprim(3,ii)*rcan(3,ia) 3451 end do 3452 ! Then to cartesian coordinates 3453 ra(:)=xred(1)*acell(1)*rprim(:,1)+xred(2)*acell(2)*rprim(:,2)+xred(3)*acell(3)*rprim(:,3) 3454 do ib=1,natom 3455 do ii=1,3 3456 xred(ii)=gprim(1,ii)*rcan(1,ib)+gprim(2,ii)*rcan(2,ib)+gprim(3,ii)*rcan(3,ib) 3457 end do 3458 do ii=1,3 3459 rb(ii)=xred(1)*acell(1)*rprim(ii,1)+xred(2)*acell(2)*rprim(ii,2)+xred(3)*acell(3)*rprim(ii,3) 3460 end do 3461 do irpt=1,nrpt 3462 ! First transform it to reduced coordinates 3463 do ii=1,3 3464 red(ii)=gprim(1,ii)*rpt(1,irpt)+gprim(2,ii)*rpt(2,irpt)+gprim(3,ii)*rpt(3,irpt) 3465 end do 3466 ! Then to cartesian coordinates 3467 do ii=1,3 3468 rptcar(ii)=red(1)*acell(1)*rprim(ii,1)+red(2)*acell(2)*rprim(ii,2)+red(3)*acell(3)*rprim(ii,3) 3469 end do 3470 do ii=1,3 3471 rdiff(ii)=-rptcar(ii)+ra(ii)-rb(ii) 3472 end do 3473 dist(ia,ib,irpt)=(rdiff(1)**2+rdiff(2)**2+rdiff(3)**2)**0.5 3474 end do 3475 end do 3476 end do 3477 3478 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
SOURCE
5337 subroutine dymfz9(dynmat,natom,nqpt,gprim,option,spqpt,trans) 5338 5339 !Arguments ------------------------------- 5340 !scalars 5341 integer,intent(in) :: natom,nqpt,option 5342 !arrays 5343 real(dp),intent(in) :: gprim(3,3),spqpt(3,nqpt),trans(3,natom) 5344 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt) 5345 5346 !Local variables ------------------------- 5347 !scalars 5348 integer :: ia,ib,iqpt,mu,nu 5349 real(dp) :: im,ktrans,re 5350 !arrays 5351 real(dp) :: kk(3) 5352 5353 ! ********************************************************************* 5354 5355 do iqpt=1,nqpt 5356 ! Definition of q in normalized reciprocal space 5357 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 5358 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 5359 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 5360 5361 if(option==1)then 5362 kk(1)=-kk(1) 5363 kk(2)=-kk(2) 5364 kk(3)=-kk(3) 5365 end if 5366 5367 do ia=1,natom 5368 do ib=1,natom 5369 ! Product of q with the differences between the two atomic translations 5370 ktrans=kk(1)*(trans(1,ia)-trans(1,ib))+kk(2)*(trans(2,ia)-trans(2,ib))+kk(3)*(trans(3,ia)-trans(3,ib)) 5371 do mu=1,3 5372 do nu=1,3 5373 re=dynmat(1,mu,ia,nu,ib,iqpt) 5374 im=dynmat(2,mu,ia,nu,ib,iqpt) 5375 ! Transformation of the Old dynamical matrices by New ones by multiplication by a phase shift 5376 dynmat(1,mu,ia,nu,ib,iqpt)=re*cos(two_pi*ktrans)-im*sin(two_pi*ktrans) 5377 dynmat(2,mu,ia,nu,ib,iqpt)=re*sin(two_pi*ktrans)+im*cos(two_pi*ktrans) 5378 end do 5379 end do 5380 end do 5381 end do 5382 end do 5383 5384 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)!!)
SOURCE
3716 subroutine dynmat_dq(qpt,natom,gprim,nrpt,rpt,atmfrc,wghatm,dddq) 3717 3718 !Arguments ------------------------------- 3719 !scalars 3720 integer,intent(in) :: natom,nrpt 3721 !arrays 3722 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt(3) 3723 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 3724 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 3725 real(dp),intent(out) :: dddq(2,3,natom,3,natom,3) 3726 3727 !Local variables ------------------------- 3728 !scalars 3729 integer :: ia,ib,irpt,mu,nu,ii 3730 real(dp) :: im,kr,re 3731 !arrays 3732 real(dp) :: kk(3),fact(2,3) 3733 3734 ! ********************************************************************* 3735 3736 dddq = zero 3737 3738 do irpt=1,nrpt 3739 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3740 kk(1) = qpt(1)*gprim(1,1)+ qpt(2)*gprim(1,2) + qpt(3)*gprim(1,3) 3741 kk(2) = qpt(1)*gprim(2,1)+ qpt(2)*gprim(2,2) + qpt(3)*gprim(2,3) 3742 kk(3) = qpt(1)*gprim(3,1)+ qpt(2)*gprim(3,2) + qpt(3)*gprim(3,3) 3743 3744 ! Product of k and r 3745 kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3746 3747 ! Get phase factor 3748 re=cos(two_pi*kr); im=sin(two_pi*kr) 3749 3750 ! Inner loop on atoms and directions 3751 do ib=1,natom 3752 do ia=1,natom 3753 if (abs(wghatm(ia,ib,irpt))>1.0d-10) then 3754 ! take into account rotation due to i. 3755 fact(1,:) = -im * wghatm(ia,ib,irpt) * rpt(:,irpt) 3756 fact(2,:) = re * wghatm(ia,ib,irpt) * rpt(:,irpt) 3757 do nu=1,3 3758 do mu=1,3 3759 ! Real and imaginary part of the dynamical matrices 3760 ! Atmfrc should be real 3761 do ii=1,3 3762 dddq(1,mu,ia,nu,ib,ii) = dddq(1,mu,ia,nu,ib,ii) + fact(1,ii) * atmfrc(mu,ia,nu,ib,irpt) 3763 dddq(2,mu,ia,nu,ib,ii) = dddq(2,mu,ia,nu,ib,ii) + fact(2,ii) * atmfrc(mu,ia,nu,ib,irpt) 3764 end do 3765 end do 3766 end do 3767 end if 3768 end do 3769 end do 3770 end do 3771 3772 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
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
6350 subroutine ftgam (wghatm,gam_qpt,gam_rpt,natom,nqpt,nrpt,qtor,coskr, sinkr) 6351 6352 !Arguments ------------------------------- 6353 !scalars 6354 integer,intent(in) :: natom,nqpt,nrpt,qtor 6355 !arrays 6356 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 6357 real(dp),intent(inout) :: gam_qpt(2,3*natom*3*natom,nqpt) 6358 real(dp),intent(inout) :: gam_rpt(2,3*natom*3*natom,nrpt) 6359 real(dp),intent(in) :: coskr(nqpt,nrpt) 6360 real(dp),intent(in) :: sinkr(nqpt,nrpt) 6361 6362 !Local variables ------------------------- 6363 !scalars 6364 integer :: iatom,idir,ip,iqpt,irpt,jatom,jdir 6365 real(dp) :: im,re 6366 character(len=500) :: msg 6367 6368 ! ********************************************************************* 6369 6370 select case (qtor) 6371 case (1) 6372 ! Recip to real space 6373 gam_rpt(:,:,:) = zero 6374 do irpt=1,nrpt 6375 do iqpt=1,nqpt 6376 ! Get the phase factor with normalization! 6377 re=coskr(iqpt,irpt) 6378 im=sinkr(iqpt,irpt) 6379 do ip=1,3*natom*3*natom 6380 ! Real and imaginary part of the real-space gam matrices 6381 gam_rpt(1,ip,irpt) = gam_rpt(1,ip,irpt) + re*gam_qpt(1,ip,iqpt) + im*gam_qpt(2,ip,iqpt) 6382 gam_rpt(2,ip,irpt) = gam_rpt(2,ip,irpt) + re*gam_qpt(2,ip,iqpt) - im*gam_qpt(1,ip,iqpt) 6383 end do 6384 end do 6385 end do 6386 gam_rpt = gam_rpt/nqpt 6387 6388 case (0) 6389 ! Recip space from real space 6390 gam_qpt(:,:,:)=zero 6391 6392 do irpt=1,nrpt 6393 do iqpt=1,nqpt 6394 6395 do iatom=1,natom 6396 do jatom=1,natom 6397 re = coskr(iqpt,irpt)*wghatm(iatom,jatom,irpt) 6398 im = sinkr(iqpt,irpt)*wghatm(iatom,jatom,irpt) 6399 6400 do idir=1,3 6401 do jdir=1,3 6402 ! Get phase factor 6403 6404 ip= jdir + (jatom-1)*3 + (idir-1)*3*natom + (iatom-1)*9*natom 6405 ! Real and imaginary part of the interatomic forces 6406 gam_qpt(1,ip,iqpt) = gam_qpt(1,ip,iqpt) + re*gam_rpt(1,ip,irpt) - im*gam_rpt(2,ip,irpt) 6407 !DEBUG 6408 gam_qpt(2,ip,iqpt) = gam_qpt(2,ip,iqpt) + im*gam_rpt(1,ip,irpt) + re*gam_rpt(2,ip,irpt) 6409 !ENDDEBUG 6410 end do ! end jdir 6411 end do ! end idir 6412 end do 6413 end do ! end iatom 6414 6415 end do ! end iqpt 6416 end do ! end irpt 6417 6418 case default 6419 write(msg,'(a,i0,a)' )'The only allowed values for qtor are 0 or 1, while qtor= ',qtor,' has been required.' 6420 ABI_BUG(msg) 6421 end select 6422 6423 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
SOURCE
6448 subroutine ftgam_init (gprim,nqpt,nrpt,qpt_full,rpt,coskr, sinkr) 6449 6450 !Arguments ------------------------------- 6451 !scalars 6452 integer,intent(in) :: nqpt,nrpt 6453 !arrays 6454 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),qpt_full(3,nqpt) 6455 real(dp),intent(out) :: coskr(nqpt,nrpt) 6456 real(dp),intent(out) :: sinkr(nqpt,nrpt) 6457 6458 !Local variables ------------------------- 6459 !scalars 6460 integer :: iqpt,irpt 6461 real(dp) :: kr 6462 !arrays 6463 real(dp) :: kk(3) 6464 6465 ! ********************************************************************* 6466 6467 ! Prepare the phase factors 6468 do iqpt=1,nqpt 6469 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 6470 kk(1) = qpt_full(1,iqpt)*gprim(1,1) + qpt_full(2,iqpt)*gprim(1,2) + qpt_full(3,iqpt)*gprim(1,3) 6471 kk(2) = qpt_full(1,iqpt)*gprim(2,1) + qpt_full(2,iqpt)*gprim(2,2) + qpt_full(3,iqpt)*gprim(2,3) 6472 kk(3) = qpt_full(1,iqpt)*gprim(3,1) + qpt_full(2,iqpt)*gprim(3,2) + qpt_full(3,iqpt)*gprim(3,3) 6473 do irpt=1,nrpt 6474 ! Product of k and r 6475 kr = kk(1)*rpt(1,irpt)+ kk(2)*rpt(2,irpt)+ kk(3)*rpt(3,irpt) 6476 coskr(iqpt,irpt)=cos(two_pi*kr) 6477 sinkr(iqpt,irpt)=sin(two_pi*kr) 6478 end do 6479 end do 6480 6481 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 comm=MPI communicator.
OUTPUT
atmfrc(3,natom,3,natom,nrpt)= Interatomic Forces in real space.
SOURCE
3508 subroutine ftifc_q2r(atmfrc,dynmat,gprim,natom,nqpt,nrpt,rpt,spqpt,comm) 3509 3510 !Arguments ------------------------------- 3511 !scalars 3512 integer,intent(in) :: natom,nqpt,nrpt,comm 3513 !arrays 3514 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt) 3515 real(dp),intent(out) :: atmfrc(3,natom,3,natom,nrpt) 3516 real(dp),intent(in) :: dynmat(2,3,natom,3,natom,nqpt) 3517 3518 !Local variables ------------------------- 3519 !scalars 3520 integer :: ia,ib,iqpt,irpt,mu,nu,nprocs,my_rank,ierr 3521 real(dp) :: im,kr,re 3522 !arrays 3523 real(dp) :: kk(3) 3524 3525 ! ********************************************************************* 3526 3527 nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 3528 3529 ! Interatomic Forces from Dynamical Matrices 3530 atmfrc = zero 3531 do irpt=1,nrpt 3532 if (mod(irpt, nprocs) /= my_rank) cycle ! mpi-parallelism 3533 do iqpt=1,nqpt 3534 3535 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3536 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)*gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 3537 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)*gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 3538 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)*gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 3539 3540 ! Product of k and r 3541 kr=kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3542 3543 ! Get the phase factor 3544 re=cos(two_pi*kr) 3545 im=sin(two_pi*kr) 3546 3547 ! Now, big inner loops on atoms and directions 3548 ! The indices are ordered to give better speed 3549 do ib=1,natom 3550 do nu=1,3 3551 do ia=1,natom 3552 do mu=1,3 3553 ! Real and imaginary part of the interatomic forces 3554 atmfrc(mu,ia,nu,ib,irpt)=atmfrc(mu,ia,nu,ib,irpt) & 3555 +re*dynmat(1,mu,ia,nu,ib,iqpt)& 3556 +im*dynmat(2,mu,ia,nu,ib,iqpt) 3557 !The imaginary part should be equal to zero !!!!!! 3558 !atmfrc(2,mu,ia,nu,ib,irpt)=atmfrc(2,mu,ia,nu,ib,irpt) & 3559 ! +re*dynmat(2,mu,ia,nu,ib,iqpt) & 3560 ! -im*dynmat(1,mu,ia,nu,ib,iqpt) 3561 end do 3562 end do 3563 end do 3564 end do 3565 3566 end do 3567 end do 3568 3569 call xmpi_sum(atmfrc, comm, ierr) 3570 !The sumifc has to be weighted by a normalization factor of 1/nqpt 3571 atmfrc = atmfrc/nqpt 3572 3573 end subroutine ftifc_q2r
m_dynmat/ftifc_r2q [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
ftifc_r2q
FUNCTION
Generates the Fourier transform of the interatomic forces to obtain dynamical matrices in reciprocal space: R --> q.
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 comm: MPI communicator
OUTPUT
dynmat(2,3,natom,3,natom,nqpt)= Dynamical matrices coming from the Derivative Data Base
SOURCE
3604 subroutine ftifc_r2q(atmfrc, dynmat, gprim, natom, nqpt, nrpt, rpt, spqpt, wghatm, comm) 3605 3606 !Arguments ------------------------------- 3607 !scalars 3608 integer,intent(in) :: natom,nqpt,nrpt,comm 3609 !arrays 3610 real(dp),intent(in) :: gprim(3,3),rpt(3,nrpt),spqpt(3,nqpt) 3611 real(dp),intent(in) :: wghatm(natom,natom,nrpt) 3612 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 3613 real(dp),intent(out) :: dynmat(2,3,natom,3,natom,nqpt) 3614 3615 !Local variables ------------------------- 3616 !scalars 3617 integer :: ia,ib,iqpt,irpt,mu,nu,cnt,my_rank,nprocs, ierr 3618 real(dp) :: facti,factr,im,kr,re 3619 !real(dp) : w(2, natom, natom) 3620 !arrays 3621 real(dp) :: kk(3) 3622 3623 ! ********************************************************************* 3624 3625 my_rank = xmpi_comm_rank(comm); nprocs = xmpi_comm_size(comm) 3626 dynmat = zero; cnt = 0 3627 3628 ! MG: This is an hotspot. I don'tknow whether one should rewrite with BLAS1 dot or not. 3629 ! Note, however, that simply removing the check on the weights inside the loop over atoms. 3630 ! leads to a non-negligible speedup with intel (~30% if dipdip -1 is used) 3631 do iqpt=1,nqpt 3632 3633 ! Calculation of the k coordinates in Normalized Reciprocal coordinates 3634 kk(1)=spqpt(1,iqpt)*gprim(1,1)+spqpt(2,iqpt)* gprim(1,2)+spqpt(3,iqpt)*gprim(1,3) 3635 kk(2)=spqpt(1,iqpt)*gprim(2,1)+spqpt(2,iqpt)* gprim(2,2)+spqpt(3,iqpt)*gprim(2,3) 3636 kk(3)=spqpt(1,iqpt)*gprim(3,1)+spqpt(2,iqpt)* gprim(3,2)+spqpt(3,iqpt)*gprim(3,3) 3637 3638 do irpt=1,nrpt 3639 cnt = cnt + 1; if (mod(cnt, nprocs) /= my_rank) cycle ! MPI parallelism. 3640 3641 ! k.R 3642 kr = kk(1)*rpt(1,irpt)+kk(2)*rpt(2,irpt)+kk(3)*rpt(3,irpt) 3643 ! Get phase factor 3644 re = cos(two_pi*kr); im = sin(two_pi*kr) 3645 3646 ! Inner loop on atoms and directions 3647 do ib=1,natom 3648 do ia=1,natom 3649 !if (abs(wghatm(ia,ib,irpt)) > tol10) then ! Commented by MG 3650 factr = re * wghatm(ia,ib,irpt) 3651 facti = im * wghatm(ia,ib,irpt) 3652 do nu=1,3 3653 do mu=1,3 3654 ! Real and imaginary part of the dynamical matrices 3655 ! Atmfrc should be real 3656 dynmat(1,mu,ia,nu,ib,iqpt) = dynmat(1,mu,ia,nu,ib,iqpt) + factr * atmfrc(mu,ia,nu,ib,irpt) 3657 dynmat(2,mu,ia,nu,ib,iqpt) = dynmat(2,mu,ia,nu,ib,iqpt) + facti * atmfrc(mu,ia,nu,ib,irpt) 3658 end do 3659 end do 3660 !end if 3661 end do 3662 end do 3663 3664 ! MG: New version: I don't know if it's faster. 3665 !w(1,:,:) = re * wghatm(:,:,irpt) 3666 !w(2,:,:) = im * wghatm(:,:,irpt) 3667 !do ib=1,natom 3668 ! do nu=1,3 3669 ! do ia=1,natom 3670 ! do mu=1,3 3671 ! ! Real and imaginary part of the dynamical matrices 3672 ! ! Atmfrc should be real 3673 ! dynmat(1,mu,ia,nu,ib,iqpt) = dynmat(1,mu,ia,nu,ib,iqpt) + w(1,ia,ib) * atmfrc(mu,ia,nu,ib,irpt) 3674 ! dynmat(2,mu,ia,nu,ib,iqpt) = dynmat(2,mu,ia,nu,ib,iqpt) + w(2,ia,ib) * atmfrc(mu,ia,nu,ib,irpt) 3675 ! end do 3676 ! end do 3677 ! end do 3678 !end do 3679 3680 end do 3681 end do 3682 3683 if (nprocs > 1) call xmpi_sum(dynmat, comm, ierr) 3684 3685 end subroutine ftifc_r2q
m_dynmat/get_bigbox_and_weights [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
get_bigbox_and_weights
FUNCTION
Compute the Big Box containing the R points in the cartesian real space needed to Fourier Transform 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.) natom= Number of atoms nqbz= Number of q-points in BZ. ngqpt(3)= Numbers used to generate the q points to sample the Brillouin zone using an homogeneous grid nqshift= number of shifts in q-mesh qshift(3, nqshift) = Q-mesh shifts rprim(3,3)= Normalized coordinates in real space. rprimd, gprimd rcan(3,natom) = Atomic position in canonical coordinates cutmode=Define the cutoff used to filter the output R-points according to their weights. 0 --> No cutoff (mainly for debugging) 1 --> Include only those R-points for which sum(abs(wg(:,:,irpt)) < tol20 This is the approach used for the dynamical matrix. 2 --> Include only those R-points for which the trace over iatom of abs(wg(iat,iat,irpt)) < tol20 This option is used for objects that depend on a single atomic index. comm= MPI communicator
OUTPUT
nrpt= Total Number of R points in the Big Box cell(3,nrpt) Give the index of the the cell and irpt rpt(3,nrpt)= Canonical coordinates of the R points in the unit cell. These coordinates are normalized (=> * acell(3)!!) r_inscribed_sphere wghatm(natom,natom,nrpt)= Weights associated to a pair of atoms and to a R vector
SOURCE
2672 subroutine get_bigbox_and_weights(brav, natom, nqbz, ngqpt, nqshift, qshift, rprim, rprimd, gprim, rcan, & 2673 cutmode, nrpt, rpt, cell, wghatm, r_inscribed_sphere, comm) 2674 2675 !Arguments ------------------------------- 2676 !scalars 2677 integer,intent(in) :: brav, natom, nqbz, nqshift, cutmode, comm 2678 integer,intent(out) :: nrpt 2679 real(dp),intent(out) :: r_inscribed_sphere 2680 !arrays 2681 integer,intent(in) :: ngqpt(3) 2682 real(dp),intent(in) :: gprim(3,3),rprim(3,3),rprimd(3,3), rcan(3, natom) 2683 real(dp),intent(in) :: qshift(3, nqshift) 2684 integer,allocatable,intent(out) :: cell(:,:) 2685 real(dp),allocatable,intent(out) :: rpt(:,:), wghatm(:,:,:) 2686 2687 !Local variables ------------------------- 2688 !scalars 2689 integer :: my_ierr, ierr, ii, irpt, all_nrpt 2690 real(dp) :: toldist 2691 integer :: ngqpt9(9) 2692 character(len=500*4) :: msg 2693 !arrays 2694 integer,allocatable :: all_cell(:,:) 2695 real(dp),allocatable :: all_rpt(:,:), all_wghatm(:,:,:) 2696 2697 ! ********************************************************************* 2698 2699 ABI_CHECK(any(cutmode == [0, 1, 2]), "cutmode should be in [0, 1, 2]") 2700 2701 ! Create the Big Box of R vectors in real space and compute the number of points (cells) in real space 2702 call make_bigbox(brav, all_cell, ngqpt, nqshift, rprim, all_nrpt, all_rpt) 2703 2704 ! Weights associated to these R points and to atomic pairs 2705 ABI_MALLOC(all_wghatm, (natom, natom, all_nrpt)) 2706 2707 ! HM: this tolerance is highly dependent on the compilation/architecture 2708 ! numeric errors in the DDB text file. Try a few tolerances and check whether all the weights are found. 2709 ngqpt9 = 0; ngqpt9(1:3) = ngqpt(1:3) 2710 toldist = tol8 2711 do while (toldist <= tol6) 2712 ! Note ngqpt(9) with intent(inout)! 2713 call wght9(brav, gprim, natom, ngqpt9, nqbz, nqshift, all_nrpt, qshift, rcan, & 2714 all_rpt, rprimd, toldist, r_inscribed_sphere, all_wghatm, my_ierr) 2715 call xmpi_max(my_ierr, ierr, comm, ii) 2716 if (ierr > 0) toldist = toldist * 10 2717 if (ierr == 0) exit 2718 end do 2719 2720 if (ierr > 0) then 2721 write(msg, '(3a,es14.4,2a,i0, 14a)' ) & 2722 'The sum of the weight is not equal to nqpt.',ch10,& 2723 'The sum of the weights is: ',sum(all_wghatm),ch10,& 2724 'The number of q points is: ',nqbz, ch10, & 2725 'This might have several sources.',ch10,& 2726 'If toldist is larger than 1.0e-8, the atom positions might be loose.',ch10,& 2727 'and the q point weights not computed properly.',ch10,& 2728 'Action: make input atomic positions more symmetric.',ch10,& 2729 'Otherwise, you might increase "buffer" in m_dynmat.F90 see bigbx9 subroutine and recompile.',ch10,& 2730 'Actually, this can also happen when ngqpt is 0 0 0,',ch10,& 2731 'if abs(brav) /= 1, in this case you should change brav to 1. If brav is already set to 1 (default) try -1.' 2732 ABI_ERROR(msg) 2733 end if 2734 2735 ! Only conserve the necessary points in rpt. 2736 nrpt = 0 2737 do irpt=1,all_nrpt 2738 if (filterw(all_wghatm(:,:,irpt))) cycle 2739 nrpt = nrpt + 1 2740 end do 2741 2742 ! Allocate output arrays and transfer data. 2743 ABI_MALLOC(rpt, (3, nrpt)) 2744 ABI_MALLOC(cell, (3, nrpt)) 2745 ABI_MALLOC(wghatm, (natom, natom, nrpt)) 2746 2747 ii = 0 2748 do irpt=1,all_nrpt 2749 if (filterw(all_wghatm(:,:,irpt))) cycle 2750 ii = ii + 1 2751 rpt(:, ii) = all_rpt(:,irpt) 2752 wghatm(:,:,ii) = all_wghatm(:,:,irpt) 2753 cell(:,ii) = all_cell(:,irpt) 2754 end do 2755 2756 ABI_FREE(all_rpt) 2757 ABI_FREE(all_wghatm) 2758 ABI_FREE(all_cell) 2759 2760 contains 2761 2762 logical pure function filterw(wg) 2763 2764 real(dp),intent(in) :: wg(natom,natom) 2765 integer :: iat 2766 real(dp) :: trace 2767 2768 select case (cutmode) 2769 case (1) 2770 filterw = sum(abs(wg)) < tol20 2771 case (2) 2772 trace = zero 2773 do iat=1,natom 2774 trace = trace + abs(wg(iat,iat)) 2775 end do 2776 filterw = trace < tol20 2777 case default 2778 filterw = .False. 2779 end select 2780 2781 end function filterw 2782 2783 end subroutine get_bigbox_and_weights
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 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 comm=MPI communicator. [dipquad] = if 1, atmfrc has been build without dipole-quadrupole part [quadquad] = if 1, atmfrc has been build without quadrupole-quadrupole part
OUTPUT
d2cart(2,3,mpert,3,mpert)=dynamical matrix obtained for the wavevector qpt (normalized using qphnrm)
SOURCE
5510 subroutine gtdyn9(acell,atmfrc,dielt,dipdip,dyewq0,d2cart,gmet,gprim,mpert,natom,& 5511 nrpt,qphnrm,qpt,rmet,rprim,rpt,trans,ucvol,wghatm,xred,zeff,qdrp_cart,ewald_option,comm,& 5512 dipquad,quadquad) ! optional 5513 5514 !Arguments ------------------------------- 5515 !scalars 5516 integer,intent(in) :: dipdip,mpert,natom,nrpt,ewald_option,comm 5517 real(dp),intent(in) :: qphnrm,ucvol 5518 integer,optional,intent(in) :: dipquad, quadquad 5519 !arrays 5520 real(dp),intent(in) :: acell(3),dielt(3,3),gmet(3,3),gprim(3,3),qpt(3) 5521 real(dp),intent(in) :: rmet(3,3),rprim(3,3),rpt(3,nrpt) 5522 real(dp),intent(in) :: trans(3,natom),wghatm(natom,natom,nrpt),xred(3,natom) 5523 real(dp),intent(in) :: zeff(3,3,natom) 5524 real(dp),intent(in) :: qdrp_cart(3,3,3,natom) 5525 real(dp),intent(in) :: atmfrc(3,natom,3,natom,nrpt) 5526 real(dp),intent(in) :: dyewq0(3,3,natom) 5527 real(dp),intent(out) :: d2cart(2,3,mpert,3,mpert) 5528 5529 !Local variables ------------------------- 5530 !scalars 5531 integer,parameter :: nqpt1 = 1, option2 = 2, sumg0 = 0, plus1 = 1, iqpt1 = 1 5532 integer :: i1, i2, ib, nsize, dipquad_, quadquad_ 5533 !arrays 5534 real(dp) :: qphon(3) !, tsec(2) 5535 real(dp),allocatable :: dq(:,:,:,:,:),dyew(:,:,:,:,:) 5536 5537 ! ********************************************************************* 5538 5539 ! Keep track of time spent in gtdyn9 5540 !call timab(1750, 1, tsec) 5541 5542 ABI_MALLOC(dq,(2,3,natom,3,natom)) 5543 5544 ! Define quadrupolar options 5545 dipquad_=0; if(present(dipquad)) dipquad_=dipquad 5546 quadquad_=0; if(present(quadquad)) quadquad_=quadquad 5547 5548 ! Get the normalized wavevector 5549 if(abs(qphnrm)<1.0d-7)then 5550 qphon(1:3)=zero 5551 else 5552 qphon(1:3)=qpt(1:3)/qphnrm 5553 end if 5554 5555 ! Generate the analytical part from the interatomic forces 5556 call ftifc_r2q(atmfrc, dq, gprim, natom, nqpt1, nrpt, rpt, qphon, wghatm, comm) 5557 5558 ! The analytical dynamical matrix dq has been generated 5559 ! in the normalized canonical coordinate system. 5560 ! Now, the phase is modified, in order to recover the usual (xred) coordinate of atoms. 5561 call dymfz9(dq,natom,nqpt1,gprim,option2,qphon,trans) 5562 5563 if (dipdip==1.or.dipquad_==1.or.quadquad_==1) then 5564 ! Add the non-analytical part 5565 ! Compute dyew(2,3,natom,3,natom)= Ewald part of the dynamical matrix, 5566 ! second energy derivative wrt xred(3,natom) in Hartrees (Denoted A-bar in the notes) 5567 ABI_MALLOC(dyew,(2,3,natom,3,natom)) 5568 5569 call ewald9(acell,dielt,dyew,gmet,gprim,natom,qphon,rmet,rprim,sumg0,ucvol,xred,zeff,& 5570 qdrp_cart,option=ewald_option,dipquad=dipquad_,quadquad=quadquad_) 5571 5572 call q0dy3_apply(natom,dyewq0,dyew) 5573 call nanal9(dyew,dq,iqpt1,natom,nqpt1,plus1) 5574 5575 ABI_FREE(dyew) 5576 end if 5577 5578 ! Copy the dynamical matrix in the proper location 5579 ! First zero all the elements 5580 nsize=2*(3*mpert)**2 5581 d2cart = zero 5582 5583 ! Copy the elements from dq to d2cart 5584 d2cart(:,:,1:natom,:,1:natom)=dq(:,:,1:natom,:,1:natom) 5585 5586 ! In case we have the gamma point, 5587 if(qphon(1)**2+qphon(2)**2+qphon(3)**2<1.d-14)then 5588 ! Copy the effective charge and dielectric constant in the final array 5589 do i1=1,3 5590 do i2=1,3 5591 d2cart(1,i1,natom+2,i2,natom+2)=dielt(i1,i2) 5592 do ib=1,natom 5593 d2cart(1,i1,natom+2,i2,ib)=zeff(i1,i2,ib) 5594 d2cart(1,i2,ib,i1,natom+2)=zeff(i1,i2,ib) 5595 end do 5596 end do 5597 end do 5598 end if 5599 5600 ABI_FREE(dq) 5601 5602 !call timab(1750, 2, tsec) 5603 5604 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
SOURCE
3797 subroutine ifclo9(ifccar,ifcloc,vect1,vect2,vect3) 3798 3799 !Arguments ------------------------------- 3800 !arrays 3801 real(dp),intent(in) :: ifccar(3,3),vect1(3),vect2(3),vect3(3) 3802 real(dp),intent(out) :: ifcloc(3,3) 3803 3804 !Local variables ------------------------- 3805 !scalars 3806 integer :: ii,jj 3807 !arrays 3808 real(dp) :: work(3,3) 3809 3810 ! ********************************************************************* 3811 3812 do jj=1,3 3813 do ii=1,3 3814 work(jj,ii)=zero 3815 end do 3816 do ii=1,3 3817 work(jj,1)=work(jj,1)+ifccar(jj,ii)*vect1(ii) 3818 work(jj,2)=work(jj,2)+ifccar(jj,ii)*vect2(ii) 3819 work(jj,3)=work(jj,3)+ifccar(jj,ii)*vect3(ii) 3820 end do 3821 end do 3822 3823 do jj=1,3 3824 do ii=1,3 3825 ifcloc(ii,jj)=zero 3826 end do 3827 do ii=1,3 3828 ifcloc(1,jj)=ifcloc(1,jj)+vect1(ii)*work(ii,jj) 3829 ifcloc(2,jj)=ifcloc(2,jj)+vect2(ii)*work(ii,jj) 3830 ifcloc(3,jj)=ifcloc(3,jj)+vect3(ii)*work(ii,jj) 3831 end do 3832 end do 3833 3834 end subroutine ifclo9
m_dynmat/make_bigbox [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
make_bigbox
FUNCTION
Helper functions to faciliate 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(3,nrpt)= integer coordinates of the cells (R points) in the rprim basis nprt= Number of cells (R points) in the Big Box rpt(3,mrpt)= canonical coordinates of the cells (R points) These coordinates are normalized (=> * acell(3)!!) The array is allocated here with the proper dimension. Client code is responsible for the deallocation.
SOURCE
2817 subroutine make_bigbox(brav, cell, ngqpt, nqshft, rprim, nrpt, rpt) 2818 2819 !Arguments ------------------------------- 2820 !scalars 2821 integer,intent(in) :: brav,nqshft 2822 integer,intent(out) :: nrpt 2823 !arrays 2824 integer,intent(in) :: ngqpt(3) 2825 real(dp),intent(in) :: rprim(3,3) 2826 real(dp),allocatable,intent(out) :: rpt(:,:) 2827 integer,allocatable,intent(out) :: cell(:,:) 2828 2829 !Local variables ------------------------- 2830 !scalars 2831 integer :: choice,mrpt 2832 !arrays 2833 real(dp) :: dummy_rpt(3,1) 2834 integer:: dummy_cell(1,3) 2835 2836 ! ********************************************************************* 2837 2838 ! Compute the number of points (cells) in real space 2839 choice=0 2840 call bigbx9(brav,dummy_cell,choice,1,ngqpt,nqshft,mrpt,rprim,dummy_rpt) 2841 2842 ! Now we can allocate and calculate the points and the weights. 2843 nrpt = mrpt 2844 ABI_MALLOC(rpt,(3,nrpt)) 2845 ABI_MALLOC(cell,(3,nrpt)) 2846 2847 choice=1 2848 call bigbx9(brav,cell,choice,mrpt,ngqpt,nqshft,nrpt,rprim,rpt) 2849 2850 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,...) [herm_opt]= 1 to hermitianize mat (default) 0 if no symmetrization should be performed
SIDE EFFECTS
mat(2*3*natom*3*natom)=Multiplies by atomic masses in output.
SOURCE
6218 subroutine massmult_and_breaksym(natom, ntypat, typat, amu, mat, & 6219 herm_opt) ! optional 6220 6221 !Arguments ------------------------------- 6222 !scalars 6223 integer,intent(in) :: natom,ntypat 6224 integer,optional,intent(in) :: herm_opt 6225 !arrays 6226 integer,intent(in) :: typat(natom) 6227 real(dp),intent(in) :: amu(ntypat) 6228 real(dp),intent(inout) :: mat(2*3*natom*3*natom) 6229 6230 !Local variables ------------------------- 6231 !scalars 6232 integer :: i1,i2,idir1,idir2,index,ipert1,ipert2, herm_opt__ 6233 real(dp),parameter :: break_symm=1.0d-12 6234 !real(dp),parameter :: break_symm=zero 6235 real(dp) :: fac 6236 !arrays 6237 real(dp) :: nearidentity(3,3) 6238 6239 ! ********************************************************************* 6240 6241 herm_opt__ = 1; if (present(herm_opt)) herm_opt__ = herm_opt 6242 6243 ! This slight breaking of the symmetry allows the results to be more portable between machines 6244 nearidentity(:,:)=one 6245 nearidentity(1,1)=one+break_symm 6246 nearidentity(3,3)=one-break_symm 6247 6248 ! Include the masses in the dynamical matrix 6249 do ipert1=1,natom 6250 do ipert2=1,natom 6251 fac=1.0_dp/sqrt(amu(typat(ipert1))*amu(typat(ipert2)))/amu_emass 6252 do idir1=1,3 6253 do idir2=1,3 6254 i1=idir1+(ipert1-1)*3 6255 i2=idir2+(ipert2-1)*3 6256 index=i1+3*natom*(i2-1) 6257 mat(2*index-1)=mat(2*index-1)*fac*nearidentity(idir1,idir2) 6258 mat(2*index )=mat(2*index )*fac*nearidentity(idir1,idir2) 6259 ! This is to break slightly the translation invariance, and make the automatic tests more portable 6260 if(ipert1==ipert2 .and. idir1==idir2)then 6261 mat(2*index-1)=mat(2*index-1)+break_symm*natom/amu_emass/idir1*0.01_dp 6262 end if 6263 end do 6264 end do 6265 end do 6266 end do 6267 6268 ! Make the dynamical matrix hermitian 6269 if (herm_opt__ == 1) call mkherm(mat,3*natom) 6270 6271 end subroutine massmult_and_breaksym
m_dynmat/massmult_and_breaksym_cplx [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
mult_masses_and_break_symms_cplx
FUNCTION
Similar to massmult_and_breaksym, the only difference is that it receives complex array.
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. For plus=0, see Eq.(76) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], get the left hand side.
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
SOURCE
5415 subroutine nanal9(dyew,dynmat,iqpt,natom,nqpt,plus) 5416 5417 !Arguments ------------------------------- 5418 !scalars 5419 integer,intent(in) :: iqpt,natom,nqpt,plus 5420 !arrays 5421 real(dp),intent(in) :: dyew(2,3,natom,3,natom) 5422 real(dp),intent(inout) :: dynmat(2,3,natom,3,natom,nqpt) 5423 5424 !Local variables ------------------------- 5425 !scalars 5426 integer :: ia,ib,mu,nu 5427 character(len=500) :: msg 5428 5429 ! ********************************************************************* 5430 5431 if (plus==0) then 5432 5433 do ia=1,natom 5434 do ib=1,natom 5435 do mu=1,3 5436 do nu=1,3 5437 ! The following four lines are OK 5438 dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) - dyew(1,mu,ia,nu,ib) 5439 dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) - dyew(2,mu,ia,nu,ib) 5440 end do 5441 end do 5442 end do 5443 end do 5444 5445 else if (plus==1) then 5446 do ia=1,natom 5447 do ib=1,natom 5448 do mu=1,3 5449 do nu=1,3 5450 dynmat(1,mu,ia,nu,ib,iqpt)=dynmat(1,mu,ia,nu,ib,iqpt) + dyew(1,mu,ia,nu,ib) 5451 dynmat(2,mu,ia,nu,ib,iqpt)=dynmat(2,mu,ia,nu,ib,iqpt) + dyew(2,mu,ia,nu,ib) 5452 end do 5453 end do 5454 end do 5455 end do 5456 5457 else 5458 write(msg,'(3a,i0,a)' )& 5459 'The argument "plus" must be equal to 0 or 1.',ch10,& 5460 'The value ',plus,' is not available.' 5461 ABI_BUG(msg) 5462 end if 5463 5464 end subroutine nanal9
m_dynmat/phdispl_from_eigvec [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
phdispl_from_eigvec
FUNCTION
Phonon displacements in cart coords from eigenvectors
INPUTS
natom: number of atoms in unit cell ntypat=number of atom types typat(natom)=integer label of each type of atom (1,2,...) amu(ntypat)=mass of the atoms (atomic mass unit) matrix (diagonal in the atoms) eigvec(2*3*natom*3*natom)= eigenvectors of the dynamical matrix in cartesian coordinates.
OUTPUT
displ(2*3*natom*3*natom)=displacements of atoms in cartesian coordinates.
SOURCE
5946 pure subroutine phdispl_from_eigvec(natom, ntypat, typat, amu, eigvec, displ) 5947 5948 !Arguments ------------------------------- 5949 !scalars 5950 integer,intent(in) :: natom, ntypat 5951 !arrays 5952 integer,intent(in) :: typat(natom) 5953 real(dp),intent(in) :: amu(ntypat) 5954 real(dp),intent(in) :: eigvec(2*3*natom*3*natom) 5955 real(dp),intent(out) :: displ(2*3*natom*3*natom) 5956 5957 !Local variables ------------------------- 5958 !scalars 5959 integer :: i1,idir1,imode,ipert1, index 5960 5961 ! ********************************************************************* 5962 5963 do imode=1,3*natom 5964 5965 do idir1=1,3 5966 do ipert1=1,natom 5967 i1=idir1+(ipert1-1)*3 5968 index=i1+3*natom*(imode-1) 5969 displ(2*index-1)=eigvec(2*index-1) / sqrt(amu(typat(ipert1))*amu_emass) 5970 displ(2*index )=eigvec(2*index ) / sqrt(amu(typat(ipert1))*amu_emass) 5971 end do 5972 end do 5973 5974 end do 5975 5976 end subroutine phdispl_from_eigvec
m_dynmat/pheigvec_normalize [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
pheigvec_normalize
FUNCTION
Normalize input eigenvectors in cartesian coordinates
INPUTS
natom: number of atoms in unit cell
SIDE EFFECTS
eigvec(2*3*natom*3*natom)=in output the normalized eigenvectors in cartesian coordinates.
SOURCE
5884 pure subroutine pheigvec_normalize(natom, eigvec) 5885 5886 !Arguments ------------------------------- 5887 !scalars 5888 integer,intent(in) :: natom 5889 !arrays 5890 real(dp),intent(inout) :: eigvec(2*3*natom*3*natom) 5891 5892 !Local variables ------------------------- 5893 !scalars 5894 integer :: i1,idir1,imode,ipert1,index 5895 real(dp) :: norm 5896 5897 ! ********************************************************************* 5898 5899 do imode=1,3*natom 5900 5901 norm=zero 5902 do idir1=1,3 5903 do ipert1=1,natom 5904 i1=idir1+(ipert1-1)*3 5905 index=i1+3*natom*(imode-1) 5906 norm=norm+eigvec(2*index-1)**2+eigvec(2*index)**2 5907 end do 5908 end do 5909 norm=sqrt(norm) 5910 5911 do idir1=1,3 5912 do ipert1=1,natom 5913 i1=idir1+(ipert1-1)*3 5914 index=i1+3*natom*(imode-1) 5915 eigvec(2*index-1)=eigvec(2*index-1)/norm 5916 eigvec(2*index)=eigvec(2*index)/norm 5917 end do 5918 end do 5919 5920 end do 5921 5922 end subroutine pheigvec_normalize
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 See Eq.(71) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], get the left hand side.
INPUTS
dyewq0(3,3,natom) = part needed to correct the dynamical matrix for atom self-interaction. natom= number of atom in the unit cell
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
SOURCE
1879 subroutine q0dy3_apply(natom,dyewq0,dyew) 1880 1881 !Arguments ------------------------------- 1882 !scalars 1883 integer,intent(in) :: natom 1884 !arrays 1885 real(dp),intent(in) :: dyewq0(3,3,natom) 1886 real(dp),intent(inout) :: dyew(2,3,natom,3,natom) 1887 1888 !Local variables ------------------------- 1889 !scalars 1890 integer :: ia,mu,nu 1891 1892 ! ********************************************************************* 1893 1894 do mu=1,3 1895 do nu=1,3 1896 do ia=1,natom 1897 dyew(1,mu,ia,nu,ia)=dyew(1,mu,ia,nu,ia)-dyewq0(mu,nu,ia) 1898 end do 1899 end do 1900 end do 1901 1902 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 See Eq.(71) in Gonze&Lee PRB 55, 10355 (1997) [[cite:Gonze1997a]], the sum over \kappa"
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
SOURCE
1949 subroutine q0dy3_calc(natom,dyewq0,dyew,option) 1950 1951 !Arguments ------------------------------- 1952 !scalars 1953 integer,intent(in) :: natom,option 1954 !arrays 1955 real(dp),intent(in) :: dyew(2,3,natom,3,natom) 1956 real(dp),intent(out) :: dyewq0(3,3,natom) 1957 1958 !Local variables ------------------------- 1959 !scalars 1960 integer :: ia,ib,mu,nu 1961 character(len=500) :: msg 1962 1963 ! ********************************************************************* 1964 1965 if(option==1.or.option==2)then 1966 do mu=1,3 1967 do nu=1,3 1968 do ia=1,natom 1969 dyewq0(mu,nu,ia)=zero 1970 do ib=1,natom 1971 dyewq0(mu,nu,ia)=dyewq0(mu,nu,ia)+dyew(1,mu,ia,nu,ib) 1972 end do 1973 end do 1974 end do 1975 end do 1976 else 1977 write (msg, '(3a)')& 1978 & 'option should be 1 or 2.',ch10,& 1979 & 'action: correct calling routine' 1980 ABI_BUG(msg) 1981 end if 1982 1983 if(option==2)then 1984 do ia=1,natom 1985 do mu=1,3 1986 do nu=mu,3 1987 dyewq0(mu,nu,ia)=(dyewq0(mu,nu,ia)+dyewq0(nu,mu,ia))/2 1988 dyewq0(nu,mu,ia)=dyewq0(mu,nu,ia) 1989 end do 1990 end do 1991 end do 1992 end if 1993 1994 end subroutine q0dy3_calc
m_dynmat/sylwtens [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
sylwtens
FUNCTION
Determines the set of irreductible elements of the non-linear optical susceptibility and Raman tensors
INPUTS
has_strain = if .true. i2pert includes strain perturbation 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 reduced space) symrel(3,3,nsym)=3x3 matrices of the group symmetries (real reduced space) symrel_cart(3,3,nsym)=3x3 matrices of the group symmetries (real cartesian 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,m_nonlinear,m_respfn_driver
CHILDREN
SOURCE
4980 !subroutine sylwtens(indsym,mpert,natom,nsym,rfpert,symrec,symrel,symrel_cart) 4981 subroutine sylwtens(indsym,mpert,natom,nsym,rfpert,symrec,symrel) 4982 4983 !Arguments ------------------------------- 4984 !scalars 4985 integer,intent(in) :: mpert,natom,nsym 4986 !arrays 4987 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4988 integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert) 4989 ! real(dp),intent(in) :: symrel_cart(3,3,nsym) 4990 4991 !Local variables ------------------------- 4992 !scalars 4993 integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_ 4994 ! integer :: i2dir_a,i2dir_b 4995 integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2 4996 integer :: ipesy3,isym 4997 ! integer :: istr,idisy2_a,idisy2_b 4998 logical :: is_strain 4999 ! real(dp) :: flag_dp 5000 !arrays 5001 ! integer,save :: idx(18)=(/1,1,2,2,3,3,3,2,3,1,2,1,2,3,1,3,1,2/) 5002 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 5003 integer,allocatable :: pertsy(:,:,:,:,:,:) 5004 5005 !*********************************************************************** 5006 5007 ABI_MALLOC(pertsy,(3,mpert,3,mpert,3,mpert)) 5008 pertsy(:,:,:,:,:,:) = 0 5009 5010 !Loop over perturbations 5011 5012 do i1pert_ = 1, mpert 5013 do i2pert_ = 1, mpert 5014 is_strain=.false. 5015 do i3pert_ = 1, mpert 5016 5017 do i1dir_ = 1, 3 5018 do i2dir_ = 1, 3 5019 do i3dir_ = 1, 3 5020 5021 i1pert = (mpert - i1pert_ + 1) 5022 if (i1pert <= natom) i1pert = natom + 1 - i1pert 5023 i2pert = (mpert - i2pert_ + 1) 5024 if (i2pert <= natom) i2pert = natom + 1 - i2pert 5025 i3pert = (mpert - i3pert_ + 1) 5026 if (i3pert <= natom) i3pert = natom + 1 - i3pert 5027 5028 if (i1pert <= natom) then 5029 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 5030 else if (i2pert <= natom) then 5031 i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_ 5032 else if (i3pert <= natom) then 5033 i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_ 5034 else 5035 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 5036 end if 5037 5038 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then 5039 5040 ! Loop over all symmetries 5041 5042 flag = 0 5043 do isym = 1, nsym 5044 5045 found = 1 5046 5047 ! Select the symmetric element of i1pert,i2pert,i3pert 5048 5049 if (i1pert <= natom) then 5050 ipesy1 = indsym(4,isym,i1pert) 5051 sym1(:,:) = symrec(:,:,isym) 5052 else if (i1pert == natom + 2) then 5053 ipesy1 = i1pert 5054 sym1(:,:) = symrel(:,:,isym) 5055 else 5056 found = 0 5057 end if 5058 5059 if (i2pert <= natom) then 5060 ipesy2 = indsym(4,isym,i2pert) 5061 sym2(:,:) = symrec(:,:,isym) 5062 else if (i2pert == natom + 2) then 5063 ipesy2 = i2pert 5064 sym2(:,:) = symrel(:,:,isym) 5065 else if (i2pert == natom + 3.or. i2pert == natom + 4) then 5066 ! !TODO: Symmetries on strain perturbation do not work yet. 5067 found = 0 5068 is_strain=.true. 5069 ! 5070 ! ipesy2 = i2pert 5071 ! sym2(:,:) = NINT(symrel_cart(:,:,isym)) 5072 ! if (i2pert==natom+3) istr=i2dir 5073 ! if (i2pert==natom+4) istr=3+i2dir 5074 ! i2dir_a=idx(2*istr-1); i2dir_b=idx(2*istr) 5075 else 5076 found = 0 5077 end if 5078 5079 if (i3pert == natom + 8) then 5080 ipesy3 = i3pert 5081 sym3(:,:) = symrel(:,:,isym) 5082 else 5083 found = 0 5084 end if 5085 5086 5087 ! See if the symmetric element is available and check if some 5088 ! of the elements may be zero. In the latter case, they do not need 5089 ! to be computed. 5090 5091 if (.not.is_strain) then 5092 if ((flag /= -1).and.& 5093 & (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then 5094 flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir) 5095 end if 5096 5097 do idisy1 = 1, 3 5098 do idisy2 = 1, 3 5099 do idisy3 = 1, 3 5100 5101 if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.& 5102 & (sym3(i3dir,idisy3) /= 0)) then 5103 if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then 5104 found = 0 5105 ! exit ! exit loop over symmetries 5106 end if 5107 end if 5108 5109 5110 if ((flag == -1).and.& 5111 & ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then 5112 if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.& 5113 & (sym3(i3dir,idisy3)/=0)) then 5114 flag = 0 5115 end if 5116 end if 5117 5118 end do 5119 end do 5120 end do 5121 ! else 5122 ! if ((flag_dp /= -1).and.& 5123 !& (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then 5124 ! flag = sym1(i1dir,i1dir)*sym2(i2dir_a,i2dir_a)* & 5125 ! & sym2(i2dir_b,i2dir_b)*sym3(i3dir,i3dir) 5126 ! end if 5127 ! 5128 ! do idisy1 = 1, 3 5129 ! do idisy2 = 1, 3 5130 ! if (ipesy2==natom+3) istr=idisy2 5131 ! if (ipesy2==natom+4) istr=3+idisy2 5132 ! idisy2_a=idx(2*istr-1); idisy2_b=idx(2*istr) 5133 ! do idisy3 = 1, 3 5134 ! 5135 ! if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir_a,idisy2_a) /= 0).and.& 5136 !& (sym2(i2dir_b,idisy2_b) /= 0).and.(sym3(i3dir,idisy3) /= 0)) then 5137 ! if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then 5138 ! found = 0 5139 !! exit ! exit loop over symmetries 5140 ! end if 5141 ! end if 5142 ! 5143 ! 5144 ! if ((flag == -1).and.& 5145 !& ((idisy1/=i1dir).or.(idisy2_a/=i2dir_a).or.(idisy2_b/=i2dir_b).or.(idisy3/=i3dir))) then 5146 ! if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir_a,idisy2_a)/=0).and.& 5147 !& (sym2(i2dir_b,idisy2_b)/=0).and.(sym3(i3dir,idisy3)/=0)) then 5148 ! flag = 0 5149 ! end if 5150 ! end if 5151 ! 5152 ! end do 5153 ! end do 5154 ! end do 5155 end if 5156 5157 if (found == 1) then 5158 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1 5159 end if 5160 5161 ! In case a symmetry operation only changes the sign of an 5162 ! element, this element has to be equal to zero 5163 5164 if (flag == -1) then 5165 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2 5166 exit 5167 end if 5168 5169 end do ! close loop on symmetries 5170 5171 ! If the elemetn i1pert,i2pert,i3pert is not symmetric 5172 ! to a basis element, it is a basis element 5173 5174 if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then 5175 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 5176 end if 5177 5178 end if ! rfpert /= 0 5179 5180 end do ! close loop over perturbations 5181 end do 5182 end do 5183 end do 5184 end do 5185 end do 5186 5187 !Now, take into account the permutation of (i1pert,i1dir) 5188 !and (i2pert,i2dir) 5189 5190 do i1pert = 1, mpert 5191 do i2pert = 1, mpert 5192 do i3pert = 1, mpert 5193 5194 do i1dir = 1, 3 5195 do i2dir = 1, 3 5196 do i3dir = 1, 3 5197 5198 if ((i1pert /= i2pert).or.(i1dir /= i2dir)) then 5199 5200 if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.& 5201 (pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) == 1)) then 5202 pertsy(i2dir,i2pert,i1dir,i1pert,i3dir,i3pert) = -1 5203 end if 5204 5205 end if 5206 5207 end do 5208 end do 5209 end do 5210 5211 end do 5212 end do 5213 end do 5214 5215 rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:) 5216 5217 ABI_FREE(pertsy) 5218 5219 end subroutine sylwtens
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) [[cite:Hendrikse1995]]
TODO
A full description of the equations should be included
SOURCE
2033 subroutine symdyma(dmati,indsym,natom,nsym,qptn,rprimd,symrel,symafm) 2034 2035 !Arguments ------------------------------- 2036 !scalars 2037 integer,intent(in) :: natom,nsym 2038 !arrays 2039 integer,intent(in) :: indsym(4,nsym,natom),symrel(3,3,nsym) 2040 integer,intent(in) :: symafm(nsym) 2041 real(dp),intent(in) :: qptn(3),rprimd(3,3) 2042 real(dp),intent(inout) :: dmati(2*3*natom*3*natom) 2043 2044 !Local variables ------------------------- 2045 !scalars 2046 integer :: i1,i2,iat,idir,ii,index,isgn,isym,itirev,jat,jdir,jj,kk,ll 2047 integer :: niat,njat,timrev 2048 real(dp) :: arg1,arg2,dmint,im,re,sumi,sumr 2049 !arrays 2050 integer :: indij(natom,natom),symq(4,2,nsym),symrec(3,3,nsym) 2051 real(dp) :: TqR(3,3),TqS_(3,3),dynmat(2,3,natom,3,natom) 2052 real(dp) :: dynmatint(2*nsym,2,3,natom,3,natom),gprimd(3,3) 2053 real(dp) :: symcart(3,3,nsym) 2054 2055 ! ********************************************************************* 2056 2057 ! 0) initializations 2058 call matr3inv(rprimd,gprimd) 2059 do isym=1,nsym 2060 call mati3inv(symrel(:,:,isym),symrec(:,:,isym)) 2061 end do 2062 2063 TqR=zero 2064 TqS_=zero 2065 dynmat=zero 2066 2067 !Note: dynmat is used as work space here 2068 i1=0 2069 do iat=1,natom 2070 do idir=1,3 2071 i1=i1+1 2072 i2=0 2073 do jat=1,natom 2074 do jdir=1,3 2075 i2=i2+1 2076 index=i1+3*natom*(i2-1) 2077 dynmat(1,idir,iat,jdir,jat)=dmati(2*index-1) 2078 dynmat(2,idir,iat,jdir,jat)=dmati(2*index ) 2079 end do 2080 end do 2081 end do 2082 end do 2083 2084 !Transform symrel to cartesian coordinates (RC coding) 2085 !do isym=1,nsym 2086 !symcart(:,:,isym)=matmul(rprimd,matmul(dble(symrel(:,:,isym)),gprimd)) 2087 !end do 2088 2089 !Coding from symdm9 2090 do isym=1,nsym 2091 do jj=1,3 2092 symcart(:,jj,isym)=zero 2093 do kk=1,3 2094 do ll=1,3 2095 symcart(:,jj,isym)=symcart(:,jj,isym)+rprimd(:,kk)*gprimd(jj,ll)*symrel(kk,ll,isym) 2096 end do 2097 end do 2098 end do 2099 end do 2100 2101 ! Get the symq of the CURRENT Q POINT 2102 ! mjv: set prtvol=0 for production runs. 2103 call littlegroup_q(nsym,qptn,symq,symrec,symafm,timrev,prtvol=0) 2104 2105 indij(:,:)=0 2106 dynmatint=zero 2107 2108 do isym=1,nsym ! loop over all the symmetries 2109 ! write(std_out,*) 'current symmetry',isym 2110 do itirev=1,2 ! loop over the time-reversal symmetry 2111 isgn=3-2*itirev 2112 ! write(std_out,*) 'timereversal',isgn 2113 2114 if (symq(4,itirev,isym)==1) then ! isym belongs to the wave vector point group 2115 ! write(std_out,*) 'isym belongs to the wave vector point group' 2116 do iat=1,natom 2117 do jat=1,natom 2118 niat=indsym(4,isym,iat) ! niat={R|t}iat 2119 njat=indsym(4,isym,jat) ! njat={R|t}jat 2120 indij(niat,njat)=indij(niat,njat)+1 2121 ! write(std_out,'(a,5i5)') 'current status:',iat,jat,niat,njat,indij(niat,njat) 2122 ! phase calculation, arg1 and arg2 because of two-atom derivative 2123 arg1=two_pi*( qptn(1)*indsym(1,isym,iat)+& 2124 qptn(2)*indsym(2,isym,iat)+& 2125 qptn(3)*indsym(3,isym,iat) ) 2126 arg2=two_pi*( qptn(1)*indsym(1,isym,jat)+& 2127 qptn(2)*indsym(2,isym,jat)+& 2128 qptn(3)*indsym(3,isym,jat) ) 2129 2130 re=cos(arg1)*cos(arg2)+sin(arg1)*sin(arg2) 2131 im=isgn*(cos(arg2)*sin(arg1)-cos(arg1)*sin(arg2)) 2132 2133 do idir=1,3 ! loop over displacements 2134 do jdir=1,3 ! loop over displacements 2135 ! we pick the (iat,jat) (3x3) block of the dyn.mat. 2136 sumr=zero 2137 sumi=zero 2138 do ii=1,3 2139 do jj=1,3 2140 sumr=sumr+symcart(idir,ii,isym)*dynmat(1,ii,niat,jj,njat)*symcart(jdir,jj,isym) 2141 sumi=sumi+symcart(idir,ii,isym)*dynmat(2,ii,niat,jj,njat)*symcart(jdir,jj,isym) 2142 end do 2143 end do 2144 sumi=isgn*sumi 2145 2146 dynmatint(nsym*(itirev-1)+isym,1,idir,iat,jdir,jat)=re*sumr-im*sumi 2147 dynmatint(nsym*(itirev-1)+isym,2,idir,iat,jdir,jat)=re*sumi+im*sumr 2148 end do 2149 end do 2150 end do 2151 end do ! end treatment of the (iat,jat) (3x3) block of dynmat 2152 end if ! symmetry check 2153 end do ! time-reversal 2154 end do ! symmetries 2155 2156 !4) make the average, get the final symmetric dynamical matrix 2157 do iat=1,natom 2158 do jat=1,natom 2159 do idir=1,3 2160 do jdir=1,3 2161 dmint=zero 2162 do isym=1,2*nsym 2163 dmint=dmint+dynmatint(isym,1,idir,iat,jdir,jat) 2164 end do 2165 dynmat(1,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat)) 2166 dmint=zero 2167 do isym=1,2*nsym 2168 dmint=dmint+dynmatint(isym,2,idir,iat,jdir,jat) 2169 end do 2170 dynmat(2,idir,iat,jdir,jat)=dmint/dble(indij(iat,jat)) 2171 end do 2172 end do 2173 end do 2174 end do 2175 2176 i1=0 2177 do iat=1,natom 2178 do idir=1,3 2179 i1=i1+1 2180 i2=0 2181 do jat=1,natom 2182 do jdir=1,3 2183 i2=i2+1 2184 index=i1+3*natom*(i2-1) 2185 dmati(2*index-1)=dynmat(1,idir,iat,jdir,jat) 2186 dmati(2*index )=dynmat(2,idir,iat,jdir,jat) 2187 end do 2188 end do 2189 end do 2190 end do 2191 2192 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
SOURCE
4745 subroutine sytens(indsym,mpert,natom,nsym,rfpert,symrec,symrel) 4746 4747 !Arguments ------------------------------- 4748 !scalars 4749 integer,intent(in) :: mpert,natom,nsym 4750 !arrays 4751 integer,intent(in) :: indsym(4,nsym,natom),symrec(3,3,nsym),symrel(3,3,nsym) 4752 integer,intent(inout) :: rfpert(3,mpert,3,mpert,3,mpert) 4753 4754 !Local variables ------------------------- 4755 !scalars 4756 integer :: flag,found,i1dir,i1dir_,i1pert,i1pert_,i2dir,i2dir_,i2pert,i2pert_ 4757 integer :: i3dir,i3dir_,i3pert,i3pert_,idisy1,idisy2,idisy3,ipesy1,ipesy2 4758 integer :: ipesy3,isym 4759 !arrays 4760 integer :: sym1(3,3),sym2(3,3),sym3(3,3) 4761 integer,allocatable :: pertsy(:,:,:,:,:,:) 4762 4763 !*********************************************************************** 4764 4765 ABI_MALLOC(pertsy,(3,mpert,3,mpert,3,mpert)) 4766 pertsy(:,:,:,:,:,:) = 0 4767 4768 !Loop over perturbations 4769 4770 do i1pert_ = 1, mpert 4771 do i2pert_ = 1, mpert 4772 do i3pert_ = 1, mpert 4773 4774 do i1dir_ = 1, 3 4775 do i2dir_ = 1, 3 4776 do i3dir_ = 1, 3 4777 4778 i1pert = (mpert - i1pert_ + 1) 4779 if (i1pert <= natom) i1pert = natom + 1 - i1pert 4780 i2pert = (mpert - i2pert_ + 1) 4781 if (i2pert <= natom) i2pert = natom + 1 - i2pert 4782 i3pert = (mpert - i3pert_ + 1) 4783 if (i3pert <= natom) i3pert = natom + 1 - i3pert 4784 4785 if (i1pert <= natom) then 4786 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 4787 else if (i2pert <= natom) then 4788 i1dir = i2dir_ ; i2dir = i1dir_ ; i3dir = i3dir_ 4789 else if (i3pert <= natom) then 4790 i1dir = i3dir_ ; i2dir = i2dir_ ; i3dir = i1dir_ 4791 else 4792 i1dir = i1dir_ ; i2dir = i2dir_ ; i3dir = i3dir_ 4793 end if 4794 4795 if (rfpert(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) /= 0) then 4796 4797 ! Loop over all symmetries 4798 4799 flag = 0 4800 do isym = 1, nsym 4801 4802 found = 1 4803 4804 ! Select the symmetric element of i1pert,i2pert,i3pert 4805 4806 if (i1pert <= natom) then 4807 ipesy1 = indsym(4,isym,i1pert) 4808 sym1(:,:) = symrec(:,:,isym) 4809 else if (i1pert == natom + 2) then 4810 ipesy1 = i1pert 4811 sym1(:,:) = symrel(:,:,isym) 4812 else 4813 found = 0 4814 end if 4815 4816 if (i2pert <= natom) then 4817 ipesy2 = indsym(4,isym,i2pert) 4818 sym2(:,:) = symrec(:,:,isym) 4819 else if (i2pert == natom + 2) then 4820 ipesy2 = i2pert 4821 sym2(:,:) = symrel(:,:,isym) 4822 else 4823 found = 0 4824 end if 4825 4826 if (i3pert <= natom) then 4827 ipesy3 = indsym(4,isym,i3pert) 4828 sym3(:,:) = symrec(:,:,isym) 4829 else if (i3pert == natom + 2) then 4830 ipesy3 = i3pert 4831 sym3(:,:) = symrel(:,:,isym) 4832 else 4833 found = 0 4834 end if 4835 4836 ! See if the symmetric element is available and check if some 4837 ! of the elements may be zeor. In the latter case, they do not need 4838 ! to be computed. 4839 4840 4841 if ((flag /= -1).and.& 4842 & (ipesy1==i1pert).and.(ipesy2==i2pert).and.(ipesy3==i3pert)) then 4843 flag = sym1(i1dir,i1dir)*sym2(i2dir,i2dir)*sym3(i3dir,i3dir) 4844 end if 4845 4846 4847 do idisy1 = 1, 3 4848 do idisy2 = 1, 3 4849 do idisy3 = 1, 3 4850 4851 if ((sym1(i1dir,idisy1) /= 0).and.(sym2(i2dir,idisy2) /= 0).and.& 4852 & (sym3(i3dir,idisy3) /= 0)) then 4853 if (pertsy(idisy1,ipesy1,idisy2,ipesy2,idisy3,ipesy3) == 0) then 4854 found = 0 4855 ! exit ! exit loop over symmetries 4856 end if 4857 end if 4858 4859 4860 if ((flag == -1).and.& 4861 & ((idisy1/=i1dir).or.(idisy2/=i2dir).or.(idisy3/=i3dir))) then 4862 if ((sym1(i1dir,idisy1)/=0).and.(sym2(i2dir,idisy2)/=0).and.& 4863 & (sym3(i3dir,idisy3)/=0)) then 4864 flag = 0 4865 end if 4866 end if 4867 4868 end do 4869 end do 4870 end do 4871 4872 if (found == 1) then 4873 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -1 4874 end if 4875 4876 ! In case a symmetry operation only changes the sign of an 4877 ! element, this element has to be equal to zero 4878 4879 if (flag == -1) then 4880 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = -2 4881 exit 4882 end if 4883 4884 end do ! close loop on symmetries 4885 4886 ! If the elemetn i1pert,i2pert,i3pert is not symmetric 4887 ! to a basis element, it is a basis element 4888 4889 if (pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) > -1) then 4890 pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) = 1 4891 end if 4892 4893 end if ! rfpert /= 0 4894 4895 end do ! close loop over perturbations 4896 end do 4897 end do 4898 end do 4899 end do 4900 end do 4901 4902 !Now, take into account the permutation of (i1pert,i1dir) 4903 !and (i3pert,i3dir) 4904 4905 do i1pert = 1, mpert 4906 do i2pert = 1, mpert 4907 do i3pert = 1, mpert 4908 4909 do i1dir = 1, 3 4910 do i2dir = 1, 3 4911 do i3dir = 1, 3 4912 4913 if ((i1pert /= i3pert).or.(i1dir /= i3dir)) then 4914 4915 if ((pertsy(i1dir,i1pert,i2dir,i2pert,i3dir,i3pert) == 1).and.& 4916 & (pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) == 1)) then 4917 pertsy(i3dir,i3pert,i2dir,i2pert,i1dir,i1pert) = -1 4918 end if 4919 4920 end if 4921 4922 end do 4923 end do 4924 end do 4925 4926 end do 4927 end do 4928 end do 4929 4930 rfpert(:,:,:,:,:,:) = pertsy(:,:,:,:,:,:) 4931 4932 ABI_FREE(pertsy) 4933 4934 end subroutine sytens
m_dynmat/wght9 [ Functions ]
[ Top ] [ m_dynmat ] [ Functions ]
NAME
wght9
FUNCTION
Generates a weight to each R point 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) toldist= Tolerance on the distance between two R points.
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
SOURCE
3874 subroutine wght9(brav,gprim,natom,ngqpt,nqpt,nqshft,nrpt,qshft,rcan,rpt,rprimd,toldist,r_inscribed_sphere,wghatm,ierr) 3875 3876 !Arguments ------------------------------- 3877 !scalars 3878 integer,intent(in) :: brav,natom,nqpt,nqshft,nrpt 3879 integer,intent(out) :: ierr 3880 real(dp),intent(out) :: r_inscribed_sphere 3881 real(dp),intent(in) :: toldist 3882 !arrays 3883 integer,intent(inout) :: ngqpt(9) 3884 real(dp),intent(in) :: gprim(3,3),qshft(3,4),rcan(3,natom),rpt(3,nrpt),rprimd(3,3) 3885 real(dp),intent(out) :: wghatm(natom,natom,nrpt) 3886 3887 !Local variables ------------------------- 3888 !scalars 3889 integer :: ia,ib,ii,jj,kk,iqshft,irpt,jqshft,nbordh,tok,nptws,nreq,idir 3890 real(dp) :: factor,sumwght,normsq,proj 3891 character(len=500) :: msg 3892 !arrays 3893 integer :: nbord(9) 3894 real(dp) :: rdiff(9),red(3,3),ptws(4, 729),pp(3),rdiff_tmp(3) 3895 3896 ! ********************************************************************* 3897 3898 ierr = 0 3899 3900 ! First analyze the vectors qshft 3901 if (nqshft /= 1) then 3902 3903 if (brav == 4) then 3904 write(msg,'(3a,i0,3a)' )& 3905 'For the time being, only nqshft=1',ch10,& 3906 'is allowed with brav=4, while it is nqshft=',nqshft,'.',ch10,& 3907 'Action: in the input file, correct either brav or nqshft.' 3908 ABI_ERROR(msg) 3909 end if 3910 3911 if (nqshft == 2) then 3912 ! Make sure that the q vectors form a BCC lattice 3913 do ii=1,3 3914 if(abs(abs(qshft(ii,1)-qshft(ii,2))-.5_dp)>1.d-10)then 3915 write(msg, '(a,a,a,a,a,a,a)' )& 3916 'The test of the q1shft vectors shows that they',ch10,& 3917 'do not generate a body-centered lattice, which',ch10,& 3918 'is mandatory for nqshft=2.',ch10,& 3919 'Action: change the q1shft vectors in your input file.' 3920 ABI_ERROR(msg) 3921 end if 3922 end do 3923 else if (nqshft == 4) then 3924 ! Make sure that the q vectors form a FCC lattice 3925 do iqshft=1,3 3926 do jqshft=iqshft+1,4 3927 tok=0 3928 do ii=1,3 3929 ! Test on the presence of a +-0.5 difference 3930 if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-.5_dp) <1.d-10) tok=tok+1 3931 ! Test on the presence of a 0 or +-1.0 difference 3932 if(abs(abs(qshft(ii,iqshft)-qshft(ii,jqshft))-1._dp) <1.d-10 .or.& 3933 abs(qshft(ii,iqshft)-qshft(ii,jqshft)) < 1.d-10) tok=tok+4 3934 end do 3935 ! Test 1 should be satisfied twice, and test 2 once 3936 if(tok/=6)then 3937 write(msg, '(7a)' )& 3938 'The test of the q1shft vectors shows that they',ch10,& 3939 'do not generate a face-centered lattice, which',ch10,& 3940 'is mandatory for nqshft=4.',ch10,& 3941 'Action: change the q1shft vectors in your input file.' 3942 ABI_ERROR(msg) 3943 end if 3944 end do 3945 end do 3946 else 3947 write(msg, '(a,i4,3a)' )& 3948 'nqshft must be 1, 2 or 4. It is nqshft=',nqshft,'.',ch10,& 3949 'Action: change nqshft in your input file.' 3950 ABI_ERROR(msg) 3951 end if 3952 end if 3953 3954 factor=0.5_dp 3955 if(brav==2 .or. brav==3) factor=0.25_dp 3956 if(nqshft/=1)factor=factor*2 3957 3958 if (brav==1) then 3959 ! Does not support multiple shifts 3960 if (nqshft/=1) then 3961 ABI_ERROR('This version of the weights does not support nqshft/=1.') 3962 end if 3963 3964 ! Find the points of the lattice given by ngqpt*acell. These are used to define 3965 ! a Wigner-Seitz cell around the origin. The origin is excluded from the list. 3966 ! TODO : in principle this should be only -1 to +1 for ii jj kk! 3967 nptws=0 3968 do ii=-2,2 3969 do jj=-2,2 3970 do kk=-2,2 3971 do idir=1,3 3972 pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3) 3973 end do 3974 normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) 3975 if (normsq > tol6) then 3976 nptws = nptws + 1 3977 ptws(:3,nptws) = pp(:) 3978 ptws(4,nptws) = half*normsq 3979 end if 3980 end do 3981 end do 3982 end do 3983 end if ! end new_wght 3984 !write(std_out,*)'factor,ngqpt',factor,ngqpt(1:3) 3985 3986 r_inscribed_sphere = sum((matmul(rprimd(:,:),ngqpt(1:3)))**2) 3987 do ii=-1,1 3988 do jj=-1,1 3989 do kk=-1,1 3990 if (ii==0 .and. jj==0 .and. kk==0) cycle 3991 do idir=1,3 3992 pp(idir)=ii*ngqpt(1)*rprimd(idir,1)+ jj*ngqpt(2)*rprimd(idir,2)+ kk*ngqpt(3)*rprimd(idir,3) 3993 end do 3994 normsq = pp(1)*pp(1)+pp(2)*pp(2)+pp(3)*pp(3) 3995 r_inscribed_sphere = min(r_inscribed_sphere, normsq) 3996 end do 3997 end do 3998 end do 3999 r_inscribed_sphere = sqrt(r_inscribed_sphere) 4000 4001 4002 !Begin the big loop on ia and ib 4003 do ia=1,natom 4004 do ib=1,natom 4005 4006 ! Simple Lattice 4007 if (abs(brav)==1) then 4008 ! In this case, it is better to work in reduced coordinates 4009 ! As rcan is in canonical coordinates, => multiplication by gprim 4010 do ii=1,3 4011 red(1,ii)= rcan(1,ia)*gprim(1,ii) +rcan(2,ia)*gprim(2,ii) +rcan(3,ia)*gprim(3,ii) 4012 red(2,ii)= rcan(1,ib)*gprim(1,ii) +rcan(2,ib)*gprim(2,ii) +rcan(3,ib)*gprim(3,ii) 4013 end do 4014 end if 4015 4016 do irpt=1,nrpt 4017 4018 ! Initialization of the weights to 1.0 4019 wghatm(ia,ib,irpt)=1.0_dp 4020 4021 ! Compute the difference vector 4022 4023 ! Simple Cubic Lattice 4024 if (abs(brav)==1) then 4025 ! Change of rpt to reduced coordinates 4026 do ii=1,3 4027 red(3,ii)= rpt(1,irpt)*gprim(1,ii) +rpt(2,irpt)*gprim(2,ii) +rpt(3,irpt)*gprim(3,ii) 4028 rdiff(ii)=red(2,ii)-red(1,ii)+red(3,ii) 4029 end do 4030 if (brav==1) then 4031 ! rdiff in cartesian coordinates 4032 do ii=1,3 4033 rdiff_tmp(ii)=rdiff(1)*rprimd(ii,1)+rdiff(2)*rprimd(ii,2)+rdiff(3)*rprimd(ii,3) 4034 end do 4035 rdiff(1:3)=rdiff_tmp(1:3) 4036 end if 4037 4038 else 4039 ! Other lattices 4040 do ii=1,3 4041 rdiff(ii)=rcan(ii,ib)-rcan(ii,ia)+rpt(ii,irpt) 4042 end do 4043 end if 4044 4045 ! Assignement of weights 4046 4047 if(nqshft==1 .and. brav/=4)then 4048 4049 if (brav/=1) then 4050 do ii=1,3 4051 ! If the rpt vector is greater than the allowed space => weight = 0.0 4052 if (abs(rdiff(ii))-tol10>factor*ngqpt(ii)) then 4053 wghatm(ia,ib,irpt)=zero 4054 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4055 ! If the point is in a boundary position => weight/2 4056 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4057 end if 4058 end do 4059 else 4060 ! new weights 4061 wghatm(ia,ib,irpt)=zero 4062 nreq = 1 4063 do ii=1,nptws 4064 proj = rdiff(1)*ptws(1,ii)+rdiff(2)*ptws(2,ii)+rdiff(3)*ptws(3,ii) 4065 ! if rdiff closer to ptws than the origin the weight is zero 4066 ! if rdiff close to the origin with respect to all the other ptws the weight is 1 4067 ! if rdiff is equidistant from the origin and N other ptws the weight is 1/(N+1) 4068 if (proj - ptws(4,ii) > toldist) then 4069 nreq = 0 4070 EXIT 4071 else if(abs(proj-ptws(4,ii)) <= toldist) then 4072 nreq=nreq+1 4073 end if 4074 end do 4075 if (nreq>0) then 4076 wghatm(ia,ib,irpt)=one/DBLE(nreq) 4077 end if 4078 end if 4079 4080 else if(brav==4)then 4081 ! Hexagonal 4082 ! Examination of the X and Y boundaries in order to form an hexagon 4083 ! First generate the relevant boundaries 4084 rdiff(4)=0.5_dp*( rdiff(1)+sqrt(3.0_dp)*rdiff(2) ) 4085 ngqpt(4)=ngqpt(1) 4086 rdiff(5)=0.5_dp*( rdiff(1)-sqrt(3.0_dp)*rdiff(2) ) 4087 ngqpt(5)=ngqpt(1) 4088 4089 ! Test the four inequalities 4090 do ii=1,5 4091 if(ii/=2)then 4092 4093 nbord(ii)=0 4094 ! If the rpt vector is greater than the allowed space => weight = 0.0 4095 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4096 wghatm(ia,ib,irpt)=zero 4097 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4098 ! If the point is in a boundary position increment nbord(ii) 4099 nbord(ii)=1 4100 end if 4101 4102 end if 4103 end do 4104 4105 ! Computation of weights 4106 nbordh=nbord(1)+nbord(4)+nbord(5) 4107 if (nbordh==1) then 4108 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4109 else if (nbordh==2) then 4110 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4111 else if (nbordh/=0) then 4112 ABI_BUG('There is a problem of borders and weights (hex).') 4113 end if 4114 if (nbord(3)==1)then 4115 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4116 end if 4117 4118 else if(nqshft==2 .and. brav/=4)then 4119 4120 ! BCC packing of k-points 4121 ! First, generate the relevant boundaries 4122 rdiff(4)= rdiff(1)+rdiff(2) 4123 rdiff(5)= rdiff(1)-rdiff(2) 4124 rdiff(6)= rdiff(1)+rdiff(3) 4125 rdiff(7)= rdiff(1)-rdiff(3) 4126 rdiff(8)= rdiff(3)+rdiff(2) 4127 rdiff(9)= rdiff(3)-rdiff(2) 4128 if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then 4129 write(msg, '(a,a,a,3i6,a,a,a,a)' )& 4130 'In the BCC case, the three ngqpt numbers ',ch10,& 4131 ' ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,& 4132 'should be equal.',ch10,& 4133 'Action: use identical ngqpt(1:3) in your input file.' 4134 ABI_ERROR(msg) 4135 end if 4136 do ii=4,9 4137 ngqpt(ii)=ngqpt(1) 4138 end do 4139 4140 ! Test the relevant inequalities 4141 nbord(1)=0 4142 do ii=4,9 4143 ! If the rpt vector is greater than the allowed space => weight = 0.0 4144 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4145 wghatm(ia,ib,irpt)=zero 4146 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4147 ! If the point is in a boundary position increment nbord(1) 4148 nbord(1)=nbord(1)+1 4149 end if 4150 end do 4151 4152 ! Computation of weights 4153 if (nbord(1)==1) then 4154 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4155 else if (nbord(1)==2) then 4156 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4157 else if (nbord(1)==3) then 4158 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4 4159 else if (nbord(1)==4) then 4160 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/6 4161 else if (nbord(1)/=0) then 4162 ABI_ERROR(' There is a problem of borders and weights (BCC).') 4163 end if 4164 4165 else if(nqshft==4 .and. brav/=4)then 4166 4167 ! FCC packing of k-points 4168 ! First, generate the relevant boundaries 4169 rdiff(4)= (rdiff(1)+rdiff(2)+rdiff(3))*2._dp/3._dp 4170 rdiff(5)= (rdiff(1)-rdiff(2)+rdiff(3))*2._dp/3._dp 4171 rdiff(6)= (rdiff(1)+rdiff(2)-rdiff(3))*2._dp/3._dp 4172 rdiff(7)= (rdiff(1)-rdiff(2)-rdiff(3))*2._dp/3._dp 4173 if(ngqpt(2)/=ngqpt(1) .or. ngqpt(3)/=ngqpt(1))then 4174 write(msg, '(a,a,a,3i6,a,a,a,a)' )& 4175 'In the FCC case, the three ngqpt numbers ',ch10,& 4176 ' ',ngqpt(1),ngqpt(2),ngqpt(3),ch10,& 4177 'should be equal.',ch10,& 4178 'Action: use identical ngqpt(1:3) in your input file.' 4179 ABI_ERROR(msg) 4180 end if 4181 do ii=4,7 4182 ngqpt(ii)=ngqpt(1) 4183 end do 4184 4185 ! Test the relevant inequalities 4186 nbord(1)=0 4187 do ii=1,7 4188 ! If the rpt vector is greater than the allowed space => weight = 0.0 4189 if (abs(rdiff(ii))-1.0d-10>factor*ngqpt(ii)) then 4190 wghatm(ia,ib,irpt)=zero 4191 ! If the point is in a boundary position increment nbord(1) 4192 else if (abs(abs(rdiff(ii))-factor*ngqpt(ii)) <=1.0d-10) then 4193 nbord(1)=nbord(1)+1 4194 end if 4195 end do 4196 4197 ! Computation of weights 4198 if (nbord(1)==1) then 4199 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/2 4200 else if (nbord(1)==2) then 4201 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/3 4202 else if (nbord(1)==3) then 4203 wghatm(ia,ib,irpt)=wghatm(ia,ib,irpt)/4 4204 else if (nbord(1)/=0 .and. wghatm(ia,ib,irpt)>1.d-10) then 4205 ! Interestingly nbord(1)==4 happens for some points outside of the volume 4206 ABI_BUG(' There is a problem of borders and weights (FCC).') 4207 end if 4208 4209 else 4210 write(msg, '(3a,i0,a)' )& 4211 'One should not arrive here ... ',ch10,& 4212 'The value nqshft ',nqshft,' is not available' 4213 ABI_BUG(msg) 4214 end if 4215 end do ! Assignement of weights is done 4216 end do ! End of the double loop on ia and ib 4217 end do 4218 4219 ! Check the results 4220 do ia=1,natom 4221 do ib=1,natom 4222 sumwght=zero 4223 do irpt=1,nrpt 4224 ! Check if the sum of the weights is equal to the number of q points 4225 sumwght=sumwght+wghatm(ia,ib,irpt) 4226 !write(std_out,'(a,3(i0,1x))' )' atom1, atom2, irpt ; rpt ; wghatm ',ia,ib,irpt 4227 !write(std_out,'(3es16.6,es18.6)' )rpt(1,irpt),rpt(2,irpt),rpt(3,irpt),wghatm(ia,ib,irpt) 4228 end do 4229 if (abs(sumwght-nqpt)>tol10) ierr = 1 4230 end do 4231 end do 4232 4233 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
SOURCE
2506 subroutine wings3(carflg,d2cart,mpert) 2507 2508 !Arguments ------------------------------- 2509 !scalars 2510 integer,intent(in) :: mpert 2511 !arrays 2512 integer,intent(inout) :: carflg(3,mpert,3,mpert) 2513 real(dp),intent(inout) :: d2cart(2,3,mpert,3,mpert) 2514 2515 !Local variables ------------------------- 2516 !scalars 2517 integer :: idir,idir1,ipert,ipert1 2518 2519 ! ********************************************************************* 2520 2521 do ipert=1,mpert 2522 do idir=1,3 2523 if(carflg(idir,ipert,idir,ipert)==0)then 2524 do ipert1=1,mpert 2525 do idir1=1,3 2526 carflg(idir,ipert,idir1,ipert1)=0 2527 carflg(idir1,ipert1,idir,ipert)=0 2528 d2cart(1,idir,ipert,idir1,ipert1)=zero 2529 d2cart(2,idir,ipert,idir1,ipert1)=zero 2530 d2cart(1,idir1,ipert1,idir,ipert)=zero 2531 d2cart(2,idir1,ipert1,idir,ipert)=zero 2532 end do 2533 end do 2534 end if 2535 end do 2536 end do 2537 2538 end subroutine wings3