TABLE OF CONTENTS
ABINIT/pred_isothermal [ Functions ]
NAME
pred_isothermal
FUNCTION
Ionmov predictors (13) Isothermal integrator
IONMOV 13:
Reversible integrator of Martyna at al.
The equation of motion of the ions in contact with a thermostat
and a barostat are solved with the algorithm proposed by Martyna,
Tuckermann Tobias and Klein [Mol. Phys., 1996, p. 1117].
Related parameters : the time step (dtion),
the initial temperature mdtemp(1), the final temperature mdtemp(2),
the number of thermostats (nnos), and the masses of thermostats (qmass).
If optcell=1 or 2, the mass of the barostat (bmass)
must be given in addition.
There are three sub cases according to the value of optcell
optcell=0: isothermal
optcell=1: homogeneous cell fluctuations
optcell=2: full cell fluctuation in addition to temperature
control.
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
PARENTS
mover
CHILDREN
dsyev,hist2var,isopress,isostress,isotemp,metric,mkrdim,var2hist,wrtout xcart2xred,xred2xcart
SOURCE
56 #if defined HAVE_CONFIG_H 57 #include "config.h" 58 #endif 59 60 #include "abi_common.h" 61 62 subroutine pred_isothermal(ab_mover,hist,itime,mttk_vars,ntime,zDEBUG,iexit) 63 64 use defs_basis 65 use m_errors 66 use m_profiling_abi 67 use m_abimover 68 use m_abihist 69 use m_linalg_interfaces 70 71 use m_numeric_tools, only : uniformrandom 72 use m_geometry, only : mkrdim, xcart2xred, xred2xcart, metric 73 74 !This section has been created automatically by the script Abilint (TD). 75 !Do not modify the following lines by hand. 76 #undef ABI_FUNC 77 #define ABI_FUNC 'pred_isothermal' 78 use interfaces_14_hidewrite 79 use interfaces_45_geomoptim, except_this_one => pred_isothermal 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 type(mttk_type),intent(inout) :: mttk_vars 89 integer,intent(in) :: itime 90 integer,intent(in) :: ntime 91 integer,intent(in) :: iexit 92 logical,intent(in) :: zDEBUG 93 94 !Local variables------------------------------- 95 !scalars 96 integer :: ii,kk,iatom,idim,idum=5,ierr 97 integer,parameter :: lwork=8 98 real(dp) :: ucvol,ucvol0,ucvol_next,mttk_aloc,mttk_aloc2,mttk_bloc,ekin 99 real(dp) :: massvol=0 100 real(dp),parameter :: esh2=one/six,esh4=esh2/20._dp,esh6=esh4/42._dp 101 real(dp),parameter :: esh8=esh6/72._dp,nosetol=tol10,v2tol=tol8 102 real(dp) :: etotal,rescale_vel,polysh,s1,s2,sigma2,v2gauss,vtest 103 real(dp),save :: ktemp,vlogv 104 character(len=5000) :: message 105 !arrays 106 real(dp),allocatable,save :: fcart_m(:,:),vel_nexthalf(:,:) 107 108 real(dp) :: mttk_alc(3),mttk_alc2(3),mttk_blc(3),mttk_psh(3) 109 real(dp) :: mttk_tv(3,3),mttk_vt(3,3),mttk_ubox(3,3) 110 real(dp) :: mttk_uu(3),mttk_uv(3),mttk_veig(3) 111 real(dp) :: acell(3),acell0(3),acell_next(3) 112 real(dp) :: rprimd(3,3),rprimd0(3,3),rprim(3,3),rprimd_next(3,3),rprim_next(3,3) 113 real(dp) :: gprimd(3,3) 114 real(dp) :: gmet(3,3) 115 real(dp) :: rmet(3,3) 116 real(dp) :: fcart(3,ab_mover%natom) 117 !real(dp) :: fred_corrected(3,ab_mover%natom) 118 real(dp) :: xcart(3,ab_mover%natom),xcart_next(3,ab_mover%natom) 119 real(dp) :: xred(3,ab_mover%natom),xred_next(3,ab_mover%natom) 120 real(dp) :: vel(3,ab_mover%natom) 121 real(dp) :: strten(6),work(lwork) 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(vel_nexthalf)) then 132 ABI_DEALLOCATE(vel_nexthalf) 133 end if 134 return 135 end if 136 137 !write(std_out,*) 'isothermal 01' 138 !########################################################## 139 !### 01. Debugging and Verbose 140 141 if(zDEBUG)then 142 write(std_out,'(a,3a,41a,36a)') ch10,('-',kk=1,3),& 143 & 'Debugging and Verbose for pred_isothermal',('-',kk=1,36) 144 write(std_out,*) 'ionmov: ',13 145 write(std_out,*) 'itime: ',itime 146 end if 147 148 !write(std_out,*) 'isothermal 01' 149 !########################################################## 150 !### 01. Allocate the vectors vin, vout and hessian matrix 151 !### These arrays could be allocated from a previus 152 !### dataset that exit before itime==ntime 153 154 if(itime==1)then 155 if (allocated(fcart_m)) then 156 ABI_DEALLOCATE(fcart_m) 157 end if 158 if (allocated(vel_nexthalf)) then 159 ABI_DEALLOCATE(vel_nexthalf) 160 end if 161 end if 162 163 if (.not.allocated(fcart_m)) then 164 ABI_ALLOCATE(fcart_m,(3,ab_mover%natom)) 165 end if 166 if (.not.allocated(vel_nexthalf)) then 167 ABI_ALLOCATE(vel_nexthalf,(3,ab_mover%natom)) 168 end if 169 170 !write(std_out,*) 'isothermal 02' 171 !########################################################## 172 !### 02. Obtain the present values from the history 173 174 call hist2var(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 175 176 fcart(:,:)=hist%fcart(:,:,hist%ihist) 177 strten(:) =hist%strten(:,hist%ihist) 178 vel(:,:) =hist%vel(:,:,hist%ihist) 179 etotal =hist%etot(hist%ihist) 180 181 do ii=1,3 182 rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3) 183 end do 184 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 185 186 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 187 188 if(zDEBUG)then 189 write (std_out,*) 'fcart:' 190 do kk=1,ab_mover%natom 191 write (std_out,*) fcart(:,kk) 192 end do 193 write (std_out,*) 'vel:' 194 do kk=1,ab_mover%natom 195 write (std_out,*) vel(:,kk) 196 end do 197 write (std_out,*) 'strten:' 198 write (std_out,*) strten(1:3),ch10,strten(4:6) 199 write (std_out,*) 'etotal:' 200 write (std_out,*) etotal 201 end if 202 203 !Save initial values 204 acell0(:)=acell(:) 205 rprimd0(:,:)=rprimd(:,:) 206 ucvol0=ucvol 207 208 !Get rid of mean force on whole unit cell, but only if no 209 !generalized constraints are in effect 210 ! call fcart2fred(hist%fcart(:,:,hist%ihist),fred_corrected,rprimd,ab_mover%natom) 211 ! if(ab_mover%nconeq==0)then 212 ! amass_tot=sum(ab_mover%amass(:)) 213 ! do ii=1,3 214 ! if (ii/=3.or.ab_mover%jellslab==0) then 215 ! favg=sum(fred_corrected(ii,:))/dble(ab_mover%natom) 216 ! fred_corrected(ii,:)=fred_corrected(ii,:)-favg*ab_mover%amass(:)/amass_tot 217 ! end if 218 ! end do 219 ! end if 220 221 !write(std_out,*) 'isothermal 03' 222 !########################################################## 223 !### 05. Seconde half velocity step 224 225 if (itime>1) then 226 227 ! Next Half velocity step 228 do idim=1,3 229 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 230 end do 231 vel(:,:)=vel_nexthalf(:,:)+ab_mover%dtion/two*fcart_m(:,:) 232 233 if (ab_mover%optcell==0) then 234 ! Update Thermostat variables and velocity 235 call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,& 236 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel) 237 else if (ab_mover%optcell==1) then 238 ! Update Thermostat variables and velocity 239 call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 240 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,& 241 & strten,ab_mover%strtarget,ucvol,vel,vlogv) 242 else if (ab_mover%optcell==2) then 243 ! Next half step for extended variables 244 call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 245 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,& 246 & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel) 247 end if 248 249 if(itime==2) massvol=ekin+etotal 250 251 if (ab_mover%optcell==2) then 252 ! Evolution of cell and volume 253 acell_next(:)=acell(:) 254 ucvol_next=ucvol 255 end if 256 257 call mkrdim(acell,rprim,rprimd) 258 call metric(gmet,gprimd,-1,rmet,rprimd,ucvol) 259 260 if(zDEBUG)then 261 write(std_out,*) 'Second half velocity step' 262 write(std_out,*) 'Cell parameters:' 263 write(std_out,*) 'rprimd:' 264 do kk=1,3 265 write(std_out,*) rprimd(:,kk) 266 end do 267 write(std_out,*) 'rprim:' 268 do kk=1,3 269 write(std_out,*) rprim(:,kk) 270 end do 271 write(std_out,*) 'acell:' 272 write(std_out,*) acell(:) 273 write(std_out,*) 'Conserved energy:',(ekin+etotal)-massvol,ekin,etotal 274 write(std_out,*) 'Volume of unitqry cell (ucvol):',ucvol 275 end if 276 277 end if ! if (itime>1) 278 279 !write(std_out,*) 'isothermal 04' 280 !########################################################## 281 !### 03. Compute the next values 282 283 !The temperature is linear between initial and final values 284 !It is here converted from Kelvin to Hartree (kb_HaK) 285 ktemp=(ab_mover%mdtemp(1)+((ab_mover%mdtemp(2)-ab_mover%mdtemp(1))/dble(ntime-1))*(itime-1))*kb_HaK 286 287 if(zDEBUG)then 288 write(std_out,*) 'Temperature in Kelvin (ktemp):',ktemp 289 write(std_out,*) 'Initial temp (mdtemp(1)):',ab_mover%mdtemp(1) 290 write(std_out,*) 'Final temp (mdtemp(2)):',ab_mover%mdtemp(2) 291 write(std_out,*) 'Delay for atom permutation (delayperm)',ab_mover%delayperm 292 write(std_out,*) 'nnos:', ab_mover%nnos 293 write(std_out,*) 'qmass', ab_mover%qmass(:) 294 write(std_out,*) 'bmass',ab_mover%bmass 295 end if 296 297 if(itime==1) then 298 mttk_vars%glogs(:)=zero; mttk_vars%vlogs(:)=zero; mttk_vars%xlogs(:)=zero 299 mttk_vars%vboxg(:,:)=zero 300 vlogv=zero 301 ! v2gauss is twice the kinetic energy 302 v2gauss=0.0_dp 303 do iatom=1,ab_mover%natom 304 do idim=1,3 305 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 306 end do 307 end do 308 309 ! If there is no kinetic energy 310 if (v2gauss<=v2tol.and.itime==1) then 311 ! Maxwell-Boltzman distribution 312 v2gauss=zero 313 vtest=zero 314 do iatom=1,ab_mover%natom 315 do idim=1,3 316 vel(idim,iatom)=sqrt(kb_HaK*ab_mover%mdtemp(1)/ab_mover%amass(iatom))*cos(two_pi*uniformrandom(idum)) 317 vel(idim,iatom)=vel(idim,iatom)*sqrt(-2._dp*log(uniformrandom(idum))) 318 end do 319 end do 320 321 ! Get rid of center-of-mass velocity 322 s1=sum(ab_mover%amass(:)) 323 do idim=1,3 324 s2=sum(ab_mover%amass(:)*vel(idim,:)) 325 vel(idim,:)=vel(idim,:)-s2/s1 326 end do 327 328 ! Recompute v2gauss 329 do iatom=1,ab_mover%natom 330 do idim=1,3 331 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 332 vtest=vtest+vel(idim,iatom)/(3._dp*ab_mover%natom) 333 end do 334 end do 335 336 ! Now rescale the velocities to give the exact temperature 337 rescale_vel=sqrt(3._dp*ab_mover%natom*kb_HaK*ab_mover%mdtemp(1)/v2gauss) 338 vel(:,:)=vel(:,:)*rescale_vel 339 340 ! Recompute v2gauss with the rescaled velocities 341 v2gauss=zero 342 do iatom=1,ab_mover%natom 343 do idim=1,3 344 v2gauss=v2gauss+vel(idim,iatom)*vel(idim,iatom)*ab_mover%amass(iatom) 345 end do 346 end do 347 348 ! Compute the variance and print 349 sigma2=(v2gauss/(3._dp*ab_mover%natom)-ab_mover%amass(1)*vtest**2)/kb_HaK 350 351 if (zDEBUG)then 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 end if !(v2gauss<=v2tol.and.itime==1) 364 end if !(itime==1) 365 366 !XG070613 : Do not take away the following line , seems needed for the pathscale compiler 367 368 if (zDEBUG) write(std_out,*) 'vboxg',mttk_vars%vboxg(:,:) 369 370 371 !write(std_out,*) 'isothermal 05' 372 !########################################################## 373 !### 03. First half velocity step 374 375 !write(std_out,*) 'FIRST HALF VELOCITY STEP',ucvol 376 !write(std_out,*) 'OPTCELL option selected:',ab_mover%optcell 377 !write(std_out,*) 'RPRIMD' 378 !do kk=1,3 379 !write(std_out,*) rprimd(:,kk) 380 !end do 381 !write(std_out,*) 'RPRIM' 382 !do kk=1,3 383 !write(std_out,*) rprim(:,kk) 384 !end do 385 !write(std_out,*) 'ACELL' 386 !write(std_out,*) acell(:) 387 388 !%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 389 !%%% BEGIN sub case optcell=0 Isothermal Ensemble 390 !%%% 391 if(ab_mover%optcell==0) then 392 ! There is no evolution of cell 393 acell_next(:)=acell(:) 394 ucvol_next=ucvol 395 rprim_next(:,:)=rprim(:,:) 396 rprimd_next(:,:)=rprimd(:,:) 397 ! Update Thermostat variables and scale velocitie 398 call isotemp(ab_mover%amass,ab_mover%dtion,ekin,ab_mover%iatfix,& 399 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,vel) 400 401 ! Half velocity step 402 do idim=1,3 403 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 404 end do 405 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 406 ! New positions 407 ! Convert input xred (reduced coordinates) to xcart (cartesian) 408 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 409 xcart_next(:,:)=xcart(:,:)+vel_nexthalf(:,:)*ab_mover%dtion 410 ! Convert back to xred (reduced coordinates) 411 call xcart2xred(ab_mover%natom,rprimd,xcart_next,xred_next) 412 ! %%% 413 ! %%% END sub case optcell=0 Isothermal Ensemble 414 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 415 416 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 417 ! %%% BEGIN sub case optcell=1 Isothermal-Isenthalpic 418 ! %%% Ensemble (homogeneous cell deformation) 419 ! %%% 420 else if (ab_mover%optcell==1) then 421 ! Only homogeneous evolution of cell 422 ! Evolution of cell we keep rprim constant 423 rprim_next(:,:)=rprim(:,:) 424 ! Update Thermostat variables and velocity 425 call isopress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 426 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,ab_mover%qmass,& 427 & strten,ab_mover%strtarget,ucvol,vel,vlogv) 428 429 ! Half velocity step 430 do idim=1,3 431 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 432 end do 433 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 434 ! New positions 435 mttk_aloc=exp(ab_mover%dtion/two*vlogv) 436 mttk_aloc2=(vlogv*ab_mover%dtion/two)**2 437 polysh=(((esh8*mttk_aloc2+esh6)*mttk_aloc2+esh4)*mttk_aloc2+esh2)*mttk_aloc2+one 438 mttk_bloc=mttk_aloc*polysh*ab_mover%dtion 439 ! Convert input xred (reduced coordinates) to xcart (cartesian) 440 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 441 xcart_next(:,:)=xcart(:,:)*mttk_aloc**2+vel_nexthalf(:,:)*mttk_bloc 442 ! Update the volume and related quantities 443 acell_next(:)=acell(:)*exp(ab_mover%dtion*vlogv) 444 ! ucvol=ucvol*exp(ab_mover%dtion*vlogv) 445 call mkrdim(acell_next,rprim,rprimd_next) 446 call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol_next) 447 ! Convert back to xred (reduced coordinates) 448 call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next) 449 ! Computation of the forces for the new positions 450 ! Compute LDA forces (big loop) 451 452 ! COMMENTED 453 ! This should be in mover.F90 454 455 ! ! If metric has changed since the initialization, update the Ylm's 456 ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then 457 ! call status(0,dtfil%filstat,iexit,level,'call initylmg ') 458 ! option=0;if (ab_mover%iscf>0) option=1 459 ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,& 460 ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr) 461 ! end if 462 463 464 ! %%% 465 ! %%% END sub case optcell=1 Isothermal-Isenthalpic 466 ! %%% Ensemble (homogeneous cell deformation) 467 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 468 469 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 470 ! %%% BEGIN sub case optcell=2 Isothermal-Isenthalpic 471 ! %%% Ensemble (full cell deformation) 472 ! %%% 473 else if (ab_mover%optcell==2) then 474 acell_next=acell 475 ! Fisrt half step for extended variables 476 call isostress(ab_mover%amass,ab_mover%bmass,ab_mover%dtion,ekin,ab_mover%iatfix,& 477 & ktemp,mttk_vars,ab_mover%natom,ab_mover%nnos,& 478 & ab_mover%qmass,strten,ab_mover%strtarget,ucvol,vel) 479 ! Half velocity step 480 do idim=1,3 481 fcart_m(idim,:)=fcart(idim,:)/ab_mover%amass(:) 482 end do 483 vel_nexthalf(:,:)=vel(:,:)+ab_mover%dtion/two*fcart_m(:,:) 484 ! Convert input xred (reduced coordinates) to xcart (cartesian) 485 call xred2xcart(ab_mover%natom,rprimd,xcart,xred) 486 ! New positions 487 mttk_vt(:,:)=mttk_vars%vboxg(:,:) 488 call dsyev('V','U',3,mttk_vt,3,mttk_veig,work,lwork,ierr) 489 mttk_tv(:,:)=transpose(mttk_vt) 490 mttk_alc(:)=exp(ab_mover%dtion/two*mttk_veig(:)) 491 mttk_alc2(:)=(mttk_veig(:)*ab_mover%dtion/two)**2 492 mttk_psh(:)=(((esh8*mttk_alc2(:)+esh6)*mttk_alc2(:)+esh4)*mttk_alc2(:)+esh2)*mttk_alc2(:)+one 493 mttk_blc(:)=mttk_alc(:)*mttk_psh(:)*ab_mover%dtion 494 ! Update the positions 495 do iatom=1,ab_mover%natom 496 mttk_uu(:)=matmul(mttk_tv,xcart(:,iatom)) 497 mttk_uv(:)=matmul(mttk_tv,vel_nexthalf(:,iatom)) 498 mttk_uu(:)=mttk_uu(:)*mttk_alc(:)**2+mttk_uv(:)*mttk_blc(:) 499 xcart_next(:,iatom)=matmul(mttk_vt,mttk_uu) 500 end do 501 ! Update the box (rprimd and rprim) 502 mttk_ubox(:,:)=matmul(mttk_tv,rprimd) 503 do idim=1,3 504 mttk_ubox(:,idim)=mttk_ubox(:,idim)*mttk_alc(:)**2 505 end do 506 rprimd_next(:,:)=matmul(mttk_vt,mttk_ubox) 507 do idim=1,3 508 rprim_next(idim,:)=rprimd_next(idim,:)/acell(:) 509 end do 510 ! Update the volume 511 call metric(gmet,gprimd,-1,rmet,rprimd_next,ucvol) 512 ! Convert back to xred (reduced coordinates) 513 call xcart2xred(ab_mover%natom,rprimd_next,xcart_next,xred_next) 514 ! Computation of the forces for the new positions 515 516 ! COMMENTED 517 ! This should be in mover.F90 518 519 ! ! If metric has changed since the initialization, update the Ylm's 520 ! if (ab_mover%optcell/=0.and.psps%useylm==1.and.itime>1)then 521 ! call status(0,dtfil%filstat,iexit,level,'call initylmg ') 522 ! option=0;if (ab_mover%iscf>0) option=1 523 ! call initylmg(gprimd,kg,ab_mover%kptns,ab_mover%mkmem,mpi_enreg,psps%mpsang,ab_mover%mpw,ab_mover%nband,ab_mover%nkpt,& 524 ! & npwarr,ab_mover%nsppol,option,rprimd_next,ylm,ylmgr) 525 ! end if 526 527 ! %%% 528 ! %%% END sub case optcell=2 Isothermal-Isenthalpic 529 ! %%% Ensemble (full cell deformation) 530 ! %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 531 else 532 write(message, '(a,i12,a,a)' )& 533 & ' Disallowed value for optcell=',ab_mover%optcell,ch10,& 534 & ' Allowed values with ionmov==13 : 0 to 2.' 535 MSG_BUG(message) 536 end if 537 538 !write(std_out,*) 'OLD PARAMETERS' 539 !write(std_out,*) 'RPRIMD' 540 !do kk=1,3 541 !write(std_out,*) rprimd(:,kk) 542 !end do 543 !write(std_out,*) 'RPRIM' 544 !do kk=1,3 545 !write(std_out,*) rprim(:,kk) 546 !end do 547 !write(std_out,*) 'ACELL' 548 !write(std_out,*) acell(:) 549 550 !write(std_out,*) 'NEXT PARAMETERS' 551 !write(std_out,*) 'RPRIMD' 552 !do kk=1,3 553 !write(std_out,*) rprimd_next(:,kk) 554 !end do 555 !write(std_out,*) 'RPRIM' 556 !do kk=1,3 557 !write(std_out,*) rprim_next(:,kk) 558 !end do 559 !write(std_out,*) 'ACELL' 560 !write(std_out,*) acell_next(:) 561 562 563 !Those are the values store into the history 564 rprim=rprim_next 565 rprimd=rprimd_next 566 xred=xred_next 567 xcart=xcart_next 568 acell=acell_next 569 570 !write(std_out,*) 'isothermal 06' 571 !########################################################## 572 !### 06. Update the history with the prediction 573 574 !Increase indexes 575 hist%ihist = abihist_findIndex(hist,+1) 576 577 !Fill the history with the variables 578 !xred, acell, rprimd, vel 579 call var2hist(acell,hist,ab_mover%natom,rprimd,xred,zDEBUG) 580 hist%vel(:,:,hist%ihist)=vel(:,:) 581 hist%time(hist%ihist)=real(itime,kind=dp)*ab_mover%dtion 582 583 if(zDEBUG)then 584 write (std_out,*) 'fcart:' 585 do kk=1,ab_mover%natom 586 write (std_out,*) fcart(:,kk) 587 end do 588 write (std_out,*) 'vel:' 589 do kk=1,ab_mover%natom 590 write (std_out,*) vel(:,kk) 591 end do 592 write (std_out,*) 'strten:' 593 write (std_out,*) strten(1:3),ch10,strten(4:6) 594 write (std_out,*) 'etotal:' 595 write (std_out,*) etotal 596 end if 597 598 end subroutine pred_isothermal