Line data Source code
1 : !!****m* ABINIT/m_ddb_interpolate
2 : !! NAME
3 : !! m_ddb_interpolate
4 : !!
5 : !! FUNCTION
6 : !! Interpolate the ddb onto a fine set of q-points using
7 : !! the interatomic force constants and write the result in a DDB file.
8 : !!
9 : !! COPYRIGHT
10 : !! Copyright (C) 2008-2026 ABINIT group (GA)
11 : !! This file is distributed under the terms of the
12 : !! GNU General Public License, see ~abinit/COPYING
13 : !! or http://www.gnu.org/copyleft/gpl.txt .
14 : !!
15 : !! SOURCE
16 :
17 : #if defined HAVE_CONFIG_H
18 : #include "config.h"
19 : #endif
20 :
21 : #include "abi_common.h"
22 :
23 : module m_ddb_interpolate
24 :
25 : use defs_basis
26 : use m_errors
27 : use m_xmpi
28 : use m_abicore
29 : use m_ddb
30 : use m_ddb_hdr
31 : use m_ifc
32 : use m_nctk
33 : use netcdf
34 :
35 : use m_anaddb_dataset, only : anaddb_dataset_type
36 : use m_bz_mesh, only : make_path
37 : use m_crystal, only : crystal_t
38 : use m_io_tools, only : get_unit
39 : use m_fstrings, only : strcat
40 : use m_dynmat, only : gtdyn9, d2cart_to_red
41 :
42 : implicit none
43 :
44 : private
45 : !!***
46 :
47 : public :: ddb_interpolate
48 : !!***
49 :
50 : contains
51 : !!***
52 :
53 : !!****f* ABINIT/ddb_interpolate
54 : !!
55 : !! NAME
56 : !! ddb_interpolate
57 : !!
58 : !! FUNCTION
59 : !! Interpolate the ddb onto a fine set of q-points using
60 : !! the interatomic force constants and write the result in a DDB file.
61 : !!
62 : !! INPUTS
63 : !!
64 : !! OUTPUT
65 : !!
66 : !! NOTES
67 : !!
68 : !! SOURCE
69 :
70 1 : subroutine ddb_interpolate(ifc, crystal, inp, ddb, ddb_hdr, asrq0, comm)
71 :
72 : !Arguments -------------------------------
73 : !scalars
74 : type(ifc_type),intent(in) :: ifc
75 : type(anaddb_dataset_type),target,intent(inout) :: inp
76 : type(crystal_t),intent(in) :: crystal
77 : type(ddb_type),intent(inout) :: ddb
78 : type(ddb_hdr_type),intent(inout) :: ddb_hdr
79 : type(asrq0_t),intent(inout) :: asrq0
80 : integer,intent(in) :: comm
81 : !arrays
82 :
83 : !Local variables -------------------------
84 : !scalars
85 : integer,parameter :: master=0
86 : integer :: nsym,natom,ntypat,mband,nqpt_fine
87 : integer :: msize,nsize,mpert,nblok
88 : integer :: rftyp
89 : integer :: ii,iblok,jblok,iqpt,ipert1,ipert2,idir1,idir2
90 : integer :: nprocs,my_rank
91 : real(dp) :: eta
92 : character(len=500) :: msg
93 : character(len=fnlen) :: ddb_out_filename, ddb_out_nc_filename
94 1 : type(ddb_type) :: ddb_new
95 : !arrays
96 : integer :: rfphon(4),rfelfd(4),rfstrs(4)
97 1 : integer,allocatable :: blkflg(:,:,:,:)
98 : real(dp) :: qpt(3), qptnrm(3), qpt_padded(3,3)
99 1 : real(dp),allocatable :: d2cart(:,:,:,:,:),d2red(:,:,:,:,:)
100 1 : real(dp),pointer :: qpt_fine(:,:)
101 1 : integer,allocatable :: ndiv(:)
102 1 : real(dp),allocatable,target :: alloc_path(:,:)
103 :
104 : ! *********************************************************************
105 :
106 :
107 : ! Only master works for the time being
108 1 : nprocs = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
109 1 : if (my_rank /= master) return
110 :
111 : ! ===================================
112 : ! Copy dimensions and allocate arrays
113 : ! ===================================
114 :
115 81 : write(msg, '(a,(80a),a,a,a,a)' ) ch10,('=',ii=1,80),ch10,ch10,' Treat the first list of vectors ',ch10
116 1 : call wrtout(std_out,msg,'COLL')
117 1 : call wrtout(ab_out,msg,'COLL')
118 :
119 1 : nullify(qpt_fine)
120 1 : nqpt_fine = inp%nph1l
121 1 : qpt_fine => inp%qph1l
122 :
123 1 : if(inp%nph1l==0) then
124 0 : if (inp%nqpath==0) then
125 : return ! if there is nothing to do, return
126 : else
127 : ! allow override of nph1l with nqpath if the former is not set
128 0 : ABI_MALLOC(ndiv,(inp%nqpath-1))
129 0 : call make_path(inp%nqpath,inp%qpath,Crystal%gmet,'G',inp%ndivsm,ndiv,nqpt_fine,alloc_path,std_out)
130 0 : ABI_FREE(ndiv)
131 0 : qpt_fine => alloc_path
132 : end if
133 : end if
134 :
135 1 : rftyp=inp%rfmeth
136 :
137 1 : nsym = Crystal%nsym
138 1 : natom = Crystal%natom
139 1 : ntypat = Crystal%ntypat
140 :
141 1 : mband = ddb_hdr%mband
142 :
143 : ! Interpolation is limited to 2nd derivatives of total energy
144 : ! GA: What??
145 1 : ddb_hdr%has_d3E_xx = .false.
146 1 : ddb_hdr%has_d3E_lw = .false.
147 1 : ddb_hdr%has_d2eig = .false.
148 :
149 1 : mpert = ddb%mpert
150 1 : msize = 3 * mpert * 3 * mpert !; if (ddb_hdr%has_d3E_xx) msize=msize*3*mpert
151 1 : nsize = 3 * mpert * 3 * mpert
152 1 : nblok = nqpt_fine
153 :
154 1 : ddb_new%nblok = nblok
155 1 : call ddb_new%malloc(msize,nblok,natom,ntypat,mpert)
156 16209 : ddb_new%flg = 0
157 3 : ddb_new%amu = ddb%amu
158 1 : if (rftyp == 1 .or. rftyp == 2) then
159 9 : ddb_new%typ = 1
160 0 : else if (rftyp == 85) then
161 0 : ddb_new%typ = 85
162 : end if
163 81 : ddb_new%qpt = zero
164 33 : ddb_new%omega = zero
165 33 : ddb_new%nrm = one
166 :
167 5 : ABI_MALLOC(d2cart,(2,3,mpert,3,mpert))
168 3 : ABI_MALLOC(d2red,(2,3,mpert,3,mpert))
169 5 : ABI_MALLOC(blkflg,(3,mpert,3,mpert))
170 :
171 2761 : blkflg = 1
172 :
173 7 : rfphon(1:2)=1; rfelfd(1:2)=0; rfstrs(1:2)=0
174 1 : qpt_padded = zero
175 :
176 :
177 1 : ddb_hdr%dscrpt = 'Interpolated DDB using interatomic force constants'
178 1 : ddb_hdr%nblok = nblok
179 :
180 : ! ================================================
181 : ! Interpolate the dynamical matrix at each q-point
182 : ! ================================================
183 :
184 4 : qptnrm = one
185 :
186 9 : do iqpt=1,nqpt_fine
187 :
188 : ! Initialisation of the phonon wavevector
189 32 : qpt(:)=qpt_fine(:,iqpt)
190 :
191 8 : if (inp%nph1l /= 0) qptnrm(1) = inp%qnrml1(iqpt)
192 :
193 : ! Look for the information in the DDB
194 32 : qpt_padded(:,1) = qpt
195 8 : call ddb%get_block(iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
196 :
197 8 : if (iblok /= 0) then
198 :
199 : ! q-point is present in the ddb. No interpolation needed.
200 :
201 8295 : d2cart(1,:,:,:,:) = reshape(ddb%val(1,:,iblok), shape = (/3,mpert,3,mpert/))
202 8295 : d2cart(2,:,:,:,:) = reshape(ddb%val(2,:,iblok), shape = (/3,mpert,3,mpert/))
203 : else
204 :
205 : ! Get d2cart using the interatomic forces and the
206 : ! long-range coulomb interaction through Ewald summation
207 : call gtdyn9(Ifc%acell,Ifc%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart, &
208 : crystal%gmet,Ifc%gprim,Ifc%mpert,natom,Ifc%nrpt,qptnrm(1), &
209 : qpt, crystal%rmet,Ifc%rprim,Ifc%rpt,Ifc%trans,crystal%ucvol, &
210 : Ifc%wghatm,crystal%xred,ifc%zeff,ifc%qdrp_cart,ifc%ewald_option,eta, xmpi_comm_self, &
211 : inp%sys_dim,dipquad=Ifc%dipquad,quadquad=Ifc%quadquad,&
212 5 : dielt_thick=inp%dielt_thick,dielt_env=inp%dielt_env)
213 :
214 : end if
215 :
216 : ! Eventually impose the acoustic sum rule based on previously calculated d2asr
217 8 : call asrq0%apply(natom, ddb%mpert, ddb%msize, qpt_padded, crystal, d2cart)
218 :
219 : ! Transform d2cart into reduced coordinates.
220 : call d2cart_to_red(d2cart,d2red,crystal%gprimd,crystal%rprimd,mpert, &
221 8 : & natom,ntypat,crystal%typat,crystal%ucvol,crystal%zion)
222 :
223 : ! TODO (GA): Should replace this with ddb_set_d2matr
224 : ! Store the dynamical matrix into a block of the new ddb
225 8 : jblok = iqpt
226 16216 : ddb_new%val(1,1:nsize,jblok) = reshape(d2red(1,:,:,:,:), shape = (/nsize/))
227 16216 : ddb_new%val(2,1:nsize,jblok) = reshape(d2red(2,:,:,:,:), shape = (/nsize/))
228 :
229 : ! Store the q-point
230 32 : ddb_new%qpt(1:3,jblok) = qpt
231 8 : ddb_new%nrm(1,jblok) = qptnrm(1)
232 :
233 : ! Set the flags
234 8 : ii=0
235 129 : do ipert2=1,mpert
236 488 : do idir2=1,3
237 5880 : do ipert1=1,mpert
238 21960 : do idir1=1,3
239 16200 : ii=ii+1
240 21600 : if (ipert1<=natom.and.ipert2<=natom) then
241 288 : ddb_new%flg(ii,jblok) = 1
242 : end if
243 : end do
244 : end do
245 : end do
246 : end do
247 :
248 : end do ! iqpt
249 :
250 : ! Copy the flags for Gamma
251 4 : qpt_padded(:,1) = zero
252 4 : qptnrm = one
253 1 : call ddb%get_block(iblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
254 1 : call ddb_new%get_block(jblok,qpt_padded,qptnrm,rfphon,rfelfd,rfstrs,rftyp)
255 :
256 1 : if (iblok /= 0 .and. jblok /= 0) then
257 :
258 : ii=0
259 16 : do ipert2=1,mpert
260 61 : do idir2=1,3
261 735 : do ipert1=1,mpert
262 2745 : do idir1=1,3
263 2025 : ii=ii+1
264 2700 : ddb_new%flg(ii,jblok) = ddb%flg(ii,iblok)
265 : end do
266 : end do
267 : end do
268 : end do
269 :
270 : end if
271 :
272 :
273 : ! =======================
274 : ! Write the new DDB files
275 : ! =======================
276 :
277 : if (my_rank == master) then
278 :
279 : ! GA: TODO choice of txt vs. nc should be set by user
280 1 : ddb_out_filename = strcat(inp%prefix_outdata, "_DDB")
281 :
282 1 : call ddb_new%write_txt(ddb_hdr, ddb_out_filename)
283 :
284 1 : ddb_out_nc_filename = strcat(inp%prefix_outdata, "_DDB.nc")
285 1 : call ddb_new%write_nc(ddb_hdr, ddb_out_nc_filename)
286 :
287 : ! Write one separate nc file for each q-point
288 : !do jblok=1,nblok
289 : ! write(ddb_out_nc_filename,'(2a,i5.5,a)') trim(prefix),'_qpt_',jblok,'_DDB.nc'
290 : ! call ddb_new%write_nc(ddb_hdr, ddb_out_nc_filename, jblok)
291 : !end do
292 :
293 : end if
294 :
295 : ! ===========
296 : ! Free memory
297 : ! ===========
298 :
299 1 : call ddb_new%free()
300 1 : ABI_FREE(d2cart)
301 1 : ABI_FREE(d2red)
302 1 : ABI_FREE(blkflg)
303 :
304 2 : end subroutine ddb_interpolate
305 : !!***
306 :
307 : end module m_ddb_interpolate
308 : !!***
|
