Branch data Line data Source code
1 : : !!****m* ABINIT/m_spin_mover
2 : : !! NAME
3 : : !! m_spin_mover
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the spin mover, which controls how the spin
7 : : !!
8 : : !!
9 : : !! Datatypes:
10 : : !!
11 : : !! * spin_mover_t
12 : : !!
13 : : !! Subroutines:
14 : : !!
15 : : !! * spin_mover_t_initialize
16 : : !! * spin_mover_t_run_one_step
17 : : !! * spin_mover_t_run_time
18 : : !! * TODO: update this when F2003 documentation format decided.
19 : : !!
20 : : !!
21 : : !! COPYRIGHT
22 : : !! Copyright (C) 2001-2026 ABINIT group (hexu)
23 : : !! This file is distributed under the terms of the
24 : : !! GNU General Public License, see ~abinit/COPYING
25 : : !! or http://www.gnu.org/copyleft/gpl.txt .
26 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
27 : : !!
28 : : !! SOURCE
29 : :
30 : :
31 : : #if defined HAVE_CONFIG_H
32 : : #include "config.h"
33 : : #endif
34 : :
35 : : #include "abi_common.h"
36 : :
37 : :
38 : : module m_spin_mover
39 : :
40 : : use defs_basis
41 : : use m_nctk
42 : : #if defined HAVE_NETCDF
43 : : use netcdf
44 : : #endif
45 : : use m_errors
46 : : use m_abicore
47 : : use m_xmpi
48 : : use m_io_tools, only : get_unit, open_file, close_unit
49 : : use m_mpi_scheduler, only: mpi_scheduler_t, init_mpi_info
50 : : use m_mathfuncs, only : cross
51 : : use m_spin_observables , only : spin_observable_t
52 : : use m_spin_potential, only: spin_potential_t
53 : : use m_spin_hist, only: spin_hist_t
54 : : use m_spin_ncfile, only: spin_ncfile_t
55 : : use m_multibinit_dataset, only: multibinit_dtset_type
56 : : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
57 : : use m_random_xoroshiro128plus, only: set_seed, rand_normal_array, rng_t
58 : : use m_abstract_potential, only: abstract_potential_t
59 : : use m_abstract_mover, only: abstract_mover_t
60 : : use m_spin_mc_mover, only : spin_mc_t
61 : : use m_hashtable_strval, only: hash_table_t
62 : : implicit none
63 : : private
64 : : !!***
65 : :
66 : :
67 : : !!****t* m_spin_mover/spin_mover_t
68 : : !! NAME
69 : : !! spin_mover_t
70 : : !!
71 : : !! FUNCTION
72 : : !! this type contains the parameters for the spin mover.
73 : : !!
74 : : !! It contains:
75 : : !! dt: time step
76 : : !! total_time
77 : : !! temperature.
78 : : !! nspin number of magnetic atoms
79 : : !! SOURCE
80 : :
81 : :
82 : : type, public, extends(abstract_mover_t) :: spin_mover_t
83 : : integer :: nspin, method=0
84 : : real(dp), allocatable :: gyro_ratio(:), damping(:), gamma_L(:), H_lang_coeff(:), ms(:), Stmp(:,:), Stmp2(:,:)
85 : : real(dp), allocatable :: Heff_tmp(:,:), Htmp(:,:), Hrotate(:,:), H_lang(:,:), buffer(:,:)
86 : : real(dp) :: init_qpoint(3), init_rotate_axis(3) ! qpoint and rotation axis to set up initial spin configuration
87 : : real(dp) :: init_orientation(3) ! spin orientation in primitive cell which is then rotated
88 : : type(spin_hist_t) :: hist
89 : : logical :: gamma_l_calculated
90 : : type(spin_mc_t) :: spin_mc
91 : : type(mpi_scheduler_t) :: mps
92 : : type(spin_observable_t) :: spin_ob
93 : : type(spin_ncfile_t) :: spin_ncfile
94 : : type(multibinit_dtset_type), pointer :: params
95 : : CONTAINS
96 : : procedure :: initialize
97 : : procedure :: finalize
98 : : procedure :: set_initial_state
99 : : procedure :: read_hist_spin_state
100 : : procedure, private :: run_one_step_DM => spin_mover_t_run_one_step_DM
101 : : procedure, private :: run_one_step_HeunP => spin_mover_t_run_one_step_HeunP
102 : : procedure, private :: run_one_step_dummy=> spin_mover_t_run_one_step_dummy
103 : : procedure, private :: run_one_step_MC=> spin_mover_t_run_one_step_MC
104 : : procedure :: run_one_step => spin_mover_t_run_one_step
105 : : procedure :: run_time => spin_mover_t_run_time
106 : : procedure :: run_MvT
107 : : procedure :: set_temperature
108 : : procedure, private :: prepare_ncfile
109 : : procedure, private ::get_Langevin_Heff
110 : : procedure :: current_spin
111 : : procedure :: set_ncfile_name
112 : : end type spin_mover_t
113 : : !!***
114 : :
115 : : contains
116 : :
117 : : !!****f* m_spin_mover/initialize
118 : : !!
119 : : !! NAME
120 : : !! initialize
121 : : !!
122 : : !! FUNCTION
123 : : !! initialize the spin mover
124 : : !!
125 : : !! INPUTS
126 : : !!
127 : : !! OUTPUT
128 : : !!
129 : : !! NOTES
130 : : !!
131 : : !! SOURCE
132 : 2 : subroutine initialize(self, params, supercell, rng)
133 : : class(spin_mover_t), intent(inout) :: self
134 : : type(multibinit_dtset_type), target :: params
135 : : type(mbsupercell_t), target :: supercell
136 : : type(rng_t), target, intent(in) :: rng
137 : : integer :: nspin
138 : :
139 : : integer :: master, my_rank, comm, nproc, ierr
140 : : logical :: iam_master
141 : 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
142 : :
143 : 2 : self%params=>params
144 : 2 : self%supercell=>supercell
145 : 2 : self%rng => rng
146 [ + - ]: 2 : if (iam_master) then
147 : 2 : nspin=supercell%spin%nspin
148 : 2 : self%nspin=nspin
149 : 2 : self%dt= params%spin_dt
150 : 2 : self%thermal_time= params%spin_ntime_pre * self%dt
151 : 2 : self%total_time= params%spin_ntime * self%dt
152 : 2 : self%temperature=params%spin_temperature
153 [ + - ]: 2 : if(params%spin_dynamics>=0) then
154 : 2 : self%method=params%spin_dynamics
155 : : endif
156 [ - + ]: 2 : if(params%spin_init_state==3) then
157 [ # # ]: 0 : self%init_qpoint = params%spin_init_qpoint
158 [ # # ]: 0 : self%init_rotate_axis = params%spin_init_rotate_axis
159 [ # # ]: 0 : self%init_orientation = params%spin_init_orientation
160 : : endif
161 : : end if
162 [ - + ]: 2 : if(params%spin_dynamics==3) then ! Monte carlo
163 : 0 : call self%spin_mc%initialize(nspin=nspin, angle=0.2_dp, temperature=params%spin_temperature)
164 : : end if
165 : 2 : call xmpi_bcast(self%nspin, master, comm, ierr)
166 : 2 : call xmpi_bcast(self%dt, master, comm, ierr)
167 : 2 : call xmpi_bcast(self%thermal_time, master, comm, ierr)
168 : 2 : call xmpi_bcast(self%total_time, master, comm, ierr)
169 : 2 : call xmpi_bcast(self%temperature, master, comm, ierr)
170 : 2 : call xmpi_bcast(self%method, master, comm, ierr)
171 : :
172 [ - + + - : 6 : ABI_MALLOC(self%ms, (self%nspin) )
+ - + - -
+ - + -
+ ]
173 [ - + + - : 6 : ABI_MALLOC(self%gyro_ratio, (self%nspin) )
+ - + - -
+ - + -
+ ]
174 [ - + + - : 6 : ABI_MALLOC(self%damping, (self%nspin) )
+ - + - -
+ - + -
+ ]
175 [ - + + - : 6 : ABI_MALLOC(self%gamma_l, (self%nspin) )
+ - + - -
+ - + -
+ ]
176 [ - + + - : 6 : ABI_MALLOC(self%H_lang_coeff, (self%nspin) )
+ - + - -
+ - + -
+ ]
177 : :
178 [ - + + - : 6 : ABI_MALLOC(self%Heff_tmp, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
179 [ - + + - : 6 : ABI_MALLOC(self%Htmp, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
180 [ - + + - : 6 : ABI_MALLOC(self%Hrotate, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
181 [ - + + - : 6 : ABI_MALLOC(self%Stmp, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
182 [ - + + - : 6 : ABI_MALLOC(self%Stmp2, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
183 [ - + + - : 6 : ABI_MALLOC(self%buffer, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
184 [ - + + - : 6 : ABI_MALLOC(self%H_lang, (3,self%nspin) )
+ - + - -
+ - + -
+ ]
185 : :
186 : 2 : self%gamma_l_calculated=.False.
187 : 2 : call self%mps%initialize(ntasks=nspin,master=master, comm=comm)
188 : :
189 : :
190 : 2 : call xmpi_bcast(params%spin_damping, master, comm, ierr)
191 : :
192 [ + - ]: 2 : if (iam_master) then
193 [ + - ]: 2 : if (params%spin_damping >=0) then
194 [ + + ]: 434 : self%damping(:)= params%spin_damping
195 : : else
196 [ # # ]: 0 : self%damping(:)=supercell%spin%gilbert_damping(:)
197 : : end if
198 : :
199 [ + + ]: 434 : self%gyro_ratio(:)=supercell%spin%gyro_ratio(:)
200 [ + + ]: 434 : self%ms(:)=supercell%spin%ms(:)
201 : : endif
202 : :
203 : 2 : call xmpi_bcast(self%damping, master, comm, ierr)
204 : 2 : call xmpi_bcast(self%gyro_ratio, master, comm, ierr)
205 : 2 : call xmpi_bcast(self%ms, master, comm, ierr)
206 : 2 : call self%set_temperature(temperature=params%spin_temperature)
207 : :
208 : : ! Hist and set initial spin state
209 [ + - ]: 2 : if(iam_master) then
210 : : call self%hist%initialize(nspin=self%nspin, &
211 : 2 : & mxhist=3, has_latt=.False.)
212 : : call self%hist%set_params(spin_nctime=params%spin_nctime, &
213 : 2 : & spin_temperature=params%spin_temperature)
214 : : endif
215 : :
216 : 2 : call self%set_initial_state(mode=params%spin_init_state)
217 : :
218 : : ! observable
219 [ + - ]: 2 : if(iam_master) then
220 : 2 : call self%spin_ob%initialize(self%supercell, params)
221 : : endif
222 : :
223 : 6 : end subroutine initialize
224 : : !!***
225 : :
226 : : !-------------------------------------------------------------------!
227 : : !set_ncfile_name :
228 : : !-------------------------------------------------------------------!
229 : 2 : subroutine set_ncfile_name(self, params, fname)
230 : : class(spin_mover_t), intent(inout) :: self
231 : : type(multibinit_dtset_type) :: params
232 : : character(len=fnlen), intent(in) :: fname
233 : : integer :: master, my_rank, comm, nproc
234 : : logical :: iam_master
235 : 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
236 [ + - ]: 2 : if (iam_master) then
237 [ + - ]: 2 : call self%prepare_ncfile(params, trim(fname)//'_spinhist.nc')
238 : 2 : call self%spin_ncfile%write_one_step(self%hist)
239 : : endif
240 : 2 : end subroutine set_ncfile_name
241 : :
242 : :
243 : : !-------------------------------------------------------------------!
244 : : ! read_hist_spin_state
245 : : ! read the last step of spin from hist file.
246 : : ! and save if to self%Stmp
247 : : !-------------------------------------------------------------------!
248 : 0 : subroutine read_hist_spin_state(self, fname)
249 : : class(spin_mover_t), intent(inout) :: self
250 : : character(len=fnlen), intent(in) :: fname
251 : : integer :: ierr, ncid, varid
252 : : integer :: nspin, ntime
253 : : character(len=500) :: msg
254 : : ! open file
255 : :
256 : : #if defined HAVE_NETCDF
257 : 0 : ierr=nf90_open(trim(fname), NF90_NOWRITE, ncid)
258 [ # # # # ]: 0 : NCF_CHECK_MSG(ierr, "The spin_init_mode is set to 4. But opening netcdf file "//trim(fname)//" Failed. ")
259 : :
260 : : ! sanity check. If the hist file is consistent with the current calculation
261 : 0 : ierr=nctk_get_dim(ncid, "nspin" , nspin)
262 [ # # ]: 0 : NCF_CHECK_MSG(ierr, "when reading nspin")
263 : :
264 : : msg="The number of spins in histfile is not equal & & to the present calculation." // &
265 : 0 : & " Please check if the file is consistent."
266 [ # # ]: 0 : if (nspin /= self%nspin) then
267 : 0 : ABI_ERROR(msg)
268 : : end if
269 : :
270 : :
271 : 0 : ierr=nctk_get_dim(ncid, "ntime", ntime)
272 [ # # ]: 0 : NCF_CHECK_MSG(ierr, "when reading ntime")
273 : :
274 : :
275 : : ! TODO: more check ???
276 : :
277 : : ! read Spin and set as initial state
278 : 0 : ierr =nf90_inq_varid(ncid, "S", varid)
279 : : msg="when reading S. Try using spin_init_state=3 option instead (specify spin_init_qpoint," // &
280 : 0 : & " spin_init_rotate_axis and spin_init_orientation as needed)."
281 [ # # ]: 0 : NCF_CHECK_MSG(ierr, msg)
282 : :
283 [ # # # # ]: 0 : ierr = nf90_get_var(ncid=ncid, varid=varid, values=self%Stmp(:,:), start=(/1, 1, ntime/), count=(/3, nspin,1/))
284 [ # # ]: 0 : NCF_CHECK_MSG(ierr, "when reading S from spin hist file")
285 : :
286 : : ! close file
287 : 0 : ierr=nf90_close(ncid)
288 [ # # ]: 0 : NCF_CHECK_MSG(ierr, "Close netcdf file")
289 : : #else
290 : : ABI_ERROR("spin_init_state set to 4 but abinit is not compiled with netcdf.")
291 : : #endif
292 : :
293 : 0 : end subroutine read_hist_spin_state
294 : :
295 : : !----------------------------------------------------------------------------!
296 : : !set_initial_state:
297 : : ! mode: which configuration to use
298 : : ! 1. Random
299 : : ! 2. reference state from potential file
300 : : ! 3. spin configuration using qpoint and rotation axis (e.g. for FM or AFM)
301 : : ! 4. Restart from last entry of hist netcdf file
302 : : !----------------------------------------------------------------------------!
303 : 2 : subroutine set_initial_state(self, mode)
304 : : class(spin_mover_t), intent(inout) :: self
305 : : integer, optional, intent(in) :: mode
306 : :
307 : : integer :: i, init_mode
308 : : character(len=500) :: msg
309 : :
310 : : integer :: master, my_rank, comm, nproc, ierr
311 : : logical :: iam_master
312 : 2 : real(dp), allocatable :: Sprim(:,:)
313 : :
314 : 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
315 : :
316 [ + - ]: 2 : if(iam_master) then
317 [ + - ]: 2 : if (present(mode)) then
318 : 2 : init_mode=mode
319 : : else
320 : : init_mode=1
321 : : end if
322 : :
323 : 2 : select case (init_mode)
324 : : case (1)
325 : : ! randomize S using uniform random number
326 : 2 : write(msg,*) "Initial spins set to random values."
327 : 2 : call wrtout(ab_out,msg,'COLL')
328 : 2 : call wrtout(std_out,msg,'COLL')
329 : : !call random_number(self%Stmp)
330 : 2 : call self%rng%rand_unif_01_array(self%Stmp, self%nspin*3 )
331 [ + + + + ]: 1730 : self%Stmp=self%Stmp-0.5
332 [ + + ]: 434 : do i=1, self%nspin
333 [ + + + + ]: 3026 : self%Stmp(:,i)=self%Stmp(:,i)/sqrt(sum(self%Stmp(:, i)**2))
334 : : end do
335 : :
336 : : case (2)
337 : : ! set spin to reference state using the reference qpoint and rotation axis from potential file
338 : 0 : write(msg,*) "Initial spins set to reference configuration."
339 : 0 : call wrtout(ab_out,msg,'COLL')
340 : 0 : call wrtout(std_out,msg,'COLL')
341 : :
342 [ # # ]: 0 : do i=1, self%nspin
343 [ # # # # ]: 0 : self%Stmp(:,:) = self%supercell%spin%Sref(:,:)
344 : : end do
345 : :
346 : : case (3)
347 : 0 : write(msg,*) "Initial spins set according to spin_init_* variables."
348 : 0 : call wrtout(ab_out,msg,'COLL')
349 : 0 : call wrtout(std_out,msg,'COLL')
350 : :
351 [ # # # # : 0 : ABI_MALLOC(Sprim, (3,self%supercell%unitcell%spin%nspin) )
# # # # #
# # # ]
352 : :
353 : : ! set inital spin state using the input variables
354 : : ! set spin to ferromagnetic along init_orientation then rotate
355 [ # # ]: 0 : do i=1, self%supercell%unitcell%spin%nspin
356 [ # # ]: 0 : Sprim(:,i)=self%init_orientation(:)
357 : : enddo
358 [ # # # # ]: 0 : self%Stmp(:,:) = 0.0d0
359 : :
360 : : call self%supercell%supercell_maker%generate_spin_wave_vectorlist(A=Sprim, &
361 : 0 : & kpoint=self%init_qpoint, axis=self%init_rotate_axis, A_sc=self%Stmp)
362 : :
363 [ # # ]: 0 : ABI_SFREE(SPrim)
364 : :
365 : : case (4)
366 : : ! read from last step of hist file
367 : 0 : write(msg,'(a,a,a)') "Initial spins set to input spin hist file ",&
368 : 0 : & trim(self%params%spin_init_hist_fname), '.'
369 : 0 : call wrtout(ab_out,msg,'COLL')
370 : 0 : call wrtout(std_out,msg,'COLL')
371 [ + - - - : 2 : call self%read_hist_spin_state(fname=self%params%spin_init_hist_fname)
- ]
372 : :
373 : : end select
374 : :
375 : : call self%hist%set_vars(S=self%Stmp, Snorm=self%supercell%spin%ms, &
376 : 2 : & time=0.0_dp, ihist_latt=0, inc=.True.)
377 : :
378 : : endif
379 : 2 : call xmpi_bcast(self%Stmp, 0, comm, ierr)
380 [ - + ]: 2 : end subroutine set_initial_state
381 : :
382 : :
383 : : !-------------------------------------------------------------------!
384 : : ! prepare_ncfile:
385 : : !-------------------------------------------------------------------!
386 : 2 : subroutine prepare_ncfile(self, params, fname)
387 : : class(spin_mover_t), intent(inout) :: self
388 : : type(multibinit_dtset_type) :: params
389 : : character(len=*), intent(in) :: fname
390 : :
391 : : integer :: master, my_rank, comm, nproc
392 : : logical :: iam_master
393 : 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
394 : :
395 [ + - ]: 2 : if(iam_master) then
396 [ + - ]: 2 : call self%spin_ncfile%initialize( trim(fname), params%spin_write_traj)
397 : 2 : call self%spin_ncfile%def_spindynamics_var(self%hist)
398 : 2 : call self%spin_ncfile%def_observable_var(self%spin_ob)
399 : 2 : call self%spin_ncfile%write_primitive_cell(self%supercell%unitcell)
400 : 2 : call self%spin_ncfile%write_supercell(self%supercell)
401 : 2 : call self%spin_ncfile%write_parameters(params)
402 : : endif
403 : 2 : end subroutine prepare_ncfile
404 : :
405 : :
406 : :
407 : 4 : subroutine set_temperature(self, temperature)
408 : : class(spin_mover_t), intent(inout) :: self
409 : : real(dp), optional, intent(in) :: temperature
410 : :
411 : : integer :: master, my_rank, comm, nproc, ierr
412 : : logical :: iam_master
413 : 2 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
414 : :
415 [ + - ]: 2 : if(present(temperature)) self%temperature=temperature
416 : 2 : call xmpi_bcast(self%temperature, master, comm, ierr)
417 [ - + ]: 2 : if(self%method==3) then
418 [ # # ]: 0 : if(iam_master) self%spin_mc%temperature = temperature
419 [ # # ]: 0 : if(iam_master) self%spin_mc%beta=1.0_dp/temperature
420 : 0 : call xmpi_bcast(self%spin_mc%temperature, master, comm, ierr)
421 : 0 : call xmpi_bcast(self%spin_mc%beta, master, comm, ierr)
422 : : end if
423 [ + + ]: 434 : self%gamma_l(:)= self%gyro_ratio(:)/(1.0_dp+ self%damping(:)**2)
424 : 2 : self%gamma_l_calculated=.True.
425 : : self%H_lang_coeff(:)=sqrt(2.0*self%damping(:)* self%temperature &
426 [ + + ]: 434 : & /(self%gyro_ratio(:)* self%dt *self%ms(:)))
427 : 2 : end subroutine set_temperature
428 : :
429 : :
430 : 2002 : subroutine get_Langevin_Heff(self, H_lang)
431 : : class(spin_mover_t), intent(inout) :: self
432 : : real(dp), intent(inout):: H_lang(3,self%nspin)
433 : : integer :: i
434 [ - + ]: 2002 : if ( self%temperature .gt. 1d-7) then
435 : 2002 : call rand_normal_array(self%rng, H_lang(:, self%mps%istart:self%mps%iend), 3*self%mps%ntask)
436 [ + + ]: 434434 : do i = self%mps%istart, self%mps%iend
437 [ + + ]: 1731730 : H_lang(:,i)= H_lang(:,i) * self%H_lang_coeff(i)
438 : : end do
439 : : else
440 [ # # # # ]: 0 : H_lang(:,:)=0.0_dp
441 : : end if
442 : 2002 : end subroutine get_Langevin_Heff
443 : :
444 : :
445 : : !!****f* m_spin_mover/spin_mover_t_run_one_step_HeunP
446 : : !!
447 : : !! NAME
448 : : !! spin_mover_t_run_one_step_HeunP
449 : : !!
450 : : !! FUNCTION
451 : : !! run one spin step using HeunP method
452 : : !!
453 : : !! INPUTS
454 : : !! effpot: abstract_potential_t type.
455 : : !! S_in : input spin. (3*nspin)
456 : : !!
457 : : !! OUTPUT
458 : : !! etot: energy (scalar)
459 : : !!
460 : : !! SOURCE
461 : 0 : subroutine spin_mover_t_run_one_step_HeunP(self, effpot, S_in, &
462 [ # # # # : 0 : & etot, displacement, strain, lwf, energy_table)
# # # # #
# # # # #
# # # # ]
463 : : !class (spin_mover_t), intent(inout):: self
464 : : class(spin_mover_t), intent(inout):: self
465 : : class(abstract_potential_t), intent(inout) :: effpot
466 : :
467 : : real(dp), optional, intent(inout):: displacement(:,:), &
468 : : strain(:,:), lwf(:)
469 : : real(dp), intent(inout) :: S_in(3,self%nspin)
470 : : real(dp), intent(out) :: etot
471 : : integer :: i
472 : : real(dp) :: dSdt(3), Htmp(3), Ri(3)
473 : : type(hash_table_t),optional, intent(inout) :: energy_table
474 : :
475 : : !integer :: master, my_rank, comm, nproc, ierr
476 : : !logical :: iam_master
477 : : !call init_mpi_info(master, iam_master, my_rank, comm, nproc)
478 : :
479 : : ! predict
480 : 0 : etot=0.0
481 [ # # # # ]: 0 : self%Heff_tmp(:,:)=0.0
482 : : call effpot%calculate(displacement=displacement, strain=strain, lwf=lwf, spin=S_in, &
483 [ # # # # : 0 : & bfield=self%Heff_tmp, energy=etot, energy_table=energy_table)
# # ]
484 : 0 : call self%get_Langevin_Heff(self%H_lang)
485 [ # # ]: 0 : do i=self%mps%istart, self%mps%iend
486 [ # # ]: 0 : Htmp=self%Heff_tmp(:,i)+self%H_lang(:,i)
487 : 0 : Ri = cross(S_in(:,i),Htmp)
488 [ # # ]: 0 : dSdt = -self%gamma_L(i)*(Ri+self%damping(i)* cross(S_in(:,i), Ri))
489 [ # # ]: 0 : Ri=S_in(:,i)+dSdt*self%dt
490 [ # # # # : 0 : Ri=Ri/norm2(Ri)
# # # # ]
491 [ # # ]: 0 : self%Stmp2(:,i)=Ri
492 : : end do
493 : 0 : call self%mps%allgatherv_dp2d(self%Stmp2, 3, buffer=self%buffer)
494 : :
495 : : ! correction
496 [ # # # # ]: 0 : self%Htmp(:,:)=0.0
497 : 0 : etot=0.0
498 : : call effpot%calculate(displacement=displacement, strain=strain, lwf=lwf,spin=self%Stmp2, &
499 [ # # # # : 0 : & bfield=self%Htmp, energy=etot, energy_table=energy_table)
# # ]
500 [ # # ]: 0 : do i=self%mps%istart, self%mps%iend
501 [ # # ]: 0 : Htmp=(self%Heff_tmp(:,i)+self%Htmp(:,i))*0.5_dp+self%H_lang(:,i)
502 : 0 : Ri = cross(S_in(:,i),Htmp)
503 [ # # ]: 0 : dSdt = -self%gamma_L(i)*(Ri+self%damping(i)* cross(S_in(:,i), Ri))
504 [ # # ]: 0 : Ri=S_in(:,i)+dSdt*self%dt
505 [ # # # # : 0 : Ri=Ri/norm2(Ri)
# # # # ]
506 [ # # ]: 0 : self%Stmp(:,i)=Ri
507 : : end do
508 : 0 : call self%mps%allgatherv_dp2d(self%Stmp, 3, buffer=self%buffer)
509 : 0 : end subroutine spin_mover_t_run_one_step_HeunP
510 : : !!***
511 : :
512 : :
513 : :
514 : : !!****f* m_spin_mover/spin_mover_t_run_one_step_dummy
515 : : !!
516 : : !! NAME
517 : : !! spin_mover_t_run_one_step_dummy
518 : : !!
519 : : !! FUNCTION
520 : : !! run one spin step using dummy method
521 : : !!
522 : : !! INPUTS
523 : : !! effpot: abstract_potential_t type.
524 : : !! S_in : input spin. (3*nspin)
525 : : !!
526 : : !! OUTPUT
527 : : !! etot: energy (scalar)
528 : : !!
529 : : !! SOURCE
530 : 0 : subroutine spin_mover_t_run_one_step_dummy(self, effpot, S_in, etot, &
531 [ # # # # : 0 : & displacement, strain, lwf, energy_table)
# # # # #
# # # # #
# # # # ]
532 : : !class (spin_mover_t), intent(inout):: self
533 : : class(spin_mover_t), intent(inout):: self
534 : : class(abstract_potential_t), intent(inout) :: effpot
535 : :
536 : : real(dp), optional, intent(inout):: displacement(:,:), &
537 : : strain(:,:), lwf(:)
538 : : real(dp), intent(inout) :: S_in(3,self%nspin)
539 : : real(dp), intent(out) :: etot
540 : : type(hash_table_t),optional, intent(inout) :: energy_table
541 : : integer :: i
542 : : real(dp) :: Htmp(3), Ri(3)
543 : :
544 : : ! predict
545 : 0 : etot=0.0
546 [ # # # # ]: 0 : self%Heff_tmp(:,:)=0.0
547 : : call effpot%calculate(displacement=displacement, strain=strain, lwf=lwf, spin=S_in, &
548 [ # # # # : 0 : & bfield=self%Heff_tmp, energy=etot, energy_table=energy_table)
# # ]
549 : 0 : call self%get_Langevin_Heff(self%H_lang)
550 [ # # ]: 0 : do i=self%mps%istart, self%mps%iend
551 : 0 : Htmp=self%Heff_tmp(:,i)+self%H_lang(:,i)
552 : : !Ri = cross(S_in(:,i),Htmp)
553 : : !dSdt = -self%gamma_L(i)*(Ri+self%damping(i)* cross(S_in(:,i), Ri))
554 [ # # ]: 0 : Ri=S_in(:,i)!+dSdt*self%dt
555 [ # # # # : 0 : Ri=Ri/norm2(Ri)
# # # # ]
556 [ # # ]: 0 : self%Stmp(:,i)=Ri
557 : : end do
558 : 0 : call self%mps%allgatherv_dp2d(self%Stmp2, 3, buffer=self%buffer)
559 : :
560 : 0 : end subroutine spin_mover_t_run_one_step_dummy
561 : : !!***
562 : :
563 : :
564 : :
565 : : !----------------------------------------------------------------------
566 : : !> @brief rotate spin with a rotation matrix
567 : : !>
568 : : !> @param[in] S_in: input spin array(3)
569 : : !> @param[in] Heff: effective field
570 : : !> @param[in] dt: time step
571 : : !> @param[out] S_out: output spin
572 : : !----------------------------------------------------------------------
573 : :
574 [ - + ]: 864864 : pure function rotate_S_DM(S_in, Heff, dt) result(S_out)
575 : : ! Depondt & Mertens method to rotate S_in
576 : : real(dp), intent(in) :: S_in(3), Heff(3), dt
577 : : real(dp) :: S_out(3)
578 : : real(dp) :: B(3) , w, u, Bnorm, R(3,3), cosw, sinw
579 [ + + ]: 3459456 : Bnorm=sqrt(sum(Heff*Heff))
580 [ + + ]: 3459456 : B(:)=Heff(:)/Bnorm ! axis of rotation
581 : 864864 : w=Bnorm*dt ! amplitude of rotation
582 : 864864 : sinw=sin(w)
583 : 864864 : cosw=cos(w)
584 : 864864 : u=1.0d0-cosw
585 : :
586 : : ! R is rotation matrix
587 : 864864 : R(1,1)=B(1)*B(1)*u+cosw
588 : 864864 : R(2,1)=B(1)*B(2)*u+B(3)*sinw
589 : 864864 : R(3,1)=B(1)*B(3)*u-B(2)*sinw
590 : :
591 : 864864 : R(1,2)=B(1)*B(2)*u-B(3)*sinw
592 : 864864 : R(2,2)=B(2)*B(2)*u+cosw
593 : 864864 : R(3,2)=B(2)*B(3)*u+B(1)*sinw
594 : :
595 : 864864 : R(1,3)=B(1)*B(3)*u+B(2)*sinw
596 : 864864 : R(2,3)=B(2)*B(3)*u-B(1)*sinw
597 : 864864 : R(3,3)=B(3)*B(3)*u+cosw
598 : : ! rotate
599 [ + + + + : 13837824 : S_out=matmul(R, S_in)
+ + ]
600 : 864864 : end function rotate_S_DM
601 : :
602 [ + - + - : 2002 : subroutine spin_mover_t_run_one_step_DM(self, effpot, S_in, etot, displacement, strain,&
- + - + -
- - - ]
603 [ - + - - : 2002 : & lwf, energy_table)
- - ]
604 : : ! Depondt & Mertens (2009) method, using a rotation matrix so length doesn't change.
605 : : !class (spin_mover_t), intent(inout):: self
606 : : class(spin_mover_t), intent(inout):: self
607 : : class(abstract_potential_t), intent(inout) :: effpot
608 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), lwf(:)
609 : : real(dp), intent(inout) :: S_in(3,self%nspin)
610 : : real(dp), intent(out) :: etot
611 : : type(hash_table_t),optional, intent(inout) :: energy_table
612 : : real(dp) :: Htmp(3)
613 : : integer :: i
614 : :
615 : : ! predict
616 [ + + + + ]: 1731730 : self%Stmp2(:,:)=0.0_dp
617 : 2002 : etot=0.0
618 [ + + + + ]: 1731730 : self%Heff_tmp(:,:)=0.0_dp
619 : : call effpot%calculate(displacement=displacement, strain=strain, lwf=lwf,spin=S_in, &
620 [ - + + - : 6006 : bfield=self%Heff_tmp, energy=etot, energy_table=energy_table)
+ - ]
621 : 2002 : call self%get_Langevin_Heff(self%H_lang)
622 [ + + ]: 434434 : do i=self%mps%istart, self%mps%iend
623 [ + + ]: 1729728 : Htmp=self%Heff_tmp(:,i)+self%H_lang(:,i)
624 : : ! Note that there is no - , because dsdt =-cross (S, Hrotate)
625 [ + + + + ]: 3027024 : self%Hrotate(:,i) = self%gamma_L(i) * (Htmp + self%damping(i)* cross(S_in(:,i), Htmp))
626 [ + + ]: 1731730 : self%Stmp2(:,i)= rotate_S_DM(S_in(:,i), self%Hrotate(:,i), self%dt)
627 : : end do
628 : 2002 : call self%mps%allgatherv_dp2d(self%Stmp2, 3, self%buffer)
629 : :
630 : : ! correction
631 [ + + + + ]: 1731730 : self%Htmp(:,:)=0.0_dp
632 : 2002 : etot=0.0
633 : : call effpot%calculate(displacement=displacement, strain=strain, lwf=lwf, spin=self%Stmp2, &
634 [ - + + - : 6006 : bfield=self%Htmp, energy=etot, energy_table=energy_table)
+ - ]
635 : :
636 [ + + ]: 434434 : do i=self%mps%istart, self%mps%iend
637 [ + + ]: 1729728 : Htmp=(self%Heff_tmp(:,i)+self%Htmp(:,i))*0.5_dp + self%H_lang(:,i)
638 [ + + + + ]: 3027024 : self%Hrotate(:,i) = self%gamma_L(i) * (Htmp + self%damping(i)* cross(S_in(:,i), Htmp))
639 [ + + ]: 1731730 : self%Stmp(:, i)= rotate_S_DM(S_in(:,i), self%Hrotate(:,i), self%dt)
640 : : end do
641 : 2002 : call self%mps%allgatherv_dp2d(self%Stmp, 3, self%buffer)
642 : 2002 : end subroutine spin_mover_t_run_one_step_DM
643 : :
644 [ # # # # : 0 : subroutine spin_mover_t_run_one_step_MC(self, effpot, S_in, etot, displacement, strain, lwf, energy_table)
# # ]
645 : : class(spin_mover_t), intent(inout) :: self
646 : : class(abstract_potential_t), intent(inout) :: effpot
647 : : real(dp), optional, intent(inout) :: displacement(:, :), strain(:,:), lwf(:)
648 : : real(dp), intent(inout) :: S_in(3,self%nspin)
649 : : real(dp), intent(out) :: etot
650 : : type(hash_table_t),optional, intent(inout) :: energy_table
651 [ # # # # ]: 0 : if(present(displacement) .or. present(lwf) .or. present(strain)) then
652 : 0 : ABI_BUG("Monte Carlo only implemented for spin.")
653 : : else
654 : 0 : call self%spin_mc%run_MC(self%rng, effpot, S_in, etot, bfield=self%Htmp)
655 : : end if
656 : 0 : call energy_table%put(self%label, etot)
657 : 0 : end subroutine spin_mover_t_run_one_step_MC
658 : :
659 [ - + + - : 2002 : subroutine spin_mover_t_run_one_step(self, effpot, displacement, strain, spin, lwf, energy_table)
+ - - + -
+ - - - -
- + - - -
- ]
660 : : class(spin_mover_t), intent(inout) :: self
661 : : class(abstract_potential_t), intent(inout) :: effpot
662 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
663 : : real(dp) :: etot
664 : : type(hash_table_t),optional, intent(inout) :: energy_table
665 : :
666 [ - + ]: 2002 : if(present(spin)) ABI_ERROR("spin should not be input for spin mover.")
667 [ - + ]: 2002 : if(self%method==1) then
668 : : call self%run_one_step_HeunP(effpot=effpot, S_in=self%Stmp, etot=etot, &
669 [ # # # # : 0 : displacement=displacement, strain=strain, lwf=lwf, energy_table=energy_table)
# # ]
670 [ + - ]: 2002 : else if (self%method==2) then
671 : : call self%run_one_step_DM(effpot=effpot, S_in=self%Stmp, etot=etot,&
672 [ - + + - : 6006 : displacement=displacement, strain=strain, lwf=lwf, energy_table=energy_table)
+ - ]
673 [ # # ]: 0 : else if (self%method==3) then
674 [ # # # # ]: 0 : if(present(displacement) .or. present(strain) .or. present(lwf)) then
675 : 0 : ABI_ERROR("Monte carlo not implemented for lattice and lwf yet.")
676 : : endif
677 : 0 : call self%run_one_step_MC(effpot, self%Stmp, etot, energy_table=energy_table)
678 [ # # ]: 0 : else if (self%method==20) then
679 : : call self%run_one_step_dummy(effpot=effpot, S_in=self%Stmp, etot=etot, &
680 [ # # # # : 0 : displacement=displacement, strain=strain, lwf=lwf, energy_table=energy_table)
# # ]
681 : : end if
682 : :
683 : : ! do not inc until time is set to hist.
684 : : ! run one step does not know about time. So it will be done in the outer loop.
685 [ + - ]: 2002 : if(self%mps%irank==0) then
686 : 2002 : call self%hist%set_vars(S=self%Stmp, Snorm=effpot%supercell%spin%ms, etot=etot, inc=.False.)
687 : : end if
688 : 2002 : end subroutine spin_mover_t_run_one_step
689 : :
690 : : !!****f* m_spin_mover/spin_mover_t_run_time
691 : : !!
692 : : !! NAME
693 : : !! spin_mover_t_run_time
694 : : !!
695 : : !! FUNCTION
696 : : !! run all spin step
697 : : !!
698 : : !! INPUTS
699 : : !!
700 : : !! OUTPUT
701 : : !!
702 : : !! NOTES
703 : : !!
704 : : !!
705 : : !! SOURCE
706 [ # # # # : 0 : subroutine spin_mover_t_run_time(self, calculator, displacement, strain, spin, lwf, energy_table)
# # # # #
# # # # #
# # # # #
# # # #
# ]
707 : :
708 : : class(spin_mover_t), intent(inout):: self
709 : : class(abstract_potential_t), intent(inout) :: calculator
710 : :
711 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), lwf(:), spin(:,:)
712 : :
713 : : type(hash_table_t),optional, intent(inout) :: energy_table
714 : : !type(spin_hist_t), intent(inout) :: hist
715 : : !type(spin_ncfile_t), intent(inout) :: ncfile
716 : : !type(spin_observable_t), intent(inout) :: ob
717 : : !real(dp) :: S(3, self%nspin)
718 : : real(dp):: t, etotal
719 : : integer :: counter, i, ii
720 : : character(len=80) :: msg, msg_empty
721 : :
722 : : integer :: master, my_rank, comm, nproc
723 : : logical :: iam_master
724 : :
725 : 0 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
726 : :
727 : 0 : t=0.0
728 : 0 : counter=0
729 [ # # ]: 0 : if(iam_master) then
730 : 0 : msg_empty=ch10
731 : :
732 : 0 : msg=repeat("=", 80)
733 : 0 : call wrtout(std_out,msg,'COLL')
734 : 0 : call wrtout(ab_out, msg, 'COLL')
735 : 0 : write(msg, '(A20)') "Spin dynamic steps:"
736 : 0 : call wrtout(std_out,msg,'COLL')
737 : 0 : call wrtout(ab_out, msg, 'COLL')
738 : 0 : msg=repeat("=", 80)
739 : 0 : call wrtout(std_out,msg,'COLL')
740 : 0 : call wrtout(ab_out, msg, 'COLL')
741 : :
742 : 0 : write(msg, "(A13, 4X, A13, 6X, A13, 4X, A13)") "Iteration", "time(s)", "Avg_Mst/Ms", "ETOT(Ha/uc)"
743 : 0 : call wrtout(std_out,msg,'COLL')
744 : 0 : call wrtout(ab_out, msg, 'COLL')
745 : :
746 : 0 : msg=repeat("-", 80)
747 : 0 : call wrtout(std_out,msg,'COLL')
748 : 0 : call wrtout(ab_out, msg, 'COLL')
749 : : end if
750 : :
751 [ # # ]: 0 : if (abs(self%thermal_time) > 1e-30) then
752 [ # # ]: 0 : if (iam_master) then
753 : 0 : msg="Thermalization run:"
754 : 0 : call wrtout(std_out,msg,'COLL')
755 : 0 : call wrtout(ab_out, msg, 'COLL')
756 : : end if
757 : :
758 [ # # ]: 0 : do while(t<self%thermal_time)
759 : 0 : counter=counter+1
760 : : call self%run_one_step(effpot=calculator, displacement=displacement, strain=strain, &
761 [ # # # # : 0 : & lwf=lwf, energy_table=energy_table)
# # ]
762 [ # # ]: 0 : if (iam_master) then
763 : 0 : call self%hist%set_vars( time=t, inc=.True.)
764 [ # # ]: 0 : if(mod(counter, self%hist%spin_nctime)==0) then
765 : : call self%spin_ob%get_observables( self%hist%S(:,:, self%hist%ihist_prev), &
766 : 0 : self%hist%Snorm(:,self%hist%ihist_prev),self%hist%etot(self%hist%ihist_prev))
767 : 0 : etotal = energy_table%sum_val()
768 : 0 : write(msg, "(A1, 1X, I13, 4X, ES13.5, 4X, ES13.5, 4X, ES13.5)") "-", counter, t*Time_Sec, &
769 : 0 : & self%spin_ob%Mst_norm_total/self%spin_ob%Snorm_total, &
770 : 0 : & etotal/self%spin_ob%nscell
771 : : ! total : 13+4+...= 64
772 : 0 : call wrtout(std_out,msg,'COLL')
773 : 0 : call wrtout(ab_out, msg, 'COLL')
774 : : endif
775 : : end if
776 : 0 : t=t+self%dt
777 : : end do
778 : :
779 : 0 : t=0.0
780 : 0 : counter=0
781 [ # # ]: 0 : if (iam_master) then
782 : 0 : call self%hist%reset(array_to_zero=.False.)
783 : : end if
784 : : endif
785 [ # # ]: 0 : if(iam_master) then
786 : 0 : call self%spin_ob%reset()
787 : : endif
788 : :
789 [ # # ]: 0 : if (iam_master) then
790 : 0 : msg="Measurement run:"
791 : 0 : call wrtout(std_out,msg,'COLL')
792 : 0 : call wrtout(ab_out, msg, 'COLL')
793 : : end if
794 : :
795 [ # # ]: 0 : do while(t<self%total_time)
796 : 0 : counter=counter+1
797 : : call self%run_one_step(effpot=calculator, displacement=displacement, strain=strain, &
798 [ # # # # : 0 : & spin=spin, lwf=lwf, energy_table=energy_table)
# # # # ]
799 [ # # ]: 0 : if (iam_master) then
800 : 0 : call self%hist%set_vars(time=t, inc=.True.)
801 : : call self%spin_ob%get_observables(self%hist%S(:,:, self%hist%ihist_prev), &
802 : 0 : self%hist%Snorm(:,self%hist%ihist_prev), self%hist%etot(self%hist%ihist_prev))
803 [ # # ]: 0 : if(modulo(counter, self%hist%spin_nctime)==0) then
804 : 0 : call self%spin_ncfile%write_one_step(self%hist)
805 : 0 : etotal = energy_table%sum_val()
806 : 0 : write(msg, "(A1, 1X, I13, 4X, ES13.5, 4X, ES13.5, 4X, ES13.5)") "-", counter, t*Time_Sec, &
807 : 0 : & self%spin_ob%Mst_norm_total/self%spin_ob%Snorm_total, &
808 : 0 : & etotal/self%spin_ob%nscell
809 : 0 : call wrtout(std_out,msg,'COLL')
810 : 0 : call wrtout(ab_out, msg, 'COLL')
811 : : endif
812 : : end if
813 : 0 : t=t+self%dt
814 : : enddo
815 : :
816 [ # # ]: 0 : if (iam_master) then
817 : 0 : msg=repeat("-", 80)
818 : 0 : call wrtout(std_out,msg,'COLL')
819 : 0 : call wrtout(ab_out, msg, 'COLL')
820 : :
821 : 0 : write(msg, "(A27)") "Summary of spin dynamics:"
822 : 0 : call wrtout(std_out,msg,'COLL')
823 : 0 : call wrtout(ab_out, msg, 'COLL')
824 : :
825 : 0 : write(msg, "(A65)") "At the end of the run, the average spin at each sublattice is"
826 : 0 : call wrtout(std_out,msg,'COLL')
827 : 0 : call wrtout(ab_out, msg, 'COLL')
828 : :
829 : 0 : write(msg, "(6X, A10, 5X, 3A10, A11)") 'Sublattice', '<M_i>(x)', '<M_i>(y)', '<M_i>(z)', '||<M_i>||'
830 : 0 : call wrtout(std_out,msg,'COLL')
831 : 0 : call wrtout(ab_out, msg, 'COLL')
832 : :
833 [ # # ]: 0 : do i =1, self%spin_ob%nsublatt
834 : 0 : write(msg, "(A1, 5X, 2X, I5.4, 8X, 4F10.5)") '-', i, &
835 [ # # ]: 0 : (self%spin_ob%Mst_sub(ii,i)/self%spin_ob%nspin_sub(i)/mu_B , ii=1, 3), &
836 [ # # # # ]: 0 : sqrt(sum((self%spin_ob%Mst_sub(:, i)/self%spin_ob%nspin_sub(i)/mu_B)**2))
837 : 0 : call wrtout(std_out,msg,'COLL')
838 : 0 : call wrtout(ab_out, msg, 'COLL')
839 : : end do
840 : :
841 : 0 : call wrtout(std_out,msg_empty,'COLL')
842 : 0 : call wrtout(ab_out, msg_empty, 'COLL')
843 : :
844 : : write(msg, "(A1, 1X, A11, 3X, A13, 3X, A13, 3X, A13, 3X, A13 )" ) &
845 : 0 : "#", "Temperature", "Cv", "chi", "BinderU4", "Mst"
846 : 0 : call wrtout(std_out, msg, "COLL")
847 : 0 : call wrtout(ab_out, msg, "COLL")
848 : : write(msg, "(2X, F11.5, 3X, ES13.5, 3X, ES13.5, 3X, E13.5, 3X, ES13.5, 3X )" ) &
849 : 0 : self%temperature*Ha_K , self%spin_ob%Cv, self%spin_ob%chi, &
850 : 0 : self%spin_ob%binderU4, self%spin_ob%Avg_Mst_norm_total/self%spin_ob%snorm_total
851 : 0 : call wrtout(std_out, msg, "COLL")
852 : 0 : call wrtout(ab_out, msg, "COLL")
853 : :
854 : 0 : msg=repeat("=", 80)
855 : 0 : call wrtout(std_out,msg,'COLL')
856 : 0 : call wrtout(ab_out, msg, 'COLL')
857 : : end if
858 : 0 : end subroutine spin_mover_t_run_time
859 : : !!***
860 : :
861 : :
862 : :
863 : :
864 : : !!****f* m_spin_mover/run_MvT
865 : : !!
866 : : !! NAME
867 : : !! run_MvT
868 : : !!
869 : : !! FUNCTION
870 : : !! run M vs Temperature
871 : : !!
872 : : !! INPUTS
873 : : !! pot: potential
874 : : !! T_start, Tend, T_nstep
875 : : !u
876 : : !! OUTPUT
877 : : !!
878 : : !! SOURCE
879 [ # # # # : 0 : subroutine run_MvT(self, pot, ncfile_prefix, displacement, strain, spin, lwf, energy_table)
# # # # #
# # # # #
# # # # #
# # # #
# ]
880 : : class(spin_mover_t), intent(inout) :: self
881 : : class(abstract_potential_t), intent(inout) :: pot
882 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), lwf(:), spin(:,:)
883 : : character(fnlen), intent(inout) :: ncfile_prefix
884 : : type(hash_table_t), optional, intent(inout) :: energy_table
885 : : real(dp) :: T_start, T_end
886 : : integer :: T_nstep
887 : : !type(spin_ncfile_t) :: spin_ncfile
888 : : character(len=4) :: post_fname
889 : : real(dp) :: T, T_step
890 : : integer :: i, ii, Tfile, iostat
891 : : character(len=90) :: msg
892 : : character(len=4200) :: Tmsg ! to write to var T file
893 : : character(len=150) :: iomsg
894 : : character(fnlen) :: Tfname ! file name for output various T calculation
895 : 0 : real(dp), allocatable :: Tlist(:), chi_list(:), Cv_list(:), binderU4_list(:)
896 : 0 : real(dp), allocatable :: Mst_sub_norm_list(:, :)
897 : 0 : real(dp), allocatable :: Mst_norm_total_list(:)
898 : :
899 : : integer :: master, my_rank, comm, nproc, ierr
900 : : logical :: iam_master
901 : 0 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
902 : :
903 [ # # ]: 0 : if (iam_master) then
904 : 0 : T_start=self%params%spin_temperature_start
905 : 0 : T_end=self%params%spin_temperature_end
906 : 0 : T_nstep=self%params%spin_temperature_nstep
907 : 0 : Tfile=get_unit()
908 [ # # ]: 0 : Tfname = trim(ncfile_prefix)//'.varT'
909 : 0 : iostat=open_file(file=Tfname, unit=Tfile, iomsg=iomsg )
910 [ # # ]: 0 : if (T_nstep<=1) then
911 : : T_step=0.0
912 : : else
913 : 0 : T_step=(T_end-T_start)/(T_nstep-1)
914 : : endif
915 : : write(msg, "(A52, ES13.5, A11, ES13.5, A1)") &
916 : 0 : & "Starting temperature dependent calculations. T from ", &
917 : 0 : & T_start*Ha_K, "K to ", T_end*Ha_K, " K."
918 : 0 : call wrtout(std_out, msg, "COLL")
919 : 0 : call wrtout(ab_out, msg, "COLL")
920 : :
921 [ # # # # : 0 : ABI_MALLOC(Tlist, (T_nstep))
# # # # #
# # # ]
922 [ # # # # : 0 : ABI_MALLOC(chi_list, (T_nstep))
# # # # ]
923 [ # # # # : 0 : ABI_MALLOC(Cv_list, (T_nstep))
# # # # ]
924 [ # # # # : 0 : ABI_MALLOC(binderU4_list, (T_nstep))
# # # # ]
925 [ # # # # : 0 : ABI_MALLOC(Mst_sub_norm_list, (self%spin_ob%nsublatt, T_nstep))
# # # # #
# # # # #
# # ]
926 [ # # # # : 0 : ABI_MALLOC( Mst_norm_total_list, (T_nstep))
# # # # ]
927 : : end if
928 : :
929 : 0 : call xmpi_bcast(T_nstep, 0, comm, ierr)
930 : :
931 : : ! write header of varT file
932 [ # # ]: 0 : if(iam_master) then
933 : : write(Tmsg, "(A1, 1X, A11, 3X, A13, 3X, A13, 3X, A13, 3X, A13, 3X, *(I13, 3X) )" ) &
934 [ # # # # ]: 0 : "#", "Temperature (K)", "Cv (1)", "chi (1)", "BinderU4 (1)", "Mst/Ms(1)", (ii, ii=1, self%spin_ob%nsublatt)
935 : 0 : call wrtout(Tfile, Tmsg, "COLL")
936 : 0 : flush(Tfile)
937 : : endif
938 : :
939 : :
940 : :
941 : :
942 [ # # ]: 0 : do i=1, T_nstep
943 [ # # ]: 0 : if(iam_master) then
944 : 0 : T=T_start+(i-1)*T_step
945 : 0 : msg=repeat("=", 79)
946 : 0 : call wrtout(std_out, msg, "COLL")
947 : 0 : call wrtout(ab_out, msg, "COLL")
948 : :
949 : 0 : write(msg, "(A13, 5X, ES13.5, A3)") "Temperature: ", T*Ha_K, " K."
950 : 0 : call wrtout(std_out, msg, "COLL")
951 : 0 : call wrtout(ab_out, msg, "COLL")
952 : :
953 : 0 : call self%hist%reset(array_to_zero=.False.)
954 : : ! set temperature
955 : : ! TODO make this into a subroutine set_params
956 : 0 : self%params%spin_temperature=T
957 : : endif
958 : 0 : call self%set_temperature(temperature=T)
959 : :
960 [ # # ]: 0 : if(iam_master) then
961 : : call self%hist%set_params(spin_nctime=self%params%spin_nctime, &
962 : 0 : & spin_temperature=T)
963 : 0 : call self%spin_ob%reset(self%params)
964 : : endif
965 : : ! uncomment if then to use spin initializer at every temperature. otherwise use last temperature
966 [ # # ]: 0 : if(i==1) then
967 : 0 : call self%set_initial_state(mode=self%params%spin_init_state)
968 : : else
969 [ # # ]: 0 : if(iam_master) then
970 : 0 : call self%hist%inc1()
971 : : endif
972 : : endif
973 : :
974 [ # # ]: 0 : if(iam_master) then
975 : 0 : write(post_fname, "(I4.4)") i
976 : : call self%prepare_ncfile( self%params, &
977 [ # # ]: 0 : & trim(ncfile_prefix)//'_T'//post_fname//'_spinhist.nc')
978 : 0 : call self%spin_ncfile%write_one_step(self%hist)
979 : : endif
980 : :
981 : : ! run in parallel
982 : : call self%run_time(pot, displacement=displacement, strain=strain, spin=spin, &
983 [ # # # # : 0 : & lwf=lwf, energy_table=energy_table)
# # # # ]
984 : :
985 [ # # ]: 0 : if(iam_master) then
986 : 0 : call self%spin_ncfile%close()
987 : : ! save observables
988 : 0 : Tlist(i)=T
989 : 0 : chi_list(i)=self%spin_ob%chi
990 : 0 : Cv_list(i)=self%spin_ob%Cv
991 : 0 : binderU4_list(i)=self%spin_ob%binderU4
992 : : !Mst_sub_list(:,:,i)=self%spin_ob%Mst_sub(:,:) ! not useful
993 [ # # ]: 0 : Mst_sub_norm_list(:,i)=self%spin_ob%Avg_Mst_sub_norm(:)
994 : 0 : Mst_norm_total_list(i)=self%spin_ob%Avg_Mst_norm_total
995 : :
996 : : ! write to varT file
997 : : write(Tmsg, "(2X, F11.5, 3X, ES13.5, 3X, ES13.5, 3X, E13.5, 3X, ES13.5, 3X, *(ES13.5, 3X) )" ) &
998 : 0 : Tlist(i)*Ha_K, Cv_list(i), chi_list(i), binderU4_list(i), Mst_norm_total_list(i)/self%spin_ob%snorm_total,&
999 [ # # # # ]: 0 : & (Mst_sub_norm_list(ii,i)/mu_B, ii=1, self%spin_ob%nsublatt)
1000 : 0 : call wrtout(Tfile, Tmsg, "COLL")
1001 : 0 : flush(Tfile)
1002 : :
1003 : : endif
1004 : :
1005 : : end do
1006 : :
1007 : :
1008 [ # # ]: 0 : if(iam_master) then
1009 : : ! write summary of MvT run
1010 : 0 : msg=repeat("=", 79)
1011 : 0 : call wrtout(std_out, msg, "COLL")
1012 : 0 : call wrtout(ab_out, msg, "COLL")
1013 : :
1014 : 0 : write(msg, *) "Summary of various T run: "
1015 : 0 : call wrtout(std_out, msg, "COLL")
1016 : 0 : call wrtout(ab_out, msg, "COLL")
1017 : :
1018 : : write(msg, "(A1, 1X, A11, 3X, A13, 3X, A13, 3X, A13, 3X, A13)" ) &
1019 : 0 : "#", "Temperature", "Cv", "chi", "BinderU4", "Mst"
1020 : 0 : call wrtout(std_out, msg, "COLL")
1021 : 0 : call wrtout(ab_out, msg, "COLL")
1022 : :
1023 [ # # ]: 0 : do i = 1, T_nstep
1024 : : write(msg, "(2X, F11.5, 3X, ES13.5, 3X, ES13.5, 3X, E13.5, 3X, ES13.5 )" ) &
1025 : 0 : Tlist(i)*Ha_K, Cv_list(i), chi_list(i), binderU4_list(i), Mst_norm_total_list(i)/self%spin_ob%snorm_total
1026 : 0 : call wrtout(std_out, msg, "COLL")
1027 : 0 : call wrtout(ab_out, msg, "COLL")
1028 : : end do
1029 : :
1030 : 0 : msg=repeat("=", 79)
1031 : 0 : call wrtout(std_out, msg, "COLL")
1032 : 0 : call wrtout(ab_out, msg, "COLL")
1033 : :
1034 : :
1035 : : ! close varT file
1036 : 0 : iostat= close_unit(unit=Tfile, iomsg=iomsg)
1037 : :
1038 [ # # ]: 0 : ABI_FREE(Tlist)
1039 [ # # ]: 0 : ABI_FREE(chi_list)
1040 [ # # ]: 0 : ABI_FREE(Cv_list)
1041 [ # # ]: 0 : ABI_FREE(binderU4_list)
1042 [ # # ]: 0 : ABI_FREE(Mst_sub_norm_list)
1043 [ # # ]: 0 : ABI_FREE( Mst_norm_total_list)
1044 : :
1045 : : endif
1046 [ # # # # : 0 : end subroutine run_MvT
# # # # #
# # # ]
1047 : : !!***
1048 : :
1049 : : !!****f* m_spin_mover/current_spin
1050 : : !!
1051 : : !! NAME
1052 : : !! current_spin
1053 : : !!
1054 : : !! FUNCTION
1055 : : !! return the current spin state
1056 : : !!
1057 : : !! INPUTS
1058 : : !!
1059 : : !!
1060 : : !! OUTPUT
1061 : : !!
1062 : : !!
1063 : : !!
1064 : : !! SOURCE
1065 : 0 : function current_spin(self) result(ret)
1066 : : class(spin_mover_t), target, intent(inout) :: self
1067 : : real(dp), pointer :: ret(:,:)
1068 : : integer :: i
1069 : 0 : i=self%hist%findIndex(step=0)
1070 : 0 : ret => self%hist%S(:,:,i)
1071 : 0 : end function current_spin
1072 : : !!***
1073 : :
1074 : :
1075 : :
1076 : :
1077 : : !!****f* m_spin_mover/finalize
1078 : : !!
1079 : : !! NAME
1080 : : !! finalize
1081 : : !!
1082 : : !! FUNCTION
1083 : : !! finalize spin mover.
1084 : : !!
1085 : : !! INPUTS
1086 : : !!
1087 : : !! OUTPUT
1088 : : !!
1089 : : !! NOTES
1090 : : !! does nothing. But it's better to preserve initialize-finalize symmetry.
1091 : : !!
1092 : : !! SOURCE
1093 : 2 : subroutine finalize(self)
1094 : :
1095 : : class(spin_mover_t), intent(inout):: self
1096 [ + - ]: 2 : if(allocated(self%gyro_ratio) ) then
1097 : 2 : ABI_FREE(self%gyro_ratio)
1098 : : end if
1099 : :
1100 [ + - ]: 2 : if(allocated(self%damping) ) then
1101 : 2 : ABI_FREE(self%damping)
1102 : : end if
1103 : :
1104 [ + - ]: 2 : if(allocated(self%gamma_l) ) then
1105 : 2 : ABI_FREE(self%gamma_l)
1106 : : end if
1107 : :
1108 [ + - ]: 2 : if(allocated(self%H_lang_coeff) ) then
1109 : 2 : ABI_FREE(self%H_lang_coeff)
1110 : : end if
1111 : :
1112 [ + - ]: 2 : if(allocated(self%ms) ) then
1113 : 2 : ABI_FREE(self%ms)
1114 : : end if
1115 : :
1116 : :
1117 : :
1118 [ - + ]: 2 : if(self%method==3) then
1119 : 0 : call self%spin_mc%finalize()
1120 : : end if
1121 : :
1122 [ + - ]: 2 : if(allocated(self%Stmp)) then
1123 : 2 : ABI_FREE(self%Stmp)
1124 : : end if
1125 : :
1126 [ + - ]: 2 : if(allocated(self%Stmp2)) then
1127 : 2 : ABI_FREE(self%Stmp2)
1128 : : end if
1129 : :
1130 : :
1131 [ + - ]: 2 : if(allocated(self%Heff_tmp)) then
1132 : 2 : ABI_FREE(self%Heff_tmp)
1133 : : end if
1134 : :
1135 [ + - ]: 2 : if(allocated(self%Htmp)) then
1136 : 2 : ABI_FREE(self%Htmp)
1137 : : end if
1138 : :
1139 [ + - ]: 2 : if(allocated(self%Hrotate)) then
1140 : 2 : ABI_FREE(self%Hrotate)
1141 : : end if
1142 : :
1143 [ + - ]: 2 : if(allocated(self%H_lang)) then
1144 : 2 : ABI_FREE(self%H_lang)
1145 : : end if
1146 : :
1147 [ + - ]: 2 : if(allocated(self%buffer)) then
1148 : 2 : ABI_FREE(self%buffer)
1149 : : end if
1150 : :
1151 : :
1152 : 2 : nullify(self%supercell)
1153 : 2 : nullify(self%params)
1154 : 2 : nullify(self%rng)
1155 : 2 : call self%mps%finalize()
1156 : 2 : call self%hist%finalize()
1157 : 2 : call self%spin_ob%finalize()
1158 : : !call self%spin_ncfile%close()
1159 : 2 : call self%spin_ob%finalize()
1160 : 2 : end subroutine finalize
1161 : : !!***
1162 : :
1163 [ # # # # : 0 : end module m_spin_mover
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # # #
# # # ]
|
