Branch data 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 : : !!***
|
