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 :
|
