Line data Source code
1 : !!****m* ABINIT/m_gstateimg
2 : !! NAME
3 : !! m_gstateimg
4 : !!
5 : !! FUNCTION
6 : !!
7 : !! COPYRIGHT
8 : !! Copyright (C) 1998-2026 ABINIT group (XG, AR, GG, MT)
9 : !! This file is distributed under the terms of the
10 : !! GNU General Public License, see ~abinit/COPYING
11 : !! or http://www.gnu.org/copyleft/gpl.txt .
12 : !!
13 : !! SOURCE
14 :
15 : #if defined HAVE_CONFIG_H
16 : #include "config.h"
17 : #endif
18 :
19 : #include "abi_common.h"
20 :
21 : module m_gstateimg
22 :
23 : use defs_basis
24 : use defs_wvltypes
25 : use defs_rectypes
26 : use m_abicore
27 : use m_abihist
28 : use m_mep
29 : use m_ga
30 : use m_use_ga
31 : use m_pimd
32 : use m_xmpi
33 : use m_errors
34 : use m_rec
35 : use m_args_gs
36 : use m_results_gs
37 : use m_results_img
38 : use m_scf_history
39 : use m_io_redirect
40 : use m_m1geo
41 : use m_abimover
42 : use m_yaml
43 : use m_dtfil
44 :
45 : use defs_datatypes, only : pseudopotential_type
46 : use defs_abitypes, only : MPI_type
47 : use m_time, only : timab
48 : use m_geometry, only : mkradim, mkrdim, fcart2gred, xred2xcart, metric
49 : use m_specialmsg, only : specialmsg_mpisum
50 : use m_libpaw_tools, only : libpaw_spmsg_mpisum
51 : use m_pawang, only : pawang_type
52 : use m_pawrad, only : pawrad_type
53 : use m_pawtab, only : pawtab_type
54 : use m_gstate, only : gstate
55 : use m_predtk, only : prtxvf
56 : use m_precpred_1geo, only : precpred_1geo
57 : use m_pred_simple, only : prec_simple
58 :
59 : #if defined HAVE_BIGDFT
60 : use BigDFT_API, only: mpi_environment_set
61 : #endif
62 :
63 : implicit none
64 :
65 : private
66 : !!***
67 :
68 : public :: gstateimg
69 : !!***
70 :
71 : contains
72 : !!***
73 :
74 : !!****f* ABINIT/gstateimg
75 : !! NAME
76 : !! gstateimg
77 : !!
78 : !! FUNCTION
79 : !! Routine for conducting DFT calculations for a set of (dynamical) images
80 : !!
81 : !! INPUTS
82 : !! codvsn=code version
83 : !! cpui=initial CPU time
84 : !! nimage=number of images of the cell (treated by current proc)
85 : !! === Optional arguments (needed when nimage>1) ===
86 : !! filnam(5)=character strings giving file names
87 : !! filstat=character strings giving name of status file
88 : !! idtset=index of the dataset
89 : !! jdtset(0:ndtset)=actual index of the datasets
90 : !! ndtset=number of datasets
91 : !!
92 : !! OUTPUT
93 : !! etotal_img=total energy, for each image
94 : !! fcart_img(3,natom,nimage)=forces, in cartesian coordinates, for each image
95 : !! gred_img(3,natom,nimage)=gradient of E wrt nuclear positions, in reduced coordinates, for each image
96 : !! intgres_img(nspden,natom,nimage)=gradient wrt constraints, for each image
97 : !! npwtot(nkpt) = total number of plane waves at each k point
98 : !! strten_img(6,nimage)=stress tensor, for each image
99 : !!
100 : !! SIDE EFFECTS
101 : !! acell_img(3,nimage)=unit cell length scales (bohr), for each image
102 : !! amu_img(ntypat,nimage)=value of mass for each atomic type, for each image
103 : !! dtfil <type(datafiles_type)>=variables related to files
104 : !! dtset <type(dataset_type)>=all input variables in this dataset
105 : !! | mband =maximum number of bands (IN)
106 : !! | mgfft =maximum single fft dimension (IN)
107 : !! | mkmem =maximum number of k points which can fit in core memory (IN)
108 : !! | mpw =maximum number of planewaves in basis sphere (large number) (IN)
109 : !! | natom =number of atoms in unit cell (IN)
110 : !! | nfft =(effective) number of FFT grid points (for this processor) (IN)
111 : !! | nkpt =number of k points (IN)
112 : !! | nspden=number of spin-density components (IN)
113 : !! | nsppol=number of channels for spin-polarization (1 or 2) (IN)
114 : !! | nsym =number of symmetry elements in space group
115 : !! iexit= exit flag
116 : !! mixalch_img(npspalch,ntypalch,nimage)=value of alchemical mixing factors,for each image
117 : !! mpi_enreg=MPI-parallelisation information (some already initialized,
118 : !! some others to be initialized here)
119 : !! occ_img(mband*nkpt*nsppol,nimage) = occupation number for each band and k, for each image
120 : !! pawang <type(pawang_type)>=paw angular mesh and related data
121 : !! pawrad(ntypat*usepaw) <type(pawrad_type)>=paw radial mesh and related data
122 : !! pawtab(ntypat*usepaw) <type(pawtab_type)>=paw tabulated starting data
123 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
124 : !! Before entering the first time in gstateimg, a significant part of
125 : !! psps has been initialized :
126 : !! the integers dimekb,lmnmax,lnmax,mpssang,mpssoang,mpsso,mgrid,
127 : !! ntypat,n1xccc,usepaw,useylm, and the arrays dimensioned to npsp
128 : !! All the remaining components of psps are to be initialized in the call
129 : !! to pspini .
130 : !! The next time the code enters gstateimg, psps might be identical to the
131 : !! one of the previous dtset, in which case, no reinitialisation is scheduled
132 : !! in pspini.f .
133 : !! rprim_img(3,3,nimage)=dimensionless real space primitive translations, for each image
134 : !! vel_cell_img(3,3,nimage)=value of cell parameters velocities, for each image
135 : !! vel_img(3,natom,nimage)=value of atomic velocities,for each image
136 : !! xred_img(3,natom,nimage) = reduced atomic coordinates, for each image
137 : !!
138 : !! NOTES
139 : !! USE OF FFT GRIDS:
140 : !! =================
141 : !! In case of PAW:
142 : !! ---------------
143 : !! Two FFT grids are used:
144 : !! - A "coarse" FFT grid (defined by ecut)
145 : !! for the application of the Hamiltonian on the plane waves basis.
146 : !! It is defined by nfft, ngfft, mgfft, ...
147 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
148 : !! are expressed on this grid.
149 : !! - A "fine" FFT grid (defined) by ecutdg)
150 : !! for the computation of the density inside PAW spheres.
151 : !! It is defined by nfftf, ngfftf, mgfftf, ...
152 : !! Total density, potentials, ...
153 : !! are expressed on this grid.
154 : !! In case of norm-conserving:
155 : !! ---------------------------
156 : !! - Only the usual FFT grid (defined by ecut) is used.
157 : !! It is defined by nfft, ngfft, mgfft, ...
158 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
159 : !! are set equal to (nfft,ngfft,mgfft) in that case.
160 : !! In case of wavelets:
161 : !! --------------------
162 : !! - Only the usual FFT grid (defined by wvl_crmult) is used.
163 : !! It is defined by nfft, ngfft, mgfft, ... This is strictly not
164 : !! an FFT grid since its dimensions are not suited for FFTs. They are
165 : !! defined by wvl_setngfft().
166 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
167 : !! are set equal to (nfft,ngfft,mgfft) in that case.
168 : !!
169 : !! TODO
170 : !! Not yet possible to use restartxf in parallel when localrdwf==0
171 : !!
172 : !! SOURCE
173 :
174 7950 : subroutine gstateimg(acell_img,amu_img,codvsn,cpui,dtfil,dtset,etotal_img,fcart_img,&
175 3975 : & gred_img,iexit,intgres_img,mixalch_img,mpi_enreg,nimage,npwtot,occ_img,&
176 3975 : & pawang,pawrad,pawtab,psps,&
177 3975 : & rprim_img,strten_img,vel_cell_img,vel_img,wvl,xred_img,&
178 3975 : & filnam,filstat,idtset,jdtset,ndtset) ! optional arguments
179 :
180 : !Arguments ------------------------------------
181 : !scalars
182 : integer,intent(in) :: nimage
183 : integer,optional,intent(in) :: idtset,ndtset
184 : integer,intent(inout) :: iexit
185 : real(dp),intent(in) :: cpui
186 : character(len=8),intent(in) :: codvsn
187 : character(len=fnlen),optional,intent(in) :: filstat
188 : type(MPI_type),intent(inout) :: mpi_enreg
189 : type(datafiles_type),target,intent(inout) :: dtfil
190 : type(dataset_type),target,intent(inout) :: dtset
191 : type(pawang_type),intent(inout) :: pawang
192 : type(pseudopotential_type),intent(inout) :: psps
193 : type(wvl_data),intent(inout) :: wvl
194 : !arrays
195 : integer,optional,intent(in) :: jdtset(:)
196 : integer,intent(out) :: npwtot(dtset%nkpt)
197 : character(len=fnlen),optional,intent(in) :: filnam(:)
198 : real(dp), intent(out) :: etotal_img(nimage),fcart_img(3,dtset%natom,nimage)
199 : real(dp), intent(out) :: gred_img(3,dtset%natom,nimage)
200 : real(dp), intent(out) :: intgres_img(dtset%nspden,dtset%natom,nimage)
201 : real(dp), intent(out) :: strten_img(6,nimage)
202 : real(dp),intent(inout) :: acell_img(3,nimage),amu_img(dtset%ntypat,nimage)
203 : real(dp),intent(inout) :: mixalch_img(dtset%npspalch,dtset%ntypalch,nimage)
204 : real(dp),intent(inout) :: occ_img(dtset%mband*dtset%nkpt*dtset%nsppol,nimage)
205 : real(dp),intent(inout) :: rprim_img(3,3,nimage),vel_cell_img(3,3,nimage),vel_img(3,dtset%natom,nimage)
206 : real(dp),intent(inout) :: xred_img(3,dtset%natom,nimage)
207 : type(pawrad_type),intent(inout) :: pawrad(psps%ntypat*psps%usepaw)
208 : type(pawtab_type),intent(inout) :: pawtab(psps%ntypat*psps%usepaw)
209 :
210 : !Local variables-------------------------------
211 : !Define file format for different type of files. Presently,
212 : !only one file format is supported for each type of files, but this might
213 : !change soon ...
214 : !2 for wavefunction file, new format (version 2.0 and after) (fform) NOT USED
215 : !52 for density rho(r) (fformr)
216 : !102 for potential V(r) file. (fformv) NOT USED
217 : !scalars
218 : integer,parameter :: formeig=0,level=100,ndtpawuj=0,response=0
219 : integer :: history_size,idelta,idynimage,ierr,ifirst
220 : integer :: ii,iimage,ih,itimimage,itimimage_eff,itimimage_gstate,itimimage_prev,ndynimage,nocc
221 : integer :: ntimimage,ntimimage_stored,ntimimage_max
222 : logical :: check_conv,compute_all_images,compute_static_images
223 : logical :: isVused,isARused,is_master,is_mep,is_pimd
224 : logical :: call_predictor,use_hist,use_hist_prev
225 : real(dp) :: delta_energy,dtion
226 : character(len=500) :: hist_filename,msg
227 : type(args_gs_type) :: args_gs
228 3975 : type(mep_type) :: mep_param
229 3975 : type(ga_type) :: ga_param
230 3975 : type(m1geo_type) :: m1geo_param
231 3975 : type(pimd_type) :: pimd_param
232 : !arrays
233 3975 : integer,allocatable :: list_dynimage(:),scf_initialized(:)
234 : character(len=60),parameter :: imagealgo_str(0:13)=(/ &
235 : & 'IMAGE COPY ',& ! 0
236 : & 'IMAGE STEEPEST DESCENT ',& ! 1
237 : & 'STRING METHOD ',& ! 2
238 : & 'METADYNAMICS ',& ! 3
239 : & 'GENETIC ALGORITHM ',& ! 4
240 : & 'NUDGED ELASTIC BAND ',& ! 5
241 : & 'LINEAR COMBINATION OF CONSTRAINED DFT ENERGIES ',& ! 6
242 : & ' ',& ! 7
243 : & ' ',& ! 8
244 : & 'PATH-INTEGRAL MOLECULAR DYNAMICS (LANGEVIN) ',& ! 9
245 : & 'PATH-INTEGRAL MOLECULAR DYNAMICS (QUANTUM THERMAL BATH) ',& ! 10
246 : & ' ',& ! 11
247 : & ' ',& ! 12
248 : & 'PATH-INTEGRAL MOLECULAR DYNAMICS (CHAIN OF THERMOSTATS) '/) ! 13
249 : character(len=24),parameter :: stgalgo_str(0:2)=(/ &
250 : & 'ORIGINAL ALGO. ',& ! 0
251 : & 'SIMPLIFIED + EQUAL ARC ',& ! 1
252 : & 'SIMPLIFIED + ENERGY-WGTH'/) ! 2
253 : character(len=20),parameter :: nebalgo_str(0:2)=(/ &
254 : & 'ORIGINAL ALGO. ',& ! 0
255 : & 'IMPROVED TANGENT ',& ! 1
256 : & 'CLIMBING IMAGE '/) ! 2
257 : character(len=20),parameter :: mepsolver_str(0:4)=(/ &
258 : & 'STEEPEST-DESCENT ',& ! 0
259 : & 'QUICK-MIN OPT. ',& ! 1
260 : & 'L-BFGS ',& ! 2
261 : & 'GL-BFGS ',& ! 3
262 : & 'ORDER 4 RUNGE-KUTTA '/) ! 4
263 : real(dp) :: acell(3),rprim(3,3),rprimd(3,3),tsec(2),vel_cell(3,3)
264 3975 : real(dp),allocatable :: amass(:,:),occ(:),vel(:,:),xred(:,:)
265 : !real(dp),pointer :: vel_prev(:,:)
266 3975 : type(abihist),allocatable :: hist(:),hist_prev(:)
267 3975 : type(results_img_type),pointer :: results_img(:,:),res_img(:)
268 3975 : type(scf_history_type),allocatable :: scf_history(:)
269 : !type(abiforstr) :: preconforstr ! Preconditioned forces and stress ... Only needed to deallocate an internal matrix in prec_simple
270 :
271 : ! ***********************************************************************
272 :
273 : DBG_ENTER("COLL")
274 :
275 3975 : call timab(1200,1,tsec)
276 3975 : call timab(1203,3,tsec)
277 :
278 : !Arguments check
279 3975 : if (dtset%nimage>1) then
280 : if ((.not.present(filnam)).or.(.not.present(filnam)).or.(.not.present(idtset)).or.&
281 50 : & (.not.present(ndtset)).or.(.not.present(jdtset))) then
282 : write(msg,'(3a)') &
283 0 : & 'When nimage>1, all the following argument should be present:',ch10,&
284 0 : & 'filnam, filstat, idtset, ndtset, jdtset !'
285 0 : ABI_BUG(msg)
286 : end if
287 : end if
288 :
289 : !Set flag for the effective computation (once) of static images
290 : !For the time being only set when parallelization is activated
291 : !Note: if you modify this flag, do not forget to change it in outvars and outvar1
292 3975 : compute_static_images=(dtset%istatimg>0)
293 :
294 : !Prepare the calculation, by computing flags and dimensions
295 3975 : is_pimd=(dtset%imgmov==9.or.dtset%imgmov==10.or.dtset%imgmov==13)
296 3975 : is_mep =(dtset%imgmov==1.or.dtset%imgmov== 2.or.dtset%imgmov== 5)
297 3975 : ntimimage=dtset%ntimimage
298 3975 : ntimimage_stored=ntimimage;if(is_pimd)ntimimage_stored=2
299 3975 : nocc=dtset%mband*dtset%nkpt*dtset%nsppol
300 3975 : is_master=(mpi_enreg%me_cell==0.and.mpi_enreg%me_img==0)
301 3975 : delta_energy=zero
302 :
303 : !Management of dynamics/relaxation history (positions, forces, stresses, ...)
304 3975 : use_hist=(dtset%imgmov/=0.and.nimage>0) ; use_hist_prev=.false.
305 3975 : isVused=is_pimd;isARused=(dtset%optcell/=0)
306 3975 : if (use_hist) then
307 : !Read history from file (and broadcast if MPI)
308 43 : use_hist_prev=(dtset%restartxf==-1.and.nimage>0)
309 43 : hist_filename=trim(dtfil%filnam_ds(4))//'_HIST.nc'
310 43 : if (use_hist_prev)then
311 8 : ABI_MALLOC(hist_prev,(nimage))
312 1 : if (mpi_enreg%me_cell==0) then
313 : call read_md_hist_img(hist_filename,hist_prev,isVused,isARused,&
314 1 : & imgtab=mpi_enreg%my_imgtab)
315 : end if
316 1 : call abihist_bcast(hist_prev,0,mpi_enreg%comm_cell)
317 1 : if (nimage>0) then
318 6 : if (any(hist_prev(:)%mxhist/=hist_prev(1)%mxhist)) then
319 0 : msg='History problem: all images should have the same number of time steps!'
320 0 : ABI_ERROR(msg)
321 : end if
322 1 : use_hist_prev=(hist_prev(1)%mxhist>0)
323 1 : if (use_hist_prev) ntimimage=ntimimage+hist_prev(1)%mxhist
324 : end if
325 1 : if (.not.use_hist_prev) then
326 0 : call abihist_free(hist_prev)
327 0 : ABI_FREE(hist_prev)
328 : end if
329 : end if
330 : !Initialize a variable to write the history
331 395 : ABI_MALLOC(hist,(nimage))
332 43 : call abihist_init(hist,dtset%natom,ntimimage,isVused,isARused)
333 : end if ! imgmov/=0
334 :
335 : !Various allocations
336 11925 : ABI_MALLOC(occ,(nocc))
337 11925 : ABI_MALLOC(vel,(3,dtset%natom))
338 7950 : ABI_MALLOC(xred,(3,dtset%natom))
339 :
340 : !Select dynamical images
341 11925 : ABI_MALLOC(list_dynimage,(dtset%ndynimage))
342 3975 : ndynimage=0
343 8210 : do iimage=1,nimage
344 4235 : ii=mpi_enreg%my_imgtab(iimage)
345 8210 : if (dtset%dynimage(ii)==1) then
346 4175 : ndynimage=ndynimage+1
347 4175 : list_dynimage(ndynimage)=iimage
348 : end if
349 : end do
350 :
351 : !Management of SCF history (density/WF predictions from one time step to another)
352 16160 : ABI_MALLOC(scf_history,(nimage))
353 11925 : ABI_MALLOC(scf_initialized,(nimage))
354 8210 : scf_initialized=0
355 3975 : history_size=-1
356 3975 : if (dtset%ntimimage<=1) then
357 3932 : if (dtset%usewvl==0.and.dtset%ionmov>0.and. (abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) then
358 98 : history_size=2
359 98 : if(dtset%extrapwf==2) history_size=3
360 : end if
361 : else
362 43 : if (abs(dtset%densfor_pred)==2.or.abs(dtset%densfor_pred)==3) history_size=0
363 43 : if (dtset%imgwfstor==1) history_size=1
364 141 : if (dtset%usewvl==0.and.(abs(dtset%densfor_pred)==5.or.abs(dtset%densfor_pred)==6)) history_size=2
365 : end if
366 8210 : do iimage=1,nimage
367 4235 : call scf_history_nullify(scf_history(iimage))
368 8210 : scf_history(iimage)%history_size=history_size
369 : end do
370 :
371 : !In some cases, need amass variable
372 3975 : if (use_hist) then
373 172 : ABI_MALLOC(amass,(dtset%natom,nimage))
374 309 : do iimage=1,nimage
375 737 : if (any(amu_img(:,iimage)/=amu_img(:,1))) then
376 0 : ABI_ERROR('Problem with variable masses!')
377 : end if
378 1713 : amass(:,iimage)=amu_emass*amu_img(dtset%typat(:),iimage)
379 : end do
380 : end if
381 :
382 : !In the case of the 4th-order Runge-Kutta solver,
383 : !one must have a number of step multiple of 4.
384 3975 : ntimimage_max=ntimimage;idelta=1
385 3975 : if (dtset%imgmov==2.and.mep_param%mep_solver==4) then
386 0 : ntimimage_max=4*(ntimimage_max/4)
387 0 : idelta=4
388 : end if
389 :
390 : !MEP search: fill in eventually the data structure mep_param
391 3975 : call mep_init(dtset,mep_param)
392 :
393 : !GA search: fill in eventually the data structure ga_param
394 3975 : call ga_init(dtset,ga_param)
395 :
396 : !Move 1GEO approach: fill the data structure m1geo_param
397 3975 : call m1geo_init(dtfil,dtset,m1geo_param)
398 :
399 : !PIMD: fill in the data structure pimd_param and init velocities
400 3975 : call pimd_init(dtset,pimd_param,is_master)
401 3975 : dtion=one;if (is_pimd) dtion=pimd_param%dtion
402 :
403 : !Initialization of main variables
404 26525 : ABI_MALLOC(results_img,(nimage,ntimimage_stored))
405 8291 : do itimimage=1,ntimimage_stored
406 4316 : res_img => results_img(:,itimimage)
407 : call init_results_img(dtset%natom,dtset%npspalch,dtset%nspden,dtset%nsppol,dtset%ntypalch,&
408 4316 : & dtset%ntypat,res_img)
409 14600 : do iimage=1,nimage
410 25236 : res_img(iimage)%acell(:) =acell_img(:,iimage)
411 15700 : res_img(iimage)%amu(:) =amu_img(:,iimage)
412 6384 : res_img(iimage)%mixalch(:,:) =mixalch_img(:,:,iimage)
413 82017 : res_img(iimage)%rprim(:,:) =rprim_img(:,:,iimage)
414 64453 : res_img(iimage)%xred(:,:) =xred_img(:,:,iimage)
415 64453 : res_img(iimage)%vel(:,:) =vel_img(:,:,iimage)
416 86333 : res_img(iimage)%vel_cell(:,:)=vel_cell_img(:,:,iimage)
417 : end do
418 : end do
419 :
420 : !Set Number of degrees Of Freedom for PIMD algorithms
421 3975 : if(use_hist) then
422 309 : hist(:)%ndof=3*dtset%natom ! Init ndof to 3N (no iatfix for PIMD)
423 43 : if(dtset%imgmov==9.or.dtset%imgmov==10) then
424 : if(pimd_param%pitransform==1.or.pimd_param%pitransform==2.or.&
425 4 : & pimd_param%constraint==1.or.pimd_param%optcell==2) then
426 21 : hist(:)%ndof=hist(:)%ndof-3
427 : end if
428 : end if
429 : end if
430 :
431 3975 : call timab(1203,2,tsec)
432 :
433 : !-----------------------------------------------------------------------------------------
434 : !Big loop on the propagation of all images
435 3975 : itimimage_eff=1 ; itimimage_gstate=1
436 4271 : do itimimage=1,ntimimage
437 :
438 4271 : res_img => results_img(:,itimimage_eff)
439 4271 : call_predictor=(ntimimage>1)
440 :
441 : ! If history is activated and if current image is inside it: do not compute anything
442 4271 : if (use_hist_prev) then
443 35 : if (all(hist_prev(:)%ihist<=hist_prev(:)%mxhist)) then
444 30 : do iimage=1,nimage
445 25 : ih=hist_prev(iimage)%ihist
446 25 : call abihist_copy(hist_prev(iimage),hist(iimage))
447 25 : call mkradim(hist_prev(iimage)%acell(:,ih),rprim,hist_prev(iimage)%rprimd(:,:,ih))
448 100 : res_img(iimage)%acell(:)=hist_prev(iimage)%acell(:,ih)
449 325 : res_img(iimage)%rprim(:,:)=rprim
450 625 : res_img(iimage)%xred(:,:)=hist_prev(iimage)%xred(:,:,ih)
451 625 : res_img(iimage)%vel(:,:)=hist_prev(iimage)%vel(:,:,ih)
452 325 : res_img(iimage)%vel_cell(:,:)=hist_prev(iimage)%vel_cell(:,:,ih)
453 625 : res_img(iimage)%results_gs%fcart(:,:)=hist_prev(iimage)%fcart(:,:,ih)
454 175 : res_img(iimage)%results_gs%strten(:)=hist_prev(iimage)%strten(:,ih)
455 25 : res_img(iimage)%results_gs%etotal=hist_prev(iimage)%etot(ih)
456 25 : res_img(iimage)%results_gs%energies%entropy=hist_prev(iimage)%entropy(ih)
457 : call fcart2gred(res_img(iimage)%results_gs%fcart,res_img(iimage)%results_gs%gred,&
458 25 : hist_prev(iimage)%rprimd(:,:,ih),dtset%natom)
459 30 : hist_prev(iimage)%ihist=hist_prev(iimage)%ihist+1
460 : end do
461 : !PI-QTB: skip a record in random force file
462 5 : if (pimd_param%use_qtb==1) call pimd_skip_qtb(pimd_param)
463 : !call_predictor=.false.
464 : goto 110 ! This is temporary
465 : end if
466 : end if
467 :
468 4266 : call timab(1204,1,tsec)
469 4266 : call localfilnam(mpi_enreg%comm_img,mpi_enreg%comm_cell,mpi_enreg%comm_world,filnam,'_IMG',dtset%nimage)
470 4266 : compute_all_images=(compute_static_images.and.itimimage==1)
471 :
472 : ! Print title for time step
473 4266 : if (dtset%nimage>1.or.dtset%ntimimage>1) then
474 342 : if (dtset%prtvolimg<2) then
475 274 : msg=ch10;if (itimimage >1) write(msg,'(2a)') ch10,ch10
476 274 : write(msg,'(5a)') trim(msg),&
477 274 : '================================================================================',&
478 548 : ch10,' ',trim(imagealgo_str(dtset%imgmov))
479 : else
480 68 : msg='';if (itimimage >1) msg=ch10
481 68 : write(msg,'(5a)') trim(msg),&
482 68 : '--------------------------------------------------------------------------------',&
483 136 : ch10,' ',trim(imagealgo_str(dtset%imgmov))
484 : end if
485 342 : if (dtset%imgmov==2) then
486 206 : write(msg,'(6a)') trim(msg),' (',trim(stgalgo_str(mep_param%string_algo)),' + ',&
487 412 : trim(mepsolver_str(mep_param%mep_solver)),')'
488 : end if
489 342 : if (dtset%imgmov==5) then
490 46 : ii=merge(mep_param%neb_algo,1,mep_param%neb_algo/=2.or.itimimage>=mep_param%cineb_start)
491 46 : write(msg,'(6a)') trim(msg),' (',trim(nebalgo_str(ii)),' + ',&
492 92 : trim(mepsolver_str(mep_param%mep_solver)),')'
493 : end if
494 342 : if (dtset%ntimimage==1) write(msg,'(2a)') trim(msg),' FOR 1 TIME STEP'
495 342 : if (dtset%ntimimage >1) write(msg,'(2a,i5)') trim(msg),' - TIME STEP ',itimimage
496 342 : if (dtset%prtvolimg<2) then
497 274 : write(msg,'(3a)') trim(msg),ch10,&
498 548 : & '================================================================================'
499 : end if
500 342 : call wrtout(ab_out ,msg,'COLL')
501 342 : call wrtout(std_out,msg,'PERS')
502 :
503 342 : call yaml_iterstart('timimage', itimimage, dev_null, dtset%use_yaml)
504 : end if
505 :
506 4266 : if (dtset%use_yaml == 1) call yaml_iterstart('timimage', itimimage, ab_out, dtset%use_yaml)
507 :
508 4266 : call timab(1204,2,tsec)
509 :
510 : ! Loop on the dynamical images
511 4266 : idynimage=0
512 10347 : do iimage=1,nimage
513 :
514 6081 : ii=mpi_enreg%my_imgtab(iimage)
515 6081 : if (dtset%dynimage(ii)==1) idynimage=idynimage+1
516 :
517 : ! Compute static image only at first time step
518 6081 : if (dtset%dynimage(ii)==1.or.compute_all_images) then
519 :
520 5629 : call timab(1205,1,tsec)
521 :
522 : ! Change file names according to image index (if nimage>1)
523 5629 : if (dtset%nimage>1) then
524 1700 : call dtfil%init(dtset,filnam,filstat,idtset,jdtset,mpi_enreg,ndtset,image_index=ii)
525 1700 : if (itimimage>1) then
526 1399 : dtfil%ireadwf=0;dtfil%ireadden=0;dtfil%ireadkden=0
527 : end if
528 1700 : call yaml_iterstart('image', iimage, dev_null, 0)
529 : end if
530 :
531 5629 : if (dtset%use_yaml == 1) call yaml_iterstart('image', iimage, ab_out, dtset%use_yaml)
532 :
533 : ! Redefine output units
534 5629 : call localwrfile(mpi_enreg%comm_cell,ii,dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg)
535 :
536 : ! Print title for image
537 5629 : if (dtset%nimage>1.and.(dtset%prtvolimg==0.or.do_write_log)) then
538 1324 : if (ii==1) write(msg,'(a)' ) ch10
539 1324 : if (ii >1) write(msg,'(2a)') ch10,ch10
540 1324 : write(msg,'(6a,i4,a,i4,3a)') trim(msg),&
541 1324 : '--------------------------------------------------------------------------------',ch10,&
542 1324 : ' ',trim(imagealgo_str(dtset%imgmov)),' - CELL # ',ii,'/',dtset%nimage,ch10,&
543 2648 : '--------------------------------------------------------------------------------',ch10
544 1324 : if (dtset%prtvolimg==0) call wrtout(ab_out ,msg,'COLL')
545 1324 : if (do_write_log) call wrtout(std_out,msg,'PERS')
546 : end if
547 :
548 22516 : acell(:) =res_img(iimage)%acell(:)
549 73177 : rprim(:,:) =res_img(iimage)%rprim(:,:)
550 58573 : vel(:,:) =res_img(iimage)%vel(:,:)
551 73177 : vel_cell(:,:)=res_img(iimage)%vel_cell(:,:)
552 58573 : xred(:,:) =res_img(iimage)%xred(:,:)
553 449988 : occ(:) =occ_img(:,iimage)
554 :
555 : call args_gs_init(args_gs, &
556 : res_img(iimage)%amu(:),dtset%cellcharge(ii),res_img(iimage)%mixalch(:,:),&
557 5629 : dtset%dmatpawu(:,:,:,:,ii),dtset%upawu(:,ii),dtset%jpawu(:,ii),dtset%rprimd_orig(:,:,ii))
558 :
559 5629 : call timab(1205,2,tsec)
560 :
561 : call gstate(args_gs,acell,codvsn,cpui,dtfil,dtset,iexit,scf_initialized(iimage),itimimage_gstate,&
562 : mpi_enreg,npwtot,occ,pawang,pawrad,pawtab,psps,&
563 : res_img(iimage)%results_gs,&
564 5629 : rprim,scf_history(iimage),vel,vel_cell,wvl,xred)
565 :
566 5629 : itimimage_gstate=itimimage_gstate+1
567 :
568 5629 : call timab(1206,1,tsec)
569 :
570 5629 : call args_gs_free(args_gs)
571 :
572 5629 : if (dtset%dynimage(ii)==1) then
573 22292 : res_img(iimage)%acell(:) =acell(:)
574 72449 : res_img(iimage)%rprim(:,:) =rprim(:,:)
575 58021 : res_img(iimage)%vel(:,:) =vel(:,:)
576 72449 : res_img(iimage)%vel_cell(:,:)=vel_cell(:,:)
577 58021 : res_img(iimage)%xred(:,:) =xred(:,:)
578 448816 : occ_img(:,iimage) =occ(:)
579 : end if
580 :
581 : ! check change of rprim and rewriting in hist
582 : ! check change of xred and rewriting in hist
583 :
584 : ! Close output units ; restore defaults
585 5629 : call localredirect(mpi_enreg%comm_cell,mpi_enreg%comm_world,dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg)
586 5629 : call timab(1206,2,tsec)
587 :
588 452 : else if (itimimage>1) then ! For static images, simply copy one time step to the other
589 448 : itimimage_prev=itimimage_eff-1;if (itimimage_prev<1) itimimage_prev=ntimimage_stored
590 : call copy_results_img(results_img(iimage,itimimage_prev), &
591 448 : & results_img(iimage,itimimage_eff ))
592 : end if
593 :
594 : ! Store results in hist datastructure
595 10347 : if (use_hist) then
596 2110 : ih=hist(iimage)%ihist
597 2110 : call mkrdim(res_img(iimage)%acell(:),res_img(iimage)%rprim(:,:),rprimd)
598 : call var2hist(res_img(iimage)%acell(:),hist(iimage),dtset%natom,&
599 2110 : & rprimd,res_img(iimage)%xred(:,:),.FALSE.)
600 2110 : call vel2hist(amass(:,iimage),hist(iimage),res_img(iimage)%vel(:,:),res_img(iimage)%vel_cell(:,:))
601 21094 : hist(iimage)%fcart(:,:,ih)=res_img(iimage)%results_gs%fcart(:,:)
602 14770 : hist(iimage)%strten(:,ih)=res_img(iimage)%results_gs%strten(:)
603 2110 : hist(iimage)%etot(ih)=res_img(iimage)%results_gs%etotal
604 2110 : hist(iimage)%entropy(ih)=res_img(iimage)%results_gs%energies%entropy
605 2110 : hist(iimage)%time(ih)=real(itimimage,kind=dp)*dtion
606 : end if
607 :
608 : end do ! iimage
609 :
610 4266 : if(mpi_enreg%paral_img==1)then
611 0 : call timab(1208,1,tsec)
612 0 : call xmpi_barrier(mpi_enreg%comm_img)
613 0 : call timab(1208,2,tsec)
614 : end if
615 :
616 4266 : call timab(1209,1,tsec)
617 :
618 : ! Output when images are used
619 4266 : if (dtset%nimage>1) then
620 : ! === 1st option: reduced outputs ===
621 337 : if (dtset%prtvolimg>0) then
622 : call prtimg(dtset%dynimage,imagealgo_str(dtset%imgmov),dtset%imgmov,ab_out,&
623 131 : & mpi_enreg,nimage,dtset%nimage,compute_all_images,dtset%prtvolimg,res_img)
624 : end if
625 : end if
626 :
627 : ! Manage log files when images are used
628 : call localrdfile(mpi_enreg%comm_img,mpi_enreg%comm_world,compute_all_images,&
629 4266 : & dtset%nimage,mpi_enreg%paral_img,dtset%prtvolimg,dyn=dtset%dynimage)
630 :
631 : ! Write hist datastructure in HIST file
632 : ! Note: for PIMD, writing is done later
633 : ! if (use_hist.and.mpi_enreg%me_cell==0) then
634 4266 : if (use_hist.and.mpi_enreg%me_cell==0.and.(.not.is_pimd)) then
635 190 : ifirst=merge(0,1,itimimage>1)
636 : call write_md_hist_img(hist,hist_filename,ifirst,itimimage,dtset%natom,dtset%ntypat,&
637 : & dtset%typat,amu_img(:,1),dtset%znucl,dtion,&
638 : & nimage=dtset%nimage,imgmov=dtset%imgmov,mdtemp=dtset%mdtemp,comm_img=mpi_enreg%comm_img,&
639 190 : & imgtab=mpi_enreg%my_imgtab)
640 : end if
641 :
642 : ! TESTS WHETHER ONE CONTINUES THE LOOP
643 : ! Here we calculate the change in energy, and exit if delta_energy < tolimg
644 35 : delta_energy=zero
645 : ! Doesn't check convergence in case of PIMD
646 4231 : check_conv=((.not.is_pimd).and.itimimage>1)
647 : ! In case of 4th-order Runge-Kutta, does check convergence every 4 steps
648 4266 : if (dtset%imgmov==2.and.mep_param%mep_solver==4) then
649 0 : check_conv=(mod(itimimage,4)==0.and.itimimage>4)
650 : end if
651 4266 : if (check_conv) then
652 1497 : do idynimage=1,ndynimage
653 1234 : iimage=list_dynimage(idynimage)
654 : delta_energy=delta_energy &
655 : & +abs(results_img(iimage,itimimage)%results_gs%etotal &
656 1497 : & -results_img(iimage,itimimage-idelta)%results_gs%etotal)
657 : end do
658 263 : if (mpi_enreg%paral_img==1) then
659 0 : call xmpi_sum(delta_energy,mpi_enreg%comm_img,ierr)
660 : end if
661 263 : delta_energy=delta_energy/dtset%ndynimage
662 263 : if (delta_energy<dtset%tolimg) then
663 33 : if (dtset%prtvolimg<2) then
664 26 : write(msg,'(5a,i5,6a,es11.3,a,es11.3,2a)') ch10,ch10,&
665 26 : & '================================================================================',ch10,&
666 26 : & ' At time step ',itimimage,ch10,&
667 26 : & ' ',trim(imagealgo_str(dtset%imgmov)),' has reached energy convergence',ch10,&
668 26 : & ' with Average[Abs(Etotal(t)-Etotal(t-dt))]=',delta_energy,'<tolimg=',dtset%tolimg,ch10,&
669 52 : & '================================================================================'
670 : else
671 7 : write(msg,'(4a,i5,6a,es11.3,a,es11.3)') ch10,&
672 7 : & '--------------------------------------------------------------------------------',ch10,&
673 7 : & ' At time step ',itimimage,ch10,&
674 7 : & ' ',trim(imagealgo_str(dtset%imgmov)),' has reached energy convergence',ch10,&
675 14 : & ' with Average[Abs(Etotal(t)-Etotal(t-dt))]=',delta_energy,'<tolimg=',dtset%tolimg
676 : end if
677 99 : call wrtout([std_out, ab_out] ,msg,'COLL')
678 33 : call timab(1209,2,tsec) ! This is the first place where counter 1209 is stopped.
679 33 : exit ! exit itimimage
680 : end if
681 : end if
682 :
683 : !Temporary statement
684 : ! 110 continue
685 :
686 : ! Dont call the predictor at last time step (except for PIMD)
687 4233 : if (itimimage>=ntimimage_max) call_predictor=(call_predictor.and.is_pimd)
688 :
689 : ! Predict the next value of the images
690 4233 : if (call_predictor) then
691 : call predictimg(delta_energy,imagealgo_str(dtset%imgmov),dtset%imgmov,itimimage,&
692 : & itimimage_eff,list_dynimage,ga_param,mep_param,mpi_enreg,m1geo_param,dtset%natom,ndynimage,&
693 298 : & nimage,dtset%nimage,ntimimage_stored,pimd_param,dtset%prtvolimg,results_img)
694 : end if
695 :
696 : ! Write hist datastructure in HIST file in case of PIMD
697 : ! Note : velocities have been updated
698 : #if defined HAVE_NETCDF
699 4233 : if (use_hist.and.is_pimd) then
700 240 : do iimage=1,nimage
701 205 : ih=hist(iimage)%ihist
702 : call vel2hist(amass(:,iimage),hist(iimage),results_img(iimage,itimimage_eff)%vel(:,:),&
703 240 : & results_img(iimage,itimimage_eff)%vel_cell(:,:))
704 : end do
705 35 : if (mpi_enreg%me_cell==0) then
706 35 : ifirst=merge(0,1,itimimage>1)
707 : call write_md_hist_img(hist,hist_filename,ifirst,itimimage,dtset%natom,dtset%ntypat,&
708 : & dtset%typat,amu_img(:,1),dtset%znucl,dtion,&
709 : & nimage=dtset%nimage,imgmov=dtset%imgmov,mdtemp=dtset%mdtemp,comm_img=mpi_enreg%comm_img,&
710 35 : & imgtab=mpi_enreg%my_imgtab)
711 : end if
712 : end if
713 : #endif
714 :
715 : !Temporary statement
716 : 110 continue
717 :
718 : ! Increment indexes
719 4238 : if (itimimage>=ntimimage_max) exit
720 296 : itimimage_eff=itimimage_eff+1;if (itimimage_eff>ntimimage_stored) itimimage_eff=1
721 296 : if (use_hist) then
722 2164 : do iimage=1,nimage
723 2164 : hist(iimage)%ihist=hist(iimage)%ihist+1
724 : end do
725 : end if
726 :
727 4238 : call timab(1209,2,tsec) ! This is the second place where counter 1209 is stopped.
728 :
729 : end do ! itimimage
730 : !-----------------------------------------------------------------------------------------
731 :
732 3975 : call timab(1210,1,tsec)
733 :
734 : !Copy the results of the computation in the appropriate arguments of the routine
735 8210 : do iimage=1,nimage
736 4235 : ii=mpi_enreg%my_imgtab(iimage)
737 8210 : if (dtset%dynimage(ii)==1) then
738 16700 : acell_img(:,iimage) =results_img(iimage,itimimage_eff)%acell(:)
739 9609 : amu_img(:,iimage) =results_img(iimage,itimimage_eff)%amu(:)
740 4250 : mixalch_img(:,:,iimage) =results_img(iimage,itimimage_eff)%mixalch(:,:)
741 54275 : rprim_img(:,:,iimage) =results_img(iimage,itimimage_eff)%rprim(:,:)
742 44015 : vel_img(:,:,iimage) =results_img(iimage,itimimage_eff)%vel(:,:)
743 54275 : vel_cell_img(:,:,iimage)=results_img(iimage,itimimage_eff)%vel_cell(:,:)
744 44015 : xred_img(:,:,iimage) =results_img(iimage,itimimage_eff)%xred(:,:)
745 4175 : etotal_img(iimage) =results_img(iimage,itimimage_eff)%results_gs%etotal
746 44015 : fcart_img(:,:,iimage) =results_img(iimage,itimimage_eff)%results_gs%fcart(:,:)
747 44015 : gred_img(:,:,iimage) =results_img(iimage,itimimage_eff)%results_gs%gred(:,:)
748 26031 : intgres_img(:,:,iimage) =results_img(iimage,itimimage_eff)%results_gs%intgres(:,:)
749 29225 : strten_img(:,iimage) =results_img(iimage,itimimage_eff)%results_gs%strten(:)
750 60 : else if (compute_static_images) then
751 56 : etotal_img(iimage) =results_img(iimage,1)%results_gs%etotal
752 552 : fcart_img(:,:,iimage) =results_img(iimage,1)%results_gs%fcart(:,:)
753 552 : gred_img(:,:,iimage) =results_img(iimage,1)%results_gs%gred(:,:)
754 304 : intgres_img(:,:,iimage)=results_img(iimage,1)%results_gs%intgres(:,:)
755 392 : strten_img(:,iimage) =results_img(iimage,1)%results_gs%strten(:)
756 : end if
757 : end do
758 :
759 : !When parallelizattion over images is activated, has to sum number of warnings
760 : !and comments written in log file
761 3975 : if (mpi_enreg%paral_img==1.and.mpi_enreg%comm_cell==0) then
762 0 : call specialmsg_mpisum(mpi_enreg%comm_img)
763 0 : call libpaw_spmsg_mpisum(mpi_enreg%comm_img)
764 : end if
765 :
766 :
767 : !Final deallocations
768 :
769 : !This call is needed to free internal storages in different routines (prec_simple, pred_bfgs ...)
770 3975 : if(dtset%imgmov==6)then
771 6 : m1geo_param%iexit=1
772 : call predictimg(delta_energy,imagealgo_str(dtset%imgmov),dtset%imgmov,itimimage,&
773 : & itimimage_eff,list_dynimage,ga_param,mep_param,mpi_enreg,m1geo_param,dtset%natom,ndynimage,&
774 6 : & nimage,dtset%nimage,ntimimage_stored,pimd_param,dtset%prtvolimg,results_img)
775 : endif
776 :
777 3975 : ABI_FREE(occ)
778 3975 : ABI_FREE(vel)
779 3975 : ABI_FREE(xred)
780 3975 : ABI_FREE(list_dynimage)
781 :
782 3975 : ABI_SFREE(amass)
783 :
784 8291 : do itimimage=1,ntimimage_stored
785 8291 : call destroy_results_img(results_img(:,itimimage))
786 : end do
787 3975 : ABI_FREE(results_img)
788 8210 : do iimage=1,nimage
789 8210 : call scf_history_free(scf_history(iimage))
790 : end do
791 8210 : ABI_FREE(scf_history)
792 3975 : ABI_FREE(scf_initialized)
793 3975 : if (allocated(hist_prev)) then
794 1 : call abihist_free(hist_prev)
795 6 : ABI_FREE(hist_prev)
796 : end if
797 3975 : if (allocated(hist)) then
798 43 : call abihist_free(hist)
799 309 : ABI_FREE(hist)
800 : end if
801 :
802 3975 : call mep_destroy(mep_param)
803 3975 : call ga_destroy(ga_param)
804 3975 : call m1geo_destroy(m1geo_param)
805 3975 : call pimd_destroy(pimd_param)
806 :
807 3975 : call timab(1210,2,tsec)
808 3975 : call timab(1200,2,tsec)
809 :
810 : DBG_EXIT("COLL")
811 :
812 7950 : end subroutine gstateimg
813 : !!***
814 :
815 : !!****f* ABINIT/prtimg
816 : !! NAME
817 : !! prtimg
818 : !!
819 : !! FUNCTION
820 : !! Print out results obtained by as ground-state calculation of
821 : !! an image of the cell. The printing format is condensed in order
822 : !! to facilitate the reading.
823 : !!
824 : !! INPUTS
825 : !! dynimage(nimagetot)=flags defining static/dynamic state of images
826 : !! imagealgo_str=name of the algorithm (with images) used
827 : !! imgmov=index of algorithm (with images) used
828 : !! iout=unit number for output
829 : !! mpi_enreg=MPI-parallelisation information
830 : !! nimage=number of images stored on current proc
831 : !! nimage_tot=total number of images (should be dtset%nimage)
832 : !! prt_all_images=true if all images have to be printed out (ignoring dynimage)
833 : !! prtvolimg=printing volume for each image
834 : !! <0 : nothing
835 : !! 0 : only a title
836 : !! 1 : energy, residuals, forces, stresses, velocities, atomic positions
837 : !! 2 : energy, residuals
838 : !! resimg(nimage) <type(results_img_type)>=results of the ground-state computations
839 : !! for all images treated by current proc
840 : !!
841 : !! OUTPUT
842 : !! (data written to unit iout)
843 : !!
844 : !! SOURCE
845 :
846 131 : subroutine prtimg(dynimage,imagealgo_str,imgmov,iout,mpi_enreg,nimage,nimage_tot,&
847 131 : & prt_all_images,prtvolimg,resimg)
848 :
849 : !Arguments ------------------------------------
850 : !scalars
851 : integer,intent(in) :: nimage_tot,dynimage(nimage_tot),imgmov,iout,nimage,prtvolimg !vz_d
852 : logical,intent(in) :: prt_all_images
853 : character(len=60),intent(in) :: imagealgo_str
854 : type(MPI_type),intent(in) :: mpi_enreg
855 : !arrays
856 : type(results_img_type),target,intent(inout) :: resimg(nimage)
857 :
858 : !Local variables-------------------------------
859 : !scalars
860 : integer :: ii,prtvel
861 : logical :: test_img
862 : real(dp) :: ucvol_img
863 : character(len=500) :: msg
864 : !arrays
865 131 : integer,allocatable :: iatfix_img(:,:)
866 131 : real(dp),allocatable :: gmet_img(:,:),gprimd_img(:,:),rmet_img(:,:),xcart_img(:,:)
867 131 : type(results_img_type),pointer :: resimg_all(:)
868 :
869 : ! ****************************************************************
870 :
871 : DBG_ENTER('COLL')
872 :
873 131 : if (prtvolimg<=0) return
874 131 : if (mpi_enreg%me_cell/=0) return
875 :
876 : !Gather data
877 59 : if (prtvolimg==1.or.prtvolimg==2) then
878 59 : test_img=(nimage_tot/=1.and.mpi_enreg%paral_img==1)
879 : if (test_img) then
880 0 : if (mpi_enreg%me==0) then
881 0 : ABI_MALLOC(resimg_all,(nimage_tot))
882 : end if
883 : call gather_results_img(mpi_enreg,resimg,resimg_all,master=0,&
884 0 : & allgather=.false.,only_one_per_img=.true.)
885 : else
886 59 : resimg_all => resimg
887 : end if
888 : end if
889 :
890 : !===== First option for the printing volume ===
891 59 : if (prtvolimg==1.and.mpi_enreg%me==0) then
892 :
893 27 : prtvel=0;if (imgmov==0.or.imgmov==9.or.imgmov==10.or.imgmov==13) prtvel=1
894 :
895 216 : do ii=1,nimage_tot
896 216 : if (dynimage(ii)==1.or.prt_all_images) then
897 :
898 : ! Title
899 141 : write(msg,'(6a,i4,a,i4,2a)') ch10,&
900 141 : '----------------------------------------------------------------------',ch10,&
901 141 : ' ',trim(imagealgo_str),' - CELL # ',ii,'/',nimage_tot,ch10,&
902 282 : '----------------------------------------------------------------------'
903 141 : call wrtout(iout,msg,'COLL')
904 :
905 : ! Total energy
906 141 : write(msg,'(2a,es20.12)') ch10,' Total energy for the cell [Ha]: ',resimg_all(ii)%results_gs%etotal
907 141 : call wrtout(iout,msg,'COLL')
908 :
909 : ! Residuals of the SCF cycle
910 141 : write(msg,'(3a,4(a,es16.8,a))') ch10,&
911 141 : ' Residuals from SCF cycle: ',ch10,&
912 141 : ' Total energy difference =',resimg_all(ii)%results_gs%deltae,ch10,&
913 141 : ' Maximal forces difference =',resimg_all(ii)%results_gs%diffor,ch10,&
914 141 : ' Max. residual of wave-functions=',resimg_all(ii)%results_gs%residm,ch10,&
915 282 : ' Density/potential residual (^2)=',resimg_all(ii)%results_gs%res2,ch10
916 141 : call wrtout(iout,msg,'COLL')
917 :
918 : ! Cell parameters
919 141 : ABI_MALLOC(rmet_img,(3,3))
920 141 : ABI_MALLOC(gmet_img,(3,3))
921 141 : ABI_MALLOC(gprimd_img,(3,3))
922 141 : call metric(gmet_img,gprimd_img,iout,rmet_img,resimg_all(ii)%rprim,ucvol_img)
923 141 : ABI_FREE(rmet_img)
924 141 : ABI_FREE(gmet_img)
925 141 : ABI_FREE(gprimd_img)
926 :
927 : ! Positions, forces and velocities
928 423 : ABI_MALLOC(iatfix_img,(3,resimg_all(ii)%natom))
929 423 : ABI_MALLOC(xcart_img,(3,resimg_all(ii)%natom))
930 1269 : iatfix_img=0
931 141 : call xred2xcart(resimg_all(ii)%natom,resimg_all(ii)%rprim,xcart_img,resimg_all(ii)%xred)
932 : call prtxvf(resimg_all(ii)%results_gs%fcart,resimg_all(ii)%results_gs%gred,&
933 : & iatfix_img,iout,resimg_all(ii)%natom,prtvel,&
934 141 : & resimg_all(ii)%vel,xcart_img,resimg_all(ii)%xred)
935 141 : ABI_FREE(iatfix_img)
936 141 : ABI_FREE(xcart_img)
937 :
938 : ! Stress tensor
939 : write(msg, '(a,es12.4,a)' ) &
940 141 : & '-Cartesian components of stress tensor (GPa) [Pressure=',&
941 : & -(resimg_all(ii)%results_gs%strten(1)+resimg_all(ii)%results_gs%strten(2) &
942 282 : & +resimg_all(ii)%results_gs%strten(3))*HaBohr3_GPa/three,' GPa]'
943 141 : call wrtout(iout,msg,'COLL')
944 141 : write(msg, '(2(a,1p,e16.8))' ) '- sigma(1 1)=',resimg_all(ii)%results_gs%strten(1)*HaBohr3_GPa,&
945 282 : & ' sigma(3 2)=',resimg_all(ii)%results_gs%strten(4)*HaBohr3_GPa
946 141 : call wrtout(iout,msg,'COLL')
947 141 : write(msg, '(2(a,1p,e16.8))' ) '- sigma(2 2)=',resimg_all(ii)%results_gs%strten(2)*HaBohr3_GPa,&
948 282 : & ' sigma(3 1)=',resimg_all(ii)%results_gs%strten(5)*HaBohr3_GPa
949 141 : call wrtout(iout,msg,'COLL')
950 141 : write(msg, '(2(a,1p,e16.8))' ) '- sigma(3 3)=',resimg_all(ii)%results_gs%strten(3)*HaBohr3_GPa,&
951 282 : & ' sigma(2 1)=',resimg_all(ii)%results_gs%strten(6)*HaBohr3_GPa
952 282 : call wrtout(iout,msg,'COLL')
953 : end if
954 : end do
955 : end if
956 :
957 :
958 : !===== 2nd option for the printing volume ===
959 59 : if (prtvolimg==2.and.mpi_enreg%me==0) then
960 32 : write(msg,'(a,1x,a)') ch10,'Cell Total_energy[Ha] deltae diffor residm res2'
961 32 : call wrtout(iout,msg,'COLL')
962 246 : do ii=1,nimage_tot
963 246 : if (dynimage(ii)==1.or.prt_all_images) then
964 : write(msg,'(1x,i4,2x,es16.8,4(1x,es13.5))') &
965 166 : & ii,resimg_all(ii)%results_gs%etotal,resimg_all(ii)%results_gs%deltae,&
966 166 : & resimg_all(ii)%results_gs%diffor,resimg_all(ii)%results_gs%residm,&
967 332 : & resimg_all(ii)%results_gs%res2
968 166 : call wrtout(iout,msg,'COLL')
969 : end if
970 : end do
971 : end if
972 :
973 : !=====
974 59 : if (prtvolimg==1.or.prtvolimg==2) then
975 59 : if (test_img.and.mpi_enreg%me==0) then
976 0 : call destroy_results_img(resimg_all)
977 0 : ABI_FREE(resimg_all)
978 : end if
979 : nullify(resimg_all)
980 : end if
981 :
982 : DBG_EXIT('COLL')
983 :
984 131 : end subroutine prtimg
985 : !!***
986 :
987 : !!****f* ABINIT/predictimg
988 : !! NAME
989 : !! predictimg
990 : !!
991 : !! FUNCTION
992 : !! Given the past history of images, predict the new set of images
993 : !!
994 : !! INPUTS
995 : !! deltae=averaged energy difference used to control convergence over images
996 : !! imagealgo_str=name of the algorithm (with images) used
997 : !! imgmov=gives the algorithm to be used for prediction of new set of images
998 : !! itimimage=time index for image propagation (itimimage+1 is to be predicted here)
999 : !! itimimage_eff=time index in the history
1000 : !! list_dynimage(nimage)=list of dynamical images. The non-dynamical ones will not change.
1001 : !! Example : in the NEB method, or in the string method, one expect the two end images to be fixed.
1002 : !! mep_param=several parameters for Minimal Energy Path (MEP) search
1003 : !! mpi_enreg=MPI-parallelisation information
1004 : !! natom= number of atoms
1005 : !! ndynimage=number of dynamical images
1006 : !! nimage=number of images (treated by current proc)
1007 : !! nimage_tot=total number of images
1008 : !! ntimimage_stored=number of time steps stored in the history
1009 : !! pimd_param=several parameters for Path-Integral MD
1010 : !! prtvolimg=printing volume
1011 : !!
1012 : !! OUTPUT
1013 : !!
1014 : !! SIDE EFFECTS
1015 : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
1016 : !! results_img(:,:)%acell(3)
1017 : !! at input, history of the values of acell for all images
1018 : !! at output, the predicted values of acell for all images
1019 : !! results_img(:,:)%results_gs
1020 : !! at input, history of the values of energies and forces for all images
1021 : !! results_img(:,:)%rprim(3,3)
1022 : !! at input, history of the values of rprim for all images
1023 : !! at output, the predicted values of rprim for all images
1024 : !! results_img(:,:)%vel(3,natom)
1025 : !! at input, history of the values of vel for all images
1026 : !! at output, the predicted values of vel for all images
1027 : !! results_img(:,:)%vel_cell(3,3)
1028 : !! at input, history of the values of vel_cell for all images
1029 : !! at output, the predicted values of vel_cell for all images
1030 : !! results_img(:,:)%xred(3,natom)
1031 : !! at input, history of the values of xred for all images
1032 : !! at output, the predicted values of xred for all images
1033 : !!
1034 : !! SOURCE
1035 :
1036 608 : subroutine predictimg(deltae,imagealgo_str,imgmov,itimimage,itimimage_eff,list_dynimage,&
1037 : & ga_param,mep_param,mpi_enreg,m1geo_param,natom,ndynimage,nimage,nimage_tot,&
1038 304 : & ntimimage_stored,pimd_param,prtvolimg,results_img)
1039 :
1040 : use m_results_gs , only : results_gs_type
1041 : use m_predict_neb, only : predict_neb
1042 : use m_predict_steepest, only : predict_steepest
1043 : use m_predict_pimd, only : predict_pimd
1044 : use m_predict_string, only : predict_string
1045 :
1046 : !Arguments ------------------------------------
1047 : !scalars
1048 : integer,intent(in) :: imgmov,itimimage,itimimage_eff,natom,ndynimage
1049 : integer,intent(in) :: nimage,nimage_tot,ntimimage_stored,prtvolimg
1050 : character(len=60),intent(in) :: imagealgo_str
1051 : real(dp),intent(in) :: deltae
1052 : type(mep_type),intent(inout) :: mep_param
1053 : type(m1geo_type),intent(inout) :: m1geo_param
1054 : type(ga_type),intent(inout) :: ga_param
1055 : type(pimd_type),intent(inout) :: pimd_param
1056 : type(MPI_type),intent(in) :: mpi_enreg
1057 : !arrays
1058 : integer,intent(in) :: list_dynimage(ndynimage)
1059 : type(results_img_type) :: results_img(nimage,ntimimage_stored)
1060 :
1061 : !Local variables-------------------------------
1062 : !scalars
1063 : integer,save :: idum=5
1064 : logical :: is_pimd
1065 : character(len=500) :: msg
1066 : ! *************************************************************************
1067 :
1068 304 : is_pimd=(imgmov==9.or.imgmov==10.or.imgmov==13)
1069 :
1070 : !Write convergence info
1071 304 : write(msg,'(3a)') ch10,&
1072 608 : & '------------------------------------------------------------',ch10
1073 304 : if (prtvolimg<2) write(msg,'(5a)') trim(msg),' ',trim(imagealgo_str),':',ch10
1074 :
1075 : !Specific case of 4th-order RK algorithm
1076 304 : if (mep_param%mep_solver==4) then
1077 0 : if (mod(itimimage,4)==0) then
1078 0 : write(msg,'(4a)') trim(msg),&
1079 0 : & ' Fourth-order Runge-Kutta algorithm - final step',ch10
1080 0 : if (itimimage>4) write(msg,'(2a,es11.3,2a)') trim(msg),&
1081 0 : & ' Average[Abs(Etotal(t)-Etotal(t-dt))]=',deltae,' Hartree',ch10
1082 0 : write(msg,'(2a)') trim(msg),' Moving images of the cell...'
1083 : else
1084 0 : write(msg,'(2a,i1,2a)') trim(msg),&
1085 0 : & ' Fourth-order Runge-Kutta algorithm - intermediate step ',mod(itimimage,4),ch10
1086 0 : write(msg,'(2a)') trim(msg),' Computing new intermediate positions...'
1087 : end if
1088 304 : else if (is_pimd) then
1089 :
1090 : ! PIMD
1091 35 : write(msg,'(2a)') trim(msg),' Moving images of the cell...'
1092 : else
1093 :
1094 : ! Other cases
1095 502 : if (itimimage>1) write(msg,'(2a,es11.3,2a)') trim(msg),&
1096 466 : & ' Average[Abs(Etotal(t)-Etotal(t-dt))]=',deltae,' Hartree',ch10
1097 269 : write(msg,'(2a)') trim(msg),' Moving images of the cell...'
1098 :
1099 : end if
1100 :
1101 : !Write the msg
1102 : !Prevent writing if iexit==1, which at present only happens for imgmov==6 algo
1103 900 : if(imgmov/=6 .or. m1geo_param%iexit==0) call wrtout([std_out, ab_out] ,msg)
1104 :
1105 1 : select case(imgmov)
1106 :
1107 : case(0)
1108 : call predict_copy(itimimage_eff,list_dynimage,ndynimage,nimage,&
1109 1 : & ntimimage_stored,results_img)
1110 :
1111 : case(1)
1112 : call predict_steepest(itimimage,itimimage_eff,list_dynimage,mep_param,natom,ndynimage,nimage,&
1113 10 : & ntimimage_stored,results_img)
1114 :
1115 : case(2)
1116 : call predict_string(itimimage,itimimage_eff,list_dynimage,mep_param,mpi_enreg,natom,&
1117 182 : & ndynimage,nimage,nimage_tot,ntimimage_stored,results_img)
1118 :
1119 : case(4)
1120 : call predict_ga(itimimage_eff,idum,ga_param,natom,nimage,&
1121 4 : & ntimimage_stored,results_img)
1122 :
1123 : case(5)
1124 : call predict_neb(itimimage,itimimage_eff,list_dynimage,mep_param,mpi_enreg,natom,&
1125 42 : & ndynimage,nimage,nimage_tot,ntimimage_stored,results_img)
1126 :
1127 : case(6)
1128 30 : call move_1geo(itimimage_eff,m1geo_param,mpi_enreg,nimage,nimage_tot,ntimimage_stored,pimd_param,results_img)
1129 :
1130 : case(9, 10, 13)
1131 : ! Path Integral Molecular Dynamics
1132 : call predict_pimd(imgmov,itimimage,itimimage_eff,mpi_enreg,natom,nimage,nimage_tot,&
1133 304 : & ntimimage_stored,pimd_param,prtvolimg,results_img)
1134 :
1135 : case default
1136 :
1137 : end select
1138 :
1139 304 : end subroutine predictimg
1140 : !!***
1141 :
1142 : !!****f* ABINIT/predict_copy
1143 : !! NAME
1144 : !! predict_copy
1145 : !!
1146 : !! FUNCTION
1147 : !! Given the past history of images, predict the new set of images.
1148 : !! Here, simple copy of the previous image.
1149 : !!
1150 : !! INPUTS
1151 : !! itimimage_eff=time index in the history
1152 : !! list_dynimage(nimage)=list of dynamical images. The non-dynamical ones will not change.
1153 : !! Example : in the NEB of string method, one expect the two end images to be fixed.
1154 : !! ndynimage=number of dynamical images
1155 : !! nimage=number of images
1156 : !! ntimimage_stored=number of time steps stored in the history
1157 : !!
1158 : !! OUTPUT
1159 : !!
1160 : !! SIDE EFFECTS
1161 : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
1162 : !! results_img(:,:)%acell(3)
1163 : !! at input, history of the values of acell for all images
1164 : !! at output, the predicted values of acell for all images
1165 : !! results_img(:,:)%results_gs
1166 : !! at input, history of the values of energies and forces for all images
1167 : !! results_img(:,:)%rprim(3,3)
1168 : !! at input, history of the values of rprim for all images
1169 : !! at output, the predicted values of rprim for all images
1170 : !! results_img(:,:)%vel(3,natom)
1171 : !! at input, history of the values of vel for all images
1172 : !! at output, the predicted values of vel for all images
1173 : !! results_img(:,:)%vel_cell(3,3)
1174 : !! at input, history of the values of vel_cell for all images
1175 : !! at output, the predicted values of vel_cell for all images
1176 : !! results_img(:,:)%xred(3,natom)
1177 : !! at input, history of the values of xred for all images
1178 : !! at output, the predicted values of xred for all images
1179 : !!
1180 : !! SOURCE
1181 :
1182 1 : subroutine predict_copy(itimimage_eff,list_dynimage,ndynimage,nimage,&
1183 1 : & ntimimage_stored,results_img)
1184 :
1185 : !Arguments ------------------------------------
1186 : !scalars
1187 : integer,intent(in) :: itimimage_eff,ndynimage,nimage,ntimimage_stored
1188 : !arrays
1189 : integer,intent(in) :: list_dynimage(ndynimage)
1190 : type(results_img_type),intent(inout) :: results_img(nimage,ntimimage_stored)
1191 :
1192 : !Local variables-------------------------------
1193 : !scalars
1194 : integer :: idynimage,iimage,next_itimimage
1195 : ! *************************************************************************
1196 :
1197 1 : next_itimimage=itimimage_eff+1
1198 1 : if (next_itimimage>ntimimage_stored) next_itimimage=1
1199 :
1200 3 : do idynimage=1,ndynimage
1201 :
1202 2 : iimage=list_dynimage(idynimage)
1203 :
1204 8 : results_img(iimage,next_itimimage)%acell(:) =results_img(iimage,itimimage_eff)%acell(:)
1205 26 : results_img(iimage,next_itimimage)%rprim(:,:) =results_img(iimage,itimimage_eff)%rprim(:,:)
1206 18 : results_img(iimage,next_itimimage)%vel(:,:) =results_img(iimage,itimimage_eff)%vel(:,:)
1207 26 : results_img(iimage,next_itimimage)%vel_cell(:,:)=results_img(iimage,itimimage_eff)%vel_cell(:,:)
1208 19 : results_img(iimage,next_itimimage)%xred(:,:) =results_img(iimage,itimimage_eff)%xred(:,:)
1209 :
1210 : end do ! idynimage
1211 :
1212 304 : end subroutine predict_copy
1213 : !!***
1214 :
1215 : !!****f* ABINIT/move_1geo
1216 : !! NAME
1217 : !! move_1geo
1218 : !!
1219 : !! FUNCTION
1220 : !! This subroutine uses the forces, stresses and other results obtained for several images with one, common, geometry,
1221 : !! weight them to deliver averaged forces, stresses, etc, and uses these to predict the next common geometry.
1222 : !! All images must be dynamical.
1223 : !! WARNING : at present, only forces are used, to change atomic positions. No change of cell geometry.
1224 : !! Since this is not the PIMD algorithm, suppose ntimimage_stored=ntimimage, and itimimage=itimimage_eff.
1225 : !!
1226 : !! INPUTS
1227 : !! itimimage_eff=time index in the history
1228 : !! nimage=number of images
1229 : !! ntimimage_stored=number of time steps stored in the history
1230 : !! mpi_enreg=MPI-parallelisation information
1231 : !! m1geo_param=parameters for the 1geo algorithms
1232 : !! pimd_param=datastructure that contains all the parameters necessary to Path-Integral MD
1233 : !!
1234 : !! OUTPUT
1235 : !!
1236 : !! SIDE EFFECTS
1237 : !! results_img(ntimimage_stored,nimage)=datastructure that holds the history of previous computations.
1238 : !! results_img(:,:)%acell(3)
1239 : !! at input, history of the values of acell for all images
1240 : !! at output, the predicted values of acell for all images
1241 : !! results_img(:,:)%results_gs
1242 : !! at input, history of the values of energies and forces for all images
1243 : !! results_img(:,:)%rprim(3,3)
1244 : !! at input, history of the values of rprim for all images
1245 : !! at output, the predicted values of rprim for all images
1246 : !! results_img(:,:)%vel(3,natom)
1247 : !! at input, history of the values of vel for all images
1248 : !! at output, the predicted values of vel for all images
1249 : !! results_img(:,:)%vel_cell(3,3)
1250 : !! at input, history of the values of vel_cell for all images
1251 : !! at output, the predicted values of vel_cell for all images
1252 : !! results_img(:,:)%xred(3,natom)
1253 : !! at input, history of the values of xred for all images
1254 : !! at output, the predicted values of xred for all images
1255 : !!
1256 : !! SOURCE
1257 :
1258 30 : subroutine move_1geo(itimimage_eff,m1geo_param,mpi_enreg,nimage,nimage_tot,ntimimage_stored,pimd_param,results_img)
1259 :
1260 : !Arguments ------------------------------------
1261 : !scalars
1262 : integer,intent(in) :: itimimage_eff,nimage,nimage_tot,ntimimage_stored
1263 : type(MPI_type),intent(in) :: mpi_enreg
1264 : type(m1geo_type),intent(inout) :: m1geo_param
1265 : type(pimd_type),intent(in) :: pimd_param
1266 : !arrays
1267 : type(results_img_type),target,intent(inout) :: results_img(nimage,ntimimage_stored)
1268 :
1269 : !Local variables-------------------------------
1270 : !scalars
1271 : integer :: ihist,iimage,natom,next_itimimage,nspden,nsppol
1272 : !integer :: iatom
1273 : real(dp) :: deltae,diffor,etotal,entropy,fermie,res2,residm
1274 : logical :: test_img
1275 30 : type(results_gs_type) :: results_gs_lincomb
1276 : !arrays
1277 : real(dp) :: acell(3),rprim(3,3),rprimd(3,3),strten(6),vel_cell(3,3)
1278 30 : real(dp),allocatable :: fcart(:,:),vel(:,:),xred(:,:)
1279 : logical :: DEBUG=.FALSE.
1280 30 : type(results_img_type),pointer :: resimg_all(:)
1281 : ! *************************************************************************
1282 :
1283 30 : natom=m1geo_param%ab_mover%natom
1284 30 : ihist=m1geo_param%hist_1geo%ihist
1285 :
1286 90 : ABI_MALLOC(fcart,(3,natom))
1287 60 : ABI_MALLOC(vel,(3,natom))
1288 60 : ABI_MALLOC(xred,(3,natom))
1289 :
1290 : !Of course, assume that the geometry parameters are the same for all images, so take them from the first one.
1291 342 : xred(:,:) =results_img(1,itimimage_eff)%xred(:,:)
1292 120 : acell(:) =results_img(1,itimimage_eff)%acell(:)
1293 390 : rprim(:,:) =results_img(1,itimimage_eff)%rprim(:,:)
1294 342 : vel(:,:) =results_img(1,itimimage_eff)%vel(:,:)
1295 390 : vel_cell(:,:)=results_img(1,itimimage_eff)%vel_cell(:,:)
1296 :
1297 30 : call mkrdim(acell,rprim,rprimd)
1298 :
1299 : !Fill history with the values of xred, acell and rprimd
1300 30 : call var2hist(acell,m1geo_param%hist_1geo,natom,rprimd,xred,DEBUG)
1301 :
1302 : !Fill history with velocities and ionic kinetic energy
1303 30 : call vel2hist(m1geo_param%ab_mover%amass,m1geo_param%hist_1geo,vel,vel_cell)
1304 30 : m1geo_param%hist_1geo%time(ihist)=zero
1305 :
1306 : !In case of image parallelism, collect results accross processors
1307 30 : test_img=(nimage_tot/=1.and.mpi_enreg%paral_img==1)
1308 : if (test_img) then
1309 0 : ABI_MALLOC(resimg_all,(nimage_tot))
1310 : call gather_results_img(mpi_enreg,results_img(1:nimage,itimimage_eff),resimg_all,&
1311 0 : & allgather=.true.,only_one_per_img=.false.)
1312 : else
1313 30 : resimg_all => results_img(:,itimimage_eff)
1314 : end if
1315 :
1316 : !Compute energy, entropy, fermie, forces and stresses for the 1geo : take the weighted average.
1317 : !Compute maximum of deltae,diffor,res2,residm
1318 30 : etotal=zero
1319 30 : entropy=zero
1320 30 : fermie=zero
1321 342 : fcart(:,:)=zero
1322 30 : strten(:)=zero
1323 30 : deltae=zero
1324 30 : diffor=zero
1325 30 : res2=zero
1326 30 : residm=zero
1327 :
1328 88 : do iimage=1,nimage_tot
1329 58 : etotal=etotal+resimg_all(iimage)%results_gs%etotal*m1geo_param%mixesimgf(iimage)
1330 58 : entropy=entropy+resimg_all(iimage)%results_gs%entropy*m1geo_param%mixesimgf(iimage)
1331 58 : fermie=fermie+resimg_all(iimage)%results_gs%fermie*m1geo_param%mixesimgf(iimage)
1332 738 : fcart(:,:)=fcart(:,:)+resimg_all(iimage)%results_gs%fcart(:,:)*m1geo_param%mixesimgf(iimage)
1333 406 : strten(:) =strten(:) +resimg_all(iimage)%results_gs%strten(:)*m1geo_param%mixesimgf(iimage)
1334 58 : if( deltae<resimg_all(iimage)%results_gs%deltae ) deltae=resimg_all(iimage)%results_gs%deltae
1335 58 : if( diffor<resimg_all(iimage)%results_gs%diffor ) diffor=resimg_all(iimage)%results_gs%diffor
1336 58 : if( res2<resimg_all(iimage)%results_gs%res2 ) res2=resimg_all(iimage)%results_gs%res2
1337 88 : if( residm<resimg_all(iimage)%results_gs%residm ) residm=resimg_all(iimage)%results_gs%residm
1338 : enddo
1339 :
1340 : !Set up a results_gs datastructure with the linear combination of images
1341 30 : nspden=resimg_all(1)%results_gs%nspden
1342 30 : nsppol=resimg_all(1)%results_gs%nsppol
1343 30 : call results_gs_lincomb%init(natom,nspden,nsppol)
1344 30 : call resimg_all(1)%results_gs%copy(results_gs_lincomb)
1345 30 : results_gs_lincomb%etotal=etotal
1346 30 : results_gs_lincomb%entropy=entropy
1347 30 : results_gs_lincomb%fermie=fermie
1348 372 : results_gs_lincomb%fcart=fcart
1349 210 : results_gs_lincomb%strten=strten
1350 30 : results_gs_lincomb%deltae=deltae
1351 30 : results_gs_lincomb%diffor=diffor
1352 30 : results_gs_lincomb%res2=res2
1353 30 : results_gs_lincomb%residm=residm
1354 :
1355 : !From now on, all procs contain the same information about the geometry, etotal, forces, stress, etc.
1356 : !Nothing more needs to be transmitted, and resimg_all is not needed anymore.
1357 30 : if (test_img) then
1358 0 : call destroy_results_img(resimg_all)
1359 0 : ABI_FREE(resimg_all)
1360 : end if
1361 30 : nullify(resimg_all)
1362 :
1363 : !Echo result_gs_lincomb
1364 30 : call results_gs_lincomb%yaml_write(ab_out, info="Linear combination of ground state results")
1365 :
1366 : !Destroy result_gs_lincomb
1367 30 : call results_gs_lincomb%free()
1368 :
1369 : !Store fcart and strten in hist_1geo
1370 342 : m1geo_param%hist_1geo%fcart(:,:,ihist)=fcart(:,:)
1371 210 : m1geo_param%hist_1geo%strten(:,ihist) =strten(:)
1372 :
1373 : !Store them in ab_xfh
1374 : !THIS IS TO BE DONE !
1375 :
1376 : !Compute new atomic positions and cell characteristics in the single geometry
1377 : call precpred_1geo(m1geo_param%ab_mover,&
1378 : & m1geo_param%ab_xfh_1geo,&
1379 : & m1geo_param%ab_mover%amu_curr,&
1380 : & m1geo_param%deloc,&
1381 : & m1geo_param%dt_chkdilatmx,&
1382 : & mpi_enreg%comm_cell,&
1383 : & m1geo_param%dilatmx,&
1384 : & m1geo_param%filnam_ds4,&
1385 : & m1geo_param%hist_1geo,&
1386 : & m1geo_param%hmctt,&
1387 : & m1geo_param%icycle,&
1388 : & m1geo_param%iexit,&
1389 : & itimimage_eff,&
1390 : & m1geo_param%mttk_vars,&
1391 : & m1geo_param%nctime,&
1392 : & m1geo_param%ncycle,&
1393 : & m1geo_param%nerr_dilatmx,&
1394 : & m1geo_param%npsp,&
1395 : & m1geo_param%ntime,&
1396 : & pimd_param,&
1397 : & m1geo_param%rprimd_orig,&
1398 : & m1geo_param%skipcycle,&
1399 30 : & m1geo_param%usewvl)
1400 :
1401 : !Retrieve the new positions, cell parameters [and velocities ?!]
1402 30 : call hist2var(acell,m1geo_param%hist_1geo,natom,rprimd,xred,DEBUG)
1403 :
1404 : !Store acell, rprim, xred and vel for the new iteration if relevant
1405 30 : if(m1geo_param%iexit==0)then
1406 24 : next_itimimage=itimimage_eff+1
1407 24 : if (next_itimimage>ntimimage_stored)then
1408 0 : ABI_ERROR('next_itimimage>ntimimage_stored')
1409 : endif
1410 70 : do iimage=1,nimage
1411 558 : results_img(iimage,next_itimimage)%xred(:,:) =xred(:,:)
1412 184 : results_img(iimage,next_itimimage)%acell(:) =acell(:)
1413 622 : results_img(iimage,next_itimimage)%rprim(:,:) =rprim(:,:)
1414 : ! WARNING : Should also store vel and vel_cell of course ...
1415 : ! results_img(iimage,next_itimimage)%vel(:,:) =vel(:,:)
1416 : ! results_img(iimage,next_itimimage)%vel_cell(:,:)=vel_cell(:,:)
1417 : end do
1418 : endif
1419 30 : ABI_FREE(fcart)
1420 30 : ABI_FREE(vel)
1421 30 : ABI_FREE(xred)
1422 :
1423 60 : end subroutine move_1geo
1424 : !!***
1425 :
1426 : end module m_gstateimg
1427 : !!***
|
