Branch data Line data Source code
1 : : !!****m* ABINIT/m_lattice_berendsen_NPT_mover
2 : : !! TODO: This is not yet implemented.
3 : : !! NAME
4 : : !! m_lattice_berendsen_NPT_mover
5 : : !!
6 : : !! FUNCTION
7 : : !! This module contains the berendsen (NPT) 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_NPT_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_NPT_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_NPT_mover_t
61 : : real(dp) :: taut ! the characteristic time of the relaxation of velocity.
62 : : ! it is usually larger than the time step.
63 : : real(dp) :: taup ! the characteristic time of the relaxation of pressure.
64 : : real(dp) :: compressibility
65 : : contains
66 : : procedure :: initialize
67 : : procedure :: finalize
68 : : procedure :: run_one_step
69 : : procedure :: scale_velocities
70 : : end type lattice_berendsen_NPT_mover_t
71 : :
72 : : contains
73 : :
74 : :
75 : :
76 : 0 : subroutine initialize(self,params, supercell, rng)
77 : : class(lattice_berendsen_NPT_mover_t), intent(inout) :: self
78 : : type(multibinit_dtset_type), target, intent(in):: params
79 : : type(mbsupercell_t), target, intent(in) :: supercell
80 : : type(rng_t), target, intent(in) :: rng
81 : 0 : self%taut = params%latt_taut
82 : : !self%taup = params%latt_taup
83 : : !self%compressibility =params%latt_compressibility
84 : 0 : call self%lattice_mover_t%initialize(params, supercell, rng)
85 : 0 : ABI_ERROR("The Berendsen NPT mover has not yet been implemented")
86 : : !TODO: Implement
87 : 0 : end subroutine initialize
88 : :
89 : :
90 : 0 : subroutine finalize(self)
91 : : class(lattice_berendsen_NPT_mover_t), intent(inout) :: self
92 : 0 : call self%lattice_mover_t%finalize()
93 : 0 : end subroutine finalize
94 : :
95 : :
96 : : !-------------------------------------------------------------------!
97 : : ! scale_velocities:
98 : : ! scale the velocities so that they get close to the required temperture
99 : : !
100 : : !-------------------------------------------------------------------!
101 : 0 : subroutine scale_velocities(self)
102 : : class(lattice_berendsen_NPT_mover_t), intent(inout) :: self
103 : : real(dp) :: tautscl, old_temperature, scale_temperature, tmp
104 : 0 : tautscl = self%dt / self%taut
105 : 0 : old_temperature=self%T_ob
106 : 0 : tmp=1.0 +(self%temperature / old_temperature - 1.0) * tautscl
107 [ # # ]: 0 : if(tmp< 0.0) then
108 : 0 : ABI_ERROR("The time scale for the Berendsen Algorithm should be at least larger than dtion.")
109 : : else
110 : 0 : scale_temperature=sqrt(tmp)
111 : : end if
112 : : ! Limit the velocity scaling to reasonable values
113 [ # # ]: 0 : if( scale_temperature > 1.1) then
114 : : scale_temperature = 1.1
115 : : elseif (scale_temperature < 0.9) then
116 : : scale_temperature = 0.9
117 : : endif
118 [ # # # # ]: 0 : self%current_vcart(:,:) = self%current_vcart(:,:) * scale_temperature
119 : 0 : end subroutine scale_velocities
120 : :
121 : :
122 : : !-------------------------------------------------------------------!
123 : : ! run_one_step.
124 : : ! The algorithm is almost the same as the velocity verlet algorithm,
125 : : ! except at the begining, the velocities are scaled so that the temperature
126 : : ! is getting closer to the required temperature.
127 : : !-------------------------------------------------------------------!
128 [ # # # # : 0 : subroutine run_one_step(self, effpot,displacement, strain, spin, lwf, energy_table)
# # # # #
# # # ]
129 : : class(lattice_berendsen_NPT_mover_t), intent(inout) :: self
130 : : class(abstract_potential_t), intent(inout) :: effpot
131 : : real(dp), optional, intent(inout) :: displacement(:,:), strain(:,:), spin(:,:), lwf(:)
132 : : type(hash_table_t), optional, intent(inout) :: energy_table
133 : : integer :: i
134 : : character(len=40) :: key
135 : :
136 : : ABI_UNUSED(displacement)
137 : : ABI_UNUSED(strain)
138 : :
139 : : ! scale the velocity.
140 : 0 : call self%scale_velocities()
141 : :
142 : 0 : self%energy=0.0
143 [ # # # # ]: 0 : self%forces(:,:) =0.0
144 : : call effpot%calculate( displacement=self%displacement, strain=self%strain, &
145 : : & spin=spin, lwf=lwf, force=self%forces, stress=self%stress, &
146 [ # # # # ]: 0 : & energy=self%energy, energy_table=energy_table)
147 [ # # ]: 0 : do i=1, self%natom
148 : : self%current_vcart(:,i) = self%current_vcart(:,i) + &
149 [ # # ]: 0 : & (0.5_dp * self%dt) * self%forces(:,i)/self%masses(i)
150 : : end do
151 : 0 : call self%force_stationary()
152 [ # # # # ]: 0 : self%displacement(:,:) = self%displacement(:,:)+self%current_vcart(:,:) * self%dt
153 : :
154 : :
155 : : ! second half of velocity update.
156 : : ! v(t+dt) = v(t + 1/2 dt) + F/m * 1/2 dt
157 : : ! NOTE: energy and forces should be initialized before every calculation!
158 : 0 : self%energy=0.0
159 [ # # # # ]: 0 : self%forces(:,:)=0.0
160 : : call effpot%calculate( displacement=self%displacement, &
161 : : & strain=self%strain, spin=spin, lwf=lwf, force=self%forces, &
162 [ # # # # ]: 0 : & stress=self%stress, energy=self%energy, energy_table=energy_table)
163 [ # # ]: 0 : do i=1, self%natom
164 : : self%current_vcart(:,i) = self%current_vcart(:,i) &
165 [ # # ]: 0 : & + (0.5_dp * self%dt) * self%forces(:,i)/self%masses(i)
166 : : end do
167 : 0 : call self%force_stationary()
168 : :
169 : 0 : call self%get_T_and_Ek()
170 [ # # ]: 0 : if (present(energy_table)) then
171 : 0 : key = 'Lattice kinetic energy'
172 : 0 : call energy_table%put(key, self%Ek)
173 : : end if
174 : :
175 : :
176 : 0 : end subroutine run_one_step
177 : :
178 [ # # # # : 0 : end module m_lattice_berendsen_NPT_mover
# # # # #
# # # ]
179 : :
|
