TABLE OF CONTENTS
ABINIT/m_pred_langevin [ Modules ]
NAME
m_pred_langevin
FUNCTION
COPYRIGHT
Copyright (C) 1998-2024 ABINIT group (DCA, XG, GMR, JCC, SE) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
16 #if defined HAVE_CONFIG_H 17 #include "config.h" 18 #endif 19 20 #include "abi_common.h" 21 22 module m_pred_langevin 23 24 use defs_basis 25 use m_abicore 26 use m_abimover 27 use m_abihist 28 29 use m_numeric_tools, only : uniformrandom 30 use m_geometry, only : xcart2xred, xred2xcart, metric 31 use m_results_gs , only : results_gs_type 32 33 implicit none 34 35 private
ABINIT/pred_langevin [ Functions ]
NAME
pred_langevin
FUNCTION
Ionmov predictors (9) Langevin dynamics algorithm IONMOV 9: Uses a Langevin dynamics algorithm : see J. Chelikowsky, J. Phys. D : Appl Phys. 33(2000)R33 [[cite:Chelikowsky2000]]
INPUTS
ab_mover <type(abimover)> : Datatype with all the information needed by the preditor itime : Index of the present iteration ntime : Maximal number of iterations icycle : Index of the present cycle ncycle : Maximal number of cycles zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
71 subroutine pred_langevin(ab_mover,hist,icycle,itime,ncycle,ntime,zDEBUG,iexit,skipcycle) 72 73 implicit none 74 75 !Arguments ------------------------------------ 76 !scalars 77 type(abimover),intent(in) :: ab_mover 78 type(abihist),intent(inout) :: hist 79 integer,intent(in) :: itime 80 integer,intent(in) :: ntime 81 integer,intent(in) :: iexit 82 integer,intent(in) :: icycle 83 integer,intent(inout) :: ncycle 84 logical,intent(in) :: zDEBUG 85 logical,intent(out) :: skipcycle 86 87 !Local variables------------------------------- 88 !scalars 89 integer :: ii,kk,iatom,idim,iatom1,iatom2,itypat,idum=5,ihist_prev,mcfac 90 real(dp) :: ucvol,ucvol_next 91 real(dp),parameter :: v2tol=tol8 92 real(dp) :: etotal,rescale_vel,ran_num1,ran_num2 93 real(dp) :: ktemp,dist,distx,disty,distz,maxp1,maxp2,v2nose 94 real(dp) :: sig_gauss,delxi 95 logical :: jump_end_of_cycle=.FALSE. 96 character(len=5000) :: message 97 !arrays 98 integer,allocatable :: imax_perm(:) 99 real(dp),allocatable,save :: max_perm(:),pot_perm(:) 100 real(dp),allocatable,save :: ran_force(:,:),lang_force(:,:) 101 real(dp),allocatable,save :: fcart_mold(:,:),fcart_m(:,:) 102 103 real(dp) :: acell(3),acell_next(3) 104 real(dp) :: rprim(3,3),rprimd(3,3),rprimd_next(3,3),rprim_next(3,3) 105 real(dp) :: gprimd(3,3),gmet(3,3),rmet(3,3) 106 real(dp) :: fcart(3,ab_mover%natom) 107 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 108 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 109 real(dp) :: vel(3,ab_mover%natom) 110 real(dp) :: strten(6) 111 112 !*************************************************************************** 113 !Beginning of executable session 114 !*************************************************************************** 115 116 if(iexit/=0)then 117 if (allocated(pot_perm)) then 118 ABI_FREE(pot_perm) 119 end if 120 if (allocated(max_perm)) then 121 ABI_FREE(max_perm) 122 end if 123 if (allocated(imax_perm)) then 124 ABI_FREE(imax_perm) 125 end if 126 if (allocated(ran_force)) then 127 ABI_FREE(ran_force) 128 end if 129 if (allocated(lang_force)) then 130 ABI_FREE(lang_force) 131 end if 132 if (allocated(fcart_mold)) then 133 ABI_FREE(fcart_mold) 134 end if 135 if (allocated(fcart_m)) then 136 ABI_FREE(fcart_m) 137 end if 138 return 139 end if 140 141 jump_end_of_cycle=.FALSE. 142 143 !write(std_out,*) 'langevin 02',jump_end_of_cycle 144 !########################################################## 145 !### 02. Allocate the arrays 146 !### These arrays could be allocated from a previus 147 !### dataset that exit before itime==ntime 148 149 if(itime==1)then 150 if (allocated(pot_perm)) then 151 ABI_FREE(pot_perm) 152 end if 153 if (allocated(max_perm)) then 154 ABI_FREE(max_perm) 155 end if 156 if (allocated(imax_perm)) then 157 ABI_FREE(imax_perm) 158 end if 159 if (allocated(ran_force)) then 160 ABI_FREE(ran_force) 161 end if 162 if (allocated(lang_force)) then 163 ABI_FREE(lang_force) 164 end if 165 if (allocated(fcart_mold)) then 166 ABI_FREE(fcart_mold) 167 end if 168 if (allocated(fcart_m)) then 169 ABI_FREE(fcart_m) 170 end if 171 end if 172 173 if (.not.allocated(pot_perm)) then 174 ABI_MALLOC(pot_perm,(ab_mover%natom)) 175 end if 176 if (.not.allocated(max_perm)) then 177 ABI_MALLOC(max_perm,(ab_mover%ntypat)) 178 end if 179 if (.not.allocated(imax_perm)) then 180 ABI_MALLOC(imax_perm,(ab_mover%ntypat)) 181 end if 182 if (.not.allocated(ran_force)) then 183 ABI_MALLOC(ran_force,(3,ab_mover%natom)) 184 end if 185 if (.not.allocated(lang_force)) then 186 ABI_MALLOC(lang_force,(3,ab_mover%natom)) 187 end if 188 if (.not.allocated(fcart_mold)) then 189 ABI_MALLOC(fcart_mold,(3,ab_mover%natom)) 190 end if 191 if (.not.allocated(fcart_m)) then 192 ABI_MALLOC(fcart_m,(3,ab_mover%natom)) 193 end if 194 195 !write(std_out,*) 'langevin 03',jump_end_of_cycle 196 !########################################################## 197 !### 03. Obtain the present values from the history 198 199 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 200 201 fcart(:,:)=hist%fcart(:,:,hist%ihist) 202 strten(:) =hist%strten(:,hist%ihist) 203 vel(:,:) =hist%vel(:,:,hist%ihist) 204 etotal =hist%etot(hist%ihist) 205 206 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 207 do ii=1,3 208 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 209 end do 210 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 211 212 !write(std_out,*) 'langevin 04',jump_end_of_cycle 213 !########################################################## 214 !### 04. Compute the next values (Only for the first cycle) 215 216 if (icycle==1) then 217 218 ! The temperature is linear between initial and final values 219 ! It is here converted from Kelvin to Hartree (kb_HaK) 220 ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK 221 ! write(std_out,*) 'KTEMP=',ktemp 222 ! write(std_out,*) 'MDITEMP=',ab_mover%mdtemp(1) 223 ! write(std_out,*) 'MDFTEMP=',ab_mover%mdtemp(2) 224 ! write(std_out,*) 'DELAYPERM=',ab_mover%delayperm 225 226 227 ! %%% LANGEVIN DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 228 229 ! 1. Next values are the present ones 230 acell_next(:)=acell(:) 231 ucvol_next=ucvol 232 rprim_next(:,:)=rprim(:,:) 233 rprimd_next(:,:)=rprimd(:,:) 234 235 if (zDEBUG) then 236 write (std_out,*) '1. Next values are the present ones' 237 write(std_out,*) 'RPRIMD' 238 do kk=1,3 239 write(std_out,*) rprimd(:,kk) 240 end do 241 write(std_out,*) 'RPRIM' 242 do kk=1,3 243 write(std_out,*) rprim(:,kk) 244 end do 245 write(std_out,*) 'ACELL' 246 write(std_out,*) acell(:) 247 end if 248 249 if(itime==1)then 250 251 ! 2. Compute twice the kinetic energy of the system, called v2nose 252 v2nose=0.0_dp 253 do iatom=1,ab_mover%natom 254 do idim=1,3 255 v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 256 end do 257 end do 258 if (zDEBUG) then 259 write (std_out,*) '2. Compute twice the kinetic energy of the system, called v2nose' 260 write (std_out,*) 'V2NOSE=',v2nose 261 end if 262 263 ! 3. If there is no kinetic energy, use random numbers 264 if (v2nose<=v2tol) then 265 v2nose=0.0_dp 266 do iatom=1,ab_mover%natom 267 do idim=1,3 268 ! uniformrandom returns a uniform random deviate between 0.0 and 1.0 269 ! if it were always 0 or 1, then the following expression 270 ! would give the requested temperature 271 vel(idim,iatom)=(1.0_dp-2.0_dp*uniformrandom(idum))*sqrt(ktemp/ab_mover%amass(iatom)) 272 ! Recompute v2nose 273 v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 274 end do 275 end do 276 end if 277 278 if (zDEBUG) then 279 write (std_out,*) '3. If there is no kinetic energy, use random numbers' 280 write (std_out,*) 'VEL' 281 do kk=1,ab_mover%natom 282 write (std_out,*) vel(:,kk) 283 end do 284 write (std_out,*) 'V2NOSE=',v2nose 285 end if 286 287 288 ! Now, rescale the velocities to give the proper temperature 289 rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose) 290 vel(:,:)=vel(:,:)*rescale_vel 291 ! Recompute v2nose with the rescaled velocities 292 v2nose=0.0_dp 293 do iatom=1,ab_mover%natom 294 do idim=1,3 295 v2nose=v2nose+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 296 end do 297 end do 298 write(message, '(a)' )& 299 & ' Rescaling or initializing velocities to initial temperature' 300 call wrtout(ab_out,message,'COLL') 301 call wrtout(std_out,message,'COLL') 302 write(message, '(a,D12.5,a,D12.5)' )& 303 & ' --- Scaling factor : ',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1) 304 call wrtout(ab_out,message,'COLL') 305 call wrtout(std_out,message,'COLL') 306 write(message, '(a,D12.5)' )& 307 & ' --- Effective temperature',v2nose/3.0_dp/(kb_HaK*ab_mover%natom) 308 call wrtout(ab_out,message,'COLL') 309 call wrtout(std_out,message,'COLL') 310 ! end if itime==0 311 end if 312 313 ! This section is devoted to the optional atom permutation (JYR 001114) 314 ! Two input variables are needed 315 ! ab_mover%delayperm : is the interval (in time steps) at which 316 ! atoms are tentatively permuted 317 ! default value could be 0 318 ! ab_mover%signperm : is the type of bias for the permutation 319 ! +1 to favor alternation of species 320 ! -1 to favor segregation 321 322 ! Force no permutation at initial step 323 if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2) then 324 if (mod(itime-1,ab_mover%delayperm)==0) then 325 ! Try commutation of atoms. 326 write(message, '(a)')' Attempt of commutation ' 327 call wrtout(ab_out,message,'COLL') 328 call wrtout(std_out,message,'COLL') 329 ! Compute a 'permutation potential' 330 do iatom=1,ab_mover%natom 331 pot_perm(iatom)=0.0_dp 332 do iatom1=1,ab_mover%natom 333 if (iatom1.ne.iatom) then 334 distx=xcart(1,iatom)-xcart(1,iatom1) 335 distx=distx-acell(1)*nint(distx/acell(1)) 336 disty=xcart(2,iatom)-xcart(2,iatom1) 337 disty=disty-acell(2)*nint(disty/acell(2)) 338 distz=xcart(3,iatom)-xcart(3,iatom1) 339 distz=distz-acell(3)*nint(distz/acell(3)) 340 ! Here we count each atom below 2 angstr as 1, could be customized 341 dist=sqrt(distx*distx+disty*disty+distz*distz)/3.7807 342 write(std_out,*) iatom,iatom1,dist 343 if (ab_mover%typat(iatom).ne.ab_mover%typat(iatom1)) then 344 mcfac=-1 345 else 346 mcfac=1 347 end if 348 if (dist<1.0_dp) dist=1.0_dp 349 pot_perm(iatom)=pot_perm(iatom)+mcfac*(ab_mover%signperm)*1.0_dp& 350 & /exp(log(dist)*6.0_dp) 351 end if 352 end do 353 end do 354 write(std_out,*) ' Perm_pot ',pot_perm(:) 355 ! write(message, '(a,10f12.5)' )' Perm_pot ',& 356 ! & (pot_perm(iatom1),iatom1=1,ab_mover%natom) 357 ! call wrtout(ab_out,message,'COLL') 358 ! call wrtout(std_out,message,'COLL') 359 360 ! Find the two atoms, of different types, with the highest perm_pot 361 max_perm(:)=-1.0d9 362 do iatom=1,ab_mover%natom 363 if (pot_perm(iatom) > max_perm(ab_mover%typat(iatom))) then 364 max_perm(ab_mover%typat(iatom))=pot_perm(iatom) 365 imax_perm(ab_mover%typat(iatom))=iatom 366 end if 367 end do 368 369 if(zDEBUG)then 370 ! write(message, '(a,10f12.5)' )' max_Perm ',& 371 ! & (max_perm(itypat),itypat=1,ab_mover%ntypat) 372 ! call wrtout(std_out,message,'COLL') 373 ! write(message, '(a,10i12)' )' imax_Perm ',& 374 ! & (imax_perm(itypat),itypat=1,ab_mover%ntypat) 375 ! call wrtout(std_out,message,'COLL') 376 write(std_out,*) 'NTYPAT',ab_mover%ntypat 377 write(message, '(a,10f12.5)' )' max_Perm ',& 378 & (max_perm(:)) 379 call wrtout(std_out,message,'COLL') 380 write(message, '(a,10i12)' )' imax_Perm ',& 381 & (imax_perm(:)) 382 call wrtout(std_out,message,'COLL') 383 end if 384 385 ! Loop and keep the 2 largest values 386 if (max_perm(1)>max_perm(2)) then 387 maxp1=max_perm(1) 388 maxp2=max_perm(2) 389 iatom1=imax_perm(1) 390 iatom2=imax_perm(2) 391 else 392 maxp1=max_perm(2) 393 maxp2=max_perm(1) 394 iatom1=imax_perm(2) 395 iatom2=imax_perm(1) 396 end if 397 398 do itypat=3,ab_mover%ntypat 399 if (max_perm(itypat)>maxp1) then 400 maxp2=maxp1 401 iatom2=iatom1 402 maxp1=max_perm(itypat) 403 iatom1=imax_perm(itypat) 404 else if (max_perm(itypat)>maxp2) then 405 maxp2=max_perm(itypat) 406 iatom2=imax_perm(itypat) 407 end if 408 end do 409 write(message, '(2(a,i5))' )' Will commute atom...',iatom1,'...of type ',& 410 & ab_mover%typat(iatom1) 411 call wrtout(ab_out,message,'COLL') 412 call wrtout(std_out,message,'COLL') 413 write(message, '(2(a,i5))' )' with atom...',iatom2,'...of type ',& 414 & ab_mover%typat(iatom2) 415 call wrtout(ab_out,message,'COLL') 416 call wrtout(std_out,message,'COLL') 417 418 ! Commute the atoms positions 419 distx=xcart(1,iatom1) 420 disty=xcart(2,iatom1) 421 distz=xcart(3,iatom1) 422 xcart(1,iatom1)=xcart(1,iatom2) 423 xcart(2,iatom1)=xcart(2,iatom2) 424 xcart(3,iatom1)=xcart(3,iatom2) 425 xcart(1,iatom2)=distx 426 xcart(2,iatom2)=disty 427 xcart(3,iatom2)=distz 428 ! Convert back to xred (reduced coordinates) 429 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 430 431 end if ! if(mod(itime,ab_mover%delayperm)==0) 432 433 else 434 ! write(std_out,*) "I will jump to the end of cycle" 435 jump_end_of_cycle=.TRUE. 436 437 end if ! if (itime/=1 .and. ab_mover%delayperm/=0 .and. ab_mover%ntypat>2) 438 ! End of the commutation section 439 440 end if ! if (icycle==1) 441 442 if (allocated(imax_perm)) then 443 ABI_FREE(imax_perm) 444 end if 445 446 !write(std_out,*) 'langevin 05',jump_end_of_cycle 447 !########################################################## 448 !### 05. Compute the next values (Only for extra cycles) 449 450 if (icycle>1) then 451 452 ! write(std_out,*) "Entering internal cycle 2",icycle 453 454 ! If the energy computed is higher than the current 455 ! (etotal_temp) we have to discard the changes 456 ! and compute again 457 458 ! write(std_out,*) ch10 459 ! write(std_out,*) 'EVALUATION FORCES',etotal,hist%etot(abihist_findIndex(hist,-1)) 460 ! write(std_out,*) ch10 461 462 ! This is the worst case (2 evaluations of SCFCV) 463 ihist_prev = abihist_findIndex(hist,-1) 464 if (etotal>hist%etot(ihist_prev).and.icycle==2) then 465 466 ! Discard the changes 467 acell(:) =hist%acell(:,ihist_prev) 468 rprimd(:,:)=hist%rprimd(:,:,ihist_prev) 469 xred(:,:) =hist%xred(:,:,ihist_prev) 470 fcart(:,:) =hist%fcart(:,:,ihist_prev) 471 strten(:) =hist%strten(:,ihist_prev) 472 vel(:,:) =hist%vel(:,:,ihist_prev) 473 etotal =hist%etot(ihist_prev) 474 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 475 476 ! distx=xcart(1,iatom1) 477 ! disty=xcart(2,iatom1) 478 ! distz=xcart(3,iatom1) 479 ! xcart(1,iatom1)=xcart(1,iatom2) 480 ! xcart(2,iatom1)=xcart(2,iatom2) 481 ! xcart(3,iatom1)=xcart(3,iatom2) 482 ! xcart(1,iatom2)=distx 483 ! xcart(2,iatom2)=disty 484 ! xcart(3,iatom2)=distz 485 486 ! Convert back to xred (reduced coordinates) 487 ! call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 488 write(message, '(a)' )' Commutation unsuccessful, recomputing the forces' 489 call wrtout(ab_out,message,'COLL') 490 call wrtout(std_out,message,'COLL') 491 492 ! This is the best case (only 1 evaluation of SCFCV) 493 else 494 495 write(message, '(a)')' Commutation successful ! Going on' 496 call wrtout(ab_out,message,'COLL') 497 call wrtout(std_out,message,'COLL') 498 499 ! In thisc case we do not need to compute again SCFCV 500 ! We avoid the second iteration on ii 501 jump_end_of_cycle=.TRUE. 502 503 end if ! etotal > etotal_temp 504 505 end if ! if (icycle>1) 506 507 !write(std_out,*) 'langevin 06',jump_end_of_cycle 508 !########################################################## 509 !### 06. Compute the next values (Only for the last cycle) 510 511 if(jump_end_of_cycle)then 512 ! icycle=ncycle 513 ! write(std_out,*) 'This is the last cycle, avoid the others and continue' 514 skipcycle=.TRUE. 515 else 516 skipcycle=.FALSE. 517 end if 518 519 if ((icycle==ncycle).OR.(skipcycle)) then 520 521 ! write(std_out,*) 'ENTERING THE FINAL PART',icycle,ncycle 522 523 ! Specific to Langevin dynamics 524 ! Initialize an array of random forces 525 ! No random force at itime=0 526 ! if (itime==0) then 527 if (itime<0) then 528 529 ran_force(:,:)=0.0_dp 530 531 else 532 533 do iatom=1,ab_mover%natom 534 ! sig_gauss is the std deviation of the random distribution 535 sig_gauss=sqrt(2.0_dp*(ab_mover%friction)*ab_mover%amass(iatom)*ktemp) 536 ! write(std_out,*) 'sig_gauss=',sig_gauss 537 ! write(std_out,*) 'friction=',ab_mover%friction 538 ! write(std_out,*) 'ktemp=',ktemp 539 do idim=1,3 540 delxi=2.0_dp 541 do while (delxi >= 1.0_dp) 542 ran_num1=2.0_dp*uniformrandom(idum)-1.0_dp 543 ran_num2=2.0_dp*uniformrandom(idum)-1.0_dp 544 delxi=ran_num1*ran_num1+ran_num2*ran_num2 545 ! write(std_out,*) delxi,ran_num1,ran_num2 546 end do 547 ran_force(idim,iatom)=ran_num1*sqrt(-2.0_dp*log(delxi)/delxi)& 548 & *sig_gauss/sqrt(ab_mover%dtion) 549 ! write(std_out,*) 'ran_force',ran_force(idim,iatom) 550 end do 551 end do 552 553 if (zDEBUG) then 554 write (std_out,*) '4. Different forces computed' 555 write (std_out,*) 'RAN_FORCE' 556 do kk=1,ab_mover%natom 557 write (std_out,*) ran_force(:,kk) 558 end do 559 end if 560 561 562 if(zDEBUG)then 563 ! The distribution should be gaussian 564 delxi=0.0_dp 565 do iatom=1,ab_mover%natom 566 do idim=1,3 567 delxi=delxi+(ran_force(idim,iatom)*ab_mover%dtion)**2 568 end do 569 end do 570 delxi=delxi/(3.0_dp*ab_mover%natom) 571 write(message, '(2(a,es22.14))' )' variance =',delxi,' asked =',& 572 & 2.0_dp*(ab_mover%friction)*ab_mover%amass(2)*ktemp*ab_mover%dtion 573 call wrtout(std_out,message,'COLL') 574 end if 575 ! end if itime\=0 576 577 end if 578 579 ! zDEBUG 580 ! write(message, '(a)' )' after initializing ran_force' 581 ! call wrtout(ab_out,message,'COLL') 582 ! call wrtout(std_out,message,'COLL') 583 ! ENDzDEBUG 584 585 do iatom=1,ab_mover%natom 586 do idim=1,3 587 fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom) 588 ran_force(idim,iatom)=ran_force(idim,iatom)/ab_mover%amass(iatom) 589 end do 590 end do 591 lang_force(:,:)=ran_force(:,:)-(ab_mover%friction)*vel(:,:)+fcart_m(:,:) 592 593 if (zDEBUG) then 594 write (std_out,*) '4. Different forces computed' 595 write (std_out,*) 'FCART_M' 596 do kk=1,ab_mover%natom 597 write (std_out,*) fcart_m(:,kk) 598 end do 599 write (std_out,*) 'RAN_FORCE' 600 do kk=1,ab_mover%natom 601 write (std_out,*) ran_force(:,kk) 602 end do 603 write (std_out,*) 'LANG_FORCE' 604 do kk=1,ab_mover%natom 605 write (std_out,*) lang_force(:,kk) 606 end do 607 end if 608 609 ! zDEBUG 610 ! write(message, '(a)' )'before verlet' 611 ! call wrtout(ab_out,message,'COLL') 612 ! call wrtout(std_out,message,'COLL') 613 ! ENDzDEBUG 614 615 ! Compute next atomic coordinates using Verlet algorithm 616 617 ! Impose no change of acell, ucvol, rprim, and rprimd 618 acell_next(:)=acell(:) 619 ucvol_next=ucvol 620 rprim_next(:,:)=rprim(:,:) 621 rprimd_next(:,:)=rprimd(:,:) 622 623 ! Convert input xred (reduced coordinates) to xcart (cartesian) 624 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 625 ! Uses the velocity 626 ! 627 ! If an atom wants to cross the walls, velocity is reversed. 628 ! 629 do iatom=1,ab_mover%natom 630 do idim=1,3 631 delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ & 632 & 0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom) 633 if ( (delxi > (rprimd(idim,idim)+(ab_mover%mdwall)) ) .or. & 634 & (delxi < - (ab_mover%mdwall) ) ) then 635 vel(idim,iatom)=-vel(idim,iatom) 636 delxi=xcart(idim,iatom)+ab_mover%dtion*vel(idim,iatom)+ & 637 & 0.5_dp*ab_mover%dtion*ab_mover%dtion*lang_force(idim,iatom) 638 end if 639 xcart_next(idim,iatom)=delxi 640 end do 641 end do 642 xcart(:,:)=xcart_next(:,:) 643 if (zDEBUG) then 644 write (std_out,*) '5. If an atom wants to cross the walls, velocity is reversed.' 645 write (std_out,*) 'XCART' 646 do kk=1,ab_mover%natom 647 write (std_out,*) xcart(:,kk) 648 end do 649 end if 650 651 ! Convert back to xred_next (reduced coordinates) 652 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 653 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 654 655 if (itime==1) then 656 ! no old forces are available at first step 657 ! Simple update of the velocity 658 ! first compute vel_nexthalf for next steps 659 vel(:,:)=vel(:,:)+ab_mover%dtion*lang_force(:,:) 660 else 661 ! case itime /= 0 normal verlet integration 662 vel(:,:)=vel(:,:)+0.5_dp*ab_mover%dtion*(fcart_mold(:,:)+lang_force(:,:)) 663 end if 664 if (zDEBUG) then 665 write (std_out,*) '5. Change velocity with verlet' 666 write (std_out,*) 'VEL' 667 do kk=1,ab_mover%natom 668 write (std_out,*) vel(:,kk) 669 end do 670 end if 671 672 ! Store 'current force' as 'old force' 673 fcart_mold(:,:)=lang_force(:,:) 674 675 end if ! if (icycle==ncycle) 676 677 !write(std_out,*) 'langevin 07',jump_end_of_cycle 678 !########################################################## 679 !### 07. Update the history with the prediction 680 681 !Increase indexes 682 hist%ihist = abihist_findIndex(hist,+1) 683 684 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 685 hist%vel(:,:,hist%ihist)=vel(:,:) 686 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 687 688 if (ab_mover%delayperm==0 .or. ab_mover%ntypat<=2) ncycle=1 689 if(itime==ntime-1) ncycle=1 690 691 end subroutine pred_langevin