Line data Source code
1 : !!****m* ABINIT/m_lwf_mc_mover
2 : !! NAME
3 : !! m_lwf_mc_mover
4 : !!
5 : !! FUNCTION
6 : !! This module contains the lwf Markov chain Monte Carlo functions for lwf mover .
7 : !!
8 : !!
9 : !! Datatypes:
10 : !!
11 : !! * lwf_mc_t : MCMC. It defines how to move lwfs in one step,
12 : !! attempt function: whether to accept move
13 : !! accecpt/reject method which define what to do if move is
14 : !! accepted or rejected!! .
15 : !!
16 : !! Subroutines:
17 : !! TODO: add this when F2003 doc style is determined.
18 : !!
19 : !!
20 : !! COPYRIGHT
21 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
22 : !! This file is distributed under the terms of the
23 : !! GNU General Public License, see ~abinit/COPYING
24 : !! or http://www.gnu.org/copyleft/gpl.txt .
25 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
26 : !!
27 : !! SOURCE
28 :
29 : #if defined HAVE_CONFIG_H
30 : #include "config.h"
31 : #endif
32 : #include "abi_common.h"
33 : module m_lwf_mc_mover
34 : use defs_basis
35 : use m_abicore
36 : use m_errors
37 : use m_abstract_potential, only: abstract_potential_t
38 : use m_random_xoroshiro128plus, only: rng_t
39 : use m_hashtable_strval, only: hash_table_t
40 : use m_multibinit_dataset, only: multibinit_dtset_type
41 : use m_lwf_mover, only: lwf_mover_t
42 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
43 : implicit none
44 : !!***
45 : private
46 :
47 : !----------------------------------------------------------------------
48 : !> @brief An helper type to run lwf dynamics
49 : ! The Metropolis-Hasting algorithm is used.
50 : !----------------------------------------------------------------------
51 : type,public, extends(lwf_mover_t) :: lwf_mc_t
52 : real(dp) :: avg_amp! an angle to rotate by average
53 : real(dp) :: deltaE ! energy and the change of energy when change one lwf
54 : real(dp) :: lwf_new, lwf_old! energy and the change of energy when change one lwf
55 : real(dp) :: beta ! 1/(kb T)
56 : integer :: nstep ! number of steps
57 : integer :: imove ! index of lwf to be moved
58 : integer :: naccept ! number of accepted steps
59 : integer :: nattempt
60 : contains
61 : procedure :: initialize
62 : procedure :: finalize
63 : procedure, private :: attempt
64 : procedure, private :: accept
65 : procedure, private :: reject
66 : procedure :: run_one_step
67 : procedure :: set_temperature
68 : procedure, private :: run_one_mc_step
69 : end type lwf_mc_t
70 :
71 :
72 : contains
73 :
74 0 : subroutine initialize(self, params, supercell, rng)
75 : class(lwf_mc_t), intent(inout) :: self
76 : type(multibinit_dtset_type),target, intent(in) :: params
77 : type(mbsupercell_t),target, intent(in) :: supercell
78 : type(rng_t), target, intent(in) :: rng
79 0 : call self%lwf_mover_t%initialize(params, supercell, rng)
80 0 : self%nstep=self%nlwf
81 0 : self%avg_amp=params%lwf_mc_avg_amp
82 0 : self%temperature=params%lwf_temperature
83 0 : self%beta=1.0/self%temperature ! Kb in a.u. is 1.
84 0 : self%lwf_new=0.0_dp
85 0 : self%lwf_old=0.0_dp
86 0 : self%nattempt=0
87 0 : self%naccept=0
88 0 : end subroutine initialize
89 :
90 : !----------------------------------------------------------------------
91 : !> @brief finalize
92 : !----------------------------------------------------------------------
93 0 : subroutine finalize(self)
94 : class(lwf_mc_t), intent(inout) :: self
95 0 : call self%lwf_mover_t%finalize()
96 0 : self%nstep=0
97 0 : self%avg_amp=0.0
98 0 : self%temperature=0.0
99 0 : self%beta=0.0
100 0 : end subroutine finalize
101 :
102 0 : subroutine set_temperature(self, temperature)
103 : class(lwf_mc_t), intent(inout) :: self
104 : real(dp), intent(in) :: temperature
105 0 : call self%lwf_mover_t%set_temperature(temperature)
106 0 : self%temperature=temperature
107 0 : self%beta=1.0/self%temperature ! Kb in a.u. is 1.
108 0 : end subroutine set_temperature
109 :
110 :
111 : !----------------------------------------------------------------------
112 : !> @brief run one monte carlo step
113 : !> @param[in] rngL rundom number generator
114 : !> @param[in] effpot: effective lwf potential
115 : !----------------------------------------------------------------------
116 0 : subroutine run_one_mc_step(self, effpot)
117 : class(lwf_mc_t) :: self
118 : class(abstract_potential_t), intent(inout) :: effpot
119 : real(dp) :: r
120 :
121 : ! try to change lwf
122 0 : r=self%attempt(self%rng, effpot)
123 : ! metropolis-hastings
124 0 : self%nattempt=self%nattempt+1
125 0 : if(self%rng%rand_unif_01()< min(1.0_dp, r) .and. abs(self%lwf_new)<0.5 ) then
126 0 : self%naccept=self%naccept+1
127 0 : call self%accept()
128 : else
129 0 : call self%reject()
130 : end if
131 0 : end subroutine run_one_mc_step
132 :
133 :
134 : !-------------------------------------------------------------------!
135 : !> @brief: Run_one_step
136 : !> effpot: the potential (which do the calculation of E and dE/dvar)
137 : !> param[in]: effpot
138 : !> param[in]: (optional) displacement
139 : !> param[in]: (optional) strain
140 : !> param[in]: (optional) spin
141 : !> param[in]: (optional) lwf
142 : ! NOTE: No need to pass the variable already saved in the mover.
143 : ! e.g. For spin mover, do NOT pass the spin to it.
144 : ! The other variables are only required if there is coupling with
145 : ! the mover variable.
146 : !-------------------------------------------------------------------!
147 0 : subroutine run_one_step(self, effpot, displacement, strain, spin, lwf, energy_table)
148 : ! run one step. (For MC also?)
149 : class(lwf_mc_t), intent(inout) :: self
150 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
151 : type(hash_table_t), optional, intent(inout) :: energy_table
152 : class(abstract_potential_t), intent(inout) :: effpot
153 : integer :: i
154 : character(len=40) :: key
155 :
156 0 : ABI_UNUSED_A(displacement)
157 0 : ABI_UNUSED_A(strain)
158 0 : ABI_UNUSED_A(spin)
159 :
160 : !self%lwf(:)=lwf(:)
161 0 : self%lwf_force(:) = 0.0_dp
162 : ! calculate energy at first step
163 0 : self%energy=0.0
164 :
165 :
166 : call effpot%calculate( displacement=displacement, strain=strain, &
167 : & spin=spin, lwf=self%lwf, lwf_force=self%lwf_force, &
168 0 : & energy=self%energy, energy_table=energy_table)
169 : !print *, "Calculated energy", self%energy/self%supercell%ncell
170 :
171 0 : do i = 1, self%nstep
172 0 : call self%run_one_mc_step(effpot)
173 : end do
174 : !print *, "Calculated MC energy", self%energy/self%supercell%ncell
175 : !print *, "Accept rate:", (1.0_dp*self%naccept)/self%nattempt
176 0 : lwf(:)=self%lwf(:)
177 0 : if (present(energy_table)) then
178 0 : key = 'LWF energy'
179 0 : call energy_table%put(key, self%energy)
180 : end if
181 :
182 0 : end subroutine run_one_step
183 :
184 :
185 : !----------------------------------------------------------------------
186 : !> @brief accept the trail step, which update the lwf and energy
187 : !----------------------------------------------------------------------
188 0 : subroutine accept(self)
189 : class(lwf_mc_t), intent(inout) :: self
190 0 : self%lwf(self%imove)=self%lwf_new
191 0 : self%energy=self%energy+self%deltaE
192 : !print *, "E:", self%energy/self%supercell%ncell
193 0 : end subroutine accept
194 :
195 : !----------------------------------------------------------------------
196 : !> @brief reject the trail step, changes nothing.
197 : !----------------------------------------------------------------------
198 0 : subroutine reject(self)
199 : class(lwf_mc_t), intent(inout) :: self
200 : ! do nothing.
201 0 : ABI_UNUSED_A(self)
202 0 : end subroutine reject
203 :
204 : !----------------------------------------------------------------------
205 : !> @brief define a trail step using Hinzke_nowak method and calculate energy difference
206 : !----------------------------------------------------------------------
207 0 : function attempt(self,rng, effpot) result(r)
208 : class(lwf_mc_t) :: self
209 : class(rng_t) :: rng
210 : class(abstract_potential_t), intent(inout) :: effpot
211 : real(dp) :: r
212 : ! choose one site
213 0 : self%imove = rng%rand_choice(self%nlwf)
214 0 : self%lwf_old= self%lwf(self%imove)
215 0 : self%deltaE=0.0
216 0 : call move(rng, self%lwf_old, self%lwf_new, self%avg_amp)
217 0 : call effpot%get_delta_E_lwf( self%lwf, self%imove, self%lwf_new, self%deltaE)
218 0 : r=exp(-self%deltaE *self%beta)
219 0 : end function attempt
220 :
221 : !----------------------------------------------------------------------
222 : !> @brief add to lwf by a random value.
223 : !----------------------------------------------------------------------
224 0 : subroutine move(rng, lwf_old, lwf_new, avg_amp)
225 : type(rng_t) :: rng
226 : real(dp), intent(in) :: lwf_old, avg_amp
227 : real(dp), intent(out) :: lwf_new
228 : real(dp):: dlwf
229 0 : dlwf=rng%rand_normal()
230 0 : lwf_new=lwf_old + dlwf*avg_amp
231 0 : end subroutine move
232 :
233 0 : end module m_lwf_mc_mover
|
