Branch data Line data Source code
1 : : !!****m* ABINIT/m_lattice_berendsen_NVT_mover
2 : : !! TODO: This is not yet implemented.
3 : : !! NAME
4 : : !! m_lattice_berendsen_NVT_mover
5 : : !!
6 : : !! FUNCTION
7 : : !! This module contains the berendsen (NVT) lattice mover.
8 : : !! The method is described in
9 : : !! H.J.C. Berendsen, J.P.M. Postma, A. DiNola, and J.R. Haak,
10 : : !! "Molecular dynamics with coupling to an external bath,"
11 : : !! J. Chem. Phys., 81 3684-3690 (1984)
12 : : !! NOTE: that this method does NOT generate properly the thermostated
13 : : !! ensemble. It does not have the correct distribution of the kinetic energy.
14 : : !! However, it approches the target temperature exponentially without oscillation,
15 : : !! for which the steps can be easily controlled.
16 : : !!
17 : : !! Datatypes:
18 : : !!
19 : : !! * lattice_berendsen_NVT_mover_t: defines the lattice movers
20 : : !!
21 : : !! Subroutines:
22 : : !! TODO: add this when F2003 doc style is determined.
23 : : !!
24 : : !!
25 : : !! COPYRIGHT
26 : : !! Copyright (C) 2001-2026 ABINIT group (hexu)
27 : : !! This file is distributed under the terms of the
28 : : !! GNU General Public License, see ~abinit/COPYING
29 : : !! or http://www.gnu.org/copyleft/gpl.txt .
30 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
31 : : !!
32 : : !! SOURCE
33 : :
34 : :
35 : :
36 : : #if defined HAVE_CONFIG_H
37 : : #include "config.h"
38 : : #endif
39 : :
40 : : #include "abi_common.h"
41 : :
42 : : module m_lattice_berendsen_NVT_mover
43 : : use defs_basis
44 : : use m_abicore
45 : : use m_errors
46 : :
47 : : use m_multibinit_dataset, only: multibinit_dtset_type
48 : : use m_abstract_potential, only: abstract_potential_t
49 : : use m_abstract_mover, only: abstract_mover_t
50 : : use m_lattice_mover, only: lattice_mover_t
51 : : use m_multibinit_cell, only: mbcell_t, mbsupercell_t
52 : : use m_random_xoroshiro128plus, only: rng_t
53 : : use m_hashtable_strval, only: hash_table_t
54 : : !!***
55 : :
56 : : implicit none
57 : :
58 : : private
59 : :
60 : : type, public, extends(lattice_mover_t) :: lattice_berendsen_NVT_mover_t
61 : : real(dp) :: taut ! the characteristic time of the relaxation of velocity.
62 : : ! usually larger than time step
63 : : contains
64 : : procedure :: initialize
65 : : procedure :: finalize
66 : : procedure :: run_one_step
67 : : procedure :: scale_velocities
68 : : end type lattice_berendsen_NVT_mover_t
69 : :
70 : : contains
71 : :
72 : :
73 : : !----------------------------------------------------------------------
74 : : !> @brief initialize
75 : : !>
76 : : !> @param[in] params: input parameters
77 : : !> @param[in] supercell: superell
78 : : !> @param[in] rng: random number generator
79 : : !----------------------------------------------------------------------
80 : 1 : subroutine initialize(self,params, supercell, rng)
81 : : class(lattice_berendsen_NVT_mover_t), intent(inout) :: self
82 : : type(multibinit_dtset_type), target, intent(in):: params
83 : : type(mbsupercell_t), target, intent(in) :: supercell
84 : : type(rng_t), target, intent(in) :: rng
85 : 1 : self%taut = params%latt_taut
86 : 1 : call self%lattice_mover_t%initialize(params, supercell, rng)
87 : 1 : end subroutine initialize
88 : :
89 : :
90 : : !----------------------------------------------------------------------
91 : : !> @brief finalize
92 : : !----------------------------------------------------------------------
93 : 1 : subroutine finalize(self)
94 : : class(lattice_berendsen_NVT_mover_t), intent(inout) :: self
95 : 1 : call self%lattice_mover_t%finalize()
96 : 1 : end subroutine finalize
97 : :
98 : :
99 : : !-------------------------------------------------------------------!
100 : : ! scale_velocities:
101 : : ! scale the velocities so that they get close to the required temperture
102 : : !
103 : : !-------------------------------------------------------------------!
104 : 1000 : subroutine scale_velocities(self)
105 : : class(lattice_berendsen_NVT_mover_t), intent(inout) :: self
106 : : real(dp) :: tautscl, old_temperature, scale_temperature, tmp
107 : 1000 : tautscl = self%dt / self%taut
108 : 1000 : old_temperature=self%T_ob
109 : 1000 : tmp=1.0 +(self%temperature / old_temperature - 1.0) * tautscl
110 [ - + ]: 1000 : if(tmp< 0.0) then
111 : 0 : ABI_ERROR("The time scale for the Berendsen algorithm should be at least larger than dtion")
112 : : else
113 : 1000 : scale_temperature=sqrt(tmp)
114 : : end if
115 : : ! Limit the velocity scaling to reasonable values
116 [ + - ]: 1000 : if( scale_temperature > 1.1) then
117 : : scale_temperature = 1.1
118 : : elseif (scale_temperature < 0.9) then
119 : : scale_temperature = 0.9
120 : : endif
121 [ + + + + ]: 4321000 : self%current_vcart(:,:) = self%current_vcart(:,:) * scale_temperature
122 : 1000 : end subroutine scale_velocities
123 : :
124 : :
125 : : !-------------------------------------------------------------------!
126 : : ! run_one_step.
127 : : ! The algorithm is almost the same as the velocity verlet algorithm,
128 : : ! except at the begining, the velocities are scaled so that the temperature
129 : : ! is getting closer to the required temperature.
130 : : !-------------------------------------------------------------------!
131 [ - + - - : 1000 : subroutine run_one_step(self, effpot,displacement, strain, spin, lwf, energy_table)
- - - + -
- - - ]
132 : : class(lattice_berendsen_NVT_mover_t), intent(inout) :: self
133 : : class(abstract_potential_t), intent(inout) :: effpot
134 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
135 : : type(hash_table_t), optional, intent(inout) :: energy_table
136 : : integer :: i
137 : : character(len=40) :: key
138 : :
139 : :
140 : : ABI_UNUSED(displacement)
141 : : ABI_UNUSED(strain)
142 : :
143 : : ! scale the velocity.
144 : 1000 : call self%scale_velocities()
145 : :
146 : 1000 : self%energy=0.0
147 [ + + + + ]: 4321000 : self%forces(:,:) =0.0
148 : : call effpot%calculate( displacement=self%displacement, strain=self%strain, &
149 : : & spin=spin, lwf=lwf, force=self%forces, stress=self%stress, &
150 [ + - + - ]: 3000 : & energy=self%energy, energy_table=energy_table)
151 [ + + ]: 1081000 : do i=1, self%natom
152 : : self%current_vcart(:,i) = self%current_vcart(:,i) + &
153 [ + + ]: 4321000 : & (0.5_dp * self%dt) * self%forces(:,i)/self%masses(i)
154 : : end do
155 : 1000 : call self%force_stationary()
156 [ + + + + ]: 4321000 : self%displacement(:,:) = self%displacement(:,:)+self%current_vcart(:,:) * self%dt
157 : :
158 : :
159 : : ! second half of velocity update.
160 : : ! v(t+dt) = v(t + 1/2 dt) + F/m * 1/2 dt
161 : : ! NOTE: energy and forces should be initialized before every calculation!
162 : 1000 : self%energy=0.0
163 [ + + + + ]: 4321000 : self%forces(:,:)=0.0
164 : : call effpot%calculate( displacement=self%displacement, &
165 : : & strain=self%strain, spin=spin, lwf=lwf, force=self%forces, &
166 [ + - + - ]: 3000 : & stress=self%stress, energy=self%energy, energy_table=energy_table)
167 [ + + ]: 1081000 : do i=1, self%natom
168 : : self%current_vcart(:,i) = self%current_vcart(:,i) &
169 [ + + ]: 4321000 : & + (0.5_dp * self%dt) * self%forces(:,i)/self%masses(i)
170 : : end do
171 : 1000 : call self%force_stationary()
172 : 1000 : call self%get_T_and_Ek()
173 [ + - ]: 1000 : if (present(energy_table)) then
174 : 1000 : key = 'Lattice kinetic energy'
175 : 1000 : call energy_table%put(key, self%Ek)
176 : : end if
177 : 1000 : end subroutine run_one_step
178 : :
179 [ + - - + : 1 : end module m_lattice_berendsen_NVT_mover
- + - + -
+ - + ]
180 : :
|
