TABLE OF CONTENTS
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
COPYRIGHT
Copyright (C) 1998-2018 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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
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
NOTES
PARENTS
mover
CHILDREN
hist2var,metric,var2hist,wrtout,xcart2xred,xred2xcart
SOURCE
59 #if defined HAVE_CONFIG_H 60 #include "config.h" 61 #endif 62 63 #include "abi_common.h" 64 65 subroutine pred_nose(ab_mover,hist,itime,ntime,zDEBUG,iexit) 66 67 use defs_basis 68 use m_profiling_abi 69 use m_abimover 70 use m_abihist 71 72 use m_numeric_tools, only : uniformrandom 73 use m_geometry, only : xcart2xred, xred2xcart, metric 74 75 !This section has been created automatically by the script Abilint (TD). 76 !Do not modify the following lines by hand. 77 #undef ABI_FUNC 78 #define ABI_FUNC 'pred_nose' 79 use interfaces_14_hidewrite 80 !End of the abilint section 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) :: fred_corrected(3,ab_mover%natom) 114 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 115 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 116 real(dp) :: vel(3,ab_mover%natom),vel_temp(3,ab_mover%natom) 117 real(dp) :: finose(3,ab_mover%natom),binose(3,ab_mover%natom) 118 real(dp) :: vonose(3,ab_mover%natom),hnose(3,ab_mover%natom) 119 real(dp),allocatable,save :: fcart_m(:,:),fcart_mold(:,:) 120 real(dp) :: strten(6) 121 character(len=500) :: message 122 123 !*************************************************************************** 124 !Beginning of executable session 125 !*************************************************************************** 126 127 if(iexit/=0)then 128 if (allocated(fcart_m)) then 129 ABI_DEALLOCATE(fcart_m) 130 end if 131 if (allocated(fcart_mold)) then 132 ABI_DEALLOCATE(fcart_mold) 133 end if 134 return 135 end if 136 137 !write(std_out,*) 'nose 01' 138 !########################################################## 139 !### 01. Allocate the arrays fcart_m and fcart_mold 140 141 if(itime==1)then 142 if (allocated(fcart_m)) then 143 ABI_DEALLOCATE(fcart_m) 144 end if 145 if (allocated(fcart_mold)) then 146 ABI_DEALLOCATE(fcart_mold) 147 end if 148 end if 149 150 if(.not.allocated(fcart_m)) then 151 ABI_ALLOCATE(fcart_m,(3,ab_mover%natom)) 152 end if 153 if(.not.allocated(fcart_mold)) then 154 ABI_ALLOCATE(fcart_mold,(3,ab_mover%natom)) 155 end if 156 157 !write(std_out,*) 'nose 02' 158 !########################################################## 159 !### 02. Obtain the present values from the history 160 161 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 162 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 163 164 fcart(:,:)=hist%fcart(:,:,hist%ihist) 165 strten(:)=hist%strten(:,hist%ihist) 166 vel(:,:)=hist%vel(:,:,hist%ihist) 167 etotal=hist%etot(hist%ihist) 168 169 write(std_out,*) 'RPRIMD' 170 do ii=1,3 171 write(std_out,*) rprimd(:,ii) 172 end do 173 write(std_out,*) 'RMET' 174 do ii=1,3 175 write(std_out,*) rmet(ii,:) 176 end do 177 178 !Get rid of mean force on whole unit cell, but only if no 179 !generalized constraints are in effect 180 ! call fcart2fred(fcart,fred_corrected,rprimd,ab_mover%natom) 181 ! if(ab_mover%nconeq==0)then 182 ! amass_tot=sum(ab_mover%amass(:)) 183 ! do ii=1,3 184 ! if (ii/=3.or.ab_mover%jellslab==0) then 185 ! favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom) 186 ! fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot 187 ! end if 188 ! end do 189 ! end if 190 191 !write(std_out,*) 'nose 03' 192 !########################################################## 193 !### 03. Fill the vectors vin and vout 194 195 !write(std_out,*) 'nose 04' 196 !########################################################## 197 !### 04. Initialize or update the hessian matrix 198 199 !write(std_out,*) 'nose 05' 200 !########################################################## 201 !### 05. Compute the next values 202 203 !The temperature is linear between initial and final values 204 !It is here converted from Kelvin to Hartree (kb_HaK) 205 ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK 206 207 !%%% NOSE DYNAMICS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 208 209 acell_next(:)=acell(:) 210 ucvol_next=ucvol 211 rprimd_next(:,:)=rprimd(:,:) 212 213 if(itime==1)then 214 snose=0.0_dp 215 xi_nose=0.0_dp 216 ! Compute twice the kinetic energy of the system, called v2nose 217 v2nose=0.0_dp 218 do kk=1,ab_mover%natom 219 do jj=1,3 220 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 221 end do 222 end do 223 if (zDEBUG)then 224 write(std_out,*) 'itime ntime KTEMP=',itime-1,ntime-1,ktemp 225 write(std_out,*) 'V2NOSE=',v2nose 226 write (std_out,*) 'VEL' 227 do kk=1,ab_mover%natom 228 write (std_out,*) vel(:,kk) 229 end do 230 end if 231 232 ! If there is no kinetic energy, use a random initial velocity 233 if (v2nose<=v2tol) then 234 v2nose=0.0_dp 235 do kk=1,ab_mover%natom 236 do jj=1,3 237 ! Uniform random returns a uniform random deviate between 0.0 238 ! and 1.0 239 ! if it were always 0 or 1, then the following expression 240 ! would give the requested temperature 241 vel(jj,kk)=(1.0_dp-2.0_dp*uniformrandom(idum))*& 242 & sqrt( (ab_mover%mdtemp(1)) * kb_HaK / ab_mover%amass(kk) ) 243 ! Recompute v2nose 244 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 245 if (zDEBUG)then 246 write(std_out,*) 'jj kk vel(jj,kk)=',jj,kk,vel(jj,kk) 247 write(std_out,*) 'jj kk V2NOSE=',jj,kk,v2nose 248 end if 249 end do 250 end do 251 end if 252 write(std_out,*) 'V2NOSE=',v2nose 253 254 ! Now, rescale the velocities to give the proper temperature 255 rescale_vel=sqrt(3.0_dp*ab_mover%natom*(ab_mover%mdtemp(1))*kb_HaK/v2nose) 256 write(std_out,*) 'RESCALE_VEL=',rescale_vel 257 vel(:,:)=vel(:,:)*rescale_vel 258 ! Recompute v2nose with the rescaled velocities 259 v2nose=0.0_dp 260 do kk=1,ab_mover%natom 261 do jj=1,3 262 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 263 end do 264 end do 265 write(message, '(a)' )& 266 & ' Rescaling or initializing velocities to initial temperature' 267 call wrtout(std_out,message,'COLL') 268 call wrtout(std_out,message,'COLL') 269 write(message, '(2(a,es22.14))' )& 270 & ' --- Scaling factor : ',rescale_vel,& 271 & ' Asked T (K) ',ab_mover%mdtemp(1) 272 call wrtout(std_out,message,'COLL') 273 call wrtout(std_out,message,'COLL') 274 write(message, '(a,es22.14)' )& 275 & ' --- Effective temperature',v2nose/(3.0_dp*ab_mover%natom*kb_HaK) 276 call wrtout(std_out,message,'COLL') 277 call wrtout(std_out,message,'COLL') 278 end if 279 280 do kk=1,ab_mover%natom 281 do jj=1,3 282 fcart_m(jj,kk)=fcart(jj,kk)/ab_mover%amass(kk) 283 end do 284 end do 285 286 !First step of velocity verlet algorithm 287 gnose=3*ab_mover%natom 288 289 !Convert input xred (reduced coordinates) to xcart (cartesian) 290 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 291 292 !Calculate nose-hoover force on atoms 293 !If first iteration, no old force are available, so use present 294 !forces 295 if (itime==1) fcart_mold(:,:)=fcart_m(:,:) 296 297 if (zDEBUG)then 298 write (std_out,*) 'FCART_MOLD' 299 do kk=1,ab_mover%natom 300 write (std_out,*) fcart_mold(:,kk) 301 end do 302 write (std_out,*) 'FCART_M' 303 do kk=1,ab_mover%natom 304 write (std_out,*) fcart_m(:,kk) 305 end do 306 end if 307 308 finose(:,:)=fcart_mold(:,:)-xi_nose*vel(:,:) 309 xcart(:,:)=xcart(:,:)+ab_mover%dtion*(vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp) 310 311 !Convert back to xred (reduced coordinates) 312 call xcart2xred(ab_mover%natom,rprimd,xcart,xred) 313 314 if (zDEBUG)then 315 write (std_out,*) 'VEL' 316 do kk=1,ab_mover%natom 317 write (std_out,*) vel(:,kk) 318 end do 319 end if 320 321 !Calculate v2nose 322 v2nose=0.0_dp 323 do kk=1,ab_mover%natom 324 do jj=1,3 325 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 326 end do 327 end do 328 vel(:,:)=vel(:,:)+ab_mover%dtion*finose(:,:)/2.0_dp 329 330 if (zDEBUG)then 331 write(std_out,*) 'NOSE BEFORE' 332 write(std_out,*) 'FSNOSE=',fsnose 333 write(std_out,*) 'SNOSE=',snose 334 write(std_out,*) 'XI_NOSE=',xi_nose 335 write (std_out,*) 'VEL' 336 do kk=1,ab_mover%natom 337 write (std_out,*) vel(:,kk) 338 end do 339 write (std_out,*) 'NOSEINERT',ab_mover%noseinert 340 end if 341 342 !Update thermostat 343 fsnose=(v2nose-gnose*ktemp)/ab_mover%noseinert 344 snose=snose+ab_mover%dtion*(xi_nose+ab_mover%dtion*fsnose/2.0_dp) 345 xi_nose=xi_nose+ab_mover%dtion*fsnose/2.0_dp 346 if (zDEBUG)then 347 write(std_out,*) 'NOSE AFTER' 348 write(std_out,*) 'FSNOSE=',fsnose 349 write(std_out,*) 'SNOSE=',snose 350 write(std_out,*) 'XI_NOSE=',xi_nose 351 write (std_out,*) 'VEL' 352 do kk=1,ab_mover%natom 353 write (std_out,*) vel(:,kk) 354 end do 355 end if 356 357 !Second step of the velocity Verlet algorithm, uses the 'new forces' 358 !Calculate v2nose 359 v2nose=0.0_dp 360 do kk=1,ab_mover%natom 361 do jj=1,3 362 v2nose=v2nose+vel(jj,kk)*vel(jj,kk)*ab_mover%amass(kk) 363 end do 364 end do 365 vel_temp(:,:)=vel(:,:) 366 367 if (zDEBUG)then 368 write(std_out,*) 'V2NOSE=',v2nose 369 write (std_out,*) 'VEL' 370 do kk=1,ab_mover%natom 371 write (std_out,*) vel(:,kk) 372 end do 373 write (std_out,*) 'Starting Newton Raphson' 374 end if 375 376 xin_nose=xi_nose 377 378 !Start Newton-Raphson loop 379 ready=.false. 380 do while (.not.ready) 381 xio=xin_nose 382 delxi=0.0D0 383 vonose(:,:)=vel_temp(:,:) 384 hnose(:,:)=-ab_mover%dtion/2.0_dp*(fcart_m(:,:)-xio*vonose(:,:))-& 385 & (vel(:,:)-vonose(:,:)) 386 do kk=1,ab_mover%natom 387 do jj=1,3 388 binose(jj,kk)=vonose(jj,kk)*ab_mover%dtion/ab_mover%noseinert*ab_mover%amass(kk) ! a verifier 389 delxi=delxi+hnose(jj,kk)*binose(jj,kk) 390 end do 391 end do 392 dnose=-(xio*ab_mover%dtion/2.0D0+1.0D0) 393 delxi=delxi-dnose*((-v2nose+gnose*ktemp)*ab_mover%dtion/2.0_dp/ & 394 & ab_mover%noseinert-(xi_nose-xio)) 395 delxi=delxi/(-ab_mover%dtion*ab_mover%dtion/2.0_dp*v2nose/ab_mover%noseinert+dnose) 396 397 ! hzeronose=-(xio-xi_nose-(v2nose-gnose*ktemp) 398 ! *dtion/(2.0_dp*ab_mover%noseinert) ) 399 ! cibinose=-v2nose*dtion*dtion/(2.0_dp*ab_mover%noseinert) 400 ! delxi=(delxi+hzeronose*dnose)/(dnose+cibinose) 401 402 ! DEBUG 403 ! write(message, '(a,es22.14)' )' after delxi',delxi 404 ! call wrtout(std_out,message,'COLL') 405 ! call wrtout(std_out,message,'COLL') 406 ! ENDDEBUG 407 v2nose=0.0_dp 408 409 vel_temp(:,:)=vel_temp(:,:)+& 410 & (hnose+ab_mover%dtion/2.0_dp*vonose(:,:)*delxi)/dnose 411 do kk=1,ab_mover%natom 412 do jj=1,3 413 v2nose=v2nose+vel_temp(jj,kk)*& 414 & vel_temp(jj,kk)*ab_mover%amass(kk) 415 end do 416 end do 417 ! New guess for xi 418 xin_nose=xio+delxi 419 420 ! zDEBUG 421 ! write(message, '(a,es22.14)' )' v2nose=',v2nose 422 ! call wrtout(std_out,message,'COLL') 423 ! call wrtout(std_out,message,'COLL') 424 ! ENDDEBUG 425 426 ready=.true. 427 ! Test for convergence 428 kk=0 429 jj=1 430 do while((kk<=ab_mover%natom).and.(jj<=3).and.ready) 431 kk=kk+1 432 if (kk>ab_mover%natom) then 433 kk=1 434 jj=jj+1 435 end if 436 if ((kk<=ab_mover%natom) .and.(jj<=3)) then 437 if (abs(vel_temp(jj,kk))<1.0d-50)& 438 & vel_temp(jj,kk)=1.0d-50 439 if (abs((vel_temp(jj,kk)-vonose(jj,kk))& 440 & /vel_temp(jj,kk))>nosetol) ready=.false. 441 else 442 if (xin_nose<1.0d-50) xin_nose=1.0d-50 443 if (abs((xin_nose-xio)/xin_nose)>nosetol) ready=.false. 444 end if 445 end do ! end of while 446 447 ! Enddo ready 448 end do 449 450 !Update velocities to converged value 451 vel(:,:)=vel_temp(:,:) 452 write(message, '(a,es14.7)' )' converged velocities for T=',ktemp 453 call wrtout(std_out,message,'COLL') 454 455 if (zDEBUG)then 456 write (std_out,*) 'Final Values for NOSE' 457 write (std_out,*) 'VEL' 458 do kk=1,ab_mover%natom 459 write (std_out,*) vel(:,kk) 460 end do 461 write (std_out,*) 'XCART' 462 do kk=1,ab_mover%natom 463 write (std_out,*) xcart(:,kk) 464 end do 465 end if 466 467 !Update thermostat 468 xi_nose=xin_nose 469 xcart_next(:,:)=xcart(:,:) 470 !Convert back to xred_next (reduced coordinates) 471 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 472 !Store 'new force' as 'old force' 473 fcart_mold(:,:)=fcart_m(:,:) 474 475 !write(std_out,*) 'nose 06' 476 !########################################################## 477 !### 06. Update the history with the prediction 478 479 !Increase indexes 480 hist%ihist=abihist_findIndex(hist,+1) 481 482 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 483 hist%vel(:,:,hist%ihist)=vel(:,:) 484 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 485 486 end subroutine pred_nose