TABLE OF CONTENTS
ABINIT/gstateimg [ Functions ]
NAME
gstateimg
FUNCTION
Routine for conducting DFT calculations for a set of (dynamical) images
COPYRIGHT
Copyright (C) 1998-2018 ABINIT group (XG, AR, GG, MT) 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
codvsn=code version
cpui=initial CPU time
nimage=number of images of the cell (treated by current proc)
=== Optional arguments (needed when nimage>1) ===
filnam(5)=character strings giving file names
filstat=character strings giving name of status file
idtset=index of the dataset
jdtset(0:ndtset)=actual index of the datasets
ndtset=number of datasets
OUTPUT
etotal_img=total energy, for each image fcart_img(3,natom,nimage)=forces, in cartesian coordinates, for each image fred_img(3,natom,nimage)=forces, in reduced coordinates, for each image npwtot(nkpt) = total number of plane waves at each k point strten_img(6,nimage)=stress tensor, for each image
SIDE EFFECTS
acell_img(3,nimage)=unit cell length scales (bohr), for each image amu_img(ntypat,nimage)=value of mass for each atomic type, for each image dtfil <type(datafiles_type)>=variables related to files dtset <type(dataset_type)>=all input variables in this dataset | mband =maximum number of bands (IN) | mgfft =maximum single fft dimension (IN) | mkmem =maximum number of k points which can fit in core memory (IN) | mpw =maximum number of planewaves in basis sphere (large number) (IN) | natom =number of atoms in unit cell (IN) | nfft =(effective) number of FFT grid points (for this processor) (IN) | nkpt =number of k points (IN) | nspden=number of spin-density components (IN) | nsppol=number of channels for spin-polarization (1 or 2) (IN) | nsym =number of symmetry elements in space group iexit= exit flag mixalch_img(npspalch,ntypalch,nimage)=value of alchemical mixing factors,for each image mpi_enreg=MPI-parallelisation information (some already initialized, some others to be initialized here) occ_img(mband*nkpt*nsppol,nimage) = occupation number for each band and k, for each image pawang <type(pawang_type)>=paw angular mesh and related data pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data psps <type(pseudopotential_type)>=variables related to pseudopotentials Before entering the first time in gstateimg, a significant part of psps has been initialized : the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid, ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp All the remaining components of psps are to be initialized in the call to pspini . The next time the code enters gstateimg, psps might be identical to the one of the previous dtset, in which case, no reinitialisation is scheduled in pspini.f . rprim_img(3,3,nimage)=dimensionless real space primitive translations, for each image vel_cell_img(3,3,nimage)=value of cell parameters velocities, for each image vel_img(3,natom,nimage)=value of atomic velocities,for each image xred_img(3,natom,nimage) = reduced atomic coordinates, for each image
NOTES
USE OF FFT GRIDS:
=================
In case of PAW:
---------------
Two FFT grids are used:
- A "coarse" FFT grid (defined by ecut)
for the application of the Hamiltonian on the plane waves basis.
It is defined by nfft, ngfft, mgfft, ...
Hamiltonian, wave-functions, density related to WFs (rhor here), ...
are expressed on this grid.
- A "fine" FFT grid (defined) by ecutdg)
for the computation of the density inside PAW spheres.
It is defined by nfftf, ngfftf, mgfftf, ...
Total density, potentials, ...
are expressed on this grid.
In case of norm-conserving:
---------------------------
- Only the usual FFT grid (defined by ecut) is used.
It is defined by nfft, ngfft, mgfft, ...
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
In case of wavelets:
--------------------
- Only the usual FFT grid (defined by wvl_crmult) is used.
It is defined by nfft, ngfft, mgfft, ... This is strictly not
an FFT grid since its dimensions are not suited for FFTs. They are
defined by wvl_setngfft().
For compatibility reasons, (nfftf,ngfftf,mgfftf)
are set equal to (nfft,ngfft,mgfft) in that case.
TODO
Not yet possible to use restartxf in parallel when localrdwf==0
PARENTS
driver,gwls_sternheimer
CHILDREN
abihist_bcast,abihist_copy,abihist_free,abihist_init,args_gs_free args_gs_init,copy_results_img,destroy_results_img,dtfil_init,fcart2fred ga_destroy,ga_init,gstate,init_results_img,libpaw_spmsg_mpisum localfilnam,localrdfile,localredirect,localwrfile,mep_destroy,mep_init mkradim,mkrdim,pimd_destroy,pimd_init,pimd_skip_qtb,predictimg,prtimg read_md_hist_img,scf_history_free,scf_history_nullify,specialmsg_mpisum status,timab,var2hist,vel2hist,write_md_hist_img,wrtout,xmpi_barrier xmpi_sum
SOURCE
122 #if defined HAVE_CONFIG_H 123 #include "config.h" 124 #endif 125 126 #include "abi_common.h" 127 128 subroutine gstateimg(acell_img,amu_img,codvsn,cpui,dtfil,dtset,etotal_img,fcart_img,& 129 & fred_img,iexit,mixalch_img,mpi_enreg,nimage,npwtot,occ_img,& 130 & pawang,pawrad,pawtab,psps,& 131 & rprim_img,strten_img,vel_cell_img,vel_img,wvl,xred_img,& 132 & filnam,filstat,idtset,jdtset,ndtset) ! optional arguments 133 134 135 use defs_basis 136 use defs_datatypes 137 use defs_abitypes 138 use defs_wvltypes 139 use defs_parameters 140 use defs_rectypes 141 use m_profiling_abi 142 use m_abihist 143 use m_mep 144 use m_ga 145 use m_pimd 146 use m_xmpi 147 use m_errors 148 use m_rec 149 use m_args_gs 150 use m_results_img 151 use m_scf_history 152 use m_io_redirect 153 154 use m_geometry, only : mkradim, mkrdim, fcart2fred 155 use m_specialmsg, only : specialmsg_mpisum 156 use m_libpaw_tools, only : libpaw_spmsg_mpisum 157 use m_pawang, only : pawang_type 158 use m_pawrad, only : pawrad_type 159 use m_pawtab, only : pawtab_type 160 161 #if defined HAVE_BIGDFT 162 use BigDFT_API, only: mpi_environment_set 163 #endif 164 165 !This section has been created automatically by the script Abilint (TD). 166 !Do not modify the following lines by hand. 167 #undef ABI_FUNC 168 #define ABI_FUNC 'gstateimg' 169 use interfaces_14_hidewrite 170 use interfaces_18_timing 171 use interfaces_32_util 172 use interfaces_45_geomoptim 173 use interfaces_67_common 174 use interfaces_95_drive, except_this_one => gstateimg 175 !End of the abilint section 176 177 implicit none 178 179 !Arguments ------------------------------------ 180 !scalars 181 integer,intent(in) :: nimage 182 integer,optional,intent(in) :: idtset,ndtset 183 integer,intent(inout) :: iexit 184 real(dp),intent(in) :: cpui 185 character(len=6),intent(in) :: codvsn 186 character(len=fnlen),optional,intent(in) :: filstat 187 type(MPI_type),intent(inout) :: mpi_enreg 188 type(datafiles_type),target,intent(inout) :: dtfil 189 type(dataset_type),target,intent(inout) :: dtset 190 type(pawang_type),intent(inout) :: pawang 191 type(pseudopotential_type),intent(inout) :: psps 192 type(wvl_data),intent(inout) :: wvl 193 !arrays 194 integer,optional,intent(in) :: jdtset(:) 195 integer,intent(out) :: npwtot(dtset%nkpt) 196 character(len=fnlen),optional,intent(in) :: filnam(:) 197 real(dp), intent(out) :: etotal_img(nimage),fcart_img(3,dtset%natom,nimage) 198 real(dp), intent(out) :: fred_img(3,dtset%natom,nimage),strten_img(6,nimage) 199 real(dp),intent(inout) :: acell_img(3,nimage),amu_img(dtset%ntypat,nimage) 200 real(dp),intent(inout) :: mixalch_img(dtset%npspalch,dtset%ntypalch,nimage) 201 real(dp),intent(inout) :: occ_img(dtset%mband*dtset%nkpt*dtset%nsppol,nimage) 202 real(dp),intent(inout) :: rprim_img(3,3,nimage),vel_cell_img(3,3,nimage),vel_img(3,dtset%natom,nimage) 203 real(dp),intent(inout) :: xred_img(3,dtset%natom,nimage) 204 type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw) 205 type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw) 206 207 !Local variables------------------------------- 208 !Define file format for different type of files. Presently, 209 !only one file format is supported for each type of files, but this might 210 !change soon ... 211 !2 for wavefunction file, new format (version 2.0 and after) (fform) NOT USED 212 !52 for density rho(r) (fformr) 213 !102 for potential V(r) file. (fformv) NOT USED 214 !scalars 215 integer,parameter :: formeig=0,level=100,ndtpawuj=0,response=0 216 integer :: history_size,idelta,idynimage,ierr,ifirst 217 integer :: ii,iimage,ih,itimimage,itimimage_eff,itimimage_prev,ndynimage,nocc 218 integer :: ntimimage,ntimimage_stored,ntimimage_max 219 logical :: check_conv,compute_all_images,compute_static_images 220 logical :: isVused,isARused,is_master,is_mep,is_pimd 221 logical :: call_predictor,use_hist,use_hist_prev 222 real(dp) :: delta_energy,dtion 223 character(len=500) :: hist_filename,msg 224 type(args_gs_type) :: args_gs 225 type(mep_type) :: mep_param 226 type(ga_type) :: ga_param 227 type(pimd_type) :: pimd_param 228 !arrays 229 integer,allocatable :: list_dynimage(:),scf_initialized(:) 230 character(len=60),parameter :: imagealgo_str(0:13)=(/ & 231 & 'IMAGE COPY ',& ! 0 232 & 'IMAGE STEEPEST DESCENT ',& ! 1 233 & 'STRING METHOD ',& ! 2 234 & 'METADYNAMICS ',& ! 3 235 & 'GENETIC ALGORITHM ',& ! 4 236 & 'NUDGED ELASTIC BAND ',& ! 5 237 & ' ',& ! 6 238 & ' ',& ! 7 239 & ' ',& ! 8 240 & 'PATH-INTEGRAL MOLECULAR DYNAMICS (LANGEVIN) ',& ! 9 241 & 'PATH-INTEGRAL MOLECULAR DYNAMICS (QUANTUM THERMAL BATH) ',& ! 10 242 & ' ',& ! 11 243 & ' ',& ! 12 244 & 'PATH-INTEGRAL MOLECULAR DYNAMICS (CHAIN OF THERMOSTATS) '/) ! 13 245 real(dp) :: acell(3),rprim(3,3),rprimd(3,3),tsec(2) 246 real(dp),allocatable :: amass(:,:),occ(:),vel(:,:),vel_cell(:,:),xred(:,:) 247 type(abihist),allocatable :: hist(:),hist_prev(:) 248 type(results_img_type),pointer :: results_img_timimage(:,:),res_img(:) 249 type(scf_history_type),allocatable :: scf_history(:) 250 251 ! *********************************************************************** 252 253 DBG_ENTER("COLL") 254 255 call timab(700,1,tsec) 256 call timab(703,3,tsec) 257 258 call status(0,dtfil%filstat,iexit,level,'enter ') 259 260 !Arguments check 261 if (dtset%nimage>1) then 262 if ((.not.present(filnam)).or.(.not.present(filnam)).or.(.not.present(idtset)).or.& 263 & (.not.present(ndtset)).or.(.not.present(jdtset))) then 264 write(msg,'(3a)') & 265 & 'When nimage>1, all the following argument should be present:',ch10,& 266 & 'filnam, filstat, idtset, ndtset, jdtset !' 267 MSG_BUG(msg) 268 end if 269 end if 270 271 !Set flag for the effective computation (once) of static images 272 !For the time being only set when parallelization is activated 273 !Note: if you modify this flag, do not forget to change it in outvars and outvar1 274 compute_static_images=(dtset%istatimg>0) 275 276 !Prepare the calculation, by computing flags and dimensions 277 is_pimd=(dtset%imgmov==9.or.dtset%imgmov==10.or.dtset%imgmov==13) 278 is_mep =(dtset%imgmov==1.or.dtset%imgmov== 2.or.dtset%imgmov== 5) 279 ntimimage=dtset%ntimimage 280 ntimimage_stored=ntimimage;if(is_pimd)ntimimage_stored=2 281 nocc=dtset%mband*dtset%nkpt*dtset%nsppol 282 is_master=(mpi_enreg%me_cell==0.and.mpi_enreg%me_img==0) 283 delta_energy=zero 284 285 !Management of dynamics/relaxation history (positions, forces, stresses, ...) 286 use_hist=(dtset%imgmov/=0.and.nimage>0) ; use_hist_prev=.false. 287 isVused=is_pimd;isARused=(dtset%optcell/=0) 288 if (use_hist) then 289 !Read history from file (and broadcast if MPI) 290 #if defined HAVE_NETCDF 291 use_hist_prev=(dtset%restartxf==-1.and.nimage>0) 292 #endif 293 hist_filename=trim(dtfil%filnam_ds(4))//'_HIST.nc' 294 if (use_hist_prev)then 295 ABI_DATATYPE_ALLOCATE(hist_prev,(nimage)) 296 if (mpi_enreg%me_cell==0) then 297 call read_md_hist_img(hist_filename,hist_prev,isVused,isARused,& 298 & imgtab=mpi_enreg%my_imgtab) 299 end if 300 call abihist_bcast(hist_prev,0,mpi_enreg%comm_cell) 301 if (nimage>0) then 302 if (any(hist_prev(:)%mxhist/=hist_prev(1)%mxhist)) then 303 msg='History problem: all images should have the same number of time steps!' 304 MSG_ERROR(msg) 305 end if 306 use_hist_prev=(hist_prev(1)%mxhist>0) 307 if (use_hist_prev) ntimimage=ntimimage+hist_prev(1)%mxhist 308 end if 309 if (.not.use_hist_prev) then 310 call abihist_free(hist_prev) 311 ABI_DATATYPE_DEALLOCATE(hist_prev) 312 end if 313 end if 314 !Initialize a variable to write the history 315 ABI_DATATYPE_ALLOCATE(hist,(nimage)) 316 call abihist_init(hist,dtset%natom,ntimimage,isVused,isARused) 317 end if ! imgmov/=0 318 319 !Allocations 320 ABI_ALLOCATE(occ,(nocc)) 321 ABI_ALLOCATE(vel,(3,dtset%natom)) 322 ABI_ALLOCATE(vel_cell,(3,3)) 323 ABI_ALLOCATE(xred,(3,dtset%natom)) 324 ABI_DATATYPE_ALLOCATE(results_img_timimage,(nimage,ntimimage_stored)) 325 ABI_ALLOCATE(list_dynimage,(dtset%ndynimage)) 326 do itimimage=1,ntimimage_stored 327 res_img => results_img_timimage(:,itimimage) 328 call init_results_img(dtset%natom,dtset%npspalch,dtset%nsppol,dtset%ntypalch,& 329 & dtset%ntypat,res_img) 330 do iimage=1,nimage 331 res_img(iimage)%acell(:) =acell_img(:,iimage) 332 res_img(iimage)%amu(:) =amu_img(:,iimage) 333 res_img(iimage)%mixalch(:,:) =mixalch_img(:,:,iimage) 334 res_img(iimage)%rprim(:,:) =rprim_img(:,:,iimage) 335 res_img(iimage)%xred(:,:) =xred_img(:,:,iimage) 336 res_img(iimage)%vel(:,:) =vel_img(:,:,iimage) 337 res_img(iimage)%vel_cell(:,:)=vel_cell_img(:,:,iimage) 338 end do 339 end do 340 ndynimage=0 341 do iimage=1,nimage 342 ii=mpi_enreg%my_imgtab(iimage) 343 if (dtset%dynimage(ii)==1) then 344 ndynimage=ndynimage+1 345 list_dynimage(ndynimage)=iimage 346 end if 347 end do 348 349 !Management of SCF history (density/WF predictions from one time step to another) 350 ABI_DATATYPE_ALLOCATE(scf_history,(nimage)) 351 ABI_ALLOCATE(scf_initialized,(nimage)) 352 scf_initialized=0 353 history_size=-1 354 if (dtset%ntimimage<=1) then 355 if (dtset%usewvl==0.and.dtset%ionmov>0.and. & 356 & (abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) history_size=2 357 else 358 if (abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==3) history_size=0 359 if (dtset%usewvl==0.and.(abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) history_size=2 360 end if 361 do iimage=1,nimage 362 call scf_history_nullify(scf_history(iimage)) 363 scf_history(iimage)%history_size=history_size 364 end do 365 366 !MEP search: fill in eventually the data structure mep_param 367 call mep_init(dtset,mep_param) 368 369 !GA search: fill in eventually the data structure ga_param 370 call ga_init(dtset,ga_param) 371 372 !PIMD: fill in eventually the data structure pimd_param 373 call pimd_init(dtset,pimd_param,is_master) 374 dtion=one;if (is_pimd) dtion=pimd_param%dtion 375 376 !In some cases, need amass variable 377 if (use_hist) then 378 ABI_ALLOCATE(amass,(dtset%natom,nimage)) 379 do iimage=1,nimage 380 if (any(amu_img(:,iimage)/=amu_img(:,1))) then 381 msg='HIST file is not compatible with variable masses!' 382 MSG_ERROR(msg) 383 end if 384 amass(:,iimage)=amu_emass*amu_img(dtset%typat(:),iimage) 385 end do 386 end if 387 388 !In the case of the 4th-order Runge-Kutta solver, 389 !one must have a number of step multiple of 4. 390 ntimimage_max=ntimimage;idelta=1 391 if (dtset%imgmov==2.and.mep_param%mep_solver==4) then 392 ntimimage_max=4*(ntimimage_max/4) 393 idelta=4 394 end if 395 396 call timab(703,2,tsec) 397 398 !----------------------------------------------------------------------------------------- 399 !Big loop on the propagation of all images 400 itimimage_eff=1 401 do itimimage=1,ntimimage 402 403 res_img => results_img_timimage(:,itimimage_eff) 404 call_predictor=(ntimimage>1) 405 406 ! If history is activated and if current image is inside it: do not compute anything 407 if (use_hist_prev) then 408 if (all(hist_prev(:)%ihist<=hist_prev(:)%mxhist)) then 409 do iimage=1,nimage 410 ih=hist_prev(iimage)%ihist 411 call abihist_copy(hist_prev(iimage),hist(iimage)) 412 call mkradim(hist_prev(iimage)%acell(:,ih),rprim,hist_prev(iimage)%rprimd(:,:,ih)) 413 res_img(iimage)%acell(:)=hist_prev(iimage)%acell(:,ih) 414 res_img(iimage)%rprim(:,:)=rprim 415 res_img(iimage)%xred(:,:)=hist_prev(iimage)%xred(:,:,ih) 416 res_img(iimage)%vel(:,:)=hist_prev(iimage)%vel(:,:,ih) 417 res_img(iimage)%vel_cell(:,:)=hist_prev(iimage)%vel_cell(:,:,ih) 418 res_img(iimage)%results_gs%fcart(:,:)=hist_prev(iimage)%fcart(:,:,ih) 419 res_img(iimage)%results_gs%strten(:)=hist_prev(iimage)%strten(:,ih) 420 res_img(iimage)%results_gs%etotal=hist_prev(iimage)%etot(ih) 421 res_img(iimage)%results_gs%energies%entropy=hist_prev(iimage)%entropy(ih) 422 call fcart2fred(res_img(iimage)%results_gs%fcart,res_img(iimage)%results_gs%fred,& 423 & hist_prev(iimage)%rprimd(:,:,ih),dtset%natom) 424 hist_prev(iimage)%ihist=hist_prev(iimage)%ihist+1 425 end do 426 !PI-QTB: skip a record in random force file 427 if (pimd_param%use_qtb==1) call pimd_skip_qtb(pimd_param) 428 !call_predictor=.false. 429 goto 110 ! This is temporary 430 end if 431 end if 432 433 call timab(704,1,tsec) 434 call localfilnam(mpi_enreg%comm_img,mpi_enreg%comm_cell,mpi_enreg%comm_world,filnam,'_IMG',dtset%nimage) 435 compute_all_images=(compute_static_images.and.itimimage==1) 436 437 ! Print title for time step 438 if (dtset%nimage>1.or.dtset%ntimimage>1) then 439 if (dtset%prtvolimg<2) then 440 msg=ch10;if (itimimage >1) write(msg,'(2a)') ch10,ch10 441 write(msg,'(5a)') trim(msg),& 442 & '================================================================================',& 443 & ch10,' ',trim(imagealgo_str(dtset%imgmov)) 444 else 445 msg='';if (itimimage >1) msg=ch10 446 write(msg,'(5a)') trim(msg),& 447 & '--------------------------------------------------------------------------------',& 448 & ch10,' ',trim(imagealgo_str(dtset%imgmov)) 449 end if 450 if (dtset%ntimimage==1) write(msg,'(2a)') trim(msg),' FOR 1 TIME STEP' 451 if (dtset%ntimimage >1) write(msg,'(2a,i5)') trim(msg),' - TIME STEP ',itimimage 452 if (dtset%prtvolimg<2) then 453 write(msg,'(3a)') trim(msg),ch10,& 454 & '================================================================================' 455 end if 456 call wrtout(ab_out ,msg,'COLL') 457 call wrtout(std_out,msg,'PERS') 458 end if 459 460 call timab(704,2,tsec) 461 462 ! Loop on the dynamical images 463 idynimage=0 464 do iimage=1,nimage 465 466 call timab(705,1,tsec) 467 468 ii=mpi_enreg%my_imgtab(iimage) 469 if (dtset%dynimage(ii)==1) idynimage=idynimage+1 470 471 ! Compute static image only at first time step 472 if (dtset%dynimage(ii)==1.or.compute_all_images) then 473 474 ! Change file names according to image index (if nimage>1) 475 if (dtset%nimage>1) then 476 call dtfil_init(dtfil,dtset,filnam,filstat,idtset,jdtset,mpi_enreg,ndtset,& 477 & image_index=ii) 478 if (itimimage>1) then 479 dtfil%ireadwf=0;dtfil%ireadden=0;dtfil%ireadkden=0 480 end if 481 end if 482 483 ! Redefine output units 484 call localwrfile(mpi_enreg%comm_cell,ii,dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg) 485 486 ! Print title for image 487 if (dtset%nimage>1.and.(dtset%prtvolimg==0.or.do_write_log)) then 488 if (ii==1) write(msg,'(a)' ) ch10 489 if (ii >1) write(msg,'(2a)') ch10,ch10 490 write(msg,'(6a,i4,a,i4,3a)') trim(msg),& 491 & '--------------------------------------------------------------------------------',ch10,& 492 & ' ',trim(imagealgo_str(dtset%imgmov)),' - CELL # ',ii,'/',dtset%nimage,ch10,& 493 & '--------------------------------------------------------------------------------',ch10 494 if (dtset%prtvolimg==0) then 495 call wrtout(ab_out ,msg,'COLL') 496 end if 497 if (do_write_log) then 498 call wrtout(std_out,msg,'PERS') 499 end if 500 end if 501 502 acell(:) =res_img(iimage)%acell(:) 503 rprim(:,:) =res_img(iimage)%rprim(:,:) 504 vel(:,:) =res_img(iimage)%vel(:,:) 505 vel_cell(:,:)=res_img(iimage)%vel_cell(:,:) 506 xred(:,:) =res_img(iimage)%xred(:,:) 507 occ(:) =occ_img(:,iimage) 508 509 call args_gs_init(args_gs, & 510 & res_img(iimage)%amu(:),res_img(iimage)%mixalch(:,:),& 511 & dtset%dmatpawu(:,:,:,:,ii),dtset%upawu(:,ii),dtset%jpawu(:,ii),& 512 & dtset%rprimd_orig(:,:,ii)) 513 514 call timab(705,2,tsec) 515 516 call status(idynimage+100*itimimage,dtfil%filstat,iexit,level,'call gstate ') 517 call gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,scf_initialized(iimage),& 518 & mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,& 519 & res_img(iimage)%results_gs,& 520 & rprim,scf_history(iimage),vel,vel_cell,wvl,xred) 521 522 call timab(706,1,tsec) 523 524 call args_gs_free(args_gs) 525 526 if (dtset%dynimage(ii)==1) then 527 res_img(iimage)%acell(:) =acell(:) 528 res_img(iimage)%rprim(:,:) =rprim(:,:) 529 res_img(iimage)%vel(:,:) =vel(:,:) 530 res_img(iimage)%vel_cell(:,:)=vel_cell(:,:) 531 res_img(iimage)%xred(:,:) =xred(:,:) 532 occ_img(:,iimage) =occ(:) 533 end if 534 535 ! check change de rprim et reecriture dans hist 536 ! check change de xred et reecriture dans hist 537 538 ! Close output units ; restore defaults 539 call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg) 540 call timab(706,2,tsec) 541 542 else if (itimimage>1) then ! For static images, simply copy one time step to the other 543 itimimage_prev=itimimage_eff-1 544 if (itimimage_prev<1) itimimage_prev=ntimimage_stored 545 call copy_results_img(results_img_timimage(iimage,itimimage_prev), & 546 & results_img_timimage(iimage,itimimage_eff )) 547 end if 548 549 ! Store results in hist datastructure 550 if (use_hist) then 551 ih=hist(iimage)%ihist 552 call mkrdim(res_img(iimage)%acell(:),res_img(iimage)%rprim(:,:),rprimd) 553 call var2hist(res_img(iimage)%acell(:),hist(iimage),dtset%natom,& 554 & rprimd,res_img(iimage)%xred(:,:),.FALSE.) 555 call vel2hist(amass(:,iimage),hist(iimage),res_img(iimage)%vel(:,:),& 556 & res_img(iimage)%vel_cell(:,:)) 557 hist(iimage)%fcart(:,:,ih)=res_img(iimage)%results_gs%fcart(:,:) 558 hist(iimage)%strten(:,ih)=res_img(iimage)%results_gs%strten(:) 559 hist(iimage)%etot(ih)=res_img(iimage)%results_gs%etotal 560 hist(iimage)%entropy(ih)=res_img(iimage)%results_gs%energies%entropy 561 hist(iimage)%time(ih)=real(itimimage,kind=dp)*dtion 562 end if 563 564 end do ! iimage 565 566 if(mpi_enreg%paral_img==1)then 567 call timab(702,1,tsec) 568 call xmpi_barrier(mpi_enreg%comm_img) 569 call timab(702,2,tsec) 570 end if 571 572 call timab(707,1,tsec) 573 574 ! Output when images are used 575 if (dtset%nimage>1) then 576 ! === 1st option: reduced outputs === 577 if (dtset%prtvolimg>0) then 578 call prtimg(dtset%dynimage,imagealgo_str(dtset%imgmov),dtset%imgmov,ab_out,& 579 & mpi_enreg,nimage,dtset%nimage,compute_all_images,dtset%prtvolimg,res_img) 580 end if 581 end if 582 583 ! Manage log files when images are used 584 call localrdfile(mpi_enreg%comm_img,mpi_enreg%comm_world,compute_all_images,& 585 & dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg,dyn=dtset%dynimage) 586 587 ! Write hist datastructure in HIST file 588 #if defined HAVE_NETCDF 589 if (use_hist.and.mpi_enreg%me_cell==0) then 590 ifirst=merge(0,1,itimimage>1) 591 call write_md_hist_img(hist,hist_filename,ifirst,itimimage,dtset%natom,dtset%ntypat,& 592 & dtset%typat,amu_img(:,1),dtset%znucl,dtion,& 593 & nimage=dtset%nimage,imgmov=dtset%imgmov,mdtemp=dtset%mdtemp,comm_img=mpi_enreg%comm_img,& 594 & imgtab=mpi_enreg%my_imgtab) 595 end if 596 #endif 597 598 ! TESTS WHETHER ONE CONTINUES THE LOOP 599 ! Here we calculate the change in energy, and exit if delta_energy < tolimg 600 delta_energy=zero 601 ! Doesn't check convergence in case of PIMD 602 check_conv=((.not.is_pimd).and.itimimage>1) 603 ! In case of 4th-order Runge-Kutta, does check convergence every 4 steps 604 if (dtset%imgmov==2.and.mep_param%mep_solver==4) then 605 check_conv=(mod(itimimage,4)==0.and.itimimage>4) 606 end if 607 if (check_conv) then 608 do idynimage=1,ndynimage 609 iimage=list_dynimage(idynimage) 610 delta_energy=delta_energy & 611 & +abs(results_img_timimage(iimage,itimimage)%results_gs%etotal & 612 & -results_img_timimage(iimage,itimimage-idelta)%results_gs%etotal) 613 end do 614 if (mpi_enreg%paral_img==1) then 615 call xmpi_sum(delta_energy,mpi_enreg%comm_img,ierr) 616 end if 617 delta_energy=delta_energy/dtset%ndynimage 618 if (delta_energy<dtset%tolimg) then 619 if (dtset%prtvolimg<2) then 620 write(msg,'(5a,i5,6a,es11.3,a,es11.3,2a)') ch10,ch10,& 621 & '================================================================================',ch10,& 622 & ' At time step ',itimimage,ch10,& 623 & ' ',trim(imagealgo_str(dtset%imgmov)),' has reached energy convergence',ch10,& 624 & ' with Average[Abs(Etotal(t)-Etotal(t-dt))]=',delta_energy,'<tolimg=',dtset%tolimg,ch10,& 625 & '================================================================================' 626 else 627 write(msg,'(4a,i5,6a,es11.3,a,es11.3)') ch10,& 628 & '--------------------------------------------------------------------------------',ch10,& 629 & ' At time step ',itimimage,ch10,& 630 & ' ',trim(imagealgo_str(dtset%imgmov)),' has reached energy convergence',ch10,& 631 & ' with Average[Abs(Etotal(t)-Etotal(t-dt))]=',delta_energy,'<tolimg=',dtset%tolimg 632 end if 633 call wrtout(ab_out ,msg,'COLL') 634 call wrtout(std_out,msg,'COLL') 635 call timab(707,2,tsec) 636 exit ! exit itimimage 637 end if 638 end if 639 640 !Temporary statement 641 110 continue 642 643 ! Dont call the predictor at last time step 644 if (itimimage>=ntimimage_max) call_predictor=(call_predictor.and.is_pimd) 645 646 ! Predict the next value of the images 647 if (call_predictor) then 648 call predictimg(delta_energy,imagealgo_str(dtset%imgmov),dtset%imgmov,itimimage,& 649 & itimimage_eff,list_dynimage,ga_param,mep_param,mpi_enreg,dtset%natom,ndynimage,& 650 & nimage,dtset%nimage,ntimimage_stored,pimd_param,dtset%prtvolimg,results_img_timimage) 651 end if 652 653 ! Increment indexes 654 if (itimimage>=ntimimage_max) exit 655 itimimage_eff=itimimage_eff+1;if (itimimage_eff>ntimimage_stored) itimimage_eff=1 656 if (use_hist) then 657 do iimage=1,nimage 658 hist(iimage)%ihist=hist(iimage)%ihist+1 659 end do 660 end if 661 662 call timab(707,2,tsec) 663 664 end do ! itimimage 665 !----------------------------------------------------------------------------------------- 666 667 call timab(708,1,tsec) 668 669 !Copy the results of the computation in the appropriate arguments of the routine 670 do iimage=1,nimage 671 ii=mpi_enreg%my_imgtab(iimage) 672 if (dtset%dynimage(ii)==1) then 673 acell_img(:,iimage) =results_img_timimage(iimage,itimimage_eff)%acell(:) 674 amu_img(:,iimage) =results_img_timimage(iimage,itimimage_eff)%amu(:) 675 mixalch_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%mixalch(:,:) 676 rprim_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%rprim(:,:) 677 vel_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%vel(:,:) 678 vel_cell_img(:,:,iimage)=results_img_timimage(iimage,itimimage_eff)%vel_cell(:,:) 679 xred_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%xred(:,:) 680 etotal_img(iimage) =results_img_timimage(iimage,itimimage_eff)%results_gs%etotal 681 fcart_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%results_gs%fcart(:,:) 682 fred_img(:,:,iimage) =results_img_timimage(iimage,itimimage_eff)%results_gs%fred(:,:) 683 strten_img(:,iimage) =results_img_timimage(iimage,itimimage_eff)%results_gs%strten(:) 684 else if (compute_static_images) then 685 etotal_img(iimage) =results_img_timimage(iimage,1)%results_gs%etotal 686 fcart_img(:,:,iimage) =results_img_timimage(iimage,1)%results_gs%fcart(:,:) 687 fred_img(:,:,iimage) =results_img_timimage(iimage,1)%results_gs%fred(:,:) 688 strten_img(:,iimage) =results_img_timimage(iimage,1)%results_gs%strten(:) 689 end if 690 end do 691 692 !When parallelizattion over images is activated, has to sum number of warnings 693 !and comments written in log file 694 if (mpi_enreg%paral_img==1.and.mpi_enreg%comm_cell==0) then 695 call specialmsg_mpisum(mpi_enreg%comm_img) 696 call libpaw_spmsg_mpisum(mpi_enreg%comm_img) 697 end if 698 !Final deallocations 699 ABI_DEALLOCATE(occ) 700 ABI_DEALLOCATE(vel) 701 ABI_DEALLOCATE(vel_cell) 702 ABI_DEALLOCATE(xred) 703 ABI_DEALLOCATE(list_dynimage) 704 if (allocated(amass)) then 705 ABI_DEALLOCATE(amass) 706 end if 707 708 do itimimage=1,ntimimage_stored 709 call destroy_results_img(results_img_timimage(:,itimimage)) 710 end do 711 ABI_DATATYPE_DEALLOCATE(results_img_timimage) 712 do iimage=1,nimage 713 call scf_history_free(scf_history(iimage)) 714 end do 715 ABI_DATATYPE_DEALLOCATE(scf_history) 716 ABI_DEALLOCATE(scf_initialized) 717 if (allocated(hist_prev)) then 718 call abihist_free(hist_prev) 719 ABI_DATATYPE_DEALLOCATE(hist_prev) 720 end if 721 if (allocated(hist)) then 722 call abihist_free(hist) 723 ABI_DATATYPE_DEALLOCATE(hist) 724 end if 725 726 call mep_destroy(mep_param) 727 call ga_destroy(ga_param) 728 call pimd_destroy(pimd_param) 729 730 call status(0,dtfil%filstat,iexit,level,'exit ') 731 732 call timab(708,2,tsec) 733 call timab(700,2,tsec) 734 735 DBG_EXIT("COLL") 736 737 end subroutine gstateimg