TABLE OF CONTENTS


ABINIT/m_pred_hmc [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_hmc

FUNCTION

COPYRIGHT

  Copyright (C) 2017-2018 ABINIT group (SPr)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

PARENTS

CHILDREN

SOURCE

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_pred_hmc
27 
28  implicit none
29 
30  private

ABINIT/pred_hmc [ Functions ]

[ Top ] [ Functions ]

NAME

  pred_hmc

FUNCTION

  Hybrid Monte Carlo simulation algorithm. The routine generates a markov
  chain of structural configurations (states characterized by ionic positions
  and lattice parameters) with probability of observing a certian state
  equal to Gibbs statistical weight (exp(-etotal/kT)/Z).

INPUTS

  ab_mover =  Data structure containing information about
              input variables related to MD, e.g dtion, masses, etc.
  hist     =  history of ionic positions, forces,
  itime    =  index of current iteration
  icycle   =  index of current cycle of the iteration
  ntime    =  total number of iterations
  ncycle   =  total number of cycles
  zDEBUG   =  flag indicating whether to print debug info
  iexit    =  flag indicating finilization of mover loop

OUTPUT

  hist  = ionic positions, lattice parameters etc. are updated

SIDE EFFECTS

NOTES

PARENTS

      mover

CHILDREN

      hist2var,pred_velverlet,var2hist,xred2xcart

SOURCE

 75 subroutine pred_hmc(ab_mover,hist,itime,icycle,ntime,ncycle,zDEBUG,iexit)
 76 
 77  use defs_basis
 78  use m_errors
 79  use m_abicore
 80  use m_abimover
 81  use m_abihist
 82  use m_io_tools
 83  use m_hmc
 84 
 85  use m_geometry,  only : xred2xcart
 86  use m_numeric_tools,  only : uniformrandom
 87  use m_pred_velverlet,     only : pred_velverlet
 88 
 89 !This section has been created automatically by the script Abilint (TD).
 90 !Do not modify the following lines by hand.
 91 #undef ABI_FUNC
 92 #define ABI_FUNC 'pred_hmc'
 93 !End of the abilint section
 94 
 95  implicit none
 96 
 97 !Arguments ------------------------------------
 98  type(abimover),intent(in)   :: ab_mover
 99  type(abihist),intent(inout) :: hist
100  integer,intent(in)          :: itime
101  integer,intent(in)          :: icycle
102  integer,intent(in)          :: ntime
103  integer,intent(in)          :: ncycle
104  integer,intent(in)          :: iexit
105  logical,intent(in)          :: zDEBUG
106 
107 !Local variables-------------------------------
108 
109  integer,save  :: seed                                   ! seed for rnd generator
110  integer       :: ii,jj,iacc                             ! dummy integers for loop indexes and acceptance decision flag
111  real(dp)      :: etotal,epot,ekin,de                    ! total, potential (electronic), kinetic (ionic) energies and energy difference
112  real(dp)      :: mv2tot,factor                          ! dummies used for rescaling of velocities
113  real(dp)      :: rnd
114  real(dp)      :: xred(3,ab_mover%natom)                 ! reduced coordinates of all ions
115  real(dp)      :: vel(3,ab_mover%natom)                  ! ionic velocities in Cartesian coordinates
116  real(dp)      :: mvtot(3)                               ! total momentum of the cell used to rescale velocities
117  real(dp)      :: mtot,kbtemp                            ! total ionic mass and target temperature in energy units
118  real(dp)      :: acell(3)                               ! lattice parameters
119  real(dp)      :: rprimd(3,3)                            ! lattice vectors
120 
121  real(dp),save :: etotal_hmc_prev,epot_hmc_prev          ! total energy of the initial state
122  real(dp),save :: acell_hmc_prev(3)                      !
123  real(dp),save :: rprimd_hmc_prev(3,3)                   !
124  real(dp),save :: strain_hmc_prev(3,3)                   !
125  real(dp),save :: strain(3,3),dstrain                    ! strain tensor
126  real(dp),save :: rprimd_original(3,3)                   ! initial lattice vectors <= itime=1,icycle=1
127  real(dp),allocatable,save :: xred_hmc_prev(:,:)         ! reduced coordinates of the ions corresponding to the initial state
128  real(dp),allocatable,save :: fcart_hmc_prev(:,:)        ! reduced coordinates of the ions corresponding to the initial state
129 
130  logical,save  :: strain_updated
131  logical,save  :: xred_updated
132  integer,save  :: strain_steps
133  logical       :: strain_sweep
134 
135  character(len=500) :: message
136 ! *************************************************************************
137 
138  DBG_ENTER("COLL")
139 
140 ! if (option/=1 .and. option/=2 ) then
141 !   write(msg,'(3a,i0)')&
142 !&   'The argument option should be 1 or 2,',ch10,&
143 !&   'however, option=',option
144 !   MSG_BUG(msg)
145 ! end if
146 !
147 ! if (sizein<1) then
148 !   write(msg,'(3a,i0)')&
149 !&   'The argument sizein should be a positive number,',ch10,&
150 !&   'however, sizein=',sizein
151 !   MSG_ERROR(msg)
152 ! end if
153 
154  DBG_EXIT("COLL")
155 
156  strain_sweep=.FALSE.
157 
158  if(iexit/=0)then  !icycle=ncycle and itime=ntime
159    if (allocated(xred_hmc_prev))  then
160      ABI_DEALLOCATE(xred_hmc_prev)
161    end if
162    if (allocated(fcart_hmc_prev))  then
163      ABI_DEALLOCATE(fcart_hmc_prev)
164    end if
165    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
166    return
167  end if
168 
169 
170  !get current values of ionic positions and cell geometry and set up the target temperature
171  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
172 
173  vel(:,:)   = hist%vel(:,:,hist%ihist)                ! velocities of all ions, not needed in reality
174  epot       = hist%etot(hist%ihist)                   ! electronic sub-system energy, not needed
175  ekin       = hist%ekin(hist%ihist)                   ! kinetic energy, not needed
176 
177  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
178 
179  if(itime==1.and.icycle==1) then
180    if (allocated(xred_hmc_prev))  then
181      ABI_DEALLOCATE(xred_hmc_prev)
182    end if
183    if (allocated(fcart_hmc_prev))  then
184      ABI_DEALLOCATE(fcart_hmc_prev)
185    end if
186 
187    ABI_ALLOCATE(xred_hmc_prev,(3,ab_mover%natom))
188    ABI_ALLOCATE(fcart_hmc_prev,(3,ab_mover%natom))
189 
190    seed=-239
191 
192    rprimd_original(:,:)=rprimd(:,:)
193    strain(:,:) = 0.0_dp
194    strain_steps=0
195    dstrain=0.001
196 
197    strain_updated=.FALSE.
198    xred_updated=.FALSE.
199  end if
200 
201 
202  !IN CASE THE SWEEP IS FOR UPDATE OF ATOMIC COORDINATES************************************************
203  !if(.NOT.strain_sweep) then
204 
205    ! *---->*
206    ! 1     n
207 
208  if (icycle==1) then
209 
210    if(itime==1) then
211      iacc=1
212      etotal = epot + ekin
213      de=zero
214    else
215      etotal = epot + ekin
216      de = etotal - etotal_hmc_prev
217      call metropolis_check(seed,de,kbtemp,iacc)
218 !DEBUG
219      write(std_out,*)' m_pred_hmc, after metropolis_check : seed,de,kbtemp,iacc=',seed,de,kbtemp,iacc
220 !ENDDEBUG
221    end if
222 
223    if(iacc==0)then  !in case the new state is not accepted, then roll back the coordinates and energies
224      xred(:,:)= xred_hmc_prev(:,:)
225      epot     = epot_hmc_prev
226      hist%fcart(:,:,hist%ihist) = fcart_hmc_prev(:,:)
227    else
228      xred_hmc_prev(:,:)=xred(:,:)
229      fcart_hmc_prev(:,:) = hist%fcart(:,:,hist%ihist)
230      epot_hmc_prev   = epot         !update reference potential energy
231    end if
232 
233    write(message,'(2a,i7,a,i2,a,E24.16,a,E24.16,a,E24.16)') ch10,' HMC Sweep => ',itime,' iacc= ', iacc,' epot= ',&
234 &   epot,' ekin=',ekin,' de=',de
235    call wrtout(ab_out,message,'COLL')
236    call wrtout(std_out,message,'COLL')
237 
238    call generate_random_velocities(ab_mover,kbtemp,seed,vel,ekin)  ! this routine also computes the new kinetic energy
239    hist%vel(:,:,hist%ihist)=vel(:,:)
240    call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
241    etotal_hmc_prev=epot+ekin ! either old or current potential energy + new kinetic energy
242 
243    call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine
244 
245  else !icycle/=1
246 
247    call pred_velverlet(ab_mover,hist,itime,ntime,zDEBUG,iexit,1,icycle,ncycle) ! 1 is indicating that velverlet is called from hmc routine
248 
249  end if
250 
251  !end if
252  !END OF ATOMIC COORDINATES SWEEP************************************************
253 
254 ! else if(icycle>ncycle) then ! strain update
255 !   strain_updated = .TRUE.
256 !   strain_steps   = strain_steps + 1
257  ! Metropolis update of lattice vectors and parameters in case optcell/=0
258 !   if(icycle==ncycle+1.and.xred_updated) then
259 !     !save rprimd_hmc_prev and total electronic energy etotal_hmc_prev
260 !     call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG)
261 !     etotal_hmc_prev = hist%etot(hist%ihist)
262 !     strain_hmc_prev(:,:) = strain(:,:)
263 
264 !     select case (ab_mover%optcell)
265 !     case (1) !volume optimization only
266 !       acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
267 !     case (2,3,7,8,9) !full geometry optimization
268 !       !suggest new strain tensor values
269 !       do ii=1,3
270 !         do jj=ii,3
271 !           strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp)
272 !           strain(jj,ii) = strain(ii,jj)
273 !         enddo
274 !       enddo
275 !       if(ab_mover%optcell==3) then !eliminate volume change if optcell==3
276 !         do ii=1,3
277 !           strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3))
278 !         enddo
279 !       endif
280 !       do jj=1,3    ! sum over three lattice vectors
281 !         do ii=1,3  ! sum over Cart components
282 !           rprimd(ii,jj)=rprimd_original(ii,jj)+&
283 !&                        rprimd_original(1,jj)*strain(ii,1)+&
284 !&                        rprimd_original(2,jj)*strain(ii,2)+&
285 !&                        rprimd_original(3,jj)*strain(ii,3)
286 !         enddo
287 !       enddo
288 !       if(ab_mover%optcell==7) then
289 !         rprimd(:,1)=rprimd_original(:,1)
290 !       else if (ab_mover%optcell==8) then
291 !         rprimd(:,2)=rprimd_original(:,2)
292 !       else if (ab_mover%optcell==9) then
293 !         rprimd(:,3)=rprimd_original(:,3)
294 !       endif
295 !     case(4)
296 !       acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
297 !     case(5)
298 !       acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
299 !     case(6)
300 !       acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
301 !     case default
302 !     !  write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell
303 !     !  MSG_ERROR(message)
304 !     end select
305 !
306 !     !update the new suggested rprimd and or acell in the history record
307 !     hist%ihist=abihist_findIndex(hist,+1)
308 !     call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
309 !   else
310 !
311 !     etotal = hist%etot(hist%ihist)
312 !     de = etotal - etotal_hmc_prev
313 !
314 !     iacc=0
315 !     rnd=uniformrandom(seed)
316 !     if(de<0)then
317 !       iacc=1
318 !     else
319 !       if(exp(-de/kbtemp)>rnd)then
320 !         iacc=1
321 !       end if
322 !     end if
323 
324 !     if(iacc==0) then
325 !      strain(:,:)=strain_hmc_prev(:,:)
326 !      acell(:)=acell_hmc_prev(:)
327 !     else
328 !      call hist2var(acell_hmc_prev,hist,ab_mover%natom,rprimd_hmc_prev,xred,zDEBUG)
329 !      strain_hmc_prev(:,:) = strain(:,:)
330 !      etotal_hmc_prev=etotal
331 !     endif
332 
333      !suggest new acell/rprimd values depending on the optcell value
334 !     select case (ab_mover%optcell)
335 !     case (1) !volume optimization only
336 !       acell(:)=acell(:)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
337 !     case (2,3,7,8,9) !full geometry optimization
338 !       !suggest new strain tensor values
339 !       do ii=1,3
340 !         do jj=ii,3
341 !           strain(ii,jj) = strain(ii,jj)+ 2.0_dp*dstrain*(uniformrandom(seed)-0.5_dp)
342 !           strain(jj,ii) = strain(ii,jj)
343 !         enddo
344 !       enddo
345 !       if(ab_mover%optcell==3) then !eliminate volume change if optcell==3
346 !         do ii=1,3
347 !           strain(ii,ii) = strain(ii,ii) -(strain(1,1)+strain(2,2)+strain(3,3))
348 !         enddo
349 !       endif
350 !       do jj=1,3    ! sum over three lattice vectors
351 !         do ii=1,3  ! sum over Cart components
352 !           rprimd(ii,jj)=rprimd_original(ii,jj)+&
353 !&                        rprimd_original(1,jj)*strain(ii,1)+&
354 !&                        rprimd_original(2,jj)*strain(ii,2)+&
355 !&                        rprimd_original(3,jj)*strain(ii,3)
356 !         enddo
357 !       enddo
358 !       if(ab_mover%optcell==7) then
359 !         rprimd(:,1)=rprimd_original(:,1)
360 !       else if (ab_mover%optcell==8) then
361 !         rprimd(:,2)=rprimd_original(:,2)
362 !       else if (ab_mover%optcell==9) then
363 !         rprimd(:,3)=rprimd_original(:,3)
364 !       endif
365 !     case(4)
366 !       acell(1)=acell(1)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
367 !     case(5)
368 !       acell(2)=acell(2)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
369 !     case(6)
370 !       acell(3)=acell(3)*(1.0_dp+dstrain*2.0_dp*(uniformrandom(seed)-0.5_dp))
371 !     case default
372 !     !  write(message,"(a,i0)") "Wrong value of optcell: ",ab_mover%optcell
373 !     !  MSG_ERROR(message)
374 !     end select
375 !
376 !     !update the new suggested rprimd/acell in the history record
377 !     hist%ihist=abihist_findIndex(hist,+1)
378 !     call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
379 !
380 !   endif
381 !
382 ! end if
383 
384 end subroutine pred_hmc