TABLE OF CONTENTS


ABINIT/m_pred_nose [ Modules ]

[ Top ] [ Modules ]

NAME

  m_pred_nose

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_nose
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 
37  implicit none
38 
39  private

ABINIT/pred_nose [ Functions ]

[ Top ] [ Functions ]

NAME

 pred_nose

FUNCTION

 Ionmov predictors (8) Verlet algorithm with a nose-hoover thermostat

 IONMOV 8:
 Given a starting point xred that is a vector of length 3*natom
 (reduced nuclei coordinates), a velocity vector (in cartesian
 coordinates), and unit cell parameters (acell and rprimd -
 without velocities in the present implementation),
 the Verlet dynamics is performed, using the gradient of the
 energy (atomic forces and stresses) as calculated by the routine scfcv.

 Some atoms can be kept fixed, while the propagation of unit cell
 parameters is only performed if optcell/=0.
 No more than "ntime" steps are performed.
 The time step is governed by dtion (contained in dtset)
 Returned quantities are xred, and eventually acell and rprimd
 (new ones!).

 See ionmov=6, but with a nose-hoover thermostat
 Velocity verlet algorithm : Swope et al JCP 76 (1982) 637

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

SIDE EFFECTS

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

PARENTS

      mover

CHILDREN

      hist2var,metric,var2hist,wrtout,xcart2xred,xred2xcart

SOURCE

 91 subroutine pred_nose(ab_mover,hist,itime,ntime,zDEBUG,iexit)
 92 
 93 
 94 !This section has been created automatically by the script Abilint (TD).
 95 !Do not modify the following lines by hand.
 96 #undef ABI_FUNC
 97 #define ABI_FUNC 'pred_nose'
 98 !End of the abilint section
 99 
100  implicit none
101 
102 !Arguments ------------------------------------
103 !scalars
104  type(abimover),intent(in)       :: ab_mover
105  type(abihist),intent(inout) :: hist
106  integer,intent(in) :: itime
107  integer,intent(in) :: ntime
108  integer,intent(in) :: iexit
109  logical,intent(in) :: zDEBUG
110 
111 !Local variables-------------------------------
112 !scalars
113  integer  :: ii,jj,kk
114  integer  :: idum=-5
115  real(dp),parameter :: v2tol=tol8,nosetol=tol10
116  real(dp) :: delxi,xio,ktemp,rescale_vel
117  real(dp) :: dnose,v2nose,xin_nose
118  real(dp),save :: xi_nose,fsnose,snose
119  real(dp) :: gnose
120  real(dp) :: ucvol,ucvol_next
121  real(dp) :: etotal
122  logical  :: ready
123 
124 !arrays
125  real(dp) :: acell(3),acell_next(3)
126  real(dp) :: rprimd(3,3),rprimd_next(3,3)
127  real(dp) :: gprimd(3,3)
128  real(dp) :: gmet(3,3)
129  real(dp) :: rmet(3,3)
130  real(dp) :: fcart(3,ab_mover%natom)
131 !real(dp) :: fred_corrected(3,ab_mover%natom)
132  real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom)
133  real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom)
134  real(dp) :: vel(3,ab_mover%natom),vel_temp(3,ab_mover%natom)
135  real(dp) :: finose(3,ab_mover%natom),binose(3,ab_mover%natom)
136  real(dp) :: vonose(3,ab_mover%natom),hnose(3,ab_mover%natom)
137  real(dp),allocatable,save :: fcart_m(:,:),fcart_mold(:,:)
138  real(dp) :: strten(6)
139  character(len=500) :: message
140 
141 !***************************************************************************
142 !Beginning of executable session
143 !***************************************************************************
144 
145  if(iexit/=0)then
146    if (allocated(fcart_m))     then
147      ABI_DEALLOCATE(fcart_m)
148    end if
149    if (allocated(fcart_mold))  then
150      ABI_DEALLOCATE(fcart_mold)
151    end if
152    return
153  end if
154 
155 !write(std_out,*) 'nose 01'
156 !##########################################################
157 !### 01. Allocate the arrays fcart_m and fcart_mold
158 
159  if(itime==1)then
160    if (allocated(fcart_m))     then
161      ABI_DEALLOCATE(fcart_m)
162    end if
163    if (allocated(fcart_mold))  then
164      ABI_DEALLOCATE(fcart_mold)
165    end if
166  end if
167 
168  if(.not.allocated(fcart_m))     then
169    ABI_ALLOCATE(fcart_m,(3,ab_mover%natom))
170  end if
171  if(.not.allocated(fcart_mold))  then
172    ABI_ALLOCATE(fcart_mold,(3,ab_mover%natom))
173  end if
174 
175 !write(std_out,*) 'nose 02'
176 !##########################################################
177 !### 02. Obtain the present values from the history
178 
179  call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
180  call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
181 
182  fcart(:,:)=hist%fcart(:,:,hist%ihist)
183  strten(:)=hist%strten(:,hist%ihist)
184  vel(:,:)=hist%vel(:,:,hist%ihist)
185  etotal=hist%etot(hist%ihist)
186 
187  write(std_out,*) 'RPRIMD'
188  do ii=1,3
189    write(std_out,*) rprimd(:,ii)
190  end do
191  write(std_out,*) 'RMET'
192  do ii=1,3
193    write(std_out,*) rmet(ii,:)
194  end do
195 
196 !Get rid of mean force on whole unit cell, but only if no
197 !generalized constraints are in effect
198 !  call fcart2fred(fcart,fred_corrected,rprimd,ab_mover%natom)
199 !  if(ab_mover%nconeq==0)then
200 !    amass_tot=sum(ab_mover%amass(:))
201 !    do ii=1,3
202 !      if (ii/=3.or.ab_mover%jellslab==0) then
203 !        favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom)
204 !        fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot
205 !      end if
206 !    end do
207 !  end if
208 
209 !write(std_out,*) 'nose 03'
210 !##########################################################
211 !### 03. Fill the vectors vin and vout
212 
213 !write(std_out,*) 'nose 04'
214 !##########################################################
215 !### 04. Initialize or update the hessian matrix
216 
217 !write(std_out,*) 'nose 05'
218 !##########################################################
219 !### 05. Compute the next values
220 
221 !The temperature is linear between initial and final values
222 !It is here converted from Kelvin to Hartree (kb_HaK)
223  ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK
224 
225 !%%% NOSE DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
226 
227  acell_next(:)=acell(:)
228  ucvol_next=ucvol
229  rprimd_next(:,:)=rprimd(:,:)
230 
231  if(itime==1)then
232    snose=0.0_dp
233    xi_nose=0.0_dp
234 !  Compute twice the kinetic energy of the system, called v2nose
235    v2nose=0.0_dp
236    do kk=1,ab_mover%natom
237      do jj=1,3
238        v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
239      end do
240    end do
241    if (zDEBUG)then
242      write(std_out,*) 'itime ntime KTEMP=',itime-1,ntime-1,ktemp
243      write(std_out,*) 'V2NOSE=',v2nose
244      write (std_out,*) 'VEL'
245      do kk=1,ab_mover%natom
246        write (std_out,*) vel(:,kk)
247      end do
248    end if
249 
250 !  If there is no kinetic energy, use a random initial velocity
251    if (v2nose<=v2tol) then
252      v2nose=0.0_dp
253      do kk=1,ab_mover%natom
254        do jj=1,3
255 !        Uniform random returns a uniform random deviate between 0.0
256 !        and 1.0
257 !        if it were always 0 or 1, then the following expression
258 !        would give the requested temperature
259          vel(jj,kk)=(1.0_dp-2.0_dp*uniformrandom(idum))*&
260 &         sqrt( (ab_mover%mdtemp(1)) * kb_HaK / ab_mover%amass(kk) )
261 !        Recompute v2nose
262          v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
263          if (zDEBUG)then
264            write(std_out,*) 'jj kk vel(jj,kk)=',jj,kk,vel(jj,kk)
265            write(std_out,*) 'jj kk V2NOSE=',jj,kk,v2nose
266          end if
267        end do
268      end do
269    end if
270    write(std_out,*) 'V2NOSE=',v2nose
271 
272 !  Now, rescale the velocities to give the proper temperature
273    rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose)
274    write(std_out,*) 'RESCALE_VEL=',rescale_vel
275    vel(:,:)=vel(:,:)*rescale_vel
276 !  Recompute v2nose with the rescaled velocities
277    v2nose=0.0_dp
278    do kk=1,ab_mover%natom
279      do jj=1,3
280        v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
281      end do
282    end do
283    write(message, '(a)' )&
284 &   ' Rescaling or initializing velocities to initial temperature'
285    call wrtout(std_out,message,'COLL')
286    call wrtout(std_out,message,'COLL')
287    write(message, '(2(a,es22.14))' )&
288 &   ' ---  Scaling factor : ',rescale_vel,&
289 &   ' Asked T (K) ',ab_mover%mdtemp(1)
290    call wrtout(std_out,message,'COLL')
291    call wrtout(std_out,message,'COLL')
292    write(message, '(a,es22.14)' )&
293 &   ' ---  Effective temperature',v2nose/(3.0_dp*ab_mover%natom*kb_HaK)
294    call wrtout(std_out,message,'COLL')
295    call wrtout(std_out,message,'COLL')
296  end if
297 
298  do kk=1,ab_mover%natom
299    do jj=1,3
300      fcart_m(jj,kk)=fcart(jj,kk)/ab_mover%amass(kk)
301    end do
302  end do
303 
304 !First step of velocity verlet algorithm
305  gnose=3*ab_mover%natom
306 
307 !Convert input xred (reduced coordinates) to xcart (cartesian)
308  call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
309 
310 !Calculate nose-hoover force on atoms
311 !If first iteration, no old force are available, so use present
312 !forces
313  if (itime==1) fcart_mold(:,:)=fcart_m(:,:)
314 
315  if (zDEBUG)then
316    write (std_out,*) 'FCART_MOLD'
317    do kk=1,ab_mover%natom
318      write (std_out,*) fcart_mold(:,kk)
319    end do
320    write (std_out,*) 'FCART_M'
321    do kk=1,ab_mover%natom
322      write (std_out,*) fcart_m(:,kk)
323    end do
324  end if
325 
326  finose(:,:)=fcart_mold(:,:)-xi_nose*vel(:,:)
327  xcart(:,:)=xcart(:,:)+ab_mover%dtion*(vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp)
328 
329 !Convert back to xred (reduced coordinates)
330  call xcart2xred(ab_mover%natom,rprimd,xcart,xred)
331 
332  if (zDEBUG)then
333    write (std_out,*) 'VEL'
334    do kk=1,ab_mover%natom
335      write (std_out,*) vel(:,kk)
336    end do
337  end if
338 
339 !Calculate v2nose
340  v2nose=0.0_dp
341  do kk=1,ab_mover%natom
342    do jj=1,3
343      v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
344    end do
345  end do
346  vel(:,:)=vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp
347 
348  if (zDEBUG)then
349    write(std_out,*) 'NOSE BEFORE'
350    write(std_out,*) 'FSNOSE=',fsnose
351    write(std_out,*) 'SNOSE=',snose
352    write(std_out,*) 'XI_NOSE=',xi_nose
353    write (std_out,*) 'VEL'
354    do kk=1,ab_mover%natom
355      write (std_out,*) vel(:,kk)
356    end do
357    write (std_out,*) 'NOSEINERT',ab_mover%noseinert
358  end if
359 
360 !Update thermostat
361  fsnose=(v2nose-gnose*ktemp)/ab_mover%noseinert
362  snose=snose+ab_mover%dtion*(xi_nose+ab_mover%dtion*fsnose/2.0_dp)
363  xi_nose=xi_nose+ab_mover%dtion*fsnose/2.0_dp
364  if (zDEBUG)then
365    write(std_out,*) 'NOSE AFTER'
366    write(std_out,*) 'FSNOSE=',fsnose
367    write(std_out,*) 'SNOSE=',snose
368    write(std_out,*) 'XI_NOSE=',xi_nose
369    write (std_out,*) 'VEL'
370    do kk=1,ab_mover%natom
371      write (std_out,*) vel(:,kk)
372    end do
373  end if
374 
375 !Second step of the velocity Verlet algorithm, uses the 'new forces'
376 !Calculate v2nose
377  v2nose=0.0_dp
378  do kk=1,ab_mover%natom
379    do jj=1,3
380      v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk)
381    end do
382  end do
383  vel_temp(:,:)=vel(:,:)
384 
385  if (zDEBUG)then
386    write(std_out,*) 'V2NOSE=',v2nose
387    write (std_out,*) 'VEL'
388    do kk=1,ab_mover%natom
389      write (std_out,*) vel(:,kk)
390    end do
391    write (std_out,*) 'Starting Newton Raphson'
392  end if
393 
394  xin_nose=xi_nose
395 
396 !Start Newton-Raphson loop
397  ready=.false.
398  do while (.not.ready)
399    xio=xin_nose
400    delxi=0.0D0
401    vonose(:,:)=vel_temp(:,:)
402    hnose(:,:)=-ab_mover%dtion/2.0_dp*(fcart_m(:,:)-xio*vonose(:,:))-&
403 &   (vel(:,:)-vonose(:,:))
404    do kk=1,ab_mover%natom
405      do jj=1,3
406        binose(jj,kk)=vonose(jj,kk)*ab_mover%dtion/ab_mover%noseinert*ab_mover%amass(kk) ! a verifier
407        delxi=delxi+hnose(jj,kk)*binose(jj,kk)
408      end do
409    end do
410    dnose=-(xio*ab_mover%dtion/2.0D0+1.0D0)
411    delxi=delxi-dnose*((-v2nose+gnose*ktemp)*ab_mover%dtion/2.0_dp/ &
412 &   ab_mover%noseinert-(xi_nose-xio))
413    delxi=delxi/(-ab_mover%dtion*ab_mover%dtion/2.0_dp*v2nose/ab_mover%noseinert+dnose)
414 
415 !  hzeronose=-(xio-xi_nose-(v2nose-gnose*ktemp)
416 !  *dtion/(2.0_dp*ab_mover%noseinert) )
417 !  cibinose=-v2nose*dtion*dtion/(2.0_dp*ab_mover%noseinert)
418 !  delxi=(delxi+hzeronose*dnose)/(dnose+cibinose)
419 
420 !  DEBUG
421 !  write(message, '(a,es22.14)' )' after delxi',delxi
422 !  call wrtout(std_out,message,'COLL')
423 !  call wrtout(std_out,message,'COLL')
424 !  ENDDEBUG
425    v2nose=0.0_dp
426 
427    vel_temp(:,:)=vel_temp(:,:)+&
428 &   (hnose+ab_mover%dtion/2.0_dp*vonose(:,:)*delxi)/dnose
429    do kk=1,ab_mover%natom
430      do jj=1,3
431        v2nose=v2nose+vel_temp(jj,kk)*&
432 &       vel_temp(jj,kk)*ab_mover%amass(kk)
433      end do
434    end do
435 !  New guess for xi
436    xin_nose=xio+delxi
437 
438 !  zDEBUG
439 !  write(message, '(a,es22.14)' )' v2nose=',v2nose
440 !  call wrtout(std_out,message,'COLL')
441 !  call wrtout(std_out,message,'COLL')
442 !  ENDDEBUG
443 
444    ready=.true.
445 !  Test for convergence
446    kk=0
447    jj=1
448    do while((kk<=ab_mover%natom).and.(jj<=3).and.ready)
449      kk=kk+1
450      if (kk>ab_mover%natom) then
451        kk=1
452        jj=jj+1
453      end if
454      if ((kk<=ab_mover%natom) .and.(jj<=3)) then
455        if (abs(vel_temp(jj,kk))<1.0d-50)&
456 &       vel_temp(jj,kk)=1.0d-50
457        if (abs((vel_temp(jj,kk)-vonose(jj,kk))&
458 &       /vel_temp(jj,kk))>nosetol) ready=.false.
459      else
460        if (xin_nose<1.0d-50) xin_nose=1.0d-50
461        if (abs((xin_nose-xio)/xin_nose)>nosetol) ready=.false.
462      end if
463    end do   ! end of while
464 
465 !  Enddo ready
466  end do
467 
468 !Update velocities to converged value
469  vel(:,:)=vel_temp(:,:)
470  write(message, '(a,es14.7)' )' converged velocities for T=',ktemp
471  call wrtout(std_out,message,'COLL')
472 
473  if (zDEBUG)then
474    write (std_out,*) 'Final Values for NOSE'
475    write (std_out,*) 'VEL'
476    do kk=1,ab_mover%natom
477      write (std_out,*) vel(:,kk)
478    end do
479    write (std_out,*) 'XCART'
480    do kk=1,ab_mover%natom
481      write (std_out,*) xcart(:,kk)
482    end do
483  end if
484 
485 !Update thermostat
486  xi_nose=xin_nose
487  xcart_next(:,:)=xcart(:,:)
488 !Convert back to xred_next (reduced coordinates)
489  call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next)
490 !Store 'new force' as 'old force'
491  fcart_mold(:,:)=fcart_m(:,:)
492 
493 !write(std_out,*) 'nose 06'
494 !##########################################################
495 !### 06. Update the history with the prediction
496 
497 !Increase indexes
498  hist%ihist=abihist_findIndex(hist,+1)
499 
500  call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG)
501  hist%vel(:,:,hist%ihist)=vel(:,:)
502  hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion
503 
504 end subroutine pred_nose