TABLE OF CONTENTS
- ABINIT/lotfpath
- lotfpath/init_lotf
- lothpath/alpha_update
- lothpath/end_lotf
- lothpath/extrapolation_loop
- lothpath/fitclus
- lothpath/force0
- lothpath/force_to_vel
- lothpath/intparms
- lothpath/lotf_extrapolation
- lothpath/lotf_interpolation
- lothpath/upd_lis
- lothpath/vel_rescale
- lothpath/vel_to_gauss
ABINIT/lotfpath [ 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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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 ]
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