TABLE OF CONTENTS
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)], as worked out by Minary et al [J. Chem. Phys. 188, 2510 (2003)]. 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)), the final temperature (mdtemp(2)), and the friction coefficient (friction).
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,var2hist,wrtout,xcart2xred,xred2xcart
SOURCE
52 #if defined HAVE_CONFIG_H 53 #include "config.h" 54 #endif 55 56 #include "abi_common.h" 57 58 subroutine pred_isokinetic(ab_mover,hist,itime,ntime,zDEBUG,iexit) 59 60 use m_profiling_abi 61 use defs_basis 62 use m_abimover 63 use m_abihist 64 65 use m_numeric_tools, only : uniformrandom 66 use m_geometry, only : xcart2xred, xred2xcart 67 68 !This section has been created automatically by the script Abilint (TD). 69 !Do not modify the following lines by hand. 70 #undef ABI_FUNC 71 #define ABI_FUNC 'pred_isokinetic' 72 use interfaces_14_hidewrite 73 !End of the abilint section 74 75 implicit none 76 77 !Arguments ------------------------------------ 78 !scalars 79 integer,intent(in) :: itime 80 integer,intent(in) :: ntime 81 integer,intent(in) :: iexit 82 logical,intent(in) :: zDEBUG 83 type(abimover),intent(in) :: ab_mover 84 type(abihist),intent(inout) :: hist 85 86 !Local variables------------------------------- 87 !scalars 88 integer :: kk,iatom,idim,idum=5,nfirst,ifirst 89 real(dp) :: a,as,b,sqb,s,s1,s2,scdot,sigma2,vtest,v2gauss 90 real(dp),parameter :: v2tol=tol8 91 real(dp) :: etotal,rescale_vel 92 character(len=5000) :: message 93 !arrays 94 real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:) 95 96 real(dp) :: acell(3),rprimd(3,3) 97 real(dp) :: fcart(3,ab_mover%natom) 98 !real(dp) :: fred_corrected(3,ab_mover%natom) 99 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 100 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 101 real(dp) :: vel(3,ab_mover%natom) 102 real(dp) :: strten(6) 103 104 !*************************************************************************** 105 !Beginning of executable session 106 !*************************************************************************** 107 108 if(iexit/=0)then 109 if (allocated(fcart_m)) then 110 ABI_DEALLOCATE(fcart_m) 111 end if 112 if (allocated(vel_nexthalf)) then 113 ABI_DEALLOCATE(vel_nexthalf) 114 end if 115 return 116 end if 117 118 !write(std_out,*) 'isokinetic 01' 119 !########################################################## 120 !### 01. Debugging and Verbose 121 122 if(zDEBUG)then 123 write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),& 124 & 'Debugging and Verbose for pred_isokinetic',('-',kk=1,37) 125 write(std_out,*) 'ionmov: ',12 126 write(std_out,*) 'itime: ',itime 127 end if 128 129 !write(std_out,*) 'isokinetic 02' 130 !########################################################## 131 !### 02. Allocate the vectors vin, vout and hessian matrix 132 !### These arrays could be allocated from a previus 133 !### dataset that exit before itime==ntime 134 135 if(itime==1)then 136 if (allocated(fcart_m)) then 137 ABI_DEALLOCATE(fcart_m) 138 end if 139 if (allocated(vel_nexthalf)) then 140 ABI_DEALLOCATE(vel_nexthalf) 141 end if 142 end if 143 144 if (.not.allocated(fcart_m)) then 145 ABI_ALLOCATE(fcart_m,(3,ab_mover%natom)) 146 end if 147 if (.not.allocated(vel_nexthalf)) then 148 ABI_ALLOCATE(vel_nexthalf,(3,ab_mover%natom)) 149 end if 150 151 !write(std_out,*) 'isokinetic 03' 152 !########################################################## 153 !### 03. Obtain the present values from the history 154 155 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 156 157 fcart(:,:)=hist%fcart(:,:,hist%ihist) 158 strten(:) =hist%strten(:,hist%ihist) 159 vel(:,:) =hist%vel(:,:,hist%ihist) 160 etotal =hist%etot(hist%ihist) 161 162 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 163 164 if(zDEBUG)then 165 write (std_out,*) 'fcart:' 166 do kk=1,ab_mover%natom 167 write (std_out,*) fcart(:,kk) 168 end do 169 write (std_out,*) 'vel:' 170 do kk=1,ab_mover%natom 171 write (std_out,*) vel(:,kk) 172 end do 173 write (std_out,*) 'strten:' 174 write (std_out,*) strten(1:3),ch10,strten(4:6) 175 write (std_out,*) 'etotal:' 176 write (std_out,*) etotal 177 end if 178 179 !Get rid of mean force on whole unit cell, but only if no 180 !generalized constraints are in effect 181 ! call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom) 182 ! if(ab_mover%nconeq==0)then 183 ! amass_tot=sum(ab_mover%amass(:)) 184 ! do ii=1,3 185 ! if (ii/=3.or.ab_mover%jellslab==0) then 186 ! favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom) 187 ! fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot 188 ! end if 189 ! end do 190 ! end if 191 192 !write(std_out,*) 'isokinetic 04' 193 !########################################################## 194 !### 04. Second half-velocity (Only after the first itime) 195 196 if(itime>1) then 197 198 do iatom=1,ab_mover%natom 199 do idim=1,3 200 fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom) 201 end do 202 end do 203 204 ! Computation of vel(:,:) at the next positions 205 ! Computation of v2gauss 206 v2gauss=0.0_dp 207 do iatom=1,ab_mover%natom 208 do idim=1,3 209 v2gauss=v2gauss+& 210 & vel_nexthalf(idim,iatom)*vel_nexthalf(idim,iatom)*& 211 & ab_mover%amass(iatom) 212 end do 213 end do 214 215 ! Calcul de a et b (4.13 de Ref.1) 216 a=0.0_dp 217 b=0.0_dp 218 do iatom=1,ab_mover%natom 219 do idim=1,3 220 a=a+fcart_m(idim,iatom)*vel_nexthalf(idim,iatom)*ab_mover%amass(iatom) 221 b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom) 222 end do 223 end do 224 a=a/v2gauss 225 b=b/v2gauss 226 227 228 ! Calcul de s et scdot 229 sqb=sqrt(b) 230 as=sqb*ab_mover%dtion/2. 231 ! jmb 232 if ( as > 300.0 ) as=300.0 233 s1=cosh(as) 234 s2=sinh(as) 235 s=a*(s1-1.)/b+s2/sqb 236 scdot=a*s2/sqb+s1 237 vel(:,:)=(vel_nexthalf(:,:)+fcart_m(:,:)*s)/scdot 238 239 if (zDEBUG)then 240 write(std_out,*) 'Computation of the second half-velocity' 241 write(std_out,*) 'Cartesian forces per atomic mass (fcart_m):' 242 do kk=1,ab_mover%natom 243 write (std_out,*) fcart_m(:,kk) 244 end do 245 write(std_out,*) 'vel:' 246 do kk=1,ab_mover%natom 247 write (std_out,*) vel(:,kk) 248 end do 249 write(std_out,*) 'v2gauss:',v2gauss 250 write(std_out,*) 'a:',a 251 write(std_out,*) 'b:',b 252 write(std_out,*) 's:',s 253 write(std_out,*) 'scdot:',scdot 254 end if 255 256 end if ! (if itime>1) 257 258 !write(std_out,*) 'isokinetic 05' 259 !########################################################## 260 !### 05. First half-time (First cycle the loop is double) 261 262 if (itime==1) then 263 nfirst=2 264 else 265 nfirst=1 266 end if 267 268 do ifirst=1,nfirst 269 270 ! Application of Gauss' principle of least constraint according to Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, p.6102) 271 272 ! v2gauss is twice the kinetic energy 273 v2gauss=0.0_dp 274 do iatom=1,ab_mover%natom 275 do idim=1,3 276 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 277 end do 278 end do 279 280 ! If there is no kinetic energy 281 if (v2gauss<=v2tol.and.itime==1) then 282 ! Maxwell-Boltzman distribution 283 v2gauss=zero 284 vtest=zero 285 do iatom=1,ab_mover%natom 286 do idim=1,3 287 vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum)) 288 vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum))) 289 end do 290 end do 291 292 ! Get rid of center-of-mass velocity 293 s1=sum(ab_mover%amass(:)) 294 do idim=1,3 295 s2=sum(ab_mover%amass(:)*vel(idim,:)) 296 vel(idim,:)=vel(idim,:)-s2/s1 297 end do 298 299 ! Recompute v2gauss 300 do iatom=1,ab_mover%natom 301 do idim=1,3 302 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 303 vtest=vtest+vel(idim,iatom)/(3._dp*ab_mover%natom) 304 end do 305 end do 306 307 ! Now rescale the velocities to give the exact temperature 308 rescale_vel=sqrt(3._dp*ab_mover%natom*kb_HaK*ab_mover%mdtemp(1)/v2gauss) 309 vel(:,:)=vel(:,:)*rescale_vel 310 311 ! Recompute v2gauss with the rescaled velocities 312 v2gauss=zero 313 do iatom=1,ab_mover%natom 314 do idim=1,3 315 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 316 end do 317 end do 318 319 ! Compute the variance and print 320 sigma2=(v2gauss/(3._dp*ab_mover%natom)-ab_mover%amass(1)*vtest**2)/kb_HaK 321 322 end if 323 324 do iatom=1,ab_mover%natom 325 do idim=1,3 326 fcart_m(idim,iatom)=fcart(idim,iatom)/ab_mover%amass(iatom) 327 end do 328 end do 329 330 if (zDEBUG)then 331 write(std_out,*) 'Calculation first half-velocity ' 332 write (std_out,*) 'vel:' 333 do kk=1,ab_mover%natom 334 write (std_out,*) vel(:,kk) 335 end do 336 write (std_out,*) 'xcart:' 337 do kk=1,ab_mover%natom 338 write (std_out,*) xcart(:,kk) 339 end do 340 write (std_out,*) 'xred:' 341 do kk=1,ab_mover%natom 342 write (std_out,*) xred(:,kk) 343 end do 344 write (std_out,*) 'fcart_m' 345 do kk=1,ab_mover%natom 346 write (std_out,*) fcart_m(:,kk) 347 end do 348 write(std_out,*) 's2',s2 349 write(std_out,*) 'v2gauss',v2gauss 350 write(std_out,*) 'sigma2',sigma2 351 352 write(message, '(a)' )& 353 & ' --- Rescaling or initializing velocities to initial temperature' 354 call wrtout(std_out,message,'COLL') 355 write(message, '(a,d12.5,a,D12.5)' )& 356 & ' --- Scaling factor :',rescale_vel,' Asked T (K) ',ab_mover%mdtemp(1) 357 call wrtout(std_out,message,'COLL') 358 write(message, '(a,d12.5,a,D12.5)' )& 359 & ' --- Effective temperature',v2gauss/(3*ab_mover%natom*kb_HaK),' From variance', sigma2 360 call wrtout(std_out,message,'COLL') 361 end if 362 363 ! Convert input xred (reduced coordinates) to xcart (cartesian) 364 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 365 366 if(itime==1.and.ifirst==1) then 367 call wrtout(std_out,'if itime==1','COLL') 368 vel_nexthalf(:,:)=vel(:,:) 369 xcart_next(:,:)=xcart(:,:) 370 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 371 xred=xred_next 372 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 373 end if 374 375 end do 376 377 !Computation of vel_nexthalf (4.16 de Ref.1) 378 !Computation of a and b (4.13 de Ref.1) 379 a=0.0_dp 380 b=0.0_dp 381 do iatom=1,ab_mover%natom 382 do idim=1,3 383 a=a+fcart_m(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 384 b=b+fcart_m(idim,iatom)*fcart_m(idim,iatom)*ab_mover%amass(iatom) 385 end do 386 end do 387 a=a/v2gauss+tol20 388 b=b/v2gauss+tol20 389 !Computation of s and scdot 390 sqb=sqrt(b)+tol20 391 as=sqb*ab_mover%dtion/2. 392 ! jmb 393 if ( as > 300.0 ) as=300.0 394 s1=cosh(as) 395 s2=sinh(as) 396 s=a*(s1-1.)/b+s2/sqb 397 scdot=a*s2/sqb+s1 398 vel_nexthalf(:,:)=(vel(:,:)+fcart_m(:,:)*s)/scdot 399 400 !Computation of the next positions 401 xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion 402 403 if (zDEBUG)then 404 write(std_out,*) 'a:',a 405 write(std_out,*) 'b:',b 406 write(std_out,*) 's:',s 407 write(std_out,*) 'scdot:',scdot 408 end if 409 410 !Convert back to xred (reduced coordinates) 411 412 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 413 414 !write(std_out,*) 'isokinetic 06' 415 !########################################################## 416 !### 06. Update the history with the prediction 417 418 xcart=xcart_next 419 xred=xred_next 420 421 !increment the ihist 422 hist%ihist = abihist_findIndex(hist,+1) 423 424 !Fill the history with the variables 425 !xred, acell, rprimd, vel 426 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 427 hist%vel(:,:,hist%ihist)=vel(:,:) 428 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 429 430 if(zDEBUG)then 431 write (std_out,*) 'fcart:' 432 do kk=1,ab_mover%natom 433 write (std_out,*) fcart(:,kk) 434 end do 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,*) 'strten:' 440 write (std_out,*) strten(1:3),ch10,strten(4:6) 441 write (std_out,*) 'etotal:' 442 write (std_out,*) etotal 443 end if 444 445 if (.false.) write(std_out,*) ntime 446 447 end subroutine pred_isokinetic