TABLE OF CONTENTS
- ABINIT/eval_lotf
- eval_lotf/calc_coord2
- eval_lotf/eval_force_devs_new_d
- eval_lotf/eval_forces_U_n
- eval_lotf/eval_forces_U_n_2
- eval_lotf/phi_n_calc
- eval_lotf/tuneparms
- eval_lotf/upd_lis0
ABINIT/eval_lotf [ Modules ]
NAME
eval_lotf
FUNCTION
COPYRIGHT
Copyright (C) 2005-2024 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
14 #if defined HAVE_CONFIG_H 15 #include "config.h" 16 #endif 17 18 #include "abi_common.h" 19 20 module eval_lotf 21 22 use defs_basis 23 use m_errors 24 25 implicit none 26 public 27 28 public :: & 29 eval_forces_u_n, & 30 phi_n_calc, & 31 calc_coord2, & 32 eval_forces_u_n_2, & 33 eval_force_devs_new_d, & 34 upd_lis0, & 35 tuneparms 36 37 contains
eval_lotf/calc_coord2 [ Functions ]
[ Top ] [ eval_lotf ] [ Functions ]
NAME
calc_coord2
FUNCTION
INPUTS
SOURCE
174 subroutine calc_coord2(nneig,r0,rv,coordatom_dum) 175 176 USE pbc_lotf,only : dist_pbc 177 USE GLUE_LOTF,only : calc_coord,rmrho 178 integer,intent(in) :: nneig 179 real(dp),intent(out) :: coordatom_dum 180 real(dp),intent(in) :: r0(3),rv(3,nneig) 181 182 !Local --------------------------- 183 integer :: ii 184 real(dp) :: r_au(nneig),RDV(3,nneig),rcut2_rho 185 186 ! ************************************************************************* 187 188 !--Cutoff radius for the density 189 rcut2_rho = rmrho**2 190 191 !--Run over neighbours 192 do ii=1,nneig 193 call dist_pbc(rv(:,ii),r0,r_au(ii),RDV(:,ii)) 194 195 !--Cutoff density 196 if (r_au(ii) < rcut2_rho) then 197 call calc_coord(r_au(ii),coordatom_dum) 198 end if 199 enddo 200 end subroutine calc_coord2
eval_lotf/eval_force_devs_new_d [ Functions ]
[ Top ] [ eval_lotf ] [ Functions ]
NAME
eval_force_devs_new_d
FUNCTION
INPUTS
SOURCE
333 subroutine eval_force_devs_new_d(alpha_dum,nneig,nlist,neig2,nlist2,& 334 r0,rv,rv2,up_list,upp_list,fact2,ffit,& 335 forc_dum2,rho_p_sum,dcost_dalpha) 336 337 use pbc_lotf,only : dist_pbc 338 use GLUE_LOTF,only : rmrho,dphi,rho_devs 339 use defs_param_lotf,only : lotfvar 340 use bond_lotf,only : nbondex,nfitmax,tafit,imat,ibmat_large 341 implicit none 342 343 !Arguments ------------------------ 344 real(dp),intent(in) :: r0(3) 345 real(dp) :: alpha_dum(3,nbondex) 346 integer :: nneig,nlist(0:lotfvar%nneigx) 347 integer :: neig2(lotfvar%nneigx),nlist2(lotfvar%nneigx,lotfvar%nneigx) 348 real(dp) :: rv(3,lotfvar%nneigx),rv2(3,lotfvar%nneigx,lotfvar%nneigx) 349 real(dp) :: up_list(lotfvar%natom),upp_list(lotfvar%natom) 350 real(dp) :: fact2(nfitmax),ffit(3,nfitmax) 351 real(dp) :: forc_dum2(3,0:nfitmax) 352 real(dp) :: dcost_dalpha(3,nbondex) 353 !Local --------------------------- 354 integer :: iat,ii,jat,i0,nbond,ixyz 355 real(dp) :: r_au,RDV(3),r_st,rcut2_rho,U_tot0 356 real(dp) :: rho_neigh_d,rho_neigh_p,rho_neigh_pd 357 real(dp) :: rho_p_sum(3,nfitmax) 358 real(dp) :: dFi_dbond_ij(3) 359 real(dp) :: r_au_ik,r_st_ik,RDV_ik(3) 360 real(dp) :: r_au_jk,r_st_jk,RDV_jk(3) 361 real(dp) :: rcut2_rho_ik 362 integer :: nbond_ik 363 real(dp) :: rho_d_ik,rho_p_ik,rho_pd_ik 364 real(dp) :: rcut2_rho_jk 365 integer :: nbond_jk 366 real(dp) :: rho_d_jk,rho_p_jk,rho_pd_jk 367 real(dp):: dFi_dbond_jk(3) 368 integer :: i2,i3,i4,kat 369 integer :: iunique(lotfvar%nneigx*lotfvar%nneigx),indunique 370 371 ! ************************************************************************* 372 373 U_tot0 = zero 374 dFi_dbond_ij(:) = zero 375 dFi_dbond_jk(:) = zero 376 iat = nlist(0) 377 378 indunique = 0 379 iunique(:) = 0 380 381 !--1st NEIGHBOURS 382 first_neighbours : do ii=1,nneig 383 jat=nlist(ii) 384 385 call dist_pbc(rv(:,ii),r0,r_au,RDV) 386 r_st = sqrt(r_au) 387 i0 = imat(iat) 388 i2 = imat(jat) 389 390 rho_neigh_d = zero 391 rho_neigh_p = zero 392 rho_neigh_pd = zero 393 394 !--Fit 395 if (tafit(jat)) then 396 397 nbond = ibmat_large(i2,i0) 398 rcut2_rho = (rmrho+alpha_dum(1,nbond)-dphi)**2 399 400 !--Cutoff 401 if (r_au < rcut2_rho) then 402 call rho_devs(r_au,alpha_dum(1,nbond),rho_neigh_d,& 403 rho_neigh_p,rho_neigh_pd) 404 405 U_tot0 = up_list(iat) + up_list(jat) 406 407 dFi_dbond_ij(:) = upp_list(iat)*rho_neigh_d*rho_p_sum(:,i0) + & 408 upp_list(jat)*rho_neigh_d*rho_neigh_p*RDV(:)/r_st+& 409 U_tot0*rho_neigh_pd*RDV(:)/r_st 410 411 do ixyz=1,3 412 !--Atom iat 413 dcost_dalpha(1,nbond) = dcost_dalpha(1,nbond) - & 414 fact2(i0) * (forc_dum2(ixyz,i0) - ffit(ixyz,i0)) *& 415 dFi_dbond_ij(ixyz) 416 enddo 417 418 endif 419 endif !--jat FIT atom 420 421 !--Debug 422 go to 1002 423 424 !--2nd NEIGHBOURS 425 second_neighbours : do i3=1,neig2(ii) 426 427 kat=nlist2(i3,ii) 428 429 ! Make sure that kat is 'NEW' (not in the neigbours of iat neither 430 ! in previously explored neigbours lists...) 431 432 if (tafit(jat).AND.tafit(kat).AND.kat /= iat) then 433 call dist_pbc(rv2(:,i3,ii),r0,r_au_ik,RDV_ik) 434 r_st_ik = sqrt(r_au_ik) 435 call dist_pbc(rv2(:,i3,ii),rv(:,ii),r_au_jk,RDV_jk) 436 r_st_jk = sqrt(r_au_jk) 437 438 i4 = imat(kat) 439 440 rho_d_ik = zero 441 rho_p_ik = zero 442 rho_pd_ik = zero 443 444 !--FIT pair iat-kat 445 if (tafit(kat)) then 446 447 nbond_ik = ibmat_large(i4,i0) 448 rcut2_rho_ik = (rmrho+alpha_dum(1,nbond_ik)-dphi)**2 449 !--Cutoff 450 if (r_au_ik < rcut2_rho_ik) then 451 call rho_devs(r_au_ik,alpha_dum(1,nbond_ik),rho_d_ik,rho_p_ik,rho_pd_ik) 452 endif ! cutoffs 453 454 endif 455 456 rho_d_jk = zero 457 rho_p_jk = zero 458 rho_pd_jk = zero 459 460 !--FIT pair jat-kat 461 if (tafit(kat).AND.tafit(jat)) then 462 463 nbond_jk = ibmat_large(i4,i0) 464 rcut2_rho_jk = (rmrho+alpha_dum(1,nbond_jk)-dphi)**2 465 466 !--Cutoff + no multiple counting 467 if (r_au_jk < rcut2_rho_jk.AND.kat > iat) then 468 call rho_devs(r_au_jk,alpha_dum(1,nbond_jk),rho_d_jk,& 469 rho_p_jk,rho_pd_jk) 470 endif 471 endif ! FIT pair jat-kat 472 473 if(r_au_jk < rcut2_rho_jk) then 474 if(r_au_ik < rcut2_rho_ik) then 475 dFi_dbond_jk(:) = dFi_dbond_jk(:) + upp_list(kat) * rho_d_jk * rho_p_ik *& 476 RDV_ik(:)/r_st_ik 477 endif ! ik cutoff 478 479 if(r_au < rcut2_rho) then 480 dFi_dbond_jk(:) = dFi_dbond_jk(:) + upp_list(jat) * rho_d_jk * rho_neigh_p *& 481 RDV(:)/r_st 482 endif ! ij cutoff 483 endif ! jk cutoff 484 485 !--Cutoff 486 if (r_au_jk < rcut2_rho_jk) then 487 do ixyz=1,3 488 !--Atom iat 489 dcost_dalpha(1,nbond_jk) = dcost_dalpha(1,nbond_jk) + & 490 fact2(i0) * (forc_dum2(ixyz,i0) - ffit(ixyz,i0)) *& 491 dFi_dbond_jk(ixyz) 492 enddo 493 endif 494 495 endif ! atoms jat and kat are FIT and kat is not iat again 496 497 enddo second_neighbours 498 499 1002 continue 500 501 enddo first_neighbours 502 503 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
SOURCE
211 subroutine eval_forces_U_n(nneig,nlist,r0,rv,up_list,forc_dum2) 212 use pbc_lotf,only : dist_pbc 213 use defs_param_lotf,only : lotfvar 214 use GLUE_LOTF, only : rmrho,calc_rhop 215 integer,intent(in) :: nneig,nlist(0:lotfvar%nneigx) 216 real(dp),intent(in) :: r0(3),rv(3,nneig),up_list(lotfvar%natom) 217 real(dp),intent(out):: forc_dum2(3) 218 219 !Local --------------------------- 220 integer :: jat,iat 221 real(dp) :: U_tot,r_au,RDV(3),r_st,rcut2_rho,rhop_dum 222 223 ! ************************************************************************* 224 225 !--Cutoffradius for the density 226 rcut2_rho = rmrho**2 227 U_tot = zero 228 forc_dum2(:) = zero 229 230 iat = nlist(0) 231 do jat=1,nneig 232 call dist_pbc(rv(:,jat),r0,r_au,RDV) 233 if(r_au < rcut2_rho) then 234 235 U_tot = up_list(iat) + up_list(nlist(jat)) 236 r_st = sqrt(r_au) 237 call calc_rhop(r_st,rhop_dum) 238 239 U_tot = U_tot * rhop_dum 240 forc_dum2(:) = forc_dum2(:) + U_tot * RDV(:) / r_st 241 endif 242 enddo 243 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
SOURCE
255 subroutine eval_forces_U_n_2(alpha_dum,nneig,nlist,& 256 r0,rv,up_list,rho_p_sum,forc_dum2) 257 258 use pbc_lotf,only : dist_pbc 259 use GLUE_LOTF,only : rmrho,dphi,rhop_value 260 use defs_param_lotf,only : lotfvar 261 use bond_lotf,only : nbondex,tafit,imat,ibmat_large 262 263 ! Input/Output variables 264 real(dp),intent(in) :: alpha_dum(3,nbondex) 265 integer,intent(in) :: nneig,nlist(0:lotfvar%nneigx) 266 real(dp),intent(in) :: r0(3),rv(3,nneig),up_list(lotfvar%natom) 267 real(dp),intent(out) :: rho_p_sum(3),forc_dum2(3) 268 269 !Local --------------------------- 270 integer :: iat,ii,jat,i0,nbond,i2 271 real(dp) :: U_tot,r_au,RDV(3),r_st,rcut2_rho 272 real(dp) :: rho_p_esc 273 274 ! ************************************************************************* 275 276 U_tot = zero 277 rho_p_sum(:) = zero 278 279 iat = nlist(0) 280 281 !--(0) RHOP_SUM(iat) and FORC_DUM2(iat) (both need to be summed up before calculating FORCE DERIVATIVES) 282 do ii=1,nneig 283 284 jat=nlist(ii) 285 call dist_pbc(rv(:,ii),r0,r_au,RDV) 286 r_st = sqrt(r_au) 287 288 if (tafit(jat)) then 289 290 !--(0a) RHOP_SUM 291 i0 = imat(iat) 292 i2 = imat(jat) 293 nbond = ibmat_large(i2,i0) 294 295 !--Cutoff radius for the density 296 rcut2_rho = (rmrho+alpha_dum(1,nbond)-dphi)**2 297 if(r_au < rcut2_rho) then 298 call rhop_value(r_st,alpha_dum(1,nbond),rho_p_esc) 299 end if 300 else 301 302 rcut2_rho = rmrho**2 303 if(r_au < rcut2_rho) then 304 call rhop_value(r_st,dphi,rho_p_esc) 305 end if 306 endif 307 308 if (r_au < rcut2_rho) then 309 310 rho_p_sum(:) = rho_p_sum(:) + rho_p_esc * RDV(:) / r_st 311 312 !--(0b) FORCES 313 U_tot = up_list(iat) + up_list(jat) 314 U_tot = U_tot * rho_p_esc 315 316 !--2-body + many-body 317 forc_dum2(:) = forc_dum2(:) + U_tot * RDV(:) / r_st 318 endif 319 enddo 320 end subroutine eval_forces_U_n_2
eval_lotf/phi_n_calc [ Functions ]
[ Top ] [ eval_lotf ] [ Functions ]
NAME
phi_n_calc
FUNCTION
INPUTS
SOURCE
48 !--Similar to SWCALC but without the triplets (used for the pair potential part and for the coordination) 49 subroutine phi_n_calc(alpha_dum,nneig,nlist,r0,rv,epot_dum,& 50 & forcv,coordatom_dum,alpha_fdum) 51 52 use defs_param_lotf,only : lotfvar 53 use bond_lotf,only : nbondex,tafit,imat,ibmat,ibmat_large 54 use glue_lotf,only : dphi,rcphi,rmrho,glue_pair,glue_pair_devs,calc_coord,calc_coord_new_d 55 use pbc_lotf,only : dist_pbc 56 57 ! INPUT/OUTPUT 58 real(dp),intent(in) :: alpha_dum(3,nbondex) 59 integer :: nneig, iat, jat 60 integer :: nlist(0:nneig) 61 real(dp),intent(in) :: r0(3) 62 real(dp) :: forcv(3,0:nneig) 63 real(dp) :: RDV(3,nneig) 64 real(dp) :: r_au(nneig) 65 real(dp) :: rv(3,nneig) 66 real(dp) :: coordatom_dum 67 real(dp) :: alpha_fdum(3,2,nbondex,1) 68 69 !Local --------------------------- 70 integer :: i1, i0, i3, nbond 71 real(dp) :: epot_dum, epot_2, drcphi 72 real(dp) :: rcut2_pair, rcut2_rho 73 real(dp) :: rref(3), fp(3), dfp(3,2) 74 75 ! ************************************************************************* 76 77 drcphi = rcphi - dphi 78 rref(:) = zero 79 forcv(:,0:nneig) = zero 80 epot_dum = zero 81 coordatom_dum = zero 82 83 iat = nlist(0) 84 85 86 run_over_neighbours: do i1 = 1, nneig ! Run over neighbours 87 !--PAIR: POTENTIAL, FORCES, DERIVATIVES 88 call dist_pbc(rv(:,i1),r0,r_au(i1),RDV(:,i1)) 89 90 jat = nlist(i1) 91 if (tafit(jat).and.(jat > iat)) then ! Only for fit pairs + NO doUBLE COUNTING 92 93 i0 = imat(iat) ! Index of atom iat in the fitting indexing (1,nfit) 94 nbond = ibmat(i1,i0) 95 96 !--Cutoff radius for the pair potential 97 rcut2_pair = (alpha_dum(1,nbond) + drcphi)**2 98 99 if ((r_au(i1) < rcut2_pair)) then 100 101 call glue_pair_devs(alpha_dum(:,nbond),RDV(:,i1),r_au(i1),epot_2,fp,dfp) 102 103 epot_dum = epot_dum + epot_2 104 105 forcv(:,0) = forcv(:,0) + fp(:) 106 forcv(:,i1) = forcv(:,i1) - fp(:) 107 108 alpha_fdum(:,1,nbond,1) = dfp(:,1) 109 alpha_fdum(:,2,nbond,1) = dfp(:,2) 110 endif 111 112 else !--We have to count the forces on the fit atoms by the non-fit neighbours 113 114 rcut2_pair = rcphi**2 115 if ( (r_au(i1) < rcut2_pair).and.(jat > iat) ) then 116 !--Cutoff pair potential + NO MULTIPLE COUNTING 117 call glue_pair(RDV(1,i1),r_au(i1),epot_2,fp) 118 epot_dum = epot_dum + epot_2 119 120 forcv(:,0) = forcv(:,0) + fp(:) 121 forcv(:,i1) = forcv(:,i1) - fp(:) 122 endif 123 endif ! Pair potential and pair forces and derivs only for fit pairs 124 125 !--COORDINATION 126 select case(lotfvar%classic) 127 case(5) 128 !--Cutoff radius for the density 129 rcut2_rho = rmrho**2 130 131 !--Cutoff density 132 if (r_au(i1) < rcut2_rho) then 133 call calc_coord(r_au(i1),coordatom_dum) 134 end if 135 136 case(6) 137 if (tafit(jat)) then 138 139 i0 = imat(iat) ! Index of atom iat in the fitting indexing (1,nfit) 140 i3 = imat(jat) 141 nbond = ibmat_large(i3,i0) 142 143 !--Cutoff radius for the density 144 rcut2_rho = (rmrho+alpha_dum(1,nbond)-dphi)**2 145 146 !--Cutoff density 147 if (r_au(i1) < rcut2_rho) then 148 call calc_coord_new_d(r_au(i1),alpha_dum(1,nbond),coordatom_dum) 149 end if 150 151 else ! tafit = .false 152 !--Cutoff radius for the density 153 rcut2_rho = rmrho**2 154 if (r_au(i1) < rcut2_rho) then 155 call calc_coord(r_au(i1),coordatom_dum) 156 end if 157 158 endif ! Fit/NonFit 159 end select 160 enddo run_over_neighbours 161 162 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..
SOURCE
758 subroutine tuneparms(tau0,tfit_int,rcf2_int) 759 760 use defs_param_lotf,only : lotfvar 761 use bond_lotf,only : ibn_tot,ibnd_mat,nbondex 762 USE pbc_lotf,only : dist_pbc 763 implicit none 764 765 !Arguments ------------------------ 766 real(dp),intent(in) :: rcf2_int 767 real(dp),intent(in) :: tau0(3,lotfvar%natom) 768 logical,intent(out) :: tfit_int(3,nbondex) 769 !Local variables ------------------------------ 770 integer :: il_fit_int,ibn,iat,i,il_fit_par,ntot 771 real(dp) :: r,r12,R1(3) 772 character(len=500) :: message 773 integer :: nbnds(lotfvar%natom),npract(lotfvar%natom) 774 775 ! ************************************************************************* 776 777 tfit_int = .false. 778 nbnds(:) = 0 779 npract(:) = 0 780 781 il_fit_int = 0 782 il_fit_par = 0 783 do ibn = 1,ibn_tot 784 iat = ibnd_mat(1,ibn) 785 i = ibnd_mat(2,ibn) 786 call dist_pbc(tau0(:,i),tau0(:,iat),r,R1) 787 r12 = sqrt(r) 788 if(r < rcf2_int) then 789 il_fit_int = il_fit_int + 1 790 nbnds(iat) = nbnds(iat) + 1 791 nbnds(i) = nbnds(i) + 1 792 793 tfit_int(:2,ibn) = .true. 794 il_fit_par = il_fit_par + 2 795 npract(iat) = npract(iat) + 2 796 npract(i) = npract(i) + 2 797 endif 798 enddo 799 800 ntot = sum(npract) 801 802 do ibn = 1,ibn_tot 803 iat = ibnd_mat(1,ibn) 804 i = ibnd_mat(2,ibn) 805 if( (npract(iat) > 50).and.(npract(i) > 50) )then 806 ABI_ERROR('LOTF: NPRACT 50 (A) ') 807 endif 808 enddo 809 810 write(message,'(2(a,i8,a))')& 811 & ' Total Active Bonds: ',il_fit_int,ch10,& 812 & ' Total Active Bondparms: ',il_fit_par,ch10 813 call wrtout(std_out,message,'COLL') 814 815 816 end subroutine tuneparms 817 818 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
SOURCE
523 subroutine upd_lis0(tau0,neighl,nneig,niter) 524 use defs_param_lotf,only : lotfvar 525 use work_var_lotf,only : rcut_nbl 526 use tools_lotf,only : icf 527 use pbc_lotf,only : dist_pbc,r2,rd,pbc_bb_proj,pbc_bb_contract,pbc_aa_contract 528 implicit none 529 530 !Arguments ------------------------ 531 integer,intent(in) :: niter 532 integer,intent(out) :: nneig(lotfvar%natom), neighl(lotfvar%nneigx,lotfvar%natom) 533 real(dp),intent(inout) :: tau0(3,lotfvar%natom) 534 !Local --------------------------- 535 integer,parameter :: icellx_a=9500 536 integer :: icell_tot,icell 537 integer :: icnr, icnr2, icn, icn2 538 integer :: ic 539 integer :: ix, iy, iz 540 integer :: n_count, iat,i 541 integer :: ica(3) 542 integer :: head(icellx_a) 543 integer :: list(lotfvar%natom) 544 integer :: neigh_cel(icellx_a,13) 545 real(dp) :: r3(3), rcut2 546 real(dp) :: rdum(3) 547 real(dp) :: raa(3) 548 character(len=500) :: message 549 550 ! ************************************************************************* 551 552 n_count = niter - lotfvar%n0 553 write(message,'(a,4i8)')& 554 & 'Checking Input',n_count,niter,mod(n_count,lotfvar%nitex),lotfvar%nitex 555 call wrtout(std_out,message,'COLL') 556 557 rcut2 = rcut_nbl*rcut_nbl 558 559 !--(1) gets all atoms within the "same" repeated cell zone 560 ! (from -0.5 to 0.5 in relative units) 561 !at this part is changed to take into account symetries other than orthorombic. 562 r3(:) = zero 563 564 !--length of the cell vectors : 565 raa(:) = pbc_aa_contract() 566 567 !--put tau0 and taum in the cell -0.5/+0.5 for tau0 568 if(lotfvar%version==1) then 569 do i =1, lotfvar%natom 570 !--compute rd 571 call dist_pbc(tau0(:,i),r3) 572 rdum(:) = rd(:) - tau0(:,i) 573 tau0(:,i) = tau0(:,i) + rdum(:) 574 ! taum(:,i) = taum(:,i) + rdum(:) 575 enddo 576 elseif(lotfvar%version==2) then 577 do i =1, lotfvar%natom 578 !--compute rd 579 call dist_pbc(tau0(:,i),r3) 580 rdum(:) = rd(:) - tau0(:,i) 581 tau0(:,i) = tau0(:,i) + rdum(:) 582 enddo 583 endif ! lotfvar%version says if taum is touched 584 585 !--clears neighbour lists 586 nneig(:) = 0 587 neighl(:,:) = 0 588 589 590 !--OK, NOW KNUTH TRICK 591 ! determines the cell partition of the run 592 ! it means we cut volumes with edges // to the cell vectors 593 ! but with length divided by an integer ica(1:3) 594 595 ica(:) = int(raa(:)/sqrt(4*rcut2)) 596 597 598 write(message,'(a,3f12.8,4a,2f12.8)')& 599 & ' raa = ', raa(:),ch10,& 600 & ' Neighbour List Cutoff: ',ch10,& 601 & ' rcut2, sqrt(4*rcut2) ',rcut2,sqrt(4*rcut2) 602 call wrtout(std_out,message,'COLL') 603 604 where(ica < 1) ica = 1 605 icell_tot = product(ica) 606 607 write(message,'(3a,2i8,2a,3i8)')& 608 & 'UPDLIS0 : SYSTEM SUBCELL DISTRIBUTION: ',ch10,& 609 & 'TOT. NO. OF CELLS (& max.) : ',icell_tot, icellx_a,ch10,& 610 & 'ICX,ICY,ICZ = ',ica(1),ica(2),ica(3) 611 call wrtout(std_out,message,'COLL') 612 613 614 if(icell_tot > icellx_a) then 615 write(message,'(a,i8,2a)')& 616 & 'ERROR IN UPD_LIS0: ICELL_TOT = ',icell_tot,ch10,& 617 & 'ERROR IN UPD_LIS0: RAISE ICELLX_A' 618 ABI_ERROR(message) 619 endif 620 621 !-- clears head vector & constructs head & list 622 head(:icellx_a) = 0 623 624 do i=1,lotfvar%natom 625 rdum = pbc_bb_proj(tau0(:,i)) 626 icell = 1 + mod(int( (rdum(1)+0.5)* float(ica(1)) ),ica(1)) & 627 + mod(int( (rdum(2)+0.5)* float(ica(2)) ),ica(2)) * ica(1) & 628 + mod(int( (rdum(3)+0.5)* float(ica(3)) ),ica(3)) * ica(1) * ica(2) 629 630 list(i) = head(icell) 631 head(icell) = i 632 enddo 633 634 !--Constructs the 13 neighbours of each cell 635 do IZ = 1, ica(3) 636 do IY = 1, ica(2) 637 do IX = 1, ica(1) 638 ic = icf(ix ,iy,iz ,ica(1),ica(2),ica(3)) 639 neigh_cel(ic,1) = icf(ix+1,iy,iz ,ica(1),ica(2),ica(3)) 640 neigh_cel(ic,2) = icf(ix+1,iy+1,iz ,ica(1),ica(2),ica(3)) 641 neigh_cel(ic,3) = icf(ix ,iy+1,iz ,ica(1),ica(2),ica(3)) 642 neigh_cel(ic,4) = icf(ix-1,iy+1,iz ,ica(1),ica(2),ica(3)) 643 neigh_cel(ic,5) = icf(ix+1,iy ,iz-1,ica(1),ica(2),ica(3)) 644 neigh_cel(ic,6) = icf(ix+1,iy+1,iz-1,ica(1),ica(2),ica(3)) 645 neigh_cel(ic,7) = icf(ix ,iy+1,iz-1,ica(1),ica(2),ica(3)) 646 neigh_cel(ic,8) = icf(ix-1,iy+1,iz-1,ica(1),ica(2),ica(3)) 647 neigh_cel(ic,9) = icf(ix+1,iy ,iz+1,ica(1),ica(2),ica(3)) 648 neigh_cel(ic,10) = icf(ix+1,iy+1,iz+1,ica(1),ica(2),ica(3)) 649 neigh_cel(ic,11) = icf(ix ,iy+1,iz+1,ica(1),ica(2),ica(3)) 650 neigh_cel(ic,12) = icf(ix-1,iy+1,iz+1,ica(1),ica(2),ica(3)) 651 neigh_cel(ic,13) = icf(ix ,iy ,iz+1,ica(1),ica(2),ica(3)) 652 enddo 653 enddo 654 enddo 655 656 !--Safety loops to avoid repetitions 657 658 !--(1) to avoid having twice the same neigh. cell 659 do icell = 1,icell_tot 660 do icnr = 1,13 661 icn = neigh_cel(icell,icnr) 662 do icnr2 = icnr+1,13 663 icn2 = neigh_cel(icell,icnr2) 664 if(icn2==icn) neigh_cel(icell,icnr2) = icell 665 enddo 666 enddo 667 enddo 668 669 !--(2) to avoid counting twice the interaction between the same 670 ! couple of neigh. cells 671 do icell = 1,icell_tot 672 do icnr = 1,13 673 icn = neigh_cel(icell,icnr) 674 do icnr2 = 1,13 675 icn2 = neigh_cel(icn,icnr2) 676 if(icn2==icell) neigh_cel(icn,icnr2) = icn 677 enddo 678 enddo 679 enddo 680 681 !--(3) constructs neighbour list looking through neighbour cells only 682 do icell = 1,icell_tot 683 iat = head(icell) 684 do while(iat > 0) 685 i = list(iat) 686 do while(i > 0) 687 if(i==iat) ABI_ERROR("i==iat") 688 call dist_pbc(tau0(:,i),tau0(:,iat)) 689 if(r2 < rcut2) then 690 nneig(iat) = nneig(iat) + 1 691 nneig(i) = nneig(i) + 1 692 neighl(nneig(iat),iat) = i 693 neighl(nneig(i) ,i ) = iat 694 endif ! distance check 695 i = list(i) 696 enddo 697 698 if(ANY(nneig(:) > lotfvar%nneigx)) then 699 write(message,'(a,i8,a)')& 700 'UPD_LIS CLASSIC: max no. of neighbours: ',lotfvar%nneigx,' is too small' 701 ABI_ERROR(message) 702 endif 703 704 705 ! mask to avoid self interaction within the same cell 706 ! considered as a neighbor and with different cells more than once 707 do icnr = 1,13 708 icn = neigh_cel(icell,icnr) 709 if(icn==icell) cycle 710 i = head(icn) 711 712 do while(i>0) 713 if(i==iat) ABI_ERROR("i==iat") 714 call dist_pbc(tau0(:,i),tau0(:,iat)) 715 if(r2 < rcut2) then 716 nneig(iat) = nneig(iat) + 1 717 nneig(i) = nneig(i) + 1 718 neighl(nneig(iat),iat) = i 719 neighl(nneig(i) ,i ) = iat 720 endif ! distance check 721 i = list(i) 722 enddo 723 if(ANY(nneig(:) > lotfvar%nneigx)) then 724 write(message,'(a,i8,a)')& 725 'UPD_LIS CLASSIC: max no. of neighbours: ',lotfvar%nneigx,' is too small' 726 ABI_ERROR(message) 727 endif 728 enddo 729 730 iat = list(iat) 731 enddo 732 enddo 733 end subroutine upd_lis0