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