File: force.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 (230 lines) | stat: -rw-r--r-- 7,973 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
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230

! Copyright (C) 2002-2005 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: force
! !INTERFACE:
subroutine force
! !USES:
use modmain
use modtest
use modmpi
use modomp
! !DESCRIPTION:
!   Computes the various contributions to the atomic forces. In principle, the
!   force acting on a nucleus is simply the gradient at that site of the
!   classical electrostatic potential from the other nuclei and the electronic
!   density. This is a result of the Hellmann-Feynman theorem. However because
!   the basis set is dependent on the nuclear coordinates and is not complete,
!   the Hellman-Feynman force is inacurate and corrections to it are required.
!   The first is the core correction which arises because the core wavefunctions
!   were determined by neglecting the non-spherical parts of the Kohn-Sham
!   potential $v_s$. Explicitly this is given by
!   $$ {\bf F}_{\rm core}^{\alpha}=\int_{\rm MT_{\alpha}} v_s({\bf r})
!    \nabla\rho_{\rm core}^{\alpha}({\bf r})\,d{\bf r} $$
!   for atom $\alpha$. The second, which is the incomplete basis set (IBS)
!   correction, is due to the position dependence of the APW functions, and is
!   derived by considering the change in total energy if the eigenvector
!   coefficients were fixed and the APW functions themselves were changed. This
!   would result in changes to the first-variational Hamiltonian and overlap
!   matrices given by
!   \begin{align*}
!    \delta H_{\bf G,G'}^{\alpha}&=i({\bf G-G'})
!    \left(H^{\alpha}_{\bf G+k,G'+k}-\frac{1}{2}({\bf G+k})\cdot({\bf G'+k})
!    \tilde{\Theta}_{\alpha}({\bf G-G'})e^{-i({\bf G-G'})\cdot{\bf r}_{\alpha}}
!    \right)\\
!    \delta O_{\bf G,G'}^{\alpha}&=i({\bf G-G'})\left(O^{\alpha}_{\bf G+k,G'+k}
!    -\tilde{\Theta}_{\alpha}({\bf G-G'})e^{-i({\bf G-G'})\cdot{\bf r}_{\alpha}}
!    \right)
!   \end{align*}
!   where both ${\bf G}$ and ${\bf G'}$ run over the APW indices;
!   $\tilde{\Theta}_{\alpha}$ is the form factor of the smooth step function for
!   muffin-tin $\alpha$; and $H^{\alpha}$ and $O^{\alpha}$ are the muffin-tin
!   Hamiltonian and overlap matrices, respectively. The APW-local-orbital part
!   is given by
!   \begin{align*}
!    \delta H_{\bf G,G'}^{\alpha}&=i({\bf G+k})H^{\alpha}_{\bf G+k,G'+k}\\
!    \delta O_{\bf G,G'}^{\alpha}&=i({\bf G+k})O^{\alpha}_{\bf G+k,G'+k}
!   \end{align*}
!   where ${\bf G}$ runs over the APW indices and ${\bf G'}$ runs over the
!   local-orbital indices. There is no contribution from the
!   local-orbital-local-orbital part of the matrices. We can now write the IBS
!   correction in terms of the basis of first-variational states as
!   \begin{align*}
!    {\bf F}_{ij}^{\alpha{\bf k}}=\sum_{\bf G,G'}
!    b^{i{\bf k}*}_{\bf G}b^{j{\bf k}}_{\bf G'}\left(
!    \delta H_{\bf G,G'}^{\alpha}-\epsilon_j\delta O_{\bf G,G'}^{\alpha}\right),
!   \end{align*}
!   where $b^{i{\bf k}}$ is the first-variational eigenvector.
!   Finally, the ${\bf F}_{ij}^{\alpha{\bf k}}$ matrix elements can be
!   multiplied by the second-variational coefficients, and contracted over all
!   indices to obtain the IBS force:
!   \begin{align*}
!    {\bf F}_{\rm IBS}^{\alpha}=\sum_{\bf k}w_{\bf k}\sum_{l\sigma}n_{l{\bf k}}
!    \sum_{ij}c_{\sigma i}^{l{\bf k}*}c_{\sigma j}^{l{\bf k}}
!    {\bf F}_{ij}^{\alpha{\bf k}}
!    +\int_{\rm MT_{\alpha}}v_s({\bf r})\nabla\left[\rho({\bf r})
!    -\rho^{\alpha}_{\rm core}({\bf r})\right]\,d{\bf r},
!   \end{align*}
!   where $c^{l{\bf k}}$ are the second-variational coefficients, $w_{\bf k}$
!   are the $k$-point weights, $n_{l{\bf k}}$ are the occupancies. See routines
!   {\tt hmlaa}, {\tt olpaa}, {\tt hmlalo}, {\tt olpalo}, {\tt energy},
!   {\tt eveqn} and {\tt gencfun}.
!
! !REVISION HISTORY:
!   Created January 2004 (JKD)
!   Fixed problem with second-variational forces, May 2008 (JKD)
!EOP
!BOC
implicit none
! local variables
integer ik,idm,ispn
integer is,ias,nr,nri
integer np,i,nthd
real(8) t1
real(8) ts0,ts1
! allocatable arrays
real(8), allocatable :: rfmt1(:,:),rfmt2(:),grfmt(:,:)
! external functions
real(8) rfmtinp
external rfmtinp
call timesec(ts0)
allocate(grfmt(npmtmax,3))
!---------------------------------!
!     Hellmann-Feynman forces     !
!---------------------------------!
do ias=1,natmtot
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  call gradrfmt(nr,nri,rsp(:,is),vclmt(:,ias),npmtmax,grfmt)
  do i=1,3
    forcehf(i,ias)=-spzn(is)*grfmt(1,i)*y00
  end do
end do
! symmetrise Hellmann-Feynman forces
call symveca(forcehf)
!----------------------------------!
!     IBS correction to forces     !
!----------------------------------!
! set the IBS forces to zero
forceibs(:,:)=0.d0
! compute k-point dependent contribution to the IBS forces
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 forcek(ik)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
! add IBS forces from each process and redistribute
if (np_mpi.gt.1) then
  call mpi_allreduce(mpi_in_place,forceibs,3*natmtot,mpi_double_precision, &
   mpi_sum,mpicom,ierror)
end if
! integral of Kohn-Sham potential with gradient of density
do ias=1,natmtot
  is=idxis(ias)
  nr=nrmt(is)
  nri=nrmti(is)
  call gradrfmt(nr,nri,rsp(:,is),rhomt(:,ias),npmtmax,grfmt)
  do i=1,3
    t1=rfmtinp(nr,nri,rsp(:,is),r2sp(:,is),vsmt(:,ias),grfmt(:,i))
    forceibs(i,ias)=forceibs(i,ias)+t1
  end do
end do
! integral of Kohn-Sham magnetic field with magnetisation gradient
if (spinpol) then
  allocate(rfmt1(npmtmax,natmtot))
  do idm=1,ndmag
    do ias=1,natmtot
      is=idxis(ias)
      call rfsht(nrcmt(is),nrcmti(is),bsmt(:,ias,idm),rfmt1(:,ias))
    end do
    call rfmtctof(rfmt1)
    do ias=1,natmtot
      is=idxis(ias)
      nr=nrmt(is)
      nri=nrmti(is)
      call gradrfmt(nr,nri,rsp(:,is),magmt(:,ias,idm),npmtmax,grfmt)
      do i=1,3
        t1=rfmtinp(nr,nri,rsp(:,is),r2sp(:,is),rfmt1(:,ias),grfmt(:,i))
        forceibs(i,ias)=forceibs(i,ias)+t1
      end do
    end do
  end do
  deallocate(rfmt1)
end if
! integral of Kohn-Sham tau-DFT potential with kinetic energy density gradient
if (xcgrad.eq.4) then
  allocate(rfmt1(npmtmax,natmtot),rfmt2(npmtmax))
  do ias=1,natmtot
    is=idxis(ias)
    call rfsht(nrcmt(is),nrcmti(is),wsmt(:,ias),rfmt1(:,ias))
  end do
  call rfmtctof(rfmt1)
  do ispn=1,nspinor
    do ias=1,natmtot
      is=idxis(ias)
      nr=nrmt(is)
      nri=nrmti(is)
      np=npmt(is)
      rfmt2(1:np)=taumt(1:np,ias,ispn)-taucr(1:np,ias,ispn)
      call gradrfmt(nr,nri,rsp(:,is),rfmt2,npmtmax,grfmt)
      do i=1,3
        t1=rfmtinp(nr,nri,rsp(:,is),r2sp(:,is),rfmt1(:,ias),grfmt(:,i))
        forceibs(i,ias)=forceibs(i,ias)+t1
      end do
    end do
  end do
  deallocate(rfmt1,rfmt2)
end if
! symmetrise IBS forces
call symveca(forceibs)
! total force on each atom
do ias=1,natmtot
  forcetot(:,ias)=forcehf(:,ias)+forceibs(:,ias)
end do
! symmetrise total forces
call symveca(forcetot)
! compute the average force
do i=1,3
  forceav(i)=0.d0
  do ias=1,natmtot
    forceav(i)=forceav(i)+forcetot(i,ias)
  end do
  forceav(i)=forceav(i)/dble(natmtot)
end do
! remove the average force, if required, to prevent translation of atomic basis
if (tfav0) then
  do ias=1,natmtot
    forcetot(:,ias)=forcetot(:,ias)-forceav(:)
  end do
end if
! zero force on atoms with negative mass
do ias=1,natmtot
  is=idxis(ias)
  if (spmass(is).le.0.d0) forcetot(:,ias)=0.d0
end do
! compute maximum force magnitude over all atoms
forcemax=0.d0
do ias=1,natmtot
  t1=sqrt(forcetot(1,ias)**2+forcetot(2,ias)**2+forcetot(3,ias)**2)
  if (t1.gt.forcemax) forcemax=t1
end do
deallocate(grfmt)
call timesec(ts1)
timefor=timefor+ts1-ts0
! write total forces to test file
call writetest(750,'total forces',nv=3*natmtot,tol=1.d-3,rva=forcetot)
return
end subroutine
!EOC