Branch data 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%atmfrc,Ifc%dielt,Ifc%dipdip,Ifc%dyewq0,d2cart, &
208 : : crystal%gmet,Ifc%gprimd,Ifc%mpert,natom,Ifc%nrpt,qptnrm(1), &
209 : : qpt, crystal%rmet,Ifc%rprimd,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 : : !!***
|
