File: curden.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (139 lines) | stat: -rw-r--r-- 3,286 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139

! Copyright (C) 2018 S. Sharma, J. K. Dewhurst and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine curden(afield)
use modmain
use modmpi
use modomp
implicit none
! arguments
real(8), intent(in) :: afield(3)
! local variables
integer ik,is,ias,ispn
integer nr,nri,iro,np
integer nrc,nrci,npc
integer ir,n,i,nthd
real(8) ca,t1
! allocatable arrays
real(8), allocatable :: rfmt(:)
! external functions
real(8) rfint
external rfint
! coupling constant of the external A-field (1/c)
ca=1.d0/solsc
! set the current density to zero
do i=1,3
  do ias=1,natmtot
    is=idxis(ias)
    cdmt(1:npcmt(is),ias,i)=0.d0
  end do
end do
cdir(:,:)=0.d0
call omp_hold(nkpt/np_mpi,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ik=1,nkpt
! distribute among MPI processes
  if (mod(ik-1,np_mpi).ne.lp_mpi) cycle
  call curdenk(ik)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! convert muffin-tin current density to spherical harmonics
call omp_hold(natmtot,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(rfmt,is,nrc,nrci,npc,i) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ias=1,natmtot
  allocate(rfmt(npcmtmax))
  is=idxis(ias)
  nrc=nrcmt(is)
  nrci=nrcmti(is)
  npc=npcmt(is)
  do i=1,3
    rfmt(1:npc)=cdmt(1:npc,ias,i)
    call rfsht(nrc,nrci,rfmt,cdmt(:,ias,i))
  end do
  deallocate(rfmt)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! symmetrise the current density
call symrvf(.false.,.true.,nrcmt,nrcmti,npcmt,npmtmax,cdmt,cdir)
! convert the current density from a coarse to a fine radial mesh
call omp_hold(3,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do i=1,3
  call rfmtctof(cdmt(:,:,i))
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! add current densities from each process and redistribute
if (np_mpi.gt.1) then
  n=npmtmax*natmtot*3
  call mpi_allreduce(mpi_in_place,cdmt,n,mpi_double_precision,mpi_sum,mpicom, &
   ierror)
  n=ngtot*3
  call mpi_allreduce(mpi_in_place,cdir,n,mpi_double_precision,mpi_sum,mpicom, &
   ierror)
end if
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
! add vector potential contribution to make current gauge invariant
call omp_hold(natmtot,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(rfmt,is,nr,nri,np) &
!$OMP PRIVATE(iro,ispn,i,ir,t1) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ias=1,natmtot
  allocate(rfmt(npmtmax))
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  np=npmt(is)
  iro=nri+1
! remove the core density from the muffin-tin density
  call dcopy(np,rhomt(:,ias),1,rfmt,1)
  do ispn=1,nspncr
    i=1
    do ir=1,nri
      rfmt(i)=rfmt(i)-rhocr(ir,ias,ispn)
      i=i+lmmaxi
    end do
    do ir=iro,nr
      rfmt(i)=rfmt(i)-rhocr(ir,ias,ispn)
      i=i+lmmaxo
    end do
  end do
  do i=1,3
    t1=-ca*afield(i)
    call daxpy(np,t1,rfmt,1,cdmt(:,ias,i),1)
  end do
  deallocate(rfmt)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
do i=1,3
  t1=-ca*afield(i)
  call daxpy(ngtot,t1,rhoir,1,cdir(:,i),1)
end do
! compute the total current in the unit cell
do i=1,3
  curtot(i)=rfint(cdmt(:,:,i),cdir(:,i))
end do
! total current magnitude
curtotm=sqrt(curtot(1)**2+curtot(2)**2+curtot(3)**2)
return
end subroutine