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