TABLE OF CONTENTS


ABINIT/m_pred_lotf [ Modules ]

[ Top ] [ Modules ]

NAME

 m_pred_lotf

FUNCTION

 Contains the predictor for LOTF (ionmov==23)

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 .
 For the initials of contributors,
 see ~abinit/doc/developers/contributors.txt .

SOURCE

18 #if defined HAVE_CONFIG_H
19 #include "config.h"
20 #endif
21 
22 #include "abi_common.h"
23 
24 
25 module m_pred_lotf
26 
27  use m_abicore
28  use defs_basis
29  use m_abimover
30  use m_abihist
31  use lotfpath
32 
33  use m_numeric_tools, only : uniformrandom
34  use m_geometry,      only : xcart2xred
35 
36  implicit none
37 
38  private
39 
40  public :: pred_lotf
41 
42 CONTAINS !===========================================================

ABINIT/m_pred_lotf/pred_lotf [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_lotf

FUNCTION

 Ionmov predictors (12) Lotf ensemble molecular dynamics

 IONMOV 23:
 Lotf ensemble molecular dynamics.

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 .
 For the initials of contributors,
 see ~abinit/doc/developers/contributors.txt .

INPUTS

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

OUTPUT

SIDE EFFECTS

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

NOTES

PARENTS

      mover

CHILDREN

      extrapolation_loop,fitclus,hist2var,init_lotf,intparms
      lotf_interpolation,var2hist,vel_rescale,vel_to_gauss,wrtout,xcart2xred
      xred2xcart

SOURCE

 88  subroutine pred_lotf(ab_mover,hist,itime,icycle,zDEBUG,iexit)
 89 
 90  use m_geometry,       only : xred2xcart
 91 
 92 !This section has been created automatically by the script Abilint (TD).
 93 !Do not modify the following lines by hand.
 94 #undef ABI_FUNC
 95 #define ABI_FUNC 'pred_lotf'
 96 !End of the abilint section
 97 
 98   implicit none
 99 
100   !Arguments ------------------------
101   type(abimover),intent(in)       :: ab_mover
102   type(abihist),intent(inout) :: hist
103   integer,intent(in) :: itime
104   integer,intent(in) :: icycle
105   integer,intent(in) :: iexit
106   logical,intent(in) :: zDEBUG
107 
108   !Local variables-------------------------------
109   !scalars
110   integer  :: kk,iatom,idim,idum=5
111   real(dp) :: v2gauss
112   real(dp),parameter :: v2tol=tol8
113   real(dp) :: etotal
114   character(len=5000) :: message
115   logical :: do_extrap,do_interp,do_first
116   !arrays
117   real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:)
118   real(dp),allocatable,save :: xcart_old(:,:),vel_old(:,:)
119 
120   real(dp) :: acell(3),rprimd(3,3),fcart(3,ab_mover%natom)
121   real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
122   real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
123   real(dp) :: vel(3,ab_mover%natom)
124   real(dp) :: strten(6)
125 
126   !***************************************************************************
127   !Beginning of executable session
128   !***************************************************************************
129 
130   if(iexit/=0)then
131     if (allocated(fcart_m))       then
132       ABI_DEALLOCATE(fcart_m)
133     end if
134     if (allocated(vel_nexthalf))  then
135       ABI_DEALLOCATE(vel_nexthalf)
136     end if
137     if (allocated(xcart_old))       then
138       ABI_DEALLOCATE(xcart_old)
139     end if
140     if (allocated(vel_old))       then
141       ABI_DEALLOCATE(vel_old)
142     end if
143     return
144   end if
145 
146   !write(std_out,*) 'lotf 01'
147   !##########################################################
148   !### 01. Debugging and Verbose
149 
150   if(zDEBUG)then
151     write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),&
152 &   'Debugging and Verbose for pred_lotf',('-',kk=1,37)
153     write(std_out,*) 'ionmov: ',12
154     write(std_out,*) 'itime:  ',itime
155   end if
156 
157   !write(std_out,*) 'lotf 02'
158   !##########################################################
159   !### 02. Allocate the vectors vin, vout and hessian matrix
160   !###     These arrays could be allocated from a previus
161   !###     dataset that exit before itime==ntime
162 
163 
164   if (.not.allocated(fcart_m))       then
165     ABI_ALLOCATE(fcart_m,(3,ab_mover%natom))
166   end if
167   if (.not.allocated(vel_nexthalf))  then
168     ABI_ALLOCATE(vel_nexthalf,(3,ab_mover%natom))
169   end if
170   if (.not.allocated(vel_old))       then
171     ABI_ALLOCATE(vel_old,(3,ab_mover%natom))
172   end if
173   if (.not.allocated(xcart_old))  then
174     ABI_ALLOCATE(xcart_old,(3,ab_mover%natom))
175   end if
176 
177 
178   !write(std_out,*) 'lotf 03'
179   !##########################################################
180   !### 03. Set what Pred_lotf has to do
181 
182   !--First step only the first time pred_lotf is called
183   do_first = (itime==1 .and. icycle==1)
184 
185   !--Extrapolation is computed only is itime is a multiple of lotfvar%nitex
186   do_extrap = lotf_extrapolation(itime)
187 
188   !--Interpolation is done at icycle==2 when extrapolation is active,
189   !--else at the first cycle
190   if(do_extrap) then
191     do_interp = (icycle==2)
192   else
193     do_interp = (icycle==1)
194   endif
195 
196   !write(std_out,*) 'lotf 04'
197   !##########################################################
198   !### 04. Obtain the present values from the history
199 
200   call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
201 
202   fcart(:,:)=hist%fcart(:,:,hist%ihist)
203   strten(:) =hist%strten(:,hist%ihist)
204   vel(:,:)  =hist%vel(:,:,hist%ihist)
205   etotal    =hist%etot(hist%ihist)
206 
207   if(zDEBUG)then
208     write (std_out,*) 'fcart:'
209     do kk=1,ab_mover%natom
210       write (std_out,*) fcart(:,kk)
211     end do
212    write (std_out,*) 'vel:'
213     do kk=1,ab_mover%natom
214       write (std_out,*) vel(:,kk)
215     end do
216     write (std_out,*) 'strten:'
217     write (std_out,*) strten(1:3),ch10,strten(4:6)
218     write (std_out,*) 'etotal:'
219     write (std_out,*) etotal
220   end if
221 
222   !write(std_out,*) 'lotf 05'
223   !##########################################################
224   !### 05. First half-time (First cycle the loop is double)
225 
226   first_step: if(do_first) then
227     write(message, '(a)' ) ' --- LOTF INITIALIZATION---'
228     call wrtout(ab_out,message,'COLL')
229     call wrtout(std_out,message,'COLL')
230 
231     !--convert from xred to xcart
232     call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
233 
234     !PRINT *,'HH-first-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
235 
236     !--call the LOTF initialization
237     call init_lotf(itime,ab_mover%natom,acell,rprimd,xcart)
238 
239     !--Application of Gauss' principle of least constraint according to 
240     ! Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, p.6102 [[cite:Zhang1997]])
241     !--v2gauss is twice the kinetic energy
242     call vel_to_gauss(vel,ab_mover%amass,v2gauss)
243 
244     !--If there is no kinetic energy
245     if(v2gauss<=v2tol) then
246       !--Maxwell-Boltzman distribution
247       v2gauss=zero
248       do iatom=1,ab_mover%natom
249         do idim=1,3
250           vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum))
251           vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum)))
252         end do
253       end do
254     endif
255 
256     !--Init alpha_end parameters
257     call fitclus(.true.,fcart(:,:),xcart,alpha_end,1,ab_mover%natom)
258 
259     !--Forces/mass for the first step
260     forall(iatom = 1:ab_mover%natom) fcart_m(:,iatom) = fcart(:,iatom)/ab_mover%amass(iatom)
261 
262     !print *,'b-end 1step',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
263     !PRINT *,'HH-first-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
264 
265   end if first_step
266 
267 
268   !--Do extrapolation:
269   !  when extrapolation is active for a fixed itime, mover procedure
270   !  will call pred_lotf twice,
271   !  1-a) The final values of alpha are set as initial value
272   !  1-b) The actual coordinates and velocities are stored in a local saved array
273   !  1-c) I the first the coordinated of the atoms at the step itime+lotfvar%nitex
274   !        will be extrapolated and alpha_in updated
275   !  2-a) The new final value of alpha_end (step itime+lotfvar%nitex) are computed
276   !  2-b) Restore velocities and positions
277   extrapolation_lotf: if(do_extrap) then
278     select case(icycle)
279     case(1)
280       !--
281       write(message, '(a)' )' LOTF : Extrapolation'
282       call wrtout(ab_out,message,'COLL')
283       call wrtout(std_out,message,'COLL')
284 
285       !--Intial alpha values to do extrpolation correspond the last one
286       alpha_in = alpha_end
287 
288       !--store the value of xcart
289       xcart_old = xcart
290       vel_old = vel
291 
292     !PRINT *,'HH-extra-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
293 
294       !--Compute new bond, and positions by extrpolation (xcartfit) at
295       !  the step=itime+lotfvar%nitex
296       call extrapolation_loop(itime,ab_mover%mdtemp(1),ab_mover%dtion,ab_mover%amass,xcart,&
297 &                     vel,xcart_next,vel_nexthalf)
298 
299       !print *,'b-end extra',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
300 
301       !--Convert back to xred (reduced coordinates)
302       call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
303 
304     !PRINT *,'HH-extra-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
305 
306     case(2)
307     !PRINT *,'HH-extra-3',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
308 
309       call fitclus(.true.,fcart(:,:),xcart,alpha_end,1,ab_mover%natom)
310       !print *,'b-end fitclus',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
311 
312       !--restore the value of xcart
313       xcart = xcart_old
314       vel = vel_old
315     !PRINT *,'HH-extra-4',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
316 
317     end select
318   endif extrapolation_lotf
319 
320 
321   !--Interpolation(pass here for extrapolation when icycle==2 else when icycle==1)
322   if(do_interp)then
323 
324     !PRINT *,'HH-norma-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
325 
326     !--Compute rescaled vel (center of mass speed and gaussian)
327     call vel_rescale(ab_mover%mdtemp(1),vel,ab_mover%amass,v2gauss)
328 
329     !PRINT *,'HH-norma-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
330 
331     !--VERLETVEL (interpolation) using the interpolated value of alpha
332     !  in the interval [alpha_in,alpha_end]
333     call lotf_interpolation(itime,ab_mover%dtion,v2gauss,ab_mover%amass,xcart,vel,fcart_m,xcart_next,vel_nexthalf)
334 
335     !PRINT *,'HH-norma-3',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
336 
337     !print *,'b-end interpo',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
338 
339     !--Convert back to xred (reduced coordinates)
340     call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
341 
342     !--Interpolate alpha_in,alpha_end to obtain: alpha for the next step
343     call intparms(itime)
344     !print *,'b-end intparm',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
345   endif
346 
347   !write(std_out,*) 'lotf 06'
348   !##########################################################
349   !### 06. Update the history with the prediction
350 
351   xcart = xcart_next
352   xred = xred_next
353 
354   !PRINT *,'HH-final-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:)))
355 
356   !print *,"XCART_final",sum(abs(xcart(1,:))),sum(abs(xcart(2,:))),sum(abs(xcart(3,:)))
357 
358   write(message, '(a,3d20.9)' )&
359 &    ' --- XCART_final',sum(abs(xcart(1,:))),sum(abs(xcart(2,:))),sum(abs(xcart(3,:)))
360   call wrtout(std_out,message,'PERS')
361 
362   !Increase indexes
363   hist%ihist = abihist_findIndex(hist,+1)
364 
365   !Fill the history with the variables
366   !xcart, xred, acell, rprimd
367   call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
368   hist%vel(:,:,hist%ihist) = vel(:,:)
369   hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
370 
371   if(zDEBUG)then
372     write (std_out,*) 'xcart:'
373     do kk=1,ab_mover%natom
374       write (std_out,*) xcart(:,kk)
375     end do
376     write (std_out,*) 'fcart:'
377     do kk=1,ab_mover%natom
378       write (std_out,*) fcart(:,kk)
379     end do
380     write (std_out,*) 'vel:'
381     do kk=1,ab_mover%natom
382       write (std_out,*) vel(:,kk)
383     end do
384     write (std_out,*) 'strten:'
385     write (std_out,*) strten(1:3),ch10,strten(4:6)
386     write (std_out,*) 'etotal:'
387     write (std_out,*) etotal
388   end if
389 
390  end subroutine pred_lotf