TABLE OF CONTENTS


ABINIT/lotfpath [ Modules ]

[ Top ] [ Modules ]

NAME

 lotfpath

FUNCTION

COPYRIGHT

 Copyright (C) 2005-2022 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

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

lotfpath/init_lotf [ Functions ]

[ Top ] [ lotfpath ] [ Functions ]

NAME

 init_lotf

FUNCTION

INPUTS

SOURCE

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

SOURCE

1076  subroutine alpha_update(dphi,jbo,alpha_dum)
1077   use bond_lotf,only : ibn_tot
1078   implicit none
1079 
1080   !Arguments ------------------------
1081   real(dp),intent(in) :: dphi
1082   integer,intent(in) :: jbo(:)
1083   real(dp),intent(inout) :: alpha_dum(:,:)
1084   !Local ---------------------------
1085   integer :: ibn,jb_old
1086   real(dp) :: alphawork(3,ibn_tot)
1087 
1088 ! *************************************************************************
1089 
1090    do ibn = 1,ibn_tot
1091      jb_old = jbo(ibn)
1092      if(jb_old /= 0) then
1093       !--swaps old parms
1094        alphawork(:,ibn) = alpha_dum(:,jb_old)
1095      else
1096       !--or reinitialise bond parms
1097        alphawork(:,ibn) = (/ dphi, one, zero /)
1098      end if
1099    end do
1100 
1101    alpha_dum(:,:ibn_tot) = alphawork(:,:)
1102 
1103  end subroutine alpha_update

lothpath/end_lotf [ Functions ]

[ Top ] [ Functions ]

NAME

 end_lotf

FUNCTION

INPUTS

SOURCE

201  subroutine end_LOTF()
202 
203   use work_var_lotf
204   use bond_lotf
205 
206   !--init_lotf
207 ! *************************************************************************
208    ABI_FREE(alpha)
209    ABI_FREE(alpha_in)
210    ABI_FREE(alpha_end)
211    ABI_FREE(nneig)
212    ABI_FREE(neighl)
213    ABI_FREE(fcartfit)
214    ABI_FREE(xcartfit)
215    ABI_FREE(velfit)
216 
217   !--deallocate LOTF internal variables
218    call work_var_dealloc()
219   !--deallocate LOTF bond variables
220    call bond_dealloc()
221 
222  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

SOURCE

1344  subroutine extrapolation_loop(itime,mditemp,dtion,amass,&
1345 &                           xcart,vel,xcart_next,vel_nexthalf)
1346 
1347   implicit none
1348 
1349   !Arguments ------------------------------------
1350   integer,intent(in) :: itime
1351   real(dp),intent(in) :: mditemp,dtion
1352   real(dp),intent(in) :: amass(:)
1353   real(dp),intent(in) :: vel(:,:)
1354   real(dp),intent(in) :: xcart(:,:)
1355   real(dp),intent(out) :: vel_nexthalf(:,:)
1356   real(dp),intent(out) :: xcart_next(:,:)
1357   !Local ---------------------------
1358   integer :: itex
1359   real(dp) :: v2gauss
1360 
1361 ! *************************************************************************
1362 
1363   !--inital values for position and forces
1364    xcartfit(:,:) = xcart(:,:)
1365    velfit(:,:) = vel(:,:)
1366 
1367   !--Update bond variables and alpha
1368    call upd_lis(xcartfit,itime,alpha_in)
1369 
1370   !print *,'b-end uplis',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2))
1371 
1372   !--Store alpha in alpha_in
1373    alpha(:,:) = alpha_in(:,:)
1374 
1375   !--compute fcartfit (fcartfit is the accelleration)
1376    call force0(xcartfit,fcartfit,alpha,amass)
1377 
1378    do_itex : do itex = 1, lotfvar%nitex
1379     !--Compute rescaled velfit (center of mass speed and gaussian)
1380      call vel_rescale(mditemp,velfit,amass,v2gauss)
1381 
1382     ! start  verletvel here
1383      if(itime==1.AND.itex==1) then  ! check this
1384        vel_nexthalf(:,:) = velfit(:,:)
1385        xcart_next(:,:) = xcartfit(:,:)
1386      else
1387       !--Computation of vel_nexthalf (4.16 de Ref.1)
1388        call force_to_vel(v2gauss,dtion,amass,velfit,fcartfit,vel_nexthalf)
1389 
1390       !--Computation of the next positions
1391        xcart_next(:,:) = xcartfit(:,:)+vel_nexthalf(:,:)*dtion
1392 
1393 
1394       !--compute fcartfit and alpha (fcartfit is the accelleration)
1395        call force0(xcart_next,fcartfit,alpha,amass)
1396 
1397 
1398       !--Computation of vel(:,:) at the next positions
1399       !--Computation of v2gauss
1400        call vel_to_gauss(vel_nexthalf,amass,v2gauss)
1401 
1402       !--Compute velocity from force
1403        call force_to_vel(v2gauss,dtion,amass,vel_nexthalf,fcartfit,velfit)
1404      end if
1405      xcartfit = xcart_next
1406 
1407    end do do_itex
1408 
1409  end subroutine extrapolation_loop

lothpath/fitclus [ Functions ]

[ Top ] [ Functions ]

NAME

 fitclus

FUNCTION

INPUTS

SOURCE

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

lothpath/force0 [ Functions ]

[ Top ] [ Functions ]

NAME

 force0

FUNCTION

INPUTS

SOURCE

 912  subroutine force0(tau0,forc0,alpha_tr,amass)
 913   ! tau0(3,natom) : atomic positions (input)
 914   ! forc0(3,natom) : atomic forces(output)
 915   ! alpha_tr(3,nbondex) : two body potential parameters
 916   ! neighl(lotfvar%nneigx,natom) : list of neighbours
 917   ! nneig(natom) : number of neighbours
 918   ! epotlotf : potential energy (parameter dependent)
 919   use glue_lotf,only : eval_U_n
 920   use tools_lotf,only : dlvsum
 921   use bond_lotf,only : ibn_tot,nbondex,tafit
 922   use eval_lotf, only : phi_n_calc,calc_coord2,eval_forces_u_n
 923   implicit none
 924 
 925   !Arguments ------------------------------------
 926   real(dp),intent(in) :: amass(lotfvar%natom)
 927   real(dp),intent(in) :: tau0(3,lotfvar%natom)
 928   real(dp),intent(out) :: forc0(3,lotfvar%natom)
 929   real(dp),intent(in) :: alpha_tr(3,nbondex)
 930   !Local ----------------------------------------
 931   integer :: i, j
 932   integer :: istride,ibmin,ibmax,iat
 933   logical  :: tcalc(3)
 934   real(dp) :: epotlotf_dum
 935   real(dp) :: epotlotf_2
 936   real(dp) :: stress(3,3)
 937 
 938   !--Glue variables
 939   integer ::  nlist(0:lotfvar%nneigx)
 940   real(dp) ::  tauv(3,lotfvar%nneigx),forcv(3,0:lotfvar%nneigx)
 941   real(dp) ::  coordatom(lotfvar%natom),alpha_fdum_v(3,2,nbondex,1)
 942   real(dp) ::  up_list(lotfvar%natom),epotlotf_dum2,epotlotf_dum3
 943   real(dp) ::  forc_dum2(3)
 944 
 945 ! *************************************************************************
 946 
 947    epotlotf_2 = zero
 948    tcalc(:) = .false.
 949 
 950    forc0(:,:) = zero    ! MODifIED FOR ABINITttime 23/07/2008
 951 
 952   !--parallel version
 953    istride = ibn_tot/lotfvar%nproc
 954    if(istride*lotfvar%nproc < ibn_tot) istride = istride + 1
 955    ibmin = (lotfvar%me-1) * istride + 1
 956    ibmax =  lotfvar%me * istride
 957    if(ibmax > ibn_tot) ibmax = ibn_tot
 958 
 959 
 960   !--All but glue forces
 961    nlist(:) = 0
 962    epotlotf_dum = zero
 963    stress(:,:) = zero
 964    coordatom(:) = zero
 965    up_list(:) = zero
 966    epotlotf_dum2 = zero
 967 
 968    do iat=1,lotfvar%natom
 969 
 970      nlist(0) = iat
 971      do j = 1,nneig(iat)
 972        i = neighl(j,iat)
 973        tauv(:,j) = tau0(:,i)
 974        nlist(j) = i
 975      end do
 976 
 977      if (tafit(iat)) then
 978        call phi_n_calc(alpha_tr,nneig(iat),nlist,tau0(1,iat),&
 979        tauv,epotlotf_dum,forcv,coordatom(iat),alpha_fdum_v)
 980       !--PAIR energy: Fit atoms and NEIGHBOURS --> OK
 981        epotlotf_2 = epotlotf_2 + epotlotf_dum
 982 
 983        forc0(:,iat) = forc0(:,iat) + forcv(:,0)
 984        do j = 1,nneig(iat)
 985          i = neighl(j,iat)
 986          forc0(:,i) = forc0(:,i) + forcv(:,j)
 987        end do
 988 
 989        call eval_U_n(coordatom(iat),epotlotf_dum2,up_list(iat))
 990        epotlotf_2 = epotlotf_2 + epotlotf_dum2 ! GLUE energy: Fit atoms ONLY --> OK
 991 
 992      else
 993 
 994       !--Evaluate coordination and up_list for ALL the non fit atoms ... CAN BE IMPROVED to include only neighbours of neighbours of fit atoms...
 995       ! For example: run over neigbours, in case one of them is FIT atom, proceed, otherwise, do nothing
 996        call calc_coord2(nneig(iat),tau0(1,iat),tauv,coordatom(iat))
 997 
 998        call eval_U_n(coordatom(iat),epotlotf_dum3,up_list(iat))
 999       !--We do NOT accumulate epotlotf_dum3
1000      end if
1001 
1002    end do
1003 
1004   !--Glue forces
1005    do iat=1,lotfvar%natom
1006      if (tafit(iat)) then
1007        nlist(0) = iat
1008        do j = 1,nneig(iat)
1009          i = neighl(j,iat)
1010          tauv(:,j) = tau0(:,i)
1011          nlist(j) = i
1012        end do
1013        call eval_forces_U_n(nneig(iat),nlist,tau0(1,iat),tauv,up_list,forc_dum2)
1014 
1015        forc0(:,iat) = forc0(:,iat) + forc_dum2(:)
1016      end if
1017 
1018     !--Renomalization of the force (to get acceleration)
1019      forc0(:,iat) = forc0(:,iat)/amass(iat)
1020    end do
1021 
1022    epotlotf = epotlotf + epotlotf_2
1023   !--ends parallelisation
1024 
1025  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

SOURCE

1156  subroutine force_to_vel(v2gauss,dtion,amass,vel_in,fcart,vel_out)
1157 
1158   implicit none
1159 
1160   !Arguments ------------------------------------
1161   real(dp),intent(in) :: v2gauss,dtion
1162   real(dp),intent(in) :: amass(:)
1163   real(dp),intent(in) :: vel_in(:,:)
1164   real(dp),intent(in) :: fcart(:,:)
1165   real(dp),intent(out) :: vel_out(:,:)
1166   !Local ---------------------------
1167   integer :: iatom,idim
1168   real(dp) :: a,b,sqb,as,s1,s2,s,scdot
1169 
1170 ! *************************************************************************
1171 
1172   !--Compute a and b (4.13 de Ref.1)
1173    a = zero
1174    b = zero
1175    do iatom=1,size(amass)
1176      do idim=1,3
1177        a=a+fcart(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1178        b=b+fcart(idim,iatom)*fcart(idim,iatom)*amass(iatom)
1179      end do
1180    end do
1181 
1182    a= a/v2gauss
1183    b= b/v2gauss
1184 
1185   !--Campute  s et scdot
1186    sqb = sqrt(b)
1187    as = sqb*dtion/2.
1188    s1 = cosh(as)
1189    s2 = sinh(as)
1190    s = a*(s1-1.)/b+s2/sqb
1191    scdot = a*s2/sqb+s1
1192    vel_out(:,:) = (vel_in(:,:)+fcart(:,:)*s)/scdot
1193  end subroutine force_to_vel

lothpath/intparms [ Functions ]

[ Top ] [ Functions ]

NAME

 intparms

FUNCTION

INPUTS

SOURCE

1037  subroutine intparms(itime)
1038   use bond_lotf,only : ibn_tot
1039   use tools_lotf,only : pinterp,pinterp_nolinear
1040   implicit none
1041   !Arguments ------------------------------------
1042   integer,intent(in) :: itime!,nitex
1043   !Local ----------------------------------------
1044   integer :: ibn,interp_type,nitdu
1045 
1046 ! *************************************************************************
1047 
1048    nitdu = mod(itime-lotfvar%n0,lotfvar%nitex) + 1
1049 
1050   !--Reset alpha
1051    alpha(:,:) = zero
1052 
1053   !--Select here the type of interpolation....................
1054    interp_type = 1
1055   !   interp_type = 2
1056   !--II body
1057    if(interp_type==1) then
1058      do ibn = 1, ibn_tot
1059        call pinterp(alpha_in(:,ibn),alpha_end(:,ibn),alpha(:,ibn),3,lotfvar%nitex,nitdu)
1060      end do
1061    end if
1062  end subroutine intparms

lothpath/lotf_extrapolation [ Functions ]

[ Top ] [ Functions ]

NAME

 lotf_extrapolation

FUNCTION

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

INPUTS

SOURCE

1115  function lotf_extrapolation(itime)
1116 
1117   implicit none
1118 
1119   !Arguments ------------------------------------
1120   integer,intent(in) :: itime
1121   logical :: lotf_extrapolation
1122 
1123 ! *************************************************************************
1124 
1125    if(itime == 1) then
1126      lotf_extrapolation = .true.
1127      return
1128    end if
1129    if(lotfvar%nitex-lotfvar%n0==0) then
1130      lotf_extrapolation = .true.
1131    else
1132      lotf_extrapolation = mod(itime-lotfvar%n0,lotfvar%nitex)==0
1133    end if
1134    return
1135  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

SOURCE

1436  subroutine lotf_interpolation(itime,dtion,v2gauss,amass,xcart,vel,&
1437 &                           fcart_m,xcart_next,vel_nexthalf)
1438 
1439   implicit none
1440 
1441   !Arguments ------------------------------------
1442   integer,intent(in) :: itime
1443   real(dp),intent(in) :: dtion
1444   real(dp),intent(inout) :: v2gauss
1445   real(dp),intent(in) :: amass(:)
1446   real(dp),intent(in) :: xcart(:,:)
1447   real(dp),intent(inout) :: vel(:,:)
1448   real(dp),intent(inout) :: fcart_m(:,:)
1449   real(dp),intent(out) :: vel_nexthalf(:,:)
1450   real(dp),intent(out) :: xcart_next(:,:)
1451 
1452 ! *************************************************************************
1453 
1454    write(message,'(a,i8)') ' ---LOTF interpolation: itime=',itime
1455    call wrtout(std_out,message,'COLL')
1456 
1457   !--VERLETVEL (interpolation)
1458    if(itime==lotfvar%n0) then  ! check this
1459     !--itime=0 fist step
1460      vel_nexthalf(:,:) = vel(:,:)
1461      xcart_next(:,:)  = xcart(:,:)
1462    else
1463     !--Computation of vel_nexthalf (4.16 de Ref.1)
1464      call force_to_vel(v2gauss,dtion,amass,vel,fcart_m,vel_nexthalf)
1465 
1466     !--Computation of the next positions
1467      xcart_next(:,:) = xcart(:,:)+vel_nexthalf(:,:)*dtion
1468 
1469     !--compute fcart_m and alpha (fcart_m is the accelleration)
1470      call force0(xcart_next,fcart_m,alpha,amass)
1471 
1472     !--Computation of vel(:,:) at the next positions
1473      call vel_to_gauss(vel_nexthalf,amass,v2gauss)
1474      call force_to_vel(v2gauss,dtion,amass,vel_nexthalf,fcart_m,vel)
1475 
1476    end if
1477 
1478  end subroutine lotf_interpolation
1479 
1480 end module LOTFPATH

lothpath/upd_lis [ Functions ]

[ Top ] [ Functions ]

NAME

 upd_lis

FUNCTION

INPUTS

SOURCE

798  subroutine  upd_lis(tau0,niter,alpha_dum)
799   ! lotfvar%nneigx : max number of neighbours
800   ! tau0(3,natom) : atomic positions
801   ! neighl(lotfvar%nneigx,natom) : list of neighbours
802   ! neighl_old(lotfvar%nneigx,natom) : old list of neighbours
803   ! nneig_old(natom) : old number of neighbours
804   ! niter  : iteration number (itime)
805 
806   !  updates the neighbour list for each atom, and updates the active
807   !  parameter lists
808 
809   USE GLUE_LOTF,only : dphi
810   use work_var_lotf,only : smallfit
811   use bond_lotf,only : ibn_tot,ibn_tot2,ibn_tots,ibmat,ibnd_mat,&
812 &                  nfitmax,imat,bond_matrix_set,bond_compute,&
813 &                  bond_fit_set,nbondex
814   use eval_lotf,only : upd_lis0
815   implicit none
816 
817   !Arguments ------------------------------------
818   integer,intent(in)   :: niter
819   real(dp),intent(inout) :: tau0(3,lotfvar%natom)
820   real(dp),intent(inout) :: alpha_dum(3,nbondex)
821   !Local ----------------------------------------
822   integer :: j,ibn_count
823   integer :: ifo,jato,jb_old
824   integer :: iat,jat,nfitdum
825   integer :: jbo(nbondex)
826   !take care nsurf_at*(nbond_at)/2
827 
828 ! *************************************************************************
829 
830   !--INITIALIZATION AND DECISION :
831    jbo(:) = 0
832 
833   !--neighbours update
834    call upd_lis0(tau0,neighl,nneig,niter)
835 
836   !--set tafit and compute nfitdum
837    call SMALLFIT(tau0,nfitdum)
838 
839   !--control nfitmax and set ifit,nfit
840    call bond_fit_set(lotfvar%natom,nfitdum)
841 
842   !--Updates bond matrix, associates the bond to the atom neighlists
843    call bond_compute(nneig,neighl)
844 
845   !--CREATES THE REORDERING ARRAY (ibmat)
846   !  ONLY variational BONDS ARE CONSIDERED IN HISTORY
847    do ibn_count = 1, ibn_tot  ! only on fitted atoms pairs
848     !--finds to which atoms it corresponds
849      iat = ibnd_mat(1,ibn_count)
850      jat = ibnd_mat(2,ibn_count)
851 
852      ifo = imat(iat)    !--imat finds the old atomic 'fitted number'
853      if(iat > jat) ABI_ERROR('UPDLIS 177')
854 
855     !--Set to 0, finds to which old bond (if any) these two correspond
856      jb_old = 0               !--atom jat is a new neighbour of atom iat
857      do j =1,nneig_old(iat)
858        jato = neighl_old(j,iat)
859        if(jato==jat) then
860          jb_old = ibmat(j,ifo)  !--atom jat was already a neighbour of atom iat
861          exit
862        end if
863      end do
864      jbo(ibn_count) = jb_old  !--index array for reordering
865    end do
866 
867   !--updates bond matrix, associates the bond to the atom neighlists
868    call bond_matrix_set(nneig,neighl)
869 
870 
871    write(message,'(2a,2(a,i8,a),(a,2i8,a),a,i8,a)') &
872 &   ' LOTF : UPD_LIS',ch10,&
873 &   ' ITERATION:  ',niter,ch10,&
874 &   ' NO. OF BONDS IN ACTIVE ZONE : ',ibn_tot,ch10,&
875 &   ' BORDER BOUNDS (&max) : ', ibn_tot2,6*nfitmax ,ch10,&
876 &   ' TOTAL N.OF BOUNDS : ', ibn_tots,ch10
877    call wrtout(std_out,message,'COLL')
878 
879   !--updates old neighbour lists
880    do iat =1,lotfvar%natom
881      nneig_old(iat) = nneig(iat)
882      neighl_old(:nneig(iat),iat) = neighl(:nneig(iat),iat)
883     !write(6,*) iat,nneig(iat),(neighl(j,iat),j=1,nneig(iat))
884    end do
885 
886   !--updates parameter list
887    if(niter /= lotfvar%n0) then
888      call alpha_update(dphi,jbo,alpha_dum)
889    end if
890 
891    if (ibn_tots > nbondex) then
892      write(message,'(2a,2(a,i8))') 'LOTF: ibn_tots > nbondex  ! ',ch10,&
893 &     'UPDLIS  stop : IBNTOTS = ',ibn_tots,' NBONDEX = ',nbondex
894      ABI_ERROR(message)
895    end if
896 
897 
898  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

SOURCE

1271  subroutine vel_rescale(mditemp,vel,amass,v2gauss)
1272 
1273   real(dp),intent(in) :: mditemp
1274   real(dp),intent(out) :: v2gauss
1275   real(dp),intent(in) :: amass(:)
1276   real(dp),intent(inout) :: vel(:,:)
1277 
1278   integer :: idim,natom
1279   real(dp) :: mass_tot,momentum,rescale_vel,vtest,sigma2
1280 
1281 ! *************************************************************************
1282 
1283    natom = size(amass)
1284 
1285   !--Get rid of center-of-mass velocity
1286    mass_tot = sum(amass(:))
1287    do idim=1,3
1288      momentum = sum(amass(:)*vel(idim,:))
1289      vel(idim,:) = vel(idim,:)-momentum/mass_tot
1290    end do
1291 
1292   !--Recompute v2gauss
1293    call vel_to_gauss(vel,amass,v2gauss,vtest)
1294 
1295   !--Now rescale the velocities to give the exact temperature
1296    rescale_vel = sqrt(three*natom*kb_HaK*mditemp/v2gauss)
1297    vel(:,:) = vel(:,:)*rescale_vel
1298 
1299   !--Recompute v2gauss
1300    call vel_to_gauss(vel,amass,v2gauss)
1301 
1302   !--Compute the variance and print
1303    sigma2 = (v2gauss/(3._dp*natom)-amass(1)*vtest**2)/kb_HaK
1304 
1305    write(message, '(a,d12.5,a,D12.5)' )&
1306 &   ' --- LOTF STEP : Effective temperature',&
1307 &   v2gauss/(3*natom*kb_HaK),' From variance', sigma2
1308    call wrtout(ab_out,message,'COLL')
1309    call wrtout(std_out,message,'COLL')
1310 
1311 
1312  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

SOURCE

1213  subroutine vel_to_gauss(vel_in,amass,v2gauss,vtest)
1214 
1215   implicit none
1216 
1217   !Arguments ------------------------------------
1218   real(dp),intent(out) :: v2gauss
1219   real(dp),intent(out),optional :: vtest
1220   real(dp),intent(in) :: vel_in(:,:)
1221   real(dp),intent(in) :: amass(:)
1222   !Local ---------------------------
1223   integer :: iatom,idim
1224 
1225 ! *************************************************************************
1226 
1227    if(present(vtest)) then
1228      v2gauss = zero
1229      vtest = zero
1230      do iatom=1,size(amass)
1231        do idim=1,3
1232          v2gauss = v2gauss+vel_in(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1233          vtest = vtest+vel_in(idim,iatom)
1234        end do
1235      end do
1236      vtest = vtest/(3._dp*size(amass))
1237    else
1238     !--Recompute v2gauss with the rescaled velocities
1239      v2gauss = zero
1240      do iatom=1,size(amass)
1241        do idim=1,3
1242          v2gauss = v2gauss+vel_in(idim,iatom)*vel_in(idim,iatom)*amass(iatom)
1243        end do
1244      end do
1245    end if
1246  end subroutine vel_to_gauss