TABLE OF CONTENTS


ABINIT/eval_lotf [ Modules ]

[ Top ] [ Modules ]

NAME

 eval_lotf

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

15 #if defined HAVE_CONFIG_H
16 #include "config.h"
17 #endif
18 
19 #include "abi_common.h"
20 
21 module eval_lotf
22 
23  use defs_basis
24  use m_errors
25 
26  implicit none
27  public  
28 
29  public ::                 &
30    eval_forces_u_n,        &
31    phi_n_calc,             &
32    calc_coord2,            &
33    eval_forces_u_n_2,      &
34    eval_force_devs_new_d,  &
35    upd_lis0,               &
36    tuneparms
37 
38 contains

eval_lotf/calc_coord2 [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 calc_coord2

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

194  subroutine calc_coord2(nneig,r0,rv,coordatom_dum)
195 
196   USE pbc_lotf,only : dist_pbc
197   USE GLUE_LOTF,only : calc_coord,rmrho
198 
199 !This section has been created automatically by the script Abilint (TD).
200 !Do not modify the following lines by hand.
201 #undef ABI_FUNC
202 #define ABI_FUNC 'calc_coord2'
203 !End of the abilint section
204 
205   integer,intent(in) :: nneig
206   real(dp),intent(out) :: coordatom_dum
207   real(dp),intent(in) :: r0(3),rv(3,nneig)
208 
209   !Local ---------------------------
210   integer :: ii
211   real(dp) :: r_au(nneig),RDV(3,nneig),rcut2_rho
212 
213 ! *************************************************************************
214 
215   !--Cutoff radius for the density
216   rcut2_rho  = rmrho**2 
217 
218   !--Run over neighbours
219   do ii=1,nneig  
220     call dist_pbc(rv(:,ii),r0,r_au(ii),RDV(:,ii))
221 
222     !--Cutoff density
223     if (r_au(ii) < rcut2_rho) then
224       call calc_coord(r_au(ii),coordatom_dum)
225     end if
226   enddo
227  end subroutine calc_coord2

eval_lotf/eval_force_devs_new_d [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 eval_force_devs_new_d

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

392  subroutine eval_force_devs_new_d(alpha_dum,nneig,nlist,neig2,nlist2,&
393    r0,rv,rv2,up_list,upp_list,fact2,ffit,&
394    forc_dum2,rho_p_sum,dcost_dalpha)
395 
396   use pbc_lotf,only : dist_pbc
397   use GLUE_LOTF,only : rmrho,dphi,rho_devs
398   use defs_param_lotf,only : lotfvar
399   use bond_lotf,only : nbondex,nfitmax,tafit,imat,ibmat_large
400 
401 !This section has been created automatically by the script Abilint (TD).
402 !Do not modify the following lines by hand.
403 #undef ABI_FUNC
404 #define ABI_FUNC 'eval_force_devs_new_d'
405 !End of the abilint section
406 
407   implicit none
408 
409   !Arguments ------------------------
410   real(dp),intent(in) ::  r0(3)
411   real(dp) ::  alpha_dum(3,nbondex)
412   integer ::  nneig,nlist(0:lotfvar%nneigx)
413   integer ::  neig2(lotfvar%nneigx),nlist2(lotfvar%nneigx,lotfvar%nneigx)
414   real(dp) ::  rv(3,lotfvar%nneigx),rv2(3,lotfvar%nneigx,lotfvar%nneigx)
415   real(dp) ::  up_list(lotfvar%natom),upp_list(lotfvar%natom)
416   real(dp) ::  fact2(nfitmax),ffit(3,nfitmax)
417   real(dp) ::  forc_dum2(3,0:nfitmax)
418   real(dp) ::  dcost_dalpha(3,nbondex)
419   !Local ---------------------------
420   integer :: iat,ii,jat,i0,nbond,ixyz
421   real(dp) :: r_au,RDV(3),r_st,rcut2_rho,U_tot0
422   real(dp) :: rho_neigh_d,rho_neigh_p,rho_neigh_pd
423   real(dp) :: rho_p_sum(3,nfitmax) 
424   real(dp) :: dFi_dbond_ij(3)
425   real(dp) :: r_au_ik,r_st_ik,RDV_ik(3)
426   real(dp) :: r_au_jk,r_st_jk,RDV_jk(3)
427   real(dp) :: rcut2_rho_ik
428   integer :: nbond_ik
429   real(dp) :: rho_d_ik,rho_p_ik,rho_pd_ik
430   real(dp) :: rcut2_rho_jk
431   integer :: nbond_jk
432   real(dp) :: rho_d_jk,rho_p_jk,rho_pd_jk
433   real(dp):: dFi_dbond_jk(3)
434   integer :: i2,i3,i4,kat
435   integer :: iunique(lotfvar%nneigx*lotfvar%nneigx),indunique
436 
437 ! *************************************************************************
438 
439   U_tot0 = zero
440   dFi_dbond_ij(:) = zero
441   dFi_dbond_jk(:) = zero
442   iat = nlist(0)
443 
444   indunique = 0
445   iunique(:) = 0
446 
447   !--1st NEIGHBOURS
448   first_neighbours : do ii=1,nneig 
449     jat=nlist(ii)
450 
451     call dist_pbc(rv(:,ii),r0,r_au,RDV)
452     r_st = sqrt(r_au)
453     i0 = imat(iat)
454     i2 = imat(jat)
455 
456     rho_neigh_d = zero
457     rho_neigh_p = zero
458     rho_neigh_pd = zero
459 
460     !--Fit
461     if (tafit(jat)) then  
462 
463       nbond = ibmat_large(i2,i0)
464       rcut2_rho  = (rmrho+alpha_dum(1,nbond)-dphi)**2 
465 
466       !--Cutoff 
467       if (r_au < rcut2_rho) then 
468         call rho_devs(r_au,alpha_dum(1,nbond),rho_neigh_d,&
469           rho_neigh_p,rho_neigh_pd)
470 
471         U_tot0 = up_list(iat) + up_list(jat)
472 
473         dFi_dbond_ij(:) = upp_list(iat)*rho_neigh_d*rho_p_sum(:,i0) + &
474           upp_list(jat)*rho_neigh_d*rho_neigh_p*RDV(:)/r_st+&
475           U_tot0*rho_neigh_pd*RDV(:)/r_st
476 
477         do ixyz=1,3
478           !--Atom iat
479           dcost_dalpha(1,nbond) = dcost_dalpha(1,nbond) - &
480             fact2(i0) * (forc_dum2(ixyz,i0) - ffit(ixyz,i0)) *&
481             dFi_dbond_ij(ixyz)
482         enddo
483 
484       endif
485     endif !--jat FIT atom
486 
487     !--Debug
488     go to 1002
489 
490     !--2nd NEIGHBOURS
491     second_neighbours : do i3=1,neig2(ii) 
492 
493       kat=nlist2(i3,ii)
494 
495       ! Make sure that kat is 'NEW' (not in the neigbours of iat neither
496       ! in previously explored neigbours lists...)
497 
498       if (tafit(jat).AND.tafit(kat).AND.kat /= iat) then
499         call dist_pbc(rv2(:,i3,ii),r0,r_au_ik,RDV_ik)
500         r_st_ik = sqrt(r_au_ik)
501         call dist_pbc(rv2(:,i3,ii),rv(:,ii),r_au_jk,RDV_jk)
502         r_st_jk = sqrt(r_au_jk)
503 
504         i4 = imat(kat)
505 
506         rho_d_ik = zero
507         rho_p_ik = zero
508         rho_pd_ik = zero
509 
510         !--FIT pair iat-kat 
511         if (tafit(kat)) then 
512 
513           nbond_ik = ibmat_large(i4,i0)
514           rcut2_rho_ik  = (rmrho+alpha_dum(1,nbond_ik)-dphi)**2
515           !--Cutoff 
516           if (r_au_ik < rcut2_rho_ik) then
517             call rho_devs(r_au_ik,alpha_dum(1,nbond_ik),rho_d_ik,rho_p_ik,rho_pd_ik)
518           endif ! cutoffs
519 
520         endif
521 
522         rho_d_jk = zero
523         rho_p_jk = zero
524         rho_pd_jk = zero
525 
526         !--FIT pair jat-kat 
527         if (tafit(kat).AND.tafit(jat)) then 
528 
529           nbond_jk = ibmat_large(i4,i0)
530           rcut2_rho_jk  = (rmrho+alpha_dum(1,nbond_jk)-dphi)**2
531 
532           !--Cutoff + no multiple counting
533           if (r_au_jk < rcut2_rho_jk.AND.kat > iat) then 
534             call rho_devs(r_au_jk,alpha_dum(1,nbond_jk),rho_d_jk,&
535               rho_p_jk,rho_pd_jk)
536           endif
537         endif ! FIT pair jat-kat
538 
539         if(r_au_jk < rcut2_rho_jk) then
540           if(r_au_ik < rcut2_rho_ik) then
541             dFi_dbond_jk(:) = dFi_dbond_jk(:) + upp_list(kat) * rho_d_jk * rho_p_ik *&
542               RDV_ik(:)/r_st_ik 
543           endif ! ik cutoff
544 
545           if(r_au < rcut2_rho) then
546             dFi_dbond_jk(:) = dFi_dbond_jk(:) + upp_list(jat) * rho_d_jk * rho_neigh_p *&
547               RDV(:)/r_st 
548           endif ! ij cutoff
549         endif ! jk cutoff 
550 
551         !--Cutoff 
552         if (r_au_jk < rcut2_rho_jk) then
553           do ixyz=1,3
554             !--Atom iat
555             dcost_dalpha(1,nbond_jk) = dcost_dalpha(1,nbond_jk) + &
556               fact2(i0) * (forc_dum2(ixyz,i0) - ffit(ixyz,i0)) *&
557               dFi_dbond_jk(ixyz)
558           enddo
559         endif
560 
561       endif ! atoms jat and kat are FIT and kat is not iat again 
562 
563     enddo second_neighbours
564 
565 1002 continue
566 
567   enddo first_neighbours
568 
569  end subroutine eval_force_devs_new_d

eval_lotf/eval_forces_U_n [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 eval_forces_U_n

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

244  subroutine eval_forces_U_n(nneig,nlist,r0,rv,up_list,forc_dum2)
245   use pbc_lotf,only : dist_pbc
246   use defs_param_lotf,only : lotfvar
247   use GLUE_LOTF, only :  rmrho,calc_rhop
248 
249 !This section has been created automatically by the script Abilint (TD).
250 !Do not modify the following lines by hand.
251 #undef ABI_FUNC
252 #define ABI_FUNC 'eval_forces_U_n'
253 !End of the abilint section
254 
255   integer,intent(in) ::  nneig,nlist(0:lotfvar%nneigx)
256   real(dp),intent(in) :: r0(3),rv(3,nneig),up_list(lotfvar%natom)
257   real(dp),intent(out):: forc_dum2(3)
258 
259   !Local ---------------------------
260   integer :: jat,iat
261   real(dp) :: U_tot,r_au,RDV(3),r_st,rcut2_rho,rhop_dum
262 
263 ! *************************************************************************
264 
265   !--Cutoffradius for the density
266   rcut2_rho  = rmrho**2 
267   U_tot = zero
268   forc_dum2(:) = zero
269 
270   iat = nlist(0)
271   do jat=1,nneig
272     call dist_pbc(rv(:,jat),r0,r_au,RDV)
273     if(r_au < rcut2_rho) then
274 
275       U_tot = up_list(iat) + up_list(nlist(jat))
276       r_st = sqrt(r_au)
277       call calc_rhop(r_st,rhop_dum)
278 
279       U_tot = U_tot * rhop_dum
280       forc_dum2(:) = forc_dum2(:) + U_tot * RDV(:) / r_st
281     endif
282   enddo
283  end subroutine eval_forces_U_n

eval_lotf/eval_forces_U_n_2 [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 eval_forces_U_n_2

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

301  subroutine eval_forces_U_n_2(alpha_dum,nneig,nlist,&
302    r0,rv,up_list,rho_p_sum,forc_dum2)
303 
304   use pbc_lotf,only : dist_pbc
305   use GLUE_LOTF,only : rmrho,dphi,rhop_value
306   use defs_param_lotf,only : lotfvar
307   use bond_lotf,only : nbondex,tafit,imat,ibmat_large
308 
309   ! Input/Output variables
310 
311 !This section has been created automatically by the script Abilint (TD).
312 !Do not modify the following lines by hand.
313 #undef ABI_FUNC
314 #define ABI_FUNC 'eval_forces_U_n_2'
315 !End of the abilint section
316 
317   real(dp),intent(in) ::  alpha_dum(3,nbondex)
318   integer,intent(in) ::  nneig,nlist(0:lotfvar%nneigx)
319   real(dp),intent(in) ::  r0(3),rv(3,nneig),up_list(lotfvar%natom)
320   real(dp),intent(out) ::  rho_p_sum(3),forc_dum2(3)
321 
322   !Local ---------------------------
323   integer :: iat,ii,jat,i0,nbond,i2
324   real(dp) :: U_tot,r_au,RDV(3),r_st,rcut2_rho
325   real(dp) :: rho_p_esc
326 
327 ! *************************************************************************
328 
329   U_tot = zero
330   rho_p_sum(:) = zero
331 
332   iat = nlist(0)
333 
334   !--(0) RHOP_SUM(iat) and FORC_DUM2(iat) (both need to be summed up before calculating FORCE DERIVATIVES)
335   do ii=1,nneig
336 
337     jat=nlist(ii)
338     call dist_pbc(rv(:,ii),r0,r_au,RDV)
339     r_st = sqrt(r_au)
340 
341     if (tafit(jat)) then
342 
343       !--(0a) RHOP_SUM
344       i0 = imat(iat) 
345       i2 = imat(jat)
346       nbond = ibmat_large(i2,i0)
347 
348       !--Cutoff radius for the density
349       rcut2_rho  = (rmrho+alpha_dum(1,nbond)-dphi)**2 
350       if(r_au < rcut2_rho) then
351         call rhop_value(r_st,alpha_dum(1,nbond),rho_p_esc)
352       end if
353     else
354 
355       rcut2_rho = rmrho**2
356       if(r_au < rcut2_rho) then
357         call rhop_value(r_st,dphi,rho_p_esc)
358       end if
359     endif
360 
361     if (r_au < rcut2_rho) then 
362 
363       rho_p_sum(:) = rho_p_sum(:) + rho_p_esc * RDV(:) / r_st 
364 
365       !--(0b) FORCES 
366       U_tot = up_list(iat) + up_list(jat)
367       U_tot = U_tot * rho_p_esc
368 
369       !--2-body + many-body
370       forc_dum2(:) = forc_dum2(:) + U_tot * RDV(:) / r_st
371     endif
372   enddo
373  end subroutine eval_forces_U_n_2

eval_lotf/phi_n_calc [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 phi_n_calc

FUNCTION

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

 55  !--Similar to SWCALC but without the triplets (used for the pair potential part and for the coordination) 
 56  subroutine phi_n_calc(alpha_dum,nneig,nlist,r0,rv,epot_dum,&
 57    &                        forcv,coordatom_dum,alpha_fdum)
 58 
 59   use defs_param_lotf,only : lotfvar
 60   use bond_lotf,only : nbondex,tafit,imat,ibmat,ibmat_large
 61   use glue_lotf,only : dphi,rcphi,rmrho,glue_pair,glue_pair_devs,calc_coord,calc_coord_new_d
 62   use pbc_lotf,only : dist_pbc
 63 
 64   ! INPUT/OUTPUT
 65 
 66 !This section has been created automatically by the script Abilint (TD).
 67 !Do not modify the following lines by hand.
 68 #undef ABI_FUNC
 69 #define ABI_FUNC 'phi_n_calc'
 70 !End of the abilint section
 71 
 72   real(dp),intent(in) :: alpha_dum(3,nbondex)
 73   integer ::  nneig,  iat, jat
 74   integer ::  nlist(0:nneig)
 75   real(dp),intent(in) :: r0(3)
 76   real(dp) ::  forcv(3,0:nneig)
 77   real(dp) ::  RDV(3,nneig)
 78   real(dp) ::  r_au(nneig)
 79   real(dp) ::  rv(3,nneig)
 80   real(dp) ::  coordatom_dum
 81   real(dp) ::  alpha_fdum(3,2,nbondex,1)
 82 
 83   !Local ---------------------------
 84   integer ::   i1,  i0, i3, nbond
 85   real(dp) ::  epot_dum, epot_2, drcphi
 86   real(dp) ::  rcut2_pair, rcut2_rho
 87   real(dp) ::  rref(3), fp(3), dfp(3,2)
 88 
 89 ! *************************************************************************
 90 
 91   drcphi = rcphi - dphi
 92   rref(:) = zero
 93   forcv(:,0:nneig) = zero
 94   epot_dum = zero
 95   coordatom_dum = zero
 96 
 97   iat  = nlist(0)
 98 
 99 
100   run_over_neighbours: do i1 = 1, nneig  ! Run over neighbours
101     !--PAIR: POTENTIAL, FORCES, DERIVATIVES
102     call dist_pbc(rv(:,i1),r0,r_au(i1),RDV(:,i1))
103 
104     jat = nlist(i1)
105     if (tafit(jat).and.(jat > iat)) then ! Only for fit pairs + NO doUBLE COUNTING
106 
107       i0 = imat(iat) ! Index of atom iat in the fitting indexing (1,nfit)
108       nbond = ibmat(i1,i0)
109 
110       !--Cutoff radius for the pair potential
111       rcut2_pair = (alpha_dum(1,nbond) + drcphi)**2 
112 
113       if ((r_au(i1) < rcut2_pair)) then 
114 
115         call glue_pair_devs(alpha_dum(:,nbond),RDV(:,i1),r_au(i1),epot_2,fp,dfp)
116 
117         epot_dum = epot_dum + epot_2
118 
119         forcv(:,0)  = forcv(:,0)  + fp(:)
120         forcv(:,i1) = forcv(:,i1) - fp(:)
121 
122         alpha_fdum(:,1,nbond,1) =  dfp(:,1)
123         alpha_fdum(:,2,nbond,1) =  dfp(:,2)
124       endif
125 
126     else !--We have to count the forces on the fit atoms by the non-fit neighbours
127 
128       rcut2_pair = rcphi**2
129       if ( (r_au(i1) < rcut2_pair).and.(jat > iat) ) then 
130         !--Cutoff pair potential + NO MULTIPLE COUNTING
131         call glue_pair(RDV(1,i1),r_au(i1),epot_2,fp)
132         epot_dum = epot_dum + epot_2
133 
134         forcv(:,0)  = forcv(:,0)  + fp(:)
135         forcv(:,i1) = forcv(:,i1) - fp(:)
136       endif
137     endif ! Pair potential and pair forces and derivs only for fit pairs
138 
139     !--COORDINATION 
140     select case(lotfvar%classic)
141     case(5)
142       !--Cutoff radius for the density 
143       rcut2_rho  = rmrho**2 
144 
145       !--Cutoff density
146       if (r_au(i1) < rcut2_rho) then
147         call calc_coord(r_au(i1),coordatom_dum)
148       end if
149 
150     case(6) 
151       if (tafit(jat)) then
152 
153         i0 = imat(iat) ! Index of atom iat in the fitting indexing (1,nfit)
154         i3 = imat(jat)
155         nbond = ibmat_large(i3,i0)
156 
157         !--Cutoff radius for the density
158         rcut2_rho  = (rmrho+alpha_dum(1,nbond)-dphi)**2 
159 
160         !--Cutoff density
161         if (r_au(i1) < rcut2_rho) then
162           call calc_coord_new_d(r_au(i1),alpha_dum(1,nbond),coordatom_dum)
163         end if
164 
165       else ! tafit = .false
166         !--Cutoff radius for the density
167         rcut2_rho  = rmrho**2 
168         if (r_au(i1) < rcut2_rho) then
169           call calc_coord(r_au(i1),coordatom_dum)
170         end if
171 
172       endif ! Fit/NonFit
173     end select
174   enddo run_over_neighbours
175 
176  end subroutine phi_n_calc

eval_lotf/tuneparms [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 tuneparms

FUNCTION

INPUTS

  tau0(3,natom)=atomic positions 
  rcf2_int=square of the cut off radius for this interaction 

OUTPUT

  tfit_int(3,nbondex)=logical that determines which parameters
  will be optimized ( if true the parameter is optimized)  

NOTES

  simple SW version (eventually modify for Tersoff or others             
  all the bonds considered here are already 100% in the fitting zone.    
  Here we just need to eliminate bonds or triplets that are too long..

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

843  subroutine tuneparms(tau0,tfit_int,rcf2_int)
844 
845   use defs_param_lotf,only : lotfvar
846   use bond_lotf,only : ibn_tot,ibnd_mat,nbondex
847   USE pbc_lotf,only : dist_pbc
848 
849 !This section has been created automatically by the script Abilint (TD).
850 !Do not modify the following lines by hand.
851 #undef ABI_FUNC
852 #define ABI_FUNC 'tuneparms'
853 !End of the abilint section
854 
855   implicit none
856 
857   !Arguments ------------------------
858   real(dp),intent(in) :: rcf2_int
859   real(dp),intent(in) :: tau0(3,lotfvar%natom)
860   logical,intent(out) :: tfit_int(3,nbondex)
861   !Local variables ------------------------------
862   integer :: il_fit_int,ibn,iat,i,il_fit_par,ntot
863   real(dp) :: r,r12,R1(3)
864   character(len=500) :: message
865   integer :: nbnds(lotfvar%natom),npract(lotfvar%natom)
866 
867 ! *************************************************************************
868 
869   tfit_int = .false.
870   nbnds(:) = 0
871   npract(:) = 0
872 
873   il_fit_int = 0
874   il_fit_par = 0
875   do ibn = 1,ibn_tot
876     iat = ibnd_mat(1,ibn)
877     i   = ibnd_mat(2,ibn)
878     call dist_pbc(tau0(:,i),tau0(:,iat),r,R1)
879     r12 = sqrt(r) 
880     if(r < rcf2_int) then
881       il_fit_int = il_fit_int + 1
882       nbnds(iat) = nbnds(iat) + 1 
883       nbnds(i)   = nbnds(i) + 1 
884 
885       tfit_int(:2,ibn) = .true.
886       il_fit_par = il_fit_par + 2
887       npract(iat) = npract(iat) + 2
888       npract(i)   = npract(i) + 2
889     endif
890   enddo
891 
892   ntot = sum(npract)
893 
894   do ibn = 1,ibn_tot
895     iat = ibnd_mat(1,ibn)
896     i   = ibnd_mat(2,ibn)
897     if( (npract(iat) > 50).and.(npract(i) > 50) )then  
898       MSG_ERROR('LOTF: NPRACT 50 (A) ') 
899     endif
900   enddo
901 
902   write(message,'(2(a,i8,a))')&
903     &    ' Total Active Bonds:          ',il_fit_int,ch10,&
904     &    ' Total Active Bondparms:      ',il_fit_par,ch10
905   call wrtout(std_out,message,'COLL')
906 
907 
908  end subroutine tuneparms
909 
910 end module eval_lotf

eval_lotf/upd_lis0 [ Functions ]

[ Top ] [ eval_lotf ] [ Functions ]

NAME

 upd_lis0

FUNCTION

  update neigbours relations

 lotfvar%nneigx : max number of neighbours
 tau0(3,natom) : atomic positions 
 neighl(lotfvar%nneigx,natom) : list of neighbours 
 nneig(natom) : number of neighbours 
 niter  : iteration number (itime) 

 upgraded to use the linked cell method (Knuth) 

INPUTS

PARENTS

      m_lotf

CHILDREN

      dist_pbc,wrtout

SOURCE

595  subroutine  upd_lis0(tau0,neighl,nneig,niter)   
596   use defs_param_lotf,only : lotfvar
597   use work_var_lotf,only : rcut_nbl
598   use tools_lotf,only : icf
599   use pbc_lotf,only : dist_pbc,r2,rd,pbc_bb_proj,pbc_bb_contract,pbc_aa_contract
600 
601 !This section has been created automatically by the script Abilint (TD).
602 !Do not modify the following lines by hand.
603 #undef ABI_FUNC
604 #define ABI_FUNC 'upd_lis0'
605 !End of the abilint section
606 
607   implicit none
608 
609   !Arguments ------------------------
610   integer,intent(in)   :: niter  
611   integer,intent(out)  :: nneig(lotfvar%natom), neighl(lotfvar%nneigx,lotfvar%natom)
612   real(dp),intent(inout) :: tau0(3,lotfvar%natom)
613   !Local --------------------------- 
614   integer,parameter :: icellx_a=9500
615   integer  :: icell_tot,icell
616   integer  :: icnr, icnr2, icn, icn2
617   integer  :: ic
618   integer  :: ix, iy, iz
619   integer  :: n_count, iat,i 
620   integer  :: ica(3) 
621   integer  :: head(icellx_a)  
622   integer  :: list(lotfvar%natom)  
623   integer  :: neigh_cel(icellx_a,13) 
624   real(dp) :: r3(3), rcut2 
625   real(dp) :: rdum(3)
626   real(dp) :: raa(3) 
627   character(len=500) :: message
628 
629 ! *************************************************************************
630 
631   n_count = niter - lotfvar%n0
632   write(message,'(a,4i8)')&
633     &  'Checking Input',n_count,niter,mod(n_count,lotfvar%nitex),lotfvar%nitex
634   call wrtout(std_out,message,'COLL')
635 
636   rcut2   = rcut_nbl*rcut_nbl     
637 
638   !--(1) gets all atoms within the "same" repeated cell zone
639   !     (from -0.5 to 0.5 in relative units) 
640   !at this part is changed to take into account symetries other than orthorombic.
641   r3(:) = zero
642 
643   !--length of the cell vectors :   
644   raa(:) = pbc_aa_contract()
645 
646   !--put tau0 and taum in the cell -0.5/+0.5 for tau0
647   if(lotfvar%version==1) then 
648     do i =1, lotfvar%natom
649       !--compute rd
650       call dist_pbc(tau0(:,i),r3)
651       rdum(:)  =  rd(:) - tau0(:,i)
652       tau0(:,i) =  tau0(:,i) + rdum(:)
653       ! taum(:,i) =  taum(:,i) + rdum(:)
654     enddo
655   elseif(lotfvar%version==2) then  
656     do i =1, lotfvar%natom
657       !--compute rd
658       call dist_pbc(tau0(:,i),r3)
659       rdum(:)  =  rd(:) - tau0(:,i)
660       tau0(:,i) =  tau0(:,i) + rdum(:)
661     enddo
662   endif ! lotfvar%version says if taum is touched
663 
664   !--clears neighbour lists
665   nneig(:) = 0 
666   neighl(:,:) = 0
667 
668 
669   !--OK, NOW KNUTH TRICK 
670   ! determines the cell partition of the run 
671   ! it means we cut volumes with edges // to the cell vectors   
672   ! but with length divided by an integer ica(1:3)
673 
674   ica(:) = int(raa(:)/sqrt(4*rcut2)) 
675 
676 
677   write(message,'(a,3f12.8,4a,2f12.8)')&
678     &     ' raa = ', raa(:),ch10,&
679     &     ' Neighbour List Cutoff: ',ch10,&
680     &     '   rcut2, sqrt(4*rcut2) ',rcut2,sqrt(4*rcut2) 
681   call wrtout(std_out,message,'COLL')
682 
683   where(ica < 1) ica = 1
684   icell_tot =  product(ica)
685 
686   write(message,'(3a,2i8,2a,3i8)')&
687     & 'UPDLIS0 : SYSTEM SUBCELL DISTRIBUTION: ',ch10,&
688     & 'TOT. NO. OF CELLS (& max.) : ',icell_tot, icellx_a,ch10,&
689     & 'ICX,ICY,ICZ = ',ica(1),ica(2),ica(3)
690   call wrtout(std_out,message,'COLL')
691 
692 
693   if(icell_tot > icellx_a) then 
694     write(message,'(a,i8,2a)')&
695       & 'ERROR IN UPD_LIS0: ICELL_TOT = ',icell_tot,ch10,&
696       & 'ERROR IN UPD_LIS0: RAISE ICELLX_A'
697     MSG_ERROR(message)
698   endif
699 
700   !-- clears head vector & constructs head & list      
701   head(:icellx_a) = 0
702 
703   do i=1,lotfvar%natom
704     rdum = pbc_bb_proj(tau0(:,i))
705     icell = 1 + mod(int( (rdum(1)+0.5)* float(ica(1)) ),ica(1)) &
706       + mod(int( (rdum(2)+0.5)* float(ica(2)) ),ica(2)) * ica(1) & 
707       + mod(int( (rdum(3)+0.5)* float(ica(3)) ),ica(3)) * ica(1) * ica(2)
708 
709     list(i) = head(icell) 
710     head(icell) = i 
711   enddo
712 
713   !--Constructs the 13 neighbours of each cell 
714   do  IZ = 1, ica(3)
715     do  IY = 1, ica(2)
716       do  IX = 1, ica(1)
717         ic               = icf(ix  ,iy,iz    ,ica(1),ica(2),ica(3))
718         neigh_cel(ic,1)  = icf(ix+1,iy,iz    ,ica(1),ica(2),ica(3)) 
719         neigh_cel(ic,2)  = icf(ix+1,iy+1,iz  ,ica(1),ica(2),ica(3)) 
720         neigh_cel(ic,3)  = icf(ix  ,iy+1,iz  ,ica(1),ica(2),ica(3)) 
721         neigh_cel(ic,4)  = icf(ix-1,iy+1,iz  ,ica(1),ica(2),ica(3)) 
722         neigh_cel(ic,5)  = icf(ix+1,iy  ,iz-1,ica(1),ica(2),ica(3)) 
723         neigh_cel(ic,6)  = icf(ix+1,iy+1,iz-1,ica(1),ica(2),ica(3)) 
724         neigh_cel(ic,7)  = icf(ix  ,iy+1,iz-1,ica(1),ica(2),ica(3)) 
725         neigh_cel(ic,8)  = icf(ix-1,iy+1,iz-1,ica(1),ica(2),ica(3)) 
726         neigh_cel(ic,9)  = icf(ix+1,iy  ,iz+1,ica(1),ica(2),ica(3)) 
727         neigh_cel(ic,10) = icf(ix+1,iy+1,iz+1,ica(1),ica(2),ica(3)) 
728         neigh_cel(ic,11) = icf(ix  ,iy+1,iz+1,ica(1),ica(2),ica(3)) 
729         neigh_cel(ic,12) = icf(ix-1,iy+1,iz+1,ica(1),ica(2),ica(3)) 
730         neigh_cel(ic,13) = icf(ix  ,iy  ,iz+1,ica(1),ica(2),ica(3)) 
731       enddo
732     enddo
733   enddo
734 
735   !--Safety loops to avoid repetitions
736 
737   !--(1) to avoid having twice the same neigh. cell 
738   do icell = 1,icell_tot
739     do icnr = 1,13
740       icn =  neigh_cel(icell,icnr)
741       do icnr2 = icnr+1,13
742         icn2 =  neigh_cel(icell,icnr2)
743         if(icn2==icn)  neigh_cel(icell,icnr2) = icell 
744       enddo
745     enddo
746   enddo
747 
748   !--(2) to avoid counting twice the interaction between the same 
749   !     couple of  neigh. cells 
750   do icell = 1,icell_tot
751     do icnr = 1,13
752       icn =  neigh_cel(icell,icnr)
753       do icnr2 = 1,13
754         icn2 =  neigh_cel(icn,icnr2)
755         if(icn2==icell) neigh_cel(icn,icnr2) = icn 
756       enddo
757     enddo
758   enddo
759 
760   !--(3) constructs neighbour list looking through neighbour cells only
761   do icell = 1,icell_tot
762     iat = head(icell) 
763     do while(iat > 0)
764       i = list(iat) 
765       do while(i  > 0) 
766         if(i==iat) MSG_ERROR("i==iat")
767         call dist_pbc(tau0(:,i),tau0(:,iat))
768         if(r2 < rcut2) then 
769           nneig(iat) = nneig(iat) + 1
770           nneig(i)   = nneig(i)   + 1
771           neighl(nneig(iat),iat) = i 
772           neighl(nneig(i)  ,i  ) = iat 
773         endif  ! distance check 
774         i = list(i) 
775       enddo
776 
777       if(ANY(nneig(:) > lotfvar%nneigx))  then 
778         write(message,'(a,i8,a)')&
779           'UPD_LIS CLASSIC: max no. of neighbours: ',lotfvar%nneigx,' is too small'
780         MSG_ERROR(message)
781       endif
782 
783 
784       !   mask to avoid self interaction within the same cell
785       !   considered as a neighbor and with different cells more than once
786       do icnr = 1,13
787         icn = neigh_cel(icell,icnr)
788         if(icn==icell) cycle
789         i = head(icn) 
790 
791         do while(i>0)
792           if(i==iat) MSG_ERROR("i==iat")
793           call dist_pbc(tau0(:,i),tau0(:,iat))
794           if(r2 < rcut2) then 
795             nneig(iat) = nneig(iat) + 1
796             nneig(i)   = nneig(i)   + 1
797             neighl(nneig(iat),iat) = i 
798             neighl(nneig(i)  ,i  ) = iat 
799           endif ! distance check  
800           i = list(i) 
801         enddo
802         if(ANY(nneig(:) > lotfvar%nneigx))  then 
803           write(message,'(a,i8,a)')&
804             'UPD_LIS CLASSIC: max no. of neighbours: ',lotfvar%nneigx,' is too small'
805           MSG_ERROR(message)
806         endif
807       enddo
808 
809       iat = list(iat) 
810     enddo
811   enddo
812  end subroutine upd_lis0