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