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>
36#include <vector>
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
48using std::endl;
49using std::setfill;
50using std::setw;
51using 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
67namespace 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
107
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>
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>
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 //*************************************************************************
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 //*************************************************************************
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 //*************************************************************************
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 {
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 {
945 const int k = FunctionDefaults<3>::get_k();
946
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();
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 //*************************************************************************
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
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 {
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
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
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:");
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
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
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 {
1841 }
1842 }
1843 //*************************************************************************
1844
1845 //*************************************************************************
1846 // Constructor
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 {
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
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
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 */
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 //*************************************************************************
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 //*************************************************************************
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 //*************************************************************************
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 //*************************************************************************
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);
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);
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 //*************************************************************************
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,
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);
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:");
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 //*************************************************************************
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);
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
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,
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:");
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);
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 //*************************************************************************
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
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 const Tensor< double > & get_cell_width()
Returns the width of each user cell dimension.
Definition funcdefaults.h:468
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
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
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
std::vector< KPoint > genkmesh(unsigned int ngridk0, unsigned ngridk1, unsigned int ngridk2, double koffset0, double koffset1, double koffset2, double R)
Definition solver.h:890
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
std::vector< poperatorT > make_bsh_operators(const std::vector< T > &eigs)
Definition solver.h:2081
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
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::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
A tensor is a multidimension array.
Definition tensor.h:317
T * ptr()
Returns a pointer to the internal data.
Definition tensor.h:1824
A simple, fixed dimension vector.
Definition vector.h:64
Definition solver.h:120
double ky
Definition solver.h:124
std::vector< coord_3d > specialpt
Definition solver.h:126
std::vector< coord_3d > special_points() const
Override this to return list of special points to be refined more deeply.
Definition solver.h:158
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.
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 rot(x, k)
Definition lookup3.c:72
#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
Namespace for all elements and tools of MADNESS.
Definition DFParameters.h:10
@ BC_PERIODIC
Definition funcdefaults.h:56
static const char * filename
Definition legendre.cc:96
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
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
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
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
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< 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
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 truncate(World &world, response_space &v, double tol, bool fence)
Definition basic_operators.cc:30
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
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
static double r2(const coord_3d &x)
Definition smooth.h:45
void compress(World &world, const std::vector< Function< T, NDIM > > &v, bool fence=true)
Compress a vector of functions.
Definition vmra.h:133
const std::vector< Function< T, NDIM > > & reconstruct(const std::vector< Function< T, NDIM > > &v)
reconstruct a vector of functions
Definition vmra.h:156
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 const Slice _(0,-1, 1)
Tensor< double_complex > tensor_complex
Definition functypedefs.h:44
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
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
response_space apply(World &world, std::vector< std::vector< std::shared_ptr< real_convolution_3d > > > &op, response_space &f)
Definition basic_operators.cc:39
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
NDIM & f
Definition mra.h:2416
std::shared_ptr< FunctionFunctorInterface< double, 3 > > functorT
Definition corepotential.cc:55
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
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:1671
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
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
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
double real(double x)
Definition complexfun.h:52
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
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 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
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 m
Definition relops.cc:9
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
static const double eshift
Definition tdse1d.cc:153
AtomicInt sum
Definition test_atomicint.cc:46
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 double alpha
Definition testcosine.cc:10
#define START_TIMER
Definition testgaxpyext.cc:9
#define END_TIMER(msg)
Definition testgaxpyext.cc:10
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