Line data Source code
1 : !!****m* ABINIT/m_spin_mc_mover
2 : !! NAME
3 : !! m_spin_mc_mover
4 : !!
5 : !! FUNCTION
6 : !! This module contains the spin Markov chain Monte Carlo functions for spin mover .
7 : !!
8 : !!
9 : !! Datatypes:
10 : !!
11 : !! * spin_mc_t : MCMC. It defines how to move spins 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_spin_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 : implicit none
41 : !!***
42 : private
43 :
44 : !----------------------------------------------------------------------
45 : !> @brief An helper type to run spin dynamics
46 : ! The Metropolis-Hasting algorithm is used.
47 : !----------------------------------------------------------------------
48 : type,public :: spin_mc_t
49 : real(dp), allocatable :: S(:,:) ! the spin for the whole structure
50 : real(dp) :: Sold(3), Snew(3) ! old and new S for one spin
51 : real(dp) :: angle ! an angle to rotate by average
52 : real(dp) :: energy, deltaE ! energy and the change of energy when change one spin
53 : real(dp) :: temperature
54 : real(dp) :: beta ! 1/(kb T)
55 : integer :: nspin ! number of spins
56 : integer :: nstep ! number of steps
57 : integer :: imove ! index of spin to be moved
58 : integer :: naccept ! number of accepted steps
59 : integer :: nattempt ! number of attempted steps
60 : contains
61 : procedure :: initialize
62 : procedure :: finalize
63 : procedure, private :: attempt
64 : procedure, private :: accept
65 : procedure, private :: reject
66 : procedure, private :: run_one_step
67 : procedure :: run_MC
68 : end type spin_mc_t
69 :
70 :
71 : contains
72 : !----------------------------------------------------------------------
73 : !> @brief initialize mc helper class
74 : !>
75 : !> @param[in] nspin: number of spins
76 : !> @param[in] angle: a angle to rotate
77 : !> @param[in] temperature: temperature
78 : !----------------------------------------------------------------------
79 0 : subroutine initialize(self, nspin, angle, temperature)
80 : class(spin_mc_t), intent(inout) :: self
81 : integer, intent(in) :: nspin
82 : real(dp), intent(in) :: angle, temperature
83 0 : self%nspin=nspin
84 0 : self%nstep=self%nspin
85 0 : ABI_MALLOC(self%S, (3, self%nspin))
86 0 : self%angle=angle
87 0 : self%temperature=temperature
88 0 : self%beta=1.0/temperature ! Kb in a.u. is 1.
89 0 : self%Sold(:)=0.0_dp
90 0 : self%Snew(:)=0.0_dp
91 0 : self%naccept=0
92 0 : self%nattempt=0
93 0 : end subroutine initialize
94 :
95 : !----------------------------------------------------------------------
96 : !> @brief finalize
97 : !----------------------------------------------------------------------
98 0 : subroutine finalize(self)
99 : class(spin_mc_t), intent(inout) :: self
100 0 : if (allocated(self%S)) then
101 0 : ABI_FREE(self%S)
102 : end if
103 0 : self%Sold=zero
104 0 : self%Snew=zero
105 0 : self%nspin=0
106 0 : self%nstep=0
107 0 : end subroutine finalize
108 :
109 :
110 : !----------------------------------------------------------------------
111 : !> @brief run one monte carlo step
112 : !> @param[in] rngL rundom number generator
113 : !> @param[in] effpot: effective spin potential
114 : !----------------------------------------------------------------------
115 0 : subroutine run_one_step(self, rng, effpot)
116 : class(spin_mc_t) :: self
117 : class(rng_t) :: rng
118 : class(abstract_potential_t), intent(inout) :: effpot
119 : real(dp) :: r
120 :
121 : ! try to change spin
122 0 : r=self%attempt(rng, effpot)
123 : ! metropolis-hastings
124 0 : self%nattempt = self%nattempt+1
125 0 : if(rng%rand_unif_01()< min(1.0_dp, r) ) then
126 0 : self%naccept=self%naccept+1
127 0 : call self%accept()
128 : !print *, "accepted"
129 : else
130 0 : call self%reject()
131 : !print *, "rejected"
132 : end if
133 0 : end subroutine run_one_step
134 :
135 : !----------------------------------------------------------------------
136 : !> @brief run a number of MC steps. Since one step only changes
137 : !> too little things, a few steps are bunched as one. Then things like
138 : !> output or calculation of observables are done after the big step.
139 : !>
140 : !> @param[in] rng: random number generator
141 : !> @param[in] effpot: the spin potential
142 : !> @param[in] S_in: the intial spin state
143 : !> @param[out] etot: the final total energy
144 : !----------------------------------------------------------------------
145 0 : subroutine run_MC(self, rng, effpot, S_in, etot, bfield)
146 : class(spin_mc_t), intent(inout) :: self
147 : type(rng_t) :: rng
148 : class(abstract_potential_t), intent(inout) :: effpot
149 : real(dp), intent(inout) :: S_in(3,self%nspin)
150 : real(dp), intent(out) :: etot
151 : real(dp), optional,intent(inout) :: bfield(:,:)
152 : real(dp) :: etmp
153 :
154 : integer :: i
155 0 : self%S(:,:)=S_in(:,:)
156 0 : call effpot%calculate(spin=S_in, energy=self%energy, bfield=bfield)
157 0 : do i = 1, self%nstep
158 0 : call self%run_one_step(rng, effpot)
159 : end do
160 0 : S_in(:, :)=self%S(:,:)
161 : !call effpot%calculate(spin=self%S, energy=self%energy, bfield=bfield)
162 : !print *, self%energy
163 0 : etot=self%energy
164 0 : call effpot%calculate(spin=S_in, energy=etmp, bfield=bfield)
165 : !print *, "energy: ", self%energy, etmp, self%energy-etmp
166 0 : end subroutine run_MC
167 :
168 : !----------------------------------------------------------------------
169 : !> @brief accept the trail step, which update the spin and energy
170 : !----------------------------------------------------------------------
171 0 : subroutine accept(self)
172 : class(spin_mc_t), intent(inout) :: self
173 0 : self%S(:,self%imove)=self%Snew(:)
174 0 : self%energy=self%energy+self%deltaE
175 0 : end subroutine accept
176 :
177 : !----------------------------------------------------------------------
178 : !> @brief reject the trail step, changes nothing.
179 : !----------------------------------------------------------------------
180 0 : subroutine reject(self)
181 : class(spin_mc_t), intent(inout) :: self
182 : ! do nothing.
183 0 : ABI_UNUSED_A(self)
184 0 : end subroutine reject
185 :
186 : !----------------------------------------------------------------------
187 : !> @brief define a trail step using Hinzke_nowak method and calculate energy difference
188 : !----------------------------------------------------------------------
189 0 : function attempt(self,rng, effpot) result(r)
190 : class(spin_mc_t) :: self
191 : class(rng_t) :: rng
192 : class(abstract_potential_t), intent(inout) :: effpot
193 : real(dp) :: r
194 : ! choose one site
195 0 : self%imove = rng%rand_choice(self%nspin)
196 0 : self%Sold(:)= self%S(:,self%imove)
197 0 : self%deltaE=0.0
198 0 : call move_hinzke_nowak(rng, self%Sold, self%Snew, self%angle)
199 0 : call effpot%get_delta_E( self%S, self%imove, self%Snew, self%deltaE)
200 : !print *, "delta E", self%deltaE
201 0 : r=exp(-self%deltaE *self%beta)
202 0 : end function attempt
203 :
204 : !----------------------------------------------------------------------
205 : !> @brief rotate the spin by the average of angle (normal distribution)
206 : !----------------------------------------------------------------------
207 0 : subroutine move_angle(rng, Sold, Snew, angle)
208 : type(rng_t) :: rng
209 : real(dp), intent(in) :: Sold(3), angle
210 : real(dp), intent(out) :: Snew(3)
211 0 : call rng%rand_normal_array(Snew, 3)
212 0 : Snew(:)=Sold(:) + Snew(:)*angle
213 0 : Snew(:)=Snew(:)/norm2(Snew)
214 0 : end subroutine move_angle
215 :
216 : !----------------------------------------------------------------------
217 : !> @brief flip one spin
218 : !----------------------------------------------------------------------
219 : subroutine move_flip(Sold, Snew)
220 : real(dp), intent(in) :: Sold(3)
221 : real(dp), intent(out) :: Snew(3)
222 0 : Snew(:)=-Sold(:)
223 : end subroutine move_flip
224 :
225 : !----------------------------------------------------------------------
226 : !> @brief set spin to random orientation
227 : !----------------------------------------------------------------------
228 0 : subroutine move_uniform(rng, Snew)
229 : type(rng_t), intent(inout) :: rng
230 : real(dp), intent(out) :: Snew(3)
231 0 : call rng%rand_normal_array(Snew, 3)
232 0 : Snew(:)=Snew(:)/norm2(Snew)
233 0 : end subroutine move_uniform
234 :
235 : !----------------------------------------------------------------------
236 : !> @brief combine rotate, flip and random set.
237 : !----------------------------------------------------------------------
238 0 : subroutine move_hinzke_nowak(rng, Sold, Snew, angle)
239 : type(rng_t), intent(inout) :: rng
240 : real(dp), intent(in) :: Sold(3), angle
241 : real(dp), intent(out) :: Snew(3)
242 : integer :: move
243 0 : move=rng%rand_choice(3)
244 0 : select case (move)
245 : case (1)
246 0 : call move_angle(rng, Sold, Snew, angle)
247 : case(2)
248 0 : call move_flip(Sold, Snew)
249 : case(3)
250 0 : call move_uniform(rng, Snew)
251 : case default
252 0 : call move_angle(rng, Sold, Snew, angle)
253 : end select
254 0 : end subroutine move_hinzke_nowak
255 :
256 0 : end module m_spin_mc_mover
|
