1
!{\src2tex{textfont=tt}}
7
!! Primary routine for conducting DFT calculations of the polarisability
8
!! within the random phase approximation (RPA)
11
!! Copyright (C) 2000-2007 ABINIT group (XG,GMR,MF)
12
!! This file is distributed under the terms of the
13
!! GNU General Public License, see ~abinit/COPYING
14
!! or http://www.gnu.org/copyleft/gpl.txt .
15
!! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
18
!! dtfil <type(datafiles_type)>=variables related to files
19
!! dtset <type(dataset_type)>=all input variables in this dataset
20
!! mband =maximum number of bands
21
!! mgfft =maximum single fft dimension
22
!! mpi_enreg=informations about MPI parallelization
23
!! mpw =maximum number of planewaves in basis sphere (large number)
24
!! natom =number of atoms in unit cell
25
!! nkpt =number of k points
26
!! nspden=number of spin-density components
27
!! nspinor=number of spinorial components of the wavefunctions
28
!! nsppol=number of channels for spin-polarization (1 or 2)
29
!! nsym=number of symmetry elements in space group
30
!! xred(3,natom)=reduced dimensionless atomic coordinates
33
!! (no direct output : results written)
36
!! mkmem =maximum number of k points which can fit in core memory
42
!! distrb2,getfreqsus,getmpw,getng,hdr_clean,inwffil3,ioarr,kpgio,metric
43
!! mkrdim,newocc,setsym,sphereboundary,status,suscep_dyn,suscep_stat,timab
44
!! wffclose,wrtout,xcacfd,xme_init
48
#if defined HAVE_CONFIG_H
52
subroutine suscep(dtfil,dtset,iexit,&
53
& mband,mgfft,mkmem,mpi_enreg,mpw,natom,nfft,nkpt,&
54
& nspden,nspinor,nsppol,nsym,occ,xred,ncdims,ngfft)
59
!This section has been created automatically by the script Abilint (TD). Do not modify these by hand.
60
#ifdef HAVE_FORTRAN_INTERFACES
61
use interfaces_01manage_mpi
64
use interfaces_12geometry
65
use interfaces_13io_mpi
66
use interfaces_13recipspace
67
use interfaces_14iowfdenpot
68
use interfaces_14occeig
69
use interfaces_17suscep, except_this_one => suscep
70
use interfaces_lib01hidempi
72
!End of the abilint section
76
!Arguments ------------------------------------
78
integer,intent(in) :: iexit,mband,mgfft,mpw,natom,nfft,nkpt,nspden,nsppol,nsym
79
integer,intent(inout) :: mkmem,nspinor
80
type(MPI_type),intent(inout) :: mpi_enreg
81
type(datafiles_type),intent(in) :: dtfil
82
type(dataset_type),intent(in) :: dtset
83
type(vardims_type),intent(inout) :: ncdims
85
integer,intent(in) :: ngfft(18)
86
real(dp),intent(in) :: xred(3,natom)
87
real(dp),intent(inout) :: occ(mband*nkpt*nsppol)
89
!Local variables-------------------------------
91
integer,parameter :: level=3
92
integer :: accessfil,dielop,fformr,ierr,ifreq,ii,ipw1,ipw2,isp,isp1,isp2
93
integer :: master,me,mgfftdiel,nband_mx,ncblocks,nfftdiel,nfreqsus,npwdiel
94
integer :: prtvol,rdwr,rdwrpaw,susopt
95
real(dp) :: diecut,dummy,ecutsus,entropy,etotal,fermie,residm,ucvol
96
character(len=500) :: message
97
character(len=fnlen) :: filkgs,filsustr,kgnam
98
type(dens_sym_operator_type) :: densymop_diel
100
type(wffile_type) :: wff1
102
integer :: ngfftdiel(18),npwarr_diel(1),npwtot_diel(1)
103
integer,allocatable :: gbound_diel(:,:),indsym(:,:,:),irrzondiel(:,:,:)
104
integer,allocatable :: kg(:,:),kg_diel(:,:),nband_dum(:),npwarr(:),npwtot(:)
105
integer,allocatable :: symrec(:,:,:)
106
real(dp) :: dielar(7),gmet(3,3),gprimd(3,3),kpt_diel(3),rmet(3,3),rprimd(3,3)
108
real(dp),allocatable :: cg(:,:),dielinv(:,:,:,:,:),doccde(:),eigen(:),freq(:)
109
real(dp),allocatable :: phnonsdiel(:,:,:),rhor(:,:),susd_non_dyn(:,:,:,:)
110
real(dp),allocatable :: susmat(:,:,:,:,:),susmat_dyn(:,:,:,:,:,:),wght_freq(:)
112
!***********************************************************************
114
call timab(84,1,tsec)
115
call status(0,dtfil%filstat,iexit,level,'enter ')
119
mpi_enreg%nproc_fft=1
120
mpi_enreg%paral_fft=0
121
mpi_enreg%paral_level=2
124
! If dtset%accesswff == 2 set all array outputs to netcdf format
127
if (dtset%accesswff == 2) then
130
if (dtset%accesswff == 3) then
136
call xme_init(mpi_enreg,me)
138
if(mpi_enreg%paral_compil_kpt==1)then
139
allocate(mpi_enreg%proc_distrb(nkpt,mband,nsppol))
141
call distrb2(mband, dtset%nband, nkpt, nsppol, mpi_enreg)
144
!Structured debugging if prtvol==-level
146
if(prtvol==-level)then
147
write(message,'(80a,a,a)') ('=',ii=1,80),ch10,&
148
& ' suscep : enter , debug mode '
149
call wrtout(06,message,'COLL')
152
!Loop input variables
153
nfreqsus=dtset%nfreqsus
155
dielar(1)=dtset%diecut
156
dielar(2)=dtset%dielng
157
dielar(3)=dtset%diemac
158
dielar(4)=dtset%diemix
159
dielar(5)=dtset%diegap
160
dielar(6)=dtset%dielam
162
!Unit numbers and file name for the _KGS file
163
filkgs=trim(dtfil%filnam_ds(5))//'_KGS'
164
filsustr=trim(dtfil%filnam_ds(5))//'_SUSTR'
166
!Impose mkmem=0 to read cg from disk
169
call status(0,dtfil%filstat,iexit,level,'call inwffil3 ')
171
!Compute different geometric tensor, as well as ucvol, from rprimd
172
call mkrdim(dtset%acell_orig,dtset%rprim_orig,rprimd)
173
call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
175
!Get diecut, and the fft grid to be used for the susceptibility computation
176
diecut=abs(dielar(1))
177
if( dielar(1)<0.0_dp )then
180
ecutsus= ( sqrt( dtset%ecut) *0.5_dp + sqrt(diecut) *0.25_dp )**2
183
ngfftdiel(1:3)=0 ; ngfftdiel(7)=101 ; ngfftdiel(8:18)=dtset%ngfft(8:18)
184
call getng(dtset%boxcutmin,ecutsus,gmet,mpi_enreg%me_fft,mgfftdiel,nfftdiel,ngfftdiel,&
185
& mpi_enreg%nproc_fft,nsym,mpi_enreg%fft_option_lob,mpi_enreg%paral_fft,dtset%symrel)
187
!Compute the size of the dielectric matrix
188
kpt_diel(1:3)=(/ 0.0_dp, 0.0_dp, 0.0_dp /)
189
call getmpw(diecut,dtset%exchn2n3,gmet,(/1/),kpt_diel,&
190
& mpi_enreg,npwdiel,1,ucvol)
192
!Now, performs allocation
193
allocate(cg(2,mpw*nspinor*mband*mkmem*nsppol))
194
allocate(eigen(mband*nkpt*nsppol))
195
allocate(kg(3,mpw*mkmem))
196
allocate(kg_diel(3,npwdiel))
197
allocate(npwarr(nkpt),npwtot(nkpt))
198
allocate(gbound_diel(2*mgfftdiel+8,2))
199
allocate(irrzondiel(nfftdiel**(1-1/nsym),2,nspden/nsppol))
200
allocate(phnonsdiel(2,nfftdiel**(1-1/nsym),nspden/nsppol))
201
allocate(nband_dum(nsppol))
203
!Then, initialize and compute the values of different arrays
204
call kpgio(dtset%ecut,dtset%exchn2n3,gmet,dtset%istwfk,kg,filkgs,dtset%kptns,mkmem,&
205
& dtset%nband,nkpt,'PERS',mpi_enreg,mpw,npwarr,npwtot,nsppol,dtfil%unkg)
206
!This kpgio call for going from the suscep FFT grid to the diel sphere
207
!Note : kgnam is dummy, npwarr_diel is dummy, npwtot_diel is dummy, nband_dum is dummy
209
call kpgio(diecut,dtset%exchn2n3,gmet,(/1/),kg_diel,kgnam,&
210
& kpt_diel,1,nband_dum,1,'COLL',mpi_enreg,npwdiel,npwarr_diel,npwtot_diel,&
213
call sphereboundary(gbound_diel,1,kg_diel,mgfftdiel,npwdiel)
215
allocate(indsym(4,nsym,natom),symrec(3,3,nsym))
217
call setsym(densymop_diel,indsym,irrzondiel,dtset%iscf,natom,&
218
& nfftdiel,ngfftdiel,nspden,nsppol,nsym,phnonsdiel,&
219
& dtset%symafm,symrec,dtset%symrel,dtset%tnons,dtset%typat,xred)
222
!Read eigenvalues from the wavefunction file
223
!Also, initialize wff1 and hdr
225
! mpi_enreg%paralbd=0
226
call inwffil3(dtset,eigen,hdr,dtset%istwfk,mband,mpi_enreg,dtset%nband,&
227
& nkpt,npwarr,nsppol,prtvol,wff1,dtfil%unwff1,dtfil%fnamewffk)
229
!Compute new occupation numbers if needed
230
allocate(doccde(mband*nkpt*nsppol))
231
if(dtset%occopt>=3.and.dtset%occopt<=7) then
232
call status(0,dtfil%filstat,iexit,level,'call newocc ')
233
call newocc(doccde,eigen,entropy,fermie,dtset%fixmom,mband,dtset%nband,&
234
& dtset%nelect,nkpt,nspinor,nsppol,occ,dtset%occopt,prtvol,dtset%stmbias,&
235
& dtset%tphysel,dtset%tsmear,dtset%wtk)
238
dielop=2 ! Immediate computation of dielectric matrix
241
! Perform allocations
242
allocate(dielinv(2,npwdiel,nspden,npwdiel,nspden))
243
allocate(susmat(2,npwdiel,nspden,npwdiel,nspden))
244
susmat(:,:,:,:,:)=0._dp
246
! Compute the static susceptibility matrix
247
call suscep_stat(cg,densymop_diel,&
248
& dielar,dielop,doccde,eigen,gbound_diel,&
249
& irrzondiel,dtset%istwfk,kg,kg_diel,&
250
& mband,mgfftdiel,mkmem,mpi_enreg,mpw,dtset%nband,nfftdiel,ngfftdiel,&
251
& nkpt,npwarr,npwdiel,nspden,nspinor,nsppol,nsym,&
252
& occ,dtset%occopt,phnonsdiel,prtvol,&
253
& susmat,dtset%symafm,dtset%symrel,dtset%tnons,ucvol,&
254
& dtfil%unkg,wff1,dtset%wtk)
256
! Print the susceptibility matrix
259
write(6,'(5x,a,2i2)') 'Susceptibility matrix for spins=',isp1,isp2
260
write(6,'(9x,a,13x,a,10x,a,10x,a)') "g","g'","real","imag"
263
write(6,'(2x,3i4,2x,3i4,2x,f12.8,2x,f12.8)') &
264
& kg_diel(1:3,ipw1),kg_diel(1:3,ipw2),&
265
& susmat(1,ipw1,isp1,ipw2,isp2),susmat(2,ipw1,isp1,ipw2,isp2)
271
! Perform deallocations
274
deallocate(nband_dum)
276
else if(nfreqsus > 0) then
278
! Perform allocations
279
allocate(freq(nfreqsus))
280
allocate(rhor(nfft,nspden))
281
allocate(wght_freq(nfreqsus))
283
! Perform initializations
288
! Read in the density, needed for ALDA kernel
289
call status(0,dtfil%filstat,iexit,level,'call ioarr ')
290
! Read rho(r) from a disk file
291
! Unit numbers and file name for the _KGS file
293
! Note : etotal is read here, and might serve in the tddft routine.
295
call ioarr(accessfil,rhor,etotal,fformr,dtfil%fildensin,hdr,&
296
& nfft,nspden,rdwr,rdwrpaw,ngfft)
297
call status(0,dtfil%filstat,iexit,level,'call fourdp ')
300
!leave in for MF to check the density
303
! dummy=dummy+rhor(ii,1)
305
! write(6,*) '%suscep: nfft=',nfft
306
! write(6,*) '%suscep: dummy=',dummy*ucvol/dble(nfft)
307
! write(6,*) '%suscep: ucvol=',ucvol
309
! Compute up+down rho(G) by fft
310
! allocate(work(nfft))
312
! call fourdp(1,rhog,work,-1,mpi_enreg,nfft,ngfft,0)
316
! Create frequency grid and weights, 2 stands for preassigned grid
317
call getfreqsus(freq,wght_freq,nfreqsus,dtset%optfreqsus,dtset%freqsuslo,dtset%freqsusin)
320
! write(6,*)' suscep : after getfreqsus '
324
call xcacfd(cg,densymop_diel,dielar,dielop,doccde,&
325
& dtfil,dtset,eigen,filsustr,freq,gbound_diel,gmet,&
326
& gprimd,irrzondiel,kg,kg_diel,mband,mgfftdiel,mkmem,&
327
& mpi_enreg,mpw,nfft,nfftdiel,nfreqsus,dtset%ngfft,&
328
& ngfftdiel,nkpt,npwarr,npwdiel,nspden,nspinor,nsppol,&
329
& nsym,occ,phnonsdiel,rhor,rmet,rprimd,ucvol,wff1,wght_freq,ncdims)
331
! Perform deallocations
334
deallocate(wght_freq)
338
! Leave intact for testing purposes
339
nfreqsus=abs(nfreqsus)
341
! Perform allocations
342
allocate(freq(nfreqsus))
343
allocate(susd_non_dyn(2,npwdiel,nspden,nfreqsus))
344
allocate(susmat_dyn(2,npwdiel,nspden,npwdiel,nspden,nfreqsus))
346
! Perform initializations
348
susmat_dyn(:,:,:,:,:,:)=0.0_dp
350
! Create a linear frequency grid
351
freq(1)=dtset%freqsuslo
353
freq(ifreq)=freq(ifreq-1)+dtset%freqsusin
357
! write(6,*) '%suscep: nband_mx=', nband_mx
360
! Compute the dynamical susceptibility matrices
361
call suscep_dyn(cg,densymop_diel,dielar,dielop,doccde,dtset,&
362
& eigen,freq,gbound_diel,&
363
& irrzondiel,dtset%istwfk,kg,kg_diel,&
364
& mband,mgfftdiel,mkmem,mpi_enreg,mpw,dtset%nband,nband_mx,nfftdiel,nfreqsus,&
365
& ngfftdiel,nkpt,npwarr,npwdiel,nspden,nspinor,nsppol,nsym,&
366
& occ,dtset%occopt,phnonsdiel,prtvol,&
367
& susopt,susd_non_dyn,susmat_dyn,dtset%symafm,&
368
& dtset%symrel,dtset%tnons,ucvol,dtfil%unkg,wff1,dtset%wtk)
370
! Print the dynamical susceptibility matrices
372
write(6,'(/,2x,a,f12.8,a)') '---Susceptibility matrices for frequency=',freq(ifreq),'i'
375
write(6,'(5x,a,2i2)') 'Susceptibility matrix for spins=',isp1,isp2
376
write(6,'(9x,a,13x,a,10x,a,10x,a)') "g","g'","real","imag"
379
write(6,'(2x,3i4,2x,3i4,2x,f12.8,2x,f12.8)') &
380
& kg_diel(1:3,ipw1),kg_diel(1:3,ipw2),&
381
& susmat_dyn(1,ipw1,isp1,ipw2,isp2,ifreq),susmat_dyn(2,ipw1,isp1,ipw2,isp2,ifreq)
388
! Perform deallocations
390
deallocate(susd_non_dyn)
391
deallocate(susmat_dyn)
393
end if !condition nfreqsus
395
!Performs deallocations
396
deallocate(cg,doccde,eigen,gbound_diel)
397
deallocate(indsym,irrzondiel)
398
deallocate(kg_diel,kg,npwarr,npwtot,phnonsdiel,symrec)
402
if(mpi_enreg%paral_compil_kpt==0)then
403
! Unit dtfil%unkg was opened in kpgio
404
close (unit=dtfil%unkg,status='delete')
407
else if(mpi_enreg%paral_compil_kpt==1)then
409
! All procs close the file dtfil%unkg
410
close(unit=dtfil%unkg)
411
if(mpi_enreg%me==0)then
412
! only proc 0 delete the file
413
open(unit=dtfil%unkg,file=filkgs,form='unformatted',status='unknown')
414
close(unit=dtfil%unkg,status='delete')
420
call WffClose(wff1,ierr)
425
if(mpi_enreg%paral_compil_kpt==1)then
426
deallocate(mpi_enreg%proc_distrb)
429
write(message, '(a,a)' ) ch10,' suscep : exiting '
430
call wrtout(06,message,'COLL')
432
call status(0,dtfil%filstat,iexit,level,'exit ')
433
call timab(84,2,tsec)
435
end subroutine suscep
added
removed
src/17suscep/suscep.F90
