MADNESS  0.10.1
solver.h
Go to the documentation of this file.
1 /*
2  This file is part of MADNESS.
3 
4  Copyright (C) 2007,2010 Oak Ridge National Laboratory
5 
6  This program is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation; either version 2 of the License, or
9  (at your option) any later version.
10 
11  This program is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with this program; if not, write to the Free Software
18  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 
20  For more information please contact:
21 
22  Robert J. Harrison
23  Oak Ridge National Laboratory
24  One Bethel Valley Road
25  P.O. Box 2008, MS-6367
26 
27  email: harrisonrj@ornl.gov
28  tel: 865-241-3937
29  fax: 865-572-0680
30 */
31 
32 #include <madness/mra/mra.h>
33 //#include <madness/mra/lbdeux.h>
34 #include <madness/world/MADworld.h>
35 #include <madness/tensor/solvers.h>
36 #include <vector>
37 #include <madness/fortran_ctypes.h>
38 #include <cmath>
39 
40 #include "poperator.h"
41 #include "libxc.h"
43 #include "complexfun.h"
44 #include "esolver.h"
45 
46 #ifndef SOLVER_H_
47 
48 using std::endl;
49 using std::setfill;
50 using std::setw;
51 using std::ostream;
52 
53 /*!
54  \ingroup applications
55  \defgroup periodic_solver Periodic Solver
56  \brief The Periodic Solver group is a group that contains the software
57  objects that are needed to solve a periodic Kohn-Sham hamiltonian.
58 */
59 
60 //*****************************************************************************
61 /* static double onesfunc(const coordT& x) */
62 /* { */
63 /* return 1.0; */
64 /* } */
65 //*****************************************************************************
66 
67 namespace madness
68 {
69  //***************************************************************************
70  /*!
71  \ingroup periodic_solver
72  \Fermi-Dirac distribution function for fixing the occupation numbers
73  */
74  template <typename T>
75  T stheta_fd(const T& x)
76  {
77  if (x > 50.0)
78  {
79  return 1.0;
80  }
81  else if (x < -50.0)
82  {
83  return 0.0;
84  }
85  else
86  {
87  return 1.0/(1.0 + exp(-x));
88  }
89  }
90  //***************************************************************************
91 
92  // NOTE: so this is totally hacked because of the requirement of using a box
93  // of L*L*L rather than the convention of using a combination of
94  // lattice vectors and reciprocal lattice vectors (FOR NOW)
95  // On second thought, it might be ok for general use if used "right".
96  //***************************************************************************
97  template <std::size_t NDIM>
98  class ComplexExp : public FunctionFunctorInterface<double_complex,NDIM> {
99  public:
104 
106  : exponent(exponent), coeff(coeff) {}
107 
108  double_complex operator()(const coordT& x) const {
109  double sum = 0.0;
110  for (std::size_t i=0; i<NDIM; ++i) {
111  sum += x[i]*exponent[i];
112  };
113  return coeff*exp(double_complex(0.0,sum));
114  }
115  };
116  //***************************************************************************
117 
118 
119  //***************************************************************************
120  class WSTAtomicBasisFunctor : public FunctionFunctorInterface<double_complex,3> {
121  private:
123  double L;
124  double kx, ky, kz;
125 
126  std::vector<coord_3d> specialpt;
127  public:
128  WSTAtomicBasisFunctor(const AtomicBasisFunction& aofunc, double L, double kx, double ky, double kz)
129  : aofunc(aofunc), L(L), kx(kx), ky(ky), kz(kz)
130  {
131  double x, y, z;
132  aofunc.get_coords(x,y,z);
133  coord_3d r;
134  r[0]=x; r[1]=y; r[2]=z;
135  specialpt=std::vector<coord_3d>(1,r);
136  }
137 
139 
141  double_complex sum = 0.0;
142  const int R = 1;
143  const double_complex I = double_complex(0.0,1.0);
144  for (int i=-R; i<=+R; i++) {
145  double xx = x[0]+i*L;
146  for (int j=-R; j<=+R; j++) {
147  double yy = x[1]+j*L;
148  for (int k=-R; k<=+R; k++) {
149  double zz = x[2]+k*L;
150  sum += -exp(-I*(kx*xx+ky*yy+kz*zz+constants::pi/2))*aofunc(xx, yy, zz);
151  //print(kx,ky,kz,exp(-I*(kx*xx+ky*yy+kz*zz+constants::pi/2)),aofunc(xx, yy, zz));
152  }
153  }
154  }
155  return sum;
156  }
157 
158  std::vector<coord_3d> special_points() const {return specialpt;}
159  };
160  //***************************************************************************
161 
162 
163  /// A MADNESS functor to compute either x, y, or z
164  //***************************************************************************
165  class DipoleFunctor : public FunctionFunctorInterface<double,3> {
166  private:
167  const int axis;
168  public:
170  double operator()(const coordT& x) const {
171  return x[axis];
172  }
173  virtual ~DipoleFunctor() {}
174  };
175  //***************************************************************************
176 
177  /*!
178  \ingroup periodic_solver
179 
180  \brief The SubspaceK class is a container class holding previous orbitals
181  and residuals.
182  \par
183  The Solver class uses the Krylov Accelerated Inexact Newton Solver (KAIN)
184  accelerate the convergence a given calculation. The KAIN solver needs to
185  store a subspace of previous orbitals and residuals. In the case is this
186  implementation, the orbitals are store according to which k-point to which
187  they belong.
188  */
189 
190  //***************************************************************************
191  template <typename T, int NDIM>
192  class SubspaceK
193  {
194  // Typedef's
195  typedef std::complex<T> valueT;
197  typedef std::vector<functionT> vecfuncT;
198  typedef std::pair<vecfuncT,vecfuncT> pairvecfuncT;
199  typedef std::vector<pairvecfuncT> subspaceT;
201  typedef std::vector<tensorT> vectensorT;
202  typedef std::vector<subspaceT> vecsubspaceT;
203 
204  //*************************************************************************
206  //*************************************************************************
207 
208  //*************************************************************************
210  //*************************************************************************
211 
212  //*************************************************************************
214  //*************************************************************************
215 
216  //*************************************************************************
217  std::vector<KPoint> _kpoints;
218  //*************************************************************************
219 
220  //*************************************************************************
222  //*************************************************************************
223 
224  //*************************************************************************
225  double _residual;
226  //*************************************************************************
227 
228  public:
229 
230  //*************************************************************************
232  const std::vector<KPoint>& kpoints) : _world(world), _kpoints(kpoints),
233  _params(params)
234  {
235  _residual = 1e6;
236  for (unsigned int kp = 0; kp < _kpoints.size(); kp++)
237  {
238  _Q.push_back(tensorT(1,1));
239  _subspace.push_back(subspaceT());
240  }
241  }
242  //*************************************************************************
243 
244  //*************************************************************************
245  void update_subspace(vecfuncT& awfs_new,
246  vecfuncT& bwfs_new,
247  const vecfuncT& awfs_old,
248  const vecfuncT& bwfs_old)
249  {
250  // compute residuals (total)
251  vecfuncT t_rm = sub(_world, awfs_old, awfs_new);
252  if (_params.spinpol)
253  {
254  vecfuncT br = sub(_world, bwfs_old, bwfs_new);
255  t_rm.insert(t_rm.end(), br.begin(), br.end());
256  }
257  std::vector<double> rnvec = norm2<valueT,NDIM>(_world, t_rm);
258  if (_world.rank() == 0)
259  {
260  double rnorm = 0.0;
261  for (unsigned int i = 0; i < rnvec.size(); i++) rnorm += rnvec[i];
262  if (_world.rank() == 0) print("residual = ", rnorm);
263  _residual = rnorm;
264  }
266 
267  for (unsigned int kp = 0; kp < _kpoints.size(); kp++)
268  {
269  KPoint kpoint = _kpoints[kp];
270 
271  vecfuncT k_phisa(awfs_old.begin() + kpoint.begin, awfs_old.begin() + kpoint.end);
272  vecfuncT k_phisb(bwfs_old.begin() + kpoint.begin, bwfs_old.begin() + kpoint.end);
273  vecfuncT k_awfs(awfs_new.begin() + kpoint.begin, awfs_new.begin() + kpoint.end);
274  vecfuncT k_bwfs(bwfs_new.begin() + kpoint.begin, bwfs_new.begin() + kpoint.end);
275 
276  // compute residuals for k-point
277  // concatentate up and down spins
278  vecfuncT k_rm = sub(_world, k_phisa, k_awfs);
279  vecfuncT k_vm(k_phisa);
280  if (_params.spinpol)
281  {
282  vecfuncT k_br = sub(_world, k_phisb, k_bwfs);
283  k_rm.insert(k_rm.end(), k_br.begin(), k_br.end());
284  k_vm.insert(k_vm.end(), k_phisb.begin(), k_phisb.end());
285  }
286 
287  // Update subspace and matrix Q
288  compress(_world, k_vm, false);
289  compress(_world, k_rm, false);
290  _world.gop.fence();
291  subspaceT k_subspace = _subspace[kp];
292  k_subspace.push_back(pairvecfuncT(k_vm,k_rm));
293 
294  int m = k_subspace.size();
295  tensorT ms(m);
296  tensorT sm(m);
297  for (int s = 0; s < m; s++)
298  {
299  const vecfuncT& k_vs = k_subspace[s].first;
300  const vecfuncT& k_rs = k_subspace[s].second;
301  for (unsigned int i = 0; i < k_vm.size(); i++)
302  {
303  ms[s] += k_vm[i].inner_local(k_rs[i]);
304  sm[s] += k_vs[i].inner_local(k_rm[i]);
305  }
306  }
307  _world.gop.sum(ms.ptr(),m);
308  _world.gop.sum(sm.ptr(),m);
309 
310  tensorT newQ(m,m);
311  if (m > 1) newQ(Slice(0,-2),Slice(0,-2)) = _Q[kp];
312  newQ(m-1,_) = ms;
313  newQ(_,m-1) = sm;
314 
315  _Q[kp] = newQ;
316  if (_world.rank() == 0) print(_Q[kp]);
317 
318  // Solve the subspace equations
319  tensorT c;
320  if (_world.rank() == 0) {
321  double rcond = 1e-12;
322  while (1) {
323  c = KAIN(_Q[kp],rcond);
324  if (abs(c[m-1]) < 3.0) {
325  break;
326  }
327  else if (rcond < 0.01) {
328  if (_world.rank() == 0)
329  print("Increasing subspace singular value threshold ", c[m-1], rcond);
330  rcond *= 100;
331  }
332  else {
333  if (_world.rank() == 0)
334  print("Forcing full step due to subspace malfunction");
335  c = 0.0;
336  c[m-1] = 1.0;
337  break;
338  }
339  }
340  }
341 
343  if (_world.rank() == 0) {
344  //print("Subspace matrix");
345  //print(Q);
346  print("Subspace solution", c);
347  }
348 
349  // Form linear combination for new solution
350  vecfuncT k_phisa_new = zero_functions_compressed<valueT,NDIM>(_world, k_phisa.size());
351  vecfuncT k_phisb_new = zero_functions_compressed<valueT,NDIM>(_world, k_phisb.size());
352  std::complex<double> one = std::complex<double>(1.0,0.0);
353  unsigned int norbitals = awfs_old.size() / _kpoints.size();
354  for (unsigned int m = 0; m < k_subspace.size(); m++)
355  {
356  const vecfuncT& k_vm = k_subspace[m].first;
357  const vecfuncT& k_rm = k_subspace[m].second;
358  // WSTHORNTON Stopped here!
359  const vecfuncT vma(k_vm.begin(),k_vm.begin() + norbitals);
360  const vecfuncT rma(k_rm.begin(),k_rm.begin() + norbitals);
361  const vecfuncT vmb(k_vm.end() - norbitals, k_vm.end());
362  const vecfuncT rmb(k_rm.end() - norbitals, k_rm.end());
363 
364  gaxpy(_world, one, k_phisa_new, c(m), vma, false);
365  gaxpy(_world, one, k_phisa_new,-c(m), rma, false);
366  gaxpy(_world, one, k_phisb_new, c(m), vmb, false);
367  gaxpy(_world, one, k_phisb_new,-c(m), rmb, false);
368  }
369  _world.gop.fence();
370 
371  if (_params.maxsub <= 1) {
372  // Clear subspace if it is not being used
373  k_subspace.clear();
374  }
375  else if (k_subspace.size() == _params.maxsub) {
376  // Truncate subspace in preparation for next iteration
377  k_subspace.erase(k_subspace.begin());
378  _Q[kp] = _Q[kp](Slice(1,-1),Slice(1,-1));
379  }
380  // Save subspace
381  _subspace[kp] = k_subspace;
382 
383  for (unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
384  wi++, fi++)
385  {
386  awfs_new[wi] = k_phisa_new[fi];
387  bwfs_new[wi] = k_phisb_new[fi];
388  }
389  }
390  }
391  //*************************************************************************
392 
393  };
394 
395  /*!
396  \ingroup periodic_solver
397 
398  \brief The SubspaceK class is a container class holding previous orbitals
399  and residuals.
400  \par
401  The Solver class uses the Krylov Accelerated Inexact Newton Solver (KAIN)
402  accelerate the convergence a given calculation. The KAIN solver needs to
403  store a subspace of previous orbitals and residuals.
404  */
405 
406  //***************************************************************************
407  template <typename T, int NDIM>
408  class Subspace
409  {
410  // Typedef's
411  typedef std::complex<T> valueT;
413  typedef std::vector<functionT> vecfuncT;
414  typedef std::pair<vecfuncT,vecfuncT> pairvecfuncT;
415  typedef std::vector<pairvecfuncT> subspaceT;
417 
418  //*************************************************************************
420  //*************************************************************************
421 
422  //*************************************************************************
424  //*************************************************************************
425 
426  //*************************************************************************
428  //*************************************************************************
429 
430  //*************************************************************************
431  std::vector<KPoint> _kpoints;
432  //*************************************************************************
433 
434  //*************************************************************************
436  //*************************************************************************
437 
438  public:
439 
440  //*************************************************************************
441  Subspace(World& world, const ElectronicStructureParams& params)
442  : _world(world), _params(params)
443  {
444  }
445  //*************************************************************************
446 
447  //*************************************************************************
448  void update_subspace(vecfuncT& awfs_new,
449  vecfuncT& bwfs_new,
450  const vecfuncT& awfs_old,
451  const vecfuncT& bwfs_old,
452  const vecfuncT& rm)
453  {
454  // concatentate up and down spins
455  vecfuncT vm = awfs_old;
456  if (_params.spinpol)
457  {
458  vm.insert(vm.end(), bwfs_old.begin(), bwfs_old.end());
459  }
460 
461  // Update subspace and matrix Q
462  compress(_world, vm, false);
463  compress(_world, rm, false);
464  _world.gop.fence();
465  _subspace.push_back(pairvecfuncT(vm,rm));
466 
467  int m = _subspace.size();
468  tensorT ms(m);
469  tensorT sm(m);
470  for (int s=0; s<m; s++)
471  {
472  const vecfuncT& vs = _subspace[s].first;
473  const vecfuncT& rs = _subspace[s].second;
474  for (unsigned int i=0; i<vm.size(); i++)
475  {
476  ms[s] += vm[i].inner_local(rs[i]);
477  sm[s] += vs[i].inner_local(rm[i]);
478  }
479  }
480  _world.gop.sum(ms.ptr(),m);
481  _world.gop.sum(sm.ptr(),m);
482 
483  tensorT newQ(m,m);
484  if (m > 1) newQ(Slice(0,-2),Slice(0,-2)) = _Q;
485  newQ(m-1,_) = ms;
486  newQ(_,m-1) = sm;
487 
488  _Q = newQ;
489  if (_world.rank() == 0) print(_Q);
490 
491  // Solve the subspace equations
492  tensorT c;
493  if (_world.rank() == 0) {
494  double rcond = 1e-12;
495  while (1) {
496  c = KAIN(_Q,rcond);
497  if (abs(c[m-1]) < 3.0) {
498  break;
499  }
500  else if (rcond < 0.01) {
501  if (_world.rank() == 0)
502  print("Increasing subspace singular value threshold ", c[m-1], rcond);
503  rcond *= 100;
504  }
505  else {
506  if (_world.rank() == 0)
507  print("Forcing full step due to subspace malfunction");
508  c = 0.0;
509  c[m-1] = 1.0;
510  break;
511  }
512  }
513  }
514 
516  if (_world.rank() == 0) {
517  //print("Subspace matrix");
518  //print(Q);
519  print("Subspace solution", c);
520  }
521 
522  // Form linear combination for new solution
523  vecfuncT phisa_new = zero_functions_compressed<valueT,NDIM>(_world, awfs_old.size());
524  vecfuncT phisb_new = zero_functions_compressed<valueT,NDIM>(_world, bwfs_old.size());
525  std::complex<double> one = std::complex<double>(1.0,0.0);
526  for (unsigned int m=0; m<_subspace.size(); m++) {
527  const vecfuncT& vm = _subspace[m].first;
528  const vecfuncT& rm = _subspace[m].second;
529  const vecfuncT vma(vm.begin(),vm.begin()+awfs_old.size());
530  const vecfuncT rma(rm.begin(),rm.begin()+awfs_old.size());
531  const vecfuncT vmb(vm.end()-bwfs_old.size(), vm.end());
532  const vecfuncT rmb(rm.end()-bwfs_old.size(), rm.end());
533 
534  gaxpy(_world, one, phisa_new, c(m), vma, false);
535  gaxpy(_world, one, phisa_new,-c(m), rma, false);
536  gaxpy(_world, one, phisb_new, c(m), vmb, false);
537  gaxpy(_world, one, phisb_new,-c(m), rmb, false);
538  }
539  _world.gop.fence();
540 
541  if (_params.maxsub <= 1) {
542  // Clear subspace if it is not being used
543  _subspace.clear();
544  }
545  else if (_subspace.size() == _params.maxsub) {
546  // Truncate subspace in preparation for next iteration
547  _subspace.erase(_subspace.begin());
548  _Q = _Q(Slice(1,-1),Slice(1,-1));
549  }
550  awfs_new = phisa_new;
551  bwfs_new = phisb_new;
552  }
553  //*************************************************************************
554 
555  //*************************************************************************
556  void reproject()
557  {
558  // //if (_world.rank() == 0)
559  // //printf("\n\nreprojecting subspace to wavelet order: %d and thresh: %.5e\n\n",
560  // //FunctionDefaults<3>::get_k(), FunctionDefaults<3>::get_thresh());
561  //
562  // unsigned int m = _subspace.size();
563  // for (unsigned int s = 0; s < m; s++)
564  // {
565  // vecfuncT& vs = _subspace[s].first;
566  // vecfuncT& rs = _subspace[s].second;
567  // reconstruct(_world, vs);
568  // reconstruct(_world, rs);
569  // unsigned int vm = vs.size();
570  // for (unsigned int i = 0; i < vm; i++)
571  // {
572  // vs[i] = madness::project(vs[i], FunctionDefaults<3>::get_k(),
573  // FunctionDefaults<3>::get_thresh(), false);
574  // rs[i] = madness::project(rs[i], FunctionDefaults<3>::get_k(),
575  // FunctionDefaults<3>::get_thresh(), false);
576  // }
577  // _world.gop.fence();
578  // truncate(_world, vs);
579  // truncate(_world, rs);
580  // normalize(_world, vs);
581  // }
582  // _world.gop.fence();
583 
584  }
585  //*************************************************************************
586 
587  };
588 
589 // template <typename T, int NDIM>
590 // struct lbcost {
591 // double leaf_value;
592 // double parent_value;
593 // lbcost(double leaf_value=1.0, double parent_value=0.0) : leaf_value(leaf_value), parent_value(parent_value) {}
594 // double operator()(const Key<NDIM>& key, const FunctionNode<T,NDIM>& node) const {
595 // if (key.level() <= 1) {
596 // return 100.0*(leaf_value+parent_value);
597 // }
598 // else if (node.is_leaf()) {
599 // return leaf_value;
600 // }
601 // else {
602 // return parent_value;
603 // }
604 // }
605 // };
606 
607  //***************************************************************************
608 
609  /*! \ingroup periodic_solver
610  \brief The main class of the periodic DFT solver
611  \f[
612  z = frac{x}{1 - y^2}
613  \f]
614  */
615 
616  //***************************************************************************
617  template <typename T, int NDIM>
618  class Solver
619  {
620  // Typedef's
621  typedef std::complex<T> valueT;
625  typedef std::vector<functionT> vecfuncT;
629  typedef std::shared_ptr<operatorT> poperatorT;
633  typedef std::pair<vecfuncT,vecfuncT> pairvecfuncT;
634  typedef std::vector<pairvecfuncT> subspaceT;
635  typedef std::vector<tensorT> vectensorT;
636  typedef std::vector<subspaceT> vecsubspaceT;
637 
638  //*************************************************************************
640  //*************************************************************************
641 
642  //*************************************************************************
643  // This variable could either be a nuclear potiential or a nuclear charge
644  // density depending on the "ispotential" boolean variable in the
645  // ElectronicStructureParams class.
647  //*************************************************************************
648 
649  //*************************************************************************
651  //*************************************************************************
652 
653  //*************************************************************************
655  //*************************************************************************
656 
657  //*************************************************************************
658  std::vector<T> _eigsa;
659  //*************************************************************************
660 
661  //*************************************************************************
662  std::vector<T> _eigsb;
663  //*************************************************************************
664 
665  //*************************************************************************
667  //*************************************************************************
668 
669  //*************************************************************************
670  std::vector<KPoint> _kpoints;
671  //*************************************************************************
672 
673  //*************************************************************************
674  std::vector<double> _occsa;
675  //*************************************************************************
676 
677  //*************************************************************************
678  std::vector<double> _occsb;
679  //*************************************************************************
680 
681  //*************************************************************************
683  //*************************************************************************
684 
685  //*************************************************************************
687  //*************************************************************************
688 
689  //*************************************************************************
691  //*************************************************************************
692 
693  //*************************************************************************
695  //*************************************************************************
696 
697  //*************************************************************************
699  //*************************************************************************
700 
701 // //*************************************************************************
702 // vecsubspaceT _subspace;
703 // //*************************************************************************
704 //
705 // //*************************************************************************
706 // vectensorT _Q;
707 // //*************************************************************************
708 
709  //*************************************************************************
711  //*************************************************************************
712 
713  //*************************************************************************
715  //*************************************************************************
716 
717  //*************************************************************************
718  double _residual;
719  //*************************************************************************
720 
721  //*************************************************************************
723  //*************************************************************************
724 
725  //*************************************************************************
726  double _maxthresh;
727  //*************************************************************************
728 
729  //*************************************************************************
730  std::ofstream _outputF;
731  //*************************************************************************
732 
733  //*************************************************************************
734  std::ofstream _matF;
735  //*************************************************************************
736 
737  //*************************************************************************
738  std::ofstream _eigF;
739  //*************************************************************************
740 
741  //*************************************************************************
742  std::ofstream _kmeshF;
743  //*************************************************************************
744 
745  //*************************************************************************
746  int _it;
747  //*************************************************************************
748 
749  //*************************************************************************
751  //*************************************************************************
752 
753  //*************************************************************************
754  int _nao;
755  //*************************************************************************
756 
757  public:
758 
759  //*************************************************************************
760  double ttt, sss;
761  void START_TIMER(World& world) {
762  world.gop.fence(); ttt=wall_time(); sss=cpu_time();
763  }
764 
765  void END_TIMER(World& world, const char* msg) {
766  ttt=wall_time()-ttt; sss=cpu_time()-sss; if (world.rank()==0) printf("timer: %20.20s %8.2fs %8.2fs\n", msg, sss, ttt);
767  }
768  //*************************************************************************
769 
770  //*************************************************************************
771  Solver(World& world, const std::string& filename) : _world(world)
772  {
773  init(filename);
775  _residual = 1e5;
777 
778  if (_params.restart==0) initial_guess();
779 // for (unsigned int kp = 0; kp < _kpoints.size(); kp++)
780 // {
781 // _Q.push_back(tensorT(1,1));
782 // _subspace.push_back(subspaceT());
783 // }
784  _subspace = new Subspace<T,NDIM>(world, _params);
785  }
786  //*************************************************************************
787 
788  //*************************************************************************
789  void init(const std::string& filename)
790  {
791 
792  // params
793  if (_world.rank() == 0)
794  {
796  //_params.fractional = false;
797  }
798  // Send params
800  if (_params.centered)
802  else
804 // FunctionDefaults<3>::set_cubic_cell(0,_params.L);
807 
808  // mentity and aobasis
809  if (_world.rank() == 0)
810  {
814  }
815  // Send mentity and aobasis
818  // set number of electrons to the total nuclear charge of the mentity
820  // total number of bands include empty
823  if ((_params.nelec % 2) == 1) _params.nelec++;
824 
825  // Open output files
826  _outputF.open("output.txt");
827  _matF.open("matrix.txt");
828  _eigF.open("eigs.txt");
829 
830  // kmesh
831  if (_params.restart == 0)
832  {
833  if (_params.periodic) // PERIODIC
834  {
835  // GAMMA POINT
836  if ((_params.ngridk0 == 1) && (_params.ngridk1 == 1) && (_params.ngridk2 == 1))
837  {
838  _kpoints.push_back(KPoint(coordT(0.0), 1.0));
839  }
840  if ((_params.ngridk0 == 0) && (_params.ngridk1 == 0) && (_params.ngridk2 == 0))
841  {
842  double TWO_PI = 2.0 * madness::constants::pi;
843  double t1 = TWO_PI/_params.L;
844  coordT c1 {0.0,0.0,0.0};
845  coordT c2 {0.5*t1,0.5*t1,0.5*t1};
846  _kpoints.push_back(KPoint(c1, 0.5));
847  _kpoints.push_back(KPoint(c2, 0.5));
848  }
849  else // NORMAL BANDSTRUCTURE
850  {
854  _params.L);
855  }
856  }
857  else // NOT-PERIODIC
858  {
859  _kpoints.push_back(KPoint(coordT(0.0), 1.0));
860  }
861  if (_world.rank() == 0)
862  {
863  _kmeshF.open("kpoints.txt");
864  _kmeshF << "kpts: " << _kpoints.size() << endl;
865  _kmeshF << "ik" << setw(10) << "kpt" << setw(30) << "weight" << endl;
866  _kmeshF << "--" << setw(10) << "---" << setw(30) << "------" << endl;
867 
868  //_kmeshF << setfill('-') << setw(55) << "-" << endl;
869  //_kmeshF << setfill(' ') << endl;
870  _kmeshF << endl;
871  for (unsigned int i = 0; i < _kpoints.size(); i++)
872  {
873  KPoint kpoint = _kpoints[i];
874  _kmeshF << i << setw(10) << kpoint.k[0];
875  _kmeshF << setw(10) << kpoint.k[1];
876  _kmeshF << setw(10) << kpoint.k[2];
877  _kmeshF << setw(10) << kpoint.weight << endl;
878  }
879  _kmeshF.close();
880  }
881  }
882  else
883  {
884  load_orbitals();
885  }
886  }
887  //*************************************************************************
888 
889  //*************************************************************************
890  std::vector<KPoint> genkmesh(unsigned int ngridk0, unsigned ngridk1, unsigned int ngridk2,
891  double koffset0, double koffset1, double koffset2, double R)
892  {
893  std::vector<KPoint> kmesh;
894  double step0 = 1.0/ngridk0;
895  double step1 = 1.0/ngridk1;
896  double step2 = 1.0/ngridk2;
897  double weight = 1.0/(ngridk0*ngridk1*ngridk2);
898  double TWO_PI = 2.0 * madness::constants::pi;
899  for (unsigned int i = 0; i < ngridk0; i++)
900  {
901  for (unsigned int j = 0; j < ngridk1; j++)
902  {
903  for (unsigned int k = 0; k < ngridk2; k++)
904  {
905  //double k0 = (i*step0 - step0/2) * TWO_PI/R;
906  //double k1 = (j*step1 - step1/2) * TWO_PI/R;
907  //double k2 = (k*step2 - step2/2) * TWO_PI/R;
908  double k0 = ((i*step0)+koffset0) * TWO_PI/R;
909  double k1 = ((j*step1)+koffset1) * TWO_PI/R;
910  double k2 = ((k*step2)+koffset2) * TWO_PI/R;
911  KPoint kpoint(k0, k1, k2, weight);
912  kmesh.push_back(kpoint);
913  }
914  }
915  }
916  return kmesh;
917  }
918  //*************************************************************************
919 
920  //*************************************************************************
922  {
924  ar & _params.spinpol;
925  ar & (unsigned int)(_kpoints.size());
926  for (unsigned int i = 0; i < _kpoints.size(); i++) ar & _kpoints[i];
927 // ar & (unsigned int)(_occsa.size());
928 // for (unsigned int i = 0; i < _occsa.size(); i++) ar & _occsa[i];
929  ar & (unsigned int)(_phisa.size());
930  for (unsigned int i = 0; i < _phisa.size(); i++) ar & _phisa[i];
931  if (_params.spinpol)
932  {
933 // ar & (unsigned int)(_occsb.size());
934 // for (unsigned int i = 0; i < _occsb.size(); i++) ar & _occsb[i];
935  ar & (unsigned int)(_phisb.size());
936  for (unsigned int i = 0; i < _phisb.size(); i++) ar & _phisb[i];
937  }
938  }
939  //*************************************************************************
940 
941  //*************************************************************************
943  {
944  const double thresh = FunctionDefaults<3>::get_thresh();
945  const int k = FunctionDefaults<3>::get_k();
946 
947  archive::ParallelInputArchive ar(_world, "orbitals");
948 
949  // spin-polarized
950  bool spinrest;
951  ar & spinrest;
952  // kpoints
953  unsigned int nkpts;
954  ar & nkpts;
955  _kpoints.clear();
956  for (unsigned int i = 0; i < nkpts; i++)
957  {
958  KPoint tkpt;
959  ar & tkpt;
960  _kpoints.push_back(tkpt);
961  }
962  // occs
963 // unsigned int noccs;
964 // ar & noccs;
965 // _occs.clear();
966 // for (unsigned int i = 0; i < noccs; i++)
967 // {
968 // double tocc;
969 // ar & tocc;
970 // _occs.push_back(tocc);
971 // }
972  // orbitals
973  unsigned int norbs;
974  ar & norbs;
975  _phisa.clear();
976  _eigsa.clear();
977  for (unsigned int i = 0; i < norbs; i++)
978  {
979  functionT tfunc;
980  ar & tfunc;
981  _phisa.push_back(tfunc);
982  _eigsa.push_back(-0.1);
983  }
984  // check for k mismatch
985  if (_phisa[0].k() != k)
986  {
988  for (unsigned int i = 0; i < _phisa.size(); i++)
989  _phisa[i] = madness::project(_phisa[i], k, thresh, false);
990  _world.gop.fence();
991  }
992  // orthonormalize
993  for (unsigned int i = 0; i < _kpoints.size(); i++)
995 
996  if (_params.spinpol)
997  {
998  _phisb.clear();
999  _eigsb.clear();
1000  for (unsigned int i = 0; i < norbs; i++)
1001  {
1002  functionT tfunc;
1003  ar & tfunc;
1004  _phisb.push_back(tfunc);
1005  _eigsb.push_back(-0.1);
1006  }
1007  // check for k mismatch
1008  if (_phisb[0].k() != k)
1009  {
1011  for (unsigned int i = 0; i < _phisb.size(); i++)
1012  _phisb[i] = madness::project(_phisb[i], k, thresh, false);
1013  _world.gop.fence();
1014  }
1015  // orthonormalize
1016  for (unsigned int i = 0; i < _kpoints.size(); i++)
1018  }
1019  else
1020  {
1021  for (unsigned int i = 0; i < norbs; i++)
1022  {
1023  _phisb.push_back(_phisa[i]);
1024  _eigsb.push_back(_eigsa[i]);
1025  }
1026  }
1027  // create vector for occupation numbers
1028  _occsa = std::vector<double>(norbs, 0.0);
1029  _occsb = std::vector<double>(norbs, 0.0);
1030  }
1031  //*************************************************************************
1032 
1033  //*************************************************************************
1035  {
1036  _vnucrhon = rfactoryT(_world).functor(rfunctorT(new MolecularPotentialFunctor(_mentity))).thresh(_params.thresh * 0.1).truncate_on_project();
1037  _vnuc = copy(_vnucrhon);
1038  }
1039  //*************************************************************************
1040 
1041  //*************************************************************************
1043  {
1044  std::vector<coordT> specialpts;
1045  for (int i = 0; i < _mentity.natom(); i++)
1046  {
1047  coordT pt(0.0);
1048  Atom atom = _mentity.get_atom(i);
1049  pt[0] = atom.x; pt[1] = atom.y; pt[2] = atom.z;
1050  specialpts.push_back(pt);
1051  if (_world.rank() == 0) print("Special point: ", pt);
1052  }
1053  double now = wall_time();
1054  std::cout << "THE FUNCTOR IS " << MolecularNuclearChargeDensityFunctor(_mentity, _params.L, _params.periodic, specialpts)(specialpts[0]) << endl;
1055  // WSTHORNTON
1056  //MADNESS_ASSERT(false);
1057 
1058 
1059  _vnucrhon = rfactoryT(_world).functor(
1061  thresh(_params.thresh).initial_level(6).truncate_on_project();
1062 
1063  if (_world.rank() == 0) printf("%f\n", wall_time() - now);
1064  if (_world.rank() == 0) print("calculating trace of rhon ..\n\n");
1065  double rtrace = _vnucrhon.trace();
1066  if (_world.rank() == 0) print("rhon trace = ", rtrace);
1067  now = wall_time();
1068  _vnucrhon.truncate();
1069  _vnuc = apply(*_cop, _vnucrhon);
1070  if (_world.rank() == 0) printf("%f\n", wall_time() - now);
1071  if (_world.rank() == 0) print("Done creating nuclear potential ..\n");
1072 
1073 
1074  // WSTHORNTON
1075  rfunctionT vnuc2 = rfactoryT(_world).functor(rfunctorT(new
1077  thresh(_params.thresh * 0.1).truncate_on_project();
1078  rfunctionT vnucdiff = _vnuc-vnuc2;
1079  double t1 = vnucdiff.trace();
1080  if (_world.rank() == 0) printf("Difference in nuclear potential: %15.8f\n\n", t1);
1081  for (int i = 0; i < 101; i++)
1082  {
1083  double dx = _params.L/100;
1084  double pt2 = -_params.L/2 + dx*i;
1085  coordT cpt(pt2);
1086  double val = vnucdiff(cpt);
1087  if (_world.rank() == 0) printf("%10.5f %15.8f\n",pt2,val);
1088  }
1089  //MADNESS_ASSERT(false);
1090  }
1091  //*************************************************************************
1092 
1093  //*************************************************************************
1095  {
1096  // Make coulomb operator
1098 
1100  if (_world.rank() == 0)
1101  {
1102  printf("Cell parameters\n");
1103  printf("------------------------------\n");
1104  print("cell(x) is ",csize(0,0), csize(0,1));
1105  print("cell(y) is ",csize(1,0), csize(1,1));
1106  print("cell(z) is ",csize(2,0), csize(2,1));
1107  printf("\n");
1108  }
1109  if (_params.ispotential) // potential
1110  {
1112  }
1113  else // charge density
1114  {
1116  }
1117  }
1118  //*************************************************************************
1119 
1120  //*************************************************************************
1121  struct GuessDensity : public FunctionFunctorInterface<double,3> {
1124  double R;
1125  const bool periodic;
1126 
1127  double operator()(const coordT& x) const
1128  {
1129  double value = 0.0;
1130  if (periodic)
1131  {
1132  for (int xr = -2; xr <= 2; xr += 1)
1133  {
1134  for (int yr = -2; yr <= 2; yr += 1)
1135  {
1136  for (int zr = -2; zr <= 2; zr += 1)
1137  {
1139  x[0]+xr*R, x[1]+yr*R, x[2]+zr*R);
1140  }
1141  }
1142  }
1143  }
1144  else
1145  {
1146  value = aobasis.eval_guess_density(mentity, x[0], x[1], x[2]);
1147  }
1148  return value;
1149  }
1150 
1152  const double& R, const bool& periodic)
1154  };
1155  //*************************************************************************
1156 
1157  //*************************************************************************
1158  template <typename Q>
1160  {
1161  int npts = 101;
1162  double eps = 1e-8;
1163  coordT r(0.0);
1164  double delta = _params.L/(npts-1);
1165  double begin = -_params.L/2;
1166  double end = _params.L/2;
1167  double tol = 1e-6;
1168  // x-axis
1169  printf("periodicity along x-axis: \n");
1170  printf("-------------------------------------------------------------\n\n");
1171  for (int i = 0; i < npts; i++)
1172  {
1173  printf("\n-------------------- edge --------------------\n");
1174  for (int j = 0; j < npts; j++)
1175  {
1176  coordT r1, r2, dr1, dr2;
1177  r1[0] = begin+eps; r1[1] = (i*delta)+begin; r1[2] = (j*delta)+begin;
1178  r2[0] = end-eps; r2[1] = (i*delta)+begin; r2[2] = (j*delta)+begin;
1179  dr1[0] = begin+2*eps; dr1[1] = (i*delta)+begin; dr1[2] = (j*delta)+begin;
1180  dr2[0] = end-2*eps; dr2[1] = (i*delta)+begin; dr2[2] = (j*delta)+begin;
1181  double val = std::abs(f(r1)-f(r2));
1182  std::string success = val < tol ? "PASS!" : "FAIL!";
1183  printf("%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1184  r1[0],r1[1],r1[2],r2[0],r2[1],r2[2],val,success.c_str());
1185  }
1186  }
1187  // y-axis
1188  printf("\nperiodicity along y-axis: \n");
1189  printf("-------------------------------------------------------------\n\n");
1190  for (int i = 0; i < npts; i++)
1191  {
1192  printf("\n-------------------- edge --------------------\n");
1193  for (int j = 0; j < npts; j++)
1194  {
1195  coordT r1, r2, dr1, dr2;
1196  r1[0] = (i*delta)+begin; r1[1] = begin+eps; r1[2] = (j*delta)+begin;
1197  r2[0] = (i*delta)+begin; r2[1] = end-eps; r2[2] = (j*delta)+begin;
1198  dr1[0] = (i*delta)+begin; dr1[1] = begin+2*eps; dr1[2] = (j*delta)+begin;
1199  dr2[0] = (i*delta)+begin; dr2[1] = end-2*eps; dr2[2] = (j*delta)+begin;
1200  double val = std::abs(f(r1)-f(r2));
1201  std::string success = val < tol ? "PASS!" : "FAIL!";
1202  printf("%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1203  r1[0],r1[1],r1[2],r2[0],r2[1],r2[2],val,success.c_str());
1204  }
1205  }
1206  // z-axis
1207  printf("\nperiodicity along z-axis: \n");
1208  printf("-------------------------------------------------------------\n\n");
1209  for (int i = 0; i < npts; i++)
1210  {
1211  printf("\n-------------------- edge --------------------\n");
1212  for (int j = 0; j < npts; j++)
1213  {
1214  coordT r1, r2, dr1, dr2;
1215  r1[0] = (i*delta)+begin; r1[1] = (j*delta)+begin; r1[2] = begin+eps;
1216  r2[0] = (i*delta)+begin; r2[1] = (j*delta)+begin; r2[2] = end-eps;
1217  dr1[0] = (i*delta)+begin; dr1[1] = (j*delta)+begin; dr1[2] = begin+2*eps;
1218  dr2[0] = (i*delta)+begin; dr2[1] = (j*delta)+begin; dr2[2] = end-2*eps;
1219  double val = std::abs(f(r1)-f(r2));
1220  std::string success = val < tol ? "PASS!" : "FAIL!";
1221  printf("%10.6f%10.6f%10.6f %10.6f%10.6f%10.6f %10.5e %s\n",
1222  r1[0],r1[1],r1[2],r2[0],r2[1],r2[2],val,success.c_str());
1223  }
1224  }
1225  }
1226  //*************************************************************************
1227 
1228  //*************************************************************************
1229  rfunctionT
1231  const rfunctionT& arho,
1232  const rfunctionT& brho,
1233  const rfunctionT& adelrhosq,
1234  const rfunctionT& bdelrhosq)
1235  {
1236  // MADNESS_ASSERT(!_params.spinpol);
1237  rfunctionT vlda = copy(arho);
1238 // vlda.unaryop(&::libxc_ldaop<double>);
1239  vlda.unaryop(&::ldaop<double>);
1240  return vlda;
1241  }
1242 
1245 
1246 // Level initial_level = 3;
1247 // for (int i=0; i < _aobasis.nbf(_mentity); i++) {
1248 // functorT aofunc(new AtomicBasisFunctor(_aobasis.get_atomic_basis_function(_mentity,i),
1249 // _params.L, _params.periodic, kpt));
1250 // ao[i] = factoryT(world).functor(aofunc).initial_level(initial_level).truncate_on_project().nofence();
1251 // }
1252 // world.gop.fence();
1253 
1254 // WSTAtomicBasisFunctor aofunc(_aobasis.get_atomic_basis_function(_mentity,4),
1255 // _params.L, kpt.k[0], kpt.k[1], kpt.k[2]);
1256 // coord_3d p(5e-2);
1257 // double_complex val = aofunc(p);
1258 // print("AO_PROJECT: val is ", val);
1259 // MADNESS_ASSERT(false);
1260 
1261  for (int i=0; i < _aobasis.nbf(_mentity); i++) {
1263  _params.L, kpt.k[0], kpt.k[1], kpt.k[2]));
1264  ao[i] = factoryT(world).functor(aofunc).truncate_on_project().truncate_mode(0);
1265  }
1266  world.gop.fence();
1267 
1268 // std::vector<double> norms;
1269 // norms = norm2s(world, ao);
1270 // for (int i=0; i<_aobasis.nbf(_mentity); i++) {
1271 // if (world.rank() == 0 && fabs(norms[i]-1.0)>1e-3) print(i," bad ao norm?", norms[i]);
1272 // norms[i] = 1.0/norms[i];
1273 // }
1274 //
1275 // scale(world, ao, norms);
1276 // norms = norm2s(world, ao);
1277 // for (int i=0; i<_aobasis.nbf(_mentity); i++) {
1278 // if (world.rank() == 0 && fabs(norms[i]-1.0)>1e-3) print(i," bad ao norm?", norms[i]);
1279 // norms[i] = 1.0/norms[i];
1280 // }
1281 // scale(world, ao, norms);
1282 
1283  return ao;
1284  }
1285  //*************************************************************************
1286 
1287  //*************************************************************************
1288  // Constructor
1289  Solver(World& world,
1290  rfunctionT vnucrhon,
1291  vecfuncT phisa,
1292  vecfuncT phisb,
1293  std::vector<T> eigsa,
1294  std::vector<T> eigsb,
1296  MolecularEntity mentity)
1297  : _world(world), _vnucrhon(vnucrhon), _phisa(phisa), _phisb(phisb),
1298  _eigsa(eigsa), _eigsb(eigsb), _params(params), _mentity(mentity)
1299  {
1300  _residual = 1e5;
1301  _cop = CoulombOperatorPtr(const_cast<World&>(world), params.lo, params.thresh * 0.1);
1302 
1303  if (params.ispotential)
1304  {
1305  _vnuc = copy(_vnucrhon);
1306  }
1307  else
1308  {
1309  _vnuc = apply(*_cop, _vnucrhon);
1310  }
1311  }
1312  //*************************************************************************
1313 
1314  /// Initializes alpha and beta mos, occupation numbers, eigenvalues
1315  //*************************************************************************
1317  {
1318  // Get initial guess for the electronic density
1319  if (_world.rank() == 0) print("Guessing rho ...\n\n");
1320 
1321  rfunctionT rho = rfactoryT(_world);
1322  if (_params.restart == 0)
1323  {
1324  rho = rfactoryT(_world).functor(rfunctorT(new GuessDensity(_mentity,
1325  _aobasis, _params.L, _params.periodic))).initial_level(3);
1326  }
1327  else
1328  {
1329  MADNESS_EXCEPTION("restart not working right now!",0);
1330  }
1331  // This is a cheat
1332  rho.scale(_params.ncharge/rho.trace());
1333 
1334  char fname3[25];
1335  coord_3d rp1(-_params.L/2);
1336  coord_3d rp2(_params.L/2);
1337  plot_line(fname3,101,rp1,rp2,rho);
1338 
1339  // load balance
1340  //if(_world.size() > 1) {
1341  // START_TIMER(_world);
1342  // LoadBalanceDeux<3> lb(_world);
1343  // lb.add_tree(_vnuc, lbcost<double,3>(1.0, 0.0), false);
1344  // lb.add_tree(rho, lbcost<double,3>(1.0, 1.0), true);
1345 
1346  // FunctionDefaults<3>::redistribute(_world, lb.load_balance(6.0));
1347  //}
1348 
1349  if (_params.restart != 1)
1350  {
1351  // build effective potential
1352  rfunctionT vlocal;
1353  // Is this a many-body system?
1354  if (_params.ncharge > 1.0)
1355  {
1356  if (_world.rank() == 0) print("Creating Coulomb op ...\n\n");
1358  // Is this system periodic?
1359  if (_params.periodic)
1360  {
1362  //op = PeriodicCoulombOpPtr<double, 3> (_world, _params.waveorder,
1363  // _params.lo, _params.thresh * 0.1, cellsize);
1365  }
1366  else
1367  {
1369  }
1370  if (_world.rank() == 0) print("Building effective potential ...\n\n");
1371  rfunctionT vc = apply(*op, rho);
1372  vlocal = _vnuc + vc; //.scale(1.0-1.0/nel); // Reduce coulomb to increase binding
1373  rho.scale(0.5);
1374  // WSTHORNTON
1375  _rhoa = rho;
1376  _rhob = rho;
1377  // Do the LDA
1378  rfunctionT vlda = make_lda_potential(_world, rho, rho, rfunctionT(), rfunctionT());
1379  vlocal = vlocal + vlda;
1380  delete op;
1381  }
1382  else
1383  {
1384  vlocal = _vnuc;
1385  }
1386 
1387  // Clear these functions
1388 // rho.clear();
1389  vlocal.reconstruct();
1390 
1391  // Get size information from k-points and ao_basis so that we can correctly size
1392  // the _orbitals data structure and the eigs tensor
1393  // number of orbitals in the basis set
1394  int _nao = _aobasis.nbf(_mentity);
1395 
1396  // number of kpoints
1397  int nkpts = _kpoints.size();
1398  // total number of orbitals to be processed (no symmetry)
1399  int norbs = _params.nbands * nkpts;
1400  // Check to see if the basis set can accomodate the number of bands
1401  if (_params.nbands > _nao)
1402  MADNESS_EXCEPTION("Error: basis not large enough to accomodate number of bands", _nao);
1403  // set the number of orbitals
1404  _eigsa = std::vector<double>(norbs, 0.0);
1405  _eigsb = std::vector<double>(norbs, 0.0);
1406  _occsa = std::vector<double>(norbs, 0.0);
1407  _occsb = std::vector<double>(norbs, 0.0);
1408  if (_world.rank() == 0) print("Building kinetic energy matrix ...\n\n");
1409  // Need to do kinetic piece for every k-point
1410  for (int ki = 0; ki < nkpts; ki++)
1411  {
1412  // These are our initial basis functions
1413  if (_world.rank() == 0) print("Projecting atomic orbitals ...\n\n");
1416  END_TIMER(_world, "projecting atomic orbital basis");
1417 // for (unsigned int iao = 0; iao < _nao; iao)
1418 // {
1419 // _aobasisf.push_back(ao[iao]);
1420 // }
1421 
1422  for (unsigned int ai = 0; ai < ao.size(); ai++)
1423  {
1424  std::vector<long> npt(3,101);
1425  char fnamedx[50];
1426  sprintf(fnamedx, "aofunc_k_%2.2d__%2.2d__.dx",ki,ai);
1427  std::vector<long> npt2(3,101);
1428  plotdx(ao[ai], fnamedx, FunctionDefaults<3>::get_cell(), npt2);
1429  }
1430 
1431  coord_3d p1(-_params.L/2);
1432  coord_3d p2(_params.L/2);
1433  for (int i = 0; i < ao.size(); i++)
1434  {
1435  char fname[25];
1436  sprintf(fname, "ao_line_%d.txt",i);
1437  plot_line(fname,101,p1,p2,ao[i]);
1438  }
1439  char fname2[25];
1440  sprintf(fname2, "vnuc_line.txt");
1441  plot_line(fname2,101,p1,p2,_vnuc);
1442  char fname3[25];
1443  sprintf(fname3, "vlocal_line.txt");
1444  plot_line(fname3,101,p1,p2,vlocal);
1445 
1446 
1447  // load balancing
1448  //if(_world.size() > 1)
1449  //{
1450  // LoadBalanceDeux<3> lb(_world);
1451  // lb.add_tree(_vnuc, lbcost<double,3>(1.0, 1.0), false);
1452  // for(unsigned int i = 0; i < ao.size(); i++)
1453  // {
1454  // lb.add_tree(ao[i], lbcost<valueT,3>(1.0, 1.0), false);
1455  // }
1456 
1457  // FunctionDefaults<3>::redistribute(_world, lb.load_balance(6.0));
1458  //}
1459 
1460  // Get k-point from list
1461  KPoint& kpoint = _kpoints[ki];
1462 
1463 
1464 
1465 
1466  // WSTHORNTON
1467 // ctensorT S = matrix_inner(_world,ao,ao,true);
1468 // print(S);
1469 // ctensorT U = csqrt(S);
1470 // ao = transform(_world, ao, U, _params.thresh, true);
1471 
1472 
1473 
1474 
1475  // Build kinetic matrx
1476  //ctensorT kinetic = ::kinetic_energy_matrix_slow<T,NDIM>(_world, ao, _params.periodic, kpt);
1477  ctensorT kinetic = ::kinetic_energy_matrix_slow<T,NDIM>(_world, ao, _params.periodic, kpoint);
1478  // Build the overlap matrix
1479  if (_world.rank() == 0) print("Building overlap matrix ...\n\n");
1480  ctensorT overlap = matrix_inner(_world, ao, ao, true);
1481  // Build the potential matrix
1482  reconstruct(_world, ao);
1483  if (_world.rank() == 0) print("Building potential energy matrix ...\n\n");
1484  //cvecfuncT vpsi = mul_sparse(_world, vlocal, ao, _params.thresh);
1485  cvecfuncT vpsi;
1486  for (int i = 0; i < ao.size(); i++)
1487  vpsi.push_back(vlocal*ao[i]);
1488  _world.gop.fence();
1489  compress(_world, vpsi);
1490  truncate(_world, vpsi);
1491  compress(_world, ao);
1492  // Build the potential matrix
1493  ctensorT potential = matrix_inner(_world, vpsi, ao, true);
1494  _world.gop.fence();
1495  // free memory
1496  vpsi.clear();
1497  _world.gop.fence();
1498 
1499  // Construct and diagonlize Fock matrix
1500  ctensorT fock = potential + kinetic;
1501  ctensorT fockzero = fock-conj_transpose(fock);
1502  if (_world.rank() == 0)
1503  {
1504  print("Kinetic:");
1505  print(kinetic);
1506  print("Potential:");
1507  print(potential);
1508  print("Fock: (pre-symmetrized)");
1509  print(fock);
1510  print("FockZero: (should be zero)");
1511  print(fockzero);
1512  print("Overlap:");
1513  print(overlap);
1514  }
1515  fock = 0.5 * (fock + conj_transpose(fock));
1516  for (unsigned int i = 0; i < fock.dim(0); i++)
1517  {
1518  fock(i,i) += i*_params.thresh*0.1;
1519  }
1520 
1521  ctensorT c; rtensorT e;
1522  sygv(fock, overlap, 1, c, e);
1523 
1524  // diagonlize kinetic
1525  ctensorT ck; rtensorT ek;
1526  sygv(kinetic, overlap, 1, ck, ek);
1527  // diagonalize potential
1528  ctensorT cp; rtensorT ep;
1529  sygv(potential, overlap, 1, cp, ep);
1530  // diagonalize overlap
1531  ctensorT co; rtensorT eo;
1532  syev(overlap, co, eo);
1533 
1534  if (_world.rank() == 0)
1535  {
1536  print("fock eigenvectors dims:",c.dim(0),c.dim(1));
1537  print("fock eigenvectors:");
1538  print(c);
1539  }
1540 
1541  if (_world.rank() == 0)
1542  {
1543  print("kinetic eigenvalues");
1544  print(ek);
1545  }
1546 
1547  if (_world.rank() == 0)
1548  {
1549  print("potential eigenvalues");
1550  print(ep);
1551  }
1552 
1553  if (_world.rank() == 0)
1554  {
1555  print("overlap eigenvalues");
1556  print(eo);
1557  }
1558  if (_world.rank() == 0)
1559  {
1560  print("fock eigenvalues");
1561  print(e);
1562  }
1563 
1564  if (_world.rank() == 0)
1565  {
1566  printf("Overlap: \n");
1567  for (int i = 0; i < overlap.dim(0); i++)
1568  {
1569  for (int j = 0; j < overlap.dim(1); j++)
1570  {
1571  printf("%10.5f", real(overlap(i,j)));
1572  }
1573  printf("\n");
1574  }
1575  printf("\n");
1576  printf("\n");
1577  }
1578 
1579  compress(_world, ao);
1580  _world.gop.fence();
1581  // Take linear combinations of the gaussian basis orbitals as the starting
1582  // orbitals for solver
1583  vecfuncT tmp_orbitals = transform(_world, ao, c(_, Slice(0, _nao - 1)));
1584 
1585 // tmp_orbitals[0] = ao[0];
1586 // tmp_orbitals[1] = ao[1];
1587 // tmp_orbitals[2] = ao[2];
1588 // tmp_orbitals[3] = ao[3];
1589 // tmp_orbitals[4] = ao[4];
1590 
1591  // WSTHORNTON
1592  // checking the symmetry
1593 // int npts = 101;
1594 // double delta = _params.L/(npts-1);
1595 // double begin = -_params.L/2;
1596 // double end = _params.L/2;
1597 // double tol = 1e-6;
1598 //
1599 // for (int i = 0; i < npts; i++)
1600 // {
1601 // double x = (i*delta)+begin;
1602 // for (int j = 0; j < npts; j++)
1603 // {
1604 // double y = (j*delta)+begin;
1605 // for (int k = 0; k < npts; k++)
1606 // {
1607 // double z = (k*delta)+begin;
1608 // coord_3d r1 {x,y,z};
1609 // coord_3d r2 {y,x,z};
1610 // coord_3d r3 {z,y,x};
1611 // coord_3d r4 {x,z,y};
1612 // double_complex pxr2 = ao[2](r2); double_complex pyr1 = ao[3](r1);
1613 // double_complex pxr3 = ao[2](r3); double_complex pzr1 = ao[4](r1);
1614 // double_complex pyr4 = ao[3](r4);
1615 // double err_xy = std::abs(pxr2-pyr1);
1616 // double err_xz = std::abs(pxr3-pzr1);
1617 // std::string success = err_xy < tol ? "PASS!" : "FAIL!";
1618 // printf("%10.4e + %10.4e %10.4e + %10.4e %15.8e\n %s\n",
1619 // std::real(pxr2), std::imag(pxr2), std::real(pyr1), std::imag(pyr1),
1620 // err_xy, success.c_str());
1621 // }
1622 // }
1623 // }
1624 
1625 // MADNESS_ASSERT(false);
1626 
1627  // WSTHORNTON
1628 // const double PI_OVER_8 = madness::constants::pi/8;
1629 // SeparatedConvolution<double,3> op = CoulombOperator(_world, _params.lo, _params.thresh * 0.1);
1630 // rfunctionT vlda = make_lda_potential(_world, rho, rho, rfunctionT(), rfunctionT());
1631 // rho.scale(2.0);
1632 // rfunctionT vc2 = apply(op,rho);
1633 // double rtr2 = rho.trace();
1634 // rfunctionT vlocal2 = vc2 + _vnuc;
1635 //
1636 // SeparatedConvolution<double,3> op2 = CoulombOperator(_world, _params.lo, _params.thresh * 0.1);
1637 // SeparatedConvolution<double,3> bshop = BSHOperator3D(_world, -0.2, _params.lo, _params.thresh * 0.1);
1638 
1639 
1640 // if (_world.rank() == 0)
1641 // {
1642 // print("trace of rho:");
1643 // print(rtr2);
1644 // }
1645 // for (int i = 0; i < 5; i++)
1646 // {
1647 // double th = i*PI_OVER_8;
1648 // tmp_orbitals[2] = ao[2];
1649 // tmp_orbitals[3] = cos(th)*ao[3] + sin(th)*ao[4];
1650 // tmp_orbitals[4] = -sin(th)*ao[3] + cos(th)*ao[4];
1651 //
1652 // ctensorT fred(3,3);
1653 // for (int ii = 2; ii < 5; ii++)
1654 // {
1655 // for (int jj = 2; jj < 5; jj++)
1656 // {
1657 // fred(ii-2,jj-2) = inner(tmp_orbitals[ii], _vnuc*tmp_orbitals[jj]);
1658 // if (std::abs(fred(ii-2,jj-2)) < 1e-6)
1659 // fred(ii-2,jj-2) = std::complex<double>(0.0,0.0);
1660 //
1661 // }
1662 // }
1663 
1664 // complex_derivative_3d Dx(_world,0);
1665 // complex_derivative_3d Dy(_world,1);
1666 // complex_derivative_3d Dz(_world,2);
1667 //
1668 // ctensorT fred(3,3);
1669 // std::complex<double> I = std::complex<double>(0.0,1.0);
1670 // std::complex<double> one = std::complex<double>(1.0,0.0);
1671 // for (int ii = 2; ii < 5; ii++)
1672 // {
1673 // for (int jj = 2; jj < 5; jj++)
1674 // {
1675 // double ksq = kpoint.k[0]*kpoint.k[0] +
1676 // kpoint.k[1]*kpoint.k[1] +
1677 // kpoint.k[2]*kpoint.k[2];
1678 //
1679 // cfunctionT ct1 = -I*kpoint.k[0]*Dx(tmp_orbitals[jj]) -
1680 // I*kpoint.k[1]*Dy(tmp_orbitals[jj]) -
1681 // I*kpoint.k[2]*Dz(tmp_orbitals[jj]);
1682 //
1683 // cfunctionT ct2 = -0.5*Dx(Dx(tmp_orbitals[jj])) +
1684 // Dy(Dy(tmp_orbitals[jj])) +
1685 // Dz(Dz(tmp_orbitals[jj]));
1686 //
1687 // cfunctionT ct3 = 0.5*ksq*tmp_orbitals[jj];
1688 //
1689 // fred(ii-2,jj-2) = inner(tmp_orbitals[ii], ct1+ct2+ct3);
1690 // if (std::abs(fred(ii-2,jj-2)) < 1e-6)
1691 // fred(ii-2,jj-2) = std::complex<double>(0.0,0.0);
1692 // }
1693 // }
1694 
1695 // ctensorT cf; rtensorT ef;
1696 // syev(fred, cf, ef);
1697 //
1698 // if (_world.rank() == 0)
1699 // {
1700 // print("Fred:", i);
1701 // print(fred);
1702 // print("eigenvalues:");
1703 // print(ef);
1704 // print("eigenvectors:");
1705 // print(cf);
1706 // }
1707 // }
1708 // MADNESS_ASSERT(false);
1709 
1710 
1711  _world.gop.fence();
1712  truncate(_world, tmp_orbitals);
1713  normalize(_world, tmp_orbitals);
1714 
1715  // Build the overlap matrix
1716  if (_world.rank() == 0) print("Building overlap matrix ...\n\n");
1717  ctensorT overlap2 = matrix_inner(_world, tmp_orbitals, tmp_orbitals, true);
1718 
1719  rtensorT tmp_eigs = e(Slice(0, _nao - 1));
1720 
1721  if (_world.rank() == 0) printf("(%8.4f,%8.4f,%8.4f)\n",kpoint.k[0], kpoint.k[1], kpoint.k[2]);
1722  if (_world.rank() == 0) print(tmp_eigs);
1723  if (_world.rank() == 0) print("\n");
1724 
1725  //if (_world.rank() == 0) print("kinetic energy for kp = ", kp);
1726  //if (_world.rank() == 0) print(kinetic);
1727  //if (_world.rank() == 0) print("\n");
1728 
1729  // DEBUG
1730  if (_world.rank() == 0) {
1731  printf("Overlap: \n");
1732  for (int i = 0; i < kinetic.dim(0); i++)
1733  {
1734  for (int j = 0; j < kinetic.dim(1); j++)
1735  {
1736  printf("%10.5f", real(overlap(i,j)));
1737  }
1738  printf("\n");
1739  }
1740  printf("\n");
1741  printf("\n");
1742 
1743  printf("Kinetic: \n");
1744  for (int i = 0; i < kinetic.dim(0); i++)
1745  {
1746  for (int j = 0; j < kinetic.dim(1); j++)
1747  {
1748  printf("%10.5f", real(kinetic(i,j)));
1749  }
1750  printf("\n");
1751  }
1752  printf("\n");
1753  printf("\n");
1754 
1755  printf("V: \n");
1756  for (int i = 0; i < potential.dim(0); i++)
1757  {
1758  for (int j = 0; j < potential.dim(1); j++)
1759  {
1760  printf("%10.5f", real(potential(i,j)));
1761  }
1762  printf("\n");
1763  }
1764  printf("\n");
1765  printf("\n");
1766 
1767  printf("Fock: \n");
1768  for (int i = 0; i < fock.dim(0); i++)
1769  {
1770  for (int j = 0; j < fock.dim(1); j++)
1771  {
1772  printf("%10.5f", real(fock(i,j)));
1773  }
1774  printf("\n");
1775  }
1776  printf("\n");
1777  printf("\n");
1778 
1779  printf("New overlap: \n");
1780  for (int i = 0; i < overlap2.dim(0); i++)
1781  {
1782  for (int j = 0; j < overlap2.dim(1); j++)
1783  {
1784  printf("%10.5f", real(overlap2(i,j)));
1785  }
1786  printf("\n");
1787  }
1788  printf("\n");
1789  printf("\n");
1790  }
1791 
1792  // Fill in orbitals and eigenvalues
1793  kpoint.begin = ki*_params.nbands;
1794  kpoint.end = (ki+1)*_params.nbands;
1795  for (unsigned int oi = kpoint.begin, ti = 0; oi < kpoint.end; oi++, ti++)
1796  {
1797  //if (_world.rank() == 0) print(oi, ti, kpt.begin, kpt.end);
1798  // normalize the orbitals
1799  //tmp_orbitals[ti].scale(1.0/tmp_orbitals[ti].norm2());
1800  _phisa.push_back(tmp_orbitals[ti]);
1801  _phisb.push_back(tmp_orbitals[ti]);
1802  _eigsa[oi] = tmp_eigs[ti];
1803  _eigsb[oi] = tmp_eigs[ti];
1804  }
1805  }
1806  }
1807  }
1808  //*************************************************************************
1809 
1810  //*************************************************************************
1811  // Constructor
1812  Solver(World& world,
1813  const rfunctionT& vnucrhon,
1814  const vecfuncT& phis,
1815  const std::vector<T>& eigs,
1816  const ElectronicStructureParams& params,
1817  MolecularEntity mentity)
1818  : _world(world), _vnucrhon(vnucrhon), _phisa(phis), _phisb(phis),
1819  _eigsa(eigs), _eigsb(eigs), _params(params), _mentity(mentity)
1820  {
1821  _residual = 1e5;
1822  if (params.periodic)
1823  {
1825  //_cop = PeriodicCoulombOpPtr<T,NDIM>(const_cast<World&>(world),
1826  // FunctionDefaults<NDIM>::get_k(), params.lo, params.thresh * 0.1, box);
1827  _cop = CoulombOperatorPtr(const_cast<World&>(world), params.lo, params.thresh * 0.1);
1828  }
1829  else
1830  {
1831  _cop = CoulombOperatorPtr(const_cast<World&>(world), params.lo, params.thresh * 0.1);
1832  }
1833 
1834  if (params.ispotential)
1835  {
1836  _vnuc = copy(_vnucrhon);
1837  }
1838  else
1839  {
1840  _vnuc = apply(*_cop, _vnucrhon);
1841  }
1842  }
1843  //*************************************************************************
1844 
1845  //*************************************************************************
1846  // Constructor
1847  Solver(World& world,
1848  rfunctionT vnucrhon,
1849  vecfuncT phis,
1850  std::vector<T> eigs,
1851  std::vector<KPoint> kpoints,
1852  std::vector<double> occs,
1854  MolecularEntity mentity)
1855  : _world(world), _vnucrhon(vnucrhon), _phisa(phis), _phisb(phis),
1856  _eigsa(eigs), _eigsb(eigs), _params(params),
1857  _kpoints(kpoints), _occsa(occs), _occsb(occs), _mentity(mentity)
1858  {
1859  _residual = 1e5;
1860  if (params.periodic)
1861  {
1863  _cop = CoulombOperatorPtr(const_cast<World&>(world), params.lo, params.thresh * 0.1);
1864  }
1865  else
1866  {
1867  _cop = CoulombOperatorPtr(const_cast<World&>(world), params.lo, params.thresh * 0.1);
1868  }
1869 
1870  if (params.ispotential)
1871  {
1872  _vnuc = copy(_vnucrhon);
1873  }
1874  else
1875  {
1876  _vnuc = apply(*_cop, _vnucrhon);
1877  }
1878  }
1879  //*************************************************************************
1880 
1881  //*************************************************************************
1882  virtual ~Solver()
1883  {
1884  _outputF.close();
1885  _matF.close();
1886  _eigF.close();
1887  delete _subspace;
1888  }
1889  //*************************************************************************
1890 
1891  //***************************************************************************
1892  // set occupation numbers (only for insulators ... no smearing)
1893  void set_occs2(const std::vector<KPoint>& kpoints,
1894  const std::vector<double>& eigsa,
1895  const std::vector<double>& eigsb,
1896  std::vector<double>& occsa,
1897  std::vector<double>& occsb)
1898  {
1899  for (unsigned int ik = 0; ik < kpoints.size(); ik++)
1900  {
1901  // get k-point
1902  KPoint k = kpoints[ik];
1903 
1904  // pull subset of data that corresponds to k
1905  const std::vector<double> k_eigsa(eigsa.begin() + k.begin, eigsa.begin() + k.end);
1906  const std::vector<double> k_eigsb(eigsb.begin() + k.begin, eigsb.begin() + k.end);
1907  std::vector<double> k_occsa(occsa.begin() + k.begin, occsa.begin() + k.end);
1908  std::vector<double> k_occsb(occsb.begin() + k.begin, occsb.begin() + k.end);
1909 
1910  // demand all vectors have the same size
1911  unsigned int sz = k_eigsa.size();
1912  MADNESS_ASSERT(k_eigsb.size() == sz);
1913  MADNESS_ASSERT(k_occsa.size() == sz);
1914  MADNESS_ASSERT(k_occsb.size() == sz);
1915  // concatenate eigenvalues
1916  std::vector<double> teigs;
1917  //std::copy(k_eigsa.begin(), k_eigsa.end(), teigs.begin());
1918  //teigs.insert(teigs.end(), k_eigsb.begin(), k_eigsb.end());
1919  //std::copy(k_eigsb.begin(), k_eigsb.end(), back_inserter(teigs));
1920  for (unsigned int ist = 0; ist < k_eigsa.size(); ist++) teigs.push_back(k_eigsa[ist]);
1921  for (unsigned int ist = 0; ist < k_eigsb.size(); ist++) teigs.push_back(k_eigsb[ist]);
1922 
1923  if (_world.rank() == 0) printf("setting occs ....\n\n");
1924  for (unsigned int ist = 0; ist < teigs.size(); ist++)
1925  {
1926  if (_world.rank() == 0)
1927  {
1928  printf("%5d %15.8f\n", ist, teigs[ist]);
1929  }
1930  }
1931 
1932  // indicies
1933  std::vector<unsigned int> inds(2*sz);
1934  for (unsigned int i = 0; i < 2*sz; i++) inds[i] = i;
1935  // sort by eigenvalue
1936  for (unsigned int i = 0; i < 2*sz; i++)
1937  {
1938  for (unsigned int j = i+1; j < 2*sz; j++)
1939  {
1940  if (teigs[j] < teigs[i])
1941  {
1942  double t1 = teigs[i];
1943  teigs[i] = teigs[j];
1944  teigs[j] = t1;
1945  int it1 = inds[i];
1946  inds[i] = inds[j];
1947  inds[j] = it1;
1948  }
1949  }
1950  }
1951 
1952  if (_world.rank() == 0)
1953  printf("\nSorted eigenvalues:\n");
1954  for (unsigned int i = 0; i < teigs.size(); i++)
1955  {
1956  if (_world.rank() == 0)
1957  printf("%10d%10d %15.8f\n",i,inds[i],teigs[i]);
1958  }
1959 
1960  // assign occupation numbers
1961  double availablecharge = _params.ncharge;
1962  for (unsigned int i = 0; (i < 2*sz) && (availablecharge > 0.0) ; i++)
1963  {
1964  unsigned int current = inds[i];
1965  if (current >= sz)
1966  {
1967  current -= sz;
1968  k_occsb[current] = 1.0;
1969  availablecharge -= 1.0;
1970  }
1971  else
1972  {
1973  k_occsa[current] = 1.0;
1974  availablecharge -= 1.0;
1975  }
1976  }
1977 
1978  for (unsigned int ik1 = k.begin, ik2 = 0; ik1 < k.end; ik1++,ik2++)
1979  {
1980  occsa[ik1] = k_occsa[ik2];
1981  occsb[ik1] = k_occsb[ik2];
1982  }
1983  }
1984 
1985  for (unsigned int ik = 0; ik < kpoints.size(); ik++)
1986  {
1987  KPoint k = kpoints[ik];
1988  if (_world.rank() == 0)
1989  {
1990  printf("k-point is: %d: \n",ik);
1991  }
1992  for (unsigned int ist = k.begin; ist < k.end; ist++)
1993  {
1994  if (_world.rank() == 0)
1995  {
1996  printf("occa: %12.5f occb: %12.5f \n",occsa[ist],occsb[ist]);
1997  }
1998  }
1999  }
2000 
2001  }
2002  //***************************************************************************
2003 
2004  //***************************************************************************
2005  /*!
2006  \ingroup periodic_solver
2007  \brief Compute the electronic density for either a molecular or periodic
2008  system.
2009  */
2010  rfunctionT compute_rho_slow(const vecfuncT& phis, std::vector<KPoint> kpoints,
2011  std::vector<double> occs)
2012  {
2013  // Electron density
2014  rfunctionT rho = rfactoryT(_world);
2015  _world.gop.fence();
2016  if (_world.rank() == 0) _outputF << "computing rho ..." << endl;
2017  // Loop over k-points
2018  for (unsigned int kp = 0; kp < kpoints.size(); kp++)
2019  {
2020  // get k-point
2021  KPoint kpoint = kpoints[kp];
2022  // loop through bands
2023  for (unsigned int j = kpoint.begin; j < kpoint.end; j++)
2024  {
2025  // Get phi(j) from iterator
2026  const functionT& phij = phis[j];
2027  // Compute the j-th density
2028  //rfunctionT prod = abs_square(phij);
2029  rfunctionT prod = abssq(phij,true);
2030  double rnrm = prod.trace();
2031  prod.scale(1/rnrm);
2032 // rho += 0.5 * _occs[j] * kpoint.weight * prod;
2033  rho += occs[j] * kpoint.weight * prod;
2034  }
2035  }
2036  rho.truncate();
2037 
2038  return rho;
2039  }
2040 
2041  // computes the electronic density for 1 spin
2042  rfunctionT compute_rho(const vecfuncT& phis, std::vector<KPoint> kpoints,
2043  std::vector<double> occs)
2044  {
2045  if (_world.rank() == 0) _outputF << "computing rho ..." << endl;
2046  rfunctionT rho = rfactoryT(_world); // Electron density
2047 
2048  reconstruct(_world, phis); // For max parallelism
2049  std::vector<rfunctionT> phisq(phis.size());
2050  for (unsigned int i=0; i<phis.size(); i++) {
2051  phisq[i] = abssq(phis[i],false);
2052  }
2053  _world.gop.fence();
2054  std::vector<double> phinorm = norm2s(_world, phis);
2055 
2056  compress(_world,phisq); // since will be using gaxpy for accumulation
2057  rho.compress();
2058 
2059  // Loop over k-points
2060  for (unsigned int kp = 0; kp < kpoints.size(); kp++)
2061  {
2062  // get k-point
2063  KPoint kpoint = kpoints[kp];
2064  // loop through bands
2065  for (unsigned int j = kpoint.begin; j < kpoint.end; j++)
2066  {
2067  rho.gaxpy(1.0, phisq[j], occs[j] * kpoint.weight / (phinorm[j]*phinorm[j]), false);
2068  }
2069  }
2070  _world.gop.fence();
2071  phisq.clear();
2072  rho.truncate();
2073 
2074  return rho;
2075  }
2076 
2077 
2078  //***************************************************************************
2079 
2080  //***************************************************************************
2081  std::vector<poperatorT> make_bsh_operators(const std::vector<T>& eigs)
2082  {
2083  // Make BSH vector
2084  std::vector<poperatorT> bops;
2085  // Get defaults
2086  double tol = FunctionDefaults<NDIM>::get_thresh();
2087  // Loop through eigenvalues, adding a BSH operator to bops
2088  // for each eigenvalue
2089  int sz = eigs.size();
2090  for (int i = 0; i < sz; i++)
2091  {
2092  T eps = eigs[i];
2093  if (eps > 0)
2094  {
2095  if (_world.rank() == 0)
2096  {
2097  std::cout << "bsh: warning: positive eigenvalue" << i << eps << endl;
2098  }
2099  eps = -0.1;
2100  }
2101 
2102  bops.push_back(poperatorT(BSHOperatorPtr3D(_world, sqrt(-2.0*eps), _params.lo, tol * 0.1)));
2103  }
2104  return bops;
2105  }
2106  //*************************************************************************
2107 
2108  //*************************************************************************
2109 // void loadbal()
2110 // {
2111 // if(_world.size() == 1)
2112 // return;
2113 //
2114 // LoadBalanceDeux<3> lb(_world);
2115 // lb.add_tree(_vnuc, lbcost<double,3>(1.0, 0.0), false);
2116 // lb.add_tree(_rhoa, lbcost<double,3>(1.0, 1.0), false);
2117 // for(unsigned int i = 0;i < _phisa.size();i++)
2118 // {
2119 // lb.add_tree(_phisa[i], lbcost<valueT,3>(1.0, 1.0), false);
2120 // }
2121 // if(_params.spinpol)
2122 // {
2123 // lb.add_tree(_rhob, lbcost<double,3>(1.0, 1.0), false);
2124 // for(unsigned int i = 0;i < _phisb.size();i++)
2125 // {
2126 // lb.add_tree(_phisa[i], lbcost<valueT,3>(1.0, 1.0), false);
2127 // }
2128 // }
2129 //
2130 // FunctionDefaults<3>::redistribute(_world, lb.load_balance(6.0));
2131 // }
2132  //*************************************************************************
2133 
2134  //*************************************************************************
2136  {
2137  double ke = 0.0;
2138  if (!_params.periodic)
2139  {
2143  for (unsigned int i = 0; i < _phisa.size(); i++)
2144  {
2145  functionT dpdx = Dx(_phisa[i]);
2146  functionT dpdy = Dy(_phisa[i]);
2147  functionT dpdz = Dz(_phisa[i]);
2148  ke += 0.5 * (real(inner(dpdx,dpdx)) + real(inner(dpdy,dpdy))
2149  + real(inner(dpdz,dpdz)));
2150  }
2151  if (_params.spinpol)
2152  {
2153  for (unsigned int i = 0; i < _phisb.size(); i++)
2154  {
2155  functionT dpdx = Dx(_phisb[i]);
2156  functionT dpdy = Dy(_phisb[i]);
2157  functionT dpdz = Dz(_phisb[i]);
2158  ke += 0.5 * (real(inner(dpdx,dpdx)) + real(inner(dpdy,dpdy))
2159  + real(inner(dpdz,dpdz)));
2160  }
2161  }
2162  else
2163  {
2164  ke *= 2.0;
2165  }
2166  }
2167  return ke;
2168  }
2169  //*************************************************************************
2170 
2171  //*************************************************************************
2172  /*!
2173  \ingroup periodic_solver
2174  \brief Applies the LDA effective potential to each orbital. Currently only
2175  lda and spin-polarized is not implemented.
2176  */
2177  void apply_potential(vecfuncT& pfuncsa,
2178  vecfuncT& pfuncsb, const vecfuncT& phisa,
2179  const vecfuncT& phisb, const rfunctionT& rhoa, const rfunctionT& rhob,
2180  const rfunctionT& rho)
2181  {
2182  // Nuclear and coulomb potentials
2183  rfunctionT vc = apply(*_cop, rho);
2184 
2185  // combined density
2186  rfunctionT rho2 = rho + _vnucrhon;
2187  double vnucrhontrace = _vnucrhon.trace();
2188  double rhotrace = rho.trace();
2189  double rho2trace = rho2.trace();
2190  if (_world.rank() == 0) printf("_vnucrhon trace: %10e\n", vnucrhontrace);
2191  if (_world.rank() == 0) printf("rho trace: %10e\n", rhotrace);
2192  if (_world.rank() == 0) printf("rho2 trace: %10e\n", rho2trace);
2193  rfunctionT vlocal = (_params.ispotential) ? _vnuc + vc : apply(*_cop, rho2);
2194  rfunctionT vlocal2 = _vnuc + vc;
2195  double vlerr = (vlocal-vlocal2).norm2();
2196  if (_world.rank() == 0) printf("vlerr trace: %10e\n\n", vlerr);
2197 
2198 
2199  // Calculate energies for Coulomb and nuclear
2200  double ce = 0.5*inner(vc,rho);
2201  double pe = inner(_vnuc,rho);
2202  double xc = 0.0;
2203  double ke = calculate_kinetic_energy();
2204  // Exchange
2205  if (_params.functional == 1)
2206  {
2207  // LDA, is calculation spin-polarized?
2208  if (_params.spinpol)
2209  {
2210  MADNESS_EXCEPTION("Spin polarized not implemented!",0);
2211 // // potential
2212 // rfunctionT vxca = binary_op(rhoa, rhob, &::libxc_ldaop_sp<double>);
2213 // rfunctionT vxcb = binary_op(rhob, rhoa, &::libxc_ldaop_sp<double>);
2214 // pfuncsa = mul_sparse(_world, vlocal + vxca, phisa, _params.thresh * 0.1);
2215 // pfuncsb = mul_sparse(_world, vlocal + vxcb, phisb, _params.thresh * 0.1);
2216 // // energy
2217 // rfunctionT fca = binary_op(rhoa, rhob, &::libxc_ldaeop_sp<double>);
2218 // rfunctionT fcb = binary_op(rhob, rhoa, &::libxc_ldaeop_sp<double>);
2219 // xc = fca.trace() + fcb.trace();
2220  }
2221  else
2222  {
2223  // potential
2224  rfunctionT vxc = copy(rhoa);
2225 // vxc.unaryop(&::libxc_ldaop<double>);
2227  vxc.unaryop(&::ldaop<double>);
2228 
2229  //test_periodicity(vc);
2230 // for (unsigned int i = 0; i < phisa.size(); i++)
2231 // {
2232 // test_periodicity(phisa[i]);
2233 // }
2234 
2235  pfuncsa = mul_sparse(_world, vlocal2+vxc, phisa, _params.thresh * 0.1);
2236 // rfunctionT vxc2 = binary_op(rhoa, rhoa, &::libxc_ldaop_sp<double>);
2237 // pfuncsa = mul_sparse(_world, vlocal + vxc2, phisa, _params.thresh * 0.1);
2238 
2239 // for (unsigned int i = 0; i < pfuncsa.size(); i++)
2240 // {
2241 // test_periodicity(pfuncsa[i]);
2242 // }
2243 
2244 
2245  END_TIMER(_world,"Applying LDA potential");
2246  // energy
2247  rfunctionT fc = copy(rhoa);
2248  fc.unaryop(&::ldaeop<double>);
2249  xc = fc.trace();
2250  }
2251  }
2252  else if (_params.functional == 2)
2253  {
2255  pfuncsa = mul_sparse(_world, vlocal, phisa, _params.thresh * 0.1);
2256  END_TIMER(_world,"Applying local potential");
2258  // gamma-point?
2259  if ((_params.ngridk0 == 1) && (_params.ngridk1 == 1) && (_params.ngridk2 == 1))
2260  {
2261  apply_hf_exchange3(_phisa, _phisb, pfuncsa, pfuncsb, xc);
2262  }
2263  else
2264  {
2265  apply_hf_exchange4(_phisa, _phisb, pfuncsa, pfuncsb, xc);
2266  }
2267  END_TIMER(_world, "Applying HF exchange");
2268  }
2269  std::cout.precision(8);
2270  if (_world.rank() == 0)
2271  {
2272  printf("Energies:\n");
2273  printf("Kinetic energy: %20.10f\n", ke);
2274  printf("Potential energy: %20.10f\n", pe);
2275  printf("Coulomb energy: %20.10f\n", ce);
2276  printf("Exchange energy: %20.10f\n", xc);
2277  printf("Total energy: %20.10f\n\n", ke + pe + ce + xc);
2278  }
2279  }
2280  //*************************************************************************
2281 
2282  // working for a molecule / atom ... need to test with periodic boundary
2283  // conditions as a gamma point only calculation
2284  //*************************************************************************
2285 // void apply_hf_exchange2(vecfuncT& phisa, vecfuncT& phisb,
2286 // vecfuncT& funcsa, vecfuncT& funcsb,
2287 // double& xc)
2288 // {
2289 // Vector<double,3> q {0.0,0.0,0.0};
2290 // SeparatedConvolution<double_complex,3> hfexop =
2291 // PeriodicHFExchangeOperator(_world, q, _params.lo,
2292 // FunctionDefaults<3>::get_thresh() * 0.1);
2293 // for (unsigned int i = 0; i < phisa.size(); i++)
2294 // {
2295 // bool isreal1 = is_real<T,NDIM>(phisa[i]);
2296 // MADNESS_ASSERT(isreal1);
2297 //// rfunctionT phi_i = real(phisa[i]);
2298 // for (unsigned int j = 0; j < phisa.size(); j++)
2299 // {
2300 // bool isreal2 = is_real<T,NDIM>(phisa[j]);
2301 // MADNESS_ASSERT(isreal2);
2302 //// rfunctionT phi_j = real(phisa[j]);
2303 //// rfunctionT prod = phi_i*phi_j;
2304 // functionT prod = phisa[i]*phisa[j];
2305 // prod.truncate();
2306 //// rfunctionT Vex = apply(*_cop,prod);
2307 // rfunctionT Vex = real(apply(hfexop,prod));
2308 // functionT tf1 = Vex*phisa[j];
2309 // funcsa[i] -= tf1;
2310 // xc -= real(inner(phisa[i],tf1));
2311 // }
2312 // }
2313 // }
2314  //*************************************************************************
2315 
2316  // working for a molecule / atom ... works for a gamma point calculation
2317  //*************************************************************************
2318  void apply_hf_exchange3(vecfuncT& phisa, vecfuncT& phisb,
2319  vecfuncT& funcsa, vecfuncT& funcsb,
2320  double& xc)
2321  {
2322  for (unsigned int j = 0; j < phisa.size(); j++)
2323  {
2324  rfunctionT phi_j = real(phisa[j]);
2325  // do diagonal piece first
2326  rfunctionT dprod = phi_j*phi_j;
2327  dprod.truncate();
2328  rfunctionT dVex = apply(*_cop,dprod);
2329  functionT tf_jjj = dVex*phisa[j];
2330  funcsa[j] -= tf_jjj;
2331  xc -= real(inner(phisa[j],tf_jjj));
2332  for (unsigned int i = j+1; i < phisa.size(); i++)
2333  {
2334  rfunctionT phi_i = real(phisa[i]);
2335  rfunctionT prod = phi_i*phi_j;
2336  prod.truncate();
2337  rfunctionT Vex = apply(*_cop,prod);
2338  // do the jij-th term
2339  functionT tf_jij = Vex*phisa[j];
2340  funcsa[i] -= tf_jij;
2341  xc -= real(inner(phisa[i],tf_jij));
2342  // do the iji-th term (in the complex case, use the complex
2343  // conjugate)
2344  functionT tf_iji = Vex*phisa[i];
2345  funcsa[j] -= tf_iji;
2346  xc -= real(inner(phisa[j],tf_iji));
2347  }
2348  }
2349  }
2350  //*************************************************************************
2351 
2352  //*************************************************************************
2353  KPoint find_kpt_from_orb(unsigned int idx)
2354  {
2355  for (unsigned int i = 0; i < _kpoints.size(); i++)
2356  {
2357  KPoint k1 = _kpoints[i];
2358  if (k1.is_orb_in_kpoint(idx)) return k1;
2359  }
2360  MADNESS_EXCEPTION("Error: find_kpt_from_orb: didn't find kpoint\n", 0);
2361  }
2362  //************************************************************************
2363 
2364  // hf exchange with k-points for periodic solid
2365  //*************************************************************************
2366  void apply_hf_exchange4(vecfuncT& phisa, vecfuncT& phisb,
2367  vecfuncT& funcsa, vecfuncT& funcsb,
2368  double& xc)
2369  {
2370  for (unsigned int i = 0; i < phisa.size(); i++)
2371  {
2372  functionT phi_i = phisa[i];
2373  KPoint k_i = find_kpt_from_orb(i);
2374  for (unsigned int j = 0; j < phisa.size(); j++)
2375  {
2376  functionT phi_j = phisa[j];
2377  KPoint k_j = find_kpt_from_orb(j);
2378  Vector<double,3> q {(k_i.k[0]-k_j.k[0])*_params.L,
2379  (k_i.k[1]-k_j.k[1])*_params.L,
2380  (k_i.k[2]-k_j.k[2])*_params.L};
2381  Vector<double,3> q2 {k_i.k[0]-k_j.k[0],
2382  k_i.k[1]-k_j.k[1],
2383  k_i.k[2]-k_j.k[2]};
2384  functionT cexp =
2385  factoryT(_world).functor(functorT(new
2386  ComplexExp<3>(q2, double_complex(1.0,0.0))));
2387  functionT cexp2 = conj(cexp);
2388  cexp.truncate();
2389  cexp2.truncate();
2390  functionT prod = phi_i*phi_j*cexp;
2394  functionT Vex = apply(hfexop,prod);
2395  functionT tf1 = Vex*phisa[j]*cexp2;
2396  funcsa[i] -= tf1;
2397  xc -= real(inner(phisa[i],tf1));
2398  }
2399  }
2400  }
2401  //*************************************************************************
2402 
2403  // not working and not tested .... supposed to be for the periodic boundary
2404  // conditions with k-points
2405  //*************************************************************************
2406  void apply_hf_exchange(vecfuncT& phisa, vecfuncT& phisb,
2407  vecfuncT& funcsa, vecfuncT& funcsb)
2408  {
2409  for (unsigned int ink1 = 0, ik1 = 0; ink1 < _phisa.size(); ++ink1)
2410  {
2411  for (unsigned int ink2 = 0, ik2 = 0; ink2 <= ink1; ink2++)
2412  {
2413  KPoint k1 = _kpoints[ik1];
2414  KPoint k2 = _kpoints[ik2];
2415 
2416  if (ink1 == k1.end) ik1++;
2417  if (ink2 == k2.end) ik2++;
2418 
2419  MADNESS_ASSERT(ik1 == 0);
2420  MADNESS_ASSERT(ik2 == 0);
2421 
2422  // no phase factor
2423  if (ik1 == ik2)
2424  {
2425 // // same state (diagonal piece)
2426 // if (ink1 == ink2)
2427 // {
2428 // rfunctionT prod = abs_square(phisa[ink1]);
2429 // rfunctionT fr = apply(*_cop,prod);
2430 // funcsa[ink1] -= funcsa[ink1]*fr;
2431 // }
2432 // else
2433  {
2434  Vector<double,3> q = 0.0;
2435  functionT f = phisa[ink1]*conj(phisa[ink2]);
2438  functionT fr = apply(hfexop,f);
2439  funcsa[ink1] -= funcsa[ink2]*fr;
2440  funcsa[ink2] -= funcsa[ink1]*conj(fr);
2441 
2442  // test code
2443  rfunctionT g1 = abs(phisa[ink1]);
2444  rfunctionT g2 = abs(phisa[ink2]);
2445  rfunctionT ff = g1*g2;
2446 
2447  printf("norm diff: %20.8f\n", abs(ff.norm2()-f.norm2()));
2448  MADNESS_ASSERT(abs(ff.norm2()-f.norm2()) <= 1e-5);
2449  rfunctionT fr2 = apply(*_cop,ff);
2450  MADNESS_ASSERT(abs(fr.norm2()-fr2.norm2()) <= 1e-5);
2451  }
2452  }
2453 // else
2454 // {
2455 // Vector<double,3> q = VectorFactory(k1.k[0]-k2.k[0],
2456 // k1.k[1]-k2.k[1],
2457 // k1.k[2]-k2.k[2]);
2458 // functionT cexp = factoryT(_world).functor(functorT(new ComplexExp<3>(q, double_complex(1.0,0.0))));
2459 // cexp.truncate();
2460 //
2461 // functionT f = phisa[ink1]*conj(phisa[ink2])*cexp;
2462 // SeparatedConvolution<double_complex,3> hfexop =
2463 // PeriodicHFExchangeOperator(_world, q, _params.lo, FunctionDefaults<3>::get_thresh() * 0.1);
2464 // functionT fr = apply(hfexop,f);
2465 // funcsa[ink1] -= phisa[ink1]*fr*conj(cexp);
2466 // funcsa[ink2] -= phisa[ink2]*conj(fr)*cexp;
2467 // }
2468  }
2469  }
2470  }
2471  //*************************************************************************
2472 
2473  //*************************************************************************
2474  void reproject()
2475  {
2476  _params.waveorder += 2;
2477  _params.thresh /= 100;
2480  if (_world.rank() == 0) _outputF << "reprojecting to wavelet order "
2481  << _params.waveorder << endl;
2483  for(unsigned int i = 0; i < _phisa.size(); i++)
2484  {
2487  }
2488  _world.gop.fence();
2491  if(_params.spinpol)
2492  {
2494  for(unsigned int i = 0; i < _phisb.size(); i++)
2495  {
2498  }
2499  _world.gop.fence();
2502  }
2503 
2504  delete _cop;
2506  //_subspace->reproject();
2507  delete _subspace;
2509  }
2510  //*************************************************************************
2511 
2512  //*************************************************************************
2513  void solve()
2514  {
2515 
2516  // WSTHORNTON (debug) Test periodicity
2517 // if (_world.rank() == 0) printf("initial orbitals (periodicity) ...\n\n");
2518 // for (unsigned int i = 0; i < _phisa.size(); i++)
2519 // {
2520 // if (_world.rank() == 0) printf("orbital %d\n",i);
2521 // test_periodicity(_phisa[i]);
2522 // if (_world.rank() == 0) printf("\n\n");
2523 // }
2524 // if (_world.rank() == 0) printf("\n\n\n\n");
2525 
2526  if (_world.rank() == 0) print("size of phisa is: ", _phisa.size());
2527  // keep track of how many iterations have gone by without reprojecting
2528  int rit = 0;
2529  int rpthresh = 20;
2530  for (_it = 0; _it < _params.maxits && _residual > _params.rcriterion; _it++, rit++)
2531  {
2532  // should we reproject?
2533  if ((_it > 0) && ((_residual < 20*_params.thresh) || (rit == rpthresh)))
2534  {
2535  // reproject orbitals and reset threshold
2536  reproject();
2537  rit = 0;
2538  }
2539 
2540  if (_world.rank() == 0) _outputF << "_it = " << _it << endl;
2541 
2542  // Set occupation numbers
2544 
2545  // Compute density
2548 // rfunctionT rhoa = _rhoa;
2549 // rfunctionT rhob = _rhob;
2550  double drhoa = (_rhoa-rhoa).trace();
2551  double drhob = (_rhob-rhob).trace();
2552  if (_world.rank() == 0)
2553  {
2554  printf("diff of alpha rho: %15.7f\n",drhoa);
2555  printf("diff of beta rho: %15.7f\n",drhob);
2556  }
2557 
2558  _rho = rhoa + rhob;
2559 
2560  _rhoa = rhoa;
2561  _rhob = rhob;
2562  double rtrace = _rho.trace();
2563  if (_world.rank() == 0) _outputF << "trace of rho" << rtrace << endl;
2564 
2565 // if(_it < 2 || (_it % 10) == 0)
2566 // {
2567 // START_TIMER(_world);
2568 // //loadbal();
2569 // END_TIMER(_world, "Load balancing");
2570 // }
2571 
2572  std::vector<functionT> pfuncsa =
2573  zero_functions<valueT,NDIM>(_world, _phisa.size());
2574  std::vector<functionT> pfuncsb =
2575  zero_functions<valueT,NDIM>(_world, _phisb.size());
2576 
2577  // Apply the potentials to the orbitals
2578  if (_world.rank() == 0) _outputF << "applying potential ...\n" << endl;
2579  //START_TIMER(_world);
2580  apply_potential(pfuncsa, pfuncsb, _phisa, _phisb, _rhoa, _rhob, _rho);
2581  //END_TIMER(_world,"apply potential");
2582 
2583  // Do right hand side for all k-points
2584  std::vector<double> alpha(pfuncsa.size(), 0.0);
2585  do_rhs_simple(_phisa, pfuncsa, _kpoints, alpha, _eigsa);
2586 
2587  if (_params.plotorbs)
2588  {
2589  std::vector<long> npt(3,101);
2590  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2591  {
2592  KPoint kpoint = _kpoints[ik];
2593  int ist = 0;
2594  for (unsigned int kst = kpoint.begin; kst < kpoint.end; kst++, ist++)
2595  {
2596  std::ostringstream strm;
2597  strm << "pre_unk_" << ik << "_" << ist << ".dx";
2598  std::string fname = strm.str();
2599  plotdx(_phisa[kst], fname.c_str(), FunctionDefaults<3>::get_cell(), npt);
2600  }
2601  }
2602  }
2603  // WSTHORNTON
2604  // DEBUG
2605  if (_world.rank() == 0)
2606  printf("\n\n\n\n------ Debugging BSH operator -----");
2607  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2608  {
2609  KPoint kpoint = _kpoints[ik];
2610  std::vector<double> k_alpha(alpha.begin() + kpoint.begin, alpha.begin() + kpoint.end);
2611  if (_world.rank() == 0)
2612  printf("alpha: (%6.4f,%6.4f,%6.4f)\n",kpoint.k[0],kpoint.k[1],kpoint.k[2]);
2613  for (unsigned int ia = 0; ia < k_alpha.size(); ia++)
2614  {
2615  if (_world.rank() == 0) printf("%15.8f\n", k_alpha[ia]);
2616  }
2617  }
2618  if (_world.rank() == 0) printf("\n\n\n\n");
2619 
2620  // WSTHORNTON (debug)
2621  if (_world.rank() == 0) printf("before BSH application ...\n\n");
2622  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2623  {
2624  KPoint kpoint = _kpoints[ik];
2625  double k0 = kpoint.k[0];
2626  double k1 = kpoint.k[1];
2627  double k2 = kpoint.k[2];
2628  if (_world.rank() == 0)
2629  printf("(%6.5f, %6.5f, %6.5f)\n",k0,k1,k2);
2630  std::vector<functionT> k_phisa(_phisa.begin() + kpoint.begin, _phisa.begin() + kpoint.end);
2631  std::vector<functionT> k_pfuncsa(pfuncsa.begin() + kpoint.begin, pfuncsa.begin() + kpoint.end);
2632  KPoint kpoint_gamma;
2633  print_fock_matrix_eigs(k_phisa, k_pfuncsa, kpoint_gamma);
2634 
2635  // diagonalize overlap
2636  tensorT overlap = matrix_inner(_world,k_pfuncsa,k_pfuncsa,true);
2637  ctensorT co; rtensorT eo;
2638  syev(overlap, co, eo);
2639 
2640  if (_world.rank() == 0) printf("Overlap eigenvalues: \n");
2641  if (_world.rank() == 0) print(overlap);
2642  for (unsigned int ie = 0; ie < eo.dim(0); ie++)
2643  {
2644  if (_world.rank() == 0) printf("%d %15.8e\n", ie, eo(ie,ie));
2645  }
2646 
2647  }
2648  if (_world.rank() == 0) printf("\n\n\n\n");
2649 
2650  // WSTHORNTON (debug) Test periodicity
2651 // if (_world.rank() == 0) printf("before BSH application (periodicity) ...\n\n");
2652 // for (unsigned int i = 0; i < _phisa.size(); i++)
2653 // {
2654 // if (_world.rank() == 0) printf("orbital %d\n",i);
2655 // test_periodicity(_phisa[i]);
2656 // if (_world.rank() == 0) printf("\n\n");
2657 // }
2658 // if (_world.rank() == 0) printf("\n\n\n\n");
2659 
2660  // Make BSH Green's function
2661  std::vector<poperatorT> bopsa = make_bsh_operators(alpha);
2662  std::vector<T> sfactor(pfuncsa.size(), -2.0);
2663  scale(_world, pfuncsa, sfactor);
2664 
2665  // Apply Green's function to orbitals
2666  if (_world.rank() == 0) std::cout << "applying BSH operator ...\n" << endl;
2667  truncate<valueT,NDIM>(_world, pfuncsa);
2669  std::vector<functionT> tmpa = apply(_world, bopsa, pfuncsa);
2670  END_TIMER(_world,"apply BSH");
2671  bopsa.clear();
2672 
2673  // WSTHORNTON
2674  // norms
2675  if (_world.rank() == 0) printf("norms of tmpa\n");
2676  std::vector<double> tmpa_norms = norm2s(_world, tmpa);
2677  for (unsigned int i = 0; i < tmpa_norms.size(); i++)
2678  {
2679  if (_world.rank() == 0) printf("%10d %15.8f\n", i, tmpa_norms[i]);
2680  }
2681  if (_world.rank() == 0) printf("\n\n\n\n");
2682 
2683  // Do other spin
2684  vecfuncT tmpb = zero_functions<valueT,NDIM>(_world, _phisb.size());
2685  if (_params.spinpol)
2686  {
2687  alpha = std::vector<double>(_phisb.size(), 0.0);
2688  do_rhs_simple(_phisb, pfuncsb, _kpoints, alpha, _eigsb);
2689  std::vector<poperatorT> bopsb = make_bsh_operators(alpha);
2690  scale(_world, pfuncsb, sfactor);
2691  truncate<valueT,NDIM>(_world, pfuncsb);
2692  tmpb = apply(_world, bopsb, pfuncsb);
2693  bopsb.clear();
2694  }
2695  else
2696  {
2697  for (unsigned int i = 0; i < _eigsa.size(); i++) _eigsb[i] = _eigsa[i];
2698  }
2699 
2700  // WSTHORNTON (debug)
2701  std::vector<functionT> pfuncsa2=
2702  zero_functions<valueT,NDIM>(_world, _phisa.size());
2703  std::vector<functionT> pfuncsb2=
2704  zero_functions<valueT,NDIM>(_world, _phisa.size());
2705 
2706  // Apply the potentials to the orbitals
2707  if (_world.rank() == 0) _outputF << "applying potential2 ...\n" << endl;
2708  apply_potential(pfuncsa2, pfuncsb2, tmpa, tmpb, _rhoa, _rhob, _rho);
2709  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2710  {
2711  KPoint kpoint = _kpoints[ik];
2712  double k0 = kpoint.k[0];
2713  double k1 = kpoint.k[1];
2714  double k2 = kpoint.k[2];
2715  if (_world.rank() == 0)
2716  printf("(%6.5f, %6.5f, %6.5f)\n",k0,k1,k2);
2717  std::vector<functionT> k_tmpa(tmpa.begin() + kpoint.begin, tmpa.begin() + kpoint.end);
2718  std::vector<functionT> k_pfuncsa2(pfuncsa2.begin() + kpoint.begin, pfuncsa2.begin() + kpoint.end);
2719  print_fock_matrix_eigs(k_tmpa, k_pfuncsa2, kpoint);
2720 
2721 
2722  // diagonalize overlap
2723 // tensorT overlap = matrix_inner(_world,k_tmpa,k_tmpa,true);
2724 // if (_world.rank() == 0) print(overlap);
2725 // ctensorT co; rtensorT eo;
2726 // syev(overlap, co, eo);
2727 //
2728 // if (_world.rank() == 0) printf("Overlap eigenvalues: \n");
2729 // for (unsigned int ie = 0; ie < eo.dim(0); ie++)
2730 // {
2731 // if (_world.rank() == 0) printf("%d %15.8e\n", ie, eo(ie,ie));
2732 // }
2733  }
2734  if (_world.rank() == 0) printf("\n\n\n\n");
2735 
2736 
2737 
2738  if (_world.rank() == 0) printf("\n\n\n\n");
2739 
2740  // Update orbitals
2741  update_orbitals(tmpa, tmpb, _kpoints);
2742  save_orbitals();
2743 
2744  // WSTHORNTON
2745  if (_world.rank() == 0) _outputF << "after update_orbitals() ...\n" << endl;
2746  pfuncsa2=zero_functions<valueT,NDIM>(_world, _phisa.size());
2747  pfuncsb2=zero_functions<valueT,NDIM>(_world, _phisa.size());
2748  apply_potential(pfuncsa2, pfuncsb2, _phisa, _phisb, _rhoa, _rhob, _rho);
2749  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2750  {
2751  KPoint kpoint = _kpoints[ik];
2752  double k0 = kpoint.k[0];
2753  double k1 = kpoint.k[1];
2754  double k2 = kpoint.k[2];
2755  if (_world.rank() == 0)
2756  printf("(%6.5f, %6.5f, %6.5f)\n",k0,k1,k2);
2757  std::vector<functionT> k_phisa(_phisa.begin() + kpoint.begin, _phisa.begin() + kpoint.end);
2758  std::vector<functionT> k_pfuncsa2(pfuncsa2.begin() + kpoint.begin, pfuncsa2.begin() + kpoint.end);
2759  print_fock_matrix_eigs(k_phisa, k_pfuncsa2, kpoint);
2760  }
2761  if (_world.rank() == 0) printf("\n\n\n\n");
2762 
2763  }
2764 
2765  if (_params.plotorbs)
2766  {
2767  std::vector<long> npt(3,101);
2768  for (unsigned int ik = 0; ik < _kpoints.size(); ik++)
2769  {
2770  KPoint kpoint = _kpoints[ik];
2771  int ist = 0;
2772  for (unsigned int kst = kpoint.begin; kst < kpoint.end; kst++, ist++)
2773  {
2774  std::ostringstream strm;
2775  strm << "unk_" << ik << "_" << ist << ".dx";
2776  std::string fname = strm.str();
2777  plotdx(_phisa[kst], fname.c_str(), FunctionDefaults<3>::get_cell(), npt);
2778  }
2779  }
2780  }
2781  save_orbitals();
2782  }
2783  //*************************************************************************
2784 
2785  //*************************************************************************
2787  const double tol = 1e-13;
2788  MADNESS_ASSERT(A.dim(0) == A.dim(1));
2789 
2790  // Scale A by a power of 2 until it is "small"
2791  double anorm = A.normf();
2792  int n = 0;
2793  double scale = 1.0;
2794  while (anorm*scale > 0.1)
2795  {
2796  n++;
2797  scale *= 0.5;
2798  }
2799  tensorT B = scale*A; // B = A*2^-n
2800 
2801  // Compute exp(B) using Taylor series
2802  ctensorT expB = ctensorT(2, B.dims());
2803  for (int i = 0; i < expB.dim(0); i++) expB(i,i) = std::complex<T>(1.0,0.0);
2804 
2805  int k = 1;
2806  ctensorT term = B;
2807  while (term.normf() > tol)
2808  {
2809  expB += term;
2810  term = inner(term,B);
2811  k++;
2812  term.scale(1.0/k);
2813  }
2814 
2815  // Repeatedly square to recover exp(A)
2816  while (n--)
2817  {
2818  expB = inner(expB,expB);
2819  }
2820 
2821  return expB;
2822  }
2823  //*************************************************************************
2824 
2825  //*************************************************************************
2826  template<typename Q>
2827  void print_tensor2d(ostream& os, Tensor<Q> t)
2828  {
2829  os.precision(5);
2830  for (int i = 0; i < t.dim(0); i++)
2831  {
2832  for (int j = 0; j < t.dim(1); j++)
2833  {
2834  os << t(i,j) << setw(12);
2835  }
2836  os << endl;
2837  }
2838  os << endl;
2839  }
2840  //*************************************************************************
2841 
2842  //*************************************************************************
2843  void print_potential_matrix_eigs(const vecfuncT& wf, const vecfuncT& vwf)
2844  {
2845  // Build the potential matrix
2847  tensorT potential = matrix_inner(_world, wf, vwf, true);
2848  _world.gop.fence();
2849  END_TIMER(_world,"potential energy matrix");
2850  if (_world.rank()==0) printf("\n");
2851  tensorT overlap = matrix_inner(_world,wf,wf,true);
2852  _world.gop.fence();
2853 
2854  // diagonalize potential
2855  ctensorT cp; rtensorT ep;
2856  sygv(potential, overlap, 1, cp, ep);
2857  if (_world.rank() == 0)
2858  {
2859  print("potential eigenvectors dims:",cp.dim(0),cp.dim(1));
2860  print("potential eigenvectors:");
2861  print(cp);
2862  printf("\n\n");
2863  print("potential eigenvalues:");
2864  print(ep);
2865  printf("\n\n");
2866  }
2867 
2868  }
2869  //*************************************************************************
2870 
2871  //*************************************************************************
2872  void print_fock_matrix_eigs(const vecfuncT& wf, const vecfuncT& vwf, KPoint kpoint)
2873  {
2874  // Build the potential matrix
2876  tensorT potential = matrix_inner(_world, wf, vwf, true);
2877  _world.gop.fence();
2878  END_TIMER(_world,"potential energy matrix");
2879  if (_world.rank()==0) printf("\n");
2880 
2882  if (_world.rank() == 0) _outputF << "Building kinetic energy matrix ...\n\n" << endl;
2883  //tensorT kinetic = ::kinetic_energy_matrix_slow<T,NDIM>(_world, psi,
2884  // _params.periodic,
2885  // kpoint);
2886  tensorT kinetic = ::kinetic_energy_matrix<T,NDIM>(_world, wf,
2887  _params.periodic,
2888  kpoint);
2889  _world.gop.fence();
2890  END_TIMER(_world,"kinetic energy matrix");
2891  if (_world.rank() == 0) printf("\n");
2892 
2893  if (_world.rank() == 0) _outputF << "Constructing Fock matrix ...\n\n" << endl;
2894  tensorT fock = potential + kinetic;
2895  fock = 0.5 * (fock + conj_transpose(fock));
2896  _world.gop.fence();
2897 
2898  // DEBUG
2899  tensorT overlap = matrix_inner(_world,wf,wf,true);
2900  _world.gop.fence();
2901 
2902  // diagonlize kinetic
2903  ctensorT ck; rtensorT ek;
2904  sygv(kinetic, overlap, 1, ck, ek);
2905  // diagonalize potential
2906  ctensorT cp; rtensorT ep;
2907  sygv(potential, overlap, 1, cp, ep);
2908  // diagonalize overlap
2909  ctensorT co; rtensorT eo;
2910  syev(overlap, co, eo);
2911  ctensorT c; rtensorT e;
2912  sygv(fock, overlap, 1, c, e);
2913 
2914  if (_world.rank() == 0)
2915  {
2916  print("kinetic matrix:");
2917  print(kinetic);
2918  print("\nkinetic eigenvalues:");
2919  print(ek);
2920  print("\n");
2921 
2922  print("potential matrix:");
2923  print(potential);
2924  print("\npotential eigenvalues:");
2925  print(ep);
2926  print("\n");
2927 
2928  print("fock matrix:");
2929  print(fock);
2930  print("\nfock eigenvalues:");
2931  print(e);
2932  print("\n");
2933  }
2934 
2935  }
2936  //*************************************************************************
2937 
2938  //*************************************************************************
2939  void do_rhs(vecfuncT& wf,
2940  vecfuncT& vwf,
2941  std::vector<KPoint> kpoints,
2942  std::vector<T>& alpha,
2943  std::vector<double>& eigs)
2944  {
2945  // tolerance
2946  double trantol = 0.1*_params.thresh/std::min(30.0,double(wf.size()));
2947  double thresh = 1e-4;
2948 
2949  if (_world.rank() == 0) _eigF << "Iteration: " << _it << endl;
2950  for (unsigned int kp = 0; kp < kpoints.size(); kp++)
2951  {
2952  // Get k-point and orbitals for this k-point
2953  KPoint kpoint = kpoints[kp];
2954  double k0 = kpoint.k[0];
2955  double k1 = kpoint.k[1];
2956  double k2 = kpoint.k[2];
2957  // Extract the relevant portion of the list of orbitals and the list of the
2958  // V times the orbitals
2959  vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
2960  vecfuncT k_vwf(vwf.begin() + kpoint.begin, vwf.begin() + kpoint.end);
2961  // Build fock matrix
2962  tensorT fock = build_fock_matrix(k_wf, k_vwf, kpoint);
2963  tensorT overlap = matrix_inner(_world, k_wf, k_wf, true);
2964 
2965 // // Do right hand side stuff for kpoint
2966 // bool isgamma = (is_equal(k0,0.0,1e-5) &&
2967 // is_equal(k1,0.0,1e-5) &&
2968 // is_equal(k2,0.0,1e-5));
2969 // if (_params.periodic && !isgamma) // Non-zero k-point
2970 // {
2971 // // Do the gradient term and k^2/2
2972 // vecfuncT d_wf = zero_functions<valueT,NDIM>(_world, k_wf.size());
2973 // complex_derivative_3d Dx(_world,0);
2974 // complex_derivative_3d Dy(_world,1);
2975 // complex_derivative_3d Dz(_world,2);
2976 // for (unsigned int i = 0; i < k_wf.size(); i++)
2977 // {
2978 // // gradient
2979 // functionT dx_wf = Dx(k_wf[i]);
2980 // functionT dy_wf = Dy(k_wf[i]);
2981 // functionT dz_wf = Dz(k_wf[i]);
2982 // d_wf[i] = std::complex<T>(0.0,k0)*dx_wf +
2983 // std::complex<T>(0.0,k1)*dy_wf +
2984 // std::complex<T>(0.0,k2)*dz_wf;
2985 // // k^/2
2986 // double ksq = k0*k0 + k1*k1 + k2*k2;
2987 // k_vwf[i] += 0.5 * ksq * k_wf[i];
2988 // k_vwf[i] -= d_wf[i];
2989 // }
2990 // }
2991 
2992  if (_params.canon) // canonical orbitals
2993  {
2994  ctensorT U; rtensorT e;
2995  sygv(fock, overlap, 1, U, e);
2996 
2997  unsigned int nmo = k_wf.size();
2998  // Fix phases.
2999  long imax;
3000  for (long j = 0; j < nmo; j++)
3001  {
3002  // Get index of largest value in column
3003  U(_,j).absmax(&imax);
3004  T ang = arg(U(imax,j));
3005  std::complex<T> phase = exp(std::complex<T>(0.0,-ang));
3006  // Loop through the rest of the column and divide by the phase
3007  for (long i = 0; i < nmo; i++)
3008  {
3009  U(i,j) *= phase;
3010  }
3011  }
3012 
3013  // Within blocks with the same occupation number attempt to
3014  // keep orbitals in the same order (to avoid confusing the
3015  // non-linear solver). Have to run the reordering multiple
3016  // times to handle multiple degeneracies.
3017  int maxpass = 5;
3018  for (int pass = 0; pass < maxpass; pass++)
3019  {
3020  long j;
3021  for (long i = 0; i < nmo; i++)
3022  {
3023  U(_, i).absmax(&j);
3024  if (i != j)
3025  {
3026  tensorT tmp = copy(U(_, i));
3027  U(_, i) = U(_, j);
3028  U(_, j) = tmp;
3029  //swap(e[i], e[j]);
3030  T ti = e[i];
3031  T tj = e[j];
3032  e[i] = tj; e[j] = ti;
3033  }
3034  }
3035  }
3036 
3037  // Rotations between effectively degenerate states confound
3038  // the non-linear equation solver ... undo these rotations
3039  long ilo = 0; // first element of cluster
3040  while (ilo < nmo-1) {
3041  long ihi = ilo;
3042  while (fabs(real(e[ilo]-e[ihi+1])) < thresh*10.0*max(fabs(real(e[ilo])),1.0)) {
3043  ihi++;
3044  if (ihi == nmo-1) break;
3045  }
3046  long nclus = ihi - ilo + 1;
3047  if (nclus > 1) {
3048  if (_world.rank() == 0) print(" found cluster", ilo, ihi);
3049  tensorT q = copy(U(Slice(ilo,ihi),Slice(ilo,ihi)));
3050  //print(q);
3051  // Special code just for nclus=2
3052  // double c = 0.5*(q(0,0) + q(1,1));
3053  // double s = 0.5*(q(0,1) - q(1,0));
3054  // double r = sqrt(c*c + s*s);
3055  // c /= r;
3056  // s /= r;
3057  // q(0,0) = q(1,1) = c;
3058  // q(0,1) = -s;
3059  // q(1,0) = s;
3060 
3061  // Iteratively construct unitary rotation by
3062  // exponentiating the antisymmetric part of the matrix
3063  // ... is quadratically convergent so just do 3
3064  // iterations
3066  q = inner(q,rot);
3067  ctensorT rot2 = matrix_exponential(-0.5*(q - conj_transpose(q)));
3068  q = inner(q,rot2);
3069  ctensorT rot3 = matrix_exponential(-0.5*(q - conj_transpose(q)));
3070  q = inner(rot,inner(rot2,rot3));
3071  U(_,Slice(ilo,ihi)) = inner(U(_,Slice(ilo,ihi)),q);
3072  }
3073  ilo = ihi+1;
3074  }
3075 
3076  // Debug output
3077  if (_params.print_matrices && _world.rank() == 0)
3078  {
3079  printf("(%10.5f, %10.5f, %10.5f)\n", k0, k1, k2);
3080  print("Overlap matrix:");
3081  print(overlap);
3082 
3083  print("Fock matrix:");
3084  print(fock);
3085 
3086  print("U matrix: (eigenvectors)");
3087  print(U);
3088 
3089  print("Fock matrix eigenvalues:");
3090  print(e);
3091  }
3092 
3093  // WSTHORNTON
3094  //print_fock_matrix_eigs(k_wf, k_vwf, kpoint);
3095 
3096  // transform orbitals and V * (orbitals)
3097  //k_vwf = transform(_world, k_vwf, U, 1e-5 / std::min(30.0, double(k_wf.size())), false);
3098  //k_wf = transform(_world, k_wf, U, FunctionDefaults<3>::get_thresh() / std::min(30.0, double(k_wf.size())), true);
3099 
3100  // WSTHORNTON
3101  //print_fock_matrix_eigs(k_wf, k_vwf, kpoint);
3102 
3103  // Do right hand side stuff for kpoint
3104  bool isgamma = (is_equal(k0,0.0,1e-5) &&
3105  is_equal(k1,0.0,1e-5) &&
3106  is_equal(k2,0.0,1e-5));
3107  if (_params.periodic && !isgamma) // Non-zero k-point
3108  {
3109  // Do the gradient term and k^2/2
3110  vecfuncT d_wf = zero_functions_compressed<valueT,NDIM>(_world, k_wf.size());
3114  for (unsigned int i = 0; i < k_wf.size(); i++)
3115  {
3116  // gradient
3117  functionT dx_wf = Dx(k_wf[i]);
3118  functionT dy_wf = Dy(k_wf[i]);
3119  functionT dz_wf = Dz(k_wf[i]);
3120  d_wf[i] = std::complex<T>(0.0,k0)*dx_wf +
3121  std::complex<T>(0.0,k1)*dy_wf +
3122  std::complex<T>(0.0,k2)*dz_wf;
3123  // k^/2
3124  double ksq = k0*k0 + k1*k1 + k2*k2;
3125  k_vwf[i] += 0.5 * ksq * k_wf[i];
3126  k_vwf[i] -= d_wf[i];
3127  }
3128  }
3129 
3130  // WSTHORNTON (new code)
3131  unsigned int eimax = -1;
3132  double eimaxval = -1e10;
3133  for (unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3134  ei++, fi++)
3135  {
3136  if ((real(e(fi,fi)) > 0.0) && (real(e(fi,fi)) > eimaxval))
3137  {
3138  eimax = fi;
3139  eimaxval = real(e(fi,fi));
3140  }
3141  }
3142 
3143  double eshift = (eimaxval > 0.0) ? eimaxval + 0.1 : 0.0;
3144  for (unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3145  ei++, fi++)
3146  {
3147  // Save the latest eigenvalues
3148  eigs[ei] = real(e(fi,fi));
3149  alpha[ei] = e(fi,fi)-eshift;
3150  k_vwf[fi] += (alpha[ei]-eigs[ei])*k_wf[fi];
3151  }
3152 
3153  if (_world.rank() == 0)
3154  {
3155  _eigF << "kpt: " << kp << endl;
3156  _eigF << setfill('-') << setw(20) << " " << endl;
3157  for (unsigned int ei = kpoint.begin; ei < kpoint.end; ei++)
3158  {
3159  char eigstr[50];
3160  printf("%3d%15.10f",ei,real(eigs[ei]));
3161 // _eigF << ei << setfill(' ') << setw(12) << real(eigs[ei]) << endl;
3162  _eigF << eigstr << endl;
3163  }
3164  _eigF << "\n\n" << endl;
3165  }
3166  }
3167  else // non-canonical orbitals
3168  {
3169  // diagonlize just to print eigenvalues
3170  tensorT overlap = matrix_inner(_world, k_wf, k_wf, true);
3171  ctensorT c; rtensorT e;
3172  sygv(fock, overlap, 1, c, e);
3173  for (unsigned int ei = 0; ei < e.dim(0); ei++)
3174  {
3175  double diffe = (ei == 0) ? 0.0 : real(e(ei,ei))-real(e(ei-1,ei-1));
3176  if (_world.rank() == 0)
3177  print("kpoint ", kp, "ei ", ei, "eps ", real(e(ei,ei)), "\tdiff\t", diffe);
3178  }
3179 
3180  for (unsigned int ei = kpoint.begin, fi = 0;
3181  ei < kpoint.end; ei++, fi++)
3182  {
3183  alpha[ei] = std::min(-0.1, real(fock(fi,fi)));
3184  fock(fi,fi) -= std::complex<T>(alpha[ei], 0.0);
3185  }
3186 
3187  std::vector<functionT> fwf = transform(_world, k_wf, fock, trantol);
3188  gaxpy(_world, 1.0, k_vwf, -1.0, fwf);
3189  fwf.clear();
3190  }
3191  for (unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3192  wi++, fi++)
3193  {
3194  wf[wi] = k_wf[fi];
3195  vwf[wi] = k_vwf[fi];
3196  }
3197  }
3198  }
3199  //*************************************************************************
3200 
3201  //*************************************************************************
3203  vecfuncT& vwf,
3204  std::vector<KPoint> kpoints,
3205  std::vector<T>& alpha,
3206  std::vector<double>& eigs)
3207  {
3208  // tolerance
3209  double trantol = 0.1*_params.thresh/std::min(30.0,double(wf.size()));
3210  double thresh = 1e-4;
3211 
3212  if (_world.rank() == 0) _eigF << "Iteration: " << _it << endl;
3213  for (unsigned int kp = 0; kp < kpoints.size(); kp++)
3214  {
3215  // Get k-point and orbitals for this k-point
3216  KPoint kpoint = kpoints[kp];
3217  double k0 = kpoint.k[0];
3218  double k1 = kpoint.k[1];
3219  double k2 = kpoint.k[2];
3220  // Extract the relevant portion of the list of orbitals and the list of the
3221  // V times the orbitals
3222  vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
3223  vecfuncT k_vwf(vwf.begin() + kpoint.begin, vwf.begin() + kpoint.end);
3224 
3225  // Build fock matrix
3226  tensorT fock = build_fock_matrix(k_wf, k_vwf, kpoint);
3227  tensorT overlap = matrix_inner(_world, k_wf, k_wf, true);
3228  for (unsigned int i = 0; i < k_wf.size(); i++)
3229  fock(i,i) += std::complex<double>(i*_params.thresh*0.1,0.0);
3230 
3231  ctensorT U; rtensorT e;
3232  sygv(fock, overlap, 1, U, e);
3233 
3234  if (_params.print_matrices && _world.rank() == 0)
3235  {
3236  printf("(%10.5f, %10.5f, %10.5f)\n", k0, k1, k2);
3237  print("Overlap matrix:");
3238  print(overlap);
3239 
3240  print("Fock matrix:");
3241  print(fock);
3242 
3243  print("U matrix: (eigenvectors)");
3244  print(U);
3245 
3246  print("Fock matrix eigenvalues:");
3247  print(e);
3248  }
3249 
3250 
3251  // this is all of the B.S. for the solver
3252  if (_params.solver == 1)
3253  {
3254  unsigned int nmo = k_wf.size();
3255  // Fix phases.
3256  long imax;
3257  for (long j = 0; j < nmo; j++)
3258  {
3259  // Get index of largest value in column
3260  U(_,j).absmax(&imax);
3261  T ang = arg(U(imax,j));
3262  std::complex<T> phase = exp(std::complex<T>(0.0,-ang));
3263  // Loop through the rest of the column and divide by the phase
3264  for (long i = 0; i < nmo; i++)
3265  {
3266  U(i,j) *= phase;
3267  }
3268  }
3269 
3270  // Within blocks with the same occupation number attempt to
3271  // keep orbitals in the same order (to avoid confusing the
3272  // non-linear solver). Have to run the reordering multiple
3273  // times to handle multiple degeneracies.
3274  int maxpass = 5;
3275  for (int pass = 0; pass < maxpass; pass++)
3276  {
3277  long j;
3278  for (long i = 0; i < nmo; i++)
3279  {
3280  U(_, i).absmax(&j);
3281  if (i != j)
3282  {
3283  tensorT tmp = copy(U(_, i));
3284  U(_, i) = U(_, j);
3285  U(_, j) = tmp;
3286  //swap(e[i], e[j]);
3287  T ti = e[i];
3288  T tj = e[j];
3289  e[i] = tj; e[j] = ti;
3290  }
3291  }
3292  }
3293 
3294  // Rotations between effectively degenerate states confound
3295  // the non-linear equation solver ... undo these rotations
3296  long ilo = 0; // first element of cluster
3297  while (ilo < nmo-1) {
3298  long ihi = ilo;
3299  while (fabs(real(e[ilo]-e[ihi+1])) < thresh*10.0*max(fabs(real(e[ilo])),1.0)) {
3300  ihi++;
3301  if (ihi == nmo-1) break;
3302  }
3303  long nclus = ihi - ilo + 1;
3304  if (nclus > 1) {
3305  if (_world.rank() == 0) print(" found cluster", ilo, ihi);
3306  tensorT q = copy(U(Slice(ilo,ihi),Slice(ilo,ihi)));
3307  //print(q);
3308  // Special code just for nclus=2
3309  // double c = 0.5*(q(0,0) + q(1,1));
3310  // double s = 0.5*(q(0,1) - q(1,0));
3311  // double r = sqrt(c*c + s*s);
3312  // c /= r;
3313  // s /= r;
3314  // q(0,0) = q(1,1) = c;
3315  // q(0,1) = -s;
3316  // q(1,0) = s;
3317 
3318  // Iteratively construct unitary rotation by
3319  // exponentiating the antisymmetric part of the matrix
3320  // ... is quadratically convergent so just do 3
3321  // iterations
3323  q = inner(q,rot);
3324  ctensorT rot2 = matrix_exponential(-0.5*(q - conj_transpose(q)));
3325  q = inner(q,rot2);
3326  ctensorT rot3 = matrix_exponential(-0.5*(q - conj_transpose(q)));
3327  q = inner(rot,inner(rot2,rot3));
3328  U(_,Slice(ilo,ihi)) = inner(U(_,Slice(ilo,ihi)),q);
3329  }
3330  ilo = ihi+1;
3331  }
3332  }
3333 
3334  // transform orbitals and V * (orbitals)
3335  k_vwf = transform(_world, k_vwf, U, 1e-5 / std::min(30.0, double(k_wf.size())), false);
3336  k_wf = transform(_world, k_wf, U, FunctionDefaults<3>::get_thresh() / std::min(30.0, double(k_wf.size())), true);
3337 
3338  // Do right hand side stuff for kpoint
3339  bool isgamma = (is_equal(k0,0.0,1e-5) &&
3340  is_equal(k1,0.0,1e-5) &&
3341  is_equal(k2,0.0,1e-5));
3342  if (_params.periodic && !isgamma) // Non-zero k-point
3343  {
3344  // Do the gradient term and k^2/2
3345  vecfuncT d_wf = zero_functions_compressed<valueT,NDIM>(_world, k_wf.size());
3349  for (unsigned int i = 0; i < k_wf.size(); i++)
3350  {
3351  // gradient
3352  functionT dx_wf = Dx(k_wf[i]);
3353  functionT dy_wf = Dy(k_wf[i]);
3354  functionT dz_wf = Dz(k_wf[i]);
3355 
3356  // WSTHORNTON
3357 // double delxx = std::abs(inner(k_wf[i],Dx(Dx(k_wf[i]))));
3358 // double delyy = std::abs(inner(k_wf[i],Dy(Dy(k_wf[i]))));
3359 // double delzz = std::abs(inner(k_wf[i],Dz(Dz(k_wf[i]))));
3360 
3361 // if (_world.rank() == 0)
3362 // printf("orb %2.2d delxx: %15.8e delyy: %15.8e delzz: %15.8e\n",i, delxx, delyy, delzz);
3363 
3364  // WSTHORNTON
3365 // std::vector<long> npt(3,101);
3366 // char fnamedx[50];
3367 // sprintf(fnamedx, "xorb%2.2d__.dx",i);
3368 // plotdx(k_wf[i], fnamedx, FunctionDefaults<3>::get_cell(), npt);
3369 //
3370 // char fnamedx_dx[50];
3371 // sprintf(fnamedx_dx, "xorb%2.2d__dx.dx",i);
3372 // plotdx(dx_wf, fnamedx_dx, FunctionDefaults<3>::get_cell(), npt);
3373 // char fnamedx_dy[50];
3374 // sprintf(fnamedx_dy, "xorb%2.2d__dy.dx",i);
3375 // plotdx(dy_wf, fnamedx_dy, FunctionDefaults<3>::get_cell(), npt);
3376 // char fnamedx_dz[50];
3377 // sprintf(fnamedx_dz, "xorb%2.2d__dz.dx",i);
3378 // plotdx(dz_wf, fnamedx_dz, FunctionDefaults<3>::get_cell(), npt);
3379 //
3380 // rfunctionT xfunc = rfactoryT(_world).functor(
3381 // rfunctorT(new DipoleFunctor(0))).
3382 // thresh(_params.thresh).truncate_on_project();
3383 // rfunctionT yfunc = rfactoryT(_world).functor(
3384 // rfunctorT(new DipoleFunctor(1))).
3385 // thresh(_params.thresh).truncate_on_project();
3386 // rfunctionT zfunc = rfactoryT(_world).functor(
3387 // rfunctorT(new DipoleFunctor(2))).
3388 // thresh(_params.thresh).truncate_on_project();
3389 //
3390 // double xdip = std::abs(inner(k_wf[i],xfunc*k_wf[i]));
3391 // double ydip = std::abs(inner(k_wf[i],yfunc*k_wf[i]));
3392 // double zdip = std::abs(inner(k_wf[i],zfunc*k_wf[i]));
3393 //
3394 // if (_world.rank() == 0)
3395 // printf("orb: %1d dipole moment: (%7.5e,%7.5e,%7.5e)\n",i,xdip,ydip,zdip);
3396 //
3397 // double xdip1 = std::real(inner(k_wf[i],xfunc*k_wf[i]));
3398 // double ydip1 = std::real(inner(k_wf[i],yfunc*k_wf[i]));
3399 // double zdip1 = std::real(inner(k_wf[i],zfunc*k_wf[i]));
3400 //
3401 // if (_world.rank() == 0)
3402 // printf("orb: %1d dipole moment: (%7.5e,%7.5e,%7.5e)\n",i,xdip1,ydip1,zdip1);
3403 //
3404 // double xdip2 = std::imag(inner(k_wf[i],xfunc*k_wf[i]));
3405 // double ydip2 = std::imag(inner(k_wf[i],yfunc*k_wf[i]));
3406 // double zdip2 = std::imag(inner(k_wf[i],zfunc*k_wf[i]));
3407 //
3408 // if (_world.rank() == 0)
3409 // printf("orb: %1d dipole moment: (%7.5e,%7.5e,%7.5e)\n\n\n",i,xdip2,ydip2,zdip2);
3410 //
3411 // // WSTHORNTON
3412 // char fname_x[50];
3413 // char fname_y[50];
3414 // char fname_z[50];
3415 //
3416 // char fnamex_x[50];
3417 // char fnamex_y[50];
3418 // char fnamex_z[50];
3419 //
3420 // char fnamey_x[50];
3421 // char fnamey_y[50];
3422 // char fnamey_z[50];
3423 //
3424 // char fnamez_x[50];
3425 // char fnamez_y[50];
3426 // char fnamez_z[50];
3427 //
3428 // sprintf(fname_x, "orb%2.2d__x.dat",i);
3429 // sprintf(fname_y, "orb%2.2d__y.dat",i);
3430 // sprintf(fname_z, "orb%2.2d__z.dat",i);
3431 //
3432 // sprintf(fnamex_x, "orb%2.2d_dx_x.dat",i);
3433 // sprintf(fnamex_y, "orb%2.2d_dx_y.dat",i);
3434 // sprintf(fnamex_z, "orb%2.2d_dx_z.dat",i);
3435 //
3436 // sprintf(fnamey_x, "orb%2.2d_dy_x.dat",i);
3437 // sprintf(fnamey_y, "orb%2.2d_dy_y.dat",i);
3438 // sprintf(fnamey_z, "orb%2.2d_dy_z.dat",i);
3439 //
3440 // sprintf(fnamez_x, "orb%2.2d_dz_x.dat",i);
3441 // sprintf(fnamez_y, "orb%2.2d_dz_y.dat",i);
3442 // sprintf(fnamez_z, "orb%2.2d_dz_z.dat",i);
3443 //
3444 // coord_3d pt1x {-_params.L/2,0.0,0.0};
3445 // coord_3d pt2x { _params.L/2,0.0,0.0};
3446 //
3447 // coord_3d pt1y {0.0,-_params.L/2,0.0};
3448 // coord_3d pt2y {0.0, _params.L/2,0.0};
3449 //
3450 // coord_3d pt1z {0.0,0.0,-_params.L/2};
3451 // coord_3d pt2z {0.0,0.0, _params.L/2};
3452 //
3453 // plot_line(fname_x,30000,pt1x,pt2x,k_wf[i]);
3454 // plot_line(fname_y,30000,pt1y,pt2y,k_wf[i]);
3455 // plot_line(fname_z,30000,pt1z,pt2z,k_wf[i]);
3456 //
3457 // plot_line(fnamex_x,30000,pt1x,pt2x,dx_wf);
3458 // plot_line(fnamex_y,30000,pt1y,pt2y,dx_wf);
3459 // plot_line(fnamex_z,30000,pt1z,pt2z,dx_wf);
3460 //
3461 // plot_line(fnamey_x,30000,pt1x,pt2x,dy_wf);
3462 // plot_line(fnamey_y,30000,pt1y,pt2y,dy_wf);
3463 // plot_line(fnamey_z,30000,pt1z,pt2z,dy_wf);
3464 //
3465 // plot_line(fnamez_x,30000,pt1x,pt2x,dz_wf);
3466 // plot_line(fnamez_y,30000,pt1y,pt2y,dz_wf);
3467 // plot_line(fnamez_z,30000,pt1z,pt2z,dz_wf);
3468 
3469 
3470 
3471  d_wf[i] = std::complex<T>(0.0,k0)*dx_wf +
3472  std::complex<T>(0.0,k1)*dy_wf +
3473  std::complex<T>(0.0,k2)*dz_wf;
3474  // k^/2
3475  double ksq = k0*k0 + k1*k1 + k2*k2;
3476  k_vwf[i] += 0.5 * ksq * k_wf[i];
3477  k_vwf[i] -= d_wf[i];
3478  }
3479  }
3480 
3481 
3482  // WSTHORNTON (new code)
3483  unsigned int eimax = -1;
3484  double eimaxval = -1e10;
3485  for (unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3486  ei++, fi++)
3487  {
3488  if ((real(e(fi,fi)) > 0.0) && (real(e(fi,fi)) > eimaxval))
3489  {
3490  eimax = fi;
3491  eimaxval = real(e(fi,fi));
3492  }
3493  }
3494 
3495  double eshift = (eimaxval > 0.0) ? eimaxval + 0.1 : 0.0;
3496  for (unsigned int ei = kpoint.begin, fi = 0; ei < kpoint.end;
3497  ei++, fi++)
3498  {
3499  // Save the latest eigenvalues
3500  eigs[ei] = real(e(fi,fi));
3501  alpha[ei] = e(fi,fi)-eshift;
3502  k_vwf[fi] += (alpha[ei]-eigs[ei])*k_wf[fi];
3503  }
3504 
3505 // if (_world.rank() == 0) printf("do_rhs_simple: (3) ...\n\n");
3506 // print_fock_matrix_eigs(k_wf, k_vwf, kpoint_gamma);
3507 
3508  if (_world.rank() == 0)
3509  {
3510  _eigF << "kpt: " << kp << endl;
3511  _eigF << setfill('-') << setw(20) << " " << endl;
3512  for (unsigned int ei = kpoint.begin; ei < kpoint.end; ei++)
3513  {
3514  char eigstr[50];
3515  sprintf(eigstr,"%3d%15.10f",ei,real(eigs[ei]));
3516  _eigF << eigstr << endl;
3517  }
3518  _eigF << "\n\n" << endl;
3519  }
3520  for (unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3521  wi++, fi++)
3522  {
3523  wf[wi] = k_wf[fi];
3524  vwf[wi] = k_vwf[fi];
3525  }
3526  }
3527  }
3528  //*************************************************************************
3529 
3530  //*************************************************************************
3532  const double_complex I = double_complex(0.0,1.0);
3533  double kx = k.k[0];
3534  double ky = k.k[1];
3535  double kz = k.k[2];
3536 
3537  complex_derivative_3d Dx(world, 0);
3538  complex_derivative_3d Dy(world, 1);
3539  complex_derivative_3d Dz(world, 2);
3540 
3541  vector_complex_function_3d dvx = apply(world, Dx, v);
3542  vector_complex_function_3d dvy = apply(world, Dy, v);
3543  vector_complex_function_3d dvz = apply(world, Dz, v);
3544 
3545  // -1/2 (del + ik)^2 = -1/2 del^2 - i k.del + 1/2 k^2
3546  // -1/2 <p|del^2|q> = +1/2 <del p | del q>
3547 
3548  tensor_complex f1 = 0.5 * (matrix_inner(world, dvx, dvx, true) +
3549  matrix_inner(world, dvy, dvy, true) +
3550  matrix_inner(world, dvz, dvz, true));
3551 
3552  tensor_complex f2 =
3553  (-I*kx)*matrix_inner(world, v, dvx, false) +
3554  (-I*ky)*matrix_inner(world, v, dvy, false) +
3555  (-I*kz)*matrix_inner(world, v, dvz, false);
3556 
3557  tensor_complex f3 = (0.5 * (kx*kx + ky*ky + kz*kz)) * matrix_inner(world, v, v, true);
3558 
3559  return f1 + f2 + f3;
3560  }
3561  //*************************************************************************
3562 
3563  //*************************************************************************
3565  vecfuncT& vpsi,
3566  KPoint kpoint)
3567  {
3568  // Build the potential matrix
3570  tensorT potential = matrix_inner(_world, psi, vpsi, true);
3571  _world.gop.fence();
3572  END_TIMER(_world,"potential energy matrix");
3573  if (_world.rank()==0) printf("\n");
3574 
3576  if (_world.rank() == 0) _outputF << "Building kinetic energy matrix ...\n\n" << endl;
3577  tensorT kinetic = ::kinetic_energy_matrix<T,NDIM>(_world, psi,
3578  _params.periodic,
3579  kpoint);
3580  tensorT kinetic2 = make_kinetic_matrix(_world, psi, kpoint);
3581  _world.gop.fence();
3582  END_TIMER(_world,"kinetic energy matrix");
3583  if (_world.rank() == 0) printf("\n");
3584 
3585  if (_world.rank() == 0) _outputF << "Constructing Fock matrix ...\n\n" << endl;
3586  tensorT fock = potential + kinetic;
3587  fock = 0.5 * (fock + conj_transpose(fock));
3588  _world.gop.fence();
3589 
3590  if (_world.rank() == 0)
3591  {
3592  print("KINETIC:");
3593  print(kinetic);
3594  print("KINETIC2:");
3595  print(kinetic2);
3596  print("POTENTIAL:");
3597  print(potential);
3598  print("FOCK:");
3599  print(fock);
3600  }
3601 
3602  return fock;
3603  }
3604  //*************************************************************************
3605 
3606  //*************************************************************************
3608  {
3609  for (unsigned int fi = kpoint.begin; fi < kpoint.end; ++fi)
3610  {
3611  // Project out the lower states
3612  for (unsigned int fj = kpoint.begin; fj < fi; ++fj)
3613  {
3614  valueT overlap = inner(f[fj], f[fi]);
3615  f[fi] -= overlap*f[fj];
3616  }
3617  f[fi].scale(1.0/f[fi].norm2());
3618  }
3619  }
3620  //*************************************************************************
3621 
3622  //*************************************************************************
3623  /// Given overlap matrix, return rotation with 3rd order error to orthonormalize the vectors
3624  tensorT Q3(const tensorT& s) {
3625  tensorT Q = inner(s,s);
3626  Q.gaxpy(0.2,s,-2.0/3.0);
3627  for (int i=0; i<s.dim(0); ++i) Q(i,i) += 1.0;
3628  return Q.scale(15.0/8.0);
3629  }
3630  //*************************************************************************
3631 
3632  //*************************************************************************
3633  /// Computes matrix square root (not used any more?)
3634  ctensorT csqrt(const ctensorT& s, double tol=1e-8) {
3635  int n=s.dim(0), m=s.dim(1);
3636  MADNESS_ASSERT(n==m);
3637  ctensorT c; rtensorT e;
3638  //s.gaxpy(0.5,conj_transpose(s),0.5); // Ensure exact symmetry
3639  syev(s, c, e);
3640  for (int i=0; i<n; ++i) {
3641  if (e(i) < -tol) {
3642  MADNESS_EXCEPTION("Matrix square root: negative eigenvalue",i);
3643  }
3644  else if (e(i) < tol) { // Ugh ..
3645  print("Matrix square root: Warning: small eigenvalue ", i, e(i));
3646  e(i) = tol;
3647  }
3648  e(i) = 1.0/sqrt(e(i));
3649  }
3650  for (int j=0; j<n; ++j) {
3651  for (int i=0; i<n; ++i) {
3652  c(j,i) *= e(i);
3653  }
3654  }
3655  return c;
3656  }
3657  //*************************************************************************
3658 
3659  //*************************************************************************
3660  void orthonormalize(vecfuncT& wf, KPoint kpoint)
3661  {
3662  // extract k-point orbitals
3663  vecfuncT k_wf(wf.begin() + kpoint.begin, wf.begin() + kpoint.end);
3664  ctensorT S = matrix_inner(_world,k_wf,k_wf,true);
3665  printf("orthonormalize: \n");
3666  printf("before matrix (after): \n");
3667  print(S);
3668  ctensorT U = csqrt(S);
3669  k_wf = transform(_world, k_wf, U, _params.thresh, true);
3670 
3671  ctensorT S2 = matrix_inner(_world,k_wf,k_wf,true);
3672  printf("overlap matrix (after): \n");
3673  print(S2);
3674  for (unsigned int wi = kpoint.begin, fi = 0; wi < kpoint.end;
3675  wi++, fi++)
3676  {
3677  wf[wi] = k_wf[fi];
3678  }
3679  }
3680  //*************************************************************************
3681 
3682  //*************************************************************************
3684  const vecfuncT& bwfs)
3685  {
3686  // vector of residual functions
3687  vecfuncT rm = sub(_world, _phisa, awfs);
3688  // if spin-polarized insert beta spin-orbital residual functions
3689  if (_params.spinpol)
3690  {
3691  vecfuncT br = sub(_world, _phisb, bwfs);
3692  rm.insert(rm.end(), br.begin(), br.end());
3693  }
3694  // scalar residual
3695  std::vector<double> rnvec = norm2s<valueT,NDIM>(_world, rm);
3696  double rnorm = 0.0;
3697  for (unsigned int i = 0; i < rnvec.size(); i++) rnorm += rnvec[i];
3698  // renormalize and print
3699  _residual = rnorm / rnvec.size();
3700  if (_world.rank() == 0) _outputF << "\nResiduals\n---------" << endl;
3701  if (_world.rank() == 0) _outputF << std::setiosflags(std::ios::scientific) << "residual = " << _residual << endl;
3702  if (_world.rank() == 0)
3703  {
3704  _outputF << endl;
3705  for (unsigned int i = 0; i < rnvec.size(); i++)
3706  {
3707  _outputF << "residual" << i << "\t" << rnvec[i] << endl;
3708  }
3709  _outputF << endl;
3710  }
3711 
3712  return rm;
3713  }
3714 
3715  //*************************************************************************
3717  vecfuncT& bwfs,
3718  std::vector<KPoint> kpoints)
3719  {
3720  // truncate before we do anyting
3721  truncate<valueT,NDIM> (_world, awfs);
3722  truncate<valueT,NDIM> (_world, _phisa);
3723  if (_params.spinpol)
3724  {
3725  truncate<valueT,NDIM> (_world, bwfs);
3726  truncate<valueT,NDIM> (_world, _phisb);
3727  }
3728  // compute residual
3729  vecfuncT rm = compute_residual(awfs, bwfs);
3730  if (_params.solver > 0 && _params.maxsub > 1)
3731  {
3732  // nonlinear solver
3733  _subspace->update_subspace(awfs, bwfs, _phisa, _phisb, rm);
3734  }
3735 
3736  // do step restriction
3737  step_restriction(_phisa, awfs, 0);
3738  if (_params.spinpol)
3739  {
3740  step_restriction(_phisb, bwfs, 1);
3741  }
3742  // do gram-schmidt
3743  for (unsigned int kp = 0; kp < kpoints.size(); kp++)
3744  {
3745  gram_schmidt(awfs, kpoints[kp]);
3746 // orthonormalize(awfs, kpoints[kp]);
3747  if (_params.spinpol)
3748  {
3749  gram_schmidt(bwfs, kpoints[kp]);
3750 // orthonormalize(bwfs, kpoints[kp]);
3751  }
3752  }
3753 
3754  // update alpha and beta orbitals
3755  truncate<valueT,NDIM>(_world, awfs);
3756  for (unsigned int ai = 0; ai < awfs.size(); ai++) {
3757  _phisa[ai] = awfs[ai].scale(1.0 / awfs[ai].norm2());
3758  }
3759  if (_params.spinpol)
3760  {
3761  truncate<valueT,NDIM>(_world, bwfs);
3762  for (unsigned int bi = 0; bi < bwfs.size(); bi++)
3763  {
3764  _phisb[bi] = bwfs[bi].scale(1.0 / bwfs[bi].norm2());
3765  }
3766  }
3767  else
3768  {
3769  for (unsigned int ia = 0; ia < awfs.size(); ia++)
3770  {
3771  _phisb[ia] = _phisa[ia];
3772  }
3773  }
3774  }
3775  //*************************************************************************
3776 
3777  //*************************************************************************
3779  vecfuncT& nwfs,
3780  int aorb)
3781  {
3782  double s = (_it < 4) ? 0.75 : _params.sd;
3783  if (_world.rank() == 0) print("damping factor: ", s);
3784  for (unsigned int i = 0; i < owfs.size(); i++)
3785  nwfs[i].gaxpy(1.0 - s, owfs[i], s, false);
3786 // std::vector<double> rnorm = norm2s(_world, sub(_world, owfs, nwfs));
3787 // // Step restriction
3788 // int nres = 0;
3789 // for (unsigned int i = 0; i < owfs.size(); i++)
3790 // {
3791 // if (rnorm[i] > _params.maxrotn)
3792 // {
3793 // double s = _params.maxrotn / rnorm[i];
3794 // nres++;
3795 // if (_world.rank() == 0)
3796 // {
3797 // if (!aorb && nres == 1) _outputF << " restricting step for alpha orbitals:" << endl;
3798 // if (aorb && nres == 1) _outputF << " restricting step for beta orbitals:" << endl;
3799 // _outputF << i;
3800 // }
3801 // nwfs[i].gaxpy(s, owfs[i], 1.0 - s, false);
3802 // }
3803 // }
3804 // if (nres > 0 && _world.rank() == 0) printf("\n");
3805 // _world.gop.fence();
3806  }
3807  //*************************************************************************
3808 
3809  //*************************************************************************
3810  void fix_occupations(const std::vector<T>& eps,
3811  std::vector<double>& occs)
3812  {
3813  // Find max/min eigenvalues
3814  double emax = eps[0];
3815  double emin = emax;
3816  for (int i = 0; i < eps.size(); i++)
3817  {
3818  emax = (eps[i] > emax) ? eps[i] : emax;
3819  emin = (eps[i] < emin) ? eps[i] : emin;
3820  }
3821 
3822  int maxits = 1000;
3823  // This is hardcoded to 2.0 (non-spinpolarized case) for now.
3824  double occmax = 2.0;
3825  // Fermi energy
3826  double efermi = 0.0;
3827  // Use bisection method to find the fermi energy and update occupation numbers
3828  bool bstop = false;
3829  // Some smoothing parameter
3830  double t1 = 1.0/_params.swidth;
3831  for (int it = 0; (it < maxits)&&(!bstop); it++)
3832  {
3833  // Proposed fermi energy
3834  efermi = 0.5 * (emax + emin);
3835  // Accumulated charge
3836  double charge = 0.0;
3837  // Loop over all eigenvalues and count the charge
3838  for (int i = 0; i < eps.size(); i++)
3839  {
3840  double x = (efermi-eps[i]) * t1;
3841  // need to add some smearing function here
3842  occs[i] = occmax*stheta_fd(x);
3843  charge += _kpoints[i].weight() * occs[i];
3844  }
3845  if (fabs(emax-emin) < 1e-5)
3846  bstop = true;
3847  else if (charge < _params.ncharge)
3848  emin = efermi;
3849  else
3850  emax = efermi;
3851  }
3852  }
3853  //*************************************************************************
3854 
3855 // //*************************************************************************
3856 // void update_eigenvalues(const vecfuncT& wavefs,
3857 // const vecfuncT& pfuncs, const vecfuncT& phis,
3858 // std::vector<T>& eigs)
3859 // {
3860 // // Update e
3861 // if (_world.rank() == 0) printf("Updating e ...\n\n");
3862 // for (unsigned int ei = 0; ei < eigs.size(); ei++)
3863 // {
3864 // functionT r = wavefs[ei] - phis[ei];
3865 // double tnorm = wavefs[ei].norm2();
3866 // // Compute correction to the eigenvalues
3867 // T ecorrection = -0.5*real(inner(pfuncs[ei], r)) / (tnorm*tnorm);
3868 // T eps_old = eigs[ei];
3869 // T eps_new = eps_old + ecorrection;
3870 //// if (_world.rank() == 0) printf("ecorrection = %.8f\n\n", ecorrection);
3871 //// if (_world.rank() == 0) printf("eps_old = %.8f eps_new = %.8f\n\n", eps_old, eps_new);
3872 // // Sometimes eps_new can go positive, THIS WILL CAUSE THE ALGORITHM TO CRASH. So,
3873 // // I bounce the new eigenvalue back into the negative side of the real axis. I
3874 // // keep doing this until it's good or I've already done it 10 times.
3875 // int counter = 50;
3876 // while (eps_new >= 0.0 && counter < 20)
3877 // {
3878 // // Split the difference between the new and old estimates of the
3879 // // pi-th eigenvalue.
3880 // eps_new = eps_old + 0.5 * (eps_new - eps_old);
3881 // counter++;
3882 // }
3883 // // Still no go, forget about it. (1$ to Donnie Brasco)
3884 // if (eps_new >= 0.0)
3885 // {
3886 // if (_world.rank() == 0) printf("FAILURE OF WST: exiting!!\n\n");
3887 // _exit(0);
3888 // }
3889 // // Set new eigenvalue
3890 // eigs[ei] = eps_new;
3891 // }
3892 // }
3893 // //*************************************************************************
3894 
3895 // //*************************************************************************
3896 // double get_eig(int indx)
3897 // {
3898 // return _solver->get_eig(indx);
3899 // }
3900 // //*************************************************************************
3901 //
3902 // //*************************************************************************
3903 // functionT get_phi(int indx)
3904 // {
3905 // return _solver->get_phi(indx);
3906 // }
3907 // //*************************************************************************
3908 //
3909 // //*************************************************************************
3910 // const std::vector<double>& eigs()
3911 // {
3912 // return _solver->eigs();
3913 // }
3914 // //*************************************************************************
3915 //
3916 // //*************************************************************************
3917 // const vecfuncT& phis()
3918 // {
3919 // return _solver->phis();
3920 // }
3921 // //*************************************************************************
3922 
3923  };
3924  //***************************************************************************
3925 
3926 }
3927 #define SOLVER_H_
3928 
3929 #endif /* SOLVER_H_ */
double potential(const coord_3d &r)
Definition: 3dharmonic.cc:132
double q(double t)
Definition: DKops.h:18
real_convolution_3d A(World &world)
Definition: DKops.h:230
This header should include pretty much everything needed for the parallel runtime.
std::complex< double > double_complex
Definition: cfft.h:14
Definition: test_ar.cc:118
Definition: test_ar.cc:141
Definition: mentity.h:95
void read_file(const std::string &filename, bool fractional)
Definition: mentity.cc:289
int natom() const
Definition: mentity.h:112
double total_nuclear_charge() const
Definition: mentity.cc:444
void center()
Moves the center of nuclear charge to the origin.
Definition: mentity.cc:413
const Atom & get_atom(unsigned int i) const
Definition: mentity.cc:369
Definition: electronicstructureapp.h:85
Definition: molecule.h:58
double y
Definition: molecule.h:60
double x
Definition: molecule.h:60
double z
Definition: molecule.h:60
Used to represent one basis function from a shell on a specific center.
Definition: madness/chem/molecularbasis.h:405
void get_coords(double &x, double &y, double &z) const
Definition: madness/chem/molecularbasis.h:450
Contracted Gaussian basis.
Definition: madness/chem/molecularbasis.h:465
double eval_guess_density(const Molecule &molecule, double x, double y, double z) const
Evaluates the guess density.
Definition: madness/chem/molecularbasis.h:642
AtomicBasisFunction get_atomic_basis_function(const Molecule &molecule, size_t ibf) const
Returns the ibf'th atomic basis function.
Definition: madness/chem/molecularbasis.h:596
void read_file(std::string filename)
read the atomic basis set from file
Definition: molecularbasis.cc:118
int nbf(const Molecule &molecule) const
Given a molecule count the number of basis functions.
Definition: madness/chem/molecularbasis.h:620
long dim(int i) const
Returns the size of dimension i.
Definition: basetensor.h:147
Definition: solver.h:98
const vec3dT exponent
Definition: solver.h:103
Vector< double, NDIM > coordT
Definition: solver.h:100
double_complex operator()(const coordT &x) const
You should implement this to return f(x)
Definition: solver.h:108
ComplexExp(vec3dT exponent, double_complex coeff)
Definition: solver.h:105
Vector< double, NDIM > vec3dT
Definition: solver.h:101
const double_complex coeff
Definition: solver.h:102
Implements derivatives operators with variety of boundary conditions on simulation domain.
Definition: derivative.h:266
A MADNESS functor to compute either x, y, or z.
Definition: solver.h:165
virtual ~DipoleFunctor()
Definition: solver.h:173
DipoleFunctor(int axis)
Definition: solver.h:169
double operator()(const coordT &x) const
Definition: solver.h:170
const int axis
Definition: solver.h:167
FunctionDefaults holds default paramaters as static class members.
Definition: funcdefaults.h:204
static int get_k()
Returns the default wavelet order.
Definition: funcdefaults.h:266
static void set_thresh(double value)
Sets the default threshold.
Definition: funcdefaults.h:286
static void set_k(int value)
Sets the default wavelet order.
Definition: funcdefaults.h:273
static const double & get_thresh()
Returns the default threshold.
Definition: funcdefaults.h:279
static void set_bc(const BoundaryConditions< NDIM > &value)
Sets the default boundary conditions.
Definition: funcdefaults.h:418
static void set_cubic_cell(double lo, double hi)
Sets the user cell to be cubic with each dimension having range [lo,hi].
Definition: funcdefaults.h:461
static const Tensor< double > & get_cell()
Gets the user cell for the simulation.
Definition: funcdefaults.h:446
static const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition: funcdefaults.h:468
FunctionFactory implements the named-parameter idiom for Function.
Definition: function_factory.h:86
Abstract base class interface required for functors used as input to Functions.
Definition: function_interface.h:68
A multiresolution adaptive numerical function.
Definition: mra.h:122
T trace() const
Returns global value of int(f(x),x) ... global comm required.
Definition: mra.h:1099
double norm2() const
Returns the 2-norm of the function ... global sum ... works in either basis.
Definition: mra.h:688
Function< T, NDIM > & scale(const Q q, bool fence=true)
Inplace, scale the function by a constant. No communication except for optional fence.
Definition: mra.h:953
Function< T, NDIM > & truncate(double tol=0.0, bool fence=true)
Truncate the function with optional fence. Compresses with fence if not compressed.
Definition: mra.h:602
const Function< T, NDIM > & reconstruct(bool fence=true) const
Reconstructs the function, transforming into scaling function basis. Possible non-blocking comm.
Definition: mra.h:775
const Function< T, NDIM > & compress(bool fence=true) const
Compresses the function, transforming into wavelet basis. Possible non-blocking comm.
Definition: mra.h:721
Function< T, NDIM > & gaxpy(const T &alpha, const Function< Q, NDIM > &other, const R &beta, bool fence=true)
Inplace, general bi-linear operation in wavelet basis. No communication except for optional fence.
Definition: mra.h:981
void unaryop(T(*f)(T))
Inplace unary operation on function values.
Definition: mra.h:895
Definition: eigsolver.h:63
double weight()
Definition: eigsolver.h:78
Definition: potentialmanager.h:55
Convolutions in separated form (including Gaussian)
Definition: operator.h:136
A slice defines a sub-range or patch of a dimension.
Definition: slice.h:103
The main class of the periodic DFT solver.
Definition: solver.h:619
Tensor< valueT > tensorT
Definition: solver.h:632
rfunctionT make_lda_potential(World &world, const rfunctionT &arho, const rfunctionT &brho, const rfunctionT &adelrhosq, const rfunctionT &bdelrhosq)
Definition: solver.h:1230
vecfuncT compute_residual(const vecfuncT &awfs, const vecfuncT &bwfs)
Definition: solver.h:3683
std::vector< functionT > vecfuncT
Definition: solver.h:625
int _it
Definition: solver.h:746
void apply_hf_exchange3(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb, double &xc)
Definition: solver.h:2318
void apply_hf_exchange4(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb, double &xc)
Definition: solver.h:2366
World & _world
Definition: solver.h:639
void do_rhs_simple(vecfuncT &wf, vecfuncT &vwf, std::vector< KPoint > kpoints, std::vector< T > &alpha, std::vector< double > &eigs)
Definition: solver.h:3202
MolecularEntity _mentity
Definition: solver.h:714
std::vector< KPoint > genkmesh(unsigned int ngridk0, unsigned ngridk1, unsigned int ngridk2, double koffset0, double koffset1, double koffset2, double R)
Definition: solver.h:890
void init(const std::string &filename)
Definition: solver.h:789
cvecfuncT _aobasisf
Definition: solver.h:750
ctensorT csqrt(const ctensorT &s, double tol=1e-8)
Computes matrix square root (not used any more?)
Definition: solver.h:3634
std::vector< subspaceT > vecsubspaceT
Definition: solver.h:636
vecfuncT _phisa
Definition: solver.h:650
Function< valueT, NDIM > functionT
Definition: solver.h:624
ctensorT matrix_exponential(const ctensorT &A)
Definition: solver.h:2786
Solver(World &world, rfunctionT vnucrhon, vecfuncT phisa, vecfuncT phisb, std::vector< T > eigsa, std::vector< T > eigsb, ElectronicStructureParams params, MolecularEntity mentity)
Definition: solver.h:1289
AtomicBasisSet _aobasis
Definition: solver.h:722
int _nao
Definition: solver.h:754
tensor_complex make_kinetic_matrix(World &world, const vector_complex_function_3d &v, const KPoint &k)
Definition: solver.h:3531
rfunctionT _vnucrhon
Definition: solver.h:646
SeparatedConvolution< T, NDIM > * _cop
Definition: solver.h:698
std::shared_ptr< operatorT> poperatorT
Definition: solver.h:629
void set_occs2(const std::vector< KPoint > &kpoints, const std::vector< double > &eigsa, const std::vector< double > &eigsb, std::vector< double > &occsa, std::vector< double > &occsb)
Definition: solver.h:1893
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:633
std::ofstream _matF
Definition: solver.h:734
std::complex< T > valueT
Definition: solver.h:621
std::vector< T > _eigsa
Definition: solver.h:658
std::vector< KPoint > _kpoints
Definition: solver.h:670
vecfuncT project_ao_basis(World &world, KPoint kpt)
Definition: solver.h:1243
std::vector< tensorT > vectensorT
Definition: solver.h:635
void print_fock_matrix_eigs(const vecfuncT &wf, const vecfuncT &vwf, KPoint kpoint)
Definition: solver.h:2872
std::vector< double > _occsb
Definition: solver.h:678
double calculate_kinetic_energy()
Definition: solver.h:2135
void END_TIMER(World &world, const char *msg)
Definition: solver.h:765
FunctionFactory< T, NDIM > rfactoryT
Definition: solver.h:623
KPoint find_kpt_from_orb(unsigned int idx)
Definition: solver.h:2353
rfunctionT _vnuc
Definition: solver.h:694
ElectronicStructureParams _params
Definition: solver.h:666
Vector< double, NDIM > kvecT
Definition: solver.h:627
std::vector< double > _occsa
Definition: solver.h:674
void print_potential_matrix_eigs(const vecfuncT &wf, const vecfuncT &vwf)
Definition: solver.h:2843
void save_orbitals()
Definition: solver.h:921
void step_restriction(vecfuncT &owfs, vecfuncT &nwfs, int aorb)
Definition: solver.h:3778
FunctionFactory< valueT, NDIM > factoryT
Definition: solver.h:626
void reproject()
Definition: solver.h:2474
void START_TIMER(World &world)
Definition: solver.h:761
std::vector< poperatorT > make_bsh_operators(const std::vector< T > &eigs)
Definition: solver.h:2081
std::ofstream _outputF
Definition: solver.h:730
double _residual
Definition: solver.h:718
void initial_guess()
Initializes alpha and beta mos, occupation numbers, eigenvalues.
Definition: solver.h:1316
virtual ~Solver()
Definition: solver.h:1882
void test_periodicity(const Function< Q, 3 > &f)
Definition: solver.h:1159
void make_nuclear_potential_impl()
Definition: solver.h:1034
Solver(World &world, rfunctionT vnucrhon, vecfuncT phis, std::vector< T > eigs, std::vector< KPoint > kpoints, std::vector< double > occs, ElectronicStructureParams params, MolecularEntity mentity)
Definition: solver.h:1847
rfunctionT _rhob
Definition: solver.h:686
std::vector< T > _eigsb
Definition: solver.h:662
tensorT build_fock_matrix(vecfuncT &psi, vecfuncT &vpsi, KPoint kpoint)
Definition: solver.h:3564
void orthonormalize(vecfuncT &wf, KPoint kpoint)
Definition: solver.h:3660
Solver(World &world, const rfunctionT &vnucrhon, const vecfuncT &phis, const std::vector< T > &eigs, const ElectronicStructureParams &params, MolecularEntity mentity)
Definition: solver.h:1812
rfunctionT compute_rho(const vecfuncT &phis, std::vector< KPoint > kpoints, std::vector< double > occs)
Definition: solver.h:2042
Tensor< double > rtensorT
Definition: solver.h:630
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:634
SeparatedConvolution< T, 3 > operatorT
Definition: solver.h:628
std::ofstream _eigF
Definition: solver.h:738
Tensor< std::complex< double > > ctensorT
Definition: solver.h:631
void apply_hf_exchange(vecfuncT &phisa, vecfuncT &phisb, vecfuncT &funcsa, vecfuncT &funcsb)
Definition: solver.h:2406
Subspace< T, NDIM > * _subspace
Definition: solver.h:710
rfunctionT _rho
Definition: solver.h:690
double sss
Definition: solver.h:760
void make_nuclear_potential()
Definition: solver.h:1094
void make_nuclear_charge_density_impl()
Definition: solver.h:1042
rfunctionT _rhoa
Definition: solver.h:682
void gram_schmidt(vecfuncT &f, KPoint kpoint)
Definition: solver.h:3607
vecfuncT _phisb
Definition: solver.h:654
double _maxthresh
Definition: solver.h:726
void solve()
Definition: solver.h:2513
Function< T, NDIM > rfunctionT
Definition: solver.h:622
void load_orbitals()
Definition: solver.h:942
void do_rhs(vecfuncT &wf, vecfuncT &vwf, std::vector< KPoint > kpoints, std::vector< T > &alpha, std::vector< double > &eigs)
Definition: solver.h:2939
tensorT Q3(const tensorT &s)
Given overlap matrix, return rotation with 3rd order error to orthonormalize the vectors.
Definition: solver.h:3624
double ttt
Definition: solver.h:760
void fix_occupations(const std::vector< T > &eps, std::vector< double > &occs)
Definition: solver.h:3810
void update_orbitals(vecfuncT &awfs, vecfuncT &bwfs, std::vector< KPoint > kpoints)
Definition: solver.h:3716
void print_tensor2d(ostream &os, Tensor< Q > t)
Definition: solver.h:2827
std::ofstream _kmeshF
Definition: solver.h:742
Solver(World &world, const std::string &filename)
Definition: solver.h:771
The SubspaceK class is a container class holding previous orbitals and residuals.
Definition: solver.h:193
Tensor< valueT > tensorT
Definition: solver.h:200
World & _world
Definition: solver.h:205
vectensorT _Q
Definition: solver.h:209
Function< valueT, NDIM > functionT
Definition: solver.h:196
void update_subspace(vecfuncT &awfs_new, vecfuncT &bwfs_new, const vecfuncT &awfs_old, const vecfuncT &bwfs_old)
Definition: solver.h:245
std::vector< KPoint > _kpoints
Definition: solver.h:217
std::complex< T > valueT
Definition: solver.h:195
vecsubspaceT _subspace
Definition: solver.h:213
std::vector< tensorT > vectensorT
Definition: solver.h:201
std::vector< subspaceT > vecsubspaceT
Definition: solver.h:202
SubspaceK(World &world, const ElectronicStructureParams &params, const std::vector< KPoint > &kpoints)
Definition: solver.h:231
std::vector< functionT > vecfuncT
Definition: solver.h:197
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:199
ElectronicStructureParams _params
Definition: solver.h:221
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:198
double _residual
Definition: solver.h:225
The SubspaceK class is a container class holding previous orbitals and residuals.
Definition: solver.h:409
tensorT _Q
Definition: solver.h:423
std::vector< KPoint > _kpoints
Definition: solver.h:431
std::pair< vecfuncT, vecfuncT > pairvecfuncT
Definition: solver.h:414
ElectronicStructureParams _params
Definition: solver.h:435
Subspace(World &world, const ElectronicStructureParams &params)
Definition: solver.h:441
std::vector< pairvecfuncT > subspaceT
Definition: solver.h:415
subspaceT _subspace
Definition: solver.h:427
World & _world
Definition: solver.h:419
std::vector< functionT > vecfuncT
Definition: solver.h:413
Tensor< valueT > tensorT
Definition: solver.h:416
void update_subspace(vecfuncT &awfs_new, vecfuncT &bwfs_new, const vecfuncT &awfs_old, const vecfuncT &bwfs_old, const vecfuncT &rm)
Definition: solver.h:448
Function< valueT, NDIM > functionT
Definition: solver.h:412
std::complex< T > valueT
Definition: solver.h:411
void reproject()
Definition: solver.h:556
scalar_type absmax(long *ind=0) const
Return the absolute maximum value (and if ind is non-null, its index) in the Tensor.
Definition: tensor.h:1754
T * ptr()
Returns a pointer to the internal data.
Definition: tensor.h:1824
float_scalar_type normf() const
Returns the Frobenius norm of the tensor.
Definition: tensor.h:1726
IsSupported< TensorTypeData< Q >, Tensor< T > & >::type scale(Q x)
Inplace multiplication by scalar of supported type (legacy name)
Definition: tensor.h:686
Definition: solver.h:120
double ky
Definition: solver.h:124
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition: solver.h:158
std::vector< coord_3d > specialpt
Definition: solver.h:126
double kx
Definition: solver.h:124
double_complex operator()(const coord_3d &x) const
Definition: solver.h:140
double L
Definition: solver.h:123
WSTAtomicBasisFunctor(const AtomicBasisFunction &aofunc, double L, double kx, double ky, double kz)
Definition: solver.h:128
const AtomicBasisFunction aofunc
Definition: solver.h:122
double kz
Definition: solver.h:124
~WSTAtomicBasisFunctor()
Definition: solver.h:138
void broadcast_serializable(objT &obj, ProcessID root)
Broadcast a serializable object.
Definition: worldgop.h:754
void fence(bool debug=false)
Synchronizes all processes in communicator AND globally ensures no pending AM or tasks.
Definition: worldgop.cc:161
void broadcast(void *buf, size_t nbyte, ProcessID root, bool dowork=true, Tag bcast_tag=-1)
Broadcasts bytes from process root while still processing AM & tasks.
Definition: worldgop.cc:173
void sum(T *buf, size_t nelem)
Inplace global sum while still processing AM & tasks.
Definition: worldgop.h:870
A parallel world class.
Definition: world.h:132
ProcessID rank() const
Returns the process rank in this World (same as MPI_Comm_rank()).
Definition: world.h:318
ProcessID size() const
Returns the number of processes in this World (same as MPI_Comm_size()).
Definition: world.h:328
WorldGopInterface & gop
Global operations.
Definition: world.h:205
An archive for storing local or parallel data, wrapping a BinaryFstreamInputArchive.
Definition: parallel_archive.h:366
An archive for storing local or parallel data wrapping a BinaryFstreamOutputArchive.
Definition: parallel_archive.h:321
static const double R
Definition: csqrt.cc:46
double(* f1)(const coord_3d &)
Definition: derivatives.cc:55
double(* f3)(const coord_3d &)
Definition: derivatives.cc:57
double(* f2)(const coord_3d &)
Definition: derivatives.cc:56
const double delta
Definition: dielectric_external_field.cc:119
bool is_equal(double val1, double val2, double eps)
Definition: esolver.h:190
std::shared_ptr< FunctionFunctorInterface< double, 3 > > rfunctorT
Definition: esolver.h:46
Correspondence between C++ and Fortran types.
const double m
Definition: gfit.cc:199
auto T(World &world, response_space &f) -> response_space
Definition: global_functions.cc:34
rfunctionT compute_rho_slow(const vecfuncT &phis, std::vector< KPoint > kpoints, std::vector< double > occs)
Compute the electronic density for either a molecular or periodic system.
Definition: solver.h:2010
T stheta_fd(const T &x)
Definition: solver.h:75
void apply_potential(vecfuncT &pfuncsa, vecfuncT &pfuncsb, const vecfuncT &phisa, const vecfuncT &phisb, const rfunctionT &rhoa, const rfunctionT &rhob, const rfunctionT &rho)
Applies the LDA effective potential to each orbital. Currently only lda and spin-polarized is not imp...
Definition: solver.h:2177
Tensor< T > conj_transpose(const Tensor< T > &t)
Returns a new deep copy of the complex conjugate transpose of the input tensor.
Definition: tensor.h:2027
Tensor< typename Tensor< T >::scalar_type > arg(const Tensor< T > &t)
Return a new tensor holding the argument of each element of t (complex types only)
Definition: tensor.h:2502
double psi(const Vector< double, 3 > &r)
Definition: hatom_energy.cc:78
static const double v
Definition: hatom_sf_dirac.cc:20
static const long oi
Definition: he.cc:16
Tensor< double > op(const Tensor< double > &x)
Definition: kain.cc:508
#define max(a, b)
Definition: lda.h:51
#define MADNESS_EXCEPTION(msg, value)
Macro for throwing a MADNESS exception.
Definition: madness_exception.h:119
#define MADNESS_ASSERT(condition)
Assert a condition that should be free of side-effects since in release builds this might be a no-op.
Definition: madness_exception.h:134
Main include file for MADNESS and defines Function interface.
const double pi
Mathematical constant .
Definition: constants.h:48
File holds all helper structures necessary for the CC_Operator and CC2 class.
Definition: DFParameters.h:10
@ BC_PERIODIC
Definition: funcdefaults.h:56
Function< TENSOR_RESULT_TYPE(L, R), NDIM > sub(const Function< L, NDIM > &left, const Function< R, NDIM > &right, bool fence=true)
Same as operator- but with optional fence and no automatic compression.
Definition: mra.h:1971
static const char * filename
Definition: legendre.cc:96
std::vector< complex_function_3d > vector_complex_function_3d
Definition: functypedefs.h:86
double abs(double x)
Definition: complexfun.h:48
static double cpu_time()
Returns the cpu time in seconds relative to an arbitrary origin.
Definition: timers.h:127
Vector< double, 3 > coordT
Definition: corepotential.cc:54
void sygv(const Tensor< T > &A, const Tensor< T > &B, int itype, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Generalized real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:1052
response_space scale(response_space a, double b)
Function< T, NDIM > conj(const Function< T, NDIM > &f, bool fence=true)
Return the complex conjugate of the input function with the same distribution and optional fence.
Definition: mra.h:2046
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d >>> &op, response_space &f)
Definition: basic_operators.cc:39
Function< double, NDIM > abssq(const Function< double_complex, NDIM > &z, bool fence=true)
Returns a new function that is the square of the absolute value of the input.
Definition: mra.h:2713
void truncate(World &world, response_space &v, double tol, bool fence)
Definition: basic_operators.cc:30
Function< T, NDIM > project(const Function< T, NDIM > &other, int k=FunctionDefaults< NDIM >::get_k(), double thresh=FunctionDefaults< NDIM >::get_thresh(), bool fence=true)
Definition: mra.h:2399
Tensor< T > KAIN(const Tensor< T > &Q, double rcond=1e-12)
Solves non-linear equation using KAIN (returns coefficients to compute next vector)
Definition: solvers.h:98
static double r2(const coord_3d &x)
Definition: smooth.h:45
Function< T, NDIM > copy(const Function< T, NDIM > &f, const std::shared_ptr< WorldDCPmapInterface< Key< NDIM > > > &pmap, bool fence=true)
Create a new copy of the function with different distribution and optional fence.
Definition: mra.h:2002
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition: vmra.h:133
void plotdx(const Function< T, NDIM > &f, const char *filename, const Tensor< double > &cell=FunctionDefaults< NDIM >::get_cell(), const std::vector< long > &npt=std::vector< long >(NDIM, 201L), bool binary=true)
Writes an OpenDX format file with a cube/slice of points on a uniform grid.
Definition: mraimpl.h:3448
double norm2(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norm of a vector of functions.
Definition: vmra.h:851
static SeparatedConvolution< double_complex, 3 > PeriodicHFExchangeOperator(World &world, Vector< double, 3 > args, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition: operator.h:1721
static const Slice _(0,-1, 1)
std::vector< Function< T, NDIM > > rot(const std::vector< Function< T, NDIM > > &v, bool do_refine=false, bool fence=true)
shorthand rot operator
Definition: vmra.h:1976
void print(const T &t, const Ts &... ts)
Print items to std::cout (items separated by spaces) and terminate with a new line.
Definition: print.h:225
Function< T, NDIM > sum(World &world, const std::vector< Function< T, NDIM > > &f, bool fence=true)
Returns new function — q = sum_i f[i].
Definition: vmra.h:1421
NDIM & f
Definition: mra.h:2416
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition: corepotential.cc:55
Function< TENSOR_RESULT_TYPE(L, R), NDIM > mul_sparse(const Function< L, NDIM > &left, const Function< R, NDIM > &right, double tol, bool fence=true)
Sparse multiplication — left and right must be reconstructed and if tol!=0 have tree of norms already...
Definition: mra.h:1742
double wall_time()
Returns the wall time in seconds relative to an arbitrary origin.
Definition: timers.cc:48
std::vector< complex_functionT > cvecfuncT
Definition: SCF.h:80
double inner(response_space &a, response_space &b)
Definition: response_functions.h:442
static SeparatedConvolution< double, 3 > * CoulombOperatorPtr(World &world, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with 1/r in 3D.
Definition: operator.h:1762
void normalize(World &world, std::vector< Function< T, NDIM > > &v, bool fence=true)
Normalizes a vector of functions — v[i] = v[i].scale(1.0/v[i].norm2())
Definition: vmra.h:1614
void plot_line(World &world, const char *filename, int npt, const Vector< double, NDIM > &lo, const Vector< double, NDIM > &hi, const opT &op)
Generates ASCII file tabulating f(r) at npoints along line r=lo,...,hi.
Definition: funcplot.h:438
double real(double x)
Definition: complexfun.h:52
static SeparatedConvolution< double, 3 > * BSHOperatorPtr3D(World &world, double mu, double lo, double eps, const BoundaryConditions< 3 > &bc=FunctionDefaults< 3 >::get_bc(), int k=FunctionDefaults< 3 >::get_k())
Factory function generating separated kernel for convolution with exp(-mu*r)/(4*pi*r) in 3D.
Definition: operator.h:2023
std::vector< Function< TENSOR_RESULT_TYPE(T, R), NDIM > > transform(World &world, const std::vector< Function< T, NDIM > > &v, const Tensor< R > &c, bool fence=true)
Transforms a vector of functions according to new[i] = sum[j] old[j]*c[j,i].
Definition: vmra.h:689
void matrix_inner(DistributedMatrix< T > &A, const std::vector< Function< T, NDIM > > &f, const std::vector< Function< T, NDIM > > &g, bool sym=false)
Definition: distpm.cc:46
std::vector< double > norm2s(World &world, const std::vector< Function< T, NDIM > > &v)
Computes the 2-norms of a vector of functions.
Definition: vmra.h:827
void syev(const Tensor< T > &A, Tensor< T > &V, Tensor< typename Tensor< T >::scalar_type > &e)
Real-symmetric or complex-Hermitian eigenproblem.
Definition: lapack.cc:969
static double phase(long i)
Definition: twoscale.cc:85
void gaxpy(const double a, ScalarResult< T > &left, const double b, const T &right, const bool fence=true)
the result type of a macrotask must implement gaxpy
Definition: macrotaskq.h:140
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition: vmra.h:156
static long abs(long a)
Definition: tensor.h:218
XCfunctional xc
Definition: newsolver_lda.cc:53
double Q(double a)
Definition: relops.cc:20
static const double c
Definition: relops.cc:10
static const double thresh
Definition: rk.cc:45
static const long k
Definition: rk.cc:44
Defines interfaces for optimization and non-linear equation solvers.
Definition: electronicstructureparams.h:46
unsigned int maxsub
Definition: electronicstructureparams.h:86
bool spinpol
Definition: electronicstructureparams.h:56
int nio
Definition: electronicstructureparams.h:98
bool fractional
Definition: electronicstructureparams.h:84
double ncharge
Definition: electronicstructureparams.h:103
double koffset2
Definition: electronicstructureparams.h:94
int ngridk1
Definition: electronicstructureparams.h:78
bool ispotential
Definition: electronicstructureparams.h:62
int restart
Definition: electronicstructureparams.h:101
bool canon
Definition: electronicstructureparams.h:90
int maxits
Definition: electronicstructureparams.h:60
double rcriterion
Definition: electronicstructureparams.h:111
int nempty
Definition: electronicstructureparams.h:72
int nelec
Definition: electronicstructureparams.h:50
double koffset1
Definition: electronicstructureparams.h:94
int nbands
Definition: electronicstructureparams.h:76
std::string basis
Definition: electronicstructureparams.h:96
void read_file(const std::string &filename)
Definition: electronicstructureparams.h:170
int waveorder
Definition: electronicstructureparams.h:66
double swidth
Definition: electronicstructureparams.h:105
double sd
Definition: electronicstructureparams.h:115
double lo
Definition: electronicstructureparams.h:54
double thresh
Definition: electronicstructureparams.h:64
bool plotorbs
Definition: electronicstructureparams.h:109
int functional
Definition: electronicstructureparams.h:52
int ngridk0
Definition: electronicstructureparams.h:78
double L
Definition: electronicstructureparams.h:48
bool centered
Definition: electronicstructureparams.h:113
bool print_matrices
Definition: electronicstructureparams.h:107
int solver
Definition: electronicstructureparams.h:92
bool periodic
Definition: electronicstructureparams.h:58
int ngridk2
Definition: electronicstructureparams.h:78
double koffset0
Definition: electronicstructureparams.h:94
Definition: solver.h:1121
double operator()(const coordT &x) const
Definition: solver.h:1127
const bool periodic
Definition: solver.h:1125
const MolecularEntity & mentity
Definition: solver.h:1122
double R
Definition: solver.h:1124
GuessDensity(const MolecularEntity &mentity, const AtomicBasisSet &aobasis, const double &R, const bool &periodic)
Definition: solver.h:1151
const AtomicBasisSet & aobasis
Definition: solver.h:1123
static const double_complex I
Definition: tdse1d.cc:164
Tensor< double > B
Definition: tdse1d.cc:166
static const double eshift
Definition: tdse1d.cc:153
const double alpha
Definition: test_coulomb.cc:54
void e()
Definition: test_sig.cc:75
static const double ky
Definition: testcosine.cc:16
static const double kz
Definition: testcosine.cc:16
static const double kx
Definition: testcosine.cc:16
static const std::size_t NDIM
Definition: testpdiff.cc:42
double k0
Definition: testperiodic.cc:66
double k2
Definition: testperiodic.cc:68
double k1
Definition: testperiodic.cc:67
FLOAT weight(const FLOAT &x)
Definition: y.cc:309