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