TABLE OF CONTENTS

ABINIT/m_paw_gaussfit [ Modules ]

[ Top ] [ Modules ]

NAME

  m_paw_gaussfit

FUNCTION

  Module to fit PAW related data to sums of gaussians

COPYRIGHT

  Copyright (C) 2012-2024 ABINIT group (T. Rangel)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .

NOTES

  FOR DEVELOPPERS: in order to preserve the portability of libPAW library,
  please consult ~abinit/src/??_libpaw/libpaw-coding-rules.txt

SOURCE

20 #include "libpaw.h"
21 
22 module m_paw_gaussfit
23 
24  USE_DEFS
25  USE_MSG_HANDLING
26  USE_MPI_WRAPPERS
27  USE_MEMORY_PROFILING
28 
29  use m_paw_numeric, only : paw_splint, paw_spline
30  use m_pawrad,      only : pawrad_type, pawrad_init, pawrad_deducer0, pawrad_free, pawrad_ifromr
31 
32  implicit none
33 
34  private
35 
36  public:: gaussfit_projector !fit non-local projectors to a sum of gaussians
37  public:: gaussfit_main !main routine to fit
38 !Routines related to MPI:
39  private:: gaussfit_mpi_set_weight
40  private:: gaussfit_mpi_remove_item
41  private:: gaussfit_mpi_add_item
42  private:: gaussfit_mpi_assign
43  private:: gaussfit_mpi_main
44  private:: gaussfit_mpi_calc_deviation
45  private:: gaussfit_mpi_swap
46 !Routines related to fitting:
47  private:: gaussfit_fit !fit a function to gaussians
48  private:: gaussfit_calc_deriv_r  !calc. derivatives for real gauss.
49  private:: gaussfit_calc_deriv_c  !calc. derivatives for cplex. gauss. of form 1
50  private:: gaussfit_calc_deriv_c2 !calc. derivatives for cplex. gauss. of form 2
51  private:: gaussfit_calc_deriv_c3 !calc. derivatives for cplex. gauss. of form 3
52  private:: gaussfit_calc_deriv_c4 !calc. derivatives for cplex. gauss. of form 4
53  private:: gaussfit_rlsf !retreined least squares fit
54  private:: gaussfit_chisq_alpha_beta
55 !set parameters for LSF (for 5 different forms of gauss. sums):
56  private:: gaussfit_set_param1
57  private:: gaussfit_set_param2
58  private:: gaussfit_set_param3
59  private:: gaussfit_set_param4
60  private:: gaussfit_set_param5
61  private:: gaussfit_constrains_init !initialize constrains
62  private:: gaussfit_apply_constrains !apply cons. to get new params.

gaussfit_set_param2/gaussfit_param2_findsign [ Functions ]

[ Top ] [ gaussfit_set_param2 ] [ Functions ]

NAME

  gaussfit_param2_findsign

FUNCTION

  Finds the value of y at a given point
  This was taken out of gaussfit_set_param2 to make the code more
  readable.

INPUTS

OUTPUT

SOURCE

1790  subroutine gaussfit_param2_findsign()
1791 
1792 !Arguments -------------------------------
1793 !Local variables-------------------------------
1794  integer::ix,minx
1795  real(dp)::dist,mindist,xx,yy
1796 
1797 ! *********************************************************************
1798 
1799    mindist=rpaw
1800    do ix=1,nx
1801      xx=x(ix)
1802      dist=abs(raux-xx)
1803      if(dist<mindist) then
1804        mindist=dist
1805        minx=ix
1806      end if
1807    end do
1808    yy=y(minx)
1809    sig=yy
1810 
1811  end subroutine gaussfit_param2_findsign
1812 
1813 end subroutine gaussfit_set_param2

m_paw_gaussfit/gaussfit_apply_constrains [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_apply_constrains

FUNCTION

  Apply constrains to get new set of parameters

INPUTS

OUTPUT

SOURCE

2028 subroutine gaussfit_apply_constrains(const,limit,nparam,ioparams)
2029 
2030 !Arguments -------------------------------
2031  integer,intent(in):: nparam
2032  integer,intent(in):: const(nparam)
2033  real(dp),intent(in):: limit(nparam)
2034  real(dp),intent(inout):: ioparams(nparam)
2035 
2036 !Local variables-------------------------------
2037  integer::ii
2038 
2039 ! *********************************************************************
2040 
2041  do ii=1,nparam
2042    if(const(ii)==restricted .or. const(ii)==restricted_and_positive) then
2043      if(ioparams(ii)<limit(ii)) ioparams(ii)=limit(ii)
2044    end if
2045    if(const(ii)==positive .or. const(ii)==restricted_and_positive ) ioparams(ii)=abs(ioparams(ii))
2046 
2047  end do
2048 
2049 end subroutine gaussfit_apply_constrains

m_paw_gaussfit/gaussfit_calc_deriv_c [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_calc_deriv_c

FUNCTION

  Calculate expressions and derivatives for Gaussians fitting.
   The Gaussians expressions are defined in the comments of "gaussfit_main"

INPUTS

OUTPUT

SOURCE

1149 subroutine gaussfit_calc_deriv_c(nparam,nterm,nx,opt,param,x,y_out,&
1150 & deriv) ! optional
1151 
1152 !Arguments -------------------------------
1153  integer,intent(in)::nparam !number of parameters
1154  integer,intent(in)::nterm !number of gaussian expressions
1155  integer,intent(in)::nx     !number of point in the x grid
1156  integer,intent(in)::opt    !option:
1157                             !1) calculate only f(x)
1158                             !2) calculate f(x) and its derivatives
1159  real(dp),intent(in)::param(nparam) !parameters
1160  real(dp),intent(in)::x(nx) !xgrid
1161  real(dp),intent(out)::y_out(nx) !f(x)
1162  real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1163 
1164 !Local variables-------------------------------
1165  integer::iexp,ii
1166  real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1167  real(dp)::alpha4(nterm),alpha5(nterm),alpha6(nterm)
1168  real(dp)::aux1(nx),aux2(nx)
1169  real(dp)::cos1(nx,nterm),cos2(nx,nterm),sin1(nx,nterm),sin2(nx,nterm)
1170  real(dp)::term1(nx,nterm),term2(nx,nterm)
1171 
1172 ! *********************************************************************
1173 
1174 !
1175 !Initialize
1176 !
1177  y_out(:)=0.d0
1178 !
1179 !Get parameters from param array:
1180 !
1181  alpha1(:)=param(1:nterm)
1182  alpha2(:)=param(nterm+1:2*nterm)
1183  alpha3(:)=param(2*nterm+1:3*nterm)
1184  alpha4(:)=param(3*nterm+1:4*nterm)
1185  alpha5(:)=param(4*nterm+1:5*nterm)
1186  alpha6(:)=param(5*nterm+1:6*nterm)
1187 !
1188 !calculate useful quantities
1189 !
1190  do iexp=1,nterm
1191    aux1(:)=-alpha2(iexp)*x(:)**2
1192    term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
1193  end do
1194 !
1195  do iexp=1,nterm
1196    aux1(:)=alpha4(iexp)*x(:)**2
1197    sin1(:,iexp)=sin(aux1(:))
1198 !
1199    aux1(:)=alpha6(iexp)*x(:)**2
1200    sin2(:,iexp)=sin(aux1(:))
1201 !
1202    aux1(:)=alpha4(iexp)*x(:)**2
1203    cos1(:,iexp)=cos(aux1(:))
1204 !
1205    aux1(:)=alpha6(iexp)*x(:)**2
1206    cos2(:,iexp)=cos(aux1(:))
1207  end do
1208 !
1209  do iexp=1,nterm
1210    aux1(:)=alpha3(iexp)*sin1(:,iexp)
1211    aux2(:)=alpha5(iexp)*cos2(:,iexp)
1212    term2(:,iexp)=aux1(:)+aux2(:)
1213    y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp)
1214  end do
1215 !
1216 !Calculate derivatives:
1217 !
1218  if(opt==2) then
1219 !
1220 !  alpha1
1221 !
1222    do iexp=1,nterm
1223      aux1(:)=term1(:,iexp)/alpha1(iexp)
1224      aux2(:)=aux1(:)*term2(:,iexp)
1225      deriv(:,iexp)=aux2(:)
1226    end do
1227 !
1228 !  alpha2
1229 !
1230    do iexp=1,nterm
1231      ii=nterm+iexp
1232      aux1(:)=-term1(:,iexp)*term2(:,iexp)
1233      aux2(:)=aux1(:)*x(:)**2
1234      deriv(:,ii)=aux2(:)
1235    end do
1236 !
1237 !  alpha3
1238 !
1239    do iexp=1,nterm
1240      ii=2*nterm+iexp
1241      aux1(:)=term1(:,iexp)*sin1(:,iexp)
1242      deriv(:,ii)=aux1(:)
1243    end do
1244 !
1245 !  alpha4
1246 !
1247    do iexp=1,nterm
1248      ii=3*nterm+iexp
1249      aux1(:)=term1(:,iexp)*alpha3(iexp)
1250      aux2(:)=cos1(:,iexp)*x(:)**2
1251      deriv(:,ii)=aux2(:)*aux1(:)
1252    end do
1253 !
1254 !  alpha5
1255 !
1256    do iexp=1,nterm
1257      ii=4*nterm+iexp
1258      aux1(:)=term1(:,iexp)*cos2(:,iexp)
1259      deriv(:,ii)=aux1(:)
1260    end do
1261 !
1262 !  alpha6
1263 !
1264    do iexp=1,nterm
1265      ii=5*nterm+iexp
1266      aux1(:)=-term1(:,iexp)*alpha5(iexp)
1267      aux2(:)=sin2(:,iexp)*x(:)**2
1268      deriv(:,ii)=aux1(:)*aux2(:)
1269    end do
1270  end if
1271 
1272 end subroutine gaussfit_calc_deriv_c

m_paw_gaussfit/gaussfit_calc_deriv_c2 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_calc_deriv_c2

FUNCTION

  Calculate expressions and derivatives for Gaussians fitting.
   The Gaussians expressions are defined in the comments of "gaussfit_main"

INPUTS

OUTPUT

SOURCE

1035 subroutine gaussfit_calc_deriv_c2(nparam,nterm,nx,opt,param,x,y_out,&
1036 & deriv) ! optional
1037 
1038 !Arguments -------------------------------
1039  integer,intent(in)::nparam !number of param
1040  integer,intent(in)::nterm  !number of gaussian expressions
1041  integer,intent(in)::nx     !number of point in the x grid
1042  integer,intent(in)::opt    !option:
1043                             !1) calculate only f(x)
1044                             !2) calculate f(x) and its derivatives
1045  real(dp),intent(in)::param(nparam) !parameters
1046  real(dp),intent(in)::x(nx) !xgrid
1047  real(dp),intent(out)::y_out(nx) !f(x)
1048  real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1049 
1050 !Local variables-------------------------------
1051  integer::iexp,ii
1052  real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1053  real(dp)::alpha4(nterm)
1054  real(dp)::term1(nx,nterm)
1055  real(dp)::sin1(nx,nterm),sin2(nx,nterm),cos1(nx,nterm),cos2(nx,nterm)
1056  real(dp)::aux1(nx),aux2(nx)
1057 
1058 ! *********************************************************************
1059 
1060 !
1061 !Initialize
1062 !
1063  y_out(:)=0.d0
1064 !
1065 !Get param from param array:
1066 !
1067  alpha1(:)=param(1:nterm)
1068  alpha2(:)=param(nterm+1:2*nterm)
1069  alpha3(:)=param(2*nterm+1:3*nterm)
1070  alpha4(:)=param(3*nterm+1:4*nterm)
1071 !
1072 !calculate useful quantities
1073 !
1074 !
1075  do iexp=1,nterm
1076    aux1(:)=alpha2(iexp)*x(:)**2
1077    sin1(:,iexp)=sin(aux1(:))
1078 !
1079    aux1(:)=alpha4(iexp)*x(:)**2
1080    sin2(:,iexp)=sin(aux1(:))
1081 !
1082    aux1(:)=alpha2(iexp)*x(:)**2
1083    cos1(:,iexp)=cos(aux1(:))
1084 !
1085    aux1(:)=alpha4(iexp)*x(:)**2
1086    cos2(:,iexp)=cos(aux1(:))
1087  end do
1088 !
1089  do iexp=1,nterm
1090    aux1(:)=alpha1(iexp)*sin1(:,iexp)
1091    aux2(:)=alpha3(iexp)*cos2(:,iexp)
1092    term1(:,iexp)=aux1(:)+aux2(:)
1093    y_out(:)=y_out(:)+term1(:,iexp)
1094  end do
1095 !
1096 !Calculate derivatives:
1097 !
1098  if(opt==2) then
1099 !
1100 !  alpha1
1101 !
1102    do iexp=1,nterm
1103      deriv(:,iexp)=sin1(:,iexp)
1104    end do
1105 !
1106 !  alpha2
1107 !
1108    do iexp=1,nterm
1109      ii=nterm+iexp
1110      aux1(:)=alpha1(iexp)*cos1(:,iexp)*x(:)**2
1111      deriv(:,ii)=aux1(:)
1112    end do
1113 !
1114 !  alpha3
1115 !
1116    do iexp=1,nterm
1117      ii=2*nterm+iexp
1118      deriv(:,ii)=cos2(:,iexp)
1119    end do
1120 !
1121 !  alpha4
1122 !
1123    do iexp=1,nterm
1124      ii=3*nterm+iexp
1125      aux1(:)=-alpha3(iexp)*sin2(:,iexp)*x(:)**2
1126      deriv(:,ii)=aux1(:)
1127    end do
1128  end if
1129 
1130 end subroutine gaussfit_calc_deriv_c2

m_paw_gaussfit/gaussfit_calc_deriv_c3 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_calc_deriv_c3

FUNCTION

  Calculate expressions and derivatives for Gaussians fitting.
   The Gaussians expressions are defined in the comments of "gaussfit_main"

INPUTS

OUTPUT

SOURCE

 941 subroutine gaussfit_calc_deriv_c3(nparam,nterm,nx,opt,param,x,y_out,&
 942 & deriv) ! optional
 943 
 944 !Arguments -------------------------------
 945  integer,intent(in)::nparam !number of parameters
 946  integer,intent(in)::nterm !number of gaussian expressions
 947  integer,intent(in)::nx     !number of point in the x grid
 948  integer,intent(in)::opt    !option:
 949                             !1) calculate only f(x)
 950                             !2) calculate f(x) and its derivatives
 951  real(dp),intent(in)::param(nparam) !parameters
 952  real(dp),intent(in)::x(nx) !xgrid
 953  real(dp),intent(out)::y_out(nx) !f(x)
 954  real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
 955 
 956 !Local variables-------------------------------
 957  integer::iexp,ii
 958  real(dp)::sep
 959  real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
 960  real(dp)::term1(nx,nterm)
 961  real(dp)::sin1(nx,nterm),cos1(nx,nterm)
 962  real(dp)::aux1(nx),aux2(nx)
 963 
 964 ! *********************************************************************
 965 
 966 !
 967 !Initialize
 968 !
 969  y_out(:)=0.d0
 970 !
 971  sep=1.2d0
 972 !
 973 !Get param from param array:
 974 !
 975  alpha1(:)=param(1:nterm)
 976  alpha2(:)=param(nterm+1:2*nterm)
 977 !
 978  do ii=1,nterm
 979    alpha3(ii)=sep**(ii)
 980  end do
 981 !
 982 !calculate useful quantities
 983 !
 984  do iexp=1,nterm
 985    aux1(:)=alpha3(iexp)*x(:)**2
 986 !
 987    sin1(:,iexp)=sin(aux1(:))
 988    cos1(:,iexp)=cos(aux1(:))
 989  end do
 990 !
 991  do iexp=1,nterm
 992    aux1(:)=alpha1(iexp)*sin1(:,iexp)
 993    aux2(:)=alpha2(iexp)*cos1(:,iexp)
 994    term1(:,iexp)=aux1(:)+aux2(:)
 995    y_out(:)=y_out(:)+term1(:,iexp)
 996  end do
 997 !
 998 !Calculate derivatives:
 999 !
1000  if(opt==2) then
1001 !
1002 !  alpha1
1003 !
1004    do iexp=1,nterm
1005      deriv(:,iexp)=sin1(:,iexp)
1006    end do
1007 !
1008 !  alpha2
1009 !
1010    do iexp=1,nterm
1011      ii=nterm+iexp
1012      deriv(:,ii)=cos1(:,iexp)
1013    end do
1014  end if
1015 
1016 end subroutine gaussfit_calc_deriv_c3

m_paw_gaussfit/gaussfit_calc_deriv_c4 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_calc_deriv_c4

FUNCTION

  Calculate expressions and derivatives for Gaussians fitting.
   The Gaussians expressions are defined in the comments of "gaussfit_main"

INPUTS

OUTPUT

SOURCE

1291 subroutine gaussfit_calc_deriv_c4(nparam,nterm,nx,opt,param,x,y_out,&
1292 & deriv) ! optional
1293 
1294 !Arguments -------------------------------
1295  integer,intent(in)::nparam !number of parameters
1296  integer,intent(in)::nterm !number of gaussian expressions
1297  integer,intent(in)::nx     !number of point in the x grid
1298  integer,intent(in)::opt    !option:
1299                             !1) calculate only f(x)
1300                             !2) calculate f(x) and its derivatives
1301  real(dp),intent(in)::param(nparam) !parameters
1302  real(dp),intent(in)::x(nx) !xgrid
1303  real(dp),intent(out)::y_out(nx) !f(x)
1304  real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
1305 
1306 !Local variables-------------------------------
1307  integer::iexp,ii
1308  real(dp)::raux,sep
1309  real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
1310  real(dp)::alpha4(nterm),alpha5(nterm)
1311  real(dp)::aux1(nx),aux2(nx)
1312  real(dp)::cos1(nx,nterm),sin1(nx,nterm)
1313  real(dp)::term1(nx,nterm),term2(nx,nterm)
1314 
1315 ! *********************************************************************
1316 
1317 !
1318 !Initialize
1319 !
1320  sep=1.1d0
1321  y_out(:)=0.d0
1322 
1323 !Get parameters from param array:
1324 !
1325  alpha1(:)=param(1:nterm)
1326  alpha2(:)=param(nterm+1:2*nterm)
1327  alpha3(:)=param(2*nterm+1:3*nterm)
1328  alpha4(:)=param(3*nterm+1:4*nterm)
1329 !
1330 !
1331  raux=(2.d0*pi)/real(nterm,dp)
1332  do ii=1,nterm
1333    alpha5(ii)=sep**(ii)
1334 !  alpha5(ii)=raux*real(ii-1,dp)
1335  end do
1336 !
1337 !calculate useful quantities
1338 !
1339  do iexp=1,nterm
1340    aux1(:)=-alpha2(iexp)*x(:)**2
1341    term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
1342  end do
1343 !
1344  do iexp=1,nterm
1345    aux1(:)=alpha5(iexp)*x(:)**2
1346 !
1347    sin1(:,iexp)=sin(aux1(:))
1348    cos1(:,iexp)=cos(aux1(:))
1349  end do
1350 !
1351  do iexp=1,nterm
1352    aux1(:)=alpha3(iexp)*sin1(:,iexp)
1353    aux2(:)=alpha4(iexp)*cos1(:,iexp)
1354    term2(:,iexp)=aux1(:)+aux2(:)
1355    y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp)
1356  end do
1357 !
1358 !Calculate derivatives:
1359 !
1360  if(opt==2) then
1361 !
1362 !  alpha1
1363 !
1364    do iexp=1,nterm
1365      aux1(:)=term1(:,iexp)/alpha1(iexp)
1366      aux2(:)=aux1(:)*term2(:,iexp)
1367      deriv(:,iexp)=aux2(:)
1368    end do
1369 !
1370 !  alpha2
1371 !
1372    do iexp=1,nterm
1373      ii=nterm+iexp
1374      aux1(:)=-term1(:,iexp)*term2(:,iexp)
1375      aux2(:)=aux1(:)*x(:)**2
1376      deriv(:,ii)=aux2(:)
1377    end do
1378 !
1379 !  alpha3
1380 !
1381    do iexp=1,nterm
1382      ii=2*nterm+iexp
1383      aux1(:)=term1(:,iexp)*sin1(:,iexp)
1384      deriv(:,ii)=aux1(:)
1385    end do
1386 !
1387 !  alpha4
1388 !
1389    do iexp=1,nterm
1390      ii=3*nterm+iexp
1391      aux1(:)=term1(:,iexp)*cos1(:,iexp)
1392      deriv(:,ii)=aux1(:)
1393    end do
1394  end if
1395 
1396 end subroutine gaussfit_calc_deriv_c4

m_paw_gaussfit/gaussfit_calc_deriv_r [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_calc_deriv_r

FUNCTION

  Calculate derivatives for Gaussians
   Only relevant for fitting Gaussians algorithm.
   The Gaussians expressions are defined in the comments of "gaussfit_main"

INPUTS

OUTPUT

SOURCE

836 subroutine gaussfit_calc_deriv_r(nterm,nparam,nx,opt,param,x,y_out,&
837 & deriv) ! optional
838 
839 !Arguments -------------------------------
840  integer,intent(in)::nx     !number of point in the x grid
841  integer,intent(in)::nparam !number of parameters
842  integer,intent(in)::nterm !number of gaussian expressions
843  integer,intent(in)::opt    !option:
844                             !1) calculate only f(x)
845                             !2) calculate f(x) and its derivatives
846  real(dp),intent(in)::param(nparam) !parameters
847  real(dp),intent(in)::x(nx) !xgrid
848  real(dp),intent(out)::y_out(nx) !f(x)
849  real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives
850 
851 !Local variables-------------------------------
852  integer::iexp,ii
853  real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm)
854  real(dp)::term1(nx,nterm)
855  real(dp)::aux1(nx)
856  !real(dp)::step
857 
858 ! *********************************************************************
859 
860 !
861 !Initialize
862 !
863  y_out(:)=0.d0
864 !
865 !Get parameters from parameters array:
866 !
867  alpha1(:)=param(1:nterm)
868  alpha2(:)=param(nterm+1:2*nterm)
869  alpha3(:)=param(2*nterm+1:3*nterm)
870 !
871 !alpha3
872 !set to constant values of x
873 !step=rpaw/real(nterm,dp)
874 !do ii=1,nterm
875 !raux=step*real(ii,dp)
876 !alpha3(ii)=raux
877 !end do
878 !
879 !
880 !
881 !calculate useful quantities
882 !
883  do iexp=1,nterm
884    aux1(:)=-alpha2(iexp)*(x(:)-alpha3(iexp))**2
885    term1(:,iexp)=alpha1(iexp)*exp(aux1(:))
886  end do
887 !
888  do iexp=1,nterm
889    y_out(:)=y_out(:)+term1(:,iexp)
890  end do
891 !
892 !Calculate derivatives:
893 !
894  if(opt==2) then
895 !
896 !  alpha1
897 !
898    do iexp=1,nterm
899      aux1(:)=term1(:,iexp)/alpha1(iexp)
900      deriv(:,iexp)=aux1(:)
901    end do
902 !
903 !  alpha2
904 !
905    do iexp=1,nterm
906      ii=nterm+iexp
907      aux1(:)=-term1(:,iexp)*(x(:)-alpha3(iexp))
908      deriv(:,ii)=aux1(:)
909      deriv(:,ii)=0.1d0
910    end do
911 !
912 !  alpha3
913 !
914    do iexp=1,nterm
915      ii=2*nterm+iexp
916      aux1(:)=term1(:,iexp)*2.d0*alpha2(iexp)
917      aux1(:)=aux1(:)*(x(:)-alpha3(iexp))
918      deriv(:,ii)=aux1(:)
919    end do
920  end if
921 
922 end subroutine gaussfit_calc_deriv_r

m_paw_gaussfit/gaussfit_chisq_alpha_beta [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_chisq_alpha_beta

FUNCTION

  Finds chisq, alpha and beta parameters for LSF using the Levenberg-Marquardt algorithm.

COPYRIGHT

  Copyright (C) 2011-2024 ABINIT group (T. Rangel)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .
  The original Levemberg Marquardt routines were written by Armando Sole
  These were modified for the ARPUS spectra in the BigDFT code by A. Mirone.
  These were re-writen in Fortran and further modified in ABINIT for our particular needs.

INPUTS

OUTPUT

SOURCE

1599 subroutine gaussfit_chisq_alpha_beta(alpha,beta,chisq,&
1600 & nparam,nterm,nx,option,parameters,x,y)
1601 
1602 !Arguments -------------------------------
1603  integer,intent(in)::nparam,nterm,nx
1604  integer,intent(in)::option
1605  real(dp),intent(out)::chisq
1606  real(dp),intent(in)::parameters(nparam),x(nx),y(nx)
1607  real(dp),intent(out)::alpha(nparam,nparam),beta(nparam)
1608 
1609 !Local variables-------------------------------
1610  integer::ii,jj,kk
1611  real(dp)::deltax,help1
1612  !arrays
1613  real(dp)::deltay(nx),deriv(nx,nparam),derivi(nx)
1614  real(dp)::yfit(nx)
1615  real(dp)::help0(nx),help2(nx),help3(nparam)
1616 
1617 ! *********************************************************************
1618 
1619  deltax=x(2)-x(1) !we assume a linear grid
1620 !
1621  if(option==1) then
1622    call gaussfit_calc_deriv_c2(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1623  elseif(option==2) then
1624    call gaussfit_calc_deriv_c(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1625  elseif(option==3) then
1626    call gaussfit_calc_deriv_c3(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1627  elseif(option==4) then
1628    call gaussfit_calc_deriv_c4(nparam,nterm,nx,2,parameters,x,yfit,deriv)
1629  end if
1630  deltay=y-yfit
1631  help0=deltay
1632 !
1633  do ii=1,nparam
1634    derivi(:)=deriv(:,ii)
1635    help1=0.d0
1636    do jj=1,nx
1637      help1=help1+help0(jj)*derivi(jj)
1638    end do
1639    beta(ii)=help1
1640 !  help1 = innerproduct(deriv,weight*derivi)
1641 !  below I use help3 instead for the array dimenstions
1642    help3=0.d0
1643    do kk=1,nparam
1644      do jj=1,nx
1645        help3(kk)=help3(kk)+deriv(jj,kk)*derivi(jj)
1646      end do
1647    end do
1648 !  !
1649    alpha(:,ii)=help3(:)
1650  end do
1651 !
1652  help2(:)=help0(:)*deltay(:)
1653  chisq=sum(help2)
1654  chisq=chisq*deltax
1655 
1656 end subroutine gaussfit_chisq_alpha_beta

m_paw_gaussfit/gaussfit_constrains_init [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_constrains_init

FUNCTION

  Initialise constrains for LSF.
  It will constrain the Gaussians width
  It will also constraint the Delta use in the LSF algorithm, to jump slowly at each step.

INPUTS

OUTPUT

SOURCE

1959 subroutine gaussfit_constrains_init(cons1,cons2,limit,nparam,nterm,nx,option,rpaw,y)
1960 
1961 !Arguments -------------------------------
1962  integer,intent(in)::nterm,option,nparam,nx
1963  integer,intent(out)::cons2(nparam)
1964  real(dp),intent(in)::rpaw
1965  real(dp),intent(in)::y(nx)
1966  real(dp),intent(out)::cons1(nparam),limit(nparam)
1967 
1968 !Local variables-------------------------------
1969  integer :: ix
1970  real(dp)::rc,a1,mm,raux
1971 
1972 ! *********************************************************************
1973 
1974 !
1975 !DEFAULT: no weight
1976  cons1(:)=1.d0
1977  limit(:)=0.d0
1978  cons2=1
1979 !
1980  if(option==4) then
1981      cons1(1:nterm)=0.2d0
1982      cons1(nterm+1:nterm*2)=0.3d0
1983  end if
1984 !
1985  if(option==4) then
1986 !  parameters
1987 
1988 !  a1=maxval( y(:),nx)/real(nterm,dp)
1989 !  maxval gives problems in some architectures:
1990    a1=-9999999
1991    do ix=1,nx
1992      if(a1<y(ix)) a1=y(ix)
1993    end do
1994    a1=a1/real(nterm,dp)
1995 
1996    mm=0.01
1997 !
1998    rc=1.7d0*rpaw
1999    raux=log(a1/mm)/rc**2
2000 !
2001 !  Constraint exponential of gaussians to be positive (multiplied by -1),
2002 !  so that it decays to zero.
2003 !  Constraint as well its value, so that it does not get too big
2004 !  and it decays soon,
2005 !  This prevents that a gaussian grows at a very large x value.
2006    cons2(nterm+1:nterm*2)=restricted_and_positive
2007    limit(nterm+1:nterm*2)=raux
2008  end if
2009 
2010 end subroutine gaussfit_constrains_init

m_paw_gaussfit/gaussfit_fit [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_fit

FUNCTION

  Fits a given input function f(r) to a sum of gaussians

INPUTS

  chisq= It measures how good is the fitting.
   It is defined here as sum_i^N_i f(x_i)-y(x_i)/N_i
  constrains= constraints for Gaussians
  limit(nparam)= it limits the widths of Gaussians
  maxiter=maximum number of iterations
  nparam= number of parameters (a constant times nterm)
         if(option==1)nparam=nterm*4
         if(option==2)nparam=nterm*6
         if(option==3)nparam=nterm*2
         if(option==4)nparam=nterm*4
  nterm= number of Gaussians
  nx=number of points along the x axis
  option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
         2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
         3 fit to a1 cos (k x^2) + a2 sin (k x^2)
         4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
  outfile= filename for output (only written if verbosity>1)
  verbosity= controls output volume
  weight(nparam)= weights for the fitting procedure
  x(nx)= points along the x axis
  y(nx)= function to be fitted
  rpaw ,optional= paw radius

OUTPUT

 y_out(nx)= fitted function

SIDE EFFECTS

 if(verbosity>1) output files are written with y(x) and y_out(x)

SOURCE

759 subroutine gaussfit_fit(chisq,constrains,&
760 & limit,maxiter,nparam,nterm,nx,option,outfile,param,&
761 & verbosity,weight,x,y,y_out)
762 
763 !Arguments ------------------------------------
764  integer, intent(in) :: maxiter,nparam
765  integer,intent(in)  :: nterm,nx,option,verbosity
766  !real(dp),optional,intent(in)::rpaw
767  !arrays
768  integer,intent(in)::constrains(nparam)
769  real(dp),intent(in)::limit(nparam),weight(nparam)
770  real(dp),intent(in)::x(nx),y(nx)
771  real(dp),intent(inout)::param(nparam)
772  real(dp),intent(out)::chisq,y_out(nx)
773  character(80),intent(in)::outfile
774 
775 !Local variables ------------------------------
776  integer, parameter :: wfn_unit=1007
777  integer::ix
778  real(dp)::rerror
779  real(dp),allocatable::sy(:)
780 
781 ! *************************************************************************
782 
783  LIBPAW_ALLOCATE(sy,(nx))
784 
785  sy(:)=1.0d0
786 
787  call gaussfit_rlsf(&
788 & chisq,constrains,limit,maxiter,&
789 & nterm,nparam,nx,option,param(1:nparam),&
790 & verbosity,weight,x,y)
791 !
792  if(verbosity>=1) then
793    if(option==1) then
794      call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,param,x,y_out)
795    elseif(option==2) then
796      call gaussfit_calc_deriv_c(nparam,nterm,nx,1,param,x,y_out)
797    elseif(option==3) then
798      call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,param,x,y_out)
799    elseif(option==4) then
800      call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,param,x,y_out)
801    end if
802 !
803 !
804    open(wfn_unit,file=outfile,form='formatted',status='unknown')
805 !  per_error=0.d0
806    do ix=1, nx
807      rerror=abs(y(ix)-y_out(ix))
808      write(wfn_unit,'(6(e20.12,1x))')x(ix),y(ix),y_out(ix),rerror
809    end do
810    close(wfn_unit)
811 
812  end if
813 
814  LIBPAW_DEALLOCATE(sy)
815 
816 end subroutine gaussfit_fit

m_paw_gaussfit/gaussfit_main [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_main

FUNCTION

  Fits a given input function f(r) to a sum of gaussians

INPUTS

  nterm_bounds= sets the minimum and maximum number of terms to be taken into account.
  mparam= maximum number of parameters
         if(option==1)mparam=nterm_bounds(2)*4
         if(option==2)mparam=nterm_bounds(2)*6
         if(option==3)mparam=nterm_bounds(2)*2
         if(option==4)mparam=nterm_bounds(2)*4
  nr= number of real space points
  pawrad= pawrad type
  option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
         2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
         3 fit to a1 cos (k x^2) + a2 sin (k x^2)
         4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
  Given these definitions, the number of complex gaussians are:
    ngauss=4*nterm (for option=1,2)
    ngauss=2*nterm (for option=3,4)
  outfile= filename to write out fitted functions.
  rpaw=paw radius
  y(nr)= function to fit

OUTPUT

 nparam_out= number of parameters found.
 param_out(nparam_out)= parameters (coefficients and factors of complex gaussians).

SOURCE

106  subroutine gaussfit_main(mparam,nparam_out,nterm_bounds,nr,&
107 &           param_out,pawrad,option,outfile,rpaw,y,comm_mpi)
108 
109 !Arguments ------------------------------------
110  integer,intent(in)::mparam,nr,nterm_bounds(2)
111  integer,intent(in)::option
112  integer,intent(in),optional:: comm_mpi
113  real(dp),intent(in)::rpaw
114  real(dp),intent(inout)::y(nr)
115  character(80),intent(in)::outfile
116  type(pawrad_type),intent(in) :: pawrad
117  integer,intent(out)::nparam_out
118  real(dp),intent(out)::param_out(mparam)
119 
120 !Local variables-------------------------------
121 !scalars
122  logical,parameter::modify_y=.false. !used only for plotting purposes
123  integer :: ichisq,ierr,ii,jj
124  integer :: master,me
125  integer :: maxiter,minterm,my_chisq_size,counts_all
126  integer :: ngauss,nparam,nproc,nterm
127  integer :: verbosity
128  real(dp) :: chisq,chisq_min
129 !real(dp) :: T1,T2 !uncomment for timming
130  !arrays
131  integer::constrains(mparam)
132  integer::proc_dist(nterm_bounds(1):nterm_bounds(2))
133  integer,allocatable::counts(:),disp(:),map_nterm(:)
134  real(dp)::limit(mparam),param_tmp(mparam),weight(mparam)
135  real(dp),allocatable::chisq_array(:),recv_buf(:),send_buf(:)
136  real(dp),allocatable::y_out(:)
137  character(len=500) :: msg
138 
139 ! *************************************************************************
140 
141 !initialize variables
142  maxiter=200
143 !
144 !initialize mpi quantities:
145  master=0; me=0; nproc=1; proc_dist=1;
146  if(present(comm_mpi)) then
147    me=xmpi_comm_rank(comm_mpi)
148    nproc=xmpi_comm_size(comm_mpi)
149  end if
150  if(nproc>1) then
151 !  Find distribution (master)
152    if(me==master) then
153      call gaussfit_mpi_main(nproc,nterm_bounds,proc_dist)
154    end if
155 !  send distribution to all processors
156    call xmpi_bcast(proc_dist(nterm_bounds(1):nterm_bounds(2)),&
157 &    master,comm_mpi,ierr)
158  end if
159 !Set size of chisq treated by each proc
160  my_chisq_size=nterm_bounds(2)-nterm_bounds(1)+1 !all terms
161  if(nproc>1 .and. .not. me==master) then
162    do nterm=nterm_bounds(1),nterm_bounds(2)
163      if(.not. proc_dist(nterm)==me+1) cycle
164      my_chisq_size=my_chisq_size+1
165    end do
166  end if
167 
168 !
169 !Allocate objects
170 !
171  LIBPAW_ALLOCATE(y_out,(nr))
172  LIBPAW_ALLOCATE(chisq_array,(my_chisq_size))
173  if(master==me ) then
174    LIBPAW_BOUND1_ALLOCATE(map_nterm,BOUNDS(nterm_bounds(1),nterm_bounds(2)))
175    jj=1
176    do ii=1,nproc
177      do nterm=nterm_bounds(1),nterm_bounds(2)
178        if(proc_dist(nterm)==ii) then
179          map_nterm(nterm)=jj
180          jj=jj+1
181        end if
182      end do
183    end do
184  end if
185 !
186 !fill with zeros
187 !
188  chisq_array=zero
189  nparam_out=0
190  y_out=0.d0
191  verbosity=1 !print the minimum at the terminal
192 !
193  ichisq=0
194  do nterm=nterm_bounds(1),nterm_bounds(2)
195 !  mpi distribution
196    if(.not. proc_dist(nterm)==me+1) cycle
197    ichisq=ichisq+1
198 
199 !   call CPU_TIME(T1)
200 
201    if(option==1) nparam=nterm*4
202    if(option==2) nparam=nterm*6
203    if(option==3) nparam=nterm*2
204    if(option==4) nparam=nterm*4
205 !  set initial guess
206 !   if(option==1) then
207 !     call gaussfit_set_param3(nterm,nparam,param_tmp(1:nparam),sep(ii))
208 !   elseif(option==2) then
209 !     call gaussfit_set_param1(nterm,nparam,nr,&
210 !&     param_tmp(1:nparam),sep(ii),pawrad%rad(1:nr),y)
211    if(option==3) then
212      call gaussfit_set_param4(nparam,param_tmp(1:nparam))
213    elseif(option==4) then
214      call gaussfit_set_param5(nterm,nparam,nr,&
215 &     param_tmp(1:nparam),rpaw,y)
216    end if
217 !
218 !
219    call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),&
220 &   limit(1:nparam),nparam,nterm,nr,option,rpaw,y)
221 !
222    call gaussfit_fit(chisq,constrains(1:nparam),&
223 &   limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,&
224 &   verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out)
225 
226 !  if there was an error, set chisq to very high
227 !  if there is an NaN, for instance:
228    if(abs(chisq+1.d0)<tol8) chisq=99999
229    if(abs(chisq)==chisq*chisq) chisq=99999
230    if(chisq .ne. chisq) chisq=99999
231 !
232    chisq_array(ichisq)=chisq
233 
234 !   call CPU_TIME(T2)
235 !   print *, 'Time for fit ', T2-T1, 'seconds.'
236 
237  end do !nterm
238 
239 !mpicast results
240 !send distribution to master
241  if(nproc>1) then
242 !  Prepare communications:
243    LIBPAW_ALLOCATE(counts,(nproc))
244    counts(:)=0
245    do nterm=nterm_bounds(1),nterm_bounds(2)
246      counts(proc_dist(nterm))=counts(proc_dist(nterm))+1
247    end do
248    counts_all=sum(counts)
249    LIBPAW_ALLOCATE(send_buf,(counts(me+1)))
250    send_buf(:)=chisq_array(1:counts(me+1))
251    if(me==master) then
252      LIBPAW_ALLOCATE(recv_buf,(counts_all))
253    else
254      LIBPAW_ALLOCATE(recv_buf,(1))
255    end if
256    LIBPAW_ALLOCATE(disp,(nproc))
257    disp(1)=0
258    do ii=2,nproc
259      disp(ii)=disp(ii-1)+counts(ii-1)
260    end do
261 !  communicate all info to master
262    call xmpi_gatherv(send_buf,counts(me+1),recv_buf,&
263 &    counts,disp,master,comm_mpi,ierr)
264 !  fill in chisq_array with all received info:
265    if(master==me) then
266      do ii=1,counts_all
267        chisq_array(ii)=recv_buf(ii)
268      end do
269    end if
270 !  Deallocate MPI arrays:
271    LIBPAW_DEALLOCATE(recv_buf)
272    LIBPAW_DEALLOCATE(counts)
273    LIBPAW_DEALLOCATE(disp)
274    LIBPAW_DEALLOCATE(send_buf)
275  end if
276 
277 !Print out info:
278  if(me==master) then
279    write(msg,'(3a)')'Preliminary results (with only 200 iter.):',ch10,'   ngauss    chisq'
280    call wrtout(std_out,msg,'COLL')
281    do nterm=nterm_bounds(1),nterm_bounds(2)
282      if(option==1) ngauss=nterm*4
283      if(option==2) ngauss=nterm*4
284      if(option==3) ngauss=nterm*2
285      if(option==4) ngauss=nterm*2
286      write(msg,'(i4,2x,e13.6,1x)')ngauss,&
287 &     chisq_array(map_nterm(nterm))
288      call wrtout(std_out,msg,'COLL')
289    end do
290  end if
291 
292 !get minterm for best accuracy:
293  if(me==master) then
294    chisq_min=999999999
295    do nterm=nterm_bounds(1),nterm_bounds(2)
296      if(chisq_array(map_nterm(nterm))<chisq_min) then
297        chisq_min=chisq_array(map_nterm(nterm))
298        minterm=nterm
299      end if
300    end do
301 
302 !run again with the best parameters
303    nterm=minterm
304    if(option==1)nparam=4*nterm
305    if(option==2)nparam=6*nterm
306    if(option==3)nparam=2*nterm
307    if(option==4)nparam=4*nterm
308 
309 !  set initial guess
310    if (option==3) then
311      call gaussfit_set_param4(nparam,param_tmp(1:nparam))
312    elseif (option==4) then
313      call gaussfit_set_param5(nterm,nparam,nr,&
314 &     param_tmp(1:nparam),rpaw,y)
315    end if
316 !
317    call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),&
318 &   limit(1:nparam),nparam,nterm,nr,option,rpaw,y)
319 !
320    verbosity=1;
321    maxiter=1000 !this time we do more iterations
322    call gaussfit_fit(chisq,constrains(1:nparam),&
323 &   limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,&
324 &   verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out)
325 
326 !  Write out best solution
327    write(msg,'(3a)')"Best solution (with more iterations):",ch10,"   ngauss    chisq"
328    call wrtout(std_out,msg,'COLL')
329    if(option==1) ngauss=nterm*4
330    if(option==2) ngauss=nterm*4
331    if(option==3) ngauss=nterm*2
332    if(option==4) ngauss=nterm*2
333    write(msg,'(i4,2x,e13.6,1x)')ngauss,chisq
334    call wrtout(std_out,msg,'COLL')
335  end if
336 
337 !Fill output variables
338 !and communicate results to all procs:
339  if(option==1) then
340    nparam_out=minterm*4
341  elseif (option==2) then
342    nparam_out=minterm*6
343  elseif (option==3) then
344    nparam_out=minterm*2
345  elseif (option==4) then
346    nparam_out=minterm*4
347  end if
348 
349 !communicate
350  if(nproc>1) then
351    call xmpi_bcast(nparam_out,master,comm_mpi,ierr)
352    call xmpi_bcast(param_tmp(1:nparam_out),master,comm_mpi,ierr)
353  end if !nproc>1
354 
355  param_out(:)=param_tmp
356 
357  if(modify_y) then
358 !   call xmpi_scatterv(y_out,nr,mpi_displ,y_out,nr,0,comm_mpi,ierr)
359    y=y_out !at output modify y for the fitted y
360  end if
361 
362  LIBPAW_DEALLOCATE(y_out)
363  LIBPAW_DEALLOCATE(chisq_array)
364  if(me==master) then
365    LIBPAW_DEALLOCATE(map_nterm)
366  end if
367 
368 end subroutine gaussfit_main

m_paw_gaussfit/gaussfit_mpi_add_item [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_add_item

FUNCTION

INPUTS

OUTPUT

SOURCE

462  subroutine gaussfit_mpi_add_item(iterm,pload)
463 
464 !Arguments ------------------------------------
465  integer,intent(in)::iterm
466  integer,intent(inout)::pload
467 
468 !Local variables ------------------------------
469  integer:: f_i
470 
471 !************************************************************************
472 
473  call gaussfit_mpi_set_weight(f_i,iterm)
474  pload=pload+f_i
475 
476  end subroutine gaussfit_mpi_add_item

m_paw_gaussfit/gaussfit_mpi_assign [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_assign

FUNCTION

  Set task to a processor

INPUTS

OUTPUT

SOURCE

592  subroutine gaussfit_mpi_assign(iterm,nproc,nterm_bounds,&
593 & proc_dist,proc_load)
594 
595 !Arguments ------------------------------------
596  integer,intent(in)::iterm,nproc,nterm_bounds(2)
597  integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc)
598 
599 !Local variables ------------------------------
600  integer:: iproc,jproc,dev
601  integer:: deviation(nproc),mindev
602  character(len=100) :: msg
603 
604 !************************************************************************
605 
606  do iproc=1,nproc
607   !add this term to iproc
608   call gaussfit_mpi_add_item(iterm,proc_load(iproc))
609   !calculate the deviation for this configuration
610   call gaussfit_mpi_calc_deviation(dev,nproc,proc_load)
611   deviation(iproc)=dev
612   !remove this term from iproc
613   call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
614  end do
615 
616 !assign to jproc, the proc with minimal deviation above:
617  jproc=-1; mindev=999999999
618  do iproc=1,nproc
619    if(deviation(iproc)<mindev) then
620      mindev=deviation(iproc)
621      jproc=iproc
622    end if
623  end do
624  if(jproc==-1) then
625 !  One should not get here!
626    msg = 'error in accomodate_mpi'
627    LIBPAW_BUG(msg)
628  end if
629 
630 !assign this term for jproc
631  proc_dist(iterm)=jproc
632  call gaussfit_mpi_add_item(iterm,proc_load(jproc))
633 
634  end subroutine gaussfit_mpi_assign

m_paw_gaussfit/gaussfit_mpi_calc_deviation [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_calc_deviation

FUNCTION

INPUTS

OUTPUT

SOURCE

493  subroutine gaussfit_mpi_calc_deviation(deviation,nproc,proc_load)
494 
495 !Arguments ------------------------------------
496  integer,intent(in)::nproc
497  integer,intent(in)::proc_load(nproc)
498  integer,intent(out)::deviation
499 
500 !Local variables ------------------------------
501  integer:: jproc,kproc,kload,jload
502 
503 !************************************************************************
504 
505 ! deviation=0
506 ! do jproc=1,nproc
507 !   deviation=deviation+abs(proc_load(jproc)-ideal)
508 ! end do
509 
510  deviation=0
511  do jproc=1,nproc
512    jload=proc_load(jproc)
513    do kproc=1,jproc
514      kload=proc_load(kproc)
515      deviation=deviation+abs(kload-jload)
516    end do
517  end do
518 
519  end subroutine gaussfit_mpi_calc_deviation

m_paw_gaussfit/gaussfit_mpi_main [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_main

FUNCTION

  Set charge for each processor

INPUTS

OUTPUT

SOURCE

652  subroutine gaussfit_mpi_main(nproc,nterm_bounds,proc_dist)
653 
654 !Arguments ------------------------------------
655  integer,intent(in)::nproc,nterm_bounds(2)
656  integer,intent(out)::proc_dist(nterm_bounds(1):nterm_bounds(2))
657 
658 !Local variables ------------------------------
659  integer:: dev1,dev2,ii
660  integer:: iproc,iterm,jterm,ngauss,weight
661  integer:: proc_load(nproc)
662  character(len=500) :: msg
663 
664 !************************************************************************
665 
666  proc_load=0; proc_dist=0 !initializations
667 
668 !1) get a first-trial distribution:
669  do iterm=nterm_bounds(2),nterm_bounds(1),-1
670      call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load)
671  end do
672 
673 !Do the following 20 times
674  do ii=1,20
675 !  Calculate initial state
676    call gaussfit_mpi_calc_deviation(dev1,nproc,proc_load)
677 !  Try to swap tasks between two processors:
678    do iterm=nterm_bounds(2),nterm_bounds(1),-1
679      do jterm=nterm_bounds(2),nterm_bounds(1),-1
680        call gaussfit_mpi_swap(iterm,jterm,&
681 &       nproc,nterm_bounds,proc_dist,proc_load)
682      end do
683    end do
684 !!  Try to reassign tasks to different processors
685    do iterm=nterm_bounds(2),nterm_bounds(1),-1
686      iproc=proc_dist(iterm)
687 !    Remove this job from this node
688      call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
689 !    Accomodate this again:
690      call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load)
691    end do
692 !  Calculate final state
693    call gaussfit_mpi_calc_deviation(dev2,nproc,proc_load)
694 !  If final state equals initial state, exit:
695    if(dev2 == dev1) exit
696 !   write(*,'(a)')'Deviation: ',dev2
697  end do
698 
699 !Write down distribution:
700  write(msg,'(3a)') 'MPI distribution',ch10,'N. gauss, iproc, weight '
701  call wrtout(std_out,msg,'COLL')
702  do iterm=nterm_bounds(2),nterm_bounds(1),-1
703      ngauss=iterm*2
704      call gaussfit_mpi_set_weight(weight,iterm)
705      write(msg,'(3(i4,1x))') ngauss,proc_dist(iterm),weight
706      call wrtout(std_out,msg,'COLL')
707  end do
708  write(msg,'(a)') 'Load per processor: '
709  call wrtout(std_out,msg,'COLL')
710  do iproc=1,nproc
711    write(msg,'(i5,1x,i10)') iproc,proc_load(iproc)
712    call wrtout(std_out,msg,'COLL')
713  end do
714 
715  end subroutine gaussfit_mpi_main

m_paw_gaussfit/gaussfit_mpi_remove_item [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_remove_item

FUNCTION

INPUTS

OUTPUT

SOURCE

431  subroutine gaussfit_mpi_remove_item(iterm,pload)
432 
433 !Arguments ------------------------------------
434  integer,intent(in)::iterm
435  integer,intent(inout)::pload
436 
437 !Local variables ------------------------------
438  integer:: f_i
439 
440 !***********************************************************************
441 
442  call gaussfit_mpi_set_weight(f_i,iterm)
443  pload=pload-f_i
444 
445  end subroutine gaussfit_mpi_remove_item

m_paw_gaussfit/gaussfit_mpi_set_weight [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_set_weight

FUNCTION

 It sets a weight to the number of
 Gaussians used.
 This was calculated by measuring the time
 it takes to fit a projector with different
 number of gaussians.

INPUTS

OUTPUT

SOURCE

390  subroutine gaussfit_mpi_set_weight(f,x)
391 
392 !Arguments ------------------------------------
393  integer,intent(in)::x
394  integer,intent(out)::f
395 
396 !Local variables ------------------------------
397  real(dp)::a,b,c,d,ff,xx
398 
399 !************************************************************************
400 
401  !The following parameters were obtained
402  !from the time (in seconds)
403  !it takes to fit a given projector
404  !using from 1 to 100 gaussians.
405  a               = -0.374137d0      !  +/- 2.013        (538.1%)
406  b               = 0.207854d0       !  +/- 0.3385       (162.8%)
407  c               = 0.0266371d0      !  +/- 0.01534      (57.59%)
408  d               = 0.000152476d0    !  +/- 0.0001978    (129.7%)
409 
410  xx=real(x,dp)
411  ff=a+b*xx+c*xx**2+d*xx**3
412  f=max(1,ceiling(ff))
413 
414  end subroutine gaussfit_mpi_set_weight

m_paw_gaussfit/gaussfit_mpi_swap [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_mpi_swap

FUNCTION

INPUTS

OUTPUT

SOURCE

536  subroutine gaussfit_mpi_swap(iterm,jterm,&
537 &           nproc,nterm_bounds,proc_dist,proc_load)
538 
539 !Arguments ------------------------------------
540  integer,intent(in)::iterm,jterm,nproc,nterm_bounds(2)
541  integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc)
542 
543 !Local variables ------------------------------
544  integer:: deviation1,deviation2
545  integer:: iproc,jproc
546 
547 !************************************************************************
548 
549 !Calculate initial state
550  call gaussfit_mpi_calc_deviation(deviation1,nproc,proc_load)
551  iproc=proc_dist(iterm)
552  jproc=proc_dist(jterm)
553 !Swap terms:
554  call gaussfit_mpi_add_item(jterm,proc_load(iproc))
555  call gaussfit_mpi_remove_item(iterm,proc_load(iproc))
556  call gaussfit_mpi_add_item(iterm,proc_load(jproc))
557  call gaussfit_mpi_remove_item(jterm,proc_load(jproc))
558 !Calculate final state
559  call gaussfit_mpi_calc_deviation(deviation2,nproc,proc_load)
560 !Swap them only if final state is better than the initial one
561  if(deviation2<deviation1) then
562    proc_dist(iterm)=jproc
563    proc_dist(jterm)=iproc
564  else
565 !  Return work load to initial state
566    call gaussfit_mpi_add_item(iterm,proc_load(iproc))
567    call gaussfit_mpi_remove_item(jterm,proc_load(iproc))
568    call gaussfit_mpi_add_item(jterm,proc_load(jproc))
569    call gaussfit_mpi_remove_item(iterm,proc_load(jproc))
570 !   write(*,*)'Back initial state'
571 !   write(*,*)'proc_load',proc_load(:)
572  end if
573 
574  end subroutine gaussfit_mpi_swap

m_paw_gaussfit/gaussfit_projector [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_projector

FUNCTION

 Fit tproj to Gaussians

INPUTS

  basis_size= size of the PAW basis
  orbitals= indicates the l quantum number for all orbitals.
  rpaw= PAW radius
  pawrad <type(pawrad_type)>=paw radial mesh and related data
  tproj= projectors
  maxterm= maximum number of terms used to fit the projectors.
  mparam= maximum number of parameters (Gaussian coefficients and factors) used.

OUTPUT

  nparam_array= number of parameters found.
  param = parameters found  (Gaussian coefficients and factors).

NOTES

 chisq=accuracy_p= sum_x abs(f(x)-y(x))/nx.
    nx is the number of points, f(x) and y(x) are the fitted and original functions.

SOURCE

2080 subroutine gaussfit_projector(basis_size,mparam,nparam_array,nterm_bounds,orbitals,param,pawrad,&
2081 & rpaw,tproj,comm_mpi)
2082 
2083 !Arguments ------------------------------------
2084  integer,intent(in) :: basis_size
2085  integer,intent(in) :: orbitals(basis_size)
2086  integer, optional,intent(in) :: comm_mpi
2087  real(dp),intent(in) :: rpaw
2088  type(pawrad_type),intent(in) :: pawrad
2089  real(dp),intent(in) :: tproj(:,:)
2090  integer,intent(in) :: mparam,nterm_bounds(2)
2091  integer,intent(out) :: nparam_array(basis_size)
2092  real(dp),intent(out) :: param(mparam,basis_size)
2093  type(pawrad_type)::mesh_tmp
2094 
2095 !Local variables ------------------------------
2096  integer :: ibasis,ierr,il,ir
2097  integer :: msz1,msz2,option
2098  real(dp) :: raux(1),rr(1)
2099  real(dp),allocatable :: d2(:),tproj_tmp1(:),tproj_tmp2(:)
2100  character(len=500) :: msg
2101  character(80) :: outfile
2102  !debug: uncomment
2103  !integer::i,nterm ,unitp
2104  !real(dp),allocatable::y(:)
2105  !end debug
2106 
2107 !************************************************************************
2108 
2109  if(size(tproj,2)<basis_size) then
2110    msg = 'wrong size for tproj in gaussfit_projector!'
2111    LIBPAW_BUG(msg)
2112  end if
2113 
2114  option=4  !see gaussfit_main
2115  nparam_array(:)=0
2116  msz1=min(pawrad_ifromr(pawrad,rpaw)+2,size(tproj,1))
2117  !msz1=pawrad%mesh_size
2118 
2119 !Augment the mesh size
2120 !this is to make the Gaussians go to zero after paw_radius
2121 !This is done by creating a new pawrad objet: mesh_tmp
2122 !Change to a linear grid:
2123  mesh_tmp%mesh_type=1 !linear grid
2124  mesh_tmp%rstep=0.0005 !very fine grid
2125  msz2=ceiling(pawrad%rmax*two/mesh_tmp%rstep)
2126  mesh_tmp%lstep=zero !only needed for log grids
2127  call pawrad_init(mesh_tmp,mesh_size=msz2,mesh_type=mesh_tmp%mesh_type,&
2128 & rstep=mesh_tmp%rstep,lstep=mesh_tmp%lstep)
2129 
2130  LIBPAW_ALLOCATE(tproj_tmp1,(msz1))
2131  LIBPAW_ALLOCATE(d2,(msz1))
2132  LIBPAW_ALLOCATE(tproj_tmp2,(msz2))
2133 
2134  do ibasis=1,basis_size
2135 
2136    write(msg,'(a," - Fitting wfn ",i4," to Gaussians")')ch10,ibasis
2137    call wrtout(std_out,  msg,'COLL')
2138 
2139    tproj_tmp1=zero; d2=zero; tproj_tmp2=zero
2140 
2141 !  take out r^il factor:
2142 !  il=psps%indlmn(1,ilmn,itypat)
2143    il=orbitals(ibasis)
2144 
2145    tproj_tmp1(2:msz1)=tproj(2:msz1,ibasis)/((pawrad%rad(2:msz1)+tol8)**(il))
2146 
2147 !  take out 1/r factor from eq.(3) of M. Torrent CMS 42, 337 (2008)
2148 !  since: <phi|proj>=1 from atompaw, and phi=phi*r, alors proj=proj/r
2149 
2150    tproj_tmp1(2:msz1)=tproj_tmp1(2:msz1)/(pawrad%rad(2:msz1))
2151    call pawrad_deducer0(tproj_tmp1(1:msz1),msz1,pawrad)
2152 
2153 !  splint to a different mesh:
2154 !  get second derivative of tproj and store it
2155    call paw_spline(pawrad%rad,tproj_tmp1(:),msz1,&
2156 &   zero,zero,d2)
2157 
2158    do ir=2,msz2
2159      rr=mesh_tmp%rad(ir)
2160      if( rr(1)-rpaw > tol8 ) then
2161        !after rpaw projectors are zero
2162        raux=zero
2163      else
2164        call paw_splint(msz1,pawrad%rad,&
2165 &       tproj_tmp1(:),d2(:),&
2166 &       1,rr,raux,ierr=ierr)
2167      end if
2168      tproj_tmp2(ir)=raux(1)
2169    end do
2170 
2171 !  Obtain the name for the output file
2172    if(ibasis<10) then
2173      write(outfile,'("wfn",i1,".fit")')ibasis
2174    elseif(ibasis<100) then
2175      write(outfile,'("wfn",i2,".fit")')ibasis
2176      write(msg,'(a,a,a,a)')ch10,&
2177 &     "ib (basis index) is too big!",ch10,&
2178 &     "Action: check your pseudopotentials"
2179      LIBPAW_BUG(msg)
2180    end if
2181 
2182    if(present(comm_mpi)) then
2183      call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,&
2184 &     param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2,comm_mpi)
2185    else
2186      call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,&
2187 &     param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2)
2188    end if
2189 
2190 !  check
2191 !  LIBPAW_ALLOCATE(y,(mesh_tmp%mesh_size))
2192 !  nterm=nparam_array(ibasis)/4
2193 !  call calcgaussc4(nparam_array(ibasis),nterm,mesh_tmp%mesh_size,1,param(:,ibasis),&
2194 !  &  mesh_tmp%rad,y)
2195 !  ! unitp=600+ibasis
2196 !  ! do ir=1,mesh_tmp%mesh_size
2197 !  !  write(unitp,'(2(f16.7,x,f16.7))')mesh_tmp%rad(ir),y(ir)
2198 !  ! end do
2199 !  LIBPAW_DEALLOCATE(y)
2200 
2201  end do
2202 
2203 !Deallocate
2204  call pawrad_free(mesh_tmp)
2205  LIBPAW_DEALLOCATE(tproj_tmp1)
2206  LIBPAW_DEALLOCATE(tproj_tmp2)
2207  LIBPAW_DEALLOCATE(d2)
2208 
2209 end subroutine gaussfit_projector

m_paw_gaussfit/gaussfit_rlsf [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_rlsf

FUNCTION

  Fits a given function to a sum of Gaussians.
  Uses the Levenberg-Marquardt algorithm.

COPYRIGHT

  Copyright (C) 2011-2024 ABINIT group (T. Rangel)
  This file is distributed under the terms of the
  GNU General Public License, see ~abinit/COPYING
  or http://www.gnu.org/copyleft/gpl.txt .
  The original Levemberg Marquardt routines were written by Armando Sole
  These were modified for the ARPUS spectra in the BigDFT code by A. Mirone.
  These were re-writen in Fortran and further modified in ABINIT for our particular needs.

INPUTS

  option=1 fit to a1 cos(a2 x^2)+ a3 sin( a4 x^2)
         2 fit to a1 exp(-a2 x^2)*(a3 cos (a4 x^2) + a5 sin (a6 x^2) )
         3 fit to a1 cos (k x^2) + a2 sin (k x^2)
         4 fit to a1 exp(-a2 x^2)* (a3 cos(k x^2)+ a4 sin (k x^2))
  if(option==1)mparam=nterm_bounds(2)*4
  if(option==2)mparam=nterm_bounds(2)*6
  if(option==3)mparam=nterm_bounds(2)*2
  if(option==4)mparam=nterm_bounds(2)*4

OUTPUT

SOURCE

1432 subroutine gaussfit_rlsf(&
1433 &chisq,constrains,limit,maxiter,&
1434 &nterm,nparam,nx,option,parameters,&
1435 &verbosity,weight,x,y)
1436 
1437 !Arguments -------------------------------
1438  real(dp),parameter::deltachi=tol10
1439  integer, intent(in) ::maxiter,nparam,nterm,nx
1440  integer, intent(in) ::option ,verbosity
1441  integer, intent(in) ::constrains(nparam)
1442  real(dp),intent(out)::chisq
1443  !arrays
1444  real(dp),intent(in)::limit(nparam),weight(nparam)
1445  real(dp),intent(inout)::parameters(nparam)
1446  real(dp),intent(in)::x(nx),y(nx)
1447 
1448 !Local variables-------------------------------
1449  integer::flag,ii,info,iter,jj,niter
1450  real(dp):: deltax
1451  real(dp)::chisq0,flambda,eta,lastdeltachi
1452  integer::ipvt(nparam)
1453  real(dp)::alpha(nparam,nparam)
1454  real(dp)::alpha0(nparam,nparam),beta(nparam)
1455  real(dp)::deltapar(nparam)
1456  real(dp)::tmp1(nparam,nparam)
1457  real(dp)::work(nparam)
1458  real(dp)::workpar(nparam)
1459  real(dp)::yfit(nx)
1460  character(len=500) :: msg
1461 
1462 ! *********************************************************************
1463 
1464 !
1465  !flambda=1e-6
1466  flambda=1e-7
1467 !iter=maxiter !later it is changed
1468  niter=0
1469  deltax=x(2)-x(1) !we assume this is a linear grid
1470 !
1471  iter_loop: do iter=1,maxiter
1472 !
1473    call gaussfit_chisq_alpha_beta(alpha0,beta,chisq0,&
1474 &   nparam,nterm,nx,option,parameters,x,y)
1475 !
1476    flag=0
1477    lastdeltachi=chisq0
1478 !
1479 !
1480    while_flag: do
1481      if(flag .ne. 0) exit while_flag
1482 !
1483      tmp1=0.d0
1484      do ii=1,nparam
1485        tmp1(ii,ii)=1.d0*flambda  !identity matrix * flambda
1486      end do
1487      alpha=alpha0+tmp1*alpha0
1488 !    Invert alpha matrix
1489      tmp1=alpha
1490      call dgetrf(nparam,nparam,tmp1,nparam,ipvt,info)
1491      if (.not.info==0) then
1492        if(verbosity>1) then
1493          write(msg,'(a)')'Matrix is singular'
1494          call wrtout(std_out,msg,'COLL')
1495        end if
1496        chisq=-1.d0
1497        exit iter_loop
1498      end if
1499      call dgetri(nparam,tmp1,nparam,ipvt,work,nparam,info)
1500      deltapar=0.d0
1501      if (.not.info==0) then
1502        if(verbosity>2) then
1503          write(msg,'(a)')'Matrix is singular'
1504          call wrtout(std_out,msg,'COLL')
1505        end if
1506        chisq=-1.d0
1507        exit iter_loop
1508      end if
1509 !
1510      if(tmp1(1,1) .ne. tmp1(1,1)) then !If is NaN
1511        chisq=-1.d0
1512        exit iter_loop
1513      end if
1514      if(abs(tmp1(1,1)) == tmp1(1,1)*tmp1(1,1)) then !If is infinity
1515        chisq=-1.d0
1516        exit iter_loop
1517      end if
1518 !
1519      do ii=1,nparam
1520        do jj=1,nparam
1521          deltapar(ii)=deltapar(ii)+beta(jj)*tmp1(jj,ii)
1522        end do
1523      end do
1524 !    apply constrains
1525      workpar(1:nparam)=parameters(1:nparam)+deltapar(1:nparam)*weight(1:nparam)
1526      call gaussfit_apply_constrains(constrains,limit,nparam,workpar)
1527 !
1528      if(option==1) then
1529        call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,workpar,x,yfit)
1530      elseif(option==2) then
1531        call gaussfit_calc_deriv_c(nparam,nterm,nx,1,workpar,x,yfit)
1532      elseif(option==3) then
1533        call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,workpar,x,yfit)
1534      elseif(option==4) then
1535        call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,workpar,x,yfit)
1536      end if
1537      chisq=0.d0
1538      do ii=1,nx
1539        chisq=chisq + ((y(ii)-yfit(ii)))**2
1540      end do
1541      chisq=chisq*deltax
1542 !
1543 !    write(*,'("chisq ",f12.5,"  chisq0 ",f12.5)')chisq,chisq0
1544 !
1545      if(chisq > chisq0) then
1546        flambda=flambda*2.0d0
1547        if( flambda > 1000.d0) then
1548          flag=1
1549 !        iter=0
1550          if(verbosity>2) then
1551            write(msg,'(a)')'flambda > 1000.d0'
1552            call wrtout(std_out,msg,'COLL')
1553          end if
1554          exit iter_loop
1555        end if
1556      else
1557        flag=1
1558        parameters=workpar
1559        eta=0.d0
1560        lastdeltachi=(chisq0-chisq) !/(chisq0+eta)
1561        if(lastdeltachi<deltachi) cycle
1562        chisq0=chisq
1563        flambda=flambda/2.d0
1564        if(verbosity>2) then
1565          write(msg,'("iter = ",i4," chisq = ",e15.6)')iter,chisq
1566          call wrtout(std_out,msg,'COLL')
1567        end if
1568      end if
1569    end do while_flag
1570  end do iter_loop
1571 
1572 end subroutine gaussfit_rlsf

m_paw_gaussfit/gaussfit_set_param1 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_set_param1

FUNCTION

  Sets parameters for LSF

INPUTS

OUTPUT

SOURCE

1674 subroutine gaussfit_set_param1(nterm,nparam,nx,param,sep,x,y)
1675 
1676 !Arguments -------------------------------
1677  integer,intent(in)::nterm,nparam,nx
1678  real(dp),intent(in)::sep
1679  real(dp),intent(in)::x(nx),y(nx)
1680  real(dp),intent(out)::param(nparam)
1681 
1682 !Local variables-------------------------------
1683  integer::ii,jj
1684  real(dp)::raux
1685 
1686 ! *********************************************************************
1687 
1688 !
1689  param(:)=1.0d0
1690 !exps=1.0/(x(nx)**2)
1691 !
1692 !alpha1
1693 
1694 !raux=maxval( y(:),nx )
1695 !maxval gives problems in some architectures:
1696  raux=-9999999
1697  do ii=1,nx
1698    if(raux<y(ii)) raux=y(ii)
1699  end do
1700  param(1:nterm)=raux
1701 !
1702 !alpha2
1703 !
1704 !y(r_c)=e^{-\alpha1 r_c}
1705  raux=-log(abs(y(nx))+tol10)/(x(nx)**2)
1706  param(nterm+1:nterm+2)=raux
1707 !
1708 !raux=0.5d0*pi/real(nterm,dp)
1709  do jj=1,nterm
1710    ii=jj+3*nterm
1711    param(ii)=sep**(jj)
1712 !  param(ii)=raux*real(jj,dp)
1713  end do
1714 !
1715  do jj=1,nterm
1716    ii=jj+5*nterm
1717    param(ii)=sep**(jj)
1718 !  param(ii)=raux*real(jj,dp)
1719  end do
1720 
1721 end subroutine gaussfit_set_param1

m_paw_gaussfit/gaussfit_set_param2 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_set_param2

FUNCTION

  Sets parameters for LSF

INPUTS

OUTPUT

SOURCE

1739 subroutine gaussfit_set_param2(nterm,nparam,nx,param,rpaw,x,y)
1740 
1741 !Arguments -------------------------------
1742  integer,intent(in)::nterm,nparam,nx
1743  real(dp),intent(in)::rpaw
1744  real(dp),intent(in)::x(nx),y(nx)
1745  real(dp),intent(out)::param(nparam)
1746 
1747 !Local variables-------------------------------
1748  integer::ii,jj
1749  real(dp)::exps,raux,sig,step
1750 
1751 ! *************************************************************************
1752 
1753  step=rpaw/real(nterm-1,dp)
1754 !exps=1.0/(rpaw/(real(nterm,dp)/2.d0))**2
1755  exps=1.0/(step*1.0d0)**2
1756 !alpha2 (width of gaussians)
1757 !Set to exps*real(nterm,dp)**2, for a good guess
1758  param(nterm+1:2*nterm)=exps !*real(nterm,dp)**2
1759 !
1760  do jj=1,nterm
1761 !  alpha3
1762 !  set to constant values of x
1763    ii=jj+2*nterm
1764    raux=step*real(jj-1,dp)
1765    param(ii)=raux
1766 !  alpha1
1767 !  set to the value of y at that point
1768    call gaussfit_param2_findsign()
1769    param(jj)=sig
1770  end do
1771 !
1772  contains

m_paw_gaussfit/gaussfit_set_param3 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_set_param3

FUNCTION

  Sets parameters for LSF

INPUTS

OUTPUT

SOURCE

1831 subroutine gaussfit_set_param3(nterm,nparam,param,sep)
1832 
1833 !Arguments -------------------------------
1834  integer,intent(in)::nterm,nparam
1835  real(dp),intent(in)::sep
1836  !real(dp),intent(in)::x(nx),y(nx)
1837  real(dp),intent(out)::param(nparam)
1838 
1839 !Local variables-------------------------------
1840  integer::ii,jj
1841 
1842 ! *********************************************************************
1843 
1844  param(:)=1.0d0
1845 !
1846  do jj=1,nterm
1847    ii=jj+nterm
1848    param(ii)=sep**(jj)
1849 !  param(ii)=raux*real(i,dp)
1850  end do
1851 !
1852  do jj=1,nterm
1853    ii=jj+3*nterm
1854    param(ii)=sep**(jj)
1855 !  param(ii)=raux*real(i,dp)
1856  end do
1857 
1858 end subroutine gaussfit_set_param3

m_paw_gaussfit/gaussfit_set_param4 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_set_param4

FUNCTION

  Sets parameters for LSF

INPUTS

OUTPUT

SOURCE

1876 subroutine gaussfit_set_param4(nparam,param)
1877 
1878 !Arguments -------------------------------
1879  integer,intent(in)::nparam
1880  real(dp),intent(out)::param(nparam)
1881 
1882 !Local variables-------------------------------
1883 
1884 ! *********************************************************************
1885 
1886  param(:)=1.0d0
1887 
1888 end subroutine gaussfit_set_param4

m_paw_gaussfit/gaussfit_set_param5 [ Functions ]

[ Top ] [ m_paw_gaussfit ] [ Functions ]

NAME

  gaussfit_set_param5

FUNCTION

  Sets parameters for LSF

INPUTS

OUTPUT

SOURCE

1906 subroutine gaussfit_set_param5(nterm,nparam,nx,param,rpaw,y)
1907 
1908 !Arguments -------------------------------
1909  integer,intent(in)::nterm,nparam,nx
1910  real(dp),intent(in)::rpaw
1911  real(dp),intent(in)::y(nx)
1912  real(dp),intent(out)::param(nparam)
1913 
1914 !Local variables-------------------------------
1915  integer::ix
1916  real(dp)::raux,a1,r_c,m
1917 
1918 ! *********************************************************************
1919 
1920  param(:)=1.0d0
1921 !
1922 !alpha1
1923 !a1=maxval( y(:),nx )
1924 !maxval gives problems in some architectures:
1925  a1=-9999999
1926  do ix=1,nx
1927    if(a1<y(ix)) a1=y(ix)
1928  end do
1929  param(1:nterm)=a1
1930 !
1931 !alpha2
1932 !
1933  r_c=rpaw+0.5d0 !paw sphere + a bit more.
1934 !this is not arbitrary since it is an initial guess
1935  m=0.01d0
1936  raux=log(a1/m)/r_c**2
1937  param(nterm+1:nterm*2)=raux
1938 
1939 end subroutine gaussfit_set_param5