Line data Source code
1 : !!****m* ABINIT/m_lwf_potential
2 : !! NAME
3 : !! m_lwf_potential
4 : !!
5 : !! FUNCTION
6 : !! This module contains an LWF potential.
7 : !!
8 : !! Datatypes:
9 : !!
10 : !! Subroutines:
11 : !! TODO: add this when F2003 doc style is determined.
12 : !!
13 : !!
14 : !! COPYRIGHT
15 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
16 : !! This file is distributed under the terms of the
17 : !! GNU General Public License, see ~abinit/COPYING
18 : !! or http://www.gnu.org/copyleft/gpl.txt .
19 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
20 : !!
21 : !! SOURCE
22 :
23 :
24 : #if defined HAVE_CONFIG_H
25 : #include "config.h"
26 : #endif
27 :
28 : #include "abi_common.h"
29 :
30 :
31 : module m_lwf_potential
32 : use defs_basis
33 : use m_errors
34 : use m_abicore
35 : use m_abstract_potential, only: abstract_potential_t
36 : use m_spmat_ndcoo, only: ndcoo_mat_t
37 : use m_spmat_coo, only: COO_mat_t
38 : use m_spmat_csr, only: CSR_mat_t
39 : use m_spmat_spvec, only: sp_real_vec
40 : use m_spmat_convert, only: spmat_convert
41 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
42 : use m_hashtable_strval, only: hash_table_t
43 : use m_twobody_interaction, only: get_twobody_dEdx, get_twobody_delta_E
44 : use m_multibinit_dataset, only: multibinit_dtset_type
45 : use m_lattice_lwf_map, only: lwf_latt_coeff_t
46 : implicit none
47 : !!***
48 :
49 : private
50 :
51 : type, public, extends(abstract_potential_t) :: lwf_potential_t
52 : integer :: nlwf ! number of lwf
53 : real(dp) :: ref_energy=0.0 ! reference energy
54 : logical :: csr_mat_ready=.False.
55 : type(COO_mat_t) :: coeff_coo ! coefficient. A COO sparse matrix (3N*3N).
56 : type(CSR_mat_t) :: coeff
57 : type(NDCOO_mat_t) :: onebody_coeff ! onebody anharmonic
58 : type(NDCOO_mat_t) :: coeff2 ! twobody anharmonic
59 :
60 : logical :: use_harmonic = .True.
61 : !logical :: has_self_bound_term = .False.
62 : !integer :: self_bound_order=0
63 : !real(dp) :: self_bound_coeff=0.0_dp
64 :
65 : real(dp) :: beta
66 : real(dp), allocatable :: coeff_diag(:)
67 :
68 : logical :: as_lattice_anharmonic=.False.
69 : ! tmp arrays for forces and lwf
70 : real(dp), allocatable :: lwf_force(:), lwf_amp(:)
71 : contains
72 : procedure :: initialize
73 : procedure :: finalize
74 : procedure :: use_as_lattice_anharmonic
75 : procedure :: set_supercell
76 : procedure :: set_ref_energy
77 : procedure :: set_params
78 : procedure :: calculate
79 : procedure :: add_term
80 : procedure :: convert_coeff_to_csr
81 : procedure :: get_delta_E_lwf
82 : procedure :: add_onebody_term
83 : !procedure :: add_self_bound_term
84 : end type lwf_potential_t
85 :
86 : contains
87 :
88 : !-------------------------------------------------------------------!
89 : ! initialize
90 : ! Input:
91 : ! natom: number of atoms
92 : !-------------------------------------------------------------------!
93 1 : subroutine initialize(self, nlwf)
94 : class(lwf_potential_t), intent(inout) :: self
95 : integer, intent(in) :: nlwf
96 1 : self%has_lwf= .True.
97 1 : self%is_null = .False.
98 1 : self%label="lwf_potential"
99 1 : self%nlwf=nlwf
100 1 : self%ref_energy=0.0_dp
101 3 : call self%coeff_coo%initialize(mshape= [self%nlwf, self%nlwf])
102 3 : call self%onebody_coeff%initialize(mshape= [self%nlwf, -1])
103 3 : call self%coeff2%initialize(mshape= [self%nlwf, self%nlwf])
104 1 : self%csr_mat_ready=.False.
105 3 : ABI_MALLOC(self%coeff_diag, (self%nlwf))
106 :
107 2 : ABI_MALLOC(self%lwf_force, (self%nlwf))
108 2 : ABI_MALLOC(self%lwf_amp, (self%nlwf))
109 :
110 1 : end subroutine initialize
111 :
112 : !-------------------------------------------------------------------!
113 : ! Finalize
114 : !-------------------------------------------------------------------!
115 1 : subroutine finalize(self)
116 : class(lwf_potential_t), intent(inout) :: self
117 1 : self%has_displacement=.False.
118 1 : call self%onebody_coeff%finalize()
119 1 : call self%coeff2%finalize()
120 1 : call self%abstract_potential_t%finalize()
121 :
122 1 : if (.not. self%csr_mat_ready) then
123 0 : call self%coeff_coo%finalize()
124 : else
125 1 : call self%coeff%finalize()
126 : end if
127 1 : ABI_SFREE(self%coeff_diag)
128 1 : self%nlwf=0
129 1 : self%is_null=.True.
130 1 : self%as_lattice_anharmonic=.False.
131 1 : ABI_SFREE(self%lwf_force)
132 1 : ABI_SFREE(self%lwf_amp)
133 1 : self%csr_mat_ready=.False.
134 1 : end subroutine finalize
135 :
136 0 : subroutine use_as_lattice_anharmonic(self)
137 : class(lwf_potential_t), intent(inout) :: self
138 0 : self%as_lattice_anharmonic=.True.
139 0 : self%use_harmonic = .False.
140 0 : self%has_displacement = .True.
141 0 : self%has_lwf= .False.
142 0 : self%is_null = .False.
143 0 : self%label="latt_lwf_potential"
144 0 : end subroutine use_as_lattice_anharmonic
145 :
146 :
147 :
148 : !===============================================================
149 : !
150 : !> @
151 : !===============================================================
152 : subroutine add_lattice_coeffs(self, ilwf, ilatt, val)
153 : class(lwf_potential_t), intent(inout) :: self
154 : integer , intent(in) :: ilwf, ilatt
155 : real(dp) , intent(in) :: val
156 : !call self%lwf_latt_coeffs(ilwf)%push(ilatt, val)
157 : call self%supercell%lwf%lwf_latt_coeffs%coeffs%add_entry( [ilatt, ilwf], val)
158 : end subroutine add_lattice_coeffs
159 :
160 : !-------------------------------------------------------------------!
161 : ! Add a term to the potential
162 : !-------------------------------------------------------------------!
163 67584 : subroutine add_term(self, i,j, val)
164 : class(lwf_potential_t), intent(inout) :: self
165 : integer, intent(in) :: i, j
166 : real(dp), intent(in) :: val
167 202752 : call self%coeff_coo%add_entry([i,j], val)
168 67584 : end subroutine add_term
169 :
170 3072 : subroutine add_onebody_term(self, i, order, val)
171 : class(lwf_potential_t), intent(inout) :: self
172 : integer, intent(in) :: i, order
173 : real(dp), intent(in) :: val
174 9216 : call self%onebody_coeff%add_entry([i, order], val)
175 3072 : end subroutine add_onebody_term
176 :
177 : subroutine add_higher_order_term(self, i, j, orderi, orderj, val)
178 : class(lwf_potential_t), intent(inout) :: self
179 : integer, intent(in) :: i, j, orderi, orderj
180 : real(dp), intent(in) :: val
181 : call self%coeff2%add_entry([i,j, orderi, orderj], val)
182 : end subroutine add_higher_order_term
183 :
184 42000 : subroutine convert_coeff_to_csr(self)
185 : class(lwf_potential_t), intent(inout) :: self
186 42000 : if (.not. self%csr_mat_ready) then
187 : !call init_mpi_info(master, iam_master, my_rank, comm, nproc)
188 : !if(iam_master) then
189 1 : call spmat_convert(self%coeff_coo, self%coeff)
190 1 : call self%coeff_coo%diag(self%coeff_diag)
191 1 : call self%coeff_coo%finalize()
192 : !endif
193 : !call self%bilinear_csr_mat%sync(master=master, comm=comm, nblock=1)
194 1 : self%csr_mat_ready=.True.
195 : !call xmpi_bcast(self%csr_mat_ready, master, comm, ierr)
196 : endif
197 :
198 42000 : end subroutine convert_coeff_to_csr
199 :
200 :
201 : !-------------------------------------------------------------------!
202 : ! Set the reference energy.
203 : !-------------------------------------------------------------------!
204 0 : subroutine set_ref_energy(self, ref_energy)
205 : class(lwf_potential_t), intent(inout) :: self
206 : real(dp), intent(in) :: ref_energy
207 0 : self%ref_energy=ref_energy
208 0 : end subroutine set_ref_energy
209 :
210 : !-------------------------------------------------------------------!
211 : ! set_supercell
212 : ! link the supercell with potential.
213 : ! Inputs:
214 : ! supercell: mbsupercell_t
215 : !-------------------------------------------------------------------!
216 2 : subroutine set_supercell(self, supercell)
217 : class(lwf_potential_t), intent(inout) :: self
218 : type(mbsupercell_t), target, intent(inout) :: supercell
219 2 : self%supercell => supercell
220 2 : end subroutine set_supercell
221 :
222 :
223 : !-------------------------------------------------------------------!
224 : ! calculate force and energy from harmonic potential
225 : ! F= - IFC .matmul. displacement
226 : ! E = 1/2 (-F) .dot. displacement = 1/2<disp|IFC|disp>
227 : ! Input:
228 : ! displacement: required.
229 : !-------------------------------------------------------------------!
230 42000 : subroutine calculate(self, displacement, strain, spin, lwf, force, stress, bfield, lwf_force, energy, energy_table)
231 : ! This function calculate the energy and its first derivative
232 : ! the inputs and outputs are optional so that each effpot can adapt to its
233 : ! own.
234 : ! In principle, the 1st derivatives are only calculated if asked to (present). However, they can be computed if it is simply convinient to do.
235 : class(lwf_potential_t), intent(inout) :: self ! the effpot may save the states.
236 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
237 : real(dp), optional, intent(inout) :: force(:,:), stress(:,:), bfield(:,:), lwf_force(:), energy
238 : type(hash_table_t),optional, intent(inout) :: energy_table
239 : real(dp) :: etmp, val
240 : integer :: i,ilwf, order
241 : ! if present in input
242 : ! calculate if required
243 42000 : ABI_UNUSED_A(strain)
244 42000 : ABI_UNUSED_A(spin)
245 42000 : ABI_UNUSED_A(stress)
246 42000 : ABI_UNUSED_A(bfield)
247 :
248 : !if(.not. present(lwf)) then
249 : ! MSG_BUG("lwf not exist")
250 : !end if
251 :
252 :
253 42000 : if(self%as_lattice_anharmonic) then
254 0 : call self%supercell%lwf%lwf_latt_coeffs%lattice_to_lwf_projection( displacement, lwf)
255 : end if
256 :
257 :
258 43050000 : self%lwf_force(:) =0.0_dp
259 42000 : etmp=0.0_dp
260 : ! Harmonic term
261 42000 : if (self%use_harmonic) then
262 42000 : call self%convert_coeff_to_csr()
263 42000 : call self%coeff%mv(lwf, self%lwf_force)
264 43050000 : etmp=etmp+0.5_dp * sum(self%lwf_force*lwf)
265 43050000 : self%lwf_force(:) = -self%lwf_force(:)
266 : end if
267 :
268 : ! self_bound_term as from the input parameters
269 : !if (self%has_self_bound_term) then
270 : ! self%lwf_force(:) = self%lwf_force(:) - &
271 : ! & self%self_bound_coeff*self%self_bound_order* lwf**(self%self_bound_order-1)
272 : ! etmp = etmp + &
273 : ! & self%self_bound_coeff*sum(lwf**(self%self_bound_order))
274 : !endif
275 :
276 42000 : if (self%onebody_coeff%nnz/= 0) then
277 129066000 : do i =1, self%onebody_coeff%nnz
278 129024000 : ilwf= self%onebody_coeff%ind%data(1, i)
279 129024000 : order= self%onebody_coeff%ind%data(2, i)
280 129024000 : val=self%onebody_coeff%val%data(i)
281 129024000 : etmp= etmp + val*lwf(ilwf)**order
282 129066000 : self%lwf_force(ilwf) =self%lwf_force(ilwf) - val*order*lwf(ilwf)**(order-1)
283 : end do
284 : end if
285 :
286 : !TODO: remove. For testing only
287 : !etmp = etmp+ self%beta*sum(lwf(::2)**2 * lwf(1::2)**2)
288 :
289 42000 : if(self%as_lattice_anharmonic) then
290 0 : call self%supercell%lwf%lwf_latt_coeffs%lwf_force_to_lattice(self%lwf_force, force)
291 : else
292 43050000 : lwf_force=lwf_force+self%lwf_force
293 : end if
294 :
295 42000 : if (present(energy)) then
296 42000 : energy=energy+etmp
297 : endif
298 42000 : if(present(energy_table)) then
299 42000 : call energy_table%put(self%label, etmp)
300 : endif
301 42000 : end subroutine calculate
302 :
303 : !----------------------------------------------------------------------
304 : !> @brief get_delta_E_lwf: calculate the energy difference when a given lwf
305 : !> is changed. This is to be used for spin Monte Carlo. Currently the
306 : !> only supported is the spin model.
307 : !>
308 : !> @param[in] lwf: lwf of full structure. array of (nlwf)
309 : !> @param[in] ilwf: the index of spin changed. integer
310 : !> @param[in] lwf_new: the new value of the changed spin.
311 : !> @param[out] deltaE: the energy difference
312 : !----------------------------------------------------------------------
313 0 : subroutine get_delta_E_lwf(self, lwf, ilwf, lwf_new, deltaE)
314 : class(lwf_potential_t), intent(inout) :: self ! the effpot may save the states.
315 : real(dp), intent(inout) :: lwf(:), lwf_new
316 : integer, intent(in) :: ilwf
317 : real(dp), intent(inout) :: deltaE
318 : real(dp) :: tmp, dlwf, lold
319 :
320 0 : lold=lwf(ilwf)
321 0 : dlwf=lwf_new-lold
322 0 : lwf(ilwf) = lwf_new
323 :
324 0 : call self%convert_coeff_to_csr()
325 :
326 0 : deltaE=0.0_dp
327 0 : tmp=0.0_dp
328 0 : call self%coeff%mv_one_row(ilwf, lwf, tmp)
329 0 : deltaE=deltaE+ tmp*dlwf - 0.5* self%coeff_diag(ilwf)*dlwf*dlwf
330 :
331 : ! bound term
332 : !if (self%has_self_bound_term) then
333 : ! deltaE= deltaE+ &
334 : ! & self%self_bound_coeff*(lwf_new**(self%self_bound_order) &
335 : ! & - lold**(self%self_bound_order))
336 : !! Adding x^6 and x^8 for VO2
337 : ! deltaE=deltaE -1.1344*(lwf_new**6-lold**6) + 0.438*(lwf_new**8-lold**8)
338 : !end if
339 :
340 : ! (Q1 Q2)^2 term
341 : ! if(modulo(ilwf, 2)==0) then
342 : ! deltaE=deltaE+ self%beta*((lwf_new**2- lold**2)*lwf(ilwf-1)**2)
343 : ! else
344 : ! deltaE=deltaE+ self%beta*((lwf_new**2- lold**2)*lwf(ilwf+1)**2)
345 : ! endif
346 :
347 : ! if (self%onebody_coeff%nnz/= 0) then
348 : ! do i =1, self%onebody_coeff%nnz
349 : ! ilwf= self%onebody_coeff%ind%data(i, 1)
350 : ! order= self%onebody_coeff%ind%data(i, 2)
351 : ! val=self%onebody_coeff%val%data(i)
352 : ! deltaE= deltaE + val*self%lwf_amp(ilwf)**order
353 : ! self%lwf_force(ilwf) =self%lwf_force(ilwf) - val*order*self%lwf_amp(ilwf)
354 : ! end do
355 : ! end if
356 :
357 0 : if(modulo(ilwf, 2)==0) then
358 0 : deltaE=deltaE+ self%beta*((lwf_new**2- lold**2)*lwf(ilwf-1)**2)
359 : else
360 0 : deltaE=deltaE+ self%beta*((lwf_new**2- lold**2)*lwf(ilwf+1)**2)
361 : endif
362 :
363 0 : lwf(ilwf)=lwf(ilwf)-dlwf
364 0 : end subroutine get_delta_E_lwf
365 :
366 : !subroutine add_self_bound_term(self, order, coeff)
367 : ! class(lwf_potential_t), intent(inout) :: self
368 : ! integer, intent(in) :: order
369 : ! real(dp), intent(in) :: coeff
370 : ! if (order /= 0) then
371 : ! self%has_self_bound_term=.True.
372 : ! self%self_bound_order=order
373 : ! self%self_bound_coeff=coeff
374 : ! end if
375 : !end subroutine add_self_bound_term
376 :
377 : !----------------------------------------------------------------------
378 : !> @brief set_params: set the parameters from input file parameters
379 : !>
380 : !> @param[in] params: multibinit_dtset_type: from input file
381 : !----------------------------------------------------------------------
382 1 : subroutine set_params(self, params)
383 : class(lwf_potential_t), intent(inout) :: self
384 : type(multibinit_dtset_type), intent(inout) :: params
385 1 : self%beta=params%spin_damping
386 1 : end subroutine set_params
387 :
388 :
389 3 : end module m_lwf_potential
|
