Line data Source code
1 : !!****m* ABINIT/m_mover
2 : !! NAME
3 : !! m_mover
4 : !!
5 : !! FUNCTION
6 : !! Move ion or change acell according to forces and stresses
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 1998-2026 ABINIT group (DCA, XG, GMR, SE, FLambert,MT)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_mover
23 :
24 : use defs_basis
25 : use m_abicore
26 : use m_errors
27 : use m_profiling_abi
28 : use m_abimover
29 : use m_abihist
30 : use m_dtset
31 : use m_pimd
32 : use m_xmpi
33 : use m_nctk
34 : use m_dtfil
35 : use m_yaml
36 : use netcdf
37 : #if defined HAVE_LOTF
38 : use lotfpath
39 : use m_pred_lotf
40 : #endif
41 :
42 : use defs_abitypes, only : MPI_type
43 : use m_fstrings, only : strcat, sjoin, indent, itoa
44 : use m_matrix, only : matr3inv
45 : use m_symtk, only : symmetrize_xred
46 : use m_geometry, only : fcart2gred, chkdilatmx, xred2xcart, metric
47 : use m_time, only : abi_wtime, sec2str
48 : use m_exit, only : get_start_time, have_timelimit_in, get_timelimit, enable_timelimit_in
49 : use m_electronpositron, only : electronpositron_type
50 : use m_scfcv, only : scfcv_t, scfcv_run
51 : use m_effective_potential,only : effective_potential_type, effective_potential_evaluate
52 : use m_initylmg, only : initylmg
53 : use m_kg, only : getcut, getph
54 : use m_xfpack, only : xfh_update
55 : use m_mkrho, only : initro
56 : use m_pawfgr, only : pawfgr_type, pawfgr_init, pawfgr_destroy
57 : use m_precpred_1geo, only : precpred_1geo
58 : use m_pred_simple, only : prec_simple
59 : !use m_generate_training_set, only : generate_training_set
60 : use m_wvl_wfsinp, only : wvl_wfsinp_reformat
61 : use m_wvl_rho, only : wvl_mkrho
62 : use m_effective_potential_file, only : effective_potential_file_mapHistToRef
63 : #if defined DEV_MS_SCALEUP
64 : use scup_global, only : global_set_parent_iter,global_set_print_parameters
65 : #endif
66 : use m_scup_dataset
67 : use m_multibinit_dataset
68 :
69 : implicit none
70 :
71 : private
72 : !!***
73 :
74 : public :: mover
75 : !!***
76 :
77 : contains
78 : !!***
79 :
80 : !!****f* ABINIT/mover
81 : !! NAME
82 : !! mover
83 : !!
84 : !! FUNCTION
85 : !! Move ion or change acell acording to forces and stresses
86 : !!
87 : !! INPUTS
88 : !! amu_curr(ntypat)=mass of each atom for the current image
89 : !! dtfil <type(datafiles_type)>=variables related to files
90 : !! dtset <type(dataset_type)>=all input variables for this dataset
91 : !! | mband=maximum number of bands
92 : !! | mgfft=maximum size of 1D FFTs
93 : !! | mkmem =number of k points treated by this node
94 : !! | angular momentum for nonlocal pseudopotential
95 : !! | mpw=maximum dimensioned size of npw.
96 : !! | natom=number of atoms in unit cell
97 : !! | except on first call (hartree/bohr); updated on output
98 : !! | nfft=(effective) number of FFT grid points (for this processor)
99 : !! | for the "coarse" grid (see NOTES below)
100 : !! | nkpt=number of k points.
101 : !! | nspden=number of spin-density components
102 : !! | nsppol=1 for unpolarized, 2 for spin-polarized
103 : !! | nsym=number of symmetry elements in space group
104 : !! itimimage_gstate= [optional] counter for the itimimage loop, in the calling routine.
105 : !! mcg=size of wave-functions array (cg) =mpw*nspinor*mband*mkmem*nsppol
106 : !! mpi_enreg=information about MPI parallelization
107 : !! nfftf=(effective) number of FFT grid points (for this processor)
108 : !! for the "fine" grid (see NOTES below)
109 : !! npwarr(nkpt)=number of planewaves in basis and boundary at this k point.
110 : !! nattyp(ntypat)= # atoms of each type.
111 : !! paw_dmft <type(paw_dmft_type)>= paw+dmft related data
112 : !! psps <type(pseudopotential_type)>=variables related to pseudopotentials
113 : !! | mpsang= 1+maximum angular momentum for nonlocal pseudopotentials
114 : !! rprimd(3,3)=dimensional primitive translations (bohr)
115 : !! scup_dtset <type(scup_dtset_type) = derived datatype holding all options
116 : !! for the evaluation of an effective electronic model using SCALE UP
117 : !!
118 : !! OUTPUT
119 : !! results_gs <type(results_gs_type)>=results (energy and its components,
120 : !! forces and its components, the stress tensor) of a ground-state computation
121 : !! eigen(mband*nkpt*nsppol)=array for holding eigenvalues (hartree)
122 : !! resid(mband*nkpt*nsppol)=residuals for each band over all k points.
123 : !!
124 : !! SIDE EFFECTS
125 : !! Rest of i/o is related to lda
126 : !! acell(3)=length scales of primitive translations (bohr)
127 : !! cg(2,mcg)=array for planewave coefficients of wavefunctions.
128 : !! electronpositron <type(electronpositron_type)>=quantities for the electron-positron annihilation
129 : !! initialized= if 0 the initialisation of the gstate run is not yet finished
130 : !! irrzon(nfft**(1-1/nsym),2,(nspden/nsppol)-3*(nspden/4))=irreducible zone data
131 : !! occ(mband*nkpt*nsppol=occupation number for each band (usually 2) at each k point.
132 : !! rhog(2,nfftf)=array for Fourier transform of electron density
133 : !! rhor(nfftf,nspden)=array for electron density in electrons/bohr**3.
134 : !! scf_history <type(scf_history_type)>=arrays obtained from previous SCF cycles
135 : !! symrec(3,3,nsym)=symmetry operations in reciprocal space
136 : !! taug(2,nfftf*dtset%usekden)=array for Fourier transform of kinetic energy density
137 : !! taur(nfftf,nspden*dtset%usekden)=array for kinetic energy density
138 : !! vel(3,natom)=old value of velocity; updated on output
139 : !! vel_cell(3,3)=old value of cell parameters velocity; updated on output
140 : !! xred(3,natom)=reduced dimensionless atomic coordinates; updated on output
141 : !! xred_old(3,natom)=work space for old xred
142 : !! eff_pot<type(effective_potential_type)> = optional,effective_potential datatype
143 : !! verbose = optional, default is true, flag to disable the verbose mode
144 : !! write_HIST = optional, default is true, flag to disble the write of the HIST file
145 : !!
146 : !! NOTES
147 : !! This subroutine uses the arguments natom, xred, vel, amu_curr,
148 : !! vis, and dtion (the last two contained in dtset) to make
149 : !! molecular dynamics updates. The rest of the lengthy
150 : !! argument list supports the underlying lda computation
151 : !! of forces, returned from subroutine scfcv
152 : !!
153 : !! USE OF FFT GRIDS:
154 : !! =================
155 : !! In case of PAW:
156 : !! ---------------
157 : !! Two FFT grids are used:
158 : !! - A "coarse" FFT grid (defined by ecut)
159 : !! for the application of the Hamiltonian on the plane waves basis.
160 : !! It is defined by nfft, ngfft, mgfft, ...
161 : !! Hamiltonian, wave-functions, density related to WFs (rhor here), ...
162 : !! are expressed on this grid.
163 : !! - A "fine" FFT grid (defined) by ecutdg)
164 : !! for the computation of the density inside PAW spheres.
165 : !! It is defined by nfftf, ngfftf, mgfftf, ...
166 : !! Total density, potentials, ...
167 : !! are expressed on this grid.
168 : !! In case of norm-conserving:
169 : !! ---------------------------
170 : !! - Only the usual FFT grid (defined by ecut) is used.
171 : !! It is defined by nfft, ngfft, mgfft, ...
172 : !! For compatibility reasons, (nfftf,ngfftf,mgfftf)
173 : !! are set equal to (nfft,ngfft,mgfft) in that case.
174 : !!
175 : !! SOURCE
176 :
177 393 : subroutine mover(scfcv_args,ab_xfh,acell,amu_curr,dtfil,&
178 393 : & electronpositron,rhog,rhor,rprimd,vel,vel_cell,xred,xred_old,&
179 : & effective_potential,filename_ddb,itimimage_gstate,verbose,verbose_cycle,writeHIST,scup_dtset,sc_size,multibinit_dtset)
180 :
181 : !Arguments ------------------------------------
182 : !scalars
183 : integer, intent(in), optional :: itimimage_gstate
184 : type(scfcv_t),intent(inout) :: scfcv_args
185 : type(datafiles_type),intent(inout),target :: dtfil
186 : type(electronpositron_type),pointer :: electronpositron
187 : type(ab_xfh_type),intent(inout) :: ab_xfh
188 : type(effective_potential_type),optional,intent(inout) :: effective_potential
189 : logical,optional,intent(in) :: verbose
190 : logical,optional,intent(in) :: verbose_cycle
191 : logical,optional,intent(in) :: writeHIST
192 : character(len=fnlen),optional,intent(in) :: filename_ddb
193 : !arrays
194 : real(dp),intent(inout) :: acell(3)
195 : real(dp), intent(in),target :: amu_curr(:) !(scfcv%dtset%ntypat)
196 : real(dp), pointer :: rhog(:,:),rhor(:,:)
197 : real(dp), intent(inout) :: xred(3,scfcv_args%dtset%natom),xred_old(3,scfcv_args%dtset%natom)
198 : real(dp), intent(inout) :: vel(3,scfcv_args%dtset%natom),vel_cell(3,3),rprimd(3,3)
199 : type(scup_dtset_type),optional, intent(inout) :: scup_dtset
200 : type(multibinit_dtset_type),optional, intent(inout) ::multibinit_dtset
201 : integer,optional,intent(in) :: sc_size(3)
202 :
203 : !Local variables-------------------------------
204 : !scalars
205 : integer,parameter :: level=102,master=0
206 393 : type(abihist) :: hist,hist_prev
207 : type(abimover) :: ab_mover
208 : type(abimover_specs) :: specs
209 393 : type(abiforstr) :: preconforstr ! Preconditioned forces and stress
210 393 : type(delocint) :: deloc
211 : type(pawfgr_type) :: pawfgr
212 393 : type(mttk_type) :: mttk_vars
213 393 : type(pimd_type) :: pimd_param
214 : integer :: itime,icycle,itime_hist,iexit=0,ifirst,ihist_prev,ihist_prev2,timelimit_exit,ncycle,nhisttot,kk,jj,me
215 : integer :: nmpi,ntime,option,comm,mgfftf,nfftf
216 : integer :: nerr_dilatmx,my_quit,ierr,quitsum_request
217 : integer ABI_ASYNC :: quitsum_async
218 : character(len=500) :: msg
219 : !character(len=500) :: dilatmx_errmsg
220 : character(len=8) :: stat4xml
221 : character(len=35) :: fmt
222 : character(len=fnlen) :: filename,fname_ddb,name_file
223 : character(len=500) :: MY_NAME = "mover"
224 : real(dp) :: gr_avg,ecut_eff,ecutdg_eff,ucvol,boxcut,gsqcut_eff
225 : logical :: DEBUG=.FALSE., effective_potential_verbose, need_verbose=.TRUE.,need_verbose_cycle=.TRUE.,need_writeHIST=.TRUE.
226 : logical :: need_scfcv_cycle = .TRUE., need_elec_eval = .FALSE.
227 : logical :: changed,useprtxfase
228 : logical :: skipcycle,force_hist_copy=.FALSE.
229 : integer :: minIndex,ii,similar,conv_retcode
230 : integer :: iapp
231 : logical :: file_exists
232 : logical :: re_init_rho
233 : #ifdef FC_NVHPC
234 : logical :: wrong=.false. !Silly trick to prevent NVHPC optimization issue
235 : #endif
236 : real(dp) :: minE,wtime_step,now,prev
237 : !arrays
238 : integer :: itimes(2),ngfft(18),ngfftf(18)
239 : real(dp) :: gprimd(3,3),rprim(3,3),rprimd_prev(3,3),gmet(3,3),rmet(3,3)
240 393 : real(dp),allocatable :: gred_corrected(:,:),xred_prev(:,:),ph1df(:,:),tmp(:,:),tmp_1d(:)
241 : real(dp) :: k0(3)
242 : ! ***************************************************************
243 393 : need_verbose=.TRUE.
244 393 : if(present(verbose)) need_verbose = verbose
245 :
246 393 : need_verbose_cycle=.TRUE.
247 393 : if(present(verbose_cycle)) need_verbose_cycle = verbose_cycle
248 :
249 393 : need_writeHIST=.TRUE.
250 393 : if(present(writeHIST)) need_writeHIST = writeHIST
251 :
252 : ! enable time limit handler if not done in callers.
253 393 : if (enable_timelimit_in(MY_NAME) == MY_NAME) then
254 393 : if (need_verbose) then
255 0 : write(std_out,*)"Enabling timelimit check in function: ",trim(MY_NAME)," with timelimit: ",trim(sec2str(get_timelimit()))
256 : end if
257 : end if
258 :
259 393 : re_init_rho = .FALSE.
260 :
261 : !Table of contents
262 : !(=>) Refers to an important call (scfcv,pred_*)
263 : !
264 : !01. Initialization of indexes and allocations of arrays
265 : !02. Particularities of each predictor
266 : !03. Set the number of iterations ntime
267 : !04. Try to read history of previous calculations
268 : !05. Allocate the hist structure
269 : !06. First output before any itime or icycle
270 : !07. Fill the history of the first SCFCV
271 : !08. Loop for itime (From 1 to ntime)
272 : !09. Loop for icycle (From 1 to ncycle)
273 : !10. Output for each icycle (and itime)
274 : !11. Symmetrize atomic coordinates over space group elements
275 : !12. => Call to SCFCV routine and fill history with forces
276 : !13. Write the history into the _HIST file
277 : !14. Output after SCFCV
278 : !15. => Test Convergence of forces and stresses
279 : !16. => Precondition forces, stress and energy
280 : !17. => Call to each predictor
281 : !18. Use the history to extract the new values
282 : !19. End loop icycle
283 : !20. End loop itime
284 : !21. Set the final values of xred
285 : !22. XML Output at the end
286 : !23. Deallocate hist and ab_mover datatypes
287 : !
288 393 : call abimover_ini(ab_mover,amu_curr,dtfil,scfcv_args%dtset,specs)
289 :
290 393 : if (ab_mover%ionmov==10 .or. ab_mover%ionmov==11) call delocint_ini(deloc)
291 :
292 393 : if (ab_mover%ionmov==13 .or. ab_mover%ionmov==25)then
293 31 : call mttk_ini(mttk_vars,ab_mover%nnos)
294 : end if
295 :
296 : !###########################################################
297 : !### 03. Set the number of iterations ntime
298 : !### By default ntime==1 but if the user enters a lower
299 : !### value, mover will execute at least one iteration
300 :
301 393 : if (scfcv_args%dtset%ntime<1)then
302 0 : ntime=1
303 : else
304 393 : ntime=scfcv_args%dtset%ntime
305 : end if
306 :
307 : !###########################################################
308 : !### 04. Try to read history of previous calculations
309 : !### It requires NetCDF library
310 :
311 : !Init MPI data
312 393 : comm=scfcv_args%mpi_enreg%comm_cell
313 393 : me=xmpi_comm_rank(comm)
314 :
315 :
316 393 : filename=trim(ab_mover%filnam_ds(4))//'_HIST.nc'
317 :
318 :
319 393 : if (ab_mover%restartxf<0)then
320 : ! Read history from file (and broadcast if MPI)
321 15 : if (me==master) then
322 15 : call read_md_hist(filename,hist_prev,specs%isVused,specs%isARused,ab_mover%restartxf==-3)
323 : endif
324 :
325 :
326 15 : call abihist_bcast(hist_prev,master,comm)
327 :
328 :
329 : ! If restartxf specifies to reconstruct the history
330 15 : if (hist_prev%mxhist>0.and.ab_mover%restartxf==-1)then
331 1 : ntime=ntime+hist_prev%mxhist
332 : end if
333 :
334 : ! If non deterministic algorithm is used, forcing reading of input hist file
335 15 : if (ab_mover%ionmov==16) then
336 0 : force_hist_copy=.TRUE.
337 : end if
338 :
339 : ! If restartxf specifies to start from the lowest energy
340 15 : if (hist_prev%mxhist>0.and.ab_mover%restartxf==-2)then
341 0 : minE=hist_prev%etot(1)
342 0 : minIndex=1
343 0 : do ii=1,hist_prev%mxhist
344 0 : if(need_verbose) write(std_out,*) 'Iteration:',ii,' Total Energy:',hist_prev%etot(ii)
345 0 : if (minE>hist_prev%etot(ii))then
346 0 : minE=hist_prev%etot(ii)
347 0 : minIndex=ii
348 : end if
349 : end do
350 0 : if(need_verbose)write(std_out,*) 'The lowest energy occurs at iteration:',minIndex,'etotal=',minE
351 0 : acell(:) =hist_prev%acell(:,minIndex)
352 0 : rprimd(:,:)=hist_prev%rprimd(:,:,minIndex)
353 0 : xred(:,:) =hist_prev%xred(:,:,minIndex)
354 0 : vel(:, :) = hist_prev%vel(:, :, minIndex)
355 0 : call abihist_free(hist_prev)
356 : end if
357 : ! If restarxf specifies to start to the last iteration
358 15 : if (hist_prev%mxhist>0.and.ab_mover%restartxf==-3)then
359 3 : if(present(effective_potential))then
360 3 : call effective_potential_file_mapHistToRef(effective_potential,hist_prev,comm,scfcv_args%dtset%iatfix,need_verbose,sc_size) ! Map Hist to Ref to order atoms
361 : !xred(:,:) = hist_prev%xred(:,:,1) ! Fill xred with new ordering
362 3 : hist%ihist = 1
363 : end if
364 12 : acell(:) =hist_prev%acell(:,hist_prev%mxhist)
365 39 : rprimd(:,:)=hist_prev%rprimd(:,:,hist_prev%mxhist)
366 323 : xred(:,:) =hist_prev%xred(:,:,hist_prev%mxhist)
367 323 : vel(:, :) = hist_prev%vel(:, :, hist_prev%mxhist)
368 3 : call abihist_free(hist_prev)
369 : end if
370 :
371 : end if !if (ab_mover%restartxf<=0)
372 :
373 : !###########################################################
374 : !### 05. Allocate the hist structure
375 :
376 393 : iexit=0; timelimit_exit=0
377 393 : ncycle=specs%ncycle
378 :
379 393 : if(ab_mover%ionmov==25.and.scfcv_args%dtset%hmctt>=0)then
380 2 : ncycle=scfcv_args%dtset%hmctt
381 2 : if(scfcv_args%dtset%hmcsst>0.and.ab_mover%optcell/=0)then
382 0 : ncycle=ncycle+scfcv_args%dtset%hmcsst
383 : endif
384 : endif
385 :
386 393 : nhisttot=ncycle*ntime;if (scfcv_args%dtset%nctime>0) nhisttot=nhisttot+1
387 : !AM_2017 New version of the hist, we just store the needed history step not all of them...
388 393 : if(specs%nhist/=-1)then
389 393 : nhisttot = specs%nhist! We don't need to store all the history
390 : endif
391 :
392 393 : call abihist_init(hist,ab_mover%natom,nhisttot,specs%isVused,specs%isARused)
393 393 : call abiforstr_ini(preconforstr,ab_mover%natom)
394 :
395 : !###########################################################
396 : !### 06. First output before any itime or icycle iteration
397 :
398 : !If effective potential is present forces will be compute with it
399 393 : if (present(effective_potential))then
400 44 : need_scfcv_cycle = .FALSE.
401 44 : if(need_verbose)then
402 : write(msg,'(2a,i2,5a,80a)')&
403 44 : & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',trim(specs%method),' with effective potential',&
404 3608 : & ch10,('=',kk=1,80)
405 132 : call wrtout([std_out, ab_out], msg)
406 : end if
407 44 : need_elec_eval = .FALSE.
408 44 : if(present(scup_dtset))then
409 44 : need_elec_eval = scup_dtset%scup_elec_model
410 : endif
411 : else
412 349 : if(need_verbose)then
413 : write(msg,'(a,a,i2,a,a,a,80a)')&
414 349 : & ch10,'=== [ionmov=',ab_mover%ionmov,'] ',trim(specs%method),&
415 28618 : & ch10,('=',kk=1,80)
416 1047 : call wrtout([std_out, ab_out], msg)
417 : end if
418 : end if
419 :
420 : !Format for printing on each cycle
421 : write(fmt,'(a6,i2,a4,i2,a4,i2,a4,i2,a9)')&
422 393 : '(a,a,i',int(log10(real(ntime))+1),&
423 393 : ',a,i',int(log10(real(ntime))+1),&
424 393 : ',a,i',int(log10(real(ncycle))+1),&
425 393 : ',a,i',int(log10(real(ncycle))+1),&
426 786 : ',a,a,80a)'
427 :
428 : !###########################################################
429 : !### 07. Fill the history of the first SCFCV
430 :
431 393 : if (ab_mover%ionmov==26)then
432 :
433 : !Tdep call need to merge with adewandre branch
434 393 : else if (ab_mover%ionmov==27)then
435 0 : if(present(filename_ddb))then
436 0 : fname_ddb = trim(filename_ddb)
437 : else
438 0 : fname_ddb = trim(ab_mover%filnam_ds(3))//'_DDB'
439 : end if
440 0 : INQUIRE(FILE=filename, EXIST=file_exists)
441 :
442 0 : ABI_ERROR("This section has been disabled, ph_freez_disp is not defined in main ABINIT")
443 :
444 : ! XG 20200322 : The input variables ph_freez_disp are not documented neither tested, so they
445 : ! have been removed from the allowed list in the parser. Also, you should not be here !
446 : ! call generate_training_set(acell,ab_mover%ph_freez_disp_addStrain==1,ab_mover%ph_freez_disp_ampl,&
447 : !& fname_ddb,hist,ab_mover%natom,ab_mover%ph_freez_disp_nampl,ntime,&
448 : !& ab_mover%ph_ngqpt,ab_mover%ph_nqshift,ab_mover%ph_freez_disp_option,&
449 : !& ab_mover%ph_qshift,scfcv_args%dtset%supercell_latt,&
450 : !& rprimd,ab_mover%mdtemp(2),xred,comm,DEBUG)
451 :
452 :
453 : !Fill history with the values of xred, acell and rprimd of the first configuration
454 0 : acell(:) =hist%acell(:,1)
455 0 : rprimd(:,:)=hist%rprimd(:,:,1)
456 0 : xred(:,:) =hist%xred(:,:,1)
457 :
458 : else
459 :
460 : ! Fill history with the values of xred, acell and rprimd
461 393 : call var2hist(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
462 :
463 :
464 : ! Fill velocities and ionic kinetic energy
465 393 : call vel2hist(ab_mover%amass,hist,vel,vel_cell)
466 393 : hist%time(hist%ihist)=zero
467 :
468 : end if
469 :
470 : !Decide if prtxfase will be called
471 393 : useprtxfase=.FALSE.
472 5276 : do ii=1,ab_mover%natom
473 5276 : if (ab_mover%prtatlist(ii)/=0)then
474 : useprtxfase=.TRUE.; exit
475 : end if
476 : end do
477 :
478 : !At beginning no error
479 393 : nerr_dilatmx = 0
480 : !Copy the number of degrees of freedom in hist structure
481 393 : hist%ndof=ab_mover%ndof
482 :
483 1179 : ABI_MALLOC(xred_prev,(3,scfcv_args%dtset%natom))
484 :
485 393 : if (present(effective_potential) .and. need_verbose .and. .not.need_verbose_cycle .and. ab_mover%ionmov==25) then
486 162 : write(msg,'(3a,80a)') ch10,'--- Initial state',ch10,('-',kk=1,80)
487 6 : call wrtout([std_out, ab_out], msg)
488 2 : if(present(multibinit_dtset))then
489 : call effective_potential_evaluate( &
490 : & effective_potential,scfcv_args%results_gs%etotal,scfcv_args%results_gs%fcart,scfcv_args%results_gs%gred,&
491 : & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=need_verbose,&
492 : & elec_eval=need_elec_eval,efield_type=multibinit_dtset%efield_type,efield=multibinit_dtset%efield,&
493 : & efield_lambda=multibinit_dtset%efield_lambda,nefield=multibinit_dtset%nefield,&
494 : & efield_period=multibinit_dtset%efield_period,efield_phase=multibinit_dtset%efield_phase,&
495 : & efield_gmean=multibinit_dtset%efield_gmean,efield_gvel=multibinit_dtset%efield_gvel,efield_sigma=multibinit_dtset%efield_sigma,&
496 2 : & efield_background=multibinit_dtset%efield_background,time=zero)
497 : else
498 : call effective_potential_evaluate( &
499 : & effective_potential,scfcv_args%results_gs%etotal,scfcv_args%results_gs%fcart,scfcv_args%results_gs%gred,&
500 : & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=need_verbose,&
501 0 : & elec_eval=need_elec_eval,time=zero)
502 : end if
503 : end if
504 :
505 : !###########################################################
506 : !### 08. Loop for itime (From 1 to ntime)
507 393 : quitsum_request = xmpi_request_null
508 :
509 3655 : do itime=1,ntime
510 :
511 3655 : call yaml_iterstart("itime", itime, dev_null, scfcv_args%dtset%use_yaml)
512 :
513 : ! Handle time limit condition.
514 3655 : if (itime == 1) prev = abi_wtime()
515 3655 : if (itime > 1) then
516 3262 : now = abi_wtime()
517 3262 : wtime_step = now - prev
518 3262 : prev = now
519 3262 : write(msg,*)sjoin("{mover_itime:", itoa(itime - 1), ", wall_time: '", sec2str(wtime_step), "'} <<< TIME")
520 3262 : if(need_verbose)call wrtout(std_out, msg)
521 3262 : if (have_timelimit_in(MY_NAME)) then
522 0 : if (itime > 2) then
523 0 : call xmpi_wait(quitsum_request,ierr)
524 0 : if (quitsum_async > 0) then
525 0 : write(msg,"(3a)")"Approaching time limit ",trim(sec2str(get_timelimit())), ". Will exit itime loop in mover."
526 0 : if(need_verbose) ABI_COMMENT(msg)
527 0 : if(need_verbose) call wrtout(ab_out, msg)
528 : timelimit_exit = 1
529 : exit
530 : end if
531 : end if
532 :
533 0 : my_quit = 0; if (now - get_start_time() + 2.15 * wtime_step > get_timelimit()) my_quit = 1
534 0 : call xmpi_isum(my_quit,quitsum_async,comm,quitsum_request,ierr)
535 : end if
536 : end if
537 :
538 3655 : skipcycle=.FALSE.
539 : #if defined HAVE_LOTF
540 : if(ab_mover%ionmov==23 .and. .not. lotf_extrapolation(itime)) skipcycle=.True.
541 : #endif
542 :
543 : ! If RMM-DIIS is used, decrease the number of NSCF steps done with wfoptalg before activating RMM-DIIS.
544 : ! In vtowfk we have the condition: istep > 3 + dtset%rmm_diis
545 : ! so setting rmm_diis = 1 gives:
546 : ! 4 NSCF iterations for itime == 1
547 : ! 1 NSCF iterations for itime >= 2.
548 3655 : if (scfcv_args%dtset%rmm_diis /= 0 .and. itime == 2) then
549 0 : scfcv_args%dtset%rmm_diis = scfcv_args%dtset%rmm_diis - 3
550 0 : if (scfcv_args%dtset%rmm_diis == 0) scfcv_args%dtset%rmm_diis = 1
551 0 : call wrtout(std_out, sjoin(" itime == 2 with RMM-DIIS --> setting rmm_diis to:", itoa(scfcv_args%dtset%rmm_diis)))
552 : end if
553 :
554 3655 : if (need_verbose .and. .not.need_verbose_cycle) then
555 16524 : write(msg,'(2a,i0,a,i0,2a,80a)') ch10,'--- Iteration: (',itime,'/',ntime,')',ch10,('-',kk=1,80)
556 612 : call wrtout([std_out, ab_out], msg)
557 : end if
558 :
559 : ! ###########################################################
560 : ! ### 09. Loop for icycle (From 1 to ncycle)
561 12503 : do icycle=1,ncycle
562 :
563 11496 : call yaml_iterstart("icycle", icycle, dev_null, scfcv_args%dtset%use_yaml)
564 :
565 11496 : itime_hist = (itime-1)*ncycle + icycle ! Store the time step in the history
566 :
567 : ! ###########################################################
568 : ! ### 10. Output for each icycle (and itime)
569 11496 : if(need_verbose.and.need_verbose_cycle)then
570 : write(msg,fmt)&
571 282204 : ch10,'--- Iteration: (',itime,'/',ntime,') Internal Cycle: (',icycle,'/',ncycle,')',ch10,('-',kk=1,80)
572 10452 : call wrtout([std_out, ab_out], msg)
573 : end if
574 11496 : if (useprtxfase) call prtxfase(ab_mover,hist,itime_hist,std_out,mover_BEFORE)
575 :
576 35517980 : xred_prev(:,:)=xred(:,:)
577 11496 : rprimd_prev(:,:)=rprimd(:,:)
578 :
579 : ! ###########################################################
580 : ! ### 11. Symmetrize atomic coordinates over space group elements
581 :
582 : call symmetrize_xred(ab_mover%natom,&
583 11496 : scfcv_args%dtset%nsym,scfcv_args%dtset%symrel,scfcv_args%dtset%tnons,xred,indsym=scfcv_args%indsym)
584 :
585 35516132 : changed = any(xred /= xred_prev)
586 11496 : if (changed)then
587 1909 : hist%xred(:,:,hist%ihist)=xred(:,:)
588 49 : if(need_verbose) then
589 49 : write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED'
590 49 : write(std_out,*) 'DIFFERENCES:'
591 514 : do kk=1,ab_mover%natom
592 1909 : write(std_out,*) xred(:,kk)-xred_prev(:,kk)
593 : end do
594 : end if
595 1909 : xred_prev(:,:)=xred(:,:)
596 : end if
597 :
598 : ! ###########################################################
599 : ! ### 12. => Call to SCFCV routine and fill history with forces
600 11496 : if (need_verbose.and.need_verbose_cycle) then
601 3484 : if (need_scfcv_cycle) then
602 : write(msg,'(a,3a,33a,44a)')&
603 78432 : ch10,('-',kk=1,3),'SELF-CONSISTENT-FIELD CONVERGENCE',('-',kk=1,44)
604 : else
605 : write(msg,'(a,3a,33a,44a)')&
606 88800 : ch10,('-',kk=1,3),'EFFECTIVE POTENTIAL CALCULATION',('-',kk=1,44)
607 : end if
608 10452 : call wrtout([std_out, ab_out], msg)
609 : end if
610 :
611 11496 : if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then
612 :
613 2 : call abihist_compare_and_copy(hist_prev,hist,ab_mover%natom,similar,tol8,specs%nhist==nhisttot,force_hist_copy)
614 2 : hist_prev%ihist=hist_prev%ihist+1
615 :
616 2 : if (hist_prev%ihist==hist_prev%mxhist) re_init_rho=.TRUE.
617 :
618 : else
619 11494 : scfcv_args%ndtpawuj=0
620 11494 : iapp=itime
621 11494 : if(icycle>1.and.icycle/=ncycle) iapp=-1
622 11494 : if(itime==1 .and. icycle/=ncycle ) iapp=-icycle-1
623 11494 : if (ab_mover%ionmov==14.and.(icycle<ncycle)) iapp=-1
624 :
625 : #if defined HAVE_LOTF
626 : if (ab_mover%ionmov/=23 .or.(lotf_extrapolation(itime).and.(icycle/=1.or.itime==1)))then
627 : #endif
628 : !call scfcv_new2(scfcv_args,electronpositron,rhog,rhor,rprimd,xred,xred_old,conv_retcode)
629 :
630 : !WVL - reformat the wavefunctions in the case of xred != xred_old
631 35517946 : if (scfcv_args%dtset%usewvl == 1 .and. maxval(xred_old - xred) > zero) then
632 : ! Before running scfcv, on non-first geometry step iterations,
633 : ! we need to reformat the wavefunctions, taking into acount the new
634 : ! coordinates. We prepare to change rhog (to be removed) and rhor.
635 0 : ABI_FREE(rhog)
636 0 : ABI_FREE(rhor)
637 : call wvl_wfsinp_reformat(scfcv_args%dtset, scfcv_args%mpi_enreg,&
638 0 : & scfcv_args%psps, rprimd, scfcv_args%wvl, xred, xred_old)
639 0 : scfcv_args%nfftf = scfcv_args%dtset%nfft
640 0 : ABI_MALLOC(rhog,(2, scfcv_args%dtset%nfft))
641 0 : ABI_MALLOC(rhor,(2, scfcv_args%dtset%nfft))
642 : call wvl_mkrho(scfcv_args%dtset, scfcv_args%irrzon, scfcv_args%mpi_enreg,&
643 0 : & scfcv_args%phnons, rhor,scfcv_args%wvl%wfs,scfcv_args%wvl%den)
644 : end if
645 :
646 : !Do another initialization of rho using the last atomic positions
647 : !to avoid potential problems during restart of MD
648 11494 : if (re_init_rho) then
649 : !Recompute some local quantities required by initro
650 1 : call pawfgr_init(pawfgr,scfcv_args%dtset,mgfftf,nfftf,ecut_eff,ecutdg_eff,ngfft,ngfftf)
651 3 : ABI_MALLOC(ph1df,(2,3*(2*mgfftf+1)*scfcv_args%dtset%natom))
652 1 : call getph(scfcv_args%atindx,scfcv_args%dtset%natom,ngfftf(1),ngfftf(2),ngfftf(3),ph1df,xred)
653 1 : call metric(gmet,gprimd,ab_out,rmet,rprimd,ucvol)
654 1 : k0(:)=0.0_dp
655 1 : call getcut(boxcut,ecutdg_eff,gmet,gsqcut_eff,scfcv_args%dtset%iboxcut,ab_out,k0,ngfftf)
656 : !Reinitialise density from current xred
657 : call initro(scfcv_args%atindx,scfcv_args%dtset%densty,gmet,gsqcut_eff,&
658 : & scfcv_args%psps%usepaw,mgfftf,scfcv_args%mpi_enreg,scfcv_args%psps%mqgrid_vl,&
659 : & scfcv_args%dtset%natom,scfcv_args%nattyp,nfftf,ngfftf,scfcv_args%dtset%nspden,&
660 : & scfcv_args%psps%ntypat,scfcv_args%psps,scfcv_args%pawtab,ph1df,&
661 : & scfcv_args%psps%qgrid_vl,rhog,rhor,scfcv_args%dtset%spinat,ucvol,&
662 1 : & scfcv_args%psps%usepaw,scfcv_args%dtset%ziontypat,scfcv_args%dtset%znucl)
663 1 : call pawfgr_destroy(pawfgr)
664 1 : ABI_FREE(ph1df)
665 1 : re_init_rho = .FALSE.
666 : ! Also update scf_history with the correct density for density prediction if required
667 3 : if (scfcv_args%scf_history%history_size>0) scfcv_args%scf_history%atmrho_last(:)=rhor(:,1)
668 : end if
669 :
670 : ! MAIN CALL TO SELF-CONSISTENT FIELD ROUTINE
671 11494 : if (need_scfcv_cycle) then
672 :
673 1632 : call dtfil%init_time(iapp)
674 1632 : itimes(1)=itime; itimes(2)=1
675 1632 : if(present(itimimage_gstate))then
676 1632 : itimes(2)=itimimage_gstate
677 : endif
678 :
679 : !DEBUG
680 : ! write(std_out,'(a,5i4)')' m_mover, before scfcv_run : itimes(1:2)=',itimes(1:2)
681 : !ENDDEBUG
682 1632 : call scfcv_run(scfcv_args, electronpositron, itimes, rhog, rhor, rprimd, xred, xred_old, conv_retcode)
683 1632 : if (conv_retcode == -1) then
684 0 : msg = "Scf cycle returned conv_retcode == -1 (timelimit is approaching), this should not happen inside mover"
685 0 : ABI_WARNING(msg)
686 : end if
687 :
688 : else
689 : ! For monte carlo don't need to recompute energy here (done in pred_montecarlo)
690 9862 : name_file='MD_anharmonic_terms_energy.dat'
691 9862 : if(itime == 1 .and. ab_mover%restartxf==-3)then
692 9 : if(icycle==1)call effective_potential_file_mapHistToRef(effective_potential,hist,comm,scfcv_args%dtset%iatfix,&
693 9 : & need_verbose.and.need_verbose_cycle,sc_size=sc_size)!Map Hist to Ref to order atoms
694 2249 : xred(:,:) = hist%xred(:,:,1) ! Fill xred with new ordering
695 9 : hist%ihist = 1
696 : end if
697 :
698 : #if defined DEV_MS_SCALEUP
699 : !If we a SCALE UP effective electron model give the iteration and set print-options
700 : if(need_elec_eval)then
701 : call global_set_parent_iter(itime)
702 : ! Set all print options to false.
703 : call global_set_print_parameters(geom=.FALSE.,eigvals=.FALSE.,eltic=.FALSE.,&
704 : & orbocc=.FALSE.,bands=.FALSE.)
705 : if(itime == 1 .or. modulo(itime,scup_dtset%scup_printniter) == 0)then
706 : call global_set_print_parameters(scup_dtset%scup_printgeom,scup_dtset%scup_printeigv,scup_dtset%scup_printeltic,&
707 : & scup_dtset%scup_printorbocc,scup_dtset%scup_printbands)
708 : end if
709 : end if
710 : #endif
711 :
712 9862 : effective_potential_verbose = need_verbose .and. (need_verbose_cycle .or. icycle == ncycle)
713 9862 : if(present(multibinit_dtset))then
714 : !TODO: use multibinit_dtset to set the parameters of the effective potential
715 : call effective_potential_evaluate( &
716 : & effective_potential,scfcv_args%results_gs%etotal,scfcv_args%results_gs%fcart,scfcv_args%results_gs%gred,&
717 : & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=effective_potential_verbose,&
718 : & filename=name_file,elec_eval=need_elec_eval,efield_type=multibinit_dtset%efield_type,efield=multibinit_dtset%efield,&
719 : & efield_lambda=multibinit_dtset%efield_lambda,nefield=multibinit_dtset%nefield,&
720 : & efield_period=multibinit_dtset%efield_period,efield_phase=multibinit_dtset%efield_phase,&
721 : & efield_gmean=multibinit_dtset%efield_gmean,efield_gvel=multibinit_dtset%efield_gvel,efield_sigma=multibinit_dtset%efield_sigma,&
722 9862 : & efield_background=multibinit_dtset%efield_background,time=itime*ab_mover%dtion)
723 : else
724 : call effective_potential_evaluate( &
725 : & effective_potential,scfcv_args%results_gs%etotal,scfcv_args%results_gs%fcart,scfcv_args%results_gs%gred,&
726 : & scfcv_args%results_gs%strten,ab_mover%natom,rprimd,xred=xred,verbose=effective_potential_verbose,&
727 0 : & filename=name_file,elec_eval=need_elec_eval,time=itime*ab_mover%dtion)
728 : end if
729 :
730 :
731 :
732 :
733 : ! Check if the simulation did not diverge...
734 9862 : if(itime > 3 .and.ABS(scfcv_args%results_gs%etotal - hist%etot(1)) > 1E5)then
735 : ! We set to false the flag corresponding to the bound
736 0 : effective_potential%anharmonics_terms%bounded = .FALSE.
737 0 : if(need_verbose.and.me==master)then
738 0 : ABI_WARNING("The simulation is diverging, please check your effective potential")
739 : end if
740 : ! Set the flag to finish the simulation
741 0 : iexit=1
742 0 : stat4xml="Failed"
743 : else
744 : ! We set to true the flag corresponding to the bound
745 9862 : effective_potential%anharmonics_terms%bounded = .TRUE.
746 : end if
747 : end if
748 : #if defined HAVE_LOTF
749 : end if
750 : #endif
751 : ! ANOMALOUS SITUATION
752 : ! This is the only case where rprimd could change inside scfcv
753 : ! It generates a weird condition, we start with a certain
754 : ! value for rprimd before scfcv and after we finish with a different value.
755 : ! Notice that normally scfcv should not change rprimd
756 : ! And even worse if optcell==0
757 : ! The solution here is to recompute acell and store these value
758 : ! in the present record even if initially it was not exactly
759 : ! the value entering in scfcv
760 : ! One test case with these condition is bigdft/t10
761 149422 : if (any(rprimd(:,:)/=rprimd_prev(:,:))) then
762 0 : hist%acell(:,hist%ihist)=acell(:)
763 0 : hist%rprimd(:,:,hist%ihist)=rprimd(:,:)
764 : end if
765 :
766 : ! ANOMALOUS SITUATIONS
767 : ! * In ionmov 4 & 5 xred could change inside SCFCV
768 : ! So we need to take the values from the output
769 : !
770 : ! * Inside scfcv_core.F90 there is a call to symmetrize_xred.F90
771 : ! for the first SCF cycle symmetrize_xred could change xred
772 11494 : if (ab_mover%ionmov<10)then
773 17439 : changed = any(xred /= xred_prev)
774 1411 : if (changed)then
775 168 : hist%xred(:,:,hist%ihist)=xred(:,:)
776 16 : if(need_verbose)then
777 16 : write(std_out,*) 'WARNING: ATOMIC COORDINATES WERE SYMMETRIZED AFTER SCFCV'
778 16 : write(std_out,*) 'DIFFERENCES:'
779 54 : do kk=1,ab_mover%natom
780 168 : write(std_out,*) xred(:,kk)-xred_prev(:,kk)
781 : end do
782 : end if
783 : end if
784 : end if
785 :
786 : ! Fill velocities and ionic kinetic energy
787 11494 : call vel2hist(ab_mover%amass,hist,vel,vel_cell)
788 :
789 35517946 : hist%fcart(:,:,hist%ihist)=scfcv_args%results_gs%fcart(:,:)
790 80458 : hist%strten(:,hist%ihist) =scfcv_args%results_gs%strten(:)
791 11494 : hist%etot(hist%ihist) =scfcv_args%results_gs%etotal
792 11494 : hist%entropy(hist%ihist) =scfcv_args%results_gs%energies%entropy
793 11494 : hist%time(hist%ihist) =real(itime,kind=dp)
794 :
795 : ! !######################################################################
796 : end if ! if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1.and.hist_prev%ihist<=hist_prev%mxhist)then
797 :
798 : ! Store trajectory in xfh file
799 11496 : if((ab_xfh%nxfh==0.or.itime/=1)) then
800 34488 : ABI_MALLOC(gred_corrected,(3,scfcv_args%dtset%natom))
801 11496 : call fcart2gred(hist%fcart(:,:,hist%ihist),gred_corrected,rprimd,ab_mover%natom)
802 : ! Get rid of mean force on whole unit cell,
803 : ! but only if no generalized constraints are in effect
804 11496 : if (ab_mover%nconeq==0)then
805 45668 : do ii=1,3
806 45668 : if (ii/=3.or.ab_mover%jellslab==0) then
807 26663404 : gr_avg=sum(gred_corrected(ii,:))/dble(ab_mover%natom)
808 26663405 : gred_corrected(ii,:)=gred_corrected(ii,:)-gr_avg
809 : end if
810 : end do
811 : end if
812 11496 : if (ncycle<10.and.ab_mover%restartxf>=0) then
813 13192 : do ii=1,3
814 42874 : rprim(ii,1:3)=rprimd(ii,1:3)/acell(1:3)
815 : end do
816 :
817 : ! The size of ab_xfh%xfhist is to big for very large supercell.
818 : ! Call it only for specific ionmov
819 14561 : if(any((/2,3,10,11,22/)==ab_mover%ionmov)) then
820 1158 : call xfh_update(ab_xfh,acell,gred_corrected,ab_mover%natom,rprim,hist%strten(:,hist%ihist),xred)
821 : end if
822 : end if
823 11496 : ABI_FREE(gred_corrected)
824 : end if
825 :
826 : ! ###########################################################
827 : ! ### 13. Write the history into the _HIST file
828 : ! ###
829 :
830 11496 : if (need_writeHIST.and.me==master) then
831 10394 : ifirst=merge(0,1,(itime>1.or.icycle>1))
832 : call write_md_hist(hist,filename,ifirst,itime_hist,ab_mover%natom,scfcv_args%dtset%nctime,&
833 10394 : & ab_mover%ntypat,ab_mover%typat,amu_curr,ab_mover%znucl,ab_mover%dtion,scfcv_args%dtset%mdtemp)
834 : end if
835 :
836 : ! ###########################################################
837 : ! ### 14. Output after SCFCV
838 11496 : if(need_verbose.and.need_verbose_cycle.and.need_scfcv_cycle)then
839 122550 : write(msg,'(a,3a,a,72a)')ch10,('-',kk=1,3),'OUTPUT',('-',kk=1,71)
840 4902 : call wrtout([std_out, ab_out], msg)
841 : end if
842 11496 : if (useprtxfase.and..not.ab_mover%use_pimd_routine) then
843 1624 : call prtxfase(ab_mover,hist,itime_hist,ab_out,mover_AFTER)
844 1624 : call prtxfase(ab_mover,hist,itime_hist,std_out,mover_AFTER)
845 : end if
846 :
847 : ! ###########################################################
848 : ! ### 15. => Test Convergence of forces and stresses
849 :
850 11496 : if (itime==ntime.and.icycle==ncycle)then
851 183 : iexit=1
852 183 : stat4xml="Failed"
853 : else
854 11313 : stat4xml="Succeded"
855 : end if
856 :
857 : ! Only if convergence is needed
858 11496 : if(specs%isFconv)then
859 1344 : if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5).or.mod(itime,2)==1)then
860 1330 : if (scfcv_args%dtset%tolmxf/=0)then
861 :
862 : call fconv(hist%fcart(:,:,hist%ihist),&
863 : & scfcv_args%dtset%iatfix, &
864 : & iexit,itime,&
865 : & ab_mover%natom,&
866 : & ntime,&
867 : & ab_mover%optcell,&
868 : & scfcv_args%dtset%strfact,&
869 : & scfcv_args%dtset%strtarget,&
870 : & hist%strten(:,hist%ihist),&
871 : & rprim,&
872 1326 : & scfcv_args%dtset%tolmxf)
873 : else
874 4 : call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde)
875 : end if
876 : end if
877 : end if
878 :
879 11496 : if (itime==ntime.and.icycle==ncycle) iexit=1
880 :
881 : ! ###########################################################
882 : ! ### 16. => Precondition forces, stress and energy
883 : ! ### 17. => Call to each predictor
884 : ! Some MOLDYN algorithms require pimd_param to be initialized
885 11496 : if(ab_mover%use_pimd_routine) then
886 10 : call pimd_init(scfcv_args%dtset,pimd_param,me==master,force_imgmov=9)
887 : end if
888 :
889 : call precpred_1geo(ab_mover,ab_xfh,amu_curr,deloc,&
890 : & scfcv_args%dtset%chkdilatmx,&
891 : & scfcv_args%mpi_enreg%comm_cell,&
892 : & scfcv_args%dtset%dilatmx,dtfil%filnam_ds(4),&
893 : & hist,scfcv_args%dtset%hmctt,&
894 : & icycle,iexit,itime,mttk_vars,&
895 : & scfcv_args%dtset%nctime,ncycle,nerr_dilatmx,scfcv_args%dtset%npsp,ntime,&
896 11496 : & pimd_param,scfcv_args%dtset%rprimd_orig,skipcycle,scfcv_args%dtset%usewvl)
897 :
898 : ! Write MOLDYN netcdf and POSABIN files (done every dtset%nctime time step)
899 : ! This file is not created for multibinit run
900 11496 : if(need_scfcv_cycle .and. (ab_mover%ionmov/=23 .or. icycle==1))then
901 1634 : if (scfcv_args%dtset%nctime>0) then
902 34 : jj=itime; if(hist_prev%mxhist>0.and.ab_mover%restartxf==-1) jj=jj-hist_prev%mxhist
903 34 : if (jj>0) then
904 32 : option=3
905 32 : ihist_prev = abihist_findIndex(hist,-1)
906 : call wrt_moldyn_netcdf(ab_mover%amass,scfcv_args%dtset,jj,option,dtfil%fnameabo_moldyn,&
907 : & scfcv_args%mpi_enreg,scfcv_args%results_gs,&
908 : & hist%rprimd(:,:,ihist_prev),dtfil%unpos,hist%vel(:,:,hist%ihist),&
909 32 : & hist%xred(:,:,ihist_prev))
910 : end if
911 34 : if (iexit==1) hist%ihist=ihist_prev
912 : end if
913 : end if
914 11496 : if(iexit/=0) exit
915 :
916 : ! ###########################################################
917 : ! ### 18. Use the history to extract the new values of acell, rprimd and xred
918 :
919 11103 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
920 : ! /!\ ---- DO NOT CHANGE THESE LINES WITHOUT CORE DEVELOPERS PERMISSION ---- /!\
921 : ! LB-03/2026:
922 : ! A noise can accumulate in acell,rprimd and xred after each iterations,
923 : ! resulting in different results for different MPI processes.
924 : ! This has been observed using threads, but could happen in other contexts.
925 : ! This slowly worsens the ionic dynamics, leading to wrong results after many iterations.
926 : ! So here we compute the mean over all MPI processes to reduce the noise.
927 : ! This error is difficult to test as it is observed in long runs only, so BE VERY CAREFUL.
928 : ! Note : the cost of these MPI communications is negligible.
929 : ! comm = comm_cell
930 11103 : nmpi = xmpi_comm_size(comm)
931 11103 : if (nmpi>1) then
932 4518 : ABI_MALLOC(tmp,(size(xred,1),size(xred,2)))
933 287690 : tmp(:,:) = xred(:,:) / nmpi
934 1506 : call xmpi_sum(tmp,comm,ierr)
935 1506 : if (ierr/=0) then
936 0 : ABI_ERROR("Error in mpi sum (tmp)")
937 : end if
938 287690 : xred(:,:) = tmp(:,:)
939 1506 : ABI_FREE(tmp)
940 1506 : ABI_MALLOC(tmp,(size(rprimd,1),size(rprimd,2)))
941 19578 : tmp(:,:) = rprimd(:,:) / nmpi
942 1506 : call xmpi_sum(tmp,comm,ierr)
943 1506 : if (ierr/=0) then
944 0 : ABI_ERROR("Error in mpi sum (tmp)")
945 : end if
946 19578 : rprimd(:,:) = tmp(:,:)
947 1506 : ABI_FREE(tmp)
948 1506 : ABI_MALLOC(tmp_1d,(size(acell)))
949 6024 : tmp_1d(:) = acell(:) / nmpi
950 1506 : call xmpi_sum(tmp_1d,comm,ierr)
951 1506 : if (ierr/=0) then
952 0 : ABI_ERROR("Error in mpi sum (tmp)")
953 : end if
954 6024 : acell(:) = tmp_1d(:)
955 6024 : ABI_FREE(tmp_1d)
956 : end if
957 : ! /!\--------------------/!\
958 :
959 11103 : if (ab_mover%optcell/=0) then
960 : ! Cell may change
961 :
962 : #ifdef FC_NVHPC
963 : ! Yet another wild NVHPC bug (only on eos_nvhpc_23.9_elpa)
964 : if(wrong) write(100,*) xred
965 : #endif
966 1295 : call matr3inv(rprimd,gprimd)
967 :
968 : ! If metric has changed since the initialization, update the Ylm's
969 1295 : if (scfcv_args%psps%useylm==1)then
970 100 : option=0
971 100 : if (scfcv_args%dtset%iscf>0) option=1
972 : call initylmg(gprimd,&
973 : & scfcv_args%kg,&
974 : & scfcv_args%dtset%kptns,&
975 : & scfcv_args%dtset%mkmem,&
976 : & scfcv_args%mpi_enreg,&
977 : & scfcv_args%psps%mpsang,&
978 : & scfcv_args%dtset%mpw,&
979 : & scfcv_args%dtset%nband,&
980 : & scfcv_args%dtset%nkpt,&
981 : & scfcv_args%npwarr,&
982 : & scfcv_args%dtset%nsppol,&
983 : & option,rprimd,&
984 : & scfcv_args%ylm,&
985 100 : & scfcv_args%ylmgr)
986 : end if
987 :
988 : end if
989 :
990 35494135 : vel(:,:)=hist%vel(:,:,hist%ihist)
991 :
992 : ! vel_cell(3,3)= velocities of cell parameters
993 : ! Not yet used here but compute it for consistency
994 144339 : vel_cell(:,:)=hist%vel_cell(:,:,hist%ihist)
995 11103 : if (ab_mover%ionmov==13 .and. hist%mxhist >= 2) then
996 1026 : if (itime_hist>2) then
997 968 : ihist_prev2 = abihist_findIndex(hist,-2)
998 12584 : vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)- hist%rprimd(:,:,ihist_prev2))/(two*ab_mover%dtion)
999 58 : else if (itime_hist>1) then
1000 29 : ihist_prev = abihist_findIndex(hist,-1)
1001 377 : vel_cell(:,:)=(hist%rprimd(:,:,hist%ihist)-hist%rprimd(:,:,ihist_prev))/(ab_mover%dtion)
1002 : end if
1003 : end if
1004 :
1005 : ! This is needed for some compilers such as
1006 : ! pathscale, g95, xlf that do not exit
1007 : ! from a loop if you change the upper limit
1008 : ! inside
1009 11103 : if (icycle>=ncycle .and. scfcv_args%mpi_enreg%me == 0) then
1010 2255 : if(need_verbose.and.need_verbose_cycle)write(std_out,*) 'EXIT:',icycle,ncycle
1011 : exit
1012 : end if
1013 :
1014 :
1015 8848 : if (need_verbose.and.need_verbose_cycle) then
1016 1040 : write(msg,*) 'ICYCLE',icycle,skipcycle
1017 1040 : call wrtout(std_out,msg)
1018 1040 : write(msg,*) 'NCYCLE',ncycle
1019 1040 : call wrtout(std_out,msg)
1020 : end if
1021 10248 : if (skipcycle) exit
1022 :
1023 : ! ###########################################################
1024 : ! ### 19. End loop icycle
1025 :
1026 : end do ! icycle
1027 :
1028 3655 : if(iexit/=0)exit
1029 :
1030 : ! ###########################################################
1031 : ! ### 20. End loop itime
1032 :
1033 : end do ! itime
1034 :
1035 : ! Call fconv here if we exited due to wall time limit.
1036 0 : if (timelimit_exit==1 .and. specs%isFconv) then
1037 0 : iexit = timelimit_exit
1038 0 : ntime = itime-1
1039 0 : ihist_prev = abihist_findIndex(hist,-1)
1040 0 : if ((ab_mover%ionmov/=4.and.ab_mover%ionmov/=5)) then
1041 0 : if (scfcv_args%dtset%tolmxf/=0)then
1042 : call fconv(hist%fcart(:,:,ihist_prev),&
1043 : & scfcv_args%dtset%iatfix, &
1044 : & iexit, itime,&
1045 : & ab_mover%natom,&
1046 : & ntime,&
1047 : & ab_mover%optcell,&
1048 : & scfcv_args%dtset%strfact,&
1049 : & scfcv_args%dtset%strtarget,&
1050 : & hist%strten(:,ihist_prev),&
1051 : & rprim,&
1052 0 : & scfcv_args%dtset%tolmxf)
1053 : else
1054 0 : call erlxconv(hist,iexit,itime,itime_hist,ntime,scfcv_args%dtset%tolmxde)
1055 : end if
1056 : end if
1057 : end if
1058 :
1059 : ! Avoid pending requests if itime == ntime.
1060 393 : call xmpi_wait(quitsum_request,ierr)
1061 :
1062 : !###########################################################
1063 : !### 21. Set the final values of xred with the last computed values (not the last predicted)
1064 :
1065 393 : call hist2var(acell,hist,ab_mover%natom,rprimd,xred,DEBUG)
1066 23845 : vel(:,:)=hist%vel(:,:,hist%ihist)
1067 :
1068 393 : if (DEBUG .and. ab_mover%ionmov==1)then
1069 0 : write (std_out,*) 'vel'
1070 0 : do kk=1,ab_mover%natom
1071 0 : write (std_out,*) hist%vel(:,kk,hist%ihist)
1072 : end do
1073 : end if
1074 :
1075 : !###########################################################
1076 : !### 22. XML Output at the end
1077 :
1078 : !XML output of the status
1079 393 : if (scfcv_args%mpi_enreg%me == 0 .and. scfcv_args%dtset%prtxml == 1) then
1080 1 : write(ab_xml_out, "(3a)") ' <geometryMinimisation type="',trim(specs%type4xml),'">'
1081 1 : write(ab_xml_out, "(5a)") ' <status cvState="',trim(stat4xml) ,'" stop-criterion="',trim(specs%crit4xml),'" />'
1082 1 : write(ab_xml_out, "(3a)") ' </geometryMinimisation>'
1083 : end if
1084 :
1085 : !###########################################################
1086 : !### 23. Deallocate hist and ab_mover datatypes
1087 :
1088 : !This call is needed to free an internal matrix. However, this is not optimal ...
1089 : !One should instead have a datastructure associated with the preconditioner.
1090 393 : if (ab_mover%goprecon>0) call prec_simple(ab_mover,preconforstr,hist,1,1,1)
1091 :
1092 393 : if (ab_mover%ionmov==13 .or. ab_mover%ionmov==25)then
1093 31 : call mttk_fin(mttk_vars)
1094 : end if
1095 :
1096 393 : if (ab_mover%ionmov==10 .or. ab_mover%ionmov==11) call delocint_fin(deloc)
1097 :
1098 393 : ABI_FREE(xred_prev)
1099 :
1100 393 : call abihist_free(hist)
1101 :
1102 393 : call abihist_free(hist_prev)
1103 393 : call abimover_destroy(ab_mover)
1104 786 : call abiforstr_fin(preconforstr)
1105 :
1106 : contains
1107 : !!***
1108 :
1109 : !!****f* ABINIT/fconv
1110 : !!
1111 : !! NAME
1112 : !! fconv
1113 : !!
1114 : !! FUNCTION
1115 : !! Check maximal absolute value of force (hartree/bohr) against input tolerance; if below tolerance, return iexit=1.
1116 : !! Takes into account the fact that the Broyden (or moldyn) step
1117 : !! might be the last one (last itime), to print eventually modified message.
1118 : !! Stresses are also included in the check, provided that optcell/=0.
1119 : !!
1120 : !! If optcell=1, takes only the trace into account
1121 : !! optcell=2, takes all components into account
1122 : !! optcell=3, takes traceless stress into account
1123 : !! optcell=4, takes sigma(1,1), sigma(2,1) and sigma(3,1) into account
1124 : !! optcell=5, takes sigma(1,2), sigma(2,2) and sigma(3,2) into account
1125 : !! optcell=6, takes sigma(1,3), sigma(2,3) and sigma(3,3) into account
1126 : !! optcell=7, takes sigma(2,2),(2,3) and (3 3) into account
1127 : !! optcell=8, takes sigma(1,1),(1,3) and (3 3) into account
1128 : !! optcell=9, takes sigma(1,1),(1,2) and (2 2) into account
1129 : !! In the case of stresses, target the tensor strtarget, and take into account the factor strfact
1130 : !!
1131 : !! INPUTS
1132 : !! fcart(3,natom)= forces on atoms in hartree/bohr in cartesian coordinates
1133 : !! iatfix(3,natom)=1 for frozen atom, 0 for unfrozen
1134 : !! itime=current number of Broyden/Moldyn iterations
1135 : !! natom=number of atoms in unit cell
1136 : !! ntime=maximum number of Broyden/Moldyn iterations allowed
1137 : !! optcell=option for taking stresses into account (see above)
1138 : !! strfact=factor that multiplies the stresses when they are compared to forces.
1139 : !! strtarget(6)=components of the target stress tensor (hartree/bohr^3)
1140 : !! strten(6)=components of the stress tensor (hartree/bohr^3)
1141 : !! tolmxf=tolerance on maximal absolute value of components of forces
1142 : !!
1143 : !! OUTPUT
1144 : !! writes to unit std_out and to ab_out, and returns
1145 : !!
1146 : !! SIDE EFFECTS
1147 : !! Input/Output
1148 : !! at input : iexit= 0 if not the last itime, 1 if the last itime
1149 : !! at output : iexit= 0 if not below tolerance, 1 if below tolerance
1150 : !!
1151 : !! SOURCE
1152 :
1153 1326 : subroutine fconv(fcart,iatfix,iexit,itime,natom,ntime,optcell,strfact,strtarget,strten,rprim,tolmxf)
1154 :
1155 : !Arguments ------------------------------------
1156 : !scalars
1157 : integer,intent(in) :: itime,natom,ntime,optcell
1158 : integer,intent(inout) :: iexit
1159 : real(dp),intent(in) :: strfact,tolmxf
1160 : !arrays
1161 : integer,intent(in) :: iatfix(3,natom)
1162 : real(dp),intent(in) :: fcart(3,natom),strtarget(6),strten(6)
1163 : real(dp), intent(in) :: rprim(3,3)
1164 :
1165 : !Local variables-------------------------------
1166 : !scalars
1167 : integer :: iatom,idir,istr
1168 : real(dp) :: fmax,strdiag,fcell
1169 : character(len=500) :: msg
1170 : !arrays
1171 : real(dp) :: dstr(6)
1172 :
1173 : ! *************************************************************************
1174 :
1175 : ABI_UNUSED(rprim)
1176 :
1177 : !Compute maximal component of forces, EXCLUDING any fixed components
1178 1326 : fmax=zero
1179 5158 : do iatom=1,natom
1180 16654 : do idir=1,3
1181 15328 : if (iatfix(idir,iatom) /= 1) then
1182 10470 : if( abs(fcart(idir,iatom)) >= fmax ) fmax=abs(fcart(idir,iatom))
1183 : end if
1184 : end do
1185 : end do
1186 :
1187 9282 : dstr(:)=strten(:)-strtarget(:)
1188 :
1189 : !Eventually take into account the stress
1190 1326 : if(optcell==1)then
1191 226 : strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp
1192 226 : if(abs(strdiag)*strfact >= fmax ) fmax=abs(strdiag)*strfact
1193 : else if(optcell==2)then
1194 938 : do istr=1,6
1195 938 : if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact
1196 : end do
1197 : else if(optcell==3)then
1198 : ! Must take away the trace from diagonal elements
1199 14 : strdiag=(dstr(1)+dstr(2)+dstr(3))/3.0_dp
1200 56 : do istr=1,3
1201 56 : if(abs(dstr(istr)-strdiag)*strfact >= fmax ) fmax=abs(dstr(istr)-strdiag)*strfact
1202 : end do
1203 56 : do istr=4,6
1204 56 : if(abs(dstr(istr))*strfact >= fmax ) fmax=abs(dstr(istr))*strfact
1205 : end do
1206 : ! else if(optcell==4 .or. optcell==5 .or. optcell==6)then
1207 : ! if(abs(dstr(optcell-3))*strfact >= fmax ) fmax=abs(dstr(optcell-3))*strfact
1208 : else if(optcell==4) then
1209 : ! Only the first lattice vector is relaxed. Check the three stress components
1210 : ! conjugate to its generalized coordinates, consistently with xfpack_f2vout.
1211 50 : fcell = maxval(abs(dstr([1,5,6]))) * strfact
1212 10 : if (fcell >= fmax) fmax=fcell
1213 : else if(optcell==5) then
1214 : ! Only the second lattice vector is relaxed.
1215 0 : fcell = maxval(abs(dstr([2,4,6]))) * strfact
1216 0 : if (fcell >= fmax) fmax=fcell
1217 : else if(optcell==6) then
1218 : ! Only the third lattice vector is relaxed.
1219 35 : fcell = maxval(abs(dstr([3,4,5]))) * strfact
1220 7 : if (fcell >= fmax) fmax=fcell
1221 : else if(optcell==7)then
1222 14 : if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact
1223 14 : if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact
1224 14 : if(abs(dstr(4))*strfact >= fmax ) fmax=abs(dstr(4))*strfact
1225 : else if(optcell==8)then
1226 0 : if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact
1227 0 : if(abs(dstr(3))*strfact >= fmax ) fmax=abs(dstr(3))*strfact
1228 0 : if(abs(dstr(5))*strfact >= fmax ) fmax=abs(dstr(5))*strfact
1229 : else if(optcell==9)then
1230 5 : if(abs(dstr(1))*strfact >= fmax ) fmax=abs(dstr(1))*strfact
1231 5 : if(abs(dstr(2))*strfact >= fmax ) fmax=abs(dstr(2))*strfact
1232 5 : if(abs(dstr(6))*strfact >= fmax ) fmax=abs(dstr(6))*strfact
1233 : end if
1234 :
1235 1326 : if (fmax<tolmxf) then
1236 220 : write(msg, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1237 220 : ' At Broyd/MD step',itime,', gradients are converged : ',ch10,&
1238 440 : ' max grad (force/stress) =',fmax,' < tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1239 660 : call wrtout([std_out, ab_out], msg)
1240 220 : iexit=1
1241 : else
1242 1106 : if(iexit==1)then
1243 96 : write(msg, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1244 96 : ' fconv : WARNING -',ch10,&
1245 96 : ' ntime=',ntime,' was not enough Broyd/MD steps to converge gradients: ',ch10,&
1246 192 : ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1247 288 : call wrtout([std_out, ab_out], msg)
1248 :
1249 96 : write(std_out,"(8a)")ch10,&
1250 96 : "--- !RelaxConvergenceWarning",ch10,&
1251 96 : "message: | ",ch10,TRIM(indent(msg)),ch10,&
1252 192 : "..."
1253 :
1254 : else
1255 : write(msg, '(a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) &
1256 1010 : ' fconv : at Broyd/MD step',itime,', gradients have not converged yet. ',ch10,&
1257 2020 : ' max grad (force/stress) =',fmax,' > tolmxf=',tolmxf,' ha/bohr (free atoms)',ch10
1258 1010 : call wrtout(std_out,msg,'COLL')
1259 : end if
1260 1106 : iexit=0
1261 : end if
1262 :
1263 1326 : end subroutine fconv
1264 : !!***
1265 :
1266 : !!****f* ABINIT/erlxconv
1267 : !! NAME
1268 : !! erlxconv
1269 : !!
1270 : !! FUNCTION
1271 : !! FIXME: add description.
1272 : !!
1273 : !! INPUTS
1274 : !!
1275 : !! OUTPUT
1276 : !!
1277 : !! SOURCE
1278 :
1279 4 : subroutine erlxconv(hist,iexit,itime,itime_hist,ntime,tolmxde)
1280 :
1281 : !Arguments ------------------------------------
1282 : !scalars
1283 : integer,intent(in) :: itime,itime_hist,ntime
1284 : integer,intent(inout) :: iexit
1285 : real(dp), intent(in) :: tolmxde
1286 : !arrays
1287 : type(abihist),intent(inout) :: hist
1288 :
1289 : !Local variables-------------------------------
1290 : integer :: ihist,ihist_prev,ihist_prev2
1291 : real(dp) :: ediff1,ediff2,maxediff
1292 : character(len=500) :: msg
1293 : ! *************************************************************************
1294 :
1295 4 : if (itime_hist<3) then
1296 2 : write(msg, '(a,a,a)' ) ch10,&
1297 4 : ' erlxconv : minimum 3 Broyd/MD steps to check convergence of energy in relaxations',ch10
1298 2 : call wrtout(std_out,msg,'COLL')
1299 : else
1300 2 : ihist = hist%ihist
1301 2 : ihist_prev = abihist_findIndex(hist,-1)
1302 2 : ihist_prev2 = abihist_findIndex(hist,-2)
1303 2 : ediff1 = hist%etot(ihist) - hist%etot(ihist_prev)
1304 2 : ediff2 = hist%etot(ihist) - hist%etot(ihist_prev2)
1305 2 : if ((abs(ediff1)<tolmxde).and.(abs(ediff2)<tolmxde)) then
1306 1 : write(msg, '(a,a,i4,a,a,a,a,a,es11.4,a,a)' ) ch10,&
1307 1 : ' At Broyd/MD step',itime,', energy is converged : ',ch10,&
1308 1 : ' the difference in energy with respect to the two ',ch10,&
1309 2 : ' previous steps is < tolmxde=',tolmxde,' ha',ch10
1310 3 : call wrtout([std_out, ab_out], msg)
1311 1 : iexit=1
1312 : else
1313 1 : maxediff = max(abs(ediff1),abs(ediff2))
1314 1 : if(iexit==1)then
1315 0 : write(msg, '(a,a,a,a,i5,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1316 0 : ' erlxconv : WARNING -',ch10,&
1317 0 : ' ntime=',ntime,' was not enough Broyd/MD steps to converge energy: ',ch10,&
1318 0 : ' max difference in energy =',maxediff,' > tolmxde=',tolmxde,' ha',ch10
1319 0 : call wrtout([std_out, ab_out], msg)
1320 :
1321 0 : write(std_out,"(8a)")ch10,&
1322 0 : "--- !RelaxConvergenceWarning",ch10,&
1323 0 : "message: | ",ch10,TRIM(indent(msg)),ch10,&
1324 0 : "..."
1325 : else
1326 1 : write(msg, '(a,a,i4,a,a,a,es11.4,a,es11.4,a,a)' ) ch10,&
1327 1 : ' erlxconv : at Broyd/MD step',itime,', energy has not converged yet. ',ch10,&
1328 2 : ' max difference in energy=',maxediff,' > tolmxde=',tolmxde,' ha',ch10
1329 1 : call wrtout(std_out,msg,'COLL')
1330 : end if
1331 : end if
1332 : end if
1333 :
1334 4 : end subroutine erlxconv
1335 : !!***
1336 :
1337 : end subroutine mover
1338 : !!***
1339 :
1340 : !!****f* ABINIT/prtxfase
1341 : !!
1342 : !! NAME
1343 : !! prtxfase
1344 : !!
1345 : !! FUNCTION
1346 : !! Print the values of xcart (X), forces (F) acell (A), Stresses (S), and energy (E)
1347 : !! All values come from the history hist
1348 : !! Also compute and print max and rms forces.
1349 : !! Also compute absolute and relative differences with previous calculation
1350 : !!
1351 : !! INPUTS
1352 : !! ab_mover<type abimover>=Subset of dtset only related with movement of ions and acell, contains:
1353 : !! | dtion: Time step
1354 : !! ! natom: Number of atoms
1355 : !! | vis: viscosity
1356 : !! | iatfix: Index of atoms and directions fixed
1357 : !! | amass: Mass of ions
1358 : !! hist<type abihist>=Historical record of positions, forces, stresses, cell and energies,
1359 : !! itime= time step
1360 : !! iout=unit number for printing
1361 : !!
1362 : !! OUTPUT
1363 : !! (only writing)
1364 : !!
1365 : !! SOURCE
1366 :
1367 4882 : subroutine prtxfase(ab_mover,hist,itime,iout,pos)
1368 :
1369 : !Arguments ------------------------------------
1370 : !scalars
1371 : type(abimover),intent(in) :: ab_mover
1372 : type(abihist),intent(in),target :: hist
1373 : integer,intent(in) :: itime,iout
1374 : integer,intent(in) :: pos
1375 : !arrays
1376 :
1377 : !Local variables-------------------------------
1378 : !scalars
1379 : integer :: jj,kk,unfixd,iprt
1380 : real(dp) :: val_max,val_rms,ucvol ! Values maximal and RMS, Volume of Unitary cell
1381 : real(dp) :: dEabs,dErel ! Diff of energy absolute and relative
1382 : real(dp) :: ekin
1383 : real(dp) :: angle(3),rmet(3,3)
1384 : !character(len=80*(max(ab_mover%natom,3)+1)) :: msg
1385 : !MGNAG: This is not very safe. One should use line-based output istead of appending chars
1386 : ! and then outputting everything! For the time being I use this temporary hack to solve the problem with NAG
1387 4882 : character(len=max(80*(max(ab_mover%natom,3)+1),50000)) :: msg
1388 : character(len=18) :: fmt1
1389 : logical :: prtallatoms
1390 : !arrays
1391 9764 : logical :: atlist(ab_mover%natom)
1392 4882 : real(dp),allocatable :: gred(:,:),xcart(:,:)
1393 4882 : real(dp),pointer :: acell(:),fcart(:,:),rprimd(:,:),strten(:),vel(:,:),xred(:,:)
1394 :
1395 : ! ***********************************************************
1396 :
1397 4882 : fmt1='(a,a,1p,3e22.14)'
1398 :
1399 : !##########################################################
1400 : !### 1. Organize list of atoms to print
1401 :
1402 4882 : prtallatoms=.TRUE.
1403 29185 : do kk=1,ab_mover%natom
1404 29185 : if (ab_mover%prtatlist(kk)/=kk) prtallatoms=.FALSE.
1405 : end do
1406 :
1407 29185 : atlist(:)=.FALSE.
1408 29185 : do iprt=1,ab_mover%natom
1409 29185 : if (ab_mover%prtatlist(iprt)>0.and.ab_mover%prtatlist(iprt)<=ab_mover%natom) atlist(ab_mover%prtatlist(iprt))=.TRUE.
1410 : end do
1411 :
1412 4882 : acell => hist%acell(:,hist%ihist)
1413 4882 : rprimd => hist%rprimd(:,:,hist%ihist)
1414 4882 : xred => hist%xred(:,:,hist%ihist)
1415 4882 : fcart => hist%fcart(:,:,hist%ihist)
1416 4882 : strten => hist%strten(:,hist%ihist)
1417 4882 : vel => hist%vel(:,:,hist%ihist)
1418 :
1419 : !###########################################################
1420 : !### 1. Positions
1421 :
1422 14646 : ABI_MALLOC(xcart,(3,ab_mover%natom))
1423 4882 : call xred2xcart(ab_mover%natom,rprimd,xcart,xred)
1424 :
1425 4882 : write(msg, '(a,a)' )ch10,' Cartesian coordinates (xcart) [bohr]'
1426 4882 : call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xcart)
1427 :
1428 4882 : write(msg, '(a)' )' Reduced coordinates (xred)'
1429 4882 : call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,xred)
1430 :
1431 4882 : ABI_FREE(xcart)
1432 :
1433 : !###########################################################
1434 : !### 2. Forces
1435 :
1436 4882 : if(pos==mover_AFTER)then
1437 :
1438 9744 : ABI_MALLOC(gred,(3,ab_mover%natom))
1439 3248 : call fcart2gred(fcart,gred,rprimd,ab_mover%natom)
1440 :
1441 : ! Compute max |f| and rms f,
1442 : ! EXCLUDING the components determined by iatfix
1443 3248 : val_max=0.0_dp
1444 3248 : val_rms=0.0_dp
1445 3248 : unfixd=0
1446 19410 : do kk=1,ab_mover%natom
1447 67896 : do jj=1,3
1448 64648 : if (ab_mover%iatfix(jj,kk) /= 1) then
1449 45690 : unfixd=unfixd+1
1450 45690 : val_rms=val_rms+fcart(jj,kk)**2
1451 45690 : val_max=max(val_max,abs(fcart(jj,kk)**2))
1452 : end if
1453 : end do
1454 : end do
1455 3248 : if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd))
1456 :
1457 3248 : write(msg, '(a,1p,2e12.5,a)' ) ' Cartesian forces (fcart) [Ha/bohr]; max,rms=',sqrt(val_max),val_rms,' (free atoms)'
1458 3248 : call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,fcart)
1459 :
1460 3248 : write(msg, '(a)' )' Gradient of E wrt nuclear positions in reduced coordinates (gred)'
1461 3248 : call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,gred)
1462 3248 : ABI_FREE(gred)
1463 : end if
1464 :
1465 : !###########################################################
1466 : !### 3. Velocities
1467 :
1468 : !Only if the velocities are being used
1469 4882 : if (hist%isVused)then
1470 : ! Only if velocities are recorded in a history
1471 1219 : if (allocated(hist%vel))then
1472 : ! Compute max |v| and rms v,
1473 : ! EXCLUDING the components determined by iatfix
1474 1219 : val_max=0.0_dp
1475 1219 : val_rms=0.0_dp
1476 1219 : unfixd=0
1477 17038 : do kk=1,ab_mover%natom
1478 64495 : do jj=1,3
1479 63276 : if (ab_mover%iatfix(jj,kk) /= 1) then
1480 45864 : unfixd=unfixd+1
1481 45864 : val_rms=val_rms+vel(jj,kk)**2
1482 45864 : val_max=max(val_max,abs(vel(jj,kk)**2))
1483 : end if
1484 : end do
1485 : end do
1486 1219 : if ( unfixd /= 0 ) val_rms=sqrt(val_rms/dble(unfixd))
1487 :
1488 : write(msg, '(a,1p,2e12.5,a)' ) &
1489 1219 : & ' Cartesian velocities (vel) [bohr*Ha/hbar]; max,rms=',sqrt(val_max),val_rms,' (free atoms)'
1490 1219 : call prtnatom(atlist,iout,msg,ab_mover%natom,prtallatoms,vel)
1491 :
1492 : ! Compute the ionic kinetic energy (no cell shape kinetic energy yet)
1493 1219 : ekin=0.0_dp
1494 17038 : do kk=1,ab_mover%natom
1495 64495 : do jj=1,3
1496 : ! Warning : the fixing of atoms is implemented in reduced
1497 : ! coordinates, so that this expression is wrong
1498 63276 : if (ab_mover%iatfix(jj,kk) == 0) then
1499 45864 : ekin=ekin+0.5_dp*ab_mover%amass(kk)*vel(jj,kk)**2
1500 : end if
1501 : end do
1502 : end do
1503 1219 : write(msg, '(a,1p,e22.14,a)' )' Kinetic energy of ions (ekin) [Ha]=',ekin
1504 1219 : call wrtout(iout,msg,'COLL')
1505 : end if
1506 : end if
1507 :
1508 : !###########################################################
1509 : !### 3. ACELL
1510 :
1511 : !Only if the acell is being used
1512 4882 : if (hist%isARused)then
1513 : ! Only if acell is recorded in a history
1514 1286 : if (allocated(hist%acell))then
1515 1286 : write(msg, '(a)' )' Scale of Primitive Cell (acell) [bohr]'
1516 1286 : write(msg,fmt1)TRIM(msg),ch10,acell(:)
1517 1286 : call wrtout(iout,msg,'COLL')
1518 : end if
1519 : end if
1520 :
1521 : !###########################################################
1522 : !### 4. RPRIMD
1523 :
1524 : !Only if the acell is being used
1525 4882 : if (hist%isARused)then
1526 : ! Only if rprimd is recorded in a history
1527 1286 : if (allocated(hist%rprimd))then
1528 1286 : write(msg, '(a)' )' Real space primitive translations (rprimd) [bohr]'
1529 5144 : do kk=1,3
1530 5144 : write(msg,fmt1)TRIM(msg),ch10,rprimd(:,kk)
1531 : end do
1532 1286 : call wrtout(iout,msg,'COLL')
1533 : end if
1534 : end if
1535 :
1536 : !###########################################################
1537 : !### 5. Unitary cell volume
1538 :
1539 4882 : if (ab_mover%optcell/=0)then
1540 :
1541 : ucvol=&
1542 : & rprimd(1,1)*(rprimd(2,2)*rprimd(3,3)-rprimd(3,2)*rprimd(2,3))+&
1543 : & rprimd(2,1)*(rprimd(3,2)*rprimd(1,3)-rprimd(1,2)*rprimd(3,3))+&
1544 1223 : & rprimd(3,1)*(rprimd(1,2)*rprimd(2,3)-rprimd(2,2)*rprimd(1,3))
1545 :
1546 1223 : write(msg, '(a,1p,e22.14)' )' Unitary Cell Volume (ucvol) [Bohr^3]=',ucvol
1547 1223 : call wrtout(iout,msg,'COLL')
1548 :
1549 : ! ###########################################################
1550 : ! ### 5. Angles and lengths
1551 :
1552 : ! Compute real space metric.
1553 81941 : rmet = MATMUL(TRANSPOSE(rprimd),rprimd)
1554 :
1555 1223 : angle(1)=acos(rmet(2,3)/sqrt(rmet(2,2)*rmet(3,3)))/two_pi*360.0d0
1556 1223 : angle(2)=acos(rmet(1,3)/sqrt(rmet(1,1)*rmet(3,3)))/two_pi*360.0d0
1557 1223 : angle(3)=acos(rmet(1,2)/sqrt(rmet(1,1)*rmet(2,2)))/two_pi*360.0d0
1558 :
1559 1223 : write(msg, '(a)' )' Angles (23,13,12)= [degrees]'
1560 1223 : write(msg,fmt1)TRIM(msg),ch10,angle(:)
1561 1223 : call wrtout(iout,msg,'COLL')
1562 :
1563 1223 : write(msg, '(a)' ) ' Lengths [Bohr]'
1564 1223 : write(msg,fmt1)TRIM(msg),ch10,sqrt(rmet(1,1)),sqrt(rmet(2,2)),sqrt(rmet(3,3))
1565 1223 : call wrtout(iout,msg,'COLL')
1566 :
1567 : ! ###########################################################
1568 : ! ### 5. Stress Tensor
1569 :
1570 1223 : if(pos==mover_AFTER)then
1571 : ! Only if strten is recorded in a history
1572 812 : if (allocated(hist%strten))then
1573 :
1574 812 : write(msg, '(a)' ) ' Stress tensor in cartesian coordinates (strten) [Ha/bohr^3]'
1575 :
1576 812 : write(msg,fmt1)TRIM(msg),ch10,strten(1),strten(6),strten(5)
1577 812 : write(msg,fmt1)TRIM(msg),ch10,strten(6),strten(2),strten(4)
1578 812 : write(msg,fmt1)TRIM(msg),ch10,strten(5),strten(4),strten(3)
1579 812 : call wrtout(iout,msg,'COLL')
1580 : end if
1581 : end if
1582 : end if
1583 :
1584 : !###########################################################
1585 : !### 6. Energy
1586 :
1587 9764 : if(pos==mover_AFTER)then
1588 3248 : write(msg, '(a,1p,e22.14)' )' Total energy (etotal) [Ha]=',hist%etot(hist%ihist)
1589 :
1590 3248 : if (itime>1)then
1591 2554 : jj = abihist_findIndex(hist,-1)
1592 2554 : dEabs=hist%etot(hist%ihist)-hist%etot(jj)
1593 2554 : dErel=2*dEabs/(abs(hist%etot(hist%ihist))+abs(hist%etot(jj)))
1594 2554 : write(msg, '(a,a,a,a)' )TRIM(msg),ch10,ch10,' Difference of energy with previous step (new-old):'
1595 : write(msg, '(a,a,10a,a,1p,e12.5,a,10a,a,1p,e12.5)')&
1596 2554 : TRIM(msg),ch10,&
1597 30648 : (' ',jj=1,10),' Absolute (Ha)=',dEabs,ch10,&
1598 58742 : (' ',jj=1,10),' Relative =',dErel
1599 : end if
1600 3248 : call wrtout(iout,msg,'COLL')
1601 : end if
1602 :
1603 : contains
1604 : !!***
1605 :
1606 : !!****f* ABINIT/gettag
1607 : !!
1608 : !! NAME
1609 : !! gettag
1610 : !!
1611 : !! FUNCTION
1612 : !! Set the tag associated to each atom,
1613 : !!
1614 : !! INPUTS
1615 : !! prtallatoms = Logical for PRTint ALL ATOMS
1616 : !! atlist = ATom LIST
1617 : !! index = index for each atom
1618 : !! natom = Number of ATOMs
1619 : !!
1620 : !! OUTPUT
1621 : !! tag = The string to put for each atom
1622 : !!
1623 : !! SOURCE
1624 :
1625 86999 : subroutine gettag(atlist,index,natom,prtallatoms,tag)
1626 :
1627 : !Arguments ------------------------------------
1628 : !scalars
1629 : logical,intent(in) :: prtallatoms
1630 : integer,intent(in) :: natom
1631 : logical,intent(in) :: atlist(natom)
1632 : integer,intent(in) :: index
1633 : character(len=7),intent(out) :: tag
1634 :
1635 : ! *********************************************************************
1636 : !The numbering will be from (1) to (9999)
1637 :
1638 86999 : if (prtallatoms)then
1639 84269 : tag=''
1640 2730 : elseif (atlist(index)) then
1641 2730 : if (natom<10) then
1642 0 : write(tag, '(a,I1.1,a)') ' (',index,')'
1643 2730 : elseif (natom<100) then
1644 2730 : write(tag, '(a,I2.2,a)') ' (',index,')'
1645 0 : elseif (natom<1000) then
1646 0 : write(tag, '(a,I3.3,a)') ' (',index,')'
1647 0 : elseif (natom<10000) then
1648 0 : write(tag, '(a,I4.4,a)') ' (',index,')'
1649 : end if
1650 : end if
1651 :
1652 86999 : end subroutine gettag
1653 : !!***
1654 :
1655 : !!****f* ABINIT/prtnatom
1656 : !!
1657 : !! NAME
1658 : !! prtnatom
1659 : !!
1660 : !! FUNCTION
1661 : !! Print information for N atoms
1662 : !!
1663 : !! INPUTS
1664 : !! prtallatoms = Logical for PRTint ALL ATOMS
1665 : !! atlist = ATom LIST
1666 : !! index = index for each atom
1667 : !! natom = Number of ATOMs
1668 : !!
1669 : !! OUTPUT
1670 : !! tag = The string to put for aech atom
1671 : !!
1672 : !! SOURCE
1673 :
1674 :
1675 17479 : subroutine prtnatom(atlist,iout,message,natom,prtallatoms,thearray)
1676 :
1677 : !Arguments ------------------------------------
1678 : !scalars
1679 : logical,intent(in) :: prtallatoms
1680 : integer,intent(in) :: natom
1681 : logical,intent(in) :: atlist(natom)
1682 : integer,intent(in) :: iout
1683 : character(len=*),intent(inout) :: message
1684 : !arrays
1685 : real(dp) :: thearray(3,natom)
1686 :
1687 : !Local variables-------------------------------
1688 : !scalars
1689 : integer :: kk
1690 : character(len=7) :: tag ! Maximal ' (9999)'
1691 : character(len=18) :: fmt
1692 :
1693 : ! *********************************************************************
1694 :
1695 17479 : fmt='(a,a,1p,3e22.14,a)'
1696 :
1697 114228 : do kk=1,natom
1698 114228 : if (atlist(kk)) then
1699 86999 : call gettag(atlist,kk,natom,prtallatoms,tag)
1700 86999 : write(message,fmt)TRIM(message),ch10,thearray(:,kk),tag
1701 : end if
1702 : end do
1703 17479 : call wrtout(iout,message,'COLL')
1704 :
1705 17479 : end subroutine prtnatom
1706 : !!***
1707 :
1708 : end subroutine prtxfase
1709 : !!***
1710 :
1711 : !!****f* ABINIT/wrt_moldyn_netcdf
1712 : !! NAME
1713 : !! wrt_moldyn_netcdf
1714 : !!
1715 : !! FUNCTION
1716 : !! Write two files for later molecular dynamics analysis:
1717 : !! - MOLDYN.nc (netcdf format) : evolution of key quantities with time (pressure, energy, ...)
1718 : !! - POSABIN : values of coordinates and velocities for the next time step
1719 : !!
1720 : !! INPUTS
1721 : !! amass(natom)=mass of each atom, in unit of electronic mass (=amu*1822...)
1722 : !! dtset <type(dataset_type)>=all input variables for this dataset
1723 : !! itime=time step index
1724 : !! option=1: write MOLDYN.nc file (netcdf format)
1725 : !! 2: write POSABIN file
1726 : !! 3: write both
1727 : !! moldyn_file=name of the MD netcdf file
1728 : !! mpi_enreg=information about MPI parallelization
1729 : !! results_gs <type(results_gs_type)>=results (energy and its components,
1730 : !! forces and its components, the stress tensor) of a ground-state computation
1731 : !! rprimd(3,3)=real space primitive translations
1732 : !! unpos=unit number for POSABIN file
1733 : !! vel(3,natom)=velocities of atoms
1734 : !! xred(3,natom)=reduced coordinates of atoms
1735 : !!
1736 : !! OUTPUT
1737 : !! -- only printing --
1738 : !!
1739 : !! SIDE EFFECTS
1740 : !!
1741 : !! SOURCE
1742 :
1743 32 : subroutine wrt_moldyn_netcdf(amass,dtset,itime,option,moldyn_file,mpi_enreg,&
1744 32 : & results_gs,rprimd,unpos,vel,xred)
1745 :
1746 : use defs_basis
1747 : use defs_abitypes
1748 : use m_results_gs
1749 : use m_abicore
1750 : use m_errors
1751 : use netcdf
1752 :
1753 : use m_io_tools, only : open_file, get_unit
1754 : use m_geometry, only : xcart2xred, xred2xcart, metric
1755 :
1756 : !Arguments ------------------------------------
1757 : !scalars
1758 : integer,intent(in) :: itime,option,unpos
1759 : character(fnlen),intent(in) :: moldyn_file
1760 : type(dataset_type),intent(in) :: dtset
1761 : type(MPI_type),intent(in) :: mpi_enreg
1762 : type(results_gs_type),intent(in) :: results_gs
1763 : !arrays
1764 : real(dp),intent(in) :: amass(dtset%natom),rprimd(3,3)
1765 : real(dp),intent(in),target :: vel(3,dtset%natom)
1766 : real(dp),intent(in) :: xred(3,dtset%natom)
1767 :
1768 : !Local variables-------------------------------
1769 : !scalars
1770 : integer,save :: ipos=0
1771 : integer :: iatom,ii
1772 : character(len=500) :: msg
1773 : integer :: AtomNumDimid,AtomNumId,CelId,CellVolumeId,DimCoordid,DimScalarid,DimVectorid
1774 : integer :: EkinDimid,EkinId,EpotDimid,EpotId,EntropyDimid,EntropyId,MassDimid,MassId,NbAtomsid
1775 : integer :: ncerr,ncid,PosId,StressDimid,StressId,TensorSymDimid
1776 : integer :: TimeDimid,TimestepDimid,TimestepId
1777 : logical :: atom_fix
1778 : real(dp) :: ekin,ucvol
1779 : character(len=fnlen) :: ficname
1780 : character(len=16) :: chain
1781 : !arrays
1782 : integer :: PrimVectId(3)
1783 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
1784 32 : real(dp),allocatable :: xcart(:,:)
1785 32 : real(dp),pointer :: vcart(:,:),vred(:,:),vtmp(:,:)
1786 : ! *************************************************************************
1787 :
1788 : !Only done by master processor, every nctime step
1789 26 : if (mpi_enreg%me==0.and.dtset%nctime>0) then
1790 :
1791 : ! Netcdf file name
1792 26 : ficname = trim(moldyn_file)//'.nc'
1793 :
1794 : ! Xcart from Xred
1795 78 : ABI_MALLOC(xcart,(3,dtset%natom))
1796 26 : call xred2xcart(dtset%natom,rprimd,xcart,xred)
1797 :
1798 : ! ==========================================================================
1799 : ! First time step: write header of netcdf file
1800 : ! ==========================================================================
1801 26 : if (itime==1.and.(option==1.or.option==3)) then
1802 :
1803 5 : ipos=0
1804 :
1805 : ! Write message
1806 5 : write(msg,'(4a)')ch10,' Open file ',trim(ficname),' to store molecular dynamics information.'
1807 5 : call wrtout(std_out,msg,'COLL')
1808 :
1809 : ! Create netcdf file
1810 5 : ncerr = nf90_create(ficname, NF90_CLOBBER , ncid)
1811 5 : NCF_CHECK_MSG(ncerr,'nf90_create')
1812 :
1813 : ! Dimension time for netcdf (time dim is unlimited)
1814 5 : ncerr = nf90_def_dim(ncid, "time", nf90_unlimited, TimeDimid)
1815 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1816 :
1817 : ! Symetric Tensor Dimension
1818 5 : ncerr = nf90_def_dim(ncid, "DimTensor", size(results_gs%strten), TensorSymDimid)
1819 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1820 :
1821 : ! Coordinates Dimension
1822 5 : ncerr = nf90_def_dim(ncid, "DimCoord", size(xcart,1), DimCoordid)
1823 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1824 :
1825 : ! Atoms Dimensions
1826 5 : ncerr = nf90_def_dim(ncid, "NbAtoms", dtset%natom, NbAtomsid)
1827 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1828 :
1829 : ! Vector Dimension
1830 5 : ncerr = nf90_def_dim(ncid, "DimVector", 3 , DimVectorid)
1831 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1832 :
1833 : ! Scalar Dimension
1834 5 : ncerr = nf90_def_dim(ncid, "DimScalar", 1 , DimScalarid)
1835 5 : NCF_CHECK_MSG(ncerr,'nf90_def_dim')
1836 :
1837 : ! Time step and time unit
1838 5 : ncerr = nf90_def_var(ncid, "Time_step", nf90_double , DimScalarid, TimestepDimid)
1839 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1840 5 : ncerr = nf90_put_att(ncid, TimestepDimid, "units", "atomic time unit")
1841 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1842 :
1843 : ! Ionic masses
1844 5 : ncerr = nf90_def_var(ncid, "Ionic_Mass", nf90_double , NbAtomsid, MassDimid)
1845 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1846 5 : ncerr = nf90_put_att(ncid, MassDimid, "units", "atomic mass unit")
1847 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1848 :
1849 : ! Ionic atomic numbers
1850 5 : ncerr = nf90_def_var(ncid, "Ionic_Atomic_Number", nf90_double , NbAtomsid, AtomNumDimid)
1851 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1852 :
1853 : ! E_pot
1854 5 : ncerr = nf90_def_var(ncid, "E_pot", nf90_double , TimeDimid, EpotDimid)
1855 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1856 5 : ncerr = nf90_put_att(ncid, EpotDimid, "units", "hartree")
1857 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1858 :
1859 : ! E_kin
1860 5 : ncerr = nf90_def_var(ncid, "E_kin", nf90_double , TimeDimid, EkinDimid)
1861 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1862 5 : ncerr = nf90_put_att(ncid, EkinDimid, "units", "hartree")
1863 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1864 :
1865 : ! Entropy
1866 5 : ncerr = nf90_def_var(ncid, "Entropy", nf90_double , TimeDimid, EntropyDimid)
1867 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1868 5 : ncerr = nf90_put_att(ncid, EntropyDimid, "units", "")
1869 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1870 :
1871 : ! Stress tensor
1872 15 : ncerr = nf90_def_var(ncid, "Stress", nf90_double , (/TensorSymDimid,TimeDimid/), StressDimid)
1873 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1874 5 : ncerr = nf90_put_att(ncid, StressDimid, "units", "hartree/bohr^3")
1875 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1876 :
1877 : ! Positions
1878 20 : ncerr = nf90_def_var(ncid, "Position", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), PosId)
1879 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1880 5 : ncerr = nf90_put_att(ncid, PosId, "units", "bohr")
1881 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1882 :
1883 : ! Celerities
1884 20 : ncerr = nf90_def_var(ncid, "Celerity", nf90_double ,(/DimCoordid,NbAtomsid,TimeDimid/), CelId)
1885 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1886 5 : ncerr = nf90_put_att(ncid, CelId, "units", "bohr/(atomic time unit)")
1887 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1888 :
1889 : ! In case of volume cell constant
1890 5 : if (dtset%optcell==0) then
1891 : ! Primitive vectors
1892 20 : do ii = 1,3
1893 15 : write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii
1894 15 : ncerr = nf90_def_var(ncid, trim(chain), nf90_double , DimVectorid, PrimVectId(ii))
1895 20 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1896 : end do
1897 : ! Cell Volume
1898 5 : ncerr = nf90_def_var(ncid, "Cell_Volume", nf90_double , DimScalarid, CellVolumeId)
1899 5 : NCF_CHECK_MSG(ncerr,'nf90_def_var')
1900 5 : ncerr = nf90_put_att(ncid, CellVolumeId, "units", "bohr^3")
1901 5 : NCF_CHECK_MSG(ncerr,'nf90_put_att')
1902 : end if
1903 :
1904 : ! Leave define mode and close file
1905 5 : ncerr = nf90_enddef(ncid)
1906 5 : NCF_CHECK_MSG(ncerr,'nf90_enddef')
1907 5 : ncerr = nf90_close(ncid)
1908 5 : NCF_CHECK_MSG(ncerr,'nf90_close')
1909 : end if
1910 :
1911 : ! ==========================================================================
1912 : ! Write data to netcdf file (every nctime time step)
1913 : ! ==========================================================================
1914 26 : if (mod(itime, dtset%nctime)==0.and.(option==1.or.option==3)) then
1915 :
1916 10 : ipos=ipos+1
1917 :
1918 : ! Write message
1919 10 : write(msg,'(3a)')ch10,' Store molecular dynamics information in file ',trim(ficname)
1920 10 : call wrtout(std_out,msg,'COLL')
1921 :
1922 : ! Open netcdf file
1923 10 : ncerr = nf90_open(ficname, nf90_write, ncid)
1924 10 : NCF_CHECK_MSG(ncerr,'nf90_open')
1925 :
1926 : ! Time step
1927 10 : ncerr = nf90_inq_varid(ncid, "Time_step", TimestepId)
1928 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1929 10 : ncerr = nf90_put_var(ncid, TimestepId, dtset%dtion)
1930 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1931 :
1932 : ! Ionic masses
1933 10 : ncerr = nf90_inq_varid(ncid, "Ionic_Mass", MassId)
1934 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1935 20 : ncerr = nf90_put_var(ncid, MassId, amass, start = (/ 1 /), count=(/dtset%natom/))
1936 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1937 :
1938 : ! Ionic atomic numbers
1939 10 : ncerr = nf90_inq_varid(ncid, "Ionic_Atomic_Number", AtomNumId)
1940 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1941 100 : ncerr = nf90_put_var(ncid, AtomNumId, dtset%znucl(dtset%typat(:)),start=(/1/),count=(/dtset%natom/))
1942 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1943 :
1944 : ! Epot
1945 10 : ncerr = nf90_inq_varid(ncid, "E_pot", EpotId)
1946 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1947 30 : ncerr = nf90_put_var(ncid, EpotId, (/results_gs%etotal/), start=(/ipos/),count=(/1/))
1948 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1949 :
1950 : ! Ekin
1951 170 : ekin=zero;atom_fix=(maxval(dtset%iatfix)>0)
1952 10 : if (dtset%ionmov==1.or.(.not.atom_fix)) then
1953 10 : vcart => vel
1954 : else
1955 0 : ABI_MALLOC(vcart,(3,dtset%natom))
1956 0 : ABI_MALLOC(vred,(3,dtset%natom))
1957 0 : vtmp => vel
1958 0 : call xcart2xred(dtset%natom,rprimd,vtmp,vred)
1959 0 : do iatom=1,dtset%natom
1960 0 : do ii=1,3
1961 0 : if (dtset%iatfix(ii,iatom)==1) vred(ii,iatom)=zero
1962 : end do
1963 : end do
1964 0 : call xred2xcart(dtset%natom,rprimd,vcart,vred)
1965 0 : ABI_FREE(vred)
1966 : end if
1967 50 : do iatom=1,dtset%natom
1968 170 : do ii=1,3
1969 160 : ekin=ekin+half*amass(iatom)*vcart(ii,iatom)**2
1970 : end do
1971 : end do
1972 10 : if (dtset%ionmov/=1.and.atom_fix) then
1973 0 : ABI_FREE(vcart)
1974 : end if
1975 10 : ncerr = nf90_inq_varid(ncid, "E_kin", EkinId)
1976 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1977 30 : ncerr = nf90_put_var(ncid, EkinId, (/ekin/), start = (/ipos/),count=(/1/))
1978 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1979 :
1980 : ! EntropyDimid
1981 10 : ncerr = nf90_inq_varid(ncid, "Entropy", EntropyId)
1982 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1983 30 : ncerr = nf90_put_var(ncid, EntropyId, (/results_gs%energies%entropy/),start = (/ipos/),count=(/1/))
1984 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1985 :
1986 : ! Stress tensor
1987 10 : ncerr = nf90_inq_varid(ncid, "Stress", StressId)
1988 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1989 : ncerr = nf90_put_var(ncid, StressId, results_gs%strten, &
1990 30 : & start=(/1,ipos/),count=(/size(results_gs%strten)/))
1991 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1992 :
1993 : ! Positions
1994 10 : ncerr = nf90_inq_varid(ncid, "Position", PosId)
1995 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
1996 : ncerr = nf90_put_var(ncid, PosId, xcart, start=(/1,1,ipos/), &
1997 70 : & count=(/size(xcart,1),dtset%natom,1/))
1998 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
1999 :
2000 : ! Celerities
2001 10 : ncerr = nf90_inq_varid(ncid, "Celerity", CelId)
2002 10 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
2003 : ncerr = nf90_put_var(ncid, CelId, vel, start=(/1,1,ipos/), &
2004 70 : & count=(/size(vel,1),dtset%natom,1/) )
2005 10 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
2006 :
2007 : ! In case of volume cell constant
2008 10 : if (dtset%optcell==0.and.ipos==1) then
2009 : ! Primitive vectors
2010 20 : do ii = 1,3
2011 15 : write(unit=chain,fmt='(a15,i1)') "PrimitiveVector",ii
2012 15 : ncerr = nf90_inq_varid(ncid, trim(chain), PrimVectId(ii) )
2013 15 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
2014 15 : ncerr = nf90_put_var(ncid, PrimVectId(ii), rprimd(:,ii))
2015 20 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
2016 : end do
2017 : ! Cell Volume
2018 5 : call metric(gmet,gprimd,-1,rmet,rprimd,ucvol)
2019 5 : ncerr = nf90_inq_varid(ncid, "Cell_Volume" , CellVolumeId)
2020 5 : NCF_CHECK_MSG(ncerr,'nf90_inq_varid')
2021 5 : ncerr = nf90_put_var(ncid, CellVolumeId, ucvol)
2022 5 : NCF_CHECK_MSG(ncerr,'nf90_put_var')
2023 : end if
2024 :
2025 : ! Close file
2026 10 : ncerr = nf90_close(ncid)
2027 10 : NCF_CHECK_MSG(ncerr,'nf90_close')
2028 : end if
2029 :
2030 : ! ==========================================================================
2031 : ! Write data to POSABIN file (every nctime time step if option=3)
2032 : ! ==========================================================================
2033 26 : if ((mod(itime, dtset%nctime)==0.and.option==3).or.(option==2)) then
2034 :
2035 : ! Open file for writing
2036 10 : if (open_file('POSABIN',msg,unit=unpos,status='replace',form='formatted') /= 0 ) then
2037 0 : ABI_ERROR(msg)
2038 : end if
2039 :
2040 : ! Write Positions
2041 10 : if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'xred ',(xred(ii,1),ii=1,3)
2042 10 : if (dtset%natom>1) then
2043 40 : do iatom=2,dtset%natom
2044 40 : write(unpos,'(7x,3d18.5)') (xred(ii,iatom),ii=1,3)
2045 : end do
2046 : end if
2047 :
2048 : ! Write Velocities
2049 10 : if (dtset%natom>=1) write(unpos,'(a7,3d18.5)') 'vel ',(vel(ii,1),ii=1,3)
2050 10 : if (dtset%natom>1) then
2051 40 : do iatom=2,dtset%natom
2052 40 : write(unpos,'(7x,3d18.5)') (vel(ii,iatom),ii=1,3)
2053 : end do
2054 : end if
2055 :
2056 : ! Close file
2057 10 : close(unpos)
2058 : end if
2059 :
2060 26 : ABI_FREE(xcart)
2061 :
2062 : ! ==========================================================================
2063 : ! End if master proc
2064 : end if
2065 :
2066 58 : end subroutine wrt_moldyn_netcdf
2067 : !!***
2068 :
2069 : end module m_mover
2070 : !!***
|
