TABLE OF CONTENTS
ABINIT/m_pred_lotf [ Modules ]
NAME
m_pred_lotf
FUNCTION
Contains the predictor for LOTF (ionmov==23)
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 . For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
SOURCE
17 #if defined HAVE_CONFIG_H 18 #include "config.h" 19 #endif 20 21 #include "abi_common.h" 22 23 24 module m_pred_lotf 25 26 use m_abicore 27 use defs_basis 28 use m_abimover 29 use m_abihist 30 use lotfpath 31 32 use m_numeric_tools, only : uniformrandom 33 use m_geometry, only : xcart2xred 34 35 implicit none 36 37 private 38 39 public :: pred_lotf 40 41 CONTAINS !===========================================================
ABINIT/m_pred_lotf/pred_lotf [ Functions ]
NAME
pred_lotf
FUNCTION
Ionmov predictors (12) Lotf ensemble molecular dynamics IONMOV 23: Lotf ensemble molecular dynamics.
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 . 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 icycle : Index of the cycle in the iteration zDEBUG : if true print some debugging information
OUTPUT
SIDE EFFECTS
hist <type(abihist)> : History of positions,forces acell, rprimd, stresses
NOTES
SOURCE
78 subroutine pred_lotf(ab_mover,hist,itime,icycle,zDEBUG,iexit) 79 80 use m_geometry, only : xred2xcart 81 implicit none 82 83 !Arguments ------------------------ 84 type(abimover),intent(in) :: ab_mover 85 type(abihist),intent(inout) :: hist 86 integer,intent(in) :: itime 87 integer,intent(in) :: icycle 88 integer,intent(in) :: iexit 89 logical,intent(in) :: zDEBUG 90 91 !Local variables------------------------------- 92 !scalars 93 integer :: kk,iatom,idim,idum=5 94 real(dp) :: v2gauss 95 real(dp),parameter :: v2tol=tol8 96 real(dp) :: etotal 97 character(len=5000) :: message 98 logical :: do_extrap,do_interp,do_first 99 !arrays 100 real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:) 101 real(dp),allocatable,save :: xcart_old(:,:),vel_old(:,:) 102 103 real(dp) :: acell(3),rprimd(3,3),fcart(3,ab_mover%natom) 104 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 105 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 106 real(dp) :: vel(3,ab_mover%natom) 107 real(dp) :: strten(6) 108 109 !*************************************************************************** 110 !Beginning of executable session 111 !*************************************************************************** 112 113 if(iexit/=0)then 114 if (allocated(fcart_m)) then 115 ABI_FREE(fcart_m) 116 end if 117 if (allocated(vel_nexthalf)) then 118 ABI_FREE(vel_nexthalf) 119 end if 120 if (allocated(xcart_old)) then 121 ABI_FREE(xcart_old) 122 end if 123 if (allocated(vel_old)) then 124 ABI_FREE(vel_old) 125 end if 126 return 127 end if 128 129 !write(std_out,*) 'lotf 01' 130 !########################################################## 131 !### 01. Debugging and Verbose 132 133 if(zDEBUG)then 134 write(std_out,'(a,3a,40a,37a)') ch10,('-',kk=1,3),& 135 & 'Debugging and Verbose for pred_lotf',('-',kk=1,37) 136 write(std_out,*) 'ionmov: ',12 137 write(std_out,*) 'itime: ',itime 138 end if 139 140 !write(std_out,*) 'lotf 02' 141 !########################################################## 142 !### 02. Allocate the vectors vin, vout and hessian matrix 143 !### These arrays could be allocated from a previus 144 !### dataset that exit before itime==ntime 145 146 147 if (.not.allocated(fcart_m)) then 148 ABI_MALLOC(fcart_m,(3,ab_mover%natom)) 149 end if 150 if (.not.allocated(vel_nexthalf)) then 151 ABI_MALLOC(vel_nexthalf,(3,ab_mover%natom)) 152 end if 153 if (.not.allocated(vel_old)) then 154 ABI_MALLOC(vel_old,(3,ab_mover%natom)) 155 end if 156 if (.not.allocated(xcart_old)) then 157 ABI_MALLOC(xcart_old,(3,ab_mover%natom)) 158 end if 159 160 161 !write(std_out,*) 'lotf 03' 162 !########################################################## 163 !### 03. Set what Pred_lotf has to do 164 165 !--First step only the first time pred_lotf is called 166 do_first = (itime==1 .and. icycle==1) 167 168 !--Extrapolation is computed only is itime is a multiple of lotfvar%nitex 169 do_extrap = lotf_extrapolation(itime) 170 171 !--Interpolation is done at icycle==2 when extrapolation is active, 172 !--else at the first cycle 173 if(do_extrap) then 174 do_interp = (icycle==2) 175 else 176 do_interp = (icycle==1) 177 endif 178 179 !write(std_out,*) 'lotf 04' 180 !########################################################## 181 !### 04. Obtain the present values from the history 182 183 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 184 185 fcart(:,:)=hist%fcart(:,:,hist%ihist) 186 strten(:) =hist%strten(:,hist%ihist) 187 vel(:,:) =hist%vel(:,:,hist%ihist) 188 etotal =hist%etot(hist%ihist) 189 190 if(zDEBUG)then 191 write (std_out,*) 'fcart:' 192 do kk=1,ab_mover%natom 193 write (std_out,*) fcart(:,kk) 194 end do 195 write (std_out,*) 'vel:' 196 do kk=1,ab_mover%natom 197 write (std_out,*) vel(:,kk) 198 end do 199 write (std_out,*) 'strten:' 200 write (std_out,*) strten(1:3),ch10,strten(4:6) 201 write (std_out,*) 'etotal:' 202 write (std_out,*) etotal 203 end if 204 205 !write(std_out,*) 'lotf 05' 206 !########################################################## 207 !### 05. First half-time (First cycle the loop is double) 208 209 first_step: if(do_first) then 210 write(message, '(a)' ) ' --- LOTF INITIALIZATION---' 211 call wrtout(ab_out,message,'COLL') 212 call wrtout(std_out,message,'COLL') 213 214 !--convert from xred to xcart 215 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 216 217 !PRINT *,'HH-first-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 218 219 !--call the LOTF initialization 220 call init_lotf(itime,ab_mover%natom,acell,rprimd,xcart) 221 222 !--Application of Gauss' principle of least constraint according to 223 ! Fei Zhang's algorithm (J. Chem. Phys. 106, 1997, p.6102 [[cite:Zhang1997]]) 224 !--v2gauss is twice the kinetic energy 225 call vel_to_gauss(vel,ab_mover%amass,v2gauss) 226 227 !--If there is no kinetic energy 228 if(v2gauss<=v2tol) then 229 !--Maxwell-Boltzman distribution 230 v2gauss=zero 231 do iatom=1,ab_mover%natom 232 do idim=1,3 233 vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum)) 234 vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum))) 235 end do 236 end do 237 endif 238 239 !--Init alpha_end parameters 240 call fitclus(.true.,fcart(:,:),xcart,alpha_end,1,ab_mover%natom) 241 242 !--Forces/mass for the first step 243 forall(iatom = 1:ab_mover%natom) fcart_m(:,iatom) = fcart(:,iatom)/ab_mover%amass(iatom) 244 245 !print *,'b-end 1step',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2)) 246 !PRINT *,'HH-first-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 247 248 end if first_step 249 250 251 !--Do extrapolation: 252 ! when extrapolation is active for a fixed itime, mover procedure 253 ! will call pred_lotf twice, 254 ! 1-a) The final values of alpha are set as initial value 255 ! 1-b) The actual coordinates and velocities are stored in a local saved array 256 ! 1-c) I the first the coordinated of the atoms at the step itime+lotfvar%nitex 257 ! will be extrapolated and alpha_in updated 258 ! 2-a) The new final value of alpha_end (step itime+lotfvar%nitex) are computed 259 ! 2-b) Restore velocities and positions 260 extrapolation_lotf: if(do_extrap) then 261 select case(icycle) 262 case(1) 263 !-- 264 write(message, '(a)' )' LOTF : Extrapolation' 265 call wrtout(ab_out,message,'COLL') 266 call wrtout(std_out,message,'COLL') 267 268 !--Intial alpha values to do extrpolation correspond the last one 269 alpha_in = alpha_end 270 271 !--store the value of xcart 272 xcart_old = xcart 273 vel_old = vel 274 275 !PRINT *,'HH-extra-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 276 277 !--Compute new bond, and positions by extrpolation (xcartfit) at 278 ! the step=itime+lotfvar%nitex 279 call extrapolation_loop(itime,ab_mover%mdtemp(1),ab_mover%dtion,ab_mover%amass,xcart,& 280 & vel,xcart_next,vel_nexthalf) 281 282 !print *,'b-end extra',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2)) 283 284 !--Convert back to xred (reduced coordinates) 285 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 286 287 !PRINT *,'HH-extra-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 288 289 case(2) 290 !PRINT *,'HH-extra-3',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 291 292 call fitclus(.true.,fcart(:,:),xcart,alpha_end,1,ab_mover%natom) 293 !print *,'b-end fitclus',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2)) 294 295 !--restore the value of xcart 296 xcart = xcart_old 297 vel = vel_old 298 !PRINT *,'HH-extra-4',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 299 300 end select 301 endif extrapolation_lotf 302 303 304 !--Interpolation(pass here for extrapolation when icycle==2 else when icycle==1) 305 if(do_interp)then 306 307 !PRINT *,'HH-norma-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 308 309 !--Compute rescaled vel (center of mass speed and gaussian) 310 call vel_rescale(ab_mover%mdtemp(1),vel,ab_mover%amass,v2gauss) 311 312 !PRINT *,'HH-norma-2',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 313 314 !--VERLETVEL (interpolation) using the interpolated value of alpha 315 ! in the interval [alpha_in,alpha_end] 316 call lotf_interpolation(itime,ab_mover%dtion,v2gauss,ab_mover%amass,xcart,vel,fcart_m,xcart_next,vel_nexthalf) 317 318 !PRINT *,'HH-norma-3',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 319 320 !print *,'b-end interpo',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2)) 321 322 !--Convert back to xred (reduced coordinates) 323 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 324 325 !--Interpolate alpha_in,alpha_end to obtain: alpha for the next step 326 call intparms(itime) 327 !print *,'b-end intparm',itime,sum(sum(alpha,dim=2)),sum(sum(alpha_in,dim=2)),sum(sum(alpha_end,dim=2)) 328 endif 329 330 !write(std_out,*) 'lotf 06' 331 !########################################################## 332 !### 06. Update the history with the prediction 333 334 xcart = xcart_next 335 xred = xred_next 336 337 !PRINT *,'HH-final-1',itime,icycle,sum(abs(xcart(1,:))),sum(abs(fcart(1,:))),sum(abs(vel(1,:))) 338 339 !print *,"XCART_final",sum(abs(xcart(1,:))),sum(abs(xcart(2,:))),sum(abs(xcart(3,:))) 340 341 write(message, '(a,3d20.9)' )& 342 & ' --- XCART_final',sum(abs(xcart(1,:))),sum(abs(xcart(2,:))),sum(abs(xcart(3,:))) 343 call wrtout(std_out,message,'PERS') 344 345 !Increase indexes 346 hist%ihist = abihist_findIndex(hist,+1) 347 348 !Fill the history with the variables 349 !xcart, xred, acell, rprimd 350 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 351 hist%vel(:,:,hist%ihist) = vel(:,:) 352 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 353 354 if(zDEBUG)then 355 write (std_out,*) 'xcart:' 356 do kk=1,ab_mover%natom 357 write (std_out,*) xcart(:,kk) 358 end do 359 write (std_out,*) 'fcart:' 360 do kk=1,ab_mover%natom 361 write (std_out,*) fcart(:,kk) 362 end do 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,*) 'strten:' 368 write (std_out,*) strten(1:3),ch10,strten(4:6) 369 write (std_out,*) 'etotal:' 370 write (std_out,*) etotal 371 end if 372 373 end subroutine pred_lotf