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