Line data Source code
1 : !!****m*ABINIT/m_lwf_primitive_potential
2 : !! NAME
3 : !! m_lwf_primitive_potential
4 : !!
5 : !! FUNCTION
6 : !! This module contains the example lwf primitive potential. It is not the
7 : !! one really used. The purpose of this is to show how to extend multibinit
8 : !! to new degree of freedom
9 : !!
10 : !! Datatypes:
11 : !! LWF_primitive_potential_t
12 : !!
13 : !! Subroutines:
14 : !!
15 : !! COPYRIGHT
16 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
17 : !! This file is distributed under the terms of the
18 : !! GNU General Public License, see ~abinit/COPYING
19 : !! or http://www.gnu.org/copyleft/gpl.txt .
20 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
21 : !!
22 : !! SOURCE
23 :
24 :
25 : #if defined HAVE_CONFIG_H
26 : #include "config.h"
27 : #endif
28 : #include "abi_common.h"
29 :
30 : #define HAVE_NETCDF 1
31 :
32 : module m_lwf_primitive_potential
33 : use, intrinsic :: iso_c_binding
34 : use defs_basis
35 : use m_abicore
36 : use m_errors
37 : use m_nctk
38 : #if defined HAVE_NETCDF
39 : use netcdf
40 : #endif
41 : use m_xmpi
42 : use m_mathfuncs, only: eigensh
43 : use m_multibinit_dataset, only: multibinit_dtset_type
44 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
45 : use m_primitive_potential, only: primitive_potential_t
46 : use m_abstract_potential, only: abstract_potential_t
47 : use m_dynamic_array, only: int2d_array_type
48 : use m_supercell_maker, only: supercell_maker_t
49 : use m_spmat_ndcoo, only: ndcoo_mat_t
50 : use m_lwf_potential, only: lwf_potential_t
51 : implicit none
52 : private
53 : !!***
54 :
55 : !-------------------------------------------------------------------!
56 : ! An harmonic potential which has only the IFC (with no dipole-dipole)
57 : !
58 : ! IFC is written in a coefficient matrix M(R, i, j)=val
59 : !-------------------------------------------------------------------!
60 : type, public, extends(primitive_potential_t):: lwf_primitive_potential_t
61 : integer:: nlwf = 0, natom = 0, nR = 0 ! number of LWF
62 : type(ndcoo_mat_t):: coeff ! A N-dimensional COO matrix.
63 : ! The indices are (ind_R, i, j). Note that xyz is included in i and j.
64 : real(dp), allocatable:: lwf_masses(:)
65 :
66 : integer:: onebody_nterm = 0
67 :
68 : real(dp), allocatable:: lattice_coeffs(:,:, :) ! (natom*3, nlwf, nR)
69 :
70 : integer, allocatable:: onebody_i(:), onebody_order(:)
71 : real(dp), allocatable:: onebody_val(:)
72 :
73 : type(ndcoo_mat_t):: coeff_twobody ! Two body interaction parameters of higher order.
74 : ! The indices are (indR, i, j, orderi, orderj), H(R, i, j)= val A_i^(orderi) A_j^(orderj)
75 : integer, allocatable:: Rlist(:,:) ! The list of R points (3, number of R-points)
76 :
77 : ! For testing the bounding.
78 : !logical:: has_self_bound_term= .False.
79 : !integer:: self_bound_order = 0
80 : !real(dp):: self_bound_coeff = 0.0_dp
81 :
82 : !. Rlist(:, ind_R) is a R-vector.
83 : real(dp):: ref_energy = 0.0 ! reference energy
84 :
85 : logical:: as_lattice_anharmonic=.False.
86 :
87 : contains
88 : procedure:: initialize
89 : procedure:: finalize
90 : procedure:: load_from_files ! load potential from files listed in files file
91 : procedure:: load_from_netcdf ! load potential from a netcdf file
92 : procedure:: fill_supercell ! fill a supercell potential
93 : procedure:: get_hamk ! generate hamiltonian for one k point.
94 : procedure:: get_eigen ! eigen values and eigen vectors
95 : !procedure:: add_self_bound_term ! add a self bound term
96 : end type lwf_primitive_potential_t
97 :
98 :
99 : contains
100 :
101 : !-------------------------------------------------------------------!
102 : ! Initialize
103 : ! Input:
104 : ! nlwf: number of Wannier functions
105 : ! primcell: the reference primitive cell.
106 : !-------------------------------------------------------------------!
107 1 : subroutine initialize(self, primcell)
108 : class(lwf_primitive_potential_t), intent(inout):: self
109 : type(mbcell_t), target, intent(inout):: primcell
110 : !integer, intent(in):: nlwf
111 1 : self%primcell=>primcell
112 1 : self%label="lwf_primitive_potential"
113 1 : self%has_spin=.False.
114 1 : self%has_displacement=.False.
115 1 : self%has_strain=.False.
116 1 : self%has_lwf=.True.
117 1 : end subroutine initialize
118 :
119 :
120 : !-------------------------------------------------------------------!
121 : ! Finalize
122 : !-------------------------------------------------------------------!
123 1 : subroutine finalize(self)
124 : class(lwf_primitive_potential_t), intent(inout):: self
125 1 : call self%coeff%finalize()
126 1 : call self%coeff_twobody%finalize()
127 1 : ABI_SFREE(self%Rlist)
128 1 : nullify(self%primcell)
129 1 : self%nlwf = 0
130 1 : self%label="Destroyed lwf_primitive_potential"
131 1 : call self%primitive_potential_t%finalize()
132 :
133 1 : ABI_SFREE(self%lattice_coeffs)
134 1 : ABI_SFREE(self%lwf_masses)
135 1 : self%natom = 0
136 1 : self%nR = 0
137 :
138 1 : if (self%onebody_nterm /= 0) then
139 1 : ABI_SFREE(self%onebody_i)
140 1 : ABI_SFREE(self%onebody_val)
141 1 : ABI_SFREE(self%onebody_order)
142 : end if
143 1 : self%onebody_nterm = 0
144 1 : end subroutine finalize
145 :
146 :
147 : !-------------------------------------------------------------------!
148 : ! Load from files:
149 : !> params: parameters
150 : !> fnames: file names from files file
151 : !-------------------------------------------------------------------!
152 1 : subroutine load_from_files(self, params, fnames)
153 : class(lwf_primitive_potential_t), intent(inout):: self
154 : type(multibinit_dtset_type), intent(in):: params
155 : character(len = fnlen), intent(in):: fnames(:)
156 1 : call self%load_from_netcdf(fnames(1))
157 : !call self%add_self_bound_term(params%lwf_self_bound_order, &
158 : ! & params%lwf_self_bound_coeff)
159 1 : self%as_lattice_anharmonic= (params%latt_lwf_anharmonic == 1)
160 1 : end subroutine load_from_files
161 :
162 : !-------------------------------------------------------------------!
163 : ! load potential from netcdf file
164 : ! Note that the lattic part of the primitive cell is also loaded.
165 : ! Input:
166 : ! fname: filename
167 : !-------------------------------------------------------------------!
168 1 : subroutine load_from_netcdf(self, fname)
169 : class(lwf_primitive_potential_t), intent(inout):: self
170 : character(len = fnlen), intent(in):: fname
171 : integer:: ncid, ierr
172 : integer:: iR, nR, nlwf, natom, twobody_nterm, onebody_nterm
173 : real(dp):: cell(3, 3)
174 1 : real(dp), allocatable:: ifc(:, :, :), xcart(:,:), masses(:)
175 1 : real(dp), allocatable :: twobody_val(:)
176 1 : integer, allocatable :: zion(:), twobody_iR(:), twobody_i(:), twobody_j(:), twobody_orderi(:), twobody_orderj(:)
177 : integer:: varid, i, j
178 : #if defined HAVE_NETCDF
179 1 : ierr = nf90_open(trim(fname), NF90_NOWRITE, ncid)
180 1 : NCF_CHECK_MSG(ierr, "Open netcdf file")
181 :
182 : ! read primcell info
183 1 : ierr = nctk_get_dim(ncid, "wann_natom", natom)
184 1 : NCF_CHECK_MSG(ierr, "getting natom in lwf potential file")
185 :
186 1 : ierr = nctk_get_dim(ncid, "nR" , nR)
187 1 : NCF_CHECK_MSG(ierr, "getting nR in lwf potential file")
188 1 : self%nR = nR
189 :
190 1 : ierr = nctk_get_dim(ncid, "wann_nwann", nlwf)
191 1 : NCF_CHECK_MSG(ierr, "getting wann_nwann in lwf potential file")
192 :
193 : ! TODO : add lattice
194 : !call self%primcell%set_lattice()
195 :
196 1 : call self%primcell%set_lwf(nlwf)
197 1 : self%nlwf = nlwf
198 4 : call self%coeff%initialize([ nR, nlwf, nlwf])
199 6 : call self%coeff_twobody%initialize([nR, nlwf, nlwf, -1, -1])
200 :
201 5 : ABI_MALLOC(ifc, (nlwf, nlwf, nR))
202 3 : ABI_MALLOC(self%Rlist, (3, nR))
203 3 : ABI_MALLOC(xcart, (3, natom))
204 3 : ABI_MALLOC(zion, (natom))
205 3 : ABI_MALLOC(masses, (natom))
206 :
207 3 : ABI_MALLOC(self%lwf_masses, (nlwf))
208 :
209 1 : ierr = nf90_inq_varid(ncid, "wann_lwf_masses", varid)
210 1 : NCF_CHECK_MSG(ierr, "lwf_masses")
211 1 : ierr = nf90_get_var(ncid, varid, self%lwf_masses)
212 1 : NCF_CHECK_MSG(ierr, "lwf_masses")
213 3 : self%lwf_masses(:) = self%lwf_masses(:) * amu_emass
214 4 : self%primcell%lwf%lwf_masses = self%lwf_masses
215 :
216 1 : ierr = nf90_inq_dimid(ncid, "wann_onebody_nterm", onebody_nterm)
217 1 : if (ierr /= nf90_noerr) then
218 0 : onebody_nterm = 0
219 : else
220 1 : ierr = nctk_get_dim(ncid, "wann_onebody_nterm", onebody_nterm)
221 1 : NCF_CHECK_MSG(ierr, "getting wann_onebody_nterm in lwf potential file")
222 : end if
223 1 : self%onebody_nterm = onebody_nterm
224 :
225 :
226 1 : ierr = nf90_inq_dimid(ncid, "wann_twobody_nterm", twobody_nterm)
227 1 : if (ierr /= nf90_noerr) then
228 1 : twobody_nterm = 0
229 : else
230 0 : ierr = nctk_get_dim(ncid, "wann_twobody_nterm", twobody_nterm)
231 0 : NCF_CHECK_MSG(ierr, "getting wann_twobody_nterm in lwf potential file")
232 : end if
233 :
234 1 : self%natom = natom
235 :
236 :
237 1 : ierr = nf90_inq_varid(ncid, "wann_cell", varid)
238 1 : NCF_CHECK_MSG(ierr, "wann_cell")
239 1 : ierr = nf90_get_var(ncid, varid, cell)
240 1 : NCF_CHECK_MSG(ierr, "wann_cell")
241 13 : cell(:,:)=cell(:,:)/ Bohr_Ang
242 :
243 :
244 1 : ierr = nf90_inq_varid(ncid, "wann_atomic_xcart", varid)
245 1 : NCF_CHECK_MSG(ierr, "wann_atomic_xcart")
246 1 : ierr = nf90_get_var(ncid, varid, xcart)
247 1 : NCF_CHECK_MSG(ierr, "wann_atomic_xcart")
248 25 : xcart(:,:)=xcart(:,:)/ Bohr_Ang
249 :
250 1 : ierr = nf90_inq_varid(ncid, "wann_atomic_numbers", varid)
251 1 : NCF_CHECK_MSG(ierr, "wann_atomic_numbers")
252 1 : ierr = nf90_get_var(ncid, varid, zion)
253 1 : NCF_CHECK_MSG(ierr, "wann_atomic_numbers")
254 :
255 1 : ierr = nf90_inq_varid(ncid, "wann_atomic_masses", varid)
256 1 : NCF_CHECK_MSG(ierr, "wann_atomic_masses")
257 1 : ierr = nf90_get_var(ncid, varid, masses)
258 1 : NCF_CHECK_MSG(ierr, "wann_atomic_masses")
259 :
260 1 : call self%primcell%set_lattice(natom, cell, xcart, masses, zion)
261 :
262 :
263 5 : ABI_MALLOC(self%lattice_coeffs, (nlwf, natom*3, nR))
264 :
265 1 : ierr = nf90_inq_varid(ncid, "wann_wannier_function_real", varid)
266 1 : NCF_CHECK_MSG(ierr, "wann_wannier_function_real")
267 1 : ierr = nf90_get_var(ncid, varid, self%lattice_coeffs)
268 1 : NCF_CHECK_MSG(ierr, "wann_wannier_function_real")
269 :
270 1 : ierr = nf90_inq_varid(ncid, "wann_Rlist", varid)
271 1 : NCF_CHECK_MSG(ierr, "wann_Rlist")
272 1 : ierr = nf90_get_var(ncid, varid, self%Rlist)
273 1 : NCF_CHECK_MSG(ierr, "wann_Rlist")
274 :
275 1 : ierr = nf90_inq_varid(ncid, "wann_HamR_real", varid)
276 1 : NCF_CHECK_MSG(ierr, "wann_HamR_real")
277 1 : ierr = nf90_get_var(ncid, varid, ifc)
278 1 : NCF_CHECK_MSG(ierr, "wann_HamR_real")
279 :
280 :
281 316 : ifc(:,:,:) = ifc(:,:,:) * eV_Ha * (Bohr_Ang*Bohr_Ang)
282 :
283 46 : do iR = 1, nR
284 136 : do i = 1, nlwf
285 315 : do j = 1, nlwf
286 270 : if (abs(ifc(j, i, iR))>1e-4) then
287 : ! NOTE: in fortran the order of index in reversed when reading netcdf array.
288 528 : call self%coeff%add_entry([iR, i, j], ifc(j, i, iR))
289 : end if
290 : end do
291 : end do
292 : end do
293 :
294 :
295 1 : if(self%onebody_nterm /= 0) then
296 3 : ABI_MALLOC(self%onebody_i, (self%onebody_nterm))
297 2 : ABI_MALLOC(self%onebody_order, (self%onebody_nterm))
298 3 : ABI_MALLOC(self%onebody_val, (self%onebody_nterm))
299 :
300 1 : ierr = nf90_inq_varid(ncid, "wann_onebody_i", varid)
301 1 : NCF_CHECK_MSG(ierr, "wann_onebody_i")
302 1 : ierr = nf90_get_var(ncid, varid, self%onebody_i)
303 1 : NCF_CHECK_MSG(ierr, "wann_onebody_i")
304 :
305 7 : self%onebody_i = self%onebody_i+1
306 :
307 1 : ierr = nf90_inq_varid(ncid, "wann_onebody_val", varid)
308 1 : NCF_CHECK_MSG(ierr, "wann_onebody_val")
309 1 : ierr = nf90_get_var(ncid, varid, self%onebody_val)
310 1 : NCF_CHECK_MSG(ierr, "wann_onebody_val")
311 :
312 1 : ierr = nf90_inq_varid(ncid, "wann_onebody_order", varid)
313 1 : NCF_CHECK_MSG(ierr, "wann_onebody_order")
314 1 : ierr = nf90_get_var(ncid, varid, self%onebody_order)
315 1 : NCF_CHECK_MSG(ierr, "wann_onebody_order")
316 :
317 7 : do i = 1, self%onebody_nterm
318 7 : self%onebody_val(i)=self%onebody_val(i) * eV_Ha * (Bohr_Ang**self%onebody_order(i))
319 : end do
320 : end if
321 :
322 1 : if(twobody_nterm /= 0) then
323 0 : ABI_MALLOC(twobody_iR, (twobody_nterm))
324 0 : ABI_MALLOC(twobody_i, (twobody_nterm))
325 0 : ABI_MALLOC(twobody_j, (twobody_nterm))
326 0 : ABI_MALLOC(twobody_orderi, (twobody_nterm))
327 0 : ABI_MALLOC(twobody_orderj, (twobody_nterm))
328 0 : ABI_MALLOC(twobody_val, (twobody_nterm))
329 :
330 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_iR", varid)
331 0 : NCF_CHECK_MSG(ierr, "wann_twobody_iR")
332 0 : ierr = nf90_get_var(ncid, varid, twobody_iR)
333 0 : NCF_CHECK_MSG(ierr, "wann_twobody_iR")
334 :
335 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_i", varid)
336 0 : NCF_CHECK_MSG(ierr, "wann_twobody_i")
337 0 : ierr = nf90_get_var(ncid, varid, twobody_i)
338 0 : NCF_CHECK_MSG(ierr, "wann_twobody_i")
339 :
340 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_j", varid)
341 0 : NCF_CHECK_MSG(ierr, "wann_twobody_j")
342 0 : ierr = nf90_get_var(ncid, varid, twobody_j)
343 0 : NCF_CHECK_MSG(ierr, "wann_twobody_j")
344 :
345 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_orderi", varid)
346 0 : NCF_CHECK_MSG(ierr, "wann_twobody_orderi")
347 0 : ierr = nf90_get_var(ncid, varid, twobody_orderi)
348 0 : NCF_CHECK_MSG(ierr, "wann_twobody_i")
349 :
350 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_orderj", varid)
351 0 : NCF_CHECK_MSG(ierr, "wann_twobody_orderj")
352 0 : ierr = nf90_get_var(ncid, varid, twobody_orderj)
353 0 : NCF_CHECK_MSG(ierr, "wann_twobody_orderj")
354 :
355 0 : ierr = nf90_inq_varid(ncid, "wann_twobody_val", varid)
356 0 : NCF_CHECK_MSG(ierr, "wann_twobody_val")
357 0 : ierr = nf90_get_var(ncid, varid, twobody_val)
358 0 : NCF_CHECK_MSG(ierr, "wann_twobody_val")
359 :
360 0 : ierr = nf90_close(ncid)
361 0 : NCF_CHECK_MSG(ierr, "Close netcdf file")
362 :
363 0 : do i = 1, twobody_nterm
364 : twobody_val(i)=twobody_val(i) * eV_Ha * &
365 0 : & (Bohr_Ang**(twobody_orderi(i)+ twobody_orderj(i)))
366 : call self%coeff_twobody%add_entry([twobody_iR(i), twobody_i(i), &
367 0 : & twobody_j(i), twobody_orderi(i), twobody_orderj(i)], twobody_val(i))
368 : end do
369 : end if
370 :
371 1 : ABI_SFREE(xcart)
372 1 : ABI_SFREE(zion)
373 1 : ABI_SFREE(masses)
374 1 : ABI_SFREE(ifc)
375 :
376 1 : if(twobody_nterm /= 0) then
377 0 : ABI_SFREE(twobody_i)
378 0 : ABI_SFREE(twobody_j)
379 0 : ABI_SFREE(twobody_orderi)
380 0 : ABI_SFREE(twobody_orderj)
381 0 : ABI_SFREE(twobody_val)
382 : endif
383 : #else
384 : NETCDF_NOTENABLED_ERROR()
385 : #endif
386 :
387 1 : end subroutine load_from_netcdf
388 :
389 : !-------------------------------------------------------------------!
390 : !Fill supercell
391 : ! Inputs:
392 : ! scmaker : supercell_maker_t
393 : ! Output:
394 : ! scpot: a class pointer to an ABSTRACT potential.
395 : ! Not a potential because this is inherited from
396 : ! an abstract_primitive_potential_t, which doesn't know
397 : ! the type of the supercell potential.
398 : !
399 : !-------------------------------------------------------------------!
400 1 : subroutine fill_supercell(self, scmaker, params, scpot, supercell)
401 : use m_spmat_convert, only: COO_to_dense
402 :
403 : class(lwf_primitive_potential_t), intent(inout):: self
404 : type(supercell_maker_t), intent(inout):: scmaker
405 : type(multibinit_dtset_type), intent(inout):: params
406 : class(abstract_potential_t), pointer, intent(inout):: scpot
407 : type(mbsupercell_t), target:: supercell
408 :
409 : integer:: nlwf, sc_nlwf
410 : integer:: inz, iR, R(3), i, j, icell
411 1 : integer, allocatable:: ilist_sc(:), jlist_sc(:), Rlist_sc(:,:)
412 : real(dp):: val
413 :
414 1 : ABI_UNUSED_A(params)
415 :
416 :
417 1 : nlwf = self%nlwf
418 1 : sc_nlwf = nlwf*scmaker%ncells
419 :
420 : ! Step 1: allocate the scpot as a corresponding supercell potential
421 1 : ABI_MALLOC_TYPE_SCALAR(lwf_potential_t, scpot)
422 : ! Fortran does not know the functions specific to the derived class pointer.
423 : ! Only the ones inheritated from abstract class,
424 : ! unless select type is used:
425 : select type(scpot)
426 : type is (lwf_potential_t)
427 1 : call scpot%initialize(sc_nlwf)
428 1 : call scpot%set_supercell(supercell)
429 : ! IFC is an COO_mat_t, which has the index of R1, R2, R3, i, j and the value of val
430 : ! list of index R: coeff%ind%data(1, 1:coeff%nnz)
431 : ! list of i: coeff%ind%data(2, 1:coeff%nnz)
432 : ! list of j: coeff%ind%data(3, 1:coeff%nnz)
433 : ! IFC: (ind) = val
434 1 : if(self%as_lattice_anharmonic) then
435 0 : call scpot%use_as_lattice_anharmonic()
436 : else
437 : ! harmonic terms
438 133 : do inz = 1, self%coeff%nnz
439 : ! For each non-zero entry in the coeff matrix
440 : ! get the R, i, and j, val
441 132 : iR = self%coeff%ind%data(1, inz)
442 528 : R(:) = self%Rlist(:, iR)
443 132 : i = self%coeff%ind%data(2, inz)
444 132 : j = self%coeff%ind%data(3, inz)
445 132 : val = self%coeff%val%data(inz)
446 : ! translate i to i in supercell.
447 : ! No need to allocate, it is done by trans_i . but remember to deallocate!
448 : ! nbasis is the number in one primitive cell.
449 : ! e.g. there are 3*natom possible i (3: x, y, z) in each primitive cell.
450 132 : call scmaker%trans_i(nbasis = self%nlwf, i = i, i_sc = ilist_sc )
451 : ! translate j, Rj to supercell.
452 132 : call scmaker%trans_j_and_Rj(nbasis = self%nlwf, j = j, Rj = R, j_sc = jlist_sc, Rj_sc = Rlist_sc)
453 : ! values are repeated in cells
454 67716 : do icell = 1, scmaker%ncells
455 67716 : call scpot%add_term(ilist_sc(icell), jlist_sc(icell), val )
456 : end do
457 132 : ABI_SFREE(ilist_sc)
458 132 : ABI_SFREE(jlist_sc)
459 133 : ABI_SFREE(Rlist_sc)
460 : end do
461 : end if
462 :
463 : ! coefficients of atomic displacements in supercell
464 1 : call scpot%supercell%lwf%lwf_latt_coeffs%initialize(self%nlwf*scmaker%ncells, self%natom*3*scmaker%ncells )
465 3 : do i = 1, self%nlwf
466 2 : call scmaker%trans_i(nbasis = nlwf, i = i, i_sc = ilist_sc)
467 :
468 92 : do iR = 1, self%nR
469 360 : R(:) = self%Rlist(:, iR)
470 1712 : do j = 1, self%natom*3
471 1620 : val = self%lattice_coeffs(i, j, iR) ! ilwf, iatom3, iR
472 1710 : if (abs(val) > 1e-4) then
473 1348 : call scmaker%trans_j_and_Rj(nbasis = self%natom*3, j = j, Rj = R, j_sc = jlist_sc, Rj_sc = Rlist_sc)
474 691524 : do icell = 1, scmaker%ncells
475 : !call scpot%supercell%lwf%lwf_latt_coeffs(ilist_sc(icell))%push(jlist_sc(icell), val)
476 2071876 : call scpot%supercell%lwf%lwf_latt_coeffs%coeffs%add_entry([jlist_sc(icell), ilist_sc(icell)], val)
477 : end do
478 1348 : ABI_SFREE(jlist_sc)
479 1348 : ABI_SFREE(Rlist_sc)
480 : end if
481 : end do
482 : end do
483 3 : ABI_SFREE(ilist_sc)
484 : end do
485 :
486 : ! anharmonic terms
487 : !call scpot%set_ref_energy(self%ref_energy*scmaker%ncells)
488 : !if (self%has_self_bound_term) then
489 : ! call scpot%add_self_bound_term(self%self_bound_order, self%self_bound_coeff)
490 : !end if
491 :
492 2 : if (self%onebody_nterm /= 0) then
493 7 : do i = 1, self%onebody_nterm
494 6 : call scmaker%trans_i(nbasis = self%nlwf, i = self%onebody_i(i), i_sc = ilist_sc)
495 3078 : do icell = 1, scmaker%ncells
496 3078 : call scpot%add_onebody_term(i = ilist_sc(icell), order = self%onebody_order(i), val = self%onebody_val(i))
497 : end do
498 7 : ABI_SFREE(ilist_sc)
499 : end do
500 : endif
501 : ! Test the phonon energy
502 : !call COO_to_dense(scpot%coeff, real_sc_evecs)
503 : !sc_evecs(:,:) = real_sc_evecs(:,:)
504 : !call eigensh(sc_evals, sc_evecs)
505 :
506 : end select
507 :
508 :
509 1 : end subroutine fill_supercell
510 :
511 : !-------------------------------------------------------------------!
512 : ! calculate hamiltonian at k point (or more precisely, q point)
513 : ! H_ij(k) = sum_R H_ij(R) exp( i2pi k.dot.R)
514 : !-------------------------------------------------------------------!
515 0 : subroutine get_hamk(self, kpoint, hamk)
516 : class(lwf_primitive_potential_t), intent(in):: self
517 : real(dp), intent(in):: kpoint(3)
518 : complex(dp), intent(inout):: hamk(:,:)
519 : integer:: inz, iR, R(3), i, j
520 : real(dp):: val
521 0 : hamk(:,:) = cmplx(0.0, 0.0)
522 0 : do inz = 1, self%coeff%nnz
523 0 : iR = self%coeff%ind%data(1, inz)
524 0 : R(:) = self%Rlist(:, iR)
525 0 : i = self%coeff%ind%data(2, inz)
526 0 : j = self%coeff%ind%data(3, inz)
527 0 : val = self%coeff%val%data(inz)
528 0 : hamk(i, j) =hamk(i, j) + val*exp(cmplx(0.0, 2.0) *pi*dot_product(kpoint, R))
529 : end do
530 0 : end subroutine get_hamk
531 :
532 : !-------------------------------------------------------------------!
533 : !calculate eigenvalue and eigen vector
534 : !-------------------------------------------------------------------!
535 0 : subroutine get_eigen(self, kpoint, evals, evecs)
536 : class(lwf_primitive_potential_t), intent(in):: self
537 : real(dp), intent(in):: kpoint(3)
538 : real(dp), intent(inout):: evals(:)
539 : complex(dp), intent(inout):: evecs(:,:)
540 0 : call self%get_hamk(kpoint, evecs)
541 : ! The evecs array is reused both as the matrix and eigenvectors.
542 0 : call eigensh(evals, evecs)
543 0 : end subroutine get_eigen
544 :
545 : !subroutine add_self_bound_term(self, order, coeff)
546 : ! class(lwf_primitive_potential_t), intent(inout):: self
547 : ! integer, intent(in):: order
548 : ! real(dp), intent(in):: coeff
549 : ! if (order /= 0) then
550 : ! self%has_self_bound_term=.True.
551 : ! self%self_bound_order = order
552 : ! self%self_bound_coeff = coeff
553 : ! end if
554 :
555 : !end subroutine add_self_bound_term
556 :
557 :
558 4 : end module m_lwf_primitive_potential
|
