TABLE OF CONTENTS


ABINIT/m_pred_langevin [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_langevin

FUNCTION

COPYRIGHT

  Copyright (C) 1998-2018 ABINIT group (DCA, XG, GMR, JCC, SE)
  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

21 #if defined HAVE_CONFIG_H
22 #include "config.h"
23 #endif
24 
25 #include "abi_common.h"
26 
27 module m_pred_langevin
28 
29  use defs_basis
30  use m_abicore
31  use m_abimover
32  use m_abihist
33 
34  use m_numeric_tools,  only : uniformrandom
35  use m_geometry,    only : xcart2xred, xred2xcart, metric
36  use m_results_gs , only : results_gs_type
37 
38  implicit none
39 
40  private

ABINIT/pred_langevin [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_langevin

FUNCTION

 Ionmov predictors (9) Langevin dynamics algorithm

 IONMOV 9:
 Uses a Langevin dynamics algorithm :
 see J. Chelikowsky, J. Phys. D : Appl Phys. 33(2000)R33 [[cite:Chelikowsky2000]]

INPUTS

 ab_mover <type(abimover)> : Datatype with all the information
                                needed by the preditor
 itime  : Index of the present iteration
 ntime  : Maximal number of iterations
 icycle : Index of the present cycle
 ncycle : Maximal number of cycles
 zDEBUG : if true print some debugging information

OUTPUT

SIDE EFFECTS

 hist <type(abihist)> : History of positions,forces acell, rprimd, stresses

PARENTS

      mover

CHILDREN

      hist2var,metric,var2hist,wrtout,xcart2xred,xred2xcart

SOURCE

 82 subroutine pred_langevin(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit,skipcycle)
 83 
 84 
 85 !This section has been created automatically by the script Abilint (TD).
 86 !Do not modify the following lines by hand.
 87 #undef ABI_FUNC
 88 #define ABI_FUNC 'pred_langevin'
 89 !End of the abilint section
 90 
 91  implicit none
 92 
 93 !Arguments ------------------------------------
 94 !scalars
 95  type(abimover),intent(in)       :: ab_mover
 96  type(abihist),intent(inout) :: hist
 97  integer,intent(in)    :: itime
 98  integer,intent(in)    :: ntime
 99  integer,intent(in)    :: iexit
100  integer,intent(in)    :: icycle
101  integer,intent(inout) :: ncycle
102  logical,intent(in)    :: zDEBUG
103  logical,intent(out)   :: skipcycle
104 
105 !Local variables-------------------------------
106 !scalars
107  integer  :: ii,kk,iatom,idim,iatom1,iatom2,itypat,idum=5,ihist_prev,mcfac
108  real(dp) :: ucvol,ucvol_next
109  real(dp),parameter :: v2tol=tol8
110  real(dp) :: etotal,rescale_vel,ran_num1,ran_num2
111  real(dp) :: ktemp,dist,distx,disty,distz,maxp1,maxp2,v2nose
112  real(dp) :: sig_gauss,delxi
113  logical  :: jump_end_of_cycle=.FALSE.
114  character(len=5000) :: message
115 !arrays
116  integer,allocatable :: imax_perm(:)
117  real(dp),allocatable,save :: max_perm(:),pot_perm(:)
118  real(dp),allocatable,save :: ran_force(:,:),lang_force(:,:)
119  real(dp),allocatable,save :: fcart_mold(:,:),fcart_m(:,:)
120 
121  real(dp) :: acell(3),acell_next(3)
122  real(dp) :: rprim(3,3),rprimd(3,3),rprimd_next(3,3),rprim_next(3,3)
123  real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3)
124  real(dp) :: fcart(3,ab_mover%natom)
125 !real(dp) :: fred_corrected(3,ab_mover%natom)
126  real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
127  real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
128  real(dp) :: vel(3,ab_mover%natom)
129  real(dp) :: strten(6)
130 
131 !***************************************************************************
132 !Beginning of executable session
133 !***************************************************************************
134 
135  if(iexit/=0)then
136    if (allocated(pot_perm))    then
137      ABI_DEALLOCATE(pot_perm)
138    end if
139    if (allocated(max_perm))    then
140      ABI_DEALLOCATE(max_perm)
141    end if
142    if (allocated(imax_perm))   then
143      ABI_DEALLOCATE(imax_perm)
144    end if
145    if (allocated(ran_force))   then
146      ABI_DEALLOCATE(ran_force)
147    end if
148    if (allocated(lang_force))  then
149      ABI_DEALLOCATE(lang_force)
150    end if
151    if (allocated(fcart_mold))  then
152      ABI_DEALLOCATE(fcart_mold)
153    end if
154    if (allocated(fcart_m))     then
155      ABI_DEALLOCATE(fcart_m)
156    end if
157    return
158  end if
159 
160  jump_end_of_cycle=.FALSE.
161 
162 !write(std_out,*) 'langevin 02',jump_end_of_cycle
163 !##########################################################
164 !### 02. Allocate the arrays
165 !###     These arrays could be allocated from a previus
166 !###     dataset that exit before itime==ntime
167 
168  if(itime==1)then
169    if (allocated(pot_perm))    then
170      ABI_DEALLOCATE(pot_perm)
171    end if
172    if (allocated(max_perm))    then
173      ABI_DEALLOCATE(max_perm)
174    end if
175    if (allocated(imax_perm))   then
176      ABI_DEALLOCATE(imax_perm)
177    end if
178    if (allocated(ran_force))   then
179      ABI_DEALLOCATE(ran_force)
180    end if
181    if (allocated(lang_force))  then
182      ABI_DEALLOCATE(lang_force)
183    end if
184    if (allocated(fcart_mold))  then
185      ABI_DEALLOCATE(fcart_mold)
186    end if
187    if (allocated(fcart_m))     then
188      ABI_DEALLOCATE(fcart_m)
189    end if
190  end if
191 
192  if (.not.allocated(pot_perm))    then
193    ABI_ALLOCATE(pot_perm,(ab_mover%natom))
194  end if
195  if (.not.allocated(max_perm))    then
196    ABI_ALLOCATE(max_perm,(ab_mover%ntypat))
197  end if
198  if (.not.allocated(imax_perm))   then
199    ABI_ALLOCATE(imax_perm,(ab_mover%ntypat))
200  end if
201  if (.not.allocated(ran_force))   then
202    ABI_ALLOCATE(ran_force,(3,ab_mover%natom))
203  end if
204  if (.not.allocated(lang_force))  then
205    ABI_ALLOCATE(lang_force,(3,ab_mover%natom))
206  end if
207  if (.not.allocated(fcart_mold))  then
208    ABI_ALLOCATE(fcart_mold,(3,ab_mover%natom))
209  end if
210  if (.not.allocated(fcart_m))     then
211    ABI_ALLOCATE(fcart_m,(3,ab_mover%natom))
212  end if
213 
214 !write(std_out,*) 'langevin 03',jump_end_of_cycle
215 !##########################################################
216 !### 03. Obtain the present values from the history
217 
218  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
219 
220  fcart(:,:)=hist%fcart(:,:,hist%ihist)
221  strten(:) =hist%strten(:,hist%ihist)
222  vel(:,:)  =hist%vel(:,:,hist%ihist)
223  etotal    =hist%etot(hist%ihist)
224 
225  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
226  do ii=1,3
227    rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
228  end do
229  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
230 
231 !Get rid of mean force on whole unit cell, but only if no
232 !generalized constraints are in effect
233 !  call fcart2fred(fcart,fred_corrected,rprimd,ab_mover%natom)
234 !  if(ab_mover%nconeq==0)then
235 !    amass_tot=sum(ab_mover%amass(:))
236 !    do ii=1,3
237 !      if (ii/=3.or.ab_mover%jellslab==0) then
238 !        favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
239 !        fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot
240 !      end if
241 !    end do
242 !  end if
243 
244 !write(std_out,*) 'langevin 04',jump_end_of_cycle
245 !##########################################################
246 !### 04. Compute the next values (Only for the first cycle)
247 
248  if (icycle==1) then
249 
250 !  The temperature is linear between initial and final values
251 !  It is here converted from Kelvin to Hartree (kb_HaK)
252    ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
253 !  write(std_out,*) 'KTEMP=',ktemp
254 !  write(std_out,*) 'MDITEMP=',ab_mover%mdtemp(1)
255 !  write(std_out,*) 'MDFTEMP=',ab_mover%mdtemp(2)
256 !  write(std_out,*) 'DELAYPERM=',ab_mover%delayperm
257 
258 
259 !  %%% LANGEVIN DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
260 
261 !  1. Next values are the present ones
262    acell_next(:)=acell(:)
263    ucvol_next=ucvol
264    rprim_next(:,:)=rprim(:,:)
265    rprimd_next(:,:)=rprimd(:,:)
266 
267    if (zDEBUG) then
268      write (std_out,*) '1. Next values are the present ones'
269      write(std_out,*) 'RPRIMD'
270      do kk=1,3
271        write(std_out,*) rprimd(:,kk)
272      end do
273      write(std_out,*) 'RPRIM'
274      do kk=1,3
275        write(std_out,*) rprim(:,kk)
276      end do
277      write(std_out,*) 'ACELL'
278      write(std_out,*) acell(:)
279    end if
280 
281    if(itime==1)then
282 
283 !    2.   Compute twice the kinetic energy of the system, called v2nose
284      v2nose=0.0_dp
285      do iatom=1,ab_mover%natom
286        do idim=1,3
287          v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
288        end do
289      end do
290      if (zDEBUG) then
291        write (std_out,*) '2.   Compute twice the kinetic energy of the system, called v2nose'
292        write (std_out,*) 'V2NOSE=',v2nose
293      end if
294 
295 !    3.   If there is no kinetic energy, use random numbers
296      if (v2nose<=v2tol) then
297        v2nose=0.0_dp
298        do iatom=1,ab_mover%natom
299          do idim=1,3
300 !          uniformrandom returns a uniform random deviate between 0.0 and 1.0
301 !          if it were always 0 or 1, then the following expression
302 !          would give the requested temperature
303            vel(idim,iatom)=(1.0_dp-2.0_dp*uniformrandom(idum))*sqrt(ktemp/ab_mover%amass(iatom))
304 !          Recompute v2nose
305            v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
306          end do
307        end do
308      end if
309 
310      if (zDEBUG) then
311        write (std_out,*) '3.   If there is no kinetic energy, use random numbers'
312        write (std_out,*) 'VEL'
313        do kk=1,ab_mover%natom
314          write (std_out,*) vel(:,kk)
315        end do
316        write (std_out,*) 'V2NOSE=',v2nose
317      end if
318 
319 
320 !    Now, rescale the velocities to give the proper temperature
321      rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose)
322      vel(:,:)=vel(:,:)*rescale_vel
323 !    Recompute v2nose with the rescaled velocities
324      v2nose=0.0_dp
325      do iatom=1,ab_mover%natom
326        do idim=1,3
327          v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
328        end do
329      end do
330      write(message, '(a)' )&
331 &     ' Rescaling or initializing velocities to initial temperature'
332      call wrtout(ab_out,message,'COLL')
333      call wrtout(std_out,message,'COLL')
334      write(message, '(a,D12.5,a,D12.5)' )&
335 &     ' ---  Scaling factor : ',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
336      call wrtout(ab_out,message,'COLL')
337      call wrtout(std_out,message,'COLL')
338      write(message, '(a,D12.5)' )&
339 &     ' ---  Effective temperature',v2nose/3.0_dp/(kb_HaK*ab_mover%natom)
340      call wrtout(ab_out,message,'COLL')
341      call wrtout(std_out,message,'COLL')
342 !    end if itime==0
343    end if
344 
345 !  This section is devoted to the optional atom permutation (JYR 001114)
346 !  Two input variables are needed
347 !  ab_mover%delayperm : is the interval (in time steps) at which
348 !  atoms are tentatively permuted
349 !  default value could be 0
350 !  ab_mover%signperm  : is the type of bias for the permutation
351 !  +1  to favor alternation of species
352 !  -1  to favor segregation
353 
354 !  Force no permutation at initial step
355    if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2) then
356      if (mod(itime-1,ab_mover%delayperm)==0) then
357 !      Try commutation of atoms.
358        write(message, '(a)')' Attempt of commutation '
359        call wrtout(ab_out,message,'COLL')
360        call wrtout(std_out,message,'COLL')
361 !      Compute a 'permutation potential'
362        do iatom=1,ab_mover%natom
363          pot_perm(iatom)=0.0_dp
364          do iatom1=1,ab_mover%natom
365            if (iatom1.ne.iatom) then
366              distx=xcart(1,iatom)-xcart(1,iatom1)
367              distx=distx-acell(1)*nint(distx/acell(1))
368              disty=xcart(2,iatom)-xcart(2,iatom1)
369              disty=disty-acell(2)*nint(disty/acell(2))
370              distz=xcart(3,iatom)-xcart(3,iatom1)
371              distz=distz-acell(3)*nint(distz/acell(3))
372 !            Here we count each atom below 2 angstr as 1, could be customized
373              dist=sqrt(distx*distx+disty*disty+distz*distz)/3.7807
374              write(std_out,*) iatom,iatom1,dist
375              if (ab_mover%typat(iatom).ne.ab_mover%typat(iatom1)) then
376                mcfac=-1
377              else
378                mcfac=1
379              end if
380              if (dist<1.0_dp)  dist=1.0_dp
381              pot_perm(iatom)=pot_perm(iatom)+mcfac*(ab_mover%signperm)*1.0_dp&
382 &             /exp(log(dist)*6.0_dp)
383            end if
384          end do
385        end do
386        write(std_out,*) ' Perm_pot ',pot_perm(:)
387 !      write(message, '(a,10f12.5)' )' Perm_pot ',&
388 !      &         (pot_perm(iatom1),iatom1=1,ab_mover%natom)
389 !      call wrtout(ab_out,message,'COLL')
390 !      call wrtout(std_out,message,'COLL')
391 
392 !      Find the two atoms, of different types, with the highest perm_pot
393        max_perm(:)=-1.0d9
394        do iatom=1,ab_mover%natom
395          if (pot_perm(iatom) > max_perm(ab_mover%typat(iatom))) then
396            max_perm(ab_mover%typat(iatom))=pot_perm(iatom)
397            imax_perm(ab_mover%typat(iatom))=iatom
398          end if
399        end do
400 
401        if(zDEBUG)then
402 !        write(message, '(a,10f12.5)' )' max_Perm ',&
403 !        &      (max_perm(itypat),itypat=1,ab_mover%ntypat)
404 !        call wrtout(std_out,message,'COLL')
405 !        write(message, '(a,10i12)' )' imax_Perm ',&
406 !        &      (imax_perm(itypat),itypat=1,ab_mover%ntypat)
407 !        call wrtout(std_out,message,'COLL')
408          write(std_out,*) 'NTYPAT',ab_mover%ntypat
409          write(message, '(a,10f12.5)' )' max_Perm ',&
410 &         (max_perm(:))
411          call wrtout(std_out,message,'COLL')
412          write(message, '(a,10i12)' )' imax_Perm ',&
413 &         (imax_perm(:))
414          call wrtout(std_out,message,'COLL')
415        end if
416 
417 !      Loop and keep the 2 largest values
418        if (max_perm(1)>max_perm(2)) then
419          maxp1=max_perm(1)
420          maxp2=max_perm(2)
421          iatom1=imax_perm(1)
422          iatom2=imax_perm(2)
423        else
424          maxp1=max_perm(2)
425          maxp2=max_perm(1)
426          iatom1=imax_perm(2)
427          iatom2=imax_perm(1)
428        end if
429 
430        do itypat=3,ab_mover%ntypat
431          if (max_perm(itypat)>maxp1) then
432            maxp2=maxp1
433            iatom2=iatom1
434            maxp1=max_perm(itypat)
435            iatom1=imax_perm(itypat)
436          else if (max_perm(itypat)>maxp2) then
437            maxp2=max_perm(itypat)
438            iatom2=imax_perm(itypat)
439          end if
440        end do
441        write(message, '(2(a,i5))' )' Will commute atom...',iatom1,'...of type ',&
442 &       ab_mover%typat(iatom1)
443        call wrtout(ab_out,message,'COLL')
444        call wrtout(std_out,message,'COLL')
445        write(message, '(2(a,i5))' )'         with atom...',iatom2,'...of type ',&
446 &       ab_mover%typat(iatom2)
447        call wrtout(ab_out,message,'COLL')
448        call wrtout(std_out,message,'COLL')
449 
450 !      Commute the atoms positions
451        distx=xcart(1,iatom1)
452        disty=xcart(2,iatom1)
453        distz=xcart(3,iatom1)
454        xcart(1,iatom1)=xcart(1,iatom2)
455        xcart(2,iatom1)=xcart(2,iatom2)
456        xcart(3,iatom1)=xcart(3,iatom2)
457        xcart(1,iatom2)=distx
458        xcart(2,iatom2)=disty
459        xcart(3,iatom2)=distz
460 !      Convert back to xred (reduced coordinates)
461        call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
462 
463      end if ! if(mod(itime,ab_mover%delayperm)==0)
464 
465    else
466 !    write(std_out,*) "I will jump to the end of cycle"
467      jump_end_of_cycle=.TRUE.
468 
469    end if ! if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2)
470 !  End of the commutation section
471 
472  end if ! if (icycle==1)
473 
474  if (allocated(imax_perm))   then
475    ABI_DEALLOCATE(imax_perm)
476  end if
477 
478 !write(std_out,*) 'langevin 05',jump_end_of_cycle
479 !##########################################################
480 !### 05. Compute the next values (Only for extra cycles)
481 
482  if (icycle>1) then
483 
484 !  write(std_out,*) "Entering internal cycle 2",icycle
485 
486 !  If the energy computed is higher than the current
487 !  (etotal_temp) we have to discard the changes
488 !  and compute again
489 
490 !  write(std_out,*) ch10
491 !  write(std_out,*) 'EVALUATION FORCES',etotal,hist%etot(abihist_findIndex(hist,-1))
492 !  write(std_out,*) ch10
493 
494 !  This is the worst case (2 evaluations of SCFCV)
495    ihist_prev = abihist_findIndex(hist,-1)
496    if (etotal>hist%etot(ihist_prev).and.icycle==2) then
497 
498 !    Discard the changes
499      acell(:)   =hist%acell(:,ihist_prev)
500      rprimd(:,:)=hist%rprimd(:,:,ihist_prev)
501      xred(:,:)  =hist%xred(:,:,ihist_prev)
502      fcart(:,:) =hist%fcart(:,:,ihist_prev)
503      strten(:)  =hist%strten(:,ihist_prev)
504      vel(:,:)   =hist%vel(:,:,ihist_prev)
505      etotal     =hist%etot(ihist_prev)
506      call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
507 
508 !    distx=xcart(1,iatom1)
509 !    disty=xcart(2,iatom1)
510 !    distz=xcart(3,iatom1)
511 !    xcart(1,iatom1)=xcart(1,iatom2)
512 !    xcart(2,iatom1)=xcart(2,iatom2)
513 !    xcart(3,iatom1)=xcart(3,iatom2)
514 !    xcart(1,iatom2)=distx
515 !    xcart(2,iatom2)=disty
516 !    xcart(3,iatom2)=distz
517 
518 !    Convert back to xred (reduced coordinates)
519 !    call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
520      write(message, '(a)' )' Commutation unsuccessful, recomputing the forces'
521      call wrtout(ab_out,message,'COLL')
522      call wrtout(std_out,message,'COLL')
523 
524 !    This is the best case (only 1 evaluation of SCFCV)
525    else
526 
527      write(message, '(a)')' Commutation successful ! Going on'
528      call wrtout(ab_out,message,'COLL')
529      call wrtout(std_out,message,'COLL')
530 
531 !    Get rid of mean force on whole unit cell, but only if no generalized
532 !    constraints are in effect
533 !    call fcart2fred(fcart,fred_corrected,rprimd,ab_mover%natom)
534 !    if(ab_mover%nconeq==0)then
535 !      amass_tot=sum(ab_mover%amass(:))
536 !      do ii=1,3
537 !        if (ii/=3.or.ab_mover%jellslab==0) then
538 !          favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
539 !          fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot
540 !        end if
541 !      end do
542 !    end if
543 
544 !    In thisc case we do not need to compute again SCFCV
545 !    We avoid the second iteration on ii
546      jump_end_of_cycle=.TRUE.
547 
548    end if ! etotal > etotal_temp
549 
550  end if ! if (icycle>1)
551 
552 !write(std_out,*) 'langevin 06',jump_end_of_cycle
553 !##########################################################
554 !### 06. Compute the next values (Only for the last cycle)
555 
556  if(jump_end_of_cycle)then
557 !  icycle=ncycle
558 !  write(std_out,*) 'This is the last cycle, avoid the others and continue'
559    skipcycle=.TRUE.
560  else
561    skipcycle=.FALSE.
562  end if
563 
564  if ((icycle==ncycle).OR.(skipcycle)) then
565 
566 !  write(std_out,*) 'ENTERING THE FINAL PART',icycle,ncycle
567 
568 !  Specific to Langevin dynamics
569 !  Initialize an array of random forces
570 !  No random force at itime=0
571 !  if (itime==0) then
572    if (itime<0) then
573 
574      ran_force(:,:)=0.0_dp
575 
576    else
577 
578      do iatom=1,ab_mover%natom
579 !      sig_gauss is the std deviation of the random distribution
580        sig_gauss=sqrt(2.0_dp*(ab_mover%friction)*ab_mover%amass(iatom)*ktemp)
581 !      write(std_out,*) 'sig_gauss=',sig_gauss
582 !      write(std_out,*) 'friction=',ab_mover%friction
583 !      write(std_out,*) 'ktemp=',ktemp
584        do idim=1,3
585          delxi=2.0_dp
586          do while (delxi >= 1.0_dp)
587            ran_num1=2.0_dp*uniformrandom(idum)-1.0_dp
588            ran_num2=2.0_dp*uniformrandom(idum)-1.0_dp
589            delxi=ran_num1*ran_num1+ran_num2*ran_num2
590 !          write(std_out,*) delxi,ran_num1,ran_num2
591          end do
592          ran_force(idim,iatom)=ran_num1*sqrt(-2.0_dp*log(delxi)/delxi)&
593 &         *sig_gauss/sqrt(ab_mover%dtion)
594 !        write(std_out,*) 'ran_force',ran_force(idim,iatom)
595        end do
596      end do
597 
598      if (zDEBUG) then
599        write (std_out,*) '4. Different forces computed'
600        write (std_out,*) 'RAN_FORCE'
601        do kk=1,ab_mover%natom
602          write (std_out,*) ran_force(:,kk)
603        end do
604      end if
605 
606 
607      if(zDEBUG)then
608 !      The distribution should be gaussian
609        delxi=0.0_dp
610        do iatom=1,ab_mover%natom
611          do idim=1,3
612            delxi=delxi+(ran_force(idim,iatom)*ab_mover%dtion)**2
613          end do
614        end do
615        delxi=delxi/(3.0_dp*ab_mover%natom)
616        write(message, '(2(a,es22.14))' )' variance =',delxi,'  asked =',&
617 &       2.0_dp*(ab_mover%friction)*ab_mover%amass(2)*ktemp*ab_mover%dtion
618        call wrtout(std_out,message,'COLL')
619      end if
620 !    end if itime\=0
621 
622    end if
623 
624 !  zDEBUG
625 !  write(message, '(a)' )' after initializing ran_force'
626 !  call wrtout(ab_out,message,'COLL')
627 !  call wrtout(std_out,message,'COLL')
628 !  ENDzDEBUG
629 
630    do iatom=1,ab_mover%natom
631      do idim=1,3
632        fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
633        ran_force(idim,iatom)=ran_force(idim,iatom)/ab_mover%amass(iatom)
634      end do
635    end do
636    lang_force(:,:)=ran_force(:,:)-(ab_mover%friction)*vel(:,:)+fcart_m(:,:)
637 
638    if (zDEBUG) then
639      write (std_out,*) '4. Different forces computed'
640      write (std_out,*) 'FCART_M'
641      do kk=1,ab_mover%natom
642        write (std_out,*) fcart_m(:,kk)
643      end do
644      write (std_out,*) 'RAN_FORCE'
645      do kk=1,ab_mover%natom
646        write (std_out,*) ran_force(:,kk)
647      end do
648      write (std_out,*) 'LANG_FORCE'
649      do kk=1,ab_mover%natom
650        write (std_out,*) lang_force(:,kk)
651      end do
652    end if
653 
654 !  zDEBUG
655 !  write(message, '(a)' )'before verlet'
656 !  call wrtout(ab_out,message,'COLL')
657 !  call wrtout(std_out,message,'COLL')
658 !  ENDzDEBUG
659 
660 !  Compute next atomic coordinates using Verlet algorithm
661 
662 !  Impose no change of acell, ucvol, rprim, and rprimd
663    acell_next(:)=acell(:)
664    ucvol_next=ucvol
665    rprim_next(:,:)=rprim(:,:)
666    rprimd_next(:,:)=rprimd(:,:)
667 
668 !  Convert input xred (reduced coordinates) to xcart (cartesian)
669    call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
670 !  Uses the velocity
671 !
672 !  If an atom wants to cross the walls, velocity is reversed.
673 !
674    do iatom=1,ab_mover%natom
675      do idim=1,3
676        delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ &
677 &       0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom)
678        if ( (delxi > (rprimd(idim,idim)+(ab_mover%mdwall)) ) .or. &
679 &       (delxi < - (ab_mover%mdwall)                   )       ) then
680          vel(idim,iatom)=-vel(idim,iatom)
681          delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ &
682 &         0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom)
683        end if
684        xcart_next(idim,iatom)=delxi
685      end do
686    end do
687    xcart(:,:)=xcart_next(:,:)
688    if (zDEBUG) then
689      write (std_out,*) '5. If an atom wants to cross the walls, velocity is reversed.'
690      write (std_out,*) 'XCART'
691      do kk=1,ab_mover%natom
692        write (std_out,*) xcart(:,kk)
693      end do
694    end if
695 
696 !  Convert back to xred_next (reduced coordinates)
697    call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
698    call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
699 
700    if (itime==1) then
701 !    no old forces are available at first step
702 !    Simple update of the velocity
703 !    first compute vel_nexthalf for next steps
704      vel(:,:)=vel(:,:)+ab_mover%dtion*lang_force(:,:)
705    else
706 !    case itime /= 0 normal verlet integration
707      vel(:,:)=vel(:,:)+0.5_dp*ab_mover%dtion*(fcart_mold(:,:)+lang_force(:,:))
708    end if
709    if (zDEBUG) then
710      write (std_out,*) '5. Change velocity with verlet'
711      write (std_out,*) 'VEL'
712      do kk=1,ab_mover%natom
713        write (std_out,*) vel(:,kk)
714      end do
715    end if
716 
717 !  Store 'current force' as 'old force'
718    fcart_mold(:,:)=lang_force(:,:)
719 
720  end if ! if (icycle==ncycle)
721 
722 !write(std_out,*) 'langevin 07',jump_end_of_cycle
723 !##########################################################
724 !### 07. Update the history with the prediction
725 
726 !Increase indexes
727  hist%ihist = abihist_findIndex(hist,+1)
728 
729  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
730  hist%vel(:,:,hist%ihist)=vel(:,:)
731  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
732 
733  if (ab_mover%delayperm==0 .or. ab_mover%ntypat<=2) ncycle=1
734  if(itime==ntime-1) ncycle=1
735 
736 end subroutine pred_langevin