Line data Source code
1 : !!****m* ABINIT/m_pred_hmc
2 : !! NAME
3 : !! m_pred_hmc
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 2017-2026 ABINIT group (SPr)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_pred_hmc
22 :
23 : implicit none
24 :
25 : private
26 : !!***
27 :
28 : public :: pred_hmc
29 : !!***
30 :
31 : contains
32 : !!***
33 :
34 : !!****f* ABINIT/pred_hmc
35 : !! NAME
36 : !! pred_hmc
37 : !!
38 : !! FUNCTION
39 : !! Hybrid Monte Carlo simulation algorithm. The routine generates a markov
40 : !! chain of structural configurations (states characterized by ionic positions
41 : !! and lattice parameters) with probability of observing a certian state
42 : !! equal to Gibbs statistical weight (exp(-etotal/kT)/Z).
43 : !!
44 : !! INPUTS
45 : !! ab_mover = Data structure containing information about
46 : !! input variables related to MD, e.g dtion, masses, etc.
47 : !! hist = history of ionic positions, forces,
48 : !! itime = index of current iteration
49 : !! icycle = index of current cycle of the iteration
50 : !! ntime = total number of iterations
51 : !! ncycle = total number of cycles
52 : !! zDEBUG = flag indicating whether to print debug info
53 : !! iexit = flag indicating finilization of mover loop
54 : !!
55 : !! OUTPUT
56 : !! hist = ionic positions, lattice parameters etc. are updated
57 : !!
58 : !! SIDE EFFECTS
59 : !!
60 : !! NOTES
61 : !!
62 : !! SOURCE
63 :
64 8014 : subroutine pred_hmc(ab_mover,hist,itime,icycle,ntime,ncycle,mttk_vars,zDEBUG,iexit)
65 :
66 : use defs_basis
67 : use m_errors
68 : use m_abicore
69 : use m_abimover
70 : use m_abihist
71 : use m_io_tools
72 : use m_hmc
73 :
74 : use m_geometry, only : xred2xcart
75 : use m_numeric_tools, only : uniformrandom
76 : use m_pred_velverlet, only : pred_velverlet
77 : use m_pred_isothermal, only : pred_isothermal
78 :
79 : !Arguments ------------------------------------
80 : type(abimover),intent(in) :: ab_mover
81 : type(abihist),intent(inout) :: hist
82 : type(mttk_type),intent(inout) :: mttk_vars
83 : integer,intent(in) :: itime
84 : integer,intent(in) :: icycle
85 : integer,intent(in) :: ntime
86 : integer,intent(in) :: ncycle
87 : integer,intent(in) :: iexit
88 : logical,intent(in) :: zDEBUG
89 :
90 : !Local variables-------------------------------
91 : integer,save :: seed ! seed for rnd generator
92 : integer :: iacc ! dummy integers for loop indexes and acceptance decision flag
93 : real(dp) :: etotal,epot,ekin,de ! total, potential (electronic), kinetic (ionic) energies and energy difference
94 : !real(dp) :: mv2tot,factor ! dummies used for rescaling of velocities
95 16028 : real(dp) :: xred(3,ab_mover%natom) ! reduced coordinates of all ions
96 16028 : real(dp) :: vel(3,ab_mover%natom) ! ionic velocities in Cartesian coordinates
97 : !real(dp) :: mvtot(3) ! total momentum of the cell used to rescale velocities
98 : real(dp) :: kbtemp !mtot, ! total ionic mass and target temperature in energy units
99 : real(dp) :: acell(3) ! lattice parameters
100 : real(dp) :: rprimd(3,3) ! lattice vectors
101 :
102 : real(dp),save :: etotal_hmc_prev,epot_hmc_prev ! total energy of the initial state
103 : real(dp),save :: strain(3,3),dstrain ! strain tensor
104 : real(dp),save :: rprimd_original(3,3) ! initial lattice vectors <= itime=1,icycle=1
105 : real(dp),allocatable,save :: xred_hmc_prev(:,:) ! reduced coordinates of the ions corresponding to the initial state
106 : real(dp),allocatable,save :: fcart_hmc_prev(:,:) ! reduced coordinates of the ions corresponding to the initial state
107 :
108 : logical,save :: strain_updated
109 : logical,save :: xred_updated
110 : integer,save :: strain_steps
111 : logical :: strain_sweep
112 :
113 : ! character(len=500) :: message
114 : ! *************************************************************************
115 :
116 : DBG_ENTER("COLL")
117 :
118 : ! if (option/=1 .and. option/=2 ) then
119 : ! write(msg,'(3a,i0)')&
120 : !& 'The argument option should be 1 or 2,',ch10,&
121 : !& 'however, option=',option
122 : ! ABI_BUG(msg)
123 : ! end if
124 : !
125 : ! if (sizein<1) then
126 : ! write(msg,'(3a,i0)')&
127 : !& 'The argument sizein should be a positive number,',ch10,&
128 : !& 'however, sizein=',sizein
129 : ! ABI_ERROR(msg)
130 : ! end if
131 :
132 : DBG_EXIT("COLL")
133 :
134 8014 : strain_sweep=.FALSE.
135 :
136 8014 : if(iexit/=0)then !icycle=ncycle and itime=ntime
137 2 : ABI_SFREE(xred_hmc_prev)
138 2 : ABI_SFREE(fcart_hmc_prev)
139 : !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! this is needed to deallocate vel_prev array allocated in pred_velverlet
140 2 : call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit)
141 : return
142 : end if
143 :
144 :
145 : !get current values of ionic positions and cell geometry and set up the target temperature
146 8012 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
147 :
148 34569932 : vel(:,:) = hist%vel(:,:,hist%ihist) ! velocities of all ions, not needed in reality
149 8012 : epot = hist%etot(hist%ihist) ! electronic sub-system energy, not needed
150 8012 : ekin = hist%ekin(hist%ihist) ! kinetic energy, not needed
151 :
152 8012 : kbtemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK ! correct temperature taking into account the possible heating/cooling
153 :
154 8012 : if(itime==1.and.icycle==1) then
155 2 : ABI_SFREE(xred_hmc_prev)
156 2 : ABI_SFREE(fcart_hmc_prev)
157 :
158 6 : ABI_MALLOC(xred_hmc_prev,(3,ab_mover%natom))
159 4 : ABI_MALLOC(fcart_hmc_prev,(3,ab_mover%natom))
160 :
161 2 : seed=-239
162 :
163 2 : rprimd_original(:,:)=rprimd(:,:)
164 2 : strain(:,:) = 0.0_dp
165 2 : strain_steps=0
166 2 : dstrain=0.001
167 :
168 2 : strain_updated=.FALSE.
169 2 : xred_updated=.FALSE.
170 :
171 2 : etotal_hmc_prev = zero
172 2 : epot_hmc_prev = zero
173 :
174 : end if
175 :
176 :
177 : !IN CASE THE SWEEP IS FOR UPDATE OF ATOMIC COORDINATES************************************************
178 : !if(.NOT.strain_sweep) then
179 :
180 : ! *---->*
181 : ! 1 n
182 :
183 8012 : if (icycle==1) then
184 :
185 204 : if(itime==1) then
186 2 : iacc=1
187 2 : etotal = epot + ekin
188 2 : de=zero
189 : else
190 202 : etotal = epot + ekin
191 202 : de = etotal - etotal_hmc_prev
192 202 : call metropolis_check(seed,de,kbtemp,iacc)
193 : !DEBUG
194 : ! write(std_out,*)' m_pred_hmc, after metropolis_check : seed,de,kbtemp,iacc=',seed,de,kbtemp,iacc
195 : !ENDDEBUG
196 : end if
197 :
198 204 : if(iacc==0)then !in case the new state is not accepted, then roll back the coordinates and energies
199 0 : xred(:,:)= xred_hmc_prev(:,:)
200 0 : epot = epot_hmc_prev
201 0 : hist%fcart(:,:,hist%ihist) = fcart_hmc_prev(:,:)
202 : else
203 864844 : xred_hmc_prev(:,:)=xred(:,:)
204 864844 : fcart_hmc_prev(:,:) = hist%fcart(:,:,hist%ihist)
205 204 : epot_hmc_prev = epot !update reference potential energy
206 : end if
207 :
208 : ! write(message,'(2a,i7,a,i2,a,E24.16,a,E24.16,a,E24.16)') ch10,' HMC Sweep => ',itime,' iacc= ', iacc,' epot= ',&
209 : !& epot,' ekin=',ekin,' de=',de
210 : ! call wrtout(ab_out,message,'COLL')
211 : ! call wrtout(std_out,message,'COLL')
212 :
213 : !call generate_random_velocities(ab_mover,kbtemp,seed,vel,ekin) ! this routine also computes the new kinetic energy
214 : !hist%vel(:,:,hist%ihist)=vel(:,:)
215 864844 : hist%vel(:,:,hist%ihist)=0
216 : !call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
217 : !etotal_hmc_prev=epot+ekin ! either old or current potential energy + new kinetic energy
218 :
219 : !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine
220 204 : call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit)
221 :
222 7808 : elseif(icycle > 1 .and. icycle <= ncycle)then !icycle/=1
223 :
224 : !call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine
225 7808 : call pred_isothermal(ab_mover,hist,icycle,mttk_vars,ncycle,zDEBUG,iexit)
226 :
227 : !end if
228 : !END OF ATOMIC COORDINATES SWEEP************************************************
229 :
230 : ! else if(icycle>ncycle) then ! strain update
231 : ! strain_updated = .TRUE.
232 : ! strain_steps = strain_steps + 1
233 : !! Metropolis update of lattice vectors and parameters in case optcell/=0
234 : ! if(icycle==ncycle+1.and.xred_updated) then
235 : ! !save rprimd_hmc_prev and total electronic energy etotal_hmc_prev
236 : ! call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG)
237 : ! etotal_hmc_prev = hist%etot(hist%ihist)
238 : ! strain_hmc_prev(:,:) = strain(:,:)
239 : !
240 : ! select case (ab_mover%optcell)
241 : ! case (1) !volume optimization only
242 : ! acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
243 : ! case (2,3,7,8,9) !full geometry optimization
244 : ! !suggest new strain tensor values
245 : ! do ii=1,3
246 : ! do jj=ii,3
247 : ! strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp)
248 : ! strain(jj,ii) = strain(ii,jj)
249 : ! enddo
250 : ! enddo
251 : ! if(ab_mover%optcell==3) then !eliminate volume change if optcell==3
252 : ! do ii=1,3
253 : ! strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3))
254 : ! enddo
255 : ! endif
256 : ! do jj=1,3 ! sum over three lattice vectors
257 : ! do ii=1,3 ! sum over Cart components
258 : ! rprimd(ii,jj)=rprimd_original(ii,jj)+&
259 : !& rprimd_original(1,jj)*strain(ii,1)+&
260 : !& rprimd_original(2,jj)*strain(ii,2)+&
261 : !& rprimd_original(3,jj)*strain(ii,3)
262 : ! enddo
263 : ! enddo
264 : ! if(ab_mover%optcell==7) then
265 : ! rprimd(:,1)=rprimd_original(:,1)
266 : ! else if (ab_mover%optcell==8) then
267 : ! rprimd(:,2)=rprimd_original(:,2)
268 : ! else if (ab_mover%optcell==9) then
269 : ! rprimd(:,3)=rprimd_original(:,3)
270 : ! endif
271 : ! case(4)
272 : ! acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
273 : ! case(5)
274 : ! acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
275 : ! case(6)
276 : ! acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
277 : ! !case default
278 : ! ! write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell
279 : ! ! ABI_ERROR(message)
280 : ! end select
281 : !
282 : ! !update the new suggested rprimd and or acell in the history record
283 : ! hist%ihist=abihist_findIndex(hist,+1)
284 : ! call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
285 : ! else
286 : !
287 : ! etotal = hist%etot(hist%ihist)
288 : ! de = etotal - etotal_hmc_prev
289 : !
290 : ! iacc=0
291 : ! rnd=uniformrandom(seed)
292 : ! if(de<0)then
293 : ! iacc=1
294 : ! else
295 : ! if(exp(-de/kbtemp)>rnd)then
296 : ! iacc=1
297 : ! end if
298 : ! end if
299 : !
300 : ! if(iacc==0) then
301 : ! strain(:,:)=strain_hmc_prev(:,:)
302 : ! acell(:)=acell_hmc_prev(:)
303 : ! else
304 : ! call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG)
305 : ! strain_hmc_prev(:,:) = strain(:,:)
306 : ! etotal_hmc_prev=etotal
307 : ! endif
308 : !
309 : ! !suggest new acell/rprimd values depending on the optcell value
310 : ! select case (ab_mover%optcell)
311 : ! case (1) !volume optimization only
312 : ! acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
313 : ! case (2,3,7,8,9) !full geometry optimization
314 : ! !suggest new strain tensor values
315 : ! do ii=1,3
316 : ! do jj=ii,3
317 : ! strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp)
318 : ! strain(jj,ii) = strain(ii,jj)
319 : ! enddo
320 : ! enddo
321 : ! if(ab_mover%optcell==3) then !eliminate volume change if optcell==3
322 : ! do ii=1,3
323 : ! strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3))
324 : ! enddo
325 : ! endif
326 : ! do jj=1,3 ! sum over three lattice vectors
327 : ! do ii=1,3 ! sum over Cart components
328 : ! rprimd(ii,jj)=rprimd_original(ii,jj)+&
329 : !& rprimd_original(1,jj)*strain(ii,1)+&
330 : !& rprimd_original(2,jj)*strain(ii,2)+&
331 : !& rprimd_original(3,jj)*strain(ii,3)
332 : ! enddo
333 : ! enddo
334 : ! if(ab_mover%optcell==7) then
335 : ! rprimd(:,1)=rprimd_original(:,1)
336 : ! else if (ab_mover%optcell==8) then
337 : ! rprimd(:,2)=rprimd_original(:,2)
338 : ! else if (ab_mover%optcell==9) then
339 : ! rprimd(:,3)=rprimd_original(:,3)
340 : ! endif
341 : ! case(4)
342 : ! acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
343 : ! case(5)
344 : ! acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
345 : ! case(6)
346 : ! acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
347 : ! case default
348 : ! ! write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell
349 : ! ! ABI_ERROR(message)
350 : ! end select
351 : !
352 : ! !update the new suggested rprimd/acell in the history record
353 : ! hist%ihist=abihist_findIndex(hist,+1)
354 : ! call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
355 : !
356 : ! endif
357 :
358 : end if
359 :
360 8014 : end subroutine pred_hmc
361 : !!***
362 :
363 : end module m_pred_hmc
364 : !!***
|
