Branch data Line data Source code
1 : :
2 : : !!****m* ABINIT/m_opt_effpot
3 : : !!
4 : : !! NAME
5 : : !! m_opt_effpot
6 : : !!
7 : : !! FUNCTION
8 : : !!
9 : : !!
10 : : !!
11 : : !! COPYRIGHT
12 : : !! Copyright (C) 2010-2026 ABINIT group (AM)
13 : : !! This file is distributed under the terms of the
14 : : !! GNU General Public Licence, see ~abinit/COPYING
15 : : !! or http://www.gnu.org/copyleft/gpl.txt .
16 : : !! For the initials of contributors, see ~abinit/doc/developers/contributors.txt .
17 : : !!
18 : : !! SOURCE
19 : :
20 : :
21 : :
22 : : #if defined HAVE_CONFIG_H
23 : : #include "config.h"
24 : : #endif
25 : :
26 : : #include "abi_common.h"
27 : : !#include "trace.h"
28 : :
29 : :
30 : : module m_opt_effpot
31 : :
32 : : use defs_basis
33 : : use defs_datatypes
34 : : use defs_abitypes
35 : : use m_errors
36 : : use m_abicore
37 : : use m_xmpi
38 : : use m_effective_potential
39 : : use m_effective_potential_file, only : effective_potential_file_mapHistToRef
40 : : use m_fit_data
41 : : use m_fit_polynomial_coeff
42 : : use m_polynomial_coeff
43 : : use m_polynomial_term
44 : : use m_crystal,only : symbols_crystal
45 : :
46 : : implicit none
47 : :
48 : : public :: opt_effpot
49 : : public :: opt_effpotbound
50 : : public :: opt_getHOforterm
51 : : public :: opt_getCombisforterm
52 : : public :: opt_getHoTerms
53 : : public :: opt_getHOstrain
54 : : public :: opt_getHOcrossdisp
55 : : public :: opt_filterdisp
56 : : public :: opt_getSingleDispTerms
57 : : public :: opt_getHOSingleDispTerms
58 : : private :: opt_boundcoeff
59 : : private :: check_to_skip
60 : : !!****
61 : : CONTAINS
62 : :
63 : : !!****f* m_opt_effpot/opt_effpot
64 : : !!
65 : : !! NAME
66 : : !! opt_effpot
67 : : !!
68 : : !! FUNCTION
69 : : !! Optimize Effective Potential by fitting the value of certain
70 : : !! coefficients while keeping the values of the others
71 : : !!
72 : : !! INPUTS
73 : : !! eff_pot<type(effective_potential)> = effective potential
74 : : !!
75 : : !! opt_coeff(opt_ncoeff) = list of terms whose coefficients are to be
76 : : !! optimized
77 : : !!
78 : : !! hist<type(abihist)> = Training set Data(or snapshot of DFT)
79 : : !! comm = MPI communicator
80 : : !!
81 : : !! OUTPUT
82 : : !! eff_pot<type(effective_potential)> = effective potential datatype with new fitted coefficients
83 : : !!
84 : : !! SOURCE
85 : :
86 : 1 : subroutine opt_effpot(eff_pot,opt_ncoeff,opt_coeff,hist,opt_on,opt_factors,comm, fit_weight_T, print_anh)
87 : :
88 : : !Arguments ------------------------------------
89 : : !scalars
90 : : integer,intent(in) :: comm,opt_ncoeff
91 : : type(effective_potential_type),intent(inout) :: eff_pot
92 : : type(abihist),intent(inout) :: hist
93 : : !arrays
94 : : integer,intent(in) :: opt_coeff(opt_ncoeff)
95 : : real(dp),intent(in) :: opt_factors(3)
96 : : real(dp), intent(in):: fit_weight_T
97 : : !Logicals
98 : : logical,intent(in) :: opt_on(3)
99 : : logical,optional,intent(in) :: print_anh
100 : : !Strings
101 : : !Local variables ------------------------------
102 : : !scalars
103 : : integer :: ii, info,natom_sc,ntime,unit_anh1,unit_anh2
104 : : integer :: master,nproc,my_rank
105 : : real(dp) :: factor,mse,msef,mses
106 : : real(dp),parameter :: HaBohr_eVAng = Ha_eV /Bohr_Ang
107 : : !arrays
108 : : integer :: sc_size(3)
109 : 2 : integer :: coeff_inds(opt_ncoeff)
110 : 1 : type(fit_data_type) :: fit_data
111 [ + + + + : 11 : type(polynomial_coeff_type) :: my_coeffs(opt_ncoeff)
+ + + - -
- - - - -
- - - - -
- - - -
- ]
112 : 2 : real(dp) :: coeff_values(opt_ncoeff), coeff_init_values(opt_ncoeff)
113 : 1 : real(dp), allocatable :: energy_coeffs(:,:),fcart_coeffs(:,:,:,:)
114 : 1 : real(dp), allocatable :: strten_coeffs(:,:,:)
115 : 1 : real(dp), allocatable :: my_weights(:)
116 : : !Logicals
117 : : logical :: need_print_anh,file_opened,iam_master
118 : : !Strings
119 : : character(len=1000) :: message
120 : : character(len=1000) :: frmt
121 : : character(len=fnlen) :: fn_bf='before_opt_diff', fn_af='after_opt_diff'
122 : : ! *************************************************************************
123 : : !MPI
124 : 1 : master = 0
125 : 1 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
126 : : iam_master = (my_rank == master)
127 : :
128 : : !Setting/Initializing Variables
129 : 1 : ntime = hist%mxhist
130 : 1 : natom_sc = size(hist%xred,2)
131 : 1 : factor = 1._dp/natom_sc
132 : 1 : need_print_anh =.False.
133 : :
134 : 1 : call get_weight_from_hist(hist, fit_weight_T, ntime, natom_sc, my_weights, comm)
135 : :
136 : :
137 [ - + ]: 1 : if(present(print_anh)) then
138 [ + - ]: 1 : if(print_anh) need_print_anh=.True.
139 : : end if
140 : : !if the number of atoms in reference supercell into effpot is not correct,
141 : : !wrt to the number of atom in the hist, we set map the hist and set the good supercell
142 [ - + ]: 1 : if (natom_sc /= eff_pot%supercell%natom) then
143 : 0 : call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=.TRUE.)
144 : : end if
145 : :
146 : : !we get the size of the supercell in the hist file
147 : : ! This is only valid for diagonal sc_size.
148 : : ! F08: sc_size(:) = int(anint( norm2(eff_pot%supercell%rprimd, dim=2) / &
149 : : ! & norm2(eff_pot%crystal%rprimd, dim=2) ))
150 [ + + ]: 4 : do ii=1,3
151 : : sc_size(ii) = int(anint(sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
152 : : & eff_pot%supercell%rprimd(ii,2)**2+&
153 : : & eff_pot%supercell%rprimd(ii,3)**2) / &
154 : : & sqrt(eff_pot%crystal%rprimd(ii,1)**2+&
155 : : & eff_pot%crystal%rprimd(ii,2)**2+&
156 : 4 : & eff_pot%crystal%rprimd(ii,3)**2)))
157 : : end do
158 : :
159 : :
160 : : !Before the fit, compute constants with fit_data_compute.
161 : : !Compute the strain of each configuration.
162 : : !Compute the displacmeent of each configuration.
163 : : !Compute the variation of the displacement due to strain of each configuration.
164 : : !Compute fixed forces and stresse and get the standard deviation.
165 : : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
166 : 1 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=.FALSE.)
167 : :
168 : :
169 [ - + ]: 1 : if(need_print_anh) call effective_potential_writeAnhHead(eff_pot%anharmonics_terms%ncoeff,&
170 : 0 : & fn_bf,eff_pot%anharmonics_terms)
171 : :
172 : : !Before deleting coefficients calculate MSD of initial model
173 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
174 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
175 : 1 : & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=fn_bf)
176 : :
177 : :
178 : : ! Print the standard devition of initial model
179 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
180 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
181 : 1 : & ' with respect to the training-set before optimization (meV^2/atm):',&
182 : 1 : & ch10,' Energy : ',&
183 : 1 : & mse* (Ha_EV*1000)**2 *factor ,ch10,&
184 : 1 : & ' Goal function values of the effective.potential',ch10,&
185 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
186 : 1 : & ' Forces+Stresses : ',&
187 : 1 : & (msef+mses)*(HaBohr_eVAng)**2,ch10,&
188 : 1 : & ' Forces : ',&
189 : 1 : & msef*(HaBohr_eVAng)**2,ch10,&
190 : 1 : & ' Stresses : ',&
191 : 2 : & mses*(HaBohr_eVAng)**2,ch10
192 : 1 : call wrtout(ab_out,message,'COLL')
193 : 1 : call wrtout(std_out,message,'COLL')
194 : :
195 : :
196 : : ! Write terms to my_coeffs(ii) and zero them in eff_pot
197 [ + + ]: 3 : do ii=1,opt_ncoeff
198 : : !Store indices for later
199 : 2 : coeff_inds(ii) = ii
200 : : !Initialize coefficients for optimizing
201 : : call polynomial_coeff_init(coeff_values(ii),eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%nterm,&
202 : : & my_coeffs(ii), eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%terms, &
203 : 2 : & check=.TRUE.)
204 : :
205 : : !DMSG(my_coeffs(ii)%debug_str)
206 : :
207 : : !Store initial values of coefficients
208 : 2 : coeff_init_values(ii) = eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%coefficient
209 : : !Put them temporarely to zero
210 : 3 : eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%coefficient = zero
211 : : end do
212 : :
213 : : !Before the fit, compute constants with fit_data_compute.
214 : : !And coefficients to be optimized put to zero
215 : : !Conpute the strain of each configuration.
216 : : !Compute the displacmeent of each configuration.
217 : : !Compute the variation of the displacement due to strain of each configuration.
218 : : !Compute fixed forces and stresse and get the standard deviation.
219 : : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
220 : 1 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=.TRUE.)
221 : :
222 : : !After deleting coefficients calculate MSD
223 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
224 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
225 : 1 : & compute_anharmonic=.TRUE.)
226 : :
227 : :
228 : : ! Print the standard deviation after deleting
229 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
230 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
231 : 1 : & ' with respect to the training-set after deleting selected terms (meV^2/atm):',&
232 : 1 : & ch10,' Energy : ',&
233 : 1 : & mse* (Ha_EV*1000)**2 *factor ,ch10,&
234 : 1 : & ' Goal function values of the effective.potential',ch10,&
235 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
236 : 1 : & ' Forces+Stresses : ',&
237 : 1 : & (msef+mses)*(HaBohr_eVAng)**2,ch10,&
238 : 1 : & ' Forces : ',&
239 : 1 : & msef*(HaBohr_eVAng)**2,ch10,&
240 : 1 : & ' Stresses : ',&
241 : 2 : & mses*(HaBohr_eVAng)**2,ch10
242 : 1 : call wrtout(ab_out,message,'COLL')
243 : 1 : call wrtout(std_out,message,'COLL')
244 : :
245 : :
246 : :
247 : : ! Allocate necessary arrays for the fit-data
248 [ - + + - : 4 : ABI_MALLOC(energy_coeffs,(opt_ncoeff,ntime))
- + + - +
- + - - +
- + ]
249 [ - + + - : 5 : ABI_MALLOC(fcart_coeffs,(3,natom_sc,opt_ncoeff,ntime))
- + + - -
+ + - + -
+ - - + -
+ ]
250 [ - + + - : 4 : ABI_MALLOC(strten_coeffs,(6,ntime,opt_ncoeff))
- + + - +
- + - - +
- + ]
251 : : ! Calculate forces and stresses per coefficient, which are to be optimized
252 : : call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,&
253 : : & fit_data%training_set%displacement,&
254 : : & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,&
255 : : & opt_ncoeff,ntime,sc_size,fit_data%training_set%strain,&
256 : 1 : & strten_coeffs,fit_data%training_set%ucvol,coeff_inds,opt_ncoeff)
257 : :
258 : :
259 : : ! call the fit process routine
260 : : ! This routine solves the linear system proposed
261 : : ! by C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
262 : : call fit_polynomial_coeff_solve(coeff_values(1:opt_ncoeff),fcart_coeffs,fit_data%fcart_diff,&
263 : : & energy_coeffs,fit_data%energy_diff,info,&
264 : : & coeff_inds,natom_sc,opt_ncoeff,opt_ncoeff,ntime,&
265 : : & strten_coeffs,fit_data%strten_diff,&
266 : : & fit_data%training_set%sqomega,opt_on,opt_factors, &
267 : 1 : & weights=my_weights)
268 [ + - - + : 1 : if (info /= 0 .and. all(coeff_values < tol16))then
- + - - ]
269 : 0 : write(frmt,*) opt_ncoeff
270 : 0 : write(message, '(2a,'//ADJUSTR(frmt)//'I4,8a)' ) ch10,&
271 : 0 : & ' The attempt to optimize the terms: ', opt_coeff ,ch10,&
272 : 0 : & ' , returned a singular solution', ch10,&
273 : 0 : & ' The terms could not be optimized ',ch10,&
274 : 0 : & ' and the effective potential has not been altered.', ch10,&
275 : 0 : & ' Action: Change training set or coefficients to be optimized.'
276 : 0 : ABI_WARNING(message)
277 [ # # ]: 0 : do ii=1,opt_ncoeff
278 : 0 : eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%coefficient = coeff_init_values(ii)
279 : 0 : call polynomial_coeff_free(my_coeffs(ii))
280 : : end do
281 : : else
282 : : ! Transfer new fitted values to coefficients and write them into effective potential
283 : : ! Deallcoate temporary coefficients my_coeffs
284 [ + + ]: 3 : do ii=1,opt_ncoeff
285 : 2 : eff_pot%anharmonics_terms%coefficients(opt_coeff(ii))%coefficient = coeff_values(ii)
286 : 3 : call polynomial_coeff_free(my_coeffs(ii))
287 : : end do
288 : : !Recalculate MSD of Final Model
289 : :
290 : : !Conpute the strain of each configuration.
291 : : !Compute the displacmeent of each configuration.
292 : : !Compute the variation of the displacement due to strain of each configuration.
293 : : !Compute fixed forces and stresse and get the standard deviation.
294 : : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
295 : 1 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=.TRUE.)
296 : :
297 : :
298 [ - + ]: 1 : if(need_print_anh) call effective_potential_writeAnhHead(eff_pot%anharmonics_terms%ncoeff,&
299 : 0 : & fn_af,eff_pot%anharmonics_terms)
300 : :
301 : : !After optimization of coefficients opt_coeff recalculate MSD
302 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
303 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
304 : 1 : & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=fn_af)
305 : :
306 : :
307 : : ! Print the standard deviation after optimization
308 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
309 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
310 : 1 : & ' with respect to the training-set after optimizing selected terms (meV^2/atm):',&
311 : 1 : & ch10,' Energy : ',&
312 : 1 : & mse* (Ha_EV*1000)**2 *factor ,ch10,&
313 : 1 : & ' Goal function values of the effective.potential',ch10,&
314 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
315 : 1 : & ' Forces+Stresses : ',&
316 : 1 : & (msef+mses)*(HaBohr_eVAng)**2,ch10,&
317 : 1 : & ' Forces : ',&
318 : 1 : & msef*(HaBohr_eVAng)**2,ch10,&
319 : 1 : & ' Stresses : ',&
320 : 2 : & mses*(HaBohr_eVAng)**2,ch10
321 : 1 : call wrtout(ab_out,message,'COLL')
322 : 1 : call wrtout(std_out,message,'COLL')
323 : : end if
324 : :
325 : : !Deallocation of fitting variables
326 [ + - ]: 1 : ABI_SFREE(energy_coeffs)
327 [ + - ]: 1 : ABI_SFREE(fcart_coeffs)
328 [ + - ]: 1 : ABI_SFREE(strten_coeffs)
329 [ + - ]: 1 : ABI_SFREE(my_weights)
330 : :
331 [ - + ]: 1 : if(need_print_anh)then
332 : 0 : INQUIRE(FILE='before_opt_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh1)
333 [ # # ]: 0 : if(file_opened) close(unit_anh1)
334 : 0 : INQUIRE(FILE='after_opt_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh2)
335 [ # # ]: 0 : if(file_opened) close(unit_anh2)
336 : : end if
337 : : ! Deallocate and delete the fit-date
338 : 1 : call fit_data_free(fit_data)
339 [ - + - + : 7 : end subroutine opt_effpot
- + - + -
+ - + - +
- + - + ]
340 : : !!***
341 : :
342 : : !!****f* m_opt_effpot/opt_effpotbound
343 : : !!
344 : : !! NAME
345 : : !! opt_effpotbound
346 : : !!
347 : : !! FUNCTION
348 : : !! Compute and add high order terms to existing odd or negative even anharmonic terms
349 : : !! Fix the coefficient of the added new high order terms to a value such that it
350 : : !! doesn't influence the precision of the existing anharmonic potential with respect
351 : : !! to a relevent training set (ATTENTIION: A user must know what a relevant training set
352 : : !! is for the system he want's to study. Typically something oscillating around its ground-state.)
353 : : !! Finally optimize the coefficients of the orignal anharmonic terms under the presence of the
354 : : !! added high order terms.
355 : : !!
356 : : !! INPUTS
357 : : !! eff_pot: existing effective potential
358 : : !! order: order for which bounding terms are generated
359 : : !! order_ran: ?
360 : : !! bound_EFS:
361 : : !! bound_factors:
362 : : !! bound_penalty:
363 : : !! comm: MPI communicator
364 : : !! fit_weight_T: ?
365 : : !! print_anh: ?
366 : : !! bound_option: 1 or 2. Whether the bounding term should be the same nbody with the term to be bounded. 1. Yes (default), 2. No
367 : : !!
368 : : !!
369 : : !! OUTPUT
370 : : !! eff_pot new effective potential
371 : : !!
372 : : !!
373 : : !! SOURCE
374 : :
375 : 1 : subroutine opt_effpotbound(eff_pot,order_ran,hist,bound_EFS,bound_factors,bound_penalty,comm, fit_weight_T,bound_option, print_anh)
376 : :
377 : : !Arguments ------------------------------------
378 : : !scalars
379 : : integer,intent(in) :: comm
380 : : type(effective_potential_type),target,intent(inout) :: eff_pot
381 : : type(abihist),intent(inout) :: hist
382 : : real(dp) :: bound_penalty
383 : : !arrays
384 : : integer,intent(in) :: order_ran(2),bound_EFS(3)
385 : : integer, intent(in) :: bound_option
386 : : real(dp),intent(in) :: bound_factors(3)
387 : : real(dp), intent(in) :: fit_weight_T
388 : : !Logicals
389 : : logical,optional,intent(in) :: print_anh
390 : : !Strings
391 : : !Local variables ------------------------------
392 : : !scalars
393 : : integer :: i,ii,natom_sc,ntime,iterm,nterm, icombi3, nbody_term
394 : : integer :: jterm, ncombi,ncombi1,ncombi2, ncombi1_real, ncombi2_real
395 : : integer :: icombi
396 : : integer :: nterm_start,nterm2
397 : : integer :: nproc,my_rank,master
398 : : !1406
399 : : real(dp) :: factor,mse_ini,msef_ini,mses_ini,mse,msef,mses
400 : : real(dp) :: coeff_tmp
401 : : real(dp),parameter :: HaBohr_eVAng = Ha_eV /Bohr_Ang
402 : : !arrays
403 : : integer :: sc_size(3),temp_cntr
404 : 1 : integer,allocatable :: terms(:)
405 : 1 : logical,allocatable :: exists(:)
406 : : logical :: any_exists
407 : 1 : type(fit_data_type) :: fit_data
408 : : real(dp) :: GF_arr(2),coeff_opt(2)
409 : 1 : real(dp), allocatable :: my_weights(:)
410 : : !real(dp), allocatable :: energy_coeffs(:,:),fcart_coeffs(:,:,:,:)
411 : : !real(dp), allocatable :: strten_coeffs(:,:,:)
412 : : !1406 strain_temrs_tmp
413 : 1 : type(polynomial_coeff_type),target,allocatable :: my_coeffs(:),my_coeffs_tmp(:)
414 : 1 : type(polynomial_coeff_type),allocatable :: singledisp_terms(:),HOsingledisp_terms(:)
415 : 1 : type(polynomial_coeff_type),allocatable :: HOcrossdisp_terms(:)
416 : : !Logicals
417 : : logical :: need_print_anh ! MARCUS FOR THE MOMENT PRINT NO FILES
418 : : logical :: to_skip,iam_master
419 : : !Strings
420 : : character(len=5),allocatable :: symbols(:)
421 : : character(len=200):: name
422 : : character(len=1000) :: message
423 : : character(len=fnlen) :: fn_bf='before_opt_diff'!, fn_af='after_opt_diff'
424 : : ! types
425 : 1 : type(SymPairs_t) :: sympairs
426 : : !*************************************************************************
427 : : !MPI variables
428 : 1 : master = 0
429 : 1 : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
430 : : iam_master = (my_rank == master)
431 : :
432 : : ! Say hello to the world!
433 : 1 : write(message, '(3a)' )'-Start Bound optimization of Anharmonic Potential ',ch10
434 : 1 : call wrtout(ab_out,message,'COLL')
435 : 1 : call wrtout(std_out,message,'COLL')
436 : :
437 : : !Setting/Initializing Variables
438 : 1 : ntime = hist%mxhist
439 : :
440 : :
441 : 1 : natom_sc = size(hist%xred,2)
442 : 1 : factor = 1._dp/natom_sc
443 : 1 : nterm =eff_pot%anharmonics_terms%ncoeff
444 : 1 : need_print_anh=.False.
445 [ + - ]: 1 : if(present(print_anh))then
446 [ # # ]: 0 : if(print_anh) need_print_anh = .True.
447 : : endif
448 : 1 : call get_weight_from_hist(hist, fit_weight_T, ntime, natom_sc, my_weights, comm )
449 : :
450 [ - + + - : 3 : ABI_MALLOC(symbols,(eff_pot%crystal%natom))
+ - + - -
+ - + ]
451 [ - + + - : 3 : ABI_MALLOC(terms,(nterm))
+ - + - -
+ - + ]
452 : : call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,&
453 : 1 : & symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl)
454 : :
455 : :
456 : : !if the number of atoms in reference supercell into effpot is not correct,
457 : : !wrt to the number of atom in the hist, we set map the hist and set the good supercell
458 [ - + ]: 1 : if (natom_sc /= eff_pot%supercell%natom) then
459 : 0 : call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=.TRUE.)
460 : : end if
461 : :
462 : : !Check if input of order is correct
463 : : ! TODO write error message here
464 [ + + + - : 3 : if (any(mod(order_ran,2) /= 0)) return
+ - ]
465 : :
466 : : !we get the size of the supercell in the hist file
467 [ + + ]: 4 : do ii=1,3
468 : : sc_size(ii) = int(anint(sqrt(eff_pot%supercell%rprimd(ii,1)**2+&
469 : : & eff_pot%supercell%rprimd(ii,2)**2+&
470 : : & eff_pot%supercell%rprimd(ii,3)**2) / &
471 : : & sqrt(eff_pot%crystal%rprimd(ii,1)**2+&
472 : : & eff_pot%crystal%rprimd(ii,2)**2+&
473 : 4 : & eff_pot%crystal%rprimd(ii,3)**2)))
474 : : end do
475 : :
476 : :
477 : : !Before the fit, compute constants with fit_data_compute.
478 : : !Conpute the strain of each configuration.
479 : : !Compute the displacmeent of each configuration.
480 : : !Compute the variation of the displacement due to strain of each configuration.
481 : : !Compute fixed forces and stresse and get the standard deviation.
482 : : !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]].
483 : : !call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=.FALSE.)
484 : 1 : call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=.FALSE.)
485 : :
486 [ - + ]: 1 : if(need_print_anh) call effective_potential_writeAnhHead(eff_pot%anharmonics_terms%ncoeff,&
487 : 0 : & fn_bf,eff_pot%anharmonics_terms)
488 : :
489 : : !Before adding bound coefficients calculate MSD of initial model
490 : : !MS FOR THE MOMENT PRINT NO FILE
491 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse_ini,msef_ini,mses_ini,&
492 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
493 : 1 : & compute_anharmonic=.TRUE.,print_file=.FALSE.)
494 : :
495 : :
496 : : ! Print the standard devition of initial model
497 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
498 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
499 : 1 : & ' with respect to the training-set before attempted bounding (meV^2/atm):',&
500 : 1 : & ch10,' Energy : ',&
501 : 1 : & mse_ini* (Ha_EV*1000)**2 *factor ,ch10,&
502 : 1 : & ' Goal function values of the effective.potential',ch10,&
503 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
504 : 1 : & ' Forces+Stresses : ',&
505 : 1 : & (msef_ini+mses_ini)*(HaBohr_eVAng)**2,ch10,&
506 : 1 : & ' Forces : ',&
507 : 1 : & msef_ini*(HaBohr_eVAng)**2,ch10,&
508 : 1 : & ' Stresses : ',&
509 : 2 : & mses_ini*(HaBohr_eVAng)**2,ch10
510 : 1 : call wrtout(ab_out,message,'COLL')
511 : 1 : call wrtout(std_out,message,'COLL')
512 : :
513 : : !MS DEV
514 : : ! single displacement terms with second order only.
515 : 1 : call opt_getSingleDispTerms(singledisp_terms,eff_pot%crystal, sc_size,comm)
516 : :
517 : :
518 : : ! create pair list:
519 : 1 : call sympairs%init(eff_pot%crystal, sc_size)
520 : :
521 : : !For the moment order loop commented
522 : : !do iorder=order(1),order(2),2, Order will be done per term
523 : : !Loop over all original terms + 1
524 : : ! + 1 to bound pure strain
525 [ + + ]: 3 : do iterm =1,nterm +1
526 [ + + ]: 2 : if(iterm <=nterm)then
527 : 1 : ncombi1=0
528 : 1 : ncombi2=0
529 : 1 : ncombi1_real=0
530 : 1 : ncombi2_real=0
531 : : !Store for optimization
532 : 1 : terms(iterm) = iterm
533 : : !Message: The world wants to know where we stand Batman
534 : 1 : write(message, '(a,(80a),a)' ) ch10,&
535 [ + - + + ]: 81 : & ('_',ii=1,80),ch10
536 : 1 : call wrtout(ab_out,message,'COLL')
537 : 1 : call wrtout(std_out,message,'COLL')
538 : 1 : write(message,'(2a,I3,a,I3,3a)' )ch10,&
539 [ + - ]: 2 : & ' Check term (',iterm,'/',nterm,'): ', trim(eff_pot%anharmonics_terms%coefficients(iterm)%name),ch10
540 : 1 : call wrtout(ab_out,message,'COLL')
541 : 1 : call wrtout(std_out,message,'COLL')
542 : 1 : to_skip = .FALSE.
543 : 1 : to_skip = check_to_skip(eff_pot%anharmonics_terms%coefficients(iterm))
544 : : !Skip term if it doesn't need bounding
545 [ + - ]: 1 : if(.not. to_skip)then
546 : : associate(term1=>eff_pot%anharmonics_terms%coefficients(iterm)%terms(1))
547 : 1 : nbody_term = term1%get_nbody()
548 : :
549 : : !Get List of high order single Terms for terms
550 : : call opt_getHOSingleDispTerms(eff_pot%anharmonics_terms%coefficients(iterm),&
551 [ - - - - : 1 : & HOsingledisp_terms,symbols,singledisp_terms,order_ran,ncombi1)
- - - - -
- - - - -
- - - - -
- - + ]
552 : :
553 [ - + + - : 4 : ABI_MALLOC(my_coeffs,(size(eff_pot%anharmonics_terms%coefficients)))
+ - + - -
+ - + - +
+ + ]
554 : 1 : my_coeffs=eff_pot%anharmonics_terms%coefficients
555 : 1 : do icombi3=1, size(my_coeffs)
556 : : !DMSG(my_coeffs(icombi3)%debug_str)
557 : : end do
558 : :
559 : : ! then add the single disp terms.
560 [ + - ]: 1 : if (ncombi1>0) then
561 [ + + ]: 3 : do icombi3=1,ncombi1
562 [ + - - + ]: 3 : if (HOsingledisp_terms(icombi3)%terms(1)%get_nbody() == nbody_term .or. bound_option/=1) then
563 : 0 : call coeffs_list_append(my_coeffs,HOsingledisp_terms(icombi3), check=.TRUE.)
564 : : !DMSG(my_coeffs(icombi3)%debug_str)
565 : 0 : ncombi1_real = ncombi1_real + 1
566 : : endif
567 : : end do
568 : : endif
569 : :
570 [ + - ]: 1 : if(allocated(HOsingledisp_terms)) call polynomial_coeff_list_free(HOsingledisp_terms)
571 : :
572 : :
573 : : !Get List of high order cross Terms for term if ndisp > 1
574 [ - + - - : 1 : if(term1%ndisp>1 .or. &
- - ]
575 : : & term1%ndisp /= 0 .and. & ! why >1 or /=0?
576 : : & term1%nstrain /= 0)then
577 : : ! output HOcrossdisp_terms, ncombi2.
578 : 1 : call opt_getHOcrossdisp(HOcrossdisp_terms,ncombi2,eff_pot%anharmonics_terms%coefficients(iterm),order_ran)
579 : : endif
580 : : ! then add the crossdisp terms.
581 [ + - ]: 1 : if(ncombi2 > 0)then
582 [ + + ]: 4 : do icombi3=1,ncombi2
583 : : !call coeffs_list_conc_onsite(my_coeffs, HOcrossdisp_terms(icombi3))
584 [ + + - + ]: 4 : if (HOcrossdisp_terms(icombi3)%terms(1)%get_nbody() == nbody_term .or. bound_option/=1) then
585 : 1 : call coeffs_list_append(my_coeffs,HOcrossdisp_terms(icombi3), check=.TRUE.)
586 : 1 : ncombi2_real = ncombi2_real + 1
587 : : !DMSG(my_coeffs(ncombi2_real)%debug_str)
588 : : endif
589 : : end do
590 : : endif
591 [ - + + - ]: 2 : if(allocated(HOcrossdisp_terms)) call polynomial_coeff_list_free(HOcrossdisp_terms)
592 : : end associate
593 : : else ! to_skip
594 : : ncombi2=0
595 : : ncombi1=0
596 [ # # # # : 0 : ABI_MALLOC(my_coeffs,(size(eff_pot%anharmonics_terms%coefficients)))
# # # # #
# # # # #
# # ]
597 : 0 : my_coeffs = eff_pot%anharmonics_terms%coefficients
598 : : endif
599 : 1 : ncombi = ncombi1_real + ncombi2_real
600 : 1 : nterm_start = eff_pot%anharmonics_terms%ncoeff
601 : : else ! if iterm = nterm + 1 => Take care about strain
602 : : block
603 : 1 : integer :: max_nbody_tmp(order_ran(2))
604 [ + - ]: 1 : if (bound_option==1) then
605 [ + + ]: 7 : max_nbody_tmp(:) = 1
606 : : else
607 [ # # ]: 0 : max_nbody_tmp(:) = 888
608 : : endif
609 : 1 : call opt_getHOstrain(my_coeffs,ncombi,nterm_start,eff_pot,order_ran,comm, max_nbody=max_nbody_tmp)
610 : : end block
611 : : endif !
612 : :
613 : :
614 : : ! Modification of Alireza.
615 : : !---------------------------------------------
616 : : !Here my_coeffs contains the previous coeff and the new bounding ones. wight is allways +1 so the power in SAT generation does not matter!!
617 : :
618 : : ! 1. We need to use generateTermsFromList module to create terms.
619 : : ! 2. The generateTermsFromList module accepts list combinations of terms: (\1,1,1,1,7,7\)
620 : : ! 3. We need to convert the terms from bounding process to a list like (\1,1,1,1,7,7\)
621 : : ! 4. use polynomial_coeff_getList to get pairs.
622 : : ! 5. compare th list and the terms and find out what is the number of diplacement in list for a term.
623 : : ! 6. convert term to list.
624 : : ! 7. ask for the temrs from generateTermsFromList
625 : : ! 8. check and see if the coeff is already considered or not. use coeffs_compare function
626 : : !! strain the same thing should be done for the strain
627 : : !! find the list of all the list_str,
628 : : !! fond what is number of strain in the displacemt comapring it to
629 : : !! the list_str and give this to
630 : : ! get coeff from > generateTermsFromList(list_disp,****)
631 : : ! setcoeff_to_tmp_coeffs
632 : : ! generateTermsFromList(cell,index_coeff,list_coeff,list_str,ncoeff,ndisp_max,nrpt,nstr,nsym,nterm,terms)
633 : :
634 : : ! call the new function to get the data required to call polynomial_coeff_getList
635 : :
636 : :
637 : 2 : call generate_bounding_term_and_add_to_list( sympairs, nterm_start, ncombi, my_coeffs, temp_cntr)
638 : :
639 : :
640 : :
641 [ + - ]: 2 : if (temp_cntr>0) then
642 [ + + ]: 5 : do icombi=1,temp_cntr
643 : : ! Copy all the terms in eff pot
644 : : call polynomial_coeff_getName(name, &
645 : 3 : & my_coeffs(nterm_start+icombi),symbols,recompute=.TRUE.)
646 : 3 : call polynomial_coeff_SetName(name,my_coeffs(nterm_start+icombi))
647 : :
648 : : ! Set dimensions of temporary my_coeffs array
649 : 3 : nterm2 = eff_pot%anharmonics_terms%ncoeff + 1
650 [ - + + - : 18 : ABI_MALLOC(my_coeffs_tmp,(nterm2))
+ - + - -
+ - + - +
+ + ]
651 : : ! Copy terms of previous cycle
652 : 3 : my_coeffs_tmp(1:nterm2-1) = eff_pot%anharmonics_terms%coefficients
653 : : !Put new term to my_coeffs_tmp
654 [ + - + - : 33 : my_coeffs_tmp(nterm2) = my_coeffs(nterm_start+icombi)
+ + + - +
- + - + -
+ - + - +
- + - - +
- - - - -
- - - - -
- - - - -
- ]
655 : :
656 : : !Do not use an intrinsic assignment here: the defined assignment (coeffs_list_copy)
657 : : !only applies to arrays, so the scalar assignment would deep-copy the allocatable
658 : : !components without ABI_MALLOC, and the memory tracking would detect a mismatch
659 : : !when they are freed with polynomial_coeff_list_free.
660 : : call polynomial_coeff_init(my_coeffs(nterm_start+icombi)%coefficient, &
661 : : & my_coeffs(nterm_start+icombi)%nterm, &
662 : : & my_coeffs_tmp(nterm2),my_coeffs(nterm_start+icombi)%terms, &
663 : 3 : & my_coeffs(nterm_start+icombi)%name)
664 : :
665 : : ! If order is greater than specified cycle
666 [ + + + + : 14 : if(sum(my_coeffs_tmp(nterm2)%terms(1)%power_disp) &
+ + - + ]
667 : : & +sum(my_coeffs_tmp(nterm2)%terms(1)%power_strain) > maxval(order_ran))then
668 : 0 : call polynomial_coeff_list_free(my_coeffs_tmp)
669 : 0 : cycle
670 : : endif
671 : : ! Message to Output
672 : 3 : write(message,'(5a)' )ch10,&
673 [ + - ]: 6 : & ' ==> high order term: ', trim(my_coeffs_tmp(nterm2)%name),' created',ch10
674 : 3 : call wrtout(ab_out,message,'COLL')
675 : 3 : call wrtout(std_out,message,'COLL')
676 : : ! Check if generated term is not already contained in effpot
677 : : ! If yes cycle
678 [ - + + - : 9 : ABI_MALLOC(exists, (nterm2))
+ - + - -
+ - + -
+ ]
679 [ + + ]: 12 : exists=.FALSE.
680 [ + + ]: 9 : do jterm=1,nterm2-1
681 : 9 : exists(jterm) = coeffs_compare(my_coeffs_tmp(jterm),my_coeffs_tmp(nterm2))
682 : : enddo !jterm
683 [ + + + - ]: 12 : any_exists=any(exists)
684 [ + - ]: 3 : ABI_SFREE(exists)
685 [ - + ]: 3 : if(any_exists)then
686 : 0 : write(message,'(3a)' )ch10,&
687 : 0 : & ' ==> Term exists already. We cycle',ch10
688 : 0 : call wrtout(ab_out,message,'COLL')
689 : 0 : call wrtout(std_out,message,'COLL')
690 : 0 : call polynomial_coeff_list_free(my_coeffs_tmp)
691 : 0 : cycle
692 : : endif
693 : :
694 : :
695 : : ! Set new term into effective potential
696 : 3 : call effective_potential_setCoeffs(my_coeffs_tmp,eff_pot,nterm2)
697 : :
698 : : ! Tell the world what we do, They want to know.
699 : 3 : write(message,'(3a)' )ch10,&
700 : 6 : & ' ==> Optimizing coefficient',ch10
701 : 3 : call wrtout(ab_out,message,'COLL')
702 : 3 : call wrtout(std_out,message,'COLL')
703 : :
704 : : ! Deallocation in loop
705 : 3 : call polynomial_coeff_list_free(my_coeffs_tmp)
706 : : !Optimizing coefficient old style
707 : :
708 [ + + ]: 5 : if(iterm>nterm)then
709 : : ! MS 2006 Decomment for old style optimization
710 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
711 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
712 : 2 : & compute_anharmonic=.TRUE.,print_file=.FALSE., weights=my_weights)
713 : 2 : i = 0
714 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses)/(msef_ini+mses_ini): ", (msef+mses)/(msef_ini+mses_ini)
715 : 2 : call wrtout(std_out,message,'COLL')
716 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses): ", (msef+mses)
717 : 2 : call wrtout(std_out,message,'COLL')
718 [ - + ]: 2 : do while((msef+mses)/(msef_ini+mses_ini) >= 1.001)
719 : 0 : i = i + 1
720 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = &!coeff_ini / 2**i
721 : 0 : & eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient / 2**1
722 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
723 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
724 : 0 : & compute_anharmonic=.TRUE.,print_file=.FALSE., weights=my_weights)
725 : :
726 : 0 : write(message,'(a,I2,a,ES24.16)') "cycle ",i," (msef+mses)/(msef_ini+mses_ini): ",(msef+mses)/(msef_ini+mses_ini)
727 : 0 : call wrtout(std_out,message,'COLL')
728 : 0 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses): ", (msef+mses)
729 : 0 : call wrtout(std_out,message,'COLL')
730 : : enddo ! while mse/mse_ini>1.0001
731 : 2 : write(message,'(a,ES24.16)') "coeff after opt:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
732 : 2 : call wrtout(std_out,message,'COLL')
733 : 2 : msef_ini = msef
734 : 2 : mses_ini = mses
735 : : ! !Optimize coefficient with opt routine
736 : : ! optterm(1)= nterm2
737 : : ! nterm_opt = 1
738 : : ! call opt_effpot(eff_pot,nterm_opt,optterm,hist,comm,print_anh=.FALSE.)
739 : : ! write(*,*) "coeff after opt:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
740 : : ! !Store new "inital precision for next coefficient
741 : : ! msef_ini = msef
742 : : ! mses_ini = mses
743 : :
744 : : else
745 : : !Optimizing coefficient with GF criterion
746 : 1 : coeff_opt = 0
747 : 1 : GF_arr = 0
748 : 1 : i = 1
749 [ + + ]: 3 : do while(i<=2)
750 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = &
751 [ + - ]: 2 : & eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient/ 2**(i-1)
752 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
753 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
754 : 2 : & compute_anharmonic=.TRUE.,print_file=.FALSE.)
755 : : GF_arr(i) = (bound_factors(1)*bound_EFS(1)*mse+bound_factors(2)*bound_EFS(2)*msef&
756 : : & +bound_factors(3)*bound_EFS(3)*mses) / &
757 : : & (bound_factors(1)*bound_EFS(1)*mse_ini+bound_factors(2)*bound_EFS(2)*msef_ini&
758 : 2 : & +bound_factors(3)*bound_EFS(3)*mses_ini)
759 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ",i," GF/GF_ini: ",GF_arr(i)
760 : 2 : call wrtout(std_out,message,'COLL')
761 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," GF: ",(bound_factors(1)*bound_EFS(1)*mse&
762 : : & +bound_factors(2)*bound_EFS(2)*msef&
763 : 4 : & +bound_factors(3)*bound_EFS(3)*mses)
764 : 2 : call wrtout(std_out,message,'COLL')
765 : 2 : coeff_opt(i) = eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
766 [ + + - + ]: 3 : if(i==2 .and. abs(GF_arr(1)-GF_arr(2)) < tol8)then
767 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient =&
768 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient*10d5
769 : 0 : write(message,'(5a)') ch10,"Differences between test-cycles to small increase",ch10, &
770 : 0 : & "test coefficient value by factor 1000",ch10
771 : 0 : call wrtout(std_out,message,'COLL')
772 : 0 : i = i+1
773 : : ! This is weird, it keeps on increasing the coefficient by 1000, until a NaN is reached.
774 : : else
775 : 2 : i=i+1
776 : : end if
777 : : enddo ! while mse/mse_ini>10
778 [ - + ]: 1 : if (eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient /= eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient) then
779 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = 0.0_dp
780 : : end if
781 : :
782 : 1 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = opt_boundcoeff(GF_arr,coeff_opt,bound_penalty)
783 [ - + ]: 1 : if (eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient /= eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient) then
784 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = 0.0_dp
785 : : end if
786 : :
787 : 1 : write(message,'(a,ES24.16)') "coeff after opt1:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
788 : 1 : call wrtout(std_out,message,'COLL')
789 : 1 : coeff_tmp = ANINT(eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient*10d10)
790 : 1 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = coeff_tmp/10d10
791 : 1 : write(message,'(a,ES24.16)') "coeff after opt2:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
792 : 1 : call wrtout(std_out,message,'COLL')
793 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
794 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
795 : 1 : & compute_anharmonic=.TRUE.,print_file=.FALSE.)
796 : 1 : write(message,'(a,ES24.16)') "GF/GF_ini after_opt: ", (bound_factors(1)*bound_EFS(1)*mse&
797 : : & +bound_factors(2)*bound_EFS(2)*msef&
798 : : & +bound_factors(3)*bound_EFS(3)*mses) / &
799 : : & (bound_factors(1)*bound_EFS(1)*mse_ini&
800 : : & +bound_factors(2)*bound_EFS(2)*msef_ini&
801 : 2 : & +bound_factors(3)*bound_EFS(3)*mses_ini)
802 : 1 : call wrtout(std_out,message,'COLL')
803 : 1 : mse_ini = mse
804 : 1 : msef_ini = msef
805 : 1 : mses_ini = mses
806 : : endif
807 : : !DEALLOCATION
808 : : !ABI_FREE(exists)
809 : : enddo ! icombi
810 : : end if
811 : :
812 : 5 : call polynomial_coeff_list_free(my_coeffs)
813 : : end do !iterm
814 : : !enddo ! order
815 : :
816 [ + - ]: 1 : if(allocated(singledisp_terms)) call polynomial_coeff_list_free(singledisp_terms)
817 : :
818 : 1 : write(message, '(a,(80a),a)' ) ch10,&
819 [ + - + + ]: 81 : &('_',ii=1,80),ch10
820 : 1 : call wrtout(ab_out,message,'COLL')
821 : 1 : call wrtout(std_out,message,'COLL')
822 : :
823 : 1 : write(message,'(3a)' )ch10,&
824 : 2 : & ' Finished creating high-order terms',ch10 !,&
825 : : !& ' Optimize initial anharmonic terms !NOT IS COMMENTED NOW!',ch10
826 : 1 : call wrtout(ab_out,message,'COLL')
827 : 1 : call wrtout(std_out,message,'COLL')
828 : : ! call opt_effpot(eff_pot,nterm,terms,hist,comm,print_anh=.FALSE.)
829 : :
830 : : ! Print the standard devition of final model
831 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
832 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
833 : 1 : & ' with respect to the training-set after attempted bounding (meV^2/atm):',&
834 : 1 : & ch10,' Energy : ',&
835 : 1 : & mse_ini* (Ha_EV*1000)**2 *factor ,ch10,&
836 : 1 : & ' Goal function values of the effective.potential',ch10,&
837 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
838 : 1 : & ' Forces+Stresses : ',&
839 : 1 : & (msef_ini+mses_ini)*(HaBohr_eVAng)**2,ch10,&
840 : 1 : & ' Forces : ',&
841 : 1 : & msef_ini*(HaBohr_eVAng)**2,ch10,&
842 : 1 : & ' Stresses : ',&
843 : 2 : & mses_ini*(HaBohr_eVAng)**2,ch10
844 : 1 : call wrtout(ab_out,message,'COLL')
845 : 1 : call wrtout(std_out,message,'COLL')
846 : :
847 : :
848 : : !DEALLOCATION
849 [ + - ]: 1 : ABI_SFREE(symbols)
850 [ + - ]: 1 : ABI_SFREE(terms)
851 [ + - ]: 1 : ABI_SFREE(my_weights)
852 : :
853 : : !ABI_FREE(my_coeffs)
854 : 1 : call fit_data_free(fit_data)
855 : 1 : call sympairs%free()
856 : :
857 : :
858 [ - + - + : 1 : end subroutine opt_effpotbound
- + - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - + - -
- - - - -
- - - - -
- - - - -
- - - - +
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + -
+ ]
859 : : !!***
860 : :
861 : :
862 : : !!****f* m_opt_effpot/opt_getHOforterm
863 : : !!
864 : : !! NAME
865 : : !! opt_effpotbound
866 : : !!
867 : : !! FUNCTION
868 : : !! Compute possible high orders for a given anharmonic term.
869 : : !! In the range of order_start,order_stop
870 : : !!
871 : : !!
872 : : !! INPUTS
873 : : !! term<polynomial_coeff_type>:anharmonic term
874 : : !! order_range(2):start and stop order desired by user
875 : : !!
876 : : !!
877 : : !! OUTPUT
878 : : !! order_start:possible start order
879 : : !! order_stop: possible stop order
880 : : !!
881 : : !! SOURCE
882 : :
883 : 2 : subroutine opt_getHOforterm(term,order_range,order_start,order_stop)
884 : :
885 : : !Arguments ------------------------------------
886 : : !scalars
887 : : type(polynomial_coeff_type),intent(in) :: term
888 : : !arrays
889 : : integer,intent(in) :: order_range(2)
890 : : integer,intent(out) :: order_start, order_stop
891 : : !Logicals
892 : : !Strings
893 : : !Local variables ------------------------------
894 : : !scalars
895 : : integer :: idisp,ndisp,nstrain,nterm_of_term,power_tot
896 : : integer :: nbody_tot
897 : : !arrays
898 : 2 : integer,allocatable :: powers(:)
899 : : !Logicals
900 : : !Strings
901 : : !*************************************************************************
902 : :
903 : : !Get/Initialize variables
904 : 2 : ndisp = term%terms(1)%ndisp
905 : 2 : nstrain = term%terms(1)%nstrain
906 : 2 : nbody_tot = ndisp + nstrain
907 : 2 : nterm_of_term = term%nterm
908 : :
909 [ - + + - : 6 : ABI_MALLOC(powers,(nbody_tot))
+ - + - -
+ - + ]
910 [ + + ]: 6 : powers(:ndisp) = term%terms(1)%power_disp
911 [ + + ]: 3 : powers(ndisp+1:) = term%terms(1)%power_strain
912 : 7 : power_tot = 0
913 : : !write(std_out,*) "powers in getHOforterm", powers
914 : :
915 : : !Get rid off odd displacements
916 [ + + ]: 7 : do idisp=1,nbody_tot
917 [ + + + + ]: 7 : if(idisp <= ndisp .and. mod(powers(idisp),2) == 1)then
918 : 2 : powers(idisp) = powers(idisp) + 1
919 [ + + ]: 3 : else if(mod(powers(idisp),2) == 1)then
920 : 1 : powers(idisp) = powers(idisp) + 1
921 : : endif
922 : : enddo !idisp
923 : : ! Count order
924 [ + + ]: 7 : do idisp=1,nbody_tot
925 : 7 : power_tot = power_tot + powers(idisp)
926 : : enddo
927 : : ! Get start and stop order for this term
928 : : ! If term doesn't fit in order range give back order_start = order_stop = 0
929 [ + + + - ]: 2 : if(power_tot >= order_range(1) .and. power_tot <=order_range(2))then
930 : 1 : order_start = power_tot
931 : 1 : order_stop = order_range(2)
932 [ + - ]: 1 : elseif(power_tot < order_range(1))then
933 : 1 : order_start = order_range(1)
934 : 1 : order_stop = order_range(2)
935 [ # # ]: 0 : elseif(power_tot > order_range(2))then
936 : 0 : order_start = 0
937 : 0 : order_stop = 0
938 : : endif
939 : :
940 [ + - ]: 2 : ABI_SFREE(powers)
941 : :
942 : 2 : end subroutine opt_getHOforterm
943 : : !!***
944 : :
945 : :
946 : : !!****f* m_opt_effpot/opt_getCombisforterm
947 : : !!
948 : : !! NAME
949 : : !! opt_getCombisforterm
950 : : !!
951 : : !! FUNCTION
952 : : !! For a given order range: order_start, order_stop
953 : : !! calculate number of total possible combinations ncombi
954 : : !! and calculat combinations per order ncombi_order(i)
955 : : !!
956 : : !!
957 : : !! INPUTS
958 : : !! order_start: start order for bounding terms
959 : : !! order_end: end order for bounding terms
960 : : !! ndisp: number of displacements a given term contains
961 : : !!
962 : : !!
963 : : !! OUTPUT
964 : : !! ncombi: total number of combinations
965 : : !! ncombi_order: array with number of combinations per order
966 : : !!
967 : : !!
968 : : !! SOURCE
969 : :
970 [ - + ]: 2 : subroutine opt_getCombisforterm(order_start,order_end,ndisp,ncombi,ncombi_order)
971 : :
972 : : !Arguments ------------------------------------
973 : : !scalars
974 : : integer,intent(in) :: ndisp
975 : : integer,intent(in) :: order_start, order_end
976 : : !arrays
977 : : integer,intent(out) :: ncombi
978 : : integer,intent(out) :: ncombi_order(:)
979 : : !Logicals
980 : : !Strings
981 : : !Local variables ------------------------------
982 : : !scalars
983 : : integer :: i
984 : : integer :: order,iorder1,iorder2
985 : : !arrays
986 : : !integer
987 : : !Logicals
988 : : !Strings
989 : : character(len=1000) :: message
990 : : !*************************************************************************
991 : :
992 : : !Test
993 [ + - - + ]: 2 : if(mod(order_start,2) /= 0 .or. mod(order_end,2) /= 0)then
994 : : ! Message to Output
995 : 0 : write(message,'(4a)' )ch10,&
996 : 0 : & 'Either start or stop order are not even numbers',ch10,&
997 : 0 : & 'Action: change bound_range in input',ch10
998 : 0 : ABI_ERROR(message)
999 : : endif
1000 : :
1001 : : !Initialize Variables
1002 : 2 : i = 0
1003 : 2 : ncombi = 0
1004 [ + + ]: 4 : ncombi_order = 0
1005 : :
1006 : : !Calculate Combinations
1007 [ + - + - ]: 4 : do order=order_start,order_end,2
1008 : 2 : i = i+1
1009 [ - + ]: 4 : if(ndisp == 1)then
1010 : 0 : ncombi = ncombi + 1
1011 : 0 : ncombi_order(i) = 1
1012 : : else
1013 [ + - + + ]: 5 : do iorder1 = 2,order-2*(ndisp-1),2
1014 [ + + ]: 3 : if(ndisp*iorder1 == order)then
1015 : 1 : ncombi = ncombi + 1
1016 : 1 : ncombi_order(i) = ncombi_order(i) + 1
1017 : 1 : cycle
1018 : : endif
1019 [ + + + - ]: 5 : do iorder2=iorder1+2,order-2*(ndisp-1),2
1020 [ + - ]: 3 : if( iorder1 + (ndisp-1)*iorder2 == order)then
1021 : 1 : ncombi = ncombi + ndisp
1022 : 1 : ncombi_order(i) = ncombi_order(i) + ndisp
1023 [ # # ]: 0 : elseif(iorder1 * (ndisp-1) + iorder2 == order)then
1024 : 0 : ncombi = ncombi + ndisp
1025 : 0 : ncombi_order(i) = ncombi_order(i) + ndisp
1026 : : endif
1027 : : enddo !iorder2
1028 : : enddo !iorder1 !
1029 : : endif
1030 : : !write(*,*) 'ncombi(',i,') for order',order,' is:', ncombi_order(i), 'are we happy?'
1031 : : enddo !order
1032 : :
1033 : : !write(*,*) ncombi_order(:)
1034 : : !write(*,*) 'ncombi for term is:', ncombi, 'are we happy?'
1035 : :
1036 : :
1037 : 2 : end subroutine opt_getCombisforterm
1038 : : !!***
1039 : :
1040 : : !!****f* m_opt_effpot/opt_getHoTerms
1041 : : !!
1042 : : !! NAME
1043 : : !! opt_geHoTerms
1044 : : !!
1045 : : !! FUNCTION
1046 : : !! For a term give all possible high order terms
1047 : : !! Attention order_start, order_stop, ncombi and
1048 : : !! ncombi_order have to be calculated before!
1049 : : !!
1050 : : !!
1051 : : !!
1052 : : !! INPUTS
1053 : : !! eff_pot: existing effective potential
1054 : : !! order: order for which bounding terms are generated
1055 : : !!
1056 : : !!
1057 : : !! OUTPUT
1058 : : !! eff_pot new effective potential
1059 : : !!
1060 : : !!
1061 : : !! SOURCE
1062 : :
1063 [ - + - + ]: 2 : subroutine opt_getHoTerms(terms,order_start,order_stop,ndisp,ncombi_order)
1064 : :
1065 : : !Arguments ------------------------------------
1066 : : !scalars
1067 : : integer,intent(in) :: ndisp
1068 : : integer,intent(in) :: order_start, order_stop
1069 : : !arrays
1070 : : integer,intent(in) :: ncombi_order(:)
1071 : : !Logicals
1072 : : type(polynomial_coeff_type),intent(inout) :: terms(:)
1073 : : !Strings
1074 : : !Local variables ------------------------------
1075 : : !scalars
1076 : : integer :: i,icombi,icombi2,icombi_start,icombi_stop,idisp,nterm_of_term
1077 : : integer :: order,iterm_of_term,jdisp,power_tot
1078 : : integer :: jdisp1,jdisp2,sec,nstrain,nbody_tot
1079 : : real(sp) :: to_divide,divider1,divider2,divided
1080 : : !arrays
1081 : : !integer
1082 : : !Logicals
1083 : : logical :: equal_term_done
1084 : : !Strings
1085 : : character(len=1000) :: message
1086 : : !*************************************************************************
1087 : : !Get Variables
1088 : 2 : nterm_of_term = terms(1)%nterm
1089 : 2 : nstrain = terms(1)%terms(1)%nstrain
1090 : 2 : nbody_tot = ndisp + nstrain
1091 : :
1092 : : ! Create all possible combinaions for the specified orders
1093 : : ! If anybody has ever, ever to read and understand this I'm terribly sorry
1094 : : ! ---this will be quite hell---
1095 : 2 : power_tot = 0
1096 : 2 : icombi_start = 1
1097 : 2 : i = 0
1098 : : !write(*,*) "what is ncombi?", ncombi
1099 : : !write(*,*) "what is ncmobi_order", ncombi_order
1100 : :
1101 : : !write(*,*) "order_start", order_start
1102 : : !write(*,*) "order_stop", order_stop
1103 : 2 : order = order_start
1104 : : ! TODO work here icombi start and order counting does not work yet.
1105 [ - + - + ]: 4 : do order=order_start,order_stop,2
1106 : : !write(std_out,*) 'Was I now here?!(order-loop)'
1107 : 2 : i = i + 1
1108 : 2 : icombi_stop = icombi_start + ncombi_order(i) - 1
1109 : 2 : equal_term_done = .FALSE.
1110 : : !write(std_out,*) 'order', order
1111 : : !write(std_out,*) 'icombi_start', icombi_start
1112 : : !write(std_out,*) 'icombi_stop', icombi_stop
1113 : 2 : icombi=icombi_start
1114 [ + + ]: 6 : do while (icombi<=icombi_stop)
1115 : : power_tot = 0
1116 [ + + ]: 7 : do idisp=1,nbody_tot
1117 : : !write(*,*) "what is icombi here?", icombi
1118 : : !write(*,*) "what is icombi_stop here?", icombi_stop
1119 [ + + ]: 7 : if(idisp<=ndisp)then
1120 : 4 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(idisp)
1121 : : else
1122 : 1 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(idisp-ndisp)
1123 : : endif
1124 : : enddo
1125 : : ! Probably have to increase order at same time
1126 : : ! If the term already has the right order we cycle
1127 : : ! we increase icombi_start and icombi to go to the next term in the array
1128 [ + + ]: 4 : if(power_tot == order)then
1129 : 1 : icombi_start = icombi_start + 1
1130 : 1 : icombi = icombi + 1
1131 : : !write(*,*) 'icombi-start in if', icombi_start
1132 : 1 : cycle
1133 : : ! If the term is not already in the right order, we manipulate it until
1134 : : ! it is
1135 : : else
1136 : : jdisp1 = 1
1137 : 3 : sec = 1
1138 : : !write(*,*) 'what is ndisp actually', ndisp
1139 : : ! Treat single permutations from the bottom of the order
1140 : : ! so start from ^2^2^2 and get ^6^2^2,^2^6^2, and ^2^2^6 f.E.
1141 : : ! do loop over displacements
1142 [ + + ]: 3 : do while(jdisp1<=nbody_tot .and. sec < 100)
1143 : : !write(*,*) "what is jdisp1 here?", jdisp1
1144 : : !write(*,*) "what is icombi here?", icombi
1145 : : !write(*,*) "what is icombi_sotp?", icombi_stop
1146 : : ! Increase the order of the displacements
1147 [ + + ]: 50 : do iterm_of_term=1,nterm_of_term
1148 [ + - ]: 50 : if(jdisp1 <= ndisp)then
1149 : : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) =&
1150 : 48 : & terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) + 2
1151 : : else
1152 : : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) =&
1153 : 0 : & terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) + 2
1154 : :
1155 : : endif
1156 : : enddo
1157 : : ! Check total order of term after increase
1158 : : power_tot=0
1159 [ + + ]: 6 : do jdisp2=1,nbody_tot
1160 [ + - ]: 6 : if(jdisp2 <= ndisp)then
1161 : 4 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(jdisp2)
1162 : : else
1163 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(jdisp2-ndisp)
1164 : : endif
1165 : : enddo
1166 : : ! If the term is at the right order do next displacement
1167 : : ! increase icombi and icombi_start to go to next term in array
1168 [ - + ]: 3 : if(power_tot == order)then
1169 : 2 : icombi = icombi + 1
1170 : 2 : icombi_start = icombi_start +1
1171 : : !write(*,*) 'what is icombi_start here', icombi_start
1172 : 2 : jdisp1 = jdisp1 + 1
1173 : : !write(*,*) 'and what is jdisp1?', jdisp1
1174 : :
1175 : : endif
1176 : : enddo!jdisp1
1177 : : !Treat permutations in terms ndisp > 2 and order >= 10
1178 : : !Start f.E. from ^4^4^4 and get ^4^2^4, ^2^4^4,^4^4^2
1179 [ - + - - ]: 1 : if(icombi_stop - icombi > 1 .and. nbody_tot >2 )then
1180 [ # # ]: 0 : do icombi2=icombi,icombi+nbody_tot-1
1181 [ # # ]: 0 : do jdisp=1,nbody_tot
1182 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1183 [ # # ]: 0 : if(jdisp <= ndisp)then
1184 : : terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) = &
1185 : 0 : & terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) +2
1186 : : !write(*,*) "What's the power now?", terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp)
1187 : : else
1188 : : terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = &
1189 : 0 : & terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) +2
1190 : : endif
1191 : : enddo
1192 : : enddo !jdisp
1193 : : enddo
1194 : : jdisp1 = 1
1195 [ # # ]: 0 : do while(jdisp1<=nbody_tot .and. sec < 100)
1196 : 0 : sec = sec + 1
1197 : : !write(*,*) "how often did I go here, hu ?"
1198 : : !write(*,*) "I did at least one displacement"
1199 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1200 [ # # ]: 0 : if(jdisp1 <= ndisp)then
1201 : : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) =&
1202 : 0 : & terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) + 2
1203 : : else
1204 : : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) =&
1205 : 0 : & terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) + 2
1206 : :
1207 : : endif
1208 : : enddo
1209 : : power_tot=0
1210 [ # # ]: 0 : do jdisp2=1,nbody_tot
1211 [ # # ]: 0 : if(jdisp2 <= ndisp)then
1212 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(jdisp2)
1213 : : else
1214 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(jdisp2-ndisp)
1215 : : endif
1216 : : enddo
1217 [ # # ]: 0 : if(power_tot == order)then
1218 : 0 : icombi = icombi + 1
1219 : 0 : icombi_start = icombi_start +1
1220 : : !write(*,*) 'what is icombi_start here', icombi_start
1221 : 0 : jdisp1 = jdisp1 + 1
1222 : : !write(*,*) 'and what is jdisp1?', jdisp1
1223 : : endif
1224 : : enddo!jdisp1
1225 : : ! Message to Output
1226 : : if(sec>100)then
1227 : : write(message,'(4a)' )ch10,&
1228 : : & "You're stuck in a while loop.",ch10,&
1229 : : & 'Action: Contact Abinit Group',ch10
1230 : : ABI_ERROR(message)
1231 : : endif
1232 : : endif! (icombi_stop - icombi)
1233 : : !write(*,*) 'I was here!'
1234 : 1 : to_divide = real(order)
1235 : 1 : divider1 = real(nbody_tot)
1236 : 1 : divided = real(to_divide/divider1)
1237 : 1 : divider2 = real(2)
1238 : : !write(*,*) 'divided', divided, 'divider2', divider2
1239 : : !Treat terms with even power f.E. ^2^2^2^2, ^4^4 etc...
1240 [ + - - - ]: 1 : if(mod(divided,divider2) == 0 .and. .not. equal_term_done .and. nbody_tot > 1)then
1241 : : !write(*,*) "Sometimes I should be here sometimes I shouldn't"
1242 [ # # ]: 0 : do jdisp=1,nbody_tot
1243 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1244 [ # # ]: 0 : if(jdisp<=ndisp)then
1245 : 0 : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp) = order/ndisp
1246 : : else
1247 : 0 : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = order/ndisp
1248 : : endif
1249 : : enddo
1250 : : enddo !jdisp
1251 [ # # ]: 0 : if(order < order_stop)then
1252 [ # # ]: 0 : do icombi2=icombi+1,icombi+nbody_tot
1253 [ # # ]: 0 : do jdisp=1,nbody_tot
1254 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1255 [ # # ]: 0 : if(jdisp<=ndisp)then
1256 : 0 : terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) = order/ndisp
1257 : : else
1258 : 0 : terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = order/ndisp
1259 : : endif
1260 : : enddo
1261 : : enddo !jdisp
1262 : : enddo
1263 : : endif
1264 : 0 : equal_term_done = .TRUE.
1265 : 0 : icombi = icombi + 1
1266 : 0 : icombi_start = icombi_start +1
1267 : : endif ! equal term if
1268 : : endif ! power_tot == order
1269 : : enddo !icombination
1270 : : enddo !order
1271 : :
1272 : 2 : end subroutine opt_getHoTerms
1273 : : !!***
1274 : :
1275 : : !!****f* m_opt_effpot/opt_filterdisp
1276 : : !!
1277 : : !! NAME
1278 : : !! opt_opt_filterdisp
1279 : : !!
1280 : : !! FUNCTION
1281 : : !! If a anharmonic term represents a strain-phonon coupling
1282 : : !! delete the strain and only keep the displacement part.
1283 : : !!
1284 : : !!
1285 : : !! INPUTS
1286 : : !! term<polynomial_coeff_type>: anharmonic term to check
1287 : : !! nterm_of_term: number of symmetry equivalent terms for term
1288 : : !!
1289 : : !! OUTPUT
1290 : : !! term<polynomial_coeff_type>: only the displacement part of original term
1291 : : !!
1292 : : !! SOURCE
1293 : :
1294 : 1 : subroutine opt_filterdisp(term,nterm_of_term)
1295 : :
1296 : : !Arguments ------------------------------------
1297 : : !scalars
1298 : : type(polynomial_coeff_type),intent(inout) :: term
1299 : : integer,intent(in) :: nterm_of_term
1300 : : !arrays
1301 : : !Logicals
1302 : : !Strings
1303 : : !Local variables ------------------------------
1304 : : !scalars
1305 : : integer :: iterm_of_term
1306 : : !reals
1307 : : real(dp) :: coeff
1308 : : !arrays
1309 [ + + + + : 242 : type(polynomial_term_type) :: terms(nterm_of_term)
+ + + - +
- + - + -
+ - + - +
- ]
1310 : : !Logicals
1311 : : !Strings
1312 : : !*************************************************************************
1313 : :
1314 [ + + ]: 25 : do iterm_of_term=1,nterm_of_term
1315 : : !terms(iterm_of_term) = term%terms(iterm_of_term)
1316 : :
1317 : : call polynomial_term_init(term%terms(iterm_of_term)%atindx, &
1318 : : & term%terms(iterm_of_term)%cell,&
1319 : : & term%terms(iterm_of_term)%direction,&
1320 : : & term%terms(iterm_of_term)%ndisp, &
1321 : : & term%terms(iterm_of_term)%nstrain, &
1322 : : & terms(iterm_of_term),&
1323 : : & term%terms(iterm_of_term)%power_disp, &
1324 : : & term%terms(iterm_of_term)%power_strain,&
1325 : : & term%terms(iterm_of_term)%strain, &
1326 : : & term%terms(iterm_of_term)%weight, &
1327 : 24 : & check=.TRUE.)
1328 : 24 : terms(iterm_of_term)%nstrain = 0
1329 [ + + ]: 48 : terms(iterm_of_term)%power_strain = 0
1330 [ + + ]: 49 : terms(iterm_of_term)%strain = 0
1331 : : !DMSG(terms(iterm_of_term)%debug_str)
1332 : : enddo
1333 : :
1334 : 1 : call polynomial_coeff_free(term)
1335 : : !Reinitial term
1336 : : !check=.TRUE. checks for duplicate terms
1337 : 1 : call polynomial_coeff_init(coeff,nterm_of_term,term,terms, check=.TRUE.)
1338 : : !DMSG(term%debug_str)
1339 : :
1340 : :
1341 [ + + ]: 25 : do iterm_of_term=1,nterm_of_term
1342 : 25 : call polynomial_term_free(terms(iterm_of_term))
1343 : : enddo
1344 : : !if(nterm_of_term /= term%nterm)then
1345 : : ! write(*,*) "nterm_of_term changed after deleting strain"
1346 : : !endif
1347 : :
1348 : :
1349 : 192 : end subroutine opt_filterdisp
1350 : : !!***
1351 : :
1352 : : !!****f* m_opt_effpot/opt_getHOstrain
1353 : : !!
1354 : : !! NAME
1355 : : !! opt_getHOstrain
1356 : : !!
1357 : : !! FUNCTION
1358 : : !! Get HO anharmnonic strain terms for bounding and add them to list
1359 : : !! of existing terms in effective potential
1360 : : !!
1361 : : !! INPUTS
1362 : : !! eff_pot<effective_potential_type>: datatype with all the information
1363 : : !! about the effective potential
1364 : : !! power_strain(2): start and stop order for strain terms
1365 : : !! comm: mpi communicator (at the moment only sequential tested)
1366 : : !!
1367 : : !! OUTPUT
1368 : : !! terms<polynomial_coeff_type>: list with original terms in effective
1369 : : !! potential + HO even strain terms
1370 : : !!
1371 : : !! SOURCE
1372 : :
1373 [ - + ]: 1 : subroutine opt_getHOstrain(terms,ncombi,nterm_start,eff_pot,power_strain,comm, max_nbody)
1374 : :
1375 : : !Arguments ------------------------------------
1376 : : !scalars
1377 : : integer,intent(in) :: comm
1378 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms(:)
1379 : : type(effective_potential_type), intent(in) :: eff_pot
1380 : : integer,intent(in) :: power_strain(2), max_nbody(:)
1381 : : integer,intent(out) :: ncombi,nterm_start
1382 : : !arrays
1383 : : !Logicals
1384 : : !Strings
1385 : : !Local variables ------------------------------
1386 : : !scalars
1387 : : integer :: nterm_tot_tmp
1388 : : integer :: i,ii
1389 : : !integer :: nbody, ref_nbody
1390 : : real(dp) :: coeff_ini
1391 : : !reals
1392 [ + + + + : 56 : type(crystal_t) :: crystal
+ + + + +
+ + + + +
+ + + + ]
1393 : : !arrays
1394 : 1 : type(polynomial_coeff_type),allocatable :: strain_terms_tmp(:)
1395 : : !Logicals
1396 : : !Strings
1397 : : character(len=1000) :: message
1398 : : !*************************************************************************
1399 : : !Get variables
1400 [ + - + - : 1 : crystal = eff_pot%crystal
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - ]
1401 : 1 : coeff_ini = 1000000
1402 : :
1403 : 1 : write(message, '(a,(80a),a)' ) ch10,&
1404 [ + - + + ]: 81 : & ('_',ii=1,80),ch10
1405 : 1 : call wrtout(ab_out,message,'COLL')
1406 : 1 : call wrtout(std_out,message,'COLL')
1407 : 1 : write(message,'(3a)' )ch10,&
1408 : 2 : & ' Chreate high order strain terms ',ch10
1409 : 1 : call wrtout(ab_out,message,'COLL')
1410 : 1 : call wrtout(std_out,message,'COLL')
1411 : :
1412 : :
1413 : : !1406 get count of high order even anharmonic strain terms and the strain terms itself
1414 : 1 : call polynomial_coeff_getEvenAnhaStrain(strain_terms_tmp,crystal,ncombi,power_strain,comm, max_nbody)
1415 : : ! Allocate my_coeffs with ncombi free space to work with
1416 : :
1417 : :
1418 : :
1419 : 1 : nterm_start = eff_pot%anharmonics_terms%ncoeff
1420 : 1 : nterm_tot_tmp = eff_pot%anharmonics_terms%ncoeff + ncombi
1421 [ - + + - : 7 : ABI_MALLOC(terms,(nterm_tot_tmp))
+ - + - -
+ - + - +
+ + + - ]
1422 [ + + ]: 5 : do i=1,nterm_tot_tmp
1423 [ + + ]: 5 : if(i<=nterm_start)then
1424 : : call polynomial_coeff_init(coeff_ini,eff_pot%anharmonics_terms%coefficients(i)%nterm,terms(i),&
1425 : : & eff_pot%anharmonics_terms%coefficients(i)%terms,eff_pot%anharmonics_terms%coefficients(i)%name,&
1426 : 2 : & check=.TRUE.)
1427 : : !DMSG(terms(i)%debug_str)
1428 : : else
1429 [ - + ]: 2 : if (.not.allocated(strain_terms_tmp(i-nterm_start)%terms)) then
1430 : 0 : write(message,'(4a)' )ch10,&
1431 : 0 : & 'Error in opt_getHOstrain: strain_terms_tmp(i-nterm_start)%terms not allocated',ch10,&
1432 : 0 : & 'Action: Contact Abinit Group',ch10
1433 : 0 : ABI_ERROR(message)
1434 : : endif
1435 : : call polynomial_coeff_init(coeff_ini,strain_terms_tmp(i-nterm_start)%nterm,terms(i),&
1436 : : & strain_terms_tmp(i-nterm_start)%terms,strain_terms_tmp(i-nterm_start)%name,&
1437 : 2 : & check=.TRUE.)
1438 : : !DMSG(strain_terms_tmp(i-nterm_start)%debug_str)
1439 : : endif
1440 : : enddo
1441 : :
1442 : 1 : call polynomial_coeff_list_free(strain_terms_tmp)
1443 : :
1444 [ - + - - : 1 : end subroutine opt_getHOstrain
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - ]
1445 : : !!***
1446 : :
1447 : : !!****f* m_opt_effpot/opt_getHOcrossdisp
1448 : : !!
1449 : : !! NAME
1450 : : !! opt_getHOcrossdisp
1451 : : !!
1452 : : !! FUNCTION
1453 : : !! Get even high order displacement terms for a given input term and
1454 : : !! add them to an excisting list of terms. If the term is strain phonon type
1455 : : !! the strain part gets deleted and the high order terms for the displacement
1456 : : !! part are computed.
1457 : : !! Example: for a tree linear term with three displacements x*y*z the even high order
1458 : : !! possibilites are computed and stored.
1459 : : !! For range 6 to 8: x^2*y^2*z^2,x^4*y^2*z^2,x^2*y^4*z^2,x^2*y^2*z^4.
1460 : : !!
1461 : : !! INPUTS
1462 : : !! eff_pot<effective_potential_type>: datatype with all the information
1463 : : !! about the effective potential
1464 : : !! power_disp(2): start and stop order for disp terms
1465 : : !! comm: mpi communicator (at the moment only sequential tested)
1466 : : !!
1467 : : !! OUTPUT
1468 : : !! terms<polynomial_coeff_type>: list with original terms in effective
1469 : : !! potential + HO even disp terms
1470 : : !!
1471 : : !! SOURCE
1472 : :
1473 : 1 : subroutine opt_getHOcrossdisp(terms_out,ncombi,term_in,power_disp)
1474 : :
1475 : : !Arguments ------------------------------------
1476 : : !scalars
1477 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms_out(:)
1478 : : type(polynomial_coeff_type),intent(inout) :: term_in
1479 : : integer,intent(in) :: power_disp(2)
1480 : : integer,intent(out) :: ncombi
1481 : : !arrays
1482 : : !Logicals
1483 : : !Strings
1484 : : !Local variables ------------------------------
1485 : : !scalars
1486 : : integer :: ndisp,nterm_of_term,nstrain,nbody_tot
1487 : : integer :: order_start,order_stop,norder
1488 : : integer :: order_start_str,order_stop_str
1489 : : integer :: icombi,idisp,iterm_of_term
1490 : : integer :: ncombi_tot,ncombi_str
1491 : : !reals
1492 : : real(dp) :: coeff_ini=1
1493 : : !arrays
1494 : 1 : type(polynomial_coeff_type) :: term
1495 : : integer,allocatable :: ncombi_order(:),ncombi_order_str(:)
1496 : : !Logicals
1497 : : logical :: had_strain
1498 : : !Strings
1499 : : character(len=1000) :: message
1500 : : !*************************************************************************
1501 : : !Get/Set Variables
1502 : 1 : norder = abs(((power_disp(2)-power_disp(1))/2)) + 1
1503 [ - + + - : 3 : ABI_MALLOC(ncombi_order,(norder))
+ - - + -
+ ]
1504 [ - + + - : 2 : ABI_MALLOC(ncombi_order_str,(norder))
- + - + ]
1505 [ + + ]: 2 : ncombi_order = 0
1506 [ + + ]: 2 : ncombi_order_str = 0
1507 : 1 : order_start_str = 0
1508 : 1 : order_stop_str = 0
1509 : :
1510 : 1 : ncombi = 0
1511 : : !Get this term (iterm) and infromations about it
1512 : : !Get number of displacements and equivalent terms for this term
1513 : : !Chose term one to get ndisp. ndisp is equal for all terms of the term
1514 : : !Get minimum oder for this term
1515 : : !Get total number of terms in effpot for message
1516 : 1 : ndisp = term_in%terms(1)%ndisp
1517 : 1 : nstrain = term_in%terms(1)%nstrain
1518 : 1 : nbody_tot = ndisp + nstrain
1519 : 1 : nterm_of_term = term_in%nterm
1520 : 1 : call polynomial_coeff_init(coeff_ini,nterm_of_term,term,term_in%terms,term_in%name, check=.true.)
1521 : : !DMSG(term%debug_str)
1522 : : ! Check if term has strain component.
1523 : : ! If yes filter strain and fit high order atomic displacement terms
1524 : 1 : had_strain = .FALSE.
1525 [ + - ]: 1 : if(term%terms(1)%nstrain /= 0)then
1526 : : ! Message to Output
1527 : 1 : write(message,'(5a)' )ch10,&
1528 : 1 : & '- Term has strain compenent',ch10,&
1529 : 2 : & ' -> Filter Displacement',ch10
1530 : 1 : call wrtout(ab_out,message,'COLL')
1531 : 1 : call wrtout(std_out,message,'COLL')
1532 : 1 : call opt_filterdisp(term,nterm_of_term)
1533 : : !Get new value of symmetry equivalent term nterm_of_term
1534 : 1 : nterm_of_term = term%nterm
1535 : : !Remember if this term had strain
1536 : 1 : had_strain = .TRUE.
1537 : : !cycle
1538 : : endif
1539 : : ! Ok we want it. Let's go.
1540 : :
1541 : : ! get start and stop order for this term
1542 : 1 : call opt_getHOforterm(term,power_disp,order_start,order_stop)
1543 [ + - ]: 1 : if(had_strain) call opt_getHOforterm(term_in,power_disp,order_start_str,order_stop_str)
1544 [ - + ]: 1 : if(order_start == 0)then
1545 : : ! Message to Output
1546 : 0 : write(message,'(5a,I2,a,I2,3a)' )ch10,&
1547 : 0 : & " ==> High order cross product terms for term ", trim(term%name),ch10,&
1548 : 0 : & " ==> do not fit into specified order range from ", power_disp(1),' to ',power_disp(2),ch10,&
1549 : 0 : & " ==> Can not construct high order cross product bounding term",ch10
1550 : 0 : call wrtout(ab_out,message,'COLL')
1551 : 0 : call wrtout(std_out,message,'COLL')
1552 [ # # ]: 0 : ABI_SFREE(ncombi_order)
1553 [ # # ]: 0 : ABI_SFREE(ncombi_order_str)
1554 : 0 : return
1555 : : end if
1556 : :
1557 [ - + - - ]: 1 : if(order_start_str == 0 .and. had_strain)then
1558 : : ! Message to Output
1559 : 0 : write(message,'(5a,I2,a,I2,3a)' )ch10,&
1560 : 0 : & " ==> High order cross product terms for term ", trim(term_in%name),ch10,&
1561 : 0 : & " ==> do not fit into specified order range from ", power_disp(1),' to ',power_disp(2),ch10,&
1562 : 0 : & " ==> Can not construct high order cross product bounding term",ch10
1563 : 0 : call wrtout(ab_out,message,'COLL')
1564 : 0 : call wrtout(std_out,message,'COLL')
1565 : 0 : had_strain = .FALSE.
1566 : : end if
1567 : : ! get total amount of combinations and combinations per order for the term
1568 : 1 : call opt_getCombisforterm(order_start,order_stop,ndisp,ncombi,ncombi_order)
1569 [ + - ]: 1 : if(had_strain)then
1570 : 1 : call opt_getCombisforterm(order_start_str,order_stop_str,nbody_tot,ncombi_str,ncombi_order_str)
1571 : : !write(std_out,*) "I was here ncombi_str is: ", ncombi_str
1572 [ - + + - : 6 : ABI_MALLOC(terms_out,(ncombi+ncombi_str))
+ - + - -
+ - + - +
+ + + - ]
1573 : 2 : ncombi_tot = ncombi+ncombi_str
1574 : : else
1575 [ # # # # : 0 : ABI_MALLOC(terms_out,(ncombi))
# # # # #
# # # # #
# # # # ]
1576 : 0 : ncombi_tot = ncombi
1577 : : endif
1578 : : ! Copy current term to the ncombination elemenst a the end in array terms
1579 : : ! change the value of their coefficient to a start value
1580 : : ! The start is estimed to not be larger then half the initial term's value
1581 : : ! This is because higher order terms should have smaller coefficients
1582 [ + + ]: 4 : do icombi=1,ncombi_tot
1583 [ + + ]: 3 : if(icombi <= ncombi)then
1584 : 2 : coeff_ini = 1 !abs(terms_out(icombi)%coefficient / 2)
1585 : 2 : nterm_of_term = term%nterm
1586 : 2 : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out(icombi),term%terms, check=.true.)
1587 : : !DMSG(terms_out(icombi)%debug_str)
1588 : : else
1589 : 1 : coeff_ini = 10d3
1590 : 1 : nterm_of_term = term_in%nterm
1591 : 1 : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out(icombi),term_in%terms, check=.true.)
1592 : : !DMSG(terms_out(icombi)%debug_str)
1593 : : endif
1594 : : ! Set the power of all terms we want to add to two. We find the correct power later
1595 : : ! Change the weight of the term to 1 (even terms have allways weight=1)
1596 [ + + ]: 76 : do iterm_of_term=1,nterm_of_term
1597 : 72 : terms_out(icombi)%terms(iterm_of_term)%weight = 1
1598 [ + + ]: 72 : if(icombi > ncombi)then
1599 [ + + ]: 48 : terms_out(icombi)%terms(iterm_of_term)%power_strain = 2
1600 : 24 : coeff_ini = 10d3 !abs(terms_out(icombi)%coefficient / 2)
1601 : 24 : terms_out(icombi)%coefficient = coeff_ini
1602 : : ! elseif(icombi>2*ncombi)then
1603 : : ! terms_out(icombi)%terms(iterm_of_term)%power_strain = 4
1604 : : ! coeff_ini = 10d5 !abs(terms_out(icombi)%coefficient / 2)
1605 : : ! terms_out(icombi)%coefficient = coeff_ini
1606 : : endif
1607 [ + + ]: 219 : do idisp=1,ndisp
1608 : 216 : terms_out(icombi)%terms(iterm_of_term)%power_disp(idisp) = 2
1609 : : enddo !idisp
1610 : : enddo !iterm_of_term
1611 : : enddo !icombi
1612 : :
1613 : : !If term had strain we had to reinitialize it in the process
1614 : : !Refree memory
1615 : : !if(had_strain)
1616 : 1 : call polynomial_coeff_free(term)
1617 : :
1618 : : ! Get high order combinations
1619 [ + - ]: 1 : if(had_strain)then
1620 : 1 : call opt_getHoTerms(terms_out(:ncombi),order_start,order_stop,ndisp,ncombi_order)
1621 : 1 : call opt_getHoTerms(terms_out(ncombi+1:),order_start_str,order_stop_str,ndisp,ncombi_order_str)
1622 : 1 : ncombi = ncombi_tot
1623 : : else
1624 : 0 : call opt_getHoTerms(terms_out,order_start,order_stop,ndisp,ncombi_order)
1625 : : endif
1626 : : !DEALLOCATION
1627 [ + - ]: 1 : ABI_SFREE(ncombi_order)
1628 [ + - ]: 1 : ABI_SFREE(ncombi_order_str)
1629 : :
1630 [ - + - - : 3 : end subroutine opt_getHOcrossdisp
- - - - -
- - - - -
- - - - -
+ - + ]
1631 : : !!***
1632 : :
1633 : : !!****f* m_opt_effpot/opt_getSingleDispTerms
1634 : : !!
1635 : : !! NAME
1636 : : !! opt_getSingleDispTerms
1637 : : !!
1638 : : !! FUNCTION
1639 : : !! Get polynomial terms (<polynomial_coeff_type>) with single displacements
1640 : : !! at second order inside a given range defined by the supercell size
1641 : : !!
1642 : : !! INPUTS
1643 : : !! sc_size(3): supercell size
1644 : : !! crystal<type(crystal_t)>: all information about the crystal
1645 : : !! comm: mpi communicator (at the moment only sequential tested)
1646 : : !!
1647 : : !! OUTPUT
1648 : : !! terms<polynomial_coeff_type>: list single displacement polynomial_coeffs
1649 : : !!
1650 : : !! SOURCE
1651 : :
1652 : 1 : subroutine opt_getSingleDispTerms(terms,crystal, sc_size,comm)
1653 : :
1654 : : !Arguments ------------------------------------
1655 : : !scalars
1656 : : integer,intent(in) :: comm
1657 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms(:)
1658 : : type(crystal_t),intent(inout) :: crystal
1659 : : real(dp) :: cutoff
1660 : : !arrays
1661 : : integer :: sc_size(3)
1662 : :
1663 : : !Logicals
1664 : : !Strings
1665 : : !Local variables ------------------------------
1666 : : !scalars
1667 : : integer :: natom,nsym,nrpt,ncoeff_sym,nstr_sym
1668 : : integer :: ncoeff,ncoeff_out,power_strph,option_GN,option
1669 : : integer :: nterms_out,nterm1,iterm1,iterm2,ind,iatom,i
1670 : : integer :: ncopy
1671 : : integer :: ii !,ia,ib,r1,r2,r3
1672 : : !integer :: irpt,irpt_ref
1673 : : integer :: master,nproc,my_rank
1674 : : !arrays
1675 : : integer :: power_disp(2)
1676 : 1 : integer,allocatable :: cell(:,:)
1677 : 1 : integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:)
1678 : 1 : logical,allocatable :: terms_to_copy(:)
1679 : 1 : type(polynomial_coeff_type),allocatable :: terms_tmp(:),terms_tmp2(:)
1680 : : !type(polynomial_coeff_type),allocatable :: terms(:)
1681 : : !real(dp),allocatable :: xcart(:,:),xred(:,:),rpt(:,:)
1682 : : real(dp) :: rprimd(3,3),range_ifc(3)
1683 : 1 : real(dp),allocatable :: dist(:,:,:,:)
1684 : 1 : character(len=5),allocatable :: symbols(:)
1685 : : !Logicals
1686 : : logical :: iam_master,need_verbose
1687 : : !Strings
1688 : : character(len=1000) :: message
1689 : : !*************************************************************************
1690 : :
1691 : 1 : option = 2
1692 : :
1693 : : if(option == 1)then
1694 : : !MPI variables
1695 : : master = 0
1696 : : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1697 : : iam_master = (my_rank == master)
1698 : :
1699 : : !Set/Get other variables
1700 : : !natom = crystal%natom
1701 : : !nsym = crystal%nsym
1702 : : rprimd = crystal%rprimd
1703 : : need_verbose = .TRUE.
1704 : :
1705 : :
1706 : : call prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc , symbols)
1707 : :
1708 : : ! FIXME: this is dangerous!
1709 : : ! It assumes that the structure is cubic in a usual axis setting.
1710 : : cutoff = rprimd(1,1)
1711 : :
1712 : : if(need_verbose)then
1713 : : write(message,'(1a)')' Generation of the list of all the possible coefficients'
1714 : : call wrtout(std_out,message,'COLL')
1715 : : end if
1716 : : call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1717 : : & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size)
1718 : :
1719 : : ABI_SFREE(dist)
1720 : : !ABI_FREE(rpt)
1721 : :
1722 : : !write(*,*) "polynomial_getList worked"
1723 : :
1724 : : !Get Order1 term s
1725 : : !Check difference between ncoeff_out, ncoeff
1726 : : call polynomial_coeff_getOrder1(cell,terms,list_symcoeff,natom,nterms_out,ncoeff_sym,nrpt,nsym,symbols)
1727 : :
1728 : : !do ii=1,nterms_out
1729 : : ! write(*,*) "Term(",ii,"/",nterms_out,"): ", terms(ii)%name
1730 : : !enddo
1731 : :
1732 : : elseif(option == 2)then
1733 : : !Get/set variables
1734 : : ! FIXME: hexu: check these values of 2. Why are they are hard coded here.
1735 : : ! Does this mean we only generate high-order terms in 2x2x2 cell?
1736 : 1 : power_disp = (/2,2/)
1737 : 1 : power_strph = zero
1738 : 1 : option_GN = 0
1739 : 1 : sc_size = (/2,2,2/)
1740 : 1 : cutoff = 0
1741 : 1 : natom = crystal%natom
1742 : : ! XXX: hexu: why cutoff = sum of axis lengths?
1743 [ + + ]: 4 : do ii=1,3
1744 : : cutoff = cutoff + sqrt(crystal%rprimd(ii,1)**2 + &
1745 : : & crystal%rprimd(ii,2)**2 + &
1746 : 4 : & crystal%rprimd(ii,3)**2)
1747 : : enddo
1748 : :
1749 : :
1750 : : !do iatom=1,3
1751 [ + + ]: 6 : do iatom=1,natom !Subhadeep
1752 : : ! TODO: to be validated.
1753 : : !do iatom=1,crystal%nirredat
1754 : : call polynomial_coeff_getNorder(terms_tmp,crystal,cutoff,ncoeff,ncoeff_out,power_disp,&
1755 : : & power_strph,option_GN,sc_size,comm,anharmstr=.false.,spcoupling=.false.,&
1756 : : & only_odd_power=.false.,only_even_power=.true.,verbose=.false.,&
1757 : : & compute_symmetric=.false.,fit_iatom=iatom, &
1758 : 5 : & max_nbody=[999,999,999,999,999,999, 999,999,999,999])
1759 : : !TEST MS
1760 : : ! write(std_out,*) "behind call getNorder"
1761 : : ! write(std_out,*) "ncoeff_out: ", ncoeff_out
1762 : : ! do ii=1,ncoeff_out
1763 : : ! write(*,*) "Term(",ii,"/",ncoeff_out,"): ", terms_tmp(ii)%name
1764 : : ! enddo
1765 : : !TEST MS
1766 [ + + ]: 6 : if(iatom == 1)then
1767 [ - + + - : 13 : ABI_MALLOC(terms,(size(terms_tmp)))
+ - + - -
+ - + - +
+ + + - ]
1768 : 1 : call coeffs_list_copy(terms,terms_tmp)
1769 : : else
1770 [ + + + - : 12 : ABI_MALLOC(terms_to_copy,(size(terms_tmp)))
+ - + + -
+ - + -
+ ]
1771 : : ! note: when iatom>1, terms is already initialized.
1772 : 4 : nterm1 = size(terms)
1773 [ - + + - : 145 : ABI_MALLOC(terms_tmp2,(nterm1))
+ - + - -
+ - + - +
+ + ]
1774 : 4 : terms_tmp2 = terms
1775 [ + + ]: 44 : terms_to_copy = .TRUE.
1776 [ + + ]: 44 : do iterm1=1,size(terms_tmp)
1777 [ + + ]: 795 : do iterm2=1,size(terms)
1778 [ - + ]: 791 : if(terms_tmp(iterm1) == terms(iterm2))then
1779 : 0 : terms_to_copy(iterm1) = .FALSE.
1780 : 0 : exit
1781 : : endif
1782 : : enddo!iterm1
1783 : : enddo!iterm2
1784 [ + + + - ]: 44 : ncopy = count(terms_to_copy)
1785 : 4 : call polynomial_coeff_list_free(terms)
1786 [ - + + - : 185 : ABI_MALLOC(terms,(nterm1+ncopy))
+ - + - -
+ - + - +
+ + + - ]
1787 : 4 : call coeffs_list_copy(terms(:nterm1),terms_tmp2)
1788 : 4 : ind = 0
1789 [ + + ]: 44 : do i =1,size(terms_tmp)
1790 [ + - ]: 44 : if(terms_to_copy(i))then
1791 : 40 : ind=ind+1
1792 : : call polynomial_coeff_init(terms_tmp(i)%coefficient,terms_tmp(i)%nterm,terms(nterm1+ind),terms_tmp(i)%terms,&
1793 : 40 : & terms_tmp(i)%name,check=.TRUE.)
1794 : : !DMSG(terms(nterm1+ind)%debug_str)
1795 : : endif
1796 : : enddo!i=1,nterm1+ncopy
1797 : 4 : call polynomial_coeff_list_free(terms_tmp)
1798 : 4 : call polynomial_coeff_list_free(terms_tmp2)
1799 [ - + ]: 4 : ABI_SFREE(terms_to_copy)
1800 : : endif!iatom==1
1801 : : enddo !iatom=1,natom
1802 : : ! call polynomial_coeff_getNorder(terms,crystal,cutoff,ncoeff,ncoeff_out,power_disp,&
1803 : : !& power_strph,option_GN,sc_size,comm,anharmstr=.false.,spcoupling=.false.,&
1804 : : !& only_odd_power=.false.,only_even_power=.true.,verbose=.false.,&
1805 : : !& compute_symmetric=.false.)
1806 : :
1807 : : endif !option
1808 : : !TEST MS
1809 : : ! write(std_out,*) "behind call getNorder"
1810 : : ! write(std_out,*) "ncoeff_out: ", size(terms)
1811 : : ! do ii=1,size(terms)
1812 : : ! write(*,*) "Term(",ii,"/",size(terms),"): ", terms(ii)%name
1813 : : ! enddo
1814 : : !TEST MS
1815 : :
1816 [ - + - + : 1 : end subroutine opt_getSingleDispTerms
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - ]
1817 : : !!***
1818 : :
1819 : : !!****f* m_opt_effpot/opt_getHOSingleDispTerms
1820 : : !!
1821 : : !! NAME
1822 : : !! opt_getHOSingleDispTerms
1823 : : !!
1824 : : !! FUNCTION
1825 : : !! For a given anharmonic term that might consits of product of
1826 : : !! terms find all even high order single displacement terms
1827 : : !! Example: Input term: x*y*z
1828 : : !! Input HO-range: 6 8
1829 : : !! Output terms: x^6,x^8,y^6,y^8,z^6,z^8
1830 : : !! Caution only atomic displacements are taken into account
1831 : : !!
1832 : : !! INPUTS
1833 : : !! term_in<polynomial_coeff_type>: input anharmonic term
1834 : : !! crystal<type(crystal_t)>: all information about the crystal
1835 : : !! single_disp_terms<polynomial_coeff_out>: list of single disp terms at
1836 : : !! second order to select terms from
1837 : : !! power_disp(2): Start and stop power for HO terms
1838 : : !! comm: mpi communicator (at the moment only sequential tested)
1839 : : !!
1840 : : !! OUTPUT
1841 : : !! terms_out<polynomial_coeff_out>: output high order even terms
1842 : : !! ncoeff: number of coefficients
1843 : : !!
1844 : : !! SOURCE
1845 : :
1846 [ - + - + ]: 1 : subroutine opt_getHOSingleDispTerms(term_in,terms_out,symbols,single_disp_terms,power_disp,ncoeff)
1847 : :
1848 : : !Arguments ------------------------------------
1849 : : !scalars
1850 : : integer,intent(out) :: ncoeff
1851 : : type(polynomial_coeff_type),intent(in) :: term_in
1852 : : type(polynomial_coeff_type),intent(in) :: single_disp_terms(:)
1853 : : type(polynomial_coeff_type),allocatable,intent(out) :: terms_out(:)
1854 : : !type(crystal_t),intent(inout) :: crystal
1855 : : !arrays
1856 : : integer, intent(in) :: power_disp(2)
1857 : : character(len=5),intent(in) :: symbols(:)
1858 : : !Logicals
1859 : : !Strings
1860 : : !Local variables ------------------------------
1861 : : !scalars
1862 : : integer :: ndisp,norder, nterm_of_term
1863 : : integer :: icoeff,iorder,idisp, iterm1,iterm2,iterm3
1864 : : real(dp) :: coeff_ini = 1
1865 : : !Strings
1866 : : character(len=200):: name
1867 : : !arrays
1868 : : type(polynomial_coeff_type),allocatable :: terms_out_tmp(:)
1869 : : !Logicals
1870 : 1 : logical,allocatable :: found(:)
1871 : : !*************************************************************************
1872 : : !Get/Set Variables
1873 : : !Number of output terms
1874 : 1 : ndisp = term_in%terms(1)%ndisp
1875 : 1 : norder = abs(power_disp(2)-power_disp(1))/2 + 1
1876 : 1 : ncoeff = norder * ndisp
1877 : : !Allocate output terms
1878 [ - + + - : 5 : ABI_MALLOC(terms_out_tmp,(ncoeff))
+ - + - -
+ - + +
+ ]
1879 : : !find equivalent second order terms in list of single disp terms
1880 : : !for each displacement in input term
1881 : : !Transfer to output term and increase order
1882 : 1 : icoeff = 0
1883 [ + + ]: 3 : do idisp=1,ndisp
1884 [ + + ]: 103 : do iterm1=1,size(single_disp_terms)
1885 [ + + ]: 2346 : do iterm2=1,single_disp_terms(iterm1)%nterm
1886 [ + + + + : 2980 : if(all(term_in%terms(1)%atindx(:,idisp) == single_disp_terms(iterm1)%terms(iterm2)%atindx(:,1)))then
+ + ]
1887 [ + + + + : 606 : if(all(term_in%terms(1)%cell(:,:,idisp) == single_disp_terms(iterm1)%terms(iterm2)%cell(:,:,1)))then
+ + + + ]
1888 [ + + ]: 6 : if(term_in%terms(1)%direction(idisp) == single_disp_terms(iterm1)%terms(iterm2)%direction(1))then
1889 [ + + ]: 4 : do iorder=1,norder
1890 : 2 : icoeff = icoeff + 1
1891 [ + - + - : 32 : terms_out_tmp(icoeff) = single_disp_terms(iterm1)
+ + + - +
- + - + -
+ - + - +
- + - - +
- - - - -
- - - - -
- - - - -
- ]
1892 : 2 : nterm_of_term = single_disp_terms(iterm1)%nterm
1893 : : !Change order of term
1894 : : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out_tmp(icoeff),&
1895 : 2 : & single_disp_terms(iterm1)%terms(:), check=.true.)
1896 : : !DMSG(terms_out_tmp(icoeff)%debug_str)
1897 : :
1898 [ + + ]: 34 : do iterm3=1,nterm_of_term
1899 [ + + ]: 62 : terms_out_tmp(icoeff)%terms(iterm3)%power_disp = power_disp(1) + (iorder-1)*2
1900 : : enddo !iterm3
1901 : : enddo !iorder
1902 : : endif
1903 : : endif
1904 : : endif
1905 : : enddo !iterm2
1906 : : enddo !iterm1
1907 : : enddo!idisp
1908 : : !Change Name
1909 [ + + ]: 3 : do icoeff=1,ncoeff
1910 : 2 : call polynomial_coeff_getName(name,terms_out_tmp(icoeff),symbols,recompute=.TRUE.)
1911 : 3 : call polynomial_coeff_SetName(name,terms_out_tmp(icoeff))
1912 : : enddo
1913 : :
1914 : : !Check for doubles and delete them
1915 : : !First count irreducible terms
1916 [ - + + - : 3 : ABI_MALLOC(found,(ncoeff))
+ - + - -
+ - + ]
1917 [ + + ]: 3 : found = .FALSE.
1918 : : iterm3 = 0
1919 [ + + ]: 3 : do iterm1=1,ncoeff
1920 [ + + ]: 4 : do iterm2=iterm1+1,ncoeff
1921 [ - + - - ]: 3 : if(terms_out_tmp(iterm1) == terms_out_tmp(iterm2) .and. .not. found(iterm2))then
1922 : 0 : found(iterm2) = .TRUE.
1923 : 0 : iterm3 = iterm3 + 1
1924 : : endif
1925 : : enddo
1926 : : enddo
1927 : 1 : iterm3 = ncoeff - iterm3
1928 [ - + + - : 5 : ABI_MALLOC(terms_out,(iterm3))
+ - + - -
+ - + - +
+ + + - ]
1929 : 1 : iterm3 = 0
1930 [ + + ]: 3 : do iterm1=1,ncoeff
1931 [ + - ]: 3 : if(.not. found(iterm1))then
1932 : 2 : iterm3 = iterm3 + 1
1933 : : call polynomial_coeff_init(coeff_ini,terms_out_tmp(iterm1)%nterm,terms_out(iterm3),&
1934 : 2 : & terms_out_tmp(iterm1)%terms,terms_out_tmp(iterm1)%name, check=.TRUE.)
1935 : : !DMSG(terms_out(iterm3)%debug_str)
1936 : : endif
1937 : : enddo
1938 [ + - ]: 1 : ABI_SFREE(found)
1939 : 1 : call polynomial_coeff_list_free(terms_out_tmp)
1940 : 1 : ncoeff = iterm3
1941 [ - + - - : 1 : end subroutine opt_getHOSingleDispTerms
- - - - -
- - - - -
- - - - -
- - - ]
1942 : : !!***
1943 : :
1944 : :
1945 : :
1946 : : !!****f* m_opt_effpot/opt_boundcoeff
1947 : : !! NAME
1948 : : !! opt_boundcoeff
1949 : : !!
1950 : : !! FUNCTION
1951 : : !!
1952 : : !! optimize a bound coefficient if optimized value is negative
1953 : : !! put an positive value that respects a precision penalty
1954 : : !!
1955 : : !! INPUTS
1956 : : !!
1957 : : !! OUTPUT
1958 : : !!
1959 : : !! SOURCE
1960 : :
1961 : 1 : function opt_boundcoeff(yvalues,cvalues,penalty_in) result (coeff)
1962 : : !Arguments ------------------------------------
1963 : :
1964 : : !Arguments ------------------------------------
1965 : : real(dp),intent(in) :: yvalues(2),cvalues(2),penalty_in
1966 : : real(dp) :: coeff
1967 : : !local
1968 : : !variable
1969 : : real(dp) :: a,b,coeff_tmp,x1,x2,penalty
1970 : : !array
1971 : : ! *************************************************************************
1972 : :
1973 : 1 : a = ( (yvalues(1) - 1) - (yvalues(2)-1)*(cvalues(1)/cvalues(2))) / (cvalues(1)**2 - cvalues(1)*cvalues(2))
1974 : :
1975 : : b = ( (yvalues(2) - 1)/cvalues(2) ) - ( (yvalues(1) -1)*cvalues(2) - (yvalues(2) - 1)*cvalues(1) )&
1976 : 1 : & / (cvalues(1)**2 - cvalues(1)*cvalues(2))
1977 : :
1978 : 1 : penalty = penalty_in - 1
1979 [ + - + - ]: 1 : if (a/=a .or. b/=b .or. a==0) then
1980 : 1 : coeff = 0.0_dp
1981 : : return
1982 : : end if
1983 : 1 : coeff_tmp = -b/(2*a)
1984 : : !write(*,*) "coeff_tmp", coeff_tmp
1985 [ + - ]: 1 : if(coeff_tmp > 0)then
1986 : : coeff = coeff_tmp
1987 : : else
1988 : 1 : x1 = (-b + sqrt(b**2 + 4*a*penalty)) / (2*a) ! 1.001 penalty value
1989 : 1 : x2 = (-b - sqrt(b**2 + 4*a*penalty)) / (2*a)
1990 [ - + ]: 1 : if(x1>0)then
1991 : : coeff = x1
1992 : : else
1993 : 0 : coeff = x2
1994 : : endif
1995 : : endif
1996 : :
1997 : : end function opt_boundcoeff
1998 : : !!***
1999 : :
2000 : : !!****f* m_opt_effpot/check_to_skip
2001 : : !! NAME
2002 : : !! check_to_skip
2003 : : !!
2004 : : !! FUNCTION
2005 : : !!
2006 : : !! Check if term contains only bodies with even power
2007 : : !! and has a positive coefficient. If yes term doesn't need
2008 : : !! a bounding high order equivalent and we can skip it.
2009 : : !! Function retursn logical to_skip.
2010 : : !!
2011 : : !! INPUTS
2012 : : !!
2013 : : !! term<polynomial_coeff_type>:anharmonic term
2014 : : !!
2015 : : !! OUTPUT
2016 : : !!
2017 : : !! logical: to_skip
2018 : : !!
2019 : : !! SOURCE
2020 : :
2021 : 1 : function check_to_skip(term) result (to_skip)
2022 : : !Arguments ------------------------------------
2023 : : type(polynomial_coeff_type),intent(in) :: term
2024 : : logical :: to_skip
2025 : : ! ------------------------------------
2026 : : !local
2027 : : !variable
2028 : : character(len=1000) :: message
2029 : : !array
2030 : : ! *************************************************************************
2031 : :
2032 : 1 : to_skip = .FALSE.
2033 : : ! Let's check if we really want all this mess
2034 : : ! If the term is even and its coefficient positive we skip it. Also here we take terms(1) as example for all equivalent terms of term
2035 [ + - + + : 2 : if(term%coefficient > 0 .and. .not. any(mod(term%terms(1)%power_disp(:),2) /= 0))then
- + - - ]
2036 [ # # # # : 0 : if(.not. any(mod(term%terms(1)%power_strain(:),2) /= 0))then
# # ]
2037 : : ! Message to Output
2038 : 0 : write(message,'(3a)' )ch10,&
2039 : 0 : & ' ==> No need for high order bounding term',ch10
2040 : 0 : call wrtout(ab_out,message,'COLL')
2041 : 0 : call wrtout(std_out,message,'COLL')
2042 : 0 : to_skip = .TRUE.
2043 : 0 : return
2044 : : end if
2045 : : end if
2046 : :
2047 : : end function check_to_skip
2048 : : !!***
2049 : :
2050 : 2 : subroutine generate_bounding_term_and_add_to_list(sympairs, nterm_start, ncombi, my_coeffs, temp_cntr)
2051 : : ! check if myterm
2052 : : !type(polynomial_coeff_type),target, intent(in) :: terms(:)
2053 : :
2054 : : type(polynomial_coeff_type), allocatable, intent(inout) :: my_coeffs(:)
2055 : : type(SymPairs_t), intent(inout) :: sympairs
2056 : : integer, intent(in) :: nterm_start, ncombi
2057 : : integer, intent(out) :: temp_cntr
2058 : :
2059 : 2 : type(polynomial_term_type), allocatable :: terms2(:)
2060 : : type(polynomial_term_type), pointer :: myterm=> null()
2061 : : type(polynomial_coeff_type), target, allocatable :: my_coeffs_tmp(:)
2062 : : integer :: idisp, my_nrpt, list_cntr, number_coeff, icombi
2063 : : integer :: jterm
2064 : : integer :: tmp_list(5), pwr
2065 : : logical :: found
2066 : : logical :: reverse_i, any_exists
2067 : 2 : logical, allocatable :: reverse(:), false_reverse(:), exists(:)
2068 : : integer :: counter, nterm, tot_power
2069 : 2 : type(polynomial_coeff_type):: temp_coeff
2070 : 2 : integer, allocatable :: list_disp(:)
2071 : : character(len=200):: name
2072 : : integer :: ncoeff
2073 : :
2074 : : ! copy the coefficients to a temporary array my_coeffs_tmp
2075 : 2 : ncoeff=size(my_coeffs)
2076 : :
2077 : :
2078 [ - + + - : 12 : ABI_MALLOC(my_coeffs_tmp,(ncoeff))
+ - + - -
+ - + +
+ ]
2079 : :
2080 : : ! check if my_coeffs i allocated and has the right size
2081 [ - + ]: 2 : if(.not. allocated(my_coeffs)) then
2082 : 0 : ABI_BUG("my_coeffs is not allocated")
2083 : : end if
2084 [ - + ]: 2 : if(size(my_coeffs) /= ncoeff) then
2085 : 0 : ABI_BUG("my_coeffs has the wrong size")
2086 : : end if
2087 [ - + ]: 2 : if(size(my_coeffs_tmp) /= ncoeff) then
2088 : 0 : ABI_BUG("my_coeffs_tmp has the wrong size")
2089 : : end if
2090 : :
2091 : 2 : my_coeffs_tmp=my_coeffs
2092 : :
2093 : : !call coeffs_list_copy(my_coeffs_tmp, my_coeffs)
2094 : :
2095 : : ! free the original array and allocate a new one with the correct size
2096 : 2 : call polynomial_coeff_list_free(my_coeffs)
2097 [ - + + - : 10 : ABI_MALLOC(my_coeffs, (ncoeff))
- + - + -
+ + + +
- ]
2098 : : !my_coeffs(1:nterm_start) = my_coeffs_tmp(1:nterm_start)
2099 : 2 : call coeffs_list_copy(my_coeffs(1:nterm_start), my_coeffs_tmp(1:nterm_start))
2100 : :
2101 : :
2102 : :
2103 : 2 : temp_cntr =0
2104 [ + + ]: 5 : do icombi=1, ncombi
2105 [ + - ]: 3 : if(size( my_coeffs_tmp(nterm_start+icombi)%terms) == 0) cycle
2106 : 3 : myterm=> my_coeffs_tmp(nterm_start+icombi)%terms(1)
2107 [ + + + + ]: 8 : tot_power=sum(myterm%power_disp) + sum(myterm%power_strain)
2108 [ - + + - : 9 : ABI_MALLOC(list_disp, (tot_power))
+ - + - -
+ - + -
+ ]
2109 [ - + + - : 6 : ABI_MALLOC(reverse, (tot_power))
- + - + -
+ ]
2110 [ - + + - : 6 : ABI_MALLOC(false_reverse, (tot_power))
- + - + -
+ ]
2111 : :
2112 [ + + ]: 21 : list_disp(:)=0
2113 [ + + ]: 21 : reverse(:)=.False.
2114 [ + + ]: 21 : false_reverse(:)=.False.
2115 : 3 : counter=0
2116 : : ! put the indices of the disp terms into a list
2117 [ + + ]: 5 : do idisp=1,myterm%ndisp
2118 : : ! for each disp, find its index in the pair list.
2119 : 2 : my_nrpt=find_irpt(cells=sympairs%cell, cell=myterm%cell(:, 2, idisp))
2120 : : tmp_list(:) = (/myterm%direction(idisp),myterm%atindx(1,idisp), &
2121 [ + + ]: 12 : & myterm%atindx(2,idisp),my_nrpt,int(myterm%weight) /)
2122 : 2 : found=.False. ! to find if the symmetry adapted term already exist.
2123 : 2 : reverse_i=.False.
2124 [ + + ]: 1518 : do list_cntr=1,size(sympairs%list_symcoeff, 2)
2125 : : !do isym=1, sympairs%nsym
2126 [ + + + + : 2220 : if (all(tmp_list(:4)==sympairs%list_symcoeff(:4,list_cntr,1))) then
+ + ]
2127 : 2 : number_coeff = sympairs%list_symcoeff(6,list_cntr,1)
2128 : 2 : found = .True.
2129 : : !if (found) exit
2130 : : end if
2131 : : !end do
2132 : : !if(found) exit
2133 : : end do
2134 [ - + ]: 2 : if(.not. found) then
2135 : 0 : reverse_i=.True.
2136 [ # # # # ]: 0 : my_nrpt=find_irpt(cells=sympairs%cell, cell=-myterm%cell(:, 2, idisp))
2137 : : tmp_list(:) = (/myterm%direction(idisp),myterm%atindx(2,idisp), &
2138 [ # # ]: 0 : & myterm%atindx(1,idisp),my_nrpt,-int(myterm%weight) /)
2139 [ # # ]: 0 : do list_cntr=1,size(sympairs%list_symcoeff, 2)
2140 : : !do isym=1, sympairs%nsym
2141 [ # # # # : 0 : if (all(tmp_list(:4)==sympairs%list_symcoeff(:4,list_cntr,1))) then
# # ]
2142 : 0 : number_coeff = sympairs%list_symcoeff(6,list_cntr,1)
2143 : 0 : found = .True.
2144 : : !if (found) exit
2145 : : end if
2146 : : !end do
2147 : : !if(found) exit
2148 : : end do
2149 : : end if
2150 : :
2151 [ - - ]: 3 : if (found )then
2152 [ + + ]: 6 : do pwr=1,myterm%power_disp(idisp)
2153 : 4 : counter=counter+1
2154 : 4 : list_disp(counter) = number_coeff
2155 : 6 : reverse(counter) = reverse_i
2156 : : end do
2157 : : else
2158 : 0 : ABI_BUG("The pair is not found in generated bounding terms.")
2159 : : end if
2160 : : end do
2161 : :
2162 [ + + ]: 6 : do idisp=1,myterm%nstrain
2163 [ + + ]: 20 : do pwr=1,myterm%power_strain(idisp)
2164 : 14 : counter = counter+1
2165 : 17 : list_disp(counter) = myterm%strain(idisp)+size(sympairs%list_symcoeff,2)
2166 : : end do
2167 : : end do
2168 : 3 : nullify(myterm)
2169 : :
2170 : : ! generate the terms with false reverse, and save to the temp_coeff.
2171 : 3 : nterm=sympairs%nsym
2172 [ - + + - : 153 : ABI_MALLOC(terms2, (nterm))
+ - + - -
+ - + - +
+ + ]
2173 : 3 : call sympairs%generateTerms(list_disp, tot_power, nterm, terms2, reverse=false_reverse)
2174 : : call polynomial_coeff_init(one,nterm,temp_coeff, &
2175 : 3 : &terms2(1:nterm), check=.true.)
2176 : 3 : call polynomial_term_list_free(terms2)
2177 : : !DMSG(temp_coeff%debug_str)
2178 : :
2179 : : call polynomial_coeff_getName(name, &
2180 : 3 : & temp_coeff,sympairs%symbols,recompute=.TRUE.)
2181 : :
2182 : : ! check if temp coeff already exists in the list of coeffs my_coeffs
2183 [ + - - + : 9 : ABI_MALLOC(exists, (nterm_start+icombi-1))
+ - + - +
- - + -
+ ]
2184 [ + + ]: 9 : exists=.False.
2185 [ + + ]: 9 : do jterm=1,(nterm_start+icombi-1)
2186 : 9 : exists(jterm) = coeffs_compare(my_coeffs(jterm),temp_coeff)
2187 : : enddo !jterm
2188 : 3 : call polynomial_coeff_free(temp_coeff)
2189 : :
2190 [ + + + - ]: 9 : any_exists=any(exists)
2191 [ - + ]: 3 : ABI_FREE(exists)
2192 : :
2193 : : ! if the term already exists, we skip it. Otherwise we add also the reverse.
2194 [ + - ]: 3 : if (.not. any_exists) then
2195 : 3 : temp_cntr = temp_cntr+1
2196 [ - + + - : 153 : ABI_MALLOC(terms2, (nterm))
+ - + - -
+ - + - +
+ + ]
2197 : 3 : call sympairs%generateTerms(list_disp, tot_power, nterm, terms2, reverse=reverse)
2198 : : call polynomial_coeff_init(one,nterm, &
2199 : 3 : & my_coeffs(nterm_start+temp_cntr),terms2, check=.true.)
2200 : : !DMSG(my_coeffs(nterm_start+temp_cntr)%debug_str)
2201 : 3 : call polynomial_term_list_free(terms2)
2202 : : end if
2203 : :
2204 [ + - ]: 3 : ABI_SFREE(list_disp)
2205 [ + - ]: 3 : ABI_SFREE(reverse)
2206 [ - + ]: 5 : ABI_SFREE(false_reverse)
2207 : : end do
2208 : :
2209 : 2 : call polynomial_coeff_list_free(my_coeffs_tmp)
2210 : :
2211 : : ! truncate my_coeffs.
2212 : 2 : call coeffs_list_truncate(my_coeffs, nterm_start+temp_cntr)
2213 [ - + - - : 2 : end subroutine generate_bounding_term_and_add_to_list
- - - - -
- - - - -
- - - - -
+ - - - -
- - - - -
- - - - -
- - - + -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + ]
2214 : :
2215 : :
2216 : : end module m_opt_effpot
2217 : : !!***
|
