TABLE OF CONTENTS
- ABINIT/m_fit_polynomial_coeff
- m_fit_polynomial_coeff/fit_polynomial_coeff_computeGF
- m_fit_polynomial_coeff/fit_polynomial_coeff_computeMSD
- m_fit_polynomial_coeff/fit_polynomial_coeff_fit
- m_fit_polynomial_coeff/fit_polynomial_coeff_getCoeffBound
- m_fit_polynomial_coeff/fit_polynomial_coeff_getFS
- m_fit_polynomial_coeff/fit_polynomial_coeff_getPositive
- m_fit_polynomial_coeff/fit_polynomial_coeff_solve
- m_fit_polynomial_coeff/fit_polynomial_printSystemFiles
- m_fit_polynomiaL_coeff/testEffPot
ABINIT/m_fit_polynomial_coeff [ Modules ]
NAME
m_fit_polynomial_coeff
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
18 #if defined HAVE_CONFIG_H 19 #include "config.h" 20 #endif 21 22 #include "abi_common.h" 23 24 module m_fit_polynomial_coeff 25 26 use defs_basis 27 use m_errors 28 use m_abicore 29 use m_polynomial_coeff 30 use m_atomdata 31 use m_xmpi 32 use m_supercell 33 34 use m_special_funcs,only : factorial 35 use m_geometry, only : xred2xcart 36 use m_crystal,only : symbols_crystal 37 use m_strain,only : strain_type,strain_get 38 use m_effective_potential,only : effective_potential_type, effective_potential_evaluate 39 use m_effective_potential,only : effective_potential_freeCoeffs,effective_potential_setCoeffs 40 use m_effective_potential,only : effective_potential_getDisp, effective_potential_writeAnhHead 41 use m_effective_potential,only : effective_potential_copy,effective_potential_free,effective_potential_init 42 use m_effective_potential_file, only : effective_potential_file_mapHistToRef 43 use m_io_tools, only : open_file,get_unit 44 use m_abihist, only : abihist,abihist_free,abihist_init,abihist_copy,write_md_hist,var2hist 45 use m_random_zbq 46 use m_fit_data 47 use m_geometry, only: metric 48 use m_scup_dataset 49 #if defined DEV_MS_SCALEUP 50 use scup_global, only : global_set_parent_iter,global_set_print_parameters 51 #endif 52 53 implicit none 54 55 public :: fit_polynomial_coeff_computeGF 56 public :: fit_polynomial_coeff_computeMSD 57 public :: fit_polynomial_coeff_fit 58 public :: fit_polynomial_coeff_getFS 59 public :: fit_polynomial_coeff_getPositive 60 public :: fit_polynomial_coeff_getCoeffBound 61 public :: fit_polynomial_coeff_solve 62 public :: fit_polynomial_coeff_testEffPot 63 public :: fit_polynomial_printSystemFiles 64 public :: genereList
m_fit_polynomial_coeff/fit_polynomial_coeff_computeGF [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_computeGF
FUNCTION
Compute the values of the goal function (Mean squared error) for gf_value(1) = stresses+forces (Ha/Bohr)**2 gf_value(2) = forces (Ha/Bohr)**2 gf_value(3) = stresses (Ha/Bohr)**2 gf_value(4) = energy (Ha)**2
INPUTS
coefficients(ncoeff) = type(polynomial_coeff_type) energy_coeffs(ncoeff,ntime) = value of the energy for each coefficient (Ha) energy_diff(ntime) = Difference of energ ybetween DFT calculation and fixed part of the model (more often harmonic part) fixed part of the model (more often harmonic part) fcart_coeffs(ncoeff,3,natom,ntime) = value of the forces for each coefficient (-1 factor is taking into acount) (Ha/Bohr) fcart_diff(3,natom,ntime) = Difference of cartesian forces between DFT calculation and fixed part of the model (more often harmonic part) list_coeffs(ncoeff_fit) = List with the indexes of the coefficients used for this model natom = Number of atoms ncoeff_fit = Number of coefficients fitted ncoeff_max = Maximum number of coeff in the list ntime = Number of time in the history strten_coeffs(ncoeff,3,natom,ntime)= value of the stresses for each coefficient (1/ucvol factor is taking into acount) (Ha/Bohr^3) strten_diff(6,natom) = Difference of stress tensor between DFT calculation and fixed part of the model (more often harmonic part) sqomega = Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]
OUTPUT
gf_value(4) = Goal function
SOURCE
2115 subroutine fit_polynomial_coeff_computeGF(coefficients,energy_coeffs,energy_diff,& 2116 & fcart_coeffs,fcart_diff,gf_value,list_coeffs,& 2117 & natom,ncoeff_fit,ncoeff_max,ntime,strten_coeffs,& 2118 & strten_diff,sqomega) 2119 2120 implicit none 2121 2122 !Arguments ------------------------------------ 2123 !scalars 2124 integer,intent(in) :: natom,ncoeff_fit,ncoeff_max,ntime 2125 !arrays 2126 integer,intent(in) :: list_coeffs(ncoeff_fit) 2127 real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime) 2128 real(dp),intent(in) :: energy_diff(ntime) 2129 real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime) 2130 real(dp),intent(in) :: fcart_diff(3,natom,ntime) 2131 real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max) 2132 real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime) 2133 real(dp),intent(in) :: coefficients(ncoeff_fit) 2134 real(dp),intent(out) :: gf_value(4) 2135 !Local variables------------------------------- 2136 !scalar 2137 integer :: ia,icoeff,icoeff_tmp,itime,mu 2138 real(dp):: etmp,emu,fmu,ftmp,smu,stmp 2139 real(dp) :: ffact,sfact,efact 2140 !arrays 2141 ! ************************************************************************* 2142 2143 !1-Compute the value of the goal function 2144 ! see equation 9 of PRB 95 094115(2017) [[cite:Escorihuela-Sayalero2017]] 2145 gf_value = zero 2146 etmp = zero 2147 ftmp = zero 2148 stmp = zero 2149 2150 !Compute factors 2151 ffact = one/(3*natom*ntime) 2152 sfact = one/(6*ntime) 2153 efact = one/(ntime) 2154 2155 ! loop over the configuration 2156 do itime=1,ntime 2157 ! Fill energy 2158 emu = zero 2159 do icoeff=1,ncoeff_fit 2160 icoeff_tmp = list_coeffs(icoeff) 2161 emu = emu + coefficients(icoeff)*energy_coeffs(icoeff_tmp,itime) 2162 end do 2163 ! uncomment the next line to be consistent with the definition of the goal function 2164 etmp = etmp + (energy_diff(itime)-emu)**2/(sqomega(itime)**(1.0/2.0)) 2165 ! uncomment the next get a measure in Ha instead of Ha^2 2166 ! etmp = etmp + abs(energy_diff(itime)-emu) 2167 ! Fill forces 2168 do ia=1,natom 2169 do mu=1,3 2170 fmu = zero 2171 do icoeff=1,ncoeff_fit 2172 icoeff_tmp = list_coeffs(icoeff) 2173 fmu = fmu + coefficients(icoeff)*fcart_coeffs(mu,ia,icoeff_tmp,itime) 2174 end do 2175 ftmp = ftmp + (fcart_diff(mu,ia,itime)-fmu)**2 2176 end do !End loop dir 2177 end do !End loop natom 2178 do mu=1,6 2179 smu = zero 2180 do icoeff=1,ncoeff_fit 2181 icoeff_tmp = list_coeffs(icoeff) 2182 smu = smu + coefficients(icoeff)*strten_coeffs(mu,itime,icoeff_tmp) 2183 end do 2184 stmp = stmp + sqomega(itime)*(strten_diff(mu,itime)-smu)**2 2185 end do !End loop stress dir 2186 end do ! End loop time 2187 2188 gf_value(1) = ffact*ftmp + sfact*stmp !+ efact*etmp !Stresses + Forces 2189 gf_value(2) = ffact*ftmp ! only Forces 2190 gf_value(3) = sfact*stmp ! only Stresses 2191 gf_value(4) = efact*etmp !abs(Energy) 2192 2193 end subroutine fit_polynomial_coeff_computeGF
m_fit_polynomial_coeff/fit_polynomial_coeff_computeMSD [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_computeMSD
FUNCTION
Compute the Mean square error of the energy, forces and stresses
INPUTS
eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = The history of the MD natom = number of atom ntime = number of time in the hist sqomega = Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]] compute_anharmonic = TRUE if the anharmonic part of the effective potential has to be taking into acount print_file = if True, a ASCII file with the difference in energy will be print
OUTPUT
mse = Mean square error of the energy (Hatree) msef = Mean square error of the forces (Hatree/Bohr)**2 mses = Mean square error of the stresses (Hatree/Bohr)**2
SOURCE
2507 subroutine fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,sqomega,comm,& 2508 & compute_anharmonic,print_file,filename,scup_dtset,prt_ph) 2509 2510 implicit none 2511 2512 !Arguments ------------------------------------ 2513 !scalars 2514 integer, intent(in) :: natom,ntime,comm 2515 real(dp),intent(out):: mse,msef,mses 2516 logical,optional,intent(in) :: compute_anharmonic,print_file,prt_ph 2517 !arrays 2518 real(dp),intent(in) :: sqomega(ntime) 2519 type(effective_potential_type),intent(in) :: eff_pot 2520 type(abihist),intent(in) :: hist 2521 !Strings/Characters 2522 character(len=fnlen),optional,intent(in) :: filename 2523 type(scup_dtset_type),optional,intent(inout) :: scup_dtset 2524 !Local variables------------------------------- 2525 !scalar 2526 integer :: ii,ia,mu,unit_energy,unit_stress,itime,master,nproc,my_rank,i 2527 !Uncommend for dipdip test 2528 integer :: ifirst 2529 real(dp):: energy,energy_harm 2530 logical :: need_anharmonic,need_print,need_elec_eval,iam_master 2531 logical :: need_prt_ph 2532 !arrays 2533 real(dp):: fcart(3,natom),gred(3,natom),strten(6),rprimd(3,3),xred(3,natom) 2534 !Strings/Characters 2535 character(len=fnlen) :: file_energy, file_stress, file_anh, name_file 2536 character(len=500) :: msg 2537 !Uncommend for dipdip test 2538 type(abihist) :: hist_out 2539 character(len=200) :: filename_hist 2540 2541 ! ************************************************************************* 2542 !MS Hide SCALE-UP variables 2543 ABI_UNUSED(itime) 2544 2545 !MPI 2546 master = 0 2547 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 2548 iam_master = (my_rank == master) 2549 2550 !Do some checks 2551 if(ntime /= hist%mxhist)then 2552 write(msg,'(a)')'ntime is not correct' 2553 ABI_BUG(msg) 2554 end if 2555 2556 if(natom /= size(hist%xred,2)) then 2557 write(msg,'(a)')'natom is not correct' 2558 ABI_BUG(msg) 2559 end if 2560 2561 need_anharmonic = .TRUE. 2562 if(present(compute_anharmonic))then 2563 need_anharmonic = compute_anharmonic 2564 end if 2565 2566 name_file='' 2567 if(present(filename))name_file = filename 2568 2569 need_print=.FALSE. 2570 if(present(print_file))need_print=print_file 2571 2572 need_elec_eval = .FALSE. 2573 if(present(scup_dtset))need_elec_eval=scup_dtset%scup_elec_model 2574 2575 need_prt_ph=.FALSE. 2576 if(present(prt_ph))need_prt_ph=prt_ph 2577 2578 2579 if(need_print .and. present(filename))then 2580 !MS hist out uncommented for PHONOPY test 2581 call abihist_init(hist_out,natom,ntime,.false.,.false.) 2582 file_energy=trim(name_file)//'_energy.dat' 2583 unit_energy = get_unit() 2584 if (open_file(file_energy,msg,unit=unit_energy,form="formatted",& 2585 & status="unknown",action="write") /= 0) then 2586 ABI_ERROR(msg) 2587 end if 2588 unit_stress = get_unit() 2589 file_stress=trim(name_file)//'_stress.dat' 2590 if (open_file(file_stress,msg,unit=unit_stress,form="formatted",& 2591 & status="unknown",action="write") /= 0) then 2592 ABI_ERROR(msg) 2593 end if 2594 else if(need_print .and. .not. present(filename))then 2595 write(msg,'(3a)')' You asked for printing of the MSD-values',ch10,& 2596 & ' without specifying a filename' 2597 ABI_ERROR(msg) 2598 end if 2599 2600 file_anh=trim(name_file)//'_anharmonic_terms_energy.dat' 2601 2602 mse = zero 2603 msef = zero 2604 mses = zero 2605 do ii=1,ntime ! Loop over configurations 2606 xred(:,:) = hist%xred(:,:,ii) 2607 rprimd(:,:) = hist%rprimd(:,:,ii) 2608 itime = ii 2609 #if defined DEV_MS_SCALEUP 2610 !Pass print options to scale-up 2611 itime = ii 2612 if(need_elec_eval)then 2613 call global_set_parent_iter(itime) 2614 ! Set all print options to false. 2615 call global_set_print_parameters(geom=.FALSE.,eigvals=.FALSE.,eltic=.FALSE.,& 2616 & orbocc=.FALSE.,bands=.FALSE.) 2617 if(ii == 1 .or. modulo(ii,scup_dtset%scup_printniter) == 0)then 2618 call global_set_print_parameters(scup_dtset%scup_printgeom,scup_dtset%scup_printeigv,scup_dtset%scup_printeltic,& 2619 & scup_dtset%scup_printorbocc,scup_dtset%scup_printbands) 2620 end if 2621 end if 2622 #endif 2623 call effective_potential_evaluate(eff_pot,energy_harm,fcart,gred,strten,natom,rprimd,& 2624 & xred=xred,compute_anharmonic=.False.,verbose=.false.,& 2625 & elec_eval=need_elec_eval) 2626 2627 call effective_potential_evaluate(eff_pot,energy,fcart,gred,strten,natom,rprimd,& 2628 & xred=xred,compute_anharmonic=need_anharmonic,verbose=.false.,& 2629 & filename=file_anh,elec_eval=need_elec_eval) 2630 2631 if(need_print .and. iam_master)then 2632 WRITE(unit_energy ,'(I10,7(F23.14))') ii,hist%etot(ii),energy_harm,energy,& 2633 & abs(hist%etot(ii) - energy_harm),abs(hist%etot(ii) - energy) 2634 WRITE(unit_stress,'(I10,12(F23.14))') ii,hist%strten(:,ii),strten(:) 2635 end if 2636 2637 !MS Uncommented for abihist test 2638 if(need_prt_ph)then 2639 if(ii == 1)then 2640 write(msg,'(a,(80a))') ch10,('-',i=1,80) 2641 call wrtout(ab_out,msg,'COLL') 2642 call wrtout(std_out,msg,'COLL') 2643 write(msg,'(3a)') ch10,'test_prt_ph == 1, write evulation of Model on the TEST-set into ph_test.nc',ch10 2644 call wrtout(ab_out,msg,'COLL') 2645 call wrtout(std_out,msg,'COLL') 2646 endif 2647 ifirst=merge(0,1,(ii>1)) 2648 filename_hist = trim("ph_test.nc") 2649 hist_out%fcart(:,:,hist_out%ihist) = fcart(:,:) 2650 hist_out%strten(:,hist_out%ihist) = strten(:) 2651 hist_out%etot(hist_out%ihist) = energy 2652 hist_out%entropy(hist_out%ihist) = hist%entropy(ii) 2653 hist_out%time(hist_out%ihist) = real(ii,kind=dp) 2654 ! call vel2hist(ab_mover%amass,hist,vel,vel_cell) 2655 call var2hist(hist%acell(:,ii),hist_out,natom,hist%rprimd(:,:,ii),hist%xred(:,:,ii),.false.) 2656 if(iam_master)then 2657 call write_md_hist(hist_out,filename_hist,ifirst,ii,natom,1,eff_pot%crystal%ntypat,& 2658 & eff_pot%supercell%typat,eff_pot%crystal%amu,eff_pot%crystal%znucl,& 2659 & real(100,dp),(/real(100,dp),real(100,dp)/)) 2660 endif 2661 endif!(need_prt_ph) 2662 2663 mse = mse + ((hist%etot(ii) - energy))**2/(sqomega(ii)**(1.0/2.0)) !+abs(hist$etot(ii) - energy) 2664 do ia=1,natom ! Loop over atoms 2665 do mu=1,3 ! Loop over cartesian directions 2666 msef = msef + (hist%fcart(mu,ia,ii) - fcart(mu,ia))**2 2667 end do 2668 end do 2669 do mu=1,6 ! Loop over stresses 2670 mses = mses + sqomega(ii)*(hist%strten(mu,ii) - strten(mu))**2 2671 end do 2672 end do ! End loop itime 2673 if(need_prt_ph)then 2674 write(msg,'(a,(80a))') ch10,('-',i=1,80) 2675 call wrtout(ab_out,msg,'COLL') 2676 call wrtout(std_out,msg,'COLL') 2677 endif 2678 2679 mse = mse / ntime 2680 msef = msef / (3*natom*ntime) 2681 mses = mses / (6*ntime) 2682 2683 if(need_print)then 2684 close(unit_energy) 2685 close(unit_stress) 2686 end if 2687 2688 !MS uncommented for PHONOPY TEST 2689 call abihist_free(hist_out) 2690 2691 end subroutine fit_polynomial_coeff_computeMSD
m_fit_polynomial_coeff/fit_polynomial_coeff_fit [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_fit
FUNCTION
Fit the list of coefficients included in eff_pot, if the coefficients are not set in eff_pot, this routine will genenerate a list of coefficients by taking into acount the symmetries of the system and the cutoff
INPUTS
eff_pot<type(effective_potential)> = effective potential
bancoeff(nbancoeff) = list of bannned coeffcients, these coefficients will NOT be
used during the fit process
fixcoeff(nfixcoeff) = list of fixed coefficient, these coefficients will be
imposed during the fit process
hist<type(abihist)> = The history of the MD (or snapshot of DFT)
generateterm = term to activate the generation of the term set
power_disps(2) = array with the minimal and maximal power_disp to be computed
nbancoeff = number of banned coeffcients
ncycle_in = number of maximum cycle (maximum coefficient to be fitted)
nfixcoeff = Number of coefficients imposed during the fit process
option = option of the fit process : 1 - selection of the coefficient one by one
2 - selection of the coefficients with Monte Carlo(testversion)
comm = MPI communicator
cutoff_in = optional,cut off to apply to the range of interation if
the coefficient are genereted in this routine
max_power_strain = maximum order of the strain of the strain phonon coupling
fit_initializeData = optional, logical !If true, we store all the information for the fit,
it will reduce the computation time but increase a lot the memory...
fit_tolMSDF = optional, tolerance in eV^2/A^2 on the Forces for the fit process
fit_tolMSDS = optional, tolerance in eV^2/A^2 on the Stresses for the fit process
fit_tolMSDE = optional, tolerance in meV^2/A^2 on the Energy for the fit process
fit_tolMSDFS= optional, tolerance in eV^2/A^2 on the Forces+stresses for the fit process
positive = optional, TRUE will return only positive coefficients
FALSE, default
verbose = optional, flag for the verbose mode
anhstr = logical, optional : TRUE, the anharmonic strain are computed
FALSE, (default) the anharmonic strain are not computed
only_odd_power = logical, optional : if TRUE generate only odd power
only_even_power= logical, optional : if TRUE generate only even power
OUTPUT
eff_pot<type(effective_potential)> = effective potential datatype with new fitted coefficients
SOURCE
119 subroutine fit_polynomial_coeff_fit(eff_pot,bancoeff,fixcoeff,hist,generateterm,power_disps,& 120 & nbancoeff,ncycle_in,nfixcoeff,nimposecoeff,imposecoeff,& 121 & option,comm,cutoff_in,max_power_strain,initialize_data,& 122 & fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS,fit_tolGF,& 123 & positive,verbose,anharmstr,spcoupling,& 124 & only_odd_power,only_even_power,prt_anh,& 125 & fit_iatom,prt_files,fit_on,sel_on,fit_factors,prt_GF_csv,& 126 & dispterms) 127 128 implicit none 129 130 !Arguments ------------------------------------ 131 !scalars 132 integer,intent(in) :: ncycle_in,nfixcoeff,nimposecoeff,comm 133 integer,intent(in) :: generateterm,nbancoeff,option 134 !arrays 135 integer,intent(in) :: fixcoeff(nfixcoeff), bancoeff(nbancoeff),imposecoeff(nimposecoeff) 136 integer,intent(in) :: power_disps(2) 137 type(effective_potential_type),target,intent(inout) :: eff_pot 138 type(abihist),intent(inout) :: hist 139 integer,optional,intent(in) :: max_power_strain,prt_anh,fit_iatom 140 real(dp),optional,intent(in) :: cutoff_in,fit_tolMSDF,fit_tolMSDS,fit_tolMSDE,fit_tolMSDFS 141 real(dp),optional,intent(in) :: fit_tolGF 142 logical,optional,intent(in) :: verbose,positive,anharmstr,spcoupling 143 logical,optional,intent(in) :: only_odd_power,only_even_power 144 logical,optional,intent(in) :: initialize_data,prt_files,prt_GF_csv 145 logical,optional,intent(in) :: fit_on(3), sel_on(3),dispterms 146 real(dp),optional,intent(in) :: fit_factors(3) 147 !Local variables------------------------------- 148 !scalar 149 integer :: ii,icoeff,my_icoeff,icycle,icycle_tmp,ierr,info,index_min,iproc,isweep,jcoeff,ia 150 integer :: master,max_power_strain_in,my_rank,my_ncoeff,ncoeff_model,ncoeff_tot,natom_sc,ncell,ncycle 151 integer :: ncycle_tot,ncycle_max,nproc,ntime,nsweep,size_mpi,ncoeff_fix,ncoeff_out 152 integer :: rank_to_send,unit_anh,fit_iatom_in,unit_GF_val,nfix_and_impose,nfixcoeff_corr 153 integer :: ncopy_terms 154 real(dp) :: cutoff,factor,time,tolMSDF,tolMSDS,tolMSDE,tolMSDFS,tolGF,check_value 155 real(dp),parameter :: HaBohr_eVAng = 27.21138386d0 / 0.529177249d0 156 type(effective_potential_type) :: eff_pot_fixed 157 logical :: iam_master,need_verbose,need_positive,converge,file_opened 158 logical :: need_anharmstr,need_spcoupling,ditributed_coefficients,need_prt_anh 159 logical :: need_only_odd_power,need_only_even_power,need_initialize_data 160 logical :: need_prt_files,need_prt_GF_csv,need_disp 161 !arrays 162 real(dp) :: mingf(4),int_fit_factors(3) 163 integer :: sc_size(3) 164 logical,allocatable :: fix_and_impose(:) 165 integer,allocatable :: buffsize(:),buffdisp(:),buffin(:),fixcoeff_corr(:) 166 integer,allocatable :: list_coeffs(:),list_coeffs_tmp(:),list_coeffs_tmp2(:) 167 integer,allocatable :: my_coeffindexes(:),singular_coeffs(:) 168 integer,allocatable :: my_coefflist(:) ,stat_coeff(:),list_coeffs_copy(:) 169 real(dp),allocatable :: gf_values_iter(:,:) 170 real(dp),allocatable :: buffGF(:,:),coeff_values(:),energy_coeffs(:,:) 171 real(dp),allocatable :: energy_coeffs_tmp(:,:) 172 real(dp),allocatable :: fcart_coeffs(:,:,:,:),gf_values(:,:),gf_mpi(:,:) 173 real(dp),allocatable :: fcart_coeffs_tmp(:,:,:,:),strten_coeffs_tmp(:,:,:) 174 real(dp),allocatable :: strten_coeffs(:,:,:) 175 type(polynomial_coeff_type),allocatable :: my_coeffs(:) 176 type(polynomial_coeff_type),allocatable :: coeffs_out(:) 177 type(polynomial_coeff_type),target,allocatable :: coeffs_tmp(:) 178 type(polynomial_coeff_type),pointer :: coeffs_in(:) 179 type(fit_data_type) :: fit_data 180 character(len=1000) :: message,message2 181 character(len=fnlen) :: filename 182 character(len=3) :: i_char 183 character(len=7) :: j_char 184 character(len=5),allocatable :: symbols(:) 185 ! ************************************************************************* 186 187 !MPI variables 188 master = 0 189 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 190 iam_master = (my_rank == master) 191 192 !Initialisation of optional arguments 193 need_verbose = .TRUE. 194 if(present(verbose)) need_verbose = verbose 195 need_initialize_data = .TRUE. 196 if(present(initialize_data)) need_initialize_data = initialize_data 197 need_positive = .FALSE. 198 if(present(positive)) need_positive = positive 199 need_anharmstr = .FALSE. 200 if(present(anharmstr)) need_anharmstr = anharmstr 201 need_disp = .TRUE. 202 if(present(dispterms)) need_disp = dispterms 203 need_spcoupling = .TRUE. 204 if(present(spcoupling)) need_spcoupling = spcoupling 205 need_only_odd_power = .FALSE. 206 if(present(only_odd_power)) need_only_odd_power = only_odd_power 207 need_prt_anh = .FALSE. 208 if(present(prt_anh))then 209 if(prt_anh == 1) need_prt_anh = .TRUE. 210 end if 211 need_prt_GF_csv = .FALSE. 212 if(present(prt_GF_csv)) need_prt_GF_csv = prt_GF_csv 213 need_prt_files = .TRUE. 214 if(present(prt_files))need_prt_files=prt_files 215 need_only_even_power = .FALSE. 216 if(present(only_even_power)) need_only_even_power = only_even_power 217 if(need_only_odd_power.and.need_only_even_power)then 218 write(message, '(3a)' )& 219 & 'need_only_odd_power and need_only_even_power are both true',ch10,& 220 & 'Action: contact abinit group' 221 ABI_ERROR(message) 222 end if 223 max_power_strain_in = 1 224 if(present(max_power_strain))then 225 max_power_strain_in = max_power_strain 226 end if 227 if(max_power_strain_in <= 0)then 228 write(message, '(3a)' )& 229 & 'max_power_strain can not be inferior or equal to zero',ch10,& 230 & 'Action: contact abinit group' 231 ABI_ERROR(message) 232 end if 233 !Check which atom to fit, if not present do all atoms 234 if(present(fit_iatom))then 235 fit_iatom_in = fit_iatom 236 else 237 fit_iatom_in = -1 238 endif 239 !Set int fit factors to default value if fit factors not present 240 int_fit_factors = (/1,1,1/) 241 if (present(fit_factors)) int_fit_factors = fit_factors 242 243 !Set the tolerance for the fit 244 tolMSDF=zero;tolMSDS=zero;tolMSDE=zero;tolMSDFS=zero;tolGF=zero 245 if(present(fit_tolMSDF)) tolMSDF = fit_tolMSDF 246 if(present(fit_tolMSDS)) tolMSDS = fit_tolMSDS 247 if(present(fit_tolMSDE)) tolMSDE = fit_tolMSDE 248 if(present(fit_tolMSDFS))tolMSDFS = fit_tolMSDFS 249 if(present(fit_tolGF)) tolGF = fit_tolGF 250 251 !Copy the input effective potential eff_pot to eff_pot fixed 252 !If nimposecoeff=0 the fixed potential is the harmonic potential 253 ncoeff_model = eff_pot%anharmonics_terms%ncoeff 254 if (nimposecoeff > ncoeff_model)then 255 write(message,'(2a)') "fit_nimposecoeff is greater then the number of anharmonic terms& 256 & provided by input effective potential."& 257 & ,"Action -> Change fit_nimposecoeff, and fit_imposecoeff in the input" 258 ABI_ERROR(message) 259 endif 260 !Impose some coefficients of the input potential 261 if (nimposecoeff > 0)then 262 if (any(imposecoeff > ncoeff_model))then 263 write(message,'(2a)') "At least one value in fit_imposeccoeff is greater then the number of anharmonic terms& 264 & provided by input effective potential."& 265 & ,"Action -> Change fit_imposecoeff in the input" 266 ABI_ERROR(message) 267 endif 268 ABI_MALLOC(coeffs_tmp,(nimposecoeff)) 269 do ia = 1,nimposecoeff 270 ii = imposecoeff(ia) 271 call polynomial_coeff_init(eff_pot%anharmonics_terms%coefficients(ii)%coefficient,& 272 & eff_pot%anharmonics_terms%coefficients(ii)%nterm,coeffs_tmp(ia),& 273 & eff_pot%anharmonics_terms%coefficients(ii)%terms,& 274 & eff_pot%anharmonics_terms%coefficients(ii)%name,& 275 & check = .TRUE.) 276 enddo 277 !Copy the input eff pot, free the coeffs and set the ones who shall be imposed to fixed 278 call effective_potential_copy(eff_pot_fixed,eff_pot,comm) 279 call effective_potential_freeCoeffs(eff_pot_fixed) 280 call effective_potential_setCoeffs(coeffs_tmp,eff_pot_fixed,nimposecoeff) 281 !Deallocate coeffs tmp 282 do ii = 1,nimposecoeff 283 call polynomial_coeff_free(coeffs_tmp(ii)) 284 enddo 285 ABI_FREE(coeffs_tmp) 286 !Fix the whole input potential, copy the input potential to eff_pot 287 elseif (nimposecoeff == -1)then 288 call effective_potential_copy(eff_pot_fixed,eff_pot,comm) 289 !Fix only the harmonic potential 290 else 291 call effective_potential_copy(eff_pot_fixed,eff_pot,comm) 292 call effective_potential_freeCoeffs(eff_pot_fixed) 293 endif 294 295 ncopy_terms = 0 296 !Set consistency between fixcoeff and imposecoeff. 297 if ( nfixcoeff > 0 .and. nimposecoeff >0)then 298 ABI_MALLOC(fix_and_impose,(nfixcoeff)) 299 fix_and_impose = .FALSE. 300 do ii = 1,nfixcoeff 301 if (any(imposecoeff == fixcoeff(ii)))then 302 fix_and_impose(ii) = .TRUE. 303 endif 304 enddo 305 nfix_and_impose = count(fix_and_impose) 306 nfixcoeff_corr = nfixcoeff - nfix_and_impose 307 ABI_MALLOC(fixcoeff_corr,(nfixcoeff_corr)) 308 ia = 1 309 do ii = 1,nfixcoeff 310 if (.not. fix_and_impose(ii))then 311 fixcoeff_corr(ia) = fixcoeff(ii) 312 ia = ia + 1 313 endif 314 enddo 315 ncopy_terms = ncoeff_model - nimposecoeff 316 ABI_MALLOC(list_coeffs_copy,(ncopy_terms)) 317 ia = 1 318 do ii = 1,ncoeff_model 319 if( .not. any(imposecoeff == ii))then 320 list_coeffs_copy(ia) = ii 321 ia = ia + 1 322 endif 323 enddo 324 do ii = 1,ncopy_terms 325 do ia = 1,nfixcoeff_corr 326 if (list_coeffs_copy(ii) == fixcoeff_corr(ia))then 327 fixcoeff_corr(ia) = ii 328 endif 329 enddo 330 enddo 331 elseif (nfixcoeff == -1 .and. nimposecoeff > 0)then 332 ABI_MALLOC(fix_and_impose,(ncoeff_model)) 333 fix_and_impose = .FALSE. 334 do ii = 1,ncoeff_model 335 if (any( imposecoeff == ii))then 336 fix_and_impose(ii) = .TRUE. 337 endif 338 enddo 339 nfix_and_impose = nimposecoeff 340 ncopy_terms = ncoeff_model - nimposecoeff 341 ABI_MALLOC(fixcoeff_corr,(ncopy_terms)) 342 ABI_MALLOC(list_coeffs_copy,(ncopy_terms)) 343 ia = 1 344 do ii = 1,ncoeff_model 345 if (.not. fix_and_impose(ii))then 346 fixcoeff_corr(ia) = ia 347 list_coeffs_copy(ia) = ii 348 ia = ia + 1 349 endif 350 enddo 351 nfixcoeff_corr = ncoeff_model- nfix_and_impose 352 elseif (nfixcoeff == -1 .and. nimposecoeff ==-1)then 353 nfixcoeff_corr = 0 354 write(message,'(3a)') "nfixcoeff and nimposecoeff are set to -1.",ch10,& 355 & "This does not make sense. nfixcoeff will be set to 0." 356 if(iam_master) ABI_WARNING(message) 357 ncopy_terms = 0 358 ABI_MALLOC(list_coeffs_copy,(ncopy_terms)) 359 elseif (nfixcoeff >0 .and. nimposecoeff ==-1)then 360 nfixcoeff_corr = 0 361 write(message,'(3a)') "nfixcoeff is > 0 and nimposecoeff is set to -1.",ch10,& 362 & "This does not make sense. nfixcoeff will be set to 0." 363 if(iam_master) ABI_WARNING(message) 364 ncopy_terms = 0 365 else 366 nfixcoeff_corr = nfixcoeff 367 ABI_MALLOC(fixcoeff_corr,(nfixcoeff_corr)) 368 !nimposecoeff or nfixcoeff always 0 here so fix_and_impose is empty 369 ABI_MALLOC(fix_and_impose,(0)) 370 fixcoeff_corr = fixcoeff 371 ncopy_terms = ncoeff_model 372 ABI_MALLOC(list_coeffs_copy,(ncopy_terms)) 373 do ii = 1,ncopy_terms 374 list_coeffs_copy(ii) = ii 375 enddo 376 endif 377 378 if(need_verbose) then 379 write(message,'(a,(80a))') ch10,('-',ii=1,80) 380 call wrtout(ab_out,message,'COLL') 381 call wrtout(std_out,message,'COLL') 382 write(message,'(2a)') ch10,' Starting Fit process' 383 call wrtout(ab_out,message,'COLL') 384 call wrtout(std_out,message,'COLL') 385 write(message,'(a,(80a))') ch10,('-',ii=1,80) 386 call wrtout(ab_out,message,'COLL') 387 call wrtout(std_out,message,'COLL') 388 end if 389 390 ditributed_coefficients = .true. 391 if(option==2) ditributed_coefficients = .false. 392 393 !if the number of atoms in reference supercell into effpot is not correct, 394 !wrt to the number of atom in the hist, we set map the hist and set the good supercell 395 if (size(hist%xred,2) /= eff_pot%supercell%natom) then 396 call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose) 397 end if 398 399 !Set the cut off 400 cutoff = zero 401 if(present(cutoff_in))then 402 cutoff = cutoff_in 403 end if 404 !If the cutoff is set to zero, we define a default value 405 if(abs(cutoff)<tol16)then 406 do ii=1,3 407 cutoff = cutoff + sqrt(eff_pot%supercell%rprimd(ii,1)**2+& 408 & eff_pot%supercell%rprimd(ii,2)**2+& 409 & eff_pot%supercell%rprimd(ii,3)**2) 410 end do 411 cutoff = cutoff / 3.0_dp 412 end if 413 !we get the size of the supercell in the hist file 414 do ii=1,3 415 sc_size(ii) = int(anint(sqrt(eff_pot%supercell%rprimd(ii,1)**2+& 416 & eff_pot%supercell%rprimd(ii,2)**2+& 417 & eff_pot%supercell%rprimd(ii,3)**2) / & 418 & sqrt(eff_pot%crystal%rprimd(ii,1)**2+& 419 & eff_pot%crystal%rprimd(ii,2)**2+& 420 & eff_pot%crystal%rprimd(ii,3)**2))) 421 end do 422 423 424 !Get the list of coefficients to fit: 425 !get from the eff_pot type (from the input) 426 !or 427 !regenerate the list 428 my_ncoeff = 0 429 ncoeff_tot = 0 430 431 !Reset ncoeff_tot 432 if(ncoeff_model > 0)then 433 if(need_verbose)then 434 write(message, '(4a)' )ch10,' The coefficients present in the effective',& 435 & ' potential will be used for the fit' 436 call wrtout(std_out,message,'COLL') 437 call wrtout(ab_out,message,'COLL') 438 end if 439 end if 440 441 ABI_MALLOC(symbols,(eff_pot%crystal%natom)) 442 call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,& 443 & symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl) 444 445 if(generateterm == 1)then 446 ! we need to regenerate them 447 if(need_verbose)then 448 if(fit_iatom_in > 0)then 449 write(message, '(2a,I3,4a)' )ch10,' The coefficients for the fit around atom', fit_iatom_in,': ',& 450 & trim(symbols(fit_iatom)),', will be generated',ch10 451 call wrtout(std_out,message,'COLL') 452 call wrtout(ab_out,message,'COLL') 453 else 454 write(message, '(3a)' )ch10,' The coefficients for the fit will be generated with all cross terms',ch10 455 call wrtout(std_out,message,'COLL') 456 call wrtout(ab_out,message,'COLL') 457 endif 458 write(message,'(a,F6.3,a)') " Cutoff of ",cutoff," Bohr is imposed" 459 call wrtout(std_out,message,'COLL') 460 end if 461 462 call polynomial_coeff_getNorder(coeffs_tmp,eff_pot%crystal,cutoff,my_ncoeff,ncoeff_tot,power_disps,& 463 & max_power_strain_in,0,sc_size,comm,anharmstr=need_anharmstr,& 464 & spcoupling=need_spcoupling,distributed=.true.,& 465 & only_odd_power=need_only_odd_power,& 466 & only_even_power=need_only_even_power,& 467 & fit_iatom=fit_iatom_in,dispterms=need_disp) 468 end if 469 470 ABI_FREE(symbols) 471 472 !Copy the initial coefficients from the model on the CPU 0 473 ncoeff_tot = ncoeff_tot + ncoeff_model 474 if((iam_master .and. ncopy_terms > 0)) my_ncoeff = my_ncoeff + ncopy_terms 475 476 !Get number of fixed coeff 477 ncoeff_fix = 0 478 if(nfixcoeff_corr /=0) then 479 if(nfixcoeff_corr == -1)then 480 ncoeff_fix = ncoeff_model 481 else 482 ncoeff_fix = nfixcoeff_corr 483 endif 484 endif 485 486 487 !Get the list with the number of coeff on each CPU 488 !In order to be abble to compute the my_coeffindexes array which is for example: 489 ! if CPU0 has 200 Coeff and CPU1 has 203 Coeff then 490 ! for CPU0:my_coeffindexes=>1-200 and for CPU1:my_coeffindexes=>201-403 491 !Also fill the my_coeffs array with the generated coefficients and/or the coefficient from the input xml 492 ABI_MALLOC(buffin,(nproc)) 493 buffin = 0 494 buffin(my_rank+1) = my_ncoeff 495 call xmpi_sum(buffin,comm,ierr) 496 ABI_MALLOC(my_coeffindexes,(my_ncoeff)) 497 ABI_MALLOC(my_coefflist,(my_ncoeff)) 498 ABI_MALLOC(my_coeffs,(my_ncoeff)) 499 do icoeff=1,my_ncoeff 500 501 jcoeff = icoeff 502 my_coefflist(icoeff) = icoeff 503 504 if(my_rank==0) then 505 my_coeffindexes(icoeff) = icoeff 506 else 507 my_coeffindexes(icoeff) = sum(buffin(1:my_rank)) + icoeff 508 end if 509 510 ! Only copy the input coefficients on the CPU0 511 if(my_rank==0) then 512 if(icoeff <= ncopy_terms)then 513 coeffs_in => eff_pot%anharmonics_terms%coefficients 514 jcoeff = list_coeffs_copy(icoeff) 515 else 516 coeffs_in => coeffs_tmp 517 jcoeff = jcoeff - ncopy_terms 518 end if 519 else 520 coeffs_in => coeffs_tmp 521 end if 522 call polynomial_coeff_init(one,coeffs_in(jcoeff)%nterm,& 523 & my_coeffs(icoeff),coeffs_in(jcoeff)%terms,& 524 & coeffs_in(jcoeff)%name,& 525 & check=.true.) 526 call polynomial_coeff_free(coeffs_in(jcoeff)) 527 end do 528 529 !Deallocation 530 ABI_SFREE(coeffs_tmp) 531 NULLIFY(coeffs_in) 532 ABI_FREE(buffin) 533 534 !wait everybody 535 call xmpi_barrier(comm) 536 537 !Reset the output (we free the memory) 538 call effective_potential_freeCoeffs(eff_pot) 539 540 541 !Check if ncycle_in is not zero or superior to ncoeff_tot 542 if(need_verbose.and.(ncycle_in > ncoeff_tot).or.(ncycle_in<0.and.nfixcoeff_corr /= -1)) then 543 write(message, '(6a,I0,3a)' )ch10,& 544 & ' --- !WARNING',ch10,& 545 & ' The number of cycle requested in the input is not correct.',ch10,& 546 & ' This number will be set to the maximum of coefficients: ',ncoeff_tot,ch10,& 547 & ' ---',ch10 548 call wrtout(std_out,message,"COLL") 549 end if 550 551 !Use fixcoeff 552 !ncycle_tot store the curent number of coefficient in the model 553 !Do not reset this variable... 554 ncycle_tot = 0 555 if (nfixcoeff_corr == -1)then 556 write(message, '(3a)')' nfixcoeff is set to -1, the coefficients present in the model',& 557 & ' are imposed.',ch10 558 ncycle_tot = ncycle_tot + ncoeff_model 559 else 560 if (nfixcoeff_corr > 0)then 561 if(maxval(fixcoeff_corr(:)) > ncoeff_tot) then 562 write(message, '(4a,I0,6a)' )ch10,& 563 & ' --- !WARNING',ch10,& 564 & ' The value ',maxval(fixcoeff_corr(:)),' is not in the list.',ch10,& 565 & ' Start from scratch...',ch10,& 566 & ' ---',ch10 567 else 568 ncycle_tot = ncycle_tot + nfixcoeff_corr 569 write(message, '(2a)')' Some coefficients are imposed from the input.',ch10 570 end if 571 else 572 write(message, '(4a)')' There is no coefficient imposed from the input.',ch10,& 573 & ' Start from scratch',ch10 574 end if 575 end if 576 577 if(need_verbose) call wrtout(std_out,message,'COLL') 578 579 !Compute the number of cycle: 580 ncycle = ncycle_in 581 !Compute the maximum number of cycle 582 ncycle_max = ncycle_in + ncycle_tot 583 584 !Check if the number of request cycle + the initial number of coeff is superior to 585 !the maximum number of coefficient allowed 586 if(ncycle_max > ncoeff_tot) then 587 ncycle = ncoeff_tot - ncycle_tot 588 ncycle_max = ncoeff_tot 589 write(message, '(4a,I0,2a,I0,2a,I0,3a)' )ch10,& 590 & ' --- !WARNING',ch10,& 591 & ' The number of cycle + the number of imposed coefficients: ',ncycle_max,ch10,& 592 & ' is superior to the maximum number of coefficients in the initial list: ',ncoeff_tot,ch10,& 593 & ' The number of cycle is set to ',ncycle,ch10,& 594 & ' ---',ch10 595 if(need_verbose) call wrtout(std_out,message,'COLL') 596 else if (option==2)then 597 ! Always set to the maximum 598 ncycle_max = ncoeff_tot 599 end if 600 601 !Initialisation of constants 602 ntime = hist%mxhist 603 natom_sc = eff_pot%supercell%natom 604 ncell = eff_pot%supercell%ncells 605 factor = 1._dp/natom_sc 606 607 !Initialisation of arrays: 608 ABI_MALLOC(energy_coeffs_tmp,(ncycle_max,ntime)) 609 ABI_MALLOC(list_coeffs,(ncycle_max)) 610 ABI_MALLOC(fcart_coeffs_tmp,(3,natom_sc,ncycle_max,ntime)) 611 ABI_MALLOC(strten_coeffs_tmp,(6,ntime,ncycle_max)) 612 list_coeffs = 0 613 614 !if ncycle_tot > 0 fill list_coeffs with the fixed coefficients 615 if(ncycle_tot > 0)then 616 do ii = 1,ncycle_tot 617 if(nfixcoeff_corr == -1)then 618 if(ii <= ncoeff_model)then 619 list_coeffs(ii) = ii 620 end if 621 else 622 list_coeffs(ii) = fixcoeff_corr(ii) 623 end if 624 end do 625 end if 626 627 !Get the decomposition for each coefficients of the forces and stresses for 628 !each atoms and each step equations 11 & 12 of PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]] 629 if(need_verbose)then 630 write(message, '(a)' ) ' Initialisation of the fit process...' 631 call wrtout(std_out,message,'COLL') 632 end if 633 !Before the fit, compute constants with fit_data_compute. 634 !Conpute the strain of each configuration. 635 !Compute the displacmeent of each configuration. 636 !Compute the variation of the displacement due to strain of each configuration. 637 !Compute fixed forces and stresse and get the standard deviation. 638 !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]. 639 call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose) 640 641 !Get the decomposition for each coefficients of the forces,stresses and energy for 642 !each atoms and each step (see equations 11 & 12 of 643 !PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]]+ allocation 644 !If the user does not turn off this initialization, we store all the information for the fit, 645 !it will reduce the computation time but increase a lot the memory... 646 if(need_initialize_data)then 647 ABI_MALLOC(energy_coeffs,(my_ncoeff,ntime)) 648 ABI_MALLOC(fcart_coeffs,(3,natom_sc,my_ncoeff,ntime)) 649 ABI_MALLOC(strten_coeffs,(6,ntime,my_ncoeff)) 650 call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,& 651 & fit_data%training_set%displacement,& 652 & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,& 653 & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,& 654 & strten_coeffs,fit_data%training_set%ucvol,my_coefflist,my_ncoeff) 655 else 656 ! Allocate just 1 dimension ! Save MEMORY ! 657 ABI_MALLOC(energy_coeffs,(1,ntime)) 658 ABI_MALLOC(fcart_coeffs,(3,natom_sc,1,ntime)) 659 ABI_MALLOC(strten_coeffs,(6,ntime,1)) 660 end if 661 662 !Allocation of arrays 663 ABI_MALLOC(coeffs_tmp,(ncycle_max)) 664 ABI_MALLOC(singular_coeffs,(max(1,my_ncoeff))) 665 ABI_MALLOC(coeff_values,(ncycle_max)) 666 ABI_MALLOC(gf_values,(4,max(1,my_ncoeff))) 667 ABI_MALLOC(list_coeffs_tmp,(ncycle_max)) 668 ABI_MALLOC(list_coeffs_tmp2,(ncycle_max)) 669 ABI_MALLOC(stat_coeff,(ncoeff_tot)) 670 coeff_values = zero 671 singular_coeffs = 0 672 stat_coeff = 0 673 !Set mpi buffer 674 !Set the bufsize for mpi allgather 675 ABI_MALLOC(buffsize,(nproc)) 676 ABI_MALLOC(buffdisp,(nproc)) 677 ABI_MALLOC(buffGF,(5,1)) 678 ABI_MALLOC(gf_mpi,(5,nproc)) 679 buffsize(:) = 0 680 buffdisp(1) = 0 681 do ii= 1,nproc 682 buffsize(ii) = 5 683 end do 684 do ii = 2,nproc 685 buffdisp(ii) = buffdisp(ii-1) + buffsize(ii-1) 686 end do 687 size_mpi = 5*nproc 688 !If some coeff are imposed by the input, we need to fill the arrays 689 !with this coeffs and broadcast to the others CPUs : 690 if(ncycle_tot>=1)then 691 do icycle = 1,ncycle_tot 692 list_coeffs_tmp(icycle) = icycle 693 rank_to_send = 0 694 do icoeff=1,my_ncoeff 695 if((my_coeffindexes(icoeff)==list_coeffs(icycle)))then 696 697 if(need_initialize_data)then 698 my_icoeff = icoeff 699 else 700 my_icoeff = 1 701 ! Need to initialized the data for the fit for this coefficient 702 call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,& 703 & fit_data%training_set%displacement,& 704 & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,& 705 & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,& 706 & strten_coeffs,fit_data%training_set%ucvol,& 707 & my_coefflist(icoeff),1) 708 end if 709 710 energy_coeffs_tmp(icycle,:) = energy_coeffs(my_icoeff,:) 711 fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:) 712 strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,my_icoeff) 713 rank_to_send = my_rank 714 call polynomial_coeff_free(coeffs_tmp(icycle)) 715 call polynomial_coeff_init(coeff_values(icycle),my_coeffs(icoeff)%nterm,& 716 & coeffs_tmp(icycle),my_coeffs(icoeff)%terms,& 717 & my_coeffs(icoeff)%name,& 718 & check=.false.) 719 exit 720 end if 721 end do 722 ! Need to send the rank with the chosen coefficient 723 call xmpi_sum(rank_to_send, comm, ierr) 724 ! Boadcast the coefficient 725 call xmpi_bcast(energy_coeffs_tmp(icycle,:), rank_to_send, comm, ierr) 726 call xmpi_bcast(fcart_coeffs_tmp(:,:,icycle,:) , rank_to_send, comm, ierr) 727 call xmpi_bcast(strten_coeffs_tmp(:,:,icycle), rank_to_send, comm, ierr) 728 call polynomial_coeff_broadcast(coeffs_tmp(icycle), rank_to_send, comm) 729 end do 730 end if 731 732 !Waiting for all 733 if(nproc > 1) then 734 if(need_verbose)then 735 write(message, '(a)') ' Initialisation done... waiting for all the CPU' 736 call wrtout(std_out,message,'COLL') 737 end if 738 call xmpi_barrier(comm) 739 end if 740 741 !Compute GF, coeff_values,strten_coeffs and fcart_coeffs are set to zero 742 !it means that only the harmonic part wiil be computed 743 coeff_values = zero 744 745 call fit_polynomial_coeff_computeGF(coeff_values,energy_coeffs,fit_data%energy_diff,fcart_coeffs,& 746 & fit_data%fcart_diff,gf_values(:,1),int((/1/)),natom_sc,& 747 & 0,my_ncoeff,ntime,strten_coeffs,fit_data%strten_diff,& 748 & fit_data%training_set%sqomega) 749 750 !Print the standard deviation before the fit 751 write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' ) ch10,& 752 !& ' Mean Standard Deviation values at the begining of the fit process (meV**2/atm):',& 753 !& ch10,' Energy : ',& 754 !& gf_values(4,1)*factor*(Ha_EV*1000)**2 ,ch10,& 755 & ' Goal function values at the begining of the fit process (eV^2/A^2):',ch10,& 756 & ' Energy : ',& 757 & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,& 758 & ' Forces+Stresses : ',& 759 & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,& 760 & ' Forces : ',& 761 & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,& 762 & ' Stresses : ',& 763 & gf_values(3,1)*(HaBohr_eVAng)**2,ch10 764 if(need_verbose)then 765 call wrtout(ab_out,message,'COLL') 766 call wrtout(std_out,message,'COLL') 767 end if 768 769 ABI_MALLOC(gf_values_iter,(4,ncycle+1)) 770 gf_values_iter(:,:) = zero 771 !Store initial gf_values as first value in gf_values_iter 772 gf_values_iter(:,1) = gf_values(:,1) 773 774 select case(option) 775 776 case(1) 777 !Option 1, we select the coefficients one by one 778 if(need_verbose.and.ncycle > 0)then 779 write(message,'(a,3x,a,10x,a,14x,a,14x,a,14x,a)') " N","Selecting","MSDE","MSDFS","MSDF","MSDS" 780 call wrtout(ab_out,message,'COLL') 781 write(message,'(4x,a,6x,a,8x,a,8x,a,8x,a)') "Coefficient","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",& 782 & "(eV^2/A^2)" 783 call wrtout(ab_out,message,'COLL') 784 end if 785 786 ! Start fit process 787 do icycle_tmp = 1,ncycle 788 icycle = ncycle_tot + 1 789 list_coeffs_tmp(icycle) = icycle 790 if(need_verbose)then 791 write(message, '(4a,I0,a)')ch10,'--',ch10,' Try to find the best model with ',& 792 & icycle,' coefficient' 793 if(icycle > 1) write(message, '(2a)') trim(message),'s' 794 if(nproc > 1) then 795 if(my_ncoeff>=1) then 796 write(message, '(2a,I0,a)')trim(message), ' (only the ',my_ncoeff,& 797 & ' first are printed for this CPU)' 798 else 799 write(message, '(2a)')trim(message), ' (no coefficient treated by this CPU)' 800 end if 801 end if 802 call wrtout(std_out,message,'COLL') 803 if(icycle>1 .or. any(list_coeffs(:) > zero))then 804 write(message, '(3a)') ' The coefficient numbers from the previous cycle are:',ch10,' [' 805 do ii=1,icycle-1 806 if(ii<icycle-1)then 807 write(message, '(a,I0,a)') trim(message),list_coeffs(ii),',' 808 else 809 write(message, '(a,I0)') trim(message),list_coeffs(ii) 810 end if 811 end do 812 write(message, '(3a)') trim(message),']',ch10 813 call wrtout(std_out,message,'COLL') 814 end if 815 816 write(message,'(2x,a,12x,a,14x,a,13x,a,14x,a)') " Testing","MSDE","MSDFS","MSDF","MSDS" 817 call wrtout(std_out,message,'COLL') 818 write(message,'(a,7x,a,8x,a,8x,a,8x,a)') " Coefficient","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",& 819 & "(eV^2/A^2)" 820 call wrtout(std_out,message,'COLL') 821 end if!End if verbose 822 823 !Print all GF VALUES in CSV if wanted 824 !Open *csv file for storing GF values of all cores for this iteration 825 if(need_prt_GF_csv)then 826 write(filename,'(a,I1,a,I3.3,a,I3.3,a)') "GF_values_iatom",fit_iatom,"_proc",my_rank,"_iter",icycle,".csv" 827 unit_GF_val = get_unit() 828 if (open_file(filename,message,unit=unit_GF_val,form="formatted",& 829 & status="unknown",action="write") /= 0) then 830 ABI_ERROR(message) 831 end if 832 end if 833 ! Reset gf_values 834 gf_values(:,:) = zero 835 do icoeff=1,my_ncoeff 836 ! cycle if this coefficient is not allowed 837 if(any(list_coeffs==my_coeffindexes(icoeff)).or.singular_coeffs(icoeff) == 1)then 838 gf_values(:,icoeff) = zero 839 cycle 840 endif 841 if(nbancoeff >= 1)then 842 if(any(bancoeff==my_coeffindexes(icoeff)))then 843 gf_values(:,icoeff) = zero 844 cycle 845 endif 846 end if 847 list_coeffs(icycle) = my_coeffindexes(icoeff) 848 849 if(need_initialize_data)then 850 my_icoeff = icoeff 851 else 852 ! Need to initialized the data for the fit for this coefficient 853 my_icoeff = 1 854 call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,& 855 & fit_data%training_set%displacement,& 856 & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,& 857 & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,& 858 & strten_coeffs,fit_data%training_set%ucvol,& 859 & my_coefflist(icoeff),1) 860 end if 861 862 ! Fill the temporary arrays 863 energy_coeffs_tmp(icycle,:) = energy_coeffs(my_icoeff,:) 864 fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:) 865 strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,my_icoeff) 866 867 ! call the fit process routine 868 ! This routine solves the linear system proposed 869 ! by C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]] 870 call fit_polynomial_coeff_solve(coeff_values(1:icycle),fcart_coeffs_tmp,fit_data%fcart_diff,& 871 & energy_coeffs_tmp,fit_data%energy_diff,info,& 872 & list_coeffs_tmp(1:icycle),natom_sc,icycle,ncycle_max,ntime,& 873 & strten_coeffs_tmp,fit_data%strten_diff,& 874 & fit_data%training_set%sqomega,fit_on,int_fit_factors) 875 876 if(info==0)then 877 if (need_positive.and.any(coeff_values(ncoeff_fix+1:icycle) < zero)) then 878 write(message, '(a)') ' Negative value detected...' 879 gf_values(:,icoeff) = zero 880 coeff_values = zero 881 else 882 call fit_polynomial_coeff_computeGF(coeff_values(1:icycle),energy_coeffs_tmp,& 883 & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,& 884 & gf_values(:,icoeff),list_coeffs_tmp(1:icycle),natom_sc,& 885 & icycle,ncycle_max,ntime,strten_coeffs_tmp,& 886 & fit_data%strten_diff,fit_data%training_set%sqomega) 887 888 889 write (j_char, '(i7)') my_coeffindexes(icoeff) 890 write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char),& 891 !& gf_values(4,icoeff)*factor*(1000*Ha_ev)**2 ,& 892 & gf_values(4,icoeff)*HaBohr_eVAng**2,& 893 & gf_values(1,icoeff)*HaBohr_eVAng**2,& 894 & gf_values(2,icoeff)*HaBohr_eVAng**2,& 895 & gf_values(3,icoeff)*HaBohr_eVAng**2 896 if(need_prt_GF_csv)then 897 write(message2, '(I7.7,3a,ES18.10,a,ES18.10,a,ES18.10,a,ES18.10)') my_coeffindexes(icoeff),",",& 898 !& gf_values(4,icoeff)*factor*(1000*Ha_ev)**2,",",& 899 & trim(my_coeffs(icoeff)%name),",",& 900 & gf_values(4,icoeff)*HaBohr_eVAng**2,",",& 901 & gf_values(1,icoeff)*HaBohr_eVAng**2,",",& 902 & gf_values(2,icoeff)*HaBohr_eVAng**2,",",& 903 & gf_values(3,icoeff)*HaBohr_eVAng**2 904 end if 905 end if 906 else!In this case the matrix is singular 907 gf_values(:,icoeff) = zero 908 singular_coeffs(icoeff) = 1 909 write(message, '(a)') ' The matrix is singular...' 910 if(need_prt_GF_csv)then 911 write(message2, '(I7.7,10a)') my_coeffindexes(icoeff),",",& 912 !& gf_values(4,icoeff)*factor*(1000*Ha_ev)**2,",",& 913 & trim(my_coeffs(icoeff)%name),",",& 914 & "None",",",& 915 & "None",",",& 916 & "None",",",& 917 & "None" 918 endif 919 end if 920 if(need_verbose)then 921 call wrtout(std_out,message,'COLL') 922 if(need_prt_GF_csv)then 923 call wrtout(unit_GF_val,message2,'PERS',do_flush=.TRUE.) 924 end if 925 endif 926 end do !icoeff=1,my_ncoeff 927 928 !Close *csv file for GF values of this iteration 929 if(need_prt_GF_csv)close(unit_GF_val) 930 931 ! find the best coeff on each CPU 932 mingf(:) = 9D99 933 index_min = 0 934 do icoeff=1,my_ncoeff 935 if(gf_values(1,icoeff) < zero) cycle 936 if(abs(gf_values(1,icoeff)) <tol16) cycle 937 if(sum(gf_values(2:4,icoeff),MASK=sel_on) < sum(mingf(2:4),MASK=sel_on))then 938 mingf(:) = gf_values(:,icoeff) 939 index_min = my_coeffindexes(icoeff) 940 end if 941 end do 942 943 ! MPI GATHER THE BEST COEFF ON EACH CPU 944 if(nproc > 1)then 945 buffGF(1,1) = index_min 946 buffGF(2:5,1) = mingf(:) 947 948 call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr) 949 ! find the best coeff 950 mingf(:) = 9D99 951 index_min= 0 952 do icoeff=1,nproc 953 if(gf_mpi(2,icoeff) < zero) cycle 954 if(abs(gf_mpi(2,icoeff)) < tol16) cycle 955 if(sum(gf_mpi(3:5,icoeff),MASK=sel_on) < sum(mingf(2:4),MASK=sel_on))then 956 mingf(:) = gf_mpi(2:5,icoeff) 957 index_min = int(gf_mpi(1,icoeff)) 958 end if 959 end do 960 end if 961 962 ! Check if there is still coefficient 963 if(index_min==0) then 964 exit 965 else 966 list_coeffs(icycle) = index_min 967 end if 968 969 ! Check if this coeff is treat by this cpu and fill the 970 ! temporary array before broadcast 971 rank_to_send = 0 972 do icoeff=1,my_ncoeff 973 974 975 if((my_coeffindexes(icoeff)==list_coeffs(icycle)))then 976 977 if(need_initialize_data)then 978 my_icoeff = icoeff 979 else 980 ! Need to initialized the data for the fit for this coefficient 981 my_icoeff = 1 982 call fit_polynomial_coeff_getFS(my_coeffs,fit_data%training_set%du_delta,& 983 & fit_data%training_set%displacement,& 984 & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,& 985 & my_ncoeff,ntime,sc_size,fit_data%training_set%strain,& 986 & strten_coeffs,fit_data%training_set%ucvol,& 987 & my_coefflist(icoeff),1) 988 end if 989 990 energy_coeffs_tmp(icycle,:) = energy_coeffs(my_icoeff,:) 991 fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,my_icoeff,:) 992 strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,my_icoeff) 993 call polynomial_coeff_free(coeffs_tmp(icycle)) 994 call polynomial_coeff_init(coeff_values(icycle),my_coeffs(icoeff)%nterm,& 995 & coeffs_tmp(icycle),my_coeffs(icoeff)%terms,& 996 & my_coeffs(icoeff)%name,& 997 & check=.false.) 998 999 rank_to_send = my_rank 1000 exit 1001 end if 1002 end do 1003 ! Need to send the rank with the chosen coefficient 1004 call xmpi_sum(rank_to_send, comm, ierr) 1005 1006 ! Boadcast the coefficient 1007 call xmpi_bcast(energy_coeffs_tmp(icycle,:), rank_to_send, comm, ierr) 1008 call xmpi_bcast(fcart_coeffs_tmp(:,:,icycle,:) , rank_to_send, comm, ierr) 1009 call xmpi_bcast(strten_coeffs_tmp(:,:,icycle), rank_to_send, comm, ierr) 1010 call polynomial_coeff_broadcast(coeffs_tmp(icycle), rank_to_send, comm) 1011 1012 if(need_verbose) then 1013 write(message, '(a,I0,2a)' )' Selecting the coefficient number ',list_coeffs(icycle),& 1014 & ' ===> ',trim(coeffs_tmp(icycle)%name) 1015 call wrtout(std_out,message,'COLL') 1016 1017 ! write(message, '(2a,I0,a,ES24.16)' )' Standard deviation of the energy for',& 1018 !& ' the iteration ',icycle_tmp,' (meV^2/atm): ',& 1019 !& mingf(4)* factor * (Ha_eV *1000)**2 1020 ! call wrtout(std_out,message,'COLL') 1021 1022 write (i_char, '(i3)') icycle 1023 write (j_char, '(i7)') list_coeffs(icycle) 1024 write(message, '(a,a,3x,a,3x,4ES18.10)') " ",adjustl(i_char),adjustl(j_char),& 1025 !& mingf(4)* factor * (Ha_eV *1000)**2,& 1026 & mingf(4)*HaBohr_eVAng**2,& 1027 & mingf(1)*HaBohr_eVAng**2,& 1028 & mingf(2)*HaBohr_eVAng**2,& 1029 & mingf(3)*HaBohr_eVAng**2 1030 call wrtout(ab_out,message,'COLL') 1031 end if 1032 1033 ncycle_tot = ncycle_tot + 1 1034 1035 !Store GF Values of this iteration 1036 gf_values_iter(:,icycle_tmp+1) = mingf(:) 1037 1038 ! Check the stopping criterion 1039 converge = .false. 1040 if(tolGF > zero)then 1041 check_value = (sum(gf_values_iter(2:4,icycle_tmp+1),MASK=sel_on) - sum(gf_values_iter(2:4,icycle_tmp),MASK=sel_on)) & 1042 & /(sum(gf_values_iter(2:4,icycle_tmp+1),MASK=sel_on) - sum(gf_values_iter(2:4,1),MASK=sel_on)) 1043 if(check_value < tolGF)then 1044 write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",& 1045 & check_value ," < ",tolGF,& 1046 & ' Goal Function is converged' 1047 converge = .true. 1048 end if 1049 endif 1050 if(tolMSDE > zero)then 1051 if(abs(tolMSDE) > abs(mingf(4)* (Ha_eV *1000)**2 *factor))then 1052 write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",& 1053 & mingf(4)* (Ha_eV *1000)**2 * factor ," < ",tolMSDE,& 1054 & ' for MSDE' 1055 converge = .true. 1056 end if 1057 end if 1058 if(tolMSDF > zero) then 1059 if(abs(tolMSDF) > abs(mingf(2)*HaBohr_eVAng**2))then 1060 write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",& 1061 & mingf(2)*HaBohr_eVAng**2 ," < ",tolMSDF,& 1062 & ' for MSDF' 1063 converge = .true. 1064 end if 1065 end if 1066 if(tolMSDS > zero) then 1067 if(abs(tolMSDS) > abs(mingf(3)*HaBohr_eVAng**2))then 1068 write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",& 1069 & mingf(3)*HaBohr_eVAng**2 ," < ",tolMSDS,& 1070 & ' for MSDS' 1071 converge = .true. 1072 end if 1073 end if 1074 if(tolMSDFS > zero)then 1075 if(abs(tolMSDFS) > abs(mingf(1)*HaBohr_eVAng**2))then 1076 write(message,'(2a,ES18.10,a,ES18.10,a)') ch10," Fit process complete =>",& 1077 & mingf(1)*HaBohr_eVAng**2 ," < ",tolMSDFS,& 1078 & ' for MSDFS' 1079 converge = .true. 1080 end if 1081 end if 1082 if(converge)then 1083 call wrtout(ab_out,message,'COLL') 1084 call wrtout(std_out,message,'COLL') 1085 exit 1086 else 1087 if(any((/abs(tolMSDE),abs(tolMSDF),abs(tolMSDS),abs(tolMSDFS)/) > tol20) .and.& 1088 & icycle_tmp == ncycle)then 1089 write(message,'(2a,I0,a)') ch10," WARNING: ",ncycle,& 1090 & " cycles was not enougth to converge the fit process" 1091 call wrtout(ab_out,message,'COLL') 1092 call wrtout(std_out,message,'COLL') 1093 end if 1094 end if 1095 1096 end do !icycle_tmp=1,ncycle 1097 1098 case(2) 1099 1100 ! Monte Carlo selection 1101 nsweep = 10000 1102 ! If no coefficient imposed in the inputs we reset the goal function 1103 if (ncycle_tot == 0) then 1104 gf_values(:,:) = zero 1105 mingf(:) = 9D99 1106 else 1107 mingf = gf_values(:,1) 1108 end if 1109 call cpu_time(time) 1110 call ZBQLINI(int(time*1000000/(my_rank+1))) 1111 if(need_verbose)then 1112 write(message,'(a,I0,a)') " Start Monte Carlo simulations on ", nproc," CPU" 1113 if(nproc>1) write(message,'(2a)') trim(message)," (only print result of the master)" 1114 call wrtout(std_out,message,'COLL') 1115 call wrtout(ab_out,message,'COLL') 1116 write(message,'(a,2x,a,9x,a,14x,a,13x,a,14x,a)') ch10," Iteration ","MSDE","MSDFS","MSDF","MSdS" 1117 call wrtout(std_out,message,'COLL') 1118 write(message,'(a,5x,a,8x,a,8x,a,8x,a)') " ","(eV^2/A^2)","(eV^2/A^2)","(eV^2/A^2)",& 1119 & "(eV^2/A^2)" 1120 call wrtout(std_out,message,'COLL') 1121 1122 end if 1123 1124 do ii = 1,1!nyccle 1125 do isweep =1,nsweep 1126 write (j_char, '(i7)') isweep 1127 !TEST_AM 1128 icycle_tmp = int(ZBQLU01(zero)*(ncycle+1-1))+1 1129 icycle_tmp = ncycle 1130 do icycle=1,icycle_tmp 1131 icoeff = int(ZBQLU01(zero)*(my_ncoeff))+1 1132 ! icycle = int(ZBQLU01(zero)*(ncycle))+1 1133 list_coeffs_tmp2(icycle) = icoeff 1134 list_coeffs_tmp(icycle)= icycle 1135 ! Fill the temporary arrays 1136 energy_coeffs_tmp(icycle,:) = energy_coeffs(icoeff,:) 1137 fcart_coeffs_tmp(:,:,icycle,:) = fcart_coeffs(:,:,icoeff,:) 1138 strten_coeffs_tmp(:,:,icycle) = strten_coeffs(:,:,icoeff) 1139 end do 1140 !TEST_AM 1141 1142 ! call the fit process routine 1143 ! This routine solves the linear system proposed by 1144 ! C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]] 1145 call fit_polynomial_coeff_solve(coeff_values(1:icycle_tmp),fcart_coeffs_tmp,fit_data%fcart_diff,& 1146 & energy_coeffs_tmp,fit_data%energy_diff,info,& 1147 & list_coeffs_tmp(1:icycle_tmp),natom_sc,icycle_tmp,ncycle_max,& 1148 & ntime,strten_coeffs_tmp,fit_data%strten_diff,& 1149 & fit_data%training_set%sqomega,fit_on,int_fit_factors) 1150 if(info==0)then 1151 call fit_polynomial_coeff_computeGF(coeff_values(1:icycle_tmp),energy_coeffs_tmp,& 1152 & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,& 1153 & gf_values(:,1),list_coeffs_tmp(1:icycle_tmp),natom_sc,& 1154 & icycle_tmp,ncycle_max,ntime,strten_coeffs_tmp,& 1155 & fit_data%strten_diff,fit_data%training_set%sqomega) 1156 1157 else!In this case the matrix is singular 1158 gf_values(:,icoeff) = zero 1159 singular_coeffs(icoeff) = 1 1160 end if 1161 1162 if(gf_values(1,1) > zero.and.abs(gf_values(1,1))>tol16.and.& 1163 & gf_values(1,1) < mingf(1) ) then 1164 mingf = gf_values(:,1) 1165 list_coeffs(1:icycle_tmp) = list_coeffs_tmp2(1:icycle_tmp) 1166 ncycle_tot = icycle_tmp 1167 1168 write(message, '(4x,a,3x,4ES18.10)') adjustl(j_char),& 1169 & gf_values(4,1)* (1000*Ha_ev)**2 *factor,& 1170 & gf_values(1,1)*HaBohr_eVAng**2,& 1171 & gf_values(2,1)*HaBohr_eVAng**2,& 1172 & gf_values(3,1)*HaBohr_eVAng**2 1173 if(need_verbose) call wrtout(std_out,message,'COLL') 1174 else 1175 list_coeffs_tmp2(1:icycle_tmp) = list_coeffs(1:icycle_tmp) 1176 end if 1177 end do 1178 1179 if(nproc > 1) then 1180 !TEST_AM 1181 do iproc=1,ncycle_tot 1182 stat_coeff(list_coeffs(iproc)) = stat_coeff(list_coeffs(iproc)) + 1 1183 end do 1184 !TEST_AM 1185 1186 ! Find the best model on all the CPUs 1187 buffGF(1,1) = zero 1188 buffGF(2:5,1) = mingf(:) 1189 call xmpi_allgatherv(buffGF,5,gf_mpi,buffsize,buffdisp, comm, ierr) 1190 ! find the best coeff 1191 mingf(:) = 9D99 1192 index_min= 0 1193 do iproc=1,nproc 1194 if(gf_mpi(2,iproc) < zero) cycle 1195 if(abs(gf_mpi(2,iproc)) <tol16) cycle 1196 if(gf_mpi(2,iproc) < mingf(1) ) then 1197 mingf(:) = gf_mpi(2:5,iproc) 1198 index_min = int(gf_mpi(1,iproc)) 1199 rank_to_send = iproc-1 1200 end if 1201 end do 1202 write(message, '(2a,I0)') ch10,' Best model found on the CPU: ', rank_to_send 1203 call wrtout(std_out,message,'COLL') 1204 end if 1205 end do 1206 1207 !TEST_AM 1208 ! call xmpi_sum(stat_coeff, comm, ierr) 1209 ! do ii=1,ncoeff_tot 1210 ! write(100,*) ii,stat_coeff(ii) 1211 ! end do 1212 ! close(100) 1213 !TEST_AM 1214 1215 ! Transfer final model 1216 if(nproc>1)then 1217 call xmpi_bcast(ncycle_tot,rank_to_send,comm,ierr) 1218 call xmpi_bcast(list_coeffs(1:ncycle_tot),rank_to_send,comm,ierr) 1219 end if 1220 do ii=1,ncycle_tot 1221 icoeff = list_coeffs(ii) 1222 list_coeffs_tmp(ii) = ii 1223 ! Fill the temporary arrays 1224 energy_coeffs_tmp(ii,:) = energy_coeffs(icoeff,:) 1225 fcart_coeffs_tmp(:,:,ii,:) = fcart_coeffs(:,:,icoeff,:) 1226 strten_coeffs_tmp(:,:,ii) = strten_coeffs(:,:,icoeff) 1227 call polynomial_coeff_free(coeffs_tmp(ii)) 1228 call polynomial_coeff_init(one,my_coeffs(icoeff)%nterm,& 1229 & coeffs_tmp(ii),my_coeffs(icoeff)%terms,& 1230 & my_coeffs(icoeff)%name,& 1231 & check=.false.) 1232 end do 1233 end select 1234 1235 !This routine solves the linear system proposed by 1236 ! C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]] 1237 if(ncycle_tot > 0)then 1238 1239 call fit_polynomial_coeff_solve(coeff_values(1:ncycle_tot),fcart_coeffs_tmp,fit_data%fcart_diff,& 1240 & energy_coeffs_tmp,fit_data%energy_diff,info,& 1241 & list_coeffs_tmp(1:ncycle_tot),natom_sc,& 1242 & ncycle_tot,ncycle_max,ntime,strten_coeffs_tmp,& 1243 & fit_data%strten_diff,fit_data%training_set%sqomega,fit_on,int_fit_factors) 1244 1245 if(need_verbose) then 1246 write(message, '(3a)') ch10,' Fitted coefficients at the end of the fit process: ' 1247 call wrtout(ab_out,message,'COLL') 1248 call wrtout(std_out,message,'COLL') 1249 end if 1250 do ii = 1,ncycle_tot 1251 if(list_coeffs(ii) ==0) cycle 1252 ! Set the value of the coefficient 1253 coeffs_tmp(ii)%coefficient = coeff_values(ii) 1254 if(need_verbose) then 1255 write(message, '(a,I0,a,ES19.10,2a)') " ",list_coeffs(ii)," =>",coeff_values(ii),& 1256 & " ",trim(coeffs_tmp(ii)%name) 1257 call wrtout(ab_out,message,'COLL') 1258 call wrtout(std_out,message,'COLL') 1259 end if 1260 end do 1261 1262 call fit_polynomial_coeff_computeGF(coeff_values(1:ncycle_tot),energy_coeffs_tmp,& 1263 & fit_data%energy_diff,fcart_coeffs_tmp,fit_data%fcart_diff,& 1264 & gf_values(:,1),list_coeffs_tmp(1:ncycle_tot),natom_sc,& 1265 & ncycle_tot,ncycle_max,ntime,strten_coeffs_tmp,& 1266 & fit_data%strten_diff,fit_data%training_set%sqomega) 1267 1268 if(need_verbose) then 1269 ! Print the standard deviation after the fit 1270 write(message,'(4a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,& 1271 !& ' Mean Standard Deviation values at the end of the fit process (meV^2/atm): ',ch10,& 1272 !& ' Energy : ',& 1273 !& gf_values(4,1)*(Ha_EV*1000)**2 *factor ,ch10,& 1274 & ' Goal function values at the end of the fit process (eV^2/A^2):',ch10,& 1275 & ' Energy : ',& 1276 & gf_values(4,1)*(HaBohr_eVAng)**2,ch10,& 1277 & ' Forces+Stresses : ',& 1278 & gf_values(1,1)*(HaBohr_eVAng)**2,ch10,& 1279 & ' Forces : ',& 1280 & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,& 1281 & ' Stresses : ',& 1282 & gf_values(3,1)*(HaBohr_eVAng)**2,ch10 1283 call wrtout(ab_out,message,'COLL') 1284 call wrtout(std_out,message,'COLL') 1285 end if 1286 1287 !Allocate output coeffs -> selected plus fixed ones 1288 ncoeff_out = ncycle_tot+eff_pot_fixed%anharmonics_terms%ncoeff 1289 ! DEBUG MS 1290 ! write(*,*) ncoeff_out,ncycle_tot 1291 ABI_MALLOC(coeffs_out,(ncoeff_out)) 1292 do ii = 1,ncoeff_out 1293 if(ii <= eff_pot_fixed%anharmonics_terms%ncoeff)then 1294 call polynomial_coeff_init(eff_pot_fixed%anharmonics_terms%coefficients(ii)%coefficient,& 1295 & eff_pot_fixed%anharmonics_terms%coefficients(ii)%nterm,coeffs_out(ii),& 1296 & eff_pot_fixed%anharmonics_terms%coefficients(ii)%terms,& 1297 & eff_pot_fixed%anharmonics_terms%coefficients(ii)%name,& 1298 & check = .TRUE.) 1299 else 1300 ia = ii - eff_pot_fixed%anharmonics_terms%ncoeff 1301 call polynomial_coeff_init(coeffs_tmp(ia)%coefficient,coeffs_tmp(ia)%nterm,& 1302 & coeffs_out(ii),coeffs_tmp(ia)%terms,& 1303 & coeffs_tmp(ia)%name,& 1304 & check=.true.) 1305 1306 endif 1307 ! DEBUG MS 1308 ! write(*,*) "coeffs_out(", ii,")%name:",coeffs_out(ii)%name 1309 enddo 1310 1311 1312 ! Set the final set of coefficients into the eff_pot type 1313 call effective_potential_setCoeffs(coeffs_out(:),eff_pot,ncoeff_out) 1314 1315 ! If Wanted open the anharmonic_terms_file and write header 1316 filename = "TRS_fit_diff" 1317 ncoeff_model = eff_pot%anharmonics_terms%ncoeff 1318 if(need_prt_anh .and. ncoeff_model > 0 )then 1319 call effective_potential_writeAnhHead(ncoeff_model,filename,& 1320 & eff_pot%anharmonics_terms) 1321 else if (need_prt_anh)then 1322 write(message, '(6a,I3,3a)' )ch10,& 1323 & ' --- !WARNING',ch10,& 1324 & ' Printing of anharmonic terms has been asked,but',ch10,& 1325 & ' there are',ncoeff_model,'anharmonic terms in the potential',ch10,& 1326 & ' ---',ch10 1327 call wrtout(ab_out,message,'COLL') 1328 call wrtout(std_out,message,'COLL') 1329 end if 1330 1331 ! Calculate MSD values for final model 1332 if(need_prt_files)call fit_polynomial_coeff_computeMSD(eff_pot,hist,gf_values(4,1),gf_values(2,1),gf_values(1,1),& 1333 & natom_sc,ntime,fit_data%training_set%sqomega,comm,& 1334 & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=filename) 1335 1336 1337 INQUIRE(FILE='TRS_fit_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh) 1338 if(file_opened) close(unit_anh) 1339 1340 ! if(need_verbose) then 1341 ! ! Print the standard deviation after the fit 1342 ! write(message,'(4a,ES24.16,4a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,& 1343 ! & ' Mean Standard Deviation values at the end of the fit process (meV/f.u.):',& 1344 ! & ch10,' Energy : ',& 1345 ! & gf_values(4,1)*Ha_EV*1000/ ncell ,ch10,& 1346 ! & ' Goal function values at the end of the fit process (eV^2/A^2):',ch10,& 1347 ! & ' Forces+Stresses : ',& 1348 ! & (gf_values(1,1)+gf_values(2,1))*(HaBohr_eVAng)**2,ch10,& 1349 ! & ' Forces : ',& 1350 ! & gf_values(2,1)*(HaBohr_eVAng)**2,ch10,& 1351 ! & ' Stresses : ',& 1352 ! & gf_values(3,1)*(HaBohr_eVAng)**2,ch10 1353 ! call wrtout(ab_out,message,'COLL') 1354 ! call wrtout(std_out,message,'COLL') 1355 ! end if 1356 1357 else 1358 ncoeff_out = 0 1359 if(need_verbose) then 1360 write(message, '(9a)' )ch10,& 1361 & ' --- !WARNING',ch10,& 1362 & ' The fit process does not provide possible terms.',ch10,& 1363 & ' Please make sure that the terms set is correct',ch10,& 1364 & ' ---',ch10 1365 call wrtout(ab_out,message,'COLL') 1366 call wrtout(std_out,message,'COLL') 1367 end if 1368 end if 1369 !Deallocation of arrays 1370 call fit_data_free(fit_data) 1371 1372 !Deallocate the temporary coefficient 1373 do ii=1,ncycle_max 1374 call polynomial_coeff_free(coeffs_tmp(ii)) 1375 end do 1376 ABI_FREE(coeffs_tmp) 1377 do ii=1,my_ncoeff 1378 call polynomial_coeff_free(my_coeffs(ii)) 1379 end do 1380 ABI_FREE(my_coeffs) 1381 do ii=1,ncoeff_out 1382 call polynomial_coeff_free(coeffs_out(ii)) 1383 end do 1384 ABI_FREE(coeffs_out) 1385 1386 1387 !Deallocate fixed eff_pot 1388 call effective_potential_free(eff_pot_fixed) 1389 1390 1391 !Other deallocations 1392 ABI_FREE(gf_values_iter) 1393 ABI_FREE(buffsize) 1394 ABI_FREE(buffdisp) 1395 ABI_FREE(buffGF) 1396 ABI_FREE(coeff_values) 1397 ABI_FREE(energy_coeffs) 1398 ABI_FREE(energy_coeffs_tmp) 1399 ABI_FREE(fcart_coeffs) 1400 ABI_FREE(fcart_coeffs_tmp) 1401 ABI_FREE(gf_mpi) 1402 ABI_FREE(gf_values) 1403 ABI_FREE(list_coeffs) 1404 ABI_FREE(list_coeffs_tmp) 1405 ABI_FREE(list_coeffs_tmp2) 1406 ABI_FREE(my_coeffindexes) 1407 ABI_FREE(my_coefflist) 1408 ABI_FREE(singular_coeffs) 1409 ABI_FREE(strten_coeffs) 1410 ABI_FREE(strten_coeffs_tmp) 1411 ABI_FREE(stat_coeff) 1412 1413 ABI_SFREE(fixcoeff_corr) 1414 ABI_SFREE(fix_and_impose) 1415 ABI_SFREE(list_coeffs_copy) 1416 1417 end subroutine fit_polynomial_coeff_fit
m_fit_polynomial_coeff/fit_polynomial_coeff_getCoeffBound [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_getCoeffBound
FUNCTION
This routine fit a list of possible model. Return in the isPositive array:
INPUTS
NEED TO UPDATE eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = The history of the MD (or snapshot of DFT comm = MPI communicator verbose = optional, flag for the verbose mode
OUTPUT
SOURCE
1650 subroutine fit_polynomial_coeff_getCoeffBound(eff_pot,coeffs_out,hist,ncoeff_bound,comm,verbose) 1651 1652 implicit none 1653 1654 !Arguments ------------------------------------ 1655 !scalars 1656 integer,intent(in) :: comm 1657 integer,intent(out) :: ncoeff_bound 1658 logical,optional,intent(in) :: verbose 1659 !arrays 1660 type(abihist),intent(inout) :: hist 1661 type(effective_potential_type),target,intent(inout) :: eff_pot 1662 type(polynomial_coeff_type),allocatable,intent(out) :: coeffs_out(:) 1663 !Local variables------------------------------- 1664 !scalar 1665 integer :: counter,icoeff,icoeff_bound,idisp,istrain,ii 1666 integer :: istart,iterm,ndisp,nstrain,nterm,ncoeff_model,ncoeff_in,ncoeff_max 1667 real(dp):: weight 1668 logical :: need_verbose 1669 !arrays 1670 integer,allocatable :: atindx(:,:),cells(:,:,:),direction(:) 1671 integer,allocatable :: power_disps(:),power_strain(:),strain(:) 1672 type(polynomial_term_type),dimension(:),allocatable :: terms 1673 integer,allocatable :: odd_coeff(:),need_bound(:) 1674 type(polynomial_coeff_type),pointer :: coeffs_in(:) 1675 type(polynomial_coeff_type),allocatable :: coeffs_test(:) 1676 character(len=5),allocatable :: symbols(:) 1677 character(len=200):: name 1678 character(len=500) :: msg 1679 ! ************************************************************************* 1680 1681 1682 !set the inputs varaibles 1683 ncoeff_model = eff_pot%anharmonics_terms%ncoeff 1684 coeffs_in => eff_pot%anharmonics_terms%coefficients 1685 1686 !Do check 1687 if(ncoeff_model == 0)then 1688 write(msg,'(a)')'ncoeff_model must be different to 0' 1689 ABI_BUG(msg) 1690 end if 1691 1692 !Map the hist in order to be consistent with the supercell into reference_effective_potential 1693 call effective_potential_file_mapHistToRef(eff_pot,hist,comm) 1694 1695 !Initialisation of optional arguments 1696 need_verbose = .TRUE. 1697 if(present(verbose)) need_verbose = verbose 1698 1699 write(msg, '(a)' ) ' Detection of the unbound coefficients' 1700 if(need_verbose)call wrtout(std_out,msg,'COLL') 1701 1702 !Allocation 1703 ncoeff_max = 2 * ncoeff_model 1704 ABI_MALLOC(odd_coeff,(ncoeff_max)) 1705 ABI_MALLOC(need_bound,(ncoeff_max)) 1706 1707 ABI_MALLOC(symbols,(eff_pot%crystal%natom)) 1708 call symbols_crystal(eff_pot%crystal%natom,eff_pot%crystal%ntypat,eff_pot%crystal%npsp,& 1709 & symbols,eff_pot%crystal%typat,eff_pot%crystal%znucl) 1710 1711 1712 ABI_MALLOC(coeffs_test,(ncoeff_max)) 1713 1714 do icoeff=1,ncoeff_model 1715 call polynomial_coeff_init(coeffs_in(icoeff)%coefficient,coeffs_in(icoeff)%nterm,& 1716 & coeffs_test(icoeff),coeffs_in(icoeff)%terms,& 1717 & coeffs_in(icoeff)%name,check=.false.) 1718 end do 1719 1720 !array to know which coeff has to be bound 1721 need_bound(:) = 1 1722 counter = 0 1723 ncoeff_in = ncoeff_model 1724 1725 do while(.not.all(need_bound == 0).and.counter<1) 1726 ! Get the coefficients with odd coefficient 1727 odd_coeff = 0 1728 if(counter>0) then 1729 need_bound(1:ncoeff_in) = 0 1730 icoeff_bound = ncoeff_in 1731 else 1732 icoeff_bound = 1 1733 end if 1734 1735 do icoeff=icoeff_bound,ncoeff_model 1736 if(any(mod(coeffs_in(icoeff)%terms(1)%power_disp(:),2)/=0))then 1737 odd_coeff(icoeff) = 1 1738 end if 1739 if(any(mod(coeffs_in(icoeff)%terms(1)%power_strain(:),2)/=0))then 1740 odd_coeff(icoeff) = 1 1741 end if 1742 if(odd_coeff(icoeff) == 0 .and. coeffs_in(icoeff)%coefficient > zero) then 1743 need_bound(icoeff) = 0 1744 else 1745 need_bound(icoeff) = 1 1746 end if 1747 end do 1748 if(need_verbose)then 1749 write(msg, '(a)' ) ' The following coefficients need to be bound:' 1750 call wrtout(std_out,msg,'COLL') 1751 do icoeff=1,ncoeff_model 1752 if(need_bound(icoeff) == 1)then 1753 write(msg, '(2a)' ) ' =>',trim(coeffs_in(icoeff)%name) 1754 call wrtout(std_out,msg,'COLL') 1755 end if 1756 end do 1757 end if 1758 1759 1760 icoeff_bound = ncoeff_in + 1 1761 if(counter==0)then 1762 istart = 1 1763 else 1764 istart = ncoeff_in 1765 end if 1766 1767 ncoeff_bound = count(need_bound(istart:ncoeff_model)==1) 1768 1769 do icoeff=istart,ncoeff_model 1770 if(need_bound(icoeff)==1)then 1771 1772 nterm = coeffs_in(icoeff)%nterm 1773 ndisp = coeffs_in(icoeff)%terms(1)%ndisp 1774 nstrain = coeffs_in(icoeff)%terms(1)%nstrain 1775 1776 ABI_MALLOC(terms,(nterm)) 1777 ABI_MALLOC(atindx,(2,ndisp)) 1778 ABI_MALLOC(cells,(3,2,ndisp)) 1779 ABI_MALLOC(direction,(ndisp)) 1780 ABI_MALLOC(power_disps,(ndisp)) 1781 ABI_MALLOC(power_strain,(nstrain)) 1782 ABI_MALLOC(strain,(nstrain)) 1783 1784 do iterm=1,coeffs_in(icoeff)%nterm 1785 atindx(:,:) = coeffs_in(icoeff)%terms(iterm)%atindx(:,:) 1786 cells(:,:,:) = coeffs_in(icoeff)%terms(iterm)%cell(:,:,:) 1787 direction(:) = coeffs_in(icoeff)%terms(iterm)%direction(:) 1788 power_strain(:) = coeffs_in(icoeff)%terms(iterm)%power_strain(:) 1789 power_disps(:) = coeffs_in(icoeff)%terms(iterm)%power_disp(:) 1790 strain(:) = coeffs_in(icoeff)%terms(iterm)%strain(:) 1791 weight = 1 1792 do idisp=1,ndisp 1793 if(mod(power_disps(idisp),2) /= 0) then 1794 power_disps(idisp) = power_disps(idisp) + 1 1795 else 1796 power_disps(idisp) = power_disps(idisp) + 2 1797 end if 1798 end do 1799 do istrain=1,nstrain 1800 if(mod(power_strain(istrain),2) /= 0)then 1801 power_strain(istrain) = power_strain(istrain) + 1 1802 else 1803 if(power_strain(istrain) < 4 ) power_strain(istrain) = power_strain(istrain) + 2 1804 end if 1805 end do 1806 1807 call polynomial_term_init(atindx,cells,direction,ndisp,nstrain,terms(iterm),& 1808 & power_disps,power_strain,strain,weight,check=.true.) 1809 end do 1810 1811 name = "" 1812 call polynomial_coeff_init(one,nterm,coeffs_test(icoeff_bound),terms,name,check=.true.) 1813 call polynomial_coeff_getName(name,coeffs_test(icoeff_bound),symbols,recompute=.TRUE.) 1814 call polynomial_coeff_SetName(name,coeffs_test(icoeff_bound)) 1815 1816 ! Deallocate the terms 1817 do iterm=1,nterm 1818 call polynomial_term_free(terms(iterm)) 1819 end do 1820 ABI_FREE(terms) 1821 ABI_FREE(atindx) 1822 ABI_FREE(cells) 1823 ABI_FREE(direction) 1824 ABI_FREE(power_disps) 1825 ABI_FREE(power_strain) 1826 ABI_FREE(strain) 1827 1828 icoeff_bound = icoeff_bound + 1 1829 1830 end if 1831 end do 1832 1833 1834 if(counter==0)ncoeff_model = ncoeff_model + ncoeff_bound 1835 ! call effective_potential_setCoeffs(coeffs_test,eff_pot,ncoeff_model) 1836 call fit_polynomial_coeff_fit(eff_pot,(/0/),(/0/),hist,0,(/0,0/),1,0,& 1837 & -1,0,(/0/),1,comm,verbose=.true.,positive=.false.) 1838 1839 coeffs_in => eff_pot%anharmonics_terms%coefficients 1840 1841 counter = counter + 1 1842 end do 1843 1844 ABI_MALLOC(coeffs_out,(ncoeff_bound)) 1845 do ii=1,ncoeff_bound 1846 icoeff_bound = ncoeff_in + ii 1847 call polynomial_coeff_init(one,coeffs_test(icoeff_bound)%nterm,coeffs_out(ii),& 1848 & coeffs_test(icoeff_bound)%terms,coeffs_test(icoeff_bound)%name,check=.true.) 1849 end do 1850 !Deallocation 1851 do ii=ncoeff_model,ncoeff_max 1852 call polynomial_coeff_free(coeffs_test(ii)) 1853 end do 1854 1855 !Deallocation 1856 do icoeff=1,ncoeff_max 1857 call polynomial_coeff_free(coeffs_test(icoeff)) 1858 end do 1859 ABI_FREE(coeffs_test) 1860 ABI_FREE(odd_coeff) 1861 ABI_FREE(need_bound) 1862 ABI_FREE(symbols) 1863 1864 1865 end subroutine fit_polynomial_coeff_getCoeffBound
m_fit_polynomial_coeff/fit_polynomial_coeff_getFS [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_getFS
FUNCTION
Compute all the matrix elements of eq.11 and 12 in PRB95,094115 (2017) [[cite:Escorihuela-Sayalero2017]]
INPUTS
coefficients(ncoeff) = type(polynomial_coeff_type) du_delta(6,3,natom_sc,ntime) = Variation to displacements wrt to the strain (Bohr) displacement(3,natom_sc,ntime)= Atomic displacement wrt to the reference (Bohr) natom_sc = Number of atoms in the supercell natom_uc = Number of atoms in the unit cell ncoeff = Number of coefficients ntime = Number of time in the history sc_size(3) = Size of the supercell strain(6,ntime) = Strain ucvol(ntime) = Volume of the supercell for each time (Bohr^3) cells(ncell) = Indexes of the cell treat by this CPU ncell = Number of cell treat by this CPU index_cells(ncell,3) = Indexes of the cells (1 1 1, 0 0 0 for instance) treat by this CPU comm = MPI communicator
OUTPUT
fcart_out(ncoeff,3,natom,ntime) = value of the forces for each coefficient (-1 factor is taking into acount) (Ha/Bohr) strten_out(ncoeff,3,natom,ntime)= value of the stresses for each coefficient (-1/ucvol factor is taking into acount) (Ha/Bohr^3) energy_out(ncoeff,ntime) = value of the energy for each coefficient (Ha)
SOURCE
2230 subroutine fit_polynomial_coeff_getFS(coefficients,du_delta,displacement,energy_out,fcart_out,& 2231 & natom_sc,natom_uc,ncoeff_max,ntime,sc_size,strain,strten_out,& 2232 & ucvol,coeffs,ncoeff) 2233 2234 implicit none 2235 2236 !Arguments ------------------------------------ 2237 !scalars 2238 integer,intent(in) :: natom_sc,natom_uc,ncoeff_max,ntime 2239 integer,intent(in) :: ncoeff 2240 !arrays 2241 integer,intent(in) :: sc_size(3) 2242 integer,intent(in) :: coeffs(ncoeff_max) 2243 real(dp),intent(in) :: du_delta(6,3,natom_sc,ntime) 2244 real(dp),intent(in) :: displacement(3,natom_sc,ntime) 2245 real(dp),intent(in) :: strain(6,ntime),ucvol(ntime) 2246 real(dp),intent(out):: energy_out(ncoeff,ntime) 2247 real(dp),intent(out) :: fcart_out(3,natom_sc,ncoeff,ntime) 2248 real(dp),intent(out) :: strten_out(6,ntime,ncoeff) 2249 type(polynomial_coeff_type), intent(in) :: coefficients(ncoeff_max) 2250 !Local variables------------------------------- 2251 !scalar 2252 integer :: i1,i2,i3,ia1,ia2,ib1,ib2,ii,icell,icoeff,icoeff_tmp 2253 integer :: idir1,idir2,idisp1,idisp2,idisp1_strain,idisp2_strain 2254 integer :: iterm,itime,ndisp,ndisp_tot,nstrain,power_disp,power_strain 2255 real(dp):: disp1,disp2,tmp1,tmp2,tmp3,weight 2256 !arrays 2257 integer :: cell_atoma1(3),cell_atoma2(3) 2258 integer :: cell_atomb1(3),cell_atomb2(3) 2259 2260 ! ************************************************************************* 2261 2262 2263 !1-Get forces and stresses from the model 2264 ! Initialisation of variables 2265 fcart_out(:,:,:,:) = zero 2266 strten_out(:,:,:) = zero 2267 energy_out(:,:) = zero 2268 2269 icell = 0; ib1=0; ia1=0 2270 do i1=1,sc_size(1) 2271 do i2=1,sc_size(2) 2272 do i3=1,sc_size(3) 2273 ii = icell*natom_uc 2274 icell = icell + 1 2275 ! Loop over configurations 2276 do itime=1,ntime 2277 ! Loop over coefficients 2278 do icoeff_tmp=1,ncoeff 2279 icoeff = coeffs(icoeff_tmp) 2280 ! Loop over terms of this coefficient 2281 do iterm=1,coefficients(icoeff)%nterm 2282 ndisp = coefficients(icoeff)%terms(iterm)%ndisp 2283 nstrain = coefficients(icoeff)%terms(iterm)%nstrain 2284 ndisp_tot = ndisp + nstrain 2285 ! Set the weight of this term 2286 weight =coefficients(icoeff)%terms(iterm)%weight 2287 tmp1 = one 2288 ! Loop over displacement and strain 2289 do idisp1=1,ndisp_tot 2290 2291 ! Set to one the acculation of forces and strain 2292 tmp2 = one 2293 tmp3 = one 2294 ! Strain case idir => -6, -5, -4, -3, -2 or -1 2295 if (idisp1 > ndisp)then 2296 idisp1_strain = idisp1 - ndisp 2297 power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp1_strain) 2298 ! Get the direction of the displacement or strain 2299 idir1 = coefficients(icoeff)%terms(iterm)%strain(idisp1_strain) 2300 if(abs(strain(idir1,itime)) > tol10)then 2301 ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z)) 2302 tmp1 = tmp1 * (strain(idir1,itime))**power_strain 2303 if(power_strain > 1) then 2304 ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...)) 2305 tmp3 = tmp3 * power_strain*(strain(idir1,itime))**(power_strain-1) 2306 end if 2307 else 2308 tmp1 = zero 2309 if(power_strain > 1) then 2310 tmp3 = zero 2311 end if 2312 end if 2313 else 2314 ! Set the power_disp of the displacement: 2315 power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp1) 2316 ! Get the direction of the displacement or strain 2317 idir1 = coefficients(icoeff)%terms(iterm)%direction(idisp1) 2318 ! Displacement case idir = 1, 2 or 3 2319 ! indexes of the cell of the atom a 2320 cell_atoma1 = coefficients(icoeff)%terms(iterm)%cell(:,1,idisp1) 2321 if(cell_atoma1(1)/=0.or.cell_atoma1(2)/=0.or.cell_atoma1(3)/=0) then 2322 ! if the cell is not 0 0 0 we apply PBC: 2323 cell_atoma1(1) = i1 + cell_atoma1(1) 2324 cell_atoma1(2) = i2 + cell_atoma1(2) 2325 cell_atoma1(3) = i3 + cell_atoma1(3) 2326 call getPBCIndexes_supercell(cell_atoma1(1:3),sc_size(1:3)) 2327 ! index of the first atom (position in the supercell if the cell is not 0 0 0) 2328 ia1 = (cell_atoma1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 2329 & (cell_atoma1(2)-1)*sc_size(3)*natom_uc+& 2330 & (cell_atoma1(3)-1)*natom_uc+& 2331 & coefficients(icoeff)%terms(iterm)%atindx(1,idisp1) 2332 else 2333 ! index of the first atom (position in the supercell if the cell is 0 0 0) 2334 ia1 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp1) 2335 end if 2336 2337 ! indexes of the cell of the atom b (with PBC) same as ia1 2338 cell_atomb1 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp1) 2339 if(cell_atomb1(1)/=0.or.cell_atomb1(2)/=0.or.cell_atomb1(3)/=0) then 2340 cell_atomb1(1) = i1 + cell_atomb1(1) 2341 cell_atomb1(2) = i2 + cell_atomb1(2) 2342 cell_atomb1(3) = i3 + cell_atomb1(3) 2343 call getPBCIndexes_supercell(cell_atomb1(1:3),sc_size(1:3)) 2344 2345 ! index of the second atom in the (position in the supercell if the cell is not 0 0 0) 2346 ib1 = (cell_atomb1(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 2347 & (cell_atomb1(2)-1)*sc_size(3)*natom_uc+& 2348 & (cell_atomb1(3)-1)*natom_uc+& 2349 & coefficients(icoeff)%terms(iterm)%atindx(2,idisp1) 2350 else 2351 ! index of the first atom (position in the supercell if the cell is 0 0 0) 2352 ib1 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp1) 2353 end if 2354 2355 ! Get the displacement for the both atoms 2356 disp1 = displacement(idir1,ia1,itime) 2357 disp2 = displacement(idir1,ib1,itime) 2358 2359 if(abs(disp1) > tol10 .or. abs(disp2)> tol10)then 2360 ! Accumulate energy fo each displacement (\sum ((A_x-O_x)^Y(A_y-O_c)^Z)) 2361 tmp1 = tmp1 * (disp1-disp2)**power_disp 2362 if(power_disp > 1) then 2363 ! Accumulate forces for each displacement (\sum (Y(A_x-O_x)^Y-1(A_y-O_c)^Z+...)) 2364 tmp2 = tmp2 * power_disp*(disp1-disp2)**(power_disp-1) 2365 end if 2366 else 2367 tmp1 = zero 2368 if(power_disp > 1) then 2369 tmp2 = zero 2370 end if 2371 end if 2372 end if 2373 2374 do idisp2=1,ndisp_tot 2375 if(idisp2 /= idisp1) then 2376 2377 ! Strain case 2378 if (idisp2 > ndisp)then 2379 idisp2_strain = idisp2 - ndisp 2380 idir2 = coefficients(icoeff)%terms(iterm)%strain(idisp2_strain) 2381 ! Set the power_strain of the strain: 2382 power_strain = coefficients(icoeff)%terms(iterm)%power_strain(idisp2_strain) 2383 ! Accumulate energy forces 2384 tmp2 = tmp2 * (strain(idir2,itime))**power_strain 2385 ! Accumulate stress for each strain (\sum (Y(eta_2)^Y-1(eta_2)^Z+...)) 2386 tmp3 = tmp3 * (strain(idir2,itime))**power_strain 2387 ! Atomic displacement case 2388 else 2389 ! Set the power_disp of the displacement: 2390 power_disp = coefficients(icoeff)%terms(iterm)%power_disp(idisp2) 2391 ! Set the direction of the displacement: 2392 idir2 = coefficients(icoeff)%terms(iterm)%direction(idisp2) 2393 2394 cell_atoma2=coefficients(icoeff)%terms(iterm)%cell(:,1,idisp2) 2395 if(cell_atoma2(1)/=0.or.cell_atoma2(2)/=0.or.cell_atoma2(3)/=0) then 2396 cell_atoma2(1) = i1 + cell_atoma2(1) 2397 cell_atoma2(2) = i2 + cell_atoma2(2) 2398 cell_atoma2(3) = i3 + cell_atoma2(3) 2399 call getPBCIndexes_supercell(cell_atoma2(1:3),sc_size(1:3)) 2400 ! index of the first atom (position in the supercell and direction) 2401 ! if the cell of the atom a is not 0 0 0 (may happen) 2402 ia2 = (cell_atoma2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 2403 & (cell_atoma2(2)-1)*sc_size(3)*natom_uc+& 2404 & (cell_atoma2(3)-1)*natom_uc+& 2405 & coefficients(icoeff)%terms(iterm)%atindx(1,idisp2) 2406 else 2407 ! index of the first atom (position in the supercell and direction) 2408 ia2 = ii + coefficients(icoeff)%terms(iterm)%atindx(1,idisp2) 2409 end if 2410 2411 cell_atomb2 = coefficients(icoeff)%terms(iterm)%cell(:,2,idisp2) 2412 2413 if(cell_atomb2(1)/=0.or.cell_atomb2(2)/=0.or.cell_atomb2(3)/=0) then 2414 ! indexes of the cell2 (with PBC) 2415 cell_atomb2(1) = i1 + cell_atomb2(1) 2416 cell_atomb2(2) = i2 + cell_atomb2(2) 2417 cell_atomb2(3) = i3 + cell_atomb2(3) 2418 call getPBCIndexes_supercell(cell_atomb2(1:3),sc_size(1:3)) 2419 2420 ! index of the second atom in the (position in the supercell) 2421 ib2 = (cell_atomb2(1)-1)*sc_size(2)*sc_size(3)*natom_uc+& 2422 & (cell_atomb2(2)-1)*sc_size(3)*natom_uc+& 2423 & (cell_atomb2(3)-1)*natom_uc+& 2424 & coefficients(icoeff)%terms(iterm)%atindx(2,idisp2) 2425 else 2426 ib2 = ii + coefficients(icoeff)%terms(iterm)%atindx(2,idisp2) 2427 end if 2428 2429 disp1 = displacement(idir2,ia2,itime) 2430 disp2 = displacement(idir2,ib2,itime) 2431 2432 tmp2 = tmp2 * (disp1-disp2)**power_disp 2433 tmp3 = tmp3 * (disp1-disp2)**power_disp 2434 2435 end if 2436 end if 2437 end do 2438 2439 if(idisp1 > ndisp)then 2440 ! Accumule stress tensor 2441 strten_out(idir1,itime,icoeff_tmp) = strten_out(idir1,itime,icoeff_tmp) + & 2442 & weight * tmp3 / ucvol(itime) 2443 else 2444 ! Accumule forces 2445 fcart_out(idir1,ia1,icoeff_tmp,itime)=fcart_out(idir1,ia1,icoeff_tmp,itime)+weight*tmp2 2446 fcart_out(idir1,ib1,icoeff_tmp,itime)=fcart_out(idir1,ib1,icoeff_tmp,itime)-weight*tmp2 2447 end if 2448 end do 2449 2450 ! accumule energy 2451 energy_out(icoeff_tmp,itime) = energy_out(icoeff_tmp,itime) + weight * tmp1 2452 2453 end do!End do iterm 2454 end do!End do coeff 2455 end do!End time 2456 end do!End do i3 2457 end do!End do i2 2458 end do!End do i1 2459 2460 ! multiply by -1 2461 fcart_out(:,:,:,:) = -1 * fcart_out(:,:,:,:) 2462 2463 !ADD stress due to forces on atoms and variation of disp with strain 2464 do icoeff=1,ncoeff 2465 do itime=1,ntime 2466 do ia1=1,natom_sc 2467 do idir1=1,3 2468 do idir2=1,6 2469 strten_out(idir2,itime,icoeff) = strten_out(idir2,itime,icoeff) + & 2470 & du_delta(idir2,idir1,ia1,itime)*fcart_out(idir1,ia1,icoeff,itime)/ucvol(itime) 2471 end do 2472 end do 2473 end do 2474 end do 2475 end do 2476 2477 2478 end subroutine fit_polynomial_coeff_getFS
m_fit_polynomial_coeff/fit_polynomial_coeff_getPositive [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_getPositive
FUNCTION
This routine fit a list of possible model. Return in the isPositive array: 0 if the model ii does not contain possive coefficients 1 if the model ii contain possive coefficients
INPUTS
eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = The history of the MD (or snapshot of DFT coeff_values(nmodel,ncoeff) = values of the coefficients for each model isPositive(nmodel) = see description below list_coeff(nmodel,ncoeff) = list of the models ncoeff = number of coeff per model nfixcoeff = will not test the nfixcoeff first coeffcients nmodel = number of model comm = MPI communicator verbose = optional, flag for the verbose mode
OUTPUT
eff_pot = effective potential datatype with new fitted coefficients
SOURCE
1448 subroutine fit_polynomial_coeff_getPositive(eff_pot,hist,coeff_values,isPositive,list_coeff,ncoeff,& 1449 & nfixcoeff,nmodel,comm,verbose) 1450 1451 implicit none 1452 1453 !Arguments ------------------------------------ 1454 !scalars 1455 integer,intent(in) :: ncoeff,nfixcoeff,nmodel,comm 1456 !arrays 1457 integer,intent(in) :: list_coeff(nmodel,ncoeff) 1458 integer,intent(out) :: isPositive(nmodel) 1459 real(dp),intent(out) :: coeff_values(nmodel,ncoeff) 1460 type(effective_potential_type),intent(inout) :: eff_pot 1461 type(abihist),intent(inout) :: hist 1462 logical,optional,intent(in) :: verbose 1463 !Local variables------------------------------- 1464 !scalar 1465 integer :: ierr,ii,info,imodel,my_nmodel,nmodel_alone 1466 integer :: master,my_rank,ncoeff_tot,natom_sc,ncell 1467 integer :: nproc,ntime 1468 logical :: iam_master,need_verbose 1469 !arrays 1470 integer :: sc_size(3) 1471 integer,allocatable :: list_coeffs(:),my_modelindexes(:),my_modellist(:) 1472 real(dp),allocatable :: energy_coeffs(:,:),fcart_coeffs(:,:,:,:), strten_coeffs(:,:,:) 1473 type(polynomial_coeff_type),allocatable :: coeffs_in(:) 1474 type(fit_data_type) :: fit_data 1475 character(len=500) :: message 1476 logical :: fit_on(3) 1477 real(dp) :: fit_factors(3) 1478 ! ************************************************************************* 1479 1480 !MPI variables 1481 master = 0 1482 nproc = xmpi_comm_size(comm); my_rank = xmpi_comm_rank(comm) 1483 iam_master = (my_rank == master) 1484 1485 !Initialisation of optional arguments 1486 need_verbose = .TRUE. 1487 if(present(verbose)) need_verbose = verbose 1488 1489 fit_on(1) = .FALSE. 1490 fit_on(2) = .TRUE. 1491 fit_on(3) = .TRUE. 1492 1493 fit_factors = (/1,1,1/) 1494 1495 !Get the list of coefficients from the eff_pot 1496 if(eff_pot%anharmonics_terms%ncoeff > 0)then 1497 ! Copy the initial coefficients array 1498 ncoeff_tot = eff_pot%anharmonics_terms%ncoeff 1499 ABI_MALLOC(coeffs_in,(ncoeff_tot)) 1500 do ii=1,ncoeff_tot 1501 call polynomial_coeff_init(eff_pot%anharmonics_terms%coefficients(ii)%coefficient,& 1502 & eff_pot%anharmonics_terms%coefficients(ii)%nterm,& 1503 & coeffs_in(ii),& 1504 & eff_pot%anharmonics_terms%coefficients(ii)%terms,& 1505 & eff_pot%anharmonics_terms%coefficients(ii)%name,& 1506 & check=.false.) 1507 end do 1508 end if 1509 1510 !Reset the output (we free the memory) 1511 call effective_potential_freeCoeffs(eff_pot) 1512 1513 !if the number of atoms in reference supercell into effpot is not corret, 1514 !wrt to the number of atom in the hist, we set map the hist and set the good 1515 !supercell 1516 if (size(hist%xred,2) /= eff_pot%supercell%natom) then 1517 call effective_potential_file_mapHistToRef(eff_pot,hist,comm,verbose=need_verbose) 1518 end if 1519 1520 !Initialisation of constants 1521 natom_sc = eff_pot%supercell%natom 1522 ncell = eff_pot%supercell%ncells 1523 ntime = hist%mxhist 1524 do ii = 1, 3 1525 sc_size(ii) = eff_pot%supercell%rlatt(ii,ii) 1526 end do 1527 1528 !Initialisation of arrays: 1529 ABI_MALLOC(list_coeffs,(ncoeff_tot)) 1530 list_coeffs = 0 1531 do ii = 1,ncoeff_tot 1532 list_coeffs(ii) = ii 1533 end do 1534 1535 !Get the decomposition for each coefficients of the forces and stresses for 1536 !each atoms and each step equations 11 & 12 of PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]] 1537 if(need_verbose)then 1538 write(message, '(a)' ) ' Initialisation of the fit process...' 1539 call wrtout(std_out,message,'COLL') 1540 end if 1541 !Before the fit, compute constants with fit_data_compute. 1542 !Conpute the strain of each configuration. 1543 !Compute the displacmeent of each configuration. 1544 !Compute the variation of the displacement due to strain of each configuration. 1545 !Compute fixed forces and stresse and get the standard deviation. 1546 !Compute Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]. 1547 call fit_data_compute(fit_data,eff_pot,hist,comm,verbose=need_verbose) 1548 1549 !Get the decomposition for each coefficients of the forces,stresses and energy for 1550 !each atoms and each step (see equations 11 & 12 of 1551 ! PRB95,094115(2017)) [[cite:Escorihuela-Sayalero2017]] + allocation 1552 ABI_MALLOC(energy_coeffs,(ncoeff_tot,ntime)) 1553 ABI_MALLOC(fcart_coeffs,(3,natom_sc,ncoeff_tot,ntime)) 1554 ABI_MALLOC(strten_coeffs,(6,ntime,ncoeff_tot)) 1555 1556 call fit_polynomial_coeff_getFS(coeffs_in,fit_data%training_set%du_delta,& 1557 & fit_data%training_set%displacement,& 1558 & energy_coeffs,fcart_coeffs,natom_sc,eff_pot%crystal%natom,& 1559 & ncoeff_tot,ntime,sc_size,& 1560 & fit_data%training_set%strain,strten_coeffs,& 1561 & fit_data%training_set%ucvol,list_coeffs,ncoeff_tot) 1562 1563 1564 !set MPI, really basic stuff... 1565 nmodel_alone = mod(nmodel,nproc) 1566 my_nmodel = int(aint(real(nmodel,sp)/(nproc))) 1567 1568 if(my_rank >= (nproc-nmodel_alone)) then 1569 my_nmodel = my_nmodel + 1 1570 end if 1571 1572 ABI_MALLOC(my_modelindexes,(my_nmodel)) 1573 ABI_MALLOC(my_modellist,(my_nmodel)) 1574 1575 !2:compute the number of model and the list of the corresponding for each CPU. 1576 do imodel=1,my_nmodel 1577 if(my_rank >= (nproc-nmodel_alone))then 1578 my_modelindexes(imodel)=(int(aint(real(nmodel,sp)/nproc)))*(my_rank)+& 1579 & (my_rank - (nproc-nmodel_alone)) + imodel 1580 my_modellist(imodel) = imodel 1581 else 1582 my_modelindexes(imodel)=(my_nmodel)*(my_rank) + imodel 1583 my_modellist(imodel) = imodel 1584 end if 1585 end do 1586 1587 1588 !Start fit process 1589 isPositive = 0 1590 coeff_values = zero 1591 do ii=1,my_nmodel 1592 imodel = my_modelindexes(ii) 1593 call fit_polynomial_coeff_solve(coeff_values(imodel,1:ncoeff),fcart_coeffs,fit_data%fcart_diff,& 1594 & energy_coeffs,fit_data%energy_diff,info,& 1595 & list_coeff(imodel,1:ncoeff),natom_sc,ncoeff,& 1596 & ncoeff_tot,ntime,strten_coeffs,fit_data%strten_diff,& 1597 & fit_data%training_set%sqomega,fit_on,fit_factors) 1598 1599 if(info==0)then 1600 1601 if (any(coeff_values(imodel,nfixcoeff+1:ncoeff) < zero))then 1602 ! coeff_values(imodel,:) = zero 1603 isPositive(imodel) = 0 1604 else 1605 isPositive(imodel) = 1 1606 end if 1607 end if 1608 end do 1609 1610 call xmpi_sum(isPositive, comm, ierr) 1611 call xmpi_sum(coeff_values, comm, ierr) 1612 1613 !Deallocation of arrays 1614 do ii=1,ncoeff_tot 1615 call polynomial_coeff_free(coeffs_in(ii)) 1616 end do 1617 call fit_data_free(fit_data) 1618 ABI_FREE(coeffs_in) 1619 ABI_FREE(energy_coeffs) 1620 ABI_FREE(fcart_coeffs) 1621 ABI_FREE(list_coeffs) 1622 ABI_FREE(my_modelindexes) 1623 ABI_FREE(my_modellist) 1624 ABI_FREE(strten_coeffs) 1625 1626 end subroutine fit_polynomial_coeff_getPositive
m_fit_polynomial_coeff/fit_polynomial_coeff_solve [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_coeff_solve
FUNCTION
Build and the solve the system to get the values of the coefficients This routine solves the linear system proposed by C.Escorihuela-Sayalero see PRB95,094115(2017) [[cite:Escorihuela-Sayalero2017]]
INPUTS
fcart_coeffs(3,natom_sc,ncoeff_max,ntime) = List of the values of the contribution to the cartesian forces for all coefficients for each direction and each time fcart_diff(3,natom,ntime) = Difference of cartesian forces between DFT calculation and fixed part of the model (more often harmonic part) energy_coeffs(ncoeff,ntime) = value of the energy for each coefficient (Ha) energy_diff(ntime) = Difference of energ ybetween DFT calculation and fixed part of the model (more often harmonic part) list_coeffs(ncoeff_fit) = List with the index of the coefficients used for this model natom = Number of atoms ncoeff_fit = Number of coeff for the fit (dimension of the system) ncoeff_max = Maximum number of coeff in the list ntime = Number of time (number of snapshot, number of md step...) strten_coeffs(6,ntime,ncoeff_max) = List of the values of the contribution to the stress tensor of the coefficients for each direction,time strten_diff(6,natom) = Difference of stress tensor between DFT calculation and fixed part of the model (more often harmonic part) sqomega(ntime) = Sheppard and al Factors \Omega^{2} see J.Chem Phys 136, 074103 (2012) [[cite:Sheppard2012]]
OUTPUT
coefficients(ncoeff_fit) = Values of the coefficients
info_out = 0: successful exit
< 0: if INFO = -i, the i-th argument had an illegal value
> 0: if INFO = i, U(i,i) computed in DOUBLE PRECISION is
exactly zero. The factorization has been completed,
but the factor U is exactly singular, so the solution
could not be computed. = 0: successful exit
information from the subroutine dsgesv in LAPACK
SOURCE
1911 subroutine fit_polynomial_coeff_solve(coefficients,fcart_coeffs,fcart_diff,energy_coeffs,energy_diff,& 1912 & info_out,list_coeffs,natom,ncoeff_fit,ncoeff_max,ntime,& 1913 & strten_coeffs,strten_diff,sqomega,fit_on,fit_factors) 1914 1915 implicit none 1916 1917 !Arguments ------------------------------------ 1918 !scalars 1919 integer,intent(in) :: natom,ncoeff_fit,ncoeff_max,ntime 1920 integer,intent(out) :: info_out 1921 !arrays 1922 real(dp),intent(in) :: energy_coeffs(ncoeff_max,ntime) 1923 real(dp),intent(in) :: energy_diff(ntime) 1924 integer,intent(in) :: list_coeffs(ncoeff_fit) 1925 real(dp),intent(in) :: fcart_coeffs(3,natom,ncoeff_max,ntime) 1926 real(dp),intent(in) :: fcart_diff(3,natom,ntime) 1927 real(dp),intent(in) :: strten_coeffs(6,ntime,ncoeff_max) 1928 real(dp),intent(in) :: strten_diff(6,ntime),sqomega(ntime) 1929 real(dp),intent(out):: coefficients(ncoeff_fit) 1930 real(dp),intent(in) :: fit_factors(3) 1931 logical,intent(in) :: fit_on(3) 1932 !Local variables------------------------------- 1933 !scalar 1934 integer :: ia,itime,icoeff,jcoeff,icoeff_tmp,jcoeff_tmp,mu,LDA,LDB,LDX,LDAF,N,NRHS 1935 real(dp):: efact,ffact,sfact,ftmpA,stmpA,ftmpB,stmpB,etmpA,etmpB,fmu,fnu,smu,snu,emu,enu 1936 integer :: INFO,ITER 1937 real(dp):: RCOND 1938 real(dp):: fcart_coeffs_tmp(3,natom,ntime) 1939 real(dp),allocatable:: AF(:,:),BERR(:),FERR(:),WORK(:),C(:),R(:) 1940 integer,allocatable :: IPIV(:),IWORK(:),SWORK(:) 1941 !arrays 1942 real(dp),allocatable :: A(:,:),B(:,:) 1943 character(len=1) :: FACT,EQUED,TRANS 1944 ! character(len=500) :: message 1945 ! ************************************************************************* 1946 1947 !0-Set variables for the 1948 N = ncoeff_fit; NRHS = 1; LDA = ncoeff_fit; LDB = ncoeff_fit; LDX = ncoeff_fit 1949 LDAF = ncoeff_fit; RCOND = zero; INFO = 0; TRANS='N'; EQUED='N'; FACT='N' 1950 1951 !Set the factors 1952 efact = fit_factors(1)*one/(ntime) 1953 ffact = fit_factors(2)*one/(3*natom*ntime) 1954 sfact = fit_factors(3)*one/(6*ntime) 1955 1956 !0-Allocation 1957 ABI_MALLOC(A,(LDA,N)) 1958 ABI_MALLOC(B,(LDB,NRHS)) 1959 ABI_MALLOC(AF,(LDAF,N)) 1960 ABI_MALLOC(IPIV,(N)) 1961 ABI_MALLOC(R,(N)) 1962 ABI_MALLOC(C,(N)) 1963 ABI_MALLOC(FERR,(NRHS)) 1964 ABI_MALLOC(BERR,(NRHS)) 1965 ABI_MALLOC(WORK,(4*N)) 1966 ABI_MALLOC(IWORK,(N)) 1967 ABI_MALLOC(SWORK,(N*(N+NRHS))) 1968 A=zero; B=zero; 1969 AF = zero; IPIV = 1; 1970 R = one; C = one; 1971 FERR = zero; BERR = zero 1972 IWORK = 0; WORK = 0 1973 1974 !1-Get forces and stresses from the model and fill A 1975 ! Fill alsor B with the forces and stresses from 1976 ! the DFT snapshot and the model 1977 ! See equation 17 of PRB95 094115 (2017) [[cite:Escorihuela-Sayalero2017]] 1978 do icoeff=1,ncoeff_fit 1979 icoeff_tmp = list_coeffs(icoeff) 1980 fcart_coeffs_tmp(:,:,:) = fcart_coeffs(:,:,icoeff_tmp,:) 1981 ftmpA= zero; ftmpB = zero 1982 stmpA= zero; stmpB = zero 1983 etmpA= zero; etmpB = zero 1984 ! loop over the configuration 1985 do itime=1,ntime 1986 ! Fill energy 1987 emu = energy_coeffs(icoeff_tmp,itime) 1988 do jcoeff=1,ncoeff_fit 1989 jcoeff_tmp = list_coeffs(jcoeff) 1990 enu = energy_coeffs(jcoeff_tmp,itime) 1991 if(fit_on(3))then 1992 etmpA = emu*enu/(sqomega(itime)**(1.0/2.0)) 1993 A(icoeff,jcoeff) = A(icoeff,jcoeff) + efact*etmpA 1994 endif 1995 end do 1996 if(fit_on(3))then 1997 etmpB = etmpB + energy_diff(itime)*emu/(sqomega(itime)**(1.0/2.0)) !/ (sqomega(itime)**3) 1998 else 1999 etmpB = zero ! REMOVE THIS LINE TO TAKE INTO ACOUNT THE ENERGY 2000 endif 2001 ! Fill forces 2002 do ia=1,natom 2003 do mu=1,3 2004 fmu = fcart_coeffs_tmp(mu,ia,itime) 2005 do jcoeff=1,ncoeff_fit 2006 jcoeff_tmp = list_coeffs(jcoeff) 2007 fnu = fcart_coeffs(mu,ia,jcoeff_tmp,itime) 2008 ftmpA = fmu*fnu 2009 if(fit_on(1))A(icoeff,jcoeff) = A(icoeff,jcoeff) + ffact*ftmpA 2010 end do 2011 if(fit_on(1))then 2012 ftmpB = ftmpB + fcart_diff(mu,ia,itime)*fmu 2013 else 2014 ftmpB = zero 2015 endif 2016 end do !End loop dir 2017 end do !End loop natom 2018 ! Fill stresses 2019 do mu=1,6 2020 smu = strten_coeffs(mu,itime,icoeff_tmp) 2021 do jcoeff=1,ncoeff_fit 2022 jcoeff_tmp = list_coeffs(jcoeff) 2023 snu = strten_coeffs(mu,itime,jcoeff_tmp) 2024 stmpA = sqomega(itime)*smu*snu 2025 if(fit_on(2))A(icoeff,jcoeff) = A(icoeff,jcoeff) + sfact*stmpA 2026 end do 2027 if(fit_on(2))then 2028 stmpB = stmpB + sqomega(itime)*strten_diff(mu,itime)*smu 2029 else 2030 stmpB = zero 2031 endif 2032 end do !End loop stress dir 2033 end do ! End loop time 2034 B(icoeff,1) = B(icoeff,1) + ffact*ftmpB + sfact*stmpB + efact*etmpB 2035 end do ! End loop icoeff 2036 2037 !2-Solve Ax=B 2038 !OLD VERSION.. 2039 ! call dgesvx(FACT,TRANS,N,NRHS,A,LDA,AF,LDAF,IPIV,EQUED,R,C,B,LDB,coefficients,LDX,& 2040 ! RCOND,FERR,BERR,WORK,IWORK,INFO) 2041 !U is nonsingular 2042 ! if (INFO==N+1) then 2043 ! coefficients = zero 2044 ! end if 2045 call DSGESV(N,NRHS,A,LDA,IPIV,B,LDB,coefficients,LDX,WORK,SWORK,ITER,INFO) 2046 2047 !other routine 2048 ! call dgesv(N,NRHS,A,LDA,IPIV,B,LDB,INFO) 2049 ! coefficients = B(:,NRHS) 2050 !U is nonsingular 2051 if (INFO==N+2) then 2052 coefficients = zero 2053 end if 2054 2055 if(any(abs(coefficients)>1.0E10))then 2056 INFO = 1 2057 coefficients = zero 2058 end if 2059 2060 info_out = INFO 2061 2062 ABI_FREE(AF) 2063 ABI_FREE(IPIV) 2064 ABI_FREE(R) 2065 ABI_FREE(C) 2066 ABI_FREE(FERR) 2067 ABI_FREE(BERR) 2068 ABI_FREE(WORK) 2069 ABI_FREE(IWORK) 2070 ABI_FREE(SWORK) 2071 ABI_FREE(A) 2072 ABI_FREE(B) 2073 2074 end subroutine fit_polynomial_coeff_solve
m_fit_polynomial_coeff/fit_polynomial_printSystemFiles [ Functions ]
[ Top ] [ m_fit_polynomial_coeff ] [ Functions ]
NAME
fit_polynomial_printSystemFiles
FUNCTION
Print the files for the fitting script
INPUTS
eff_pot<type(effective_potential)> = effective potential hist<type(abihist)> = datatype with the history of the MD
OUTPUT
SOURCE
2838 subroutine fit_polynomial_printSystemFiles(eff_pot,hist) 2839 2840 implicit none 2841 2842 !Arguments ------------------------------------ 2843 !scalars 2844 !arrays 2845 type(effective_potential_type), intent(in) :: eff_pot 2846 type(abihist),intent(in) :: hist 2847 !Local variables------------------------------- 2848 !scalar 2849 integer :: ia,ib,ib1,ii,jj,irpt,kk,ll,mu,nu,nstep,nshift 2850 integer :: natom_uc 2851 integer :: unit_born=22,unit_epsiloninf=23,unit_md=24 2852 integer :: unit_harmonic=25,unit_ref=26,unit_strain=27,unit_sym=28 2853 !arrays 2854 integer,allocatable :: typat_order(:),typat_order_uc(:) 2855 integer, dimension(3) :: A,ncell 2856 real(dp), allocatable :: xcart(:,:),fcart(:,:) 2857 character(len=500) :: msg 2858 type(supercell_type) :: supercell 2859 ! ************************************************************************* 2860 2861 !Create new supercell corresponding to the MD 2862 ncell = (/2,2,2/) 2863 call init_supercell(eff_pot%crystal%natom, (/ncell(1),0,0, 0,ncell(2),0, 0,0,ncell(3)/),& 2864 & eff_pot%crystal%rprimd,eff_pot%crystal%typat,& 2865 & eff_pot%crystal%xcart,eff_pot%crystal%znucl, supercell) 2866 2867 !allocation of array 2868 ABI_MALLOC(xcart,(3,supercell%natom)) 2869 ABI_MALLOC(fcart,(3,supercell%natom)) 2870 ABI_MALLOC(typat_order,(supercell%natom)) 2871 ABI_MALLOC(typat_order_uc,(eff_pot%crystal%natom)) 2872 2873 A = (/ 2, 3, 1/) 2874 2875 nshift = product(ncell) 2876 natom_uc = eff_pot%crystal%natom 2877 !Fill the typat_order array: 2878 !In the fit script the atom must be in the order 11111 222222 33333 .. 2879 !and the order of the atom can not be change in the fit script, 2880 !we transform into the format of the script 2881 ib = 1 2882 ib1= 1 2883 do ii=1,eff_pot%crystal%ntypat 2884 jj = A(ii) 2885 do kk=1,natom_uc 2886 if(supercell%typat(kk)==jj)then 2887 typat_order_uc(ib1) = kk 2888 ib1 = ib1 + 1 2889 do ll=1,nshift 2890 ia = (ll-1)*natom_uc + kk 2891 typat_order(ib) = ia 2892 ib = ib + 1 2893 end do 2894 end if 2895 end do 2896 end do 2897 2898 ! BORN CHARGES FILE 2899 if (open_file('system/Born_Charges',msg,unit=unit_born,form="formatted",& 2900 & status="replace",action="write") /= 0) then 2901 ABI_ERROR(msg) 2902 end if 2903 do ii=1,eff_pot%crystal%ntypat 2904 jj = A(ii) 2905 do ia=1,eff_pot%crystal%natom 2906 if(eff_pot%crystal%typat(ia)==jj)then 2907 write(unit_born,'(i2,a,1F10.5)') ia," ",eff_pot%crystal%amu(eff_pot%crystal%typat(ia)) 2908 do mu=1,3 2909 WRITE(unit_born,'(a,3(F23.14))') " ",eff_pot%harmonics_terms%zeff(:,mu,ia) 2910 end do 2911 end if 2912 end do 2913 end do 2914 2915 !DIELECTRIC TENSOR FILE 2916 if (open_file('system/Dielectric_Tensor',msg,unit=unit_epsiloninf,form="formatted",& 2917 & status="replace",action="write") /= 0) then 2918 ABI_ERROR(msg) 2919 end if 2920 do mu=1,3 2921 WRITE(unit_epsiloninf,'(3(F23.14))') eff_pot%harmonics_terms%epsilon_inf(:,mu) 2922 end do 2923 2924 2925 !REFERENCE STRUCTURE FILE 2926 if (open_file('system/Reference_structure',msg,unit=unit_ref,form="formatted",& 2927 & status="replace",action="write") /= 0) then 2928 ABI_ERROR(msg) 2929 end if 2930 2931 write(unit_ref,'("Energy (Hartree)")') 2932 write(unit_ref,'("================")') 2933 write(unit_ref,'(F23.14)') (hist%etot(1)/nshift) 2934 write(unit_ref,'("")') 2935 write(unit_ref,'("Cell vectors")') 2936 write(unit_ref,'("============")') 2937 do jj=1,3 2938 write(unit_ref,'(3(F22.14))') (supercell%rprimd(:,jj)) 2939 end do 2940 2941 write(unit_ref,'("")') 2942 write(unit_ref,'("Atomic positions (Bohr radius)")') 2943 write(unit_ref,'("==============================")') 2944 2945 do ia=1,supercell%natom 2946 write(unit_ref,'(3(F23.14))') supercell%xcart(:,typat_order(ia)) 2947 end do 2948 2949 !Harmonic XML file 2950 if (open_file('system/harmonic.xml',msg,unit=unit_harmonic,form="formatted",& 2951 & status="replace",action="write") /= 0) then 2952 ABI_ERROR(msg) 2953 end if 2954 2955 !Write header 2956 write(unit_harmonic,'("<?xml version=""1.0"" ?>")') 2957 write(unit_harmonic,'("<name>")') 2958 2959 do irpt=1,eff_pot%harmonics_terms%ifcs%nrpt 2960 if(any(abs(eff_pot%harmonics_terms%ifcs%short_atmfrc(:,:,:,:,irpt))>tol9)) then 2961 write(unit_harmonic,'(" <local_force_constant units=""hartree/bohrradius**2"">")') 2962 write(unit_harmonic,'(" <data>")') 2963 do ia=1,eff_pot%crystal%natom 2964 do mu=1,3 2965 do ib=1,eff_pot%crystal%natom 2966 do nu=1,3 2967 write(unit_harmonic,'(F22.14)', advance="no")& 2968 & (eff_pot%harmonics_terms%ifcs%short_atmfrc(mu,typat_order_uc(ia),& 2969 & nu,typat_order_uc(ib),irpt)) 2970 end do 2971 end do 2972 write(unit_harmonic,'(a)')'' 2973 end do 2974 end do 2975 write(unit_harmonic,'(" </data>")') 2976 write(unit_harmonic,'(" <cell>")') 2977 write(unit_harmonic,'(3(I4))') (eff_pot%harmonics_terms%ifcs%cell(:,irpt)) 2978 write(unit_harmonic,'(" </cell>")') 2979 write(unit_harmonic,'(" </local_force_constant>")') 2980 end if 2981 end do 2982 write(unit_harmonic,'("</name>")') 2983 2984 !STRAIN FILE 2985 if (open_file('system/Strain_Tensor',msg,unit=unit_strain,form="formatted",& 2986 & status="replace",action="write") /= 0) then 2987 ABI_ERROR(msg) 2988 end if 2989 write(unit_strain,'(6(F23.14))') (eff_pot%harmonics_terms%elastic_constants) 2990 2991 ! SYM FILE 2992 if (open_file('system/symmetry_operations',msg,unit=unit_sym,form="formatted",& 2993 & status="replace",action="write") /= 0) then 2994 ABI_ERROR(msg) 2995 end if 2996 write(unit_sym,'("(x,y,z) (y,-x,z) (z,x,y) (y,z,x) (x,z,y) (y,x,z) (z,y,x) (x,-y,-z) (z,-x,-y)",& 2997 & " (y,-z,-x) (x,-z,-y) (y,-x,-z) (z,-y,-x) (-x,y,-z) (-z,x,-y) (-y,z,-x) (-x,z,-y)",& 2998 & " (-y,x,-z) (-z,y,-x) (-x,-y,z) (-z,-x,y) (-y,-z,x) (-x,-z,y) (-y,-x,z) (-z,-y,x)",& 2999 & " (-x,-y,-z) (-z,-x,-y) (-y,-z,-x) (-x,-z,-y) (-y,-x,-z) (-z,-y,-x) (-x,y,z)",& 3000 & " (-z,x,y) (-y,z,x) (-x,z,y) (-y,x,z) (-z,y,x) (x,-y,z) (z,-x,y) (y,-z,x) (x,-z,y)",& 3001 & " (z,-y,x) (x,y,-z) (z,x,-y) (y,z,-x) (x,z,-y) (y,x,-z) (z,y,-x)")') 3002 3003 3004 !MD file 3005 nstep = hist%mxhist 3006 if (open_file('system/Molecular_dynamic',msg,unit=unit_md,form="formatted",& 3007 & status="replace",action="write") /= 0) then 3008 ABI_ERROR(msg) 3009 end if 3010 do ii=1,nstep 3011 write(unit_md,'(I5)') ii-1 3012 write(unit_md,'(F22.14)') hist%etot(ii)/nshift 3013 do jj=1,3 3014 write(unit_md,'(3(F22.14))') (hist%rprimd(:,jj,ii)) 3015 end do 3016 ! Set xcart and fcart for this step 3017 call xred2xcart(supercell%natom,hist%rprimd(:,:,ii),& 3018 & xcart,hist%xred(:,:,ii)) 3019 3020 fcart(:,:) = hist%fcart(:,:,ii) 3021 3022 do ia=1,supercell%natom 3023 write(unit_md,'(3(E22.14),3(E22.14))') xcart(:,typat_order(ia)),fcart(:,typat_order(ia)) 3024 end do 3025 write(unit_md,'(6(E22.14))') hist%strten(:,ii) 3026 end do 3027 3028 !Close files 3029 close(unit_ref) 3030 close(unit_born) 3031 close(unit_harmonic) 3032 close(unit_epsiloninf) 3033 close(unit_md) 3034 close(unit_strain) 3035 close(unit_sym) 3036 3037 !Deallocation array 3038 ABI_FREE(typat_order) 3039 ABI_FREE(typat_order_uc) 3040 ABI_FREE(xcart) 3041 ABI_FREE(fcart) 3042 call destroy_supercell(supercell) 3043 3044 end subroutine fit_polynomial_printSystemFiles
m_fit_polynomiaL_coeff/testEffPot [ Functions ]
NAME
testEffPot
FUNCTION
Calculate the energy, forces for displacements provided in an test-set (input:hist) within a given effective potential (input: eff_pot) If the test set is from DFT and contains DFT energies and forces calculate the Goal Function values and the MSD of the Energy with respect to the DFT energies
INPUTS
eff_pot = effective_potential datatype hist = abihist datatype
OUTPUT
SOURCE
2715 subroutine fit_polynomial_coeff_testEffPot(eff_pot,hist,master,comm,print_anharmonic,scup_dtset,prt_ph) 2716 2717 2718 implicit none 2719 2720 !Arguments ------------------------------------ 2721 !scalars 2722 integer,intent(in) :: master,comm 2723 integer,optional,intent(in) :: prt_ph 2724 !logicals 2725 logical,optional,intent(in) :: print_anharmonic 2726 !array 2727 type(effective_potential_type),intent(inout) :: eff_pot 2728 type(abihist),intent(in) :: hist 2729 type(scup_dtset_type),optional,intent(inout) :: scup_dtset 2730 !Local variables------------------------------- 2731 !reals 2732 real(dp) :: factor,mse,msef,mses 2733 real(dp),allocatable :: sqomega(:),ucvol(:) 2734 real(dp),parameter :: HaBohr_eVAng = 27.21138386d0 / 0.529177249d0 2735 !scalar 2736 integer :: itime,unit_anh 2737 integer :: natom,ntime,ncoeff,my_rank 2738 !logicals 2739 logical :: iam_master, need_print_anharmonic,file_opened,need_prt_ph 2740 !strings/characters 2741 character(len=fnlen) :: filename 2742 character(len=1000) :: message 2743 !arrays 2744 real(dp) :: gmet(3,3),gprimd(3,3),rmet(3,3) 2745 ! ************************************************************************* 2746 2747 !MPI variables 2748 my_rank=xmpi_comm_rank(comm) 2749 iam_master = (my_rank == master) 2750 2751 !Initialisation of optional arguments 2752 need_print_anharmonic = .FALSE. 2753 if(present(print_anharmonic)) need_print_anharmonic = print_anharmonic 2754 2755 need_prt_ph = .FALSE. 2756 if(present(prt_ph))then 2757 if(prt_ph==1) need_prt_ph=.TRUE. 2758 endif 2759 2760 !Setting/Allocating other Variables 2761 natom = size(hist%xred,2) 2762 factor = 1._dp/natom 2763 ntime = hist%mxhist 2764 ABI_MALLOC(sqomega,(ntime)) 2765 ABI_MALLOC(ucvol,(ntime)) 2766 sqomega = zero 2767 filename = 'TES_fit_diff' 2768 ncoeff = eff_pot%anharmonics_terms%ncoeff 2769 2770 do itime=1,ntime 2771 ! Compute \Omega^{2} and ucvol for each time 2772 call metric(gmet,gprimd,-1,rmet,hist%rprimd(:,:,itime),ucvol(itime)) 2773 ! Formula: sqomega(itime) = (((ucvol(itime)**(-2.))* ((natom)**(0.5)))**(-1.0/3.0))**2 2774 ! Compact form: 2775 sqomega(itime) = ((ucvol(itime)**(4.0/3.0)) / ((natom)**(1/3.0))) 2776 end do 2777 2778 2779 if(need_print_anharmonic) call effective_potential_writeAnhHead(ncoeff,& 2780 & filename,eff_pot%anharmonics_terms) 2781 2782 call fit_polynomial_coeff_computeMSD(eff_pot,hist,mse,msef,mses,natom,ntime,& 2783 & sqomega,comm,& 2784 & compute_anharmonic=.TRUE.,print_file=.TRUE.,filename=filename,scup_dtset=scup_dtset,prt_ph=need_prt_ph) 2785 2786 2787 ! Print the standard deviation after the fit 2788 write(message,'(6a,ES24.16,2a,ES24.16,2a,ES24.16,2a,ES24.16,a)' )ch10,& 2789 !& ' Mean Standard Deviation values of the effective-potential',ch10,& 2790 !& ' with respect to the test-set (meV^2/atm):',& 2791 !& ch10,' Energy : ',& 2792 !& mse* (Ha_EV*1000)**2 *factor ,ch10,& 2793 & ' Goal function values of the effective.potential',ch10,& 2794 & ' with respect to the test-set (eV^2/A^2):',ch10,& 2795 & ' Energy : ',& 2796 & (mse)*(HaBohr_eVAng)**2,ch10,& 2797 & ' Forces+Stresses : ',& 2798 & (msef+mses)*(HaBohr_eVAng)**2,ch10,& 2799 & ' Forces : ',& 2800 & msef*(HaBohr_eVAng)**2,ch10,& 2801 & ' Stresses : ',& 2802 & mses*(HaBohr_eVAng)**2,ch10 2803 call wrtout(ab_out,message,'COLL') 2804 call wrtout(std_out,message,'COLL') 2805 2806 2807 !Deallocating 2808 ABI_FREE(sqomega) 2809 ABI_FREE(ucvol) 2810 2811 INQUIRE(FILE='TES_fit_diff_anharmonic_terms_energy.dat',OPENED=file_opened,number=unit_anh) 2812 if(file_opened) close(unit_anh) 2813 2814 2815 end subroutine fit_polynomial_coeff_testEffPot