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
|
! Copyright (C) 2002-2008 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.
subroutine phononsc
use modmain
use modphonon
use modmpi
use moddelf
implicit none
! local variables
integer is,ia,ja,ias,jas
integer ip,nph,p,n,i
real(8) a,b,t1
real(8) ft(3,maxatoms*maxspecies)
complex(8) z1,z2
! allocatable arrays
real(8), allocatable :: vsmt0(:,:),vsir0(:)
complex(8), allocatable :: dyn(:,:)
! store original parameters
avec0(:,:)=avec(:,:)
natoms0(:)=natoms(:)
atposl0(:,:,:)=atposl(:,:,:)
bfcmt00(:,:,:)=bfcmt0(:,:,:)
mommtfix0(:,:,:)=mommtfix(:,:,:)
tshift0=tshift
tforce0=tforce
autokpt0=autokpt
primcell0=primcell
ngridk0(:)=ngridk(:)
! no shifting of atomic basis allowed
tshift=.false.
! require forces
tforce=.true.
! determine k-point grid size from radkpt
autokpt=.true.
! no reduction to primitive cell
primcell=.false.
! initialise universal variables
call init0
! initialise q-point dependent variables
call init2
! store original parameters
natmtot0=natmtot
idxis0(1:natmtot)=idxis(1:natmtot)
bvec0(:,:)=bvec(:,:)
binv0(:,:)=binv(:,:)
atposc0(:,:,:)=atposc(:,:,:)
ngridg0(:)=ngridg(:)
ngtot0=ngtot
if (allocated(ivg0)) deallocate(ivg0)
allocate(ivg0(3,ngtot0))
ivg0(:,:)=ivg(:,:)
if (allocated(igfft0)) deallocate(igfft0)
allocate(igfft0(ngtot0))
igfft0(:)=igfft(:)
! allocate the Kohn-Sham potential derivative arrays
if (allocated(dvsbs)) deallocate(dvsbs)
n=npmtmax*natmtot
allocate(dvsbs(n))
dvsmt(1:npmtmax,1:natmtot) => dvsbs(1:)
if (allocated(dvsir)) deallocate(dvsir)
allocate(dvsir(ngtot))
! allocate supercell offset vector array
if (allocated(vscph)) deallocate(vscph)
allocate(vscph(3,nqptnr))
! allocate dynamical matrix column
allocate(dyn(3,natmtot))
! begin new phonon task
10 continue
natoms(:)=natoms0(:)
! find a dynamical matrix to calculate
call dyntask(80,filext)
! if nothing more to do then restore original input parameters and return
if (iqph == 0) then
filext='.OUT'
natoms(:)=natoms0(:)
avec(:,:)=avec0(:,:)
atposl(:,:,:)=atposl0(:,:,:)
bfcmt0(:,:,:)=bfcmt00(:,:,:)
mommtfix(:,:,:)=mommtfix0(:,:,:)
tshift=tshift0
tforce=tforce0
autokpt=autokpt0
primcell=primcell0
ngridk(:)=ngridk0(:)
deallocate(ivg0,igfft0)
return
end if
if (mp_mpi) write(*,'("Info(phononsc): working on ",A)') 'DYN'//trim(filext)
! dry run: just generate empty DYN files
if (task == 202) goto 10
! zero the dynamical matrix row
dyn(:,:)=0.d0
! zero the Kohn-Sham potential derivative
dvsmt(:,:)=0.d0
dvsir(:)=0.d0
! check to see if mass is considered infinite
if (spmass(isph) <= 0.d0) goto 20
! loop over phases: 0 = cos and 1 = sin displacements
if ((ivq(1,iqph) == 0).and.(ivq(2,iqph) == 0).and.(ivq(3,iqph) == 0)) then
nph=0
else
nph=1
end if
! initialise or read the charge density and potentials from file
trdstate=(task == 201)
! loop over cos and sin displacements
do p=0,nph
! generate the supercell with negative displacement
call genscph(p,-0.5d0*deltaph)
! run the ground-state calculation
call gndstate
! subsequent calculations will read in this supercell potential
trdstate=.true.
! store the total force for this displacement
do ias=1,natmtot
ft(:,ias)=forcetot(:,ias)
end do
! store the Kohn-Sham potential for this displacement
allocate(vsmt0(npmtmax,natmtot),vsir0(ngtot))
vsmt0(:,:)=vsmt(:,:)
vsir0(:)=vsir(:)
! generate the supercell again with positive displacement
call genscph(p,0.5d0*deltaph)
! run the ground-state calculation again
call gndstate
! compute the complex Kohn-Sham potential derivative with implicit q-phase
call phscdvs(p,vsmt0,vsir0)
deallocate(vsmt0,vsir0)
! Fourier transform the force differences to obtain the dynamical matrix
z1=1.d0/(dble(nscph)*deltaph)
! multiply by i for sin-like displacement
if (p == 1) z1=z1*zi
ias=0
jas=0
do is=1,nspecies
ja=0
do ia=1,natoms0(is)
ias=ias+1
do i=1,nscph
ja=ja+1
jas=jas+1
t1=-dot_product(vqc(:,iqph),vscph(:,i))
z2=z1*cmplx(cos(t1),sin(t1),8)
do ip=1,3
t1=-(forcetot(ip,jas)-ft(ip,jas))
dyn(ip,ias)=dyn(ip,ias)+z2*t1
end do
end do
end do
end do
end do
20 continue
! write dynamical matrix row to file
if (mp_mpi) then
ias=0
do is=1,nspecies
do ia=1,natoms0(is)
ias=ias+1
do ip=1,3
a=dble(dyn(ip,ias))
b=aimag(dyn(ip,ias))
if (abs(a) < 1.d-12) a=0.d0
if (abs(b) < 1.d-12) b=0.d0
write(80,'(2G18.10," : is = ",I4,", ia = ",I4,", ip = ",I4)') a,b,is, &
ia,ip
end do
end do
end do
close(80)
! write the complex Kohn-Sham potential derivative to file
natoms(:)=natoms0(:)
natmtot=natmtot0
idxis(1:natmtot)=idxis0(1:natmtot)
ngridg(:)=ngridg0(:)
call writedvs(filext)
end if
! delete the non-essential files
call delfiles(evec=.true.,eval=.true.,occ=.true.)
! synchronise MPI processes
call mpi_barrier(mpicom,ierror)
goto 10
end subroutine
|
