TABLE OF CONTENTS
ABINIT/m_pred_nose [ Modules ]
NAME
m_pred_nose
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_nose 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 32 implicit none 33 34 private
ABINIT/pred_nose [ Functions ]
NAME
pred_nose
FUNCTION
Ionmov predictors (8) Verlet algorithm with a nose-hoover thermostat IONMOV 8: Given a starting point xred that is a vector of length 3*natom (reduced nuclei coordinates), a velocity vector (in cartesian coordinates), and unit cell parameters (acell and rprimd - without velocities in the present implementation), the Verlet dynamics is performed, using the gradient of the energy (atomic forces and stresses) as calculated by the routine scfcv. Some atoms can be kept fixed, while the propagation of unit cell parameters is only performed if optcell/=0. No more than "ntime" steps are performed. The time step is governed by dtion (contained in dtset) Returned quantities are xred, and eventually acell and rprimd (new ones!). See ionmov=6, but with a nose-hoover thermostat Velocity verlet algorithm : Swope et al JCP 76 (1982) 637
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
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
SOURCE
80 subroutine pred_nose(ab_mover,hist,itime,ntime,zDEBUG,iexit) 81 82 implicit none 83 84 !Arguments ------------------------------------ 85 !scalars 86 type(abimover),intent(in) :: ab_mover 87 type(abihist),intent(inout) :: hist 88 integer,intent(in) :: itime 89 integer,intent(in) :: ntime 90 integer,intent(in) :: iexit 91 logical,intent(in) :: zDEBUG 92 93 !Local variables------------------------------- 94 !scalars 95 integer :: ii,jj,kk 96 integer :: idum=-5 97 real(dp),parameter :: v2tol=tol8,nosetol=tol10 98 real(dp) :: delxi,xio,ktemp,rescale_vel 99 real(dp) :: dnose,v2nose,xin_nose 100 real(dp),save :: xi_nose,fsnose,snose 101 real(dp) :: gnose 102 real(dp) :: ucvol,ucvol_next 103 real(dp) :: etotal 104 logical :: ready 105 106 !arrays 107 real(dp) :: acell(3),acell_next(3) 108 real(dp) :: rprimd(3,3),rprimd_next(3,3) 109 real(dp) :: gprimd(3,3) 110 real(dp) :: gmet(3,3) 111 real(dp) :: rmet(3,3) 112 real(dp) :: fcart(3,ab_mover%natom) 113 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 114 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 115 real(dp) :: vel(3,ab_mover%natom),vel_temp(3,ab_mover%natom) 116 real(dp) :: finose(3,ab_mover%natom),binose(3,ab_mover%natom) 117 real(dp) :: vonose(3,ab_mover%natom),hnose(3,ab_mover%natom) 118 real(dp),allocatable,save :: fcart_m(:,:),fcart_mold(:,:) 119 real(dp) :: strten(6) 120 character(len=500) :: message 121 122 !*************************************************************************** 123 !Beginning of executable session 124 !*************************************************************************** 125 126 if(iexit/=0)then 127 ABI_SFREE(fcart_m) 128 ABI_SFREE(fcart_mold) 129 return 130 end if 131 132 !write(std_out,*) 'nose 01' 133 !########################################################## 134 !### 01. Allocate the arrays fcart_m and fcart_mold 135 136 if(itime==1)then 137 ABI_SFREE(fcart_m) 138 ABI_SFREE(fcart_mold) 139 end if 140 141 if(.not.allocated(fcart_m)) then 142 ABI_MALLOC(fcart_m,(3,ab_mover%natom)) 143 end if 144 if(.not.allocated(fcart_mold)) then 145 ABI_MALLOC(fcart_mold,(3,ab_mover%natom)) 146 end if 147 148 !write(std_out,*) 'nose 02' 149 !########################################################## 150 !### 02. Obtain the present values from the history 151 152 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 153 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 154 155 fcart(:,:)=hist%fcart(:,:,hist%ihist) 156 strten(:)=hist%strten(:,hist%ihist) 157 vel(:,:)=hist%vel(:,:,hist%ihist) 158 etotal=hist%etot(hist%ihist) 159 160 write(std_out,*) 'RPRIMD' 161 do ii=1,3 162 write(std_out,*) rprimd(:,ii) 163 end do 164 write(std_out,*) 'RMET' 165 do ii=1,3 166 write(std_out,*) rmet(ii,:) 167 end do 168 169 !write(std_out,*) 'nose 03' 170 !########################################################## 171 !### 03. Fill the vectors vin and vout 172 173 !write(std_out,*) 'nose 04' 174 !########################################################## 175 !### 04. Initialize or update the hessian matrix 176 177 !write(std_out,*) 'nose 05' 178 !########################################################## 179 !### 05. Compute the next values 180 181 !The temperature is linear between initial and final values 182 !It is here converted from Kelvin to Hartree (kb_HaK) 183 ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK 184 185 !%%% NOSE DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 186 187 acell_next(:)=acell(:) 188 ucvol_next=ucvol 189 rprimd_next(:,:)=rprimd(:,:) 190 191 if(itime==1)then 192 snose=0.0_dp 193 xi_nose=0.0_dp 194 ! Compute twice the kinetic energy of the system, called v2nose 195 v2nose=0.0_dp 196 do kk=1,ab_mover%natom 197 do jj=1,3 198 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 199 end do 200 end do 201 if (zDEBUG)then 202 write(std_out,*) 'itime ntime KTEMP=',itime-1,ntime-1,ktemp 203 write(std_out,*) 'V2NOSE=',v2nose 204 write (std_out,*) 'VEL' 205 do kk=1,ab_mover%natom 206 write (std_out,*) vel(:,kk) 207 end do 208 end if 209 210 ! If there is no kinetic energy, use a random initial velocity 211 if (v2nose<=v2tol) then 212 v2nose=0.0_dp 213 do kk=1,ab_mover%natom 214 do jj=1,3 215 ! Uniform random returns a uniform random deviate between 0.0 216 ! and 1.0 217 ! if it were always 0 or 1, then the following expression 218 ! would give the requested temperature 219 vel(jj,kk)=(1.0_dp-2.0_dp*uniformrandom(idum))*& 220 & sqrt( (ab_mover%mdtemp(1)) * kb_HaK / ab_mover%amass(kk) ) 221 ! Recompute v2nose 222 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 223 if (zDEBUG)then 224 write(std_out,*) 'jj kk vel(jj,kk)=',jj,kk,vel(jj,kk) 225 write(std_out,*) 'jj kk V2NOSE=',jj,kk,v2nose 226 end if 227 end do 228 end do 229 end if 230 write(std_out,*) 'V2NOSE=',v2nose 231 232 ! Now, rescale the velocities to give the proper temperature 233 rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose) 234 write(std_out,*) 'RESCALE_VEL=',rescale_vel 235 vel(:,:)=vel(:,:)*rescale_vel 236 ! Recompute v2nose with the rescaled velocities 237 v2nose=0.0_dp 238 do kk=1,ab_mover%natom 239 do jj=1,3 240 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 241 end do 242 end do 243 write(message, '(a)' )& 244 & ' Rescaling or initializing velocities to initial temperature' 245 call wrtout(std_out,message,'COLL') 246 call wrtout(std_out,message,'COLL') 247 write(message, '(2(a,es22.14))' )& 248 & ' --- Scaling factor : ',rescale_vel,& 249 & ' Asked T (K) ',ab_mover%mdtemp(1) 250 call wrtout(std_out,message,'COLL') 251 call wrtout(std_out,message,'COLL') 252 write(message, '(a,es22.14)' )& 253 & ' --- Effective temperature',v2nose/(3.0_dp*ab_mover%natom*kb_HaK) 254 call wrtout(std_out,message,'COLL') 255 call wrtout(std_out,message,'COLL') 256 end if 257 258 do kk=1,ab_mover%natom 259 do jj=1,3 260 fcart_m(jj,kk)=fcart(jj,kk)/ab_mover%amass(kk) 261 end do 262 end do 263 264 !First step of velocity verlet algorithm 265 gnose=3*ab_mover%natom 266 267 !Convert input xred (reduced coordinates) to xcart (cartesian) 268 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 269 270 !Calculate nose-hoover force on atoms 271 !If first iteration, no old force are available, so use present 272 !forces 273 if (itime==1) fcart_mold(:,:)=fcart_m(:,:) 274 275 if (zDEBUG)then 276 write (std_out,*) 'FCART_MOLD' 277 do kk=1,ab_mover%natom 278 write (std_out,*) fcart_mold(:,kk) 279 end do 280 write (std_out,*) 'FCART_M' 281 do kk=1,ab_mover%natom 282 write (std_out,*) fcart_m(:,kk) 283 end do 284 end if 285 286 finose(:,:)=fcart_mold(:,:)-xi_nose*vel(:,:) 287 xcart(:,:)=xcart(:,:)+ab_mover%dtion*(vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp) 288 289 !Convert back to xred (reduced coordinates) 290 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 291 292 if (zDEBUG)then 293 write (std_out,*) 'VEL' 294 do kk=1,ab_mover%natom 295 write (std_out,*) vel(:,kk) 296 end do 297 end if 298 299 !Calculate v2nose 300 v2nose=0.0_dp 301 do kk=1,ab_mover%natom 302 do jj=1,3 303 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 304 end do 305 end do 306 vel(:,:)=vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp 307 308 if (zDEBUG)then 309 write(std_out,*) 'NOSE BEFORE' 310 write(std_out,*) 'FSNOSE=',fsnose 311 write(std_out,*) 'SNOSE=',snose 312 write(std_out,*) 'XI_NOSE=',xi_nose 313 write (std_out,*) 'VEL' 314 do kk=1,ab_mover%natom 315 write (std_out,*) vel(:,kk) 316 end do 317 write (std_out,*) 'NOSEINERT',ab_mover%noseinert 318 end if 319 320 !Update thermostat 321 fsnose=(v2nose-gnose*ktemp)/ab_mover%noseinert 322 snose=snose+ab_mover%dtion*(xi_nose+ab_mover%dtion*fsnose/2.0_dp) 323 xi_nose=xi_nose+ab_mover%dtion*fsnose/2.0_dp 324 if (zDEBUG)then 325 write(std_out,*) 'NOSE AFTER' 326 write(std_out,*) 'FSNOSE=',fsnose 327 write(std_out,*) 'SNOSE=',snose 328 write(std_out,*) 'XI_NOSE=',xi_nose 329 write (std_out,*) 'VEL' 330 do kk=1,ab_mover%natom 331 write (std_out,*) vel(:,kk) 332 end do 333 end if 334 335 !Second step of the velocity Verlet algorithm, uses the 'new forces' 336 !Calculate v2nose 337 v2nose=0.0_dp 338 do kk=1,ab_mover%natom 339 do jj=1,3 340 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 341 end do 342 end do 343 vel_temp(:,:)=vel(:,:) 344 345 if (zDEBUG)then 346 write(std_out,*) 'V2NOSE=',v2nose 347 write (std_out,*) 'VEL' 348 do kk=1,ab_mover%natom 349 write (std_out,*) vel(:,kk) 350 end do 351 write (std_out,*) 'Starting Newton Raphson' 352 end if 353 354 xin_nose=xi_nose 355 356 !Start Newton-Raphson loop 357 ready=.false. 358 do while (.not.ready) 359 xio=xin_nose 360 delxi=0.0D0 361 vonose(:,:)=vel_temp(:,:) 362 hnose(:,:)=-ab_mover%dtion/2.0_dp*(fcart_m(:,:)-xio*vonose(:,:))-& 363 & (vel(:,:)-vonose(:,:)) 364 do kk=1,ab_mover%natom 365 do jj=1,3 366 binose(jj,kk)=vonose(jj,kk)*ab_mover%dtion/ab_mover%noseinert*ab_mover%amass(kk) ! a verifier 367 delxi=delxi+hnose(jj,kk)*binose(jj,kk) 368 end do 369 end do 370 dnose=-(xio*ab_mover%dtion/2.0D0+1.0D0) 371 delxi=delxi-dnose*((-v2nose+gnose*ktemp)*ab_mover%dtion/2.0_dp/ & 372 & ab_mover%noseinert-(xi_nose-xio)) 373 delxi=delxi/(-ab_mover%dtion*ab_mover%dtion/2.0_dp*v2nose/ab_mover%noseinert+dnose) 374 375 ! hzeronose=-(xio-xi_nose-(v2nose-gnose*ktemp) 376 ! *dtion/(2.0_dp*ab_mover%noseinert) ) 377 ! cibinose=-v2nose*dtion*dtion/(2.0_dp*ab_mover%noseinert) 378 ! delxi=(delxi+hzeronose*dnose)/(dnose+cibinose) 379 380 ! DEBUG 381 ! write(message, '(a,es22.14)' )' after delxi',delxi 382 ! call wrtout(std_out,message,'COLL') 383 ! call wrtout(std_out,message,'COLL') 384 ! ENDDEBUG 385 v2nose=0.0_dp 386 387 vel_temp(:,:)=vel_temp(:,:)+& 388 & (hnose+ab_mover%dtion/2.0_dp*vonose(:,:)*delxi)/dnose 389 do kk=1,ab_mover%natom 390 do jj=1,3 391 v2nose=v2nose+vel_temp(jj,kk)*& 392 & vel_temp(jj,kk)*ab_mover%amass(kk) 393 end do 394 end do 395 ! New guess for xi 396 xin_nose=xio+delxi 397 398 ! zDEBUG 399 ! write(message, '(a,es22.14)' )' v2nose=',v2nose 400 ! call wrtout(std_out,message,'COLL') 401 ! call wrtout(std_out,message,'COLL') 402 ! ENDDEBUG 403 404 ready=.true. 405 ! Test for convergence 406 kk=0 407 jj=1 408 do while((kk<=ab_mover%natom).and.(jj<=3).and.ready) 409 kk=kk+1 410 if (kk>ab_mover%natom) then 411 kk=1 412 jj=jj+1 413 end if 414 if ((kk<=ab_mover%natom) .and.(jj<=3)) then 415 if (abs(vel_temp(jj,kk))<1.0d-50)& 416 & vel_temp(jj,kk)=1.0d-50 417 if (abs((vel_temp(jj,kk)-vonose(jj,kk))& 418 & /vel_temp(jj,kk))>nosetol) ready=.false. 419 else 420 if (xin_nose<1.0d-50) xin_nose=1.0d-50 421 if (abs((xin_nose-xio)/xin_nose)>nosetol) ready=.false. 422 end if 423 end do ! end of while 424 425 ! Enddo ready 426 end do 427 428 !Update velocities to converged value 429 vel(:,:)=vel_temp(:,:) 430 write(message, '(a,es14.7)' )' converged velocities for T=',ktemp 431 call wrtout(std_out,message,'COLL') 432 433 if (zDEBUG)then 434 write (std_out,*) 'Final Values for NOSE' 435 write (std_out,*) 'VEL' 436 do kk=1,ab_mover%natom 437 write (std_out,*) vel(:,kk) 438 end do 439 write (std_out,*) 'XCART' 440 do kk=1,ab_mover%natom 441 write (std_out,*) xcart(:,kk) 442 end do 443 end if 444 445 !Update thermostat 446 xi_nose=xin_nose 447 xcart_next(:,:)=xcart(:,:) 448 !Convert back to xred_next (reduced coordinates) 449 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 450 !Store 'new force' as 'old force' 451 fcart_mold(:,:)=fcart_m(:,:) 452 453 !write(std_out,*) 'nose 06' 454 !########################################################## 455 !### 06. Update the history with the prediction 456 457 !Increase indexes 458 hist%ihist=abihist_findIndex(hist,+1) 459 460 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 461 hist%vel(:,:,hist%ihist)=vel(:,:) 462 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 463 464 end subroutine pred_nose