TABLE OF CONTENTS
ABINIT/m_pred_isokinetic [ Modules ]
NAME
m_pred_isokinetic
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_isokinetic 23 24 use m_abicore 25 use defs_basis 26 use m_abimover 27 use m_abihist 28 29 use m_numeric_tools, only : uniformrandom 30 use m_geometry, only : xcart2xred, xred2xcart 31 32 33 implicit none 34 35 private
ABINIT/pred_isokinetic [ Functions ]
NAME
pred_isokinetic
FUNCTION
Ionmov predictors (12) Isokinetic ensemble molecular dynamics IONMOV 12: Isokinetic ensemble molecular dynamics. The equation of motion of the ions in contact with a thermostat are solved with the algorithm proposed by Zhang [J. Chem. Phys. 106, 6102 (1997)] [[cite:Zhang1997]], as worked out by Minary et al, J. Chem. Phys. 188, 2510 (2003) [[cite:Minary2003]]. The conservation of the kinetic energy is obtained within machine precision, at each step. Related parameters: the time step (dtion), the initial temperature (mdtemp(1)) if the velocities are not defined to start with.
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 zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
72 subroutine pred_isokinetic(ab_mover,hist,itime,ntime,zDEBUG,iexit) 73 74 implicit none 75 76 !Arguments ------------------------------------ 77 !scalars 78 integer,intent(in) :: itime 79 integer,intent(in) :: ntime 80 integer,intent(in) :: iexit 81 logical,intent(in) :: zDEBUG 82 type(abimover),intent(in) :: ab_mover 83 type(abihist),intent(inout) :: hist 84 85 !Local variables------------------------------- 86 !scalars 87 integer :: kk,iatom,idim,idum=5,nxyzatfree,ndegfreedom,nfirst,ifirst 88 real(dp) :: a,as,b,sqb,s,s1,s2,scdot,sigma2,vtest,v2gauss 89 real(dp),parameter :: v2tol=tol8 90 real(dp) :: etotal,rescale_vel 91 character(len=5000) :: message 92 !arrays 93 real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:) 94 95 real(dp) :: acell(3),rprimd(3,3) 96 real(dp) :: fcart(3,ab_mover%natom) 97 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 98 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 99 real(dp) :: vel(3,ab_mover%natom) 100 real(dp) :: strten(6) 101 102 !*************************************************************************** 103 !Beginning of executable session 104 !*************************************************************************** 105 106 !DEBUG 107 !write(std_out,*)' pred_isokinetic : enter ' 108 !stop 109 !ENDDEBUG 110 111 if(iexit/=0)then 112 ABI_SFREE(fcart_m) 113 ABI_SFREE(vel_nexthalf) 114 return 115 end if 116 117 !write(std_out,*) 'isokinetic 01' 118 !########################################################## 119 !### 01. Debugging and Verbose 120 121 if(zDEBUG)then 122 write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),& 123 & 'Debugging and Verbose for pred_isokinetic',('-',kk=1,37) 124 write(std_out,*) 'ionmov: ',12 125 write(std_out,*) 'itime: ',itime 126 end if 127 128 !write(std_out,*) 'isokinetic 02' 129 !########################################################## 130 !### 02. Allocate the vectors vin, vout and hessian matrix 131 !### These arrays could be allocated from a previous 132 !### dataset that exit before itime==ntime 133 134 if(itime==1)then 135 ABI_SFREE(fcart_m) 136 ABI_SFREE(vel_nexthalf) 137 end if 138 139 if (.not.allocated(fcart_m)) then 140 ABI_MALLOC(fcart_m,(3,ab_mover%natom)) 141 end if 142 if (.not.allocated(vel_nexthalf)) then 143 ABI_MALLOC(vel_nexthalf,(3,ab_mover%natom)) 144 end if 145 146 !write(std_out,*) 'isokinetic 03' 147 !########################################################## 148 !### 03. Obtain the present values from the history 149 150 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 151 152 fcart(:,:)=hist%fcart(:,:,hist%ihist) 153 strten(:) =hist%strten(:,hist%ihist) 154 vel(:,:) =hist%vel(:,:,hist%ihist) 155 etotal =hist%etot(hist%ihist) 156 157 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 158 159 if(zDEBUG)then 160 write (std_out,*) 'fcart:' 161 do kk=1,ab_mover%natom 162 write (std_out,*) fcart(:,kk) 163 end do 164 write (std_out,*) 'vel:' 165 do kk=1,ab_mover%natom 166 write (std_out,*) vel(:,kk) 167 end do 168 write (std_out,*) 'strten:' 169 write (std_out,*) strten(1:3),ch10,strten(4:6) 170 write (std_out,*) 'etotal:' 171 write (std_out,*) etotal 172 end if 173 174 !Count the number of degrees of freedom, taking into account iatfix. 175 !Also fix the velocity to zero for the fixed atoms 176 nxyzatfree=0 177 do iatom=1,ab_mover%natom 178 do idim=1,3 179 if(ab_mover%iatfix(idim,iatom)==0)then 180 nxyzatfree=nxyzatfree+1 181 else 182 vel(idim,iatom)=zero 183 endif 184 enddo 185 enddo 186 187 !Now, the number of degrees of freedom is reduced by four because of the kinetic energy conservation 188 !and because of the conservation of the total momentum for each dimension, in case no atom position is fixed for that dimension 189 !(in the latter case, one degree of freedom has already been taken away) 190 !This was not done until v8.9 of ABINIT ... 191 ndegfreedom=nxyzatfree 192 ndegfreedom=nxyzatfree-1 ! Kinetic energy conservation 193 do idim=1,3 194 if(sum(ab_mover%iatfix(idim,:))==0)then 195 ndegfreedom=ndegfreedom-1 ! Macroscopic momentum 196 endif 197 enddo 198 199 !write(std_out,*) 'isokinetic 04' 200 !########################################################## 201 !### 04. Second half-velocity (Only after the first itime) 202 203 if(itime>1) then 204 205 do iatom=1,ab_mover%natom 206 do idim=1,3 207 if(ab_mover%iatfix(idim,iatom)==0)then 208 fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom) 209 else 210 fcart_m(idim,iatom)=zero 211 endif 212 end do 213 end do 214 215 ! Computation of vel(:,:) at the next positions 216 ! Computation of v2gauss, actually twice the kinetic energy. 217 ! Called 2K, cf Eq. (A13) of [[cite:Minary2003]]. 218 v2gauss=0.0_dp 219 do iatom=1,ab_mover%natom 220 do idim=1,3 221 v2gauss=v2gauss+& 222 & vel_nexthalf(idim,iatom)*vel_nexthalf(idim,iatom)*& 223 & ab_mover%amass(iatom) 224 end do 225 end do 226 227 ! Computation of a and b (4.13 of [[cite:Minary2003]]) 228 a=0.0_dp 229 b=0.0_dp 230 do iatom=1,ab_mover%natom 231 do idim=1,3 232 a=a+fcart_m(idim,iatom)*vel_nexthalf(idim,iatom)*ab_mover%amass(iatom) 233 b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom) 234 end do 235 end do 236 a=a/v2gauss 237 b=b/v2gauss 238 239 240 ! Computation of s and scdot 241 sqb=sqrt(b) 242 as=sqb*ab_mover%dtion/2. 243 ! jmb 244 if ( as > 300.0 ) as=300.0 245 s1=cosh(as) 246 s2=sinh(as) 247 s=a*(s1-1.)/b+s2/sqb 248 scdot=a*s2/sqb+s1 249 250 do iatom=1,ab_mover%natom 251 do idim=1,3 252 if(ab_mover%iatfix(idim,iatom)==0)then 253 vel(idim,iatom)=(vel_nexthalf(idim,iatom)+fcart_m(idim,iatom)*s)/scdot 254 else 255 vel(idim,iatom)=zero 256 endif 257 enddo 258 enddo 259 260 if (zDEBUG)then 261 write(std_out,*) 'Computation of the second half-velocity' 262 write(std_out,*) 'Cartesian forces per atomic mass (fcart_m):' 263 do kk=1,ab_mover%natom 264 write (std_out,*) fcart_m(:,kk) 265 end do 266 write(std_out,*) 'vel:' 267 do kk=1,ab_mover%natom 268 write (std_out,*) vel(:,kk) 269 end do 270 write(std_out,*) 'v2gauss:',v2gauss 271 write(std_out,*) 'a:',a 272 write(std_out,*) 'b:',b 273 write(std_out,*) 's:',s 274 write(std_out,*) 'scdot:',scdot 275 end if 276 277 end if ! (if itime>1) 278 279 !write(std_out,*) 'isokinetic 05' 280 !########################################################## 281 !### 05. First half-time (First cycle the loop is double) 282 283 if (itime==1) then 284 nfirst=2 285 else 286 nfirst=1 287 end if 288 289 do ifirst=1,nfirst 290 291 ! Application of Gauss' principle of least constraint according to Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, [[cite:Zhang1997]] p.6102) 292 293 ! v2gauss is twice the kinetic energy 294 v2gauss=0.0_dp 295 do iatom=1,ab_mover%natom 296 do idim=1,3 297 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 298 end do 299 end do 300 301 ! If there is no kinetic energy to start with ... 302 if (v2gauss<=v2tol.and.itime==1) then 303 ! Maxwell-Boltzman distribution 304 v2gauss=zero 305 vtest=zero 306 do iatom=1,ab_mover%natom 307 do idim=1,3 308 if(ab_mover%iatfix(idim,iatom)==0)then 309 vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum)) 310 vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum))) 311 else 312 vel(idim,iatom)=zero 313 endif 314 end do 315 end do 316 317 ! Get rid of center-of-mass velocity 318 s1=sum(ab_mover%amass(:)) 319 do idim=1,3 320 if(sum(ab_mover%iatfix(idim,:))==0)then 321 s2=sum(ab_mover%amass(:)*vel(idim,:)) 322 vel(idim,:)=vel(idim,:)-s2/s1 323 endif 324 end do 325 326 ! Recompute v2gauss 327 do iatom=1,ab_mover%natom 328 do idim=1,3 329 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 330 vtest=vtest+vel(idim,iatom)/ndegfreedom 331 end do 332 end do 333 334 ! Now rescale the velocities to give the exact temperature 335 rescale_vel=sqrt(ndegfreedom*kb_HaK*ab_mover%mdtemp(1)/v2gauss) 336 vel(:,:)=vel(:,:)*rescale_vel 337 338 ! Recompute v2gauss with the rescaled velocities 339 v2gauss=zero 340 do iatom=1,ab_mover%natom 341 do idim=1,3 342 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 343 end do 344 end do 345 346 ! Compute the variance and print 347 sigma2=(v2gauss/ndegfreedom-ab_mover%amass(1)*vtest**2)/kb_HaK 348 349 end if 350 351 do iatom=1,ab_mover%natom 352 do idim=1,3 353 if(ab_mover%iatfix(idim,iatom)==0)then 354 fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom) 355 else 356 fcart_m(idim,iatom)=zero 357 endif 358 end do 359 end do 360 361 if (zDEBUG)then 362 write(std_out,*) 'Calculation first half-velocity ' 363 write (std_out,*) 'vel:' 364 do kk=1,ab_mover%natom 365 write (std_out,*) vel(:,kk) 366 end do 367 write (std_out,*) 'xcart:' 368 do kk=1,ab_mover%natom 369 write (std_out,*) xcart(:,kk) 370 end do 371 write (std_out,*) 'xred:' 372 do kk=1,ab_mover%natom 373 write (std_out,*) xred(:,kk) 374 end do 375 write (std_out,*) 'fcart_m' 376 do kk=1,ab_mover%natom 377 write (std_out,*) fcart_m(:,kk) 378 end do 379 write(std_out,*) 's2',s2 380 write(std_out,*) 'v2gauss',v2gauss 381 write(std_out,*) 'sigma2',sigma2 382 383 write(message, '(a)' )& 384 & ' --- Rescaling or initializing velocities to initial temperature' 385 call wrtout(std_out,message,'COLL') 386 write(message, '(a,d12.5,a,D12.5)' )& 387 & ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1) 388 call wrtout(std_out,message,'COLL') 389 write(message, '(a,d12.5,a,D12.5)' )& 390 & ' --- Effective temperature',v2gauss/(ndegfreedom*kb_HaK),' From variance', sigma2 391 call wrtout(std_out,message,'COLL') 392 end if 393 394 ! Convert input xred (reduced coordinates) to xcart (cartesian) 395 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 396 397 if(itime==1.and.ifirst==1) then 398 call wrtout(std_out,'if itime==1','COLL') 399 vel_nexthalf(:,:)=vel(:,:) 400 xcart_next(:,:)=xcart(:,:) 401 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 402 xred=xred_next 403 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 404 end if 405 406 end do 407 408 !Computation of vel_nexthalf (4.16 of [[cite:Minary2003]]) 409 !Computation of a and b (4.13 of [[cite:Minary2003]]) 410 a=0.0_dp 411 b=0.0_dp 412 do iatom=1,ab_mover%natom 413 do idim=1,3 414 a=a+fcart_m(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 415 b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom) 416 end do 417 end do 418 a=a/v2gauss+tol20 419 b=b/v2gauss+tol20 420 !Computation of s and scdot 421 sqb=sqrt(b)+tol20 422 as=sqb*ab_mover%dtion/2. 423 ! jmb 424 if ( as > 300.0 ) as=300.0 425 s1=cosh(as) 426 s2=sinh(as) 427 s=a*(s1-1.)/b+s2/sqb 428 scdot=a*s2/sqb+s1 429 do iatom=1,ab_mover%natom 430 do idim=1,3 431 if(ab_mover%iatfix(idim,iatom)==0)then 432 vel_nexthalf(idim,iatom)=(vel(idim,iatom)+fcart_m(idim,iatom)*s)/scdot 433 else 434 vel_nexthalf(idim,iatom)=zero 435 endif 436 enddo 437 enddo 438 439 !Computation of the next positions 440 xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion 441 442 if (zDEBUG)then 443 write(std_out,*) 'a:',a 444 write(std_out,*) 'b:',b 445 write(std_out,*) 's:',s 446 write(std_out,*) 'scdot:',scdot 447 end if 448 449 !Convert back to xred (reduced coordinates) 450 451 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 452 453 !write(std_out,*) 'isokinetic 06' 454 !########################################################## 455 !### 06. Update the history with the prediction 456 457 xcart=xcart_next 458 xred=xred_next 459 460 !increment the ihist 461 hist%ihist = abihist_findIndex(hist,+1) 462 463 !Fill the history with the variables 464 !xred, acell, rprimd, vel 465 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 466 hist%vel(:,:,hist%ihist)=vel(:,:) 467 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 468 469 if(zDEBUG)then 470 write (std_out,*) 'fcart:' 471 do kk=1,ab_mover%natom 472 write (std_out,*) fcart(:,kk) 473 end do 474 write (std_out,*) 'vel:' 475 do kk=1,ab_mover%natom 476 write (std_out,*) vel(:,kk) 477 end do 478 write (std_out,*) 'strten:' 479 write (std_out,*) strten(1:3),ch10,strten(4:6) 480 write (std_out,*) 'etotal:' 481 write (std_out,*) etotal 482 end if 483 484 if (.false.) write(std_out,*) ntime 485 486 end subroutine pred_isokinetic