Line data Source code
1 : !!****m* ABINIT/m_fit_polynomial_coeff
2 : !!
3 : !! NAME
4 : !! m_fit_polynomial_coeff
5 : !!
6 : !! FUNCTION
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2010-2026 ABINIT group (AM)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public Licence, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
14 : !!
15 : !! SOURCE
16 :
17 :
18 : #if defined HAVE_CONFIG_H
19 : #include "config.h"
20 : #endif
21 :
22 : #include "abi_common.h"
23 :
24 : module m_fit_polynomial_coeff
25 :
26 : use defs_basis
27 : use m_errors
28 : use m_abicore
29 : use m_polynomial_term
30 : use m_polynomial_coeff
31 : use m_atomdata
32 : use m_xmpi
33 : use m_supercell
34 : use m_fstrings, only : itoa, ftoa
35 :
36 : use m_hashtable_strval, only: hash_table_t
37 : use m_mergesort, only: mergesort
38 : use m_dynamic_array, only : int_array_type, real_array_type
39 : use m_special_funcs,only : factorial
40 : use m_geometry, only : xred2xcart
41 : use m_crystal,only : symbols_crystal
42 : use m_strain,only : strain_type,strain_get
43 : use m_effective_potential,only : effective_potential_type, effective_potential_evaluate
44 : use m_effective_potential,only : effective_potential_freeCoeffs,effective_potential_setCoeffs
45 : use m_effective_potential,only : effective_potential_getDisp, effective_potential_writeAnhHead
46 : use m_effective_potential,only : effective_potential_copy,effective_potential_free,effective_potential_init
47 : use m_effective_potential_file, only : effective_potential_file_mapHistToRef
48 : use m_io_tools, only : open_file,get_unit
49 : use m_abihist, only : abihist,abihist_free,abihist_init,abihist_copy,write_md_hist,var2hist
50 : use m_random_zbq
51 : use m_fit_data
52 : use m_geometry, only: metric
53 : use m_scup_dataset
54 : #if defined DEV_MS_SCALEUP
55 : use scup_global, only : global_set_parent_iter,global_set_print_parameters
56 : #endif
57 :
58 : implicit none
59 :
60 : public :: fit_polynomial_coeff_computeGF
61 : public :: fit_polynomial_coeff_computeMSD
62 : public :: fit_polynomial_coeff_fit
63 : public :: fit_polynomial_coeff_getFS
64 : public :: fit_polynomial_coeff_getPositive
65 : public :: fit_polynomial_coeff_getCoeffBound
66 : public :: fit_polynomial_coeff_solve
67 : public :: fit_polynomial_coeff_testEffPot
68 : public :: fit_polynomial_printSystemFiles
69 : public :: genereList
70 : !!***
71 :
72 : CONTAINS
73 : !===========================================================================================
74 :
75 :
76 : !!****f* m_fit_polynomial_coeff/fit_polynomial_coeff_fit
77 : !!
78 : !! NAME
79 : !! fit_polynomial_coeff_fit
80 : !!
81 : !! FUNCTION
82 : !! Fit the list of coefficients included in eff_pot,
83 : !! if the coefficients are not set in eff_pot, this routine will genenerate
84 : !! a list of coefficients by taking into acount the symmetries of the system
85 : !! and the cutoff
86 : !!
87 : !! INPUTS
88 : !! eff_pot<type(effective_potential)> = effective potential
89 : !! bancoeff(nbancoeff) = list of bannned coeffcients, these coefficients will NOT be
90 : !! used during the fit process
91 : !! fixcoeff(nfixcoeff) = list of fixed coefficient, these coefficients will be
92 : !! imposed during the fit process
93 : !! hist<type(abihist)> = The history of the MD (or snapshot of DFT)
94 : !! generateterm = term to activate the generation of the term set
95 : !! power_disps(2) = array with the minimal and maximal power_disp to be computed
96 : !! nbancoeff = number of banned coeffcients
97 : !! ncycle_in = number of maximum cycle (maximum coefficient to be fitted)
98 : !! nfixcoeff = Number of coefficients imposed during the fit process
99 : !! option = option of the fit process : 1 - selection of the coefficient one by one
100 : !! 2 - selection of the coefficients with Monte Carlo(testversion)
101 : !! comm = MPI communicator
102 : !! cutoff_in = optional,cut off to apply to the range of interation if
103 : !! the coefficient are genereted in this routine
104 : !! max_power_strain = maximum order of the strain of the strain phonon coupling
105 : !! fit_initializeData = optional, logical !If true, we store all the information for the fit,
106 : !! it will reduce the computation time but increase a lot the memory...
107 : !! fit_tolMSDF = optional, tolerance in eV^2/A^2 on the Forces for the fit process
108 : !! fit_tolMSDS = optional, tolerance in eV^2/A^2 on the Stresses for the fit process
109 : !! fit_tolMSDE = optional, tolerance in meV^2/A^2 on the Energy for the fit process
110 : !! fit_tolMSDFS= optional, tolerance in eV^2/A^2 on the Forces+stresses for the fit process
111 : !! positive = optional, TRUE will return only positive coefficients
112 : !! FALSE, default
113 : !! verbose = optional, flag for the verbose mode
114 : !! anhstr = logical, optional : TRUE, the anharmonic strain are computed
115 : !! FALSE, (default) the anharmonic strain are not computed
116 : !! only_odd_power = logical, optional : if TRUE generate only odd power
117 : !! only_even_power= logical, optional : if TRUE generate only even power
118 : !!
119 : !! OUTPUT
120 : !! eff_pot<type(effective_potential)> = effective potential datatype with new fitted coefficients
121 : !!
122 : !! SOURCE
123 :
124 76 : subroutine fit_polynomial_coeff_fit(eff_pot,bancoeff,fixcoeff,hist,generateterm,power_disps,&
125 38 : & nbancoeff,ncoeff_in,nfixcoeff,nimposecoeff,imposecoeff,&
126 : & option,comm,cutoff_in,max_power_strain,initialize_data,&
127 : & fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS,fit_tolGF,&
128 : & positive,verbose,anharmstr,spcoupling,&
129 : & only_odd_power,only_even_power,prt_anh,&
130 : & fit_iatom,prt_files,fit_on,sel_on,fit_factors,prt_GF_csv,&
131 38 : & dispterms,coeff_file_rw,read_effective_potential, max_nbody, &
132 : & drop_rate, ncoeff_per_cycle, fit_weight_T)
133 :
134 : implicit none
135 :
136 : !Arguments ------------------------------------
137 : !scalars
138 : integer,intent(in) :: ncoeff_in,nfixcoeff,nimposecoeff,comm
139 : integer,intent(in) :: generateterm,nbancoeff,option, ncoeff_per_cycle
140 : !arrays
141 : integer,intent(in) :: fixcoeff(nfixcoeff), bancoeff(nbancoeff),imposecoeff(nimposecoeff)
142 : integer,intent(in) :: power_disps(2)
143 : type(effective_potential_type),target,intent(inout) :: eff_pot
144 : type(effective_potential_type),optional,intent(inout) :: read_effective_potential
145 : type(abihist),intent(inout) :: hist
146 : integer,optional,intent(in) :: max_power_strain,prt_anh,fit_iatom, coeff_file_rw, max_nbody(:)
147 : real(dp),optional,intent(in) :: cutoff_in,fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS
148 : real(dp),optional,intent(in) :: fit_tolGF
149 : logical,optional,intent(in) :: verbose,positive,anharmstr,spcoupling
150 : logical,optional,intent(in) :: only_odd_power,only_even_power
151 : logical,optional,intent(in) :: initialize_data,prt_files,prt_GF_csv
152 : logical,optional,intent(in) :: fit_on(3), sel_on(3),dispterms
153 : real(dp),optional,intent(in) :: fit_factors(3)
154 : real(dp), optional, intent(in) :: drop_rate
155 : real(dp), intent(in) :: fit_weight_T
156 : !Local variables-------------------------------
157 : !scalar
158 : integer :: nfix
159 : integer :: ii,icoeff,my_icoeff,icycle,icycle_tmp,ierr,info,index_min,iproc,isweep,jcoeff,ia,generateterm_in
160 : integer :: master,max_power_strain_in,my_rank,my_ncoeff,ncoeff_model,ncoeff_tot,natom_sc,ncell,ncoeff_to_select,ncoeff_tot_tmp
161 : integer :: ncoeff_preselected, ncoeff_selected, ncoeff_to_fit,nproc,ntime,nsweep,size_mpi,ncoeff_fix,ncoeff_out
162 : integer :: my_ncoeff_start,my_ncoeff_end,my_ncoeff_simple,ncoeff_alone
163 : integer :: rank_to_send,unit_anh,fit_iatom_in,unit_GF_val,nfix_and_impose,nfixcoeff_corr,atom_start,atom_end
164 : integer :: ncopy_terms
165 : real(dp) :: cutoff,factor,time,tolMSDF,tolMSDS,tolMSDE,tolMSDFS,tolGF,check_value
166 : real(dp),parameter :: HaBohr_eVAng = Ha_eV /Bohr_Ang
167 2394 : type(effective_potential_type) :: eff_pot_fixed
168 : logical :: iam_master,need_verbose,need_positive,converge,file_opened
169 : logical :: need_anharmstr,need_spcoupling,ditributed_coefficients,need_prt_anh
170 : logical :: need_only_odd_power,need_only_even_power,need_initialize_data
171 : logical :: need_prt_files,need_prt_GF_csv,need_disp
172 : logical :: fit_iatom_all
173 : !arrays
174 : real(dp) :: mingf(4),int_fit_factors(3)
175 : real(dp) :: remaining_rate
176 : integer :: sc_size(3)
177 38 : logical, allocatable :: isbanned(:)
178 38 : logical, allocatable :: isselected(:)
179 38 : logical,allocatable :: fix_and_impose(:)
180 38 : integer,allocatable :: buffsize(:),buffdisp(:),buffin(:),fixcoeff_corr(:)
181 76 : integer,allocatable :: list_coeffs(:),list_coeffs_tmp(:),list_coeffs_tmp2(:)
182 76 : integer,allocatable :: my_coeffindexes(:),singular_coeffs(:)
183 38 : integer,allocatable :: my_coefflist(:) ,stat_coeff(:),list_coeffs_copy(:)
184 38 : real(dp),allocatable :: gf_values_iter(:,:)
185 76 : real(dp),allocatable :: buffGF(:,:),coeff_values(:),energy_coeffs(:,:)
186 38 : real(dp),allocatable :: energy_coeffs_tmp(:,:)
187 76 : real(dp),allocatable :: fcart_coeffs(:,:,:,:),gf_values(:,:),gf_mpi(:,:)
188 38 : real(dp),allocatable :: fcart_coeffs_tmp(:,:,:,:),strten_coeffs_tmp(:,:,:)
189 38 : real(dp),allocatable :: strten_coeffs(:,:,:)
190 38 : real(dp), allocatable :: weights(:)
191 38 : type(polynomial_coeff_type),allocatable :: my_coeffs(:),coeffs_iatom(:)
192 38 : type(polynomial_coeff_type),allocatable :: coeffs_out(:)
193 38 : type(polynomial_coeff_type),target,allocatable :: coeffs_tmp(:)
194 38 : type(polynomial_coeff_type),pointer :: coeffs_in(:)
195 38 : type(fit_data_type) :: fit_data
196 :
197 : !type(real_array_type) :: gf_array
198 : character(len=20000) :: message,message2
199 : character(len=fnlen) :: filename
200 : character(len=3) :: i_char
201 : character(len=7) :: j_char
202 38 : character(len=5),allocatable :: symbols(:)
203 38 : integer, allocatable :: list_bound(:), list_fix(:)
204 : integer :: n_remaining
205 :
206 : ! *************************************************************************
207 :
208 : !MPI variables
209 38 : master = 0
210 38 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
211 38 : iam_master = (my_rank == master)
212 :
213 : #if defined CC_INTEL_ONEAPI
214 : if (iam_master) then
215 : write(message,'(3a)') ch10,' This part of code does not compile with intel IFX compiler'
216 : ABI_ERROR(message)
217 : end if
218 : #else
219 38 : call initialize_parameters()
220 38 : call copy_eff_pot_to_eff_pot_fixed()
221 38 : ncopy_terms = 0
222 38 : call combine_fixcoeff_and_imposecoeff()
223 38 : call print_start_fitting()
224 38 : call map_hist_to_supercell()
225 38 : call generate_list_of_coefficients_to_fit()
226 38 : call read_or_write_coeffs_to_xml()
227 38 : ncoeff_tot = ncoeff_tot + ncoeff_model
228 38 : if((iam_master .and. ncopy_terms > 0)) my_ncoeff = my_ncoeff + ncopy_terms
229 : call get_num_fixed_coeff()
230 38 : call assign_coeffs_to_cpu()
231 38 : call check_sanity_ncoeff_in()
232 38 : call initialize_tags_for_coeffs()
233 : !call select_bounding_terms()
234 38 : call select_fix_terms()
235 38 : call get_ncoeff_preselected()
236 38 : call get_ncoeff_to_select_and_ncoeff_to_fit()
237 38 : call initialize_fitting_parameters_for_coefficients()
238 38 : call get_weight_from_hist(hist, fit_weight_T, ntime, eff_pot%supercell%natom, weights, comm )
239 38 : call initialize_gf()
240 76 : select case(option)
241 : case(1)
242 38 : call select_one_by_one()
243 : case(2)
244 38 : call select_with_monte_carlo()
245 : end select
246 38 : call fit_all_selected_coefficients()
247 38 : call deallocate_arrays()
248 :
249 : contains
250 38 : subroutine initialize_parameters()
251 38 : need_verbose = .TRUE.
252 38 : if(present(verbose)) need_verbose = verbose
253 38 : need_initialize_data = .TRUE.
254 38 : if(present(initialize_data)) need_initialize_data = initialize_data
255 38 : need_positive = .FALSE.
256 38 : if(present(positive)) need_positive = positive
257 38 : need_anharmstr = .FALSE.
258 38 : if(present(anharmstr)) need_anharmstr = anharmstr
259 38 : need_disp = .TRUE.
260 38 : if(present(dispterms)) need_disp = dispterms
261 38 : need_spcoupling = .TRUE.
262 38 : if(present(spcoupling)) need_spcoupling = spcoupling
263 38 : need_only_odd_power = .FALSE.
264 : !----------------
265 38 : generateterm_in = generateterm
266 38 : if(present(coeff_file_rw))then
267 38 : if (coeff_file_rw==2) generateterm_in = 0
268 : end if
269 : !----------------
270 38 : if(present(only_odd_power)) need_only_odd_power = only_odd_power
271 38 : need_prt_anh = .FALSE.
272 38 : if(present(prt_anh))then
273 38 : if(prt_anh == 1) need_prt_anh = .TRUE.
274 : end if
275 38 : need_prt_GF_csv = .FALSE.
276 38 : if(present(prt_GF_csv)) need_prt_GF_csv = prt_GF_csv
277 38 : need_prt_files = .TRUE.
278 38 : if(present(prt_files))need_prt_files=prt_files
279 38 : need_only_even_power = .FALSE.
280 38 : if(present(only_even_power)) need_only_even_power = only_even_power
281 38 : if(need_only_odd_power.and.need_only_even_power)then
282 : write(message, '(3a)' )&
283 0 : & 'need_only_odd_power and need_only_even_power are both true',ch10,&
284 0 : & 'Action: contact abinit group'
285 0 : ABI_ERROR(message)
286 : end if
287 38 : max_power_strain_in = 1
288 38 : if(present(max_power_strain))then
289 38 : max_power_strain_in = max_power_strain
290 : end if
291 38 : if(max_power_strain_in <= 0)then
292 : write(message, '(3a)' )&
293 0 : & 'max_power_strain can not be inferior or equal to zero',ch10,&
294 0 : & 'Action: contact abinit group'
295 0 : ABI_ERROR(message)
296 : end if
297 : !Check which atom to fit, if not present do all atoms
298 38 : if(present(fit_iatom))then
299 38 : fit_iatom_in = fit_iatom
300 : else
301 0 : fit_iatom_in = -1
302 : endif
303 :
304 38 : fit_iatom_all=.False.
305 38 : if(present(fit_iatom)) then
306 38 : fit_iatom_all=(fit_iatom==-2)
307 : end if
308 :
309 : !Set int fit factors to default value if fit factors not present
310 38 : int_fit_factors = (/1,1,1/)
311 38 : if (present(fit_factors)) int_fit_factors = fit_factors
312 :
313 38 : ditributed_coefficients = .true.
314 38 : if(option==2) ditributed_coefficients = .false.
315 :
316 :
317 : !Set the tolerance for the fit
318 38 : tolMSDF=zero;tolMSDS=zero;tolMSDE=zero;tolMSDFS=zero;tolGF=zero
319 38 : if(present(fit_tolMSDF)) tolMSDF = fit_tolMSDF
320 38 : if(present(fit_tolMSDS)) tolMSDS = fit_tolMSDS
321 38 : if(present(fit_tolMSDE)) tolMSDE = fit_tolMSDE
322 38 : if(present(fit_tolMSDFS))tolMSDFS = fit_tolMSDFS
323 38 : if(present(fit_tolGF)) tolGF = fit_tolGF
324 :
325 38 : remaining_rate=0.5
326 38 : if(present(drop_rate)) remaining_rate=1.0_dp-drop_rate
327 38 : ncoeff_selected=0
328 :
329 153 : if(nfixcoeff==-1 .and. any(max_nbody==-1)) then
330 0 : ABI_ERROR("nfixcoeff=-1 and max_nbody -1 cannot be used together in the current version of Multibinit.")
331 : end if
332 38 : end subroutine initialize_parameters
333 :
334 :
335 38 : subroutine deallocate_arrays()
336 : !Deallocation of arrays
337 38 : call fit_data_free(fit_data)
338 : !Deallocate the temporary coefficient
339 :
340 198 : do ii=1,size(coeffs_tmp)
341 198 : call polynomial_coeff_free(coeffs_tmp(ii))
342 : end do
343 198 : ABI_SFREE(coeffs_tmp)
344 :
345 : !do ii =1, my_ncoeff
346 5043 : do ii=1,size(my_coeffs)
347 5043 : call polynomial_coeff_free(my_coeffs(ii))
348 : end do
349 5043 : ABI_SFREE(my_coeffs)
350 :
351 198 : do ii=1,size(coeffs_out)
352 198 : call polynomial_coeff_free(coeffs_out(ii))
353 : end do
354 198 : ABI_SFREE(coeffs_out)
355 :
356 : !Deallocate fixed eff_pot
357 38 : call effective_potential_free(eff_pot_fixed)
358 : !Other deallocations
359 38 : ABI_SFREE(list_bound)
360 38 : ABI_SFREE(list_fix)
361 38 : ABI_SFREE(gf_values_iter)
362 38 : ABI_SFREE(buffsize)
363 38 : ABI_SFREE(buffdisp)
364 38 : ABI_SFREE(buffGF)
365 38 : ABI_SFREE(coeff_values)
366 38 : ABI_SFREE(energy_coeffs)
367 38 : ABI_SFREE(energy_coeffs_tmp)
368 38 : ABI_SFREE(fcart_coeffs)
369 38 : ABI_SFREE(fcart_coeffs_tmp)
370 38 : ABI_SFREE(gf_mpi)
371 38 : ABI_SFREE(gf_values)
372 38 : ABI_SFREE(list_coeffs)
373 38 : ABI_SFREE(list_coeffs_tmp)
374 38 : ABI_SFREE(list_coeffs_tmp2)
375 38 : ABI_SFREE(my_coeffindexes)
376 38 : ABI_SFREE(my_coefflist)
377 38 : ABI_SFREE(singular_coeffs)
378 38 : ABI_SFREE(strten_coeffs)
379 38 : ABI_SFREE(strten_coeffs_tmp)
380 38 : ABI_SFREE(stat_coeff)
381 :
382 38 : ABI_FREE(isselected)
383 38 : ABI_FREE(isbanned)
384 38 : ABI_SFREE(fixcoeff_corr)
385 38 : ABI_SFREE(fix_and_impose)
386 38 : ABI_SFREE(list_coeffs_copy)
387 38 : ABI_FREE(weights)
388 38 : end subroutine deallocate_arrays
389 :
390 : !> Copy effective potential parameters to a fixed version
391 : !>
392 : !> If there are no imposed coefficients, copies the potential without coefficients.
393 : !> If there are imposed coefficients, copies only those coefficients.
394 : !> If nimposecoeff=-1, imposes all coefficients.
395 : !>
396 : !> Sets ncoeff_model based on number of anharmonic terms.
397 : !> Validates that imposed coefficients are within valid range.
398 38 : subroutine copy_eff_pot_to_eff_pot_fixed()
399 : ! if there is no imposed coefficient, just copy the effpot without any coefficient.
400 : ! if impose, copy all the imposed coefficient.
401 : ! if nimposedcoeff=-1, impose all the coefficients.
402 :
403 : !Copy the input effective potential eff_pot to eff_pot fixed
404 : !If nimposecoeff=0 the fixed potential is the harmonic potential
405 38 : ncoeff_model = eff_pot%anharmonics_terms%ncoeff
406 38 : if (nimposecoeff > ncoeff_model)then
407 : write(message,'(2a)') "fit_nimposecoeff is greater then the number of anharmonic terms&
408 0 : & provided by input effective potential."&
409 0 : & ,"Action -> Change fit_nimposecoeff, and fit_imposecoeff in the input"
410 0 : ABI_ERROR(message)
411 : endif
412 : !Impose some coefficients of the input potential
413 38 : if (nimposecoeff > 0)then
414 0 : if (any(imposecoeff > ncoeff_model))then
415 : write(message,'(2a)') "At least one value in fit_imposeccoeff is greater then the number of anharmonic terms&
416 0 : & provided by input effective potential."&
417 0 : & ,"Action -> Change fit_imposecoeff in the input"
418 0 : ABI_ERROR(message)
419 : endif
420 :
421 0 : ABI_MALLOC(coeffs_tmp,(nimposecoeff))
422 : ! Copy the imposed coefficients to coeffs_tmp, then set them in eff_pot_fixed
423 0 : do ia = 1,nimposecoeff
424 0 : ii = imposecoeff(ia)
425 : ! set imposed coefficient to coeff_tmp
426 : call polynomial_coeff_init(eff_pot%anharmonics_terms%coefficients(ii)%coefficient,&
427 : & eff_pot%anharmonics_terms%coefficients(ii)%nterm,coeffs_tmp(ia),&
428 : & eff_pot%anharmonics_terms%coefficients(ii)%terms,&
429 : & eff_pot%anharmonics_terms%coefficients(ii)%name,&
430 0 : & check = .TRUE.)
431 : enddo
432 : !Copy the input eff pot, free the coeffs and set the ones who shall be imposed to fixed
433 0 : call effective_potential_copy(eff_pot_fixed,eff_pot,comm)
434 0 : call effective_potential_freeCoeffs(eff_pot_fixed)
435 0 : call effective_potential_setCoeffs(coeffs_tmp,eff_pot_fixed,nimposecoeff)
436 : ! set the eff_pot_fixed with coeffs_tmp with imposed coefficients
437 : !Deallocate coeffs tmp
438 0 : do ii = 1,nimposecoeff
439 0 : call polynomial_coeff_free(coeffs_tmp(ii))
440 : enddo
441 0 : ABI_SFREE(coeffs_tmp)
442 38 : elseif (nimposecoeff == -1)then
443 0 : call effective_potential_copy(eff_pot_fixed,eff_pot,comm)
444 : else
445 38 : call effective_potential_copy(eff_pot_fixed,eff_pot,comm)
446 38 : call effective_potential_freeCoeffs(eff_pot_fixed)
447 : endif
448 38 : end subroutine copy_eff_pot_to_eff_pot_fixed
449 :
450 :
451 :
452 : !> Combine fixed and imposed coefficients for consistency
453 : !>
454 : !> This routine:
455 : !> - Checks for overlap between fixed and imposed coefficients
456 : !> - Corrects the fixed coefficient list to avoid duplicates
457 : !> - Handles special cases like nfixcoeff=-1
458 : !> - Creates a mapping between original and corrected coefficient indices
459 : !> - Sets up arrays for tracking fixed/imposed status
460 : !> - Updates nfixcoeff_corr and ncopy_terms accordingly
461 38 : subroutine combine_fixcoeff_and_imposecoeff()
462 : !Set consistency between fixcoeff and imposecoeff.
463 38 : if ( nfixcoeff > 0 .and. nimposecoeff >0)then
464 0 : ABI_MALLOC(fix_and_impose,(nfixcoeff))
465 0 : fix_and_impose = .FALSE.
466 0 : do ii = 1,nfixcoeff
467 0 : if (any(imposecoeff == fixcoeff(ii)))then
468 0 : fix_and_impose(ii) = .TRUE.
469 : endif
470 : enddo
471 0 : nfix_and_impose = count(fix_and_impose)
472 0 : nfixcoeff_corr = nfixcoeff - nfix_and_impose
473 0 : ABI_MALLOC(fixcoeff_corr,(nfixcoeff_corr))
474 0 : ia = 1
475 0 : do ii = 1,nfixcoeff
476 0 : if (.not. fix_and_impose(ii))then
477 0 : fixcoeff_corr(ia) = fixcoeff(ii)
478 0 : ia = ia + 1
479 : endif
480 : enddo
481 0 : ncopy_terms = ncoeff_model - nimposecoeff
482 0 : ABI_MALLOC(list_coeffs_copy,(ncopy_terms))
483 0 : ia = 1
484 0 : do ii = 1,ncoeff_model
485 0 : if( .not. any(imposecoeff == ii))then
486 0 : list_coeffs_copy(ia) = ii
487 0 : ia = ia + 1
488 : endif
489 : enddo
490 0 : do ii = 1,ncopy_terms
491 0 : do ia = 1,nfixcoeff_corr
492 0 : if (list_coeffs_copy(ii) == fixcoeff_corr(ia))then
493 0 : fixcoeff_corr(ia) = ii
494 : endif
495 : enddo
496 : enddo
497 38 : elseif (nfixcoeff == -1 .and. nimposecoeff > 0)then
498 0 : ABI_MALLOC(fix_and_impose,(ncoeff_model))
499 0 : fix_and_impose = .FALSE.
500 0 : do ii = 1,ncoeff_model
501 0 : if (any( imposecoeff == ii))then
502 0 : fix_and_impose(ii) = .TRUE.
503 : endif
504 : enddo
505 0 : nfix_and_impose = nimposecoeff
506 0 : ncopy_terms = ncoeff_model - nimposecoeff
507 0 : ABI_MALLOC(fixcoeff_corr,(ncopy_terms))
508 0 : ABI_MALLOC(list_coeffs_copy,(ncopy_terms))
509 0 : ia = 1
510 0 : do ii = 1,ncoeff_model
511 0 : if (.not. fix_and_impose(ii))then
512 0 : fixcoeff_corr(ia) = ia
513 0 : list_coeffs_copy(ia) = ii
514 0 : ia = ia + 1
515 : endif
516 : enddo
517 0 : nfixcoeff_corr = ncoeff_model- nfix_and_impose
518 38 : elseif (nfixcoeff == -1 .and. nimposecoeff ==-1)then
519 0 : nfixcoeff_corr = 0
520 0 : write(message,'(3a)') "nfixcoeff and nimposecoeff are set to -1.",ch10,&
521 0 : & "This does not make sense. nfixcoeff will be set to 0."
522 0 : if(iam_master) ABI_WARNING(message)
523 0 : ncopy_terms = 0
524 0 : ABI_MALLOC(list_coeffs_copy,(ncopy_terms))
525 38 : elseif (nfixcoeff >0 .and. nimposecoeff ==-1)then
526 0 : nfixcoeff_corr = 0
527 0 : write(message,'(3a)') "nfixcoeff is > 0 and nimposecoeff is set to -1.",ch10,&
528 0 : & "This does not make sense. nfixcoeff will be set to 0."
529 0 : if(iam_master) ABI_WARNING(message)
530 0 : ncopy_terms = 0
531 : else
532 38 : nfixcoeff_corr = nfixcoeff
533 114 : ABI_MALLOC(fixcoeff_corr,(nfixcoeff_corr))
534 : !nimposecoeff or nfixcoeff always 0 here so fix_and_impose is empty
535 38 : ABI_MALLOC(fix_and_impose,(0))
536 90 : fixcoeff_corr = fixcoeff
537 38 : ncopy_terms = ncoeff_model
538 114 : ABI_MALLOC(list_coeffs_copy,(ncopy_terms))
539 126 : do ii = 1,ncopy_terms
540 126 : list_coeffs_copy(ii) = ii
541 : enddo
542 : endif
543 38 : end subroutine combine_fixcoeff_and_imposecoeff
544 :
545 :
546 : !> Print the header message for starting the fitting process
547 38 : subroutine print_start_fitting()
548 38 : if(need_verbose) then
549 3078 : write(message,'(a,(80a))') ch10,('-',ii=1,80)
550 38 : call wrtout(ab_out,message,'COLL')
551 38 : call wrtout(std_out,message,'COLL')
552 38 : write(message,'(2a)') ch10,' Starting Fit process'
553 38 : call wrtout(ab_out,message,'COLL')
554 38 : call wrtout(std_out,message,'COLL')
555 3078 : write(message,'(a,(80a))') ch10,('-',ii=1,80)
556 38 : call wrtout(ab_out,message,'COLL')
557 38 : call wrtout(std_out,message,'COLL')
558 : end if
559 38 : end subroutine print_start_fitting
560 :
561 : !> Map history data to supercell and set up cutoff parameters
562 : !>
563 : !> - Maps history data to supercell if sizes don't match
564 : !> - Sets cutoff distance for interactions
565 : !> - If cutoff is 0, uses default of cell diagonal/3
566 : !> - Calculates supercell size based on cell dimensions
567 38 : subroutine map_hist_to_supercell()
568 : !if the number of atoms in reference supercell into effpot is not correct,
569 : !wrt to the number of atom in the hist, we set map the hist and set the good supercell
570 38 : if (size(hist%xred,2) /= eff_pot%supercell%natom) then
571 0 : call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose)
572 : end if
573 : !Set the cut off
574 38 : cutoff = zero
575 38 : if(present(cutoff_in))then
576 38 : cutoff = cutoff_in
577 : end if
578 : !If the cutoff is set to zero, we define a default value
579 : ! which is the diagonal of the cell/3.
580 38 : if(abs(cutoff)<tol16)then
581 0 : do ii=1,3
582 : cutoff = cutoff + sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
583 : & eff_pot%supercell%rprimd(ii,2)**2+&
584 0 : & eff_pot%supercell%rprimd(ii,3)**2)
585 : end do
586 0 : cutoff = cutoff / 3.0_dp
587 : end if
588 : !we get the size of the supercell in the hist file
589 152 : do ii=1,3
590 : sc_size(ii) = int(anint(sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
591 : & eff_pot%supercell%rprimd(ii,2)**2+&
592 : & eff_pot%supercell%rprimd(ii,3)**2) / &
593 : & sqrt(eff_pot%crystal%rprimd(ii,1)**2+&
594 : & eff_pot%crystal%rprimd(ii,2)**2+&
595 152 : & eff_pot%crystal%rprimd(ii,3)**2)))
596 : end do
597 38 : end subroutine map_hist_to_supercell
598 :
599 : !> Generate or retrieve list of coefficients to be fitted
600 : !>
601 : !> This routine either:
602 : !> - Uses coefficients present in the input effective potential, or
603 : !> - Generates a new list of coefficients based on symmetries and cutoffs
604 : !>
605 : !> For each atom to fit:
606 : !> - Generates coefficients within cutoff radius
607 : !> - Handles special cases like fit_iatom=-2 (fit all atoms)
608 : !> - Validates coefficient generation parameters
609 : !> - Sets up coefficient arrays and counters
610 38 : subroutine generate_list_of_coefficients_to_fit()
611 : integer :: max_nbody_tmp(20)
612 : !Get the list of coefficients to fit:
613 : !get from the eff_pot type (from the input)
614 : !or
615 : !regenerate the list
616 38 : my_ncoeff = 0
617 38 : ncoeff_tot = 0
618 :
619 : !Reset ncoeff_tot
620 38 : if(ncoeff_model > 0)then
621 27 : if(need_verbose)then
622 27 : write(message, '(4a)' )ch10,' The coefficients present in the effective',&
623 54 : & ' potential will be used for the fit'
624 27 : call wrtout(std_out,message,'COLL')
625 27 : call wrtout(ab_out,message,'COLL')
626 : end if
627 : end if
628 :
629 114 : ABI_MALLOC(symbols,(eff_pot%crystal%natom))
630 : call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,&
631 38 : & symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
632 :
633 38 : if (fit_iatom_all .and. generateterm==1) then
634 0 : atom_start = 1
635 0 : atom_end = eff_pot%crystal%nirredat
636 : else
637 38 : atom_start = fit_iatom_in
638 38 : atom_end = fit_iatom_in
639 : end if
640 :
641 38 : if(generateterm == 1)then
642 : ! we need to regenerate them
643 38 : ncoeff_tot = 0
644 76 : do ii = atom_start, atom_end
645 38 : if (fit_iatom_all) then ! fit_iatom=-2
646 0 : fit_iatom_in=eff_pot%crystal%irredatindx(ii)
647 : end if
648 38 : if(need_verbose)then
649 38 : if(fit_iatom_in > 0)then
650 31 : write(message, '(2a,I3,4a)' )ch10,' The coefficients for the fit around atom', fit_iatom_in,': ',&
651 62 : & trim(symbols(fit_iatom_in)),', will be generated',ch10
652 31 : call wrtout(std_out,message,'COLL')
653 31 : call wrtout(ab_out,message,'COLL')
654 : else
655 7 : write(message, '(3a)' )ch10,' The coefficients for the fit will be generated with all cross terms',ch10
656 7 : call wrtout(std_out,message,'COLL')
657 7 : call wrtout(ab_out,message,'COLL')
658 : endif
659 38 : write(message,'(a,F6.3,a)') " Cutoff of ",cutoff," Bohr is imposed"
660 38 : call wrtout(std_out,message,'COLL')
661 : end if !need_verbose
662 :
663 798 : max_nbody_tmp(:) = 20
664 153 : max_nbody_tmp(1: size(max_nbody)) = max_nbody
665 : call polynomial_coeff_getNorder(coeffs_iatom,eff_pot%crystal,cutoff,my_ncoeff,ncoeff_tot_tmp,power_disps,&
666 : & max_power_strain_in,0,sc_size,comm,anharmstr=(ii==1 .and. need_anharmstr),&
667 : & spcoupling=need_spcoupling,distributed=.true.,&
668 : & only_odd_power=need_only_odd_power,&
669 : & only_even_power=need_only_even_power,&
670 : & fit_iatom=fit_iatom_in,dispterms=need_disp, &
671 73 : & max_nbody=max_nbody_tmp)
672 :
673 :
674 76 : if (.not. fit_iatom_all) then
675 38 : call polynomial_coeff_list_free(coeffs_tmp)
676 38 : ncoeff_tot = ncoeff_tot_tmp
677 5071 : ABI_MALLOC(coeffs_tmp,(my_ncoeff))
678 : !if(allocated(coeffs_iatom)) then
679 38 : call coeffs_list_copy(coeffs_tmp,coeffs_iatom)
680 38 : call polynomial_coeff_list_free(coeffs_iatom)
681 : !end if
682 : else
683 0 : ncoeff_tot = ncoeff_tot+ncoeff_tot_tmp
684 0 : if(.not.(allocated(coeffs_tmp))) then
685 0 : ABI_MALLOC(coeffs_tmp,(size(coeffs_iatom)))
686 0 : call coeffs_list_copy(coeffs_tmp,coeffs_iatom)
687 : else
688 : ! FIXME: this does not work on ubu_intel
689 0 : if(allocated(coeffs_iatom)) then
690 0 : call coeffs_list_conc_onsite(coeffs_tmp,coeffs_iatom)
691 : endif
692 : end if ! not allocate coeffs_tmp
693 0 : call polynomial_coeff_list_free(coeffs_iatom)
694 : end if !fit_iatom/=-2
695 : end do ! ii = atom_start, atom_end
696 :
697 : !call coeffs_list_reduce_duplicate(coeffs_tmp, eff_pot%crystal, sc_size, fit_iatom_in, cutoff , power_disps(2))
698 38 : my_ncoeff = size(coeffs_tmp)
699 38 : if (fit_iatom_all .and. iam_master) then
700 : ! FIXME: this does not compile on alpa intel 2025 elpa!
701 0 : message=ch10 // ' fit_iatom = -2 : The total number of coefficients for all atoms are'// itoa(my_ncoeff) // ch10
702 0 : call wrtout(std_out,message,'COLL')
703 0 : call wrtout(ab_out,message,'COLL')
704 : end if
705 : end if ! generateterm == 1
706 38 : ABI_FREE(symbols)
707 38 : end subroutine generate_list_of_coefficients_to_fit
708 :
709 : !> Check if a coefficient is a duplicate of any in the given list
710 : !>
711 : !> This function:
712 : !> - Compares a coefficient against a list of existing coefficients
713 : !> - Checks both direct and inverse term comparisons
714 : !> - Returns true if a match is found, false otherwise
715 : !>
716 : !> @param coeff The coefficient to check
717 : !> @param lcoeffs List of coefficients to compare against
718 : !> @param ncoeff_to_compare Number of coefficients to compare
719 : !> @return found True if duplicate found, false otherwise
720 42 : function is_duplicate_coeff(coeff, lcoeffs, ncoeff_to_compare) result(found)
721 : type(polynomial_coeff_type), intent(in):: coeff, lcoeffs(:)
722 : integer, intent(in) :: ncoeff_to_compare
723 : logical :: found
724 42 : found =.False.
725 42 : do ii=1, ncoeff_to_compare
726 0 : if(terms_compare(coeff%terms(1), lcoeffs(ii)%terms(1))) then
727 : found=.True.
728 : cycle
729 : end if
730 42 : if(terms_compare_inverse(coeff%terms(1), lcoeffs(ii)%terms(1))) then
731 0 : found=.True.
732 : cycle
733 : end if
734 : end do
735 42 : end function is_duplicate_coeff
736 :
737 :
738 : !> Read or write coefficients to XML file based on coeff_file_rw setting
739 : !>
740 : !> This routine:
741 : !> - If coeff_file_rw=1: Writes generated coefficients to XML file
742 : !> - If coeff_file_rw=2: Reads coefficients from XML file
743 : !> - Redistributes coefficients across CPUs for parallel processing
744 : !> - Updates coefficient counts and data structures accordingly
745 : !>
746 : !> Handles MPI distribution of coefficients when reading from file
747 38 : subroutine read_or_write_coeffs_to_xml()
748 :
749 38 : if(present(coeff_file_rw)) then
750 : ! TODO hexu: This functionality needs more test and is now temporarily disabled
751 : ! If test ok, remove setting it to 0.
752 :
753 38 : if (coeff_file_rw==1 .and. generateterm_in == 1 ) then
754 0 : write (filename, "(A9,I2,I2,A4)") "TEST_TERMS", my_rank+1,fit_iatom_in,".xml"
755 0 : call polynomial_coeff_writeXML(coeffs_tmp,my_ncoeff,filename=filename)
756 : end if
757 :
758 38 : if (coeff_file_rw==2) then
759 0 : ncoeff_tot = read_effective_potential%anharmonics_terms%ncoeff
760 0 : if (iam_master) then
761 0 : write(message,'(1a,I7,1a)')' Reading coefficients from file',ncoeff_tot,' Coefficients read'
762 0 : call wrtout(std_out,message,'COLL')
763 0 : write(message,'(1a)')' Redistributing the coefficients over CPUs'
764 0 : call wrtout(std_out,message,'COLL')
765 : end if
766 :
767 0 : ncoeff_alone = mod(ncoeff_tot,nproc)
768 0 : my_ncoeff_simple = int(aint(real(ncoeff_tot,sp)/(nproc)))
769 0 : if(ncoeff_alone == 0 .and. ncoeff_tot >= nproc)then
770 0 : my_ncoeff_start = (my_ncoeff_simple * my_rank) + 1
771 0 : my_ncoeff_end = my_ncoeff_start + my_ncoeff_simple - 1
772 0 : else if(ncoeff_tot < nproc)then
773 0 : if(my_rank + 1 <= ncoeff_tot)then !myrank smaller than ncombi
774 0 : my_ncoeff_start = my_rank + 1
775 0 : my_ncoeff_end = my_ncoeff_start
776 : else
777 0 : my_ncoeff_start = ncoeff_tot + 1 !myrank bigger than ncombi
778 0 : my_ncoeff_end = ncoeff_tot + 1
779 : endif
780 0 : else if(ncoeff_tot > nproc .and. ncoeff_alone /= 0)then
781 0 : if(my_rank >= (nproc-ncoeff_alone)) then
782 0 : my_ncoeff_start = (my_ncoeff_simple * my_rank) + 1 + (my_rank - nproc + ncoeff_alone)
783 0 : my_ncoeff_end = my_ncoeff_start + my_ncoeff_simple
784 : else
785 0 : my_ncoeff_start = (my_ncoeff_simple * my_rank) + 1
786 0 : my_ncoeff_end = my_ncoeff_start + my_ncoeff_simple - 1
787 : endif
788 : end if
789 0 : if(my_ncoeff_end <= ncoeff_tot)then
790 0 : my_ncoeff = my_ncoeff_end+1-my_ncoeff_start
791 : else
792 0 : my_ncoeff = 0
793 : endif
794 0 : call polynomial_coeff_list_free(coeffs_tmp)
795 0 : ABI_MALLOC(coeffs_tmp,(my_ncoeff))
796 0 : if(my_ncoeff /= 0)coeffs_tmp = read_effective_potential%anharmonics_terms%coefficients(my_ncoeff_start:my_ncoeff_end)
797 : end if
798 : end if
799 38 : end subroutine read_or_write_coeffs_to_xml
800 :
801 : !> Distribute coefficients across CPUs for parallel processing
802 : !>
803 : !> This routine:
804 : !> - Determines coefficient distribution across MPI ranks
805 : !> - Creates mapping arrays for coefficient indexing
806 : !> - Copies/initializes coefficients on each CPU
807 : !> - Handles special case for CPU0 which gets input coefficients
808 : !> - Sets up MPI communication patterns
809 : !> - Updates global coefficient arrays and indexing
810 : !> - Frees temporary storage after distribution
811 38 : subroutine assign_coeffs_to_cpu()
812 : !Get the list with the number of coeff on each CPU
813 : !In order to be abble to compute the my_coeffindexes array which is for example:
814 : ! if CPU0 has 200 Coeff and CPU1 has 203 Coeff then
815 : ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403
816 : !Also fill the my_coeffs array with the generated coefficients and/or the coefficient from the input xml
817 114 : ABI_MALLOC(buffin,(nproc))
818 132 : buffin = 0
819 38 : buffin(my_rank+1) = my_ncoeff
820 : call xmpi_sum(buffin,comm,ierr)
821 114 : ABI_MALLOC(my_coeffindexes,(my_ncoeff))
822 76 : ABI_MALLOC(my_coefflist,(my_ncoeff))
823 5119 : ABI_MALLOC(my_coeffs,(my_ncoeff))
824 5043 : do icoeff=1,my_ncoeff
825 5005 : jcoeff = icoeff
826 5005 : my_coefflist(icoeff) = icoeff
827 :
828 5005 : if(my_rank==0) then
829 3850 : my_coeffindexes(icoeff) = icoeff
830 : else
831 3003 : my_coeffindexes(icoeff) = sum(buffin(1:my_rank)) + icoeff
832 : end if
833 :
834 : ! Only copy the input coefficients on the CPU0
835 5005 : if(my_rank==0) then
836 3850 : if(icoeff <= ncopy_terms)then
837 48 : coeffs_in => eff_pot%anharmonics_terms%coefficients
838 48 : jcoeff = list_coeffs_copy(icoeff)
839 : else
840 3802 : coeffs_in => coeffs_tmp
841 3802 : jcoeff = jcoeff - ncopy_terms
842 : end if
843 : else
844 1155 : coeffs_in => coeffs_tmp
845 : end if
846 : ! set the my_coeffs.
847 : call polynomial_coeff_init(zero,coeffs_in(jcoeff)%nterm,&
848 : & my_coeffs(icoeff),coeffs_in(jcoeff)%terms,&
849 : & coeffs_in(jcoeff)%name,&
850 5005 : & check=.true.)
851 5043 : call polynomial_coeff_free(coeffs_in(jcoeff))
852 : end do
853 :
854 : !Deallocation
855 4995 : ABI_SFREE(coeffs_tmp)
856 38 : NULLIFY(coeffs_in)
857 38 : ABI_FREE(buffin)
858 :
859 : !wait everybody
860 38 : call xmpi_barrier(comm)
861 :
862 : !Reset the output (we free the memory)
863 38 : call effective_potential_freeCoeffs(eff_pot)
864 38 : end subroutine assign_coeffs_to_cpu
865 :
866 38 : subroutine get_num_fixed_coeff()
867 : !Get number of fixed coeff
868 38 : ncoeff_fix = 0
869 38 : if(nfixcoeff_corr /=0) then
870 27 : if(nfixcoeff_corr == -1)then
871 20 : ncoeff_fix = ncoeff_model
872 : else
873 7 : ncoeff_fix = nfixcoeff_corr
874 : endif
875 : endif
876 : end subroutine get_num_fixed_coeff
877 :
878 : !> Check if ncoeff_in value is valid and adjust if needed
879 : !>
880 : !> This routine:
881 : !> - Verifies ncoeff_in is not zero or larger than ncoeff_tot
882 : !> - Displays warning if ncoeff_in needs adjustment
883 : !> - Sets ncoeff_in to ncoeff_tot if invalid
884 : !> - Handles special case when nfixcoeff_corr = -1
885 38 : subroutine check_sanity_ncoeff_in()
886 : !Check if ncycle_in is not zero or superior to ncoeff_tot
887 38 : if(need_verbose.and.(ncoeff_in > ncoeff_tot).or.(ncoeff_in<0.and.nfixcoeff_corr /= -1)) then
888 0 : write(message, '(6a,I0,3a)' )ch10,&
889 0 : & ' --- !WARNING',ch10,&
890 0 : & ' The number of coeff requested in the input is not correct.',ch10,&
891 0 : & ' This number will be set to the maximum of coefficients: ',ncoeff_tot,ch10,&
892 0 : & ' ---',ch10
893 0 : call wrtout(std_out,message,"COLL")
894 : end if
895 38 : end subroutine check_sanity_ncoeff_in
896 : !> Initialize coefficient selection status arrays
897 : !>
898 : !> This routine:
899 : !> - Allocates arrays to track coefficient status
900 : !> - Initializes isselected array for tracking selected coefficients
901 : !> - Initializes isbanned array for tracking banned coefficients
902 : !> - Sets initial remaining coefficient count
903 38 : subroutine initialize_tags_for_coeffs()
904 : ! initialize isselected and isbanned
905 114 : ABI_MALLOC(isselected, (ncoeff_tot))
906 76 : ABI_MALLOC(isbanned, (ncoeff_tot))
907 38 : n_remaining=ncoeff_tot
908 8779 : isselected(:)=.False.
909 8779 : isbanned(:)=.False.
910 38 : end subroutine initialize_tags_for_coeffs
911 :
912 :
913 : ! subroutine select_bounding_terms()
914 : ! integer :: ico
915 : ! type(int_array_type) :: ind_bound
916 : ! nbound=0
917 : ! do ico=1, size(my_coeffs)
918 : ! if ( my_coeffs(ico)%isbound==1) then
919 : ! ! If it is already imposed/fixed, it needs not to be added.
920 : ! !if (.not.(any(fixcoeff_corr==my_coeffindexes(ico)))) then
921 : ! nbound =nbound +1
922 : ! call ind_bound%push(my_coeffindexes(ico))
923 : ! isselected(my_coeffindexes(ico))=.True.
924 : ! !end if
925 : !end if
926 : ! end do
927 : ! call xmpi_sum(nbound, comm, ierr)
928 : ! call ind_bound%allgatherv(list_bound, comm, nproc)
929 : ! call xmpi_bcast(list_bound, master, comm, ierr)
930 : ! call ind_bound%finalize()
931 : ! call xmpi_lor(isselected, comm)
932 : ! !if (nbound /=count(isselected)) then
933 : ! ! ABI_ERROR("The number of the bounding term is not consistent with the selected terms.")
934 : ! !end if
935 : ! end subroutine select_bounding_terms
936 :
937 : !> Select and mark fixed coefficient terms
938 : !>
939 : !> This routine:
940 : !> - Handles fixed coefficient selection based on nfixcoeff_corr
941 : !> - For nfixcoeff_corr > 0: Selects specific fixed coefficients
942 : !> - For nfixcoeff_corr = -1: Selects all coefficients up to ncoeff_model
943 : !> - Updates isselected array to mark fixed terms
944 : !> - Validates selection count matches expected number
945 : !> - Uses MPI to synchronize selection across processes
946 38 : subroutine select_fix_terms()
947 : integer :: ico, ifix
948 38 : type(int_array_type) :: ind_fix
949 38 : nfix=0
950 :
951 52 : if(nfixcoeff_corr>0) then
952 163 : do ico=1, size(my_coeffs)
953 156 : ifix=my_coeffindexes(ico)
954 466 : if( any(fixcoeff_corr==ifix) ) then
955 6 : nfix=nfix+1
956 6 : call ind_fix%push(ifix)
957 6 : isselected(ifix)=.True.
958 : end if
959 : end do
960 : call xmpi_sum(nfix, comm, ierr)
961 7 : call ind_fix%allgatherv(list_fix, comm, nproc)
962 : call xmpi_bcast(list_fix, master, comm, ierr)
963 7 : call ind_fix%finalize()
964 7 : call xmpi_lor(isselected, comm)
965 371 : if (nfix/=count(isselected)) then
966 0 : ABI_ERROR("The number of the fixed term is not consistent with the selected terms.")
967 : end if
968 71 : else if(nfixcoeff_corr==-1) then
969 3114 : do ico=1, size(my_coeffs)
970 3094 : ifix=my_coeffindexes(ico)
971 3114 : if( ifix<=ncoeff_model) then
972 36 : nfix=nfix+1
973 36 : call ind_fix%push(ifix)
974 36 : isselected(ifix)=.True.
975 : end if
976 : end do
977 : call xmpi_sum(nfix, comm, ierr)
978 20 : call ind_fix%allgatherv(list_fix, comm, nproc)
979 : call xmpi_bcast(list_fix, master, comm, ierr)
980 20 : call ind_fix%finalize()
981 20 : call xmpi_lor(isselected, comm)
982 5506 : if (nfix/=count(isselected)) then
983 0 : ABI_ERROR("The number of the fixed term is not consistent with the selected terms.")
984 : end if
985 : end if
986 :
987 38 : end subroutine select_fix_terms
988 :
989 :
990 :
991 : !> Calculate number of preselected coefficients and handle banned terms
992 : !>
993 : !> This routine:
994 : !> - Updates ncoeff_model and ncoeff_tot
995 : !> - Initializes ncoeff_preselected counter
996 : !> - Marks banned coefficients in isbanned array
997 : !> - Handles special case when nfixcoeff_corr = -1
998 : !> - Outputs informative messages about coefficient selection
999 : !> - Synchronizes banned coefficient status across MPI processes
1000 38 : subroutine get_ncoeff_preselected()
1001 : !Use fixcoeff
1002 : !ncoeff_preselected store the curent number of coefficient in the model
1003 : !Do not reset this variable...
1004 :
1005 38 : ncoeff_model = eff_pot%anharmonics_terms%ncoeff
1006 38 : ncoeff_tot = ncoeff_tot + ncoeff_model
1007 38 : ncoeff_preselected = 0
1008 :
1009 38 : isbanned(bancoeff(1: nbancoeff))=.True.
1010 38 : call xmpi_lor(isbanned, comm)
1011 :
1012 :
1013 38 : if (nfixcoeff_corr == -1)then
1014 20 : write(message, '(3a)')' nfixcoeff is set to -1, the coefficients present in the model',&
1015 40 : & ' are imposed.',ch10
1016 : else
1017 18 : if (nfixcoeff_corr > 0)then
1018 21 : if(maxval(fixcoeff_corr(:)) > ncoeff_tot) then
1019 0 : write(message, '(4a,I0,6a)' )ch10,&
1020 0 : & ' --- !WARNING',ch10,&
1021 0 : & ' The value ',maxval(fixcoeff_corr(:)),' is not in the list.',ch10,&
1022 0 : & ' Start from scratch...',ch10,&
1023 0 : & ' ---',ch10
1024 : else
1025 7 : write(message, '(2a)')' Some coefficients are imposed from the input.',ch10
1026 : end if
1027 : else
1028 11 : write(message, '(4a)')' There is no coefficient imposed from the input.',ch10,&
1029 22 : & ' Start from scratch',ch10
1030 : end if
1031 : end if
1032 :
1033 38 : ncoeff_preselected = ncoeff_preselected + nfix
1034 :
1035 38 : if(need_verbose) call wrtout(std_out,message,'COLL')
1036 38 : end subroutine get_ncoeff_preselected
1037 :
1038 : !> Calculate number of coefficients to select and fit
1039 : !>
1040 : !> This routine:
1041 : !> - Computes ncoeff_to_select based on ncoeff_in
1042 : !> - Computes ncoeff_to_fit by adding preselected coefficients
1043 : !> - Validates that total coefficients don't exceed maximum allowed
1044 : !> - Adjusts counts if needed and displays warnings
1045 : !> - Handles special case for option 2 (Monte Carlo mode)
1046 : !> - Updates ncoeff_to_select and ncoeff_to_fit accordingly
1047 38 : subroutine get_ncoeff_to_select_and_ncoeff_to_fit()
1048 : !Compute the number of cycle:
1049 38 : ncoeff_to_select = ncoeff_in
1050 38 : if (ncoeff_in + nimposecoeff + nfix > ncoeff_tot .and. generateterm_in==0) then
1051 0 : ncoeff_to_select = ncoeff_tot - nimposecoeff - nfix
1052 0 : write(message, '(4a,I0,2a,I0,2a,I0,3a)' )ch10,&
1053 0 : & ' --- !WARNING',ch10,&
1054 0 : & ' The number of cycle + the number of imposed coefficients: ',ncoeff_to_select,ch10,&
1055 0 : & ' is superior to the maximum number of coefficients in the initial list: ',ncoeff_tot,ch10,&
1056 0 : & ' The number of cycle is set to ',ncoeff_to_select,ch10,&
1057 0 : & ' ---',ch10
1058 0 : if(need_verbose) call wrtout(std_out,message,'COLL')
1059 : end if
1060 : !Compute the maximum number of cycle
1061 38 : ncoeff_to_fit = ncoeff_to_select + ncoeff_preselected
1062 :
1063 :
1064 : !Check if the number of request cycle + the initial number of coeff is superior to
1065 : !the maximum number of coefficient allowed
1066 38 : if(ncoeff_to_fit > ncoeff_tot .and. generateterm_in==0) then
1067 0 : ncoeff_to_select = ncoeff_tot - ncoeff_preselected -nimposecoeff
1068 0 : ncoeff_to_fit = ncoeff_tot-nimposecoeff
1069 :
1070 0 : write(message, '(4a,I0,2a,I0,2a,I0,3a)' )ch10,&
1071 0 : & ' --- !WARNING',ch10,&
1072 0 : & ' The number of cycle + the number of imposed coefficients: ',ncoeff_to_fit,ch10,&
1073 0 : & ' is superior to the maximum number of coefficients in the initial list: ',ncoeff_tot,ch10,&
1074 0 : & ' The number of cycle is set to ',ncoeff_to_select,ch10,&
1075 0 : & ' ---',ch10
1076 0 : if(need_verbose) call wrtout(std_out,message,'COLL')
1077 38 : else if (option==2)then
1078 : ! Always set to the maximum
1079 0 : ncoeff_to_fit = ncoeff_tot
1080 : end if
1081 :
1082 : ! Debug the number of coefficients to fit
1083 :
1084 38 : end subroutine get_ncoeff_to_select_and_ncoeff_to_fit
1085 :
1086 : !> Initialize parameters and arrays needed for coefficient fitting
1087 : !>
1088 : !> This routine:
1089 : !> - Sets up constants for fitting process (ntime, natom_sc, ncell, factor)
1090 : !> - Allocates arrays for energy, forces and stress coefficients
1091 : !> - Initializes list_coeffs array for tracking coefficients
1092 : !> - Handles fixed coefficients from previous selections
1093 : !> - Sets up data structures for decomposing forces and stresses
1094 : !> - Computes constants through fit_data_compute
1095 : !> - Initializes storage for fit results based on memory settings
1096 38 : subroutine initialize_fitting_parameters_for_coefficients()
1097 : !Initialisation of constants
1098 38 : ntime = hist%mxhist
1099 38 : natom_sc = eff_pot%supercell%natom
1100 38 : ncell = eff_pot%supercell%ncells
1101 38 : factor = 1._dp/natom_sc
1102 :
1103 : !Initialisation of arrays:
1104 152 : ABI_MALLOC(energy_coeffs_tmp,(ncoeff_to_fit,ntime))
1105 114 : ABI_MALLOC(list_coeffs,(ncoeff_to_fit))
1106 190 : ABI_MALLOC(fcart_coeffs_tmp,(3,natom_sc,ncoeff_to_fit,ntime))
1107 152 : ABI_MALLOC(strten_coeffs_tmp,(6,ntime,ncoeff_to_fit))
1108 : !ABI_MALLOC(weights, (ntime))
1109 198 : list_coeffs = 0
1110 :
1111 :
1112 :
1113 38 : if(nfix > 0) then
1114 101 : do ii = 1,nfix
1115 101 : list_coeffs(ii) = list_fix(ii)
1116 : end do
1117 : end if
1118 :
1119 : !if ncoeff_preselected > 0 fill list_coeffs with the fixed coefficients
1120 : ! TODO: check if the bounding terms are already in the list of fixed.
1121 : !if(ncoeff_preselected > nbound)then
1122 : ! do ii = nbound+1,ncoeff_preselected
1123 : ! if(nfixcoeff_corr == -1)then
1124 : ! if(ii <= ncoeff_model)then
1125 : ! list_coeffs(ii) = ii
1126 : ! isselected(ii) = .True.
1127 : ! end if
1128 : ! else
1129 : ! list_coeffs(ii) = fixcoeff_corr(ii-nbound)
1130 : ! isselected(fixcoeff_corr(ii-nbound)) = .True.
1131 : ! end if
1132 : ! end do
1133 : !end if
1134 : !call xmpi_lor(isselected, comm)
1135 :
1136 :
1137 : !Get the decomposition for each coefficients of the forces and stresses for
1138 : !each atoms and each step equations 11 & 12 of PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
1139 38 : if(need_verbose)then
1140 38 : write(message, '(a)' ) ' Initialisation of the fit process...'
1141 38 : call wrtout(std_out,message,'COLL')
1142 : end if
1143 : !Before the fit, compute constants with fit_data_compute.
1144 : !Conpute the strain of each configuration.
1145 : !Compute the displacmeent of each configuration.
1146 : !Compute the variation of the displacement due to strain of each configuration.
1147 : !Compute fixed forces and stresse and get the standard deviation.
1148 : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
1149 38 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose)
1150 :
1151 : !Get the decomposition for each coefficients of the forces,stresses and energy for
1152 : !each atoms and each step (see equations 11 & 12 of
1153 : !PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]]+ allocation
1154 : !If the user does not turn off this initialization, we store all the information for the fit,
1155 : !it will reduce the computation time but increase a lot the memory...
1156 38 : if(need_initialize_data)then
1157 0 : ABI_MALLOC(energy_coeffs,(my_ncoeff,ntime))
1158 0 : ABI_MALLOC(fcart_coeffs,(3,natom_sc,my_ncoeff,ntime))
1159 0 : ABI_MALLOC(strten_coeffs,(6,ntime,my_ncoeff))
1160 : call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
1161 : & fit_data%training_set%displacement,&
1162 : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
1163 : & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
1164 0 : & strten_coeffs,fit_data%training_set%ucvol,my_coefflist,my_ncoeff)
1165 : else
1166 : ! Allocate just 1 dimension ! Save MEMORY !
1167 114 : ABI_MALLOC(energy_coeffs,(1,ntime))
1168 152 : ABI_MALLOC(fcart_coeffs,(3,natom_sc,1,ntime))
1169 114 : ABI_MALLOC(strten_coeffs,(6,ntime,1))
1170 : end if
1171 38 : end subroutine initialize_fitting_parameters_for_coefficients
1172 :
1173 : !> Initialize goal function parameters and arrays
1174 : !>
1175 : !> This routine:
1176 : !> - Allocates arrays for tracking coefficients and goal function values
1177 : !> - Sets up MPI communication buffers and patterns
1178 : !> - Handles previously imposed coefficients from input
1179 : !> - Computes initial goal function values for fixed terms
1180 : !> - Initializes arrays for tracking coefficient values and statistics
1181 : !> - Sets up tracking of goal function values during iterations
1182 : !> - Prints initial goal function values if in verbose mode
1183 : !> - Synchronizes data across MPI processes
1184 38 : subroutine initialize_gf()
1185 : integer :: ipre, duplicated, ipre_real
1186 38 : real(dp) , allocatable :: weights1(:)
1187 : !Allocation of arrays
1188 274 : ABI_MALLOC(coeffs_tmp,(ncoeff_to_fit))
1189 114 : ABI_MALLOC(singular_coeffs,(max(1,my_ncoeff)))
1190 114 : ABI_MALLOC(coeff_values,(ncoeff_to_fit))
1191 114 : ABI_MALLOC(gf_values,(4,max(1,my_ncoeff)))
1192 114 : ABI_MALLOC(list_coeffs_tmp,(ncoeff_to_fit))
1193 76 : ABI_MALLOC(list_coeffs_tmp2,(ncoeff_to_fit))
1194 114 : ABI_MALLOC(stat_coeff,(ncoeff_tot))
1195 198 : coeff_values = zero
1196 5043 : singular_coeffs = 0
1197 8779 : stat_coeff = 0
1198 : !Set mpi buffer
1199 : !Set the bufsize for mpi allgather
1200 114 : ABI_MALLOC(buffsize,(nproc))
1201 76 : ABI_MALLOC(buffdisp,(nproc))
1202 38 : ABI_MALLOC(buffGF,(5,1))
1203 114 : ABI_MALLOC(gf_mpi,(5,nproc))
1204 132 : buffsize(:) = 0
1205 38 : buffdisp(1) = 0
1206 :
1207 132 : do ii= 1,nproc
1208 132 : buffsize(ii) = 5
1209 : end do
1210 94 : do ii = 2,nproc
1211 94 : buffdisp(ii) = buffdisp(ii-1) + buffsize(ii-1)
1212 : end do
1213 38 : size_mpi = 5*nproc
1214 : !If some coeff are imposed by the input, we need to fill the arrays
1215 : !with this coeffs and broadcast to the others CPUs :
1216 38 : if(ncoeff_preselected>=1)then
1217 : ipre_real=0
1218 101 : do ipre=1, ncoeff_preselected
1219 74 : list_coeffs_tmp(ipre) = ipre
1220 74 : rank_to_send = 0
1221 74 : duplicated= 0
1222 12226 : do icoeff=1,my_ncoeff
1223 12226 : if((my_coeffindexes(icoeff)==list_coeffs(ipre)))then
1224 42 : if (is_duplicate_coeff(my_coeffs(icoeff), coeffs_tmp, ncoeff_selected)) duplicated=1
1225 : endif
1226 : end do
1227 : call xmpi_sum(duplicated, comm, ierr)
1228 175 : if(duplicated==0) then
1229 : ipre_real=ipre_real+1
1230 3079 : do icoeff=1,my_ncoeff
1231 3079 : if((my_coeffindexes(icoeff)==list_coeffs(ipre)))then
1232 42 : if(need_initialize_data)then
1233 0 : my_icoeff = icoeff
1234 : else
1235 42 : my_icoeff = 1
1236 : ! Need to initialized the data for the fit for this coefficient
1237 : call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
1238 : & fit_data%training_set%displacement,&
1239 : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
1240 : & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
1241 : & strten_coeffs,fit_data%training_set%ucvol,&
1242 42 : & my_coefflist(icoeff),1)
1243 : end if
1244 504 : energy_coeffs_tmp(ipre,:) = energy_coeffs(my_icoeff,:)
1245 74424 : fcart_coeffs_tmp(:,:,ipre,:) = fcart_coeffs(:,:,my_icoeff,:)
1246 3276 : strten_coeffs_tmp(:,:,ipre) = strten_coeffs(:,:,my_icoeff)
1247 :
1248 42 : rank_to_send = my_rank
1249 42 : call polynomial_coeff_free(coeffs_tmp(ipre))
1250 : call polynomial_coeff_init(coeff_values(ipre),my_coeffs(icoeff)%nterm,&
1251 : & coeffs_tmp(ipre),my_coeffs(icoeff)%terms,&
1252 : & my_coeffs(icoeff)%name,&
1253 42 : & check=.false.)
1254 42 : exit
1255 : end if
1256 : end do
1257 : ! Boadcast the coefficient
1258 : ! Need to send the rank with the chosen coefficient
1259 : call xmpi_sum(rank_to_send, comm, ierr)
1260 74 : call xmpi_bcast(energy_coeffs_tmp(ipre,:), rank_to_send, comm, ierr)
1261 74 : call xmpi_bcast(fcart_coeffs_tmp(:,:,ipre,:) , rank_to_send, comm, ierr)
1262 74 : call xmpi_bcast(strten_coeffs_tmp(:,:,ipre), rank_to_send, comm, ierr)
1263 74 : call polynomial_coeff_broadcast(coeffs_tmp(ipre), rank_to_send, comm)
1264 : end if
1265 : end do
1266 : end if
1267 : !Waiting for all
1268 38 : if(nproc > 1) then
1269 24 : if(need_verbose)then
1270 24 : write(message, '(a)') ' Initialisation done... waiting for all the CPU'
1271 24 : call wrtout(std_out,message,'COLL')
1272 : end if
1273 24 : call xmpi_barrier(comm)
1274 : end if
1275 :
1276 : !Compute GF, coeff_values,strten_coeffs and fcart_coeffs are set to zero
1277 : !it means that only the harmonic part wiil be computed
1278 38 : if(fit_weight_T>0.0_dp) then
1279 15 : coeff_values = zero
1280 : call fit_polynomial_coeff_computeGF(coeff_values,energy_coeffs,fit_data%energy_diff,fcart_coeffs,&
1281 : & fit_data%fcart_diff,gf_values(:,1),int((/1/)),natom_sc,&
1282 : & 0,my_ncoeff,ntime,strten_coeffs,fit_data%strten_diff,&
1283 3 : & fit_data%training_set%sqomega, weights=weights)
1284 :
1285 : !Print the standard deviation before the fit
1286 : ! FIXME: on builder ALPS_INTEL_2025_elpa, compile fails.
1287 : ! This is probably due to the fact that the format is too long.
1288 : ! We should split the write in two writes.
1289 3 : write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' ) ch10,&
1290 3 : & ' Weighted Goal function values at the begining of the fit process (eV^2/A^2):',ch10, &
1291 3 : & ' Energy : ',&
1292 3 : & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,&
1293 3 : & ' Forces+Stresses : ',&
1294 3 : & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,&
1295 3 : & ' Forces : ',&
1296 3 : & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,&
1297 3 : & ' Stresses : ',&
1298 6 : & gf_values(3,1)*(HaBohr_eVAng)**2,ch10
1299 3 : if(need_verbose)then
1300 3 : call wrtout(ab_out,message,'COLL')
1301 3 : call wrtout(std_out,message,'COLL')
1302 : end if
1303 : end if
1304 :
1305 114 : ABI_MALLOC(weights1,(ntime))
1306 595 : weights1=one
1307 198 : coeff_values = zero
1308 : call fit_polynomial_coeff_computeGF(coeff_values,energy_coeffs,fit_data%energy_diff,fcart_coeffs,&
1309 : & fit_data%fcart_diff,gf_values(:,1),int((/1/)),natom_sc,&
1310 : & 0,my_ncoeff,ntime,strten_coeffs,fit_data%strten_diff,&
1311 38 : & fit_data%training_set%sqomega, weights=weights1)
1312 :
1313 : !Print the standard deviation before the fit
1314 : ! FIXME: on builder ALPS_INTEL_2025_elpa, compile fails.
1315 : ! This is probably due to the fact that the format is too long.
1316 : ! We should split the write in two writes.
1317 38 : write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' ) ch10,&
1318 38 : & ' Goal function values at the begining of the fit process (eV^2/A^2):',ch10,&
1319 38 : & ' Energy : ',&
1320 38 : & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,&
1321 38 : & ' Forces+Stresses : ',&
1322 38 : & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,&
1323 38 : & ' Forces : ',&
1324 38 : & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,&
1325 38 : & ' Stresses : ',&
1326 76 : & gf_values(3,1)*(HaBohr_eVAng)**2,ch10
1327 : !message = ' Goal function values at the begining of the fit process (eV^2/A^2):'//ch10
1328 : !message = trim(message) // ' Energy : ' // trim(ftoa(gf_values(4,1)*(HaBohr_eVAng)**2))//ch10
1329 : !message = trim(message) // ' Forces+Stresses : ' // trim(ftoa(gf_values(1,1)*(HaBohr_eVAng)**2))//ch10
1330 : !message = trim(message) // ' Forces : ' // trim(ftoa(gf_values(2,1)*(HaBohr_eVAng)**2))//ch10
1331 : !message = trim(message) // ' Stresses : ' // trim(ftoa(gf_values(3,1)*(HaBohr_eVAng)**2))//ch10
1332 38 : if(need_verbose)then
1333 38 : call wrtout(ab_out,message,'COLL')
1334 38 : call wrtout(std_out,message,'COLL')
1335 : end if
1336 38 : ABI_FREE(weights1)
1337 :
1338 :
1339 114 : ABI_MALLOC(gf_values_iter,(4,ncoeff_to_select+1))
1340 658 : gf_values_iter(:,:) = zero
1341 : !Store initial gf_values as first value in gf_values_iter
1342 190 : gf_values_iter(:,1) = gf_values(:,1)
1343 38 : end subroutine initialize_gf
1344 :
1345 : !> Selects coefficients one by one based on goal function improvement
1346 : !>
1347 : !> This routine:
1348 : !> - Iteratively selects coefficients that minimize the goal function
1349 : !> - Handles coefficient selection in cycles with ncoeff_per_cycle terms per cycle
1350 : !> - Performs coefficient fitting and evaluation for each candidate
1351 : !> - Tracks selected coefficients and updates goal function values
1352 : !> - Applies convergence criteria based on goal function improvement
1353 : !> - Handles MPI parallelization for distributed coefficient selection
1354 : !> - Outputs progress information in verbose mode
1355 : !> - Maintains CSV log of goal function values if requested
1356 :
1357 38 : subroutine select_one_by_one()
1358 : integer :: ncycle_select ! number of cycles for the selection
1359 : integer :: ncoeff_this_cycle
1360 : integer :: ncoeff_selected
1361 : integer :: ic
1362 :
1363 38 : ncycle_select=ceiling(real(ncoeff_to_select)/real(ncoeff_per_cycle))
1364 38 : ncoeff_selected = ncoeff_preselected
1365 : !Option 1, we select the coefficients one by one
1366 38 : if(need_verbose.and.ncoeff_to_select > 0)then
1367 38 : write(message,'(a,3x,a,10x,a,14x,a,14x,a,14x,a)') " N","Selecting","MSDE","MSDFS","MSDF","MSDS"
1368 : !message = " N Selecting MSDE MSDFS MSDF MSDS"
1369 38 : call wrtout(ab_out,message,'COLL')
1370 38 : write(message,'(4x,a,6x,a,8x,a,8x,a,8x,a)') "Coefficient","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",&
1371 76 : & "(eV^2/A^2)"
1372 : !message = " Coefficient (eV^2/A^2) (eV^2/A^2) (eV^2/A^2) (eV^2/A^2)"
1373 38 : call wrtout(ab_out,message,'COLL')
1374 : end if
1375 :
1376 : ! Start fit process
1377 124 : do icycle_tmp = 1,ncycle_select
1378 86 : ncoeff_this_cycle= min(ncoeff_per_cycle, ncoeff_to_select+ncoeff_preselected-ncoeff_selected)
1379 86 : icycle=ncoeff_selected+1
1380 86 : list_coeffs_tmp(icycle)= icycle
1381 86 : if(need_verbose)then
1382 86 : write(message, '(4a,I0,a)')ch10,'--',ch10,' Try to find the best model with ',&
1383 172 : & ncoeff_selected+ncoeff_this_cycle,' coefficient'
1384 : !message = '--' // ch10 // ' Try to find the best model with '
1385 : !message = message // trim(itoa(ncoeff_selected+ncoeff_this_cycle)) // ' coefficient'
1386 :
1387 86 : if(ncoeff_selected+ncoeff_this_cycle> 1) message = trim(message)// 's'
1388 : !write(message, '(2a)') trim(message),'s'
1389 86 : if(nproc > 1) then
1390 48 : if(my_ncoeff>=1) then
1391 48 : write(message, '(2a,I0,a)')trim(message), ' (only the ',my_ncoeff,&
1392 96 : & ' first are printed for this CPU)'
1393 : ! message = trim(message) // ' (only the ' // trim(itoa(my_ncoeff))
1394 : ! message = trim(message) // ' first are printed for this CPU)'
1395 :
1396 : else
1397 0 : write(message, '(2a)')trim(message), ' (no coefficient treated by this CPU)'
1398 : !message = trim(message) // ' (no coefficient treated by this CPU)'
1399 : end if
1400 : end if
1401 86 : call wrtout(std_out,message,'COLL')
1402 118 : if(ncoeff_selected>0 .or. any(list_coeffs(:) > zero))then
1403 : !write(message, '(3a)') ' The coefficient numbers from the previous cycle are:',ch10,' ['
1404 75 : message = ' The coefficient numbers from the previous cycle are:' // ch10 // ' ['
1405 326 : do ii=1,ncoeff_selected
1406 326 : if(ii<ncoeff_selected)then
1407 176 : write(message, '(a,I0,a)') trim(message),list_coeffs(ii),','
1408 : !message = trim(message) // trim(itoa(list_coeffs(ii))) // ','
1409 : else
1410 75 : write(message, '(a,I0)') trim(message),list_coeffs(ii)
1411 : !message = trim(message) // trim(itoa(list_coeffs(ii)))
1412 : end if
1413 : end do
1414 : !write(message, '(3a)') trim(message),']',ch10
1415 75 : message = trim(message) // ']' // ch10
1416 75 : call wrtout(std_out,message,'COLL')
1417 : end if
1418 :
1419 86 : write(message,'(2x,a,12x,a,14x,a,13x,a,14x,a)') " Testing","MSDE","MSDFS","MSDF","MSDS"
1420 : !message = ' Testing MSDE MSDFS MSDF MSDS'
1421 86 : call wrtout(std_out,message,'COLL')
1422 86 : write(message,'(a,7x,a,8x,a,8x,a,8x,a)') " Coefficient","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",&
1423 172 : & "(eV^2/A^2)"
1424 : ! message = ' Coefficient (eV^2/A^2) (eV^2/A^2) (eV^2/A^2) (eV^2/A^2)'
1425 86 : call wrtout(std_out,message,'COLL')
1426 : end if!End if verbose
1427 :
1428 : !Print all GF VALUES in CSV if wanted
1429 : !Open *csv file for storing GF values of all cores for this iteration
1430 : !TODO: Should move inside the loop over all selected?
1431 86 : if(need_prt_GF_csv)then
1432 6 : write(filename,'(a,I1,a,I3.3,a,I3.3,a)') "GF_values_iatom",fit_iatom_in,"_proc",my_rank,"_iter",icycle,".csv"
1433 : !filename = "GF_values_iatom" // trim(itoa(fit_iatom_in)) // "_proc"
1434 : !filename = trim(filename) // trim(itoa(my_rank)) // "_iter" // trim(itoa(icycle)) // ".csv"
1435 :
1436 6 : unit_GF_val = get_unit()
1437 6 : if (open_file(filename,message,unit=unit_GF_val,form="formatted",&
1438 : & status="unknown",action="write") /= 0) then
1439 0 : ABI_ERROR(message)
1440 : end if
1441 : end if
1442 : ! Reset gf_values
1443 63086 : gf_values(:,:) = huge(0.0_dp)/5.0_dp
1444 :
1445 12686 : do icoeff=1,my_ncoeff
1446 12600 : if(isbanned(my_coeffindexes(icoeff)) .or. &
1447 : & isselected(my_coeffindexes(icoeff)) )then
1448 4495 : gf_values(:,icoeff) = huge(0.0_dp)/5.0_dp
1449 : cycle
1450 : endif
1451 87395 : if(any(list_coeffs==my_coeffindexes(icoeff)) .or. singular_coeffs(icoeff) == 1)then
1452 0 : gf_values(:,icoeff) = huge(0.0_dp)/5.0_dp
1453 : cycle
1454 : endif
1455 : !if(nbancoeff >= 1)then
1456 : !end if
1457 : !list_coeffs(icycle) = my_coeffindexes(icoeff)
1458 : !my_coeffindexes(my_coeffindexes(icoeff)) = .True.
1459 :
1460 11701 : if(need_initialize_data)then
1461 0 : my_icoeff = icoeff
1462 : else
1463 : ! Need to initialized the data for the fit for this coefficient
1464 11701 : my_icoeff = 1
1465 : call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
1466 : & fit_data%training_set%displacement,&
1467 : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
1468 : & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
1469 : & strten_coeffs,fit_data%training_set%ucvol,&
1470 11701 : my_coefflist(icoeff),1)
1471 : end if
1472 :
1473 : ! Fill the temporary arrays
1474 558385 : energy_coeffs_tmp(icycle,:) = energy_coeffs(my_icoeff,:)
1475 88027825 : fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:)
1476 3838489 : strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,my_icoeff)
1477 :
1478 :
1479 : ! call the fit process routine
1480 : ! This routine solves the linear system proposed
1481 : ! by C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
1482 65983 : do ic =1, icycle
1483 65983 : list_coeffs_tmp(ic) =ic
1484 : end do
1485 : call fit_polynomial_coeff_solve(coeff_values(1:icycle),fcart_coeffs_tmp,fit_data%fcart_diff,&
1486 : & energy_coeffs_tmp,fit_data%energy_diff,info,&
1487 : & list_coeffs_tmp(1:icycle),natom_sc,icycle,ncoeff_to_fit,ntime,&
1488 : & strten_coeffs_tmp,fit_data%strten_diff,&
1489 : & fit_data%training_set%sqomega,fit_on,int_fit_factors, &
1490 11701 : & weights=weights)
1491 :
1492 11701 : if(info==0)then
1493 17911 : if (need_positive.and.any(coeff_values(ncoeff_fix+1:icycle) < zero)) then
1494 : !write(message, '(a)') ' Negative value detected...'
1495 0 : message = 'Negative value detected...'
1496 0 : gf_values(:,icoeff) = zero
1497 0 : coeff_values = zero
1498 : else
1499 : call fit_polynomial_coeff_computeGF(coeff_values(1:icycle),energy_coeffs_tmp,&
1500 : & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
1501 : & gf_values(:,icoeff),list_coeffs_tmp(1:icycle),natom_sc,&
1502 : & icycle,ncoeff_to_fit,ntime,strten_coeffs_tmp,&
1503 11632 : & fit_data%strten_diff,fit_data%training_set%sqomega, weights=weights)
1504 11632 : write(j_char, '(i7)') my_coeffindexes(icoeff)
1505 : !j_char = itoa(my_coeffindexes(icoeff))
1506 :
1507 11632 : write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char), &
1508 : !& gf_values(4,icoeff)*factor*(1000*Ha_ev)**2 ,&
1509 11632 : & gf_values(4,icoeff)*HaBohr_eVAng**2, &
1510 11632 : & gf_values(1,icoeff)*HaBohr_eVAng**2, &
1511 11632 : & gf_values(2,icoeff)*HaBohr_eVAng**2, &
1512 23264 : & gf_values(3,icoeff)*HaBohr_eVAng**2
1513 : !message = trim(adjustl(j_char)) // " " // &
1514 : ! ftoa(gf_values(4,icoeff)*HaBohr_eVAng**2) // " " // &
1515 : ! ftoa(gf_values(1,icoeff)*HaBohr_eVAng**2) // " " // &
1516 : ! ftoa(gf_values(2,icoeff)*HaBohr_eVAng**2) // " " // &
1517 : ! ftoa(gf_values(3,icoeff)*HaBohr_eVAng**2)
1518 :
1519 :
1520 11632 : if(need_prt_GF_csv)then
1521 1635 : write(message2, '(I7.7,3a,ES18.10,a,ES18.10,a,ES18.10,a,ES18.10)') my_coeffindexes(icoeff),",",&
1522 1635 : & trim(my_coeffs(icoeff)%name),",", &
1523 1635 : & gf_values(4,icoeff)*HaBohr_eVAng**2,",", &
1524 1635 : & gf_values(1,icoeff)*HaBohr_eVAng**2,",", &
1525 1635 : & gf_values(2,icoeff)*HaBohr_eVAng**2,",", &
1526 3270 : & gf_values(3,icoeff)*HaBohr_eVAng**2
1527 : !message2 = itoa(my_coeffindexes(icoeff)) // "," // &
1528 : ! trim(my_coeffs(icoeff)%name) // "," // &
1529 : ! ftoa(gf_values(4,icoeff)*HaBohr_eVAng**2) // "," // &
1530 : ! ftoa(gf_values(1,icoeff)*HaBohr_eVAng**2) // "," // &
1531 : ! ftoa(gf_values(2,icoeff)*HaBohr_eVAng**2) // "," // &
1532 : ! ftoa(gf_values(3,icoeff)*HaBohr_eVAng**2)
1533 : end if
1534 : end if
1535 : else!In this case the matrix is singular.
1536 345 : gf_values(:,icoeff) = huge(0.0_dp)/5.0_dp
1537 69 : isbanned(my_coeffindexes(icoeff))=.True.
1538 69 : singular_coeffs(icoeff) = 1
1539 : !write(message, '(a)') ' The matrix is singular...'
1540 69 : message = 'The matrix is singular...'
1541 69 : if(need_prt_GF_csv)then
1542 30 : write(message2, '(I7.7,10a)') my_coeffindexes(icoeff),",", &
1543 30 : & trim(my_coeffs(icoeff)%name),",",&
1544 30 : & "None",",",&
1545 30 : & "None",",",&
1546 30 : & "None",",",&
1547 60 : & "None"
1548 : !message2 = itoa(my_coeffindexes(icoeff)) // "," // &
1549 : ! trim(my_coeffs(icoeff)%name) // "," // &
1550 : ! "None,None,None,None"
1551 : endif
1552 : end if
1553 11787 : if(need_verbose)then
1554 11701 : call wrtout(std_out,message,'COLL')
1555 11701 : if(need_prt_GF_csv)then
1556 1665 : call wrtout(unit_GF_val,message2,'PERS',do_flush=.TRUE.)
1557 : end if
1558 : endif
1559 : end do !icoeff=1,my_ncoeff
1560 86 : call xmpi_lor(isbanned, comm)
1561 :
1562 : !Close *csv file for GF values of this iteration
1563 86 : if(need_prt_GF_csv)close(unit_GF_val)
1564 :
1565 : ! find the best coeff on each CPU
1566 430 : mingf(:) = huge(0.0_dp)/5.0_dp
1567 86 : index_min = 0
1568 12686 : do icoeff=1,my_ncoeff
1569 : !if(gf_values(1,icoeff) > huge(0.0_dp)/5-1.0) cycle
1570 : !if(abs(gf_values(1,icoeff)) <tol16) cycle
1571 12600 : if(isbanned(my_coeffindexes(icoeff)) .or. isselected(my_coeffindexes(icoeff)) ) cycle
1572 81510 : if(sum(gf_values(2:4,icoeff),MASK=sel_on) < sum(mingf(2:4),MASK=sel_on))then
1573 1670 : mingf(:) = gf_values(:,icoeff)
1574 334 : index_min = my_coeffindexes(icoeff)
1575 : end if
1576 : end do
1577 :
1578 :
1579 : ! MPI GATHER THE BEST COEFF ON EACH CPU
1580 86 : if(nproc > 1)then
1581 48 : buffGF(1,1) = index_min
1582 240 : buffGF(2:5,1) = mingf(:)
1583 48 : call xmpi_barrier(comm)
1584 48 : call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr)
1585 240 : mingf(:) = huge(0.0_dp)/5.0_dp
1586 48 : index_min= 0
1587 208 : do icoeff=1,nproc
1588 160 : if(gf_mpi(2,icoeff) > huge(0.0_dp)/5-1) cycle
1589 : !if(abs(gf_mpi(2,icoeff)) < tol16) cycle
1590 1168 : if(sum(gf_mpi(3:5,icoeff),MASK=sel_on) < sum(mingf(2:4),MASK=sel_on))then
1591 310 : mingf(:) = gf_mpi(2:5,icoeff)
1592 62 : index_min = int(gf_mpi(1,icoeff))
1593 : end if
1594 : end do
1595 : end if
1596 :
1597 124 : BLOCK ! sort the coeff on each CPU
1598 172 : real(dp) :: mygf(my_ncoeff)
1599 172 : integer :: myorder(my_ncoeff), ntot
1600 86 : real(dp), allocatable :: allgf(:)
1601 86 : integer, allocatable :: allorder(:)
1602 12686 : do icoeff=1,my_ncoeff
1603 12600 : if(gf_values(1,icoeff) < zero) then
1604 0 : mygf(icoeff)=huge(0.0_dp)/5.0_dp
1605 12600 : else if(abs(gf_values(1,icoeff)) <tol16) then
1606 0 : mygf(icoeff)=huge(0.0_dp)/5.0_dp
1607 : else
1608 50400 : mygf(icoeff) = sum(gf_values(2:4,icoeff),MASK=sel_on)
1609 : end if
1610 12686 : myorder(icoeff) = my_coeffindexes(icoeff)
1611 : end do
1612 86 : call mpigatherv(mygf,myorder, my_ncoeff, allgf, allorder, ntot, comm, nproc)
1613 : BLOCK
1614 172 : real(dp) :: allgf_copy(size(allgf))
1615 172 : real(dp):: work((ntot+1)/2)
1616 172 : integer:: worder((ntot+1)/2)
1617 : integer :: i
1618 86 : logical :: ideg(ntot)
1619 20158 : ideg=.False.
1620 :
1621 20158 : allgf_copy(:)=allgf(:)
1622 86 : call MergeSort(allgf_copy, work, allorder, worder)
1623 20072 : do i=2, ntot
1624 20072 : if(abs(allgf(allorder(i))-allgf(allorder(i-1)))< 1e-16) then
1625 2592 : ideg(allorder(i))=.True.
1626 2592 : isbanned(allorder(i))=.True.
1627 : end if
1628 : end do
1629 20072 : do i=2, ntot
1630 20072 : if(ideg(allorder(i))) then
1631 : !print *, "new:", allgf(allorder(i-1)), allgf(allorder(i))
1632 2592 : allgf(allorder(i))=9D99
1633 : endif
1634 : end do
1635 86 : call MergeSort(allgf, work, allorder, worder)
1636 :
1637 86 : n_remaining = max(ceiling(n_remaining * remaining_rate), min(ncoeff_to_select*40, ncoeff_tot) )
1638 86 : do i=n_remaining+1, ncoeff_tot
1639 86 : isbanned(allorder(i))=.True.
1640 : end do
1641 :
1642 :
1643 : end BLOCK
1644 :
1645 : BLOCK ! add selected terms
1646 :
1647 : integer :: i
1648 : integer :: ind_select
1649 : integer :: nselected_this_cycle
1650 86 : nselected_this_cycle=0
1651 86 : i=0
1652 258 : do while(nselected_this_cycle<ncoeff_this_cycle)
1653 86 : i=i+1
1654 : !do while( isbanned(allorder(i)) .or. isselected(allorder(i)))
1655 : ! i=i+1
1656 : !enddo
1657 86 : index_min = allorder(i)
1658 : ! Check if there is still coefficient
1659 : !if(i>size(allorder))then
1660 : ! exit
1661 : !end if
1662 :
1663 86 : if(index_min==0) then
1664 : exit
1665 : ! TODO : check if this is necessary
1666 : !else if(is_duplicate_coeff(index_min)) then
1667 : ! cycle
1668 86 : else if ( isbanned(allorder(i)) .or. isselected(allorder(i))) then
1669 : cycle
1670 : else
1671 86 : ind_select=ncoeff_selected+1
1672 86 : list_coeffs(ind_select) = index_min
1673 86 : isselected(index_min)=.True.
1674 86 : nselected_this_cycle=nselected_this_cycle+1
1675 : end if
1676 : ! Check if this coeff is treat by this cpu and fill the
1677 : ! temporary array before broadcast
1678 86 : rank_to_send = 0
1679 6094 : do icoeff=1,my_ncoeff
1680 6094 : if((my_coeffindexes(icoeff)==list_coeffs(ind_select)))then
1681 54 : if(need_initialize_data)then
1682 0 : my_icoeff = icoeff
1683 : else
1684 : ! Need to initialized the data for the fit for this coefficient
1685 54 : my_icoeff = 1
1686 : call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
1687 : & fit_data%training_set%displacement,&
1688 : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
1689 : & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
1690 : & strten_coeffs,fit_data%training_set%ucvol,&
1691 54 : & my_coefflist(icoeff),1)
1692 : end if
1693 :
1694 2316 : energy_coeffs_tmp(ind_select,:) = energy_coeffs(my_icoeff,:)
1695 364236 : fcart_coeffs_tmp(:,:,ind_select,:) = fcart_coeffs(:,:,my_icoeff,:)
1696 15888 : strten_coeffs_tmp(:,:,ind_select) = strten_coeffs(:,:,my_icoeff)
1697 :
1698 54 : call polynomial_coeff_free(coeffs_tmp(ind_select))
1699 :
1700 :
1701 : call polynomial_coeff_init(coeff_values(ind_select),my_coeffs(icoeff)%nterm,&
1702 : & coeffs_tmp(ind_select),my_coeffs(icoeff)%terms,&
1703 : & my_coeffs(icoeff)%name,&
1704 54 : & check=.false.)
1705 54 : rank_to_send = my_rank
1706 54 : exit
1707 : end if
1708 : end do
1709 : ! Need to send the rank with the chosen coefficient
1710 : call xmpi_sum(rank_to_send, comm, ierr)
1711 : ! Boadcast the coefficient
1712 86 : call xmpi_bcast(energy_coeffs_tmp(ind_select,:), rank_to_send, comm, ierr)
1713 86 : call xmpi_bcast(fcart_coeffs_tmp(:,:,ind_select,:) , rank_to_send, comm, ierr)
1714 86 : call xmpi_bcast(strten_coeffs_tmp(:,:,ind_select), rank_to_send, comm, ierr)
1715 86 : call polynomial_coeff_broadcast(coeffs_tmp(ind_select), rank_to_send, comm)
1716 86 : if(need_verbose) then
1717 86 : write(message, '(a,I0,2a)' )' Selecting the coefficient number ',list_coeffs(ind_select),&
1718 172 : & ' ===> ',trim(coeffs_tmp(ind_select)%name)
1719 86 : call wrtout(std_out,message,'COLL')
1720 :
1721 86 : write (i_char, '(i3)') ind_select
1722 86 : write (j_char, '(i7)') list_coeffs(ind_select)
1723 86 : write(message, '(a,a,3x,a,3x,4ES18.10)') " ",adjustl(i_char),adjustl(j_char),&
1724 : !& mingf(4)* factor * (Ha_eV *1000)**2,&
1725 86 : & mingf(4)*HaBohr_eVAng**2,&
1726 86 : & mingf(1)*HaBohr_eVAng**2,&
1727 86 : & mingf(2)*HaBohr_eVAng**2,&
1728 172 : & mingf(3)*HaBohr_eVAng**2
1729 86 : call wrtout(ab_out,message,'COLL')
1730 : end if
1731 430 : gf_values_iter(:,icycle_tmp+1) = mingf(:)
1732 : !Store GF Values of this iteration
1733 : ! Check the stopping criterion
1734 86 : converge = .false.
1735 86 : if(tolGF > zero)then
1736 : check_value = (sum(gf_values_iter(2:4,icycle_tmp+1),MASK=sel_on) - &
1737 : & sum(gf_values_iter(2:4,icycle_tmp),MASK=sel_on)) &
1738 : & /(sum(gf_values_iter(2:4,icycle_tmp+1),MASK=sel_on) &
1739 0 : & - sum(gf_values_iter(2:4,1),MASK=sel_on))
1740 0 : if(check_value < tolGF)then
1741 0 : write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
1742 0 : & check_value ," < ",tolGF,&
1743 0 : & ' Goal Function is converged'
1744 0 : converge = .true.
1745 : end if
1746 : endif
1747 86 : if(tolMSDE > zero)then
1748 0 : if(abs(tolMSDE) > abs(mingf(4)* (Ha_eV *1000)**2 *factor))then
1749 0 : write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
1750 0 : & mingf(4)* (Ha_eV *1000)**2 * factor ," < ",tolMSDE,&
1751 0 : & ' for MSDE'
1752 0 : converge = .true.
1753 : end if
1754 : end if
1755 86 : if(tolMSDF > zero) then
1756 0 : if(abs(tolMSDF) > abs(mingf(2)*HaBohr_eVAng**2))then
1757 0 : write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
1758 0 : & mingf(2)*HaBohr_eVAng**2 ," < ",tolMSDF,&
1759 0 : & ' for MSDF'
1760 0 : converge = .true.
1761 : end if
1762 : end if
1763 86 : if(tolMSDS > zero) then
1764 0 : if(abs(tolMSDS) > abs(mingf(3)*HaBohr_eVAng**2))then
1765 0 : write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
1766 0 : & mingf(3)*HaBohr_eVAng**2 ," < ",tolMSDS,&
1767 0 : & ' for MSDS'
1768 0 : converge = .true.
1769 : end if
1770 : end if
1771 86 : if(tolMSDFS > zero)then
1772 0 : if(abs(tolMSDFS) > abs(mingf(1)*HaBohr_eVAng**2))then
1773 0 : write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",&
1774 0 : & mingf(1)*HaBohr_eVAng**2 ," < ",tolMSDFS,&
1775 0 : & ' for MSDFS'
1776 0 : converge = .true.
1777 : end if
1778 : end if
1779 86 : if(converge)then
1780 0 : call wrtout(ab_out,message,'COLL')
1781 0 : call wrtout(std_out,message,'COLL')
1782 0 : exit
1783 : else
1784 86 : if(any((/abs(tolMSDE),abs(tolMSDF),abs(tolMSDS),abs(tolMSDFS)/) > tol20) .and.&
1785 : & icycle_tmp == ncoeff_to_select)then
1786 0 : write(message,'(2a,I0,a)') ch10," WARNING: ",ncoeff_to_select,&
1787 0 : & " cycles was not enougth to converge the fit process"
1788 0 : call wrtout(ab_out,message,'COLL')
1789 0 : call wrtout(std_out,message,'COLL')
1790 : end if
1791 : end if
1792 86 : ncoeff_selected = ncoeff_selected + 1
1793 : !list_coeffs_tmp(ncoeff_selected) = ncoeff_selected
1794 : end do
1795 : end BLOCK
1796 86 : ABI_FREE(allgf)
1797 86 : ABI_FREE(allorder)
1798 : end BLOCK
1799 : end do !icycle_tmp=1,ncoeff_to_select
1800 38 : end subroutine select_one_by_one
1801 :
1802 :
1803 : !> Selects coefficients using Monte Carlo sampling approach
1804 : !>
1805 : !> This routine:
1806 : !> - Performs Monte Carlo sampling over coefficient combinations
1807 : !> - Runs for a fixed number of sweeps (default 10000)
1808 : !> - Randomly selects coefficients and evaluates goal function
1809 : !> - Keeps track of best coefficient combination found
1810 : !> - Handles MPI parallelization for distributed sampling
1811 : !> - Updates goal function values and coefficient statistics
1812 : !> - Transfers final model from best performing CPU
1813 : !> - Outputs progress and final results in verbose mode
1814 : !> - Optional validation of coefficient selection
1815 0 : subroutine select_with_monte_carlo()
1816 : ! Monte Carlo selection
1817 0 : nsweep = 10000
1818 : ! If no coefficient imposed in the inputs we reset the goal function
1819 0 : if (ncoeff_preselected == 0) then
1820 0 : gf_values(:,:) = zero
1821 0 : mingf(:) = 9D99
1822 : else
1823 0 : mingf = gf_values(:,1)
1824 : end if
1825 0 : call cpu_time(time)
1826 0 : call ZBQLINI(int(time*1000000/(my_rank+1)))
1827 0 : if(need_verbose)then
1828 0 : write(message,'(a,I0,a)') " Start Monte Carlo simulations on ", nproc," CPU"
1829 0 : if(nproc>1) write(message,'(2a)') trim(message)," (only print result of the master)"
1830 0 : call wrtout(std_out,message,'COLL')
1831 0 : call wrtout(ab_out,message,'COLL')
1832 0 : write(message,'(a,2x,a,9x,a,14x,a,13x,a,14x,a)') ch10," Iteration ","MSDE","MSDFS","MSDF","MSdS"
1833 0 : call wrtout(std_out,message,'COLL')
1834 0 : write(message,'(a,5x,a,8x,a,8x,a,8x,a)') " ","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",&
1835 0 : & "(eV^2/A^2)"
1836 0 : call wrtout(std_out,message,'COLL')
1837 :
1838 : end if
1839 :
1840 0 : do ii = 1,1!nyccle
1841 0 : do isweep =1,nsweep
1842 0 : write (j_char, '(i7)') isweep
1843 : !TEST_AM
1844 0 : icycle_tmp = int(ZBQLU01(zero)*(ncoeff_to_select+1-1))+1
1845 0 : icycle_tmp = ncoeff_to_select
1846 0 : do icycle=1,icycle_tmp
1847 0 : icoeff = int(ZBQLU01(zero)*(my_ncoeff))+1
1848 : ! icycle = int(ZBQLU01(zero)*(ncoeff_to_select))+1
1849 0 : list_coeffs_tmp2(icycle) = icoeff
1850 0 : list_coeffs_tmp(icycle)= icycle
1851 : ! Fill the temporary arrays
1852 0 : energy_coeffs_tmp(icycle,:) = energy_coeffs(icoeff,:)
1853 0 : fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,icoeff,:)
1854 0 : strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,icoeff)
1855 : end do
1856 : !TEST_AM
1857 :
1858 : ! call the fit process routine
1859 : ! This routine solves the linear system proposed by
1860 : ! C.Escorihuela-Sayalero see PRB95,094115(2017)
1861 : ! [[cite:Escorihuela-Sayalero2017]]
1862 : call fit_polynomial_coeff_solve(coeff_values(1:icycle_tmp)&
1863 : &,fcart_coeffs_tmp,fit_data%fcart_diff, energy_coeffs_tmp,fit_data&
1864 : &%energy_diff,info, list_coeffs_tmp(1:icycle_tmp),natom_sc,icycle_tmp&
1865 : &,ncoeff_to_fit, ntime,strten_coeffs_tmp,fit_data%strten_diff, fit_data&
1866 0 : &%training_set%sqomega,fit_on,int_fit_factors, weights=weights )
1867 0 : if(info==0)then
1868 : call fit_polynomial_coeff_computeGF(coeff_values(1:icycle_tmp),energy_coeffs_tmp,&
1869 : & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
1870 : & gf_values(:,1),list_coeffs_tmp(1:icycle_tmp),natom_sc,&
1871 : & icycle_tmp,ncoeff_to_fit,ntime,strten_coeffs_tmp,&
1872 0 : & fit_data%strten_diff,fit_data%training_set%sqomega, weights=weights)
1873 :
1874 : else!In this case the matrix is singular
1875 0 : gf_values(:,icoeff) = zero
1876 0 : singular_coeffs(icoeff) = 1
1877 : end if
1878 :
1879 0 : if(gf_values(1,1) > zero.and.abs(gf_values(1,1))>tol16.and.&
1880 0 : & gf_values(1,1) < mingf(1) ) then
1881 0 : mingf = gf_values(:,1)
1882 0 : list_coeffs(1:icycle_tmp) = list_coeffs_tmp2(1:icycle_tmp)
1883 0 : ncoeff_preselected = icycle_tmp
1884 :
1885 0 : write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char),&
1886 0 : & gf_values(4,1)* (1000*Ha_ev)**2 *factor,&
1887 0 : & gf_values(1,1)*HaBohr_eVAng**2,&
1888 0 : & gf_values(2,1)*HaBohr_eVAng**2,&
1889 0 : & gf_values(3,1)*HaBohr_eVAng**2
1890 0 : if(need_verbose) call wrtout(std_out,message,'COLL')
1891 : else
1892 0 : list_coeffs_tmp2(1:icycle_tmp) = list_coeffs(1:icycle_tmp)
1893 : end if
1894 : end do
1895 :
1896 0 : if(nproc > 1) then
1897 : !TEST_AM
1898 0 : do iproc=1,ncoeff_preselected
1899 0 : stat_coeff(list_coeffs(iproc)) = stat_coeff(list_coeffs(iproc)) + 1
1900 : end do
1901 : !TEST_AM
1902 :
1903 : ! Find the best model on all the CPUs
1904 0 : buffGF(1,1) = zero
1905 0 : buffGF(2:5,1) = mingf(:)
1906 0 : call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr)
1907 : ! find the best coeff
1908 0 : mingf(:) = huge(0.0_dp)/5.0_dp
1909 0 : index_min= 0
1910 0 : do iproc=1,nproc
1911 0 : if(gf_mpi(2,iproc) < zero) cycle
1912 0 : if(abs(gf_mpi(2,iproc)) <tol16) cycle
1913 0 : if(gf_mpi(2,iproc) < mingf(1) ) then
1914 0 : mingf(:) = gf_mpi(2:5,iproc)
1915 0 : index_min = int(gf_mpi(1,iproc))
1916 0 : rank_to_send = iproc-1
1917 : end if
1918 : end do
1919 0 : write(message, '(2a,I0)') ch10,' Best model found on the CPU: ', rank_to_send
1920 0 : call wrtout(std_out,message,'COLL')
1921 : end if
1922 : end do
1923 :
1924 : !TEST_AM
1925 : ! call xmpi_sum(stat_coeff, comm, ierr)
1926 : ! do ii=1,ncoeff_tot
1927 : ! write(100,*) ii,stat_coeff(ii)
1928 : ! end do
1929 : ! close(100)
1930 : !TEST_AM
1931 :
1932 : ! Transfer final model
1933 0 : if(nproc>1)then
1934 : call xmpi_bcast(ncoeff_preselected,rank_to_send,comm,ierr)
1935 0 : call xmpi_bcast(list_coeffs(1:ncoeff_preselected),rank_to_send,comm,ierr)
1936 : end if
1937 0 : do ii=1,ncoeff_preselected
1938 0 : icoeff = list_coeffs(ii)
1939 0 : list_coeffs_tmp(ii) = ii
1940 : ! Fill the temporary arrays
1941 0 : energy_coeffs_tmp(ii,:) = energy_coeffs(icoeff,:)
1942 0 : fcart_coeffs_tmp(:,:,ii,:) = fcart_coeffs(:,:,icoeff,:)
1943 0 : strten_coeffs_tmp(:,:,ii) = strten_coeffs(:,:,icoeff)
1944 0 : call polynomial_coeff_free(coeffs_tmp(ii))
1945 : call polynomial_coeff_init(one,my_coeffs(icoeff)%nterm,&
1946 : & coeffs_tmp(ii),my_coeffs(icoeff)%terms,&
1947 : & my_coeffs(icoeff)%name,&
1948 0 : & check=.false.)
1949 : end do
1950 :
1951 0 : end subroutine select_with_monte_carlo
1952 :
1953 : !> Fits all selected coefficients to find their optimal values
1954 : !>
1955 : !> This routine:
1956 : !> - Performs final fitting of all selected coefficients together
1957 : !> - Solves the complete linear system for all coefficients
1958 : !> - Updates coefficient values in the effective potential
1959 : !> - Computes final goal function values for the fit
1960 : !> - Handles weighted and unweighted goal function calculations
1961 : !> - Outputs detailed fitting results in verbose mode
1962 : !> - Updates the effective potential with fitted coefficients
1963 : !> - Optionally writes anharmonic terms to output files
1964 38 : subroutine fit_all_selected_coefficients()
1965 : !This routine solves the linear system proposed by
1966 : ! C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
1967 38 : real(dp), allocatable :: weights1(:)
1968 38 : if(ncoeff_to_fit > 0)then
1969 : block
1970 : integer :: ic
1971 198 : do ic =1, ncoeff_to_fit
1972 198 : list_coeffs_tmp(ic) =ic
1973 : end do
1974 : end block
1975 : call fit_polynomial_coeff_solve(coeff_values(1:ncoeff_to_fit),fcart_coeffs_tmp,fit_data%fcart_diff,&
1976 : & energy_coeffs_tmp,fit_data%energy_diff,info,&
1977 : & list_coeffs_tmp(1:ncoeff_to_fit),natom_sc,&
1978 : & ncoeff_to_fit,ncoeff_to_fit,ntime,strten_coeffs_tmp,&
1979 : & fit_data%strten_diff,fit_data%training_set%sqomega,fit_on,int_fit_factors, &
1980 38 : & weights=weights)
1981 :
1982 :
1983 38 : if(need_verbose) then
1984 38 : write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit process: '
1985 38 : call wrtout(ab_out,message,'COLL')
1986 38 : call wrtout(std_out,message,'COLL')
1987 : end if
1988 198 : do ii = 1,ncoeff_to_fit
1989 160 : if(list_coeffs(ii) ==0) cycle
1990 : ! Set the value of the coefficient
1991 160 : coeffs_tmp(ii)%coefficient = coeff_values(ii)
1992 198 : if(need_verbose) then
1993 160 : write(message, '(a,I0,a,ES19.10,2a)') " ",list_coeffs(ii)," =>",coeff_values(ii),&
1994 320 : & " ",trim(coeffs_tmp(ii)%name)
1995 160 : call wrtout(ab_out,message,'COLL')
1996 160 : call wrtout(std_out,message,'COLL')
1997 : end if
1998 : end do
1999 38 : if(fit_weight_T>0.0_dp) then
2000 : call fit_polynomial_coeff_computeGF(coeff_values(1:ncoeff_to_fit),energy_coeffs_tmp,&
2001 : & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
2002 : & gf_values(:,1),list_coeffs_tmp(1:ncoeff_to_fit),natom_sc,&
2003 : & ncoeff_to_fit,ncoeff_to_fit,ntime,strten_coeffs_tmp,&
2004 3 : & fit_data%strten_diff,fit_data%training_set%sqomega, weights=weights)
2005 :
2006 3 : if(need_verbose) then
2007 : ! Print the standard deviation after the fit
2008 3 : write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
2009 : !& ' Mean Standard Deviation values at the end of the fit process (meV^2/atm): ',ch10,&
2010 : !& ' Energy : ',&
2011 : !& gf_values(4,1)*(Ha_EV*1000)**2 *factor ,ch10,&
2012 3 : & 'Weighted goal function values at the end of the fit process (eV^2/A^2):',ch10,&
2013 3 : & ' Energy : ',&
2014 3 : & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,&
2015 3 : & ' Forces+Stresses : ',&
2016 3 : & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,&
2017 3 : & ' Forces : ',&
2018 3 : & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,&
2019 3 : & ' Stresses : ',&
2020 6 : & gf_values(3,1)*(HaBohr_eVAng)**2,ch10
2021 3 : call wrtout(ab_out,message,'COLL')
2022 3 : call wrtout(std_out,message,'COLL')
2023 : end if
2024 : endif
2025 :
2026 114 : ABI_MALLOC(weights1,(ntime))
2027 595 : weights1(:) = one
2028 : call fit_polynomial_coeff_computeGF(coeff_values(1:ncoeff_to_fit),energy_coeffs_tmp,&
2029 : & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,&
2030 : & gf_values(:,1),list_coeffs_tmp(1:ncoeff_to_fit),natom_sc,&
2031 : & ncoeff_to_fit,ncoeff_to_fit,ntime,strten_coeffs_tmp,&
2032 38 : & fit_data%strten_diff,fit_data%training_set%sqomega, weights=weights1)
2033 :
2034 38 : if(need_verbose) then
2035 : ! Print the standard deviation after the fit
2036 38 : write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
2037 : !& ' Mean Standard Deviation values at the end of the fit process (meV^2/atm): ',ch10,&
2038 : !& ' Energy : ',&
2039 : !& gf_values(4,1)*(Ha_EV*1000)**2 *factor ,ch10,&
2040 38 : & ' Goal function values at the end of the fit process (eV^2/A^2):',ch10,&
2041 38 : & ' Energy : ',&
2042 38 : & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,&
2043 38 : & ' Forces+Stresses : ',&
2044 38 : & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,&
2045 38 : & ' Forces : ',&
2046 38 : & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,&
2047 38 : & ' Stresses : ',&
2048 76 : & gf_values(3,1)*(HaBohr_eVAng)**2,ch10
2049 38 : call wrtout(ab_out,message,'COLL')
2050 38 : call wrtout(std_out,message,'COLL')
2051 : end if
2052 38 : ABI_FREE(weights1)
2053 :
2054 :
2055 : !Allocate output coeffs -> selected plus fixed ones
2056 38 : ncoeff_out = ncoeff_to_fit+eff_pot_fixed%anharmonics_terms%ncoeff
2057 :
2058 274 : ABI_MALLOC(coeffs_out,(ncoeff_out))
2059 198 : do ii = 1,ncoeff_out
2060 198 : if(ii <= eff_pot_fixed%anharmonics_terms%ncoeff)then
2061 : call polynomial_coeff_init(eff_pot_fixed%anharmonics_terms%coefficients(ii)%coefficient,&
2062 : & eff_pot_fixed%anharmonics_terms%coefficients(ii)%nterm,coeffs_out(ii),&
2063 : & eff_pot_fixed%anharmonics_terms%coefficients(ii)%terms,&
2064 : & eff_pot_fixed%anharmonics_terms%coefficients(ii)%name,&
2065 0 : & check = .TRUE.)
2066 : else
2067 160 : ia = ii - eff_pot_fixed%anharmonics_terms%ncoeff
2068 : call polynomial_coeff_init( &
2069 : coeffs_tmp(ia)%coefficient, & ! Coefficient value
2070 : coeffs_tmp(ia)%nterm, & ! Number of terms
2071 : coeffs_out(ii), & ! Output coefficient structure
2072 : coeffs_tmp(ia)%terms, & ! Terms array
2073 : coeffs_tmp(ia)%name, & ! Name of the coefficient
2074 : check = .true. & ! Validation flag
2075 160 : )
2076 :
2077 : !call polynomial_coeff_free(coeffs_tmp(ia))
2078 : endif
2079 :
2080 : enddo
2081 : !ABI_SFREE(coeffs_tmp)
2082 :
2083 :
2084 :
2085 :
2086 : ! Set the final set of coefficients into the eff_pot type
2087 38 : call effective_potential_setCoeffs(coeffs_out(:),eff_pot,ncoeff_out)
2088 :
2089 : ! If Wanted open the anharmonic_terms_file and write header
2090 38 : filename = "TRS_fit_diff"
2091 38 : if(need_prt_anh .and. ncoeff_model > 0 )then
2092 : call effective_potential_writeAnhHead(ncoeff_model,filename,&
2093 0 : & eff_pot%anharmonics_terms)
2094 38 : else if (need_prt_anh)then
2095 0 : write(message, '(6a,I3,3a)' )ch10,&
2096 0 : & ' --- !WARNING',ch10,&
2097 0 : & ' Printing of anharmonic terms has been asked,but',ch10,&
2098 0 : & ' there are',ncoeff_model,'anharmonic terms in the potential',ch10,&
2099 0 : & ' ---',ch10
2100 0 : call wrtout(ab_out,message,'COLL')
2101 0 : call wrtout(std_out,message,'COLL')
2102 : end if
2103 :
2104 : ! Calculate MSD values for final model
2105 38 : if(need_prt_files)call fit_polynomial_coeff_computeMSD(eff_pot,hist,gf_values(4,1),gf_values(2,1),gf_values(1,1),&
2106 : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
2107 18 : & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=filename, weights=weights)
2108 :
2109 :
2110 : INQUIRE(FILE='TRS_fit_diff_anharmonic_terms_energy.dat',OPENED=file_opened&
2111 38 : &,number=unit_anh)
2112 38 : if(file_opened) close(unit_anh)
2113 : else
2114 0 : ncoeff_out = 0
2115 0 : if(need_verbose) then
2116 0 : write(message, '(9a)' )ch10,&
2117 0 : & ' --- !WARNING',ch10,&
2118 0 : & ' The fit process does not provide possible terms.',ch10,&
2119 0 : & ' Please make sure that the terms set is correct',ch10,&
2120 0 : & ' ---',ch10
2121 0 : call wrtout(ab_out,message,'COLL')
2122 0 : call wrtout(std_out,message,'COLL')
2123 : end if
2124 : end if
2125 38 : end subroutine fit_all_selected_coefficients
2126 :
2127 : !> Selects a single coefficient to include in the model
2128 : !>
2129 : !> This routine:
2130 : !> - Takes a coefficient index and marks it as selected
2131 : !> - Updates global selection status arrays
2132 : !> - Updates coefficient counts and remaining coefficient pool
2133 : !> - Handles MPI synchronization of selection status
2134 : !> - Updates list of selected coefficients
2135 : !>
2136 : !> @param ind Index of coefficient to select
2137 : subroutine select_one_coeff(ind)
2138 : ! Note ind is my_coeffindexes(icoeff)
2139 : integer :: ind
2140 : if (.not. isselected(ind)) then
2141 : rank_to_send= my_rank
2142 : isselected(ind)=.True.
2143 : call xmpi_lor(isselected, comm)
2144 : end if
2145 : !ncoeff_selected=count(isselected)
2146 : if(isselected(ind)) then
2147 : ncoeff_selected = ncoeff_selected+1
2148 : n_remaining=n_remaining-1
2149 : call xmpi_bcast(n_remaining, rank_to_send, comm, ierr)
2150 : end if
2151 : list_coeffs(ncoeff_selected) = ind
2152 : end subroutine select_one_coeff
2153 :
2154 : !> Bans a coefficient from being selected in future iterations
2155 : !>
2156 : !> This routine:
2157 : !> - Marks a coefficient as banned in the global banned array
2158 : !> - Updates the count of remaining available coefficients
2159 : !> - Synchronizes banned status across MPI processes
2160 : !> - Only bans if coefficient wasn't already banned
2161 : !>
2162 : !> @param icoeff Index of coefficient to ban
2163 : subroutine ban_one_term(icoeff)
2164 : integer :: icoeff
2165 : if (.not. isbanned(my_coeffindexes(icoeff))) then
2166 : isbanned(my_coeffindexes(icoeff))=.True.
2167 : call xmpi_lor(isbanned, comm)
2168 : n_remaining=n_remaining-1
2169 : call xmpi_bcast(n_remaining, rank_to_send, comm, ierr)
2170 : end if
2171 : end subroutine ban_one_term
2172 :
2173 :
2174 : #endif
2175 :
2176 : end subroutine fit_polynomial_coeff_fit
2177 : !!***
2178 :
2179 :
2180 :
2181 : !!****f* m_fit_polynomial_coeff/fit_polynomial_coeff_getPositive
2182 : !!
2183 : !! NAME
2184 : !! fit_polynomial_coeff_getPositive
2185 : !!
2186 : !! FUNCTION
2187 : !! This routine fit a list of possible model.
2188 : !! Return in the isPositive array:
2189 : !! 0 if the model ii does not contain possive coefficients
2190 : !! 1 if the model ii contain possive coefficients
2191 : !!
2192 : !! INPUTS
2193 : !! eff_pot<type(effective_potential)> = effective potential
2194 : !! hist<type(abihist)> = The history of the MD (or snapshot of DFT
2195 : !! coeff_values(nmodel,ncoeff) = values of the coefficients for each model
2196 : !! isPositive(nmodel) = see description below
2197 : !! list_coeff(nmodel,ncoeff) = list of the models
2198 : !! ncoeff = number of coeff per model
2199 : !! nfixcoeff = will not test the nfixcoeff first coeffcients
2200 : !! nmodel = number of model
2201 : !! comm = MPI communicator
2202 : !! verbose = optional, flag for the verbose mode
2203 : !!
2204 : !! OUTPUT
2205 : !! eff_pot = effective potential datatype with new fitted coefficients
2206 : !!
2207 : !! SOURCE
2208 :
2209 0 : subroutine fit_polynomial_coeff_getPositive(eff_pot,hist,coeff_values,isPositive,list_coeff,ncoeff,&
2210 : & nfixcoeff,nmodel,comm,verbose, fit_weight_T)
2211 :
2212 : implicit none
2213 :
2214 : !Arguments ------------------------------------
2215 : !scalars
2216 : integer,intent(in) :: ncoeff,nfixcoeff,nmodel,comm
2217 : !arrays
2218 : integer,intent(in) :: list_coeff(nmodel,ncoeff)
2219 : integer,intent(out) :: isPositive(nmodel)
2220 : real(dp),intent(out) :: coeff_values(nmodel,ncoeff)
2221 : type(effective_potential_type),intent(inout) :: eff_pot
2222 : type(abihist),intent(inout) :: hist
2223 : logical,optional,intent(in) :: verbose
2224 : real(dp) :: fit_weight_T
2225 : !Local variables-------------------------------
2226 : !scalar
2227 : integer :: ierr,ii,info,imodel,my_nmodel,nmodel_alone
2228 : integer :: master,my_rank,ncoeff_tot,natom_sc,ncell
2229 : integer :: nproc,ntime
2230 : logical :: iam_master,need_verbose
2231 : !arrays
2232 0 : real(dp), allocatable :: weights(:)
2233 : integer :: sc_size(3)
2234 0 : integer,allocatable :: list_coeffs(:),my_modelindexes(:),my_modellist(:)
2235 0 : real(dp),allocatable :: energy_coeffs(:,:),fcart_coeffs(:,:,:,:), strten_coeffs(:,:,:)
2236 0 : type(polynomial_coeff_type),allocatable :: coeffs_in(:)
2237 0 : type(fit_data_type) :: fit_data
2238 : character(len=500) :: message
2239 : logical :: fit_on(3)
2240 : real(dp) :: fit_factors(3)
2241 : ! *************************************************************************
2242 :
2243 : !MPI variables
2244 0 : master = 0
2245 0 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
2246 0 : iam_master = (my_rank == master)
2247 :
2248 : !Initialisation of optional arguments
2249 0 : need_verbose = .TRUE.
2250 0 : if(present(verbose)) need_verbose = verbose
2251 :
2252 0 : fit_on(1) = .FALSE.
2253 0 : fit_on(2) = .TRUE.
2254 0 : fit_on(3) = .TRUE.
2255 :
2256 0 : fit_factors = (/1,1,1/)
2257 :
2258 : !Get the list of coefficients from the eff_pot
2259 0 : if(eff_pot%anharmonics_terms%ncoeff > 0)then
2260 : ! Copy the initial coefficients array
2261 0 : ncoeff_tot = eff_pot%anharmonics_terms%ncoeff
2262 0 : ABI_MALLOC(coeffs_in,(ncoeff_tot))
2263 0 : do ii=1,ncoeff_tot
2264 : call polynomial_coeff_init(eff_pot%anharmonics_terms%coefficients(ii)%coefficient,&
2265 : & eff_pot%anharmonics_terms%coefficients(ii)%nterm,&
2266 : & coeffs_in(ii),&
2267 : & eff_pot%anharmonics_terms%coefficients(ii)%terms,&
2268 : & eff_pot%anharmonics_terms%coefficients(ii)%name,&
2269 0 : & check=.false.)
2270 : end do
2271 : end if
2272 :
2273 : !Reset the output (we free the memory)
2274 0 : call effective_potential_freeCoeffs(eff_pot)
2275 :
2276 : !if the number of atoms in reference supercell into effpot is not corret,
2277 : !wrt to the number of atom in the hist, we set map the hist and set the good
2278 : !supercell
2279 0 : if (size(hist%xred,2) /= eff_pot%supercell%natom) then
2280 0 : call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose)
2281 : end if
2282 :
2283 : !Initialisation of constants
2284 0 : natom_sc = eff_pot%supercell%natom
2285 0 : ncell = eff_pot%supercell%ncells
2286 0 : ntime = hist%mxhist
2287 0 : do ii = 1, 3
2288 0 : sc_size(ii) = eff_pot%supercell%rlatt(ii,ii)
2289 : end do
2290 :
2291 : !Initialisation of arrays:
2292 0 : ABI_MALLOC(list_coeffs,(ncoeff_tot))
2293 0 : list_coeffs = 0
2294 0 : do ii = 1,ncoeff_tot
2295 0 : list_coeffs(ii) = ii
2296 : end do
2297 :
2298 : !Get the decomposition for each coefficients of the forces and stresses for
2299 : !each atoms and each step equations 11 & 12 of PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
2300 0 : if(need_verbose)then
2301 0 : write(message, '(a)' ) ' Initialisation of the fit process...'
2302 0 : call wrtout(std_out,message,'COLL')
2303 : end if
2304 : !Before the fit, compute constants with fit_data_compute.
2305 : !Conpute the strain of each configuration.
2306 : !Compute the displacmeent of each configuration.
2307 : !Compute the variation of the displacement due to strain of each configuration.
2308 : !Compute fixed forces and stresse and get the standard deviation.
2309 : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
2310 0 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose)
2311 :
2312 :
2313 0 : call get_weight_from_hist(hist, fit_weight_T, ntime, eff_pot%supercell%natom, weights, comm )
2314 :
2315 : !Get the decomposition for each coefficients of the forces,stresses and energy for
2316 : !each atoms and each step (see equations 11 & 12 of
2317 : ! PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]] + allocation
2318 0 : ABI_MALLOC(energy_coeffs,(ncoeff_tot,ntime))
2319 0 : ABI_MALLOC(fcart_coeffs,(3,natom_sc,ncoeff_tot,ntime))
2320 0 : ABI_MALLOC(strten_coeffs,(6,ntime,ncoeff_tot))
2321 :
2322 : call fit_polynomial_coeff_getFS(coeffs_in,fit_data%training_set%du_delta,&
2323 : & fit_data%training_set%displacement,&
2324 : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
2325 : & ncoeff_tot,ntime,sc_size,&
2326 : & fit_data%training_set%strain,strten_coeffs,&
2327 0 : & fit_data%training_set%ucvol,list_coeffs,ncoeff_tot)
2328 :
2329 :
2330 :
2331 : !set MPI, really basic stuff...
2332 0 : nmodel_alone = mod(nmodel,nproc)
2333 0 : my_nmodel = int(aint(real(nmodel,sp)/(nproc)))
2334 :
2335 0 : if(my_rank >= (nproc-nmodel_alone)) then
2336 0 : my_nmodel = my_nmodel + 1
2337 : end if
2338 :
2339 0 : ABI_MALLOC(my_modelindexes,(my_nmodel))
2340 0 : ABI_MALLOC(my_modellist,(my_nmodel))
2341 :
2342 :
2343 :
2344 : !2:compute the number of model and the list of the corresponding for each CPU.
2345 0 : do imodel=1,my_nmodel
2346 0 : if(my_rank >= (nproc-nmodel_alone))then
2347 : my_modelindexes(imodel)=(int(aint(real(nmodel,sp)/nproc)))*(my_rank)+&
2348 0 : & (my_rank - (nproc-nmodel_alone)) + imodel
2349 0 : my_modellist(imodel) = imodel
2350 : else
2351 0 : my_modelindexes(imodel)=(my_nmodel)*(my_rank) + imodel
2352 0 : my_modellist(imodel) = imodel
2353 : end if
2354 : end do
2355 :
2356 :
2357 :
2358 : !Start fit process
2359 0 : isPositive = 0
2360 0 : coeff_values = zero
2361 0 : do ii=1,my_nmodel
2362 0 : imodel = my_modelindexes(ii)
2363 : call fit_polynomial_coeff_solve(coeff_values(imodel,1:ncoeff),fcart_coeffs&
2364 : &,fit_data%fcart_diff, energy_coeffs,fit_data%energy_diff,info,&
2365 : & list_coeff(imodel,1:ncoeff),natom_sc,ncoeff, ncoeff_tot,ntime&
2366 : &,strten_coeffs,fit_data%strten_diff, fit_data%training_set%sqomega,fit_on&
2367 0 : &,fit_factors, weights=weights)
2368 :
2369 0 : if(info==0)then
2370 :
2371 0 : if (any(coeff_values(imodel,nfixcoeff+1:ncoeff) < zero))then
2372 : ! coeff_values(imodel,:) = zero
2373 0 : isPositive(imodel) = 0
2374 : else
2375 0 : isPositive(imodel) = 1
2376 : end if
2377 : end if
2378 : end do
2379 :
2380 0 : call xmpi_sum(isPositive, comm, ierr)
2381 0 : call xmpi_sum(coeff_values, comm, ierr)
2382 :
2383 :
2384 :
2385 : !Deallocation of arrays
2386 0 : do ii=1,ncoeff_tot
2387 0 : call polynomial_coeff_free(coeffs_in(ii))
2388 : end do
2389 0 : call fit_data_free(fit_data)
2390 :
2391 0 : ABI_FREE(coeffs_in)
2392 0 : ABI_FREE(energy_coeffs)
2393 0 : ABI_FREE(fcart_coeffs)
2394 0 : ABI_FREE(list_coeffs)
2395 0 : ABI_FREE(my_modelindexes)
2396 0 : ABI_FREE(my_modellist)
2397 0 : ABI_FREE(strten_coeffs)
2398 0 : ABI_FREE(weights)
2399 0 : end subroutine fit_polynomial_coeff_getPositive
2400 : !!***
2401 :
2402 : !!****f* m_fit_polynomial_coeff/fit_polynomial_coeff_getCoeffBound
2403 : !!
2404 : !! NAME
2405 : !! fit_polynomial_coeff_getCoeffBound
2406 : !!
2407 : !! FUNCTION
2408 : !! This routine fit a list of possible model.
2409 : !! Return in the isPositive array:
2410 : !!
2411 : !! INPUTS
2412 : !! NEED TO UPDATE
2413 : !! eff_pot<type(effective_potential)> = effective potential
2414 : !! hist<type(abihist)> = The history of the MD (or snapshot of DFT
2415 : !! comm = MPI communicator
2416 : !! verbose = optional, flag for the verbose mode
2417 : !!
2418 : !! OUTPUT
2419 : !!
2420 : !!
2421 : !! SOURCE
2422 :
2423 0 : subroutine fit_polynomial_coeff_getCoeffBound(eff_pot,coeffs_out,hist,ncoeff_bound,comm,verbose, max_nbody)
2424 :
2425 : implicit none
2426 :
2427 : !Arguments ------------------------------------
2428 : !scalars
2429 : integer,intent(in) :: comm, max_nbody(:)
2430 : integer,intent(out) :: ncoeff_bound
2431 : logical,optional,intent(in) :: verbose
2432 : !arrays
2433 : type(abihist),intent(inout) :: hist
2434 : type(effective_potential_type),target,intent(inout) :: eff_pot
2435 : type(polynomial_coeff_type),allocatable,intent(out) :: coeffs_out(:)
2436 : !Local variables-------------------------------
2437 : !scalar
2438 : integer :: counter,icoeff,icoeff_bound,idisp,istrain,ii
2439 : integer :: istart,iterm,ndisp,nstrain,nterm,ncoeff_model,ncoeff_in,ncoeff_max
2440 : real(dp):: weight
2441 : logical :: need_verbose
2442 : !arrays
2443 0 : integer,allocatable :: atindx(:,:),cells(:,:,:),direction(:)
2444 0 : integer,allocatable :: power_disps(:),power_strain(:),strain(:)
2445 0 : type(polynomial_term_type),dimension(:),allocatable :: terms
2446 0 : integer,allocatable :: odd_coeff(:),need_bound(:)
2447 0 : type(polynomial_coeff_type),pointer :: coeffs_in(:)
2448 0 : type(polynomial_coeff_type),allocatable :: coeffs_test(:)
2449 0 : character(len=5),allocatable :: symbols(:)
2450 : character(len=200):: name
2451 : character(len=500) :: msg
2452 : ! *************************************************************************
2453 :
2454 :
2455 : !set the inputs varaibles
2456 0 : ncoeff_model = eff_pot%anharmonics_terms%ncoeff
2457 0 : coeffs_in => eff_pot%anharmonics_terms%coefficients
2458 :
2459 : !Do check
2460 0 : if(ncoeff_model == 0)then
2461 0 : write(msg,'(a)')'ncoeff_model must be different to 0'
2462 0 : ABI_BUG(msg)
2463 : end if
2464 :
2465 : !Map the hist in order to be consistent with the supercell into reference_effective_potential
2466 0 : call effective_potential_file_mapHistToRef(eff_pot,hist,comm)
2467 :
2468 : !Initialisation of optional arguments
2469 0 : need_verbose = .TRUE.
2470 0 : if(present(verbose)) need_verbose = verbose
2471 :
2472 0 : write(msg, '(a)' ) ' Detection of the unbound coefficients'
2473 0 : if(need_verbose)call wrtout(std_out,msg,'COLL')
2474 :
2475 : !Allocation
2476 0 : ncoeff_max = 2 * ncoeff_model
2477 0 : ABI_MALLOC(odd_coeff,(ncoeff_max))
2478 0 : ABI_MALLOC(need_bound,(ncoeff_max))
2479 :
2480 0 : ABI_MALLOC(symbols,(eff_pot%crystal%natom))
2481 : call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,&
2482 0 : & symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
2483 :
2484 :
2485 0 : ABI_MALLOC(coeffs_test,(ncoeff_max))
2486 :
2487 0 : do icoeff=1,ncoeff_model
2488 : call polynomial_coeff_init(coeffs_in(icoeff)%coefficient,coeffs_in(icoeff)%nterm,&
2489 : & coeffs_test(icoeff),coeffs_in(icoeff)%terms,&
2490 : & coeffs_in(icoeff)%name, &
2491 0 : & check=.false.)
2492 : end do
2493 :
2494 : !array to know which coeff has to be bound
2495 0 : need_bound(:) = 1
2496 : counter = 0
2497 : ncoeff_in = ncoeff_model
2498 :
2499 0 : do while(.not.all(need_bound == 0).and.counter<1)
2500 : ! Get the coefficients with odd coefficient
2501 0 : odd_coeff = 0
2502 : if(counter>0) then
2503 : need_bound(1:ncoeff_in) = 0
2504 : icoeff_bound = ncoeff_in
2505 : else
2506 : icoeff_bound = 1
2507 : end if
2508 :
2509 0 : do icoeff=icoeff_bound,ncoeff_model
2510 0 : if(any(mod(coeffs_in(icoeff)%terms(1)%power_disp(:),2)/=0))then
2511 0 : odd_coeff(icoeff) = 1
2512 : end if
2513 0 : if(any(mod(coeffs_in(icoeff)%terms(1)%power_strain(:),2)/=0))then
2514 0 : odd_coeff(icoeff) = 1
2515 : end if
2516 0 : if(odd_coeff(icoeff) == 0 .and. coeffs_in(icoeff)%coefficient > zero) then
2517 0 : need_bound(icoeff) = 0
2518 : else
2519 0 : need_bound(icoeff) = 1
2520 : end if
2521 : end do
2522 0 : if(need_verbose)then
2523 0 : write(msg, '(a)' ) ' The following coefficients need to be bound:'
2524 0 : call wrtout(std_out,msg,'COLL')
2525 0 : do icoeff=1,ncoeff_model
2526 0 : if(need_bound(icoeff) == 1)then
2527 0 : write(msg, '(2a)' ) ' =>',trim(coeffs_in(icoeff)%name)
2528 0 : call wrtout(std_out,msg,'COLL')
2529 : end if
2530 : end do
2531 : end if
2532 :
2533 :
2534 0 : icoeff_bound = ncoeff_in + 1
2535 : if(counter==0)then
2536 : istart = 1
2537 : else
2538 : istart = ncoeff_in
2539 : end if
2540 :
2541 0 : ncoeff_bound = count(need_bound(istart:ncoeff_model)==1)
2542 :
2543 0 : do icoeff=istart,ncoeff_model
2544 0 : if(need_bound(icoeff)==1)then
2545 :
2546 0 : nterm = coeffs_in(icoeff)%nterm
2547 0 : ndisp = coeffs_in(icoeff)%terms(1)%ndisp
2548 0 : nstrain = coeffs_in(icoeff)%terms(1)%nstrain
2549 :
2550 0 : ABI_MALLOC(terms,(nterm))
2551 0 : ABI_MALLOC(atindx,(2,ndisp))
2552 0 : ABI_MALLOC(cells,(3,2,ndisp))
2553 0 : ABI_MALLOC(direction,(ndisp))
2554 0 : ABI_MALLOC(power_disps,(ndisp))
2555 0 : ABI_MALLOC(power_strain,(nstrain))
2556 0 : ABI_MALLOC(strain,(nstrain))
2557 :
2558 0 : do iterm=1,coeffs_in(icoeff)%nterm
2559 0 : atindx(:,:) = coeffs_in(icoeff)%terms(iterm)%atindx(:,:)
2560 0 : cells(:,:,:) = coeffs_in(icoeff)%terms(iterm)%cell(:,:,:)
2561 0 : direction(:) = coeffs_in(icoeff)%terms(iterm)%direction(:)
2562 0 : power_strain(:) = coeffs_in(icoeff)%terms(iterm)%power_strain(:)
2563 0 : power_disps(:) = coeffs_in(icoeff)%terms(iterm)%power_disp(:)
2564 0 : strain(:) = coeffs_in(icoeff)%terms(iterm)%strain(:)
2565 0 : weight = 1
2566 0 : do idisp=1,ndisp
2567 0 : if(mod(power_disps(idisp),2) /= 0) then
2568 0 : power_disps(idisp) = power_disps(idisp) + 1
2569 : else
2570 0 : power_disps(idisp) = power_disps(idisp) + 2
2571 : end if
2572 : end do
2573 0 : do istrain=1,nstrain
2574 0 : if(mod(power_strain(istrain),2) /= 0)then
2575 0 : power_strain(istrain) = power_strain(istrain) + 1
2576 : else
2577 0 : if(power_strain(istrain) < 4 ) power_strain(istrain) = power_strain(istrain) + 2
2578 : end if
2579 : end do
2580 :
2581 : call polynomial_term_init(atindx,cells,direction,ndisp,nstrain,terms(iterm),&
2582 0 : & power_disps,power_strain,strain,weight,check=.true.)
2583 : end do
2584 :
2585 0 : name = ""
2586 0 : call polynomial_coeff_init(one,nterm,coeffs_test(icoeff_bound),terms,name, check=.true.)
2587 0 : call polynomial_coeff_getName(name,coeffs_test(icoeff_bound),symbols,recompute=.TRUE.)
2588 0 : call polynomial_coeff_SetName(name,coeffs_test(icoeff_bound))
2589 :
2590 : ! Deallocate the terms
2591 0 : do iterm=1,nterm
2592 0 : call polynomial_term_free(terms(iterm))
2593 : end do
2594 0 : ABI_FREE(terms)
2595 0 : ABI_FREE(atindx)
2596 0 : ABI_FREE(cells)
2597 0 : ABI_FREE(direction)
2598 0 : ABI_FREE(power_disps)
2599 0 : ABI_FREE(power_strain)
2600 0 : ABI_FREE(strain)
2601 :
2602 0 : icoeff_bound = icoeff_bound + 1
2603 :
2604 : end if
2605 : end do
2606 :
2607 :
2608 0 : if(counter==0)ncoeff_model = ncoeff_model + ncoeff_bound
2609 : ! call effective_potential_setCoeffs(coeffs_test,eff_pot,ncoeff_model)
2610 : call fit_polynomial_coeff_fit(eff_pot,(/0/),(/0/),hist,0,(/0,0/),1,0,&
2611 : & -1,0,(/0/),1,comm,verbose=.true.,positive=.false., max_nbody=max_nbody, &
2612 0 : & drop_rate=0.0_dp, ncoeff_per_cycle=1, fit_weight_T=-0.1_dp)
2613 :
2614 0 : coeffs_in => eff_pot%anharmonics_terms%coefficients
2615 :
2616 0 : counter = counter + 1
2617 : end do
2618 :
2619 0 : ABI_MALLOC(coeffs_out,(ncoeff_bound))
2620 0 : do ii=1,ncoeff_bound
2621 0 : icoeff_bound = ncoeff_in + ii
2622 : call polynomial_coeff_init(one,coeffs_test(icoeff_bound)%nterm,coeffs_out(ii),&
2623 : & coeffs_test(icoeff_bound)%terms,coeffs_test(icoeff_bound)%name,&
2624 0 : & check=.true.)
2625 : end do
2626 : !Deallocation
2627 0 : do ii=ncoeff_model,ncoeff_max
2628 0 : call polynomial_coeff_free(coeffs_test(ii))
2629 : end do
2630 :
2631 : !Deallocation
2632 0 : do icoeff=1,ncoeff_max
2633 0 : call polynomial_coeff_free(coeffs_test(icoeff))
2634 : end do
2635 0 : ABI_FREE(coeffs_test)
2636 0 : ABI_FREE(odd_coeff)
2637 0 : ABI_FREE(need_bound)
2638 0 : ABI_FREE(symbols)
2639 :
2640 :
2641 0 : end subroutine fit_polynomial_coeff_getCoeffBound
2642 : !!***
2643 :
2644 :
2645 : !!****f* m_fit_polynomial_coeff/fit_polynomial_coeff_solve
2646 : !!
2647 : !! NAME
2648 : !! fit_polynomial_coeff_solve
2649 : !!
2650 : !! FUNCTION
2651 : !! Build and the solve the system to get the values of the coefficients
2652 : !! This routine solves the linear system proposed by
2653 : !! C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
2654 : !!
2655 : !! INPUTS
2656 : !! fcart_coeffs(3,natom_sc,ncoeff_max,ntime) = List of the values of the contribution to the
2657 : !! cartesian forces for all coefficients
2658 : !! for each direction and each time
2659 : !! fcart_diff(3,natom,ntime) = Difference of cartesian forces between DFT calculation and
2660 : !! fixed part of the model (more often harmonic part)
2661 : !! energy_coeffs(ncoeff,ntime) = value of the energy for each coefficient (Ha)
2662 : !! energy_diff(ntime) = Difference of energ ybetween DFT calculation and fixed part
2663 : !! of the model (more often harmonic part)
2664 : !! list_coeffs(ncoeff_fit) = List with the index of the coefficients used for this model
2665 : !! natom = Number of atoms
2666 : !! ncoeff_fit = Number of coeff for the fit (dimension of the system)
2667 : !! ncoeff_max = Maximum number of coeff in the list
2668 : !! ntime = Number of time (number of snapshot, number of md step...)
2669 : !! strten_coeffs(6,ntime,ncoeff_max) = List of the values of the contribution to the stress tensor
2670 : !! of the coefficients for each direction,time
2671 : !! strten_diff(6,natom) = Difference of stress tensor between DFT calculation and
2672 : !! fixed part of the model (more often harmonic part)
2673 : !! sqomega(ntime) = Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]
2674 : !! weights(ntime) = weight of each configuration
2675 : !! fit_on(3) = Flag to know if we fit the energy, forces and stresses
2676 : !! fit_factors(3) = Factors to apply to the energy, forces and stresses
2677 : !!
2678 : !! OUTPUT
2679 : !! coefficients(ncoeff_fit) = Values of the coefficients
2680 : !! info_out = 0: successful exit
2681 : !! < 0: if INFO = -i, the i-th argument had an illegal value
2682 : !! > 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is
2683 : !! exactly zero. The factorization has been completed,
2684 : !! but the factor U is exactly singular, so the solution
2685 : !! could not be computed. = 0: successful exit
2686 : !! information from the subroutine dsgesv in LAPACK
2687 : !!
2688 : !! SOURCE
2689 :
2690 11740 : subroutine fit_polynomial_coeff_solve(coefficients,fcart_coeffs,fcart_diff,energy_coeffs,energy_diff,&
2691 11740 : & info_out,list_coeffs,natom,ncoeff_fit,ncoeff_max,ntime,&
2692 11740 : & strten_coeffs,strten_diff,sqomega,fit_on,fit_factors, weights)
2693 :
2694 : implicit none
2695 :
2696 : !Arguments ------------------------------------
2697 : !scalars
2698 : integer,intent(in) :: natom,ncoeff_fit,ncoeff_max,ntime
2699 : integer,intent(out) :: info_out
2700 : !arrays
2701 : real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime)
2702 : real(dp),intent(in) :: energy_diff(ntime)
2703 : integer,intent(in) :: list_coeffs(ncoeff_fit)
2704 : real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime)
2705 : real(dp),intent(in) :: fcart_diff(3,natom,ntime)
2706 : real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max)
2707 : real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime)
2708 : real(dp),intent(out):: coefficients(ncoeff_fit)
2709 : real(dp),intent(in) :: fit_factors(3)
2710 : logical,intent(in) :: fit_on(3)
2711 : real(dp), intent(in) :: weights(ntime)
2712 :
2713 : !Local variables-------------------------------
2714 : !scalar
2715 : integer :: ia,itime,icoeff,jcoeff,icoeff_tmp,jcoeff_tmp,mu,LDA,LDB,LDX,LDAF,N,NRHS
2716 : real(dp):: efact,ffact,sfact,ftmpA,stmpA,ftmpB,stmpB,etmpA,etmpB,fmu,fnu,smu,snu,emu,enu
2717 : integer :: INFO
2718 : !integer :: ITER ! Only needed if DSGESV is used
2719 : real(dp):: RCOND
2720 23480 : real(dp):: fcart_coeffs_tmp(3,natom,ntime)
2721 11740 : real(dp),allocatable:: AF(:,:),BERR(:),FERR(:),WORK(:),C(:),R(:)
2722 11740 : integer,allocatable :: IPIV(:),IWORK(:),SWORK(:)
2723 : !arrays
2724 11740 : real(dp),allocatable :: A(:,:),B(:,:)
2725 : character(len=1) :: FACT,EQUED,TRANS
2726 : ! character(len=500) :: message
2727 : ! *************************************************************************
2728 :
2729 : !0-Set variables for the
2730 11740 : N = ncoeff_fit; NRHS = 1; LDA = ncoeff_fit; LDB = ncoeff_fit; LDX = ncoeff_fit
2731 11740 : LDAF = ncoeff_fit; RCOND = zero; INFO = 0; TRANS='N'; EQUED='N'; FACT='N'
2732 :
2733 : !Set the factors
2734 11740 : efact = fit_factors(1)*one/(ntime)
2735 11740 : ffact = fit_factors(2)*one/(3*natom*ntime)
2736 11740 : sfact = fit_factors(3)*one/(6*ntime)
2737 :
2738 : !0-Allocation
2739 46960 : ABI_MALLOC(A,(LDA,N))
2740 35220 : ABI_MALLOC(B,(LDB,NRHS))
2741 35220 : ABI_MALLOC(AF,(LDAF,N))
2742 35220 : ABI_MALLOC(IPIV,(N))
2743 35220 : ABI_MALLOC(R,(N))
2744 23480 : ABI_MALLOC(C,(N))
2745 11740 : ABI_MALLOC(FERR,(NRHS))
2746 11740 : ABI_MALLOC(BERR,(NRHS))
2747 35220 : ABI_MALLOC(WORK,(4*N))
2748 23480 : ABI_MALLOC(IWORK,(N))
2749 35220 : ABI_MALLOC(SWORK,(N*(N+NRHS)))
2750 465112 : A=zero; B=zero;
2751 453372 : AF = zero; IPIV = 1;
2752 120628 : R = one; C = one;
2753 35220 : FERR = zero; BERR = zero
2754 283960 : IWORK = 0; WORK = 0
2755 :
2756 : !1-Get forces and stresses from the model and fill A
2757 : ! Fill alsor B with the forces and stresses from
2758 : ! the DFT snapshot and the model
2759 : ! See equation 17 of PRB95 094115 (2017) [[cite:Escorihuela-Sayalero2017]]
2760 66184 : do icoeff=1,ncoeff_fit
2761 54444 : icoeff_tmp = list_coeffs(icoeff)
2762 530387800 : fcart_coeffs_tmp(:,:,:) = fcart_coeffs(:,:,icoeff_tmp,:)
2763 3348440 : ftmpA= zero; ftmpB = zero
2764 3348440 : stmpA= zero; stmpB = zero
2765 3348440 : etmpA= zero; etmpB = zero
2766 : ! loop over the configuration
2767 3348440 : do itime=1,ntime
2768 : ! Fill energy
2769 3293996 : emu = energy_coeffs(icoeff_tmp,itime)
2770 29312412 : do jcoeff=1,ncoeff_fit
2771 26018416 : jcoeff_tmp = list_coeffs(jcoeff)
2772 26018416 : enu = energy_coeffs(jcoeff_tmp,itime)
2773 29312412 : if(fit_on(3))then
2774 361889 : etmpA = emu*enu/(sqomega(itime)**(1.0/2.0))
2775 361889 : A(icoeff,jcoeff) = A(icoeff,jcoeff) + efact*etmpA* weights(itime)
2776 : endif
2777 : end do
2778 3293996 : if(fit_on(3))then
2779 73403 : etmpB = etmpB + energy_diff(itime)*emu/(sqomega(itime)**(1.0/2.0)) * weights(itime)!/ (sqomega(itime)**3)
2780 : else
2781 : etmpB = zero ! REMOVE THIS LINE TO TAKE INTO ACOUNT THE ENERGY
2782 : endif
2783 : ! Fill forces
2784 135053836 : do ia=1,natom
2785 530333356 : do mu=1,3
2786 395279520 : fmu = fcart_coeffs_tmp(mu,ia,itime)
2787 3517489440 : do jcoeff=1,ncoeff_fit
2788 3122209920 : jcoeff_tmp = list_coeffs(jcoeff)
2789 3122209920 : fnu = fcart_coeffs(mu,ia,jcoeff_tmp,itime)
2790 3122209920 : ftmpA = fmu*fnu
2791 3517489440 : if(fit_on(1))A(icoeff,jcoeff) = A(icoeff,jcoeff) + ffact*ftmpA* weights(itime)
2792 : end do
2793 527039360 : if(fit_on(1))then
2794 395279520 : ftmpB = ftmpB + fcart_diff(mu,ia,itime)*fmu* weights(itime)
2795 : else
2796 : ftmpB = zero
2797 : endif
2798 : end do !End loop dir
2799 : end do !End loop natom
2800 : ! Fill stresses
2801 23112416 : do mu=1,6
2802 19763976 : smu = strten_coeffs(mu,itime,icoeff_tmp)
2803 175874472 : do jcoeff=1,ncoeff_fit
2804 156110496 : jcoeff_tmp = list_coeffs(jcoeff)
2805 156110496 : snu = strten_coeffs(mu,itime,jcoeff_tmp)
2806 156110496 : stmpA = sqomega(itime)*smu*snu
2807 175874472 : if(fit_on(2))A(icoeff,jcoeff) = A(icoeff,jcoeff) + sfact*stmpA *weights(itime)
2808 : end do
2809 23057972 : if(fit_on(2))then
2810 19763976 : stmpB = stmpB + sqomega(itime)*strten_diff(mu,itime)*smu * weights(itime)
2811 : else
2812 : stmpB = zero
2813 : endif
2814 : end do !End loop stress dir
2815 : end do ! End loop time
2816 66184 : B(icoeff,1) = B(icoeff,1) + ffact*ftmpB + sfact*stmpB + efact*etmpB
2817 : end do ! End loop icoeff
2818 :
2819 : !2-Solve Ax=B
2820 :
2821 : !== Three different LAPACK routines are possible ==
2822 : !OLD VERSION..
2823 : ! call dgesvx(FACT,TRANS,N,NRHS,A,LDA,AF,LDAF,IPIV,EQUED,R,C,B,LDB,coefficients,LDX,&
2824 : !RCOND,FERR,BERR,WORK,IWORK,INFO)
2825 : !U is nonsingular
2826 : ! if (INFO==N+1) then
2827 : ! coefficients = zero
2828 : ! end if
2829 :
2830 : !VERSION PRIOR TO 20240817. However test multibinit multi_l_6_1 fails on the new reference machine (eos).
2831 : !call DSGESV(N,NRHS,A,LDA,IPIV,B,LDB,coefficients,LDX,WORK,SWORK,ITER,INFO)
2832 :
2833 : !Other routine, activated on 20240817, together with the change of reference machine
2834 : ! Check matrix conditioning
2835 : !if (any(isnan(A))) then
2836 : ! INFO = -2
2837 : ! print *, "NaN in A"
2838 : ! coefficients = zero
2839 : !endif
2840 :
2841 : !if (any(isnan(B))) then
2842 : ! INFO = -2
2843 : ! !print *, "NaN in B"
2844 : ! coefficients = zero
2845 : !endif
2846 : ! Solve the system
2847 : ! call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
2848 : !if (INFO == 0) then
2849 : ! Check solution for NaN or very large values
2850 : ! if (any(isnan(B(:,NRHS))) .or. any(abs(B(:,NRHS)) > 1.0e10_dp)) then
2851 : ! INFO = -3
2852 : ! coefficients = zero
2853 : ! !print *, "NaN or very large values in coefficients"
2854 : ! else
2855 : ! coefficients = B(:,NRHS)
2856 : ! print *, "INFO=", INFO
2857 : ! print *, "Good solution", coefficients
2858 : ! endif
2859 : !else
2860 : ! print*, "INFO=", INFO, "setting coefficients to zero"
2861 : ! coefficients = zero
2862 : !end if
2863 :
2864 : !==================================================
2865 :
2866 :
2867 11740 : call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
2868 66184 : coefficients = B(:,NRHS)
2869 :
2870 : !U is nonsingular
2871 11740 : if (INFO==N+2) then
2872 0 : coefficients = zero
2873 : end if
2874 :
2875 132299 : if(any(abs(coefficients)>1.0E10) .or. any(coefficients /= coefficients))then
2876 69 : INFO = 1
2877 218 : coefficients = zero
2878 : end if
2879 :
2880 11740 : info_out = INFO
2881 :
2882 11740 : ABI_FREE(AF)
2883 11740 : ABI_FREE(IPIV)
2884 11740 : ABI_FREE(R)
2885 11740 : ABI_FREE(C)
2886 11740 : ABI_FREE(FERR)
2887 11740 : ABI_FREE(BERR)
2888 11740 : ABI_FREE(WORK)
2889 11740 : ABI_FREE(IWORK)
2890 11740 : ABI_FREE(SWORK)
2891 11740 : ABI_FREE(A)
2892 11740 : ABI_FREE(B)
2893 :
2894 11740 : end subroutine fit_polynomial_coeff_solve
2895 : !!***
2896 :
2897 :
2898 :
2899 : !> Computes weights for history configurations based on forces and temperature
2900 : !>
2901 : !> This routine:
2902 : !> - Calculates average force magnitudes for each configuration
2903 : !> - Applies Boltzmann weighting based on input temperature
2904 : !> - Uses atomic forces to determine configuration importance
2905 : !> - Normalizes weights to sum to ntime
2906 : !> - Outputs detailed statistics about forces and weights
2907 : !> - Handles both temperature-dependent and uniform weighting
2908 : !>
2909 : !> @param hist Input history containing force data
2910 : !> @param temperature Temperature in K for Boltzmann weighting
2911 : !> @param ntime Number of timesteps in history
2912 : !> @param natom Number of atoms in system
2913 : !> @param weights Output array of weights for each configuration
2914 : !> @param comm MPI communicator
2915 40 : subroutine get_weight_from_hist(hist, temperature, ntime, natom, weights, comm)
2916 : type(abihist), intent(in) :: hist
2917 : !real(dp), intent(in) :: DFT_forces(3,natom,ntime)
2918 : real(dp), intent(in) :: temperature
2919 : integer, intent(in) :: ntime, natom
2920 : integer, intent(in) :: comm
2921 : real(dp), allocatable, intent(out) :: weights(:)
2922 : integer :: itime, iatom
2923 : ! Boltzmann constant in Ha/K
2924 : real(dp), parameter :: kb = 3.166815d-6
2925 80 : real(dp) :: average_forces(ntime)
2926 40 : character(len=ntime*100) :: msg
2927 : integer :: master, my_rank
2928 : logical :: iam_master
2929 :
2930 40 : master = 0
2931 40 : my_rank = xmpi_comm_rank(comm)
2932 40 : iam_master = (my_rank == master)
2933 120 : ABI_MALLOC(weights, (ntime))
2934 799 : weights=1.0_dp
2935 :
2936 : ! compute average forces for each time step
2937 799 : average_forces = 0.0_dp
2938 799 : do itime=1,ntime
2939 31119 : do iatom=1,natom
2940 122199 : average_forces(itime) = average_forces(itime) + norm2(hist%fcart(:,iatom,itime))
2941 : end do
2942 799 : average_forces(itime) = average_forces(itime)/natom
2943 : end do
2944 :
2945 40 : if (temperature > 0.0_dp) then
2946 : ! compute the average norm of force for each atom
2947 : ! compute the weight
2948 36 : weights = 0.0_dp
2949 36 : do itime=1,ntime
2950 33 : weights(itime) = exp(-average_forces(itime)/(kb*temperature))
2951 36 : if (weights(itime)< 0.001) then
2952 21 : weights(itime) = 0.001
2953 : endif
2954 : end do
2955 : ! normalize the weight
2956 69 : weights = weights/sum(weights)* ntime
2957 : end if
2958 : ! For each time, print the average forces, the energy, the weight,
2959 : ! the maximum of axial stress, and maximum of shear stress
2960 : ! first write the header of the table
2961 : ! write a separator
2962 :
2963 40 : if(iam_master)then
2964 24 : write(msg, '(a)') '-----------------------------------------------------------------------------------'
2965 24 : write(msg, '(a)') 'Time Average_forces Etot-min(Etot) Weight Max(axial_stress) Max(shear_stress)'
2966 24 : call wrtout(std_out,msg,'COLL')
2967 : ! then write the data
2968 607 : do itime=1,ntime
2969 583 : write(msg, '(i5, 3e15.6, 2e15.6)') itime, average_forces(itime), &
2970 46609 : & hist%etot(itime)-minval(hist%etot(:)), weights(itime), &
2971 6413 : & maxval(abs(hist%strten(1:3,itime))), maxval(abs(hist%strten(4:6,itime)))
2972 607 : call wrtout(std_out,msg,'COLL')
2973 : end do
2974 24 : write(msg, '(a)') '-----------------------------------------------------------------------------------'
2975 : endif
2976 :
2977 :
2978 40 : end subroutine get_weight_from_hist
2979 :
2980 :
2981 :
2982 :
2983 : !> Computes goal function values comparing model predictions to reference data
2984 : !>
2985 : !> This routine calculates:
2986 : !> - Forces + stresses combined goal function (gf_value(1))
2987 : !> - Forces only goal function (gf_value(2))
2988 : !> - Stresses only goal function (gf_value(3))
2989 : !> - Energy goal function (gf_value(4))
2990 : !>
2991 : !> The calculations:
2992 : !> - Use supplied coefficient values to evaluate model
2993 : !> - Compare to reference energy/force/stress differences
2994 : !> - Apply configuration weights if provided
2995 : !> - Handle volume/supercell normalization
2996 : !> - Support parallel evaluation across processes
2997 : !>
2998 : !> @param coefficients Model coefficient values
2999 : !> @param energy_coeffs Energy terms for each coefficient
3000 : !> @param energy_diff Reference energy differences
3001 : !> @param fcart_coeffs Force terms for each coefficient
3002 : !> @param fcart_diff Reference force differences
3003 : !> @param gf_value Output goal function values
3004 : !> @param list_coeffs List of coefficients to include
3005 : !> @param natom Number of atoms
3006 : !> @param ncoeff_fit Number of coefficients being fit
3007 : !> @param ncoeff_max Maximum number of coefficients
3008 : !> @param ntime Number of configurations
3009 : !> @param strten_coeffs Stress terms for each coefficient
3010 : !> @param strten_diff Reference stress differences
3011 : !> @param sqomega Configuration weighting factors
3012 : !> @param weights Optional weights for each configuration
3013 11714 : subroutine fit_polynomial_coeff_computeGF(coefficients,energy_coeffs,energy_diff,&
3014 11714 : & fcart_coeffs,fcart_diff,gf_value,list_coeffs,&
3015 11714 : & natom,ncoeff_fit,ncoeff_max,ntime,strten_coeffs,&
3016 11714 : & strten_diff,sqomega, weights)
3017 :
3018 : implicit none
3019 :
3020 : !Arguments ------------------------------------
3021 : !scalars
3022 : integer,intent(in) :: natom,ncoeff_fit,ncoeff_max,ntime
3023 : !arrays
3024 : integer,intent(in) :: list_coeffs(ncoeff_fit)
3025 : real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime)
3026 : real(dp),intent(in) :: energy_diff(ntime)
3027 : real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime)
3028 : real(dp),intent(in) :: fcart_diff(3,natom,ntime)
3029 : real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max)
3030 : real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime)
3031 : real(dp),intent(in) :: coefficients(ncoeff_fit)
3032 : real(dp), intent(in) :: weights(ntime)
3033 : real(dp),intent(out) :: gf_value(4)
3034 : !Local variables-------------------------------
3035 : !scalar
3036 : integer :: ia,icoeff,icoeff_tmp,itime,mu
3037 : real(dp):: etmp,emu,fmu,ftmp,smu,stmp
3038 : real(dp) :: ffact,sfact,efact
3039 : !arrays
3040 : ! *************************************************************************
3041 :
3042 : !1-Compute the value of the goal function
3043 : ! see equation 9 of PRB 95 094115(2017) [[cite:Escorihuela-Sayalero2017]]
3044 : gf_value = zero
3045 11714 : etmp = zero
3046 11714 : ftmp = zero
3047 11714 : stmp = zero
3048 :
3049 : !Compute factors
3050 11714 : ffact = one/(3*natom*ntime)
3051 11714 : sfact = one/(6*ntime)
3052 11714 : efact = one/(ntime)
3053 :
3054 : ! loop over the configuration
3055 558819 : do itime=1,ntime
3056 : ! Fill energy
3057 : emu = zero
3058 3839392 : do icoeff=1,ncoeff_fit
3059 3292287 : icoeff_tmp = list_coeffs(icoeff)
3060 3839392 : emu = emu + coefficients(icoeff)*energy_coeffs(icoeff_tmp,itime)
3061 : end do
3062 : ! uncomment the next line to be consistent with the definition of the goal function
3063 547105 : etmp = etmp + (energy_diff(itime)-emu)**2* weights(itime)/(sqomega(itime)**(1.0/2.0))
3064 : ! uncomment the next get a measure in Ha instead of Ha^2
3065 : ! etmp = etmp + abs(energy_diff(itime)-emu)
3066 : ! Fill forces
3067 22431305 : do ia=1,natom
3068 88083905 : do mu=1,3
3069 : fmu = zero
3070 460727040 : do icoeff=1,ncoeff_fit
3071 395074440 : icoeff_tmp = list_coeffs(icoeff)
3072 460727040 : fmu = fmu + coefficients(icoeff)*fcart_coeffs(mu,ia,icoeff_tmp,itime)
3073 : end do
3074 87536800 : ftmp = ftmp + (fcart_diff(mu,ia,itime)-fmu)**2* weights(itime)
3075 : end do !End loop dir
3076 : end do !End loop natom
3077 3841449 : do mu=1,6
3078 : smu = zero
3079 23036352 : do icoeff=1,ncoeff_fit
3080 19753722 : icoeff_tmp = list_coeffs(icoeff)
3081 23036352 : smu = smu + coefficients(icoeff)*strten_coeffs(mu,itime,icoeff_tmp)
3082 : end do
3083 3829735 : stmp = stmp + sqomega(itime)* weights(itime) *(strten_diff(mu,itime)-smu)**2
3084 : end do !End loop stress dir
3085 : end do ! End loop time
3086 :
3087 11714 : gf_value(1) = ffact*ftmp + sfact*stmp !+ efact*etmp !Stresses + Forces
3088 11714 : gf_value(2) = ffact*ftmp ! only Forces
3089 11714 : gf_value(3) = sfact*stmp ! only Stresses
3090 11714 : gf_value(4) = efact*etmp !abs(Energy)
3091 :
3092 11714 : end subroutine fit_polynomial_coeff_computeGF
3093 : !!***
3094 :
3095 :
3096 : !> Computes force and stress contributions for all polynomial coefficients
3097 : !>
3098 : !> This routine:
3099 : !> - Evaluates energy, force and stress terms for each coefficient
3100 : !> - Handles displacement-displacement and displacement-strain couplings
3101 : !> - Applies periodic boundary conditions for supercell calculations
3102 : !> - Computes derivatives needed for force and stress evaluations
3103 : !> - Includes strain-dependent volume corrections
3104 : !> - Supports parallel evaluation across processes
3105 : !>
3106 : !> Based on equations 11 & 12 in PRB95,094115 (2017) [[cite:Escorihuela-Sayalero2017]]
3107 : !>
3108 : !> @param coefficients Input polynomial coefficients
3109 : !> @param du_delta Changes in displacements due to strain
3110 : !> @param displacement Atomic displacements
3111 : !> @param energy_out Output energy terms for each coefficient
3112 : !> @param fcart_out Output force terms for each coefficient
3113 : !> @param natom_sc Number of atoms in supercell
3114 : !> @param natom_uc Number of atoms in unit cell
3115 : !> @param ncoeff_max Maximum number of coefficients
3116 : !> @param ntime Number of configurations
3117 : !> @param sc_size Supercell dimensions
3118 : !> @param strain Strain values
3119 : !> @param strten_out Output stress terms for each coefficient
3120 : !> @param ucvol Unit cell volumes
3121 : !> @param coeffs List of coefficients to evaluate
3122 : !> @param ncoeff Number of coefficients to evaluate
3123 11798 : subroutine fit_polynomial_coeff_getFS(coefficients,du_delta,displacement,energy_out,fcart_out,&
3124 11798 : & natom_sc,natom_uc,ncoeff_max,ntime,sc_size,strain,strten_out,&
3125 11798 : & ucvol,coeffs,ncoeff)
3126 :
3127 : implicit none
3128 :
3129 : !Arguments ------------------------------------
3130 : !scalars
3131 : integer,intent(in) :: natom_sc,natom_uc,ncoeff_max,ntime
3132 : integer,intent(in) :: ncoeff
3133 : !arrays
3134 : integer,intent(in) :: sc_size(3)
3135 : integer,intent(in) :: coeffs(ncoeff_max)
3136 : real(dp),intent(in) :: du_delta(6,3,natom_sc,ntime)
3137 : real(dp),intent(in) :: displacement(3,natom_sc,ntime)
3138 : real(dp),intent(in) :: strain(6,ntime),ucvol(ntime)
3139 : real(dp),intent(out):: energy_out(ncoeff,ntime)
3140 : real(dp),intent(out) :: fcart_out(3,natom_sc,ncoeff,ntime)
3141 : real(dp),intent(out) :: strten_out(6,ntime,ncoeff)
3142 : type(polynomial_coeff_type), intent(in) :: coefficients(ncoeff_max)
3143 : !Local variables-------------------------------
3144 : !scalar
3145 : integer :: i1,i2,i3,ia1,ia2,ib1,ib2,ii,icell,icoeff,icoeff_tmp
3146 : integer :: idir1,idir2,idisp1,idisp2,idisp1_strain,idisp2_strain
3147 : integer :: iterm,itime,ndisp,ndisp_tot,nstrain,power_disp,power_strain
3148 : real(dp):: disp1,disp2,tmp1,tmp2,tmp3,weight
3149 : !arrays
3150 : integer :: cell_atoma1(3),cell_atoma2(3)
3151 : integer :: cell_atomb1(3),cell_atomb2(3)
3152 :
3153 : ! *************************************************************************
3154 :
3155 :
3156 : !1-Get forces and stresses from the model
3157 : ! Initialisation of variables
3158 89048517 : fcart_out(:,:,:,:) = zero
3159 3870867 : strten_out(:,:,:) = zero
3160 1110917 : energy_out(:,:) = zero
3161 11798 : icell = 0; ib1=0; ia1=0
3162 35394 : do i1=1,sc_size(1)
3163 82586 : do i2=1,sc_size(2)
3164 165172 : do i3=1,sc_size(3)
3165 94384 : ii = icell*natom_uc
3166 94384 : icell = icell + 1
3167 : ! Loop over configurations
3168 4537648 : do itime=1,ntime
3169 : ! Loop over coefficients
3170 8887336 : do icoeff_tmp=1,ncoeff
3171 4396880 : icoeff = coeffs(icoeff_tmp)
3172 : ! Loop over terms of this coefficient
3173 129038056 : do iterm=1,coefficients(icoeff)%nterm
3174 120245104 : ndisp = coefficients(icoeff)%terms(iterm)%ndisp
3175 120245104 : nstrain = coefficients(icoeff)%terms(iterm)%nstrain
3176 120245104 : ndisp_tot = ndisp + nstrain
3177 : ! Set the weight of this term
3178 120245104 : weight =coefficients(icoeff)%terms(iterm)%weight
3179 120245104 : tmp1 = one
3180 : ! Loop over displacement and strain
3181 527428432 : do idisp1=1,ndisp_tot
3182 :
3183 : ! Set to one the acculation of forces and strain
3184 407183328 : tmp2 = one
3185 407183328 : tmp3 = one
3186 : ! Strain case idir => -6, -5, -4, -3, -2 or -1
3187 407183328 : if (idisp1 > ndisp)then
3188 48745920 : idisp1_strain = idisp1 - ndisp
3189 48745920 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain)
3190 : ! Get the direction of the displacement or strain
3191 48745920 : idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain)
3192 48745920 : if(abs(strain(idir1,itime)) > tol10)then
3193 : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
3194 29382624 : tmp1 = tmp1 * (strain(idir1,itime))**power_strain
3195 29382624 : if(power_strain > 1) then
3196 : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
3197 45360 : tmp3 = tmp3 * power_strain*(strain(idir1,itime))**(power_strain-1)
3198 : end if
3199 : else
3200 19363296 : tmp1 = zero
3201 19363296 : if(power_strain > 1) then
3202 2160 : tmp3 = zero
3203 : end if
3204 : end if
3205 : else
3206 : ! Set the power_disp of the displacement:
3207 358437408 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1)
3208 : ! Get the direction of the displacement or strain
3209 358437408 : idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1)
3210 : ! Displacement case idir = 1, 2 or 3
3211 : ! indexes of the cell of the atom a
3212 1433749632 : cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1)
3213 358437408 : if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then
3214 : ! if the cell is not 0 0 0 we apply PBC:
3215 0 : cell_atoma1(1) = i1 + cell_atoma1(1)
3216 0 : cell_atoma1(2) = i2 + cell_atoma1(2)
3217 0 : cell_atoma1(3) = i3 + cell_atoma1(3)
3218 0 : call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3))
3219 : ! index of the first atom (position in the supercell if the cell is not 0 0 0)
3220 : ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3221 : & (cell_atoma1(2)-1)*sc_size(3)*natom_uc+&
3222 : & (cell_atoma1(3)-1)*natom_uc+&
3223 0 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
3224 : else
3225 : ! index of the first atom (position in the supercell if the cell is 0 0 0)
3226 358437408 : ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
3227 : end if
3228 :
3229 : ! indexes of the cell of the atom b (with PBC) same as ia1
3230 1433749632 : cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1)
3231 358437408 : if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then
3232 195199584 : cell_atomb1(1) = i1 + cell_atomb1(1)
3233 195199584 : cell_atomb1(2) = i2 + cell_atomb1(2)
3234 195199584 : cell_atomb1(3) = i3 + cell_atomb1(3)
3235 195199584 : call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3))
3236 :
3237 : ! index of the second atom in the (position in the supercell if the cell is not 0 0 0)
3238 : ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3239 : & (cell_atomb1(2)-1)*sc_size(3)*natom_uc+&
3240 : & (cell_atomb1(3)-1)*natom_uc+&
3241 195199584 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
3242 : else
3243 : ! index of the first atom (position in the supercell if the cell is 0 0 0)
3244 163237824 : ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
3245 : end if
3246 :
3247 : ! Get the displacement for the both atoms
3248 358437408 : disp1 = displacement(idir1,ia1,itime)
3249 358437408 : disp2 = displacement(idir1,ib1,itime)
3250 :
3251 358437408 : if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then
3252 : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
3253 351505920 : tmp1 = tmp1 * (disp1-disp2)**power_disp
3254 351505920 : if(power_disp > 1) then
3255 : ! Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...))
3256 32850240 : tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1)
3257 : end if
3258 : else
3259 6931488 : tmp1 = zero
3260 6931488 : if(power_disp > 1) then
3261 368736 : tmp2 = zero
3262 : end if
3263 : end if
3264 : end if
3265 :
3266 1831660032 : do idisp2=1,ndisp_tot
3267 1831660032 : if(idisp2 /= idisp1) then
3268 :
3269 : ! Strain case
3270 1017293376 : if (idisp2 > ndisp)then
3271 126657792 : idisp2_strain = idisp2 - ndisp
3272 126657792 : idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain)
3273 : ! Set the power_strain of the strain:
3274 126657792 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain)
3275 : ! Accumulate energy forces
3276 126657792 : tmp2 = tmp2 * (strain(idir2,itime))**power_strain
3277 : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
3278 126657792 : tmp3 = tmp3 * (strain(idir2,itime))**power_strain
3279 : ! Atomic displacement case
3280 : else
3281 : ! Set the power_disp of the displacement:
3282 890635584 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2)
3283 : ! Set the direction of the displacement:
3284 890635584 : idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2)
3285 :
3286 3562542336 : cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2)
3287 890635584 : if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then
3288 0 : cell_atoma2(1) = i1 + cell_atoma2(1)
3289 0 : cell_atoma2(2) = i2 + cell_atoma2(2)
3290 0 : cell_atoma2(3) = i3 + cell_atoma2(3)
3291 0 : call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3))
3292 : ! index of the first atom (position in the supercell and direction)
3293 : ! if the cell of the atom a is not 0 0 0 (may happen)
3294 : ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3295 : & (cell_atoma2(2)-1)*sc_size(3)*natom_uc+&
3296 : & (cell_atoma2(3)-1)*natom_uc+&
3297 0 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
3298 : else
3299 : ! index of the first atom (position in the supercell and direction)
3300 890635584 : ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
3301 : end if
3302 :
3303 3562542336 : cell_atomb2 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2)
3304 :
3305 890635584 : if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then
3306 : ! indexes of the cell2 (with PBC)
3307 475655760 : cell_atomb2(1) = i1 + cell_atomb2(1)
3308 475655760 : cell_atomb2(2) = i2 + cell_atomb2(2)
3309 475655760 : cell_atomb2(3) = i3 + cell_atomb2(3)
3310 475655760 : call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3))
3311 :
3312 : ! index of the second atom in the (position in the supercell)
3313 : ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3314 : & (cell_atomb2(2)-1)*sc_size(3)*natom_uc+&
3315 : & (cell_atomb2(3)-1)*natom_uc+&
3316 475655760 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
3317 : else
3318 414979824 : ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
3319 : end if
3320 :
3321 890635584 : disp1 = displacement(idir2,ia2,itime)
3322 890635584 : disp2 = displacement(idir2,ib2,itime)
3323 :
3324 890635584 : tmp2 = tmp2 * (disp1-disp2)**power_disp
3325 890635584 : tmp3 = tmp3 * (disp1-disp2)**power_disp
3326 :
3327 : end if
3328 : end if
3329 : end do
3330 :
3331 527428432 : if(idisp1 > ndisp)then
3332 : ! Accumule stress tensor
3333 : strten_out(idir1,itime,icoeff_tmp) = strten_out(idir1,itime,icoeff_tmp) + &
3334 48745920 : & weight * tmp3
3335 : else
3336 : ! Accumule forces
3337 358437408 : fcart_out(idir1,ia1,icoeff_tmp,itime)=fcart_out(idir1,ia1,icoeff_tmp,itime)+weight*tmp2
3338 358437408 : fcart_out(idir1,ib1,icoeff_tmp,itime)=fcart_out(idir1,ib1,icoeff_tmp,itime)-weight*tmp2
3339 : end if
3340 : end do
3341 :
3342 : ! accumule energy
3343 124641984 : energy_out(icoeff_tmp,itime) = energy_out(icoeff_tmp,itime) + weight * tmp1
3344 :
3345 : end do!End do iterm
3346 : end do!End do coeff
3347 : end do!End time
3348 : end do!End do i3
3349 : end do!End do i2
3350 : end do!End do i1
3351 :
3352 : ! multiply by -1
3353 89048517 : fcart_out(:,:,:,:) = -1 * fcart_out(:,:,:,:)
3354 :
3355 : !ADD stress due to forces on atoms and variation of disp with strain
3356 23597 : do icoeff=1,ncoeff
3357 573207 : do itime=1,ntime
3358 22545809 : do ia1=1,natom_sc
3359 88487210 : do idir1=1,3
3360 483656800 : do idir2=1,6
3361 : strten_out(idir2,itime,icoeff) = strten_out(idir2,itime,icoeff) - &
3362 461672400 : & du_delta(idir2,idir1,ia1,itime)*fcart_out(idir1,ia1,icoeff,itime)
3363 : end do
3364 : end do
3365 : end do
3366 : end do
3367 : end do
3368 :
3369 561307 : do itime=1, ntime
3370 3858361 : do idir1=1,6
3371 3846563 : if (idir1 < 4) then
3372 3297357 : strten_out(idir1,itime,:) = strten_out(idir1,itime,:) * (1 + strain(idir1,itime))/ucvol(itime)
3373 : else
3374 3297357 : strten_out(idir1,itime,:) = strten_out(idir1,itime,:) * (1 - strain(idir1,itime)**2)/ucvol(itime)
3375 : end if
3376 : end do
3377 : end do
3378 :
3379 11798 : end subroutine fit_polynomial_coeff_getFS
3380 :
3381 :
3382 :
3383 :
3384 :
3385 : !> Computes mean squared deviations between model and reference data
3386 : !>
3387 : !> This routine:
3388 : !> - Evaluates the effective potential for all configurations
3389 : !> - Computes energy, force and stress differences vs reference
3390 : !> - Calculates weighted mean squared deviations
3391 : !> - Handles optional printing of detailed comparisons
3392 : !> - Supports electronic model evaluations if requested
3393 : !> - Can output phonon-specific analysis data
3394 : !>
3395 : !> @param eff_pot Effective potential to evaluate
3396 : !> @param hist Reference history with DFT data
3397 : !> @param mse Output mean squared error for energies
3398 : !> @param msef Output mean squared error for forces
3399 : !> @param mses Output mean squared error for stresses
3400 : !> @param natom Number of atoms in system
3401 : !> @param ntime Number of configurations
3402 : !> @param sqomega Sheppard-like weighting factors
3403 : !> @param comm MPI communicator
3404 : !> @param compute_anharmonic Flag to include anharmonic terms
3405 : !> @param print_file Flag to print detailed comparison files
3406 : !> @param filename Base name for output files
3407 : !> @param scup_dtset Optional SCUP dataset parameters
3408 : !> @param prt_ph Flag to print phonon analysis
3409 : !> @param weights Optional weights for each configuration
3410 56 : subroutine fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,sqomega,comm,&
3411 : & compute_anharmonic,print_file,filename,scup_dtset,&
3412 20 : & prt_ph, weights)
3413 :
3414 : implicit none
3415 :
3416 : !Arguments ------------------------------------
3417 : !scalars
3418 : integer, intent(in) :: natom,ntime,comm
3419 : real(dp),intent(out):: mse,msef,mses
3420 : logical,optional,intent(in) :: compute_anharmonic,print_file,prt_ph
3421 : !arrays
3422 : real(dp),intent(in) :: sqomega(ntime)
3423 : type(effective_potential_type),intent(in) :: eff_pot
3424 : type(abihist),intent(in) :: hist
3425 : !Strings/Characters
3426 : character(len=fnlen),optional,intent(in) :: filename
3427 : type(scup_dtset_type),optional,intent(inout) :: scup_dtset
3428 : real(dp),optional, intent(in) :: weights(ntime)
3429 : !Local variables-------------------------------
3430 : !scalar
3431 : integer :: ii,ia,mu,unit_energy,unit_stress,itime,master,nproc,my_rank,i
3432 : !Uncommend for dipdip test
3433 : integer :: ifirst
3434 : real(dp):: energy,energy_harm
3435 : logical :: need_anharmonic,need_print,need_elec_eval,iam_master
3436 : logical :: need_prt_ph
3437 : !arrays
3438 84 : real(dp):: fcart(3,natom),gred(3,natom),strten(6),rprimd(3,3),xred(3,natom)
3439 : !Strings/Characters
3440 : character(len=fnlen) :: file_energy, file_stress, file_anh, name_file
3441 : character(len=500) :: msg
3442 : !Uncommend for dipdip test
3443 28 : type(abihist) :: hist_out
3444 : character(len=200) :: filename_hist
3445 :
3446 56 : real(dp) :: weights_tmp(ntime)
3447 : ! *************************************************************************
3448 : !MS Hide SCALE-UP variables
3449 : ABI_UNUSED(itime)
3450 :
3451 : !MPI
3452 28 : master = 0
3453 28 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3454 28 : iam_master = (my_rank == master)
3455 :
3456 : !Do some checks
3457 28 : if(ntime /= hist%mxhist)then
3458 0 : write(msg,'(a)')'ntime is not correct'
3459 0 : ABI_BUG(msg)
3460 : end if
3461 :
3462 28 : if(natom /= size(hist%xred,2)) then
3463 0 : write(msg,'(a)')'natom is not correct'
3464 0 : ABI_BUG(msg)
3465 : end if
3466 :
3467 28 : need_anharmonic = .TRUE.
3468 28 : if(present(compute_anharmonic))then
3469 28 : need_anharmonic = compute_anharmonic
3470 : end if
3471 :
3472 28 : name_file=''
3473 28 : if(present(filename))name_file = filename
3474 :
3475 28 : need_print=.FALSE.
3476 28 : if(present(print_file))need_print=print_file
3477 :
3478 28 : need_elec_eval = .FALSE.
3479 28 : if(present(scup_dtset))need_elec_eval=scup_dtset%scup_elec_model
3480 :
3481 28 : need_prt_ph=.FALSE.
3482 28 : if(present(prt_ph))need_prt_ph=prt_ph
3483 :
3484 :
3485 28 : if(need_print .and. present(filename))then
3486 : !MS hist out uncommented for PHONOPY test
3487 21 : call abihist_init(hist_out,natom,ntime,.false.,.false.)
3488 21 : file_energy=trim(name_file)//'_energy.dat'
3489 21 : unit_energy = get_unit()
3490 21 : if (open_file(file_energy,msg,unit=unit_energy,form="formatted",&
3491 : & status="unknown",action="write") /= 0) then
3492 0 : ABI_ERROR(msg)
3493 : end if
3494 21 : unit_stress = get_unit()
3495 21 : file_stress=trim(name_file)//'_stress.dat'
3496 21 : if (open_file(file_stress,msg,unit=unit_stress,form="formatted",&
3497 : & status="unknown",action="write") /= 0) then
3498 0 : ABI_ERROR(msg)
3499 : end if
3500 7 : else if(need_print .and. .not. present(filename))then
3501 0 : write(msg,'(3a)')' You asked for printing of the MSD-values',ch10,&
3502 0 : & ' without specifying a filename'
3503 0 : ABI_ERROR(msg)
3504 : end if
3505 :
3506 28 : if(present(weights))then
3507 559 : weights_tmp = weights
3508 : else
3509 816 : weights_tmp = one
3510 : end if
3511 :
3512 28 : file_anh=trim(name_file)//'_anharmonic_terms_energy.dat'
3513 :
3514 28 : mse = zero
3515 28 : msef = zero
3516 28 : mses = zero
3517 1375 : do ii=1,ntime ! Loop over configurations
3518 216867 : xred(:,:) = hist%xred(:,:,ii)
3519 17511 : rprimd(:,:) = hist%rprimd(:,:,ii)
3520 1347 : itime = ii
3521 : #if defined DEV_MS_SCALEUP
3522 : !Pass print options to scale-up
3523 : itime = ii
3524 : if(need_elec_eval)then
3525 : call global_set_parent_iter(itime)
3526 : ! Set all print options to false.
3527 : call global_set_print_parameters(geom=.FALSE.,eigvals=.FALSE.,eltic=.FALSE.,&
3528 : & orbocc=.FALSE.,bands=.FALSE.)
3529 : if(ii == 1 .or. modulo(ii,scup_dtset%scup_printniter) == 0)then
3530 : call global_set_print_parameters(scup_dtset%scup_printgeom,scup_dtset%scup_printeigv,scup_dtset%scup_printeltic,&
3531 : & scup_dtset%scup_printorbocc,scup_dtset%scup_printbands)
3532 : end if
3533 : end if
3534 : #endif
3535 : call effective_potential_evaluate(eff_pot,energy_harm,fcart,gred,strten,natom,rprimd,&
3536 : & xred=xred,compute_anharmonic=.False.,verbose=.false.,&
3537 1347 : & elec_eval=need_elec_eval)
3538 :
3539 : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3540 : & xred=xred,compute_anharmonic=need_anharmonic,verbose=.false.,&
3541 1347 : & filename=file_anh,elec_eval=need_elec_eval)
3542 :
3543 1347 : if(need_print .and. iam_master)then
3544 552 : WRITE(unit_energy ,'(I10,7(F23.14))') ii,hist%etot(ii),energy_harm,energy,&
3545 1104 : & abs(hist%etot(ii) - energy_harm),abs(hist%etot(ii) - energy)
3546 552 : WRITE(unit_stress,'(I10,12(F23.14))') ii,hist%strten(:,ii),strten(:)
3547 : end if
3548 :
3549 : !MS Uncommented for abihist test
3550 1347 : if(need_prt_ph)then
3551 0 : if(ii == 1)then
3552 0 : write(msg,'(a,(80a))') ch10,('-',i=1,80)
3553 0 : call wrtout(ab_out,msg,'COLL')
3554 0 : call wrtout(std_out,msg,'COLL')
3555 0 : write(msg,'(3a)') ch10,'test_prt_ph == 1, write evulation of Model on the TEST-set into ph_test.nc',ch10
3556 0 : call wrtout(ab_out,msg,'COLL')
3557 0 : call wrtout(std_out,msg,'COLL')
3558 : endif
3559 0 : ifirst=merge(0,1,(ii>1))
3560 0 : filename_hist = trim("ph_test.nc")
3561 0 : hist_out%fcart(:,:,hist_out%ihist) = fcart(:,:)
3562 0 : hist_out%strten(:,hist_out%ihist) = strten(:)
3563 0 : hist_out%etot(hist_out%ihist) = energy
3564 0 : hist_out%entropy(hist_out%ihist) = hist%entropy(ii)
3565 0 : hist_out%time(hist_out%ihist) = real(ii,kind=dp)
3566 : ! call vel2hist(ab_mover%amass,hist,vel,vel_cell)
3567 0 : call var2hist(hist%acell(:,ii),hist_out,natom,hist%rprimd(:,:,ii),hist%xred(:,:,ii),.false.)
3568 0 : if(iam_master)then
3569 : call write_md_hist(hist_out,filename_hist,ifirst,ii,natom,1,eff_pot%crystal%ntypat,&
3570 : & eff_pot%supercell%typat,eff_pot%crystal%amu,eff_pot%crystal%znucl,&
3571 0 : & real(100,dp),(/real(100,dp),real(100,dp)/))
3572 : endif
3573 : endif!(need_prt_ph)
3574 :
3575 1347 : mse = mse + weights_tmp(ii) * ((hist%etot(ii) - energy))**2/(sqomega(ii)**(1.0/2.0)) !+abs(hist$etot(ii) - energy)
3576 55227 : do ia=1,natom ! Loop over atoms
3577 216867 : do mu=1,3 ! Loop over cartesian directions
3578 215520 : msef = msef + weights_tmp(ii) * (hist%fcart(mu,ia,ii) - fcart(mu,ia))**2
3579 : end do
3580 : end do
3581 9457 : do mu=1,6 ! Loop over stresses
3582 9429 : mses = mses + weights_tmp(ii) * sqomega(ii)*(hist%strten(mu,ii) - strten(mu))**2
3583 : end do
3584 : end do ! End loop itime
3585 28 : if(need_prt_ph)then
3586 0 : write(msg,'(a,(80a))') ch10,('-',i=1,80)
3587 0 : call wrtout(ab_out,msg,'COLL')
3588 0 : call wrtout(std_out,msg,'COLL')
3589 : endif
3590 :
3591 28 : mse = mse / ntime
3592 28 : msef = msef / (3*natom*ntime)
3593 28 : mses = mses / (6*ntime)
3594 :
3595 28 : if(need_print)then
3596 21 : close(unit_energy)
3597 21 : close(unit_stress)
3598 : end if
3599 :
3600 : !MS uncommented for PHONOPY TEST
3601 28 : call abihist_free(hist_out)
3602 :
3603 56 : end subroutine fit_polynomial_coeff_computeMSD
3604 :
3605 :
3606 :
3607 : !MARCUS_EXPERIMENTAL_SECTION
3608 : !!****f* m_fit_polynomiaL_coeff/testEffPot
3609 : !! NAME
3610 : !! testEffPot
3611 : !!
3612 : !! FUNCTION
3613 : !! Calculate the energy, forces for displacements provided
3614 : !! in an test-set (input:hist) within a given effective potential
3615 : !! (input: eff_pot)
3616 : !! If the test set is from DFT and contains DFT energies and forces
3617 : !! calculate the Goal Function values and the MSD of the Energy with
3618 : !! respect to the DFT energies
3619 : !!
3620 : !! INPUTS
3621 : !! eff_pot = effective_potential datatype
3622 : !! hist = abihist datatype
3623 : !!
3624 : !! OUTPUT
3625 : !!
3626 : !! SOURCE
3627 1 : subroutine fit_polynomial_coeff_testEffPot(eff_pot,hist,master,comm,print_anharmonic,scup_dtset,prt_ph)
3628 :
3629 :
3630 : implicit none
3631 :
3632 : !Arguments ------------------------------------
3633 : !scalars
3634 : integer,intent(in) :: master,comm
3635 : integer,optional,intent(in) :: prt_ph
3636 : !logicals
3637 : logical,optional,intent(in) :: print_anharmonic
3638 : !array
3639 : type(effective_potential_type),intent(inout) :: eff_pot
3640 : type(abihist),intent(in) :: hist
3641 : type(scup_dtset_type),optional,intent(inout) :: scup_dtset
3642 : !Local variables-------------------------------
3643 : !reals
3644 : real(dp) :: factor,mse,msef,mses
3645 1 : real(dp),allocatable :: sqomega(:),ucvol(:)
3646 : real(dp),parameter :: HaBohr_eVAng = Ha_eV /Bohr_Ang
3647 : !scalar
3648 : integer :: itime,unit_anh
3649 : integer :: natom,ntime,ncoeff,my_rank
3650 : !logicals
3651 : logical :: iam_master, need_print_anharmonic,file_opened,need_prt_ph
3652 : !strings/characters
3653 : character(len=fnlen) :: filename
3654 : character(len=1000) :: message
3655 : !arrays
3656 : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
3657 : ! *************************************************************************
3658 :
3659 : !MPI variables
3660 1 : my_rank=xmpi_comm_rank(comm)
3661 1 : iam_master = (my_rank == master)
3662 :
3663 : !Initialisation of optional arguments
3664 1 : need_print_anharmonic = .FALSE.
3665 1 : if(present(print_anharmonic)) need_print_anharmonic = print_anharmonic
3666 :
3667 1 : need_prt_ph = .FALSE.
3668 1 : if(present(prt_ph))then
3669 1 : if(prt_ph==1) need_prt_ph=.TRUE.
3670 : endif
3671 :
3672 : !Setting/Allocating other Variables
3673 1 : natom = size(hist%xred,2)
3674 1 : factor = 1._dp/natom
3675 1 : ntime = hist%mxhist
3676 3 : ABI_MALLOC(sqomega,(ntime))
3677 2 : ABI_MALLOC(ucvol,(ntime))
3678 102 : sqomega = zero
3679 1 : filename = 'TES_fit_diff'
3680 1 : ncoeff = eff_pot%anharmonics_terms%ncoeff
3681 :
3682 102 : do itime=1,ntime
3683 : ! Compute \Omega^{2} and ucvol for each time
3684 101 : call metric(gmet,gprimd,-1,rmet,hist%rprimd(:,:,itime),ucvol(itime))
3685 : ! Formula: sqomega(itime) = (((ucvol(itime)**(-2.))* ((natom)**(0.5)))**(-1.0/3.0))**2
3686 : ! Compact form:
3687 102 : sqomega(itime) = ((ucvol(itime)**(4.0/3.0)) / ((natom)**(1/3.0)))
3688 : end do
3689 :
3690 :
3691 1 : if(need_print_anharmonic) call effective_potential_writeAnhHead(ncoeff,&
3692 1 : & filename,eff_pot%anharmonics_terms)
3693 :
3694 : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,&
3695 : & sqomega,comm,&
3696 : & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=filename,&
3697 1 : & scup_dtset=scup_dtset,prt_ph=need_prt_ph)
3698 :
3699 :
3700 : ! Print the standard deviation after the fit
3701 1 : write(message,'(6a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
3702 : !& ' Mean Standard Deviation values of the effective-potential',ch10,&
3703 : !& ' with respect to the test-set (meV^2/atm):',&
3704 : !& ch10,' Energy : ',&
3705 : !& mse* (Ha_EV*1000)**2 *factor ,ch10,&
3706 1 : & ' Goal function values of the effective.potential',ch10,&
3707 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
3708 1 : & ' Energy : ',&
3709 1 : & (mse)*(HaBohr_eVAng)**2,ch10,&
3710 1 : & ' Forces+Stresses : ',&
3711 1 : & (msef+mses)*(HaBohr_eVAng)**2,ch10,&
3712 1 : & ' Forces : ',&
3713 1 : & msef*(HaBohr_eVAng)**2,ch10,&
3714 1 : & ' Stresses : ',&
3715 2 : & mses*(HaBohr_eVAng)**2,ch10
3716 1 : call wrtout(ab_out,message,'COLL')
3717 1 : call wrtout(std_out,message,'COLL')
3718 :
3719 :
3720 : !Deallocating
3721 1 : ABI_FREE(sqomega)
3722 1 : ABI_FREE(ucvol)
3723 :
3724 1 : INQUIRE(FILE='TES_fit_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh)
3725 1 : if(file_opened) close(unit_anh)
3726 :
3727 :
3728 1 : end subroutine fit_polynomial_coeff_testEffPot
3729 : !!***
3730 :
3731 :
3732 :
3733 : !> Prints system files needed for external fitting scripts
3734 : !>
3735 : !> This routine writes:
3736 : !> - Born effective charges
3737 : !> - Dielectric tensor
3738 : !> - Reference structure and energy
3739 : !> - Harmonic force constants in XML format
3740 : !> - Strain tensors
3741 : !> - Symmetry operations
3742 : !> - Molecular dynamics trajectory data
3743 : !>
3744 : !> The files are organized in a standard format expected by external fitting tools.
3745 : !> Special care is taken to properly order atoms and handle supercell transformations.
3746 : !>
3747 : !> @param eff_pot Effective potential containing system parameters
3748 : !> @param hist History containing MD trajectory data
3749 0 : subroutine fit_polynomial_printSystemFiles(eff_pot,hist)
3750 :
3751 : implicit none
3752 :
3753 : !Arguments ------------------------------------
3754 : !scalars
3755 : !arrays
3756 : type(effective_potential_type), intent(in) :: eff_pot
3757 : type(abihist),intent(in) :: hist
3758 : !Local variables-------------------------------
3759 : !scalar
3760 : integer :: ia,ib,ib1,ii,jj,irpt,kk,ll,mu,nu,nstep,nshift
3761 : integer :: natom_uc
3762 : integer :: unit_born=22,unit_epsiloninf=23,unit_md=24
3763 : integer :: unit_harmonic=25,unit_ref=26,unit_strain=27,unit_sym=28
3764 : !arrays
3765 0 : integer,allocatable :: typat_order(:),typat_order_uc(:)
3766 : integer, dimension(3) :: A,ncell
3767 0 : real(dp), allocatable :: xcart(:,:),fcart(:,:)
3768 : character(len=500) :: msg
3769 0 : type(supercell_type) :: supercell
3770 : ! *************************************************************************
3771 :
3772 : !Create new supercell corresponding to the MD
3773 0 : ncell = (/2,2,2/)
3774 : call supercell%init(eff_pot%crystal%natom, (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/),&
3775 : eff_pot%crystal%rprimd,eff_pot%crystal%typat,&
3776 0 : eff_pot%crystal%xcart,eff_pot%crystal%znucl)
3777 :
3778 : !allocation of array
3779 0 : ABI_MALLOC(xcart,(3,supercell%natom))
3780 0 : ABI_MALLOC(fcart,(3,supercell%natom))
3781 0 : ABI_MALLOC(typat_order,(supercell%natom))
3782 0 : ABI_MALLOC(typat_order_uc,(eff_pot%crystal%natom))
3783 :
3784 0 : A = (/ 2, 3, 1/)
3785 :
3786 0 : nshift = product(ncell)
3787 0 : natom_uc = eff_pot%crystal%natom
3788 : !Fill the typat_order array:
3789 : !In the fit script the atom must be in the order 11111 222222 33333 ..
3790 : !and the order of the atom can not be change in the fit script,
3791 : !we transform into the format of the script
3792 0 : ib = 1
3793 0 : ib1= 1
3794 0 : do ii=1,eff_pot%crystal%ntypat
3795 0 : jj = A(ii)
3796 0 : do kk=1,natom_uc
3797 0 : if(supercell%typat(kk)==jj)then
3798 0 : typat_order_uc(ib1) = kk
3799 0 : ib1 = ib1 + 1
3800 0 : do ll=1,nshift
3801 0 : ia = (ll-1)*natom_uc + kk
3802 0 : typat_order(ib) = ia
3803 0 : ib = ib + 1
3804 : end do
3805 : end if
3806 : end do
3807 : end do
3808 :
3809 : ! BORN CHARGES FILE
3810 0 : if (open_file('system/Born_Charges',msg,unit=unit_born,form="formatted",&
3811 : & status="replace",action="write") /= 0) then
3812 0 : ABI_ERROR(msg)
3813 : end if
3814 0 : do ii=1,eff_pot%crystal%ntypat
3815 0 : jj = A(ii)
3816 0 : do ia=1,eff_pot%crystal%natom
3817 0 : if(eff_pot%crystal%typat(ia)==jj)then
3818 0 : write(unit_born,'(i2,a,1F10.5)') ia," ",eff_pot%crystal%amu(eff_pot%crystal%typat(ia))
3819 0 : do mu=1,3
3820 0 : WRITE(unit_born,'(a,3(F23.14))') " ",eff_pot%harmonics_terms%zeff(:,mu,ia)
3821 : end do
3822 : end if
3823 : end do
3824 : end do
3825 :
3826 : !DIELECTRIC TENSOR FILE
3827 0 : if (open_file('system/Dielectric_Tensor',msg,unit=unit_epsiloninf,form="formatted",&
3828 : & status="replace",action="write") /= 0) then
3829 0 : ABI_ERROR(msg)
3830 : end if
3831 0 : do mu=1,3
3832 0 : WRITE(unit_epsiloninf,'(3(F23.14))') eff_pot%harmonics_terms%epsilon_inf(:,mu)
3833 : end do
3834 :
3835 :
3836 : !REFERENCE STRUCTURE FILE
3837 0 : if (open_file('system/Reference_structure',msg,unit=unit_ref,form="formatted",&
3838 : & status="replace",action="write") /= 0) then
3839 0 : ABI_ERROR(msg)
3840 : end if
3841 :
3842 0 : write(unit_ref,'("Energy (Hartree)")')
3843 0 : write(unit_ref,'("================")')
3844 0 : write(unit_ref,'(F23.14)') (hist%etot(1)/nshift)
3845 0 : write(unit_ref,'("")')
3846 0 : write(unit_ref,'("Cell vectors")')
3847 0 : write(unit_ref,'("============")')
3848 0 : do jj=1,3
3849 0 : write(unit_ref,'(3(F22.14))') (supercell%rprimd(:,jj))
3850 : end do
3851 :
3852 0 : write(unit_ref,'("")')
3853 0 : write(unit_ref,'("Atomic positions (Bohr radius)")')
3854 0 : write(unit_ref,'("==============================")')
3855 :
3856 0 : do ia=1,supercell%natom
3857 0 : write(unit_ref,'(3(F23.14))') supercell%xcart(:,typat_order(ia))
3858 : end do
3859 :
3860 : !Harmonic XML file
3861 0 : if (open_file('system/harmonic.xml',msg,unit=unit_harmonic,form="formatted",&
3862 : & status="replace",action="write") /= 0) then
3863 0 : ABI_ERROR(msg)
3864 : end if
3865 :
3866 : !Write header
3867 0 : write(unit_harmonic,'("<?xml version=""1.0"" ?>")')
3868 0 : write(unit_harmonic,'("<name>")')
3869 :
3870 0 : do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt
3871 0 : if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol9)) then
3872 0 : write(unit_harmonic,'(" <local_force_constant units=""hartree/bohrradius**2"">")')
3873 0 : write(unit_harmonic,'(" <data>")')
3874 0 : do ia=1,eff_pot%crystal%natom
3875 0 : do mu=1,3
3876 0 : do ib=1,eff_pot%crystal%natom
3877 0 : do nu=1,3
3878 : write(unit_harmonic,'(F22.14)', advance="no")&
3879 : & (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,typat_order_uc(ia),&
3880 0 : & nu,typat_order_uc(ib),irpt))
3881 : end do
3882 : end do
3883 0 : write(unit_harmonic,'(a)')''
3884 : end do
3885 : end do
3886 0 : write(unit_harmonic,'(" </data>")')
3887 0 : write(unit_harmonic,'(" <cell>")')
3888 0 : write(unit_harmonic,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
3889 0 : write(unit_harmonic,'(" </cell>")')
3890 0 : write(unit_harmonic,'(" </local_force_constant>")')
3891 : end if
3892 : end do
3893 0 : write(unit_harmonic,'("</name>")')
3894 :
3895 : !STRAIN FILE
3896 0 : if (open_file('system/Strain_Tensor',msg,unit=unit_strain,form="formatted",&
3897 : & status="replace",action="write") /= 0) then
3898 0 : ABI_ERROR(msg)
3899 : end if
3900 0 : write(unit_strain,'(6(F23.14))') (eff_pot%harmonics_terms%elastic_constants)
3901 :
3902 : ! SYM FILE
3903 0 : if (open_file('system/symmetry_operations',msg,unit=unit_sym,form="formatted",&
3904 : & status="replace",action="write") /= 0) then
3905 0 : ABI_ERROR(msg)
3906 : end if
3907 : write(unit_sym,'("(x,y,z) (y,-x,z) (z,x,y) (y,z,x) (x,z,y) (y,x,z) (z,y,x) (x,-y,-z) (z,-x,-y)",&
3908 : & " (y,-z,-x) (x,-z,-y) (y,-x,-z) (z,-y,-x) (-x,y,-z) (-z,x,-y) (-y,z,-x) (-x,z,-y)",&
3909 : & " (-y,x,-z) (-z,y,-x) (-x,-y,z) (-z,-x,y) (-y,-z,x) (-x,-z,y) (-y,-x,z) (-z,-y,x)",&
3910 : & " (-x,-y,-z) (-z,-x,-y) (-y,-z,-x) (-x,-z,-y) (-y,-x,-z) (-z,-y,-x) (-x,y,z)",&
3911 : & " (-z,x,y) (-y,z,x) (-x,z,y) (-y,x,z) (-z,y,x) (x,-y,z) (z,-x,y) (y,-z,x) (x,-z,y)",&
3912 0 : & " (z,-y,x) (x,y,-z) (z,x,-y) (y,z,-x) (x,z,-y) (y,x,-z) (z,y,-x)")')
3913 :
3914 :
3915 : !MD file
3916 0 : nstep = hist%mxhist
3917 0 : if (open_file('system/Molecular_dynamic',msg,unit=unit_md,form="formatted",&
3918 : & status="replace",action="write") /= 0) then
3919 0 : ABI_ERROR(msg)
3920 : end if
3921 0 : do ii=1,nstep
3922 0 : write(unit_md,'(I5)') ii-1
3923 0 : write(unit_md,'(F22.14)') hist%etot(ii)/nshift
3924 0 : do jj=1,3
3925 0 : write(unit_md,'(3(F22.14))') (hist%rprimd(:,jj,ii))
3926 : end do
3927 : ! Set xcart and fcart for this step
3928 : call xred2xcart(supercell%natom,hist%rprimd(:,:,ii),&
3929 0 : & xcart,hist%xred(:,:,ii))
3930 :
3931 0 : fcart(:,:) = hist%fcart(:,:,ii)
3932 :
3933 0 : do ia=1,supercell%natom
3934 0 : write(unit_md,'(3(E22.14),3(E22.14))') xcart(:,typat_order(ia)),fcart(:,typat_order(ia))
3935 : end do
3936 0 : write(unit_md,'(6(E22.14))') hist%strten(:,ii)
3937 : end do
3938 :
3939 : !Close files
3940 0 : close(unit_ref)
3941 0 : close(unit_born)
3942 0 : close(unit_harmonic)
3943 0 : close(unit_epsiloninf)
3944 0 : close(unit_md)
3945 0 : close(unit_strain)
3946 0 : close(unit_sym)
3947 :
3948 : !Deallocation array
3949 0 : ABI_FREE(typat_order)
3950 0 : ABI_FREE(typat_order_uc)
3951 0 : ABI_FREE(xcart)
3952 0 : ABI_FREE(fcart)
3953 0 : call supercell%free()
3954 :
3955 0 : end subroutine fit_polynomial_printSystemFiles
3956 :
3957 :
3958 : !> Recursively generates combinations of numbers with specific conditions
3959 : !>
3960 : !> This routine:
3961 : !> - Generates combinations of numbers from 1 to n_max of length m_max
3962 : !> - Ensures each number in a combination is greater than previous one
3963 : !> - Can either count combinations or store them based on compute flag
3964 : !> - Uses recursive approach to build combinations
3965 : !>
3966 : !> @param i Current index/counter for output combinations
3967 : !> @param m Current position in combination being built
3968 : !> @param m_max Maximum length of combinations
3969 : !> @param n_max Maximum number to use in combinations
3970 : !> @param list Working array for current combination
3971 : !> @param list_out Output array for storing combinations
3972 : !> @param size Size of output array
3973 : !> @param compute Whether to store combinations or just count them
3974 0 : recursive subroutine genereList(i,m,m_max,n_max,list,list_out,size,compute)
3975 :
3976 : implicit none
3977 :
3978 : !Arguments ---------------------------------------------
3979 : !scalar
3980 : integer, intent(in) :: m_max,n_max,m,size
3981 : integer, intent(inout) :: i
3982 : logical,intent(in) :: compute
3983 : !arrays
3984 : integer, intent(out) :: list(m_max),list_out(size,m_max)
3985 : !Local variables ---------------------------------------
3986 : !scalar
3987 : integer n
3988 : !arrays
3989 :
3990 : ! *************************************************************************
3991 0 : if (m > m_max) then
3992 0 : i = i + 1
3993 0 : if(compute)list_out(i,:) = list(:)
3994 : else
3995 0 : do n = 1, n_max
3996 0 : if (m == 1)then
3997 0 : list(m) = n
3998 0 : call genereList (i, m + 1,m_max,n_max,list,list_out,size,compute)
3999 0 : else if (n > list(m - 1)) then
4000 0 : list(m) = n
4001 0 : call genereList (i, m + 1,m_max,n_max,list,list_out,size,compute)
4002 : end if
4003 : end do
4004 : end if
4005 :
4006 0 : end subroutine genereList
4007 : !!***
4008 :
4009 : !> Gathers arrays from all MPI processes with varying lengths
4010 : !>
4011 : !> This routine:
4012 : !> - Gathers real values and their corresponding ordering indices
4013 : !> - Handles arrays of different sizes from different processes
4014 : !> - Allocates output buffers based on total gathered size
4015 : !> - Sets up displacement arrays for MPI communication
4016 : !> - Performs AllGatherv for both values and indices
4017 : !>
4018 : !> @param A Input array of values from this process
4019 : !> @param order Input array of indices from this process
4020 : !> @param n Size of input arrays from this process
4021 : !> @param buff_A Output buffer for gathered values
4022 : !> @param buff_order Output buffer for gathered indices
4023 : !> @param totsize Total size across all processes
4024 : !> @param comm MPI communicator
4025 : !> @param nproc Number of MPI processes
4026 86 : subroutine mpigatherv(A, order, n ,buff_A, buff_order, totsize, comm, nproc)
4027 : integer, intent(inout) :: n, totsize
4028 : real(dp), intent(inout) :: A(n)
4029 : integer, intent(inout):: order(n)
4030 : integer, intent(in) :: comm, nproc
4031 : real(dp), allocatable, intent(out) :: buff_A(:)
4032 : integer, allocatable, intent(out) :: buff_order(:)
4033 172 : integer :: disps(nproc), sizes(nproc)
4034 : integer :: ierr, i
4035 86 : totsize=n
4036 86 : call xmpi_sum(totsize, comm, ierr)
4037 258 : ABI_MALLOC(buff_A, (totsize))
4038 258 : ABI_MALLOC(buff_order, (totsize))
4039 86 : call xmpi_allgather(n, sizes, comm, ierr)
4040 86 : disps(1)=0
4041 198 : do i=2, nproc
4042 198 : disps(i)=disps(i-1)+sizes(i-1)
4043 : end do
4044 86 : call xmpi_allgatherv(A, n, buff_A, sizes, disps, comm, ierr )
4045 86 : call xmpi_allgatherv(order, n, buff_order, sizes, disps, comm, ierr )
4046 86 : end subroutine mpigatherv
4047 :
4048 :
4049 :
4050 : end module m_fit_polynomial_coeff
4051 : !!***
|
