TABLE OF CONTENTS


ABINIT/m_pred_moldyn [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_moldyn

FUNCTION

COPYRIGHT

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

20 #if defined HAVE_CONFIG_H
21 #include "config.h"
22 #endif
23 
24 #include "abi_common.h"
25 
26 module m_pred_moldyn
27 
28  use defs_basis
29  use m_abicore
30  use m_abimover
31  use m_abihist
32 
33  use m_geometry,  only : xcart2xred, xred2xcart
34  use m_predtk,    only : fdtion
35 
36  implicit none
37 
38  private

ABINIT/pred_moldyn [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_moldyn

FUNCTION

 Ionmov predictor (1) Molecular dynamics

 Molecular dynamics, with or without viscous damping
 This function should be after the call to scfcv
 Updates positions, velocities and forces

INPUTS

 ab_mover<type abimover>=Subset of dtset only related with
          |                 movement of ions and acell, contains:
          | dtion:  Time step
          ! natom:  Number of atoms
          | vis:    viscosity
          | iatfix: Index of atoms and directions fixed
          | amass:  Mass of ions
 icycle: Index of the internal cycle inside a time step (itime)
 itime: Index of time iteration
 zDEBUG : if true print some debugging information

OUTPUT

SIDE EFFECTS

 hist<type abihist>=Historical record of positions, forces,
                               stresses, cell and energies,

 ncycle: Number of cycles of a particular time step

NOTES

 * This routine is a predictor, it only produces new positions
   to be computed in the next iteration, this routine should
   produce not output at all
 * ncycle changes from 4 for the first iteration (itime==1) to 1 for (itime>1)
 * The arrays vec_tmp1 and vec_tmp2 are triky, they are use with
   different meanings, during the initialization they contains
   working positions and velocities that acumulated produce the
   first positions of itime=1, for itime>1 they will contain
   positions in 2 previous steps, those values are different
   from the values store in the history, thats the reason why
   we cannot simply use hist%xred to obtain those positions.

PARENTS

      mover

CHILDREN

      hist2var,var2hist,xcart2xred,xred2xcart

SOURCE

 99 subroutine pred_moldyn(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit)
100 
101 
102 !This section has been created automatically by the script Abilint (TD).
103 !Do not modify the following lines by hand.
104 #undef ABI_FUNC
105 #define ABI_FUNC 'pred_moldyn'
106 !End of the abilint section
107 
108 implicit none
109 
110 !Arguments ------------------------------------
111 !scalars
112 type(abimover),intent(in)       :: ab_mover
113 type(abihist),intent(inout),target :: hist
114 integer,intent(in)    :: icycle
115 integer,intent(inout) :: ncycle
116 integer,intent(in)    :: itime
117 integer,intent(in)    :: ntime
118 integer,intent(in)    :: iexit
119 logical,intent(in)    :: zDEBUG
120 
121 !Local variables-------------------------------
122 !scalars
123 integer  :: kk,jj,ihist,ihist_next,ihist_prev,ihist_prev2
124 integer  :: ihist_prev4,ihist_prev5
125 real(dp) :: aa,alfa,bb,cc,x0,xm,em,vis,dx,dv
126 real(dp) :: fcart,fprev,fprev2
127 real(dp) :: xc
128 real(dp) :: vel,vnow,xnow,vprev
129 real(dp),save :: hh,time
130 !arrays
131 real(dp) :: acell(3),rprimd(3,3)
132 real(dp) :: xred(3,ab_mover%natom)
133 real(dp),allocatable :: xcart(:,:),xcart_prev(:,:)
134 real(dp),save,allocatable :: vec_tmp1(:,:)
135 real(dp),save,allocatable :: vec_tmp2(:,:)
136 real(dp), ABI_CONTIGUOUS pointer :: vel_cur(:,:),vel_next(:,:)
137 real(dp),pointer :: fcart_cur(:,:),fcart_prev(:,:),fcart_prev2(:,:)
138 
139 !***************************************************************************
140 !Beginning of executable session
141 !***************************************************************************
142 
143  if(iexit/=0)then
144    if(allocated(vec_tmp1))  then
145      ABI_DEALLOCATE(vec_tmp1)
146    end if
147    if(allocated(vec_tmp2))  then
148      ABI_DEALLOCATE(vec_tmp2)
149    end if
150    return
151  end if
152 
153  vis= ab_mover%vis
154 !Just to avoid warnings of uninitialized variables
155  fprev=0.0_dp
156  fprev=0.0_dp
157  fprev2=0.0_dp
158  vnow=0.0_dp
159  vprev=0.0_dp
160  xnow=0.0_dp
161 
162 !Those arrays contains intermediary results used with
163 !different meanings during the different time steps
164 !We need to preserv the allocation status, this is the
165 !reason to be 'SAVE'
166  if (itime==1.and.icycle==1)then
167    if(allocated(vec_tmp1))  then
168      ABI_DEALLOCATE(vec_tmp1)
169    end if
170    if(allocated(vec_tmp2))  then
171      ABI_DEALLOCATE(vec_tmp2)
172    end if
173    ABI_ALLOCATE(vec_tmp1,(3,ab_mover%natom))
174    ABI_ALLOCATE(vec_tmp2,(3,ab_mover%natom))
175  end if
176 
177 !write(std_out,*) '00'
178 !##########################################################
179 !### 00. Copy from the history to the variables
180 
181  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
182 
183  ABI_ALLOCATE(xcart,(3,ab_mover%natom))
184  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
185 
186  if (itime==1.or.itime==2)then
187    ABI_ALLOCATE(xcart_prev,(3,ab_mover%natom))
188    call xred2xcart(ab_mover%natom,rprimd,xcart_prev,hist%xred(:,:,1))
189  end if
190 
191  ihist = abihist_findIndex(hist, 0)
192  ihist_prev  = abihist_findIndex(hist,-1)
193  ihist_prev2 = abihist_findIndex(hist,-2)
194  ihist_prev4 = abihist_findIndex(hist,-4)
195  ihist_prev5 = abihist_findIndex(hist,-5)
196  ihist_next  = abihist_findIndex(hist,+1)
197 
198  fcart_cur => hist%fcart(:,:,ihist)
199  if (itime==2) fcart_prev  => hist%fcart(:,:,ihist_prev4)
200  if (itime==3) fcart_prev2 => hist%fcart(:,:,ihist_prev5)
201  if (itime >2.or. icycle>=2)fcart_prev  => hist%fcart(:,:,ihist_prev)
202  if (itime >3.or. icycle>=3)fcart_prev2 => hist%fcart(:,:,ihist_prev2)
203 
204  vel_cur  => hist%vel(:,:,ihist)
205  vel_next => hist%vel(:,:,ihist_next)
206 
207 !write(std_out,*) '01'
208 !##########################################################
209 !### 01. Get or compute de time step dtion
210 
211  if (ab_mover%dtion>0)then
212    hh = ab_mover%dtion
213  else
214    hh=fdtion(ab_mover,itime,xcart,fcart_cur,vel_cur)
215  end if
216 
217 !write(std_out,*) '02'
218 !##########################################################
219 !### 02. For all atoms and directions
220  do kk=1,ab_mover%natom
221    em=ab_mover%amass(kk)
222    do jj=1,3
223 
224 !    write(std_out,*) '03'
225 !    ##########################################################
226 !    ### 03. Filling other values from history (forces and vel)
227      fcart=fcart_cur(jj,kk)
228      xc=xcart(jj,kk)
229      vel=hist%vel(jj,kk,1)
230 
231 !    Previous values only after first iteration
232      if (itime>=2.or.icycle>=2) then
233        fprev=fcart_prev(jj,kk)
234        vprev=hist%vel(jj,kk,hist%ihist)
235      end if
236      if (itime>=3.or.icycle>=3) then
237        fprev2=fcart_prev2(jj,kk)
238      end if
239 
240      if (itime==2)then
241        vec_tmp1(jj,kk)=xcart_prev(jj,kk)
242        vec_tmp2(jj,kk)=xcart(jj,kk)
243      end if
244 
245 !    write(std_out,*) '04'
246 !    ##########################################################
247 !    ### 04. Take first the atoms that are not allowed to move along
248 !    ###     this direction
249 !    ###     Warning : implemented in cartesian coordinates
250      if (ab_mover%iatfix(jj,kk)==1) then
251 !      Their positions will be the same as xcart
252        xnow=xcart(jj,kk)
253 !      Their velocities are zero
254        vnow=0.0_dp
255      else
256 
257 !      write(std_out,*) '05'
258 !      ##########################################################
259 !      ### 05. Initialization (itime==1):
260 !      ###     4 calls to obtain the forces are neeeded
261 !      ###     The variables vec_tmp2 and vec_tmp1 from previous
262 !      ###     calls are used in the following ones.
263        if(itime==1)then
264          x0=xcart_prev(jj,kk)
265 
266 !        Prepare the second cycle
267          if(icycle==1)then
268            dx=hh*vel
269            dv=hh/em*(fcart-vis*vel)
270            xnow=x0+.5_dp*dx
271            vnow=vel+.5_dp*dv
272            vec_tmp2(jj,kk)=xc+sixth*dx
273            vec_tmp1(jj,kk)=vel+sixth*dv
274          else if(icycle==2)then
275            dx=hh*vprev
276            dv=hh/em*(fcart-vis*vprev)
277            xnow=x0+.5_dp*dx
278            vnow=vel+.5_dp*dv
279            vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx
280            vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv
281          else if(icycle==3)then
282            dx=hh*vprev
283            dv=hh/em*(fcart-vis*vprev)
284            xnow=x0+dx
285            vnow=vel+dv
286            vec_tmp2(jj,kk)=vec_tmp2(jj,kk)+third*dx
287            vec_tmp1(jj,kk)=vec_tmp1(jj,kk)+third*dv
288          else if(icycle==4)then
289            dx=hh*vprev
290            dv=hh/em*(fcart-vis*vprev)
291            xnow=vec_tmp2(jj,kk)+sixth*dx
292            vnow=vec_tmp1(jj,kk)+sixth*dv
293          end if
294        else !(itime/=1)
295 
296 !        write(std_out,*) '06'
297 !        ##########################################################
298 !        ### 06. Change positions and velocities
299 !        ###     These changes only applies for itime>2
300          if (itime>2)then
301 !          Uses a corrector to have better value of xnow, and
302 !          derive vnow. Only update atoms position and
303 !          velocity along its allowed directions
304            aa=fprev
305            bb=(fcart-fprev2)/(2._dp*hh)
306            cc=(fcart+fprev2-2._dp*fprev)/(2._dp*hh*hh)
307            x0=vec_tmp2(jj,kk)
308            xm=vec_tmp1(jj,kk)
309            if(abs(vis)<=1.d-8)then
310 !            NON-DAMPED DYNAMICS (Post-Code)
311              xnow=2._dp*x0-xm+hh**2/em/12._dp*&
312 &             (fprev2+10._dp*fprev+fcart)
313              vnow=(bb*hh**2)/(3._dp*em)&
314 &             +1.5_dp*aa*hh/em+&
315 &             (5._dp/12._dp)*cc*hh**3/em&
316 &             +x0/hh-xm/hh
317            else
318 !            DAMPED DYNAMICS (Post-Code)
319              alfa=exp(-vis*hh/em)
320              xnow=((-aa*hh*vis**2+0.5_dp*bb*hh**2*vis**2&
321 &             -third*cc*hh**3*vis**2+em*bb*hh*vis&
322 &             -em*cc*hh**2*vis-2._dp*em**2*cc*hh+x0*vis**3-xm*vis**3)*alfa&
323 &             +aa*hh*vis**2-em*bb*hh*vis+third*cc*hh**3*vis**2&
324 &             +2._dp*em**2*cc*hh+0.5D0*bb*hh**2*vis**2-em*cc*hh**2*vis+x0*vis**3)&
325 &             /vis**3
326              vnow=(em*aa*vis**2*alfa-em*aa*vis**2+bb*hh*vis**2*em*alfa&
327 &             -bb*hh*vis**2*em+cc*hh**2*vis**2*em*alfa-cc*hh**2*vis**2*em&
328 &             -em**2*bb*vis*alfa+em**2*bb*vis-2._dp*em**2*cc*hh*vis*alfa+&
329 &             2._dp*em**2*cc*hh*vis+2._dp*em**3*cc*alfa-2._dp*em**3*cc+&
330 &             vis**3*alfa**2*aa*hh-0.5_dp*vis**3*alfa**2*bb*hh**2+&
331 &             third*vis**3*alfa**2*cc*hh**3-vis**2*&
332 &             alfa**2*em*bb*hh+vis**2*alfa**2*em*cc*hh**2+&
333 &             2._dp*vis*alfa**2*em**2*cc*hh-vis**4*alfa**2*x0+&
334 &             vis**4*alfa**2*xm)/vis**3/(alfa-1._dp)/em
335 
336            end if !if(abs(vis)<=1.d-8)
337 
338            xc=xnow
339            vec_tmp1(jj,kk)=vec_tmp2(jj,kk)
340            vec_tmp2(jj,kk)=xnow
341          else
342            vnow=vprev
343          end if !if(itime>2)
344 
345 !        write(std_out,*) '07'
346 !        ##########################################################
347 !        ### 07. Correct positions
348 !        ###     These changes only applies for itime>1
349 
350          if(abs(vis)<=1.d-8)then
351 !          NON-DAMPED DYNAMICS (Pre-Code)
352 !          If the viscosity is too small, the equations become
353 !          ill conditioned due to rounding error so do regular
354 !          Verlet predictor Numerov corrector.
355            x0=vec_tmp2(jj,kk)
356            xm=vec_tmp1(jj,kk)
357            xnow=2._dp*x0-xm&
358 &           + hh**2/em*fcart
359          else
360 !          DAMPED DYNAMICS (Pre-Code)
361 !          These equations come from solving
362 !          m*d2x/dt2+vis*dx/dt=a+b*t+c*t**2
363 !          analytically under the boundary conditions that
364 !          x(0)=x0 and x(-h)=xm, and the following is the
365 !          expression for x(h). a, b and c are determined
366 !          from our knowledge of the driving forces.
367            aa=fcart
368            bb=(fcart-fprev)/hh
369            x0=vec_tmp2(jj,kk)
370            xm=vec_tmp1(jj,kk)
371            alfa=exp(-vis*hh/em)
372            xnow=( (-aa*hh*vis**2 +0.5_dp*bb*hh**2*vis**2&
373 &           +em*bb*hh*vis +x0*vis**3 -xm*vis**3)*alfa&
374 &           +aa*hh*vis**2 -em*bb*hh*vis&
375 &           +0.5_dp*bb*hh**2*vis**2 +x0*vis**3)/vis**3
376 !          End of choice between initialisation, damped
377 !          dynamics and non-damped dynamics
378          end if
379 
380        end if !if(itime==1)
381 
382      end if !if(ab_mover%iatfix(jj,kk)==1)
383 
384 !    write(std_out,*) '08'
385 !    ##########################################################
386 !    ### 08. Update history
387 
388      xcart(jj,kk)=xnow
389      vel_next(jj,kk)=vnow
390 
391 !    write(std_out,*) '09'
392 !    ##########################################################
393 !    ### 09. End loops of atoms and directions
394    end do ! jj=1,3
395  end do ! kk=1,ab_mover%natom
396 
397 !write(std_out,*) '10'
398 !##########################################################
399 !### 10. Filling history with the new values
400 
401  hist%ihist = abihist_findIndex(hist,+1)
402 
403  call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
404  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
405 
406 !Change ncycle for itime>1
407  if (icycle==4)  ncycle=1
408 
409  if (itime==1)then
410    time=0.0_dp
411    if (ab_mover%dtion<0)then
412      write(std_out,*) 'Time=',time
413    end if
414  end if
415  time=time+hh
416  hist%time(hist%ihist)=time
417 
418  if (allocated(xcart)) then
419    ABI_DEALLOCATE(xcart)
420  end if
421  if (allocated(xcart_prev)) then
422    ABI_DEALLOCATE(xcart_prev)
423  end if
424 
425  if (.false.) write(std_out,*) ntime
426 
427 end subroutine pred_moldyn