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