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 : :
657 : : !call polynomial_coeff_init(my_coeffs(nterm_start+icombi)%coefficient, &
658 : : ! & my_coeffs(nterm_start+icombi)%nterm, &
659 : : ! & my_coeffs_tmp(nterm2),my_coeffs(nterm_start+icombi)%terms, &
660 : : ! &my_coeffs(nterm_start+icombi)%name)
661 : :
662 : : ! DMSG(my_coeffs_tmp(nterm2)%debug_str)
663 : :
664 : : ! If order is greater than specified cycle
665 [ + + + + : 14 : if(sum(my_coeffs_tmp(nterm2)%terms(1)%power_disp) &
+ + - + ]
666 : : & +sum(my_coeffs_tmp(nterm2)%terms(1)%power_strain) > maxval(order_ran))then
667 : 0 : call polynomial_coeff_list_free(my_coeffs_tmp)
668 : 0 : cycle
669 : : endif
670 : : ! Message to Output
671 : 3 : write(message,'(5a)' )ch10,&
672 [ + - ]: 6 : & ' ==> high order term: ', trim(my_coeffs_tmp(nterm2)%name),' created',ch10
673 : 3 : call wrtout(ab_out,message,'COLL')
674 : 3 : call wrtout(std_out,message,'COLL')
675 : : ! Check if generated term is not already contained in effpot
676 : : ! If yes cycle
677 [ - + + - : 9 : ABI_MALLOC(exists, (nterm2))
+ - + - -
+ - + -
+ ]
678 [ + + ]: 12 : exists=.FALSE.
679 [ + + ]: 9 : do jterm=1,nterm2-1
680 : 9 : exists(jterm) = coeffs_compare(my_coeffs_tmp(jterm),my_coeffs_tmp(nterm2))
681 : : enddo !jterm
682 [ + + + - ]: 12 : any_exists=any(exists)
683 [ + - ]: 3 : ABI_SFREE(exists)
684 [ - + ]: 3 : if(any_exists)then
685 : 0 : write(message,'(3a)' )ch10,&
686 : 0 : & ' ==> Term exists already. We cycle',ch10
687 : 0 : call wrtout(ab_out,message,'COLL')
688 : 0 : call wrtout(std_out,message,'COLL')
689 : 0 : call polynomial_coeff_list_free(my_coeffs_tmp)
690 : 0 : cycle
691 : : endif
692 : :
693 : :
694 : : ! Set new term into effective potential
695 : 3 : call effective_potential_setCoeffs(my_coeffs_tmp,eff_pot,nterm2)
696 : :
697 : : ! Tell the world what we do, They want to know.
698 : 3 : write(message,'(3a)' )ch10,&
699 : 6 : & ' ==> Optimizing coefficient',ch10
700 : 3 : call wrtout(ab_out,message,'COLL')
701 : 3 : call wrtout(std_out,message,'COLL')
702 : :
703 : : ! Deallocation in loop
704 : 3 : call polynomial_coeff_list_free(my_coeffs_tmp)
705 : : !Optimizing coefficient old style
706 : :
707 [ + + ]: 5 : if(iterm>nterm)then
708 : : ! MS 2006 Decomment for old style optimization
709 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
710 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
711 : 2 : & compute_anharmonic=.TRUE.,print_file=.FALSE., weights=my_weights)
712 : 2 : i = 0
713 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses)/(msef_ini+mses_ini): ", (msef+mses)/(msef_ini+mses_ini)
714 : 2 : call wrtout(std_out,message,'COLL')
715 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses): ", (msef+mses)
716 : 2 : call wrtout(std_out,message,'COLL')
717 [ - + ]: 2 : do while((msef+mses)/(msef_ini+mses_ini) >= 1.001)
718 : 0 : i = i + 1
719 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = &!coeff_ini / 2**i
720 : 0 : & eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient / 2**1
721 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
722 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
723 : 0 : & compute_anharmonic=.TRUE.,print_file=.FALSE., weights=my_weights)
724 : :
725 : 0 : write(message,'(a,I2,a,ES24.16)') "cycle ",i," (msef+mses)/(msef_ini+mses_ini): ",(msef+mses)/(msef_ini+mses_ini)
726 : 0 : call wrtout(std_out,message,'COLL')
727 : 0 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," (msef+mses): ", (msef+mses)
728 : 0 : call wrtout(std_out,message,'COLL')
729 : : enddo ! while mse/mse_ini>1.0001
730 : 2 : write(message,'(a,ES24.16)') "coeff after opt:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
731 : 2 : call wrtout(std_out,message,'COLL')
732 : 2 : msef_ini = msef
733 : 2 : mses_ini = mses
734 : : ! !Optimize coefficient with opt routine
735 : : ! optterm(1)= nterm2
736 : : ! nterm_opt = 1
737 : : ! call opt_effpot(eff_pot,nterm_opt,optterm,hist,comm,print_anh=.FALSE.)
738 : : ! write(*,*) "coeff after opt:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
739 : : ! !Store new "inital precision for next coefficient
740 : : ! msef_ini = msef
741 : : ! mses_ini = mses
742 : :
743 : : else
744 : : !Optimizing coefficient with GF criterion
745 : 1 : coeff_opt = 0
746 : 1 : GF_arr = 0
747 : 1 : i = 1
748 [ + + ]: 3 : do while(i<=2)
749 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = &
750 [ + - ]: 2 : & eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient/ 2**(i-1)
751 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
752 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
753 : 2 : & compute_anharmonic=.TRUE.,print_file=.FALSE.)
754 : : GF_arr(i) = (bound_factors(1)*bound_EFS(1)*mse+bound_factors(2)*bound_EFS(2)*msef&
755 : : & +bound_factors(3)*bound_EFS(3)*mses) / &
756 : : & (bound_factors(1)*bound_EFS(1)*mse_ini+bound_factors(2)*bound_EFS(2)*msef_ini&
757 : 2 : & +bound_factors(3)*bound_EFS(3)*mses_ini)
758 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ",i," GF/GF_ini: ",GF_arr(i)
759 : 2 : call wrtout(std_out,message,'COLL')
760 : 2 : write(message,'(a,I2,a,ES24.16)') "cycle ", i ," GF: ",(bound_factors(1)*bound_EFS(1)*mse&
761 : : & +bound_factors(2)*bound_EFS(2)*msef&
762 : 4 : & +bound_factors(3)*bound_EFS(3)*mses)
763 : 2 : call wrtout(std_out,message,'COLL')
764 : 2 : coeff_opt(i) = eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
765 [ + + - + ]: 3 : if(i==2 .and. abs(GF_arr(1)-GF_arr(2)) < tol8)then
766 : : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient =&
767 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient*10d5
768 : 0 : write(message,'(5a)') ch10,"Differences between test-cycles to small increase",ch10, &
769 : 0 : & "test coefficient value by factor 1000",ch10
770 : 0 : call wrtout(std_out,message,'COLL')
771 : 0 : i = i+1
772 : : ! This is weird, it keeps on increasing the coefficient by 1000, until a NaN is reached.
773 : : else
774 : 2 : i=i+1
775 : : end if
776 : : enddo ! while mse/mse_ini>10
777 [ - + ]: 1 : if (eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient /= eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient) then
778 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = 0.0_dp
779 : : end if
780 : :
781 : 1 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = opt_boundcoeff(GF_arr,coeff_opt,bound_penalty)
782 [ - + ]: 1 : if (eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient /= eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient) then
783 : 0 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = 0.0_dp
784 : : end if
785 : :
786 : 1 : write(message,'(a,ES24.16)') "coeff after opt1:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
787 : 1 : call wrtout(std_out,message,'COLL')
788 : 1 : coeff_tmp = ANINT(eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient*10d10)
789 : 1 : eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient = coeff_tmp/10d10
790 : 1 : write(message,'(a,ES24.16)') "coeff after opt2:", eff_pot%anharmonics_terms%coefficients(nterm2)%coefficient
791 : 1 : call wrtout(std_out,message,'COLL')
792 : : call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,&
793 : : & natom_sc,ntime,fit_data%training_set%sqomega,comm,&
794 : 1 : & compute_anharmonic=.TRUE.,print_file=.FALSE.)
795 : 1 : write(message,'(a,ES24.16)') "GF/GF_ini after_opt: ", (bound_factors(1)*bound_EFS(1)*mse&
796 : : & +bound_factors(2)*bound_EFS(2)*msef&
797 : : & +bound_factors(3)*bound_EFS(3)*mses) / &
798 : : & (bound_factors(1)*bound_EFS(1)*mse_ini&
799 : : & +bound_factors(2)*bound_EFS(2)*msef_ini&
800 : 2 : & +bound_factors(3)*bound_EFS(3)*mses_ini)
801 : 1 : call wrtout(std_out,message,'COLL')
802 : 1 : mse_ini = mse
803 : 1 : msef_ini = msef
804 : 1 : mses_ini = mses
805 : : endif
806 : : !DEALLOCATION
807 : : !ABI_FREE(exists)
808 : : enddo ! icombi
809 : : end if
810 : :
811 : 5 : call polynomial_coeff_list_free(my_coeffs)
812 : : end do !iterm
813 : : !enddo ! order
814 : :
815 [ + - ]: 1 : if(allocated(singledisp_terms)) call polynomial_coeff_list_free(singledisp_terms)
816 : :
817 : 1 : write(message, '(a,(80a),a)' ) ch10,&
818 [ + - + + ]: 81 : &('_',ii=1,80),ch10
819 : 1 : call wrtout(ab_out,message,'COLL')
820 : 1 : call wrtout(std_out,message,'COLL')
821 : :
822 : 1 : write(message,'(3a)' )ch10,&
823 : 2 : & ' Finished creating high-order terms',ch10 !,&
824 : : !& ' Optimize initial anharmonic terms !NOT IS COMMENTED NOW!',ch10
825 : 1 : call wrtout(ab_out,message,'COLL')
826 : 1 : call wrtout(std_out,message,'COLL')
827 : : ! call opt_effpot(eff_pot,nterm,terms,hist,comm,print_anh=.FALSE.)
828 : :
829 : : ! Print the standard devition of final model
830 : 1 : write(message,'(6a,ES24.16,6a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,&
831 : 1 : & ' Mean Standard Deviation values of the effective-potential',ch10,&
832 : 1 : & ' with respect to the training-set after attempted bounding (meV^2/atm):',&
833 : 1 : & ch10,' Energy : ',&
834 : 1 : & mse_ini* (Ha_EV*1000)**2 *factor ,ch10,&
835 : 1 : & ' Goal function values of the effective.potential',ch10,&
836 : 1 : & ' with respect to the test-set (eV^2/A^2):',ch10,&
837 : 1 : & ' Forces+Stresses : ',&
838 : 1 : & (msef_ini+mses_ini)*(HaBohr_eVAng)**2,ch10,&
839 : 1 : & ' Forces : ',&
840 : 1 : & msef_ini*(HaBohr_eVAng)**2,ch10,&
841 : 1 : & ' Stresses : ',&
842 : 2 : & mses_ini*(HaBohr_eVAng)**2,ch10
843 : 1 : call wrtout(ab_out,message,'COLL')
844 : 1 : call wrtout(std_out,message,'COLL')
845 : :
846 : :
847 : : !DEALLOCATION
848 [ + - ]: 1 : ABI_SFREE(symbols)
849 [ + - ]: 1 : ABI_SFREE(terms)
850 [ + - ]: 1 : ABI_SFREE(my_weights)
851 : :
852 : : !ABI_FREE(my_coeffs)
853 : 1 : call fit_data_free(fit_data)
854 : 1 : call sympairs%free()
855 : :
856 : :
857 [ - + - + : 1 : end subroutine opt_effpotbound
- + - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - - - -
- - - - -
- - - - -
- - - - -
- - + - -
- - - - -
- - - - -
- - - - -
- - - - +
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + -
+ ]
858 : : !!***
859 : :
860 : :
861 : : !!****f* m_opt_effpot/opt_getHOforterm
862 : : !!
863 : : !! NAME
864 : : !! opt_effpotbound
865 : : !!
866 : : !! FUNCTION
867 : : !! Compute possible high orders for a given anharmonic term.
868 : : !! In the range of order_start,order_stop
869 : : !!
870 : : !!
871 : : !! INPUTS
872 : : !! term<polynomial_coeff_type>:anharmonic term
873 : : !! order_range(2):start and stop order desired by user
874 : : !!
875 : : !!
876 : : !! OUTPUT
877 : : !! order_start:possible start order
878 : : !! order_stop: possible stop order
879 : : !!
880 : : !! SOURCE
881 : :
882 : 2 : subroutine opt_getHOforterm(term,order_range,order_start,order_stop)
883 : :
884 : : !Arguments ------------------------------------
885 : : !scalars
886 : : type(polynomial_coeff_type),intent(in) :: term
887 : : !arrays
888 : : integer,intent(in) :: order_range(2)
889 : : integer,intent(out) :: order_start, order_stop
890 : : !Logicals
891 : : !Strings
892 : : !Local variables ------------------------------
893 : : !scalars
894 : : integer :: idisp,ndisp,nstrain,nterm_of_term,power_tot
895 : : integer :: nbody_tot
896 : : !arrays
897 : 2 : integer,allocatable :: powers(:)
898 : : !Logicals
899 : : !Strings
900 : : !*************************************************************************
901 : :
902 : : !Get/Initialize variables
903 : 2 : ndisp = term%terms(1)%ndisp
904 : 2 : nstrain = term%terms(1)%nstrain
905 : 2 : nbody_tot = ndisp + nstrain
906 : 2 : nterm_of_term = term%nterm
907 : :
908 [ - + + - : 6 : ABI_MALLOC(powers,(nbody_tot))
+ - + - -
+ - + ]
909 [ + + ]: 6 : powers(:ndisp) = term%terms(1)%power_disp
910 [ + + ]: 3 : powers(ndisp+1:) = term%terms(1)%power_strain
911 : 7 : power_tot = 0
912 : : !write(std_out,*) "powers in getHOforterm", powers
913 : :
914 : : !Get rid off odd displacements
915 [ + + ]: 7 : do idisp=1,nbody_tot
916 [ + + + + ]: 7 : if(idisp <= ndisp .and. mod(powers(idisp),2) == 1)then
917 : 2 : powers(idisp) = powers(idisp) + 1
918 [ + + ]: 3 : else if(mod(powers(idisp),2) == 1)then
919 : 1 : powers(idisp) = powers(idisp) + 1
920 : : endif
921 : : enddo !idisp
922 : : ! Count order
923 [ + + ]: 7 : do idisp=1,nbody_tot
924 : 7 : power_tot = power_tot + powers(idisp)
925 : : enddo
926 : : ! Get start and stop order for this term
927 : : ! If term doesn't fit in order range give back order_start = order_stop = 0
928 [ + + + - ]: 2 : if(power_tot >= order_range(1) .and. power_tot <=order_range(2))then
929 : 1 : order_start = power_tot
930 : 1 : order_stop = order_range(2)
931 [ + - ]: 1 : elseif(power_tot < order_range(1))then
932 : 1 : order_start = order_range(1)
933 : 1 : order_stop = order_range(2)
934 [ # # ]: 0 : elseif(power_tot > order_range(2))then
935 : 0 : order_start = 0
936 : 0 : order_stop = 0
937 : : endif
938 : :
939 [ + - ]: 2 : ABI_SFREE(powers)
940 : :
941 : 2 : end subroutine opt_getHOforterm
942 : : !!***
943 : :
944 : :
945 : : !!****f* m_opt_effpot/opt_getCombisforterm
946 : : !!
947 : : !! NAME
948 : : !! opt_getCombisforterm
949 : : !!
950 : : !! FUNCTION
951 : : !! For a given order range: order_start, order_stop
952 : : !! calculate number of total possible combinations ncombi
953 : : !! and calculat combinations per order ncombi_order(i)
954 : : !!
955 : : !!
956 : : !! INPUTS
957 : : !! order_start: start order for bounding terms
958 : : !! order_end: end order for bounding terms
959 : : !! ndisp: number of displacements a given term contains
960 : : !!
961 : : !!
962 : : !! OUTPUT
963 : : !! ncombi: total number of combinations
964 : : !! ncombi_order: array with number of combinations per order
965 : : !!
966 : : !!
967 : : !! SOURCE
968 : :
969 [ - + ]: 2 : subroutine opt_getCombisforterm(order_start,order_end,ndisp,ncombi,ncombi_order)
970 : :
971 : : !Arguments ------------------------------------
972 : : !scalars
973 : : integer,intent(in) :: ndisp
974 : : integer,intent(in) :: order_start, order_end
975 : : !arrays
976 : : integer,intent(out) :: ncombi
977 : : integer,intent(out) :: ncombi_order(:)
978 : : !Logicals
979 : : !Strings
980 : : !Local variables ------------------------------
981 : : !scalars
982 : : integer :: i
983 : : integer :: order,iorder1,iorder2
984 : : !arrays
985 : : !integer
986 : : !Logicals
987 : : !Strings
988 : : character(len=1000) :: message
989 : : !*************************************************************************
990 : :
991 : : !Test
992 [ + - - + ]: 2 : if(mod(order_start,2) /= 0 .or. mod(order_end,2) /= 0)then
993 : : ! Message to Output
994 : 0 : write(message,'(4a)' )ch10,&
995 : 0 : & 'Either start or stop order are not even numbers',ch10,&
996 : 0 : & 'Action: change bound_range in input',ch10
997 : 0 : ABI_ERROR(message)
998 : : endif
999 : :
1000 : : !Initialize Variables
1001 : 2 : i = 0
1002 : 2 : ncombi = 0
1003 [ + + ]: 4 : ncombi_order = 0
1004 : :
1005 : : !Calculate Combinations
1006 [ + - + - ]: 4 : do order=order_start,order_end,2
1007 : 2 : i = i+1
1008 [ - + ]: 4 : if(ndisp == 1)then
1009 : 0 : ncombi = ncombi + 1
1010 : 0 : ncombi_order(i) = 1
1011 : : else
1012 [ + - + + ]: 5 : do iorder1 = 2,order-2*(ndisp-1),2
1013 [ + + ]: 3 : if(ndisp*iorder1 == order)then
1014 : 1 : ncombi = ncombi + 1
1015 : 1 : ncombi_order(i) = ncombi_order(i) + 1
1016 : 1 : cycle
1017 : : endif
1018 [ + + + - ]: 5 : do iorder2=iorder1+2,order-2*(ndisp-1),2
1019 [ + - ]: 3 : if( iorder1 + (ndisp-1)*iorder2 == order)then
1020 : 1 : ncombi = ncombi + ndisp
1021 : 1 : ncombi_order(i) = ncombi_order(i) + ndisp
1022 [ # # ]: 0 : elseif(iorder1 * (ndisp-1) + iorder2 == order)then
1023 : 0 : ncombi = ncombi + ndisp
1024 : 0 : ncombi_order(i) = ncombi_order(i) + ndisp
1025 : : endif
1026 : : enddo !iorder2
1027 : : enddo !iorder1 !
1028 : : endif
1029 : : !write(*,*) 'ncombi(',i,') for order',order,' is:', ncombi_order(i), 'are we happy?'
1030 : : enddo !order
1031 : :
1032 : : !write(*,*) ncombi_order(:)
1033 : : !write(*,*) 'ncombi for term is:', ncombi, 'are we happy?'
1034 : :
1035 : :
1036 : 2 : end subroutine opt_getCombisforterm
1037 : : !!***
1038 : :
1039 : : !!****f* m_opt_effpot/opt_getHoTerms
1040 : : !!
1041 : : !! NAME
1042 : : !! opt_geHoTerms
1043 : : !!
1044 : : !! FUNCTION
1045 : : !! For a term give all possible high order terms
1046 : : !! Attention order_start, order_stop, ncombi and
1047 : : !! ncombi_order have to be calculated before!
1048 : : !!
1049 : : !!
1050 : : !!
1051 : : !! INPUTS
1052 : : !! eff_pot: existing effective potential
1053 : : !! order: order for which bounding terms are generated
1054 : : !!
1055 : : !!
1056 : : !! OUTPUT
1057 : : !! eff_pot new effective potential
1058 : : !!
1059 : : !!
1060 : : !! SOURCE
1061 : :
1062 [ - + - + ]: 2 : subroutine opt_getHoTerms(terms,order_start,order_stop,ndisp,ncombi_order)
1063 : :
1064 : : !Arguments ------------------------------------
1065 : : !scalars
1066 : : integer,intent(in) :: ndisp
1067 : : integer,intent(in) :: order_start, order_stop
1068 : : !arrays
1069 : : integer,intent(in) :: ncombi_order(:)
1070 : : !Logicals
1071 : : type(polynomial_coeff_type),intent(inout) :: terms(:)
1072 : : !Strings
1073 : : !Local variables ------------------------------
1074 : : !scalars
1075 : : integer :: i,icombi,icombi2,icombi_start,icombi_stop,idisp,nterm_of_term
1076 : : integer :: order,iterm_of_term,jdisp,power_tot
1077 : : integer :: jdisp1,jdisp2,sec,nstrain,nbody_tot
1078 : : real(sp) :: to_divide,divider1,divider2,divided
1079 : : !arrays
1080 : : !integer
1081 : : !Logicals
1082 : : logical :: equal_term_done
1083 : : !Strings
1084 : : character(len=1000) :: message
1085 : : !*************************************************************************
1086 : : !Get Variables
1087 : 2 : nterm_of_term = terms(1)%nterm
1088 : 2 : nstrain = terms(1)%terms(1)%nstrain
1089 : 2 : nbody_tot = ndisp + nstrain
1090 : :
1091 : : ! Create all possible combinaions for the specified orders
1092 : : ! If anybody has ever, ever to read and understand this I'm terribly sorry
1093 : : ! ---this will be quite hell---
1094 : 2 : power_tot = 0
1095 : 2 : icombi_start = 1
1096 : 2 : i = 0
1097 : : !write(*,*) "what is ncombi?", ncombi
1098 : : !write(*,*) "what is ncmobi_order", ncombi_order
1099 : :
1100 : : !write(*,*) "order_start", order_start
1101 : : !write(*,*) "order_stop", order_stop
1102 : 2 : order = order_start
1103 : : ! TODO work here icombi start and order counting does not work yet.
1104 [ - + - + ]: 4 : do order=order_start,order_stop,2
1105 : : !write(std_out,*) 'Was I now here?!(order-loop)'
1106 : 2 : i = i + 1
1107 : 2 : icombi_stop = icombi_start + ncombi_order(i) - 1
1108 : 2 : equal_term_done = .FALSE.
1109 : : !write(std_out,*) 'order', order
1110 : : !write(std_out,*) 'icombi_start', icombi_start
1111 : : !write(std_out,*) 'icombi_stop', icombi_stop
1112 : 2 : icombi=icombi_start
1113 [ + + ]: 6 : do while (icombi<=icombi_stop)
1114 : : power_tot = 0
1115 [ + + ]: 7 : do idisp=1,nbody_tot
1116 : : !write(*,*) "what is icombi here?", icombi
1117 : : !write(*,*) "what is icombi_stop here?", icombi_stop
1118 [ + + ]: 7 : if(idisp<=ndisp)then
1119 : 4 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(idisp)
1120 : : else
1121 : 1 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(idisp-ndisp)
1122 : : endif
1123 : : enddo
1124 : : ! Probably have to increase order at same time
1125 : : ! If the term already has the right order we cycle
1126 : : ! we increase icombi_start and icombi to go to the next term in the array
1127 [ + + ]: 4 : if(power_tot == order)then
1128 : 1 : icombi_start = icombi_start + 1
1129 : 1 : icombi = icombi + 1
1130 : : !write(*,*) 'icombi-start in if', icombi_start
1131 : 1 : cycle
1132 : : ! If the term is not already in the right order, we manipulate it until
1133 : : ! it is
1134 : : else
1135 : : jdisp1 = 1
1136 : 3 : sec = 1
1137 : : !write(*,*) 'what is ndisp actually', ndisp
1138 : : ! Treat single permutations from the bottom of the order
1139 : : ! so start from ^2^2^2 and get ^6^2^2,^2^6^2, and ^2^2^6 f.E.
1140 : : ! do loop over displacements
1141 [ + + ]: 3 : do while(jdisp1<=nbody_tot .and. sec < 100)
1142 : : !write(*,*) "what is jdisp1 here?", jdisp1
1143 : : !write(*,*) "what is icombi here?", icombi
1144 : : !write(*,*) "what is icombi_sotp?", icombi_stop
1145 : : ! Increase the order of the displacements
1146 [ + + ]: 50 : do iterm_of_term=1,nterm_of_term
1147 [ + - ]: 50 : if(jdisp1 <= ndisp)then
1148 : : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) =&
1149 : 48 : & terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) + 2
1150 : : else
1151 : : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) =&
1152 : 0 : & terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) + 2
1153 : :
1154 : : endif
1155 : : enddo
1156 : : ! Check total order of term after increase
1157 : : power_tot=0
1158 [ + + ]: 6 : do jdisp2=1,nbody_tot
1159 [ + - ]: 6 : if(jdisp2 <= ndisp)then
1160 : 4 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(jdisp2)
1161 : : else
1162 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(jdisp2-ndisp)
1163 : : endif
1164 : : enddo
1165 : : ! If the term is at the right order do next displacement
1166 : : ! increase icombi and icombi_start to go to next term in array
1167 [ - + ]: 3 : if(power_tot == order)then
1168 : 2 : icombi = icombi + 1
1169 : 2 : icombi_start = icombi_start +1
1170 : : !write(*,*) 'what is icombi_start here', icombi_start
1171 : 2 : jdisp1 = jdisp1 + 1
1172 : : !write(*,*) 'and what is jdisp1?', jdisp1
1173 : :
1174 : : endif
1175 : : enddo!jdisp1
1176 : : !Treat permutations in terms ndisp > 2 and order >= 10
1177 : : !Start f.E. from ^4^4^4 and get ^4^2^4, ^2^4^4,^4^4^2
1178 [ - + - - ]: 1 : if(icombi_stop - icombi > 1 .and. nbody_tot >2 )then
1179 [ # # ]: 0 : do icombi2=icombi,icombi+nbody_tot-1
1180 [ # # ]: 0 : do jdisp=1,nbody_tot
1181 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1182 [ # # ]: 0 : if(jdisp <= ndisp)then
1183 : : terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) = &
1184 : 0 : & terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) +2
1185 : : !write(*,*) "What's the power now?", terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp)
1186 : : else
1187 : : terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = &
1188 : 0 : & terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) +2
1189 : : endif
1190 : : enddo
1191 : : enddo !jdisp
1192 : : enddo
1193 : : jdisp1 = 1
1194 [ # # ]: 0 : do while(jdisp1<=nbody_tot .and. sec < 100)
1195 : 0 : sec = sec + 1
1196 : : !write(*,*) "how often did I go here, hu ?"
1197 : : !write(*,*) "I did at least one displacement"
1198 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1199 [ # # ]: 0 : if(jdisp1 <= ndisp)then
1200 : : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) =&
1201 : 0 : & terms(icombi)%terms(iterm_of_term)%power_disp(jdisp1) + 2
1202 : : else
1203 : : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) =&
1204 : 0 : & terms(icombi)%terms(iterm_of_term)%power_strain(jdisp1-ndisp) + 2
1205 : :
1206 : : endif
1207 : : enddo
1208 : : power_tot=0
1209 [ # # ]: 0 : do jdisp2=1,nbody_tot
1210 [ # # ]: 0 : if(jdisp2 <= ndisp)then
1211 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_disp(jdisp2)
1212 : : else
1213 : 0 : power_tot = power_tot + terms(icombi)%terms(1)%power_strain(jdisp2-ndisp)
1214 : : endif
1215 : : enddo
1216 [ # # ]: 0 : if(power_tot == order)then
1217 : 0 : icombi = icombi + 1
1218 : 0 : icombi_start = icombi_start +1
1219 : : !write(*,*) 'what is icombi_start here', icombi_start
1220 : 0 : jdisp1 = jdisp1 + 1
1221 : : !write(*,*) 'and what is jdisp1?', jdisp1
1222 : : endif
1223 : : enddo!jdisp1
1224 : : ! Message to Output
1225 : : if(sec>100)then
1226 : : write(message,'(4a)' )ch10,&
1227 : : & "You're stuck in a while loop.",ch10,&
1228 : : & 'Action: Contact Abinit Group',ch10
1229 : : ABI_ERROR(message)
1230 : : endif
1231 : : endif! (icombi_stop - icombi)
1232 : : !write(*,*) 'I was here!'
1233 : 1 : to_divide = real(order)
1234 : 1 : divider1 = real(nbody_tot)
1235 : 1 : divided = real(to_divide/divider1)
1236 : 1 : divider2 = real(2)
1237 : : !write(*,*) 'divided', divided, 'divider2', divider2
1238 : : !Treat terms with even power f.E. ^2^2^2^2, ^4^4 etc...
1239 [ + - - - ]: 1 : if(mod(divided,divider2) == 0 .and. .not. equal_term_done .and. nbody_tot > 1)then
1240 : : !write(*,*) "Sometimes I should be here sometimes I shouldn't"
1241 [ # # ]: 0 : do jdisp=1,nbody_tot
1242 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1243 [ # # ]: 0 : if(jdisp<=ndisp)then
1244 : 0 : terms(icombi)%terms(iterm_of_term)%power_disp(jdisp) = order/ndisp
1245 : : else
1246 : 0 : terms(icombi)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = order/ndisp
1247 : : endif
1248 : : enddo
1249 : : enddo !jdisp
1250 [ # # ]: 0 : if(order < order_stop)then
1251 [ # # ]: 0 : do icombi2=icombi+1,icombi+nbody_tot
1252 [ # # ]: 0 : do jdisp=1,nbody_tot
1253 [ # # ]: 0 : do iterm_of_term=1,nterm_of_term
1254 [ # # ]: 0 : if(jdisp<=ndisp)then
1255 : 0 : terms(icombi2)%terms(iterm_of_term)%power_disp(jdisp) = order/ndisp
1256 : : else
1257 : 0 : terms(icombi2)%terms(iterm_of_term)%power_strain(jdisp-ndisp) = order/ndisp
1258 : : endif
1259 : : enddo
1260 : : enddo !jdisp
1261 : : enddo
1262 : : endif
1263 : 0 : equal_term_done = .TRUE.
1264 : 0 : icombi = icombi + 1
1265 : 0 : icombi_start = icombi_start +1
1266 : : endif ! equal term if
1267 : : endif ! power_tot == order
1268 : : enddo !icombination
1269 : : enddo !order
1270 : :
1271 : 2 : end subroutine opt_getHoTerms
1272 : : !!***
1273 : :
1274 : : !!****f* m_opt_effpot/opt_filterdisp
1275 : : !!
1276 : : !! NAME
1277 : : !! opt_opt_filterdisp
1278 : : !!
1279 : : !! FUNCTION
1280 : : !! If a anharmonic term represents a strain-phonon coupling
1281 : : !! delete the strain and only keep the displacement part.
1282 : : !!
1283 : : !!
1284 : : !! INPUTS
1285 : : !! term<polynomial_coeff_type>: anharmonic term to check
1286 : : !! nterm_of_term: number of symmetry equivalent terms for term
1287 : : !!
1288 : : !! OUTPUT
1289 : : !! term<polynomial_coeff_type>: only the displacement part of original term
1290 : : !!
1291 : : !! SOURCE
1292 : :
1293 : 1 : subroutine opt_filterdisp(term,nterm_of_term)
1294 : :
1295 : : !Arguments ------------------------------------
1296 : : !scalars
1297 : : type(polynomial_coeff_type),intent(inout) :: term
1298 : : integer,intent(in) :: nterm_of_term
1299 : : !arrays
1300 : : !Logicals
1301 : : !Strings
1302 : : !Local variables ------------------------------
1303 : : !scalars
1304 : : integer :: iterm_of_term
1305 : : !reals
1306 : : real(dp) :: coeff
1307 : : !arrays
1308 [ + + + + : 242 : type(polynomial_term_type) :: terms(nterm_of_term)
+ + + - +
- + - + -
+ - + - +
- ]
1309 : : !Logicals
1310 : : !Strings
1311 : : !*************************************************************************
1312 : :
1313 [ + + ]: 25 : do iterm_of_term=1,nterm_of_term
1314 : : !terms(iterm_of_term) = term%terms(iterm_of_term)
1315 : :
1316 : : call polynomial_term_init(term%terms(iterm_of_term)%atindx, &
1317 : : & term%terms(iterm_of_term)%cell,&
1318 : : & term%terms(iterm_of_term)%direction,&
1319 : : & term%terms(iterm_of_term)%ndisp, &
1320 : : & term%terms(iterm_of_term)%nstrain, &
1321 : : & terms(iterm_of_term),&
1322 : : & term%terms(iterm_of_term)%power_disp, &
1323 : : & term%terms(iterm_of_term)%power_strain,&
1324 : : & term%terms(iterm_of_term)%strain, &
1325 : : & term%terms(iterm_of_term)%weight, &
1326 : 24 : & check=.TRUE.)
1327 : 24 : terms(iterm_of_term)%nstrain = 0
1328 [ + + ]: 48 : terms(iterm_of_term)%power_strain = 0
1329 [ + + ]: 49 : terms(iterm_of_term)%strain = 0
1330 : : !DMSG(terms(iterm_of_term)%debug_str)
1331 : : enddo
1332 : :
1333 : 1 : call polynomial_coeff_free(term)
1334 : : !Reinitial term
1335 : : !check=.TRUE. checks for duplicate terms
1336 : 1 : call polynomial_coeff_init(coeff,nterm_of_term,term,terms, check=.TRUE.)
1337 : : !DMSG(term%debug_str)
1338 : :
1339 : :
1340 [ + + ]: 25 : do iterm_of_term=1,nterm_of_term
1341 : 25 : call polynomial_term_free(terms(iterm_of_term))
1342 : : enddo
1343 : : !if(nterm_of_term /= term%nterm)then
1344 : : ! write(*,*) "nterm_of_term changed after deleting strain"
1345 : : !endif
1346 : :
1347 : :
1348 : 192 : end subroutine opt_filterdisp
1349 : : !!***
1350 : :
1351 : : !!****f* m_opt_effpot/opt_getHOstrain
1352 : : !!
1353 : : !! NAME
1354 : : !! opt_getHOstrain
1355 : : !!
1356 : : !! FUNCTION
1357 : : !! Get HO anharmnonic strain terms for bounding and add them to list
1358 : : !! of existing terms in effective potential
1359 : : !!
1360 : : !! INPUTS
1361 : : !! eff_pot<effective_potential_type>: datatype with all the information
1362 : : !! about the effective potential
1363 : : !! power_strain(2): start and stop order for strain terms
1364 : : !! comm: mpi communicator (at the moment only sequential tested)
1365 : : !!
1366 : : !! OUTPUT
1367 : : !! terms<polynomial_coeff_type>: list with original terms in effective
1368 : : !! potential + HO even strain terms
1369 : : !!
1370 : : !! SOURCE
1371 : :
1372 [ - + ]: 1 : subroutine opt_getHOstrain(terms,ncombi,nterm_start,eff_pot,power_strain,comm, max_nbody)
1373 : :
1374 : : !Arguments ------------------------------------
1375 : : !scalars
1376 : : integer,intent(in) :: comm
1377 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms(:)
1378 : : type(effective_potential_type), intent(in) :: eff_pot
1379 : : integer,intent(in) :: power_strain(2), max_nbody(:)
1380 : : integer,intent(out) :: ncombi,nterm_start
1381 : : !arrays
1382 : : !Logicals
1383 : : !Strings
1384 : : !Local variables ------------------------------
1385 : : !scalars
1386 : : integer :: nterm_tot_tmp
1387 : : integer :: i,ii
1388 : : !integer :: nbody, ref_nbody
1389 : : real(dp) :: coeff_ini
1390 : : !reals
1391 [ + + + + : 68 : type(crystal_t) :: crystal
+ + + + +
+ + + + +
+ + + + +
+ ]
1392 : : !arrays
1393 : 1 : type(polynomial_coeff_type),allocatable :: strain_terms_tmp(:)
1394 : : !Logicals
1395 : : !Strings
1396 : : character(len=1000) :: message
1397 : : !*************************************************************************
1398 : : !Get variables
1399 [ + - + - : 1 : crystal = eff_pot%crystal
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - ]
1400 : 1 : coeff_ini = 1000000
1401 : :
1402 : 1 : write(message, '(a,(80a),a)' ) ch10,&
1403 [ + - + + ]: 81 : & ('_',ii=1,80),ch10
1404 : 1 : call wrtout(ab_out,message,'COLL')
1405 : 1 : call wrtout(std_out,message,'COLL')
1406 : 1 : write(message,'(3a)' )ch10,&
1407 : 2 : & ' Chreate high order strain terms ',ch10
1408 : 1 : call wrtout(ab_out,message,'COLL')
1409 : 1 : call wrtout(std_out,message,'COLL')
1410 : :
1411 : :
1412 : : !1406 get count of high order even anharmonic strain terms and the strain terms itself
1413 : 1 : call polynomial_coeff_getEvenAnhaStrain(strain_terms_tmp,crystal,ncombi,power_strain,comm, max_nbody)
1414 : : ! Allocate my_coeffs with ncombi free space to work with
1415 : :
1416 : :
1417 : :
1418 : 1 : nterm_start = eff_pot%anharmonics_terms%ncoeff
1419 : 1 : nterm_tot_tmp = eff_pot%anharmonics_terms%ncoeff + ncombi
1420 [ - + + - : 7 : ABI_MALLOC(terms,(nterm_tot_tmp))
+ - + - -
+ - + - +
+ + + - ]
1421 [ + + ]: 5 : do i=1,nterm_tot_tmp
1422 [ + + ]: 5 : if(i<=nterm_start)then
1423 : : call polynomial_coeff_init(coeff_ini,eff_pot%anharmonics_terms%coefficients(i)%nterm,terms(i),&
1424 : : & eff_pot%anharmonics_terms%coefficients(i)%terms,eff_pot%anharmonics_terms%coefficients(i)%name,&
1425 : 2 : & check=.TRUE.)
1426 : : !DMSG(terms(i)%debug_str)
1427 : : else
1428 [ - + ]: 2 : if (.not.allocated(strain_terms_tmp(i-nterm_start)%terms)) then
1429 : 0 : write(message,'(4a)' )ch10,&
1430 : 0 : & 'Error in opt_getHOstrain: strain_terms_tmp(i-nterm_start)%terms not allocated',ch10,&
1431 : 0 : & 'Action: Contact Abinit Group',ch10
1432 : 0 : ABI_ERROR(message)
1433 : : endif
1434 : : call polynomial_coeff_init(coeff_ini,strain_terms_tmp(i-nterm_start)%nterm,terms(i),&
1435 : : & strain_terms_tmp(i-nterm_start)%terms,strain_terms_tmp(i-nterm_start)%name,&
1436 : 2 : & check=.TRUE.)
1437 : : !DMSG(strain_terms_tmp(i-nterm_start)%debug_str)
1438 : : endif
1439 : : enddo
1440 : :
1441 : 1 : call polynomial_coeff_list_free(strain_terms_tmp)
1442 : :
1443 [ - + - - : 1 : end subroutine opt_getHOstrain
- - - - -
- - - - -
- - - - -
- - - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - +
- + - + -
+ - + - ]
1444 : : !!***
1445 : :
1446 : : !!****f* m_opt_effpot/opt_getHOcrossdisp
1447 : : !!
1448 : : !! NAME
1449 : : !! opt_getHOcrossdisp
1450 : : !!
1451 : : !! FUNCTION
1452 : : !! Get even high order displacement terms for a given input term and
1453 : : !! add them to an excisting list of terms. If the term is strain phonon type
1454 : : !! the strain part gets deleted and the high order terms for the displacement
1455 : : !! part are computed.
1456 : : !! Example: for a tree linear term with three displacements x*y*z the even high order
1457 : : !! possibilites are computed and stored.
1458 : : !! 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.
1459 : : !!
1460 : : !! INPUTS
1461 : : !! eff_pot<effective_potential_type>: datatype with all the information
1462 : : !! about the effective potential
1463 : : !! power_disp(2): start and stop order for disp terms
1464 : : !! comm: mpi communicator (at the moment only sequential tested)
1465 : : !!
1466 : : !! OUTPUT
1467 : : !! terms<polynomial_coeff_type>: list with original terms in effective
1468 : : !! potential + HO even disp terms
1469 : : !!
1470 : : !! SOURCE
1471 : :
1472 : 1 : subroutine opt_getHOcrossdisp(terms_out,ncombi,term_in,power_disp)
1473 : :
1474 : : !Arguments ------------------------------------
1475 : : !scalars
1476 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms_out(:)
1477 : : type(polynomial_coeff_type),intent(inout) :: term_in
1478 : : integer,intent(in) :: power_disp(2)
1479 : : integer,intent(out) :: ncombi
1480 : : !arrays
1481 : : !Logicals
1482 : : !Strings
1483 : : !Local variables ------------------------------
1484 : : !scalars
1485 : : integer :: ndisp,nterm_of_term,nstrain,nbody_tot
1486 : : integer :: order_start,order_stop,norder
1487 : : integer :: order_start_str,order_stop_str
1488 : : integer :: icombi,idisp,iterm_of_term
1489 : : integer :: ncombi_tot,ncombi_str
1490 : : !reals
1491 : : real(dp) :: coeff_ini=1
1492 : : !arrays
1493 : 1 : type(polynomial_coeff_type) :: term
1494 : : integer,allocatable :: ncombi_order(:),ncombi_order_str(:)
1495 : : !Logicals
1496 : : logical :: had_strain
1497 : : !Strings
1498 : : character(len=1000) :: message
1499 : : !*************************************************************************
1500 : : !Get/Set Variables
1501 : 1 : norder = abs(((power_disp(2)-power_disp(1))/2)) + 1
1502 [ - + + - : 3 : ABI_MALLOC(ncombi_order,(norder))
+ - - + -
+ ]
1503 [ - + + - : 2 : ABI_MALLOC(ncombi_order_str,(norder))
- + - + ]
1504 [ + + ]: 2 : ncombi_order = 0
1505 [ + + ]: 2 : ncombi_order_str = 0
1506 : 1 : order_start_str = 0
1507 : 1 : order_stop_str = 0
1508 : :
1509 : 1 : ncombi = 0
1510 : : !Get this term (iterm) and infromations about it
1511 : : !Get number of displacements and equivalent terms for this term
1512 : : !Chose term one to get ndisp. ndisp is equal for all terms of the term
1513 : : !Get minimum oder for this term
1514 : : !Get total number of terms in effpot for message
1515 : 1 : ndisp = term_in%terms(1)%ndisp
1516 : 1 : nstrain = term_in%terms(1)%nstrain
1517 : 1 : nbody_tot = ndisp + nstrain
1518 : 1 : nterm_of_term = term_in%nterm
1519 : 1 : call polynomial_coeff_init(coeff_ini,nterm_of_term,term,term_in%terms,term_in%name, check=.true.)
1520 : : !DMSG(term%debug_str)
1521 : : ! Check if term has strain component.
1522 : : ! If yes filter strain and fit high order atomic displacement terms
1523 : 1 : had_strain = .FALSE.
1524 [ + - ]: 1 : if(term%terms(1)%nstrain /= 0)then
1525 : : ! Message to Output
1526 : 1 : write(message,'(5a)' )ch10,&
1527 : 1 : & '- Term has strain compenent',ch10,&
1528 : 2 : & ' -> Filter Displacement',ch10
1529 : 1 : call wrtout(ab_out,message,'COLL')
1530 : 1 : call wrtout(std_out,message,'COLL')
1531 : 1 : call opt_filterdisp(term,nterm_of_term)
1532 : : !Get new value of symmetry equivalent term nterm_of_term
1533 : 1 : nterm_of_term = term%nterm
1534 : : !Remember if this term had strain
1535 : 1 : had_strain = .TRUE.
1536 : : !cycle
1537 : : endif
1538 : : ! Ok we want it. Let's go.
1539 : :
1540 : : ! get start and stop order for this term
1541 : 1 : call opt_getHOforterm(term,power_disp,order_start,order_stop)
1542 [ + - ]: 1 : if(had_strain) call opt_getHOforterm(term_in,power_disp,order_start_str,order_stop_str)
1543 [ - + ]: 1 : if(order_start == 0)then
1544 : : ! Message to Output
1545 : 0 : write(message,'(5a,I2,a,I2,3a)' )ch10,&
1546 : 0 : & " ==> High order cross product terms for term ", trim(term%name),ch10,&
1547 : 0 : & " ==> do not fit into specified order range from ", power_disp(1),' to ',power_disp(2),ch10,&
1548 : 0 : & " ==> Can not construct high order cross product bounding term",ch10
1549 : 0 : call wrtout(ab_out,message,'COLL')
1550 : 0 : call wrtout(std_out,message,'COLL')
1551 [ # # ]: 0 : ABI_SFREE(ncombi_order)
1552 [ # # ]: 0 : ABI_SFREE(ncombi_order_str)
1553 : 0 : return
1554 : : end if
1555 : :
1556 [ - + - - ]: 1 : if(order_start_str == 0 .and. had_strain)then
1557 : : ! Message to Output
1558 : 0 : write(message,'(5a,I2,a,I2,3a)' )ch10,&
1559 : 0 : & " ==> High order cross product terms for term ", trim(term_in%name),ch10,&
1560 : 0 : & " ==> do not fit into specified order range from ", power_disp(1),' to ',power_disp(2),ch10,&
1561 : 0 : & " ==> Can not construct high order cross product bounding term",ch10
1562 : 0 : call wrtout(ab_out,message,'COLL')
1563 : 0 : call wrtout(std_out,message,'COLL')
1564 : 0 : had_strain = .FALSE.
1565 : : end if
1566 : : ! get total amount of combinations and combinations per order for the term
1567 : 1 : call opt_getCombisforterm(order_start,order_stop,ndisp,ncombi,ncombi_order)
1568 [ + - ]: 1 : if(had_strain)then
1569 : 1 : call opt_getCombisforterm(order_start_str,order_stop_str,nbody_tot,ncombi_str,ncombi_order_str)
1570 : : !write(std_out,*) "I was here ncombi_str is: ", ncombi_str
1571 [ - + + - : 6 : ABI_MALLOC(terms_out,(ncombi+ncombi_str))
+ - + - -
+ - + - +
+ + + - ]
1572 : 2 : ncombi_tot = ncombi+ncombi_str
1573 : : else
1574 [ # # # # : 0 : ABI_MALLOC(terms_out,(ncombi))
# # # # #
# # # # #
# # # # ]
1575 : 0 : ncombi_tot = ncombi
1576 : : endif
1577 : : ! Copy current term to the ncombination elemenst a the end in array terms
1578 : : ! change the value of their coefficient to a start value
1579 : : ! The start is estimed to not be larger then half the initial term's value
1580 : : ! This is because higher order terms should have smaller coefficients
1581 [ + + ]: 4 : do icombi=1,ncombi_tot
1582 [ + + ]: 3 : if(icombi <= ncombi)then
1583 : 2 : coeff_ini = 1 !abs(terms_out(icombi)%coefficient / 2)
1584 : 2 : nterm_of_term = term%nterm
1585 : 2 : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out(icombi),term%terms, check=.true.)
1586 : : !DMSG(terms_out(icombi)%debug_str)
1587 : : else
1588 : 1 : coeff_ini = 10d3
1589 : 1 : nterm_of_term = term_in%nterm
1590 : 1 : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out(icombi),term_in%terms, check=.true.)
1591 : : !DMSG(terms_out(icombi)%debug_str)
1592 : : endif
1593 : : ! Set the power of all terms we want to add to two. We find the correct power later
1594 : : ! Change the weight of the term to 1 (even terms have allways weight=1)
1595 [ + + ]: 76 : do iterm_of_term=1,nterm_of_term
1596 : 72 : terms_out(icombi)%terms(iterm_of_term)%weight = 1
1597 [ + + ]: 72 : if(icombi > ncombi)then
1598 [ + + ]: 48 : terms_out(icombi)%terms(iterm_of_term)%power_strain = 2
1599 : 24 : coeff_ini = 10d3 !abs(terms_out(icombi)%coefficient / 2)
1600 : 24 : terms_out(icombi)%coefficient = coeff_ini
1601 : : ! elseif(icombi>2*ncombi)then
1602 : : ! terms_out(icombi)%terms(iterm_of_term)%power_strain = 4
1603 : : ! coeff_ini = 10d5 !abs(terms_out(icombi)%coefficient / 2)
1604 : : ! terms_out(icombi)%coefficient = coeff_ini
1605 : : endif
1606 [ + + ]: 219 : do idisp=1,ndisp
1607 : 216 : terms_out(icombi)%terms(iterm_of_term)%power_disp(idisp) = 2
1608 : : enddo !idisp
1609 : : enddo !iterm_of_term
1610 : : enddo !icombi
1611 : :
1612 : : !If term had strain we had to reinitialize it in the process
1613 : : !Refree memory
1614 : : !if(had_strain)
1615 : 1 : call polynomial_coeff_free(term)
1616 : :
1617 : : ! Get high order combinations
1618 [ + - ]: 1 : if(had_strain)then
1619 : 1 : call opt_getHoTerms(terms_out(:ncombi),order_start,order_stop,ndisp,ncombi_order)
1620 : 1 : call opt_getHoTerms(terms_out(ncombi+1:),order_start_str,order_stop_str,ndisp,ncombi_order_str)
1621 : 1 : ncombi = ncombi_tot
1622 : : else
1623 : 0 : call opt_getHoTerms(terms_out,order_start,order_stop,ndisp,ncombi_order)
1624 : : endif
1625 : : !DEALLOCATION
1626 [ + - ]: 1 : ABI_SFREE(ncombi_order)
1627 [ + - ]: 1 : ABI_SFREE(ncombi_order_str)
1628 : :
1629 [ - + - - : 3 : end subroutine opt_getHOcrossdisp
- - - - -
- - - - -
- - - - -
+ - + ]
1630 : : !!***
1631 : :
1632 : : !!****f* m_opt_effpot/opt_getSingleDispTerms
1633 : : !!
1634 : : !! NAME
1635 : : !! opt_getSingleDispTerms
1636 : : !!
1637 : : !! FUNCTION
1638 : : !! Get polynomial terms (<polynomial_coeff_type>) with single displacements
1639 : : !! at second order inside a given range defined by the supercell size
1640 : : !!
1641 : : !! INPUTS
1642 : : !! sc_size(3): supercell size
1643 : : !! crystal<type(crystal_t)>: all information about the crystal
1644 : : !! comm: mpi communicator (at the moment only sequential tested)
1645 : : !!
1646 : : !! OUTPUT
1647 : : !! terms<polynomial_coeff_type>: list single displacement polynomial_coeffs
1648 : : !!
1649 : : !! SOURCE
1650 : :
1651 : 1 : subroutine opt_getSingleDispTerms(terms,crystal, sc_size,comm)
1652 : :
1653 : : !Arguments ------------------------------------
1654 : : !scalars
1655 : : integer,intent(in) :: comm
1656 : : type(polynomial_coeff_type),allocatable,intent(inout) :: terms(:)
1657 : : type(crystal_t),intent(inout) :: crystal
1658 : : real(dp) :: cutoff
1659 : : !arrays
1660 : : integer :: sc_size(3)
1661 : :
1662 : : !Logicals
1663 : : !Strings
1664 : : !Local variables ------------------------------
1665 : : !scalars
1666 : : integer :: natom,nsym,nrpt,ncoeff_sym,nstr_sym
1667 : : integer :: ncoeff,ncoeff_out,power_strph,option_GN,option
1668 : : integer :: nterms_out,nterm1,iterm1,iterm2,ind,iatom,i
1669 : : integer :: ncopy
1670 : : integer :: ii !,ia,ib,r1,r2,r3
1671 : : !integer :: irpt,irpt_ref
1672 : : integer :: master,nproc,my_rank
1673 : : !arrays
1674 : : integer :: power_disp(2)
1675 : 1 : integer,allocatable :: cell(:,:)
1676 : 1 : integer,allocatable :: list_symcoeff(:,:,:),list_symstr(:,:,:)
1677 : 1 : logical,allocatable :: terms_to_copy(:)
1678 : 1 : type(polynomial_coeff_type),allocatable :: terms_tmp(:),terms_tmp2(:)
1679 : : !type(polynomial_coeff_type),allocatable :: terms(:)
1680 : : !real(dp),allocatable :: xcart(:,:),xred(:,:),rpt(:,:)
1681 : : real(dp) :: rprimd(3,3),range_ifc(3)
1682 : 1 : real(dp),allocatable :: dist(:,:,:,:)
1683 : 1 : character(len=5),allocatable :: symbols(:)
1684 : : !Logicals
1685 : : logical :: iam_master,need_verbose
1686 : : !Strings
1687 : : character(len=1000) :: message
1688 : : !*************************************************************************
1689 : :
1690 : 1 : option = 2
1691 : :
1692 : : if(option == 1)then
1693 : : !MPI variables
1694 : : master = 0
1695 : : nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm)
1696 : : iam_master = (my_rank == master)
1697 : :
1698 : : !Set/Get other variables
1699 : : !natom = crystal%natom
1700 : : !nsym = crystal%nsym
1701 : : rprimd = crystal%rprimd
1702 : : need_verbose = .TRUE.
1703 : :
1704 : :
1705 : : call prepare_for_getList(crystal,sc_size, dist, cell, natom, nsym, nrpt, range_ifc , symbols)
1706 : :
1707 : : ! FIXME: this is dangerous!
1708 : : ! It assumes that the structure is cubic in a usual axis setting.
1709 : : cutoff = rprimd(1,1)
1710 : :
1711 : : if(need_verbose)then
1712 : : write(message,'(1a)')' Generation of the list of all the possible coefficients'
1713 : : call wrtout(std_out,message,'COLL')
1714 : : end if
1715 : : call polynomial_coeff_getList(cell,crystal,dist,list_symcoeff,list_symstr,&
1716 : : & natom,nstr_sym,ncoeff_sym,nrpt,range_ifc,cutoff,sc_size=sc_size)
1717 : :
1718 : : ABI_SFREE(dist)
1719 : : !ABI_FREE(rpt)
1720 : :
1721 : : !write(*,*) "polynomial_getList worked"
1722 : :
1723 : : !Get Order1 term s
1724 : : !Check difference between ncoeff_out, ncoeff
1725 : : call polynomial_coeff_getOrder1(cell,terms,list_symcoeff,natom,nterms_out,ncoeff_sym,nrpt,nsym,symbols)
1726 : :
1727 : : !do ii=1,nterms_out
1728 : : ! write(*,*) "Term(",ii,"/",nterms_out,"): ", terms(ii)%name
1729 : : !enddo
1730 : :
1731 : : elseif(option == 2)then
1732 : : !Get/set variables
1733 : : ! FIXME: hexu: check these values of 2. Why are they are hard coded here.
1734 : : ! Does this mean we only generate high-order terms in 2x2x2 cell?
1735 : 1 : power_disp = (/2,2/)
1736 : 1 : power_strph = zero
1737 : 1 : option_GN = 0
1738 : 1 : sc_size = (/2,2,2/)
1739 : 1 : cutoff = 0
1740 : 1 : natom = crystal%natom
1741 : : ! XXX: hexu: why cutoff = sum of axis lengths?
1742 [ + + ]: 4 : do ii=1,3
1743 : : cutoff = cutoff + sqrt(crystal%rprimd(ii,1)**2 + &
1744 : : & crystal%rprimd(ii,2)**2 + &
1745 : 4 : & crystal%rprimd(ii,3)**2)
1746 : : enddo
1747 : :
1748 : :
1749 : : !do iatom=1,3
1750 [ + + ]: 6 : do iatom=1,natom !Subhadeep
1751 : : ! TODO: to be validated.
1752 : : !do iatom=1,crystal%nirredat
1753 : : call polynomial_coeff_getNorder(terms_tmp,crystal,cutoff,ncoeff,ncoeff_out,power_disp,&
1754 : : & power_strph,option_GN,sc_size,comm,anharmstr=.false.,spcoupling=.false.,&
1755 : : & only_odd_power=.false.,only_even_power=.true.,verbose=.false.,&
1756 : : & compute_symmetric=.false.,fit_iatom=iatom, &
1757 : 5 : & max_nbody=[999,999,999,999,999,999, 999,999,999,999])
1758 : : !TEST MS
1759 : : ! write(std_out,*) "behind call getNorder"
1760 : : ! write(std_out,*) "ncoeff_out: ", ncoeff_out
1761 : : ! do ii=1,ncoeff_out
1762 : : ! write(*,*) "Term(",ii,"/",ncoeff_out,"): ", terms_tmp(ii)%name
1763 : : ! enddo
1764 : : !TEST MS
1765 [ + + ]: 6 : if(iatom == 1)then
1766 [ - + + - : 13 : ABI_MALLOC(terms,(size(terms_tmp)))
+ - + - -
+ - + - +
+ + + - ]
1767 : 1 : call coeffs_list_copy(terms,terms_tmp)
1768 : : else
1769 [ + + + - : 12 : ABI_MALLOC(terms_to_copy,(size(terms_tmp)))
+ - + + -
+ - + -
+ ]
1770 : : ! note: when iatom>1, terms is already initialized.
1771 : 4 : nterm1 = size(terms)
1772 [ - + + - : 145 : ABI_MALLOC(terms_tmp2,(nterm1))
+ - + - -
+ - + - +
+ + ]
1773 : 4 : terms_tmp2 = terms
1774 [ + + ]: 44 : terms_to_copy = .TRUE.
1775 [ + + ]: 44 : do iterm1=1,size(terms_tmp)
1776 [ + + ]: 795 : do iterm2=1,size(terms)
1777 [ - + ]: 791 : if(terms_tmp(iterm1) == terms(iterm2))then
1778 : 0 : terms_to_copy(iterm1) = .FALSE.
1779 : 0 : exit
1780 : : endif
1781 : : enddo!iterm1
1782 : : enddo!iterm2
1783 [ + + + - ]: 44 : ncopy = count(terms_to_copy)
1784 : 4 : call polynomial_coeff_list_free(terms)
1785 [ - + + - : 185 : ABI_MALLOC(terms,(nterm1+ncopy))
+ - + - -
+ - + - +
+ + + - ]
1786 : 4 : call coeffs_list_copy(terms(:nterm1),terms_tmp2)
1787 : 4 : ind = 0
1788 [ + + ]: 44 : do i =1,size(terms_tmp)
1789 [ + - ]: 44 : if(terms_to_copy(i))then
1790 : 40 : ind=ind+1
1791 : : call polynomial_coeff_init(terms_tmp(i)%coefficient,terms_tmp(i)%nterm,terms(nterm1+ind),terms_tmp(i)%terms,&
1792 : 40 : & terms_tmp(i)%name,check=.TRUE.)
1793 : : !DMSG(terms(nterm1+ind)%debug_str)
1794 : : endif
1795 : : enddo!i=1,nterm1+ncopy
1796 : 4 : call polynomial_coeff_list_free(terms_tmp)
1797 : 4 : call polynomial_coeff_list_free(terms_tmp2)
1798 [ - + ]: 4 : ABI_SFREE(terms_to_copy)
1799 : : endif!iatom==1
1800 : : enddo !iatom=1,natom
1801 : : ! call polynomial_coeff_getNorder(terms,crystal,cutoff,ncoeff,ncoeff_out,power_disp,&
1802 : : !& power_strph,option_GN,sc_size,comm,anharmstr=.false.,spcoupling=.false.,&
1803 : : !& only_odd_power=.false.,only_even_power=.true.,verbose=.false.,&
1804 : : !& compute_symmetric=.false.)
1805 : :
1806 : : endif !option
1807 : : !TEST MS
1808 : : ! write(std_out,*) "behind call getNorder"
1809 : : ! write(std_out,*) "ncoeff_out: ", size(terms)
1810 : : ! do ii=1,size(terms)
1811 : : ! write(*,*) "Term(",ii,"/",size(terms),"): ", terms(ii)%name
1812 : : ! enddo
1813 : : !TEST MS
1814 : :
1815 [ - + - + : 1 : end subroutine opt_getSingleDispTerms
- - - - -
- - - - -
- - - - -
- - - - -
- + - - -
- - - - -
- - - - -
- - - - -
- - ]
1816 : : !!***
1817 : :
1818 : : !!****f* m_opt_effpot/opt_getHOSingleDispTerms
1819 : : !!
1820 : : !! NAME
1821 : : !! opt_getHOSingleDispTerms
1822 : : !!
1823 : : !! FUNCTION
1824 : : !! For a given anharmonic term that might consits of product of
1825 : : !! terms find all even high order single displacement terms
1826 : : !! Example: Input term: x*y*z
1827 : : !! Input HO-range: 6 8
1828 : : !! Output terms: x^6,x^8,y^6,y^8,z^6,z^8
1829 : : !! Caution only atomic displacements are taken into account
1830 : : !!
1831 : : !! INPUTS
1832 : : !! term_in<polynomial_coeff_type>: input anharmonic term
1833 : : !! crystal<type(crystal_t)>: all information about the crystal
1834 : : !! single_disp_terms<polynomial_coeff_out>: list of single disp terms at
1835 : : !! second order to select terms from
1836 : : !! power_disp(2): Start and stop power for HO terms
1837 : : !! comm: mpi communicator (at the moment only sequential tested)
1838 : : !!
1839 : : !! OUTPUT
1840 : : !! terms_out<polynomial_coeff_out>: output high order even terms
1841 : : !! ncoeff: number of coefficients
1842 : : !!
1843 : : !! SOURCE
1844 : :
1845 [ - + - + ]: 1 : subroutine opt_getHOSingleDispTerms(term_in,terms_out,symbols,single_disp_terms,power_disp,ncoeff)
1846 : :
1847 : : !Arguments ------------------------------------
1848 : : !scalars
1849 : : integer,intent(out) :: ncoeff
1850 : : type(polynomial_coeff_type),intent(in) :: term_in
1851 : : type(polynomial_coeff_type),intent(in) :: single_disp_terms(:)
1852 : : type(polynomial_coeff_type),allocatable,intent(out) :: terms_out(:)
1853 : : !type(crystal_t),intent(inout) :: crystal
1854 : : !arrays
1855 : : integer, intent(in) :: power_disp(2)
1856 : : character(len=5),intent(in) :: symbols(:)
1857 : : !Logicals
1858 : : !Strings
1859 : : !Local variables ------------------------------
1860 : : !scalars
1861 : : integer :: ndisp,norder, nterm_of_term
1862 : : integer :: icoeff,iorder,idisp, iterm1,iterm2,iterm3
1863 : : real(dp) :: coeff_ini = 1
1864 : : !Strings
1865 : : character(len=200):: name
1866 : : !arrays
1867 : : type(polynomial_coeff_type),allocatable :: terms_out_tmp(:)
1868 : : !Logicals
1869 : 1 : logical,allocatable :: found(:)
1870 : : !*************************************************************************
1871 : : !Get/Set Variables
1872 : : !Number of output terms
1873 : 1 : ndisp = term_in%terms(1)%ndisp
1874 : 1 : norder = abs(power_disp(2)-power_disp(1))/2 + 1
1875 : 1 : ncoeff = norder * ndisp
1876 : : !Allocate output terms
1877 [ - + + - : 5 : ABI_MALLOC(terms_out_tmp,(ncoeff))
+ - + - -
+ - + +
+ ]
1878 : : !find equivalent second order terms in list of single disp terms
1879 : : !for each displacement in input term
1880 : : !Transfer to output term and increase order
1881 : 1 : icoeff = 0
1882 [ + + ]: 3 : do idisp=1,ndisp
1883 [ + + ]: 103 : do iterm1=1,size(single_disp_terms)
1884 [ + + ]: 2346 : do iterm2=1,single_disp_terms(iterm1)%nterm
1885 [ + + + + : 2980 : if(all(term_in%terms(1)%atindx(:,idisp) == single_disp_terms(iterm1)%terms(iterm2)%atindx(:,1)))then
+ + ]
1886 [ + + + + : 606 : if(all(term_in%terms(1)%cell(:,:,idisp) == single_disp_terms(iterm1)%terms(iterm2)%cell(:,:,1)))then
+ + + + ]
1887 [ + + ]: 6 : if(term_in%terms(1)%direction(idisp) == single_disp_terms(iterm1)%terms(iterm2)%direction(1))then
1888 [ + + ]: 4 : do iorder=1,norder
1889 : 2 : icoeff = icoeff + 1
1890 [ + - + - : 32 : terms_out_tmp(icoeff) = single_disp_terms(iterm1)
+ + + - +
- + - + -
+ - + - +
- + - - +
- - - - -
- - - - -
- - - - -
- ]
1891 : 2 : nterm_of_term = single_disp_terms(iterm1)%nterm
1892 : : !Change order of term
1893 : : call polynomial_coeff_init(coeff_ini,nterm_of_term,terms_out_tmp(icoeff),&
1894 : 2 : & single_disp_terms(iterm1)%terms(:), check=.true.)
1895 : : !DMSG(terms_out_tmp(icoeff)%debug_str)
1896 : :
1897 [ + + ]: 34 : do iterm3=1,nterm_of_term
1898 [ + + ]: 62 : terms_out_tmp(icoeff)%terms(iterm3)%power_disp = power_disp(1) + (iorder-1)*2
1899 : : enddo !iterm3
1900 : : enddo !iorder
1901 : : endif
1902 : : endif
1903 : : endif
1904 : : enddo !iterm2
1905 : : enddo !iterm1
1906 : : enddo!idisp
1907 : : !Change Name
1908 [ + + ]: 3 : do icoeff=1,ncoeff
1909 : 2 : call polynomial_coeff_getName(name,terms_out_tmp(icoeff),symbols,recompute=.TRUE.)
1910 : 3 : call polynomial_coeff_SetName(name,terms_out_tmp(icoeff))
1911 : : enddo
1912 : :
1913 : : !Check for doubles and delete them
1914 : : !First count irreducible terms
1915 [ - + + - : 3 : ABI_MALLOC(found,(ncoeff))
+ - + - -
+ - + ]
1916 [ + + ]: 3 : found = .FALSE.
1917 : : iterm3 = 0
1918 [ + + ]: 3 : do iterm1=1,ncoeff
1919 [ + + ]: 4 : do iterm2=iterm1+1,ncoeff
1920 [ - + - - ]: 3 : if(terms_out_tmp(iterm1) == terms_out_tmp(iterm2) .and. .not. found(iterm2))then
1921 : 0 : found(iterm2) = .TRUE.
1922 : 0 : iterm3 = iterm3 + 1
1923 : : endif
1924 : : enddo
1925 : : enddo
1926 : 1 : iterm3 = ncoeff - iterm3
1927 [ - + + - : 5 : ABI_MALLOC(terms_out,(iterm3))
+ - + - -
+ - + - +
+ + + - ]
1928 : 1 : iterm3 = 0
1929 [ + + ]: 3 : do iterm1=1,ncoeff
1930 [ + - ]: 3 : if(.not. found(iterm1))then
1931 : 2 : iterm3 = iterm3 + 1
1932 : : call polynomial_coeff_init(coeff_ini,terms_out_tmp(iterm1)%nterm,terms_out(iterm3),&
1933 : 2 : & terms_out_tmp(iterm1)%terms,terms_out_tmp(iterm1)%name, check=.TRUE.)
1934 : : !DMSG(terms_out(iterm3)%debug_str)
1935 : : endif
1936 : : enddo
1937 [ + - ]: 1 : ABI_SFREE(found)
1938 : 1 : call polynomial_coeff_list_free(terms_out_tmp)
1939 : 1 : ncoeff = iterm3
1940 [ - + - - : 1 : end subroutine opt_getHOSingleDispTerms
- - - - -
- - - - -
- - - - -
- - - ]
1941 : : !!***
1942 : :
1943 : :
1944 : :
1945 : : !!****f* m_opt_effpot/opt_boundcoeff
1946 : : !! NAME
1947 : : !! opt_boundcoeff
1948 : : !!
1949 : : !! FUNCTION
1950 : : !!
1951 : : !! optimize a bound coefficient if optimized value is negative
1952 : : !! put an positive value that respects a precision penalty
1953 : : !!
1954 : : !! INPUTS
1955 : : !!
1956 : : !! OUTPUT
1957 : : !!
1958 : : !! SOURCE
1959 : :
1960 : 1 : function opt_boundcoeff(yvalues,cvalues,penalty_in) result (coeff)
1961 : : !Arguments ------------------------------------
1962 : :
1963 : : !Arguments ------------------------------------
1964 : : real(dp),intent(in) :: yvalues(2),cvalues(2),penalty_in
1965 : : real(dp) :: coeff
1966 : : !local
1967 : : !variable
1968 : : real(dp) :: a,b,coeff_tmp,x1,x2,penalty
1969 : : !array
1970 : : ! *************************************************************************
1971 : :
1972 : 1 : a = ( (yvalues(1) - 1) - (yvalues(2)-1)*(cvalues(1)/cvalues(2))) / (cvalues(1)**2 - cvalues(1)*cvalues(2))
1973 : :
1974 : : b = ( (yvalues(2) - 1)/cvalues(2) ) - ( (yvalues(1) -1)*cvalues(2) - (yvalues(2) - 1)*cvalues(1) )&
1975 : 1 : & / (cvalues(1)**2 - cvalues(1)*cvalues(2))
1976 : :
1977 : 1 : penalty = penalty_in - 1
1978 [ + - + - ]: 1 : if (a/=a .or. b/=b .or. a==0) then
1979 : 1 : coeff = 0.0_dp
1980 : : return
1981 : : end if
1982 : 1 : coeff_tmp = -b/(2*a)
1983 : : !write(*,*) "coeff_tmp", coeff_tmp
1984 [ + - ]: 1 : if(coeff_tmp > 0)then
1985 : : coeff = coeff_tmp
1986 : : else
1987 : 1 : x1 = (-b + sqrt(b**2 + 4*a*penalty)) / (2*a) ! 1.001 penalty value
1988 : 1 : x2 = (-b - sqrt(b**2 + 4*a*penalty)) / (2*a)
1989 [ - + ]: 1 : if(x1>0)then
1990 : : coeff = x1
1991 : : else
1992 : 0 : coeff = x2
1993 : : endif
1994 : : endif
1995 : :
1996 : : end function opt_boundcoeff
1997 : : !!***
1998 : :
1999 : : !!****f* m_opt_effpot/check_to_skip
2000 : : !! NAME
2001 : : !! check_to_skip
2002 : : !!
2003 : : !! FUNCTION
2004 : : !!
2005 : : !! Check if term contains only bodies with even power
2006 : : !! and has a positive coefficient. If yes term doesn't need
2007 : : !! a bounding high order equivalent and we can skip it.
2008 : : !! Function retursn logical to_skip.
2009 : : !!
2010 : : !! INPUTS
2011 : : !!
2012 : : !! term<polynomial_coeff_type>:anharmonic term
2013 : : !!
2014 : : !! OUTPUT
2015 : : !!
2016 : : !! logical: to_skip
2017 : : !!
2018 : : !! SOURCE
2019 : :
2020 : 1 : function check_to_skip(term) result (to_skip)
2021 : : !Arguments ------------------------------------
2022 : : type(polynomial_coeff_type),intent(in) :: term
2023 : : logical :: to_skip
2024 : : ! ------------------------------------
2025 : : !local
2026 : : !variable
2027 : : character(len=1000) :: message
2028 : : !array
2029 : : ! *************************************************************************
2030 : :
2031 : 1 : to_skip = .FALSE.
2032 : : ! Let's check if we really want all this mess
2033 : : ! 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
2034 [ + - + + : 2 : if(term%coefficient > 0 .and. .not. any(mod(term%terms(1)%power_disp(:),2) /= 0))then
- + - - ]
2035 [ # # # # : 0 : if(.not. any(mod(term%terms(1)%power_strain(:),2) /= 0))then
# # ]
2036 : : ! Message to Output
2037 : 0 : write(message,'(3a)' )ch10,&
2038 : 0 : & ' ==> No need for high order bounding term',ch10
2039 : 0 : call wrtout(ab_out,message,'COLL')
2040 : 0 : call wrtout(std_out,message,'COLL')
2041 : 0 : to_skip = .TRUE.
2042 : 0 : return
2043 : : end if
2044 : : end if
2045 : :
2046 : : end function check_to_skip
2047 : : !!***
2048 : :
2049 : 2 : subroutine generate_bounding_term_and_add_to_list(sympairs, nterm_start, ncombi, my_coeffs, temp_cntr)
2050 : : ! check if myterm
2051 : : !type(polynomial_coeff_type),target, intent(in) :: terms(:)
2052 : :
2053 : : type(polynomial_coeff_type), allocatable, intent(inout) :: my_coeffs(:)
2054 : : type(SymPairs_t), intent(inout) :: sympairs
2055 : : integer, intent(in) :: nterm_start, ncombi
2056 : : integer, intent(out) :: temp_cntr
2057 : :
2058 : 2 : type(polynomial_term_type), allocatable :: terms2(:)
2059 : : type(polynomial_term_type), pointer :: myterm=> null()
2060 : : type(polynomial_coeff_type), target, allocatable :: my_coeffs_tmp(:)
2061 : : integer :: idisp, my_nrpt, list_cntr, number_coeff, icombi
2062 : : integer :: jterm
2063 : : integer :: tmp_list(5), pwr
2064 : : logical :: found
2065 : : logical :: reverse_i, any_exists
2066 : 2 : logical, allocatable :: reverse(:), false_reverse(:), exists(:)
2067 : : integer :: counter, nterm, tot_power
2068 : 2 : type(polynomial_coeff_type):: temp_coeff
2069 : 2 : integer, allocatable :: list_disp(:)
2070 : : character(len=200):: name
2071 : : integer :: ncoeff
2072 : :
2073 : : ! copy the coefficients to a temporary array my_coeffs_tmp
2074 : 2 : ncoeff=size(my_coeffs)
2075 : :
2076 : :
2077 [ - + + - : 12 : ABI_MALLOC(my_coeffs_tmp,(ncoeff))
+ - + - -
+ - + +
+ ]
2078 : :
2079 : : ! check if my_coeffs i allocated and has the right size
2080 [ - + ]: 2 : if(.not. allocated(my_coeffs)) then
2081 : 0 : ABI_BUG("my_coeffs is not allocated")
2082 : : end if
2083 [ - + ]: 2 : if(size(my_coeffs) /= ncoeff) then
2084 : 0 : ABI_BUG("my_coeffs has the wrong size")
2085 : : end if
2086 [ - + ]: 2 : if(size(my_coeffs_tmp) /= ncoeff) then
2087 : 0 : ABI_BUG("my_coeffs_tmp has the wrong size")
2088 : : end if
2089 : :
2090 : 2 : my_coeffs_tmp=my_coeffs
2091 : :
2092 : : !call coeffs_list_copy(my_coeffs_tmp, my_coeffs)
2093 : :
2094 : : ! free the original array and allocate a new one with the correct size
2095 : 2 : call polynomial_coeff_list_free(my_coeffs)
2096 [ - + + - : 10 : ABI_MALLOC(my_coeffs, (ncoeff))
- + - + -
+ + + +
- ]
2097 : : !my_coeffs(1:nterm_start) = my_coeffs_tmp(1:nterm_start)
2098 : 2 : call coeffs_list_copy(my_coeffs(1:nterm_start), my_coeffs_tmp(1:nterm_start))
2099 : :
2100 : :
2101 : :
2102 : 2 : temp_cntr =0
2103 [ + + ]: 5 : do icombi=1, ncombi
2104 [ + - ]: 3 : if(size( my_coeffs_tmp(nterm_start+icombi)%terms) == 0) cycle
2105 : 3 : myterm=> my_coeffs_tmp(nterm_start+icombi)%terms(1)
2106 [ + + + + ]: 8 : tot_power=sum(myterm%power_disp) + sum(myterm%power_strain)
2107 [ - + + - : 9 : ABI_MALLOC(list_disp, (tot_power))
+ - + - -
+ - + -
+ ]
2108 [ - + + - : 6 : ABI_MALLOC(reverse, (tot_power))
- + - + -
+ ]
2109 [ - + + - : 6 : ABI_MALLOC(false_reverse, (tot_power))
- + - + -
+ ]
2110 : :
2111 [ + + ]: 21 : list_disp(:)=0
2112 [ + + ]: 21 : reverse(:)=.False.
2113 [ + + ]: 21 : false_reverse(:)=.False.
2114 : 3 : counter=0
2115 : : ! put the indices of the disp terms into a list
2116 [ + + ]: 5 : do idisp=1,myterm%ndisp
2117 : : ! for each disp, find its index in the pair list.
2118 : 2 : my_nrpt=find_irpt(cells=sympairs%cell, cell=myterm%cell(:, 2, idisp))
2119 : : tmp_list(:) = (/myterm%direction(idisp),myterm%atindx(1,idisp), &
2120 [ + + ]: 12 : & myterm%atindx(2,idisp),my_nrpt,int(myterm%weight) /)
2121 : 2 : found=.False. ! to find if the symmetry adapted term already exist.
2122 : 2 : reverse_i=.False.
2123 [ + + ]: 1518 : do list_cntr=1,size(sympairs%list_symcoeff, 2)
2124 : : !do isym=1, sympairs%nsym
2125 [ + + + + : 2220 : if (all(tmp_list(:4)==sympairs%list_symcoeff(:4,list_cntr,1))) then
+ + ]
2126 : 2 : number_coeff = sympairs%list_symcoeff(6,list_cntr,1)
2127 : 2 : found = .True.
2128 : : !if (found) exit
2129 : : end if
2130 : : !end do
2131 : : !if(found) exit
2132 : : end do
2133 [ - + ]: 2 : if(.not. found) then
2134 : 0 : reverse_i=.True.
2135 [ # # # # ]: 0 : my_nrpt=find_irpt(cells=sympairs%cell, cell=-myterm%cell(:, 2, idisp))
2136 : : tmp_list(:) = (/myterm%direction(idisp),myterm%atindx(2,idisp), &
2137 [ # # ]: 0 : & myterm%atindx(1,idisp),my_nrpt,-int(myterm%weight) /)
2138 [ # # ]: 0 : do list_cntr=1,size(sympairs%list_symcoeff, 2)
2139 : : !do isym=1, sympairs%nsym
2140 [ # # # # : 0 : if (all(tmp_list(:4)==sympairs%list_symcoeff(:4,list_cntr,1))) then
# # ]
2141 : 0 : number_coeff = sympairs%list_symcoeff(6,list_cntr,1)
2142 : 0 : found = .True.
2143 : : !if (found) exit
2144 : : end if
2145 : : !end do
2146 : : !if(found) exit
2147 : : end do
2148 : : end if
2149 : :
2150 [ - - ]: 3 : if (found )then
2151 [ + + ]: 6 : do pwr=1,myterm%power_disp(idisp)
2152 : 4 : counter=counter+1
2153 : 4 : list_disp(counter) = number_coeff
2154 : 6 : reverse(counter) = reverse_i
2155 : : end do
2156 : : else
2157 : 0 : ABI_BUG("The pair is not found in generated bounding terms.")
2158 : : end if
2159 : : end do
2160 : :
2161 [ + + ]: 6 : do idisp=1,myterm%nstrain
2162 [ + + ]: 20 : do pwr=1,myterm%power_strain(idisp)
2163 : 14 : counter = counter+1
2164 : 17 : list_disp(counter) = myterm%strain(idisp)+size(sympairs%list_symcoeff,2)
2165 : : end do
2166 : : end do
2167 : 3 : nullify(myterm)
2168 : :
2169 : : ! generate the terms with false reverse, and save to the temp_coeff.
2170 : 3 : nterm=sympairs%nsym
2171 [ - + + - : 153 : ABI_MALLOC(terms2, (nterm))
+ - + - -
+ - + - +
+ + ]
2172 : 3 : call sympairs%generateTerms(list_disp, tot_power, nterm, terms2, reverse=false_reverse)
2173 : : call polynomial_coeff_init(one,nterm,temp_coeff, &
2174 : 3 : &terms2(1:nterm), check=.true.)
2175 : 3 : call polynomial_term_list_free(terms2)
2176 : : !DMSG(temp_coeff%debug_str)
2177 : :
2178 : : call polynomial_coeff_getName(name, &
2179 : 3 : & temp_coeff,sympairs%symbols,recompute=.TRUE.)
2180 : :
2181 : : ! check if temp coeff already exists in the list of coeffs my_coeffs
2182 [ + - - + : 9 : ABI_MALLOC(exists, (nterm_start+icombi-1))
+ - + - +
- - + -
+ ]
2183 [ + + ]: 9 : exists=.False.
2184 [ + + ]: 9 : do jterm=1,(nterm_start+icombi-1)
2185 : 9 : exists(jterm) = coeffs_compare(my_coeffs(jterm),temp_coeff)
2186 : : enddo !jterm
2187 : 3 : call polynomial_coeff_free(temp_coeff)
2188 : :
2189 [ + + + - ]: 9 : any_exists=any(exists)
2190 [ - + ]: 3 : ABI_FREE(exists)
2191 : :
2192 : : ! if the term already exists, we skip it. Otherwise we add also the reverse.
2193 [ + - ]: 3 : if (.not. any_exists) then
2194 : 3 : temp_cntr = temp_cntr+1
2195 [ - + + - : 153 : ABI_MALLOC(terms2, (nterm))
+ - + - -
+ - + - +
+ + ]
2196 : 3 : call sympairs%generateTerms(list_disp, tot_power, nterm, terms2, reverse=reverse)
2197 : : call polynomial_coeff_init(one,nterm, &
2198 : 3 : & my_coeffs(nterm_start+temp_cntr),terms2, check=.true.)
2199 : : !DMSG(my_coeffs(nterm_start+temp_cntr)%debug_str)
2200 : 3 : call polynomial_term_list_free(terms2)
2201 : : end if
2202 : :
2203 [ + - ]: 3 : ABI_SFREE(list_disp)
2204 [ + - ]: 3 : ABI_SFREE(reverse)
2205 [ - + ]: 5 : ABI_SFREE(false_reverse)
2206 : : end do
2207 : :
2208 : 2 : call polynomial_coeff_list_free(my_coeffs_tmp)
2209 : :
2210 : : ! truncate my_coeffs.
2211 : 2 : call coeffs_list_truncate(my_coeffs, nterm_start+temp_cntr)
2212 [ - + - - : 2 : end subroutine generate_bounding_term_and_add_to_list
- - - - -
- - - - -
- - - - -
+ - - - -
- - - - -
- - - - -
- - - + -
- - - - -
- - - - -
- - - - -
- - - - -
+ - + ]
2213 : :
2214 : :
2215 : : end module m_opt_effpot
2216 : : !!***
|
