File: forcek.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 (168 lines) | stat: -rw-r--r-- 4,833 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
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168

! Copyright (C) 2002-2006 J. K. Dewhurst, S. Sharma and C. Ambrosch-Draxl.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

!BOP
! !ROUTINE: forcek
! !INTERFACE:
subroutine forcek(ik)
! !USES:
use modmain
! !INPUT/OUTPUT PARAMETERS:
!   ik : reduced k-point number (in,integer)
! !DESCRIPTION:
!   Computes the {\bf k}-dependent contribution to the incomplete basis set
!   (IBS) force. See the calling routine {\tt force} for a full description.
!
! !REVISION HISTORY:
!   Created June 2006 (JKD)
!   Updated for spin-spiral case, May 2007 (Francesco Cricchio and JKD)
!EOP
!BOC
implicit none
! arguments
integer, intent(in) :: ik
! local variables
integer ispn0,ispn1,ispn,jspn
integer n,nm,nm2,is,ias,ist,jst
integer iv(3),jv(3),ig,i,j,k,l
real(8) vj(3),sum,t1
complex(8) z1,z2
! allocatable arrays
integer, allocatable :: ijg(:)
real(8), allocatable :: dp(:),evalfv(:,:)
complex(8), allocatable :: apwalm(:,:,:,:)
complex(8), allocatable :: evecfv(:,:,:),evecsv(:,:)
complex(8), allocatable :: h(:),o(:),dlh(:),dlo(:)
complex(8), allocatable :: vh(:),vo(:),ffv(:,:),y(:)
! external functions
complex(8) zdotc
external zdotc
nm2=nmatmax**2
! allocate local arrays
allocate(ijg(nm2),dp(nm2))
allocate(evalfv(nstfv,nspnfv))
allocate(apwalm(ngkmax,apwordmax,lmmaxapw,natmtot))
allocate(evecfv(nmatmax,nstfv,nspnfv))
allocate(h(nm2),o(nm2),dlh(nm2),dlo(nm2))
allocate(vh(nmatmax),vo(nmatmax))
allocate(ffv(nstfv,nstfv),y(nstfv))
! get the eigenvalues/vectors from file
call getevalfv(filext,ik,vkl(:,ik),evalfv)
call getevecfv(filext,ik,vkl(:,ik),vgkl(:,:,:,ik),evecfv)
if (tevecsv) then
  allocate(evecsv(nstsv,nstsv))
  call getevecsv(filext,ik,vkl(:,ik),evecsv)
end if
! loop over first-variational spin components
do jspn=1,nspnfv
  if (spinsprl) then
    ispn0=jspn; ispn1=jspn
  else
    ispn0=1; ispn1=nspinor
  end if
  n=ngk(jspn,ik)
  nm=nmat(jspn,ik)
  do j=1,n
    k=(j-1)*nm
    jv(:)=ivg(:,igkig(j,jspn,ik))
    vj(:)=0.5d0*vgkc(:,j,jspn,ik)
    do i=1,j
      k=k+1
      iv(:)=ivg(:,igkig(i,jspn,ik))-jv(:)
      ijg(k)=ivgig(iv(1),iv(2),iv(3))
      dp(k)=dot_product(vgkc(:,i,jspn,ik),vj(:))
    end do
  end do
! find the matching coefficients
  call match(n,gkc(:,jspn,ik),tpgkc(:,:,jspn,ik),sfacgk(:,:,jspn,ik),apwalm)
! loop over species and atoms
  do ias=1,natmtot
    is=idxis(ias)
! Hamiltonian and overlap matrices
    h(:)=0.d0
    call hmlaa(ias,n,apwalm(:,:,:,ias),nm,h)
    call hmlalo(ias,n,apwalm(:,:,:,ias),nm,h)
    o(:)=0.d0
    call olpaa(ias,n,apwalm(:,:,:,ias),nm,o)
    call olpalo(ias,n,apwalm(:,:,:,ias),nm,o)
! loop over Cartesian directions
    do l=1,3
! APW-APW contribution
      do j=1,n
        k=(j-1)*nm
        do i=1,j
          k=k+1
          ig=ijg(k)
          t1=vgc(l,ig)
          z1=-ffacg(ig,is)*conjg(sfacg(ig,ias))
          z2=t1*(dp(k)*z1+h(k))
          dlh(k)=cmplx(-aimag(z2),dble(z2),8)
          z2=t1*(z1+o(k))
          dlo(k)=cmplx(-aimag(z2),dble(z2),8)
        end do
      end do
      do j=n+1,nm
        k=(j-1)*nm
! APW-local-orbital contribution
        do i=1,n
          k=k+1
          t1=vgkc(l,i,jspn,ik)
          z1=t1*h(k)
          dlh(k)=cmplx(-aimag(z1),dble(z1),8)
          z1=t1*o(k)
          dlo(k)=cmplx(-aimag(z1),dble(z1),8)
        end do
! zero the local-orbital-local-orbital contribution
        do i=n+1,j
          k=k+1
          dlh(k)=0.d0
          dlo(k)=0.d0
        end do
      end do
! compute the force matrix elements in the first-variational basis
      do jst=1,nstfv
        call zhemv('U',nm,zone,dlh,nm,evecfv(:,jst,jspn),1,zzero,vh,1)
        call zhemv('U',nm,zone,dlo,nm,evecfv(:,jst,jspn),1,zzero,vo,1)
        t1=evalfv(jst,jspn)
        do ist=1,nstfv
          z1=zdotc(nm,evecfv(:,ist,jspn),1,vh,1)
          z2=zdotc(nm,evecfv(:,ist,jspn),1,vo,1)
          ffv(ist,jst)=z1-t1*z2
        end do
      end do
! compute the force using the second-variational coefficients if required
      sum=0.d0
      if (tevecsv) then
! spin-polarised case
        do j=1,nstsv
          do ispn=ispn0,ispn1
            i=(ispn-1)*nstfv+1
            call zgemv('N',nstfv,nstfv,zone,ffv,nstfv,evecsv(i,j),1,zzero,y,1)
            z1=zdotc(nstfv,evecsv(i,j),1,y,1)
            sum=sum+occsv(j,ik)*dble(z1)
          end do
        end do
      else
! spin-unpolarised case
        do j=1,nstsv
          sum=sum+occsv(j,ik)*dble(ffv(j,j))
        end do
      end if
!$OMP CRITICAL(forcek_)
      forceibs(l,ias)=forceibs(l,ias)+wkpt(ik)*sum
!$OMP END CRITICAL(forcek_)
! end loop over Cartesian components
    end do
! end loop over atoms and species
  end do
! end loop over first-variational spins
end do
deallocate(ijg,dp,evalfv,apwalm,evecfv)
if (tevecsv) deallocate(evecsv)
deallocate(h,o,dlh,dlo,vh,vo,ffv,y)
return
end subroutine
!EOC