TABLE OF CONTENTS
ABINIT/alignph [ Functions ]
NAME
alignph
FUNCTION
Construct linear combinations of the phonon eigendisplacements of degenerate modes in order to align the mode effective charges along the axes of the cartesian frame.
COPYRIGHT
Copyright (C) 1999-2018 ABINIT group (MVeithen) This file is distributed under the terms of the GNU General Public Licence, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
INPUTS
amu(ntypat)=mass of the atoms (atomic mass unit) displ(2,3*natom,3*natom)= 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. d2cart(2,3,mpert,3,mpert)= dynamical matrix, effective charges, dielectric tensor,.... all in cartesian coordinates mpert =maximum number of ipert natom=number of atoms in unit cell ntypat=number of types of atoms 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. typat(natom)=integer label of each type of atom (1,2,...)
OUTPUT
displ(2,3*natom,3*natom)= the displacements of atoms in cartesian coordinates. The eigendisplacements of degenerate modes have been aligned along the cartesian axes.
PARENTS
ddb_diel
CHILDREN
SOURCE
51 #if defined HAVE_CONFIG_H 52 #include "config.h" 53 #endif 54 55 #include "abi_common.h" 56 57 58 subroutine alignph(amu,displ,d2cart,mpert,natom,ntypat,phfrq,typat) 59 60 use defs_basis 61 use m_profiling_abi 62 63 !This section has been created automatically by the script Abilint (TD). 64 !Do not modify the following lines by hand. 65 #undef ABI_FUNC 66 #define ABI_FUNC 'alignph' 67 !End of the abilint section 68 69 implicit none 70 71 !Arguments ------------------------------- 72 !scalars 73 integer,intent(in) :: mpert,natom,ntypat 74 !arrays 75 integer,intent(in) :: typat(natom) 76 real(dp),intent(in) :: amu(ntypat),d2cart(2,3,mpert,3,mpert),phfrq(3*natom) 77 real(dp),intent(inout) :: displ(2,3*natom,3*natom) 78 79 !Local variables ------------------------- 80 !scalars 81 integer :: i1,idir1,idir2,ii,imode,imodex,imodey,imodez,ipert1 82 real(dp) :: theta 83 !arrays 84 integer,allocatable :: deg(:) 85 real(dp) :: zvec(3,3),zvect(3,3) 86 real(dp),allocatable :: modez(:,:,:),modezabs(:),oscstr(:,:,:),vec(:,:),vect(:,:) 87 88 ! ********************************************************************* 89 90 !DEBUG 91 ! write(std_out,*)'alignph : enter' 92 !ENDDEBUG 93 94 !Get the oscillator strength and mode effective charge for each mode 95 ABI_ALLOCATE(oscstr,(2,3,3*natom)) 96 ABI_ALLOCATE(modez,(2,3,3*natom)) 97 ABI_ALLOCATE(modezabs,(3*natom)) 98 ABI_ALLOCATE(vec,(3*natom,3)) 99 ABI_ALLOCATE(vect,(3*natom,3)) 100 ABI_ALLOCATE(deg,(3*natom)) 101 102 write(std_out,'(a,a)')ch10,' alignph : before modifying the eigenvectors, mode number and mode effective charges :' 103 do imode=1,3*natom 104 modezabs(imode)=zero 105 do ii=1,2 106 do idir2=1,3 107 oscstr(ii,idir2,imode)=zero 108 modez(ii,idir2,imode)=zero 109 do idir1=1,3 110 do ipert1=1,natom 111 i1=idir1+(ipert1-1)*3 112 oscstr(ii,idir2,imode)=oscstr(ii,idir2,imode)+& 113 & displ(ii,i1,imode)*& 114 & d2cart(1,idir1,ipert1,idir2,natom+2) 115 modez(ii,idir2,imode)=modez(ii,idir2,imode)+& 116 & displ(ii,i1,imode)*& 117 & d2cart(1,idir1,ipert1,idir2,natom+2)*& 118 & sqrt(amu(typat(ipert1))*amu_emass) 119 end do 120 end do 121 if(abs(modez(ii,idir2,imode))>modezabs(imode))modezabs(imode)=abs(modez(ii,idir2,imode)) 122 end do 123 end do 124 write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode) 125 end do 126 127 !Find degenerate modes with non-zero mode effective charge 128 imode = 0 129 do while (imode < 3*natom) 130 imode = imode + 1 131 if (imode == 3*natom) then 132 deg(imode) = 1 133 else if (abs(phfrq(imode) - phfrq(imode+1)) > tol6 .or. modezabs(imode)<tol8 .or. modezabs(imode+1)<tol8) then 134 ! Differ by phonon frequency or zero mode effective charge 135 deg(imode) = 1 136 else 137 deg(imode) = 2 138 if (imode < 3*natom - 1) then 139 if (abs(phfrq(imode) - phfrq(imode+2)) < tol6 .and. modezabs(imode+2)>tol8) then 140 deg(imode) = 3 141 imode = imode + 1 142 end if 143 end if 144 imode = imode + 1 145 end if 146 end do 147 148 149 !In case of a degenerate mode, with non-zero mode effective charge, align the mode effective charge vector along 150 !the axes of the cartesian frame 151 imode = 1 152 do while (imode <= 3*natom) 153 154 write(std_out,'(a,a,i4,a,i2)')ch10,' Mode number ',imode,' has degeneracy ',deg(imode) 155 write(std_out,'(a,3es16.6)') ' Mode effective charge of this mode =',modez(1,:,imode) 156 157 if (deg(imode) == 2) then 158 159 ! Optimize on the x direction 160 write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1) 161 if (abs(modez(1,1,imode)) > tol8) then 162 theta = atan(-modez(1,1,imode+1)/modez(1,1,imode)) 163 vec(:,1) = displ(1,:,imode) 164 vec(:,2) = displ(1,:,imode+1) 165 displ(1,:,imode) = cos(theta)*vec(:,1) - sin(theta)*vec(:,2) 166 displ(1,:,imode+1) = sin(theta)*vec(:,1) + cos(theta)*vec(:,2) 167 end if 168 169 else if (deg(imode) == 3) then 170 171 write(std_out,'(a,3es16.6)') ' Mode effective charge of next mode =',modez(1,:,imode+1) 172 write(std_out,'(a,3es16.6)') ' Mode effective charge of next-next mode =',modez(1,:,imode+2) 173 174 ! Before mixing them, select the mode-effective charge vectors as being predominently "x", "y" or "z" type. 175 if(abs(modez(1,1,imode))>abs(modez(1,2,imode))-tol12 .and. & 176 & abs(modez(1,1,imode))>abs(modez(1,3,imode))-tol12) then 177 imodex=imode 178 if(abs(modez(1,2,imode+1))>abs(modez(1,3,imode+1))-tol12)then 179 imodey=imode+1 ; imodez=imode+2 180 else 181 imodez=imode+1 ; imodey=imode+2 182 end if 183 else if(abs(modez(1,2,imode))>abs(modez(1,1,imode))-tol12 .and. & 184 & abs(modez(1,2,imode))>abs(modez(1,3,imode))-tol12) then 185 imodey=imode 186 if(abs(modez(1,1,imode+1))>abs(modez(1,3,imode+1))-tol12)then 187 imodex=imode+1 ; imodez=imode+2 188 else 189 imodez=imode+1 ; imodex=imode+2 190 end if 191 else 192 imodez=imode 193 if(abs(modez(1,1,imode+1))>abs(modez(1,2,imode+1))-tol12)then 194 imodex=imode+1 ; imodey=imode+2 195 else 196 imodey=imode+1 ; imodex=imode+2 197 end if 198 end if 199 vec(:,1)=displ(1,:,imodex) 200 vec(:,2)=displ(1,:,imodey) 201 vec(:,3)=displ(1,:,imodez) 202 zvec(:,1)=modez(1,:,imodex) 203 zvec(:,2)=modez(1,:,imodey) 204 zvec(:,3)=modez(1,:,imodez) 205 206 207 ! Optimize along x : does the first vector has a component along x ? 208 if (abs(zvec(1,1)) > tol8) then 209 ! Optimize on the (1,2) pair of modes along x 210 theta = atan(-zvec(1,2)/zvec(1,1)) 211 zvect(:,:)=zvec(:,:) 212 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,2) 213 zvec(:,2) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,2) 214 vect(:,:)=vec(:,:) 215 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,2) 216 vec(:,2) = sin(theta)*vect(:,1) + cos(theta)*vect(:,2) 217 ! Optimize on the (1,3) pair of modes along x 218 theta = atan(-zvec(1,3)/zvec(1,1)) 219 zvect(:,:)=zvec(:,:) 220 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3) 221 zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3) 222 vect(:,:)=vec(:,:) 223 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3) 224 vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3) 225 if (abs(zvec(2,2)) > tol8) then 226 ! Optimize on the (2,3) pair of modes along y 227 theta = atan(-zvec(2,3)/zvec(2,2)) 228 zvect(:,:)=zvec(:,:) 229 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 230 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 231 vect(:,:)=vec(:,:) 232 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 233 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 234 end if 235 ! Likely, the remaining is not needed ... because the vectors have been ordered in x, y, and z major component ... 236 ! Optimize along x : does the second vector has a component along x ? 237 else if(abs(zvec(1,2)) > tol8) then 238 ! Optimize on the (2,3) pair of modes along x 239 theta = atan(-zvec(1,3)/zvec(1,2)) 240 zvect(:,:)=zvec(:,:) 241 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 242 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 243 vect(:,:)=vec(:,:) 244 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 245 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 246 ! Optimize on the (1,3) pair of modes along y 247 if (abs(zvec(2,1)) > tol8) then 248 theta = atan(-zvec(2,3)/zvec(2,1)) 249 zvect(:,:)=zvec(:,:) 250 zvec(:,1) = cos(theta)*zvect(:,1) - sin(theta)*zvect(:,3) 251 zvec(:,3) = sin(theta)*zvect(:,1) + cos(theta)*zvect(:,3) 252 vect(:,:)=vec(:,:) 253 vec(:,1) = cos(theta)*vect(:,1) - sin(theta)*vect(:,3) 254 vec(:,3) = sin(theta)*vect(:,1) + cos(theta)*vect(:,3) 255 end if 256 ! We are left with the pair of vectors (2,3) 257 else if (abs(zvec(2,2)) > tol8) then 258 ! Optimize on the (2,3) pair of modes along y 259 theta = atan(-zvec(2,3)/zvec(2,2)) 260 zvect(:,:)=zvec(:,:) 261 zvec(:,2) = cos(theta)*zvect(:,2) - sin(theta)*zvect(:,3) 262 zvec(:,3) = sin(theta)*zvect(:,2) + cos(theta)*zvect(:,3) 263 vect(:,:)=vec(:,:) 264 vec(:,2) = cos(theta)*vect(:,2) - sin(theta)*vect(:,3) 265 vec(:,3) = sin(theta)*vect(:,2) + cos(theta)*vect(:,3) 266 end if 267 268 displ(1,:,imodex)=vec(:,1) 269 displ(1,:,imodey)=vec(:,2) 270 displ(1,:,imodez)=vec(:,3) 271 272 ! Previous coding, from Marek. Apparently, break the orthogonalization of vectors ... 273 ! do ii = 1,3 274 ! coeff(:) = 0._dp 275 ! if (ii == 1) then 276 ! jj = 2 ; kk = 3 277 ! else if (ii == 2) then 278 ! jj = 1 ; kk = 3 279 ! else 280 ! jj = 1 ; kk = 2 281 ! end if 282 ! coeff(ii) = 1._dp 283 ! c1 = modez(1,jj,imode+ii-1) 284 ! c2 = modez(1,jj,imode+jj-1) 285 ! c3 = modez(1,jj,imode+kk-1) 286 ! c4 = modez(1,kk,imode+ii-1) 287 ! c5 = modez(1,kk,imode+jj-1) 288 ! c6 = modez(1,kk,imode+kk-1) 289 ! dtm = c2*c6 - c3*c5 290 ! if (abs(dtm) > tol8) then 291 ! coeff(jj) = (c3*c4 - c1*c6)/dtm 292 ! coeff(kk) = (c1*c5 - c2*c4)/dtm 293 ! end if 294 ! mod_ = sqrt(1._dp + coeff(jj)*coeff(jj) + coeff(kk)*coeff(kk)) 295 ! coeff(:) = coeff(:)/mod_ 296 ! displ(1,:,imode+ii-1) = coeff(1)*vec(1,:) + coeff(2)*vec(2,:) + & 297 !& coeff(3)*vec(3,:) 298 ! end do 299 300 end if ! if deg mode 301 302 imode = imode + deg(imode) 303 304 end do 305 306 write(std_out,'(a,a)')ch10,' alignph : after modifying the eigenvectors, mode number and mode effective charges :' 307 do imode=1,3*natom 308 do ii=1,2 309 do idir2=1,3 310 modez(ii,idir2,imode)=zero 311 do idir1=1,3 312 do ipert1=1,natom 313 i1=idir1+(ipert1-1)*3 314 modez(ii,idir2,imode)=modez(ii,idir2,imode)+& 315 & displ(ii,i1,imode)*& 316 & d2cart(1,idir1,ipert1,idir2,natom+2)*& 317 & sqrt(amu(typat(ipert1))*amu_emass) 318 end do 319 end do 320 end do 321 end do 322 write(std_out,'(i4,3f16.6)')imode,modez(1,:,imode) 323 end do 324 325 ABI_DEALLOCATE(deg) 326 ABI_DEALLOCATE(oscstr) 327 ABI_DEALLOCATE(modez) 328 ABI_DEALLOCATE(modezabs) 329 ABI_DEALLOCATE(vec) 330 ABI_DEALLOCATE(vect) 331 332 end subroutine alignph