TABLE OF CONTENTS
- ABINIT/m_paw_gaussfit
- gaussfit_set_param2/gaussfit_param2_findsign
- m_paw_gaussfit/gaussfit_apply_constrains
- m_paw_gaussfit/gaussfit_calc_deriv_c
- m_paw_gaussfit/gaussfit_calc_deriv_c2
- m_paw_gaussfit/gaussfit_calc_deriv_c3
- m_paw_gaussfit/gaussfit_calc_deriv_c4
- m_paw_gaussfit/gaussfit_calc_deriv_r
- m_paw_gaussfit/gaussfit_chisq_alpha_beta
- m_paw_gaussfit/gaussfit_constrains_init
- m_paw_gaussfit/gaussfit_fit
- m_paw_gaussfit/gaussfit_main
- m_paw_gaussfit/gaussfit_mpi_add_item
- m_paw_gaussfit/gaussfit_mpi_assign
- m_paw_gaussfit/gaussfit_mpi_calc_deviation
- m_paw_gaussfit/gaussfit_mpi_main
- m_paw_gaussfit/gaussfit_mpi_remove_item
- m_paw_gaussfit/gaussfit_mpi_set_weight
- m_paw_gaussfit/gaussfit_mpi_swap
- m_paw_gaussfit/gaussfit_projector
- m_paw_gaussfit/gaussfit_rlsf
- m_paw_gaussfit/gaussfit_set_param1
- m_paw_gaussfit/gaussfit_set_param2
- m_paw_gaussfit/gaussfit_set_param3
- m_paw_gaussfit/gaussfit_set_param4
- m_paw_gaussfit/gaussfit_set_param5
ABINIT/m_paw_gaussfit [ 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