Branch data 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
+ + - + -
+ - + - +
- + + - -
+ ]
|
