QMCPACK
Public Types | Public Member Functions | Public Attributes | List of all members
LRBreakup< BreakupBasis > Struct Template Reference
+ Collaboration diagram for LRBreakup< BreakupBasis >:

Public Types

using ParticleLayout = ParticleSet::ParticleLayout
 

Public Member Functions

void AddKToList (mRealType k, mRealType degeneracy=1.0)
 
int SetupKVecs (mRealType kc, mRealType kcont, mRealType kmax)
 setup KList More...
 
mRealType DoBreakup (mRealType *Fk, mRealType *t, mRealType *adjust)
 
mRealType DoGradBreakup (mRealType *Fk, mRealType *t, mRealType *adjust)
 
mRealType DoStrainBreakup (mRealType *Fk, mRealType *dFk, mRealType *t, mRealType *adjust)
 
void DoAllBreakup (mRealType *chisqr, mRealType *Fk, mRealType *dFk, mRealType *t, mRealType *gt, mRealType *dt, mRealType *adjust)
 
mRealType DoBreakup (mRealType *Fk, mRealType *t)
 
 LRBreakup (BreakupBasis &bref)
 
mRealType DoGradBreakup (mRealType *Fk, mRealType *t)
 

Public Attributes

BreakupBasis & Basis
 The basis to be used for breakup. More...
 
std::vector< TinyVector< mRealType, 2 > > KList
 For each k, KList[k][0] = |k| and KList[k][1] = degeneracy. More...
 

Detailed Description

template<class BreakupBasis>
struct qmcplusplus::LRBreakup< BreakupBasis >

Definition at line 29 of file LRBreakup.h.

Member Typedef Documentation

◆ ParticleLayout

using ParticleLayout = ParticleSet::ParticleLayout

Definition at line 34 of file LRBreakup.h.

Constructor & Destructor Documentation

◆ LRBreakup()

LRBreakup ( BreakupBasis &  bref)
inline

Definition at line 191 of file LRBreakup.h.

191  : Basis(bref)
192  { /*Do Nothing*/
193  }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42

Member Function Documentation

◆ AddKToList()

void AddKToList ( mRealType  k,
mRealType  degeneracy = 1.0 
)

Definition at line 316 of file LRBreakup.h.

References qmcplusplus::abs().

317 {
318  //Search for this k already in list
319  int ki = 0;
320  while ((ki < KList.size()) && (std::abs(k - KList[ki][0]) > 1.0e-12))
321  ki++;
322  if (ki == KList.size())
323  {
324  TinyVector<mRealType, 2> temp(k, degeneracy);
325  KList.push_back(temp);
326  }
327  else
328  KList[ki][1] += degeneracy;
329 }
qmcplusplus::abs
MakeReturn< UnaryNode< FnFabs, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t abs(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:88
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoAllBreakup()

void DoAllBreakup ( mRealType *  chisqr,
mRealType *  Fk,
mRealType *  dFk,
mRealType *  t,
mRealType *  gt,
mRealType *  dt,
mRealType *  adjust 
)

Definition at line 894 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), qmcplusplus::Units::distance::m, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

901 {
902  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
903  //t and adjust must be allocated up to Basis.NumBasisElem();
904  //Fk must be allocated and filled up to KList.size();
905  // assert(t.size()==adjust.size());
906  // assert(t.size()==Basis.NumBasisElem());
907  Matrix<mRealType> A;
908  Matrix<mRealType> Af;
909  Matrix<mRealType> As;
910 
911  std::vector<mRealType> b;
912  std::vector<mRealType> bf;
913  std::vector<mRealType> bs;
914 
915  Matrix<mRealType> cnk;
916  Matrix<mRealType> dcnk;
917 
918  int N = Basis.NumBasisElem(); //t.size();
919 
920  A.resize(N, N);
921  Af.resize(N, N);
922  As.resize(N, N);
923 
924  b.resize(N, 0.0);
925  bf.resize(N, 0.0);
926  bs.resize(N, 0.0);
927 
928  cnk.resize(N, KList.size());
929  dcnk.resize(N, KList.size());
930  //Fill in cnk.
931  for (int n = 0; n < N; n++)
932  {
933  for (int ki = 0; ki < KList.size(); ki++)
934  {
935  mRealType k = KList[ki][0];
936  cnk(n, ki) = Basis.c(n, k);
937  dcnk(n, ki) = Basis.dc_dk(n, k); //-Basis.c(n,k);
938  }
939  }
940  //Fill in A and b
941  A = 0.0;
942  Af = 0.0;
943  As = 0.0;
944 
945  for (int l = 0; l < N; l++)
946  {
947  for (int ki = 0; ki < KList.size(); ki++)
948  {
949  mRealType k2 = KList[ki][0] * KList[ki][0];
950 
951  mRealType temp = KList[ki][1] * Fk[ki] * cnk(l, ki);
952  // b[l] += k2*KList[ki][1]*(dFk[ki]-Fk[ki]) * dcnk(l, ki);
953  b[l] += temp;
954  bf[l] += k2 * temp;
955  bs[l] += k2 * KList[ki][1] * dFk[ki] * dcnk(l, ki);
956 
957  for (int n = 0; n < N; n++)
958  {
959  temp = KList[ki][1] * cnk(l, ki) * cnk(n, ki);
960  A(l, n) += temp;
961  Af(l, n) += k2 * temp;
962  As(l, n) += k2 * KList[ki][1] * dcnk(l, ki) * dcnk(n, ki);
963  }
964  }
965  }
966  //************************************
967  //FOR POTENTIAL AND FORCE
968  //************************************
969 
970  //Reduce for constraints
971  int M = N;
972  for (int i = 0; i < N; i++)
973  if (!adjust[i])
974  M--;
975  //The c is for "constrained"
976  Matrix<mRealType> Ac;
977  Matrix<mRealType> Afc;
978  Matrix<mRealType> Asc;
979  Ac.resize(M, M);
980  Afc.resize(M, M);
981  Asc.resize(M, M);
982  std::vector<mRealType> bc(M, 0.0), bfc(M, 0.0), bsc(M, 0.0), tc(M, 0.0), tfc(M, 0.0), tsc(M, 0.0);
983  //Build constrained Ac and bc
984  int j = 0;
985  for (int col = 0; col < N; col++)
986  {
987  if (adjust[col])
988  {
989  // Copy column a A to Ac
990  int i = 0;
991  for (int row = 0; row < N; row++)
992  if (adjust[row])
993  {
994  Ac(i, j) = A(row, col);
995  Afc(i, j) = Af(row, col);
996  Asc(i, j) = As(row, col);
997  i++;
998  }
999  j++;
1000  }
1001  else
1002  {
1003  // Otherwise, subtract t(col)*A(:,col) from bc
1004  for (int row = 0; row < N; row++)
1005  {
1006  b[row] -= A(row, col) * t[col];
1007  bf[row] -= Af(row, col) * gt[col];
1008  bs[row] -= As(row, col) * dt[col];
1009  }
1010  }
1011  }
1012  j = 0;
1013  for (int row = 0; row < N; row++)
1014  if (adjust[row])
1015  {
1016  bc[j] = b[row];
1017  bfc[j] = bf[row];
1018  bsc[j] = bs[row];
1019  j++;
1020  }
1021  // Do SVD:
1022  // -------
1023  // Matrix<mRealType> U(M, M), V(M, M);
1024  // std::vector<mRealType> S(M), Sinv(M);
1025  // SVdecomp(Ac, U, S, V);
1026  ////////////////////////////////
1027  int m = Ac.rows();
1028  int n = Ac.cols();
1029  Matrix<mRealType> A_trans(n, m);
1030  Matrix<mRealType> Af_trans(n, m);
1031  Matrix<mRealType> As_trans(n, m);
1032  Matrix<mRealType> U, V;
1033  Matrix<mRealType> Uf, Vf;
1034  Matrix<mRealType> Us, Vs;
1035 
1036  U.resize(std::min(m, n), m);
1037  V.resize(n, std::min(m, n));
1038 
1039  Uf.resize(std::min(m, n), m);
1040  Vf.resize(n, std::min(m, n));
1041 
1042  Us.resize(std::min(m, n), m);
1043  Vs.resize(n, std::min(m, n));
1044 
1045  std::vector<mRealType> S, Sinv;
1046  S.resize(std::min(n, m));
1047 
1048  std::vector<mRealType> Sf, Sfinv;
1049  Sf.resize(std::min(n, m));
1050 
1051  std::vector<mRealType> Ss, Ssinv;
1052  Ss.resize(std::min(n, m));
1053  //do the transpose
1054  for (int i = 0; i < m; i++)
1055  {
1056  for (int j = 0; j < n; j++)
1057  {
1058  A_trans(j, i) = Ac(i, j);
1059  Af_trans(j, i) = Afc(i, j);
1060  As_trans(j, i) = Asc(i, j);
1061  }
1062  }
1063  char JOBU = 'S';
1064  char JOBVT = 'S';
1065  int LDA = m;
1066  int LDU = m;
1067  int LDVT = std::min(m, n);
1068  int LWORK = 10 * std::max(3 * std::min(n, m) + std::max(m, n), 5 * std::min(m, n));
1069  std::vector<mRealType> WORK(LWORK);
1070  int INFO;
1071  LAPACK::gesvd(JOBU, JOBVT, m, n, A_trans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
1072  assert(INFO == 0);
1073 
1074  LAPACK::gesvd(JOBU, JOBVT, m, n, Af_trans.data(), LDA, &Sf[0], Uf.data(), LDU, Vf.data(), LDVT, &WORK[0], LWORK,
1075  INFO);
1076  assert(INFO == 0);
1077 
1078  LAPACK::gesvd(JOBU, JOBVT, m, n, As_trans.data(), LDA, &Ss[0], Us.data(), LDU, Vs.data(), LDVT, &WORK[0], LWORK,
1079  INFO);
1080  assert(INFO == 0);
1081 
1082  int ur = U.rows();
1083  int uc = U.cols();
1084  Matrix<mRealType> U_trans(uc, ur);
1085  Matrix<mRealType> Uf_trans(uc, ur);
1086  Matrix<mRealType> Us_trans(uc, ur);
1087  for (int i = 0; i < ur; i++)
1088  {
1089  for (int j = 0; j < uc; j++)
1090  {
1091  U_trans(j, i) = U(i, j);
1092  Uf_trans(j, i) = Uf(i, j);
1093  Us_trans(j, i) = Us(i, j);
1094  }
1095  }
1096  U.resize(uc, ur);
1097  U = U_trans;
1098 
1099  Uf.resize(uc, ur);
1100  Uf = Uf_trans;
1101 
1102  Us.resize(uc, ur);
1103  Us = Us_trans;
1104  //////////////////////////////////
1105  // Zero out near-singular values
1106 
1107  //First, do normal breakup.
1108  mRealType Smax = S[0];
1109  for (int i = 1; i < M; i++)
1110  Smax = std::max(S[i], Smax);
1111  for (int i = 0; i < M; i++)
1112  if (S[i] < 0.0)
1113  std::cout << "negative singlar value.\n";
1114  // perr << "Smax = " << Smax << std::endl;
1115  Sinv.resize(S.size());
1116  for (int i = 0; i < M; i++)
1117  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
1118  int numSingular = 0;
1119  for (int i = 0; i < Sinv.size(); i++)
1120  if (Sinv[i] == 0.0)
1121  numSingular++;
1122  if (numSingular > 0)
1123  std::cout << "There were " << numSingular << " singular values in energy breakup.\n";
1124  // Compute t_n, removing singular values
1125 
1126  //Second, do force.
1127  Smax = Sf[0];
1128  for (int i = 1; i < M; i++)
1129  Smax = std::max(Sf[i], Smax);
1130  for (int i = 0; i < M; i++)
1131  if (Sf[i] < 0.0)
1132  std::cout << "negative singlar value.\n";
1133  // perr << "Smax = " << Smax << std::endl;
1134  Sfinv.resize(Sf.size());
1135  for (int i = 0; i < M; i++)
1136  Sfinv[i] = (Sf[i] < (tolerance * Smax)) ? 0.0 : (1.0 / Sf[i]);
1137  numSingular = 0;
1138  for (int i = 0; i < Sfinv.size(); i++)
1139  if (Sfinv[i] == 0.0)
1140  numSingular++;
1141  if (numSingular > 0)
1142  std::cout << "There were " << numSingular << " singular values in force breakup.\n";
1143  // Compute t_n, removing singular values
1144 
1145  //First, do strain.
1146  Smax = Ss[0];
1147  for (int i = 1; i < M; i++)
1148  Smax = std::max(Ss[i], Smax);
1149  for (int i = 0; i < M; i++)
1150  if (Ss[i] < 0.0)
1151  std::cout << "negative singlar value.\n";
1152  // perr << "Smax = " << Smax << std::endl;
1153  Ssinv.resize(Ss.size());
1154  for (int i = 0; i < M; i++)
1155  Ssinv[i] = (Ss[i] < (tolerance * Smax)) ? 0.0 : (1.0 / Ss[i]);
1156  numSingular = 0;
1157  for (int i = 0; i < Ssinv.size(); i++)
1158  if (Ssinv[i] == 0.0)
1159  numSingular++;
1160  if (numSingular > 0)
1161  std::cout << "There were " << numSingular << " singular values in strain breakup.\n";
1162  // Compute t_n, removing singular values
1163 
1164  for (int i = 0; i < M; i++)
1165  {
1166  mRealType coef = 0.0;
1167  mRealType coef_f = 0.0;
1168  mRealType coef_s = 0.0;
1169 
1170  for (int j = 0; j < M; j++)
1171  {
1172  coef += U(j, i) * bc[j];
1173  coef_f += Uf(j, i) * bfc[j];
1174  coef_s += Us(j, i) * bsc[j];
1175  }
1176 
1177  coef *= Sinv[i];
1178  coef_f *= Sfinv[i];
1179  coef_s *= Ssinv[i];
1180 
1181  for (int k = 0; k < M; k++)
1182  {
1183  tc[k] += coef * V(k, i);
1184  tfc[k] += coef_f * Vf(k, i);
1185  tsc[k] += coef_s * Vs(k, i);
1186  }
1187  }
1188  // Now copy tc values into t
1189  j = 0;
1190  for (int i = 0; i < N; i++)
1191  if (adjust[i])
1192  {
1193  t[i] = tc[j];
1194  gt[i] = tfc[j];
1195  dt[i] = tsc[j];
1196  j++;
1197  }
1198  // Calculate chi-squared
1199  mRealType Yk(0.0), chi2(0.0);
1200  mRealType Yk_f(0.0), chi2_f(0.0);
1201  mRealType Yk_s(0.0), chi2_s(0.0);
1202 
1203  for (int ki = 0; ki < KList.size(); ki++)
1204  {
1205  Yk = Fk[ki]; //-Fk[ki];
1206  Yk_f = Fk[ki];
1207  Yk_s = dFk[ki];
1208 
1209  for (int n = 0; n < N; n++)
1210  {
1211  Yk -= cnk(n, ki) * t[n];
1212  Yk_f -= cnk(n, ki) * gt[n];
1213  Yk_s -= dcnk(n, ki) * dt[n];
1214  }
1215  chi2 += KList[ki][1] * Yk * Yk;
1216  chi2_f += KList[ki][0] * KList[ki][0] * KList[ki][1] * Yk_f * Yk_f;
1217  chi2_s += KList[ki][0] * KList[ki][0] * KList[ki][1] * Yk_s * Yk_s;
1218  }
1219  // std::vector<mRealType> chisqrtmp(3);
1220 
1221  chisqrlist[0] = chi2;
1222  chisqrlist[1] = chi2_f;
1223  chisqrlist[2] = chi2_s;
1224 
1225  //chisqrlist=chisqrtmp;
1226 }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoBreakup() [1/2]

LRBreakup< BreakupBasis >::mRealType DoBreakup ( mRealType *  Fk,
mRealType *  t,
mRealType *  adjust 
)

Definition at line 389 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), qmcplusplus::Units::distance::m, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

Referenced by LRRPABFeeHandlerTemp< Func, BreakupBasis >::InitBreakup(), LRRPAHandlerTemp< Func, BreakupBasis >::InitBreakup(), and LRHandlerTemp< Func, BreakupBasis >::InitBreakup().

392 {
393  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
394  //t and adjust must be allocated up to Basis.NumBasisElem();
395  //Fk must be allocated and filled up to KList.size();
396  // assert(t.size()==adjust.size());
397  // assert(t.size()==Basis.NumBasisElem());
398  Matrix<mRealType> A;
399  std::vector<mRealType> b;
400  Matrix<mRealType> cnk;
401  int N = Basis.NumBasisElem(); //t.size();
402  A.resize(N, N);
403  b.resize(N, 0.0);
404  cnk.resize(N, KList.size());
405  //Fill in cnk.
406  for (int n = 0; n < N; n++)
407  {
408  for (int ki = 0; ki < KList.size(); ki++)
409  {
410  mRealType k = KList[ki][0];
411  cnk(n, ki) = Basis.c(n, k);
412  }
413  }
414  //Fill in A and b
415  A = 0.0;
416  for (int l = 0; l < N; l++)
417  {
418  for (int ki = 0; ki < KList.size(); ki++)
419  {
420  b[l] += KList[ki][1] * Fk[ki] * cnk(l, ki);
421  for (int n = 0; n < N; n++)
422  A(l, n) += KList[ki][1] * cnk(l, ki) * cnk(n, ki);
423  }
424  }
425  //Reduce for constraints
426  int M = N;
427  for (int i = 0; i < N; i++)
428  if (!adjust[i])
429  M--;
430  //The c is for "constrained"
431  Matrix<mRealType> Ac;
432  Ac.resize(M, M);
433  std::vector<mRealType> bc(M, 0.0), tc(M, 0.0);
434  //Build constrained Ac and bc
435  int j = 0;
436  for (int col = 0; col < N; col++)
437  {
438  if (adjust[col])
439  {
440  // Copy column a A to Ac
441  int i = 0;
442  for (int row = 0; row < N; row++)
443  if (adjust[row])
444  {
445  Ac(i, j) = A(row, col);
446  i++;
447  }
448  j++;
449  }
450  else
451  {
452  // Otherwise, subtract t(col)*A(:,col) from bc
453  for (int row = 0; row < N; row++)
454  b[row] -= A(row, col) * t[col];
455  }
456  }
457  j = 0;
458  for (int row = 0; row < N; row++)
459  if (adjust[row])
460  {
461  bc[j] = b[row];
462  j++;
463  }
464  // Do SVD:
465  // -------
466  // Matrix<mRealType> U(M, M), V(M, M);
467  // std::vector<mRealType> S(M), Sinv(M);
468  // SVdecomp(Ac, U, S, V);
469  ////////////////////////////////
470  int m = Ac.rows();
471  int n = Ac.cols();
472  Matrix<mRealType> Atrans(n, m);
473  Matrix<mRealType> U, V;
474  U.resize(std::min(m, n), m);
475  V.resize(n, std::min(m, n));
476  std::vector<mRealType> S, Sinv;
477  S.resize(std::min(n, m));
478  //do the transpose
479  for (int i = 0; i < m; i++)
480  {
481  for (int j = 0; j < n; j++)
482  Atrans(j, i) = Ac(i, j);
483  }
484  char JOBU = 'S';
485  char JOBVT = 'S';
486  int LDA = m;
487  int LDU = m;
488  int LDVT = std::min(m, n);
489  int LWORK = 10 * std::max(3 * std::min(n, m) + std::max(m, n), 5 * std::min(m, n));
490  std::vector<mRealType> WORK(LWORK);
491  int INFO;
492  LAPACK::gesvd(JOBU, JOBVT, m, n, Atrans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
493  assert(INFO == 0);
494  int ur = U.rows();
495  int uc = U.cols();
496  Matrix<mRealType> Utrans(uc, ur);
497  for (int i = 0; i < ur; i++)
498  {
499  for (int j = 0; j < uc; j++)
500  Utrans(j, i) = U(i, j);
501  }
502  U.resize(uc, ur);
503  U = Utrans;
504  //////////////////////////////////
505  // Zero out near-singular values
506  mRealType Smax = S[0];
507  for (int i = 1; i < M; i++)
508  Smax = std::max(S[i], Smax);
509  for (int i = 0; i < M; i++)
510  if (S[i] < 0.0)
511  std::cout << "negative singlar value.\n";
512  // perr << "Smax = " << Smax << std::endl;
513  Sinv.resize(S.size());
514  for (int i = 0; i < M; i++)
515  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
516  int numSingular = 0;
517  for (int i = 0; i < Sinv.size(); i++)
518  if (Sinv[i] == 0.0)
519  numSingular++;
520  if (numSingular > 0)
521  std::cout << "There were " << numSingular << " singular values in breakup.\n";
522  // Compute t_n, removing singular values
523  for (int i = 0; i < M; i++)
524  {
525  mRealType coef = 0.0;
526  for (int j = 0; j < M; j++)
527  coef += U(j, i) * bc[j];
528  coef *= Sinv[i];
529  for (int k = 0; k < M; k++)
530  tc[k] += coef * V(k, i);
531  }
532  // Now copy tc values into t
533  j = 0;
534  for (int i = 0; i < N; i++)
535  if (adjust[i])
536  {
537  t[i] = tc[j];
538  j++;
539  }
540  // Calculate chi-squared
541  mRealType Yk, chi2;
542  chi2 = 0.0;
543  for (int ki = 0; ki < KList.size(); ki++)
544  {
545  Yk = Fk[ki];
546  for (int n = 0; n < N; n++)
547  {
548  Yk -= cnk(n, ki) * t[n];
549  }
550  chi2 += KList[ki][1] * Yk * Yk;
551  }
552  return (chi2);
553 }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoBreakup() [2/2]

mRealType DoBreakup ( mRealType *  Fk,
mRealType *  t 
)
inline

Definition at line 71 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, LRBreakup< BreakupBasis >::Basis, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), LRBreakup< BreakupBasis >::KList, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

72  {
73  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
74  //t must be allocated up to Basis.NumBasisElem();
75  //Fk must be allocated and filled up to KList.size();
76  // assert(t.size()==Basis.NumBasisElem());
77  Matrix<mRealType> A;
78  std::vector<mRealType> b;
79  Matrix<mRealType> cnk;
80  int numElem = Basis.NumBasisElem(); //t.size();
81  A.resize(numElem, numElem);
82  b.resize(numElem, 0.0);
83  cnk.resize(numElem, KList.size());
84 // Fill in cnk.
85 // app_log() << "Check OMP size : numElem, KList.size : " << numElem << " , " << KList.size() << std::endl;
86 #pragma omp parallel for shared(cnk)
87  for (int n = 0; n < numElem; n++)
88  {
89  for (int ki = 0; ki < KList.size(); ki++)
90  {
91  mRealType k = KList[ki][0];
92  cnk(n, ki) = Basis.c(n, k);
93  }
94  }
95  // Now, fill in A and b
96  A = 0.0;
97  for (int l = 0; l < numElem; l++)
98  {
99  for (int ki = 0; ki < KList.size(); ki++)
100  {
101  b[l] += KList[ki][1] * Fk[ki] * cnk(l, ki);
102  for (int n = 0; n < numElem; n++)
103  A(l, n) += KList[ki][1] * cnk(l, ki) * cnk(n, ki);
104  }
105  }
106  //////////////////////////
107  //Do the SVD:
108  // Matrix<mRealType> U(numElem, numElem), V(numElem, numElem);
109  // std::vector<mRealType> S(numElem), Sinv(numElem);
110  //////////////////////////
111  // SVdecomp(A, U, S, V);
112  //////////////////////////
113  int M = A.rows();
114  int N = A.cols();
115  Matrix<mRealType> Atrans(N, M);
116  Matrix<mRealType> U, V;
117  U.resize(std::min(M, N), M);
118  V.resize(N, std::min(M, N));
119  std::vector<mRealType> S, Sinv;
120  S.resize(std::min(N, M));
121  //Do the transpose
122  for (int i = 0; i < M; i++)
123  {
124  for (int j = 0; j < N; j++)
125  Atrans(j, i) = A(i, j);
126  }
127  char JOBU = 'S';
128  char JOBVT = 'S';
129  int LDA = M;
130  int LDU = M;
131  int LDVT = std::min(M, N);
132  int LWORK = 10 * std::max(3 * std::min(N, M) + std::max(M, N), 5 * std::min(M, N));
133  std::vector<mRealType> WORK(LWORK);
134  int INFO;
135  LAPACK::gesvd(JOBU, JOBVT, M, N, Atrans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
136  assert(INFO == 0);
137  int ur = U.rows();
138  int uc = U.cols();
139  Matrix<mRealType> Utrans(uc, ur);
140  for (int i = 0; i < ur; i++)
141  {
142  for (int j = 0; j < uc; j++)
143  Utrans(j, i) = U(i, j);
144  }
145  U.resize(uc, ur);
146  U = Utrans;
147  ///////////////////////////////////
148  // Zero out near-singular values
149  mRealType Smax = S[0];
150  for (int i = 1; i < S.size(); i++)
151  Smax = std::max(S[i], Smax);
152  Sinv.resize(S.size());
153  for (int i = 0; i < S.size(); i++)
154  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
155  int numSingular = 0;
156  for (int i = 0; i < Sinv.size(); i++)
157  if (Sinv[i] == 0.0)
158  numSingular++;
159  if (numSingular > 0)
160  std::cout << "There were " << numSingular << " singular values in breakup.\n";
161  for (int i = 0; i < numElem; i++)
162  t[i] = 0.0;
163  // Compute t_n, removing singular values
164  for (int i = 0; i < numElem; i++)
165  {
166  mRealType coef = 0.0;
167  for (int j = 0; j < numElem; j++)
168  coef += U(j, i) * b[j];
169  coef *= Sinv[i];
170  for (int k = 0; k < numElem; k++)
171  t[k] += coef * V(k, i);
172  }
173  // Calculate chi-squared
174  mRealType Yk, chi2;
175  chi2 = 0.0;
176  for (int ki = 0; ki < KList.size(); ki++)
177  {
178  Yk = Fk[ki];
179  for (int n = 0; n < numElem; n++)
180  {
181  Yk -= cnk(n, ki) * t[n];
182  }
183  chi2 += KList[ki][1] * Yk * Yk;
184  }
185  return (chi2);
186  }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoGradBreakup() [1/2]

LRBreakup< BreakupBasis >::mRealType DoGradBreakup ( mRealType *  Fk,
mRealType *  t,
mRealType *  adjust 
)

Definition at line 557 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), qmcplusplus::Units::distance::m, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

Referenced by LRHandlerSRCoulomb< Func, BreakupBasis >::InitBreakup().

560 {
561  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
562  //t and adjust must be allocated up to Basis.NumBasisElem();
563  //Fk must be allocated and filled up to KList.size();
564  // assert(t.size()==adjust.size());
565  // assert(t.size()==Basis.NumBasisElem());
566  Matrix<mRealType> A;
567  std::vector<mRealType> b;
568  Matrix<mRealType> cnk;
569  int N = Basis.NumBasisElem(); //t.size();
570  A.resize(N, N);
571  b.resize(N, 0.0);
572  cnk.resize(N, KList.size());
573  //Fill in cnk.
574  for (int n = 0; n < N; n++)
575  {
576  for (int ki = 0; ki < KList.size(); ki++)
577  {
578  mRealType k = KList[ki][0];
579  cnk(n, ki) = Basis.c(n, k);
580  }
581  }
582  //Fill in A and b
583  A = 0.0;
584  for (int l = 0; l < N; l++)
585  {
586  for (int ki = 0; ki < KList.size(); ki++)
587  {
588  mRealType k2 = KList[ki][0] * KList[ki][0];
589  b[l] += k2 * KList[ki][1] * Fk[ki] * cnk(l, ki);
590  for (int n = 0; n < N; n++)
591  A(l, n) += k2 * KList[ki][1] * cnk(l, ki) * cnk(n, ki);
592  }
593  }
594  //Reduce for constraints
595  int M = N;
596  for (int i = 0; i < N; i++)
597  if (!adjust[i])
598  M--;
599  //The c is for "constrained"
600  Matrix<mRealType> Ac;
601  Ac.resize(M, M);
602  std::vector<mRealType> bc(M, 0.0), tc(M, 0.0);
603  //Build constrained Ac and bc
604  int j = 0;
605  for (int col = 0; col < N; col++)
606  {
607  if (adjust[col])
608  {
609  // Copy column a A to Ac
610  int i = 0;
611  for (int row = 0; row < N; row++)
612  if (adjust[row])
613  {
614  Ac(i, j) = A(row, col);
615  i++;
616  }
617  j++;
618  }
619  else
620  {
621  // Otherwise, subtract t(col)*A(:,col) from bc
622  for (int row = 0; row < N; row++)
623  b[row] -= A(row, col) * t[col];
624  }
625  }
626  j = 0;
627  for (int row = 0; row < N; row++)
628  if (adjust[row])
629  {
630  bc[j] = b[row];
631  j++;
632  }
633  // Do SVD:
634  // -------
635  // Matrix<mRealType> U(M, M), V(M, M);
636  // std::vector<mRealType> S(M), Sinv(M);
637  // SVdecomp(Ac, U, S, V);
638  ////////////////////////////////
639  int m = Ac.rows();
640  int n = Ac.cols();
641  Matrix<mRealType> Atrans(n, m);
642  Matrix<mRealType> U, V;
643  U.resize(std::min(m, n), m);
644  V.resize(n, std::min(m, n));
645  std::vector<mRealType> S, Sinv;
646  S.resize(std::min(n, m));
647  //do the transpose
648  for (int i = 0; i < m; i++)
649  {
650  for (int j = 0; j < n; j++)
651  Atrans(j, i) = Ac(i, j);
652  }
653  char JOBU = 'S';
654  char JOBVT = 'S';
655  int LDA = m;
656  int LDU = m;
657  int LDVT = std::min(m, n);
658  int LWORK = 10 * std::max(3 * std::min(n, m) + std::max(m, n), 5 * std::min(m, n));
659  std::vector<mRealType> WORK(LWORK);
660  int INFO;
661  LAPACK::gesvd(JOBU, JOBVT, m, n, Atrans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
662  assert(INFO == 0);
663  int ur = U.rows();
664  int uc = U.cols();
665  Matrix<mRealType> Utrans(uc, ur);
666  for (int i = 0; i < ur; i++)
667  {
668  for (int j = 0; j < uc; j++)
669  Utrans(j, i) = U(i, j);
670  }
671  U.resize(uc, ur);
672  U = Utrans;
673  //////////////////////////////////
674  // Zero out near-singular values
675  mRealType Smax = S[0];
676  for (int i = 1; i < M; i++)
677  Smax = std::max(S[i], Smax);
678  for (int i = 0; i < M; i++)
679  if (S[i] < 0.0)
680  std::cout << "negative singlar value.\n";
681  // perr << "Smax = " << Smax << std::endl;
682  Sinv.resize(S.size());
683  for (int i = 0; i < M; i++)
684  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
685  int numSingular = 0;
686  for (int i = 0; i < Sinv.size(); i++)
687  if (Sinv[i] == 0.0)
688  numSingular++;
689  if (numSingular > 0)
690  std::cout << "There were " << numSingular << " singular values in breakup.\n";
691  // Compute t_n, removing singular values
692  for (int i = 0; i < M; i++)
693  {
694  mRealType coef = 0.0;
695  for (int j = 0; j < M; j++)
696  coef += U(j, i) * bc[j];
697  coef *= Sinv[i];
698  for (int k = 0; k < M; k++)
699  tc[k] += coef * V(k, i);
700  }
701  // Now copy tc values into t
702  j = 0;
703  for (int i = 0; i < N; i++)
704  if (adjust[i])
705  {
706  t[i] = tc[j];
707  j++;
708  }
709  // Calculate chi-squared
710  mRealType Yk, chi2;
711  chi2 = 0.0;
712  for (int ki = 0; ki < KList.size(); ki++)
713  {
714  Yk = Fk[ki];
715  for (int n = 0; n < N; n++)
716  {
717  Yk -= cnk(n, ki) * t[n];
718  }
719  chi2 += KList[ki][0] * KList[ki][0] * KList[ki][1] * Yk * Yk;
720  }
721  return (chi2);
722 }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoGradBreakup() [2/2]

mRealType DoGradBreakup ( mRealType *  Fk,
mRealType *  t 
)
inline

Definition at line 196 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, LRBreakup< BreakupBasis >::Basis, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), LRBreakup< BreakupBasis >::KList, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

197  {
198  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
199  //t must be allocated up to Basis.NumBasisElem();
200  //Fk must be allocated and filled up to KList.size();
201  // assert(t.size()==Basis.NumBasisElem());
202  Matrix<mRealType> A;
203  std::vector<mRealType> b;
204  Matrix<mRealType> cnk;
205  int numElem = Basis.NumBasisElem(); //t.size();
206  A.resize(numElem, numElem);
207  b.resize(numElem, 0.0);
208  cnk.resize(numElem, KList.size());
209  // Fill in cnk.
210  for (int n = 0; n < numElem; n++)
211  {
212  for (int ki = 0; ki < KList.size(); ki++)
213  {
214  mRealType k = KList[ki][0];
215  cnk(n, ki) = Basis.c(n, k);
216  }
217  }
218  // Now, fill in A and b
219 
220  A = 0.0;
221  for (int l = 0; l < numElem; l++)
222  {
223  for (int ki = 0; ki < KList.size(); ki++)
224  {
225  mRealType k2 = KList[ki][0] * KList[ki][0];
226  b[l] += k2 * KList[ki][1] * Fk[ki] * cnk(l, ki);
227  for (int n = 0; n < numElem; n++)
228  A(l, n) += k2 * KList[ki][1] * cnk(l, ki) * cnk(n, ki);
229  }
230  }
231  //////////////////////////
232  //Do the SVD:
233  // Matrix<mRealType> U(numElem, numElem), V(numElem, numElem);
234  // std::vector<mRealType> S(numElem), Sinv(numElem);
235  //////////////////////////
236  // SVdecomp(A, U, S, V);
237  //////////////////////////
238  int M = A.rows();
239  int N = A.cols();
240  Matrix<mRealType> Atrans(N, M);
241  Matrix<mRealType> U, V;
242  U.resize(std::min(M, N), M);
243  V.resize(N, std::min(M, N));
244  std::vector<mRealType> S, Sinv;
245  S.resize(std::min(N, M));
246  //Do the transpose
247  for (int i = 0; i < M; i++)
248  {
249  for (int j = 0; j < N; j++)
250  Atrans(j, i) = A(i, j);
251  }
252  char JOBU = 'S';
253  char JOBVT = 'S';
254  int LDA = M;
255  int LDU = M;
256  int LDVT = std::min(M, N);
257  int LWORK = 10 * std::max(3 * std::min(N, M) + std::max(M, N), 5 * std::min(M, N));
258  std::vector<mRealType> WORK(LWORK);
259  int INFO;
260  LAPACK::gesvd(JOBU, JOBVT, M, N, Atrans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
261  assert(INFO == 0);
262  int ur = U.rows();
263  int uc = U.cols();
264  Matrix<mRealType> Utrans(uc, ur);
265  for (int i = 0; i < ur; i++)
266  {
267  for (int j = 0; j < uc; j++)
268  Utrans(j, i) = U(i, j);
269  }
270  U.resize(uc, ur);
271  U = Utrans;
272  ///////////////////////////////////
273  // Zero out near-singular values
274  mRealType Smax = S[0];
275  for (int i = 1; i < S.size(); i++)
276  Smax = std::max(S[i], Smax);
277  Sinv.resize(S.size());
278  for (int i = 0; i < S.size(); i++)
279  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
280  int numSingular = 0;
281  for (int i = 0; i < Sinv.size(); i++)
282  if (Sinv[i] == 0.0)
283  numSingular++;
284  if (numSingular > 0)
285  std::cout << "There were " << numSingular << " singular values in breakup.\n";
286  for (int i = 0; i < numElem; i++)
287  t[i] = 0.0;
288  // Compute t_n, removing singular values
289  for (int i = 0; i < numElem; i++)
290  {
291  mRealType coef = 0.0;
292  for (int j = 0; j < numElem; j++)
293  coef += U(j, i) * b[j];
294  coef *= Sinv[i];
295  for (int k = 0; k < numElem; k++)
296  t[k] += coef * V(k, i);
297  }
298  // Calculate chi-squared
299  mRealType Yk, chi2;
300  chi2 = 0.0;
301  for (int ki = 0; ki < KList.size(); ki++)
302  {
303  mRealType k2 = KList[ki][0] * KList[ki][0];
304  Yk = Fk[ki];
305  for (int n = 0; n < numElem; n++)
306  {
307  Yk -= cnk(n, ki) * t[n];
308  }
309  chi2 += k2 * KList[ki][1] * Yk * Yk;
310  }
311  return (chi2);
312  }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ DoStrainBreakup()

LRBreakup< BreakupBasis >::mRealType DoStrainBreakup ( mRealType *  Fk,
mRealType *  dFk,
mRealType *  t,
mRealType *  adjust 
)

Definition at line 724 of file LRBreakup.h.

References qmcplusplus::Units::distance::A, Matrix< T, Alloc >::cols(), Matrix< T, Alloc >::data(), LAPACK::gesvd(), qmcplusplus::Units::distance::m, omptarget::min(), qmcplusplus::Units::force::N, qmcplusplus::n, Matrix< T, Alloc >::resize(), and Matrix< T, Alloc >::rows().

728 {
729  const mRealType tolerance = std::numeric_limits<mRealType>::epsilon();
730  //t and adjust must be allocated up to Basis.NumBasisElem();
731  //Fk must be allocated and filled up to KList.size();
732  // assert(t.size()==adjust.size());
733  // assert(t.size()==Basis.NumBasisElem());
734  Matrix<mRealType> A;
735  std::vector<mRealType> b;
736  Matrix<mRealType> dcnk;
737  int N = Basis.NumBasisElem(); //t.size();
738  A.resize(N, N);
739  b.resize(N, 0.0);
740  dcnk.resize(N, KList.size());
741  //Fill in cnk.
742  for (int n = 0; n < N; n++)
743  {
744  for (int ki = 0; ki < KList.size(); ki++)
745  {
746  mRealType k = KList[ki][0];
747  dcnk(n, ki) = Basis.dc_dk(n, k); //-Basis.c(n,k);
748  }
749  }
750  //Fill in A and b
751  A = 0.0;
752  for (int l = 0; l < N; l++)
753  {
754  for (int ki = 0; ki < KList.size(); ki++)
755  {
756  mRealType k2 = KList[ki][0] * KList[ki][0];
757  // b[l] += k2*KList[ki][1]*(dFk[ki]-Fk[ki]) * dcnk(l, ki);
758  b[l] += k2 * KList[ki][1] * (dFk[ki]) * dcnk(l, ki);
759  for (int n = 0; n < N; n++)
760  A(l, n) += k2 * KList[ki][1] * dcnk(l, ki) * dcnk(n, ki);
761  }
762  }
763  //Reduce for constraints
764  int M = N;
765  for (int i = 0; i < N; i++)
766  if (!adjust[i])
767  M--;
768  //The c is for "constrained"
769  Matrix<mRealType> Ac;
770  Ac.resize(M, M);
771  std::vector<mRealType> bc(M, 0.0), tc(M, 0.0);
772  //Build constrained Ac and bc
773  int j = 0;
774  for (int col = 0; col < N; col++)
775  {
776  if (adjust[col])
777  {
778  // Copy column a A to Ac
779  int i = 0;
780  for (int row = 0; row < N; row++)
781  if (adjust[row])
782  {
783  Ac(i, j) = A(row, col);
784  i++;
785  }
786  j++;
787  }
788  else
789  {
790  // Otherwise, subtract t(col)*A(:,col) from bc
791  for (int row = 0; row < N; row++)
792  b[row] -= A(row, col) * t[col];
793  }
794  }
795  j = 0;
796  for (int row = 0; row < N; row++)
797  if (adjust[row])
798  {
799  bc[j] = b[row];
800  j++;
801  }
802  // Do SVD:
803  // -------
804  // Matrix<mRealType> U(M, M), V(M, M);
805  // std::vector<mRealType> S(M), Sinv(M);
806  // SVdecomp(Ac, U, S, V);
807  ////////////////////////////////
808  int m = Ac.rows();
809  int n = Ac.cols();
810  Matrix<mRealType> Atrans(n, m);
811  Matrix<mRealType> U, V;
812  U.resize(std::min(m, n), m);
813  V.resize(n, std::min(m, n));
814  std::vector<mRealType> S, Sinv;
815  S.resize(std::min(n, m));
816  //do the transpose
817  for (int i = 0; i < m; i++)
818  {
819  for (int j = 0; j < n; j++)
820  Atrans(j, i) = Ac(i, j);
821  }
822  char JOBU = 'S';
823  char JOBVT = 'S';
824  int LDA = m;
825  int LDU = m;
826  int LDVT = std::min(m, n);
827  int LWORK = 10 * std::max(3 * std::min(n, m) + std::max(m, n), 5 * std::min(m, n));
828  std::vector<mRealType> WORK(LWORK);
829  int INFO;
830  LAPACK::gesvd(JOBU, JOBVT, m, n, Atrans.data(), LDA, &S[0], U.data(), LDU, V.data(), LDVT, &WORK[0], LWORK, INFO);
831  assert(INFO == 0);
832  int ur = U.rows();
833  int uc = U.cols();
834  Matrix<mRealType> Utrans(uc, ur);
835  for (int i = 0; i < ur; i++)
836  {
837  for (int j = 0; j < uc; j++)
838  Utrans(j, i) = U(i, j);
839  }
840  U.resize(uc, ur);
841  U = Utrans;
842  //////////////////////////////////
843  // Zero out near-singular values
844  mRealType Smax = S[0];
845  for (int i = 1; i < M; i++)
846  Smax = std::max(S[i], Smax);
847  for (int i = 0; i < M; i++)
848  if (S[i] < 0.0)
849  std::cout << "negative singlar value.\n";
850  // perr << "Smax = " << Smax << std::endl;
851  Sinv.resize(S.size());
852  for (int i = 0; i < M; i++)
853  Sinv[i] = (S[i] < (tolerance * Smax)) ? 0.0 : (1.0 / S[i]);
854  int numSingular = 0;
855  for (int i = 0; i < Sinv.size(); i++)
856  if (Sinv[i] == 0.0)
857  numSingular++;
858  if (numSingular > 0)
859  std::cout << "There were " << numSingular << " singular values in breakup.\n";
860  // Compute t_n, removing singular values
861  for (int i = 0; i < M; i++)
862  {
863  mRealType coef = 0.0;
864  for (int j = 0; j < M; j++)
865  coef += U(j, i) * bc[j];
866  coef *= Sinv[i];
867  for (int k = 0; k < M; k++)
868  tc[k] += coef * V(k, i);
869  }
870  // Now copy tc values into t
871  j = 0;
872  for (int i = 0; i < N; i++)
873  if (adjust[i])
874  {
875  t[i] = tc[j];
876  j++;
877  }
878  // Calculate chi-squared
879  mRealType Yk, chi2;
880  chi2 = 0.0;
881  for (int ki = 0; ki < KList.size(); ki++)
882  {
883  Yk = dFk[ki]; //-Fk[ki];
884  for (int n = 0; n < N; n++)
885  {
886  Yk -= dcnk(n, ki) * t[n];
887  }
888  chi2 += KList[ki][0] * KList[ki][0] * KList[ki][1] * Yk * Yk;
889  }
890  return (chi2);
891 }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
omptarget::min
T min(T a, T b)
Definition: OMPTargetMath.hpp:36
LAPACK::gesvd
static void gesvd(const char &jobu, const char &jobvt, const int &m, const int &n, float *a, const int &lda, float *s, float *u, const int &ldu, float *vt, const int &ldvt, float *work, const int &lwork, int &info)
Definition: BLAS.hpp:520
qmcplusplus::LRBreakup::KList
std::vector< TinyVector< mRealType, 2 > > KList
For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.
Definition: LRBreakup.h:44

◆ SetupKVecs()

int SetupKVecs ( mRealType  kc,
mRealType  kcont,
mRealType  kmax 
)

setup KList

Parameters
kck-space cutoff for long-range sums
kcontk at which approximate (spherical shell) degeneracies are used.
kmaxlargest k used for performing the breakup
Returns
the maximum kshell for the given kc

Definition at line 333 of file LRBreakup.h.

References qmcplusplus::app_log(), KContainer::kshell, KContainer::ksq, qmcplusplus::Units::force::N, KContainer::numk, qmcplusplus::sqrt(), and KContainer::updateKLists().

Referenced by LRRPABFeeHandlerTemp< Func, BreakupBasis >::InitBreakup(), LRRPAHandlerTemp< Func, BreakupBasis >::InitBreakup(), LRHandlerTemp< Func, BreakupBasis >::InitBreakup(), and LRHandlerSRCoulomb< Func, BreakupBasis >::InitBreakup().

334 {
335  //Add low |k| ( < kcont) k-points with exact degeneracy
336  KContainer kexact;
337  kexact.updateKLists(Basis.get_Lattice(), kcont, Basis.get_Lattice().ndim);
338  bool findK = true;
339  mRealType kc2 = kc * kc;
340  //use at least one shell
341  size_t ks = 0;
342  kc2 = std::max(kc2, static_cast<mRealType>(kexact.ksq[kexact.kshell[ks]]));
343  while (findK)
344  {
345  if (kexact.ksq[kexact.kshell[ks]] > kc2)
346  findK = false;
347  else
348  ks++;
349  }
350  size_t maxkshell = ks;
351  size_t numk = kexact.numk - kexact.kshell[ks];
352  for (; ks < kexact.kshell.size() - 1; ks++)
353  AddKToList(std::sqrt(kexact.ksq[kexact.kshell[ks]]), kexact.kshell[ks + 1] - kexact.kshell[ks]);
354  ////Add these vectors to the internal list
355  //int numk=0;
356  //mRealType modk2;
357  //for(int ki=0; ki<kexact.numk; ki++) {
358  // modk2 = dot(kexact.kpts_cart[ki],kexact.kpts_cart[ki]);
359  // if(modk2 > (kc*kc)) { //Breakup needs kc < k < kcont.
360  // AddKToList(std::sqrt(modk2));
361  // numk++;
362  // }
363  //}
364  //Add high |k| ( >kcont, <kmax) k-points with approximate degeneracy
365  //Volume of 1 K-point is (2pi)^3/(a1.a2^a3)
366  mRealType kelemvol = 8 * M_PI * M_PI * M_PI / Basis.get_CellVolume();
367  //Generate 4000 shells:
368  const int N = 4000;
369  mRealType deltak = (kmax - kcont) / N;
370  for (int i = 0; i < N; i++)
371  {
372  mRealType k1 = kcont + deltak * i;
373  mRealType k2 = k1 + deltak;
374  mRealType kmid = 0.5 * (k1 + k2);
375  mRealType shellvol = 4.0 * M_PI * (k2 * k2 * k2 - k1 * k1 * k1) / 3.0;
376  mRealType degeneracy = shellvol / kelemvol;
377  AddKToList(kmid, degeneracy);
378  numk += static_cast<int>(degeneracy);
379  }
380  app_log() << " NUMBER OF OPT_BREAK KVECS = " << numk << std::endl;
381 
382  return maxkshell;
383  //numk now contains the total number of vectors.
384  //this->klist.size() contains the number of unique vectors.
385 }
qmcplusplus::LRBreakup::Basis
BreakupBasis & Basis
The basis to be used for breakup.
Definition: LRBreakup.h:42
qmcplusplus::app_log
std::ostream & app_log()
Definition: OutputManager.h:65
qmcplusplus::mRealType
EwaldHandler3D::mRealType mRealType
Definition: test_ewald3d.cpp:21
qmcplusplus::LRBreakup::AddKToList
void AddKToList(mRealType k, mRealType degeneracy=1.0)
Definition: LRBreakup.h:316
qmcplusplus::sqrt
MakeReturn< UnaryNode< FnSqrt, typename CreateLeaf< Vector< T1, C1 > >::Leaf_t > >::Expression_t sqrt(const Vector< T1, C1 > &l)
Definition: OhmmsVectorOperators.h:136

Member Data Documentation

◆ Basis

BreakupBasis& Basis

The basis to be used for breakup.

Definition at line 42 of file LRBreakup.h.

Referenced by LRBreakup< BreakupBasis >::DoBreakup(), and LRBreakup< BreakupBasis >::DoGradBreakup().

◆ KList

std::vector<TinyVector<mRealType, 2> > KList

For each k, KList[k][0] = |k| and KList[k][1] = degeneracy.

Definition at line 44 of file LRBreakup.h.

Referenced by LRBreakup< BreakupBasis >::DoBreakup(), LRBreakup< BreakupBasis >::DoGradBreakup(), LRRPABFeeHandlerTemp< Func, BreakupBasis >::InitBreakup(), LRRPAHandlerTemp< Func, BreakupBasis >::InitBreakup(), LRHandlerTemp< Func, BreakupBasis >::InitBreakup(), and LRHandlerSRCoulomb< Func, BreakupBasis >::InitBreakup().

The documentation for this struct was generated from the following file: