TABLE OF CONTENTS


ABINIT/lotfpath [ Modules ]

[ Top ] [ Modules ]

NAME

 lotfpath

FUNCTION

COPYRIGHT

 Copyright (C) 2005-2018 ABINIT group (MMancini)
 This file is distributed under the terms of the
 GNU General Public License, see ~abinit/COPYING
 or http://www.gnu.org/copyleft/gpl.txt .

SOURCE

16 #if defined HAVE_CONFIG_H
17 #include "config.h"
18 #endif
19 
20 #include "abi_common.h"
21 
22 module LOTFPATH
23 
24  use defs_basis
25  use m_errors
26  use m_abicore
27 
28  use defs_param_lotf, only : lotfvar
29 
30  implicit none
31  public
32 
33  !--Lotf variables used in non-lotf procedure (ex.moldyn)
34  real(dp),allocatable,dimension(:,:) :: alpha,alpha_in,alpha_end
35 
36  real(dp),private :: epotlotf
37  integer,private,allocatable,dimension(:)   :: nneig,nneig_old
38  integer,private,allocatable,dimension(:,:) :: neighl,neighl_old
39  real(dp),private,allocatable,dimension(:,:) :: xcartfit,velfit,fcartfit
40  character(len=500),private :: message
41 
42  public  :: init_lotf
43  public  :: end_lotf
44  public  :: fitclus
45  public  :: upd_lis
46  public  :: force0
47  public  :: intparms
48  public  :: lotf_extrapolation
49  public  :: vel_to_gauss
50  public  :: force_to_vel
51  public  :: vel_rescale
52  public  :: extrapolation_loop
53 
54  private :: alpha_update
55 
56 contains

lotfpath/init_lotf [ Functions ]

[ Top ] [ lotfpath ] [ Functions ]

NAME

 init_lotf

FUNCTION

INPUTS

PARENTS

      m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

 75  subroutine init_lotf(itime,natom,acell,rprimd,xcart)
 76   ! natom : number of atoms, abinit notation
 77   ! acell : cell length, abinit notation
 78 
 79   use glue_lotf,only  : glue_init
 80   use pbc_lotf,only   : pbc_init
 81   use eval_lotf,only  : upd_lis0
 82   use work_var_lotf,only : work_var_set,smallfit,cutoff_init
 83 
 84   use bond_lotf,only : ibn_tot,ibn_tots,nbondex,ibn_tot2,nfitmax,&
 85 &                  bond_fit_set,bond_tafit_init,&
 86 &                  bond_atom_init,     &
 87 &                  bond_matrix_alloc,  &
 88 &                  bond_matrix_set
 89 
 90 !This section has been created automatically by the script Abilint (TD).
 91 !Do not modify the following lines by hand.
 92 #undef ABI_FUNC
 93 #define ABI_FUNC 'init_lotf'
 94 !End of the abilint section
 95 
 96   integer,intent(in) :: itime
 97   integer,intent(in) :: natom
 98   real(dp),intent(in) :: acell(3),rprimd(3,3)
 99   real(dp),intent(out) :: xcart(3,natom)
100 
101 
102   integer :: nfitdum
103 
104 ! *************************************************************************
105 
106   ! I should modify units : LOTF uses atomic units
107   !                         ABINIT uses  ???
108 
109   ! !-----------------------------------------
110   ! Transfered in gstate
111   ! call lotfvar_init(natom,&
112   !   &               2, & !--version: set type of MD algo
113   !   &           itime, & !--nstart: initial step
114   !   &           nitex,& !--nitex: number of LOTF steps
115   !   &               40,& !--nneigx: roughly decide the number of neighbours
116   !   &               5, & !--classic: stick with the adaptable Glue model (rough version)
117   !   &               1,1& !--me,nproc : disabled parallel LOTF
118   !   )
119 
120   !-----------------------------------------
121   !--Prepare LOTF variables used in moldyn.f90
122 
123   !--initialize potential energy :
124    epotlotf = zero
125 
126    ABI_ALLOCATE(fcartfit,(3,natom))
127    ABI_ALLOCATE(xcartfit,(3,natom))
128    ABI_ALLOCATE(velfit,(3,natom))
129 
130    ABI_ALLOCATE(neighl,(lotfvar%nneigx,natom))
131    ABI_ALLOCATE(nneig,(lotfvar%nneigx))
132    ABI_ALLOCATE(neighl_old,(lotfvar%nneigx,natom))
133    ABI_ALLOCATE(nneig_old,(lotfvar%nneigx))
134    neighl = 0
135    nneig = 0
136    neighl_old = 0
137    nneig_old = 0
138 
139   !--set work variables of LOTF
140    call work_var_set()
141 
142   !--control on lotfvar%classic
143    if(lotfvar%classic/=5 .and. lotfvar%classic/=6) then
144      write(message,'(3a,3f12.6,a)')&
145 &     'LOTF: INIT_LIST: wrong value for lotfvar%classic = ',&
146 &     lotfvar%classic,ch10,&
147 &     'change lotfvar%classic 5 or 6 '
148      MSG_ERROR(message)
149    end if
150 
151   !--Init cell and pbc_lotf
152    call pbc_init(rprimd)
153 
154   !--Initializes Glue, here we suppose we have atomic units.
155    call Glue_INIT()
156 
157   !--Initializes cut-off radius
158    call cutoff_init()
159 
160   !--last argument is to force the search for neighbors in upd_lis0
161    call upd_lis0(xcart,neighl,nneig,itime)
162 
163   ! SMALL_FIT FINDS A RESTRICTED REGION WHERE THE FIT WILL BE
164   ! PERFOMED. dimensionS WILL BE SET ACCORDING THE NUMBER OF ATOMS
165   ! IN THIS SMALL REGION
166   ! if niter=lotfvar%n0 it will set some dimensions for later use
167 
168   !--set tafit
169    call bond_tafit_init(lotfvar%natom)
170    call SMALLFIT(xcart,nfitdum)
171 
172   !--set nfitmax,ifit,nfit
173    call bond_fit_set(lotfvar%natom,nfitdum)
174 
175   !--Initialize nbondex,neighl,neeig
176    call bond_atom_init(lotfvar%nneigx,nneig,neighl)
177 
178   !---------------------------------------------
179   !--Initialize/creates arrays needed by updlis :
180 
181   !--Allocate bond matrix
182    call bond_matrix_alloc(lotfvar%natom,lotfvar%nneigx)
183 
184   !--updates bond matrix, associates the bond to the atom neighlists
185    call bond_matrix_set(nneig,neighl)
186 
187   !--initialize bond parms alpha
188    ABI_ALLOCATE(alpha,(3,nbondex))
189    ABI_ALLOCATE(alpha_in,(3,nbondex))
190    ABI_ALLOCATE(alpha_end,(3,nbondex))
191    alpha = zero
192    alpha_in = zero
193    alpha_end = zero
194 
195    write(message,'(2a,i6,a,i8,2a,i8,2a,i8,i8)')ch10,&
196 &   'ITERATION:  ',itime,&
197 &   ' NO. OF BONDS BETWEEN FITTED ATOMS: ', ibn_tot,ch10,&
198 &   'TOTAL N.OF BOUNDS : ', ibn_tots,ch10,&
199 &   'BORDER BOUNDS (&max) : ', ibn_tot2,6*nfitmax
200    call wrtout(std_out,message,'COLL')
201 
202  end subroutine init_lotf

lothpath/alpha_update [ Functions ]

[ Top ] [ Functions ]

NAME

 alpha_update

FUNCTION

  updates parameter list

INPUTS

  dphi=parameter to reinatialise bond parameters
  jbo= index array for reordering

PARENTS

      m_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1160  subroutine alpha_update(dphi,jbo,alpha_dum)
1161   use bond_lotf,only : ibn_tot
1162 
1163 !This section has been created automatically by the script Abilint (TD).
1164 !Do not modify the following lines by hand.
1165 #undef ABI_FUNC
1166 #define ABI_FUNC 'alpha_update'
1167 !End of the abilint section
1168 
1169   implicit none
1170 
1171   !Arguments ------------------------
1172   real(dp),intent(in) :: dphi
1173   integer,intent(in) :: jbo(:)
1174   real(dp),intent(inout) :: alpha_dum(:,:)
1175   !Local ---------------------------
1176   integer :: ibn,jb_old
1177   real(dp) :: alphawork(3,ibn_tot)
1178 
1179 ! *************************************************************************
1180 
1181    do ibn = 1,ibn_tot
1182      jb_old = jbo(ibn)
1183      if(jb_old /= 0) then
1184       !--swaps old parms
1185        alphawork(:,ibn) = alpha_dum(:,jb_old)
1186      else
1187       !--or reinitialise bond parms
1188        alphawork(:,ibn) = (/ dphi, one, zero /)
1189      end if
1190    end do
1191 
1192    alpha_dum(:,:ibn_tot) = alphawork(:,:)
1193 
1194  end subroutine alpha_update

lothpath/end_lotf [ Functions ]

[ Top ] [ Functions ]

NAME

 end_lotf

FUNCTION

INPUTS

PARENTS

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

220  subroutine end_LOTF()
221 
222   use work_var_lotf
223   use bond_lotf
224 
225   !--init_lotf
226 ! *************************************************************************
227 
228 !This section has been created automatically by the script Abilint (TD).
229 !Do not modify the following lines by hand.
230 #undef ABI_FUNC
231 #define ABI_FUNC 'end_LOTF'
232 !End of the abilint section
233 
234    ABI_DEALLOCATE(alpha)
235    ABI_DEALLOCATE(alpha_in)
236    ABI_DEALLOCATE(alpha_end)
237    ABI_DEALLOCATE(nneig)
238    ABI_DEALLOCATE(neighl)
239    ABI_DEALLOCATE(fcartfit)
240    ABI_DEALLOCATE(xcartfit)
241    ABI_DEALLOCATE(velfit)
242 
243   !--deallocate LOTF internal variables
244    call work_var_dealloc()
245   !--deallocate LOTF bond variables
246    call bond_dealloc()
247 
248  end subroutine end_LOTF

lothpath/extrapolation_loop [ Functions ]

[ Top ] [ Functions ]

NAME

 extrapolation_loop

FUNCTION

  Compute the LOTF extrapolation:
  Starting from xcart_0,vel_0,fcart_0, it computes in first the new
  bond, then the forces on the atoms.
  At this point if compute the position,speed and forces on atoms
  upto the step nitex:  xcart_nitex, vel_nitex, fcart_nitex
  This will be used in a SCF to compute the forces in the final
  point of the LOTF approximation.
  In output the positions (extrapoled) in the step ntitex.

INPUTS

  itime=numeber of MD step
  mditemp=temperature of ions
  dtion=time step for Molecular Dynamics
  amass(natom)=masse of the ions
  xcart(3,natom)=position of the ions initial
  vel(3,natom)=speed of the ions initial

 OUT
  xcart_next(3,natom)=positions of the ions final step
  the following variable are used as work variables:
  vel_nexthalf(3,natom)=velocity of ions in next half step

PARENTS

      m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1489  subroutine extrapolation_loop(itime,mditemp,dtion,amass,&
1490 &                           xcart,vel,xcart_next,vel_nexthalf)
1491 
1492 
1493 !This section has been created automatically by the script Abilint (TD).
1494 !Do not modify the following lines by hand.
1495 #undef ABI_FUNC
1496 #define ABI_FUNC 'extrapolation_loop'
1497 !End of the abilint section
1498 
1499   implicit none
1500 
1501   !Arguments ------------------------------------
1502   integer,intent(in) :: itime
1503   real(dp),intent(in) :: mditemp,dtion
1504   real(dp),intent(in) :: amass(:)
1505   real(dp),intent(in) :: vel(:,:)
1506   real(dp),intent(in) :: xcart(:,:)
1507   real(dp),intent(out) :: vel_nexthalf(:,:)
1508   real(dp),intent(out) :: xcart_next(:,:)
1509   !Local ---------------------------
1510   integer :: itex
1511   real(dp) :: v2gauss
1512 
1513 ! *************************************************************************
1514 
1515   !--inital values for position and forces
1516    xcartfit(:,:) = xcart(:,:)
1517    velfit(:,:) = vel(:,:)
1518 
1519   !--Update bond variables and alpha
1520    call upd_lis(xcartfit,itime,alpha_in)
1521 
1522   !print *,'b-end uplis',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
1523 
1524   !--Store alpha in alpha_in
1525    alpha(:,:) = alpha_in(:,:)
1526 
1527   !--compute fcartfit (fcartfit is the accelleration)
1528    call force0(xcartfit,fcartfit,alpha,amass)
1529 
1530    do_itex : do itex = 1, lotfvar%nitex
1531     !--Compute rescaled velfit (center of mass speed and gaussian)
1532      call vel_rescale(mditemp,velfit,amass,v2gauss)
1533 
1534     ! start  verletvel here
1535      if(itime==1.AND.itex==1) then  ! check this
1536        vel_nexthalf(:,:) = velfit(:,:)
1537        xcart_next(:,:) = xcartfit(:,:)
1538      else
1539       !--Computation of vel_nexthalf (4.16 de Ref.1)
1540        call force_to_vel(v2gauss,dtion,amass,velfit,fcartfit,vel_nexthalf)
1541 
1542       !--Computation of the next positions
1543        xcart_next(:,:) = xcartfit(:,:)+vel_nexthalf(:,:)*dtion
1544 
1545 
1546       !--compute fcartfit and alpha (fcartfit is the accelleration)
1547        call force0(xcart_next,fcartfit,alpha,amass)
1548 
1549 
1550       !--Computation of vel(:,:) at the next positions
1551       !--Computation of v2gauss
1552        call vel_to_gauss(vel_nexthalf,amass,v2gauss)
1553 
1554       !--Compute velocity from force
1555        call force_to_vel(v2gauss,dtion,amass,vel_nexthalf,fcartfit,velfit)
1556      end if
1557      xcartfit = xcart_next
1558 
1559    end do do_itex
1560 
1561  end subroutine extrapolation_loop

lothpath/fitclus [ Functions ]

[ Top ] [ Functions ]

NAME

 fitclus

FUNCTION

INPUTS

PARENTS

      m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

268  subroutine fitclus(tfor,forc_in,tau0,alpha_tr,nqmin,nqmax)
269 
270   ! tfor : true if the forces are given in input
271   !        false if fitclus calculates the forces
272   ! forc_in(3,natom) : input forces to be fitted
273   ! tau0(3,natom) : atomic positions
274   ! alpha_tr(3,nbondex) : two body potential parameters
275   ! neighl(lotfvar%nneigx,natom) : list of neighbours
276   ! nneig(natom) : number of neighbours
277   ! nqmin : lowest  index for the quantum atoms (nqmin = 1 ) !!!!NOT USED!!!
278   ! nqmax : maximum index for the quantum atoms ( nqmax = natom)!!!!NOT USED!!!
279   ! niter  : iteration number (itime)
280 
281   use work_var_lotf,only : rcut,ifixed
282   use bond_lotf,only : nbondex,ifit,ibn_tots,nfitmax,imat,tafit,&
283 &                  ibnd_mat,ibn_tot,nfit
284   use tools_lotf,only : dlvsum
285   use glue_lotf, only : dphi,calc_coord_new_d,eval_u_n,eval_Upp_n
286   use eval_lotf,only : tuneparms,phi_n_calc,calc_coord2,eval_forces_U_n,eval_forces_U_n_2,eval_force_devs_new_d
287 
288   !--Evaluates "real" forces from external source, and computes
289   ! the istantaneous best fit of the current model.
290 
291 !This section has been created automatically by the script Abilint (TD).
292 !Do not modify the following lines by hand.
293 #undef ABI_FUNC
294 #define ABI_FUNC 'fitclus'
295 !End of the abilint section
296 
297   implicit none
298   !Arguments ------------------------------------
299   logical,intent(in) :: tfor
300   integer,intent(in) :: nqmin, nqmax
301   real(dp),intent(in) :: forc_in(3,lotfvar%natom)
302   real(dp) :: tau0(3,lotfvar%natom)
303   real(dp),intent(out) :: alpha_tr(3,nbondex)
304   !Local ----------------------------------------
305   real(dp), allocatable :: fspare(:,:),alpha_fce(:,:,:,:)
306   real(dp)              :: rcut_fit_int, rcf2_int
307   integer  :: i,j,k, iprecmass, iat
308   integer  :: jat
309   real(dp)  :: dcost,dcost_rms,dcost_old
310   integer :: ibn,nitmax,icheck,nwarn,ibn_count,ifiat
311   integer :: ifi
312   real(dp)  :: epotlotf_dum,f_dum,dampfact,dcost_rm0,dcost_rm1
313   real(dp)  :: d_dcost,dcost_check,dummy,vel,force
314   real(dp)  :: dtm,toeva,dtest
315   integer :: n,iem
316 
317   logical    :: tfit(3,nbondex),t_2nd(nbondex)
318   real(dp)   :: fact(nfitmax),fact2(nfitmax),ffit(3,nfitmax)
319 
320   real(dp), allocatable :: alpha_dum(:,:),alpha_old(:,:)
321   real(dp)   :: dcost_dalpha(3,nbondex),fcold(3,nbondex)
322 
323   real(dp) :: forc0_dum(3,0:nfitmax)
324   integer  ::  iwrdri
325   real(dp) :: dt_par(3),dmaspars(3)
326   real(dp),parameter :: prec_lotf = 1.0d-5
327 
328   !--Glue variables
329   integer :: istride, imin,imax,kat
330   integer :: nlist(0:lotfvar%nneigx)
331   integer :: neig2(lotfvar%nneigx),nlist2(lotfvar%nneigx,lotfvar%nneigx)
332   real(dp) :: tauv(3,lotfvar%nneigx),forcv(3,0:lotfvar%nneigx)
333   real(dp) :: stress(3,3)
334   real(dp) :: coordatom(lotfvar%natom),up_list(lotfvar%natom),upp_list(lotfvar%natom),forc_dum2(3)
335   real(dp) :: alpha_avg(3), alpha_disp(3)
336   real(dp) :: tauv2(3,lotfvar%nneigx,lotfvar%nneigx)
337   real(dp) :: rho_p_sum(3,nfitmax)
338 
339 ! *************************************************************************
340 
341   ! #####################  end DECLARATION #########################
342    call wrtout(std_out,' LOTF : FITCLUS','COLL')
343 
344    iwrdri = 0
345 
346    ABI_ALLOCATE(fspare,(3,nbondex))
347    ABI_ALLOCATE(alpha_fce,(3,2,nbondex,1))
348 
349   !--(0,-2) initialises a few parameters
350    rcut_fit_int = rcut
351    rcf2_int     = rcut_fit_int * rcut_fit_int
352 
353    dmaspars = (/ 0.01_dp, 0.1_dp,  0.1_dp /)
354    dt_par   = (/ 0.4_dp, 0.4_dp, 0.0_dp /)
355 
356   !--Not needed in this case... see if we can safely get rid of it
357   !     dt_ang          = zero
358   ! (0,-1) sets (optional: preconditions) masses of bond parameters
359 
360    iprecmass = 0
361 
362   !--Iitialize variables  :
363    t_2nd = .true.
364    tfit = .false.
365    ffit = zero
366    dcost_dalpha = zero
367 
368   !--set constant array for speed-up : fact(i)
369    do i = 1,nfit
370      iat = ifit(i)
371      fact2(i) =        real(ifixed(iat),dp)
372      fact(i)  = half * real(ifixed(iat),dp)
373    end do
374 
375   !--(0,0) decides which parameters will be fitted
376    call tuneparms(tau0,tfit,rcf2_int)
377 
378   !--(1,0) computes true forces: needs the external "expensive" routine.
379    fcold = zero
380 
381    if (tfor) then ! tfor = true <-> WE CALCULATE THE FORCES HERE !
382      do i =1, nfit
383        iat = ifit(i)
384       !write(97,*) iat,forc_in(1,iat)
385        ffit(:,i) = forc_in(:,iat)
386      end do
387    else
388      MSG_ERROR('LOTF : HERE WE SHOULD HAVE THE FORCES ALREADY !! ')
389    end if ! TFOR
390 
391   !--THE REST OF THE ROUTINE IS THE FIT
392    write(message,'(2(a,i8,a))')&
393    ' nbondex',  nbondex,ch10,' ibn_tots : ', ibn_tots,ch10
394    call wrtout(std_out,message,'COLL')
395 
396 
397    ABI_ALLOCATE(alpha_dum,(3,nbondex))
398    ABI_ALLOCATE(alpha_old,(3,nbondex))
399    alpha_dum = zero
400    alpha_old = zero
401 
402   !--(1,2) initialises the parameter set Alpha_dum, which IS modified
403   !       in the optimisation section
404    forall(ibn=1:ibn_tots)  alpha_dum(:,ibn) = (/ dphi,one, zero /)
405 
406   !--(2,0) computes derivatives of the ionic forces wrt parameters
407   !       at fixed ionic positions.
408 
409   !--Debug.......................
410    nitmax = 1000000
411 
412    icheck     = 0
413    nwarn      = 0
414    dcost      = zero
415    dcost_old  = zero
416    alpha_fce  = zero
417    dcost_rms  = one
418 
419   !###################################
420   !--MAIN MINIMISATION LOOP BEGINS
421   !###################################
422    main_minimization: do while(dcost_rms >  prec_lotf)
423     ! if(mod(icheck,20)==1.and.(lotfvar%me==1)) write(*,*)  icheck&
424     !   &,dcost,dcost_rms , prec_lotf,dcost_rms > prec_lotf
425 
426     !---------
427     !       derivatives of two-body forces w.r.t. parameters
428     !---------
429 
430     !--I choose to duplicate the code, can be made better :
431      if(lotfvar%classic==5) then
432        forc0_dum = zero
433       !--Testing........................................
434        alpha_fce = zero
435       !................................................
436       !--Test_5
437        dcost_dalpha = zero
438       ! End Test_5
439 
440       !--I cannot run only over fit pairs because of the term U(n)
441       ! so I run over fit atoms and ALL its neighbours twice, as I did in Force_Glue...
442 
443       !--Also: I only consider 1 parameter as variable, alpha(1,:)
444       ! therefore I don't use t_2nd nor tfit(:,ibn_count)
445 
446       !--(I) Pair potential (and its forces), coordination,  U, and  Up
447        nlist(:) = 0
448        stress(:,:) = zero
449        coordatom(:) = zero
450        up_list(:) = zero
451 
452       !--Parallel version
453        istride = lotfvar%natom/lotfvar%nproc
454        if(istride*lotfvar%nproc < lotfvar%natom) istride = istride + 1
455        imin = (lotfvar%me-1)*istride + 1
456        imax = lotfvar%me*istride
457        if(imax > lotfvar%natom) imax = lotfvar%natom
458 
459       !--First run over pairs: pair potential, its forces, coordination, U, total E, and Up
460        overatom1 : do iat = imin,imax
461          nlist(0) = iat
462         !   write(*,*) 'number of neighbours',nneig(iat)
463 
464          do j = 1,nneig(iat)
465            i = neighl(j,iat)
466            tauv(:,j) = tau0(:,i)
467            nlist(j) = i
468          end do
469 
470          if (tafit(iat)) then !--Only for fit atoms
471 
472           !--Pair potential (and its forces) & coordination(i)
473            call phi_n_calc(alpha_dum,nneig(iat),nlist,tau0(:,iat),&
474 &           tauv,epotlotf_dum,forcv,coordatom(iat),alpha_fce)
475 
476            forc0_dum(:,imat(iat)) = forc0_dum(:,imat(iat)) + forcv(:,0)
477 
478            do j = 1,nneig(iat)
479              i = neighl(j,iat)
480              forc0_dum(:,imat(i)) = forc0_dum(:,imat(i)) + forcv(:,j)
481            end do
482 
483            call eval_U_n(coordatom(iat),epotlotf_dum,up_list(iat))
484          else   !--non fit atoms
485 
486           !--Evaluate coordination for ALL the non fit atoms ... SHOULD BE IMPROVED to include only neighbours of neighbours of fit atoms...
487            call calc_coord2(nneig(iat),tau0(1,iat),tauv,coordatom(iat))
488 
489            call eval_U_n(coordatom(iat),epotlotf_dum,up_list(iat))
490 
491          end if !--fit / non fit atoms
492        end do overatom1
493 
494        call dlvsum(lotfvar%me-1,lotfvar%nproc,alpha_fce(1,1,1,1),6*nbondex)
495        call dlvsum(lotfvar%me-1,lotfvar%nproc,up_list(1),lotfvar%natom)
496 
497 
498       !--(II) Glue forces (coming only from the embedding function)
499        do iat = imin,imax
500          if (tafit(iat)) then
501            nlist(0) = iat
502            do j = 1,nneig(iat)
503              i = neighl(j,iat)
504              tauv(:,j) = tau0(:,i)
505              nlist(j) = i
506            end do
507 
508            call eval_forces_U_n(nneig(iat),nlist,tau0(1,iat),tauv,up_list,&
509 &           forc_dum2)
510            forc0_dum(:,imat(iat)) = forc0_dum(:,imat(iat)) + forc_dum2(:)
511          end if
512        end do
513      elseif (lotfvar%classic==6) then !-----------------------------------------------------------
514        forc0_dum = zero
515       !--Testing
516        alpha_fce = zero
517 
518       !--Test_5
519        dcost_dalpha = zero
520       !--End Test_5
521 
522       !--(I) Pair potential (and its forces), coordination,  U, and  Up
523        nlist(:) = 0
524        stress(:,:) = zero
525        coordatom(:) = zero
526        up_list(:) = zero
527        upp_list(:) = zero
528 
529       !--Parallel version
530        istride = lotfvar%natom/lotfvar%nproc
531        if(istride*lotfvar%nproc < lotfvar%natom) istride = istride + 1
532        imin = (lotfvar%me-1)*istride + 1
533        imax = lotfvar%me*istride
534        if(imax > lotfvar%natom) imax = lotfvar%natom
535 
536       !--First run over pairs: pair potential, its forces, coordination, U, total E, and Up
537        overatom2 : do iat = imin,imax
538          nlist(0) = iat
539          do j = 1,nneig(iat)
540            i = neighl(j,iat)
541            tauv(:,j) = tau0(:,i)
542            nlist(j) = i
543          end do
544 
545          if(tafit(iat)) then ! Only for fit atoms
546 
547           !--Pair potential (and its forces) & coordination(i)
548            call phi_n_calc(alpha_dum,nneig(iat),nlist,tau0(:,iat),&
549 &           tauv,epotlotf_dum,forcv,coordatom(iat),alpha_fce) ! COORDATOM is OK (lotfvar%classic dependence)
550            forc0_dum(:,imat(iat)) = forc0_dum(:,imat(iat)) + forcv(:,0)
551            do j = 1,nneig(iat)
552              i = neighl(j,iat)
553              forc0_dum(:,imat(i)) = forc0_dum(:,imat(i)) + forcv(:,j)
554            end do
555 
556           !--Up(coord(i)) & Upp(coord(i))
557            call eval_Upp_n(coordatom(iat),up_list(iat),upp_list(iat)) ! Up and Upp
558          else   ! non fit atoms
559 
560           !--Evaluate coordination for ALL the non fit atoms ... CAN BE IMPROVED to include only neighbours of neighbours of fit atoms...
561            call calc_coord2(nneig(iat),tau0(:,iat),tauv,&
562 &           coordatom(iat))                   ! OK (non-fit)
563 
564            call eval_Upp_n(coordatom(iat),up_list(iat),upp_list(iat)) ! OK (non-fit)
565 
566          end if ! fit / non fit atoms
567        end do overatom2
568 
569        call dlvsum(lotfvar%me-1,lotfvar%nproc,up_list(1),lotfvar%natom)
570        call dlvsum(lotfvar%me-1,lotfvar%nproc,upp_list(1),lotfvar%natom)
571        call dlvsum(lotfvar%me-1,lotfvar%nproc,alpha_fce(1,1,1,1),6*nbondex)
572 
573       !--(II) Glue forces (coming only from the embedding function)
574        do iat = imin,imax
575          if (tafit(iat)) then
576            nlist(0) = iat
577            do j = 1,nneig(iat)
578              jat = neighl(j,iat)
579              tauv(:,j) = tau0(:,jat)
580              nlist(j) = jat
581            end do
582 
583            call eval_forces_U_n_2(alpha_dum,nneig(iat),nlist,&
584 &           tau0(:,iat),tauv,up_list,&
585 &           rho_p_sum(:,imat(iat)),&
586 &           forc0_dum(:,imat(iat)))
587          end if
588        end do
589 
590        call dlvsum(lotfvar%me-1,lotfvar%nproc,rho_p_sum,3*nfitmax)
591        call dlvsum(lotfvar%me-1,lotfvar%nproc,forc0_dum,(nfitmax+1)*3)
592 
593       !--(III) Derivatives of the forces (inside DCOST_DALPHA)
594        neig2(:) = 0
595        nlist2(:,:) = 0
596        tauv(:,:) = zero
597        tauv2(:,:,:) = zero
598 
599        do iat = imin,imax
600          if (tafit(iat)) then
601            nlist(0) = iat
602            do j = 1,nneig(iat)
603              jat = neighl(j,iat)
604              tauv(:,j) = tau0(:,jat)
605              nlist(j) = jat
606              do k=1,nneig(jat)
607                neig2(j)=nneig(jat)
608                kat = neighl(k,jat)
609                nlist2(k,j)=kat
610                tauv2(:,k,j)=tau0(:,kat)
611              end do
612            end do
613 
614            call eval_force_devs_new_d(alpha_dum,nneig(iat),nlist,neig2,nlist2,&
615 &           tau0(:,iat),tauv,tauv2,up_list,upp_list,&
616 &           fact2,ffit,forc0_dum,rho_p_sum,dcost_dalpha)
617          end if
618        end do
619 
620       ! In paralell all this is going to cause trouble...........................
621       !      call dlvsum(lotfvar%me-1,lotfvar%nproc,dcost_dalpha,3*nbondex)
622       !......... check later more operations over dcost_dalpha...................
623 
624       !--------------------------------------------------------------------------------------------------
625 
626      end if ! lotfvar%classic
627 
628     ! if(lotfvar%classic /= 5) then
629     !   call dlvsum(lotfvar%me-1,lotfvar%nproc,forc0_dum,(nfitmax+1)*3)
630     ! end if
631 
632      if (lotfvar%classic /= 6)  then
633        call dlvsum(lotfvar%me-1,lotfvar%nproc,forc0_dum,(nfitmax+1)*3)
634      end if
635 
636     !--Check1
637     !--(1) evaluate cost function: quadratic deviation from true forces
638      dcost = zero
639      do i = 1,nfit
640        dcost = dcost+fact(i)*sum((forc0_dum(:,i)-ffit(:,i))**2,dim=1)
641      end do
642      dcost_rms = (sqrt((two*dcost)/nfit))*two*13.6058d0/0.529177d0
643 
644     !--minimisation is achived
645      if(dcost_rms < prec_lotf) exit
646 
647     !--(2) evaluate its derivative downhill w.r.t two body parms
648      do ibn_count = 1,ibn_tot
649 
650        if( (lotfvar%me == (mod(ibn_count-1,lotfvar%nproc)+1)) ) then
651 
652          iat   = ibnd_mat(1,ibn_count)
653          i     = ibnd_mat(2,ibn_count)
654          ifiat = imat(iat)
655          ifi   = imat(i)
656 
657          do n=1,2
658            f_dum = zero
659            if(tfit(n,ibn_count)) then
660             !--note : reduction here does not give good results so the
661             !loop in k has to be left
662              do k=1,3
663                f_dum = f_dum +  fact2(ifi)&
664 &               *  (forc0_dum(k,ifi)   - ffit(k,ifi))&
665 &               *  alpha_fce(k,n,ibn_count,1)
666 
667                f_dum = f_dum -  fact2(ifiat)&
668 &               *  (forc0_dum(k,ifiat) - ffit(k,ifiat))&
669 &               *  alpha_fce(k,n,ibn_count,1)
670              end do
671            end if
672 
673           ! Debug
674            dcost_dalpha(n,ibn_count) = dcost_dalpha(n,ibn_count) + f_dum
675          end do
676        end if
677      end do
678 
679 
680     !--parameter OPTIMIZATION STEP --------------------------- start
681      icheck = icheck + 1
682 
683      if(icheck ==1) alpha_old(:,:ibn_tot) = alpha_dum(:,:ibn_tot)
684 
685     !-----------------------------------------
686     !--damped dynamics
687 
688     !--(1) damping factor
689      if(icheck > 40) then
690        dampfact = 0.6
691        if(icheck > 200)  dampfact = 0.6
692      else
693        dampfact = 0.6
694      end if
695 
696      if(dcost  >  dcost_old)  nwarn = nwarn + 1
697 
698      dcost_rm0 = (sqrt((two*dcost_old)/nfit))*two*13.6058d0/0.529177d0
699      dcost_rm1 = (sqrt((two*dcost)/nfit))*two*13.6058d0/0.529177d0
700      d_dcost = dcost_rm1 - dcost_rm0
701      dcost_old = dcost
702 
703 
704      if(icheck < nitmax) then
705        if(icheck==nitmax/2) dcost_check = dcost_rms
706 
707       !--II BODY
708        bodyii: do  ibn_count = 1, ibn_tot
709          if( (lotfvar%me == (mod(ibn_count-1,lotfvar%nproc)+1)) ) then
710            if(.not.t_2nd(ibn_count)) cycle
711 
712           !--(2) updating strategy
713            do n=1,3
714              dummy = alpha_dum(n,ibn_count)
715              if(tfit(n,ibn_count)) then
716                vel   = (alpha_dum(n,ibn_count)-alpha_old(n,ibn_count))
717                force = dcost_dalpha(n,ibn_count)
718 
719                fcold(n,ibn_count) = force
720 
721               !--(3) Verlet step
722                alpha_dum(n,ibn_count) = alpha_dum(n,ibn_count)&
723 &               + dampfact * vel + dt_par(n) * dt_par(n) * force/dmaspars(n)
724 
725               !--(4) bound control for the parameters
726                if(alpha_dum(1,ibn_count)  >  6.1d0) then
727                  alpha_dum(1,ibn_count) = 6.1d0
728                elseif(alpha_dum(1,ibn_count)  <  two ) then
729                  alpha_dum(1,ibn_count) = two
730                 !MSG_ERROR('LOTF: Alpha1 reaches 2 au... Too small value!')
731                end if
732 
733              end if   ! tfit(n)
734 
735              alpha_old(n,ibn_count) = dummy
736            end do
737 
738          else
739            alpha_old(:,ibn_count) = zero
740            alpha_dum(:,ibn_count) = zero
741          end if
742 
743        end do bodyii
744 
745        call dlvsum(lotfvar%me-1,lotfvar%nproc,alpha_dum,3*ibn_tot)
746 
747        dcost_rms=(sqrt((two*dcost)/nfit))*two*13.6058d0/0.529177d0
748 
749       !--Minimization is not achived
750       !if(dcost_rms >  prec_lotf)  cycle
751      end if   ! if icheck < ##
752     !!qui!!
753     !--Minimization is not achived (dcost_rms > prec_lotf)
754    end do main_minimization
755 
756    if(dcost_rms >  prec_lotf) then
757      MSG_ERROR('LOTF: ACHTUNG: REQD.TOLERANCE NOT ACHIEVED IN THE FIT')
758    end if
759 
760    iwrdri = iwrdri + 1
761 
762   !--To have a look at what is going on with alpha...
763    if(lotfvar%me==1) then
764 
765      alpha_avg(:) = sum(alpha_dum(:,:ibn_tots),dim=2)
766      alpha_avg(:) = (one/ibn_tots)*alpha_avg(:)
767 
768      alpha_disp(:) = zero
769      do i=1,ibn_tots
770        alpha_disp(:) = alpha_disp(:) + (alpha_dum(:,i) - alpha_avg(:))**2
771      end do
772      alpha_disp(:) = sqrt((one/ibn_tots)*alpha_disp(:))
773 
774      write(message,'(2(2a,2f12.8))')ch10,&
775 &     'Alpha_avg =',(alpha_avg(i),i=1,2),ch10,&
776 &     'Alpha_disp =',(alpha_disp(i),i=1,2)
777      call wrtout(std_out,message,'COLL')
778 
779 
780      dtm = zero
781 
782      write(message,'(a,2f12.8,a)')'FIT.PREC. : ',dcost, dcost_rms,' EV/A '
783      call wrtout(std_out,message,'COLL')
784 
785      toeva = 2.d0*13.6058d0/0.529177d0
786      do i = 1,nfit
787        dtest = zero
788        do k=1,3
789          dtest = dtest+sqrt(fact2(i)*(forc0_dum(k,i)-ffit(k,i))**2)*toeva
790        end do
791        if(dtest > dtm) then
792          dtm = dtest
793          iem = ifit(i)
794        end if
795      end do
796      write(message,'(a,f12.8,a,i8)')'Max Error : ',dtm, ' on atom : ',iem
797      call wrtout(std_out,message,'COLL')
798 
799    end if
800 
801   !--PARAMETER OPTIMIZATION STEP ---------------------------end
802   !if (lotfvar%classic /= 5 .AND. lotfvar%classic /= 6) then
803   !  call dlvsum(lotfvar%me-1,lotfvar%nproc,alpha_dum,3*ibn_tot)
804   !end if
805 
806   !--Final result for alpha............................
807    if (lotfvar%me==1) then
808      call wrtout(std_out,'alpha(1)   alpha(2)','COLL')
809      do ibn_count = 1,ibn_tots
810        write(message,'(2f12.8)') (alpha_dum(n,ibn_count),n=1,2)
811        call wrtout(std_out,message,'COLL')
812      end do
813    end if
814 
815   !--prepares "true" updated parameters
816    alpha_tr(:,:ibn_tots) = alpha_dum(:,:ibn_tots)
817 
818    ABI_DEALLOCATE(alpha_fce)
819    ABI_DEALLOCATE(alpha_dum)
820    ABI_DEALLOCATE(alpha_old)
821   !-----------------------------------------------------------
822 
823  end subroutine fitclus

lothpath/force0 [ Functions ]

[ Top ] [ Functions ]

NAME

 force0

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

 970  subroutine force0(tau0,forc0,alpha_tr,amass)
 971   ! tau0(3,natom) : atomic positions (input)
 972   ! forc0(3,natom) : atomic forces(output)
 973   ! alpha_tr(3,nbondex) : two body potential parameters
 974   ! neighl(lotfvar%nneigx,natom) : list of neighbours
 975   ! nneig(natom) : number of neighbours
 976   ! epotlotf : potential energy (parameter dependent)
 977   use glue_lotf,only : eval_U_n
 978   use tools_lotf,only : dlvsum
 979   use bond_lotf,only : ibn_tot,nbondex,tafit
 980   use eval_lotf, only : phi_n_calc,calc_coord2,eval_forces_u_n
 981 
 982 !This section has been created automatically by the script Abilint (TD).
 983 !Do not modify the following lines by hand.
 984 #undef ABI_FUNC
 985 #define ABI_FUNC 'force0'
 986 !End of the abilint section
 987 
 988   implicit none
 989 
 990   !Arguments ------------------------------------
 991   real(dp),intent(in) :: amass(lotfvar%natom)
 992   real(dp),intent(in) :: tau0(3,lotfvar%natom)
 993   real(dp),intent(out) :: forc0(3,lotfvar%natom)
 994   real(dp),intent(in) :: alpha_tr(3,nbondex)
 995   !Local ----------------------------------------
 996   integer :: i, j
 997   integer :: istride,ibmin,ibmax,iat
 998   logical  :: tcalc(3)
 999   real(dp) :: epotlotf_dum
1000   real(dp) :: epotlotf_2
1001   real(dp) :: stress(3,3)
1002 
1003   !--Glue variables
1004   integer ::  nlist(0:lotfvar%nneigx)
1005   real(dp) ::  tauv(3,lotfvar%nneigx),forcv(3,0:lotfvar%nneigx)
1006   real(dp) ::  coordatom(lotfvar%natom),alpha_fdum_v(3,2,nbondex,1)
1007   real(dp) ::  up_list(lotfvar%natom),epotlotf_dum2,epotlotf_dum3
1008   real(dp) ::  forc_dum2(3)
1009 
1010 ! *************************************************************************
1011 
1012    epotlotf_2 = zero
1013    tcalc(:) = .false.
1014 
1015    forc0(:,:) = zero    ! MODifIED FOR ABINITttime 23/07/2008
1016 
1017   !--parallel version
1018    istride = ibn_tot/lotfvar%nproc
1019    if(istride*lotfvar%nproc < ibn_tot) istride = istride + 1
1020    ibmin = (lotfvar%me-1) * istride + 1
1021    ibmax =  lotfvar%me * istride
1022    if(ibmax > ibn_tot) ibmax = ibn_tot
1023 
1024 
1025   !--All but glue forces
1026    nlist(:) = 0
1027    epotlotf_dum = zero
1028    stress(:,:) = zero
1029    coordatom(:) = zero
1030    up_list(:) = zero
1031    epotlotf_dum2 = zero
1032 
1033    do iat=1,lotfvar%natom
1034 
1035      nlist(0) = iat
1036      do j = 1,nneig(iat)
1037        i = neighl(j,iat)
1038        tauv(:,j) = tau0(:,i)
1039        nlist(j) = i
1040      end do
1041 
1042      if (tafit(iat)) then
1043        call phi_n_calc(alpha_tr,nneig(iat),nlist,tau0(1,iat),&
1044        tauv,epotlotf_dum,forcv,coordatom(iat),alpha_fdum_v)
1045       !--PAIR energy: Fit atoms and NEIGHBOURS --> OK
1046        epotlotf_2 = epotlotf_2 + epotlotf_dum
1047 
1048        forc0(:,iat) = forc0(:,iat) + forcv(:,0)
1049        do j = 1,nneig(iat)
1050          i = neighl(j,iat)
1051          forc0(:,i) = forc0(:,i) + forcv(:,j)
1052        end do
1053 
1054        call eval_U_n(coordatom(iat),epotlotf_dum2,up_list(iat))
1055        epotlotf_2 = epotlotf_2 + epotlotf_dum2 ! GLUE energy: Fit atoms ONLY --> OK
1056 
1057      else
1058 
1059       !--Evaluate coordination and up_list for ALL the non fit atoms ... CAN BE IMPROVED to include only neighbours of neighbours of fit atoms...
1060       ! For example: run over neigbours, in case one of them is FIT atom, proceed, otherwise, do nothing
1061        call calc_coord2(nneig(iat),tau0(1,iat),tauv,coordatom(iat))
1062 
1063        call eval_U_n(coordatom(iat),epotlotf_dum3,up_list(iat))
1064       !--We do NOT accumulate epotlotf_dum3
1065      end if
1066 
1067    end do
1068 
1069   !--Glue forces
1070    do iat=1,lotfvar%natom
1071      if (tafit(iat)) then
1072        nlist(0) = iat
1073        do j = 1,nneig(iat)
1074          i = neighl(j,iat)
1075          tauv(:,j) = tau0(:,i)
1076          nlist(j) = i
1077        end do
1078        call eval_forces_U_n(nneig(iat),nlist,tau0(1,iat),tauv,up_list,forc_dum2)
1079 
1080        forc0(:,iat) = forc0(:,iat) + forc_dum2(:)
1081      end if
1082 
1083     !--Renomalization of the force (to get acceleration)
1084      forc0(:,iat) = forc0(:,iat)/amass(iat)
1085    end do
1086 
1087    epotlotf = epotlotf + epotlotf_2
1088   !--ends parallelisation
1089 
1090  end subroutine force0

lothpath/force_to_vel [ Functions ]

[ Top ] [ Functions ]

NAME

 force_to_vel

FUNCTION

  Compute velocity starting from : vel_in,v2gauss and forces

INPUTS

  v2gauss=gauss factor (2*kinetic energy)
  dtion=time step for Molecular Dynamics
  amass(natom)=masse of the ions
  vel_in(3,natom)=initial velocity of ions
  fcart(3,natom)=force on ions

 OUTPUTS
  vel_out(3,natom)=new velocity of ions

PARENTS

      m_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1262  subroutine force_to_vel(v2gauss,dtion,amass,vel_in,fcart,vel_out)
1263 
1264 
1265 !This section has been created automatically by the script Abilint (TD).
1266 !Do not modify the following lines by hand.
1267 #undef ABI_FUNC
1268 #define ABI_FUNC 'force_to_vel'
1269 !End of the abilint section
1270 
1271   implicit none
1272 
1273   !Arguments ------------------------------------
1274   real(dp),intent(in) :: v2gauss,dtion
1275   real(dp),intent(in) :: amass(:)
1276   real(dp),intent(in) :: vel_in(:,:)
1277   real(dp),intent(in) :: fcart(:,:)
1278   real(dp),intent(out) :: vel_out(:,:)
1279   !Local ---------------------------
1280   integer :: iatom,idim
1281   real(dp) :: a,b,sqb,as,s1,s2,s,scdot
1282 
1283 ! *************************************************************************
1284 
1285   !--Compute a and b (4.13 de Ref.1)
1286    a = zero
1287    b = zero
1288    do iatom=1,size(amass)
1289      do idim=1,3
1290        a=a+fcart(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1291        b=b+fcart(idim,iatom)*fcart(idim,iatom)*amass(iatom)
1292      end do
1293    end do
1294 
1295    a= a/v2gauss
1296    b= b/v2gauss
1297 
1298   !--Campute  s et scdot
1299    sqb = sqrt(b)
1300    as = sqb*dtion/2.
1301    s1 = cosh(as)
1302    s2 = sinh(as)
1303    s = a*(s1-1.)/b+s2/sqb
1304    scdot = a*s2/sqb+s1
1305    vel_out(:,:) = (vel_in(:,:)+fcart(:,:)*s)/scdot
1306  end subroutine force_to_vel

lothpath/intparms [ Functions ]

[ Top ] [ Functions ]

NAME

 intparms

FUNCTION

INPUTS

PARENTS

      m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1108  subroutine intparms(itime)
1109   use bond_lotf,only : ibn_tot
1110   use tools_lotf,only : pinterp,pinterp_nolinear
1111 
1112 !This section has been created automatically by the script Abilint (TD).
1113 !Do not modify the following lines by hand.
1114 #undef ABI_FUNC
1115 #define ABI_FUNC 'intparms'
1116 !End of the abilint section
1117 
1118   implicit none
1119   !Arguments ------------------------------------
1120   integer,intent(in) :: itime!,nitex
1121   !Local ----------------------------------------
1122   integer :: ibn,interp_type,nitdu
1123 
1124 ! *************************************************************************
1125 
1126    nitdu = mod(itime-lotfvar%n0,lotfvar%nitex) + 1
1127 
1128   !--Reset alpha
1129    alpha(:,:) = zero
1130 
1131   !--Select here the type of interpolation....................
1132    interp_type = 1
1133   !   interp_type = 2
1134   !--II body
1135    if(interp_type==1) then
1136      do ibn = 1, ibn_tot
1137        call pinterp(alpha_in(:,ibn),alpha_end(:,ibn),alpha(:,ibn),3,lotfvar%nitex,nitdu)
1138      end do
1139    end if
1140  end subroutine intparms

lothpath/lotf_extrapolation [ Functions ]

[ Top ] [ Functions ]

NAME

 lotf_extrapolation

FUNCTION

  return true if mod(itime,nitex) == 0

INPUTS

CHILDREN

SOURCE

1208  function lotf_extrapolation(itime)
1209 
1210 
1211 !This section has been created automatically by the script Abilint (TD).
1212 !Do not modify the following lines by hand.
1213 #undef ABI_FUNC
1214 #define ABI_FUNC 'lotf_extrapolation'
1215 !End of the abilint section
1216 
1217   implicit none
1218 
1219   !Arguments ------------------------------------
1220   integer,intent(in) :: itime
1221   logical :: lotf_extrapolation
1222 
1223 ! *************************************************************************
1224 
1225    if(itime == 1) then
1226      lotf_extrapolation = .true.
1227      return
1228    end if
1229    if(lotfvar%nitex-lotfvar%n0==0) then
1230      lotf_extrapolation = .true.
1231    else
1232      lotf_extrapolation = mod(itime-lotfvar%n0,lotfvar%nitex)==0
1233    end if
1234    return
1235  end function lotf_extrapolation

lothpath/lotf_interpolation [ Functions ]

[ Top ] [ Functions ]

NAME

 lotf_interpolation

FUNCTION

  Compute the LOTF interpolation:

INPUTS

  itime=numeber of MD step
  dtion=time step for Molecular Dynamics
  amass(natom)=masse of the ions
  xcart(3,natom)=position of the ions initial

 OUTPUTS
  xcart_next(3,natom)=positions of the ions final step
  vel_nexthalf(3,natom)=velocity of ions in next half step

SIDE EFFECTS

  v2gauss=gauss factor (twice the kinetic energy) (initial and final)
  vel(3,natom)=velocity of ions
  fcart_m(3,natom)=forces on ions in

PARENTS

      m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1594  subroutine lotf_interpolation(itime,dtion,v2gauss,amass,xcart,vel,&
1595 &                           fcart_m,xcart_next,vel_nexthalf)
1596 
1597 
1598 !This section has been created automatically by the script Abilint (TD).
1599 !Do not modify the following lines by hand.
1600 #undef ABI_FUNC
1601 #define ABI_FUNC 'lotf_interpolation'
1602 !End of the abilint section
1603 
1604   implicit none
1605 
1606   !Arguments ------------------------------------
1607   integer,intent(in) :: itime
1608   real(dp),intent(in) :: dtion
1609   real(dp),intent(inout) :: v2gauss
1610   real(dp),intent(in) :: amass(:)
1611   real(dp),intent(in) :: xcart(:,:)
1612   real(dp),intent(inout) :: vel(:,:)
1613   real(dp),intent(inout) :: fcart_m(:,:)
1614   real(dp),intent(out) :: vel_nexthalf(:,:)
1615   real(dp),intent(out) :: xcart_next(:,:)
1616 
1617 ! *************************************************************************
1618 
1619    write(message,'(a,i8)') ' ---LOTF interpolation: itime=',itime
1620    call wrtout(std_out,message,'COLL')
1621 
1622   !--VERLETVEL (interpolation)
1623    if(itime==lotfvar%n0) then  ! check this
1624     !--itime=0 fist step
1625      vel_nexthalf(:,:) = vel(:,:)
1626      xcart_next(:,:)  = xcart(:,:)
1627    else
1628     !--Computation of vel_nexthalf (4.16 de Ref.1)
1629      call force_to_vel(v2gauss,dtion,amass,vel,fcart_m,vel_nexthalf)
1630 
1631     !--Computation of the next positions
1632      xcart_next(:,:) = xcart(:,:)+vel_nexthalf(:,:)*dtion
1633 
1634     !--compute fcart_m and alpha (fcart_m is the accelleration)
1635      call force0(xcart_next,fcart_m,alpha,amass)
1636 
1637     !--Computation of vel(:,:) at the next positions
1638      call vel_to_gauss(vel_nexthalf,amass,v2gauss)
1639      call force_to_vel(v2gauss,dtion,amass,vel_nexthalf,fcart_m,vel)
1640 
1641    end if
1642 
1643  end subroutine lotf_interpolation
1644 
1645 end module LOTFPATH

lothpath/upd_lis [ Functions ]

[ Top ] [ Functions ]

NAME

 upd_lis

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

843  subroutine  upd_lis(tau0,niter,alpha_dum)
844   ! lotfvar%nneigx : max number of neighbours
845   ! tau0(3,natom) : atomic positions
846   ! neighl(lotfvar%nneigx,natom) : list of neighbours
847   ! neighl_old(lotfvar%nneigx,natom) : old list of neighbours
848   ! nneig_old(natom) : old number of neighbours
849   ! niter  : iteration number (itime)
850 
851   !  updates the neighbour list for each atom, and updates the active
852   !  parameter lists
853 
854   USE GLUE_LOTF,only : dphi
855   use work_var_lotf,only : smallfit
856   use bond_lotf,only : ibn_tot,ibn_tot2,ibn_tots,ibmat,ibnd_mat,&
857 &                  nfitmax,imat,bond_matrix_set,bond_compute,&
858 &                  bond_fit_set,nbondex
859   use eval_lotf,only : upd_lis0
860 
861 !This section has been created automatically by the script Abilint (TD).
862 !Do not modify the following lines by hand.
863 #undef ABI_FUNC
864 #define ABI_FUNC 'upd_lis'
865 !End of the abilint section
866 
867   implicit none
868 
869   !Arguments ------------------------------------
870   integer,intent(in)   :: niter
871   real(dp),intent(inout) :: tau0(3,lotfvar%natom)
872   real(dp),intent(inout) :: alpha_dum(3,nbondex)
873   !Local ----------------------------------------
874   integer :: j,ibn_count
875   integer :: ifo,jato,jb_old
876   integer :: iat,jat,nfitdum
877   integer :: jbo(nbondex)
878   !take care nsurf_at*(nbond_at)/2
879 
880 ! *************************************************************************
881 
882   !--INITIALIZATION AND DECISION :
883    jbo(:) = 0
884 
885   !--neighbours update
886    call upd_lis0(tau0,neighl,nneig,niter)
887 
888   !--set tafit and compute nfitdum
889    call SMALLFIT(tau0,nfitdum)
890 
891   !--control nfitmax and set ifit,nfit
892    call bond_fit_set(lotfvar%natom,nfitdum)
893 
894   !--Updates bond matrix, associates the bond to the atom neighlists
895    call bond_compute(nneig,neighl)
896 
897   !--CREATES THE REORDERING ARRAY (ibmat)
898   !  ONLY variational BONDS ARE CONSIDERED IN HISTORY
899    do ibn_count = 1, ibn_tot  ! only on fitted atoms pairs
900     !--finds to which atoms it corresponds
901      iat = ibnd_mat(1,ibn_count)
902      jat = ibnd_mat(2,ibn_count)
903 
904      ifo = imat(iat)    !--imat finds the old atomic 'fitted number'
905      if(iat > jat) MSG_ERROR('UPDLIS 177')
906 
907     !--Set to 0, finds to which old bond (if any) these two correspond
908      jb_old = 0               !--atom jat is a new neighbour of atom iat
909      do j =1,nneig_old(iat)
910        jato = neighl_old(j,iat)
911        if(jato==jat) then
912          jb_old = ibmat(j,ifo)  !--atom jat was already a neighbour of atom iat
913          exit
914        end if
915      end do
916      jbo(ibn_count) = jb_old  !--index array for reordering
917    end do
918 
919   !--updates bond matrix, associates the bond to the atom neighlists
920    call bond_matrix_set(nneig,neighl)
921 
922 
923    write(message,'(2a,2(a,i8,a),(a,2i8,a),a,i8,a)') &
924 &   ' LOTF : UPD_LIS',ch10,&
925 &   ' ITERATION:  ',niter,ch10,&
926 &   ' NO. OF BONDS IN ACTIVE ZONE : ',ibn_tot,ch10,&
927 &   ' BORDER BOUNDS (&max) : ', ibn_tot2,6*nfitmax ,ch10,&
928 &   ' TOTAL N.OF BOUNDS : ', ibn_tots,ch10
929    call wrtout(std_out,message,'COLL')
930 
931   !--updates old neighbour lists
932    do iat =1,lotfvar%natom
933      nneig_old(iat) = nneig(iat)
934      neighl_old(:nneig(iat),iat) = neighl(:nneig(iat),iat)
935     !write(6,*) iat,nneig(iat),(neighl(j,iat),j=1,nneig(iat))
936    end do
937 
938   !--updates parameter list
939    if(niter /= lotfvar%n0) then
940      call alpha_update(dphi,jbo,alpha_dum)
941    end if
942 
943    if (ibn_tots > nbondex) then
944      write(message,'(2a,2(a,i8))') 'LOTF: ibn_tots > nbondex  ! ',ch10,&
945 &     'UPDLIS  stop : IBNTOTS = ',ibn_tots,' NBONDEX = ',nbondex
946      MSG_ERROR(message)
947    end if
948 
949 
950  end subroutine upd_lis

lothpath/vel_rescale [ Functions ]

[ Top ] [ Functions ]

NAME

 vel_rescale

FUNCTION

  Starting from the velocity, it recompute the velocities in the
  center of mass. Then compute the gauss factor and renormalises
  the velocities with respect the gauss distribution. Then is
  recompute the gauss factor

INPUTS

  mditemp=temperature of ions
  amass(natom)=masse of the ions

 OUTPUTS
  v2gauss=2*kinetic energy

SIDE EFFECTS

  vel(3,natom)=velocity of ions

PARENTS

      m_lotf,m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1403  subroutine vel_rescale(mditemp,vel,amass,v2gauss)
1404 
1405 
1406 !This section has been created automatically by the script Abilint (TD).
1407 !Do not modify the following lines by hand.
1408 #undef ABI_FUNC
1409 #define ABI_FUNC 'vel_rescale'
1410 !End of the abilint section
1411 
1412   real(dp),intent(in) :: mditemp
1413   real(dp),intent(out) :: v2gauss
1414   real(dp),intent(in) :: amass(:)
1415   real(dp),intent(inout) :: vel(:,:)
1416 
1417   integer :: idim,natom
1418   real(dp) :: mass_tot,momentum,rescale_vel,vtest,sigma2
1419 
1420 ! *************************************************************************
1421 
1422    natom = size(amass)
1423 
1424   !--Get rid of center-of-mass velocity
1425    mass_tot = sum(amass(:))
1426    do idim=1,3
1427      momentum = sum(amass(:)*vel(idim,:))
1428      vel(idim,:) = vel(idim,:)-momentum/mass_tot
1429    end do
1430 
1431   !--Recompute v2gauss
1432    call vel_to_gauss(vel,amass,v2gauss,vtest)
1433 
1434   !--Now rescale the velocities to give the exact temperature
1435    rescale_vel = sqrt(three*natom*kb_HaK*mditemp/v2gauss)
1436    vel(:,:) = vel(:,:)*rescale_vel
1437 
1438   !--Recompute v2gauss
1439    call vel_to_gauss(vel,amass,v2gauss)
1440 
1441   !--Compute the variance and print
1442    sigma2 = (v2gauss/(3._dp*natom)-amass(1)*vtest**2)/kb_HaK
1443 
1444    write(message, '(a,d12.5,a,D12.5)' )&
1445 &   ' --- LOTF STEP : Effective temperature',&
1446 &   v2gauss/(3*natom*kb_HaK),' From variance', sigma2
1447    call wrtout(ab_out,message,'COLL')
1448    call wrtout(std_out,message,'COLL')
1449 
1450 
1451  end subroutine vel_rescale

lothpath/vel_to_gauss [ Functions ]

[ Top ] [ Functions ]

NAME

 vel_to_gauss

FUNCTION

  Compute gauss factor sum(v**2*m) the double of kinetic energy
  If present vtest compute also sum(v)/(3*sum(m))

INPUTS

  vel_in(3,natom)=velocity to use
  amass(natom)=masse of the ions

OUTPUT

  v2gauss=2*kinetic energy
  vtest=pick velocity

PARENTS

      m_lotf,m_pred_lotf

CHILDREN

      force0,force_to_vel,vel_to_gauss,wrtout

SOURCE

1332  subroutine vel_to_gauss(vel_in,amass,v2gauss,vtest)
1333 
1334 
1335 !This section has been created automatically by the script Abilint (TD).
1336 !Do not modify the following lines by hand.
1337 #undef ABI_FUNC
1338 #define ABI_FUNC 'vel_to_gauss'
1339 !End of the abilint section
1340 
1341   implicit none
1342 
1343   !Arguments ------------------------------------
1344   real(dp),intent(out) :: v2gauss
1345   real(dp),intent(out),optional :: vtest
1346   real(dp),intent(in) :: vel_in(:,:)
1347   real(dp),intent(in) :: amass(:)
1348   !Local ---------------------------
1349   integer :: iatom,idim
1350 
1351 ! *************************************************************************
1352 
1353    if(present(vtest)) then
1354      v2gauss = zero
1355      vtest = zero
1356      do iatom=1,size(amass)
1357        do idim=1,3
1358          v2gauss = v2gauss+vel_in(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1359          vtest = vtest+vel_in(idim,iatom)
1360        end do
1361      end do
1362      vtest = vtest/(3._dp*size(amass))
1363    else
1364     !--Recompute v2gauss with the rescaled velocities
1365      v2gauss = zero
1366      do iatom=1,size(amass)
1367        do idim=1,3
1368          v2gauss = v2gauss+vel_in(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1369        end do
1370      end do
1371    end if
1372  end subroutine vel_to_gauss