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 56 : 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 : !!***
|
