Branch data 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 : 11740 : real(dp),allocatable:: fcart_coeffs_tmp(:,:,:)
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 [ - + + - : 46960 : ABI_MALLOC(fcart_coeffs_tmp,(3,natom,ntime))
- + + - +
- + - - +
- + ]
2751 [ + + + + : 465112 : A=zero; B=zero;
+ + + + ]
2752 [ + + + + : 453372 : AF = zero; IPIV = 1;
+ + ]
2753 [ + + + + ]: 120628 : R = one; C = one;
2754 [ + + + + ]: 35220 : FERR = zero; BERR = zero
2755 [ + + + + ]: 283960 : IWORK = 0; WORK = 0
2756 : :
2757 : : !1-Get forces and stresses from the model and fill A
2758 : : ! Fill alsor B with the forces and stresses from
2759 : : ! the DFT snapshot and the model
2760 : : ! See equation 17 of PRB95 094115 (2017) [[cite:Escorihuela-Sayalero2017]]
2761 [ + + ]: 66184 : do icoeff=1,ncoeff_fit
2762 : 54444 : icoeff_tmp = list_coeffs(icoeff)
2763 [ + + + + : 530387800 : fcart_coeffs_tmp(:,:,:) = fcart_coeffs(:,:,icoeff_tmp,:)
+ + ]
2764 : 3348440 : ftmpA= zero; ftmpB = zero
2765 : 3348440 : stmpA= zero; stmpB = zero
2766 : 3348440 : etmpA= zero; etmpB = zero
2767 : : ! loop over the configuration
2768 [ + + ]: 3348440 : do itime=1,ntime
2769 : : ! Fill energy
2770 : 3293996 : emu = energy_coeffs(icoeff_tmp,itime)
2771 [ + + ]: 29312412 : do jcoeff=1,ncoeff_fit
2772 : 26018416 : jcoeff_tmp = list_coeffs(jcoeff)
2773 : 26018416 : enu = energy_coeffs(jcoeff_tmp,itime)
2774 [ + + ]: 29312412 : if(fit_on(3))then
2775 : 361889 : etmpA = emu*enu/(sqomega(itime)**(1.0/2.0))
2776 : 361889 : A(icoeff,jcoeff) = A(icoeff,jcoeff) + efact*etmpA* weights(itime)
2777 : : endif
2778 : : end do
2779 [ + + ]: 3293996 : if(fit_on(3))then
2780 : 73403 : etmpB = etmpB + energy_diff(itime)*emu/(sqomega(itime)**(1.0/2.0)) * weights(itime)!/ (sqomega(itime)**3)
2781 : : else
2782 : : etmpB = zero ! REMOVE THIS LINE TO TAKE INTO ACOUNT THE ENERGY
2783 : : endif
2784 : : ! Fill forces
2785 [ + + ]: 135053836 : do ia=1,natom
2786 [ + + ]: 530333356 : do mu=1,3
2787 : 395279520 : fmu = fcart_coeffs_tmp(mu,ia,itime)
2788 [ + + ]: 3517489440 : do jcoeff=1,ncoeff_fit
2789 : 3122209920 : jcoeff_tmp = list_coeffs(jcoeff)
2790 : 3122209920 : fnu = fcart_coeffs(mu,ia,jcoeff_tmp,itime)
2791 : 3122209920 : ftmpA = fmu*fnu
2792 [ + - ]: 3517489440 : if(fit_on(1))A(icoeff,jcoeff) = A(icoeff,jcoeff) + ffact*ftmpA* weights(itime)
2793 : : end do
2794 [ + - ]: 527039360 : if(fit_on(1))then
2795 : 395279520 : ftmpB = ftmpB + fcart_diff(mu,ia,itime)*fmu* weights(itime)
2796 : : else
2797 : : ftmpB = zero
2798 : : endif
2799 : : end do !End loop dir
2800 : : end do !End loop natom
2801 : : ! Fill stresses
2802 [ + + ]: 23112416 : do mu=1,6
2803 : 19763976 : smu = strten_coeffs(mu,itime,icoeff_tmp)
2804 [ + + ]: 175874472 : do jcoeff=1,ncoeff_fit
2805 : 156110496 : jcoeff_tmp = list_coeffs(jcoeff)
2806 : 156110496 : snu = strten_coeffs(mu,itime,jcoeff_tmp)
2807 : 156110496 : stmpA = sqomega(itime)*smu*snu
2808 [ + - ]: 175874472 : if(fit_on(2))A(icoeff,jcoeff) = A(icoeff,jcoeff) + sfact*stmpA *weights(itime)
2809 : : end do
2810 [ + - ]: 23057972 : if(fit_on(2))then
2811 : 19763976 : stmpB = stmpB + sqomega(itime)*strten_diff(mu,itime)*smu * weights(itime)
2812 : : else
2813 : : stmpB = zero
2814 : : endif
2815 : : end do !End loop stress dir
2816 : : end do ! End loop time
2817 : 66184 : B(icoeff,1) = B(icoeff,1) + ffact*ftmpB + sfact*stmpB + efact*etmpB
2818 : : end do ! End loop icoeff
2819 : :
2820 : : !2-Solve Ax=B
2821 : :
2822 : : !== Three different LAPACK routines are possible ==
2823 : : !OLD VERSION..
2824 : : ! call dgesvx(FACT,TRANS,N,NRHS,A,LDA,AF,LDAF,IPIV,EQUED,R,C,B,LDB,coefficients,LDX,&
2825 : : !RCOND,FERR,BERR,WORK,IWORK,INFO)
2826 : : !U is nonsingular
2827 : : ! if (INFO==N+1) then
2828 : : ! coefficients = zero
2829 : : ! end if
2830 : :
2831 : : !VERSION PRIOR TO 20240817. However test multibinit multi_l_6_1 fails on the new reference machine (eos).
2832 : : !call DSGESV(N,NRHS,A,LDA,IPIV,B,LDB,coefficients,LDX,WORK,SWORK,ITER,INFO)
2833 : :
2834 : : !Other routine, activated on 20240817, together with the change of reference machine
2835 : : ! Check matrix conditioning
2836 : : !if (any(isnan(A))) then
2837 : : ! INFO = -2
2838 : : ! print *, "NaN in A"
2839 : : ! coefficients = zero
2840 : : !endif
2841 : :
2842 : : !if (any(isnan(B))) then
2843 : : ! INFO = -2
2844 : : ! !print *, "NaN in B"
2845 : : ! coefficients = zero
2846 : : !endif
2847 : : ! Solve the system
2848 : : ! call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
2849 : : !if (INFO == 0) then
2850 : : ! Check solution for NaN or very large values
2851 : : ! if (any(isnan(B(:,NRHS))) .or. any(abs(B(:,NRHS)) > 1.0e10_dp)) then
2852 : : ! INFO = -3
2853 : : ! coefficients = zero
2854 : : ! !print *, "NaN or very large values in coefficients"
2855 : : ! else
2856 : : ! coefficients = B(:,NRHS)
2857 : : ! print *, "INFO=", INFO
2858 : : ! print *, "Good solution", coefficients
2859 : : ! endif
2860 : : !else
2861 : : ! print*, "INFO=", INFO, "setting coefficients to zero"
2862 : : ! coefficients = zero
2863 : : !end if
2864 : :
2865 : : !==================================================
2866 : :
2867 : :
2868 : 11740 : call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO)
2869 [ + + ]: 66184 : coefficients = B(:,NRHS)
2870 : :
2871 : : !U is nonsingular
2872 [ - + ]: 11740 : if (INFO==N+2) then
2873 [ # # ]: 0 : coefficients = zero
2874 : : end if
2875 : :
2876 [ + + + + : 132299 : if(any(abs(coefficients)>1.0E10) .or. any(coefficients /= coefficients))then
+ + + - +
+ ]
2877 : 69 : INFO = 1
2878 [ + + ]: 218 : coefficients = zero
2879 : : end if
2880 : :
2881 : 11740 : info_out = INFO
2882 : :
2883 [ - + ]: 11740 : ABI_FREE(AF)
2884 [ - + ]: 11740 : ABI_FREE(IPIV)
2885 [ - + ]: 11740 : ABI_FREE(R)
2886 [ - + ]: 11740 : ABI_FREE(C)
2887 [ - + ]: 11740 : ABI_FREE(FERR)
2888 [ - + ]: 11740 : ABI_FREE(BERR)
2889 [ - + ]: 11740 : ABI_FREE(WORK)
2890 [ - + ]: 11740 : ABI_FREE(IWORK)
2891 [ - + ]: 11740 : ABI_FREE(SWORK)
2892 [ - + ]: 11740 : ABI_FREE(A)
2893 [ - + ]: 11740 : ABI_FREE(B)
2894 [ - + ]: 11740 : ABI_FREE(fcart_coeffs_tmp)
2895 : :
2896 : 11740 : end subroutine fit_polynomial_coeff_solve
2897 : : !!***
2898 : :
2899 : :
2900 : :
2901 : : !> Computes weights for history configurations based on forces and temperature
2902 : : !>
2903 : : !> This routine:
2904 : : !> - Calculates average force magnitudes for each configuration
2905 : : !> - Applies Boltzmann weighting based on input temperature
2906 : : !> - Uses atomic forces to determine configuration importance
2907 : : !> - Normalizes weights to sum to ntime
2908 : : !> - Outputs detailed statistics about forces and weights
2909 : : !> - Handles both temperature-dependent and uniform weighting
2910 : : !>
2911 : : !> @param hist Input history containing force data
2912 : : !> @param temperature Temperature in K for Boltzmann weighting
2913 : : !> @param ntime Number of timesteps in history
2914 : : !> @param natom Number of atoms in system
2915 : : !> @param weights Output array of weights for each configuration
2916 : : !> @param comm MPI communicator
2917 : 40 : subroutine get_weight_from_hist(hist, temperature, ntime, natom, weights, comm)
2918 : : type(abihist), intent(in) :: hist
2919 : : !real(dp), intent(in) :: DFT_forces(3,natom,ntime)
2920 : : real(dp), intent(in) :: temperature
2921 : : integer, intent(in) :: ntime, natom
2922 : : integer, intent(in) :: comm
2923 : : real(dp), allocatable, intent(out) :: weights(:)
2924 : : integer :: itime, iatom
2925 : : ! Boltzmann constant in Ha/K
2926 : : real(dp), parameter :: kb = 3.166815d-6
2927 : 80 : real(dp) :: average_forces(ntime)
2928 : 40 : character(len=ntime*100) :: msg
2929 : : integer :: master, my_rank
2930 : : logical :: iam_master
2931 : :
2932 : 40 : master = 0
2933 : 40 : my_rank = xmpi_comm_rank(comm)
2934 : 40 : iam_master = (my_rank == master)
2935 [ - + + - : 120 : ABI_MALLOC(weights, (ntime))
+ - + - -
+ - + -
+ ]
2936 [ + + ]: 799 : weights=1.0_dp
2937 : :
2938 : : ! compute average forces for each time step
2939 [ + + ]: 799 : average_forces = 0.0_dp
2940 [ + + ]: 799 : do itime=1,ntime
2941 [ + + ]: 31119 : do iatom=1,natom
2942 [ + + + - : 122199 : average_forces(itime) = average_forces(itime) + norm2(hist%fcart(:,iatom,itime))
- + ]
2943 : : end do
2944 : 799 : average_forces(itime) = average_forces(itime)/natom
2945 : : end do
2946 : :
2947 [ + + ]: 40 : if (temperature > 0.0_dp) then
2948 : : ! compute the average norm of force for each atom
2949 : : ! compute the weight
2950 [ + + ]: 36 : weights = 0.0_dp
2951 [ + + ]: 36 : do itime=1,ntime
2952 : 33 : weights(itime) = exp(-average_forces(itime)/(kb*temperature))
2953 [ + + ]: 36 : if (weights(itime)< 0.001) then
2954 : 21 : weights(itime) = 0.001
2955 : : endif
2956 : : end do
2957 : : ! normalize the weight
2958 [ + + + + ]: 69 : weights = weights/sum(weights)* ntime
2959 : : end if
2960 : : ! For each time, print the average forces, the energy, the weight,
2961 : : ! the maximum of axial stress, and maximum of shear stress
2962 : : ! first write the header of the table
2963 : : ! write a separator
2964 : :
2965 [ + + ]: 40 : if(iam_master)then
2966 : 24 : write(msg, '(a)') '-----------------------------------------------------------------------------------'
2967 : 24 : write(msg, '(a)') 'Time Average_forces Etot-min(Etot) Weight Max(axial_stress) Max(shear_stress)'
2968 : 24 : call wrtout(std_out,msg,'COLL')
2969 : : ! then write the data
2970 [ + + ]: 607 : do itime=1,ntime
2971 : 583 : write(msg, '(i5, 3e15.6, 2e15.6)') itime, average_forces(itime), &
2972 [ + - - + : 46609 : & hist%etot(itime)-minval(hist%etot(:)), weights(itime), &
- - + + ]
2973 [ - + + - : 6413 : & maxval(abs(hist%strten(1:3,itime))), maxval(abs(hist%strten(4:6,itime)))
+ + - + +
- + + ]
2974 : 607 : call wrtout(std_out,msg,'COLL')
2975 : : end do
2976 : 24 : write(msg, '(a)') '-----------------------------------------------------------------------------------'
2977 : : endif
2978 : :
2979 : :
2980 : 40 : end subroutine get_weight_from_hist
2981 : :
2982 : :
2983 : :
2984 : :
2985 : : !> Computes goal function values comparing model predictions to reference data
2986 : : !>
2987 : : !> This routine calculates:
2988 : : !> - Forces + stresses combined goal function (gf_value(1))
2989 : : !> - Forces only goal function (gf_value(2))
2990 : : !> - Stresses only goal function (gf_value(3))
2991 : : !> - Energy goal function (gf_value(4))
2992 : : !>
2993 : : !> The calculations:
2994 : : !> - Use supplied coefficient values to evaluate model
2995 : : !> - Compare to reference energy/force/stress differences
2996 : : !> - Apply configuration weights if provided
2997 : : !> - Handle volume/supercell normalization
2998 : : !> - Support parallel evaluation across processes
2999 : : !>
3000 : : !> @param coefficients Model coefficient values
3001 : : !> @param energy_coeffs Energy terms for each coefficient
3002 : : !> @param energy_diff Reference energy differences
3003 : : !> @param fcart_coeffs Force terms for each coefficient
3004 : : !> @param fcart_diff Reference force differences
3005 : : !> @param gf_value Output goal function values
3006 : : !> @param list_coeffs List of coefficients to include
3007 : : !> @param natom Number of atoms
3008 : : !> @param ncoeff_fit Number of coefficients being fit
3009 : : !> @param ncoeff_max Maximum number of coefficients
3010 : : !> @param ntime Number of configurations
3011 : : !> @param strten_coeffs Stress terms for each coefficient
3012 : : !> @param strten_diff Reference stress differences
3013 : : !> @param sqomega Configuration weighting factors
3014 : : !> @param weights Optional weights for each configuration
3015 : 11714 : subroutine fit_polynomial_coeff_computeGF(coefficients,energy_coeffs,energy_diff,&
3016 : 11714 : & fcart_coeffs,fcart_diff,gf_value,list_coeffs,&
3017 : 11714 : & natom,ncoeff_fit,ncoeff_max,ntime,strten_coeffs,&
3018 : 11714 : & strten_diff,sqomega, weights)
3019 : :
3020 : : implicit none
3021 : :
3022 : : !Arguments ------------------------------------
3023 : : !scalars
3024 : : integer,intent(in) :: natom,ncoeff_fit,ncoeff_max,ntime
3025 : : !arrays
3026 : : integer,intent(in) :: list_coeffs(ncoeff_fit)
3027 : : real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime)
3028 : : real(dp),intent(in) :: energy_diff(ntime)
3029 : : real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime)
3030 : : real(dp),intent(in) :: fcart_diff(3,natom,ntime)
3031 : : real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max)
3032 : : real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime)
3033 : : real(dp),intent(in) :: coefficients(ncoeff_fit)
3034 : : real(dp), intent(in) :: weights(ntime)
3035 : : real(dp),intent(out) :: gf_value(4)
3036 : : !Local variables-------------------------------
3037 : : !scalar
3038 : : integer :: ia,icoeff,icoeff_tmp,itime,mu
3039 : : real(dp):: etmp,emu,fmu,ftmp,smu,stmp
3040 : : real(dp) :: ffact,sfact,efact
3041 : : !arrays
3042 : : ! *************************************************************************
3043 : :
3044 : : !1-Compute the value of the goal function
3045 : : ! see equation 9 of PRB 95 094115(2017) [[cite:Escorihuela-Sayalero2017]]
3046 : : gf_value = zero
3047 : 11714 : etmp = zero
3048 : 11714 : ftmp = zero
3049 : 11714 : stmp = zero
3050 : :
3051 : : !Compute factors
3052 : 11714 : ffact = one/(3*natom*ntime)
3053 : 11714 : sfact = one/(6*ntime)
3054 : 11714 : efact = one/(ntime)
3055 : :
3056 : : ! loop over the configuration
3057 [ + + ]: 558819 : do itime=1,ntime
3058 : : ! Fill energy
3059 : : emu = zero
3060 [ + + ]: 3839392 : do icoeff=1,ncoeff_fit
3061 : 3292287 : icoeff_tmp = list_coeffs(icoeff)
3062 : 3839392 : emu = emu + coefficients(icoeff)*energy_coeffs(icoeff_tmp,itime)
3063 : : end do
3064 : : ! uncomment the next line to be consistent with the definition of the goal function
3065 : 547105 : etmp = etmp + (energy_diff(itime)-emu)**2* weights(itime)/(sqomega(itime)**(1.0/2.0))
3066 : : ! uncomment the next get a measure in Ha instead of Ha^2
3067 : : ! etmp = etmp + abs(energy_diff(itime)-emu)
3068 : : ! Fill forces
3069 [ + + ]: 22431305 : do ia=1,natom
3070 [ + + ]: 88083905 : do mu=1,3
3071 : : fmu = zero
3072 [ + + ]: 460727040 : do icoeff=1,ncoeff_fit
3073 : 395074440 : icoeff_tmp = list_coeffs(icoeff)
3074 : 460727040 : fmu = fmu + coefficients(icoeff)*fcart_coeffs(mu,ia,icoeff_tmp,itime)
3075 : : end do
3076 : 87536800 : ftmp = ftmp + (fcart_diff(mu,ia,itime)-fmu)**2* weights(itime)
3077 : : end do !End loop dir
3078 : : end do !End loop natom
3079 [ + + ]: 3841449 : do mu=1,6
3080 : : smu = zero
3081 [ + + ]: 23036352 : do icoeff=1,ncoeff_fit
3082 : 19753722 : icoeff_tmp = list_coeffs(icoeff)
3083 : 23036352 : smu = smu + coefficients(icoeff)*strten_coeffs(mu,itime,icoeff_tmp)
3084 : : end do
3085 : 3829735 : stmp = stmp + sqomega(itime)* weights(itime) *(strten_diff(mu,itime)-smu)**2
3086 : : end do !End loop stress dir
3087 : : end do ! End loop time
3088 : :
3089 : 11714 : gf_value(1) = ffact*ftmp + sfact*stmp !+ efact*etmp !Stresses + Forces
3090 : 11714 : gf_value(2) = ffact*ftmp ! only Forces
3091 : 11714 : gf_value(3) = sfact*stmp ! only Stresses
3092 : 11714 : gf_value(4) = efact*etmp !abs(Energy)
3093 : :
3094 : 11714 : end subroutine fit_polynomial_coeff_computeGF
3095 : : !!***
3096 : :
3097 : :
3098 : : !> Computes force and stress contributions for all polynomial coefficients
3099 : : !>
3100 : : !> This routine:
3101 : : !> - Evaluates energy, force and stress terms for each coefficient
3102 : : !> - Handles displacement-displacement and displacement-strain couplings
3103 : : !> - Applies periodic boundary conditions for supercell calculations
3104 : : !> - Computes derivatives needed for force and stress evaluations
3105 : : !> - Includes strain-dependent volume corrections
3106 : : !> - Supports parallel evaluation across processes
3107 : : !>
3108 : : !> Based on equations 11 & 12 in PRB95,094115 (2017) [[cite:Escorihuela-Sayalero2017]]
3109 : : !>
3110 : : !> @param coefficients Input polynomial coefficients
3111 : : !> @param du_delta Changes in displacements due to strain
3112 : : !> @param displacement Atomic displacements
3113 : : !> @param energy_out Output energy terms for each coefficient
3114 : : !> @param fcart_out Output force terms for each coefficient
3115 : : !> @param natom_sc Number of atoms in supercell
3116 : : !> @param natom_uc Number of atoms in unit cell
3117 : : !> @param ncoeff_max Maximum number of coefficients
3118 : : !> @param ntime Number of configurations
3119 : : !> @param sc_size Supercell dimensions
3120 : : !> @param strain Strain values
3121 : : !> @param strten_out Output stress terms for each coefficient
3122 : : !> @param ucvol Unit cell volumes
3123 : : !> @param coeffs List of coefficients to evaluate
3124 : : !> @param ncoeff Number of coefficients to evaluate
3125 : 11798 : subroutine fit_polynomial_coeff_getFS(coefficients,du_delta,displacement,energy_out,fcart_out,&
3126 : 11798 : & natom_sc,natom_uc,ncoeff_max,ntime,sc_size,strain,strten_out,&
3127 : 11798 : & ucvol,coeffs,ncoeff)
3128 : :
3129 : : implicit none
3130 : :
3131 : : !Arguments ------------------------------------
3132 : : !scalars
3133 : : integer,intent(in) :: natom_sc,natom_uc,ncoeff_max,ntime
3134 : : integer,intent(in) :: ncoeff
3135 : : !arrays
3136 : : integer,intent(in) :: sc_size(3)
3137 : : integer,intent(in) :: coeffs(ncoeff_max)
3138 : : real(dp),intent(in) :: du_delta(6,3,natom_sc,ntime)
3139 : : real(dp),intent(in) :: displacement(3,natom_sc,ntime)
3140 : : real(dp),intent(in) :: strain(6,ntime),ucvol(ntime)
3141 : : real(dp),intent(out):: energy_out(ncoeff,ntime)
3142 : : real(dp),intent(out) :: fcart_out(3,natom_sc,ncoeff,ntime)
3143 : : real(dp),intent(out) :: strten_out(6,ntime,ncoeff)
3144 : : type(polynomial_coeff_type), intent(in) :: coefficients(ncoeff_max)
3145 : : !Local variables-------------------------------
3146 : : !scalar
3147 : : integer :: i1,i2,i3,ia1,ia2,ib1,ib2,ii,icell,icoeff,icoeff_tmp
3148 : : integer :: idir1,idir2,idisp1,idisp2,idisp1_strain,idisp2_strain
3149 : : integer :: iterm,itime,ndisp,ndisp_tot,nstrain,power_disp,power_strain
3150 : : real(dp):: disp1,disp2,tmp1,tmp2,tmp3,weight
3151 : : !arrays
3152 : : integer :: cell_atoma1(3),cell_atoma2(3)
3153 : : integer :: cell_atomb1(3),cell_atomb2(3)
3154 : :
3155 : : ! *************************************************************************
3156 : :
3157 : :
3158 : : !1-Get forces and stresses from the model
3159 : : ! Initialisation of variables
3160 [ + + + + : 89048517 : fcart_out(:,:,:,:) = zero
+ + + + ]
3161 [ + + + + : 3870867 : strten_out(:,:,:) = zero
+ + ]
3162 [ + + + + ]: 1110917 : energy_out(:,:) = zero
3163 : 11798 : icell = 0; ib1=0; ia1=0
3164 [ + + ]: 35394 : do i1=1,sc_size(1)
3165 [ + + ]: 82586 : do i2=1,sc_size(2)
3166 [ + + ]: 165172 : do i3=1,sc_size(3)
3167 : 94384 : ii = icell*natom_uc
3168 : 94384 : icell = icell + 1
3169 : : ! Loop over configurations
3170 [ + + ]: 4537648 : do itime=1,ntime
3171 : : ! Loop over coefficients
3172 [ + + ]: 8887336 : do icoeff_tmp=1,ncoeff
3173 : 4396880 : icoeff = coeffs(icoeff_tmp)
3174 : : ! Loop over terms of this coefficient
3175 [ + + ]: 129038056 : do iterm=1,coefficients(icoeff)%nterm
3176 : 120245104 : ndisp = coefficients(icoeff)%terms(iterm)%ndisp
3177 : 120245104 : nstrain = coefficients(icoeff)%terms(iterm)%nstrain
3178 : 120245104 : ndisp_tot = ndisp + nstrain
3179 : : ! Set the weight of this term
3180 : 120245104 : weight =coefficients(icoeff)%terms(iterm)%weight
3181 : 120245104 : tmp1 = one
3182 : : ! Loop over displacement and strain
3183 [ + + ]: 527428432 : do idisp1=1,ndisp_tot
3184 : :
3185 : : ! Set to one the acculation of forces and strain
3186 : 407183328 : tmp2 = one
3187 : 407183328 : tmp3 = one
3188 : : ! Strain case idir => -6, -5, -4, -3, -2 or -1
3189 [ + + ]: 407183328 : if (idisp1 > ndisp)then
3190 : 48745920 : idisp1_strain = idisp1 - ndisp
3191 : 48745920 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain)
3192 : : ! Get the direction of the displacement or strain
3193 : 48745920 : idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain)
3194 [ + + ]: 48745920 : if(abs(strain(idir1,itime)) > tol10)then
3195 : : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
3196 : 29382624 : tmp1 = tmp1 * (strain(idir1,itime))**power_strain
3197 [ + + ]: 29382624 : if(power_strain > 1) then
3198 : : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
3199 : 45360 : tmp3 = tmp3 * power_strain*(strain(idir1,itime))**(power_strain-1)
3200 : : end if
3201 : : else
3202 : 19363296 : tmp1 = zero
3203 [ + + ]: 19363296 : if(power_strain > 1) then
3204 : 2160 : tmp3 = zero
3205 : : end if
3206 : : end if
3207 : : else
3208 : : ! Set the power_disp of the displacement:
3209 : 358437408 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1)
3210 : : ! Get the direction of the displacement or strain
3211 : 358437408 : idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1)
3212 : : ! Displacement case idir = 1, 2 or 3
3213 : : ! indexes of the cell of the atom a
3214 [ + + ]: 1433749632 : cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1)
3215 [ + - + - : 358437408 : if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then
- + ]
3216 : : ! if the cell is not 0 0 0 we apply PBC:
3217 : 0 : cell_atoma1(1) = i1 + cell_atoma1(1)
3218 : 0 : cell_atoma1(2) = i2 + cell_atoma1(2)
3219 : 0 : cell_atoma1(3) = i3 + cell_atoma1(3)
3220 : 0 : call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3))
3221 : : ! index of the first atom (position in the supercell if the cell is not 0 0 0)
3222 : : ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3223 : : & (cell_atoma1(2)-1)*sc_size(3)*natom_uc+&
3224 : : & (cell_atoma1(3)-1)*natom_uc+&
3225 : 0 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
3226 : : else
3227 : : ! index of the first atom (position in the supercell if the cell is 0 0 0)
3228 : 358437408 : ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1)
3229 : : end if
3230 : :
3231 : : ! indexes of the cell of the atom b (with PBC) same as ia1
3232 [ + + ]: 1433749632 : cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1)
3233 [ + + + + : 358437408 : if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then
+ + ]
3234 : 195199584 : cell_atomb1(1) = i1 + cell_atomb1(1)
3235 : 195199584 : cell_atomb1(2) = i2 + cell_atomb1(2)
3236 : 195199584 : cell_atomb1(3) = i3 + cell_atomb1(3)
3237 : 195199584 : call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3))
3238 : :
3239 : : ! index of the second atom in the (position in the supercell if the cell is not 0 0 0)
3240 : : ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3241 : : & (cell_atomb1(2)-1)*sc_size(3)*natom_uc+&
3242 : : & (cell_atomb1(3)-1)*natom_uc+&
3243 : 195199584 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
3244 : : else
3245 : : ! index of the first atom (position in the supercell if the cell is 0 0 0)
3246 : 163237824 : ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1)
3247 : : end if
3248 : :
3249 : : ! Get the displacement for the both atoms
3250 : 358437408 : disp1 = displacement(idir1,ia1,itime)
3251 : 358437408 : disp2 = displacement(idir1,ib1,itime)
3252 : :
3253 [ + + - + ]: 358437408 : if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then
3254 : : ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z))
3255 : 351505920 : tmp1 = tmp1 * (disp1-disp2)**power_disp
3256 [ + + ]: 351505920 : if(power_disp > 1) then
3257 : : ! Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...))
3258 : 32850240 : tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1)
3259 : : end if
3260 : : else
3261 : 6931488 : tmp1 = zero
3262 [ + + ]: 6931488 : if(power_disp > 1) then
3263 : 368736 : tmp2 = zero
3264 : : end if
3265 : : end if
3266 : : end if
3267 : :
3268 [ + + ]: 1831660032 : do idisp2=1,ndisp_tot
3269 [ + + ]: 1831660032 : if(idisp2 /= idisp1) then
3270 : :
3271 : : ! Strain case
3272 [ + + ]: 1017293376 : if (idisp2 > ndisp)then
3273 : 126657792 : idisp2_strain = idisp2 - ndisp
3274 : 126657792 : idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain)
3275 : : ! Set the power_strain of the strain:
3276 : 126657792 : power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain)
3277 : : ! Accumulate energy forces
3278 : 126657792 : tmp2 = tmp2 * (strain(idir2,itime))**power_strain
3279 : : ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...))
3280 : 126657792 : tmp3 = tmp3 * (strain(idir2,itime))**power_strain
3281 : : ! Atomic displacement case
3282 : : else
3283 : : ! Set the power_disp of the displacement:
3284 : 890635584 : power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2)
3285 : : ! Set the direction of the displacement:
3286 : 890635584 : idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2)
3287 : :
3288 [ + + ]: 3562542336 : cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2)
3289 [ + - + - : 890635584 : if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then
- + ]
3290 : 0 : cell_atoma2(1) = i1 + cell_atoma2(1)
3291 : 0 : cell_atoma2(2) = i2 + cell_atoma2(2)
3292 : 0 : cell_atoma2(3) = i3 + cell_atoma2(3)
3293 : 0 : call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3))
3294 : : ! index of the first atom (position in the supercell and direction)
3295 : : ! if the cell of the atom a is not 0 0 0 (may happen)
3296 : : ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3297 : : & (cell_atoma2(2)-1)*sc_size(3)*natom_uc+&
3298 : : & (cell_atoma2(3)-1)*natom_uc+&
3299 : 0 : & coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
3300 : : else
3301 : : ! index of the first atom (position in the supercell and direction)
3302 : 890635584 : ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2)
3303 : : end if
3304 : :
3305 [ + + ]: 3562542336 : cell_atomb2 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2)
3306 : :
3307 [ + + + + : 890635584 : if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then
+ + ]
3308 : : ! indexes of the cell2 (with PBC)
3309 : 475655760 : cell_atomb2(1) = i1 + cell_atomb2(1)
3310 : 475655760 : cell_atomb2(2) = i2 + cell_atomb2(2)
3311 : 475655760 : cell_atomb2(3) = i3 + cell_atomb2(3)
3312 : 475655760 : call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3))
3313 : :
3314 : : ! index of the second atom in the (position in the supercell)
3315 : : ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+&
3316 : : & (cell_atomb2(2)-1)*sc_size(3)*natom_uc+&
3317 : : & (cell_atomb2(3)-1)*natom_uc+&
3318 : 475655760 : & coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
3319 : : else
3320 : 414979824 : ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2)
3321 : : end if
3322 : :
3323 : 890635584 : disp1 = displacement(idir2,ia2,itime)
3324 : 890635584 : disp2 = displacement(idir2,ib2,itime)
3325 : :
3326 : 890635584 : tmp2 = tmp2 * (disp1-disp2)**power_disp
3327 : 890635584 : tmp3 = tmp3 * (disp1-disp2)**power_disp
3328 : :
3329 : : end if
3330 : : end if
3331 : : end do
3332 : :
3333 [ + + ]: 527428432 : if(idisp1 > ndisp)then
3334 : : ! Accumule stress tensor
3335 : : strten_out(idir1,itime,icoeff_tmp) = strten_out(idir1,itime,icoeff_tmp) + &
3336 : 48745920 : & weight * tmp3
3337 : : else
3338 : : ! Accumule forces
3339 : 358437408 : fcart_out(idir1,ia1,icoeff_tmp,itime)=fcart_out(idir1,ia1,icoeff_tmp,itime)+weight*tmp2
3340 : 358437408 : fcart_out(idir1,ib1,icoeff_tmp,itime)=fcart_out(idir1,ib1,icoeff_tmp,itime)-weight*tmp2
3341 : : end if
3342 : : end do
3343 : :
3344 : : ! accumule energy
3345 : 124641984 : energy_out(icoeff_tmp,itime) = energy_out(icoeff_tmp,itime) + weight * tmp1
3346 : :
3347 : : end do!End do iterm
3348 : : end do!End do coeff
3349 : : end do!End time
3350 : : end do!End do i3
3351 : : end do!End do i2
3352 : : end do!End do i1
3353 : :
3354 : : ! multiply by -1
3355 [ + + + + : 89048517 : fcart_out(:,:,:,:) = -1 * fcart_out(:,:,:,:)
+ + + + ]
3356 : :
3357 : : !ADD stress due to forces on atoms and variation of disp with strain
3358 [ + + ]: 23597 : do icoeff=1,ncoeff
3359 [ + + ]: 573207 : do itime=1,ntime
3360 [ + + ]: 22545809 : do ia1=1,natom_sc
3361 [ + + ]: 88487210 : do idir1=1,3
3362 [ + + ]: 483656800 : do idir2=1,6
3363 : : strten_out(idir2,itime,icoeff) = strten_out(idir2,itime,icoeff) - &
3364 : 461672400 : & du_delta(idir2,idir1,ia1,itime)*fcart_out(idir1,ia1,icoeff,itime)
3365 : : end do
3366 : : end do
3367 : : end do
3368 : : end do
3369 : : end do
3370 : :
3371 [ + + ]: 561307 : do itime=1, ntime
3372 [ + + ]: 3858361 : do idir1=1,6
3373 [ + + ]: 3846563 : if (idir1 < 4) then
3374 [ + + ]: 3297357 : strten_out(idir1,itime,:) = strten_out(idir1,itime,:) * (1 + strain(idir1,itime))/ucvol(itime)
3375 : : else
3376 [ + + ]: 3297357 : strten_out(idir1,itime,:) = strten_out(idir1,itime,:) * (1 - strain(idir1,itime)**2)/ucvol(itime)
3377 : : end if
3378 : : end do
3379 : : end do
3380 : :
3381 : 11798 : end subroutine fit_polynomial_coeff_getFS
3382 : :
3383 : :
3384 : :
3385 : :
3386 : :
3387 : : !> Computes mean squared deviations between model and reference data
3388 : : !>
3389 : : !> This routine:
3390 : : !> - Evaluates the effective potential for all configurations
3391 : : !> - Computes energy, force and stress differences vs reference
3392 : : !> - Calculates weighted mean squared deviations
3393 : : !> - Handles optional printing of detailed comparisons
3394 : : !> - Supports electronic model evaluations if requested
3395 : : !> - Can output phonon-specific analysis data
3396 : : !>
3397 : : !> @param eff_pot Effective potential to evaluate
3398 : : !> @param hist Reference history with DFT data
3399 : : !> @param mse Output mean squared error for energies
3400 : : !> @param msef Output mean squared error for forces
3401 : : !> @param mses Output mean squared error for stresses
3402 : : !> @param natom Number of atoms in system
3403 : : !> @param ntime Number of configurations
3404 : : !> @param sqomega Sheppard-like weighting factors
3405 : : !> @param comm MPI communicator
3406 : : !> @param compute_anharmonic Flag to include anharmonic terms
3407 : : !> @param print_file Flag to print detailed comparison files
3408 : : !> @param filename Base name for output files
3409 : : !> @param scup_dtset Optional SCUP dataset parameters
3410 : : !> @param prt_ph Flag to print phonon analysis
3411 : : !> @param weights Optional weights for each configuration
3412 : 56 : subroutine fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,sqomega,comm,&
3413 : : & compute_anharmonic,print_file,filename,scup_dtset,&
3414 : 20 : & prt_ph, weights)
3415 : :
3416 : : implicit none
3417 : :
3418 : : !Arguments ------------------------------------
3419 : : !scalars
3420 : : integer, intent(in) :: natom,ntime,comm
3421 : : real(dp),intent(out):: mse,msef,mses
3422 : : logical,optional,intent(in) :: compute_anharmonic,print_file,prt_ph
3423 : : !arrays
3424 : : real(dp),intent(in) :: sqomega(ntime)
3425 : : type(effective_potential_type),intent(in) :: eff_pot
3426 : : type(abihist),intent(in) :: hist
3427 : : !Strings/Characters
3428 : : character(len=fnlen),optional,intent(in) :: filename
3429 : : type(scup_dtset_type),optional,intent(inout) :: scup_dtset
3430 : : real(dp),optional, intent(in) :: weights(ntime)
3431 : : !Local variables-------------------------------
3432 : : !scalar
3433 : : integer :: ii,ia,mu,unit_energy,unit_stress,itime,master,nproc,my_rank,i
3434 : : !Uncommend for dipdip test
3435 : : integer :: ifirst
3436 : : real(dp):: energy,energy_harm
3437 : : logical :: need_anharmonic,need_print,need_elec_eval,iam_master
3438 : : logical :: need_prt_ph
3439 : : !arrays
3440 : : real(dp):: strten(6),rprimd(3,3)
3441 : 28 : real(dp),allocatable :: fcart(:,:),gred(:,:),xred(:,:)
3442 : : !Strings/Characters
3443 : : character(len=fnlen) :: file_energy, file_stress, file_anh, name_file
3444 : : character(len=500) :: msg
3445 : : !Uncommend for dipdip test
3446 : 28 : type(abihist) :: hist_out
3447 : : character(len=200) :: filename_hist
3448 : :
3449 : 28 : real(dp),allocatable :: weights_tmp(:)
3450 : : ! *************************************************************************
3451 : : !MS Hide SCALE-UP variables
3452 : : ABI_UNUSED(itime)
3453 : :
3454 [ - + + - : 84 : ABI_MALLOC(fcart,(3,natom))
+ - + - -
+ - + ]
3455 [ - + + - : 56 : ABI_MALLOC(gred,(3,natom))
- + - + ]
3456 [ - + + - : 56 : ABI_MALLOC(xred,(3,natom))
- + - + ]
3457 [ - + + - : 84 : ABI_MALLOC(weights_tmp,(ntime))
+ - + - -
+ - + ]
3458 : :
3459 : : !MPI
3460 : 28 : master = 0
3461 : 28 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
3462 : 28 : iam_master = (my_rank == master)
3463 : :
3464 : : !Do some checks
3465 [ - + ]: 28 : if(ntime /= hist%mxhist)then
3466 : 0 : write(msg,'(a)')'ntime is not correct'
3467 : 0 : ABI_BUG(msg)
3468 : : end if
3469 : :
3470 [ - + ]: 28 : if(natom /= size(hist%xred,2)) then
3471 : 0 : write(msg,'(a)')'natom is not correct'
3472 : 0 : ABI_BUG(msg)
3473 : : end if
3474 : :
3475 : 28 : need_anharmonic = .TRUE.
3476 [ + - ]: 28 : if(present(compute_anharmonic))then
3477 : 28 : need_anharmonic = compute_anharmonic
3478 : : end if
3479 : :
3480 : 28 : name_file=''
3481 [ + + ]: 28 : if(present(filename))name_file = filename
3482 : :
3483 : 28 : need_print=.FALSE.
3484 [ + + ]: 28 : if(present(print_file))need_print=print_file
3485 : :
3486 : 28 : need_elec_eval = .FALSE.
3487 [ + + ]: 28 : if(present(scup_dtset))need_elec_eval=scup_dtset%scup_elec_model
3488 : :
3489 : 28 : need_prt_ph=.FALSE.
3490 [ + + ]: 28 : if(present(prt_ph))need_prt_ph=prt_ph
3491 : :
3492 : :
3493 [ + + ]: 28 : if(need_print .and. present(filename))then
3494 : : !MS hist out uncommented for PHONOPY test
3495 : 21 : call abihist_init(hist_out,natom,ntime,.false.,.false.)
3496 [ + - ]: 21 : file_energy=trim(name_file)//'_energy.dat'
3497 : 21 : unit_energy = get_unit()
3498 [ - + ]: 21 : if (open_file(file_energy,msg,unit=unit_energy,form="formatted",&
3499 : : & status="unknown",action="write") /= 0) then
3500 : 0 : ABI_ERROR(msg)
3501 : : end if
3502 : 21 : unit_stress = get_unit()
3503 [ + - ]: 21 : file_stress=trim(name_file)//'_stress.dat'
3504 [ - + ]: 21 : if (open_file(file_stress,msg,unit=unit_stress,form="formatted",&
3505 : : & status="unknown",action="write") /= 0) then
3506 : 0 : ABI_ERROR(msg)
3507 : : end if
3508 [ - + ]: 7 : else if(need_print .and. .not. present(filename))then
3509 : 0 : write(msg,'(3a)')' You asked for printing of the MSD-values',ch10,&
3510 : 0 : & ' without specifying a filename'
3511 : 0 : ABI_ERROR(msg)
3512 : : end if
3513 : :
3514 [ + + ]: 28 : if(present(weights))then
3515 [ + - - + : 579 : weights_tmp = weights
- - + + ]
3516 : : else
3517 [ + + ]: 816 : weights_tmp = one
3518 : : end if
3519 : :
3520 [ + - ]: 28 : file_anh=trim(name_file)//'_anharmonic_terms_energy.dat'
3521 : :
3522 : 28 : mse = zero
3523 : 28 : msef = zero
3524 : 28 : mses = zero
3525 [ + + ]: 1375 : do ii=1,ntime ! Loop over configurations
3526 [ + + + + ]: 216867 : xred(:,:) = hist%xred(:,:,ii)
3527 [ + + + + ]: 17511 : rprimd(:,:) = hist%rprimd(:,:,ii)
3528 : 1347 : itime = ii
3529 : : #if defined DEV_MS_SCALEUP
3530 : : !Pass print options to scale-up
3531 : : itime = ii
3532 : : if(need_elec_eval)then
3533 : : call global_set_parent_iter(itime)
3534 : : ! Set all print options to false.
3535 : : call global_set_print_parameters(geom=.FALSE.,eigvals=.FALSE.,eltic=.FALSE.,&
3536 : : & orbocc=.FALSE.,bands=.FALSE.)
3537 : : if(ii == 1 .or. modulo(ii,scup_dtset%scup_printniter) == 0)then
3538 : : call global_set_print_parameters(scup_dtset%scup_printgeom,scup_dtset%scup_printeigv,scup_dtset%scup_printeltic,&
3539 : : & scup_dtset%scup_printorbocc,scup_dtset%scup_printbands)
3540 : : end if
3541 : : end if
3542 : : #endif
3543 : : call effective_potential_evaluate(eff_pot,energy_harm,fcart,gred,strten,natom,rprimd,&
3544 : : & xred=xred,compute_anharmonic=.False.,verbose=.false.,&
3545 : 1347 : & elec_eval=need_elec_eval)
3546 : :
3547 : : call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,&
3548 : : & xred=xred,compute_anharmonic=need_anharmonic,verbose=.false.,&
3549 : 1347 : & filename=file_anh,elec_eval=need_elec_eval)
3550 : :
3551 [ + + ]: 1347 : if(need_print .and. iam_master)then
3552 : 552 : WRITE(unit_energy ,'(I10,7(F23.14))') ii,hist%etot(ii),energy_harm,energy,&
3553 : 1104 : & abs(hist%etot(ii) - energy_harm),abs(hist%etot(ii) - energy)
3554 : 552 : WRITE(unit_stress,'(I10,12(F23.14))') ii,hist%strten(:,ii),strten(:)
3555 : : end if
3556 : :
3557 : : !MS Uncommented for abihist test
3558 [ - + ]: 1347 : if(need_prt_ph)then
3559 [ # # ]: 0 : if(ii == 1)then
3560 [ # # # # ]: 0 : write(msg,'(a,(80a))') ch10,('-',i=1,80)
3561 : 0 : call wrtout(ab_out,msg,'COLL')
3562 : 0 : call wrtout(std_out,msg,'COLL')
3563 : 0 : write(msg,'(3a)') ch10,'test_prt_ph == 1, write evulation of Model on the TEST-set into ph_test.nc',ch10
3564 : 0 : call wrtout(ab_out,msg,'COLL')
3565 : 0 : call wrtout(std_out,msg,'COLL')
3566 : : endif
3567 [ # # ]: 0 : ifirst=merge(0,1,(ii>1))
3568 : 0 : filename_hist = trim("ph_test.nc")
3569 [ # # # # ]: 0 : hist_out%fcart(:,:,hist_out%ihist) = fcart(:,:)
3570 [ # # ]: 0 : hist_out%strten(:,hist_out%ihist) = strten(:)
3571 : 0 : hist_out%etot(hist_out%ihist) = energy
3572 : 0 : hist_out%entropy(hist_out%ihist) = hist%entropy(ii)
3573 : 0 : hist_out%time(hist_out%ihist) = real(ii,kind=dp)
3574 : : ! call vel2hist(ab_mover%amass,hist,vel,vel_cell)
3575 : 0 : call var2hist(hist%acell(:,ii),hist_out,natom,hist%rprimd(:,:,ii),hist%xred(:,:,ii),.false.)
3576 [ # # ]: 0 : if(iam_master)then
3577 : : call write_md_hist(hist_out,filename_hist,ifirst,ii,natom,1,eff_pot%crystal%ntypat,&
3578 : : & eff_pot%supercell%typat,eff_pot%crystal%amu,eff_pot%crystal%znucl,&
3579 : 0 : & real(100,dp),(/real(100,dp),real(100,dp)/))
3580 : : endif
3581 : : endif!(need_prt_ph)
3582 : :
3583 : 1347 : mse = mse + weights_tmp(ii) * ((hist%etot(ii) - energy))**2/(sqomega(ii)**(1.0/2.0)) !+abs(hist$etot(ii) - energy)
3584 [ + + ]: 55227 : do ia=1,natom ! Loop over atoms
3585 [ + + ]: 216867 : do mu=1,3 ! Loop over cartesian directions
3586 : 215520 : msef = msef + weights_tmp(ii) * (hist%fcart(mu,ia,ii) - fcart(mu,ia))**2
3587 : : end do
3588 : : end do
3589 [ + + ]: 9457 : do mu=1,6 ! Loop over stresses
3590 : 9429 : mses = mses + weights_tmp(ii) * sqomega(ii)*(hist%strten(mu,ii) - strten(mu))**2
3591 : : end do
3592 : : end do ! End loop itime
3593 [ - + ]: 28 : if(need_prt_ph)then
3594 [ # # # # ]: 0 : write(msg,'(a,(80a))') ch10,('-',i=1,80)
3595 : 0 : call wrtout(ab_out,msg,'COLL')
3596 : 0 : call wrtout(std_out,msg,'COLL')
3597 : : endif
3598 : :
3599 : 28 : mse = mse / ntime
3600 : 28 : msef = msef / (3*natom*ntime)
3601 : 28 : mses = mses / (6*ntime)
3602 : :
3603 [ + + ]: 28 : if(need_print)then
3604 : 21 : close(unit_energy)
3605 : 21 : close(unit_stress)
3606 : : end if
3607 : :
3608 : : !MS uncommented for PHONOPY TEST
3609 : 28 : call abihist_free(hist_out)
3610 : :
3611 [ - + ]: 28 : ABI_FREE(fcart)
3612 [ - + ]: 28 : ABI_FREE(gred)
3613 [ - + ]: 28 : ABI_FREE(xred)
3614 [ - + ]: 28 : ABI_FREE(weights_tmp)
3615 : :
3616 [ + + - + : 56 : end subroutine fit_polynomial_coeff_computeMSD
- + - + -
+ - + - +
- + - + -
+ - + -
+ ]
3617 : :
3618 : :
3619 : :
3620 : : !MARCUS_EXPERIMENTAL_SECTION
3621 : : !!****f* m_fit_polynomiaL_coeff/testEffPot
3622 : : !! NAME
3623 : : !! testEffPot
3624 : : !!
3625 : : !! FUNCTION
3626 : : !! Calculate the energy, forces for displacements provided
3627 : : !! in an test-set (input:hist) within a given effective potential
3628 : : !! (input: eff_pot)
3629 : : !! If the test set is from DFT and contains DFT energies and forces
3630 : : !! calculate the Goal Function values and the MSD of the Energy with
3631 : : !! respect to the DFT energies
3632 : : !!
3633 : : !! INPUTS
3634 : : !! eff_pot = effective_potential datatype
3635 : : !! hist = abihist datatype
3636 : : !!
3637 : : !! OUTPUT
3638 : : !!
3639 : : !! SOURCE
3640 : 1 : subroutine fit_polynomial_coeff_testEffPot(eff_pot,hist,master,comm,print_anharmonic,scup_dtset,prt_ph)
3641 : :
3642 : :
3643 : : implicit none
3644 : :
3645 : : !Arguments ------------------------------------
3646 : : !scalars
3647 : : integer,intent(in) :: master,comm
3648 : : integer,optional,intent(in) :: prt_ph
3649 : : !logicals
3650 : : logical,optional,intent(in) :: print_anharmonic
3651 : : !array
3652 : : type(effective_potential_type),intent(inout) :: eff_pot
3653 : : type(abihist),intent(in) :: hist
3654 : : type(scup_dtset_type),optional,intent(inout) :: scup_dtset
3655 : : !Local variables-------------------------------
3656 : : !reals
3657 : : real(dp) :: factor,mse,msef,mses
3658 : 1 : real(dp),allocatable :: sqomega(:),ucvol(:)
3659 : : real(dp),parameter :: HaBohr_eVAng = Ha_eV /Bohr_Ang
3660 : : !scalar
3661 : : integer :: itime,unit_anh
3662 : : integer :: natom,ntime,ncoeff,my_rank
3663 : : !logicals
3664 : : logical :: iam_master, need_print_anharmonic,file_opened,need_prt_ph
3665 : : !strings/characters
3666 : : character(len=fnlen) :: filename
3667 : : character(len=1000) :: message
3668 : : !arrays
3669 : : real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3)
3670 : : ! *************************************************************************
3671 : :
3672 : : !MPI variables
3673 : 1 : my_rank=xmpi_comm_rank(comm)
3674 : 1 : iam_master = (my_rank == master)
3675 : :
3676 : : !Initialisation of optional arguments
3677 : 1 : need_print_anharmonic = .FALSE.
3678 [ + - ]: 1 : if(present(print_anharmonic)) need_print_anharmonic = print_anharmonic
3679 : :
3680 : 1 : need_prt_ph = .FALSE.
3681 [ + - ]: 1 : if(present(prt_ph))then
3682 [ - + ]: 1 : if(prt_ph==1) need_prt_ph=.TRUE.
3683 : : endif
3684 : :
3685 : : !Setting/Allocating other Variables
3686 : 1 : natom = size(hist%xred,2)
3687 : 1 : factor = 1._dp/natom
3688 : 1 : ntime = hist%mxhist
3689 [ - + + - : 3 : ABI_MALLOC(sqomega,(ntime))
+ - + - -
+ - + ]
3690 [ - + + - : 2 : ABI_MALLOC(ucvol,(ntime))
- + - + ]
3691 [ + + ]: 102 : sqomega = zero
3692 : 1 : filename = 'TES_fit_diff'
3693 : 1 : ncoeff = eff_pot%anharmonics_terms%ncoeff
3694 : :
3695 [ + + ]: 102 : do itime=1,ntime
3696 : : ! Compute \Omega^{2} and ucvol for each time
3697 : 101 : call metric(gmet,gprimd,-1,rmet,hist%rprimd(:,:,itime),ucvol(itime))
3698 : : ! Formula: sqomega(itime) = (((ucvol(itime)**(-2.))* ((natom)**(0.5)))**(-1.0/3.0))**2
3699 : : ! Compact form:
3700 : 102 : sqomega(itime) = ((ucvol(itime)**(4.0/3.0)) / ((natom)**(1/3.0)))
3701 : : end do
3702 : :
3703 : :
3704 [ + - ]: 1 : if(need_print_anharmonic) call effective_potential_writeAnhHead(ncoeff,&
3705 : 1 : & filename,eff_pot%anharmonics_terms)
3706 : :
3707 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,&
3708 : : & sqomega,comm,&
3709 : : & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=filename,&
3710 : 1 : & scup_dtset=scup_dtset,prt_ph=need_prt_ph)
3711 : :
3712 : :
3713 : : ! Print the standard deviation after the fit
3714 : 1 : write(message,'(6a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
3715 : : !& ' Mean Standard Deviation values of the effective-potential',ch10,&
3716 : : !& ' with respect to the test-set (meV^2/atm):',&
3717 : : !& ch10,' Energy : ',&
3718 : : !& mse* (Ha_EV*1000)**2 *factor ,ch10,&
3719 : 1 : & ' Goal function values of the effective.potential',ch10,&
3720 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
3721 : 1 : & ' Energy : ',&
3722 : 1 : & (mse)*(HaBohr_eVAng)**2,ch10,&
3723 : 1 : & ' Forces+Stresses : ',&
3724 : 1 : & (msef+mses)*(HaBohr_eVAng)**2,ch10,&
3725 : 1 : & ' Forces : ',&
3726 : 1 : & msef*(HaBohr_eVAng)**2,ch10,&
3727 : 1 : & ' Stresses : ',&
3728 : 2 : & mses*(HaBohr_eVAng)**2,ch10
3729 : 1 : call wrtout(ab_out,message,'COLL')
3730 : 1 : call wrtout(std_out,message,'COLL')
3731 : :
3732 : :
3733 : : !Deallocating
3734 [ - + ]: 1 : ABI_FREE(sqomega)
3735 [ - + ]: 1 : ABI_FREE(ucvol)
3736 : :
3737 : 1 : INQUIRE(FILE='TES_fit_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh)
3738 [ + - ]: 1 : if(file_opened) close(unit_anh)
3739 : :
3740 : :
3741 : 1 : end subroutine fit_polynomial_coeff_testEffPot
3742 : : !!***
3743 : :
3744 : :
3745 : :
3746 : : !> Prints system files needed for external fitting scripts
3747 : : !>
3748 : : !> This routine writes:
3749 : : !> - Born effective charges
3750 : : !> - Dielectric tensor
3751 : : !> - Reference structure and energy
3752 : : !> - Harmonic force constants in XML format
3753 : : !> - Strain tensors
3754 : : !> - Symmetry operations
3755 : : !> - Molecular dynamics trajectory data
3756 : : !>
3757 : : !> The files are organized in a standard format expected by external fitting tools.
3758 : : !> Special care is taken to properly order atoms and handle supercell transformations.
3759 : : !>
3760 : : !> @param eff_pot Effective potential containing system parameters
3761 : : !> @param hist History containing MD trajectory data
3762 : 0 : subroutine fit_polynomial_printSystemFiles(eff_pot,hist)
3763 : :
3764 : : implicit none
3765 : :
3766 : : !Arguments ------------------------------------
3767 : : !scalars
3768 : : !arrays
3769 : : type(effective_potential_type), intent(in) :: eff_pot
3770 : : type(abihist),intent(in) :: hist
3771 : : !Local variables-------------------------------
3772 : : !scalar
3773 : : integer :: ia,ib,ib1,ii,jj,irpt,kk,ll,mu,nu,nstep,nshift
3774 : : integer :: natom_uc
3775 : : integer :: unit_born=22,unit_epsiloninf=23,unit_md=24
3776 : : integer :: unit_harmonic=25,unit_ref=26,unit_strain=27,unit_sym=28
3777 : : !arrays
3778 : 0 : integer,allocatable :: typat_order(:),typat_order_uc(:)
3779 : : integer, dimension(3) :: A,ncell
3780 : 0 : real(dp), allocatable :: xcart(:,:),fcart(:,:)
3781 : : character(len=500) :: msg
3782 : 0 : type(supercell_type) :: supercell
3783 : : ! *************************************************************************
3784 : :
3785 : : !Create new supercell corresponding to the MD
3786 : 0 : ncell = (/2,2,2/)
3787 : : call supercell%init(eff_pot%crystal%natom, (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/),&
3788 : : eff_pot%crystal%rprimd,eff_pot%crystal%typat,&
3789 [ # # ]: 0 : eff_pot%crystal%xcart,eff_pot%crystal%znucl)
3790 : :
3791 : : !allocation of array
3792 [ # # # # : 0 : ABI_MALLOC(xcart,(3,supercell%natom))
# # # # #
# # # ]
3793 [ # # # # : 0 : ABI_MALLOC(fcart,(3,supercell%natom))
# # # # ]
3794 [ # # # # : 0 : ABI_MALLOC(typat_order,(supercell%natom))
# # # # #
# # # ]
3795 [ # # # # : 0 : ABI_MALLOC(typat_order_uc,(eff_pot%crystal%natom))
# # # # #
# # # ]
3796 : :
3797 : 0 : A = (/ 2, 3, 1/)
3798 : :
3799 [ # # ]: 0 : nshift = product(ncell)
3800 : 0 : natom_uc = eff_pot%crystal%natom
3801 : : !Fill the typat_order array:
3802 : : !In the fit script the atom must be in the order 11111 222222 33333 ..
3803 : : !and the order of the atom can not be change in the fit script,
3804 : : !we transform into the format of the script
3805 : 0 : ib = 1
3806 : 0 : ib1= 1
3807 [ # # ]: 0 : do ii=1,eff_pot%crystal%ntypat
3808 : 0 : jj = A(ii)
3809 [ # # ]: 0 : do kk=1,natom_uc
3810 [ # # ]: 0 : if(supercell%typat(kk)==jj)then
3811 : 0 : typat_order_uc(ib1) = kk
3812 : 0 : ib1 = ib1 + 1
3813 [ # # ]: 0 : do ll=1,nshift
3814 : 0 : ia = (ll-1)*natom_uc + kk
3815 : 0 : typat_order(ib) = ia
3816 : 0 : ib = ib + 1
3817 : : end do
3818 : : end if
3819 : : end do
3820 : : end do
3821 : :
3822 : : ! BORN CHARGES FILE
3823 [ # # ]: 0 : if (open_file('system/Born_Charges',msg,unit=unit_born,form="formatted",&
3824 : : & status="replace",action="write") /= 0) then
3825 : 0 : ABI_ERROR(msg)
3826 : : end if
3827 [ # # ]: 0 : do ii=1,eff_pot%crystal%ntypat
3828 : 0 : jj = A(ii)
3829 [ # # ]: 0 : do ia=1,eff_pot%crystal%natom
3830 [ # # ]: 0 : if(eff_pot%crystal%typat(ia)==jj)then
3831 : 0 : write(unit_born,'(i2,a,1F10.5)') ia," ",eff_pot%crystal%amu(eff_pot%crystal%typat(ia))
3832 [ # # ]: 0 : do mu=1,3
3833 : 0 : WRITE(unit_born,'(a,3(F23.14))') " ",eff_pot%harmonics_terms%zeff(:,mu,ia)
3834 : : end do
3835 : : end if
3836 : : end do
3837 : : end do
3838 : :
3839 : : !DIELECTRIC TENSOR FILE
3840 [ # # ]: 0 : if (open_file('system/Dielectric_Tensor',msg,unit=unit_epsiloninf,form="formatted",&
3841 : : & status="replace",action="write") /= 0) then
3842 : 0 : ABI_ERROR(msg)
3843 : : end if
3844 [ # # ]: 0 : do mu=1,3
3845 : 0 : WRITE(unit_epsiloninf,'(3(F23.14))') eff_pot%harmonics_terms%epsilon_inf(:,mu)
3846 : : end do
3847 : :
3848 : :
3849 : : !REFERENCE STRUCTURE FILE
3850 [ # # ]: 0 : if (open_file('system/Reference_structure',msg,unit=unit_ref,form="formatted",&
3851 : : & status="replace",action="write") /= 0) then
3852 : 0 : ABI_ERROR(msg)
3853 : : end if
3854 : :
3855 : 0 : write(unit_ref,'("Energy (Hartree)")')
3856 : 0 : write(unit_ref,'("================")')
3857 : 0 : write(unit_ref,'(F23.14)') (hist%etot(1)/nshift)
3858 : 0 : write(unit_ref,'("")')
3859 : 0 : write(unit_ref,'("Cell vectors")')
3860 : 0 : write(unit_ref,'("============")')
3861 [ # # ]: 0 : do jj=1,3
3862 [ # # ]: 0 : write(unit_ref,'(3(F22.14))') (supercell%rprimd(:,jj))
3863 : : end do
3864 : :
3865 : 0 : write(unit_ref,'("")')
3866 : 0 : write(unit_ref,'("Atomic positions (Bohr radius)")')
3867 : 0 : write(unit_ref,'("==============================")')
3868 : :
3869 [ # # ]: 0 : do ia=1,supercell%natom
3870 : 0 : write(unit_ref,'(3(F23.14))') supercell%xcart(:,typat_order(ia))
3871 : : end do
3872 : :
3873 : : !Harmonic XML file
3874 [ # # ]: 0 : if (open_file('system/harmonic.xml',msg,unit=unit_harmonic,form="formatted",&
3875 : : & status="replace",action="write") /= 0) then
3876 : 0 : ABI_ERROR(msg)
3877 : : end if
3878 : :
3879 : : !Write header
3880 : 0 : write(unit_harmonic,'("<?xml version=""1.0"" ?>")')
3881 : 0 : write(unit_harmonic,'("<name>")')
3882 : :
3883 [ # # ]: 0 : do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt
3884 [ # # # # : 0 : if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol9)) then
# # # # #
# # # ]
3885 : 0 : write(unit_harmonic,'(" <local_force_constant units=""hartree/bohrradius**2"">")')
3886 : 0 : write(unit_harmonic,'(" <data>")')
3887 [ # # ]: 0 : do ia=1,eff_pot%crystal%natom
3888 [ # # ]: 0 : do mu=1,3
3889 [ # # ]: 0 : do ib=1,eff_pot%crystal%natom
3890 [ # # ]: 0 : do nu=1,3
3891 : : write(unit_harmonic,'(F22.14)', advance="no")&
3892 : : & (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,typat_order_uc(ia),&
3893 : 0 : & nu,typat_order_uc(ib),irpt))
3894 : : end do
3895 : : end do
3896 : 0 : write(unit_harmonic,'(a)')''
3897 : : end do
3898 : : end do
3899 : 0 : write(unit_harmonic,'(" </data>")')
3900 : 0 : write(unit_harmonic,'(" <cell>")')
3901 [ # # ]: 0 : write(unit_harmonic,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt))
3902 : 0 : write(unit_harmonic,'(" </cell>")')
3903 : 0 : write(unit_harmonic,'(" </local_force_constant>")')
3904 : : end if
3905 : : end do
3906 : 0 : write(unit_harmonic,'("</name>")')
3907 : :
3908 : : !STRAIN FILE
3909 [ # # ]: 0 : if (open_file('system/Strain_Tensor',msg,unit=unit_strain,form="formatted",&
3910 : : & status="replace",action="write") /= 0) then
3911 : 0 : ABI_ERROR(msg)
3912 : : end if
3913 [ # # # # ]: 0 : write(unit_strain,'(6(F23.14))') (eff_pot%harmonics_terms%elastic_constants)
3914 : :
3915 : : ! SYM FILE
3916 [ # # ]: 0 : if (open_file('system/symmetry_operations',msg,unit=unit_sym,form="formatted",&
3917 : : & status="replace",action="write") /= 0) then
3918 : 0 : ABI_ERROR(msg)
3919 : : end if
3920 : : 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)",&
3921 : : & " (y,-z,-x) (x,-z,-y) (y,-x,-z) (z,-y,-x) (-x,y,-z) (-z,x,-y) (-y,z,-x) (-x,z,-y)",&
3922 : : & " (-y,x,-z) (-z,y,-x) (-x,-y,z) (-z,-x,y) (-y,-z,x) (-x,-z,y) (-y,-x,z) (-z,-y,x)",&
3923 : : & " (-x,-y,-z) (-z,-x,-y) (-y,-z,-x) (-x,-z,-y) (-y,-x,-z) (-z,-y,-x) (-x,y,z)",&
3924 : : & " (-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)",&
3925 : 0 : & " (z,-y,x) (x,y,-z) (z,x,-y) (y,z,-x) (x,z,-y) (y,x,-z) (z,y,-x)")')
3926 : :
3927 : :
3928 : : !MD file
3929 : 0 : nstep = hist%mxhist
3930 [ # # ]: 0 : if (open_file('system/Molecular_dynamic',msg,unit=unit_md,form="formatted",&
3931 : : & status="replace",action="write") /= 0) then
3932 : 0 : ABI_ERROR(msg)
3933 : : end if
3934 [ # # ]: 0 : do ii=1,nstep
3935 : 0 : write(unit_md,'(I5)') ii-1
3936 : 0 : write(unit_md,'(F22.14)') hist%etot(ii)/nshift
3937 [ # # ]: 0 : do jj=1,3
3938 [ # # ]: 0 : write(unit_md,'(3(F22.14))') (hist%rprimd(:,jj,ii))
3939 : : end do
3940 : : ! Set xcart and fcart for this step
3941 : : call xred2xcart(supercell%natom,hist%rprimd(:,:,ii),&
3942 : 0 : & xcart,hist%xred(:,:,ii))
3943 : :
3944 [ # # # # ]: 0 : fcart(:,:) = hist%fcart(:,:,ii)
3945 : :
3946 [ # # ]: 0 : do ia=1,supercell%natom
3947 : 0 : write(unit_md,'(3(E22.14),3(E22.14))') xcart(:,typat_order(ia)),fcart(:,typat_order(ia))
3948 : : end do
3949 : 0 : write(unit_md,'(6(E22.14))') hist%strten(:,ii)
3950 : : end do
3951 : :
3952 : : !Close files
3953 : 0 : close(unit_ref)
3954 : 0 : close(unit_born)
3955 : 0 : close(unit_harmonic)
3956 : 0 : close(unit_epsiloninf)
3957 : 0 : close(unit_md)
3958 : 0 : close(unit_strain)
3959 : 0 : close(unit_sym)
3960 : :
3961 : : !Deallocation array
3962 [ # # ]: 0 : ABI_FREE(typat_order)
3963 [ # # ]: 0 : ABI_FREE(typat_order_uc)
3964 [ # # ]: 0 : ABI_FREE(xcart)
3965 [ # # ]: 0 : ABI_FREE(fcart)
3966 : 0 : call supercell%free()
3967 : :
3968 [ # # # # : 0 : end subroutine fit_polynomial_printSystemFiles
# # # # #
# # # #
# ]
3969 : :
3970 : :
3971 : : !> Recursively generates combinations of numbers with specific conditions
3972 : : !>
3973 : : !> This routine:
3974 : : !> - Generates combinations of numbers from 1 to n_max of length m_max
3975 : : !> - Ensures each number in a combination is greater than previous one
3976 : : !> - Can either count combinations or store them based on compute flag
3977 : : !> - Uses recursive approach to build combinations
3978 : : !>
3979 : : !> @param i Current index/counter for output combinations
3980 : : !> @param m Current position in combination being built
3981 : : !> @param m_max Maximum length of combinations
3982 : : !> @param n_max Maximum number to use in combinations
3983 : : !> @param list Working array for current combination
3984 : : !> @param list_out Output array for storing combinations
3985 : : !> @param size Size of output array
3986 : : !> @param compute Whether to store combinations or just count them
3987 : 0 : recursive subroutine genereList(i,m,m_max,n_max,list,list_out,size,compute)
3988 : :
3989 : : implicit none
3990 : :
3991 : : !Arguments ---------------------------------------------
3992 : : !scalar
3993 : : integer, intent(in) :: m_max,n_max,m,size
3994 : : integer, intent(inout) :: i
3995 : : logical,intent(in) :: compute
3996 : : !arrays
3997 : : integer, intent(out) :: list(m_max),list_out(size,m_max)
3998 : : !Local variables ---------------------------------------
3999 : : !scalar
4000 : : integer n
4001 : : !arrays
4002 : :
4003 : : ! *************************************************************************
4004 [ # # ]: 0 : if (m > m_max) then
4005 : 0 : i = i + 1
4006 [ # # # # ]: 0 : if(compute)list_out(i,:) = list(:)
4007 : : else
4008 [ # # ]: 0 : do n = 1, n_max
4009 [ # # ]: 0 : if (m == 1)then
4010 : 0 : list(m) = n
4011 : 0 : call genereList (i, m + 1,m_max,n_max,list,list_out,size,compute)
4012 [ # # ]: 0 : else if (n > list(m - 1)) then
4013 : 0 : list(m) = n
4014 : 0 : call genereList (i, m + 1,m_max,n_max,list,list_out,size,compute)
4015 : : end if
4016 : : end do
4017 : : end if
4018 : :
4019 : 0 : end subroutine genereList
4020 : : !!***
4021 : :
4022 : : !> Gathers arrays from all MPI processes with varying lengths
4023 : : !>
4024 : : !> This routine:
4025 : : !> - Gathers real values and their corresponding ordering indices
4026 : : !> - Handles arrays of different sizes from different processes
4027 : : !> - Allocates output buffers based on total gathered size
4028 : : !> - Sets up displacement arrays for MPI communication
4029 : : !> - Performs AllGatherv for both values and indices
4030 : : !>
4031 : : !> @param A Input array of values from this process
4032 : : !> @param order Input array of indices from this process
4033 : : !> @param n Size of input arrays from this process
4034 : : !> @param buff_A Output buffer for gathered values
4035 : : !> @param buff_order Output buffer for gathered indices
4036 : : !> @param totsize Total size across all processes
4037 : : !> @param comm MPI communicator
4038 : : !> @param nproc Number of MPI processes
4039 : 86 : subroutine mpigatherv(A, order, n ,buff_A, buff_order, totsize, comm, nproc)
4040 : : integer, intent(inout) :: n, totsize
4041 : : real(dp), intent(inout) :: A(n)
4042 : : integer, intent(inout):: order(n)
4043 : : integer, intent(in) :: comm, nproc
4044 : : real(dp), allocatable, intent(out) :: buff_A(:)
4045 : : integer, allocatable, intent(out) :: buff_order(:)
4046 : 172 : integer :: disps(nproc), sizes(nproc)
4047 : : integer :: ierr, i
4048 : 86 : totsize=n
4049 : 86 : call xmpi_sum(totsize, comm, ierr)
4050 [ - + + - : 258 : ABI_MALLOC(buff_A, (totsize))
+ - + - -
+ - + -
+ ]
4051 [ - + + - : 258 : ABI_MALLOC(buff_order, (totsize))
+ - + - -
+ - + -
+ ]
4052 : 86 : call xmpi_allgather(n, sizes, comm, ierr)
4053 : 86 : disps(1)=0
4054 [ + + ]: 198 : do i=2, nproc
4055 : 198 : disps(i)=disps(i-1)+sizes(i-1)
4056 : : end do
4057 : 86 : call xmpi_allgatherv(A, n, buff_A, sizes, disps, comm, ierr )
4058 : 86 : call xmpi_allgatherv(order, n, buff_order, sizes, disps, comm, ierr )
4059 : 86 : end subroutine mpigatherv
4060 : :
4061 : :
4062 : :
4063 : : end module m_fit_polynomial_coeff
4064 : : !!***
|
