1! Copyright (C) 2016-2020 Marios Zacharias, Feliciano Giustino 2! 3! This file is distributed under the terms of the GNU General Public 4! License. See the file `LICENSE' in the root directory of the 5! present distribution, or http://www.gnu.org/copyleft.gpl.txt . 6! 7Module ifconstants 8 ! This code generates ZG displacements 9 ! 10 ! All variables read from file that need dynamical allocation 11 ! 12 USE kinds, ONLY: DP 13 REAL(DP), ALLOCATABLE :: frc(:,:,:,:,:,:,:), tau_blk(:,:), zeu(:,:,:), & 14 m_loc(:,:) 15 ! frc : interatomic force constants in REAL space 16 ! tau_blk : atomic positions for the original cell 17 ! zeu : effective charges for the original cell 18 ! m_loc: the magnetic moments of each atom 19 INTEGER, ALLOCATABLE :: ityp_blk(:) 20 ! ityp_blk : atomic types for each atom of the original cell 21 ! 22 CHARACTER(LEN=3), ALLOCATABLE :: atm(:) 23end Module ifconstants 24! 25!------------------------------------------------------------------------- 26PROGRAM ZG 27 !----------------------------------------------------------------------- 28 !----------------------------------------------------------------------- 29 !! authors: Marios Zacharias, Feliciano Giustino 30 !! acknowledgement: Hyungjun Lee for help packaging this release 31 !! version: v0.1 32 !! license: GNU 33 ! 34 ! This program generates the ZG_displacement for a list of 35 ! q vectors that are comensurate to the supercell size to be employed for 36 ! the special displacement method. The first part of the code is obtained by 37 ! modifying the matdyn.f90 subroutine of QE. The program starts from the 38 ! interatomic force constants generated from the DFPT phonon code through 39 ! the companion program q2r. 40 ! 41 ! ZG_displacement generates a supercell of the original cell with the atoms 42 ! displaced via Eq. (2) of https://arxiv.org/ABS/1912.10929. 43 ! Data required for the ZG_displacement are read from the force constant 44 ! file "*.fc" and. 45 ! 46 ! Input cards for ZG.in: namelist &input (first six as for matdyn.f90) 47 ! flfrc file produced by q2r containing force constants (needed) 48 ! It is the same as in the input of q2r.x (+ the .xml extension 49 ! IF the dynamical matrices produced by ph.x were in xml 50 ! format). No default value: must be specified. 51 ! asr (character) indicates the type of Acoustic Sum Rule imposed 52 ! - 'no': no Acoustic Sum Rules imposed (default) 53 ! - 'simple': previous implementation of the asr used 54 ! (3 translational asr imposed by correction of 55 ! the diagonal elements of the force constants matrix) 56 ! - 'crystal': 3 translational asr imposed by optimized 57 ! correction of the force constants (projection). 58 ! - 'one-dim': 3 translational asr + 1 rotational asr 59 ! imposed by optimized correction of the force constants 60 ! (the rotation axis is the direction of periodicity; 61 ! it will work only IF this axis considered is one of 62 ! the cartesian axis). 63 ! - 'zero-dim': 3 translational asr + 3 rotational asr 64 ! imposed by optimized correction of the force constants 65 ! Note that in certain cases, not all the rotational asr 66 ! can be applied (e.g. IF there are only 2 atoms in a 67 ! molecule or IF all the atoms are aligned, etc.). 68 ! In these cases the supplementary asr are cancelled 69 ! during the orthonormalization procedure (see below). 70 ! using tetrahedra and a uniform q-point grid (see below) 71 ! NB: may not work properly in noncubic materials 72 ! IF .FALSE. calculate phonon bands from the list of q-points 73 ! supplied in input (default) 74 ! amass masses of atoms in the supercell (a.m.u.), one per atom type 75 ! (default: use masses read from file flfrc) 76 ! "q_in_band_form" and "q_in_cryst_coord" meaningful if "q_external" 77 ! (see below) is set to .true. 78 ! q_in_band_form IF .TRUE. the q points are given in band form: 79 ! Only the first and last point of one or more lines 80 ! are given. See below. (default: .FALSE.). 81 ! q_in_cryst_coord IF .TRUE. input q points are in crystalline 82 ! coordinates (default: .FALSE.) 83 ! loto_2d set to .true. to activate two-dimensional treatment of LO-TO 84 ! siplitting. 85 ! 86 ! IF (q_in_band_form) THEN 87 ! nq ! number of q points 88 ! (q(i, n), i = 1, 3), nptq nptq is the number of points between this point 89 ! and the next. These points are automatically 90 ! generated. the q points are given in Cartesian 91 ! coordinates, 2pi/a units (a=lattice parameters) 92 ! ELSE : 93 ! nq number of q-points 94 ! (q(i, n), i = 1, 3) nq q-points in cartesian coordinates, 2pi/a units 95 ! If q = 0, the direction qhat (q=>0) for the non-analytic part 96 ! is extracted from the sequence of q-points as follows: 97 ! qhat = q(n) - q(n- 1) or qhat = q(n) - q(n+ 1) 98 ! depending on which one is available and nonzero. 99 ! For low-symmetry crystals, specify twice q = 0 in the list 100 ! IF you want to have q = 0 results for two different directions 101 ! ---------------------------------------------------------------------------------- 102 ! 103 ! Input cards to control "ZG_configuration" subroutine: 104 ! 105 ! "ZG_conf" : Logical flag that enables the creation of the ZG-displacement. 106 ! (default .true.) 107 ! "T" : Real number indicating the temperature at which the calculations will be performed. 108 ! "T" essentially defines the amplitude of the normal coordinates. 109 ! (default 0.00) 110 ! "dimx","dimy","dimz" : Integers corresponding to the dimensionality of the supercell. 111 ! (default 0,0,0) 112 ! "atm_zg(1), etc.." : String describing the element of each atomic species 113 ! (default "Element") 114 ! "synch" : Logical flag that enables the synchronization of the modes. 115 ! (default .false.) 116 ! "nloops" : Integer for the number of loops the algorithm needs to 117 ! go through for finding the optimum configuration. The algorithm 118 ! generates a set of "+,-,+,-" signs and its possible permutations, 119 ! trying to minimize the error coming from the coupling of modes with 120 ! the same q-wavevector but at different branch. For a finite supercell 121 ! size the order of using the "+,-,+,-" set and its permutations is 122 ! important giving different results. Therefore the algorithm checks 123 ! the combination that brings the error lower than a threshold. 124 ! (default 15000) 125 ! "compute_error" : Logical flag: if set to .true. allows the code to find the optimal ZG configuration 126 ! by minimizing the error based on the "threshold" flag (see below). Set it 127 ! to .false. if speed up is required. Setting it to .false. is useful when 128 ! (i) large supercell sizes are considered for which the error is minimized by 129 ! the first set of signs, and (ii) only single phonon displacements are of interest (see below) 130 ! (default .true.) 131 ! "threshold" : Real number indicating the error at which the algorithm stops while it's 132 ! looking for possible combinations of signs. Once this limit is reached 133 ! the ZG-displacement is constructed. The threshold is usually chosen 134 ! to be less than 5% of the diagonal terms, i.e. those terms that contribute 135 ! to the calculation of temperature-dependent properties. 136 ! (default 0.05) 137 ! "incl_qA" : Logical flag, to decide whether to include phonon modes in set A or not. 138 ! (default .true.) 139 ! "single_phonon_displ": Logical flag that allows to displace the nuclei along single phonon modes. 140 ! Use output configurations to compute electron-phonon matrix elements with a direct 141 ! supercell calculation. Set the displacement to the zero point by "T = 0". 142 ! This finite displacement should carry precisely the effect of diagonal elements of [g(q)+g(-q)]. 143 ! Output files: "single_phonon-displacements.dat" and 144 ! "single_phonon-velocities.dat". 145 ! (default .false.) 146 ! "q_external" : Logical flag that allows the use of a q-point list specified by the user in the input file. 147 ! If .false. the q-point list is specified by the supercell dimensions dimx, dimy, and dimz. 148 ! If .false. any q-point list after the input flags is ignored. 149 ! If .true. the q-point list must be provided by the user (see "qlist_AB.txt"). 150 ! (default .false.) 151 ! "qlist_AB.txt" : This file contains the external q-list in crystal coordinates as in the "ZG_444.in" example, 152 ! after the input flags. It corresponds to the q-points commensurate to the supercell size. 153 ! Only one of the q-point time-reversal partners is kept for the construction of the 154 ! ZG-displacement. The calculations, for the moment, assume systems with time-reversal symmetry. 155 ! For the generation of the "qlist_AB.txt" set the q-gird in file 156 ! "example/silicon/input/qlist.in" and run "../../../src/create_qlist.x < qlist.in > qlist.out". 157 ! One can modify the "create_qlist.f90" to generate a different path for consecutive q-points. 158 ! Paste the output of "qlist_AB.txt" to "ZG.in" after namelist &input. Set the flag 159 ! q_external = .true. for the code to read the list. 160 ! 161 USE kinds, ONLY : DP 162 USE mp, ONLY : mp_bcast 163 USE mp_world, ONLY : world_comm 164 USE mp_global, ONLY : mp_startup, mp_global_end 165 USE environment, ONLY : environment_start, environment_end 166 USE io_global, ONLY : ionode, ionode_id, stdout 167 USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, & 168 read_ifc_param, read_ifc 169 USE cell_base, ONLY : at, bg, celldm 170 USE constants, ONLY : RY_TO_THZ, RY_TO_CMM1, amu_ry 171 USE symm_base, ONLY : set_sym 172 USE rap_point_group, ONLY : code_group 173 USE bz_form, ONLY : transform_label_coord 174 USE parser, ONLY : read_line 175 USE rigid, ONLY : dyndiag, nonanal, nonanal_ifc 176 177 USE ifconstants, ONLY : frc, atm, zeu, tau_blk, ityp_blk, m_loc 178 ! 179 IMPLICIT NONE 180 ! 181 ! 182 ! variables *_blk refer to the original cell, other variables 183 ! to the (super)cell (which may coincide with the original cell) 184 ! 185 INTEGER, PARAMETER:: ntypx= 10, nrwsx=200 186 REAL(DP), PARAMETER :: eps = 1.0d-6 187 INTEGER :: nr1, nr2, nr3, nsc, ibrav 188 CHARACTER(LEN=256) :: flfrc, filename 189 CHARACTER(LEN= 10) :: asr 190 LOGICAL :: has_zstar, q_in_cryst_coord 191 COMPLEX(DP), ALLOCATABLE :: dyn(:,:,:,:), dyn_blk(:,:,:,:), frc_ifc(:,:,:,:) 192 COMPLEX(DP), ALLOCATABLE :: z(:,:) 193 REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), wq(:) 194 INTEGER, ALLOCATABLE:: ityp(:), itau_blk(:) 195 REAL(DP) :: omega, alat, &! cell parameters and volume 196 at_blk(3, 3), bg_blk(3, 3), &! original cell 197 omega_blk, &! original cell volume 198 epsil(3, 3), &! dielectric tensor 199 amass(ntypx), &! atomic masses 200 amass_blk(ntypx), &! original atomic masses 201 atws(3, 3), &! lattice vector for WS initialization 202 rws(0:3, nrwsx) ! nearest neighbor list, rws(0,*) = norm^2 203 ! 204 INTEGER :: nat, nat_blk, ntyp, ntyp_blk, & 205 l1, l2, l3, &! supercell dimensions 206 nrws, &! number of nearest neighbor 207 code_group_old 208 209 INTEGER :: nspin_mag, nqs, ios 210 ! 211 LOGICAL :: xmlifc, lo_to_split, loto_2d 212 ! 213 REAL(DP) :: qhat(3), qh, E 214 REAL(DP) :: delta 215 REAL(DP), ALLOCATABLE :: xqaux(:,:) 216 INTEGER, ALLOCATABLE :: nqb(:) 217 INTEGER :: n, i, j, it, nq, nqx, na, nb, nqtot 218 LOGICAL, EXTERNAL :: has_xml 219 INTEGER, ALLOCATABLE :: num_rap_mode(:,:) 220 LOGICAL, ALLOCATABLE :: high_sym(:) 221 LOGICAL :: q_in_band_form 222 ! .... variables for band plotting based on similarity of eigenvalues 223 COMPLEX(DP), ALLOCATABLE :: f_of_q(:,:,:,:) 224 INTEGER :: location(1), isig 225 CHARACTER(LEN=6) :: int_to_char 226 INTEGER :: npk_label, nch 227 CHARACTER(LEN=3), ALLOCATABLE :: letter(:) 228 INTEGER, ALLOCATABLE :: label_list(:) 229 LOGICAL :: tend, terr 230 CHARACTER(LEN=256) :: input_line, buffer 231 CHARACTER(LEN= 10) :: point_label_type 232 CHARACTER(len=80) :: k_points = 'tpiba' 233 ! 234 COMPLEX(DP), ALLOCATABLE :: z_nq_zg(:,:,:) ! nomdes, nmodes, nq 235 REAL(DP), ALLOCATABLE :: q_nq_zg(:,:) ! 3, nq 236 LOGICAL :: ZG_conf, synch, incl_qA, q_external 237 LOGICAL :: compute_error, single_phonon_displ 238 INTEGER :: dimx, dimy, dimz, nloops 239 REAL(DP) :: error_thresh, T 240 CHARACTER(LEN=3) :: atm_zg(ntypx) 241 ! 242 ! 243 NAMELIST /input/ flfrc, amass, asr, at, ntyp, loto_2d, & 244 & q_in_band_form, q_in_cryst_coord, point_label_type, & 245! we add the inputs for generating the ZG-configuration 246 & ZG_conf, dimx, dimy, dimz, nloops, error_thresh, q_external, & 247 & compute_error, synch, atm_zg, T, incl_qA, single_phonon_displ 248! ZG_conf --> IF TRUE compute the ZG_configuration 249 ! 250 CALL mp_startup() 251 CALL environment_start('ZG') 252 ! 253 l1= 1 254 l2= 1 255 l3= 1 256 IF (ionode) CALL input_from_file ( ) 257 ! 258 ! ... all calculations are done by the first cpu 259 ! 260 ! set namelist default 261 ! 262 asr ='no' 263 flfrc=' ' 264 amass(:) = 0.d0 265 amass_blk(:) = 0.d0 266 at(:,:) = 0.d0 267 ntyp = 0 268 q_in_band_form=.FALSE. 269 q_in_cryst_coord = .FALSE. 270 point_label_type='SC' 271 loto_2d=.FALSE. 272 ! 273 ZG_conf = .TRUE. 274 compute_error = .TRUE. 275 synch = .FALSE. 276 q_external = .FALSE. 277 incl_qA = .TRUE. 278 single_phonon_displ = .FALSE. 279 T = 0 280 error_thresh = 5.0E-02 281 dimx = 0 282 dimy = 0 283 dimz = 0 284 nloops = 15000 285 atm_zg = "Element" 286 ! 287 ! 288 ! 289 IF (ionode) READ (5, input,IOSTAT=ios) 290 CALL mp_bcast(ios, ionode_id, world_comm) 291 CALL errore('ZG', 'reading input namelist', ABS(ios)) 292 CALL mp_bcast(asr, ionode_id, world_comm) 293 CALL mp_bcast(flfrc, ionode_id, world_comm) 294 CALL mp_bcast(amass, ionode_id, world_comm) 295 CALL mp_bcast(amass_blk, ionode_id, world_comm) 296 CALL mp_bcast(at, ionode_id, world_comm) 297 CALL mp_bcast(ntyp, ionode_id, world_comm) 298 CALL mp_bcast(q_in_band_form, ionode_id, world_comm) 299 CALL mp_bcast(q_in_cryst_coord, ionode_id, world_comm) 300 CALL mp_bcast(point_label_type, ionode_id, world_comm) 301 CALL mp_bcast(loto_2d,ionode_id, world_comm) 302 ! 303 CALL mp_bcast(ZG_conf, ionode_id, world_comm) 304 CALL mp_bcast(compute_error, ionode_id, world_comm) 305 CALL mp_bcast(synch, ionode_id, world_comm) 306 CALL mp_bcast(q_external, ionode_id, world_comm) 307 CALL mp_bcast(incl_qA, ionode_id, world_comm) 308 CALL mp_bcast(single_phonon_displ, ionode_id, world_comm) 309 CALL mp_bcast(T, ionode_id, world_comm) 310 CALL mp_bcast(error_thresh, ionode_id, world_comm) 311 CALL mp_bcast(dimx, ionode_id, world_comm) 312 CALL mp_bcast(dimy, ionode_id, world_comm) 313 CALL mp_bcast(dimz, ionode_id, world_comm) 314 CALL mp_bcast(nloops, ionode_id, world_comm) 315 CALL mp_bcast(atm_zg, ionode_id, world_comm) 316 ! 317 ! To check that use specify supercell dimensions 318 IF (ZG_conf) THEN 319 IF ((dimx < 1) .OR. (dimy < 1) .OR. (dimz < 1)) CALL errore('ZG', 'reading supercell size', dimx) 320 ENDIF 321 ! 322 ! 323 ! read force constants 324 ! 325 ntyp_blk = ntypx ! avoids fake out-of-bound error 326 xmlifc=has_xml(flfrc) 327 IF (xmlifc) THEN 328 CALL read_dyn_mat_param(flfrc, ntyp_blk, nat_blk) 329 ALLOCATE (m_loc(3, nat_blk)) 330 ALLOCATE (tau_blk(3, nat_blk)) 331 ALLOCATE (ityp_blk(nat_blk)) 332 ALLOCATE (atm(ntyp_blk)) 333 ALLOCATE (zeu(3, 3, nat_blk)) 334 CALL read_dyn_mat_header(ntyp_blk, nat_blk, ibrav, nspin_mag, & 335 celldm, at_blk, bg_blk, omega_blk, atm, amass_blk, & 336 tau_blk, ityp_blk, m_loc, nqs, has_zstar, epsil, zeu ) 337 alat= celldm(1) 338 call volume(alat, at_blk(1, 1), at_blk(1, 2), at_blk(1, 3), omega_blk) 339 CALL read_ifc_param(nr1, nr2, nr3) 340 ALLOCATE(frc(nr1, nr2, nr3, 3, 3, nat_blk, nat_blk)) 341 CALL read_ifc(nr1, nr2, nr3, nat_blk,frc) 342 ELSE 343 CALL readfc ( flfrc, nr1, nr2, nr3, epsil, nat_blk, & 344 ibrav, alat, at_blk, ntyp_blk, & 345 amass_blk, omega_blk, has_zstar) 346 ENDIF 347 ! 348 CALL recips ( at_blk(1, 1), at_blk(1, 2), at_blk(1, 3), & 349 bg_blk(1, 1), bg_blk(1, 2), bg_blk(1, 3) ) 350 ! 351 ! set up (super)cell 352 ! 353 IF (ntyp < 0) THEN 354 call errore ('ZG','wrong ntyp ', ABS(ntyp)) 355 ELSE IF (ntyp == 0) THEN 356 ntyp =ntyp_blk 357 ENDIF 358 ! 359 ! masses (for mass approximation) 360 ! 361 DO it= 1, ntyp 362 IF (amass(it) < 0.d0) THEN 363 CALL errore ('ZG','wrong mass in the namelist', it) 364 ELSE IF (amass(it) == 0.d0) THEN 365 IF (it.LE.ntyp_blk) THEN 366 WRITE (stdout,'(a, i3, a, a)') ' mass for atomic type ', it, & 367 & ' not given; uses mass from file ',flfrc 368 amass(it) = amass_blk(it) 369 ELSE 370 CALL errore ('ZG','missing mass in the namelist', it) 371 ENDIF 372 ENDIF 373 ENDDO 374 ! 375 ! lattice vectors 376 ! 377 IF (SUM(ABS(at(:,:))) == 0.d0) THEN 378 IF (l1.LE.0 .OR. l2.LE.0 .OR. l3.LE.0) CALL & 379 & errore ('ZG',' wrong l1,l2 or l3', 1) 380 at(:, 1) = at_blk(:, 1) * DBLE(l1) 381 at(:, 2) = at_blk(:, 2) * DBLE(l2) 382 at(:, 3) = at_blk(:, 3) * DBLE(l3) 383 ENDIF 384 ! 385 CALL check_at(at, bg_blk, alat, omega) 386 ! 387 ! the supercell contains "nsc" times the original unit cell 388 ! 389 nsc = NINT(omega / omega_blk) 390 IF (ABS(omega / omega_blk-nsc) > eps) & 391 CALL errore ('ZG', 'volume ratio not integer', 1) 392 ! 393 ! read/generate atomic positions of the (super)cell 394 ! 395 nat = nat_blk * nsc 396 ! 397 ALLOCATE ( tau (3, nat), ityp(nat), itau_blk(nat) ) 398 ! 399 ! read atomic positions from IFC file 400 CALL set_tau & 401 (nat, nat_blk, at, at_blk, tau, tau_blk, ityp, ityp_blk, itau_blk) 402 ! 403 ! 404 ! reciprocal lattice vectors 405 ! 406 CALL recips (at(1, 1), at(1, 2), at(1, 3), bg(1, 1), bg(1, 2), bg(1, 3)) 407 ! 408 ! build the WS cell corresponding to the force constant grid 409 ! 410 atws(:, 1) = at_blk(:, 1) * DBLE(nr1) 411 atws(:, 2) = at_blk(:, 2) * DBLE(nr2) 412 atws(:, 3) = at_blk(:, 3) * DBLE(nr3) 413 ! initialize WS r-vectors 414 CALL wsinit(rws, nrwsx, nrws, atws) 415 ! 416 ! end of (super)cell setup 417 ! 418 ! 419 ! read q-point list 420 ! 421 ! 422 IF (.NOT. q_external) THEN 423 CALL qpoint_gen1(dimx, dimy, dimz, nq) 424 ! nq = ctrAB 425 CALL mp_bcast(nq, ionode_id, world_comm) 426 ALLOCATE ( q(3, nq) ) 427 CALL qpoint_gen2(dimx, dimy, dimz, nq, q) 428 ! 429 CALL mp_bcast(q, ionode_id, world_comm) 430 ! 431 CALL cryst_to_cart(nq, q, bg, +1) ! convert them to Cartesian 432 ELSE 433 ! 434 IF (ionode) READ (5, *) nq 435 CALL mp_bcast(nq, ionode_id, world_comm) 436 ALLOCATE ( q(3, nq) ) 437 IF (.NOT.q_in_band_form) THEN 438 DO n = 1, nq 439 IF (ionode) READ (5, *) (q(i, n), i = 1, 3) 440 ! IF (ionode) READ (5,'(3F10.6)') q(:, n) 441 ENDDO 442 CALL mp_bcast(q, ionode_id, world_comm) 443 ! 444 IF (q_in_cryst_coord) CALL cryst_to_cart(nq, q, bg, +1) 445 ELSE 446 ALLOCATE( nqb(nq) ) 447 ALLOCATE( xqaux(3, nq) ) 448 ALLOCATE( letter(nq) ) 449 ALLOCATE( label_list(nq) ) 450 npk_label= 0 451 DO n = 1, nq 452 CALL read_line( input_line, end_of_file = tend, error = terr ) 453 IF (tend) CALL errore('ZG','Missing lines', 1) 454 IF (terr) CALL errore('ZG','Error reading q points', 1) 455 DO j = 1, 256 ! loop over all characters of input_line 456 IF ( (ICHAR(input_line(j:j)) < 58 .AND. & ! a digit 457 ICHAR(input_line(j:j)) > 47) & 458 .OR.ICHAR(input_line(j:j)) == 43 .OR. & ! the + sign 459 ICHAR(input_line(j:j)) == 45 .OR. & ! the - sign 460 ICHAR(input_line(j:j)) == 46 ) THEN ! a dot . 461 ! 462 ! This is a digit, therefore this line contains the coordinates of the 463 ! k point. We read it and EXIT from the loop on characters 464 ! 465 READ(input_line,*) xqaux(1, n), xqaux(2, n), xqaux(3, n), & 466 nqb(n) 467 EXIT 468 ELSEIF ((ICHAR(input_line(j:j)) < 123 .AND. & 469 ICHAR(input_line(j:j)) > 64)) THEN 470 ! 471 ! This is a letter, not a space character. We read the next three 472 ! characters and save them in the letter array, save also which k point 473 ! it is 474 ! 475 npk_label=npk_label+ 1 476 READ(input_line(j:),'(a3)') letter(npk_label) 477 label_list(npk_label) =n 478 ! 479 ! now we remove the letters from input_line and read the number of points 480 ! of the line. The next two line should account for the case in which 481 ! there is only one space between the letter and the number of points. 482 ! 483 nch=3 484 IF ( ICHAR(input_line(j+ 1:j+ 1)) ==32 .OR. & 485 ICHAR(input_line(j+2:j+2)) ==32 ) nch=2 486 buffer =input_line(j+nch:) 487 READ(buffer,*, err =20, iostat=ios) nqb(n) 488 20 IF (ios /= 0) CALL errore('ZG',& 489 'problem reading number of points', 1) 490 EXIT 491 ENDIF 492 ENDDO 493 ENDDO 494 IF (q_in_cryst_coord) k_points ='crystal' 495 IF ( npk_label > 0 ) & 496 CALL transform_label_coord(ibrav, celldm, xqaux, letter, & 497 label_list, npk_label, nq, k_points, point_label_type ) 498 499 DEALLOCATE(letter) 500 DEALLOCATE(label_list) 501 502 CALL mp_bcast(xqaux, ionode_id, world_comm) 503 CALL mp_bcast(nqb, ionode_id, world_comm) 504 IF (q_in_cryst_coord) CALL cryst_to_cart(nq,xqaux, bg,+ 1) 505 nqtot=SUM(nqb(1 : nq - 1)) + 1 506 DO i = 1, nq - 1 507 IF (nqb(i) == 0) nqtot=nqtot+ 1 508 ENDDO 509 DEALLOCATE(q) 510 ALLOCATE(q(3, nqtot)) 511 ALLOCATE(wq(nqtot)) 512 CALL generate_k_along_lines(nq, xqaux, nqb, q, wq, nqtot) 513 nq = nqtot 514 DEALLOCATE(xqaux) 515 DEALLOCATE(nqb) 516 ENDIF 517 ! 518 ENDIF ! q_external, q-list 519 ! 520 IF (asr /= 'no') THEN 521 CALL set_asr (asr, nr1, nr2, nr3, frc, zeu, & 522 nat_blk, ibrav, tau_blk) 523 ENDIF 524 ! 525 ALLOCATE ( dyn(3, 3, nat, nat), dyn_blk(3, 3, nat_blk, nat_blk) ) 526 ALLOCATE ( z(3 * nat, 3 * nat), w2(3*nat, nq), f_of_q(3, 3, nat, nat) ) 527 ! 528 IF (ionode .AND. ZG_conf) THEN 529 ALLOCATE ( z_nq_zg(3 * nat, 3 * nat, nq), q_nq_zg(3, nq)) 530 z_nq_zg(:, :, :) = (0.d0, 0.d0) 531 q_nq_zg(:, :) = 0.d0 532 ENDIF 533 ! 534 535 IF (xmlifc) CALL set_sym(nat, tau, ityp, nspin_mag, m_loc ) 536 537 ALLOCATE(num_rap_mode(3 * nat, nq)) 538 ALLOCATE(high_sym(nq)) 539 num_rap_mode=- 1 540 high_sym=.TRUE. 541 542 DO n = 1, nq 543 dyn(:,:,:,:) = (0.d0, 0.d0) 544 545 lo_to_split = .FALSE. 546 f_of_q(:,:,:,:) = (0.d0, 0.d0) 547 548 CALL setupmat (q(1,n), dyn, nat, at, bg, tau, itau_blk, nsc, alat, & 549 dyn_blk, nat_blk, at_blk, bg_blk, tau_blk, omega_blk, & 550 loto_2d, & 551 epsil, zeu, frc, nr1,nr2,nr3, has_zstar, rws, nrws, f_of_q) 552 553 IF (.not.loto_2d) THEN 554 qhat(1) = q(1, n) * at(1, 1) + q(2, n) * at(2, 1) + q(3, n) * at(3, 1) 555 qhat(2) = q(1, n) * at(1, 2) + q(2, n) * at(2, 2) + q(3, n) * at(3, 2) 556 qhat(3) = q(1, n) * at(1, 3) + q(2, n) * at(2, 3) + q(3, n) * at(3, 3) 557 IF ( ABS( qhat(1) - NINT (qhat(1) ) ) <= eps .AND. & 558 ABS( qhat(2) - NINT (qhat(2) ) ) <= eps .AND. & 559 ABS( qhat(3) - NINT (qhat(3) ) ) <= eps ) THEN 560 ! 561 ! q = 0 : we need the direction q => 0 for the non-analytic part 562 ! 563 IF ( n == 1 ) THEN 564 ! IF q is the first point in the list 565 IF ( nq > 1 ) THEN 566 ! one more point 567 qhat(:) = q(:, n) - q(:, n+ 1) 568 ELSE 569 ! no more points 570 qhat(:) = 0.d0 571 ENDIF 572 ELSE IF ( n > 1 ) THEN 573 ! IF q is not the first point in the list 574 IF ( q(1, n- 1) == 0.d0 .AND. & 575 q(2, n- 1) == 0.d0 .AND. & 576 q(3, n- 1) == 0.d0 .AND. n < nq ) THEN 577 ! IF the preceding q is also 0 : 578 qhat(:) = q(:, n) - q(:, n+ 1) 579 ELSE 580 ! IF the preceding q is npt 0 : 581 qhat(:) = q(:, n) - q(:, n- 1) 582 ENDIF 583 ENDIF 584 qh = SQRT(qhat(1)**2 +qhat(2)**2+qhat(3) **2) 585 ! WRITE(*,*) ' qh, has_zstar ',qh, has_zstar 586 IF (qh /= 0.d0) qhat(:) = qhat(:) / qh 587 IF (qh /= 0.d0 .AND. .NOT. has_zstar) THEN 588 CALL infomsg & 589 ('ZG','Z* not found in file '//TRIM(flfrc)// & 590 ', TO-LO splitting at q = 0 will be ABSent!') 591 ELSE 592 lo_to_split=.TRUE. 593 ENDIF 594 ! 595 CALL nonanal (nat, nat_blk, itau_blk, epsil, qhat, zeu, omega, dyn) 596 ! 597 ENDIF 598 ! 599 END IF 600 601 CALL dyndiag(nat, ntyp, amass, ityp, dyn, w2(1, n), z) 602 ! fill a 3D matrix with all eigenvectors 603 CALL mp_bcast(z, ionode_id, world_comm) 604 IF (ionode .AND. ZG_conf) THEN 605 z_nq_zg(:, :, n) = z(:, :) 606 q_nq_zg(:, n) = q(:, n) 607 ENDIF 608 ! 609 ! Cannot use the small group of \Gamma to analize the symmetry 610 ! of the mode IF there is an electric field. 611 ! 612 IF (xmlifc.AND..NOT.lo_to_split) THEN 613 WRITE(stdout,'(10x,"xq=", 3F8.4)') q(:, n) 614 CALL find_representations_mode_q(nat, ntyp, q(:, n), & 615 w2(:, n), z,tau, ityp, amass, num_rap_mode(:, n), nspin_mag) 616 IF (code_group == code_group_old.OR.high_sym(n- 1)) high_sym(n) =.FALSE. 617 code_group_old= code_group 618 ENDIF 619 ! 620 ! 621 ! 622 ENDDO !nq 623 ! 624 ! 625 ! 626 ! 627 ! If the force constants are in the xml format we WRITE also 628 ! the file with the representations of each mode 629 ! 630 ! 631 ! 632 ! 633 CALL mp_bcast(w2, ionode_id, world_comm) 634 IF ( ionode .AND. ZG_conf ) call ZG_configuration(nq, nat, ntyp, amass, & 635 ityp, q_nq_zg, w2, z_nq_zg, ios, & 636 dimx, dimy, dimz, nloops, error_thresh, synch, tau, alat, atm_zg, & 637 ntypx, at, q_in_cryst_coord, q_external, T, incl_qA, & 638 compute_error, single_phonon_displ) 639 ! 640 ! 641 DEALLOCATE (z, w2, dyn, dyn_blk) 642 ! 643 IF (ionode .AND. ZG_conf) DEALLOCATE (z_nq_zg, q_nq_zg) 644 ! 645 ! 646 ! for a2F 647 ! 648 DEALLOCATE(num_rap_mode) 649 DEALLOCATE(high_sym) 650 ! 651 652 CALL environment_end('ZG') 653 ! 654 CALL mp_global_end() 655 ! 656 STOP 657 ! 658END PROGRAM ZG 659! 660!----------------------------------------------------------------------- 661SUBROUTINE readfc ( flfrc, nr1, nr2, nr3, epsil, nat, & 662 ibrav, alat, at, ntyp, amass, omega, has_zstar ) 663 !----------------------------------------------------------------------- 664 ! 665 USE kinds, ONLY : DP 666 USE ifconstants,ONLY : tau => tau_blk, ityp => ityp_blk, frc, zeu 667 USE cell_base, ONLY : celldm 668 USE io_global, ONLY : ionode, ionode_id, stdout 669 USE mp, ONLY : mp_bcast 670 USE mp_world, ONLY : world_comm 671 USE constants, ONLY : amu_ry 672 ! 673 IMPLICIT NONE 674 ! I/O variable 675 CHARACTER(LEN=256) :: flfrc 676 INTEGER :: ibrav, nr1, nr2, nr3, nat, ntyp 677 REAL(DP) :: alat, at(3, 3), epsil(3, 3) 678 LOGICAL :: has_zstar 679 ! local variables 680 INTEGER :: i, j, na, nb, m1,m2,m3 681 INTEGER :: ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid 682 REAL(DP) :: amass(ntyp), amass_from_file, omega 683 INTEGER :: nt 684 CHARACTER(LEN=3) :: atm 685 ! 686 ! 687 IF (ionode) OPEN (unit= 1,file=flfrc, status ='old',form='formatted') 688 ! 689 ! read cell data 690 ! 691 IF (ionode)THEN 692 READ(1,*) ntyp, nat, ibrav,(celldm(i), i = 1,6) 693 IF (ibrav== 0) THEN 694 read(1,*) ((at(i, j), i = 1, 3), j =1, 3) 695 ENDIF 696 ENDIF 697 CALL mp_bcast(ntyp, ionode_id, world_comm) 698 CALL mp_bcast(nat, ionode_id, world_comm) 699 CALL mp_bcast(ibrav, ionode_id, world_comm) 700 CALL mp_bcast(celldm, ionode_id, world_comm) 701 IF (ibrav== 0) THEN 702 CALL mp_bcast(at, ionode_id, world_comm) 703 ENDIF 704 ! 705 CALL latgen(ibrav, celldm, at(1, 1), at(1, 2), at(1, 3), omega) 706 alat = celldm(1) 707 at = at / alat ! bring at in units of alat 708 CALL volume(alat, at(1, 1), at(1, 2), at(1, 3), omega) 709 ! 710 ! read atomic types, positions and masses 711 ! 712 DO nt = 1, ntyp 713 IF (ionode) READ(1,*) i, atm, amass_from_file 714 CALL mp_bcast(i, ionode_id, world_comm) 715 CALL mp_bcast(atm, ionode_id, world_comm) 716 CALL mp_bcast(amass_from_file, ionode_id, world_comm) 717 IF (i.NE.nt) CALL errore ('readfc','wrong data read', nt) 718 IF (amass(nt).EQ.0.d0) THEN 719 amass(nt) = amass_from_file/amu_ry 720 ELSE 721 WRITE(stdout,*) 'for atomic type', nt,' mass from file not used' 722 ENDIF 723 ENDDO 724 ! 725 ALLOCATE (tau(3, nat), ityp(nat), zeu(3, 3, nat)) 726 ! 727 DO na= 1, nat 728 IF (ionode) READ(1,*) i, ityp(na),(tau(j, na), j = 1, 3) 729 CALL mp_bcast(i, ionode_id, world_comm) 730 IF (i.NE.na) CALL errore ('readfc','wrong data read', na) 731 ENDDO 732 CALL mp_bcast(ityp, ionode_id, world_comm) 733 CALL mp_bcast(tau, ionode_id, world_comm) 734 ! 735 ! read macroscopic variable 736 ! 737 IF (ionode) READ (1,*) has_zstar 738 CALL mp_bcast(has_zstar, ionode_id, world_comm) 739 IF (has_zstar) THEN 740 IF (ionode) READ(1,*) ((epsil(i, j), j = 1, 3), i =1, 3) 741 CALL mp_bcast(epsil, ionode_id, world_comm) 742 IF (ionode) THEN 743 DO na= 1, nat 744 READ(1,*) 745 READ(1,*) ((zeu(i, j, na), j = 1, 3), i =1, 3) 746 ENDDO 747 ENDIF 748 CALL mp_bcast(zeu, ionode_id, world_comm) 749 ELSE 750 zeu (:,:,:) = 0.d0 751 epsil(:,:) = 0.d0 752 ENDIF 753 ! 754 IF (ionode) READ (1,*) nr1, nr2, nr3 755 CALL mp_bcast(nr1, ionode_id, world_comm) 756 CALL mp_bcast(nr2, ionode_id, world_comm) 757 CALL mp_bcast(nr3, ionode_id, world_comm) 758 ! 759 ! read REAL-space interatomic force constants 760 ! 761 ALLOCATE ( frc(nr1, nr2, nr3, 3, 3, nat, nat) ) 762 frc(:,:,:,:,:,:,:) = 0.d0 763 DO i = 1, 3 764 DO j = 1, 3 765 DO na= 1, nat 766 DO nb= 1, nat 767 IF (ionode) READ (1,*) ibid, jbid, nabid, nbbid 768 CALL mp_bcast(ibid, ionode_id, world_comm) 769 CALL mp_bcast(jbid, ionode_id, world_comm) 770 CALL mp_bcast(nabid, ionode_id, world_comm) 771 CALL mp_bcast(nbbid, ionode_id, world_comm) 772 IF(i .NE.ibid .OR. j .NE.jbid .OR. & 773 na.NE.nabid .OR. nb.NE.nbbid) & 774 CALL errore ('readfc','error in reading', 1) 775 IF (ionode) READ (1,*) (((m1bid, m2bid, m3bid, & 776 frc(m1,m2,m3, i, j, na, nb), & 777 m1= 1, nr1),m2 =1, nr2),m3=1, nr3) 778 779 CALL mp_bcast(frc(:,:,:, i, j, na, nb), ionode_id, world_comm) 780 ENDDO 781 ENDDO 782 ENDDO 783 ENDDO 784 ! 785 IF (ionode) CLOSE(unit= 1) 786 ! 787 RETURN 788END SUBROUTINE readfc 789! 790!----------------------------------------------------------------------- 791SUBROUTINE frc_blk(dyn,q,tau, nat, nr1, nr2, nr3,frc, at, bg, rws, nrws,f_of_q) 792 !----------------------------------------------------------------------- 793 ! calculates the dynamical matrix at q from the (short-range part of the) 794 ! force constants 795 ! 796 USE kinds, ONLY : DP 797 USE constants, ONLY : tpi 798 USE io_global, ONLY : stdout 799 ! 800 IMPLICIT NONE 801 INTEGER nr1, nr2, nr3, nat, n1, n2, n3, nr1_, nr2_, nr3_, & 802 ipol, jpol, na, nb, m1, m2, m3, NINT, i, j, nrws 803 COMPLEX(DP) dyn(3, 3, nat, nat), f_of_q(3, 3, nat, nat) 804 REAL(DP) frc(nr1, nr2, nr3, 3, 3, nat, nat), tau(3, nat), q(3), arg, & 805 at(3, 3), bg(3, 3), r(3), weight, r_ws(3), & 806 total_weight, rws(0:3, nrws), alat 807 REAL(DP), EXTERNAL :: wsweight 808 REAL(DP),SAVE,ALLOCATABLE :: wscache(:,:,:,:,:) 809 REAL(DP), ALLOCATABLE :: ttt(:,:,:,:,:), tttx(:,:) 810 LOGICAL,SAVE :: first=.TRUE. 811 ! 812 nr1_=2*nr1 813 nr2_=2*nr2 814 nr3_=2*nr3 815 FIRST_TIME : IF (first) THEN 816 first=.FALSE. 817 ALLOCATE( wscache(-nr3_:nr3_, -nr2_:nr2_, -nr1_:nr1_, nat, nat) ) 818 DO na= 1, nat 819 DO nb= 1, nat 820 total_weight= 0.0d0 821 ! 822 DO n1=-nr1_, nr1_ 823 DO n2 =-nr2_, nr2_ 824 DO n3=-nr3_, nr3_ 825 DO i = 1, 3 826 r(i) = n1*at(i, 1) +n2*at(i, 2) +n3*at(i, 3) 827 r_ws(i) = r(i) + tau(i, na) -tau(i, nb) 828 ENDDO 829 wscache(n3, n2, n1, nb, na) = wsweight(r_ws, rws, nrws) 830 ENDDO 831 ENDDO 832 ENDDO 833 ENDDO 834 ENDDO 835 ENDIF FIRST_TIME 836 ! 837 ALLOCATE(ttt(3, nat, nr1, nr2, nr3)) 838 ALLOCATE(tttx(3, nat*nr1*nr2*nr3)) 839 ttt(:,:,:,:,:) = 0.d0 840 841 DO na= 1, nat 842 DO nb= 1, nat 843 total_weight= 0.0d0 844 DO n1=-nr1_, nr1_ 845 DO n2 =-nr2_, nr2_ 846 DO n3=-nr3_, nr3_ 847 ! 848 ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY SAFE RANGE! 849 ! 850 DO i = 1, 3 851 r(i) = n1*at(i, 1) +n2*at(i, 2) +n3*at(i, 3) 852 ENDDO 853 854 weight = wscache(n3, n2, n1, nb, na) 855 IF (weight > 0.0d0) THEN 856 ! 857 ! FIND THE VECTOR CORRESPONDING TO R IN THE ORIGINAL CELL 858 ! 859 m1 = MOD(n1 + 1, nr1) 860 IF(m1.LE.0) m1=m1 +nr1 861 m2 = MOD(n2 + 1, nr2) 862 IF(m2.LE.0) m2 =m2 +nr2 863 m3 = MOD(n3 + 1, nr3) 864 IF(m3.LE.0) m3=m3 +nr3 865 ! WRITE(*,'(6i4)') n1, n2, n3,m1,m2,m3 866 ! 867 ! FOURIER TRANSFORM 868 ! 869 DO i = 1, 3 870 ttt(i, na,m1,m2,m3) =tau(i, na) +m1*at(i, 1) +m2*at(i, 2) +m3*at(i, 3) 871 ttt(i, nb,m1,m2,m3) =tau(i, nb) +m1*at(i, 1) +m2*at(i, 2) +m3*at(i, 3) 872 ENDDO 873 874 arg = tpi* (q(1) *r(1) + q(2) *r(2) + q(3) *r(3)) 875 DO ipol= 1, 3 876 DO jpol= 1, 3 877 dyn(ipol, jpol, na, nb) = & 878 dyn(ipol, jpol, na, nb) + & 879 (frc(m1,m2,m3, ipol, jpol, na, nb) +f_of_q(ipol, jpol, na, nb)) & 880 *CMPLX(COS(arg),-SIN(arg), kind= DP) *weight 881 ENDDO 882 ENDDO 883 ENDIF 884 total_weight=total_weight + weight 885 ENDDO 886 ENDDO 887 ENDDO 888 IF (ABS(total_weight-nr1*nr2*nr3).GT.1.0d-8) THEN 889 WRITE(stdout,*) total_weight 890 CALL errore ('frc_blk','wrong total_weight', 1) 891 ENDIF 892 ENDDO 893 ENDDO 894 ! 895 RETURN 896END SUBROUTINE frc_blk 897! 898!----------------------------------------------------------------------- 899SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, & 900 & dyn_blk,nat_blk,at_blk,bg_blk,tau_blk,omega_blk, & 901 & loto_2d, & 902 & epsil,zeu,frc,nr1,nr2,nr3,has_zstar,rws,nrws,f_of_q) 903 !----------------------------------------------------------------------- 904 ! compute the dynamical matrix (the analytic part only) 905 ! 906 USE kinds, ONLY : DP 907 USE constants, ONLY : tpi 908 USE cell_base, ONLY : celldm 909 USE rigid, ONLY : rgd_blk 910 ! 911 IMPLICIT NONE 912 ! 913 ! I/O variables 914 ! 915 INTEGER:: nr1, nr2, nr3, nat, nat_blk, nsc, nrws, itau_blk(nat) 916 REAL(DP) :: q(3), tau(3, nat), at(3, 3), bg(3, 3), alat, & 917 epsil(3, 3), zeu(3, 3, nat_blk), rws(0:3, nrws), & 918 frc(nr1, nr2, nr3, 3, 3, nat_blk, nat_blk) 919 REAL(DP) :: tau_blk(3, nat_blk), at_blk(3, 3), bg_blk(3, 3), omega_blk 920 COMPLEX(DP) dyn_blk(3, 3, nat_blk, nat_blk), f_of_q(3, 3, nat, nat) 921 COMPLEX(DP) :: dyn(3, 3, nat, nat) 922 LOGICAL :: has_zstar 923 ! 924 ! local variables 925 ! 926 REAL(DP) :: arg 927 COMPLEX(DP) :: cfac(nat) 928 INTEGER :: i, j, k, na, nb, na_blk, nb_blk, iq 929 REAL(DP) :: qp(3), qbid(3, nsc) ! automatic array 930 LOGICAL :: loto_2d 931 ! 932 ! 933 CALL q_gen(nsc,qbid, at_blk, bg_blk, at, bg) 934 ! 935 DO iq= 1, nsc 936 ! 937 DO k = 1, 3 938 qp(k) = q(k) + qbid(k, iq) 939 ENDDO 940 ! 941 dyn_blk(:,:,:,:) = (0.d0,0.d0) 942 CALL frc_blk (dyn_blk, qp,tau_blk, nat_blk, & 943 & nr1, nr2, nr3,frc, at_blk, bg_blk, rws, nrws,f_of_q) 944 IF (has_zstar) & 945 CALL rgd_blk(nr1, nr2, nr3, nat_blk, dyn_blk, qp,tau_blk, & 946 epsil, zeu, bg_blk, omega_blk, celldm(1), loto_2d, +1.d0) 947 ! 948 DO na= 1, nat 949 na_blk = itau_blk(na) 950 DO nb= 1, nat 951 nb_blk = itau_blk(nb) 952 ! 953 arg=tpi* ( qp(1) * ( (tau(1, na) -tau_blk(1, na_blk)) - & 954 (tau(1, nb) -tau_blk(1, nb_blk)) ) + & 955 qp(2) * ( (tau(2, na) -tau_blk(2, na_blk)) - & 956 (tau(2, nb) -tau_blk(2, nb_blk)) ) + & 957 qp(3) * ( (tau(3, na) -tau_blk(3, na_blk)) - & 958 (tau(3, nb) -tau_blk(3, nb_blk)) ) ) 959 ! 960 cfac(nb) = CMPLX(COS(arg),SIN(arg), kind= DP)/nsc 961 ! 962 ENDDO ! nb 963 ! 964 DO i = 1, 3 965 DO j = 1, 3 966 ! 967 DO nb= 1, nat 968 nb_blk = itau_blk(nb) 969 dyn(i, j, na, nb) = dyn(i, j, na, nb) + cfac(nb) * & 970 dyn_blk(i, j, na_blk, nb_blk) 971 ENDDO ! nb 972 ! 973 ENDDO ! j 974 ENDDO ! i 975 ENDDO ! na 976 ! 977 ENDDO ! iq 978 ! 979 RETURN 980END SUBROUTINE setupmat 981! 982! 983!---------------------------------------------------------------------- 984SUBROUTINE set_asr (asr, nr1, nr2, nr3, frc, zeu, nat, ibrav, tau) 985 !----------------------------------------------------------------------- 986 ! 987 USE kinds, ONLY : DP 988 USE io_global, ONLY : stdout 989 ! 990 IMPLICIT NONE 991 CHARACTER (LEN= 10), intent(in) :: asr 992 INTEGER, intent(in) :: nr1, nr2, nr3, nat, ibrav 993 REAL(DP), intent(in) :: tau(3, nat) 994 REAL(DP), intent(inout) :: frc(nr1, nr2, nr3, 3, 3, nat, nat), zeu(3,3, nat) 995 ! 996 INTEGER :: axis, n, i, j, na, nb, n1, n2, n3, m, p, k,l,q, r, i1, j1, na1 997 REAL(DP) :: zeu_new(3, 3, nat) 998 REAL(DP), ALLOCATABLE :: frc_new(:,:,:,:,:,:,:) 999 type vector 1000 REAL(DP), pointer :: vec(:,:,:,:,:,:,:) 1001 end type vector 1002 ! 1003 type (vector) u(6*3*nat) 1004 ! These are the "vectors" associated with the sum rules on force-constants 1005 ! 1006 integer :: u_less(6*3*nat), n_less, i_less 1007 ! indices of the vectors u that are not independent to the preceding ones, 1008 ! n_less = number of such vectors, i_less = temporary parameter 1009 ! 1010 integer, allocatable :: ind_v(:,:,:) 1011 REAL(DP), allocatable :: v(:,:) 1012 ! These are the "vectors" associated with symmetry conditions, coded by 1013 ! indicating the positions (i.e. the seven indices) of the non-zero elements (there 1014 ! should be only 2 of them) and the value of that element. We DO so in order 1015 ! to limit the amount of memory used. 1016 ! 1017 REAL(DP), allocatable :: w(:,:,:,:,:,:,:), x(:,:,:,:,:,:,:) 1018 ! temporary vectors and parameters 1019 REAL(DP) :: scal, norm2, sum 1020 ! 1021 REAL(DP) :: zeu_u(6*3, 3, 3, nat) 1022 ! These are the "vectors" associated with the sum rules on effective charges 1023 ! 1024 integer :: zeu_less(6*3), nzeu_less, izeu_less 1025 ! indices of the vectors zeu_u that are not independent to the preceding ones, 1026 ! nzeu_less = number of such vectors, izeu_less = temporary parameter 1027 ! 1028 REAL(DP) :: zeu_w(3, 3, nat), zeu_x(3, 3, nat) 1029 ! temporary vectors 1030 1031 ! Initialization. n is the number of sum rules to be considered (IF asr.ne.'simple') 1032 ! and 'axis' is the rotation axis in the case of a 1D system 1033 ! (i.e. the rotation axis is (Ox) IF axis ='1', (Oy) if axis='2' and (Oz) if axis='3') 1034 ! 1035 if((asr.ne.'simple').and.(asr.ne.'crystal').and.(asr.ne.'one-dim') & 1036 .and.(asr.ne.'zero-dim')) THEN 1037 call errore('set_asr','invalid Acoustic Sum Rule:' // asr, 1) 1038 ENDIF 1039 ! 1040 if(asr.eq.'simple') THEN 1041 ! 1042 ! Simple Acoustic Sum Rule on effective charges 1043 ! 1044 DO i = 1, 3 1045 DO j = 1, 3 1046 sum= 0.0d0 1047 DO na= 1, nat 1048 sum = sum + zeu(i, j, na) 1049 ENDDO 1050 DO na= 1, nat 1051 zeu(i, j, na) = zeu(i, j, na) - sum/nat 1052 ENDDO 1053 ENDDO 1054 ENDDO 1055 ! 1056 ! Simple Acoustic Sum Rule on force constants in REAL space 1057 ! 1058 DO i = 1, 3 1059 DO j = 1, 3 1060 DO na= 1, nat 1061 sum= 0.0d0 1062 DO nb= 1, nat 1063 DO n1= 1, nr1 1064 DO n2 = 1, nr2 1065 DO n3= 1, nr3 1066 sum=sum+frc(n1, n2, n3, i, j, na, nb) 1067 ENDDO 1068 ENDDO 1069 ENDDO 1070 ENDDO 1071 frc(1, 1, 1, i, j, na, na) = frc(1, 1, 1, i, j, na, na) - sum 1072 ! WRITE(6,*) ' na, i, j, sum = ', na, i, j, sum 1073 ENDDO 1074 ENDDO 1075 ENDDO 1076 ! 1077 return 1078 ! 1079 ENDIF 1080 if(asr.eq.'crystal') n=3 1081 if(asr.eq.'one-dim') THEN 1082 ! the direction of periodicity is the rotation axis 1083 ! It will work only IF the crystal axis considered is one of 1084 ! the cartesian axis (typically, ibrav= 1, 6 or 8, or 4 along the 1085 ! z-direction) 1086 IF (nr1*nr2*nr3.eq.1) axis =3 1087 IF ((nr1.ne.1).and.(nr2*nr3.eq.1)) axis = 1 1088 IF ((nr2.ne.1).and.(nr1*nr3.eq.1)) axis =2 1089 IF ((nr3.ne.1).and.(nr1*nr2.eq.1)) axis =3 1090 IF (((nr1.ne.1).and.(nr2.ne.1)).or.((nr2.ne.1).and. & 1091 (nr3.ne.1)).or.((nr1.ne.1).and.(nr3.ne.1))) THEN 1092 call errore('set_asr','too many directions of & 1093 & periodicity in 1D system', axis) 1094 ENDIF 1095 IF ((ibrav.ne.1).and.(ibrav.ne.6).and.(ibrav.ne.8).and. & 1096 ((ibrav.ne.4).or.(axis.ne.3)) ) THEN 1097 WRITE(stdout,*) 'asr: rotational axis may be wrong' 1098 ENDIF 1099 WRITE(stdout,'("asr rotation axis in 1D system= ",I4)') axis 1100 n=4 1101 ENDIF 1102 if(asr.eq.'zero-dim') n=6 1103 ! 1104 ! Acoustic Sum Rule on effective charges 1105 ! 1106 ! generating the vectors of the orthogonal of the subspace to project 1107 ! the effective charges matrix on 1108 ! 1109 zeu_u(:,:,:,:) = 0.0d0 1110 DO i = 1, 3 1111 DO j = 1, 3 1112 DO na= 1, nat 1113 zeu_new(i, j, na) =zeu(i, j, na) 1114 ENDDO 1115 ENDDO 1116 ENDDO 1117 ! 1118 p = 0 1119 DO i = 1, 3 1120 DO j = 1, 3 1121 ! These are the 3*3 vectors associated with the 1122 ! translational acoustic sum rules 1123 p = p + 1 1124 zeu_u(p, i, j,:) = 1.0d0 1125 ! 1126 ENDDO 1127 ENDDO 1128 ! 1129 IF (n.eq.4) THEN 1130 DO i = 1, 3 1131 ! These are the 3 vectors associated with the 1132 ! single rotational sum rule (1D system) 1133 p = p + 1 1134 DO na= 1, nat 1135 zeu_u(p, i, MOD(axis, 3) + 1, na) =-tau(MOD(axis+ 1, 3) + 1, na) 1136 zeu_u(p, i, MOD(axis+ 1, 3) + 1, na) =tau(MOD(axis, 3) + 1, na) 1137 ENDDO 1138 ! 1139 ENDDO 1140 ENDIF 1141 ! 1142 IF (n.eq.6) THEN 1143 DO i = 1, 3 1144 DO j = 1, 3 1145 ! These are the 3*3 vectors associated with the 1146 ! three rotational sum rules (0D system - typ. molecule) 1147 p = p + 1 1148 DO na= 1, nat 1149 zeu_u(p, i, MOD(j, 3) + 1, na) =-tau(MOD(j+ 1, 3) + 1, na) 1150 zeu_u(p, i, MOD(j+ 1, 3) + 1, na) =tau(MOD(j, 3) + 1, na) 1151 ENDDO 1152 ! 1153 ENDDO 1154 ENDDO 1155 ENDIF 1156 ! 1157 ! Gram-Schmidt orthonormalization of the set of vectors created. 1158 ! 1159 nzeu_less = 0 1160 DO k = 1, p 1161 zeu_w(:,:,:) =zeu_u(k,:,:,:) 1162 zeu_x(:,:,:) =zeu_u(k,:,:,:) 1163 DO q= 1, k- 1 1164 r = 1 1165 DO izeu_less = 1, nzeu_less 1166 IF (zeu_less(izeu_less).eq.q) r = 0 1167 ENDDO 1168 IF (r.ne.0) THEN 1169 call sp_zeu(zeu_x, zeu_u(q,:,:,:), nat, scal) 1170 zeu_w(:,:,:) = zeu_w(:,:,:) - scal* zeu_u(q,:,:,:) 1171 ENDIF 1172 ENDDO 1173 call sp_zeu(zeu_w, zeu_w, nat, norm2) 1174 IF (norm2.gt.1.0d-16) THEN 1175 zeu_u(k,:,:,:) = zeu_w(:,:,:) / DSQRT(norm2) 1176 ELSE 1177 nzeu_less =nzeu_less+ 1 1178 zeu_less(nzeu_less) =k 1179 ENDIF 1180 ENDDO 1181 ! 1182 ! Projection of the effective charge "vector" on the orthogonal of the 1183 ! subspace of the vectors verifying the sum rules 1184 ! 1185 zeu_w(:,:,:) = 0.0d0 1186 DO k = 1, p 1187 r = 1 1188 DO izeu_less = 1, nzeu_less 1189 IF (zeu_less(izeu_less).eq.k) r = 0 1190 ENDDO 1191 IF (r.ne.0) THEN 1192 zeu_x(:,:,:) =zeu_u(k,:,:,:) 1193 call sp_zeu(zeu_x, zeu_new, nat, scal) 1194 zeu_w(:,:,:) = zeu_w(:,:,:) + scal*zeu_u(k,:,:,:) 1195 ENDIF 1196 ENDDO 1197 ! 1198 ! Final substraction of the former projection to the initial zeu, to get 1199 ! the new "projected" zeu 1200 ! 1201 zeu_new(:,:,:) =zeu_new(:,:,:) - zeu_w(:,:,:) 1202 call sp_zeu(zeu_w, zeu_w, nat, norm2) 1203 WRITE(stdout,'("Norm of the difference between old and new effective ", & 1204 & "charges: ", F25.20)') SQRT(norm2) 1205 ! 1206 ! Check projection 1207 ! 1208 !WRITE(6,'("Check projection of zeu")') 1209 !DO k = 1, p 1210 ! zeu_x(:,:,:) =zeu_u(k,:,:,:) 1211 ! call sp_zeu(zeu_x, zeu_new, nat, scal) 1212 ! IF (DABS(scal).gt.1d- 10) WRITE(6,'("k = ",I8," zeu_new|zeu_u(k) = ", F15.10)') k, scal 1213 !ENDDO 1214 ! 1215 DO i = 1, 3 1216 DO j = 1, 3 1217 DO na= 1, nat 1218 zeu(i, j, na) =zeu_new(i, j, na) 1219 ENDDO 1220 ENDDO 1221 ENDDO 1222 ! 1223 ! Acoustic Sum Rule on force constants 1224 ! 1225 ! 1226 ! generating the vectors of the orthogonal of the subspace to project 1227 ! the force-constants matrix on 1228 ! 1229 DO k = 1, 18*nat 1230 ALLOCATE(u(k) % vec(nr1, nr2, nr3, 3, 3, nat, nat)) 1231 u(k) % vec (:,:,:,:,:,:,:) = 0.0d0 1232 ENDDO 1233 ALLOCATE (frc_new(nr1, nr2, nr3, 3, 3, nat, nat)) 1234 DO i = 1, 3 1235 DO j = 1, 3 1236 DO na= 1, nat 1237 DO nb= 1, nat 1238 DO n1= 1, nr1 1239 DO n2 = 1, nr2 1240 DO n3= 1, nr3 1241 frc_new(n1, n2, n3, i, j, na, nb) =frc(n1, n2, n3, i, j, na, nb) 1242 ENDDO 1243 ENDDO 1244 ENDDO 1245 ENDDO 1246 ENDDO 1247 ENDDO 1248 ENDDO 1249 ! 1250 p = 0 1251 DO i = 1, 3 1252 DO j = 1, 3 1253 DO na= 1, nat 1254 ! These are the 3*3*nat vectors associated with the 1255 ! translational acoustic sum rules 1256 p = p + 1 1257 u(p) % vec (:,:,:, i, j, na,:) = 1.0d0 1258 ! 1259 ENDDO 1260 ENDDO 1261 ENDDO 1262 ! 1263 IF (n.eq.4) THEN 1264 DO i = 1, 3 1265 DO na= 1, nat 1266 ! These are the 3*nat vectors associated with the 1267 ! single rotational sum rule (1D system) 1268 p = p + 1 1269 DO nb= 1, nat 1270 u(p) % vec (:,:,:, i, MOD(axis, 3) + 1, na, nb) =-tau(MOD(axis+ 1, 3) + 1, nb) 1271 u(p) % vec (:,:,:, i, MOD(axis+ 1, 3) + 1, na, nb) =tau(MOD(axis, 3) + 1, nb) 1272 ENDDO 1273 ! 1274 ENDDO 1275 ENDDO 1276 ENDIF 1277 ! 1278 IF (n.eq.6) THEN 1279 DO i = 1, 3 1280 DO j = 1, 3 1281 DO na= 1, nat 1282 ! These are the 3*3*nat vectors associated with the 1283 ! three rotational sum rules (0D system - typ. molecule) 1284 p = p + 1 1285 DO nb= 1, nat 1286 u(p) % vec (:,:,:, i, MOD(j, 3) + 1, na, nb) =-tau(MOD(j+ 1, 3) + 1, nb) 1287 u(p) % vec (:,:,:, i, MOD(j+ 1, 3) + 1, na, nb) =tau(MOD(j, 3) + 1, nb) 1288 ENDDO 1289 ! 1290 ENDDO 1291 ENDDO 1292 ENDDO 1293 ENDIF 1294 ! 1295 ALLOCATE (ind_v(9*nat*nat*nr1*nr2*nr3, 2,7), v(9*nat*nat*nr1*nr2*nr3, 2) ) 1296 m= 0 1297 DO i = 1, 3 1298 DO j = 1, 3 1299 DO na= 1, nat 1300 DO nb= 1, nat 1301 DO n1= 1, nr1 1302 DO n2 = 1, nr2 1303 DO n3= 1, nr3 1304 ! These are the vectors associated with the symmetry constraints 1305 q= 1 1306 l= 1 1307 DO while((l.le.m).and.(q.ne.0)) 1308 IF ((ind_v(l, 1, 1).eq.n1).and.(ind_v(l, 1, 2).eq.n2).and. & 1309 (ind_v(l, 1, 3).eq.n3).and.(ind_v(l, 1,4).eq.i).and. & 1310 (ind_v(l, 1,5).eq.j).and.(ind_v(l, 1,6).eq.na).and. & 1311 (ind_v(l, 1,7).eq.nb)) q= 0 1312 IF ((ind_v(l, 2, 1).eq.n1).and.(ind_v(l, 2, 2).eq.n2).and. & 1313 (ind_v(l, 2, 3).eq.n3).and.(ind_v(l, 2,4).eq.i).and. & 1314 (ind_v(l, 2,5).eq.j).and.(ind_v(l, 2,6).eq.na).and. & 1315 (ind_v(l, 2,7).eq.nb)) q= 0 1316 l=l+ 1 1317 ENDDO 1318 IF ((n1.eq.MOD(nr1 + 1-n1, nr1) + 1).and.(n2.eq.MOD(nr2 + 1-n2, nr2) + 1) & 1319 .and.(n3.eq.MOD(nr3 + 1-n3, nr3) + 1).and.(i.eq.j).and.(na.eq.nb)) q= 0 1320 IF (q.ne.0) THEN 1321 m=m+ 1 1322 ind_v(m, 1, 1) =n1 1323 ind_v(m, 1, 2) =n2 1324 ind_v(m, 1, 3) =n3 1325 ind_v(m, 1,4) =i 1326 ind_v(m, 1,5) =j 1327 ind_v(m, 1,6) =na 1328 ind_v(m, 1,7) =nb 1329 v(m, 1) = 1.0d0/DSQRT(2.0d0) 1330 ind_v(m, 2, 1) =MOD(nr1 + 1-n1, nr1) + 1 1331 ind_v(m, 2, 2) =MOD(nr2 + 1-n2, nr2) + 1 1332 ind_v(m, 2, 3) =MOD(nr3 + 1-n3, nr3) + 1 1333 ind_v(m, 2,4) =j 1334 ind_v(m, 2,5) =i 1335 ind_v(m, 2,6) =nb 1336 ind_v(m, 2,7) =na 1337 v(m, 2) =- 1.0d0/DSQRT(2.0d0) 1338 ENDIF 1339 ENDDO 1340 ENDDO 1341 ENDDO 1342 ENDDO 1343 ENDDO 1344 ENDDO 1345 ENDDO 1346 ! 1347 ! Gram-Schmidt orthonormalization of the set of vectors created. 1348 ! Note that the vectors corresponding to symmetry constraints are already 1349 ! orthonormalized by construction. 1350 ! 1351 n_less = 0 1352 ALLOCATE (w(nr1, nr2, nr3, 3, 3, nat, nat), x(nr1, nr2, nr3,3,3, nat, nat)) 1353 DO k = 1, p 1354 w(:,:,:,:,:,:,:) =u(k) % vec (:,:,:,:,:,:,:) 1355 x(:,:,:,:,:,:,:) =u(k) % vec (:,:,:,:,:,:,:) 1356 DO l= 1,m 1357 ! 1358 call sp2(x,v(l,:), ind_v(l,:,:), nr1, nr2, nr3, nat, scal) 1359 DO r = 1, 2 1360 n1=ind_v(l, r, 1) 1361 n2 =ind_v(l, r, 2) 1362 n3=ind_v(l, r, 3) 1363 i =ind_v(l, r,4) 1364 j =ind_v(l, r,5) 1365 na=ind_v(l, r,6) 1366 nb=ind_v(l, r,7) 1367 w(n1, n2, n3, i, j, na, nb) =w(n1, n2, n3, i, j, na, nb) -scal*v(l, r) 1368 ENDDO 1369 ENDDO 1370 IF (k.le.(9*nat)) THEN 1371 na1=MOD(k, nat) 1372 IF (na1.eq.0) na1=nat 1373 j1=MOD((k-na1)/nat, 3) + 1 1374 i1=MOD((((k-na1)/nat) -j1 + 1)/3, 3) + 1 1375 ELSE 1376 q=k-9*nat 1377 IF (n.eq.4) THEN 1378 na1=MOD(q, nat) 1379 IF (na1.eq.0) na1=nat 1380 i1=MOD((q-na1)/nat, 3) + 1 1381 ELSE 1382 na1=MOD(q, nat) 1383 IF (na1.eq.0) na1=nat 1384 j1=MOD((q-na1)/nat, 3) + 1 1385 i1=MOD((((q-na1)/nat) -j1 + 1)/3, 3) + 1 1386 ENDIF 1387 ENDIF 1388 DO q= 1, k- 1 1389 r = 1 1390 DO i_less = 1, n_less 1391 IF (u_less(i_less).eq.q) r = 0 1392 ENDDO 1393 IF (r.ne.0) THEN 1394 call sp3(x,u(q) % vec (:,:,:,:,:,:,:), i1, na1, nr1, nr2, nr3, nat, scal) 1395 w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) - scal* u(q) % vec (:,:,:,:,:,:,:) 1396 ENDIF 1397 ENDDO 1398 call sp1(w,w, nr1, nr2, nr3, nat, norm2) 1399 IF (norm2.gt.1.0d-16) THEN 1400 u(k) % vec (:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) / DSQRT(norm2) 1401 ELSE 1402 n_less =n_less+ 1 1403 u_less(n_less) =k 1404 ENDIF 1405 ENDDO 1406 ! 1407 ! Projection of the force-constants "vector" on the orthogonal of the 1408 ! subspace of the vectors verifying the sum rules and symmetry contraints 1409 ! 1410 w(:,:,:,:,:,:,:) = 0.0d0 1411 DO l= 1,m 1412 call sp2(frc_new,v(l,:), ind_v(l,:,:), nr1, nr2, nr3, nat, scal) 1413 DO r = 1, 2 1414 n1=ind_v(l, r, 1) 1415 n2 =ind_v(l, r, 2) 1416 n3=ind_v(l, r, 3) 1417 i =ind_v(l, r,4) 1418 j =ind_v(l, r,5) 1419 na=ind_v(l, r,6) 1420 nb=ind_v(l, r,7) 1421 w(n1, n2, n3, i, j, na, nb) =w(n1, n2, n3, i, j, na, nb) +scal*v(l, r) 1422 ENDDO 1423 ENDDO 1424 DO k = 1, p 1425 r = 1 1426 DO i_less = 1, n_less 1427 IF (u_less(i_less).eq.k) r = 0 1428 ENDDO 1429 IF (r.ne.0) THEN 1430 x(:,:,:,:,:,:,:) =u(k) % vec (:,:,:,:,:,:,:) 1431 call sp1(x,frc_new, nr1, nr2, nr3, nat, scal) 1432 w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) + scal*u(k)%vec(:,:,:,:,:,:,:) 1433 ENDIF 1434 DEALLOCATE(u(k) % vec) 1435 ENDDO 1436 ! 1437 ! Final substraction of the former projection to the initial frc, to get 1438 ! the new "projected" frc 1439 ! 1440 frc_new(:,:,:,:,:,:,:) =frc_new(:,:,:,:,:,:,:) - w(:,:,:,:,:,:,:) 1441 call sp1(w,w, nr1, nr2, nr3, nat, norm2) 1442 WRITE(stdout,'("Norm of the difference between old and new force-constants:",& 1443 & F25.20)') SQRT(norm2) 1444 ! 1445 ! Check projection 1446 ! 1447 !WRITE(6,'("Check projection IFC")') 1448 !DO l= 1,m 1449 ! call sp2(frc_new,v(l,:), ind_v(l,:,:), nr1, nr2, nr3, nat, scal) 1450 ! IF (DABS(scal).gt.1d- 10) WRITE(6,'("l= ",I8," frc_new|v(l) = ", F15.10)') l, scal 1451 !ENDDO 1452 !DO k = 1, p 1453 ! x(:,:,:,:,:,:,:) =u(k) % vec (:,:,:,:,:,:,:) 1454 ! call sp1(x,frc_new, nr1, nr2, nr3, nat, scal) 1455 ! IF (DABS(scal).gt.1d- 10) WRITE(6,'("k = ",I8," frc_new|u(k) = ", F15.10)') k, scal 1456 ! DEALLOCATE(u(k) % vec) 1457 !ENDDO 1458 ! 1459 DO i = 1, 3 1460 DO j = 1, 3 1461 DO na= 1, nat 1462 DO nb= 1, nat 1463 DO n1= 1, nr1 1464 DO n2 = 1, nr2 1465 DO n3= 1, nr3 1466 frc(n1, n2, n3, i, j, na, nb) =frc_new(n1, n2, n3, i, j, na, nb) 1467 ENDDO 1468 ENDDO 1469 ENDDO 1470 ENDDO 1471 ENDDO 1472 ENDDO 1473 ENDDO 1474 DEALLOCATE (x, w) 1475 DEALLOCATE (v, ind_v) 1476 DEALLOCATE (frc_new) 1477 ! 1478 return 1479end subroutine set_asr 1480! 1481!---------------------------------------------------------------------- 1482subroutine sp_zeu(zeu_u, zeu_v, nat, scal) 1483 !----------------------------------------------------------------------- 1484 ! 1485 ! does the scalar product of two effective charges matrices zeu_u and zeu_v 1486 ! (considered as vectors in the R^(3*3*nat) space, and coded in the usual way) 1487 ! 1488 USE kinds, ONLY: DP 1489 implicit none 1490 integer i, j, na, nat 1491 REAL(DP) zeu_u(3, 3, nat) 1492 REAL(DP) zeu_v(3, 3, nat) 1493 REAL(DP) scal 1494 ! 1495 ! 1496 scal= 0.0d0 1497 DO i = 1, 3 1498 DO j = 1, 3 1499 DO na= 1, nat 1500 scal=scal+zeu_u(i, j, na) *zeu_v(i, j, na) 1501 ENDDO 1502 ENDDO 1503 ENDDO 1504 ! 1505 return 1506 ! 1507end subroutine sp_zeu 1508! 1509! 1510!---------------------------------------------------------------------- 1511subroutine sp1(u,v, nr1, nr2, nr3, nat, scal) 1512 !----------------------------------------------------------------------- 1513 ! 1514 ! does the scalar product of two force-constants matrices u and v (considered as 1515 ! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space, and coded in the usual way) 1516 ! 1517 USE kinds, ONLY: DP 1518 implicit none 1519 integer nr1, nr2, nr3, i, j, na, nb, n1, n2, n3, nat 1520 REAL(DP) u(nr1, nr2, nr3, 3, 3, nat, nat) 1521 REAL(DP) v(nr1, nr2, nr3, 3, 3, nat, nat) 1522 REAL(DP) scal 1523 ! 1524 ! 1525 scal= 0.0d0 1526 DO i = 1, 3 1527 DO j = 1, 3 1528 DO na= 1, nat 1529 DO nb= 1, nat 1530 DO n1= 1, nr1 1531 DO n2 = 1, nr2 1532 DO n3= 1, nr3 1533 scal=scal+u(n1, n2, n3, i, j, na, nb) *v(n1, n2, n3, i, j, na, nb) 1534 ENDDO 1535 ENDDO 1536 ENDDO 1537 ENDDO 1538 ENDDO 1539 ENDDO 1540 ENDDO 1541 ! 1542 return 1543 ! 1544end subroutine sp1 1545! 1546!---------------------------------------------------------------------- 1547subroutine sp2(u,v, ind_v, nr1, nr2, nr3, nat, scal) 1548 !----------------------------------------------------------------------- 1549 ! 1550 ! does the scalar product of two force-constants matrices u and v (considered as 1551 ! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space). u is coded in the usual way 1552 ! but v is coded as EXPlained when defining the vectors corresponding to the 1553 ! symmetry constraints 1554 ! 1555 USE kinds, ONLY: DP 1556 implicit none 1557 integer nr1, nr2, nr3, i, nat 1558 REAL(DP) u(nr1, nr2, nr3, 3, 3, nat, nat) 1559 integer ind_v(2,7) 1560 REAL(DP) v(2) 1561 REAL(DP) scal 1562 ! 1563 ! 1564 scal= 0.0d0 1565 DO i = 1, 2 1566 scal=scal+u(ind_v(i, 1), ind_v(i, 2), ind_v(i, 3), ind_v(i,4), ind_v(i,5), ind_v(i,6), & 1567 ind_v(i,7)) *v(i) 1568 ENDDO 1569 ! 1570 return 1571 ! 1572end subroutine sp2 1573! 1574!---------------------------------------------------------------------- 1575subroutine sp3(u,v, i, na, nr1, nr2, nr3, nat, scal) 1576 !----------------------------------------------------------------------- 1577 ! 1578 ! like sp1, but in the particular case when u is one of the u(k)%vec 1579 ! defined in set_asr (before orthonormalization). In this case most of the 1580 ! terms are zero (the ones that are not are characterized by i and na), so 1581 ! that a lot of computer time can be saved (during Gram-Schmidt). 1582 ! 1583 USE kinds, ONLY: DP 1584 implicit none 1585 integer nr1, nr2, nr3, i, j, na, nb, n1, n2, n3, nat 1586 REAL(DP) u(nr1, nr2, nr3, 3, 3, nat, nat) 1587 REAL(DP) v(nr1, nr2, nr3, 3, 3, nat, nat) 1588 REAL(DP) scal 1589 ! 1590 ! 1591 scal= 0.0d0 1592 DO j = 1, 3 1593 DO nb= 1, nat 1594 DO n1= 1, nr1 1595 DO n2 = 1, nr2 1596 DO n3= 1, nr3 1597 scal=scal+u(n1, n2, n3, i, j, na, nb) *v(n1, n2, n3, i, j, na, nb) 1598 ENDDO 1599 ENDDO 1600 ENDDO 1601 ENDDO 1602 ENDDO 1603 ! 1604 return 1605 ! 1606end subroutine sp3 1607! 1608!----------------------------------------------------------------------- 1609SUBROUTINE q_gen(nsc,qbid, at_blk, bg_blk, at, bg) 1610 !----------------------------------------------------------------------- 1611 ! generate list of q (qbid) that are G-vectors of the supercell 1612 ! but not of the bulk 1613 ! 1614 USE kinds, ONLY : DP 1615 ! 1616 IMPLICIT NONE 1617 INTEGER :: nsc 1618 REAL(DP) qbid(3, nsc), at_blk(3, 3), bg_blk(3, 3), at(3,3), bg(3,3) 1619 ! 1620 INTEGER, PARAMETER:: nr1=4, nr2 =4, nr3=4, & 1621 nrm=(2*nr1 + 1) * (2*nr2 + 1) * (2*nr3 + 1) 1622 REAL(DP), PARAMETER:: eps = 1.0d-7 1623 INTEGER :: i, j, k, i1, i2, i3, idum(nrm), iq 1624 REAL(DP) :: qnorm(nrm), qbd(3, nrm) ,qwork(3), delta 1625 LOGICAL lbho 1626 ! 1627 i = 0 1628 DO i1=-nr1, nr1 1629 DO i2 =-nr2, nr2 1630 DO i3=-nr3, nr3 1631 i = i + 1 1632 DO j = 1, 3 1633 qwork(j) = i1*bg(j, 1) + i2*bg(j, 2) + i3*bg(j, 3) 1634 ENDDO ! j 1635 ! 1636 qnorm(i) = qwork(1)**2 + qwork(2)**2 + qwork(3) **2 1637 ! 1638 DO j = 1, 3 1639 ! 1640 qbd(j, i) = at_blk(1, j) *qwork(1) + & 1641 at_blk(2, j) *qwork(2) + & 1642 at_blk(3, j) *qwork(3) 1643 ENDDO ! j 1644 ! 1645 idum(i) = 1 1646 ! 1647 ENDDO ! i3 1648 ENDDO ! i2 1649 ENDDO ! i1 1650 ! 1651 DO i = 1, nrm- 1 1652 IF (idum(i).EQ.1) THEN 1653 DO j =i + 1, nrm 1654 IF (idum(j).EQ.1) THEN 1655 lbho=.TRUE. 1656 DO k = 1, 3 1657 delta = qbd(k, i) -qbd(k, j) 1658 lbho = lbho.AND. (ABS(NINT(delta) -delta)<eps) 1659 ENDDO ! k 1660 IF (lbho) THEN 1661 IF(qnorm(i).GT.qnorm(j)) THEN 1662 qbd(1, i) = qbd(1, j) 1663 qbd(2, i) = qbd(2, j) 1664 qbd(3, i) = qbd(3, j) 1665 qnorm(i) = qnorm(j) 1666 ENDIF 1667 idum(j) = 0 1668 ENDIF 1669 ENDIF 1670 ENDDO ! j 1671 ENDIF 1672 ENDDO ! i 1673 ! 1674 iq = 0 1675 DO i = 1, nrm 1676 IF (idum(i).EQ.1) THEN 1677 iq=iq+ 1 1678 qbid(1, iq) = bg_blk(1, 1) *qbd(1, i) + & 1679 bg_blk(1, 2) *qbd(2, i) + & 1680 bg_blk(1, 3) *qbd(3, i) 1681 qbid(2, iq) = bg_blk(2, 1) *qbd(1, i) + & 1682 bg_blk(2, 2) *qbd(2, i) + & 1683 bg_blk(2, 3) *qbd(3, i) 1684 qbid(3, iq) = bg_blk(3, 1) *qbd(1, i) + & 1685 bg_blk(3, 2) *qbd(2, i) + & 1686 bg_blk(3, 3) *qbd(3, i) 1687 ENDIF 1688 ENDDO ! i 1689 ! 1690 IF (iq.NE.nsc) CALL errore('q_gen',' probably nr1, nr2, nr3 too small ', iq) 1691 RETURN 1692END SUBROUTINE q_gen 1693! 1694!----------------------------------------------------------------------- 1695SUBROUTINE check_at(at, bg_blk, alat, omega) 1696 !----------------------------------------------------------------------- 1697 ! 1698 USE kinds, ONLY : DP 1699 USE io_global, ONLY : stdout 1700 ! 1701 IMPLICIT NONE 1702 ! 1703 REAL(DP) :: at(3, 3), bg_blk(3, 3), alat, omega 1704 REAL(DP) :: work(3, 3) 1705 INTEGER :: i, j 1706 REAL(DP), PARAMETER :: small= 1.d-6 1707 ! 1708 work(:,:) = at(:,:) 1709 CALL cryst_to_cart(3,work, bg_blk,- 1) 1710 ! 1711 DO j = 1, 3 1712 DO i = 1, 3 1713 IF ( ABS(work(i, j) -NINT(work(i, j))) > small) THEN 1714 WRITE (stdout,'(3f9.4)') work(:,:) 1715 CALL errore ('check_at','at not multiple of at_blk', 1) 1716 ENDIF 1717 ENDDO 1718 ENDDO 1719 ! 1720 omega =alat**3 * ABS(at(1, 1) * (at(2, 2) *at(3, 3) -at(3, 2) *at(2, 3)) - & 1721 at(1, 2) * (at(2, 1) *at(3, 3) -at(2, 3) *at(3, 1)) + & 1722 at(1, 3) * (at(2, 1) *at(3, 2) -at(2, 2) *at(3, 1))) 1723 ! 1724 RETURN 1725END SUBROUTINE check_at 1726! 1727!----------------------------------------------------------------------- 1728SUBROUTINE set_tau (nat, nat_blk, at, at_blk, tau, tau_blk, & 1729 ityp, ityp_blk, itau_blk) 1730 !----------------------------------------------------------------------- 1731 ! 1732 USE kinds, ONLY : DP 1733 ! 1734 IMPLICIT NONE 1735 INTEGER nat, nat_blk, ityp(nat), ityp_blk(nat_blk), itau_blk(nat) 1736 REAL(DP) at(3, 3), at_blk(3, 3),tau(3, nat),tau_blk(3, nat_blk) 1737 ! 1738 REAL(DP) bg(3, 3), r(3) ! work vectors 1739 INTEGER i, i1, i2, i3, na, na_blk 1740 REAL(DP) small 1741 INTEGER NN1,NN2,NN3 1742 PARAMETER (NN1=8, NN2 =8, NN3=8, small= 1.d-8) 1743 ! 1744 CALL recips (at(1, 1), at(1, 2), at(1, 3), bg(1, 1), bg(1, 2), bg(1, 3)) 1745 ! 1746 na = 0 1747 ! 1748 DO i1 = -NN1,NN1 1749 DO i2 = -NN2,NN2 1750 DO i3 = -NN3,NN3 1751 r(1) = i1*at_blk(1, 1) + i2*at_blk(1, 2) + i3*at_blk(1, 3) 1752 r(2) = i1*at_blk(2, 1) + i2*at_blk(2, 2) + i3*at_blk(2, 3) 1753 r(3) = i1*at_blk(3, 1) + i2*at_blk(3, 2) + i3*at_blk(3, 3) 1754 CALL cryst_to_cart(1, r, bg,- 1) 1755 ! 1756 IF ( r(1).GT.-small .AND. r(1)<1.d0-small .AND. & 1757 r(2).GT.-small .AND. r(2)<1.d0-small .AND. & 1758 r(3).GT.-small .AND. r(3)<1.d0-small ) THEN 1759 CALL cryst_to_cart(1, r, at,+ 1) 1760 ! 1761 DO na_blk = 1, nat_blk 1762 na = na + 1 1763 IF (na.GT.nat) CALL errore('set_tau','too many atoms', na) 1764 tau(1, na) = tau_blk(1, na_blk) + r(1) 1765 tau(2, na) = tau_blk(2, na_blk) + r(2) 1766 tau(3, na) = tau_blk(3, na_blk) + r(3) 1767 ityp(na) = ityp_blk(na_blk) 1768 itau_blk(na) = na_blk 1769 ENDDO 1770 ! 1771 ENDIF 1772 ! 1773 ENDDO 1774 ENDDO 1775 ENDDO 1776 ! 1777 IF (na.NE.nat) CALL errore('set_tau','too few atoms: increase NNs', na) 1778 ! 1779 RETURN 1780END SUBROUTINE set_tau 1781! 1782!----------------------------------------------------------------------- 1783SUBROUTINE read_tau & 1784 (nat, nat_blk, ntyp, bg_blk, tau, tau_blk, ityp, itau_blk) 1785 !--------------------------------------------------------------------- 1786 ! 1787 USE kinds, ONLY : DP 1788 USE io_global, ONLY : ionode_id, ionode 1789 USE mp, ONLY : mp_bcast 1790 USE mp_world, ONLY : world_comm 1791 ! 1792 IMPLICIT NONE 1793 ! 1794 INTEGER nat, nat_blk, ntyp, ityp(nat), itau_blk(nat) 1795 REAL(DP) bg_blk(3, 3),tau(3, nat),tau_blk(3, nat_blk) 1796 ! 1797 REAL(DP) r(3) ! work vectors 1798 INTEGER i, na, na_blk 1799 ! 1800 REAL(DP) small 1801 PARAMETER ( small = 1.d-6 ) 1802 ! 1803 DO na= 1, nat 1804 IF (ionode) READ(5,*) (tau(i, na), i = 1, 3), ityp(na) 1805 CALL mp_bcast(tau(:, na), ionode_id, world_comm) 1806 CALL mp_bcast(ityp(na), ionode_id, world_comm) 1807 IF (ityp(na).LE.0 .OR. ityp(na) > ntyp) & 1808 CALL errore('read_tau',' wrong atomic type', na) 1809 DO na_blk = 1, nat_blk 1810 r(1) = tau(1, na) - tau_blk(1, na_blk) 1811 r(2) = tau(2, na) - tau_blk(2, na_blk) 1812 r(3) = tau(3, na) - tau_blk(3, na_blk) 1813 CALL cryst_to_cart(1, r, bg_blk,- 1) 1814 IF (ABS( r(1) -NINT(r(1)) ) < small .AND. & 1815 ABS( r(2) -NINT(r(2)) ) < small .AND. & 1816 ABS( r(3) -NINT(r(3)) ) < small ) THEN 1817 itau_blk(na) = na_blk 1818 go to 999 1819 ENDIF 1820 ENDDO 1821 CALL errore ('read_tau',' wrong atomic position ', na) 1822999 CONTINUE 1823 ENDDO 1824 ! 1825 RETURN 1826END SUBROUTINE read_tau 1827! 1828!----------------------------------------------------------------------- 1829! 1830!--------------------------------------------------------------------- 1831! 1832!----------------------------------------------------------------------- 1833subroutine setgam (q, gam, nat, at, bg,tau, itau_blk, nsc, alat, & 1834 & gam_blk, nat_blk, at_blk, bg_blk,tau_blk, omega_blk, & 1835 & frcg, nr1, nr2, nr3, rws, nrws) 1836 !----------------------------------------------------------------------- 1837 ! compute the dynamical matrix (the analytic part only) 1838 ! 1839 USE kinds, ONLY : DP 1840 USE constants, ONLY : tpi 1841 implicit none 1842 ! 1843 ! I/O variables 1844 ! 1845 integer :: nr1, nr2, nr3, nat, nat_blk, & 1846 nsc, nrws, itau_blk(nat) 1847 REAL(DP) :: q(3), tau(3, nat), at(3, 3), bg(3, 3), alat, rws(0:3, nrws) 1848 REAL(DP) :: tau_blk(3, nat_blk), at_blk(3, 3), bg_blk(3, 3), omega_blk, & 1849 frcg(nr1, nr2, nr3, 3, 3, nat_blk, nat_blk) 1850 COMPLEX(DP) :: gam_blk(3, 3, nat_blk, nat_blk),f_of_q(3, 3, nat, nat) 1851 COMPLEX(DP) :: gam(3, 3, nat, nat) 1852 ! 1853 ! local variables 1854 ! 1855 REAL(DP) :: arg 1856 complex(DP) :: cfac(nat) 1857 integer :: i, j, k, na, nb, na_blk, nb_blk, iq 1858 REAL(DP) :: qp(3), qbid(3, nsc) ! automatic array 1859 ! 1860 ! 1861 call q_gen(nsc,qbid, at_blk, bg_blk, at, bg) 1862 ! 1863 f_of_q=(0.0_DP,0.0_DP) 1864 DO iq= 1, nsc 1865 ! 1866 DO k = 1, 3 1867 qp(k) = q(k) + qbid(k, iq) 1868 ENDDO 1869 ! 1870 gam_blk(:,:,:,:) = (0.d0,0.d0) 1871 CALL frc_blk (gam_blk, qp,tau_blk, nat_blk, & 1872 nr1, nr2, nr3,frcg, at_blk, bg_blk, rws, nrws,f_of_q) 1873 ! 1874 DO na= 1, nat 1875 na_blk = itau_blk(na) 1876 DO nb= 1, nat 1877 nb_blk = itau_blk(nb) 1878 ! 1879 arg = tpi * ( qp(1) * ( (tau(1, na) -tau_blk(1, na_blk)) - & 1880 (tau(1, nb) -tau_blk(1, nb_blk)) ) + & 1881 qp(2) * ( (tau(2, na) -tau_blk(2, na_blk)) - & 1882 (tau(2, nb) -tau_blk(2, nb_blk)) ) + & 1883 qp(3) * ( (tau(3, na) -tau_blk(3, na_blk)) - & 1884 (tau(3, nb) -tau_blk(3, nb_blk)) ) ) 1885 ! 1886 cfac(nb) = CMPLX(cos(arg), sin(arg), kind=dp)/nsc 1887 ! 1888 ENDDO ! nb 1889 DO nb= 1, nat 1890 DO i = 1, 3 1891 DO j = 1, 3 1892 nb_blk = itau_blk(nb) 1893 gam(i, j, na, nb) = gam(i, j, na, nb) + cfac(nb) * & 1894 gam_blk(i, j, na_blk, nb_blk) 1895 ENDDO ! j 1896 ENDDO ! i 1897 ENDDO ! nb 1898 ENDDO ! na 1899 ! 1900 ENDDO ! iq 1901 ! 1902 return 1903end subroutine setgam 1904! 1905!-------------------------------------------------------------------- 1906function dos_gam (nbndx, nq, jbnd, gamma, et, ef) 1907 !-------------------------------------------------------------------- 1908 ! calculates weights with the tetrahedron method (Bloechl version) 1909 ! this subroutine is based on tweights.f90 belonging to PW 1910 ! it calculates a2F on the surface of given frequency <=> histogram 1911 ! Band index means the frequency mode here 1912 ! and "et" means the frequency(mode,q-point) 1913 ! 1914 USE kinds, ONLY: DP 1915 USE parameters 1916 USE ktetra, ONLY : ntetra, tetra 1917 implicit none 1918 ! 1919 integer :: nq, nbndx, jbnd 1920 REAL(DP) :: et(nbndx, nq), gamma(nbndx, nq), func 1921 1922 REAL(DP) :: ef 1923 REAL(DP) :: e1, e2, e3, e4, c1, c2, c3, c4, etetra(4) 1924 integer :: ik, ibnd, nt, nk, ns, i, ik1, ik2, ik3, ik4, itetra(4) 1925 1926 REAL(DP) :: f12,f13,f14,f23,f24,f34, f21,f31,f41,f42,f32,f43 1927 REAL(DP) :: P1,P2,P3,P4, G, o13, Y1,Y2,Y3,Y4, eps,vol, Tint 1928 REAL(DP) :: dos_gam 1929 1930 Tint = 0.0d0 1931 o13 = 1.0_dp/3.0_dp 1932 eps = 1.0d-14 1933 vol = 1.0d0/ntetra 1934 P1 = 0.0_dp 1935 P2 = 0.0_dp 1936 P3 = 0.0_dp 1937 P4 = 0.0_dp 1938 DO nt = 1, ntetra 1939 ibnd = jbnd 1940 ! 1941 ! etetra are the energies at the vertexes of the nt-th tetrahedron 1942 ! 1943 DO i = 1, 4 1944 etetra(i) = et(ibnd, tetra(i, nt)) 1945 ENDDO 1946 itetra(1) = 0 1947 call hpsort (4, etetra, itetra) 1948 ! 1949 ! ...sort in ascending order: e1 < e2 < e3 < e4 1950 ! 1951 e1 = etetra (1) 1952 e2 = etetra (2) 1953 e3 = etetra (3) 1954 e4 = etetra (4) 1955 ! 1956 ! kp1-kp4 are the irreducible k-points corresponding to e1- e4 1957 ! 1958 ik1 = tetra(itetra(1), nt) 1959 ik2 = tetra(itetra(2), nt) 1960 ik3 = tetra(itetra(3), nt) 1961 ik4 = tetra(itetra(4), nt) 1962 Y1 = gamma(ibnd, ik1)/et(ibnd, ik1) 1963 Y2 = gamma(ibnd, ik2)/et(ibnd, ik2) 1964 Y3 = gamma(ibnd, ik3)/et(ibnd, ik3) 1965 Y4 = gamma(ibnd, ik4)/et(ibnd, ik4) 1966 1967 IF ( e3 < ef .and. ef < e4) THEN 1968 1969 f14 = (ef- e4)/(e1- e4) 1970 f24 = (ef- e4)/(e2- e4) 1971 f34 = (ef- e4)/(e3- e4) 1972 1973 G = 3.0_dp * f14 * f24 * f34 / (e4- ef) 1974 P1 = f14 * o13 1975 P2 = f24 * o13 1976 P3 = f34 * o13 1977 P4 = (3.0_dp - f14 - f24 - f34 ) * o13 1978 1979 ELSE IF ( e2 < ef .and. ef < e3 ) THEN 1980 1981 f13 = (ef- e3)/(e1- e3) 1982 f31 = 1.0_dp - f13 1983 f14 = (ef- e4)/(e1- e4) 1984 f41 = 1.0_dp-f14 1985 f23 = (ef- e3)/(e2- e3) 1986 f32 = 1.0_dp - f23 1987 f24 = (ef- e4)/(e2- e4) 1988 f42 = 1.0_dp - f24 1989 1990 G = 3.0_dp * (f23*f31 + f32*f24) 1991 P1 = f14 * o13 + f13*f31*f23 / G 1992 P2 = f23 * o13 + f24*f24*f32 / G 1993 P3 = f32 * o13 + f31*f31*f23 / G 1994 P4 = f41 * o13 + f42*f24*f32 / G 1995 G = G / (e4- e1) 1996 1997 ELSE IF ( e1 < ef .and. ef < e2 ) THEN 1998 1999 f12 = (ef- e2)/(e1- e2) 2000 f21 = 1.0_dp - f12 2001 f13 = (ef- e3)/(e1- e3) 2002 f31 = 1.0_dp - f13 2003 f14 = (ef- e4)/(e1- e4) 2004 f41 = 1.0_dp - f14 2005 2006 G = 3.0_dp * f21 * f31 * f41 / (ef- e1) 2007 P1 = o13 * (f12 + f13 + f14) 2008 P2 = o13 * f21 2009 P3 = o13 * f31 2010 P4 = o13 * f41 2011 2012 ELSE 2013 2014 G = 0.0_dp 2015 2016 ENDIF 2017 2018 Tint = Tint + G * (Y1*P1 + Y2*P2 + Y3*P3 + Y4*P4) * vol 2019 2020 ENDDO ! ntetra 2021 2022 2023 dos_gam = Tint !2 because DOS_ee is per 1 spin 2024 2025 return 2026end function dos_gam 2027! 2028! 2029!----------------------------------------------------------------------- 2030subroutine readfg ( ifn, nr1, nr2, nr3, nat, frcg ) 2031 !----------------------------------------------------------------------- 2032 ! 2033 USE kinds, ONLY : DP 2034 USE io_global, ONLY : ionode, ionode_id, stdout 2035 USE mp, ONLY : mp_bcast 2036 USE mp_world, ONLY : world_comm 2037 implicit none 2038 ! I/O variable 2039 integer, intent(in) :: nr1, nr2, nr3, nat 2040 REAL(DP), intent(out) :: frcg(nr1, nr2, nr3, 3, 3, nat, nat) 2041 ! local variables 2042 integer i, j, na, nb, m1,m2,m3, ifn 2043 integer ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid 2044 ! 2045 ! 2046 IF (ionode) READ (ifn,*) m1, m2, m3 2047 CALL mp_bcast(m1, ionode_id, world_comm) 2048 CALL mp_bcast(m2, ionode_id, world_comm) 2049 CALL mp_bcast(m3, ionode_id, world_comm) 2050 IF ( m1 /= nr1 .or. m2 /= nr2 .or. m3 /= nr3) & 2051 call errore('readfg','inconsistent nr1, nr2, nr3 read', 1) 2052 DO i = 1, 3 2053 DO j = 1, 3 2054 DO na= 1, nat 2055 DO nb= 1, nat 2056 IF (ionode) read (ifn,*) ibid, jbid, nabid, nbbid 2057 CALL mp_bcast(ibid, ionode_id, world_comm) 2058 CALL mp_bcast(jbid, ionode_id, world_comm) 2059 CALL mp_bcast(nabid, ionode_id, world_comm) 2060 CALL mp_bcast(nbbid, ionode_id, world_comm) 2061 2062 if(i.ne.ibid.or.j.ne.jbid.or.na.ne.nabid.or.nb.ne.nbbid) THEN 2063 WRITE(stdout,*) i, j, na, nb,' <> ', ibid, jbid, nabid, nbbid 2064 call errore ('readfG','error in reading', 1) 2065 ELSE 2066 IF (ionode) read (ifn,*) (((m1bid, m2bid, m3bid, & 2067 frcg(m1,m2,m3, i, j, na, nb), & 2068 m1= 1, nr1),m2 =1, nr2),m3=1, nr3) 2069 ENDIF 2070 CALL mp_bcast(frcg(:,:,:, i, j, na, nb), ionode_id, world_comm) 2071 ENDDO 2072 ENDDO 2073 ENDDO 2074 ENDDO 2075 ! 2076 IF (ionode) CLOSE(ifn) 2077 ! 2078 return 2079end subroutine readfg 2080! 2081! 2082SUBROUTINE find_representations_mode_q ( nat, ntyp, xq, w2, u, tau, ityp, & 2083 amass, num_rap_mode, nspin_mag ) 2084 2085 USE kinds, ONLY : DP 2086 USE cell_base, ONLY : at, bg 2087 USE symm_base, ONLY : s, sr, ft, irt, nsym, nrot, t_rev, time_reversal,& 2088 sname, copy_sym, s_axis_to_cart 2089 2090 IMPLICIT NONE 2091 INTEGER, INTENT(IN) :: nat, ntyp, nspin_mag 2092 REAL(DP), INTENT(IN) :: xq(3), amass(ntyp), tau(3, nat) 2093 REAL(DP), INTENT(IN) :: w2(3*nat) 2094 INTEGER, INTENT(IN) :: ityp(nat) 2095 COMPLEX(DP), INTENT(IN) :: u(3*nat, 3*nat) 2096 INTEGER, INTENT(OUT) :: num_rap_mode(3*nat) 2097 REAL(DP) :: gi (3, 48), gimq (3), sr_is(3, 3,48), rtau(3,48, nat) 2098 INTEGER :: irotmq, nsymq, nsym_is, isym, i, ierr 2099 LOGICAL :: minus_q, search_sym, sym(48), magnetic_sym 2100! 2101! find the small group of q 2102! 2103 time_reversal=.TRUE. 2104 IF (.NOT.time_reversal) minus_q=.FALSE. 2105 2106 sym(1:nsym) =.TRUE. 2107 call smallg_q (xq, 0, at, bg, nsym, s, ft, sym, minus_q) 2108 nsymq= copy_sym(nsym, sym ) 2109 call s_axis_to_cart () 2110 CALL set_giq (xq, s, nsymq, nsym, irotmq,minus_q,gi,gimq) 2111! 2112! IF the small group of q is non symmorphic, 2113! search the symmetries only IF there are no G such that Sq -> q+G 2114! 2115 search_sym=.TRUE. 2116 IF ( ANY ( ABS(ft(:,1:nsymq)) > 1.0d-8 ) ) THEN 2117 DO isym= 1, nsymq 2118 search_sym=( search_sym.and.(ABS(gi(1, isym))<1.d-8).and. & 2119 (ABS(gi(2, isym))<1.d-8).and. & 2120 (ABS(gi(3, isym))<1.d-8) ) 2121 ENDDO 2122 ENDIF 2123! 2124! Set the representations tables of the small group of q and 2125! find the mode symmetry 2126! 2127 IF (search_sym) THEN 2128 magnetic_sym=(nspin_mag==4) 2129 CALL prepare_sym_analysis(nsymq, sr,t_rev,magnetic_sym) 2130 sym (1:nsym) = .TRUE. 2131 CALL sgam_lr (at, bg, nsym, s, irt, tau, rtau, nat) 2132 CALL find_mode_sym_new (u, w2, tau, nat, nsymq, s, sr, irt, xq, & 2133 rtau, amass, ntyp, ityp, 1, .FALSE., .FALSE., num_rap_mode, ierr) 2134 2135 ENDIF 2136 RETURN 2137 END SUBROUTINE find_representations_mode_q 2138 2139SUBROUTINE qpoint_gen1(dimx, dimy, dimz, ctrAB) 2140! 2141 use kinds, only: dp 2142 2143 IMPLICIT NONE 2144 ! input 2145 INTEGER, intent(in) :: dimx, dimy, dimz 2146 INTEGER, intent(out) :: ctrAB 2147!! REAL(DP), intent(out) :: q_AB(:,:) 2148 ! local 2149 INTEGER :: i, j, k, n, ctr, nqs 2150 REAL(DP), ALLOCATABLE :: q_all(:,:) 2151 REAL(DP) :: q_B(3), q_A(3), eps 2152 ! 2153 nqs = dimx * dimy * dimz 2154 eps = 1.0E-06 2155 ! 2156 ALLOCATE(q_all(3, nqs)) 2157 ! 2158 DO i = 1, dimx 2159 DO j = 1, dimy 2160 DO k = 1, dimz 2161 ! this is nothing but consecutive ordering 2162 n = (k - 1) + (j - 1) * dimz + (i - 1) * dimy * dimz + 1 2163 ! q_all are the components of the complete grid in crystal axis 2164 q_all(1, n) = dble(i - 1) / dimx ! + dble(k1)/2/dimx 2165 q_all(2, n) = dble(j - 1) / dimy ! + dble(k2)/2/dimy 2166 q_all(3, n) = dble(k - 1) / dimz ! + dble(k3)/2/dimz ! k1 , k2 , k3 is for the shift 2167 ENDDO 2168 ENDDO 2169 ENDDO 2170 ! 2171 ctr = 0 2172 ctrAB = 0 2173 DO i = 1, nqs 2174 q_A = q_all(:, i) + q_all(:, i) ! q_A to find if q belongs in A 2175 IF (((ABS(q_A(1)) .LT. eps) .OR. (abs(abs(q_A(1)) - 1) .LT. eps)) .AND. & 2176 ((ABS(q_A(2)) .LT. eps) .OR. (abs(abs(q_A(2)) - 1) .LT. eps)) .AND. & 2177 ((ABS(q_A(3)) .LT. eps) .OR. (abs(abs(q_A(3)) - 1) .LT. eps))) THEN 2178 ctrAB = ctrAB + 1 2179 ELSE 2180 DO j = i + 1, nqs 2181 q_B = q_all(:, i) + q_all(:, j) 2182 IF (((ABS(q_B(1)) .LT. eps) .OR. (abs(abs(q_B(1)) - 1) .LT. eps)) .AND. & 2183 ((ABS(q_B(2)) .LT. eps) .OR. (abs(abs(q_B(2)) - 1) .LT. eps)) .AND. & 2184 ((ABS(q_B(3)) .LT. eps) .OR. (abs(abs(q_B(3)) - 1) .LT. eps))) THEN 2185 ctr = ctr + 1 2186 ctrAB = ctrAB + 1 2187 END IF 2188 END DO 2189 END IF 2190 END DO 2191 ! 2192 DEALLOCATE(q_all) 2193 ! 2194 ! 2195END SUBROUTINE qpoint_gen1 2196 2197SUBROUTINE qpoint_gen2(dimx, dimy, dimz, ctrAB, q_AB) 2198! 2199 use kinds, only: dp 2200 2201 IMPLICIT NONE 2202 ! input 2203 INTEGER, intent(in) :: dimx, dimy, dimz, ctrAB 2204 REAL(DP), intent(out) :: q_AB(3, ctrAB) 2205 ! local 2206 INTEGER :: i, j, k, n, ctr, nqs 2207 REAL(DP), ALLOCATABLE :: q_all(:, :) 2208 REAL(DP) :: q_B(3), q_A(3), eps 2209 ! 2210 nqs = dimx * dimy * dimz 2211 eps = 1.0E-06 2212 ! 2213 ALLOCATE(q_all(3, nqs)) 2214 DO i = 1, dimx 2215 DO j = 1, dimy 2216 DO k = 1, dimz 2217 ! this is nothing but consecutive ordering 2218 n = (k - 1) + (j - 1) * dimz + (i - 1) * dimy * dimz + 1 2219 ! q_all are the components of the complete grid in crystal axis 2220 q_all(1, n) = dble(i - 1) / dimx ! + dble(k1)/2/dimx 2221 q_all(2, n) = dble(j - 1) / dimy ! + dble(k2)/2/dimy 2222 q_all(3, n) = dble(k - 1) / dimz ! + dble(k3)/2/dimz ! k1 , k2 , k3 is for the shift 2223 ENDDO 2224 ENDDO 2225 ENDDO 2226 ! 2227 ctr = 0 2228 DO i = 1, nqs 2229 q_A = q_all(:, i) + q_all(:, i) ! q_A to find if q belongs in A 2230 IF (((ABS(q_A(1)) .LT. eps) .OR. (abs(abs(q_A(1)) - 1) .LT. eps)) .AND. & 2231 ((ABS(q_A(2)) .LT. eps) .OR. (abs(abs(q_A(2)) - 1) .LT. eps)) .AND. & 2232 ((ABS(q_A(3)) .LT. eps) .OR. (abs(abs(q_A(3)) - 1) .LT. eps))) THEN 2233 ctr = ctr + 1 2234 q_AB(:, ctr) = q_all(:, i) 2235 ! write(*,*) "A", q_AB(:, ctr) 2236 ELSE 2237 DO j = i + 1, nqs 2238 q_B = q_all(:, i) + q_all(:, j) 2239 IF (((ABS(q_B(1)) .LT. eps) .OR. (abs(abs(q_B(1)) - 1) .LT. eps)) .AND. & 2240 ((ABS(q_B(2)) .LT. eps) .OR. (abs(abs(q_B(2)) - 1) .LT. eps)) .AND. & 2241 ((ABS(q_B(3)) .LT. eps) .OR. (abs(abs(q_B(3)) - 1) .LT. eps))) THEN 2242 ctr = ctr + 1 2243 q_AB(:, ctr) = q_all(:, i) 2244 ! write(*,*) q_AB(:, ctr) 2245 END IF 2246 END DO 2247 END IF 2248 END DO 2249 ! 2250 DEALLOCATE(q_all) 2251 2252END SUBROUTINE qpoint_gen2 2253 2254SUBROUTINE ZG_configuration(nq, nat, ntyp, amass, ityp, q, w2, z_nq_zg, ios, & 2255 dimx, dimy, dimz, nloops, error_thresh, synch, tau, alat, atm, & 2256 ntypx, at, q_in_cryst_coord, q_external, T, incl_qA, & 2257 compute_error, single_phonon_displ) 2258! 2259 use kinds, only: dp 2260 use constants, only: amu_ry, ry_to_thz, ry_to_cmm1, H_PLANCK_SI, & 2261 K_BOLTZMANN_SI, AMU_SI, pi 2262 USE cell_base, ONLY : bg 2263 USE io_global, ONLY : ionode 2264 IMPLICIT NONE 2265 ! input 2266 CHARACTER(LEN=3), intent(in) :: atm(ntypx) 2267 LOGICAL, intent(in) :: synch, q_in_cryst_coord, q_external 2268 LOGICAL, intent(in) :: incl_qA, compute_error, single_phonon_displ 2269 INTEGER, intent(in) :: dimx, dimy, dimz, nloops 2270 INTEGER, intent(in) :: nq, nat, ntyp, ios, ntypx 2271 INTEGER, intent(in) :: ityp(nat) 2272 ! nq is the number of qpoints in sets A and B 2273 REAL(DP), intent(in) :: error_thresh, alat, T 2274 REAL(DP), intent(in) :: at(3, 3) 2275 REAL(DP), intent(in) :: q(3, nq), w2(3 * nat, nq), amass(ntyp), tau(3, nat) 2276 COMPLEX(DP), intent(in) :: z_nq_zg(3 * nat, 3 * nat, nq) 2277 ! 2278 ! local 2279 CHARACTER(len=256) :: filename 2280 CHARACTER(len=256) :: pointer_etta 2281 ! 2282 INTEGER :: nat3, na, nta, ipol, i, j, k, qp, ii, p, kk 2283 INTEGER :: nq_tot, pn, combs, combs_all, sum_zg 2284 INTEGER, ALLOCATABLE :: Mx_mat(:, :), Mx_mat_or(:, :), M_mat(:, :), V_mat(:) 2285 INTEGER, ALLOCATABLE :: Rlist(:, :) 2286 ! nq_tot total number of q-points (including sets A, B, C) 2287 ! pn combinations 2288 INTEGER :: ctr, ctr2, ctrA, ctrB, ctrAB 2289 ! 2290 REAL(DP) :: freq(3 * nat, nq), ph_w(3 * nat, nq), l_q(3 * nat, nq) 2291 REAL(DP) :: q_A(3), q_B(3), p_q(3 * nat, nq), e_nq(3 * nat, nq) 2292 ! l_q is the amplitude \sigma at temperature T, 2293 ! e_nq --> to calculate total vibrational energy 2294 ! PE_nq --> Potential enrgy: 1/2 Mp \omega_\nu^2 x_\nu^2 2295 ! p_q is the momentum on the nuclei \hbar\2\l_\bq\nu \SQRT(n_{q\nu, T}+ 1/2) 2296 ! 2297 ! 2298 REAL(DP), PARAMETER :: eps = 1.0d-6 2299 REAL(DP) :: hbar, ang, JOULE_TO_RY, u_rand, dotp, PE_nq, KE_nq 2300 ! ALLOCATE TABLES 2301 REAL(DP), ALLOCATABLE :: equil_p(:, :, :), T_fact(:, :), DW_fact(:, :), qA(:, :), qB(:, :), DWp_fact(:, :), Tp_fact(:, :) 2302 ! for displacements 2303 REAL(DP), ALLOCATABLE :: Cx_matA(:, :), Cx_matB(:, :), Cx_matAB(:, :), Bx_vect(:) 2304 ! for momenta/velocities 2305 REAL(DP), ALLOCATABLE :: Cpx_matA(:,:), Cpx_matB(:,:), Cpx_matAB(:,:) 2306 ! matrices to account for the coupling terms between different phonon branches ! 2307 REAL(DP), ALLOCATABLE :: sum_error_D(:, :), sum_diag_D(:, :), sum_error_B(:), sum_diag_B(:), sum_error_B2(:), sum_diag_B2(:) 2308 REAL(DP), ALLOCATABLE :: D_tau(:, :, :), P_tau(:, :, :), ratio_zg(:)! displacements and velocities 2309 REAL(DP), ALLOCATABLE :: R_mat(:, :), E_vect(:, :), D_vect(:, :), F_vect(:, :) 2310 ! D_tau : atomic displacements 2311 ! z_nq_A : eigenvectors for q-points in set A 2312 ! z_nq_B : eigenvectors for q-points in set B 2313 ! M matrices : sign matrices 2314 ! R_mat, E_vect, D_vect, F_vect : are used to compute the minimization of the 2315 ! error coming from the off diagonal terms --> sum_error_D ! 2316 ! sum_diag_D : the sum of diagonal terms contributing to the T-dependent properties 2317 ! 2318 COMPLEX(DP) :: z_zg(3 * nat, 3 * nat, nq) 2319 COMPLEX(DP) :: imagi 2320 COMPLEX(DP), ALLOCATABLE :: z_nq_synch(:, :, :), z_nq_A(:, :, :), z_nq_B(:, :, :) 2321 ! singular value decomposition matrices U = R*conj(L) 2322 ! 2323 INTEGER :: INFO, N_dim, M_dim, K_dim, L_dim , LWORK 2324 REAL(DP), ALLOCATABLE :: RWORK(:), S_svd(:) 2325 COMPLEX(DP), ALLOCATABLE :: M_over(:, :, :), U_svd(:, :, :), U_svd_d(:, :), dotp_mat(:, :) 2326 COMPLEX(DP), ALLOCATABLE :: L_svd(:, :), R_svd(:, :), WORK(:), U_svd_d_new(:, :) 2327 COMPLEX*16 dum( 1 ) ! for the ZGEEV 2328 ! 2329 ! 2330 ! 2331 ! constants to be used 2332 hbar = 0.5 * H_PLANCK_SI / pi ! reduce Plnack constant 2333 JOULE_TO_RY = 2.1798741E-18 ! joule to rydberg 2.1798741E-18 2334 ang = 1.0E-10 ! angstrom units 2335 imagi = (0.0d0, 1.0d0) !imaginary unit 2336 ! Inititialize eigenvectors matrix 2337 z_zg = (0.d0, 0.d0) 2338 ! Set intitial values 2339 nq_tot = dimx * dimy * dimz 2340 nat3 = 3 * nat 2341 pn = 2**(nat3 - 1) 2342 ! it's pointless to ALLOCATE more signs a very large number of branches 2343 IF ( nat3 > 12) pn = 2**(12 - 1) 2344 ! 2345 ! SVD parameters 2346 M_dim = nat3 2347 N_dim = nat3 2348 K_dim = MIN(M_dim, N_dim) 2349 L_dim = MAX(M_dim, N_dim) 2350 ! 2351 ! create equilibrium configuration 2352 ALLOCATE(equil_p(nq_tot, nat, 3)) 2353 call create_supercell(at, tau, alat, dimx, dimy, dimz, nat, equil_p) 2354 filename = 'equil_pos.txt' 2355 OPEN (unit = 70, file = filename, status = 'unknown', form = 'formatted') 2356 WRITE(70,*) "Number of atoms", nat * dimx * dimy * dimz 2357 WRITE(70,*) 'equilibrium positions, (Ang):' 2358 DO i = 1, nq_tot 2359 DO k = 1, nat 2360 WRITE(70,'(A6, 3F13.8)') atm(ityp(k)), equil_p(i, k, :) 2361 ENDDO 2362 ENDDO 2363 CLOSE(70) 2364 ! 2365 ! convert eigenvectors to mass-unscalled 2366 DO i = 1, nat3 2367 DO na = 1, nat 2368 nta = ityp(na) 2369 DO ipol = 1, 3 2370 DO qp = 1, nq 2371 z_zg((na - 1) * 3 + ipol, i, qp) = z_nq_zg((na - 1) * 3 + ipol, i, qp) * SQRT(amu_ry*amass(nta)) 2372 ENDDO 2373 ENDDO 2374 ENDDO 2375 ENDDO 2376! 2377! 2378! Frequency check 2379 WRITE(*,*) 2380 freq = 0.0d0 2381 DO qp = 1, nq 2382 DO i = 1, nat3 2383 ! IF (w2(i, qp) .lt. 1.0E-8) THEN 2384 IF (w2(i, qp) .lt. 0.0d0) THEN 2385 WRITE(*,*) "WARNING: Negative ph. freqs:", w2(i, qp), i, qp 2386 freq(i, qp) = SQRT(ABS(w2(i, qp))) 2387 ELSE 2388 freq(i, qp) = SQRT(ABS(w2(i, qp))) 2389 ENDIF 2390 ENDDO 2391 ENDDO 2392 ! 2393 ph_w = freq * ry_to_thz * (1.0E12) * 2 * pi ! correct frequency for phonons in SI 2394 ! 2395 ! 2396 ! set amplitudes of displacements l_q = \sigma_\bq\nu and momenta 2397 p_q = 0.0d0 2398 DO qp = 1, nq 2399 DO i = 1, nat3 2400 l_q(i, qp) = SQRT(hbar / ph_w(i, qp) / 2.0d0 / AMU_SI / ang**2.0d0) * & 2401 SQRT(dble(1.0d0 + 2.0d0 / (EXP(hbar * ph_w(i, qp) / (K_BOLTZMANN_SI * T)) - 1))) 2402 p_q(i, qp) = hbar / SQRT(2.0d0) / (SQRT(hbar / ph_w(i, qp) / 2.0d0 /AMU_SI)) / ang * & !*1.0E-12& 2403 SQRT(dble(0.5d0 + 1.0d0 / (EXP(hbar * ph_w(i, qp) / (K_BOLTZMANN_SI * T)) - 1))) 2404 ! we can multiply by 1.0E-12 to get 'picos' 2405 !! WRITE(*,*) p_q(i, qp),l_q(i, qp), amass(1) 2406 ENDDO 2407 ENDDO 2408 ! 2409! WRITE(*,*) "total vibrational energy per cell", 2*dotp/dimx/dimy/dimz, "Ry" 2410 IF (q_external) THEN 2411 IF (q_in_cryst_coord .EQV. .FALSE.) THEN 2412 ! in both cases convert them to crystal 2413 CALL cryst_to_cart(nq, q, at, -1) 2414 ELSE 2415 CALL cryst_to_cart(nq, q, at, -1) 2416 ENDIF 2417 ELSE 2418 CALL cryst_to_cart(nq, q, at, -1) 2419 ENDIF 2420 ! To distinguish between different sets of qpoints, A, B, C 2421 ! to find how many points belong to set A and then allocate matrix accordingly 2422 ! NOTE that we want the qpoints always in crystal coordinates 2423 ! 2424 ctrA = 0 2425 ctrAB = 0 2426 dotp = 0.0d0 2427 PE_nq = 0.0d0 2428 KE_nq = 0.0d0 2429 DO qp = 1, nq 2430 q_A = q(:, qp) + q(:, qp) ! q_A to find IF q belongs in A 2431 IF (((ABS(q_A(1)) < eps) .OR. (ABS(ABS(q_A(1)) - 1) < eps)) .AND. & 2432 ((ABS(q_A(2)) < eps) .OR. (ABS(ABS(q_A(2)) - 1) < eps)) .AND. & 2433 ((ABS(q_A(3)) < eps) .OR. (ABS(ABS(q_A(3)) - 1) < eps))) THEN 2434 ! WRITE(*,*) "set A", qp, q(:, qp) 2435 ctrA = ctrA + 1 2436 ctrAB = ctrAB + 1 2437 DO i = 1, nat3 2438 e_nq(i, qp) = hbar * ph_w(i, qp) * (0.5d0 + 1.0d0 / (EXP(hbar*ph_w(i, qp) / (K_BOLTZMANN_SI * T)) - 1)) 2439 dotp = dotp + e_nq(i, qp) / JOULE_TO_RY ! joule to rydberg 2.1798741E-18 2440 PE_nq = PE_nq + 0.5d0 * AMU_SI * ph_w(i, qp)**2 * l_q(i, qp)**2 * ang**2.0d0 / JOULE_TO_RY 2441 KE_nq = KE_nq + 0.5d0 / AMU_SI * p_q(i, qp)**2 * ang**2 / JOULE_TO_RY 2442 ENDDO 2443 ELSE 2444 ctrAB = ctrAB + 1 2445 DO i = 1, nat3 2446 e_nq(i, qp) = hbar * ph_w(i, qp) * (0.5d0 + 1.0d0 / (EXP(hbar*ph_w(i, qp) / (K_BOLTZMANN_SI * T)) - 1)) 2447 ! Factor of two because I account half of the q-points (set B equiv. to set C) 2448 dotp = dotp + 2.0d0 * e_nq(i, qp) / JOULE_TO_RY ! joule to rydberg 2.1798741E-18 2449 PE_nq = PE_nq + 2.0d0 * 0.5d0 * AMU_SI * ph_w(i, qp)**2 * l_q(i, qp)**2 * ang**2.0d0 / JOULE_TO_RY 2450 KE_nq = KE_nq + 2.0d0 * 0.5d0 / AMU_SI * p_q(i, qp)**2 * ang**2 / JOULE_TO_RY 2451 ENDDO 2452 ENDIF 2453 ENDDO 2454 ! 2455 WRITE(*,*) 2456 WRITE(*,'(A30, 1F13.8, A4)') "Total vibrational energy: ", dotp, "Ry" 2457 WRITE(*,*) 2458 WRITE(*,'(A30, 1F13.8, A4)') "Potential energy: ", PE_nq, "Ry" 2459 WRITE(*,*) 2460 WRITE(*,'(A30, 1F13.8, A4)') "Kinetic energy: ", KE_nq, "Ry" 2461 WRITE(*,*) 2462 WRITE(*,*) "Note that the total energy output from a DFT-ZG & 2463 calculation accounts for half the total vibrational energy" 2464 WRITE(*,*) 2465 ! 2466 ctrB = ctrAB - ctrA 2467 WRITE(*,*) "Points in sets AB, A, B :", ctrAB, ctrA, ctrB 2468 WRITE(*,*) 2469 ! 2470 ALLOCATE(qA(ctrA, 3), qB(ctrB, 3), z_nq_A(nat3, nat3, ctrA), z_nq_B(nat3, nat3, ctrB)) 2471 ALLOCATE(Cx_matAB(nat3, ctrAB), Cx_matA(nat3, ctrA), Cx_matB(nat3, ctrB)) 2472 ALLOCATE(Cpx_matAB(nat3, ctrAB), Cpx_matA(nat3, ctrA), Cpx_matB(nat3, ctrB)) 2473 ALLOCATE(D_tau(nq_tot, nat, 3), P_tau(nq_tot, nat, 3), Rlist(nq_tot, 3)) 2474 ALLOCATE(T_fact(nat, 3), DW_fact(nat, 3), DWp_fact(nat, 3), Tp_fact(nat, 3)) 2475 Cx_matAB = 0 2476 Cpx_matAB = 0 2477 ! 2478 ! Generate lattice vectors in crystal coordinates 2479 ctr2 = 1 2480 DO i = 0, dimz - 1 2481 DO j = 0, dimy - 1 2482 DO k = 0, dimx - 1 2483 Rlist(ctr2, 1) = k 2484 Rlist(ctr2, 2) = j 2485 Rlist(ctr2, 3) = i !(/ k, j, i /) 2486 !WRITE(*,*) Rlist(ctr2,:) 2487 ctr2 = ctr2 + 1 2488 ENDDO 2489 ENDDO 2490 ENDDO 2491 ! 2492 ! for accoustic modes put l_q\nu = 0and p_q\nu = 0 so we freeze them 2493 ! 2494 DO qp = 1, ctrAB 2495 q_A = q(:, qp) ! q_A to find IF q belongs in A 2496 IF (((ABS(q_A(1)) < eps)) .AND. ((ABS(q_A(2)) < eps)) .AND. & 2497 ((ABS(q_A(3)) < eps))) THEN 2498 ! WRITE(*,*) "accoustic modes at Gamma", q_A 2499 l_q(1, qp) = 0.0d0 2500 l_q(2, qp) = 0.0d0 2501 l_q(3, qp) = 0.0d0 2502 p_q(1, qp) = 0.0d0 2503 p_q(2, qp) = 0.0d0 2504 p_q(3, qp) = 0.0d0 2505 ENDIF 2506 ENDDO 2507 ! 2508 ! 2509 ! First step of synchronization to make all eigenvectors "positive" based on their first entry. 2510 !DO i = 1, nq 2511 ! DO j = 1, nat3 2512 ! IF (REAL(z_zg(1, j, i)) < 0.0d0) THEN 2513 ! z_zg(:, j, i) = -z_zg(:, j, i) 2514 ! ENDIF 2515 ! ENDDO 2516 !ENDDO 2517 ! Now main synch proc as in paper. Do this procedure for every pn. 2518 ! 2519 IF (synch) THEN 2520 ! 2521 ALLOCATE(M_over(nat3, nat3, pn - 1), U_svd(nat3, nat3, pn - 1), z_nq_synch(nat3, nat3, ctrAB)) 2522 ALLOCATE(U_svd_d(nat3, pn - 1), dotp_mat(nat3, nat3), U_svd_d_new(nat3, pn-1)) 2523 ALLOCATE(L_svd(M_dim,K_dim), R_svd(K_dim,N_dim),S_svd(K_dim)) 2524 z_nq_synch = (0.0d0 , 0.0d0) 2525 ! query workspace 2526 ! 2527 LWORK = 5 * nat3 !MAX(1, 2*K_dim+L_dim) 2528 ! 2529 ALLOCATE( RWORK( LWORK ) ) 2530 ALLOCATE( WORK( LWORK ) ) 2531 LWORK = - 1 2532 2533 call ZGESVD('A','A', nat3, nat3, M_over(:,:, 1), nat3, S_svd, L_svd, & 2534 nat3, R_svd, nat3, WORK, LWORK, RWORK, INFO) 2535 ! 2536 LWORK = INT(WORK(1)) + 1 2537 ! 2538 IF( LWORK > SIZE( WORK ) ) THEN 2539 DEALLOCATE( WORK ) 2540 ALLOCATE( WORK( LWORK ) ) 2541 ENDIF 2542 ! 2543 ! 2544 DO i = 0, ctrAB - pn, pn 2545 z_nq_synch(:,:, i + 1) = z_zg(:,:, i + 1) 2546 ! z_nq_synch(:,:, ctrAB-pn -i + 1) = z_zg(:,:, i + 1) 2547 DO ii = 1, pn - 1 2548 M_over = 0.d0 2549 ! Construct the overlap matrix M_{\nu,\nu'} 2550 S_svd = 0.d0 2551 DO p = 1, nat3 2552 DO j = 1, nat3 2553 DO k = 1, nat3 ! sum over \k,\a 2554 M_over(j, p, ii) = M_over(j, p, ii) + (z_zg(k, j, ii + i + 1) * CONJG(z_nq_synch(k, p, ii + i))) 2555 ENDDO ! k-loop 2556 ENDDO ! j-loop 2557 ENDDO ! p-loop 2558 ! perform singular value decomposition 2559 call ZGESVD('A', 'A', nat3, nat3, M_over(:,:, ii), nat3, S_svd, L_svd, & 2560 nat3, R_svd, nat3, WORK, LWORK, RWORK,INFO) 2561 U_svd(:,:, ii) = MATMUL(TRANSPOSE(CONJG(R_svd)),TRANSPOSE(CONJG(L_svd))) 2562 call ZGEEV( 'N', 'N', nat3, U_svd(:,:, ii), nat3, U_svd_d(:, ii), dum, 1, dum, 1, & 2563 WORK, LWORK, RWORK, INFO ) 2564 M_over = 0.0d0 2565 DO p = 1, nat3 2566 DO j = 1, nat3 2567 DO k = 1, nat3 ! sum over \k,\a 2568 M_over(j, p, ii) = M_over(j, p, ii) + (z_zg(k, j, ii + i + 1) * CONJG(z_nq_synch(k, p, ii + i))) 2569 ENDDO ! k-loop 2570 ENDDO ! j-loop 2571 ENDDO ! p-loop 2572 DO qp = 1, nat3 2573 DO k = 1, nat3 2574 dotp_mat(qp, k) = CONJG(M_over(qp, qp, ii)) * CONJG(U_svd_d(k, ii)) 2575 ENDDO 2576 ENDDO 2577 dotp_mat = ABS(REAL(dotp_mat)) 2578 DO qp = 1, nat3 2579 p = MAXLOC(REAL(dotp_mat(qp,:)), 1) 2580 U_svd_d_new(qp, ii) = U_svd_d(p, ii) 2581 ENDDO 2582 DO qp = 1, nat3 2583 DO k = 1, nat3 2584 z_nq_synch(k, qp, ii + i + 1) = U_svd_d_new(qp, ii) * z_zg(k, qp, ii + i + 1) 2585 ENDDO 2586 ENDDO 2587 z_zg(:, :, ii + i + 1) = z_nq_synch(:, :, ii + i + 1) 2588 ENDDO ! ii-loop 2589 ENDDO ! i-loop 2590 ! Here we synchronize the remaining eigenvectors IF ctrAB is not divided by pn 2591 IF (mod(ctrAB, pn) > 0) THEN 2592 ctr = ctrAB - mod(ctrAB, pn) 2593 z_nq_synch(:, :, ctr + 1) = z_zg(:, :, ctr + 1) 2594 DO ii = 1, mod(ctrAB, pn) - 1 2595 M_over = 0.0d0 2596 ! Construct the overlap matrix M_{\nu,\nu'} 2597 S_svd = 0.0d0 2598 DO p = 1, nat3 2599 DO j = 1, nat3 2600 DO k = 1, nat3 2601 M_over(j, p, ii) = M_over(j, p, ii) + (z_zg(k, j, ii + ctr + 1) * CONJG(z_nq_synch(k, p, ii + ctr))) 2602 ENDDO ! k-loop 2603 ENDDO ! j-loop 2604 ENDDO ! p-loop 2605 ! perform singular value decomposition 2606 call ZGESVD('S','S', nat3, nat3, M_over(:,:, ii), nat3, S_svd, L_svd, & 2607 nat3, R_svd, nat3, WORK, LWORK, RWORK,INFO) 2608 ! ZGESVD returns R_svd**H (the hermitian TRANSPOSE of R_svd) 2609 U_svd(:,:, ii) = MATMUL(TRANSPOSE(CONJG(R_svd)),TRANSPOSE(CONJG(L_svd))) 2610 call ZGEEV( 'N', 'N', nat3, U_svd(:, :, ii), nat3, U_svd_d(:, ii), dum, 1, dum, 1, & 2611 WORK, LWORK, RWORK, INFO ) 2612 M_over = 0.0d0 2613 DO p = 1, nat3 2614 DO j = 1, nat3 2615 DO k = 1, nat3 ! sum over \k,\a 2616 M_over(j, p, ii) = M_over(j, p, ii) + (z_zg(k, j, ii + ctr + 1) * CONJG(z_nq_synch(k, p, ii + ctr))) 2617 ENDDO ! k-loop 2618 ENDDO ! j-loop 2619 ENDDO ! p-loop 2620 DO qp = 1, nat3 2621 DO k = 1, nat3 2622 dotp_mat(qp, k) = CONJG(M_over(qp, qp, ii)) * CONJG(U_svd_d(k, ii)) 2623 ENDDO 2624 ENDDO 2625 dotp_mat = ABS(DBLE(dotp_mat)) 2626 DO qp = 1, nat3 2627 p = MAXLOC(DBLE(dotp_mat(qp, :)), 1) 2628 U_svd_d_new(qp, ii) = U_svd_d(p, ii) 2629 ENDDO 2630 ! 2631 DO qp = 1, nat3 2632 DO k = 1, nat3 2633 z_nq_synch(k, qp, ii + ctr + 1) = U_svd_d_new(qp, ii) * z_zg(k, qp, ii + ctr + 1) 2634 ENDDO 2635 ENDDO 2636 ! overwrite z_zg 2637 z_zg(:, :, ii + ctr + 1) = z_nq_synch(:, :, ii + ctr + 1) 2638 ENDDO ! ii-loop 2639 ENDIF ! mod(ctrAB, pn) 2640 DEALLOCATE(WORK, RWORK) 2641 ! 2642 ! DEALLOCATE matrices for synch proc. 2643 DEALLOCATE(M_over, U_svd, z_nq_synch, U_svd_d, U_svd_d_new, dotp_mat) 2644 DEALLOCATE(L_svd, R_svd,S_svd) 2645 ! 2646 ENDIF ! end IF synch is TRUE 2647 ! 2648 ! Initialize sign matrices: 2649 ! sign matrices Mx and MY . Total entries pf Mx should be 2^nmodes/2! I divide by two so we get only independent entries. 2650 ! i.e. [1 1 1 1 1 1] gives same result to [- 1 -1 -1 -1 -1 -1]. 2651 ! 2652 ! 2653 ALLOCATE(Mx_mat_or(pn, nat3)) 2654 ALLOCATE(Mx_mat(pn, nat3), M_mat(2 * pn, nat3), V_mat(2)) 2655 V_mat = (/ 1, -1/) ! initialize V_mat whose entries will generate the sign matrices 2656 DO i = 1, nat3 2657 ctr = 1 2658 DO p = 1, 2**(i - 1) 2659 DO qp = 1, 2 2660 DO k = 1, 2**(nat3 - i) 2661 IF (ctr > 2 * pn) EXIT ! I DO this in case there many branches in the system and 2662 ! in that case we DO not need to ALLOCATE more signs 2663 M_mat(ctr, i) = V_mat(qp) 2664 ctr = ctr + 1 2665 IF (ctr > 2 * pn) EXIT 2666 ENDDO 2667 ENDDO 2668 ENDDO 2669 ENDDO 2670 ! NOTE: In M_mat the first half entries are the independent set of signs (2** (nat3- 1)) ! 2671 ! The rest entries are the antithetics !! 2672 ! checks 2673 WRITE(*,*) "Sign matrix" 2674 WRITE(*,*) "-----------" 2675 DO j = 1, 2 * pn !2** (nat3) 2676 IF (MOD(j, 2) == 0) M_mat(j,:) = -1 * M_mat(j,:) 2677 WRITE(*,'(100i3)') M_mat(j,:) 2678 ENDDO 2679 ! checks_done 2680 ! 2681 combs = 0! how many unique pairs x1x2, x1x3,... we have. Note without diagonal terms x1^2, x2^2, ... ! So we define R matrix 2682 DO i = 1, nat3 - 1 2683 combs = combs + i 2684 ENDDO 2685 combs_all = 2 * combs + nat3 !; % with x1^2, x2^2 ... 2686 ! 2687 ! combs_all refere also to all possible pais ({\k,\a}, {\k' \a'}) 2688 ! 2689 ALLOCATE(ratio_zg(combs_all)) 2690 ratio_zg = 0.d0 2691 IF (compute_error) THEN 2692 ALLOCATE(sum_error_D(combs_all, INT(ctrAB / pn) + 1), sum_error_B(combs_all)) 2693 ALLOCATE(sum_error_B2(combs_all * NINT(nq_tot / 2.0d0)), sum_diag_B2(combs_all * NINT(nq_tot / 2.0d0))) 2694 ! I add one because of the reminder when ctrAB is not divided by pn 2695 ALLOCATE(sum_diag_D(combs_all, INT(ctrAB / pn) + 1), sum_diag_B(combs_all)) 2696 ALLOCATE(R_mat(combs_all, combs_all), D_vect(combs_all, ctrAB), F_vect(combs_all, ctrAB)) 2697 ALLOCATE(Bx_vect(combs_all), E_vect(combs_all, ctrAB)) 2698 ENDIF 2699 ! 2700 ! Instead of taking all possible permutations which are pn=2** (nmodes- 1)! 2701 ! we just select possible permutations until the error is lower than a 2702 ! threshold. The lower the threshold the longer the algorithm can take. 2703 ! filename = 'ZG-configuration.txt' 2704 WRITE(pointer_etta,'(f5.3)') error_thresh 2705 filename = 'ZG-configuration_' // TRIM( pointer_etta ) // '.dat' !'.fp' 2706 OPEN (unit = 80, file = filename, status = 'unknown', form = 'formatted') 2707 filename = 'ZG-velocities_' // TRIM( pointer_etta ) // '.dat' !'.fp' 2708 OPEN (unit = 81, file = filename, status = 'unknown', form = 'formatted') 2709 DO kk = 1, nloops 2710 ! Allocate original matrices ! half the entries of M_mat 2711 ! We also make the inherent choice that each column of Mx_mat_or has the same number of positive and negative signs 2712 Mx_mat_or = 0 2713 DO i = 1, 2 * pn / 4, 2 2714 Mx_mat_or(i, :) = M_mat(i, :) 2715 ENDDO 2716 ctr = 1 2717 DO i = 2, 2 * pn / 4, 2 2718 Mx_mat_or(i, :) = M_mat(2 * pn + 1 - i, :) 2719 ctr = ctr + 1 2720 ENDDO 2721 DO i = 2 * pn / 4 + 1, 2 * pn / 2, 2 2722 Mx_mat_or(i, :) = M_mat(i + 1, :) 2723 ctr = ctr + 1 2724 ENDDO 2725 DO i = 2 * pn / 4 + 2, 2 * pn / 2, 2 2726 Mx_mat_or(i, :) = M_mat(2 * pn + 2 - i, :) 2727 ctr = ctr + 1 2728 ENDDO 2729 ! 2730 ! 2731 ! 2732 DO i = pn, 1, -1 2733 ! To generate integer numbers from 1 to pn 2734 ! and ALLOCATE the M_x,y matrices ! so we DO not have specific order ! 2735 call random_number(u_rand) 2736 ii = 1 + FLOOR(i * u_rand) 2737 Mx_mat(i, :) = Mx_mat_or(ii, :) 2738 Mx_mat_or(ii, :) = Mx_mat_or(i, :) ! so I DO not repeat this entry 2739 ENDDO 2740 ! DO q-points in sets A n B 2741 ! based on the sets of signs we genereated from the above loop 2742 DO ii = 1, ctrAB ! loop over qpoints 2743 ! change the signs of Mx in every 2^(nmodes-2) entries 2744 ! I use Mx_mat_or to apply concatenate matrices and obtain set E 2745 IF (mod(ii, pn) .EQ. 1) THEN 2746 Mx_mat_or(1 : pn - 1, :) = Mx_mat(2 : pn, :) ! second element goes to the top 2747 Mx_mat_or(pn, :) = Mx_mat(1, :) ! first element goes to the bottom 2748 ENDIF 2749 Mx_mat = Mx_mat_or 2750 ! to take antithetics every pn 2751 ! Remember always the error goes to very small as long we have equal number of + and - signs 2752 ! and displacements remain around equilibrium ! 2753 !IF (mod(ii, pn).EQ.1) THEN 2754 ! Mx_mat=-Mx_mat 2755 !ENDIF 2756 ! 2757 ! 2758 ! Cx_matAB contains all the sigmas with the appropriates signs 2759 IF (MOD(ii, pn) > 0) THEN 2760 DO k = 1, nat3 2761 Cx_matAB(k, ii) = l_q(k, ii) * Mx_mat(mod(ii, pn), k) ! mod so values of every pn q-points are repeated 2762 Cpx_matAB(k, ii) = p_q(k, ii) * Mx_mat(mod(ii, pn), k) ! mod so values of every pn q-points are repeated 2763 ENDDO 2764 ELSE 2765 DO k = 1, nat3 2766 Cx_matAB(k, ii) = l_q(k, ii) * Mx_mat(pn, k) 2767 Cpx_matAB(k, ii) = p_q(k, ii) * Mx_mat(pn, k) 2768 ENDDO 2769 ENDIF 2770 ! create R matrix that contains Re[e_ka^nu(q) * e_k'a'^nu'* (q)] 2771 IF (compute_error) THEN 2772 R_mat = 0.0d0 2773 D_vect = 0.0d0 2774 ctr = 1 2775 DO p = 1, nat3 2776 DO qp = 1, nat3 !those are for \k \a, \k' \a' 2777 ctr2 = 1 2778 !------------------------------ 2779 DO i = 1, nat3 2780 DO j = 1, nat3 2781 D_vect(ctr, ii) = D_vect(ctr, ii) + DBLE(z_zg(p, i, ii) * CONJG(z_zg(qp, j, ii))) * & 2782 Cx_matAB(i, ii) * Cx_matAB(j, ii) 2783 R_mat(ctr, ctr2) = DBLE(z_zg(p, i, ii) * CONJG(z_zg(qp, j, ii))) 2784 ctr2 = ctr2 + 1 2785 ENDDO 2786 ENDDO 2787 !---------------------------- 2788 ctr = ctr + 1 2789 ENDDO 2790 ENDDO ! end p loop ! R_mat is filled 2791 ! 2792 Bx_vect = 0.d0 2793 ! Bx_vect will contain all the cross terms v .neq. v' and diagonal for each q 2794 ctr = 1 2795 DO i = 1, nat3 2796 DO j = 1, nat3 2797 Bx_vect(ctr) = Cx_matAB(i, ii) * Cx_matAB(j, ii) 2798 ctr = ctr + 1 2799 ENDDO 2800 ENDDO 2801 ! 2802 E_vect(:, ii) = 0.d0 2803 ! E_vect contains only the diagonal terms (i.e. v = v') 2804 ctr = 1 2805 DO p = 1, nat3 ! those are for \k \a, \k' \a' 2806 DO qp = 1, nat3 2807 DO i = 1, nat3 2808 DO j = 1, nat3 2809 IF (j == i) THEN 2810 E_vect(ctr, ii) = E_vect(ctr, ii) + DBLE(z_zg(p, i, ii) * CONJG(z_zg(qp, j, ii))) * Cx_matAB(i, ii)**2 2811 ENDIF 2812 ENDDO 2813 ENDDO 2814 ctr = ctr + 1 2815 ENDDO 2816 ENDDO ! p loop 2817 ! D_vect(:, ii) = MATMUL(R_mat, Bx_vect) 2818 ! D_vect contains the diagonal and the non-diagonal terms (i.e. v = v' and v .neq. v') 2819 ! E_vect contains only the diagonal terms (i.e. v = v') 2820 ! F_vect contains the error (i.e.each entry is the contribution from v \neq v') at each q point to minimize 2821 ! checks 2822 F_vect(:, ii) = D_vect(:, ii) - E_vect(:, ii) 2823 ! 2824 ENDIF ! compute_error 2825 ENDDO ! ii loop over qpoints 2826 ! 2827 ! 2828 !Compute error 2829 ! 2830 ! 2831 IF (compute_error) THEN 2832 sum_error_D = 0.0d0 2833 sum_error_B = 0.0d0 2834 ! sum_error_D : contains the error from \nu and \nu' every pn 2835 !!!!!!! 2836 WRITE(*,*) 2837 WRITE(*,*) "Minimize error based on threshold" 2838 DO p = 1, combs_all 2839 ctr = 1 2840 DO i = 0, INT(ctrAB / pn) - 1 2841 sum_error_D(p, ctr) =SUM(F_vect(p, pn * i + 1 : pn * (i + 1))) ! pn) is the length of Mx 2842 ctr = ctr + 1 2843 ENDDO 2844 ! Here we add the reminder IF ctrAB is not divided exactly by pn 2845 IF (mod(ctrAB, pn) > 0) THEN 2846 sum_error_D(p, ctr) = SUM(F_vect(p, ctrAB - mod(ctrAB, pn) + 1 : ctrAB)) ! add the remaining terms 2847 ENDIF 2848 ! evaluate also error from all q-points in B 2849 DO i = 1, ctrAB 2850 sum_error_B(p) = sum_error_B(p) + F_vect(p, i) ! 2851 ENDDO 2852 ! 2853 ENDDO ! end p-loop 2854 ! 2855 sum_error_B2 = 0.0d0 2856 ctr2 = 1 2857 DO j = 1, NINT(nq_tot / 2.0d0) 2858 DO p = 1, combs_all ! for ever \k,\a,\k',\a' 2859 DO i = 1, ctrAB ! over qpoints 2860 dotp = 0.0d0 2861 DO ii = 1, 3 2862 dotp = dotp + q(i, ii) * Rlist(j, ii)! 2863 ENDDO ! ii 2864 sum_error_B2(ctr2) = sum_error_B2(ctr2) + cos(2.0d0 * pi * dotp) * F_vect(p, i) 2865 ! 2866 ENDDO ! i 2867 ctr2 = ctr2 + 1 2868 ENDDO ! p 2869 ENDDO ! j 2870 ! 2871 sum_diag_D = 0.0d0 2872 sum_diag_B = 0.0d0 2873 ! sum_diag_D : contains the diagonal terms 2874 DO p = 1, combs_all 2875 ctr = 1 2876 DO i = 0, INT(ctrAB / pn) - 1 2877 sum_diag_D(p, ctr) =SUM(E_vect(p, pn * i + 1 : pn * (i + 1))) ! pn) is the length of Mx 2878 ctr = ctr + 1 2879 ENDDO 2880 ! Here we add the reminder IF ctrAB is not divided exactly by pn 2881 IF (mod(ctrAB, pn) > 0) THEN 2882 sum_diag_D(p, ctr) = SUM(E_vect(p, ctrAB - mod(ctrAB, pn) + 1 : ctrAB)) ! add the remaining terms 2883 ENDIF 2884 DO i = 1, ctrAB 2885 sum_diag_B(p) = sum_diag_B(p) + E_vect(p, i) ! pn) is the length of Mx 2886 ENDDO 2887 ENDDO ! end p-loop 2888 sum_diag_B2 = 0.0d0 2889 ctr = 1 2890 sum_zg = 0 2891 ratio_zg = 0.0 2892 ! 2893 DO p = 1, nat3 ! those are for \k \a, \k' \a' 2894 DO qp = 1, nat3 2895 IF (p == qp) THEN 2896 ratio_zg(ctr) = sum_error_B(ctr) / sum_diag_B(ctr) 2897 !!! WRITE(*,*) "Error from each branch", sum_diag_B(ctr), sum_error_B(ctr), p , qp, ratio_zg(ctr) 2898 IF (ABS(ratio_zg(ctr)) < error_thresh) THEN 2899 sum_zg = sum_zg + 1 2900 ENDIF 2901 ENDIF 2902 ctr = ctr + 1 2903 ENDDO 2904 ENDDO 2905 ! 2906 WRITE(*,*) "Total error:", SUM(ABS(ratio_zg)) / nat3 2907 ENDIF ! compute_error 2908 ! 2909 !IF (sum_zg == nat3) THEN 2910 IF (SUM(ABS(ratio_zg)) / nat3 < error_thresh ) THEN 2911 ctrA = 0 2912 ctrB = 0 2913 ! here we distinguish between q-points in A and B to perform the appropriate summations for the atomic displacements 2914 DO qp = 1, ctrAB 2915 q_A = q(:, qp) + q(:, qp) ! q_A to find IF q belongs in A 2916 IF (((ABS(q_A(1)) < eps) .OR. (ABS(ABS(q_A(1)) - 1) < eps)) .AND. & 2917 ((ABS(q_A(2)) < eps) .OR. (ABS(ABS(q_A(2)) - 1) < eps)) .AND. & 2918 ((ABS(q_A(3)) < eps) .OR. (ABS(ABS(q_A(3)) - 1) < eps))) THEN 2919 ctrA = ctrA + 1 2920 Cx_matA(:, ctrA) = Cx_matAB(:, qp) 2921 Cpx_matA(:, ctrA) = Cpx_matAB(:, qp) 2922 z_nq_A(:, :, ctrA) = z_zg(:, :, qp) 2923 qA(ctrA, :) = q(:, qp) 2924 IF (ABS(qA(ctrA, 1)) < eps) qA(ctrA, 1) = 0.0 2925 IF (ABS(qA(ctrA, 2)) < eps) qA(ctrA, 2) = 0.0 2926 IF (ABS(qA(ctrA, 3)) < eps) qA(ctrA, 3) = 0.0 2927 ELSE 2928 ctrB = ctrB + 1 2929 Cx_matB(:, ctrB) = Cx_matAB(:, qp) 2930 Cpx_matB(:, ctrB) = Cpx_matAB(:, qp) 2931 z_nq_B(:, :, ctrB) = z_zg(:, :, qp) 2932 qB(ctrB,:) = q(:, qp) 2933 IF (ABS(qB(ctrB, 1)) < eps) qB(ctrB, 1) = 0.0 2934 IF (ABS(qB(ctrB, 2)) < eps) qB(ctrB, 2) = 0.0 2935 IF (ABS(qB(ctrB, 3)) < eps) qB(ctrB, 3) = 0.0 2936 ! 2937 ENDIF 2938 ENDDO 2939 ! 2940 IF (single_phonon_displ) THEN 2941 WRITE(*,*) "Print single phonon displacements" 2942 CALL single_phonon(nq_tot, nat, ctrB, ctrA, nat3, ityp, ntyp, & 2943 ntypx, qA, qB, amass, atm, equil_p, & 2944 Rlist, z_nq_B, z_nq_A, Cx_matB, & 2945 Cx_matA, Cpx_matB, Cpx_matA) 2946 ENDIF 2947 ! 2948 WRITE(*,*) 2949 WRITE(*,*) "Print ZG configuration" 2950 IF (compute_error) THEN 2951 WRITE(80,*) "Sum of diagonal terms per q-point:", DBLE(SUM(sum_diag_B) / ctrAB) 2952 WRITE(80,*) "Error and loop index:", SUM(ABS(ratio_zg)) / nat3, kk ! 2953 ENDIF 2954 !WRITE(80,*) "Sum of error per q-point and loop index:", SUM(sum_error_B)/ctrAB, kk ! 2955 WRITE(80,'(A20, 1F6.2,A2)') 'Temperature is: ' , T ,' K' 2956 WRITE(80,*) "Atomic positions", nat * nq_tot 2957 WRITE(81,*) "ZG-Velocities (Ang/ps)" 2958 ! Generate displacements and velocities. 2959 ! Remember nq_tot is also equal to the number of cells 2960 ! Here the displacements are generated according to 2961 ! Np^(- 1/2)(Mo/Mk)^(1/2)[\sum_{q \in B} e^{1qR_p}e^v_{ka}(q)(x_{qv}+y_{q\nu}) 2962 ! z_zg(nat3, nat3, nq)) 2963 ! 2964 D_tau = 0.0d0 2965 P_tau = 0.0d0 2966 ! 2967 ! Main loop to construct ZG configuration 2968 DO p = 1, nq_tot 2969 ctr = 1 2970 DO k = 1, nat ! k represents the atom 2971 nta = ityp(k) 2972 DO i = 1, 3 ! i is for cart directions 2973 ! 2974 DO qp = 1, ctrB 2975 dotp = 0.0d0 2976 DO ii = 1, 3 2977 dotp = dotp + qB(qp, ii) * Rlist(p, ii)! dot product between q and R 2978 ENDDO 2979 DO j = 1, nat3 2980 D_tau(p, k, i) = D_tau(p, k, i) + SQRT(2.0d0 / nq_tot / amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 2981 * z_nq_B(ctr, j, qp) * (1.d0 + imagi) * Cx_matB(j, qp)) 2982 P_tau(p, k, i) = P_tau(p, k, i) + SQRT(2.0d0 / nq_tot * amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 2983 * z_nq_B(ctr, j, qp) * (1.d0 + imagi) * Cpx_matB(j, qp)) / (amass(nta) * AMU_SI) 2984 ! Here we calculate the momenta of the nuclei and finally 2985 !we divide by (amass(nta) *AMU_SI) to get the velocities. 2986 ENDDO 2987 ENDDO ! qp loop 2988 ! 2989 IF (incl_qA) THEN ! If we want to include modes in set A. Those modes are known 2990 ! to break the degeneracy for finite size systems 2991 DO qp = 1, ctrA 2992 dotp = 0.0d0 2993 DO ii = 1, 3 2994 dotp = dotp + qA(qp, ii) * Rlist(p, ii)! 2995 ENDDO 2996 DO j = 1, nat3 2997 D_tau(p, k, i) = D_tau(p, k, i) + SQRT(1.0d0 / nq_tot / amass(nta)) * cos(2.0d0 * pi * dotp) & 2998 * DBLE(z_nq_A(ctr, j, qp) * Cx_matA(j, qp)) ! 2999 P_tau(p, k, i) = P_tau(p, k, i) + SQRT(1.0d0 / nq_tot * amass(nta)) * cos(2.0d0 * pi * dotp) & 3000 * DBLE(z_nq_A(ctr, j, qp) * Cpx_matA(j, qp)) / (amass(nta) * AMU_SI) ! 3001 ENDDO 3002 ENDDO 3003 ENDIF ! IF incl_qA 3004 ! 3005 ctr = ctr + 1 ! for k and i 3006 IF (ABS(D_tau(p, k, i)) .GT. 5) CALL errore('ZG', 'Displacement very large', D_tau(p, k, i) ) 3007 D_tau(p, k, i) = equil_p(p, k, i) + D_tau(p, k, i) ! add equil structure 3008 ENDDO ! end i for cart directions 3009 ENDDO ! end k loop over nat 3010 ENDDO ! end p loop over unit cells 3011 ! print displacements 3012 IF (ionode) THEN 3013 DO p = 1, nq_tot 3014 DO k = 1, nat ! k represents the atom 3015 WRITE(80,'(A6, 3F13.8)') atm(ityp(k)), D_tau(p, k, :) 3016 WRITE(*,'(A10, A6, 3F13.8)') "ZG_conf:", atm(ityp(k)), D_tau(p, k, :) 3017 WRITE(81,'(A6, 3F15.8)') atm(ityp(k)), P_tau(p, k, :) * 1.0E-12 ! multiply to obtain picoseconds 3018 ENDDO 3019 ENDDO 3020 !WRITE(80,*) "sign matrices" 3021 !DO i = 1, pn 3022 ! WRITE(80,*) Mx_mat(i,:) 3023 !ENDDO 3024 WRITE(80,*) 'Anisotropic displacement tensor vs exact values:' 3025 WRITE(81,*) 'ZG-velocities vs exact velocities from momentum operator in second quantization:' 3026 ! 3027 ! Exact anisotropic displacement parameter 3028 DW_fact = 0.0d0 3029 DWp_fact = 0.0d0 3030 ctr = 1 3031 DO k = 1, nat ! k represents the atom 3032 nta = ityp(k) 3033 DO i = 1, 3 ! i is for cart directions 3034 ! 3035 DO qp = 1, ctrB 3036 DO j = 1, nat3 3037 DW_fact(k, i) = DW_fact(k, i) + 2.0d0 / DBLE(nq_tot * amass(nta)) * z_nq_B(ctr, j, qp) & 3038 * CONJG(z_nq_B(ctr, j, qp)) * Cx_matB(j, qp)**2 3039 DWp_fact(k, i) = DWp_fact(k, i) + 2.0d0 * amass(nta) / DBLE(nq_tot) * z_nq_B(ctr, j, qp) & 3040 * CONJG(z_nq_B(ctr, j, qp)) * Cpx_matB(j, qp)**2 & 3041 / ((amass(nta) * AMU_SI)**2) * 1.0E-24 3042 ENDDO 3043 ENDDO 3044 ! 3045 IF (incl_qA) THEN ! 3046 DO qp = 1, ctrA 3047 DO j = 1, nat3 3048 DW_fact(k, i) = DW_fact(k, i) + 1.0d0 / DBLE(nq_tot * amass(nta)) * z_nq_A(ctr, j, qp) & 3049 * CONJG(z_nq_A(ctr, j, qp)) * Cx_matA(j, qp)**2 3050 DWp_fact(k, i) = DWp_fact(k, i) + amass(nta) / DBLE(nq_tot) * z_nq_A(ctr, j, qp) * CONJG(z_nq_A(ctr, j, qp)) & 3051 * Cpx_matA(j, qp)**2 / ((amass(nta) * AMU_SI)**2) * 1.0E-24 3052 ENDDO 3053 ENDDO 3054 ENDIF 3055 ! 3056 ctr = ctr + 1 ! for k and i 3057 ENDDO ! end i for cart directions 3058 ENDDO ! end k loop over nat 3059 ! 3060 T_fact(:,:) = 0.d0 3061 Tp_fact(:,:) = 0.d0 3062 DO k = 1, nat 3063 nta = ityp(k) 3064 DO p = 1, nq_tot 3065 T_fact(k, 1) = T_fact(k, 1) + (D_tau(p, k, 1) - equil_p(p, k, 1))**2 / nq_tot 3066 T_fact(k, 2) = T_fact(k, 2) + (D_tau(p, k, 2) - equil_p(p, k, 2))**2 / nq_tot 3067 T_fact(k, 3) = T_fact(k, 3) + (D_tau(p, k, 3) - equil_p(p, k, 3))**2 / nq_tot 3068 Tp_fact(k, 1) = Tp_fact(k, 1) + (P_tau(p, k, 1))**2 / nq_tot * 1.0E-24 3069 Tp_fact(k, 2) = Tp_fact(k, 2) + (P_tau(p, k, 2))**2 / nq_tot * 1.0E-24 3070 Tp_fact(k, 3) = Tp_fact(k, 3) + (P_tau(p, k, 3))**2 / nq_tot * 1.0E-24 3071 ENDDO 3072 WRITE(80,'(A6, 2i1)') "Atom: ", k 3073 WRITE(80,'(A6, 3F11.6)') atm(nta), T_fact(k, 1) * 8 * pi**2, T_fact(k, 2) * 8 * pi**2, T_fact(k, 3) * 8 * pi**2 3074 WRITE(80,'(A20, 3F11.6)') "Exact values" 3075 WRITE(80,'(A6, 3F11.6)') atm(nta), DW_fact(k, 1) * 8 * pi**2, DW_fact(k, 2) * 8 * pi**2, DW_fact(k, 3) * 8 * pi**2 3076 !Here we print the DW velocities, in the same spirit that with DW factor. see p.237 of Maradudin Book 3077 WRITE(81,'(A6, 3F12.8)') atm(nta), SQRT(Tp_fact(k, 1)), SQRT(Tp_fact(k, 2)), SQRT(Tp_fact(k, 3)) 3078 WRITE(81,'(A6, 3F12.8)') atm(nta), SQRT(DWp_fact(k, 1)), SQRT(DWp_fact(k, 2)), SQRT(DWp_fact(k, 3)) 3079 ENDDO 3080 ! 3081 ! off-diagonal terms of tensor 3082 WRITE(80,*) "off-diagonal terms" 3083 T_fact(:,:) = 0.d0 3084 Tp_fact(:,:) = 0.d0 3085 DO k = 1, nat 3086 nta = ityp(k) 3087 DO p = 1, nq_tot 3088 T_fact(k, 1) = T_fact(k, 1) + (D_tau(p, k, 1) - equil_p(p, k, 1)) * & 3089 (D_tau(p, k, 2) - equil_p(p, k, 2)) / nq_tot 3090 T_fact(k, 2) = T_fact(k, 2) + (D_tau(p, k, 1) - equil_p(p, k, 1)) * & 3091 (D_tau(p, k, 3) - equil_p(p, k, 3)) / nq_tot 3092 T_fact(k, 3) = T_fact(k, 3) + (D_tau(p, k, 2) - equil_p(p, k, 2)) * & 3093 (D_tau(p, k, 3) - equil_p(p, k, 3)) / nq_tot 3094 ENDDO 3095 WRITE(80,'(A6, 3F11.6)') atm(nta), T_fact(k, 1) * 8 * pi**2, T_fact(k, 2) * 8 * pi**2, T_fact(k, 3) * 8 * pi**2 3096 !! WRITE(80,'(A6, 3F11.6)') atm(nta), DW_fact(k, 1) *8*pi**2, DW_fact(k, 2) *8*pi**2, DW_fact(k, 3) *8*pi**2 3097 ENDDO 3098 ENDIF ! (ionode) 3099 EXIT ! exit kk-loop if the error is less than a threshold 3100 ENDIF 3101!!!!! 3102 ENDDO ! end kk for nloops 3103 CLOSE(80) ! close ZG-configuration file 3104 CLOSE(81) ! close ZG-velocities file 3105 ! 3106 ! 3107 DEALLOCATE(T_fact, Tp_fact, DW_fact, DWp_fact) 3108 DEALLOCATE(equil_p, Rlist, D_tau, qA, qB, z_nq_A, z_nq_B) 3109 DEALLOCATE(Cx_matA, Cx_matB, Cx_matAB) 3110 DEALLOCATE(Cpx_matA, Cpx_matB, Cpx_matAB) 3111 DEALLOCATE(Mx_mat_or, Mx_mat, M_mat, V_mat, ratio_zg) 3112 IF (compute_error) THEN 3113 DEALLOCATE(R_mat, D_vect, F_vect, Bx_vect, E_vect) 3114 DEALLOCATE(sum_error_D, sum_diag_D, sum_error_B, sum_diag_B, sum_error_B2, sum_diag_B2) 3115 ENDIF 3116 ! 3117 ! 3118END SUBROUTINE ZG_configuration 3119 3120SUBROUTINE create_supercell(at,tau, alat, dimx, dimy, dimz, nat, equil_p) 3121! 3122 USE kinds, ONLY : DP 3123 USE constants, ONLY : BOHR_RADIUS_ANGS 3124 USE cell_base, ONLY : bg 3125 IMPLICIT NONE 3126! 3127! 3128 INTEGER, intent(in) :: dimx, dimy, dimz, nat 3129 REAL(DP), intent(in) :: tau(3, nat), at(3, 3) 3130 REAL(DP), intent(out) :: equil_p(dimx * dimy * dimz, nat, 3) 3131 INTEGER :: i, j, k, ctr, p 3132 REAL(DP) :: alat, crystal_pos(dimx * dimy * dimz, nat, 3), abc(3, nat) 3133 alat = alat * BOHR_RADIUS_ANGS !bohr_to_angst ! to convert them in angstrom ! 3134 abc = tau 3135 ! to convert tau/abc in crystal coordinates 3136 call cryst_to_cart(nat, abc, bg, - 1) 3137 ! 3138 ! 3139 ! 3140 crystal_pos = 0 3141 ctr = 1 3142 ! I put dimz loop first so that I am consistent with espresso, but It doesn't 3143 ! matter 3144 DO i = 0, dimz - 1 3145 DO j = 0, dimy - 1 3146 DO k = 0, dimx - 1 3147 DO p = 1, nat 3148 crystal_pos(ctr, p, 1) = (abc(1, p) + k) / float(dimx) 3149 crystal_pos(ctr, p, 2) = (abc(2, p) + j) / float(dimy) 3150 crystal_pos(ctr, p, 3) = (abc(3, p) + i) / float(dimz) 3151 ENDDO 3152 ctr = ctr + 1 3153 ENDDO 3154 ENDDO 3155 ENDDO 3156 ! 3157 ! 3158 equil_p = 0.d0 3159 DO i = 1, dimx*dimy*dimz 3160 DO p = 1, nat 3161 ! matrix maltiplication to cnvert crystaL COORDINATES to cartesian 3162 equil_p(i, p, 1) = (at(1, 1) * crystal_pos(i, p, 1) + & 3163 at(1, 2) * crystal_pos(i, p, 2) + & 3164 at(1, 3) * crystal_pos(i, p, 3)) * alat * dimx 3165 equil_p(i, p, 2) = (at(2, 1) * crystal_pos(i, p, 1) + & 3166 at(2, 2) * crystal_pos(i, p, 2) + & 3167 at(2, 3) * crystal_pos(i, p, 3)) * alat * dimy 3168 equil_p(i, p, 3) = (at(3, 1) * crystal_pos(i, p, 1) + & 3169 at(3, 2) * crystal_pos(i, p, 2) + & 3170 at(3, 3) * crystal_pos(i, p, 3)) * alat * dimz 3171 ENDDO 3172 ENDDO 3173 ! 3174 ! 3175 END SUBROUTINE create_supercell 3176 3177 3178SUBROUTINE single_phonon(nq_tot, nat, ctrB, ctrA, nat3, ityp, ntyp, & 3179 ntypx, qA, qB, amass, atm, equil_p, & 3180 Rlist, z_nq_B, z_nq_A, Cx_matB, & 3181 Cx_matA, Cpx_matB, Cpx_matA) 3182! 3183 USE kinds, ONLY : DP 3184 USE constants, ONLY : AMU_SI, pi 3185 IMPLICIT NONE 3186! 3187! 3188 CHARACTER(LEN=3), intent(in) :: atm(ntypx) 3189 INTEGER, intent(in) :: nq_tot, nat, ctrB, ctrA, nat3, ntyp, ntypx 3190 INTEGER, intent(in) :: Rlist(nq_tot, 3) 3191 INTEGER, intent(in) :: ityp(nat) 3192 REAL(DP), intent(in) :: qA(ctrA, 3), qB(ctrB, 3), amass(ntyp), equil_p(nq_tot, nat, 3) 3193 REAL(DP), intent(in) :: Cx_matB(nat3, ctrB), Cx_matA(nat3, ctrA), Cpx_matA(nat3, ctrA), Cpx_matB(nat3, ctrB) 3194 COMPLEX(DP), intent(in) :: z_nq_A(nat3, nat3, ctrA), z_nq_B(nat3, nat3, ctrB) 3195 CHARACTER(len=256) :: filename 3196! 3197 INTEGER :: i, j, k, p, ii, ctr, nta, qp 3198 REAL(DP) :: dotp 3199 REAL(DP), ALLOCATABLE :: D_tau(:, :, :), P_tau(:, :, :) 3200 COMPLEX(DP) :: imagi 3201 ! 3202 ! 3203 imagi = (0.0d0, 1.0d0) !imaginary unit 3204 ! 3205 ALLOCATE(D_tau(nq_tot, nat, 3), P_tau(nq_tot, nat, 3)) 3206 ! 3207 ! 3208 filename = 'single_phonon-displacements.dat' !'.fp' 3209 OPEN (unit = 85, file = filename, status = 'unknown', form = 'formatted') 3210 ! 3211 filename = 'single_phonon-velocities.dat' !'.fp' 3212 OPEN (unit = 86, file = filename, status = 'unknown', form = 'formatted') 3213 ! Main loop to give single phonon displacements / velocities 3214 WRITE(85,'(A50)') "Displaced positions along phonon modes in set B" 3215 DO qp = 1, ctrB 3216 DO j = 1, nat3 3217 WRITE(85,'(A25, 3F8.4, A15, i)') "Phonon mode at q-point", qB(qp, :), " and branch:", j 3218 WRITE(86,'(A25, 3F8.4, A15, i)') "Phonon mode at q-point", qB(qp, :), " and branch:", j 3219 D_tau = 0.0d0 3220 P_tau = 0.0d0 3221 DO p = 1, nq_tot 3222 dotp = 0.0d0 3223 DO ii = 1, 3 3224 dotp = dotp + qB(qp, ii) * Rlist(p, ii)! dot product between q and R 3225 ENDDO 3226 ctr = 1 3227 DO k = 1, nat ! k represents the atom 3228 nta = ityp(k) 3229 DO i = 1, 3 ! i is for cart directions 3230 D_tau(p, k, i) = D_tau(p, k, i) + SQRT(2.0d0 / nq_tot / amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 3231 * z_nq_B(ctr, j, qp) * (1.d0 + imagi) * ABS(Cx_matB(j, qp))) 3232 P_tau(p, k, i) = P_tau(p, k, i) + SQRT(2.0d0 / nq_tot * amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 3233 * z_nq_B(ctr, j, qp) * (1.d0 + imagi) * ABS(Cpx_matB(j, qp))) / (amass(nta) * AMU_SI) 3234 ! Here we calculate the momenta of the nuclei and finally 3235 !we divide by (amass(nta) *AMU_SI) to get the velocities. 3236 ctr = ctr + 1 ! for k and i 3237 IF (ABS(D_tau(p, k, i)) .GT. 5) CALL errore('ZG', 'Displacement very large', D_tau(p, k, i) ) 3238 D_tau(p, k, i) = equil_p(p, k, i) + D_tau(p, k, i) ! add equil structure 3239 ENDDO ! i loop 3240 ! write output data 3241 WRITE(85,'(A6, 3F13.8)') atm(ityp(k)), D_tau(p, k, :) 3242 WRITE(86,'(A6, 3F15.8)') atm(ityp(k)), P_tau(p, k,:) * 1.0E-12 ! multiply to obtain picoseconds 3243 ! 3244 ENDDO ! k loop 3245 ENDDO ! p loop 3246 ENDDO ! j loop 3247 ENDDO ! qp loop 3248 ! 3249 WRITE(85,'(A50)') "Displaced positions along phonon modes in set A" 3250 DO qp = 1, ctrA 3251 DO j = 1, nat3 3252 WRITE(85,'(A25, 3F8.4, A15, i)') "Phonon mode at q-point", qA(qp, :), " and branch:", j 3253 WRITE(86,'(A25, 3F8.4, A15, i)') "Phonon mode at q-point", qA(qp, :), " and branch:", j 3254 D_tau = 0.0d0 3255 P_tau = 0.0d0 3256 DO p = 1, nq_tot 3257 dotp = 0.0d0 3258 DO ii = 1, 3 3259 dotp = dotp + qA(qp, ii) * Rlist(p, ii)! dot product between q and R 3260 ENDDO 3261 ctr = 1 3262 DO k = 1, nat ! k represents the atom 3263 nta = ityp(k) 3264 DO i = 1, 3 ! i is for cart directions 3265 D_tau(p, k, i) = D_tau(p, k, i) + SQRT(2.0d0 / nq_tot / amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 3266 * z_nq_A(ctr, j, qp) * (1.d0 + imagi) * ABS(Cx_matA(j, qp))) 3267 P_tau(p, k, i) = P_tau(p, k, i) + SQRT(2.0d0 / nq_tot * amass(nta)) * DBLE(EXP(imagi * 2.0d0 * pi * dotp) & 3268 * z_nq_A(ctr, j, qp) * (1.d0 + imagi) * ABS(Cpx_matA(j, qp))) / (amass(nta) * AMU_SI) 3269 ! Here we calculate the momenta of the nuclei and finally 3270 !we divide by (amass(nta) *AMU_SI) to get the velocities. 3271 ctr = ctr + 1 ! for k and i 3272 IF (ABS(D_tau(p, k, i)) .GT. 5) CALL errore('ZG', 'Displacement very large', D_tau(p, k, i) ) 3273 D_tau(p, k, i) = equil_p(p, k, i) + D_tau(p, k, i) ! add equil structure 3274 ENDDO ! i loop 3275 ! write output data 3276 WRITE(85,'(A6, 3F13.8)') atm(ityp(k)), D_tau(p, k, :) 3277 WRITE(86,'(A6, 3F15.8)') atm(ityp(k)), P_tau(p, k,:) * 1.0E-12 ! multiply to obtain picoseconds 3278 ! 3279 ENDDO ! k loop 3280 ENDDO ! p loop 3281 ENDDO ! j loop 3282 ENDDO ! qp loop 3283 ! 3284! 3285! 3286 DEALLOCATE(D_tau, P_tau) 3287 CLOSE(85) 3288 CLOSE(86) 3289! 3290! 3291END SUBROUTINE single_phonon 3292