Branch data 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 : : !!***
|
