1 2! Copyright (C) 2009 T. McQueen and J. K. Dewhurst. 3! This file is distributed under the terms of the GNU General Public License. 4! See the file COPYING for license details. 5 6!BOP 7! !ROUTINE: ggamt_sp_2b 8! !INTERFACE: 9subroutine ggamt_sp_2b(is,np,g2up,g2dn,gvup,gvdn,vxup,vxdn,vcup,vcdn,dxdgu2, & 10 dxdgd2,dxdgud,dcdgu2,dcdgd2,dcdgud) 11! !USES: 12use modmain 13! !DESCRIPTION: 14! Post processing step of muffin-tin gradients for GGA type 2. See routine 15! {\tt ggamt\_sp\_2a} for full details. 16! 17! !REVISION HISTORY: 18! Created November 2009 (JKD and TMcQ) 19!EOP 20!BOC 21implicit none 22! arguments 23integer, intent(in) :: is,np 24real(8), intent(in) :: g2up(np),g2dn(np) 25real(8), intent(in) :: gvup(np,3),gvdn(np,3) 26real(8), intent(inout) :: vxup(np),vxdn(np) 27real(8), intent(inout) :: vcup(np),vcdn(np) 28real(8), intent(in) :: dxdgu2(np),dxdgd2(np),dxdgud(np) 29real(8), intent(in) :: dcdgu2(np),dcdgd2(np),dcdgud(np) 30! local variables 31integer nr,nri,i 32! allocatable arrays 33real(8), allocatable :: rfmt1(:),rfmt2(:),grfmt(:,:) 34allocate(rfmt1(np),rfmt2(np),grfmt(np,3)) 35nr=nrmt(is) 36nri=nrmti(is) 37!------------------! 38! exchange ! 39!------------------! 40! convert dxdgu2 to spherical harmonics 41call rfsht(nr,nri,dxdgu2,rfmt1) 42! compute grad dxdgu2 43call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 44! (grad dxdgu2).(grad rhoup) in spherical coordinates 45rfmt1(1:np)=0.d0 46do i=1,3 47 call rbsht(nr,nri,grfmt(:,i),rfmt2) 48 rfmt1(1:np)=rfmt1(1:np)+rfmt2(1:np)*gvup(1:np,i) 49end do 50vxup(1:np)=vxup(1:np)-2.d0*(rfmt1(1:np)+dxdgu2(1:np)*g2up(1:np)) & 51 -dxdgud(1:np)*g2dn(1:np) 52! convert dxdgd2 to spherical harmonics 53call rfsht(nr,nri,dxdgd2,rfmt1) 54! compute grad dxdgd2 55call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 56! (grad dxdgd2).(grad rhodn) in spherical coordinates 57rfmt1(1:np)=0.d0 58do i=1,3 59 call rbsht(nr,nri,grfmt(:,i),rfmt2) 60 rfmt1(1:np)=rfmt1(1:np)+rfmt2(1:np)*gvdn(1:np,i) 61end do 62vxdn(1:np)=vxdn(1:np)-2.d0*(rfmt1(1:np)+dxdgd2(1:np)*g2dn(1:np)) & 63 -dxdgud(1:np)*g2up(1:np) 64! convert dxdgud to spherical harmonics 65call rfsht(nr,nri,dxdgud,rfmt1) 66! compute grad dxdgud 67call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 68! (grad dxdgud).(grad rhodn) and (grad dxdgud).(grad rhoup) 69do i=1,3 70 call rbsht(nr,nri,grfmt(:,i),rfmt1) 71 vxup(1:np)=vxup(1:np)-rfmt1(1:np)*gvdn(1:np,i) 72 vxdn(1:np)=vxdn(1:np)-rfmt1(1:np)*gvup(1:np,i) 73end do 74!---------------------! 75! correlation ! 76!---------------------! 77! convert dcdgu2 to spherical harmonics 78call rfsht(nr,nri,dcdgu2,rfmt1) 79! compute grad dcdgu2 80call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 81! (grad dcdgu2).(grad rhoup) in spherical coordinates 82rfmt1(1:np)=0.d0 83do i=1,3 84 call rbsht(nr,nri,grfmt(:,i),rfmt2) 85 rfmt1(1:np)=rfmt1(1:np)+rfmt2(1:np)*gvup(1:np,i) 86end do 87vcup(1:np)=vcup(1:np)-2.d0*(rfmt1(1:np)+dcdgu2(1:np)*g2up(1:np)) & 88 -dcdgud(1:np)*g2dn(1:np) 89! convert dcdgd2 to spherical harmonics 90call rfsht(nr,nri,dcdgd2,rfmt1) 91! compute grad dcdgd2 92call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 93! (grad dcdgd2).(grad rhodn) in spherical coordinates 94rfmt1(1:np)=0.d0 95do i=1,3 96 call rbsht(nr,nri,grfmt(:,i),rfmt2) 97 rfmt1(1:np)=rfmt1(1:np)+rfmt2(1:np)*gvdn(1:np,i) 98end do 99vcdn(1:np)=vcdn(1:np)-2.d0*(rfmt1(1:np)+dcdgd2(1:np)*g2dn(1:np)) & 100 -dcdgud(1:np)*g2up(1:np) 101! convert dcdgud to spherical harmonics 102call rfsht(nr,nri,dcdgud,rfmt1) 103! compute grad dcdgud 104call gradrfmt(nr,nri,rlmt(:,-1,is),wcrmt(:,:,is),rfmt1,np,grfmt) 105! (grad dcdgud).(grad rhodn) and (grad dcdgud).(grad rhoup) 106do i=1,3 107 call rbsht(nr,nri,grfmt(:,i),rfmt1) 108 vcup(1:np)=vcup(1:np)-rfmt1(1:np)*gvdn(1:np,i) 109 vcdn(1:np)=vcdn(1:np)-rfmt1(1:np)*gvup(1:np,i) 110end do 111deallocate(rfmt1,rfmt2,grfmt) 112end subroutine 113!EOC 114 115