TABLE OF CONTENTS


ABINIT/m_pred_isokinetic [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_isokinetic

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_isokinetic
28 
29  use m_abicore
30  use defs_basis
31  use m_abimover
32  use m_abihist
33 
34  use m_numeric_tools,  only : uniformrandom
35  use m_geometry,  only : xcart2xred, xred2xcart
36 
37 
38  implicit none
39 
40  private

ABINIT/pred_isokinetic [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_isokinetic

FUNCTION

 Ionmov predictors (12) Isokinetic ensemble molecular dynamics

 IONMOV 12:
 Isokinetic ensemble molecular dynamics.
 The equation of motion of the ions in contact with a thermostat
 are solved with the algorithm proposed by Zhang [J. Chem. Phys. 106, 6102 (1997)] [[cite:Zhang1997]],
 as worked out by Minary et al, J. Chem. Phys. 188, 2510 (2003) [[cite:Minary2003]].
 The conservation of the kinetic energy is obtained within machine precision, at each step.
 Related parameters : the time step (dtion), the initial temperature (mdtemp(1)) if the velocities are not defined to start with.

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
 zDEBUG : if true print some debugging information

OUTPUT

SIDE EFFECTS

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

PARENTS

      mover

CHILDREN

      hist2var,var2hist,wrtout,xcart2xred,xred2xcart

SOURCE

 83 subroutine pred_isokinetic(ab_mover,hist,itime,ntime,zDEBUG,iexit)
 84 
 85 
 86 !This section has been created automatically by the script Abilint (TD).
 87 !Do not modify the following lines by hand.
 88 #undef ABI_FUNC
 89 #define ABI_FUNC 'pred_isokinetic'
 90 !End of the abilint section
 91 
 92  implicit none
 93 
 94 !Arguments ------------------------------------
 95 !scalars
 96  integer,intent(in) :: itime
 97  integer,intent(in) :: ntime
 98  integer,intent(in) :: iexit
 99  logical,intent(in) :: zDEBUG
100  type(abimover),intent(in)       :: ab_mover
101  type(abihist),intent(inout) :: hist
102 
103 !Local variables-------------------------------
104 !scalars
105  integer  :: kk,iatom,idim,idum=5,nxyzatfree,ndegfreedom,nfirst,ifirst
106  real(dp) :: a,as,b,sqb,s,s1,s2,scdot,sigma2,vtest,v2gauss
107  real(dp),parameter :: v2tol=tol8
108  real(dp) :: etotal,rescale_vel
109  character(len=5000) :: message
110 !arrays
111  real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:)
112 
113  real(dp) :: acell(3),rprimd(3,3)
114  real(dp) :: fcart(3,ab_mover%natom)
115 !real(dp) :: fred_corrected(3,ab_mover%natom)
116  real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
117  real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
118  real(dp) :: vel(3,ab_mover%natom)
119  real(dp) :: strten(6)
120 
121 !***************************************************************************
122 !Beginning of executable session
123 !***************************************************************************
124 
125 !DEBUG
126 !write(std_out,*)' pred_isokinetic : enter '
127 !stop
128 !ENDDEBUG
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    return
138  end if
139 
140 !write(std_out,*) 'isokinetic 01'
141 !##########################################################
142 !### 01. Debugging and Verbose
143 
144  if(zDEBUG)then
145    write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),&
146 &   'Debugging and Verbose for pred_isokinetic',('-',kk=1,37)
147    write(std_out,*) 'ionmov: ',12
148    write(std_out,*) 'itime:  ',itime
149  end if
150 
151 !write(std_out,*) 'isokinetic 02'
152 !##########################################################
153 !### 02. Allocate the vectors vin, vout and hessian matrix
154 !###     These arrays could be allocated from a previous
155 !###     dataset that exit before itime==ntime
156 
157  if(itime==1)then
158    if (allocated(fcart_m))       then
159      ABI_DEALLOCATE(fcart_m)
160    end if
161    if (allocated(vel_nexthalf))  then
162      ABI_DEALLOCATE(vel_nexthalf)
163    end if
164  end if
165 
166  if (.not.allocated(fcart_m))       then
167    ABI_ALLOCATE(fcart_m,(3,ab_mover%natom))
168  end if
169  if (.not.allocated(vel_nexthalf))  then
170    ABI_ALLOCATE(vel_nexthalf,(3,ab_mover%natom))
171  end if
172 
173 !write(std_out,*) 'isokinetic 03'
174 !##########################################################
175 !### 03. Obtain the present values from the history
176 
177  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
178 
179  fcart(:,:)=hist%fcart(:,:,hist%ihist)
180  strten(:) =hist%strten(:,hist%ihist)
181  vel(:,:)  =hist%vel(:,:,hist%ihist)
182  etotal    =hist%etot(hist%ihist)
183 
184  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
185 
186  if(zDEBUG)then
187    write (std_out,*) 'fcart:'
188    do kk=1,ab_mover%natom
189      write (std_out,*) fcart(:,kk)
190    end do
191    write (std_out,*) 'vel:'
192    do kk=1,ab_mover%natom
193      write (std_out,*) vel(:,kk)
194    end do
195    write (std_out,*) 'strten:'
196    write (std_out,*) strten(1:3),ch10,strten(4:6)
197    write (std_out,*) 'etotal:'
198    write (std_out,*) etotal
199  end if
200 
201 !Get rid of mean force on whole unit cell, but only if no
202 !generalized constraints are in effect
203 !  call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom)
204 !  if(ab_mover%nconeq==0)then
205 !    amass_tot=sum(ab_mover%amass(:))
206 !    do ii=1,3
207 !      if (ii/=3.or.ab_mover%jellslab==0) then
208 !        favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
209 !        fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot
210 !      end if
211 !    end do
212 !  end if
213 
214 !Count the number of degrees of freedom, taking into account iatfix.
215 !Also fix the velocity to zero for the fixed atoms
216  nxyzatfree=0
217  do iatom=1,ab_mover%natom
218    do idim=1,3
219      if(ab_mover%iatfix(idim,iatom)==0)then
220        nxyzatfree=nxyzatfree+1
221      else
222        vel(idim,iatom)=zero
223      endif
224    enddo
225  enddo
226 
227 !Now, the number of degrees of freedom is reduced by four because of the kinetic energy conservation
228 !and because of the conservation of the total momentum for each dimension, in case no atom position is fixed for that dimension 
229 !(in the latter case, one degree of freedom has already been taken away)
230 !This was not done until v8.9 of ABINIT ...
231  ndegfreedom=nxyzatfree
232  ndegfreedom=nxyzatfree-1 ! Kinetic energy conservation
233  do idim=1,3
234    if(sum(ab_mover%iatfix(idim,:))==0)then
235      ndegfreedom=ndegfreedom-1 ! Macroscopic momentum
236    endif
237  enddo
238 
239 !write(std_out,*) 'isokinetic 04'
240 !##########################################################
241 !### 04. Second half-velocity (Only after the first itime)
242 
243  if(itime>1) then
244 
245    do iatom=1,ab_mover%natom
246      do idim=1,3
247        if(ab_mover%iatfix(idim,iatom)==0)then
248          fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
249        else
250          fcart_m(idim,iatom)=zero
251        endif
252      end do
253    end do
254 
255 !  Computation of vel(:,:) at the next positions
256 !  Computation of v2gauss, actually twice the kinetic energy. 
257 !  Called 2K, cf Eq. (A13) of [[cite:Minary2003]].
258    v2gauss=0.0_dp
259    do iatom=1,ab_mover%natom
260      do idim=1,3
261        v2gauss=v2gauss+&
262 &       vel_nexthalf(idim,iatom)*vel_nexthalf(idim,iatom)*&
263 &       ab_mover%amass(iatom)
264      end do
265    end do
266 
267 !  Computation of a and b (4.13 of [[cite:Minary2003]])
268    a=0.0_dp
269    b=0.0_dp
270    do iatom=1,ab_mover%natom
271      do idim=1,3
272        a=a+fcart_m(idim,iatom)*vel_nexthalf(idim,iatom)*ab_mover%amass(iatom)
273        b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom)
274      end do
275    end do
276    a=a/v2gauss
277    b=b/v2gauss
278 
279 
280 !  Computation of s and scdot
281    sqb=sqrt(b)
282    as=sqb*ab_mover%dtion/2.
283 ! jmb
284    if ( as > 300.0 ) as=300.0
285    s1=cosh(as)
286    s2=sinh(as)
287    s=a*(s1-1.)/b+s2/sqb
288    scdot=a*s2/sqb+s1
289 
290    do iatom=1,ab_mover%natom
291      do idim=1,3
292        if(ab_mover%iatfix(idim,iatom)==0)then
293          vel(idim,iatom)=(vel_nexthalf(idim,iatom)+fcart_m(idim,iatom)*s)/scdot
294        else
295          vel(idim,iatom)=zero
296        endif
297      enddo
298    enddo  
299 
300    if (zDEBUG)then
301      write(std_out,*) 'Computation of the second half-velocity'
302      write(std_out,*) 'Cartesian forces per atomic mass (fcart_m):'
303      do kk=1,ab_mover%natom
304        write (std_out,*) fcart_m(:,kk)
305      end do
306      write(std_out,*) 'vel:'
307      do kk=1,ab_mover%natom
308        write (std_out,*) vel(:,kk)
309      end do
310      write(std_out,*) 'v2gauss:',v2gauss
311      write(std_out,*) 'a:',a
312      write(std_out,*) 'b:',b
313      write(std_out,*) 's:',s
314      write(std_out,*) 'scdot:',scdot
315    end if
316 
317  end if ! (if itime>1)
318 
319 !write(std_out,*) 'isokinetic 05'
320 !##########################################################
321 !### 05. First half-time (First cycle the loop is double)
322 
323  if (itime==1) then
324    nfirst=2
325  else
326    nfirst=1
327  end if
328 
329  do ifirst=1,nfirst
330 
331 !  Application of Gauss' principle of least constraint according to Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, [[cite:Zhang1997]] p.6102)
332 
333 !  v2gauss is twice the kinetic energy
334    v2gauss=0.0_dp
335    do iatom=1,ab_mover%natom
336      do idim=1,3
337        v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
338      end do
339    end do
340 
341 !  If there is no kinetic energy to start with ...
342    if (v2gauss<=v2tol.and.itime==1) then
343 !    Maxwell-Boltzman distribution
344      v2gauss=zero
345      vtest=zero
346      do iatom=1,ab_mover%natom
347        do idim=1,3
348          if(ab_mover%iatfix(idim,iatom)==0)then
349            vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum))
350            vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum)))
351          else
352            vel(idim,iatom)=zero
353          endif
354        end do
355      end do
356 
357 !    Get rid of center-of-mass velocity
358      s1=sum(ab_mover%amass(:))
359      do idim=1,3
360        if(sum(ab_mover%iatfix(idim,:))==0)then
361          s2=sum(ab_mover%amass(:)*vel(idim,:))
362          vel(idim,:)=vel(idim,:)-s2/s1
363        endif
364      end do
365 
366 !    Recompute v2gauss
367      do iatom=1,ab_mover%natom
368        do idim=1,3
369          v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
370          vtest=vtest+vel(idim,iatom)/ndegfreedom
371        end do
372      end do
373 
374 !    Now rescale the velocities to give the exact temperature
375      rescale_vel=sqrt(ndegfreedom*kb_HaK*ab_mover%mdtemp(1)/v2gauss)
376      vel(:,:)=vel(:,:)*rescale_vel
377 
378 !    Recompute v2gauss with the rescaled velocities
379      v2gauss=zero
380      do iatom=1,ab_mover%natom
381        do idim=1,3
382          v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
383        end do
384      end do
385 
386 !    Compute the variance and print
387      sigma2=(v2gauss/ndegfreedom-ab_mover%amass(1)*vtest**2)/kb_HaK
388 
389    end if
390 
391    do iatom=1,ab_mover%natom
392      do idim=1,3
393        if(ab_mover%iatfix(idim,iatom)==0)then
394          fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom)
395        else
396          fcart_m(idim,iatom)=zero
397        endif
398      end do
399    end do
400 
401    if (zDEBUG)then
402      write(std_out,*) 'Calculation first half-velocity '
403      write (std_out,*) 'vel:'
404      do kk=1,ab_mover%natom
405        write (std_out,*) vel(:,kk)
406      end do
407      write (std_out,*) 'xcart:'
408      do kk=1,ab_mover%natom
409        write (std_out,*) xcart(:,kk)
410      end do
411      write (std_out,*) 'xred:'
412      do kk=1,ab_mover%natom
413        write (std_out,*) xred(:,kk)
414      end do
415      write (std_out,*) 'fcart_m'
416      do kk=1,ab_mover%natom
417        write (std_out,*) fcart_m(:,kk)
418      end do
419      write(std_out,*) 's2',s2
420      write(std_out,*) 'v2gauss',v2gauss
421      write(std_out,*) 'sigma2',sigma2
422 
423      write(message, '(a)' )&
424 &     ' --- Rescaling or initializing velocities to initial temperature'
425      call wrtout(std_out,message,'COLL')
426      write(message, '(a,d12.5,a,D12.5)' )&
427 &     ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1)
428      call wrtout(std_out,message,'COLL')
429      write(message, '(a,d12.5,a,D12.5)' )&
430 &     ' --- Effective temperature',v2gauss/(ndegfreedom*kb_HaK),' From variance', sigma2
431      call wrtout(std_out,message,'COLL')
432    end if
433 
434 !  Convert input xred (reduced coordinates) to xcart (cartesian)
435    call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
436 
437    if(itime==1.and.ifirst==1) then
438      call wrtout(std_out,'if itime==1','COLL')
439      vel_nexthalf(:,:)=vel(:,:)
440      xcart_next(:,:)=xcart(:,:)
441      call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
442      xred=xred_next
443      call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
444    end if
445 
446  end do
447 
448 !Computation of vel_nexthalf (4.16 of [[cite:Minary2003]])
449 !Computation of a and b (4.13 of [[cite:Minary2003]])
450  a=0.0_dp
451  b=0.0_dp
452  do iatom=1,ab_mover%natom
453    do idim=1,3
454      a=a+fcart_m(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom)
455      b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom)
456    end do
457  end do
458  a=a/v2gauss+tol20
459  b=b/v2gauss+tol20
460 !Computation of s and scdot
461  sqb=sqrt(b)+tol20
462  as=sqb*ab_mover%dtion/2.
463 ! jmb
464  if ( as > 300.0 ) as=300.0
465  s1=cosh(as)
466  s2=sinh(as)
467  s=a*(s1-1.)/b+s2/sqb
468  scdot=a*s2/sqb+s1
469  do iatom=1,ab_mover%natom
470    do idim=1,3
471      if(ab_mover%iatfix(idim,iatom)==0)then
472        vel_nexthalf(idim,iatom)=(vel(idim,iatom)+fcart_m(idim,iatom)*s)/scdot
473      else
474        vel_nexthalf(idim,iatom)=zero
475      endif
476    enddo
477  enddo
478 
479 !Computation of the next positions
480  xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion
481 
482  if (zDEBUG)then
483    write(std_out,*) 'a:',a
484    write(std_out,*) 'b:',b
485    write(std_out,*) 's:',s
486    write(std_out,*) 'scdot:',scdot
487  end if
488 
489 !Convert back to xred (reduced coordinates)
490 
491  call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
492 
493 !write(std_out,*) 'isokinetic 06'
494 !##########################################################
495 !### 06. Update the history with the prediction
496 
497  xcart=xcart_next
498  xred=xred_next
499 
500 !increment the ihist
501  hist%ihist = abihist_findIndex(hist,+1)
502 
503 !Fill the history with the variables
504 !xred, acell, rprimd, vel
505  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
506  hist%vel(:,:,hist%ihist)=vel(:,:)
507  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
508 
509  if(zDEBUG)then
510    write (std_out,*) 'fcart:'
511    do kk=1,ab_mover%natom
512      write (std_out,*) fcart(:,kk)
513    end do
514    write (std_out,*) 'vel:'
515    do kk=1,ab_mover%natom
516      write (std_out,*) vel(:,kk)
517    end do
518    write (std_out,*) 'strten:'
519    write (std_out,*) strten(1:3),ch10,strten(4:6)
520    write (std_out,*) 'etotal:'
521    write (std_out,*) etotal
522  end if
523 
524  if (.false.) write(std_out,*) ntime
525 
526 end subroutine pred_isokinetic