Line data Source code
1 : !!****m* ABINIT/m_lattice_harmonic_primitive_potential
2 : !! NAME
3 : !! m_lattice_harmonic_primitive_potential
4 : !!
5 : !! FUNCTION
6 : !! This module contains the example lattice 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 : !! lattice_harmonic_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 : module m_lattice_harmonic_primitive_potential
31 : use, intrinsic :: iso_c_binding
32 : !use m_dynamic_array, only: int_array_type, real_array_type, int2d_array_type
33 : use defs_basis
34 : use m_abicore
35 : use m_errors
36 : use m_nctk
37 : #if defined HAVE_NETCDF
38 : use netcdf
39 : #endif
40 : use m_xmpi
41 : use m_mathfuncs, only: eigensh
42 : use m_multibinit_dataset, only: multibinit_dtset_type
43 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
44 : use m_primitive_potential, only: primitive_potential_t
45 : use m_abstract_potential, only: abstract_potential_t
46 : use m_dynamic_array, only: int2d_array_type
47 : use m_supercell_maker, only: supercell_maker_t
48 : use m_spmat_ndcoo, only: ndcoo_mat_t
49 : use m_lattice_harmonic_potential, only: lattice_harmonic_potential_t
50 : implicit none
51 : private
52 : !!***
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) :: lattice_harmonic_primitive_potential_t
61 : integer :: natom ! number of atoms
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 : integer, allocatable :: Rlist(:,:) ! The list of R points (3, number of R-points)
65 : !. Rlist(:, ind_R) is a R-vector.
66 : real(dp) :: ref_energy=0.0 ! reference energy
67 : contains
68 : procedure:: initialize
69 : procedure:: finalize
70 : procedure :: load_from_files ! load potential from files listed in files file
71 : procedure :: load_from_netcdf ! load potential from a netcdf file
72 : procedure:: fill_supercell ! fill a supercell potential
73 : procedure :: get_hamk ! generate hamiltonian for one k point.
74 : procedure :: get_eigen ! eigen values and eigen vectors
75 : end type lattice_harmonic_primitive_potential_t
76 :
77 : contains
78 :
79 : !-------------------------------------------------------------------!
80 : ! Initialize
81 : ! Input:
82 : ! primcell: the reference primitive cell.
83 : !-------------------------------------------------------------------!
84 5 : subroutine initialize(self, primcell)
85 : class(lattice_harmonic_primitive_potential_t), intent(inout) :: self
86 : type(mbcell_t), target, intent(inout) :: primcell
87 : !integer, intent(in) :: nspin
88 5 : self%primcell=>primcell
89 5 : self%label="lattice_harmonic_primitive_potential"
90 5 : self%has_spin=.False.
91 5 : self%has_displacement=.True.
92 5 : self%has_strain=.False.
93 5 : self%has_lwf=.False.
94 5 : end subroutine initialize
95 :
96 :
97 : !-------------------------------------------------------------------!
98 : ! Finalize
99 : !-------------------------------------------------------------------!
100 5 : subroutine finalize(self)
101 : class(lattice_harmonic_primitive_potential_t), intent(inout) :: self
102 5 : call self%coeff%finalize()
103 5 : ABI_FREE(self%Rlist)
104 5 : nullify(self%primcell)
105 5 : self%natom=0
106 5 : self%label="Destroyed lattice_harmonic_primitive_potential"
107 5 : call self%primitive_potential_t%finalize()
108 5 : end subroutine finalize
109 :
110 :
111 : !-------------------------------------------------------------------!
112 : ! Load from files:
113 : !> params: parameters
114 : !> fnames: file names from files file
115 : !-------------------------------------------------------------------!
116 5 : subroutine load_from_files(self, params, fnames)
117 : class(lattice_harmonic_primitive_potential_t), intent(inout) :: self
118 : type(multibinit_dtset_type), intent(in) :: params
119 : character(len=fnlen), intent(in) :: fnames(:)
120 : ABI_UNUSED(fnames)
121 5 : if(trim(params%latt_pot_fname)=='') then
122 0 : call self%load_from_netcdf(fnames(3))
123 : else
124 5 : call self%load_from_netcdf(params%latt_pot_fname)
125 : end if
126 5 : end subroutine load_from_files
127 :
128 : !-------------------------------------------------------------------!
129 : ! load potential from netcdf file
130 : ! Note that the lattic part of the primitive cell is also loaded.
131 : ! Input:
132 : ! fname: filename
133 : !-------------------------------------------------------------------!
134 5 : subroutine load_from_netcdf(self, fname)
135 : class(lattice_harmonic_primitive_potential_t), intent(inout) :: self
136 : character(len=fnlen), intent(in) :: fname
137 : integer :: ncid, ierr
138 : integer :: iR, nR, natom, natom3
139 : real(dp) :: ref_energy
140 : real(dp) :: cell(3,3)
141 5 : real(dp), allocatable :: xcart(:,:), masses(:)
142 5 : integer, allocatable :: zion(:)
143 5 : real(dp), allocatable :: ifc_vallist(:,:,:)
144 : integer :: varid, i, j
145 : #if defined HAVE_NETCDF
146 5 : ierr=nf90_open(trim(fname), NF90_NOWRITE, ncid)
147 5 : NCF_CHECK_MSG(ierr, "Open netcdf file "//trim(fname))
148 :
149 5 : ierr=nctk_get_dim(ncid, "ifc_nR" , nR)
150 5 : ierr=nctk_get_dim(ncid, "natom", natom)
151 5 : ierr=nctk_get_dim(ncid, "natom3", natom3)
152 :
153 :
154 15 : ABI_MALLOC(masses, (natom))
155 15 : ABI_MALLOC(xcart, (3, natom))
156 15 : ABI_MALLOC(zion,(natom))
157 :
158 25 : ABI_MALLOC(ifc_vallist, (natom3, natom3, nR))
159 15 : ABI_MALLOC(self%Rlist,(3, nR))
160 :
161 :
162 5 : ierr =nf90_inq_varid(ncid, "ref_energy", varid)
163 5 : NCF_CHECK_MSG(ierr, "ref_energy")
164 5 : ierr = nf90_get_var(ncid, varid, ref_energy)
165 5 : NCF_CHECK_MSG(ierr, "ref_energy")
166 :
167 :
168 5 : ierr =nf90_inq_varid(ncid, "ref_masses", varid)
169 5 : NCF_CHECK_MSG(ierr, "ref_masses")
170 5 : ierr = nf90_get_var(ncid, varid, masses)
171 5 : NCF_CHECK_MSG(ierr, "ref_masses")
172 :
173 5 : ierr =nf90_inq_varid(ncid, "ref_xcart", varid)
174 5 : NCF_CHECK_MSG(ierr, "ref_xcart")
175 5 : ierr = nf90_get_var(ncid, varid, xcart)
176 5 : NCF_CHECK_MSG(ierr, "ref_xcart")
177 :
178 5 : ierr =nf90_inq_varid(ncid, "ref_cell", varid)
179 5 : NCF_CHECK_MSG(ierr, "rec_cell")
180 5 : ierr = nf90_get_var(ncid, varid, cell)
181 5 : NCF_CHECK_MSG(ierr, "ref_cell")
182 :
183 5 : ierr =nf90_inq_varid(ncid, "ref_zion", varid)
184 5 : NCF_CHECK_MSG(ierr, "ref_zion")
185 5 : ierr = nf90_get_var(ncid, varid, zion)
186 5 : NCF_CHECK_MSG(ierr, "ref_zion")
187 :
188 : ! Unit conversions
189 5 : ref_energy = ref_energy * eV_Ha
190 30 : masses(:) = masses(:) * amu_emass
191 65 : cell(:,:) = cell(:,:) / Bohr_Ang
192 105 : xcart(:,:) = xcart(:,:) /Bohr_Ang
193 :
194 5 : self%natom=natom
195 5 : self%ref_energy=ref_energy
196 5 : call self%primcell%set_lattice(natom, cell, xcart, masses, zion)
197 :
198 20 : call self%coeff%initialize([ nR, natom3, natom3 ])
199 5 : ierr =nf90_inq_varid(ncid, "ifc_Rlist", varid)
200 5 : NCF_CHECK_MSG(ierr, "ifc_Rlist")
201 5 : ierr = nf90_get_var(ncid, varid, self%Rlist)
202 5 : NCF_CHECK_MSG(ierr, "ifc_Rlist")
203 :
204 5 : ierr =nf90_inq_varid(ncid, "ifc_vallist", varid)
205 5 : NCF_CHECK_MSG(ierr, "ifc_vallist")
206 5 : ierr = nf90_get_var(ncid, varid, ifc_vallist)
207 5 : NCF_CHECK_MSG(ierr, "ifc_vallist")
208 :
209 32540 : ifc_vallist(:,:,:) = ifc_vallist(:,:,:) * eV_Ha * (Bohr_Ang * Bohr_Ang)
210 :
211 140 : do iR =1, nR
212 2165 : do i=1 , natom3
213 32535 : do j=1, natom3
214 32400 : if (abs(ifc_vallist(j, i, iR))>1e-3) then
215 : ! NOTE: in fortran the order of index in reversed when reading netcdf array.
216 11460 : call self%coeff%add_entry([iR, i, j], ifc_vallist(j, i, iR))
217 : end if
218 : end do
219 : end do
220 : end do
221 :
222 5 : ierr=nf90_close(ncid)
223 5 : NCF_CHECK_MSG(ierr, "Close netcdf file")
224 :
225 5 : ABI_FREE(masses)
226 5 : ABI_FREE(xcart)
227 5 : ABI_FREE(zion)
228 5 : ABI_FREE(ifc_vallist)
229 : #else
230 : NETCDF_NOTENABLED_ERROR()
231 : #endif
232 :
233 5 : end subroutine load_from_netcdf
234 :
235 :
236 : !-------------------------------------------------------------------!
237 : !Fill supercell
238 : ! Inputs:
239 : ! scmaker : supercell_maker_t
240 : ! Output:
241 : ! scpot: a class pointer to an ABSTRACT potential.
242 : ! Not a harmonic potential because this is inherited from
243 : ! an abstract_primitive_potential_t, which doesn't know
244 : ! the type of the supercell potential.
245 : !
246 : !-------------------------------------------------------------------!
247 5 : subroutine fill_supercell(self, scmaker, params, scpot, supercell)
248 : use m_spmat_convert, only: COO_to_dense
249 :
250 : class(lattice_harmonic_primitive_potential_t) , intent(inout) :: self
251 : type(supercell_maker_t), intent(inout) :: scmaker
252 : type(multibinit_dtset_type), intent(inout) :: params
253 : class(abstract_potential_t), pointer, intent(inout) :: scpot
254 : type(mbsupercell_t), target :: supercell
255 :
256 : integer :: natom, sc_natom
257 : integer :: inz, iR, R(3), i, j, icell
258 5 : integer, allocatable :: ilist_sc(:), jlist_sc(:), Rlist_sc(:,:)
259 : real(dp):: val
260 :
261 5 : ABI_UNUSED_A(params)
262 :
263 :
264 5 : natom=self%natom
265 5 : sc_natom= natom* scmaker%ncells
266 :
267 : !! NOTE: the code below can be used as a pattern to build supercell from primitivecell.
268 : ! Step 1: allocate the scpot as a corresponding supercell potential
269 5 : ABI_MALLOC_TYPE_SCALAR(lattice_harmonic_potential_t, scpot)
270 : ! Fortran does not know the functions specific to the derived class pointer.
271 : ! Only the ones inheritated from abstract class,
272 : ! unless select type is used:
273 : select type(scpot)
274 : type is (lattice_harmonic_potential_t)
275 5 : call scpot%initialize(sc_natom)
276 5 : call scpot%set_supercell(supercell)
277 : ! IFC is an COO_mat_t, which has the index of R1, R2, R3, i, j and the value of val
278 : ! list of index R: coeff%ind%data(1, 1:coeff%nnz)
279 : ! list of i: coeff%ind%data(2, 1:coeff%nnz)
280 : ! list of j: coeff%ind%data(3, 1:coeff%nnz)
281 : ! IFC: (ind) = val
282 2870 : do inz =1 , self%coeff%nnz
283 : ! For each non-zero entry in the coeff matrix
284 : ! get the R, i, and j, val
285 2865 : iR=self%coeff%ind%data(1,inz)
286 11460 : R(:) = self%Rlist(:, iR)
287 2865 : i=self%coeff%ind%data(2, inz)
288 2865 : j=self%coeff%ind%data(3, inz)
289 2865 : val=self%coeff%val%data(inz)
290 : ! translate i to i in supercell.
291 : ! No need to allocate, it is done by trans_i . but remember to deallocate!
292 : ! nbasis is the number in one primitive cell.
293 : ! e.g. there are 3*natom possible i (3: x, y, z) in each primitive cell.
294 2865 : call scmaker%trans_i(nbasis=self%natom*3, i=i, i_sc=ilist_sc )
295 : ! translate j, Rj to supercell.
296 2865 : call scmaker%trans_j_and_Rj(nbasis=self%natom*3, j=j, Rj=R, j_sc=jlist_sc, Rj_sc=Rlist_sc)
297 : ! values are repeated in cells
298 534614 : do icell=1, scmaker%ncells
299 534609 : call scpot%add_term(ilist_sc(icell), jlist_sc(icell), val )
300 : end do
301 : end do
302 :
303 10 : call scpot%set_ref_energy(self%ref_energy * scmaker%ncells)
304 :
305 : ! Test the phonon energy
306 : !call COO_to_dense(scpot%coeff, real_sc_evecs)
307 : !sc_evecs(:,:) = real_sc_evecs(:,:)
308 : !call eigensh(sc_evals, sc_evecs)
309 :
310 : end select
311 :
312 :
313 5 : ABI_SFREE(ilist_sc)
314 5 : ABI_SFREE(jlist_sc)
315 5 : ABI_SFREE(Rlist_sc)
316 :
317 5 : ABI_UNUSED_A(params)
318 :
319 5 : end subroutine fill_supercell
320 :
321 :
322 : !-------------------------------------------------------------------!
323 : ! calculate hamiltonian at k point (or more precisely, q point)
324 : ! H_ij(k) = sum_R H_ij(R) exp( i2pi k.dot.R)
325 : !-------------------------------------------------------------------!
326 0 : subroutine get_hamk(self, kpoint, hamk)
327 : class(lattice_harmonic_primitive_potential_t), intent(in) :: self
328 : real(dp), intent(in) :: kpoint(3)
329 : complex(dp), intent(inout) :: hamk(:,:)
330 : integer :: inz, iR, R(3), i, j
331 : real(dp) :: val
332 0 : hamk(:,:) = cmplx(0.0,0.0)
333 0 : do inz =1, self%coeff%nnz
334 0 : iR=self%coeff%ind%data(1,inz)
335 0 : R(:) = self%Rlist(:, iR)
336 0 : i=self%coeff%ind%data(2, inz)
337 0 : j=self%coeff%ind%data(3, inz)
338 0 : val=self%coeff%val%data(inz)
339 0 : hamk(i, j) =hamk(i, j) + val*exp(cmplx(0.0,2.0) *pi * dot_product(kpoint, R))
340 : end do
341 0 : end subroutine get_hamk
342 :
343 : !-------------------------------------------------------------------!
344 : !calculate eigenvalue and eigen vector
345 : !-------------------------------------------------------------------!
346 0 : subroutine get_eigen(self, kpoint, evals, evecs)
347 : class(lattice_harmonic_primitive_potential_t), intent(in) :: self
348 : real(dp), intent(in) :: kpoint(3)
349 : real(dp), intent(inout) :: evals(:)
350 : complex(dp), intent(inout) :: evecs(:,:)
351 0 : call self%get_hamk(kpoint, evecs)
352 : ! The evecs array is reused both as the matrix and eigenvectors.
353 0 : call eigensh(evals, evecs)
354 0 : end subroutine get_eigen
355 :
356 20 : end module m_lattice_harmonic_primitive_potential
|
