Branch data Line data Source code
1 : : !!****m* ABINIT/m_lattice_langevin_mover
2 : : !! NAME
3 : : !! m_lattice_langevin_mover
4 : : !!
5 : : !! FUNCTION
6 : : !! This module contains the langevin (NVT) lattice mover.
7 : : !! It is a translation from the ASE (GPL licenced) Langevin mover python code to fortran.
8 : : !! The original code can be found at
9 : : !! https://gitlab.com/ase/ase/blob/master/ase/md/langevin.py
10 : : !! The method is described in
11 : : !! E. V.-Eijnden, and G. Ciccotti, Chem. Phys. Lett. 429, 310 (2006)
12 : : !! https://doi.org/10.1016/j.cplett.2006.07.086
13 : : !!
14 : : !!
15 : : !! Datatypes:
16 : : !!
17 : : !! * lattice_langevin_mover_t: defines the lattice movers
18 : : !!
19 : : !! Subroutines:
20 : : !! TODO: add this when F2003 doc style is determined.
21 : : !!
22 : : !!
23 : : !! COPYRIGHT
24 : : !! Copyright (C) 2001-2026 ABINIT group (hexu)
25 : : !! This file is distributed under the terms of the
26 : : !! GNU General Public License, see ~abinit/COPYING
27 : : !! or http://www.gnu.org/copyleft/gpl.txt .
28 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
29 : : !!
30 : : !! SOURCE
31 : :
32 : :
33 : :
34 : : #if defined HAVE_CONFIG_H
35 : : #include "config.h"
36 : : #endif
37 : :
38 : : #include "abi_common.h"
39 : :
40 : : module m_lattice_langevin_mover
41 : : use defs_basis
42 : : use m_abicore
43 : : use m_errors
44 : :
45 : : use m_multibinit_dataset, only: multibinit_dtset_type
46 : : use m_abstract_potential, only: abstract_potential_t
47 : : use m_abstract_mover, only: abstract_mover_t
48 : : use m_lattice_mover, only: lattice_mover_t
49 : : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
50 : : use m_random_xoroshiro128plus, only: rng_t
51 : : use m_hashtable_strval, only: hash_table_t
52 : : !!***
53 : :
54 : : implicit none
55 : :
56 : : private
57 : :
58 : : !-------------------------------------------------------------------!
59 : : ! Lattice_langevin_mover_t
60 : : !
61 : : !-------------------------------------------------------------------!
62 : : type, public, extends(lattice_mover_t) :: lattice_langevin_mover_t
63 : : ! c1 to c5 are a constants (temperature and mass dependent)
64 : : ! c3 c4 c5 has dimension of natom.
65 : : real(dp) :: c1, c2
66 : : real(dp), allocatable :: c3(:), c4(:), c5(:)
67 : : real(dp) :: fr =1e-4 ! friction, usually 1e-2~1e-4
68 : : ! xi and eta: random numbers of dimension (3, natom)
69 : : real(dp), allocatable :: xi(:,:), eta(:,:)
70 : : contains
71 : : procedure :: initialize
72 : : procedure :: finalize
73 : : procedure :: run_one_step
74 : : procedure :: update_vars
75 : : end type lattice_langevin_mover_t
76 : :
77 : : contains
78 : :
79 : : !-------------------------------------------------------------------!
80 : : ! Initialize:
81 : : ! read parameters
82 : : ! allocate memory and call update_vars()
83 : : !-------------------------------------------------------------------!
84 : 1 : subroutine initialize(self,params, supercell, rng)
85 : : class(lattice_langevin_mover_t), intent(inout) :: self
86 : : type(multibinit_dtset_type), target, intent(in):: params
87 : : type(mbsupercell_t), target, intent(in) :: supercell
88 : : type(rng_t), target, intent(in) :: rng
89 : 1 : call self%lattice_mover_t%initialize(params, supercell, rng)
90 : :
91 : : ! TODO: add friction
92 : 1 : self%fr = params%latt_friction
93 : :
94 [ - + + - : 3 : ABI_MALLOC(self%c3, (self%natom))
+ - + - -
+ - + -
+ ]
95 [ - + + - : 3 : ABI_MALLOC(self%c4, (self%natom))
+ - + - -
+ - + -
+ ]
96 [ - + + - : 3 : ABI_MALLOC(self%c5, (self%natom))
+ - + - -
+ - + -
+ ]
97 [ - + + - : 3 : ABI_MALLOC(self%xi, (3,self%natom))
+ - + - -
+ - + -
+ ]
98 [ - + + - : 3 : ABI_MALLOC(self%eta, (3,self%natom))
+ - + - -
+ - + -
+ ]
99 : :
100 : 1 : call self%update_vars()
101 : 1 : end subroutine initialize
102 : :
103 : :
104 : : !-------------------------------------------------------------------!
105 : : ! Finalize
106 : : !-------------------------------------------------------------------!
107 : 1 : subroutine finalize(self)
108 : : class(lattice_langevin_mover_t), intent(inout) :: self
109 [ - + ]: 1 : ABI_FREE(self%c3)
110 [ - + ]: 1 : ABI_FREE(self%c4)
111 [ - + ]: 1 : ABI_FREE(self%c5)
112 [ - + ]: 1 : ABI_FREE(self%xi)
113 [ - + ]: 1 : ABI_FREE(self%eta)
114 : 1 : call self%lattice_mover_t%finalize()
115 : 1 : end subroutine finalize
116 : :
117 : : !-------------------------------------------------------------------!
118 : : ! update_vars:
119 : : ! calculate c1 to c5 from dt, masses and temperature
120 : : ! It is called by the initialization function.
121 : : !-------------------------------------------------------------------!
122 : 1 : subroutine update_vars(self)
123 : : class(lattice_langevin_mover_t), intent(inout) :: self
124 : 2 : real(dp) :: dt, T, fr, sigma(self%natom)
125 : :
126 : 1 : dt=self%dt
127 : 1 : T= self%temperature
128 : 1 : fr=self%fr
129 [ + + ]: 321 : sigma(:) = sqrt(2.0*T*fr/self%masses(:))
130 : 1 : self%c1 = dt / 2.0 - dt * dt * fr / 8.0
131 : 1 : self%c2 = dt * fr / 2.0 - dt * dt * fr * fr / 8.0
132 [ + - - + : 322 : self%c3 = sqrt(dt) * sigma / 2.0 - dt**1.5 * fr * sigma / 8.0
- - + + ]
133 [ + - - + : 322 : self%c5 = dt**1.5 * sigma / (2.0 * sqrt(3.0))
- - + + ]
134 [ + - - + : 322 : self%c4 = fr / 2. * self%c5
- - + + ]
135 : 1 : end subroutine update_vars
136 : :
137 : : !-------------------------------------------------------------------!
138 : : !> @brief: Run_one_step using a Langevin heat bath.
139 : : !> effpot: the potential (which do the calculation of E and dE/dvar)
140 : : !> param[in]: effpot
141 : : !> param[in]: (optional) displacement
142 : : !> param[in]: (optional) strain
143 : : !> param[in]: (optional) spin
144 : : !> param[in]: (optional) lwf
145 : : ! NOTE: No need to pass the variable already saved in the mover.
146 : : ! e.g. For spin mover, do NOT pass the spin to it.
147 : : ! The other variables are only required if there is coupling with
148 : : ! the mover variable.
149 : : !-------------------------------------------------------------------!
150 [ - + - - : 10 : subroutine run_one_step(self, effpot,displacement, strain, spin, lwf , energy_table)
- - - + -
- - - ]
151 : : class(lattice_langevin_mover_t), intent(inout) :: self
152 : : class(abstract_potential_t), intent(inout) :: effpot
153 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
154 : : type(hash_table_t), optional, intent(inout) :: energy_table
155 : : integer :: i
156 : : character(len=40) :: key
157 : :
158 : : ! do not use displacement and strain because they are stored in the mover.
159 : 10 : ABI_UNUSED_A(displacement)
160 : 10 : ABI_UNUSED_A(strain)
161 : :
162 : 10 : self%energy = 0.0
163 [ + + + + ]: 12810 : self%forces(:,:) =0.0
164 : : call effpot%calculate( displacement=self%displacement, strain=self%strain, &
165 : : & spin=spin, lwf=lwf, force=self%forces, stress=self%stress, &
166 [ + - + - ]: 30 : & energy=self%energy, energy_table=energy_table)
167 : 10 : call self%rng%rand_normal_array(self%xi, 3*self%natom)
168 : 10 : call self%rng%rand_normal_array(self%eta, 3*self%natom)
169 : :
170 : :
171 : :
172 : : ! First half of velocity update
173 [ + + ]: 3210 : do i =1, self%natom
174 : : self%current_vcart(:,i) = self%current_vcart(:,i) + &
175 : : & self%c1 * self%forces(:,i) / self%masses(i) - &
176 : : & self%c2 * self%current_vcart(:,i) + &
177 [ + + ]: 12800 : & self%c3(i) * self%xi(:, i) - self%c4(i) * self%eta(:,i)
178 : :
179 : : self%displacement(:, i) = self%displacement(:, i) &
180 [ + + ]: 12810 : & + self%dt * self%current_vcart(:, i) + self%c5( i) *self%eta(:,i)
181 : : end do
182 : :
183 : : ! second half, update the velocity but not the displacement.
184 : 10 : self%energy=0.0
185 [ + + + + ]: 12810 : self%forces=0.0
186 : : call effpot%calculate( displacement=self%displacement, strain=self%strain, &
187 : : & spin=spin, lwf=lwf, force=self%forces, stress=self%stress, &
188 [ + - + - ]: 30 : & energy=self%energy, energy_table=energy_table)
189 [ + + ]: 3210 : do i =1, self%natom
190 : : self%current_vcart(:,i) = self%current_vcart(:,i) + &
191 : : &self%c1 * self%forces(:,i) / self%masses(i) - &
192 : : & self%c2 * self%current_vcart(:,i) + &
193 [ + + ]: 12810 : & self%c3(i) * self%xi(:, i) - self%c4(i) * self%eta(:,i)
194 : : end do
195 : :
196 : 10 : call self%get_T_and_Ek()
197 [ + - ]: 10 : if (present(energy_table)) then
198 : 10 : key = 'Lattice kinetic energy'
199 : 10 : call energy_table%put(key, self%Ek)
200 : : end if
201 : :
202 : 10 : end subroutine run_one_step
203 : :
204 [ + + - + : 4 : end module m_lattice_langevin_mover
+ + - + +
+ - + - +
+ + - + -
+ - + ]
205 : :
|
