Line data Source code
1 : !!****m* ABINIT/m_lattice_mover
2 : !! NAME
3 : !! m_lattice_mover
4 : !!
5 : !! FUNCTION
6 : !! This module contains the lattice mover.
7 : !!
8 : !!
9 : !! Datatypes:
10 : !!
11 : !! * lattice_mover_t: defines the lattice movers
12 : !!
13 : !! Subroutines:
14 : !! TODO: add this when F2003 doc style is determined.
15 : !!
16 : !!
17 : !! COPYRIGHT
18 : !! Copyright (C) 2001-2026 ABINIT group (hexu)
19 : !! This file is distributed under the terms of the
20 : !! GNU General Public License, see ~abinit/COPYING
21 : !! or http://www.gnu.org/copyleft/gpl.txt .
22 : !! For the initials of contributorsi see ~abinit/doc/developers/contributors.txt .
23 : !!
24 : !! SOURCE
25 :
26 :
27 :
28 : #if defined HAVE_CONFIG_H
29 : #include "config.h"
30 : #endif
31 :
32 : #include "abi_common.h"
33 :
34 : module m_lattice_mover
35 : use defs_basis
36 : use m_abicore
37 : use m_errors
38 : use m_xmpi
39 : use m_multibinit_dataset, only: multibinit_dtset_type
40 : use m_abstract_potential, only: abstract_potential_t
41 : use m_abstract_mover, only: abstract_mover_t
42 : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
43 : use m_random_xoroshiro128plus, only: rng_t
44 : use m_hashtable_strval, only: hash_table_t
45 : use m_lattice_ncfile, only: lattice_ncfile_t
46 : use m_mpi_scheduler, only: init_mpi_info
47 : !!***
48 :
49 : implicit none
50 : private
51 :
52 : type ,public, extends(abstract_mover_t) :: lattice_mover_t
53 : !> This is the abstract lattice mover
54 :
55 : type(multibinit_dtset_type), pointer :: params=>null() ! input parameters
56 : integer :: natom=0 ! number of atoms
57 : real(dp) :: stress(3,3), strain(3,3) ! stress and strain
58 : real(dp), allocatable :: masses(:) ! masses
59 : integer :: latt_dynamics=0 ! type of lattice dynamics
60 :
61 : ! hexu: is xcart needed?
62 : real(dp), allocatable :: current_xcart(:,:) ! xcart of current step
63 : real(dp), allocatable :: current_vcart(:,:) ! vcart of current step
64 : real(dp), allocatable :: forces(:,:) ! forces
65 : real(dp), allocatable :: displacement(:,:) ! displacement
66 : real(dp) :: energy ! total energy
67 : real(dp) :: Ek ! kinetic energy
68 : real(dp) :: T_ob ! observed temperature
69 : logical :: is_null = .True.
70 : !> TODO: hist
71 : !type(lattice_hist_t) :: hist
72 :
73 : real(dp) :: mass_total
74 : type(lattice_ncfile_t) :: ncfile
75 : contains
76 : procedure:: initialize ! perhaps each effpot type should have own
77 : procedure :: finalize
78 : procedure :: set_params
79 : procedure :: prepare_ncfile
80 : procedure :: set_ncfile_name
81 : procedure :: set_initial_state ! initial state
82 : procedure :: set_temperature
83 : procedure :: force_stationary
84 : procedure :: run_one_step
85 : procedure :: run_time
86 : procedure :: run_varT
87 : procedure :: reset ! reset the mover
88 : procedure :: get_T_and_Ek ! calculate temperature and kinetic energy.
89 : procedure :: calc_observables ! call functions to calculate observables
90 : procedure :: write_hist ! write hist file
91 : end type lattice_mover_t
92 :
93 :
94 : contains
95 :
96 : !-------------------------------------------------------------------!
97 : ! Initialize:
98 : !
99 : ! Inputs:
100 : !> params: input parameters
101 : !> supercell: supercell.
102 : !> rng: random number generator
103 : !-------------------------------------------------------------------!
104 5 : subroutine initialize(self, params, supercell, rng)
105 : class(lattice_mover_t), intent(inout) :: self
106 : type(multibinit_dtset_type),target, intent(in) :: params
107 : type(mbsupercell_t),target, intent(in) :: supercell
108 : type(rng_t), target, intent(in) :: rng
109 5 : self%params=>params
110 5 : self%supercell=>supercell
111 5 : self%label="Lattice Mover"
112 5 : self%natom = supercell%lattice%natom
113 15 : ABI_MALLOC(self%masses, (self%natom))
114 15 : ABI_MALLOC(self%displacement, (3, self%natom))
115 15 : ABI_MALLOC(self%current_xcart, (3, self%natom))
116 15 : ABI_MALLOC(self%current_vcart, (3, self%natom))
117 15 : ABI_MALLOC(self%forces, (3,self%natom))
118 5 : self%is_null=.False.
119 65 : self%strain(:,:) = 0.0
120 65 : self%stress(:,:) = 0.0
121 18565 : self%forces(:,:) = 0.0
122 18565 : self%displacement(:,:) = 0.0
123 18565 : self%current_vcart(:,:) = 0.0
124 5 : call self%set_params(params)
125 5 : call self%set_rng(rng)
126 5 : end subroutine initialize
127 :
128 : !-------------------------------------------------------------------!
129 : ! Finalize:
130 : !-------------------------------------------------------------------!
131 5 : subroutine finalize(self)
132 : class(lattice_mover_t), intent(inout) :: self
133 5 : nullify(self%supercell)
134 5 : nullify(self%params)
135 5 : self%label="Destroyed lattice mover"
136 5 : if (.not.self%is_null) then
137 5 : ABI_FREE(self%masses)
138 5 : ABI_FREE(self%current_xcart)
139 5 : ABI_FREE(self%current_vcart)
140 5 : ABI_FREE(self%forces)
141 5 : ABI_FREE(self%displacement)
142 : endif
143 5 : self%is_null=.True.
144 5 : end subroutine finalize
145 :
146 : !-------------------------------------------------------------------!
147 : ! Set the mover using the input parameters
148 : !-------------------------------------------------------------------!
149 5 : subroutine set_params(self, params)
150 : ! set parameters from input file. (something else, like temperature for MvT calculation?)
151 : class(lattice_mover_t), intent(inout) :: self
152 : type(multibinit_dtset_type) :: params
153 5 : self%temperature = params%temperature !TODO: to Hartree ??
154 5 : self%dt =params%dtion
155 4645 : self%masses(:)=self%supercell%lattice%masses(:)
156 4645 : self%mass_total = sum(self%masses)
157 5 : self%total_time = self%dt * params%ntime
158 5 : self%latt_dynamics = params%dynamics
159 5 : end subroutine set_params
160 :
161 : !-------------------------------------------------------------------!
162 : ! Set the mover temperature
163 : !-------------------------------------------------------------------!
164 0 : subroutine set_temperature(self, temperature)
165 : class(lattice_mover_t), intent(inout) :: self
166 : real(dp),intent(in) :: temperature
167 0 : self%temperature = temperature !TODO: to Hartree ??
168 0 : end subroutine set_temperature
169 :
170 :
171 : !-------------------------------------------------------------------!
172 : ! set initial state:
173 : ! Inputs:
174 : ! mode: integer
175 : ! if mode=1, use a Boltzman distribution to init the velocities.
176 : ! if mode=2, ...
177 : !-------------------------------------------------------------------!
178 5 : subroutine set_initial_state(self, mode)
179 : ! set initial positions, spin, etc
180 : class(lattice_mover_t), intent(inout) :: self
181 : integer, optional, intent(in) :: mode
182 10 : real(dp) :: xi(3, self%natom)
183 : integer :: i
184 :
185 :
186 5 : if(mode==1) then ! using a boltzmann distribution.
187 : ! Should only be used for a constant Temperature mover
188 : ! which includes:
189 : ! 102: Langevin
190 : ! 103: Brendesen
191 5 : if (.not.( &
192 : self%latt_dynamics==101 .or. & ! TODO remove
193 : self%latt_dynamics==102 .or. self%latt_dynamics==103 ) ) then
194 0 : ABI_ERROR("Only set lattice initial state with a Boltzmann distribution in a constant T mover.")
195 : end if
196 5 : call self%rng%rand_normal_array(xi, 3*self%natom)
197 4645 : do i=1, self%natom
198 18565 : self%current_vcart(:,i) = xi(:, i) *sqrt(self%temperature/self%masses(i))
199 : end do
200 5 : call self%force_stationary()
201 5 : call self%get_T_and_Ek()
202 18565 : self%current_xcart(:, :) = self%supercell%lattice%xcart(:,:)
203 0 : else if(mode==2) then ! Use reference structure and 0 velocity.
204 : ! other modes.
205 0 : if(self%latt_dynamics==102 .or. self%latt_dynamics==103 ) then
206 0 : ABI_ERROR("Displacement and velocity set to zero in a NVT mover.")
207 : end if
208 0 : do i=1, self%natom
209 0 : self%current_vcart(:,i) = 0.0
210 : end do
211 0 : self%current_xcart(:, :) = self%supercell%lattice%xcart(:,:)
212 0 : call self%get_T_and_Ek()
213 : end if
214 :
215 :
216 5 : end subroutine set_initial_state
217 :
218 3 : subroutine prepare_ncfile(self, params, fname)
219 : class(lattice_mover_t), intent(inout) :: self
220 : type(multibinit_dtset_type) :: params
221 : character(len=*), intent(in) :: fname
222 : integer :: master, my_rank, comm, nproc
223 : logical :: iam_master
224 3 : ABI_UNUSED_A(params)
225 3 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
226 3 : if(iam_master) then
227 3 : call self%ncfile%initialize( trim(fname), 1)
228 3 : call self%ncfile%write_cell(self%supercell)
229 3 : call self%ncfile%def_lattice_var()
230 : end if
231 3 : end subroutine prepare_ncfile
232 :
233 :
234 : !-------------------------------------------------------------------!
235 : !set_ncfile_name :
236 : !-------------------------------------------------------------------!
237 3 : subroutine set_ncfile_name(self, params, fname)
238 : class(lattice_mover_t), intent(inout) :: self
239 : type(multibinit_dtset_type) :: params
240 : character(len=fnlen), intent(in) :: fname
241 : integer :: master, my_rank, comm, nproc
242 : logical :: iam_master
243 3 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
244 3 : if (iam_master) then
245 3 : call self%prepare_ncfile(params, trim(fname)//'_latthist.nc')
246 3 : call self%ncfile%write_one_step(self%current_xcart, self%current_vcart, self%energy, self%Ek )
247 : endif
248 3 : end subroutine set_ncfile_name
249 :
250 :
251 :
252 : !-------------------------------------------------------------------!
253 : ! Make sure the mass center does not move.
254 : !-------------------------------------------------------------------!
255 2005 : subroutine force_stationary(self)
256 : class(lattice_mover_t), intent(inout) :: self
257 : integer :: i
258 : real(dp) :: p(3), pavg(3)
259 2005 : p(:)=0.0
260 2166645 : do i = 1, self%natom
261 8660565 : p(:)=p(:)+self%current_vcart(:,i) * self%masses(i)
262 : end do
263 8020 : pavg=p/self%mass_total
264 2166645 : do i = 1, self%natom
265 8660565 : self%current_vcart(:, i) = self%current_vcart(:, i) - pavg(:)
266 : end do
267 2005 : end subroutine force_stationary
268 :
269 :
270 : !-------------------------------------------------------------------!
271 : ! Force the temperature strictly.
272 : ! Since the boltzman distribution has some fluctuation
273 : !-------------------------------------------------------------------!
274 : subroutine force_temperature(self)
275 : class(lattice_mover_t), intent(inout) :: self
276 : ABI_UNUSED_A(self)
277 : end subroutine force_temperature
278 :
279 :
280 : !-------------------------------------------------------------------!
281 : ! run_one_step
282 : ! run one step of dynamics.
283 : ! Should be overrided.
284 : ! Inputs:
285 : !> effpot: effective potential
286 : !> displacement: should NOT be provided, since it is already stored.
287 : !> strain: Also should NOT be provided.
288 : !> spin: should be provided only if there is spin-lattice coupling
289 : !> lwf : should be provided only if there is lattice-lwf coupling (unlikely)
290 : !> energy_table: energy_table.
291 : !-------------------------------------------------------------------!
292 0 : subroutine run_one_step(self, effpot, displacement, strain, spin, lwf, energy_table)
293 : ! run one step. (For MC also?)
294 : class(lattice_mover_t), intent(inout) :: self ! array of effective potentials so that there can be multiple of them.
295 : class(abstract_potential_t), intent(inout) :: effpot
296 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
297 : type(hash_table_t), optional, intent(inout) :: energy_table
298 :
299 : character(len=40) :: key
300 :
301 0 : if(present(displacement) .or. present(strain)) then
302 0 : ABI_ERROR("displacement and strain should not be input for lattice mover")
303 : end if
304 :
305 0 : ABI_BUG("The abstract lattice mover is used, which should be a bug.")
306 :
307 0 : ABI_UNUSED_A(self)
308 0 : ABI_UNUSED_A(effpot)
309 0 : ABI_UNUSED_A(displacement)
310 0 : ABI_UNUSED_A(strain)
311 0 : ABI_UNUSED_A(spin)
312 0 : ABI_UNUSED_A(lwf)
313 0 : ABI_UNUSED_A(energy_table)
314 :
315 0 : call self%get_T_and_Ek()
316 0 : if (present(energy_table)) then
317 0 : key = 'Lattice kinetic energy'
318 0 : call energy_table%put(key, self%Ek)
319 : end if
320 0 : end subroutine run_one_step
321 :
322 :
323 : !-------------------------------------------------------------------!
324 : !get_temperature_and_kinetic_energy
325 : ! Ek = 1/2 \sum m_i vi^2
326 : ! T = 2/3 Ek/natom (in a.u.)
327 : !-------------------------------------------------------------------!
328 4017 : subroutine get_T_and_Ek(self)
329 : class(lattice_mover_t), intent(inout) :: self
330 : integer :: i
331 4017 : self%Ek=0.0
332 4334017 : do i =1, self%natom
333 : self%Ek = self%Ek+ 0.5* self%masses(i) * &
334 17324017 : &sum(self%current_vcart(:,i)*self%current_vcart(:,i))
335 : end do
336 : ! temperature
337 4017 : self%T_ob = 2.0*self%Ek/(3*self%natom)
338 4017 : end subroutine get_T_and_Ek
339 :
340 :
341 : !-------------------------------------------------------------------!
342 : ! run from begining to end.
343 : !-------------------------------------------------------------------!
344 3 : subroutine run_time(self, effpot, displacement, strain, spin, lwf, energy_table)
345 : class(lattice_mover_t), intent(inout) :: self
346 : ! array of effective potentials so that there can be multiple of them.
347 : class(abstract_potential_t), intent(inout) :: effpot
348 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
349 : type(hash_table_t), optional, intent(inout) :: energy_table
350 : integer :: i, nstep
351 : character(len=90) :: msg
352 3 : if(present(displacement) .or. present(strain)) then
353 0 : ABI_ERROR("displacement and strain should not be input for lattice mover")
354 : end if
355 3 : ABI_UNUSED_A(self)
356 3 : ABI_UNUSED_A(effpot)
357 3 : ABI_UNUSED_A(spin)
358 3 : ABI_UNUSED_A(lwf)
359 3 : ABI_UNUSED_A(energy_table)
360 :
361 :
362 3 : msg=repeat("=", 90)
363 3 : call wrtout(std_out,msg,'COLL')
364 3 : call wrtout(ab_out, msg, 'COLL')
365 3 : write(msg, '(A22)') "Lattice dynamic steps:"
366 3 : call wrtout(std_out,msg,'COLL')
367 3 : call wrtout(ab_out, msg, 'COLL')
368 3 : msg=repeat("=", 90)
369 3 : call wrtout(std_out,msg,'COLL')
370 3 : call wrtout(ab_out, msg, 'COLL')
371 :
372 : write(msg, "(A13, 4X, A15, 4X, A15, 4X, A15, 4X, A15)") &
373 3 : & "Iteration", "temperature(K)", "Ekin(Ha/uc)", &
374 6 : & "Epot(Ha/uc)", "ETOT(Ha/uc)"
375 3 : call wrtout(std_out,msg,'COLL')
376 3 : call wrtout(ab_out, msg, 'COLL')
377 :
378 3 : nstep=floor(self%thermal_time/self%dt)
379 3 : do i =1, nstep
380 3 : call self%run_one_step(effpot=effpot, spin=spin, lwf=lwf, energy_table=energy_table)
381 : end do
382 :
383 3 : nstep=floor(self%total_time/self%dt)
384 2013 : do i =1, nstep
385 6030 : call self%run_one_step(effpot=effpot, spin=spin, lwf=lwf, energy_table=energy_table)
386 2013 : if(modulo(i, self%params%nctime)==0) then
387 210 : write(msg, "(I13, 4X, F15.5, 4X, ES15.5, 4X, ES15.5, 4X, ES15.5)") i, self%T_ob*Ha_K, &
388 210 : & self%Ek/self%supercell%ncell, self%energy/self%supercell%ncell, &
389 420 : & (self%Ek+self%energy)/self%supercell%ncell
390 210 : call wrtout(std_out,msg,'COLL')
391 210 : call wrtout(ab_out, msg, 'COLL')
392 877220 : self%current_xcart = self%supercell%lattice%xcart+self%displacement
393 210 : call self%ncfile%write_one_step(self%current_xcart, self%current_vcart, self%energy, self%Ek)
394 : end if
395 : !TODO: output, observables
396 : end do
397 :
398 3 : msg=repeat("=", 90)
399 3 : call wrtout(std_out,msg,'COLL')
400 3 : call wrtout(ab_out, msg, 'COLL')
401 :
402 3 : end subroutine run_time
403 :
404 :
405 : !-------------------------------------------------------------------!
406 : ! Reset:
407 : ! It reset the counter of steps, but does not set the initial state
408 : ! again.
409 : !-------------------------------------------------------------------!
410 0 : subroutine reset(self)
411 : ! reset the state of mover (e.g. counter->0)
412 : ! so it can be reused.
413 : class(lattice_mover_t), intent(inout) :: self
414 0 : ABI_UNUSED_A(self)
415 0 : end subroutine reset
416 :
417 : !-------------------------------------------------------------------!
418 : !Calc_observables
419 : !
420 : !-------------------------------------------------------------------!
421 0 : subroutine calc_observables(self)
422 : ! call functions to calculate observables.
423 : class(lattice_mover_t), intent(inout) :: self
424 0 : ABI_UNUSED_A(self)
425 0 : end subroutine calc_observables
426 :
427 : !-------------------------------------------------------------------!
428 : ! Write_hist: write to hist file
429 : !-------------------------------------------------------------------!
430 0 : subroutine write_hist(self)
431 : ! write to hist file
432 : class(lattice_mover_t), intent(inout) :: self
433 0 : ABI_UNUSED_A(self)
434 :
435 0 : end subroutine write_hist
436 :
437 : !-------------------------------------------------------------------!
438 : ! Get_state: get the current state
439 : !-------------------------------------------------------------------!
440 : subroutine get_state(self, displacement, strain, spin, lwf, ihist)
441 : ! get the state of the ihist(th) step. ihist can be 0 (current), -1 (last), ... -maxhist..
442 : class(lattice_mover_t), intent(in):: self
443 : real(dp), optional, intent(inout) :: displacement, strain, spin, lwf
444 : integer, optional, intent(in):: ihist
445 : ABI_UNUSED_A(self)
446 : ABI_UNUSED_A(displacement)
447 : ABI_UNUSED_A(strain)
448 : ABI_UNUSED_A(spin)
449 : ABI_UNUSED_A(lwf)
450 : ABI_UNUSED_A(ihist)
451 : end subroutine get_state
452 :
453 :
454 : !!****f* m_lwf_mover/run_varT
455 : !!
456 : !! NAME
457 : !! run_varT
458 : !!
459 : !! FUNCTION
460 : !! run M vs Temperature
461 : !!
462 : !! INPUTS
463 : !! pot: potential
464 : !! T_start, Tend, T_nstep
465 : !u
466 : !! OUTPUT
467 : !!
468 : !! SOURCE
469 0 : subroutine run_varT(self, pot, ncfile_prefix, displacement, strain, spin, lwf, energy_table)
470 : class(lattice_mover_t), intent(inout) :: self
471 : class(abstract_potential_t), intent(inout) :: pot
472 : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), lwf(:), spin(:,:)
473 : character(fnlen), intent(inout) :: ncfile_prefix
474 : type(hash_table_t), optional, intent(inout) :: energy_table
475 : real(dp) :: T_start, T_end
476 : integer :: T_nstep
477 : !type(lwf_ncfile_t) :: lwf_ncfile
478 : character(len=4) :: post_fname
479 : real(dp) :: T, T_step
480 : integer :: i
481 : !integer :: Tfile, iostat
482 : character(len=90) :: msg
483 : !character(len=4200) :: Tmsg ! to write to var T file
484 : !character(len=150) :: iomsg
485 : !character(fnlen) :: Tfname ! file name for output various T calculation
486 : !real(dp), allocatable :: Tlist(:), chi_list(:), Cv_list(:), binderU4_list(:)
487 : !real(dp), allocatable :: Mst_sub_norm_list(:, :)
488 : !real(dp), allocatable :: Mst_norm_total_list(:)
489 :
490 : integer :: master, my_rank, comm, nproc, ierr
491 : logical :: iam_master
492 0 : call init_mpi_info(master, iam_master, my_rank, comm, nproc)
493 :
494 0 : ABI_UNUSED_A(displacement)
495 0 : ABI_UNUSED_A(strain)
496 :
497 0 : if (iam_master) then
498 0 : T_start=self%params%latt_temperature_start
499 0 : T_end=self%params%latt_temperature_end
500 0 : T_nstep=self%params%latt_temperature_nstep
501 : !Tfile=get_unit()
502 : !Tfname = trim(ncfile_prefix)//'.varT'
503 : !iostat=open_file(file=Tfname, unit=Tfile, iomsg=iomsg )
504 0 : if (T_nstep<=1) then
505 : T_step=0.0
506 : else
507 0 : T_step=(T_end-T_start)/(T_nstep-1)
508 : endif
509 : write(msg, "(A52, ES13.5, A11, ES13.5, A1)") &
510 0 : & "Starting temperature dependent calculations. T from ", &
511 0 : & T_start*Ha_K, "K to ", T_end*Ha_K, " K."
512 0 : call wrtout(std_out, msg, "COLL")
513 0 : call wrtout(ab_out, msg, "COLL")
514 : end if
515 :
516 0 : call xmpi_bcast(T_nstep, 0, comm, ierr)
517 0 : do i=1, T_nstep
518 0 : if(iam_master) then
519 0 : T=T_start+(i-1)*T_step
520 0 : msg=repeat("=", 79)
521 0 : call wrtout(std_out, msg, "COLL")
522 0 : call wrtout(ab_out, msg, "COLL")
523 :
524 0 : write(msg, "(A13, 5X, ES13.5, A3)") "Temperature: ", T*Ha_K, " K."
525 0 : call wrtout(std_out, msg, "COLL")
526 0 : call wrtout(ab_out, msg, "COLL")
527 :
528 : ! set temperature
529 : ! TODO make this into a subroutine set_params
530 : endif
531 0 : call self%set_temperature(temperature=T)
532 0 : if(iam_master) then
533 0 : if(i==1) then
534 0 : call self%set_initial_state(mode=1)
535 : endif
536 :
537 0 : write(post_fname, "(I4.4)") i
538 : call self%prepare_ncfile( self%params, &
539 0 : & trim(ncfile_prefix)//'_T'//post_fname//'_latthist.nc')
540 0 : call self%ncfile%write_one_step(self%current_xcart, self%current_vcart, self%energy, self%Ek )
541 : endif
542 :
543 : call self%run_time(pot, spin=spin, &
544 0 : & lwf=lwf, energy_table=energy_table)
545 :
546 0 : if(iam_master) then
547 0 : call self%ncfile%finalize()
548 : endif
549 : end do
550 :
551 0 : end subroutine run_varT
552 : !!***
553 :
554 :
555 :
556 15 : end module m_lattice_mover
557 :
|
