1! 2! Copyright (C) 1996-2016 The SIESTA group 3! This file is distributed under the terms of the 4! GNU General Public License: see COPYING in the top directory 5! or http://www.gnu.org/copyleft/gpl.txt. 6! See Docs/Contributors.txt for a list of contributors. 7! 8 subroutine obc( polR, ucell, dx, nspin) 9C ******************************************************************* 10C Writes the born effective charges to file 11C by Tom Archer 12C Modified by Alberto Garcia, April 2007 13C ********* INPUT *************************************************** 14C real*8 dx : atomic displacements (in Bohr) 15C real*8 polR(3,nspin) : polarization along lattice vectors(Bohr) 16C real*8 ucell(3,3) : cell vectors 17C integer nspin : spin polarized calculation flag 18C ********** BEHAVIOUR ********************************************** 19C On the first call (undisplaced coordinates), initial polarization 20C is calculated. 21C The born charges are calculated as the difference between the 22C polarization of the undisplaced and the displaced coordinates 23C divided by the displacement 24C 25C a phase is introduced in the KSV calculation which is removed before 26C the BC matrix is saved. 27C ******************************************************************* 28 29 use precision, only: dp 30 use fdf, only: fdf_string 31 use alloc, only: re_alloc 32 33 implicit none 34 35 integer, intent(in) :: nspin 36 real(dp), intent(in) :: polR(3,nspin), 37 $ ucell(3,3), dx 38 external io_assign, io_close 39 40c Internal 41 42 integer igrd, ispin, ix, unit1 43 44 character(len=33), save :: fname 45 logical, save :: first_time = .true. 46 47 real(dp) :: dmod(3), phaseR(3), uR(3,3), pR(3), pxyz(3) 48 49 double precision, dimension(:,:), pointer, save :: pres 50 51!--------- 52 53 do ispin=1,nspin 54 if (nspin.gt.1) then 55 if (ispin.eq.1) write(6,'(/,a)') 56 . 'obc: Macroscopic polarization for spin Up:' 57 if (ispin.eq.2) write(6,'(/,a)') 58 . 'obc: Macroscopic polarization for spin Down:' 59 endif 60 write(6,'(/,a,a)') 61 . 'obc: Macroscopic polarization per unit cell', 62 . ' along the lattice vectors (Bohr): ' 63 write(6,'(a,3f12.6)') 'obc:',(polR(ix,ispin),ix=1,3) 64 enddo 65 66 if (first_time) then 67 68 fname = trim(fdf_string('SystemLabel','siesta')) // '.BC' 69 70 nullify( pres ) 71 call re_alloc( pres, 1, 3, 1, nspin, name='pres', 72 & routine='obc' ) 73 74 call io_assign(unit1) 75 open( unit1, file=fname, status='unknown' ) 76 rewind(unit1) 77 if (nspin.eq.1) then 78 write(unit1,'(a)') 'BC matrix' 79 else 80 write(unit1,'(a)') 'BC matrix for spin-up + spin-down' 81 endif 82 ! Set values of residual polarization 83 do ispin=1,nspin 84 do igrd=1,3 85 pres(igrd,ispin) = polR(igrd,ispin) 86 enddo 87 enddo 88 call io_close(unit1) 89 90 first_time = .false. 91 92 else 93 94 !matrix to convert from lattice vectors to xyz 95 do igrd=1,3 96 dmod(igrd)=dot_product(ucell(:,igrd),ucell(:,igrd)) 97 dmod(igrd)=sqrt(dmod(igrd)) 98 do ix=1,3 99 uR(ix,igrd)=ucell(ix,igrd)/dmod(igrd) 100 enddo 101 enddo 102 103 ! write BC matrix in lattice vectors 104 pR(:)=0.0_dp 105 pxyz(:)=0.0_dp 106 do ispin=1,nspin 107 do igrd=1,3 108 pR(igrd)=pR(igrd)+(polR(igrd,ispin)-pres(igrd,ispin)) 109 enddo 110 enddo 111 112 !remove extra phase 113 do igrd=1,3 114 phaseR(igrd) = nint((pR(igrd)/dmod(igrd))) * dmod(igrd) 115 pR(igrd) = pR(igrd) - phaseR(igrd) 116 enddo 117 118 !get polarization in cartesian coordinates 119 do ix=1,3 120 pxyz(ix) = uR(ix,1)*pR(1)+ uR(ix,2)*pR(2)+uR(ix,3)*pR(3) 121 enddo 122 123 call io_assign(unit1) 124 open( unit1, file=fname, status='old',position="append", 125 $ action="write") 126 write(unit1,'(3f15.7)') (pxyz(ix)/dx,ix=1,3) 127 call io_close(unit1) 128 129 endif 130 131 end 132