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-2018 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
21 #include "libpaw.h" 22 23 module m_paw_gaussfit 24 25 USE_DEFS 26 USE_MSG_HANDLING 27 USE_MPI_WRAPPERS 28 USE_MEMORY_PROFILING 29 30 use m_paw_numeric, only : paw_splint, paw_spline 31 use m_pawrad, only : pawrad_type, pawrad_init, pawrad_deducer0, pawrad_free, pawrad_ifromr 32 33 implicit none 34 35 private 36 37 public:: gaussfit_projector !fit non-local projectors to a sum of gaussians 38 public:: gaussfit_main !main routine to fit 39 !Routines related to MPI: 40 private:: gaussfit_mpi_set_weight 41 private:: gaussfit_mpi_remove_item 42 private:: gaussfit_mpi_add_item 43 private:: gaussfit_mpi_assign 44 private:: gaussfit_mpi_main 45 private:: gaussfit_mpi_calc_deviation 46 private:: gaussfit_mpi_swap 47 !Routines related to fitting: 48 private:: gaussfit_fit !fit a function to gaussians 49 private:: gaussfit_calc_deriv_r !calc. derivatives for real gauss. 50 private:: gaussfit_calc_deriv_c !calc. derivatives for cplex. gauss. of form 1 51 private:: gaussfit_calc_deriv_c2 !calc. derivatives for cplex. gauss. of form 2 52 private:: gaussfit_calc_deriv_c3 !calc. derivatives for cplex. gauss. of form 3 53 private:: gaussfit_calc_deriv_c4 !calc. derivatives for cplex. gauss. of form 4 54 private:: gaussfit_rlsf !retreined least squares fit 55 private:: gaussfit_chisq_alpha_beta 56 !set parameters for LSF (for 5 different forms of gauss. sums): 57 private:: gaussfit_set_param1 58 private:: gaussfit_set_param2 59 private:: gaussfit_set_param3 60 private:: gaussfit_set_param4 61 private:: gaussfit_set_param5 62 private:: gaussfit_constrains_init !initialize constrains 63 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2084 subroutine gaussfit_param2_findsign() 2085 2086 !Arguments ------------------------------- 2087 !Local variables------------------------------- 2088 2089 !This section has been created automatically by the script Abilint (TD). 2090 !Do not modify the following lines by hand. 2091 #undef ABI_FUNC 2092 #define ABI_FUNC 'gaussfit_param2_findsign' 2093 !End of the abilint section 2094 2095 integer::ix,minx 2096 real(dp)::dist,mindist,xx,yy 2097 2098 ! ********************************************************************* 2099 2100 mindist=rpaw 2101 do ix=1,nx 2102 xx=x(ix) 2103 dist=abs(raux-xx) 2104 if(dist<mindist) then 2105 mindist=dist 2106 minx=ix 2107 end if 2108 end do 2109 yy=y(minx) 2110 sig=yy 2111 2112 end subroutine gaussfit_param2_findsign 2113 2114 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2399 subroutine gaussfit_apply_constrains(const,limit,nparam,ioparams) 2400 2401 2402 !This section has been created automatically by the script Abilint (TD). 2403 !Do not modify the following lines by hand. 2404 #undef ABI_FUNC 2405 #define ABI_FUNC 'gaussfit_apply_constrains' 2406 !End of the abilint section 2407 2408 implicit none 2409 2410 !Arguments ------------------------------- 2411 integer,intent(in):: nparam 2412 integer,intent(in):: const(nparam) 2413 real(dp),intent(in):: limit(nparam) 2414 real(dp),intent(inout):: ioparams(nparam) 2415 2416 !Local variables------------------------------- 2417 integer::ii 2418 2419 ! ********************************************************************* 2420 2421 do ii=1,nparam 2422 if(const(ii)==restricted .or. const(ii)==restricted_and_positive) then 2423 if(ioparams(ii)<limit(ii)) ioparams(ii)=limit(ii) 2424 end if 2425 if(const(ii)==positive .or. const(ii)==restricted_and_positive ) ioparams(ii)=abs(ioparams(ii)) 2426 2427 end do 2428 2429 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1349 subroutine gaussfit_calc_deriv_c(nparam,nterm,nx,opt,param,x,y_out,& 1350 & deriv) ! optional 1351 1352 1353 !This section has been created automatically by the script Abilint (TD). 1354 !Do not modify the following lines by hand. 1355 #undef ABI_FUNC 1356 #define ABI_FUNC 'gaussfit_calc_deriv_c' 1357 !End of the abilint section 1358 1359 implicit none 1360 1361 !Arguments ------------------------------- 1362 integer,intent(in)::nparam !number of parameters 1363 integer,intent(in)::nterm !number of gaussian expressions 1364 integer,intent(in)::nx !number of point in the x grid 1365 integer,intent(in)::opt !option: 1366 !1) calculate only f(x) 1367 !2) calculate f(x) and its derivatives 1368 real(dp),intent(in)::param(nparam) !parameters 1369 real(dp),intent(in)::x(nx) !xgrid 1370 real(dp),intent(out)::y_out(nx) !f(x) 1371 real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives 1372 1373 !Local variables------------------------------- 1374 integer::iexp,ii 1375 real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm) 1376 real(dp)::alpha4(nterm),alpha5(nterm),alpha6(nterm) 1377 real(dp)::aux1(nx),aux2(nx) 1378 real(dp)::cos1(nx,nterm),cos2(nx,nterm),sin1(nx,nterm),sin2(nx,nterm) 1379 real(dp)::term1(nx,nterm),term2(nx,nterm) 1380 1381 ! ********************************************************************* 1382 1383 ! 1384 !Initialize 1385 ! 1386 y_out(:)=0.d0 1387 ! 1388 !Get parameters from param array: 1389 ! 1390 alpha1(:)=param(1:nterm) 1391 alpha2(:)=param(nterm+1:2*nterm) 1392 alpha3(:)=param(2*nterm+1:3*nterm) 1393 alpha4(:)=param(3*nterm+1:4*nterm) 1394 alpha5(:)=param(4*nterm+1:5*nterm) 1395 alpha6(:)=param(5*nterm+1:6*nterm) 1396 ! 1397 !calculate useful quantities 1398 ! 1399 do iexp=1,nterm 1400 aux1(:)=-alpha2(iexp)*x(:)**2 1401 term1(:,iexp)=alpha1(iexp)*exp(aux1(:)) 1402 end do 1403 ! 1404 do iexp=1,nterm 1405 aux1(:)=alpha4(iexp)*x(:)**2 1406 sin1(:,iexp)=sin(aux1(:)) 1407 ! 1408 aux1(:)=alpha6(iexp)*x(:)**2 1409 sin2(:,iexp)=sin(aux1(:)) 1410 ! 1411 aux1(:)=alpha4(iexp)*x(:)**2 1412 cos1(:,iexp)=cos(aux1(:)) 1413 ! 1414 aux1(:)=alpha6(iexp)*x(:)**2 1415 cos2(:,iexp)=cos(aux1(:)) 1416 end do 1417 ! 1418 do iexp=1,nterm 1419 aux1(:)=alpha3(iexp)*sin1(:,iexp) 1420 aux2(:)=alpha5(iexp)*cos2(:,iexp) 1421 term2(:,iexp)=aux1(:)+aux2(:) 1422 y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp) 1423 end do 1424 ! 1425 !Calculate derivatives: 1426 ! 1427 if(opt==2) then 1428 ! 1429 ! alpha1 1430 ! 1431 do iexp=1,nterm 1432 aux1(:)=term1(:,iexp)/alpha1(iexp) 1433 aux2(:)=aux1(:)*term2(:,iexp) 1434 deriv(:,iexp)=aux2(:) 1435 end do 1436 ! 1437 ! alpha2 1438 ! 1439 do iexp=1,nterm 1440 ii=nterm+iexp 1441 aux1(:)=-term1(:,iexp)*term2(:,iexp) 1442 aux2(:)=aux1(:)*x(:)**2 1443 deriv(:,ii)=aux2(:) 1444 end do 1445 ! 1446 ! alpha3 1447 ! 1448 do iexp=1,nterm 1449 ii=2*nterm+iexp 1450 aux1(:)=term1(:,iexp)*sin1(:,iexp) 1451 deriv(:,ii)=aux1(:) 1452 end do 1453 ! 1454 ! alpha4 1455 ! 1456 do iexp=1,nterm 1457 ii=3*nterm+iexp 1458 aux1(:)=term1(:,iexp)*alpha3(iexp) 1459 aux2(:)=cos1(:,iexp)*x(:)**2 1460 deriv(:,ii)=aux2(:)*aux1(:) 1461 end do 1462 ! 1463 ! alpha5 1464 ! 1465 do iexp=1,nterm 1466 ii=4*nterm+iexp 1467 aux1(:)=term1(:,iexp)*cos2(:,iexp) 1468 deriv(:,ii)=aux1(:) 1469 end do 1470 ! 1471 ! alpha6 1472 ! 1473 do iexp=1,nterm 1474 ii=5*nterm+iexp 1475 aux1(:)=-term1(:,iexp)*alpha5(iexp) 1476 aux2(:)=sin2(:,iexp)*x(:)**2 1477 deriv(:,ii)=aux1(:)*aux2(:) 1478 end do 1479 end if 1480 1481 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1219 subroutine gaussfit_calc_deriv_c2(nparam,nterm,nx,opt,param,x,y_out,& 1220 & deriv) ! optional 1221 1222 1223 !This section has been created automatically by the script Abilint (TD). 1224 !Do not modify the following lines by hand. 1225 #undef ABI_FUNC 1226 #define ABI_FUNC 'gaussfit_calc_deriv_c2' 1227 !End of the abilint section 1228 1229 implicit none 1230 1231 !Arguments ------------------------------- 1232 integer,intent(in)::nparam !number of param 1233 integer,intent(in)::nterm !number of gaussian expressions 1234 integer,intent(in)::nx !number of point in the x grid 1235 integer,intent(in)::opt !option: 1236 !1) calculate only f(x) 1237 !2) calculate f(x) and its derivatives 1238 real(dp),intent(in)::param(nparam) !parameters 1239 real(dp),intent(in)::x(nx) !xgrid 1240 real(dp),intent(out)::y_out(nx) !f(x) 1241 real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives 1242 1243 !Local variables------------------------------- 1244 integer::iexp,ii 1245 real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm) 1246 real(dp)::alpha4(nterm) 1247 real(dp)::term1(nx,nterm) 1248 real(dp)::sin1(nx,nterm),sin2(nx,nterm),cos1(nx,nterm),cos2(nx,nterm) 1249 real(dp)::aux1(nx),aux2(nx) 1250 1251 ! ********************************************************************* 1252 1253 ! 1254 !Initialize 1255 ! 1256 y_out(:)=0.d0 1257 ! 1258 !Get param from param array: 1259 ! 1260 alpha1(:)=param(1:nterm) 1261 alpha2(:)=param(nterm+1:2*nterm) 1262 alpha3(:)=param(2*nterm+1:3*nterm) 1263 alpha4(:)=param(3*nterm+1:4*nterm) 1264 ! 1265 !calculate useful quantities 1266 ! 1267 ! 1268 do iexp=1,nterm 1269 aux1(:)=alpha2(iexp)*x(:)**2 1270 sin1(:,iexp)=sin(aux1(:)) 1271 ! 1272 aux1(:)=alpha4(iexp)*x(:)**2 1273 sin2(:,iexp)=sin(aux1(:)) 1274 ! 1275 aux1(:)=alpha2(iexp)*x(:)**2 1276 cos1(:,iexp)=cos(aux1(:)) 1277 ! 1278 aux1(:)=alpha4(iexp)*x(:)**2 1279 cos2(:,iexp)=cos(aux1(:)) 1280 end do 1281 ! 1282 do iexp=1,nterm 1283 aux1(:)=alpha1(iexp)*sin1(:,iexp) 1284 aux2(:)=alpha3(iexp)*cos2(:,iexp) 1285 term1(:,iexp)=aux1(:)+aux2(:) 1286 y_out(:)=y_out(:)+term1(:,iexp) 1287 end do 1288 ! 1289 !Calculate derivatives: 1290 ! 1291 if(opt==2) then 1292 ! 1293 ! alpha1 1294 ! 1295 do iexp=1,nterm 1296 deriv(:,iexp)=sin1(:,iexp) 1297 end do 1298 ! 1299 ! alpha2 1300 ! 1301 do iexp=1,nterm 1302 ii=nterm+iexp 1303 aux1(:)=alpha1(iexp)*cos1(:,iexp)*x(:)**2 1304 deriv(:,ii)=aux1(:) 1305 end do 1306 ! 1307 ! alpha3 1308 ! 1309 do iexp=1,nterm 1310 ii=2*nterm+iexp 1311 deriv(:,ii)=cos2(:,iexp) 1312 end do 1313 ! 1314 ! alpha4 1315 ! 1316 do iexp=1,nterm 1317 ii=3*nterm+iexp 1318 aux1(:)=-alpha3(iexp)*sin2(:,iexp)*x(:)**2 1319 deriv(:,ii)=aux1(:) 1320 end do 1321 end if 1322 1323 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1109 subroutine gaussfit_calc_deriv_c3(nparam,nterm,nx,opt,param,x,y_out,& 1110 & deriv) ! optional 1111 1112 1113 !This section has been created automatically by the script Abilint (TD). 1114 !Do not modify the following lines by hand. 1115 #undef ABI_FUNC 1116 #define ABI_FUNC 'gaussfit_calc_deriv_c3' 1117 !End of the abilint section 1118 1119 implicit none 1120 1121 !Arguments ------------------------------- 1122 integer,intent(in)::nparam !number of parameters 1123 integer,intent(in)::nterm !number of gaussian expressions 1124 integer,intent(in)::nx !number of point in the x grid 1125 integer,intent(in)::opt !option: 1126 !1) calculate only f(x) 1127 !2) calculate f(x) and its derivatives 1128 real(dp),intent(in)::param(nparam) !parameters 1129 real(dp),intent(in)::x(nx) !xgrid 1130 real(dp),intent(out)::y_out(nx) !f(x) 1131 real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives 1132 1133 !Local variables------------------------------- 1134 integer::iexp,ii 1135 real(dp)::sep 1136 real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm) 1137 real(dp)::term1(nx,nterm) 1138 real(dp)::sin1(nx,nterm),cos1(nx,nterm) 1139 real(dp)::aux1(nx),aux2(nx) 1140 1141 ! ********************************************************************* 1142 1143 ! 1144 !Initialize 1145 ! 1146 y_out(:)=0.d0 1147 ! 1148 sep=1.2d0 1149 ! 1150 !Get param from param array: 1151 ! 1152 alpha1(:)=param(1:nterm) 1153 alpha2(:)=param(nterm+1:2*nterm) 1154 ! 1155 do ii=1,nterm 1156 alpha3(ii)=sep**(ii) 1157 end do 1158 ! 1159 !calculate useful quantities 1160 ! 1161 do iexp=1,nterm 1162 aux1(:)=alpha3(iexp)*x(:)**2 1163 ! 1164 sin1(:,iexp)=sin(aux1(:)) 1165 cos1(:,iexp)=cos(aux1(:)) 1166 end do 1167 ! 1168 do iexp=1,nterm 1169 aux1(:)=alpha1(iexp)*sin1(:,iexp) 1170 aux2(:)=alpha2(iexp)*cos1(:,iexp) 1171 term1(:,iexp)=aux1(:)+aux2(:) 1172 y_out(:)=y_out(:)+term1(:,iexp) 1173 end do 1174 ! 1175 !Calculate derivatives: 1176 ! 1177 if(opt==2) then 1178 ! 1179 ! alpha1 1180 ! 1181 do iexp=1,nterm 1182 deriv(:,iexp)=sin1(:,iexp) 1183 end do 1184 ! 1185 ! alpha2 1186 ! 1187 do iexp=1,nterm 1188 ii=nterm+iexp 1189 deriv(:,ii)=cos1(:,iexp) 1190 end do 1191 end if 1192 1193 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1507 subroutine gaussfit_calc_deriv_c4(nparam,nterm,nx,opt,param,x,y_out,& 1508 & deriv) ! optional 1509 1510 1511 !This section has been created automatically by the script Abilint (TD). 1512 !Do not modify the following lines by hand. 1513 #undef ABI_FUNC 1514 #define ABI_FUNC 'gaussfit_calc_deriv_c4' 1515 !End of the abilint section 1516 1517 implicit none 1518 1519 !Arguments ------------------------------- 1520 integer,intent(in)::nparam !number of parameters 1521 integer,intent(in)::nterm !number of gaussian expressions 1522 integer,intent(in)::nx !number of point in the x grid 1523 integer,intent(in)::opt !option: 1524 !1) calculate only f(x) 1525 !2) calculate f(x) and its derivatives 1526 real(dp),intent(in)::param(nparam) !parameters 1527 real(dp),intent(in)::x(nx) !xgrid 1528 real(dp),intent(out)::y_out(nx) !f(x) 1529 real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives 1530 1531 !Local variables------------------------------- 1532 integer::iexp,ii 1533 real(dp)::raux,sep 1534 real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm) 1535 real(dp)::alpha4(nterm),alpha5(nterm) 1536 real(dp)::aux1(nx),aux2(nx) 1537 real(dp)::cos1(nx,nterm),sin1(nx,nterm) 1538 real(dp)::term1(nx,nterm),term2(nx,nterm) 1539 1540 ! ********************************************************************* 1541 1542 ! 1543 !Initialize 1544 ! 1545 sep=1.1d0 1546 y_out(:)=0.d0 1547 1548 !Get parameters from param array: 1549 ! 1550 alpha1(:)=param(1:nterm) 1551 alpha2(:)=param(nterm+1:2*nterm) 1552 alpha3(:)=param(2*nterm+1:3*nterm) 1553 alpha4(:)=param(3*nterm+1:4*nterm) 1554 ! 1555 ! 1556 raux=(2.d0*pi)/real(nterm,dp) 1557 do ii=1,nterm 1558 alpha5(ii)=sep**(ii) 1559 ! alpha5(ii)=raux*real(ii-1,dp) 1560 end do 1561 ! 1562 !calculate useful quantities 1563 ! 1564 do iexp=1,nterm 1565 aux1(:)=-alpha2(iexp)*x(:)**2 1566 term1(:,iexp)=alpha1(iexp)*exp(aux1(:)) 1567 end do 1568 ! 1569 do iexp=1,nterm 1570 aux1(:)=alpha5(iexp)*x(:)**2 1571 ! 1572 sin1(:,iexp)=sin(aux1(:)) 1573 cos1(:,iexp)=cos(aux1(:)) 1574 end do 1575 ! 1576 do iexp=1,nterm 1577 aux1(:)=alpha3(iexp)*sin1(:,iexp) 1578 aux2(:)=alpha4(iexp)*cos1(:,iexp) 1579 term2(:,iexp)=aux1(:)+aux2(:) 1580 y_out(:)=y_out(:)+term1(:,iexp)*term2(:,iexp) 1581 end do 1582 ! 1583 !Calculate derivatives: 1584 ! 1585 if(opt==2) then 1586 ! 1587 ! alpha1 1588 ! 1589 do iexp=1,nterm 1590 aux1(:)=term1(:,iexp)/alpha1(iexp) 1591 aux2(:)=aux1(:)*term2(:,iexp) 1592 deriv(:,iexp)=aux2(:) 1593 end do 1594 ! 1595 ! alpha2 1596 ! 1597 do iexp=1,nterm 1598 ii=nterm+iexp 1599 aux1(:)=-term1(:,iexp)*term2(:,iexp) 1600 aux2(:)=aux1(:)*x(:)**2 1601 deriv(:,ii)=aux2(:) 1602 end do 1603 ! 1604 ! alpha3 1605 ! 1606 do iexp=1,nterm 1607 ii=2*nterm+iexp 1608 aux1(:)=term1(:,iexp)*sin1(:,iexp) 1609 deriv(:,ii)=aux1(:) 1610 end do 1611 ! 1612 ! alpha4 1613 ! 1614 do iexp=1,nterm 1615 ii=3*nterm+iexp 1616 aux1(:)=term1(:,iexp)*cos1(:,iexp) 1617 deriv(:,ii)=aux1(:) 1618 end do 1619 end if 1620 1621 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
PARENTS
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
988 subroutine gaussfit_calc_deriv_r(nterm,nparam,nx,opt,param,x,y_out,& 989 & deriv) ! optional 990 991 992 !This section has been created automatically by the script Abilint (TD). 993 !Do not modify the following lines by hand. 994 #undef ABI_FUNC 995 #define ABI_FUNC 'gaussfit_calc_deriv_r' 996 !End of the abilint section 997 998 implicit none 999 1000 !Arguments ------------------------------- 1001 integer,intent(in)::nx !number of point in the x grid 1002 integer,intent(in)::nparam !number of parameters 1003 integer,intent(in)::nterm !number of gaussian expressions 1004 integer,intent(in)::opt !option: 1005 !1) calculate only f(x) 1006 !2) calculate f(x) and its derivatives 1007 real(dp),intent(in)::param(nparam) !parameters 1008 real(dp),intent(in)::x(nx) !xgrid 1009 real(dp),intent(out)::y_out(nx) !f(x) 1010 real(dp),optional,intent(out)::deriv(nx,nparam) !derivatives 1011 1012 !Local variables------------------------------- 1013 integer::iexp,ii 1014 real(dp)::alpha1(nterm),alpha2(nterm),alpha3(nterm) 1015 real(dp)::term1(nx,nterm) 1016 real(dp)::aux1(nx) 1017 !real(dp)::step 1018 1019 ! ********************************************************************* 1020 1021 ! 1022 !Initialize 1023 ! 1024 y_out(:)=0.d0 1025 ! 1026 !Get parameters from parameters array: 1027 ! 1028 alpha1(:)=param(1:nterm) 1029 alpha2(:)=param(nterm+1:2*nterm) 1030 alpha3(:)=param(2*nterm+1:3*nterm) 1031 ! 1032 !alpha3 1033 !set to constant values of x 1034 !step=rpaw/real(nterm,dp) 1035 !do ii=1,nterm 1036 !raux=step*real(ii,dp) 1037 !alpha3(ii)=raux 1038 !end do 1039 ! 1040 ! 1041 ! 1042 !calculate useful quantities 1043 ! 1044 do iexp=1,nterm 1045 aux1(:)=-alpha2(iexp)*(x(:)-alpha3(iexp))**2 1046 term1(:,iexp)=alpha1(iexp)*exp(aux1(:)) 1047 end do 1048 ! 1049 do iexp=1,nterm 1050 y_out(:)=y_out(:)+term1(:,iexp) 1051 end do 1052 ! 1053 !Calculate derivatives: 1054 ! 1055 if(opt==2) then 1056 ! 1057 ! alpha1 1058 ! 1059 do iexp=1,nterm 1060 aux1(:)=term1(:,iexp)/alpha1(iexp) 1061 deriv(:,iexp)=aux1(:) 1062 end do 1063 ! 1064 ! alpha2 1065 ! 1066 do iexp=1,nterm 1067 ii=nterm+iexp 1068 aux1(:)=-term1(:,iexp)*(x(:)-alpha3(iexp)) 1069 deriv(:,ii)=aux1(:) 1070 deriv(:,ii)=0.1d0 1071 end do 1072 ! 1073 ! alpha3 1074 ! 1075 do iexp=1,nterm 1076 ii=2*nterm+iexp 1077 aux1(:)=term1(:,iexp)*2.d0*alpha2(iexp) 1078 aux1(:)=aux1(:)*(x(:)-alpha3(iexp)) 1079 deriv(:,ii)=aux1(:) 1080 end do 1081 end if 1082 1083 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-2018 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1847 subroutine gaussfit_chisq_alpha_beta(alpha,beta,chisq,& 1848 & nparam,nterm,nx,option,parameters,x,y) 1849 1850 1851 !This section has been created automatically by the script Abilint (TD). 1852 !Do not modify the following lines by hand. 1853 #undef ABI_FUNC 1854 #define ABI_FUNC 'gaussfit_chisq_alpha_beta' 1855 !End of the abilint section 1856 1857 implicit none 1858 1859 !Arguments ------------------------------- 1860 integer,intent(in)::nparam,nterm,nx 1861 integer,intent(in)::option 1862 real(dp),intent(out)::chisq 1863 real(dp),intent(in)::parameters(nparam),x(nx),y(nx) 1864 real(dp),intent(out)::alpha(nparam,nparam),beta(nparam) 1865 1866 !Local variables------------------------------- 1867 integer::ii,jj,kk 1868 real(dp)::deltax,help1 1869 !arrays 1870 real(dp)::deltay(nx),deriv(nx,nparam),derivi(nx) 1871 real(dp)::yfit(nx) 1872 real(dp)::help0(nx),help2(nx),help3(nparam) 1873 1874 ! ********************************************************************* 1875 1876 deltax=x(2)-x(1) !we assume a linear grid 1877 ! 1878 if(option==1) then 1879 call gaussfit_calc_deriv_c2(nparam,nterm,nx,2,parameters,x,yfit,deriv) 1880 elseif(option==2) then 1881 call gaussfit_calc_deriv_c(nparam,nterm,nx,2,parameters,x,yfit,deriv) 1882 elseif(option==3) then 1883 call gaussfit_calc_deriv_c3(nparam,nterm,nx,2,parameters,x,yfit,deriv) 1884 elseif(option==4) then 1885 call gaussfit_calc_deriv_c4(nparam,nterm,nx,2,parameters,x,yfit,deriv) 1886 end if 1887 deltay=y-yfit 1888 help0=deltay 1889 ! 1890 do ii=1,nparam 1891 derivi(:)=deriv(:,ii) 1892 help1=0.d0 1893 do jj=1,nx 1894 help1=help1+help0(jj)*derivi(jj) 1895 end do 1896 beta(ii)=help1 1897 ! help1 = innerproduct(deriv,weight*derivi) 1898 ! below I use help3 instead for the array dimenstions 1899 help3=0.d0 1900 do kk=1,nparam 1901 do jj=1,nx 1902 help3(kk)=help3(kk)+deriv(jj,kk)*derivi(jj) 1903 end do 1904 end do 1905 ! ! 1906 alpha(:,ii)=help3(:) 1907 end do 1908 ! 1909 help2(:)=help0(:)*deltay(:) 1910 chisq=sum(help2) 1911 chisq=chisq*deltax 1912 1913 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2314 subroutine gaussfit_constrains_init(cons1,cons2,limit,nparam,nterm,nx,option,rpaw,y) 2315 2316 2317 !This section has been created automatically by the script Abilint (TD). 2318 !Do not modify the following lines by hand. 2319 #undef ABI_FUNC 2320 #define ABI_FUNC 'gaussfit_constrains_init' 2321 !End of the abilint section 2322 2323 implicit none 2324 2325 !Arguments ------------------------------- 2326 integer,intent(in)::nterm,option,nparam,nx 2327 integer,intent(out)::cons2(nparam) 2328 real(dp),intent(in)::rpaw 2329 real(dp),intent(in)::y(nx) 2330 real(dp),intent(out)::cons1(nparam),limit(nparam) 2331 2332 !Local variables------------------------------- 2333 integer :: ix 2334 real(dp)::rc,a1,mm,raux 2335 2336 ! ********************************************************************* 2337 2338 ! 2339 !DEFAULT: no weight 2340 cons1(:)=1.d0 2341 limit(:)=0.d0 2342 cons2=1 2343 ! 2344 if(option==4) then 2345 cons1(1:nterm)=0.2d0 2346 cons1(nterm+1:nterm*2)=0.3d0 2347 end if 2348 ! 2349 if(option==4) then 2350 ! parameters 2351 2352 ! a1=maxval( y(:),nx)/real(nterm,dp) 2353 ! maxval gives problems in some architectures: 2354 a1=-9999999 2355 do ix=1,nx 2356 if(a1<y(ix)) a1=y(ix) 2357 end do 2358 a1=a1/real(nterm,dp) 2359 2360 mm=0.01 2361 ! 2362 rc=1.7d0*rpaw 2363 raux=log(a1/mm)/rc**2 2364 ! 2365 ! Constraint exponential of gaussians to be positive (multiplied by -1), 2366 ! so that it decays to zero. 2367 ! Constraint as well its value, so that it does not get too big 2368 ! and it decays soon, 2369 ! This prevents that a gaussian grows at a very large x value. 2370 cons2(nterm+1:nterm*2)=restricted_and_positive 2371 limit(nterm+1:nterm*2)=raux 2372 end if 2373 2374 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)
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
896 subroutine gaussfit_fit(chisq,constrains,& 897 & limit,maxiter,nparam,nterm,nx,option,outfile,param,& 898 & verbosity,weight,x,y,y_out) 899 900 901 !This section has been created automatically by the script Abilint (TD). 902 !Do not modify the following lines by hand. 903 #undef ABI_FUNC 904 #define ABI_FUNC 'gaussfit_fit' 905 !End of the abilint section 906 907 implicit none 908 909 !Arguments ------------------------------------ 910 integer, intent(in) :: maxiter,nparam 911 integer,intent(in) :: nterm,nx,option,verbosity 912 !real(dp),optional,intent(in)::rpaw 913 !arrays 914 integer,intent(in)::constrains(nparam) 915 real(dp),intent(in)::limit(nparam),weight(nparam) 916 real(dp),intent(in)::x(nx),y(nx) 917 real(dp),intent(inout)::param(nparam) 918 real(dp),intent(out)::chisq,y_out(nx) 919 character(80),intent(in)::outfile 920 921 !Local variables ------------------------------ 922 integer, parameter :: wfn_unit=1007 923 integer::ix 924 real(dp)::rerror 925 real(dp),allocatable::sy(:) 926 927 ! ************************************************************************* 928 929 LIBPAW_ALLOCATE(sy,(nx)) 930 931 sy(:)=1.0d0 932 933 call gaussfit_rlsf(& 934 & chisq,constrains,limit,maxiter,& 935 & nterm,nparam,nx,option,param(1:nparam),& 936 & verbosity,weight,x,y) 937 ! 938 if(verbosity>=1) then 939 if(option==1) then 940 call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,param,x,y_out) 941 elseif(option==2) then 942 call gaussfit_calc_deriv_c(nparam,nterm,nx,1,param,x,y_out) 943 elseif(option==3) then 944 call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,param,x,y_out) 945 elseif(option==4) then 946 call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,param,x,y_out) 947 end if 948 ! 949 ! 950 open(wfn_unit,file=outfile,form='formatted',status='unknown') 951 ! per_error=0.d0 952 do ix=1, nx 953 rerror=abs(y(ix)-y_out(ix)) 954 write(wfn_unit,'(6(e20.12,1x))')x(ix),y(ix),y_out(ix),rerror 955 end do 956 close(wfn_unit) 957 958 end if 959 960 LIBPAW_DEALLOCATE(sy) 961 962 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).
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
115 subroutine gaussfit_main(mparam,nparam_out,nterm_bounds,nr,& 116 & param_out,pawrad,option,outfile,rpaw,y,comm_mpi) 117 118 119 !This section has been created automatically by the script Abilint (TD). 120 !Do not modify the following lines by hand. 121 #undef ABI_FUNC 122 #define ABI_FUNC 'gaussfit_main' 123 !End of the abilint section 124 125 implicit none 126 127 !Arguments ------------------------------------ 128 integer,intent(in)::mparam,nr,nterm_bounds(2) 129 integer,intent(in)::option 130 integer,intent(in),optional:: comm_mpi 131 real(dp),intent(in)::rpaw 132 real(dp),intent(inout)::y(nr) 133 character(80),intent(in)::outfile 134 type(pawrad_type),intent(in) :: pawrad 135 integer,intent(out)::nparam_out 136 real(dp),intent(out)::param_out(mparam) 137 138 !Local variables------------------------------- 139 !scalars 140 logical,parameter::modify_y=.false. !used only for plotting purposes 141 integer :: ichisq,ierr,ii,jj 142 integer :: master,me 143 integer :: maxiter,minterm,my_chisq_size,counts_all 144 integer :: ngauss,nparam,nproc,nterm 145 integer :: verbosity 146 real(dp) :: chisq,chisq_min 147 !real(dp) :: T1,T2 !uncomment for timming 148 !arrays 149 integer::constrains(mparam) 150 integer::proc_dist(nterm_bounds(1):nterm_bounds(2)) 151 integer,allocatable::counts(:),disp(:),map_nterm(:) 152 real(dp)::limit(mparam),param_tmp(mparam),weight(mparam) 153 real(dp),allocatable::chisq_array(:),recv_buf(:),send_buf(:) 154 real(dp),allocatable::y_out(:) 155 character(len=500) :: msg 156 157 ! ************************************************************************* 158 159 !initialize variables 160 maxiter=200 161 ! 162 !initialize mpi quantities: 163 master=0; me=0; nproc=1; proc_dist=1; 164 if(present(comm_mpi)) then 165 me=xmpi_comm_rank(comm_mpi) 166 nproc=xmpi_comm_size(comm_mpi) 167 end if 168 if(nproc>1) then 169 ! Find distribution (master) 170 if(me==master) then 171 call gaussfit_mpi_main(nproc,nterm_bounds,proc_dist) 172 end if 173 ! send distribution to all processors 174 call xmpi_bcast(proc_dist(nterm_bounds(1):nterm_bounds(2)),& 175 & master,comm_mpi,ierr) 176 end if 177 !Set size of chisq treated by each proc 178 my_chisq_size=nterm_bounds(2)-nterm_bounds(1)+1 !all terms 179 if(nproc>1 .and. .not. me==master) then 180 do nterm=nterm_bounds(1),nterm_bounds(2) 181 if(.not. proc_dist(nterm)==me+1) cycle 182 my_chisq_size=my_chisq_size+1 183 end do 184 end if 185 186 ! 187 !Allocate objects 188 ! 189 LIBPAW_ALLOCATE(y_out,(nr)) 190 LIBPAW_ALLOCATE(chisq_array,(my_chisq_size)) 191 if(master==me ) then 192 LIBPAW_BOUND1_ALLOCATE(map_nterm,BOUNDS(nterm_bounds(1),nterm_bounds(2))) 193 jj=1 194 do ii=1,nproc 195 do nterm=nterm_bounds(1),nterm_bounds(2) 196 if(proc_dist(nterm)==ii) then 197 map_nterm(nterm)=jj 198 jj=jj+1 199 end if 200 end do 201 end do 202 end if 203 ! 204 !fill with zeros 205 ! 206 chisq_array=zero 207 nparam_out=0 208 y_out=0.d0 209 verbosity=1 !print the minimum at the terminal 210 ! 211 ichisq=0 212 do nterm=nterm_bounds(1),nterm_bounds(2) 213 ! mpi distribution 214 if(.not. proc_dist(nterm)==me+1) cycle 215 ichisq=ichisq+1 216 217 ! call CPU_TIME(T1) 218 219 if(option==1) nparam=nterm*4 220 if(option==2) nparam=nterm*6 221 if(option==3) nparam=nterm*2 222 if(option==4) nparam=nterm*4 223 ! set initial guess 224 ! if(option==1) then 225 ! call gaussfit_set_param3(nterm,nparam,param_tmp(1:nparam),sep(ii)) 226 ! elseif(option==2) then 227 ! call gaussfit_set_param1(nterm,nparam,nr,& 228 !& param_tmp(1:nparam),sep(ii),pawrad%rad(1:nr),y) 229 if(option==3) then 230 call gaussfit_set_param4(nparam,param_tmp(1:nparam)) 231 elseif(option==4) then 232 call gaussfit_set_param5(nterm,nparam,nr,& 233 & param_tmp(1:nparam),rpaw,y) 234 end if 235 ! 236 ! 237 call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),& 238 & limit(1:nparam),nparam,nterm,nr,option,rpaw,y) 239 ! 240 call gaussfit_fit(chisq,constrains(1:nparam),& 241 & limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,& 242 & verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out) 243 244 ! if there was an error, set chisq to very high 245 ! if there is an NaN, for instance: 246 if(abs(chisq+1.d0)<tol8) chisq=99999 247 if(abs(chisq)==chisq*chisq) chisq=99999 248 if(chisq .ne. chisq) chisq=99999 249 ! 250 chisq_array(ichisq)=chisq 251 252 ! call CPU_TIME(T2) 253 ! print *, 'Time for fit ', T2-T1, 'seconds.' 254 255 end do !nterm 256 257 !mpicast results 258 !send distribution to master 259 if(nproc>1) then 260 ! Prepare communications: 261 LIBPAW_ALLOCATE(counts,(nproc)) 262 counts(:)=0 263 do nterm=nterm_bounds(1),nterm_bounds(2) 264 counts(proc_dist(nterm))=counts(proc_dist(nterm))+1 265 end do 266 counts_all=sum(counts) 267 LIBPAW_ALLOCATE(send_buf,(counts(me+1))) 268 send_buf(:)=chisq_array(1:counts(me+1)) 269 if(me==master) then 270 LIBPAW_ALLOCATE(recv_buf,(counts_all)) 271 else 272 LIBPAW_ALLOCATE(recv_buf,(1)) 273 end if 274 LIBPAW_ALLOCATE(disp,(nproc)) 275 disp(1)=0 276 do ii=2,nproc 277 disp(ii)=disp(ii-1)+counts(ii-1) 278 end do 279 ! communicate all info to master 280 call xmpi_gatherv(send_buf,counts(me+1),recv_buf,& 281 & counts,disp,master,comm_mpi,ierr) 282 ! fill in chisq_array with all received info: 283 if(master==me) then 284 do ii=1,counts_all 285 chisq_array(ii)=recv_buf(ii) 286 end do 287 end if 288 ! Deallocate MPI arrays: 289 LIBPAW_DEALLOCATE(recv_buf) 290 LIBPAW_DEALLOCATE(counts) 291 LIBPAW_DEALLOCATE(disp) 292 LIBPAW_DEALLOCATE(send_buf) 293 end if 294 295 !Print out info: 296 if(me==master) then 297 write(msg,'(3a)')'Preliminary results (with only 200 iter.):',ch10,' ngauss chisq' 298 call wrtout(std_out,msg,'COLL') 299 do nterm=nterm_bounds(1),nterm_bounds(2) 300 if(option==1) ngauss=nterm*4 301 if(option==2) ngauss=nterm*4 302 if(option==3) ngauss=nterm*2 303 if(option==4) ngauss=nterm*2 304 write(msg,'(i4,2x,e13.6,1x)')ngauss,& 305 & chisq_array(map_nterm(nterm)) 306 call wrtout(std_out,msg,'COLL') 307 end do 308 end if 309 310 !get minterm for best accuracy: 311 if(me==master) then 312 chisq_min=999999999 313 do nterm=nterm_bounds(1),nterm_bounds(2) 314 if(chisq_array(map_nterm(nterm))<chisq_min) then 315 chisq_min=chisq_array(map_nterm(nterm)) 316 minterm=nterm 317 end if 318 end do 319 320 !run again with the best parameters 321 nterm=minterm 322 if(option==1)nparam=4*nterm 323 if(option==2)nparam=6*nterm 324 if(option==3)nparam=2*nterm 325 if(option==4)nparam=4*nterm 326 327 ! set initial guess 328 if (option==3) then 329 call gaussfit_set_param4(nparam,param_tmp(1:nparam)) 330 elseif (option==4) then 331 call gaussfit_set_param5(nterm,nparam,nr,& 332 & param_tmp(1:nparam),rpaw,y) 333 end if 334 ! 335 call gaussfit_constrains_init(weight(1:nparam),constrains(1:nparam),& 336 & limit(1:nparam),nparam,nterm,nr,option,rpaw,y) 337 ! 338 verbosity=1; 339 maxiter=1000 !this time we do more iterations 340 call gaussfit_fit(chisq,constrains(1:nparam),& 341 & limit(1:nparam),maxiter,nparam,nterm,nr,option,outfile,param_tmp,& 342 & verbosity,weight(1:nparam),pawrad%rad(1:nr),y,y_out) 343 344 ! Write out best solution 345 write(msg,'(3a)')"Best solution (with more iterations):",ch10," ngauss chisq" 346 call wrtout(std_out,msg,'COLL') 347 if(option==1) ngauss=nterm*4 348 if(option==2) ngauss=nterm*4 349 if(option==3) ngauss=nterm*2 350 if(option==4) ngauss=nterm*2 351 write(msg,'(i4,2x,e13.6,1x)')ngauss,chisq 352 call wrtout(std_out,msg,'COLL') 353 end if 354 355 !Fill output variables 356 !and communicate results to all procs: 357 if(option==1) then 358 nparam_out=minterm*4 359 elseif (option==2) then 360 nparam_out=minterm*6 361 elseif (option==3) then 362 nparam_out=minterm*2 363 elseif (option==4) then 364 nparam_out=minterm*4 365 end if 366 367 !communicate 368 if(nproc>1) then 369 call xmpi_bcast(nparam_out,master,comm_mpi,ierr) 370 call xmpi_bcast(param_tmp(1:nparam_out),master,comm_mpi,ierr) 371 end if !nproc>1 372 373 param_out(:)=param_tmp 374 375 if(modify_y) then 376 ! call xmpi_scatterv(y_out,nr,mpi_displ,y_out,nr,0,comm_mpi,ierr) 377 y=y_out !at output modify y for the fitted y 378 end if 379 380 LIBPAW_DEALLOCATE(y_out) 381 LIBPAW_DEALLOCATE(chisq_array) 382 if(me==master) then 383 LIBPAW_DEALLOCATE(map_nterm) 384 end if 385 386 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
519 subroutine gaussfit_mpi_add_item(iterm,pload) 520 521 522 !This section has been created automatically by the script Abilint (TD). 523 !Do not modify the following lines by hand. 524 #undef ABI_FUNC 525 #define ABI_FUNC 'gaussfit_mpi_add_item' 526 !End of the abilint section 527 528 implicit none 529 530 !Arguments ------------------------------------ 531 integer,intent(in)::iterm 532 integer,intent(inout)::pload 533 534 !Local variables ------------------------------ 535 integer:: f_i 536 537 !************************************************************************ 538 539 call gaussfit_mpi_set_weight(f_i,iterm) 540 pload=pload+f_i 541 542 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
697 subroutine gaussfit_mpi_assign(iterm,nproc,nterm_bounds,& 698 & proc_dist,proc_load) 699 700 701 !This section has been created automatically by the script Abilint (TD). 702 !Do not modify the following lines by hand. 703 #undef ABI_FUNC 704 #define ABI_FUNC 'gaussfit_mpi_assign' 705 !End of the abilint section 706 707 implicit none 708 709 !Arguments ------------------------------------ 710 integer,intent(in)::iterm,nproc,nterm_bounds(2) 711 integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc) 712 713 !Local variables ------------------------------ 714 integer:: iproc,jproc,dev 715 integer:: deviation(nproc),mindev 716 character(len=100) :: msg 717 718 !************************************************************************ 719 720 do iproc=1,nproc 721 !add this term to iproc 722 call gaussfit_mpi_add_item(iterm,proc_load(iproc)) 723 !calculate the deviation for this configuration 724 call gaussfit_mpi_calc_deviation(dev,nproc,proc_load) 725 deviation(iproc)=dev 726 !remove this term from iproc 727 call gaussfit_mpi_remove_item(iterm,proc_load(iproc)) 728 end do 729 730 !assign to jproc, the proc with minimal deviation above: 731 jproc=-1; mindev=999999999 732 do iproc=1,nproc 733 if(deviation(iproc)<mindev) then 734 mindev=deviation(iproc) 735 jproc=iproc 736 end if 737 end do 738 if(jproc==-1) then 739 ! One should not get here! 740 msg = 'error in accomodate_mpi' 741 MSG_BUG(msg) 742 end if 743 744 !assign this term for jproc 745 proc_dist(iterm)=jproc 746 call gaussfit_mpi_add_item(iterm,proc_load(jproc)) 747 748 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
566 subroutine gaussfit_mpi_calc_deviation(deviation,nproc,proc_load) 567 568 569 !This section has been created automatically by the script Abilint (TD). 570 !Do not modify the following lines by hand. 571 #undef ABI_FUNC 572 #define ABI_FUNC 'gaussfit_mpi_calc_deviation' 573 !End of the abilint section 574 575 implicit none 576 577 !Arguments ------------------------------------ 578 integer,intent(in)::nproc 579 integer,intent(in)::proc_load(nproc) 580 integer,intent(out)::deviation 581 582 !Local variables ------------------------------ 583 integer:: jproc,kproc,kload,jload 584 585 !************************************************************************ 586 587 ! deviation=0 588 ! do jproc=1,nproc 589 ! deviation=deviation+abs(proc_load(jproc)-ideal) 590 ! end do 591 592 deviation=0 593 do jproc=1,nproc 594 jload=proc_load(jproc) 595 do kproc=1,jproc 596 kload=proc_load(kproc) 597 deviation=deviation+abs(kload-jload) 598 end do 599 end do 600 601 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
773 subroutine gaussfit_mpi_main(nproc,nterm_bounds,proc_dist) 774 775 776 !This section has been created automatically by the script Abilint (TD). 777 !Do not modify the following lines by hand. 778 #undef ABI_FUNC 779 #define ABI_FUNC 'gaussfit_mpi_main' 780 !End of the abilint section 781 782 implicit none 783 784 !Arguments ------------------------------------ 785 integer,intent(in)::nproc,nterm_bounds(2) 786 integer,intent(out)::proc_dist(nterm_bounds(1):nterm_bounds(2)) 787 788 !Local variables ------------------------------ 789 integer:: dev1,dev2,ii 790 integer:: iproc,iterm,jterm,ngauss,weight 791 integer:: proc_load(nproc) 792 character(len=500) :: msg 793 794 !************************************************************************ 795 796 proc_load=0; proc_dist=0 !initializations 797 798 !1) get a first-trial distribution: 799 do iterm=nterm_bounds(2),nterm_bounds(1),-1 800 call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load) 801 end do 802 803 !Do the following 20 times 804 do ii=1,20 805 ! Calculate initial state 806 call gaussfit_mpi_calc_deviation(dev1,nproc,proc_load) 807 ! Try to swap tasks between two processors: 808 do iterm=nterm_bounds(2),nterm_bounds(1),-1 809 do jterm=nterm_bounds(2),nterm_bounds(1),-1 810 call gaussfit_mpi_swap(iterm,jterm,& 811 & nproc,nterm_bounds,proc_dist,proc_load) 812 end do 813 end do 814 !! Try to reassign tasks to different processors 815 do iterm=nterm_bounds(2),nterm_bounds(1),-1 816 iproc=proc_dist(iterm) 817 ! Remove this job from this node 818 call gaussfit_mpi_remove_item(iterm,proc_load(iproc)) 819 ! Accomodate this again: 820 call gaussfit_mpi_assign(iterm,nproc,nterm_bounds,proc_dist,proc_load) 821 end do 822 ! Calculate final state 823 call gaussfit_mpi_calc_deviation(dev2,nproc,proc_load) 824 ! If final state equals initial state, exit: 825 if(dev2 == dev1) exit 826 ! write(*,'(a)')'Deviation: ',dev2 827 end do 828 829 !Write down distribution: 830 write(msg,'(3a)') 'MPI distribution',ch10,'N. gauss, iproc, weight ' 831 call wrtout(std_out,msg,'COLL') 832 do iterm=nterm_bounds(2),nterm_bounds(1),-1 833 ngauss=iterm*2 834 call gaussfit_mpi_set_weight(weight,iterm) 835 write(msg,'(3(i4,1x))') ngauss,proc_dist(iterm),weight 836 call wrtout(std_out,msg,'COLL') 837 end do 838 write(msg,'(a)') 'Load per processor: ' 839 call wrtout(std_out,msg,'COLL') 840 do iproc=1,nproc 841 write(msg,'(i5,1x,i10)') iproc,proc_load(iproc) 842 call wrtout(std_out,msg,'COLL') 843 end do 844 845 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
472 subroutine gaussfit_mpi_remove_item(iterm,pload) 473 474 475 !This section has been created automatically by the script Abilint (TD). 476 !Do not modify the following lines by hand. 477 #undef ABI_FUNC 478 #define ABI_FUNC 'gaussfit_mpi_remove_item' 479 !End of the abilint section 480 481 implicit none 482 483 !Arguments ------------------------------------ 484 integer,intent(in)::iterm 485 integer,intent(inout)::pload 486 487 !Local variables ------------------------------ 488 integer:: f_i 489 490 !*********************************************************************** 491 492 call gaussfit_mpi_set_weight(f_i,iterm) 493 pload=pload-f_i 494 495 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
415 subroutine gaussfit_mpi_set_weight(f,x) 416 417 418 !This section has been created automatically by the script Abilint (TD). 419 !Do not modify the following lines by hand. 420 #undef ABI_FUNC 421 #define ABI_FUNC 'gaussfit_mpi_set_weight' 422 !End of the abilint section 423 424 implicit none 425 426 !Arguments ------------------------------------ 427 integer,intent(in)::x 428 integer,intent(out)::f 429 430 !Local variables ------------------------------ 431 real(dp)::a,b,c,d,ff,xx 432 433 !************************************************************************ 434 435 !The following parameters were obtained 436 !from the time (in seconds) 437 !it takes to fit a given projector 438 !using from 1 to 100 gaussians. 439 a = -0.374137d0 ! +/- 2.013 (538.1%) 440 b = 0.207854d0 ! +/- 0.3385 (162.8%) 441 c = 0.0266371d0 ! +/- 0.01534 (57.59%) 442 d = 0.000152476d0 ! +/- 0.0001978 (129.7%) 443 444 xx=real(x,dp) 445 ff=a+b*xx+c*xx**2+d*xx**3 446 f=max(1,ceiling(ff)) 447 448 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
625 subroutine gaussfit_mpi_swap(iterm,jterm,& 626 & nproc,nterm_bounds,proc_dist,proc_load) 627 628 629 !This section has been created automatically by the script Abilint (TD). 630 !Do not modify the following lines by hand. 631 #undef ABI_FUNC 632 #define ABI_FUNC 'gaussfit_mpi_swap' 633 !End of the abilint section 634 635 implicit none 636 637 !Arguments ------------------------------------ 638 integer,intent(in)::iterm,jterm,nproc,nterm_bounds(2) 639 integer,intent(inout)::proc_dist(nterm_bounds(1):nterm_bounds(2)),proc_load(nproc) 640 641 !Local variables ------------------------------ 642 integer:: deviation1,deviation2 643 integer:: iproc,jproc 644 645 !************************************************************************ 646 647 !Calculate initial state 648 call gaussfit_mpi_calc_deviation(deviation1,nproc,proc_load) 649 iproc=proc_dist(iterm) 650 jproc=proc_dist(jterm) 651 !Swap terms: 652 call gaussfit_mpi_add_item(jterm,proc_load(iproc)) 653 call gaussfit_mpi_remove_item(iterm,proc_load(iproc)) 654 call gaussfit_mpi_add_item(iterm,proc_load(jproc)) 655 call gaussfit_mpi_remove_item(jterm,proc_load(jproc)) 656 !Calculate final state 657 call gaussfit_mpi_calc_deviation(deviation2,nproc,proc_load) 658 !Swap them only if final state is better than the initial one 659 if(deviation2<deviation1) then 660 proc_dist(iterm)=jproc 661 proc_dist(jterm)=iproc 662 else 663 ! Return work load to initial state 664 call gaussfit_mpi_add_item(iterm,proc_load(iproc)) 665 call gaussfit_mpi_remove_item(jterm,proc_load(iproc)) 666 call gaussfit_mpi_add_item(jterm,proc_load(jproc)) 667 call gaussfit_mpi_remove_item(iterm,proc_load(jproc)) 668 ! write(*,*)'Back initial state' 669 ! write(*,*)'proc_load',proc_load(:) 670 end if 671 672 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.
PARENTS
m_pawpsp
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2467 subroutine gaussfit_projector(basis_size,mparam,nparam_array,nterm_bounds,orbitals,param,pawrad,& 2468 & rpaw,tproj,comm_mpi) 2469 2470 2471 !This section has been created automatically by the script Abilint (TD). 2472 !Do not modify the following lines by hand. 2473 #undef ABI_FUNC 2474 #define ABI_FUNC 'gaussfit_projector' 2475 !End of the abilint section 2476 2477 implicit none 2478 2479 !Arguments ------------------------------------ 2480 integer,intent(in) :: basis_size 2481 integer,intent(in) :: orbitals(basis_size) 2482 integer, optional,intent(in) :: comm_mpi 2483 real(dp),intent(in) :: rpaw 2484 type(pawrad_type),intent(in) :: pawrad 2485 real(dp),intent(in) :: tproj(:,:) 2486 integer,intent(in) :: mparam,nterm_bounds(2) 2487 integer,intent(out) :: nparam_array(basis_size) 2488 real(dp),intent(out) :: param(mparam,basis_size) 2489 type(pawrad_type)::mesh_tmp 2490 2491 !Local variables ------------------------------ 2492 integer :: ibasis,ierr,il,ir 2493 integer :: msz1,msz2,option 2494 real(dp) :: raux(1),rr(1) 2495 real(dp),allocatable :: d2(:),tproj_tmp1(:),tproj_tmp2(:) 2496 character(len=500) :: msg 2497 character(80) :: outfile 2498 !debug: uncomment 2499 !integer::i,nterm ,unitp 2500 !real(dp),allocatable::y(:) 2501 !end debug 2502 2503 !************************************************************************ 2504 2505 if(size(tproj,2)<basis_size) then 2506 msg = 'wrong size for tproj in gaussfit_projector!' 2507 MSG_BUG(msg) 2508 end if 2509 2510 option=4 !see gaussfit_main 2511 nparam_array(:)=0 2512 msz1=min(pawrad_ifromr(pawrad,rpaw)+2,size(tproj,1)) 2513 !msz1=pawrad%mesh_size 2514 2515 !Augment the mesh size 2516 !this is to make the Gaussians go to zero after paw_radius 2517 !This is done by creating a new pawrad objet: mesh_tmp 2518 !Change to a linear grid: 2519 mesh_tmp%mesh_type=1 !linear grid 2520 mesh_tmp%rstep=0.0005 !very fine grid 2521 msz2=ceiling(pawrad%rmax*two/mesh_tmp%rstep) 2522 mesh_tmp%lstep=zero !only needed for log grids 2523 call pawrad_init(mesh_tmp,mesh_size=msz2,mesh_type=mesh_tmp%mesh_type,& 2524 & rstep=mesh_tmp%rstep,lstep=mesh_tmp%lstep) 2525 2526 LIBPAW_ALLOCATE(tproj_tmp1,(msz1)) 2527 LIBPAW_ALLOCATE(d2,(msz1)) 2528 LIBPAW_ALLOCATE(tproj_tmp2,(msz2)) 2529 2530 do ibasis=1,basis_size 2531 2532 write(msg,'(a," - Fitting wfn ",i4," to Gaussians")')ch10,ibasis 2533 call wrtout(std_out, msg,'COLL') 2534 2535 tproj_tmp1=zero; d2=zero; tproj_tmp2=zero 2536 2537 ! take out r^il factor: 2538 ! il=psps%indlmn(1,ilmn,itypat) 2539 il=orbitals(ibasis) 2540 2541 tproj_tmp1(2:msz1)=tproj(2:msz1,ibasis)/((pawrad%rad(2:msz1)+tol8)**(il)) 2542 2543 ! take out 1/r factor from eq.(3) of M. Torrent CMS 42, 337 (2008) 2544 ! since: <phi|proj>=1 from atompaw, and phi=phi*r, alors proj=proj/r 2545 2546 tproj_tmp1(2:msz1)=tproj_tmp1(2:msz1)/(pawrad%rad(2:msz1)) 2547 call pawrad_deducer0(tproj_tmp1(1:msz1),msz1,pawrad) 2548 2549 ! splint to a different mesh: 2550 ! get second derivative of tproj and store it 2551 call paw_spline(pawrad%rad,tproj_tmp1(:),msz1,& 2552 & zero,zero,d2) 2553 2554 do ir=2,msz2 2555 rr=mesh_tmp%rad(ir) 2556 if( rr(1)-rpaw > tol8 ) then 2557 !after rpaw projectors are zero 2558 raux=zero 2559 else 2560 call paw_splint(msz1,pawrad%rad,& 2561 & tproj_tmp1(:),d2(:),& 2562 & 1,rr,raux,ierr=ierr) 2563 end if 2564 tproj_tmp2(ir)=raux(1) 2565 end do 2566 2567 ! Obtain the name for the output file 2568 if(ibasis<10) then 2569 write(outfile,'("wfn",i1,".fit")')ibasis 2570 elseif(ibasis<100) then 2571 write(outfile,'("wfn",i2,".fit")')ibasis 2572 write(msg,'(a,a,a,a)')ch10,& 2573 & "ib (basis index) is too big!",ch10,& 2574 & "Action: check your pseudopotentials" 2575 MSG_BUG(msg) 2576 end if 2577 2578 if(present(comm_mpi)) then 2579 call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,& 2580 & param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2,comm_mpi) 2581 else 2582 call gaussfit_main(mparam,nparam_array(ibasis),nterm_bounds,msz2,& 2583 & param(:,ibasis),mesh_tmp,option,outfile,rpaw,tproj_tmp2) 2584 end if 2585 2586 ! check 2587 ! LIBPAW_ALLOCATE(y,(mesh_tmp%mesh_size)) 2588 ! nterm=nparam_array(ibasis)/4 2589 ! call calcgaussc4(nparam_array(ibasis),nterm,mesh_tmp%mesh_size,1,param(:,ibasis),& 2590 ! & mesh_tmp%rad,y) 2591 ! ! unitp=600+ibasis 2592 ! ! do ir=1,mesh_tmp%mesh_size 2593 ! ! write(unitp,'(2(f16.7,x,f16.7))')mesh_tmp%rad(ir),y(ir) 2594 ! ! end do 2595 ! LIBPAW_DEALLOCATE(y) 2596 2597 end do 2598 2599 !Deallocate 2600 call pawrad_free(mesh_tmp) 2601 LIBPAW_DEALLOCATE(tproj_tmp1) 2602 LIBPAW_DEALLOCATE(tproj_tmp2) 2603 LIBPAW_DEALLOCATE(d2) 2604 2605 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-2018 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1664 subroutine gaussfit_rlsf(& 1665 &chisq,constrains,limit,maxiter,& 1666 &nterm,nparam,nx,option,parameters,& 1667 &verbosity,weight,x,y) 1668 1669 1670 !This section has been created automatically by the script Abilint (TD). 1671 !Do not modify the following lines by hand. 1672 #undef ABI_FUNC 1673 #define ABI_FUNC 'gaussfit_rlsf' 1674 !End of the abilint section 1675 1676 implicit none 1677 1678 !Arguments ------------------------------- 1679 real(dp),parameter::deltachi=tol10 1680 integer, intent(in) ::maxiter,nparam,nterm,nx 1681 integer, intent(in) ::option ,verbosity 1682 integer, intent(in) ::constrains(nparam) 1683 real(dp),intent(out)::chisq 1684 !arrays 1685 real(dp),intent(in)::limit(nparam),weight(nparam) 1686 real(dp),intent(inout)::parameters(nparam) 1687 real(dp),intent(in)::x(nx),y(nx) 1688 1689 !Local variables------------------------------- 1690 integer::flag,ii,info,iter,jj,niter 1691 real(dp):: deltax 1692 real(dp)::chisq0,flambda,eta,lastdeltachi 1693 integer::ipvt(nparam) 1694 real(dp)::alpha(nparam,nparam) 1695 real(dp)::alpha0(nparam,nparam),beta(nparam) 1696 real(dp)::deltapar(nparam) 1697 real(dp)::tmp1(nparam,nparam) 1698 real(dp)::work(nparam) 1699 real(dp)::workpar(nparam) 1700 real(dp)::yfit(nx) 1701 character(len=500) :: msg 1702 1703 ! ********************************************************************* 1704 1705 ! 1706 !flambda=1e-6 1707 flambda=1e-7 1708 !iter=maxiter !later it is changed 1709 niter=0 1710 deltax=x(2)-x(1) !we assume this is a linear grid 1711 ! 1712 iter_loop: do iter=1,maxiter 1713 ! 1714 call gaussfit_chisq_alpha_beta(alpha0,beta,chisq0,& 1715 & nparam,nterm,nx,option,parameters,x,y) 1716 ! 1717 flag=0 1718 lastdeltachi=chisq0 1719 ! 1720 ! 1721 while_flag: do 1722 if(flag .ne. 0) exit while_flag 1723 ! 1724 tmp1=0.d0 1725 do ii=1,nparam 1726 tmp1(ii,ii)=1.d0*flambda !identity matrix * flambda 1727 end do 1728 alpha=alpha0+tmp1*alpha0 1729 ! Invert alpha matrix 1730 tmp1=alpha 1731 call dgetrf(nparam,nparam,tmp1,nparam,ipvt,info) 1732 if (.not.info==0) then 1733 if(verbosity>1) then 1734 write(msg,'(a)')'Matrix is singular' 1735 call wrtout(std_out,msg,'COLL') 1736 end if 1737 chisq=-1.d0 1738 exit iter_loop 1739 end if 1740 call dgetri(nparam,tmp1,nparam,ipvt,work,nparam,info) 1741 deltapar=0.d0 1742 if (.not.info==0) then 1743 if(verbosity>2) then 1744 write(msg,'(a)')'Matrix is singular' 1745 call wrtout(std_out,msg,'COLL') 1746 end if 1747 chisq=-1.d0 1748 exit iter_loop 1749 end if 1750 ! 1751 if(tmp1(1,1) .ne. tmp1(1,1)) then !If is NaN 1752 chisq=-1.d0 1753 exit iter_loop 1754 end if 1755 if(abs(tmp1(1,1)) == tmp1(1,1)*tmp1(1,1)) then !If is infinity 1756 chisq=-1.d0 1757 exit iter_loop 1758 end if 1759 ! 1760 do ii=1,nparam 1761 do jj=1,nparam 1762 deltapar(ii)=deltapar(ii)+beta(jj)*tmp1(jj,ii) 1763 end do 1764 end do 1765 ! apply constrains 1766 workpar(1:nparam)=parameters(1:nparam)+deltapar(1:nparam)*weight(1:nparam) 1767 call gaussfit_apply_constrains(constrains,limit,nparam,workpar) 1768 ! 1769 if(option==1) then 1770 call gaussfit_calc_deriv_c2(nparam,nterm,nx,1,workpar,x,yfit) 1771 elseif(option==2) then 1772 call gaussfit_calc_deriv_c(nparam,nterm,nx,1,workpar,x,yfit) 1773 elseif(option==3) then 1774 call gaussfit_calc_deriv_c3(nparam,nterm,nx,1,workpar,x,yfit) 1775 elseif(option==4) then 1776 call gaussfit_calc_deriv_c4(nparam,nterm,nx,1,workpar,x,yfit) 1777 end if 1778 chisq=0.d0 1779 do ii=1,nx 1780 chisq=chisq + ((y(ii)-yfit(ii)))**2 1781 end do 1782 chisq=chisq*deltax 1783 ! 1784 ! write(*,'("chisq ",f12.5," chisq0 ",f12.5)')chisq,chisq0 1785 ! 1786 if(chisq > chisq0) then 1787 flambda=flambda*2.0d0 1788 if( flambda > 1000.d0) then 1789 flag=1 1790 ! iter=0 1791 if(verbosity>2) then 1792 write(msg,'(a)')'flambda > 1000.d0' 1793 call wrtout(std_out,msg,'COLL') 1794 end if 1795 exit iter_loop 1796 end if 1797 else 1798 flag=1 1799 parameters=workpar 1800 eta=0.d0 1801 lastdeltachi=(chisq0-chisq) !/(chisq0+eta) 1802 if(lastdeltachi<deltachi) cycle 1803 chisq0=chisq 1804 flambda=flambda/2.d0 1805 if(verbosity>2) then 1806 write(msg,'("iter = ",i4," chisq = ",e15.6)')iter,chisq 1807 call wrtout(std_out,msg,'COLL') 1808 end if 1809 end if 1810 end do while_flag 1811 end do iter_loop 1812 1813 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
PARENTS
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
1937 subroutine gaussfit_set_param1(nterm,nparam,nx,param,sep,x,y) 1938 1939 1940 !This section has been created automatically by the script Abilint (TD). 1941 !Do not modify the following lines by hand. 1942 #undef ABI_FUNC 1943 #define ABI_FUNC 'gaussfit_set_param1' 1944 !End of the abilint section 1945 1946 implicit none 1947 1948 !Arguments ------------------------------- 1949 integer,intent(in)::nterm,nparam,nx 1950 real(dp),intent(in)::sep 1951 real(dp),intent(in)::x(nx),y(nx) 1952 real(dp),intent(out)::param(nparam) 1953 1954 !Local variables------------------------------- 1955 integer::ii,jj 1956 real(dp)::raux 1957 1958 ! ********************************************************************* 1959 1960 ! 1961 param(:)=1.0d0 1962 !exps=1.0/(x(nx)**2) 1963 ! 1964 !alpha1 1965 1966 !raux=maxval( y(:),nx ) 1967 !maxval gives problems in some architectures: 1968 raux=-9999999 1969 do ii=1,nx 1970 if(raux<y(ii)) raux=y(ii) 1971 end do 1972 param(1:nterm)=raux 1973 ! 1974 !alpha2 1975 ! 1976 !y(r_c)=e^{-\alpha1 r_c} 1977 raux=-log(abs(y(nx))+tol10)/(x(nx)**2) 1978 param(nterm+1:nterm+2)=raux 1979 ! 1980 !raux=0.5d0*pi/real(nterm,dp) 1981 do jj=1,nterm 1982 ii=jj+3*nterm 1983 param(ii)=sep**(jj) 1984 ! param(ii)=raux*real(jj,dp) 1985 end do 1986 ! 1987 do jj=1,nterm 1988 ii=jj+5*nterm 1989 param(ii)=sep**(jj) 1990 ! param(ii)=raux*real(jj,dp) 1991 end do 1992 1993 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
PARENTS
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2017 subroutine gaussfit_set_param2(nterm,nparam,nx,param,rpaw,x,y) 2018 2019 2020 !This section has been created automatically by the script Abilint (TD). 2021 !Do not modify the following lines by hand. 2022 #undef ABI_FUNC 2023 #define ABI_FUNC 'gaussfit_set_param2' 2024 !End of the abilint section 2025 2026 implicit none 2027 2028 !Arguments ------------------------------- 2029 integer,intent(in)::nterm,nparam,nx 2030 real(dp),intent(in)::rpaw 2031 real(dp),intent(in)::x(nx),y(nx) 2032 real(dp),intent(out)::param(nparam) 2033 2034 !Local variables------------------------------- 2035 integer::ii,jj 2036 real(dp)::exps,raux,sig,step 2037 2038 ! ************************************************************************* 2039 2040 step=rpaw/real(nterm-1,dp) 2041 !exps=1.0/(rpaw/(real(nterm,dp)/2.d0))**2 2042 exps=1.0/(step*1.0d0)**2 2043 !alpha2 (width of gaussians) 2044 !Set to exps*real(nterm,dp)**2, for a good guess 2045 param(nterm+1:2*nterm)=exps !*real(nterm,dp)**2 2046 ! 2047 do jj=1,nterm 2048 ! alpha3 2049 ! set to constant values of x 2050 ii=jj+2*nterm 2051 raux=step*real(jj-1,dp) 2052 param(ii)=raux 2053 ! alpha1 2054 ! set to the value of y at that point 2055 call gaussfit_param2_findsign() 2056 param(jj)=sig 2057 end do 2058 ! 2059 contains
m_paw_gaussfit/gaussfit_set_param3 [ Functions ]
[ Top ] [ m_paw_gaussfit ] [ Functions ]
NAME
gaussfit_set_param3
FUNCTION
Sets parameters for LSF
INPUTS
OUTPUT
PARENTS
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2138 subroutine gaussfit_set_param3(nterm,nparam,param,sep) 2139 2140 2141 !This section has been created automatically by the script Abilint (TD). 2142 !Do not modify the following lines by hand. 2143 #undef ABI_FUNC 2144 #define ABI_FUNC 'gaussfit_set_param3' 2145 !End of the abilint section 2146 2147 implicit none 2148 2149 !Arguments ------------------------------- 2150 integer,intent(in)::nterm,nparam 2151 real(dp),intent(in)::sep 2152 !real(dp),intent(in)::x(nx),y(nx) 2153 real(dp),intent(out)::param(nparam) 2154 2155 !Local variables------------------------------- 2156 integer::ii,jj 2157 2158 ! ********************************************************************* 2159 2160 param(:)=1.0d0 2161 ! 2162 do jj=1,nterm 2163 ii=jj+nterm 2164 param(ii)=sep**(jj) 2165 ! param(ii)=raux*real(i,dp) 2166 end do 2167 ! 2168 do jj=1,nterm 2169 ii=jj+3*nterm 2170 param(ii)=sep**(jj) 2171 ! param(ii)=raux*real(i,dp) 2172 end do 2173 2174 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2199 subroutine gaussfit_set_param4(nparam,param) 2200 2201 2202 !This section has been created automatically by the script Abilint (TD). 2203 !Do not modify the following lines by hand. 2204 #undef ABI_FUNC 2205 #define ABI_FUNC 'gaussfit_set_param4' 2206 !End of the abilint section 2207 2208 implicit none 2209 2210 !Arguments ------------------------------- 2211 integer,intent(in)::nparam 2212 real(dp),intent(out)::param(nparam) 2213 2214 !Local variables------------------------------- 2215 2216 ! ********************************************************************* 2217 2218 param(:)=1.0d0 2219 2220 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
PARENTS
m_paw_gaussfit
CHILDREN
gaussfit_main,paw_spline,paw_splint,pawrad_deducer0,pawrad_free pawrad_init,wrtout
SOURCE
2245 subroutine gaussfit_set_param5(nterm,nparam,nx,param,rpaw,y) 2246 2247 2248 !This section has been created automatically by the script Abilint (TD). 2249 !Do not modify the following lines by hand. 2250 #undef ABI_FUNC 2251 #define ABI_FUNC 'gaussfit_set_param5' 2252 !End of the abilint section 2253 2254 implicit none 2255 2256 !Arguments ------------------------------- 2257 integer,intent(in)::nterm,nparam,nx 2258 real(dp),intent(in)::rpaw 2259 real(dp),intent(in)::y(nx) 2260 real(dp),intent(out)::param(nparam) 2261 2262 !Local variables------------------------------- 2263 integer::ix 2264 real(dp)::raux,a1,r_c,m 2265 2266 ! ********************************************************************* 2267 2268 param(:)=1.0d0 2269 ! 2270 !alpha1 2271 !a1=maxval( y(:),nx ) 2272 !maxval gives problems in some architectures: 2273 a1=-9999999 2274 do ix=1,nx 2275 if(a1<y(ix)) a1=y(ix) 2276 end do 2277 param(1:nterm)=a1 2278 ! 2279 !alpha2 2280 ! 2281 r_c=rpaw+0.5d0 !paw sphere + a bit more. 2282 !this is not arbitrary since it is an initial guess 2283 m=0.01d0 2284 raux=log(a1/m)/r_c**2 2285 param(nterm+1:nterm*2)=raux 2286 2287 end subroutine gaussfit_set_param5