File: dforce.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 (165 lines) | stat: -rw-r--r-- 5,111 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

! Copyright (C) 2012 J. K. Dewhurst, S. Sharma 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 dforce(dyn)
use modmain
use modphonon
use modomp
implicit none
! arguments
complex(8), intent(out) :: dyn(3,natmtot)
! local variables
integer ik,is,ias,nthd
integer nr,nri,ir,np,i
complex(8) zsum,z1
! allocatable arrays
complex(8), allocatable :: zrhomt(:,:),zrhoir(:)
complex(8), allocatable :: grhomt(:,:,:),grhoir(:,:)
complex(8), allocatable :: zvclmt(:,:),zvclir(:)
complex(8), allocatable :: gvclmt(:,:,:),gvclir(:,:)
complex(8), allocatable :: zfmt(:),gzfmt(:,:)
! external functions
complex(8) zfmtinp
external zfmtinp
allocate(zrhomt(npmtmax,natmtot),zrhoir(ngtot))
allocate(grhomt(npmtmax,natmtot,3),grhoir(ngtot,3))
allocate(zvclmt(npmtmax,natmtot),zvclir(ngtot))
allocate(gvclmt(npmtmax,natmtot,3),gvclir(ngtot,3))
allocate(zfmt(npmtmax),gzfmt(npmtmax,3))
! make complex copy of the density
do ias=1,natmtot
  is=idxis(ias)
  call rtozfmt(nrmt(is),nrmti(is),rhomt(:,ias),zrhomt(:,ias))
end do
zrhoir(:)=rhoir(:)
! compute the gradient of the density
call gradzf(zrhomt,zrhoir,grhomt,grhoir)
!--------------------------------------------------------------!
!     Hellmann-Feynman force derivative for displaced atom     !
!--------------------------------------------------------------!
nr=nrmt(isph)
nri=nrmti(isph)
np=npmt(isph)
! zero the interstitial density
zrhoir(:)=0.d0
zfmt(1:np)=0.d0
i=1
do ir=1,nri
  zfmt(i)=vcln(ir,isph)
  i=i+lmmaxi
end do
do ir=nri+1,nr
  zfmt(i)=vcln(ir,isph)
  i=i+lmmaxo
end do
call gradzfmt(nr,nri,rsp(:,isph),zfmt,npmtmax,gzfmt)
! compute the q-dependent nuclear Coulomb potential derivative
zvclmt(:,:)=0.d0
zvclmt(1:np,iasph)=gzfmt(1:np,ipph)
tphdyn=.true.
call zpotcoul(nrmt,nrmti,npmt,npmti,nrspmax,rsp,ngridg,igfft,ngvec,gqc,gclgq, &
 ngvec,jlgqrmt,ylmgq,sfacgq,zrhoir,npmtmax,zvclmt,zvclir)
zfmt(1:np)=zvnmt(1:np)
! multiply with density derivative and integrate
zsum=0.d0
do ir=1,ngtot
  zsum=zsum+cfunir(ir)*conjg(zvclir(ir))*drhoir(ir)
end do
zsum=zsum*omega/dble(ngtot)
do ias=1,natmtot
  is=idxis(ias)
  zsum=zsum+zfmtinp(nrmt(is),nrmti(is),rsp(:,is),r2sp(:,is),zvclmt(:,ias), &
   drhomt(:,ias))
end do
dyn(ipph,iasph)=-zsum
! compute the lattice-periodic nuclear Coulomb potential derivative
zvclmt(:,:)=0.d0
zvclmt(1:np,iasph)=gzfmt(1:np,ipph)
call zpotcoul(nrmt,nrmti,npmt,npmti,nrspmax,rsp,ngridg,igfft,ngvec,gc,gclg, &
 ngvec,jlgrmt,ylmg,sfacg,zrhoir,npmtmax,zvclmt,zvclir)
tphdyn=.false.
! multiply with density gradient and integrate
zsum=0.d0
do ir=1,ngtot
  zsum=zsum+cfunir(ir)*zvclir(ir)*grhoir(ir,ipph)
end do
zsum=zsum*omega/dble(ngtot)
do ias=1,natmtot
  is=idxis(ias)
  zsum=zsum+zfmtinp(nrmt(is),nrmti(is),rsp(:,is),r2sp(:,is),zvclmt(:,ias), &
   grhomt(:,ias,ipph))
end do
dyn(ipph,iasph)=dyn(ipph,iasph)-zsum
! nuclear-nuclear term
zvclmt(1:np,iasph)=zvnmt(1:np)-zfmt(1:np)
call gradzf(zvclmt,zvclir,gvclmt,gvclir)
do ias=1,natmtot
  is=idxis(ias)
  z1=spzn(is)*gvclmt(1,ias,ipph)*y00
  dyn(ipph,iasph)=dyn(ipph,iasph)+z1
end do
!-------------------------------------------------------------------!
!     Hellmann-Feynman force derivative for non-displaced atoms     !
!-------------------------------------------------------------------!
do ias=1,natmtot
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  np=npmt(is)
! remove the gradient part of the Coulomb potential for displaced muffin-tin
  if (ias.eq.iasph) then
    call rtozfmt(nr,nri,vclmt(:,iasph),zfmt)
    call gradzfmt(nr,nri,rsp(:,isph),zfmt,npmtmax,gzfmt)
    dvclmt(1:np,ias)=dvclmt(1:np,ias)+gzfmt(1:np,ipph)
  end if
! compute the gradient of the Coulomb potential derivative at the nucleus
  call gradzfmt(nr,nri,rsp(:,is),dvclmt(:,ias),npmtmax,gzfmt)
  do i=1,3
    if ((ias.eq.iasph).and.(i.eq.ipph)) cycle
    dyn(i,ias)=spzn(is)*gzfmt(1,i)*y00
  end do
end do
!--------------------------------------------!
!     IBS correction to force derivative     !
!--------------------------------------------!
! k-point dependent part
call omp_hold(nkptnr,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP NUM_THREADS(nthd)
!$OMP DO
do ik=1,nkptnr
  call dforcek(ik,dyn)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! k-point independent part
do ias=1,natmtot
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  np=npmt(is)
  do i=1,3
    z1=zfmtinp(nr,nri,rsp(:,is),r2sp(:,is),grhomt(:,ias,i),dvsmt(:,ias))
    dyn(i,ias)=dyn(i,ias)-z1
  end do
! convert Kohn-Sham potential to complex spherical harmonics
  call rtozfmt(nr,nri,vsmt(:,ias),zfmt)
! remove the gradient part from the density derivative for displaced muffin-tin
  if (ias.eq.iasph) then
    drhomt(1:np,ias)=drhomt(1:np,ias)+grhomt(1:np,ias,ipph)
  end if
! compute the gradient of the density derivative
  call gradzfmt(nr,nri,rsp(:,is),drhomt(:,ias),npmtmax,gzfmt)
  do i=1,3
    z1=zfmtinp(nr,nri,rsp(:,is),r2sp(:,is),zfmt,gzfmt(:,i))
    dyn(i,ias)=dyn(i,ias)-z1
  end do
end do
deallocate(zrhomt,zrhoir,grhomt,grhoir)
deallocate(zvclmt,zvclir,gvclmt,gvclir,zfmt,gzfmt)
return
end subroutine